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

SCM Repository

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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 5146 - (download) (annotate)
Fri Jul 21 05:54:36 2017 UTC (22 months, 3 weeks ago) by glk
File size: 13567 byte(s)
adding pbfs1 pbfs2
input image(3)[] img ("data to analyze") = image("pbfs-vol.nrrd");
input real isoval ("isovalue of isosurface to sample") = 0;
input vec3{} ipos ("initial point positions") = load("pbfs-pos.txt");

field#2(3)[] F = bspln3 ⊛ clamp(img);
//field#2(3)[] F = c4hexic ⊛ clamp(img);

/*
function vec3 featureStep(vec3 x) =
  (isoval - F(x))*∇F(x)/(∇F(x)•∇F(x));
function tensor[3,3] featureProj(vec3 x) {
  vec3 norm = normalize(∇F(x));
  return identity[3] - norm⊗norm;
}
function bool featureOut(vec3 x) =
  |∇F(x)| == 0;
*/
/*
function vec3 featureStep(vec3 x) {
  vec3{3} E = evecs(∇⊗∇F(x));
  real{3} L = evals(∇⊗∇F(x));
  return -(1/L{2})*E{2}⊗E{2}•∇F(x);
}
function tensor[3,3] featureProj(vec3 x) {
  vec3{3} E = evecs(∇⊗∇F(x));
  return E{0}⊗E{0} + E{1}⊗E{1};
}
function bool featureOut(vec3 x) {
  real{3} L = evals(∇⊗∇F(x));
  return -L{2}/(|∇F(x)| + 2) < 10;
}
*/
function vec3 featureStep(vec3 x) {
  vec3{3} E = evecs(∇⊗∇F(x));
  real{3} L = evals(∇⊗∇F(x));
  return -(E{1}⊗E{1}/L{1} + E{2}⊗E{2}/L{2})•∇F(x);
}
function tensor[3,3] featureProj(vec3 x) {
  vec3{3} E = evecs(∇⊗∇F(x));
  return E{0}⊗E{0};
}
function bool featureOut(vec3 x) {
  real{3} L = evals(∇⊗∇F(x));
  return -L{1}/(∇F(x)•∇F(x) + 240) < 0.001 /* 0.1 */;
}

input real tipd ("target inter-particle distance") = 0.1;
input real travMax ("max travel allowed to seek feature") = 2;
input int seekMax ("max # steps allowed to seek feature") = 10;
input real seekEps ("conv. threshold on feature seeking") = 0.001;
input real mvmtEps ("conv. threshold on point movement") = 0.01;
input real geoEps ("conv. threshold on system geometry") = 0.05;
input int pcp ("periodicity of population control (PC)") = 5;
input real hhInit ("initial step size for energy descent") = 1;
input bool verb ("verbose") = false; // HIDE

/* Each particle wants between nnmm[0] and nnmm[1] neighbors */
input vec2 nnmm ("desired # neighbors min max") = [6.0, 8.0];

/* Potential function phi(r) has phi(0)=1, phi(r)=0 for r >= 1,
   and a minima at phi'(0.666)=0 and phi(0.666)=-0.005, with C^3
   continuity across the well and with 0 for r >= 1. The
   potential well enables using inter-particle energy to create
   the desired packing density */
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;
// phi'(r) is d phi(r) / dr
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;    // radius of potential well
real rad = tipd/phiWellRad; // actual radius of potential support
function real enr(vec3 x) = phi(|x|/rad); // energy at x
function vec3 frc(vec3 x) = // force on x (via chain rule)
   phi'(|x|/rad) * (1/rad) * x/|x|;
real newDist = tipd;        // how far away to birth
real stepMax = rad;         // limit on travel per iter
int iter = 0;               // current iteration of system
/* how much history matters in averaging movement; higher
   means more memory and more stringent test of convergence */
real mvlerp = 0.5;

/* a number in [0,1) roughly based on low 32 bits of significand
   of given real x. Useful only when compiling with --double */
function real urnd(real x) {
  if (x==0) return 0;
  real l2 = log2(|x|);
  real frxp = 2^(l2 - floor(l2) - 1); // like frexp(x)
  real tmp = fmod(frxp*1048576, 1);
  return fmod(2*tmp, 1);
//  return fmod((frxp-0.5)*1048576,1); // 1048576==2^20; 2097152 == 2^21
}

// generate random-ish value uniformly sampling [0,1)
function real posrnd(vec3 v) {
  //real ret = fmod(urnd(v[0]) + urnd(v[1]) + urnd(v[2]), 1);
  real ret = urnd(v[0]);
  print("66666666 ", v[0], " ", v[1], " ", v[2], " ", ret, "\n");
  return ret;
}

