testSimpleMatrix.cpp

#include <iostream>
#include "DGtal/base/Common.h"
#include "DGtal/kernel/SimpleMatrix.h"

using namespace std;
using namespace DGtal;

// Functions for testing class SimpleMatrix.

bool testSimpleMatrix()
{
  unsigned int nbok = 0;
  unsigned int nb = 0;
  
  trace.beginBlock ( "Testing create ..." );
  
  typedef SimpleMatrix<double,3,4> M34d;
  
  M34d m34d;
  trace.info() << m34d<<std::endl;
  
  m34d.setComponent(1,2, 0.5);
  trace.info() << m34d<<std::endl;
  
  nbok += (m34d(1,2) == 0.5) ? 1 : 0; 
  nb++;
  trace.info() << "(" << nbok << "/" << nb << ") "
               << "true == true" << std::endl;
  
  M34d matrix;
  bool res=true;

  matrix.constant(12.3);
  trace.info() << matrix;
  for(DGtal::Dimension i = 0; i< 3; ++i)
    for(DGtal::Dimension j = 0; j< 4; ++j)
        res = res && (matrix(i,j) == 12.3);
  nbok += res ? 1 : 0;
  nb++;
  trace.info() << "(" << nbok << "/" << nb << ") "
               << "all equals to 12.3" << std::endl;
  

  trace.endBlock();
  
  return nbok == nb;
}

bool testArithm()
{
  unsigned int nbok = 0;
  unsigned int nb = 0;
  
  
  typedef SimpleMatrix<double,3,4> M34d;
  typedef SimpleMatrix<double,4,3> M43d;
  typedef SimpleMatrix<double,3,3> M33d;
  
  M34d m34d, two,four;
  M34d m34dbis, resadd, ressub;
  
  two.constant(2);
  four.constant(4);
  
  for(DGtal::Dimension i = 0; i< 3; ++i)
    for(DGtal::Dimension j = 0; j< 4; ++j)
      {
        m34d.setComponent(i,j,i*j);
        m34dbis.setComponent(i,j,i+j);
        resadd.setComponent(i,j,i*j+i+j);
        ressub.setComponent(i,j,i*j-(double)i-(double)j);
      }
  

  trace.info() << m34d <<std::endl;
  trace.info() << m34dbis<<std::endl;
 
  trace.beginBlock ( "Testing add ..." ); 
  nbok += ((m34d + m34dbis) == resadd) ? 1 : 0; 
  nb++;
  trace.info() << "(" << nbok << "/" << nb << ") "
               << "ok" << std::endl;
  nbok += ((m34dbis + m34d) == resadd) ? 1 : 0; 
  nb++;
  trace.info() << "(" << nbok << "/" << nb << ") "
               << "ok commutative" << std::endl;
  
  M34d other;
  other += m34d;
  nbok += (other == m34d) ? 1 : 0; 
  nb++;
  trace.info() << "(" << nbok << "/" << nb << ") "
               << "ok +=" << std::endl;
  
  trace.endBlock();
  
  trace.beginBlock ( "Testing substraction ..." ); 
  nbok += ((m34d - m34dbis) == ressub) ? 1 : 0; 
  nb++;
  trace.info()<<ressub<<std::endl;
  trace.info()<<m34d - m34dbis<<std::endl;
  
  trace.info() << "(" << nbok << "/" << nb << ") "
               << "ok simple" << std::endl;
  trace.endBlock();
  
  trace.beginBlock ( "Testing scalar product/divide ..." ); 
  nbok += ( (two*2.0) == four) ? 1 : 0; 
  nb++;
  trace.info()<<ressub<<std::endl;
  trace.info() << "(" << nbok << "/" << nb << ") "
               << " [2]*2 == [4]" << std::endl;
 
  nbok += ( two == four/2.0) ? 1 : 0; 
  nb++;
  trace.info()<<ressub<<std::endl;
  trace.info() << "(" << nbok << "/" << nb << ") "
               << " [2]= [4]/2" << std::endl;
  trace.endBlock();
  

  trace.beginBlock ( "Testing transpose ..." ); 
  M43d transp = m34d.transpose();
  nbok += (transp.transpose() == m34d) ? 1 : 0; 
  nb++;
  trace.info() << "(" << nbok << "/" << nb << ") "
               << "ok idem potent" << std::endl;
  trace.endBlock();

  trace.beginBlock ( "Testing product ..." ); 
 
  M43d one;
  M33d eight33;

  one.constant(1);
  eight33.constant(8);
  trace.info() << two * one<<std::endl;
  nbok += (two * one  == eight33) ? 1 : 0; 
  nb++;
  trace.info() << "(" << nbok << "/" << nb << ") "
               << " [2]*[1] = [8]" << std::endl;
  trace.endBlock();


  
  return nbok == nb;

}

