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], _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},

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.


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.

  • $\begingroup$ If I copy your code as given, then trying to compile the very first function fails with an error "Sequence specification...". Can you please confirm that you have all the code necessary to reproduce the problem? $\endgroup$ – dr.blochwave Sep 15 '14 at 8:57
  • $\begingroup$ Yes that's right, what I meant was that if you copy the first code, and then by hand replace the i with 1 or 2 or 3, then it does work, but if you do these replacements in a loop, it doesn't work. $\endgroup$ – Jansen Sep 15 '14 at 9:10
  • $\begingroup$ I found a way to do it in a loop just now. If in the last loopComp you replace Module[{da=mat.a},exprIndexed[i]] with Module[{da=mat.a,expr=exprIndexed[i]},expr], it works. I would still like to understand why though. $\endgroup$ – Jansen Sep 15 '14 at 9:12

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.