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