source: trunk/test/normalization.cc @ 3342

// $Id: normalization.cc 3342 2014-11-06 05:26:24Z peter $

/*
  Copyright (C) 2008, 2009 Jari Häkkinen, Peter Johansson
  Copyright (C) 2010, 2012, 2014 Peter Johansson

  This file is part of the yat library, http://dev.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 3 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 yat. If not, see <http://www.gnu.org/licenses/>.
*/

#include <config.h>

#include "Suite.h"

#include "yat/normalizer/Centralizer.h"
#include "yat/normalizer/ColumnNormalizer.h"
#include "yat/normalizer/Gauss.h"
#include "yat/normalizer/qQuantileNormalizer.h"
#include "yat/normalizer/QuantileNormalizer.h"
#include "yat/normalizer/RowNormalizer.h"
#include "yat/normalizer/Spearman.h"
#include "yat/normalizer/Zscore.h"

#include "yat/utility/DataIterator.h"
#include "yat/utility/FileUtil.h"
#include "yat/utility/Matrix.h"
#include "yat/utility/MatrixWeighted.h"
#include "yat/utility/WeightIterator.h"

#include <boost/concept_archetype.hpp>
#include <boost/iterator/iterator_archetypes.hpp>

#include <climits>
#include <fstream>
#include <iterator>
#include <limits>
#include <vector>

using namespace theplu::yat;
void test_centralizer(test::Suite&);
void test_column_normalize(test::Suite&);
void test_gauss_normalize(test::Suite&);
void test_qquantile_normalize(test::Suite&);
void test_qquantile_normalize_weighted(test::Suite&);
void test_quantile_normalize(test::Suite&);
void test_row_normalize(test::Suite&);
void test_spearman(test::Suite&);
void test_spearman_weighted(test::Suite&);
void test_z_score(test::Suite&);

int main(int argc, char* argv[])
{
  test::Suite suite(argc, argv);
  suite.err() << "testing normalizations ... " << std::endl;

  test_centralizer(suite);
  test_column_normalize(suite);
  test_qquantile_normalize(suite);
  test_qquantile_normalize_weighted(suite);
  test_quantile_normalize(suite);
  test_gauss_normalize(suite);
  test_row_normalize(suite);
  test_spearman(suite);
  test_z_score(suite);

  return suite.return_value();
}


void test_centralizer(test::Suite& suite)
{
  suite.err() << "Testing Centralizer\n";
  std::vector<double> vec;
  vec.push_back(1);
  vec.push_back(2);
  vec.push_back(3);
  normalizer::Centralizer<> c;
  c(vec.begin(), vec.end(), vec.begin());
  for (size_t i=0; i<vec.size(); ++i)
    suite.add(suite.equal(vec[i], static_cast<double>(i)-1.0));

  std::vector<utility::DataWeight> vec2;
  vec2.push_back(utility::DataWeight(1,1));
  vec2.push_back(utility::DataWeight(2,0.5));
  vec2.push_back(utility::DataWeight(2,0.5));
  std::vector<utility::DataWeight> vec3(vec2.size());
  c(vec2.begin(), vec2.end(), vec3.begin());
  for (size_t i=0; i<vec2.size(); ++i)
    suite.add(suite.equal(vec3[i].weight(), vec2[i].weight()));
  suite.add(suite.equal(vec3[0].data(), -0.5));
  suite.add(suite.equal(vec3[1].data(), 0.5));
  suite.add(suite.equal(vec3[2].data(), 0.5));

  // compile test should not be run
  if (false) {

    using boost::detail::dummy_constructor;
    c(boost::forward_iterator_archetype<double>(),
      boost::forward_iterator_archetype<double>(),
      boost::mutable_forward_iterator_archetype<double>());

    c(boost::forward_iterator_archetype<double>(),
      boost::forward_iterator_archetype<double>(),
      boost::mutable_forward_iterator_archetype<utility::DataWeight>());

    c(boost::forward_iterator_archetype<utility::DataWeight>(),
      boost::forward_iterator_archetype<utility::DataWeight>(),
      boost::mutable_forward_iterator_archetype<double>());

    c(boost::forward_iterator_archetype<utility::DataWeight>(),
      boost::forward_iterator_archetype<utility::DataWeight>(),
      boost::mutable_forward_iterator_archetype<utility::DataWeight>());
  }
}


