User:Devon McCormick/DLA/DiffLimAgg
Jump to navigation
Jump to search
NB.* DiffLimAgg.ijs: diffusion-limited aggregation of random particles. NB. Adapted from http://www.jsoftware.com/jwiki/Studio/Gallery/DLA load 'coutil' cocurrent 'dla' load '~Code/math.ijs mystats viewmat' init=: 3 : 0 D=: y NB. How many dimensions? RNI=: D dimnbr 1 NB. Relative neighbor indexes PTS=: ,:y$0 NB. List of points ) NB.EG init 2 NB. Initialize 2-dimensional space. NB.* aggStruc: make max rand walk y tries to aggregate point to structure, NB. (>0{x) away from random point in >1{x. aggStruc=: 4 : 0"(_ 0) point=. (>0{x) release >1{x [ ctr=. _1 while. y>ctr=. >:ctr do. point=. walk point if. PTS check_collision point do. y=. ctr [ addPt point end. end. ctr NB.EG (5;PTS) aggStruc 1e2$1e2 ) NB.* check_collision: 1 iff x in neighborhood of points y. check_collision=: 4 : '1 e. x e. RNI+"1 y' NB.* dimnbr: x-dimensional y-nearest neighbor offsets. dimnbr=: 13 : '(x$0)-.~,/>{x$<y-~i.>:+:y' NB. RNI=: (D$0)-.~(_1 0 1){~(D$3)#:i.3^D NB. Relative neighbor offsets NB.* addPt: add point to cluster, track indexes of points added. addPt=: 3 : 'PTS=: PTS,y' NB.* walk: walk one random step. walk=: 3 : 'y+RNI{~?#RNI' NB.* releaseIfOpen: release point x away from random point in y if NB. open (empty) neighborhood there. releaseIfOpen=: 4 : 0 while. 1 e. PTS e. RNI+"1 newpt=. ((]{~[:?#)y)+x*_1 1{~?D$2 do. end. newpt NB.EG dis releaseIfOpen PERIMPTS ) NB.* release: release new particle x away to find open neighborhood. release=: releaseIfOpen NB.* bordFill: fill 0-mat w/1 according to y, having border of x cells. bordFill=: 4 : '(1)(<"1 y-"1 x-~<./y)}0$~2$(>:+:x)+(>./y)-<./y' NB.EG viewmat 1 bordFill PTS neigh2=: 13 : 'x-.~~.,/y+"1/RNI' NB. Empty neighbors of these points. NB.EG PTS neigh2 PTS calcPerim=: 3 : 0 NP=: PTS neigh2^:(0{y)]PTS edges=. ~.PTS-.~NP -.~ ,/NP+"1/RNI edges=. ~.PTS neigh2^:(1{y)]edges NB.EG PPTS_dla_=: calcPerim 5 2 ) do1=: 3 : '(sv-~#PTS)%tm=. 6!:2 ''(2;PTS) aggStruc y'' [ sv=. #PTS' do1_dla_=: 4 : 0 if. 0=?10 do. PPTS=: calcPerim >1{x end. ((>0{x);PPTS) aggStruc y ) NB.* countNeighbors: count neighbors x in each direction from point y. countNeighbors=: 4 : '+/PTS e. y+"1 D dimnbr x'"1 NB.* pickGTMeanDis: pick points > mean dis from center of cluster. pickGTMeanDis=: [:p2c"1 [:(]#~[:(]>mean)0{"1])c2p"1 NB.* estPerimeter: estimate peripheral points by who has fewer neighbors. estPerimeter=: 3 : 'pickGTMeanDis y#~2><.-:1 countNeighbors"1 y' NB.* pushOutPerimeter: push out perimeter x cells from center of points y. pushOutPerimeter=: 4 : 0 y#~ad>x+<./ad=. 0{"1 c2p"1 y-"1 mean y ) alt0release=: 4 : 0 if. 0=?1000 do. PERIMPTS=: estPerimeter PTS end. x releaseIfOpen PERIMPTS ) alt1release=: 4 : 0 if. 0=?10000 do. RELEASEDIS=: >:RELEASEDIS end. <.y+RELEASEDIS r. 2p1*?0 ) NB.** Experiment with releasing near the perimeter of the cluster. egShowPerimeter=: 0 : 0 pp=. c2p"1 PERIM#PTS-"1]500 500 mpc=: (mean 0{"1 pp),.2p1*-:>:i:1j999 NB. Mean perimeter circle: 1000 points mpc=: p2c"1 mpc ofst=. (mean mpc)-~mean 340-~PTS opc=: (25+mean 0{"1 pp),.2p1*-:>:i:1j999 NB. Outer perimeter circle: 1000 points NB. PERIMPTS=: <.0.5+opc=: p2c"1 opc knownspace=: 1 bordFill PTS oofst=. (mean opc)-~mean 340-~PTS mpcmat=: (10)(<"1]<.0.5+ofst+"1 mpc)}(2><.-:PTN1) (<"1]340-~PTS)}knownspace (<./,>./)<.0.5+oofst+"1 opc viewmat (20)(<"1]<.0.5+oofst+"1 opc)}mpcmat ) NB.** Apply perimeter-release with periodic perimeter adjustment. reduceToIn=: 13 : 'x#~x e. y' adjPts=: 13 : '<"1 y-x-~<./,y' growOut=: 3 : 0 np=. 10 [ po=. 10 while. 0<:y=. <:y do. tm=. 6!:2 'nn=: (np;PERIMPTS) aggStruc 1e3$10+*:np' smoutput (qts ''),tm,(#PTS),(nn+/ . =>./nn),usus nn PERIMPTS=: estPerimeter PTS PERIMPTS=: PTS reduceToIn ~.,/PERIMPTS+"1/D dimnbr 2 PERIMPTS=: PTS reduceToIn ~.,/ PERIMPTS+"1/D dimnbr 2 PERIMPTS=: po pushOutPerimeter PERIMPTS if. 0=?5 do. np=. >:np [ po=. >:po end. end. viewmat 2 (<"1]1+PERIMPTS-"1 <./PTS)}1 bordFill PTS ) NB.* bumpAggParms: aggregate with changing parameters. bumpAggParms=: 3 : 0 'bgn end inc np npp'=. y [ stats=. 0 7$0 while. end>:bgn do. tm=. 6!:2 'nn=. (bgn;npp{.PTS) aggStruc np$*:bgn' stats=. stats,tm,(#PTS),(usus nn),nn+/ . = >./nn bgn=. bgn+inc end. stats NB.EG bumpAggParms 4 34 2 100 _30 NB. 4->34 by 2; 100 pts at a time; last 30. )