Re: efficiently adding many 2D Gaussians
- To: mathgroup at smc.vnet.net
- Subject: [mg91572] Re: [mg91562] efficiently adding many 2D Gaussians
- From: Carl Woll <carlw at wolfram.com>
- Date: Thu, 28 Aug 2008 03:15:28 -0400 (EDT)
- References: <200808271043.GAA21263@smc.vnet.net>
M.Roellig 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
>
>
Perhaps the following will help:
SeedRandom[1111111];
posTab = RandomReal[NormalDistribution[0, 1], {2000, 2}];
In[3]:= MemoryInUse[]
MaxMemoryUsed[]
Out[3]= 5898936
Out[4]= 7225904
In[5]:= Table[Total[Exp[-((x - posTab[[All, 1]])^2 + (y - posTab[[All,
2]])^2)/.125]], {x, -5, 5, .1}, {y, -5, 5, .1}]; // Timing
Out[5]= {10.031,Null}
In[6]:= MemoryInUse[]
MaxMemoryUsed[]
Out[6]= 7098408
Out[7]= 7237848
Carl Woll
Wolfram Research
- References:
- efficiently adding many 2D Gaussians
- From: "M.Roellig" <markus.roellig@googlemail.com>
- efficiently adding many 2D Gaussians