Mathematica быстрый алгоритм 2D-биннинга
У меня возникли проблемы с разработкой подходящего быстрого алгоритма биннинга в Mathematica. У меня есть большой (~100 тыс. Элементов) набор данных в форме T={{x1,y1,z1},{x2,y2,z2},....}, и я хочу объединить его в двумерный массив размером около 100x100 бинов, причем значение бина определяется суммой значений Z, попадающих в каждый бин.
В настоящее время я перебираю каждый элемент таблицы, использую Select, чтобы выбрать, в каком бине он должен быть, основываясь на списках границ бина, и добавляю значение z в список значений, занимающих этот бин. В конце я отображаю Total в список ячеек, суммируя их содержимое (я делаю это, потому что иногда я хочу делать другие вещи, такие как максимизация).
Я пытался использовать Gather и другие подобные функции, чтобы сделать это, но вышеупомянутый метод был смехотворно быстрее, хотя, возможно, я использую Gather плохо. В любом случае, сортировка по моему методу все еще занимает несколько минут, и я чувствую, что Mathematica может добиться большего успеха. У кого-нибудь есть хороший эффективный алгоритм?
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}]]]];
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, ноль}
res3 == res2
True
"TreatRepeatedEntries" -> 1 добавляет дубликаты вверх.
Я намерен переписать код ниже из-за проблем с читабельностью Сабольча. До этого, знайте, что если ваши мусорные ведра регулярны, и вы можете использовать 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
Refactored:
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]
Вот мой подход:
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
) может обрабатывать 100 тыс. и 200 тыс. элементов в секунду соответственно на этой машине. Это достаточно быстро, или вам нужно запустить всю эту программу в цикле?
Вот мой подход с использованием функции SelectEquivalents, определенной в разделе Что находится в вашей сумке с инструментами Mathematica? который идеально подходит для такой проблемы, как эта.
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]]];