original_ifgt.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 ORIGINAL_IFGT_H
#define ORIGINAL_IFGT_H

#include "fastlib/fastlib.h"

class OriginalIFGT {

 private:


  struct datanode *module_;

  int dim_;

  int num_reference_points_;

  Matrix query_set_;

  Matrix reference_set_;

  Vector reference_weights_;

  double bandwidth_;

  double bandwidth_sq_;

  double bandwidth_factor_;

  double epsilon_;

  int pterms_;

  int total_num_coeffs_;
  
  Matrix weighted_coeffs_;

  Matrix unweighted_coeffs_;

  int num_cluster_desired_;

  double cut_off_radius_;

  double max_radius_cluster_;

  Matrix cluster_centers_;

  ArrayList<int> index_during_clustering_;

  ArrayList<int> cluster_index_;

  Vector cluster_radii_;

  ArrayList<int> num_reference_points_in_cluster_;

  Vector densities_;

  void TaylorExpansion() {
      
    Vector tmp_coeffs;
    tmp_coeffs.Init(total_num_coeffs_);
    
    ComputeUnweightedCoeffs_(tmp_coeffs);
    
    ComputeWeightedCoeffs_(tmp_coeffs);
    
    return;     
  }

  void ComputeUnweightedCoeffs_(Vector &taylor_coeffs) {
        
    ArrayList<int> heads;
    heads.Init(dim_ + 1);
    ArrayList<int> cinds;
    cinds.Init(total_num_coeffs_);
    
    for (int i = 0; i < dim_; i++) {
      heads[i] = 0;
    }
    heads[dim_] = INT_MAX;
    
    cinds[0] = 0;
    taylor_coeffs[0] = 1.0;
    
    for(int k = 1, t = 1, tail = 1; k < pterms_; k++, tail = t) {
      for(int i = 0; i < dim_; i++) {
        int head = heads[i];
        heads[i] = t;
        for(int j = head; j < tail; j++, t++) {
          cinds[t] = (j < heads[i+1]) ? cinds[j] + 1 : 1;
          taylor_coeffs[t] = 2.0 * taylor_coeffs[j];
          taylor_coeffs[t] /= (double) cinds[t];
        }
      }
    }
    return;    
  }

  void ComputeWeightedCoeffs_(Vector &taylor_coeffs) {
    
    Vector dx;
    dx.Init(dim_);
    Vector prods;
    prods.Init(total_num_coeffs_);
    ArrayList<int> heads;
    heads.Init(dim_);
    
    // initialize coefficients for all clusters to be zero.
    weighted_coeffs_.SetZero();
    unweighted_coeffs_.SetZero();
    
    for(int n = 0; n < num_reference_points_; n++) {
      
      int ix2c = cluster_index_[n];
      double sum = 0.0;
      
      for(int i = 0; i < dim_; i++) {
        dx[i] = (reference_set_.get(i, n) - cluster_centers_.get(i, ix2c)) / 
          bandwidth_factor_;

        sum -= dx[i] * dx[i];
        heads[i] = 0;
      }
      
      prods[0] = exp(sum);
      for(int k = 1, t = 1, tail = 1; k < pterms_; k++, tail = t) {
        
        for (int i = 0; i < dim_; i++) {
          int head = heads[i];
          heads[i] = t;
          for(int j = head; j < tail; j++, t++)
            prods[t] = dx[i] * prods[j];
        } // for i
      } // for k
      
      // compute the weighted coefficients and unweighted coefficients.
      for(int i = 0; i < total_num_coeffs_; i++) {
        weighted_coeffs_.set(i, ix2c, weighted_coeffs_.get(i, ix2c) +
                             reference_weights_[n] * prods[i]);
        unweighted_coeffs_.set(i, ix2c, unweighted_coeffs_.get(i, ix2c) +
                               prods[i]);
      }
      
    }// for n
    
    // normalize by the Taylor coefficients.
    for(int k = 0; k < num_cluster_desired_; k++) {
      for(int i = 0; i < total_num_coeffs_; i++) {
        weighted_coeffs_.set(i, k, weighted_coeffs_.get(i, k) * 
                             taylor_coeffs[i]);
        unweighted_coeffs_.set(i, k, unweighted_coeffs_.get(i, k) * 
                               taylor_coeffs[i]);
      }
    }
    return;
  }

