source: src/parser.cpp@ 178f92

/** \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 =============================================