User:Devon McCormick/DLA/SeekingRelease

From J Wiki
Jump to: navigation, search


We explore how the "release policy" shapes the resulting shape of a diffusion-limited aggregation (DLA).

Basic Code

We start with this basic code, all of which can be found here. First, load some sub-routines, define our own namespace dla, and create an initialization routine which defines a few globals:

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: 1-away in D-space
   PTS=: ,:y$0      NB. List of points: start with one at origin.
NB.EG  init 2       NB. Initialize 2-dimensional space.

Here's the main routine, which aggregates points to the global PTS:

NB.* aggStruc: randomly walk a point y times to aggregate it to cluster:
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.
NB.EG (5;PTS) aggStruc 1e2$1e2

The basic sub-routines:

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.
addPt=: 3 : 'PTS=: PTS,y'

NB.* walk: walk one random step.
walk=: 3 : 'y+RNI{~?#RNI'

Finally, an early version of the nub of the problem: the function that "releases" a new particle.

NB.* releaseIfOpen: release new particle x away to find open neighborhood.
releaseIfOpen=: 4 : 0
   while. 1 e. PTS e. RNI+"1 newpt=. ((]{~[:?#)y)+x*_1 1{~?D$2 do. end.

Make this our release policy (for now):

release=: releaseIfOpen

The main routine randomly moves the newly-released point until it either encounters the existing cluster and sticks, or finishes a pre-set number of moves without encountering the cluster, at which point it disappears.

We'll explain in more detail how the release policy works shortly but first let's see how to run the code. If we run the above lines of J, then move back into the base namespace, and make the dla namespace transparently available here, we can initialize our cluster.

   coclass 'base'
   coinsert 'dla'
   init 2
0 0
1 2

So, we start with one point at the origin; its shape is 1x2 because it's one point in two dimensions: we plan to concatenate new points to the bottom of this matrix (as can be seen by the definition of addPt above).

The arguments to the main routine are a parameter of 3 - to be explained shortly - the existing list of points and a right argument of how long our random walks should be before we give up. We'll arbitrarily try 20 walks of length 10 each:

   (3;PTS) aggStruc 20$10
4 10 10 10 10 10 10 1 10 10 10 10 3 10 10 7 10 9 10 10

The result is how long each walk went on. We see that most of them went the maximum - which means they didn't encounter the cluster (or encountered it only on the last step) - but a few ended early, hence became part of the cluster.

So, how many points do we now have and what are they?

7 2
0 _1 _2 _2  0 _1 _2
0 _1 _2  0 _2 _3  1

We have 7 points which we show transposed simply because it's a more efficient use of display space.

Let's try the same thing again:

   (3;PTS) aggStruc 20$10
10 10 1 10 6 10 10 10 3 10 10 10 2 3 10 1 10 4 4 10
|16 2|0 _1 _2 _2  0 _1 _2 _3  1 _4 _1 0 _2 1  2 _3|
|    |0 _1 _2  0 _2 _3  1  2 _2  1  2 3 _4 1 _1 _3|

This time we see that more points stuck. Again, we display the shape and the values as a vector of boxed items just to make more efficient use of our display.

How about a picture of these points? To do this, re-enter the namespace and define a new utility function bordFill, then return to the base namespace. This new function takes a left argument of how many empty cells we want to border our points and a right argument of our point co-ordinates.

   coclass 'dla'

   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

   coclass 'base'

We use this in conjunction with the general "viewmat" utility to display the points:

   viewmat 1 bordFill PTS

Which looks like this:


Explaining First Release Policy

Our first release function is very crude and its limitations will become apparent as the cluster grows but it's fine for a start. Here's most of the code:

while. 1 e. PTS e. RNI+"1 newpt=. ((]{~[:?#)y)+x*_1 1{~?D$2 do. end.

This is a while loop where all the work gets done in the conditional part. Here, we check if any of the immediate neighbors of the new, randomly chosen point newpt, are in the existing set of points. If so, we do nothing and try again until this is false. Once we exit the loop, we return newpt: a neighborless point near the cluster.

Let's examine this line in detail. The first part after while. is

1 e. PTS e. RNI+"1 newpt

which evaluates true if one or more points are in the set formed by RNI+"1 newpt. RNI is one of the three globals from our initialization: it looks like this (transposed for more efficient display):

_1 _1 _1  0 0  1 1 1
_1  0  1 _1 1 _1 0 1

Arranging these eight points around an origin,

   3 3$<"1]1 1 1 1 0 1 1 1 1#^:_1]RNI
|_1 _1|_1 0|_1 1|
|0 _1 |0 0 |0 1 |
|1 _1 |1 0 |1 1 |

