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

SCM Repository

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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4704 - (view) (download)

1 : glk 4704
2 :     vec2{} ipos = {[0,0],[0,0.1]};
3 :     input vec2 radmm ("particle minimum > 0 and maximum radius") = [0.01, 0.1];
4 :     input real eps ("system convergence threshold, computed as the coefficient-of-variation of normalized distances to nearest neighbors") = 0.05;
5 :     input int pcp ("periodicity of population control (if greater than zero)") = 2;
6 :     input real hhInit ("initial step size for potential energy gradient descent") = 1;
7 :    
8 :     real newDist = 0.55; // how far away to birth, as multiple of radius
9 :     real stepMax = 1; // limit on travel per iter, as multiple of radius
10 :     int iter = 0; // which iteration we're on
11 :    
12 :     // Univariate energy functions; see ../circle/circle.diderot for alternatives
13 :     // well radius 2/3, depth -0.02
14 :     function real phi(real r) =
15 :     1 + r*(-3.87 + r*(4.725 - r*1.8225)) if r < 0.666666666666 else
16 :     2.7 + r*(-12.96 + r*(22.68 + r*(-17.28 + r*4.86)));
17 :     function real phi'(real r) =
18 :     -0.0225*(-2 + 3*r)*(-86 + 81*r) if r < 0.666666666666 else
19 :     6.48*(r - 1)*(r - 1)*(-2 + r*3);
20 :     // Energy and force from particle (with radius rad) at vec2 x
21 :     function real enr(vec2 x, real rad) = phi(|x|/rad);
22 :     function vec2 frc(vec2 x, real rad) = phi'(|x|/rad) * (1/rad) * x/|x|;
23 :    
24 :     // From a given vec2, find a random-ish value uniformly sampling [0,1)
25 :     function real posrnd(vec2 v, real rad) {
26 :     vec2 p = 10000*v/rad;
27 :     return fmod(|fmod(p[0],1)| + |fmod(p[1],1)|, 1);
28 :     }
29 :    
30 :     /* this one compiles okay
31 :     function int pcIter() = 1 if (pcp > 0 && iter > 0 && 0 == iter % pcp) else 0;
32 :     */
33 :     /* BUG this one crashes the compiler */
34 :     function int pcIter() {
35 :     if (!(pcp > 0 && iter > 0 && 0 == iter % pcp)) {
36 :     return 0;
37 :     }
38 :     return 1 if (iter/pcp) % 2 == 0 else -1;
39 :     }
40 :    
41 :    
42 :     /* The particle is initialized at position pos0, with initial stepsize hh0.
43 :     The first set of particles gets hhInit for the initial stepsize, but new
44 :     particles created by population control benefit from getting the stepsize
45 :     that was adaptively learned by the parent. */
46 :     strand particle (vec2 pos0, real hh0) {
47 :     output vec2 pos = pos0;
48 :     real hh = hh0;
49 :     vec2 step = [0,0]; // step along force
50 :     real rad = 0; // zero isn't a valid value
51 :     real closest = 1; // distance to closest neighbor (as mult of radius)
52 :     int ncount = 0; // how many neighbors did we have
53 :     real undone = 1; // same as in ../sphere
54 :     update {
55 :     if (|pos[0]| > 1 || |pos[1]| > 1) { die; }
56 :     rad = clamp(radmm[0], radmm[1], lerp(radmm[0], radmm[1], lerp(0, 1, -1, pos[0], 1)));
57 :     // compute energy and forces on us from neighbors
58 :     real energyLast = 0;
59 :     vec2 force = [0,0];
60 :     ncount = 0;
61 :     foreach (particle P in sphere(rad)) {
62 :     vec2 x = P.pos - pos;
63 :     if (|x| == 0) {
64 :     die; // same rationale as in ../sphere
65 :     }
66 :     energyLast += enr(x, rad);
67 :     force += frc(x, rad);
68 :     ncount += 1;
69 :     }
70 :     if (0 == ncount) {
71 :     if (pcIter() > 0) { // no neighbors, so let's make one
72 :     real a = 2*π*posrnd(pos, rad);
73 :     vec2 npos = pos + newDist*rad*[cos(a),sin(a)];
74 :     new particle(npos, hh);
75 :     undone = 1;
76 :     } else {
77 :     undone *= pow(0.5, 1.0/(2*(pcp if pcp > 0 else 1)));
78 :     }
79 :     // set closest to something in case used in global update
80 :     closest = newDist;
81 :     continue;
82 :     }
83 :    
84 :     /* Else have interacting neighbors; find step, limit step size */
85 :     step = hh*force;
86 :     if (|step| > stepMax*rad) {
87 :     hh *= stepMax*rad/|step|;
88 :     step = hh*force;
89 :     }
90 :    
91 :     // Take step, re-apply constraint, find new energy
92 :     vec2 posLast = pos;
93 :     pos += step;
94 :     real energy = 0;
95 :     closest = 1;
96 :     ncount = 0;
97 :     vec2 mnd = [0,0]; // mean neighbor difference
98 :     foreach (particle P in sphere(rad)) {
99 :     vec2 nd = P.pos - pos;
100 :     energy += enr(nd, rad);
101 :     closest = min(closest, |nd|/rad);
102 :     mnd += nd;
103 :     ncount += 1;
104 :     }
105 :    
106 :     // Armijo-Goldstein sufficient decrease condition
107 :     if (energy - energyLast > 0.5*(pos - posLast)•(-force)) {
108 :     hh *= 0.5; // energy didn't decrease as expected: backtrack
109 :     pos = posLast; // try again next iteration
110 :     // no progress, so decrease of undone
111 :     } else {
112 :     hh *= 1.02; // opportunistically increase step size
113 :     // indicate progress; may be over-written below
114 :     undone *= pow(0.5, 1.0/(2*(pcp if pcp > 0 else 1)));
115 :     // try to have between 5 and 8 neighbors
116 :     if (pcIter() != 0) {
117 :     if (ncount < 6) {
118 :     if (pcIter() > 0) {
119 :     new particle(pos + newDist*rad*normalize(-mnd), hh);
120 :     }
121 :     undone = 1;
122 :     } else if (ncount > 8 && energy > 0 && pcIter() < 0) {
123 :     /* the logic for using posrnd() this way is the same as in
124 :     ../sphere, but here we can exploit how the potential
125 :     function has a (negative) well, so energy > 0 means that the
126 :     total system energy would be lower without this particle */
127 :     if (posrnd(pos, rad) < (ncount - 8.0)/ncount) {
128 :     die;
129 :     }
130 :     // else not done if too many neighbors, w/ population control
131 :     undone = 1;
132 :     }
133 :     }
134 :     }
135 :    
136 :     // Record actual step taken, in case used in global update
137 :     step = pos - posLast;
138 :     }
139 :     }
140 :    
141 :     global {
142 :     /* Compute coefficient-of-variation of distance to closest neighbor */
143 :     real meancl = mean { P.closest | P in particle.all};
144 :     real varicl = mean { (P.closest - meancl)^2 | P in particle.all};
145 :     real covcl = sqrt(varicl) / meancl;
146 :     real meanncount = mean { real(P.ncount) | P in particle.all};
147 :     real maxundone = max { P.undone | P in particle.all};
148 :     print("(iter ", iter, "; pcIter ", pcIter(), ") COV(cl)=", covcl,
149 :     "; mean(cl)=", meancl,
150 :     "; mean(ncount)=", meanncount,
151 :     "; max(undone)=", maxundone, "\n");
152 :    
153 :     if (covcl < eps // seem to be geometrically uniform
154 :     && maxundone < 0.5) { // and no particle recently set undone=1
155 :     print("Stabilizing ", numActive(), " points with COV(closest) ", covcl,
156 :     " < ", eps, " and maxundone ", maxundone, " < 0.5 (iter ", iter, ")\n");
157 :     stabilize;
158 :     }
159 :     iter += 1;
160 :     }
161 :    
162 :     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