source: src/Fragmentation/MatrixContainer.cpp@ dcbb5d

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

/*
 * MatrixContainer.cpp
 *
 *  Created on: Sep 15, 2011
 *      Author: heber
 */

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

#include "CodePatterns/MemDebug.hpp"

#include <cstring>
#include <fstream>
#include <iomanip>

#include "CodePatterns/Log.hpp"
#include "KeySetsContainer.hpp"

#include "Fragmentation/helpers.hpp"
#include "Helpers/defs.hpp"
#include "Helpers/helpers.hpp"

#include "MatrixContainer.hpp"

/** Constructor of MatrixContainer class.
 */
MatrixContainer::MatrixContainer()
{
  Header.resize(1);
  RowCounter.resize(1);
  ColumnCounter.resize(1);
  ColumnCounter[0] = -1;
  MatrixCounter = 0;
};

/** Destructor of MatrixContainer class.
 */
MatrixContainer::~MatrixContainer()
{}

/** Either copies index matrix from another MatrixContainer or initialises with trivial mapping if NULL.
 * This either copies the index matrix or just maps 1 to 1, 2 to 2 and so on for all fragments.
 * \param *Matrix pointer to other MatrixContainer
 * \return true - copy/initialisation sucessful, false - dimension false for copying
 */
bool MatrixContainer::InitialiseIndices(class MatrixContainer *_container)
{
  if (_container == NULL) {
    LOG(3, "INFO: Initialising indices with trivial mapping.");
    Indices.resize(MatrixCounter + 1);
    for(int i=MatrixCounter+1;i--;) {
      Indices[i].resize(RowCounter[i]);
      for(int j=RowCounter[i];j--;)
        Indices[i][j] = j;
    }
  } else {
    std::stringstream output;
    if (MatrixCounter != _container->MatrixCounter)
      return false;
    Indices.resize(MatrixCounter + 1);
    for(int i=MatrixCounter+1;i--;) {
      if (RowCounter[i] != _container->RowCounter[i])
        return false;
      Indices[i].resize(_container->RowCounter[i]);
      for(int j=_container->RowCounter[i];j--;) {
        Indices[i][j] = _container->Indices[i][j];
        output << Indices[i][j] << "\t";
      }
    }
    LOG(3, "INFO: Initialising indices from other MatrixContainer: " << output.str());
  }
  return true;
};

/** Parsing a number of matrices.
 *    -# open the matrix file
 *    -# skip some lines (\a skiplines)
 *    -# scan header lines for number of columns
 *    -# scan lines for number of rows
 *    -# allocate a temporary matrix
 *    -# loop over found column and row counts and parse in each entry
 *    -# use MatrixContainer::AddMatrix() to add the parsed matrix to the internal
 * \param &input input stream
 * \param skiplines number of inital lines to skip
 * \param skiplines number of inital columns to skip
 * \param MatrixNr index number in Matrix array to parse into
 * \return parsing successful
 */
bool MatrixContainer::ParseMatrix(std::istream &input, int skiplines, int skipcolumns, size_t MatrixNr)
{
  stringstream line;
  string token;
  char filename[1023];

  if (input.fail()) {
    ELOG(1, endl << "MatrixContainer::ParseMatrix: Unable to parse istream.");
    //performCriticalExit();
    return false;
  }

  // parse header
  std::string header;
  char dummy[1024];
  for (int m=skiplines+1;m--;)
    input.getline(dummy, 1023);
  line.str(dummy);
  for(int k=skipcolumns;k--;)
    line >> header;
  LOG(3, "INFO: Header of Matrix " << MatrixNr << " :" << line.str());

  // scan header for number of columns
  size_t ColumnCounter = 0;
  while ( getline(line,token, '\t') ) {
    if (token.length() > 0)
      ColumnCounter++;
  }
  LOG(3, "INFO: "+line.str());

  // scan rest for number of rows/lines
  size_t RowCounter = -1;
  while (!input.eof()) {
    input.getline(filename, 1023);
    LOG(3, "INFO: Comparing: " << strncmp(filename,"MeanForce",9));
    RowCounter++; // then line was not last MeanForce
    if (strncmp(filename,"MeanForce", 9) == 0) {
      break;
    }
  }

  // allocate temporary matrix
  MatrixArray temp_matrix;
  temp_matrix.resize(RowCounter);
  for(MatrixArray::iterator iter = temp_matrix.begin(); iter != temp_matrix.end(); ++iter)
    (*iter).resize(ColumnCounter);

  // parse in each entry for this matrix
  input.clear();
  input.seekg(ios::beg);
  for (int m=skiplines+1;m--;)
    input.getline(dummy, 1023);    // skip header
  line.str(dummy);
  LOG(3, "INFO: Header: " << line.str());
  for(int k=skipcolumns;k--;)  // skip columns in header too
    line >> filename;
  header = line.str();
  for(size_t j=0;j<RowCounter;j++) {
    input.getline(filename, 1023);
    std::stringstream lines(filename);
    std::stringstream output;
    output << "INFO: Matrix at level " << j << ":";// << filename << endl;
    for(int k=skipcolumns;k--;)
      lines >> filename;
    for(size_t k=0;(k<ColumnCounter) && (!lines.eof());k++) {
      lines >> temp_matrix[j][k];
      output << " " << std::setprecision(2) << temp_matrix[j][k] << endl;
    }
    LOG(3, output.str());
  }

  // finally add the matrix
  return AddMatrix(header, temp_matrix, MatrixNr);
}

