User:Devon McCormick/DLA03

Here's some work being done by Ric Sherlock and Devon McCormick to come up with a good Rosetta Code entry for the "Brownian Tree" problem - another name for Diffusion-Limited Aggregation.

Here's an initial attempt with a few helper functions:

NB.* brownianTree3.ijs: make diffusion-limited aggregation.

NB.* brownianTree3: aggregate points in cluster by Brownian motion.
brownianTree3=: 3 : 0
   'world npts'=. y           NB. Boolean mat world, # pts to insert.
   if. 0=#world do. sz=. 100  NB. Default to 100x100 world
       peri=. 2#3-~sz         NB. Track periphery of cluster
       world=. 1 (<<.-:peri)} (2#sz)$0 end.  NB. Seed center
   peri=. 2#3-~sz=. #world
   for_pt. i. npts do.        NB. Random pt on periphery,
       xy=. >:?peri           NB.  not on populated cell or edge.
       while. -.+./world{~<"1 nn=. getNbrs xy do. NB. No populated neighbor?
           xy=. nn {~ ?#nn                        NB. Move to random neighbor
           if. +./ xy (1&>@[ , >) peri do.        NB.  near periphery.
               xy=. >:?peri
           end.
       end.
       world=. (1) (<xy)}world
  end.
NB.EG brownianTree3 '';100
)

NB.* getNbrs: get neighbors of a cell.
getNbrs=: (0 0 -.~ <: 3 3#: i.9) +"1 ]
NB.* pad0s: pad rows and cols of numeric mat with x zeroes in each direction.
pad0s=: 4 : 'y (<x+&.>i.&.>$y)}0$~(2*x)+$y'
NB.* minmaxrow: indexes of min and max occupied (=1) cell of boolean mat.
minmaxrow=: [:(([:<./0{"1]) , [:>./1{"1]) [:>(<_ __)-.~ [:(<./,>./) &.>[:I.&.><"1
NB.* minmaxrc: min and max occupied (=1) row and column of boolean mat.
minmaxrc=: minmaxrow,minmaxrow@:|:
NB.* trim0s: trim all-zero rows and columns from boolean mat.
trim0s=: [: (] #"1~ 0 +./ .~:])] #~ 0 +./ .~:"1 ]
NB.* density: portion of occupied cells in boolean mat.
density=: ([: +/,) % [:#,

One problem is that this is very inefficient code for a large, empty matrix. That's why the arguments are an existing matrix and a number of points: this allows us to start with a small world so that short random walks will encounter the cluster, then increase the size of the matrix by padding the borders as the cluster grows. This means a typical use of this code will be something like the following.

First, start with the default 100x100 for 100 points:

   sz=. +/,w3 [ tm=. 6!:2 'w3=. brownianTree3 '''';100'  NB. Start tracking number of cells and run-times.
   sz,:tm
       99
2.6285197
   viewmat w3

This gives something like this:

Next, aggregate more cells to the existing cluster:

   sz=. sz,+/,w3 [ tm=. tm,6!:2 'w3=. brownianTree3 w3;500'
   sz,:tm
       99       568
2.4815702 3.0286292

Notice how much more efficient is the second round since we have a larger target, hence shorter random walks. This gives us something like this:

Are we close to the borders of our world yet?

   minmaxrc w3
21 82
19 81

No, this looks like a sufficient margin to continue adding more points. Let's incorporate this check as part of our routine:

   minmaxrc w3 [ sz=. sz,+/,w3 [ tm=. tm,6!:2 'w3=. brownianTree3 w3;1000'
 8 88
10 90
   sz,:tm
       99       568      1392
2.4815702 3.0286292 1.0529608

Now it looks like this:

Since we're approaching the borders of our world, let's make it larger the next time we add more points:

   minmaxrc w3 [ sz=. sz,+/,w3 [ tm=. tm,6!:2 'w3=. brownianTree3 (50 pad0s w3);2000'
   sz,:tm
33 165
29 168
      sz,:tm
       99       568      1392      3187
2.4815702 3.0286292 1.0529608 21.477968

We notice a drop in efficiency as we increased the empty area by padding the borders to give us a 200x200 world.

We have enough free border space to add another bunch of points without further padding for now:

   minmaxrc w3 [ sz=. sz,+/,w3 [ tm=. tm,6!:2 'w3=. brownianTree3 w3;2000'
21 178
20 177
      sz,:tm
       99       568      1392      3187      4882
2.4815702 3.0286292 1.0529608 21.477968 6.5040638

The aggregate now looks like this:

If we continue growing this cluster in this fashion, we'll see some aesthetic flaws in our eventual result. Here's a sample of a 500x500 world:

This looks too "filled-in" in the center. Also, it's assuming a suspiciously square outline - possibly because random, independent (X,Y) selection fails to pick points uniformly on a circle? Or maybe it has do with some other aspect of how we choose the random "release points" for the new particles.