fn_fqmi.hpp

Go to the documentation of this file.

// Copyright (C) 2011 the authors listed below
// http://ecocpak.sourceforge.net
// 
// Authors:
// - Dimitrios Bouzas (bouzas at ieee dot org)
// - Nikolaos Arvanitopoulos (niarvani at ieee dot org)
// - Anastasios Tefas (tefas at aiia dot csd dot auth dot gr)
// 
// This file is part of the ECOC PAK C++ library. It is 
// provided without any warranty of fitness for any purpose.
//
// You can redistribute this file and/or modify it under 
// the terms of the GNU Lesser General Public License (LGPL) 
// as published by the Free Software Foundation, either 
// version 3 of the License or (at your option) any later 
// version.
// (see http://www.opensource.org/licenses for more info)





double 
compute_sigma
  (
  const vector<ClassData>& classes_vector
  )
  {  
  double maximum = 0.0;
  double d = 0.0;
  
  // create samples matrix
  mat samples = create_samples_matrix(classes_vector);
  
  const u32 n_rows_j = samples.n_rows;
  const u32 n_rows_i = samples.n_rows - 1;

        // find the minimum distance between all the samples in the set
  for(u32 i  = 0; i < n_rows_i; i++)
    {
    for(u32 j = i + 1; j < n_rows_j; j++)
      {
      // compute distance as the Euclidean norm
      d = norm(samples.row(i) - samples.row(j), 2);
      
      if(maximum < d)
        {
        maximum = d;
        }

      }
    
    } 

  return maximum / 2.0;
  }



template<typename T1, typename T2>
typename T1::elem_type
fqmi
  (
  const Base<typename T1::elem_type, T1>& X, 
  const Base<typename T2::elem_type, T2>& Y,
  const      typename T1::elem_type       sigma
  )
  {
  arma_extra_debug_sigprint();
  
  // alias to element type of samples' matrix
  typedef typename T1::elem_type eT1;
  
  // alias to element type of labels' vector 
  typedef typename T2::elem_type eT2;
  
  // unwrap samples matrix
  const unwrap<T1>             tmp_X(X.get_ref());
  const   Mat<eT1>& samples =  tmp_X.M;
  
  // unwrap labels vector
  const unwrap<T2>        tmp_Y(Y.get_ref());
  const Mat<eT2> labels = tmp_Y.M;
  
  // number of samples
  const u32 n_samples = samples.n_rows;
  
  // number of attributes
  const u32 n_attributes = samples.n_cols;
  
  // check whether input labels is a vector
  arma_debug_check
    (
    (labels.is_vec() == false), 
    "fqmi(): Input class labels should be a vector"
    );
  
  // check whether input labels number of elements is equal to the 
  // number of samples
  arma_debug_check
    (
    labels.n_elem != n_samples, 
    "fqmi(): Input class labels must have the same number of elements with the number of input samples"
    );
  
  // Initialize VIN, VALL and VBTW terms to zero.
  
  // VIN can be seen as interactions between pairs of
  // samples inside each class.
  eT1 vin  = 0;
  
  // VBTW consists of interactions between samples of
  // each class against all other samples.
  eT1 vbtw = 0;
  
  // VALL consists of interactions between all pairs of
  // samples, regardless of class membership.
  eT1 vall = 0;

  // find number of classes
  // find number of samples in each class
  // create new labels that start from one if necessary
  u32 n_classes = 0;   
  Col<eT1> n_samples_each_class;
  Col<eT2> new_labels = 
    process_labels(labels, n_samples_each_class, n_classes);

  // compute normalization constant for later use
  eT1 n_samples_square = 1.0 / eT1(n_samples * n_samples);

  // compute normalization constant for later use
  eT1 v = 2 * sigma * sigma;
  
  // compact expression used in the kernel for later use
  eT1 exp_sigma = eT1(-1) / eT1(2 * v);
  
  // compute normalization constant for later use
  eT1 nc = 
    eT1(1) /  eT1
                (
                std::sqrt
                  (
                  std::pow(2.0 * math::pi(), 
                  n_attributes) * std::pow(v, n_attributes)
                  )
                );
  
  // initialize temporary vector VBTW_tmp to hold intermediate 
  // computations of VBTW term actually each position of vbtw_tmp 
  // corresponds to a class label and holds the sums of the VBTW term 
  // for each of the available classes
  Col<eT1> vbtw_tmp = zeros< Col<eT1> >(n_classes);
  
  // initialize temporary scalar to hold intermediate computation of the 
  // kernel
  eT1 tmp = 0;
  
  // double for to compute the interactions of all samples with each 
  // other it needs N^2 / 2 iterations, where N is the number of samples
  
  // iterate through samples
  for(u32 i = 0; i < n_samples; i++)
    {
    // compute self interaction of the ith sample
    // exponential becomes zero and consequently 
    // interaction is a constant as such update
    // VIN, VALL and VBTW accordingly
    vall += nc;
    vbtw_tmp[labels[i] - 1] += nc;
    vin  += nc;
    
    // compute interaction of ith sample with the rest of the samples
    for(u32 j = i + 1; j < n_samples; j++)
      {
      // compute the difference between ith and jth vectors
      Row<eT1> centered_sample = samples.row(i) - samples.row(j);
  
      // compute kernel
      tmp = 
        nc * std::exp
               (
               exp_sigma * accu(centered_sample % centered_sample)
               );

      // update each of term VALL, VBTW and VIN accordingly
      
      // multiply the kernel result by 2 because k(i,j) = k(j,i)
      vall += (2 * tmp);
      
      // update the sum jth sample's class
      vbtw_tmp[labels[j] - 1] += tmp;
      
      // update the sum jth sample's class
      vbtw_tmp[labels[i] - 1] += tmp;
      
      // kernel result contributes to the VIN sum if ith and jth
      // samples belong to the same class (i.e., their respective
      // class labels are equal)
      if(labels[i] == labels[j])
        {
        // multiply the kernel result by 2 because k(i,j) = k(j,i)
        vin += (2 * tmp);
        }
   
      }
    
    }

  // normalization of VIN, VBTW and VALL terms
  vall = 
    accu
      (
        (
        n_samples_each_class % n_samples_each_class
        ) * 
        n_samples_square * vall
      ) * 
      n_samples_square;
      
  vbtw = 
    n_samples_square * 
    accu(vbtw_tmp % n_samples_each_class) / eT1(n_samples);
    
  vin *= n_samples_square;
  
  // compute and return MI value from vin, vall and vbtw terms
  
  // equation is obtained by using a Parzen estimate with a
  // symmetrical Gaussian kernel of width sigma
  return (vin + vall) - (2.0 * vbtw);
  }