/*
function int posid(vec3 v) {
  vec3 p = 10000000*(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 (PC)?
   If not, pcIter() returns 0. Otherwise, returns 1 when should
   birth new particles, and -1 when should kill then off.  This
   alternation is not due to any language limitations; just a
   strategy for aiding the PC heuristics. */
function int pcIter() = ((iter/pcp)%2)*2 - 1
                        if (pcp > 0 && iter > 0 && 0 == iter % pcp)
                        else 0;

/* each strand seeks the feature (if !found), else interacts
   with other strands to uniformly sample the feature */
strand point (vec3 pos0, real hh0, int pnn) {
  output vec4 posn = [pos0[0], pos0[1], pos0[2], 0];
  vec3 pos = pos0;
  vec3 posLast = pos0; // previous position
  real hh = hh0;       // integration step size
  vec3 step = [0,0,0]; // step towards feature or lower energy
  real closest = rad;  // distance to closest neighbor
  int nn = 0;          // how many neighbors do I have
  bool found = false;  // whether or not I've found the feature
  real trav = 0;       // distance traveled looking for feature
  int steps = 0;       // # steps taken looking for feature
  real mvmt = 1;       // recent movement avg, to test convergence
  int fresh = 1;       // 1: just created, else 0
  update {
    bool vv = false; // pos•[cos(1.658),sin(1.658),0] > 0.82; // HIDE
    real tmp = posrnd(pos0);
    stabilize;
    if (verb && vv) print("hello from ", pos, " (found ", found, ")--------------------------------------\n"); // HIDE
    // can't proceed if outside field or feature isn't there
    if (!inside(pos, F) || featureOut(pos)) {
      if (verb && vv) print("!inside || featureOut DIE\n"); // HIDE
      die;
    }
    posLast = pos;
    if (!found) {
      // various reasons to call it quits
      if ((seekMax > 0 && steps > seekMax)     // too many steps
         || (travMax > 0 && trav > travMax)) { // too much travel
        if (1==pnn /* verb && vv */) { // HIDE
          print("   steps = ", steps, " (vs seekMax ", seekMax, ")\n"); // HIDE
          print("   trav = ", trav, " (vs travMax ", travMax, ")\n"); // HIDE
          print("   featureOut = ", featureOut(pos), "\n"); // HIDE
        } // HIDE
        die;
      }
      if (verb && vv) print("   still here; F=", F(pos), "\n"); // HIDE
      // Take one step towards feature of interest
      step = featureStep(pos);
      pos += step;
      posn = [pos[0], pos[1], pos[2], 0];          
      if (verb && vv) print("  step ", step, " (length ", |step|, ") to ", pos, " where F=", F(pos), "\n"); // HIDE
      // have converged if step is small enough
      if (|step|/rad < seekEps) {
        found = true;
        if (verb && vv) print("   and I have found!\n"); // HIDE
      } else {
        trav += |step|;
        steps += 1;
      }
      // leave mvmt at 1; only change with interactions below
    } else { // feature has been found
      // every path through following should update mvmt
      real energyLast = 0;  // energy at current location
      vec3 force = [0,0,0]; // force on me from neighbors
      nn = 0;
      foreach (point P in sphere(rad)) {
        /* interact with point P even when !P.found to
           avoid disruption from change in P.found */
        vec3 off = P.pos - pos;
        //if (|off| == 0) { die; } // something has gone wrong                       
        energyLast += enr(off);
        force += frc(off);
        nn += 1;
      }
      if (verb && vv) print(" (feat) at ", pos, " where F=", F(pos), ", nn = ", nn, "\n"); // HIDE
      if (0 == nn) {
        /* No neighbors. Can just wait with "continue", or, can
           make a neighbor, but it takes some work to make sure
           that the new position is well-defined for all
           feature codimensions and orientations. */
        vec3 noff0 = featureProj(pos)•[newDist,0,0];
        vec3 noff1 = featureProj(pos)•[0,newDist,0];
        vec3 noff2 = featureProj(pos)•[0,0,newDist];
        vec3 noff = noff0;
        noff = noff if |noff| > |noff1| else noff1;
        noff = noff if |noff| > |noff2| else noff2;
        if (verb && vv) print(" (feat) at ", pos, ": creating new at ", newDist*normalize(noff) + pos, "\n"); // HIDE
        new point(newDist*normalize(noff) + pos, hh, nn);
        mvmt = 1; // reset convergence
        continue;
      }
      /* Else I have interacting neighbors */
      if (verb && vv) print(" (feat) force = ", force); // HIDE
      force = featureProj(pos)•force;
      step = hh*force;
      if (verb && vv) print(" ---project--> ", force, "; step = ", step); // HIDE
      if (|step| > stepMax) { // rescale hh to limit step taken
        hh *= stepMax/|step|;
        step = hh*force;
      }
      if (verb && vv) print(" ---limit--> ", step, "\n"); // HIDE

      // Take step, re-apply constraint, find new energy
      pos += step;
      pos += featureStep(pos);
      pos += featureStep(pos);
      real energy = 0;
      closest = rad;
      nn = 0;
      /* find mean neighbor offset to know (opposite) direction
         in which to add new particles with PC */
      vec3 mno = [0,0,0];
      foreach (point P in sphere(rad)) {
        vec3 off = P.pos - pos;
        energy += enr(off);
        closest = min(closest, |off|);
        mno += off;
        nn += 1;
      }
      mno /= nn;
      posn = [pos[0], pos[1], pos[2], nn];          
      if (verb && vv) print(" (feat) moved to pos ", pos, "; nn = ", nn, "\n"); // HIDE
      if (verb && vv) print(" (feat) energy - energyLast == ", energy, " - ", energyLast, " == ", energy - energyLast, "\n"); // HIDE

      /* Armijo-Goldstein sufficient decrease condition: change
         in energy should roughly match change predicted by force.
         If not, backtrack and try again next iter with smaller hh */
      step = pos - posLast;
      if (energy - energyLast > 0.5*step•(-force)) {
        hh *= 0.5;
        if (verb && vv) print("      OOPS backtrack!  hh == ", hh, "\n"); //HIDE
        pos = posLast;  // try again next iteration
        posn = [pos[0], pos[1], pos[2], nn];          
      } else {
        hh *= 1.1; // opportunistic increase in step size
        if (verb && vv) print("(", pos, ")   progress  hh == ", hh, "; nn = ", nn, "; pcIter = ", pcIter(), "\n"); // HIDE
      }
      if (|step|/rad < seekEps) {
        /* even after backtracking, want to permit PC as long as
           there is small (possibly tentative) motion */
        if (verb && vv) print("(", pos, ") |step|/rad=", |step|/rad, " < ", seekEps, "=seekEps: maybe PC iter=", pcIter(), ", nn=", nn, ", energy=", energy, "\n"); // HIDE
        if (pcIter() > 0 && energy < 0 && nn < nnmm[0]) {
          /* the lower nn is compared to nnmm[0], the more likely
             I'll try to birth a new particle */
          vec3 npos = pos + newDist*normalize(-mno);
          /* don't want new point being killed right away */
          if (inside(npos, F) && !featureOut(npos)) {
            if (verb && vv) print("(", pos, ")   npos=", npos, "; in=", inside(npos,F), "; !fout=", !featureOut(npos), "; posrnd=", posrnd(pos), "\n"); // HIDE
            if (posrnd(pos) < (nnmm[0] - nn)/nnmm[0]) {
              new point(npos, hh, nn);
            }
          }
        } else if (pcIter() < 0 && energy > 0 && nn > nnmm[1]) {
          /* If this particle has nn neighbors, then all those
             neighbors probably have a similar # neighbors. To
             get down to nnmm[1] neighbors, all die with chance
             of nn-nnmm[1] out of nn. */
          if (posrnd(pos) < (nn - nnmm[1])/nn) {
            if (1==pnn) { print("BBB pnn=1 dying\n"); }
            die;
          }
        }
      } // successfully moved to lower energy
      // update movement history
      mvmt = lerp(|pos - posLast|/rad, mvmt, mvlerp);
    } // else found
  } // update
}

