kdtree.w,v

Lin-Kernighan heuristic for the TSP and minimum weight perfect matching

Because this process is naturally recursive, we use a separate recursive
procedure to do the dirty work.  But we must prime it with the right
arguments.  

It turns out that the range for the entire data set has already been
computed for us in the instance-reading function.  We just access those
values.

The set of values to be partitioned is specified using a closed-on-the-left
and open-on-the-right index range.

@@<Build the tree@@>=
coord = tsp->coord;
errorif(n<=0,"Need at least one point; instance has %d points",n);
E2_root = E2_build_helper(NULL,0,0,0,n,
	tsp->xmin,tsp->xmax,tsp->ymin,tsp->ymax);

@@ We need to declare the tree variable and the coordinate array variable.
The coordinates will be needed for queries later, which is why need to
keep it around on a more permanent basis.
@@<Module variables@@>=
static E2_node_t *E2_root;
static coord_2d *coord;

@@ The helper routine returns the
root of the subtree that it creates.

We are guaranteed to be given a positive number of points to partition,
so there will be at least one bucket created at this level.  So we create
it right away and fill in the common fields.

We increment the level number before testing to see whether to include
a bounding box 
because there is no sense in having a bounding box test at the root.

This invocation creates a bucket if there are few enough elements
(no more than |kd_bucket_cutoff|), or when both dimensions have already
been flattened.  The dimensions which have been flattened are recorded
in the bit vector |flat_dimen|.  Bit 0, with value $2^0=1$ is turned
on when the $x$ dimension is flattened, \ie, when one of our ancestors,
possibly this node, is an |eq_child| of a node whose cutting dimension
is the $x$ dimension.  Similarly, bit 1, with value $2^1=2$ is turned
on if the $y$ dimension is flattened.

@@<Module subroutines@@>=
static E2_node_t *
E2_build_helper(E2_node_t *parent,int flat_dimens,int level,
	int lo, int hi, 
	double xmin, double xmax, double ymin, double ymax) {
	E2_node_t *node = new_node();
	node->parent = parent;
#if !defined(KD_NO_HIDDEN_BIT)
	node->hidden = lo>=hi;
#endif
	if ( (++level % kd_bbox_skip)==0 ) {
		node->bbox = new_box();
		node->bbox->xmin = xmin;
		node->bbox->xmax = xmax;
		node->bbox->ymin = ymin;
		node->bbox->ymax = ymax;
		@@<Paranoid: check bounds@@>@@;
	} else {
		node->bbox = NULL;
	}
	if ( hi-lo <= kd_bucket_cutoff || flat_dimens == 0x03 ) {
		node->is_bucket = 1;
		node->f.e.lo = lo;
		node->f.e.hi = hi;
		@@<Fill other bucket fields@@>@@;
		@@<Make elements point to this bucket@@>@@;
	} else {
		node->is_bucket = 0;
		@@<Create three children@@>@@;
	}
	return node;
}

@@
It is possible for a node have no elements in it.  Such nodes must be buckets,
and have |hi==lo|.   The bounding box of such a bucket should not be trusted.

I could have chosen to avoid creating a node for an empty bucket, in which
case returning |NULL| would be the right answer.  However, I didn't do so
for the following reasons.

First, I didn't think of the possibility of empty buckets
before I wrote the code.  That's an honest answer.

Second, empty buckets are likely to be rare.  (I should measure this.)
So checking for a |NULL| pointer {\it every\/} time I want to visit a node
is likely to take more time overall than visiting an empty bucket.    
We can use the same simple code to cover both cases---empty and non-empty
buckets---it's just the termination test of the |for| loop.   
This shows again the benefit of treating zero as a first-class citizen
among numbers.

@@  Bottom-up searching only works if we can instantly jump to 
the right bucket.  We define an array that does this mapping for us.
@@<Make elements point to this bucket@@>=
{ int i;
for (i=lo;i<hi;i++) E2_point_to_bucket[perm[i]] = node;
}

@@ Of course, we need to declare this array.
@@<Module variables@@>=
static E2_node_t **E2_point_to_bucket;

@@ And we must allocate space for it.
@@<Allocate structures@@>=
E2_point_to_bucket = new_arr_of(E2_node_t *,n);

@@ And finally free it.
@@<Other deallocation@@>=
free_mem(E2_point_to_bucket);

@@ The recursive case is much more interesting.  It picks a cutting dimension,
finds a pivot, and then partitions.

