dualtree_kde_cv.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_CV_H
#define DUALTREE_KDE_CV_H

#define INSIDE_DUALTREE_KDE_CV_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_cv_common.h"
#include "kde_cv_stat.h"

template<typename TKernelAux>
class DualtreeKdeCV {

  friend class DualtreeKdeCommon;
  
  friend class DualtreeKdeCVCommon;

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


  static const int num_initial_samples_per_query_ = 25;

  static const int sample_multiple_ = 1;


  struct datanode *module_;

  TKernelAux first_ka_;

  TKernelAux second_ka_;

  Matrix rset_;
  
  Tree *rroot_;

  Vector rset_weights_;

  double first_sum_l_;
  
  double first_sum_e_;

  double first_sum_u_;

  double second_sum_l_;

  double second_sum_e_;

  double second_sum_u_;

  double first_mult_const_;
  
  double second_mult_const_;

  double first_used_error_;
  
  double second_used_error_;
  
  double 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_references_;


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

  bool PrunableEnhanced_(Tree *qnode, Tree *rnode, double probability, 
                         DRange &dsqd_range, DRange &first_kernel_value_range,
                         DRange &second_kernel_value_range, double &first_dl,
                         double &first_de, double &first_du, 
                         double &first_used_error, int &first_order,
                         double &second_dl, double &second_de,
                         double &second_du, double &second_used_error, 
                         int &second_order, double &n_pruned);
  
  void EvalUnnormOnSq_(index_t reference_point_index, double squared_distance,
                       double *first_kernel_value,
                       double *second_kernel_value);

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

  void PreProcess(Tree *node);

 public:


  DualtreeKdeCV() {
    rroot_ = NULL;
  }

  ~DualtreeKdeCV() {    
    delete rroot_;
  }




  double Compute() {

    // Compute normalization constant.
    first_mult_const_ = 1.0 / 
      (pow(sqrt(2), rset_.n_rows()) * 
       second_ka_.kernel_.CalcNormConstant(rset_.n_rows()));
    second_mult_const_ = 1.0 /
      second_ka_.kernel_.CalcNormConstant(rset_.n_rows());

    // Set accuracy parameters.
    relative_error_ = fx_param_double(module_, "relative_error", 0.1);
    threshold_ = fx_param_double(module_, "threshold", 0) *
      first_ka_.kernel_.CalcNormConstant(rset_.n_rows());

    // 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(second_ka_.kernel_.bandwidth_sq()));
    fx_timer_start(NULL, "fast_kde_compute");

    // Reset the accumulated sum...
    first_sum_l_ = first_sum_e_ = 0;
    first_sum_u_ = rset_weight_sum_ * rroot_->count();
    second_sum_l_ = second_sum_e_ = 0;
    second_sum_u_ = rset_weight_sum_ * rroot_->count();
    first_used_error_ = second_used_error_ = 0;
    n_pruned_ = 0;

    // Preprocessing step for initializing series expansion objects
    PreProcess(rroot_);
        
    // Get the required probability guarantee for each query and call
    // the main routine.
    double probability = fx_param_double(module_, "probability", 1);
    DualtreeKdeCVCanonical_(rroot_, rroot_, probability);
    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_);

    // Normalize accordingly.
    first_sum_e_ *= (first_mult_const_ / rset_weight_sum_);
    second_sum_e_ *= (second_mult_const_ / rset_weight_sum_);

    // Return the sum of the two sums.
    double lscv_score = 
      (first_sum_e_ - 2.0 * second_sum_e_ +
       2.0 * second_ka_.kernel_.EvalUnnormOnSq(0.0) / 
       second_ka_.kernel_.CalcNormConstant(rset_.n_rows())) /
      ((double) rset_.n_cols());
    return lscv_score;
  }

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

    // point to the incoming module
    module_ = module_in;

    // 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];
    }

    // 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_);
    fx_timer_stop(NULL, "tree_d");

    // Initialize the kernel.
    double bandwidth = fx_param_double_req(module_, "bandwidth");

    // Initialize the series expansion object. I should think about
    // whether this is true for kernels other than Gaussian.
    if(rset_.n_rows() <= 2) {
      first_ka_.Init(sqrt(2) * bandwidth, fx_param_int(module_, "order", 7), 
                     rset_.n_rows());
      second_ka_.Init(bandwidth, fx_param_int(module_, "order", 7), 
                      rset_.n_rows());
    }
    else if(rset_.n_rows() <= 3) {
      first_ka_.Init(sqrt(2) * bandwidth, fx_param_int(module_, "order", 5), 
                     rset_.n_rows());
      second_ka_.Init(bandwidth, fx_param_int(module_, "order", 5), 
                      rset_.n_rows());
    }
    else if(rset_.n_rows() <= 5) {
      first_ka_.Init(sqrt(2) * bandwidth, fx_param_int(module_, "order", 3), 
                     rset_.n_rows());
      second_ka_.Init(bandwidth, fx_param_int(module_, "order", 3), 
                      rset_.n_rows());
    }
    else if(rset_.n_rows() <= 6) {
      first_ka_.Init(sqrt(2) * bandwidth, fx_param_int(module_, "order", 1), 
                     rset_.n_rows());
      second_ka_.Init(bandwidth, fx_param_int(module_, "order", 1), 
                      rset_.n_rows());
    }
    else {
      first_ka_.Init(sqrt(2) * bandwidth, fx_param_int(module_, "order", 0), 
                     rset_.n_rows());
      second_ka_.Init(bandwidth, fx_param_int(module_, "order", 0), 
                      rset_.n_rows()); 
    }
  }
};

#include "dualtree_kde_cv_impl.h"
#undef INSIDE_DUALTREE_KDE_CV_H

#endif