Services & ResourcesWolfram ForumsMathGroup Archive

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

efficiently adding many 2D Gaussians

• To: mathgroup at smc.vnet.net
• Subject: [mg91562] efficiently adding many 2D Gaussians
• From: "M.Roellig" <markus.roellig at googlemail.com>
• Date: Wed, 27 Aug 2008 06:43:36 -0400 (EDT)

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

• Follow-Ups:
  • Re: efficiently adding many 2D Gaussians
    • From: Carl Woll <carlw@wolfram.com>