source: trunk/test/matrix_test.cc @ 616

// $Id: matrix_test.cc 616 2006-08-31 08:52:02Z jari $

#include <c++_tools/utility/matrix.h>

#include <unistd.h>
#include <fstream>
#include <iostream>

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

  inline theplu::utility::vector
  row(const size_t& i) { return theplu::utility::vector(m_,i); }

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

private:
  theplu::utility::matrix m_;
};


int main(const int argc,const char* argv[])
{
  using namespace theplu;
  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 << "matrix_test -v : for printing extra information\n";
  }

  *error << "Testing matrix class" << std::endl;
  bool ok = true;

  *error << "\tcopy constructor and operator!=" << std::endl;
  utility::matrix m(3,3,9);
  utility::matrix m2(m);
  if (m2!=m)
    ok=false;

  *error << "\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("data/tmp_test_matrix.txt");
  my_out << m2;
  my_out.close();
  std::ifstream is("data/tmp_test_matrix.txt");
  utility::matrix m3(is);
  is.close();
  if (m3!=m2)
    ok=false;
  unlink("data/tmp_test_matrix.txt");

  *error << "\toperator*(double)" << std::endl;
  utility::matrix m4(3,3,1);
  m4 *= 9;
  if (m4!=m) {
    ok=false;
    *error << "error operator*=(double)" << std::endl;
  }

  *error << "\tsub-matrix" << std::endl;
  // Checking that sub-matrices work, i.e. mutation to the view are
  // reflected in the viewed object.
  is.open("data/knni_matrix.data");
  // 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);
  utility::matrix* m5sub=new utility::matrix(m5,3,3,3,3);
  double m5sub_sum=0;
  for (size_t i=0; i<m5sub->rows(); ++i)
    for (size_t j=0; j<m5sub->columns(); ++j) {
      m5sub_sum+=(*m5sub)(i,j);
      (*m5sub)(i,j)=1;
    }
  delete m5sub;
  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);
  if (m5_sum2-m5_sum-9+m5sub_sum>1e-13) {
    ok=false;
    *error << "error sub-matrix test" << std::endl;
  }

  *error << "\tsub-(row)vector" << std::endl;
  // Checking that the row view works, i.e. mutation to the view are
  // reflected in the viewed object.
  utility::vector 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) {
    ok=false;
    *error << "error sub-vector test 1" << std::endl;
  }

  *error << "\tsub-(column)vector" << std::endl;
  // Checking that the column view works, i.e. mutation to the view
  // are reflected in the viewed object.
  utility::vector v5subcolumn(m5,0,false);
  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) {
    ok=false;
    *error << "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) {
    ok=false;
    *error << "error sub-vector test 3" << std::endl;
  }

  *error << "\tsub-vector and vector copying" << std::endl;
  // Checking that a view is not inherited through the copy
  // contructor.
  utility::vector v6(v5subrow);
  v6.set_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) {
    ok=false;
    *error << "error sub-vector test 4" << std::endl;
  }
  // Checking that values in a vector is copied into a viewed matrix.
  v5subrow.set(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) {
    ok=false;
    *error << "error sub-vector test 5" << std::endl;
  }
  // Checking that vector::operator= indeed makes a view to become a
  // "normal" vector.
  v5subrow=utility::vector(23,22);
  double m5_sum6=0;
  for (size_t i=0; i<m5.rows(); ++i)
    for (size_t j=0; j<m5.columns(); ++j)
      m5_sum6+=m5(i,j);
  if (m5_sum6-m5_sum5>1e-13) {
    ok=false;
    *error << "error sub-vector test 6" << std::endl;
  }

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

  *error << "\tsub-matrix of a sub-matrix" << std::endl;
  // Checking that a sub-matrix of a sub-matrix can be created.
  utility::matrix sub(m5,3,3,3,3);
  utility::matrix subsub(sub,2,1,1,2);
  subsub(0,0)=23221121;
  if (&sub(2,1)!=&subsub(0,0) || subsub(0,0)!=m5(5,4) || &subsub(0,0)!=&m5(5,4)){
    ok=false;
    *error << "error sub-matrix test 1" << std::endl;
  }

  *error << "\tmatrix::nan()" << std::endl;
  is.open("data/sorlie_centroids.txt");
  utility::matrix* m_nan = new utility::matrix(is,'\t');
  utility::matrix m_weight;
  m_nan->nan(m_weight);
  is.close();
  if (m_weight(0,0)!=1){
    *error << "error in matrix::nan(): element(0,0) is " << m_weight(0,0)
           << " expected 1." << std::endl;
    ok=false;
  }
  if (m_weight(3,2)!=0){
    *error << "error in matrix::nan(): element(3,2) is " << m_weight(3,2)
           << " expected 0." << std::endl;
    ok=false;
  }
  delete m_nan;

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

  return (ok ? 0 : -1);
}