source: src/Fragmentation/MatrixContainer.cpp@ a9b86d

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

/*
 * 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() :
  Indices(NULL)
{
  Header = new char*[1];
  Matrix = new double**[1]; // one more each for the total molecule
  RowCounter = new int[1];
  ColumnCounter = new int[1];
  Header[0] = NULL;
  Matrix[0] = NULL;
  RowCounter[0] = -1;
  MatrixCounter = 0;
  ColumnCounter[0] = -1;
};

/** Destructor of MatrixContainer class.
 */
MatrixContainer::~MatrixContainer() {
  if (Matrix != NULL) {
    // free Matrix[MatrixCounter]
    if ((ColumnCounter != NULL) && (RowCounter != NULL) && (Matrix[MatrixCounter] != NULL))
      for(int j=RowCounter[MatrixCounter]+1;j--;)
        delete[](Matrix[MatrixCounter][j]);
    //if (MatrixCounter != 0)
    delete[](Matrix[MatrixCounter]);
    // free all matrices but ultimate one
    for(int i=MatrixCounter;i--;) {
      if ((ColumnCounter != NULL) && (RowCounter != NULL)) {
          for(int j=RowCounter[i]+1;j--;)
            delete[](Matrix[i][j]);
          delete[](Matrix[i]);
      }
    }
    delete[](Matrix);
  }
  // free indices
  if (Indices != NULL)
    for(int i=MatrixCounter+1;i--;) {
      delete[](Indices[i]);
    }
  delete[](Indices);

  // free header and counters
  if (Header != NULL)
    for(int i=MatrixCounter+1;i--;)
      delete[](Header[i]);
  delete[](Header);
  delete[](RowCounter);
  delete[](ColumnCounter);
};

