Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

Mathematica large table periodic interpolation

I have a very large table in Mathematica ((dimcub-1)^3 elements) coming from an inverse FFT. I need to use periodic interpolation on this table. Since periodic interpolation requires that the first elements and last elements are equal, I create a new table of dim^3 elements manually and use that in my interpolation. It works but it is ugly/slow and due to my superfluous intermediate table, I hit the memory barrier sooner. Can any one tell me either how to turn my old table into a periodic one somehow by appending elements or use my non periodic table to make a periodic interpolation function? Here is my current piece of code:

mr 1 is the new table:

mr1 = Table[  0. , {i, 1, dimcub}, {j, 1, dimcub}, {k, 1, dimcub}];

Do[Do[  Do[   
      mr1[[m, n, k]] = oldtable[[m, n, k]] ;  , {m, 1, 
       dimcub - 1}]; , {n, 1, dimcub - 1}]; , {k, 1, dimcub - 1}]; 
Do[Do[     mr1[[m, n, dimcub]] =  mr1[[m, n, 1]]; 
  mr1[[m, dimcub, n]] =  mr1[[m, 1, n]];  
  mr1[[dimcub, m, n]] =  mr1[[1, m, n]];     , {m, 1, dimcub - 1}];  
 mr1[[n, dimcub, dimcub]] =  mr1[[n, 1, 1]]; 
 mr1[[dimcub, n, dimcub]] =  mr1[[1, n, 1]];  
 mr1[[dimcub, dimcub, n]] =  mr1[[1, 1, n]]; , {n, 1, dimcub - 1}]; 
mr1[[dimcub, dimcub, dimcub]] = mr1[[1, 1, 1]]; 

Remove[oldtable]; 

myinterpolatingfunction = 
 ListInterpolation[mr1, {{1, dimcub}, {1, dimcub}, {1, dimcub}}, 
  InterpolationOrder -> 1, 
  PeriodicInterpolation -> True];

 Remove[mr1];

myinterpolatingfunction takes much less memory and works perfectly once i remove the older tables. Any help will be greatly appreciated.

like image 732
Hsn Avatar asked May 31 '11 22:05

Hsn


2 Answers

Both Leonid's and Mr.Wizard's answers do too much work. In Leonid's case only the first three lines are necessary. To show this I'll change the last 4 Sets to Equals:

In[65]:= len = 4; oldtable = 
 Partition[Partition[Range[len^3], len], len]

Out[65]= {{{1, 2, 3, 4}, {5, 6, 7, 8}, {9, 10, 11, 12}, {13, 14, 15, 
   16}}, {{17, 18, 19, 20}, {21, 22, 23, 24}, {25, 26, 27, 28}, {29, 
   30, 31, 32}}, {{33, 34, 35, 36}, {37, 38, 39, 40}, {41, 42, 43, 
   44}, {45, 46, 47, 48}}, {{49, 50, 51, 52}, {53, 54, 55, 56}, {57, 
   58, 59, 60}, {61, 62, 63, 64}}}

In[66]:= oldtable[[All, All, -1]] = oldtable[[All, All, 1]];
oldtable[[All, -1, All]] = oldtable[[All, 1, All]];
oldtable[[-1, All, All]] = oldtable[[1, All, All]];
oldtable[[All, -1, -1]] == oldtable[[All, 1, 1]]
oldtable[[-1, All, -1]] == oldtable[[1, All, 1]]
oldtable[[-1, -1, All]] == oldtable[[1, 1, All]]
oldtable[[-1, -1, -1]] == oldtable[[1, 1, 1]]

Out[69]= True

Out[70]= True

Out[71]= True

Out[72]= True

What Leonid does is illustrated in the figures below. Lines 4-6 of his code do something as illustrated in the left-hand panel: copying a line (darker color) of a plane already copied (light colors). Line 7 is illustrated by the right-hand panel. This is a cell to cell copy of diagonally opposing positions, and its operation is not included in any of the first three copy actions separately, but is a result of their successive operation.

enter image description here

like image 160
Sjoerd C. de Vries Avatar answered Dec 13 '22 16:12

Sjoerd C. de Vries


You can get it much faster and more memory-efficiently by modifying the original table as follows:

oldtable[[All, All, -1]] = oldtable[[All, All, 1]];
oldtable[[All, -1, All]] = oldtable[[All, 1, All]];
oldtable[[-1, All, All]] = oldtable[[1, All, All]];
oldtable[[All, -1, -1]] = oldtable[[All, 1, 1]];
oldtable[[-1, All, -1]] = oldtable[[1, All, 1]];
oldtable[[-1, -1, All]] = oldtable[[1, 1, All]];
oldtable[[-1, -1, -1]] = oldtable[[1, 1, 1]];

These assignments replace nested loops and are much faster, plus you don't need memory to store the copy. This is based on an extended vectorized functionality of the Part command (array and general expression indexing), particularly vectorized assignments. It also matters to have your numerical array in a Packed array form, which is often the case.

like image 42
Leonid Shifrin Avatar answered Dec 13 '22 14:12

Leonid Shifrin