4
$\begingroup$

This question already has an answer here:

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.

$\endgroup$

marked as duplicate by Alexey Popkov, MarcoB, user9660, Yves Klett, Mr.Wizard Jul 13 '16 at 6:45

This question has been asked before and already has an answer. If those answers do not fully address your question, please ask a new question.

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

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.

$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – xslittlegrass Jul 10 '13 at 15:51
  • $\begingroup$ @xslittlegrass, you always can get the content of the Piecewise by using List @@ sol[t] $\endgroup$ – user21 Jul 10 '13 at 15:57
  • $\begingroup$ 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. $\endgroup$ – xslittlegrass Jul 10 '13 at 16:03
  • $\begingroup$ @xslittlegrass, look at the CompilePrint@f3 to see the difference. $\endgroup$ – user21 Jul 10 '13 at 16:08
  • $\begingroup$ 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? $\endgroup$ – xslittlegrass Jul 10 '13 at 21:51
4
$\begingroup$

I know I'm late, but was too curious. The problem really stems from _InterpolationFunction in

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

from the question not matching s[t]. First, I'm not sure if it's a typo or InterpolationFunction has been renamed as InterpolatingFunction in Mathematica 10. Second, _InterpolatingFunction only matches the head of s[t]. What is needed is really _InterpolatingFunction[_] or InterpolatingFunction[__][_] to specify that s[t] is complex. Using this pattern in Compile instead of _InterpolatingFunction yields a compiled function which correctly accepts a real and returns a complex.

$\endgroup$
  • 1
    $\begingroup$ Good catch (+1)! It was indeed a typo since InterpolatingFunction has this name starting from version 2 where it was introduced. The form _InterpolatingFunction[_] works nicely, tested with versions 5.2, 8.0.4 and 10.4.1. It's a shame that no one noticed this before... $\endgroup$ – Alexey Popkov Jul 6 '16 at 11:49
  • 1
    $\begingroup$ I've found this method doesn't work for functions the compiler can create instructions which do not call the Wolfram Language evaluator, such as Sqrt and Power. For such functions, we would need to force a call to the Wolfram Language evaluator, by e.g. Compile[{{t, _Real}}, sqrt[-t], {{_sqrt, _Complex}}] /. sqrt -> Sqrt; $\endgroup$ – obsolesced Jul 7 '16 at 6:01
  • $\begingroup$ I'm not an expert in Compile. I'm surprised that simple f = Compile[{{t, _Complex}}, Sqrt[-t]]; f[-4] returns 2. + 0. I instead of 2 I. I advise you to create a question on these issues. $\endgroup$ – Alexey Popkov Jul 7 '16 at 8:37
  • $\begingroup$ @AlexeyPopkov as you have defined, f[-4] should rightfully evaluate to 2 because of the minus sign in Sqrt[-t]. f[4] does give 0. + 2. I. $\endgroup$ – obsolesced Jul 7 '16 at 15:29
  • 1
    $\begingroup$ Apparently the OP posted another question in which the same mistake was found. $\endgroup$ – obsolesced Jul 7 '16 at 16:44
1
$\begingroup$

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
$\endgroup$

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