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];
Do[
boxIndx = 1 + Floor[(coordlist[[particleIndx, 1 ;; 3]] - xmin)/d].{1, ng, ng^2};

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

{particleIndx, Length[particlelist]}];
Return[ParticleListByBox];]

xyztst = RandomReal[{-10, 10}, {10^6, 3}];
free = RandomSample[Range[10^6], 10^4];
BoxEmUp3[xyztst[[free]], free, 20, 1., -10] // Timing

• take a look at BinLists[] Apr 3, 2014 at 21:37
• 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.
– JEP
Apr 3, 2014 at 21:49
• 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.
– JEP
Apr 3, 2014 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;
ParticleListByBox]

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} *)

ans1==ans2

(* True *)