kmrand.cpp

高效的k-means算法实现

//----------------------------------------------------------------------
//	File:		KMrand.cpp
//	Programmer:	Sunil Arya and David Mount
//	Last modified:	05/14/04
//	Description:	Routines for random point generation
//----------------------------------------------------------------------
// 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 "KMrand.h"			// random generator declarations

#ifdef WIN32				// Visual C++ (no srandom/random)
void srandom(unsigned int seed) { srand(seed); }
long random(void) { return long(rand()); }
#endif

//----------------------------------------------------------------------
//  Globals
//----------------------------------------------------------------------
int	kmIdum = 0;			// used for random number generation

//------------------------------------------------------------------------
//	kmRan0 - (safer) uniform random number generator
//
//	The code given here is taken from "Numerical Recipes in C" by
//	William Press, Brian Flannery, Saul Teukolsky, and William
//	Vetterling. The task of the code is to do an additional randomizing
//	shuffle on the system-supplied random number generator to make it
//	safer to use. 
//
//	Returns a uniform deviate between 0.0 and 1.0 using the
//	system-supplied routine "random()". Set kmIdum to any negative value
//	to initialise or reinitialise the sequence.
//------------------------------------------------------------------------

static double kmRan0()
{
    int j;

    static double y, maxran, v[98];	// The exact number 98 is unimportant
    static int iff = 0;

    // As a precaution against misuse, we will always initialize on the first
    // call, even if "kmIdum" is not set negative. Determine "maxran", the
    // next integer after the largest representable value of type int. We
    // assume this is a factor of 2 smaller than the corresponding value of
    // type unsigned int. 

    if (kmIdum < 0 || iff == 0) {	// initialize
		/* compute maximum random number */
#ifdef WIN32				// Microsoft Visual C++
	maxran = RAND_MAX;
#else
	unsigned i, k;
	i = 2;
	do {
	    k = i;
	    i <<= 1;
	} while (i);
	maxran = (double) k;
#endif
 	iff = 1;
  
	srandom(kmIdum);
	kmIdum = 1;

	for (j = 1; j <= 97; j++)	// exercise the system routine
	    random();			// (value intentionally ignored)

	for (j = 1; j <= 97; j++)	// Then save 97 values and a 98th
	    v[j] = random();
	y = random();
     }

    // This is where we start if not initializing. Use the previously saved
    // random number y to get an index j between 1 and 97. Then use the
    // corresponding v[j] for both the next j and as the output number. */

    j = 1 + (int) (97.0 * (y / maxran));
    y = v[j];
    v[j] = random();			// Finally, refill the table entry
					// with the next random number from
					// "random()" 
    return(y / maxran);
}

//------------------------------------------------------------------------
//  kmRanInt - generate a random integer from {0,1,...,n-1}
//
//	If n == 0, then -1 is returned.
//------------------------------------------------------------------------

int kmRanInt(
    int                 n)
{
    int r = (int) (kmRan0()*n);
    if (r == n) r--;			// (in case kmRan0() == 1 or n == 0)
    return r;
}

//------------------------------------------------------------------------
//  kmRanUnif - generate a random uniform in [lo,hi]
//------------------------------------------------------------------------

double kmRanUnif(
    double		lo,
    double		hi)
{
    return kmRan0()*(hi-lo) + lo;
}

//------------------------------------------------------------------------
//  kmRanGauss - Gaussian random number generator
//	Returns a normally distributed deviate with zero mean and unit
//	variance, using kmRan0() as the source of uniform deviates.
//------------------------------------------------------------------------

static double kmRanGauss()
{
    static int iset=0;
    static double gset;

    if (iset == 0) {			// we don't have a deviate handy
	double v1, v2;
	double r = 2.0;
	while (r >= 1.0) {
	    //------------------------------------------------------------
	    // Pick two uniform numbers in the square extending from -1 to
	    // +1 in each direction, see if they are in the circle of radius
	    // 1.  If not, try again 
	    //------------------------------------------------------------
	    v1 = kmRanUnif(-1, 1);
	    v2 = kmRanUnif(-1, 1);
	    r = v1 * v1 + v2 * v2;
	}
        double fac = sqrt(-2.0 * log(r) / r);
	//-----------------------------------------------------------------
	// Now make the Box-Muller transformation to get two normal
	// deviates.  Return one and save the other for next time.
	//-----------------------------------------------------------------
	gset = v1 * fac;
	iset = 1;		    	// set flag
	return v2 * fac;
    }
    else {				// we have an extra deviate handy
	iset = 0;			// so unset the flag
	return gset;			// and return it
    }
}

