/** \file parsing.cpp
 *
 * Declarations of various class functions for the parsing of value files.
 *    
 */

// ======================================= INCLUDES ==========================================

#include "helpers.hpp" 
#include "parser.hpp"

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

// ======================================= FUNCTIONS ==========================================

/** Test if given filename can be opened.
 * \param filename name of file
 * \param test true - no error message, false - print error
 * \return given file exists
 */
bool FilePresent(const char *filename, bool test)
{
  ifstream input;
  
  input.open(filename, ios::in);
  if (input == NULL) {
    if (!test)
      cout << endl << "Unable to open " << filename << ", is the directory correct?" << endl;
    return false;
  }
  input.close();
  return true;
};

/** Test the given parameters.
 * \param argc argument count
 * \param **argv arguments array
 * \return given inputdir is valid
 */
bool TestParams(int argc, char **argv)
{
  ifstream input;
  stringstream line;

  line << argv[1] << FRAGMENTPREFIX << KEYSETFILE;
  return FilePresent(line.str().c_str(), false);
};

// ======================================= CLASS MatrixContainer =============================

/** Constructor of MatrixContainer class.
 */
MatrixContainer::MatrixContainer() {
  Indices = NULL;
  Header = (char *) Malloc(sizeof(char)*1023, "MatrixContainer::MatrixContainer: *Header");
  Matrix = (double ***) Malloc(sizeof(double **)*(1), "MatrixContainer::MatrixContainer: ***Matrix"); // one more each for the total molecule
  RowCounter = (int *) Malloc(sizeof(int)*(1), "MatrixContainer::MatrixContainer: *RowCounter");
  Matrix[0] = NULL;
  RowCounter[0] = -1;
  MatrixCounter = 0;
  ColumnCounter = 0;
};

/** Destructor of MatrixContainer class.
 */
MatrixContainer::~MatrixContainer() {
  if (Matrix != NULL) {
    for(int i=MatrixCounter;i--;) {
      if (RowCounter != NULL) {
          for(int j=RowCounter[i]+1;j--;)
            Free((void **)&Matrix[i][j], "MatrixContainer::~MatrixContainer: *Matrix[][]");
        Free((void **)&Matrix[i], "MatrixContainer::~MatrixContainer: **Matrix[]");
      }
    }
    if ((RowCounter != NULL) && (Matrix[MatrixCounter] != NULL))
      for(int j=RowCounter[MatrixCounter]+1;j--;)
        Free((void **)&Matrix[MatrixCounter][j], "MatrixContainer::~MatrixContainer: *Matrix[MatrixCounter][]");
    if (MatrixCounter != 0)
      Free((void **)&Matrix[MatrixCounter], "MatrixContainer::~MatrixContainer: **Matrix[MatrixCounter]");
    Free((void **)&Matrix, "MatrixContainer::~MatrixContainer: ***Matrix");
  }
  if (Indices != NULL)
    for(int i=MatrixCounter+1;i--;) {
      Free((void **)&Indices[i], "MatrixContainer::~MatrixContainer: *Indices[]");
    }
  Free((void **)&Indices, "MatrixContainer::~MatrixContainer: **Indices"); 
  
  Free((void **)&Header, "MatrixContainer::~MatrixContainer: *Header");
  Free((void **)&RowCounter, "MatrixContainer::~MatrixContainer: *RowCounter");
};


