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[1111111];
> posTab =
> RandomReal[NormalDistribution[0, 1], {1000, 2}];
>
> 1. idea (fastest I could get):
>
> gsum = Compile[{{x, _Real, 1}},
> Evaluate[
> Total[Map[(1.*
> Exp[-(((x[[1]] - #[[1]])^2 + (x[[2]] - #[[2]])^2)/(
> 2. 0.25^2.))]) &, posTab]]]];
>
> (* Warning Output:
> -----------------------
> Part::partd: Part specification x[[1]] is longer than depth of \
> object. >>
>
> Part::partd: Part specification x[[2]] is longer than depth of \
> object. >>
> ------------------------*)
>
> xyArray = Array[{#1, #2} &, {100, 100}] // MatrixForm;
> map1e4 = Map[gsum, xyArray, {3}]; // Timing
>
> Out[22]= {11.843, Null}
>
> 2. idea (trying to use mapping, and came up with a very inefficient
> solution :-)
>
> Total[
> Table[Exp[-(((x - #[[1]])^2 + (y - #[[2]])^2)/0.125)], {x, -5,
> 5, .1}, {y, -5, 5, .1}] & /@ posTab]; // Timing
>
> Out[9]= {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
> before adding them up.
>
> 3. completely different approach hoping to take advantage of internal
> mathematica optimization:
>
> SeedRandom[1111111]; posTab =
> RandomReal[NormalDistribution[0, 100], {1000, 2}];
> clumpSum =
> Total[Map[(1.*
> Exp[-(((x - #[[1]])^2 + (y - #[[2]])^2)/(2. 10.^2.))]) &,
> posTab]];
> clumpEm100 =
> Plot3D[clumpSum, {x, -350, 350}, {y, -350, 350}, PlotRange -> All,
> PlotPoints -> 100, MaxRecursion -> 0]; // Timing
>
> mapEm100 = clumpEm100[[1]] /. GraphicsComplex[pts_, __] -> pts;
>
> Out[25]= {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[1111111]; 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 - #[[1]])^2 + (y - #[[2]])^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[[1]] /. GraphicsComplex[pts_, __] -> pts;
> ListPlot3D[ensMap100, PlotRange -> All]
>
> Hope somebody can help me.
>
> Cheers,
>
> Markus