Ignore:
File:
1 edited

Legend:

Unmodified
Added
Removed

  • TabularUnified src/parser.cpp

    r6ac7ee rce5ac3  
    22 *
    33 * Declarations of various class functions for the parsing of value files.
    4  *
     4 *   
    55 */
    66
    77// ======================================= INCLUDES ==========================================
    88
    9 #include "helpers.hpp"
     9#include "helpers.hpp" 
    1010#include "parser.hpp"
    1111
    … …  
    2424bool FilePresent(const char *filename, bool test)
    2525{
    26         ifstream input;
    27 
    28         input.open(filename, ios::in);
    29         if (input == NULL) {
    30                 if (!test)
    31                         cout << endl << "Unable to open " << filename << ", is the directory correct?" << endl;
    32                 return false;
    33         }
    34         input.close();
    35         return true;
     26  ifstream input;
     27 
     28  input.open(filename, ios::in);
     29  if (input == NULL) {
     30    if (!test)
     31      cout << endl << "Unable to open " << filename << ", is the directory correct?" << endl;
     32    return false;
     33  }
     34  input.close();
     35  return true;
    3636};
    3737
    … …  
    4343bool TestParams(int argc, char **argv)
    4444{
    45         ifstream input;
    46         stringstream line;
    47 
    48         line << argv[1] << FRAGMENTPREFIX << KEYSETFILE;
    49         return FilePresent(line.str().c_str(), false);
     45  ifstream input;
     46  stringstream line;
     47
     48  line << argv[1] << FRAGMENTPREFIX << KEYSETFILE;
     49  return FilePresent(line.str().c_str(), false);
    5050};
    5151
    … …  
    5555 */
    5656MatrixContainer::MatrixContainer() {
    57         Indices = NULL;
    58         Header = (char *) Malloc(sizeof(char)*1023, "MatrixContainer::MatrixContainer: *Header");
    59         Matrix = (double ***) Malloc(sizeof(double **)*(1), "MatrixContainer::MatrixContainer: ***Matrix"); // one more each for the total molecule
    60         RowCounter = (int *) Malloc(sizeof(int)*(1), "MatrixContainer::MatrixContainer: *RowCounter");
    61         Matrix[0] = NULL;
    62         RowCounter[0] = -1;
    63         MatrixCounter = 0;
    64         ColumnCounter = 0;
     57  Indices = NULL;
     58  Header = (char **) Malloc(sizeof(char)*1, "MatrixContainer::MatrixContainer: **Header");
     59  Matrix = (double ***) Malloc(sizeof(double **)*(1), "MatrixContainer::MatrixContainer: ***Matrix"); // one more each for the total molecule
     60  RowCounter = (int *) Malloc(sizeof(int)*(1), "MatrixContainer::MatrixContainer: *RowCounter");
     61  ColumnCounter = (int *) Malloc(sizeof(int)*(1), "MatrixContainer::MatrixContainer: *ColumnCounter");
     62  Header[0] = NULL;
     63  Matrix[0] = NULL;
     64  RowCounter[0] = -1;
     65  MatrixCounter = 0;
     66  ColumnCounter[0] = -1;
    6567};
    6668
    … …  
    6870 */
    6971MatrixContainer::~MatrixContainer() {
    70         if (Matrix != NULL) {
    71                 for(int i=MatrixCounter;i--;) {
    72                         if (RowCounter != NULL) {
    73                                         for(int j=RowCounter[i]+1;j--;)
    74                                                 Free((void **)&Matrix[i][j], "MatrixContainer::~MatrixContainer: *Matrix[][]");
    75                                 Free((void **)&Matrix[i], "MatrixContainer::~MatrixContainer: **Matrix[]");
    76                         }
    77                 }
    78                 if ((RowCounter != NULL) && (Matrix[MatrixCounter] != NULL))
    79                         for(int j=RowCounter[MatrixCounter]+1;j--;)
    80                                 Free((void **)&Matrix[MatrixCounter][j], "MatrixContainer::~MatrixContainer: *Matrix[MatrixCounter][]");
    81                 if (MatrixCounter != 0)
    82                         Free((void **)&Matrix[MatrixCounter], "MatrixContainer::~MatrixContainer: **Matrix[MatrixCounter]");
    83                 Free((void **)&Matrix, "MatrixContainer::~MatrixContainer: ***Matrix");
    84         }
    85         if (Indices != NULL)
    86                 for(int i=MatrixCounter+1;i--;) {
    87                         Free((void **)&Indices[i], "MatrixContainer::~MatrixContainer: *Indices[]");
    88                 }
    89         Free((void **)&Indices, "MatrixContainer::~MatrixContainer: **Indices");
    90 
    91         Free((void **)&Header, "MatrixContainer::~MatrixContainer: *Header");
    92         Free((void **)&RowCounter, "MatrixContainer::~MatrixContainer: *RowCounter");
    93 };
    94 
     72  if (Matrix != NULL) {
     73    for(int i=MatrixCounter;i--;) {
     74      if ((ColumnCounter != NULL) && (RowCounter != NULL)) {
     75          for(int j=RowCounter[i]+1;j--;)
     76            Free((void **)&Matrix[i][j], "MatrixContainer::~MatrixContainer: *Matrix[][]");
     77        Free((void **)&Matrix[i], "MatrixContainer::~MatrixContainer: **Matrix[]");
     78      }
     79    }
     80    if ((ColumnCounter != NULL) && (RowCounter != NULL) && (Matrix[MatrixCounter] != NULL))
     81      for(int j=RowCounter[MatrixCounter]+1;j--;)
     82        Free((void **)&Matrix[MatrixCounter][j], "MatrixContainer::~MatrixContainer: *Matrix[MatrixCounter][]");
     83    if (MatrixCounter != 0)
     84      Free((void **)&Matrix[MatrixCounter], "MatrixContainer::~MatrixContainer: **Matrix[MatrixCounter]");
     85    Free((void **)&Matrix, "MatrixContainer::~MatrixContainer: ***Matrix");
     86  }
     87  if (Indices != NULL)
     88    for(int i=MatrixCounter+1;i--;) {
     89      Free((void **)&Indices[i], "MatrixContainer::~MatrixContainer: *Indices[]");
     90    }
     91  Free((void **)&Indices, "MatrixContainer::~MatrixContainer: **Indices");
     92 
     93  if (Header != NULL)
     94    for(int i=MatrixCounter+1;i--;)
     95      Free((void **)&Header[i], "MatrixContainer::~MatrixContainer: *Header[]");
     96  Free((void **)&Header, "MatrixContainer::~MatrixContainer: **Header");
     97  Free((void **)&RowCounter, "MatrixContainer::~MatrixContainer: *RowCounter");
     98  Free((void **)&ColumnCounter, "MatrixContainer::~MatrixContainer: *RowCounter");
     99};
     100
     101/** Either copies index matrix from another MatrixContainer or initialises with trivial mapping if NULL.
     102 * This either copies the index matrix or just maps 1 to 1, 2 to 2 and so on for all fragments.
     103 * \param *Matrix pointer to other MatrixContainer
     104 * \return true - copy/initialisation sucessful, false - dimension false for copying
     105 */
     106bool MatrixContainer::InitialiseIndices(class MatrixContainer *Matrix)
     107{
     108  cout << "Initialising indices";
     109  if (Matrix == NULL) {
     110    cout << " with trivial mapping." << endl;
     111    Indices = (int **) Malloc(sizeof(int *)*(MatrixCounter+1), "MatrixContainer::InitialiseIndices: **Indices");
     112    for(int i=MatrixCounter+1;i--;) {
     113      Indices[i] = (int *) Malloc(sizeof(int)*RowCounter[i], "MatrixContainer::InitialiseIndices: *Indices[]");
     114      for(int j=RowCounter[i];j--;)
     115        Indices[i][j] = j;
     116    }
     117  } else {
     118    cout << " from other MatrixContainer." << endl;
     119    if (MatrixCounter != Matrix->MatrixCounter)
     120      return false;
     121    Indices = (int **) Malloc(sizeof(int *)*(MatrixCounter+1), "MatrixContainer::InitialiseIndices: **Indices");
     122    for(int i=MatrixCounter+1;i--;) {
     123      if (RowCounter[i] != Matrix->RowCounter[i])
     124        return false;
     125      Indices[i] = (int *) Malloc(sizeof(int)*Matrix->RowCounter[i], "MatrixContainer::InitialiseIndices: *Indices[]");
     126      for(int j=Matrix->RowCounter[i];j--;) {
     127        Indices[i][j] = Matrix->Indices[i][j];
     128        //cout << Indices[i][j] << "\t";
     129      }
     130      //cout << endl;
     131    }
     132  }
     133  return true;
     134};
    95135
    96136/** Parsing a number of matrices.
    97  *              -# open the matrix file
    98  *              -# skip some lines (\a skiplines)
    99  *              -# scan header lines for number of columns
    100  *              -# scan lines for number of rows
    101  *              -# allocate matrix
    102  *              -# loop over found column and row counts and parse in each entry
     137 *    -# open the matrix file
     138 *    -# skip some lines (\a skiplines)
     139 *    -# scan header lines for number of columns
     140 *    -# scan lines for number of rows
     141 *    -# allocate matrix
     142 *    -# loop over found column and row counts and parse in each entry
    103143 * \param *name directory with files
    104144 * \param skiplines number of inital lines to skip
    105145 * \param skiplines number of inital columns to skip
    106146 * \param MatrixNr index number in Matrix array to parse into
    107  * \return parsing successful
     147 * \return parsing successful 
    108148 */
    109149bool MatrixContainer::ParseMatrix(const char *name, int skiplines, int skipcolumns, int MatrixNr)
    110150{
    111         ifstream input;
    112         stringstream line;
    113         string token;
    114         char filename[1023];
    115 
    116         input.open(name, ios::in);
    117         //cout << "Opening " << name << " ... " << input << endl;
    118         if (input == NULL) {
    119                 cerr << endl << "Unable to open " << name << ", is the directory correct?" << endl;
    120                 return false;
    121         }
    122 
    123         // skip some initial lines
    124         for (int m=skiplines+1;m--;)
    125                 input.getline(Header, 1023);
    126 
    127         // scan header for number of columns
    128         line.str(Header);
    129         for(int k=skipcolumns;k--;)
    130                 line >> Header;
    131         //cout << line.str() << endl;
    132         ColumnCounter=0;
    133         while ( getline(line,token, '\t') ) {
    134                 if (token.length() > 0)
    135                         ColumnCounter++;
    136         }
    137         //cout << line.str() << endl;
    138         //cout << "ColumnCounter: " << ColumnCounter << "." << endl;
    139         if (ColumnCounter == 0)
    140                 cerr << "ColumnCounter: " << ColumnCounter << " from file " << name << ", this is probably an error!" << endl;
    141 
    142         // scan rest for number of rows/lines
    143         RowCounter[MatrixNr]=-1;                // counts one line too much
    144         while (!input.eof()) {
    145                 input.getline(filename, 1023);
    146                 //cout << "Comparing: " << strncmp(filename,"MeanForce",9) << endl;
    147                 RowCounter[MatrixNr]++; // then line was not last MeanForce
    148                 if (strncmp(filename,"MeanForce", 9) == 0) {
    149                         break;
    150                 }
    151         }
    152         //cout << "RowCounter[" << MatrixNr << "]: " << RowCounter[MatrixNr] << " from file " << name << "." << endl;
    153         if (RowCounter[MatrixNr] == 0)
    154                 cerr << "RowCounter[" << MatrixNr << "]: " << RowCounter[MatrixNr] << " from file " << name << ", this is probably an error!" << endl;
    155 
    156         // allocate matrix if it's not zero dimension in one direction
    157         if ((ColumnCounter > 0) && (RowCounter[MatrixNr] > -1)) {
    158                 Matrix[MatrixNr] = (double **) Malloc(sizeof(double *)*(RowCounter[MatrixNr]+1), "MatrixContainer::ParseFragmentMatrix: **Matrix[]");
    159 
    160                 // parse in each entry for this matrix
    161                 input.clear();
    162                 input.seekg(ios::beg);
    163                 for (int m=skiplines+1;m--;)
    164                         input.getline(Header, 1023);            // skip header
    165                 line.str(Header);
    166                 for(int k=skipcolumns;k--;)     // skip columns in header too
    167                         line >> filename;
    168                 strncpy(Header, line.str().c_str(), 1023);
    169                 for(int j=0;j<RowCounter[MatrixNr];j++) {
    170                         Matrix[MatrixNr][j] = (double *) Malloc(sizeof(double)*ColumnCounter, "MatrixContainer::ParseFragmentMatrix: *Matrix[][]");
    171                         input.getline(filename, 1023);
    172                         stringstream lines(filename);
    173                         //cout << "Matrix at level " << j << ":";// << filename << endl;
    174                         for(int k=skipcolumns;k--;)
    175                                 lines >> filename;
    176                         for(int k=0;(k<ColumnCounter) && (!lines.eof());k++) {
    177                                 lines >> Matrix[MatrixNr][j][k];
    178                                 //cout << " " << setprecision(2) << Matrix[MatrixNr][j][k];;
    179                         }
    180                         //cout << endl;
    181                         Matrix[MatrixNr][ RowCounter[MatrixNr] ] = (double *) Malloc(sizeof(double)*ColumnCounter, "MatrixContainer::ParseFragmentMatrix: *Matrix[RowCounter[MatrixNr]][]");
    182                         for(int j=ColumnCounter;j--;)
    183                                 Matrix[MatrixNr][ RowCounter[MatrixNr] ][j] = 0.;
    184                 }
    185         } else {
    186                 cerr << "ERROR: Matrix nr. " << MatrixNr << " has column and row count of (" << ColumnCounter << "," << RowCounter[MatrixNr] << "), could not allocate nor parse!" << endl;
    187         }
    188         input.close();
    189         return true;
     151  ifstream input;
     152  stringstream line;
     153  string token;
     154  char filename[1023];
     155
     156  input.open(name, ios::in);
     157  //cout << "Opening " << name << " ... "  << input << endl;
     158  if (input == NULL) {
     159    cerr << endl << "Unable to open " << name << ", is the directory correct?" << endl;
     160    return false;
     161  }
     162
     163  // parse header
     164  Header[MatrixNr] = (char *) Malloc(sizeof(char)*1024, "MatrixContainer::ParseMatrix: *Header[]");
     165  for (int m=skiplines+1;m--;)
     166    input.getline(Header[MatrixNr], 1023);
     167 
     168  // scan header for number of columns
     169  line.str(Header[MatrixNr]);
     170  for(int k=skipcolumns;k--;)
     171    line >> Header[MatrixNr];
     172  //cout << line.str() << endl;
     173  ColumnCounter[MatrixNr]=0;
     174  while ( getline(line,token, '\t') ) {
     175    if (token.length() > 0)
     176      ColumnCounter[MatrixNr]++;
     177  }
     178  //cout << line.str() << endl;
     179  //cout << "ColumnCounter[" << MatrixNr << "]: " << ColumnCounter[MatrixNr] << "." << endl;
     180  if (ColumnCounter[MatrixNr] == 0)
     181    cerr << "ColumnCounter[" << MatrixNr << "]: " << ColumnCounter[MatrixNr] << " from file " << name << ", this is probably an error!" << endl;
     182 
     183  // scan rest for number of rows/lines
     184  RowCounter[MatrixNr]=-1;    // counts one line too much
     185  while (!input.eof()) {
     186    input.getline(filename, 1023);
     187    //cout << "Comparing: " << strncmp(filename,"MeanForce",9) << endl;
     188    RowCounter[MatrixNr]++; // then line was not last MeanForce
     189    if (strncmp(filename,"MeanForce", 9) == 0) {
     190      break;
     191    }
     192  }
     193  //cout << "RowCounter[" << MatrixNr << "]: " << RowCounter[MatrixNr] << " from file " << name << "." << endl;
     194  if (RowCounter[MatrixNr] == 0)
     195    cerr << "RowCounter[" << MatrixNr << "]: " << RowCounter[MatrixNr] << " from file " << name << ", this is probably an error!" << endl;
     196 
     197  // allocate matrix if it's not zero dimension in one direction
     198  if ((ColumnCounter[MatrixNr] > 0) && (RowCounter[MatrixNr] > -1)) {
     199    Matrix[MatrixNr] = (double **) Malloc(sizeof(double *)*(RowCounter[MatrixNr]+1), "MatrixContainer::ParseMatrix: **Matrix[]");
     200 
     201    // parse in each entry for this matrix
     202    input.clear();
     203    input.seekg(ios::beg);
     204    for (int m=skiplines+1;m--;)
     205      input.getline(Header[MatrixNr], 1023);    // skip header
     206    line.str(Header[MatrixNr]);
     207    for(int k=skipcolumns;k--;)  // skip columns in header too
     208      line >> filename;
     209    strncpy(Header[MatrixNr], line.str().c_str(), 1023); 
     210    for(int j=0;j<RowCounter[MatrixNr];j++) {
     211      Matrix[MatrixNr][j] = (double *) Malloc(sizeof(double)*ColumnCounter[MatrixNr], "MatrixContainer::ParseMatrix: *Matrix[][]");
     212      input.getline(filename, 1023);
     213      stringstream lines(filename);
     214      //cout << "Matrix at level " << j << ":";// << filename << endl;
     215      for(int k=skipcolumns;k--;)
     216        lines >> filename;
     217      for(int k=0;(k<ColumnCounter[MatrixNr]) && (!lines.eof());k++) {
     218        lines >> Matrix[MatrixNr][j][k];
     219        //cout << " " << setprecision(2) << Matrix[MatrixNr][j][k];;
     220      }
     221      //cout << endl;
     222      Matrix[MatrixNr][ RowCounter[MatrixNr] ] = (double *) Malloc(sizeof(double)*ColumnCounter[MatrixNr], "MatrixContainer::ParseMatrix: *Matrix[RowCounter[MatrixNr]][]");
     223      for(int j=ColumnCounter[MatrixNr];j--;)
     224        Matrix[MatrixNr][ RowCounter[MatrixNr] ][j] = 0.;
     225    }
     226  } else {
     227    cerr << "ERROR: Matrix nr. " << MatrixNr << " has column and row count of (" << ColumnCounter[MatrixNr] << "," << RowCounter[MatrixNr] << "), could not allocate nor parse!" << endl;
     228  }
     229  input.close();
     230  return true;
    190231};
    191232
    192233/** Parsing a number of matrices.
    193234 * -# First, count the number of matrices by counting lines in KEYSETFILE
    194  * -# Then,
    195  *              -# construct the fragment number
    196  *              -# open the matrix file
    197  *              -# skip some lines (\a skiplines)
    198  *              -# scan header lines for number of columns
    199  *              -# scan lines for number of rows
    200  *              -# allocate matrix
    201  *              -# loop over found column and row counts and parse in each entry
    202  * -# Finally, allocate one additional matrix (\a MatrixCounter) containing combined or temporary values
     235 * -# Then, 
     236 *    -# construct the fragment number
     237 *    -# open the matrix file
     238 *    -# skip some lines (\a skiplines)
     239 *    -# scan header lines for number of columns
     240 *    -# scan lines for number of rows
     241 *    -# allocate matrix
     242 *    -# loop over found column and row counts and parse in each entry
     243 * -# Finally, allocate one additional matrix (\a MatrixCounter) containing combined or temporary values 
    203244 * \param *name directory with files
    204245 * \param *prefix prefix of each matrix file
    … …  
    206247 * \param skiplines number of inital lines to skip
    207248 * \param skiplines number of inital columns to skip
    208  * \return parsing successful
     249 * \return parsing successful 
    209250 */
    210251bool MatrixContainer::ParseFragmentMatrix(char *name, char *prefix, string suffix, int skiplines, int skipcolumns)
    211252{
    212         char filename[1023];
    213         ifstream input;
    214         char *FragmentNumber = NULL;
    215         stringstream file;
    216         string token;
    217 
    218         // count the number of matrices
    219         MatrixCounter = -1; // we count one too much
    220         file << name << FRAGMENTPREFIX << KEYSETFILE;
    221         input.open(file.str().c_str(), ios::in);
    222         if (input == NULL) {
    223                 cout << endl << "Unable to open " << file.str() << ", is the directory correct?" << endl;
    224                 return false;
    225         }
    226         while (!input.eof()) {
    227                 input.getline(filename, 1023);
    228                 stringstream zeile(filename);
    229                 MatrixCounter++;
    230         }
    231         input.close();
    232         cout << "Determined " << MatrixCounter << " fragments." << endl;
    233 
    234         cout << "Parsing through each fragment and retrieving " << prefix << suffix << "." << endl;
    235         Matrix = (double ***) ReAlloc(Matrix, sizeof(double **)*(MatrixCounter+1), "MatrixContainer::ParseFragmentMatrix: ***Matrix"); // one more each for the total molecule
    236         RowCounter = (int *) ReAlloc(RowCounter, sizeof(int)*(MatrixCounter+1), "MatrixContainer::ParseFragmentMatrix: *RowCounter");
    237         for(int i=MatrixCounter+1;i--;) {
    238                 Matrix[i] = NULL;
    239                 RowCounter[i] = -1;
    240         }
    241         for(int i=0; i < MatrixCounter;i++) {
    242                 // open matrix file
    243                 FragmentNumber = FixedDigitNumber(MatrixCounter, i);
    244                 file.str(" ");
    245                 file << name << FRAGMENTPREFIX << FragmentNumber << prefix << suffix;
    246                 if (!ParseMatrix(file.str().c_str(), skiplines, skipcolumns, i))
    247                         return false;
    248                 Free((void **)&FragmentNumber, "MatrixContainer::ParseFragmentMatrix: *FragmentNumber");
    249         }
    250         return true;
     253  char filename[1023];
     254  ifstream input;
     255  char *FragmentNumber = NULL;
     256  stringstream file;
     257  string token;
     258 
     259  // count the number of matrices
     260  MatrixCounter = -1; // we count one too much
     261  file << name << FRAGMENTPREFIX << KEYSETFILE;
     262  input.open(file.str().c_str(), ios::in);
     263  if (input == NULL) {
     264    cout << endl << "Unable to open " << file.str() << ", is the directory correct?" << endl;
     265    return false;
     266  }
     267  while (!input.eof()) {
     268    input.getline(filename, 1023);
     269    stringstream zeile(filename);
     270    MatrixCounter++;
     271  }
     272  input.close(); 
     273  cout << "Determined " << MatrixCounter << " fragments." << endl;
     274 
     275  cout << "Parsing through each fragment and retrieving " << prefix << suffix << "." << endl;
     276  Header = (char **) ReAlloc(Header, sizeof(char *)*(MatrixCounter+1), "MatrixContainer::ParseFragmentMatrix: **Header"); // one more each for the total molecule
     277  Matrix = (double ***) ReAlloc(Matrix, sizeof(double **)*(MatrixCounter+1), "MatrixContainer::ParseFragmentMatrix: ***Matrix"); // one more each for the total molecule
     278  RowCounter = (int *) ReAlloc(RowCounter, sizeof(int)*(MatrixCounter+1), "MatrixContainer::ParseFragmentMatrix: *RowCounter");
     279  ColumnCounter = (int *) ReAlloc(ColumnCounter, sizeof(int)*(MatrixCounter+1), "MatrixContainer::ParseFragmentMatrix: *ColumnCounter");
     280  for(int i=MatrixCounter+1;i--;) {
     281    Matrix[i] = NULL;
     282    Header[i] = NULL;
     283    RowCounter[i] = -1;
     284    ColumnCounter[i] = -1;
     285  }
     286  for(int i=0; i < MatrixCounter;i++) {
     287    // open matrix file
     288    FragmentNumber = FixedDigitNumber(MatrixCounter, i);
     289    file.str(" ");
     290    file << name << FRAGMENTPREFIX << FragmentNumber << prefix << suffix;
     291    if (!ParseMatrix(file.str().c_str(), skiplines, skipcolumns, i))
     292      return false;
     293    Free((void **)&FragmentNumber, "MatrixContainer::ParseFragmentMatrix: *FragmentNumber");
     294  }
     295  return true;
    251296};
    252297
    253298/** Allocates and resets the memory for a number \a MCounter of matrices.
    254  * \param *GivenHeader Header line
     299 * \param **GivenHeader Header line for each matrix
    255300 * \param MCounter number of matrices
    256301 * \param *RCounter number of rows for each matrix
    257  * \param CCounter number of columns for all matrices
    258  * \return Allocation successful
    259  */
    260 bool MatrixContainer::AllocateMatrix(char *GivenHeader, int MCounter, int *RCounter, int CCounter)
    261 {
    262         Header = (char *) Malloc(sizeof(char)*1024, "MatrixContainer::ParseFragmentMatrix: *EnergyHeader");
    263         strncpy(Header, GivenHeader, 1023);
    264 
    265         MatrixCounter = MCounter;
    266         ColumnCounter = CCounter;
    267         Matrix = (double ***) Malloc(sizeof(double **)*(MatrixCounter+1), "MatrixContainer::ParseFragmentMatrix: ***Matrix"); // one more each for the total molecule
    268         RowCounter = (int *) Malloc(sizeof(int)*(MatrixCounter+1), "MatrixContainer::ParseFragmentMatrix: *RowCounter");
    269         for(int i=MatrixCounter+1;i--;) {
    270                 RowCounter[i] = RCounter[i];
    271                 Matrix[i] = (double **) Malloc(sizeof(double *)*(RowCounter[i]+1), "MatrixContainer::ParseFragmentMatrix: **Matrix[]");
    272                 for(int j=RowCounter[i]+1;j--;) {
    273                         Matrix[i][j] = (double *) Malloc(sizeof(double)*ColumnCounter, "MatrixContainer::ParseFragmentMatrix: *Matrix[][]");
    274                         for(int k=ColumnCounter;k--;)
    275                                 Matrix[i][j][k] = 0.;
    276                 }
    277         }
    278         return true;
     302 * \param *CCounter number of columns for each matrix
     303 * \return Allocation successful
     304 */
     305bool MatrixContainer::AllocateMatrix(char **GivenHeader, int MCounter, int *RCounter, int *CCounter)
     306{
     307  MatrixCounter = MCounter;
     308  Header = (char **) Malloc(sizeof(char *)*(MatrixCounter+1), "MatrixContainer::AllocateMatrix: *Header");
     309  Matrix = (double ***) Malloc(sizeof(double **)*(MatrixCounter+1), "MatrixContainer::AllocateMatrix: ***Matrix"); // one more each for the total molecule
     310  RowCounter = (int *) Malloc(sizeof(int)*(MatrixCounter+1), "MatrixContainer::AllocateMatrix: *RowCounter");
     311  ColumnCounter = (int *) Malloc(sizeof(int)*(MatrixCounter+1), "MatrixContainer::AllocateMatrix: *ColumnCounter");
     312  for(int i=MatrixCounter+1;i--;) {
     313    Header[i] = (char *) Malloc(sizeof(char)*1024, "MatrixContainer::AllocateMatrix: *Header[i]");
     314    strncpy(Header[i], GivenHeader[i], 1023);
     315    RowCounter[i] = RCounter[i];
     316    ColumnCounter[i] = CCounter[i];
     317    Matrix[i] = (double **) Malloc(sizeof(double *)*(RowCounter[i]+1), "MatrixContainer::AllocateMatrix: **Matrix[]"); 
     318    for(int j=RowCounter[i]+1;j--;) {
     319      Matrix[i][j] = (double *) Malloc(sizeof(double)*ColumnCounter[i], "MatrixContainer::AllocateMatrix: *Matrix[][]");
     320      for(int k=ColumnCounter[i];k--;)
     321        Matrix[i][j][k] = 0.;
     322    }
     323  }
     324  return true;
    279325};
    280326
    … …  
    284330bool MatrixContainer::ResetMatrix()
    285331{
    286         for(int i=MatrixCounter+1;i--;)
    287                 for(int j=RowCounter[i]+1;j--;)
    288                         for(int k=ColumnCounter;k--;)
    289                                 Matrix[i][j][k] = 0.;
    290         return true;
     332  for(int i=MatrixCounter+1;i--;)
     333    for(int j=RowCounter[i]+1;j--;)
     334      for(int k=ColumnCounter[i];k--;)
     335        Matrix[i][j][k] = 0.;
     336  return true;
    291337};
    292338
    … …  
    296342double MatrixContainer::FindMaxValue()
    297343{
    298         double max = Matrix[0][0][0];
    299         for(int i=MatrixCounter+1;i--;)
    300                 for(int j=RowCounter[i]+1;j--;)
    301                         for(int k=ColumnCounter;k--;)
    302                                 if (fabs(Matrix[i][j][k]) > max)
    303                                         max = fabs(Matrix[i][j][k]);
    304         if (fabs(max) < MYEPSILON)
    305                 max += MYEPSILON;
     344  double max = Matrix[0][0][0];
     345  for(int i=MatrixCounter+1;i--;)
     346    for(int j=RowCounter[i]+1;j--;)
     347      for(int k=ColumnCounter[i];k--;)
     348        if (fabs(Matrix[i][j][k]) > max)
     349          max = fabs(Matrix[i][j][k]);
     350  if (fabs(max) < MYEPSILON)
     351    max += MYEPSILON;
    306352 return max;
    307353};
    … …  
    312358double MatrixContainer::FindMinValue()
    313359{
    314         double min = Matrix[0][0][0];
    315         for(int i=MatrixCounter+1;i--;)
    316                 for(int j=RowCounter[i]+1;j--;)
    317                         for(int k=ColumnCounter;k--;)
    318                                 if (fabs(Matrix[i][j][k]) < min)
    319                                         min = fabs(Matrix[i][j][k]);
    320         if (fabs(min) < MYEPSILON)
    321                 min += MYEPSILON;
    322         return min;
     360  double min = Matrix[0][0][0];
     361  for(int i=MatrixCounter+1;i--;)
     362    for(int j=RowCounter[i]+1;j--;)
     363      for(int k=ColumnCounter[i];k--;)
     364        if (fabs(Matrix[i][j][k]) < min)
     365          min = fabs(Matrix[i][j][k]);
     366  if (fabs(min) < MYEPSILON)
     367    min += MYEPSILON;
     368  return min;
    323369};
    324370
    … …  
    330376bool MatrixContainer::SetLastMatrix(double value, int skipcolumns)
    331377{
    332         for(int j=RowCounter[MatrixCounter]+1;j--;)
    333                 for(int k=skipcolumns;k<ColumnCounter;k++)
    334                         Matrix[MatrixCounter][j][k] = value;
    335         return true;
     378  for(int j=RowCounter[MatrixCounter]+1;j--;)
     379    for(int k=skipcolumns;k<ColumnCounter[MatrixCounter];k++)
     380      Matrix[MatrixCounter][j][k] = value;
     381  return true;
    336382};
    337383
    … …  
    343389bool MatrixContainer::SetLastMatrix(double **values, int skipcolumns)
    344390{
    345         for(int j=RowCounter[MatrixCounter]+1;j--;)
    346                 for(int k=skipcolumns;k<ColumnCounter;k++)
    347                         Matrix[MatrixCounter][j][k] = values[j][k];
    348         return true;
    349 };
    350 
    351 /** Sums the energy with each factor and put into last element of \a ***Energies.
     391  for(int j=RowCounter[MatrixCounter]+1;j--;)
     392    for(int k=skipcolumns;k<ColumnCounter[MatrixCounter];k++)
     393      Matrix[MatrixCounter][j][k] = values[j][k];
     394  return true;
     395};
     396
     397/** Sums the entries with each factor and put into last element of \a ***Matrix.
    352398 * Sums over "E"-terms to create the "F"-terms
    353  * \param Matrix MatrixContainer with matrices (LevelCounter by ColumnCounter) with all the energies.
     399 * \param Matrix MatrixContainer with matrices (LevelCounter by *ColumnCounter) with all the energies.
    354400 * \param KeySet KeySetContainer with bond Order and association mapping of each fragment to an order
    355401 * \param Order bond order
    … …  
    358404bool MatrixContainer::SumSubManyBodyTerms(class MatrixContainer &MatrixValues, class KeySetsContainer &KeySet, int Order)
    359405{
    360         // go through each order
    361         for (int CurrentFragment=0;CurrentFragment<KeySet.FragmentsPerOrder[Order];CurrentFragment++) {
    362                 //cout << "Current Fragment is " << CurrentFragment << "/" << KeySet.OrderSet[Order][CurrentFragment] << "." << endl;
    363                 // then go per order through each suborder and pick together all the terms that contain this fragment
    364                 for(int SubOrder=0;SubOrder<=Order;SubOrder++) { // go through all suborders up to the desired order
    365                         for (int j=0;j<KeySet.FragmentsPerOrder[SubOrder];j++) { // go through all possible fragments of size suborder
    366                                 if (KeySet.Contains(KeySet.OrderSet[Order][CurrentFragment], KeySet.OrderSet[SubOrder][j])) {
    367                                         //cout << "Current other fragment is " << j << "/" << KeySet.OrderSet[SubOrder][j] << "." << endl;
    368                                         // if the fragment's indices are all in the current fragment
    369                                         for(int k=0;k<RowCounter[ KeySet.OrderSet[SubOrder][j] ];k++) { // go through all atoms in this fragment
    370                                                 int m = MatrixValues.Indices[ KeySet.OrderSet[SubOrder][j] ][k];
    371                                                 //cout << "Current index is " << k << "/" << m << "." << endl;
    372                                                 if (m != -1) { // if it's not an added hydrogen
    373                                                         for (int l=0;l<RowCounter[ KeySet.OrderSet[Order][CurrentFragment] ];l++) { // look for the corresponding index in the current fragment
    374                                                                 //cout << "Comparing " << m << " with " << MatrixValues.Indices[ KeySet.OrderSet[Order][CurrentFragment] ][l] << "." << endl;
    375                                                                 if (m == MatrixValues.Indices[ KeySet.OrderSet[Order][CurrentFragment] ][l]) {
    376                                                                         m = l;
    377                                                                         break;
    378                                                                 }
    379                                                         }
    380                                                         //cout << "Corresponding index in CurrentFragment is " << m << "." << endl;
    381                                                         if (m > RowCounter[ KeySet.OrderSet[Order][CurrentFragment] ]) {
    382                                                                 cerr << "In fragment No. " << KeySet.OrderSet[Order][CurrentFragment]    << " current force index " << m << " is greater than " << RowCounter[ KeySet.OrderSet[Order][CurrentFragment] ] << "!" << endl;
    383                                                                 return false;
    384                                                         }
    385                                                         if (Order == SubOrder) { // equal order is always copy from Energies
    386                                                                 for(int l=ColumnCounter;l--;) // then adds/subtract each column
    387                                                                         Matrix[ KeySet.OrderSet[Order][CurrentFragment] ][m][l] += MatrixValues.Matrix[ KeySet.OrderSet[SubOrder][j] ][k][l];
    388                                                         } else {
    389                                                                 for(int l=ColumnCounter;l--;)
    390                                                                         Matrix[ KeySet.OrderSet[Order][CurrentFragment] ][m][l] -= Matrix[ KeySet.OrderSet[SubOrder][j] ][k][l];
    391                                                         }
    392                                                 }
    393                                                 //if ((ColumnCounter>1) && (RowCounter[0]-1 >= 1))
    394                                                  //cout << "Fragments[ KeySet.OrderSet[" << Order << "][" << CurrentFragment << "]=" << KeySet.OrderSet[Order][CurrentFragment] << " ][" << RowCounter[0]-1 << "][" << 1 << "] = " <<   Matrix[ KeySet.OrderSet[Order][CurrentFragment] ][RowCounter[0]-1][1] << endl;
    395                                         }
    396                                 } else {
    397                                         //cout << "Fragment " << KeySet.OrderSet[SubOrder][j] << " is not contained in fragment " << KeySet.OrderSet[Order][CurrentFragment] << "." << endl;
    398                                 }
    399                         }
    400                 }
    401          //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;
    402         }
    403 
    404         return true;
     406  // go through each order
     407  for (int CurrentFragment=0;CurrentFragment<KeySet.FragmentsPerOrder[Order];CurrentFragment++) {
     408    //cout << "Current Fragment is " << CurrentFragment << "/" << KeySet.OrderSet[Order][CurrentFragment] << "." << endl;
     409    // then go per order through each suborder and pick together all the terms that contain this fragment
     410    for(int SubOrder=0;SubOrder<=Order;SubOrder++) { // go through all suborders up to the desired order
     411      for (int j=0;j<KeySet.FragmentsPerOrder[SubOrder];j++) { // go through all possible fragments of size suborder
     412        if (KeySet.Contains(KeySet.OrderSet[Order][CurrentFragment], KeySet.OrderSet[SubOrder][j])) {
     413          //cout << "Current other fragment is " << j << "/" << KeySet.OrderSet[SubOrder][j] << "." << endl;
     414          // if the fragment's indices are all in the current fragment
     415          for(int k=0;k<RowCounter[ KeySet.OrderSet[SubOrder][j] ];k++) { // go through all atoms in this fragment
     416            int m = MatrixValues.Indices[ KeySet.OrderSet[SubOrder][j] ][k];
     417            //cout << "Current index is " << k << "/" << m << "." << endl;
     418            if (m != -1) { // if it's not an added hydrogen
     419              for (int l=0;l<RowCounter[ KeySet.OrderSet[Order][CurrentFragment] ];l++) { // look for the corresponding index in the current fragment
     420                //cout << "Comparing " << m << " with " << MatrixValues.Indices[ KeySet.OrderSet[Order][CurrentFragment] ][l] << "." << endl;
     421                if (m == MatrixValues.Indices[ KeySet.OrderSet[Order][CurrentFragment] ][l]) {
     422                  m = l;
     423                  break; 
     424                }
     425              }
     426              //cout << "Corresponding index in CurrentFragment is " << m << "." << endl;
     427              if (m > RowCounter[ KeySet.OrderSet[Order][CurrentFragment] ]) {
     428                cerr << "In fragment No. " << KeySet.OrderSet[Order][CurrentFragment]  << " current force index " << m << " is greater than " << RowCounter[ KeySet.OrderSet[Order][CurrentFragment] ] << "!" << endl;
     429                return false;
     430              }
     431              if (Order == SubOrder) { // equal order is always copy from Energies
     432                for(int l=ColumnCounter[ KeySet.OrderSet[SubOrder][j] ];l--;) // then adds/subtract each column
     433                  Matrix[ KeySet.OrderSet[Order][CurrentFragment] ][m][l] += MatrixValues.Matrix[ KeySet.OrderSet[SubOrder][j] ][k][l];
     434              } else {
     435                for(int l=ColumnCounter[ KeySet.OrderSet[SubOrder][j] ];l--;)
     436                  Matrix[ KeySet.OrderSet[Order][CurrentFragment] ][m][l] -= Matrix[ KeySet.OrderSet[SubOrder][j] ][k][l];
     437              }
     438            }
     439            //if ((ColumnCounter[ KeySet.OrderSet[SubOrder][j] ]>1) && (RowCounter[0]-1 >= 1))
     440             //cout << "Fragments[ KeySet.OrderSet[" << Order << "][" << CurrentFragment << "]=" << KeySet.OrderSet[Order][CurrentFragment] << " ][" << RowCounter[0]-1 << "][" << 1 << "] = " <<  Matrix[ KeySet.OrderSet[Order][CurrentFragment] ][RowCounter[0]-1][1] << endl;
     441          }
     442        } else {
     443          //cout << "Fragment " << KeySet.OrderSet[SubOrder][j] << " is not contained in fragment " << KeySet.OrderSet[Order][CurrentFragment] << "." << endl;
     444        }
     445      }
     446    }
     447   //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;
     448  }
     449 
     450  return true;
    405451};
    406452
    … …  
    412458bool MatrixContainer::WriteTotalFragments(const char *name, const char *prefix)
    413459{
    414         ofstream output;
    415         char *FragmentNumber = NULL;
    416 
    417         cout << "Writing fragment files." << endl;
    418         for(int i=0;i<MatrixCounter;i++) {
    419                 stringstream line;
    420                 FragmentNumber = FixedDigitNumber(MatrixCounter, i);
    421                 line << name << FRAGMENTPREFIX << FragmentNumber << "/" << prefix;
    422                 Free((void **)&FragmentNumber, "*FragmentNumber");
    423                 output.open(line.str().c_str(), ios::out);
    424                 if (output == NULL) {
    425                         cerr << "Unable to open output energy file " << line.str() << "!" << endl;
    426                         return false;
    427                 }
    428                 output << Header << endl;
    429                 for(int j=0;j<RowCounter[i];j++) {
    430                         for(int k=0;k<ColumnCounter;k++)
    431                                 output << scientific << Matrix[i][j][k] << "\t";
    432                         output << endl;
    433                 }
    434                 output.close();
    435         }
    436         return true;
     460  ofstream output;
     461  char *FragmentNumber = NULL;
     462
     463  cout << "Writing fragment files." << endl;
     464  for(int i=0;i<MatrixCounter;i++) {
     465    stringstream line;
     466    FragmentNumber = FixedDigitNumber(MatrixCounter, i);
     467    line << name << FRAGMENTPREFIX << FragmentNumber << "/" << prefix;
     468    Free((void **)&FragmentNumber, "*FragmentNumber");
     469    output.open(line.str().c_str(), ios::out);
     470    if (output == NULL) {
     471      cerr << "Unable to open output energy file " << line.str() << "!" << endl;
     472      return false;
     473    }
     474    output << Header[i] << endl;
     475    for(int j=0;j<RowCounter[i];j++) {
     476      for(int k=0;k<ColumnCounter[i];k++)
     477        output << scientific << Matrix[i][j][k] << "\t";
     478      output << endl;
     479    }
     480    output.close();
     481  }
     482  return true;
    437483};
    438484
    … …  
    445491bool MatrixContainer::WriteLastMatrix(const char *name, const char *prefix, const char *suffix)
    446492{
    447         ofstream output;
    448         stringstream line;
    449 
    450         cout << "Writing matrix values of " << suffix << "." << endl;
    451         line << name << prefix << suffix;
    452         output.open(line.str().c_str(), ios::out);
    453         if (output == NULL) {
    454                 cerr << "Unable to open output matrix file " << line.str() << "!" << endl;
    455                 return false;
    456         }
    457         output << Header << endl;
    458         for(int j=0;j<RowCounter[MatrixCounter];j++) {
    459                 for(int k=0;k<ColumnCounter;k++)
    460                         output << scientific << Matrix[MatrixCounter][j][k] << "\t";
    461                 output << endl;
    462         }
    463         output.close();
    464         return true;
     493  ofstream output;
     494  stringstream line;
     495
     496  cout << "Writing matrix values of " << suffix << "." << endl;
     497  line << name << prefix << suffix;
     498  output.open(line.str().c_str(), ios::out);
     499  if (output == NULL) {
     500    cerr << "Unable to open output matrix file " << line.str() << "!" << endl;
     501    return false;
     502  }
     503  output << Header[MatrixCounter] << endl;
     504  for(int j=0;j<RowCounter[MatrixCounter];j++) {
     505    for(int k=0;k<ColumnCounter[MatrixCounter];k++)
     506      output << scientific << Matrix[MatrixCounter][j][k] << "\t";
     507    output << endl;
     508  }
     509  output.close();
     510  return true;
    465511};
    466512
    467513// ======================================= CLASS EnergyMatrix =============================
    468 
    469 /** Create a trivial energy index mapping.
    470  * This just maps 1 to 1, 2 to 2 and so on for all fragments.
    471  * \return creation sucessful
    472  */
    473 bool EnergyMatrix::ParseIndices()
    474 {
    475         cout << "Parsing energy indices." << endl;
    476         Indices = (int **) Malloc(sizeof(int *)*(MatrixCounter+1), "EnergyMatrix::ParseIndices: **Indices");
    477         for(int i=MatrixCounter+1;i--;) {
    478                 Indices[i] = (int *) Malloc(sizeof(int)*RowCounter[i], "EnergyMatrix::ParseIndices: *Indices[]");
    479                 for(int j=RowCounter[i];j--;)
    480                         Indices[i][j] = j;
    481         }
    482         return true;
    483 };
    484514
    485515/** Sums the energy with each factor and put into last element of \a EnergyMatrix::Matrix.
    486516 * Sums over the "F"-terms in ANOVA decomposition.
    487  * \param Matrix MatrixContainer with matrices (LevelCounter by ColumnCounter) with all the energies.
     517 * \param Matrix MatrixContainer with matrices (LevelCounter by *ColumnCounter) with all the energies.
    488518 * \param CorrectionFragments MatrixContainer with hydrogen saturation correction per fragments
    489519 * \param KeySet KeySetContainer with bond Order and association mapping of each fragment to an order
    … …  
    494524bool EnergyMatrix::SumSubEnergy(class EnergyMatrix &Fragments, class EnergyMatrix *CorrectionFragments, class KeySetsContainer &KeySet, int Order, double sign)
    495525{
    496         // sum energy
    497         if (CorrectionFragments == NULL)
    498                 for(int i=KeySet.FragmentsPerOrder[Order];i--;)
    499                         for(int j=RowCounter[ KeySet.OrderSet[Order][i] ];j--;)
    500                                 for(int k=ColumnCounter;k--;)
    501                                         Matrix[MatrixCounter][j][k] += sign*Fragments.Matrix[ KeySet.OrderSet[Order][i] ][j][k];
    502         else
    503                 for(int i=KeySet.FragmentsPerOrder[Order];i--;)
    504                         for(int j=RowCounter[ KeySet.OrderSet[Order][i] ];j--;)
    505                                 for(int k=ColumnCounter;k--;)
    506                                         Matrix[MatrixCounter][j][k] += sign*(Fragments.Matrix[ KeySet.OrderSet[Order][i] ][j][k] + CorrectionFragments->Matrix[ KeySet.OrderSet[Order][i] ][j][k]);
    507         return true;
     526  // sum energy
     527  if (CorrectionFragments == NULL)
     528    for(int i=KeySet.FragmentsPerOrder[Order];i--;)
     529      for(int j=RowCounter[ KeySet.OrderSet[Order][i] ];j--;)
     530        for(int k=ColumnCounter[ KeySet.OrderSet[Order][i] ];k--;)
     531          Matrix[MatrixCounter][j][k] += sign*Fragments.Matrix[ KeySet.OrderSet[Order][i] ][j][k];
     532  else
     533    for(int i=KeySet.FragmentsPerOrder[Order];i--;)
     534      for(int j=RowCounter[ KeySet.OrderSet[Order][i] ];j--;)
     535        for(int k=ColumnCounter[ KeySet.OrderSet[Order][i] ];k--;)
     536          Matrix[MatrixCounter][j][k] += sign*(Fragments.Matrix[ KeySet.OrderSet[Order][i] ][j][k] + CorrectionFragments->Matrix[ KeySet.OrderSet[Order][i] ][j][k]);
     537  return true;
    508538};
    509539
    … …  
    515545 * \param skiplines number of inital columns to skip
    516546 * \return parsing successful
    517  */
     547 */ 
    518548bool EnergyMatrix::ParseFragmentMatrix(char *name, char *prefix, string suffix, int skiplines, int skipcolumns)
    519549{
    520         char filename[1024];
    521         bool status = MatrixContainer::ParseFragmentMatrix(name, prefix, suffix, skiplines, skipcolumns);
    522 
    523         if (status) {
    524                 // count maximum of columns
    525                 RowCounter[MatrixCounter] = 0;
    526                 for(int j=0; j < MatrixCounter;j++) // (energy matrix might be bigger than number of atoms in terms of rows)
    527                         if (RowCounter[j] > RowCounter[MatrixCounter])
    528                                 RowCounter[MatrixCounter] = RowCounter[j];
    529                 // allocate last plus one matrix
    530                 cout << "Allocating last plus one matrix with " << (RowCounter[MatrixCounter]+1) << " rows and " << ColumnCounter << " columns." << endl;
    531                 Matrix[MatrixCounter] = (double **) Malloc(sizeof(double *)*(RowCounter[MatrixCounter]+1), "MatrixContainer::ParseFragmentMatrix: **Matrix[]");
    532                 for(int j=0;j<=RowCounter[MatrixCounter];j++)
    533                         Matrix[MatrixCounter][j] = (double *) Malloc(sizeof(double)*ColumnCounter, "MatrixContainer::ParseFragmentMatrix: *Matrix[][]");
    534 
    535                 // try independently to parse global energysuffix file if present
    536                 strncpy(filename, name, 1023);
    537                 strncat(filename, prefix, 1023-strlen(filename));
    538                 strncat(filename, suffix.c_str(), 1023-strlen(filename));
    539                 ParseMatrix(filename, skiplines, skipcolumns, MatrixCounter);
    540         }
    541         return status;
     550  char filename[1024];
     551  bool status = MatrixContainer::ParseFragmentMatrix(name, prefix, suffix, skiplines, skipcolumns);
     552
     553  if (status) {
     554    // count maximum of columns
     555    RowCounter[MatrixCounter] = 0;
     556    ColumnCounter[MatrixCounter] = 0;
     557    for(int j=0; j < MatrixCounter;j++) { // (energy matrix might be bigger than number of atoms in terms of rows)
     558      if (RowCounter[j] > RowCounter[MatrixCounter])
     559        RowCounter[MatrixCounter] = RowCounter[j];
     560      if (ColumnCounter[j] > ColumnCounter[MatrixCounter])  // take maximum of all for last matrix
     561        ColumnCounter[MatrixCounter] = ColumnCounter[j];
     562    }
     563    // allocate last plus one matrix
     564    cout << "Allocating last plus one matrix with " << (RowCounter[MatrixCounter]+1) << " rows and " << ColumnCounter[MatrixCounter] << " columns." << endl;
     565    Matrix[MatrixCounter] = (double **) Malloc(sizeof(double *)*(RowCounter[MatrixCounter]+1), "MatrixContainer::ParseFragmentMatrix: **Matrix[]");
     566    for(int j=0;j<=RowCounter[MatrixCounter];j++)
     567      Matrix[MatrixCounter][j] = (double *) Malloc(sizeof(double)*ColumnCounter[MatrixCounter], "MatrixContainer::ParseFragmentMatrix: *Matrix[][]");
     568   
     569    // try independently to parse global energysuffix file if present
     570    strncpy(filename, name, 1023);
     571    strncat(filename, prefix, 1023-strlen(filename));
     572    strncat(filename, suffix.c_str(), 1023-strlen(filename));
     573    ParseMatrix(filename, skiplines, skipcolumns, MatrixCounter);
     574  }
     575  return status;
    542576};
    543577
    … …  
    546580/** Parsing force Indices of each fragment
    547581 * \param *name directory with \a ForcesFile
    548  * \return parsing successful
    549  */
    550 bool ForceMatrix::ParseIndices(char *name)
    551 {
    552         ifstream input;
    553         char *FragmentNumber = NULL;
    554         char filename[1023];
    555         stringstream line;
    556 
    557         cout << "Parsing force indices for " << MatrixCounter << " matrices." << endl;
    558         Indices = (int **) Malloc(sizeof(int *)*(MatrixCounter+1), "ForceMatrix::ParseIndices: **Indices");
    559         line << name << FRAGMENTPREFIX << FORCESFILE;
    560         input.open(line.str().c_str(), ios::in);
    561         //cout << "Opening " << line.str() << " ... "   << input << endl;
    562         if (input == NULL) {
    563                 cout << endl << "Unable to open " << line.str() << ", is the directory correct?" << endl;
    564                 return false;
    565         }
    566         for (int i=0;(i<MatrixCounter) && (!input.eof());i++) {
    567                 // get the number of atoms for this fragment
    568                 input.getline(filename, 1023);
    569                 line.str(filename);
    570                 // parse the values
    571                 Indices[i] = (int *) Malloc(sizeof(int)*RowCounter[i], "ForceMatrix::ParseIndices: *Indices[]");
    572                 FragmentNumber = FixedDigitNumber(MatrixCounter, i);
    573                 //cout << FRAGMENTPREFIX << FragmentNumber << "[" << RowCounter[i] << "]:";
    574                 Free((void **)&FragmentNumber, "ForceMatrix::ParseIndices: *FragmentNumber");
    575                 for(int j=0;(j<RowCounter[i]) && (!line.eof());j++) {
    576                         line >> Indices[i][j];
    577                         //cout << " " << Indices[i][j];
    578                 }
    579                 //cout << endl;
    580         }
    581         Indices[MatrixCounter] = (int *) Malloc(sizeof(int)*RowCounter[MatrixCounter], "ForceMatrix::ParseIndices: *Indices[]");
    582         for(int j=RowCounter[MatrixCounter];j--;) {
    583                 Indices[MatrixCounter][j] = j;
    584         }
    585         input.close();
    586         return true;
     582 * \return parsing successful 
     583 */
     584bool ForceMatrix::ParseIndices(char *name) 
     585{
     586  ifstream input;
     587  char *FragmentNumber = NULL;
     588  char filename[1023];
     589  stringstream line;
     590 
     591  cout << "Parsing force indices for " << MatrixCounter << " matrices." << endl;
     592  Indices = (int **) Malloc(sizeof(int *)*(MatrixCounter+1), "ForceMatrix::ParseIndices: **Indices");
     593  line << name << FRAGMENTPREFIX << FORCESFILE;
     594  input.open(line.str().c_str(), ios::in);
     595  //cout << "Opening " << line.str() << " ... "  << input << endl;
     596  if (input == NULL) {
     597    cout << endl << "Unable to open " << line.str() << ", is the directory correct?" << endl;
     598    return false;
     599  }
     600  for (int i=0;(i<MatrixCounter) && (!input.eof());i++) {
     601    // get the number of atoms for this fragment
     602    input.getline(filename, 1023);
     603    line.str(filename);
     604    // parse the values
     605    Indices[i] = (int *) Malloc(sizeof(int)*RowCounter[i], "ForceMatrix::ParseIndices: *Indices[]");
     606    FragmentNumber = FixedDigitNumber(MatrixCounter, i);
     607    //cout << FRAGMENTPREFIX << FragmentNumber << "[" << RowCounter[i] << "]:";
     608    Free((void **)&FragmentNumber, "ForceMatrix::ParseIndices: *FragmentNumber");
     609    for(int j=0;(j<RowCounter[i]) && (!line.eof());j++) {
     610      line >> Indices[i][j];
     611      //cout << " " << Indices[i][j];
     612    }
     613    //cout << endl;
     614  }
     615  Indices[MatrixCounter] = (int *) Malloc(sizeof(int)*RowCounter[MatrixCounter], "ForceMatrix::ParseIndices: *Indices[]");
     616  for(int j=RowCounter[MatrixCounter];j--;) {
     617    Indices[MatrixCounter][j] = j;
     618  }
     619  input.close();
     620  return true;
    587621};
    588622
    589623
    590624/** Sums the forces and puts into last element of \a ForceMatrix::Matrix.
    591  * \param Matrix MatrixContainer with matrices (LevelCounter by ColumnCounter) with all the energies.
     625 * \param Matrix MatrixContainer with matrices (LevelCounter by *ColumnCounter) with all the energies.
    592626 * \param KeySet KeySetContainer with bond Order and association mapping of each fragment to an order
    593627 * \param Order bond order
    594  *      \param sign +1 or -1
     628 *  \param sign +1 or -1
    595629 * \return true if summing was successful
    596630 */
    597631bool ForceMatrix::SumSubForces(class ForceMatrix &Fragments, class KeySetsContainer &KeySet, int Order, double sign)
    598632{
    599         int FragmentNr;
    600         // sum forces
    601         for(int i=0;i<KeySet.FragmentsPerOrder[Order];i++) {
    602                 FragmentNr = KeySet.OrderSet[Order][i];
    603                 for(int l=0;l<RowCounter[ FragmentNr ];l++) {
    604                         int j = Indices[ FragmentNr ][l];
    605                         if (j > RowCounter[MatrixCounter]) {
    606                                 cerr << "Current force index " << j << " is greater than " << RowCounter[MatrixCounter] << "!" << endl;
    607                                 return false;
    608                         }
    609                         if (j != -1) {
    610                                 //if (j == 0) cout << "Summing onto ion 0, type 0 from fragment " << FragmentNr << ", ion " << l << "." << endl;
    611                                 for(int k=2;k<ColumnCounter;k++)
    612                                         Matrix[MatrixCounter][j][k] += sign*Fragments.Matrix[ FragmentNr ][l][k];
    613                         }
    614                 }
    615         }
    616         return true;
     633  int FragmentNr;
     634  // sum forces
     635  for(int i=0;i<KeySet.FragmentsPerOrder[Order];i++) {
     636    FragmentNr = KeySet.OrderSet[Order][i];
     637    for(int l=0;l<RowCounter[ FragmentNr ];l++) {
     638      int j = Indices[ FragmentNr ][l];
     639      if (j > RowCounter[MatrixCounter]) {
     640        cerr << "Current force index " << j << " is greater than " << RowCounter[MatrixCounter] << "!" << endl;
     641        return false;
     642      }
     643      if (j != -1) {
     644        //if (j == 0) cout << "Summing onto ion 0, type 0 from fragment " << FragmentNr << ", ion " << l << "." << endl;
     645        for(int k=2;k<ColumnCounter[MatrixCounter];k++)
     646          Matrix[MatrixCounter][j][k] += sign*Fragments.Matrix[ FragmentNr ][l][k];
     647      }
     648    }
     649  }
     650  return true;
    617651};
    618652
    … …  
    625659 * \param skiplines number of inital columns to skip
    626660 * \return parsing successful
    627  */
     661 */ 
    628662bool ForceMatrix::ParseFragmentMatrix(char *name, char *prefix, string suffix, int skiplines, int skipcolumns)
    629663{
    630         char filename[1023];
    631         ifstream input;
    632         stringstream file;
    633         int nr;
    634         bool status = MatrixContainer::ParseFragmentMatrix(name, prefix, suffix, skiplines, skipcolumns);
    635 
    636         if (status) {
    637                 // count number of atoms for last plus one matrix
    638                 file << name << FRAGMENTPREFIX << KEYSETFILE;
    639                 input.open(file.str().c_str(), ios::in);
    640                 if (input == NULL) {
    641                         cout << endl << "Unable to open " << file.str() << ", is the directory correct?" << endl;
    642                         return false;
    643                 }
    644                 RowCounter[MatrixCounter] = 0;
    645                 while (!input.eof()) {
    646                         input.getline(filename, 1023);
    647                         stringstream zeile(filename);
    648                         while (!zeile.eof()) {
    649                                 zeile >> nr;
    650                                 //cout << "Current index: " << nr << "." << endl;
    651                                 if (nr > RowCounter[MatrixCounter])
    652                                         RowCounter[MatrixCounter] = nr;
    653                         }
    654                 }
    655                 RowCounter[MatrixCounter]++;            // nr start at 0, count starts at 1
    656                 input.close();
    657 
    658                 // allocate last plus one matrix
    659                 cout << "Allocating last plus one matrix with " << (RowCounter[MatrixCounter]+1) << " rows and " << ColumnCounter << " columns." << endl;
    660                 Matrix[MatrixCounter] = (double **) Malloc(sizeof(double *)*(RowCounter[MatrixCounter]+1), "MatrixContainer::ParseFragmentMatrix: **Matrix[]");
    661                 for(int j=0;j<=RowCounter[MatrixCounter];j++)
    662                         Matrix[MatrixCounter][j] = (double *) Malloc(sizeof(double)*ColumnCounter, "MatrixContainer::ParseFragmentMatrix: *Matrix[][]");
    663 
    664                 // try independently to parse global forcesuffix file if present
    665                 strncpy(filename, name, 1023);
    666                 strncat(filename, prefix, 1023-strlen(filename));
    667                 strncat(filename, suffix.c_str(), 1023-strlen(filename));
    668                 ParseMatrix(filename, skiplines, skipcolumns, MatrixCounter);
    669         }
    670 
    671 
    672         return status;
     664  char filename[1023];
     665  ifstream input;
     666  stringstream file;
     667  int nr;
     668  bool status = MatrixContainer::ParseFragmentMatrix(name, prefix, suffix, skiplines, skipcolumns);
     669
     670  if (status) {
     671    // count number of atoms for last plus one matrix
     672    file << name << FRAGMENTPREFIX << KEYSETFILE;
     673    input.open(file.str().c_str(), ios::in);
     674    if (input == NULL) {
     675      cout << endl << "Unable to open " << file.str() << ", is the directory correct?" << endl;
     676      return false;
     677    }
     678    RowCounter[MatrixCounter] = 0;
     679    while (!input.eof()) {
     680      input.getline(filename, 1023);
     681      stringstream zeile(filename);
     682      while (!zeile.eof()) {
     683        zeile >> nr;
     684        //cout << "Current index: " << nr << "." << endl;
     685        if (nr > RowCounter[MatrixCounter])
     686          RowCounter[MatrixCounter] = nr;
     687      }
     688    }
     689    RowCounter[MatrixCounter]++;    // nr start at 0, count starts at 1
     690    input.close();
     691
     692    ColumnCounter[MatrixCounter] = 0;
     693    for(int j=0; j < MatrixCounter;j++) { // (energy matrix might be bigger than number of atoms in terms of rows)
     694      if (ColumnCounter[j] > ColumnCounter[MatrixCounter])  // take maximum of all for last matrix
     695        ColumnCounter[MatrixCounter] = ColumnCounter[j];
     696    }
     697 
     698    // allocate last plus one matrix
     699    cout << "Allocating last plus one matrix with " << (RowCounter[MatrixCounter]+1) << " rows and " << ColumnCounter[MatrixCounter] << " columns." << endl;
     700    Matrix[MatrixCounter] = (double **) Malloc(sizeof(double *)*(RowCounter[MatrixCounter]+1), "MatrixContainer::ParseFragmentMatrix: **Matrix[]");
     701    for(int j=0;j<=RowCounter[MatrixCounter];j++)
     702      Matrix[MatrixCounter][j] = (double *) Malloc(sizeof(double)*ColumnCounter[MatrixCounter], "MatrixContainer::ParseFragmentMatrix: *Matrix[][]");
     703
     704    // try independently to parse global forcesuffix file if present
     705    strncpy(filename, name, 1023);
     706    strncat(filename, prefix, 1023-strlen(filename));
     707    strncat(filename, suffix.c_str(), 1023-strlen(filename));
     708    ParseMatrix(filename, skiplines, skipcolumns, MatrixCounter);
     709  }
     710 
     711
     712  return status;
     713};
     714
     715// ======================================= CLASS HessianMatrix =============================
     716
     717/** Parsing force Indices of each fragment
     718 * \param *name directory with \a ForcesFile
     719 * \return parsing successful
     720 */
     721bool HessianMatrix::ParseIndices(char *name)
     722{
     723  ifstream input;
     724  char *FragmentNumber = NULL;
     725  char filename[1023];
     726  stringstream line;
     727 
     728  cout << "Parsing hessian indices for " << MatrixCounter << " matrices." << endl;
     729  Indices = (int **) Malloc(sizeof(int *)*(MatrixCounter+1), "HessianMatrix::ParseIndices: **Indices");
     730  line << name << FRAGMENTPREFIX << FORCESFILE;
     731  input.open(line.str().c_str(), ios::in);
     732  //cout << "Opening " << line.str() << " ... "  << input << endl;
     733  if (input == NULL) {
     734    cout << endl << "Unable to open " << line.str() << ", is the directory correct?" << endl;
     735    return false;
     736  }
     737  for (int i=0;(i<MatrixCounter) && (!input.eof());i++) {
     738    // get the number of atoms for this fragment
     739    input.getline(filename, 1023);
     740    line.str(filename);
     741    // parse the values
     742    Indices[i] = (int *) Malloc(sizeof(int)*RowCounter[i], "HessianMatrix::ParseIndices: *Indices[]");
     743    FragmentNumber = FixedDigitNumber(MatrixCounter, i);
     744    //cout << FRAGMENTPREFIX << FragmentNumber << "[" << RowCounter[i] << "]:";
     745    Free((void **)&FragmentNumber, "HessianMatrix::ParseIndices: *FragmentNumber");
     746    for(int j=0;(j<RowCounter[i]) && (!line.eof());j++) {
     747      line >> Indices[i][j];
     748      //cout << " " << Indices[i][j];
     749    }
     750    //cout << endl;
     751  }
     752  Indices[MatrixCounter] = (int *) Malloc(sizeof(int)*RowCounter[MatrixCounter], "HessianMatrix::ParseIndices: *Indices[]");
     753  for(int j=RowCounter[MatrixCounter];j--;) {
     754    Indices[MatrixCounter][j] = j;
     755  }
     756  input.close();
     757  return true;
     758};
     759
     760
     761/** Sums the hessian entries and puts into last element of \a HessianMatrix::Matrix.
     762 * \param Matrix MatrixContainer with matrices (LevelCounter by *ColumnCounter) with all the energies.
     763 * \param KeySet KeySetContainer with bond Order and association mapping of each fragment to an order
     764 * \param Order bond order
     765 *  \param sign +1 or -1
     766 * \return true if summing was successful
     767 */
     768bool HessianMatrix::SumSubHessians(class HessianMatrix &Fragments, class KeySetsContainer &KeySet, int Order, double sign)
     769{
     770  int FragmentNr;
     771  // sum forces
     772  for(int i=0;i<KeySet.FragmentsPerOrder[Order];i++) {
     773    FragmentNr = KeySet.OrderSet[Order][i];
     774    for(int l=0;l<RowCounter[ FragmentNr ];l++) {
     775      int j = Indices[ FragmentNr ][l];
     776      if (j > RowCounter[MatrixCounter]) {
     777        cerr << "Current hessian index " << j << " is greater than " << RowCounter[MatrixCounter] << ", where i=" << i << ", Order=" << Order << ", l=" << l << " and FragmentNr=" << FragmentNr << "!" << endl;
     778        return false;
     779      }
     780      if (j != -1) {
     781        for(int m=0;m<ColumnCounter[ FragmentNr ];m++) {
     782          int k = Indices[ FragmentNr ][m];
     783          if (k > ColumnCounter[MatrixCounter]) {
     784            cerr << "Current hessian index " << k << " is greater than " << ColumnCounter[MatrixCounter] << ", where m=" << m << ", j=" << j << ", i=" << i << ", Order=" << Order << ", l=" << l << " and FragmentNr=" << FragmentNr << "!" << endl;
     785            return false;
     786          }
     787          if (k != -1) {
     788            //cout << "Adding " << sign*Fragments.Matrix[ FragmentNr ][l][m] << " from [" << l << "][" << m << "] onto [" << j << "][" << k << "]." << endl;
     789            Matrix[MatrixCounter][j][k] += sign*Fragments.Matrix[ FragmentNr ][l][m];
     790          }
     791        }
     792      }
     793    }
     794  }
     795  return true;
     796};
     797
     798/** Constructor for class HessianMatrix.
     799 */
     800HessianMatrix::HessianMatrix() : MatrixContainer()
     801{
     802   IsSymmetric = true;
     803}
     804
     805/** Sums the hessian entries with each factor and put into last element of \a ***Matrix.
     806 * Sums over "E"-terms to create the "F"-terms
     807 * \param Matrix MatrixContainer with matrices (LevelCounter by *ColumnCounter) with all the energies.
     808 * \param KeySet KeySetContainer with bond Order and association mapping of each fragment to an order
     809 * \param Order bond order
     810 * \return true if summing was successful
     811 */
     812bool HessianMatrix::SumSubManyBodyTerms(class MatrixContainer &MatrixValues, class KeySetsContainer &KeySet, int Order)
     813{
     814  // go through each order
     815  for (int CurrentFragment=0;CurrentFragment<KeySet.FragmentsPerOrder[Order];CurrentFragment++) {
     816    //cout << "Current Fragment is " << CurrentFragment << "/" << KeySet.OrderSet[Order][CurrentFragment] << "." << endl;
     817    // then go per order through each suborder and pick together all the terms that contain this fragment
     818    for(int SubOrder=0;SubOrder<=Order;SubOrder++) { // go through all suborders up to the desired order
     819      for (int j=0;j<KeySet.FragmentsPerOrder[SubOrder];j++) { // go through all possible fragments of size suborder
     820        if (KeySet.Contains(KeySet.OrderSet[Order][CurrentFragment], KeySet.OrderSet[SubOrder][j])) {
     821          //cout << "Current other fragment is " << j << "/" << KeySet.OrderSet[SubOrder][j] << "." << endl;
     822          // if the fragment's indices are all in the current fragment
     823          for(int k=0;k<RowCounter[ KeySet.OrderSet[SubOrder][j] ];k++) { // go through all atoms in this fragment
     824            int m = MatrixValues.Indices[ KeySet.OrderSet[SubOrder][j] ][k];
     825            //cout << "Current row index is " << k << "/" << m << "." << endl;
     826            if (m != -1) { // if it's not an added hydrogen
     827              for (int l=0;l<RowCounter[ KeySet.OrderSet[Order][CurrentFragment] ];l++) { // look for the corresponding index in the current fragment
     828                //cout << "Comparing " << m << " with " << MatrixValues.Indices[ KeySet.OrderSet[Order][CurrentFragment] ][l] << "." << endl;
     829                if (m == MatrixValues.Indices[ KeySet.OrderSet[Order][CurrentFragment] ][l]) {
     830                  m = l;
     831                  break; 
     832                }
     833              }
     834              //cout << "Corresponding row index for " << k << " in CurrentFragment is " << m << "." << endl;
     835              if (m > RowCounter[ KeySet.OrderSet[Order][CurrentFragment] ]) {
     836                cerr << "In fragment No. " << KeySet.OrderSet[Order][CurrentFragment]   << " current row index " << m << " is greater than " << RowCounter[ KeySet.OrderSet[Order][CurrentFragment] ] << "!" << endl;
     837                return false;
     838              }
     839             
     840              for(int l=0;l<ColumnCounter[ KeySet.OrderSet[SubOrder][j] ];l++) {
     841                int n = MatrixValues.Indices[ KeySet.OrderSet[SubOrder][j] ][l];
     842                //cout << "Current column index is " << l << "/" << n << "." << endl;
     843                if (n != -1) { // if it's not an added hydrogen
     844                  for (int p=0;p<ColumnCounter[ KeySet.OrderSet[Order][CurrentFragment] ];p++) { // look for the corresponding index in the current fragment
     845                    //cout << "Comparing " << n << " with " << MatrixValues.Indices[ KeySet.OrderSet[Order][CurrentFragment] ][p] << "." << endl;
     846                    if (n == MatrixValues.Indices[ KeySet.OrderSet[Order][CurrentFragment] ][p]) {
     847                      n = p;
     848                      break; 
     849                    }
     850                  }
     851                  //cout << "Corresponding column index for " << l << " in CurrentFragment is " << n << "." << endl;
     852                  if (n > ColumnCounter[ KeySet.OrderSet[Order][CurrentFragment] ]) {
     853                    cerr << "In fragment No. " << KeySet.OrderSet[Order][CurrentFragment]   << " current column index " << n << " is greater than " << ColumnCounter[ KeySet.OrderSet[Order][CurrentFragment] ] << "!" << endl;
     854                    return false;
     855                  }
     856                  if (Order == SubOrder) { // equal order is always copy from Energies
     857                    //cout << "Adding " << MatrixValues.Matrix[ KeySet.OrderSet[SubOrder][j] ][k][l] << " from [" << k << "][" << l << "] onto [" << m << "][" << n << "]." << endl;
     858                    Matrix[ KeySet.OrderSet[Order][CurrentFragment] ][m][n] += MatrixValues.Matrix[ KeySet.OrderSet[SubOrder][j] ][k][l];
     859                  } else {
     860                    //cout << "Subtracting " << Matrix[ KeySet.OrderSet[SubOrder][j] ][k][l] << " from [" << k << "][" << l << "] onto [" << m << "][" << n << "]." << endl;
     861                    Matrix[ KeySet.OrderSet[Order][CurrentFragment] ][m][n] -= Matrix[ KeySet.OrderSet[SubOrder][j] ][k][l];
     862                  }
     863                }
     864              }
     865            }
     866            //if ((ColumnCounter[ KeySet.OrderSet[SubOrder][j] ]>1) && (RowCounter[0]-1 >= 1))
     867             //cout << "Fragments[ KeySet.OrderSet[" << Order << "][" << CurrentFragment << "]=" << KeySet.OrderSet[Order][CurrentFragment] << " ][" << RowCounter[0]-1 << "][" << 1 << "] = " <<  Matrix[ KeySet.OrderSet[Order][CurrentFragment] ][RowCounter[0]-1][1] << endl;
     868          }
     869        } else {
     870          //cout << "Fragment " << KeySet.OrderSet[SubOrder][j] << " is not contained in fragment " << KeySet.OrderSet[Order][CurrentFragment] << "." << endl;
     871        }
     872      }
     873    }
     874   //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;
     875  }
     876 
     877  return true;
     878};
     879
     880/** Calls MatrixContainer::ParseFragmentMatrix() and additionally allocates last plus one matrix.
     881 * \param *name directory with files
     882 * \param *prefix prefix of each matrix file
     883 * \param *suffix suffix of each matrix file
     884 * \param skiplines number of inital lines to skip
     885 * \param skiplines number of inital columns to skip
     886 * \return parsing successful
     887 */
     888bool HessianMatrix::ParseFragmentMatrix(char *name, char *prefix, string suffix, int skiplines, int skipcolumns)
     889{
     890  char filename[1023];
     891  ifstream input;
     892  stringstream file;
     893  int nr;
     894  bool status = MatrixContainer::ParseFragmentMatrix(name, prefix, suffix, skiplines, skipcolumns);
     895
     896  if (status) {
     897    // count number of atoms for last plus one matrix
     898    file << name << FRAGMENTPREFIX << KEYSETFILE;
     899    input.open(file.str().c_str(), ios::in);
     900    if (input == NULL) {
     901      cout << endl << "Unable to open " << file.str() << ", is the directory correct?" << endl;
     902      return false;
     903    }
     904    RowCounter[MatrixCounter] = 0;
     905    ColumnCounter[MatrixCounter] = 0;
     906    while (!input.eof()) {
     907      input.getline(filename, 1023);
     908      stringstream zeile(filename);
     909      while (!zeile.eof()) {
     910        zeile >> nr;
     911        //cout << "Current index: " << nr << "." << endl;
     912        if (nr > RowCounter[MatrixCounter]) {
     913          RowCounter[MatrixCounter] = nr;
     914          ColumnCounter[MatrixCounter] = nr;
     915        }
     916      }
     917    }
     918    RowCounter[MatrixCounter]++;    // nr start at 0, count starts at 1
     919    ColumnCounter[MatrixCounter]++;    // nr start at 0, count starts at 1
     920    input.close();
     921 
     922    // allocate last plus one matrix
     923    cout << "Allocating last plus one matrix with " << (RowCounter[MatrixCounter]+1) << " rows and " << ColumnCounter[MatrixCounter] << " columns." << endl;
     924    Matrix[MatrixCounter] = (double **) Malloc(sizeof(double *)*(RowCounter[MatrixCounter]+1), "MatrixContainer::ParseFragmentMatrix: **Matrix[]");
     925    for(int j=0;j<=RowCounter[MatrixCounter];j++)
     926      Matrix[MatrixCounter][j] = (double *) Malloc(sizeof(double)*ColumnCounter[MatrixCounter], "MatrixContainer::ParseFragmentMatrix: *Matrix[][]");
     927
     928    // try independently to parse global forcesuffix file if present
     929    strncpy(filename, name, 1023);
     930    strncat(filename, prefix, 1023-strlen(filename));
     931    strncat(filename, suffix.c_str(), 1023-strlen(filename));
     932    ParseMatrix(filename, skiplines, skipcolumns, MatrixCounter);
     933  }
     934 
     935
     936  return status;
    673937};
    674938
    … …  
    678942 */
    679943KeySetsContainer::KeySetsContainer() {
    680         KeySets = NULL;
    681         AtomCounter = NULL;
    682         FragmentCounter = 0;
    683         Order = 0;
    684         FragmentsPerOrder = 0;
    685         OrderSet = NULL;
     944  KeySets = NULL;
     945  AtomCounter = NULL;
     946  FragmentCounter = 0;
     947  Order = 0;
     948  FragmentsPerOrder = 0;
     949  OrderSet = NULL;
    686950};
    687951
    … …  
    689953 */
    690954KeySetsContainer::~KeySetsContainer() {
    691         for(int i=FragmentCounter;i--;)
    692                 Free((void **)&KeySets[i], "KeySetsContainer::~KeySetsContainer: *KeySets[]");
    693         for(int i=Order;i--;)
    694                 Free((void **)&OrderSet[i], "KeySetsContainer::~KeySetsContainer: *OrderSet[]");
    695         Free((void **)&KeySets, "KeySetsContainer::~KeySetsContainer: **KeySets");
    696         Free((void **)&OrderSet, "KeySetsContainer::~KeySetsContainer: **OrderSet");
    697         Free((void **)&AtomCounter, "KeySetsContainer::~KeySetsContainer: *AtomCounter");
    698         Free((void **)&FragmentsPerOrder, "KeySetsContainer::~KeySetsContainer: *FragmentsPerOrder");
     955  for(int i=FragmentCounter;i--;)
     956    Free((void **)&KeySets[i], "KeySetsContainer::~KeySetsContainer: *KeySets[]");
     957  for(int i=Order;i--;)
     958    Free((void **)&OrderSet[i], "KeySetsContainer::~KeySetsContainer: *OrderSet[]");
     959  Free((void **)&KeySets, "KeySetsContainer::~KeySetsContainer: **KeySets");
     960  Free((void **)&OrderSet, "KeySetsContainer::~KeySetsContainer: **OrderSet");
     961  Free((void **)&AtomCounter, "KeySetsContainer::~KeySetsContainer: *AtomCounter");
     962  Free((void **)&FragmentsPerOrder, "KeySetsContainer::~KeySetsContainer: *FragmentsPerOrder");
    699963};
    700964
    … …  
    706970 */
    707971bool KeySetsContainer::ParseKeySets(const char *name, const int *ACounter, const int FCounter) {
    708         ifstream input;
    709         char *FragmentNumber = NULL;
    710         stringstream file;
    711         char filename[1023];
    712 
    713         FragmentCounter = FCounter;
    714         cout << "Parsing key sets." << endl;
    715         KeySets = (int **) Malloc(sizeof(int *)*FragmentCounter, "KeySetsContainer::ParseKeySets: **KeySets");
    716         for(int i=FragmentCounter;i--;)
    717                 KeySets[i] = NULL;
    718         file << name << FRAGMENTPREFIX << KEYSETFILE;
    719         input.open(file.str().c_str(), ios::in);
    720         if (input == NULL) {
    721                 cout << endl << "Unable to open " << file.str() << ", is the directory correct?" << endl;
    722                 return false;
    723         }
    724 
    725         AtomCounter = (int *) Malloc(sizeof(int)*FragmentCounter, "KeySetsContainer::ParseKeySets: *RowCounter");
    726         for(int i=0;(i<FragmentCounter) && (!input.eof());i++) {
    727                 stringstream line;
    728                 AtomCounter[i] = ACounter[i];
    729                 // parse the values
    730                 KeySets[i] = (int *) Malloc(sizeof(int)*AtomCounter[i], "KeySetsContainer::ParseKeySets: *KeySets[]");
    731                 for(int j=AtomCounter[i];j--;)
    732                         KeySets[i][j] = -1;
    733                 FragmentNumber = FixedDigitNumber(FragmentCounter, i);
    734                 //cout << FRAGMENTPREFIX << FragmentNumber << "[" << AtomCounter[i] << "]:";
    735                 Free((void **)&FragmentNumber, "KeySetsContainer::ParseKeySets: *FragmentNumber");
    736                 input.getline(filename, 1023);
    737                 line.str(filename);
    738                 for(int j=0;(j<AtomCounter[i]) && (!line.eof());j++) {
    739                         line >> KeySets[i][j];
    740                         //cout << " " << KeySets[i][j];
    741                 }
    742                 //cout << endl;
    743         }
    744         input.close();
    745         return true;
     972  ifstream input;
     973  char *FragmentNumber = NULL;
     974  stringstream file;
     975  char filename[1023];
     976 
     977  FragmentCounter = FCounter;
     978  cout << "Parsing key sets." << endl;
     979  KeySets = (int **) Malloc(sizeof(int *)*FragmentCounter, "KeySetsContainer::ParseKeySets: **KeySets");
     980  for(int i=FragmentCounter;i--;)
     981    KeySets[i] = NULL;
     982  file << name << FRAGMENTPREFIX << KEYSETFILE;
     983  input.open(file.str().c_str(), ios::in);
     984  if (input == NULL) {
     985    cout << endl << "Unable to open " << file.str() << ", is the directory correct?" << endl;
     986    return false;
     987  }
     988 
     989  AtomCounter = (int *) Malloc(sizeof(int)*FragmentCounter, "KeySetsContainer::ParseKeySets: *RowCounter");
     990  for(int i=0;(i<FragmentCounter) && (!input.eof());i++) {
     991    stringstream line;
     992    AtomCounter[i] = ACounter[i];
     993    // parse the values
     994    KeySets[i] = (int *) Malloc(sizeof(int)*AtomCounter[i], "KeySetsContainer::ParseKeySets: *KeySets[]");
     995    for(int j=AtomCounter[i];j--;)
     996      KeySets[i][j] = -1;
     997    FragmentNumber = FixedDigitNumber(FragmentCounter, i);
     998    //cout << FRAGMENTPREFIX << FragmentNumber << "[" << AtomCounter[i] << "]:";
     999    Free((void **)&FragmentNumber, "KeySetsContainer::ParseKeySets: *FragmentNumber");
     1000    input.getline(filename, 1023);
     1001    line.str(filename);
     1002    for(int j=0;(j<AtomCounter[i]) && (!line.eof());j++) {
     1003      line >> KeySets[i][j];
     1004      //cout << " " << KeySets[i][j];
     1005    }
     1006    //cout << endl;
     1007  }
     1008  input.close();
     1009  return true;
    7461010};
    7471011
    … …  
    7511015bool KeySetsContainer::ParseManyBodyTerms()
    7521016{
    753         int Counter;
    754 
    755         cout << "Creating Fragment terms." << endl;
    756         // scan through all to determine maximum order
    757         Order=0;
    758         for(int i=FragmentCounter;i--;) {
    759                 Counter=0;
    760                 for(int j=AtomCounter[i];j--;)
    761                         if (KeySets[i][j] != -1)
    762                                 Counter++;
    763                 if (Counter > Order)
    764                         Order = Counter;
    765         }
    766         cout << "Found Order is " << Order << "." << endl;
    767 
    768         // scan through all to determine fragments per order
    769         FragmentsPerOrder = (int *) Malloc(sizeof(int)*Order, "KeySetsContainer::ParseManyBodyTerms: *FragmentsPerOrder");
    770         for(int i=Order;i--;)
    771                 FragmentsPerOrder[i] = 0;
    772         for(int i=FragmentCounter;i--;) {
    773                 Counter=0;
    774                 for(int j=AtomCounter[i];j--;)
    775                         if (KeySets[i][j] != -1)
    776                                 Counter++;
    777                 FragmentsPerOrder[Counter-1]++;
    778         }
    779         for(int i=0;i<Order;i++)
    780                 cout << "Found No. of Fragments of Order " << i+1 << " is " << FragmentsPerOrder[i] << "." << endl;
    781 
    782         // scan through all to gather indices to each order set
    783         OrderSet = (int **) Malloc(sizeof(int *)*Order, "KeySetsContainer::ParseManyBodyTerms: **OrderSet");
    784         for(int i=Order;i--;)
    785                 OrderSet[i] = (int *) Malloc(sizeof(int)*FragmentsPerOrder[i], "KeySetsContainer::ParseManyBodyTermsKeySetsContainer::ParseManyBodyTerms: *OrderSet[]");
    786         for(int i=Order;i--;)
    787                 FragmentsPerOrder[i] = 0;
    788         for(int i=FragmentCounter;i--;) {
    789                 Counter=0;
    790                 for(int j=AtomCounter[i];j--;)
    791                         if (KeySets[i][j] != -1)
    792                                 Counter++;
    793                 OrderSet[Counter-1][FragmentsPerOrder[Counter-1]] = i;
    794                 FragmentsPerOrder[Counter-1]++;
    795         }
    796         cout << "Printing OrderSet." << endl;
    797         for(int i=0;i<Order;i++) {
    798                 for (int j=0;j<FragmentsPerOrder[i];j++) {
    799                         cout << " " << OrderSet[i][j];
    800                 }
    801                 cout << endl;
    802         }
    803         cout << endl;
    804 
    805 
    806         return true;
     1017  int Counter;
     1018 
     1019  cout << "Creating Fragment terms." << endl;
     1020  // scan through all to determine maximum order
     1021  Order=0;
     1022  for(int i=FragmentCounter;i--;) {
     1023    Counter=0;
     1024    for(int j=AtomCounter[i];j--;)
     1025      if (KeySets[i][j] != -1)
     1026        Counter++;
     1027    if (Counter > Order)
     1028      Order = Counter;
     1029  }
     1030  cout << "Found Order is " << Order << "." << endl;
     1031 
     1032  // scan through all to determine fragments per order
     1033  FragmentsPerOrder = (int *) Malloc(sizeof(int)*Order, "KeySetsContainer::ParseManyBodyTerms: *FragmentsPerOrder");
     1034  for(int i=Order;i--;)
     1035    FragmentsPerOrder[i] = 0;
     1036  for(int i=FragmentCounter;i--;) {
     1037    Counter=0;
     1038    for(int j=AtomCounter[i];j--;)
     1039      if (KeySets[i][j] != -1)
     1040        Counter++;
     1041    FragmentsPerOrder[Counter-1]++;
     1042  }
     1043  for(int i=0;i<Order;i++)
     1044    cout << "Found No. of Fragments of Order " << i+1 << " is " << FragmentsPerOrder[i] << "." << endl;
     1045   
     1046  // scan through all to gather indices to each order set
     1047  OrderSet = (int **) Malloc(sizeof(int *)*Order, "KeySetsContainer::ParseManyBodyTerms: **OrderSet");
     1048  for(int i=Order;i--;)
     1049    OrderSet[i] = (int *) Malloc(sizeof(int)*FragmentsPerOrder[i], "KeySetsContainer::ParseManyBodyTermsKeySetsContainer::ParseManyBodyTerms: *OrderSet[]");
     1050  for(int i=Order;i--;)
     1051    FragmentsPerOrder[i] = 0;
     1052  for(int i=FragmentCounter;i--;) {
     1053    Counter=0;
     1054    for(int j=AtomCounter[i];j--;)
     1055      if (KeySets[i][j] != -1)
     1056        Counter++;
     1057    OrderSet[Counter-1][FragmentsPerOrder[Counter-1]] = i;
     1058    FragmentsPerOrder[Counter-1]++;
     1059  }
     1060  cout << "Printing OrderSet." << endl;
     1061  for(int i=0;i<Order;i++) {
     1062    for (int j=0;j<FragmentsPerOrder[i];j++) {
     1063      cout << " " << OrderSet[i][j];
     1064    }
     1065    cout << endl;
     1066  }
     1067  cout << endl;
     1068 
     1069   
     1070  return true;
    8071071};
    8081072
    … …  
    8141078bool KeySetsContainer::Contains(const int GreaterSet, const int SmallerSet)
    8151079{
    816         bool result = true;
    817         bool intermediate;
    818         if ((GreaterSet < 0) || (SmallerSet < 0) || (GreaterSet > FragmentCounter) || (SmallerSet > FragmentCounter)) // index out of bounds
    819                 return false;
    820         for(int i=AtomCounter[SmallerSet];i--;) {
    821                 intermediate = false;
    822                 for (int j=AtomCounter[GreaterSet];j--;)
    823                         intermediate = (intermediate || ((KeySets[SmallerSet][i] == KeySets[GreaterSet][j]) || (KeySets[SmallerSet][i] == -1)));
    824                 result = result && intermediate;
    825         }
    826 
    827         return result;
     1080  bool result = true;
     1081  bool intermediate;
     1082  if ((GreaterSet < 0) || (SmallerSet < 0) || (GreaterSet > FragmentCounter) || (SmallerSet > FragmentCounter)) // index out of bounds
     1083    return false;
     1084  for(int i=AtomCounter[SmallerSet];i--;) {
     1085    intermediate = false;
     1086    for (int j=AtomCounter[GreaterSet];j--;)
     1087      intermediate = (intermediate || ((KeySets[SmallerSet][i] == KeySets[GreaterSet][j]) || (KeySets[SmallerSet][i] == -1)));
     1088    result = result && intermediate;
     1089  }
     1090
     1091  return result;
    8281092};
    8291093
Note: See TracChangeset for help on using the changeset viewer.