3
$\begingroup$

I have developed (of course with a lot of help from experts in this forum) the following code to speed up the calculations for large digraphs. However, I ended up with much slower calculations compared to Do loop formulation.

ClearAll[n, d, G, aG, aGsym, Xt, Xt1, alfa, eta, gama, tao, teta, one,
   onePlusEta, Aterm, Bterm, Cterm, Dterm, eqs, modelFull];

(* generating a digraph *)
SeedRandom[14];
n = 10; 
d = 0.5; 
G = RandomGraph[{Round[n], Round[n*(n - 1)*d]}, DirectedEdges -> True]; 
aG = AdjacencyMatrix[G]; 
aGsym = BitAnd[aG, Transpose[aG]]; 

(* preliminaries of matrix inputs for the model *)
Xt = DiagonalMatrix[Array[Subscript[x, #][t] &, n]]; 
Xt1 = Array[Subscript[x, #]'[t] &, n];
alfa = Transpose@Array[Subscript[\[Alpha], #] &, {n, n}];
eta = Array[Subscript[\[Eta], #] &, n];
gama = Transpose[Array[Subscript[\[Gamma], #] &, {n, n}]]*aG;
tao = Array[Subscript[\[Tau], #] &, {n, n}];
teta = Transpose[Array[Subscript[\[Theta], #] &, {n, n}]*aG];
one = Array[1 &, n];

(* building the model *)
onePlusEta = DiagonalMatrix[one + eta]*Xt; Aterm = (gama + teta).Xt;  
Bterm = Xt.((alfa*tao)*aG); 
Cterm = ((alfa*alfa)*aGsym)*(Xt.Array[Diagonal[Xt] &, {n}])*aGsym;
Dterm = ((gama*gama)*aGsym)*
  Transpose[Array[(Subscript[x, #][t])^2 &, {n, n}]*aGsym];
eqs = Total[onePlusEta - Aterm + Bterm + Cterm - Dterm]//FullSimplify;   
modelFull = Thread[Xt1 == eqs - Array[Subscript[x, #][t] &, n]];     

(* parameter values for numeric solutions *)
ClearAll[params, parameters, initialVals, sol];
SeedRandom[14];
params = {abs, random, absCost, trsCost, trs} = {
    RandomReal[{0.01, 1}, n], RandomReal[{-0.5, 0.5}, n],
    RandomReal[{0.01, 0.5}, n], RandomReal[{0.01, 1}, n],  
    RandomReal[{0.01, 0.5}, n] 
    };
parameters = Flatten[Thread[#[[1]] -> #[[2]]] & /@Thread[Table[Subscript[#, n], {n, n}] & /@ {\[Alpha], \[Eta], \[Gamma], \[Tau], \[Theta]} -> params]];

initialVals = Array[Subscript[x, #][0] == RandomReal[{0.01, 1}] &, n];
sol = NDSolve[{modelFull /. parameters, initialVals} // Flatten, 
    Array[Subscript[x, #][t] &, n], {t, 50},
    Method-> "StiffnessSwitching"] // Flatten;

Plot[Evaluate[Array[Subscript[x, #][t] &, n] /. sol], {t, 0, 50}, 
    PlotRange -> All]

Is there any specific reason for the slow process? How can I really improve the performance of the code with larger digraphs?

OR Am I on the WRONG track for speeding up the calculations?

Please give me some ideas about improving the calculations.

Note:

I can also send the other code using Do and If if somebody is interested in.

$\endgroup$
  • 2
    $\begingroup$ Code works fast. And what is calculated? $\endgroup$ – Alex Trounev Oct 30 '18 at 21:19
  • $\begingroup$ @AlexTrounev: Thank you very much for checking the speed of the code. The code uses a digraph G to create a system of ODEs and then solve the system sol using the created system. I would appreciate if you run the code with n=50 and n=100 and give me the timing for each one. MyMathematica version is 10.0.2.0 and use Windows 7 (with RAM 4). Also tell me on which machine you run the model. Thank you again. $\endgroup$ – Tugrul Temel Oct 30 '18 at 21:58
  • 2
    $\begingroup$ NDSolve spits out ndsz warning and fails at about t==5.57, tested in v9.0.1 and v11.2. $\endgroup$ – xzczd Oct 31 '18 at 6:47
3
$\begingroup$

I found that the main loss of time and resources occurs while simplifying the equations in line eqs = Total[onePlusEta - Aterm + Bterm + Cterm - Dterm]//FullSimplify;. Moreover, this simplification with a large number of equations leads to errors. I replaced FullSimplify with Simplify[eqs = Total[onePlusEta - Aterm + Bterm + Cterm - Dterm], TimeConstraint -> ts]; with the option TimeConstraint and added to the code a convenient analysis of the influence of the number of equations, the time of symbolic and numerical calculations.

AbsoluteTiming[
 ClearAll[n, d, G, aG, aGsym, Xt, Xt1, alfa, eta, gama, tao, teta, 
  one, onePlusEta, Aterm, Bterm, Cterm, Dterm, eqs, modelFull];
 SeedRandom[14];

 (*generating a digraph*)
 n = 100; ts = 1; tm = 2; d = 0.5; 
 G = RandomGraph[{Round[n], Round[n*(n - 1)*d]}, 
   DirectedEdges -> True]; aG = AdjacencyMatrix[G]; 
 aGsym = BitAnd[aG, 
   Transpose[
    aG]];(*picks symmetric cells from "aG"*)(*preliminaries of matrix \
inputs for the model*)
 Xt = DiagonalMatrix[
   Array[Subscript[x, #][t] &, 
    n]];(*Subscript[x,i][t] is vertex i at time t*)
 Xt1 = Array[Subscript[x, #]'[t] &, 
   n];(*time derivative of Subscript[x,i][t]*)
 alfa = Transpose@Array[Subscript[\[Alpha], #] &, {n, n}];
 eta = Array[Subscript[\[Eta], #] &, n];
 gama = Transpose[Array[Subscript[\[Gamma], #] &, {n, n}]]*aG;
 tao = Array[Subscript[\[Tau], #] &, {n, n}];
 teta = Transpose[Array[Subscript[\[Theta], #] &, {n, n}]*aG];
 one = Array[1 &, n];
 (*building the model*)onePlusEta = DiagonalMatrix[one + eta]*Xt; 
 Aterm = (gama + teta).Xt;
 Bterm = Xt.((alfa*tao)*aG);
 Cterm = ((alfa*alfa)*aGsym)*(Xt.Array[Diagonal[Xt] &, {n}])*aGsym;
 Dterm = ((gama*gama)*aGsym)*
   Transpose[Array[(Subscript[x, #][t])^2 &, {n, n}]*aGsym];
 Simplify[eqs = Total[onePlusEta - Aterm + Bterm + Cterm - Dterm], 
  TimeConstraint -> ts];
 modelFull = Thread[Xt1 == eqs - Array[Subscript[x, #][t] &, n]];

 (*parameter values for numeric solutions*)
 ClearAll[params, parameters, initialVals, sol];
 SeedRandom[14];
 params = {abs, random, absCost, trsCost, 
    trs} = {RandomReal[{0.01, 1}, n], RandomReal[{-0.5, 0.5}, n], 
    RandomReal[{0.01, 0.5}, n], RandomReal[{0.01, 1}, n], 
    RandomReal[{0.01, 0.5}, n]};
 parameters = 
  Flatten[Thread[#[[1]] -> #[[2]]] & /@ 
    Thread[Table[
         Subscript[#, n], {n, 
          n}] & /@ {\[Alpha], \[Eta], \[Gamma], \[Tau], \[Theta]} -> 
      params]];

 initialVals = Array[Subscript[x, #][0] == RandomReal[{0.01, 1}] &, n];
 sol = NDSolve[{modelFull /. parameters, initialVals} // Flatten, 
    Array[Subscript[x, #][t] &, n], {t, tm}, 
    Method -> "StiffnessSwitching"] // Flatten;


 Plot[Evaluate[Array[Subscript[x, #][t] &, n] /. sol], {t, 0, tm}, 
  PlotRange -> All, 
  PlotLabel -> Grid[{{"n", "ts", "tm"}, {n, ts, tm}}, Frame -> All]]]

fig1 In my version, simplification does not affect the result, but increases the total time by 7-8 times.

$Version

Out[]= "11.3.0 for Microsoft Windows (64-bit) (March 7, 2018)"
$\endgroup$
  • $\begingroup$ Thank you so much for the neat analysis and result. With your fix, now the model behaves itself and works smoothly, and takes only 1 minute for n=100. I can live with this speed. $\endgroup$ – Tugrul Temel Oct 31 '18 at 9:40
  • 2
    $\begingroup$ You're welcome! You can do this analysis by simply putting on each suspicious line //AbsoluteTiming $\endgroup$ – Alex Trounev Oct 31 '18 at 9:55

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.