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

SCM Repository

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

View of /tests/vis15-bugs/posh.diderot

Parent Directory Parent Directory | Revision Log Revision Log


Revision 5111 - (download) (annotate)
Tue Jul 18 19:00:35 2017 UTC (22 months, 1 week ago) by jhr
File size: 9312 byte(s)
  added #version directive
#version 1

input image(2)[] img ("field in which to disperse particles") = image("posh-img.nrrd");
field#0(2)[] F = tent ⊛ clamp(img);

input vec3{} poshInit ("initial positions and stepsize for all particles");
input vec2 radmm ("particle minimum > 0 and maximum radius") = [0.01, 0.1];
input real eps ("system convergence threshold, computed as mean movement (normalized by particle radius)") = 0.005;
input int pcp ("periodicity of population control (if greater than zero)") = 2;
input real hhInit ("initial step size for potential energy gradient descent") = 1;
input real well ("strength of parabolic well") = 0;
input int ISZ ("size of later debugging image") = 700;

/* A potential function phi(r) found via Mathematica (see findingPhi.nb in
   this directory): this is the lowest-degree piecewise polynomial phi(r) with
   phi(0)=1 and a potential well at phi(0.666)=-0.005, with C^3 continuity
   across the well and with 0 for r >= 1. Having a slight potential well is
   essential for the strategy, used by this program, of using inter-particle
   energy to create the desired packing density, rather than the interaction
   of particles and some ambient image-based potential field. */
function real phi(real r) =
1 + r*(-4.248582222661734 + r*(5.991287388535527 + (-2.864766048887071 + 0.06790595504541913*r)*r))
if r < 0.666 else
-0.005 + ((0.4482053856358993 + (-2.6838645846460842 + (6.026642031390917 - 4.811690244623482*(-0.666 + r))*(-0.666 + r))*(-0.666 + r))*(-0.666 + r))*(-0.666 + r)
if r < 1 else 0;
function real phi'(real r) =
-4.248582222661734 + r*(11.982574777071054 + (-8.594298146661213 + 0.2716238201816765*r)*r)
if r < 0.666 else
-15.42591894092919 + 0.4482053856358993*(-1.332 + 2*r) + r*(-2.6838645846460842*(-3.996 + 3*r) + 6.026642031390917*(5.322672 + r*(-7.992 + 4*r)) - 4.811690244623482*(-5.90816592 + r*(13.30668 + r*(-13.32 + 5*r))))
if r < 1 else 0;
real phiWellRad = 0.666;

real newDist = phiWellRad; // 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

/* Energy and force from particle (with radius rad) at vec2 x. Unlike the
   sphere and circle examples, this takes the radius as an argument. */
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);
}

/* Is this an iteration in which to do population control?  If not, pcIter()
   returns 0. Otherwise, it returns 1 if we should birth new particles, and -1
   if we should kill off particles. This alternation is not because of any
   language limitations; it is just a strategy for aiding the population
   control heuristics. */
function int pcIter() = ((iter/pcp)%2)*2 - 1
                        if (pcp > 0 && iter > 0 && 0 == iter % pcp)
                        else 0;

function real radius(real ff) {
   real ret = min(radmm[1], radmm[0]*sqrt(1/clamp(0.000000001, 1, ff)))/phiWellRad;
   real rrr = 0.450450417;
   if (ret != rrr) {
      print("radius HEYHEY: ", ff, " -> ", ret, (" > " if (ret > rrr) else " < " if (ret < rrr) else " == "), rrr, "\n");
   }
   /* there is probably another BUG in here:
      the print above never happens, but replacing the next line with
      "return rrr" changes the program behavior */
   return ret;
}

