User:Devon McCormick/DLA/DiffLimAgg

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.
)