source: trunk/test/ncc_test.cc @ 865

// $Id: ncc_test.cc 865 2007-09-10 19:41:04Z peter $

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

  This file is part of the yat library, http://trac.thep.lu.se/trac/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/MatrixLookup.h"
#include "yat/classifier/MatrixLookupWeighted.h"
#include "yat/classifier/NCC.h"
#include "yat/classifier/Target.h"
#include "yat/utility/matrix.h"
#include "yat/statistics/PearsonDistance.h"
#include "yat/utility/utility.h"

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

using namespace theplu::yat;

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;

  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);
  for(size_t i=0;i<data.rows();++i)
    for(size_t j=0;j<data.columns();++j)
      if(!std::isnan(data(i,j)))
        weights(i,j)=1.0;
      
  classifier::MatrixLookupWeighted dataviewweighted(data,weights);
  statistics::PearsonDistance pearson;  
  classifier::NCC ncc(dataviewweighted,targets,pearson);
  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;
  }

  double slack = 0;
  for (size_t i=0; i<centroids.rows(); i++){
    for (size_t j=0; j<centroids.columns(); j++){
      slack += fabs(centroids(i,j)-ncc.centroids()(i,j));
    }
  }
  slack /= (centroids.columns()*centroids.rows());
  double slack_bound=2e-7;
  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;
  }

  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 = 0;
  for (size_t i=0; i<result.rows(); i++){
    for (size_t j=0; j<result.columns(); j++){
        slack += fabs(result(i,j)-prediction(i,j));
    }
  }
  slack /= (result.columns()*result.rows());
  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;
  }
  

  // Testing IGP  
  classifier:: MatrixLookup dataview(data);
  *error << "testing igp" << std::endl;
  classifier::IGP igp(dataview,targets,pearson);
  *error << igp.score() << std::endl;

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


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

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