/*
 * Project: MoleCuilder
 * Description: creates and alters molecular systems
 * Copyright (C)  2010 University of Bonn. All rights reserved.
 * Please see the LICENSE file or "Copyright notice" in builder.cpp for details.
 */

/*
 * SubspaceFactorizerUnittest.cpp
 *
 *  Created on: Nov 13, 2010
 *      Author: heber
 */

// 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 <gsl/gsl_vector.h>
#include <boost/shared_ptr.hpp>

#include "Helpers/Assert.hpp"
#include "Helpers/Log.hpp"
#include "Helpers/toString.hpp"
#include "Helpers/Verbose.hpp"
#include "LinearAlgebra/VectorContent.hpp"
#include "LinearAlgebra/MatrixContent.hpp"

#include "SubspaceFactorizerUnittest.hpp"

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

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

void SubspaceFactorizerUnittest::setUp(){
  fourbyfour = new MatrixContent(4,4);
  fourbyfour->setZero();
  for (int i=0; i<4 ; i++) {
    fourbyfour->set(i,i, 2.);
    if (i < (4-1)) {
      fourbyfour->set(i+1,i, 1.);
      fourbyfour->set(i,i+1, 1.);
    }
  }
  transformation = new MatrixContent**[3];

  // 1d subspace
  transformation[0] = new MatrixContent*[4];
  for(size_t i=0; i<4;++i) {
    transformation[0][i] = new MatrixContent(4,4);
    transformation[0][i]->setZero();
    for (size_t j=0; j<1; ++j)
      transformation[0][i]->set(i+j,i+j, 1.);
//    std::cout << i << "th transformation matrix, " << 1 << "d subspace is "
//        << *transformation[0][i] << std::endl;
  }

  // 2d subspace
  transformation[1] = new MatrixContent*[3];
  for(size_t i=0; i<3;++i) {
    transformation[1][i] = new MatrixContent(4,4);
    transformation[1][i]->setZero();
    for (size_t j=0; j<2; ++j)
      transformation[1][i]->set(i+j,i+j, 1.);
//    std::cout << i << "th transformation matrix, " << 2 << "d subspace is "
//        << *transformation[1][i] << std::endl;
  }

  // 3d subspace
  transformation[2] = new MatrixContent*[2];
  for(size_t i=0; i<2;++i) {
    transformation[2][i] = new MatrixContent(4,4);
    transformation[2][i]->setZero();
    for (size_t j=0; j<3; ++j)
      transformation[2][i]->set(i+j,i+j, 1.);
//    std::cout << i << "th transformation matrix, " << 3 << "d subspace is "
//        << *transformation[2][i] << std::endl;
  }
}
void SubspaceFactorizerUnittest::tearDown(){
  // delete test matrix
  delete fourbyfour;

  // delete all transformations
  for(size_t i=0; i<3;++i)
    delete transformation[0][i];
  delete[] transformation[0];
  for(size_t i=0; i<3;++i)
    delete transformation[1][i];
  delete[] transformation[1];
  for(size_t i=0; i<2;++i)
    delete transformation[2][i];
  delete[] transformation[2];
  delete[] transformation;
}

void SubspaceFactorizerUnittest::BlockTest()
{
  MatrixContent temp((*fourbyfour)&(*transformation[0][0]));
  std::cout << "Our matrix is " << *fourbyfour << "." << std::endl;

  std::cout << "Hadamard product of " << *fourbyfour << " with " << *transformation[0][0] << " is: " << std::endl;
  std::cout << temp << std::endl;

  gsl_vector *eigenvalues = temp.transformToEigenbasis();
  VectorContent *eigenvaluesView = new VectorViewContent(gsl_vector_subvector(eigenvalues, 0, eigenvalues->size));
  std::cout << "The resulting eigenbasis is " << temp
      << "\n\t with eigenvalues " << *eigenvaluesView << std::endl;
  delete eigenvaluesView;
  gsl_vector_free(eigenvalues);


  CPPUNIT_ASSERT_EQUAL(0,0);
}

