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

regression bugs still not fixed

Browse Source
This commit is contained in:
gmravi2003
2007-12-02 23:26:50 +00:00
parent d5aab3d9a5
commit ef10bd0322
26 changed files with 74624 additions and 574 deletions
Download Patch File Download Diff File Expand all files Collapse all files
fastlib/u/gmravi/allknn/allknn.h
+767 -473
View File
File diff suppressed because it is too large Load Diff
fastlib/u/gmravi/allknn/allnn.h
+27 -23
View File
@@ -12,10 +12,20 @@ class AllNNNaive
class Stat
{
//This is not used in case of single tree operations. used only in case of dual tree operations
//This is justa dummy class
private:
double distance_max; //This gives the max distance within which all neighbours should be found for all points
public:
void Init() {}
void set_maximum_distance(double distance)
{
this->distance_max=distance;
}
double get_maximum_distance()
{
return distance_max;
}
void Init() {set_maximum_distance(32768.0);}
void Init(const Matrix& dataset, index_t start, index_t count)
{
@@ -86,7 +96,7 @@ class AllKNNSingleTreeResults
AllKNNSingleTreeResults()
{
}
};
@@ -102,41 +112,35 @@ class AllKNNSingleTree
// results_matrix.astr.Init();
printf("Initialized properly\n");
}
friend bool same_points(double *point1,double *point2, int len);
//This is the function which will perform the actual single tree algorithm and spew the results
void ComputeAllKNNSingleTree(Matrix*,Matrix*,int); // A 2-D matrix is returned
int FindKNearestNeighbours(Tree*,double*,Matrix*,SingleTreeResults*,int,int); //A 1-D matrix of the k nearest neighbours is returned
int FindKNearestNeighbours(Tree*,double*,Matrix*,ArrayList<SingleTreeResults>&,int,int); //A 1-D matrix of the k nearest neighbours is returned
};
class AllKNNDualTreeResults //This is just the same as AllKNNSingleTreeResults
{
public:
ArrayList<ArrayList<SingleTreeResults> > astr;
};
classAllKNNDualTree
class AllKNNDualTree
{
public:
AllKNNDualTreeResults results_matrix;
//constructor
//cosntructor definition
public:
AllKNNDualTreeResults results_matrix;
AllKNNDualTree()
{
//Initialize results_matrix_ This will later be reinitialized to the number of columns in the function ComputeAllKNNSingleTree
// results_matrix.astr.Init();
printf("Initialized properly\n");
}
//This is the function which will perform the actual single tree algorithm and spew the results
void ComputeAllKNNSingleTree(Matrix*,Matrix*,int); // A 2-D matrix is returned
int FindKNearestNeighbours(Tree*,double*,Matrix*,SingleTreeResults*,int,int); //A 1-D matrix of the k nearest neighbours is returned
}
friend bool same_points(double*, double*);
//This is the function which will perform the actual dual tree algorithm and spew the results
void ComputeAllKNNDualTree(Matrix*,Matrix*,int); // A 2-D matrix is returned
double FindKNearestNeighboursDualTree(Tree *q_tree,Tree *r_tree,Matrix*,Matrix*,int k); //A 1-D matrix of the k nearest neighbours is returned
};
fastlib/u/gmravi/allknn/build.py
+7 -7
View File
@@ -1,16 +1,16 @@
librule(
name = "allknn_mine", # this line can be safely omitted
sources = [], # files that must be compiled
headers = ["allnn.h"], # include files part of the 'lib'
name = "allknn", # this line can be safely omitted
sources = [""], # files that must be compiled
headers = ["allknn.h"], # include files part of the 'lib'
deplibs = ["fastlib:fastlib_int"] # depends on fastlib core
)
binrule(
name = "main", # the executable name
sources = ["main.cc"], # compile main.cc
headers = [], # no extra headers
deplibs = [":allknn_mine"] #
name = "allknn_main", # the executable name
sources = ["allknn_main.cc"], # compile main.cc
headers = ["allknn.h"], # no extra headers
deplibs = ["fastlib:fastlib_int"]
)
# to build:
fastlib/u/gmravi/allknn/main.cc
+302 -71
View File
@@ -7,10 +7,10 @@ bool check_if_equal(ArrayList <double> arr1, ArrayList <double> arr2)
for(i=0;i<arr1.size();i++)
{
if(arr1[i]!=arr2[i])
return 0;
return false;
}
return 1;
}
return true;
}
void AllNNNaive::Compute(Matrix *q_matrix, Matrix *r_matrix)
{
//for each column vector of q_matrix find the nearest neighbour in r_matrix
@@ -175,18 +175,24 @@ void AllNNSingleTree::ComputeAllNNSingleTree(Matrix* q_matrix, Matrix *r_matrix
//return str_;
}
int find_index(SingleTreeResults *str,double dist,int start,int end)
int find_index(ArrayList<SingleTreeResults> str,double dist,int start,int end)
{
printf("Start is %d and end is %d\n",start,end);
printf("Searching for distance %f\n",dist);
printf("In function find index str is...\n");
for(int j=start;j<=end;j++)
{
printf("distance is %f\n",str[j].get_distance());
}
if(start>end)
return 0;
if(start==end)
{
//printf("str[start].get_distance() is %f",str[start].get_distance());
//printf("dist is %f\n",dist);
printf("str[start].get_distance() is %f",str[start].get_distance());
printf("dist is %f\n",dist);
if(dist>str[start].get_distance())
{
// printf("will return %d\n",end+1);
printf("will return %d\n",end+1);
return end+1; //returned if the element to be added is at the end of the list
}
else return start;
@@ -200,30 +206,32 @@ int find_index(SingleTreeResults *str,double dist,int start,int end)
if(dist<str[(end+start)/2].get_distance())
{
//go left
//printf("left\n");
printf("left\n");
return find_index(str,dist,start,(start+end)/2);
}
else
{
//go right
// printf("right\n");
//printf("New start is %d and new end is %d\n",(start+end)/2+1,end);
printf("right\n");
printf("New start is %d and new end is %d\n",(start+end)/2+1,end);
return find_index(str,dist,(start+end)/2+1,end);
}
}
}
int push_into_array(SingleTreeResults *str,int position, double dist,int length, int index,int k)
int push_into_array(ArrayList<SingleTreeResults> &str,int position, double dist,int length, int index,int k)
{
//printf("Length is %d\n",length);
//printf("Index to be inserted is %d\n",index);
//printf("Distance to be inserted is %f\n",dist);
//printf("The position is %d\n",position);
printf("Length is %d\n",length);
printf("Index to be inserted is %d\n",index);
printf("Distance to be inserted is %f\n",dist);
printf("The position is %d\n",position);
if(position==length) //that means add the element to the end of list
{
if(length==k)
{
//cannot push into array
for(int l=0;l<length;l++)
printf("distance:%f and index:%d\n",str[l].get_distance(),str[l].get_index());
return 0;
}
else
@@ -231,21 +239,29 @@ int push_into_array(SingleTreeResults *str,int position, double dist,int length,
//add it to the end of the array
str[length].set_result(dist,index);
for(int l=0;l<length+1;l++)
printf("distance:%f and index:%d\n",str[l].get_distance(),str[l].get_index());
return 1;
}
}
ArrayList <SingleTreeResults> temp;
temp.Init(k);
// printf("temp created\n");
printf("temp created\n");
//the element will be added in the middle of the array
//printf("The element will be added to the middle of the array\n");
printf("The element will be added to the middle of the array\n");
for(int j=0;j<position;j++)
temp[j].set_result(str[j].get_distance(),str[j].get_index());
printf("Distance is %f\n",dist);
printf("index is %d\n",index);
temp[position].set_result(dist,index);
//printf("Length is %d and k=%d\n",length,k);
printf("Length is %d and k=%d\n",length,k);
printf("temp at position is set up\n");
if(length==k)
{
@@ -255,22 +271,57 @@ int push_into_array(SingleTreeResults *str,int position, double dist,int length,
for(int j=0;j<length;j++)
str[j].set_result(temp[j].get_distance(),temp[j].get_index());
// delete(temp);
printf("Before returning i have str as\n");
for(int l=0;l<length;l++)
printf("distance:%f and index:%d\n",str[l].get_distance(),str[l].get_index());
return 0;
}
else
{
for(int t=length;t>position;t--)
temp[t].set_result(str[t-1].get_distance(),str[t-1].get_index());
// printf("adjusted.....");
printf("since length is not equal to k\n");
for(int t=length;t>position;t-=1)
{
printf("came here\n");
printf("Length is %d\n",length);
printf("position is %d\n",position);
printf("will copy to temp\n");
printf("t is %d\n",t);
printf("distance is %f\n",str[t-1].get_distance());
printf("index is %d\n",str[t-1].get_index());
temp[t].set_result(str[t-1].get_distance(),str[t-1].get_index());
printf("Set up..\n");
printf("t is %d\n",t);
printf("wi.ll loop...\n");
}
printf("adjusrwrewerwerwerwerwe.....");
for(int j=0;j<length+1;j++)
str[j].set_result(temp[j].get_distance(),temp[j].get_index());
{
printf("copying element by element...\n");
printf("woll copy distance %f\n",temp[j].get_distance());
printf("Will copy index =%d\n",temp[j].get_index());
printf("j is %d\n",j);
str[j].set_result(temp[j].get_distance(),temp[j].get_index());
printf("copied....\n");
}
printf("copied\n");
// delete(temp);
// printf("deleted temp");
printf("str before leaving is..\n");
for(int l=0;l<length+1;l++)
printf("distance:%f and index:%d\n",str[l].get_distance(),str[l].get_index());
return 1;
}
}
bool same_points(double *point1,double *point2, int len)
bool same_points(double *point1,double *point2,int len)
{
for(int i=0;i<len;i++)
{
@@ -286,7 +337,7 @@ bool same_points(double *point1,double *point2, int len)
return true;
}
int AllKNNSingleTree::FindKNearestNeighbours(Tree *root, double *point,Matrix *r_matrix,SingleTreeResults *str,int length,int k)
int AllKNNSingleTree::FindKNearestNeighbours(Tree *root, double *point,Matrix *r_matrix,ArrayList<SingleTreeResults> &str,int length,int k)
{
//Base case is that the node is a leaf
@@ -308,18 +359,21 @@ int AllKNNSingleTree::FindKNearestNeighbours(Tree *root, double *point,Matrix *r
if(!same_points(point,r_matrix->GetColumnPtr(i),r_matrix->n_rows()))
{
dist=la::DistanceSqEuclidean(r_matrix->n_rows(),point,r_matrix->GetColumnPtr(i));
//printf("Distance is %f\n",dist);
printf("Distance is %f\n",dist);
//find where the new element should be pushed
int start=0;
printf("Before going to funcction find index ..\n");
for(int l=0;l<length;l++)
printf("distance is %f\n",str[l].get_distance());
position=find_index(str,dist,start,end);
//printf("The position where this willbe inserted is \n");
//printf("%d\n",position);
printf("The position where this willbe inserted is \n");
printf("%d\n",position);
length+=push_into_array(str,position,dist,length,i,k);
end=length-1;
/*printf("After push...\n");
printf("After pushAFTWER PUSH AFTER PUSH...\n");
for(int count=0;count<length;count++)
printf("distance:%f\n",str[count].get_distance());*/
printf("distance:%f\n",str[count].get_distance());
//printf("length is now %d\n",length);
if(length>k)
@@ -404,7 +458,8 @@ int AllKNNSingleTree::FindKNearestNeighbours(Tree *root, double *point,Matrix *r
//Now for each query point find the k nearest neighbours
SingleTreeResults *str=new SingleTreeResults[k];
ArrayList<SingleTreeResults> str;
str.Init(k);
//printf("Will start knn calculations\n");
@@ -413,11 +468,8 @@ int AllKNNSingleTree::FindKNearestNeighbours(Tree *root, double *point,Matrix *r
for(int i=0;i<r_matrix->n_cols();i++)
{
int length=0;
// printf("came here...\n");
//printf("root=%X\n",root);
//printf("the number of points covered are %d\n",i);
//printf("Total number of points left to be covered are %d\n",r_matrix->n_cols()-i);
length=FindKNearestNeighbours(root,q_matrix->GetColumnPtr(i),r_matrix,str,length,k);
length+=FindKNearestNeighbours(root,q_matrix->GetColumnPtr(i),r_matrix,str,length,k);
printf("The nearest neighbour distances are\n");
@@ -431,12 +483,7 @@ int AllKNNSingleTree::FindKNearestNeighbours(Tree *root, double *point,Matrix *r
results_matrix.astr[i][j].set_result(str[j].get_distance(),str[j].get_index());
}
/* printf("knn for this point have been added...........\n");
for(int j=0;j<length;j++)
{
printf("Index:%d distance: %f\n",str[j].get_index(),str[j].get_distance());
}*/
printf("The origianl point is\n");
@@ -446,30 +493,215 @@ int AllKNNSingleTree::FindKNearestNeighbours(Tree *root, double *point,Matrix *r
printf("Will clear str now........\n");
printf("Length is %d\n",length);
/* for(int l=0;l<length;l++)
str[l].set_result(32768.0,-1);*/
}
printf("The knn distances are\n");
printf("The number of points are %d\n",r_matrix->n_cols());
char fname[40];
strcpy(fname,"allknn");
char query_file[40];
strcat(fname,"_colors50k.ds");
strcat(fname,query_file);
printf("The name of query_file is %s\n",query_file);
FILE *fp=fopen(fname,"w");
for(int t=0;t<r_matrix->n_cols();t++)
{
for(int l=0;l<k;l++)
fprintf(fp,"%f ",results_matrix.astr[t][l].get_distance());
fprintf(fp,"\n");
}
}
void test_module()
double AllKNNDualTree::FindKNearestNeighboursDualTree(Tree *q_tree,Tree *r_tree,Matrix *q_matrix,Matrix *r_matrix,int k)
{
//if distance between the two boxes is larger than the max _distance then return
double distance_between_boxes=q_tree->bound().MinDistanceSq (r_tree->bound()); //base case
if(q_tree->is_leaf()&& r_tree->is_leaf())
{
int start,end;
//check if pruneable
if(q_tree->stat().get_maximum_distance() < distance_between_boxes)
{
//then there is no need to go further and hence we can return
return 32768.0;
}
else
{
//not purneable. therefore carry out exhaustive point-to-point computations
double max_dist=0.0;
double distance;
int position;
int count=0;
for(int i=q_tree->begin();i<q_tree->end();i++)
{
printf("for this point count is %d\n",count);
printf("TAKING NEW I=%d............................................................\n",i);
count=0;
for(int j=r_tree->begin();j<r_tree->end();j++)
{
if(!same_points(q_matrix->GetColumnPtr(i),r_matrix->GetColumnPtr(j),q_matrix->n_rows()))
{
distance=la::DistanceSqEuclidean (r_matrix->n_rows(),q_matrix->GetColumnPtr(i),r_matrix->GetColumnPtr(j));
//we would like to find index of this point into str. this will enter into str[i]. This function takes in as argument an array list of single tree results
printf("distance is %f\n",distance);
printf("size is %d\n",results_matrix.astr[i].size());
if(results_matrix.astr[i].size()==0)
{
//initialize start and end
printf("Initialized start and end\n");
start=0;
end=-1;
}
int length=results_matrix.astr[i].size(); //this calculates the old length
printf("at the moment length is %d\n",length);
position=find_index(results_matrix.astr[i],distance,start,results_matrix.astr[i].size()-1);
if(length<k)
{
printf("since the size of array is still lesser than k will increase the size\n");
printf("Before expanding length is %d \n",length);
results_matrix.astr[i].AddBack(1); //increase the length of the arraylist, only if the number of elements as of now are lesser than k
printf("new length after expanding is %d\n",results_matrix.astr[i].size());
}
push_into_array(results_matrix.astr[i],position,distance,length,j,k);
printf("pushed into array..\n");
count++;
}
}
if(results_matrix.astr[i].size()<k)
max_dist=32768.0;
else //all k nn have been found
max_dist=max_dist > results_matrix.astr[i][k-1].get_distance()?max_dist:results_matrix.astr[i][k-1].get_distance(); //this is an update after every point IN THE QUERY NODE
for(int z=0;z<results_matrix.astr[i].size();z++)
printf("distance=%f and index=%d\n",results_matrix.astr[i][z].get_distance(),results_matrix.astr[i][z].get_index());
}
printf("Will return from this node a value of %f%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\n",max_dist);
return max_dist;
}
}
else
{
//Not base case. PRUNING CASE
if(q_tree->stat().get_maximum_distance() < distance_between_boxes)
{
//then there is no need to go further and henraghavce we can return
return 32768.0;
}
//NOT PRUNEABLE
//both are not leafs
if(!q_tree->is_leaf() && !r_tree->is_leaf())
{
printf("both are not leafs......\n");
double max_dist_q_left_r_left=FindKNearestNeighboursDualTree(q_tree->left(),r_tree->left(),q_matrix,r_matrix,k);
/* if(max_dist < q_tree->stat().get_maximum_distance())
{
q_tree->stat().set_maximum_distance(max_dist);
}*/
double max_dist_q_left_r_right=FindKNearestNeighboursDualTree(q_tree->left(),r_tree->right(),q_matrix,r_matrix,k);
double max_dist_q_left= max_dist_q_left_r_left< max_dist_q_left_r_right? max_dist_q_left_r_left: max_dist_q_left_r_right;//take minimum
/* if(max_dist_q_left< q_tree->stat().get_maximum_distance())
{
q_tree->stat().set_maximum_distance(max_dist);
}*/
double max_dist_q_right_r_left=FindKNearestNeighboursDualTree(q_tree->right(),r_tree->left(),q_matrix,r_matrix,k);
/*if(max_dist<q_tree->stat().get_maximum_distance())
{
q_tree->stat().set_maximum_distance(max_dist);
}*/
double max_dist_q_right_r_right=FindKNearestNeighboursDualTree(q_tree->right(),r_tree->right(),q_matrix,r_matrix,k);
double max_dist_q_right= max_dist_q_right_r_left< max_dist_q_right_r_right? max_dist_q_right_r_left: max_dist_q_right_r_right;
double max_dist_q=max_dist_q_left>max_dist_q_right?max_dist_q_left:max_dist_q_right;
q_tree->stat().set_maximum_distance(max_dist_q);
return q_tree->stat().get_maximum_distance();
}
else
{
printf("q is leaf and r is not...\n");
//q_tree is leaf and r_tree is not
if(q_tree->is_leaf()&&!r_tree->is_leaf())
{
double max_dist_q_r_left=FindKNearestNeighboursDualTree(q_tree,r_tree->left(),q_matrix,r_matrix,k);
/*if(max_dist < q_tree->stat().get_maximum_distance())
{
q_tree->stat().set_maximum_distance(max_dist);
}*/
double max_dist_q_r_right=FindKNearestNeighboursDualTree(q_tree,r_tree->right(),q_matrix,r_matrix,k);
/*if(max_dist<q_tree->stat().get_maximum_distance())
{
q_tree->stat().set_maximum_distance(max_dist);
}*/
double max_dist_q = max_dist_q_r_left > max_dist_q_r_right?max_dist_q_r_right:max_dist_q_r_left; //take minimum
q_tree->stat().set_maximum_distance(max_dist_q);
return q_tree->stat().get_maximum_distance();
}
else
{
printf("q is not a leaf and r is \n");
//q_tree is not a leaf and r_tree is
if(!q_tree->is_leaf()&&r_tree->is_leaf())
{
double max_dist_q_left_r=FindKNearestNeighboursDualTree(q_tree->left(),r_tree,q_matrix,r_matrix,k);
//printf("After getting back to the function the value is %f\n",max_dist);
printf("will compare this vlaue with %f\n",q_tree->stat().get_maximum_distance());
/* if(max_dist < q_tree->stat().get_maximum_distance())
{
q_tree->stat().set_maximum_distance(max_dist);
}*/
double max_dist_q_right_r=FindKNearestNeighboursDualTree(q_tree->right(),r_tree,q_matrix,r_matrix,k);
double max_dist_q= max_dist_q_left_r > max_dist_q_right_r? max_dist_q_left_r: max_dist_q_right_r;
q_tree->stat().set_maximum_distance(max_dist_q);
return max_dist_q;
}
}
}
}
}
void AllKNNDualTree::ComputeAllKNNDualTree(Matrix *q_matrix,Matrix *r_matrix,int k)
{
//printf("results_matrix initialized\n");
// Build the kdtree form the reference matrix and the query matrix
using namespace tree;
Tree *r_tree=tree::MakeKdTreeMidpoint<Tree>(*r_matrix,3,NULL,NULL); //here the leaf size is set to 3
printf("ref tree built\n");
Tree *q_tree=tree::MakeKdTreeMidpoint<Tree>(*r_matrix,3,NULL,NULL); //here the leaf size is set to 3
printf("query tree built\n");
printf("Trees built\n");
//****************************************************************
results_matrix.astr.Init(r_matrix->n_cols()); //the arraylist ahs been completely initialized
printf("size has been set to %d\n",r_matrix->n_cols());
for(int i=0;i<r_matrix->n_cols();i++)
results_matrix.astr[i].Init();
FindKNearestNeighboursDualTree(q_tree,r_tree,q_matrix,r_matrix,k);
printf("completed the algo. Printing results\n");
printf("the k nearest neighbours are............\n");
for(int l=0;l<r_matrix->n_cols();l++)
{
for(int t=0;t<k;t++)
{printf("distance is %f ",results_matrix.astr[l][t].get_distance());}
printf("\n");
}
//****************************************************************
}
int main(int argc, char *argv[])
{
@@ -490,8 +722,8 @@ int main(int argc, char *argv[])
//strcpy(ref_file,fx_param_str_req(NULL,"ref_file"));
//strcpy(method,fx_param_str_req(NULL,"method"));
strcpy(query_file,"colors50k.ds");
strcpy(ref_file,"colors50k.ds");
strcpy(query_file,"dummy.ds");
strcpy(ref_file,"dummy.ds");
int k=3;
@@ -507,7 +739,7 @@ int main(int argc, char *argv[])
//if method is allnnnaive. AllNNNaive computations follow
//printf("String compare with allknnsingletree is %d\n",strcmp(method,"allknnsingletree"));
if(strcmp(method,"allnnnaive")==0)
if(strcmp(method,"allnnnaive")==0)
{
printf("In allnaive");
AllNNNaive naive;
@@ -541,15 +773,14 @@ int main(int argc, char *argv[])
akst->ComputeAllKNNSingleTree(q_matrix,r_matrix,k);
}
//the code that follows is just a test module
if(strcmp(method,"test"))
/*
if(strcmp(method,"allknndualtree")==0)
{
printf("came here");
test_module();
}
if(strcmp(method,"allnndual")==0)
{
}
printf("in allknndualtree\n");
//AllKNNSDual tree computations
//int k=fx_param_int_req(NULL,"k");
AllKNNDualTree *akdt= new AllKNNDualTree();
akdt->ComputeAllKNNDualTree(q_matrix,r_matrix,k);
}*/
fx_done();
}
fastlib/u/gmravi/regression/Makefile
+203
View File
@@ -0,0 +1,203 @@
/net/hc295/gmravi/home/fastlib/fastlib/u/gmravi/regression/main: /net/hc295/gmravi/home/fastlib/fastlib/u/gmravi/regression/Makefile /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/u/gmravi/regression/main
@echo '... Making /net/hc295/gmravi/home/fastlib/fastlib/u/gmravi/regression/main'
@rm -f /net/hc295/gmravi/home/fastlib/fastlib/u/gmravi/regression/main
@ln -s -f /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/u/gmravi/regression/main /net/hc295/gmravi/home/fastlib/fastlib/u/gmravi/regression/main
@echo '*** Created 1 symlinks in /net/hc295/gmravi/home/fastlib/fastlib/u/gmravi/regression/.'
/net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/u/gmravi/regression/main: /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/libcol_col.a /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/libdata_data.a /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/libu_gmravi_regression_main__auto.a /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/libbase_base.a /net/hc295/gmravi/home/fastlib/fastlib/bin_keep/x86_64_Linux_COMMON_gcc_COMMON/libblaspack.a /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/libfx_fx.a /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/libla_la.a /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/libpar_par.a /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/libmath_math.a /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/libfile_file.a
@echo '... Making bin/x86_64_Linux_debug_gcc_/u/gmravi/regression/main'
@mkdir -p /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/u/gmravi/regression
@g++ -o /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/u/gmravi/regression/main -g3 -DDEBUG -O0 /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/libu_gmravi_regression_main__auto.a /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/libpar_par.a /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/libmath_math.a /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/libdata_data.a /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/libfile_file.a /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/libfx_fx.a /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/libla_la.a /net/hc295/gmravi/home/fastlib/fastlib/bin_keep/x86_64_Linux_COMMON_gcc_COMMON/libblaspack.a /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/libcol_col.a /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/libbase_base.a -lm -lpthread -lg2c
pseudo_4: /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/libcol_col.a /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/libdata_data.a /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/libu_gmravi_regression_main__auto.a /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/libbase_base.a /net/hc295/gmravi/home/fastlib/fastlib/bin_keep/x86_64_Linux_COMMON_gcc_COMMON/libblaspack.a /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/libfx_fx.a /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/libla_la.a /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/libpar_par.a /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/libmath_math.a /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/libfile_file.a
@echo '*** Done with /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/libcol_col.a /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/libdata_data.a /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/libu_gmravi_regression_main__auto.a /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/libbase_base.a /net/hc295/gmravi/home/fastlib/fastlib/bin_keep/x86_64_Linux_COMMON_gcc_COMMON/libblaspack.a /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/libfx_fx.a /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/libla_la.a /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/libpar_par.a /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/libmath_math.a /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/libfile_file.a'
/net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/libu_gmravi_regression_main__auto.a: /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/obj/u_gmravi_regression_main.o
@echo '... Making bin/x86_64_Linux_debug_gcc_/libu_gmravi_regression_main__auto.a'
@mkdir -p /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_
@ar r /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/libu_gmravi_regression_main__auto.a /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/obj/u_gmravi_regression_main.o
/net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/obj/u_gmravi_regression_main.o: /net/hc295/gmravi/home/fastlib/fastlib/u/gmravi/regression/main.cc /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_COMMON_gcc_COMMON/libfastlib_fastlib_int.h /net/hc295/gmravi/home/fastlib/fastlib/u/gmravi/regression/regression.h
@echo '... Making bin/x86_64_Linux_debug_gcc_/obj/u_gmravi_regression_main.o'
@mkdir -p /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/obj
@cd /net/hc295/gmravi/home/fastlib/fastlib/u/gmravi/regression && g++ -Wall -Woverloaded-virtual -fno-exceptions -Wparentheses -fno-exceptions -I/net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_COMMON_gcc_COMMON -I/net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_COMMON_gcc_ -I/net/hc295/gmravi/home/fastlib/fastlib -g3 -DDEBUG -O0 -c main.cc -o /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/obj/u_gmravi_regression_main.o
/net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_COMMON_gcc_COMMON/libfastlib_fastlib_int.h: /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_COMMON_gcc_COMMON/libdata_data.h /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_COMMON_gcc_COMMON/libtree_tree.h /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_COMMON_gcc_COMMON/libbase_base.h /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_COMMON_gcc_COMMON/libfx_fx.h /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_COMMON_gcc_COMMON/libmath_math.h /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_COMMON_gcc_COMMON/libcol_col.h /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_COMMON_gcc_COMMON/libfile_file_int.h /net/hc295/gmravi/home/fastlib/fastlib/fastlib/fastlib.h /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_COMMON_gcc_COMMON/libpar_par.h /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_COMMON_gcc_COMMON/libla_la.h
@echo '... Making bin/x86_64_Linux_COMMON_gcc_COMMON/libfastlib_fastlib_int.h'
@mkdir -p /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_COMMON_gcc_COMMON
@touch /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_COMMON_gcc_COMMON/libfastlib_fastlib_int.h
/net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_COMMON_gcc_COMMON/libpar_par.h: /net/hc295/gmravi/home/fastlib/fastlib/par/grain.h /net/hc295/gmravi/home/fastlib/fastlib/par/task.h /net/hc295/gmravi/home/fastlib/fastlib/par/thread.h /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_COMMON_gcc_COMMON/libcol_col.h /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_COMMON_gcc_COMMON/libbase_base.h
@echo '... Making bin/x86_64_Linux_COMMON_gcc_COMMON/libpar_par.h'
@mkdir -p /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_COMMON_gcc_COMMON
@touch /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_COMMON_gcc_COMMON/libpar_par.h
/net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_COMMON_gcc_COMMON/libtree_tree.h: /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_COMMON_gcc_COMMON/libdata_data.h /net/hc295/gmravi/home/fastlib/fastlib/tree/kdtree_impl.h /net/hc295/gmravi/home/fastlib/fastlib/tree/spacetree.h /net/hc295/gmravi/home/fastlib/fastlib/tree/kdtree.h /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_COMMON_gcc_COMMON/libbase_base.h /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_COMMON_gcc_COMMON/libfx_fx.h /net/hc295/gmravi/home/fastlib/fastlib/tree/bounds.h /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_COMMON_gcc_COMMON/libcol_col.h /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_COMMON_gcc_COMMON/libfile_file_int.h /net/hc295/gmravi/home/fastlib/fastlib/tree/statistic.h /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_COMMON_gcc_COMMON/libla_la.h
@echo '... Making bin/x86_64_Linux_COMMON_gcc_COMMON/libtree_tree.h'
@mkdir -p /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_COMMON_gcc_COMMON
@touch /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_COMMON_gcc_COMMON/libtree_tree.h
/net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_COMMON_gcc_COMMON/libmath_math.h: /net/hc295/gmravi/home/fastlib/fastlib/math/math.h /net/hc295/gmravi/home/fastlib/fastlib/math/discrete.h /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_COMMON_gcc_COMMON/libcol_col.h /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_COMMON_gcc_COMMON/libbase_base.h /net/hc295/gmravi/home/fastlib/fastlib/math/kernel.h /net/hc295/gmravi/home/fastlib/fastlib/math/geometry.h
@echo '... Making bin/x86_64_Linux_COMMON_gcc_COMMON/libmath_math.h'
@mkdir -p /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_COMMON_gcc_COMMON
@touch /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_COMMON_gcc_COMMON/libmath_math.h
/net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_COMMON_gcc_COMMON/libdata_data.h: /net/hc295/gmravi/home/fastlib/fastlib/data/dataset.h /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_COMMON_gcc_COMMON/libfile_file.h /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_COMMON_gcc_COMMON/libla_la.h /net/hc295/gmravi/home/fastlib/fastlib/data/crossvalidation.h /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_COMMON_gcc_COMMON/libbase_base.h
@echo '... Making bin/x86_64_Linux_COMMON_gcc_COMMON/libdata_data.h'
@mkdir -p /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_COMMON_gcc_COMMON
@touch /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_COMMON_gcc_COMMON/libdata_data.h
/net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_COMMON_gcc_COMMON/libfile_file_int.h: /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_COMMON_gcc_COMMON/libfile_file.h
@echo '... Making bin/x86_64_Linux_COMMON_gcc_COMMON/libfile_file_int.h'
@mkdir -p /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_COMMON_gcc_COMMON
@touch /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_COMMON_gcc_COMMON/libfile_file_int.h
/net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_COMMON_gcc_COMMON/libfx_fx.h: /net/hc295/gmravi/home/fastlib/fastlib/fx/fx.h /net/hc295/gmravi/home/fastlib/fastlib/fx/timer.h /net/hc295/gmravi/home/fastlib/fastlib/fx/datastore.h /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_COMMON_gcc_COMMON/libbase_base.h
@echo '... Making bin/x86_64_Linux_COMMON_gcc_COMMON/libfx_fx.h'
@mkdir -p /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_COMMON_gcc_COMMON
@touch /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_COMMON_gcc_COMMON/libfx_fx.h
pseudo_22: /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/libcol_col.a /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/libdata_data.a /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/libbase_base.a /net/hc295/gmravi/home/fastlib/fastlib/bin_keep/x86_64_Linux_COMMON_gcc_COMMON/libblaspack.a /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/libfx_fx.a /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/libla_la.a /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/libpar_par.a /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/libmath_math.a /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/libfile_file.a
@echo '*** Done with /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/libcol_col.a /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/libdata_data.a /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/libbase_base.a /net/hc295/gmravi/home/fastlib/fastlib/bin_keep/x86_64_Linux_COMMON_gcc_COMMON/libblaspack.a /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/libfx_fx.a /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/libla_la.a /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/libpar_par.a /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/libmath_math.a /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/libfile_file.a'
pseudo_24: /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/libcol_col.a /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/libpar_par.a /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/libbase_base.a
@echo '*** Done with /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/libcol_col.a /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/libpar_par.a /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/libbase_base.a'
/net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/libpar_par.a: /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/obj/par_thread.o
@echo '... Making bin/x86_64_Linux_debug_gcc_/libpar_par.a'
@mkdir -p /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_
@ar r /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/libpar_par.a /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/obj/par_thread.o
/net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/obj/par_thread.o: /net/hc295/gmravi/home/fastlib/fastlib/par/task.h /net/hc295/gmravi/home/fastlib/fastlib/par/thread.h /net/hc295/gmravi/home/fastlib/fastlib/par/thread.cc /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_COMMON_gcc_COMMON/libcol_col.h /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_COMMON_gcc_COMMON/libbase_base.h /net/hc295/gmravi/home/fastlib/fastlib/par/grain.h
@echo '... Making bin/x86_64_Linux_debug_gcc_/obj/par_thread.o'
@mkdir -p /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/obj
@cd /net/hc295/gmravi/home/fastlib/fastlib/par && g++ -Wall -Woverloaded-virtual -fno-exceptions -Wparentheses -fno-exceptions -I/net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_COMMON_gcc_COMMON -I/net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_COMMON_gcc_ -I/net/hc295/gmravi/home/fastlib/fastlib -g3 -DDEBUG -O0 -c thread.cc -o /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/obj/par_thread.o
pseudo_31: /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/libcol_col.a /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/libdata_data.a /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/libbase_base.a /net/hc295/gmravi/home/fastlib/fastlib/bin_keep/x86_64_Linux_COMMON_gcc_COMMON/libblaspack.a /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/libfx_fx.a /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/libla_la.a /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/libfile_file.a
@echo '*** Done with /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/libcol_col.a /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/libdata_data.a /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/libbase_base.a /net/hc295/gmravi/home/fastlib/fastlib/bin_keep/x86_64_Linux_COMMON_gcc_COMMON/libblaspack.a /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/libfx_fx.a /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/libla_la.a /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/libfile_file.a'
pseudo_33: /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/libcol_col.a /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/libmath_math.a /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/libbase_base.a
@echo '*** Done with /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/libcol_col.a /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/libmath_math.a /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/libbase_base.a'
/net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/libmath_math.a: /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/obj/math_geometry.o /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/obj/math_discrete.o
@echo '... Making bin/x86_64_Linux_debug_gcc_/libmath_math.a'
@mkdir -p /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_
@ar r /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/libmath_math.a /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/obj/math_discrete.o /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/obj/math_geometry.o
/net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/obj/math_geometry.o: /net/hc295/gmravi/home/fastlib/fastlib/math/math.h /net/hc295/gmravi/home/fastlib/fastlib/math/discrete.h /net/hc295/gmravi/home/fastlib/fastlib/math/geometry.cc /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_COMMON_gcc_COMMON/libcol_col.h /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_COMMON_gcc_COMMON/libbase_base.h /net/hc295/gmravi/home/fastlib/fastlib/math/kernel.h /net/hc295/gmravi/home/fastlib/fastlib/math/geometry.h
@echo '... Making bin/x86_64_Linux_debug_gcc_/obj/math_geometry.o'
@mkdir -p /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/obj
@cd /net/hc295/gmravi/home/fastlib/fastlib/math && g++ -Wall -Woverloaded-virtual -fno-exceptions -Wparentheses -fno-exceptions -I/net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_COMMON_gcc_COMMON -I/net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_COMMON_gcc_ -I/net/hc295/gmravi/home/fastlib/fastlib -g3 -DDEBUG -O0 -c geometry.cc -o /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/obj/math_geometry.o
/net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/obj/math_discrete.o: /net/hc295/gmravi/home/fastlib/fastlib/math/math.h /net/hc295/gmravi/home/fastlib/fastlib/math/discrete.h /net/hc295/gmravi/home/fastlib/fastlib/math/discrete.cc /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_COMMON_gcc_COMMON/libcol_col.h /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_COMMON_gcc_COMMON/libbase_base.h /net/hc295/gmravi/home/fastlib/fastlib/math/kernel.h /net/hc295/gmravi/home/fastlib/fastlib/math/geometry.h
@echo '... Making bin/x86_64_Linux_debug_gcc_/obj/math_discrete.o'
@mkdir -p /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/obj
@cd /net/hc295/gmravi/home/fastlib/fastlib/math && g++ -Wall -Woverloaded-virtual -fno-exceptions -Wparentheses -fno-exceptions -I/net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_COMMON_gcc_COMMON -I/net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_COMMON_gcc_ -I/net/hc295/gmravi/home/fastlib/fastlib -g3 -DDEBUG -O0 -c discrete.cc -o /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/obj/math_discrete.o
pseudo_43: /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/libcol_col.a /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/libdata_data.a /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/libbase_base.a /net/hc295/gmravi/home/fastlib/fastlib/bin_keep/x86_64_Linux_COMMON_gcc_COMMON/libblaspack.a /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/libla_la.a /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/libfile_file.a
@echo '*** Done with /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/libcol_col.a /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/libdata_data.a /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/libbase_base.a /net/hc295/gmravi/home/fastlib/fastlib/bin_keep/x86_64_Linux_COMMON_gcc_COMMON/libblaspack.a /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/libla_la.a /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/libfile_file.a'
/net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/libdata_data.a: /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/obj/data_dataset.o
@echo '... Making bin/x86_64_Linux_debug_gcc_/libdata_data.a'
@mkdir -p /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_
@ar r /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/libdata_data.a /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/obj/data_dataset.o
/net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/obj/data_dataset.o: /net/hc295/gmravi/home/fastlib/fastlib/data/dataset.h /net/hc295/gmravi/home/fastlib/fastlib/data/crossvalidation.h /net/hc295/gmravi/home/fastlib/fastlib/data/dataset.cc /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_COMMON_gcc_COMMON/libfile_file.h /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_COMMON_gcc_COMMON/libbase_base.h /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_COMMON_gcc_COMMON/libla_la.h
@echo '... Making bin/x86_64_Linux_debug_gcc_/obj/data_dataset.o'
@mkdir -p /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/obj
@cd /net/hc295/gmravi/home/fastlib/fastlib/data && g++ -Wall -Woverloaded-virtual -fno-exceptions -Wparentheses -fno-exceptions -I/net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_COMMON_gcc_COMMON -I/net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_COMMON_gcc_ -I/net/hc295/gmravi/home/fastlib/fastlib -g3 -DDEBUG -O0 -c dataset.cc -o /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/obj/data_dataset.o
/net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_COMMON_gcc_COMMON/libla_la.h: /net/hc295/gmravi/home/fastlib/fastlib/la/blas.h /net/hc295/gmravi/home/fastlib/fastlib/la/matrix.h /net/hc295/gmravi/home/fastlib/fastlib/la/uselapack.h /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_COMMON_gcc_COMMON/libcol_col.h /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_COMMON_gcc_COMMON/libbase_base.h /net/hc295/gmravi/home/fastlib/fastlib/la/clapack.h /net/hc295/gmravi/home/fastlib/fastlib/la/la.h
@echo '... Making bin/x86_64_Linux_COMMON_gcc_COMMON/libla_la.h'
@mkdir -p /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_COMMON_gcc_COMMON
@touch /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_COMMON_gcc_COMMON/libla_la.h
/net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_COMMON_gcc_COMMON/libfile_file.h: /net/hc295/gmravi/home/fastlib/fastlib/file/textfile.h /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_COMMON_gcc_COMMON/libcol_col.h /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_COMMON_gcc_COMMON/libbase_base.h
@echo '... Making bin/x86_64_Linux_COMMON_gcc_COMMON/libfile_file.h'
@mkdir -p /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_COMMON_gcc_COMMON
@touch /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_COMMON_gcc_COMMON/libfile_file.h
pseudo_51: /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/libcol_col.a /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/libfile_file.a /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/libbase_base.a
@echo '*** Done with /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/libcol_col.a /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/libfile_file.a /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/libbase_base.a'
pseudo_53: /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/libcol_col.a /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/libfile_file.a /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/libbase_base.a
@echo '*** Done with /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/libcol_col.a /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/libfile_file.a /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/libbase_base.a'
/net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/libfile_file.a: /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/obj/file_textfile.o
@echo '... Making bin/x86_64_Linux_debug_gcc_/libfile_file.a'
@mkdir -p /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_
@ar r /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/libfile_file.a /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/obj/file_textfile.o
/net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/obj/file_textfile.o: /net/hc295/gmravi/home/fastlib/fastlib/file/textfile.cc /net/hc295/gmravi/home/fastlib/fastlib/file/textfile.h /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_COMMON_gcc_COMMON/libcol_col.h /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_COMMON_gcc_COMMON/libbase_base.h
@echo '... Making bin/x86_64_Linux_debug_gcc_/obj/file_textfile.o'
@mkdir -p /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/obj
@cd /net/hc295/gmravi/home/fastlib/fastlib/file && g++ -Wall -Woverloaded-virtual -fno-exceptions -Wparentheses -fno-exceptions -I/net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_COMMON_gcc_COMMON -I/net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_COMMON_gcc_ -I/net/hc295/gmravi/home/fastlib/fastlib -g3 -DDEBUG -O0 -c textfile.cc -o /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/obj/file_textfile.o
pseudo_58: /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/libfx_fx.a /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/libbase_base.a
@echo '*** Done with /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/libfx_fx.a /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/libbase_base.a'
/net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/libfx_fx.a: /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/obj/fx_fx.o /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/obj/fx_timer.o /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/obj/fx_datastore.o
@echo '... Making bin/x86_64_Linux_debug_gcc_/libfx_fx.a'
@mkdir -p /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_
@ar r /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/libfx_fx.a /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/obj/fx_datastore.o /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/obj/fx_timer.o /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/obj/fx_fx.o
/net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/obj/fx_fx.o: /net/hc295/gmravi/home/fastlib/fastlib/fx/fx.h /net/hc295/gmravi/home/fastlib/fastlib/fx/fx.c /net/hc295/gmravi/home/fastlib/fastlib/fx/timer.h /net/hc295/gmravi/home/fastlib/fastlib/fx/datastore.h /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_COMMON_gcc_COMMON/libbase_base.h
@echo '... Making bin/x86_64_Linux_debug_gcc_/obj/fx_fx.o'
@mkdir -p /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/obj
@cd /net/hc295/gmravi/home/fastlib/fastlib/fx && gcc -I/net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_COMMON_gcc_COMMON -I/net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_COMMON_gcc_ -I/net/hc295/gmravi/home/fastlib/fastlib -g3 -DDEBUG -O0 -c fx.c -o /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/obj/fx_fx.o -Wall
/net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/obj/fx_timer.o: /net/hc295/gmravi/home/fastlib/fastlib/fx/fx.h /net/hc295/gmravi/home/fastlib/fastlib/fx/timer.h /net/hc295/gmravi/home/fastlib/fastlib/fx/datastore.h /net/hc295/gmravi/home/fastlib/fastlib/fx/timer.c /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_COMMON_gcc_COMMON/libbase_base.h
@echo '... Making bin/x86_64_Linux_debug_gcc_/obj/fx_timer.o'
@mkdir -p /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/obj
@cd /net/hc295/gmravi/home/fastlib/fastlib/fx && gcc -I/net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_COMMON_gcc_COMMON -I/net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_COMMON_gcc_ -I/net/hc295/gmravi/home/fastlib/fastlib -g3 -DDEBUG -O0 -c timer.c -o /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/obj/fx_timer.o -Wall
/net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/obj/fx_datastore.o: /net/hc295/gmravi/home/fastlib/fastlib/fx/fx.h /net/hc295/gmravi/home/fastlib/fastlib/fx/timer.h /net/hc295/gmravi/home/fastlib/fastlib/fx/datastore.h /net/hc295/gmravi/home/fastlib/fastlib/fx/datastore.c /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_COMMON_gcc_COMMON/libbase_base.h
@echo '... Making bin/x86_64_Linux_debug_gcc_/obj/fx_datastore.o'
@mkdir -p /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/obj
@cd /net/hc295/gmravi/home/fastlib/fastlib/fx && gcc -I/net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_COMMON_gcc_COMMON -I/net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_COMMON_gcc_ -I/net/hc295/gmravi/home/fastlib/fastlib -g3 -DDEBUG -O0 -c datastore.c -o /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/obj/fx_datastore.o -Wall
pseudo_69: /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/libcol_col.a /net/hc295/gmravi/home/fastlib/fastlib/bin_keep/x86_64_Linux_COMMON_gcc_COMMON/libblaspack.a /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/libla_la.a /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/libbase_base.a
@echo '*** Done with /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/libcol_col.a /net/hc295/gmravi/home/fastlib/fastlib/bin_keep/x86_64_Linux_COMMON_gcc_COMMON/libblaspack.a /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/libla_la.a /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/libbase_base.a'
/net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/libla_la.a: /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/obj/la_uselapack.o
@echo '... Making bin/x86_64_Linux_debug_gcc_/libla_la.a'
@mkdir -p /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_
@ar r /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/libla_la.a /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/obj/la_uselapack.o
/net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/obj/la_uselapack.o: /net/hc295/gmravi/home/fastlib/fastlib/la/blas.h /net/hc295/gmravi/home/fastlib/fastlib/la/matrix.h /net/hc295/gmravi/home/fastlib/fastlib/la/uselapack.cc /net/hc295/gmravi/home/fastlib/fastlib/la/uselapack.h /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_COMMON_gcc_COMMON/libcol_col.h /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_COMMON_gcc_COMMON/libbase_base.h /net/hc295/gmravi/home/fastlib/fastlib/la/clapack.h /net/hc295/gmravi/home/fastlib/fastlib/la/la.h
@echo '... Making bin/x86_64_Linux_debug_gcc_/obj/la_uselapack.o'
@mkdir -p /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/obj
@cd /net/hc295/gmravi/home/fastlib/fastlib/la && g++ -Wall -Woverloaded-virtual -fno-exceptions -Wparentheses -fno-exceptions -I/net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_COMMON_gcc_COMMON -I/net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_COMMON_gcc_ -I/net/hc295/gmravi/home/fastlib/fastlib -g3 -DDEBUG -O0 -c uselapack.cc -o /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/obj/la_uselapack.o
/net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_COMMON_gcc_COMMON/libcol_col.h: /net/hc295/gmravi/home/fastlib/fastlib/col/fastalloc.h /net/hc295/gmravi/home/fastlib/fastlib/col/queue.h /net/hc295/gmravi/home/fastlib/fastlib/col/heap.h /net/hc295/gmravi/home/fastlib/fastlib/col/intmap.h /net/hc295/gmravi/home/fastlib/fastlib/col/rangeset.h /net/hc295/gmravi/home/fastlib/fastlib/col/arraylist.h /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_COMMON_gcc_COMMON/libbase_base.h /net/hc295/gmravi/home/fastlib/fastlib/col/string.h
@echo '... Making bin/x86_64_Linux_COMMON_gcc_COMMON/libcol_col.h'
@mkdir -p /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_COMMON_gcc_COMMON
@touch /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_COMMON_gcc_COMMON/libcol_col.h
/net/hc295/gmravi/home/fastlib/fastlib/bin_keep/x86_64_Linux_COMMON_gcc_COMMON/libblaspack.a: /net/hc295/gmravi/home/fastlib/fastlib/la/blaspack.tgz
@echo '... Making /net/hc295/gmravi/home/fastlib/fastlib/bin_keep/x86_64_Linux_COMMON_gcc_COMMON/libblaspack.a'
@mkdir -p /net/hc295/gmravi/home/fastlib/fastlib/bin_keep/x86_64_Linux_COMMON_gcc_COMMON
@echo '... Extracting FORTRAN files...'
@mkdir -p /net/hc295/gmravi/home/fastlib/fastlib/bin_keep/x86_64_Linux_COMMON_gcc_COMMON/libblaspack_workspace
@cd /net/hc295/gmravi/home/fastlib/fastlib/bin_keep/x86_64_Linux_COMMON_gcc_COMMON/libblaspack_workspace && tar -xzf /net/hc295/gmravi/home/fastlib/fastlib/la/blaspack.tgz
@echo '*** Compiling LAPACK with BLAS reference implementation.'
@echo '!!! LAPACK WARNING: For better performance, install ATLAS or Intel MKL.'
@echo '... Our compilation differs slightly from regular LAPACK/BLAS:'
@echo '... NOTE 2: We omit complex-number routines (halves compile time).'
@echo '... NOTE 3: We require case sensitivity for LAPACK/BLAS string parameters.'
@echo '... This may take several minutes (about 800 FORTRAN files).'
@cd /net/hc295/gmravi/home/fastlib/fastlib/bin_keep/x86_64_Linux_COMMON_gcc_COMMON/libblaspack_workspace && g77 -O2 -c src/*.f
@echo '... Almost done with LAPACK/BLAS...'
@cd /net/hc295/gmravi/home/fastlib/fastlib/bin_keep/x86_64_Linux_COMMON_gcc_COMMON/libblaspack_workspace && g77 -O0 -c src/dlamch.f src/slamch.f
@cd /net/hc295/gmravi/home/fastlib/fastlib/bin_keep/x86_64_Linux_COMMON_gcc_COMMON/libblaspack_workspace && ar r /net/hc295/gmravi/home/fastlib/fastlib/bin_keep/x86_64_Linux_COMMON_gcc_COMMON/libblaspack.a *.o
@echo '... Created archive, cleaning up.'
@rm -rf /net/hc295/gmravi/home/fastlib/fastlib/bin_keep/x86_64_Linux_COMMON_gcc_COMMON/libblaspack_workspace
@echo '*** Done with LAPACK and BLAS!'
/net/hc295/gmravi/home/fastlib/fastlib/la/blaspack.tgz:
@echo '... Making /net/hc295/gmravi/home/fastlib/fastlib/la/blaspack.tgz'
@echo 'Downloading the file using curl...'
@cd /net/hc295/gmravi/home/fastlib/fastlib/la && curl -L -o blaspack.tgz "http://www.cc.gatech.edu/~garryb/fastlib/blaspack.tgz"
pseudo_81: /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/libcol_col.a /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/libbase_base.a
@echo '*** Done with /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/libcol_col.a /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/libbase_base.a'
/net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/libcol_col.a: /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/obj/col_col.o
@echo '... Making bin/x86_64_Linux_debug_gcc_/libcol_col.a'
@mkdir -p /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_
@ar r /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/libcol_col.a /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/obj/col_col.o
/net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/obj/col_col.o: /net/hc295/gmravi/home/fastlib/fastlib/col/fastalloc.h /net/hc295/gmravi/home/fastlib/fastlib/col/queue.h /net/hc295/gmravi/home/fastlib/fastlib/col/heap.h /net/hc295/gmravi/home/fastlib/fastlib/col/intmap.h /net/hc295/gmravi/home/fastlib/fastlib/col/rangeset.h /net/hc295/gmravi/home/fastlib/fastlib/col/arraylist.h /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_COMMON_gcc_COMMON/libbase_base.h /net/hc295/gmravi/home/fastlib/fastlib/col/col.cc /net/hc295/gmravi/home/fastlib/fastlib/col/string.h
@echo '... Making bin/x86_64_Linux_debug_gcc_/obj/col_col.o'
@mkdir -p /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/obj
@cd /net/hc295/gmravi/home/fastlib/fastlib/col && g++ -Wall -Woverloaded-virtual -fno-exceptions -Wparentheses -fno-exceptions -I/net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_COMMON_gcc_COMMON -I/net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_COMMON_gcc_ -I/net/hc295/gmravi/home/fastlib/fastlib -g3 -DDEBUG -O0 -c col.cc -o /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/obj/col_col.o
/net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_COMMON_gcc_COMMON/libbase_base.h: /net/hc295/gmravi/home/fastlib/fastlib/base/ccmem.h /net/hc295/gmravi/home/fastlib/fastlib/base/scale.h /net/hc295/gmravi/home/fastlib/fastlib/base/cc.h /net/hc295/gmravi/home/fastlib/fastlib/base/common.h /net/hc295/gmravi/home/fastlib/fastlib/base/otrav_impl.h /net/hc295/gmravi/home/fastlib/fastlib/base/test.h /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_COMMON_gcc_COMMON/base/basic_types.h /net/hc295/gmravi/home/fastlib/fastlib/base/compiler_impl.h /net/hc295/gmravi/home/fastlib/fastlib/base/compiler.h /net/hc295/gmravi/home/fastlib/fastlib/base/fortran.h /net/hc295/gmravi/home/fastlib/fastlib/base/debug.h /net/hc295/gmravi/home/fastlib/fastlib/base/otrav.h
@echo '... Making bin/x86_64_Linux_COMMON_gcc_COMMON/libbase_base.h'
@mkdir -p /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_COMMON_gcc_COMMON
@touch /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_COMMON_gcc_COMMON/libbase_base.h
pseudo_93: /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/libbase_base.a
@echo '*** Done with /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/libbase_base.a'
/net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/libbase_base.a: /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/obj/base_common.o /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/obj/base_cc.o /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/obj/base_otrav.o /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/obj/base_ccmem.o
@echo '... Making bin/x86_64_Linux_debug_gcc_/libbase_base.a'
@mkdir -p /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_
@ar r /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/libbase_base.a /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/obj/base_common.o /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/obj/base_cc.o /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/obj/base_ccmem.o /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/obj/base_otrav.o
/net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/obj/base_otrav.o: /net/hc295/gmravi/home/fastlib/fastlib/base/ccmem.h /net/hc295/gmravi/home/fastlib/fastlib/base/otrav.cc /net/hc295/gmravi/home/fastlib/fastlib/base/scale.h /net/hc295/gmravi/home/fastlib/fastlib/base/cc.h /net/hc295/gmravi/home/fastlib/fastlib/base/common.h /net/hc295/gmravi/home/fastlib/fastlib/base/otrav_impl.h /net/hc295/gmravi/home/fastlib/fastlib/base/test.h /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_COMMON_gcc_COMMON/base/basic_types.h /net/hc295/gmravi/home/fastlib/fastlib/base/compiler_impl.h /net/hc295/gmravi/home/fastlib/fastlib/base/compiler.h /net/hc295/gmravi/home/fastlib/fastlib/base/fortran.h /net/hc295/gmravi/home/fastlib/fastlib/base/debug.h /net/hc295/gmravi/home/fastlib/fastlib/base/otrav.h
@echo '... Making bin/x86_64_Linux_debug_gcc_/obj/base_otrav.o'
@mkdir -p /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/obj
@cd /net/hc295/gmravi/home/fastlib/fastlib/base && g++ -Wall -Woverloaded-virtual -fno-exceptions -Wparentheses -fno-exceptions -I/net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_COMMON_gcc_COMMON -I/net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_COMMON_gcc_ -I/net/hc295/gmravi/home/fastlib/fastlib -g3 -DDEBUG -O0 -c otrav.cc -o /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/obj/base_otrav.o
/net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/obj/base_ccmem.o: /net/hc295/gmravi/home/fastlib/fastlib/base/ccmem.h /net/hc295/gmravi/home/fastlib/fastlib/base/scale.h /net/hc295/gmravi/home/fastlib/fastlib/base/cc.h /net/hc295/gmravi/home/fastlib/fastlib/base/common.h /net/hc295/gmravi/home/fastlib/fastlib/base/otrav_impl.h /net/hc295/gmravi/home/fastlib/fastlib/base/test.h /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_COMMON_gcc_COMMON/base/basic_types.h /net/hc295/gmravi/home/fastlib/fastlib/base/compiler_impl.h /net/hc295/gmravi/home/fastlib/fastlib/base/compiler.h /net/hc295/gmravi/home/fastlib/fastlib/base/fortran.h /net/hc295/gmravi/home/fastlib/fastlib/base/debug.h /net/hc295/gmravi/home/fastlib/fastlib/base/ccmem.cc /net/hc295/gmravi/home/fastlib/fastlib/base/otrav.h
@echo '... Making bin/x86_64_Linux_debug_gcc_/obj/base_ccmem.o'
@mkdir -p /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/obj
@cd /net/hc295/gmravi/home/fastlib/fastlib/base && g++ -Wall -Woverloaded-virtual -fno-exceptions -Wparentheses -fno-exceptions -I/net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_COMMON_gcc_COMMON -I/net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_COMMON_gcc_ -I/net/hc295/gmravi/home/fastlib/fastlib -g3 -DDEBUG -O0 -c ccmem.cc -o /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/obj/base_ccmem.o
/net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/obj/base_cc.o: /net/hc295/gmravi/home/fastlib/fastlib/base/ccmem.h /net/hc295/gmravi/home/fastlib/fastlib/base/scale.h /net/hc295/gmravi/home/fastlib/fastlib/base/cc.h /net/hc295/gmravi/home/fastlib/fastlib/base/common.h /net/hc295/gmravi/home/fastlib/fastlib/base/otrav_impl.h /net/hc295/gmravi/home/fastlib/fastlib/base/test.h /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_COMMON_gcc_COMMON/base/basic_types.h /net/hc295/gmravi/home/fastlib/fastlib/base/cc.cc /net/hc295/gmravi/home/fastlib/fastlib/base/compiler_impl.h /net/hc295/gmravi/home/fastlib/fastlib/base/compiler.h /net/hc295/gmravi/home/fastlib/fastlib/base/fortran.h /net/hc295/gmravi/home/fastlib/fastlib/base/debug.h /net/hc295/gmravi/home/fastlib/fastlib/base/otrav.h
@echo '... Making bin/x86_64_Linux_debug_gcc_/obj/base_cc.o'
@mkdir -p /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/obj
@cd /net/hc295/gmravi/home/fastlib/fastlib/base && g++ -Wall -Woverloaded-virtual -fno-exceptions -Wparentheses -fno-exceptions -I/net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_COMMON_gcc_COMMON -I/net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_COMMON_gcc_ -I/net/hc295/gmravi/home/fastlib/fastlib -g3 -DDEBUG -O0 -c cc.cc -o /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/obj/base_cc.o
/net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/obj/base_common.o: /net/hc295/gmravi/home/fastlib/fastlib/base/ccmem.h /net/hc295/gmravi/home/fastlib/fastlib/base/scale.h /net/hc295/gmravi/home/fastlib/fastlib/base/cc.h /net/hc295/gmravi/home/fastlib/fastlib/base/common.h /net/hc295/gmravi/home/fastlib/fastlib/base/otrav_impl.h /net/hc295/gmravi/home/fastlib/fastlib/base/test.h /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_COMMON_gcc_COMMON/base/basic_types.h /net/hc295/gmravi/home/fastlib/fastlib/base/common.c /net/hc295/gmravi/home/fastlib/fastlib/base/compiler_impl.h /net/hc295/gmravi/home/fastlib/fastlib/base/compiler.h /net/hc295/gmravi/home/fastlib/fastlib/base/fortran.h /net/hc295/gmravi/home/fastlib/fastlib/base/debug.h /net/hc295/gmravi/home/fastlib/fastlib/base/otrav.h
@echo '... Making bin/x86_64_Linux_debug_gcc_/obj/base_common.o'
@mkdir -p /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/obj
@cd /net/hc295/gmravi/home/fastlib/fastlib/base && gcc -I/net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_COMMON_gcc_COMMON -I/net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_COMMON_gcc_ -I/net/hc295/gmravi/home/fastlib/fastlib -g3 -DDEBUG -O0 -c common.c -o /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/obj/base_common.o -Wall
/net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_COMMON_gcc_COMMON/base/basic_types.h: /net/hc295/gmravi/home/fastlib/fastlib/script/config.py /net/hc295/gmravi/home/fastlib/fastlib/base/config/template_types.c
@echo '... Making bin/x86_64_Linux_COMMON_gcc_COMMON/base/basic_types.h'
@/net/hc295/gmravi/home/fastlib/fastlib/script/config.py --genfiles_dir=/net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_COMMON_gcc_COMMON --source_dir=/net/hc295/gmravi/home/fastlib/fastlib
@mkdir -p /net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_COMMON_gcc_COMMON/base
clean:
@echo 'Removing the bin/ directory.'
rm -rf /net/hc295/gmravi/home/fastlib/fastlib/bin/*
fastlib/u/gmravi/regression/build.py
+18
View File
@@ -0,0 +1,18 @@
binrule(
name = "main", # the executable name
sources = ["main.cc"], #
headers = ["regression.h"], # no extra headers
deplibs = ["fastlib:fastlib_int"]
)
# to build:
# 1. make sure have environment variables set up:
# $ source /full/path/to/fastlib/script/fl-env /full/path/to/fastlib
# (you might want to put this in bashrc)
# 2. fl-build main
# - this automatically will assume --mode=check, the default
# - type fl-build --help for help
# 3. ./main
# - to build same target again, type: make
# - to force recompilation, type: make clean
fastlib/u/gmravi/regression/cmdline_params
+18
View File
@@ -0,0 +1,18 @@
These are the command line parameters which the code regression.h uses
do_naive=true/false
method=fast
fast_kde_output: file name to which the output of fast kde should be written to
naive_kde_output:file name to which the navie kde should outptut its results
kernel: the kernel u want to use
data: the file name from which the data is read
query: the file name from which the query is read
fastlib/u/gmravi/regression/cmdline_params.txt
+31
View File
@@ -0,0 +1,31 @@
These are the command line parameters which the code regression.h uses
do_naive=true/false
method=fast
fast_kde_output: file name to which the output of fast kde should be written to
naive_kde_output:file name to which the navie kde should outptut its results
kernel: the kernel u want to use
data: the file name from which the data is read
query: the file name from which the query is read
bandwidth: The bandwidth of the density estimation
weights: this file is the weights assigned to the dofferent points
the command should be:
./main --method=fast --do_naive --fast_kde_output=fast_kde_output.txt --naive_kde_output=naive_kde_output.txt --kernel=gaussian --data=query_file.csv --query=query_file.csv --scaling=none --bandwidth=0.1
./main --method=fast --do_naive --fast_kde_output=fast_kde_output.txt --naive_kde_output=naive_kde_output.txt --kernel=gaussian --data=colors50k.ds --query=colors50k.ds --scaling=none --bandwidth=0.1
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
+20
View File
@@ -0,0 +1,20 @@
0.012931
0.004784
0.008621
0.017845
0.012751
0.018197
0.021199
0.034751
0.012699
0.031747
0.033863
0.051949
0.012699
-0.004535
-0.008466
-0.069266
0.012983
0.013763
0.008707
0.029319
fastlib/u/gmravi/regression/fft_kde.h
+883
View File
@@ -0,0 +1,883 @@
#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
@@ -0,0 +1,346 @@
//-------------------------------------------------------------------
// 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
@@ -0,0 +1,139 @@
//-------------------------------------------------------------------
// 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
@@ -0,0 +1,193 @@
#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
@@ -0,0 +1,317 @@
//-------------------------------------------------------------------
// 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
@@ -0,0 +1,118 @@
//-------------------------------------------------------------------
// 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
@@ -0,0 +1,973 @@
#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
@@ -0,0 +1,3 @@
function k=kg(x)
k=exp(-x .* x ./2)/sqrt(2*pi);
end
fastlib/u/gmravi/regression/lscvscore.m
+16
View File
@@ -0,0 +1,16 @@
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/main
Symlink
+1
View File
@@ -0,0 +1 @@
/net/hc295/gmravi/home/fastlib/fastlib/bin/x86_64_Linux_debug_gcc_/u/gmravi/regression/main
fastlib/u/gmravi/regression/main.cc
+103
View File
@@ -0,0 +1,103 @@
#include "fastlib/fastlib_int.h"
#include "regression.h"
#define MAXDOUBLE 32768.0
int
main (int argc, char *argv[])
{
fx_init (argc, argv);
const char *algorithm = fx_param_str_req (NULL, "method");
bool do_naive = fx_param_exists (NULL, "do_naive");
Matrix query_dataset;
Matrix reference_dataset;
Vector rset_weights;
ArrayList<Vector> fast_kde_results;
FastKde <GaussianKernel> fast_kde;
if (!strcmp (algorithm, "fast"))
{
if (!strcmp (fx_param_str (NULL, "kernel", "gaussian"), "gaussian"))
{
// printf("will do fast kde with gaussian kernel....\n");
ArrayList<Vector> fast_kde_results;
printf ("O(p^D) expansion KDE\n");
fast_kde.Init ();
fast_kde.Compute (fx_param_double (NULL, "tau", 0.1));
//default value is 0.1
printf("all computations done......");
if (fx_param_exists (NULL, "fast_kde_output")) {
fast_kde.PrintDebug ();
}
printf("Size of reference datset is %d\n",fast_kde.get_reference_dataset().n_cols());
fast_kde_results.Init(fast_kde.get_query_dataset().n_cols());
for(int i=0;i<fast_kde.get_query_dataset().n_cols();i++)
fast_kde_results[i].Alias(fast_kde.get_density_estimates(i));
query_dataset.Copy (fast_kde.get_query_dataset());
reference_dataset.Copy (fast_kde.get_reference_dataset());
rset_weights.Copy(fast_kde.get_reference_weights());
printf("all values copied for naive kde...\n");
}
if (do_naive)
{
//Vector fast_kde_results;
NaiveKde <GaussianKernel> naive_kde;
naive_kde.Init (query_dataset, reference_dataset);
naive_kde.Compute ();
if (fx_param_exists (NULL, "naive_kde_output"))
{
naive_kde.PrintDebug ();
}
printf("Naive KDE done...\n");
printf("Fast kde results length is %d\n",fast_kde_results.size());
naive_kde.ComputeMaximumRelativeError(fast_kde_results);
}
}
//Do fastkde wih epanechnikov kernel
else
{
/* if (!strcmp (fx_param_str (NULL, "kernel", "epan"), "epan"))
{
FastKde < EpanKernel > fast_kde;
fast_kde.Init ();
fast_kde.Compute (fx_param_double (NULL, "tau", 0.1));
if (fx_param_exists (NULL, "fast_kde_output"))
{
fast_kde.PrintDebug ();
}
Vector fast_kde_results;
fast_kde_results.Alias (fast_kde.get_density_estimates ());
if (do_naive)
{
NaiveKde < EpanKernel > naive_kde;
naive_kde.Init (fast_kde.get_query_dataset (),
fast_kde.get_reference_dataset ());
naive_kde.Compute ();
if (fx_param_exists (NULL, "naive_kde_output"))
{
naive_kde.PrintDebug ();
}
naive_kde.ComputeMaximumRelativeError (fast_kde_results);
}
}*/
printf("do nothing..............\n");
fx_done ();
return 0;
}
}
fastlib/u/gmravi/regression/naive_kde_output.txt
+20
View File
@@ -0,0 +1,20 @@
0.002586
0.002679
0.002586
0.002617
0.002550
0.002548
0.002544
0.002548
0.002540
0.002540
0.002540
0.002540
0.002540
0.002540
0.002540
0.002540
0.002597
0.002569
0.002612
0.002580
fastlib/u/gmravi/regression/query_file.csv
+5
View File
@@ -0,0 +1,5 @@
1.0, 2.0, 3.0
4.0, 5.0, 6.0
7.0, 8.0, 9.0
-1.0,-2.0,-12.0
3.0,2.0,5.0
1 1.0 2.0 3.0
2 4.0 5.0 6.0
3 7.0 8.0 9.0
4 -1.0 -2.0 -12.0
5 3.0 2.0 5.0
fastlib/u/gmravi/regression/regression.h
+1191
View File
File diff suppressed because it is too large Load Diff