       Re: efficiently adding many 2D Gaussians

• To: mathgroup at smc.vnet.net
• Subject: [mg91586] Re: efficiently adding many 2D Gaussians
• From: Ray Koopman <koopman at sfu.ca>
• Date: Thu, 28 Aug 2008 03:18:05 -0400 (EDT)
• References: <g93b4e\$koc\$1@smc.vnet.net>

If the j'th Gaussian is of the form

f[j][x,y] = amplitude[[j]] * Exp[-.5 *
( (x - xmean[[j]])^2 + (y - ymean[[j]])^2 )/variance[[j]] ]

then try  f[x_,y_] := a.Exp[w((x+u)^2 + (y+v)^2)],

where a, u, v, and w are lists:

a = amplitudes, u = -xmeans, v = -ymeans, w = -.5/variances.

On Aug 27, 3:43 am, "M.Roellig" <markus.roellig at googlemail.com> wrote:
> Dear group,
>
> I need some help on reducing run time and memory usage for the
> following problem:
>
> I want to place a large number N (10000+)  of 2D Gaussians in a plane
> (on a large array, m x m ) and add them up.
> I came up with a few approaches but they all fail, or become very slow
> once N >10000 or m>100.
> Here are my ideas so far (not bothering about dimensions, units, etc.
> and using 1000 positions only):
>
> SeedRandom;
> posTab =
>  RandomReal[NormalDistribution[0, 1], {1000, 2}];
>
> 1. idea (fastest I could get):
>
> gsum = Compile[{{x, _Real, 1}},
>    Evaluate[
>     Total[Map[(1.*
>          Exp[-(((x[] - #[])^2 + (x[] - #[])^2)/(
>             2. 0.25^2.))]) &, posTab]]]];
>
> (* Warning Output:
> -----------------------
> Part::partd: Part specification x[] is longer than depth of \
> object. >>
>
> Part::partd: Part specification x[] is longer than depth of \
> object. >>
> ------------------------*)
>
> xyArray = Array[{#1, #2} &, {100, 100}] // MatrixForm;
> map1e4 = Map[gsum, xyArray, {3}]; // Timing
>
> Out= {11.843, Null}
>
> 2. idea (trying to use mapping, and came up with a very inefficient
> solution :-)
>
> Total[
>    Table[Exp[-(((x - #[])^2 + (y - #[])^2)/0.125)], {x, -5,
>        5, .1}, {y, -5, 5, .1}] & /@ posTab]; // Timing
>
> Out= {85.375, Null}
>
> Using more than 1000 Gaussians doesn't fit into my memory and makes
> WinXp swapping memory to disk
> basically making the computation time unbearable. Of course this is no
> big surprise since I create one array for each individual Gaussian
>
>  3. completely different approach hoping to take advantage of internal
> mathematica optimization:
>
> SeedRandom; posTab =
>  RandomReal[NormalDistribution[0, 100], {1000, 2}];
> clumpSum =
>   Total[Map[(1.*
>        Exp[-(((x - #[])^2 + (y - #[])^2)/(2. 10.^2.))]) &,
>     posTab]];
> clumpEm100 =
>    Plot3D[clumpSum, {x, -350, 350}, {y, -350, 350}, PlotRange -> All,
>     PlotPoints -> 100, MaxRecursion -> 0]; // Timing
>
> mapEm100 = clumpEm100[] /. GraphicsComplex[pts_, __] -> pts;
>
> Out= {53.485, Null}
>
> I am grateful for any advice on how to do this faster and more
> efficiently.
> I guess a clever usage of Fold/Nest could do it, but I failed to get
> it done.
>
> The code I use now is a little more complex, since I need to add up
> Gaussians with
> different sizes and amplitudes, which I take from a list:
>
> --------------------
>
> fullEns ={{1/1000, 63095, 374775., 0.00237989, 112.101, 0.818147},
> {1/100,
>   10000, 185962., 0.00647649, 20.143, 2.22645}, {1/10, 1584, 92273.2,
>   0.0176247, 5.91138, 6.05894}, {1, 251, 45785.5, 0.0479629, 0.837552,
>    16.4884}, {10, 39, 22718.5, 0.130523, 0.156455, 44.8706}, {100, 6,
>   11272.8, 0.355198, 0.220978, 122.108}, {1000, 1, 5593.51, 0.966614,
>   0.240513, 332.297}}
>
> (*In each sublist element 2 is the number, element 5 the amplitude,
> element 6 the stand.dev.
> I am discarding the first to rows to reduce total numbers. *)
>
> SeedRandom; posEnsTab =
>  Append[Table[
>    RandomReal[NormalDistribution[0, 100], {fullEns[[j, 2]], 2}], {j,
>     1, Length[fullEns] - 1}], {{0., 0.}}];
> ensSum = Plus @@
>    Table[Total[
>      Map[( fullEns[[j, 5]]*
>          Exp[-(((x - #[])^2 + (y - #[])^2)/(
>             2. fullEns[[j, 6]]^2.))]) &, posEnsTab[[j]]]], {j, 3,
>      Length[fullEns]}];
> ensSum100 =
>    Plot3D[ensSum, {x, -350, 350}, {y, -350, 350}, PlotRange -> All,
>     PlotPoints -> 100, MaxRecursion -> 0]; // Timing
>
> (* Out: {116.453, Null} *)
>
> ensMap100 = ensSum100[] /. GraphicsComplex[pts_, __] -> pts;
> ListPlot3D[ensMap100, PlotRange -> All]
>
> Hope somebody can help me.
>
> Cheers,
>
> Markus

• Prev by Date: Re: Re: Integration Program of Burr Distribution
• Next by Date: Re: How to get result in shape of fraction in Mathematica