Комбинации элементов - перебрать по 3 [Matlab]

[d_integral]

Не знаю, куда правильнее отнести такой вопрос, к математике, программированию, или еще чему-то, поэтому, если что, прошу прощения за неправильное размещение.
Задача состоит в том, что нужно составить алгоритм, который бы перебирал все комбинации по 3 из 232-х элементов. То есть у меня есть 232 объекта. Для каждой "тройки" (трех не одинаковых) можно посчитать некое число.
У каждого объекта есть индекс ( i=1:232 ).
И вот как бы мне их так перебрать, чтобы сконструировать массив (набор) этих самых чисел, не важно в каком порядке они будут записаны в этот массив, главное, чтобы туда вошли числа, соответствующие всем возможным "тройкам".

Используется синтаксис Matlab M

m = 1;
for i = 1:232
for j = 1:232
for k = 1:232
if i ~= j && i ~= k && j ~= k 
...G(m) = ... ;
m=m+1;
end
end
end
end
 

Это первое, что пришло в голову, но туда войдут повторяющиеся тройки (то есть как 2-3-5, так и 3-2-5), а как сделать так, чтобы этого не было?
Если вопрос совсем дурацкий, ткните меня в матчасть, или посоветуйте, что поискать в гугле :)
Заранее благодарю за совет.

[Pphantom]

Все банально:

Используется синтаксис Matlab M

for i = 1:230
for j = (i+1):231
for k = (j+1):232
... (и никаких if'ов)
 

[arseniiv]

Кстати!
d_integral, перебрать подобным способом $n$-ки сможете? Число $n$ не известно заранее.

[d_integral]

Pphantom
да, спасибо большое! Действительно все просто :) Как раз пришла в голову такая идея, и здорово, что оно без особенных трудов реализуется. Эх, не было у нас в школе комбинаторики..

Цитата:

перебрать подобным способом $n$-ки сможете? Число $n$ не известно заранее.

То есть когда не только "тройки" нужно взять, но и "пары", "четверки" и прочие вплоть до n=N=232 ?
Ну, наверное, n вложенных циклов? что-то вроде

Используется синтаксис Matlab M

for i1 = 1:(232-n)
for i2 = (i1+1):(231-n+1)
...

не?
Правда, писать много будет...

[Pphantom]

d_integral в сообщении #932741 писал(а):

Правда, писать много будет...

Да, особенно 232 вложенных цикла.

[Alexu007]

1. Отсортируйте свой массив из 232-х элементов.
2. Выкиньте повторяющиеся элементы.
3. Переберите тройки:

0-1-2, 0-1-3, 0-1-4 ...
0-2-3, 0-2-4, 0-2-5 ...
...
1-2-3, 1-2-4 ...

и так далее до конца. Это можно сделать с помощью вложенных циклов.

[d_integral]

Alexu007, к сожалению, не могу сказать, что полностью понимаю алгоритм, но примерно ясно..

А пока, если можно, я тут же, чтобы новую тему не создавать, задам короткий вопросик по поводу этой же программы на счет ее возможной оптимизации - просто у меня как-то уж больно медленно считает. Вот такой код

Используется синтаксис Matlab M

m = 1;
for i1 = 1:230
    for i2 = (i1+1):231
        for i3 = (i2+1):232
            R12 = ((Hy(i1).X-Hy(i2).X)^2 + (Hy(i1).Y-Hy(i2).Y)^2 + (Hy(i1).Z-Hy(i2).Z)^2)^0.5;
            R23 = ((Hy(i2).X-Hy(i3).X)^2 + (Hy(i2).Y-Hy(i3).Y)^2 + (Hy(i2).Z-Hy(i3).Z)^2)^0.5;
            R13 = ((Hy(i1).X-Hy(i3).X)^2 + (Hy(i1).Y-Hy(i3).Y)^2 + (Hy(i1).Z-Hy(i3).Z)^2)^0.5;
            if R12 >10 && R23 >10 && R13 >10
