naive_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 NAIVE_KDE_H
#define NAIVE_KDE_H

#include "mlpack/allknn/allknn.h"

template<typename TKernel>
class NaiveKde {
  
  FORBID_ACCIDENTAL_COPIES(NaiveKde);

 private:


  struct datanode *module_;

  Matrix qset_;
  
  Matrix rset_;

  Vector rset_weights_;

  ArrayList<TKernel> kernels_;

  Vector densities_;

  double norm_const_;

 public:


  NaiveKde() {    
  }

  ~NaiveKde() {
  }


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


  void Compute(Vector *results) {

    printf("\nStarting naive KDE...\n");
    fx_timer_start(module_, "naive_kde_compute");

    for(index_t q = 0; q < qset_.n_cols(); q++) {
      
      const double *q_col = qset_.GetColumnPtr(q);

      // Compute unnormalized sum first.
      for(index_t r = 0; r < rset_.n_cols(); r++) {
        const double *r_col = rset_.GetColumnPtr(r);
        double dsqd = la::DistanceSqEuclidean(qset_.n_rows(), q_col, r_col);

        densities_[q] += rset_weights_[r] * kernels_[r].EvalUnnormOnSq(dsqd);
      }

      // Then normalize it.
      densities_[q] /= norm_const_;
    }
    fx_timer_stop(module_, "naive_kde_compute");
    printf("\nNaive KDE completed...\n");

    // retrieve density estimates
    get_density_estimates(results);
  }

  void Compute() {

    printf("\nStarting naive KDE...\n");
    fx_timer_start(module_, "naive_kde_compute");

    for(index_t q = 0; q < qset_.n_cols(); q++) {
      
      const double *q_col = qset_.GetColumnPtr(q);
      
      // Compute unnormalized sum.
      for(index_t r = 0; r < rset_.n_cols(); r++) {
        const double *r_col = rset_.GetColumnPtr(r);
        double dsqd = la::DistanceSqEuclidean(qset_.n_rows(), q_col, r_col);
        
        densities_[q] += rset_weights_[r] * kernels_[r].EvalUnnormOnSq(dsqd);
      }
      // Then, normalize it.
      densities_[q] /= norm_const_;
    }
    fx_timer_stop(module_, "naive_kde_compute");
    printf("\nNaive KDE completed...\n");
  }

  void Init(Matrix &qset, Matrix &rset, struct datanode *module_in) {

    // Use the uniform weights for a moment.
    Matrix uniform_weights;
    uniform_weights.Init(1, rset.n_cols());
    uniform_weights.SetAll(1.0);

    Init(qset, rset, uniform_weights, module_in);
  }

  void Init(Matrix &qset, Matrix &rset, Matrix &reference_weights,
            struct datanode *module_in) {

    // Set the datanode module to be the incoming one.
    module_ = module_in;

    // Get datasets.
    qset_.Copy(qset);
    rset_.Copy(rset);
    rset_weights_.Init(reference_weights.n_cols());
    for(index_t i = 0; i < rset_weights_.length(); i++) {
      rset_weights_[i] = reference_weights.get(0, i);
    }    

    // Compute the normalizing constant.
    double weight_sum = 0;
    for(index_t i = 0; i < rset_weights_.length(); i++) {
      weight_sum += rset_weights_[i];
    }
    
    // Get bandwidth and compute the normalizing constant.
    kernels_.Init(rset_.n_cols());
    if(!strcmp(fx_param_str(module_, "mode", "variablebw"), "variablebw")) {

      // Initialize the kernels for each reference point.
      int knns = fx_param_int_req(module_, "knn");
      AllkNN all_knn;
      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.
      norm_const_ = weight_sum * min_norm_const;
    }
    else {
      for(index_t i = 0; i < kernels_.size(); i++) {
        kernels_[i].Init(fx_param_double_req(module_, "bandwidth"));
      }
      norm_const_ = kernels_[0].CalcNormConstant(qset_.n_rows()) * weight_sum;
    }

    // Allocate density storage.
    densities_.Init(qset.n_cols());
    densities_.SetZero();
  }

  void PrintDebug() {

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

    {
      fname = fx_param_str(module_, "naive_kde_output", 
                           "naive_kde_output.txt");
      stream = fopen(fname, "w+");
    }
    for(index_t q = 0; q < qset_.n_cols(); q++) {
      fprintf(stream, "%g\n", densities_[q]);
    }
    
    if(stream != stdout) {
      fclose(stream);
    }    
  }
  
  void ComputeMaximumRelativeError(const Vector &density_estimates) {
    
    double max_rel_err = 0;
    FILE *stream = fopen("relative_error_output.txt", "w+");
    int within_limit = 0;

    for(index_t q = 0; q < densities_.length(); q++) {
      double rel_err = 
        (fabs(density_estimates[q] - densities_[q]) < DBL_EPSILON) ?
        0 : fabs(density_estimates[q] - densities_[q]) / densities_[q];

      if(isnan(density_estimates[q]) || isinf(density_estimates[q]) || 
         isnan(densities_[q]) || isinf(densities_[q])) {
        printf("Warning: Got infs or nans!\n");
      }

      if(rel_err > max_rel_err) {
        max_rel_err = rel_err;
      }
      if(rel_err <= fx_param_double(module_, "relative_error", 0.01)) {
        within_limit++;
      }

      fprintf(stream, "%g\n", rel_err);
    }
    
    fclose(stream);
    fx_format_result(module_, "maximum_relative_error_for_fast_KDE", "%g", 
                     max_rel_err);
    fx_format_result(module_, "under_relative_error_limit", "%d",
                     within_limit);
  }

};

#endif