source: trunk/test/ncc_test.cc @ 1142

// $Id: ncc_test.cc 1142 2008-02-25 14:32:35Z markus $

/*
  Copyright (C) 2006 Jari Häkkinen, Markus Ringnér

  This file is part of the yat library, http://trac.thep.lu.se/yat

  The yat library 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.

  The yat library 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., 59 Temple Place - Suite 330, Boston, MA
  02111-1307, USA.
*/

#include "yat/classifier/IGP.h"
#include "yat/classifier/Kernel_MEV.h"
#include "yat/classifier/KernelLookup.h"
#include "yat/classifier/MatrixLookup.h"
#include "yat/classifier/MatrixLookupWeighted.h"
#include "yat/classifier/NCC.h"
#include "yat/classifier/PolynomialKernelFunction.h"
#include "yat/classifier/Target.h"
#include "yat/utility/Matrix.h"
#include "yat/statistics/EuclideanDistance.h"
#include "yat/statistics/PearsonDistance.h"
#include "yat/utility/utility.h"

#include <cassert>
#include <fstream>
#include <iostream>
#include <stdexcept>
#include <sstream>
#include <string>
#include <limits>
#include <cmath>

using namespace theplu::yat;

double deviation(const utility::Matrix& a, const utility::Matrix& b) {
  double sl=0;
  for (size_t i=0; i<a.rows(); i++){
    for (size_t j=0; j<a.columns(); j++){
      sl += fabs(a(i,j)-b(i,j));
    }
  }
  sl /= (a.columns()*a.rows());
  return sl;
}

