diff --git a/fastlib2/contrib/niche/build.py b/fastlib2/contrib/niche/build.py index d1958ed1ec..5f27a66bcc 100644 --- a/fastlib2/contrib/niche/build.py +++ b/fastlib2/contrib/niche/build.py @@ -1,6 +1,6 @@ binrule( - name = "test2", - sources = ["test2.cc"], + name = "solvelinsys", + sources = ["solvelinsys.cc"], headers = ["kernel_vector_mult.h"], deplibs = ["contrib/dongryel/regression:krylov_lpr", "fastlib:fastlib_int", diff --git a/fastlib2/contrib/niche/kernel_vector_mult.h b/fastlib2/contrib/niche/kernel_vector_mult.h index e837295d03..7e8d985016 100644 --- a/fastlib2/contrib/niche/kernel_vector_mult.h +++ b/fastlib2/contrib/niche/kernel_vector_mult.h @@ -407,7 +407,7 @@ class KernelVectorMult { /** * Read parameters, copy data into the class, and build the trees. */ - void Init(const Matrix& references_in, struct datanode* module_in) { + void Init(const Matrix& references_in, double bandwidth, struct datanode* module_in) { // It's a good idea to make sure the object isn't initialized a // second time, as this is almost certainly mistaken. @@ -462,15 +462,14 @@ class KernelVectorMult { number_of_prunes_ = 0; - double bandwidth = fx_param_int(module_, "bandwidth", 1); DEBUG_ASSERT(bandwidth > 0); // kernel function epan_kernel.Init(bandwidth); - epan_cutoff_dist_ = epan_kernel.MaxUnnormValue(); + epan_cutoff_dist_ = bandwidth*bandwidth; - printf("epan_cutoff_dist = %f\n", epan_cutoff_dist_); + DEBUG_ONLY(printf("epan_cutoff_dist = %f\n", epan_cutoff_dist_)); } /* Init */ @@ -502,7 +501,7 @@ class KernelVectorMult { /* Start recursion on the roots of either tree */ GNPRecursion_(query_tree_, reference_tree_, MinNodeDistSq_(query_tree_, reference_tree_)); - + printf("queries_.n_rows() = %d\n", queries_.n_rows()); la::Scale(1 / epan_kernel.CalcNormConstant(queries_.n_rows()), &weighted_sums_); fx_timer_stop(module_, "dual_tree_computation"); diff --git a/fastlib2/contrib/niche/test2.cc b/fastlib2/contrib/niche/test2.cc index bb8b8086a6..c6d550d4f0 100644 --- a/fastlib2/contrib/niche/test2.cc +++ b/fastlib2/contrib/niche/test2.cc @@ -26,12 +26,12 @@ private: // used for a simple linear transformation operator data_^T * data Matrix data_; - // Matrix K_; + //Matrix K_; - EpanKernel epan_kernel_; - double norm_constant_; + //EpanKernel epan_kernel_; + //double norm_constant_; - double sigma_; + double sigma_squared_; KernelVectorMult* kernel_vector_mult_; @@ -40,8 +40,7 @@ private: public: - SimpleLinearOperator(int n_points_in, Matrix data_in, - double bandwidth_in, double sigma_in, + SimpleLinearOperator(int n_points_in, Matrix data_in, double sigma_squared_in, KernelVectorMult* kernel_vector_mult_in) { n_points_ = n_points_in; @@ -50,8 +49,12 @@ public: data_.Copy(data_in); - //la::MulTransBInit(data_in, data_in, &K_); - //data::Load("K.txt", &K_); + + + + + + n_dims_ = data_in.n_rows(); //n_points_ = data_in.n_cols(); @@ -60,10 +63,56 @@ public: map = new Epetra_Map(n_points_, 0, comm); - epan_kernel_.Init(bandwidth_in, n_dims_); - norm_constant_ = epan_kernel_.CalcNormConstant(n_dims_); + //epan_kernel_.Init(bandwidth_in, n_dims_); + //norm_constant_ = epan_kernel_.CalcNormConstant(n_dims_); + //printf("norm_constant = %f\n", norm_constant_); - sigma_ = sigma_in; + sigma_squared_ = sigma_squared_in; + + /* + // for debugging purposes, explicitly represent kernel matrix K + K_.Init(n_points_, n_points_); + for(int i = 0; i < n_points_; i++) { + Vector v_i; + data_.MakeColumnVector(i, &v_i); + + for(int j = 0; j < n_points_; j++) { + Vector v_j; + data_.MakeColumnVector(j, &v_j); + + //double dist = la::Dot(v_i, v_j); + //K_i[j] = dist; + + double dist = la::DistanceSqEuclidean(v_i, v_j); + K_.set(i, j, epan_kernel_.EvalUnnormOnSq(dist)); + } + } + + la::Scale(1 / norm_constant_, &K_); + for(int i = 0; i < n_points_; i++) { + K_.set(i, i, K_.get(i, i) + sigma_squared_); + } + + + const char *K_file_name = "K.txt"; + data::Save(K_file_name, K_); + + +// int errors = 0; +// for(int i = 0; i < n_points_; i++) { +// for(int j = 0; j < n_points_; j++) { +// if(K_.get(i,j) != K_.get(j,i)) { +// errors++; +// } +// } +// } + +// printf("sigma_squared_ = %f\n", sigma_squared_); +// printf("errors = %d\n", errors); + + + // end debugging explicit representation of K + */ } @@ -79,7 +128,8 @@ public: int Apply (const Epetra_MultiVector &X, Epetra_MultiVector &Y) const { - /* fast summation code */ + + // FAST SUMMATION CODE Vector weights_vector; Vector results; weights_vector.Init(n_points_); @@ -93,19 +143,24 @@ public: kernel_vector_mult_ -> Reset(); kernel_vector_mult_ -> ComputeKernelMatrixVectorMultiplication(weights_vector, &results); - + printf("\n"); for(int i = 0; i < n_points_; i++) { - Y.Pointers()[0][i] = results[i] + (sigma_ * weights_vector[i]); + Y.Pointers()[0][i] = results[i] + (sigma_squared_ * weights_vector[i]); + //results[i] += sigma_squared_ * weights_vector[i]; + //printf("%f ", results[i]); } - - /* end fast summation code */ - + printf("\n"); + + // END FAST SUMMATION CODE + + printf("data.n_rows() = %d\ndata.n_cols() = %d\n", + data_.n_rows(), + data_.n_cols()); + /* - // LET'S ASSUME THAT DATA IS A D X N MATRIX - Vector K_i; K_i.Init(n_points_); @@ -118,28 +173,54 @@ public: Vector v_j; data_.MakeColumnVector(j, &v_j); + //double dist = la::Dot(v_i, v_j); + //K_i[j] = dist; + double dist = la::DistanceSqEuclidean(v_i, v_j); - - K_i[j] = epan_kernel_.EvalUnnormOnSq(dist) / norm_constant_; + K_i[j] = epan_kernel_.EvalUnnormOnSq(dist); } + + la::Scale(1 / norm_constant_, &K_i); + - //K_i[i] += sigma_; // add (sigma * I) term + K_i[i] += sigma_squared_; // add (sigma^2 * I) term double sum = 0; for(int j = 0; j < n_points_; j++) { - sum += K_i[j] * X[0][j]; + sum += (K_i[j] * X.Pointers()[0][j]); } Y.Pointers()[0][i] = sum; } + */ + /* + // explicit K + Vector x_vec; + x_vec.Init(n_points_); + for(int i = 0; i < n_points_; i++) { + x_vec[i] = X.Pointers()[0][i]; + } + + Vector y_vec; + la::MulInit(K_, x_vec, &y_vec); + + for(int i = 0; i < n_points_; i++) { + Y.Pointers()[0][i] = y_vec[i]; + } + + // end explicit K + + //y_vec.PrintDebug("y_vec"); + */ + + /* double my_squared_error = 0; for(int i = 0; i < n_points_; i++) { my_squared_error += pow(Y.Pointers()[0][i] - results[i], 2); } - printf("my_squared_error = %f\n", my_squared_error); - */ + */ @@ -190,81 +271,54 @@ public: +// A(r) x = rhs +// A(r) is some linear operator derived from reference points r +void SolveLinearSystem(Matrix references, Vector rhs, double bandwidth, double sigma_squared, Vector* solution) { - -int main(int argc, char *argv[]) { - - // Initialize FastExec... - fx_init(argc, argv); - - - - /* begin code from KernelVectorMult */ + DEBUG_ONLY(printf("references.n_rows() = %d\nreferences.n_cols() = %d\n", + references.n_rows(), references.n_cols())); - // The reference data file is a required parameter. - const char* references_file_name = fx_param_str_req(NULL, "r"); - - Matrix references; - data::Load(references_file_name, &references); - - printf("references.n_rows() = %d\nreferences.n_cols() = %d\n", - references.n_rows(), references.n_cols()); - KernelVectorMult kernel_vector_mult; struct datanode* kernel_vector_mult_module = fx_submodule(NULL, "kernel_vector_mult", "kernel_vector_mult_module"); + + + DEBUG_ONLY(printf("sigma_squared = %f\n", sigma_squared)); + + kernel_vector_mult.Init(references, bandwidth, kernel_vector_mult_module); - kernel_vector_mult.Init(references, kernel_vector_mult_module); - - /* end code from KernelVectorMult */ - - - - - Matrix data;//, data_transpose; - Matrix right_hand_side_e; - - data::Load("refined_astroset.ds", &data); - data::Load("alldata_zs", &right_hand_side_e); - /* - if(data_transpose.n_cols() < data_transpose.n_rows()) { - la::TransposeInit(data_transpose, &data); - } - else { - data = data_transpose; - } - */ + // Communication stuff? Epetra_SerialComm comm; - int row_length = data.n_cols(); // this should be correct despite obvious semantic conflicts + int row_length = references.n_cols(); // this should be correct despite obvious semantic conflicts // Required map for multivector stuff Epetra_BlockMap blockmap(1, row_length, 0, comm); - printf("row_length = %d\n", row_length); + DEBUG_ONLY(printf("row_length = %d\n", row_length)); // Define the linear problem. - Epetra_MultiVector solution(blockmap, 1, true); - Epetra_MultiVector right_hand_side(blockmap, 1, false); + Epetra_MultiVector solution_e(blockmap, 1, true); + Epetra_MultiVector right_hand_side_e(blockmap, 1, false); for(index_t j = 0; j < row_length; j++) { - (*(right_hand_side(0)))[j] = right_hand_side_e.get(0, j); + (*(right_hand_side_e(0)))[j] = rhs[j]; } - SimpleLinearOperator simple_linear_operator(row_length, data, 1, 1, - &kernel_vector_mult); + SimpleLinearOperator simple_linear_operator(row_length, references, + sigma_squared, &kernel_vector_mult); Epetra_LinearProblem linear_problem(&simple_linear_operator, - &solution, &right_hand_side); + &solution_e, &right_hand_side_e); // Declare the iterative solver. AztecOO iterative_solver; @@ -288,22 +342,58 @@ int main(int argc, char *argv[]) { // Use Conjugate Gradient //iterative_solver.SetAztecOption(AZ_solver, AZ_cg); iterative_solver.SetAztecOption(AZ_solver, AZ_gmres); + //iterative_solver.SetAztecOption(AZ_solver, AZ_cg_condnum); // Use modified Gram-Schmidt. - //iterative_solver.SetAztecOption(AZ_orthog, AZ_modified); + iterative_solver.SetAztecOption(AZ_orthog, AZ_modified); // No output. //iterative_solver.SetAztecOption(AZ_diagnostics, AZ_none); //iterative_solver.SetAztecOption(AZ_output, AZ_none); // Solve the linear system. - iterative_solver.Iterate(row_length, 1.0E-9); + iterative_solver.Iterate(row_length, 1.0E-6); + + solution -> Init(row_length); + for(index_t j = 0; j < row_length; j++) { + (*solution)[j] = solution_e.Pointers()[0][j]; + } - //cout<