dualtree_kde.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_KDE_H
#define DUALTREE_KDE_H

#define INSIDE_DUALTREE_KDE_H

#include "fastlib/fastlib.h"
#include "mlpack/series_expansion/farfield_expansion.h"
#include "mlpack/series_expansion/local_expansion.h"
#include "mlpack/series_expansion/mult_farfield_expansion.h"
#include "mlpack/series_expansion/mult_local_expansion.h"
#include "mlpack/series_expansion/kernel_aux.h"
#include "gen_metric_tree.h"
#include "dualtree_kde_common.h"
#include "kde_stat.h"

const fx_entry_doc kde_main_entries[] = {
  {"data", FX_REQUIRED, FX_STR, NULL,
   "  A file containing reference data.\n"},
  {"query", FX_PARAM, FX_STR, NULL,
   "  A file containing query data (defaults to data).\n"},
  FX_ENTRY_DOC_DONE
};

const fx_entry_doc kde_entries[] = {
  {"bandwidth", FX_PARAM, FX_DOUBLE, NULL,
   "  The bandwidth parameter.\n"},
  {"do_naive", FX_PARAM, FX_BOOL, NULL,
   "  Whether to perform naive computation as well.\n"},
  {"dwgts", FX_PARAM, FX_STR, NULL,
   "  A file that contains the weight of each point. If missing, will\
 assume uniform weight\n"},
  {"fast_kde_output", FX_PARAM, FX_STR, NULL,
   "  A file to receive the results of computation.\n"},
  {"kernel", FX_PARAM, FX_STR, NULL,
   "  The type of kernel to use.\n"},
  {"knn", FX_PARAM, FX_INT, NULL,
   "  The number of k-nearest neighbor to use for variable bandwidth.\n"},
  {"loo", FX_PARAM, FX_BOOL, NULL,
   "  Whether to output the density estimates using leave-one-out.\n"},
  {"mode", FX_PARAM, FX_STR, NULL,
   "  Fixed bandwidth or variable bandwidth mode.\n"},
  {"multiplicative_expansion", FX_PARAM, FX_BOOL, NULL,
   "  Whether to do O(p^D) kernel expansion instead of O(D^p).\n"},
  {"probability", FX_PARAM, FX_DOUBLE, NULL,
   "  The probability guarantee that the relative error accuracy holds.\n"},
  {"relative_error", FX_PARAM, FX_DOUBLE, NULL,
   "  The required relative error accuracy.\n"},
  {"threshold", FX_PARAM, FX_DOUBLE, NULL,
   "  If less than this value, then absolute error bound.\n"},
  {"scaling", FX_PARAM, FX_STR, NULL,
   "  The scaling option.\n"},
  FX_ENTRY_DOC_DONE
};

const fx_module_doc kde_doc = {
  kde_entries, NULL,
  "Performs dual-tree kernel density estimate computation.\n"
};

const fx_submodule_doc kde_main_submodules[] = {
  {"kde", &kde_doc,
   "  Responsible for dual-tree kernel density estimate computation.\n"},
  FX_SUBMODULE_DOC_DONE
};

const fx_module_doc kde_main_doc = {
  kde_main_entries, kde_main_submodules,
  "This is the driver for the kernel density estimator.\n"
};



template<typename TKernelAux>
class DualtreeKde {

  friend class DualtreeKdeCommon;

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


  static const int num_initial_samples_per_query_ = 25;

  static const int sample_multiple_ = 1;


  struct datanode *module_;

  bool leave_one_out_;

  double mult_const_;

  TKernelAux ka_;

  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_farfield_to_local_prunes_;

  int num_farfield_prunes_;
  
  int num_local_prunes_;
  
  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 RefineBoundStatistics_(Tree *destination);

  void DualtreeKdeBase_(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 DualtreeKdeCanonical_(Tree *qnode, Tree *rnode, double probability);

  void PreProcess(Tree *node);

  void PostProcess(Tree *qnode);

 public:


  DualtreeKde() {
    qroot_ = rroot_ = NULL;
  }

  ~DualtreeKde() { 
    
    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) {

    // compute normalization constant
    mult_const_ = 1.0 / ka_.kernel_.CalcNormConstant(qset_.n_rows());

    // Set accuracy parameters.
    relative_error_ = fx_param_double(module_, "relative_error", 0.1);
    threshold_ = fx_param_double(module_, "threshold", 0) *
      ka_.kernel_.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_ =
      num_farfield_to_local_prunes_ = num_farfield_prunes_ = 
      num_local_prunes_ = 0;

    printf("\nStarting fast KDE on bandwidth value of %g...\n",
           sqrt(ka_.kernel_.bandwidth_sq()));
    fx_timer_start(NULL, "fast_kde_compute");

    // Preprocessing step for initializing series expansion objects
    PreProcess(rroot_);
    if(qroot_ != rroot_) {
      PreProcess(qroot_);
    }
    
    // Get the required probability guarantee for each query and call
    // the main routine.
    double probability = fx_param_double(module_, "probability", 1);
    DualtreeKdeCanonical_(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_);
    printf("F2L prunes: %d\n", num_farfield_to_local_prunes_);
    printf("F prunes: %d\n", num_farfield_prunes_);
    printf("L prunes: %d\n", num_local_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_.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 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 kernel.
    double bandwidth = fx_param_double_req(module_, "bandwidth");

    // initialize the series expansion object
    if(qset_.n_rows() <= 2) {
      ka_.Init(bandwidth, fx_param_int(module_, "order", 7), qset_.n_rows());
    }
    else if(qset_.n_rows() <= 3) {
      ka_.Init(bandwidth, fx_param_int(module_, "order", 5), qset_.n_rows());
    }
    else if(qset_.n_rows() <= 5) {
      ka_.Init(bandwidth, fx_param_int(module_, "order", 3), qset_.n_rows());
    }
    else if(qset_.n_rows() <= 6) {
      ka_.Init(bandwidth, fx_param_int(module_, "order", 1), qset_.n_rows());
    }
    else {
      ka_.Init(bandwidth, fx_param_int(module_, "order", 0), qset_.n_rows());
    }
  }

  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_kde_impl.h"
#undef INSIDE_DUALTREE_KDE_H

#endif