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

/** \file datacreator.cpp
 *
 * Declarations of assisting functions in creating data and plot files.
 *
 */

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

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

#include "CodePatterns/MemDebug.hpp"

#include <fstream>
#include <sstream>
#include <iomanip>

#include "CodePatterns/Log.hpp"
#include "CodePatterns/Verbose.hpp"
#include "LinearAlgebra/defs.hpp"
#include "Helpers/defs.hpp"
#include "datacreator.hpp"
#include "Fragmentation/EnergyMatrix.hpp"
#include "Fragmentation/ForceMatrix.hpp"
#include "Fragmentation/HessianMatrix.hpp"
#include "Fragmentation/KeySetsContainer.hpp"

using namespace std;

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

/** Opens a file with \a *filename in \a *dir.
 * \param output file handle on return
 * \param *dir directory
 * \param *filename name of file
 * \return true if file has been opened
 */
bool OpenOutputFile(ofstream &output, const char *dir, const char *filename)
{
  stringstream name;
  name << dir << "/" << filename;
  output.open(name.str().c_str(), ios::out);
  if (output == NULL) {
    LOG(0, "Unable to open " << name.str() << " for writing, is directory correct?");
    return false;
  }
  return true;
};

/** Opens a file for appending with \a *filename in \a *dir.
 * \param output file handle on return
 * \param *dir directory
 * \param *filename name of file
 * \return true if file has been opened
 */
bool AppendOutputFile(ofstream &output, const char *dir, const char *filename)
{
  stringstream name;
  name << dir << "/" << filename;
  output.open(name.str().c_str(), ios::app);
  if (output == NULL) {
    LOG(0, "Unable to open " << name.str() << " for writing, is directory correct?");
    return false;
  }
  return true;
};

/** Plots an energy vs. order.
 * \param &Fragments EnergyMatrix class containing matrix values
 * \param KeySets KeySetContainer class holding bond KeySetContainer::Order
 * \param *prefix prefix in filename (without ending)
 * \param *msg message to be place in first line as a comment
 * \return true if file was written successfully
 */
bool CreateDataEnergyOrder(class EnergyMatrix &Fragments, class KeySetsContainer &KeySets, const char *dir, const char *prefix, const char *msg, const char *datum)
{
  stringstream filename;
  ofstream output;

  filename << prefix << ".dat";
  if (!OpenOutputFile(output, dir, filename.str().c_str())) return false;
  LOG(0, msg);
  output << "# " << msg << ", created on " << datum;
  output << "#Order\tFrag.No.\t" << Fragments.Header[Fragments.MatrixCounter] << endl;
  for (int BondOrder=0;BondOrder<KeySets.Order;BondOrder++) {
    for(int i=KeySets.FragmentsPerOrder[BondOrder];i--;) {
      for(int j=Fragments.RowCounter[ KeySets.OrderSet[BondOrder][i] ];j--;)
        for(int k=Fragments.ColumnCounter[ KeySets.OrderSet[BondOrder][i] ];k--;)
          Fragments.Matrix[Fragments.MatrixCounter][j][k] += Fragments.Matrix[ KeySets.OrderSet[BondOrder][i] ][j][k];
    }
    output << BondOrder+1 << "\t" << KeySets.FragmentsPerOrder[BondOrder];
    for (int l=0;l<Fragments.ColumnCounter[Fragments.MatrixCounter];l++)
      output << scientific << "\t" << Fragments.Matrix[Fragments.MatrixCounter][ Fragments.RowCounter[Fragments.MatrixCounter]-1 ][l];
    output << endl;
  }
  output.close();
  return true;
};

/** Plots an energy error vs. order.
 * \param &Energy EnergyMatrix class containing reference values (in MatrixCounter matrix)
 * \param &Fragments EnergyMatrix class containing matrix values
 * \param KeySets KeySetContainer class holding bond KeySetContainer::Order
 * \param *prefix prefix in filename (without ending)
 * \param *msg message to be place in first line as a comment
 * \return true if file was written successfully
 */
bool CreateDataDeltaEnergyOrder(class EnergyMatrix &Energy, class EnergyMatrix &Fragments, class KeySetsContainer &KeySets, const char *dir, const char *prefix, const char *msg, const char *datum)
{
  stringstream filename;
  ofstream output;

  filename << prefix << ".dat";
  if (!OpenOutputFile(output, dir, filename.str().c_str())) return false;
  LOG(0, msg);
  output << "# " << msg << ", created on " << datum;
  output << "#Order\tFrag.No.\t" << Fragments.Header[Fragments.MatrixCounter] << endl;
  Fragments.SetLastMatrix(Energy.Matrix[Energy.MatrixCounter],0);
  for (int BondOrder=0;BondOrder<KeySets.Order;BondOrder++) {
    for(int i=KeySets.FragmentsPerOrder[BondOrder];i--;) {
      for(int j=Fragments.RowCounter[ KeySets.OrderSet[BondOrder][i] ];j--;)
        for(int k=Fragments.ColumnCounter[ KeySets.OrderSet[BondOrder][i] ];k--;)
          Fragments.Matrix[Fragments.MatrixCounter][j][k] -= Fragments.Matrix[ KeySets.OrderSet[BondOrder][i] ][j][k];
    }
    output << BondOrder+1 << "\t" << KeySets.FragmentsPerOrder[BondOrder];
    for (int l=0;l<Fragments.ColumnCounter[Energy.MatrixCounter];l++)
      if (fabs(Energy.Matrix[Energy.MatrixCounter][ Energy.RowCounter[Energy.MatrixCounter]-1 ][l]) < MYEPSILON)
        output << scientific << "\t" << Fragments.Matrix[Fragments.MatrixCounter][ Fragments.RowCounter[Fragments.MatrixCounter]-1 ][l];
      else
        output << scientific << "\t" << (Fragments.Matrix[Fragments.MatrixCounter][ Fragments.RowCounter[Fragments.MatrixCounter]-1 ][l] / Energy.Matrix[Energy.MatrixCounter][ Energy.RowCounter[Energy.MatrixCounter]-1 ][l]);
    output << endl;
  }
  output.close();
  return true;
};

