Compile with array operations

Question

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,Compile`FunctionVariable$266]. >>

Compile::cprank: Compile cannot determine the rank of the result tensor. The length of CompileFunctionVariable$266 in MapThread[Plus,Compile`FunctionVariable$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 Compile`GetElement[qsum,Compile`Variable$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: Compile`GetElement[qsum,System`Private`CompileSymbol[0]] should be a tensor of type Integer, Real, or Complex; evaluation will use the uncompiled function. >> Compile::cplist: Compile`GetElement[qsum,System`Private`CompileSymbol[0]] should be a tensor of type Integer, Real, or Complex; evaluation will use the uncompiled function. >> Compile::cplist: Compile`GetElement[qsum,System`Private`CompileSymbol[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,Compile`FunctionVariable$3193]. >> Compile::cprank: Compile cannot determine the rank of the result tensor. The length of CompileFunctionVariable$3193 in MapThread[Plus,Compile`FunctionVariable$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.

Answer 1

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 = 
       MapThread[
        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 = 
       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], ϕ}];
      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 = 
         MapThread[
          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 = 
         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], ϕ}];
        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.