Home My Page Projects Code Snippets Project Openings diderot

# SCM Repository

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

Tue Jul 18 19:00:35 2017 UTC (2 years 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 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. */

// From a given vec2, find a random-ish value uniformly sampling [0,1)
function real posrnd(vec2 v, real rad) {
return fmod(|fmod(p,1)| + |fmod(p,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;

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, 1),
lerp(-0.5, ISZ-0.5, -1, pp, 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, pos, 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;
}
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;
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,-mno);
a += lerp(-1, 1, gen % 2)*π/6.0;
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
//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, pos, 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], poshInit[ii]], poshInit[ii], ii) | ii in 0 .. length(poshInit)-1 };
```