Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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.

share|improve this question
    
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
add comment

4 Answers

up vote 10 down vote accepted

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

share|improve this answer
    
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
4  
@rcollyer SparseArray[...]["Properties"] gives {"AdjacencyLists", "Background", "NonzeroPositions", "NonzeroValues", "PatternArray", "Properties"}. That's the complete list. –  Oleksandr R. Apr 4 '12 at 18:09
5  
@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
show 15 more comments

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}} *)
share|improve this answer
add comment

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
share|improve this answer
add comment

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} *) 
share|improve this answer
4  
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
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.