added dataset classification program, added easy option if hyperparam are set manually

FMKGPHyperparameterOptimization.cpp

@@ -261,7 +261,6 @@ void FMKGPHyperparameterOptimization::updateEigenVectors()
   if ( verbose )
   {
     std::cerr << "FMKGPHyperparameterOptimization::updateEigenVectors -- size of ikmsums: " << ikmsums.size() << std::endl;
-    std::cerr << "class of first object: " << ikmsums.begin()->first << std::endl;
   }
   
   if ( learnBalanced )
@@ -280,15 +279,12 @@ void FMKGPHyperparameterOptimization::updateEigenVectors()
   }
   else
   {
-    std::cerr << "not balanced, considere for VarApprox: " << nrOfEigenvaluesToConsiderForVarApprox << " eigenvalues" << std::endl;
-    std::cerr << "and for simple: " << nrOfEigenvaluesToConsider << std::endl;
-    if (nrOfEigenvaluesToConsiderForVarApprox > 1)
-      nrOfEigenvaluesToConsiderForVarApprox = 1;
     //compute the largest eigenvalue of K + noise
     eigenMax.resize(1);
     eigenMaxVectors.resize(1);    
     
-    eig->getEigenvalues ( * ( ikmsums.begin()->second ),  eigenMax[0], eigenMaxVectors[0], nrOfEigenvaluesToConsiderForVarApprox );
+    //TODO check why we are only interested in the largest EW!
+    eig->getEigenvalues ( * ( ikmsums.begin()->second ),  eigenMax[0], eigenMaxVectors[0], 1 /* we are only interested in the largest eigenvalue here*/ );
   }
 }
 
@@ -430,7 +426,8 @@ void FMKGPHyperparameterOptimization::performOptimization ( std::map<int, GPLike
       OPTIMIZATION::matrix_type hyperp ( 1, 1, value );
       gplikes.begin()->second->setParameterLowerBound ( value );
       gplikes.begin()->second->setParameterUpperBound ( value );
-      gplikes.begin()->second->evaluate ( hyperp );
+      //we do not need to compute the likelihood here - we are only interested in directly obtaining alpha vectors
+      gplikes.begin()->second->computeAlphaDirect( hyperp );
     }
   }
 
