2014 dxdy logo

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

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




Начать новую тему Ответить на тему
 
 портировать interpn() из Матлаб в С++
Сообщение07.09.2015, 19:38 
Экс-модератор
Аватара пользователя


23/12/05
12046
Еще одна задачка портирования ))

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

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 
Заслуженный участник


15/05/05
3445
USA
photon в сообщении #1051353 писал(а):
Еще одна задачка портирования ))
В матлабе есть: ...
Вы не заглядывали в исходные коды Scilab или Octave?
Это, конечно, не Matlab, но достаточно близко.

 Профиль  
                  
 
 Re: портировать interpn() из Матлаб в С++
Сообщение08.09.2015, 19:20 
Экс-модератор
Аватара пользователя


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

 Профиль  
                  
 
 Re: портировать interpn() из Матлаб в С++
Сообщение08.09.2015, 19:44 
Заслуженный участник
Аватара пользователя


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

 Профиль  
                  
 
 Re: портировать interpn() из Матлаб в С++
Сообщение08.09.2015, 21:41 
Заслуженный участник


15/05/05
3445
USA
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 
Экс-модератор
Аватара пользователя


23/12/05
12046
спасибо

 Профиль  
                  
 
 Re: портировать interpn() из Матлаб в С++
Сообщение09.09.2015, 04:31 
Заслуженный участник


12/07/07
4438
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 
Экс-модератор
Аватара пользователя


23/12/05
12046
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 
Заслуженный участник


12/07/07
4438
Поясню своё предыдущее сообщение кодом.
Скрипт
Используется синтаксис 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 
Заслуженный участник


12/07/07
4438
Перечитал в свободное время.
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 
Экс-модератор
Аватара пользователя


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

 Профиль  
                  
 
 Re: портировать interpn() из Матлаб в С++
Сообщение16.09.2015, 18:04 
Экс-модератор
Аватара пользователя


23/12/05
12046
Всем спасибо, работает. Матлаб обогнали )

 Профиль  
                  
 
 Re: портировать interpn() из Матлаб в С++
Сообщение16.09.2015, 18:16 
Заслуженный участник
Аватара пользователя


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

 Профиль  
                  
 
 Re: портировать interpn() из Матлаб в С++
Сообщение16.09.2015, 18:34 
Экс-модератор
Аватара пользователя


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

 Профиль  
                  
Показать сообщения за:  Поле сортировки  
Начать новую тему Ответить на тему  [ Сообщений: 14 ] 

Модераторы: Karan, Toucan, PAV, maxal, Супермодераторы



Кто сейчас на конференции

Сейчас этот форум просматривают: нет зарегистрированных пользователей


Вы не можете начинать темы
Вы не можете отвечать на сообщения
Вы не можете редактировать свои сообщения
Вы не можете удалять свои сообщения
Вы не можете добавлять вложения

Найти:
Powered by phpBB © 2000, 2002, 2005, 2007 phpBB Group