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

# SCM Repository

[diderot] Annotation of /tests/lamont-tests/unit-circle-pc3.diderot
 [diderot] / tests / lamont-tests / unit-circle-pc3.diderot

# Annotation of /tests/lamont-tests/unit-circle-pc3.diderot

Revision 4701 - (view) (download)

 1 : glk 4701 #version 2 2 : glk 4640 // unit-circle 3 : // 4 : // Demo of distributing particles on the unit circle 5 : // 6 : glk 4701 input vec2[] initPos = load_sequence("vec2.nrrd"); 7 : glk 4640 //input vec2[] initPos = {[0.15606558, -0.56803983]}; // BUG 8 : input int iterMax ("maximum number of iterations to run") = 1; 9 : input real rad ("radius of particle's energy support") = 0.2; 10 : input int pcp ("periodicity of population control") = 5; 11 : input real mvmtEps ("system convergence criterion on particle movement ") = 0.001; 12 : 13 : real hhInit = 10; // initial step size, can be big 14 : real newDist = 0.5*rad; // how far away to birth new particles 15 : real stepMax = 2*rad; // limit on distance to travel per iter 16 : 17 : int iter = 0; // which iteration we're on 18 : 19 : // inter-particle energy, and its derivative 20 : function real φ(real r) = (1 - r)^4; 21 : function real φ'(real r) = -4*(1 - r)^3; 22 : //function real φ(real r) = (1/r)*(1-r)^3; 23 : //function real φ'(real r) = 3 - 1/(r^2) - 2*r; 24 : function real enr(vec2 x) = φ(|x|/rad); 25 : function vec2 frc(vec2 x) = φ'(|x|/rad) * (1/rad) * x/|x|; // chain rule 26 : 27 : strand Particle (int parent_, int ii, vec2 pos0, real hh0) { 28 : int parent = parent_; 29 : vec2 pos = pos0/|pos0|; 30 : real hh = hh0; 31 : output vec2 outPos = pos; 32 : vec2 delta = [0,0]; 33 : int id = ii; 34 : int born = 0; 35 : real energyLast = ii; // HEY: can do convergence test based on variance of energies 36 : vec2 force = [0,0]; // or can test convergence based on sum of |force| 37 : bool first = true; 38 : update { 39 : if (first) { 40 : //print("(iter ", iter, ") hello from ", ii, "\n"); 41 : first = false; // HEY remove me 42 : } 43 : // compute energy and forces on us 44 : energyLast = 0; 45 : force = [0,0]; 46 : int ncount = 0; 47 : foreach (Particle p_j in sphere(rad)) { 48 : ncount += 1; 49 : vec2 r_ij = p_j.pos - pos; 50 : energyLast += enr(r_ij); 51 : force += frc(r_ij); 52 : } 53 : vec2 norm = normalize(pos); // surface normal (for unit circle) 54 : if (ncount == 0 && (pcp > 0 && 0 == iter % pcp)) { 55 : vec2 npos = pos + newDist*rad*[norm[1],-norm[0]]; 56 : new Particle(ii, born, npos, hh); 57 : born += 1; 58 : continue; 59 : } 60 : // else update position based on force from neighbors. 61 : // first, project force onto tangent surface 62 : force = (identity[2] - norm⊗norm)•force; 63 : delta = hh*force; 64 : if (|delta| > stepMax) { 65 : // decrease hh by factor by which step was too big 66 : hh *= stepMax/|delta|; 67 : // and find smaller step 68 : delta = hh*force; // == stepMax*normalize(force); 69 : } 70 : vec2 posLast = pos; 71 : // take step and re-apply constraint 72 : pos = normalize(pos + delta); 73 : real energy = 0; 74 : ncount = 0; 75 : foreach (Particle p_j in sphere(rad)) { 76 : energy += enr(p_j.pos - pos); 77 : ncount += 1; 78 : } 79 : if (energy - energyLast > 0.5*(pos - posLast)•(-force)) { 80 : hh *= 0.5; // backtrack 81 : pos = posLast; 82 : } else { 83 : hh *= 1.1; 84 : /* putting new particle creation here avoids possibility of 85 : starting multiple particles at same position while not 86 : moving due to backtracking */ 87 : if (ncount < 2 && (pcp > 0 && 0 == iter % pcp)) { 88 : vec2 npos = pos + newDist*rad*normalize(force); 89 : new Particle(ii, born, npos, hh); 90 : born += 1; 91 : } 92 : } 93 : outPos = pos; 94 : } 95 : } 96 : 97 : glk 4701 update { 98 : glk 4640 //print("(iter ", iter, ") hello from global\n"); 99 : real mvmt = max { |P.delta|/rad | P in Particle.all}; 100 : glk 4648 if (numActive() > 1 && mvmt < mvmtEps) { 101 : glk 4640 print("stabilizing with movment ", mvmt, " at iter ", iter, "\n"); 102 : stabilize; 103 : } 104 : iter += 1; 105 : if (iter >= iterMax) { 106 : print("stopping with movment ", mvmt, " at iter ", iter, "\n"); 107 : die; 108 : } 109 : } 110 : 111 : glk 4701 create_collection {Particle(-1, ii, initPos[ii], hhInit) 112 : | ii in 0 .. length(initPos)-1 }

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