kdtree.w,v

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

Revision 1.158  1998/02/13 21:16:55  neto
Made prototypes and function declarations match w.r.t. const.
,

Revision 1.157  1997/10/25 21:44:06  neto
Finished coding quadrant stuff (I think). Must compile and test.

Revision 1.156  1997/10/25 19:49:11  neto
Compact the code.

Revision 1.155  1997/10/18 15:23:37  neto
I had missed a swapint conversion.  Now it compiles.

Revision 1.154  1997/10/18 14:40:25  neto
Added E2 supports()

Revision 1.153  1997/10/18 14:31:05  neto
Made all floats into doubles.
Checked for CEIL 2D compliance, and added comments in key places.
Minor editing.

Revision 1.152  1997/10/17  21:48:51  neto
Removed CUTDIMEN and THISDIMEN kludge by recasting the coord 2d structure
as an array instead of two named fields.

Revision 1.151  1997/10/17  20:47:35  neto
Move short fields to end of E2 node struct.

Revision 1.150  1997/09/27 18:06:44  neto
Fixed RCS log behaviour.

Revision 1.149  1997/08/15  20:18:25  neto
Added Index major section.

Revision 1.148  1997/05/16  18:11:41  neto
Break locks by david and neto.
Include <config.h> and "lkconfig.h"

Revision 1.147  1997/01/21  21:55:55  david
Added standard copyright notice by including copyrt.w

Revision 1.146  1997/01/21  17:42:14  david
Make showing of partitioning optional.

Revision 1.145  1997/01/17  21:15:44  neto
Fixed a printing leak.

Revision 1.144  96/12/13  15:08:17  neto
Mention hypot

Revision 1.143  96/12/02  15:26:06  neto
Added copyright notice.

Revision 1.142  96/11/11  16:36:44  neto
fixed a spelling mistake.

Revision 1.141  96/09/03  15:14:03  neto
Macro for ML programming language name.

Revision 1.140  96/08/16  13:04:53  neto
Added fixincludes.

Revision 1.139  96/08/16  12:42:26  neto
Converted putchar to printf.  Otherwise, I'd never get a prototype
for SunOS' \_flusbuf.

Revision 1.138  96/08/15  14:20:27  neto
Fixed the prototype for E2\_rnn

Revision 1.137  96/08/15  12:54:30  neto
Make it pass -Wall

Revision 1.136  96/07/30  17:15:34  neto
Fixed editing bug with respect to recurse if not hidden in E2\_rnn.

Revision 1.135  96/07/30  13:03:33  neto
Fixed to compile.

Revision 1.134  96/07/30  12:59:52  neto
Added KD\_BUILD\_SMALLEST\_SEGMENT\_FIRST

Revision 1.133  96/07/29  17:08:56  neto
Fixed to compile.

Revision 1.132  96/07/29  16:19:56  neto
Added *\_rcs\_id.
Made sure RCS log is activated within this file.

Revision 1.131  96/07/26  16:19:54  neto
Report experiments with hidden bit.

Revision 1.130  96/07/26  15:41:44  neto
Added KD\_NO\_HIDDEN\_RNN\_TEST

Revision 1.129  96/07/26  13:44:55  neto
Added compile-time-option KD\_NO\_HIDDEN\_BIT

Revision 1.128  96/07/26  13:21:36  neto
Recurse in rnn ony if not hidden.

Revision 1.127  96/07/26  12:45:03  neto
Don't be preachy about recursing on largest last.

Revision 1.126  96/07/25  14:02:54  neto
Fixed some minor editing of english.
Replaced INFINITY by a computed (floating point!) upper bound.
Isolated debugging output and made it subject to verbose value and
compile-time choices.

Revision 1.125  96/07/25  12:09:27  neto
More info to catch dsj1000 bug.

Revision 1.124  96/07/24  15:57:27  neto
Removed 117 debugging output.
Made bbox checkig optional.

Revision 1.123  96/07/24  15:48:46  neto
Fixed the bbox bug.  I wasn't updating the boboounding box
on the b and c elemnts that were swapped.

Revision 1.122  96/07/24  13:51:38  neto
More debugging outupt for failure kcase on lin318

Revision 1.121  1996/07/15  22:28:48  david
Simplified the code a bit.
It fails on lin318.tsp

Revision 1.120  1996/07/15  21:18:21  david
Bounding box every 1/kd\_bucket\_skip fraction, not 1-that fraction.
Big bug removed: bounding box initialization during Dutch National Flag
problem is now fixed.  (Had max's at infinity instead of effectively
negative infinity.)
Factored the nn search in a bucket.  The general case now is free
from checking that it doesn't return the seed.  It should be faster.
Added PostScript debugging output to draw the buckets.

Revision 1.119  1996/07/15  18:33:24  david
Fixed a name of a section.
Fixed bucket search for nearer neighbour.

Revision 1.118  1996/07/15  16:00:08  david
Compiles now.  I had to add |f.i.| to many of the field names.
I also had to add $E2\_nn\_dist$ back in for the bounding box test.

Revision 1.117  1996/07/15  15:23:55  david
Give up on FRNN, ball searching for now.
We don't need them for Greedy / MF heuristic.

