source: molecuilder/src/datacreator.cpp@ 0f44198

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

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

#include "datacreator.hpp"

//=========================== 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) {
    cout << "Unable to open " << name.str() << " for writing, is directory correct?" << endl;
    return false;
  }
  return true;
}; 

/** Plots an energy vs. order.
 * \param Energy EnergyMatrix class containing matrix values (last matrix is created by summing over Fragments)
 * \param EnergyFragments EnergyMatrix class containing matrix values
 * \param KeySet 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 MatrixContainer &Energy, class MatrixContainer &EnergyFragments, class KeySetsContainer &KeySet, char *dir, char *prefix, char *msg, char *datum) 
{
  stringstream filename;
  ofstream output;

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

/** Plot forces vs. order.
 */
bool CreateDataForcesOrder(class MatrixContainer &Force, class MatrixContainer &ForceFragments, class KeySetsContainer &KeySet, char *dir, char *prefix, char *msg, char *datum, void (*CreateForce)(class MatrixContainer &, int))
{
  stringstream filename;
  ofstream output;

  filename << prefix << ".dat";
  if (!OpenOutputFile(output, dir, filename.str().c_str())) return false; 
  cout << msg << endl;
  output << "# " << msg << ", created on " << datum;
  output << "# Order\tFrag.No.\t" << Force.Header << endl;
  for(int j=0;j<=Force.RowCounter[ Force.MatrixCounter ];j++)
    for(int k=0;k<Force.ColumnCounter;k++)
       Force.Matrix[Force.MatrixCounter][j][k] = 0.;
  for (int BondOrder=0;BondOrder<KeySet.Order;BondOrder++) {
    for(int i=0;i<KeySet.FragmentsPerOrder[BondOrder];i++) {
      for(int l=0;l<Force.RowCounter[ KeySet.OrderSet[BondOrder][i] ];l++) {
        int j = Force.Indices[ KeySet.OrderSet[BondOrder][i] ][l];
        if (j > Force.RowCounter[Force.MatrixCounter]) {
          cerr << "Current force index " << j << " is greater than " << Force.RowCounter[Force.MatrixCounter] << "!" << endl;
          return 1;
        }
        if (j != -1)
          for(int k=Force.ColumnCounter;k--;) {
            Force.Matrix[Force.MatrixCounter][j][k] += ForceFragments.Matrix[ KeySet.OrderSet[BondOrder][i] ][l][k];
          }
      }
    }
    output << BondOrder+1 << "\t" << KeySet.FragmentsPerOrder[BondOrder];
    CreateForce(Force, Force.MatrixCounter);
    for (int l=0;l<Force.ColumnCounter;l++)
       output << scientific << "\t" << Force.Matrix[Force.MatrixCounter][ Force.RowCounter[Force.MatrixCounter] ][l];
    output << endl;
  }
  output.close();
  return true;
};

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

  filename << prefix << ".dat";
  if (!OpenOutputFile(output, dir, filename.str().c_str())) return false; 
  cout << msg << endl;
  output << "# " << msg << ", created on " << datum << endl;
  output << "#Order\tFrag.No.\t" << Fragment.Header << endl;
  for (int BondOrder=0;BondOrder<KeySet.Order;BondOrder++) {
    for(int i=0;i<KeySet.FragmentsPerOrder[BondOrder];i++) {
      output << BondOrder+1 << "\t" << KeySet.OrderSet[BondOrder][i]+1;
      CreateFragment(Fragment, KeySet.OrderSet[BondOrder][i]);
      for (int l=0;l<Fragment.ColumnCounter;l++)
        output << scientific << "\t" << Fragment.Matrix[ KeySet.OrderSet[BondOrder][i] ][ Fragment.RowCounter[ KeySet.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 &KeySet KeySetsContainer with associations of each fragment to a bond order
 * \param BondOrder current bond order
 */ 
void CreateMaxFragmentOrder(class MatrixContainer &Fragments, class KeySetsContainer &KeySet, int BondOrder)
{
  for(int j=Fragments.RowCounter[ Fragments.MatrixCounter ];j--;) {
    for(int i=KeySet.FragmentsPerOrder[BondOrder];i--;) {
      if (fabs(Fragments.Matrix[ Fragments.MatrixCounter ][j][1]) < fabs(Fragments.Matrix[ KeySet.OrderSet[BondOrder][i] ][j][1])) {
        for (int k=Fragments.ColumnCounter;k--;)
          Fragments.Matrix[ Fragments.MatrixCounter ][j][k] = Fragments.Matrix[ KeySet.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 &KeySet KeySetsContainer with associations of each fragment to a bond order
 * \param BondOrder current bond order
 */ 
void CreateMinFragmentOrder(class MatrixContainer &Fragments, class KeySetsContainer &KeySet, 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;k--;)
        Fragments.Matrix[ Fragments.MatrixCounter ][j][k] = Fragments.Matrix[ KeySet.OrderSet[BondOrder][i] ][j][k];
      i++;
    } while ((fabs(Fragments.Matrix[ Fragments.MatrixCounter ][j][1]) < MYEPSILON) && (i<KeySet.FragmentsPerOrder[BondOrder]));
    for(;i<KeySet.FragmentsPerOrder[BondOrder];i++) { // then find lowest
      if (fabs(Fragments.Matrix[ Fragments.MatrixCounter ][j][1]) > fabs(Fragments.Matrix[ KeySet.OrderSet[BondOrder][i] ][j][1])) {
        for (int k=Fragments.ColumnCounter;k--;)
          Fragments.Matrix[ Fragments.MatrixCounter ][j][k] = Fragments.Matrix[ KeySet.OrderSet[BondOrder][i] ][j][k];
      }
    }      
  }
};

/** Plot matrix vs. fragment.
 */
bool CreateDataFragmentOrder(class MatrixContainer &Fragment, class KeySetsContainer &KeySet, char *dir, char *prefix, char *msg, 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; 
  cout << msg << endl;
  output << "# " << msg << ", created on " << datum;
  output << "#Order\tFrag.No.\t" << Fragment.Header << endl;
  // max
  for (int BondOrder=0;BondOrder<KeySet.Order;BondOrder++) {
    Fragment.SetLastMatrix(0.,0);
    CreateFragmentOrder(Fragment, KeySet, BondOrder);
    output << BondOrder+1 << "\t" << KeySet.FragmentsPerOrder[BondOrder];
    for (int l=0;l<Fragment.ColumnCounter;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;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;l--;)
    Force.Matrix[MatrixNumber][ Force.RowCounter[MatrixNumber] ][l] = 0.;
  for (int l=5;l<Force.ColumnCounter;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;l--;)
    Force.Matrix[MatrixNumber][ Force.RowCounter[MatrixNumber] ][l] = 0.;
  for (int l=5;l<Force.ColumnCounter;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;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;
      }
    }
  }
};

/** 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;l--;)
    Force.Matrix[MatrixNumber][ Force.RowCounter[MatrixNumber] ][l] = 0.;
  for (int l=5;l<Force.ColumnCounter;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 KeySet 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 &KeySet, 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)
{
  char item[1024];
  stringstream line(Energy.Header);

  line >> item;
  for (int i=3; i< Energy.ColumnCounter;i++) {
    line >> item;
    output << "'" << prefix << ".dat' title '" << item << "' using " << xargument << ":(abs($" << i+1 << ")) " << uses;
    if (i != (Energy.ColumnCounter-1))
      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)
{
  char item[1024];
  stringstream line(Energy.Header);

  line >> item;
  for (int i=3; i< Energy.ColumnCounter;i++) {
    line >> item;
    output << "'" << prefix << ".dat' title '" << item << "' using " << xargument << ":" << i+1 << " " << uses;
    if (i != (Energy.ColumnCounter-1))
      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)
{
  char item[1024];
  stringstream line(Force.Header);

  line >> item;
  line >> item;
  line >> item;
  line >> item;
  line >> item;
  for (int i=7; i< Force.ColumnCounter;i+=NDIM) {
    line >> item;
    item[strlen(item)-1] = '\0';  // kill residual index char (the '0')
    output << "'" << prefix << ".dat' title '" << item << "' using " << xargument << ":(sqrt($" << i+1 << "*$" << i+1 << "+$" << i+2 << "*$" << i+2 << "+$" << i+3 << "*$" << i+3 << ")) " << uses;
    if (i != (Force.ColumnCounter-1))
      output << ", \\";
    output << endl;
    line >> item;
    line >> item;
  }
};

/** 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)
{
  char item[1024];
  stringstream line(Force.Header);

  line >> item;
  line >> item;
  line >> item;
  line >> item;
  line >> item;
  for (int i=7; i< Force.ColumnCounter;i+=NDIM) {
    line >> item;
    item[strlen(item)-1] = '\0';  // kill residual index char (the '0')
    output << "'" << prefix << ".dat' title '" << item << "' using " << xargument << ":(abs($" << i+1 << ")) " << uses;
    if (i != (Force.ColumnCounter-1))
      output << ", \\";
    output << endl;
    line >> item;
    line >> item;
  }
};

/** 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)
{
  char item[1024];
  stringstream line(Force.Header);
  char *fillcolor[5] = {"black", "red", "blue", "green", "cyan"};

  line >> item;
  line >> item;
  line >> item;
  line >> item;
  line >> item;
  for (int i=7; i< Force.ColumnCounter;i+=NDIM) {
    line >> item;
    item[strlen(item)-1] = '\0';  // kill residual index char (the '0')
    output << "'" << prefix << ".dat' title '" << item << "' 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-1))
      output << ", \\";
    output << endl;
    line >> item;
    line >> item;
  }
};

/** 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)
{
  char item[1024];
  stringstream line(Force.Header);
  char *fillcolor[5] = {"black", "red", "blue", "green", "cyan"};

  line >> item;
  line >> item;
  line >> item;
  line >> item;
  line >> item;
  for (int i=7; i< Force.ColumnCounter;i+=NDIM) {
    line >> item;
    item[strlen(item)-1] = '\0';  // kill residual index char (the '0')
    output << "'" << prefix << ".dat' title '" << item << "' using ($" << xargument << "+" << fixed << setprecision(1) << (double)((i-7)/3)*0.2 << "):" << i+1 << " " << uses << " " << fillcolor[(i-7)/3];
    if (i != (Force.ColumnCounter-1))
      output << ", \\";
    output << endl;
    line >> item;
    line >> item;
  }
};