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

SCM Repository

[diderot] View of /tests/vis15-bugs/halftone-bug.diderot
ViewVC logotype

View of /tests/vis15-bugs/halftone-bug.diderot

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4705 - (download) (annotate)
Thu Oct 6 18:34:41 2016 UTC (3 years ago) by glk
File size: 5934 byte(s)
simpler to run now
vec2{} ipos = {[0,0],[0,0.1]};
input vec2 radmm ("particle minimum > 0 and maximum radius") = [0.01, 0.1];
input int pcp ("periodicity of population control (if greater than zero)") = 2;
input real hhInit ("initial step size for potential energy gradient descent") = 1;

real newDist = 0.55;  // how far away to birth, as multiple of radius
real stepMax = 1;     // limit on travel per iter, as multiple of radius
int iter = 0;         // which iteration we're on

// Univariate energy functions; see ../circle/circle.diderot for alternatives
// well radius 2/3, depth -0.02
function real phi(real r) =
   1 + r*(-3.87 + r*(4.725 - r*1.8225))  if r < 0.666666666666 else
   2.7 +  r*(-12.96 + r*(22.68 + r*(-17.28 + r*4.86)));
function real phi'(real r) =
   -0.0225*(-2 + 3*r)*(-86 + 81*r)  if r < 0.666666666666 else
   6.48*(r - 1)*(r - 1)*(-2 + r*3);
// Energy and force from particle (with radius rad) at vec2 x
function real enr(vec2 x, real rad) = phi(|x|/rad);
function vec2 frc(vec2 x, real rad) = phi'(|x|/rad) * (1/rad) * x/|x|;

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

/* this one compiles okay
function int pcIter() = 1 if (pcp > 0 && iter > 0 && 0 == iter % pcp) else 0;
*/
/* BUG this one crashes the compiler */
function int pcIter() {
   if (!(pcp > 0 && iter > 0 && 0 == iter % pcp)) {
      return 0;
   }
   return 1 if (iter/pcp) % 2 == 0 else -1;
}

/* The particle is initialized at position pos0, with initial stepsize hh0.
   The first set of particles gets hhInit for the initial stepsize, but new
   particles created by population control benefit from getting the stepsize
   that was adaptively learned by the parent. */
strand particle (vec2 pos0, real hh0) {
   output vec2 pos = pos0;
   real hh = hh0;
   vec2 step = [0,0];   // step along force
   real rad = 0;        // zero isn't a valid value
   real closest = 1;    // distance to closest neighbor (as mult of radius)
   int ncount = 0;      // how many neighbors did we have
   real undone = 1;     // same as in ../sphere
   update {
      if (|pos[0]| > 1 || |pos[1]| > 1) { die; }
      rad = clamp(radmm[0], radmm[1], lerp(radmm[0], radmm[1], lerp(0, 1, -1, pos[0], 1)));
      // compute energy and forces on us from neighbors
      real energyLast = 0;
      vec2 force = [0,0];
      ncount = 0;
      foreach (particle P in sphere(rad)) {
         vec2 x = P.pos - pos;
         if (|x| == 0) {
            die; // same rationale as in ../sphere
         }
         energyLast += enr(x, rad);
         force += frc(x, rad);
         ncount += 1;
      }
      if (0 == ncount) {
         if (pcIter() > 0) {  // no neighbors, so let's make one
            real a = 2*π*posrnd(pos, rad);
            vec2 npos = pos + newDist*rad*[cos(a),sin(a)];
            new particle(npos, hh);
            undone = 1;
         } else {
            undone *= pow(0.5, 1.0/(2*(pcp if pcp > 0 else 1)));
         }
         // set closest to something in case used in global update
         closest = newDist;
         continue;
      }

      /* Else have interacting neighbors; find step, limit step size */
      step = hh*force;
      if (|step| > stepMax*rad) {
         hh *= stepMax*rad/|step|;
         step = hh*force;
      }

      // Take step, re-apply constraint, find new energy
      vec2 posLast = pos;
      pos += step;
      real energy = 0;
      closest = 1;
      ncount = 0;
      vec2 mnd = [0,0]; // mean neighbor difference
      foreach (particle P in sphere(rad)) {
         vec2 nd = P.pos - pos;
         energy += enr(nd, rad);
         closest = min(closest, |nd|/rad);
         mnd += nd;
         ncount += 1;
      }

      // 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;  // try again next iteration
         // no progress, so decrease of undone
      } else {
         hh *= 1.02; // opportunistically increase step size
         // indicate progress; may be over-written below
         undone *= pow(0.5, 1.0/(2*(pcp if pcp > 0 else 1)));
         // try to have between 5 and 8 neighbors
         if (pcIter() != 0) {
            if (ncount < 6) {
               if (pcIter() > 0) {
                  new particle(pos + newDist*rad*normalize(-mnd), hh);
               }
               undone = 1;
            } else if (ncount > 8 && energy > 0 && pcIter() < 0) {
               /* the logic for using posrnd() this way is the same as in
                  ../sphere, but here we can exploit how the potential
                  function has a (negative) well, so energy > 0 means that the
                  total system energy would be lower without this particle */
               if (posrnd(pos, rad) < (ncount - 8.0)/ncount) {
                  die;
               }
               // else not done if too many neighbors, w/ population control
               undone = 1;
            }
         }
      }

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

global {
   /* Compute coefficient-of-variation of 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;
   real meanncount = mean { real(P.ncount) | P in particle.all};
   real maxundone = max { P.undone | P in particle.all};
   print("(iter ", iter, "; pcIter ", pcIter(), ") COV(cl)=", covcl,
         "; mean(cl)=", meancl,
         "; mean(ncount)=", meanncount,
         "; max(undone)=", maxundone, "\n");

   if (iter > 10) {
      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