/*
 * 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/>.
 */

/*
 * RealSpaceMatrixUnitTest.cpp
 *
 *  Created on: Jul 7, 2010
 *      Author: crueger
 */

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

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

#include <cmath>
#include <limits>

#include "RealSpaceMatrixUnitTest.hpp"

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

#include "CodePatterns/Assert.hpp"

#include "defs.hpp"
#include "Exceptions.hpp"
#include "MatrixContent.hpp"
#include "RealSpaceMatrix.hpp"
#include "Vector.hpp"

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

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

void RealSpaceMatrixTest::setUp(){
  // failing asserts should be thrown
  ASSERT_DO(Assert::Throw);

  zero = new RealSpaceMatrix();
  for(int i =NDIM;i--;) {
    for(int j =NDIM;j--;) {
      zero->at(i,j)=0.;
    }
  }
  one = new RealSpaceMatrix();
  for(int i =NDIM;i--;){
    one->at(i,i)=1.;
  }
  full=new RealSpaceMatrix();
  for(int i=NDIM;i--;){
    for(int j=NDIM;j--;){
      full->at(i,j)=1.;
    }
  }
  diagonal = new RealSpaceMatrix();
  for(int i=NDIM;i--;){
    diagonal->at(i,i)=i+1.;
  }
  perm1 = new RealSpaceMatrix();
  perm1->column(0) = unitVec[0];
  perm1->column(1) = unitVec[2];
  perm1->column(2) = unitVec[1];


  perm2 = new RealSpaceMatrix();
  perm2->column(0) = unitVec[1];
  perm2->column(1) = unitVec[0];
  perm2->column(2) = unitVec[2];

  perm3 = new RealSpaceMatrix();
  perm3->column(0) = unitVec[1];
  perm3->column(1) = unitVec[2];
  perm3->column(2) = unitVec[0];

  perm4 = new RealSpaceMatrix();
  perm4->column(0) = unitVec[2];
  perm4->column(1) = unitVec[1];
  perm4->column(2) = unitVec[0];

  perm5 = new RealSpaceMatrix();
  perm5->column(0) = unitVec[2];
  perm5->column(1) = unitVec[0];
  perm5->column(2) = unitVec[1];

}
void RealSpaceMatrixTest::tearDown(){
  delete zero;
  delete one;
  delete full;
  delete diagonal;
  delete perm1;
  delete perm2;
  delete perm3;
  delete perm4;
  delete perm5;
}

void RealSpaceMatrixTest::AccessTest(){
  RealSpaceMatrix mat;
  for(int i=NDIM;i--;){
    for(int j=NDIM;j--;){
      CPPUNIT_ASSERT_EQUAL(mat.at(i,j),0.);
    }
  }
  int k=1;
  for(int i=NDIM;i--;){
    for(int j=NDIM;j--;){
      mat.at(i,j)=k++;
    }
  }
  k=1;
  for(int i=NDIM;i--;){
    for(int j=NDIM;j--;){
      CPPUNIT_ASSERT_EQUAL(mat.at(i,j),(double)k);
      ++k;
    }
  }
}

