dualtree_vkde.h

Go to the documentation of this file.

/* MLPACK 0.2
 *
 * Copyright (c) 2008, 2009 Alexander Gray,
 *                          Garry Boyer,
 *                          Ryan Riegel,
 *                          Nikolaos Vasiloglou,
 *                          Dongryeol Lee,
 *                          Chip Mappus, 
 *                          Nishant Mehta,
 *                          Hua Ouyang,
 *                          Parikshit Ram,
 *                          Long Tran,
 *                          Wee Chin Wong
 *
 * Copyright (c) 2008, 2009 Georgia Institute of Technology
 *
 * This program is free software; you can redistribute it and/or
 * modify it under the terms of the GNU General Public License as
 * published by the Free Software Foundation; either version 2 of the
 * License, or (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
 * 02110-1301, USA.
 */
#ifndef DUALTREE_VKDE_H
#define DUALTREE_VKDE_H

#define INSIDE_DUALTREE_VKDE_H

#include "fastlib/fastlib.h"
#include "gen_metric_tree.h"
#include "dualtree_kde_common.h"
#include "kde_stat.h"
#include "mlpack/allknn/allknn.h"

template<typename TKernel>
class DualtreeVKde {
  
  friend class DualtreeKdeCommon;

 public:
  
  // our tree type using the VKdeStat
  typedef GeneralBinarySpaceTree<DBallBound < LMetric<2>, Vector>, Matrix, VKdeStat<TKernel> > Tree;
    
 private:


  static const int num_initial_samples_per_query_ = 25;

  static const int sample_multiple_ = 10;


  struct datanode *module_;

  bool leave_one_out_;

  double mult_const_;

  ArrayList<TKernel> kernels_;

  Matrix qset_;

  Tree *qroot_;

  Matrix rset_;
  
  Tree *rroot_;

  Vector rset_weights_;

  Vector densities_l_;

  Vector densities_e_;

  Vector densities_u_;

  Vector used_error_;

  Vector n_pruned_;

  double rset_weight_sum_;

  double relative_error_;

  double threshold_;
  
  int num_finite_difference_prunes_;

  int num_monte_carlo_prunes_;

  ArrayList<index_t> old_from_new_queries_;
  
  ArrayList<index_t> old_from_new_references_;


  void DualtreeVKdeBase_(Tree *qnode, Tree *rnode, double probability);

  bool PrunableEnhanced_(Tree *qnode, Tree *rnode, double probability,
                         DRange &dsqd_range, DRange &kernel_value_range, 
                         double &dl, double &du,
                         double &used_error, double &n_pruned,
                         int &order_farfield_to_local,
                         int &order_farfield, int &order_local);

  double EvalUnnormOnSq_(index_t reference_point_index,
                         double squared_distance);

  bool DualtreeVKdeCanonical_(Tree *qnode, Tree *rnode, double probability);

  void PreProcess(Tree *node, bool reference_side);

  void PostProcess(Tree *qnode);
    
 public:


  DualtreeVKde() {
    qroot_ = rroot_ = NULL;
  }

  ~DualtreeVKde() { 
    
    if(qroot_ != rroot_ ) {
      delete qroot_; 
      delete rroot_; 
    } 
    else {
      delete rroot_;
    }

  }


  void get_density_estimates(Vector *results) { 
    results->Init(densities_e_.length());
    
    for(index_t i = 0; i < densities_e_.length(); i++) {
      (*results)[i] = densities_e_[i];
    }
  }


  void Compute(Vector *results) {

    // Set accuracy parameters.
    relative_error_ = fx_param_double(module_, "relative_error", 0.1);
    threshold_ = fx_param_double(module_, "threshold", 0) *
      kernels_[0].CalcNormConstant(qset_.n_rows());
    
    // initialize the lower and upper bound densities
    densities_l_.SetZero();
    densities_e_.SetZero();
    densities_u_.SetAll(rset_weight_sum_);

    // Set zero for error accounting stuff.
    used_error_.SetZero();
    n_pruned_.SetZero();

    // Reset prune statistics.
    num_finite_difference_prunes_ = num_monte_carlo_prunes_ = 0;

    printf("\nStarting variable KDE using %d neighbors...\n",
           (int) fx_param_int_req(module_, "knn"));

    fx_timer_start(NULL, "fast_kde_compute");

    // Preprocessing step for initializing series expansion objects
    PreProcess(rroot_, true);
    if(qroot_ != rroot_) {
      PreProcess(qroot_, false);
    }
    
    // Get the required probability guarantee for each query and call
    // the main routine.
    double probability = fx_param_double(module_, "probability", 1);
    DualtreeVKdeCanonical_(qroot_, rroot_, probability);

    // Postprocessing step for finalizing the sums.
    PostProcess(qroot_);
    fx_timer_stop(NULL, "fast_kde_compute");
    printf("\nFast KDE completed...\n");
    printf("Finite difference prunes: %d\n", num_finite_difference_prunes_);
    printf("Monte Carlo prunes: %d\n", num_monte_carlo_prunes_);

    // Reshuffle the results to account for dataset reshuffling
    // resulted from tree constructions.
    Vector tmp_q_results;
    tmp_q_results.Init(densities_e_.length());
    
    for(index_t i = 0; i < tmp_q_results.length(); i++) {
      tmp_q_results[old_from_new_queries_[i]] =
        densities_e_[i];
    }
    for(index_t i = 0; i < tmp_q_results.length(); i++) {
      densities_e_[i] = tmp_q_results[i];
    }

    // Retrieve density estimates.
    get_density_estimates(results);
  }

