source: src/helpers.cpp@ fc6053

/** \file helpers.cpp
 *
 * Implementation of some auxiliary functions for memory dis-/allocation and so on
 */

#include "Helpers/MemDebug.hpp"

#include "helpers.hpp"
#include "Helpers/fast_functions.hpp"
#include "log.hpp"
#include "memoryusageobserver.hpp"

/********************************************** helpful functions *********************************/


/** Asks for a double value and checks input
 * \param *text question
 */
double ask_value(const char *text)
{
  double test = 0.1439851348959832147598734598273456723948652983045928346598365;
  do {
    DoLog(0) && (Log() << Verbose(0) << text);
    cin >> test;
  } while (test == 0.1439851348959832147598734598273456723948652983045928346598365);
  return test;
};

/** Output of a debug message to stderr.
 * \param *P Problem at hand, points to ParallelSimulationData#me
 * \param output output string
 */
#ifdef HAVE_DEBUG
void debug_in(const char *output, const char *file, const int line) {
  if (output) fprintf(stderr,"DEBUG: in %s at line %i: %s\n", file, line, output);
}
#else
void debug_in(const char *output, const char *file, const int line) {}  // print nothing
#endif

/** modulo operator for doubles.
 * \param *b pointer to double
 * \param lower_bound lower bound
 * \param upper_bound upper bound
 */
void bound(double *b, double lower_bound, double upper_bound)
{
  double step = (upper_bound - lower_bound);
  while (*b >= upper_bound)
    *b -= step;
  while (*b < lower_bound)
    *b += step;
};

/** Counts lines in file.
 * Note we are scanning lines from current position, not from beginning.
 * \param InputFile file to be scanned.
 */
int CountLinesinFile(ifstream &InputFile)
{
  char *buffer = new char[MAXSTRINGSIZE];
  int lines=0;

  int PositionMarker = InputFile.tellg();  // not needed as Inputfile is copied, given by value, not by ref
  // count the number of lines, i.e. the number of fragments
  InputFile.getline(buffer, MAXSTRINGSIZE); // skip comment lines
  InputFile.getline(buffer, MAXSTRINGSIZE);
  while(!InputFile.eof()) {
    InputFile.getline(buffer, MAXSTRINGSIZE);
    lines++;
  }
  InputFile.seekg(PositionMarker, ios::beg);
  delete[](buffer);
  return lines;
};

/** Returns a string with \a i prefixed with 0s to match order of total number of molecules in digits.
 * \param FragmentNumber total number of fragments to determine necessary number of digits
 * \param digits number to create with 0 prefixed
 * \return allocated(!) char array with number in digits, ten base.
 */
char *FixedDigitNumber(const int FragmentNumber, const int digits)
{
  char *returnstring;
  int number = FragmentNumber;
  int order = 0;
  while (number != 0) { // determine number of digits needed
    number = (int)floor(((double)number / 10.));
    order++;
    //Log() << Verbose(0) << "Number is " << number << ", order is " << order << "." << endl;
  }
  // allocate string
  returnstring = new char[order + 2];
  // terminate  and fill string array from end backward
  returnstring[order] = '\0';
  number = digits;
  for (int i=order;i--;) {
    returnstring[i] = '0' + (char)(number % 10);
    number = (int)floor(((double)number / 10.));
  }
  //Log() << Verbose(0) << returnstring << endl;
  return returnstring;
};

/** Tests whether a given string contains a valid number or not.
 * \param *string
 * \return true - is a number, false - is not a valid number
 */
bool IsValidNumber( const char *string)
{
  int ptr = 0;
  if ((string[ptr] == '.') || (string[ptr] == '-')) // number may be negative or start with dot
    ptr++;
  if ((string[ptr] >= '0') && (string[ptr] <= '9'))
    return true;
  return false;
};

/** Blows the 6-dimensional \a cell_size array up to a full NDIM by NDIM matrix.
 * \param *symm 6-dim array of unique symmetric matrix components
 * \return allocated NDIM*NDIM array with the symmetric matrix
 */
double * ReturnFullMatrixforSymmetric(const double * const symm)
{
  double *matrix = new double[NDIM * NDIM];
  matrix[0] = symm[0];
  matrix[1] = symm[1];
  matrix[2] = symm[3];
  matrix[3] = symm[1];
  matrix[4] = symm[2];
  matrix[5] = symm[4];
  matrix[6] = symm[3];
  matrix[7] = symm[4];
  matrix[8] = symm[5];
  return matrix;
};

/** Calculate the inverse of a 3x3 matrix.
 * \param *matrix NDIM_NDIM array
 */
double * InverseMatrix( const double * const A)
{
  double *B = new double[NDIM * NDIM];
  double detA = RDET3(A);
  double detAReci;

  for (int i=0;i<NDIM*NDIM;++i)
    B[i] = 0.;
  // calculate the inverse B
  if (fabs(detA) > MYEPSILON) {;  // RDET3(A) yields precisely zero if A irregular
    detAReci = 1./detA;
    B[0] =  detAReci*RDET2(A[4],A[5],A[7],A[8]);    // A_11
    B[1] = -detAReci*RDET2(A[1],A[2],A[7],A[8]);    // A_12
    B[2] =  detAReci*RDET2(A[1],A[2],A[4],A[5]);    // A_13
    B[3] = -detAReci*RDET2(A[3],A[5],A[6],A[8]);    // A_21
    B[4] =  detAReci*RDET2(A[0],A[2],A[6],A[8]);    // A_22
    B[5] = -detAReci*RDET2(A[0],A[2],A[3],A[5]);    // A_23
    B[6] =  detAReci*RDET2(A[3],A[4],A[6],A[7]);    // A_31
    B[7] = -detAReci*RDET2(A[0],A[1],A[6],A[7]);    // A_32
    B[8] =  detAReci*RDET2(A[0],A[1],A[3],A[4]);    // A_33
  }
  return B;
};



/** Comparison function for GSL heapsort on distances in two molecules.
 * \param *a
 * \param *b
 * \return <0, \a *a less than \a *b, ==0 if equal, >0 \a *a greater than \a *b
 */
int CompareDoubles (const void * a, const void * b)
{
  if (*(double *)a > *(double *)b)
    return -1;
  else if (*(double *)a < *(double *)b)
    return 1;
  else
    return 0;
};


/**
 * Calls exit(255).
 */
void performCriticalExit() {
  exit(255);
}