int main(const int argc,const char* argv[])
{  

  std::ostream* error;
  if (argc>1 && argv[1]==std::string("-v"))
    error = &std::cerr;
  else {
    error = new std::ofstream("/dev/null");
    if (argc>1)
      std::cout << "ncc_test -v : for printing extra information\n";
  }
  *error << "testing ncc" << std::endl;
  bool ok = true;

  /////////////////////////////////////////////
  // First test of constructor and training  
  /////////////////////////////////////////////
  classifier::MatrixLookup ml(4,4);
  std::vector<std::string> vec(4, "pos");
  vec[3]="bjds";
  classifier::Target target(vec);
  classifier::NCC<statistics::EuclideanDistance> ncctmp(ml,target);
  *error << "training...\n";
  ncctmp.train();
  *error << "done\n";

  /////////////////////////////////////////////
  // A test of predictions using unweighted data
  /////////////////////////////////////////////
  *error << "test of predictions using unweighted test data\n";
  utility::Matrix data1(3,4);
  for(size_t i=0;i<3;i++) {
    data1(i,0)=3-i;
    data1(i,1)=5-i;
    data1(i,2)=i+1;
    data1(i,3)=i+3;
  }
  std::vector<std::string> vec1(4, "pos");
  vec1[0]="neg";
  vec1[1]="neg";

  classifier::MatrixLookup ml1(data1);
  classifier::Target target1(vec1);

  classifier::NCC<statistics::EuclideanDistance> ncc1(ml1,target1);
  ncc1.train();
  utility::Matrix prediction1;
  ncc1.predict(ml1,prediction1);
  double slack_bound=2e-7;
  utility::Matrix result1(2,4);
  result1(0,0)=result1(0,1)=result1(1,2)=result1(1,3)=sqrt(3.0);
  result1(0,2)=result1(0,3)=result1(1,0)=result1(1,1)=sqrt(11.0);
  double slack = deviation(prediction1,result1); 
  if (slack > slack_bound || std::isnan(slack)){
    *error << "Difference to expected prediction too large\n";
    *error << "slack: " << slack << std::endl;
    *error << "expected less than " << slack_bound << std::endl;
    ok = false;
  }

  //////////////////////////////////////////////////////////////////////////
  // A test of predictions using unweighted training and weighted test data
  //////////////////////////////////////////////////////////////////////////
  *error << "test of predictions using unweighted training and weighted test data\n";
  utility::Matrix weights1(3,4,1.0);
  weights1(0,0)=weights1(1,1)=weights1(2,2)=weights1(1,3)=0.0;
  classifier::MatrixLookupWeighted mlw1(data1,weights1);
  ncc1.predict(mlw1,prediction1); 
  result1(0,2)=result1(0,3)=result1(1,0)=result1(1,1)=sqrt(15.0);
  slack = deviation(prediction1,result1);
  if (slack > slack_bound || std::isnan(slack)){
    *error << "Difference to expected prediction too large\n";
    *error << "slack: " << slack << std::endl;
    *error << "expected less than " << slack_bound << std::endl;
    ok = false;
  }

  //////////////////////////////////////////////////////////////////////////
  // A test of predictions using weighted training resulting in NaN's
  // in centroids and unweighted test data
  //////////////////////////////////////////////////////////////////////////
  *error << "test of predictions using nan centroids and unweighted test data\n";
  utility::Matrix weights2(3,4,1.0);
  weights2(1,0)=weights2(1,1)=0.0;
  classifier::MatrixLookupWeighted mlw2(data1,weights2);
  classifier::NCC<statistics::EuclideanDistance> ncc2(mlw2,target1);
  ncc2.train();
  ncc2.predict(ml1,prediction1);  
  result1(0,0)=result1(0,1)=result1(1,2)=result1(1,3)=sqrt(3.0);
  result1(1,0)=result1(1,1)=sqrt(11.0);
  result1(0,2)=result1(0,3)=sqrt(15.0);
  slack = deviation(prediction1,result1);
  if(!std::isnan(ncc2.centroids()(1,0))) 
    ok=false;
  if (slack > slack_bound || std::isnan(slack)){
    *error << "Difference to expected prediction too large\n";
    *error << "slack: " << slack << std::endl;
    *error << "expected less than " << slack_bound << std::endl;
    ok = false;
  }

  //////////////////////////////////////////////////////////////////////////
  // A test of predictions when a centroid has nan for all variables that a
  // test sample has non-zero weights for.
  //////////////////////////////////////////////////////////////////////////
  *error << "test of predictions using nan centroids and weighted test data\n";
  *error << "... using EuclideanDistance" << std::endl;
  weights1(0,0)=weights1(2,0)=0;
  classifier::NCC<statistics::EuclideanDistance> ncc3(mlw2,target1);
  ncc3.train();
  ncc3.predict(mlw1,prediction1); 
  if(!std::isnan(ncc3.centroids()(1,0))) {
    ok=false;
    *error << "Training failed: expected nan in centroid" << std::endl;
  }
  if(!(std::isnan(prediction1(0,0)) &&
       fabs(prediction1(1,0)-sqrt(3.0))<slack_bound &&
       fabs(prediction1(0,1)-sqrt(3.0))<slack_bound &&
       fabs(prediction1(1,1)-sqrt(15.0))<slack_bound &&
       fabs(prediction1(0,2)-sqrt(27.0))<slack_bound)) {  
    ok=false; 
    *error << "Test failed: predictions incorrect" << std::endl;
  }
  *error << "... using PearsonDistance" << std::endl;;
  classifier::NCC<statistics::PearsonDistance> ncc4(mlw2,target1);
  ncc4.train();
  ncc4.predict(mlw1,prediction1); 
  if(!std::isnan(ncc4.centroids()(1,0))) {
    ok=false;
    *error << "Training failed: expected nan in centroid" << std::endl;
  }
  if(!(std::isnan(prediction1(0,0)) &&
       std::isnan(prediction1(0,2)) &&
       std::isnan(prediction1(1,0)) &&
       fabs(prediction1(0,1))<slack_bound &&
       fabs(prediction1(1,2))<slack_bound &&
       fabs(prediction1(1,3))<slack_bound &&  
       fabs(prediction1(0,3)-2.0)<slack_bound &&
       fabs(prediction1(1,1)-2.0)<slack_bound)) {
    ok=false; 
    *error << "Test failed: predictions incorrect" << std::endl;
  }


  //////////////////////////////////////////////////////////////////////////
  // A test of predictions using Sorlie data
  //////////////////////////////////////////////////////////////////////////
  *error << "test with Sorlie data\n";
  std::ifstream is("data/sorlie_centroid_data.txt");
  utility::Matrix data(is,'\t');
  is.close();

  is.open("data/sorlie_centroid_classes.txt");
  classifier::Target targets(is);
  is.close();

  // Generate weight matrix with 0 for missing values and 1 for others.
  utility::Matrix weights(data.rows(),data.columns(),0.0);
  utility::nan(data,weights);
      
  classifier::MatrixLookupWeighted dataviewweighted(data,weights);
  classifier::NCC<statistics::PearsonDistance> ncc(dataviewweighted,targets);
  *error << "training...\n";
  ncc.train();

  // Comparing the centroids to stored result
  is.open("data/sorlie_centroids.txt");
  utility::Matrix centroids(is);
  is.close();

  if(centroids.rows() != ncc.centroids().rows() ||
     centroids.columns() != ncc.centroids().columns()) {
    *error << "Error in the dimensionality of centroids\n";
    *error << "Nof rows: " << centroids.rows() << " expected: " 
           << ncc.centroids().rows() << std::endl;
    *error << "Nof columns: " << centroids.columns() << " expected: " 
           << ncc.centroids().columns() << std::endl;
  }

  slack = deviation(centroids,ncc.centroids());
  if (slack > slack_bound || std::isnan(slack)){
    *error << "Difference to stored centroids too large\n";
    *error << "slack: " << slack << std::endl;
    *error << "expected less than " << slack_bound << std::endl;
    ok = false;
  }

  *error << "...predicting...\n";
  utility::Matrix prediction;
  ncc.predict(dataviewweighted,prediction);
  
  // Comparing the prediction to stored result
  is.open("data/sorlie_centroid_predictions.txt");
  utility::Matrix result(is,'\t');
  is.close();

  slack = deviation(result,prediction);
  if (slack > slack_bound || std::isnan(slack)){
    *error << "Difference to stored prediction too large\n";
    *error << "slack: " << slack << std::endl;
    *error << "expected less than " << slack_bound << std::endl;
    ok = false;
  }
  *error << "done\n";

  if(ok)
    *error << "OK" << std::endl;
  else
    *error << "FAILED" << std::endl;

  if (error!=&std::cerr)
    delete error;

  if(ok) 
    return 0;
  return -1;  
}