Home My Page Projects Code Snippets Project Openings diderot

# SCM Repository

[diderot] Diff of /tests/examples/sphere/sphere.diderot
 [diderot] / tests / examples / sphere / sphere.diderot

# Diff of /tests/examples/sphere/sphere.diderot

revision 4685, Mon Oct 3 16:46:17 2016 UTC revision 4686, Mon Oct 3 16:46:39 2016 UTC
# Line 1  Line 1
1  /* ==========================================  /* ==========================================
2  Mutually-repelling particles on a unit sphere  Mutually-repelling particles on a unit sphere
3
4  This is based on the unit circle example; see that program for more  This is based on the unit circle example; see that program for more detailed
6    documented in comments below. The most significant new capability is
7    "population control", whereby particles create new particles if there are too
8    few (using "new"), or die if there are too may (using "die").
9  ========================================== */  ========================================== */
10
11  input vec3{} ipos ("initial positions for all particles") = load("vec3.nrrd");  input vec3{} ipos ("initial positions for all particles") = load("vec3.nrrd");
# Line 63  Line 66
66  function bool pcIter() = (pcp > 0 && iter > 0 && 0 == iter % pcp);  function bool pcIter() = (pcp > 0 && iter > 0 && 0 == iter % pcp);
67
68  /* The particle is initialized at position pos0, with initial stepsize hh0.  /* The particle is initialized at position pos0, with initial stepsize hh0.
69     The initial set of particles is passed hhInit for initial stepsize, but     The first set of particles gets hhInit for the initial stepsize, but new
70     new particles created by population control benefit from the stepsize that     particles created by population control benefit from getting the stepsize
71     was adaptively learned by the parent */     that was adaptively learned by the parent. */
72  strand particle (vec3 pos0, real hh0) {  strand particle (vec3 pos0, real hh0) {
73     output vec3 pos = pos0/|pos0|;     output vec3 pos = pos0/|pos0|;
74     real hh = hh0;     real hh = hh0;
75     vec3 step = [0,0,0]; // step along force     vec3 step = [0,0,0]; // step along force
76     real closest = 0;    // distance to closest neighbor     real closest = 0;    // distance to closest neighbor
77     int ncount = 0;     int ncount = 0;      // how many neighbors did we have
78     /* This "done" variable is a simple way to signal to the global     /* This "done" variable is a simple way to signal to the global update that
79        update that something has just changed in a way that should        something is happening or just changed in a way that should prevent
80        prevent convergence; in this program it is reset to 0 when new        convergence. In this program it is reset to 0 when new particles are
81        particles are created, and otherwise incrementally increased */        created and when there are too many neighbors; otherwise it is
82          incrementally increased towards 1. */
83     real done = 0;     real done = 0;
84     update {     update {
85        // compute energy and forces on us from neighbors        // compute energy and forces on us from neighbors
# Line 92  Line 96
96           // no neighbors, so let's make one           // no neighbors, so let's make one
97           vec3 npos = normalize(pos + newDist*normalize(perp3(norm)));           vec3 npos = normalize(pos + newDist*normalize(perp3(norm)));
98           new particle(npos, hh);           new particle(npos, hh);
99           // have to set closest to something; used in global update           // set closest to something in case used in global update
100           closest = |npos - pos|;           closest = |npos - pos|;
101           done = 0;           done = 0;
102           continue;           continue;
# Line 122  Line 126
126
127        // Armijo-Goldstein sufficient decrease condition        // Armijo-Goldstein sufficient decrease condition
128        if (energy - energyLast > 0.5*(pos - posLast)•(-force)) {        if (energy - energyLast > 0.5*(pos - posLast)•(-force)) {
129           hh *= 0.5;  // energy didn't decrease as expected, backtrack           hh *= 0.5;  // energy didn't decrease as expected, so backtrack
130           pos = posLast;           pos = posLast;
131        } else {        } else {
132           hh *= 1.02; // opportunistically increase step size           hh *= 1.02; // opportunistically increase step size
# Line 132  Line 136
136                 new particle(pos + newDist*normalize(force), hh);                 new particle(pos + newDist*normalize(force), hh);
137                 done = 0;                 done = 0;
138              } else if (ncount > 8) {              } else if (ncount > 8) {
139                 /* If this particles has ncount neighbors, then all of those                 /* If this particle has ncount neighbors, then all of those
140                    neighbors probably have a similar number of neighbors. To                    neighbors probably have a similar number of neighbors. To
141                    get down to having about 8 neighbors, all of them should                    get down to having about 8 neighbors, all of them should die
142                    die with a chance of ncount-8 out of ncount. */                    with a chance of ncount-8 out of ncount. */
143                 if (posrnd(pos) < (ncount - 8.0)/ncount) {                 if (posrnd(pos) < (ncount - 8.0)/ncount) {
144                    die;                    die;
145                 }                 }
146                 done = 0;                 done = 0; // else not done if too many neighbors
147              }              }
148           }           }
149        }        }
# Line 150  Line 154
154  }  }
155
156  global {  global {
157     /* could test convergence based on movement */     /* Compute coefficient-of-variation of distance to closest neighbor
//real mvmt = max { |P.step|/rad | P in particle.all};
/* This tests convergence based on distance to closest neighbor */
158     real meancl = mean { P.closest | P in particle.all};     real meancl = mean { P.closest | P in particle.all};
159     real varicl = mean { (P.closest - meancl)^2 | P in particle.all};     real varicl = mean { (P.closest - meancl)^2 | P in particle.all};
160     real covcl = sqrt(varicl) / meancl; // coefficient of variation     real covcl = sqrt(varicl) / meancl;
161     real meanncount = mean  { real(P.ncount) | P in particle.all};     real meanncount = mean  { real(P.ncount) | P in particle.all};
162     real mindone = min  { P.done | P in particle.all};     real mindone = min  { P.done | P in particle.all};
163     print("(iter ", iter, ") COV(cl)=", covcl,     print("(iter ", iter, ") COV(cl)=", covcl,
164           "; mean(ncount)=", meanncount, "; min(done)=", mindone, "\n");           "; mean(ncount)=", meanncount, "; min(done)=", mindone, "\n");
165
166     if (covcl < eps          // seem to be geometrically uniform     if (covcl < eps          // seem to be geometrically uniform
167         && mindone > 0.9) {  // and no particle recently set done=0         && mindone > 0.9) {  // and no particle recently set done=0
168        print("Stabilizing ", numActive(), " points with COV(closest) ", covcl,        print("Stabilizing ", numActive(), " points with COV(closest) ", covcl,

Legend:
 Removed from v.4685 changed lines Added in v.4686