utilities.h

/*
 * utilities.h
 *
 */
 
#ifndef UTILITIES_H_
#define UTILITIES_H_
 
#include <vector>
#include <math.h>
#include <algorithm>
#include <functional>
#include <string>
#include <stdexcept>
#include <iostream>
#include <sstream>
 
 
#ifndef PosType_declare
#define PosType_declare
typedef double PosType;
typedef unsigned long ULONG;
#endif
 
template <typename T> int sgn(T val) {
  return (T(0) < val) - (val < T(0));
}
 
template <class T>
inline T min(T x,T y){
  return (x < y) ? x : y;
};
template <class T>
inline T max(T x,T y){
  return (x > y) ? x : y;
};
 
 
namespace Utilities{
  
  inline void print_date(){
    time_t now = time(0);
    struct tm date = *localtime(&now);
    std::cout << date.tm_hour << ":" << date.tm_min << "   " << date.tm_mday << "/" << date.tm_mon << "/" << date.tm_year + 1900 << std::endl;
  }
 
  
  const double nXbin=64.;
  
  template <class T>
  void fill_linear ( std::vector<T> &v, size_t n, T min, T max ){
    v.resize ( n );
    for ( size_t i = 0; i < n; i++ )
      v[i] = min + (max - min) * T (i)/T (n-1);
  }
  
  template <class T>
  void fill_logarithmic ( std::vector<T> &v, size_t n, T min, T max ){
    v.resize (n);
    for ( size_t i = 0; i < n; i++ )
      v[i] = exp ( log (min) + ( log (max) - log (min) ) * T (i)/T (n-1) );
  }
  
  template <class T>
  long locate (const std::vector<T> &v, const T x)
  {
    long n = v.size ();
    if(n==0) return -1;
    
    if (x == v[0])
      return 0;
    if (x == v[n-1])
      return n-1;
    
    long jl = -1;
    long ju = n;
    bool as = (v[n-1] >= v[0]);
    while (ju-jl > 1)
    {
      long jm = (ju+jl)/2;
      if ((x >= v[jm]) == as)
        jl=jm;
      else
        ju=jm;
    }
    
    return jl;
  }
  template <class T,class F>
  long locate (const std::vector<T> &v,F x,std::function<bool (F ,const T &)> less_than)
  {
    long n = v.size ();
    
    if (less_than(x,v[0]))
      return -1;
    if (!less_than(x,v[n-1]))
      return n-1;
    
    long jl = -1;
    long ju = n;
    
    
    while (ju-jl > 1)
    {
      long jm = (ju+jl)/2;
      if ( less_than(x,v[jm]) )
        ju=jm;
      else
        jl=jm;
    }
    
    return jl;
  }
  
  template <class T>
  void locate (T *v, unsigned long n, T x, unsigned long *index)
  {
    if (x == v[0]){
      *index = 0;
      return;
    }
    if (x == v[n-1]){
      *index = n-1;
      return;
    }
    
    long jl = -1;
    long ju = n;
    bool as = (v[n-1] >= v[0]);
    while (ju-jl > 1)
    {
      long jm = (ju+jl)/2;
      if ((x >= v[jm]) == as)
        jl=jm;
      else
        ju=jm;
    }
    
    *index = jl;
    return;
  }
    template<typename T>
    size_t locate(const std::vector<T> &v                      
                  ,const std::vector<size_t> &sorted_index     
                  ,T value                                     
                  ,size_t &rank                                
    ){
        size_t n = v.size();
        
        if(n<=1) return 0;
        
        if(n != sorted_index.size() ){
            std::cerr << "sorted_index is the wrong size." << std::endl;
            throw std::invalid_argument("wrong size");
        }
        
        long imin=0,imax=v.size()-1;
        double fmax = v[sorted_index[imax]];
        double fmin = v[sorted_index[imin]];
        if(fmin > fmax){
            std::swap(imin,imax);
            std::swap(fmin,fmax);
        }
        
        if(value <= fmin) return sorted_index[imin];
        if(value >= fmax) return sorted_index[imax];
        
        size_t icurrent = (imax+imin)/2;
        double fcurrent = v[sorted_index[icurrent]];
        
        while ( abs(imax-imin) > 1 ) {
            if(fcurrent > value){
                fmax = fcurrent;
                imax = icurrent;
            }else{
                fmin = fcurrent;
                imin = icurrent;
            }
            icurrent = (imax+imin)/2;
            fcurrent = v[sorted_index[icurrent]];
        }
        
        rank=imin;
        return sorted_index[imin];
    }
 
  template <class T>
  size_t closest (const std::vector<T> &v, const T x)
  {
    
    int index = locate(v,x);
    if(index == -1) return 0;
    if(index == v.size()-1) return index;
    
    return (x - v[index]) < (v[index+1] - x) ? index : index+1;
  }
  
  
  PosType InterpolateYvec(std:: vector<PosType> &x, std:: vector<PosType> &y,PosType xi);
  
