在 Mathematica 中的其他列表上生成具有基于元素的条件的列表答案

【问题标题】：Generating a list with elements based conditionals on other lists in Mathematica在 Mathematica 中的其他列表上生成具有基于元素的条件的列表
【发布时间】：2012-05-30 13:14:36
【问题描述】：

我想计算一个二维数组“tocalc”，其中的元素是根据对其他三个列表 (z,b1,b2) 的测试来计算的。

(*example data*)
z = Range[0, 100, 10];
x = Range[10];
b1 = ConstantArray[0., Length[x]];
tocalc = ConstantArray[0, {Length[x], Length[z]}];
b2 = {0, 20, 30, 40, 50, 40, 30, 20, 10, 0};

一个解决方案是

(*simple but slow solution*)
Do[
 Do[
   If[z[[k]] <= b2[[i]] && z[[k]] >= b1[[i]], 
    tocalc[[i, k]] = (b2[[i + 1]] - b2[[i - 1]])],
   {k, 1, Length[z]}];,
 {i, 2, Length[x] - 1}]

结果

{{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {30, 30, 30, 0, 0, 0, 0, 0, 0, 0, 
  0}, {20, 20, 20, 20, 0, 0, 0, 0, 0, 0, 0}, {20, 20, 20, 20, 20, 0, 
  0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
  0}, {-20, -20, -20, -20, -20, 0, 0, 0, 0, 0, 
  0}, {-20, -20, -20, -20, 0, 0, 0, 0, 0, 0, 0}, {-20, -20, -20, 0, 0,
   0, 0, 0, 0, 0, 0}, {-20, -20, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0,
   0, 0, 0, 0, 0, 0, 0, 0}}

问题：如何在 Mathematica 中有效地做到这一点？

如果评估 10000 次，则需要 3.66 秒。而在 Matlab 中，这需要 0.04 秒，所以 Matlab 快了近 100 倍。我知道使用两个 Do 循环的解决方案对于 Mathematica 来说并不完美，因此我尝试了其他几种解决方案，例如 MapIndexed、Table、函数、条件等。但是这两个 Do 循环并不是真的更快，甚至可能更慢。这是 MapIndexed 的一个示例：

tocalc = ConstantArray[0, {Length[x], Length[z]}];
MapIndexed[
  If[z[[Part[#2, 2]]] <= b2[[Part[#2, 1]]] && 
     z[[Part[#2, 2]]] >= b1[[Part[#2, 1]]] && Part[#2, 1] >= 2 && 
     Part[#2, 1] <= Length[x] - 1, 
    tocalc[[Part[#2, 1], Part[#2, 2]]] = (b2[[Part[#2, 1] + 1]] - 
       b2[[Part[#2, 1] - 1]]), 0.] &, tocalc, {2}];

理想的解决方案应该适用于更大的矩阵和实数以及更复杂的条件。

---编辑：

由于在我的实际问题中看起来一些解决方案甚至更慢，所以这里是它的一个例子：

现实世界的问题

b2 = {0.`, 0.`, 0.`, 990.3440201085594`, 1525.7589030785484`, 
   1897.6531659202747`, 2191.6073263357594`, 2433.0441988616717`, 
   2630.6658409463894`, 2799.347578394955`, 2944.656306810331`, 
   3070.718467691769`, 3179.485627984329`, 3272.3788096129415`, 
   3346.199103579602`, 3405.384848015466`, 3346.199103579602`, 
   3272.3788096129415`, 3179.485627984329`, 3070.718467691769`, 
   2944.656306810331`, 2799.347578394955`, 2630.6658409463894`, 
   2433.0441988616717`, 2191.6073263357594`, 1897.6531659202747`, 
   1525.7589030785484`, 990.3440201085594`, 0.`, 0.`, 0.`};
z = {0.`, 250.`, 500.`, 750.`, 1000.`, 1250.`, 1500.`, 1750.`, 2000.`,
   2250.`, 2500.`, 2750.`, 3000.`, 3250.`, 
  3500.`}; (*z(k)*)
imax = 31; (*number of x(i)*)
b1 = ConstantArray[0., imax]; (*lower boundary, can be different form 0*)
deltax = 50000.`;
mmax = 10000.; (*number of calculations*)
A00 = 1.127190283243198`*^-12; (*somefactor*)
n = 3;

一个解决方案：

f2C = Compile[{{b2, _Real, 1}, {z, _Real, 1}, {b1, _Real, 1}},
   With[{zeros = {ConstantArray[0., Length[z]]}},
    Join[zeros, 
     Table[If[
       b1[[i]] <= z[[k]] <= 
        b2[[i]], -(A00*(Abs[(b2[[i + 1]] - b2[[i - 1]])/(2.0*
                 deltax)])^(n - 
              1.0)*(b2[[i]]^(n + 1.) - (b2[[i]] - z[[k]])^(n + 
                1.)))*((b2[[i + 1]] - b2[[i - 1]])/(2.0*deltax))
       , 0.],
      {i, 2, Length[b2] - 1}, {k, Length[z]}
      ], zeros]]
   , CompilationTarget -> "C"];

结果是

Timing[Do[f2C[b2, z, b1];, {mmax}]]
Out[85]= {81.3544, Null}

谢谢！

【问题讨论】：

当你说“如果这被评估了 10000 次，那么它将花费超过 30 秒。”，你的意思是把你的前两个代码块包装在 Do[ ... , {10000}] // Timing 中吗？这在 6.5 年前的移动 CPU 上运行 7 秒。或者你的意思是100000 次？只是想确保我的基准测试是正确的。
好的。谢谢你是对的，我实际上使用了 100000 次。我编辑了文本。
在 MATLAB 上需要多长时间？就这样我知道什么时候停止尝试优化它:) “矢量化”它将时间从 7 秒推到 0.5 到 1 秒之间（即大约 10 倍加速），而无需编译。但是代码并不漂亮且不直观。对于更大的矩阵，加速会更大。
我现在在 Matlab 和 10^4 评估中做了同样的例子。耗时 0.04 秒，在 Mathematica 中为 3.66。所以 Matlab 快了大约 91.5 倍。（在另一台机器上，但 Mathematica 甚至更新）。

标签： performance wolfram-mathematica

【解决方案1】：

您可以执行以下操作。不过，您需要弄清楚如何处理边界（未定义 b2[[i+1]] 或 b2[[i-1]]）。

f[x_, y_] := If[x[[1]] <= y <= x[[2]], x[[4]] - x[[3]], 0]

在这里我限制了 Outer 的水平，这样我就不需要改变头部（就像我在原始响应中所做的那样）。

In[1309]:= Outer[f, 
 Transpose[{b1, b2, RotateRight[b2], RotateLeft[b2]}], z, 1]

Out[1309]= {{20, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {30, 30, 30, 0, 0, 0, 
  0, 0, 0, 0, 0}, {20, 20, 20, 20, 0, 0, 0, 0, 0, 0, 0}, {20, 20, 20, 
  20, 20, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
  0}, {-20, -20, -20, -20, -20, 0, 0, 0, 0, 0, 
  0}, {-20, -20, -20, -20, 0, 0, 0, 0, 0, 0, 0}, {-20, -20, -20, 0, 0,
   0, 0, 0, 0, 0, 0}, {-20, -20, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {-10, 0, 
  0, 0, 0, 0, 0, 0, 0, 0, 0}}

速度检查：

In[1298]:= Timing[
 Do[Outer[f, 
   Apply[list, 
    Transpose[{b1, b2, RotateRight[b2], RotateLeft[b2]}], {1}], 
   z], {10^4}]]

Out[1298]= {2.68, Null}

我们可以编译函数以获得更好的速度。

fC = Compile[{{x, _Integer, 1}, {y, _Integer}}, 
   If[x[[1]] <= y <= x[[2]], x[[4]] - x[[3]], 0]];

In[1306]:= Timing[
 Do[Outer[fC, Transpose[{b1, b2, RotateRight[b2], RotateLeft[b2]}], z,
    1], {10^4}]]

Out[1306]= {0.8, Null}

--- 编辑---

变体包括编译整个例程。这是一个这样的。

ff = Compile[{{b1, _Integer, 1}, {b2, _Integer, 1}, {z, _Integer, 
     1}},
   With[{lc = 
      RotateRight[ListConvolve[{1, 0, -1}, b2, {-1, -1}, 0]]},
    Table[
     If[b1[[i]] <= z[[k]] <= b2[[i]], lc[[i]], 0], {i, 
      Length[b2]}, {k, Length[z]}
     ]]];
In[385]:= Timing[Do[ff[b1, b2, z], {10^4}]]

Out[385]= {0.24, Null}

如果我添加CompilationTarget -> "C"，那么它的速度会快两倍。

另一个变体，在 C 代码中，不到 0.1 秒。

In[441]:= 
ff2C = Compile[{{b1, _Integer, 1}, {b2, _Integer, 1}, {z, _Integer, 
     1}},
   With[{zeros = {ConstantArray[0, Length[z]]}},
    Join[zeros, Table[
      If[b1[[i]] <= z[[k]] <= b2[[i]], b2[[i + 1]] - b2[[i - 1]], 
       0], {i, 2, Length[b2] - 1}, {k, Length[z]}
      ], zeros]], CompilationTarget -> "C"];

In[442]:= Timing[Do[ff2C[b1, b2, z], {10^4}]]

Out[442]= {0.04, Null}

In[443]:= ff2C[b1, b2, z]

Out[443]= {{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {30, 30, 30, 0, 0, 0, 0,
   0, 0, 0, 0}, {20, 20, 20, 20, 0, 0, 0, 0, 0, 0, 0}, {20, 20, 20, 
  20, 20, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
  0}, {-20, -20, -20, -20, -20, 0, 0, 0, 0, 0, 
  0}, {-20, -20, -20, -20, 0, 0, 0, 0, 0, 0, 0}, {-20, -20, -20, 0, 0,
   0, 0, 0, 0, 0, 0}, {-20, -20, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0,
   0, 0, 0, 0, 0, 0, 0, 0}}

我猜有些变种可能会更快。

---结束编辑---

--- 编辑 2 ---

当然，如果您有全局变量（即在编译之外定义），那么还有更多工作要做。我知道两种可能性。在版本 8 之前，人们会在 Compile 周围使用 With[] 来吸收常量，如下所示。

f2C = With[{n = n, deltax = deltax, A00 = A00}, 编译[{{b2, _Real, 1}, {z, _Real, 1}, {b1, _Real, 1}}, 用[{zeros = {ConstantArray[0., Length[z]]}}，加入[零，表[如果[ b1[[i]] (b2[[i]]^(n + 1.) - (b2[[i]] - z[[k]])^(n + 1.)))((b2[[i + 1]] - b2[[i - 1]])/(2.0*deltax)), 0.], {i, 2, Length[b2] - 1}, {k, Length[z]}], zeros]], 编译目标 -> "C"]];

在版本 8 中，以下实现了相同的效果。

f2Cb = 编译[{{b2, _Real, 1}, {z, _Real, 1}, {b1, _Real, 1}}, 用[{zeros = {ConstantArray[0., Length[z]]}}，加入[零，表[如果[ b1[[i]] (b2[[i]]^(n + 1.) - (b2[[i]] - z[[k]])^(n + 1.)))((b2[[i + 1]] - b2[[i - 1]])/(2.0*deltax)), 0.], {i, 2, Length[b2] - 1}, {k, Length[z]}], zeros]], 编译目标 -> "C", CompilationOptions -> {"InlineExternalDefinitions" -> True}];

无论哪种情况，我都会在大约 0.7 秒内得到更实际示例的结果，而我的机器如果没有在 Compile 中定义这些全局变量，则需要 100 多秒。

更通用的方法可能是将它们作为参数传递（如果它们可能会更改而不是常量）。不过，这会导致运行时间稍慢。

关于该选项方法，您可以查看 Cocumentation Center 中的 ref/CompilationOptions

--- 结束编辑 2 ---

【讨论】：

好的，谢谢。这也适用于现实世界的问题，并且需要大约一半的时间作为我的其他解决方案。但它仍然比 Matlab 中的相同内容慢很多。所以我担心它永远不会像在 Matlab 中那样快。如果没有这个计算，10000 次评估需要 1.6 秒，而使用它需要 49 秒。
好的，再次感谢，我试过了，它更快，但仍然不够快。作为参考，现实世界的问题仍然比 Matlab 等效问题慢 10 倍以上。也许这里可以使用 ListCorrelate 或 ListConvolute？
我在现实世界的问题上尝试了您的新解决方案 (In[441])，但 10^4 次评估仍需要 81 秒。也许是因为我使用 Reals。我是新来的，我应该发布现实世界的问题吗？（约 1400 个字符）
我认为这不到 100 行？如果是这样，那可以在这里发帖。也就是说，我不确定我能否解决它。我可以说一件事：假设它们是机器浮点数，使用实数而不是整数不太可能成为瓶颈。
请查看我对问题的编辑。我添加了现实世界的问题。