User:Devon McCormick/DLA/DiffLimAgg02

Here's the penultimate version of our re-vamped DLA code as explained in some detail here.

NB.* DiffLimAgg02.ijs: diffusion-limited aggregation of random particles.
NB. Adapted from http://www.jsoftware.com/jwiki/Studio/Gallery/DLA;
NB. explained in the following sections as shown in "DLA_TOC" (following).

DLA_TOC=: 0 : 0
(All the links following are rooted in the J wiki at http://www.jsoftware.com/jwiki/)
The project starts here with an introduction to diffusion-limited aggregation and an implementation already on the J wiki:
    NYCJUG/2009-10-13#DLA.3AInitialJImplementation

It continues here with a superficial re-vamp that aims to preserve how the original code works but make it superficially more J-like:
    DevonMcCormick/DLA00

We look at how this revised version performs by devising some performance metrics:
    DevonMcCormick/DLA0Runs

Then we attempt a deeper re-vamp to make the code more J-like here:
    DevonMcCormick/DLA01
This section also makes use of the J profiler.

We address a basic flaw in the original implementation here:
    DevonMcCormick/DLA00/BuildingACloseNeighborhood

We explore the most crucial determinant of what shapes a DLA cluster by developing the "release policy" here:
    DevonMcCormick/DLA/SeekingRelease
This also summarizes the basic, re-vamped code that works with an explicit, multi-dimensional list of points but releases only a single particle at a time.
)

load 'coutil'
cocurrent 'dla'
NB. load 'mystats'
load '~Code/math.ijs 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
)