开发者

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]; 

myinterpo开发者_运维知识库latingfunction = 
 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.


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.

Mathematica large table periodic interpolation


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.


Just because I'm a bit of a loon, Leonid's solution can be written as:

a = {0, 1}~Tuples~3~SortBy~Tr // Rest;

MapThread[
  (oldtable[[Sequence @@ #]] = oldtable[[Sequence @@ #2]]) &,
  {-a, a} /. 0 -> All
];


Thanks for all the answers. I tried the suggestion by leonid but when I print my oldtable, it was still (dimcub -1)^3 dimensional. New elements were defined and I can see them individually but they do not show up as part of the oldtable when I print the whole table. So I ended up with something similar which is doing exactly what I needed:

oldtable= PadRight[oldtable, {dimcub, dimcub, dimcub}];
oldtable[[All, All, dimcub]] = oldtable[[All, All, 1]];
oldtable[[All, dimcub, All]] = oldtable[[All, 1, All]];
oldtable[[dimcub, All, All]] = oldtable[[1, All, All]];
oldtable[[All, dimcub, dimcub]] = oldtable[[All, 1, 1]];
oldtable[[dimcub, All, dimcub]] = oldtable[[1, All, 1]];
oldtable[[dimcub, dimcub, All]] = oldtable[[1, 1, All]];
oldtable[[dimcub, dimcub, dimcub]] = oldtable[[1, 1, 1]];

The answer by Wizard is far too advanced for my level of mathematica..

0

上一篇:

下一篇:

精彩评论

暂无评论...
验证码 换一张
取 消

最新问答

问答排行榜