How to force Compile to return multiple results?

26

14

I would like to return two separate values from a compiled function, but Mathematica refuses to use the compiled version. Here is a narrowed-down example:

    cFunc = Compile[{{a, _Integer, 1}}, {a, a.a}];

(* 
  ==> CompiledFunction::cfex: Could not complete external evaluation at instruction 3; proceeding with uncompiled evaluation. >>

  ==> {{1, 1, 1}, 3}
*)

Here a and a.a are hypothetic intermediate results of some long computation, that must be externalized in order for computation to go on outside the compiled function. I assume it is not possible to return lists that are not uniform (data type and dimensions). Is there any way to return more than one results from a compiled computation?

István Zachar

Posted 2012-02-15T14:13:27.780

Reputation: 44 135

Answers

16

What about making the result uniform inside Compile and constructing things back afterwards? Like

   cFunc = Compile[
       {{a, _Integer, 1}}, 
       Join @@ {a, {a . a}}]; 
cFunc[{1, 3}]
Function[z, {Most[z], 
       Last[z]}][%]

which has the nice feature of not calling MainEvaluate:

    Needs["CompiledFunctionTools`"]; cFunc // CompilePrint

(* 
==>
    1 argument
        2 Integer registers
        3 Tensor registers
        Underflow checking off
        Overflow checking off
        Integer overflow checking on
        RuntimeAttributes -> {}

        T(I1)0 = A1
        I0 = 4
        Result = T(I1)2

1   I1 = DotVV[ T(I1)0, T(I1)0, I0]]
2   T(I1)1 ={ I1 }
3   T(I1)2 = Join[ T(I1)0, T(I1)1]]
4   Return

*)

Rolf Mertig

Posted 2012-02-15T14:13:27.780

Reputation: 16 237

I couldn't be more specific because your method is the one that I'm still using regularly, 4 years later : ) – István Zachar – 2016-07-21T17:23:57.687

2I think your code is going to break in this case: cFunc = Compile[{{a, _Integer, 1}}, Join @@ {1.*a, {a.a}}];It makes some assuptions about a, that may not be true. Also, a single call to MainEvaluate is not so much of an issue, as long as it is not in a loop. – None – 2012-02-15T15:04:54.803

@ruebenko what is your analysis of my methods? – Mr.Wizard – 2012-02-15T15:11:56.867

@ruebenko: Definitely it is in a loop... But then the problem with Rolf's code is that joining and retrieving partial results from z might not be that easy. – István Zachar – 2012-02-15T23:55:47.373

why not? can you be more specific? – Rolf Mertig – 2012-02-16T00:16:14.657

@RolfMertig I am just learning how to use Compile. Why does the code f:=Compile[{{M,_Integer,2},{m,_Integer}}, Module[{l},l=ConstantArray[{},m]; l[[M[[1]]]]=Transpose@{M[[2]]}; l] ]; f[{{8,9,10},{2,3,1}},10] return CompiledFunction::cfn: Numerical error encountered at instruction 8; proceeding with uncompiled evaluation.? – Leo – 2020-11-21T16:20:10.997

ConstantArray[{0},m] will work. But please ask a new question in general if you do have a new, unrelated, question. – Rolf Mertig – 2020-11-22T21:38:29.500

14

You can use a technique like this:

cf = Compile[{{n, _Integer}, {m, _Integer}}, 
   Module[{cpos = RandomSample[Range[n], m]},
    Set[pos, cpos];
    RandomReal[1, Length[cpos]]]];

{valres, posres} = Block[{pos}, {cf[10, 5], pos}]

The set works on a variable which is scoped with Block. Since this is done only once, the call to MainEvaluate is not a problem.

user21

Posted 2012-02-15T14:13:27.780

Reputation:

Am I correct in thinking that any solution must force a (one-time) call to MainEvaluate because compiled functions can only work with non-ragged arrays? – Szabolcs – 2012-02-15T14:20:20.050

I think there are special cases where this is not the case (e.g. Listable compiled functions that returns arrays of different length), but to the best of my knowledge, yes, you need a call to MainEvaluate. – None – 2012-02-15T14:24:44.287