@@ -538,6 +535,7 @@ void FMKGPHyperparameterOptimization::computeMatricesAndLUTs ( const std::map<in
       fmk->hik_prepare_alpha_multiplications ( i->second, A, B );
       A.setIoUntilEndOfFile ( false );
       B.setIoUntilEndOfFile ( false );
+
       precomputedA[ i->first ] = A;
       precomputedB[ i->first ] = B;
 

GPHIKClassifier.cpp

@@ -196,6 +196,7 @@ void GPHIKClassifier::train ( const std::vector< NICE::SparseVector *> & example
   Timer t;
   t.start();
   FastMinKernel *fmk = new FastMinKernel ( examples, noise, this->debug );
+  
   t.stop();
   if (verbose)
     std::cerr << "Time used for setting up the fmk object: " << t.getLast() << std::endl;  

GPLikelihoodApprox.cpp

@@ -125,6 +125,83 @@ void GPLikelihoodApprox::calculateLikelihood ( double mypara, const FeatureMatri
   cerr << "OPTGT: " << mypara << " " << gt_nlikelihood << " " << gt_logdet << " " << gt_dataterm << endl;
 }
 
+void GPLikelihoodApprox::computeAlphaDirect(const OPTIMIZATION::matrix_type & x)
+{
+  Timer t;
+//   NICE::Vector diagonalElements;
+  
+//   ikm->getDiagonalElements ( diagonalElements );
+
+  // set simple jacobi pre-conditioning
+  ILSConjugateGradients *linsolver_cg = dynamic_cast<ILSConjugateGradients *> ( linsolver );
+
+//   //TODO why do we need this?  
+//   if ( linsolver_cg != NULL )
+//     linsolver_cg->setJacobiPreconditioner ( diagonalElements );
+  
+
+  // all alpha vectors will be stored!
+  std::map<int, NICE::Vector> alphas;
+
+  // This has to be done m times for the multi-class case
+  if (verbose)
+    std::cerr << "run ILS for every bin label. binaryLabels.size(): " << binaryLabels.size() << std::endl;
+  for ( std::map<int, NICE::Vector>::const_iterator j = binaryLabels.begin(); j != binaryLabels.end() ; j++)
+  {
+    // (b) y^T (K+sI)^{-1} y
+    int classCnt = j->first;
+    if (verbose)
+    {
+      std::cerr << "Solving linear equation system for class " << classCnt << " ..." << std::endl;
+      std::cerr << "Size of the kernel matrix " << ikm->rows() << std::endl;
+      std::cerr << "binary label: " << j->second << std::endl;
+    }
+
+    /** About finding a good initial solution
+     * K~ = K + sigma^2 I
+     *
+     * (0) we have already estimated alpha for a previous hyperparameter, then
+     *     we should use this as an initial estimate. According to my quick
+     *     tests this really helps!
+     * (1) K~ \approx lambda_max v v^T
+     * \lambda_max v v^T * alpha = y     | multiply with v^T from left
+     * => \lambda_max v^T alpha = v^T y
+     * => alpha = y / lambda_max could be a good initial start
+     * If we put everything in the first equation this gives us
+     * v = y ....which is somehow a weird assumption (cf Kernel PCA)
+     *  This reduces the number of iterations by 5 or 8
+     */
+    Vector alpha;
+    
+    if ( (usePreviousAlphas) && (lastAlphas != NULL) )
+    {
+      std::map<int, NICE::Vector>::iterator alphaIt = lastAlphas->begin();
+      alpha = (*lastAlphas)[classCnt];
+    }
+    else  
+    {
+      //TODO hand over the eigenmax
+      alpha = (binaryLabels[classCnt] ); //* (1.0 / eigenmax[0]) );
+    }
+    
+    NICE::Vector initialAlpha;
+    if ( verbose )
+     initialAlpha = alpha;
+
+    if ( verbose )
+      std::cerr << "Using the standard solver ..." << std::endl;
+
+    t.start();
+    linsolver->solveLin ( *ikm, binaryLabels[classCnt], alpha );
+    t.stop();
+   
+    alphas.insert( std::pair<int, NICE::Vector> ( classCnt, alpha) );
+  }  
+  
+  // save the parameter value and alpha vectors
+  ikm->getParameters ( min_parameter );
+  this->min_alphas = alphas;
+}
 
 double GPLikelihoodApprox::evaluate(const OPTIMIZATION::matrix_type & x)
 {

GPLikelihoodApprox.h

@@ -103,6 +103,16 @@ class GPLikelihoodApprox : public OPTIMIZATION::CostFunction
     virtual ~GPLikelihoodApprox();
      
     // ------ main methods ------
+    
+    /**
+    * @brief Compute alpha vectors for given hyperparameters
+    *
+    * @param x vector with specified hyperparameters to evaluate their likelihood
+    *
+    * @return void
+    */    
+    void computeAlphaDirect(const OPTIMIZATION::matrix_type & x);
+    
     /**
     * @brief Evaluate the likelihood for given hyperparameters
     *

progs/classifyDatasetGPHIK.cpp

@@ -0,0 +1,185 @@
+/** 
+* @file classifyDatasetGPHIK.cpp
+* @brief Perform classification on an arbitrary dataset with the GPHIK-Classifier
+* @author Alexander Freytag
+* @date 16-09-2013
+*/
+
+// STL-includes
+#include <iostream>
+#include <vector>
+
+// NICE-core includes
+#include <core/basics/Config.h>
+#include <core/basics/Exception.h>
+#include <core/vector/MatrixT.h>
+#include <core/vector/SparseVectorT.h>
+
+// gp-hik-core includes
+#include "gp-hik-core/GPHIKClassifier.h"
+
+
+void readSparseExamples ( const std::string & fn,  std::vector< NICE::SparseVector * > & examples, NICE::Vector & labels )
+{
+  // initially cleaning of variables
+  examples.clear();
+  labels.clear();
+  
+  std::vector<double> labels_std;
+  labels_std.clear();
+  
+  std::cerr << "Reading " << fn << std::endl;
+  std::ifstream ifs ( fn.c_str(), std::ios::in );
+  if ( ! ifs.good() )
+  {
+      std::cerr <<  "Unable to read " << fn << std::endl;
+      return;
+  }
+  
+  // read until no new line is in the file anymore
+  while ( !ifs.eof() )
+  {
+    int classno;
+    if ( !(ifs >> classno) )
+      break;
+    
+    labels_std.push_back( classno );
+    
+    NICE::SparseVector *v = new NICE::SparseVector; 
+    /* needed format in every line: 
+     * SVECTOR dimension size index value index value ... END
+     * with 
+     *      SVECTOR   -- starting flag
+     *      dimension -- overall feature dimension and
+     *      size      -- number of non-zero entries for the current feature vector 
+     *      index     -- integer value specifying a non-zero dimension
+     *      value     -- double value specifying the value for the corresp. non-zero dimension
+     *      END       -- ending flag
+     */
+    try
+    {
+      v->restore ( ifs, NICE::SparseVector::FORMAT_INDEX );    
+    }
+    catch ( NICE::Exception excep)
+    {
+      std::cerr << "Error while reading features. Error message: " << excep.what() << std::endl;
+      break;
+    }
+          
+    
+    examples.push_back ( v );
+  }
+  ifs.close();
+  
+  labels = NICE::Vector( labels_std );
+}
+
+void mapClassNumbersToIndices( const NICE::Vector & labels, std::map<int,int> & mapClassNoToIdx )
+{
+  mapClassNoToIdx.clear();
+  int classCnt ( 0 );
+  
+  for ( NICE::Vector::const_iterator it_labels = labels.begin(); it_labels != labels.end(); it_labels++ )
+  {
+    if ( mapClassNoToIdx.find( *it_labels ) == mapClassNoToIdx.end() )
+    {
+      mapClassNoToIdx.insert( std::pair<int,int>( (int) round(*it_labels), classCnt ) );
+      classCnt++;
+    }
+  }
+}
+
+
+int main (int argc, char* argv[])
+{  
+  std::set_terminate(__gnu_cxx::__verbose_terminate_handler);
+
+  NICE::Config conf ( argc, argv );
+ 
+  NICE::GPHIKClassifier classifier ( &conf, "GPHIKClassifier" );  
+  
+  // ========================================================================
+  //                            TRAINING STEP
+  // ========================================================================  
+   
+  // read training data
+  std::vector< NICE::SparseVector * > examplesTrain;
+  NICE::Vector labelsTrain;
+  
+  std::string s_fn_trainingSet = conf.gS("main", "trainset");
+  readSparseExamples ( s_fn_trainingSet, examplesTrain, labelsTrain );
+
+  //map the occuring classes to a minimal set of indices
+  std::map<int,int> map_classNoToClassIdx_train; // < classNo, Idx>
+  
+  mapClassNumbersToIndices( labelsTrain, map_classNoToClassIdx_train );
+  
+  //how many different classes do we have in the training set?
+  int i_noClassesTrain ( map_classNoToClassIdx_train.size() );
+
+  // train GPHIK classifier
+  classifier.train ( examplesTrain, labelsTrain );
+  
+  // ========================================================================
+  //                            TEST STEP
+  // ========================================================================
+  
+  // read test data
+  std::vector< NICE::SparseVector * > examplesTest;
+  NICE::Vector labelsTest;
+  
+  std::string s_fn_testSet = conf.gS("main", "testset");
+  readSparseExamples ( s_fn_testSet, examplesTest, labelsTest );
+  
+  //map the occuring classes to a minimal set of indices
+  std::map<int,int> map_classNoToClassIdx_test; // < classNo, Idx>
+  
+  mapClassNumbersToIndices( labelsTest, map_classNoToClassIdx_test );
+
+  //how many different classes do we have in the test set?  
+  int i_noClassesTest ( map_classNoToClassIdx_test.size() );
+  
+  // evaluate GPHIK classifier on unseen test data
+  int idx ( 0 );
+  NICE::SparseVector scores;  /* not needed in this evaluation, so we just declare it ones */
+  
+  NICE::Matrix confusion ( i_noClassesTest, i_noClassesTrain, 0.0 );
+  
+  for (std::vector< NICE::SparseVector *>::const_iterator itTestExamples = examplesTest.begin(); itTestExamples != examplesTest.end(); itTestExamples++, idx++)
+  {
+    int classno_groundtruth = labelsTest( idx );
+    int classno_predicted;
+
+    classifier.classify ( *itTestExamples, classno_predicted, scores /* not needed anyway in that evaluation*/ );
+    
+    
+    int idx_classno_groundtruth ( map_classNoToClassIdx_test[ classno_groundtruth ] );
+    int idx_classno_predicted ( map_classNoToClassIdx_train[ classno_predicted ] );
+        
+    confusion( idx_classno_groundtruth, idx_classno_predicted ) += 1;
+  }
+
+
+  confusion.normalizeRowsL1();
+  std::cerr << confusion << std::endl;
+
+  std::cerr << "average recognition rate: " << confusion.trace()/confusion.rows() << std::endl;
+  
+  
+  // ========================================================================
+  //                           clean up memory
+  // ========================================================================
+  
+  // release memore of feature vectors from training set
+  for (std::vector< NICE::SparseVector *>::const_iterator itTrainExamples = examplesTrain.begin(); itTrainExamples != examplesTrain.end(); itTrainExamples++ )
+  {
+    delete *itTrainExamples;
+  }
+  
+  // release memore of feature vectors from test set
+  for (std::vector< NICE::SparseVector *>::const_iterator itTestExamples = examplesTest.begin(); itTestExamples != examplesTest.end(); itTestExamples++ )
+  {
+    delete *itTestExamples;
+  }
+  return 0;
+}