# Compile with array operations

When I use the function f with only the Module

f[a1_, a2_, a3_]:= Module[.....]
f[0,0,1200]


works fine but gives errors when I use compile as follows.

demand[n_, k_, qr_] := Min[k Vf, n capacity - qr]; supply[n_, k_, qr_] := Min[(n Kj - k) w, n capacity - qr];
Nsupply[n_, k_, qsum_, qrsum_, qr_] := Min[(n Kj - k) Vf - qsum - qrsum, n capacity - qr];
floExct[n_, ku_, k_, qsum_, qrsum_, qr_] := Min[demand[n, ku, qr], Nsupply[n, k, qsum, qrsum, qr]];
gamma[ku_, kd_] := Min[1, supply[L, kd, 0]/(demand[L, ku, 0] + 0.001)];
inflow[phi_, qin_] := (phi - \[Beta] qin) dx;
density[qin_, qout_, qr_] := (qin - qout + qr)/Vf;

Kj = 150.; w = 20.; Vf = 100.; capacity = 2500.; n = 30; m = 100; p = 27; RML = 18; \[Alpha] = 1200.; \[Beta] = 0.; L = 1.; delta = 1.; dt =4./3600.; dx = 1./9.;

f= Compile[{{a1, _Integer}, {a2, _Integer}, {a3, _Integer}},
Module[{k0 = ConstantArray[0., n],
Fk = Table[Table[0., {i1, n}], {j, 5}],
kr = Table[Table[0., {i1, 1, p}], {i2, 1, n - 2}],
Rk = Table[Table[Table[0., {i1, 1, p}], {i2, 1, n - 2}], {i3, 4}], \[Phi],
Fq = Table[Table[0., {i1, n - 1}], {j, 4}], qin,
Fin = Table[Table[0., {i1, n - 2}], {j, 4}], qr,
Rq = Table[Table[Table[0., {i1, 1, p}], {i2, 1, n - 2}], {i3, 4}],
RMori = Table[0., {i1, 1, n - 2}],
RM = Table[0., {i1, 1, n - 2}], TT = 0., qf, qsum, qrsum},
RMori = Table[a1 (i1 dx)^2 + a2 (i1 dx) + a3, {i1, 1, n - 2}];
kr = ReplacePart[#, (\[Alpha] delta/Vf), 1] & /@ kr; AppendTo[Rk, kr];

\[Phi] = demand[1, #[[-1]], 0] & /@ kr MapThread[gamma, {Take[k0, {2, -2}], Take[k0, {3, -1}]}];
qsum = N[Plus @@ Fq];
qrsum = Join[Plus @@ Fin, {0}];
qf = MapThread[floExct[L, #1, #2, #3, #4, #5] &, {Most@k0, Rest@(Fk[[1]]), qsum,qrsum, Join[\[Phi], {0}]}];
qin = inflow[\[Phi], Most@qf];
k0 += Join[{0}, density[Most@qf, Rest@qf, qin], {0}];
qsum = MapThread[Plus, #] &@Rq;
RM = MapThread[Min[#1, #2] &, {RMori, MapThread[floExct[1, #1[[RML]], #2[[RML + 1]], #3[[RML]], 0, 0] &, {kr,Rk[[1]], qsum}]}];
qr = MapThread[Join[#1, {#2}, #3, {#4}] &, {
MapThread[floExct[1, #1, #2, #3, 0, 0] &, {Take[#, {1, RML - 1}] & /@ kr,Take[#, {2, RML}] & /@ Rk[[1]], Take[#, {1, RML - 1}] & /@ qsum}, 2], RM,
MapThread[floExct[1, #1, #2, #3, 0, 0] &, {Take[#, {RML + 1, -2}] & /@ kr, Take[#, {RML + 2, -1}] & /@ Rk[[1]], Take[#, {RML + 1, -2}] & /@ qsum}, 2], \[Phi]}];
kr = MapThread[(#1 + #2) &, {Join[{0}, #] & /@ (density[Most@#, Rest@#, 0] & /@ qr), kr}];

TT += Plus @@ k0; TT += Plus @@ (Plus @@ kr);

Fq = Delete[Fq, {1}]; Fin = Delete[Fin, {1}]; Rq = Delete[Rq, {1}]; Fk = Delete[Fk, {1}]; Rk = Delete[Rk, {1}];
AppendTo[Fq, qf]; AppendTo[Fin, qin]; AppendTo[Rq, qr]; AppendTo[Fk, k0]; AppendTo[Rk, kr];
TT]]


all the errors (as shown below) are related to the qsum variable. Can someone tell me what I am missing?

MapThread::list: List expected at position 2 in MapThread[Plus,CompileFunctionVariable$266]. >> Compile::cprank: Compile cannot determine the rank of the result tensor. The length of CompileFunctionVariable$266 in MapThread[Plus,CompileFunctionVariable$266] is unknown at compiling time; evaluation will use the uncompiled function. >> Compile::cset: Variable qsum of type {_Real,1} encountered in assignment of type {_Integer,0}. >> Compile::part: Part specification CompileGetElement[qsum,CompileVariable$912][[18]] cannot be compiled since the argument is not a tensor of sufficient rank. Evaluation will use the uncompiled function. >>
Compile::part: Part specification CompileGetElement[qsum,CompileVariable$912][[18]] cannot be compiled since the argument is not a tensor of sufficient rank. Evaluation will use the uncompiled function. >> Compile::part: Part specification CompileGetElement[qsum,CompileVariable$1852][[18]] cannot be compiled since the argument is not a tensor of sufficient rank. Evaluation will use the uncompiled function. >>
General::stop: Further output of Compile::part will be suppressed during this calculation. >>
Compile::cplist: CompileGetElement[qsum,SystemPrivateCompileSymbol[0]] should be a tensor of type Integer, Real, or Complex; evaluation will use the uncompiled function. >>
Compile::cplist: CompileGetElement[qsum,SystemPrivateCompileSymbol[0]] should be a tensor of type Integer, Real, or Complex; evaluation will use the uncompiled function. >>
Compile::cplist: CompileGetElement[qsum,SystemPrivateCompileSymbol[0]] should be a tensor of type Integer, Real, or Complex; evaluation will use the uncompiled function. >>
General::stop: Further output of Compile::cplist will be suppressed during this calculation. >>
MapThread::list: List expected at position 2 in MapThread[Plus,CompileFunctionVariable$3193]. >> Compile::cprank: Compile cannot determine the rank of the result tensor. The length of CompileFunctionVariable$3193 in MapThread[Plus,CompileFunctionVariable\$3193] is unknown at compiling time; evaluation will use the uncompiled function. >>
Compile::cset: Variable qsum of type {_Real,1} encountered in assignment of type {_Integer,0}. >>


Edit:

Please note that the above code is a simplified/partial piece of a larger code. I presented this to show the problem area. The original code has several hundred iterations (till a condition is met) of the above critical steps and the compiled function itself will be called several thousands of times (for direct search-based optimization).

Edit 2:

Thanks to xzczd, the following code compiles without any errors. However, it gives error/warnings during run-time. The code may look cumbersome, but the problem area is very specific.

f = With[{Kj = 150., w = 20., Vf = 100., capacity = 2500., n = 15, m = 30, p = 10, RML = 5, \[Alpha] = 1200., \[Beta] = 0., L = 1., delta = 1., dt = 4./3600., dx = 1./9.},
Module[{demand, supply, Nsupply, floExct, gamma, inflow, density}, demand = Min[#2 Vf, # capacity - #3] &; supply = Min[(# Kj - #2) w, # capacity - #3] &; Nsupply = Min[(# Kj - #2) Vf - #3 - #4, # capacity - #5] &;
floExct = Evaluate@Min[demand[#, #2, #6], Nsupply[#, #3, #4, #5, #6]] &; gamma = Evaluate@Min[1, supply[L, #2, 0]/(demand[L, #, 0] + 0.001)] &; inflow = (# - \[Beta] #2) dx &; density = (# - #2 + #3)/Vf &;
Compile[{{a1, _Integer}, {a2, _Integer}, {a3, _Integer}},
Module[{k0 = Table[0., {n}], Fk = Table[0., {5}, {n}], kr = Table[0., {n - 2}, {p}], Rk = Table[0., {4}, {n - 2}, {p}], Fq = Table[0., {4}, {n - 1}], Fin = Table[0., {4}, {n - 2}], Rq = Table[0., {4}, {n - 2}, {p}],
RMori = Table[a1 (i1 dx)^2 + a2 (i1 dx) + a3, {i1, 1, n - 2}], Shutoff = False, j = 0, \[Phi], qin, qr, qf, qsum, qrsum, RM, TT, NtwrkTT},
If[300 > Min[RMori] && Max[RMori] > 1200, Print["MR not within bounds"]; NtwrkTT = 10^20,
kr = ReplacePart[#, (\[Alpha] delta/Vf), 1] & /@ kr;
NtwrkTT = TT = Total@(Total@kr); AppendTo[Rk, kr];
While[TT > 0, TT = 0;
\[Phi] = demand[1, #[[-1]], 0] & /@ kr MapThread[gamma, {Take[k0, {2, -2}], Take[k0, {3, -1}]}];
qrsum = Join[Total@Fin, {0}];
qf = MapThread[floExct[L, #1, #2, #3, #4, #5] &, {Most@k0, Rest@(Fk[[1]]), Total@Fq, qrsum, Join[\[Phi], {0}]}];
qin = inflow[\[Phi], Most@qf];
k0 += Join[{0}, density[Most@qf, Rest@qf, qin], {0}];
qsum = Total@Rq;
RM = MapThread[Min[#1, #2] &, {RMori, MapThread[floExct[1, #1[[RML]], #2[[RML + 1]], #3[[RML]], 0, 0] &, {kr, Rk[[1]], qsum}]}];
qr = MapThread[Join[#1, {#2}, #3, {#4}] &, {MapThread[floExct[1, #1, #2, #3, 0, 0] &, {Take[#, {1, RML - 1}] & /@ kr, Take[#, {2, RML}] & /@ Rk[[1]], Take[#, {1, RML - 1}] & /@ qsum}, 2], RM,
MapThread[floExct[1, #1, #2, #3, 0, 0] &, {Take[#, {RML + 1, -2}] & /@ kr, Take[#, {RML + 2, -1}] & /@ Rk[[1]], Take[#, {RML + 1, -2}] & /@ qsum}, 2], \[Phi]}];
kr = MapThread[(#1 + #2) &, {Join[{0}, #] & /@ (density[Most@#, Rest@#, 0] & /@ qr), kr}];
TT += Total@k0; TT += Total@(Total@kr); NtwrkTT += TT; j++;
Fq = Delete[Fq, {1}]; Fin = Delete[Fin, {1}]; Rq = Delete[Rq, {1}]; Fk = Delete[Fk, {1}]; Rk = Delete[Rk, {1}];
AppendTo[Fq, qf]; AppendTo[Fin, qin]; AppendTo[Rq, qr]; AppendTo[Fk, k0]; AppendTo[Rk, kr];
Print["j: ", j];
If[#[[2]] > Kj - \[Alpha]/w, Print["Ramp Spillback"]; NtwrkTT = 10^20; Break[]] & /@ kr;
If[(Shutoff == False) && j > m, Shutoff = True; kr = ReplacePart[#, 0, 1] & /@ kr];
];
]; NtwrkTT dt], CompilationOptions -> {"InlineExternalDefinitions" -> True}]]];


when I run the compiled function with f[0,0,1100], it runs the While loop only once, gives the following error:

CompiledFunction::cfse: Compiled expression {Null,Null,Null,Null,Null,Null,Null,Null,Null,Null,Null,Null,Null} should be a machine-size real number. >>
CompiledFunction::cfex: Could not complete external evaluation at instruction 491; proceeding with uncompiled evaluation. >>


Then, it proceeds to run the program, I guess in an uncompiled fashion.

One possible variable that refers to {Null,Null,Null,Null,Null,Null,Null,Null,Null,Null,Null,Null,Null} may be RM and I was wondering, if it is the source of problem.

-
Can you make some effort on reproducing your problem in a minimal code sample? BTW, I doubt even if you fix your code, it'll be slow since it has made so many external calls in Compile. – xzczd Apr 15 '14 at 3:19
Several observations: 1. qsum = MapThread[Plus, #] &@Rq can be changed to qsum=Total@Rq, 2. qsum = N[Plus @@ Fq] can be embedded into the definition of qf, 3. kr = ReplacePart[#, (\[Alpha] delta/Vf), 1] & /@ kr can be changed to kr[[All, 1]] = \[Alpha] delta/Vf. There still exists warnings but I'd like to stop here, as I mentioned above, I don't think compiling these code makes much sense. – xzczd Apr 15 '14 at 4:48
@xzczd Thank you for your comments and suggestions. As I have mentioned in the edit section above, what I presented is only a simplified/partial code. Are you suggesting that I move the external functions to inside compile or compile the external functions (demand,supply,Nsupply,floExct..etc) and call them from the main compile function f? – brama Apr 15 '14 at 10:33
Precisely speaking, I suggest you to at least avoid function definition based on pattern-matching. (see this post for details. ) – xzczd Apr 15 '14 at 11:01

OK, as an after-meal exercise, I fixed your code. The main modifications I've done are:

1. Based on trick mentioned in this post, change your pattern-matching function into pure function.

2. Change qsum = MapThread[Plus, #] &@Rq(for some unclear reason it can't be compiled, maybe it's because Rq isn't a "explicit" list? ) into qsum=Total@Rq.

3. Introduce constants with With to avoid external calls in Compile caused by unlocalized variables. (Explanation for this involves explanation for how With works, which I'd like to omit here. There already exists some good posts for this issue in this site, for example this. )

4. Embed qsum = N[Plus @@ Fq] into the definition of qf to avoid error (This seems to be because the dimension of qsum will be changed in latter code if we don't eliminate it here ).

Note: I've discard my suggestion in the comment above of using kr[[All, 1]] = α delta/Vf instead of kr = ReplacePart[#, (α delta/Vf), 1] & /@ kr because I found the latter can be compiled with the modifications mentioned above, what's more, kr[[All, 1]] = α delta/Vf can't be compiled, kr[[All, 1]] = Table[α delta/Vf, {n - 2}] can, though.

The following will generate no error and be about 20 times faster than the corresponding uncompiled version:

f = With[{Kj = 150., w = 20., Vf = 100., capacity = 2500., n = 30,
m = 100, p = 27, RML = 18, α = 1200., β = 0., L = 1.,
delta = 1., dt = 4./3600., dx = 1./9.},
Module[{demand, supply, Nsupply, floExct, gamma, inflow, density},
demand = Min[#2 Vf, # capacity - #3] &;
supply = Min[(# Kj - #2) w, # capacity - #3] &;
Nsupply = Min[(# Kj - #2) Vf - #3 - #4, # capacity - #5] &;
floExct = Evaluate@Min[demand[#, #2, #5], Nsupply[#, #3, #4, #5, #6]] &;
gamma = Evaluate@Min[1, supply[L, #2, 0]/(demand[L, #, 0] + 0.001)] &;
inflow = (# - β #2) dx &;
density = (# - #2 + #3)/Vf &;
Compile[{{a1, _Integer}, {a2, _Integer}, {a3, _Integer}},
Module[{k0 = Table[0., {n}],
Fk = Table[0., {5}, {n}],
kr = Table[0., {n - 2}, {p}],
Rk = Table[0., {4}, {n - 2}, {p}], ϕ,
Fq = Table[0., {4}, {n - 1}], qin,
Fin = Table[0., {4}, {n - 2}], qr,
Rq = Table[0., {4}, {n - 2}, {p}],
RMori = Table[a1 (i1 dx)^2 + a2 (i1 dx) + a3, {i1, 1, n - 2}],
RM = Table[0., {n - 2}], TT = 0., qf, qsum, qrsum},

kr = ReplacePart[#, (α delta/Vf), 1] & /@ kr;
AppendTo[Rk, kr];
ϕ =
demand[1, #[[-1]], 0] & /@ kr MapThread[
gamma, {Take[k0, {2, -2}], Take[k0, {3, -1}]}];
qrsum = Join[Plus @@ Fin, {0}];
qf =
floExct[L, #1, #2, #3, #4, #5] &, {Most@k0, Rest@(Fk[[1]]),
Plus @@ Fq, qrsum, Join[ϕ, {0}]}];
qin = inflow[ϕ, Most@qf];
k0 += Join[{0}, density[Most@qf, Rest@qf, qin], {0}];
qsum = Total@Rq;
RM =
Min[#1, #2] &, {RMori,
floExct[1, #1[[RML]], #2[[RML + 1]], #3[[RML]], 0,
0] &, {kr, Rk[[1]], qsum}]}];
qr =
Join[#1, {#2}, #3, {#4}] &, {MapThread[
floExct[1, #1, #2, #3, 0, 0] &, {Take[#, {1, RML - 1}] & /@
kr, Take[#, {2, RML}] & /@ Rk[[1]],
Take[#, {1, RML - 1}] & /@ qsum}, 2], RM,
floExct[1, #1, #2, #3, 0, 0] &, {Take[#, {RML + 1, -2}] & /@
kr, Take[#, {RML + 2, -1}] & /@ Rk[[1]],
Take[#, {RML + 1, -2}] & /@ qsum}, 2], ϕ}];
kr =
MapThread[(#1 + #2) &, {Join[{0}, #] & /@ (density[Most@#,
Rest@#, 0] & /@ qr), kr}];
TT += Plus @@ k0; TT += Plus @@ (Plus @@ kr);
Fq = Delete[Fq, {1}]; Fin = Delete[Fin, {1}];
Rq = Delete[Rq, {1}]; Fk = Delete[Fk, {1}]; Rk = Delete[Rk, {1}];
AppendTo[Fq, qf]; AppendTo[Fin, qin]; AppendTo[Rq, qr];
AppendTo[Fk, k0]; AppendTo[Rk, kr];
TT],
CompilationOptions -> {"InlineExternalDefinitions" -> True}]]];


I believe this function can still be optimized because many CopyTensor, which is said to severely influence the performance of compiled function, exists in the output of CompilePrint, but this time I'll really stop and leave it to you. As for the resource of this issue, I think you've already found this post, there exist discussion for CopyTensor, just search CopyTensor inside the page, there exist other posts about the issue in this site too, for example this, you can have a search.

Update:

To fix your final code, I made the following modifications:

1. Change all the NtwrkTT = 10^20 into NtwrkTT = 10^20.. (Though you mentioned the code is compiled without error in your computer, it fails without this modification in my v8.0.4, Windows Vista Home Basic 32bit. )

2. Change

If[#[[2]] > Kj - α/w, Print["Ramp Spillback"]; NtwrkTT = 10^20; Break[]] & /@ kr;


into

Do[If[i[[2]] > Kj - α/w, Print["Ramp Spillback"]; NtwrkTT = 10^20; Break[]], {i, kr}];


Using Map and Apply inside Compile needs corresponding mental preparation: It can easily cause failure. (At one time I even thought that "The only function arguments supported are Times, Plus, or List", but your code has already shown it's not correct.)

The fixed code is as following:

f = With[{Kj = 150., w = 20., Vf = 100., capacity = 2500., n = 15,
m = 30, p = 10, RML = 5, α = 1200., β = 0., L = 1.,
delta = 1., dt = 4./3600., dx = 1./9.},
Module[{demand, supply, Nsupply, floExct, gamma, inflow, density},
demand = Min[#2 Vf, # capacity - #3] &;
supply = Min[(# Kj - #2) w, # capacity - #3] &;
Nsupply = Min[(# Kj - #2) Vf - #3 - #4, # capacity - #5] &;
floExct = Evaluate@Min[demand[#, #2, #6], Nsupply[#, #3, #4, #5, #6]] &;
gamma = Evaluate@Min[1, supply[L, #2, 0]/(demand[L, #, 0] + 0.001)] &;
inflow = (# - β #2) dx &;
density = (# - #2 + #3)/Vf &;
Compile[{{a1, _Integer}, {a2, _Integer}, {a3, _Integer}},
Module[{
k0 = Table[0., {n}],
Fk = Table[0., {5}, {n}],
kr = Table[0., {n - 2}, {p}],
Rk = Table[0., {4}, {n - 2}, {p}],
Fq = Table[0., {4}, {n - 1}],
Fin = Table[0., {4}, {n - 2}],
Rq = Table[0., {4}, {n - 2}, {p}],
RMori = Table[a1 (i1 dx)^2 + a2 (i1 dx) + a3, {i1, 1, n - 2}],
Shutoff = False, j = 0, ϕ, qin, qr, qf, qsum, qrsum, RM, TT, NtwrkTT},
If[300 > Min[RMori] && Max[RMori] > 1200,
Print["MR not within bounds"]; NtwrkTT = 10^20.,
kr = ReplacePart[#, (α delta/Vf), 1] & /@ kr;
NtwrkTT = TT = Total@(Total@kr); AppendTo[Rk, kr];
While[TT > 0, TT = 0;
ϕ =
demand[1, #[[-1]], 0] & /@ kr MapThread[
gamma, {Take[k0, {2, -2}], Take[k0, {3, -1}]}];
qrsum = Join[Total@Fin, {0}];
qf =
floExct[L, #1, #2, #3, #4, #5] &, {Most@k0, Rest@(Fk[[1]]),
Total@Fq, qrsum, Join[ϕ, {0}]}];
qin = inflow[ϕ, Most@qf];
k0 += Join[{0}, density[Most@qf, Rest@qf, qin], {0}];
qsum = Total@Rq;
RM =
Min[#1, #2] &, {RMori,
floExct[1, #1[[RML]], #2[[RML + 1]], #3[[RML]], 0,
0] &, {kr, Rk[[1]], qsum}]}];
qr =
Join[#1, {#2}, #3, {#4}] &, {MapThread[

floExct[1, #1, #2, #3, 0, 0] &, {Take[#, {1, RML - 1}] & /@
kr, Take[#, {2, RML}] & /@ Rk[[1]],
Take[#, {1, RML - 1}] & /@ qsum}, 2], RM,
floExct[1, #1, #2, #3, 0,
0] &, {Take[#, {RML + 1, -2}] & /@ kr,
Take[#, {RML + 2, -1}] & /@ Rk[[1]],
Take[#, {RML + 1, -2}] & /@ qsum}, 2], ϕ}];
kr =
MapThread[(#1 + #2) &, {Join[{0}, #] & /@ (density[Most@#,
Rest@#, 0] & /@ qr), kr}];
TT += Total@k0; TT += Total@(Total@kr); NtwrkTT += TT; j++;
Fq = Delete[Fq, {1}]; Fin = Delete[Fin, {1}];
Rq = Delete[Rq, {1}]; Fk = Delete[Fk, {1}];
Rk = Delete[Rk, {1}];
AppendTo[Fq, qf]; AppendTo[Fin, qin]; AppendTo[Rq, qr];
AppendTo[Fk, k0]; AppendTo[Rk, kr];
Print["j: ", j];
Do[If[i[[2]] > Kj - α/w, Print["Ramp Spillback"];
NtwrkTT = 10^20; Break[]], {i, kr}];
If[(Shutoff == False) && j > m, Shutoff = True;
kr = ReplacePart[#, 0, 1] & /@ kr];];]; NtwrkTT dt],
CompilationOptions -> {"InlineExternalDefinitions" -> True}]]];


At the end of this post, I'd like to say something about Compile. Generally speaking, though I inevitably use Compile in my task, I think it's really a black sheep inside Mathematica, not because its limitation for one's programming ideas, but its elusiveness. Your code is finally speed up by Compile, but I believe a better design can make it faster and more elegant and keep it away from the endless grief looking for the part that isn't compiled. Since your original code has become quite large, this suggestion might be a little late, but I really recommend writing your program in a more… mathematica-like way in your next task.

-
Thank you for walking through this. I am not sure what a copyTensor is. Can you please explain or point me to a resource? Also, As I have posted the final code in the edit2 section. Can you please look at it and let me know what the problem could be? – brama Apr 16 '14 at 17:53
@brama See my edit. (Links for CopyTensor are added above Update.) – xzczd Apr 17 '14 at 4:12
Thank you for your elaborate comments and suggestions. I have started to use Mathematica in the last 6 months and really learnt a lot during the course. This is one of the many pieces of code I wrote during this period. I will keep in mind your suggestions, but I have one question. Can you elaborate on what you mean by "I really recommend writing your program in a more… mathematica-like way"? – brama Apr 17 '14 at 5:01
@brama Er……for example, you've used a lot of AppendTo, I'm not sure if it's because you're following the rule of thumb "the more my code inside Compile looks like C code I'd write otherwise, the more I benefit from it", but if it's what your original code looks like, you'd better set some effort on avoiding it. BTW, AppendTo is even disliked inside Compile`, see this post for details. A complete discussion for this issue might need a whole book 囧, have you already read Leonid's excellent book? – xzczd Apr 17 '14 at 5:17
Thank you so much!! :) – brama Apr 18 '14 at 13:22