Helpers.h

/** \file Helpers.h */

#ifndef MY_HELPERS_H
#define MY_HELPERS_H

#include "AMDiS.h"

#if HAVE_PARALLEL_DOMAIN_AMDIS
#include "parallel/StdMpi.h"
#endif

#ifdef _WIN32
#include <io.h>
#else
#include <unistd.h>
#endif
#include <ios>
#include <iostream>
#include <fstream>
#include <string>

#include <boost/numeric/mtl/mtl.hpp>
#include <boost/numeric/itl/itl.hpp>
#include "boost/lexical_cast.hpp"
#include "boost/date_time/posix_time/posix_time.hpp"
#include <boost/math/special_functions/fpclassify.hpp>
#include <boost/math/special_functions/atanh.hpp>

#include <time.h>
#include <stdarg.h>

#include "VectorOperations.h"
#include "Views.h"

#define TEST_WARNING(test) if ((test));else WARNING

using namespace std;
using namespace AMDiS; 

static long random_seed_initial_value = 0;

namespace Helpers {

  // math routines
  // ==============
  
  inline double cint(double x) {
    double intpart; 
    if (modf(x, &intpart) >= 0.5)
      return x >=0 ? ceil(x) : floor(x);
    else
      return x < 0 ? ceil(x) : floor(x);
  }
  
  inline double round (double r, double places) {
    double off = pow(10.0, places);
    return cint(r*off)/off;
  }
  
  inline double signum(double val, double tol = DBL_TOL) { return val > tol ? 1.0 : (val < -tol ? -1.0 : 0.0); }
  inline int signum(int val) { return val > 0 ? 1 : (val < 0 ? -1 : 0); }
  inline int signum(bool val) { return val ? 1 : -1; }
  inline int toInt(bool val) { return val ? 1 : 0; }
  inline bool isNumber(double val) { return !boost::math::isnan(val) && !boost::math::isinf(val); }
  inline bool toBool(double val, double tol = DBL_TOL) { return val > tol || val < -tol ? true : false; }
  inline bool isPositiv(double val, double tol = DBL_TOL) { return val > -tol; }
  inline bool isStrictPositiv(double val, double tol = DBL_TOL) { return val > tol; }

  template<typename T>
  T atanh(T x)
  {
    T result;
    try {
      result = boost::math::atanh(x);
    } catch(...) { // domain_error, or overflow_error
      if (x > 0.0)
	result = 10.0;
      else
	result = -10.0;
    }
    return result;
  }
  
  // routines for conversion to string
  // =================================
  
  template<typename T>
  inline std::string toString(const T &value, 
			      ios_base::fmtflags flag = ios_base::scientific)
  {
    std::stringstream ss;
    ss.setf(flag);
    ss.precision(6);
    if (!(ss<<value))
      throw(std::runtime_error("toString() not possible"));
    return ss.str();
  }
  
  inline std::string toString(const int &value)
  {
    return boost::lexical_cast<std::string>(value);
  }
  
  template<typename T>
  inline std::string toString(const WorldVector<T> &value, 
			      bool brackets = true, 
			      ios_base::fmtflags flag = ios_base::scientific)
  {
    std::stringstream ss;
    if (brackets) ss<<"[";
    ss.setf(flag);
    ss.precision(6);
    if (!(ss<<value))
      throw(std::runtime_error("toString() not possible"));
    
    if (brackets) ss<<"]";
    return ss.str();
  }
  
  template<typename T>
  inline std::string toString(const std::vector<T> &vec, 
			      bool brackets = true, 
			      ios_base::fmtflags flag = ios_base::scientific)
  {
    std::stringstream ss;
    if (brackets) ss<<"[";
    ss.setf(flag);
    ss.precision(6);
    for (size_t i = 0; i < vec.size(); ++i) {
      if (!(ss<<(i>0?", ":"")<<vec[i]))
	throw(std::runtime_error("toString() not possible"));
    }
    if (brackets) ss<<"]";
    return ss.str();
  }
  
  template<class T>
  inline T fromString(const std::string& s)
  {
    std::istringstream stream(s);
    T t;
    stream >> t;
    return t;
  }

  inline std::string fillString(int length, char c, int numValues, ...)
  {
    va_list values;
    int value;
    va_start(values, numValues);

    int len = length;
    for(int i=0; i<numValues; i++) {
      value = va_arg(values, int);
      double tmp = static_cast<double>(value);
      while (true) {
	len --;
	tmp /= 10.0;
	if (tmp < 1.0 - DBL_TOL)
	  break;
      }
    }
    va_end(values);

    return std::string(len, c);
  }

  /// produce printable string of memory
  inline std::string memoryToString(double mem, int startIdx = 0) {
    int idx = startIdx;
    double mem_ = mem;
    while (mem_/1024.0 > 1.0) {
      mem_ /= 1024.0;
      idx++;
    }
    std::string result = toString(mem_)+" "+(idx==0?"B":(idx==1?"KB":(idx==2?"MB":"GB")));
    return result;
  }
  
  
  // some mesh routines
  // ===========================
  
