source: src/LinearAlgebra/Subspace.cpp@ 286af5f

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

/*
 * Subspace.cpp
 *
 *  Created on: Nov 22, 2010
 *      Author: heber
 */

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

#include "Helpers/MemDebug.hpp"

#include <boost/shared_ptr.hpp>
#include <boost/foreach.hpp>

#include "Helpers/Assert.hpp"
#include "Helpers/Log.hpp"
#include "Helpers/toString.hpp"
#include "Helpers/Verbose.hpp"

#include "Eigenspace.hpp"
#include "Subspace.hpp"


/** Constructor for class Subspace.
 *
 * @param _s indices that make up this subspace
 * @param _FullSpace reference to the full eigenspace
 */
Subspace::Subspace(indexset & _s, Eigenspace &_FullSpace) :
  Eigenspace(_s),
  ProjectToSubspace(_FullSpace.getIndices().size(), _s.size()),
  ProjectFromSubspace(_s.size(), _FullSpace.getIndices().size()),
  FullSpace(_FullSpace),
  ZeroVector(_FullSpace.getIndices().size())
{
  // TODO: away with both of this when calculateEigenSubspace() works
  // create projection matrices
  createProjectionMatrices();

  // create eigenspace matrices
  projectFullSpaceMatrixToSubspace();
}


/** Destructor for class Subspace.
 *
 */
Subspace::~Subspace()
{}


/** Diagonalizes the subspace matrix.
 *
 * @param fullmatrix full matrix to project into subspace and solve
 */
void Subspace::calculateEigenSubspace()
{
  // project down
  createProjectionMatrices();
  // remove subsubspace directions
  correctEigenvectorsFromSubIndices();
  // solve
  projectFullSpaceMatrixToSubspace();
  EigenvectorMatrix = EigenspaceMatrix;
  VectorContent *EigenvalueVector = new VectorContent(EigenvectorMatrix.transformToEigenbasis());
  Eigenvalues.clear();
  for(size_t i = 0; i< EigenvalueVector->getDimension(); ++i) {
    Eigenvalues.push_back( EigenvalueVector->at(i) );
  }
  delete EigenvalueVector;
  Log() << Verbose(2) << "Eigenvector matrix is " << EigenvectorMatrix << std::endl;
  Log() << Verbose(2) << "Eigenvalues are " << Eigenvalues << std::endl;
  // create mapping
  createLocalMapping();
}


/** Projects a given full matrix down into the subspace of this class.
 *
 * @param bigmatrix full matrix to project into subspace
 */
void Subspace::getSubspacematrixFromBigmatrix(const MatrixContent & bigmatrix)
{
  // check whether subsystem is big enough for both index sets
  ASSERT(Indices.size() <= bigmatrix.getRows(),
      "embedEigenspaceMatrix() - bigmatrix has less rows "+toString(bigmatrix.getRows())
      +" than needed by index set "
      +toString(Indices.size())+"!");
  ASSERT(Indices.size() <= bigmatrix.getColumns(),
      "embedEigenspaceMatrix() - bigmatrix has less columns "+toString(bigmatrix.getColumns())
      +" than needed by index set "
      +toString(Indices.size())+"!");

  // construct small matrix
  MatrixContent *subsystem =  new MatrixContent(Indices.size(), Indices.size());
  size_t localrow = 0; // local row indices for the subsystem
  size_t localcolumn = 0;
  BOOST_FOREACH( size_t rowindex, Indices) {
    localcolumn = 0;
    BOOST_FOREACH( size_t columnindex, Indices) {
      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) = (*Eigenvectors[rowindex]) * (bigmatrix * (*Eigenvectors[columnindex]));
      localcolumn++;
    }
    localrow++;
  }
}


void Subspace::correctEigenvectorsFromSubIndices()
{
}


/** Creates the inverse of LocalToGlobal.
 * Mapping is one-to-one and onto, hence it's simply turning around the map.
 */
void Subspace::invertLocalToGlobalMapping()
{
  GlobalToLocal.clear();
  for (mapping::const_iterator iter = LocalToGlobal.begin();
      iter != LocalToGlobal.end();
      ++iter) {
    GlobalToLocal.insert( make_pair(iter->second, iter->first) );
  }
}


