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

SCM Repository

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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4918 - (view) (download)

1 : glk 4918 /* ==========================================
2 :     ## circle.diderot: Mutually-repelling particles on a unit circle
3 : glk 4729
4 : glk 4918 This example computes a system of mutually-repelling particles constrained to
5 :     the unit circle. The particles are initialized from a given sequence of a
6 :     2-vectors, and each update tries to move each particle to a lower potential
7 :     energy, as determined by looking at nearby particles. As in the
8 :     [`life.diderot`](../life) example that introduced strand communication, the
9 :     special strand variable named `pos` indicates the particle position, and the
10 :     computation is monitored via a `global` update method. Like
11 :     [`sieve.diderot`](../sieve), this program computes global reductions in the
12 :     global update, used here to compute the coefficient-of-variation of nearest
13 :     inter-particle distances, as the basis of a geometric test of convergence.
14 : glk 4729
15 : glk 4918 Like the [`life.diderot`](../life) example, we use snapshots to inspect
16 :     the system state as it progresses, so we compile with:
17 :    
18 :     diderotc --snapshot --double --exec circle.diderot
19 :    
20 :     When compiling with `--double`, Diderot represents the `real` type
21 :     as a `double` in C++, instead of the default `float`, which produces
22 :     more accurate results.
23 :    
24 :     The following creates a list of `N` randomly oriented unit-length 2-vectors in
25 :     `vec2.nrrd` which will be the initial input. Setting the random number seed
26 :     of `unu 1op nrand` with `-s RNG` ensures reproducible results, or you can
27 :     leave that out to give different results each time.
28 :    
29 :     N=60
30 :     RNG=42
31 :     eval echo {1..$[2*N]} | unu reshape -s 2 $N | unu 1op nrand -s $RNG -o vec2.nrrd
32 :     unu project -i vec2.nrrd -a 0 -m l2 | unu axinsert -a 0 -s 2 | unu 2op / vec2.nrrd - -o vec2.nrrd
33 :    
34 :     Then to run with snapshots saved every iteration (`-s 1`), but limiting the program
35 :     to 500 iterations (`-l 500`), as well as cleaning results from previous run:
36 :    
37 :     rm -f pos-????.nrrd pos.nrrd
38 :     ./circle -s 1 -l 500
39 :    
40 :     If this fails with `ERROR: unexpected arg (or unrecognized flag): "-s"`, it means that
41 :     `circle` wasn't compiled with `--snapshot`, as above. It should converge in under
42 :     500 iterations, producing many `pos-????.nrrd` files, which we can post-process with:
43 :    
44 :     export NRRD_STATE_VERBOSE_IO=0
45 :     for PIIN in pos-????.nrrd; do
46 :     IIN=${PIIN#*-}
47 :     II=${IIN%.*}
48 :     echo "post-processing snapshot $II ... "
49 :     unu dice -i $PIIN -a 0 -o ./
50 :     unu 2op atan2 0.nrrd 1.nrrd | unu histax -a 0 -min -pi -max pi -b 800 -t float -o angle-$II.nrrd
51 :     done
52 :     rm -f 0.nrrd 1.nrrd
53 :     unu join -i angle-*.nrrd -a 1 |
54 :     unu quantize -b 8 -min 0 -max 1 -o angles.png
55 :     rm -f angle-????.nrrd
56 :    
57 :     This produces `angles.png` by unrolling (with atan2) positions along the circle,
58 :     and laying these out along horizontal scanlines, one per iteration. The result
59 :     should look something like [angles-ref.png](angles-ref.png). It is clear that the
60 :     particle interactions made a roughly uniform distribution early on, but
61 :     subsequent refinements took longer. Scrutinizing the top of the image shows
62 :     where strands were stationary either because they repeatedly backtracked in the
63 :     Armijo line search (until they found a step size that led to a sufficient
64 :     energy decrease), or because they had no other particles to interact with.
65 :    
66 :     Some experiments that can be tried with this example:
67 :     * If `--double` is removed from the compilation, the system may not converge to the default `eps` before the iteration limit, due to the decrease in numerical accuracy.
68 :     * There are different possible inter-particle potentials `phi` (uncomment one at a time), which can lead to different dynamics and different convergence speed towards a uniform distribution.
69 :     * If the interaction radius `rad` is decreased while maintining the number of particles, gaps may form.
70 :     * Increasing radius `rad` will tend to force convergence faster, at the expense of computing more particle interactions.
71 :     * There are constants associated with Armijo line search, changing them may lead to faster convergence.
72 :    
73 :     ========================================== */
74 :    
75 :     /* The file containing a sequence (read via "load") currently undergoes no
76 :     compile-time analysis, but there is a run-time check on the shape of the
77 :     array. A sequence of scalars should come from a 1-D array, and a sequence
78 :     of any other kind of value should come from a 2-D array: the faster axis
79 :     should contain the components per value, and the sequence of values is
80 :     along the slower axis. The filename "-" may be used to say that the
81 :     sequence should be read from standard in. A text file (one value per line)
82 :     may be used for non-scalar sequences. */
83 : glk 4917 input vec2{} ipos ("initial positions for all particles") = load("circle-vec2.nrrd");
84 : glk 4729 input real rad ("radius of particle's potential energy support") = 0.25;
85 : glk 4918 input real eps ("system convergence threshold, computed as the coefficient-of-variation of distances to nearest neighbors") = 0.01;
86 : glk 4729
87 :     real hhInit = 1; // initial step size
88 :     real stepMax = rad; // limit on distance to travel per iter
89 :     int iter = 0; // which iteration we're on
90 :    
91 : glk 4918 /* phi defines the potential energy around a particle, as a function of
92 :     distance from particle r > 0. By design, the potential function smoothly
93 :     falls to 0 at r == 1, so that the potential has compact support, which
94 :     simplifies computing the inter-particle interactions. phi'(r) is the
95 :     derivative of phi(r). The ' in phi' is merely suggestive; ' is not acting
96 :     as an operator. Below are some different possibilities for phi and phi';
97 :     only one should be uncommented.
98 :    
99 :     In a conservative force field, force is the *negative* gradient of some
100 :     potential. While phi(r) is the potential at r > 0, due to a particle at 0,
101 :     it can also be considered as the potential at 0, due to a particle at
102 :     r. Moving closer to the particle at r by a small dr > 0, the change in
103 :     potential is dphi = -phi'(r)*dr, which is positive for a monotonically
104 :     decreasing phi(r). The force at 0, or negative gradient of the potential,
105 :     is then phi'(r), which justifies how no negation is used in the frc()
106 :     definition below, based on phi'(). */
107 :    
108 : glk 4729 // phi(0) and phi'(0) are bounded
109 :     function real phi(real r) = (1 - r)^4;
110 :     function real phi'(real r) = -4*(1 - r)^3;
111 :    
112 : glk 4918 /*
113 :     // electrostatic potential 1/r, scaled to be C^2 continuous with 0 at r==1
114 :     function real phi(real r) = (1/r)*(1 - r)^3;
115 :     function real phi'(real r) = 3 - 1/(r^2) - 2*r;
116 :     */
117 :     /*
118 :     // Cotangent-based potential from Meyer et al. SMI'05
119 :     function real phi(real r) = 1/tan(r*π/2) + r*π/2 - π/2;
120 :     function real phi'(real r) = (π/2)*(1 - (1/sin(r*π/2))^2);
121 :     */
122 :    
123 :     /* The enr and frc functions use phi and phi' to define the potential due to,
124 :     and force from, a particle at 2-D offset x with a potential field extending
125 :     to radius rad. */
126 : glk 4729 function real enr(vec2 x) = phi(|x|/rad);
127 :     function vec2 frc(vec2 x) = phi'(|x|/rad) * (1/rad) * x/|x|; // chain rule
128 :    
129 :     /* The particle is initialized at position pos0 */
130 :     strand particle (vec2 pos0) {
131 :     /* These strand variables are visible to the global update and to other
132 :     strands making spatial queries. Any variable inside the scope of
133 : glk 4918 update{} will not be visible in this way. NOTE: "pos" is the special
134 : glk 4729 variable name that *must* be used to enable spatial queries of
135 :     neighboring particles via sphere(). */
136 :     output vec2 pos = pos0/|pos0|;
137 :     real hh = hhInit;
138 :     vec2 step = [0,0]; // step along force
139 :     real closest = 0; // distance to closest neighbor
140 :     int ncount = 0; // how many neighbors did we have
141 :    
142 :     update {
143 :     // Compute energy and forces on us from neighbors
144 :     real energyLast = 0;
145 :     vec2 force = [0,0];
146 :     ncount = 0;
147 :     foreach (particle P in sphere(rad)) {
148 :     energyLast += enr(P.pos - pos);
149 :     force += frc(P.pos - pos);
150 :     ncount += 1;
151 :     }
152 :     if (0 == ncount) {
153 :     /* no neighbors; do nothing to do but set closest to big value */
154 :     closest = rad;
155 :     continue;
156 :     }
157 :    
158 :     // Else we have interacting neighbors
159 :     vec2 norm = normalize(pos); // surface normal for unit circle
160 :     // project force onto tangent surface
161 :     force = (identity[2] - norm⊗norm)•force;
162 :    
163 :     /* Limiting the step size (even before testing for sufficient decrease,
164 :     below) helps keep particles near the surface they are supposed to be
165 :     sampling; this precaution matters more with a non-trivial
166 :     (data-driven) constraint surface */
167 :     step = hh*force; // compute step along force
168 :     if (|step| > stepMax) {
169 :     // decrease hh by factor by which step was too big
170 :     hh *= stepMax/|step|;
171 :     // and find smaller step (of length stepMax)
172 :     step = hh*force;
173 :     }
174 :    
175 :     // take step and re-apply constraint
176 :     vec2 posLast = pos;
177 :     pos = normalize(pos + step);
178 :     // find energy at new location, and distance to closest neighbor
179 :     real energy = 0;
180 :     closest = rad;
181 :     foreach (particle P in sphere(rad)) {
182 :     energy += enr(P.pos - pos);
183 :     closest = min(closest, |P.pos - pos|);
184 :     }
185 :     // save actual step taken
186 :     step = pos - posLast;
187 :    
188 : glk 4918 /* The Armijo-Goldstein sufficient decrease condition: The potential
189 :     gradient, -force, dotted with change in position (pos - posLast),
190 :     should predict the change in potential (energy - energyLast), which
191 :     should be negative. If the potential change is too large (not
192 :     negative enough), decrease step size and try again on the next
193 :     iteration. */
194 : glk 4729 if (energy - energyLast > 0.5*(pos - posLast)•(-force)) {
195 :     hh *= 0.5; // backtrack
196 :     pos = posLast;
197 :     } else {
198 :     hh *= 1.02; // opportunistically increase step size
199 :     }
200 :     }
201 :     }
202 :    
203 : glk 4918 /* In some programs, like ../heron/heron.diderot, strands can individually
204 :     determine when they have converged, and then execute "stabilize". In this
205 :     kind of particle system, the judgement of convergence has to come from a
206 :     more global view of the state of computation. This global update monitors
207 :     the particle system by measuring how uniformly the points are distributed
208 :     along the circle, based on the coefficient of variation of distances to
209 :     nearest neighbors. If the result has converged sufficiently, the
210 :     "stabilize" here results in all active strands being stabilized (and all
211 :     their stabilize{} methods would be called, if they had them). On the other
212 :     hand, if it is clear that the computation cannot converge, global update
213 :     can also call "die". */
214 : glk 4729 global {
215 : glk 4918 // Convergence test could be based on movement
216 :     //real mvmt = max { |P.step|/rad | P in particle.all};
217 :     // This tests convergence based on distances to closest neighbor
218 : glk 4729 real meancl = mean { P.closest | P in particle.all};
219 :     real varicl = mean { (P.closest - meancl)^2 | P in particle.all};
220 :     real covcl = sqrt(varicl) / meancl;
221 :     print("(iter ", iter, ") mean(hh)=", mean { P.hh | P in particle.all},
222 :     "; COV(closest) = ", covcl, "\n");
223 :     if (covcl < eps) {
224 :     print("Stabilizing with COV(closest) ", covcl, " < ", eps, " (iter ", iter, ")\n");
225 :     stabilize;
226 :     }
227 :     iter += 1;
228 :     }
229 :    
230 : glk 4918 /* In this example, there is no "population control" (with "new" or "die"), so
231 :     the set of strands is constant throughout the program execution. It thus
232 :     doesn't matter if the strands are created via "initially {}" (a collection
233 :     of strands) or "initially []" (an array of strands). We create a collection
234 :     of strands here to be consistent with future particle system examples.*/
235 : glk 4729 initially { particle(ipos[ii]) | ii in 0 .. length(ipos)-1 };

root@smlnj-gforge.cs.uchicago.edu
ViewVC Help
Powered by ViewVC 1.0.0