error testing stuff
This commit is contained in:
Bill March
@@ -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<AbsoluteErrorStat> absolute;
|
||||
GaussianKernelErrorTester<RelativeErrorStat> relative;
|
||||
GaussianKernelErrorTester<ExponentialErrorStat> exponential;
|
||||
GaussianKernelErrorTester<GaussianErrorStat> gaussian;
|
||||
GaussianKernelErrorTester<HybridErrorStat> 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;
|
||||
|
||||
}
|
||||
@@ -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 <typename TErrorStat>
|
||||
class GaussianKernelErrorTester {
|
||||
|
||||
//FORBID_ACCIDENTAL_COPY(GaussianKernelErrorTester);
|
||||
public:
|
||||
|
||||
GaussianKernelErrorTester() {}
|
||||
|
||||
~GaussianKernelErrorTester() {}
|
||||
|
||||
typedef BinarySpaceTree<DHrectBound<2>, 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<index_t> 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<ErrorTree>(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
|
||||
@@ -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
|
||||
@@ -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
|
||||
@@ -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
|
||||
Reference in New Issue
Block a user