/** Adds a matrix at position \a MatrixNr to MatrixContainer::Matrix.
 *
 * @param header header to add for this matrix
 * @param matrix to add
 * @param MatrixNr position in MatrixContainer::Matrix.
 * @return true - insertion ok, false - invalid matrix
 */
bool MatrixContainer::AddMatrix(const std::string &header, const MatrixArray &matrix, size_t MatrixNr)
{
  // make some pre-checks
  if (header.size() == 0)
    ELOG(2, "The given header of the matrix to add is empty.");
  if (matrix.size() == 0) {
    ELOG(1, "RowCounter[" << MatrixNr << "]: " << RowCounter[MatrixNr] << " from input stream.");
    return false;
  }
  if (matrix[0].size() == 0) {
    ELOG(1, "ColumnCounter[" << MatrixNr << "]: " << ColumnCounter[MatrixNr] << " from ostream.");
    return false;
  }

  // add header
  if (Header.size() <= MatrixNr)
    Header.resize(MatrixNr+1);
  Header[MatrixNr] = header;

  // row count
  if (RowCounter.size() <= MatrixNr)
    RowCounter.resize(MatrixNr+1);
  RowCounter[MatrixNr] = matrix.size();
  LOG(4, "INFO: RowCounter[" << MatrixNr << "]: " << RowCounter[MatrixNr] << " from input stream.");

  // column count
  if (ColumnCounter.size() <= MatrixNr)
    ColumnCounter.resize(MatrixNr+1);
  ColumnCounter[MatrixNr] = matrix[0].size();
  LOG(4, "INFO: ColumnCounter[" << MatrixNr << "]: " << ColumnCounter[MatrixNr] << ".");

  // allocate matrix ...
  if (Matrix.size() <= MatrixNr)
    Matrix.resize(MatrixNr+1);
  MatrixCounter = Matrix.size()-1;
  Matrix[MatrixNr].resize(RowCounter[MatrixNr] + 1);
  for(int j=0;j<=RowCounter[MatrixNr];++j)
    Matrix[MatrixNr][j].resize(ColumnCounter[MatrixNr]+1);

  // .. and copy values
  for(int j=0;j<RowCounter[MatrixNr];++j)
    for(int k=0;k<ColumnCounter[MatrixNr];++k)
      Matrix[MatrixNr][j][k] = matrix[j][k];
  // reset last column
  for(int j=0;j<RowCounter[MatrixNr];++j)
    Matrix[MatrixNr][j][ ColumnCounter[MatrixNr] ] = 0.;
  // reset last row
  for(int k=0;k<=ColumnCounter[MatrixNr];++k)
    Matrix[MatrixNr][ RowCounter[MatrixNr] ][ k ] = 0.;

  return true;
}

/** Parsing a number of matrices.
 * -# First, count the number of matrices by counting lines in KEYSETFILE
 * -# Then,
 *    -# construct the fragment number
 *    -# open the matrix file
 *    -# skip some lines (\a skiplines)
 *    -# scan header lines for number of columns
 *    -# scan lines for number of rows
 *    -# allocate matrix
 *    -# loop over found column and row counts and parse in each entry
 * -# Finally, allocate one additional matrix (\a MatrixCounter) containing combined or temporary values
 * \param *name directory with files
 * \param *prefix prefix of each matrix file
 * \param *suffix suffix of each matrix file
 * \param skiplines number of inital lines to skip
 * \param skiplines number of inital columns to skip
 * \return parsing successful
 */
