kcenter.cpp

快速高斯变换的程序

#include <stdlib.h>
#include <math.h>
#include <time.h>
//#include <mkl.h>
#include "FastGauss.h"

// kcenter.cpp -  K-center algorithm, also minmax algorithm
//
// Using Gonzalez's algorithm to get a 2-approximation of the 
// optimal solution.
//
// This is a cpp file for fast Gauss transform algorithm.
// Author: Changjiang Yang.
//

// $Revision: 1.60 $
// $Date: Mon Jul  8 18:11:31 EDT 2002 $

// Add another overload function which generates the backward mapping
// from centers to points. Feb. 10, 2003

// sdist is the square of the distance of two vectors(float)
inline float FastGauss::sdist(const int d, const float *x, const float *y)
{
	float t, s = 0.0;
	for (int i = d; i != 0; i--)
	{
		t = *x++ - *y++;
		s += t * t;
	}
	
	return s;
}


// ddist is the square of the distance of two vectors(double)
inline double FastGauss::ddist(const int d, const double *x, const double *y)
{
	double t, s = 0.0;
	for (int i = d; i != 0; i--)
	{
		t = *x++ - *y++;
		s += t * t;
	}
	return s;
}



//copy vector to another vector
void FastGauss::dcopy(const int d, const double *x, double *y)
{
	for (int i = d; i != 0; i--, x++, y++)
		*y = *x;
	return;
}

//find the largest element from a vector, in case of no intel mkl library.
int FastGauss::idmax(int n, double *x)
{
	int k = 0;
	double t = -1.0;
	for (int i = 0; i < n; i++, x++)
		if( t < *x )
		{
			t = *x;
			k = i;
		}
	return k;

}


//K is the number of centers.
double FastGauss::dkcenter(const int dim, const int n, const int K,
				const double *x, double *cn, int *cind)
//dimension,number of nodes, number of centers, 
//nodes(n), centers(K), index to center(n, range:K)
{
	
	const int incu = 1; //unit increment
	double *dist_C = new double[n]; //distances to the center

	// randomly pick one node as the first center.
	srand( (unsigned)time( NULL ) );
	int ind = rand() % n;
	
	// copy the ind-th node to the first center.
#ifdef _MKL_H_
	cblas_dcopy(dim, x + ind*dim, incu, cn, incu);
#else
	dcopy(dim, x + ind*dim, cn);
#endif
	
	// compute the distances from each node to the first center.
	for (int i = 0; i < n; i++)
	{
		dist_C[i] = (i==ind)? 0.0:ddist(dim, cn, x+i*dim);
		cind[i] = ind;
	}
	
	for(int k = 1; k < K; k++)
	{	 
		//find the maximum of vector dist_C, i.e., find the node
		//that is farthest away from C
#ifdef _MKL_H_
		ind = cblas_idamax(n, dist_C, incu) - 1;
#else
		ind = idmax(n,dist_C);
#endif
		
		cn += dim; //current center.
		//copy the ind-th node to the current center.
#ifdef _MKL_H_
		cblas_dcopy(dim, x + ind*dim, incu, cn, incu);
#else
		dcopy(dim, x + ind*dim, cn);
#endif
		//update the distances from each point to the current center.
		for (int j = 0; j < n; j++)
		{
			double d = (j==ind)? 0.0:ddist(dim, cn, x+j*dim);
			if (d < dist_C[j])
			{
				dist_C[j] = d;
				cind[j] = ind;
			}
		}
	}

	//find the radius of the k-center algorithm.
#ifdef _MKL_H_
	ind = cblas_idamax(n, dist_C, incu) - 1;
#else
	ind = idmax(n,dist_C);
#endif
	
	double radius = dist_C[ind];
	
	delete []dist_C;
	
	return sqrt(radius);
}



// Overloaded function without copying the centers, instead 
// using indices to the centers. Include the numbers of nodes
// belonged in each center.
double FastGauss::dkcenter(const int dim, const int n, const int K,
				const double *x, int *cn, int *cind, int *cnn)
