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.

Improve this question
$\endgroup$
7
  • 2
    $\begingroup$ Related: Randomly packing spheres of fixed radius within a cube $\endgroup$ Commented Sep 30, 2017 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$ Commented Sep 30, 2017 at 14:09
  • $\begingroup$ @becko Note also that none of the proposed solutions seem to properly sample the configurations $\endgroup$ Commented Oct 2, 2017 at 10:58
  • $\begingroup$ What code have you tried? What specific probability distribution and length and number of spheres are of interest? $\endgroup$ Commented Oct 6, 2017 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$ Commented Oct 6, 2017 at 15:26

1 Answer 1

Reset to default
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$

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.