/** Plot forces vs. order.
 * \param &Fragments ForceMatrix class containing matrix values
 * \param KeySets KeySetContainer class holding bond KeySetContainer::Order
 * \param *prefix prefix in filename (without ending)
 * \param *msg message to be place in first line as a comment
 * \param *datum current date and time
 * \return true if file was written successfully
 */
bool CreateDataForcesOrder(class ForceMatrix &Fragments, class KeySetsContainer &KeySets, const char *dir, const char *prefix,const  char *msg, const char *datum, void (*CreateForce)(class MatrixContainer &, int))
{
  stringstream filename;
  ofstream output;

  filename << prefix << ".dat";
  if (!OpenOutputFile(output, dir, filename.str().c_str())) return false;
  LOG(0, msg);
  output << "# " << msg << ", created on " << datum;
  output << "# Order\tFrag.No.\t" << Fragments.Header[Fragments.MatrixCounter] << endl;
  Fragments.SetLastMatrix(0.,0);
  for (int BondOrder=0;BondOrder<KeySets.Order;BondOrder++) {
    Fragments.SumSubForces(Fragments, KeySets, BondOrder, 1.);
    output << BondOrder+1 << "\t" << KeySets.FragmentsPerOrder[BondOrder];
    CreateForce(Fragments, Fragments.MatrixCounter);
    for (int l=0;l<Fragments.ColumnCounter[Fragments.MatrixCounter];l++)
       output << scientific << "\t" << Fragments.Matrix[Fragments.MatrixCounter][ Fragments.RowCounter[Fragments.MatrixCounter] ][l];
    output << endl;
  }
  output.close();
  return true;
};

/** Plot forces error vs. order.
 * \param &Force ForceMatrix containing reference values (in MatrixCounter matrix)
 * \param &Fragments ForceMatrix class containing matrix values
 * \param KeySets KeySetContainer class holding bond KeySetContainer::Order
 * \param *prefix prefix in filename (without ending)
 * \param *msg message to be place in first line as a comment
 * \param *datum current date and time
 * \return true if file was written successfully
 */
bool CreateDataDeltaForcesOrder(class ForceMatrix &Force, class ForceMatrix &Fragments, class KeySetsContainer &KeySets, const char *dir, const char *prefix, const char *msg, const char *datum, void (*CreateForce)(class MatrixContainer &, int))
{
  stringstream filename;
  ofstream output;

  filename << prefix << ".dat";
  if (!OpenOutputFile(output, dir, filename.str().c_str())) return false;
  LOG(0, msg);
  output << "# " << msg << ", created on " << datum;
  output << "# Order\tFrag.No.\t" << Fragments.Header[Fragments.MatrixCounter] << endl;
  Fragments.SetLastMatrix(Force.Matrix[Force.MatrixCounter],0);
  for (int BondOrder=0;BondOrder<KeySets.Order;BondOrder++) {
    Fragments.SumSubForces(Fragments, KeySets, BondOrder, -1.);
    output << BondOrder+1 << "\t" << KeySets.FragmentsPerOrder[BondOrder];
    CreateForce(Fragments, Fragments.MatrixCounter);
    for (int l=0;l<Fragments.ColumnCounter[Fragments.MatrixCounter];l++)
       output << scientific << "\t" << Fragments.Matrix[Fragments.MatrixCounter][ Fragments.RowCounter[Fragments.MatrixCounter] ][l];
    output << endl;
  }
  output.close();
  return true;
};

/** Plot forces error vs. vs atom vs. order.
 * \param &Force ForceMatrix containing reference values (in MatrixCounter matrix)
 * \param &Fragments ForceMatrix class containing matrix values
 * \param KeySets KeySetContainer class holding bond KeySetContainer::Order
 * \param *prefix prefix in filename (without ending)
 * \param *msg message to be place in first line as a comment
 * \param *datum current date and time
 * \return true if file was written successfully
 */
