apaam/mlpack
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Browse Source
This commit is contained in:
gmravi2003
2008-01-21 05:40:06 +00:00
parent 0e0e13734c
commit 80940aa0ae
17 changed files with 27 additions and 71970 deletions
Download Patch File Download Diff File Expand all files Collapse all files
fastlib/u/gmravi/regression/build.py
+3 -3
View File
@@ -1,7 +1,7 @@
binrule(
name = "main_final", # the executable name
sources = ["main_final.cc"],
headers = ["regression_new.h","regression2.h"], # no extra headers
name = "main_new", # the executable name
sources = ["main_new.cc"],
headers = ["regression_new1.h"], # no extra headers
deplibs = ["fastlib:fastlib_int"]
)
fastlib/u/gmravi/regression/colors50k.ds
-50002
View File
File diff suppressed because it is too large Load Diff
fastlib/u/gmravi/regression/fast_kde_output.txt
-25
View File
@@ -1,25 +0,0 @@
Point:1.300000 1.000000 2.000000
2.671527
3.385464
4.383671
5.647777
Point:1.100000 1.800000 1.900000
3.408173
4.289175
6.271699
7.380817
Point:2.100000 2.200000 1.500000
2.294841
3.540102
4.491305
4.387145
Point:1.000000 2.000000 3.000000
2.910201
3.260119
5.737805
7.191543
Point:1.000000 2.400000 2.600000
3.094637
3.571118
6.217940
7.430222
fastlib/u/gmravi/regression/fft_kde.h
-883
View File
@@ -1,883 +0,0 @@
#ifndef FFT_KDE_H
#define FFT_KDE_H
#include <math.h>
#include <values.h>
/**
* computing kernel estimate using Fast Fourier Transform: I have
* used multidimensional fast fourier transform called ffteasy
*/
class FFTKde {
private:
/** constant TAU */
static const double TAU = 4.0;
/** query dataset */
Matrix qset_;
/** reference dataset */
Matrix rset_;
/** kernel */
GaussianKernel kernel_;
/** computed densities */
Vector densities_;
/** number of grid points along each dimension */
int m_;
/** number of points along each dimension in the zero padded */
ArrayList<int> size_;
/** minimum coordinate along each dimension */
Vector mincoords_;
ArrayList<int> minindices_;
/** maximum coordinate along each dimension */
Vector maxcoords_;
/** difference between min and max along each dimension */
Vector diffcoords_;
/** size of grid along each dimension */
Vector gridsizes_;
/** kernel weights along each dimension */
ArrayList<int> kernelweights_dims_;
/** total number of grid points */
int numgridpts_;
/** grid box volume */
double gridbinvolume_;
/** discretized dataset storing the assigned kernel weights */
Vector discretized_;
int nyquistnum_;
Vector d_fnyquist_;
Vector k_fnyquist_;
Vector kernelweights_;
// preprocessing: scaling the dataset; this has to be moved to the dataset
// module
/* scales each attribute to 0-1 using the min/max values */
void scale_data_by_minmax() {
int num_dims = rset_.n_rows();
DHrectBound<2> qset_bound;
DHrectBound<2> rset_bound;
qset_bound.Init(qset_.n_rows());
rset_bound.Init(qset_.n_rows());
// go through each query/reference point to find out the bounds
for(index_t r = 0; r < rset_.n_cols(); r++) {
Vector ref_vector;
rset_.MakeColumnVector(r, &ref_vector);
rset_bound |= ref_vector;
}
for(index_t q = 0; q < qset_.n_cols(); q++) {
Vector query_vector;
qset_.MakeColumnVector(q, &query_vector);
qset_bound |= query_vector;
}
for(index_t i = 0; i < num_dims; i++) {
DRange qset_range = qset_bound.get(i);
DRange rset_range = rset_bound.get(i);
double min_coord = min(qset_range.lo, rset_range.lo);
double max_coord = max(qset_range.hi, rset_range.hi);
double width = max_coord - min_coord;
for(index_t j = 0; j < rset_.n_cols(); j++) {
rset_.set(i, j, (rset_.get(i, j) - min_coord) / width);
}
if(fx_param_str(NULL, "query", NULL) != NULL) {
for(index_t j = 0; j < qset_.n_cols(); j++) {
qset_.set(i, j, (qset_.get(i, j) - min_coord) / width);
}
}
}
}
/**
* Do a Fourier transform of an array of N complex numbers separated by
* steps of (complex) size skip. The array f should be of length 2N*skip
* and N must be a power of 2. Forward determines whether to do a
* forward transform (1) or an inverse one (-1)
*/
void fftc1(double *f, int N, int skip, int forward) {
int b, index1, index2, trans_size, trans;
double pi2 = 4. * asin(1.);
// used in recursive formula for Re(W^b) and Im(W^b)
double pi2n, cospi2n, sinpi2n;
// wk = W^k = e^(2 pi i b/N) in the Danielson-Lanczos formula for a
// transform of length N
struct complex wb;
// buffers for implementing recursive formulas
struct complex temp1, temp2;
// treat f as an array of N complex numbers
struct complex *c = (struct complex *)f;
// Place the elements of the array c in bit-reversed order
for(index1 = 1, index2 = 0; index1 < N; index1++) {
// to find the next bit reversed array index subtract leading 1's from
// index2
for(b = N / 2; index2 >= b; b /= 2) {
index2 -= b;
}
// Next replace the first 0 in index2 with a 1 and this gives the
// correct next value
index2 += b;
// swap each pair only the first time it is found
if(index2 > index1) {
temp1 = c[index2 * skip];
c[index2 * skip] = c[index1 * skip];
c[index1 * skip] = temp1;
}
}
// Next perform successive transforms of length 2,4,...,N using the
// Danielson-Lanczos formula
// trans_size = size of transform being computed
for(trans_size = 2; trans_size <= N; trans_size *= 2) {
// +- 2 pi/trans_size
pi2n = forward * pi2 / (double)trans_size;
// Used to calculate W^k in D-L formula
cospi2n = cos(pi2n);
sinpi2n = sin(pi2n);
// Initialize W^b for b=0
wb.real = 1.;
wb.imag = 0.;
// Step over half of the elements in the transform
for(b = 0; b < trans_size / 2; b++) {
// Iterate over all transforms of size trans_size to be computed
for(trans = 0; trans < N / trans_size; trans++) {
// Index of element in first half of transform being computed
index1 = (trans * trans_size + b) * skip;
// Index of element in second half of transform being computed
index2 = index1 + trans_size / 2 * skip;
temp1 = c[index1];
temp2 = c[index2];
// implement D-L formula
c[index1].real = temp1.real + wb.real * temp2.real -
wb.imag * temp2.imag;
c[index1].imag = temp1.imag + wb.real * temp2.imag +
wb.imag * temp2.real;
c[index2].real = temp1.real - wb.real * temp2.real +
wb.imag * temp2.imag;
c[index2].imag = temp1.imag - wb.real * temp2.imag -
wb.imag * temp2.real;
}
temp1 = wb;
// Real part of e^(2 pi i b/trans_size) used in D-L formula
wb.real = cospi2n * temp1.real - sinpi2n * temp1.imag;
// Imaginary part of e^(2 pi i b/trans_size) used in D-L formula
wb.imag = cospi2n*temp1.imag + sinpi2n*temp1.real;
}
}
// For an inverse transform divide by the number of grid points
if(forward<0) {
for(index1 = 0; index1 < skip * N; index1 += skip) {
c[index1].real /= N;
c[index1].imag /= N;
}
}
}
/**
* Do a Fourier transform of an ndims dimensional array of complex numbers
* Array dimensions are given by size[0],...,size[ndims-1]. Note that these
* are sizes of complex arrays. The array f should be of length
* 2*size[0]*...*size[ndims-1] and all sizes must be powers of 2.
* Forward determines whether to do a forward transform (1) or an inverse
* one(-1)
*/
void fftcn(double *f, int ndims, int *size, int forward) {
// These determine where to begin successive transforms and the skip
// between their elements (see below)
int planesize = 1, skip = 1;
// Total size of the ndims dimensional array
int totalsize = 1;
// determine total size of array
for(index_t dim = 0; dim < ndims; dim++) {
totalsize *= size[dim];
}
// loop over dimensions
for(index_t dim = ndims - 1; dim >= 0; dim--) {
// planesize = Product of all sizes up to and including size[dim]
planesize *= size[dim];
// Take big steps to begin loops of transforms
for(index_t i = 0; i < totalsize; i += planesize) {
// Skip sets the number of transforms in between big steps as well as
// the skip between elements
for(index_t j = 0; j < skip; j++) {
// 1-D Fourier transform. (Factor of two converts complex index to
// double index.)
fftc1(f + 2 * (i + j), size[dim], skip, forward);
}
}
// Skip = Product of all sizes up to (but not including) size[dim]
skip *= size[dim];
}
}
/**
* Do a Fourier transform of an array of N real numbers
* N must be a power of 2
* Forward determines whether to do a forward transform (>=0) or an inverse
* one(<0)
*/
void fftr1(double *f, int N, int forward) {
int b;
// pi2n = 2 Pi/N
double pi2n = 4. * asin(1.) / N, cospi2n = cos(pi2n), sinpi2n = sin(pi2n);
// wb = W^b = e^(2 pi i b/N) in the Danielson-Lanczos formula for a
// transform of length N
struct complex wb;
// Buffers for implementing recursive formulas
struct complex temp1, temp2;
// Treat f as an array of N/2 complex numbers
struct complex *c = (struct complex *)f;
// Do a transform of f as if it were N/2 complex points
if(forward == 1) {
fftc1(f, N / 2, 1, 1);
}
// initialize W^b for b = 0
wb.real = 1.;
wb.imag = 0.;
// Loop over elements of transform. See documentation for these formulas
for(b = 1; b < N / 4; b++) {
temp1 = wb;
// Real part of e^(2 pi i b/N) used in D-L formula
wb.real = cospi2n * temp1.real - sinpi2n * temp1.imag;
// Imaginary part of e^(2 pi i b/N) used in D-L formula
wb.imag = cospi2n * temp1.imag + sinpi2n * temp1.real;
temp1 = c[b];
temp2 = c[N / 2 - b];
c[b].real = .5 * (temp1.real + temp2.real + forward * wb.real *
(temp1.imag + temp2.imag) + wb.imag *
(temp1.real - temp2.real));
c[b].imag = .5 * (temp1.imag-temp2.imag - forward * wb.real *
(temp1.real - temp2.real) + wb.imag *
(temp1.imag + temp2.imag));
c[N/2-b].real = .5 * (temp1.real + temp2.real - forward * wb.real *
(temp1.imag + temp2.imag) - wb.imag *
(temp1.real - temp2.real));
c[N/2-b].imag = .5 * (-temp1.imag + temp2.imag - forward * wb.real *
(temp1.real - temp2.real) + wb.imag *
(temp1.imag + temp2.imag));
}
temp1 = c[0];
// set b = 0 term in transform
c[0].real = temp1.real+temp1.imag;
// put b = N / 2 term in imaginary part of first term
c[0].imag = temp1.real-temp1.imag;
if(forward == -1) {
c[0].real *= .5;
c[0].imag *= .5;
fftc1(f, N / 2, 1, -1);
}
}
/**
* Do a Fourier transform of an ndims dimensional array of real numbers
* Array dimensions are given by size[0],...,size[ndims-1]. All sizes must
* be powers of 2. The (complex) nyquist frequency components are stored in
* fnyquist[size[0]][size[1]]...[2*size[ndims-2]]
* Forward determines whether to do a forward transform (1) or an inverse
* one (-1)
*/
void fftrn(double *f, double *fnyquist, int ndims, int *size, int forward) {
int i, j, b;
// Positions in the 1-d arrays of points labeled by indices
// (i0,i1,...,i(ndims-1)); indexneg gives the position in the array of
// the corresponding negative frequency
int index,indexneg = 0;
int stepsize; // Used in calculating indexneg
// The size of the last dimension is used often enough to merit its own
// name.
int N = size[ndims - 1];
// pi2n = 2Pi / N
double pi2n = 4. * asin(1.) / N, cospi2n = cos(pi2n), sinpi2n = sin(pi2n);
// wb = W^b = e^(2 pi i b/N) in the Danielson-Lanczos formula for a
// transform of length N
struct complex wb;
// Buffers for implementing recursive formulas
struct complex temp1, temp2;
// Treat f and fnyquist as arrays of complex numbers
struct complex *c = (struct complex *)f,
*cnyquist = (struct complex *)fnyquist;
// Total number of complex points in array
int totalsize = 1;
// Indices for looping through array
ArrayList<int> indices;
indices.Init(ndims);
// Set size[] to be the sizes of f viewed as a complex array
size[ndims - 1] /= 2;
for(i = 0; i < ndims; i++) {
totalsize *= size[i];
indices[i] = 0;
}
// forward transform
if(forward == 1) {
// Do a transform of f as if it were N/2 complex points
fftcn(f, ndims, size, 1);
// Copy b=0 data into cnyquist so the recursion formulas below for b=0
// and cnyquist don't overwrite data they later need
for(i = 0; i < totalsize / size[ndims - 1]; i++) {
// Only copy points where last array index for c is 0
cnyquist[i] = c[i * size[ndims - 1]];
}
}
// Loop over all but last array index
for(index = 0; index < totalsize; index += size[ndims-1]) {
wb.real = 1.; /* Initialize W^b for b=0 */
wb.imag = 0.;
// Loop over elements of transform. See documentation for these formulas
for(b = 1; b < N / 4; b++) {
temp1 = wb;
// Real part of e^(2 pi i b/N_real) used in D-L formula
wb.real = cospi2n*temp1.real - sinpi2n*temp1.imag;
// Imaginary part of e^(2 pi i b/N_real) used in D-L formula
wb.imag = cospi2n*temp1.imag + sinpi2n*temp1.real;
temp1 = c[index + b];
// Note that N-b is NOT the negative frequency for b. Only
// nonnegative b momenta are stored.
temp2 = c[indexneg + N / 2 - b];
c[index + b].real = .5 * (temp1.real + temp2.real + forward * wb.real *
(temp1.imag + temp2.imag) + wb.imag *
(temp1.real - temp2.real));
c[index + b].imag = .5 * (temp1.imag - temp2.imag - forward * wb.real *
(temp1.real - temp2.real) + wb.imag *
(temp1.imag + temp2.imag));
c[indexneg + N / 2 - b].real = .5 * (temp1.real + temp2.real -
forward *
wb.real *
(temp1.imag + temp2.imag) -
wb.imag *
(temp1.real - temp2.real));
c[indexneg + N / 2 - b].imag = .5 * (-temp1.imag + temp2.imag -
forward * wb.real *
(temp1.real - temp2.real) +
wb.imag *
(temp1.imag + temp2.imag));
}
temp1 = c[index];
// Index is smaller for cnyquist because it doesn't have the last
// dimension
temp2 = cnyquist[indexneg / size[ndims - 1]];
// Set b=0 term in transform
c[index].real = .5 * (temp1.real + temp2.real + forward *
(temp1.imag + temp2.imag));
c[index].imag = .5 * (temp1.imag - temp2.imag - forward *
(temp1.real - temp2.real));
// Set b=N/2 transform.
cnyquist[indexneg / size[ndims - 1]].real =
.5 * (temp1.real + temp2.real - forward * (temp1.imag + temp2.imag));
cnyquist[indexneg / size[ndims - 1]].imag =
.5 * (-temp1.imag + temp2.imag - forward * (temp1.real - temp2.real));
// Find indices for positive and single index for negative frequency.
// In each dimension indexneg[j]=0 if index[j]=0,
// indexneg[j]=size[j]-index[j] otherwise.
// amount to increment indexneg by as each individual index is
// incremented
stepsize = size[ndims - 1];
// If the rightmost indices are maximal reset them to 0. Indexneg goes
// from 1 to 0 in these dimensions
for(j = ndims - 2; j >= 0 && indices[j] == size[j] - 1; j--) {
indices[j] = 0;
indexneg -= stepsize;
stepsize *= size[j];
}
// If index[j] goes from 0 to 1 indexneg[j] goes from 0 to size[j]-1
if(j >= 0 && indices[j] == 0) {
indexneg += stepsize * (size[j] - 1);
}
// Otherwise increasing index[j] decreases indexneg by one unit.
else {
indexneg -= stepsize;
}
// This avoids writing outside the array bounds on the last pass
// through the array loop
if(j >= 0) {
indices[j]++;
}
} // End of i loop (over total array)
// inverse transform
if(forward == -1) {
fftcn(f, ndims, size, -1);
}
// Give the user back the array size[] in its original condition
size[ndims - 1] *= 2;
}
void assign_weights(int reference_pt_num, int level, double volume, int pos,
int skip) {
if(level == -1) {
discretized_[pos] += volume;
}
else {
// Recurse in the right direction
double coord = rset_.get(level, reference_pt_num);
double leftgridcoord = mincoords_[level] + minindices_[level] *
gridsizes_[level];
double rightgridcoord = leftgridcoord + gridsizes_[level];
double leftvolume = volume * (rightgridcoord - coord);
double rightvolume = volume * (coord - leftgridcoord);
int nextskip = size_[level] * skip;
int nextleftpos = pos + skip * minindices_[level];
if(leftvolume > 0.0) {
assign_weights(reference_pt_num, level - 1, leftvolume, nextleftpos,
nextskip);
}
if(rightvolume > 0.0) {
assign_weights(reference_pt_num, level - 1, rightvolume,
nextleftpos + skip, nextskip);
}
}
}
void retrieve_weights(int query_pt_num, double volume, int level, int pos,
int skip, double divfactor) {
if(level == -1) {
densities_[query_pt_num] += discretized_[pos] * volume / divfactor;
}
else {
// Recurse in the right direction
double coord = qset_.get(level, query_pt_num);
double leftgridcoord = mincoords_[level] + minindices_[level] *
gridsizes_[level];
double rightgridcoord = leftgridcoord + gridsizes_[level];
double leftvolume = volume * (rightgridcoord - coord);
double rightvolume = volume * (coord - leftgridcoord);
int nextskip = size_[level] * skip;
int nextleftpos = pos + skip * minindices_[level];
if(leftvolume > 0.0) {
retrieve_weights(query_pt_num, leftvolume, level - 1, nextleftpos,
nextskip, divfactor);
}
if(rightvolume > 0.0) {
retrieve_weights(query_pt_num, rightvolume, level - 1,
nextleftpos + skip, nextskip, divfactor);
}
}
}
/**
* Query the normalized density for each query point.
*/
void RetrieveDensities() {
double normc =
(kernel_.CalcNormConstant(rset_.n_rows()) * rset_.n_cols());
for(index_t r = 0; r < qset_.n_cols(); r++) {
densities_[r] = 0.0;
for(index_t d = 0; d < qset_.n_rows(); d++) {
minindices_[d] = (int) floor((qset_.get(d, r) - mincoords_[d])/
gridsizes_[d]);
}
retrieve_weights(r, 1.0, qset_.n_rows() - 1, 0, 1,
gridbinvolume_ * normc);
}
}
void discretize_dataset() {
// Temporary used to count the number of elements in the enlarged
// matrices for the kernel weights and bin counts. Also calculate the
// volume of each grid bin.
numgridpts_ = 1;
gridbinvolume_ = 1.0;
double min, max;
// Find the min/max in each coordinate direction, and calculate the grid
// size in each dimension.
for(index_t d = 0; d < qset_.n_rows(); d++) {
int possiblesample;
min = MAXDOUBLE;
max = -MAXDOUBLE;
for(index_t r = 0; r < rset_.n_cols(); r++) {
double coord = rset_.get(d, r);
if(coord > max)
max = coord;
if(coord < min)
min = coord;
}
// Following Silverman's advice here
mincoords_[d] = min;
maxcoords_[d] = max;
diffcoords_[d] = maxcoords_[d] - mincoords_[d];
gridsizes_[d] = diffcoords_[d] / ((double) m_ - 1);
gridbinvolume_ *= gridsizes_[d];
// Determine how many kernel weight calculation to do for this
// dimension.
kernelweights_dims_[d] = m_ - 1;
possiblesample = (int) floor(TAU * sqrt(kernel_.bandwidth_sq()) /
gridsizes_[d]);
if(kernelweights_dims_[d] > possiblesample) {
if(possiblesample == 0) {
possiblesample = 1;
}
kernelweights_dims_[d] = possiblesample;
}
// Wand p440: Need to calculate the actual dimension of the matrix
// after the necessary 0 padding of the kernel weight matrix and the
// bin count matrix.
size_[d] = (int) ceil(log(m_ + kernelweights_dims_[d]) / log(2));
size_[d] = 1 << size_[d];
numgridpts_ *= size_[d];
}
// Allocate the memory for discretized grid count matrix and initialize
// it.
discretized_.Init(numgridpts_);
discretized_.SetZero();
double inv_gvolume = 1.0 / gridbinvolume_;
// Now loop over each data and calculate the weights at each grid point.
for(index_t r = 0; r < rset_.n_cols(); r++) {
// First locate the bin the data point falls into and identify it by
// the lower grid coordinates.
for(index_t d = 0; d < rset_.n_rows(); d++) {
minindices_[d] = (int) floor((rset_.get(d, r) - mincoords_[d])/
gridsizes_[d]);
}
// Assign the weights around the neighboring grid points due to this
// data point. This results in 2^num_dims number of recursion per data
// point.
assign_weights(r, qset_.n_rows() - 1, inv_gvolume, 0, 1);
}
}
void gaussify(double acc, double precalc, int level, int pos, int skip) {
if(level == -1) {
kernelweights_[pos] = exp(precalc * acc);
}
else {
int half = kernelweights_dims_[level];
int g;
for(g = 0; g <= half; g++) {
double addThis = g * gridsizes_[level];
double newacc = acc + addThis * addThis;
int newskip = skip * size_[level];
gaussify(newacc, precalc, level - 1, pos + skip * g, newskip);
// If this is not the 0th frequency, then do the mirror image thingie.
if(g != 0) {
gaussify(newacc, precalc, level - 1,
pos + skip * (size_[level] - g), newskip);
}
}
}
}
public:
struct complex {
double real;
double imag;
};
FFTKde() {}
~FFTKde() {}
// getters and setters
/** get the reference dataset */
Matrix &get_reference_dataset() { return rset_; }
/** get the query dataset */
Matrix &get_query_dataset() { return qset_; }
/** get the density estimate */
const Vector &get_density_estimates() { return densities_; }
void Init(Matrix &qset, Matrix &rset) {
printf("Initializing FFT KDE...\n");
fx_timer_start(NULL, "fft_kde_init");
// initialize the kernel and read in the number of grid points
kernel_.Init(fx_param_double_req(NULL, "bandwidth"));
m_ = fx_param_int(NULL, "num_grid_pts_per_dim", 128);
// set aliases to the query and reference datasets and initialize
// query density sets
qset_.Alias(qset);
densities_.Init(qset_.n_cols());
rset_.Alias(rset);
// scale dataset if the user wants to
if(!strcmp(fx_param_str(NULL, "scaling", NULL), "range")) {
scale_data_by_minmax();
}
// initialize member variables.
size_.Init(qset_.n_rows());
minindices_.Init(rset_.n_rows());
mincoords_.Init(qset_.n_rows());
maxcoords_.Init(qset_.n_rows());
diffcoords_.Init(qset_.n_rows());
gridsizes_.Init(qset_.n_rows());
kernelweights_dims_.Init(qset_.n_rows());
// set up the discretized grid for the reference dataset
discretize_dataset();
nyquistnum_ = 2 * numgridpts_ / size_[rset_.n_rows() - 1];
d_fnyquist_.Init(nyquistnum_);
k_fnyquist_.Init(nyquistnum_);
kernelweights_.Init(numgridpts_);
fx_timer_stop(NULL, "fft_kde_init");
printf("FFT KDE initialization completed...\n");
}
void Init() {
const char *rfname = fx_param_str_req(NULL, "data");
const char *qfname = fx_param_str(NULL, "query", rfname);
// initialize the kernel and read in the number of grid points
kernel_.Init(fx_param_double_req(NULL, "bandwidth"));
m_ = fx_param_int(NULL, "num_grid_pts_per_dim", 128);
// read reference dataset
Dataset ref_dataset;
ref_dataset.InitFromFile(rfname);
rset_.Own(&(ref_dataset.matrix()));
// read query dataset if different
if(!strcmp(qfname, rfname)) {
qset_.Alias(rset_);
}
else {
Dataset query_dataset;
query_dataset.InitFromFile(qfname);
qset_.Own(&(query_dataset.matrix()));
}
// scale dataset if the user wants to
if(!strcmp(fx_param_str(NULL, "scaling", NULL), "range")) {
scale_data_by_minmax();
}
printf("Initializing FFT KDE...\n");
fx_timer_start(NULL, "fft_kde_init");
// initialize member variables.
size_.Init(qset_.n_rows());
densities_.Init(qset_.n_cols());
minindices_.Init(rset_.n_rows());
mincoords_.Init(qset_.n_rows());
maxcoords_.Init(qset_.n_rows());
diffcoords_.Init(qset_.n_rows());
gridsizes_.Init(qset_.n_rows());
kernelweights_dims_.Init(qset_.n_rows());
// set up the discretized grid for the reference dataset
discretize_dataset();
nyquistnum_ = 2 * numgridpts_ / size_[rset_.n_rows() - 1];
d_fnyquist_.Init(nyquistnum_);
k_fnyquist_.Init(nyquistnum_);
kernelweights_.Init(numgridpts_);
fx_timer_stop(NULL, "fft_kde_init");
printf("FFT KDE initialization completed...\n");
}
void Compute() {
printf("Computing FFT KDE...\n");
fx_timer_start(NULL, "fft_kde");
// FFT the discretized bin count matrix.
d_fnyquist_.SetZero();
k_fnyquist_.SetZero();
kernelweights_.SetZero();
fftrn(discretized_.ptr(), d_fnyquist_.ptr(), rset_.n_rows(),
size_.begin(), 1);
// Calculate the required kernel weights at each grid point. This matrix
// will be convolved with fourier transformed data set.
double precalc = -0.5 / kernel_.bandwidth_sq();
gaussify(0.0, precalc, rset_.n_rows() - 1, 0, 1);
// FFT the kernel weight matrix.
fftrn(kernelweights_.ptr(), k_fnyquist_.ptr(),
rset_.n_rows(), size_.begin(), 1);
// We need to invoke the convolution theorem for FFT here. Take each
// corresponding complex number in kernelweights and discretized and do
// an element-wise multiplication. Later, pass it to inverse fft function,
// and we have our answer!
for(index_t d = 0; d < numgridpts_; d += 2) {
double real1 = discretized_[d];
double complex1 = discretized_[d + 1];
double real2 = kernelweights_[d];
double complex2 = kernelweights_[d + 1];
discretized_[d] = real1 * real2 - complex1 * complex2;
discretized_[d + 1] = real1 * complex2 + complex1 * real2;
}
for(index_t d = 0; d < nyquistnum_; d += 2) {
double real1 = d_fnyquist_[d];
double complex1 = d_fnyquist_[d + 1];
double real2 = k_fnyquist_[d];
double complex2 = k_fnyquist_[d + 1];
d_fnyquist_[d] = real1 * real2 - complex1 * complex2;
d_fnyquist_[d + 1] = real1 * complex2 + complex1 * real2;
}
// Inverse FFT the elementwise multiplied matrix.
fftrn(discretized_.ptr(), d_fnyquist_.ptr(),
rset_.n_rows(), size_.begin(), -1);
// Retrieve the densities of each data point.
RetrieveDensities();
fx_timer_stop(NULL, "fft_kde");
printf("FFT KDE completed...\n");
}
void NormalizeDensities() {
double norm_const = kernel_.CalcNormConstant(qset_.n_rows()) *
rset_.n_cols();
for(index_t q = 0; q < qset_.n_cols(); q++) {
densities_[q] /= norm_const;
}
}
void PrintDebug() {
FILE *stream = stdout;
const char *fname = NULL;
if((fname = fx_param_str(NULL, "fft_kde_output", NULL)) != NULL) {
stream = fopen(fname, "w+");
}
for(index_t q = 0; q < qset_.n_cols(); q++) {
fprintf(stream, "%g\n", densities_[q]);
}
if(stream != stdout) {
fclose(stream);
}
}
};
#endif
fastlib/u/gmravi/regression/fgt_kde.h
-1040
View File
File diff suppressed because it is too large Load Diff
fastlib/u/gmravi/regression/ifgt_kde.cc
-346
View File
@@ -1,346 +0,0 @@
//-------------------------------------------------------------------
// The code was written by Vikas Raykar and Changjiang Yang
// and is copyrighted under the Lesser GPL:
//
// Copyright (C) 2006 Vikas Raykar and Changjiang Yang
//
// This program is free software; you can redistribute it and/or modify
// it under the terms of the GNU Lesser General Public License as
// published by the Free Software Foundation; version 2.1 or later.
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
// See the GNU Lesser General Public License for more details.
// You should have received a copy of the GNU Lesser General Public
// License along with this program; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place - Suite 330, Boston,
// MA 02111-1307, USA.
//
// The author may be contacted via email at:
// vikas(at)umiacs(.)umd(.)edu, cyang(at)sarnoff(.)com
//-------------------------------------------------------------------
//-------------------------------------------------------------------
// File : ImprovedFastGaussTransform.cpp
// Purpose : Implementation for the Improved Fast Gauss Transform
// Author : Vikas C. Raykar (vikas@cs.umd.edu)
// Date : July 15 2005
//-------------------------------------------------------------------
#include "ifgt_kde.h"
#include "fastlib/fastlib_int.h"
#include <math.h>
#include <values.h>
//-------------------------------------------------------------------
// Constructor
//
// PURPOSE
// -------
// Initialize the class.
// Read the parameters.
// Allocate memory.
//-------------------------------------------------------------------
ImprovedFastGaussTransform::ImprovedFastGaussTransform(int Dim,
int NSources,
int MTargets,
double *pSources,
double Bandwidth,
double *pWeights,
double *pTargets,
int MaxTruncNumber,
int NumClusters,
int *pClusterIndex,
double *pClusterCenter,
double *pClusterRadii,
double CutoffRadius,
double epsilon,
double *pGaussTransform,
int *pTruncNumber
)
{
//Read the parameters
d=Dim;
N=NSources;
M=MTargets;
px=pSources;
h=Bandwidth * sqrt(2);
pq=pWeights;
py=pTargets;
p_max=MaxTruncNumber;
K=NumClusters;
pci=pClusterIndex;
pcc=pClusterCenter;
pcr=pClusterRadii;
r=CutoffRadius;
pG=pGaussTransform;
pT=pTruncNumber;
eps=epsilon;
//Memory allocation
p_max_total=nchoosek(p_max-1+d,d);
constant_series=new double[p_max_total];
source_center_monomials = new double[p_max_total];
target_center_monomials = new double[p_max_total];
dx = new double[d];
dy = new double[d];
heads = new int[d];
C=new double[K*p_max_total];
h_square=h*h;
ry=new double[K];
ry_square=new double[K];
for(int i=0; i<K; i++)
{
ry[i]=r+pcr[i];
ry_square[i]=ry[i]*ry[i];
}
}
//-------------------------------------------------------------------
// Destructor
//-------------------------------------------------------------------
ImprovedFastGaussTransform::~ImprovedFastGaussTransform()
{
}
//-------------------------------------------------------------------
// Compute the combinatorial number nchoosek.
//-------------------------------------------------------------------
int
ImprovedFastGaussTransform::nchoosek(int n, int k){
int n_k = n - k;
if (k < n_k)
{
k = n_k;
n_k = n - k;
}
int nchsk = 1;
for ( int i = 1; i <= n_k; i++)
{
nchsk *= (++k);
nchsk /= i;
}
return nchsk;
}
//-------------------------------------------------------------------
//Computes p_i such error(a,p_i,h) <= q_i epsilon.
//-------------------------------------------------------------------
int
ImprovedFastGaussTransform::return_p(double a_square, int cluster_index)
{
double a=sqrt(a_square);
double b,c;
double error=1;
double temp=1;
int p=1;
while((error > eps) & (p <= p_max))
{
b=min(((a+sqrt((a_square)+(2*p*h_square)))/2),ry[cluster_index]);
c=a-b;
temp=temp*(((2*a*b)/h_square)/p);
error=temp*(exp(-(c*c)/h_square));
p++;
}
return p-1;
}
//-------------------------------------------------------------------
// This function computes the constants 2^alpha/alpha!.
//-------------------------------------------------------------------
void
ImprovedFastGaussTransform::compute_constant_series(){
int *heads = new int[d+1];
int *cinds = new int[p_max_total];
for (int i = 0; i < d; i++)
heads[i] = 0;
heads[d] = MAXINT;
cinds[0] = 0;
constant_series[0] = 1.0;
for (int k=1, t=1, tail=1; k < p_max; k++, tail=t)
{
for (int i = 0; i < d; i++)
{
int head = heads[i];
heads[i] = t;
for ( int j = head; j < tail; j++, t++)
{
cinds[t] = (j < heads[i+1])? cinds[j] + 1 : 1;
constant_series[t] = 2.0 * constant_series[j];
constant_series[t] /= (double) cinds[t];
}
}
}
delete []cinds;
delete []heads;
}
//-------------------------------------------------------------------
// This function computes the monomials [(x_i-c_k)/h]^{alpha}
// and norm([(x_i-c_k)/h])^2
//-------------------------------------------------------------------
void
ImprovedFastGaussTransform::compute_source_center_monomials(int p)
{
for (int i = 0; i < d; i++){
dx[i]=dx[i]/h;
heads[i] = 0;
}
source_center_monomials[0] = 1.0;
for (int k=1, t=1, tail=1; k < p; k++, tail=t){
for (int i = 0; i < d; i++){
int head = heads[i];
heads[i] = t;
for ( int j = head; j < tail; j++, t++)
source_center_monomials[t] = dx[i] * source_center_monomials[j];
}
}
}
//-------------------------------------------------------------------
// This function computes the monomials [(y_j-c_k)/h]^{alpha}
//-------------------------------------------------------------------
void
ImprovedFastGaussTransform::compute_target_center_monomials()
{
for (int i = 0; i < d; i++){
dy[i]=dy[i]/h;
heads[i] = 0;
}
target_center_monomials[0] = 1.0;
for (int k=1, t=1, tail=1; k < p_max; k++, tail=t){
for (int i = 0; i < d; i++){
int head = heads[i];
heads[i] = t;
for ( int j = head; j < tail; j++, t++)
target_center_monomials[t] = dy[i] * target_center_monomials[j];
}
}
}
//-------------------------------------------------------------------
// This function computes the coeffeicients C_k for all clusters.
//-------------------------------------------------------------------
void
ImprovedFastGaussTransform::compute_C()
{
for (int i = 0; i < K*p_max_total; i++){
C[i]=0.0;
}
p_max_actual=-1;
for(int i=0; i<N; i++){
int k=pci[i];
int source_base=i*d;
int center_base=k*d;
source_center_distance_square=0.0;
for (int j = 0; j < d; j++){
dx[j]=(px[source_base+j]-pcc[center_base+j]);
source_center_distance_square += (dx[j]*dx[j]);
}
pT[i]=return_p(source_center_distance_square,k);
if (pT[i]>p_max_actual){
p_max_actual=pT[i];
}
compute_source_center_monomials(pT[i]);
double f=pq[i]*exp(-source_center_distance_square/h_square);
for(int alpha=0; alpha<nchoosek(pT[i]-1+d,d); alpha++){
C[k*p_max_total+alpha]+=(f*source_center_monomials[alpha]);
}
}
p_max_actual_total=nchoosek(p_max_actual-1+d,d);
compute_constant_series();
for(int k=0; k<K; k++){
for(int alpha=0; alpha<p_max_total; alpha++){
C[k*p_max_total+alpha]*=constant_series[alpha];
}
}
}
//-------------------------------------------------------------------
// Actual function to evaluate the Gauss Transform.
//-------------------------------------------------------------------
void
ImprovedFastGaussTransform::Evaluate()
{
compute_C();
for(int j=0; j < M; j++)
{
pG[j]=0.0;
int target_base=j*d;
for(int k=0; k<K; k++){
int center_base=k*d;
double target_center_distance_square=0.0;
for(int i=0; i<d; i++){
dy[i]=py[target_base+i]-pcc[center_base+i];
target_center_distance_square += dy[i]*dy[i];
if (target_center_distance_square > ry_square[k]) break;
}
if (target_center_distance_square <= ry_square[k]){
compute_target_center_monomials();
double g=exp(-target_center_distance_square/h_square);
for(int alpha=0; alpha<p_max_actual_total; alpha++){
pG[j]+=(C[k*p_max_total+alpha]*g*target_center_monomials[alpha]);
}
}
}
}
}
fastlib/u/gmravi/regression/ifgt_kde.h
-139
View File
@@ -1,139 +0,0 @@
//-------------------------------------------------------------------
// The code was written by Vikas Raykar and Changjiang Yang
// and is copyrighted under the Lesser GPL:
//
// Copyright (C) 2006 Vikas Raykar and Changjiang Yang
//
// This program is free software; you can redistribute it and/or modify
// it under the terms of the GNU Lesser General Public License as
// published by the Free Software Foundation; version 2.1 or later.
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
// See the GNU Lesser General Public License for more details.
// You should have received a copy of the GNU Lesser General Public
// License along with this program; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place - Suite 330, Boston,
// MA 02111-1307, USA.
//
// The author may be contacted via email at:
// vikas(at)umiacs(.)umd(.)edu, cyang(at)sarnoff(.)com
//-------------------------------------------------------------------
//-------------------------------------------------------------
// File : ImprovedFastGaussTransform.h
// Purpose : Interface for
// Data Adaptive Improved Fast Gauss Transform.
// Author : Vikas C. Raykar (vikas@cs.umd.edu)
// Date : July 15 2005
//-------------------------------------------------------------
// Data Adaptive Improved Fast Gauss Transform (IFGT).
// All Sources have the same scales 'h'.
//
//
// A new version of the IFGT where the parameters are chosen
// based on the acutal distribution of the source points.
// The truncation number for each source point is chosen based
// on its distance to the cluster center.
//
// Advantages:
// -----------
// 1. Better Speedup.
// 2. Choice of parameters is fully automatic taking into
// consideration the actual distribtion of the data points.
// 3. Uses more tight pointwise error bounds in choosing the
// parameters.
//
// Implementation based on:
//
// Fast computation of sums of Gaussians in high dimensions.
// Vikas C. Raykar, C. Yang, R. Duraiswami, and N. Gumerov,
// CS-TR-4767, Department of computer science,
// University of Maryland, Collegepark.
// ------------------------------------------------------------
#ifndef IMPROVED_FAST_GAUSS_TRANSFORM_H
#define IMPROVED_FAST_GAUSS_TRANSFORM_H
class ImprovedFastGaussTransform{
public:
//constructor
ImprovedFastGaussTransform(int Dim,
int NSources,
int MTargets,
double *pSources,
double Bandwidth,
double *pWeights,
double *pTargets,
int MaxTruncNumber,
int NumClusters,
int *pClusterIndex,
double *pClusterCenter,
double *pClusterRadii,
double CutoffRadius,
double epsilon,
double *pGaussTransform,
int *pTruncNumber
);
//destructor
~ImprovedFastGaussTransform();
//function to evaluate the Gauss Transform.
void Evaluate();
private:
//Parameters
int d;
int N;
int M;
double *px;
double h;
double *pq;
double *py;
int p_max;
int K;
int *pci;
double *pcc;
double *pcr;
double r;
double eps;
double *pG;
int *pT;
//
int p_max_total;
int p_max_actual;
int p_max_actual_total;
double *constant_series;
double *source_center_monomials;
double source_center_distance_square;
double *target_center_monomials;
double target_center_distance_square;
double *dx;
double *dy;
int *heads;
double *C;
double h_square;
double *ry;
double *ry_square;
//Functions
int nchoosek(int n, int k);
int return_p(double a_square, int cluster_index);
void compute_constant_series();
void compute_source_center_monomials(int p);
void compute_target_center_monomials();
void compute_C();
};
#endif
fastlib/u/gmravi/regression/ifgt_main.cc
-193
View File
@@ -1,193 +0,0 @@
#include "fastlib/fastlib_int.h"
#include "ifgt_kde.h"
#include "ifgt_choose_parameters.h"
#include "kde.h"
#include "ifgt_choose_truncation_number.h"
#include "kcenter_clustering.h"
void concatenate_vectors(Matrix &source, Vector &dest) {
for(index_t i = 0; i < source.n_cols(); i++) {
for(index_t j = 0; j < source.n_rows(); j++) {
dest[i * source.n_rows() + j] = source.get(j, i);
}
}
}
// preprocessing: scaling the dataset; this has to be moved to the dataset
// module
/* scales each attribute to 0-1 using the min/max values */
void scale_data_by_minmax(Matrix &qset_, Matrix &rset_) {
int num_dims = rset_.n_rows();
DHrectBound<2> qset_bound;
DHrectBound<2> rset_bound;
qset_bound.Init(qset_.n_rows());
rset_bound.Init(qset_.n_rows());
// go through each query/reference point to find out the bounds
for(index_t r = 0; r < rset_.n_cols(); r++) {
Vector ref_vector;
rset_.MakeColumnVector(r, &ref_vector);
rset_bound |= ref_vector;
}
for(index_t q = 0; q < qset_.n_cols(); q++) {
Vector query_vector;
qset_.MakeColumnVector(q, &query_vector);
qset_bound |= query_vector;
}
for(index_t i = 0; i < num_dims; i++) {
DRange qset_range = qset_bound.get(i);
DRange rset_range = rset_bound.get(i);
double min_coord = min(qset_range.lo, rset_range.lo);
double max_coord = max(qset_range.hi, rset_range.hi);
double width = max_coord - min_coord;
for(index_t j = 0; j < rset_.n_cols(); j++) {
rset_.set(i, j, (rset_.get(i, j) - min_coord) / width);
}
if(fx_param_str(NULL, "query", NULL) != NULL) {
for(index_t j = 0; j < qset_.n_cols(); j++) {
qset_.set(i, j, (qset_.get(i, j) - min_coord) / width);
}
}
}
}
int main(int argc, char *argv[]) {
fx_init(argc, argv);
// read the datasets and do k-center clustering
Dataset ref_dataset;
Matrix qset_;
Matrix rset_;
Vector pWeights;
// read the datasets
const char *rfname = fx_param_str_req(NULL, "data");
const char *qfname = fx_param_str(NULL, "query", rfname);
// read reference dataset
ref_dataset.InitFromFile(rfname);
rset_.Own(&(ref_dataset.matrix()));
// read the reference weights
char *rwfname = NULL;
if(fx_param_exists(NULL, "dwgts")) {
rwfname = (char *)fx_param_str(NULL, "dwgts", NULL);
}
if(rwfname != NULL) {
Dataset ref_weights;
ref_weights.InitFromFile(rwfname);
pWeights.Copy(ref_weights.matrix().GetColumnPtr(0),
ref_weights.matrix().n_rows());
}
else {
pWeights.Init(rset_.n_cols());
pWeights.SetAll(1);
}
if(!strcmp(qfname, rfname)) {
qset_.Alias(rset_);
}
else {
Dataset query_dataset;
query_dataset.InitFromFile(qfname);
qset_.Own(&(query_dataset.matrix()));
}
// scale dataset if the user wants to
if(!strcmp(fx_param_str(NULL, "scaling", NULL), "range")) {
scale_data_by_minmax(qset_, rset_);
}
Vector pSources;
pSources.Init(rset_.n_rows() * rset_.n_cols());
concatenate_vectors(rset_, pSources);
Vector pTargets;
pTargets.Init(qset_.n_rows() * qset_.n_cols());
concatenate_vectors(qset_, pTargets);
double Bandwidth = fx_param_double_req(NULL, "bandwidth");
Vector pGaussTransform;
pGaussTransform.Init(qset_.n_cols());
double epsilon = fx_param_double(NULL, "tau", 0.1);
// choose parameters
fx_timer_start(NULL, "ifgt_kde_compute");
ImprovedFastGaussTransformChooseParameters cp(rset_.n_rows(), Bandwidth,
epsilon,
(int) ceil(0.2 * 100 /
sqrt
(2 * Bandwidth *
Bandwidth)));
printf("Number of clusters chosen: %d\n", cp.K);
printf("Maximum truncation number: %d\n", cp.p_max);
printf("Maximum cutoff radius: %g\n", cp.r);
// run k-center clustering
ArrayList<int> pClusterIndex;
pClusterIndex.Init(rset_.n_cols());
KCenterClustering kc(rset_.n_rows(), rset_.n_cols(), pSources.ptr(),
pClusterIndex.begin(), cp.K);
kc.Cluster();
ArrayList<int> pNumPoints;
pNumPoints.Init(cp.K);
Vector pClusterCenter;
pClusterCenter.Init(rset_.n_rows() * cp.K);
Vector pClusterRadii;
pClusterRadii.Init(cp.K);
kc.ComputeClusterCenters(cp.K, pClusterCenter.ptr(), pNumPoints.begin(),
pClusterRadii.ptr());
// update truncation number
ImprovedFastGaussTransformChooseTruncationNumber ct(rset_.n_rows(),
Bandwidth, epsilon,
kc.MaxClusterRadius);
// initialize IFGT instance
ArrayList<int> pTruncNumber;
pTruncNumber.Init(rset_.n_cols());
for(index_t i = 0; i < rset_.n_cols(); i++) {
pTruncNumber[i] = 0;
}
ImprovedFastGaussTransform* pIFGT = new
ImprovedFastGaussTransform(rset_.n_rows(), rset_.n_cols(), qset_.n_cols(),
pSources.ptr(), Bandwidth, pWeights.ptr(),
pTargets.ptr(), ct.p_max,
cp.K, pClusterIndex.begin(),
pClusterCenter.ptr(), pClusterRadii.ptr(),
cp.r, epsilon, pGaussTransform.ptr(),
pTruncNumber.begin());
// run IFGT
pIFGT->Evaluate();
GaussianKernel kernel;
kernel.Init(Bandwidth);
double norm_const = kernel.CalcNormConstant(qset_.n_rows()) *
rset_.n_cols();
// normalize density estimates
for(index_t q = 0; q < qset_.n_cols(); q++) {
pGaussTransform[q] /= norm_const;
}
fx_timer_stop(NULL, "ifgt_kde_compute");
// check answer with naive
NaiveKde<GaussianKernel> naive_kde;
naive_kde.Init(qset_, rset_);
naive_kde.Compute();
naive_kde.ComputeMaximumRelativeError(pGaussTransform);
delete pIFGT;
fx_done();
return 0;
}
fastlib/u/gmravi/regression/kcenter_clustering.cc
-317
View File
@@ -1,317 +0,0 @@
//-------------------------------------------------------------------
// The code was written by Changjiang Yang and Vikas Raykar
// and is copyrighted under the Lesser GPL:
//
// Copyright (C) 2006 Changjiang Yang and Vikas Raykar
//
// This program is free software; you can redistribute it and/or modify
// it under the terms of the GNU Lesser General Public License as
// published by the Free Software Foundation; version 2.1 or later.
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
// See the GNU Lesser General Public License for more details.
// You should have received a copy of the GNU Lesser General Public
// License along with this program; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place - Suite 330, Boston,
// MA 02111-1307, USA.
//
// The author may be contacted via email at:cyang(at)sarnoff(.)com
// vikas(at)umiacs(.)umd(.)edu
//-------------------------------------------------------------------
//-------------------------------------------------------------------
// File : KCenterClustering.cpp
// Purpose : Implementation for the k-center clustering algorithm.
// Author : Vikas C. Raykar (vikas@cs.umd.edu)
// Date : April 25 2005, June 10 2005, August 23, 2005
//-------------------------------------------------------------------
#include "kcenter_clustering.h"
#include <math.h>
#include <time.h>
#include <stdlib.h>
#include "fastlib/fastlib_int.h"
//-------------------------------------------------------------------
// Constructor
//
// PURPOSE
// -------
// Initialize the class.
// Read the parameters.
//
// INPUT
// ----------
// Dim --> dimension of the points.
// NSources --> number of sources.
// pSources --> pointer to sources, (d*N).
// pClusterIndex --> pointer to a vector of length N where the
// i th element is the cluster number to
// which the i th point belongs.
//
//-------------------------------------------------------------------
KCenterClustering::KCenterClustering(int Dim,
int NSources,
double *pSources,
int *pClusterIndex,
int NumClusters
)
{
//Read the parameters
d=Dim;
N=NSources;
px=pSources;
pci=pClusterIndex;
K=NumClusters;
dist_C = new double[N]; //distances to the center.
r=new double[K];
}
//-------------------------------------------------------------------
// Destructor
//-------------------------------------------------------------------
KCenterClustering::~KCenterClustering()
{
delete [] dist_C;
delete [] r;
}
//-------------------------------------------------------------------
// ddist is the square of the distance of two vectors(double)
//-------------------------------------------------------------------
double
KCenterClustering::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;
}
//-------------------------------------------------------------------
// Find the largest element from a vector
//-------------------------------------------------------------------
int
KCenterClustering::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-center Clustering.
//-------------------------------------------------------------------
//
// Gonzalez's farthest-point clustering algorithm.
//
// OUTPUT
// ----------------
//
// MaxClusterRadius --> maximum radius of the clusters, (rx).
// pci --> vector of length N where the i th element is the
// cluster number to which the i th point belongs.
// pci[i] varies between 0 to K-1.
//-------------------------------------------------------------------
void
KCenterClustering::Cluster()
{
int *pCenters = new int[K]; //indices of the centers.
int *cprev = new int[N]; // index to the previous node
int *cnext = new int[N]; // index to the next node
int *far2c = new int[K]; // farthest node to the center
// randomly pick one node as the first center.
srand( (unsigned)time( NULL ) );
int nc = rand() % N; // new center
// add the ind-th node to the first center.
pCenters[0] = nc;
// compute the distances from each node to the first center.
// initialize the circular linked list, the center is the
// sentinel node.
const double *x_nc, *x_j;
x_nc = px + nc*d;
x_j = px;
for (int j = 0; j < N; x_j += d, j++)
{
dist_C[j] = (j==nc)? 0.0:ddist(d, x_j, x_nc);
cnext[j] = j+1;
cprev[j] = j-1;
// my fix (by Dongryeol Lee)
pci[j] = 0;
}
cnext[N-1] = 0; // link the tail to the head.
cprev[0] = N-1; // link the head to the tail.
// compute the radius of the first cluster and the farthest
// node to the center.
nc = idmax(N,dist_C);
far2c[0] = nc;
r[0] = dist_C[nc];
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. It is a new center.
nc = idmax(i,r);
nc = far2c[nc];
pCenters[i] = nc; //add the ind-th node to the current center.
r[i] = dist_C[nc] = 0.0;pci[nc]=i;
far2c[i] = nc;
cnext[cprev[nc]] = cnext[nc]; // delete nc
cprev[cnext[nc]] = cprev[nc];
cnext[nc] = cprev[nc] = nc; //self-loop
//update the distances from each point to the current center.
x_nc = px + nc*d;
for (int j = 0; j < i; j++)
{
int ct_j = pCenters[j];
x_j = px + ct_j*d;
double dc2cq = ddist(d, x_j, x_nc) / 4;
if (dc2cq < r[j]) // neighbor cluster
{
r[j] = 0.0;
far2c[j] = ct_j;
int k = cnext[ct_j];
while (k != ct_j) // visit the circular linked list
{
int nextk = cnext[k];
//compare the distances from new center
//and from current center.
double dist2c_k = dist_C[k];
if ( dc2cq < dist2c_k )
{
x_j = px + k*d;
double dd = ddist(d, x_j, x_nc);
if ( dd < dist2c_k )
{
dist_C[k] = dd; // update distances to center
pci[k]=i;
if (r[i] < dd) // find max r
{
r[i] = dd;
far2c[i] = k;
}
cnext[cprev[k]] = nextk; // delete nextk from ct_j
cprev[nextk] = cprev[k];
cnext[k] = cnext[nc]; // insert nextk to nc
cprev[cnext[nc]] = k;
cnext[nc] = k;
cprev[k] = nc;
}
else if ( r[j] < dist2c_k )
{
r[j] = dist2c_k;
far2c[j] = k;
}
}
else if ( r[j] < dist2c_k )
{
r[j] = dist2c_k;
far2c[j] = k;
} // if d < 2 r_k
k = nextk;
} // while k
} // if d < 2 r
} // for j
} // for i
nc = idmax(K,r);
MaxClusterRadius=sqrt(r[nc]);
}
//------------------------------------------------------------------------
// Computes
// [1] the cluster centers by taking the mean of all the points
// belonging to a cluster.
// [2] the number of points in each cluster.
// [3] the radius of each cluster.
//------------------------------------------------------------------------
// NumClusters --> number of clusters
// pClusterCenters --> pointer to the cluster centers, (d*K),
// pNumPoints --> pointer to the num of points in each cluster, (K).
// pClusterRadii --> pointer to the radius of each cluster, (K).
//------------------------------------------------------------------------
void
KCenterClustering::ComputeClusterCenters(
int NumClusters,
double *pClusterCenters,
int *pNumPoints,
double *pClusterRadii
)
{
int K=NumClusters;
for(int k=0; k<K; k++)
{
pNumPoints[k]=0;
pClusterRadii[k]=sqrt(r[k]);
for(int dim=0; dim<d; dim++)
{
pClusterCenters[(k*d)+dim]=0.0;
}
}
for(int i=0; i<N; i++)
{
pNumPoints[pci[i]] += 1;
for(int dim=0; dim<d; dim++)
{
pClusterCenters[(pci[i]*d)+dim] += px[(i*d)+dim];
}
}
for(int k=0; k<K; k++)
{
for(int dim=0; dim<d; dim++)
{
pClusterCenters[(k*d)+dim]=pClusterCenters[(k*d)+dim]/pNumPoints[k];
}
}
}
fastlib/u/gmravi/regression/kcenter_clustering.h
-118
View File
@@ -1,118 +0,0 @@
//-------------------------------------------------------------------
// The code was written by Changjiang Yang and Vikas Raykar
// and is copyrighted under the Lesser GPL:
//
// Copyright (C) 2006 Changjiang Yang and Vikas Raykar
//
// This program is free software; you can redistribute it and/or modify
// it under the terms of the GNU Lesser General Public License as
// published by the Free Software Foundation; version 2.1 or later.
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
// See the GNU Lesser General Public License for more details.
// You should have received a copy of the GNU Lesser General Public
// License along with this program; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place - Suite 330, Boston,
// MA 02111-1307, USA.
//
// The author may be contacted via email at:cyang(at)sarnoff(.)com
// vikas(at)umiacs(.)umd(.)edu
//-------------------------------------------------------------------
//----------------------------------------------------------------------------
// File : KCenterClustering.h
// Purpose : Interface for the k-center clustering algorithm.
// Author : Vikas C. Raykar (vikas@cs.umd.edu)
// Date : April 25 2005, June 10 2005, August 23, 2005
//
//----------------------------------------------------------------------------
// Gonzalez's farthest-point clustering algorithm.
//
// June 10, 2005:
// This version now returns the number points and the radius of each cluster.
//
// August 23, 2005:
// Speed up using the doubly circular list.
// The clusters far away are trimmed. The nodes inside the neighboring
// clusters which are within half sphere are trimmed.
// The computational complexity is reduced to O(n log k).
//
//----------------------------------------------------------------------------
//
// INPUT
// ----------------
//
// Dim --> dimension of the points.
// NSources --> number of sources.
// pSources --> pointer to sources, (d*N).
// NumClusters --> number of clusters.
//
// OUTPUT
// ----------------
//
// MaxClusterRadius --> maximum radius of the clusters, (rx).
// pClusterIndex --> vector of length N where the i th element is the
// cluster number to which the i th point belongs.
// pClusterIndex[i] varies between 0 to K-1.
// pClusterCenters --> pointer to the cluster centers, (d*K).
// pNumPoints --> pointer to the number of points in each cluster, (K).
// pClusterRadii --> pointer to the radius of each cluster, (K).
//----------------------------------------------------------------------------
#ifndef K_CENTER_CLUSTERING_H
#define K_CENTER_CLUSTERING_H
class KCenterClustering{
public:
//Output parameters
double MaxClusterRadius; //maximum cluster radius
//Functions
//constructor
KCenterClustering(int Dim,
int NSources,
double *pSources,
int *pClusterIndex,
int NumClusters
);
//destructor
~KCenterClustering();
//k-center clustering
void Cluster();
//Compute cluster centers and the number of points in each cluster
//and the radius of each cluster.
void ComputeClusterCenters(int NumClusters,
double *pClusterCenters,
int *pNumPoints,
double *pClusterRadii);
private:
//Input Parameters
int d; //dimension of the points.
int N; //number of sources.
double *px; //pointer to sources, (d*N).
int K; //number of clusters
int *pci; //pointer to a vector of length N where the i th element is the
//cluster number to which the i th point belongs.
double *dist_C; //distances to the center.
double *r;
//Functions
double ddist(const int d, const double *x, const double *y);
int idmax(int n, double *x);
};
#endif
fastlib/u/gmravi/regression/kde.h
-973
View File
@@ -1,973 +0,0 @@
#ifndef KDE_H
#define KDE_H
#include "fastlib/fastlib_int.h"
#include "u/dongryel/series_expansion/farfield_expansion.h"
#include "u/dongryel/series_expansion/local_expansion.h"
#include "u/dongryel/series_expansion/mult_farfield_expansion.h"
#include "u/dongryel/series_expansion/mult_local_expansion.h"
#include "u/dongryel/series_expansion/kernel_aux.h"
template<typename TKernel>
class NaiveKde {
private:
/** query dataset */
Matrix qset_;
/** reference dataset */
Matrix rset_;
/** kernel */
TKernel kernel_;
/** computed densities */
Vector densities_;
public:
void Compute() {
printf("\nStarting naive KDE...\n");
fx_timer_start(NULL, "naive_kde_compute");
// compute unnormalized sum
for(index_t q = 0; q < qset_.n_cols(); q++) {
const double *q_col = qset_.GetColumnPtr(q);
for(index_t r = 0; r < rset_.n_cols(); r++) {
const double *r_col = rset_.GetColumnPtr(r);
double dsqd = la::DistanceSqEuclidean(qset_.n_rows(), q_col, r_col);
densities_[q] += kernel_.EvalUnnormOnSq(dsqd);
}
}
// then normalize it
double norm_const = kernel_.CalcNormConstant(qset_.n_rows()) *
rset_.n_cols();
for(index_t q = 0; q < qset_.n_cols(); q++) {
densities_[q] /= norm_const;
}
fx_timer_stop(NULL, "naive_kde_compute");
printf("\nNaive KDE completed...\n");
}
void Init() {
densities_.SetZero();
}
void Init(Matrix &qset, Matrix &rset) {
// get datasets
qset_.Alias(qset);
rset_.Alias(rset);
// get bandwidth
kernel_.Init(fx_param_double_req(NULL, "bandwidth"));
// allocate density storage
densities_.Init(qset.n_cols());
densities_.SetZero();
}
void PrintDebug() {
FILE *stream = stdout;
const char *fname = NULL;
if(fx_param_exists(NULL, "naive_kde_output")) {
fname = fx_param_str(NULL, "naive_kde_output", NULL);
stream = fopen(fname, "w+");
}
for(index_t q = 0; q < qset_.n_cols(); q++) {
fprintf(stream, "%g\n", densities_[q]);
}
if(stream != stdout) {
fclose(stream);
}
}
void ComputeMaximumRelativeError(const Vector &density_estimate) {
double max_rel_err = 0;
for(index_t q = 0; q < densities_.length(); q++) {
double rel_err = fabs(density_estimate[q] - densities_[q]) /
densities_[q];
if(rel_err > max_rel_err) {
max_rel_err = rel_err;
}
}
fx_format_result(NULL, "maxium_relative_error_for_fast_KDE", "%g",
max_rel_err);
}
};
template<typename TKernel, typename TKernelAux>
class FastKde {
public:
// forward declaration of KdeStat class
class KdeStat;
// our tree type using the KdeStat
typedef BinarySpaceTree<DHrectBound<2>, Matrix, KdeStat > Tree;
class KdeStat {
public:
/** lower bound on the densities for the query points owned by this node
*/
double mass_l_;
/**
* additional offset for the lower bound on the densities for the query
* points owned by this node (for leaf nodes only).
*/
double more_l_;
/**
* lower bound offset passed from above
*/
double owed_l_;
/** stores the portion pruned by finite difference
*/
double mass_e_;
/** upper bound on the densities for the query points owned by this node
*/
double mass_u_;
/**
* additional offset for the upper bound on the densities for the query
* points owned by this node (for leaf nodes only)
*/
double more_u_;
/**
* upper bound offset passed from above
*/
double owed_u_;
/** extra error that can be used for the query points in this node */
double mass_t_;
/**
* Far field expansion created by the reference points in this node.
*/
typename TKernelAux::TFarFieldExpansion farfield_expansion_;
/**
* Local expansion stored in this node.
*/
typename TKernelAux::TLocalExpansion local_expansion_;
/** Initialize the statistics */
void Init() {
mass_l_ = 0;
more_l_ = 0;
owed_l_ = 0;
mass_e_ = 0;
mass_u_ = 0;
more_u_ = 0;
owed_u_ = 0;
mass_t_ = 0;
}
void Init(double bandwidth,
typename TKernelAux::TSeriesExpansionAux *sea) {
farfield_expansion_.Init(bandwidth, sea);
local_expansion_.Init(bandwidth, sea);
}
void Init(const Matrix& dataset, index_t &start, index_t &count) {
Init();
}
void Init(const Matrix& dataset, index_t &start, index_t &count,
const KdeStat& left_stat,
const KdeStat& right_stat) {
Init();
}
void Init(double bandwidth, const Vector& center,
typename TKernelAux::TSeriesExpansionAux *sea) {
farfield_expansion_.Init(bandwidth, center, sea);
local_expansion_.Init(bandwidth, center, sea);
Init();
}
void MergeChildBounds(KdeStat &left_stat, KdeStat &right_stat) {
// steal left and right children's tokens
double min_mass_t = min(left_stat.mass_t_, right_stat.mass_t_);
// improve lower and upper bound
mass_l_ = max(mass_l_, min(left_stat.mass_l_, right_stat.mass_l_));
mass_u_ = min(mass_u_, max(left_stat.mass_u_, right_stat.mass_u_));
mass_t_ += min_mass_t;
left_stat.mass_t_ -= min_mass_t;
right_stat.mass_t_ -= min_mass_t;
}
void PushDownTokens
(KdeStat &left_stat, KdeStat &right_stat, double *de,
typename TKernelAux::TLocalExpansion *local_expansion, double *dt) {
if(de != NULL) {
double de_ref = *de;
left_stat.mass_e_ += de_ref;
right_stat.mass_e_ += de_ref;
*de = 0;
}
if(local_expansion != NULL) {
local_expansion->TranslateToLocal(left_stat.local_expansion_);
local_expansion->TranslateToLocal(right_stat.local_expansion_);
}
if(dt != NULL) {
double dt_ref = *dt;
left_stat.mass_t_ += dt_ref;
right_stat.mass_t_ += dt_ref;
*dt = 0;
}
}
KdeStat() { }
~KdeStat() {}
};
private:
/** series expansion auxililary object */
typename TKernelAux::TSeriesExpansionAux sea_;
/** query dataset */
Matrix qset_;
/** query tree */
Tree *qroot_;
/** reference dataset */
Matrix rset_;
/** reference tree */
Tree *rroot_;
/** reference weights */
Vector rset_weights_;
/** list of kernels to evaluate */
TKernel kernel_;
/** lower bound on the densities */
Vector densities_l_;
/** densities computed */
Vector densities_e_;
/** upper bound on the densities */
Vector densities_u_;
/** accuracy parameter */
double tau_;
int num_farfield_to_local_prunes_;
int num_farfield_prunes_;
int num_local_prunes_;
int num_finite_difference_prunes_;
// preprocessing: scaling the dataset; this has to be moved to the dataset
// module
/* scales each attribute to 0-1 using the min/max values */
void scale_data_by_minmax() {
int num_dims = rset_.n_rows();
DHrectBound<2> qset_bound;
DHrectBound<2> rset_bound;
qset_bound.Init(qset_.n_rows());
rset_bound.Init(qset_.n_rows());
// go through each query/reference point to find out the bounds
for(index_t r = 0; r < rset_.n_cols(); r++) {
Vector ref_vector;
rset_.MakeColumnVector(r, &ref_vector);
rset_bound |= ref_vector;
}
for(index_t q = 0; q < qset_.n_cols(); q++) {
Vector query_vector;
qset_.MakeColumnVector(q, &query_vector);
qset_bound |= query_vector;
}
for(index_t i = 0; i < num_dims; i++) {
DRange qset_range = qset_bound.get(i);
DRange rset_range = rset_bound.get(i);
double min_coord = min(qset_range.lo, rset_range.lo);
double max_coord = max(qset_range.hi, rset_range.hi);
double width = max_coord - min_coord;
printf("Dimension %d range: [%g, %g]\n", i, min_coord, max_coord);
for(index_t j = 0; j < rset_.n_cols(); j++) {
rset_.set(i, j, (rset_.get(i, j) - min_coord) / width);
}
if(strcmp(fx_param_str(NULL, "query", NULL),
fx_param_str_req(NULL, "data"))) {
for(index_t j = 0; j < qset_.n_cols(); j++) {
qset_.set(i, j, (qset_.get(i, j) - min_coord) / width);
}
}
}
}
// member functions
void UpdateBounds(Tree *qnode, Tree *rnode,
double *dl, double *de, double *du, double *dt,
int *order_farfield_to_local, int *order_farfield,
int *order_local) {
// query self statistics
KdeStat &qstat = qnode->stat();
// reference node statistics
KdeStat &rstat = rnode->stat();
// incorporate into the self
double dl_ref = *dl;
double du_ref = *du;
qstat.mass_l_ += dl_ref;
qstat.mass_u_ += du_ref;
// incorporate finite difference pruning if available
if(de != NULL) {
qstat.mass_e_ += (*de);
}
// incorporate token change
if(dt != NULL) {
qstat.mass_t_ += (*dt);
}
// incorporate series approximation into the self
// far field to local translation
if(order_farfield_to_local != NULL && *order_farfield_to_local >= 0) {
rstat.farfield_expansion_.TranslateToLocal(qstat.local_expansion_,
*order_farfield_to_local);
}
// far field pruning
else if(order_farfield != NULL && *order_farfield >= 0) {
for(index_t q = qnode->begin(); q < qnode->end(); q++) {
densities_e_[q] +=
rstat.farfield_expansion_.EvaluateField(&qset_, q, NULL,
*order_farfield);
}
}
// local accumulation pruning
else if(order_local != NULL && *order_local >= 0) {
qstat.local_expansion_.AccumulateCoeffs(rset_, rset_weights_,
rnode->begin(), rnode->end(),
*order_local);
}
// for a leaf node, incorporate the lower and upper bound changes into
// its additional offset
if(qnode->is_leaf())
{
qstat.more_l_ += dl_ref;
qstat.more_u_ += du_ref;
}
// otherwise, incorporate the bound changes into the owed slots of
// the immediate descendants
else {
qnode->left()->stat().owed_l_ += dl_ref; //transmission of the owed valuesto the children
qnode->left()->stat().owed_u_ += du_ref;
qnode->right()->stat().owed_l_ += dl_ref;
qnode->right()->stat().owed_u_ += du_ref;
}
}
/** exhaustive base KDE case */
void FKdeBase(Tree *qnode, Tree *rnode) {
// compute unnormalized sum
for(index_t q = qnode->begin(); q < qnode->end(); q++) {
// get query point
const double *q_col = qset_.GetColumnPtr(q);
for(index_t r = rnode->begin(); r < rnode->end(); r++) {
// get reference point
const double *r_col = rset_.GetColumnPtr(r);
// pairwise distance and kernel value
double dsqd = la::DistanceSqEuclidean(qset_.n_rows(), q_col, r_col);
double ker_value = kernel_.EvalUnnormOnSq(dsqd);
densities_l_[q] += ker_value;
densities_e_[q] += ker_value;
densities_u_[q] += ker_value;
}
}
// tally up the unused error components due to exhaustive computation
qnode->stat().mass_t_ += rnode->count();
// get a tighter lower and upper bound by looping over each query point
// in the current query leaf node
double min_l = MAXDOUBLE;
double max_u = -MAXDOUBLE;
for(index_t q = qnode->begin(); q < qnode->end(); q++) {
if(densities_l_[q] < min_l) {
min_l = densities_l_[q];
}
if(densities_u_[q] > max_u) {
max_u = densities_u_[q];
}
}
// subtract the contribution accounted by the exhaustive computation
qnode->stat().more_u_ -= rnode->count();
// tighten lower and upper bound
qnode->stat().mass_l_ = min_l + qnode->stat().more_l_;
qnode->stat().mass_u_ = max_u + qnode->stat().more_u_;
}
/**
* checking for prunability of the query and the reference pair using
* four types of pruning methods
*/
int PrunableEnhanced(Tree *qnode, Tree *rnode, DRange &dsqd_range,
DRange &kernel_value_range, double &dl, double &du,
double &dt, int &order_farfield_to_local,
int &order_farfield, int &order_local) {
int dim = rset_.n_rows();
// actual amount of error incurred per each query/ref pair
double actual_err_farfield_to_local = 0;
double actual_err_farfield = 0;
double actual_err_local = 0;
// estimated computational cost
int cost_farfield_to_local = MAXINT;
int cost_farfield = MAXINT;
int cost_local = MAXINT;
int cost_exhaustive = (qnode->count()) * (rnode->count()) * dim;
int min_cost = 0;
// query node and reference node statistics
KdeStat &qstat = qnode->stat();
KdeStat &rstat = rnode->stat();
// expansion objects
typename TKernelAux::TFarFieldExpansion &farfield_expansion =
rstat.farfield_expansion_;
typename TKernelAux::TLocalExpansion &local_expansion =
qstat.local_expansion_;
// number of reference points
int num_references = rnode->count();
// try pruning after bound refinement:
// the new lower bound after incorporating new info
dl = kernel_value_range.lo * num_references;
du = -kernel_value_range.hi * num_references;
// refine the lower bound using the new lower bound info
double new_mass_l = qstat.mass_l_ + dl;
double allowed_err = tau_ * new_mass_l *
((double)(num_references + qstat.mass_t_)) /
((double) rroot_->count() * num_references);
// get the order of approximations
order_farfield_to_local =
farfield_expansion.OrderForConvertingToLocal
(rnode->bound(), qnode->bound(), dsqd_range.lo, dsqd_range.hi,
allowed_err, &actual_err_farfield_to_local);
order_farfield =
farfield_expansion.OrderForEvaluating(rnode->bound(), qnode->bound(),
dsqd_range.lo, dsqd_range.hi,
allowed_err, &actual_err_farfield);
order_local =
local_expansion.OrderForEvaluating(rnode->bound(), qnode->bound(),
dsqd_range.lo, dsqd_range.hi,
allowed_err, &actual_err_local);
// update computational cost and compute the minimum
if(order_farfield_to_local >= 0) {
cost_farfield_to_local = (int) pow(order_farfield_to_local + 1,
2 * dim);
}
if(order_farfield >= 0) {
cost_farfield = (int) pow(order_farfield + 1, dim) * (qnode->count());
}
if(order_local >= 0) {
cost_local = (int) pow(order_local + 1, dim) * (rnode->count());
}
min_cost = min(cost_farfield_to_local,
min(cost_farfield, min(cost_local, cost_exhaustive)));
if(cost_farfield_to_local == min_cost) {
dt = num_references *
(1.0 - (rroot_->count()) * actual_err_farfield_to_local /
(new_mass_l * tau_));
order_farfield = order_local = -1;
num_farfield_to_local_prunes_++;
return 1;
}
if(cost_farfield == min_cost) {
dt = num_references *
(1.0 - (rroot_->count()) * actual_err_farfield / (new_mass_l * tau_));
order_farfield_to_local = order_local = -1;
num_farfield_prunes_++;
return 1;
}
if(cost_local == min_cost) {
dt = num_references *
(1.0 - (rroot_->count()) * actual_err_local / (new_mass_l * tau_));
order_farfield_to_local = order_farfield = -1;
num_local_prunes_++;
return 1;
}
order_farfield_to_local = order_farfield = order_local = -1;
dl = du = dt = 0;
return 0;
}
/** checking for prunability of the query and the reference pair */
int Prunable(Tree *qnode, Tree *rnode, DRange &dsqd_range,
DRange &kernel_value_range, double &dl, double &de,
double &du, double &dt) {
// query node stat
KdeStat &stat = qnode->stat();
// number of reference points
int num_references = rnode->count();
// try pruning after bound refinement: first compute distance/kernel
// value bounds
dsqd_range.lo = qnode->bound().MinDistanceSq(rnode->bound());
dsqd_range.hi = qnode->bound().MaxDistanceSq(rnode->bound());
kernel_value_range = kernel_.RangeUnnormOnSq(dsqd_range);
// the new lower bound after incorporating new info
dl = kernel_value_range.lo * num_references;
de = 0.5 * num_references *
(kernel_value_range.lo + kernel_value_range.hi);
du = -kernel_value_range.hi * num_references;
// refine the lower bound using the new lower bound info
double new_mass_l = stat.mass_l_ + dl;
double allowed_err = tau_ * new_mass_l *
((double)(num_references + stat.mass_t_)) / ((double) rroot_->count());
// this is error per each query/reference pair for a fixed query
double m = 0.5 * (kernel_value_range.hi - kernel_value_range.lo);
// this is total error for each query point
double error = m * num_references;
// check pruning condition
if(error <= allowed_err) {
dt = num_references *
(1.0 - (rroot_->count()) * m / (new_mass_l * tau_));
num_finite_difference_prunes_++;
return 1;
}
else {
dl = de = du = dt = 0;
return 0;
}
}
/** determine which of the node to expand first */
void BestNodePartners(Tree *nd, Tree *nd1, Tree *nd2, Tree **partner1,
Tree **partner2) {
double d1 = nd->bound().MinDistanceSq(nd1->bound());
double d2 = nd->bound().MinDistanceSq(nd2->bound());
if(d1 <= d2) {
*partner1 = nd1;
*partner2 = nd2;
}
else {
*partner1 = nd2;
*partner2 = nd1;
}
}
/** canonical fast KDE case */
void FKde(Tree *qnode, Tree *rnode) {
/** temporary variable for storing lower bound change */
double dl = 0, de = 0, du = 0, dt = 0;
int order_farfield_to_local = -1, order_farfield = -1, order_local = -1;
// temporary variable for holding distance/kernel value bounds
DRange dsqd_range;
DRange kernel_value_range;
// query node statistics
KdeStat &stat = qnode->stat();
// left child and right child of query node statistics
KdeStat *left_stat = NULL;
KdeStat *right_stat = NULL;
// process density bound changes sent from the ancestor query nodes,
UpdateBounds(qnode, rnode, &stat.owed_l_, NULL, &stat.owed_u_, NULL,
NULL, NULL, NULL);
// for non-leaf query node, tighten lower/upper bounds and the
// reclaim tokens unused by the children.
if(!qnode->is_leaf()) {
left_stat = &(qnode->left()->stat());
right_stat = &(qnode->right()->stat());
stat.MergeChildBounds(*left_stat, *right_stat);
}
// try finite difference pruning first
if(Prunable(qnode, rnode, dsqd_range, kernel_value_range,
dl, de, du, dt)) {
UpdateBounds(qnode, rnode, &dl, &de, &du, &dt, NULL, NULL, NULL);
return;
}
// try series-expansion pruning
else if(PrunableEnhanced(qnode, rnode, dsqd_range, kernel_value_range,
dl, du, dt, order_farfield_to_local,
order_farfield, order_local)) {
UpdateBounds(qnode, rnode, &dl, NULL, &du, &dt,
&order_farfield_to_local, &order_farfield,
&order_local);
return;
}
// for leaf query node
if(qnode->is_leaf()) {
// for leaf pairs, go exhaustive
if(rnode->is_leaf()) {
FKdeBase(qnode, rnode);
return;
}
// for non-leaf reference, expand reference node
else {
Tree *rnode_first = NULL, *rnode_second = NULL;
BestNodePartners(qnode, rnode->left(), rnode->right(), &rnode_first,
&rnode_second);
FKde(qnode, rnode_first);
FKde(qnode, rnode_second);
return;
}
}
// for non-leaf query node
else {
// for a leaf reference node, expand query node
if(rnode->is_leaf()) {
Tree *qnode_first = NULL, *qnode_second = NULL;
stat.PushDownTokens(*left_stat, *right_stat, NULL, NULL,
&stat.mass_t_);
BestNodePartners(rnode, qnode->left(), qnode->right(), &qnode_first,
&qnode_second);
FKde(qnode_first, rnode);
FKde(qnode_second, rnode);
return;
}
// for non-leaf reference node, expand both query and reference nodes
else {
Tree *rnode_first = NULL, *rnode_second = NULL;
stat.PushDownTokens(*left_stat, *right_stat, NULL, NULL,
&stat.mass_t_);
BestNodePartners(qnode->left(), rnode->left(), rnode->right(),
&rnode_first, &rnode_second);
FKde(qnode->left(), rnode_first);
FKde(qnode->left(), rnode_second);
BestNodePartners(qnode->right(), rnode->left(), rnode->right(),
&rnode_first, &rnode_second);
FKde(qnode->right(), rnode_first);
FKde(qnode->right(), rnode_second);
return;
}
}
}
/**
* pre-processing step - this wouldn't be necessary if the core
* fastlib supported a Init function for Stat objects that take
* more arguments.
*/
void PreProcess(Tree *node) {
// initialize the center of expansions and bandwidth for
// series expansion
node->stat().Init(sqrt(kernel_.bandwidth_sq()), &sea_);
node->bound().CalculateMidpoint
(node->stat().farfield_expansion_.get_center());
node->bound().CalculateMidpoint
(node->stat().local_expansion_.get_center());
// initialize lower bound to 0
node->stat().mass_l_ = 0;
// set the finite difference approximated amounts to 0
node->stat().mass_e_ = 0;
// set the upper bound to the number of reference points
node->stat().mass_u_ = rset_.n_cols();
// set the number of tokens to 0
node->stat().mass_t_ = 0;
// for non-leaf node, recurse
if(!node->is_leaf()) {
node->stat().owed_l_ = node->stat().owed_u_ = 0;
PreProcess(node->left());
PreProcess(node->right());
// translate multipole moments
node->stat().farfield_expansion_.TranslateFromFarField
(node->left()->stat().farfield_expansion_);
node->stat().farfield_expansion_.TranslateFromFarField
(node->right()->stat().farfield_expansion_);
}
else {
node->stat().more_l_ = node->stat().more_u_ = 0;
// exhaustively compute multipole moments
node->stat().farfield_expansion_.RefineCoeffs(rset_, rset_weights_,
node->begin(), node->end(),
sea_.get_max_order());
}
}
/** post processing step */
void PostProcess(Tree *qnode) {
KdeStat &stat = qnode->stat();
// for leaf query node
if(qnode->is_leaf()) {
for(index_t q = qnode->begin(); q < qnode->end(); q++) {
densities_e_[q] +=
stat.local_expansion_.EvaluateField(&qset_, q, NULL) +
stat.mass_e_;
}
}
else {
// push down approximations
stat.PushDownTokens(qnode->left()->stat(), qnode->right()->stat(),
&stat.mass_e_, &stat.local_expansion_, NULL);
PostProcess(qnode->left());
PostProcess(qnode->right());
}
}
void NormalizeDensities() {
double norm_const = kernel_.CalcNormConstant(qset_.n_rows()) *
rset_.n_cols();
for(index_t q = 0; q < qset_.n_cols(); q++) {
densities_l_[q] /= norm_const;
densities_e_[q] /= norm_const;
densities_u_[q] /= norm_const;
}
}
public:
// constructor/destructor
FastKde() {}
~FastKde() {
if(qroot_ != rroot_ ) {
delete qroot_;
delete rroot_;
}
else {
delete rroot_;
}
}
// getters and setters
/** get the reference dataset */
Matrix &get_reference_dataset() { return rset_; }
/** get the query dataset */
Matrix &get_query_dataset() { return qset_; }
/** get the density estimate */
const Vector &get_density_estimates() { return densities_e_; }
// interesting functions...
void Compute(double tau) {
// set accuracy parameter
tau_ = tau;
// initialize the lower and upper bound densities
densities_l_.SetZero();
densities_e_.SetZero();
densities_u_.SetAll(rset_.n_cols());
num_finite_difference_prunes_ = num_farfield_to_local_prunes_ =
num_farfield_prunes_ = num_local_prunes_ = 0;
printf("\nStarting fast KDE...\n");
fx_timer_start(NULL, "fast_kde_compute");
// preprocessing step for initializing series expansion objects
PreProcess(rroot_);
if(qroot_ != rroot_) {
PreProcess(qroot_);
}
// call main routine
FKde(qroot_, rroot_);
// postprocessing step for finalizing the sums
PostProcess(qroot_);
// normalize densities
NormalizeDensities();
fx_timer_stop(NULL, "fast_kde_compute");
printf("\nFast KDE completed...\n");
printf("Finite difference prunes: %d\n", num_finite_difference_prunes_);
printf("F2L prunes: %d\n", num_farfield_to_local_prunes_);
printf("F prunes: %d\n", num_farfield_prunes_);
printf("L prunes: %d\n", num_local_prunes_);
}
void Init() {
Dataset ref_dataset;
// read in the number of points owned by a leaf
int leaflen = fx_param_int(NULL, "leaflen", 20);
// read the datasets
const char *rfname = fx_param_str_req(NULL, "data");
const char *qfname = fx_param_str(NULL, "query", rfname);
// read reference dataset
ref_dataset.InitFromFile(rfname);
rset_.Own(&(ref_dataset.matrix()));
// read the reference weights
char *rwfname = NULL;
if(fx_param_exists(NULL, "dwgts")) {
rwfname = (char *)fx_param_str(NULL, "dwgts", NULL);
}
if(rwfname != NULL) {
Dataset ref_weights;
ref_weights.InitFromFile(rwfname);
rset_weights_.Copy(ref_weights.matrix().GetColumnPtr(0),
ref_weights.matrix().n_rows());
}
else {
rset_weights_.Init(rset_.n_cols());
rset_weights_.SetAll(1);
}
if(!strcmp(qfname, rfname)) {
qset_.Alias(rset_);
}
else {
Dataset query_dataset;
query_dataset.InitFromFile(qfname);
qset_.Own(&(query_dataset.matrix()));
}
// scale dataset if the user wants to
if(!strcmp(fx_param_str(NULL, "scaling", NULL), "range")) {
scale_data_by_minmax();
}
// construct query and reference trees
fx_timer_start(NULL, "tree_d");
rroot_ = tree::MakeKdTreeMidpoint<Tree>(rset_, leaflen);
if(!strcmp(qfname, rfname)) {
qroot_ = rroot_;
}
else {
qroot_ = tree::MakeKdTreeMidpoint<Tree>(qset_, leaflen);
}
fx_timer_stop(NULL, "tree_d");
// initialize the density lists
densities_l_.Init(qset_.n_cols());
densities_e_.Init(qset_.n_cols());
densities_u_.Init(qset_.n_cols());
// initialize the kernel
kernel_.Init(fx_param_double_req(NULL, "bandwidth"));
// initialize the series expansion object
if(qset_.n_rows() <= 2) {
sea_.Init(fx_param_int(NULL, "order", 5), qset_.n_rows());
}
else {
sea_.Init(fx_param_int(NULL, "order", 0), qset_.n_rows());
}
}
void PrintDebug() {
FILE *stream = stdout;
const char *fname = NULL;
if((fname = fx_param_str(NULL, "fast_kde_output", NULL)) != NULL) {
stream = fopen(fname, "w+");
}
for(index_t q = 0; q < qset_.n_cols(); q++) {
fprintf(stream, "%g\n", densities_e_[q]);
}
if(stream != stdout) {
fclose(stream);
}
}
};
#endif
fastlib/u/gmravi/regression/kde_test_to_dongryeol.nb
-17881
View File
File diff suppressed because it is too large Load Diff
fastlib/u/gmravi/regression/kg.m
-3
View File
@@ -1,3 +0,0 @@
function k=kg(x)
k=exp(-x .* x ./2)/sqrt(2*pi);
end
fastlib/u/gmravi/regression/lscvscore.m
-16
View File
@@ -1,16 +0,0 @@
function moh=lscvscore(x,h)
% x is the input data
% h is the bandwidth
n=length(x);
sqterm=0;
xterm=0;
for i=1:n
sqterm=sqterm+sum(kg((x-x(i))/(sqrt(2)*h)))/sqrt(2);
% The sqrt(2) factors are for the K(2) term which is
% equivalent to a kernel with variance of 2 that is the convolution of 2
% gaussian kernels.
xterm=xterm+sum(kg((x-x(i))/h));
end;
sqterm=sqterm/(n*n*h);
xterm=2*(xterm/(n*n)-kg(0)/n)/h;
moh=sqterm-xterm;
fastlib/u/gmravi/regression/main2.cc
+20 -1
View File
@@ -5,10 +5,29 @@
int main (int argc, char *argv[]){
fx_init (argc, argv);
char *rfname=(char*)malloc(40);
char *qfname=(char*)malloc(40);
strcpy(rfname,fx_param_str_req (NULL, "data"));
strcpy(qfname,fx_param_str_req (NULL,"query"));
Regression2 <GaussianKernel> reg2;
printf("going to initialization function...\n");
reg2.Init();
Dataset ref_dataset ;
Dataset q_dataset;
Matrix query_dataset;
Matrix reference_dataset;
ref_dataset.InitFromFile(rfname);
reference_dataset.Own(&(ref_dataset.matrix()));
q_dataset.InitFromFile(qfname);
query_dataset.Own(&(q_dataset.matrix()));
reg2.Init(query_dataset,reference_dataset);
printf("Initializations done..\n");
// reg2.Compute(fx_param_double (NULL, "tau", 0.1));
reg2.Compute(0.1);
ArrayList<Matrix> wfkde_results;
fastlib/u/gmravi/regression/naive_kde_output.txt
-25
View File
@@ -1,25 +0,0 @@
Point:1.300000 1.000000 2.000000
2.671527
3.385464
4.383671
5.647777
Point:1.100000 1.800000 1.900000
3.408173
4.289175
6.271699
7.380817
Point:2.100000 2.200000 1.500000
2.294841
3.540102
4.491305
4.387145
Point:1.000000 2.000000 3.000000
2.910201
3.260119
5.737805
7.191543
Point:1.000000 2.400000 2.600000
3.094637
3.571118
6.217940
7.430222
fastlib/u/gmravi/regression/query_file.csv
+4 -5
View File
@@ -1,5 +1,4 @@
1.0, 2.0, 3.0
2.1, 2.2, 1.5
1.10, 1.8, 1.9
1.30, 1.0, 2.0
1.0, 2.4, 2.6
1.0, 2.0, 3.2
1.3, 2.4, 3.6
1.5, 4.5, 7.6
2.3, 3.6, 8.9
1 1.0 2.0 3.0 3.2
2 2.1 1.3 2.2 2.4 1.5 3.6
3 1.10 1.5 1.8 4.5 1.9 7.6
4 1.30 2.3 1.0 3.6 2.0 8.9
1.0 2.4 2.6