diff --git a/fastlib2/contrib/niche/kernel_vector_mult.h b/fastlib2/contrib/niche/kernel_vector_mult.h index c180763aae..20b5ce750b 100644 --- a/fastlib2/contrib/niche/kernel_vector_mult.h +++ b/fastlib2/contrib/niche/kernel_vector_mult.h @@ -131,9 +131,9 @@ class KernelVectorMult { // kernel function - EpanKernel epan_kernel; + GaussianKernel kernel; - double epan_cutoff_dist_; + double cutoff_dist_; @@ -276,7 +276,7 @@ class KernelVectorMult { la::DistanceSqEuclidean(query_point, reference_point); weighted_sums_[query_index] += - weights_[reference_index] * epan_kernel.EvalUnnormOnSq(distance); + weights_[reference_index] * kernel.EvalUnnormOnSq(distance); } /* for reference_index */ @@ -305,7 +305,7 @@ class KernelVectorMult { DEBUG_SAME_DOUBLE(lower_bound_distance, MinNodeDistSq_(query_node, reference_node)); - if (lower_bound_distance > epan_cutoff_dist_) { + if (0) {//lower_bound_distance > cutoff_dist_) { // execute prune @@ -465,11 +465,11 @@ class KernelVectorMult { DEBUG_ASSERT(bandwidth > 0); // kernel function - epan_kernel.Init(bandwidth); + kernel.Init(bandwidth); - epan_cutoff_dist_ = bandwidth*bandwidth; + cutoff_dist_ = bandwidth*bandwidth; - DEBUG_ONLY(printf("epan_cutoff_dist = %f\n", epan_cutoff_dist_)); + DEBUG_ONLY(printf("cutoff_dist = %f\n", cutoff_dist_)); } /* Init */ @@ -502,7 +502,7 @@ class KernelVectorMult { 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_); + la::Scale(1 / kernel.CalcNormConstant(queries_.n_rows()), &weighted_sums_); fx_timer_stop(module_, "dual_tree_computation"); diff --git a/fastlib2/contrib/niche/solvelinsys.h b/fastlib2/contrib/niche/solvelinsys.h index e58ce7c211..dcdcf01fac 100644 --- a/fastlib2/contrib/niche/solvelinsys.h +++ b/fastlib2/contrib/niche/solvelinsys.h @@ -30,6 +30,7 @@ private: // used for a simple linear transformation operator data_^T * data Matrix data_; + float** K_; //Matrix K_; //EpanKernel epan_kernel_; @@ -55,6 +56,7 @@ public: + @@ -73,9 +75,30 @@ public: sigma_squared_ = sigma_squared_in; - /* + // for debugging purposes, explicitly represent kernel matrix K - K_.Init(n_points_, n_points_); + + if((K_ = (float**) malloc(40000 * sizeof(float*))) == NULL) { + printf("failed\n"); + exit(1); + } + + for(int i = 0; i < 40000; i++) { + if((K_[i] = (float*) malloc(40000 * sizeof(float))) == NULL) { + printf("failed on %d\n", i); + exit(1); + } + } + + printf("start writing\n"); + for(int i = 0; i < 40000; i++) { + for(int j = 0; j < 40000; j++) { + K_[i][j] = 0; + } + } + printf("done writing\n"); + + for(int i = 0; i < n_points_; i++) { Vector v_i; data_.MakeColumnVector(i, &v_i); @@ -84,22 +107,18 @@ 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_.set(i, j, epan_kernel_.EvalUnnormOnSq(dist)); + K_[i][j] = kernel_.EvalUnnormOnSq(dist) / norm_constant_; } } - la::Scale(1 / norm_constant_, &K_); for(int i = 0; i < n_points_; i++) { - K_.set(i, i, K_.get(i, i) + sigma_squared_); + K_[i][i] += sigma_squared_; } - const char *K_file_name = "K.txt"; - data::Save(K_file_name, K_); + //const char *K_file_name = "K.txt"; + //data::Save(K_file_name, K_); // int errors = 0; @@ -116,7 +135,7 @@ public: // end debugging explicit representation of K - */ + } @@ -344,8 +363,8 @@ void SolveLinearSystem(Matrix references, Vector rhs, double bandwidth, double s //iterative_solver.SetAztecOption(AZ_precond, AZ_Jacobi); // Use Conjugate Gradient - //iterative_solver.SetAztecOption(AZ_solver, AZ_cg); - iterative_solver.SetAztecOption(AZ_solver, AZ_gmres); + 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.