bool testColRow()
{
  unsigned int nbok = 0;
  unsigned int nb = 0;
  
  DGtal::SimpleMatrix<double,3,4> mat;
  for(DGtal::Dimension i = 0; i< 3; ++i)
    for(DGtal::Dimension j = 0; j< 4; ++j)
      mat.setComponent(i,j,i+j);
  
  trace.beginBlock("Get Row");
  trace.info() << mat <<std::endl;
  DGtal::SimpleMatrix<double,3,4>::RowVector row;
  row = mat.row(1);
  trace.info() << row << std::endl;
  nbok += (row[1] == 2 ) ? 1 : 0; 
  nb++;
  trace.info() << "(" << nbok << "/" << nb << ") "
               << " row value" << std::endl;
  trace.endBlock();

  trace.beginBlock("Get Col");
  DGtal::SimpleMatrix<double,3,4>::ColumnVector col;
  col = mat.column(1);
  trace.info() << row << std::endl;
  nbok += (col[1] == 2 ) ? 1 : 0; 
  nb++;
  trace.info() << "(" << nbok << "/" << nb << ") "
               << " col value" << std::endl;
  trace.endBlock();

  

  trace.beginBlock("Prod Matrix x Row^t");
  //Row vector is a dim 4 vector
  DGtal::SimpleMatrix<double,3,4>::RowVector r(1,2,3,4);
  DGtal::SimpleMatrix<double,3,4>::ColumnVector c;
  c = mat*r;
  DGtal::SimpleMatrix<double,3,4>::ColumnVector expected(20,30,40);
  
  trace.info() << c << std::endl;
  nbok += (c == expected) ? 1 : 0; 
  nb++;
  trace.info() << "(" << nbok << "/" << nb << ") "
               << " mat*row^t" << std::endl;
  trace.endBlock();

  return nbok == nb; 
}

bool testDetCofactor()
{
  unsigned int nbok = 0;
  unsigned int nb = 0;
  
  typedef DGtal::SimpleMatrix<double,2,2> MAT2;
  MAT2 mat2;
  mat2.setComponent(0,0,1);
  mat2.setComponent(1,1,2);

  trace.beginBlock("det2x2 tests...");
  trace.info() << mat2<<std::endl;
  trace.info() << mat2.determinant() << std::endl;
  nbok += (mat2.determinant() == 2) ? 1 : 0; 
  nb++;
  trace.info() << "(" << nbok << "/" << nb << ") "
               << " 2" << std::endl;
  trace.endBlock();

  typedef DGtal::SimpleMatrix<double,3,3> MAT;
  MAT mat;
  mat.setComponent(0,0,1);
  mat.setComponent(1,1,2);
  mat.setComponent(2,2,4);
  
  trace.beginBlock("det3x3 tests...");
  trace.info() << mat<<std::endl;
  nbok += (mat.determinant() == 8) ? 1 : 0; 
   nb++;
  trace.info() << "(" << nbok << "/" << nb << ") "
               << " 8" << std::endl;
  trace.endBlock();


  typedef DGtal::SimpleMatrix<double,4,4> MAT44;
  MAT44 mat44;
  mat44.setComponent(0,0,1);
  mat44.setComponent(1,1,2);
  mat44.setComponent(2,2,4);
  mat44.setComponent(3,3,4);
  
  trace.beginBlock("det4x4 tests...");
  trace.info() << mat44 <<std::endl;
  trace.info() << mat44.determinant() << std::endl;
  nbok += (mat44.determinant() == 32) ? 1 : 0; 
  nb++;
  trace.info() << "(" << nbok << "/" << nb << ") "
               << " 32" << std::endl;
  trace.endBlock();


  return nbok == nb;
}

bool testM1Matrix()
{
  trace.beginBlock("Mx1 matrix test");
  SimpleMatrix<double, 3,1> mat;
  trace.info() << mat<<std::endl;
  trace.endBlock();
  return true;
}

bool testInverse()
{
  unsigned int nbok = 0;
  unsigned int nb = 0;
 
  trace.beginBlock("Inverse tests 2X2...");
  
  typedef DGtal::SimpleMatrix<double,2,2> MAT2;
  MAT2 mat2;
  mat2.setComponent(0,0,1);
  mat2.setComponent(1,1,2);

  MAT2 Id2;
  Id2.identity();

  trace.info() << mat2<<std::endl;
  trace.info() << mat2.inverse() << std::endl;
  nbok += (( mat2 * mat2.inverse() )== Id2 ) ? 1 : 0; 
  nb++;
  trace.info() << "(" << nbok << "/" << nb << ") "
               << " M*M^-1=Id" << std::endl;
 
  trace.endBlock();

  trace.beginBlock("Inverse tests 6x6 random...");
  
  typedef DGtal::SimpleMatrix<double,6,6> MAT6;
  MAT6 mat;
  
  for(unsigned int i=0; i< 6; i++)
    for(unsigned int j=0; j< 6; j++)
      mat.setComponent(i,j, rand() % 10);
  
  MAT6 Id6;
  Id6.identity();
  
  trace.info() << "M= "<<mat<<std::endl;
  trace.info() << "M^-1=" <<mat.inverse() << std::endl;
  trace.info() << "det(M)= "<<mat.determinant() <<std::endl;
  trace.info() << "M*M^-1= "<<mat.inverse()*mat << std::endl;

  
 
  trace.endBlock();


  return nbok == nb; 
}

// Standard services - public :

int main( int argc, char** argv )
{
  trace.beginBlock ( "Testing class SimpleMatrix" );
  trace.info() << "Args:";
  for ( int i = 0; i < argc; ++i )
    trace.info() << " " << argv[ i ];
  trace.info() << endl;

  bool res = testSimpleMatrix() && testArithm() && testColRow()
    && testDetCofactor() && testM1Matrix()
    && testInverse(); // && ... other tests
  trace.emphase() << ( res ? "Passed." : "Error." ) << endl;
  trace.endBlock();
  return res ? 0 : 1;
}
//                                                                           //