I have a few questions about compiling functions, which I think are all related to scoping and order of evaluation.
I will illustrate them by a minimal example of the problem I have. I'm sorry for the lengthy question.
The setup code is this,
mat = {{0, 1}, {1, 0}};
expr = {a^2 - ma + 3 x^3, a mb^2 + 2 b x, mc a - mb ma x^6};
Do[exprIndexed[i] = expr[[i]], {i, 1, 3}]
x = {5, -2};
atest = {1, 2}; btest = {2, 1}; ctest = {0, 0};
I have a list of expressions expr, each of which represents a vector.
The first depends only on a
, ma
stands for the matrix mat
times the vector a
, and x
is a constant known vector. The second depends on a
and b
and the last on a, b, c
.
I want to write for each of these expressions a compiled code that evaluates them. In particular I want to be able to define these compiled functions inside a do loop, I want them to really fully compile, and I want it to evaluate as much as it can without knowing the variables (a
, {a, b}
and {a, b, c}
respectively).
One attempt is this:
(NOTE: i'm omitting everywhere the last argument CompilationOptions -> {"InlineExternalDefinitions" -> True}, CompilationTarget -> "C"
for brevity)
Compile[Evaluate[({#, _Real, 1} &) /@ Take[{a, b, c}, i]],
ma = mat.a;
If[i > 1, mb = mat.b];
If[i > 2, dc = mat.c];
exprIndexed[i],
{{exprIndexed[i], _Real, 1}, {ma, _Real, 1}, {mb, _Real,1}, {mat, _Real, 2}}]
(A not so important question I have here is: when it is and when it's not necessary to include these type specifications at the end)
Plugging in by hand i = 1
, and copying for i = 2
, etc., this code works.
Wrapping it in a Do
, e.g.
Do[comp[i]= (expr above)]
does not work. Mathematica gives
CompiledFunction::cfse: Compiled expression {373. +a^2,-25.+a^2} should be a machine-size integer. >>
However, writing
loopCode[i_] := (code above)
Do[comp[i] = loopCode[i]]
does work. I wonder why this is so?
Although the code works, it does not compile completely (I checked this by doing Export["test.c", comp[1]]
). The reason is that the assignments set the global value of e.g. da, so must talk to Mathematica.
This can be fixed by changing the body of the compile to
Module[{ma = mat.a},
exprIndexed[i]]
Now it works and compiles, but it cannot be defined in a loop.
Explicitly, I tried,
loopComp[i_] :=
Compile[Evaluate[({#, _Real, 1} &) /@ Take[{a, b, c}, i]],
Module[{da = mat.a}, exprIndexed[i]],
{{exprIndexed[i], _Real, 1}, {da, _Real, 1}, {db, _Real, 1}, {mat, _Real, 2}}]
Do[comp[i] = loopComp[i], {i, 1, 3}]
I haven't found a way to automate the assignment of ma
, mb
and mc
, i.e. to get Module[{ma=mat.a},...]
for i=1
, Module[{ma=mat.a, mb=mat.b},...}]
for i=2
etc.,
so comp[2]
and comp[3]
shouldn't work, but at least comp[1]
should, but it does not.
comp[1][atest]
gives the error
CompiledFunction::cfte: Compiled expression da should be a rank 1 tensor of machine-size real numbers. >>
and the output
{376 - da, -20 - da}
I think the issue is in the scoping of Do
, Module
and Compile
, but I don't know how to fix it.
The final issue is this. Even in the working versions of this code (just by explicitly copy-pasting and inserting the i = 1, 2, 3
) each time it runs everything is computed, i.e. the x^6
is not computed by compile once and for all, but every time the code is run.
Note that if x
is defined before expr is defined, then this does not happen, but then there is another problem in that i.e. x + a
will become {5 + a, -2 + a}
and when a
is plugged in we get a matrix instead of a vector.
So I want to evaluate as much as possible, while still getting a vector as output.