If no dimensions have been flattened, then the
cut dimension is the one with the greatest range.
Otherwise, pick the nonflat dimension.

Exercise: Show that I could have replaced the entire |switch| statement
by the statement
|cutdimen = flat_dimens ? (xmax-xmin < ymax-ymin) : (2-flat_dimens);|.
Would it have been worth it?

@@<Create three children@@>=
{ int cutdimen=0; /* Satisfy GCC's dataflow analysis. */
switch ( flat_dimens ) {
case 0:
	if ( xmax-xmin > ymax-ymin ) cutdimen = 0;
	else cutdimen = 1;
	break;
case 1:
	cutdimen = 1;
	break;
case 2:
	cutdimen = 0;
	break;
case 3:	/* A valid |flat_dimens| value, but not here, so fall through. */@@;
default:
	errorif(1,"Invalid flat_dimens: %d", flat_dimens);
}
@@<Partition along dimension |cutdimen|@@>@@;
}

@@ Partitioning involves finding a pivot and then splitting the points into
three classes.

@@<Partition along dimension |cutdimen|@@>=
{ int p;	/* The pivot's index in |perm|. */
int a,b,c,d;	/* Cursors for partitioning. */
double exl,exh,lxl,lxh,gxl,gxh;	/* |x| extremes. */
double eyl,eyh,lyl,lyh,gyl,gyh;	/* |y| extremes. */
@@<Find the pivot@@>@@;
@@<Partition around |perm[p]|@@>@@;
@@<Build the subtrees@@>@@;
}

@@
In the Quicksort application, Bentley and McIlroy settled on the following
@@^Bentley, Jon Louis@@>
algorithm for finding a pivot.  For large inputs, it uses Tukey's `ninther':
the median of the three medians of three.  For smaller inputs, a median
of three is used.  

Now, 
their pivot-picking algorithm
may not be suitable here
because their application differs from ours in the following ways.
First, our cost model differs: we compare coordinates, which is fast,
while their comparison is a function-call away.  
Second, they are doing a pure sort; we too are sorting, but in $k$ dimensions.
Does this change matters?
Third, their machine is a VAX-8550; ours is likely to be a modern RISC machine.
Fourth, theirs is a library function, so side-effecting the random number
generator is a bad thing; we don't care about that: how much better or
slower is randomization? 

However, until someone gives convincing evidence that there are better choices
for the pivot, I will stick to the Bentley and McIlroy solution.
Before you run off and try to find a better solution, here's how
their algorithm stacks up:
to sort $n$ random 30-bit integers, they required $1.094n\lg n-0.74n$
comparisons.  (This is for $n$ ranging over each power of two from 128 to
65536. The function was found by a weighted least-squares regression fit.)

If you're thinking `pick the true median', then your thoughts
have already been anticipated.  In the Quicksort application, finding
the true median, although taking only linear time, is too costly.

@@d val(a) (coord[perm[(a)]].x[cutdimen])
@@d valx(a) (coord[perm[(a)]].x[0])
@@d valy(a) (coord[perm[(a)]].x[1])
@@d med3(a,b,c) (val(a)<val(b) ? (val(b)<val(c) ? b : val(a)<val(c) ? (c):(a)) @@;
							  : (val(b)>val(c) ? b : val(a)>val(c) ? (c):(a)) )

@@<Find the pivot@@>=
p = (lo+hi)/2;	/* Small arrays, middle element. */
if ( hi-lo > 7 ) {
	int p1=lo,pn=hi-1;
	if ( hi-lo > 40 ) {	/* Big arrays, pseudomedian of 9. */
		int s = (hi-lo)/8;
		p1 =  med3(p1,p1+s,p1+s+s);
		p =  med3(p-s,p,p+s);
		pn =  med3(pn-s-s,pn-s,pn);
	}
	p = med3(p1,p,pn);	/* Mid-size, median of 3. */
}

