# Getting lengths of sublists that sum to more than one

I have a long list of numbers, and I'm interested in finding the length of each sublist that totals to more than 1, so for the list

{0.423768, 0.157558, 0.675251, 0.685209, 0.580772, 0.0230333,
0.927156, 0.506085, 0.0516773, 0.485349}


I would get the lengths {3, 2, 3}, as

Total[{0.423768, 0.157558, 0.675251}] == 1.25658
Total[{0.685209, 0.580772}] == 1.26598
Total[{0.0230333, 0.927156, 0.506085}] == 1.45627


The last two numbers total to 0.537026, so we ignore them. This has an easy imperative solution:

thresholdFor[nums_List] :=
Module[{i = 0, sum = 0, k = 0},
Reap[
For[i = 1, i <= Length@nums, ++i,
sum += nums[[i]];
++k;
If[sum > 1,
Sow[k];
k = 0;
sum = 0]]][[-1, 1]]];
thresholdFor[{}] = {};


In addition to being kind of ugly and needing the For, it's also pretty slow, as functions that index into lists so often are. It takes about 0.04 seconds to process a list of 10000 random numbers (picked uniformly between 0 and 1). I futzed around with LengthWhile, TakeWhile and the like before deciding tha I really needed the full generality of Fold to accomplish what I needed to do:

