Home My Page Projects Code Snippets Project Openings diderot

# SCM Repository

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

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

revision 4695, Tue Oct 4 18:40:43 2016 UTC revision 4696, Tue Oct 4 19:26:05 2016 UTC
# Line 71  Line 71
71  // Is this an iteration in which to do population control?  // Is this an iteration in which to do population control?
72  function bool pcIter() = (pcp > 0 && iter > 0 && 0 == iter % pcp);  function bool pcIter() = (pcp > 0 && iter > 0 && 0 == iter % pcp);
73
74  strand particle (int gen, vec3 pos0, real hh0) {  function vec3 featureStep(vec3 x) = -normalize(∇F(x)) * F(x)/|∇F(x)|;
75
76    strand particle (vec3 pos0, real hh0) {
77     output vec3 pos = pos0;     output vec3 pos = pos0;
78     vec3 posLast = pos0;     vec3 posLast = pos0;
79     real hh = hh0;     real hh = hh0;
# Line 84  Line 86
86     real undone = 1;       // same as in ../sphere     real undone = 1;       // same as in ../sphere
87     update {     update {
88        posLast = pos;        posLast = pos;
89        if (gen > 0) { print("hello from ", pos, " (isoFound ", isoFound, ")--------------\n"); }        if (false) print("hello from ", pos, " (isoFound ", isoFound, ")--------------\n");
90        if (!isoFound) {        if (!isoFound) {
91           // various reasons to call it quits           // various reasons to call it quits
92           if (gen > 0) {           if (false) {
93           print("   inside = ", inside(pos, F), "\n");           print("   inside = ", inside(pos, F), "\n");
94           print("   isoSteps = ", isoSteps, " (vs isoStepsMax ", isoStepsMax, ")\n");           print("   isoSteps = ", isoSteps, " (vs isoStepsMax ", isoStepsMax, ")\n");
95           print("   isoTravel = ", isoTravel, " (vs isoTravelMax ", isoTravelMax, ")\n");           print("   isoTravel = ", isoTravel, " (vs isoTravelMax ", isoTravelMax, ")\n");
# Line 101  Line 103
103               || |∇F(pos)| == 0) {             // can't compute step               || |∇F(pos)| == 0) {             // can't compute step
104              die;              die;
105           }           }
106           if (gen > 0) print("   still here; F=", F(pos), "\n");           if (false) print("   still here; F=", F(pos), "\n");
107           // The Newton-Raphson step           // Take one step towards feature of interest
108           step = -normalize(∇F(pos)) * F(pos)/|∇F(pos)|;           step = featureStep(pos);
109           pos += step;           pos += step;
110           if (gen > 0) print("   step ", step, " (length ", |step|, ") to ", pos, " where F=", F(pos), "\n");           if (false) print("   step ", step, " (length ", |step|, ") to ", pos, " where F=", F(pos), "\n");
111           // We've converged if step is small enough           // We've converged if step is small enough
112           if (|step| < isoEps) {           if (|step| < isoEps) {
113              isoFound = true;              isoFound = true;
114              if (gen > 0) print("   and we have isoFound!\n");              if (false) print("   and we have isoFound!\n");
115           } else {           } else {
117              isoSteps += 1;              isoSteps += 1;
# Line 126  Line 128
128              force += frc(x);              force += frc(x);
129              ncount += 1;              ncount += 1;
130           }           }
131           if (gen > 0) print(" (iso) at ", pos, " where F=", F(pos), ", ncount = ", ncount, "\n");           if (false) print(" (iso) at ", pos, " where F=", F(pos), ", ncount = ", ncount, "\n");
132           vec3 norm = normalize(pos); // surface normal for unit circle           vec3 norm = normalize(∇F(pos));
133           if (0 == ncount) {           if (0 == ncount) {
134              if (pcIter()) { // no neighbors, so let's make one              if (pcIter()) { // no neighbors, so let's make one
135                 vec3 npos = pos + newDist*normalize(perp3(norm));                 vec3 npos = pos + newDist*normalize(perp3(norm));
136                 print(" (iso) from pos ", pos, " where F=", F(pos), ", delF=", ∇F(pos), ": NEW at ", npos, "\n");                 if (false) print(" (iso) from pos ", pos, " where F=", F(pos), ", delF=", ∇F(pos), ": NEW at ", npos, "\n");
137                 new particle(gen+1, npos, hh);                 new particle(npos, hh);
138                 undone = 1;                 undone = 1;
139              } else {              } else {
140                 undone *= pow(0.5, 1.0/(2*(pcp if pcp > 0 else 1)));                 undone *= pow(0.5, 1.0/(2*(pcp if pcp > 0 else 1)));
# Line 144  Line 146
146
147           /* Else we have interacting neighbors; project force onto           /* Else we have interacting neighbors; project force onto
148              tangent surface, find step, limit step size */              tangent surface, find step, limit step size */
149           if (gen > 0) print(" (iso) force = ", force);           if (false) print(" (iso) force = ", force);
150           force = (identity[3] - norm⊗norm)•force;           force = (identity[3] - norm⊗norm)•force;
151           step = hh*force;           step = hh*force;
152           if (gen > 0) print(" ---project--> ", force, "; step = ", step);           if (false) print(" ---project--> ", force, "; step = ", step);
153           if (|step| > stepMax) {           if (|step| > stepMax) {
154              hh *= stepMax/|step|;              hh *= stepMax/|step|;
155              step = hh*force;              step = hh*force;
156           }           }
157           if (gen > 0) print(" ---limit--> ", step, "\n");           if (false) print(" ---limit--> ", step, "\n");
158
159           // Take step, re-apply constraint, find new energy           // Take step, re-apply constraint, find new energy
160           pos = normalize(pos + step);           pos += step;
161             pos += featureStep(pos);
162             pos += featureStep(pos);
163             pos += featureStep(pos);
164           real energy = 0;           real energy = 0;
166           ncount = 0;           ncount = 0;
# Line 164  Line 169
169              closest = min(closest, |P.pos - pos|);              closest = min(closest, |P.pos - pos|);
170              ncount += 1;              ncount += 1;
171           }           }
172            print(" (iso) moved to pos ", pos, "; ncount = ", ncount, "\n");           if (false) print(" (iso) moved to pos ", pos, "; ncount = ", ncount, "\n");
173             if (false) print(" (iso) energy - energyLast == ", energy, " - ", energyLast, " == ", energy - energyLast, "\n");
174
175           // Armijo-Goldstein sufficient decrease condition           // Armijo-Goldstein sufficient decrease condition
176           if (energy - energyLast > 0.3*(pos - posLast)•(-force)) {           if (energy - energyLast > 0.5*(pos - posLast)•(-force)) {
177              hh *= 0.5;  // energy didn't decrease as expected: backtrack              hh *= 0.5;  // energy didn't decrease as expected: backtrack
178                if (false) print("      OOPS backtrack!  hh == ", hh, "\n");
179              pos = posLast;  // try again next iteration              pos = posLast;  // try again next iteration
180              // no progress, so decrease of undone              // no progress, so decrease of undone
181           } else {           } else {
182              hh *= 1.02; // opportunistically increase step size              hh *= 1.02; // opportunistically increase step size
183              // indicate progress; may be over-written below              // indicate progress; may be over-written below
184              undone *= pow(0.5, 1.0/(2*(pcp if pcp > 0 else 1)));              undone *= pow(0.5, 1.0/(2*(pcp if pcp > 0 else 1)));
185                if (false) print("      progress  hh == ", hh, "; undone = ", undone, "; ncount = ", ncount, "; pcIter = ", pcIter(), "\n");
186              // try to have between 5 and 8 neighbors              // try to have between 5 and 8 neighbors
187              if (pcIter()) {              if (pcIter()) {
188                   if (false) print("      considering pop cntl (iter ", iter, ", pcp ", pcp, ") with ncount ", ncount, "\n");
189                 if (ncount < 5) {                 if (ncount < 5) {
190                    // new particle(gen +1, pos + newDist*normalize(force), hh);                    new particle(pos + newDist*normalize(force), hh);
191                    undone = 1;                    undone = 1;
192                 } else if (ncount > 8) {                 } else if (ncount > 8) {
193                    if (posrnd(pos) < (ncount - 8.0)/ncount) {                    if (posrnd(pos) < (ncount - 8.0)/ncount) {
# Line 203  Line 212
212     real meanncount = mean { real(P.ncount) | P in particle.all};     real meanncount = mean { real(P.ncount) | P in particle.all};
213     real maxundone = max { P.undone | P in particle.all};     real maxundone = max { P.undone | P in particle.all};
214     bool allfound = (min { 1 if P.isoFound else 0 | P in particle.all} == 1);     bool allfound = (min { 1 if P.isoFound else 0 | P in particle.all} == 1);
215     print("\n\n= = = = (iter ", iter, " w/ ", numStrands(), " or ", numActive(), ") mean(isofound) = ",     if (false) print("\n\n\n\n");
216       print("(iter ", iter, " w/ ", numStrands(), " or ", numActive(), ") mean(isofound) = ",
217           mean { 1.0 if P.isoFound else 0.0 | P in particle.all },           mean { 1.0 if P.isoFound else 0.0 | P in particle.all },
218           "; allfound=", allfound,           "; allfound=", allfound,
219           "; mean(hh)=", mean { P.hh | P in particle.all},           "; mean(hh)=", mean { P.hh | P in particle.all},
220           "; COV(cl)=", covcl,           "; COV(cl)=", covcl,
221           "; mean(ncount)=", meanncount, "; max(undone)=", maxundone, "\n\n\n\n");           "; mean(ncount)=", meanncount, "; max(undone)=", maxundone, "\n");
222       if (false) print("\n\n\n\n");
223
224     if (covcl < geoEps        // seem to be geometrically uniform     if (covcl < geoEps        // seem to be geometrically uniform
225         && allfound           // and all found the isosurface         && allfound           // and all found the isosurface
# Line 220  Line 231
231     iter += 1;     iter += 1;
232  }  }
233
234  initially { particle(1,[cent[0] + lerp(-extent[0], extent[0], 0, ii, sz0-1)/2,  initially { particle([cent[0] + lerp(-extent[0], extent[0], 0, ii, sz0-1)/2,
235                        cent[1] + lerp(-extent[1], extent[1], 0, jj, sz1-1)/2,                        cent[1] + lerp(-extent[1], extent[1], 0, jj, sz1-1)/2,
236                        cent[2] + lerp(-extent[2], extent[2], 0, kk, sz2-1)/2],                        cent[2] + lerp(-extent[2], extent[2], 0, kk, sz2-1)/2],
237                       hhInit)                       hhInit)

Legend:
 Removed from v.4695 changed lines Added in v.4696