6
$\begingroup$

I have a code that I need to compile. The code contains a loop in which I need to gather data into a list. The question is how to best gather data so that the code works inside Compile. The following pieces of code illustrate the problem:

1. gathering just one list (tbl1) into the final list (result) - WORKS

Needs["CCompilerDriver`"];
dir = "c:\\temp";
$CCompilerInternalDirectory = dir;
func = Compile[{},
   Module[{t, tmax, dt, n, eps, tbl1, tbl2, i, result},
    n = 10;
    t = 0.0;
    tmax = 10.0;
    eps = 10.0^-2;
    tbl1 = Table[0.0, {n}];
    tbl2 = Table[0.0, {n}];
    result = {tbl1};
    While[TrueQ[t <= tmax],
     dt = RandomReal[{0.1, 0.5}];
     t = t + dt;
     For[i = 1, i <= n, i++,
      tbl1[[i]] = t^0.5;
      tbl2[[i]] = t^2;
      ];
     result = Append[result, tbl1];
     ];
    result
    ],
   CompilationTarget -> "C", 
   CompilationOptions -> {"InlineExternalDefinitions" -> True}
   ];

2. gathering two lists (tbl1 and tbl1) - NOT Working:

Needs["CCompilerDriver`"];
dir = "c:\\temp";
$CCompilerInternalDirectory = dir;
func = Compile[{},
   Module[{t, tmax, dt, n, eps, tbl1, tbl2, i, result},
    n = 10;
    t = 0.0;
    tmax = 10.0;
    eps = 10.0^-2;
    tbl1 = Table[0.0, {n}];
    tbl2 = Table[0.0, {n}];
    result = {tbl1, tbl2};
    While[TrueQ[t <= tmax],
     dt = RandomReal[{0.1, 0.5}];
     t = t + dt;
     For[i = 1, i <= n, i++,
      tbl1[[i]] = t^0.5;
      tbl2[[i]] = t^2;
      ];
     result = Append[result, {tbl1, tbl2}];
     ];
    result
    ],
   CompilationTarget -> "C", 
   CompilationOptions -> {"InlineExternalDefinitions" -> True}
   ];

Compile::cpts: The result after evaluating Insert[result,{tbl1,tbl2},-1] should be a tensor. Nontensor lists are not supported at present; evaluation will proceed with the uncompiled function. >>

Compile::cpts: The result after evaluating Insert[result,{tbl1,tbl2},-1] should be a tensor. Nontensor lists are not supported at present; evaluation will proceed with the uncompiled function. >>

3. And the same with linked lists - NOT working either:

Needs["CCompilerDriver`"];
dir = "c:\\temp";
$CCompilerInternalDirectory = dir;
func = Compile[{},
   Module[{t, tmax, dt, n, eps, tbl1, tbl2, i, result},
    n = 10;
    t = 0.0;
    tmax = 10.0;
    eps = 10.0^-2;
    tbl1 = Table[0.0, {n}];
    tbl2 = Table[0.0, {n}];
    result = {tbl1, tbl2};
    While[TrueQ[t <= tmax],
     dt = RandomReal[{0.1, 0.5}];
     t = t + dt;
     For[i = 1, i <= n, i++,
      tbl1[[i]] = t^0.5;
      tbl2[[i]] = t^2;
      ];
     result = {result, {tbl1, tbl2}};
     ];
    result
    ],
   CompilationTarget -> "C", 
   CompilationOptions -> {"InlineExternalDefinitions" -> True}
   ];

Compile::cset: Variable result of type {_Real,2} encountered in assignment of type {_Real,3}. >>

Compile::cset: Variable result of type {_Real,2} encountered in assignment of type {_Real,3}. >>

4. Trying Sow and Reap with no effect either (although these functions are just a huge mystery for me):

Needs["CCompilerDriver`"];
dir = "c:\\temp";
$CCompilerInternalDirectory = dir;
func = Compile[{},
   Module[{t, tmax, dt, n, eps, tbl1, tbl2, i, result},
    n = 10;
    t = 0.0;
    tmax = 10.0;
    eps = 10.0^-2;
    tbl1 = Table[0.0, {n}];
    tbl2 = Table[0.0, {n}];
    result = {tbl1, tbl2};
    result = Reap[
      While[TrueQ[t <= tmax],
        dt = RandomReal[{0.1, 0.5}];
        t = t + dt;
        For[i = 1, i <= n, i++,
         tbl1[[i]] = t^0.5;
         tbl2[[i]] = t^2;
         ];
        Sow[tbl1];
        Sow[tbl2];
        ];
      ];
    result
    ],
   CompilationTarget -> "C", 
   CompilationOptions -> {"InlineExternalDefinitions" -> True}
   ];

Compile::cset: Variable result of type {_Real,2} encountered in assignment of type {_Real,0}. >>

Compile::cset: Variable result of type {_Real,2} encountered in assignment of type {_Real,0}. >>

