2014 dxdy logo

Научный форум dxdy

Математика, Физика, Computer Science, Machine Learning, LaTeX, Механика и Техника, Химия,
Биология и Медицина, Экономика и Финансовая Математика, Гуманитарные науки




 
 портировать interpn() из Матлаб в С++
Сообщение07.09.2015, 19:38 
Аватара пользователя
Еще одна задачка портирования ))

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

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 в данном случае?

 
 
 
 Re: портировать interpn() из Матлаб в С++
Сообщение08.09.2015, 16:46 
photon в сообщении #1051353 писал(а):
Еще одна задачка портирования ))
В матлабе есть: ...
Вы не заглядывали в исходные коды Scilab или Octave?
Это, конечно, не Matlab, но достаточно близко.

 
 
 
 Re: портировать interpn() из Матлаб в С++
Сообщение08.09.2015, 19:20 
Аватара пользователя
со Scilab не сталкивался, а Octave... ставить виртуалку, ставить Octave, чтобы посмотреть пару функций - как-то не очень хочется. Тем более, что эти функции получаются с дополнительными ненужными мне наворотами под разные опции, разные размерности массивов, от которых надо будет избавляться.

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

 
 
 
 Re: портировать interpn() из Матлаб в С++
Сообщение08.09.2015, 21:41 
photon в сообщении #1051645 писал(а):
...ставить виртуалку, ставить Octave, чтобы посмотреть пару функций...
Я имел в виду только исходники. Их можно скачать как отдельные архивы:
octave-4.0.0.tar.gz (22 МБ)
scilab-5.5.2-src.tar.gz (66 МБ)

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

 
 
 
 Re: портировать interpn() из Матлаб в С++
Сообщение08.09.2015, 22:59 
Аватара пользователя
спасибо

 
 
 
 Re: портировать interpn() из Матлаб в С++
Сообщение09.09.2015, 04:31 
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) выполняется линейная интерполяция не обязательно для трехмерного случая. При таком раскладе уже не получится явно выписать все слагаемые. Это приводит к необходимости организации вычислений в виде двойного цикла.

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

 
 
 
 Re: портировать interpn() из Матлаб в С++
Сообщение09.09.2015, 15:33 
Аватара пользователя
GAA в сообщении #1051754 писал(а):
В трехмерном случае координаты точки, для которой необходимо вычислить значение функции, обозначим через $(x_1, x_2, x_3)$.


Чегой-то не то...
В моем случае
photon в сообщении #1051353 писал(а):
A=intepn(B,x,y,z);

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

 
 
 
 Re: портировать interpn() из Матлаб в С++
Сообщение10.09.2015, 00:33 
Поясню своё предыдущее сообщение кодом.
Скрипт
Используется синтаксис 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 к каждой точке.

 
 
 
 Re: портировать interpn() из Матлаб в С++
Сообщение15.09.2015, 15:32 
Перечитал в свободное время.
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. Прим. Значения элементов матриц взяты произвольные, лишь бы внутри сетки находились.

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

 
 
 
 Re: портировать interpn() из Матлаб в С++
Сообщение16.09.2015, 13:56 
Аватара пользователя
Да, сделал, спасибо... Но работает существенно медленнее, чем в матлабе - надо теперь думать, как ускорить

 
 
 
 Re: портировать interpn() из Матлаб в С++
Сообщение16.09.2015, 18:04 
Аватара пользователя
Всем спасибо, работает. Матлаб обогнали )

 
 
 
 Re: портировать interpn() из Матлаб в С++
Сообщение16.09.2015, 18:16 
Аватара пользователя
В чём состояло ускорение, можно поинтересоваться?

 
 
 
 Re: портировать interpn() из Матлаб в С++
Сообщение16.09.2015, 18:34 
Аватара пользователя
во-первых, сначала я смотрел время в дебаг-версии, где нет оптимизации и сам по себе код выполняется медленнее, что было в корне неверно для оценки скорости, во-вторых, распараллелил вычисления между ядрами процессора.

 
 
 [ Сообщений: 14 ] 


Powered by phpBB © 2000, 2002, 2005, 2007 phpBB Group