/** 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 *name directory with files
 * \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(const char *name, int skiplines, int skipcolumns, int MatrixNr)
{
  ifstream input;
  stringstream line;
  string token;
  char filename[1023];

  input.open(name, ios::in);
  //cout << "Opening " << name << " ... "  << input << endl;
  if (input == NULL) {
    cerr << endl << "Unable to open " << name << ", is the directory correct?" << endl;
    return false;
  }
  
  // skip some initial lines
  for (int m=skiplines+1;m--;)
    input.getline(Header, 1023);
  
  // scan header for number of columns
  line.str(Header);
  for(int k=skipcolumns;k--;)
    line >> Header;
  //cout << line.str() << endl;
  ColumnCounter=0;
  while ( getline(line,token, '\t') ) {
    if (token.length() > 0)
      ColumnCounter++;
  }
  //cout << line.str() << endl;
  //cout << "ColumnCounter: " << ColumnCounter << "." << endl;
  if (ColumnCounter == 0)
    cerr << "ColumnCounter: " << ColumnCounter << " from file " << name << ", this is probably an error!" << endl;
  
  // scan rest for number of rows/lines
  RowCounter[MatrixNr]=-1;    // counts one line too much
  while (!input.eof()) { 
    input.getline(filename, 1023);
    //cout << "Comparing: " << strncmp(filename,"MeanForce",9) << endl;
    RowCounter[MatrixNr]++; // then line was not last MeanForce
    if (strncmp(filename,"MeanForce", 9) == 0) {
      break;
    }
  }
  //cout << "RowCounter[" << MatrixNr << "]: " << RowCounter[MatrixNr] << " from file " << name << "." << endl;
  if (RowCounter[MatrixNr] == 0)
    cerr << "RowCounter[" << MatrixNr << "]: " << RowCounter[MatrixNr] << " from file " << name << ", this is probably an error!" << endl;
  
  // allocate matrix if it's not zero dimension in one direction
  if ((ColumnCounter > 0) && (RowCounter[MatrixNr] > -1)) {
    Matrix[MatrixNr] = (double **) Malloc(sizeof(double *)*(RowCounter[MatrixNr]+1), "MatrixContainer::ParseFragmentMatrix: **Matrix[]");
  
    // parse in each entry for this matrix
    input.clear();
    input.seekg(ios::beg);
    for (int m=skiplines+1;m--;)
      input.getline(Header, 1023);    // skip header
    line.str(Header);
    for(int k=skipcolumns;k--;)  // skip columns in header too
      line >> filename;
    strncpy(Header, line.str().c_str(), 1023);  
    for(int j=0;j<RowCounter[MatrixNr];j++) {
      Matrix[MatrixNr][j] = (double *) Malloc(sizeof(double)*ColumnCounter, "MatrixContainer::ParseFragmentMatrix: *Matrix[][]");
      input.getline(filename, 1023);
      stringstream lines(filename);
      //cout << "Matrix at level " << j << ":";// << filename << endl;
      for(int k=skipcolumns;k--;)
        lines >> filename;
      for(int k=0;(k<ColumnCounter) && (!lines.eof());k++) {
        lines >> Matrix[MatrixNr][j][k];
        //cout << " " << setprecision(2) << Matrix[MatrixNr][j][k];;
      }
      //cout << endl;
      Matrix[MatrixNr][ RowCounter[MatrixNr] ] = (double *) Malloc(sizeof(double)*ColumnCounter, "MatrixContainer::ParseFragmentMatrix: *Matrix[RowCounter[MatrixNr]][]");
      for(int j=ColumnCounter;j--;)
        Matrix[MatrixNr][ RowCounter[MatrixNr] ][j] = 0.;
    }
  } else {
    cerr << "ERROR: Matrix nr. " << MatrixNr << " has column and row count of (" << ColumnCounter << "," << RowCounter[MatrixNr] << "), could not allocate nor parse!" << endl; 
  }
  input.close();
  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(char *name, 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 == NULL) {
    cout << endl << "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();  
  cout << "Determined " << MatrixCounter << " fragments." << endl;
  
  cout << "Parsing through each fragment and retrieving " << prefix << suffix << "." << endl;
  Matrix = (double ***) ReAlloc(Matrix, sizeof(double **)*(MatrixCounter+1), "MatrixContainer::ParseFragmentMatrix: ***Matrix"); // one more each for the total molecule
  RowCounter = (int *) ReAlloc(RowCounter, sizeof(int)*(MatrixCounter+1), "MatrixContainer::ParseFragmentMatrix: *RowCounter");
  for(int i=MatrixCounter+1;i--;) {
    Matrix[i] = NULL;
    RowCounter[i] = -1;
  }
  for(int i=0; i < MatrixCounter;i++) {
    // open matrix file
    FragmentNumber = FixedDigitNumber(MatrixCounter, i);
    file.str(" ");
    file << name << FRAGMENTPREFIX << FragmentNumber << prefix << suffix;
    if (!ParseMatrix(file.str().c_str(), skiplines, skipcolumns, i))
      return false;
    Free((void **)&FragmentNumber, "MatrixContainer::ParseFragmentMatrix: *FragmentNumber");
  }
  return true;
};

/** Allocates and resets the memory for a number \a MCounter of matrices.
 * \param *GivenHeader Header line
 * \param MCounter number of matrices
 * \param *RCounter number of rows for each matrix
 * \param CCounter number of columns for all matrices
 * \return Allocation successful 
 */