void test_column_normalize(test::Suite& suite)
{
  using namespace normalizer;
  suite.err() << "Testing ColumnNormalizer\n";

  utility::Matrix m(2,2);
  m(0,0) = 0;
  m(0,1) = 10;
  m(1,0) = 2;
  m(1,1) = 4;
  ColumnNormalizer<Centralizer<> > qn;
  qn(m, m);
  suite.err() << "Testing m(0,0)\n";
  suite.add(suite.equal(m(0,0), -1));
  suite.err() << "Testing m(0,1)\n";
  suite.add(suite.equal(m(0,1), 3));
  suite.err() << "Testing m(1,0)\n";
  suite.add(suite.equal(m(1,0), 1));
  suite.err() << "Testing m(1,1)\n";
  suite.add(suite.equal(m(1,1), -3));

  if (false) { // do not run compile tests
    test::container2d_archetype<double> container2d; 
    test::mutable_container2d_archetype<double> mutable_container2d; 
    qn(container2d, mutable_container2d);
  }
}


void test_qquantile_normalize(test::Suite& suite)
{
  using namespace normalizer;

  suite.err() << "Testing qQuantileNormalizer\n";
  std::string data(test::filename("data/normalization_test.data"));
  if (utility::FileUtil(data.c_str()).permissions("r")) {
    suite.add(false);
    suite.err() << "Cannot access file " << data << '\n';
    return;
  }
  std::ifstream data_stream(data.c_str());

  utility::Matrix m(data_stream);

  suite.err() << "testing number of parts (Q) boundary conditions\n";
  qQuantileNormalizer(m.begin_column(0), m.end_column(0), m.rows());
  qQuantileNormalizer(m.begin_column(0), m.end_column(0), 3);

  // first column as target
  qQuantileNormalizer qqn(m.begin_column(0), m.end_column(0) ,9);  
  ColumnNormalizer<qQuantileNormalizer> cn(qqn);
  utility::Matrix result(m.rows(),m.columns());
  cn(m, result);

  suite.err() << "test that result can be stored in the source matrix...";
  cn(m,m);
  if (suite.add(result==m))
    suite.err() << " ok.\n";
  else 
    suite.err() << " failed.\n";

  // Enough iteration will make all columns to have the same values as
  // the target.
  suite.err() << "Testing that q=matrix rows gives QuantileNormalization\n";
  utility::Matrix m2(4,2);
  m2(0,0) = 0; m2(0,1) = 10;
  m2(1,0) = 2; m2(1,1) = 4;
  m2(2,0) = 1; m2(2,1) = 0;
  m2(3,0) = 3; m2(3,1) = 7;
  qQuantileNormalizer qqn2(m2.begin_column(0), m2.end_column(0), m2.rows());
  ColumnNormalizer<qQuantileNormalizer> cn2(qqn2);
  utility::Matrix result2(m2.rows(),m2.columns());
  cn2(m2,result2);
  suite.add( suite.equal_fix(m2(0,0),result2(2,1),1.0e-12) &&
             suite.equal_fix(m2(1,0),result2(3,1),1.0e-12) &&
             suite.equal_fix(m2(2,0),result2(1,1),1.0e-12) &&
             suite.equal_fix(m2(3,0),result2(0,1),1.0e-12) );
  // compile test should not be run
  if (false) {
    qQuantileNormalizer qqn3(boost::forward_iterator_archetype<double>(),
                             boost::forward_iterator_archetype<double>(),
                             100);
    qqn3(boost::random_access_iterator_archetype<double>(), 
         boost::random_access_iterator_archetype<double>(),
         boost::mutable_random_access_iterator_archetype<double>());
    using utility::DataWeight;
    qQuantileNormalizer qqn4(boost::forward_iterator_archetype<DataWeight>(),
                             boost::forward_iterator_archetype<DataWeight>(),
                             100);
    qqn4(boost::random_access_iterator_archetype<DataWeight>(), 
         boost::random_access_iterator_archetype<DataWeight>(),
         boost::mutable_random_access_iterator_archetype<DataWeight>());
    qqn4(boost::random_access_iterator_archetype<DataWeight>(), 
         boost::random_access_iterator_archetype<DataWeight>(),
         boost::mutable_random_access_iterator_archetype<double>());
  }
}