void RealSpaceMatrixTest::VectorTest(){
  RealSpaceMatrix mat;
  for(int i=NDIM;i--;){
    CPPUNIT_ASSERT_EQUAL(mat.row(i),zeroVec);
    CPPUNIT_ASSERT_EQUAL(mat.column(i),zeroVec);
  }
  CPPUNIT_ASSERT_EQUAL(mat.diagonal(),zeroVec);

  mat.setIdentity();
  CPPUNIT_ASSERT_EQUAL(mat.row(0),unitVec[0]);
  CPPUNIT_ASSERT_EQUAL(mat.row(1),unitVec[1]);
  CPPUNIT_ASSERT_EQUAL(mat.row(2),unitVec[2]);
  CPPUNIT_ASSERT_EQUAL(mat.column(0),unitVec[0]);
  CPPUNIT_ASSERT_EQUAL(mat.column(1),unitVec[1]);
  CPPUNIT_ASSERT_EQUAL(mat.column(2),unitVec[2]);

  Vector t1=Vector(1.,1.,1.);
  Vector t2=Vector(2.,2.,2.);
  Vector t3=Vector(3.,3.,3.);
  Vector t4=Vector(1.,2.,3.);

  mat.row(0)=t1;
  mat.row(1)=t2;
  mat.row(2)=t3;
  CPPUNIT_ASSERT_EQUAL(mat.row(0),t1);
  CPPUNIT_ASSERT_EQUAL(mat.row(1),t2);
  CPPUNIT_ASSERT_EQUAL(mat.row(2),t3);
  CPPUNIT_ASSERT_EQUAL(mat.column(0),t4);
  CPPUNIT_ASSERT_EQUAL(mat.column(1),t4);
  CPPUNIT_ASSERT_EQUAL(mat.column(2),t4);
  CPPUNIT_ASSERT_EQUAL(mat.diagonal(),t4);
  for(int i=NDIM;i--;){
    for(int j=NDIM;j--;){
      CPPUNIT_ASSERT_EQUAL(mat.at(i,j),i+1.);
    }
  }

  mat.column(0)=t1;
  mat.column(1)=t2;
  mat.column(2)=t3;
  CPPUNIT_ASSERT_EQUAL(mat.column(0),t1);
  CPPUNIT_ASSERT_EQUAL(mat.column(1),t2);
  CPPUNIT_ASSERT_EQUAL(mat.column(2),t3);
  CPPUNIT_ASSERT_EQUAL(mat.row(0),t4);
  CPPUNIT_ASSERT_EQUAL(mat.row(1),t4);
  CPPUNIT_ASSERT_EQUAL(mat.row(2),t4);
  CPPUNIT_ASSERT_EQUAL(mat.diagonal(),t4);
  for(int i=NDIM;i--;){
    for(int j=NDIM;j--;){
      CPPUNIT_ASSERT_EQUAL(mat.at(i,j),j+1.);
    }
  }
}

void RealSpaceMatrixTest::TransposeTest(){
  RealSpaceMatrix res;

  // transpose of unit is unit
  res.setIdentity();
  res.transpose();
  CPPUNIT_ASSERT_EQUAL(res,*one);

  // transpose of transpose is same matrix
  res.setZero();
  res.set(2,2, 1.);
  CPPUNIT_ASSERT_EQUAL(res.transpose().transpose(),res);
}