plane1 = Hy(i1).Y*(Hy(i2).Z - Hy(i3).Z) + Hy(i2).Y*(Hy(i3).Z - Hy(i1).Z) + Hy(i3).Y*(Hy(i1).Z - Hy(i2).Z) ;
plane2 = Hy(i1).Z*(Hy(i2).X - Hy(i3).X) + Hy(i2).Z*(Hy(i3).X - Hy(i1).X) + Hy(i3).Z*(Hy(i1).X - Hy(i2).X) ;
plane3 = Hy(i1).X*(Hy(i2).Y - Hy(i3).Y) + Hy(i2).X*(Hy(i3).Y - Hy(i1).Y) + Hy(i3).X*(Hy(i1).Y - Hy(i2).Y) ;
% plane4 = -(Hy(i1).X*(Hy(i2).Y*Hy(i3).Z - Hy(i3).Y*Hy(i2).Z) + Hy(i2).X*(Hy(i3).Y*Hy(i1).Z - Hy(i1).Y*Hy(i3).Z) + Hy(i3).X*(Hy(i1).Y*Hy(i2).Z - Hy(i2).Y*Hy(i1).Z));
    
                ang = plane3/((plane1)^2 + (plane2)^2 + (plane3)^2)^0.5;
                if ang > 0
                    ang_zenith(m) = acos(ang);
                    m = m+1;
                end
            end
        end
    end
end

 

То есть Hy - это массив точек с координатами, для каждых трех я провожу плоскость через них, и определяю вектор нормали. Если он направлен вверх, то записываю этот угол в массив.
Но как-то оно плохо считается и получается NaN..
Может, опечатка где?

[Alexu007]

А чё тут непонятного?

Возведение в степень 0,5 - это квадратный корень, если не ошибаюсь? Квадратный корень из отрицательного числа и будет NaN - т.е. ошибка. А число может получиться отрицательным? Да легко, если у вас присутствуют и сложения и вычитания.

plane3/((plane1)^2 - plane1 может принимать значение 0? Наверняка может, там у вас куча сложений и вычитаний, при каких то значениях будет получаться 0 - получите деление на 0 и второй NaN.

[Nemiroff]

Alexu007 в сообщении #932858 писал(а):

Квадратный корень из отрицательного числа и будет NaN - т.е. ошибка.

Не будет.

А чего, структура --- это реально нужно?

[12d3]

Одна оптимизация пришла пока в голову - перед стартом цикла по i3 вычислить R12 и если оно меньше десяти, то третий цикл вообще не стартовать.

-- Вт ноя 18, 2014 15:48:28 --

И еще, упрощение вычисления нормали. Нормаль (ненормированная) - это будет векторное произведение $\overrightarrow{R_{12}} \times \overrightarrow{R_{13}}$, с точностью до знака (это смотря в каком порядке вы обходите точки). То есть, вам нужно вычислить всего одну компоненту векторного произведения и определить, положительна она или отрицательна.

[d_integral]

Alexu007 в сообщении #932858 писал(а):

А число может получиться отрицательным? Да легко, если у вас присутствуют и сложения и вычитания.

Не, отрицательным не может - там везде, где берется корень, он берется из суммы квадратов. А вот насчет того, может ли оказаться равным нулю... Думаю, что тоже вряд ли. Не могут оказаться равными нулю одновременно все компоненты вектора нормали.

Nemiroff в сообщении #932859 писал(а):

А чего, структура --- это реально нужно?

Не знаю, если честно :) Это у меня после С++, наверное, удобно, когда есть объект и у него поля (характеристики)..
12d3, спасибо за совет, да, это наверняка упростит дело, сейчас попробую!
Вообще, если поставить не до 232 цикл, а меньше, то оно считает. И даже выдает значения, и, возможно, верные. Просто очень уж долго.

12d3 в сообщении #932862 писал(а):

И еще, упрощение вычисления нормали. Нормаль (ненормированная) - это будет векторное произведение $\overrightarrow{R_{12}} \times \overrightarrow{R_{13}}$, с точностью до знака (это смотря в каком порядке вы обходите точки). То есть, вам нужно вычислить всего одну компоненту векторного произведения и определить, положительна она или отрицательна.

То есть... В принципе -да, направления по Х и Y меня не интересуют, но нормировать нужно - мне ведь в конечном итоге нужен косинус угла между этой нормалью и осью Z. Причем нужен только в тех случаях, когда он положителен, то есть при углах от 0 до $\pi/2$. Кстати, может быть, по этому критерию тоже можно сократить часть циклов?

Пока программа выглядит вот так

Используется синтаксис Matlab M

m = 1;
for i1 = 1:29
    for i2 = 20:50
        R12 = ((Hy(i1).X-Hy(i2).X)^2 + (Hy(i1).Y-Hy(i2).Y)^2 + (Hy(i1).Z-Hy(i2).Z)^2)^0.5;
        if R12 >10
        for i3 = (i2+1):60
            R23 = ((Hy(i2).X-Hy(i3).X)^2 + (Hy(i2).Y-Hy(i3).Y)^2 + (Hy(i2).Z-Hy(i3).Z)^2)^0.5;
            R13 = ((Hy(i1).X-Hy(i3).X)^2 + (Hy(i1).Y-Hy(i3).Y)^2 + (Hy(i1).Z-Hy(i3).Z)^2)^0.5;
            if  R23 >10 && R13 >10
