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

using namespace std;

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

#include "gslmatrixsymmetricunittest.hpp"

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

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

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


void GSLMatrixSymmetricTest::setUp()
{
  m = new GSLMatrix(3,3);
};

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

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

  // set
  for (int i=0;i<3;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<3;i++)
    for (int j=0;j<3;j++) {
      CPPUNIT_ASSERT_EQUAL( (double)(i*3+j), m->Get(i,j) );
      ptr = m->Pointer(i,j);
      CPPUNIT_ASSERT_EQUAL( (double)(i*3+j), *ptr );
    }

  // assignment
  for (int i=0;i<3;i++)
    for (int j=0;j<3;j++)
      m->Set(i,j, i*3+j );
  GSLMatrix *dest = new GSLMatrix(3,3);
  *dest = *m;
  for (int i=0;i<3;i++)
    for (int j=0;j<3;j++)
      CPPUNIT_ASSERT_EQUAL( dest->Get(i,j), m->Get(i,j) );
  delete(dest);

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

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

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

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

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

};

/** Unit Test for copying matrices.
 *
 */
void GSLMatrixSymmetricTest::CopyTest()
{
  // set basis
  GSLMatrix *dest = NULL;
  for (int i=0;i<3;i++) {
    m->SetAll(i);
    dest = new GSLMatrix(m);
    for (int j=0;j<3;j++)
      CPPUNIT_ASSERT_EQUAL( m->Get(i,j) , dest->Get(i,j) );

    delete(dest);
  }
};

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

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

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

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

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

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


  // set to 1,2,3,4, ...
  for (int i=0;i<3;i++)
    for (int j=0;j<3;j++)
      m->Set(i,j, 3*i+j+1 );
  // transpose
  m->Transpose();
  for (int i=0;i<3;i++)
    for (int j=0;j<3;j++)
      CPPUNIT_ASSERT_EQUAL( (double)(3*i+j+1), m->Get(j,i) );
  // second transpose
  m->Transpose();
  for (int i=0;i<3;i++)
    for (int j=0;j<3;j++)
      CPPUNIT_ASSERT_EQUAL( (double)(3*i+j+1), m->Get(i,j) );
};

/** Unit Test for matrix properties.
 *
 */
void GSLMatrixSymmetricTest::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->SetAll(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->SetAll(-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
  double array[] = { 1., 0., 0.,
                     0., 1., 1.,
                     0., 0., 1. };
  m->SetFromDoubleArray(array);
  CPPUNIT_ASSERT_EQUAL( true, m->IsPositiveDefinite() );

  //determinant
  m->SetIdentity();
  CPPUNIT_ASSERT_EQUAL( 1., m->Determinant() );

  m->SetZero();
  CPPUNIT_ASSERT_EQUAL( 0., m->Determinant() );

  m->Set( 0, 0, 1.);
  m->Set( 1, 1, 1.);
  m->Set( 2, 1, 1.);
  CPPUNIT_ASSERT_EQUAL( 0., m->Determinant() );
  
  double array2[] = { 2., 0., 1.,
                     -3., 1., 1.,
                     1., 5.5, 1. };
  m->SetFromDoubleArray(array2);
  CPPUNIT_ASSERT_EQUAL( -26.5, m->Determinant() );
};