Mathematica 9 is now available
Services & Resources / Wolfram Forums
-----
 /
MathGroup Archive
2004
*January
*February
*March
*April
*May
*June
*July
*August
*September
*October
*November
*December
*Archive Index
*Ask about this page
*Print this page
*Give us feedback
*Sign up for the Wolfram Insider

MathGroup Archive 2004

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

Search the Archive

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


  • Prev by Date: Re: Do-loop conversion
  • Next by Date: Selecting cubic roots in functional form
  • Previous by thread: Do-loop conversion
  • Next by thread: Re: Do-loop conversion