/** Project calculated Eigenvectors into full space.
 *
 * @return set of projected eigenvectors.
 */
const Eigenspace::eigenvectorset & Subspace::getEigenvectorsInFullSpace()
{
  // check whether bigmatrix is at least as big as EigenspaceMatrix
  ASSERT(Eigenvectors.size() > 0,
      "embedEigenspaceMatrix() - no Eigenvectors given!");
  ASSERT(EigenspaceMatrix.getRows() <= Eigenvectors[0]->getDimension(),
      "embedEigenspaceMatrix() - EigenspaceMatrix has more rows "
      +toString(EigenspaceMatrix.getRows())+" than eigenvectors "
      +toString(Eigenvectors[0]->getDimension())+"!");
  ASSERT(EigenspaceMatrix.getColumns() <= Eigenvectors.size(),
      "embedEigenspaceMatrix() - EigenspaceMatrix has more columns "
      +toString(EigenspaceMatrix.getColumns())+" than eigenvectors "
      +toString(Eigenvectors.size())+"!");
  // check whether EigenspaceMatrix is big enough for both index sets
  ASSERT(EigenspaceMatrix.getColumns() == EigenspaceMatrix.getRows(),
      "embedEigenspaceMatrix() - EigenspaceMatrix is not square "
      +toString(EigenspaceMatrix.getRows())+" than needed by index set "
      +toString(EigenspaceMatrix.getColumns())+"!");
  ASSERT(Indices.size() == EigenspaceMatrix.getColumns(),
      "embedEigenspaceMatrix() - EigenspaceMatrix has not the same number of columns "
      +toString(EigenspaceMatrix.getColumns())+" compared to the index set "
      +toString(Indices.size())+"!");

  // construct intermediate matrix
  MatrixContent *intermediatematrix = new MatrixContent(Eigenvectors[0]->getDimension(), Indices.size());
  size_t localcolumn = 0;
  BOOST_FOREACH(size_t columnindex, Indices) {
    // create two vectors from each row and copy assign them
    boost::shared_ptr<VectorContent> srceigenvector(Eigenvectors[columnindex]);
    boost::shared_ptr<VectorContent> desteigenvector(intermediatematrix->getColumnVector(localcolumn));
    *desteigenvector = *srceigenvector;
    localcolumn++;
  }
  // matrix product with eigenbasis EigenspaceMatrix matrix
  *intermediatematrix *= EigenspaceMatrix;

  // and place at right columns into bigmatrix
  MatrixContent *bigmatrix = new MatrixContent(Eigenvectors[0]->getDimension(), Eigenvectors.size());
  bigmatrix->setZero();
  localcolumn = 0;
  BOOST_FOREACH(size_t columnindex, Indices) {
    // create two vectors from each row and copy assign them
    boost::shared_ptr<VectorContent> srceigenvector(intermediatematrix->getColumnVector(localcolumn));
    boost::shared_ptr<VectorContent> desteigenvector(bigmatrix->getColumnVector(columnindex));
    *desteigenvector = *srceigenvector;
    localcolumn++;
  }

  return Eigenvectors;
}


/** Remove a subset of indices from the SubIndices.
 *
 * @param _s subset to remove
 * @return true - removed, false - not found
 */
bool Subspace::removeSubset(boost::shared_ptr<Subspace> &_s)
{
  if (SubIndices.count(_s) != 0) {
    SubIndices.erase(_s);
    return true;
  } else {
    return false;
  }
}

/** Add a subspace set to SubIndices.
 *
 * @param _s subset to add
 * @return true - not present before, false - has already been present
 */
bool Subspace::addSubset(boost::shared_ptr<Subspace> & _s)
{
  if (SubIndices.count(_s) != 0)
    return false;
  else {
    SubIndices.insert(_s);
    return true;
  }
}


/** Simply a getter for Eigenvectors, as they are stored as subspace eigenvectors.
 *
 * @return set of eigenvectors in subspace
 */
