10

我在 Mathematica 中开发适当快速的分箱算法时遇到了一些麻烦。我有一个 T={{x1,y1,z1},{x2,y2,z2},....} 形式的大型（约 100k 个元素）数据集，我想将它放入一个 2D 数组中100x100 个 bin，bin 值由落入每个 bin 的 Z 值的总和给出。

目前，我正在遍历表的每个元素，使用 Select 根据 bin 边界列表选择它应该在哪个 bin 中，并将 z 值添加到占据该 bin 的值列表中。最后，我将 Total 映射到 bin 列表中，对它们的内容求和（我这样做是因为我有时想做其他事情，比如最大化）。

我曾尝试使用 Gather 和其他此类函数来执行此操作，但上述方法速度快得离谱，尽管我可能使用 Gather 很差。无论如何，按照我的方法进行排序仍然需要几分钟，我觉得 Mathematica 可以做得更好。有没有人有一个很好的高效算法方便？

4

4 回答 4

12

这是一种基于 Szabolcs 帖子的方法，速度大约快一个数量级。

data = RandomReal[5, {500000, 3}];
(*500k values*)
zvalues = data[[All, 3]];

epsilon = 1*^-10;(*prevent 101 index*)
(*rescale and round (x,y) coordinates to index pairs in the 1..100 range*)
indexes = 1 + Floor[(1 - epsilon) 100 Rescale[data[[All, {1, 2}]]]];

res2 = Module[{gb = GatherBy[Transpose[{indexes, zvalues}], First]}, 
    SparseArray[
     gb[[All, 1, 1]] -> 
      Total[gb[[All, All, 2]], {2}]]]; // AbsoluteTiming

给出大约 {2.012217, Null}

AbsoluteTiming[
 System`SetSystemOptions[ 
  "SparseArrayOptions" -> {"TreatRepeatedEntries" -> 1}];
 res3 = SparseArray[indexes -> zvalues];
 System`SetSystemOptions[ 
  "SparseArrayOptions" -> {"TreatRepeatedEntries" -> 0}];
 ]

给出大约 {0.195228, Null}

res3 == res2
True

"TreatRepeatedEntries" -> 1 增加重复的位置。

于 2011-11-20T17:37:38.310 回答
5

由于 Szabolcs 的可读性问题，我打算重写下面的代码。在那之前，要知道如果你的垃圾箱是常规的，并且你可以使用Round、Floor或Ceiling（带有第二个参数）代替Nearest，下面的代码会快得多。在我的系统上，它的测试速度比GatherBy发布的解决方案还要快。

假设我了解您的要求，我建议：

data = RandomReal[100, {75, 3}];

bins = {0, 20, 40, 60, 80, 100};

Reap[
  Sow[{#3, #2}, bins ~Nearest~ #] & @@@ data,
  bins,
  Reap[Sow[#, bins ~Nearest~ #2] & @@@ #2, bins, Tr@#2 &][[2]] &
][[2]] ~Flatten~ 1 ~Total~ {3} // MatrixForm

重构：

f[bins_] := Reap[Sow[{##2}, bins ~Nearest~ #]& @@@ #, bins, #2][[2]] &

bin2D[data_, X_, Y_] := f[X][data, f[Y][#2, #2~Total~2 &] &] ~Flatten~ 1 ~Total~ {3}

采用：

bin2D[data, xbins, ybins]
于 2011-11-18T07:21:22.797 回答
4

这是我的方法：

data = RandomReal[5, {500000, 3}]; (* 500k values *)

zvalues = data[[All, 3]];

epsilon = 1*^-10; (* prevent 101 index *)

(* rescale and round (x,y) coordinates to index pairs in the 1..100 range *)    
indexes = 1 + Floor[(1 - epsilon) 100 Rescale[data[[All, {1, 2}]]]];

(* approach 1: create bin-matrix first, then fill up elements by adding  zvalues *)
res1 = Module[
    {result = ConstantArray[0, {100, 100}]},
    Do[
      AddTo[result[[##]], zvalues[[i]]] & @@ indexes[[i]], 
      {i, Length[indexes]}
    ];
    result
    ]; // Timing

(* approach 2: gather zvalues by indexes, add them up, convert them to a matrix *)
res2 = Module[{gb = GatherBy[Transpose[{indexes, zvalues}], First]},
    SparseArray[gb[[All, 1, 1]] -> (Total /@ gb[[All, All, 2]])]
    ]; // Timing

res1 == res2

这两种方法 ( res1& res2) 可以在这台机器上分别每秒处理 100k 和 200k 个元素。这是否足够快，或者您是否需要循环运行整个程序？

于 2011-11-18T08:30:15.990 回答
3

这是我使用在您的 Mathematica 工具包中有什么中定义的函数 SelectEquivalents 的方法？这非常适合像这样的问题。

data = RandomReal[100, {75, 3}];
bins = Range[0, 100, 20];
binMiddles = (Most@bins + Rest@bins)/2;
nearest = Nearest[binMiddles];

SelectEquivalents[
   data
   ,
   TagElement -> ({First@nearest[#[[1]]], First@nearest[#[[2]]]} &)
   ,
   TransformElement -> (#[[3]] &)
   ,
   TransformResults -> (Total[#2] &)
   ,
   TagPattern -> Flatten[Outer[List, binMiddles, binMiddles], 1]
   , 
   FinalFunction -> (Partition[Flatten[# /. {} -> 0], Length[binMiddles]] &)
]

如果您想根据两个以上的维度进行分组，您可以在 FinalFunction 中使用此函数为列表结果提供所需的维度（我不记得在哪里找到它了）。

InverseFlatten[l_,dimensions_]:= Fold[Partition[#, #2] &, l, Most[Reverse[dimensions]]];
于 2011-11-18T13:40:43.820 回答