Оптимизация интерполяции в Mathematica

Question

В рамках своей работы мне часто приходится визуализировать сложные трехмерные плотности. Один набор программ, с которым я работаю, выводит радиальную составляющую плотностей в виде набора из 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.

Несколько альтернатив представили:

1. Оцените функцию плотности на фиксированной сетке и используйте вместо нее ListContourPlot.
2. Или, линейно сплайнируйте радиальную функцию и используйте Piecewise, чтобы соединить их вместе. (Это представилось само собой, так как я мог использовать упрощение, чтобы уменьшить сложность получаемой функции.)

В итоге я использовал оба, так как InterpolatingFunction дает заметную задержку в своей оценке, и с помощью до 49 интерполированных функций для оценки любая задержка может стать заметной. Кроме того, ContourPlot3D был быстрее со сплайном, но это не давало мне желаемой скорости.

Я свободно признаю, что я не пробовал Interpolation на v.7, и я не пробовал это на моем обновленном оборудовании (G4 v. Intel Core i5). Тем не менее, я ищу альтернативы моей нынешней схеме; желательно тот, где я могу использовать ContourPlot3D напрямую. Я мог бы попробовать другую форму сплайна, такую ​​как B-сплайн , и, возможно, комбинировать ее с UnitBox вместо использования Piecewise.

Редактировать: Просто чтобы прояснить, моя текущая реализация включает в себя создание сплайна первого порядка для каждой радиальной части, умножение каждого из них на их соответствующие сферические гармоники, суммирование и Simplify уравнений для каждого радиального интервала и затем с помощью Piecewise связать их в одну функцию. Итак, моя реализация является полуаналитической в ​​том смысле, что сферические гармоники являются точными, и только радиальная часть является численной. Это одна из причин, почему я хотел бы иметь возможность использовать ContourPlot3D, чтобы я мог воспользоваться полуаналитическим характером данных. Следует отметить, что радиальная сетка достаточно тонкая, чтобы получить хорошее представление радиальной части, и ее можно плавно интерполировать. Хотя это дало мне значительное ускорение, когда я писал код, он все еще замедлялся для оборудования, которое я использовал в то время.

Таким образом, вместо использования ContourPlot3D, я сначала сгенерирую функцию, как указано выше, затем я оценил бы ее на декартовой сетке 80 ³ . Это данные этого шага, которые я использовал в ListContourPlot3D. Поскольку это не адаптивная сетка, в некоторых местах это было слишком естественно, и мне не хватало возможностей.

Answer 1

Если вы можете обойтись без Mathematica, я бы посоветовал вам взглянуть на Paraview (FOSS, все платформы, финансируемые правительством США), который, как я обнаружил, превосходит все, когда дело доходит до визуализации огромных количествданных.Ядром программного обеспечения является «Visualization Toolkit» VTK , и вы можете найти / написать другие интерфейсы, если это необходимо.

VTK / Paraview может обрабатывать практически любые типы данных: скалярные ивектор на структурированных сетках или случайных точках, полигонах, данных временных рядов и т. д. Из Mathematica я часто просто выгружаю данные сетки в VTK устаревший формат , который в простейшем случае выглядит следующим образом

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

НТН!

Answer 2

Если это действительно интерполяция радиальных функций, которая замедляет вас, вы можете рассмотреть ручное кодирование этой части, основываясь на вашем знании точек выборки.Как показано ниже, это дает значительное ускорение:

Я настроил все с помощью вашей записи.lookuprvals - это список значений 100000 r для поиска времени.

Сначала рассмотрим интерполяцию акций как базовую отметку

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

Переключение на интерполяцию 0-го порядка уже выполненона порядок быстрее (первый порядок почти такой же скорости):

With[{interp=Interpolation[N@Transpose@{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} ]

и воссоздайте свою сетку, как обсуждалось здесь .