Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

I have a long list of say 1 million Uniform(0,1) random numbers, such as:

 dat = {0.71, 0.685, 0.16, 0.82, 0.73, 0.44, 0.89, 0.02, 0.47, 0.65}

I want to partition this list whenever the cumulative sum exceeds 1. For the above data, the desired output would be:

{{0.71, 0.685}, {0.16, 0.82, 0.73}, {0.44, 0.89}, {0.02, 0.47, 0.65}}

I was trying to find a neat way to do this efficiently with Split combined with say Accumulate or FoldList or Total, but my attempts with Split have not been fruitful. Any suggestions?

share|improve this question
3  
Are you most concerned with "a neat way" or absolute performance? I posted the cleanest method I could think of but I doubt it will perform as well as a simple procedural loop in Compile. – Mr.Wizard Feb 6 at 15:22
1  
I like clean methods - and yours is super! Thanks. On testing, it is not as fast as it looks ... which may be because of the AddTo (which, in the past, at least, was a slow mma function). Given dat = RandomReal[{0, 1}, 10^6];, your method takes about 1.5 seconds on my Mac ... I was hoping to find something about 3 to 5 times faster. – wolfies Feb 6 at 16:13
1  
I just updated my answer with one that's a bit faster, but IMHO not enough faster to justify the ugliness. ;^) I shall try to improve it further. – Mr.Wizard Feb 6 at 16:15
1  
Seems like a borderline dupe of this one - the only missing step from there to this one in to use Internal`PartitionRagged or Mr. Wizard's dynP on the result. In any case, very strongly related. – Leonid Shifrin Feb 7 at 14:47
1  
@Leonid Wow, that is closely related, and I seem to have forgotten it despite voting for your answer. If I had remembered that question when I first saw this one I would probably have closed it without even trying to answer. Now however my Split method would be out of place there which makes it difficult for me to be happy about a close. So I'm not sure what to do or which hat to wear (moderator or user). – Mr.Wizard Feb 7 at 14:56
up vote 21 down vote accepted
dat = {0.71, 0.685, 0.16, 0.82, 0.73, 0.44, 0.89, 0.02, 0.47, 0.65};

Module[{t = 0},
  Split[dat, (t += #) <= 1 || (t = 0) &]
]
{{0.71, 0.685}, {0.16, 0.82, 0.73}, {0.44, 0.89}, {0.02, 0.47, 0.65}}

Credit to Simon Woods for getting me to think about using Or in applications like this.


Performance

I decided to make an attempt at a higher performing solution at the cost of elegance and clarity.

f2[dat_List] := Module[{bin, lns},
   bin = 1 - Unitize @ FoldList[If[# <= 1`, #, 0`] & @ +## &, dat]; 
   lns = SparseArray[bin]["AdjacencyLists"] ~Prepend~ 0 // Differences;
   Internal`PartitionRagged[dat,
      If[# > 0, Append[lns, #], lns] &[Length @ dat - Tr @ lns]
   ]
 ]

And a second try at performance using Szabolcs's inversion:

f3[dat_List] :=
  Module[{bin},
    bin = 1 - Unitize @ FoldList[If[# <= 1`, #, 0`] & @ +## &, dat];
    bin = Reverse @ Accumulate @ Reverse @ bin;
    dat[[#]] & /@ GatherBy[Range @ Length @ dat, bin[[#]] &]
  ]

Using SplitBy seems natural here but it tested slower than GatherBy.

The other functions to compare:

f1[dat_List] := Module[{t = 0}, Split[dat, (t += #) <= 1 || (t = 0) &]]

fqwerty[dat_List] :=
  Module[{f},
    f[x_, y_] := Module[{new}, If[Total[new = Append[x, y]] >= 1, Sow[new]; {}, new]];
    Reap[Fold[f, {}, dat]][[2, 1]]
  ]

fAlgohi[dat_List] :=
 Module[{i = 0, r},
  Split[dat, (If[r, , i = 0]; i += #; r = i <= 1) &]
 ]

And a single point benchmark using "a long list of say 1 million Uniform(0,1) random numbers:"

test = RandomReal[1, 1*^6];

fqwerty[test] // Length // RepeatedTiming
fAlgohi[test] // Length // RepeatedTiming
f1[test]      // Length // RepeatedTiming
f2[test]      // Length // RepeatedTiming
f3[test]      // Length // RepeatedTiming
main[test]    // Length // RepeatedTiming    (* from LLlAMnYP's answer *)
{6.270, 367891}

{1.530, 367892}

{1.24, 367892}

{0.437, 367892}

{0.821, 367892}

{0.408, 367892}

I note that qwerty's solution has one less sublist in the output because he does not include the final trailing elements if they do not exceed one. I do not know which behavior is desired.

share|improve this answer
1  
+1, I was working on it and got another version of similar idea. i=0; Split[dat, (If[r, , i = 0]; i += #1; r = i < 1) &] – Algohi Feb 6 at 15:45
    
+1, good stuff. But. I hesitate to vote for closing, but it seems a borderline dupe of this. What do you think? – Leonid Shifrin Feb 7 at 14:46

Here's my take at making a function as fast as possible.

main = Module[{idxs = sub[Accumulate@#]}, 
    Internal`PartitionRagged[#, idxs]] &;
sub = Compile[{{list, _Real, 1}},
   Block[{i, l = Length[list], ref = 1., bag = Internal`Bag[{0}]},
    For[i = 1, i <= l, i++,
     If[list[[i]] >= ref || i == l, Internal`StuffBag[bag, i]; 
      ref = list[[i]] + 1.;]
     ];
    Differences[Internal`BagPart[bag, All]]
    ]];

It's maybe 5% faster than Mr. Wizards f2 function, but the real bottleneck is PartitionRagged which takes about 80-85% of the time. I suppose, there's not much to gain from compiling, and what's needed, is a fast ragged partition routine. Part is compilable, however Compile does not like to return ragged arrays.

EDIT
This got me thinking about proper treatment of ragged arrays. While I didn't come up with any proper solution, I did manage to construct a compilable function, that creates a rectangular array with the desired output, but padded to the right with zeros.

main1 = Function[{list}, 
   Block[{sum = Accumulate[list]}, sub2[sub1[sum], list]]];
sub1 = Compile[{{list, _Real, 1}},
   Block[{i, l = Length[list], ref = 1., bag = Internal`Bag[{0}], 
     idxs},
    For[i = 1, i <= l, i++,
     If[list[[i]] >= ref || i == l, Internal`StuffBag[bag, i]; 
      ref = list[[i]] + 1.;]
     ];
    idxs = Internal`BagPart[bag, All];
    {Most[idxs] + 1, Rest[idxs], Differences[idxs]} // Transpose
    ]];
sub2 = Compile[{{idxs, _Integer, 2}, {list, _Real, 1}},
   Block[{result = 
      ConstantArray[0., {Length[idxs], Max[idxs[[All, 3]]]}], i},
    For[i = 1, i <= Length[idxs], i++,
     result[[i, ;; idxs[[i, 3]]]] = 
      list[[idxs[[i, 1]] ;; idxs[[i, 2]]]]];
    result]];

This is about 30% faster than the previous result. One might assume, that if we're talking about running totals, more often than not we're looking at non-negative numbers, so padding with zeros (or maybe a large negative number) will not lead to ambiguity.

share|improve this answer
    
I quickly realized that Internal`PartitionRagged was a bottleneck. I then realized the foolishness of my earlier comment about "a simple procedural loop in Compile" when I tried to get around it, for the reason you mention. – Mr.Wizard Feb 7 at 13:04

Internal`PartitionRagged uses Accumulate internally to generate a list of positions from the sub-list lengths, then MapThread and Take to extract the corresponding elements from the array. You can check the internal definition with

Needs["GeneralUtilities`"];
PrintDefinitions[Internal`PartitionRagged]

The reason for pointing this out is that answers which generate a list of positions and then use Differences to convert to sub-list lengths (for input to Internal`PartitionRagged) can be marginally improved by skipping the Differences and Accumulate steps.

Here's a modified version of LLlAMnYP's code. The compiled function outputs two lists, these are the start and end points for each sub-list in the original array. The main function just MapThreads Take over these lists.

fc = Compile[{{data, _Real, 1}},
   Block[{a = Accumulate[data], i, n = Length[data], ref = 1.0, bag = Internal`Bag[{0}]},
    Do[If[a[[i]] >= ref, Internal`StuffBag[bag, i]; ref = a[[i]] + 1.0], {i, n - 1}];
    Internal`StuffBag[bag, n];
    {1 + Most@Internal`BagPart[bag, All], Rest@Internal`BagPart[bag, All]}]];

f0[data_] := Module[{p = fc@data}, MapThread[Take[data, {#1, #2}] &, p]]

In my tests this comes out a few percent faster than LLlAMnYP's. As with all the answers, the bottleneck is the unpacking of the original packed data into a ragged list.

share|improve this answer
    
Very nice observation. I was curious, why Accumulate isn't on the compilable functions' list. – LLlAMnYP Feb 7 at 22:05
f[x_, y_] := Module[{new}, If[Total[new = Append[x, y]] >= 1, Sow[new]; {}, new]]
Reap[Fold[f, {}, dat]][[2, 1]]
share|improve this answer

Borrowing heavily from other answers here, but I wanted to do as much as I could inside Compile. On my machine, this is a bit faster than LLIAMnYP's main:

runsComp = Compile[{{list, _Real, 1}},
Block[{ans = ConstantArray[{0, 0}, Length@list + 1], t = 0., j = 1,
 len = Length[list]},
Do[(t += list[[i]]) <= 1 || (ans[[j + 1]] = {ans[[j, 2]] + 1, i}; 
  j++;t = 0), {i, len}
];

If[ans[[j, 2]] != len, ans[[j + 1]] = {ans[[j, 2]] + 1, len}; 
j++];
ans[[2 ;; j]]
], CompilationTarget -> "C", RuntimeOptions -> "Speed"
];

runs[list_] := list[[#1 ;; #2]] & @@@ runsComp[list]

I'm on 10.0.2 so I don't have RepeatedTiming, but on my Macbook Pro I get

test = RandomReal[1, 1*^6];

main[test] // Length // AbsoluteTiming
runs[test] // Length // AbsoluteTiming
{0.667418, 368209}
{0.593660, 368209}
share|improve this answer
    
Not to be a sore loser, but using a C compiler makes a difference - I don't have one set up. What impresses me though, is that your roundabout implementation of PartitionRagged performs on par with the built-in. +1 – LLlAMnYP Feb 7 at 15:03
    
@LLlAMnYP, Oh, I overlooked that you didn't compile to C actually! And it indeed makes a difference, to the point that a WVM version of my function is slower than yours. I guess I have just become too used to compiling to C without even trying WVM lately... – Marius Ladegård Meyer Feb 7 at 15:11
    
The question I'm taking home from today's thread, is "what is an efficient data structure for storing and creating ragged arrays?".Internal`PartitionRagged is faster than any home-brewed functions, but still pretty slow. I tried linked lists, but they are also not compilable and Bag is, if I'm not mistaken essentially a linked list. – LLlAMnYP Feb 7 at 15:17

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.