如何加速Mathematica替换矩阵元素

Question

我有几个100x15矩阵; 其中一个是距离.当该矩阵的元素超过一个边界时,我想将这些元素重置为零,并将其他三个矩阵的相应元素重置为零.这是我的愚蠢方式(但它有效):

Do[ If[ xnow[[i, j]] > L, xnow[[i, j]] = 0.;
                  cellactvA[[i, j ]]  = 0.;
                  cellactvB[[i, j ]]  = 0.;
                  cellactvC[[i, j ]]  = 0.;   ], (* endIF  *)
   { i, 1, nstrips}, {j, 1, ncells}       ];  (* endDO *)

我试过了ReplacePart:

 xnow = ReplacePart[ xnow, Position[ xnow, x_?(# > L &) ] ]

(类似这样的东西,我没有得心应用;它已经足够正确地执行了),但是它和循环一样慢并且没有在矩阵xnow中产生正确的替换结构.请告知如何以相当快的方式执行此操作,因为此计算是在执行多次的另一个循环(随着时间的推移)内.当然,整体计算现在非常缓慢.提前致谢.

---

以下是我在R中的表现; 非常简单快捷:

    # -- find indices of cells outside window
indxoutRW  <- which( xnow > L, arr.ind=T )

    # -- reset cells outside window
cellrateA[indxoutRW] <- 0 
cellrateB[indxoutRW] <- 0 
cellrateC[indxoutRW] <- 0 

    # -- move reset cells back to left side
 xnow[indxoutRW]    <- xnow[indxoutRW] - L  

Answer 1

这个怎么样:

Timing[
 matrixMask2 = UnitStep[limit - $xnow];
 xnow = $xnow*matrixMask2;
 cellactvA2 = $a*matrixMask2;
 cellactvB2 = $b*matrixMask2;
 cellactvC2 = $c*matrixMask2;
 ]

如果你想编写快速代码一件事,以确保检查On ["Packing"]不给出消息; 或者至少你了解它们并知道它们不是问题.

编辑OP评论:

mask = UnitStep[limit - xnow];
{xnow*mask, cellactvA2*mask, cellactvB2*mask, cellactvC2*mask}

希望这有帮助,你仍然需要设置限制.

Answer 2

以下将基于SparseArrays,避免无关紧要的东西,速度非常快:

extractPositionFromSparseArray[
   HoldPattern[SparseArray[u___]]] := {u}[[4, 2, 2]];
positionExtr[x_List, n_] := 
   extractPositionFromSparseArray[
     SparseArray[Unitize[x - n], Automatic, 1]]

replaceWithZero[mat_, flatZeroPositions_List, type : (Integer | Real) : Real] :=
  Module[{copy = Flatten@mat},
     copy[[flatZeroPositions]] = If[type === Integer, 0, 0.];
     Partition[copy, Last[Dimensions[mat]]]];

getFlatZeroDistancePositions[distanceMat_, lim_] :=
  With[{flat = Flatten[distanceMat]},
     With[{originalZPos = Flatten@ positionExtr[flat , 0]},
       If[originalZPos  === {}, #, Complement[#, originalZPos ]] &@
         Flatten@positionExtr[Clip[flat , {0, lim}, {0, 0}], 0]]];

现在,我们生成矩阵,确保它们已经打包:

{xnow, cellactvA, cellactvB, cellactvC} = 
   Developer`ToPackedArray /@ RandomReal[10, {4, 100, 15}];

这是1000次这样做的基准:

In[78]:= 
Do[
  With[{L = 5},
    With[{flatzpos = getFlatZeroDistancePositions[xnow,L]},
       Map[replaceWithZero[#,flatzpos ]&,{xnow,cellactvA,cellactvB,cellactvC}]]
  ],
  {1000}]//Timing

Out[78]= {0.203,Null}

请注意,在此过程中没有解压缩,但您必须确保从一开始就打包矩阵,并为该replaceWithZero函数选择正确的类型(整数或实数).