3
$\begingroup$

My knowledge in parallel computing is very limited. I tried to figure out how to proceed using the mathematica help but did not really find what I needed. The programme looks as shown below. There are several loops, but what I would like is to parallelize the big loop (the one defined on r). I would be grateful if you had suggestions. Thanks a lot. Laurent

    ClearAll["Global`*"];
    kmax = 10000;
    m1 = 10.19357; s1 = 0.3022616; m2 = 11.795; s2 = s1;

    For[r = -0.99999999999, r <= 0.99999999999, r = r+0.09999999999,
     Clear[data];
     SeedRandom[1234];
     S = {{s1^2,r*s1*s2},{r*s1*s2,s2^2}};
     data = Sort[RandomVariate[LogMultinormalDistribution[{m1, m2}, S],kmax]];
     productivity = Table[data[[k]][[1]], {k, 1, kmax}];
     wagemin = Min[productivity];
     wagemax = Max[productivity];
     gamma = Table[data[[k]][[2]], {k, 1, kmax}];
     elasticity = Table[0.4/(wagemax - wagemin)*data[[k]][[1]] + 0.1 -
       0.4*wagemin/(wagemax - wagemin), {k, 1, kmax}];
     taxmax[r] = 0;

     For[tP = 0, tP < 1, tP = tP + 0.001,
     For[tC = 0, tC < tP, tC = tC + 0.001,
     Clear[wtilde];
     wtilde[g_, e_] := (((1 + e) g)/((1 - tC)^(1 + e) - (1 - tP)^(1+e)))^(1/(1 + e));
     For[k = 1, k <= kmax, k = k + 1,
       thresholds[k] =wtilde[gamma[[k]], elasticity[[k]]];];
     For[k = 1, k <= kmax, k = k + 1,
     If[productivity[[k]] <= thresholds[k], winf[k] = productivity[[k]], winf[k] = 0];
     If[productivity[[k]] <= thresholds[k], einf[k] = elasticity[[k]], einf[k] = 1];
     If[productivity[[k]] > thresholds[k], wmax[k] = productivity[[k]], wmax[k] = 0];
     If[productivity[[k]] > thresholds[k], emax[k] = elasticity[[k]], emax[k] = 1];];
     tax = tP*Sum[(1 - tP)^einf[k]*winf[k]^(1 + einf[k]), {k, 1, kmax}] +
       tC*Sum[(1 - tC)^emax[k]*wmax[k]^(1 + emax[k]), {k, 1, kmax}];
     gdp = Sum[(1 - tP)^einf[k]*winf[k]^(1 + einf[k]), {k, 1, kmax}] + 
       Sum[(1 - tC)^emax[k]*wmax[k]^(1 + emax[k]), {k, 1, kmax}];

      If[tax >= taxmax[r], taxmax[r] = tax; tPopt[r] = tP; tCopt[r] = tC;];
      ] (* End loop tC *);
      ] (* End loop tP*);

   Print["r =", r, " Optimum: tP =", tPopt[r], "; tC=", tCopt[r]];
$\endgroup$
  • 5
    $\begingroup$ Asking us to fix/rewrite your specific code for you won't probably be well received here, and may be closed as too localized. Perhaps you could reduce your example to a much simpler independent minimal working case, show us what you tried on that one, indicate why it didn't work, and then you will receive much better feedback. $\endgroup$ – MarcoB Feb 21 '16 at 18:04
  • $\begingroup$ One off-topic tip, For is not as optimized as Do and is (significantly) slower. I recommend using Do instead. (e.g. Do[ (*blah blah*), {r, -0.99999999999, 0.99999999999, 0.09999999999} ] and so on) $\endgroup$ – JungHwan Min Feb 21 '16 at 21:51
  • 1
    $\begingroup$ Try to rewrite your For loops into Map or Table statements. See the documentation on Map/Table for some examples. This is the hard part, but it is worth learning. Once you do this, you can replace them with ParallelMap or ParallelTable. $\endgroup$ – Searke Feb 22 '16 at 1:45
  • $\begingroup$ As a second note, your code is in one giant block. This is bad. You should try to define some functions so that your code is broken up into smaller parts. $\endgroup$ – Searke Feb 22 '16 at 1:46
  • $\begingroup$ You can gain at least a factor of two in speed, and perhaps much more, by writing the inner loop more efficiently. Concentrate of efficient code first, and worry about parallelization afterword, if you still need it. $\endgroup$ – bbgodfrey Feb 22 '16 at 5:32
2
$\begingroup$

Runtime for this code as written but with kmax reduced to 1000 and the tP and tC increments increased to 0.1, is 37 sec (AbsoluteTiming) on my four-processor PC after correcting a typo in the final line, a missing ]. (Without those parameter changes, runtime would be about 100000 times longer.)

Runtime can be reduced an order of magnitude by focusing on the two innermost For loops (i.e., For[k = 1, k <= kmax, k = k + 1, ...]) and Sums:

  • The Sums comprising tax and gdp are computed twice but need be computed only once.
  • When winf or wmax is zero, the corresponding term in the corresponding sum also is zero and need not be calculated.
  • MapThread is faster than For in this instance.

With these changes, runtime decreases to 12 sec. The outer loops also can be optimized,

  • Eliminate Sort, which is unnecessary.
  • Replace For-loops by Do-loops.
  • Replace Table by Transpose to compute productivity and gamma.
  • Replace Table by Map to compute elasticity

compressing the code to

ClearAll["Global`*"];
kmax = 1000; m1 = 10.19357; s1 = 0.3022616; m2 = 11.795; s2 = s1;

Do[
  taxmax = 0; SeedRandom[1234];
  S = {{s1^2, r*s1*s2}, {r*s1*s2, s2^2}};
  {productivity, gamma} = RandomVariate[LogMultinormalDistribution[{m1, m2}, S], kmax] 
    // Transpose;
  wagemin = Min[productivity]; wagemax = Max[productivity];
  elasticity = Map[(0.4/(wagemax - wagemin)*(# - wagemin) + 0.1) &, productivity];

  Do[
    wtilde[g_, e_] := (((1 + e) g)/((1 - tC)^(1 + e) - (1 - tP)^(1 + e)))^(1/(1 + e));
    thresholds = MapThread[wtilde, {gamma, elasticity}];
    sum1 = MapThread[If[#1 <= #2, (1 - tP)^#3*#1^(1 + #3), 0] &, 
      {productivity, thresholds, elasticity}] // Total;
    sum2 = MapThread[If[#1 > #2, (1 - tC)^#3*#1^(1 + #3), 0] &, 
      {productivity, thresholds, elasticity}] // Total;
    tax = tP*sum1 + tC*sum2; gdp = sum1 + sum2;
    If[tax >= taxmax, taxmax = tax; tPopt = tP; tCopt = tC;], 
  {tP, 0, 1 - 0.0001, 0.1}, {tC, 0, tP - 0.0001, 0.1}]

  Print["r =", r, " Optimum: tP =", tPopt, "; tC=", tCopt],
{r, -0.99999999999, 0.99999999999, 0.09999999999}]

although doing so saves negligible time. Finally, simply replacing the outermost Do by ParallelDo reduces runtime of this reduced problem to 4 sec, an overall reduction of a factor of ten.

Note that this particular problem does not require returning results from the subordinate kernels to the master kernel. To do so would require use of SetSharedVariable.

$\endgroup$
  • $\begingroup$ Thank you very much for your help. This was very useful and all the provided details helped me better understand how to write my code in a more efficient way. $\endgroup$ – Laurent Simula Feb 24 '16 at 9:14

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.