function vec2 ipos(vec2 pp) =
  [lerp(-0.5, ISZ-0.5, -1, pp[0], 1),
   lerp(-0.5, ISZ-0.5, -1, pp[1], 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 (int gen, vec2 pos0, real hh0, int idx0) {
   int iter = 0;
   int idx = idx0;
   vec2 pos = pos0;
   vec2 opos1 = [0,0];
   vec2 opos2 = [0,0];
   real hh = hh0;
   output vec4 posh = [pos[0], pos[1], hh, 0];
   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 mvmt = 1;
   real energy=0;
   bool vv = false;
   update {
      /* Limit the system to only work inside the image domain. */
      if (!inside(pos, F)) {
         die;
      }
      // rad is nominal radius for current position
      rad = radius(F(pos));
      vv = (idx==19);
      if (vv) { print("(", iter, ":", idx, ") rad=", rad, " at pos=", pos, "/ipos=", ipos(pos), "\n"); }
      real energyLast = 0;
      vec2 force = [0,0];
      ncount = 0;
      foreach (particle P in sphere(rad)) {
         vec2 x = P.pos - pos;
         real rr = radius(F(lerp(pos, P.pos, 0.5)));
         if (vv) { print("(", iter, ":", idx, ") AA seeing P.idx=", P.idx, " at P.pos=", P.pos, "/ipos=", ipos(P.pos), " --> x=", x, "; rr=", rr, ", |x|/rad=", |x|/rad, "\n"); }
         energyLast += enr(x, rr);
         force += frc(x, rr);
         ncount += 1 if |x| < rr else 0;
         if (vv) { print("              |x|=", |x|, "<" if |x| < rr else ">=", rr, " ==> ncount now ", ncount, "\n"); }
      }
      energyLast += well*(pos•pos)^2;
      force -= well*4*pos*(pos•pos);
      /* Else we have interacting neighbors; find step, limit step size */
      step = hh*force;
      if (|step| > stepMax*rad) {
         hh *= stepMax*rad/|step|;
         step = hh*force;
      }

      // Take step, find new energy
      vec2 posLast = pos;
      pos += step;
      energy = 0;
      closest = 1;
      ncount = 0;
      vec2 mno = [0,0];
      foreach (particle P in sphere(rad)) {
         real rr = radius(F(lerp(pos, P.pos, 0.5)));
         vec2 x = P.pos - pos;
         if (vv) { print("(", iter, ":", idx, ") BB seeing P.idx=", P.idx, " at P.pos=", P.pos, "/ipos=", ipos(P.pos), " --> x=", x, "; rr=", rr, ", |x|/rad=", |x|/rad, "\n"); }
         if (|x| < rr) {
            energy += enr(x, rr);
            closest = min(closest, |x|/rr);
            mno += x;
            ncount += 1;
            if (vv) { print("              |x|=", |x|, "<" if |x| < rr else ">=", rr, " ==> ncount now ", ncount, "\n"); }
         }
      }
      energy += well*(pos•pos)^2;
      mno /= ncount;

      real osci = 0.666;
      if (vv) { print("(", iter, ":", idx, ") found new energy=", energy, " with ncount=", ncount, "\n"); }
      // 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
         if (vv) { print("(", iter, ":", idx, ") OOPS backtrack \n"); }
         // don't change opos1, opos2
      } else {
         // do change opos1, opos2
         opos2 = opos1; opos1 = posLast;
         osci = |opos1 - pos|/|opos2 - pos|;
         hh *= 1.1;
         if (vv) { print("(", iter, ":", idx, ") progress with osci=", osci, " \n"); }
         if (pcIter() > 0) {
            if (energy < 0 && ncount < 6
                && posrnd(pos, rad) < (6.0 - ncount)/6.0) {
               real a = atan2(-mno[1],-mno[0]);
               a += lerp(-1, 1, gen % 2)*π/6.0;
               vec2 noff = newDist*rad*[cos(a),sin(a)];
               vec2 npos = pos + noff;
               if (inside(npos, F) && inside(pos + 2*noff, F)) {
                  new particle(gen+1, npos, hh, -1);
               }
            }
         }
      }

      /* Record information that may be used in next global update,
         or by other strands in the next iteration */
      step = pos - posLast;
      mvmt = lerp(mvmt, |step|/rad, 0.5); // should decay to 0
      rad = radius(F(pos));
      //print("iter= ", iter, " idx= ", idx, " posh= ", posh, "\n");
      if (vv) {
         print("(", iter, ":", idx, ") DONE: rad=", rad, ", osci=", osci, ", @pos=", pos, "/ipos=", ipos(pos), "]\n");
      }
      posh = [pos[0], pos[1], hh, osci];
      iter += 1;
   }
}

global {
   real totenr = sum { P.energy | P in particle.all};
   real meanmv = mean { P.mvmt | P in particle.all };
   print("=-=-= iter ", iter,
         "; tot#=", numActive(),
         "; totenr=", totenr,
         "; mean(hh)=", mean { P.hh | P in particle.all},
         "; mean(ncount)=", mean { real(P.ncount) | P in particle.all},
         "; mean(mvmt)=", meanmv,
         "\n");
   /* The variation of inter-particle distances seems to decrease the
      reliability of COV(closest) as a convergence indicator, so unlike
      ../sphere this instead uses a particle motion test */
   if (meanmv < eps) {
      print("Stabilizing\n");
      print("Stabilizing\n");
      print("Stabilizing\n");
      print("Stabilizing ", numActive(), " points with mean(movement) ", meanmv,
            " < ", eps, " (iter ", iter, ")\n");
      print("Stabilizing\n");
      print("Stabilizing\n");
      print("Stabilizing\n");
      stabilize;
   }
   iter += 1;
}

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

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