Getting lengths of sublists that sum to more than one

Question

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

Answer 1

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

Answer 2

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

Answer 3

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}} *)

Answer 4

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} *)