bool MatrixContainer::ParseFragmentMatrix(const std::string name, const std::string prefix, std::string suffix, int skiplines, int skipcolumns)
{
  char filename[1023];
  ifstream input;
  char *FragmentNumber = NULL;
  stringstream file;
  string token;

  // count the number of matrices
  MatrixCounter = -1; // we count one too much
  file << name << FRAGMENTPREFIX << KEYSETFILE;
  input.open(file.str().c_str(), ios::in);
  if (input.bad()) {
    ELOG(1, "MatrixContainer::ParseFragmentMatrix: Unable to open " << file.str() << ", is the directory correct?");
    return false;
  }
  while (!input.eof()) {
    input.getline(filename, 1023);
    stringstream zeile(filename);
    MatrixCounter++;
  }
  input.close();
  LOG(2, "INFO: Determined " << MatrixCounter << " fragments.");

  LOG(1, "STATUS: Parsing through each fragment and retrieving " << prefix << suffix << ".");
  Header.clear();
  Matrix.clear();
  RowCounter.clear();
  ColumnCounter.clear();
  Header.resize(MatrixCounter + 1); // one more each for the total molecule
  Matrix.resize(MatrixCounter + 1); // one more each for the total molecule
  RowCounter.resize(MatrixCounter + 1);
  ColumnCounter.resize(MatrixCounter + 1);
  for(int i=0; i < MatrixCounter;i++) {
    // open matrix file
    FragmentNumber = FixedDigitNumber(MatrixCounter, i);
    file.str(" ");
    file << name << FRAGMENTPREFIX << FragmentNumber << prefix << suffix;
    std::ifstream input(file.str().c_str());
    LOG(2, "INFO: Opening " << file.str() << " ... ");
    if (!ParseMatrix(input, skiplines, skipcolumns, i)) {
      input.close();
      return false;
    }
    input.close();
    delete[](FragmentNumber);
  }
  return true;
};

/** Allocates and resets the memory for a number \a MCounter of matrices.
 * \param **GivenHeader Header line for each matrix
 * \param MCounter number of matrices
 * \param *RCounter number of rows for each matrix
 * \param *CCounter number of columns for each matrix
 * \return Allocation successful
 */
bool MatrixContainer::AllocateMatrix(StringVector GivenHeader, int MCounter, IntVector RCounter, IntVector CCounter)
{
  MatrixCounter = MCounter;
  Header.resize(MatrixCounter + 1);
  Matrix.resize(MatrixCounter + 1); // one more each for the total molecule
  RowCounter.resize(MatrixCounter + 1);
  ColumnCounter.resize(MatrixCounter + 1);
  for(int i=MatrixCounter+1;i--;) {
    Header[i] = GivenHeader[i];
    RowCounter[i] = RCounter[i];
    ColumnCounter[i] = CCounter[i];
    if ((int)Matrix[i].size() <= RowCounter[i] + 2)
      Matrix[i].resize(RowCounter[i] + 1);
    for(int j=0;j<=RowCounter[i];j++)
      if ((int)Matrix[i][j].size() <= ColumnCounter[i]+1)
        Matrix[i][j].resize(ColumnCounter[i]);
    // allocation with 0 is guaranted by STL
  }
  return true;
};

/** Resets all values in MatrixContainer::Matrix.
 * \return true if successful
 */
bool MatrixContainer::ResetMatrix()
{
  for(int i=MatrixCounter+1;i--;)
    for(int j=RowCounter[i]+1;j--;)
      for(int k=ColumnCounter[i];k--;)
        Matrix[i][j][k] = 0.;
   return true;
};

/** Scans all elements of MatrixContainer::Matrix for greatest absolute value.
 * \return greatest value of MatrixContainer::Matrix
 */
double MatrixContainer::FindMaxValue()
{
  double max = Matrix[0][0][0];
  for(int i=MatrixCounter+1;i--;)
    for(int j=RowCounter[i]+1;j--;)
      for(int k=ColumnCounter[i];k--;)
        if (fabs(Matrix[i][j][k]) > max)
          max = fabs(Matrix[i][j][k]);
  if (fabs(max) < MYEPSILON)
    max += MYEPSILON;
 return max;
};

/** Scans all elements of MatrixContainer::Matrix for smallest absolute value.
 * \return smallest value of MatrixContainer::Matrix
 */
double MatrixContainer::FindMinValue()
{
  double min = Matrix[0][0][0];
  for(int i=MatrixCounter+1;i--;)
    for(int j=RowCounter[i]+1;j--;)
      for(int k=ColumnCounter[i];k--;)
        if (fabs(Matrix[i][j][k]) < min)
          min = fabs(Matrix[i][j][k]);
  if (fabs(min) < MYEPSILON)
    min += MYEPSILON;
  return min;
};

