diff --git a/src/contrib/nslagle/CMakeLists.txt b/src/contrib/nslagle/CMakeLists.txt index a79986dea2..3b028885a9 100644 --- a/src/contrib/nslagle/CMakeLists.txt +++ b/src/contrib/nslagle/CMakeLists.txt @@ -3,9 +3,9 @@ cmake_minimum_required(VERSION 2.8) # we just want to recurse into the child directories here set(DIRS myKDE - kdeDL - nested_summation_template - proximity_project + #kdeDL + #nested_summation_template + #proximity_project ) foreach(dir ${DIRS}) diff --git a/src/contrib/nslagle/kdeDL/dualtree_kde.h b/src/contrib/nslagle/kdeDL/dualtree_kde.h index 6491909cbc..9ae77314c3 100644 --- a/src/contrib/nslagle/kdeDL/dualtree_kde.h +++ b/src/contrib/nslagle/kdeDL/dualtree_kde.h @@ -65,7 +65,7 @@ #define INSIDE_DUALTREE_KDE_H #include "mlpack/core.h" -#include "mlpack/core/tree/spacetree.hpp" +#include "mlpack/core/tree/binary_space_tree.hpp" //#include "contrib/nslagle/series_expansion/farfield_expansion.h" //#include "contrib/nslagle/series_expansion/local_expansion.h" //#include "contrib/nslagle/series_expansion/mult_farfield_expansion.h" @@ -248,7 +248,7 @@ class DualtreeKde { /** @brief The permutation mapping indices of references_ to * original order. */ - arma::Col old_from_new_references_; + std::vector old_from_new_references_; ////////// Private Member Functions ////////// @@ -342,13 +342,13 @@ class DualtreeKde { // } // else { // NOTIFY("Using the default dimension of %d", qset_.n_rows); - mult_const_ = 1.0 / ka_.kernel_.CalcNormConstant(qset_.n_rows); + mult_const_ = 1.0;// / ka_.kernel_.CalcNormConstant(qset_.n_rows); // } // Set accuracy parameters. relative_error_ = CLI::GetParam("relative_error"); - threshold_ = CLI::GetParam("threshold") * - ka_.kernel_.CalcNormConstant(qset_.n_rows); + threshold_ = CLI::GetParam("threshold"); + // * ka_.kernel_.CalcNormConstant(qset_.n_rows); // initialize the lower and upper bound densities densities_l_.zeros(); @@ -365,7 +365,7 @@ class DualtreeKde { num_local_prunes_ = 0; printf("\nStarting fast KDE on bandwidth value of %g...\n", - sqrt(ka_.kernel_.bandwidth_sq())); + sqrt(ka_.Bandwidth())); CLI::StartTimer ("fast_kde_compute"); // Preprocessing step for initializing series expansion objects @@ -417,7 +417,7 @@ class DualtreeKde { (&queries == &references); // Read in the number of points owned by a leaf. - int leaflen = CLI::GetParam("leaflen"); + // TODO :int leaflen = CLI::GetParam("leaflen"); // Copy reference dataset and reference weights and compute its // sum. @@ -455,9 +455,7 @@ class DualtreeKde { // weights according to the permutation of the reference set in // the reference tree. CLI::StartTimer("tree_d"); - rroot_ = proximity::MakeGenMetricTree(rset_, leaflen, - &old_from_new_references_, - NULL); + rroot_ = new Tree(rset_, /*leaflen,*/old_from_new_references_); DualtreeKdeCommon::ShuffleAccordingToPermutation (rset_weights_, old_from_new_references_); @@ -466,12 +464,10 @@ class DualtreeKde { old_from_new_queries_ = old_from_new_references_; } else { - qroot_ = proximity::MakeGenMetricTree(qset_, leaflen, - &old_from_new_queries_, - NULL); + qroot_ = new Tree(qset_, /*leaflen,*/ old_from_new_queries_); } CLI::StopTimer("tree_d"); - + // Initialize the density lists densities_l_ = arma::vec(qset_.n_cols); densities_e_ = arma::vec(qset_.n_cols); @@ -497,35 +493,35 @@ class DualtreeKde { { order = 7; } - ka_.Init(bandwidth, order, qset_.n_rows); + ka_ = kernel::GaussianKernel(bandwidth);// order, qset_.n_rows); } else if(qset_.n_rows <= 3) { if (!hasOrder) { order = 5; } - ka_.Init(bandwidth, order, qset_.n_rows); + ka_ = kernel::GaussianKernel(bandwidth);//, order, qset_.n_rows); } else if(qset_.n_rows <= 5) { if (!hasOrder) { order = 3; } - ka_.Init(bandwidth, order, qset_.n_rows); + ka_ = kernel::GaussianKernel(bandwidth);//, order, qset_.n_rows); } else if(qset_.n_rows <= 6) { if (!hasOrder) { order = 1; } - ka_.Init(bandwidth, order, qset_.n_rows); + ka_ = kernel::GaussianKernel(bandwidth);//, order, qset_.n_rows); } else { if (!hasOrder) { order = 0; } - ka_.Init(bandwidth, order, qset_.n_rows); + ka_ = kernel::GaussianKernel(bandwidth);//, order, qset_.n_rows); } } @@ -534,7 +530,7 @@ class DualtreeKde { FILE *stream = stdout; const char *fname = NULL; - if((fname = CLI::GetParam("fast_kde_output")) != NULL) { + if((fname = CLI::GetParam("fast_kde_output").c_str()) != NULL) { stream = fopen(fname, "w+"); } for(size_t q = 0; q < qset_.n_cols; q++) { diff --git a/src/contrib/nslagle/kdeDL/dualtree_kde_common.h b/src/contrib/nslagle/kdeDL/dualtree_kde_common.h index 36c034fab1..cc852343ed 100644 --- a/src/contrib/nslagle/kdeDL/dualtree_kde_common.h +++ b/src/contrib/nslagle/kdeDL/dualtree_kde_common.h @@ -97,7 +97,7 @@ class DualtreeKdeCommon { * @param permutation The permutation. */ static void ShuffleAccordingToPermutation - (arma::vec &v, const arma::Col &permutation) { + (arma::vec &v, const std::vector &permutation) { arma::vec v_tmp(v.size()); for(size_t i = 0; i < v_tmp.size(); i++) { diff --git a/src/contrib/nslagle/kdeDL/dualtree_kde_impl.h b/src/contrib/nslagle/kdeDL/dualtree_kde_impl.h index b51d03e4ce..efa240c876 100644 --- a/src/contrib/nslagle/kdeDL/dualtree_kde_impl.h +++ b/src/contrib/nslagle/kdeDL/dualtree_kde_impl.h @@ -27,7 +27,7 @@ void DualtreeKde::DualtreeKdeBase_(Tree *qnode, Tree *rnode, // pairwise distance and kernel value double dsqd = kernel::LMetric<2,false>::Evaluate (q_col, r_col); - double kernel_value = ka_.kernel_.EvalUnnormOnSq(dsqd); + double kernel_value = ka_.kernel_.Evaluate(q_col,r_col);//EvalUnnormOnSq(dsqd); double weighted_kernel_value = rset_weights_[r] * kernel_value; densities_l_[q] += weighted_kernel_value; diff --git a/src/contrib/nslagle/kdeDL/dualtree_kde_main.cc b/src/contrib/nslagle/kdeDL/dualtree_kde_main.cc index 76bffe4304..a1f665cf0d 100644 --- a/src/contrib/nslagle/kdeDL/dualtree_kde_main.cc +++ b/src/contrib/nslagle/kdeDL/dualtree_kde_main.cc @@ -12,6 +12,7 @@ #include "dualtree_vkde.h" #include "naive_kde.h" +using namespace mlpack; using namespace mlpack::kernel; void VariableBandwidthKde(arma::mat &queries, arma::mat &references, @@ -267,18 +268,18 @@ int main(int argc, char *argv[]) { !strcmp(queries_file_name.c_str(), references_file_name.c_str()); // data::Load inits a arma::mat with the contents of a .csv or .arff. - data::Load(references_file_name.c_str(), references); + references.load(references_file_name.c_str()); if(queries_equal_references) { queries = references; } else { - data::Load(queries_file_name.c_str(), queries); + queries.load(queries_file_name.c_str()); } // If the reference weight file name is specified, then read in, // otherwise, initialize to uniform weights. if(CLI::HasParam("dwgts")) { - data::Load(CLI::GetParam("dwgts").c_str(), reference_weights); + reference_weights.load(CLI::GetParam("dwgts").c_str()); } else { reference_weights = arma::mat(1, references.n_cols); diff --git a/src/contrib/nslagle/kdeDL/dualtree_vkde.h b/src/contrib/nslagle/kdeDL/dualtree_vkde.h index d590d16913..58c7aa1f07 100644 --- a/src/contrib/nslagle/kdeDL/dualtree_vkde.h +++ b/src/contrib/nslagle/kdeDL/dualtree_vkde.h @@ -61,9 +61,9 @@ #define INSIDE_DUALTREE_VKDE_H #include "mlpack/core.h" -#include "mlpack/core/tree/spacetree.hpp" +#include "mlpack/core/tree/binary_space_tree.hpp" //#include "contrib/dongryel/proximity_project/gen_metric_tree.h" -//#include "dualtree_kde_common.h" +#include "dualtree_kde_common.h" //#include "kde_stat.h" #include "mlpack/methods/neighbor_search/neighbor_search.h" @@ -196,12 +196,12 @@ class DualtreeVKde { /** @brief The permutation mapping indices of queries_ to original * order. */ - arma::Col old_from_new_queries_; + std::vector old_from_new_queries_; /** @brief The permutation mapping indices of references_ to * original order. */ - arma::Col old_from_new_references_; + std::vector old_from_new_references_; ////////// Private Member Functions ////////// @@ -287,7 +287,7 @@ class DualtreeVKde { // Set accuracy parameters. relative_error_ = CLI::GetParam("relative_error"); threshold_ = CLI::GetParam("threshold") * - kernels_[0].CalcNormConstant(qset_.n_rows()); + kernels_[0].Normalizer();//CalcNormConstant(qset_.n_rows); // initialize the lower and upper bound densities densities_l_.zeros(); @@ -352,7 +352,7 @@ class DualtreeVKde { (&queries == &references); // read in the number of points owned by a leaf - int leaflen = CLI::GetParam("leaflen"); + //TODO int leaflen = CLI::GetParam("leaflen"); // Copy reference dataset and reference weights and compute its // sum. rset_weight_sum_ should be the raw sum of the reference @@ -392,9 +392,7 @@ class DualtreeVKde { // weights according to the permutation of the reference set in // the reference tree. CLI::StartTimer("tree_d"); - rroot_ = proximity::MakeGenMetricTree(rset_, leaflen, - old_from_new_references_, - NULL); + rroot_ = new Tree(rset_, /*leaflen,*/ old_from_new_references_); DualtreeKdeCommon::ShuffleAccordingToPermutation (rset_weights_, old_from_new_references_); @@ -403,9 +401,7 @@ class DualtreeVKde { old_from_new_queries_ = old_from_new_references_; } else { - qroot_ = proximity::MakeGenMetricTree(qset_, leaflen, - &old_from_new_queries_, - NULL); + qroot_ = new Tree(qset_, /*leaflen*/ old_from_new_queries_); } CLI::StopTimer("tree_d"); @@ -420,16 +416,16 @@ class DualtreeVKde { // Initialize the kernels for each reference point. int knns = CLI::GetParam("knn"); - AllkNN all_knn = AllkNN(rset_, 20); - kernels_.Init(rset_.n_cols()); + AllkNN all_knn = AllkNN(rset_, knns); arma::Mat resulting_neighbors; - arma::mat squared_distances; + arma::mat squared_distances; CLI::StartTimer("bandwidth_initialization"); all_knn.ComputeNeighbors(resulting_neighbors, squared_distances); - for(size_t i = 0; i < squared_distances.size(); i += knns) { - kernels_[i / knns].Init(sqrt(squared_distances[i + knns - 1])); + for(size_t i = 0; i < squared_distances.size(); i += knns) + { + kernels_.push_back(kernel::GaussianKernel(sqrt(squared_distances[i + knns - 1]))); } CLI::StopTimer("bandwidth_initialization"); @@ -437,11 +433,11 @@ class DualtreeVKde { // that have been chosen. double min_norm_const = DBL_MAX; for(size_t i = 0; i < rset_weights_.size(); i++) { - double norm_const = kernels_[i].CalcNormConstant(qset_.n_rows()); + double norm_const = kernels_[i].Normalizer();//CalcNormConstant(qset_.n_rows); min_norm_const = std::min(min_norm_const, norm_const); } for(size_t i = 0; i < rset_weights_.size(); i++) { - double norm_const = kernels_[i].CalcNormConstant(qset_.n_rows()); + double norm_const = kernels_[i].Normalizer();//CalcNormConstant(qset_.n_rows); rset_weights_[i] *= (min_norm_const / norm_const); } @@ -454,11 +450,11 @@ class DualtreeVKde { FILE *stream = stdout; const char *fname = NULL; - if((fname = CLI::GetParam("fast_kde_output")) != NULL) + if((fname = CLI::GetParam("fast_kde_output").c_str()) != NULL) { stream = fopen(fname, "w+"); } - for(size_t q = 0; q < qset_.n_cols(); q++) { + for(size_t q = 0; q < qset_.n_cols; q++) { fprintf(stream, "%g\n", densities_e_[q]); } diff --git a/src/contrib/nslagle/kdeDL/kde_stat.h b/src/contrib/nslagle/kdeDL/kde_stat.h index b213972fbc..9eac478237 100644 --- a/src/contrib/nslagle/kdeDL/kde_stat.h +++ b/src/contrib/nslagle/kdeDL/kde_stat.h @@ -200,12 +200,12 @@ class KdeStat { /** @brief The far field expansion created by the reference points * in this node. */ - typename TKernelAux::TFarFieldExpansion farfield_expansion_; - + // TODO: typename TKernelAux::TFarFieldExpansion farfield_expansion_; + /** @brief The local expansion stored in this node. */ - typename TKernelAux::TLocalExpansion local_expansion_; - + // TODO :typename TKernelAux::TLocalExpansion local_expansion_; + /** @brief The subspace associated with this node. */ //SubspaceStat subspace_; @@ -213,9 +213,9 @@ class KdeStat { /** @brief Gets the weight sum. */ double get_weight_sum() { - return farfield_expansion_.get_weight_sum(); + return 0.0;//farfield_expansion_.get_weight_sum(); } - + /** @brief Adds the other postponed contributions. */ void AddPostponed(const KdeStat& parent_stat) { @@ -281,8 +281,8 @@ class KdeStat { } void Init(const TKernelAux &ka) { - farfield_expansion_.Init(ka); - local_expansion_.Init(ka); + // TODO: farfield_expansion_.Init(ka); + // TODO: local_expansion_.Init(ka); } void Init(const arma::mat& dataset, size_t &start, size_t &count) { @@ -300,8 +300,8 @@ class KdeStat { void Init(const arma::vec& center, const TKernelAux &ka) { - farfield_expansion_.Init(center, ka); - local_expansion_.Init(center, ka); + // TODO: farfield_expansion_.Init(center, ka); + // TODO: local_expansion_.Init(center, ka); Init(); } diff --git a/src/contrib/nslagle/kdeDL/naive_kde.h b/src/contrib/nslagle/kdeDL/naive_kde.h index 8993e6ee84..ae5eaac553 100644 --- a/src/contrib/nslagle/kdeDL/naive_kde.h +++ b/src/contrib/nslagle/kdeDL/naive_kde.h @@ -98,16 +98,16 @@ class NaiveKde { printf("\nStarting naive KDE...\n"); CLI::StartTimer("naive_kde_compute"); - for(size_t q = 0; q < qset_.n_cols(); q++) { + for(size_t q = 0; q < qset_.n_cols; q++) { const arma::vec q_col = qset_.unsafe_col(q); // Compute unnormalized sum first. - for(size_t r = 0; r < rset_.n_cols(); r++) { + for(size_t r = 0; r < rset_.n_cols; r++) { const arma::vec r_col = rset_.unsafe_col(r); double dsqd = kernel::LMetric<2,false>::Evaluate(q_col, r_col); - densities_[q] += rset_weights_[r] * kernels_[r].EvalUnnormOnSq(dsqd); + densities_[q] += rset_weights_[r] * kernels_[r].Evaluate(q_col,r_col);//EvalUnnormOnSq(dsqd); } // Then normalize it. @@ -127,16 +127,16 @@ class NaiveKde { printf("\nStarting naive KDE...\n"); CLI::StartTimer("naive_kde_compute"); - for(size_t q = 0; q < qset_.n_cols(); q++) { + for(size_t q = 0; q < qset_.n_cols; q++) { const arma::vec q_col = qset_.unsafe_col(q); // Compute unnormalized sum. - for(size_t r = 0; r < rset_.n_cols(); r++) { + for(size_t r = 0; r < rset_.n_cols; r++) { const arma::vec r_col = rset_.unsafe_col(r); - double dsqd = kernel::LMetric<2,false>::Evaluate (q_col, r_col); + //double dsqd = kernel::LMetric<2,false>::Evaluate (q_col, r_col); - densities_[q] += rset_weights_[r] * kernels_[r].EvalUnnormOnSq(dsqd); + densities_[q] += rset_weights_[r] * kernels_[r].Evaluate(q_col,r_col);//EvalUnnormOnSq(dsqd); } // Then, normalize it. densities_[q] /= norm_const_; @@ -148,7 +148,7 @@ class NaiveKde { void Init(arma::mat &qset, arma::mat &rset, struct datanode *module_in) { // Use the uniform weights for a moment. - arma::mat uniform_weights(1, rset.n_cols()); + arma::mat uniform_weights(1, rset.n_cols); uniform_weights.fill(1.0); Init(qset, rset, uniform_weights, module_in); @@ -184,7 +184,7 @@ class NaiveKde { rset_(r,c) = rset(r,c); } } - rset_weights_ = arma::vec(reference_weights.n_cols()); + rset_weights_ = arma::vec(reference_weights.n_cols); for(size_t i = 0; i < rset_weights_.size(); i++) { rset_weights_[i] = reference_weights(0, i); @@ -197,20 +197,20 @@ class NaiveKde { } // Get bandwidth and compute the normalizing constant. - kernels_.Init(rset_.n_cols()); - if(!strcmp(CLI::GetParam("mode").c_str(), "variablebw")) { - + if(!strcmp(CLI::GetParam("mode").c_str(), "variablebw")) + { // Initialize the kernels for each reference point. int knns = CLI::GetParam("knn"); - AllkNN all_knn = AllkNN(rset_, 20, knns); + AllkNN all_knn = AllkNN(rset_, knns); arma::Mat resulting_neighbors; arma::mat squared_distances; CLI::StartTimer("bandwidth_initialization"); all_knn.ComputeNeighbors(resulting_neighbors, squared_distances); - for(size_t i = 0; i < squared_distances.size(); i += knns) { - kernels_[i / knns].Init(sqrt(squared_distances[i + knns - 1])); + for(size_t i = 0; i < squared_distances.size(); i += knns) + { + kernels_.push_back(kernel::GaussianKernel(sqrt(squared_distances[i + knns - 1]))); } CLI::StopTimer("bandwidth_initialization"); @@ -218,30 +218,32 @@ class NaiveKde { // that have been chosen. double min_norm_const = DBL_MAX; for(size_t i = 0; i < rset_weights_.size(); i++) { - double norm_const = kernels_[i].CalcNormConstant(qset_.n_rows()); + double norm_const = kernels_[i].Normalizer();//CalcNormConstant(qset_.n_rows()); min_norm_const = std::min(min_norm_const, norm_const); } for(size_t i = 0; i < rset_weights_.size(); i++) { - double norm_const = kernels_[i].CalcNormConstant(qset_.n_rows()); + double norm_const = kernels_[i].Normalizer();//CalcNormConstant(qset_.n_rows()); rset_weights_[i] *= (min_norm_const / norm_const); } - + // Compute normalization constant. norm_const_ = weight_sum * min_norm_const; } - else { - for(size_t i = 0; i < kernels_.size(); i++) { - kernels_[i].Init(CLI::GetParam("bandwidth")); + else + { + for(size_t i = 0; i < kernels_.size(); i++) + { + kernels_[i] = kernel::GaussianKernel(CLI::GetParam("bandwidth")); } - norm_const_ = kernels_[0].CalcNormConstant(qset_.n_rows()) * weight_sum; + norm_const_ = kernels_[0].Normalizer();//CalcNormConstant(qset_.n_rows()) * weight_sum; } // Allocate density storage. - densities_ = arma::vec(qset.n_cols()); + densities_ = arma::vec(qset.n_cols); densities_.zeros(); } - /** @brief Output KDE results to a stream + /** @brief Output KDE results to a stream * * If the user provided "--naive_kde_output=" argument, then the * output will be directed to a file whose name is provided after @@ -254,10 +256,10 @@ class NaiveKde { const char *fname = NULL; { - fname = CLI::GetParam("naive_kde_output"); + fname = CLI::GetParam("naive_kde_output").c_str(); stream = fopen(fname, "w+"); } - for(size_t q = 0; q < qset_.n_cols(); q++) { + for(size_t q = 0; q < qset_.n_cols; q++) { fprintf(stream, "%g\n", densities_[q]); } diff --git a/src/contrib/nslagle/myKDE/CMakeLists.txt b/src/contrib/nslagle/myKDE/CMakeLists.txt index cde89d749e..29cd508287 100644 --- a/src/contrib/nslagle/myKDE/CMakeLists.txt +++ b/src/contrib/nslagle/myKDE/CMakeLists.txt @@ -4,9 +4,9 @@ cmake_minimum_required(VERSION 2.8) # Anything not in this list will not be compiled into the output library # Do not include test programs here set(SOURCES - test_tree.cc kde_dual_tree.hpp kde_dual_tree_impl.hpp + test_kde_dual_tree.cpp ) # add directory name to sources @@ -18,10 +18,10 @@ endforeach() set(MLPACK_CONTRIB_SRCS ${MLPACK_CONTRIB_SRCS} ${DIR_SRCS} PARENT_SCOPE) # link dependencies of test executable -add_executable(test_tree +add_executable(test_kde_dual_tree EXCLUDE_FROM_ALL - test_tree.cc + test_kde_dual_tree.cpp ) -target_link_libraries(test_tree +target_link_libraries(test_kde_dual_tree mlpack ) diff --git a/src/contrib/nslagle/myKDE/kde_dual_tree.hpp b/src/contrib/nslagle/myKDE/kde_dual_tree.hpp index b5c8763b91..bef98f9b8f 100644 --- a/src/contrib/nslagle/myKDE/kde_dual_tree.hpp +++ b/src/contrib/nslagle/myKDE/kde_dual_tree.hpp @@ -2,7 +2,7 @@ #define KDE_DUAL_TREE_HPP #include -#include +#include #include #include @@ -17,6 +17,7 @@ namespace mlpack namespace kde { /* structure within the priority queue */ +template > > struct queueNode { TTree* T; @@ -28,13 +29,14 @@ struct queueNode size_t bLowerIndex; size_t bUpperIndex; }; +template > > class QueueNodeCompare { bool reverse; public: QueueNodeCompare(const bool& revparam=false) : reverse(revparam) {} - bool operator() (const struct queueNode& lhs, - const struct queueNode& rhs) const + bool operator() (const struct queueNode& lhs, + const struct queueNode& rhs) const { if (reverse) return (lhs.priority>rhs.priority); @@ -52,6 +54,8 @@ class KdeDualTree /* possibly, these refer to the same object */ TTree* referenceRoot; TTree* queryRoot; + std::map nodeIndices; + size_t nextAvailableNodeIndex; std::vector referenceShuffledIndices; std::vector queryShuffledIndices; arma::mat referenceData; @@ -66,17 +70,19 @@ class KdeDualTree double delta; /* relative error with respect to the density estimate */ double epsilon; - math::Range bandwidths; - std::priority_queue, - QueueNodeCompare> nodePriorityQueue; + std::priority_queue, + std::vector >, + QueueNodeCompare > nodePriorityQueue; size_t bandwidthCount; std::vector bandwidths; + std::vector inverseBandwidths; + double lowBandwidth; + double highBandwidth; size_t levelsInTree; size_t queryTreeSize; void SetDefaults(); - void MultiBandwidthDualTree(); + size_t MultiBandwidthDualTree(); void MultiBandwidthDualTreeBase(TTree* Q, TTree* T, size_t QIndex, size_t lowerBIndex, size_t upperBIndex); @@ -88,22 +94,24 @@ class KdeDualTree { return levelsInTree - node->levelsBelow(); } - void Winnow(size_t bLower, size_t bUpper, size_t* newLower, size_t* newUpper); + void Winnow(size_t level, size_t* newLower, size_t* newUpper); public: /* the two data sets are different */ KdeDualTree (arma::mat& referenceData, arma::mat& queryData); /* the reference data is also the query data */ KdeDualTree (arma::mat& referenceData); + std::vector Calculate(); /* setters and getters */ - const math::Range& BandwidthRange() const { return bandwidthRange; } const size_t& BandwidthCount() const { return bandwidthCount; } const double& Delta() const { return delta; } const double& Epsilon() const { return epsilon; } + const double& LowBandwidth() const { return lowBandwidth; } + const double& HighBandwidth() const { return highBandwidth; } - void BandwidthRange(double l, double u) { bandwidthRange = math::Range(l,u); } size_t& BandwidthCount() { return bandwidthCount; } double& Delta() { return delta; } double& Epsilon() { return epsilon; } + void SetBandwidthBounds(double l, double u); }; }; /* end namespace kde */ }; /* end namespace mlpack */ diff --git a/src/contrib/nslagle/myKDE/kde_dual_tree_impl.hpp b/src/contrib/nslagle/myKDE/kde_dual_tree_impl.hpp index 7659d1c7fa..8576f92a74 100644 --- a/src/contrib/nslagle/myKDE/kde_dual_tree_impl.hpp +++ b/src/contrib/nslagle/myKDE/kde_dual_tree_impl.hpp @@ -8,6 +8,8 @@ #endif #endif +#define MADEIT std::cout<<"made it to "<<__LINE__<<" in "<<__FILE__< KdeDualTree::KdeDualTree (arma::mat& reference, arma::mat& query) { - referenceRoot (new TTree (reference)), - queryRoot (new TTree (query)) + referenceRoot = new TTree (reference, referenceShuffledIndices), + queryRoot = new TTree (query, queryShuffledIndices); referenceData = reference; queryData = query; levelsInTree = queryRoot->levelsBelow(); @@ -37,27 +39,121 @@ KdeDualTree::KdeDualTree (arma::mat& reference) referenceRoot = new TTree (reference, referenceShuffledIndices); queryRoot = referenceRoot; + queryShuffledIndices = referenceShuffledIndices; levelsInTree = queryRoot->levelsBelow(); queryTreeSize = queryRoot->treeSize(); SetDefaults(); } template -KdeDualTree::SetDefaults() +void KdeDualTree::SetDefaults() { - BandwidthRange(0.01, 100.0); + SetBandwidthBounds(0.01, 100.0); bandwidthCount = 10; delta = epsilon = 0.05; kernel = TKernel(1.0); + nextAvailableNodeIndex = 0; } template -void KdeDualTree::MultiBandwidthDualTree() +std::vector KdeDualTree::Calculate() { + /* calculate the bandwidths */ + bandwidths.clear(); + inverseBandwidths.clear(); + + if (bandwidthCount > 1) + { + double bandwidthDelta = (highBandwidth - lowBandwidth) / (bandwidthCount - 1); + for (size_t bIndex = 0; bIndex < bandwidthCount; ++bIndex) + { + bandwidths.push_back(lowBandwidth + bandwidthDelta * bIndex); + inverseBandwidths.push_back(1.0 / bandwidths.back()); + } + } + else + { + bandwidths.push_back(lowBandwidth); + inverseBandwidths.push_back(1.0 / lowBandwidth); + } + + /* resize the critical matrices */ + upperBoundLevelByBandwidth.zeros(levelsInTree,bandwidthCount); + for (size_t bIndex = 0; bIndex < bandwidthCount; ++bIndex) + { + arma::vec col = upperBoundLevelByBandwidth.unsafe_col(bIndex); + col.fill(referenceRoot->count() * inverseBandwidths[bIndex]); + } + upperBoundLevelByBandwidth.fill(referenceRoot->count()); + lowerBoundLevelByBandwidth.zeros(levelsInTree,bandwidthCount); + upperBoundQPointByBandwidth.zeros(queryRoot->count(),bandwidthCount); + for (size_t bIndex = 0; bIndex < bandwidthCount; ++bIndex) + { + arma::vec col = upperBoundQPointByBandwidth.unsafe_col(bIndex); + col.fill(referenceRoot->count() * inverseBandwidths[bIndex]); + } + lowerBoundQPointByBandwidth.zeros(queryRoot->count(),bandwidthCount); + upperBoundQNodeByBandwidth.zeros(queryTreeSize,bandwidthCount); + for (size_t bIndex = 0; bIndex < bandwidthCount; ++bIndex) + { + arma::vec col = upperBoundQNodeByBandwidth.unsafe_col(bIndex); + col.fill(referenceRoot->count() * inverseBandwidths[bIndex]); + } + lowerBoundQNodeByBandwidth.zeros(queryTreeSize,bandwidthCount); + + arma::vec dl; + arma::vec du; + dl.zeros(bandwidthCount); + du.zeros(bandwidthCount); + double priority = pow( + queryRoot->bound().MinDistance(referenceRoot->bound()), + 0.5); + struct queueNode firstNode = + {referenceRoot,queryRoot, nextAvailableNodeIndex, dl, du, + priority, 0, bandwidthCount - 1}; + nodeIndices[queryRoot] = nextAvailableNodeIndex; + ++nextAvailableNodeIndex; + nodePriorityQueue.push(firstNode); + size_t finalLevel = MultiBandwidthDualTree(); + + size_t maxIndex = 0; + double maxLogLikelihood = (upperBoundLevelByBandwidth(finalLevel,0) + + lowerBoundLevelByBandwidth(finalLevel,0)) / 2.0; + for (size_t bIndex = 1; bIndex < bandwidthCount; ++bIndex) + { + double currentLogLikelihood = (upperBoundLevelByBandwidth(finalLevel,bIndex) + + lowerBoundLevelByBandwidth(finalLevel,bIndex)) / 2.0; + if (currentLogLikelihood > maxLogLikelihood) + { + currentLogLikelihood = maxLogLikelihood; + maxIndex = bIndex; + } + } + std::cout << upperBoundLevelByBandwidth << "\n"; + std::cout << lowerBoundLevelByBandwidth << "\n"; + std::cout << "best bandwidth " << bandwidths[maxIndex] << ";\n"; + exit(1); + std::vector densities; + for (std::vector::iterator shuffIt = queryShuffledIndices.begin(); + shuffIt != queryShuffledIndices.end(); ++shuffIt) + { + densities.push_back((upperBoundQPointByBandwidth(*shuffIt, maxIndex) + + lowerBoundQPointByBandwidth(*shuffIt, maxIndex)) / (2.0 * referenceRoot->count())); + + } + return densities; +} + +template +size_t KdeDualTree::MultiBandwidthDualTree() +{ + /* current level */ + size_t v = 0; while (!nodePriorityQueue.empty()) { /* get the first structure in the queue */ - struct queueNode queueCurrent = nodePriorityQueue.pop(); + struct queueNode queueCurrent = nodePriorityQueue.top(); + nodePriorityQueue.pop(); TTree* Q = queueCurrent.Q; TTree* T = queueCurrent.T; size_t sizeOfTNode = T->count(); @@ -66,7 +162,7 @@ void KdeDualTree::MultiBandwidthDualTree() arma::vec deltaLower = queueCurrent.deltaLower; arma::vec deltaUpper = queueCurrent.deltaUpper; /* v is the level of the Q node */ - size_t v = GetLevelOfNode(Q); + v = GetLevelOfNode(Q); size_t bUpper = queueCurrent.bUpperIndex; size_t bLower = queueCurrent.bLowerIndex; /* check to see whether we've reached the epsilon condition */ @@ -96,7 +192,7 @@ void KdeDualTree::MultiBandwidthDualTree() /* return */ if (epsilonCondition) { - return; + return v; } /* we didn't meet the criteria; let's narrow the bandwidth range */ Winnow(v, &bLower, &bUpper); @@ -109,11 +205,13 @@ void KdeDualTree::MultiBandwidthDualTree() std::vector deltaCondition; for (size_t bIndex = bLower; bIndex <= bUpper; ++bIndex) { - double bandwidth = bandwidths[bIndex]; - double dl = sizeOfTNode * kernel(dMax / bandwidth); - double du = sizeOfTNode * kernel(dMin / bandwidth); + double inverseBandwidth = inverseBandwidths[bIndex]; + double dl = sizeOfTNode * inverseBandwidth * kernel.Evaluate(dMax * inverseBandwidth); + double du = sizeOfTNode * inverseBandwidth * kernel.Evaluate(dMin * inverseBandwidth); deltaLower(bIndex) = dl; deltaUpper(bIndex) = du - sizeOfTNode; + //std::cout << "QIndex: " << QIndex << " bIndex: " << bIndex << std::endl; + //std::cout << "max QIndex: " << queryTreeSize - 1 << std::endl; if ((du - dl)/(lowerBoundQNodeByBandwidth(QIndex, bIndex) + dl) < delta) { for (size_t q = Q->begin(); q < Q->end(); ++q) @@ -151,12 +249,12 @@ void KdeDualTree::MultiBandwidthDualTree() if (meetsDeltaCondition) { /* adjust the current structure, then reinsert it into the queue */ - queueCurrent.dl = deltaLower; - queueCurrent.du = deltaUpper; + queueCurrent.deltaLower = deltaLower; + queueCurrent.deltaUpper = deltaUpper; queueCurrent.bUpperIndex = bUpper; queueCurrent.bLowerIndex = bLower; queueCurrent.priority += PRIORITY_MAX; - nodePriorityQueue.insert(queueCurrent); + nodePriorityQueue.push(queueCurrent); continue; } else @@ -208,28 +306,45 @@ void KdeDualTree::MultiBandwidthDualTree() MultiBandwidthDualTreeBase(Q, T, QIndex, bLower, bUpper); } double priority = pow(Q->bound().MinDistance(T->bound()), 0.5); - if (!Q->is_left() && !T->is_leaf()) + if (!Q->is_leaf() && !T->is_leaf()) { - struct queueNode leftLeft = - {T->left(),Q->left(), 2*QIndex + 1, arma::vec(deltaUpper), - arma::vec(deltaLower), priority, bLower, bUpper}; - struct queueNode leftRight = - {T->left(),Q->right(), 2*QIndex + 2, arma::vec(deltaUpper), - arma::vec(deltaLower), priority, bLower, bUpper}; - struct queueNode rightLeft = - {T->right(),Q->left(), 2*QIndex + 1, arma::vec(deltaUpper), - arma::vec(deltaLower), priority, bLower, bUpper}; - struct queueNode rightRight = - {T->right(),Q->right(), 2*QIndex + 2, arma::vec(deltaUpper), - arma::vec(deltaLower), priority, bLower, bUpper}; - nodePriorityQueue.insert(leftLeft); - nodePriorityQueue.insert(leftRight); - nodePriorityQueue.insert(rightLeft); - nodePriorityQueue.insert(rightRight); + //std::cout << "QIndex for the current non-leaf : " << QIndex << std::endl; + TTree* QLeft = Q->left(); + TTree* QRight = Q->right(); + if (nodeIndices.find(QLeft) == nodeIndices.end()) + { + nodeIndices[QLeft] = nextAvailableNodeIndex; + ++nextAvailableNodeIndex; + } + if (nodeIndices.find(QRight) == nodeIndices.end()) + { + nodeIndices[QRight] = nextAvailableNodeIndex; + ++nextAvailableNodeIndex; + } + size_t QLeftIndex = (*(nodeIndices.find(QLeft))).second; + size_t QRightIndex = (*(nodeIndices.find(QRight))).second; + struct queueNode leftLeft = + {T->left(),Q->left(), QLeftIndex, arma::vec(deltaLower), + arma::vec(deltaUpper), priority, bLower, bUpper}; + struct queueNode leftRight = + {T->left(),Q->right(), QRightIndex, arma::vec(deltaLower), + arma::vec(deltaUpper), priority, bLower, bUpper}; + struct queueNode rightLeft = + {T->right(),Q->left(), QLeftIndex, arma::vec(deltaLower), + arma::vec(deltaUpper), priority, bLower, bUpper}; + struct queueNode rightRight = + {T->right(),Q->right(), QRightIndex, arma::vec(deltaLower), + arma::vec(deltaUpper), priority, bLower, bUpper}; + nodePriorityQueue.push(leftLeft); + nodePriorityQueue.push(leftRight); + nodePriorityQueue.push(rightLeft); + nodePriorityQueue.push(rightRight); } } + return v; } +template void KdeDualTree::Winnow(size_t level, size_t* bLower, size_t* bUpper) @@ -310,14 +425,14 @@ void KdeDualTree::MultiBandwidthDualTreeBase(TTree* Q, { arma::vec diff = queryPoint - referenceData.unsafe_col(t); double distSquared = arma::dot(diff, diff); - size_t bandwidthIndex = upperBIndex; + size_t bandwidthIndex = upperBIndex + 1; while (bandwidthIndex > lowerBIndex) { --bandwidthIndex; - double bandwidth = bandwidths[bandwidthIndex]; - double scaledProduct = pow(distSquared, 0.5) / bandwidth; + double inverseBandwidth = inverseBandwidths[bandwidthIndex]; + double scaledProduct = pow(distSquared, 0.5) * inverseBandwidth; /* TODO: determine the power of the incoming argument */ - double contribution = kernel(scaledProduct); + double contribution = inverseBandwidth * kernel.Evaluate(scaledProduct); if (contribution > DBL_EPSILON) { upperBoundQPointByBandwidth(q, bandwidthIndex) += contribution; @@ -368,6 +483,16 @@ void KdeDualTree::MultiBandwidthDualTreeBase(TTree* Q, sizeOfQNode * log(lowerBoundQNodeByBandwidth(QIndex, bIndex)); } } +template +void KdeDualTree::SetBandwidthBounds(double l, double u) +{ + if (u <= l + DBL_EPSILON || l <= DBL_EPSILON) + { + Log::Fatal << "Incorrect bandwidth range assignment" << std::endl; + } + lowBandwidth = l; + highBandwidth = u; +} }; }; diff --git a/src/contrib/nslagle/myKDE/test_kde_dual_tree.cpp b/src/contrib/nslagle/myKDE/test_kde_dual_tree.cpp new file mode 100644 index 0000000000..88cd61cf45 --- /dev/null +++ b/src/contrib/nslagle/myKDE/test_kde_dual_tree.cpp @@ -0,0 +1,100 @@ +#include +#include "kde_dual_tree.hpp" + +PROGRAM_INFO("Kernel Density Estimation Multibandwidth Dual Tree", + "KDE multibandwidth dual tree calculates density estimates for each " + "query point, given a collection of reference points, using a collection " + "of equidistantly spaced bandwidths\n\n" + "$ kde_dual_tree --reference_file=reference.csv --query_file=query.csv\n" + " --output_file=output.csv --low_bandwidth=0.1 --high_bandwidth=100.0\n" + " --bandwidth_count=10 --epsilon=0.01 --delta=0.01", "kde_dual_tree"); + +PARAM_STRING_REQ("reference_file", "CSV file containing the reference dataset.", + ""); +PARAM_STRING("query_file", "CSV file containing query points", + "", ""); +PARAM_STRING("output_file", "File to output CSV-formatted results into.", "", + "kde_dual_tree_output.csv"); +PARAM_DOUBLE("low_bandwidth", "Low bandwidth", "", 0.1); +PARAM_DOUBLE("high_bandwidth", "Low bandwidth", "", 100.0); +PARAM_INT("bandwidth_count", "Low bandwidth", "", 10); +PARAM_DOUBLE("epsilon", "error tolerance", "", 0.01); +PARAM_DOUBLE("delta", "reversibility tolerance", "", 0.01); + +int main (int argc, char* argv[]) +{ + CLI::ParseCommandLine(argc, argv); + std::string referenceFile = CLI::GetParam("reference_file"); + std::string queryFile = CLI::GetParam("query_file"); + std::string outputFile = CLI::GetParam("output_file"); + arma::mat referenceData; + arma::mat queryData; + double epsilon = CLI::GetParam("epsilon"); + double delta = CLI::GetParam("delta"); + int bandwidthCount = CLI::GetParam("bandwidth_count"); + double lowBandwidth = CLI::GetParam("low_bandwidth"); + double highBandwidth = CLI::GetParam("high_bandwidth"); + + /* check the parameters */ + if (delta < 0.0) + { + Log::Fatal << "Improper delta: " << delta << + "; delta must be positive" << std::endl; + } + if (epsilon < 0.0) + { + Log::Fatal << "Improper epsilon: " << epsilon << + "; epsilon must be positive" << std::endl; + } + if (bandwidthCount <= 0) + { + Log::Fatal << "Improper bandwidth_count: " << bandwidthCount << + "; bandwidth_count must be positive" << std::endl; + } + if (highBandwidth <= lowBandwidth + DBL_EPSILON || lowBandwidth <= 0.0) + { + Log::Fatal << "Improper bandwidth range: " << lowBandwidth << ", " << + highBandwidth << "; bandwidth range must be a positive interval" << std::endl; + } + + if (!data::Load(referenceFile.c_str(), referenceData)) + { + Log::Fatal << "Failed to load the reference file " << referenceFile << std::endl; + } + + Log::Info << "Loaded reference data from " << referenceFile << std::endl; + + std::vector densities; + if (queryFile == "") + { + /* invoke KDE without specific query data */ + KdeDualTree<> kde = KdeDualTree<>(referenceData); + kde.Epsilon() = epsilon; + kde.Delta() = delta; + kde.BandwidthCount() = bandwidthCount; + kde.SetBandwidthBounds(lowBandwidth, highBandwidth); + densities = kde.Calculate(); + } + else + { + /* invoke KDE without specific query data */ + KdeDualTree<> kde = KdeDualTree<>(referenceData, queryData); + kde.Epsilon() = epsilon; + kde.Delta() = delta; + kde.BandwidthCount() = bandwidthCount; + kde.SetBandwidthBounds(lowBandwidth, highBandwidth); + densities = kde.Calculate(); + } + size_t index = 0; + for (std::vector::iterator dIt = densities.begin(); + dIt != densities.end(); + ++dIt) + { + if (*dIt != 0.0) + { + std::cout << "density[" << index << "]=" << *dIt << std::endl; + } + ++index; + } + return 0; +}