source: ThirdParty/LinearAlgebra/src/unittests/MatrixContentUnitTest.cpp@ 6798a5

/*
 * Project: MoleCuilder
 * Description: creates and alters molecular systems
 * Copyright (C)  2010-2012 University of Bonn. All rights reserved.
 * 
 *
 *   This file is part of MoleCuilder.
 *
 *    MoleCuilder 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.
 *
 *    MoleCuilder 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 MoleCuilder.  If not, see <http://www.gnu.org/licenses/>.
 */

/*
 * MatrixContentUnitTest.cpp
 *
 *  Created on: Jan 8, 2010
 *      Author: heber
 */

// include config.h
#ifdef HAVE_CONFIG_H
#include <config.h>
#endif

using namespace std;

#include <cppunit/CompilerOutputter.h>
#include <cppunit/extensions/TestFactoryRegistry.h>
#include <cppunit/ui/text/TestRunner.h>

#include <cmath>
#include <limits>

// include headers that implement a archive in simple text format
#include <boost/archive/text_oarchive.hpp>
#include <boost/archive/text_iarchive.hpp>

#include "MatrixContentUnitTest.hpp"

#include "MatrixContent.hpp"
#include "VectorContent.hpp"

#ifdef HAVE_TESTRUNNER
#include "UnitTestMain.hpp"
#endif /*HAVE_TESTRUNNER*/

/********************************************** Test classes **************************************/

const std::string matrixA = "   9.962848439806617e-01   5.950413028139907e-01   3.172650036465889e-01\n\
   2.070939109924951e-01   8.365712117773689e-01   3.980567886073226e-01\n\
   1.021600573071414e-01   8.061359922326598e-02   8.493537275043391e-01\n\
   6.955361514014282e-01   6.873206387346102e-02   2.649165458481005e-01";
const std::string vectorS = "1.630290761001098e+00                    7.624807744374841e-01                    6.374093742147465e-01";

// Registers the fixture into the 'registry'
CPPUNIT_TEST_SUITE_REGISTRATION( MatrixContentTest );


void MatrixContentTest::setUp()
{
  m = new MatrixContent(4,3);
};

void MatrixContentTest::tearDown()
{
  delete(m);
};

/** Unit Test for accessing matrix elements.
 *
 */
void MatrixContentTest::AccessTest()
{
  // check whether all elements are initially zero
  for (int i=0;i<4;i++)
    for (int j=0;j<3;j++)
      CPPUNIT_ASSERT_EQUAL( 0., m->at(i,j) );

  // set
  for (int i=0;i<4;i++)
    for (int j=0;j<3;j++)
      m->set(i,j, i*3+j );

  // and check
  double *ptr = NULL;
  for (int i=0;i<4;i++)
    for (int j=0;j<3;j++) {
      CPPUNIT_ASSERT_EQUAL( (double)(i*3+j), m->at(i,j) );
      ptr = m->Pointer(i,j);
      CPPUNIT_ASSERT_EQUAL( (double)(i*3+j), *ptr );
    }

  // assignment
  for (int i=0;i<4;i++)
    for (int j=0;j<3;j++)
      m->set(i,j, i*3+j );
  MatrixContent *dest = new MatrixContent(4,3);
  *dest = *m;
  for (int i=0;i<4;i++)
    for (int j=0;j<3;j++)
      CPPUNIT_ASSERT_EQUAL( dest->at(i,j), m->at(i,j) );
  delete(dest);

  // out of bounds
  //CPPUNIT_ASSERT_EQUAL(0., v->at(5,2) );
  //CPPUNIT_ASSERT_EQUAL(0., v->at(2,17) );
  //CPPUNIT_ASSERT_EQUAL(0., v->at(1024,140040) );
  //CPPUNIT_ASSERT_EQUAL(0., v->at(-1,0) );
  //CPPUNIT_ASSERT_EQUAL(0., v->at(0,-1) );
  //CPPUNIT_ASSERT_EQUAL(0., v->at(-1,-1) );
};

/** Unit Test for initializating matrices.
 *
 */
void MatrixContentTest::InitializationTest()
{
  // set zero
  m->setZero();
  for (int i=0;i<4;i++)
    for (int j=0;j<3;j++)
      CPPUNIT_ASSERT_EQUAL( 0., m->at(i,j) );

  // set all
  m->setValue(1.5);
  for (int i=0;i<4;i++)
    for (int j=0;j<3;j++)
      CPPUNIT_ASSERT_EQUAL( 1.5, m->at(i,j) );

  // set basis
  m->setIdentity();
  for (int i=0;i<4;i++)
    for (int j=0;j<3;j++)
      CPPUNIT_ASSERT_EQUAL( i == j ? 1. : 0. , m->at(i,j) );

  // set from array
  double array[] = { 1., 0., 0.,
                     0., 1., 0.,
                     0., 0., 1.,
                     0., 0., 0. };
  m->setFromDoubleArray(array);
  for (int i=0;i<4;i++)
    for (int j=0;j<3;j++)
      CPPUNIT_ASSERT_EQUAL( i == j ? 1. : 0. , m->at(i,j) );

};

