Resorting data from a multidimensional list

Question

I have a dataset that is generated in a 4D table. A snippet would look something like:

{{{{x,y},{x2,y2},...}...}...} and so on, in code form:

data=Table[{RandomReal[{0, 4}], RandomReal[{0, 1}]}, {InTh, 5}, {InE, 6}, backEn, 780}, {backAngle, 0, 89}]

Now, I want to sort/rebin the data in such a way that a program goes over the list and sums y values which lie in a certain x range, let's say my bins are generated by using:

driftBins=Table[i,{i,0,4,0.001}]

The approach I currently have is really poor and uses For loops as I am not really familiar with Mathematica to an expert level.

driftRebinned2 = {};
start = AbsoluteTime[];
For[i = 1, i < Length[driftBins], i++,
  summy = 0;
  For[InTh = 1, InTh <= 5, InTh++,
   For[InE = 1, InE < 6, InE++,
    For[backEn = 1, backEn <= 780, backEn++,
     For[backAngle = 1, 
      backAngle <= 90, backAngle++,
      If[driftBins[[i + 1]] > 
         data[[InTh, InE, backEn, backAngle]][[1]] && 
        driftBins[[i]] <= data[[InTh, InE, backEn, backAngle]][[1]], 
       summy = summy + data[[InTh, InE, backEn, backAngle]][[2]]]
      ]
     ]
    ]
   ]; AppendTo[
   driftRebinned2, {(driftBins[[i]] + driftBins[[i + 1]])/2, summy}];
  ];
AbsoluteTime[] - start

Is there a way to optimize this code or do it in a table? I tried doing it via a table but the issue I encounter there is I am not sure how to sum up y values that lie within a certain x range.

The data generation is just an example. My real data is a result of reading a data file where the values of the data file are weighted by a combination of weighing function which depend on the four iterators.

Answer 1

I'll assume here that the bins are equally spaced.

Starting like @HenrikSchumacher suggests:

data = Join[RandomReal[{0, 4}, {5, 6, 780, 90, 1}], RandomReal[{0, 1}, {5, 6, 780, 90, 1}], 5];
data1 = Flatten[data, 3];

If your data need to be weighted depending on their four indices, you could instead first do something like

weighteddata = MapIndexed[{#1[[1]], f[#1[[2]],#2]}&, data, {4}];

with f some function you define that does the weighting. Then use data1 = Flatten[weighteddata, 3] instead.

For each data point, round the $x$-coordinate down to the nearest thousandth (you may use Round or Ceiling instead of Floor to define the bins, depending on what exactly you need):

data2 = {Floor[#[[1]], 0.001], #[[2]]} & /@ data1;

Gather together all data points with the same $x$-coordinate (rounded down to the nearest thousandth):

A = GatherBy[data2, First];

From these gathered lists, calculate for each bin (i) the $x$-value of the lower bin edge (or the center of the bin, or whatever) and (ii) the sum of the $y$-values (or the mean, or length, or whatever):

B = (#[[1, 1]] -> Total[#[[All, 2]]]) & /@ A;

The results in B are in random order. Make an ordered list of all bins and their sum:

bins = Range[0, 4, 0.001];
Transpose[{bins, Lookup[B, bins, 0]}]

{{0, ...}, {0.001, ...}, {0.002, ...}, ..., {4, ...}}

(The dots stand for the sum values in each bin.)

If you prefer having the bin centers as first coordinate instead of the bin lower limits, then you could do

bins = Range[0, 3.999, 0.001];
Transpose[{bins+0.0005, Lookup[B, bins, 0]}]

{{0.0005, ...}, {0.0015, ...}, {0.0025, ...}, ..., {3.9995, ...}}

I agree with @HenrikSchumacher that BinLists is nicer than my use of GatherBy here.

Answer 2

I haven't checked this for correctness, but in runs in about a second:

data = Join[
   RandomReal[{0, 4}, {5, 6, 780, 90, 1}],
   RandomReal[{0, 1}, {5, 6, 780, 90, 1}],
   5
   ];

Because the nested structure of data is not need anywhere, I flatten it to a matrix. Afterwards I employ Nearest to find the correct bin for each entry of data1 (exploting that the bins are equally spaced). Afterwards, a compiled routine assemble adds the contributions for each bin into a large vector.

data1 = Flatten[data, 3];

driftBins = Table[i, {i, 0., 4., 0.001}];
bincenters = MovingAverage[driftBins, 2];
binradii = 0.5 (driftBins[[2]] - driftBins[[1]]);
idx = Developer`ToPackedArray[
  Nearest[bincenters -> Automatic, data1[[All, 1]], {1, binradii}]
  ];


assemble = 
 Compile[{{idx, _Integer, 1}, {y, _Real, 1}, {n, _Integer}},
  Block[{a},
   a = Table[0., {n}];
   Do[a[[idx[[i]]]] += y[[i]], {i, 1, Length[idx]}];
   a
   ]
  ];

symmylist = assemble[Flatten[idx], data1[[All, 2]], Length[bincenters]];
driftRebinned2 = Transpose[{bincenters, symmylist}];

A somewhat more flexible approach is by using BinLists:

sum = Compile[{{a, _Real, 2}},
   Total[a[[All, 2]]],
   RuntimeAttributes -> {Listable},
   Parallelization -> True
   ];

symmylist3 = sum[
   Flatten[
    BinLists[Flatten[data, 3];, {driftBins}, {{-∞, ∞}}], 
    1
    ]
   ];
driftRebinned3 = Transpose[{bincenters, symmylist3}];