       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,
>
> 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
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:= Timing[{id,res} = fTriangle[10^7];]
Out= {45.857, Null}

Daniel Lichtblau
Wolfram Research

```

• Prev by Date: Re: Mathematica in conjunction with outside program; NMinimize fails.
• Next by Date: Re: Looking for more Mathematica online user groups/forums
• Previous by thread: Help with Speeding up a For loop
• Next by thread: Re: Help with Speeding up a For loop