diff --git a/fastlib/u/gmravi/regression/main_regression_ll.cc b/fastlib/u/gmravi/regression/main_regression_ll.cc index 984fe377fc..9c30458181 100644 --- a/fastlib/u/gmravi/regression/main_regression_ll.cc +++ b/fastlib/u/gmravi/regression/main_regression_ll.cc @@ -11,6 +11,7 @@ int main(int argc, char *argv[]){ Matrix q_matrix; Matrix r_matrix; Vector rset_weights; + Vector true_regression_values; //Reading parameters and loading data @@ -18,7 +19,8 @@ int main(int argc, char *argv[]){ fx_submodule(NULL, "regression", "regression_module"); // The reference data file is a required parameter. - const char* reference_file_name = fx_param_str_req(regression_module, "data"); + const char* reference_file_name = + fx_param_str_req(regression_module, "data"); // The query data file defaults to the references. const char* query_file_name = @@ -49,8 +51,8 @@ int main(int argc, char *argv[]){ //Get the bandwidth for kernel calculations and the tolerance limit. //Both default to 0.2 - double bandwidth=fx_param_double(regression_module,"bandwidth",0.15); - double tau=fx_param_double(regression_module,"tau",0.15); + double bandwidth=fx_param_double(regression_module,"bandwidth",0.125); + double tau=fx_param_double(regression_module,"tau",0.0025); printf("tau is %f\n",tau); //Get the weights for the reference set @@ -65,7 +67,7 @@ int main(int argc, char *argv[]){ Dataset ref_weights; ref_weights.InitFromFile (rwfname); - rset_weights.Copy (ref_weights.matrix ().GetColumnPtr (0), ref_weights.matrix ().n_rows ()); //Note rset_weights_ is a vector of weights + rset_weights.Copy (ref_weights.matrix ().GetColumnPtr (0), ref_weights.matrix ().n_rows ()); //Note rset_weights is a vector of weights } @@ -76,6 +78,30 @@ int main(int argc, char *argv[]){ } + //Also get the true regression values of all the reference points + const char *are=NULL; + if (fx_param_exists (NULL, "true_reg")){ + + //are is the filename having the regression estimates of the + are = fx_param_str (NULL, "true_reg", NULL); + } + + if (are != NULL){ + + Dataset true_reg; + true_reg.InitFromFile (are); + true_regression_values.Copy (true_reg.matrix ().GetColumnPtr (0), true_reg.matrix ().n_rows ()); + + } + + else{ + + true_regression_values.Init (r_matrix.n_cols ()); + true_regression_values.SetAll (1); + + }//Hence the true regression values have been stored + + //Get the length of the leaf. Defaulted to 2 index_t leaf_length=fx_param_int(regression_module,"leaf_length",3); @@ -88,23 +114,31 @@ int main(int argc, char *argv[]){ criteria=(char*)malloc(40*sizeof(char)); strcpy(criteria,fx_param_str(regression_module,"criteria","fnorm")); - fast_regression.Init(q_matrix, r_matrix, bandwidth, tau, leaf_length, rset_weights, criteria); + fast_regression.Init(q_matrix, r_matrix, bandwidth, tau, + leaf_length, rset_weights, criteria,true_regression_values); fast_regression.Compute(); + printf("FAST CALCULATIONS ALL DONE"); + + //These are the regression estimates provided by the fast + //algorithm + + // Vector fast_regression_estimate; + //fast_regression_estimate.Copy(fast_regression.get_regression_estimate()); //Lets do naive calculations too............ - - //Lets first declare an object of the naive type - printf("FAST CALCULATIONS ALL DONE"); + ArrayList old_from_new_r; old_from_new_r.Copy(fast_regression.get_old_from_new_r()); + + + //Get the fast regression estimates.We shall use it to compare the //accuracy - Vector fast_regression_estimate; - fast_regression_estimate.Copy(fast_regression.get_regression_estimate()); + ArrayList fast_b_twy_estimates; ArrayList fast_b_twb_estimates; @@ -113,23 +147,32 @@ int main(int argc, char *argv[]){ fast_b_twb_estimates.Init(q_matrix.n_cols()); for(index_t q=0;q naive; - naive.Init(q_matrix,r_matrix,old_from_new_r,bandwidth,rset_weights); + naive.Init(q_matrix,r_matrix, + old_from_new_r,bandwidth,rset_weights); + + printf("started naive computations..\n"); naive.Compute(); + printf("completed naive calculations..\n"); naive.ComputeMaximumRelativeError - (fast_b_twy_estimates,fast_b_twb_estimates,criteria); - naive.CompareFastWithNaive(fast_regression_estimate); + (fast_b_twy_estimates,fast_b_twb_estimates,criteria); + //naive.CompareFastWithNaive(fast_regression_estimate); - //Will need to verify if this is fine to do... - - fx_done(); } } diff --git a/fastlib/u/gmravi/regression/regression_ll1.h b/fastlib/u/gmravi/regression/regression_ll1.h index e261cd9975..ac93b669a0 100644 --- a/fastlib/u/gmravi/regression/regression_ll1.h +++ b/fastlib/u/gmravi/regression/regression_ll1.h @@ -1,7 +1,7 @@ #ifndef REGRESSION_LL1_H #define REGRESSION_LL1_H #include "regression_ll2.h" -#include "pseudo_inverse.h" +//#include "pseudo_inverse.h" /** In this code we shalll evaluate (B^TWB)^-1 by first calculating B^TWB @@ -421,8 +421,8 @@ void FastRegression::UpdateBounds_(Tree *qnode, //Add dl_b_twy and du_b_twy to mass_l and mass_u of B^TWY - la::AddTo(dl_b_twy, &qnode->stat().b_twy_mass_l); - la::AddTo(du_b_twy, &qnode->stat().b_twy_mass_u); + la::AddTo(dl_b_twy, &(qnode->stat().b_twy_mass_l)); + la::AddTo(du_b_twy, &(qnode->stat().b_twy_mass_u)); if(!qnode->is_leaf()){ @@ -469,8 +469,8 @@ void FastRegression::UpdateBounds_(Tree *qnode, //Add dl_b_twb and du_b_twb to mass_l and mass_u of B^TWB - la::AddTo(dl_b_twb, &qnode->stat().b_twb_mass_l); - la::AddTo(du_b_twb, &qnode->stat().b_twb_mass_u); + la::AddTo(dl_b_twb, &(qnode->stat().b_twb_mass_l)); + la::AddTo(du_b_twb, &(qnode->stat().b_twb_mass_u)); if(!qnode->is_leaf()){ @@ -493,30 +493,30 @@ void FastRegression::UpdateBounds_(Tree *qnode, else{ if(flag==CHECK_FOR_PRUNE_B_TWY){ - //the flag is CHECK_FOR_PRUNE_B_TWY + //the flag is CHECK_FOR_PRUNE_B_TWY - //Add dl_b_twy and du_b_twy to mass_l and mass_u of B^TWY + //Add dl_b_twy and du_b_twy to mass_l and mass_u of B^TWY - la::AddTo(dl_b_twy, &qnode->stat().b_twy_mass_l); - la::AddTo(du_b_twy, &qnode->stat().b_twy_mass_u); + la::AddTo(dl_b_twy, &(qnode->stat().b_twy_mass_l)); + la::AddTo(du_b_twy, &(qnode->stat().b_twy_mass_u)); - if(!qnode->is_leaf()){ + if(!qnode->is_leaf()){ - // transmit these values to the children node + // transmit these values to the children node - la::AddTo(dl_b_twy,&(qnode->left()->stat().b_twy_owed_l)); - la::AddTo(du_b_twy,&(qnode->left()->stat().b_twy_owed_u)); + la::AddTo(dl_b_twy,&(qnode->left()->stat().b_twy_owed_l)); + la::AddTo(du_b_twy,&(qnode->left()->stat().b_twy_owed_u)); - la::AddTo(dl_b_twy,&(qnode->right()->stat().b_twy_owed_l)); - la::AddTo(du_b_twy,&(qnode->right()->stat().b_twy_owed_u)); - } + la::AddTo(dl_b_twy,&(qnode->right()->stat().b_twy_owed_l)); + la::AddTo(du_b_twy,&(qnode->right()->stat().b_twy_owed_u)); + } - else{ - - /* in case of leaf nodes add these values to more_l and more_u */ - la::AddTo(dl_b_twy,&(qnode->stat().b_twy_more_l)); - la::AddTo(du_b_twy,&(qnode->stat().b_twy_more_u)); - } + else{ + + /* in case of leaf nodes add these values to more_l and more_u */ + la::AddTo(dl_b_twy,&(qnode->stat().b_twy_more_l)); + la::AddTo(du_b_twy,&(qnode->stat().b_twy_more_u)); + } } } @@ -565,9 +565,10 @@ FastRegression::Prunable_(Tree *qnode, Tree *rnode, Matrix &dl_b_twb, /** this means both the quantities are prunable */ - //For a moment lets makc a change - - return PRUNE_BOTH; + + //FOR DEBUGGING + //return PRUNE_NOT; + return PRUNE_BOTH; } else{ if(flag1==1){ @@ -587,6 +588,7 @@ FastRegression::Prunable_(Tree *qnode, Tree *rnode, Matrix &dl_b_twb, else{ /** None of the quantities is prunable */ + printf("Couldnt prune anything..\n"); return PRUNE_NOT; } @@ -615,7 +617,8 @@ FastRegression::Prunable_(Tree *qnode, Tree *rnode, Matrix &dl_b_twb, } //flag==CHECK_FOR_PRUNE_B_TWY - index_t flag2=FastRegression::PrunableB_TWY_(qnode,rnode,dl_b_twy,du_b_twy); + index_t flag2=FastRegression:: + PrunableB_TWY_(qnode,rnode,dl_b_twy,du_b_twy); if(flag2==1){ //this means B_TWY is prunable @@ -624,6 +627,7 @@ FastRegression::Prunable_(Tree *qnode, Tree *rnode, Matrix &dl_b_twb, } else{ //it is not prunable. + printf("Could not prune..\n"); return PRUNE_NOT; } } @@ -713,9 +717,11 @@ void FastRegression::FRegression_(Tree *qnode, Tree *rnode, if(what_is_prunable==PRUNE_B_TWB){ + printf("BTWB pruned..\n"); + UpdateBoundsForPruningB_TWB_(qnode, dl_b_twb, du_b_twb); - MergeChildBounds_(qnode,flag); + // MergeChildBounds_(qnode,what_is_prunable); if((index_t)what_is_prunable==(index_t)flag){ //this means the job is done and we may return now return; @@ -731,9 +737,10 @@ void FastRegression::FRegression_(Tree *qnode, Tree *rnode, else{ if(what_is_prunable==PRUNE_B_TWY){ + printf("BTWY pruned..\n"); UpdateBoundsForPruningB_TWY_(qnode, dl_b_twy, du_b_twy); - MergeChildBounds_(qnode,flag); + //MergeChildBounds_(qnode,what_is_prunable); if((index_t)what_is_prunable==(index_t)flag){ //this means the job is done and we may return now @@ -748,10 +755,10 @@ void FastRegression::FRegression_(Tree *qnode, Tree *rnode, } else{ //both are prunable - + printf("Both pruned...\n"); UpdateBoundsForPruningB_TWB_(qnode, dl_b_twb, du_b_twb); UpdateBoundsForPruningB_TWY_(qnode, dl_b_twy, du_b_twy); - MergeChildBounds_(qnode,flag); + MergeChildBounds_(qnode,CHECK_FOR_PRUNE_BOTH); return; } @@ -832,9 +839,13 @@ void FastRegression::Compute(){ check_for_prune_t flag=CHECK_FOR_PRUNE_BOTH; fx_timer_start(NULL,"fast_timer"); + printf("entered recursive function..\n"); FRegression_(qroot_,rroot_,flag); + printf("Completed recursion..\n"); PostProcess_(qroot_); - ObtainRegressionEstimate_(); + printf("Will postprocess now..\n"); + //ObtainRegressionEstimate_(); + //CompareWithTrueValues_(); fx_timer_stop(NULL,"fast_timer"); //This will print the matrices B^TWB and B^TWY to an output file @@ -916,15 +927,22 @@ void FastRegression::SetUpperBounds_(Tree *node){ if(node->is_leaf()){ - node->stat().b_twy_mass_u.Copy(rroot_->stat().b_ty); + + node->stat().b_twy_mass_u.Init(rset_.n_rows()+1,1); + node->stat().b_twy_mass_u.CopyValues(rroot_->stat().b_ty); - node->stat().b_twb_mass_u.Copy(rroot_->stat().b_tb); + node->stat().b_twb_mass_u.Init(rset_.n_rows()+1,rset_.n_rows()+1); + node->stat().b_twb_mass_u.CopyValues(rroot_->stat().b_tb); } else { - node->stat().b_twy_mass_u .Copy(rroot_->stat().b_ty); - node->stat().b_twb_mass_u.Copy(rroot_->stat().b_tb); + node->stat().b_twy_mass_u.Init(rset_.n_rows()+1,1); + node->stat().b_twb_mass_u.Init(rset_.n_rows()+1,rset_.n_rows()+1); + + node->stat().b_twy_mass_u .CopyValues(rroot_->stat().b_ty); + node->stat().b_twb_mass_u.CopyValues(rroot_->stat().b_tb); + SetUpperBounds_(node->left()); SetUpperBounds_(node->right()); } @@ -936,18 +954,25 @@ void FastRegression::Init(Matrix &q_matrix, Matrix &r_matrix, double bandwidth, double tau, index_t leaf_length, - Vector &rset_weights, char *criteria){ + Vector &rset_weights, + char *criteria, + Vector &true_regression_values){ - //Set up value of tau,qset_,rset_ + + + //Set up value of tau,qset_,rset_, true_regression_values qset_.Alias(q_matrix); rset_.Alias(r_matrix); + true_regression_values_.Alias(true_regression_values); tau_=tau; index_t leaflen=leaf_length; rset_weights_.Alias(rset_weights); + /* Construct Query and Reference trees */ + fx_timer_start(NULL,"tree_create"); rroot_ = tree::MakeKdTreeMidpoint < Tree > (rset_, leaflen, &old_from_new_r_, &new_from_old_r_); diff --git a/fastlib/u/gmravi/regression/regression_ll2.h b/fastlib/u/gmravi/regression/regression_ll2.h index 0d8843f2ae..d5ea049e15 100644 --- a/fastlib/u/gmravi/regression/regression_ll2.h +++ b/fastlib/u/gmravi/regression/regression_ll2.h @@ -8,9 +8,9 @@ double FastRegression::Compute1NormLike_(Matrix &a){ //This function computes treats the natrix as a vector and computes //it's 1-norm - double value=0; - for(index_t col=0;col double FastRegression::SquaredFrobeniusNorm_(Matrix &a){ double sqd_frobenius_norm=0.0; - for(index_t col=0;col::PrunableB_TWY_(Tree *qnode, Tree *rnode, //THIS Means we are interested in the frobenius norm pruning................ - //Maximum allowed error is the sqaured frobenius norm of the matrix max_error + //Maximum error is the sqaured frobenius norm of the matrix max_error double squared_frobenius_norm_of_max_error=SquaredFrobeniusNorm_(max_error); double ratio_of_1norms= - Compute1NormLike_(rnode->stat().b_ty)/Compute1NormLike_(rroot_->stat().b_ty); + (double) Compute1NormLike_(rnode->stat().b_ty)/Compute1NormLike_(rroot_->stat().b_ty); //Now lets calculate the squared frobenius norm of new_mass_l double sqd_frobenius_norm_of_new_mass_l=SquaredFrobeniusNorm_(new_mass_l); double allowed_error= ratio_of_1norms*tau_*sqd_frobenius_norm_of_new_mass_l; - if(squared_frobenius_norm_of_max_error<=allowed_error) + if(squared_frobenius_norm_of_max_error:: PrunableB_TWB_(Tree *qnode, Tree *rnode, Matri //Now lets calculate the allowed error. - //allowed_error=tau*(rnode->stat().b_tb/rrot->stat().b_tb)(b_twb_mass_l+dl) + //allowed_error=tau*(rnode->stat().b_tb/rroot->stat().b_tb)(b_twb_mass_l+dl) //Lets define a matrix new_mass_l as shown below.This matrix takes //into account the addditional dl that will be added if pruning takes //place @@ -228,7 +237,7 @@ index_t FastRegression:: PrunableB_TWB_(Tree *qnode, Tree *rnode, Matri double alpha= tau_*(double)(rnode->stat().b_tb.get(row,col))/ (rroot_->stat().b_tb.get(row,col)); - if(rroot_->stat().b_tb.get(row,col)==0){ + if(rroot_->stat().b_tb.get(row,col)==0){ printf("btb.is not properly defined\n"); } @@ -249,6 +258,7 @@ index_t FastRegression:: PrunableB_TWB_(Tree *qnode, Tree *rnode, Matri } } //This means the matrix is compoenent wise prunable. hence return 1 + return 1; } //THIS MEANS WE ARE INTERESTED IN PRUNING BY FROBENIUS NORM @@ -256,14 +266,15 @@ index_t FastRegression:: PrunableB_TWB_(Tree *qnode, Tree *rnode, Matri double squared_frobenius_norm_of_max_error=SquaredFrobeniusNorm_(max_error); double ratio_of_1norms= - Compute1NormLike_(rnode->stat().b_tb)/Compute1NormLike_(rroot_->stat().b_tb); + (double)Compute1NormLike_(rnode->stat().b_tb)/Compute1NormLike_(rroot_->stat().b_tb); //Now lets calculate the squared frobenius norm of new_mass_l double squared_frobenius_norm_of_new_mass_l=SquaredFrobeniusNorm_(new_mass_l); double allowed_error=tau_*ratio_of_1norms*squared_frobenius_norm_of_new_mass_l; - if(squared_frobenius_norm_of_max_error > allowed_error){ + + if(squared_frobenius_norm_of_max_error >= allowed_error){ //then this matrix is NOT runable dl.SetAll(0); @@ -271,12 +282,19 @@ index_t FastRegression:: PrunableB_TWB_(Tree *qnode, Tree *rnode, Matri return 0; } - //This means that the quantity is prunable - return 1; + else{ + //This means that the quantity is prunable + printf("squared_frobenius_norm_of_max_error was %f\n",squared_frobenius_norm_of_max_error); + printf("allowed error is %f\n",allowed_error); + printf("kernel value lo is %f\n",kernel_value_range.lo); + printf("kernel value hi is %f\n",kernel_value_range.hi); + printf("\n\n"); + return 1; + } } - -/* The Update Boundws function is independent of the pruning criteria */ + +/* The Update Bounds function is independent of the pruning criteria */ template void FastRegression:: @@ -290,8 +308,8 @@ UpdateBoundsForPruningB_TWY_(Tree *qnode, Matrix &dl_b_twy, Matrix &du_b_twy){ //b_twy_mass_u <- b_twy_mass_u+du - la::AddTo(dl_b_twy, &qnode->stat().b_twy_mass_l); - la::AddTo(du_b_twy, &qnode->stat().b_twy_mass_u); + la::AddTo(dl_b_twy, &(qnode->stat().b_twy_mass_l)); + la::AddTo(du_b_twy, &(qnode->stat().b_twy_mass_u)); // for a leaf node, incorporate the lower and upper bound changes into @@ -396,31 +414,31 @@ MergeChildBoundsB_TWB_( FastRegression:: // parent //So lets find the componentwise minimum and maximum + if(pruning_criteria==CRITERIA_FOR_PRUNE_COMPONENT){ + Matrix max_children; + Matrix min_children; + max_children.Init(parent_stat.b_twb_mass_l.n_rows(), + parent_stat.b_twb_mass_l.n_cols()); + + min_children.Init(parent_stat.b_twb_mass_l.n_rows(), + parent_stat.b_twb_mass_l.n_cols()); - Matrix max_children; - Matrix min_children; - max_children.Init(parent_stat.b_twb_mass_l.n_rows(), - parent_stat.b_twb_mass_l.n_cols()); - - min_children.Init(parent_stat.b_twb_mass_l.n_rows(), - parent_stat.b_twb_mass_l.n_cols()); - - if(pruning_criteria==CRITERIA_FOR_PRUNE_COMPONENT){ - for(index_t col=0;colb_twb_mass_l.get(row,col) <= - right_stat->b_twb_mass_l.get(row,col)){ - - //left child has lesser mass_l value - - min_children.set(row,col, - left_stat->b_twb_mass_l.get(row,col)); - - max_children.set(row,col, - right_stat->b_twb_mass_l.get(row,col)); - } - + + for(index_t col=0;colb_twb_mass_l.get(row,col) <= + right_stat->b_twb_mass_l.get(row,col)){ + + //left child has lesser mass_l value + + min_children.set(row,col, + left_stat->b_twb_mass_l.get(row,col)); + + max_children.set(row,col, + right_stat->b_twb_mass_l.get(row,col)); + } + else{ //right child has lesser mass-l value @@ -432,9 +450,9 @@ MergeChildBoundsB_TWB_( FastRegression:: } - } + } - } + } //Now compare with parent... for(index_t col=0;col:: } } } - } + return; + } //THe pruning criteria is Frobenius norm pruning... @@ -538,6 +557,7 @@ MergeChildBoundsB_TWY_( FastRegression:: if(pruning_criteria==CRITERIA_FOR_PRUNE_COMPONENT){ + //This means we want to merge the bounds of b_twb //b_twb_mass_l_parent= //max(min(b_twy_mass_l,left_child,b_twy_mass_l_right_child), @@ -601,7 +621,9 @@ MergeChildBoundsB_TWY_( FastRegression:: } } } + return; } + //IF Pruning criteria is Frobenius Norm Pruning criteria..... else{ @@ -609,17 +631,20 @@ MergeChildBoundsB_TWY_( FastRegression:: // lower_of_parent=max (parent, min (children)) // upper_of_parent=min(parent, max(children)) - //lets first frind out the squared frobenius norms of the lower + //lets first find out the squared frobenius norms of the lower //bound masses of both the children - double sqd_frobenius_norm_of_left_child=SquaredFrobeniusNorm_(left_stat->b_twy_mass_l); - double sqd_frobenius_norm_of_right_child=SquaredFrobeniusNorm_(right_stat->b_twy_mass_l); + double sqd_frobenius_norm_of_left_child= + SquaredFrobeniusNorm_(left_stat->b_twy_mass_l); + double sqd_frobenius_norm_of_right_child= + SquaredFrobeniusNorm_(right_stat->b_twy_mass_l); //Now compare it with that of the parent - double sqd_frobenius_norm_of_parent=SquaredFrobeniusNorm_(parent_stat.b_twy_mass_l); + double sqd_frobenius_norm_of_parent= + SquaredFrobeniusNorm_(parent_stat.b_twy_mass_l); - if(sqd_frobenius_norm_of_left_child :: } //A similar logic for mass_u values - sqd_frobenius_norm_of_left_child=SquaredFrobeniusNorm_(left_stat->b_twy_mass_u); - sqd_frobenius_norm_of_right_child=SquaredFrobeniusNorm_(right_stat->b_twy_mass_u); + sqd_frobenius_norm_of_left_child= + SquaredFrobeniusNorm_(left_stat->b_twy_mass_u); - //Now compute it with that of the parent + sqd_frobenius_norm_of_right_child= + SquaredFrobeniusNorm_(right_stat->b_twy_mass_u); + + //Now compare it with that of the parent sqd_frobenius_norm_of_parent=SquaredFrobeniusNorm_(parent_stat.b_twy_mass_u); @@ -680,6 +708,43 @@ MergeChildBoundsB_TWY_( FastRegression:: * and the children */ +template +void FastRegression:: +MergeChildBounds_( Tree *qnode,prune_t flag){ + + //We will merge bounds depending on the value of the flag + + //But firstly we will check if the qnode is a leaf node. if it is + //then we will simply return + + if(qnode->is_leaf()){ + return; + } + + //The statistics of the parent and the children node + FastRegressionStat &parent_stat=qnode->stat(); + FastRegressionStat *left_stat=&(qnode->left()->stat()); + FastRegressionStat *right_stat=&(qnode->right()->stat()); + if(flag==PRUNE_B_TWB){ + + MergeChildBoundsB_TWB_(left_stat,right_stat,parent_stat); + + } + else{ + if(flag==PRUNE_B_TWY){ + + MergeChildBoundsB_TWY_(left_stat,right_stat,parent_stat); + } + else{ + if(flag==PRUNE_BOTH){ + MergeChildBoundsB_TWB_(left_stat,right_stat,parent_stat); + MergeChildBoundsB_TWY_(left_stat,right_stat,parent_stat); + } + } + } +} + + template void FastRegression:: MergeChildBounds_( Tree *qnode,check_for_prune_t flag){ @@ -708,10 +773,10 @@ MergeChildBounds_( Tree *qnode,check_for_prune_t flag){ MergeChildBoundsB_TWY_(left_stat,right_stat,parent_stat); } else{ - //printf("will merge both the bounds...\n"); - MergeChildBoundsB_TWB_(left_stat,right_stat,parent_stat); - MergeChildBoundsB_TWY_(left_stat,right_stat,parent_stat); - + if(flag==CHECK_FOR_PRUNE_BOTH){ + MergeChildBoundsB_TWB_(left_stat,right_stat,parent_stat); + MergeChildBoundsB_TWY_(left_stat,right_stat,parent_stat); + } } } } @@ -725,13 +790,15 @@ template void FastRegression::FRegressionBaseB_TWY_(Tree *qnode, Tree *rnode){ + //printf("Hit the base case for BTWY..\n"); + //subtract along each dimension as we are now doing exhaustive calculations //more_u <- more_u - rnode->stat().b_ty la::SubFrom (rnode->stat().b_ty, &qnode->stat().b_twy_more_u); - //Having subtracted calculate B^TWY exhaustively + //Having subtracted, calculate B^TWY exhaustively //One can do this by using linear algebra routines available //in LaPack. however we shall not use them because B can be a //very large matrix @@ -759,27 +826,27 @@ void FastRegression::FRegressionBaseB_TWY_(Tree *qnode, Tree *rnode){ //This is nothing but B^TWY being evaluated dimension wise if (row != 0){ - double val=b_twy_l_estimate_[q].get(row,col)+ + double val1=b_twy_l_estimate_[q].get(row,col)+ ker_value * rset_weights_[old_from_new_r_[r]] * rset_.get (row- 1, r); - b_twy_l_estimate_[q].set(row,col,val); + b_twy_l_estimate_[q].set(row,col,val1); - val=b_twy_u_estimate_[q].get(row,col)+ + double val2 =b_twy_u_estimate_[q].get(row,col)+ ker_value * rset_weights_[old_from_new_r_[r]] * rset_.get (row- 1, r); - b_twy_u_estimate_[q].set(row,col,val); + b_twy_u_estimate_[q].set(row,col,val2); } else{ - double val= b_twy_l_estimate_[q].get(row,col) + + double val1= b_twy_l_estimate_[q].get(row,col) + ker_value * rset_weights_[old_from_new_r_[r]]; - b_twy_l_estimate_[q].set(row,col,val); + b_twy_l_estimate_[q].set(row,col,val1); - val= b_twy_u_estimate_[q].get(row,col) + + double val2= b_twy_u_estimate_[q].get(row,col) + ker_value * rset_weights_[old_from_new_r_[r]]; - b_twy_u_estimate_[q].set(row,col,val); + b_twy_u_estimate_[q].set(row,col,val2); } @@ -833,13 +900,15 @@ void FastRegression::FRegressionBaseB_TWY_(Tree *qnode, Tree *rnode){ //having looped over each point qnode->stat().b_twy_mass_u.CopyValues(max_u); - qnode->stat().b_twy_mass_l.CopyValues(min_l); - + qnode->stat().b_twy_mass_l.CopyValues(min_l); + return; + } + else{ //We are doing frobenius norm pruning - /** get a tighter lower and upper boiounf by looping over each query + /** get a tighter lower and upper bound by looping over each query * point to find that particular query point that has the least sqd * frobenius norm of the matrix b_twy_l_estimate_[q]+qnode->stat().b_twy_more_l */ @@ -850,8 +919,8 @@ void FastRegression::FRegressionBaseB_TWY_(Tree *qnode, Tree *rnode){ //The min_pointer and the max_pointer store the index number of the //query point which possibly has the least and the highest squared frobenius norm - index_t min_pointer; - index_t max_pointer; + index_t min_pointer=qnode->begin(); + index_t max_pointer=qnode->begin(); temp.Init(rset_.n_rows()+1,1); for(index_t i=qnode->begin();iend();i++){ @@ -879,9 +948,13 @@ void FastRegression::FRegressionBaseB_TWY_(Tree *qnode, Tree *rnode){ //Once done with the looping process set up the mass_l and mass_u values - la::AddOverwrite(b_twy_l_estimate_[min_pointer],qnode->stat().b_twy_more_l,&qnode->stat().b_twy_mass_l); - la::AddOverwrite(b_twy_u_estimate_[max_pointer],qnode->stat().b_twy_more_u,&qnode->stat().b_twy_mass_u); + la::AddOverwrite(b_twy_l_estimate_[min_pointer], + qnode->stat().b_twy_more_l,&qnode->stat().b_twy_mass_l); + + la::AddOverwrite(b_twy_u_estimate_[max_pointer], + qnode->stat().b_twy_more_u,&qnode->stat().b_twy_mass_u); + } } @@ -890,6 +963,7 @@ void FastRegression::FRegressionBaseB_TWY_(Tree *qnode, Tree *rnode){ template void FastRegression::FRegressionBaseB_TWB_(Tree *qnode,Tree *rnode){ + printf("Hit the base case for BTWB\n"); la::SubFrom (rnode->stat().b_tb, &qnode->stat().b_twb_more_u); @@ -973,7 +1047,7 @@ void FastRegression::FRegressionBaseB_TWB_(Tree *qnode,Tree *rnode){ } } } - } + } if(pruning_criteria==CRITERIA_FOR_PRUNE_COMPONENT){ //Loop over each point and set the max and min @@ -1022,35 +1096,37 @@ void FastRegression::FRegressionBaseB_TWB_(Tree *qnode,Tree *rnode){ qnode->stat().b_twb_mass_u.CopyValues(max_u); qnode->stat().b_twb_mass_l.CopyValues(min_l); + return; } //We are interested in frobenius norm pruning + else{ - - //We are doing frobenius norm pruning /** get a tighter lower and upper boiounf by looping over each query * point to find that particular query point that has the least sqd - * frobenius norm of the matrix b_twy_l_estimate_[q]+qnode->stat().b_twy_more_l + * frobenius norm of the matrix b_twb_l_estimate_[q]+qnode->stat().b_twb_more_l */ double min_norm=DBL_MAX; double max_norm=DBL_MIN; - Matrix temp; + Matrix temp1; + Matrix temp2; //The min_pointer and the max_pointer store the index number of the //query point which possibly has the least and the highest squared frobenius norm - index_t min_pointer; - index_t max_pointer; + index_t min_pointer=qnode->begin(); + index_t max_pointer=qnode->begin(); - temp.Init(rset_.n_rows()+1,rset_.n_rows()+1); + temp1.Init(rset_.n_rows()+1,rset_.n_rows()+1); + temp2.Init(rset_.n_rows()+1,rset_.n_rows()+1); for(index_t i=qnode->begin();iend();i++){ //look for the lower bound - la::AddOverwrite(b_twb_l_estimate_[i],qnode->stat().b_twb_more_l,&temp); - - double var=SquaredFrobeniusNorm_(temp); + + la::AddOverwrite(b_twb_l_estimate_[i],qnode->stat().b_twb_more_l,&temp1); + double var=SquaredFrobeniusNorm_(temp1); if(var< min_norm){ min_pointer=i; @@ -1058,20 +1134,22 @@ void FastRegression::FRegressionBaseB_TWB_(Tree *qnode,Tree *rnode){ } //Look for the upper bound - la::AddOverwrite(b_twb_u_estimate_[i],qnode->stat().b_twb_more_u,&temp); - var=SquaredFrobeniusNorm_(temp); - if(var>max_norm){ + la::AddOverwrite(b_twb_u_estimate_[i],qnode->stat().b_twb_more_u,&temp2); + var=SquaredFrobeniusNorm_(temp2); + if(var > max_norm){ max_pointer=i; max_norm=var; } - } //Once done with the looping process set up the mass_l and mass_u values - la::AddOverwrite(b_twb_l_estimate_[min_pointer],qnode->stat().b_twb_more_l,&qnode->stat().b_twb_mass_l); - la::AddOverwrite(b_twb_u_estimate_[max_pointer],qnode->stat().b_twb_more_u,&qnode->stat().b_twb_mass_u); + la::AddOverwrite(b_twb_l_estimate_[min_pointer], + qnode->stat().b_twb_more_l,&qnode->stat().b_twb_mass_l); + + la::AddOverwrite(b_twb_u_estimate_[max_pointer], + qnode->stat().b_twb_more_u,&qnode->stat().b_twb_mass_u); } } @@ -1118,31 +1196,28 @@ void FastRegression::PostProcess_(Tree *qnode){ //mass_u <- owed_u+mass_u - - + + UpdateBounds_(qnode, qnode->stat().b_twb_owed_l, + qnode->stat().b_twb_owed_u, + qnode->stat().b_twy_owed_l, + qnode->stat().b_twy_owed_u, CHECK_FOR_PRUNE_BOTH); + if(!qnode->is_leaf()){ //Having completed all the tree calculations we now update bounds for both //the qunatities. this can be done by calling the UpdateBounds_ function //with the flag set to CHECK_FOR_PRUNE_BOTH - - UpdateBounds_(qnode, qnode->stat().b_twb_owed_l, - qnode->stat().b_twb_owed_u, - qnode->stat().b_twy_owed_l, - qnode->stat().b_twy_owed_u, CHECK_FOR_PRUNE_BOTH); + PostProcess_(qnode->left()); PostProcess_(qnode->right()); } else{ - UpdateBounds_(qnode, qnode->stat().b_twb_owed_l, - qnode->stat().b_twb_owed_u, - qnode->stat().b_twy_owed_l, - qnode->stat().b_twy_owed_u, CHECK_FOR_PRUNE_BOTH); + // b_twb_e_estimate= - //0.5*(b_twy_l_estimate+b_twy_u_estimate+b_twy_more_l+b_twy_more_u) + //0.5*(b_twb_l_estimate+b_twb_u_estimate+b_twb_more_l+b_twb_more_u) for(index_t q=qnode->begin();qend();q++){ FastRegressionStat qstat=qnode->stat(); @@ -1150,16 +1225,16 @@ void FastRegression::PostProcess_(Tree *qnode){ la::AddInit(b_twb_l_estimate_[q], b_twb_u_estimate_[q],&estimate_mean1); - la::Scale(0.50,&estimate_mean1); + // la::Scale(0.50,&estimate_mean1); Matrix more_mean1; la::AddInit(qstat.b_twb_more_l,qstat.b_twb_more_u,&more_mean1); - la::Scale(0.50,&more_mean1); + //la::Scale(0.50,&more_mean1); //Final estimate is just the sum of estimate_mean and more_mean la::AddOverwrite(estimate_mean1,more_mean1,&b_twb_e_estimate_[q]); - + la::Scale(0.50,&b_twb_e_estimate_[q]); // b_twy_e_estimate= //0.5*(b_twy_l_estimate+b_twy_u_estimate+b_twy_more_l+b_twy_more_u) @@ -1168,15 +1243,17 @@ void FastRegression::PostProcess_(Tree *qnode){ la::AddInit(b_twy_l_estimate_[q], b_twy_u_estimate_[q],&estimate_mean2); - la::Scale(0.50,&estimate_mean2); + // la::Scale(0.50,&estimate_mean2); Matrix more_mean2; la::AddInit(qstat.b_twy_more_l,qstat.b_twy_more_u,&more_mean2); - la::Scale(0.50,&more_mean2); + // la::Scale(0.50,&more_mean2); //Final estimate is just the sum of estimate_mean and more_mean la::AddOverwrite(estimate_mean2,more_mean2,&b_twy_e_estimate_[q]); + la::Scale(0.50,&b_twy_e_estimate_[q]); + } } } @@ -1298,6 +1375,20 @@ void FastRegression::PrintRegressionEstimate_(){ fclose(gp); } +template +void FastRegression::CompareWithTrueValues_(){ + //We have regression_estimate_ and true_regression_values_.Lets get + //mean squared error + + double error=0.0; + for(index_t q=0;q::PrintRegressionEstimate_(){ - - - diff --git a/fastlib/u/gmravi/regression/regression_ll_naive.h b/fastlib/u/gmravi/regression/regression_ll_naive.h index 708604ea1c..457c5ab819 100644 --- a/fastlib/u/gmravi/regression/regression_ll_naive.h +++ b/fastlib/u/gmravi/regression/regression_ll_naive.h @@ -1,7 +1,7 @@ #ifndef REGRESSION_LL_NAIVE_H #define REGRESSION_LL_NAIVE_H #include "fastlib/fastlib_int.h" -#include "pseudo_inverse.h" +//#include "pseudo_inverse.h" template double NaiveCalculation::SquaredFrobeniusNorm_(Matrix &a){ @@ -9,9 +9,9 @@ double NaiveCalculation::SquaredFrobeniusNorm_(Matrix &a){ //This function computes treats the natrix as a vector and computes //it's 1-norm double value=0; - for(index_t col=0;col::Compute (){ //Fill B^TWY naive double val=b_twy_naive_estimate_[q].get(row,col)+ - ker_value * rset_weights_[old_from_new_r_[r]] * rset_.get (row- 1, r); + ker_value * rset_weights_[old_from_new_r_[r]] * rset_.get (row-1, r); b_twy_naive_estimate_[q].set(row,col,val); @@ -75,6 +75,7 @@ void NaiveCalculation::Compute (){ } else{ + //this is row =0 //Fill B^TWY naive @@ -120,7 +121,7 @@ void NaiveCalculation::Compute (){ Print_(); //Having done this get regression estimates by calling the function //ObtainRegressionEstimates - ObtainRegressionEstimate_(); + //ObtainRegressionEstimate_(); //PrintRegressionEstimate_(fast_regression_estimate); } @@ -181,10 +182,16 @@ void NaiveCalculation::ObtainRegressionEstimate_(){ template -void NaiveCalculation:: ComputeMaximumRelativeError(ArrayList &fast_b_twy_estimate, ArrayList & fast_b_twb_estimate, char *pruning_criteria){ +void NaiveCalculation:: +ComputeMaximumRelativeError(ArrayList &fast_b_twy_estimate, + ArrayList & fast_b_twb_estimate, + char *pruning_criteria){ if(!strcmp(pruning_criteria,"fnorm")){ + FILE *fp; + fp=fopen("errors.txt","w+"); + double max_frobenius_error_b_twb=0; double max_frobenius_error_b_twy=0; @@ -195,27 +202,49 @@ void NaiveCalculation:: ComputeMaximumRelativeError(ArrayList & la::SubInit (b_twy_naive_estimate_[q],fast_b_twy_estimate[q] , &temp1); double f_norm=SquaredFrobeniusNorm_(temp1); double f_norm_naive=SquaredFrobeniusNorm_(b_twy_naive_estimate_[q]); - double rel_error_b_twy=fabs(f_norm)/f_norm_naive; - //printf("relative frobenius norm error for BTWY is %f\n",rel_error_b_twy); + double rel_error_b_twy; + if(f_norm==0.0 && f_norm_naive==0.0){ + rel_error_b_twy=0; + } + else{ + rel_error_b_twy=fabs(f_norm)/f_norm_naive; + } + fprintf(fp,"\n"); + fprintf(fp,"diff norm is %f\n",f_norm); + fprintf(fp,"naive norm is %f\n",f_norm_naive); + fprintf(fp,"relative frobenius norm error for BTWY is %f\n",rel_error_b_twy); if(max_frobenius_error_b_twy:: ComputeMaximumRelativeError(ArrayList & fabs(fast_b_twb_estimate[q].get(row,col)-b_twb_naive_estimate_[q].get(row,col)); double error=diff/b_twb_naive_estimate_[q].get(row,col); + + printf("diff =%f\n",diff); + printf("Being compared against %f\n",b_twb_naive_estimate_[q].get(row,col)); if(error>max_error_for_this_point){ max_error_for_this_point=error; @@ -245,8 +277,8 @@ void NaiveCalculation:: ComputeMaximumRelativeError(ArrayList & } }//Calculations complete for this point - printf("Maximum error for this point is %f\n",max_error_for_this_point); - if(max_error_for_this_pointmax_error_on_the_whole){ max_error_on_the_whole=max_error_for_this_point; } @@ -266,7 +298,10 @@ void NaiveCalculation:: ComputeMaximumRelativeError(ArrayList & double diff= fabs(fast_b_twy_estimate[q].get(row,col)-b_twy_naive_estimate_[q].get(row,col)); - double error=diff/b_twb_naive_estimate_[q].get(row,col); + double error=diff/b_twy_naive_estimate_[q].get(row,col); + + printf("diff =%f\n",diff); + printf("Being compared against %f\n",b_twy_naive_estimate_[q].get(row,col)); if(error>max_error_for_this_point){ max_error_for_this_point=error; @@ -274,9 +309,9 @@ void NaiveCalculation:: ComputeMaximumRelativeError(ArrayList & } } //Calculations complete for this point - printf("Maximum error for this point is %f\n",max_error_for_this_point); + printf("Maximum error for this point BTWY is %f\n",max_error_for_this_point); - if(max_error_for_this_pointmax_error_on_the_whole){ max_error_on_the_whole=max_error_for_this_point; } } diff --git a/fastlib/u/gmravi/regression/regression_parent.h b/fastlib/u/gmravi/regression/regression_parent.h index 54130bba7b..6df46da5bb 100644 --- a/fastlib/u/gmravi/regression/regression_parent.h +++ b/fastlib/u/gmravi/regression/regression_parent.h @@ -207,6 +207,9 @@ public: /** Regression estimates of the query points */ Vector regression_estimate_; + /* True regression values of the query points */ + Vector true_regression_values_; + /** Mappings from old dataset to new dataset * Remember that when we build the query tree and the * reference tree out of the query and reference @@ -246,8 +249,9 @@ public: void FRegressionBase_(Tree *qnode, Tree *rnode, check_for_prune_t flag); - void MergeChildBounds_(Tree *,check_for_prune_t flag); + void MergeChildBounds_(Tree *,prune_t flag); + void MergeChildBounds_(Tree *,check_for_prune_t flag); void CallRecursively_(Tree *qnode, Tree *rnode,check_for_prune_t flag); @@ -283,6 +287,8 @@ public: void PrintRegressionEstimate_(); double SquaredFrobeniusNorm_(Matrix &); double Compute1NormLike_(Matrix &); + + void CompareWithTrueValues_(); public: @@ -305,7 +311,7 @@ public: void Compute(); void Init(Matrix &q_matrix, Matrix &r_matrix, double bandwidth, - double tau,index_t leaf_length,Vector &rset_weights, char *); + double tau,index_t leaf_length,Vector &rset_weights, char *, Vector &); };