# How to force evaluation of an expression with Part left unnevaluated?

My goal is to minimize a real valued multi-variable function averageFreeEnergyDensity. The function accepts an integer argument n, a set of real parameters j,d,a,h,dx and two real matrices th,ph. I want to find matrices th,ph which minimize the function. FindMinimum minimizes the function well enough but it is too slow at the moment. I am trying to compile the function to accelerate FindMinimum. If I compile the function for arbitrary n it is at least an order of magnitude more slow than if I make a set of functions each of which is defined for a particular n. My current fastest solution is


averageFreeEnergyDensityCompiled["square",n_]:=
averageFreeEnergyDensityCompiled["square", n]=
Compile[{{j, _Real}, {d, _Real}, {a, _Real}, {h, _Real}, {dx,
_Real}, {th, _Real, 2}, {ph, _Real, 2}}, Evaluate[
(*apply boundary conditions*)
s[i1_, i2_] := Which[
i1 == n + 1 && i2 == n + 1, s[1, 1],
i1 == n + 1, s[1, i2],
i2 == n + 1, s[i1, 1],
True, {Cos[ph[[i1, i2]]] Sin[th[[i1, i2]]],
Sin[ph[[i1, i2]]] Sin[th[[i1, i2]]], Cos[th[[i1, i2]]]}
];
(*define finite difference functions*)
dsdx[i1_, i2_] := (s[i1 + 1, i2] - s[i1, i2])/dx;
dsdy[i1_, i2_] := (s[i1, i2 + 1] - s[i1, i2])/dx;
(*return average free energy density function*)
Sum[j (dsdx[i1, i2].dsdx[i1, i2] + dsdy[i1, i2].dsdy[i1, i2]) +
a s[i1, i2][]^2 - h s[i1, i2][], {i1, 1, n}, {i2, 1, n}]/
n^2
]]

I need to delay evaluation because the expression being compiled can only be evaluated for particular n. I store the function because I will evaluate it for many values of j,d,a,h,dx. This solution is fast but it throws an error the first time the function is evaluated

Block[{n = 100},
AbsoluteTiming[
averageFreeEnergyDensityCompiled["square", n][1, 1, 1, 0, 0.1,
Array[0 &, {n, n}], RandomReal[2 Pi, {n, n}]]]]
Part::partd: Part specification ph$[[2,1]] is longer than depth of object. >> Part::partd: Part specification th$[[2,1]] is longer than depth of object. >>
Part::partd: Part specification ph\$[[2,1]] is longer than depth of object. >>
General::stop: Further output of Part::partd will be suppressed during this calculation. >>

After several attempts to evaluate the Part function Mathematica gives up and evaluates the expression ignoring Part. The final expression can then be compiled. I want to manually tell Mathematica to ignore Part. I have tried using Inactive,Activate but this slows down the final compiled function by orders of magnitude. I would ideally give an option to Evaluate like Evaluate["insert expression here",Holding->Part] but Evaluate has no options that I know of.

Note: The boundary condition I have used here is simple and can be taken out of compile, but I will be using more complicated boundary conditions once I solve this problem. I do not need to use finite difference steps of order greater than 1. I will also be using this code on a triangular grid which will further change the finite difference.

• Welcome to Mathematica.SE! I hope you will become a regular contributor. To get started, 1) take the introductory Tour now, 2) when you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge, 3) remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign, and 4) give help too, by answering questions in your areas of expertise. – bbgodfrey May 23 '15 at 1:05

I believe Block will work for you here:

mem : averageFreeEnergyDensityCompiled["square", n_] := mem =
Block[{Part},
Compile[{{j, _Real}, {d, _Real}, {a, _Real}, {h, _Real}, {dx, _Real}, {th, _Real,
2}, {ph, _Real, 2}},
Evaluate[(*apply boundary conditions*)
s[i1_, i2_] :=
Which[i1 == n + 1 && i2 == n + 1, s[1, 1], i1 == n + 1, s[1, i2], i2 == n + 1,
s[i1, 1], True, {Cos[ph[[i1, i2]]] Sin[th[[i1, i2]]],
Sin[ph[[i1, i2]]] Sin[th[[i1, i2]]], Cos[th[[i1, i2]]]}];
(*define finite difference functions*)
dsdx[i1_, i2_] := (s[i1 + 1, i2] - s[i1, i2])/dx;
dsdy[i1_, i2_] := (s[i1, i2 + 1] - s[i1, i2])/dx;
(*return average free energy density function*)
Sum[j (dsdx[i1, i2].dsdx[i1, i2] + dsdy[i1, i2].dsdy[i1, i2]) + a s[i1, i2][]^2 -
h s[i1, i2][], {i1, 1, n}, {i2, 1, n}]/n^2]]
]


However I believe CompileGetElement will also work and perform better too. Search for "GetElement" (on this site) for more information.

mem : averageFreeEnergyDensityCompiled["square", n_] := mem =
With[{Part = CompileGetElement},
Compile[{{j, _Real}, {d, _Real}, {a, _Real}, {h, _Real}, {dx, _Real}, {th, _Real,
2}, {ph, _Real, 2}},
Evaluate[(*apply boundary conditions*)
s[i1_, i2_] :=
Which[i1 == n + 1 && i2 == n + 1, s[1, 1], i1 == n + 1, s[1, i2], i2 == n + 1,
s[i1, 1], True, {Cos[ph[[i1, i2]]] Sin[th[[i1, i2]]],
Sin[ph[[i1, i2]]] Sin[th[[i1, i2]]], Cos[th[[i1, i2]]]}];
(*define finite difference functions*)
dsdx[i1_, i2_] := (s[i1 + 1, i2] - s[i1, i2])/dx;
dsdy[i1_, i2_] := (s[i1, i2 + 1] - s[i1, i2])/dx;
(*return average free energy density function*)
Sum[j (dsdx[i1, i2].dsdx[i1, i2] + dsdy[i1, i2].dsdy[i1, i2]) + a s[i1, i2][]^2 -
h s[i1, i2][], {i1, 1, n}, {i2, 1, n}]/n^2]]
]


(See the note at the bottom of this post for explanation of the mem thing.)

• I see no speed improvements with Compile`GetElement but it is more aesthetically pleasing for me. There is still an error when using compiled functions in FindMinimum; it attempts to evaluate the function symbolically. The solution to this problem is to define a function which is only defined for numeric arguments. – James Rowland May 26 '15 at 19:52