source: trunk/test/matrix_test.cc @ 2119

// $Id: matrix_test.cc 2119 2009-12-12 23:11:43Z peter $

/*
  Copyright (C) 2005 Jari Häkkinen, Peter Johansson
  Copyright (C) 2006 Jari Häkkinen, Peter Johansson, Markus Ringnér
  Copyright (C) 2007, 2008 Jari Häkkinen, 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 "Suite.h"

#include "yat/utility/Matrix.h"

#include <cstdio>
#include <fstream>
#include <iostream>

class matrixwrapper
{
public:
  matrixwrapper(size_t i,size_t j,double value=7)
    : m_(i,j,value) {}

  inline theplu::yat::utility::VectorView
  row(const size_t& i) { return m_.row_view(i); }

  inline const theplu::yat::utility::Matrix& matrix(void) const { return m_; }

private:
  theplu::yat::utility::Matrix m_;
};


int main(int argc, char* argv[])
{
  using namespace theplu::yat;
  test::Suite suite(argc, argv);

  suite.err() << "Testing matrix class" << std::endl;
  utility::Matrix unit3x3(3,3);
  for (size_t i=0; i<unit3x3.rows(); ++i)
    unit3x3(i,i)=1;

  suite.add(suite.test_stream(unit3x3));

  suite.err() << "\tcopy constructor and operator!=" << std::endl;
  utility::Matrix m(3,3,9);
  utility::Matrix m2(m);
  if (m2!=m)
    suite.add(false);

  suite.err() << "\toutput operator and istream constructor" << std::endl;
  // Checking that the matrix output operator writes a file that the
  // input operator can read.
  std::ofstream my_out("/tmp/yat_test_matrix.txt");
  my_out << m2;
  my_out.close();
  std::ifstream is("/tmp/yat_test_matrix.txt");
  utility::Matrix m3(is);
  is.close();
  if (m3!=m2)
    suite.add(false);
  std::remove("/tmp/yat_test_matrix.txt");

  suite.err() << "\toperator*=(double)" << std::endl;
  utility::Matrix m4(3,3,1);
  m4 *= 9;
  if (m4!=m) {
    suite.add(false);
    suite.err() << "error operator*=(double)" << std::endl;
  }

  is.open(test::filename("data/knni_matrix.data").c_str());
  // The stream input is a proper matrix file, with some stray empty
  // lines and other whitespaces. The file is not expected to break
  // things.
  utility::Matrix m5(is);
  is.close();
  double m5_sum=0;
  for (size_t i=0; i<m5.rows(); ++i)
    for (size_t j=0; j<m5.columns(); ++j)
      m5_sum+=m5(i,j);

  // checking that copy constructor creates an independent object when
  // a non-view matrix is copied
  {
    suite.err() << "\tcopy constructor" << std::endl;
    utility::Matrix m2(m5);
    suite.add(m2.rows()==m5.rows());
    suite.add(m2.columns()==m5.columns());
    suite.add(m2==m5);
    suite.add(&m2 != &m5);
  }

  suite.err() << "\tsub-(row)vector" << std::endl;
  // Checking that the row view works, i.e. mutation to the view are
  // reflected in the viewed object.
  double m5_sum2=0;
  for (size_t i=0; i<m5.rows(); ++i)
    for (size_t j=0; j<m5.columns(); ++j)
      m5_sum2+=m5(i,j);
  utility::VectorView v5subrow(m5,3);
  double v5subrow_sum=0;
  for (size_t i=0; i<v5subrow.size(); ++i) {
    v5subrow_sum+=v5subrow(i);
    v5subrow(i)=0;
  }
  double m5_sum3=0;
  for (size_t i=0; i<m5.rows(); ++i)
    for (size_t j=0; j<m5.columns(); ++j)
      m5_sum3+=m5(i,j);
  if (m5_sum3-m5_sum2+v5subrow_sum>1e-13) {
    suite.add(false);
    suite.err() << "error sub-vector test 1" << std::endl;
  }

  suite.err() << "\tsub-(column)vector" << std::endl;
  // Checking that the column view works, i.e. mutation to the view
  // are reflected in the viewed object.
  m5_sum3=0;
  for (size_t i=0; i<m5.rows(); ++i)
    for (size_t j=0; j<m5.columns(); ++j)
      m5_sum3+=m5(i,j);
  utility::VectorView v5subcolumn = m5.column_view(0);
  double v5subcolumn_sum=0;
  for (size_t i=0; i<v5subcolumn.size(); ++i) {
    v5subcolumn_sum+=v5subcolumn(i);
    v5subcolumn(i)=1;
  }
  double m5_sum4=0;
  for (size_t i=0; i<m5.rows(); ++i)
    for (size_t j=0; j<m5.columns(); ++j)
      m5_sum4+=m5(i,j);
  if (m5_sum4-m5_sum3-v5subcolumn.size()+v5subcolumn_sum>1e-13) {
    suite.add(false);
    suite.err() << "error sub-vector test 2" << std::endl;
  }

  // Checking that the column view above mutates the values in the row
  // view.
  double v5subrow_sum2=0;
  for (size_t i=0; i<v5subrow.size(); ++i)
    v5subrow_sum2+=v5subrow(i);
  if (v5subrow_sum2-v5subcolumn(3)>1e-13) {
    suite.add(false);
    suite.err() << "error sub-vector test 3" << std::endl;
  }

  suite.err() << "\tsub-vector and vector copying" << std::endl;
  // Checking that a view is not inherited through the copy
  // contructor.
  utility::Vector v6(v5subrow);
  v6.all(2);
  double v5subrow_sum3=0;
  for (size_t i=0; i<v5subrow.size(); ++i)
    v5subrow_sum3+=v5subrow(i);
  if (v5subrow_sum3-v5subrow_sum2>1e-13) {
    suite.add(false);
    suite.err() << "error sub-vector test 4" << std::endl;
  }
  // Checking that values in a vector is copied into a viewed matrix.
  v5subrow = v6;
  double m5_sum5=0;
  for (size_t i=0; i<m5.rows(); ++i)
    for (size_t j=0; j<m5.columns(); ++j)
      m5_sum5+=m5(i,j);
  if (m5_sum5-m5_sum4-v5subrow.size()*2+v5subrow_sum3>1e-13) {
    suite.add(false);
    suite.err() << "error sub-vector test 5" << std::endl;
  }

  // Checking that the memberfunction matrixwrapper::row() returns a
  // view
  suite.err() << "\tthat class member returns a view" << std::endl;
  matrixwrapper mw(5,2);
  utility::VectorView mwrow=mw.row(2);
  if (mwrow.gsl_vector_p()->data != &(mw.matrix()(2,0))) {
    suite.add(false);
    suite.err() << "error sub-vector test 7" << std::endl;
  }

  suite.err() << "\tmatrix::nan()" << std::endl;
  is.open(test::filename("data/sorlie_centroid_data.txt").c_str());
  utility::Matrix* m_nan = new utility::Matrix(is,'\t');
  utility::Matrix m_weight;
  utility::nan(*m_nan,m_weight);
  is.close();
  if (m_weight(0,0)!=1){
    suite.err() << "error in matrix::nan(): element(0,0) is " << m_weight(0,0)
           << " expected 1." << std::endl;
    suite.add(false);
  }
  if (m_weight(3,2)!=0){
    suite.err() << "error in matrix::nan(): element(3,2) is " << m_weight(3,2)
           << " expected 0." << std::endl;
    suite.add(false);
  }
  delete m_nan;

  suite.err() << "\toperator*=(matrix&)" << std::endl;
  utility::Matrix m6(unit3x3);
  m6 *= m;
  if (m6!=m) {
    suite.add(false);
    suite.err() << "error operator*=(matrix) 1" << std::endl;
  }
  m6 *= unit3x3;
  if (m6!=m) {
    suite.add(false);
    suite.err() << "error operator*=(matrix) 2" << std::endl;
  }
  m6*= utility::Matrix(3,4,1.0);
  m6*= utility::Matrix(4,3,1.0);
  m6*= utility::Matrix(3,5,2.0);
  m6*= utility::Matrix(5,5,2.0);
  m6*= utility::Matrix(5,3,2.0);
  m6*= unit3x3;

  utility::Matrix m7(0,0);
  m7.resize(1,1);
  m7 = utility::Matrix(0,0);
  m7.resize(0,0);
  m7 = utility::Matrix(0,0);

  suite.return_value();
}