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.

When we don't know the number of results that will be generated, the usual way to collect results is Reap/Sow. Another alternative is linked lists. Neither of these are available in compiled functions.

AppendTo does work, but it has an $O(n^2)$ complexity, so it will be unacceptably slow for long result lists.

There was a very clever suggestion to use Internal`Bag in these situations:

Unfortunately a compiled Bag will only hold scalars. Let's try to use vectors:

cf = Compile[{}, 
  Module[{bag = Internal`Bag[{{0, 0, 0}}]}, 
   Do[Internal`StuffBag[bag, {i, i, i}], {i, {1, 2, 3}}]; 
   Internal`BagPart[bag, All]
  ]]

CompilePrint shows that this calls MainEvaluate, so this does not work.

What is the best way to collect a large number of results in a compiled function when the number of results is not known before the computation and the result type is not a (fixed size) vector or matrix?


Benchmarking the answers

  • Andy's answer

    cf = Compile[{len},
           Module[{bag = Internal`Bag[Most[{0}]]}, 
             Do[Internal`StuffBag[bag, {i, i, i}, 1], {i, len}];
             Partition[Internal`BagPart[bag, All], 3]
           ]
         ];
    
    Do[cf[500000], {100}]; // Timing
    
    (* ==> {2.87, Null} *)
    

    I needed to initialize the Bag using bag = Internal`Bag[Most[{0}]] to let Compile know that it is holding integers, not reals (see here).

  • Leonid's answer

    cf2 = 
      Compile[{len},
       Module[{arr, lim, ctr},
        arr = ConstantArray[{0, 0, 0}, 10];
        lim = Length[arr];
        ctr = 1;
        Do[
         If[ctr == lim,
          arr = Join[arr, Table[{0, 0, 0}, {lim}]];
          lim = Length[arr]];
         arr[[ctr++]] = {i, i, i},
    
         {i, len}
        ];
        Take[arr, ctr - 1]
       ]
      ];
    
    Do[cf2[500000], {100}]; // Timing
    
    (* ==> {16.474, Null} *)
    

Comparing the computational complexity of the two solutions by direct measurement:

data = Table[
   {Round[2^k], First@AbsoluteTiming@Do[cf[Round[2^k]], {100}]}, 
   {k, 13, 19, 1/2}];
data2 = Table[
   {Round[2^k], First@AbsoluteTiming@Do[cf2[Round[2^k]], {100}]}, 
   {k, 13, 19, 1/2}];

ListLogLogPlot[{data, data2}]

Mathematica graphics

(They're the same.)

share|improve this question
    
Do you mean when the number of results is unknown and the type of each result is a fixed dimensional list? Or you're asking for the more general case of different lengths and dimensions of each result to be collected? –  Rojo Feb 17 '12 at 20:17
    
@Rojo I am asking for the case when each result is the same type (a fixed dimensional tensor) but the number of results is not known beforehand. –  Szabolcs Feb 17 '12 at 21:02

2 Answers 2

up vote 12 down vote accepted

Assuming you know the dimensions of the pieces that you want to come out you can always add a second argument to Internal`StuffBag that indicates the rank of the elements going in. The result is still flat so you have to partition after the fact.

cf = Compile[{}, Module[{bag = Internal`Bag[]},
    Do[Internal`StuffBag[bag, {i, i, i}, 1], {i, {0, 1, 2, 3}}];
    Partition[Internal`BagPart[bag, All], 3]]];

Here I'm indicating that vectors will be going in to the bag. I would specify 2 for a matrix. Notice that it no longer calls MainEvaluate.

In[25]:= StringFreeQ[CompilePrint[cf], "MainEvaluate"]

Out[25]= True

In[20]:= cf[]

Out[20]= {{0., 0., 0.}, {1., 1., 1.}, {2., 2., 2.}, {3., 3., 3.}}
share|improve this answer
1  
very useful information, +1. Cheating a bit though :) –  acl Feb 18 '12 at 0:07
2  
Notice that this returns reals rather than integers. What is the best way to get it to work with integers (and also run faster)? @halirutan's trick of initializing the Bag as bag = Internal`Bag[Most[{0}]] does work, but I find it both ugly and confusing, so I was wondering if there is a better way (maybe through Compile's 3rd argument). –  Szabolcs Feb 18 '12 at 8:10
1  
@Szabolcs there may well be another way but I'm unaware of it. I typically initialize the bag as you have if I want to avoid coercion to reals. –  Andy Ross Feb 18 '12 at 17:53

As I mentioned in my recent answer, another alternative is to implement a version of a dynamic array inside Compile, as say arr = Table[{0,0},{10}] (collecting vectors of length 2, in this example) . Set up a variable (say lim) which gives the current size limit, initialize it to the initial size of the allocated array, and another one (say ctr) which counts the current maximal used position. Then, you can do something like

If[ctr==lim,
  arr = Join[arr,Table[{0,0},{lim}]];
  lim = Length[arr]
]

or, instead of doubling, use some other array expansion policy (other exponent, or additive, may depend on the problem). This adds one extra instruction (check) in your inner loop, but potentially saves you memory, and, unlike Internal`Bag, can be returned from Compile. This won't work when your results are lists of different lengths (dimensions), however.

share|improve this answer
    
This is what I would do in C too. –  Szabolcs Feb 17 '12 at 21:10
    
Can you please take a quick look at my update and verify that I did the benchmarking correctly? (Just to make sure the difference is not because of something I overlooked.) –  Szabolcs Feb 18 '12 at 8:31
1  
Leonid, congratulations on 6K rep. :-) –  Mr.Wizard Feb 19 '12 at 4:45
1  
@Mr.Wizard No, this is correct. The rep measures the level of participation on the site, as well as the quality of the answers. There are people who are many times more knowledgable than me for example, but can'y afford to participate a lot, and have a lower rep. I don't think we should pay much attention to the rep - what matters is that the questions and answers are useful to others. You bring a lot to this site, so your rep is more than well-deserved. –  Leonid Shifrin Feb 19 '12 at 13:59
1  
I really appreciate that. I still would like to meet one of these "many times more knowledgeable" people you speak of. Or perhaps not. –  Mr.Wizard Feb 19 '12 at 14:09

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.