crossvalidation.h

/* 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 DATA_CROSSVALIDATION
#define DATA_CROSSVALIDATION

#include "fastlib/data/dataset.h"
//#include "dataset.h"

#include "fastlib/la/matrix.h"
#include "fastlib/fx/fx.h"

template<class TClassifier>
class SimpleCrossValidator {
  FORBID_ACCIDENTAL_COPIES(SimpleCrossValidator);
  
 public:
  typedef TClassifier Classifier;
  
 private:
  const Dataset *data_;
  datanode *root_module_;
  datanode *kfold_module_;
  int n_folds_;
  int n_classes_;
  const char *classifier_fx_name_;
  index_t n_correct_;
  Matrix confusion_matrix_;
  
 public:
  SimpleCrossValidator() {}
  ~SimpleCrossValidator() {}
  
  void Init(
      const Dataset *data_with_labels,
      int n_labels,
      int default_k,
      struct datanode *module_root,
      const char *classifier_fx_name,
      const char *kfold_fx_name = "kfold");
  
  void Run(bool randomized = false);
  
  
  index_t n_correct() {
    return n_correct_;
  }
  
  index_t n_incorrect() {
    return data_->n_points() - n_correct_;
  }
  
  double portion_correct() {
    return n_correct_ * 1.0 / data_->n_points();
  }
  
  const Matrix& confusion_matrix() const {
    return confusion_matrix_;
  }
  
  const Dataset& data() const {
    return *data_;
  }

 private:
  void SaveTrainTest_(int i_fold,
      const Dataset& train, const Dataset& test) const;
};

template<class TClassifier>
void SimpleCrossValidator<TClassifier>::SaveTrainTest_(
    int i_fold,
    const Dataset& train, const Dataset& test) const {
  String train_name;
  String test_name;
  
  train_name.InitSprintf("train_%d.csv", i_fold);
  test_name.InitSprintf("test_%d.csv", i_fold);
  
  train.WriteCsv(train_name);
  test.WriteCsv(test_name);
}


template<class TClassifier>
void SimpleCrossValidator<TClassifier>::Init(
    const Dataset *data_with_labels,
    int n_labels,
    int default_k,
    struct datanode *module_root,
    const char *classifier_fx_name,
    const char *kfold_fx_name) {
  data_ = data_with_labels;
  
  if (n_labels <= 0) {
    const DatasetFeature *feature =
        &data_->info().feature(data_->n_features() - 1);
    DEBUG_ASSERT_MSG(feature->type() == DatasetFeature::NOMINAL,
        "Must specify number of classes/labels if the feature is not nominal.");
    n_classes_ = feature->n_values();
  } else {
    n_classes_ = n_labels;
  }
  
  root_module_ = module_root;
  kfold_module_ = fx_submodule(module_root, kfold_fx_name);
  classifier_fx_name_ = classifier_fx_name;
  
  n_folds_ = fx_param_int(kfold_module_, "k", default_k);
  
  DEBUG_ONLY(n_correct_ = BIG_BAD_NUMBER);
  
  confusion_matrix_.Init(n_classes_, n_classes_);
  confusion_matrix_.SetZero();
}

template<class TClassifier>
void SimpleCrossValidator<TClassifier>::Run(bool randomized) {
  ArrayList<index_t> permutation;
  
  if (randomized) {
    math::MakeRandomPermutation(data_->n_points(), &permutation);
  } else {
    math::MakeIdentityPermutation(data_->n_points(), &permutation);
  }
  
  n_correct_ = 0;
  
  fx_timer_start(kfold_module_, "total");
  
  for (int i_fold = 0; i_fold < n_folds_; i_fold++) {
    Classifier classifier;
    Dataset test;
    Dataset train;
    index_t local_n_correct = 0;
    datanode *classifier_module = fx_copy_module(root_module_,
        classifier_fx_name_, "%s/%d/%s",
        kfold_module_->key, i_fold, classifier_fx_name_);
    datanode *foldmodule = fx_submodule(classifier_module, "..");

    data_->SplitTrainTest(n_folds_, i_fold, permutation, &train, &test);
    
    if (fx_param_bool(kfold_module_, "save", 0)) {
      SaveTrainTest_(i_fold, train, test);
    }
  
    VERBOSE_MSG(1, "cross: Training fold %d", i_fold);
    fx_timer_start(foldmodule, "train");
    classifier.InitTrain(train, n_classes_, classifier_module);
    fx_timer_stop(foldmodule, "train");
    
    fx_timer_start(foldmodule, "test");
    VERBOSE_MSG(1, "cross: Testing fold %d", i_fold);
    for (index_t i = 0; i < test.n_points(); i++) {
      Vector test_vector_with_label;
      Vector test_vector;
      
      test.matrix().MakeColumnVector(i, &test_vector_with_label);
      test_vector_with_label.MakeSubvector(
          0, test.n_features()-1, &test_vector);
      
      int label_predict = classifier.Classify(test_vector);
      double label_expect_dbl = test_vector_with_label[test.n_features()-1];
      int label_expect = int(label_expect_dbl);
      
      DEBUG_ASSERT(double(label_expect) == label_expect_dbl);
      DEBUG_ASSERT(label_expect < n_classes_);
      DEBUG_ASSERT(label_expect >= 0);
      DEBUG_ASSERT(label_predict < n_classes_);
      DEBUG_ASSERT(label_predict >= 0);
      
      if (label_expect == label_predict) {
        local_n_correct++;
      }
      
      confusion_matrix_.ref(label_expect, label_predict) += 1;
    }
    fx_timer_stop(foldmodule, "test");
    
    fx_format_result(foldmodule, "n_correct", "%"LI"d",
        local_n_correct);
    fx_format_result(foldmodule, "n_incorrect", "%"LI"d",
        test.n_points() - local_n_correct);
    fx_format_result(foldmodule, "p_correct", "%.03f",
        local_n_correct * 1.0 / test.n_points());
    
    n_correct_ += local_n_correct;
  }
  fx_timer_stop(kfold_module_, "total");

  fx_format_result(kfold_module_, "n_points", "%"LI"d",
      data_->n_points());
  fx_format_result(kfold_module_, "n_correct", "%"LI"d",
      n_correct());
  fx_format_result(kfold_module_, "n_incorrect", "%"LI"d",
      n_incorrect());
  fx_format_result(kfold_module_, "p_correct", "%.03f",
      1.0 * portion_correct());
}






template<class TLearner>
class GeneralCrossValidator {
  FORBID_ACCIDENTAL_COPIES(GeneralCrossValidator);
  
 public:
  typedef TLearner Learner;
  
 private:

  int learner_typeid_;
  int n_folds_;
  const Dataset *data_;
  index_t num_data_points_;
  datanode *root_module_;
  datanode *kfold_module_;
  const char *learner_fx_name_;

  
  int clsf_n_classes_;
  index_t clsf_n_correct_;
  Matrix clsf_confusion_matrix_;

  double msq_err_all_folds_;


 public:
  GeneralCrossValidator() {}
  ~GeneralCrossValidator() {}
  void Init(int learner_typeid,
            int default_k,
            const Dataset *data_input,
            struct datanode *module_root,
            const char *learner_fx_name,
            const char *kfold_fx_name = "kfold");

  const Dataset& data() const {
    return *data_;
  }

  void Run(bool randomized);

  index_t clsf_n_correct() {
    return clsf_n_correct_;
  }
  index_t clsf_n_incorrect() {
    return data_->n_points() - clsf_n_correct_;
  }
  double clsf_portion_correct() {
    return clsf_n_correct_ * 1.0 / data_->n_points();
  }
  const Matrix& clsf_confusion_matrix() const {
    return clsf_confusion_matrix_;
  }


 private:
  void SaveTrainValidationSet_(int i_fold,
      const Dataset& train, const Dataset& validation) const;

  void StratifiedSplitCVSet_(int i_fold, index_t num_classes, ArrayList<index_t>& cv_labels_ct, 
                             ArrayList<index_t>& cv_labels_startpos, const ArrayList<index_t>& permutation, Dataset *train, Dataset *validation){
    // Begin stratified splitting for the i-th fold stratified CV
    index_t n_cv_features = data_->n_features();
    
    // detemine the number of data samples for training and validation according to i_fold
    index_t n_cv_validation, i_validation, i_train;
    n_cv_validation = 0;
    for (index_t i_classes=0; i_classes<num_classes; i_classes++) {
      i_validation = 0;
      for (index_t j=0; j<cv_labels_ct[i_classes]; j++) {
        if ((j - i_fold) % n_folds_ == 0) { // point for validation
          i_validation++;
        }
      }
      n_cv_validation = n_cv_validation + i_validation;
    }
    index_t n_cv_train = num_data_points_ - n_cv_validation;
    train->InitBlank();
    train->info().InitContinuous(n_cv_features);
    train->matrix().Init(n_cv_features, n_cv_train);

    validation->InitBlank();
    validation->info().InitContinuous(n_cv_features);
    validation->matrix().Init(n_cv_features, n_cv_validation);

    // make training set and vaidation set by concatenation
    i_train = 0;
    i_validation = 0;
    for (index_t i_classes=0; i_classes<num_classes; i_classes++) {
      for (index_t j=0; j<cv_labels_ct[i_classes]; j++) {
        Vector source, dest;
        if ((j - i_fold) % n_folds_ != 0) { // add to training set
          train->matrix().MakeColumnVector(i_train, &dest);
          i_train++;
        }
        else { // add to validation set
          validation->matrix().MakeColumnVector(i_validation, &dest);
          i_validation++;
        }
        data_->matrix().MakeColumnVector(cv_labels_startpos[i_classes]+j, &source);
        dest.CopyValues(source);
      }
    }
  }

};
  
template<class TLearner>
void GeneralCrossValidator<TLearner>::SaveTrainValidationSet_(
    int i_fold, const Dataset& train, const Dataset& validation) const {
  String train_name;
  String validation_name;
  
  // save training and validation sets for this fold
  train_name.InitSprintf("cv_train_%d.csv", i_fold);
  validation_name.InitSprintf("cv_validation_%d.csv", i_fold);
  
  train.WriteCsv(train_name);
  validation.WriteCsv(validation_name);
}

template<class TLearner>
void GeneralCrossValidator<TLearner>::Init(
    int learner_typeid,
    int default_k,
    const Dataset *data_input,
    struct datanode *module_root,
    const char *learner_fx_name,
    const char *kfold_fx_name) {
  learner_typeid_ = learner_typeid;
  data_ = data_input;

  root_module_ = module_root;
  kfold_module_ = fx_submodule(module_root, kfold_fx_name);
  n_folds_ = fx_param_int(kfold_module_, "k", default_k);
  learner_fx_name_ = learner_fx_name;

  if(learner_typeid_ == 0) {
    // get the number of classes
    clsf_n_classes_ = data_->n_labels();
    clsf_n_correct_ = 0;
    // initialize confusion matrix
    clsf_confusion_matrix_.Init(clsf_n_classes_, clsf_n_classes_);
    clsf_confusion_matrix_.SetZero();
  }
  else if (learner_typeid_ == 1 || learner_typeid_ == 2) {
    clsf_confusion_matrix_.Init(1,1);
    // initialize mean squared error over all folds
    msq_err_all_folds_ = 0.0;
  }

}

template<class TLearner>
void GeneralCrossValidator<TLearner>::Run(bool randomized) {  
  fx_timer_start(kfold_module_, "total");
  num_data_points_ = data_->n_points();

  if (learner_typeid_ == 0) {
    // get label information
    /* list of labels, need to be integers. e.g. [0,1,2] for a 3-class dataset */
    ArrayList<double> cv_labels_list;
    /* array of label indices, after grouping. e.g. [c1[0,5,6,7,10,13,17],c2[1,2,4,8,9],c3[...]]*/
    ArrayList<index_t> cv_labels_index;
    /* counted number of label for each class. e.g. [7,5,8]*/
    ArrayList<index_t> cv_labels_ct;
    /* start positions of each classes in the cv label list. e.g. [0,7,12] */
    ArrayList<index_t> cv_labels_startpos;
    // Get label list and label indices from the cross validation data set
    index_t num_classes = data_->n_labels();

    cv_labels_list.Init();
    cv_labels_index.Init();
    cv_labels_ct.Init();
    cv_labels_startpos.Init();
    data_->GetLabels(cv_labels_list, cv_labels_index, cv_labels_ct, cv_labels_startpos);

    // randomize the original data set within each class if necessary
    ArrayList<index_t> permutation;

    if (randomized) {
      permutation.Init(num_data_points_);
      for (index_t i_classes=0; i_classes<num_classes; i_classes++) {
        ArrayList<index_t> sub_permutation; // within class permut indices
        math::MakeRandomPermutation(cv_labels_ct[i_classes], &sub_permutation);
        // use sub-permutation indicies to form the whole permutation
        if (i_classes==0){
          for (index_t j=0; j<cv_labels_ct[i_classes]; j++)
            permutation[cv_labels_startpos[i_classes]+j] = cv_labels_index[ sub_permutation[j] ];
        }
        else {
          for (index_t j=0; j<cv_labels_ct[i_classes]; j++)
            permutation[cv_labels_startpos[i_classes]+j] = cv_labels_index[ cv_labels_ct[i_classes-1]+sub_permutation[j] ];
        }
        sub_permutation.Clear();
      }
    } // e.g. [10,13,5,17,0,6,7,,4,9,8,1,2,,...]
    else {
      permutation.InitCopy(cv_labels_index, cv_labels_index.size()); // e.g. [0,5,6,7,10,13,17,,1,2,4,8,9,,...]
    }
    // begin CV
    for (int i_fold = 0; i_fold < n_folds_; i_fold++) {
      Learner classifier;
      Dataset train;
      Dataset validation;
      
      index_t local_n_correct = 0;
      datanode *learner_module = fx_copy_module(root_module_,
          learner_fx_name_, "%s/%d/%s",
          kfold_module_->key, i_fold, learner_fx_name_);
      datanode *foldmodule = fx_submodule(learner_module, "..");

      // Split labeled data sets according to i_fold. Use Stratified Cross-Validation to ensure 
      // that approximately the same portion of data (training/validation) are used for each class.
      StratifiedSplitCVSet_(i_fold, num_classes, cv_labels_ct, cv_labels_startpos, permutation, &train, &validation);
      if (fx_param_bool(kfold_module_, "save", 0)) {
        SaveTrainValidationSet_(i_fold, train, validation);
      }
      
      VERBOSE_MSG(1, "cross: Training fold %d", i_fold);
      fx_timer_start(foldmodule, "train");
      // training
      classifier.InitTrain(learner_typeid_, train, learner_module);
      fx_timer_stop(foldmodule, "train");

      // validation; measure method: percent of correctly classified validation samples
      fx_timer_start(foldmodule, "validation");
      VERBOSE_MSG(1, "cross: Validation fold %d", i_fold);

      for (index_t i = 0; i < validation.n_points(); i++) {
        Vector validation_vector_with_label;
        Vector validation_vector;
        
        validation.matrix().MakeColumnVector(i, &validation_vector_with_label);
        validation_vector_with_label.MakeSubvector(0, validation.n_features()-1, &validation_vector);
        // testing (classification)
        int label_predict = int(classifier.Predict(learner_typeid_, validation_vector));
        double label_expect_dbl = validation_vector_with_label[validation.n_features()-1];
        int label_expect = int(label_expect_dbl);

        DEBUG_ASSERT(double(label_expect) == label_expect_dbl);
        DEBUG_ASSERT(label_expect < clsf_n_classes_);
        DEBUG_ASSERT(label_expect >= 0);
        DEBUG_ASSERT(label_predict < clsf_n_classes_);
        DEBUG_ASSERT(label_predict >= 0);
        
        if (label_expect == label_predict) {
          local_n_correct++;
        }
        clsf_confusion_matrix_.ref(label_expect, label_predict) += 1;
      }
      fx_timer_stop(foldmodule, "validation");

      fx_format_result(foldmodule, "local_n_correct", "%"LI"d", local_n_correct);
      fx_format_result(foldmodule, "local_n_incorrect", "%"LI"d", validation.n_points() - local_n_correct);
      fx_format_result(foldmodule, "local_p_correct", "%.03f", local_n_correct * 1.0 / validation.n_points());

      clsf_n_correct_ += local_n_correct;
    }
    fx_timer_stop(kfold_module_, "total");
    
    fx_format_result(kfold_module_, "n_points", "%"LI"d", num_data_points_);
    fx_format_result(kfold_module_, "n_correct", "%"LI"d", clsf_n_correct());
    fx_format_result(kfold_module_, "n_incorrect", "%"LI"d", clsf_n_incorrect());
    fx_format_result(kfold_module_, "p_correct", "%.03f", 1.0 * clsf_portion_correct());
  }
  else if (learner_typeid_ == 1 || learner_typeid_ == 2) {
    double accu_msq_err_all_folds = 0.0;

    // randomize the original data set if necessary
    ArrayList<index_t> permutation;  
    if (randomized) {
      math::MakeRandomPermutation(num_data_points_, &permutation);
    } else {
      math::MakeIdentityPermutation(num_data_points_, &permutation);
    }
    // begin CV
    for (int i_fold = 0; i_fold < n_folds_; i_fold++) {
      Learner learner;
      Dataset train;
      Dataset validation;
      
      double msq_err_local = 0.0;
      double accu_sq_err_local = 0.0;
      datanode *learner_module = fx_copy_module(root_module_,
          learner_fx_name_, "%s/%d/%s",
          kfold_module_->key, i_fold, learner_fx_name_);
      datanode *foldmodule = fx_submodule(learner_module, "..");
      
      // Split general data sets according to i_fold
      data_->SplitTrainTest(n_folds_, i_fold, permutation, &train, &validation);
      
      if (fx_param_bool(kfold_module_, "save", 0)) {
        SaveTrainValidationSet_(i_fold, train, validation);
      }
      
      VERBOSE_MSG(1, "cross: Training fold %d", i_fold);
      fx_timer_start(foldmodule, "train");
      // training
      learner.InitTrain(learner_typeid_, train, learner_module); // 0: dummy number of classes
      fx_timer_stop(foldmodule, "train");
      
      // validation
      fx_timer_start(foldmodule, "validation");
      VERBOSE_MSG(1, "cross: Validation fold %d", i_fold);
      for (index_t i = 0; i < validation.n_points(); i++) {
        Vector validation_vector_with_label;
        Vector validation_vector;
        
        validation.matrix().MakeColumnVector(i, &validation_vector_with_label);
        validation_vector_with_label.MakeSubvector(
                                             0, validation.n_features()-1, &validation_vector);
        // testing
        double value_predict = learner.Predict(learner_typeid_, validation_vector);
        double value_true = validation_vector_with_label[validation.n_features()-1];
        double value_err = value_predict - value_true;
        
        // Calculate squared error: sublevel
        accu_sq_err_local  += pow(value_err, 2);
      }
      fx_timer_stop(foldmodule, "validation");
      
      msq_err_local = accu_sq_err_local / validation.n_points();
      fx_format_result(foldmodule, "local_msq_err", "%f", msq_err_local);

      accu_msq_err_all_folds += msq_err_local;
    }
    fx_timer_stop(kfold_module_, "total");
    
    // Calculate mean squared error: over all folds
    msq_err_all_folds_ = accu_msq_err_all_folds / n_folds_;
    fx_format_result(kfold_module_, "msq_err_all_folds", "%f", msq_err_all_folds_);
  }
  else {
    fprintf(stderr, "Other learner types or Unknown learner type id! Cross validation stops!\n");
    return;
  }
}


#endif