bool CreateDataDeltaForcesOrderPerAtom(class ForceMatrix &Force, class ForceMatrix &Fragments, class KeySetsContainer &KeySets, const char *dir, const char *prefix, const char *msg, const char *datum)
{
  stringstream filename;
  ofstream output;
  double norm = 0.;

  filename << prefix << ".dat";
  if (!OpenOutputFile(output, dir, filename.str().c_str())) return false;
  LOG(0, msg);
  output << "# " << msg << ", created on " << datum;
  output << "# AtomNo\t" << Fragments.Header[Fragments.MatrixCounter] << endl;
  Fragments.SetLastMatrix(Force.Matrix[Force.MatrixCounter], 0);
  for (int BondOrder=0;BondOrder<KeySets.Order;BondOrder++) {
    //LOG(0, "Current order is " << BondOrder << ".");
    Fragments.SumSubForces(Fragments, KeySets, BondOrder, -1.);
    // errors per atom
    output << endl << "#Order\t" << BondOrder+1 << endl;
    for(int j=0;j<Fragments.RowCounter[ Fragments.MatrixCounter ];j++) {
      output << Fragments.Indices[Fragments.MatrixCounter][j] << "\t";
      for (int l=0;l<Fragments.ColumnCounter[ Fragments.MatrixCounter ];l++) {
        if (((l+1) % 3) == 0) {
          norm = 0.;
          for (int m=0;m<NDIM;m++)
            norm += Force.Matrix[Force.MatrixCounter][ j ][l+m]*Force.Matrix[Force.MatrixCounter][ j ][l+m];
          norm = sqrt(norm);
        }                                                                                                    
//        if (norm < MYEPSILON)
          output << scientific << Fragments.Matrix[Fragments.MatrixCounter][ j ][l] << "\t";
//        else
//          output << scientific << (Fragments.Matrix[Fragments.MatrixCounter][ j ][l] / norm) << "\t";
      }
      output << endl;
    }
    output << endl;
  }
  output.close();
  return true;
};

/** Plot forces error vs. vs atom vs. order.
 * \param &Fragments ForceMatrix class containing matrix values
 * \param KeySets KeySetContainer class holding bond KeySetContainer::Order
 * \param *prefix prefix in filename (without ending)
 * \param *msg message to be place in first line as a comment
 * \param *datum current date and time
 * \return true if file was written successfully
 */
bool CreateDataForcesOrderPerAtom(class ForceMatrix &Fragments, class KeySetsContainer &KeySets, const char *dir, const char *prefix, const char *msg, const char *datum)
{
  stringstream filename;
  ofstream output;

  filename << prefix << ".dat";
  if (!OpenOutputFile(output, dir, filename.str().c_str())) return false;
  LOG(0, msg);
  output << "# " << msg << ", created on " << datum;
  output << "# AtomNo\t" << Fragments.Header[Fragments.MatrixCounter] << endl;
  for (int BondOrder=0;BondOrder<KeySets.Order;BondOrder++) {
    //LOG(0, "Current order is " << BondOrder << ".");
    Fragments.SumSubForces(Fragments, KeySets, BondOrder, 1.);
    // errors per atom
    output << endl << "#Order\t" << BondOrder+1 << endl;
    for(int j=0;j<Fragments.RowCounter[ Fragments.MatrixCounter ];j++) {
      output << Fragments.Indices[Fragments.MatrixCounter][j] << "\t";
      for (int l=0;l<Fragments.ColumnCounter[ Fragments.MatrixCounter ];l++)
        output << scientific << Fragments.Matrix[Fragments.MatrixCounter][ j ][l] << "\t";
      output << endl;
    }
    output << endl;
  }
  output.close();
  return true;
};


/** Plot hessian error vs. vs atom vs. order.
 * \param &Hessian HessianMatrix containing reference values (in MatrixCounter matrix)
 * \param &Fragments HessianMatrix class containing matrix values
 * \param KeySets KeySetContainer class holding bond KeySetContainer::Order
 * \param *prefix prefix in filename (without ending)
 * \param *msg message to be place in first line as a comment
 * \param *datum current date and time
 * \return true if file was written successfully
 */
bool CreateDataDeltaHessianOrderPerAtom(class HessianMatrix &Hessian, class HessianMatrix &Fragments, class KeySetsContainer &KeySets, const char *dir, const char *prefix, const char *msg, const char *datum)
{
  stringstream filename;
  ofstream output;

  filename << prefix << ".dat";
  if (!OpenOutputFile(output, dir, filename.str().c_str())) return false; 
  LOG(0, msg);
  output << "# " << msg << ", created on " << datum;
  output << "# AtomNo\t" << Fragments.Header[Fragments.MatrixCounter] << endl;
  Fragments.SetLastMatrix(Hessian.Matrix[Hessian.MatrixCounter], 0);
  for (int BondOrder=0;BondOrder<KeySets.Order;BondOrder++) {
    //LOG(0, "Current order is " << BondOrder << ".");
    Fragments.SumSubHessians(Fragments, KeySets, BondOrder, -1.);
    // errors per atom
    output << endl << "#Order\t" << BondOrder+1 << endl;
    for(int j=0;j<Fragments.RowCounter[ Fragments.MatrixCounter ];j++) { 
      output << Fragments.Indices[Fragments.MatrixCounter][j] << "\t";
      for (int l=0;l<Fragments.ColumnCounter[ Fragments.MatrixCounter ];l++) {
        output << scientific << Fragments.Matrix[Fragments.MatrixCounter][ j ][l] << "\t";
      }
      output << endl;
    }
    output << endl;
  }
  output.close();
  return true;
};

