From f7fb517061b849c2fbf45c02ea999dc98ec0b4af Mon Sep 17 00:00:00 2001 From: Bill March Date: Sun, 20 Apr 2008 18:21:42 +0000 Subject: [PATCH] error testing stuff --- fastlib2/contrib/march/error/hybrid_error.cc | 70 ++++ fastlib2/contrib/march/error/hybrid_error.h | 258 +++++++++++++++ .../march/error/hybrid_error_analysis.h | 122 +++++++ .../contrib/march/error/hybrid_error_stat.h | 302 ++++++++++++++++++ .../contrib/march/error/naive_kernel_sum.h | 92 ++++++ 5 files changed, 844 insertions(+) create mode 100644 fastlib2/contrib/march/error/hybrid_error.cc create mode 100644 fastlib2/contrib/march/error/hybrid_error.h create mode 100644 fastlib2/contrib/march/error/hybrid_error_analysis.h create mode 100644 fastlib2/contrib/march/error/hybrid_error_stat.h create mode 100644 fastlib2/contrib/march/error/naive_kernel_sum.h diff --git a/fastlib2/contrib/march/error/hybrid_error.cc b/fastlib2/contrib/march/error/hybrid_error.cc new file mode 100644 index 0000000000..dddccccfc8 --- /dev/null +++ b/fastlib2/contrib/march/error/hybrid_error.cc @@ -0,0 +1,70 @@ +#include "hybrid_error.h" +#include "hybrid_error_stat.h" +#include "hybrid_error_analysis.h" +#include "naive_kernel_sum.h" +#include "fastlib/fastlib.h" + +int main(int argc, char* argv[]) { + + fx_init(argc, argv); + + GaussianKernelErrorTester absolute; + GaussianKernelErrorTester relative; + GaussianKernelErrorTester exponential; + GaussianKernelErrorTester gaussian; + GaussianKernelErrorTester hybrid; + + Matrix centers; + + const char* dataset = fx_param_str(NULL, "data", "test_data.csv"); + data::Load(dataset, ¢ers); + + double bandwidth = fx_param_double(NULL, "bandwidth", 0.1); + DEBUG_ASSERT(bandwidth > 0.0); + + Vector abs_results; + struct datanode* abs_mod = fx_submodule(NULL, "abs", "absolute"); + absolute.Init(abs_mod, centers, bandwidth); + absolute.ComputeTotalSum(&abs_results); + + Vector rel_results; + struct datanode* rel_mod = fx_submodule(NULL, "rel", "relative"); + relative.Init(rel_mod, centers, bandwidth); + relative.ComputeTotalSum(&rel_results); + + Vector exp_results; + struct datanode* exp_mod = fx_submodule(NULL, "exp", "exponential"); + exponential.Init(exp_mod, centers, bandwidth); + exponential.ComputeTotalSum(&exp_results); + + Vector gauss_results; + struct datanode* gauss_mod = fx_submodule(NULL, "gauss", "gaussian"); + gaussian.Init(gauss_mod, centers, bandwidth); + gaussian.ComputeTotalSum(&gauss_results); + + Vector hybrid_results; + struct datanode* hybrid_mod = fx_submodule(NULL, "hybrid", "hybrid"); + hybrid.Init(hybrid_mod, centers, bandwidth); + hybrid.ComputeTotalSum(&hybrid_results); + + NaiveKernelSum naive; + Vector naive_results; + struct datanode* naive_mod = fx_submodule(NULL, "naive", "naive"); + naive.Init(naive_mod, centers, bandwidth); + naive.ComputeTotalSum(&naive_results); + + + ErrorAnalysis analysis; + analysis.Init(abs_results, rel_results, exp_results, gauss_results, + hybrid_results, naive_results, abs_mod, rel_mod, exp_mod, + gauss_mod, hybrid_mod, naive_mod); + + analysis.ComputeResults(); + + + + fx_done(); + + return 0; + +} \ No newline at end of file diff --git a/fastlib2/contrib/march/error/hybrid_error.h b/fastlib2/contrib/march/error/hybrid_error.h new file mode 100644 index 0000000000..7bb3cd08ce --- /dev/null +++ b/fastlib2/contrib/march/error/hybrid_error.h @@ -0,0 +1,258 @@ +/** + * @file hybrid_error.h + * + * @author Bill March (march@gatech.edu) + * + * Contains a class for running timing tests for hybrid error schemes. + */ + +#ifndef HYBRID_ERROR_H +#define HYBRID_ERROR_H + +#include "fastlib/fastlib.h" + + +/* + * I need a templated error class. It needs to: + * - return the correct error bound - relative, absolute, hybrid + * - Be able to inform the tree building so I get the right stats for each kind + * of error + * + * It might be possible to just templatize with the Stat class. It can have a + * function to compute the allowed error + * + * TODO: write a naive version too, to evaluate the approximations against + */ + +/** + * Dual-tree implementation of a Gaussian kernel summation with monopole + * error bounds. + * + * TODO: I think this should be templatized to allow for different error styles + * I could templatize with the stat class, which has a function to determine if + * pruning is possible. + */ +template +class GaussianKernelErrorTester { + + //FORBID_ACCIDENTAL_COPY(GaussianKernelErrorTester); +public: + + GaussianKernelErrorTester() {} + + ~GaussianKernelErrorTester() {} + +typedef BinarySpaceTree, Matrix, TErrorStat> ErrorTree; + + private: + + // The common, global bandwidth of the Gaussians + double bandwidth_; + + // The centers of the Gaussians + Matrix centers_; + + // The result of the sum + Vector results_; + + // The fx module for timing + struct datanode* module_; + + // The total number of points + index_t num_points_; + + // The dimension + index_t dimension_; + + ErrorTree* tree_; + + ArrayList old_from_new_; + + + /** + * Computes the value of the gaussian centered at r at the point referred to + * by q. + * + * I don't think I'll worry about normalization for now. + */ + double ComputeGaussian_(index_t q, index_t r) { + + Vector q_vec; + centers_.MakeColumnVector(q, &q_vec); + + Vector r_vec; + centers_.MakeColumnVector(r, &r_vec); + + double dist = la::DistanceSqEuclidean(q_vec, r_vec); + + return(exp(-bandwidth_ * dist)); + + } // ComputeGaussian_() + + double ComputeGaussian_(double dist_sq) { + + return (exp(-bandwidth_ * dist_sq)); + + } // ComputeGaussian_ + + // This function needs to decrease the stat's query_count, since I don't need + // to allocate any more error to those query points + void ComputeSumBaseCase_(ErrorTree* query, ErrorTree* reference) { + + for (index_t query_index = query->begin(); query_index < query->end(); + query_index++) { + + double query_value = results_[query_index]; + + for (index_t ref_index = reference->begin(); ref_index < reference->end(); + ref_index++) { + + query_value = query_value + ComputeGaussian_(query_index, ref_index); + + } + + DEBUG_ASSERT(query_value >= 0.0); + + results_[query_index] = query_value; + + } + + index_t query_count = query->stat().query_count(); + index_t reference_count = reference->count(); + query->stat().set_query_count(query_count - reference_count); + + } // ComputeSumBaseCase_() + + void ComputeSumRecursion_(ErrorTree* query, ErrorTree* reference) { + + double q_min_dist = query->bound().MinDistanceSq(reference->bound()); + double q_max_dist = query->bound().MaxDistanceSq(reference->bound()); + + double q_upper_bound = ComputeGaussian_(q_min_dist) * reference->count(); + double q_lower_bound = ComputeGaussian_(q_max_dist) * reference->count(); + + if (query->is_leaf() && reference->is_leaf()) { + + ComputeSumBaseCase_(query, reference); + + } // Base case + // Maybe I should put the prune check before the base case + else if(query->stat().CanPrune( + q_upper_bound, q_lower_bound, reference->count())) { + + double approximate_result = 0.5 * (q_upper_bound + q_lower_bound); + DEBUG_ASSERT(approximate_result >= 0.0); + + Vector subvec; + results_.MakeSubvector(query->begin(), query->count(), &subvec); + + Vector approx; + approx.Init(query->count()); + approx.SetAll(approximate_result); + // I'm pretty sure this will work, but I should check + la::AddOverwrite(subvec, approx, &subvec); + + } // Pruning case + else if(query->is_leaf()) { + + ComputeSumRecursion_(query, reference->left()); + ComputeSumRecursion_(query, reference->right()); + + } // only split references + else if(reference->is_leaf()) { + + query->left()->stat().set_query_count(query->stat().query_count()); + query->right()->stat().set_query_count(query->stat().query_count()); + + ComputeSumRecursion_(query->left(), reference); + ComputeSumRecursion_(query->right(), reference); + + index_t left_count = query->left()->stat().query_count(); + // This gives an unused variable warning in fast mode + index_t right_count = query->right()->stat().query_count(); + DEBUG_ASSERT(left_count == right_count); + + query->stat().set_query_count(left_count); + + } // only split queries + else { + + // Should consider some kind of priority here + // which side could use the error more effectively? + + query->left()->stat().set_query_count(query->stat().query_count()); + query->right()->stat().set_query_count(query->stat().query_count()); + + ComputeSumRecursion_(query->left(), reference->left()); + ComputeSumRecursion_(query->left(), reference->right()); + + ComputeSumRecursion_(query->right(), reference->left()); + ComputeSumRecursion_(query->right(), reference->right()); + + index_t left_count = query->left()->stat().query_count(); + index_t right_count = query->right()->stat().query_count(); + DEBUG_ASSERT(left_count == right_count); + + query->stat().set_query_count(left_count); + + } // four-way + + } // ComputeSumRecursion_ + +public: + + void InitStats(ErrorTree* node) { + + if (!(node->is_leaf())) { + + InitStats(node->left()); + InitStats(node->right()); + + } + + node->stat().SetParams(module_); + + } // InitStats() + + + void Init(struct datanode* mod, const Matrix& cent, double band) { + + module_ = mod; + + centers_.Copy(cent); + + bandwidth_ = band; + DEBUG_ASSERT(bandwidth_ > 0.0); + + num_points_ = centers_.n_cols(); + + dimension_ = centers_.n_rows(); + + results_.Init(num_points_); + results_.SetZero(); + + tree_ = tree::MakeKdTreeMidpoint(centers_, + fx_param_int(mod, "leaf_size", 20), &old_from_new_, NULL); + + tree_->stat().set_query_count(num_points_); + + InitStats(tree_); + + } // Init() + + void ComputeTotalSum(Vector* results_vec) { + + fx_timer_start(module_, "timer"); + ComputeSumRecursion_(tree_, tree_); + fx_timer_stop(module_, "timer"); + + results_vec->Copy(results_); + + } // ComputeTotalSum() + + +}; // class GaussianKernelErrorTester + + + +#endif \ No newline at end of file diff --git a/fastlib2/contrib/march/error/hybrid_error_analysis.h b/fastlib2/contrib/march/error/hybrid_error_analysis.h new file mode 100644 index 0000000000..65d5be8d00 --- /dev/null +++ b/fastlib2/contrib/march/error/hybrid_error_analysis.h @@ -0,0 +1,122 @@ +#ifndef HYBRID_ERROR_ANALYSIS_H +#define HYBRID_ERROR_ANALYSIS_H + +#include "fastlib/fastlib.h" + +class ErrorAnalysis { + +private: + + Vector abs_vec; + Vector rel_vec; + Vector exp_vec; + Vector gauss_vec; + Vector hybrid_vec; + Vector naive_vec; + + struct datanode* abs_mod_; + struct datanode* rel_mod_; + struct datanode* exp_mod_; + struct datanode* gauss_mod_; + struct datanode* hybrid_mod_; + struct datanode* naive_mod_; + + index_t num_points_; + + void TotalError_() { + + double abs_error = 0.0; + double rel_error = 0.0; + double exp_error = 0.0; + double gauss_error = 0.0; + double hybrid_error = 0.0; + + double abs_rel = 0.0; + double rel_rel = 0.0; + double exp_rel = 0.0; + double gauss_rel = 0.0; + double hybrid_rel = 0.0; + + for (index_t i = 0; i < num_points_; i++) { + + double naive_val = naive_vec[i]; + + double this_abs = fabs(naive_val - abs_vec[i]); + double this_rel = fabs(naive_val - rel_vec[i]); + double this_exp = fabs(naive_val - exp_vec[i]); + double this_gauss = fabs(naive_val - gauss_vec[i]); + double this_hybrid = fabs(naive_val - hybrid_vec[i]); + + abs_error = abs_error + this_abs; + rel_error = rel_error + this_rel; + exp_error = exp_error + this_exp; + gauss_error = gauss_error + this_gauss; + hybrid_error = hybrid_error + this_hybrid; + + abs_rel = abs_rel + (this_abs/naive_val); + rel_rel = rel_rel + (this_rel/naive_val); + exp_rel = exp_rel + (this_exp/naive_val); + gauss_rel = gauss_rel + (this_gauss/naive_val); + hybrid_rel = hybrid_rel + (this_gauss/naive_val); + + } // i + + fx_format_result(abs_mod_, "total_absolute_error", "%g", abs_error); + fx_format_result(rel_mod_, "total_absolute_error", "%g", rel_error); + fx_format_result(exp_mod_, "total_absolute_error", "%g", exp_error); + fx_format_result(gauss_mod_, "total_absolute_error", "%g", gauss_error); + fx_format_result(hybrid_mod_, "total_absolute_error", "%g", hybrid_error); + + fx_format_result(abs_mod_, "total_relative_error", "%g", abs_rel); + fx_format_result(rel_mod_, "total_relative_error", "%g", rel_rel); + fx_format_result(exp_mod_, "total_relative_error", "%g", exp_rel); + fx_format_result(gauss_mod_, "total_relative_error", "%g", gauss_rel); + fx_format_result(hybrid_mod_, "total_relative_error", "%g", hybrid_rel); + + + } // TotalAbsoluteError_() + + + public: + + ErrorAnalysis() {} + + ~ErrorAnalysis() {} + + void Init(const Vector& abs, const Vector& rel, const Vector& exp, + const Vector& gauss, const Vector& hybrid, const Vector& naive, + struct datanode* abs_m, struct datanode* rel_m, + struct datanode* exp_m, struct datanode* gauss_m, + struct datanode* hybrid_m, struct datanode* naive_m) { + + abs_vec.Copy(abs); + rel_vec.Copy(rel); + exp_vec.Copy(exp); + gauss_vec.Copy(gauss); + hybrid_vec.Copy(hybrid); + naive_vec.Copy(naive); + + num_points_ = abs_vec.length(); + + abs_mod_ = abs_m; + rel_mod_ = rel_m; + exp_mod_ = exp_m; + gauss_mod_ = gauss_m; + hybrid_mod_ = hybrid_m; + naive_mod_ = naive_m; + + } // Init() + + + void ComputeResults() { + + TotalError_(); + + } // ComputeResults() + + + +}; // class ErrorAnalysis + + +#endif \ No newline at end of file diff --git a/fastlib2/contrib/march/error/hybrid_error_stat.h b/fastlib2/contrib/march/error/hybrid_error_stat.h new file mode 100644 index 0000000000..f05e6004fa --- /dev/null +++ b/fastlib2/contrib/march/error/hybrid_error_stat.h @@ -0,0 +1,302 @@ +/** + * @file hybrid_error_stat.h + * + * @author Bill March (march@gatech.edu) + * + * Defines the stat classes for the different kinds of error + */ + +#ifndef HYBRID_ERROR_STAT_H +#define HYBRID_ERROR_STAT_H + +// I should make all of these inherit the basic stat stuff +// The classes should be able to overlap entirely except for the Epsilon_ +// function and some private variables for the complicated functions + +class GenericErrorStat { + + protected: + + index_t query_count_; + + double epsilon_; + + virtual double Epsilon_(double upper_bound, double lower_bound) = 0; + +public: + + GenericErrorStat() {} + + virtual ~GenericErrorStat() {} + + void Init(const Matrix& matrix, index_t start, index_t count) { + + query_count_ = count; + + } // Init() (leaves) + + void Init(const Matrix& matrix, index_t start, index_t count, + const GenericErrorStat& left, const GenericErrorStat& right) { + + query_count_ = count; + + } // Init() (non-leaves) + + bool CanPrune(double q_upper_bound, double q_lower_bound, + index_t reference_count) { + + bool prune = false; + + double max_error_incurred = 0.5 * (q_upper_bound - q_lower_bound); + DEBUG_ASSERT(max_error_incurred >= 0.0); + + double allowed_error = q_lower_bound * + Epsilon_(q_upper_bound, q_lower_bound) * reference_count / query_count_; + DEBUG_ASSERT(allowed_error >= 0.0); + + if (max_error_incurred < allowed_error) { + + prune = true; + + epsilon_ = epsilon_ - max_error_incurred; + + query_count_ = query_count_ - reference_count; + DEBUG_ASSERT(query_count_ >= 0.0); + + } + + return prune; + + } // CanPrune() + + void set_query_count(index_t new_count) { + + query_count_ = new_count; + + DEBUG_ASSERT(query_count_ >= 0); + + } // set_query_count() + + index_t query_count() { + + return query_count_; + + } // query_count() + +}; // GenericErrorStat + +/** + * Prunes with an absolute error criterion + */ +class AbsoluteErrorStat : public GenericErrorStat { + + + + protected: + + /** + * Returns the error tolerance as a function of the bounds on Q. In this case + * we divide by the lower bound to get absolute error. + */ + double Epsilon_(double upper_bound, double lower_bound) { + + double eps = epsilon_ / lower_bound; + DEBUG_ASSERT(eps >= 0.0); + + return (eps); + + } // Epsilon_ + + + public: + + AbsoluteErrorStat() {} + + ~AbsoluteErrorStat() {} + + void SetParams(struct datanode* mod) { + + epsilon_ = fx_param_double_req(mod, "epsilon"); + + } // SetParams() + +}; // class AbsoluteErrorStat + + +/** + * Prunes with a relative error criterion + */ +class RelativeErrorStat : public GenericErrorStat { + + private: + + index_t query_count_; + +protected: + /** + * Relative error just depends on epsilon_ + */ + double Epsilon_(double upper_bound, double lower_bound) { + + DEBUG_ASSERT(epsilon_ >= 0.0); + return epsilon_; + + } // Epsilon_ + + public: + + RelativeErrorStat() {} + + ~RelativeErrorStat() {} + + void SetParams(struct datanode* mod) { + + epsilon_ = fx_param_double_req(mod, "epsilon"); + DEBUG_ASSERT(epsilon_ >= 0.0); + + + } // SetParams + + +}; // class RelativeErrorStat + + +/** + * Prunes with the hybrid exponential error criterion + */ +class ExponentialErrorStat : public GenericErrorStat { + + private: + + double max_error_; + + double steepness_; + + double min_error_; + + +protected: + /** + * Hybrid error using the exponential criterion + */ + double Epsilon_(double upper_bound, double lower_bound) { + + double eps = (max_error_ * exp(-steepness_ * upper_bound)) + min_error_ + + epsilon_; + + DEBUG_ASSERT(eps >= 0.0); + + return (eps); + + } // Epsilon_() + + public: + + ExponentialErrorStat() {} + + ~ExponentialErrorStat() {} + + void SetParams(struct datanode* mod) { + + max_error_ = fx_param_double_req(mod, "max_error"); + + steepness_ = fx_param_double_req(mod, "steepness"); + + min_error_ = fx_param_double_req(mod, "min_error"); + + epsilon_ = 0.0; + + } // SetParams() + + +}; // class ExponentialErrorStat + +/** + * Uses a Gaussian hybrid error criterion + */ +class GaussianErrorStat : public GenericErrorStat { + +private: + + double max_error_; + + double steepness_; + + double min_error_; + + +protected: + /** + * Hybrid error using the gaussian criterion + */ + double Epsilon_(double upper_bound, double lower_bound) { + + double eps = (max_error_ * exp(-steepness_ * upper_bound * upper_bound)) + + min_error_ + epsilon_; + + DEBUG_ASSERT(eps >= 0.0); + + return (eps); + + } // Epsilon_() + + public: + + GaussianErrorStat() {} + + ~GaussianErrorStat() {} + + void SetParams(struct datanode* mod) { + + max_error_ = fx_param_double_req(mod, "max_error"); + + steepness_ = fx_param_double_req(mod, "steepness"); + + min_error_ = fx_param_double_req(mod, "min_error"); + + epsilon_ = 0.0; + + } // SetParams + +}; // class GaussianErrorStat + +class HybridErrorStat : public GenericErrorStat { + +private: + + double steepness_; + + +protected: + + double Epsilon_(double upper_bound, double lower_bound) { + + double eps = (1 - exp(-steepness_ * lower_bound)) * epsilon_; + + eps = eps + (exp(-steepness_ * upper_bound) * epsilon_ / lower_bound); + + DEBUG_ASSERT(eps >= 0.0); + + return eps; + + } // Epsilon_() + +public: + + HybridErrorStat() {} + + ~HybridErrorStat() {} + + void SetParams(struct datanode* mod) { + + steepness_ = fx_param_double_req(mod, "steepness"); + + epsilon_ = fx_param_double_req(mod, "epsilon"); + + } // SetParams() + +}; // class HybridErrorStat + + + +#endif \ No newline at end of file diff --git a/fastlib2/contrib/march/error/naive_kernel_sum.h b/fastlib2/contrib/march/error/naive_kernel_sum.h new file mode 100644 index 0000000000..5bb24d404e --- /dev/null +++ b/fastlib2/contrib/march/error/naive_kernel_sum.h @@ -0,0 +1,92 @@ +#ifndef NAIVE_KERNEL_SUM_H +#define NAIVE_KERNEL_SUM_H + + +#include "fastlib/fastlib.h" + + +class NaiveKernelSum { + + private: + + Vector results_; + + Matrix centers_; + + double bandwidth_; + + index_t num_points_; + + index_t dimension_; + + struct datanode* module_; + + double ComputeGaussian_(index_t i, index_t j) { + + Vector i_vec; + centers_.MakeColumnVector(i, &i_vec); + + Vector j_vec; + centers_.MakeColumnVector(j, &j_vec); + + double dist_sq = la::DistanceSqEuclidean(i_vec, j_vec); + + return (exp(-bandwidth_ * dist_sq)); + + } // ComputeGaussian_() + + public: + + NaiveKernelSum() {} + + ~NaiveKernelSum() {} + + void Init(struct datanode* mod, const Matrix& cent, double band) { + + centers_.Copy(cent); + + bandwidth_ = band; + DEBUG_ASSERT(bandwidth_ > 0.0); + + num_points_ = centers_.n_cols(); + + dimension_ = centers_.n_rows(); + + results_.Init(num_points_); + results_.SetZero(); + + module_ = mod; + + } // Init() + + void ComputeTotalSum(Vector* return_results) { + + fx_timer_start(module_, "timing"); + + for (index_t i = 0; i < num_points_; i++) { + + double this_result = results_[i]; + + for (index_t j = i; j < num_points_; j++) { // for symmetry + + double this_kernel = ComputeGaussian_(i, j); + this_result = this_result + this_kernel; + + results_[j] = results_[j] + this_kernel; // take advantage of symmetry + + } // j + + results_[i] = this_result; + + } // i + + fx_timer_stop(module_, "timing"); + + return_results->Copy(results_); + + } // ComputeTotalSum() + +}; // NaiveKernelSum + + +#endif