Note: this is a clarification of 2 earlier unanswered questions, Compile issues, scoping and order of evaluation and https://mathematica.stackexchange.com/questions/60549/optimising-a-simple-computation-with-vectors .
The problem
I have a list of expressions, e.g.
expressions = {a b + x^7 a^2 + x^30 Ma^2, a b c - 2 Mb^2 x + Mb Mc, a b^2 - Mb Md - x^10 d^2};
where x is a known fixed vector, M is a known fixed matrix, Ma stands for M.a, etc. and the rest are unknown vectors.
x = {1, 2};
M = {{0, 1}, {1, 0}};
The first is a function of a,b , the second of a,b,c and the last of a,b,d.
Note: this is a minimal example, in the actual problem each of these expressions is a list of expressions like it (of fixed length), they are longer and there's more of them, and they need to be evaluated thousands to millions of times.
A code to evaluate these quickly is something like
eval1 = Compile[{{a, _Real, 1}, {b, _Real, 1}},
Module[{Ma = M.a}, a b + x^7 a^2 + x^30 Ma^2], CompilationTarget -> "C"];
eval2 = Compile[{{a, _Real, 1}, {b, _Real, 1}, {c, _Real, 1}},
Module[{Ma = M.a, Mb = M.b, Mc = M.c}, a b c - 2 Mb^2 x + Mb Mc],
CompilationTarget -> "C"];
What I want is to automate the writing of these functions, i.e. I want to write something like
Do[eval[i] = Compile[variables[[i]], Module[Mvariables[[i]], expressions[[i]]],
CompilationTarget -> "C"],{i, 1, Length[expressions]}]
where variables is the argument list, Mvariables is the assignment of M.a to Ma etc.
It is also important that all that can be evaluated in expression beforehand, i.e. without knowing the arguments, is evaluated once and for all and not every time. (for example the x^30 shouldn't be computed every time the function is evaluate but only when defining the function)
Partial solution
variables = Union@Cases[#, x_ /; Not[FreeQ[{a, b, c, d}, x]], Infinity] & /@ expressions
Do[exprs[i] = expressions[[i]], {i, 1, Length[expressions]}]
loopCode[i_] := Compile[Evaluate@({#, _Real, 1} & /@ variables[[i]]),
Module[{Ma = M.a, expr = exprs[i]}, expr], CompilationTarget -> "C",
CompilationOptions -> "InlineExternalDefinitions" -> True]
Do[finalCode[i] = loopCode[i], {i, 1, Length[expressions]}]
This gives,
finalCode[1][{1, 2}, {2, 2}]
{{7., 9.}, {36., 20.}}
This solution has 2 problems.
- I cheated with the initialisation in module.
I want it to figure out for itself which ones are needed. Doing something like
Mvariables = Union@Cases[#, x_ /; Not[FreeQ[{Ma, Mb, Mc, Md}, x]], Infinity] &/@expressions;
Mvariables2 = Map[ToExpression[StringTake[ToString[#], -1]] &, Mvariables, {2}];
Module[Evaluate@(Mvariables[[i]]=(M*#&)/@variables[[i]]),...]
doesn't work as it still isn't evaluated.
- The expression (
exprs[i]) isn't evaluated as much as possible.
To see this we can check:
M = IdentityMatrix[30];
x= 1/Range[30]//N;
c1=x^7;
c2=x^30;
newCode=Compile[{{a,_Real,1},{b,_Real,1}},
Module[{Ma=M.a}, a b + c1 a^2 + c2 Ma^2, CompilationTarget -> "C",
CompilationOptions -> "InlineExternalDefinitions" -> True];
atest = RandomReal[{0, 1}, 30];
btest = RandomReal[{0, 1}, 30];
Now timing,
Do[finalCode[1][atest, btest], {10^5}] // Timing
{0.565053, Null}
Do[newCode[atest, btest], {10^5}] // Timing
{0.365492, Null}
I think both issues have to do with evaluation inside module, but I don't see how to fix it.