/** Plot hessian error vs. vs atom vs. order in the frobenius norm.
 * \param &Hessian HessianMatrix containing reference values (in MatrixCounter matrix)
 * \param &Fragments HessianMatrix class containing matrix values
 * \param KeySets KeySetContainer class holding bond KeySetContainer::Order
 * \param *prefix prefix in filename (without ending)
 * \param *msg message to be place in first line as a comment
 * \param *datum current date and time
 * \return true if file was written successfully
 */
bool CreateDataDeltaFrobeniusOrderPerAtom(class HessianMatrix &Hessian, class HessianMatrix &Fragments, class KeySetsContainer &KeySets, const char *dir, const char *prefix, const char *msg, const char *datum)
{
  stringstream filename;
  ofstream output;
  double norm = 0;
  double tmp;

  filename << prefix << ".dat";
  if (!OpenOutputFile(output, dir, filename.str().c_str())) return false; 
  LOG(0, msg);
  output << "# " << msg << ", created on " << datum;
  output << "# AtomNo\t";
  Fragments.SetLastMatrix(Hessian.Matrix[Hessian.MatrixCounter], 0);
  for (int BondOrder=0;BondOrder<KeySets.Order;BondOrder++) {
    output << "Order" << BondOrder+1 << "\t";
  }
  output << endl;
  output << Fragments.RowCounter[ Fragments.MatrixCounter ] << "\t";
  for (int BondOrder=0;BondOrder<KeySets.Order;BondOrder++) {
    //LOG(0, "Current order is " << BondOrder << ".");
    Fragments.SumSubHessians(Fragments, KeySets, BondOrder, -1.);
    // frobenius norm of errors per atom
    norm = 0.;
    for(int j=0;j<Fragments.RowCounter[ Fragments.MatrixCounter ];j++) { 
      for (int l=0;l<Fragments.ColumnCounter[ Fragments.MatrixCounter ];l++) {
        tmp = Fragments.Matrix[Fragments.MatrixCounter][ j ][l];
        norm += tmp*tmp;
      }
    }
    output << scientific << sqrt(norm)/(Fragments.RowCounter[ Fragments.MatrixCounter ]*Fragments.ColumnCounter[ Fragments.MatrixCounter] ) << "\t";
  }
  output << endl;
  output.close();
  return true;
};

/** Plot hessian error vs. vs atom vs. order.
 * \param &Fragments HessianMatrix class containing matrix values
 * \param KeySets KeySetContainer class holding bond KeySetContainer::Order
 * \param *prefix prefix in filename (without ending)
 * \param *msg message to be place in first line as a comment
 * \param *datum current date and time
 * \return true if file was written successfully
 */
bool CreateDataHessianOrderPerAtom(class HessianMatrix &Fragments, class KeySetsContainer &KeySets, const char *dir, const char *prefix, const char *msg, const char *datum)
{
  stringstream filename;
  ofstream output;

  filename << prefix << ".dat";
  if (!OpenOutputFile(output, dir, filename.str().c_str())) return false; 
  LOG(0, msg);
  output << "# " << msg << ", created on " << datum;
  output << "# AtomNo\t" << Fragments.Header[ Fragments.MatrixCounter ] << endl;
  Fragments.SetLastMatrix(0., 0);
  for (int BondOrder=0;BondOrder<KeySets.Order;BondOrder++) {
    //LOG(0, "Current order is " << BondOrder << ".");
    Fragments.SumSubHessians(Fragments, KeySets, BondOrder, 1.);
    // errors per atom
    output << endl << "#Order\t" << BondOrder+1 << endl;
    for(int j=0;j<Fragments.RowCounter[ Fragments.MatrixCounter ];j++) { 
      output << Fragments.Indices[Fragments.MatrixCounter][j] << "\t";
      for (int l=0;l<Fragments.ColumnCounter[ Fragments.MatrixCounter ];l++)
        output << scientific << Fragments.Matrix[Fragments.MatrixCounter][ j ][l] << "\t";
      output << endl;
    }
    output << endl;
  }
  output.close();
  return true;
};

/** Plot matrix vs. fragment.
 */
bool CreateDataFragment(class MatrixContainer &Fragment, class KeySetsContainer &KeySets, const char *dir, const char *prefix, const char *msg, const char *datum, void (*CreateFragment)(class MatrixContainer &, int))
{
  stringstream filename;
  ofstream output;

  filename << prefix << ".dat";
  if (!OpenOutputFile(output, dir, filename.str().c_str())) return false;
  LOG(0, msg);
  output << "# " << msg << ", created on " << datum << endl;
  output << "#Order\tFrag.No.\t" << Fragment.Header[ Fragment.MatrixCounter ] << endl;
  for (int BondOrder=0;BondOrder<KeySets.Order;BondOrder++) {
    for(int i=0;i<KeySets.FragmentsPerOrder[BondOrder];i++) {
      output << BondOrder+1 << "\t" << KeySets.OrderSet[BondOrder][i]+1;
      CreateFragment(Fragment, KeySets.OrderSet[BondOrder][i]);
      for (int l=0;l<Fragment.ColumnCounter[ KeySets.OrderSet[BondOrder][i] ];l++)
        output << scientific << "\t" << Fragment.Matrix[ KeySets.OrderSet[BondOrder][i] ][ Fragment.RowCounter[ KeySets.OrderSet[BondOrder][i] ] ][l];
      output << endl;
    }
  }
  output.close();
  return true;
};