@Szabolcs I guess it might be possible without the MainEvaluate call if in future support for passing Internal`Bags between the Mathematica VM and the interpreter is implemented. But for now, this is indeed the only way to do it as far as I know. – Oleksandr R. – 2012-02-15T14:27:05.487

@Oleksandr Actually I think the OP's version works correctly (i.e. uses the compiled code) as well, except it issues a warning message when it gets to the last step. Wouldn't it be better to just suppress that message by wrapping the call in Quiet[..., {CompiledFunction::cfex}]? Perhaps slightly dangerous, as one never knows if the same message would be triggered by an unexpected input ... – Szabolcs – 2012-02-15T14:30:35.247

3@Szabolcs the problem there is that the error encountered by the VM, even though it occurs at the very end of the compiled code, forces a complete re-evaluation in the interpreter. Error control for compiled code isn't very sophisticated, unfortunately, so evaluation won't pick up in the interpreter where it left off in the VM. – Oleksandr R. – 2012-02-15T14:39:46.483

2@Oleksandr I see. I wasn't sure if that was the case. cf = Compile[{a}, Print[a]]; cf[2] prints twice and shows that this is indeed what happens. But now it is not clear to me when MainEvaluate will force a complete re-evaluation and when it won't. – Szabolcs – 2012-02-15T14:42:27.370

Let me correct myself. If I run cf = Compile[{a}, Print[a]]; cf[2] in a fresh kernel, it prints twice and gives a warning. If I switch away from the Mathematica window (!!!) and back, it only prints once and does not give a warning. @ruebenko Do you know what is going on? I'm on Win7. – Szabolcs – 2012-02-15T14:44:39.723

@Szabolcs MainEvaluate is fine when the VM expects it, such as when the compiler has emitted code to evaluate an uncompilable function (e.g. Sow; Set for a non-localized variable). But when errors are encountered that were not anticipated by the compiler, this message is produced and complete re-evaluation occurs. – Oleksandr R. – 2012-02-15T14:50:14.490

What is the a reason you use the FullForm Set in your code? – Mr.Wizard – 2012-02-15T15:38:49.023

@Mr.Wizard, the sole reason for this: When I woke up this morning and planed my day (LOL), was such that you could ask about it ;-) no, seriously, none, whatsoever. – None – 2012-02-15T16:09:55.963

8

A couple of ideas:

cfun = Compile[{{a, _Integer, 1}}, x[1] = a; x[2] = a.a;]
fun[a_?VectorQ] :=
 Block[{x},
  cfun[a];
  {x[1], x[2]}
 ]

fun[{1,2,3,4,5}]
{{1, 2, 3, 4, 5}, 55}
cf2 = Compile[{{a, _Integer, 1}}, Sow[a]; Sow[a.a];];

Reap[cf2@{1, 2, 3, 4, 5}][[2, 1]]
{{1, 2, 3, 4, 5}, 55}

Mr.Wizard

Posted 2012-02-15T14:13:27.780

Reputation: 259 163

2which also calls MainEvaluate – Rolf Mertig – 2012-02-15T14:37:21.777

@Rolf but it does not result in CompiledFunction::cfex isn't that the point? By the way I voted for yours; slick. – Mr.Wizard – 2012-02-15T14:40:11.583

@Mr.Wizard: nice ideas from you, too! Voted for it. Just try to understand how things work and what are the differences. – Rolf Mertig – 2012-02-15T14:42:06.090

@Rolf good point. I don't have CompilePrint in v7. I think your method will be the best if the function is low complexity. On a computationally complex function I think the overhead that my methods have will not matter, but I haven't tested that assertion. – Mr.Wizard – 2012-02-15T14:45:55.523

@Mr.Wizard, I think your methods are fine, you need to call MainEvaluate one more time than with my approach. (That's the only thing I could think of...) – None – 2012-02-15T15:36:40.930

@ruebenko thank you. – Mr.Wizard – 2012-02-15T15:38:43.143

@Mr.Wizard, do you have no access to v8? – Rolf Mertig – 2012-02-16T00:14:15.627

@Rolf I'd have to buy it and I am too cheap. I am waiting for v9. – Mr.Wizard – 2012-02-16T04:12:51.253

@Mr.Wizard why don't you ask WRI if you get a moderator-discount? They probably also should add you to the v9-tester-list (I hope I am still on that list ...) – Rolf Mertig – 2012-02-16T09:04:19.667

6

I have been encountering this problem a lot recently. Since many related questions are linked to this post and I will share my solution here.

Advantages

  • Can return multiple results (not necessarily have the same shape)

  • Can return ragged lists

  • Does not call MainEvaluate

Basic idea

Flatten all the tensors into one single list, and include enough information to reconstruct them.

The first element of my list is the number of tensors / variables to return.

The $2$nd to $2 + var - 1$ th element corresponds to the rank of each tensor

The $2 + var$ to the $2 + var + rank_i -1$ th elements corresponds to the dimension of each tensor

Construction inside Compile

1. Multiple return with different-dimension tensors (of different types)

Note:In my example there is no Complex or True|False, but since Re, Im and Boole are all compilable, they can be transformed to a real tensor and a integer tensor respectively.

This example illustrates returning 3 tensors with different dimensions.

cf1=Compile[{},
  Module[{
   m={{0,8,1,7},{1,9,2,6}},
   n={0.301,0.98},
   p={{{1,0},{2,7}},{{2,0},{0,0}}}},
  Join[{3},
   {TensorRank[m]},{TensorRank[n]},{TensorRank[p]},
   Dimensions[m],Dimensions[n],Dimensions[p],
   Flatten@m,Flatten@n,Flatten@p]
 ]
]

2. Return a ragged list (of arbitrary length)

This is a common case when a collection of Positions should be returned. This example illustrates adding 1D list of arbitrary length to the result programmatically.

cf2=Compile[{},Module[{var=0,rank={},dim={},res={},temp},
  Do[temp=RandomReal[{0,1},RandomInteger[{1,10}]];
   var++;
   AppendTo[rank,TensorRank[temp]];
   dim=Join[dim,Flatten@Dimensions[temp]];
   res=Join[res,Flatten@temp];
  ,{i,1,3}];
 Join[{var},rank,dim,res]]]

Neither of the examples have MainEvaluate when examining with CompilePrint.

Extracting the lists

extractLists[list_?VectorQ] := 
 Module[{vars = Round@First@list, rank, dim}, 
  rank = Round@list[[2 ;; 1 + vars]]; 
  dim = Round@
    Internal`PartitionRagged[
     list[[2 + vars ;; 1 + vars + Total@rank]], rank]; 
  MapThread[
   ArrayReshape, {Internal`PartitionRagged[
     list[[2 + vars + Total@rank ;;]], Times @@@ dim], dim}]]

The results (the result of cf2 is random):

extractLists[cf1[]]
(*{{{0., 8., 1., 7.}, {1., 9., 2., 6.}}, {0.301, 
  0.98}, {{{1., 0.}, {2., 7.}}, {{2., 0.}, {0., 0.}}}}*)
extractLists[cf2[]]
(*{{0.895086, 0.716247, 0.626751, 0.457065, 0.709812, 0.118539, 
  0.504491, 0.40369}, {0.2376}, {0.159539, 0.398285, 0.0233042, 
  0.246191, 0.351316, 0.580408}}*)

Notes

The type of the result is not conserved (Integer is converted to Real). This can be implemented by adding extra parameter before the rank info (I did not include it because it's not useful for my cases). Also I am not sure whether the performance of the code inside Compile is optimal. Feel free to edit if there are improvements.

vapor

Posted 2012-02-15T14:13:27.780

Reputation: 7 611

Hmm... Isn't this the same method Rolf has suggested in his answer? I.e. flatten internal results and reconstruct outside the Compile. – István Zachar – 2016-07-21T17:24:34.157

@IstvánZachar true for your original problem. I meant to solve a general one. – vapor – 2016-07-21T17:27:10.960

1Yes, but still, the general idea was quite clear from Rolf's post. Well, it's nevertheless nice to have a more general approach, so here goes my +1. – István Zachar – 2016-07-21T17:34:05.040