# Efficient conditional Mean[] on a large data set

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?

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

-
I compressed and re-uploaded the file---it was way too big at 35 MB. –  Szabolcs Feb 12 '12 at 14:23
@Szabolcs, how did you compress the file? Worth a new question ? –  500 Feb 12 '12 at 14:42
Export["data.dat.gz", Compress[AK6], "String"]. Compress gives an ASCII (i.e. 7-bit) string, so I did an additional gzip compression in Export to shave off some 10-12%. Ideally one could just Export to WDX (which I think is already compressed) but with this amount of data it is just too slow. –  Szabolcs Feb 12 '12 at 15:02
@Szabolcs, thank you, how do you then import the data in the correct format ? –  500 Feb 12 '12 at 15:04
You can use data = Uncompress@Import["data.dat.gz", "String"];. Import and Export will automatically handle compressing to or uncompressing from gz if they detect that the file is gzip-compressed. –  Szabolcs Feb 12 '12 at 15:08

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).

-
Thank You very much. I am in shock by how fast this was. Could you explain what makes this so fast? –  500 Feb 12 '12 at 14:13
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 gazNo (so 3000 times instead of only once). –  Heike Feb 12 '12 at 14:21
@Mr.Wizard Just replacing _ with Range[1000,4000] should do it I suppose (and probably checking for an empty list when applying Mean). I agree that it's a powerful technique, so I thank you just in case you are the one who introduced it to me. –  Heike Feb 12 '12 at 14:40
Yes, exactly. There is no need to check for empty lists as the function is not applied of nothing is Reaped. You can also nest this construct to make it multidimensional. –  Mr.Wizard Feb 12 '12 at 14:45
@Mr.Wizard and that form is eerily similar to my general one. –  rcollyer Feb 12 '12 at 18:45

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[DeveloperToPackedArray, 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[DeveloperToPackedArray, 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).

-
When you want to time especially c-compiled stuff, you should always use AbsoluteTiming, everything else is not what you want. –  halirutan Feb 12 '12 at 19:55
@halirutan That's a good point, but in my experience, this is more relevant for running external programs (e.g. Mathlink). C-compiled function becomes a part of the kernel process, so there should be no real difference between using Timing on it and some built-in function. In particular, I re-checked with AbsoluteTiming and got the same results. But I agree with you in general, and will edit accordingly. –  Leonid Shifrin Feb 12 '12 at 20:26
You have to be careful when using the parallelization in compiled functions, because it can happen that the cpu-times of each processors are added up. Let me try an example code: f = Compile[{{i, _Integer}}, i^2, CompilationTarget -> "C", Parallelization -> True, RuntimeAttributes -> {Listable}]; With[{arg = RandomInteger[{0, 100}, 10^7]}, First[#[f[arg]]] & /@ {Timing, AbsoluteTiming}] –  halirutan Feb 12 '12 at 21:54
@halirutan Yes, I agree - for this case, probably because of threading. I do know about this, and had I used parallel, I would have used AbsoluteTiming from the start. –  Leonid Shifrin Feb 12 '12 at 22:03

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.

-
Wizard, thank you for your attention ! –  500 Feb 12 '12 at 14:37

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).

-
Actually, it is somewhat more economical than mine, since I have to maintain one more list and do a few more things, although the idea is similar. Note that most of the time is spend in Transpose-Flatten in your method, since the sublists of AK6 are unpacked. –  Leonid Shifrin Feb 12 '12 at 17:25
@Leonid I know, but this was the fastest I could get with a reasonable effort (although what "reasonable effort" achieves depend on the person exerting the effort!) –  acl Feb 12 '12 at 17:27
+1. But, if you compare, you should also include compilation time in your case! (I only have to generate the function once). You code is somewhat faster, but for a different reason, which I outlined above, and to see that, you have to work on packed list, e.g. like I did in my answer. –  Leonid Shifrin Feb 12 '12 at 17:27
@acl Are aware of the fact, that calls to ConstantArray are never compiled? –  halirutan Feb 12 '12 at 19:49
@halirutan yes, but ConstantArray here takes something like $10^{-6}$ seconds, while in the body of the function we iterate over something $10^6$ times. So I thought the main bottleneck would be the loop, not the initialization of the array. –  acl Feb 12 '12 at 20:19

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[DeveloperToPackedArray, 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.

-
This seems a bit slower than the other compiled solutions, which is strange (maybe because of the CopyTensors?). –  acl Feb 12 '12 at 23:15
Yes, probably. The nice thing here is IMO only, that a lot of highlevel functions are used and everything is done in compiled code. –  halirutan Feb 13 '12 at 0:19
+1 for the idea. Note however that your code, when uncompiled (top-level), runs only twice slower. This is because top-level Sort with default comparison function is so efficient that does not benefit from compilation at all. The two-fold speedup is because here Table does benefit from compilation. In this case, you can take the line Transpose[Sort[Reverse /@ lst]] to the top-level, compile the rest, and get identical performance. My point is that compilation is most effective when you know which pieces will benefit from it the most. One big Compile should not be the goal IMO. –  Leonid Shifrin Feb 13 '12 at 6:37
@LeonidShifrin Yes, the problem is not one of the nicest. It's badly parallelizable and can be implemented with Tensor-transformations which are fast even without Compile. One advantage of having it all inside compiled code is, that you can now easily call it on several data-sets in parallel. –  halirutan Feb 13 '12 at 8:29
@halirutan Yes, that would make sense. Also, you may just use it on its own (take the C code), if you link against Wolfram RTL - so in some cases having everything inside Compile is good. –  Leonid Shifrin Feb 13 '12 at 9:28

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[DeveloperToPackedArray, AK6];
meanFirst /@
Gather[Gather[AK6P, Equal[#1[[2]], #2[[2]]] &],
Length[#1] === Length[#2] &]];

-
I have a 6-core machine and your code runs roughly with the same speed as other non-compiled solutions. A little profiling tells us that most time is spent in the nested Gather- function (because the list unpacks), and this sort of destroys the purpose of both compilation and parallelization, in my view. In other words, the main culprit here is not a single core, but unpacking, so if you allow unpacking, no amount of parallelization (unless you have a 100 CPU machine) would give you a speed-up w.r.t the single-core compiled code that does not unpack. –  Leonid Shifrin Feb 12 '12 at 20:41
you can see that it unpacks by doing SetSystemOptions[ "PackedArrayOptions" -> "UnpackMessage" -> True]. I like the idea though, so you have my vote. –  acl Feb 12 '12 at 20:53