plane1 = Hy(i1).Y*(Hy(i2).Z - Hy(i3).Z) + Hy(i2).Y*(Hy(i3).Z - Hy(i1).Z) + Hy(i3).Y*(Hy(i1).Z - Hy(i2).Z) ;
plane2 = Hy(i1).Z*(Hy(i2).X - Hy(i3).X) + Hy(i2).Z*(Hy(i3).X - Hy(i1).X) + Hy(i3).Z*(Hy(i1).X - Hy(i2).X) ;
plane3 = Hy(i1).X*(Hy(i2).Y - Hy(i3).Y) + Hy(i2).X*(Hy(i3).Y - Hy(i1).Y) + Hy(i3).X*(Hy(i1).Y - Hy(i2).Y) ;
    if plane3 > 0
                cosin(m) = plane3/((plane1)^2 + (plane2)^2 + (plane3)^2)^0.5;
                    m = m+1;
                end
            end
        end
    end
end

 

Сейчас буду еще думать, как его улучшить.

А, да. Может, оно станет быстрее, если для R12 (и других двух расстояний) не считать их самих, а только квадраты и сравнивать с числом 100 ?

[Nemiroff]

d_integral в сообщении #932950 писал(а):

Это у меня после С++, наверное, удобно, когда есть объект и у него поля (характеристики)..

То бишь, пользоваться преимуществами Матлаба в Матлабе --- нет?

[d_integral]

Nemiroff в сообщении #932958 писал(а):

То бишь, пользоваться преимуществами Матлаба в Матлабе --- нет?

А что, 3 отдельных массива будет лучше, чем массив структур?

-- 18.11.2014, 20:06 --

На самом деле мне, оказывается, нужен даже не косинус, а квадрат тангенса. В таком случае алгоритм приобрел следующий вид:

Используется синтаксис Matlab M

m = 1;
for i1 = 1:230
    for i2 = (i1+1):231
        R12 = (Hy(i1).X-Hy(i2).X)^2 + (Hy(i1).Y-Hy(i2).Y)^2 + (Hy(i1).Z-Hy(i2).Z)^2;
        if R12>100
        for i3 = (i2+1):232
            R23 = (Hy(i2).X-Hy(i3).X)^2 + (Hy(i2).Y-Hy(i3).Y)^2 + (Hy(i2).Z-Hy(i3).Z)^2;
            R13 = (Hy(i1).X-Hy(i3).X)^2 + (Hy(i1).Y-Hy(i3).Y)^2 + (Hy(i1).Z-Hy(i3).Z)^2;
            if  R23 >100 && R13 >100
                Norm3 = (Hy(i2).X - Hy(i1).X)*(Hy(i3).Y - Hy(i1).Y) - (Hy(i3).Y - Hy(i1).Y)*(Hy(i2).X - Hy(i1).X);
                if Norm3 > 0

Norm1 = (Hy(i2).Y - Hy(i1).Y)*(Hy(i3).Z - Hy(i1).Z) - (Hy(i3).Y - Hy(i1).Y)*(Hy(i2).Z - Hy(i1).Z);
Norm2 = -(Hy(i2).X - Hy(i1).X)*(Hy(i3).Z - Hy(i1).Z) + (Hy(i3).X - Hy(i1).X)*(Hy(i2).Z - Hy(i1).Z);
                    cos2 = Norm3^2/(Norm1^2 + Norm2^2 + Norm3^2);
                    tg2(m) = 1/cos2 - 1;
                    m = m+1;
                end
            end
        end
        end
    end
end
 

[Nemiroff]

d_integral в сообщении #932989 писал(а):

А что, 3 отдельных массива будет лучше, чем массив структур?

Не знаю. А зачем три — у вас же вроде бы один массив. Можно же в матрицу писать, почему нет?

[d_integral]

Найдена опечатка:

Используется синтаксис Matlab M

 Norm3 = (Hy(i2).X - Hy(i1).X)*(Hy(i3).Y - Hy(i1).Y) - (Hy(i2).Y - Hy(i1).Y)*(Hy(i3).X - Hy(i1).X);
 

Но теперь он почему-то не продвигается дальше m=1. То есть вот посчитался Norm3, оказался отрицательным, поэтому в условные циклы не пойдем, а берем следующую тройку... Почему же он останавливается?

Nemiroff в сообщении #933015 писал(а):

Можно же в матрицу писать, почему нет?

Да в общем-то можно, конечно, только я боюсь совсем запутаться в индексах )