# Optimising 2D binning code

I have a set of (x,y,z) data, 45,000 to be precise and I want to bin the z values in 256 equidistant bins based on their (x,y) values. The final array should be a set of 256x256 array with each slot containing an average of binned z values.

Being new to mathematica, I came up with the following code:

data = RandomReal[{12000, 35000}, {45000, 3}];
data1 = data[[All, {1, 2}]];(*strip the zvalues from the set*)
xValues = data[[All, 1]];
yValues = data[[All, 2]];
zValues = data[[All, 3]];
(*Compute maximum/minimum of x values*)
maxXvalue = Max[xValues];
minXvalue = Min[xValues];

(*Compute maximum/minimum of y values*)
maxYvalue = Max[yValues];
minYvalue = Min [yValues];

(*Compute maximum/minimum of z values*)
maxZvalue = Max[zValues];
minZValue = Min[zValues];

bbx = {Floor[minXvalue], Floor[maxXvalue],
Floor[((maxXvalue - minXvalue)/256)]}; (* equidistant x bins*)
bby = {Floor[minYvalue], Floor[maxYvalue],
Floor[((maxYvalue - minYvalue)/256)]};(* equidistant y bins*)
bList = BinLists[data1, {bbx}, {bby}];
bCount = BinCounts[data1, {bbx}, {bby}];(*Gives a count of the number of items in \
each bins*)

(*Defining array to contain final z average values*)
meanZValues = Table[0, {Length[bList]}, {Length[bList]}];

i = 0; (*initialising loop variables*)
j = 0;
k = 0;

f[x_] := zValues[[x]];(*Defining function to get z values back*)

For[i = 1, i <= Length[bList], i++,
For [j = 1, j <= Length[bList], j++, m1 = {};    (*Re-empty m1 list*)
For [k = 1, k <= Length[bList[[i, j]]], k++,
AppendTo[m1, Position[data1, bList[[i, j]][[k]]] (*accessing only the x-
coordinate index of the position on original matrix*)
];
(*Getting the indices of the binned values*)
indices = Flatten[DeleteDuplicates[Take[m1, All]]]; (*Position command above gives multiple indices if  these values occur more than once, hence deleting the duplicate ones*)

meanZValues[[i, j]] =  Mean[Map[f,indices]];  (*Compute average values of Z by accessing the original array, getting the z values  *)
]
]
]
meanZValues

It gives an output in a reasonable amount of time for up to couple of thousand values, however, it lags and maybe crashes without any output for 45,000 set of data.

How do I make this code more efficient? Thank you

-
– faysou Jan 14 '13 at 6:51
"bList = BinLists[data1, {bbx}, {bby}]; bCount = BinCounts[data1, {bbx}, {bby}];"should be "bList = BinLists[data1, bbx, bby]; bCount = BinCounts[data1, bbx, bby];"? – kptnw Jan 14 '13 at 8:54

Modifying @ruebenko's answer in the StackOverflow Q/A linked in Faysal's comment (Mathematica fast 2D binning algorithm) to get the means of z-values for each bin (using yet another undocumented setting for the option "TreatRepeatedEntries" that works in version 9 only):

