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.

I would like to learn writting more efficiently. I am working on problem of simulating random walks (abnormal diffusion). I believe this is good example of inefficient code. How can I make it faster?

xvv[n_, a_] := (
  xx = #1 Cos[#2] &;
  yy = #1 Sin[#2] &;
  lval = RandomVariate[ParetoDistribution[0.01, a], n]; 
  phival = RandomReal[{0, 2 Pi}, n];
  xv = Thread[xx[lval, phival]];
  yv = Thread[yy[lval, phival]];
  {xv, yv}

mm[n_, a_] := (
  laz = xvv[n, a];
  {Sqrt[laz[[1, All]]^2 + laz[[2, All]]^2], laz[[1, All]]}

xm2[tmax_, ni_, a_] := (
  t = 0; x = 0;j=1; values = Transpose[mm[10000, a]];
  While[t < tmax,
   {t, x} = {t + 0.1 RandomVariate[ParetoDistribution[0.1, ni]], x} + values[[j]];j++;]; x

tab2[tmax_, ni_, a_] := Table[xm2[tmax, ni, a], {i, 1, 100}];

valt = Table[10^(j), {j, 1, 1.7, .2}]
Ngamm[mi_, ni_] := (
  a = mi - 1;
  a3 = Table[{i, MedianDeviation[tab2[i, ni - 1, a]]^2}, {i, valt}]; 
  model = FindFit[Log[a3], aa xxx + b, {aa, b}, xxx];
  aa /. model

Edited (I accidentaly deleted part of xm2 (j=1; and j++;)) (now it should work). mi and ni in Ngamm should be larger than 1. (range that is interesting is 1 to 5).

Thank you for responses.

**Edit: Now codel looks a bit better, but it is not much faster.

xm3[tmax_, ni_, a_] := Module[{augvalues =
      0.1 RandomVariate[ParetoDistribution[0.1, ni], 10000]], 
      0.1 RandomVariate[ParetoDistribution[0.1, a], 10000] Cos[
        RandomReal[{0, 3.1415}, 10000]]]}, j = 1},
  While[augvalues[[1, j]] < tmax, j++]; augvalues[[2, j]]]

tab2[tmax_, ni_, a_] := Table[xm3[tmax, ni, a], {i, 1, 50}];

valt = Table[10^(j), {j, 1.0, 2.8, .1}]
  ni_] := (a3 = 
   Table[{i, MedianDeviation[tab2[i, ni - 1, mi - 1]]^2}, {i, valt}];
  model = FindFit[Log[a3], aa xxx + b, {aa, b}, xxx];
  aa /. model)
share|improve this question

closed as off-topic by Yves Klett, ciao, bobthechemist, m_goldberg, Michael E2 Mar 16 '14 at 16:10

This question appears to be off-topic. The users who voted to close gave this specific reason:

  • "The question is out of scope for this site. The answer to this question requires either advice from Wolfram support or the services of a professional consultant." – Yves Klett, ciao, bobthechemist, m_goldberg, Michael E2
If this question can be reworded to fit the rules in the help center, please edit the question.

This will be useful to read through: mathematica.stackexchange.com/questions/29349/… –  Szabolcs Mar 16 '14 at 2:18
Your code is giving me errors. What sort of input are you expecting for Ngamm? I'm guessing positive integers, but it's not clear. –  Verbeia Mar 16 '14 at 3:27
"values[[j]]".what is j? –  Apple Mar 16 '14 at 3:56
j is equal to the i in tab2? –  Apple Mar 16 '14 at 3:59

2 Answers 2

up vote 3 down vote accepted

Regarding your edit:

There are two thing you can do to speed up xm3. Firstly, the While loop can be replaced with faster code. (By the way, you should be aware of LengthWhile, which does what you want without having to explicity increment a counter.) However, for this case, we can use fast numerical functions to accomplish the task.

j = Total[UnitStep[tmax - augvalues[[1]]]]

Secondly, xm3 returns the value of augvalues[[2, j]]], with a single value of j. And yet you have computed it for every value of j up to 10000. It is far more efficient to just create a list of j values of the Cos term and total them.

The re-written xm3 looks like:

xm3[tmax_, ni_, a_] := Module[{augvalues1, j},
  augvalues1 = Accumulate[0.1 RandomVariate[ParetoDistribution[0.1, ni], 10000]];
  j = Total[UnitStep[tmax - augvalues1]];
  0.1 Total[RandomVariate[ParetoDistribution[0.1, a], j] Cos[RandomReal[{0, 3.1415}, j]]]]