/** Unit Test for copying matrices.
 *
 */
void MatrixContentTest::CopyTest()
{
  // set basis
  MatrixContent *dest = NULL;
  for (int i=0;i<4;i++) {
    m->setValue(i);
    dest = new MatrixContent(m);
    for (int j=0;j<3;j++)
      CPPUNIT_ASSERT_EQUAL( m->at(i,j) , dest->at(i,j) );

    delete(dest);
  }
};

/** Unit Test for exchanging rows and columns.
 *
 */
void MatrixContentTest::ExchangeTest()
{
  // set to 1,1,1,2, ...
  for (int i=0;i<4;i++)
    for (int j=0;j<3;j++)
      m->set(i,j, i+1 );

  // swap such that nothing happens
  CPPUNIT_ASSERT_EQUAL( true, m->SwapColumns(1,2) );
  for (int i=0;i<4;i++)
    for (int j=0;j<3;j++)
      CPPUNIT_ASSERT_EQUAL( (double)i+1., m->at(i,j) );

  // swap two rows
  CPPUNIT_ASSERT_EQUAL( true, m->SwapRows(1,2) );
  for (int i=0;i<4;i++)
    for (int j=0;j<3;j++)
      switch (i) {
        case 0:
          CPPUNIT_ASSERT_EQUAL( 1., m->at(i,j) );
          break;
        case 1:
          CPPUNIT_ASSERT_EQUAL( 3., m->at(i,j) );
          break;
        case 2:
          CPPUNIT_ASSERT_EQUAL( 2., m->at(i,j) );
          break;
        case 3:
          CPPUNIT_ASSERT_EQUAL( 4., m->at(i,j) );
          break;
        default:
          CPPUNIT_ASSERT_EQUAL( -1., m->at(i,j) );
          break;
      }
  // check that op is reversable
  CPPUNIT_ASSERT_EQUAL( true, m->SwapRows(1,2) );
  for (int i=0;i<4;i++)
    for (int j=0;j<3;j++)
      CPPUNIT_ASSERT_EQUAL( (double)i+1., m->at(i,j) );

  // set to 1,2,3,1, ...
  for (int i=0;i<4;i++)
    for (int j=0;j<3;j++)
      m->set(i,j, j+1. );

  // swap such that nothing happens
  CPPUNIT_ASSERT_EQUAL( true, m->SwapRows(0,2) );
  for (int i=0;i<4;i++)
    for (int j=0;j<3;j++)
      CPPUNIT_ASSERT_EQUAL( (double)j+1., m->at(i,j) );

  // swap two columns
  CPPUNIT_ASSERT_EQUAL( true, m->SwapColumns(0,2) );
  for (int i=0;i<4;i++)
    for (int j=0;j<3;j++)
      switch (j) {
        case 0:
          CPPUNIT_ASSERT_EQUAL( 3., m->at(i,j) );
          break;
        case 1:
          CPPUNIT_ASSERT_EQUAL( 2., m->at(i,j) );
          break;
        case 2:
          CPPUNIT_ASSERT_EQUAL( 1., m->at(i,j) );
          break;
        default:
          CPPUNIT_ASSERT_EQUAL( -1., m->at(i,j) );
          break;
      }
  // check that op is reversable
  CPPUNIT_ASSERT_EQUAL( true, m->SwapColumns(0,2) );
  for (int i=0;i<4;i++)
    for (int j=0;j<3;j++)
      CPPUNIT_ASSERT_EQUAL( (double)j+1., m->at(i,j) );


  // set to 1,2,3,4, ...
  MatrixContent *n =  new MatrixContent(3,4);
  for (int i=0;i<4;i++)
    for (int j=0;j<3;j++) {
      m->set(i,j, 3*i+j+1 );
      n->set(j,i, 3*i+j+1 );
    }
  // transpose
  MatrixContent res = (*m).transposed();
  CPPUNIT_ASSERT( *n == res );
  // second transpose
  MatrixContent res2 = res.transposed();
  CPPUNIT_ASSERT( *m == res2 );
  delete n;
};

/** Unit Test for matrix properties.
 *
 */