/** For given set of row and column indices, we extract the small block matrix.
 * @param bigmatrix big matrix to extract from
 * @param rowindexset row index set
 * @param columnindexset column index set
 * @return small matrix with dimension equal to the number of indices for row and column.
 */
MatrixContent * getSubspaceMatrix(
    MatrixContent &bigmatrix,
    const IndexSet &rowindexset,
    const IndexSet &columnindexset)
{
  // check whether subsystem is big enough for both index sets
  ASSERT(rowindexset.size() <= bigmatrix.getRows(),
      "embedSubspaceMatrix() - bigmatrix has less rows "+toString(bigmatrix.getRows())
      +" than needed by index set "
      +toString(rowindexset.size())+"!");
  ASSERT(columnindexset.size() <= bigmatrix.getColumns(),
      "embedSubspaceMatrix() - bigmatrix has less columns "+toString(bigmatrix.getColumns())
      +" than needed by index set "
      +toString(columnindexset.size())+"!");
  MatrixContent *subsystem =  new MatrixContent(rowindexset.size(), columnindexset.size());
  size_t localrow = 0; // local row indices for the subsystem
  size_t localcolumn = 0;
  for (IndexSet::const_iterator rowindex = rowindexset.begin();
      rowindex != rowindexset.end();
      ++rowindex) {
    localcolumn = 0;
    for (IndexSet::const_iterator columnindex = columnindexset.begin();
        columnindex != columnindexset.end();
        ++columnindex) {
      ASSERT((*rowindex < bigmatrix.getRows()) && (*columnindex < bigmatrix.getColumns()),
          "current index pair ("
          +toString(*rowindex)+","+toString(*columnindex)
          +") is out of bounds of bigmatrix ("
          +toString(bigmatrix.getRows())+","+toString(bigmatrix.getColumns())
          +")");
      subsystem->at(localrow,localcolumn) = bigmatrix.at(*rowindex, *columnindex);
      localcolumn++;
    }
    localrow++;
  }
  return subsystem;
}

/** For a given set of row and columns indices, we embed a small block matrix into a bigger space.
 *
 * @param bigmatrix big matrix to place submatrix into
 * @param rowindexset
 * @param columnindexset
 * @return
 */
void embedSubspaceMatrix(
    MatrixContent &bigmatrix,
    MatrixContent &subsystem,
    const IndexSet &rowindexset,
    const IndexSet &columnindexset)
{
  // check whether bigmatrix is at least as big as subsystem
  ASSERT(subsystem.getRows() <= bigmatrix.getRows(),
      "embedSubspaceMatrix() - subsystem has more rows "+toString(subsystem.getRows())
      +" than bigmatrix "
      +toString(bigmatrix.getRows())+"!");
  ASSERT(subsystem.getColumns() <= bigmatrix.getColumns(),
      "embedSubspaceMatrix() - subsystem has more columns "+toString(subsystem.getColumns())
      +" than bigmatrix "
      +toString(bigmatrix.getColumns())+"!");
  // check whether subsystem is big enough for both index sets
  ASSERT(rowindexset.size() <= subsystem.getRows(),
      "embedSubspaceMatrix() - subsystem has less rows "+toString(subsystem.getRows())
      +" than needed by index set "
      +toString(rowindexset.size())+"!");
  ASSERT(columnindexset.size() <= subsystem.getColumns(),
      "embedSubspaceMatrix() - subsystem has less columns "+toString(subsystem.getColumns())
      +" than needed by index set "
      +toString(columnindexset.size())+"!");
  size_t localrow = 0; // local row indices for the subsystem
  size_t localcolumn = 0;
  for (IndexSet::const_iterator rowindex = rowindexset.begin();
      rowindex != rowindexset.end();
      ++rowindex) {
    localcolumn = 0;
    for (IndexSet::const_iterator columnindex = columnindexset.begin();
        columnindex != columnindexset.end();
        ++columnindex) {
      bigmatrix.at(*rowindex, *columnindex) = subsystem.at(localrow,localcolumn);
      localcolumn++;
    }
    localrow++;
  }
}