  void Init(const Matrix &queries, const Matrix &references,
            const Matrix &rset_weights, bool queries_equal_references, 
            struct datanode *module_in) {

    // point to the incoming module
    module_ = module_in;

    // Set the flag for whether to perform leave-one-out computation.
    leave_one_out_ = fx_param_exists(module_in, "loo") &&
      (queries.ptr() == references.ptr());

    // read in the number of points owned by a leaf
    int leaflen = fx_param_int(module_in, "leaflen", 20);

    // Copy reference dataset and reference weights and compute its
    // sum. rset_weight_sum_ should be the raw sum of the reference
    // weights, ignoring the possibly different normalizing constants
    // in the case of variable-bandwidth case.
    rset_.Copy(references);
    rset_weights_.Init(rset_weights.n_cols());
    rset_weight_sum_ = 0;
    for(index_t i = 0; i < rset_weights.n_cols(); i++) {
      rset_weights_[i] = rset_weights.get(0, i);
      rset_weight_sum_ += rset_weights_[i];
    }

    // Copy the query dataset.
    if(queries_equal_references) {
      qset_.Alias(rset_);
    }
    else {
      qset_.Copy(queries);
    }

    // Construct query and reference trees. Shuffle the reference
    // weights according to the permutation of the reference set in
    // the reference tree.
    fx_timer_start(NULL, "tree_d");
    rroot_ = proximity::MakeGenMetricTree<Tree>(rset_, leaflen,
                                                &old_from_new_references_, 
                                                NULL);
    DualtreeKdeCommon::ShuffleAccordingToPermutation
      (rset_weights_, old_from_new_references_);

    if(queries_equal_references) {
      qroot_ = rroot_;
      old_from_new_queries_.InitCopy(old_from_new_references_);
    }
    else {
      qroot_ = proximity::MakeGenMetricTree<Tree>(qset_, leaflen,
                                                  &old_from_new_queries_, 
                                                  NULL);
    }
    fx_timer_stop(NULL, "tree_d");
    
    // Initialize the density lists
    densities_l_.Init(qset_.n_cols());
    densities_e_.Init(qset_.n_cols());
    densities_u_.Init(qset_.n_cols());

    // Initialize the error accounting stuff.
    used_error_.Init(qset_.n_cols());
    n_pruned_.Init(qset_.n_cols());

    // Initialize the kernels for each reference point.
    int knns = fx_param_int_req(module_, "knn");
    AllkNN all_knn;
    kernels_.Init(rset_.n_cols());
    all_knn.Init(rset_, 20, knns);
    ArrayList<index_t> resulting_neighbors;
    ArrayList<double> squared_distances;    

    fx_timer_start(fx_root, "bandwidth_initialization");
    all_knn.ComputeNeighbors(&resulting_neighbors, &squared_distances);

    for(index_t i = 0; i < squared_distances.size(); i += knns) {
      kernels_[i / knns].Init(sqrt(squared_distances[i + knns - 1]));
    }
    fx_timer_stop(fx_root, "bandwidth_initialization");

    // Renormalize the reference weights according to the bandwidths
    // that have been chosen.
    double min_norm_const = DBL_MAX;
    for(index_t i = 0; i < rset_weights_.length(); i++) {
      double norm_const = kernels_[i].CalcNormConstant(qset_.n_rows());
      min_norm_const = std::min(min_norm_const, norm_const);
    }
    for(index_t i = 0; i < rset_weights_.length(); i++) {
      double norm_const = kernels_[i].CalcNormConstant(qset_.n_rows());
      rset_weights_[i] *= (min_norm_const / norm_const);
    }

    // Compute normalization constant.
    mult_const_ = 1.0 / min_norm_const;
  }

  void PrintDebug() {

    FILE *stream = stdout;
    const char *fname = NULL;

    if((fname = fx_param_str(module_, "fast_kde_output", 
                             "fast_kde_output.txt")) != NULL) {
      stream = fopen(fname, "w+");
    }
    for(index_t q = 0; q < qset_.n_cols(); q++) {
      fprintf(stream, "%g\n", densities_e_[q]);
    }
    
    if(stream != stdout) {
      fclose(stream);
    }
  }

};

#include "dualtree_vkde_impl.h"
#undef INSIDE_DUALTREE_VKDE_H

#endif