0

他の 3 つのリスト (z、b1、b2) のテストに基づいて要素が計算される 2D 配列 "tocalc" を計算したいと考えています。

(*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};

これに対する1つの解決策は

(*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 倍高速です。2 つの Do ループを使用したソリューションが Mathematica に最適ではないことはわかっているので、MapIndexed、Table、関数、Conditional などを使用した他のソリューションをいくつか試しました。しかし、すべてが本当に速いわけではなく、2 つの 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;

1 つの解決策:

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"と、約2倍速くなります。

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

もちろん、グローバル変数がある場合 (つまり、Compile の外部で定義されている場合) は、もう少し作業が必要です。私は2つの可能性を認識しています。バージョン 8 より前では、次のように、Compile の周りで With[] を使用して定数を吸い込んでいました。

f2C = With[{n = n, deltax = deltax, A00 = A00}, 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* デルタ)])^(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 = Compile[{{b2, _Real, 1}, {z, _Real, 1}, {b1, _Real, 1}}, With[{zeros = {ConstantArray[0., 長さ[z]]}}, Join [zeros, Table[If[ b1[[i]] <= z[[k]]] <= b2[[i]], -(A00*(Abs[(b2[[i + 1]]] - b2[[ i - 1]])/(2.0* デルタ)])^(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, 長さ[b2] - 1}, {k, Length[z]}], zeros]], CompilationTarget -> "C", CompilationOptions -> {"InlineExternalDefinitions" -> True}];

どちらを使用しても、より現実的な例では約 0.7 秒で結果が得られますが、Compile 内でこれらのグローバルが定義されていない場合、私のマシンでは 100 秒以上かかります。

より一般的なアプローチは、それらをパラメーターとして渡すことです (定数ではなく変更される可能性が高い場合)。ただし、実行時間が少し遅くなります。

そのオプションのアプローチに関しては、Cocumentation Center の ref/CompilationOptions を参照してください。

--- 編集終了 2 ---

于 2012-05-30T13:45:31.483 に答える