/** Sets all values in the last of MatrixContainer::Matrix to \a value.
 * \param value reset value
 * \param skipcolumns skip initial columns
 * \return true if successful
 */
bool MatrixContainer::SetLastMatrix(double value, int skipcolumns)
{
  for(int j=RowCounter[MatrixCounter]+1;j--;)
    for(int k=skipcolumns;k<ColumnCounter[MatrixCounter];k++)
      Matrix[MatrixCounter][j][k] = value;
   return true;
};

/** Sets all values in the last of MatrixContainer::Matrix to \a value.
 * \param **values matrix with each value (must have at least same dimensions!)
 * \param skipcolumns skip initial columns
 * \return true if successful
 */
bool MatrixContainer::SetLastMatrix(const MatrixArray &values, int skipcolumns)
{
  for(int j=RowCounter[MatrixCounter]+1;j--;)
    for(int k=skipcolumns;k<ColumnCounter[MatrixCounter];k++)
      Matrix[MatrixCounter][j][k] = values[j][k];
   return true;
};

/** Sums the entries with each factor and put into last element of \a ***Matrix.
 * Sums over "E"-terms to create the "F"-terms
 * \param Matrix MatrixContainer with matrices (LevelCounter by *ColumnCounter) with all the energies.
 * \param KeySets KeySetContainer with bond Order and association mapping of each fragment to an order
 * \param Order bond order
 * \return true if summing was successful
 */
bool MatrixContainer::SumSubManyBodyTerms(class MatrixContainer &MatrixValues, class KeySetsContainer &KeySets, int Order)
{
  // go through each order
  for (int CurrentFragment=0;CurrentFragment<KeySets.FragmentsPerOrder[Order];CurrentFragment++) {
    //LOG(0, "Current Fragment is " << CurrentFragment << "/" << KeySets.OrderSet[Order][CurrentFragment] << ".");
    // then go per order through each suborder and pick together all the terms that contain this fragment
    for(int SubOrder=0;SubOrder<=Order;SubOrder++) { // go through all suborders up to the desired order
      for (int j=0;j<KeySets.FragmentsPerOrder[SubOrder];j++) { // go through all possible fragments of size suborder
        if (KeySets.Contains(KeySets.OrderSet[Order][CurrentFragment], KeySets.OrderSet[SubOrder][j])) {
          //LOG(0, "Current other fragment is " << j << "/" << KeySets.OrderSet[SubOrder][j] << ".");
          // if the fragment's indices are all in the current fragment
          for(int k=0;k<RowCounter[ KeySets.OrderSet[SubOrder][j] ];k++) { // go through all atoms in this fragment
            int m = MatrixValues.Indices[ KeySets.OrderSet[SubOrder][j] ][k];
            //LOG(0, "Current index is " << k << "/" << m << ".");
            if (m != -1) { // if it's not an added hydrogen
              for (int l=0;l<RowCounter[ KeySets.OrderSet[Order][CurrentFragment] ];l++) { // look for the corresponding index in the current fragment
                //LOG(0, "Comparing " << m << " with " << MatrixValues.Indices[ KeySets.OrderSet[Order][CurrentFragment] ][l] << ".");
                if (m == MatrixValues.Indices[ KeySets.OrderSet[Order][CurrentFragment] ][l]) {
                  m = l;
                  break;
                }
              }
              //LOG(0, "Corresponding index in CurrentFragment is " << m << ".");
              if (m > RowCounter[ KeySets.OrderSet[Order][CurrentFragment] ]) {
                ELOG(0, "In fragment No. " << KeySets.OrderSet[Order][CurrentFragment]   << " current force index " << m << " is greater than " << RowCounter[ KeySets.OrderSet[Order][CurrentFragment] ] << "!");
                performCriticalExit();
                return false;
              }
              if (Order == SubOrder) { // equal order is always copy from Energies
                for(int l=ColumnCounter[ KeySets.OrderSet[SubOrder][j] ];l--;) // then adds/subtract each column
                  Matrix[ KeySets.OrderSet[Order][CurrentFragment] ][m][l] += MatrixValues.Matrix[ KeySets.OrderSet[SubOrder][j] ][k][l];
              } else {
                for(int l=ColumnCounter[ KeySets.OrderSet[SubOrder][j] ];l--;)
                  Matrix[ KeySets.OrderSet[Order][CurrentFragment] ][m][l] -= Matrix[ KeySets.OrderSet[SubOrder][j] ][k][l];
              }
            }
            //if ((ColumnCounter[ KeySets.OrderSet[SubOrder][j] ]>1) && (RowCounter[0]-1 >= 1))
             //LOG(0, "Fragments[ KeySets.OrderSet[" << Order << "][" << CurrentFragment << "]=" << KeySets.OrderSet[Order][CurrentFragment] << " ][" << RowCounter[0]-1 << "][" << 1 << "] = " <<  Matrix[ KeySets.OrderSet[Order][CurrentFragment] ][RowCounter[0]-1][1]);
          }
        } else {
          //LOG(0, "Fragment " << KeySets.OrderSet[SubOrder][j] << " is not contained in fragment " << KeySets.OrderSet[Order][CurrentFragment] << ".");
        }
      }
    }
   //LOG(0, "Final Fragments[ KeySets.OrderSet[" << Order << "][" << CurrentFragment << "]=" << KeySets.OrderSet[Order][CurrentFragment] << " ][" << KeySets.AtomCounter[0]-1 << "][" << 1 << "] = " <<  Matrix[ KeySets.OrderSet[Order][CurrentFragment] ][KeySets.AtomCounter[0]-1][1]);
  }

  return true;
};

