Оптимизация алгоритма вычисления расстояний между ячейками в частицах

#c #алгоритм #cpp

Есть вопрос оптимизации алгоритма.

Это своеобразное представление кристалической решётки. Прямоугольные параллелепипеды
(частицы), состоящие из точек - ячеек (27 штук), образуют кристалличекую решётку. Задаются
расстояния по x,y,z между частицами (a,b,c) и между ячейками(a',b',c') в частице. Необходимо
раcсчитать расстояние между каждыми двумя ячейками, не входящими в одну частицу.
Реализовал расчёт расстояний через длину вектора, поскольку, поместив начало координат
в какую-либо ячейку, можно вычислить координаты всех ячеек. Но расчёт расстояний по
всем ячейкам занимает очень много времени для образца размером, например, 10*100*100
частиц. Просьба предложить алгоритм, позволяющий производить рассчёт более эффективно.
Кусок кода для рассчёта расстояний с учётом следующего алгоритма:

Расчёт расстояний внутри 1 слоя ZY
Расчёт расстояний от каждой ячейки внутри слоя до каждой ячейки в следующем слое ZY.
И так далее до последнего слоя. 
Потом необходимо сделать комбинацию расстояний, чтобы вычисить полные расстояния.


// В плоскости YZ вычисляем энергию зваимодействия ячеек в одном слое
for (int i = 0; i < 27 * Y * Z; i++) {
    if (i % 27 == 0) {
        part_num++;
    }

for (int j = 0; j < (part_num - 1) * 27; j++) {
        distance =
            sqrt(pow((x_vect.at(i) - x_vect.at(j)), 2) +
            pow((y_vect.at(i) - y_vect.at(j)), 2) + 
            pow((z_vect.at(i) - z_vect.at(j)), 2));
    }

for (int j = (part_num) * 27; j < 27 * Y * Z; j++) {
        distance =
            sqrt(pow((x_vect.at(i) - x_vect.at(j)), 2) +
            pow((y_vect.at(i) - y_vect.at(j)), 2) + 
            pow((z_vect.at(i) - z_vect.at(j)), 2));
    }
}

// Выбираем слой
for (int layer = 1; layer < X; layer++) {
    // В плоскости YZ вычисляем энергию зваимодействия ячеек в разных слоях
    for (int i = 0; i < 27 * Y * Z; i++) {
        for (int j = 27 * Y * Z * layer; j < Z * Y * 27 * (layer + 1); j++) {
            distance =
                sqrt(pow((x_vect.at(i) - x_vect.at(j)), 2) +
                pow((y_vect.at(i) - y_vect.at(j)), 2) + 
                pow((z_vect.at(i) - z_vect.at(j)), 2));
        }
    }

cout << "Layer: " << layer << endl;
}


Прошу помочь с алгоритмом.    

Ответы

Ответ 1

Дистанция между двумя ячейками может быть представлена как
d(i,p,j,q,k,r) = 
    sqrt( (i*a+p*a')*(i*a+p*a') + (j*b+q*b')*(j*b+q*b') + (k*c+r*c')*(kc+rc')),

где 

-n < i < n,  -2<= p<= 2,  -n < j < n,  -2<= q<= 2,  -n < k < n,  -2<= r<= 2;
i, j, k - разность индексов двух частиц по x, y, z соответственно;

p, q, r - разность индексов ячеек двух частиц по x, y, z соответственно при условии
совмещения центров частиц 

При этом очевидно выполняются условия
d(i,p,j,q,k,r) = 
    d(i,p,k,r,j,q) = 
        d(j,q,i,p,k,r) = 
            d(j,q,k,r,i,p) =
                d(k,r,i,p,j,q) = 
                    d(k,r,j,q,i,p)

d(-i,p,j,q,k,r) = d(i,-p,j,q,k,r)
d(i,p,-j,q,k,r) = d(i,p,j,-q,k,r)
d(i,p,j,q,-k,r) = d(i,p,j,q,k,-r)

Последнее означает, что i, j, k можно всегда считать неотрицательными в диапазоне
от 0 до n-1.
Таким образом, для случая решетки n*n*n частиц число всевозможных расстояний
d(i,p,j,q,k,r) = m*(m-1)*(m-2)/6 + m*(m-1) + m = m*(m*m +3*m +2)/6,

где m = 5*n. Если n=50 достаточно вычислить и сохранить 2635500 различных расстояний.
И вместо расчёта расстояний просто в зависимости от разности индексов ячеек в переменных
i, j, k, p, q, r брать нужное расстояние. Правда потребуется определенной порядок хранения
d(i,p,j,q,k,r) и обращения к ним. Если же не заморачиваться с этим, то можно хранить
все m*mm расстояний, т.е. при n=50 это 15625000 чисел, что также терпимо при современных
объемах оперативной памяти.


Ответ 2

К сожалению, приведённая модель вряд ли оптимизируется без дальнейших физических
аппроксимаций. Время обработки (и объём ответа) можно оценить как (27*X*Y*Z)^2 (более
строго, Theta((X*Y*Z)^2)). Рассмотрим образец 10x100x100 (X*Y*Z = 10^5). Тогда объём
ответа (т.е. сумма занимаемой памяти по всем distance), считая в single-precision,
будет примерно 29160 ГБ. Понятно, что одновременно эти данные нужны не будут, так что
проблема только во времени их обработки. На достаточно быстром сервере/кластере это
можно сделать за десяток минут (или меньше), плюс задача хорошо распараллеливается.
Вопрос в том, есть ли у Вас доступ к таким мощностям, или Вы можете себе позволить
дальнейшие аппроксимации, уменьшив количество значений distance, входящих в ответ.


Ответ 3

На мой взгляд, не нужно рассчитывать все расстояния в данной модели, т.к. они будут
дублироваться. Например, диагональных векторов будет 4.


Ответ 4

В данном случае уместно определить минимальное расстояние на котором взаимодействие
частиц приравнивается к нулю.т.к.
"При больших расстояниях между
частицами температура перехода близка к температуре
фазового перехода в изолированной частице, смещен-
ной за счет размерных эффектов"
В этом случае можно посчитать минимальное(27) и максимальное количество частиц, что
даст и ответ на вопрос нужно ли оптимизировать алгоритм.


Ответ 5

Если решетка состоит из строго ортогональных отрезков (как на рисунке) и если длина
отрезков фиксирована, то, если разумно пронумеровать ячейки, можно существенно оптимизировать.
Тогда заранее считаете стандартные длины отрезков и пишите отображение

(индексЯчейки, индексЯчейки) --> расстояние

к которому обращаетесь по мере надобности (оно будет нужное количество раз суммировать
фиксированные заранее рассчитанные длины). 
P.S. Подозреваю, что на рисунке всё сильно упрощено.


Ответ 6

Допустим есть 3х3х20 частиц по 27 ячеек.

Вычисляем расстояния первого слоя.
Потом от первого до второго.
От первого к третьему (от второго к третьему = от  первого к второму).
От первого к четвертому (от второго к четвертому = от первого к третьему, от третьего
к четвертому = от первого до второго).
...

Эффект турбо-топора. =)


Ответ 7

Если Вам нужно сложить энергии, зависящие только от расстояния, то даже собственно
расстояние  можно подсчитать потом вместе с энергией. Просто в цикле перебора всех ячеек
в ячейке памяти, соответствующей d(i,p,j,q,k,r) подсчитывать количество раз, когда
был востребован индекс i,p,j,q,k,r и затем уже, когда весь цикл по парам ячеек будет
пройден, подсчитать и дистанцию и элементарную энергию, соответствующую индексу i,p,j,q,k,r
умножив на количество раз, когда этот индекс был востребован.
Добавление.
Если Вам нужно сложить энергии, зависящие только от расстояния между ячейками, то
послойный вариант также хорош, так как в этом случае общая энергия
E = n*E(0) + (n-1)*E(1) + (n-2)*E(2) +...+E(n-1)

где E(k) - вклад в энергию от слоев, отстоящих на k единиц шага решетки между частицами
по соответствующей координате.