Mathematica快速2D分库算法

Question

我在Mathematica中开发一个合适的快速装箱算法时遇到了一些麻烦。我有一个很大的(大约100k个元素)数据集，格式为T={{x1，y1，z1}，{x2，y2，z2}，...}，我想将它放入一个大约有100x100个bin的2D数组中，其中bin值由落入每个bin的Z值之和给出。

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

我尝试过使用Gather和其他类似的函数来做这件事，但上面的方法快得离谱，尽管我可能用得不好。无论如何，按照我的方法进行排序还是需要几分钟的时间，我觉得Mathematica可以做得更好。有没有人手头有一个很好的高效算法？

Answer 1

这是一个基于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，空}

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

给出了大约{0.195228，空}

res3 == res2
True

"TreatRepeatedEntries“-> 1将重复位置相加。

Answer 2

出于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]

Answer 3

以下是我的方法：

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元素。这是否足够快，或者您是否需要在循环中运行整个程序？