/*
 * Project: MoleCuilder
 * Description: creates and alters molecular systems
 * Copyright (C)  2010-2012 University of Bonn. All rights reserved.
 * 
 *
 *   This file is part of MoleCuilder.
 *
 *    MoleCuilder is free software: you can redistribute it and/or modify
 *    it under the terms of the GNU General Public License as published by
 *    the Free Software Foundation, either version 2 of the License, or
 *    (at your option) any later version.
 *
 *    MoleCuilder is distributed in the hope that it will be useful,
 *    but WITHOUT ANY WARRANTY; without even the implied warranty of
 *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *    GNU General Public License for more details.
 *
 *    You should have received a copy of the GNU General Public License
 *    along with MoleCuilder.  If not, see <http://www.gnu.org/licenses/>.
 */

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

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

#include "CodePatterns/MemDebug.hpp"

#include <cstring>
#include <fstream>

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

#include "EnergyMatrix.hpp"

/** 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()
{
  LOG(0, "Parsing energy indices.");
  Indices.resize(MatrixCounter + 1);
  for(int i=MatrixCounter+1;i--;) {
    Indices[i].resize(RowCounter[i]);
    for(int j=RowCounter[i];j--;)
      Indices[i][j] = j;
  }
  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 KeySets 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 &KeySets, int Order, double sign)
{
  // sum energy
  if (CorrectionFragments == NULL)
    for(int i=KeySets.FragmentsPerOrder[Order];i--;)
      for(int j=RowCounter[ KeySets.OrderSet[Order][i] ];j--;)
        for(int k=ColumnCounter[ KeySets.OrderSet[Order][i] ];k--;)
          Matrix[MatrixCounter][j][k] += sign*Fragments.Matrix[ KeySets.OrderSet[Order][i] ][j][k];
  else
    for(int i=KeySets.FragmentsPerOrder[Order];i--;)
      for(int j=RowCounter[ KeySets.OrderSet[Order][i] ];j--;)
        for(int k=ColumnCounter[ KeySets.OrderSet[Order][i] ];k--;)
          Matrix[MatrixCounter][j][k] += sign*(Fragments.Matrix[ KeySets.OrderSet[Order][i] ][j][k] + CorrectionFragments->Matrix[ KeySets.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(const char *name, const char *prefix, std::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;
    ColumnCounter[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];
      if (ColumnCounter[j] > ColumnCounter[MatrixCounter])  // take maximum of all for last matrix
        ColumnCounter[MatrixCounter] = ColumnCounter[j];
    }
    // allocate last plus one matrix
    if ((int)Matrix[MatrixCounter].size() <= RowCounter[MatrixCounter] + 2)
      Matrix[MatrixCounter].resize(RowCounter[MatrixCounter] + 1);
    for(int j=0;j<=RowCounter[MatrixCounter];j++)
      if ((int)Matrix[MatrixCounter][j].size() <= ColumnCounter[MatrixCounter]+1)
        Matrix[MatrixCounter][j].resize(ColumnCounter[MatrixCounter]);

    // 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));
    std::ifstream input(filename);
    ParseMatrix(input, skiplines, skipcolumns, MatrixCounter);
    input.close();
  }
  return status;
};