Re: Do-loop conversion
- To: mathgroup at smc.vnet.net
- Subject: [mg50247] Re: [mg50211] Do-loop conversion
- From: DrBob <drbob at bigfoot.com>
- Date: Mon, 23 Aug 2004 06:34:08 -0400 (EDT)
- References: <200408210704.DAA24758@smc.vnet.net>
- Reply-to: drbob at bigfoot.com
- Sender: owner-wri-mathgroup at wolfram.com
Do loops are not always the wrong way to go and, in this case, the real problem could be this statement:
newc = AppendTo[newc, newdata]
In the first place, it can be replaced with simply
AppendTo[newc, newdata]
The Set statement isn't necessary, since AppendTo already changes the value of its first argument.
Secondly, AppendTo is slow. Consider these three timings:
n=10000
SeedRandom[1]
new1={};
Timing@Do[AppendTo[new1,Random[]],{i,1,n}]
10000
{0.938 Second,Null}
SeedRandom[1]
new2={};
Timing[Do[new2={new2,Random[]},{i,1,n}];new2=Flatten@new2;]
{0.031 Second,Null}
SeedRandom[1]
Timing[new3=Reap[Do[Sow@Random[],{i,1,n}]][[-1,1]];]
{0.016 Second,Null}
new1 == new2 == new2
True
The second method only works for you if "newdata" is a scalar (not a List).
The speed difference gets much MORE dramatic as n gets larger. But, if n=200, there may or may not be a noticeable difference. In that case, other Mathgroupers may have ideas for you. Offhand, I doubt that getting rid of the Do loop will speed things up for you in this case.
Bobby
> Mathematica colleagues,
>
> I have a problem that involves converting a Do-loop structure into
> faster-executing version using functional programming. This is an optics
> problem -- specifically, calculating the angle dependence of the X-ray
> reflectivity of an n-layered medium on an infinitely thick substrate. The
> calculation relies on a well-known recursion relation:
>
> X(j) = [ r(j) + X(j+1) p^2(j+1) ]/[1+ r(j) X(j+1) p^2(j+1)]
>
> where X(j) is the ratio of the reflected and transmitted X-ray amplitudes
> at the bottom of layer j, p(j+1) is a phase factor (given below), and r(j)
> is the Fresnel coefficient for reflection at the interface between layers j
> and j+1:
>
> r(j) = [k(z,j) - k(z,j+1)]/[k(z,j) + k(z,j+1)]
>
> where the wavevector k(z,j) = (2*Pi/wavelength)Sqrt[n^2(j) - Cos(phi)^2],
> n(j) is the complex index of refraction for X-rays for the layer n(j) = 1 -
> delta(j) + I*beta(j), and phi is incident angle of the X-rays. The phase
> factor mentioned above is given by
>
> p(j) = r(j) exp[-2 (k(j) k(j+1)) s(j+1)] with s being the roughness at
> interface j+1.
>
> The recursion layer works because with an infinitely thick substrate, the
> reflected amplitude coming up from the substrate = 0, so at a given angle
> of incidence, you work from the bottom up starting with X(bottom) = 0 and
> add the amplitudes at each interface until you get to the top surface. The
> Mathematica code that works just fine as a Do-loop is:
>
> Do[
> wavevectmap := wavetabc[revdelta, revbeta, phi];
> leftvectmap := leftvectmapc[wavevectmap];
> thtest := thtestcomp[revrough, leftvectmap, wavevectmap];
> maproughrev := roughrevcomp[thtest];
> coeffmap := coeffcomp[wavevectmap, leftvectmap, maproughrev];
> phasemap := phasecomp[wavevectmap, revthick];
> coeffs := fcoeffcomp[coeffmap, phasemap];
> recurr[prior_, {fc_, pf_}] := (fc + (prior*pf^2))/(1 +
> (fc*prior*pf^2));
> intensity := intensitycomp[Fold[recurr, 0, coeffs ]];
> newdata = newdatacomp[phi, intensity, params];
> newc = AppendTo[newc, newdata],
> {phi, phistart, phiend, phiinc}]
>
> The various functions above -- wavetabc, leftvectmapc, thtestcomp,
> roughrevcomp, coeffcomp, phasecomp, fcoeffcomp, intensitycomp, and
> newdatacomp -- are complied functions that take structural parameters for
> each layer (index of refraction parameters delta and beta, layer thickness,
> interfacial roughness) by reading them out of a list {params} that is input
> at the start of the program. I've left out these functions for brevity; if
> needed, I'll post them as well.
>
> As I said, the above Do-loop works just fine by calculating at each angle
> phi the final intensity at that angle, normalizing it with a background and
> maximum peak factor (also in the {params} list and applied by the
> newdatacomp line) and appending it to the intensity datafile
> newc. Executing this loop for a five-layer system and a couple of hundred
> angular points takes around 0.2 seconds -- not bad, except that I need to
> use this in a genetic algorithm minimization routine which needs many
> hundreds of these curves to be calculated with GA- dictated changes to the
> initial parameter matrix. So: if anyone can give me some guidance on how to
> eliminate the Do-loop over the angular range so I can speed things up, I
> would be very grateful.
>
> Thanks in advance,
> Rich Matyi
> rjmatyi at comcast.net
>
>
>
>
>
--
DrBob at bigfoot.com
www.eclecticdreams.net
- References:
- Do-loop conversion
- From: "Rich Matyi" <rjmatyi@comcast.net>
- Do-loop conversion