在Mathematica中优化插值

作为我工作的一部分，我经常需要可视化复杂的三维密度。与我一起工作的一个程序套件将密度的径向分量输出为对数网格上的一组781个点，即球面谐波的数量，即ri = (Rmax/Rstep)^((i-1)/(pts-1)。对于低对称系统，球谐函数的数量可能相当大以确保精度，例如，一个系统需要49个对应于lmax = 6的谐波。因此，要在Mathematica中使用这些数据，我将得到总计49个内插函数，每个函数乘以不同的球谐函数。在使用v.6并使用Interpolation和设置r = Sqrt(x^2 + y^2 + z^2)构建插值径向函数时，我会在一个多小时后停止ContourPlot3D而不显示任何内容。这包括减少两者InterpolationOrder和MaxRecursion为1在Mathematica中优化插值

几种替代品呈现自己：

评估固定网格密度函数，并使用ListContourPlot代替。
或者，线性旋转径向函数并使用Piecewise将它们缝合在一起。（此提出了自己，因为我可以用简化，以帮助减少所产生的功能的复杂性。）

我结束了使用这两种，作为InterpolatingFunction给出的评价明显的延迟，以及高达49个插补功能评估，任何延迟都会变得明显。此外，ContourPlot3D与样条更快，但它并没有给我加速我想要的。

我会毫不犹豫地承认，我还没有尝试012.v.7，我也没有尝试过我的升级硬件（G4 v。英特尔酷睿i5）。但是，我正在寻找替代我目前的方案;最好是我可以直接使用ContourPlot3D的地方。我可以尝试其他形式的样条，如B-spline，并可能将其与UnitBox结合使用，而不是使用Piecewise。

编辑：只是为了澄清，我的当前的实现包括创建用于每个径向部分的一阶样条，乘以每一个由各自的球谐，总结和Simplify荷兰国际集团上每一径向间隔的公式，然后使用Piecewise将它们绑定到一个函数中。所以，我的实现是半分析的，因为球谐函数是精确的，只有径向部分是数值的。这是我希望能够使用ContourPlot3D的原因之一，以便我可以利用数据的半分析性质。作为一个注意点，径向网格足够精细，可以生成径向部分的良好表示，并且可以平滑插值。虽然这给我一个显着的提速，但是当我编写代码时，对于我当时使用的硬件来说，速度仍然很慢。

因此，我不会使用ContourPlot3D，而是首先生成如上所述的函数，然后我将在笛卡尔网格上对其进行评估。这是我在ListContourPlot3D中使用的这一步的数据。由于这不是一个自适应网格，在某些地方，这太过于自然，而且我缺少特征。

来源

2010-10-25 rcollyer

只是为了确保我正确理解你：每个径向函数上有781个数据点，总共有781 * 49个标量？ – Janus 2010-10-26 05:51:09

您是否在'ListContourPlot'中使用'DataRange'和数据 - 否则它看起来像有一个内插层，请参阅http://stackoverflow.com/questions/2497517/mathematica-listcontourplot3d – Janus 2010-10-26 07:26:46

@Janus，是的，你阅读正确。有高达781 * 49的标量。正如我上面指出的，数据的角度部分是分析性的，只有径向部分是数字的。另外，在使用'ListContourPlot3D'之前，我已经遇到了用于'{x，y，z，f}'形式的低分辨率错误，但是在这里我生成了一个3D数组并指定了'DataRange'。此外，虽然速度较慢，但“ListContourPlot3D”中的内插层仍比“ContourPlot3D”快。 – rcollyer 2010-10-26 13:21:13

如果它确实是使放慢速度的径向函数的插值，那么可以考虑根据您对采样点的了解手动编码该部分。如下所示，这给出了显着的加速：

我使用符号进行设置。 lookuprvals是查找时间的100000 r值的列表。

首先，看股票插值为basemark

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

切换到0级插补已经是一个数量级的速度更快（第一顺序是几乎相同的速度）：

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

我们通过直接计算指数可以得到另一个1.5的数量级：

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

作为中间地带，请手动进行对数线性插值。这比上述溶液慢，但还是比股票插值快得多：

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

如果你有这方面的记忆，我会缓存在满格的球谐和半径（甚至半径指数）。然后展平网格缓存，以便您可以执行

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

并按照here所述重新创建网格。

来源

2010-10-26 06:41:13 Janus

对不起，我花了这么长时间才看到这个。在你的第三个例子中，你计算移动平均值，但不要对它做任何事情。此外，计算指标和使用'部分'是聪明的，基本上是我想要做的。我已经编译了这些函数，但是也可以用调度表合理完成。我不确定第0顺序对于我想要做的事情有多好，但如果它有效，那会提供很大的加速。很多想法，因为这是目前的一个副项目，在我能够实施它之前，这将是一段时间。感谢您的想法。 – rcollyer 2010-10-29 14:55:19

将示例3更改为我的想法：返回每个间隔中端点值的平均值。再看一遍，在我看来，它可能实际上是值得对径向数据进行过采样（通过样条或其他），以便您在网格上进行采样时可以执行快速0阶插值？ – Janus 2010-10-30 14:41:41

如果没有Mathematica，我建议你看看Paraview（美国政府资助的FOSS，所有平台），我发现它在可视化海量数据方面优于所有产品。该软件的核心是“可视化工具包”VTK，如果需要，您可以找到/写入其他前端。

VTK/Paraview几乎可以处理任何数据类型：结构化网格或随机点，多边形，时间序列数据等上的标量和矢量。从Mathematica我经常只将网格数据转储到VTK legacy format，然后在最简单的情况下看起来像这样

# 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！

来源

2010-10-26 04:05:35 Janus

+1，我曾经遇到过这些，但那时学习曲线已经足够让我避开它们。他们绝对是工业实力的解决方案，我可能不得不最终诉诸于他们。但是，我希望得到一个基于Mathematica的解决方案，因为这是我最了解的，实施的时间最好是几个小时。 – rcollyer 2010-10-26 04:37:22

在Mathematica中优化插值

回答

相关问题