  double ComputeCenters() {

    // set cluster max radius to zero
    max_radius_cluster_ = 0;

    // clear all centers.
    cluster_centers_.SetZero();
    
    // Compute the weighted centroid for each cluster.
    for(int j = 0; j < dim_; j++) {
      for(int i = 0; i < num_reference_points_; i++) {
        
        cluster_centers_.set(j, cluster_index_[i],
                             cluster_centers_.get(j, cluster_index_[i]) +
                             reference_set_.get(j, i));
      }
    }
    
    for(int j = 0; j < dim_; j++) {
      for(int i = 0; i < num_cluster_desired_; i++) {
        cluster_centers_.set(j, i, cluster_centers_.get(j, i) /
                             num_reference_points_in_cluster_[i]);
      }
    }
    
    // Now loop through and compute the radius of each cluster.
    cluster_radii_.SetZero();
    for(int i = 0; i < num_reference_points_; i++) {
      Vector reference_pt;
      reference_set_.MakeColumnVector(i, &reference_pt);

      // the index of the cluster this reference point belongs to.
      int cluster_id = cluster_index_[i];
      Vector center;
      cluster_centers_.MakeColumnVector(cluster_id, &center);
      cluster_radii_[cluster_id] = 
        std::max(cluster_radii_[cluster_id], 
                 sqrt(la::DistanceSqEuclidean(reference_pt, center)));
      max_radius_cluster_ =
        std::max(max_radius_cluster_, cluster_radii_[cluster_id]);
    }

    return max_radius_cluster_;
  }

  double KCenterClustering() {
    
    Vector distances_to_center;
    distances_to_center.Init(num_reference_points_);
    
    // randomly pick one node as the first center.
    srand( (unsigned)time( NULL ) );
    int ind = rand() % num_reference_points_;
    
    // add the ind-th node to the first center.
    index_during_clustering_[0] = ind;
    Vector first_center;
    reference_set_.MakeColumnVector(ind, &first_center);
    
    // compute the distances from each node to the first center and
    // initialize the index of the cluster ID to zero for all
    // reference points.
    for(int j = 0; j < num_reference_points_; j++) {
      Vector reference_point;
      reference_set_.MakeColumnVector(j, &reference_point);
      
      distances_to_center[j] = (j == ind) ? 
        0.0:la::DistanceSqEuclidean(reference_point, first_center);
      cluster_index_[j] = 0;
    }
    
    // repeat until the desired number of clusters is reached.
    for(int i = 1; i < num_cluster_desired_; i++) {
      
      // Find the reference point that is farthest away from the
      // current center.
      ind = IndexOfLargestElement(distances_to_center);
      
      // Add the ind-th node to the centroid list.
      index_during_clustering_[i] = ind;
      
      // Update the distances from each point to the current center.
      Vector center;
      reference_set_.MakeColumnVector(ind, &center);
      
      for (int j = 0; j < num_reference_points_; j++) {
        Vector reference_point;
        reference_set_.MakeColumnVector(j, &reference_point);
        double d = (j == ind)? 
          0.0:la::DistanceSqEuclidean(reference_point, center);
        
        if (d < distances_to_center[j]) {
          distances_to_center[j] = d;
          cluster_index_[j] = i;
        }
      }
    }
    
    // Find the maximum radius of the k-center algorithm.
    ind = IndexOfLargestElement(distances_to_center);
    
    double radius = distances_to_center[ind];
    
    
    for(int i = 0; i < num_cluster_desired_; i++) {
      num_reference_points_in_cluster_[i] = 0;
    }
    // tally up the number of reference points for each cluster.
    for (int i = 0; i < num_reference_points_; i++) {
      num_reference_points_in_cluster_[cluster_index_[i]]++;
    }
    
    return sqrt(radius);  
  }

  int IndexOfLargestElement(const Vector &x) {
    
    int largest_index = 0;
    double largest_quantity = -DBL_MAX;
    
    for(int i = 0; i < x.length(); i++) {
      if(largest_quantity < x[i]) {
        largest_quantity = x[i];
        largest_index = i;
      }
    }
    return largest_index;
  }

  void NormalizeDensities() {
    double norm_const = pow(2 * math::PI * bandwidth_sq_, 
                            query_set_.n_rows() / 2.0) *
      reference_set_.n_cols();

    for(index_t q = 0; q < query_set_.n_cols(); q++) {
      densities_[q] /= norm_const;
    }
  }

