# Error when compiling function containing InterpolationFunction

I'm trying to solve a coupled ODE system using NDSolve. The coupling can turn on and off in some regions, and when the coupling is off, the solution can be written down directly in terms of simple algebraic functions. So I only need to solve the regions when the coupling is on using NDSolve. Then the solutions in different regions are put together to form a Piecewise function as a whole solution. This Piecewise function now contains several InterpolationFunction objects and other simple algebraic functions. Since the InterpolationFunction only takes a small part of the whole solution region, I'm trying to compile this Piecewise function so that those simple algebraic functions would benefit from compiling.

Here is the simplified version of the problem:

suppose we have the solution in one region as an interpolation function

s = C1 /.
First@NDSolve[{I C1'[t] == C2[t] E^(-I t),
I C2'[t] == C1[t] E^(I t), C1[0.] == 1., C2[0.] == 0.}, {C1,
C2}, {t, 0., 10.}]
(*
InterpolatingFunction[{{0.,10.}},<>]
*)


and the solution in full region is constructed as

sol[t_] :=
Piecewise[{{s[t], t <= 10.}, {Sin[t], t <= 20.}, {Cos[t], t <= 30}}, 0.]


it can compile without problem, but when executed gives error

f2 = Compile[{{t, _Real}}, Evaluate@sol[t], {{_InterpolationFunction, _Complex}}];
f2[5.]


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

(*
-0.787325-0.231525 I
*)


if we compile without Evaluate, then there is no error

f1=Compile[{{t,_Real}},sol[t],{{_sol,_Complex}}];
f1[5.]
(*
-0.787325-0.231525 I
*)


but this seems to send the whole sol function to the MainEvaluator and thus those simple algebraic functions can't benefit from the compiling. Indeed if we print out the content we can see this difference.

CompilePrint@f1

"
1 argument
1 Real register
1 Complex register
Underflow checking off
Overflow checking off
Integer overflow checking on
RuntimeAttributes -> {}

R0 = A1
Result = C0

1   C0 = MainEvaluate[ Hold[sol][ R0]]
2   Return
"

CompilePrint@f2

"
1 argument
3 Boolean registers
1 Integer register
9 Real registers
Underflow checking off
Overflow checking off
Integer overflow checking on
RuntimeAttributes -> {}

R0 = A1
R2 = 7.
R1 = 10.
R4 = 20.
I0 = 30
R7 = 0.
Result = R5

1   B0 = R0 <= R1 (tol R2)
2   if[ !B0] goto 6
3   R3 = MainEvaluate[ Hold[InterpolatingFunction[{{0., 10.}}, <>]][ R0]]
4   R5 = R3
5   goto 20
6   B1 = R0 <= R4 (tol R2)
7   if[ !B1] goto 11
8   R5 = Sin[ R0]
9   R6 = R5
10  goto 19
11  R6 = I0
12  B2 = R0 <= R6 (tol R2)
13  if[ !B2] goto 17
14  R6 = Cos[ R0]
15  R8 = R6
16  goto 18
17  R8 = R7
18  R6 = R8
19  R5 = R6
20  Return
"


So how can I solve this problem? Thanks.

-
I've never seen the pattern you used; f2 = Compile[{{t, _Real}}, sol[t], {{sol[_], _Complex}}] seems to work. – b.gatessucks Jul 10 '13 at 6:24
@b.gatessucks, if you make this a question I could talk a bit about this ;-) – user21 Jul 10 '13 at 7:13
@ruebenko Do consider it a question then ! – b.gatessucks Jul 10 '13 at 7:21

Here is one way to do it:

You could compile the components one by one and inline those in a final compiled function.

c1 = With[{help = s}, Compile[{{t, _Complex}}, help[t]]];
c2 = Compile[{{t, _Complex}}, Sin[t]];
c3 = Compile[{{t, _Complex}}, Cos[t]];
f3 = Compile[{{t, _Complex}},
Piecewise[{{c1[t], t <= 10.}, {c2[t], t <= 20.}, {c3[t], t <= 30}},
0.]
, CompilationOptions -> {"InlineExternalDefinitions" -> True}]


Note, that I used _Complex as an input argument. That will automatically call the correct type signature and a complex type will be returned. Here then the conversion to complex is done when during the compiled function call.

Also note that the f3 has the compiled functions c2 and c3 optimized out:

