Re: Help with Speeding up a For loop
- To: mathgroup at smc.vnet.net
- Subject: [mg98933] Re: Help with Speeding up a For loop
- From: "Sjoerd C. de Vries" <sjoerd.c.devries at gmail.com>
- Date: Wed, 22 Apr 2009 05:11:03 -0400 (EDT)
- References: <gsivhv$9f6$1@smc.vnet.net>
The following code does 10 million iterations in 4 minutes:
Ideal = ConstantArray[0, n + 1];
Resolution = ConstantArray[0, n + 1];
For[i = 0, i < n, i++,
If[k[[i]] < Re[(E0 - Eb1[[i]])^2*Sqrt[1 - m^2/(E0 - Eb1[[i]])^2]],
Eb1[[i]] =
RandomReal[NormalDistribution[(Ideal[[i + 1]] = Eb1[[i]]),
Sigma]];
Resolution[[i + 1]] = Eb1[[i]]
];
]
It can be sped up even more, but I think this meets your requirements.
Cheers -- Sjoerd
On Apr 21, 1:13 am, Adam Dally <ada... at wisc.edu> 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 basic=
ally
> a "Monte Carlo Integration" technique. Ideal should simulate the curve.
> Resolution should simulate the curve convoluted with a normal distributio=
n.
> I want to do this for n=10 000 000 or more, but it takes far too long r=
ight
> now. I can do n=100 000 in about 1 minute, but 1 000 000 takes more tha=
n 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; S=
igma = 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]], (*c=
heck
> if the {k,Eb1} value is under the curve*)
> AppendTo[Ideal, Eb1[[i]]]; (*Keep events under curve in 'Id=
eal'*)
> y3 = Eb1[[i]]; (*cast element to a number*)
> Eb1[[i]] = RandomReal[NormalDistribution[y3, Sigma], 1]; (*c=
hoose a
> random number from a normal distribution about that point*)
> AppendTo[Resolution, Eb1[[i]]]; ]] (*Keep that event in 'Resol=
ution'*)