  /// rotate scale and shift the mesh coordinates
  void transformMesh(Mesh *mesh, WorldVector<double> scale, WorldVector<double> shift, WorldVector<double> rotate);

  /// scale a vector of meshes by different values in all directions
  void scaleMesh(std::vector<Mesh*> meshes, WorldVector<double> scale);
  
  /// scale and shift by different values in all directions
  void scaleMesh(Mesh *mesh, WorldVector<double> shift, WorldVector<double> scale);
  
  /// scale by different values in all directions
  void scaleMesh(Mesh *mesh, WorldVector<double> scale);
  
  /// scale and shift by the same values in all directions
  void scaleMesh(Mesh *mesh, double shift=0.0, double scale=1.0);
  
  /**
   * calculate the dimension of a mesh, by mesh traversal.
   * The mthod determines the maximal mesh coordinates and assumes symmetry of the
   * mesh around the origin, i.e. mesh = [-rsult, result]
   **/
  WorldVector<double> getMeshDimension(Mesh *mesh);
  
  /// calculate the dimension of a mesh, by mesh traversal.
  void getMeshDimension(Mesh *mesh, WorldVector<double> &min_corner, WorldVector<double> &max_corner);
  
  /// read DOFVector from AMDiS .dat-files
  void read_dof_vector(const std::string file, DOFVector<double> *dofvec, long size);

  // some linear algebra methods
  // ===========================

  /// calculate the determant of a 3x3 matrix of type dense2D
  inline double determinant(mtl::dense2D<double> &A) // only for 3x3 - matrices
  {
    if(num_rows(A)==3 && num_cols(A)==3) {
      double det = A(0,0)*A(1,1)*A(2,2)+A(0,1)*A(1,2)*A(2,0)+A(0,2)*A(1,0)*A(2,1);
      return det-(A(0,2)*A(1,1)*A(2,0)+A(0,1)*A(1,0)*A(2,2)+A(0,0)*A(1,2)*A(2,1));
    }
    return 1.0;
  }
  
  /// calculate the determant of a 3x3 matrix of type WorldMatrix
  inline double determinant(WorldMatrix<double> &A) // only for 3x3 - matrices
  {
    if(A.getNumRows()==3 && A.getNumCols()==3) {
      double det = A[0][0]*A[1][1]*A[2][2]+A[0][1]*A[1][2]*A[2][0]+A[0][2]*A[1][0]*A[2][1];
      return det-(A[0][2]*A[1][1]*A[2][0]+A[0][1]*A[1][0]*A[2][2]+A[0][0]*A[1][2]*A[2][1]);
    }
    return 1.0;
  }
  
  /// inverse power method to find minimal eigenvalue of matrix A and appropriate eigenvector x, using m iterations
  template <class LinearSolver, class Matrix, class EigenVector>
  inline double inverse_iteration(LinearSolver& solver, const Matrix& A, EigenVector& x, int m)
  { FUNCNAME("Helpers::inverse_iteration()");

    EigenVector y( size(x) ), res( size(x) );
    double lambda = 0.0, lambda_old = 0.0, relErr;
    mtl::seed<double> seed;
    random(x, seed);
    x /= two_norm(x);

    solver.setMultipleRhs(true);
    for (int i = 0; i < m; ++i) {
      solver.solveSystem(A, y, x ); // solve Ay=x --> y
      lambda_old = lambda;
      lambda = dot(y,x) / dot(y,y);
      relErr = abs((lambda - lambda_old) / lambda);
      if (relErr < 1.e-5)
	break;
      x = y / two_norm(y);
    }
    solver.setMultipleRhs(false);

    return lambda;
  }

  /// calculate approximation to condition number of matrix A using the OEMSolver solver
  template <class LinearSolver, class Matrix>
  inline double condition (LinearSolver& solver, const Matrix& A, int m=10)
  { FUNCNAME("Helpers::condition()");

    mtl::dense_vector<double>   x(num_rows(A)), y(num_rows(A));
    mtl::seed<double>           seed;
    random(x, seed);

    double lambda = 0.0, lambda_old = 0.0, relErr, result1;
    x /= two_norm(x);
    for (int i = 0; i < m; ++i) {
      y = A * x;
      lambda_old = lambda;
      lambda = mtl::dot(y,x);
      relErr = abs((lambda-lambda_old)/lambda);
      if (relErr < 1.e-5)
	break;
      x = y / two_norm(y) ;
    }
    result1 = std::abs(lambda / inverse_iteration(solver, A, x, m));
    return result1;
  }
  