So what is he best and Compile safe approach of accumulating data to a list if one does not know the list size beforehand? (I can see that making a sufficiently large list to store the data would solve the problem, but I consider this as a last resort, and rather bad approach; Append is a compilable function after all.

$\endgroup$
  • 3
    $\begingroup$ "Internal`Bag inside Compile" is a good start. $\endgroup$ – C. E. Mar 22 '15 at 21:45
  • 1
    $\begingroup$ There is also this question $\endgroup$ – Leonid Shifrin Mar 22 '15 at 21:49
  • 1
    $\begingroup$ I looked at that, ok, it might be a solution. There is still a question of returning multiple lists from a compiled function. Let's say I want to return a scalar and an array, so I return {scalar,list} but this does not work, Mathematica complains Compiled expression ... should be a rank 2 tensor of machine-size real numbers. So how do I do that? $\endgroup$ – leosenko Mar 22 '15 at 23:42
1
$\begingroup$

When Compile is used the list returned must have a consistent structure.
A version of your 2. code block that does compile without errors:

func = Compile[{}, Module[{t, tmax, dt, n, eps, tbl1, tbl2, i, result},
 n = 10;
 t = 0.0;
 tmax = 10.0;
 eps = 10.0^-2;
 tbl1 = Table[0.0, {n}];
 tbl2 = Table[0.0, {n}];
 result = {{tbl1, tbl2}}; (* <== extra {} added here *)
 While[t <= tmax,
  dt = RandomReal[{0.1, 0.5}];
  t = t + dt;
  For[i = 1, i <= n, i++, tbl1[[i]] = t^0.5;
   tbl2[[i]] = t^2;];
  result = Append[result, {tbl1, tbl2}];];
 result], CompilationTarget -> "C", 
CompilationOptions -> {"InlineExternalDefinitions" -> True}];

This compiled function returns a tensor of rank 3.
If you evaluate the Module of your function you'll see that the first two entries are a single list, whereas all other entries and the entries of func as defined above are lists of lists.


Regarding your comment, if there is no "natural" way (like there is in your example above) to make get a regular array as the output, you'll have to make the output a regular array.
For a scalar and a list one could use something like

 testF1 = Compile[{},
           {RandomReal[{1, 10}, 1]~Join~Table[0, {5 - 1}],
            RandomReal[{1, 10}, 5]}]

{scalar, list} = MapAt[First@# &, testF1[], 1]

{8.44659, {8.53816, 3.61359, 8.87295, 8.27914, 7.3135}}

and for two lists of different length

 testF2 = Compile[{},
           {RandomReal[{1, 10}, 2]~Join~Table[0, {5 - 2}],
            RandomReal[{1, 10}, 5]}]

{firstList, lastList} = MapAt[#[[;; 2]] &, testF2[], 1]

{{5.61392, 1.83448}, {6.85195, 8.52678, 4.18129, 2.10639, 9.0375}}

$\endgroup$
  • $\begingroup$ Thank you @Karsten 7. However, I am not sure what it means consistent structure. From what I have seen it seems, that if I have two lists (arrays) which do not have the same dimensions, there is no way to return those. So let's say, if I want to return a list and a scalar or a list that has dimensions 1x10 and another one that has dimensions 100x10, how would I do that? $\endgroup$ – leosenko Mar 25 '15 at 20:31
  • $\begingroup$ @leosenko Please see my edit. To the best of my current knowledge (and I'd like to be proven wrong on this) the only way possible is to make the array a regular one (meaning all sublists having the same length). $\endgroup$ – Karsten 7. Mar 25 '15 at 21:08
  • 1
    $\begingroup$ @leosenko Check How to force Compile to return multiple results? for other possible implementations. $\endgroup$ – Karsten 7. Mar 25 '15 at 21:15
1
$\begingroup$

This variant works. But it will of necessity do some evaluation outside of Compile so it may not offer much in the way of a speed gain (I do not have time to check right now). Another drawback is I have not succeeded in getting it to cooperate with compiling to C.

SetAttributes[myreap, HoldFirst];
myreap[a__] := Reap[a][[2, 1]]

func = Compile[{}, Module[{t, tmax, dt, n, eps, tbl1, tbl2, i, result},
    n = 10;
    t = 0.0;
    tmax = 10.0;
    eps = 10.0^-2;
    tbl1 = Table[0.0, {n}];
    tbl2 = Table[0.0, {n}];
    result = {tbl1, tbl2};
    result = 
     myreap[While[TrueQ[t <= tmax], dt = RandomReal[{0.1, 0.5}];
       t = t + dt;
       For[i = 1, i <= n, i++, tbl1[[i]] = t^0.5;
        tbl2[[i]] = t^2;];
       Sow[tbl1];
       Sow[tbl2]]];
    result], {{myreap[__], _Real, 2}},
   CompilationOptions -> {"InlineExternalDefinitions" -> True}];

func[]

(* {{0.408651552874, 0.408651552874, 0.408651552874, 0.408651552874, 
  0.408651552874, 0.408651552874, 0.408651552874, 0.408651552874, 
  0.408651552874, 0.408651552874}, 
...
 {102.197228542, 102.197228542, 102.197228542, 
  102.197228542, 102.197228542, 102.197228542, 102.197228542, 
  102.197228542, 102.197228542, 102.197228542}} *)
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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