Re: efficiently adding many 2D Gaussians
- To: mathgroup at smc.vnet.net
- Subject: [mg91590] Re: efficiently adding many 2D Gaussians
- From: Jens-Peer Kuska <kuska at informatik.uni-leipzig.de>
- Date: Thu, 28 Aug 2008 03:18:49 -0400 (EDT)
- References: <g93b4e$koc$1@smc.vnet.net>
Hi, look at http://www.umiacs.umd.edu/~yangcj/papers/siam_fgt_v11.pdf or search "Fast Gauss transform" with Google Regards Jens 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 >