//dim--dimension, n--number of nodes, K--number of centers, 
//x--nodes(n), cnt--centers(K), cind--index to center(n)
//cnn--numbers of nodes(K)
{
	const int incu = 1; //unit increment
	double *dist_C = new double[n]; //distances to the center

	// randomly pick one node as the first center.
	srand( (unsigned)time( NULL ) );
	int ind = rand() % n;
	
	// add the ind-th node to the first center.
	*cn++ = ind;
	
	// compute the distances from each node to the first center.
	const double *x_ind, *x_j;
	x_ind = x + ind*dim;
	x_j = x;
	for (int j = 0; j < n; x_j += dim, j++)
	{
		dist_C[j] = (j==ind)? 0.0:ddist(dim, x_j, x_ind);
		cind[j] = 0;
	}
	
	for(int i = 1; i < K; i++)
	{	 
		//find the maximum of vector dist_C, i.e., find the node
		//that is farthest away from C
#ifdef _MKL_H_
		ind = cblas_idamax(n, dist_C, incu) - 1;
#else
		ind = idmax(n,dist_C);
#endif
		*cn++ = ind; //add the ind-th node to the current center.

		//update the distances from each point to the current center.
		x_ind = x + ind*dim;
		x_j = x;
		for (int j = 0; j < n; x_j += dim, j++)
		{
			double d = (j==ind)? 0.0:ddist(dim, x_j, x_ind);
			if (d < dist_C[j])
			{
				dist_C[j] = d;
				cind[j] = i;
			}
		}
	}

	//find the radius of the k-center algorithm.
#ifdef _MKL_H_
	ind = cblas_idamax(n, dist_C, incu) - 1;
#else
	ind = idmax(n,dist_C);
#endif

	double radius = dist_C[ind];

	delete []dist_C;

	for (int i = 0; i < K; i++)
		cnn[i] = 0;

	for (int i = 0; i < n; i++)
		cnn[cind[i]]++;
	
	return sqrt(radius);

}


// Overloaded function without copying the centers, instead 
// using indices to the centers.
double FastGauss::dkcenter(const int dim, const int n, const int K,
				const double *x, int *cn, int *cind)
//dim--dimension, n--number of nodes, K--number of centers, 
//x--nodes(n), cnt--centers(K), cind--index to center(n)
{
	
	const int incu = 1; //unit increment
	double *dist_C = new double[n]; //distances to the center

	// randomly pick one node as the first center.
	srand( (unsigned)time( NULL ) );
	int ind = rand() % n;
	
	// add the ind-th node to the first center.
	*cn++ = ind;
	
	// compute the distances from each node to the first center.
	const double *x_ind, *x_j;
	x_ind = x + ind*dim;
	x_j = x;
	for (int j = 0; j < n; x_j += dim, j++)
	{
		dist_C[j] = (j==ind)? 0.0:ddist(dim, x_j, x_ind);
		cind[j] = 0;
	}
	
	for(int i = 1; i < K; i++)
	{	 
		//find the maximum of vector dist_C, i.e., find the node
		//that is farthest away from C
#ifdef _MKL_H_
		ind = cblas_idamax(n, dist_C, incu) - 1;
#else
		ind = idmax(n,dist_C);
#endif
		*cn++ = ind; //add the ind-th node to the current center.

		//update the distances from each point to the current center.
		x_ind = x + ind*dim;
		x_j = x;
		for (int j = 0; j < n; x_j += dim, j++)
		{
			double d = (j==ind)? 0.0:ddist(dim, x_j, x_ind);
			if (d < dist_C[j])
			{
				dist_C[j] = d;
				cind[j] = i;
			}
		}
	}

	//find the radius of the k-center algorithm.
#ifdef _MKL_H_
	ind = cblas_idamax(n, dist_C, incu) - 1;
#else
	ind = idmax(n,dist_C);
#endif

	double radius = dist_C[ind];

	delete []dist_C;
	
	return sqrt(radius);
}


// Overloaded function without copying the centers, instead 
// using indices to the centers, and get the hist information
// and the indices from centers to nodes.
double FastGauss::dkcenter(const int dim, const int n, const int K,
				const double *x, int *cn, int *cind, 
				int *heads, int *boxsz, int *nind)
//dim--dimension, n--number of nodes, K--number of centers, 
//x--nodes(n), cnt--centers(K), cind--index to center(n)
//heads--(K), boxsz--(K), nind--index to point(n)
{
	double r = dkcenter(dim,n,K,x,cn,cind);

	for (int i = 0; i < K; i++)
		boxsz[i] = 0;

	for (int i = 0; i < n; i++)
		boxsz[cind[i]]++;

	
	int ct = 0;
	for (int i = 0; i < K; i++)
	{
		ct += boxsz[i];
		heads[i] = ct;
	}
	
	for (int i = n-1; i >= 0; i--)
		nind[--heads[cind[i]]] = i;

	return r;
}