bool MatrixContainer::AllocateMatrix(char *GivenHeader, int MCounter, int *RCounter, int CCounter)
{
  Header = (char *) Malloc(sizeof(char)*1024, "MatrixContainer::ParseFragmentMatrix: *EnergyHeader");
  strncpy(Header, GivenHeader, 1023);

  MatrixCounter = MCounter;
  ColumnCounter = CCounter;
  Matrix = (double ***) Malloc(sizeof(double **)*(MatrixCounter+1), "MatrixContainer::ParseFragmentMatrix: ***Matrix"); // one more each for the total molecule
  RowCounter = (int *) Malloc(sizeof(int)*(MatrixCounter+1), "MatrixContainer::ParseFragmentMatrix: *RowCounter");
  for(int i=MatrixCounter+1;i--;) {
    RowCounter[i] = RCounter[i];
    Matrix[i] = (double **) Malloc(sizeof(double *)*(RowCounter[i]+1), "MatrixContainer::ParseFragmentMatrix: **Matrix[]");  
    for(int j=RowCounter[i]+1;j--;) {
      Matrix[i][j] = (double *) Malloc(sizeof(double)*ColumnCounter, "MatrixContainer::ParseFragmentMatrix: *Matrix[][]");
      for(int k=ColumnCounter;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;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;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;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;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;k++)
      Matrix[MatrixCounter][j][k] = values[j][k];
   return true;
};

/** Sums the energy with each factor and put into last element of \a ***Energies.
 * Sums over "E"-terms to create the "F"-terms
 * \param Matrix MatrixContainer with matrices (LevelCounter by ColumnCounter) with all the energies.
 * \param KeySet 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 &KeySet, int Order)
{
  // go through each order
  for (int CurrentFragment=0;CurrentFragment<KeySet.FragmentsPerOrder[Order];CurrentFragment++) {
    //cout << "Current Fragment is " << CurrentFragment << "/" << KeySet.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<KeySet.FragmentsPerOrder[SubOrder];j++) { // go through all possible fragments of size suborder
        if (KeySet.Contains(KeySet.OrderSet[Order][CurrentFragment], KeySet.OrderSet[SubOrder][j])) {
          //cout << "Current other fragment is " << j << "/" << KeySet.OrderSet[SubOrder][j] << "." << endl;
          // if the fragment's indices are all in the current fragment
          for(int k=0;k<RowCounter[ KeySet.OrderSet[SubOrder][j] ];k++) { // go through all atoms in this fragment
            int m = MatrixValues.Indices[ KeySet.OrderSet[SubOrder][j] ][k];
            //cout << "Current index is " << k << "/" << m << "." << endl;
            if (m != -1) { // if it's not an added hydrogen
              for (int l=0;l<RowCounter[ KeySet.OrderSet[Order][CurrentFragment] ];l++) { // look for the corresponding index in the current fragment
                //cout << "Comparing " << m << " with " << MatrixValues.Indices[ KeySet.OrderSet[Order][CurrentFragment] ][l] << "." << endl;
                if (m == MatrixValues.Indices[ KeySet.OrderSet[Order][CurrentFragment] ][l]) {
                  m = l;
                  break;  
                }
              }
              //cout << "Corresponding index in CurrentFragment is " << m << "." << endl;
              if (m > RowCounter[ KeySet.OrderSet[Order][CurrentFragment] ]) {
                cerr << "In fragment No. " << KeySet.OrderSet[Order][CurrentFragment]   << " current force index " << m << " is greater than " << RowCounter[ KeySet.OrderSet[Order][CurrentFragment] ] << "!" << endl;
                return false;
              }
              if (Order == SubOrder) { // equal order is always copy from Energies
                for(int l=ColumnCounter;l--;) // then adds/subtract each column 
                  Matrix[ KeySet.OrderSet[Order][CurrentFragment] ][m][l] += MatrixValues.Matrix[ KeySet.OrderSet[SubOrder][j] ][k][l];
              } else {
                for(int l=ColumnCounter;l--;)
                  Matrix[ KeySet.OrderSet[Order][CurrentFragment] ][m][l] -= Matrix[ KeySet.OrderSet[SubOrder][j] ][k][l];
              }
            }
            //if ((ColumnCounter>1) && (RowCounter[0]-1 >= 1))
             //cout << "Fragments[ KeySet.OrderSet[" << Order << "][" << CurrentFragment << "]=" << KeySet.OrderSet[Order][CurrentFragment] << " ][" << RowCounter[0]-1 << "][" << 1 << "] = " <<  Matrix[ KeySet.OrderSet[Order][CurrentFragment] ][RowCounter[0]-1][1] << endl;
          }
        } else {
          //cout << "Fragment " << KeySet.OrderSet[SubOrder][j] << " is not contained in fragment " << KeySet.OrderSet[Order][CurrentFragment] << "." << endl;
        }
      }
    }
   //cout << "Final Fragments[ KeySet.OrderSet[" << Order << "][" << CurrentFragment << "]=" << KeySet.OrderSet[Order][CurrentFragment] << " ][" << KeySet.AtomCounter[0]-1 << "][" << 1 << "] = " <<  Matrix[ KeySet.OrderSet[Order][CurrentFragment] ][KeySet.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;

  cout << "Writing fragment files." << endl;
  for(int i=0;i<MatrixCounter;i++) {
    stringstream line;
    FragmentNumber = FixedDigitNumber(MatrixCounter, i);
    line << name << FRAGMENTPREFIX << FragmentNumber << "/" << prefix;
    Free((void **)&FragmentNumber, "*FragmentNumber");
    output.open(line.str().c_str(), ios::out);
    if (output == NULL) {
      cerr << "Unable to open output energy file " << line.str() << "!" << endl;
      return false;
    }
    output << Header << endl;
    for(int j=0;j<RowCounter[i];j++) {
      for(int k=0;k<ColumnCounter;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;

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

// ======================================= CLASS EnergyMatrix =============================

/** Create a trivial energy index mapping.
 * This just maps 1 to 1, 2 to 2 and so on for all fragments.
 * \return creation sucessful
 */
bool EnergyMatrix::ParseIndices() 
{
  cout << "Parsing energy indices." << endl;
  Indices = (int **) Malloc(sizeof(int *)*(MatrixCounter+1), "EnergyMatrix::ParseIndices: **Indices");
  for(int i=MatrixCounter+1;i--;) {
    Indices[i] = (int *) Malloc(sizeof(int)*RowCounter[i], "EnergyMatrix::ParseIndices: *Indices[]");
    for(int j=RowCounter[i];j--;)
      Indices[i][j] = j;
  } 
  return true; 
};

/** Sums the energy with each factor and put into last element of \a EnergyMatrix::Matrix.
 * Sums over the "F"-terms in ANOVA decomposition.
 * \param Matrix MatrixContainer with matrices (LevelCounter by ColumnCounter) with all the energies.
 * \param CorrectionFragments MatrixContainer with hydrogen saturation correction per fragments
 * \param KeySet KeySetContainer with bond Order and association mapping of each fragment to an order
 * \param Order bond order
 * \parsm sign +1 or -1
 * \return true if summing was successful
 */
bool EnergyMatrix::SumSubEnergy(class EnergyMatrix &Fragments, class EnergyMatrix *CorrectionFragments, class KeySetsContainer &KeySet, int Order, double sign)
{
  // sum energy
  if (CorrectionFragments == NULL)
    for(int i=KeySet.FragmentsPerOrder[Order];i--;)
      for(int j=RowCounter[ KeySet.OrderSet[Order][i] ];j--;)
        for(int k=ColumnCounter;k--;)
          Matrix[MatrixCounter][j][k] += sign*Fragments.Matrix[ KeySet.OrderSet[Order][i] ][j][k];
  else
    for(int i=KeySet.FragmentsPerOrder[Order];i--;)
      for(int j=RowCounter[ KeySet.OrderSet[Order][i] ];j--;)
        for(int k=ColumnCounter;k--;)
          Matrix[MatrixCounter][j][k] += sign*(Fragments.Matrix[ KeySet.OrderSet[Order][i] ][j][k] + CorrectionFragments->Matrix[ KeySet.OrderSet[Order][i] ][j][k]);
  return true;
};

/** Calls MatrixContainer::ParseFragmentMatrix() and additionally allocates last plus one matrix.
 * \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 EnergyMatrix::ParseFragmentMatrix(char *name, char *prefix, string suffix, int skiplines, int skipcolumns)
{
  char filename[1024];
  bool status = MatrixContainer::ParseFragmentMatrix(name, prefix, suffix, skiplines, skipcolumns);

  if (status) {
    // count maximum of columns
    RowCounter[MatrixCounter] = 0;
    for(int j=0; j < MatrixCounter;j++) // (energy matrix might be bigger than number of atoms in terms of rows)
      if (RowCounter[j] > RowCounter[MatrixCounter])
        RowCounter[MatrixCounter] = RowCounter[j];
    // allocate last plus one matrix
    cout << "Allocating last plus one matrix with " << (RowCounter[MatrixCounter]+1) << " rows and " << ColumnCounter << " columns." << endl;
    Matrix[MatrixCounter] = (double **) Malloc(sizeof(double *)*(RowCounter[MatrixCounter]+1), "MatrixContainer::ParseFragmentMatrix: **Matrix[]");
    for(int j=0;j<=RowCounter[MatrixCounter];j++)
      Matrix[MatrixCounter][j] = (double *) Malloc(sizeof(double)*ColumnCounter, "MatrixContainer::ParseFragmentMatrix: *Matrix[][]");
    
    // try independently to parse global energysuffix file if present
    strncpy(filename, name, 1023);
    strncat(filename, prefix, 1023-strlen(filename));
    strncat(filename, suffix.c_str(), 1023-strlen(filename));
    ParseMatrix(filename, skiplines, skipcolumns, MatrixCounter);
  }
  return status;
};

// ======================================= CLASS ForceMatrix =============================

/** Parsing force Indices of each fragment
 * \param *name directory with \a ForcesFile
 * \return parsing successful 
 */
bool ForceMatrix::ParseIndices(char *name) 
{
  ifstream input;
  char *FragmentNumber = NULL;
  char filename[1023];
  stringstream line;
  
  cout << "Parsing force indices for " << MatrixCounter << " matrices." << endl;
  Indices = (int **) Malloc(sizeof(int *)*(MatrixCounter+1), "ForceMatrix::ParseIndices: **Indices");
  line << name << FRAGMENTPREFIX << FORCESFILE;
  input.open(line.str().c_str(), ios::in);
  //cout << "Opening " << line.str() << " ... "  << input << endl;
  if (input == NULL) {
    cout << endl << "Unable to open " << line.str() << ", is the directory correct?" << endl;
    return false;
  }
  for (int i=0;(i<MatrixCounter) && (!input.eof());i++) {
    // get the number of atoms for this fragment
    input.getline(filename, 1023);
    line.str(filename);
    // parse the values
    Indices[i] = (int *) Malloc(sizeof(int)*RowCounter[i], "ForceMatrix::ParseIndices: *Indices[]");
    FragmentNumber = FixedDigitNumber(MatrixCounter, i);
    //cout << FRAGMENTPREFIX << FragmentNumber << "[" << RowCounter[i] << "]:";
    Free((void **)&FragmentNumber, "ForceMatrix::ParseIndices: *FragmentNumber");
    for(int j=0;(j<RowCounter[i]) && (!line.eof());j++) {
      line >> Indices[i][j];
      //cout << " " << Indices[i][j];
    }
    //cout << endl;
  }
  Indices[MatrixCounter] = (int *) Malloc(sizeof(int)*RowCounter[MatrixCounter], "ForceMatrix::ParseIndices: *Indices[]");
  for(int j=RowCounter[MatrixCounter];j--;) {
    Indices[MatrixCounter][j] = j;
  }
  input.close();
  return true;
};


/** Sums the forces and puts into last element of \a ForceMatrix::Matrix.
 * \param Matrix MatrixContainer with matrices (LevelCounter by ColumnCounter) with all the energies.
 * \param KeySet KeySetContainer with bond Order and association mapping of each fragment to an order
 * \param Order bond order
 *  \param sign +1 or -1
 * \return true if summing was successful
 */
bool ForceMatrix::SumSubForces(class ForceMatrix &Fragments, class KeySetsContainer &KeySet, int Order, double sign)
{
  int FragmentNr;
  // sum forces
  for(int i=0;i<KeySet.FragmentsPerOrder[Order];i++) {
    FragmentNr = KeySet.OrderSet[Order][i];
    for(int l=0;l<RowCounter[ FragmentNr ];l++) {
      int j = Indices[ FragmentNr ][l];
      if (j > RowCounter[MatrixCounter]) {
        cerr << "Current force index " << j << " is greater than " << RowCounter[MatrixCounter] << "!" << endl;
        return false;
      }
      if (j != -1) {
        //if (j == 0) cout << "Summing onto ion 0, type 0 from fragment " << FragmentNr << ", ion " << l << "." << endl;
        for(int k=2;k<ColumnCounter;k++)
          Matrix[MatrixCounter][j][k] += sign*Fragments.Matrix[ FragmentNr ][l][k];
      }
    }
  }
  return true;
};


/** Calls MatrixContainer::ParseFragmentMatrix() and additionally allocates last plus one matrix.
 * \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 ForceMatrix::ParseFragmentMatrix(char *name, char *prefix, string suffix, int skiplines, int skipcolumns)
{
  char filename[1023];
  ifstream input;
  stringstream file;
  int nr;
  bool status = MatrixContainer::ParseFragmentMatrix(name, prefix, suffix, skiplines, skipcolumns);

  if (status) {
    // count number of atoms for last plus one matrix
    file << name << FRAGMENTPREFIX << KEYSETFILE;
    input.open(file.str().c_str(), ios::in);
    if (input == NULL) {
      cout << endl << "Unable to open " << file.str() << ", is the directory correct?" << endl;
      return false;
    }
    RowCounter[MatrixCounter] = 0;
    while (!input.eof()) {
      input.getline(filename, 1023);
      stringstream zeile(filename);
      while (!zeile.eof()) {
        zeile >> nr;
        //cout << "Current index: " << nr << "." << endl;
        if (nr > RowCounter[MatrixCounter])
          RowCounter[MatrixCounter] = nr;
      }
    }
    RowCounter[MatrixCounter]++;    // nr start at 0, count starts at 1
    input.close();
  
    // allocate last plus one matrix
    cout << "Allocating last plus one matrix with " << (RowCounter[MatrixCounter]+1) << " rows and " << ColumnCounter << " columns." << endl;
    Matrix[MatrixCounter] = (double **) Malloc(sizeof(double *)*(RowCounter[MatrixCounter]+1), "MatrixContainer::ParseFragmentMatrix: **Matrix[]");
    for(int j=0;j<=RowCounter[MatrixCounter];j++)
      Matrix[MatrixCounter][j] = (double *) Malloc(sizeof(double)*ColumnCounter, "MatrixContainer::ParseFragmentMatrix: *Matrix[][]");

    // try independently to parse global forcesuffix file if present
    strncpy(filename, name, 1023);
    strncat(filename, prefix, 1023-strlen(filename));
    strncat(filename, suffix.c_str(), 1023-strlen(filename));
    ParseMatrix(filename, skiplines, skipcolumns, MatrixCounter);
  }
 

  return status;
};

// ======================================= CLASS KeySetsContainer =============================

/** Constructor of KeySetsContainer class.
 */
KeySetsContainer::KeySetsContainer() {
  KeySets = NULL;
  AtomCounter = NULL;
  FragmentCounter = 0;
  Order = 0;
  FragmentsPerOrder = 0;
  OrderSet = NULL;
};

/** Destructor of KeySetsContainer class.
 */
KeySetsContainer::~KeySetsContainer() {
  for(int i=FragmentCounter;i--;)
    Free((void **)&KeySets[i], "KeySetsContainer::~KeySetsContainer: *KeySets[]");
  for(int i=Order;i--;)
    Free((void **)&OrderSet[i], "KeySetsContainer::~KeySetsContainer: *OrderSet[]");
  Free((void **)&KeySets, "KeySetsContainer::~KeySetsContainer: **KeySets");
  Free((void **)&OrderSet, "KeySetsContainer::~KeySetsContainer: **OrderSet");
  Free((void **)&AtomCounter, "KeySetsContainer::~KeySetsContainer: *AtomCounter");
  Free((void **)&FragmentsPerOrder, "KeySetsContainer::~KeySetsContainer: *FragmentsPerOrder");
};

/** Parsing KeySets into array.
 * \param *name directory with keyset file
 * \param *ACounter number of atoms per fragment
 * \param FCounter number of fragments
 * \return parsing succesful
 */
bool KeySetsContainer::ParseKeySets(const char *name, const int *ACounter, const int FCounter) {
  ifstream input;
  char *FragmentNumber = NULL;
  stringstream file;
  char filename[1023];
  
  FragmentCounter = FCounter;
  cout << "Parsing key sets." << endl;
  KeySets = (int **) Malloc(sizeof(int *)*FragmentCounter, "KeySetsContainer::ParseKeySets: **KeySets");
  for(int i=FragmentCounter;i--;)
    KeySets[i] = NULL;
  file << name << FRAGMENTPREFIX << KEYSETFILE;
  input.open(file.str().c_str(), ios::in);
  if (input == NULL) {
    cout << endl << "Unable to open " << file.str() << ", is the directory correct?" << endl;
    return false;
  }
  
  AtomCounter = (int *) Malloc(sizeof(int)*FragmentCounter, "KeySetsContainer::ParseKeySets: *RowCounter");
  for(int i=0;(i<FragmentCounter) && (!input.eof());i++) {
    stringstream line;
    AtomCounter[i] = ACounter[i];
    // parse the values
    KeySets[i] = (int *) Malloc(sizeof(int)*AtomCounter[i], "KeySetsContainer::ParseKeySets: *KeySets[]");
    for(int j=AtomCounter[i];j--;)
      KeySets[i][j] = -1;
    FragmentNumber = FixedDigitNumber(FragmentCounter, i);
    //cout << FRAGMENTPREFIX << FragmentNumber << "[" << AtomCounter[i] << "]:";
    Free((void **)&FragmentNumber, "KeySetsContainer::ParseKeySets: *FragmentNumber");
    input.getline(filename, 1023);
    line.str(filename);
    for(int j=0;(j<AtomCounter[i]) && (!line.eof());j++) {
      line >> KeySets[i][j];
      //cout << " " << KeySets[i][j];
    }
    //cout << endl;
  }
  input.close();
  return true;
};

/** Parse many body terms, associating each fragment to a certain bond order.
 * \return parsing succesful
 */
bool KeySetsContainer::ParseManyBodyTerms()
{
  int Counter;
  
  cout << "Creating Fragment terms." << endl;
  // scan through all to determine maximum order
  Order=0;
  for(int i=FragmentCounter;i--;) {
    Counter=0;
    for(int j=AtomCounter[i];j--;)
      if (KeySets[i][j] != -1)
        Counter++;
    if (Counter > Order)
      Order = Counter;
  }
  cout << "Found Order is " << Order << "." << endl;
  
  // scan through all to determine fragments per order
  FragmentsPerOrder = (int *) Malloc(sizeof(int)*Order, "KeySetsContainer::ParseManyBodyTerms: *FragmentsPerOrder");
  for(int i=Order;i--;)
    FragmentsPerOrder[i] = 0;
  for(int i=FragmentCounter;i--;) {
    Counter=0;
    for(int j=AtomCounter[i];j--;)
      if (KeySets[i][j] != -1)
        Counter++;
    FragmentsPerOrder[Counter-1]++;
  }
  for(int i=0;i<Order;i++)
    cout << "Found No. of Fragments of Order " << i+1 << " is " << FragmentsPerOrder[i] << "." << endl;
    
  // scan through all to gather indices to each order set
  OrderSet = (int **) Malloc(sizeof(int *)*Order, "KeySetsContainer::ParseManyBodyTerms: **OrderSet");
  for(int i=Order;i--;)
    OrderSet[i] = (int *) Malloc(sizeof(int)*FragmentsPerOrder[i], "KeySetsContainer::ParseManyBodyTermsKeySetsContainer::ParseManyBodyTerms: *OrderSet[]");
  for(int i=Order;i--;)
    FragmentsPerOrder[i] = 0;
  for(int i=FragmentCounter;i--;) {
    Counter=0;
    for(int j=AtomCounter[i];j--;)
      if (KeySets[i][j] != -1)
        Counter++;
    OrderSet[Counter-1][FragmentsPerOrder[Counter-1]] = i;
    FragmentsPerOrder[Counter-1]++;
  }
  cout << "Printing OrderSet." << endl;
  for(int i=0;i<Order;i++) {
    for (int j=0;j<FragmentsPerOrder[i];j++) {
      cout << " " << OrderSet[i][j];
    }
    cout << endl;
  }
  cout << endl;
  
    
  return true;
};

/** Compares each entry in \a *SmallerSet if it is containted in \a *GreaterSet.
 * \param GreaterSet index to greater set
 * \param SmallerSet index to smaller set
 * \return true if all keys of SmallerSet contained in GreaterSet
 */
bool KeySetsContainer::Contains(const int GreaterSet, const int SmallerSet)
{
  bool result = true;
  bool intermediate;
  if ((GreaterSet < 0) || (SmallerSet < 0) || (GreaterSet > FragmentCounter) || (SmallerSet > FragmentCounter)) // index out of bounds
    return false;
  for(int i=AtomCounter[SmallerSet];i--;) {
    intermediate = false;
    for (int j=AtomCounter[GreaterSet];j--;)
      intermediate = (intermediate || ((KeySets[SmallerSet][i] == KeySets[GreaterSet][j]) || (KeySets[SmallerSet][i] == -1)));
    result = result && intermediate;
  }

  return result;
};


// ======================================= END =============================================