一种有效的数据结构或方法来管理随时间增长的绘图数据 [英] An efficient data structure or method to manage plotting data that grow with time

问题描述

我想问一下，我管理模拟绘图结果的以下方式是否有效地使用了 Mathematica，以及是否有更实用"的方式来做到这一点.(可能会使用 Sow、Reap 等).

I'd like to ask if the following way I manage plotting result of simulation is efficient use of Mathematica and if there is a more 'functional' way to do it. (may be using Sow, Reap and such).

这是一个基本问题.假设您想模拟一个物理过程，比如钟摆，并且想要绘制解决方案运行时的时间序列(即时间与角度)(或任何其他类型的结果).

The problem is basic one. Suppose you want to simulate a physical process, say a pendulum, and want to plot the time-series of the solution (i.e. time vs. angle) as it runs (or any other type of result).

为了能够显示图表，需要在运行时保留数据点.

To be able to show the plot, one needs to keep the data points as it runs.

下面是一个简单的例子，它绘制了解决方案，但只绘制了当前点，而不是完整的时间序列:

The following is a simple example, that plots the solution, but only the current point, and not the full time-series:

Manipulate[
 sol = First@NDSolve[{y''[t] + 0.1 y'[t] + Sin[y[t]] == 0, y[0] == Pi/4, y'[0] == 0}, 
                      y, {t, time, time + 1}];

 With[{angle = y /. sol},
  (
   ListPlot[{{time, angle[time]}}, AxesLabel -> {"time", "angle"}, 
            PlotRange -> {{0, max}, {-Pi, Pi}}]
   )
  ],

 {{time, 0, "run"}, 0, max, Dynamic@delT, ControlType -> Trigger},
 {{delT, 0.1, "delT"}, 0.1, 1, 0.1, Appearance -> "Labeled"},

 TrackedSymbols :> {time},
 Initialization :> (max = 10)
 ]

以上并不有趣，因为人们只看到一个点在移动，而不是完整的解决方案路径.

The above is not interesting, as one only sees a point moving, and not the full solution path.

我目前处理这个问题的方法是使用 Table[] 分配一个足够大的缓冲区，以容纳可以生成的最大可能的时间序列大小.

The way currently I handle this, is allocate, using Table[], a buffer large enough to hold the largest possible time-series size that can be generated.

问题是时间步长是可以变化的，时间步长越小，生成的数据就越多.

The issue is that the time-step can change, and the smaller it is, the more data will be generated.

但由于我知道可能的最小时间步长(在此示例中为 0.1 秒)，并且我知道运行的总时间(此处为 10 秒)，因此我知道要分配多少.

But since I know the smallest possible time-step (which is 0.1 seconds in this example), and I know the total time to run (which is 10 seconds here), then I know how much to allocate.

我还需要一个索引"来跟踪缓冲区.使用这种方法，这里有一个方法:

I also need an 'index' to keep track of the buffer. Using this method, here is a way to do it:

Manipulate[

 If[time == 0, index = 0];

 sol = First@NDSolve[{y''[t] + 0.1 y'[t] + Sin[y[t]] == 0, y[0] == Pi/4,y'[0] == 0}, 
                      y, {t, time, time + 1}];

 With[{angle = y /. sol},
  (
   index += 1;
   buffer[[index]] = {time, angle[time]};

   ListPlot[buffer[[1 ;; index]], Joined -> True, AxesLabel -> {"time", "angle"}, 
            PlotRange -> {{0, 10}, {-Pi, Pi}}]
   )
  ],

 {{time, 0, "run"}, 0, 10, Dynamic@delT, AnimationRate -> 1, ControlType -> Trigger},
 {{delT, 0.1, "delT"}, 0.1, 1, 0.1, Appearance -> "Labeled"},

 {{buffer, Table[{0, 0}, {(max + 1)*10}]}, None},
 {{index, 0}, None},

 TrackedSymbols :> {time},
 Initialization :> (max = 10)
 ]

作为参考，当我在 Matlab 中做类似上面的事情时，它有一个很好的绘图工具，称为保持".这样就可以绘制一个点，然后说等一下"，这意味着下一个绘图不会删除该绘图上已有的内容，而是添加它.

For reference, when I do something like the above in Matlab, it has a nice facility for plotting, called 'hold on'. So that one can plot a point, then say 'hold on' which means that the next plot will not erase what is already on the plot, but will add it.

我在 Mathematica 中没有找到类似的东西，即动态更新当前绘图.

I did not find something like this in Mathematica, i.e. update a current plot on the fly.

我也不想在缓冲区运行时使用 Append[] 和 AppendTo[] 来构建缓冲区，因为这会很慢且效率不高.

I also did not want to use Append[] and AppendTo[] to build the buffer as it runs, as that will be slow and not efficient.

