0

我想计算一个二维数组“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}

谢谢！

4

1 回答 1

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}}, 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"]];

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

f2Cb = 编译[{{b2, _Real, 1}, {z, _Real, 1}, {b1, _Real, 1}}, With[{zeros = {ConstantArray[0., Length[z]]}}, 加入[零，表[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", CompilationOptions -> {"InlineExternalDefinitions" -> True}];

无论使用哪种方法，我都会在大约 0.7 秒内得到更实际示例的结果，而我的机器将花费 100 多秒，而这些全局变量在 Compile.

A more general approach might be to pass them as parameters (if they were likely to change rather than be constants). That would lead to a slightly slower run time though.

Regarding that option approach, you might have a look at ref/CompilationOptions in the Cocumentation Center

--- end edit 2 ---

于 2012-05-30T13:45:31.483 回答