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

SCM Repository

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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 5111 - (view) (download)

1 : jhr 5111 #version 1
2 :    
3 : glk 5106 input image(2)[] img ("field in which to disperse particles") = image("posh-img.nrrd");
4 :     field#0(2)[] F = tent ⊛ clamp(img);
5 :    
6 :     input vec3{} poshInit ("initial positions and stepsize for all particles");
7 :     input vec2 radmm ("particle minimum > 0 and maximum radius") = [0.01, 0.1];
8 :     input real eps ("system convergence threshold, computed as mean movement (normalized by particle radius)") = 0.005;
9 :     input int pcp ("periodicity of population control (if greater than zero)") = 2;
10 :     input real hhInit ("initial step size for potential energy gradient descent") = 1;
11 :     input real well ("strength of parabolic well") = 0;
12 :     input int ISZ ("size of later debugging image") = 700;
13 :    
14 :     /* A potential function phi(r) found via Mathematica (see findingPhi.nb in
15 :     this directory): this is the lowest-degree piecewise polynomial phi(r) with
16 :     phi(0)=1 and a potential well at phi(0.666)=-0.005, with C^3 continuity
17 :     across the well and with 0 for r >= 1. Having a slight potential well is
18 :     essential for the strategy, used by this program, of using inter-particle
19 :     energy to create the desired packing density, rather than the interaction
20 :     of particles and some ambient image-based potential field. */
21 :     function real phi(real r) =
22 :     1 + r*(-4.248582222661734 + r*(5.991287388535527 + (-2.864766048887071 + 0.06790595504541913*r)*r))
23 :     if r < 0.666 else
24 :     -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)
25 :     if r < 1 else 0;
26 :     function real phi'(real r) =
27 :     -4.248582222661734 + r*(11.982574777071054 + (-8.594298146661213 + 0.2716238201816765*r)*r)
28 :     if r < 0.666 else
29 :     -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))))
30 :     if r < 1 else 0;
31 :     real phiWellRad = 0.666;
32 :    
33 :     real newDist = phiWellRad; // how far away to birth, as multiple of radius
34 :     real stepMax = 1; // limit on travel per iter, as multiple of radius
35 :     int iter = 0; // which iteration we're on
36 :    
37 :     /* Energy and force from particle (with radius rad) at vec2 x. Unlike the
38 :     sphere and circle examples, this takes the radius as an argument. */
39 :     function real enr(vec2 x, real rad) = phi(|x|/rad);
40 :     function vec2 frc(vec2 x, real rad) = phi'(|x|/rad) * (1/rad) * x/|x|;
41 :    
42 :     // From a given vec2, find a random-ish value uniformly sampling [0,1)
43 :     function real posrnd(vec2 v, real rad) {
44 :     vec2 p = 10000*v/rad;
45 :     return fmod(|fmod(p[0],1)| + |fmod(p[1],1)|, 1);
46 :     }
47 :    
48 :     /* Is this an iteration in which to do population control? If not, pcIter()
49 :     returns 0. Otherwise, it returns 1 if we should birth new particles, and -1
50 :     if we should kill off particles. This alternation is not because of any
51 :     language limitations; it is just a strategy for aiding the population
52 :     control heuristics. */
53 :     function int pcIter() = ((iter/pcp)%2)*2 - 1
54 :     if (pcp > 0 && iter > 0 && 0 == iter % pcp)
55 :     else 0;
56 :    
57 :     function real radius(real ff) {
58 :     real ret = min(radmm[1], radmm[0]*sqrt(1/clamp(0.000000001, 1, ff)))/phiWellRad;
59 :     real rrr = 0.450450417;
60 :     if (ret != rrr) {
61 :     print("radius HEYHEY: ", ff, " -> ", ret, (" > " if (ret > rrr) else " < " if (ret < rrr) else " == "), rrr, "\n");
62 :     }
63 :     /* there is probably another BUG in here:
64 :     the print above never happens, but replacing the next line with
65 :     "return rrr" changes the program behavior */
66 :     return ret;
67 :     }
68 :    
69 :     function vec2 ipos(vec2 pp) =
70 :     [lerp(-0.5, ISZ-0.5, -1, pp[0], 1),
71 :     lerp(-0.5, ISZ-0.5, -1, pp[1], 1)];
72 :    
73 :     /* The particle is initialized at position pos0, with initial stepsize hh0.
74 :     The first set of particles gets hhInit for the initial stepsize, but new
75 :     particles created by population control benefit from getting the stepsize
76 :     that was adaptively learned by the parent. */
77 :     strand particle (int gen, vec2 pos0, real hh0, int idx0) {
78 :     int iter = 0;
79 :     int idx = idx0;
80 :     vec2 pos = pos0;
81 :     vec2 opos1 = [0,0];
82 :     vec2 opos2 = [0,0];
83 :     real hh = hh0;
84 :     output vec4 posh = [pos[0], pos[1], hh, 0];
85 :     vec2 step = [0,0]; // step along force
86 :     real rad = 0; // zero isn't a valid value
87 :     real closest = 1; // distance to closest neighbor (as mult of radius)
88 :     int ncount = 0; // how many neighbors did we have
89 :     real mvmt = 1;
90 :     real energy=0;
91 :     bool vv = false;
92 :     update {
93 :     /* Limit the system to only work inside the image domain. */
94 :     if (!inside(pos, F)) {
95 :     die;
96 :     }
97 :     // rad is nominal radius for current position
98 :     rad = radius(F(pos));
99 :     vv = (idx==19);
100 :     if (vv) { print("(", iter, ":", idx, ") rad=", rad, " at pos=", pos, "/ipos=", ipos(pos), "\n"); }
101 :     real energyLast = 0;
102 :     vec2 force = [0,0];
103 :     ncount = 0;
104 :     foreach (particle P in sphere(rad)) {
105 :     vec2 x = P.pos - pos;
106 :     real rr = radius(F(lerp(pos, P.pos, 0.5)));
107 :     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"); }
108 :     energyLast += enr(x, rr);
109 :     force += frc(x, rr);
110 :     ncount += 1 if |x| < rr else 0;
111 :     if (vv) { print(" |x|=", |x|, "<" if |x| < rr else ">=", rr, " ==> ncount now ", ncount, "\n"); }
112 :     }
113 :     energyLast += well*(pos•pos)^2;
114 :     force -= well*4*pos*(pos•pos);
115 :     /* Else we have interacting neighbors; find step, limit step size */
116 :     step = hh*force;
117 :     if (|step| > stepMax*rad) {
118 :     hh *= stepMax*rad/|step|;
119 :     step = hh*force;
120 :     }
121 :    
122 :     // Take step, find new energy
123 :     vec2 posLast = pos;
124 :     pos += step;
125 :     energy = 0;
126 :     closest = 1;
127 :     ncount = 0;
128 :     vec2 mno = [0,0];
129 :     foreach (particle P in sphere(rad)) {
130 :     real rr = radius(F(lerp(pos, P.pos, 0.5)));
131 :     vec2 x = P.pos - pos;
132 :     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"); }
133 :     if (|x| < rr) {
134 :     energy += enr(x, rr);
135 :     closest = min(closest, |x|/rr);
136 :     mno += x;
137 :     ncount += 1;
138 :     if (vv) { print(" |x|=", |x|, "<" if |x| < rr else ">=", rr, " ==> ncount now ", ncount, "\n"); }
139 :     }
140 :     }
141 :     energy += well*(pos•pos)^2;
142 :     mno /= ncount;
143 :    
144 :     real osci = 0.666;
145 :     if (vv) { print("(", iter, ":", idx, ") found new energy=", energy, " with ncount=", ncount, "\n"); }
146 :     // Armijo-Goldstein sufficient decrease condition
147 :     if (energy - energyLast > 0.5*(pos - posLast)•(-force)) {
148 :     hh *= 0.5; // energy didn't decrease as expected: backtrack
149 :     pos = posLast; // try again next iteration
150 :     if (vv) { print("(", iter, ":", idx, ") OOPS backtrack \n"); }
151 :     // don't change opos1, opos2
152 :     } else {
153 :     // do change opos1, opos2
154 :     opos2 = opos1; opos1 = posLast;
155 :     osci = |opos1 - pos|/|opos2 - pos|;
156 :     hh *= 1.1;
157 :     if (vv) { print("(", iter, ":", idx, ") progress with osci=", osci, " \n"); }
158 :     if (pcIter() > 0) {
159 :     if (energy < 0 && ncount < 6
160 :     && posrnd(pos, rad) < (6.0 - ncount)/6.0) {
161 :     real a = atan2(-mno[1],-mno[0]);
162 :     a += lerp(-1, 1, gen % 2)*π/6.0;
163 :     vec2 noff = newDist*rad*[cos(a),sin(a)];
164 :     vec2 npos = pos + noff;
165 :     if (inside(npos, F) && inside(pos + 2*noff, F)) {
166 :     new particle(gen+1, npos, hh, -1);
167 :     }
168 :     }
169 :     }
170 :     }
171 :    
172 :     /* Record information that may be used in next global update,
173 :     or by other strands in the next iteration */
174 :     step = pos - posLast;
175 :     mvmt = lerp(mvmt, |step|/rad, 0.5); // should decay to 0
176 :     rad = radius(F(pos));
177 :     //print("iter= ", iter, " idx= ", idx, " posh= ", posh, "\n");
178 :     if (vv) {
179 :     print("(", iter, ":", idx, ") DONE: rad=", rad, ", osci=", osci, ", @pos=", pos, "/ipos=", ipos(pos), "]\n");
180 :     }
181 :     posh = [pos[0], pos[1], hh, osci];
182 :     iter += 1;
183 :     }
184 :     }
185 :    
186 :     global {
187 :     real totenr = sum { P.energy | P in particle.all};
188 :     real meanmv = mean { P.mvmt | P in particle.all };
189 :     print("=-=-= iter ", iter,
190 :     "; tot#=", numActive(),
191 :     "; totenr=", totenr,
192 :     "; mean(hh)=", mean { P.hh | P in particle.all},
193 :     "; mean(ncount)=", mean { real(P.ncount) | P in particle.all},
194 :     "; mean(mvmt)=", meanmv,
195 :     "\n");
196 :     /* The variation of inter-particle distances seems to decrease the
197 :     reliability of COV(closest) as a convergence indicator, so unlike
198 :     ../sphere this instead uses a particle motion test */
199 :     if (meanmv < eps) {
200 :     print("Stabilizing\n");
201 :     print("Stabilizing\n");
202 :     print("Stabilizing\n");
203 :     print("Stabilizing ", numActive(), " points with mean(movement) ", meanmv,
204 :     " < ", eps, " (iter ", iter, ")\n");
205 :     print("Stabilizing\n");
206 :     print("Stabilizing\n");
207 :     print("Stabilizing\n");
208 :     stabilize;
209 :     }
210 :     iter += 1;
211 :     }
212 :    
213 :     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