void MatrixContentTest::PropertiesTest()
{
  // is zero
  m->setZero();
  CPPUNIT_ASSERT_EQUAL( true, m->IsNull() );
  CPPUNIT_ASSERT_EQUAL( false, m->IsPositive() );
  CPPUNIT_ASSERT_EQUAL( false, m->IsNegative() );
  CPPUNIT_ASSERT_EQUAL( true, m->IsNonNegative() );

  // is positive
  m->setValue(0.5);
  CPPUNIT_ASSERT_EQUAL( false, m->IsNull() );
  CPPUNIT_ASSERT_EQUAL( true, m->IsPositive() );
  CPPUNIT_ASSERT_EQUAL( false, m->IsNegative() );
  CPPUNIT_ASSERT_EQUAL( true, m->IsNonNegative() );

  // is negative
  m->setValue(-0.1);
  CPPUNIT_ASSERT_EQUAL( false, m->IsNull() );
  CPPUNIT_ASSERT_EQUAL( false, m->IsPositive() );
  CPPUNIT_ASSERT_EQUAL( true, m->IsNegative() );
  CPPUNIT_ASSERT_EQUAL( false, m->IsNonNegative() );

  // is positive definite
  CPPUNIT_ASSERT_EQUAL( false, m->IsPositiveDefinite() );
};


/** Unit test for reading from and writing matrix to stream
 *
 */
void MatrixContentTest::ReadWriteTest()
{
  // set up matrix
  for (int i=0;i<4;i++)
    for (int j=0;j<3;j++)
      m->set(i,j, i*3+j );

  // write to stream
  std::stringstream matrixstream;
  m->write(matrixstream);

  // parse in dimensions and check
  std::pair<size_t, size_t> matrixdimensions = MatrixContent::preparseMatrixDimensions(matrixstream);
  CPPUNIT_ASSERT_EQUAL( (size_t)4, matrixdimensions.first );
  CPPUNIT_ASSERT_EQUAL( (size_t)3, matrixdimensions.second );
  // parse in matrix and check
  MatrixContent* n = new MatrixContent(4,3, matrixstream);
  CPPUNIT_ASSERT_EQUAL( *m, *n );

  // free matrix
  delete n;
}

/** Unit test for MatrixContent::performSingularValueDecomposition().
 *
 */
void MatrixContentTest::SVDTest()
{
  // setup "A", U,S,V
  std::stringstream Astream(matrixA);
  std::stringstream Sstream(vectorS);
  MatrixContent A(m->getRows(), m->getColumns(), Astream);
  VectorContent S_expected(m->getColumns(), Sstream);
  *m = A;

  // check SVD
  VectorContent S(m->getColumns());
  MatrixContent V(m->getColumns(), m->getColumns());
  A.performSingularValueDecomposition(V,S);
  MatrixContent S_diag(m->getColumns(), m->getColumns());
  for (size_t index = 0; index < m->getColumns(); ++index)
    S_diag.set(index, index, S[index]);
  CPPUNIT_ASSERT_EQUAL(*m,A*S_diag*(V.transpose()));
  for (size_t index = 0; index < m->getColumns(); ++index) {
    //CPPUNIT_ASSERT_EQUAL(S_expected[index], S[index]);
    CPPUNIT_ASSERT(fabs(S_expected[index] - S[index]) < numeric_limits<double>::epsilon()*10.);
  }

  // set up right-hand side
  VectorContent b(4);
  b[0] = 1.;
  b[1] = 0.;
  b[2] = 0.;
  b[3] = 0.;

  // SVD
  VectorContent x_expected(3);
  x_expected[0] = -0.00209169;
  x_expected[1] = -0.325399;
  x_expected[2] = 0.628004;
  const VectorContent x_result(A.solveOverdeterminedLinearEquation(V,S,b));

  // check result
  for (size_t index = 0; index < m->getColumns(); ++index) {
    CPPUNIT_ASSERT(fabs(x_expected[index] - x_result[index]) < 1.e-6);
  }
}

/** Unit test for serialization
 *
 */
void MatrixContentTest::SerializationTest()
{
  m->setZero();
  int count=0;
  for (size_t x = 0; x < 4 ; ++x)
    for (size_t y = 0; y < 3 ; ++y)
      m->at(x,y) = count++;
  // write element to stream
  std::stringstream stream;
  boost::archive::text_oarchive oa(stream);
  oa << m;

  //std::cout << "Contents of archive is " << stream.str() << std::endl;

  // create and open an archive for input
  boost::archive::text_iarchive ia(stream);
  // read class state from archive
  MatrixContent *newm;

  ia >> newm;

  CPPUNIT_ASSERT (*m == *newm);
}