void test_qquantile_normalize_weighted(test::Suite& suite)
{
  using namespace normalizer;

  suite.err() << "Testing qQuantileNormalizer weighted\n";

  // test with unweighted target and source
  std::vector<double> target;
  target.reserve(1000);
  while (target.size()<1000)
    target.push_back(target.size());
  qQuantileNormalizer qQN(target.begin(), target.end(), 4);
  std::vector<double> source;
  while (source.size()<10)
    source.push_back(source.size()*10);
  std::vector<double> result(source.size());
  
  qQN(source.begin(), source.end(), result.begin());
  
  using utility::DataWeight;
  suite.err() << "Testing with unweighted target and weighted source\n";
  std::vector<utility::DataWeight> src_w(source.size(), DataWeight(0, 1));
  std::copy(source.begin(), source.end(),
            utility::data_iterator(src_w.begin()));

  std::vector<utility::DataWeight> result_w(src_w.size());
  qQN(src_w.begin(), src_w.end(), result_w.begin());
  suite.add(suite.equal_range(result.begin(), result.end(),
                              utility::data_iterator(result_w.begin())));

  suite.err() << "Testing with missing value in source\n";
  // adding a missing value
  std::vector<utility::DataWeight>::iterator MWi=src_w.begin();
  MWi+=5;
  src_w.insert(MWi, DataWeight(std::numeric_limits<double>::quiet_NaN(), 0.0));
  std::vector<utility::DataWeight> result_w2(src_w.size());
  qQN(src_w.begin(), src_w.end(), result_w2.begin());
  // excluding missing value from comparison in suite.equal_range
  MWi=result_w2.begin();
  MWi+=5;
  result_w2.erase(MWi);
  suite.add(suite.equal_range(utility::data_iterator(result_w.begin()), 
                              utility::data_iterator(result_w.end()),
                              utility::data_iterator(result_w2.begin())));

  suite.err() << "testing with weighted target" << std::endl;
  std::vector<utility::DataWeight> target_w(target.size()+1, DataWeight(0, 1));
  target_w[0] = DataWeight(5.3, 0);
  std::copy(target.begin(), target.end(),
            utility::data_iterator(target_w.begin()+1));
  qQuantileNormalizer qQNw(target_w.begin(), target_w.end(), 4);
  std::vector<utility::DataWeight> result_w3(src_w.size());
  qQNw(src_w.begin(), src_w.end(), result_w3.begin());
  // excluding missing value from comparison in suite.equal_range
  MWi=result_w3.begin();
  MWi+=5;
  result_w3.erase(MWi);
  suite.add(suite.equal_range(utility::data_iterator(result_w3.begin()), 
                              utility::data_iterator(result_w3.end()),
                              utility::data_iterator(result_w2.begin())));
  
}


void test_quantile_normalize(test::Suite& suite)
{
  suite.err() << "Testing QuantileNormalizer\n";
  utility::Matrix m(2,2);
  m(0,0) = 0;
  m(0,1) = 10;
  m(1,0) = 2;
  m(1,1) = 4;
  normalizer::QuantileNormalizer qn;
  qn(m, m);
  suite.err() << "Testing m(0,0)\n";
  suite.add(suite.equal(m(0,0), 2));
  suite.err() << "Testing m(0,1)\n";
  suite.add(suite.equal(m(0,1), 6));
  suite.err() << "Testing m(1,0)\n";
  suite.add(suite.equal(m(1,0), 6));
  suite.err() << "Testing m(1,1)\n";
  suite.add(suite.equal(m(1,1), 2));

  // testing with landscape-shaped data matrix
  m.resize(2, 10);
  qn(m, m);
}

void test_row_normalize(test::Suite& suite)
{
  using namespace normalizer;
  suite.err() << "Testing RowNormalizer\n";
  
  utility::Matrix m(2,3);
  m(0,0) = 0;
  m(0,1) = 10;
  m(1,0) = 2;
  m(1,1) = 4;
  utility::Matrix m2(m);
  m2.transpose();
  ColumnNormalizer<Centralizer<> > cn;
  RowNormalizer<Centralizer<> > rn;
  cn(m, m);
  rn(m2, m2);
  m2.transpose();
  suite.equal_range(m.begin(), m.end(), m2.begin());
  if (false) { // do not run compile tests
    test::container2d_archetype<double> container2d; 
    test::mutable_container2d_archetype<double> mutable_container2d; 
    rn(container2d, mutable_container2d);
  }
}