void RealSpaceMatrixTest::OperationTest(){
  RealSpaceMatrix res;

  res =(*zero) *(*zero);
  //std::cout << *zero << " times " << *zero << " is " << res << std::endl;
  CPPUNIT_ASSERT_EQUAL(res,*zero);
  res =(*zero) *(*one);
  CPPUNIT_ASSERT_EQUAL(res,*zero);
  res =(*zero) *(*full);
  CPPUNIT_ASSERT_EQUAL(res,*zero);
  res =(*zero) *(*diagonal);
  CPPUNIT_ASSERT_EQUAL(res,*zero);
  res =(*zero) *(*perm1);
  CPPUNIT_ASSERT_EQUAL(res,*zero);
  res =(*zero) *(*perm2);
  CPPUNIT_ASSERT_EQUAL(res,*zero);
  res =(*zero) *(*perm3);
  CPPUNIT_ASSERT_EQUAL(res,*zero);
  res =(*zero) *(*perm4);
  CPPUNIT_ASSERT_EQUAL(res,*zero);
  res =(*zero) *(*perm5);
  CPPUNIT_ASSERT_EQUAL(res,*zero);

  res =(*one)*(*one);
  CPPUNIT_ASSERT_EQUAL(res,*one);
  res =(*one)*(*full);
  CPPUNIT_ASSERT_EQUAL(res,*full);
  res =(*one)*(*diagonal);
  CPPUNIT_ASSERT_EQUAL(res,*diagonal);
  res =(*one)*(*perm1);
  CPPUNIT_ASSERT_EQUAL(res,*perm1);
  res =(*one)*(*perm2);
  CPPUNIT_ASSERT_EQUAL(res,*perm2);
  res =(*one)*(*perm3);
  CPPUNIT_ASSERT_EQUAL(res,*perm3);
  res =(*one)*(*perm4);
  CPPUNIT_ASSERT_EQUAL(res,*perm4);
  res =(*one)*(*perm5);
  CPPUNIT_ASSERT_EQUAL(res,*perm5);

  res = (*full)*(*perm1);
  CPPUNIT_ASSERT_EQUAL(res,*full);
  res = (*full)*(*perm2);
  CPPUNIT_ASSERT_EQUAL(res,*full);
  res = (*full)*(*perm3);
  CPPUNIT_ASSERT_EQUAL(res,*full);
  res = (*full)*(*perm4);
  CPPUNIT_ASSERT_EQUAL(res,*full);
  res = (*full)*(*perm5);
  CPPUNIT_ASSERT_EQUAL(res,*full);

  res = (*diagonal)*(*perm1);
  CPPUNIT_ASSERT_EQUAL(res.column(0),unitVec[0]);
  CPPUNIT_ASSERT_EQUAL(res.column(1),3*unitVec[2]);
  CPPUNIT_ASSERT_EQUAL(res.column(2),2*unitVec[1]);
  res = (*diagonal)*(*perm2);
  CPPUNIT_ASSERT_EQUAL(res.column(0),2*unitVec[1]);
  CPPUNIT_ASSERT_EQUAL(res.column(1),unitVec[0]);
  CPPUNIT_ASSERT_EQUAL(res.column(2),3*unitVec[2]);
  res = (*diagonal)*(*perm3);
  CPPUNIT_ASSERT_EQUAL(res.column(0),2*unitVec[1]);
  CPPUNIT_ASSERT_EQUAL(res.column(1),3*unitVec[2]);
  CPPUNIT_ASSERT_EQUAL(res.column(2),unitVec[0]);
  res = (*diagonal)*(*perm4);
  CPPUNIT_ASSERT_EQUAL(res.column(0),3*unitVec[2]);
  CPPUNIT_ASSERT_EQUAL(res.column(1),2*unitVec[1]);
  CPPUNIT_ASSERT_EQUAL(res.column(2),unitVec[0]);
  res = (*diagonal)*(*perm5);
  CPPUNIT_ASSERT_EQUAL(res.column(0),3*unitVec[2]);
  CPPUNIT_ASSERT_EQUAL(res.column(1),unitVec[0]);
  CPPUNIT_ASSERT_EQUAL(res.column(2),2*unitVec[1]);
}