we can see that they are the two-dimensional offsets around a given point. So, if we add each of these (one-dimensional) elements of RNI to a single point like newpt, we get the co-ordinates of all of the neighbors of newpt, e.g. for (5,5) shown here arranged two-dimensionally with the point itself inserted into the middle:

   (13 : '(<y)(<1 1)}3 3$<"1]1 1 1 1 0 1 1 1 1#^:_1]RNI+"1 y') 5 5
|4 4|4 5|4 6|
|5 4|5 5|5 6|
|6 4|6 5|6 6|

The first part of the assigment of newpt:  newpt=. ((]{~[:?#)y)+x*_1 1{~?D$2 uses a tacit function  (]{~[:?#) which generates a random number ( ? ) based on the length ( # ) of its right argument and uses this to select ( {~ ) an element (row). This is applied to the right argument y which is the list of points in the cluster, so we randomly choose one of our existing points.

This random point is added to the latter part of this expression -  x*_1 1{~?D$2 in which we apply an offset determined by the left argument x multiplied by a random combination of ones and negative ones set by _1 1{~?D$2. This multiplication essentially moves x points in a random direction from the original point. The D is in there to support generalization to more than two dimensions.

So, to return to our long-awaited explanation for the magic number 3 in our invocation of the main function ( (3;PTS) aggStruc 20$10 ), this is the random distance we travel to find an empty neighborhood in which to release our new point. So, if we pick a larger number for this argument, we'll pick release points, on average, further away from the cluster. This means the random walk will have to go on longer to have good chance of hitting the cluster.

So, if we choose a much bigger number than 3, we might see something like this:

   (20;PTS) aggStruc 20$10
10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10
16 2

We added no new points because a random walk of 10 steps is not long enough to travel (more or less) 20 points to where the cluster is. If we are willing to walk a lot more, say on the order of 20^2 steps, it's a different story:

   (20;PTS) aggStruc 20$400
400 288 400 400 400 400 309 400 400 400 400 400 400 247 400 400 400 400 400 400
19 2

However, as you might guess, it takes more time to walk further, i.e. starting with the original 16 points and timing 20 walks of 10 versus 20 walks of 400:

   6!:2 '(3;PTS) aggStruc 20$10'
30 2
   $PTS_dla_=: 16{.PTS_dla_
16 2
   4!:55 <'PTS'   NB. Ensure local copy doesn't shadow one in "dla" namespace.
   6!:2 '(20;PTS) aggStruc 20$400'
18 2

So, we see that the longer walks take many times as long and add fewer points. In fact, this set of parameters will often add no points.

Tricky Parameter

It would seem that we would want to make the initial offset parameter as small as possible for more efficient code. However, how small can this value be? If it's zero, we'll never find an empty neighborhood because we'll always choose a point on the cluster. If it's one, it fails for the same reason. How about two?

   6!:2 '(2;PTS) aggStruc 20$10'
30 2

This seems to work fine. How many points per second are we adding?


That's over 1800 points/second - what if we continue doing this?

   6!:2 '(2;PTS) aggStruc 20$10'
44 2
   6!:2 '(2;PTS) aggStruc 20$10'
57 2

It's a variable amount but still thousands per second. To reduce our typing, let's combine these separate lines into one. We'll have to save the starting number of points to calculate the number per second.

   (sv-~#PTS)%tm=. 6!:2 '(2;PTS) aggStruc 20$10' [ sv=. #PTS
   (sv-~#PTS)%tm=. 6!:2 '(2;PTS) aggStruc 20$10' [ sv=. #PTS
   (sv-~#PTS)%tm=. 6!:2 '(2;PTS) aggStruc 20$10' [ sv=. #PTS
   (sv-~#PTS)%tm=. 6!:2 '(2;PTS) aggStruc 20$10' [ sv=. #PTS
   (sv-~#PTS)%tm=. 6!:2 '(2;PTS) aggStruc 20$10' [ sv=. #PTS
123 2

This seems to be going well - what's it look like?

   viewmat 1 bordFill PTS
   savemat_jviewmat_ 'showUse01.bmp'


Let's continue but with 30 releases each time instead of 20 but, first, we'll save even more typing by defining a calling function to encapsulate what we're doing repetitively:

   do1_dla_=: 3 : '(sv-~#PTS)%tm=. 6!:2 ''(2;PTS) aggStruc y'' [ sv=. #PTS'

Since this is only one line long, we'll create it directly in the namespace instead of entering and exiting.

Now use this to run 3 batches of 30 walks:

1400.4003 931.17767 1604.4576
190 2

Now try 40 at a time:

1476.1391 822.55663 789.23452

The number per second seems to be declining. Let's try 50 at a time:

754.1739 855.81707 584.68598 669.20533

The rate of increase has definitely dropped. What if we try 100 at a time?

332.41549 394.68151 306.51331
643 2

Is it slowing because we're increasing the number of points we try at once? This seems counter-intuitive but let's test it.

146.76733 229.80636 277.27224 333.09454 379.94374 473.80991 247.69881 318.09695 179.85791
780 2
   viewmat 1 bordFill PTS
   savemat_jviewmat_ 'showUse04.bmp'


It's definitely slowing. A little reflection gives us a clue: since we're moving a small distance away from points in the cluster to find an empty neighborhood, we must make more attempts each time as the cluster grows since only points from the perimeter are likely to give us good hits.

Concentrating on the Perimeter

Fortunately, it's easy to see how to address this: instead of giving all the points as an argument to aggStruc, why don't we give only the perimeter points? This will reduce the number of false hits. How do we know which points are on the perimeter? One simple way is to look at only the most recently added points since we are, by the nature of our aggregation, adding to the perimeter.

Guessing that the most recent 10% of points is a reasonable bunch to try, let's extract those:

   $PERIMPTS_dla_=: PTS{.~->.10%~#PTS
78 2

Looking at how these compare to the whole set:

   viewmat 2 (<"1]1+PERIMPTS-"1 <./PTS)}1 bordFill PTS

The manipulation we do to PERIMPTS is necessary to map the origin-centered co-ordinates to positive only co-ordinates. This is a replication of what bordFill does internally as well.


This looks good as the perimeter is indicated by the pink points around the edge. Let's try using this as our argument to aggStruc and see how it affects our timings.

1005.5137 995.09202 564.62669
|attention interrupt: releaseIfOpen
|   point=.(>0{x)    release>1{x[ctr=._1
831 2

The points/sec initially goes back up but quickly declines, then the code freezes up and has to be interrupted. A look at the points and their perimeter makes the problem clear:

   viewmat 2 (<"1]1+PERIMPTS-"1 <./PTS)}1 bordFill PTS


We see that we quickly lose our original perimeter. Let's see how costly it is to re-calculate the perimeter with each call to the main routine:

   do1_dla_=: 3 : '(sv-~#PTS)%tm=. 6!:2 ''(2;PTS{.~->.10%~#PTS) aggStruc y'' [ sv=. #PTS'
   NB.  Inserted the perimeter calculation   ^^^^^ here ^^^^^^
1009.626 1297.3126 1809.8097 1466.0253 1293.6868

Increasing the number of attempts per call:

1578.9208 857.49015 1561.9932 1310.9005

Increasing the number of attempts again:

1033.4018 1404.6138 1208.4397 1019.4336

Increasing again:

997.21846 893.30982 1098.0115 820.66574 1135.0023 2141.7607 807.80606

This seems to be working fine though there does seem to be a gradual decrease in efficiency.

1358 2
892.7828 859.53069 827.92105 809.63688 852.42751

This looks pretty good but there is an aesthetic problem:


This DLA looks a bit "chunky" compared to other examples we've seen. Our policy of releasing very near the existing cluster has a tendency to fill in the gaps between branches of the cluster. Also, since we pick release points randomly in relation to the perimeter, sometimes we release within the existing cluster as indicated by the interior pink points on the lower part of this picture. What we'd really like to do is to pick release points only outside the perimeter. One way to do this is shown here as well in the following section.

A simple measure of how densely packed the above-pictured DLA is this:

   (13 : '(#y)%~*/$0 bordFill y') PTS

This shows us that 1 in every 2.5 points in the bounding rectangle is part of the cluster. If we add more points, getting up to almost 62,000 using these parameters, this sparsity ratio goes up to 4.07.

Same Method, Different Parameters

Before we show a new and different release policy, let's build a cluster with our existing policy using slightly different parameters. First, initialize a new set of points:

   init 2
   do1 10$10
6 2

Then we'll vary the minimal offset distance by parameterizing our do1 function:

do1_dla_=: 4 : '(sv-~#PTS)%tm=. 6!:2 ''(x;PTS) aggStruc y'' [ sv=. #PTS'

and use it with a larger offset distance than the "2" we used for the previous cluster:

   5 do1&>10$<10$30
72.617714 607.66382 155.76784 157.00792 281.65434 334.79505 543.02723 342.61396 477.80781 496.56591

Since we'll be doing this a lot, let's summarize the "points-per-second" result we get from do1. The verb usus shows the usual statistics in which we're interested: minimum, maximum, mean, and standard deviation of the points-per-second results generated by the do1 verb.

   usus 10 do1&>100$<10$100
29.740074 1414.9642 320.73076 187.40271
   usus 10 do1&>100$<10$100
50.720809 720.61703 243.40365 126.31516
   usus 10 do1&>100$<10$100
37.930561 2157.7823 222.73005 227.97863
   usus 10 do1&>100$<10$100
25.294278 441.61624 160.07025 77.476433

From these statistics, we see the gradual decline we expect as the cluster grows, giving us more fruitless picks from the interior proportional to fruitful ones from the perimeter. We continue until we have 11,652 points.

   usus 10 do1&>100$<20$100
41.947957 316.08637 115.16387 46.413008
   usus 10 do1&>100$<40$100
36.47076 130.91957 71.789116 17.662557
   usus 10 do1&>100$<80$100
24.403297 74.751489 43.270225 10.192174
11652 2
   $PPTS=: PTS{.~->.10%~#PTS
1166 2
   viewmat 2 (<"1]1+PPTS-"1 <./PTS)}1 bordFill PTS

Taking a look at the perimeter points we'd get at this point, Show2 10 04.png we see the stray interior ones that tend to fill up the cluster. This one has a sparsity slightly greater than that (4.07) of our previous cluster generated using this same method but with a smaller offset.

   (13 : '(#y)%~*/$0 bordFill y') show2pts03
5.0468589    NB. Bounding rectangle points/cluster points

New Release Policy

(The following method is also elaborated here in the context of an earlier version of the DLA code which works on a fixed-size, boolean matrix rather than an explicit set of points.)

First, we define a verb to generate all the immediate neighbors of a point set not in that set.

neigh2=: 13 : 'x-.~~.,/y+"1/RNI'   NB. Empty neighbors of these

We use this to get the nearest neighbors:

   $NP=: PTS neigh2 PTS
26386 2

Then we expand the neighborhood by finding the next closest set of neighbors to these, then their neighbors, and so on five times:

   6!:2 'NP=: PTS neigh2^:5]NP'

This gives us an envelope six deep around our original points. We find the edge of this envelope by isolating those outermost neighbors (those with neighbors not in the set of neighbors):

   $edges=. ~.PTS-.~NP -.~ ,/NP+"1/RNI
1029 2
   viewmat 2 (<"1]1+edges-"1 <./PTS,edges)}1 bordFill PTS,edges

Here's what the edge of the envelope looks like in relation to our existing cluster:

Show2 edges00.png

This looks like a good perimeter set on which to release points but it's a bit thin. Let's thicken it by finding its neighbors two levels deep:

   $edges=. ~.PTS neigh2^:2]edges
5138 2

Giving us this fatter release area:

Show2 edges01.png

Adjust our top-level calling routine to use this perimeter points global:

   do1_dla_=: 4 : '(sv-~#PTS)%tm=. 6!:2 ''(x;PPTS) aggStruc y'' [ sv=. #PTS'
   PPTS_dla_=: edges
   1 do1&>10$<10$25
12.281323 21.910126 23.042773 79.635403 23.813149 22.408149 34.180425 48.340221 21.374501 50.795483

Our initial points/second are a bit low but seem to be improving. Let's try some more:

12590 2
   viewmat 2 (<"1]1+PPTS-"1 <./PTS,PPTS)}1 bordFill PTS,PPTS
   1 do1&>20$<20$25
24.342477 61.743025 41.15663 27.630435 34.954674 34.740216 60.580126 52.822625 57.748258 22.060313 22.321549 52.2986 21.649241 40.913895 67.019374 34.597901 40.404442 52.186518 39.882113 42.753015

One problem we'll have if we try to generate too many points using one perimeter is that the cluster will eventually "invade" the perimeter. If we don't do something to adjust it, this will lead to unsightly clumps of aggregation within the perimeter band. First, we'll make a new verb to create a perimeter:

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
   $PPTS_dla_=: calcPerim 5 2
5314 2

Running this gives satisfactory statistics (min, max, mean, SD of points/second):

   usus 1 do1&>50$<20$25
5.0021381 49.32456 27.373553 10.544996
12969 2

We can adjust our main routine to incorporate the perimeter calculation directly:

   do1_dla_=: 4 : '(sv-~#PTS)%tm=. 6!:2 ''(x;calcPerim 3 1) aggStruc y'' [ sv=. #PTS'
   usus 1 do1&>50$<20$25
12.040371 61.52522 30.262127 10.363449

After trying out this version, it seems that the perimeter calculation, being somewhat expensive, should be done only every so often with the result assigned to a global like this:

   do1_dla_=: 4 : 0
   if. 0=?10 do. PPTS=: calcPerim 5 2 end.
   sv=. #PTS
   tm=. 6!:2 '(x;PPTS) aggStruc y'
   usus 1 do1&>50$<20$25
0 62.962738 22.075551 13.289167

After running this enough times to add a significant number of points, we look at the cluster with its perimeter: Show2 edges02.png

This is interesting because we can plainly see the effect of the new release policy on the appearance of the cluster. The inner, denser portion was formed using our flawed "releaseIfOpen" method. The outer, sparser part of the cluster was formed using our new perimeter-based release policy.

If we generate a cluster solely based on this newer release policy, we get something like the following with a sparsity measure of 6.7, which makes this less dense than the clusters generated by our earlier policies.