source: src/helpers.cpp@ ffe885

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


#include "helpers.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 {
    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;
};

/** Returns the power of \a n with respect to \a base.
 * \param base basis
 * \param n power
 * \return \f$base^n\f$
 */
int pot(int base, int n)
{
  int res = 1;
  int j;
  for (j=n;j--;)
    res *= base;
  return res;
};

/** 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 = Malloc<char>(MAXSTRINGSIZE, "CountLinesinFile: *buffer");
  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);
  Free(&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 = Malloc<char>(order + 2, "FixedDigitNumber: *returnstring");
  // 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 = Malloc<double>(NDIM * NDIM, "molecule::ReturnFullMatrixforSymmetric: *matrix");
  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 = Malloc<double>(NDIM * NDIM, "Vector::InverseMatrix: *B");
  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;
};

/** hard-coded determinant of a 3x3 matrix.
 * \param a[9] matrix
 * \return \f$det(a)\f$
 */
double RDET3(const double a[NDIM*NDIM])
{
  return ((a)[0]*(a)[4]*(a)[8] + (a)[3]*(a)[7]*(a)[2] + (a)[6]*(a)[1]*(a)[5] - (a)[2]*(a)[4]*(a)[6] - (a)[5]*(a)[7]*(a)[0] - (a)[8]*(a)[1]*(a)[3]);
};

/** hard-coded determinant of a 2x2 matrix.
 * \param a[4] matrix
 * \return \f$det(a)\f$
 */
double RDET2(const double a[4])
{
  return ((a[0])*(a[3])-(a[1])*(a[2]));
};

/** hard-coded determinant of a 2x2 matrix.
 * \param a0 (0,0) entry of matrix
 * \param a1 (0,1) entry of matrix
 * \param a2 (1,0) entry of matrix
 * \param a3 (1,1) entry of matrix
 * \return \f$det(a)\f$
 */
double RDET2(const double a0, const double a1, const double a2, const double a3)
{
  return ((a0)*(a3)-(a1)*(a2));
};

/** 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;
};


/** Allocates a memory range using malloc().
 * Prints the provided error message in case of a failure.
 *
 * \param number of memory slices of type X to allocate
 * \param failure message which is printed if the allocation fails
 * \return pointer to the allocated memory range, will be NULL if a failure occurred
 */
template <> char* Malloc<char>(size_t size, const char* output)
{
  char* buffer = NULL;
  buffer = (char*) malloc(sizeof(char) * (size + 1));
  for (size_t i = size; i--;)
    buffer[i] = (i % 2 == 0) ? 'p': 'c';
  buffer[size] = '\0';

  if (buffer != NULL) {
    MemoryUsageObserver::getInstance()->addMemory(buffer, size);
  } else {
    Log() << Verbose(0) << "Malloc for datatype " << typeid(char).name()
      << " failed - pointer is NULL: " << output << endl;
  }

  return buffer;
};

/**
 * Frees all memory registered by the memory observer and calls exit(225) afterwards.
 */
void performCriticalExit() {
  map<void*, size_t> pointers = MemoryUsageObserver::getInstance()->getPointersToAllocatedMemory();
  for (map<void*, size_t>::iterator runner = pointers.begin(); runner != pointers.end(); runner++) {
    Free(((void**) &runner->first));
  }

  exit(255);
}