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

SCM Repository

[diderot] Annotation of /tests/examples/sphere/sphere.diderot
ViewVC logotype

Annotation of /tests/examples/sphere/sphere.diderot

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4694 - (view) (download)

1 : glk 4675 /* ==========================================
2 : glk 4692 Mutually-repelling particles populating a unit sphere
3 : glk 4675
4 : glk 4692 This is heavily based on the [`circle.diderot`](../circle) example; see that
5 :     program for more detailed and explanatory comments. The new things added with
6 :     this example are documented in comments below. A significant capability
7 :     demonstrated here is "population control", whereby particles create new
8 :     particles (using `new`) if there seem to be too few, or die if there are too
9 :     may (using `die`). Along with this there is a new variable `undone` that acts
10 :     as a signal (to the global update method monitoring the computation) that
11 :     things are unfinished because the system population is changing.
12 : glk 4689
13 :     ... TODO: example of running program and looking at output
14 :    
15 : glk 4675 ========================================== */
16 :    
17 :     input vec3{} ipos ("initial positions for all particles") = load("vec3.nrrd");
18 : glk 4689 input real rad ("radius of particle's potential energy support") = 0.1;
19 : glk 4675 input real eps ("system convergence threshold, computed as the coefficient-of-variation of distances to nearest neighbors") = 0.05;
20 : glk 4692 input int pcp ("periodicity of population control (if greater than zero)") = 2;
21 : glk 4683 input real hhInit ("initial step size for potential energy gradient descent") = 1;
22 : glk 4675
23 : glk 4683 real newDist = 0.5*rad; // how far away to birth new particles
24 : glk 4675 real stepMax = rad; // limit on distance to travel per iter
25 :     int iter = 0; // which iteration we're on
26 :    
27 : glk 4692 // Univariate energy functions; see ../circle/circle.diderot for alternatives
28 : glk 4675 function real phi(real r) = (1 - r)^4;
29 :     function real phi'(real r) = -4*(1 - r)^3;
30 : glk 4692 // Energy and force from particle (with radius rad) at vec3 x
31 : glk 4675 function real enr(vec3 x) = phi(|x|/rad);
32 :     function vec3 frc(vec3 x) = phi'(|x|/rad) * (1/rad) * x/|x|; // chain rule
33 :    
34 : glk 4683 // Returns a non-zero vector perpendicular to given non-zero vector v
35 :     function vec3 perp3(vec3 v) {
36 :     int c = 0;
37 :     if (|v[0]| < |v[1]|) {
38 :     c = 1;
39 :     }
40 :     // not v[c] because tensors can only be indexed by constants
41 :     if (|v[1] if 1==c else v[0]| < |v[2]|) {
42 :     c = 2;
43 :     }
44 :     // now c is index of v component with largest absolute value
45 :     vec3 ret = ([v[1] - v[2], -v[0], v[0]] if (c == 0) else
46 :     [-v[1], v[0] - v[2], v[1]] if (c == 1) else
47 :     [-v[2], v[2], v[0] - v[1]]);
48 :     return ret;
49 :     }
50 :    
51 :     // From a given vec3, find a random-ish value uniformly sampling [0,1)
52 :     function real posrnd(vec3 v) {
53 :     vec3 p = 10000*v/rad;
54 :     return fmod(|fmod(p[0],1)| + |fmod(p[1],1)| + |fmod(p[2],1)|, 1);
55 :     }
56 :    
57 :     // Is this an iteration in which to do population control?
58 :     function bool pcIter() = (pcp > 0 && iter > 0 && 0 == iter % pcp);
59 :    
60 :     /* The particle is initialized at position pos0, with initial stepsize hh0.
61 : glk 4686 The first set of particles gets hhInit for the initial stepsize, but new
62 :     particles created by population control benefit from getting the stepsize
63 :     that was adaptively learned by the parent. */
64 : glk 4675 strand particle (vec3 pos0, real hh0) {
65 :     output vec3 pos = pos0/|pos0|;
66 :     real hh = hh0;
67 :     vec3 step = [0,0,0]; // step along force
68 : glk 4692 real closest = rad; // distance to closest neighbor
69 : glk 4686 int ncount = 0; // how many neighbors did we have
70 : glk 4692 /* This "undone" variable signals to global update that something is
71 :     happening or just changed that should delay convergence. In this program
72 :     it is reset to 1 when new particles are created and when there are too
73 :     many neighbors; otherwise it is slowly decreased towards 0. */
74 :     real undone = 1;
75 : glk 4675 update {
76 :     // compute energy and forces on us from neighbors
77 :     real energyLast = 0;
78 :     vec3 force = [0,0,0];
79 : glk 4683 ncount = 0;
80 : glk 4675 foreach (particle P in sphere(rad)) {
81 : glk 4692 vec3 x = P.pos - pos;
82 :     if (|x| == 0) {
83 :     /* we're exactly overlapping with another particle; would be
84 :     nice to have exactly one strand persist and kill the others;
85 :     but simpler to have all overlap-ees die here and let
86 :     population control fill in the hole as needed later */
87 :     die;
88 :     }
89 :     energyLast += enr(x);
90 :     force += frc(x);
91 : glk 4675 ncount += 1;
92 :     }
93 : glk 4683 vec3 norm = normalize(pos); // surface normal for unit circle
94 : glk 4692 if (0 == ncount) {
95 :     if (pcIter()) { // no neighbors, so let's make one
96 :     vec3 npos = pos + newDist*normalize(perp3(norm));
97 :     new particle(npos, hh);
98 :     undone = 1;
99 : glk 4693 } else {
100 :     undone *= pow(0.5, 1.0/(2*(pcp if pcp > 0 else 1)));
101 : glk 4692 }
102 : glk 4686 // set closest to something in case used in global update
103 : glk 4692 closest = newDist;
104 : glk 4675 continue;
105 :     }
106 :    
107 : glk 4683 /* Else we have interacting neighbors; project force onto
108 :     tangent surface, find step, limit step size */
109 : glk 4675 force = (identity[3] - norm⊗norm)•force;
110 :     step = hh*force;
111 :     if (|step| > stepMax) {
112 :     hh *= stepMax/|step|;
113 :     step = hh*force;
114 :     }
115 :    
116 : glk 4683 // Take step, re-apply constraint, find new energy
117 : glk 4675 vec3 posLast = pos;
118 :     pos = normalize(pos + step);
119 :     real energy = 0;
120 :     closest = rad;
121 : glk 4683 ncount = 0;
122 : glk 4675 foreach (particle P in sphere(rad)) {
123 :     energy += enr(P.pos - pos);
124 :     closest = min(closest, |P.pos - pos|);
125 : glk 4683 ncount += 1;
126 : glk 4675 }
127 :    
128 : glk 4683 // Armijo-Goldstein sufficient decrease condition
129 : glk 4675 if (energy - energyLast > 0.5*(pos - posLast)•(-force)) {
130 : glk 4692 hh *= 0.5; // energy didn't decrease as expected: backtrack
131 :     pos = posLast; // try again next iteration
132 :     // no progress, so decrease of undone
133 : glk 4675 } else {
134 :     hh *= 1.02; // opportunistically increase step size
135 : glk 4692 // indicate progress; may be over-written below
136 :     undone *= pow(0.5, 1.0/(2*(pcp if pcp > 0 else 1)));
137 :     // try to have between 5 and 8 neighbors
138 : glk 4683 if (pcIter()) {
139 :     if (ncount < 5) {
140 :     new particle(pos + newDist*normalize(force), hh);
141 : glk 4692 undone = 1;
142 : glk 4683 } else if (ncount > 8) {
143 : glk 4686 /* If this particle has ncount neighbors, then all of those
144 : glk 4683 neighbors probably have a similar number of neighbors. To
145 : glk 4686 get down to having about 8 neighbors, all of them should die
146 :     with a chance of ncount-8 out of ncount. */
147 : glk 4683 if (posrnd(pos) < (ncount - 8.0)/ncount) {
148 :     die;
149 :     }
150 : glk 4692 // else not done if too many neighbors, w/ population control
151 :     undone = 1;
152 : glk 4683 }
153 :     }
154 : glk 4675 }
155 : glk 4683
156 :     // Record actual step taken, in case used in global update
157 :     step = pos - posLast;
158 : glk 4675 }
159 :     }
160 :    
161 :     global {
162 : glk 4687 /* Compute coefficient-of-variation of distance to closest neighbor */
163 : glk 4675 real meancl = mean { P.closest | P in particle.all};
164 :     real varicl = mean { (P.closest - meancl)^2 | P in particle.all};
165 : glk 4686 real covcl = sqrt(varicl) / meancl;
166 : glk 4694 real meanncount = mean { real(P.ncount) | P in particle.all};
167 :     real maxundone = max { P.undone | P in particle.all};
168 : glk 4683 print("(iter ", iter, ") COV(cl)=", covcl,
169 : glk 4692 "; mean(ncount)=", meanncount, "; max(undone)=", maxundone, "\n");
170 : glk 4686
171 : glk 4692 if (covcl < eps // seem to be geometrically uniform
172 :     && maxundone < 0.5) { // and no particle recently set undone=1
173 : glk 4683 print("Stabilizing ", numActive(), " points with COV(closest) ", covcl,
174 : glk 4694 " < ", eps, " and maxundone ", maxundone, " < 0.5 (iter ", iter, ")\n");
175 : glk 4675 stabilize;
176 :     }
177 :     iter += 1;
178 :     }
179 :    
180 :     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