Home My Page Projects Code Snippets Project Openings diderot
Summary Activity Tracker Tasks SCM

SCM Repository

[diderot] View of /tests/examples/sphere/sphere.diderot
ViewVC logotype

View of /tests/examples/sphere/sphere.diderot

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4683 - (download) (annotate)
Mon Oct 3 16:27:34 2016 UTC (2 years, 5 months ago) by glk
File size: 6608 byte(s)
more stable example of particle system with population control, with documentation
/* ==========================================
Mutually-repelling particles on a unit sphere

This is based on the unit circle example; see that program for more
detailed and explanatory comments.
========================================== */

input vec3{} ipos ("initial positions for all particles") = load("vec3.nrrd");
input real rad ("radius of particle's energy support") = 0.1;
input real eps ("system convergence threshold, computed as the coefficient-of-variation of distances to nearest neighbors") = 0.05;
input int pcp ("periodicity of population control (if non-zero)") = 5;
input real hhInit ("initial step size for potential energy gradient descent") = 1;

real newDist = 0.5*rad; // how far away to birth new particles
real stepMax = rad;     // limit on distance to travel per iter
int iter = 0;           // which iteration we're on

/* energy functions */

// phi(0) and phi'(0) are bounded
function real  phi(real r) = (1 - r)^4;
function real phi'(real r) = -4*(1 - r)^3;

/*
// electrostatic potential 1/r, scaled to be C^2 continuous with 0 at r==1
function real  phi(real r) =  (1/r)*(1 - r)^3;
function real phi'(real r) = 3 - 1/(r^2) - 2*r;
*/
/*
// Cotangent-based potential from Meyer et al. SMI'05
function real  phi(real r) = 1/tan(r*π/2) + r*π/2 - π/2;
function real phi'(real r) = (π/2)*(1 - (1/sin(r*π/2))^2);
*/

// Energy and force from particle at vec3 x
function real enr(vec3 x) = phi(|x|/rad);
function vec3 frc(vec3 x) = phi'(|x|/rad) * (1/rad) * x/|x|; // chain rule

// Returns a non-zero vector perpendicular to given non-zero vector v
function vec3 perp3(vec3 v) {
   int c = 0;
   if (|v[0]| < |v[1]|) {
      c = 1;
   }
   // not v[c] because tensors can only be indexed by constants
   if (|v[1] if 1==c else v[0]| < |v[2]|) {
      c = 2;
   }
   // now c is index of v component with largest absolute value
   vec3 ret = ([v[1] - v[2], -v[0], v[0]] if (c == 0) else
               [-v[1], v[0] - v[2], v[1]] if (c == 1) else
               [-v[2], v[2], v[0] - v[1]]);
   return ret;
}

// From a given vec3, find a random-ish value uniformly sampling [0,1)
function real posrnd(vec3 v) {
   vec3 p = 10000*v/rad;
   return fmod(|fmod(p[0],1)| + |fmod(p[1],1)| + |fmod(p[2],1)|, 1);
}

// Is this an iteration in which to do population control?
function bool pcIter() = (pcp > 0 && iter > 0 && 0 == iter % pcp);

/* The particle is initialized at position pos0, with initial stepsize hh0.
   The initial set of particles is passed hhInit for initial stepsize, but
   new particles created by population control benefit from the stepsize that
   was adaptively learned by the parent */
strand particle (vec3 pos0, real hh0) {
   output vec3 pos = pos0/|pos0|;
   real hh = hh0;
   vec3 step = [0,0,0]; // step along force
   real closest = 0;    // distance to closest neighbor
   int ncount = 0;
   /* This "done" variable is a simple way to signal to the global
      update that something has just changed in a way that should
      prevent convergence; in this program it is reset to 0 when new
      particles are created, and otherwise incrementally increased */
   real done = 0;
   update {
      // compute energy and forces on us from neighbors
      real energyLast = 0;
      vec3 force = [0,0,0];
      ncount = 0;
      foreach (particle P in sphere(rad)) {
         energyLast += enr(P.pos - pos);
         force += frc(P.pos - pos);
         ncount += 1;
      }
      vec3 norm = normalize(pos); // surface normal for unit circle
      if (0 == ncount && pcIter()) {
         // no neighbors, so let's make one
         vec3 npos = normalize(pos + newDist*normalize(perp3(norm)));
         new particle(npos, hh);
         // have to set closest to something; used in global update
         closest = |npos - pos|;
         done = 0;
         continue;
      }

      /* Else we have interacting neighbors; project force onto
         tangent surface, find step, limit step size */
      force = (identity[3] - norm⊗norm)•force;
      step = hh*force;
      if (|step| > stepMax) {
         hh *= stepMax/|step|;
         step = hh*force;
      }

      // Take step, re-apply constraint, find new energy
      vec3 posLast = pos;
      pos = normalize(pos + step);
      real energy = 0;
      closest = rad;
      ncount = 0;
      foreach (particle P in sphere(rad)) {
         energy += enr(P.pos - pos);
         closest = min(closest, |P.pos - pos|);
         ncount += 1;
      }
      done = lerp(done, 1, 0.5);  // may be reset below

      // Armijo-Goldstein sufficient decrease condition
      if (energy - energyLast > 0.5*(pos - posLast)•(-force)) {
         hh *= 0.5;  // energy didn't decrease as expected, backtrack
         pos = posLast;
      } else {
         hh *= 1.02; // opportunistically increase step size
         /* try to have between 5 and 8 neighbors */
         if (pcIter()) {
            if (ncount < 5) {
               new particle(pos + newDist*normalize(force), hh);
               done = 0;
            } else if (ncount > 8) {
               /* If this particles has ncount neighbors, then all of those
                  neighbors probably have a similar number of neighbors. To
                  get down to having about 8 neighbors, all of them should
                  die with a chance of ncount-8 out of ncount. */
               if (posrnd(pos) < (ncount - 8.0)/ncount) {
                  die;
               }
               done = 0;
            }
         }
      }

      // Record actual step taken, in case used in global update
      step = pos - posLast;
   }
}

global {
   /* could test convergence based on movement */
   //real mvmt = max { |P.step|/rad | P in particle.all};
   /* This tests convergence based on distance to closest neighbor */
   real meancl = mean { P.closest | P in particle.all};
   real varicl = mean { (P.closest - meancl)^2 | P in particle.all};
   real covcl = sqrt(varicl) / meancl; // coefficient of variation
   real meanncount = mean  { real(P.ncount) | P in particle.all};
   real mindone = min  { P.done | P in particle.all};
   print("(iter ", iter, ") COV(cl)=", covcl,
         "; mean(ncount)=", meanncount, "; min(done)=", mindone, "\n");
   if (covcl < eps          // seem to be geometrically uniform
       && mindone > 0.9) {  // and no particle recently set done=0
      print("Stabilizing ", numActive(), " points with COV(closest) ", covcl,
            " < ", eps, " and mindone > ", mindone, " (iter ", iter, ")\n");
      stabilize;
   }
   iter += 1;
}

initially { particle(ipos[ii], hhInit) | ii in 0 .. length(ipos)-1 };

root@smlnj-gforge.cs.uchicago.edu
ViewVC Help
Powered by ViewVC 1.0.0