void RealSpaceMatrixTest::RotationTest(){
  RealSpaceMatrix res;
  RealSpaceMatrix inverse;
  Vector fixture;

  // zero rotation angles yields unity matrix
  res.setRotation(0,0,0);
  CPPUNIT_ASSERT_EQUAL(*one, res);

  // arbitrary rotation matrix has det = 1
  res.setRotation(M_PI/3.,1.,M_PI/7.);
  CPPUNIT_ASSERT(fabs(res.determinant() - 1.) <= LINALG_MYEPSILON());

  // inverse is rotation matrix with negative angles
  res.setRotation(M_PI/3.,0.,0.);
  inverse.setRotation(-M_PI/3.,0.,0.);
  CPPUNIT_ASSERT_EQUAL(*one, res * inverse);

  // ... or transposed
  res.setRotation(M_PI/3.,0.,0.);
  CPPUNIT_ASSERT_EQUAL(inverse, res.transposed());

  // check arbitrary axis
  res.setRotation(unitVec[0], M_PI/3.);
  inverse.setRotation(M_PI/3.,0.,0.);
  CPPUNIT_ASSERT_EQUAL( res, inverse);
  res.setRotation(unitVec[1], M_PI/3.);
  inverse.setRotation(0., M_PI/3.,0.);
  CPPUNIT_ASSERT_EQUAL( res, inverse);
  res.setRotation(unitVec[2], M_PI/3.);
  inverse.setRotation(0.,0.,M_PI/3.);
  CPPUNIT_ASSERT_EQUAL( res, inverse);

  // check axis not normalized throws
#ifndef NDEBUG
  CPPUNIT_ASSERT_THROW( res.setRotation(unitVec[0]+unitVec[1], M_PI/3.),Assert::AssertionFailure);
#endif

  // check arbitrary axis and det=1
  fixture = (unitVec[0]+unitVec[1]).getNormalized();
  res.setRotation(fixture, M_PI/3.);
  CPPUNIT_ASSERT(fabs(res.determinant() - 1.) <= LINALG_MYEPSILON());
  fixture = (unitVec[1]+unitVec[2]).getNormalized();
  res.setRotation(fixture, M_PI/3.);
  CPPUNIT_ASSERT(fabs(res.determinant() - 1.) <= LINALG_MYEPSILON());
  fixture = (unitVec[1]+unitVec[2]).getNormalized();
  res.setRotation(fixture, M_PI/3.);
  CPPUNIT_ASSERT(fabs(res.determinant() - 1.) <= LINALG_MYEPSILON());
}

void RealSpaceMatrixTest::InvertTest(){
  CPPUNIT_ASSERT_THROW(zero->invert(),NotInvertibleException);
  CPPUNIT_ASSERT_THROW(full->invert(),NotInvertibleException);

  RealSpaceMatrix res;
  res = (*one)*one->invert();
  CPPUNIT_ASSERT_EQUAL(res,*one);
  res = (*diagonal)*diagonal->invert();
  CPPUNIT_ASSERT_EQUAL(res,*one);
  res = (*perm1)*perm1->invert();
  CPPUNIT_ASSERT_EQUAL(res,*one);
  res = (*perm2)*perm2->invert();
  CPPUNIT_ASSERT_EQUAL(res,*one);
  res = (*perm3)*perm3->invert();
  CPPUNIT_ASSERT_EQUAL(res,*one);
  res = (*perm4)*perm4->invert();
  CPPUNIT_ASSERT_EQUAL(res,*one);
  res = (*perm5)*perm5->invert();
  CPPUNIT_ASSERT_EQUAL(res,*one);
}


void RealSpaceMatrixTest::DeterminantTest(){
  CPPUNIT_ASSERT_EQUAL(zero->determinant(),0.);
  CPPUNIT_ASSERT_EQUAL(one->determinant(),1.);
  CPPUNIT_ASSERT_EQUAL(diagonal->determinant(),6.);
  CPPUNIT_ASSERT_EQUAL(full->determinant(),0.);
  CPPUNIT_ASSERT_EQUAL(perm1->determinant(),-1.);
  CPPUNIT_ASSERT_EQUAL(perm2->determinant(),-1.);
  CPPUNIT_ASSERT_EQUAL(perm3->determinant(),1.);
  CPPUNIT_ASSERT_EQUAL(perm4->determinant(),-1.);
  CPPUNIT_ASSERT_EQUAL(perm5->determinant(),1.);
}