Revision 1.116  1996/07/05  21:02:33  david
Continued a bit on fixed radius nearest neighbours.
Added more about stack overflow.  Referenced Peyton Jones' Implementing
Functional Languages.

Revision 1.115  1996/07/05  20:34:02  david
I forgot to mention the move to using only squared distances in
nearest neighbour search.

Revision 1.114  1996/07/05  20:31:31  david
Changed $E2\_tree$ to $E2\_root$.
Implemented $unhide\_all$ and $hide\_all$.
Got started on comments about fixed radius nearest neighbour.

Revision 1.113  1996/07/05  17:13:32  david
Added some index references to people.
changed is\_internal to is\_bucket.
Worked on nearest neighbours computation.

Revision 1.112  1996/07/04  20:07:38  david
Clarify comments in first part of the document.
Make the recursive building more CWEB-like.
Finished the hiding/unhiding part, but I haven't tried compiling it.

Revision 1.111  96/06/28  14:00:54  neto
Added comments about persistence.
Fixed too-long-line for cweave.
Worked on hiding/unhiding
Added to nn.

Revision 1.110  96/06/28  12:58:29  neto
Ensure theat recursing won't break the stack.
Start work on nearest neighbour query.

Revision 1.109  96/06/26  12:15:42  neto
Made coord a static variable.

Revision 1.108  96/06/26  12:08:19  neto
Be careful with macro parameters.  Parenthesize them.

Revision 1.107  96/06/26  12:05:44  neto
This file was added after I went to 100 and above.

Revision 1.6  96/06/26  12:02:28  neto
Fixed my in-file log.  (I edited, saved, and committed the revision,
but didn't reload before the next commit...)

Revision 1.5  96/06/26  11:38:41  neto
Removed some printing clutter.

Revision 1.4  96/06/26  11:35:41  neto
Most importantly, I fixed the construction bug.  The code that was
moving equal elements to the middle did not handle the left end properly.
It was 0-based instead of being based at the real lower end: lo.

I added comments about what to do with an empty bucket, and why.

I added better paranoid checking diagnostics.

Revision 1.3  96/06/25  17:10:15  neto
Added missing break.
Fixed a swap.  Duh.
More precise output.

Revision 1.2  96/06/25  16:41:37  neto
Fixed typo in med3 computation.
Added better debugging output for paranoid.
Fails paranoid check on lin105.tsp

Revision 1.1  96/06/24  17:16:26  neto
Initial revision

}

@@s near void
@@s far int

@@*$k$-dimensional search structures.
Many geometrical algorithms can be made to run faster in practice if we
use a $k$-dimensional search structure.  Bentley (INSERT REFERENCE) has
@@^Bentley, Jon Louis@@>
shown this to be true of many algorithms related to the TSP.  

Bentley (INSERT REFERENCE) uses a $k$-d search tree over 
\term{semi-dynamic point sets}@@^semi-dynamic point sets@@>:
these are sets of points where all the points are known in advance,
and individual points may be hidden and unhidden from the queries.  

Bentley actually describes the hiding and unhiding as \term{deletion} and 
\term{undeletion}.
But as the points still remain in the structure and may be recovered, I
prefer the concept of hiding; it seems a more accurate description of
what is going on.  This is related in some ways to the notion of
persistency.
In general, we say a data structure is \term{(fully) persistent}
@@^persistence, full@@>
if we can completely recover the visible part of any of 
its old states; we say it is 
\term{partially persistent} if we may query any of its old states, but
@@^persistence, partial@@>
may update only its latest state.  Persistency 
is a hot topic in data structures, and has led to some interesting
and intuitive new algorithms (INSERT REFERENCE).
@@^Tarjan, Robert Endre@@>
@@^Sleator@@>

This module implements $k$-d trees for semi-dynamic point sets.

Actually, only one tree at a time is supported by this module.  
This is all I need (so far)
in the TSP application.  This scheme is also simpler to implement,
and it saves a lot of
parameter passing time (or indirection time) as compared with fully 
object-oriented code.  Besides, the code that Bentley describes
for his experiments in (INSERT REFERENCE) only supports one tree at a time.
%Some of his parameters and maybe even conclusions might change if he
%wrote fully object-oriented code.
%%%% Is that true?  Check C++ meaning for |static| storage class specifier
%%%% for a class member.

The outline of this module is as follows:

@@c
#include <config.h>
#include "lkconfig.h"
@@<System headers@@>@@;
@@<Early module headers@@>@@;
@@<Module headers@@>@@;

@@<Early module type definitions@@>@@;
@@<Module type definitions@@>@@;
@@<Module variables@@>@@;
@@<Global variables@@>@@;
@@<Module subroutines@@>@@;
@@<Subroutines@@>@@;
const char *kdtree_rcs_id = "$Id: kdtree.w,v 1.164 1998/10/09 21:34:29 neto Exp neto $";

@@ The exported interface is contained in the \file{kdtree.h} header file,
which has the following form.