  void IFGTChooseTruncationNumber_() {  

    double rx = max_radius_cluster_;
    double max_diameter_of_the_datasets = sqrt(dim_);
    
    double two_h_square = bandwidth_factor_ * bandwidth_factor_;
    
    double r = min(max_diameter_of_the_datasets, 
                   bandwidth_factor_ * sqrt(log(1 / epsilon_)));
    
    int p_ul=300;
    
    double rx_square = rx * rx;
    
    double error = 1;
    double temp = 1;
    int p = 0;
    while((error > epsilon_) & (p <= p_ul)) {
      p++;
      double b = min(((rx + sqrt((rx_square) + (2 * p * two_h_square))) / 2),
                     rx + r);
      double c = rx - b;
      temp = temp * (((2 * rx * b) / two_h_square) / p);
      error = temp * (exp(-(c * c) / two_h_square));                    
    }   
    
    // update the truncation order.
    pterms_ = p;

    // update the cut-off radius
    cut_off_radius_ = r;
    
  }
  
  void IFGTChooseParameters_(int max_num_clusters) {
    
    // for references and queries that fit in the unit hypercube, this
    // assumption is true, but for general case it is not.
    double max_diamater_of_the_datasets = sqrt(dim_);
    
    double two_h_square = bandwidth_factor_ * bandwidth_factor_;

    // The cut-off radius.
    double r = min(max_diamater_of_the_datasets, 
                   bandwidth_factor_ * sqrt(log(1 / epsilon_)));
    
    // Upper limit on the truncation number.
    int p_ul=200; 
    
    num_cluster_desired_ = 1;
    
    double complexity_min=1e16;
    double rx;

    for(int i = 0; i < max_num_clusters; i++){
     
      // Compute an estimate of the maximum cluster radius.
      rx = pow((double) i + 1, -1.0 / (double) dim_);
      double rx_square = rx * rx;

      // An estimate of the number of neighbors.
      double n = std::min(i + 1.0, pow(r / rx, (double) dim_));
      double error = 1;
      double temp = 1;
      int p = 0;

      // Choose the truncation order.
      while((error > epsilon_) & (p <= p_ul)) {
        p++;
        double b = 
          std::min(((rx + sqrt((rx_square) + (2 * p * two_h_square))) / 2.0),
                   rx + r);
        double c = rx - b;
        temp = temp * (((2 * rx * b) / two_h_square) / p);
        error = temp * (exp(-(c * c) / two_h_square));
      }
      double complexity = (i + 1) + log((double) i + 1) + 
        ((1 + n) * math::BinomialCoefficient(p - 1 + dim_, dim_));
        
      if(complexity < complexity_min) {
        complexity_min = complexity;
        num_cluster_desired_ = i + 1;
        pterms_ = p;
      }
    }    
  }

 public:


  OriginalIFGT() {}

  ~OriginalIFGT() {}

  

  void get_density_estimates(Vector *results) {
    results->Init(densities_.length());

    for(index_t i = 0; i < densities_.length(); i++) {
      (*results)[i] = densities_[i];
    }
  }