thresholdFold[nums_] :=
Flatten@Last@Fold[
With[{sum = #[[1]] + #2, length = #[[2]] + 1, acc = #[[3]]},
If[sum > 1,
{0, 0, {acc, length}},
{sum, length, acc}]] &,
{0, 0, {}},
nums]


This is arguably more idiomatic, but it's even a bit slower than thresholdFor. I can speed thresholdFor up a lot by compilation (with a suitable adaptation to get rid of the Reap/Sow pair) and a little wrapper to handle the empty list properly:

compiledBody =
Compile[{{nums, _Real, 1}},
With[{n = Length@nums},
Module[{result = ConstantArray[0, n], sum = 0.0, k = 0, fill = 0, i},
For[i = 1, i <= n, ++i,
sum += nums[[i]];
++k;
If[sum > 1.0,
result[[++fill]] = k;
k = 0; sum = 0]];
Take[result, fill]]]];

thresholdCompiled[nums_List] := compiledBody[nums];
thresholdCompiled[{}] = {};


This is dramatically faster after compilation--it runs in about 10 ms for the list of 100000 numbers, and is only about 10 times slower than Mean or Total on the same data. Still, I always think I'm going down a bit of a blind alley when I start using Compile to make imperative list-processing algorithms fast enough to use.

Finally, I did come up with something more Mathematica-esque, but it's actually much slower than the compiled solution (taking about 130 ms):

thresholdClip[nums_] :=
Differences@Flatten@Position[
FoldList[Clip[Plus[##], {0.0, 1.0}, {0.0, 0.0}] &, 0.0, randoms],
0.0]


I've tried a couple tweaks for speeding up thresholdClip, including one which replaces Position with Pick so I can compile the whole thing, but that didn't seem to do a lot of good.

-
The problem looks naggingly familiar. Have we had a question like this before? – Sjoerd C. de Vries Apr 4 '12 at 16:42
@SjoerdC.deVries It felt naggingly familiar to me, too, but searching didn't turn up anything helpful-looking. – Pillsy Apr 4 '12 at 16:43

Very good question / problem. Generally, this problem seems to belong to the class of problems where Compile is the best choice if maximum efficiency is looked for, since it is, by its nature, not a good fit for the Mathematica paradigm of working with lots of data at once. However, your last solution can be, in a somewhat modified form, brought to the same performance level while only partly compiled:

Clear[split, thresholdSemiCompiled];
split = Compile[{{nums, _Real, 1}},
FoldList[If[# > 1, 0, #] &[#1 + #2] &, 0.0, nums],
CompilationTarget -> "C"];

thresholdSemiCompiled[l_List] :=
Flatten@Differences[SparseArray[Unitize[split[l]] - 1]["NonzeroPositions"]]


It turns out that Clip was a major bottleneck in the compiled split, and once I replaced it with If, it became an order of magnitude faster. This one is on par with your fastest fully compiled solution, in terms of performance.

### Edit (by Oleksandr R., at Leonid's request)

In version 8, provided that an integer argument is used rather than a pattern, Position can be compiled down into a call to a new kernel/VM function "Position". This allows the above to be implemented in purely compiled code and avoids the need to use SparseArray in an undocumented and perhaps unintuitive way to get the positions of the zeros without unpacking the input array:

thresholdCompiled = Compile[{{nums, _Real, 1}},
Flatten@Differences@Position[
Unitize@FoldList[If[# > 1, 0, #] &[#1 + #2] &, 0.0, nums],
0
], CompilationTarget -> "C"
];


The performance gained by doing this is, perhaps surprisingly, relatively minor: this version returns timings of about 7/8ths those of the semi-compiled approach, irrespective of compilation to MVM bytecode or to C, probably due to the ability to locate zeros directly rather than having to subtract 1 from the array and find nonzero elements, and by avoiding the need to return a large intermediate array. (Unitize, Differences, Position, and Flatten are all compilable, but in version 8 they are compiled into simple function calls, so the only situation in which we would expect a performance gain through compilation is when, as here, it allows a larger problem to be dealt with entirely in compiled code.)

-
Last question (probably), I've never encountered SparseArray[...]["NonzeroPositions"] before, where did you find that, and are there other methods available on it? – rcollyer Apr 4 '12 at 17:46
@rcollyer I learned this from Oliver, who published some code using this a few times on Mathgroup, e.g. here. There are also "NonzeroValues", I don't know what else. But note, all these shortcuts are not strictly necessary, just convenient. One can use the SparseArray "API" like the one I posted here to achieve the same effect. – Leonid Shifrin Apr 4 '12 at 17:53
Convenient, yes. Undocumented, definitely. They don't even show up via SubValues[SparseArray]. :P – rcollyer Apr 4 '12 at 17:58
@rcollyer SparseArray[...]["Properties"] gives {"AdjacencyLists", "Background", "NonzeroPositions", "NonzeroValues", "PatternArray", "Properties"}. That's the complete list. – Oleksandr R. Apr 4 '12 at 18:09
@rcollyer yes, I think the lack of sensible output for InterpolatingFunction[...]["Properties"] is an omission. Anyway, the list of possible properties for interpolating functions is {"Domain", "Coordinates", "DerivativeOrder", "InterpolationOrder", "Grid", "ValuesOnGrid"}. – Oleksandr R. Apr 4 '12 at 18:29

Possibly not fastest, but concise and using a simple iteration over the list. Or two iterations, if we do a cleanup step as in the Compile'd variants.

subseqlens[ll_] := Reap[Module[{n = 0, tot = 0.},
Map[(n++; tot += #; If[tot > 1, tot = 0.; Sow[n]; n = 0]) &,
ll]]][[-1, 1]]

subseqlensC = Compile[{{ll, _Real, 1}}, Module[{n = 0, tot = 0., m},
Select[
Map[(n++; tot += #; If[tot > 1, tot = 0.; m = n; n = 0; m, 0]) &,
ll], # != 0 &]]];

subseqlensCC = Compile[{{ll, _Real, 1}}, Module[{n = 0, tot = 0., m},
Select[
Map[(n++; tot += #; If[tot > 1, tot = 0.; m = n; n = 0; m, 0]) &,
ll], # != 0 &]], CompilationTarget -> "C"];


Here is a fairly big example.

In[461]:= biglist = RandomReal[1, 10^6];
Timing[ss1 = subseqlens[biglist];]
Timing[ss2 = subseqlensC[biglist];]
Timing[ss3 = subseqlensCC[biglist];]
ss1 === ss2 === ss3

Out[462]= {2.92, Null}
Out[463]= {0.26, Null}
Out[464]= {0.04, Null}
Out[465]= True

-
"but concise" - as Fortran is to assembly language. What century are we in? – alancalvitti Feb 10 at 5:59
@alancalvitti Sometimes Fortran-like code is the way to proceed. I mean, if it was good enough for Michaelangelo, then who am I to snub it? – Daniel Lichtblau Feb 10 at 15:54
Why not assembly, or punched cards? – alancalvitti Feb 11 at 19:38
@alancalvitti We only developed Compile when we started running low on punched cards. – Daniel Lichtblau Feb 11 at 20:51
But now you're running low on time - should develop a RTOS or RTVM running on FPGAs. You see how long it takes to dev JSF (F35) blocks? – alancalvitti Feb 12 at 5:31

Using Fold with Reap and Sow might be cleaner:

ClearAll[f];
f = If[Total[{##}, Infinity] > 1, Sow[Length[Flatten[{##}]]]; ## &[], {##}] &;
Fold[f, First@list, Rest@list]; // Reap // Last

(* Out[1]= {{3, 2, 3}} *)

-

You can use NestWhileList together with Accumulate :

alist = {0.423768, 0.157558, 0.675251, 0.685209, 0.580772,
0.0230333, 0.927156, 0.506085, 0.0516773, 0.485349} ;


This will return the position where you need to split the input list :

NestWhileList[
{pos = Position[Accumulate[#[[2]]], _?(# >= 1 &)][[1, 1]];
pos + #[[1]],Drop[#[[2]], pos]} &,
{pos = Position[Accumulate[alist], _?(# >= 1 &)][[1, 1]], Drop[alist, pos]},
(Length[#[[2]]] > 1 && Total[#[[2]]] >= 1) &
][[All, 1]]

(* {3, 5, 8} *)

-
The performance on this solution is quadratic, because you keep re-accumulating and re-searching alist after droping smallish chunks off the front of it each time. Also, if you're going to be tracking state in NestWhileList by using a list argument, you don't need to mutate a variable like pos` as you go. – Pillsy Apr 4 '12 at 17:02