@@ Partitioning the elements into three classes means more than just finding
which elements belong to which class---it also means moving the elements
so that the members of each class is contiguous.  In addition, we'll
move all the lesser elements to the front, the greater elements to the
back, and the equal elements to the middle.  
This problem is equivalent to Dijkstra's `Dutch National Flag' problem,
@@^Dijkstra, Edsger W.@@>
after the Dutch flag, which is three horizontal stripes of red, white, and
blue.  (See {\tt http://www.cesi.it/flags/nl.html} for a rendition
of the Dutch National Flag.)

The invariant maintained is: equal, lesser, unknown, greater, equal.
The boundary indices are |lo|, |a|, |b|, |c|, |d|, and |hi-1|.

We're also careful to find the ranges over each set of elements.

@@d swapint(J,K)@@+ {@@+ int t=J;@@+ J=K;@@+ K=t;@@+ }
		/* Assumes |t| is not free in |J| or |K| */

@@<Partition around |perm[p]|@@>= 
@@<Initialize bounds@@>@@;
a=b=lo; c=d=hi-1; 
{double v = val(p), diff; 
node->f.i.cutdimen = cutdimen;  
node->f.i.cutvalue = v;
for(;;) {
	while( b<=c && (diff=val(b)-v)<=0.0 ) {
		if ( diff==0.0 ) {
			@@<Update bounds for equal elements at |b|@@>@@;
			swapint(perm[a],perm[b]);
			a++;
		} @@+ else @@+ {@@+  
			@@<Update bounds for lesser elements@@>@@+
		}
		b++;
	}
	while( c>=b && (diff=val(c)-v)>=0.0 ) {
		if ( diff==0.0 ) {
			@@<Update bounds for equal elements at |c|@@>@@;
			swapint(perm[d],perm[c]);
			d--;
		} @@+ else @@+ {@@+  
			@@<Update bounds for greater elements@@>@@+
		}
		c--;
	}
	if ( b>c ) break;
	swapint(perm[b],perm[c]);
	@@<Update bounds for lesser elements@@>@@;
	@@<Update bounds for greater elements@@>@@;
	b++; c--;
}
@@<Move equals to middle@@>@@;
}
@@<Paranoid: check the partitioning@@>@@;

@@
@@<Initialize bounds@@>=
exl=lxl=gxl=xmax;	
exh=lxh=gxh=xmin;	
eyl=lyl=gyl=ymax;	
eyh=lyh=gyh=ymin;	

@@
@@d min(x,y) ((x)<(y)?(x):(y))
@@d max(x,y) ((x)>(y)?(x):(y))
@@<Update bounds for equal elements at |b|@@>=
			exl = min(exl,valx(b));
			exh = max(exh,valx(b));
			eyl = min(eyl,valy(b));
			eyh = max(eyh,valy(b));

@@
			@@<Update bounds for lesser elements@@>=
			lxl = min(lxl,valx(b));
			lxh = max(lxh,valx(b));
			lyl = min(lyl,valy(b));
			lyh = max(lyh,valy(b));

@@
			@@<Update bounds for equal elements at |c|@@>=
			exl = min(exl,valx(c));
			exh = max(exh,valx(c));
			eyl = min(eyl,valy(c));
			eyh = max(eyh,valy(c));

@@
			@@<Update bounds for greater elements@@>=
			gxl = min(gxl,valx(c));
			gxh = max(gxh,valx(c));
			gyl = min(gyl,valy(c));
			gyh = max(gyh,valy(c));

@@
@@<Move equals to middle@@>=
{
int s,l,h;
s=min(a-lo,b-a);	/* Move first set of equals to middle. */
for(l=lo,h=b-s;s;s--) {swapint(perm[l],perm[h]);@@+l++;h++;}
s=min(d-c,hi-1-d);	/* Move last set of equals to middle. */
for(l=b,h=hi-s;s;s--) {swapint(perm[l],perm[h]);@@+l++;h++;}
}

@@ Now that the elements are partitioned, we recurse to build the subtrees.

It is possible to pick a sequence of bad pivots, causing many of the
segments to be large.  This would not only slow things down, but it would
also force us to use a lot of stack space.  

There are two ways to get around this.
First, one might always choose a pivot that guarantees a good partition,
for example the median; medians can be found in linear time (INSERT REFERENCE).
Second, one might always process the smaller segments first and
use tail recursion to process the remaining largest segment.  This guarantees
that at most $\log_2 n$ activation records are present on the stack
@@^tail recusion@@>
@@^stack overflow@@>
at any one time.

However, until this problem shows up,  
I will use the simplest solution: blind recursion.
Note that Quicksort itself runs in $O(n\log n)$ time on average.
For example, see C.~A.~R.~Hoare's INSERT REFERENCE, or Robert Sedgewick's
PhD thesis (INSERT REFERENCE), or the perhaps more readily available
Knuth (INSERT REFERENCE), section ??.
@@^Hoare, Charles Anthony Richard@@>
@@^Knuth, Donald Ervin@@>
@@^Sedgewick, Robert@@>

%%%
%%%A common trick in Quicksort implementations is to recurse on the 
%%%smaller-size array segment first.  This leaves the larger segment waiting
%%%in the activation record on the stack.   When the smaller segment is done,
%%%then we sort the larger segment.
%%%This strategy ensures that
%%%there are at most $\log_2 n$ activation records on the stack, protecting
%%%us from stack overflow.
%%%@@^stack overflow@@>
%%%
%%%Interestingly, many people know this trick, but a significant fraction
%%%of them get it backwards: they recurse on the {\it larger\/} segment
%%%first.  Knuth (CHECK THIS) comments on this phenomenon, and wags his finger
%%%@@^Knuth, Donald Ervin@@>
%%%appropriately.
%%%
%%%Here, we have three segments on which we recurse.  A moment's thought tells
%%%us 
%%%to recurse on the smallest first, the next smallest second, and the largest
%%%last.  This still guarantees that there are at most $\log_{2} n$ activation
%%%records on the stack: because the smallest segment is at most one third
%%%of the remaining elements, and the middle segment is at most half the
%%%remaining elements.
%%%
%%%We only need the largest segment last in order to guarantee this logarithmic
%%%bound.  Arranging to meet  
%%%this simpler condition isn't much less work: one fewer
%%%integer comparison in the worst case.
%%%I didn't realize this until after I coded the following, which arranges
%%%the
%%%stronger small-medium-large order, so I left it in.
%%%
%%%

@@<Build the subtrees@@>=
#if defined(KD_BUILD_SMALLEST_SEGMENT_FIRST)
{ int l = b-a, m = hi-(d-c)-lo+b-a, h = d-c;
if ( l <= m ) {
	if ( m <= h ) { @@<Build $<$ tree@@> @@+ @@<Build $=$ tree@@> @@+ 
		@@<Build $>$ tree@@> }
	else if ( l <= h ) { @@<Build $<$ tree@@> @@+ @@<Build $>$ tree@@> @@+ @@<Build $=$ tree@@> }
	else { @@<Build $>$ tree@@> @@+ @@<Build $<$ tree@@> @@+ @@<Build $=$ tree@@> }
} else {
	if ( l <= h ) { @@<Build $=$ tree@@> @@+ @@<Build $<$ tree@@> @@+ @@<Build $>$ tree@@> }
	else if ( m <= h ) { @@<Build $=$ tree@@> @@+ @@<Build $>$ tree@@> @@+ @@<Build $<$ tree@@> }
	else { @@<Build $>$ tree@@> @@+ @@<Build $=$ tree@@> @@+ @@<Build $<$ tree@@> }
}
}
#else
@@<Build $<$ tree@@>@@;
@@<Build $=$ tree@@>@@;
@@<Build $>$ tree@@>@@;
#endif

@@ To build the subtree for the smaller elements, we call ourselves 
recursively.  We  pass this node as the parent; we propagate the information
on which dimensions have been flattened, and what level the subtree is at; 
we also pass the array bounds for the segment and the spatial 
bounding box of those
points.

@@<Build $<$ tree@@>=
node->f.i.lo_child =
	E2_build_helper(node,flat_dimens,level,lo,lo+b-a,lxl,lxh,lyl,lyh);

@@ @@<Build $>$ tree@@>=
node->f.i.hi_child =
	E2_build_helper(node,flat_dimens,level,hi-(d-c),hi,gxl,gxh,gyl,gyh);

@@ When we recurse on the equal elements, we additionally 
tell the new invocation
that this dimension has already been flattened.
For $k\in\{0,1\}$, $2^k=k+1$, so we use |cutdimen+1| in place of an 
exponentiation.

@@<Build $=$ tree@@>=
node->f.i.eq_child = 
E2_build_helper(node,flat_dimens|(cutdimen+1),level,lo+b-a,hi-(d-c),exl,exh,eyl,eyh);


@@*Hiding and unhiding.
Part of the meaning of a semidynamic point set is that we may hide and unhide
points at will.  
A point becomes invisible to queries when it is hidden.  It becomes
visible again when it is unhidden.

First we will describe how to hide and unhide single points.
Later we will show how to hide and unhide all the points at once.

@@ Inside a bucket, we hide a node by moving it to the end of the bucket, and