/** 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 *Matrix)
{
  DoLog(0) && (Log() << Verbose(0) << "Initialising indices");
  if (Matrix == NULL) {
    DoLog(0) && (Log() << Verbose(0) << " with trivial mapping." << endl);
    Indices = new int*[MatrixCounter + 1];
    for(int i=MatrixCounter+1;i--;) {
      Indices[i] = new int[RowCounter[i]];
      for(int j=RowCounter[i];j--;)
        Indices[i][j] = j;
    }
  } else {
    DoLog(0) && (Log() << Verbose(0) << " from other MatrixContainer." << endl);
    if (MatrixCounter != Matrix->MatrixCounter)
      return false;
    Indices = new int*[MatrixCounter + 1];
    for(int i=MatrixCounter+1;i--;) {
      if (RowCounter[i] != Matrix->RowCounter[i])
        return false;
      Indices[i] = new int[Matrix->RowCounter[i]];
      for(int j=Matrix->RowCounter[i];j--;) {
        Indices[i][j] = Matrix->Indices[i][j];
        //Log() << Verbose(0) << Indices[i][j] << "\t";
      }
      //Log() << Verbose(0) << endl;
    }
  }
  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 matrix
 *    -# loop over found column and row counts and parse in each entry
 * \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, int MatrixNr)
{
  stringstream line;
  string token;
  char filename[1023];

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

  // parse header
  if (Header[MatrixNr] != NULL)
    delete[] Header[MatrixNr];
  Header[MatrixNr] = new char[1024];
  for (int m=skiplines+1;m--;)
    input.getline(Header[MatrixNr], 1023);

  // scan header for number of columns
  line.str(Header[MatrixNr]);
  for(int k=skipcolumns;k--;)
    line >> Header[MatrixNr];
  Log() << Verbose(0) << line.str() << endl;
  ColumnCounter[MatrixNr]=0;
  while ( getline(line,token, '\t') ) {
    if (token.length() > 0)
      ColumnCounter[MatrixNr]++;
  }
  Log() << Verbose(0) << line.str() << endl;
  Log() << Verbose(1) << "ColumnCounter[" << MatrixNr << "]: " << ColumnCounter[MatrixNr] << "." << endl;
  if (ColumnCounter[MatrixNr] == 0) {
    DoeLog(0) && (eLog()<< Verbose(0) << "ColumnCounter[" << MatrixNr << "]: " << ColumnCounter[MatrixNr] << " from ostream." << endl);
    performCriticalExit();
  }

  // scan rest for number of rows/lines
  RowCounter[MatrixNr]=-1;    // counts one line too much
  while (!input.eof()) {
    input.getline(filename, 1023);
    Log() << Verbose(0) << "Comparing: " << strncmp(filename,"MeanForce",9) << endl;
    RowCounter[MatrixNr]++; // then line was not last MeanForce
    if (strncmp(filename,"MeanForce", 9) == 0) {
      break;
    }
  }
  Log() << Verbose(1) << "RowCounter[" << MatrixNr << "]: " << RowCounter[MatrixNr] << " from input stream." << endl;
  if (RowCounter[MatrixNr] == 0) {
    DoeLog(0) && (eLog()<< Verbose(0) << "RowCounter[" << MatrixNr << "]: " << RowCounter[MatrixNr] << " from input stream." << endl);
    performCriticalExit();
  }

  // allocate matrix if it's not zero dimension in one direction
  if ((ColumnCounter[MatrixNr] > 0) && (RowCounter[MatrixNr] > -1)) {
    if (Matrix[MatrixNr] != NULL)
      delete[] Matrix[MatrixNr];
    Matrix[MatrixNr] = new double*[RowCounter[MatrixNr] + 1];
    for(int j=0;j<RowCounter[MatrixNr]+1;j++)
      Matrix[MatrixNr][j] = 0;

    // parse in each entry for this matrix
    input.clear();
    input.seekg(ios::beg);
    for (int m=skiplines+1;m--;)
      input.getline(Header[MatrixNr], 1023);    // skip header
    line.str(Header[MatrixNr]);
    Log() << Verbose(0) << "Header: " << line.str() << endl;
    for(int k=skipcolumns;k--;)  // skip columns in header too
      line >> filename;
    strncpy(Header[MatrixNr], line.str().c_str(), 1023);
    for(int j=0;j<RowCounter[MatrixNr];j++) {
      if (Matrix[MatrixNr][j] != NULL)
        delete[] Matrix[MatrixNr][j];
      Matrix[MatrixNr][j] = new double[ColumnCounter[MatrixNr]];
      for(int k=0;k<ColumnCounter[MatrixNr];k++)
        Matrix[MatrixNr][j][k] = 0;

      input.getline(filename, 1023);
      stringstream lines(filename);
      Log() << Verbose(2) << "Matrix at level " << j << ":";// << filename << endl;
      for(int k=skipcolumns;k--;)
        lines >> filename;
      for(int k=0;(k<ColumnCounter[MatrixNr]) && (!lines.eof());k++) {
        lines >> Matrix[MatrixNr][j][k];
        Log() << Verbose(1) << " " << std::setprecision(2) << Matrix[MatrixNr][j][k] << endl;
      }
      if (Matrix[MatrixNr][ RowCounter[MatrixNr] ] != NULL)
        delete[] Matrix[MatrixNr][ RowCounter[MatrixNr] ];
      Matrix[MatrixNr][ RowCounter[MatrixNr] ] = new double[ColumnCounter[MatrixNr]];
      for(int j=ColumnCounter[MatrixNr];j--;)
        Matrix[MatrixNr][ RowCounter[MatrixNr] ][j] = 0.;
    }
  } else {
    DoeLog(1) && (eLog()<< Verbose(1) << "Matrix nr. " << MatrixNr << " has column and row count of (" << ColumnCounter[MatrixNr] << "," << RowCounter[MatrixNr] << "), could not allocate nor parse!" << endl);
  }
  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 char *name, const char *prefix, 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()) {
    DoLog(0) && (Log() << Verbose(0) << endl << "MatrixContainer::ParseFragmentMatrix: Unable to open " << file.str() << ", is the directory correct?" << endl);
    return false;
  }
  while (!input.eof()) {
    input.getline(filename, 1023);
    stringstream zeile(filename);
    MatrixCounter++;
  }
  input.close();
  DoLog(0) && (Log() << Verbose(0) << "Determined " << MatrixCounter << " fragments." << endl);

  DoLog(0) && (Log() << Verbose(0) << "Parsing through each fragment and retrieving " << prefix << suffix << "." << endl);
  delete[](Header);
  delete[](Matrix);
  delete[](RowCounter);
  delete[](ColumnCounter);
  Header = new char*[MatrixCounter + 1]; // one more each for the total molecule
  Matrix = new double**[MatrixCounter + 1]; // one more each for the total molecule
  RowCounter = new int[MatrixCounter + 1];
  ColumnCounter = new int[MatrixCounter + 1];
  for(int i=MatrixCounter+1;i--;) {
    Matrix[i] = NULL;
    Header[i] = NULL;
    RowCounter[i] = -1;
    ColumnCounter[i] = -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());
    DoLog(0) &&( Log() << Verbose(0) << "Opening " << file.str() << " ... " << endl);
    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(char **GivenHeader, int MCounter, int *RCounter, int *CCounter)
{
  MatrixCounter = MCounter;
  Header = new char*[MatrixCounter + 1];
  Matrix = new double**[MatrixCounter + 1]; // one more each for the total molecule
  RowCounter = new int[MatrixCounter + 1];
  ColumnCounter = new int[MatrixCounter + 1];
  for(int i=MatrixCounter+1;i--;) {
    Header[i] = new char[1024];
    strncpy(Header[i], GivenHeader[i], 1023);
    RowCounter[i] = RCounter[i];
    ColumnCounter[i] = CCounter[i];
    Matrix[i] = new double*[RowCounter[i] + 1];
    for(int j=RowCounter[i]+1;j--;) {
      Matrix[i][j] = new double[ColumnCounter[i]];
      for(int k=ColumnCounter[i];k--;)
        Matrix[i][j][k] = 0.;
    }
  }
  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(double **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() << Verbose(0) << "Current Fragment is " << CurrentFragment << "/" << KeySets.OrderSet[Order][CurrentFragment] << "." << endl;
    // 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() << Verbose(0) << "Current other fragment is " << j << "/" << KeySets.OrderSet[SubOrder][j] << "." << endl;
          // 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() << Verbose(0) << "Current index is " << k << "/" << m << "." << endl;
            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() << Verbose(0) << "Comparing " << m << " with " << MatrixValues.Indices[ KeySets.OrderSet[Order][CurrentFragment] ][l] << "." << endl;
                if (m == MatrixValues.Indices[ KeySets.OrderSet[Order][CurrentFragment] ][l]) {
                  m = l;
                  break;
                }
              }
              //Log() << Verbose(0) << "Corresponding index in CurrentFragment is " << m << "." << endl;
              if (m > RowCounter[ KeySets.OrderSet[Order][CurrentFragment] ]) {
                DoeLog(0) && (eLog()<< Verbose(0) << "In fragment No. " << KeySets.OrderSet[Order][CurrentFragment]   << " current force index " << m << " is greater than " << RowCounter[ KeySets.OrderSet[Order][CurrentFragment] ] << "!" << endl);
                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() << Verbose(0) << "Fragments[ KeySets.OrderSet[" << Order << "][" << CurrentFragment << "]=" << KeySets.OrderSet[Order][CurrentFragment] << " ][" << RowCounter[0]-1 << "][" << 1 << "] = " <<  Matrix[ KeySets.OrderSet[Order][CurrentFragment] ][RowCounter[0]-1][1] << endl;
          }
        } else {
          //Log() << Verbose(0) << "Fragment " << KeySets.OrderSet[SubOrder][j] << " is not contained in fragment " << KeySets.OrderSet[Order][CurrentFragment] << "." << endl;
        }
      }
    }
   //Log() << Verbose(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] << endl;
  }

  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 char *name, const char *prefix)
{
  ofstream output;
  char *FragmentNumber = NULL;

  DoLog(0) && (Log() << Verbose(0) << "Writing fragment files." << endl);
  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) {
      DoeLog(0) && (eLog()<< Verbose(0) << "MatrixContainer::WriteTotalFragments: Unable to open output energy file " << line.str() << "!" << endl);
      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 char *name, const char *prefix, const char *suffix)
{
  ofstream output;
  stringstream line;

  DoLog(0) && (Log() << Verbose(0) << "Writing matrix values of " << suffix << "." << endl);
  line << name << prefix << suffix;
  output.open(line.str().c_str(), ios::out);
  if (output == NULL) {
    DoeLog(0) && (eLog()<< Verbose(0) << "MatrixContainer::WriteLastMatrix: Unable to open output matrix file " << line.str() << "!" << endl);
    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;
};

std::ostream & operator << (std::ostream &ost, const MatrixContainer &m)
{
  for (int i=0;i<=m.MatrixCounter;++i) {
    ost << "Printing matrix " << i << ":" << 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;
}