Services & ResourcesWolfram ForumsMathGroup Archive

[Date Index] [Thread Index] [Author Index]

Re: Help with Speeding up a For loop

• To: mathgroup at smc.vnet.net
• Subject: [mg98927] Re: [mg98881] Help with Speeding up a For loop
• From: Daniel Lichtblau <danl at wolfram.com>
• Date: Wed, 22 Apr 2009 05:09:56 -0400 (EDT)
• References: <200904202313.TAA09729@smc.vnet.net>

Adam Dally wrote:
> I am using an Intel MacBook with OS X 10.5.6.
> 
> I am trying to create 2 lists: "Ideal" and "Resolution".  This is basically
> a "Monte Carlo Integration" technique. Ideal should simulate the curve.
> Resolution should simulate the curve convoluted with a normal distribution.
> I want to do this for n=10 000 000 or more, but it takes far too long right
> now. I can do n=100 000 in about 1 minute, but 1 000 000 takes more than an
> hour. I haven't waited long enough for 10 000 000 to finish (it has been 5
> days).
> 
> Thank you,
> Adam Dally
> 
> Here is the code:
> 
> ClearAll[E0, Eb1, m, DeltaE, Sigma, k, n, y3, Ideal, Resolution, i,
> normalizer, maxE, minE]
> Eb1 = 0; k = 0; n = 10000; E0 = 2470; m = 0.2; DeltaE = 50; Sigma = 5; maxE
> = E0 - m; minE = E0 - DeltaE; Resolution = {Eb1}; Ideal = {Eb1};  (*Setup
> all constants, lists and ranges*)
> 
> Eb1 = RandomReal[{minE, maxE}, n];  (*create a list of 'n' random Eb1
> values*)
> k = -RandomReal[TriangularDistribution[{-2470, 0}, -0.1], n]; (*create a
> list of 'n' random k values; triangle distribution gives more successful
> results*)
> 
> For[i = 0, i < n, i++,
> 
>  If[k[[i]] < Re[(E0 - Eb1[[i]])^2*Sqrt[1 - m^2/(E0 - Eb1[[i]])^2]], (*check
> if the {k,Eb1} value is under the curve*)
>      AppendTo[Ideal, Eb1[[i]]];  (*Keep events under curve in 'Ideal'*)
>      y3 = Eb1[[i]]; (*cast element to a number*)
>      Eb1[[i]] = RandomReal[NormalDistribution[y3, Sigma], 1]; (*choose a
> random number from a normal distribution about that point*)
>      AppendTo[Resolution, Eb1[[i]]]; ]] (*Keep that event in 'Resolution'*)
> 

Best to avoid AppendTo for large list constructions (it will give 
quadratic complexity). Instead use either Sow/Reap or list nesting 
followed by Flatten.

In this particular example, one can also just preinitialize arrays to 
full length, change elements as needed, and remove the extras before 
returning. The code below does this and uses Compile, hence is fairly fast.

randomRealNormal[eb1j_,sigma_] :=
   RandomReal[NormalDistribution[eb1j,sigma]]

randomRealTriangular[lo_,hi_,max_,n_] :=
   RandomReal[TriangularDistribution[{lo,hi},max], n]

fTriangle = Compile[{{n,_Integer}}, Module[
   {ideal, resolution, eb1, e0=2470., m=.2, deltaE=50.,
     sigma=5., maxE, minE, kvec, i=1, eb1j, eminus},
	{maxE,minE} = e0-{m,deltaE};
	ideal = Table[0.,{n}];
	resolution = Table[0.,{n}];
	eb1 = RandomReal[{minE,maxE}, n];
	kvec = -randomRealTriangular[-e0, 0., -.1, n];
	Do [
	  eb1j = eb1[[j]];
	  eminus = (e0 - eb1j)^2;
	  If [kvec[[j]] < eminus*Sqrt[1-m^2/eminus],
	    i++;
		ideal[[i]] = eb1j;
		eb1j = randomRealNormal[eb1j,sigma];
		resolution[[i]] = eb1j;
		];
	  , {j,n}];
	ideal = Take[ideal,i];
	resolution = Take[resolution,i];
	{ideal,resolution}
   ],
   {{randomRealTriangular[___], _Real, 1},
   {randomRealNormal[___], _Real}
   }];

In[9]:= Timing[{id,res} = fTriangle[10^7];]
Out[9]= {45.857, Null}


Daniel Lichtblau
Wolfram Research

• References:
  • Help with Speeding up a For loop
    • From: Adam Dally <adally@wisc.edu>