kctree.cpp

高效的k-means算法实现

//----------------------------------------------------------------------
//	File:		KCtree.cc
//	Programmer:	David Mount
//	Last modified:	08/10/2005
//	Description:	Basic methods for kc-trees.
//----------------------------------------------------------------------
// Copyright (C) 2004-2005 David M. Mount and University of Maryland
// All Rights Reserved.
// 
// This program is free software; you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation; either version 2 of the License, or (at
// your option) any later version.  See the file Copyright.txt in the
// main directory.
// 
// The University of Maryland and the authors make no representations
// about the suitability or fitness of this software for any purpose.
// It is provided "as is" without express or implied warranty.
//----------------------------------------------------------------------

#include "KCtree.h"			// kc-tree declarations
#include "KMfilterCenters.h"		// center set structure
#include "KMrand.h"			// random number includes

//----------------------------------------------------------------------
//  Declaration of local utilities.  These are used in getNeighbors().
//----------------------------------------------------------------------
static int closestToBox(		// get closest point to box center
    KMctrIdxArray	cands,			// candidates for closest
    int			kCands,			// number of candidates
    KMorthRect		&bnd_box);		// bounding box of cell

static bool pruneTest(			// test whether to prune candidate
    KMcenter		cand,			// candidate to test
    KMcenter		closeCand,		// closest candidate
    KMorthRect		&bnd_box);		// bounding box

static void postNeigh(			// assign neighbors to center
    KCptr		p,			// the node posting
    KMpoint		sum,			// the sum of coordinates
    double		sumSq,			// the sum of squares
    int			n_data,			// number of points
    KMctrIdx		ctrIdx);		// center index

//----------------------------------------------------------------------
//  KCtree constructors
//	There is a skeleton kc-tree constructor which does (almost)
//	all the initialization, without actually building the tree.
//	The arguments are essentially the same as the constructor
//	for the kc-tree.  The last argument (pi) is an optional array
//	for the point indices.  Normally, this will be NULL (meaning
//	that the constructor should initialize the array of indices),
//	but for the load constructor the point indices will be set
//	created and passesd through this argument.
//----------------------------------------------------------------------

void KCtree::skeletonTree(		// construct skeleton tree
    KMdataArray		pa,		// point array (with at least n pts)
    int			n,		// number of points
    int			dd,		// dimension
    int			n_max,		// maximum number of points (optional)
    KMpoint		bb_lo,		// bounding box low point (optional)
    KMpoint		bb_hi,		// bounding box high point (optional)
    KMdatIdxArray	pi)		// point indices (optional)
{
    					// initialize basic elements
    dim = dd;				// dimension
    n_pts = n;				// number of points
    if (n_max < n) n_max = n;		// max_pts must be >= n
    max_pts = n_max;			// set max_pts

    if (pa == NULL) {			// no points supplied?
    	kmError("Points must be supplied to construct tree.", KMabort);
    }
    pts = pa;				// initialize point array

    if (pi == NULL) {			// point indices provided?
	pidx = new KMdataIdx[max_pts];	// no, allocate them
	for (int i = 0; i < n; i++)	// initialize to identity
	    pidx[i] = i;
    }
    else pidx = pi;			// yes, just use them

    			
    if (bb_lo == NULL || bb_hi == NULL) // boundng box fully specified?
	kmEnclRect(pa, pidx, n, dd, bnd_box);	// no, construct from points
					// save bounding box
    if (bb_lo != NULL)			// if lower point given, then use it
    	bnd_box.lo = kmAllocCopyPt(dd, bb_lo);

    if (bb_hi != NULL)			// same for upper point
    	bnd_box.hi = kmAllocCopyPt(dd, bb_hi);

    root = NULL;			// no associated tree yet
}

//----------------------------------------------------------------------
//  BasicGlobals - global variables
// 	To prevent long argument lists in a number of the tree traversal
// 	programs, we store a number of common global variables here.
// 	In the case of construction, these are initialized before
// 	calling buildKcTree.  They are used in getNeighbors() and by
// 	sampleCtr().
//----------------------------------------------------------------------

int		kcDim;			// dimension of space
int		kcDataSize;		// number of data points
KMdataArray	kcPoints;		// data points

//----------------------------------------------------------------------
//  initBasicGlobals - initialize basic globals
//----------------------------------------------------------------------

static void initBasicGlobals(		// initialize basic globals
    int			dim,			// dimension
    int			data_size,		// number of data points
    KMdataArray		data_pts)		// data points
{
    kcDim = dim;
    kcDataSize = data_size;
    kcPoints = data_pts;
}

//----------------------------------------------------------------------
// kc-tree constructor
//	This is the main constructor for kc-trees given a set of
//	points.  It first builds a skeleton tree, then computes the
//	bounding box, and then invokes buildKcTree() to actually
//	build the tree.  It passes in the appropriate splitting
//	routine.
//
//	The constructor has a number of optional arguments.
//
//	n_max		Max number of points.  (default: n).
//	bb_lo, bb_hi	Bounding box low and high points (default:
//			compute bounding box from points).
//
//	As long as the number of points is nonzero, or if a bounding
//	box is provided, then the constructor will build a tree with
//	at least one node (even an empty leaf).  Otherwise, it returns
//	with a null tree.
//
//	Under the current implementation, point insertion generates
//	leaf nodes with a bucket size of 1.  If the requested bucket
//	size is higher, and it looks like insertion is requested, then
//	we generate a warning message.
//----------------------------------------------------------------------