CompilePrint@f3

1 argument
3 Boolean registers
1 Integer register
5 Real registers
6 Complex registers
Underflow checking off
Overflow checking off
Integer overflow checking on
RuntimeAttributes -> {}

C0 = A1
R2 = 7.
R0 = 10.
R3 = 20.
I0 = 30
R1 = 0.
Result = C3

1   C1 = R0 + R1 I
2   B0 = C0 <= C1 (tol R2)
3   if[ !B0] goto 7
4   C1 = MainEvaluate[ Function[{t}, CompiledFunction[{t}, \
InterpolatingFunction[{{0., 10.}}, <>][t], -CompiledCode-][t]][ C0]]
5   C3 = C1
6   goto 26
7   C2 = R3 + R1 I
8   B1 = C0 <= C2 (tol R2)
9   if[ !B1] goto 14
10  C2 = C0
11  C3 = Sin[ C2]
12  C4 = C3
13  goto 25
14  R4 = I0
15  C2 = R4 + R1 I
16  B2 = C0 <= C2 (tol R2)
17  if[ !B2] goto 22
18  C2 = C0
19  C4 = Cos[ C2]
20  C2 = C4
21  goto 24
22  C5 = R1 + R1 I
23  C2 = C5
24  C4 = C2
25  C3 = C4
26  Return


To simplify this a bit you could use:

c1 = With[{help = s}, Compile[{{t, _Complex}}, help[t]]];
f3 = Compile[{{t, _Complex}},
Piecewise[{{c1[t], t <= 10.}, {Sin[t], t <= 20.}, {Cos[t],
t <= 30}}, 0.]
, CompilationOptions -> {"InlineExternalDefinitions" -> True}]


Hope this helps.

-
Thanks, but I think maybe I didn't make my self clear. I was trying to ask how to compile the whole solution sol using only sol, suppose we only have access to sol but not s. – xslittlegrass Jul 10 '13 at 15:51
@xslittlegrass, you always can get the content of the Piecewise by using List @@ sol[t] – user21 Jul 10 '13 at 15:57
yes, that's true, but in my real sol function, there are lots of different s's it would much easier to be able to compile the whole thing without sorting out each s and assign new symbols. By the way, why do you use "InlineExternalDefinitions" option? It seems it doesn't matter whether we use that option or not. And do you know why there is the error "complete external evaluation at instruction 3" in my question? I guess I'm also very interested in why it doesn't work if we have InterpolationFunction in compile. – xslittlegrass Jul 10 '13 at 16:03
@xslittlegrass, look at the CompilePrint@f3 to see the difference. – user21 Jul 10 '13 at 16:08
Thanks, but I still get an error "Compile::argcomp: The comparison, LessEqual, is invalid for arguments of type Complex. >>". The kernel version is "8.0 for Linux x86 (64-bit) (October 10, 2011)" Do you know what's wrong? – xslittlegrass Jul 10 '13 at 21:51

This works(version 9 mac)

f3 = Compile[{{t, _Complex}}, Evaluate@sol[t]]

In[4]:= f3[5.]
Out[4]= -0.787325-0.231525 I

CompilePrint@f3

1 argument
3 Boolean registers
1 Integer register
5 Real registers
6 Complex registers
Underflow checking off
Overflow checking off
Integer overflow checking on
RuntimeAttributes -> {}

C0 = A1
R2 = 7.
R0 = 10.
R3 = 20.
I0 = 30
R1 = 0.
Result = C2

1   C1 = R0 + R1 I
2   B0 = C0 <= C1 (tol R2)
3   if[ !B0] goto 7
4   C1 = MainEvaluate[ Hold[InterpolatingFunction[{{0., 10.}}, <>]][ C0]]
5   C2 = C1
6   goto 24
7   C2 = R3 + R1 I
8   B1 = C0 <= C2 (tol R2)
9   if[ !B1] goto 13
10  C2 = Sin[ C0]
11  C3 = C2
12  goto 23
13  R4 = I0
14  C3 = R4 + R1 I
15  B2 = C0 <= C3 (tol R2)
16  if[ !B2] goto 20
17  C3 = Cos[ C0]
18  C4 = C3
19  goto 22
20  C5 = R1 + R1 I
21  C4 = C5
22  C3 = C4
23  C2 = C3
24  Return

-