@@(kdtree.h@@>=
#if !defined(_KDTREE_H_)
#define _KDTREE_H_
extern const char *kdtree_rcs_id;
@@<Exported type definitions@@>@@;
@@<Exported variables@@>@@;
@@<Exported macros@@>@@;
@@<Exported subroutines@@>@@;
#endif

@@ To ensure consistency between the interface and the implementation,
we include our own header.
@@<Module headers@@>=
#include "kdtree.h"

@@*Interface.  
This section presents the interface to the $k$-d tree as specified by
Bentley. I have altered the names of the procedures to fit within my
naming scheme for modules.  
@@^Bentley, Jon Louis@@>

Currently this module only supports 2-dimensional point
sets with the Euclidean metric.
The names of the procedures all end with $E2$ to reflect this fact.

The interface is as follows:

|void E2_create(tsp_instance_t *tsp)| 
creates a tree based on the TSP instance |tsp|.  
The structure is defined
in the \module{READ} module.

|void E2_destroy()| destroys the search tree.

|void E2_hide(int i)| hides point |i| from future queries (until the
next |unhide| operation).

|void E2_unhide(int i)| makes point |i| visible to future queries (until the
next |hide| operation).

|void E2_hide_all(void)| hides all points.  

|void E2_unhide_all(void)| hides all points.

|int E2_nn(int i)| returns the index of an unhidden point that is nearest
to point |i|, but is not |i| itself. If there is none, $-1$ is returned.

|int E2_nn_quadrant(int i, int q)| returns the index of an unhidden point 
in relative quadrant $q$ that is nearest
to point |i|, and is not |i| itself. If there is none, $-1$ is returned.
For example, quadrant 1 relative to point |i| is the set of all points
to the northeast of |i|.

For efficiency reasons, we also have 
bulk versions of the nearest neighbour searches.
Function |int E2_nn_bulk(int i, int num_wanted, kd_bin_t *buffer)|
tries to find the nearest |num_wanted| neighbours to city |i|, and
writes them into |buffer|.
It returns the number of vertices actually found.
Naturally, the number found is as close as possible but no
greater than |num_wanted|.
Function
|int E2_nn_quadrant_bulk(int i, int num_wanted, kd_bin_t *buffer, int qmask)| 
does the same for quadrant-based searches.
In neither case are the buffers ordered in any particular way;
at least we don't promise so.
To help sort out the buffer we also supply comparison functions 
|kd_bin_cmp_increasing|
and
|kd_bin_cmp_decreasing|
suitable for passing to |qsort| and related functions.

To access the city field within a bin, we export the macro |kd_bin_city|.
To access the field within a bin that is monotonic with the distance to the
city, we export the macro |kd_bin_monotonic_len|.

|void E2_frnn(int i, double rad, void (*proc)(int j))| is 
a 
\term{fixed-radius nearest neighbour search}@@^fixed-radius nearest neighbour search@@>.
It calls |proc| on all points |j| that lie within |rad| units of point |i|.

The following two procedures are used to implement 
\term{sphere of influence}@@^sphere of influence@@>
computations.

|void E2_set_radius(int i, double r)| sets the radius of the 
active sphere around
point |i| to |r| units.

|void E2_ball_search(int i, void (*proc)(int j))| calls |proc| on every
point |j| in whose active sphere point |i| lies.
% The girl the boy the dog bit hit cried.
Get that?  Let's try it a different way:  this procedure determines which
points |j| have point |i| somewhere within their own (|j|'s) 
active sphere; it calls
|proc| on every such |j|.


@@<Exported subroutines@@>=
void E2_create(tsp_instance_t *tsp);
void E2_destroy(void);
void E2_hide(int i);
void E2_unhide(int i);
void E2_hide_all(void);
void E2_unhide_all(void);
int E2_nn_func(int i);
int E2_nn_quadrant(int i, const int qmask);
int E2_nn_bulk_func(int i, int num_wanted, kd_bin_t *buffer);
int E2_nn_quadrant_bulk(int i, int num_wanted, kd_bin_t *buffer, int qmask);
int kd_bin_cmp_increasing(const void *a, const void *b);
int kd_bin_cmp_decreasing(const void *a, const void *b);
void E2_frnn(int i, double rad, void (*proc)(int j));
void E2_set_radius(int i, double r);
void E2_ball_search(int i, void (*proc)(int j));


@@ It is also useful to know when a 2-d trees may be used with a particular
TSP instance.

@@<Exported subroutines@@>=
int E2_supports(tsp_instance_t *tsp);


@@  These routines support TSP instances with |EUC_2D|, |CEIL_2D|, and
|ATT| cost functions.

Now, |EUC_2D| and |CEIL_2D| are obvious choices, since they use essentially
just 2-d Euclidean distance.  (|EUC_2D| rounds to the nearest integer,
and |CEIL_2D| always rounds toward postivie infinity.)

|ATT| uses a scaled Euclidean distance; it still has 2-d coordinates.

Now, the nearest neighbour searches use only 
relative ordering based on raw Euclidean distance.  
Since the scaling and rounding are monotonic functions, we can support
all three.

@@<Subroutines@@>=
int 
E2_supports(tsp_instance_t *tsp) 
{
	switch(tsp->edge_weight_type) {