/** Copies fragment energy values into last matrix of \a Matrix with greatest total energy.
 * \param &Matrix MatrixContainer with all fragment energy values
 * \param &KeySets KeySetsContainer with associations of each fragment to a bond order
 * \param BondOrder current bond order
 */
void CreateMaxFragmentOrder(class MatrixContainer &Fragments, class KeySetsContainer &KeySets, int BondOrder)
{
  for(int j=Fragments.RowCounter[ Fragments.MatrixCounter ];j--;) {
    for(int i=KeySets.FragmentsPerOrder[BondOrder];i--;) {
      if (fabs(Fragments.Matrix[ Fragments.MatrixCounter ][j][1]) < fabs(Fragments.Matrix[ KeySets.OrderSet[BondOrder][i] ][j][1])) {
        for (int k=Fragments.ColumnCounter[ Fragments.MatrixCounter ];k--;)
          Fragments.Matrix[ Fragments.MatrixCounter ][j][k] = Fragments.Matrix[ KeySets.OrderSet[BondOrder][i] ][j][k];
      }
    }
  }
};

/** Copies fragment energy values into last matrix of \a Matrix with smallest total energy.
 * \param &Matrix MatrixContainer with all fragment energy values
 * \param &KeySets KeySetsContainer with associations of each fragment to a bond order
 * \param BondOrder current bond order
 */
void CreateMinFragmentOrder(class MatrixContainer &Fragments, class KeySetsContainer &KeySets, int BondOrder)
{
  for(int j=0;j<Fragments.RowCounter[ Fragments.MatrixCounter ];j++) {
    int i=0;
    do {  // first get a minimum value unequal to 0
      for (int k=Fragments.ColumnCounter[ Fragments.MatrixCounter ];k--;)
        Fragments.Matrix[ Fragments.MatrixCounter ][j][k] = Fragments.Matrix[ KeySets.OrderSet[BondOrder][i] ][j][k];
      i++;
    } while ((fabs(Fragments.Matrix[ Fragments.MatrixCounter ][j][1]) < MYEPSILON) && (i<KeySets.FragmentsPerOrder[BondOrder]));
    for(;i<KeySets.FragmentsPerOrder[BondOrder];i++) { // then find lowest
      if (fabs(Fragments.Matrix[ Fragments.MatrixCounter ][j][1]) > fabs(Fragments.Matrix[ KeySets.OrderSet[BondOrder][i] ][j][1])) {
        for (int k=Fragments.ColumnCounter[ Fragments.MatrixCounter ];k--;)
          Fragments.Matrix[ Fragments.MatrixCounter ][j][k] = Fragments.Matrix[ KeySets.OrderSet[BondOrder][i] ][j][k];
      }
    }
  }
};

/** Plot matrix vs. fragment.
 */
bool CreateDataFragmentOrder(class MatrixContainer &Fragment, class KeySetsContainer &KeySets, const char *dir, const char *prefix, const char *msg, const char *datum, void (*CreateFragmentOrder)(class MatrixContainer &, class KeySetsContainer &, int))
{
  stringstream filename;
  ofstream output;

  filename << prefix << ".dat";
  if (!OpenOutputFile(output, dir, filename.str().c_str())) return false;
  LOG(0, msg);
  output << "# " << msg << ", created on " << datum;
  output << "#Order\tFrag.No.\t" << Fragment.Header[ Fragment.MatrixCounter ] << endl;
  // max
  for (int BondOrder=0;BondOrder<KeySets.Order;BondOrder++) {
    Fragment.SetLastMatrix(0.,0);
    CreateFragmentOrder(Fragment, KeySets, BondOrder);
    output << BondOrder+1 << "\t" << KeySets.FragmentsPerOrder[BondOrder];
    for (int l=0;l<Fragment.ColumnCounter[ Fragment.MatrixCounter ];l++)
      output << scientific << "\t" << Fragment.Matrix[ Fragment.MatrixCounter ][ Fragment.RowCounter[ Fragment.MatrixCounter ]-1 ][l];
    output << endl;
  }
  output.close();
  return true;
};

/** Takes last but one row and copies into final row.
 * \param Energy MatrixContainer with matrix values
 * \param MatrixNumber the index for the ForceMatrix::matrix array
 */
void CreateEnergy(class MatrixContainer &Energy, int MatrixNumber)
{
  for(int k=0;k<Energy.ColumnCounter[MatrixNumber];k++)
    Energy.Matrix[MatrixNumber][ Energy.RowCounter[MatrixNumber] ] [k] =  Energy.Matrix[MatrixNumber][ Energy.RowCounter[MatrixNumber]-1 ] [k];
};

/** Scans forces for the minimum in magnitude.
 * Results are stored in the matrix ForceMatrix::MatrixCounter of \a Force.
 * \param Force ForceMatrix class containing matrix values
 * \param MatrixNumber the index for the ForceMatrix::matrix array
 */
