0
$\begingroup$

I'm trying to construct a large table where new entries depend on the most recent entries. It's quite slow: takes 2 mins and I'm embedding this into a function that does it repeatedly so any reduction in time is a big savings. I think this issue is with "AppendTo" and I've seen repeated references on here to use Reap/Sow, but I'm having trouble implementing that because the next entries in the table are conditional on the previous. Here is an example of the code I'm working with right now.

SeedRandom[1234];

TransMark = Table[RandomReal[], {i, 7}, {j, 7}]; (*Just an example of a similar matrix to what I have*)

RandDraw = Table[RandomReal[], {i, 10000}, {j, 576}];

WorkerMat = Table[{{2, 1}}, {10000}]; 

Timing[Do[
   If[RandDraw[[i, j]] <= 
     TransMark[[All, 1]][[Last[WorkerMat[[i]]][[1]]]],
    If[Last[WorkerMat[[i]]][[1]] == 1, 
     Tenure = Last[WorkerMat[[i]]][[2]] + 1, Tenure = 1];
    AppendTo[WorkerMat[[i]], {1, Tenure}],
    If[RandDraw[[i, j]] <= 
      TransMark[[All, 2]][[Last[WorkerMat[[i]]][[1]]]],
     If[Last[WorkerMat[[i]]][[1]] == 2, 
      Tenure = Last[WorkerMat[[i]]][[2]] + 1, Tenure = 1];
     AppendTo[WorkerMat[[i]], {2, Tenure}],
     If[RandDraw[[i, j]] <= 
       TransMark[[All, 3]][[Last[WorkerMat[[i]]][[1]]]],
      If[Last[WorkerMat[[i]]][[1]] == 3, 
       Tenure = Last[WorkerMat[[i]]][[2]] + 1, Tenure = 1];
      AppendTo[WorkerMat[[i]], {3, Tenure}],
      If[RandDraw[[i, j]] <= 
        TransMark[[All, 4]][[Last[WorkerMat[[i]]][[1]]]],
       If[Last[WorkerMat[[i]]][[1]] == 4, 
        Tenure = Last[WorkerMat[[i]]][[2]] + 1, Tenure = 1];
       AppendTo[WorkerMat[[i]], {4, Tenure}],
       If[
        RandDraw[[i, j]] <= 
         TransMark[[All, 5]][[Last[WorkerMat[[i]]][[1]]]],
        If[Last[WorkerMat[[i]]][[1]] == 5, 
         Tenure = Last[WorkerMat[[i]]][[2]] + 1, Tenure = 1];
        AppendTo[WorkerMat[[i]], {5, Tenure}],
        If[
         RandDraw[[i, j]] <= 
          TransMark[[All, 6]][[Last[WorkerMat[[i]]][[1]]]],
         If[Last[WorkerMat[[i]]][[1]] == 6, 
          Tenure = Last[WorkerMat[[i]]][[2]] + 1, Tenure = 1];
         AppendTo[WorkerMat[[i]], {6, Tenure}],
         If[Last[WorkerMat[[i]]][[1]] == 7, 
          Tenure = Last[WorkerMat[[i]]][[2]] + 1, Tenure = 1];
         AppendTo[WorkerMat[[i]], {7, Tenure}];]]]]]] , {j, 1, 
    576}, {i, 1, 10000}];
 ];

Which takes over 2 mins on my machine. My apologies if this is a repeat question, but I've been looking for a similar post for many hours and have not discovered one. Thanks in advance!

$\endgroup$
  • 2
    $\begingroup$ I'm voting to close this question as off-topic because it is too localized; i.e, it applies only to the local situation and needs of its poster and answers will not benefit others. $\endgroup$ – m_goldberg Dec 14 '17 at 23:34
  • 4
    $\begingroup$ Can you construct a minimal working example that illustrates your issue? I suspect the essence of the problem can be described with much less code. $\endgroup$ – jjc385 Dec 15 '17 at 1:30
  • $\begingroup$ NestList, FoldList and their "See Also"s are probably good places to start looking. $\endgroup$ – aardvark2012 Dec 15 '17 at 3:00