void test_spearman(test::Suite& suite)
{
  suite.err() << "Testing Spearman\n";
  normalizer::Spearman spearman;
  std::vector<double> vec;
  vec.push_back(0);
  vec.push_back(2);
  vec.push_back(3);
  vec.push_back(1);
  spearman(vec.begin(), vec.end(), vec.begin());
  std::vector<double> correct;
  correct.push_back(1.0/8);
  correct.push_back(5.0/8);
  correct.push_back(7.0/8);
  correct.push_back(3.0/8);
  suite.add(suite.equal_range(vec.begin(), vec.end(), correct.begin()));
  suite.err() << "Testing Spearman with ties\n";
  vec[1]=vec[2];
  correct[1] = correct[2] = (correct[1]+correct[2])/2;
  spearman(vec.begin(), vec.end(), vec.begin());
  suite.add(suite.equal_range(vec.begin(), vec.end(), correct.begin()));
  test_spearman_weighted(suite);
}


void test_gauss_normalize(test::Suite& suite)
{
  suite.err() << "Testing Gauss\n";
  normalizer::Gauss gauss;
  std::vector<double> vec;
  vec.push_back(1);
  gauss(vec.begin(), vec.end(), vec.begin());
  suite.add(suite.equal(vec.front(), 0));
  vec.push_back(1);
  gauss(vec.begin(), vec.end(), vec.begin());
  suite.add(suite.equal(vec.front(), -vec.back()));
  // compile test should not be run
  if (false) {
    using boost::iterator_archetype;
    iterator_archetype<double,
                       boost::iterator_archetypes::readable_iterator_t,
                       boost::random_access_traversal_tag>
      readable_iterator;

    iterator_archetype<double,
                       boost::iterator_archetypes::readable_writable_iterator_t,
                       boost::random_access_traversal_tag>
      writable_iterator;

    // We have to use lvalue here because otherwise proxy classes
    // provided by boost kick in and they do not provide the needed
    // data() and weight() functions that e.g. DataWeightProxy does.
    iterator_archetype<utility::DataWeight,
                       boost::iterator_archetypes::readable_lvalue_iterator_t,
                       boost::random_access_traversal_tag>
      weighted_readable_iterator;

    iterator_archetype<utility::DataWeight,
                       boost::iterator_archetypes::writable_lvalue_iterator_t,
                       boost::random_access_traversal_tag>
      weighted_writable_iterator;

    gauss(readable_iterator, readable_iterator, writable_iterator);
    gauss(readable_iterator, readable_iterator, weighted_writable_iterator);
    gauss(weighted_readable_iterator, weighted_readable_iterator,
          writable_iterator);
    gauss(weighted_readable_iterator, weighted_readable_iterator,
          weighted_writable_iterator);
  }
}