void CreateMinimumForce(class MatrixContainer &Force, int MatrixNumber)
{
  for (int l=Force.ColumnCounter[MatrixNumber];l--;)
    Force.Matrix[MatrixNumber][ Force.RowCounter[MatrixNumber] ][l] = 0.;
  for (int l=5;l<Force.ColumnCounter[MatrixNumber];l+=3) {
    double stored = 0;
    int k=0;
    do {
      for (int m=NDIM;m--;) {
        stored += Force.Matrix[MatrixNumber][ k ][l+m]
              * Force.Matrix[MatrixNumber][ k ][l+m];
        Force.Matrix[MatrixNumber][ Force.RowCounter[MatrixNumber] ][l+m]  = Force.Matrix[MatrixNumber][ k ][l+m];
      }
      k++;
    } while ((fabs(stored) < MYEPSILON) && (k<Force.RowCounter[MatrixNumber]));
    for (;k<Force.RowCounter[MatrixNumber];k++) {
      double tmp = 0;
      for (int m=NDIM;m--;)
        tmp += Force.Matrix[MatrixNumber][ k ][l+m]
              * Force.Matrix[MatrixNumber][ k ][l+m];
      if ((fabs(tmp) > MYEPSILON) && (tmp < stored)) {  // current force is greater than stored
        for (int m=NDIM;m--;)
          Force.Matrix[MatrixNumber][ Force.RowCounter[MatrixNumber] ][l+m]  = Force.Matrix[MatrixNumber][ k ][l+m];
        stored = tmp;
      }
    }
  }
};

/** Scans forces for the mean in magnitude.
 * Results are stored in the matrix ForceMatrix::MatrixCounter of \a Force.
 * \param Force ForceMatrix class containing matrix values
  * \param MatrixNumber the index for the ForceMatrix::matrix array
 */
void CreateMeanForce(class MatrixContainer &Force, int MatrixNumber)
{
  int divisor = 0;
  for (int l=Force.ColumnCounter[MatrixNumber];l--;)
    Force.Matrix[MatrixNumber][ Force.RowCounter[MatrixNumber] ][l] = 0.;
  for (int l=5;l<Force.ColumnCounter[MatrixNumber];l+=3) {
    double tmp = 0;
    for (int k=Force.RowCounter[MatrixNumber];k--;) {
      double norm = 0.;
      for (int m=NDIM;m--;)
        norm += Force.Matrix[MatrixNumber][ k ][l+m]
              * Force.Matrix[MatrixNumber][ k ][l+m];
      tmp += sqrt(norm);
      if (fabs(norm) > MYEPSILON) divisor++;
    }
    tmp /= (double)divisor;
    Force.Matrix[MatrixNumber][ Force.RowCounter[MatrixNumber] ][l] = tmp;
  }
};

/** Scans forces for the maximum in magnitude.
 * Results are stored in the matrix ForceMatrix::MatrixCounter of \a Force.
 * \param Force ForceMatrix class containing matrix values
 * \param MatrixNumber the index for the ForceMatrix::matrix array
 */
void CreateMaximumForce(class MatrixContainer &Force, int MatrixNumber)
{
  for (int l=5;l<Force.ColumnCounter[MatrixNumber];l+=3) {
    double stored = 0;
    for (int k=Force.RowCounter[MatrixNumber];k--;) {
      double tmp = 0;
      for (int m=NDIM;m--;)
        tmp += Force.Matrix[MatrixNumber][ k ][l+m]
              * Force.Matrix[MatrixNumber][ k ][l+m];
      if (tmp > stored) {  // current force is greater than stored
        for (int m=NDIM;m--;)
          Force.Matrix[MatrixNumber][ Force.RowCounter[MatrixNumber] ][l+m]  = Force.Matrix[MatrixNumber][ k ][l+m];
        stored = tmp;
      }
    }
  }
};

/** Leaves the Force.Matrix as it is.
 * \param Force ForceMatrix class containing matrix values
 * \param MatrixNumber the index for the ForceMatrix::matrix array
*/
void CreateSameForce(class MatrixContainer &Force, int MatrixNumber)
{
  // does nothing
};

/** Adds vectorwise all forces.
 * Results are stored in the matrix ForceMatrix::MatrixCounter of \a Force.
 * \param Force ForceMatrix class containing matrix values
 * \param MatrixNumber the index for the ForceMatrix::matrix array
 */
void CreateVectorSumForce(class MatrixContainer &Force, int MatrixNumber)
{
  for (int l=Force.ColumnCounter[MatrixNumber];l--;)
    Force.Matrix[MatrixNumber][ Force.RowCounter[MatrixNumber] ][l] = 0.;
  for (int l=0;l<Force.ColumnCounter[MatrixNumber];l++) {
    for (int k=Force.RowCounter[MatrixNumber];k--;)
      Force.Matrix[MatrixNumber][ Force.RowCounter[MatrixNumber] ][l] += Force.Matrix[MatrixNumber][k][l];
  }
};