This gives about a factor of 5 increase in speed.

share|improve this answer

Your code is throwing errors for me, but that is probably because it isn't clear what input you require to Ngamm, and because there is a values[[j]] element in the definition of xm2 but the j doesn't seem to come from anywhere.

Some general comments about your code follow.

There are some elements that you can pull out into separate definitions instead of defining them in a function that is itself defined using SetDelayed (:=). Doing it the way you are means that you are defining these things every time you call the function. This is not efficient, but it is probably not the cause of the slowdown.

You can use MapThread instead of Thread in the definition of xvv. I suspect that will be a bit more efficient.

Also in xvv and some of your other functions, you are using local variables or constants but they are not localised to these functions. Have a look at the documentation for Module, With and Block, and at this question.

So so far, this results in

xx = #1 Cos[#2] &;
yy = #1 Sin[#2] &;

xvv[n_, a_] := With[{
   lval = RandomVariate[ParetoDistribution[0.01, a], n],
   phival = RandomReal[{0, 2 Pi}, n]},
  {MapThread[xx, {lval, phival}], MapThread[yy, {lval, phival}]}]

You can simplify mm as well, and modularise, though on my testing I don't think it really speeds things up that much. It might scale a bit better. In particular, you don't need to specify All as the second dimension Part specification - it's assumed. Here's my version. Notice how I only call the squaring once, by adding up the two rows of laz^2 directly, by applying plus.

mm2[n_, a_] := With[{laz = xvv[n, a]},
  {Sqrt[Plus @@ laz^2], laz[[1]]}]

Anyway on my machine, mm or mm2 only take about 0.03 seconds for n->10000, so this isn't where the speedup in evaluation can be found. Instead, it's in the next sub function that your troubles arise. In the version above, j isn't defined anywhere. You can fix this up like so:

xm2[tmax_, ni_, a_] := (t = 0; x = 0; j = 1; 
  values = Transpose[mm[10000, a]];
  While[t < tmax, {t, x} = {t + 0.1 RandomVariate[ParetoDistribution[0.1, ni]], x} + 
     values[[j]]; j++]; x)

Now, what I don't understand about this function is that in values you already have a random ParetoDistribution vector being cumulated into itself, but you're also adding another one. Are you sure that's what you want? Assuming it really is, it occurs to me that you could use one of the definitions of FoldWhileList provided in the answers to this question. I am going to use Mr.Wizard's version, since it seems to scale well.

 FoldWhileList[f_, start_, rest_, test_] :=
 Module[{bag = Internal`Bag[start], g},
  g[x_?test] := (Internal`StuffBag[bag, x]; x);
  g[else_] := Return[Null, Fold];
  Fold[g @ f @ ## &, start, rest];
  Internal`BagPart[bag, All]

xm4[tmax_, ni_, a_] := 
 With[{addvalues = Transpose[mm[10000, a]] + 
   ({#, 0} & /@ (0.1 RandomVariate[ParetoDistribution[0.1, ni], 10000]))}, 
  FoldWhileList[Plus, {0, 0}, addvalues, First[#] < tmax &][[-1, 2]]]

The advantage with this approach is that you don't need to track the length of the resulting accumulated vector. The only issue here is that this will give you the last value of values before the test is failed, rather than the first position where it is.

An alternative approach would be to accumulate the random numbers first, which is fast, and then find the first position where the first row is bigger. To be honest this doesn't speed up much.

xm3[tmax_, ni_, a_] := 
 Module[{augvalues = Accumulate[Transpose[mm[10000, a]] + ({#,0} & /@ (0.1 RandomVariate[ParetoDistribution[0.1, ni], 10000]))], j=1}, 
 While[augvalues[[j, 1]] < tmax, j++];  augvalues[[j, 2]] ]

All of this cleans up your code but none of it speeds it up. Most of the time is taken in Ngamm, specifically coming up with 100 different cases of xm2, for each of the values of valt. I can't help wondering if your problem would make more sense to have a single set of random accumulations, and then testing each of the valt values on that same set. Because if the first time 10. is exceeded is say, point 100, then the first time 10^1.2 is exceeded is later, and you can use that fact to speed things up.

Also it's worth testing whether you need an 10000-length sample for each of the 100 samples. If the condition is fulfilled much earlier, you could trim that down, which would definitely save time.

share|improve this answer
Since xx and yy are listable, just using xx[lval, phival] will probably be the fastest method. –  Simon Woods Mar 16 '14 at 13:01
Duh! Of course. Baby brain... thanks @SimonWoods. –  Verbeia Mar 16 '14 at 20:19

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