9
$\begingroup$

There are two issues with your code: First, it is very procedural with lots of If statements. Mathematica is slow at processing procedural code. Good news: Procedural code in Mahematica can usually be compiled with Compile in order to speed up the execution.

The second and major flaw of your code however is that you use AppendTo quite often. Every time you use it, a list has to be copied which makes the code very memory bound (and thus slow).

Fortunately, the size of WorkingMat is known right from the start and Last[WorkerMat[[i]]] refers to WorkerMat[[i,j]] and Append[Last[WorkerMat[[i]]],value] can be replaced by WorkerMat[[i,j+1]] = value. Moreover, one can make the j-loop the outer loop and parallelize it. Afterwards, one obtains a code which is nicely compilable. Replacing the nested If construction by a lookup for the first entry in a certain list which is greater or equal a certain value, the result is the CompiledFunction cf below (I am skipping a significant amount of refactorization steps here.)

cf = Compile[{{rand, _Real, 1}, {TransMark, _Real, 2}},
  Block[{r, a, b, k},
   a = Table[0, {i, 1, Length[rand] + 1}, {j, 1, 2}];
   a[[1, 1]] = 2;
   a[[1, 2]] = 1;
   Do[
    r = Compile`GetElement[rand, j];
    b = Compile`GetElement[a, j, 1];
    k = 1;
    While[r > Compile`GetElement[TransMark, b, k] && k < 7, k++];
    a[[j + 1, 1]] = k;
    a[[j + 1, 2]] = If[b == k, Compile`GetElement[a, j, 2] + 1, 1];
    , {j, 1, Length[rand]}];
   a
   ],
  RuntimeAttributes -> {Listable},
  Parallelization -> True,
  CompilationTarget -> "C",
  RuntimeOptions -> "Speed"
  ]

Btw, the Compile`GetElement statements in the code above are read-only instructions that replace Part. I think it instructs the C code generator to skip certain bound checks. This can improve the performance considerably but makes it less "safe".

Generating some data...

SeedRandom[1234];
TransMark = RandomReal[{0, 1}, {7, 7}];
n = 10000;
m = 576;
RandDraw = RandomReal[{0, 1}, {n, m}];

... and running the code on the data

WorkerMat2 = cf[RandDraw, TransMark]; // AbsoluteTiming//First

0.05216

Your original code took about 86 seconds on my machine, so that's a speed-up of factor 1600 or so.

$\endgroup$
0
$\begingroup$

This is substantially faster and might be equivalent. If it is not, I may well have made mistakes in translating and, as was noted, a minimal working example is a useful commodity.

SeedRandom[1234];
TransMark = RandomReal[{0, 1}, {7, 7}];
maxi = 10000;
maxj = 576;
RandDraw = RandomReal[{0, 1}, {maxi, maxj}];
WorkerMat = Table[{{2, 1}}, {maxi}];

AbsoluteTiming[
 wm2 = Table[
   last = Last[WorkerMat[[i]]];
   rndi = RandDraw[[i]];
   ltable = Table[
     rndij = rndi[[j]];
     {l1, l2} = last;
     tm1 = TransMark[[l1]];
     tenure = 1;
     Which[rndij <= tm1[[1]], k = 1,
      rndij <= tm1[[2]], k = 2,
      rndij <= tm1[[3]], k = 3,
      rndij <= tm1[[4]], k = 4,
      rndij <= tm1[[5]], k = 5,
      rndij <= tm1[[6]], k = 6,
      True, k = 7
      ];
     If[l1 == k, tenure = l2 + 1];
     last = {k, tenure}
     , {j, maxj}];
   Developer`ToPackedArray[ltable]
   , {i, maxi}];]

(* Out[374]= {27.892029, Null} *)

wm3 = Join[WorkerMat, wm2, 2];

That last should give the original result.

I'm sure a further speed gain could be had using Compile.

$\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.