портировать interpn() из Матлаб в С++

[photon]

Еще одна задачка портирования ))

В матлабе есть:

A=intepn(B,x,y,z);

B - трехмерный массив (например 4x5x26, но могут быть и другие), x,y,z - двумерные массивы (в моем случае, например, 801х1001).

Такой вызов функции для указанного размера равносилен

A=intepn(1:4,1:5,1:26, B, x,y,z);

x,y - получены при помощи meshgrid() с некоторой последующей обработкой, которая сохранила линейность содержимого в этих массивах, z - некие данные.

Поскольку метод интерполяции не указан - по умолчанию, он линейный

Как конкретно получается каждая из точек A в данном случае?

[Yuri Gendelman]

photon в сообщении #1051353 писал(а):

Еще одна задачка портирования ))
В матлабе есть: ...

Вы не заглядывали в исходные коды Scilab или Octave?
Это, конечно, не Matlab, но достаточно близко.

[photon]

со Scilab не сталкивался, а Octave... ставить виртуалку, ставить Octave, чтобы посмотреть пару функций - как-то не очень хочется. Тем более, что эти функции получаются с дополнительными ненужными мне наворотами под разные опции, разные размерности массивов, от которых надо будет избавляться.

[Munin]

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

[Yuri Gendelman]

photon в сообщении #1051645 писал(а):

...ставить виртуалку, ставить Octave, чтобы посмотреть пару функций...

Я имел в виду только исходники. Их можно скачать как отдельные архивы:
octave-4.0.0.tar.gz (22 МБ)
scilab-5.5.2-src.tar.gz (66 МБ)

P.S. И почему виртуалка? Обе системы есть на Windows и на Linux.

[photon]

спасибо

[GAA]

photon в сообщении #1051353 писал(а):

Как конкретно получается каждая из точек A в данном случае?

Для простоты я буду рассматривать скалярный вариант. Т.е. вычислять значение для одной точки.

Предварительно напомню себе двумерный случай. Координаты точки, для которой необходимо вычислить значение функции, обозначим через $(x_1, x_2)$ [этому соответствуют обозначения x и y начального сообщения ветки].
$j_i = \lfloor x_i \rfloor$, $\Delta_i = x_i - j_i$.
$A = \left[B_{j_1, j_2}(1 - \Delta_1) + B_{j_1+1, j_2} \Delta_1\right] (1- \Delta_2)+\left [B_{j_1, j_2+1}(1 - \Delta_1) + B_{j_1+1, j_2+1} \Delta_1 \right] \Delta_2$.


В трехмерном случае координаты точки, для которой необходимо вычислить значение функции, обозначим через $(x_1, x_2, x_3)$.

$A =$

$ \left\{ \left[B_{j_1, j_2, j_3}(1 - \Delta_1) + B_{j_1+1, j_2, j_3} \Delta_1 \right] (1- \Delta_2) +\left[ B_{j_1, j_2+1, j_3}(1 - \Delta_1) + B_{j_1+1, j_2+1, j_3} \Delta_1] \Delta_2 \right\} (1-\Delta_3) + $

$\left\{ \left[B_{j_1, j_2, j_3+1}(1 - \Delta_1) + B_{j_1+1, j_2, j_3+1} \Delta_1 \right] (1- \Delta_2) + \left[ B_{j_1, j_2+1, j_3+1}(1 - \Delta_1) + B_{j_1+1, j_2+1, j_3+1} \Delta_1\right] \Delta_2 \right\} \Delta_3$.


Ну а в подфункции function F = linear(varargin) функции interpn (MatLab) выполняется линейная интерполяция не обязательно для трехмерного случая. При таком раскладе уже не получится явно выписать все слагаемые. Это приводит к необходимости организации вычислений в виде двойного цикла.

Надеюсь, об этом Вы спрашивали.

[photon]

GAA в сообщении #1051754 писал(а):

В трехмерном случае координаты точки, для которой необходимо вычислить значение функции, обозначим через $(x_1, x_2, x_3)$.

Чегой-то не то...
В моем случае

photon в сообщении #1051353 писал(а):

A=intepn(B,x,y,z);

результат $A$ - двумерный массив, такой же по размеру как массивы $x$, $y$, и $z$.

[GAA]

Поясню своё предыдущее сообщение кодом.
Скрипт

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

B(1,1,1) = 0; B(1, 2, 1) = 1; B(2, 1, 1) = 1; B(2, 2, 1) = 2;
B(1,1,2) = 0; B(1, 2, 2) = 1; B(2, 1, 2) = 1; B(2, 2, 2) = 2;
X1 = 1.1; X2 = 1.2; X3 = 1.3;

A1 = MyInterp3(B, X1, X2, X3);
A2 = interpn(B, X1, X2, X3);
Err = A1 - A2 

m-файл функции MyInterp3

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