/** Writes the standard pyxplot header info.
 * \param keycolumns number of columns of the key
 * \param *key position of key
 * \param *logscale axis for logscale
 * \param *extraline extra set lines if desired
 * \param mxtics small tics at ...
 * \param xtics large tics at ...
 * \param *xlabel label for x axis
 * \param *ylabel label for y axis
 */
void CreatePlotHeader(ofstream &output, const char *prefix, const int keycolumns, const char *key, const char *logscale, const char *extraline, const int mxtics, const int xtics, const char *xlabel, const char *ylabel)
{
  //output << "#!/home/heber/build/pyxplot/pyxplot" << endl << endl;
  output << "reset" << endl;
  output << "set keycolumns "<< keycolumns << endl;
  output << "set key " << key << endl;
  output << "set mxtics "<< mxtics << endl;
  output << "set xtics "<< xtics << endl;
  if (logscale != NULL)
    output << "set logscale " << logscale << endl;
  if (extraline != NULL)
    output << extraline << endl;
  output << "set xlabel '" << xlabel << "'" << endl;
  output << "set ylabel '" << ylabel << "'" << endl;
  output << "set terminal eps color" << endl;
  output << "set output '"<< prefix << ".eps'" << endl;
};

/** Creates the pyxplotfile for energy data.
 * \param Matrix MatrixContainer with matrix values
 * \param KeySets contains bond order
 * \param *dir directory
 * \param *prefix prefix for all filenames (without ending)
 * \param keycolumns number of columns of the key
 * \param *key position of key
 * \param logscale axis for logscale
 * \param mxtics small tics at ...
 * \param xtics large tics at ...
 * \param xlabel label for x axis
 * \param xlabel label for x axis
 * \param *xrange xrange
 * \param *yrange yrange
 * \param *xargument number of column to plot against (x axis)
 * \param uses using string for plot command
 * \param (*CreatePlotLines) function reference that writes a single plot line
 * \return true if file was written successfully
 */
bool CreatePlotOrder(class MatrixContainer &Matrix, const class KeySetsContainer &KeySets, const char *dir, const char *prefix, const int keycolumns, const char *key, const char *logscale, const char *extraline, const int mxtics, const int xtics, const char *xlabel, const char *ylabel, const char *xrange, const char *yrange, const char *xargument, const char *uses, void (*CreatePlotLines)(ofstream &, class MatrixContainer &, const char *, const char *, const char *))
{
  stringstream filename;
  ofstream output;

  filename << prefix << ".pyx";
  if (!OpenOutputFile(output, dir, filename.str().c_str())) return false;
  CreatePlotHeader(output, prefix, keycolumns, key, logscale, extraline, mxtics, xtics, xlabel, ylabel);
  output << "plot " << xrange << " " << yrange << " \\" << endl;
  CreatePlotLines(output, Matrix, prefix, xargument, uses);
  output.close();
  return true;
};

/** Writes plot lines for absolute energies.
 * \param output file handler
 * \param Energy MatrixContainer with matrix values
 * \param *prefix prefix of data file
 * \param *xargument number of column to plot against (x axis)
 * \param *uses uses command
 */
void AbsEnergyPlotLine(ofstream &output, class MatrixContainer &Energy, const char *prefix, const char *xargument, const char *uses)
{
  stringstream line(Energy.Header[ Energy.MatrixCounter ]);
  string token;

  getline(line, token, '\t');
  for (int i=2; i<= Energy.ColumnCounter[Energy.MatrixCounter];i++) {
    getline(line, token, '\t');
    while (token[0] == ' ') // remove leading white spaces
      token.erase(0,1);
    output << "'" << prefix << ".dat' title '" << token << "' using " << xargument << ":(abs($" << i+2 << ")) " << uses;
    if (i != (Energy.ColumnCounter[Energy.MatrixCounter]))
      output << ", \\";
    output << endl;
  }
};

/** Writes plot lines for energies.
 * \param output file handler
 * \param Energy MatrixContainer with matrix values
 * \param *prefix prefix of data file
 * \param *xargument number of column to plot against (x axis)
 * \param *uses uses command
 */
void EnergyPlotLine(ofstream &output, class MatrixContainer &Energy, const char *prefix, const char *xargument, const char *uses)
{
  stringstream line(Energy.Header[Energy.MatrixCounter]);
  string token;

  getline(line, token, '\t');
  for (int i=1; i<= Energy.ColumnCounter[Energy.MatrixCounter];i++) {
    getline(line, token, '\t');
    while (token[0] == ' ') // remove leading white spaces
      token.erase(0,1);
    output << "'" << prefix << ".dat' title '" << token << "' using " << xargument << ":" << i+2 << " " << uses;
    if (i != (Energy.ColumnCounter[Energy.MatrixCounter]))
      output << ", \\";
    output << endl;
  }
};

/** Writes plot lines for absolute force magnitudes.
 * \param output file handler
 * \param Force MatrixContainer with matrix values
 * \param *prefix prefix of data file
 * \param *xargument number of column to plot against (x axis)
 * \param *uses uses command
 */