//------------------------------------------------------------------------
//  kmRanLaplace - Laplacian random number generator
//	Returns a Laplacian distributed deviate with zero mean and
//	unit variance, using kmRan0() as the source of uniform deviates. 
//
//		prob(x) = b/2 * exp(-b * |x|).
//
//	b is chosen to be sqrt(2.0) so that the variance of the Laplacian
//	distribution [2/(b^2)] becomes 1. 
//------------------------------------------------------------------------

static double kmRanLaplace() 
{
    const double b = 1.4142136;

    double laprand = -log(kmRan0()) / b;
    double sign = kmRan0();
    if (sign < 0.5) laprand = -laprand;
    return(laprand);
}

//----------------------------------------------------------------------
//  kmUniformPts - Generate uniformly distributed points
//	A uniform distribution over [-1,1].
//----------------------------------------------------------------------

void kmUniformPts(		// uniform distribution
	KMpointArray	pa,		// point array (modified)
	int		n,		// number of points
	int		dim)		// dimension
{
    for (int i = 0; i < n; i++) {
	for (int d = 0; d < dim; d++) {
	    pa[i][d] = (KMcoord) (kmRanUnif(-1,1));
	}
    }
}

//----------------------------------------------------------------------
//  kmGaussPts - Generate Gaussian distributed points
//	A Gaussian distribution with zero mean and the given standard
//	deviation.
//----------------------------------------------------------------------

void kmGaussPts(			// Gaussian distribution
	KMpointArray	pa,		// point array (modified)
	int		n,		// number of points
	int		dim,		// dimension
	double		std_dev)	// standard deviation
{
    for (int i = 0; i < n; i++) {
	for (int d = 0; d < dim; d++) {
	    pa[i][d] = (KMcoord) (kmRanGauss() * std_dev);
	}
    }
}

//----------------------------------------------------------------------
//  kmLaplacePts - Generate Laplacian distributed points
//	Generates a Laplacian distribution (zero mean and unit variance).
//----------------------------------------------------------------------

void kmLaplacePts(		// Laplacian distribution
	KMpointArray	pa,		// point array (modified)
	int		n,		// number of points
	int		dim)		// dimension
{
    for (int i = 0; i < n; i++) {
	for (int d = 0; d < dim; d++) {
            pa[i][d] = (KMcoord) kmRanLaplace();
	}
    }
}

//----------------------------------------------------------------------
//  kmCoGaussPts - Generate correlated Gaussian distributed points
//	Generates a Gauss-Markov distribution of zero mean and unit
//	variance.
//----------------------------------------------------------------------

void kmCoGaussPts(		// correlated-Gaussian distribution
	KMpointArray	pa,		// point array (modified)
	int		n,		// number of points
	int		dim,		// dimension
	double		correlation)	// correlation
{
    double std_dev_w = sqrt(1.0 - correlation * correlation);
    for (int i = 0; i < n; i++) {
	double previous = kmRanGauss();
	pa[i][0] = (KMcoord) previous;
	for (int d = 1; d < dim; d++) {
	    previous = correlation*previous + std_dev_w*kmRanGauss();
	    pa[i][d] = (KMcoord) previous;
	} 
    }
}

//----------------------------------------------------------------------
//  kmCoLaplacePts - Generate correlated Laplacian distributed points
//	Generates a Laplacian-Markov distribution of zero mean and unit
//	variance.
//----------------------------------------------------------------------

void kmCoLaplacePts(		// correlated-Laplacian distribution
	KMpointArray	pa,		// point array (modified)
	int		n,		// number of points
	int		dim,		// dimension
	double		correlation)	// correlation
{
    double wn;
    double corr_sq = correlation * correlation;

    for (int i = 0; i < n; i++) {
	double previous = kmRanLaplace();
	pa[i][0] = (KMcoord) previous;
	for (int d = 1; d < dim; d++) {
	    double temp = kmRan0();
	    if (temp < corr_sq)
		wn = 0.0;
	    else
		wn = kmRanLaplace();
	    previous = correlation * previous + wn;
	    pa[i][d] = (KMcoord) previous;
        } 
    }
}

//----------------------------------------------------------------------
//  kmClusGaussPts - Generate clusters of Gaussian distributed points
//	Cluster centers are uniformly distributed over [-1,1], and the
//	standard deviation within each cluster is fixed.
//
//	Note: Once cluster centers have been set, they are not changed,
//	unless new_clust = true.  This is so that subsequent calls generate
//	points from the same distribution.  It follows, of course, that any
//	attempt to change the dimension or number of clusters without
//	generating new clusters is asking for trouble.
//
//	Note: Cluster centers are not generated by a call to uniformPts().
//	Although this could be done, it has been omitted for
//	compatibility with kmClusGaussPts() in the colored version,
//	rand_c.cc.
//
//	As a side effect we compute the cluster separation. This is