function C = MyInterp3(B, x1, x2, x3)
j1 = floor(x1);    j2 = floor(x2);   j3 = floor(x3);
Delta1 = x1 - j1; Delta2 = x2 - j2; Delta3 = x3 - j3;
C11 = (  B(j1, j2, j3)*(1-Delta1)  +  B(j1+1, j2, j3)*Delta1 ) * ( 1 - Delta2 );
C12 = (B(j1, j2+1, j3)*(1-Delta1)  + B(j1+1, j2+1, j3)*Delta1) * Delta2;
C21 = (B(j1, j2, j3+1)*(1-Delta1)  + B(j1+1, j2, j3+1)*Delta1 ) * ( 1 - Delta2 );
C22 = (B(j1, j2+1, j3+1)*(1-Delta1)  + B(j1+1, j2+1, j3+1 )*Delta1 ) * Delta2;
C = (C11 + C12) * (1-Delta3) + (C21 + C22) * Delta3;
 

Результат выполнения: 0. Т.е. функция interpn Matlab и моя функция MyInterp3 вернули одинаковый результат.
Можно попробовать изменить значения x1, x2 и x3; потестить. Возможно в тексте опечатки. Буду признателен за их указание.

photon в сообщении #1051901 писал(а):

результат $A$ - двумерный массив, такой же по размеру как массивы $x$, $y$, и $z$.

Если X1, X2, X3 — одномерные массивы, то результат interpn — одномерный массив. Например, для того же значения B (что и в скрипте выше) и X1 = [1.1, 1.2]; X2 = [1.2, 1.3]; X3 = [1.3, 1.4];
>> A = interpn(B, X1, X2, X3);
вернёт вектор
A = 0.3000 0.5000

Если X1, X2, X3 — матрицы, то результатом будет трехмерный массив. Например
>> [X, Y, Z] = meshgrid(1:0.1:2, 1:0.1:2, 1:0.1:2);
>> A = interpn(B, X, Y, Z);
вернёт трехмерный массив.

Чтобы получить функцию, возвращающую такой результат, достаточно применить в цикле MyInterp3 к каждой точке.

[GAA]

Перечитал в свободное время.

GAA в сообщении #1052138 писал(а):

Если X1, X2, X3 — матрицы, то результатом будет трехмерный массив.

Спать по ночам надо. :) Поправлюсь.
Если X1, X2, X3 — трехмерные массивы, то результатом будет трехмерный массив.
(Аналогично, если X1, X2, X3 — двумерные массивы, то результатом будет двумерный массив. Т.е. размерность результата определяется размерностью массивов задающих сетку.)

Функция для иллюстрации идеи интерполяции (двумерная сетка). (Проверка выхода на границу сетки не делается, проверка выхода за границу не делается. Для Matlab код абсурден, поэтому никакие проверки на входе процедуры не делаются. Короче, это идея для перевода на C. Не стал загромождать код. Все проверки имеет смысл писать уже в коде на C.)

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

function A = MyInterp3(B, X, Y, Z)
[m, n] = size(X)
for i=1:m
 for j=1:n
  A(i, j) = f(B, X(i, j), Y(i, j), Z(i, j));
 end
end

function C = f(B, x1, x2, x3)
j1 = floor(x1);    j2 = floor(x2);   j3 = floor(x3);
Delta1 = x1 - j1; Delta2 = x2 - j2; Delta3 = x3 - j3;
C11 = (  B(j1, j2, j3)*(1-Delta1)  +  B(j1+1, j2, j3)*Delta1 ) * ( 1 - Delta2 );
C12 = (B(j1, j2+1, j3)*(1-Delta1)  + B(j1+1, j2+1, j3)*Delta1) * Delta2;
C21 = (B(j1, j2, j3+1)*(1-Delta1)  + B(j1+1, j2, j3+1)*Delta1 ) * ( 1 - Delta2 );
C22 = (B(j1, j2+1, j3+1)*(1-Delta1)  + B(j1+1, j2+1, j3+1 )*Delta1 ) * Delta2;
C = (C11 + C12) * (1-Delta3) + (C21 + C22) * Delta3;
 

Скрипт для тестирования. :)

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

B(1,1,1) = 0; B(1, 2, 1) = 1; B(2, 1, 1) = 1; B(2, 2, 1) = 2;
B(1,1,2) = 0; B(1, 2, 2) = 1; B(2, 1, 2) = 1; B(2, 2, 2) = 2;

X = [1.1 1.8; 1.2 1.9]; Y=[1.05 1.095; 1.5 1.55]; Z = [1.01 1.15;  1.25 1.35];

format long
A1 = MyInterp3(B, X, Y, Z);
A2 = interpn(B, X, Y, Z);
Err = A1 - A2

Если запустить, то получим совпадение, в пределах погрешности операций с Double. Прим. Значения элементов матриц взяты произвольные, лишь бы внутри сетки находились.

(Если массивы, задающие сетку, будут трехмерные, то можно нарисовать тройной цикл.)

[photon]

Да, сделал, спасибо... Но работает существенно медленнее, чем в матлабе - надо теперь думать, как ускорить

[photon]

Всем спасибо, работает. Матлаб обогнали )

[Munin]

В чём состояло ускорение, можно поинтересоваться?

[photon]

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