void test_spearman_weighted(test::Suite& suite)
{
  suite.err() << "Testing Weighted Spearman\n";
  normalizer::Spearman spearman;

  suite.err() << "Testing that unity weights reproduces unweighted case\n";
  utility::MatrixWeighted m(1,4,0,1);
  utility::MatrixWeighted res(m.rows(), m.columns(),3.14,0);
  m(0,0).data()=0;
  m(0,1).data()=2;
  m(0,2).data()=3;
  m(0,3).data()=1;
  std::vector<double> correct(m.columns());
  std::vector<double> correct_w(m.columns(), 1.0);
  std::copy(utility::data_iterator(m.begin_row(0)),
            utility::data_iterator(m.end_row(0)),
            correct.begin());
  spearman(correct.begin(), correct.end(), correct.begin());
  spearman(m.begin_row(0), m.end_row(0), res.begin_row(0));

  using utility::data_iterator;
  suite.add(suite.equal_range(data_iterator(res.begin_row(0)),
                               data_iterator(res.end_row(0)),
                               correct.begin()));
  using utility::weight_iterator;
  suite.add(suite.equal_range(weight_iterator(res.begin_row(0)),
                               weight_iterator(res.end_row(0)),
                               correct_w.begin()));

  suite.err() << "Testing rescaling of weights\n";
  for (size_t i=0; i<m.columns(); ++i) {
    m(0,i).weight() *= 2;
    correct_w[i] *= 2;
  }   
  spearman(m.begin_row(0), m.end_row(0), res.begin_row(0));
  suite.add(suite.equal_range(data_iterator(res.begin_row(0)),
                               data_iterator(res.end_row(0)),
                               correct.begin()));
  suite.add(suite.equal_range(weight_iterator(res.begin_row(0)),
                               weight_iterator(res.end_row(0)),
                               correct_w.begin()));

  
  suite.err() << "Testing case with a zero weight\n";
  m(0,1).data() = std::numeric_limits<double>::quiet_NaN();
  m(0,1).weight() = 0.0;
  spearman(m.begin_row(0), m.end_row(0), res.begin_row(0));
  suite.add(suite.equal(res(0,0).data(), 0.5/3)); 
  suite.add(suite.equal(res(0,2).data(), 2.5/3)); 
  suite.add(suite.equal(res(0,3).data(), 1.5/3)); 

  suite.err() << "Testing case with ties\n";
  m(0,0).data() = m(0,2).data();
  spearman(m.begin_row(0), m.end_row(0), res.begin_row(0));
  suite.add(suite.equal(res(0,0).data(), 2.0/3)); 
  suite.add(suite.equal(res(0,2).data(), 2.0/3)); 
  suite.add(suite.equal(res(0,3).data(), 0.5/3)); 
  // compile test should not be run
  if (false) {
    using boost::iterator_archetype;
    iterator_archetype<double,
                       boost::iterator_archetypes::readable_iterator_t,
                       boost::random_access_traversal_tag>
      readable_iterator;

    iterator_archetype<double,
                       boost::iterator_archetypes::readable_writable_iterator_t,
                       boost::random_access_traversal_tag>
      writable_iterator;

    // We have to use lvalue here because otherwise proxy classes
    // provided by boost kick in and they do not provide the needed
    // data() and weight() functions that e.g. DataWeightProxy does.
    iterator_archetype<utility::DataWeight,
                       boost::iterator_archetypes::readable_lvalue_iterator_t,
                       boost::random_access_traversal_tag>
      weighted_readable_iterator;

    iterator_archetype<utility::DataWeight,
                       boost::iterator_archetypes::writable_lvalue_iterator_t,
                       boost::random_access_traversal_tag>
      weighted_writable_iterator;

    spearman(readable_iterator, readable_iterator, writable_iterator);
    spearman(readable_iterator, readable_iterator, weighted_writable_iterator);
    spearman(weighted_readable_iterator, weighted_readable_iterator,
             writable_iterator);
    spearman(weighted_readable_iterator, weighted_readable_iterator,
             weighted_writable_iterator);
  }

}

void test_z_score(test::Suite& suite)
{
  suite.err() << "Testing Zscore\n";
  std::vector<double> vec;
  vec.push_back(0);
  vec.push_back(3.14);
  normalizer::Zscore zscore;
  zscore(vec.begin(), vec.end(), vec.begin());
  for (size_t i=0; i<vec.size(); ++i)
    suite.add(suite.equal(vec[i], 2.0*i-1.0));

  std::vector<utility::DataWeight> vec2;
  vec2.push_back(utility::DataWeight(1,1));
  vec2.push_back(utility::DataWeight(2.13,0.5));
  vec2.push_back(utility::DataWeight(2.13,0.5));
  std::vector<utility::DataWeight> vec3(vec2.size());
  zscore(vec2.begin(), vec2.end(), vec3.begin());
  for (size_t i=0; i<vec2.size(); ++i)
    suite.add(suite.equal(vec3[i].weight(), vec2[i].weight()));
  suite.add(suite.equal(vec3[0].data(), -1.0));
  suite.add(suite.equal(vec3[1].data(), 1.0));
  suite.add(suite.equal(vec3[2].data(), 1.0));
  // compile test should not be run
  if (false) {
    boost::detail::dummy_constructor dummy_cons;
    zscore(boost::forward_iterator_archetype<double>(), 
           boost::forward_iterator_archetype<double>(),
           boost::output_iterator_archetype<double>(dummy_cons));
    using utility::DataWeight;
    zscore(boost::forward_iterator_archetype<DataWeight>(), 
           boost::forward_iterator_archetype<DataWeight>(),
           boost::mutable_forward_iterator_archetype<DataWeight>());
    zscore(boost::random_access_iterator_archetype<DataWeight>(), 
           boost::random_access_iterator_archetype<DataWeight>(),
           boost::mutable_random_access_iterator_archetype<double>());
    zscore(boost::random_access_iterator_archetype<double>(), 
           boost::random_access_iterator_archetype<double>(),
           boost::mutable_random_access_iterator_archetype<DataWeight>());
  }
}