const MatrixContent & Subspace::getEigenvectorMatrixInFullSpace()
{
  // check whether bigmatrix is at least as big as EigenspaceMatrix
  ASSERT(EigenspaceMatrix.getRows() <= Eigenvectors[0]->getDimension(),
      "embedEigenspaceMatrix() - EigenspaceMatrix has more rows "
      +toString(EigenspaceMatrix.getRows())+" than eigenvectors "
      +toString(Eigenvectors[0]->getDimension())+"!");
  ASSERT(EigenspaceMatrix.getColumns() <= Eigenvectors.size(),
      "embedEigenspaceMatrix() - EigenspaceMatrix has more columns "
      +toString(EigenspaceMatrix.getColumns())+" than eigenvectors "
      +toString(Eigenvectors.size())+"!");
  // check whether EigenspaceMatrix is big enough for both index sets
  ASSERT(EigenspaceMatrix.getColumns() == EigenspaceMatrix.getRows(),
      "embedEigenspaceMatrix() - EigenspaceMatrix is not square "
      +toString(EigenspaceMatrix.getRows())+" than needed by index set "
      +toString(EigenspaceMatrix.getColumns())+"!");
  ASSERT(Indices.size() == EigenspaceMatrix.getColumns(),
      "embedEigenspaceMatrix() - EigenspaceMatrix has not the same number of columns "
      +toString(EigenspaceMatrix.getColumns())+" compared to the index set "
      +toString(Indices.size())+"!");

  // construct intermediate matrix
  MatrixContent *intermediatematrix = new MatrixContent(Eigenvectors[0]->getDimension(), Indices.size());
  size_t localcolumn = 0;
  BOOST_FOREACH(size_t columnindex, Indices) {
    // create two vectors from each row and copy assign them
    boost::shared_ptr<VectorContent> srceigenvector(Eigenvectors[columnindex]);
    boost::shared_ptr<VectorContent> desteigenvector(intermediatematrix->getColumnVector(localcolumn));
    *desteigenvector = *srceigenvector;
    localcolumn++;
  }
  // matrix product with eigenbasis EigenspaceMatrix matrix
  *intermediatematrix *= EigenspaceMatrix;

  // and place at right columns into bigmatrix
  MatrixContent *bigmatrix = new MatrixContent(Eigenvectors[0]->getDimension(), Eigenvectors.size());
  bigmatrix->setZero();
  localcolumn = 0;
  BOOST_FOREACH(size_t columnindex, Indices) {
    // create two vectors from each row and copy assign them
    boost::shared_ptr<VectorContent> srceigenvector(intermediatematrix->getColumnVector(localcolumn));
    boost::shared_ptr<VectorContent> desteigenvector(bigmatrix->getColumnVector(columnindex));
    *desteigenvector = *srceigenvector;
    localcolumn++;
  }

  return *bigmatrix;
}


/** Creates the projection matrix from Subspace::FullSpace::EigenvectorMatrix.
 * We simply copy the respective eigenvectors from Subspace::FullSpace into
 * Subspace::Eigenvectors.
 */
void Subspace::createProjectionMatrices()
{
  // first ProjectToSubspace

  // generate column vectors
  std::vector< boost::shared_ptr<VectorContent> > ColumnVectors;
  for (size_t i = 0; i < Indices.size(); ++i) {
    boost::shared_ptr<VectorContent> p(ProjectToSubspace.getColumnVector(i));
    ColumnVectors.push_back( p );
  }

  // then copy from real eigenvectors
  size_t localindex = 0;
  BOOST_FOREACH(size_t iter, Indices) {
    *ColumnVectors[localindex] = FullSpace.getEigenvector(iter);
    localindex++;
  }
  Log() << Verbose(2) << "ProjectionToSubspace matrix is " << ProjectToSubspace << std::endl;

  // then ProjectFromSubspace is just transposed
  ProjectFromSubspace = ((const MatrixContent)ProjectToSubspace).transpose();
  Log() << Verbose(2) << "ProjectFromSubspace matrix is " << ProjectFromSubspace << std::endl;
}


/** Creates the subspace matrix by projecting down the FullSpace::EigenspaceMatrix.
 *  We just perform \f$M = P^t \cdot A \cdot P\f$, when P are the projection matrix,
 *  A the full matrix and M the subspace matrix.
 */
