Efficient conditional Mean[] on a large data set

Question

Consider a list, AK6, that has 382 sublists of length varying from 2500 to 3000. Each "subsublist" is as such : AK6[[1,1]]={5.5,1001}

With AK6[[All,1]] going from 1001 to 4001 with gaps (missing points)

I now need to compute the Mean[] for all the sublists with equal second value. That is Select with Flatten[AK6,1] all the list with 2005 as the second value (from 350 to 384 items approx)

I do this using :

ParallelTable[
             Mean[Select[Flatten[AK6s, 1], #[[2]] == gazeNo &][[All,1]]], 
             {gazeNo, Range[1000, 4000, 1]}]

It does what I need but this is very slow. Is there a way to do this computation faster?

---

Download the 13 MB list as

AK6 = Uncompress@Import["http://api2.ge.tt/0/9F8d6WD/0/blob/download", "String"];

Answer 1

You could do something like this

mean = Reap[Sow @@@ Flatten[cogAK6, 1];, _, {Mean[#2], #} &][[2]];

This will be a lot faster than your approach because by using Sow and Reap this code only iterates through the list of data once. In your code, you reiterate through all elements of the data list for every value of gazeNo (so 3000 times instead of only once).

Answer 2

This answer will be about efficiency, no ambitions for the beauty contest participation here. Since you mentioned speed, and only need mean values, here is a method that will be an order of magnitude faster and order of magnitude more memory-efficient than the (fine) solutions offered so far.

Code

First, we make a compiled function generator:

ClearAll[generateFastMean];
generateFastMean[maxIndex_Integer?Positive] :=
   Compile[{{data, _Real, 2}},
     Module[{means = Table[0., {maxIndex}], num = Table[0, {maxIndex}], 
      ctr = 0, i = 0, index = 0, resultIndices =  Table[0, {maxIndex}],
      indexHash = Table[0, {maxIndex}]},
     Do[
       index = IntegerPart[data[[i, 2]]];
       means[[index]] += data[[i, 1]];
       num[[index]]++;
       If[indexHash[[index]] == 0,
          indexHash[[index]] = 1;
          resultIndices[[++ctr]] = index
       ];,
       {i, Length[data]}
     ];
     resultIndices = Take[resultIndices, ctr];
     Transpose[{resultIndices, means[[resultIndices]]/num[[resultIndices]]}]
    ], (* Module *)
    CompilationTarget -> "C", RuntimeOptions -> "Speed"];

What happens here is that I exploit the fact that your indices are not too large integers, and use arrays as hash-tables to keep the data I accumulate. At the end, I extract only those indices which are present, and compute means by dividing a list of totals by a list of index frequencies. Since your indices are smaller than 5000, I will generate the function accordingly:

fastMean = generateFastMean[5000]; 

The case at hand

First, we load the data (by the way, compressing is a fine idea but one can have additional savings, see below):

AK6 = Uncompress@Import["http://api2.ge.tt/0/9F8d6WD/0/blob/download", "String"];

To reduce the memory this occupies, we Flatten the list as:

AK6PFl1 =  Join @@Map[Developer`ToPackedArray, AK6]

Just to illustrate what kind of savings I am speaking about:

ByteCount /@ {AK6PFl, Flatten[AK6, 1]}

(*
  ==> {18238440, 127667944}
*)

So, we are looking at the order of magnitude here. Ok, now computing:

fastMean[AK6PFl1] //Short[#,5]&//AbsoluteTiming

(*
  ==> {0.063,{{1000.,5.00499},{1001.,5.00541},{1002.,5.00556},<<2995>>,
       {3998.,3.72284},{3999.,3.72058},{4000.,3.71838}}}
*)

How to store the data

While Compress is a good idea, I would supplement it with converting sublists to packed arrays, as illustrated above. So, you can do

Export["C:\\Temp\\AK6.dat", Compress@Map[Developer`ToPackedArray, AK6], "String"]

Then, the import is fast, and there is no unpacking (which happens with the current version of the file):

Uncompress@Import["C:\\Temp\\AK6.dat","String"]//ByteCount//AbsoluteTiming

(*
 ==> {0.594,18305896}
*)

The idea of combining Compress and packed arrays (in a somehwat more complex setting since I made the data self-uncompressing), allowed me to get the remarkable performance in this case, and seems to represent a general method.

Remarks

I cheated because I did not include the compilation time, which is significant. However, if you have to do this many times, this may pay off. Of course, this will all make sense only when you need extreme performance, but I think it is good to understand what your performance limitations really are, for a given problem (I do not claim that the code has been fully optimized).

Answer 3

Heike already got Sow/Reap so here is GatherBy:

<< cogAK6.m;
j = Join @@ cogAK6;
j[[All, 2]] = Round@j[[All, 2]];

{#[[1, 2]], Mean[#[[All, 1]]]} & /@ GatherBy[j, Last] // Sort

Both methods are fast because you do not repeatedly rescan the data for each gazeNo.

By the way I added Round just to be sure that values were properly gathered in case there were slight variations.

Answer 4

Here is another compiled effort (this will be ugly):

ClearAll[f]
f = Compile[
   {
    {vals, _Real, 1},
    {ind, _Integer, 1}
    },
   Module[
    {
     s = ConstantArray[0., 4000],
     n = ConstantArray[0, 4000]
     },
    Do[
     s[[ind[[i]]]] += vals[[i]];
     n[[ind[[i]]]]++,
     {i, 1, Length@ind}
     ];
    Transpose[{Drop[s, 999]/N@Drop[n, 999], Range[1000, 4000]}]
    ],
   CompilationTarget \[Rule] "C"
   ];

then

({vals, ind} = Transpose@Flatten[AK6, 1];
  mean2 = f[vals, ind];) // AbsoluteTiming

takes 0.223913s.

EDIT: comparing to Leonid's code (which looks less ugly than this to me), this is a tiny bit faster (Leonid's code here takes 0.26s if I include the generateFastMean[5000] part (or maybe not--I should run it many times and compare).

Answer 5

I have another approach and since it is no update of my first answer, I post it separately. This one is short: We know from here how to get/check compile-able functions in Mathematica. Since we can compile Transpose, Sort, Reverse, Table, Mean, Tally and IntegerPart we can simply write:

Get["tmp/ak6.mx"];
AK6P = Join @@ Map[Developer`ToPackedArray, AK6];

meanFirst = Compile[{{lst, _Real, 2}}, 
Module[{
  vals = Most[{1.0}], ids = Most[{1.0}], ptr = 1},
  {ids, vals} = Transpose[Sort[Reverse /@ lst]];
  Table[
    With[{currentMean = Mean[vals[[ptr ;; ptr + elm[[2]] - 1]]]},
      ptr = ptr + elm[[2]];
      {elm[[1]], currentMean}
    ], {elm, Tally[IntegerPart[ids]]}]
  ], CompilationTarget -> "C", RuntimeOptions -> "Speed"
];

First[AbsoluteTiming[meanFirst[AK6P]]]

(*  0.199104 *)

The good thing is, that Sort works on tuples as expected, when the id is upfront. Then we use Tally to cluser the sorted ids and the we jump cluster-wise through the values and calculate the Mean.

This function has not call to MainEvaluate and runs completely in library code.

Update

As already pointed out in the comments, the code is not really superior to its non-compiled counterparts. The difference is, that even with Compile the code contains almost only high-level list-manipulations and is not broken down into c-like code, which is not often the case with Compile.

Answer 6

Let me give one more approach which kind of divides the work to do. The core will be a meanFirst compiled function which takes a list of the form {{val1, id}, {val2, id}, ...} and calculate the mean of all val and returns their mean together with the id. When we use Gather as already done so by @MrWizard to take all input-tuples and creating sublists which all have a common id, then what's left to do is to map meanFirst over all these sublists.

meanFirst = Compile[{{l, _Real, 2}},
   {Mean[First[Transpose[l]]], l[[1, 2]]}, CompilationTarget -> "C", 
   Parallelization -> True, RuntimeAttributes -> {Listable}];

At this point we reconsider the Map and remember, that compiled functions which are parallelized and are Listable can use very fast pthread-parallelization without any explicit Map. How to do that? Easy, just give the compiled function a list or matrix of the argument-type it usually would get. Our meanFirst gets a sublst of type {{val1, id}, {val2, id}, ...}, so the only thing we need to do is to supply a list {sublst1, sublst2, ...} as argument and all mean-calculations are done in parallel.

Unfortunately there is one Wermutstropfen to this approach. The supplied tensor is not allowed to be ragged, which means in this case all sublst need to have the same length. What we can do is Gather the list {sublst1, sublst2, ...} by the length of all sublst and we get {{sbl11, sbl21, ...}, {sbl12, sbl22, ...}, ...} where all sblN1 have the same length. Now we can securely call meanFirst[{sbl11, sbl21, ...}] on all gathered sublists (meaning we Map).

Unfortunately, this means we get only a partial parallelization, but it is, nevertheless reasonable fast:

Get["tmp/ak6.mx"];
AK6P = Join @@ Map[Developer`ToPackedArray, AK6];
meanFirst /@ 
   Gather[Gather[AK6P, Equal[#1[[2]], #2[[2]]] &], 
    Length[#1] === Length[#2] &]];