標籤：

Mathematica中已知位置列表,如何直接對目標元素進行操作？

02-06

假設我已經通過 Position 函數得到了一個位置列表，現在想利用此列表對一個已知列表的對應元素直接乘上一個常數，該如何操作？

例：

cn = 2;

vc = Table[

Cos[([Pi] (n1 + 0.5) k1)/cn] Cos[([Pi] (n2 + 0.5) k2)/cn], {n1, 0,

cn - 1}, {n2, 0, cn - 1},

{k1, 0, cn - 1}, {k2, 0, cn - 1}];

vp = Position[vc, _?Negative]

現在 vc 是一個2*2*2*2的高維數組，vp 是一個位置矩陣，其每一行表示 vc 中小於零的元素的位置指標。

下面，如果我想讓 vc 中這些小於零的元素都乘以10，該如何做呢？

用循環肯定是可以的，但可否利用位置矩陣 vp，直接完成操作呢，我查看了 Part 函數的 help，發現似乎沒有這種用法，而 ReplacePart 函數則只能進行替換操作，似乎不能完成元素運算……

希望能得到一些提示，謝謝！

Mathematica code:

Block[{n=40, t, vc},
t = Array[Cos[Pi/n (#+0.5) #2],{n,n},0];
vc = Outer[Times,t,t];
vc Clip[vc,{0,-1},{10,1}]
];//AbsoluteTiming

(* vc=Transpose[Outer[Times,t,t],{3,1,4,2}] *)
Clip[x, {0, -1}, {10, 1}]//Simplify

總結了一下樓上幾位高手的答案，另外對比了每個方案的速度，列於下圖：

結論：

利用模式匹配，可以避免循環，完成對數組中特定元素的操作，方法有多種，但是速度有所差別。

1. 在不用到特殊函數的場合，可以考慮使用Compile加速（有C編譯器時優先編譯為C，另外可加上並行運算指令）;

2. 不能編譯的話，優先使用Map函數；

3. 與Matlab一樣，多考慮向量化編程；

4. 盡量減少中間變數；避免引入副作用。（關於副作用的原理，感謝@章佳傑 知友的熱心解答！）

但是和Matlab相比，還是慢了不少，Matlab用的是最簡單的四層for循環實現；

只有Compile後的方法比Matlab快，不知原因何在？

原因是Matlab默認開啟了JIT加速，針對循環做了優化。（感謝@曹洪洋 知友的熱心解答！）

P.S.

Matlab代碼段附於下（在我的筆記本上跑，比Mathematica下的Map方法快9倍左右）：

clear all
N = 40;
tic
vc = zeros(N,N,N,N);
for k1 = 0:N-1
for k2 = 0:N-1
for m = 0:N-1
for n = 0:N-1
c1 = cos (pi/N * (m+1/2) * k1);
c2 = cos (pi/N * (n+1/2) * k2);
c = c1 * c2;
vc(k1+1,k2+1,m+1,n+1) = ((c&<0)*10 + (c&>0))*c;
end
end
end
end
vc;
TotalTime = toc

看到有人提到MATLAB，不過代碼寫的不理想。補充一下：

function vc = fun4d
N = 40; k = 0:N-1;
vc = bsxfun(@times, cos(pi/N*(k+0.5)*k),shiftdim(cos(pi/N*(k+0.5)*k),-2));
vc = vc + (vc&<0).*vc*9;

這個代碼在較新版本上應該不會比上述的其他語言低

如果是15b或16a，可以試用一下未來版本的特性，先執行（可以寫在MATLAB啟動文件中避免每次手動執行）：

builtin(_useSingletonExpansion,1);

之後bsxfun可以簡化為：

function vc = fun4d
N = 40; k = 0:N-1;
vc = cos(pi/N*(k+0.5)*k).*cos(pi/N*shiftdim((k+0.5)*k,-2));
vc = vc + (vc&<0).*vc*9;

兩者性能基本一致，後者可能略快

另外，其中cos的部分顯然有重複計算，不過之前其他回答也都沒有優化這部分，所以為了對比測試公平，之前也沒有考慮這部分優化。要優化只要增加個中間變數（或者寫用feval可以免去中間變數）：

function vc = fun4d
N = 40; k = 0:N-1;
vc = cos(pi/N*(k+0.5)*k);
vc = bsxfun(@times,vc,shiftdim(vc,-2)); %15b或16a可用vc = vc.*shiftdim(vc,-2);
vc = vc + (vc&<0).*vc*9;

在時間線上看到了這道題，發現還有很大的改進空間。不過為了效率，下面的方法方法利用的一些技巧。

生成數組：

With[{cn=40},vc=Outer[Times,#,#]~Transpose~{1,3,2,4}@
Cos[N[Pi/cn]Outer[Times,#+0.5,#]]@Range[0.,cn-1]];//AbsoluteTiming
(* {0.0200221,Null} *)

將負元素乘以10：

vc-9Ramp[-vc];//AbsoluteTiming
(* {0.0276519,Null} *)

基於題主的代碼，可以這麼做：

vc[[Sequence @@ #]] *= 10 /@ vp

更Mathematica的作法可以是：

vc /. x_?Negative :&> 10 x

Map[If[# &< 0, 10 #, #] , vc, {-1}]

這樣能達到你的要求不？

這裡的這個「操作」可以換成任意的其他操作

或者這樣：

Replace[vc, x_ /; x &< 0 :&> 10 x, {-1}]

或者更簡單一點：

vc /. x_ /; x &< 0 :&> 10 x

vc = Table[
Cos[([Pi] (n1 + 0.5) k1)/cn] Cos[([Pi] (n2 + 0.5) k2)/cn], {n1,
0, cn - 1}, {n2, 0, cn - 1}, {k1, 0, cn - 1}, {k2, 0, cn - 1}] /.
x_?Negative :&> 10 x

推薦閱讀：

TAG:編程 | WolframMathematica |