void RealSpaceMatrixTest::VecMultTest(){
  CPPUNIT_ASSERT_EQUAL((*zero)*unitVec[0],zeroVec);
  CPPUNIT_ASSERT_EQUAL((*zero)*unitVec[1],zeroVec);
  CPPUNIT_ASSERT_EQUAL((*zero)*unitVec[2],zeroVec);
  CPPUNIT_ASSERT_EQUAL((*zero)*zeroVec,zeroVec);

  CPPUNIT_ASSERT_EQUAL((*one)*unitVec[0],unitVec[0]);
  CPPUNIT_ASSERT_EQUAL((*one)*unitVec[1],unitVec[1]);
  CPPUNIT_ASSERT_EQUAL((*one)*unitVec[2],unitVec[2]);
  CPPUNIT_ASSERT_EQUAL((*one)*zeroVec,zeroVec);

  CPPUNIT_ASSERT_EQUAL((*diagonal)*unitVec[0],unitVec[0]);
  CPPUNIT_ASSERT_EQUAL((*diagonal)*unitVec[1],2*unitVec[1]);
  CPPUNIT_ASSERT_EQUAL((*diagonal)*unitVec[2],3*unitVec[2]);
  CPPUNIT_ASSERT_EQUAL((*diagonal)*zeroVec,zeroVec);

  CPPUNIT_ASSERT_EQUAL((*perm1)*unitVec[0],unitVec[0]);
  CPPUNIT_ASSERT_EQUAL((*perm1)*unitVec[1],unitVec[2]);
  CPPUNIT_ASSERT_EQUAL((*perm1)*unitVec[2],unitVec[1]);
  CPPUNIT_ASSERT_EQUAL((*perm1)*zeroVec,zeroVec);

  CPPUNIT_ASSERT_EQUAL((*perm2)*unitVec[0],unitVec[1]);
  CPPUNIT_ASSERT_EQUAL((*perm2)*unitVec[1],unitVec[0]);
  CPPUNIT_ASSERT_EQUAL((*perm2)*unitVec[2],unitVec[2]);
  CPPUNIT_ASSERT_EQUAL((*perm2)*zeroVec,zeroVec);

  CPPUNIT_ASSERT_EQUAL((*perm3)*unitVec[0],unitVec[1]);
  CPPUNIT_ASSERT_EQUAL((*perm3)*unitVec[1],unitVec[2]);
  CPPUNIT_ASSERT_EQUAL((*perm3)*unitVec[2],unitVec[0]);
  CPPUNIT_ASSERT_EQUAL((*perm3)*zeroVec,zeroVec);

  CPPUNIT_ASSERT_EQUAL((*perm4)*unitVec[0],unitVec[2]);
  CPPUNIT_ASSERT_EQUAL((*perm4)*unitVec[1],unitVec[1]);
  CPPUNIT_ASSERT_EQUAL((*perm4)*unitVec[2],unitVec[0]);
  CPPUNIT_ASSERT_EQUAL((*perm4)*zeroVec,zeroVec);

  CPPUNIT_ASSERT_EQUAL((*perm5)*unitVec[0],unitVec[2]);
  CPPUNIT_ASSERT_EQUAL((*perm5)*unitVec[1],unitVec[0]);
  CPPUNIT_ASSERT_EQUAL((*perm5)*unitVec[2],unitVec[1]);
  CPPUNIT_ASSERT_EQUAL((*perm5)*zeroVec,zeroVec);

  Vector t = Vector(1.,2.,3.);
  CPPUNIT_ASSERT_EQUAL((*perm1)*t,Vector(1,3,2));
  CPPUNIT_ASSERT_EQUAL((*perm2)*t,Vector(2,1,3));
  CPPUNIT_ASSERT_EQUAL((*perm3)*t,Vector(3,1,2));
  CPPUNIT_ASSERT_EQUAL((*perm4)*t,Vector(3,2,1));
  CPPUNIT_ASSERT_EQUAL((*perm5)*t,Vector(2,3,1));
}

void RealSpaceMatrixTest::SerializationTest()
{
  // write element to stream
  std::stringstream stream;
  boost::archive::text_oarchive oa(stream);
  oa << diagonal;

  //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
  RealSpaceMatrix *newm;

  ia >> newm;

  CPPUNIT_ASSERT (*diagonal == *newm);
}