Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

The following code works fine with the the function f defined as a module, however, when I convert it to be a compiled function, it gives me errors. Can some one help me identify the problem is with the following code:

demand[n_, k_] := Min[k Vf, n capacity]; 
supply[n_, k_] := Min[(n Kj - k) w, n capacity]; 
flo[n_, ku_, kd_] := Min[demand[n, ku], supply[n, kd]];
gamma[ku_, kd_] := Min[1, supply[L, kd]/(demand[L, ku] + 0.001)]; 
inflow [phi_, FQin_] := (phi - β FQin) dx; density[qin_, qout_, qr_] := 
  (qin - qout + qr)/Vf;

dx = 1/6; capacity = 7500.; Kj = 150.; w = 100.; Vf = 100.; n = 60; m = 300; p = 24; 
RML = 18;α = 4300.; β = 0.1; L = 3.; delta = 1.; dt = 4./3600.; 

f = Compile[{{a1, _Real}, {a2, _Real} {a3, _Real}, {a4, _Real}, {a5,_Real}}, 
  Module[{k0 = ConstantArray[0, n], 
          kr = Table[Table[0, {i1, 1, p}], {i2, 1, n - 2}], 
          ϕ = 0,  
          Shutoff = False},
    j = 0;
    RMori = 
      Table[
        If[i1 dx < 2, 100 a1, 
          If[i1 dx < 4, 100 a2, 
            If[i1 dx < 6, 100 a3, 
              If[i1 dx < 8, 100 a4, 100 a5]]]], 
        {i1, 1,n - 2}];
    kr = ReplacePart[#, 1 -> (α delta/Vf)] & /@ kr; 
    NtwrkTT = TT = Plus @@ (Plus @@ kr);
    While[TT > 0, 
      TT = 0;
      qf = MapThread[flo[L, ##] &, {Most@k0, Rest@k0}];
      ϕ = demand[1, #[[-1]] ] & /@ kr MapThread[gamma, {k0[[2 ;; -2]], k0[[3 ;;]]}];
      qin = inflow[ϕ, Most@qf];
      k0 += Join[{0}, density[Most@qf, Rest@qf, qin], {0}];
      RM = MapThread[Min[#1, #2] &, {RMori, flo[1, #[[RML]], #[[RML + 1]]] & /@ kr}];
      qr = 
        MapThread[Join[#1, {#2}, #3, {#4}] &, 
          {flo[1, ##] & @@@ Partition[#[[;; RML]], 2, 1] & /@ kr, 
           RM, 
           flo[1, ##] & @@@ Partition[#[[RML + 1 ;;]], 2, 1] & /@ kr, ϕ}];
      kr = 
        MapThread[(#1 + #2) &, 
          {Join[{0}, #] & /@ (density[Most@#,Rest@#, 0] & /@ qr), kr}];
      If[j > m && Shutoff == False, Shutoff = True; 
      kr = ReplacePart[#, 1 -> 0] & /@ kr];
      TT += Plus @@ k0; 
      TT += Plus @@ (Plus @@ kr); 
      NtwrkTT += TT; 
      j++;];
    NtwrkTT dt]]

When the function is evaluated at f[39, 39, 39, 39, 39], it gives the following errors

ReplacePart::argrx: ReplacePart called with 2 arguments; 3 arguments are expected. >>
ReplacePart::argrx: ReplacePart called with 2 arguments; 3 arguments are expected. >>
ReplacePart::argrx: ReplacePart called with 2 arguments; 3 arguments are expected. >>
General::stop: Further output of ReplacePart::argrx will be suppressed during this calculation. >>
Compile::cpintlt: 3;;All at position 2 of k0[[3;;All]] should be either a nonzero integer or a vector of nonzero integers; evaluation will use the uncompiled function. >>
Compile::cpintlt: 3;;All at position 2 of k0[[3;;All]] should be either a nonzero integer or a vector of nonzero integers; evaluation will use the uncompiled function. >>
Compile::cset: Variable \[Phi] of type {_Integer,0} encountered in assignment of type {_Integer,1}. >>
Compile::cplist: density[Most[qf],Rest[qf],qin] should be a tensor of type Integer, Real, or Complex; evaluation will use the uncompiled function. >>
Compile::cpapot: Compilation of (flo[1,##1]&)@@Compile`GetElement[System`Private`CompileSymbol[4],System`Private`CompileSymbol[5]] is not supported for the function argument flo[1,##1]&. The only function arguments supported are Times, Plus, or List. Evaluation will use the uncompiled function. >>
Compile::cpapot: Compilation of (flo[1,##1]&)@@Compile`GetElement[System`Private`CompileSymbol[4],System`Private`CompileSymbol[5]] is not supported for the function argument flo[1,##1]&. The only function arguments supported are Times, Plus, or List. Evaluation will use the uncompiled function. >>
Compile::cpapot: Compilation of (flo[1,##1]&)@@Compile`GetElement[System`Private`CompileSymbol[4],System`Private`CompileSymbol[5]] is not supported for the function argument flo[1,##1]&. The only function arguments supported are Times, Plus, or List. Evaluation will use the uncompiled function. >>
General::stop: Further output of Compile::cpapot will be suppressed during this calculation. >>
Compile::cpintlt: RML+1;;All at position 2 of Compile`GetElement[kr,System`Private`CompileSymbol[0]][[RML+1;;All]] should be either a nonzero integer or a vector of nonzero integers; evaluation will use the uncompiled function. >>
General::stop: Further output of Compile::cpintlt will be suppressed during this calculation. >>
Compile::cset: Variable kr of type {_Integer,2} encountered in assignment of type {_Integer,0}. >>
Compile::cset: Variable \[Phi] of type {_Integer,0} encountered in assignment of type {_Integer,1}. >>
General::stop: Further output of Compile::cset will be suppressed during this calculation. >>
Compile::cplist: density[Most[qf],Rest[qf],qin] should be a tensor of type Integer, Real, or Complex; evaluation will use the uncompiled function. >>
share|improve this question
2  
What error does it give? –  RunnyKine Mar 20 at 19:29
    
Sorry for the confusion, but please see the edited question. –  brama Mar 20 at 20:07
    
It contains kr = ReplacePart[#, 1 -> (α delta/Vf)] & /@ kr; So you quite clearly invoke ReplacePart with 2 args, and if it wants 3 then it's no surprise when it complains. –  Daniel Lichtblau Mar 21 at 0:21
    
ReplacePart only needs 2 arguments and works fine when used in the Module, please see reference.wolfram.com/mathematica/ref/ReplacePart.html, but MMA wants 3 args when I recode it in a Compile function. This is confusing!! –  brama Mar 21 at 9:36
    
The syntax of ReplacePart (documentation for V5.2) has changed. It might be that Compile hasn't learned that change (a probable reason is that newer versions with the new syntax can use patterns which Compile couldn't support anyway). In an uncompiled evalation the old syntax still seems to work... –  Albert Retey Mar 21 at 10:57

1 Answer 1

up vote 1 down vote accepted

With

1. the note about the behavior of ReplacePart inside Compile made by Daniel Lichtblau and Albert Retey in the comments above,

2. the correction for the simple mistake {a2, _Real} {a3, _Real},

3. the note about All inside Compile here,

4. a modification for the uncompiled definition of qr i.e. changing

qr =  MapThread[Join[#1, {#2}, #3, {#4}] &, 
   {flo[1, ##] & @@@ Partition[#[[;; RML]], 2, 1] & /@ kr, 
    RM, 
    flo[1, ##] & @@@ Partition[#[[RML + 1 ;;]], 2, 1] & /@ kr, ϕ}];

into

qr = MapThread[Join[#1, {#2}, #3, {#4}] &, 
   {MapThread[flo[1, #, #2] &, {#[[1 ;; RML - 1]], #[[2 ;; RML]]}] & /@ kr, 
    RM, 
    MapThread[flo[1, #, #2] &, {#[[RML + 1 ;; -2]], #[[RML + 2 ;; -1]]}] & /@ kr, ϕ}];

5. notes in my answer for your latter question,

I fixed your code, now for f[39, 39, 39, 39, 39] it's about 40 times faster than the uncompiled one:

f = With[{dx = 1/6, capacity = 7500., Kj = 150., w = 100., Vf = 100., n = 60, m = 300, 
          p = 24, RML = 18, α = 4300., β = 0.1, L = 3., delta = 1., dt = 4./3600.}, 
   Module[{demand, supply, flo, gam, inflow, density}, 
    demand = Function[{n, k}, Min[k Vf, n capacity]]; 
    supply = Function[{n, k}, Min[(n Kj - k) w, n capacity]]; 
    flo = Function[{n, ku, kd}, Evaluate[Min[demand[n, ku], supply[n, kd]]]]; 
    gamma = Function[{ku, kd}, Evaluate[Min[1, supply[L, kd]/(demand[L, ku] + 0.001)]]]; 
    inflow = Function[{phi, FQin}, (phi - β FQin) dx]; 
    density = Function[{qin, qout, qr}, (qin - qout + qr)/Vf]; 
    Compile[{{a1, _Real}, {a2, _Real} , {a3, _Real}, {a4, _Real}, {a5, _Real}}, 
     Module[{k0 = Table[0., {n}], 
       kr = Table[0., {n - 2}, {p}], 
       ϕ = {0.}, Shutoff = False, 
       j = 0, RMori, NtwrkTT, TT, qf, qin, RM, qr},
      RMori = Table[If[i1 dx < 2, 100 a1, 
         If[i1 dx < 4, 100 a2, 
          If[i1 dx < 6, 100 a3, 
           If[i1 dx < 8, 100 a4, 100 a5]]]], {i1, 1, n - 2}];
      kr = (ReplacePart[#1, (α delta)/Vf, 1] &) /@ kr;
      NtwrkTT = TT = Plus @@ Plus @@ kr;
      While[TT > 0,
       TT = 0;
       qf = MapThread[flo[L, ##] &, {Most[k0], Rest[k0]}]; 
       ϕ = (demand[1, #1[[-1]]] &) /@ kr MapThread[gamma, {k0[[2 ;; -2]], k0[[3 ;; -1]]}]; 
       qin = inflow[ϕ, Most[qf]]; 
       k0 += Join[{0}, density[Most[qf], Rest[qf], qin], {0}];        
       RM = MapThread[Min[#1, #2]&, {RMori, (flo[1, #1[[RML]], #1[[RML + 1]]] &)/@kr}];
       qr = MapThread[Join[#1, {#2}, #3, {#4}] &, 
             {MapThread[flo[1, #, #2] &, {#[[1 ;; RML - 1]], #[[2 ;; RML]]}] & /@ kr, 
              RM, 
              MapThread[flo[1, #, #2] &, {#[[RML + 1 ;; -2]], #[[RML + 2 ;; -1]]}] &/@kr, 
              ϕ}];        
       kr = MapThread[(#1 + #2) &, 
         {Join[{0}, #] & /@ (density[Most@#, Rest@#, 0] & /@ qr), kr}];
       If[j > m && Shutoff == False, 
          Shutoff = True; kr = (ReplacePart[#1, 0, 1] &) /@ kr]; TT += Plus @@ k0; 
          TT += Plus @@ Plus @@ kr; NtwrkTT += TT; j++;]; NtwrkTT dt], 
     CompilationOptions -> {"InlineExternalDefinitions" -> True}]]];
share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.