1
$\begingroup$

The following simulation takes over 20 seconds for evaluation. I have difficulty re-coding it with Nest and Map or Compile functions to make it more efficient. Can some one help with this?

The overview of the code is: It is a numerical solution for evolution of density of vehicles on a corridor that has L lanes. The corridor is split into n sections of dx length and each section has a sink and source coming from a ramp (which intern is divided into p segments each). the network is loaded from on-ramps with density A for m time periods only and simulated till the network is empty. The variable a1 is the control of inflow from the ramps on to the corridor. The outflow factor is B. The output of the simulation is the total density observed over the whole simulation.

Also, how can I use parallelization on this 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]];
dx = 1/6; n = 48; m = 300; p = 36; capacity = 2500; A = 18.; B = 0.1; L = 3.; 
RML = 30; Vf = 100; Kj = 150; w = 20; dt = 1/600;

f[a1_] := 
Module[{G,k0=ConstantArray[0, n],k, kr=Table[Table[0,{i1,1,p}],{i2,1,n}],P}, 
k = k0; RM = 100 a1; j = 0; RampSpill = 0;
For[i = 2, i < n, i++, kr[[i, 1]] = A]; 
NtwrkTT = TT = Plus @@ (Plus @@ kr);
While[TT > 0, TT = 0;
 Do[FQin = If[i == 2, Min[demand[L, k0[[i - 1]]], supply[L, k0[[i]]]], FQout];
 dem = demand[L, k0[[i]]]; dem = If[dem == 0, 0.001, dem];
 G = Min[1, supply[L, k0[[i + 1]]]/dem];
 P = G demand[1, kr[[i, p]]];
 Qr = (P - B FQin) dx;
 FQout = Min[demand[L, k0[[i]]], supply[L, k0[[i + 1]]]];
 k[[i]] = k0[[i]] + (FQin - FQout + Qr)/Vf;
 kr0 = kr[[i]];
 Do[MR = If[ir == RML + 1, RM, capacity];
  RQin = Min[MR, If[ir == 2, flo[1, kr0[[ir - 1]], kr0[[ir]]], RQout]];
  MR = If[ir == RML, RM, capacity];
  RQout =Min[MR, If[ir < p, flo[1, kr0[[ir]], kr0[[ir + 1]]], P]];
  kr[[i, ir]] = kr0[[ir]] + (RQin - RQout)/Vf, {ir, 2, p}];
 kr[[i, 1]] = If[j <= m, A, 0], {i, 2, n - 1}];
TT = Plus @@ (Plus @@ kr);
TT += Plus @@ k;
k0 = k; NtwrkTT += TT; j++];
NtwrkTT dt] // Timing
f[14]

Edit: Please find below the profile output:enter image description here

Improve this question
$\endgroup$
  • $\begingroup$ Perhaps you can break your problem down and identify the expensive parts. As it is, motivation to debug may be limited. $\endgroup$ – Yves Klett Mar 7 '14 at 19:47
  • $\begingroup$ @Yves Klett, Added the Profile output. $\endgroup$ – brama Mar 7 '14 at 21:50
  • $\begingroup$ You seem to be doing simple arithmetic operations inside your module. It is possible that this can be compiled efficiently... have you tried doing that? $\endgroup$ – rm -rf♦ Mar 7 '14 at 21:50
  • $\begingroup$ @rm -rf, I am a beginner in Mathematica. Can you describe more of what you mean by "compiled efficiently"? Do you mean, using` Compile[]` function? $\endgroup$ – brama Mar 7 '14 at 21:52
  • 1
    $\begingroup$ Is there some very specific function of Mathematica that you are using? It looks like a lot of loops with rather simple math in them. Here I would suggest not to use Mathematica if performance and scalability matters. In my experience you can speed up programs by roughly a factor of 1000 (non optimized Mathematica) if you switch to C/C++. Even when you optimize your code you wont reach that performance with Mathematica. $\endgroup$ – jens_bo Mar 8 '14 at 5:56
4
$\begingroup$

With

1. notes in the answer of your latter question,

2. initial values for FQout and RQout(the values are carelessly chosen real numbers because it's bound to be overlaid),

the modified code is about 40 times faster now:

f = With[{dx = 1/6, n = 48, m = 300, p = 36, capacity = 2500, A = 18.,
     B = 0.1, L = 3., RML = 30, Vf = 100, Kj = 150, w = 20, 
    dt = 1/600},
   Module[{demand, supply, flo},
    demand = {n, k} \[Function] Min[k Vf, n capacity];         
    supply = {n, k} \[Function] Min[(n Kj - k) w, n capacity];
    flo = {n, Ku, Kd} \[Function] Evaluate@Min[demand[n, Ku], supply[n, Kd]];

    Compile[{{a1, _Real}}, 
     Module[{G, k0 = Table[0., {n}], k, kr = Table[0., {n}, {p}], P, 
       j = 0, RM = 100 a1, RampSpill, NtwrkTT, TT, dem, FQin, Qr, 
       FQout = 0., kr0, MR, RQin, RQout = 0., i}, 
      k = k0; 
      RampSpill = 0;
      For[i = 2, i < n, i++, kr[[i, 1]] = A];
      NtwrkTT = TT = Plus @@ (Plus @@ kr);
      While[TT > 0, TT = 0;
       Do[
        FQin = If[i == 2, Min[demand[L, k0[[i - 1]]], supply[L, k0[[i]]]], FQout];
        dem = demand[L, k0[[i]]];
        dem = If[dem == 0, 0.001, dem];
        G = Min[1, supply[L, k0[[i + 1]]]/dem];
        P = G demand[1, kr[[i, p]]];
        Qr = (P - B FQin) dx;
        FQout = Min[demand[L, k0[[i]]], supply[L, k0[[i + 1]]]];
        k[[i]] = k0[[i]] + (FQin - FQout + Qr)/Vf;
        kr0 = kr[[i]];
        Do[
         MR = If[ir == RML + 1, RM, capacity];
         RQin = Min[MR, If[ir == 2, flo[1, kr0[[ir - 1]], kr0[[ir]]], RQout]];
         MR = If[ir == RML, RM, capacity];
         RQout = Min[MR, If[ir < p, flo[1, kr0[[ir]], kr0[[ir + 1]]], P]];
         kr[[i, ir]] = kr0[[ir]] + (RQin - RQout)/Vf, {ir, 2, p}];
        kr[[i, 1]] = If[j <= m, A, 0], {i, 2, n - 1}];
       TT = Plus @@ (Plus @@ kr);
       TT += Plus @@ k;
       k0 = k; NtwrkTT += TT; j++];
      NtwrkTT dt], 
     CompilationOptions -> {"InlineExternalDefinitions" -> True}]]];
Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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