  /// interpolate values of DOFVector along line from (x1,y1) -> (x2,y1) with nPoints points
  template<typename Container>
  inline void interpolOverLine(const Container &vec, 
			       double x1, double y1, double x2, double y2, 
			       int nPoints,
			       std::vector<WorldVector<double> > &x, std::vector<double> &y)
  { FUNCNAME("Helpers::interpolOverLine()");

    WorldVector<double> p;

    if (nPoints <= 1)
      throw(std::runtime_error("Zu wenig Punkte fuer eine Interpolation!"));
    
    for (int i = 0; i < nPoints; ++i) {
      double lambda = static_cast<double>(i)/static_cast<double>(nPoints - 1.0);
      p[0] = lambda*x2 + (1.0-lambda)*x1;
      p[1] = lambda*y2 + (1.0-lambda)*y1;
      
      double value = evalAtPoint(vec, p);
      x.push_back(p);
      y.push_back(value);
    }
  }

  /// calculate maxima of DOFVector along x-axis, using interpolOverLine
  static void calcMaxOnXAxis(DOFVector<double> *rho, std::vector<std::pair<WorldVector<double>, double> > &maxima);
  
  /// calculate maxima of DOFVector along y-axis, using interpolOverLine
  static void calcMaxOnYAxis(DOFVector<double> *rho, std::vector<std::pair<WorldVector<double>, double> > &maxima);

  /// calc normal vectors of surface from element normals by averaging
  static void getNormalsWeighted(FiniteElemSpace *feSpace, DOFVector<WorldVector<double> > *normals);
  
  /// calc normal vectors of surface from element normals by local approximation by quartic
  static void getNormals(FiniteElemSpace *feSpace, DOFVector<WorldVector<double> > *normals, DOFVector<WorldMatrix<double> > *gradNormals = NULL);
  
  /// calc curvature from given normal vectors
  static void getCurvature(DOFVector<WorldVector<double> >* normals, DOFVector<double>* curvature);
  
  // misc routines
  // =============
  
  void plot(std::vector<double> &values, int numRows = 7, int numCols = 20, std::string symbol = "*");
  
  inline long getMicroTimestamp() {
    using namespace boost::posix_time;
    ptime t0(min_date_time);
    ptime now = microsec_clock::local_time();
    time_duration td = now-t0;
    return td.total_microseconds();
  }
  
  inline long getRandomSeed(int param = 0) {
    using namespace boost::posix_time;
    ptime t0(min_date_time);
    ptime now = microsec_clock::local_time();
    time_duration td = now-t0;
    long value0 = td.total_microseconds();
    long value1 = clock();
    long value2 = random_seed_initial_value++;

    return value0 + value1 + value2 + param;
  }

  
} // end namespace Helpers

namespace AMDiS {

  /** \ingroup DOFAdministration
    * \brief
    * Implements a DOFIterator for a DOFIndexed<T> object
    */
  template<typename T>
  class DOFConstIterator : public DOFIteratorBase
  {
  public:
    /// Constructs a DOFIterator for cont of type t
    DOFConstIterator(const DOFIndexed<T> *obj, DOFIteratorType t)
      : DOFIteratorBase(dynamic_cast<DOFAdmin*>(obj->getFeSpace()->getAdmin()), t),
	iteratedObject(obj)
    {}

    /// Constructs a DOFIterator for cont of type t
    DOFConstIterator(DOFAdmin *admin,
		const DOFIndexed<T> *obj,
		DOFIteratorType t)
      : DOFIteratorBase(admin, t),
	iteratedObject(obj)
    {}

    /// Dereference operator
    inline const T& operator*()
    {
      return *it;
    }

    /// Dereference operator
    inline const T* operator->()
    {
      return &(*it);
    }

    inline bool operator!=(const DOFIterator<T>& rhs)
    {
      if (this->iteratedObject != rhs.iteratedObject)
	return true;

      if (this->it != rhs.it)
	return true;

      return false;
    }

    inline bool operator==(const DOFIterator<T>& rhs)
    {
      return !(this->operator==(rhs));
    }

  protected:
    /// Implementation of DOFIteratorBase::goToBeginOfIteratedObject()
    inline void goToBeginOfIteratedObject()
    {
      it = iteratedObject->begin();
    }

    /// Implementation of DOFIteratorBase::goToEndOfIteratedObject()
    inline void goToEndOfIteratedObject()
    {
      it = iteratedObject->end();
    }

    /// Implementation of DOFIteratorBase::incObjectIterator()
    inline void incObjectIterator()
    {
      ++it;
    }

    /// Implementation of DOFIteratorBase::incObjectIterator()
    inline void decObjectIterator()
    {
      --it;
    }

  protected:
    /// Object that is iterated
    const DOFIndexed<T> *iteratedObject;

    /// Iterator for \ref iteratedObject
    typename std::vector<T>::const_iterator it;
  };

} // end namespace AMDiS

#endif