  PosType arctanh(PosType x);
  
  PosType fmini(PosType a,PosType b);
  
  PosType fmaxi(PosType a,PosType b);
  
  template <typename T>
  T median(std::vector<T> vec){
    size_t size = vec.size();
    if (size == 0){
      // std:: cout << "median of an empty vector" << std:: endl;
      return 0;
    }
    sort(vec.begin(), vec.end());
    size_t mid = size/2;
    return size % 2 == 0 ? (vec[mid] + vec[mid-1]) / 2 : vec[mid];
  }
  
  template <typename T>
  void range(std::vector<T> vec,T &max,T &min){
    if(vec.size() == 0){
      max = min = NAN;
      return;
    }
    max = min = vec[0];
    for(T a : vec){
      if(a > max) max = a;
      if(a < min) min = a;
    }
  }
 
  
  template <typename T>
  void polintT(T xa[], T ya[], int n, T x, T *y, T *dy)
  {
    int i,m,ns=1;
    T den,dif,dift,ho,hp,w;
    T *c,*d;
    
    dif=fabs(x-xa[1]);
    
    c = new T[n+1];
    d = new T[n+1];
    
    for (i=1;i<=n;i++) {
      if ( (dift=fabs(x-xa[i])) < dif) {
        ns=i;
        dif=dift;
      }
      c[i]=ya[i];
      d[i]=ya[i];
    }
    *y=ya[ns--];
    for (m=1;m<n;m++) {
      for (i=1;i<=n-m;i++) {
        ho=xa[i]-x;
        hp=xa[i+m]-x;
        w=c[i+1]-d[i];
        if ( (den=ho-hp) == 0.0) throw std::runtime_error("Error in routine polint");
        den=w/den;
        d[i]=hp*den;
        c[i]=ho*den;
      }
      *y += (*dy=(2*ns < (n-m) ? c[ns+1] : d[ns--]));
    }
    delete[] d;
    delete[] c;
  }
 
  template <typename FunctorType,typename T = double>
  T trapz(FunctorType &func, T a, T b, int n, T *s2)
  {
    T x,tnm,del,sum;
    int it,j;
    
    if (n == 1) {
      return (*s2=0.5*(b-a)*( func(a) + func(b) ));
    } else {
      
      for (it=1,j=1;j<n-1;j++) it <<= 1;
      tnm=it;
      del=(b-a)/tnm;
      x=a+0.5*del;
      for (sum=0.0,j=1;j<=it;j++,x+=del) sum += func(x);
      *s2=0.5*(*s2+(b-a)*sum/tnm);
      
      return *s2;
    }
  }
 
  template <typename FunctorType,typename T = double>
  T nintegrateF(
                FunctorType func        
                ,T a      
                ,T b      
                ,T tols   
  )
  {
    const int JMAX = 34,K=6;
    const int JMAXP = JMAX+1;
    
    if(a==b) return 0;
    
    T ss,dss;
    T s2[JMAXP],h2[JMAXP+1];
    T ss2=0;
    
    h2[1]=1.0;
    for (int j=1;j<=JMAX;j++) {
      s2[j] = Utilities::trapz<FunctorType,T>(func,a,b,j,&ss2);
       
      if (j>=K) {
        Utilities::polintT<T>(&h2[j-K],&s2[j-K],K,0.0,&ss,&dss);
        if(fabs(dss) <= tols*fabs(ss)) return ss;
      }
      h2[j+1]=0.25*h2[j];
    }
    std::cout << "s2= "; for(int j=1;j<=JMAX;j++) std::cout << s2[j] << "  ";
    std::cout << std::endl << "Too many steps in routine nintegrate<>\n";
    throw std::invalid_argument("Too many steps in routine nintegrate<>");
    return 0.0;
  }
 
  template <typename FunctorType,typename T = double>
  T nintegrate(
               FunctorType &func        
               ,T a      
               ,T b      
               ,T tols   
  )
  {
    const int JMAX = 34,K=6;
    const int JMAXP = JMAX+1;
    
    T ss,dss;
    T s2[JMAXP],h2[JMAXP+1];
    T ss2=0;
    
    h2[1]=1.0;
    for (int j=1;j<=JMAX;j++) {
      s2[j] = Utilities::trapz<FunctorType,T>(func,a,b,j,&ss2);
      if (j>=K) {
        Utilities::polintT<T>(&h2[j-K],&s2[j-K],K,0.0,&ss,&dss);
        if(fabs(dss) <= tols*fabs(ss)) return ss;
      }
      h2[j+1]=0.25*h2[j];
    }
    std::cout << "s2= "; for(int j=1;j<=JMAX;j++) std::cout << s2[j] << "  ";
    std::cout << std::endl << "Too many steps in routine nintegrate<>\n";
    throw std::invalid_argument("Too many steps in routine nintegrate<>");
    return 0.0;
  }
 
