# Evaluating arguments of module (inside compile)

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, 2}, {2, 2}]
{{7., 9.}, {36., 20.}}


## This solution has 2 problems.

1. 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.

1. The expression (exprs[i]) isn't evaluated as much as possible.

To see this we can check:

M = IdentityMatrix;
x= 1/Range//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[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.

I suggest that you use local rules for your fixed values. Then, the following seems to do what you ask for:

ClearAll[dotQ, init, vars];
dotQ[s_String, MSymbol_Symbol] :=
StringMatchQ[s, ToString[MSymbol] ~~ LetterCharacter];

vars[MSymbol_] := Select[Not@*(dotQ[#, MSymbol] &)@*ToString]@*Variables;

init[expr_, MSymbol_Symbol] :=
Thread[#, Hold] & @ Union @ Cases[
expr,
s_Symbol :>
With[{eval = ToString[s]},
ToExpression[
StringDrop[eval, 1],
StandardForm,
Function[sym, Hold[s = MSymbol.sym], HoldAll]
] /; dotQ[eval, MSymbol]
],
Infinity
];

ClearAll[compile];
compile[expr_, MSymbol_, rules_] :=
With[{
i = init[expr, MSymbol] /. rules,
computed =
Hold[
Evaluate[expr /. {Plus -> Hold[Plus], Times -> Hold[Times]} /. rules]
] /. {Hold[Plus] -> Plus, Hold[Times] -> Times}
},
Prepend[
i /. Hold[x_] :> (computed /. Hold[y_] :> Hold[Hold[Module][x, y]]) /. Hold[Module] -> Module,
{#, _Real, 1} & /@ vars[MSymbol][expr /. rules]
] /. Hold[args__] :> Compile[args, "CompilationTarget" -> "MVM"]
];


You use this as:

ClearAll[M,a,b,c,d,x];
eval = Map[compile[#, M, {x -> {1, 2}, M -> {{0, 1}, {1, 0}}}] &]@expressions


Some examples to illustrate what the functions do:

vars[M]@expressions[]

(* {a, b, x} *)

init[expressions[], M]

(* Hold[{Mb = M.b, Mc = M.c}] *)


This could've been done in a less cryptic manner without using the tricks I use, but that would require some more code.

Note that, to my mind, your choice of using symbols like Ma etc needlessly complicate things. Note also that I am using operator forms in a few places, so this would have to be modified to work on V9 and earlier.

• Thanks, I have to study this to see how it works, but for me it doesn't give the result I want. When i apply the first function to the arguments {1,2},{3,4} I get a 2x2 matrix, whereas it should be a vector. Maybe I should have mentioned this: the multiplication in the expressions should be interpreted as element wise multiplication of lists. Sep 25, 2014 at 16:26
• And the reason that I'm using these Ma variables is that in the actual expressions they occur multiple times, so this limits the number of matrix multiplications to the minimum. Sep 25, 2014 at 16:37
• @Jansen Ok, the problem there was that on one hand, you want it to be partially evaluated, on another hand, Plus and Times are Listable, so one has to fight that during the partial computation of your expression. I didn't have the time to make the code pretty, but that part has now been corrected. At this point, it becomes clear that your problem is complex enough that a custom expression parser / code generator would be much cleaner than what I wrote. No time for this at the moment though. But I do realize that the code above can be hard to digest. Sep 25, 2014 at 16:53
• Thanks, it works exactly as I wanted. Don't know what you mean by custom expression parser/code generator though. Btw, also thanks for writing the book on mathematica, learned a lot from it :) Sep 26, 2014 at 8:30
• @Jansen Glad I could help, thanks for the accept. Re: book - good to know that it is useful! Oct 2, 2014 at 9:07