Is there a better way to write this gather? I have a list of particles, "particlelist", and a list of particle coordinates, "coordlist" and I want to put all the particles into a 3-dimensional ng x ng x ng array of cubic boxes of width d. Each coordinate takes values between xmin and xmin + d x ng. boxIndx contains the correct box for any particle given the particle coordinates. I can then place the appropriate particle into its proper box.

The code below is easy to read but it is a little slow. I know that AppendTo's are supposed to be no-no's but replacing it with Reap and Sow results in even slower times. The test problem typically takes .13-.17 on my Mac laptop. Thanks.

BoxEmUp3[coordlist_, particlelist_, ng_, d_, xmin_] :=
    Module[{ParticleListByBox,particleIndx, boxIndx},

  ParticleListByBox = ConstantArray[{}, ng^3];
   boxIndx = 1 + Floor[(coordlist[[particleIndx, 1 ;; 3]] - xmin)/d].{1, ng, ng^2};

   AppendTo[ParticleListByBox[[boxIndx]], particlelist[[particleIndx]]];,

   {particleIndx, Length[particlelist]}];

xyztst = RandomReal[{-10, 10}, {10^6, 3}];
free = RandomSample[Range[10^6], 10^4];
BoxEmUp3[xyztst[[free]], free, 20, 1., -10] // Timing
  • 2
    $\begingroup$ take a look at BinLists[] $\endgroup$ – Dr. belisarius Apr 3 '14 at 21:37
  • $\begingroup$ It looks to me that BinLists will place the particle coordinates into Boxes but not the particle indexes? I want to store the indexes, not the coordinates. $\endgroup$ – JEP Apr 3 '14 at 21:49
  • $\begingroup$ I had hoped I could use Gather[] for this but I want the empty boxes. A vanilla application of Gather only returns the boxes that have particles in them. I don't know if there is some way to trick Gather into returning the empty ones as well. $\endgroup$ – JEP Apr 3 '14 at 23:37

Although not so faster than yours, yet try next :

NewBoxEmUp3[coordlist_, particlelist_, ng_, d_, xmin_] := 
   Module[{ParticleListByBox, particleIndx, boxIndx, temp},
   ParticleListByBox = ConstantArray[{}, ng^3];
   boxIndx = {1, ng, ng^2}.# & /@ Quotient[coordlist - xmin, d] + 1;
   temp = MapThread[List, {boxIndx, particlelist} ];
   particleIndx = {#[[1, 1]], #[[All, 2]]} & /@ GatherBy[temp, First];
   (ParticleListByBox[[First@#]] = Last@#) & /@ particleIndx;

 xyztst = RandomReal[{-10, 10}, {10^6, 3}];
 free = RandomSample[Range[10^6], 10^4];

 ans1 = BoxEmUp3[xyztst[[free]], free, 20, 1., -10]; // Timing
 ans2 = NewBoxEmUp3[xyztst[[free]], free, 20, 1., -10]; // Timing

 (* {0.109201, Null} *)
 (* {0.031200, Null} *)


 (* True *)
|improve this answer|||||

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.