/** Writes the summed total fragment terms \f$F_{ij}\f$ to file.
 * \param *name inputdir
 * \param *prefix prefix before \a EnergySuffix
 * \return file was written
 */
bool MatrixContainer::WriteTotalFragments(const std::string name, const std::string prefix)
{
  ofstream output;
  char *FragmentNumber = NULL;

  LOG(1, "STATUS: Writing fragment files.");
  for(int i=0;i<MatrixCounter;i++) {
    stringstream line;
    FragmentNumber = FixedDigitNumber(MatrixCounter, i);
    line << name << FRAGMENTPREFIX << FragmentNumber << "/" << prefix;
    delete[](FragmentNumber);
    output.open(line.str().c_str(), ios::out);
    if (output == NULL) {
      ELOG(0, "MatrixContainer::WriteTotalFragments: Unable to open output energy file " << line.str() << "!");
      performCriticalExit();
      return false;
    }
    output << Header[i] << endl;
    for(int j=0;j<RowCounter[i];j++) {
      for(int k=0;k<ColumnCounter[i];k++)
        output << scientific << Matrix[i][j][k] << "\t";
      output << endl;
    }
    output.close();
  }
  return true;
};

/** Writes the summed total values in the last matrix to file.
 * \param *name inputdir
 * \param *prefix prefix
 * \param *suffix suffix
 * \return file was written
 */
bool MatrixContainer::WriteLastMatrix(const std::string name, const std::string prefix, const std::string suffix)
{
  ofstream output;
  stringstream line;

  LOG(1, "STATUS: Writing matrix values of " << suffix << ".");
  line << name << prefix << suffix;
  output.open(line.str().c_str(), ios::out);
  if (output == NULL) {
    ELOG(0, "MatrixContainer::WriteLastMatrix: Unable to open output matrix file " << line.str() << "!");
    performCriticalExit();
    return false;
  }
  output << Header[MatrixCounter] << endl;
  for(int j=0;j<RowCounter[MatrixCounter];j++) {
    for(int k=0;k<ColumnCounter[MatrixCounter];k++)
      output << scientific << Matrix[MatrixCounter][j][k] << "\t";
    output << endl;
  }
  output.close();
  return true;
};

/** Comparison operator for class MatrixContainer.
 *
 * @param other instance to compare to
 * @return true - both instances are the same in each member variable.
 */
bool MatrixContainer::operator==(const MatrixContainer &other) const
{
  // compare matrices
  if (Matrix != other.Matrix)
    return false;
  // compare Indices
  if (Indices != other.Indices)
    return false;
  // compare Headers
  if (Header != other.Header)
    return false;
  // compare MatrixCounter
  if (MatrixCounter != other.MatrixCounter)
    return false;
  // compare RowCounter
  if (RowCounter != other.RowCounter)
    return false;
  // compare ColumnCounter
  if (ColumnCounter != other.ColumnCounter)
    return false;
  return true;
}

std::ostream & operator << (std::ostream &ost, const MatrixContainer &m)
{
  for (int i=0;i<=m.MatrixCounter;++i) {
    ost << "Matrix " << i << " in MatrixContainer:" << std::endl;
    for (int j=0;j<=m.RowCounter[i];++j) {
      for (int k=0;k<m.ColumnCounter[i];++k) {
        ost << m.Matrix[i][j][k] << " ";
      }
      ost << std::endl;
    }
  }
  ost << std::endl;
  return ost;
}