Dynamic arrays and lists with compile (Append, Linked lists etc.)

Question

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.

Answer 1

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

Answer 2

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