  template <class F,typename T>
  double bisection_search(
    F &func      
    ,T fvalue    
    ,T x1       
    ,T x2
    ,T xacc    
    ){
    int Nmax = 100000000,count=0;
    
    T f1 = func(x1) - fvalue;
    T f2 = func(x2) - fvalue;
    if( f1/f2 > 0 ){
      throw std::invalid_argument("not bracketed");
    }
    while( fabs(x1-x2) > xacc && count < Nmax ){
      ++count;
      T xnew = (x1 + x2)/2;
      T fnew = func(xnew) - fvalue;
      if( fnew/f1 < 0 ){
        x2 = xnew;
        f2 = fnew;
      }else{
        x1 = xnew;
        f1 = fnew;
      }
    }
    
    return (x1+x2)/2;
  }
  
  namespace stats{
  //**********************************************************
  //*** some routines for statistics of vectors **************
  //**********************************************************
  template <typename T>
  void vector_maxmin(const std::vector<T> &v,T &max,T &min ){
 
      max = min = v[0];
    for(T d : v){
      max = (max > d) ? max : d;
      min = (min < d) ? min : d;
    }
  }
  
  template <typename T>
  T vector_mean(const std::vector<T> &v){
    T ans = 0;
    for(T a : v) ans += a;
    return ans/v.size();
  }
  
  template <typename T>
  T vector_variance(const std::vector<T> &v){
    T m = vector_mean(v);
    T ans = 0;
    for(T a : v) ans += (a - m)*(a - m);
    return ans/(v.size() - 1);
  }
  
  // Pearson's correlation coefficient
  template <typename T>
  T vector_correlation(const std::vector<T> &v1,const std::vector<T> &v2){
    size_t N = v1.size();
    if(v2.size() != N){
      throw std::invalid_argument("vectors are not the same size");
    }
    
    T m1 = vector_mean(v1);
    T m2 = vector_mean(v2);
    T var1 = vector_variance(v1);
    T var2 = vector_variance(v2);
    
    T ans = 0;
    for(size_t i = 0 ; i < N ; ++i){
      ans += (v1[i]-m1)*(v2[i]-m2);
    }
    ans /= (N-1)*sqrt(var1*var2);
    
    return ans;
  }
  
  template <typename T>
  class Accumulater{
  public:
    Accumulater(int size):
    sum(0),count(0),N(size)
    {
      products.resize(N*N,0);
      sums.resize(N,0);
    }
    
    void add(const std::vector<T> &v){
      if(v.size() != N ) throw std::invalid_argument("wrong size vector");
      for(T a : v) sum += a;
      
      for(size_t i = 0 ; i < N ; ++i){
        sums[i] += v[i];
        for(size_t j = 0 ; j < N ; ++j){
          products[j + i*N] += v[i]*v[j];
        }
      }
      ++count;
    }
 
    void add(T a){
      if(1 != N ) throw std::invalid_argument("wrong size vector");
      sum += a;
      
      sums[0] += a;
      products[0] += a*a;
      ++count;
    }
 
    T mean(){
      T sum = 0;
      for(auto d : sums) sum += d;
      return sum/count/N;
    }
    T mean(size_t i){return sums[i]/count;}
    void means(std::vector<T> &mm){
      mm = sums;
      for(T &m : mm) m /= count;
    }
    
    T covariance(size_t i,size_t j){
      return (products[j + i*N] - sums[i]*sums[j]/count)/(count - 1);
    }
    
    void variances(std::vector<T> &mm){
      mm.resize(N);
      for(size_t i=0 ; i < N ; ++i){
        mm[i] = covariance(i,i);
      }
    }
    
    T ave_variances(){
      T ans = 0;
      for(size_t i=0 ; i < N ; ++i){
        ans += covariance(i,i);
      }
      return ans/N;
    }
    
    void cov_matrix(std::vector<T> &mm){
      mm.resize(N*N);
      for(size_t i = 0 ; i<N ; ++i){
        T avi = sums[i]/count;
        for(size_t j = i ; j<N ; ++j){
          size_t k = i + j*N;
          mm[k] = (products[k] - avi*sums[j])/(count - 1);
          mm[j + i*N] = mm[k];
        }
      }
    }
    
    void cov_matrix_normalized(std::vector<T> &mm){
      std::vector<T> var;
      variances(var);
      
      mm.resize(N*N);
      for(size_t i = 0 ; i<N ; ++i){
        T avi = sums[i]/count;
        for(size_t j = i ; j<N ; ++j){
          size_t k = i + j*N;
          mm[k] = (products[k] - avi*sums[j])/(count - 1)/sqrt( var[i]*var[j] );
          mm[j + i*N] = mm[k];
        }
      }
    }
    
  private:
    T sum;
    std::vector<T> products;
    std::vector<T> sums;
    size_t count;
    size_t N;
  };
 
  }
}
#endif /* UTILITIES_H_ */