NYCJUG/2009-10-13
convex hull, DLA, geometry, literate programming, evangelizing J
Location::BEST, Hoboken, NJ
Meeting Summary
We looked at some different problems with possible solutions involving finding the "convex hull" of a set of points in two dimensions. We also compared some different styles of programming and program exposition by considering literate programming in contrast to the J-oriented divide between tacit and explicit styles.
Agenda for NYCJUG of 20091013
1. Beginner's regatta: improving code - looking at "diffusion-limited aggregation" - see "DLAIntroInclJavaEG.odt" and "ExamplesJDLAOutput.odt". 2. Show-and-tell: packing points problem - see "Point-Packing with Convex Hull.odt" and "GeometryFns.odt". 3. Advanced topics: Literate programming - example of convex Hulls - see "LiterateProgrammingConvexHullEGs.odt". 4. Learning, teaching and promoting J: literate versus tacit - see "EvangelismAndOtherAgendasForRosettaCode.odt".
Proceedings
We were well-supplied with food and drink by our generous hosts at Bayesian Enhanced Strategic Trading but we didn't take full advantage of their conferencing and display facilities at this meeting. We need to arrange topics to make better use of these aids.
Beginner's regatta
We began by considering some code on the J wiki which is an implementation of "diffusion-limited aggregation" - a fractal growth model of a wide range of physical phenomena including crystal growth and diffusion of liquids - such as water and oil - in confined areas - such as fractured rock or soil.
There's a good introduction to this and a Java applet demonstrating how it looks here. Also, here's a technical article, titled "Diffusion-Limited Aggregation: A Model for Pattern Formation" (Halsey), that's a readable and comprehensive introduction to the general subject.
DLA: Initial J Implementation
There's an implementation of DLA on the J wiki that was my inspiration for taking up this topic because it's an interesting problem with some working code but the code has room for improvement, both to make it work better and to make it more "J-like". It's a challenge to make this more amenable to a J solution because the problem, as stated, lends itself to a very scalar, linear solution which does not take much advantage of J's strengths. The code on the wiki reflects a straight-forward implementation of the problem.
However, even though the code could be improved, as it now stands, it illustrates one of J's strengths: clarity from tersness. As we see below, the high-level outline of the algorithm is fairly clear from the main routine "grow".
grow =: 3 : 0
c =. 0
point =: size release center
while. y > c do.
point =: walk point
if. point outOfBounds center do.
point =: size release center
else.
if. check_collision point do.
add_to_cluster point
size =: center dist point
break.
end.
end.
c =. >: c
end.
)
We see the general outline in which a particle is "released" onto a grid and allowed to wander randomly until it collides with the existing structure, at which point it sticks to the structure to make it slightly larger. This process is then repeated.
The general exercise of improving this code, at least the first set of revisions, is shown on a separate page.
However, though this code runs, I initially encountered two problems, one of which I addressed as a bug-fix.
At first, the results look reasonable.
6!:2 'grow_dla_"0]1e3$1e3' 47.170839 viewmat space_dla_
But, as we run for a longer time, we get an error.
6!:2 'grow_dla_"0 (10000#1000)'
|index error: check_collision
| +/+/(<drange y) {space
Turn on debug to see what the problem is.
13!:0]1
6!:2 'grow_dla_"0 (10000#1000)'
|index error: check_collision
| +/+/(<drange y) {space
|check_collision[0]
y
265j499
$space
500 500
drange y
+-----------+-----------+
|264 265 266|498 499 500|
+-----------+-----------+
$space
500 500
We're indexing out of bounds - we've hit the edge of our lattice. Put together a check for this and put it in the main routine.
check_oob_dla_=: 3 : '(0 0+./ . >+.y)+.($space)+./ . <: +.y' check_oob_dla_"0]0j0 _1j0 499j499 499j500 500j499 0 1 0 1 1
Subsequent errors led me to modify this to warn when the point is one cell away from the edge as the next check of the neighborhood would give an error.
NB.* check_oob: Check if point close to out-of-bounds (within 1 of any edge). check_oob=: 13 : '(2+./ . >:y)+. y+./ . >:$space [ y=. >:+.y'
Once we fix these boundary problems, by changing this line in the main routine:
if. (point outOfBounds center)+.check_oob point do.
we encounter a more serious problem as we experiment with the efficiency of different parameters.
+/,space_dla_[smoutput 6!:2 'grow_dla_"0]1e4$200' 99.611579 1895 +/,space_dla_[smoutput 6!:2 'grow_dla_"0]1e4$400' 196.99597 2138 +/,space_dla_[smoutput 6!:2 'grow_dla_"0]1e4$800' 386.72237 2495
Taking these above timings and counts of the number of points in the cluster, we can calculate how many points-per-second we're adding. This is a basic measure of efficiency.
NB. Points/sec (2-~/\1792 1895 2138 2495)%99.6 197 386.7 1.0341365 1.2335025 0.92319628 +/,space_dla_[smoutput 6!:2 'grow_dla_"0]1e4$400' 183.12599 2495
However, we soon notice that we're taking time but failing to add points, even as we change the parameters.
NB. Uh-oh – most efficient before doesn't work now +/,space_dla_[smoutput 6!:2 'grow_dla_"0]1e4$800' 365.47759 2495 +/,space_dla_[smoutput 6!:2 'grow_dla_"0]1e4$1200' 550.43562 2495
A look shows at the cluster shows why: we've painted ourselves into a corner.
Somehow this implementation favors an ever-lengthening cluster. Some examination of the code for "release" gives a reason why this might be.
NB. Release a new particle distance x + 5 from center y release=: 4 : '<.y + (x + 5) r. 2p1 * ? 0'
Looking at where this is used in the core of the main loop:
point=: walk point
if. (point outOfBounds center)+.check_oob point do.
point=: size release center
else. if. check_collision point do.
add_to_cluster point
size=: center dist point
break.
end.
end.
we see that the release point moves out from the center depending on the "size" global and this is set at the point of the previous addition to the cluster. If the cluster elongates in any direction, the end of it is most likely to contact any newly-released point.
This illustrates the importance of the algorithm for releasing a new point, which will be the subject of some consideration in future enhancements to this code. Other ways of releasing a particle achieve their own characteristic shapes:
rndInRange=: 13 : '(<./y)+?|>:-/y'
stst=: 13 : '(y i. 1),y i: 1'
circlePick=: 4 : '(wh{y)+(_1 1{~wh=. ?2)*?x'
release=: 4 : '(rndInRange _50 50+stst +./space) j. rndInRange _45 45+stst +./|:space'
[negative image]
This code offers a number of possibilities for improvement, but that's a lengthy topic on its own, so we'll deal with that in some separate pages, starting here and continuing here.
Jordan also brought up the possible usefulness of Voronoi diagrams - it would be great to have J code for this as it's a widely useful technique. One other possibility is to control the release locations by means of the "convex hull" around the cluster. This is a minimal enclosing boundary for a set of points and we have J code for calculating this in two dimensions. As it so happens, the next section of the meeting also leads us to convex hulls.
Show-and-tell
In this section, we re-visited a topic that had come up in a previous meeting: a programming contest in which the problem is to find the most closely-packed configuration of points on a lattice (integral co-ordinates) given some constraints. This contest, now over, is detailed here. Contests like these provide a good opportunity to showcase J if we can get some competitive solutions submitted before they're over.
This problem is stated as follows:
In example A, the area of the enclosing circle is 6.25 pi. In example B, it is 5 pi. Therefore, example B is the better packing.
One interesting feature (which led me off from the main problem) is the scoring system based on the "smallest enclosing circle". This is an interesting problem in its own right.
As it turns out, there is a rich literature and some fairly simple, well-defined solutions for this problem. In fact, it turns out there is a very efficient algorithm based on first finding the complex hull. This algorithm, taken from the source referenced here, is called "Applet's Algorithm" and is stated as follows:
1. For each vertex of H other than those of S, we compute the angle subtended by S. The minimum such angle occurs at vertex v, and the value of the minimum angle is "alpha".
* IF alpha >= 90 degrees THEN
o We are done!
o /* The minimal enclosing circle is determined by the diametric circle of S */
* IF alpha < 90 degrees THEN
o Goto step 2 below.
2. Since alpha < 90 degrees, check the remaining vertices of the triangle formed by S and v.
* IF no vertices are obtuse THEN
o We are done!
o /* The MEC is determined by the vertices of S and the vertex v */
3. If one of the other angles of the triangle formed by S and v is obtuse, then set S to be the side opposite the obtuse angle and repeat the main loop of the algorithm (i.e., Goto step 1 above).
I have not yet coded this up but we have the pieces to do this, not the least of which is a page of geometric functions, from Oleg Kobchenko, on the J wiki. I put together a crude demonstration of how this might be done as illustrated by the solutions shown for the two examples "A" and "B" shown in the "Point Packing Problem" statement above.
load '~Code/convexHull.ijs'
CvxHull=: (13 : 'j./"1 (],{.)~.y{~compl j./"1 y') NB. Crude approximation
pch=. (13 : '~.y{~compl j./"1 y') ptsa
cirpts=. circlePts ((xy3P pch);mean pch ptsDis xy3P pch);ptsa;15
pd 'show' [ pd (0{&.>cirpts) [ pd 'type point;pensize 4' NB. Circle
pd 'pensize 2' [ pd 'color RED'
pd 'show' [ pd (CvxHull |:1{&>cirpts) [ pd 'type line' NB. Hull
pd 'pensize 6' [ pd 'color GREEN'
pd 'show' [ pd (1{&.>cirpts) [ pd 'type point' NB. Points
Advanced topics
We studied some code on convex hulls written by Knuth in a "Literate Programming" style. However, without the proper "weaver", the examples were a bit perplexing.
Learning and teaching J
We discussed the possible benefits of a literate programming style of exposition of code and compared this with the preference of some forum members for tacit over explicit solutions for pedagogical purposes.
Scan of Meeting Notes
CategoryWorkInProgress CategoryWorkInProgress