3
$\begingroup$

Suppose I have a distribution of sphere radii, $P(r)$. I need to sample configurations of $N$ spheres in a cubic box of side $L$ (assume they fit), where the $N$ spheres have radii following the distribution $P(r)$.

I know Mathematica has some tools to simulate random processes, but I don't have much experience with them. Perhaps someone can suggest a way of doing this? I just need the idea, I can look into the details.

$\endgroup$
  • 2
    $\begingroup$ Related: Randomly packing spheres of fixed radius within a cube $\endgroup$ – C. E. Sep 30 '17 at 12:31
  • $\begingroup$ @C.E. That question is a particular case of this one, with $P(r) = \delta(r - 0.7)$. However, since that question did not any solution exploiting Mathematica functionality (they are just loops), that indicates that I'm better off doing this simulation in another language. $\endgroup$ – becko Sep 30 '17 at 14:09
  • $\begingroup$ @becko Note also that none of the proposed solutions seem to properly sample the configurations $\endgroup$ – Lukas Lang Oct 2 '17 at 10:58
  • $\begingroup$ What code have you tried? What specific probability distribution and length and number of spheres are of interest? $\endgroup$ – Daniel Lichtblau Oct 6 '17 at 15:12
  • $\begingroup$ @DanielLichtblau To be honest, after asking here and seeing the related question linked by C.E., I gave up and coded this in Julia. The probability distribution can vary. You have something in mind for a specific probabitlity distribution? $\endgroup$ – becko Oct 6 '17 at 15:26
4
$\begingroup$

Okay, here goes.

Given a distro, box length, and number of desired points one can generate a bunch more, remove the ones that are too close to the edges to fit, then iterate through what remains to discard neighbors that are too close. Only remove neighbors indexed higher (since a neighbor too close but indexed lower would have led to the removal of the higher indexed sphere already).

Here is an example. I'll start with 3x as many spheres as I want.

len = 10;
dist = GammaDistribution[.3, 2];
num = 200;
mult = 3;
SeedRandom[1111]
radii = RandomVariate[dist, mult*num];
centers = RandomReal[len, {mult*num, 3}];
spheres = Partition[Riffle[centers, radii], 2];
spheres2 = 
  Select[spheres, 
   Min[Abs[#[[1]] - #[[2]]]] > 0 && 
     Max[Abs[#[[1]] + #[[2]]]] < len &];

After removing the ones too close to the box boundary, how many remain?

numspheres = Length[spheres2]

(* Out[503]= 513 *)

What is the largest radius?

maxrad = Max[newradii]

(* Out[505]= 4.50929047049 *)

The iteration can be made fairly fast using a NearestFunction. In the search phase, for a given sphere we only need find neighbors within current sphere radius + maximal sphere radius, so we use this to make for a faster search.

This is pedestrian procedural code and maybe could be made faster with Scan or Compile or ...

AbsoluteTiming[newradii = spheres2[[All, 2]];
 newcenters = spheres2[[All, 1]];
 nf = Nearest[newcenters -> "Index"];
 present = ConstantArray[True, numspheres];
 Do[
  If[! present[[j]], Continue[]];
  nbrs = Rest[nf[newcenters[[j]], {Infinity, maxrad + newradii[[j]]}]];
  Do[
   If[! present[[k]], Continue[]];
   If[Norm[newcenters[[j]] - newcenters[[k]]] < 
     newradii[[j]] + newradii[[k]],
    present[[k]] = False];
   , {k, nbrs}
   ];
  , {j, numspheres}
  ];
 spheres3 = Pick[spheres2, present];]

(* Out[507]= {0.224918, Null} *)

Now how many do we have?

Length[spheres3]

(* Out[509]= 374 *)

Grab what we need.

finalspheres = Take[spheres3, UpTo[num]];

We'll have a look at them.

Graphics3D[Map[Sphere[#[[1]], #[[2]]] &, finalspheres]]

Same code for 500 took around 1.3 seconds. I would not expect this to work too well when the max radius gets large relative to the dimensions of the box. As to how well this gives the correct distribution, given the manner of excluding, I do not have the expertise to address that.

enter image description here

--- edit ---

Actually I did try this with `Compile and got around a 10x factor as follows. First some preprocessing, which is quite fast.

AbsoluteTiming[maxrad = Max[newradii];
 newradii = spheres2[[All, 2]];
 newcenters = spheres2[[All, 1]];
 nf = Nearest[newcenters -> "Index"];]

(* Out[571]= {0.000296, Null} *)

Now we use Compile and allow for external evaluations of the nearest function. I did not useQuiet` and got the warnings shown below, but it compiles and runs fine. I guess the type inferencing is less than ideal. Also I had to explicitly set a couple of initializations just so the inferencing would not be fooled at run time and punt to uncompiled evaluation.

removal = 
  Compile[{{newcenters, _Real, 2}, {newradii, _Real, 
     1}, {maxrad, _Real}}, Module[
    {numspheres = Length[newradii], present, nbrs = {1, 1}, nbr = 1},
    present = ConstantArray[1, numspheres];
    Do[
     If[0 == present[[j]], Continue[]];
     nbrs = 
      Rest[nf[newcenters[[j]], {Infinity, maxrad + newradii[[j]]}]];
     Do[
      nbr = nbrs[[k]];
      If[0 == present[[nbr]], Continue[]];
      If[Norm[newcenters[[j]] - newcenters[[nbr]]] < 
        newradii[[j]] + newradii[[nbr]],
       present[[nbr]] = 0];
      , {k, Length[nbrs]}
      ];
     , {j, numspheres}
     ];
    present], RuntimeOptions -> "Speed", CompilationTarget -> "C"];

During evaluation of In[576]:= Compile::cset: Variable nbrs of type {_Integer,1} encountered in assignment of type {_Real,1}.

During evaluation of In[576]:= Compile::cset: Variable nbrs of type {_Integer,1} encountered in assignment of type {_Real,1}.

Let's try it out. This is with the number of spheres set to 500 so spheres3 is not the same is in the earlier part of this response.

AbsoluteTiming[
 sph3 = Pick[spheres2, removal[newcenters, newradii, maxrad], 1];]
sph3 === spheres3

(* Out[577]= {0.144264, Null}

Out[578]= True *)

When we want 1000 spheres with these parameter settings, the uncompiled version takes around 7.7 seconds on my desktop and the compiled version is .35 for a factor of 22 improvement.

--- end edit ---

| improve this answer | |
$\endgroup$

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.