zvalues = data[[All, 3]];
epsilon = 1*^-10;
indexes = 1 + Floor[(1 - epsilon) 256 Rescale[data[[All, {1, 2}]]]];
SystemSetSystemOptions["SparseArrayOptions" -> {"TreatRepeatedEntries" -> (Mean[{##}] &)}];
binmeansZ = SparseArray[indexes -> zvalues];
SystemSetSystemOptions["SparseArrayOptions" -> {"TreatRepeatedEntries" -> First}];

A picture:

MatrixPlot[binmeansZ]

Update: Timings

Mr.Wizards's version 7 settings (also works in versions 8.0.4.0 and 9):

SetSystemOptions["SparseArrayOptions" -> {"TreatRepeatedEntries" -> 2}];
AbsoluteTiming[binmeans =  Normal[SparseArray[indexes -> zvalues]] /.
"List"[x__] :> Mean@{x};]
SetSystemOptions["SparseArrayOptions" -> {"TreatRepeatedEntries" -> 0}];
(* {0.086009, Null} *)

Version 9 settings:

SystemSetSystemOptions["SparseArrayOptions" -> {"TreatRepeatedEntries" -> (Mean[{##}] &)}];
AbsoluteTiming[binmeansZ = SparseArray[indexes -> zvalues];]
SystemSetSystemOptions["SparseArrayOptions" -> {"TreatRepeatedEntries" -> First}];
(* {0.035003, Null}*)
binmeansZ == SparseArray[binmeans]
(* True *)

Update 2: Default settings in versions 8.0.4.0 and 9:

"TreatRepeatedEntries" /. SystemOptions["SparseArrayOptions"][[1, 2]]
(* 0          (Version 8.0.4.0) *)
(* First      (Version 9)   *)
-
How did you learn of that option? – Mr.Wizard Jan 14 '13 at 9:20
@Mr.W learned about SparseArrayOptions from ruebenko's answer in the linked Q/A. That the option setting could be a pure function was a wild guess and it worked. – kglr Jan 14 '13 at 9:24
This doesn't work on version 7, so sadly I cannot vote for it, but that's fantastic if it works for others. – Mr.Wizard Jan 14 '13 at 9:27
@Mr.W, do you get anything that looks similar from SystemOptions["SparseArrayOptions"] in version 7? – kglr Jan 14 '13 at 9:29
I think I've got a workaround. I'll post an answer soon, and I hope you'll compare the performance of the two for me. – Mr.Wizard Jan 14 '13 at 9:34
data = RandomReal[{12000, 35000}, {45000, 3}];
n = 256;
dataT = Transpose@data;
r[x_, m_] :=  IntegerPart@N@Rescale[x, {Min[dataT[[m]]], Max[dataT[[m]]]}, {1, n + 1}]
Timing[(Mean /@ Transpose@#) & /@ GatherBy[
data /. {x_, y_, z_} -> {r[x, 1], r[y, 2], z}
/. {n + 1, x__} -> {n, x}
/. {x_, n + 1, z_} -> {x, n, z}, {#[[1]], #[[2]]} &]][[1]]

1.188

-

kguler's answer looks great but unfortunately it doesn't work on version 7.
However, I was able to find a similar method that does.

data = RandomReal[{12000, 35000}, {45000, 3}];

zvalues = data[[All, 3]];
epsilon = 1*^-10;
indexes = 1 + Floor[(1 - epsilon) 256 Rescale[data[[All, {1, 2}]]]];

SetSystemOptions["SparseArrayOptions" -> {"TreatRepeatedEntries" -> 2}];

AbsoluteTiming[
binmeans = Normal[SparseArray[indexes -> zvalues]] /. "List"[x__] :> Mean@{x};
]

SetSystemOptions["SparseArrayOptions" -> {"TreatRepeatedEntries" -> 0}];

MatrixPlot[binmeans]

{0.0300000, Null}

-
The AbsoluteTiming i get is 0.085000 (versus 0.030003 with the setting (Mean[{#}]&).(So, the setting 2 gives "List" of the values?) +1 of course... – kglr Jan 14 '13 at 9:54
@kguler Thank you. Setting 2 does give "List"[val1, val2, . . .] which seems like an odd format but still useful. It's not surprising the additional operation slows this down. Do you know if your method works in v8 or only v9? – Mr.Wizard Jan 14 '13 at 9:59
Just checked version 8.0.4.0: the method in my post works only for version 9. – kglr Jan 14 '13 at 10:47
@Mr.Wizard, i tried your codes and ran it with my data. I tried with 10,000 data and 100x100 bins. The final array of z values has a lot of zeros (implying empty bins). I ran the same data using the software Gwyddion which binned the data and returned 100x100 z values without a single zero in the array. Is there anything wrong with the code above? Thank you – Mun Jan 16 '13 at 5:41
@Mun if you had the same number of bins as data and no empty bins then there must have been one and only one item in each and every bin... is this realistic? While I wouldn't rule out a mistake, this remarkable coincidence makes me a little suspicious of the result given by Gwyddion. – Oleksandr R. Jan 16 '13 at 10:42