void ForceMagnitudePlotLine(ofstream &output, class MatrixContainer &Force, const char *prefix, const char *xargument, const char *uses)
{
  stringstream line(Force.Header[Force.MatrixCounter]);
  string token;

  getline(line, token, '\t');
  getline(line, token, '\t');
  getline(line, token, '\t');
  getline(line, token, '\t');
  getline(line, token, '\t');
  for (int i=7; i< Force.ColumnCounter[Force.MatrixCounter];i+=NDIM) {
    getline(line, token, '\t');
    while (token[0] == ' ') // remove leading white spaces
      token.erase(0,1);
    token.erase(token.length(), 1);  // kill residual index char (the '0')
    output << "'" << prefix << ".dat' title '" << token << "' using " << xargument << ":(sqrt($" << i+1 << "*$" << i+1 << "+$" << i+2 << "*$" << i+2 << "+$" << i+3 << "*$" << i+3 << ")) " << uses;
    if (i != (Force.ColumnCounter[Force.MatrixCounter]-1))
      output << ", \\";
    output << endl;
    getline(line, token, '\t');
    getline(line, token, '\t');
  }
};

/** Writes plot lines for first component of force vector.
 * \param output file handler
 * \param Force MatrixContainer with matrix values
 * \param *prefix prefix of data file
 * \param *xargument number of column to plot against (x axis)
 * \param *uses uses command
 */
void AbsFirstForceValuePlotLine(ofstream &output, class MatrixContainer &Force, const char *prefix, const char *xargument, const char *uses)
{
  stringstream line(Force.Header[Force.MatrixCounter]);
  string token;

  getline(line, token, '\t');
  getline(line, token, '\t');
  getline(line, token, '\t');
  getline(line, token, '\t');
  getline(line, token, '\t');
  for (int i=7; i< Force.ColumnCounter[Force.MatrixCounter];i+=NDIM) {
    getline(line, token, '\t');
    while (token[0] == ' ') // remove leading white spaces
      token.erase(0,1);
    token.erase(token.length(), 1);  // kill residual index char (the '0')
    output << "'" << prefix << ".dat' title '" << token << "' using " << xargument << ":(abs($" << i+1 << ")) " << uses;
    if (i != (Force.ColumnCounter[Force.MatrixCounter]-1))
      output << ", \\";
    output << endl;
    getline(line, token, '\t');
    getline(line, token, '\t');
  }
};

/** Writes plot lines for force vector as boxes with a fillcolor.
 * \param output file handler
 * \param Force MatrixContainer with matrix values
 * \param *prefix prefix of data file
 * \param *xargument number of column to plot against (x axis)
 * \param *uses uses command
 */
void BoxesForcePlotLine(ofstream &output, class MatrixContainer &Force, const char *prefix, const char *xargument, const char *uses)
{
  stringstream line(Force.Header[Force.MatrixCounter]);
  const char *fillcolor[5] = {"black", "red", "blue", "green", "cyan"};
  string token;

  getline(line, token, '\t');
  getline(line, token, '\t');
  getline(line, token, '\t');
  getline(line, token, '\t');
  getline(line, token, '\t');
  for (int i=7; i< Force.ColumnCounter[Force.MatrixCounter];i+=NDIM) {
    getline(line, token, '\t');
    while (token[0] == ' ') // remove leading white spaces
      token.erase(0,1);
    token.erase(token.length(), 1);  // kill residual index char (the '0')
    output << "'" << prefix << ".dat' title '" << token << "' using ($" << xargument << "+" << fixed << setprecision(1) << (double)((i-7)/3)*0.2 << "):(sqrt($" << i+1 << "*$" << i+1 << "+$" << i+2 << "*$" << i+2 << "+$" << i+3 << "*$" << i+3 << ")) " << uses << " " << fillcolor[(i-7)/3];
    if (i != (Force.ColumnCounter[Force.MatrixCounter]-1))
      output << ", \\";
    output << endl;
    getline(line, token, '\t');
    getline(line, token, '\t');
  }
};

/** Writes plot lines for first force vector component as boxes with a fillcolor.
 * \param output file handler
 * \param Force MatrixContainer with matrix values
 * \param *prefix prefix of data file
 * \param *xargument number of column to plot against (x axis)
 * \param *uses uses command
 */
void BoxesFirstForceValuePlotLine(ofstream &output, class MatrixContainer &Force, const char *prefix, const char *xargument, const char *uses)
{
  stringstream line(Force.Header[Force.MatrixCounter]);
  const char *fillcolor[5] = {"black", "red", "blue", "green", "cyan"};
  string token;

  getline(line, token, '\t');
  getline(line, token, '\t');
  getline(line, token, '\t');
  getline(line, token, '\t');
  getline(line, token, '\t');
  for (int i=7; i< Force.ColumnCounter[Force.MatrixCounter];i+=NDIM) {
    getline(line, token, '\t');
    while (token[0] == ' ') // remove leading white spaces
      token.erase(0,1);
    token.erase(token.length(), 1);  // kill residual index char (the '0')
    output << "'" << prefix << ".dat' title '" << token << "' using ($" << xargument << "+" << fixed << setprecision(1) << (double)((i-7)/3)*0.2 << "):" << i+1 << " " << uses << " " << fillcolor[(i-7)/3];
    if (i != (Force.ColumnCounter[Force.MatrixCounter]-1))
      output << ", \\";
    output << endl;
    getline(line, token, '\t');
    getline(line, token, '\t');
  }
};