void Subspace::projectFullSpaceMatrixToSubspace()
{
  // construct small matrix
  EigenspaceMatrix = ProjectFromSubspace * FullSpace.getEigenspaceMatrix() * ProjectToSubspace;
  Log() << Verbose(2) << "EigenspaceMatrix matrix is " << EigenspaceMatrix << std::endl;
}


/** We associate local Eigenvectors with ones from FullSpace by a paralellity
 * criterion.
 */
void Subspace::createLocalMapping()
{
  // first LocalToGlobal
  LocalToGlobal.clear();
  IndexSet CorrespondenceList(Indices); // is to ensure one-to-one mapping

  for (size_t localindex = 0; localindex< Indices.size(); ++localindex) {
    boost::shared_ptr<VectorContent> &CurrentEigenvector = Eigenvectors[localindex];
    Log() << Verbose(1) << "Current Eigenvector is " << *CurrentEigenvector << std::endl;

    // (for now settle with the one we are most parallel to)
    size_t mostparallel_index = FullSpace.getIndices().size();
    double mostparallel_scalarproduct = 0.;
    BOOST_FOREACH( size_t indexiter, CorrespondenceList) {
      Log() << Verbose(2) << "Comparing to old " << indexiter << "th eigenvector " << FullSpace.getEigenvector(indexiter) << std::endl;
      const double scalarproduct = (FullSpace.getEigenvector(indexiter)) * ( ProjectToSubspace * (*CurrentEigenvector));
      Log() << Verbose(2) << "SKP is " << scalarproduct << std::endl;
      if (fabs(scalarproduct) > fabs(mostparallel_scalarproduct)) {
        mostparallel_scalarproduct = scalarproduct;
        mostparallel_index = indexiter;
      }
    }
    // and make the assignment if we found one
    if (mostparallel_index != FullSpace.getIndices().size()) {
      // put into std::list for later use
      // invert if pointing in negative direction
      if (mostparallel_scalarproduct < 0) {
        *CurrentEigenvector *= -1.;
        Log() << Verbose(1) << "Associating (inverted) " << *CurrentEigenvector << " with " << mostparallel_index << "th's fullspace eigenvector." << std::endl;
      } else {
        Log() << Verbose(1) << "Associating " << *CurrentEigenvector << " with " << mostparallel_index << "th's fullspace eigenvector." << std::endl;
      }
      ASSERT( LocalToGlobal.count(localindex) == 0,
          "Subspace::createLocalMapping() - localindex "+toString(localindex)+" already assigned to "
          +toString(LocalToGlobal[localindex])+" !=? "+toString(mostparallel_index)+".");
      LocalToGlobal.insert( make_pair(localindex, mostparallel_index) );
      CorrespondenceList.erase(mostparallel_index);
    }
  }

  // then invert mapping
  GlobalToLocal.clear();
  BOOST_FOREACH(mapping::value_type iter, LocalToGlobal) {
    ASSERT(GlobalToLocal.count(iter.second) == 0,
        "Subspace::createLocalMapping() - LocalToGlobal is not bijective, key "
        +toString(iter.second)+" present more than once!");
    GlobalToLocal.insert( make_pair(iter.second, iter.first) );
  }
}


/** Get the local eigenvector that is most parallel to the \a i'th full one.
 * We just the internal mapping Subspace::GlobalToLocal.
 * @param i index of global eigenvector
 * @return local eigenvector, most parallel
 */
const VectorContent Subspace::getEigenvectorParallelToFullOne(size_t i)
{
  if (GlobalToLocal.count(i) == 0) {
    return ZeroVector;
  } else {
    return ProjectToSubspace*(*Eigenvectors[GlobalToLocal[i]]);
  }
}


/** Get the local eigenvalue of the eigenvector that is most parallel to the
 * \a i'th full one.
 * We just the internal mapping Subspace::GlobalToLocal.
 * @param i index of global eigenvector
 * @return eigenvalue of local eigenvector, most parallel
 */
const double Subspace::getEigenvalueOfEigenvectorParallelToFullOne(size_t i)
{
  if (GlobalToLocal.count(i) == 0) {
    return 0.;
  } else {
    return Eigenvalues[GlobalToLocal[i]];
  }
}