我的问题:除了我正在做的事情之外，是否有更高效的 Mathematica 方式(可以更快、更优雅)来完成上述典型任务?

My question: Is there a more efficient, Mathematica way (which can be faster and more elegent) to do a typical task such as the above, other than what I am doing?

谢谢，

更新:

关于为什么不一次解决所有 ODE 的问题.是的，这是可能的，但是出于性能原因，它可以简化很多事情，将其分块进行.这是一个带有初始条件的 ode 示例:

On the question on why not solving the ODE all at once. Yes, it is possible, but it simplifies things alot to do it in pieces, also for performance reasons. Here is an example with ode with initial conditions:

Manipulate[
 If[time == 0, index = 0];
 sol = First@
   NDSolve[{y''[t] + 0.1 y'[t] + Sin[y[t]] == 0, y[0] == y0, 
     y'[0] == yder0}, y, {t, time, time + 1}];

 With[{angle = (y /. sol)[time]},
  (
   index += 1;
   buffer[[index]] = {time, angle};

   ListPlot[buffer[[1 ;; index]], Joined -> True, 
    AxesLabel -> {"time", "angle"}, 
    PlotRange -> {{0, 10}, {-Pi, Pi}}])],

 {{time, 0, "run"}, 0, 10, Dynamic@delT, AnimationRate -> 1, 
  ControlType -> Trigger}, {{delT, 0.1, "delT"}, 0.1, 1, 0.1, 
  Appearance -> "Labeled"},
 {{y0, Pi/4, "y(0)"}, -Pi, Pi, Pi/100, Appearance -> "Labeled"},
 {{yder0, 0, "y'(0)"}, -1, 1, .1, Appearance -> "Labeled"},

 {{buffer, Table[{0, 0}, {(max + 1)*10}]}, None},
 {{index, 0}, None},

 TrackedSymbols :> {time},
 Initialization :> (max = 10)
 ]

现在，在一个之前解决系统一次，然后他们需要注意IC是否发生变化.这可以做到，但需要额外的逻辑，我之前已经做过很多次了，但它确实使事情复杂化了一点.我在这里写了一个小笔记.

Now, in one were to solve the system once before, then they need to watch out if the IC changes. This can be done, but need extra logic and I have done this before many times, but it does complicate things a bit. I wrote a small note on this here.

此外，我注意到随着时间的推移，通过在更小的时间段内解决系统问题，我可以获得更好的速度，而不是一次性解决整个问题.NDSolve 调用开销非常小.但是当 NDsolve for 的持续时间很长时，当人们要求 NDSolve 提供更高的准确度时，就会出现问题，如选项 AccuracyGoal ->, PrecisionGoal ->，我在时间间隔时无法做到这一点非常大.

Also, I noticed that I can get much better speed by solving the system for smaller time segments as time marches on, than the whole thing at once. NDSolve call overhead is very small. But when the time duration to NDsolve for is large, problems can result when one ask for higher accuracy from NDSolve, as in options AccuracyGoal ->, PrecisionGoal ->, which I could not when time interval is very large.

总的来说，与它在简化逻辑和速度方面的优势相比，为较小的段调用 NDSolve 的开销似乎要少得多(可能更准确，但我没有对此进行更多检查).我知道继续调用 NDSolve 似乎有点奇怪，但是在尝试了这两种方法(同时进行，但添加了检查其他控制变量的逻辑)与此方法之后，我现在倾向于使用此方法.

Overall, the overhead of calling NDSolve for smaller segments seems to much less compare to the advantages it makes in simplifing the logic, and speed (may be more accurate, but I have not checked on this more). I know it seems a bit strange to keep calling NDSolve, but after trying both methods (all at once, but add logic to check for other control variables) vs. this method, I am now leaning towards this one.

更新 2

我针对 2 个测试用例比较了以下 4 种方法:

I compared the following 4 methods for 2 test cases:

tangle[j][j] 方法(贝利撒留)AppendTo(由 Sjoerd 建议)动态链表 (Leonid)(有和没有 SetAttributes[linkedList, HoldAllComplete])预分配缓冲区(Nasser)

tangle[j][j] method (Belisarius) AppendTo (suggested by Sjoerd) Dynamic linked list (Leonid) (with and without SetAttributes[linkedList, HoldAllComplete]) preallocate buffer (Nasser)

我这样做的方法是将它运行在 2 个案例中，一个是 10,000 点，第二个是 20,000 点.我确实将 Plot[[] 命令留在那里，但没有在屏幕上显示它，这是为了消除实际渲染的任何开销.

The way I did this, is by running it over 2 cases, one for 10,000 points, and the second for 20,000 points. I did leave the Plot[[] command there, but do not display it on the screen, this is to eliminate any overhead of the actual rendering.

我在 Do 循环周围使用了 Timing[]，该循环迭代调用 NDSolve 的核心逻辑，并使用上述 delT 增量迭代时间跨度.没有使用 Manipulate.

I used Timing[] around a Do loop which iterate over the core logic which called NDSolve and iterate over the time span using delT increments as above. No Manipulate was used.

我在每次运行前都使用了 Quit[].

I used Quit[] before each run.

对于 Leonid 方法，我更改了他在 Do 循环中拥有的 Column[].我最后验证了，但是使用他的 getData[] 方法绘制数据，结果没问题.

For Leonid method, I changed the Column[] he had by the Do loop. I verified at the end, but plotting the data using his getData[] method, that the result is ok.

我使用的所有代码都在下面.我制作了一个表格，显示了 10,000 点和 20,000 点的结果.计时是每秒:

All the code I used is below. I made a table which shows the results for the 10,000 points and 20,000. Timing is per seconds:

 result = Grid[{
   {Text[Style["method", Bold]], 
    Text[Style["number of elements", Bold]], SpanFromLeft},
   {"", 10000, 20000},
   {"", SpanFromLeft},
   {"buffer", 129, 571},
   {"AppendTo", 128, 574},
   {"tangle[j][j]", 612, 2459},
   {"linkedList with SetAttribute", 25, 81},
   {"linkedList w/o SetAttribute", 27, 90}}
  ]

显然，除非我做错了什么，但下面的代码供任何人验证，Leonid 方法在这里很容易获胜.我也很惊讶 AppendTo 与预先分配数据的缓冲区方法一样好.

Clearly, unless I did something wrong, but code is below for anyone to verify, Leonid method wins easily here. I was also surprised that AppendTo did just as well as the buffer method which pre-allocated data.

这是我用来生成上述结果的稍微修改过的代码.

Here are the slightly modified code I used to generate the above results.

缓冲方法

delT = 0.01; max = 100; index = 0;
buffer = Table[{0, 0}, {(max + 1)*1/delT}];

Timing[
 Do[
  sol = First@
    NDSolve[{y''[t] + 0.1 y'[t] + Sin[y[t]] == 0, y[0] == Pi/4, 
      y'[0] == 0}, y, {t, time, time + 1}];

  With[{angle = y /. sol},
   (index += 1;
    buffer[[index]] = {time, angle[time]};
    foo = 
     ListPlot[buffer[[1 ;; index]], Joined -> True, 
      AxesLabel -> {"time", "angle"}, 
      PlotRange -> {{0, 10}, {-Pi, Pi}}]
    )
   ], {time, 0, max, delT}
  ]
 ]

AppendTo 方法

AppendTo method

Clear[y, t];
delT = 0.01; max = 200;
buffer = {{0, 0}};  (*just a hack to get ball rolling, would not do this in real code*)

Timing[
 Do[
  sol = First@
    NDSolve[{y''[t] + 0.1 y'[t] + Sin[y[t]] == 0, y[0] == Pi/4, 
      y'[0] == 0}, y, {t, time, time + 1}];

  With[{angle = y /. sol},
   (AppendTo[buffer, {time, angle[time]}];
    foo = 
     ListPlot[buffer, Joined -> True, AxesLabel -> {"time", "angle"}, 
      PlotRange -> {{0, 10}, {-Pi, Pi}}]
    )
   ], {time, 0, max, delT}
  ]
 ]

缠结[j][j]方法

Clear[y, t];
delT = 0.01; max = 200;
Timing[
 Do[
  sol = First@
    NDSolve[{y''[t] + 0.1 y'[t] + Sin[y[t]] == 0, y[0] == Pi/4, 
      y'[0] == 0}, y, {t, time, time + 1}];
  tangle[time] = y /. sol;
  foo = ListPlot[
    Table[{j, tangle[j][j]}, {j, .1, max, delT}],
    AxesLabel -> {"time", "angle"},
    PlotRange -> {{0, max}, {-Pi, Pi}}
    ]
  , {time, 0, max, delT}
  ]
 ]

动态链表方法

Timing[
 max = 200;

 ClearAll[linkedList, toLinkedList, fromLinkedList, addToList, pop, 
  emptyList];
 SetAttributes[linkedList, HoldAllComplete];
 toLinkedList[data_List] := Fold[linkedList, linkedList[], data];
 fromLinkedList[ll_linkedList] := 
  List @@ Flatten[ll, Infinity, linkedList];
 addToList[ll_, value_] := linkedList[ll, value];
 pop[ll_] := Last@ll;
 emptyList[] := linkedList[];

 Clear[getData];

 Module[{ll = emptyList[], time = 0, restart, plot, y},

  getData[] := fromLinkedList[ll];

  plot[] := Graphics[
    {
     Hue[0.67`, 0.6`, 0.6`],
     Line[fromLinkedList[ll]]
     },
    AspectRatio -> 1/GoldenRatio,
    Axes -> True,
    AxesLabel -> {"time", "angle"},
    PlotRange -> {{0, 10}, {-Pi, Pi}},
    PlotRangeClipping -> True
    ];

  DynamicModule[{sol, angle, llaux, delT = 0.01},

   restart[] := (time = 0; llaux = emptyList[]);
   llaux = ll;

   sol := 
    First@NDSolve[{y''[t] + 0.1 y'[t] + Sin[y[t]] == 0, y[0] == Pi/4, 
       y'[0] == 0}, y, {t, time, time + 1}];
   angle := y /. sol;

   ll := With[{res = 
       If[llaux === emptyList[] || pop[llaux][[1]] != time,
        addToList[llaux, {time, angle[time]}],
        (*else*)llaux]
      },
     llaux = res
     ];

   Do[
    time += delT;
    plot[]
    , {i, 0, max, delT}
    ]
   ]
  ]
 ]

感谢大家的帮助.

推荐答案

我不知道如何使用 Manipulate 获得您想要的东西，但我似乎已经设法通过自定义获得了一些接近的东西动态.以下代码将:使用链表以合理高效，使用按钮停止/恢复您的绘图，并在任何给定时间按需提供迄今为止收集的数据:

I don't know how to get what you want with Manipulate, but I seem to have managed getting something close with a custom Dynamic. The following code will: use linked lists to be reasonably efficient, stop / resume your plot with a button, and have the data collected so far available on demand at any given time:

ClearAll[linkedList, toLinkedList, fromLinkedList, addToList, pop, emptyList];
SetAttributes[linkedList, HoldAllComplete]; 
toLinkedList[data_List] := Fold[linkedList, linkedList[], data];
fromLinkedList[ll_linkedList] := List @@ Flatten[ll, Infinity, linkedList];
addToList[ll_, value_] := linkedList[ll, value];
pop[ll_] := Last@ll;
emptyList[] := linkedList[];


Clear[getData];
Module[{ll = emptyList[], time = 0, restart, plot, y},
   getData[] := fromLinkedList[ll];
   plot[] := 
      Graphics[{Hue[0.67`, 0.6`, 0.6`], Line[fromLinkedList[ll]]}, 
        AspectRatio -> 1/GoldenRatio, Axes -> True, 
        AxesLabel -> {"time", "angle"}, PlotRange -> {{0, 10}, {-Pi, Pi}}, 
        PlotRangeClipping -> True];
   DynamicModule[{sol, angle, llaux, delT = 0.1},
     restart[] := (time = 0; llaux = emptyList[]);
     llaux = ll;
     sol := First@
        NDSolve[{y''[t] + 0.1 y'[t] + Sin[y[t]] == 0, y[0] == Pi/4, y'[0] == 0}, 
             y, {t, time, time + 1}];
     angle := y /. sol;
     ll := With[{res =
              If[llaux === emptyList[] || pop[llaux][[1]] != time, 
                 addToList[llaux, {time, angle[time]}],
                 (* else  *)
                 llaux]},
              llaux = res];
     Column[{
        Row[{Dynamic@delT, Slider[Dynamic[delT], {0.1, 1., 0.1}]}],
        Dynamic[time, {None, Automatic, None}],
        Row[{
          Trigger[Dynamic[time], {0, 10, Dynamic@delT}, 
               AppearanceElements -> { "PlayPauseButton"}], 
          Button[Style["Restart", Small], restart[]]
        }],
        Dynamic[plot[]]
      }, Frame -> True]
  ]
]

此处的链接列表取代了您的 buffer，您无需预先分配并提前知道您将拥有多少数据点.plot[] 是一个自定义的低级绘图函数，尽管我们可能也可以使用 ListPlot.您可以使用播放"按钮停止和恢复绘图，并使用自定义的重新启动"按钮重置参数.

Linked lists here replace your buffer and you don't need to pre-allocate and to know in advance how many data points you will have. The plot[] is a custom low-level plotting function, although we probably could just as well use ListPlot. You use the "Play" button to both stop and resume plotting, and you use the custom "Restart" button to reset the parameters.

您可以在任何给定时间调用 getData[] 以获取迄今为止累积的数据列表，如下所示:

You can call getData[] at any given time to get a list of data accumulated so far, like so:

In[218]:= getData[]
Out[218]= {{0,0.785398},{0.2,0.771383},{0.3,0.754062},{0.4,0.730105},{0.5,0.699755},
{0.6,0.663304},{0.7,0.621093},{0.8,0.573517},{0.9,0.521021},{1.,0.464099},
{1.1,0.403294},{1.2,0.339193},{1.3,0.272424}}

这篇关于一种有效的数据结构或方法来管理随时间增长的绘图数据的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！