global {
  /* use coefficient-of-variation of distance to closest
     neighbor as way of measuring geometric uniformity */
  real meancl = mean { P.closest | P in point.all };
  real varicl = mean { (P.closest - meancl)^2 | P in point.all };
  real covcl = sqrt(varicl) / meancl;
  bool allfound = (min { 1 if P.found else 0 | P in point.all} == 1);
  real maxmv = max { P.mvmt | P in point.all };
  if (verb) { print("\n\n"); } // HIDE
  print("(iter ", iter, " w/ ", numActive(), ")",
        "; allfound=", allfound,
        "; mean(hh)=", mean { P.hh | P in point.all}, //HIDE
        "; mean(cl)=", meancl,
        "; COV(cl)=", covcl,
        "; range(nn)=[", min { P.nn | P in point.all }, ",",
                         max { P.nn | P in point.all }, "]",
        "; max(mvmt)=", maxmv, "\n");
  if (verb) { print("==============================================\n\n"); }  // HIDE
  if (allfound              // all particles have found the feature
      && covcl < geoEps     // is geometrically uniform
      && maxmv < mvmtEps) { // nothing's moving much
    print("Stabilizing ", numActive(), " (iter ", iter, ")",
        "; COV(cl)=", covcl, " < ", geoEps,
        "; max(mvmt)=", maxmv, " < ", mvmtEps, "\n");
      stabilize;
  }
  iter += 1;
}

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

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