26/05/12 1792 приходит весна?
|
Последний раз редактировалось B@R5uk 04.04.2025, 01:59, всего редактировалось 5 раз(а).
жёсткие рёбра (минимальной длины) принадлежат выпуклой оболочке вершин. Если обозначить длину ребра l, а расстояние от центра единичной сферы до этого ребра d, то будет выполняться  Это означает, если длина ребра минимально возможная, то расстояние максимально, что с необходимостью приводит к тому, что ребро принадлежит выпуклой оболочке. Потому что в противном случае найдётся пара точек на сфере, расстояние между которыми меньше величины l, что будет противоречить требованию её минимальности. Тут у меня возник другой вопрос: есть ли какой-нибудь алгоритмизируемый способ убедиться, что заданная конструкция жёстких рёбер сама по себе жесткая? То есть, что многогранник задаёт локальный минимум задаче. А то я довольно долго (ошибочно) полагал, что многогранник для 10 точек на рисунке слева ниже является оптимальным решением, и только недавно понял, что он не просто не является оптимальным, он даже локально оптимальным не является. Если разорвать жёсткое ребро 1-2 и, рассматривая угол 1-0-2 как параметр, начать его увеличивать, то весь многогранник поплывёт. При этом длина жёстких рёбер будет расти до тех пор, пока вершины 1 и 2 не упрутся жёсткими рёбрами в соответсвующие вершины под ними (на рисунке слева они соединены парой тонких чёрных линий). В результате получится фигура на рисунке справа. У неё уже 19 жёстких рёбер (на одно больше, чем у предыдущей) и она доставляет задаче глобальный максимум. Можно ли как-то это (локальную оптимальность) проверить с помощью систематического расчёта? Может быть через матрицу производных как-нибудь? -------------------------Код для расчёта картинок причесал, теперь не стыдно показать. Если кому вдруг тоже будет интересно поиграться. Файл расчёта координат многогранника для 10 точек (оптимальное решение) calc_10.m:
% dxdy.ru/topic160089.html
% B@R5uk, 2025, April
clc
clearvars
format compact
points_func = @ (w) [
w(1) 0
w(1) pi
w(2) pi / 2 - w(3)
w(2) pi / 2 + w(3)
pi - w(4) 0
pi - w(5) pi / 2
];
dist_select_func = @ (d) sum (([d(1, 3), d(3, 4), d(3, 5), d(3, 6), d(5, 6)] - d(1, 2)) .^ 2);
target_func = @ (w) dist_select_func (calc_dist_func (calc_sphere_to_cartesian (points_func (w))));
w0 = [0.6 1.4 0.6 1.1 0.4];
fmso = optimset ('fminsearch');
fmso = optimset (fmso, 'Display', 'final', 'TolFun', 1e-15, 'TolX', 1e-15, ...
'MaxIter', 20000, 'MaxFunEvals', 4000);
w = fminsearch (@(w) target_func (w), w0, fmso);
tp = [
w(5) pi / 2
w(5) -pi / 2
w(4) 0
w(4) pi
pi - w(2) pi / 2 - w(3)
pi - w(2) pi / 2 + w(3)
pi - w(2) -pi / 2 - w(3)
pi - w(2) -pi / 2 + w(3)
pi - w(1) 0
pi - w(1) pi
];
num = size (tp, 1);
coords = calc_sphere_to_cartesian (tp);
calc_plot_func (coords)
view ([-24 15])
calc_save_image (num)
Файл расчёта координат многогранника для 10 точек (неверное решение) calc_10_wrong.m:
% dxdy.ru/topic160089.html
% B@R5uk, 2025, April
clc
clearvars
format compact
points_func = @ (w) [
w(1) 0
w(1) 2 * pi / 3
w(2) pi / 3
pi - w(3) 0
pi 0
];
dist_select_func = @ (d) sum (([d(1, 3), d(3, 4), d(4, 5)] - d(1, 2)) .^ 2);
target_func = @ (w) dist_select_func (calc_dist_func (calc_sphere_to_cartesian (points_func (w))));
w0 = [0.6 1.4 0.8];
fmso = optimset ('fminsearch');
fmso = optimset (fmso, 'Display', 'final', 'TolFun', 1e-15, 'TolX', 1e-15, ...
'MaxIter', 20000, 'MaxFunEvals', 4000);
w = fminsearch (@(w) target_func (w), w0, fmso);
tp = [
w(1) 0
w(1) 2 * pi / 3
w(1) -2 * pi / 3
w(2) pi / 3
w(2) -pi / 3
w(2) pi
pi - w(3) 0
pi - w(3) 2 * pi / 3
pi - w(3) -2 * pi / 3
pi 0
];
num = size (tp, 1);
coords = calc_sphere_to_cartesian (tp);
calc_plot_func (coords)
view ([-24 15])
calc_save_image ('points_10_wrong')
Файл отрисовки многогранника (наиболее трудоёмкая в плане программирования часть) calc_plot_func.m:
% dxdy.ru/topic160089.html
% B@R5uk, 2025, April
function calc_plot_func (coords, pcolor, hcolor, ecolor)
subplot (111)
% Прозрачность граней
alpha_value = 0.6;
% Цвет граней
if 1 == nargin
pcolor = [0.8, 0.9, 1.0];
end
% Цвет "жёстких" рёбер
if 2 >= nargin
hcolor = 'b';
end
% Цвет "мягких" рёбер
if 3 >= nargin
ecolor = 'k';
end
% Допустимая погрешность для компланарности вершин
tol_plane = 1e-10;
% Допустимая погрешность для длины жёстких рёбер
tol_lngth = 1e-10;
% Количество вершин многогранника
num = size (coords, 1);
% Центр масс вершин для определения внутренней и внешней стороны каждой грани
center = mean (coords);
% Матрица расстояний между каждой парой вершин (бесконечность на главной диагонали матрицы)
dists = calc_dist_func (coords);
% Минимальное расстояние между точками для определения жестких рёбер
min_dist = min (dists (:));
plot3 (0, 0, 0)
hold on
% Массив для хранения индексов вершин текущей грани
points_list = zeros (num, 1);
% Список кодов троек вершин, проверку которых необходимо пропустить,
% потому что они уже входят в грань с числом вершин четыре и более
skip_list = [];
skip_code_vector = [num ^ 2 num 1];
% Массив флагов рёбер для отрисовки
edge_flags = zeros (num);
% Построение многогранника как выпуклой оболочки его вершин
for k = 1 : num - 2
% Направление на центр
cntr_vect = center - coords (k, :);
for l = k + 1 : num - 1
for m = l + 1 : num
% Пропустить тройку, если она в списке
if 0 ~= sum (num * (num * k + l) + m == skip_list)
continue
end
% Нормаль к текущей грани
norm_vect = cross_product (coords (l, :) - coords (k, :), coords (m, :) - coords (k, :))';
tmp_value = cntr_vect * norm_vect;
% Количество вершин текущей грани
points_count = 0;
% Флаг, помечающий, что текущая тройка вершин не принадлежит внешней грани
flag = 0;
% Цикл ищет вершины в плоскости текущей тройки точек и
% проверяет, что все остальные лежат по ту же сторону от неё,
% что и центр масс вершин многогранника
for n = 1 : num
tmp_vect = coords (n, :) - coords (k, :);
% Для вершин в плоскости текущей тройки точек скалярное
% произведение нормали на направление на вершину должно
% равняться почти нулю
if abs (tmp_vect * norm_vect) < tol_plane
points_count = points_count + 1;
points_list (points_count) = n;
continue
end
% Скалярные произведения нормали на вектора направлений
% на вершину и на центр масс должны быть одного знака
if 0 >= tmp_value * tmp_vect * norm_vect;
flag = 1;
break
end
end
if flag
continue
end
% Координаты вершин грани
face_coords = coords (points_list (1 : points_count), :);
% Для числа вершин грани четыре и более требуется
% нахождение выпуклой оболочки вершин в плоскости грани
if 3 == points_count
% Грань является треугольником
new_list = [1 2 3];
else
% Затравка для построения -- угол наибольшей величины
new_list = max_angle (face_coords);
new_list (points_count) = 0;
% Обход всех вершин
for n = 3 : points_count
new_list (n) = max_angle_next (face_coords, new_list (n - 2), new_list (n - 1));
end
% Пополняем список троек вершин к пропуску
skip_list = [skip_list, list_of_codes(points_list (1 : points_count), skip_code_vector, 3)]; %#ok<AGROW>
end
% Отрисовка многоугольника с вершинами в правильном порядке
fill3 (face_coords (new_list, 1), ...
face_coords (new_list, 2), ...
face_coords (new_list, 3), pcolor, 'EdgeAlpha', 0, 'FaceAlpha', alpha_value)
% Помечаем рёбра для отрисовки
for n = 2 : points_count
edge_flags (points_list (new_list (n - 1)), points_list (new_list (n))) = 1;
edge_flags (points_list (new_list (n)), points_list (new_list (n - 1))) = 1;
end
edge_flags (points_list (new_list (1)), points_list (new_list (points_count))) = 1;
edge_flags (points_list (new_list (points_count)), points_list (new_list (1))) = 1;
end
end
end
% Количество жёстких рёбер
num_hard = 0;
% Отрисовка рёбер
for k = 1 : num - 1
for l = k + 1 : num
% Жесткие рёбра
if dists (k, l) - min_dist < tol_lngth
num_hard = num_hard + 1;
plot3 (coords ([k l], 1), coords ([k l], 2), coords ([k l], 3), hcolor, 'LineWidth', 3)
continue
end
% Мягкие рёбра
if edge_flags (k, l)
plot3 (coords ([k l], 1), coords ([k l], 2), coords ([k l], 3), ecolor)
end
end
end
hold off
xlim ([-1, 1])
ylim ([-1, 1])
zlim ([-1, 1])
axis vis3d
axis off
set (gcf, 'Color', 'w')
disp (['Number of points: ', num2str(num)])
disp (['Number of hard edges: ', num2str(num_hard)])
disp (['Minimal distance: ', num2str(min_dist, 15)])
disp (['Spherical distance: ', num2str(2 * asind (min_dist / 2), 15), ' deg'])
disp (char (13))
end
function result = cross_product (a, b)
result = [a(2) * b(3) - a(3) * b(2), a(3) * b(1) - a(1) * b(3), a(1) * b(2) - a(2) * b(1)];
end
function result = max_angle (coords)
num = size (coords, 1);
min_value = 1;
result = [];
for m = 2 : num
tmp_vect_1 = coords (m, :) - coords (1, :);
tmp_vect_1 = tmp_vect_1' / sqrt (sum (tmp_vect_1 .^ 2));
for n = 2 : num
% Максимальный угол соответсвует наименьшему скалярному
% произведению нормированных векторов
tmp_vect_2 = coords (n, :) - coords (1, :);
tmp_vect_2 = tmp_vect_2 / sqrt (sum (tmp_vect_2 .^ 2));
value = tmp_vect_2 * tmp_vect_1;
if min_value > value
min_value = value;
result = [m 1 n];
end
end
end
end
function result = max_angle_next (coords, k, l)
num = size (coords, 1);
min_value = 1;
result = [];
tmp_vect_1 = coords (k, :) - coords (l, :);
tmp_vect_1 = tmp_vect_1' / sqrt (sum (tmp_vect_1 .^ 2));
for m = 1 : num
tmp_vect_2 = coords (m, :) - coords (l, :);
tmp_vect_2 = tmp_vect_2 / sqrt (sum (tmp_vect_2 .^ 2));
value = tmp_vect_2 * tmp_vect_1;
if min_value > value
min_value = value;
result = m;
end
end
end
function result = list_of_codes (list, vector, num)
result = vector * nchoosek (list, num)';
end
Файл сохраняющий картинку в файл calc_save_image.m:
% dxdy.ru/topic160089.html
% B@R5uk, 2025, April
function calc_save_image (indx)
color_map = [
255 255 255
0 0 0
0 0 255
128 144 160
208 240 255
128 144 255
] / 255;
pause (0.5)
frmdata = getframe (gcf ());
imgdata = rgb2ind (frmdata .cdata, color_map, 'nodither');
[height, width] = size (imgdata);
imgdata = imcrop (imgdata, [floor(width / 2) - 128, floor(height / 2) - 128, 256, 256]);
if isnumeric (indx)
filename = ['points_', num2str(indx), '.gif'];
else
filename = [indx, '.gif'];
end
imwrite (imgdata, color_map, filename, 'GIF', 'TransparentColor', 0)
end
Файл вспомогательной функции (этот и следующий) calc_sphere_to_cartesian .m: % dxdy.ru/topic160089.html
% B@R5uk, 2025, April
function result = calc_sphere_to_cartesian (tp)
result = [sin(tp (:, 1)) .* cos(tp (:, 2)) sin(tp (:, 1)) .* sin(tp (:, 2)) cos(tp (:, 1))];
end
Файл calc_dist_func .m: % dxdy.ru/topic160089.html
% B@R5uk, 2025, April
function result = calc_dist_func (coords)
num = size (coords, 1);
result = repmat (permute (coords, [3, 1, 2]), num, 1);
result = sqrt (sum ((permute (result, [2, 1, 3]) - result) .^ 2, 3));
result = result + diag (inf (num, 1));
end
|
|