  void Init(Matrix &queries, Matrix &references, struct datanode *module_in) {
    
    // set module to the incoming one.
    module_ = module_in;
    
    // set dimensionality
    dim_ = references.n_rows();
        
    // set up query set and reference set.
    query_set_.Copy(queries);
    reference_set_.Copy(references);
    num_reference_points_ = reference_set_.n_cols();
    
    // By default, the we do uniform weights only
    reference_weights_.Init(reference_set_.n_cols());
    reference_weights_.SetAll(1);

    // initialize density estimate vector
    densities_.Init(query_set_.n_cols());
    
    // A "hack" such that the code uses the proper Gaussian kernel.
    bandwidth_ = fx_param_double_req(module_, "bandwidth");
    bandwidth_sq_ = bandwidth_ * bandwidth_;
    bandwidth_factor_ = sqrt(2) * bandwidth_;
    
    // Read in the desired absolute error accuracy.
    epsilon_ = fx_param_double(module_, "absolute_error", 0.1);

    // This is the upper limit on the number of clusters.
    int cluster_limit = (int) ceilf(20.0 * sqrt(dim_) / sqrt(bandwidth_));
    
    printf("Automatic parameter selection phase...\n");

    printf("Preprocessing phase for the original IFGT...\n");

    fx_timer_start(module_, "ifgt_kde_preprocess");
    IFGTChooseParameters_(cluster_limit);
    printf("Chose %d clusters...\n", num_cluster_desired_);
    printf("Tentatively chose %d truncation order...\n", pterms_);

    // Allocate spaces for storing coefficients and clustering information.
    cluster_centers_.Init(dim_, num_cluster_desired_);
    index_during_clustering_.Init(num_cluster_desired_);
    cluster_index_.Init(num_reference_points_);
    cluster_radii_.Init(num_cluster_desired_);
    num_reference_points_in_cluster_.Init(num_cluster_desired_);    
    
    printf("Now clustering...\n");

    // Divide the source space into num_cluster_desired_ parts using
    // K-center algorithm
    max_radius_cluster_ = KCenterClustering();
    
    // computer the center of the sources
    ComputeCenters();

    // Readjust the truncation order based on the actual clustering result.
    IFGTChooseTruncationNumber_();
    // pd = C_dim^(dim+pterms-1)
    total_num_coeffs_ = 
      (int) math::BinomialCoefficient(pterms_ + dim_ - 1, dim_);
    weighted_coeffs_.Init(total_num_coeffs_, num_cluster_desired_);
    unweighted_coeffs_.Init(total_num_coeffs_, num_cluster_desired_);

    printf("Maximum radius generated in the cluster: %g...\n",
                max_radius_cluster_);
    printf("Truncation order updated to %d after clustering...\n", 
                pterms_);

    // Compute coefficients.    
    printf("Now computing Taylor coefficients...\n");
    TaylorExpansion();
    printf("Taylor coefficient computation finished...\n");
    fx_timer_stop(module_, "ifgt_kde_preprocess");
    printf("Preprocessing step finished...\n");
  }

  void Compute() {
    
    printf("Starting the original IFGT-based KDE computation...\n");

    fx_timer_start(module_, "original_ifgt_kde_compute");

    Vector dy;
    dy.Init(dim_);
    
    Vector tempy;
    tempy.Init(dim_);
    
    Vector prods;
    prods.Init(total_num_coeffs_);
    
    ArrayList<int> heads;
    heads.Init(dim_);
    
    // make sure the sum for each query point starts at zero.
    densities_.SetZero();
    
    for(int m = 0; m < query_set_.n_cols(); m++) {      
      
      // loop over each cluster and evaluate Taylor expansions.
      for(int kn = 0; kn < num_cluster_desired_; kn++) {
        
        double sum2 = 0.0;
        
        // compute the ratio of the squared distance between each query
        // point and each cluster center to the bandwidth factor.
        for (int i = 0; i < dim_; i++) {
          dy[i] = (query_set_.get(i, m) - cluster_centers_.get(i, kn)) / 
            bandwidth_factor_;
          sum2 += dy[i] * dy[i];
        }
        
        // If the ratio is greater than the cut-off, this cluster's
        // contribution is ignored.
        if (sum2 > (cut_off_radius_ + cluster_radii_[kn]) /
            (bandwidth_factor_ * bandwidth_factor_)) {
          continue;
        }
        
        for(int i = 0; i < dim_; i++) {
          heads[i] = 0;
        }
        
        prods[0] = exp(-sum2);          
        for(int k = 1, t = 1, tail = 1; k < pterms_; k++, tail = t) {
          for (int i = 0; i < dim_; i++) {
            int head = heads[i];
            heads[i] = t;
            for(int j = head; j < tail; j++, t++)
              prods[t] = dy[i] * prods[j];
          } // for i
        }// for k
        
        for(int i = 0; i < total_num_coeffs_; i++) {
          densities_[m] += weighted_coeffs_.get(i, kn) * prods[i];
        }
        
      } // for each cluster
    } //for each query point

    // normalize density estimates
    NormalizeDensities();

    fx_timer_stop(module_, "original_ifgt_kde_compute");
    printf("Computation finished...\n");
    return;
  }

  void PrintDebug() {
    
    FILE *stream = stdout;
    const char *fname = NULL;
    
    if((fname = fx_param_str(module_, "ifgt_kde_output", NULL)) != NULL) {
      stream = fopen(fname, "w+");
    }
    for(index_t q = 0; q < query_set_.n_cols(); q++) {
      fprintf(stream, "%g\n", densities_[q]);
    }
    
    if(stream != stdout) {
      fclose(stream);
    }
  }

};

#endif