Optimizing interpolation in Mathematica

As part of my work, I often have to visualize complex 3 dimensional densities. One program suite that I work with outputs the radial component of the densities as a set of 781 points on a logarithmic grid, ri = (Rmax/Rstep)^((i-1)/(pts-1), times a spherical harmonic. For low symmetry systems, the number of spherical harmonics can be fairly large to ensure accuracy, e.g. one system requires 49 harmonics corresponding to lmax = 6. So, to use this data within Mathematica, I would have a sum of up to 49 interpolated functions with each multiplied by a different spherical harmonic. While using v.6 and constructing the interpolated radial functions using Interpolation and 开发者_开发知识库setting r = Sqrt(x^2 + y^2 + z^2), I would stop ContourPlot3D after well over an hour without anything displayed. This included reducing both the InterpolationOrder and MaxRecursion to 1.

Several alternatives presented themselves:

1. Evaluate the density function on a fixed grid, and use ListContourPlot instead.
2. Or, linearly spline the radial function and use Piecewise to stitch them together. (This presented itself, as I could use simplify to help reduce the complexity of the resulting function.)

I ended up using both, as InterpolatingFunction gives a noticeable delay in its evaluation, and with up to 49 interpolated functions to evaluate, any delay can become noticeable. Also, ContourPlot3D was faster with the spline, but it didn't give me the speed up I desired.

I'll freely admit that I haven't tried Interpolation on v.7, nor I have tried this on my upgraded hardware (G4 v. Intel Core i5). However, I'm looking for alternatives to my current scheme; preferably, one where I can use ContourPlot3D directly. I could try some other form of spline, such as a B-spline, and possibly combine that with UnitBox instead of using Piecewise.

Edit: Just to clarify, my current implementation involves creating a first order spline for each radial part, multiplying each one by their respective spherical harmonic, summing and Simplifying the equations on each radial interval, and then using Piecewise to bind them into one function. So, my implementation is semi-analytical in that the spherical harmonics are exact, and only the radial part is numerical. This is part of the reason why I would like to be able to use ContourPlot3D, so that I can take advantage of the semi-analytical nature of the data. As a point of note, the radial grid is fine enough that a good representation of the radial part is generated and can be smoothly interpolated. While this gave me a significant speed-up, when I wrote the code, it was still to slow for the hardware I was using at the time.

So, instead of using ContourPlot3D, I would first generate the function, as above, then I would evaluate it on an 80³ Cartesian grid. It is the data from this step that I used in ListContourPlot3D. Since this is not an adaptive grid, in some places this was too course, and I was missing features.

If you can do without Mathematica, I would suggest you have a look at Paraview (US government funded FOSS, all platforms) which I have found to be superior to everything when it comes to visualizing massive amounts of data. The core of the software is the "Visualization Toolkit" VTK, and you can find/write other frontends if need be.

VTK/Paraview can handle almost any data-type: scalar and vector on structured grids or random points, polygons, time-series data, etc. From Mathematica I often just dump grid data into VTK legacy format which in then simplest case looks like this

# vtk DataFile Version 2.0
Generated by mma via vtkGridDump

ASCII

DATASET STRUCTURED_POINTS
DIMENSIONS 49 25 15
SPACING 0.125 0.125 0.0625
ORIGIN 8.5 5. 0.7124999999999999

POINT_DATA 18375
SCALARS  RF_pondpot_1V1MHz1amu  double 1
LOOKUP_TABLE default

0.04709501616121583
0.04135197485227461
... <18373 more numbers> ...

HTH!

If it really is the interpolation of the radial functions that is slowing you down, you could consider hand-coding that part based on your knowledge of the sample points. As demonstrated below, this gives a significant speedup:

I set things up with your notation. lookuprvals is a list of 100000 r values to look up for timing.

First, look at stock interpolation as a basemark

With[{interp=Interpolation[N@Transpose@{rvals,yvals}]},
  Timing[interp[lookuprvals]][[1]]]
Out[259]= 2.28466

Switching to 0th-order interpolation is already an order of magnitude faster (first order is almost same speed):

With[{interp=Interpolation[N@Transpose@{rvals,yvals},InterpolationOrder->0]},
  Timing[interp[lookuprvals]][[1]]]
Out[271]= 0.146486

We can get another 1.5 order of magnitude by calculating indices directly:

Module[{avg=MovingAverage[yvals,2],idxfact=N[(pts-1) /Log[Rmax/Rstep]]},
  Timing[res=Part[avg,Ceiling[idxfact Log[lookuprvals]]]][[1]]]
Out[272]= 0.006067

As a middle ground, do a log-linear interpolation by hand. This is slower than the above solution but still much faster than stock interpolation:

Module[{diffs=Differences[yvals],
  idxfact=N[(pts-1) /Log[Rmax/Rstep]]},
  Timing[Block[{idxraw,idxfloor,idxrel},
    idxraw=1+idxfact Log[lookuprvals];
    idxfloor=Floor[idxraw];
    idxrel=idxraw-idxfloor;
    res=Part[yvals,idxfloor]+Part[diffs,idxfloor]idxrel  
  ]][[1]]]
Out[276]= 0.026557

If you have the memory for it, I would cache the spherical harmonics and radius (or even radius-index) on the full grid. Then flatten the grid caches so you can do

 Sum[ interpolate[yvals[lm],gridrvals] gridylmvals[lm], {lm,lmvals} ]

and recreate your grid as discussed here.