Home My Page Projects Code Snippets Project Openings diderot

# SCM Repository

 [diderot] / tests / vis15-bugs / badptclC / pbfsmini.diderot

revision 5431, Wed Sep 27 07:10:41 2017 UTC revision 5432, Wed Sep 27 17:28:52 2017 UTC
# Line 13  Line 13
13    return identity[3] - norm⊗norm;    return identity[3] - norm⊗norm;
14  }  }
15  function bool fMask(vec3 x) = |∇F(x)| > 0;  function bool fMask(vec3 x) = |∇F(x)| > 0;
16  int iter = 0; // HIDE  int iter = 0;
17
18  // Potential energy and force functions  // Potential energy and force functions
19  function real phi(real r) = (1 - r)^4;  function real phi(real r) = (1 - r)^4;
# Line 21  Line 21
21  function real enr(vec3 x) = phi(|x|/rad);  function real enr(vec3 x) = phi(|x|/rad);
23
24  function vec2 jhpos(vec3 x) = [  function vec3 rndrgb(vec3 v) {
25     lerp(-0.5, 700-0.5, -1.6, x[0] /* if x[2] < 0 else x[0]+3.2 */, 1.6),     vec3 p = (100000 + iter)*v/rad;
26     lerp(-0.5, 700-0.5, -1.7, x[1], 1.5)];     return 4*[round(lerp(6, 63, |fmod(p[0],1)|)),
27                 round(lerp(6, 63, |fmod(p[1],1)|)),
28                 round(lerp(6, 63, |fmod(p[2],1)|))];
29    }
30
31    function vec2 jh(vec3 x) = [
32       round(lerp(-0.5, 800-0.5, 0.2, x[0], 0.8)),
33       round(lerp(-0.5, 800-0.5, 0.2, x[1], 0.8))];
34
35    function bool v3eq(vec3 a, vec3 b) = (a[0]==b[0] && a[1]==b[1] && a[2]==b[2]);
36
37  // Strands first find feature, then interact w/ or make neighbors  // Strands first find feature, then interact w/ or make neighbors
38  strand point (int idx, vec3 p0, real hh0, vec3 parentpos_) {  strand point (int idx, vec3 p0, real hh0, vec3 parentID) {
39    output vec3 pos = p0;      // current particle position    output vec3 pos = p0;      // current particle position
40    output vec3 parentpos = parentpos_;    output vec3 ID = rndrgb(p0);
41    real hh = hh0;             // energy gradient descent step size    real hh = hh0;             // energy gradient descent step size
42    vec3 step = [0,0,0];       // total movement this iteration    vec3 step = [0,0,0];       // total movement this iteration
43    bool found = false;        // whether feature has been found    bool found = false;        // whether feature has been found
44    int nstep = 0;             // # steps taken looking for feature    int nstep = 0;             // # steps taken looking for feature
45    int nback = 0;             // # times backtracked after found    int nback = 0;             // # times backtracked after found
46      bool first = true;
47    update {    update {
bool verb = |[468,468] - jhpos(pos)| < 5; // idx == 1287; // HIDE
48      vec3 oldpos = pos;           // save previous position      vec3 oldpos = pos;           // save previous position
49      if (!inside(pos, F)          // if not in field domain,      if (first) {
50          || !fMask(pos)) { die; } //   or close to feature, die        print(ID, "(", iter, ") HELLO at pos=", pos, ":", jh(pos), " born of ", parentID, "\n");
51          first = false;
52        }
53        if (!inside(pos, F) || !fMask(pos)) {
54          print(ID, "(", iter, ") ", "" if found else "!", "found: inside=", inside(pos, F), "; fMask=", fMask(pos), ": DIE\n");
55          die;
56        }
57      if (!found) {                // looking for feature      if (!found) {                // looking for feature
58        step = fStep(pos);         // one step towards feature        step = fStep(pos);         // one step towards feature
59        pos += step;        pos += step;
60        if (verb) print("!found: step=", step, " to pos=", pos, " where F=", F(pos), "\n"); // HIDE        print(ID, "(", iter, ")", " !found: step=", step, " to pos=", pos, ":", jh(pos), " where F=", F(pos), "\n");
61        if (|step|/rad > eps) {    // if took a big step        if (|step|/rad > eps) {    // if took a big step
62          nstep += 1;          nstep += 1;
63          if (nstep > 10) { die; } // too slow          if (nstep > 10) {
64        } else { found = true; }   // else found feature            print(ID, "(", iter, ")", " !found: I'm too slow; DIE\n");
65              die; // too slow
66            }
67          } else {
68            print(ID, "(", iter, ")", " !found: I've found it\n");
69            found = true;
70          }
71      } else {      // feature found; interact with other points      } else {      // feature found; interact with other points
72        real oldE = 0;             // energy at current location        real oldE = 0;             // energy at current location
73        vec3 force = [0,0,0];      // force on me from neighbors        vec3 force = [0,0,0];      // force on me from neighbors
# Line 54  Line 75
75        foreach (point P in sphere(rad)) {        foreach (point P in sphere(rad)) {
76          oldE += enr(P.pos - pos);          oldE += enr(P.pos - pos);
77          force += frc(P.pos - pos);          force += frc(P.pos - pos);
78            print(ID, "(", iter, ")", ": pos=", pos, ":", jh(pos), " forced from ", P.ID, ":", jh(P.pos), " at distance ", |P.pos - pos|, "\n");
79          nn += 1;          nn += 1;
80        }        }
81        if (0 == nn) {             // no neighbors, so make one        if (0 == nn) {             // no neighbors, so make one
82          new point(-1, pos + [0.5*rad,0,0], hh, pos);          print(ID, "(", iter, ")", ": pos=", pos, ":", jh(pos), " but zero neighbors, NEW at pos=", pos + [0.5*rad,0,0], ":", jh(pos + [0.5*rad,0,0]), "\n");
83            new point(-1, pos + [0.5*rad,0,0], hh, ID);
84            print(ID, "(", iter, ")", " " if found else " !", "found: and iter done B\n");
85          continue;          continue;
86        }                          // else interact w/ neighbors        }                          // else interact w/ neighbors
87        vec3 away = hh*fPerp(pos)•force; // away from neighbors        vec3 away = hh*fPerp(pos)•force; // away from neighbors
# Line 65  Line 89
90          away = hh*force;          away = hh*force;
91        }        }
92        if (verb) print("found: pos=", pos, "; F=", F(pos), "; GM=", |∇F(pos)|, "; |fStep|=", |fStep(pos)|, "nn=", nn, "; hh=", hh, "\n"); // HIDE        print(ID, "(", iter, ")", ": pos=", pos, ":", jh(pos), " had nn=", nn, "; oldE=", oldE, " hh=", hh, "\n");
93        pos += away;               // take step along force        pos += away;               // take step along force
94        if (verb) print("      after force step: pos=", pos, "\n"); // HIDE        print(ID, "(", iter, ")", "      after force step: pos=", pos, ":", jh(pos), "\n");
95        pos += fStep(pos);         // move towards feature        pos += fStep(pos);         // move towards feature
96        if (verb) print("           after fStep: pos=", pos, "\n"); // HIDE        //real newE = sum { enr(P.pos - pos) | P in sphere(rad) };
97        real newE = sum { enr(P.pos - pos) | P in sphere(rad) };        real newE = 0;
98        // Armijo condition for predictable energy decrease        foreach (point P in sphere(rad)) {
99        //if (newE < oldE) { // HIDE way too stochastic          newE += enr(P.pos - pos);
100            print(ID, "(", iter, ")", ": pos=", pos, ":", jh(pos), " new energy from ", P.ID, ":", jh(P.pos), " at distance ", |P.pos - pos|, "\n");
101          }
102          print(ID, "(", iter, ")", "      after fStep: pos=", pos, ":", jh(pos), "; newE=", newE, "\n");
103        if (newE - oldE < 0.5*(pos - oldpos)•(-force)) {        if (newE - oldE < 0.5*(pos - oldpos)•(-force)) {
104          hh *= 1.1;               // cautious step size increase          hh *= 1.1;               // cautious step size increase
105          nback = 0;               // reset backtrack counter          nback = 0;               // reset backtrack counter
106          if (nn < 5 && pos[0] > 0 && pos[1] > 0) {          if (nn < 5 && pos[0] > 0 && pos[1] > 0) {
107            if (verb) print("found: pos=", pos, "; nn=", nn, "; NEWNEW!\n"); // HIDE            print(ID, "(", iter, ")", " pos=", pos, ":", jh(pos), "; nn=", nn, "; NEEW at pos=", pos + 0.5*rad*normalize(away), ":", jh(pos + 0.5*rad*normalize(away)), "\n");
108            new point(-10 - nn, pos + 0.5*rad*normalize(away), hh, pos);            new point(-10 - nn, pos + 0.5*rad*normalize(away), hh, ID);
109          }          }
110        } else {        } else {
if (verb) print("           newE|oldE=", newE, "|", oldE, " --> backtrack\n"); // HIDE
111          pos = oldpos; nback += 1; // backtrack          pos = oldpos; nback += 1; // backtrack
112          if (nback > 10) { die; }  // seem to have gotten stuck          print(ID, "(", iter, ")", " newE|oldE=", newE, "|", oldE, " --> backtrack to pos=", pos, ":", jh(pos), "\n");
113            if (nback > 10) {
114               print(ID, "(", iter, ")", " got stuck DIE\n");
115               die;
116            }
117          hh *= 0.5;                // try again w/ smaller step          hh *= 0.5;                // try again w/ smaller step
118        }        }
119        step = pos - oldpos;        // move down force + to feature        step = pos - oldpos;        // move down force + to feature
120          print(ID, "(", iter, ")", " " if found else " !", "found: and iter done C\n");
121      } // else found      } // else found
122    } // update    } // update
123  }  }
124  global {  global {
125    bool allfound = all { P.found | P in point.all};    bool allfound = all { P.found | P in point.all};
126    real maxstep = max { |P.step| | P in point.all };    real maxstep = max { |P.step| | P in point.all };
127    print("(iter ", iter, " w/ ", numActive(), ")", // HIDE    print("================= done with iter ", iter, " w/ ", numActive(),
128          "; allfound=", allfound, // HIDE          "; allfound=", allfound,
129          "; max(mvmt)=", maxstep, // HIDE          "; max(mvmt)=", maxstep,
130          "\n"); // HIDE          "\n");
131    if (allfound && maxstep/rad < eps) { stabilize; }    iter += 1;
132  /*    print("Stabilizing ", numActive(), " (iter ", 0, ")", // HIDE    if (iter == 40) { stabilize; }
"; max(mvmt)=", maxstep, " < ", eps, "\n"); // HIDE */
iter += 1; // HIDE
133  }  }
134  initially { point(ii, ipos[ii], 1, [0,0,0]) | ii in 0 .. length(ipos)-1 };  initially { point(ii, ipos[ii], 1, [0,0,0]) | ii in 0 .. length(ipos)-1 };

Legend:
 Removed from v.5431 changed lines Added in v.5432