KCtree::KCtree(			// construct from point array
    KMdataArray		pa,		// point array (with at least n pts)
    int			n,		// number of points
    int			dd,		// dimension
    int			n_max,		// maximum number of points (optional)
    KMpoint		bb_lo,		// bounding box low point (optional)
    KMpoint		bb_hi) :	// bounding box high point (optional)
    bnd_box(dd)				// create initial bounding box
{
    					// set up the basic stuff
    skeletonTree(pa, n, dd, n_max, bb_lo, bb_hi, NULL);
    initBasicGlobals(dd, n, pa);	// initialize globals

    root = buildKcTree(pa, pidx, n, dd, bnd_box);

    int ignoreMe1;			// ignore results of call
    KMpoint ignoreMe2;
    double ignoreMe3;
    					// compute sums
    root->makeSums(ignoreMe1, ignoreMe2, ignoreMe3);
    assert(ignoreMe1 == n);		// should be all the points
}

//----------------------------------------------------------------------
//  buildKcTree - recursive procedure to build a kc-tree
//	Builds a kc-tree for points in pa as indexed through the array
//	pidx[0..n-1] (typically a subarray of the array used in the
//	top-level call).  This routine permutes the array pidx, but does
//	not alter pa[].
//
//	The construction is based on a standard algorithm for
//	constructing the kc-tree (see Friedman, Bentley, and Finkel,
//	``An algorithm for finding best matches in logarithmic expected
//	time,'' ACM Transactions on Mathematical Software, 3(3):209-226,
//	1977).  The procedure operates by a simple divide-and-conquer
//	strategy, which determines an appropriate orthogonal cutting
//	plane (see below), and splits the points.  When the number of
//	points falls below 1, we simply store the points in a leaf
//	node's bucket.
//
//	This procedure selects a cutting dimension and cutting value,
//	partitions pa about these values, and returns the number of
//	points on the low side of the cut.
//
//	Note that this procedure is not only used for constructing full
//	trees, but is also used by the insertion routine to rebuild a
//	subtree.
//	
//----------------------------------------------------------------------

KCptr KCtree::buildKcTree(	// recursive construction of kc-tree
    KMdataArray		pa,		// point array
    KMdatIdxArray	pidx,		// point indices to store in subtree
    int			n,		// number of points
    int			dim,		// dimension of space
    KMorthRect		&bnd_box)	// bounding box for current node
{
    if (n <= 1) {			// n small, make a leaf node
	return new KCleaf(dim, bnd_box, n, pidx); 
    }
    else {				// n large, make a splitting node
	int cd;				// cutting dimension
	KMcoord cv;			// cutting value
	int n_lo;			// number on low side of cut
	KCptr lo, hi;			// low and high children

					// invoke splitting procedure
	sl_midpt_split(pa, pidx, bnd_box, n, dim, cd, cv, n_lo);

	KMcoord lv = bnd_box.lo[cd];	// save bounds for cutting dimension
	KMcoord hv = bnd_box.hi[cd];

	bnd_box.hi[cd] = cv;		// modify bounds for left subtree
	lo = buildKcTree(		// build left subtree
		pa, pidx, n_lo,		// ...from pidx[0..n_lo-1]
		dim, bnd_box);
	bnd_box.hi[cd] = hv;		// restore bounds

	bnd_box.lo[cd] = cv;		// modify bounds for right subtree
	hi = buildKcTree(		// build right subtree
		pa, pidx + n_lo, n-n_lo,// ...from pidx[n_lo..n-1]
		dim, bnd_box);
	bnd_box.lo[cd] = lv;		// restore bounds

					// create the splitting node
	KCsplit *ptr = new KCsplit(dim, bnd_box, cd, cv,
						lv, hv, lo, hi);
	return ptr;			// return pointer to this node
    }
} 

//----------------------------------------------------------------------
//  cellMidpt - return bounding box midpoint
//	The result is returned by modifying the argument.
//  getPoint - return point in a leaf cell
//	If the leaf cell has no point, then NULL is returned.
//	(This cannot be inlined in KCtree.h because of reference
//	to thePoints.)
//----------------------------------------------------------------------

void KCnode::cellMidpt(	// compute cell midpoint
    KMpoint	pt)			// the midpoint (returned)
{
    for (int d = 0; d < kcDim; d++) {		// compute box midpoint
	pt[d] = (bnd_box.lo[d] + bnd_box.hi[d])/2;
    }
}

KMpoint KCleaf::getPoint()		// get data point
{  return (n_data == 1 ? kcPoints[bkt[0]] : NULL);  }


//----------------------------------------------------------------------
//  kc-tree make sums (part of constructor)
//	Computes the sums of points for each node of the kc-tree,
//	and the sums of squares (the sums of dot products of each
//	point with itself).  These values are returned through the
//	arguments.
//
//	The sum of points is assumed to have been allocated as part
//	of the constructor, and hence already exists.
//----------------------------------------------------------------------

void KCsplit::makeSums(
    int			&n,			// number of points (returned)
    KMpoint		&theSum,		// sum (returned)
    double		&theSumSq)		// sum of squares (returned)
{
    assert(sum != NULL);			// should already be allocated