/** Operator for output to std:: ostream operator of an IndexSet.
 * @param ost output stream
 * @param indexset index set to output
 * @return ost output stream
 */
std::ostream & operator<<(std::ostream &ost, const IndexSet &indexset)
{
  ost << "{ ";
  for (IndexSet::const_iterator iter = indexset.begin();
      iter != indexset.end();
      ++iter)
    ost << *iter << " ";
  ost << "}";
  return ost;
}

void SubspaceFactorizerUnittest::EigenvectorTest()
{
  VectorArray CurrentEigenvectors;
  VectorList *ParallelEigenvectorList = new VectorList[4];
  IndexSet AllIndices;

  // create the total index set
  for (size_t i=0;i<4;++i)
    AllIndices.insert(i);

  // create all consecutive index subsets for dim 1 to 3
  IndexMap Dimension_to_Indexset;
  for (size_t dim = 0; dim<3;++dim) {
    for (size_t i=0;i<4-dim;++i) {
      IndexSet *indexset = new IndexSet;
      for (size_t j=0;j<=dim;++j) {
        //std::cout << "Putting " << i+j << " into " << i << "th map " << dim << std::endl;
        CPPUNIT_ASSERT_MESSAGE("index "+toString(i+j)+" already present in "+toString(dim)+"-dim "+toString(i)+"th indexset.", indexset->count(i+j) == 0);
        indexset->insert(i+j);
      }
      Dimension_to_Indexset.insert( make_pair(dim, boost::shared_ptr<IndexSet>(indexset)) );
      // no need to free indexset, is stored in shared_ptr and
      // will get released when Dimension_to_Indexset is destroyed
    }
  }

  // set to first guess, i.e. the unit vectors of R^4
  for (IndexSet::const_iterator index = AllIndices.begin();
      index != AllIndices.end();
      ++index) {
    boost::shared_ptr<VectorContent> EV(new VectorContent(4));
    EV->setZero();
    EV->at(*index) = 1.;
    CurrentEigenvectors.push_back(EV);
  }

  // for every dimension
  for (size_t dim = 0; dim<3;++dim) {
    // for every index set of this dimension
    Log() << Verbose(0) << std::endl << std::endl;
    Log() << Verbose(0) << "Current dimension is " << dim << std::endl;
    std::pair<IndexMap::const_iterator,IndexMap::const_iterator> Bounds = Dimension_to_Indexset.equal_range(dim);
    for (IndexMap::const_iterator IndexsetIter = Bounds.first;
        IndexsetIter != Bounds.second;
        ++IndexsetIter) {
      // show the index set
      Log() << Verbose(0) << std::endl;
      Log() << Verbose(1) << "Current index set is " << *(IndexsetIter->second) << std::endl;

      // create transformation matrices from these
      MatrixContent *subsystem = getSubspaceMatrix(*fourbyfour, *(IndexsetIter->second), *(IndexsetIter->second));
      Log() << Verbose(2) << "Subsystem matrix is " << *subsystem << std::endl;

      // solve _small_ systems for eigenvalues
      VectorContent *Eigenvalues = new VectorContent(subsystem->transformToEigenbasis());
      Log() << Verbose(2) << "Eigenvector matrix is " << *subsystem << std::endl;
      Log() << Verbose(2) << "Eigenvalues are " << *Eigenvalues << std::endl;
      delete Eigenvalues;

      // blow up eigenvectors to 4dim column vector again
      MatrixContent *Eigenvectors = new MatrixContent(4,4);
      Eigenvectors->setZero();
      embedSubspaceMatrix(*Eigenvectors, *subsystem, *(IndexsetIter->second), *(IndexsetIter->second));
      Log() << Verbose(1) << "4x4 Eigenvector matrix is " << *Eigenvectors << std::endl;

      // we don't need the subsystem anymore
      delete subsystem;

      // copy the index set, we look for one new eigenvector for each old one
      IndexSet CorrespondenceList((*IndexsetIter->second));

      // go through all eigenvectors in this subspace
      for (IndexSet::const_iterator iter = (*IndexsetIter->second).begin();
          iter != (*IndexsetIter->second).end();
          ++iter) {
        // we rather copy the column vector, as the containing matrix is destroyed lateron
        VectorContent *CurrentEigenvector = new VectorContent(Eigenvectors->getColumnVector(*iter));
        Log() << Verbose(1) << "Current Eigenvector is " << *CurrentEigenvector << std::endl;

        // recognize eigenvectors parallel to existing ones
        // (for now settle with the one we are most parallel to)
        size_t mostparallel_index = 4;
        double mostparallel_scalarproduct = 0.;
        for (IndexSet::const_iterator indexiter = CorrespondenceList.begin();
            indexiter != CorrespondenceList.end();
            ++indexiter) {
          Log() << Verbose(2) << "Comparing to old " << *indexiter << "th eigenvector " << *(CurrentEigenvectors[*indexiter]) << std::endl;
          const double scalarproduct = (*(CurrentEigenvectors[*indexiter])) * (*CurrentEigenvector);
          Log() << Verbose(2) << "SKP is " << scalarproduct << std::endl;
          if (fabs(scalarproduct) > mostparallel_scalarproduct) {
            mostparallel_scalarproduct = fabs(scalarproduct);
            mostparallel_index = *indexiter;
          }
        }
        if (mostparallel_index != 4) {
          // put into std::list for later use
          Log() << Verbose(1) << "Pushing " << *CurrentEigenvector << " into parallel list [" << mostparallel_index << "]" << std::endl;
          ParallelEigenvectorList[mostparallel_index].push_back(boost::shared_ptr<VectorContent>(CurrentEigenvector));
          CorrespondenceList.erase(mostparallel_index);
        }
        // no need to delete CurrentEigenvector as is taken care of by shared_ptr
      }
      // free eigenvectors
      delete Eigenvectors;
    }
  }

  // print list of parallel eigenvectors
  for (IndexSet::const_iterator index = AllIndices.begin();
      index != AllIndices.end();
      ++index) {
    Log() << Verbose(0) << "Similar to " << *index << "th current eigenvector " << *(CurrentEigenvectors[*index]) << " are:" << std::endl;
    for (VectorList::const_iterator iter = ParallelEigenvectorList[*index].begin();
        iter != ParallelEigenvectorList[*index].end();
        ++iter) {
      Log() << Verbose(0) << **iter << std::endl;
    }
    Log() << Verbose(0) << std::endl;
  }

  // create new CurrentEigenvectors from averaging parallel ones.
  for (IndexSet::const_iterator index = AllIndices.begin();
      index != AllIndices.end();
      ++index) {
    for (VectorList::const_iterator iter = ParallelEigenvectorList[*index].begin();
        iter != ParallelEigenvectorList[*index].end();
        ++iter) {
      *CurrentEigenvectors[*index] += **iter;
    }
    *CurrentEigenvectors[*index] *= 1./(double)(ParallelEigenvectorList[*index].size()+1);
  }

  // show new ones
  Log() << Verbose(0) << "Resulting new eigenvectors are:" << std::endl;
  for (IndexSet::const_iterator index = AllIndices.begin();
      index != AllIndices.end();
      ++index) {
    Log() << Verbose(0) << *CurrentEigenvectors[*index] << std::endl;
  }


  delete[] ParallelEigenvectorList;

  CPPUNIT_ASSERT_EQUAL(0,0);
}