User:Devon McCormick/DLA/MultiRelease

Here we wrap up our DLA essays with some performance enhancements.

Multiple Releases

One last improvement to consider in the current implementation is driven by the question: why do we only work with a single particle at a time? This limitation we inherit from the original implementation - it's typical of inherently scalar languages.

This scalar orientation often hurts the performance of a naive implementation of an algorithm in J by introducing an (often) unnecessary bottleneck. In the case of DLA, there seems to be no good reason to retain this constraint. In fact, the point to which we've progressed in the evolution of our implementation lends itself well to a multi-release strategy.

Consider the perimeter we build around an existing cluster.

Let's call this RC for "release candidates":

   $RC_dla_=: calcPerim 6 3
9950 2

The arguments to calcPerim are two numbers which may be thought of as the distance and width of the perimeter. However, the details are slightly more intricate as seen by this small variation on the above arguments:

   $altrc=: calcPerim 5 3
11696 2
   viewmat 2 (<"1]1+altrc-"1 <./PTS,altrc)}1 bordFill PTS,altrc

Note the small islands of release candidates (in pink) within the cluster: these reflect sufficiently empty interior areas in which we have empty space at least five cells away from points in the cluster.

Building the Code

There's no reason to choose only a single point from this set when we could choose an arbitrary number of points from it. Let's see how this works. First, we'll initialize some variables we'll need: the number of points we'll start with, a counter to track the length of the walks, a vector of results, a maximum number of walks to try before giving up, and the random set of points from the list of release candidates.

   numpts=. 100
   ctr=. _1 [ rr=. i.0
   maxwalk=. 100
   $pts=. numpts (]{~ ([<.[:#]) ?[:#]) RC
100 2

At the top of our loop, we'll check if we've exceeded the maximum walk length or run out of points.

   (maxwalk>ctr=. >:ctr) *. 0<#pts
1

Now, we walk all the points and check for collisions. However, we encounter some problems because we have to modify some of our routines to handle multiple points at a time. Fortunately, the changes required are trivial: we simply specify the ranks to which our existing functions should apply.

   wcc=. PTS check_collision pts=. walk pts
|length error: walk
|   y    +RNI{~?#RNI
   walk
3 : 'y+RNI{~?#RNI'
   walk_dla_=: 3 : 'y+RNI{~?#RNI'"1

We have to modify walk to work on vectors to accommodate our non-scalar argument, then check again.

   wcc=. PTS check_collision pts
|length error: check_collision
|   1 e.x e.RNI    +"1 y
   check_collision
4 : '1 e. x e. RNI+"1 y'
   check_collision=: 4 : '1 e. x e. RNI+"1 y'"1
   wcc=. PTS check_collision pts
|length error: check_collision
|   wcc=.PTS     check_collision pts
   check_collision_dla_=: 4 : '1 e. x e. RNI+"1 y'"(_ 1)

Similarly, we need to modify our collision check to specify that we apply the entire left argument to vectors on the right, then we'll try again.

   wcc=. PTS check_collision pts
   1 e. wcc
0

We have a vector wcc where a "1" indicates we have a collision of the corresponding point with the cluster. This time, and the next few, we don't have any collisions, so we keep walking the points and counting the number of steps we've taken so far.

   1 e. wcc=. PTS check_collision pts=. walk pts [ ctr=. >:ctr
0
   1 e. wcc=. PTS check_collision pts=. walk pts [ ctr=. >:ctr
0
   1 e. wcc=. PTS check_collision pts=. walk pts [ ctr=. >:ctr
1
   1+/ . = wcc
1

Now we have one collision, so we need to add the colliding point to the cluster and remove it from our original list of released points.

   #PTS
12075
   #addPt wcc#pts
12076
   pts=. pts#~-.wcc [ rr=. rr,ctr$~+/wcc
   ctr
3
   rr
3
   $pts
99 2

Testing the Code

Now that we've figured out the body of the loop, we can put it all together in a function.

growMany_dla_=: 3 : 0
   'maxwalk numpts pp'=. y
   if. -.nameExists 'RC' do. RC=: calcPerim pp
   else. if. 0=?10 do. RC=: calcPerim pp end. end.
   pts=. numpts (]{~ ([<.[:#]) ?[:#]) RC
   ctr=. _1 [ rr=. i.0
   while. maxwalk>ctr=. >:ctr *. (0<#pts)
   do. wcc=. PTS check_collision pts=. walk pts
       if. 1 e. wcc do. addPt wcc#pts
           pts=. pts#~-.wcc [ rr=. rr,ctr$~+/wcc
       end.
   end.
   rr,maxwalk$~#pts
)

We've also added some setup code to ensure that we have a global RC of release candidates, creating it if it's missing. Also, since we'll need to refresh this list as the cluster grows toward overlapping this perimeter, we randomly refresh RC about every ten iterations.