utilities_slsim.h

#ifndef _UTILITIES_H
#define _UTILITIES_H
 
#include "standard.h"
#include <typeinfo>
#include <vector>
#include <map>
#include <string>
#include <stdexcept>
#include <algorithm>
#include <utility>
#include <iterator>
#include <cstdlib>
#include <random>
//#if __cplusplus >= 201103L
#include <typeindex>
#include <dirent.h>
#include <unistd.h>
#include <sys/stat.h>
#include <sys/types.h>
#include <boost/variant.hpp>
#include <set>
#include <iomanip>
#include <thread>
#include <cmath>
#include <cstdint>
#include <limits>
#include <stdexcept>
#include <fstream>
#include <sstream>
 
namespace Utilities
{
template<class T>
T to_numeric(const std::string &str) {
  return std::stoi(str);
};
template<>
inline long to_numeric<long>(const std::string &str) {
  try{
    return std::stol(str);
  }catch(...){
    return -500;
  }
};
template<>
inline int to_numeric<int>(const std::string &str) {
  return std::stoi(str);
  try{
    return std::stoi(str);
  }catch(...){
    return -500;
  }
};
template<>
inline float to_numeric<float>(const std::string &str) {
  return std::stof(str);
  try{
    return std::stof(str);
  }catch(...){
    return -500;
  }
};
template<>
inline double to_numeric<double>(const std::string &str) {
  try{
    return std::stod(str);
  }catch(...){
    return -500;
  }
};
//********************************************************
 
template <typename T>
bool AlwaysTrue(T t){return true;}
 
template <typename T>
bool AlwaysFalse(T t){return false;}
  
template <typename T>
T vec_sum(const std::vector<T> &v){
  double sum=0.0;
  for(const T &a : v) sum += a;
  return sum;
}
 
// this is not for the user
namespace detail
{
//#if __cplusplus < 201103L
//class type_index
//{
//public:
//  type_index(const std::type_info& type) : t(type) {}
//  inline bool operator<(const type_index& rhs) const { return t.before(rhs.t); }
//
//private:
//  const std::type_info& t;
//};
//#else
using std::type_index;
//#endif
}
 
template <class T>
void Matrix(T **matrix, long rows, long cols){
  matrix = new T*[rows];
  matrix[0] = new T[rows*cols];
  for (long i = 1; i < rows; ++i)
    matrix[i] = matrix[0] + i * cols;
}
 
template <class T>
void free_Matrix(T **matrix, long rows, long){
  if (rows) delete[] matrix[0];
  delete []matrix;
}
 
PosType **PosTypeMatrix(size_t rows, size_t cols);
void free_PosTypeMatrix(PosType **matrix, size_t rows, size_t cols);
PosType **PosTypeMatrix(long rows1,long rows2, long cols1, long cols2);
void free_PosTypeMatrix(PosType **matrix, long rows1,long rows2, long cols1, long cols2);
 
template<typename T>
class D3Matrix{
public:
  D3Matrix(size_t xsize,size_t ysize,size_t zsize)
  :xn(xsize),yn(ysize),zn(zsize){
    array = new T[xn*yn*zn];
  }
  ~D3Matrix(){
    delete[] array;
  }
  
  T& operator()(size_t i,size_t j,size_t k){
    return array[i + j*xn + k*xn*yn];
  }
  
  T& operator()(size_t m){
    return array[m];
  }
  
  size_t xindex(size_t m){
    return m % (xn);
  }
  size_t yindex(size_t m){
    return (m % (xn*yn) ) / xn;
  }
  
  size_t zindex(size_t m){
    return m / (xn*yn);
  }
  
private:
  size_t xn;
  size_t yn;
  size_t zn;
  T *array;
};
 
 
template<typename T>
class D2Matrix{
  
public:
  D2Matrix(size_t xsize,size_t ysize)
  :xn(xsize),yn(ysize){
    array = new T[xn*yn];
  }
  ~D2Matrix(){
    delete[] array;
  }
  
  T& operator()(size_t i,size_t j){
    return array[i + j*xn];
  }
  
  T& operator()(size_t m){
    return array[m];
  }
  
  size_t xindex(size_t m){
    return m % xn;
  }
  size_t yindex(size_t m){
    return m / xn;
  }
  
  
private:
  size_t xn;
  size_t yn;
  T *array;
};
 
template <typename T>
class SymmetricMatrix{
  std::vector<T> v;
  int n;
  int m;
public:
  SymmetricMatrix(int n):n(n){
    v.resize(n*(n+1)/2);
    m = 2*n-1;
  }
  T& operator()(int i,int j){
    //long k = j + (2*n-1-i)*i/2 ;
    //size_t k = (i <= j ) ? j + (2*n-1-i)*i/2 : i + (2*n-1-j)*j/2;
    size_t k = (i <= j ) ? j + (m-i)*i/2 : i + (m-j)*j/2;
    //assert(k>=0);
    //assert(k < v.size());
    return v[ k ];
  }
  T& operator[](size_t k){
    return v[ k ];
  }
  const T &operator[](size_t k) const {
    return v[k];
  }
 
  int size(){return n;}
  
  size_t oned_index(int i,int j){
    return (i <= j ) ? j + (m-i)*i/2 : i + (m-j)*j/2;
  }
};
 
 
template<typename BaseT>
class MixedVector
{
public: /* iterators */
  template<typename ValueT = BaseT> // TODO: needs to check subclass type
  class iterator
  {
  private:
    typedef typename std::vector<BaseT*>::iterator base_iterator;
    base_iterator it;
    
  public:
    typedef ValueT value_type;
    typedef ValueT* pointer;
    typedef const ValueT* const_pointer;
    typedef ValueT& reference;
    typedef const ValueT& const_reference;
    typedef typename base_iterator::difference_type difference_type;
    typedef typename base_iterator::iterator_category iterator_category;
    
    iterator() {}
    iterator(base_iterator i) : it(i) {}
    iterator(const iterator& other) : it(other.it) {}
    
    iterator& operator=(const iterator& rhs) { it = rhs.it; return *this; }
    
    reference operator*() { return (reference)(**it); }
    
    //const reference operator*() const { return (const reference)(**it); }
    
    pointer operator->() { return (pointer)(*it); }
    const_pointer operator->() const { return (const_pointer)(*it); }
    
    iterator& operator++() { ++it; return *this; }
    iterator operator++(int) { iterator tmp(*this); ++it; return tmp; }
    iterator& operator--() { --it; return *this; }
    iterator operator--(int) { iterator tmp(*this); --it; return tmp; }
    
    bool operator==(const iterator& rhs) const { return (it == rhs.it); }
    bool operator!=(const iterator& rhs) const { return (it != rhs.it); }
    
    iterator& operator+=(difference_type n) { it += n; return *this; }
    iterator& operator-=(difference_type n) { it -= n; return *this; }
    
    reference operator[](difference_type n) { return (reference)*it[n]; }
    
    //const reference operator[](difference_type n) const { return (const reference)*it[n]; }
    
    friend iterator operator+(const iterator& i, difference_type n) { return iterator(i.it + n); }
    friend iterator operator+(difference_type n, const iterator& i) { return iterator(i.it + n); }
    friend iterator operator-(const iterator& i, difference_type n) { return iterator(i.it - n); }
    friend iterator operator-(difference_type n, const iterator& i) { return iterator(i.it - n); }
    
    friend difference_type operator-(const iterator& b, const iterator& a) { return (b.it - a.it); }
    
    friend bool operator<(const iterator&a, const iterator& b) { return (a.it < b.it); }
    friend bool operator>(const iterator&a, const iterator& b) { return (a.it > b.it); }
    friend bool operator<=(const iterator&a, const iterator& b) { return (a.it <= b.it); }
    friend bool operator>=(const iterator&a, const iterator& b) { return (a.it >= b.it); }
  };
  
public:
  MixedVector()
  {
  }
  
  ~MixedVector()
  {
    clear();
  }
  
  friend void swap(MixedVector<BaseT>& a, MixedVector<BaseT>& b)
  {
    using std::swap;
    
    swap(a.items, b.items);
    swap(a.tmap, b.tmap);
  }
  
  MixedVector<BaseT>& operator=(MixedVector<BaseT> rhs)
  {
    // copy-and-swap idiom
    swap(*this, rhs);
    return *this;
  }
  
  template<typename SubclassT>
  void copy_back(const SubclassT& obj)
  {
    // make sure this is a subclass of BaseT
    check_subclass(obj);
    
    // copy the object
    SubclassT* copy = new SubclassT(obj);
    
    // add the copy of the object to the list of items
    items.push_back(copy);
    
    // add the copy to type map
    tmap[typeid(SubclassT)].push_back(copy);
  }
 
 
  void pop_back()
  {
    // get very last element
    BaseT* back = items.back();
    
    // remove from the vector in the type map
    tmap[back->type()].pop_back();
    
    // remove from items
    items.pop_back();
    
    // delete from memory
    delete back;
  }
  
  template<typename SubclassT>
  void pop_back()
  {
    // search for vector belonging to SubclassT
    typename tmap_t::iterator found = tmap.find(typeid(SubclassT));
    if(found == tmap.end())
      return;
    
    // pop from vector for type
    found->second.pop_back();
    
    // go through list of all items from back, until item of correct type is found
    for(typename base_container::reverse_iterator it = items.rbegin(); it != items.rend(); ++it)
    {
      if((*it)->type() == typeid(SubclassT))
      {
        // remove the item
        delete (*it);
        items.erase(it);
        return;
      }
    }
  }
  
  void clear()
  {
    for(std::size_t i = 0, n = items.size(); i < n; ++i)
      delete items[i];
    items.clear();
    tmap.clear();
  }
  
  template<typename SubclassT>
  void clear()
  {
    // search for vector belonging to SubclassT
    typename tmap_t::iterator found = tmap.find(typeid(SubclassT));
    if(found == tmap.end())
      return;
    
    base_container& titems = found.second;
    
    // delete items from vector for subclass
    for(std::size_t i = 0, n = titems.size(); i < n; ++i)
      delete titems[i];
    
    // remove items for subclass
    tmap.erase(found);
    
    // erase items matching type
    items.erase(std::remove_if(items.begin(), items.end(), is_type<SubclassT>), items.end());
  }
  
  BaseT* data()
  {
    return items.data();
  }
  
  const BaseT* data() const
  {
    return items.data();
  }
  
  template<typename SubclassT>
  bool type(std::size_t i)
  {
    return is_type<SubclassT>(items[i]);
  }
  
  BaseT& operator[](std::size_t i)
  {
    return *items[i];
  }
  
  const BaseT& operator[](std::size_t i) const
  {
    return *items[i];
  }
  
  BaseT& get(std::size_t i)
  {
    return *items[i];
  }
  
  const BaseT& get(std::size_t i) const
  {
    return *items[i];
  }
  
  template<typename SubclassT>
  SubclassT& get(std::size_t i)
  {
    typename tmap_t::iterator found = tmap.find(typeid(SubclassT));
    if(found == tmap.end())
      throw std::out_of_range(std::string() + "type " + typeid(SubclassT).name() + " not in vector");
    return (SubclassT&)*found->second[i];
  }
  
  template<typename SubclassT>
  const SubclassT& get(std::size_t i) const
  {
    typename tmap_t::const_iterator found = tmap.find(typeid(SubclassT));
    if(found == tmap.end())
      throw std::out_of_range(std::string() + "type " + typeid(SubclassT).name() + " not in vector");
    return (const SubclassT&)*found->second[i];
  }
  
  BaseT& at(std::size_t i)
  {
    return *items.at(i);
  }
  
  const BaseT& at(std::size_t i) const
  {
    return *items.at(i);
  }
  
  template<typename SubclassT>
  SubclassT& at(std::size_t i)
  {
    typename tmap_t::iterator found = tmap.find(typeid(SubclassT));
    if(found == tmap.end())
      throw std::out_of_range(std::string() + "type " + typeid(SubclassT).name() + " not in vector");
    return (SubclassT&)*found->second.at(i);
  }
  
  template<typename SubclassT>
  const SubclassT& at(std::size_t i) const
  {
    typename tmap_t::const_iterator found = tmap.find(typeid(SubclassT));
    if(found == tmap.end())
      throw std::out_of_range(std::string() + "type " + typeid(SubclassT).name() + " not in vector");
    return (const SubclassT&)*found->second.at(i);
  }
  
  std::size_t size() const
  {
    return items.size();
  }
  
  template<typename SubclassT>
  std::size_t size() const
  {
    typename tmap_t::const_iterator found = tmap.find(typeid(SubclassT));
    if(found == tmap.end())
      return 0;
    return found->second.size();
  }
  
  bool empty() const
  {
    return items.empty();
  }
  
  template<typename SubclassT>
  bool empty() const
  {
    typename tmap_t::const_iterator found = tmap.find(typeid(SubclassT));
    if(found == tmap.end())
      return true;
    
    return found->second.empty();
  }
  
  iterator<> begin()
  {
    return iterator<>(items.begin());
  }
  
  iterator<> end()
  {
    return iterator<>(items.end());
  }
  
  template<typename SubclassT>
  iterator<SubclassT> begin()
  {
    typename tmap_t::iterator found = tmap.find(typeid(SubclassT));
    if(found == tmap.end())
      return iterator<SubclassT>(items.end());
    return iterator<SubclassT>(found->second.begin());
  }
  
  template<typename SubclassT>
  iterator<SubclassT> end()
  {
    typename tmap_t::iterator found = tmap.find(typeid(SubclassT));
    if(found == tmap.end())
      return iterator<SubclassT>(items.end());
    return iterator<SubclassT>(found->second.end());
  }
  
private:
  MixedVector(const MixedVector<BaseT>& other);
  
  typedef std::vector<BaseT*> base_container;
  
  inline static void check_subclass(const BaseT&) {}
  
  template<typename SubclassT>
  inline static bool is_type(const BaseT* obj) { return (typeid(*obj) == typeid(SubclassT)); }
  
  base_container items;
  
  typedef std::map<detail::type_index, base_container> tmap_t;
  tmap_t tmap;
};
 
template<typename BaseT>
class MixedVector<BaseT*>
{
public: /* iterators */
  template<typename ValueT = BaseT*> // TODO: needs to check subclass pointer type
  class iterator
  {
  private:
    typedef typename std::vector<BaseT*>::iterator base_iterator;
    base_iterator it;
    
  public:
    typedef typename std::vector<ValueT>::iterator::value_type value_type;
    typedef typename std::vector<ValueT>::iterator::pointer pointer;
    typedef typename std::vector<ValueT>::iterator::reference reference;
    typedef typename base_iterator::difference_type difference_type;
    typedef typename base_iterator::iterator_category iterator_category;
    
    iterator() {}
    iterator(base_iterator i) : it(i) {}
    iterator(const iterator& other) : it(other.it) {}
    
    iterator& operator=(const iterator& rhs) { it = rhs.it; return *this; }
    
    reference operator*() { return (reference)*it; }
    const reference operator*() const { return (const reference)*it; }
    
    iterator& operator++() { ++it; return *this; }
    iterator operator++(int) { iterator tmp(*this); ++it; return tmp; }
    iterator& operator--() { --it; return *this; }
    iterator operator--(int) { iterator tmp(*this); --it; return tmp; }
    
    bool operator==(const iterator& rhs) const { return (it == rhs.it); }
    bool operator!=(const iterator& rhs) const { return (it != rhs.it); }
    
    iterator& operator+=(difference_type n) { it += n; return *this; }
    iterator& operator-=(difference_type n) { it -= n; return *this; }
    
    pointer operator[](difference_type n) { return (pointer)it[n]; }
    const pointer operator[](difference_type n) const { return (const pointer)it[n]; }
    
    friend iterator operator+(const iterator& i, difference_type n) { return iterator(i.it + n); }
    friend iterator operator+(difference_type n, const iterator& i) { return iterator(i.it + n); }
    friend iterator operator-(const iterator& i, difference_type n) { return iterator(i.it - n); }
    friend iterator operator-(difference_type n, const iterator& i) { return iterator(i.it - n); }
    
    friend difference_type operator-(const iterator& b, const iterator& a) { return (b.it - a.it); }
    
    friend bool operator<(const iterator&a, const iterator& b) { return (a.it < b.it); }
    friend bool operator>(const iterator&a, const iterator& b) { return (a.it > b.it); }
    friend bool operator<=(const iterator&a, const iterator& b) { return (a.it <= b.it); }
    friend bool operator>=(const iterator&a, const iterator& b) { return (a.it >= b.it); }
  };
  
public:
  MixedVector()
  {
  }
  
  MixedVector(const MixedVector<BaseT*>& other)
  : items(other.items), tmap(other.tmap)
  {
  }
  
  ~MixedVector()
  {
    clear();
  }
  
  friend void swap(MixedVector<BaseT*>& a, MixedVector<BaseT*>& b)
  {
    using std::swap;
    
    swap(a.items, b.items);
    swap(a.tmap, b.tmap);
  }
  
  MixedVector<BaseT*>& operator=(MixedVector<BaseT*> rhs)
  {
    // copy-and-swap idiom
    swap(*this, rhs);
    return *this;
  }
  
  BaseT& back()
  {
    return items.back();
  }
  
  const BaseT& back() const
  {
    return items.back();
  }
  
  void push_back(BaseT* obj)
  {
    // make sure this is a subclass of BaseT
    check_subclass(obj);
    
    // add the copy of the object pointer to the list of items
    items.push_back(obj);
    
    // add the copy to type map
    tmap[typeid(*obj)].push_back(obj);
  }
  
  void pop_back()
  {
    // get very last element
    BaseT* back = items.back();
    
    // remove from the vector in the type map
    tmap[typeid(*back)].pop_back();
    
    // remove from items
    items.pop_back();
    
    delete back; // delete object
  }
  
  template<typename SubclassT>
  void pop_back()
  {
    // search for vector belonging to SubclassT
    typename tmap_t::iterator found = tmap.find(typeid(SubclassT));
    if(found == tmap.end())
      return;
    
    // pop from vector for type
    found->second.pop_back();
    
    // go through list of all items from back, until item of correct type is found
    for(typename std::vector<BaseT*>::reverse_iterator it = items.rbegin(); it != items.rend(); ++it)
    {
      if(typeid(**it) == typeid(SubclassT))
      {
        // remove the item
        delete(*it);
        items.erase(it);
        return;
      }
    }
  }
  
  void clear()
  {
    while(!empty()) pop_back();
  }
  
  template<typename SubclassT>
  void clear()
  {
    // remove vector for subclass
    tmap.erase(typeid(SubclassT));
    
    // erase items matching type
    items.erase(std::remove_if(items.begin(), items.end(), is_type<SubclassT>), items.end());
  }
  
  BaseT** data()
  {
    return items.data();
  }
  
  const BaseT** data() const
  {
    return items.data();
  }
  
  template<typename SubclassT>
  bool type(std::size_t i)
  {
    return is_type<SubclassT>(items[i]);
  }
  
  BaseT* operator[](std::size_t i) const
  {
    return items[i];
  }
  
  BaseT* get(std::size_t i)
  {
    return items[i];
  }
  
  const BaseT* get(std::size_t i) const
  {
    return items[i];
  }
  
  template<typename SubclassT>
  SubclassT* get(std::size_t i)
  {
    typename tmap_t::iterator found = tmap.find(typeid(SubclassT));
    if(found == tmap.end())
      throw std::out_of_range(std::string() + "type " + typeid(SubclassT).name() + " not in vector");
    return (SubclassT*)found->second[i];
  }
  
  template<typename SubclassT>
  const SubclassT* get(std::size_t i) const
  {
    typename tmap_t::const_iterator found = tmap.find(typeid(SubclassT));
    if(found == tmap.end())
      throw std::out_of_range(std::string() + "type " + typeid(SubclassT).name() + " not in vector");
    return (const SubclassT*)found->second[i];
  }
  
  BaseT* at(std::size_t i)
  {
    return items.at(i);
  }
  
  const BaseT* at(std::size_t i) const
  {
    return items.at(i);
  }
  
  template<typename SubclassT>
  SubclassT* at(std::size_t i)
  {
    typename tmap_t::iterator found = tmap.find(typeid(SubclassT));
    if(found == tmap.end())
      throw std::out_of_range(std::string() + "type " + typeid(SubclassT).name() + " not in vector");
    return (SubclassT*)found->second.at(i);
  }
  
  template<typename SubclassT>
  const SubclassT* at(std::size_t i) const
  {
    typename tmap_t::const_iterator found = tmap.find(typeid(SubclassT));
    if(found == tmap.end())
      throw std::out_of_range(std::string() + "type " + typeid(SubclassT).name() + " not in vector");
    return (const SubclassT*)found->second.at(i);
  }
  
  std::size_t size() const
  {
    return items.size();
  }
  
  template<typename SubclassT>
  std::size_t size() const
  {
    typename tmap_t::const_iterator found = tmap.find(typeid(SubclassT));
    if(found == tmap.end())
      return 0;
    return found->second.size();
  }
  
  bool empty() const
  {
    return items.empty();
  }
  
  template<typename SubclassT>
  bool empty() const
  {
    typename tmap_t::const_iterator found = tmap.find(typeid(SubclassT));
    if(found == tmap.end())
      return true;
    
    return found->second.empty();
  }
  
  iterator<> begin()
  {
    return iterator<>(items.begin());
  }
  
  iterator<> end()
  {
    return iterator<>(items.end());
  }
  
  template<typename SubclassT>
  iterator<SubclassT> begin()
  {
    typename tmap_t::iterator found = tmap.find(typeid(SubclassT));
    if(found == tmap.end())
      return iterator<SubclassT>(items.end());
    return iterator<SubclassT>(found->second.begin());
  }
  
  template<typename SubclassT>
  iterator<SubclassT> end()
  {
    typename tmap_t::iterator found = tmap.find(typeid(SubclassT));
    if(found == tmap.end())
      return iterator<SubclassT>(items.end());
    return iterator<SubclassT>(found->second.end());
  }
  
private:
  inline void check_subclass(const BaseT*) {}
  
  template<typename SubclassT>
  inline static bool is_type(const BaseT* obj) { return (typeid(*obj) == typeid(SubclassT)); }
  
  typedef std::vector<BaseT*> base_container;
  
  base_container items;
  
  typedef std::map<detail::type_index, base_container> tmap_t;
  tmap_t tmap;
};
 
template<class BaseT>
std::size_t lower_bound(std::vector<BaseT*>& items, PosType target){
  std::size_t ju,jm,jl;
  
  if(items.size() == 0) return 0;
  
  jl=0;
  ju=items.size()-1;
  while (ju-jl > 1) {
    jm=(ju+jl) >> 1;
    if(items[jm]->compareZ(target))
      jl=jm;
    else
      ju=jm;
  }
  return jl;
}
 
template<typename Container>
void delete_container(Container& c) { while(!c.empty()) delete c.back(), c.pop_back(); }
 
class HilbertCurve{
public:
  HilbertCurve(
               PosType x_min       
               ,PosType y_min      
               ,PosType my_range   
               ,PosType smallsize  
  ):
  range(my_range)
  {
    xo[0] = x_min;
    xo[1] = y_min;
    n = (int)(range/smallsize+1);
  }
  
  int xy2d (int x, int y);
  int xy2d (PosType x, PosType y);
  void d2xy(int d, int *x, int *y);
  void d2xy(int d, PosType *x, PosType *y);
  
private:
  int n;
  PosType xo[2],range;
  void rot(int s,int *x, int *y, int rx, int ry);
};
 
template<typename T>
T between(const T& x, const T& l, const T& u)
{
  return std::max(l, std::min(u, x));
}
 
#ifdef ENABLE_CLANG
class RandomNumbers{
public:
  
  RandomNumbers(unsigned int seed);
  ~RandomNumbers(void);
  
  PosType operator()(void);
  PosType gauss(){return norm_dist(rand_gen);}
private:
  
  std::normal_distribution<> norm_dist;
  std::mt19937 rand_gen;
};
#endif
 
class RandomNumbers_NR{
public:
  
  RandomNumbers_NR(long seed);
  // sets seed from clock
  RandomNumbers_NR();
  
  PosType operator()(void);
  PosType gauss(){
    ++calls;
    if(count){
      do{
        u = 2*ran2() - 1;
        v = 2*ran2() - 1;
        s = u*u +v*v;
      }while( s > 1.0 || s == 0.0);
      
      s = sqrt(-2*log(s)/s);
      count = false;
      return s*u;
    }else{
      count = true;
      return s*v;
    }
  };
  
  int poisson(double lam){
    if(lam > 200){
      return std::max<long>(std::lround( gauss()*sqrt(lam) + lam ) ,0);
    }
    double L = exp(-lam),p=1;
    int k = 0;
    do{
      ++k;
      p *= operator()();
    }while(p>L);
    
    return k-1;
  }
  
//  int poisson(double lam){
//
//    int step = 500;
//    double Lam = lam,p=1;
//    int k = 0;
//    do{
//      ++k;
//      p *= operator()();
//      while(p<1 && Lam>0){
//        if(Lam > step){
//          p *= exp(step);
//          Lam -= step;
//        }else{
//          p *= exp(Lam);
//          Lam=0;
//        }
//      }
//    }while(p>1);
//
//    return k-1;
//  }
  
  size_t calls = 0;  
  
  long getseed(){return firstseed;}
private:
  long idum;
  PosType ran2(void);
  
  int IM1;
  int IM2;
  PosType AM;
  //int IMM1 = (IM1-1);
  int IA1;
  int IA2;
  int IQ1;
  int IQ2;
  int IR1;
  int IR2;
  int NDIV;
  PosType EPS;
  PosType RNMX;
  
  long idum2;
  long iy;
  long iv[32];
  bool count;
  PosType u,v,s;
  long firstseed;
  
};
 
class BetaDistribution {
public:
    using result_type = double;
 
    BetaDistribution(double alpha, double beta)
        : alpha_(alpha),
          beta_(beta),
          gamma_alpha_(alpha, 1.0),
          gamma_beta_(beta, 1.0)
    {
        if (!(alpha > 0.0) || !(beta > 0.0)) {
            throw std::invalid_argument(
                "Beta distribution parameters must be positive");
        }
    }
 
    double operator()(){
        // The loop protects against the extremely unlikely case that both
        // Gamma samples underflow to zero.
        for (;;) {
            const double x = gamma_alpha_(engine);
            const double y = gamma_beta_(engine);
            const double sum = x + y;
 
            if (sum > 0.0 && std::isfinite(sum)) {
                return x / sum;
            }
        }
    }
 
    double alpha() const noexcept { return alpha_; }
    double beta() const noexcept { return beta_; }
 
    static constexpr double min() noexcept { return 0.0; }
    static constexpr double max() noexcept { return 1.0; }
 
private:
    
    std::mt19937_64 engine{std::random_device{}()};
    double alpha_;
    double beta_;
    std::gamma_distribution<double> gamma_alpha_;
    std::gamma_distribution<double> gamma_beta_;
};
 
 
 
template <typename T, typename R>
void shuffle(
             std::vector<T> &vec   
             ,R &ran               
){
 
  size_t ran_t;
  if(vec.size() < 2) return;
  for (size_t i = vec.size()-1; i>0; --i) {
    ran_t = (size_t)(ran()*(i+1));
    //swap(tmp,vec[ran_t]);
    std::swap(vec[ran_t],vec[i]);
    //swap(vec[i],tmp);
  }
}
 
template <typename T>
void sort_indexes(const std::vector<T> &v     
                  ,std::vector<size_t> &index 
) {
  
  // initialise original index locations
  index.resize(v.size());
  for (size_t i = 0; i != index.size(); ++i) index[i] = i;
  
  // sort indexes based on comparing values in v
  std::sort(index.begin(), index.end(),[&v](size_t i1, size_t i2) {return v[i1] < v[i2];});
}
 
template <typename T>
void sort_indexes(const T *v     
                  ,std::vector<size_t> &index 
                  ,size_t N) {
  
  // initialise original index locations
  index.resize(N);
  
  for (size_t i = 0; i != index.size(); ++i) index[i] = i;
  
  // sort indexes based on comparing values in v
  std::sort(index.begin(), index.end(),
            [v](size_t i1, size_t i2) {return v[i1] < v[i2];});
}
 
template <typename T>
void sort_indexes_decending(const std::vector<T> &v     
                            ,std::vector<size_t> &index 
) {
  
  // initialise original index locations
  index.resize(v.size());
  for (size_t i = 0; i != index.size(); ++i) index[i] = i;
  
  // sort indexes based on comparing values in v
  std::sort(index.begin(), index.end(),
            [&v](size_t i1, size_t i2) {return v[i1] > v[i2];});
}
 
template <typename R>
class ShuffledIndex{
public:
  ShuffledIndex(size_t N,R &ran){
    if(N == 0) throw std::invalid_argument("N=0");
    
    for(size_t i=0 ; i<N ; ++i) index[i] = i;
    index.resize(N);
    Utilities::shuffle(index,ran);
    index_internal = 0;
  }
  
  size_t operator*(){return index[index_internal];}
  size_t operator++(int){
    size_t tmp = index[index_internal];
    if(index_internal == index.size()){
      index_internal = 0;
    }else{
      ++index_internal;
    }
    return tmp;
  }
  size_t operator++(){
    if(index_internal == index.size()){
      index_internal = 0;
    }else{
      ++index_internal;
    }
    return index[index_internal];
  }
  
  void reshuffle(R &ran){
    Utilities::shuffle(index,ran);
    index_internal = 0;
  }
  
private:
  std::vector<size_t> index;
  size_t index_internal;
};
 
 
// reorders vec according to index p
template <typename T>
void apply_permutation(
                       T *vec,
                       const std::vector<std::size_t>& p)
{
  std::vector<T> copy(p.size());
  
  for(std::size_t i = 0; i < p.size(); ++i){
    copy[i] = vec[i];
  }
  
  for(std::size_t i = 0; i < p.size(); ++i){
    vec[i] = copy[p[i]];
  }
}
 
template <typename T>
void apply_permutation(
                       std::vector<T>& vec,
                       const std::vector<std::size_t>& p)
{
  apply_permutation(vec.data(),p);
}
 
 
void powerspectrum2d(
                     std::valarray<double> const &aa      
                     ,std::valarray<double> const &bb     
                     ,long nx                       
                     ,long ny                       
                     ,double boxlx                 
                     ,double boxly                 
                     ,std::vector<double> &ll      
                     ,std::vector<double> &Pl      
                     ,double zeropaddingfactor
                     );
 
void powerspectrum2d(
                     std::valarray<double> &aa      
                     ,long nx                       
                     ,long ny                       
                     ,double boxlx                 
                     ,double boxly                 
                     ,std::vector<double> &ll      
                     ,std::vector<double> &Pl      
);
 
void powerspectrum2dprebin(
                     std::valarray<double> &aa      
                     ,long nx                       
                     ,long ny                       
                     ,double boxlx                 
                     ,double boxly                 
                     ,const std::vector<double> &ll      
                     ,std::vector<double> &Pl      
                     ,std::vector<double> &llave     
);
 
void powerspectrum2d(
                     std::valarray<float> const &aa      
                     ,std::valarray<float> const &bb     
                     ,long nx                       
                     ,long ny                       
                     ,double boxlx                 
                     ,double boxly                 
                     ,std::vector<double> &ll      
                     ,std::vector<double> &Pl      
                     ,double zeropaddingfactor
                     );
 
void powerspectrum2d(
                     std::valarray<float> &aa      
                     ,long nx                       
                     ,long ny                       
                     ,double boxlx                 
                     ,double boxly                 
                     ,std::vector<double> &ll      
                     ,std::vector<double> &Pl      
);
 
void powerspectrum2dprebin(
                     std::valarray<float> &aa      
                     ,int nx                       
                     ,int ny                       
                     ,double boxlx                 
                     ,double boxly                 
                     ,const std::vector<double> &ll      
                     ,std::vector<double> &Pl      
                     ,std::vector<double> &llave     
);
 
int GetNThreads();
 
namespace IO{  
 
inline bool file_exists (const std::string& name) {
  struct stat buffer;
  return (stat (name.c_str(), &buffer) == 0);
}
  
  template <class T>
  void read1columnfile(
                       std::string filename    
                       ,std::vector<T> &x     
                       ,int skiplines = 0
                       ,bool verbose = false
                       ){
 
    x.clear();
 
    std::ifstream file_in(filename.c_str());
    std::string myline;
    std::string space = " ";
    T myt1;
 
    std::string strg;
    std::stringstream buffer;
 
    if(!file_in){
      std::cout << "Can't open file " << filename << std::endl;
      ERROR_MESSAGE();
      throw std::runtime_error(" Cannot open file.");
    }
 
    std::cout << "Reading from " << filename << std::endl;
    size_t i=0;
    while(i < skiplines){
      getline(file_in,myline);
      ++i;
    }
    while(file_in.peek() == '#'){
      file_in.ignore(10000,'\n');
      ++i;
    }
    std::cout << "skipped "<< i << " comment lines in " << filename << std::endl;
 
    size_t pos;
    // read in data
    while(getline(file_in,myline)){
 
      if(myline[0] == '#'){
        std::cout << "skipped line " << i << std::endl;
        continue;
      }
 
      pos= myline.find_first_not_of(space);
      myline.erase(0,pos);
 
      buffer << myline;
      buffer >> myt1;
      if(verbose) std::cout << myt1 << " ";
      x.push_back(myt1);
 
      strg.clear();
      buffer.clear();
      myline.clear();
 
    }
    std::cout << "Read " << x.size() << " lines from " << filename << std::endl;
  }
 
template <class T1,class T2>
void read2columnfile(
                     std::string filename    
                     ,std::vector<T1> &x     
                     ,std::vector<T2> &y     
                     ,std::string delineator = " "  
                     ,int skiplines = 0
                     ,bool verbose = false
                     
                     ){
  
  x.clear();
  y.clear();
  
  std::ifstream file_in(filename.c_str());
  std::string myline;
  std::string space = " ";
  T1 myt1;
  T2 myt2;
  
  std::string strg;
  std::stringstream buffer;
  
  if(!file_in){
    std::cout << "Can't open file " << filename << std::endl;
    ERROR_MESSAGE();
    throw std::runtime_error(" Cannot open file.");
  }
  
  std::cout << "Reading caustic information from " << filename << std::endl;
  size_t i=0;
  while(i < skiplines){
    getline(file_in,myline);
    ++i;
  }
  while(file_in.peek() == '#'){
    file_in.ignore(10000,'\n');
    ++i;
  }
  std::cout << "skipped "<< i << " comment lines in " << filename << std::endl;
  
  size_t pos;
  // read in data
  while(getline(file_in,myline)){
    
    if(myline[0] == '#'){
      std::cout << "skipped line " << i << std::endl;
      continue;
    }
    
    pos= myline.find_first_not_of(space);
    myline.erase(0,pos);
    
    
    pos = myline.find(delineator);
    strg.assign(myline,0,pos);
    buffer << strg;
    buffer >> myt1;
    if(verbose) std::cout << myt1 << " ";
    x.push_back(myt1);
    
    myline.erase(0,pos+1);
    pos= myline.find_first_not_of(space);
    myline.erase(0,pos);
    
    strg.clear();
    buffer.clear();
    
    pos = myline.find(space);
    strg.assign(myline,0,pos);
    buffer << strg;
    buffer >> myt2;
    if(verbose)  std::cout << myt2 << std::endl;
    y.push_back(myt2);
    
    strg.clear();
    buffer.clear();
    myline.clear();
    
  }
  std::cout << "Read " << x.size() << " lines from " << filename << std::endl;
}
template <class T1,class T2,class T3>
void read3columnfile(
                     std::string filename    
                     ,std::vector<T1> &x     
                     ,std::vector<T2> &y     
                     ,std::vector<T3> &z     
                     ,std::string delineator = " "  
                     ,bool verbose = false
                     
                     ){
  
  
  assert(0); // Untested!!!!
  x.clear();
  y.clear();
  z.clear();
  
  std::ifstream file_in(filename.c_str());
  std::string myline;
  std::string space = " ";
  T1 myt1;
  T2 myt2;
  T3 myt3;
  
  std::string strg;
  std::stringstream buffer;
  
  if(!file_in){
    std::cout << "Can't open file " << filename << std::endl;
    ERROR_MESSAGE();
    throw std::runtime_error(" Cannot open file.");
  }
  
  std::cout << "Reading caustic information from " << filename << std::endl;
  size_t i=0;
  while(file_in.peek() == '#'){
    file_in.ignore(10000,'\n');
    ++i;
  }
  std::cout << "skipped "<< i << " comment lines in " << filename << std::endl;
  
  size_t pos;
  // read in data
  while(getline(file_in,myline)){
    
    if(myline[0] == '#'){
      std::cout << "skipped line " << i << std::endl;
      continue;
    }
    
    pos= myline.find_first_not_of(space);
    myline.erase(0,pos);
    
    
    pos = myline.find(delineator);
    strg.assign(myline,0,pos);
    buffer << strg;
    buffer >> myt1;
    if(verbose) std::cout << myt1 << " ";
    x.push_back(myt1);
    
    myline.erase(0,pos+1);
    pos= myline.find_first_not_of(space);
    myline.erase(0,pos);
    
    strg.clear();
    buffer.clear();
    
    // ******************
    
    
    pos = myline.find(space);
    strg.assign(myline,0,pos);
    buffer << strg;
    buffer >> myt2;
    if(verbose)  std::cout << myt2 << std::endl;
    y.push_back(myt2);
    
    myline.erase(0,pos+1);
    pos= myline.find_first_not_of(space);
    myline.erase(0,pos);
    
    strg.clear();
    buffer.clear();
    
    // ******************
    pos = myline.find(space);
    strg.assign(myline,0,pos);
    buffer << strg;
    buffer >> myt3;
    if(verbose)  std::cout << myt3 << std::endl;
    y.push_back(myt3);
    
    strg.clear();
    buffer.clear();
    myline.clear();
    
  }
  std::cout << "Read " << x.size() << " lines from " << filename << std::endl;
}
 
int NumberOfEntries(const std::string &string,char deliniator);
 
int CountColumns(std::string filename  
                 ,char comment_char = '#'  
                 ,char deliniator = ' '    
);
 
void ReadFileNames(
                   std::string dir              
                   ,const std::string filespec 
                   ,std::vector<std::string> & filenames  
                   ,bool verbose);
 
bool check_directory(std::string dir);
bool make_directories(const std::string &root_dir);
 
 
template <typename T>
void ReadASCII(std::vector<T> &data
               ,std::string filename
               ,int &columns
               ,int &rows
               ,char comment_char = '#'
               ,int skiplines = 0
               ,size_t MaxNrows = std::numeric_limits<size_t>::max()
               ,bool verbose = true){
  
  std::ifstream file(filename);
  // find number of particles
  if (!file.is_open()){
    std::cerr << "file " << filename << " cann't be opened." << std::endl;
    throw std::runtime_error("no file");
  }
  
  
  data.empty();
  
  std::string line;
  columns = 0;
  rows = 0;
  size_t first_data_line;
  
  // skip over first lines
  int i=0;
  while(i < skiplines){
    std::getline(file,line);
    if(!file) break;  // probably EOF
    ++i;
  }
  
  // read comment lines and first data line
  do{
    first_data_line = file.tellg();
    std::getline(file,line);
    if(!file) break;  // probably EOF
  }while(line[0] == comment_char);
  
  columns =  NumberOfEntries(line,' ');
  
  file.seekg(first_data_line);   // move back to first data line
  
  std::copy(std::istream_iterator<T>(file),
            std::istream_iterator<T>(),
            std::back_inserter(data));
  
  rows = data.size()/columns;
  if(verbose){
    std::cout << "Read " << rows << " rows of " << columns << " columns from file " << filename << std::endl;
  }
}
 
template <typename T>
int ReadCSVnumerical1(std::string filename                              
                      ,std::vector<std::vector<T> > &data               
                      ,std::vector<std::string> &column_names           
                       ,size_t MaxNumber = 100000000                     
                      ,char comment_char = '#'                          
                      ,char deliniator = ','                            
                      ,std::string replace = "\\N"                      
                      ,std::function<bool(std::vector<T> &)> accept = [](std::vector<T> &v){return true;}  
  ){
    
    bool verbose = false;
    
    std::ifstream file(filename);
    // find number of particles
    if (!file.is_open()){
      std::cerr << "file " << filename << " cann't be opened." << std::endl;
      throw std::runtime_error("no file");
    }
    
    int count_preamble=0;
    std::string line;
    // read comment lines and first data line
    do{
      std::getline(file,line);
      ++count_preamble;
      if(!file) break;  // probably EOF
    }while(line[0] == comment_char);
    
    // read the names
    std::stringstream          lineStream(line);
    std::string                cell;
    column_names.empty();
    while(std::getline(lineStream,cell, deliniator))
    {
      column_names.push_back(cell);
    }
    // This checks for a trailing comma with no data after it.
    if (!lineStream && cell.empty())
    {
      column_names.push_back("");
    }
    
    if(verbose){ // print column names
      int i = 0;
      for(auto st : column_names){
        std::cout << i++ << " " << st << std::endl;
      }
    }
    
    int ncolumns = NumberOfEntries(line,deliniator);
    
    // check if any of the columns are strings
    std::getline(file,line);
    {
      size_t n=0,m = 0,i=0;
      do{
        m = line.find(",",n);
        cell = line.substr(n,m-n);
        cell.erase(remove_if(cell.begin(),cell.end(), isspace), cell.end());
        n = m+1;
        ++i;
      }while(m != std::string::npos);
    }
  
  // count number of data line
  size_t number_of_data_lines = 1;
  while(std::getline(file, line) && number_of_data_lines < MaxNumber ) ++number_of_data_lines;
  
  //std::vector<std::vector<T> > tmp_data(columns,std::vector<T>(number_of_data_lines));
  std::vector<std::vector<T> > tmp_data(ncolumns,std::vector<T>(0));
  swap(data,tmp_data);
  std::vector<T> tmp_row(ncolumns);
  
  //file.seekg(std::ios::beg);
  file.clear();
  file.close();
  file.open(filename);
  for(int i=0; i < count_preamble; ++i){
    file.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
  }
  
  size_t rows = 0;
  while(std::getline(file,line) && rows < MaxNumber){
    // read the names
    std::stringstream          lineStream(line);
    std::string                cell;
    
    int i=0;
    while(std::getline(lineStream,cell,deliniator))
    {
      if(cell==replace) cell='0';
      
      cell.erase(remove_if(cell.begin(),cell.end(), isspace), cell.end());
      //data[i].push_back(to_numeric<T>(cell));
      //data[i][rows] = to_numeric<T>(cell);
      tmp_row[i] = to_numeric<T>(cell);
      i = (i+1)%ncolumns;
    }
    if(accept(tmp_row)){
      ++rows;
      for(int i=0 ; i < ncolumns ; ++i) data[i].push_back(tmp_row[i]);
    }
  }
  return 1;
}
 
template <typename T>
size_t ReadCSVrange(std::string filename                              
                      ,std::vector<std::vector<T> > &ranges               
                      ,std::vector<std::string> &column_names           
                       ,size_t MaxNumber = 100000000                     
                      ,char comment_char = '#'                          
                      ,char deliniator = ','                            
                      ,std::string replace = "\\N"                      
                      ,std::function<bool(std::vector<T> &)> accept = [](std::vector<T> &v){return true;}  
){
  
  bool verbose = false;
  
  std::ifstream file(filename);
  // find number of particles
  if (!file.is_open()){
    std::cerr << "file " << filename << " cann't be opened." << std::endl;
    throw std::runtime_error("no file");
  }
  
  int count_preamble=0;
  std::string line;
  // read comment lines and first data line
  do{
    std::getline(file,line);
    ++count_preamble;
    if(!file) break;  // probably EOF
  }while(line[0] == comment_char);
  
  // read the names
  std::stringstream          lineStream(line);
  std::string                cell;
  column_names.empty();
  while(std::getline(lineStream,cell, ','))
  {
    column_names.push_back(cell);
  }
  // This checks for a trailing comma with no data after it.
  if (!lineStream && cell.empty())
  {
    column_names.push_back("");
  }
  
  if(verbose){ // print colum names
    int i = 0;
    for(auto st : column_names){
      std::cout << i++ << " " << st << std::endl;
    }
  }
  
  int columns = NumberOfEntries(line,deliniator);
  
  std::vector<std::vector<T> > tmp_ranges(columns,std::vector<T>(2));
  swap(ranges,tmp_ranges);
  std::vector<T> tmp_row(columns);
  
  //file.seekg(std::ios::beg);
  file.clear();
  file.close();
  file.open(filename);
  for(int i=0; i < count_preamble; ++i){
    file.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
  }
  //std::getline(file,line);
  
  size_t rows = 0;
  while(std::getline(file,line) && rows < MaxNumber){
    // read the names
    std::stringstream          lineStream(line);
    std::string                cell;
    
    int i=0;
    while(std::getline(lineStream,cell,deliniator))
    {
      if(cell==replace) cell='0';
      
      cell.erase(remove_if(cell.begin(),cell.end(), isspace), cell.end());
      tmp_row[i] = to_numeric<T>(cell);
      i = (i+1)%columns;
    }
    if(accept(tmp_row)){
      if(rows==0){
        for(int j=0 ; j < columns ; ++j){
          ranges[j][0] = tmp_row[j];
          ranges[j][1] = tmp_row[j];
 
        }
      }else{
        for(int j=0 ; j < columns ; ++j){
          if(ranges[j][0] > tmp_row[j]) ranges[j][0] = tmp_row[j];
          if(ranges[j][1] < tmp_row[j]) ranges[j][1] = tmp_row[j];
        }
      }
      ++rows;
    }
  }
  return rows;
}
 
template <typename T>
int ReadCSVnumerical2(std::string filename   
                      ,std::vector<std::vector<T> > &data  
                      ,std::vector<std::string> &column_names 
                      ,size_t Nmax = 1000000
                      ,char comment_char = '#'  
                      ,char deliniator = ','    
                      ,bool header = true    
                      ,std::string reject = ""  
){
  
  
  std::ifstream file(filename);
  
  if (!file.is_open()){
    std::cerr << "file " << filename << " cann't be opened." << std::endl;
    throw std::runtime_error("no file");
  }
  std::string line;
  // read comment lines and first data line
  do{
    std::getline(file,line);
    if(!file) break;  // probably EOF
  }while(line[0] == comment_char);
  
  // read the names
  std::stringstream          lineStream(line);
  std::string                cell;
  column_names.empty();
  long ii = 0;
  
  if(header){
    while(std::getline(lineStream,cell,deliniator))
    {
      column_names.push_back(cell);
    }
    
    // This checks for a trailing comma with no data after it.
    if (!lineStream && cell.empty())
    {
      column_names.push_back("");
    }
  }else{
    
    while(std::getline(lineStream,cell,deliniator))
    {
      column_names.push_back(cell);
    }
    
    data.emplace_back(column_names.size());
    int i=0;
    for(auto a : column_names){
      cell.erase(remove_if(cell.begin(),cell.end(), isspace), cell.end());
      data.back()[i++] = to_numeric<T>(a);
    }
    ++ii;
  }
  
  int columns = NumberOfEntries(line,deliniator);
  
  for(auto &v : data ) v.empty();
  
  while(std::getline(file,line) && ii < Nmax){
    
    data.emplace_back(columns);
    
    // read the numbers
    std::stringstream          lineStream(line);
    std::string                cell;
    
    int i=0;
    while(std::getline(lineStream,cell, deliniator))
    {
      assert(i < columns);
      cell.erase(remove_if(cell.begin(),cell.end(), isspace), cell.end());
      if(cell == reject){
        data.pop_back();
        --ii;
        break;
      }else{
        data.back()[i] = to_numeric<T>(cell);
        ++i;
      }
    }
    ++ii;
  }
  
  return 1;
}
 
template<typename T>
void writeCSV(const std::string filename              
              ,const std::vector<std::string> header  
              ,std::vector<T *> &data                 
){
  
  std::ofstream s(filename + ".csv");
  
  long ncol = header.size();
  assert(ncol == data.size() );
  for(int i = 0 ; i < header.size()-1 ; ++i){
    s << header[i] << ",";
  }
  s << header.back() << std::endl;
  
  size_t nrow = data[0]->size();
  for(auto v : data) assert(nrow == v->size() );
  
  for(size_t j = 0 ; j < nrow ; ++j){
    for(int i = 0 ; i < ncol-1 ; ++i){
      s << data[i]->operator[](j) << ",";
    }
    s << data.back()->operator[](j) << "\n";
  }
}
 
/*** \brief A class for reading and then looking up objects from a CSV catalog.
 
 The constructor will read in the whole catalog and sort it into Nxbins X-bins.  Each
 X-bin is then sorted by Y.
 
 The find(x,y) function will set the current value to a galaxy with a x within the
 x-bin of x and a y close to y.  The other information in the row of the csv file for
 this entry can then be read.
 
 */
class XYcsvLookUp{
public:
  XYcsvLookUp(std::string datafile   
              ,std::string Xlabel    
              ,std::string Ylabel    
              ,int Nxbins            
              ,size_t MaxNumber = 1000000 
              ,bool verbose = false);
  XYcsvLookUp(
              std::string datafile   
              ,std::string Xlabel    
              ,std::string Ylabel    
              ,std::vector<double> Xbins  
              ,size_t MaxNumber = 1000000 
              ,bool verbose = false);
  
  std::vector<double> find(double x,double y);
  
  double operator[](std::string label) const;
  
  double operator[](int i) const{
    return (*current)[i];
  }
  double getX() const {
    return (*current)[Xindex];
  }
  double getY() const {
    return (*current)[Yindex];
  }
  std::vector<double> record() const{
    return *current;
  }
  std::vector<double> operator*() const {
    return *current;
  }
  
  void operator++(){
    if(current != data.end() ) ++current;
  }
  void operator--(){
    if(current != data.begin() ) --current;
  }
  std::vector<std::string> labels() const{
    return column_names;
  }
  
  double Xmin() const {return xmin;}
  double Xmax() const {return xmax;}
  double Ymin(double x) const;
  double Ymax(double x) const;
  
  void printcurrent(){
    int i = 0;
    for(auto label : column_names){
      std::cout << label << " : " << (*current)[i++] << std::endl;
    }
  }
  
private:
  int Xindex;
  int Yindex;
  double xmax;
  double xmin;
  std::vector<std::vector<double> >::iterator current;
  std::vector<std::vector<std::vector<double> >::iterator> borders;
  std::vector<std::vector<double> > data;
  std::vector<std::string> column_names;
  std::vector<double> Xborders;
  size_t NinXbins;
  std::string filename;
};
} // Utilities::IO
 
void splitstring(std::string &line,std::vector<std::string> &vec,const std::string &delimiter);
 
template< typename T>
class DataFrame{
public:
  DataFrame(std::string datafile   
            ,size_t MaxNumber = 1000000 
            ,char comment_char = '#'  
            ,char deliniator = ','    
            ,std::string replace = "\\N"    
            ,std::function<bool(std::vector<T> &)> accept = [](std::vector<T> &v){return true;}  
  ):filename(datafile){
    input(datafile,MaxNumber,comment_char,deliniator,replace,accept);
  };
 
  DataFrame(){};
  
  void input(std::string datafile   
            ,size_t MaxNumber = 1000000 
            ,char comment_char = '#'  
            ,char deliniator = ','    
            ,std::string replace = "\\N"    
            ,std::function<bool(std::vector<T> &)> accept = [](std::vector<T> &v){return true;}  
  ){
    filename = datafile;
    Utilities::IO::ReadCSVnumerical1(datafile,data,column_names,MaxNumber,'#',',',replace,accept);
    for(int i=0 ; i<column_names.size() ; ++i) datamap[column_names[i]] = i;
  };
 
  void pop(std::string colname){
    
    int i=datamap[colname];
 
    swap(column_names[i],column_names.back());
    column_names.pop_back();
    swap(data[i],data.back());
    data.pop_back();
    
    datamap.empty();
    for(int i=0 ; i<column_names.size() ; ++i) datamap[column_names[i]] = i;
  }
  
  std::vector<T>& operator[](const std::string &label){
    if(datamap.find(label) == datamap.end()){
      std::cerr << "No label - " << label << " - in " << filename <<std::endl;
      throw std::invalid_argument("no label");
    }
    return data[datamap[label]];
    //for(auto c : column_names ) std::cout << c << " ";
    //std::cout << std::endl;
    //throw std::invalid_argument(label + " was not one of the columns of the galaxy data file :" + filename);
  };
  const std::vector<T>& operator[](const std::string &label) const{
    if(datamap.find(label) == datamap.end()){
      std::cerr << "No label - " << label << " - in " << filename <<std::endl;
      throw std::invalid_argument("no label");
    }
    return data[datamap[label]];
  };
  
  void add_column(const std::string &name,const std::vector<T> &vec){
    if(vec.size() != data[0].size()){
      std::cerr << "Added column to DataFrame needs to have the same size." << std::endl;
      throw std::invalid_argument("wrong size");
    }
    column_names.push_back(name);
    data.push_back(vec);
    data[name] = data.size()-1;
  };
  
  std::vector<T>& operator[](int i){
    return data[i];
  };
  const std::vector<T>& operator[](int i) const{
    return data[i];
  };
  
  std::vector<std::string> labels() const{
    return column_names;
  };
  
  // sort by one of the columns
  void sortby(std::string name){
    std::vector<size_t> index(data[0].size());
    size_t N = index.size();
    for(size_t i=0 ; i<N ; ++i) index[i] = i;
    
    sort_indexes(data[datamap[name]],index);
    std::vector<T> tmp_v(N);
    
    for(size_t j=0 ; j<data.size() ; ++j){
      for(size_t i=0 ; i<N ; ++i){
        tmp_v[i] = data[j][index[i]];
      }
      swap(data[j],tmp_v);
    }
  }
  
  // shuffle order of rows
  void shuffle(Utilities::RandomNumbers_NR &ran){
    std::vector<size_t> index(data[0].size());
    size_t N = index.size();
    for(size_t i=0 ; i<N ; ++i) index[i] = i;
    Utilities::shuffle(index, ran);
    
    std::vector<T> tmp_v(N);
    
    for(size_t j=0 ; j<data.size() ; ++j){
      for(size_t i=0 ; i<N ; ++i){
        tmp_v[i] = data[j][index[i]];
      }
      swap(data[j],tmp_v);
    }
  }
 
  size_t number_of_rows(){return data[0].size();}
  size_t number_of_columns(){return data.size();}
  
  std::vector<std::vector<T> > data;
private:
  std::map<std::string,int> datamap;
  std::vector<std::string> column_names;
  std::string filename;
};
 
// this is a summation algorithm that maintains percision for large sequences of numbers
template <typename T>
double KleinSum(std::vector<T> &input){
  double s = 0.0,cs = 0.0,ccs = 0.0,c,cc;
  size_t n = input.size();
  for(size_t i = 0 ; i<n ; ++i){
    double t = s + input[i];
    if( abs(s) >= abs(input[i]) ){
      c = (s - t) + input[i];
    }else{
      c = (input[i] - t) + s;
    }
    s = t;
    t = cs + c;
    if( abs(cs) >= abs(c) ){
      cc = (cs - t) + c;
    }else{
      cc = (c - t) + cs;
    }
    cs = t;
    ccs = ccs + cc;
  }
  return s + cs + ccs;
}
 
 
template <typename T,typename F>
double PairWiseSum(T *begin,T *end,F value){
  double sum = 0;
  size_t n = end - begin;
  if(n <= 10){
    for(T* p=begin ; p != end ; ++p){
      sum += value(p);
    }
  }else{
    sum = PairWiseSum<T,F>(begin,begin + n/2,value) + PairWiseSum<T,F>(begin + n/2,end,value);
  }
  return sum;
}
 
template <typename T>
double PairWiseSum(T *begin,T *end){
  double sum = 0;
  size_t n = end - begin;
  if(n <= 10){
    for(T* p=begin ; p != end ; ++p){
      sum += *p;
    }
  }else{
    sum = PairWiseSum(begin,begin + n/2) + PairWiseSum(begin + n/2,end);
  }
  return sum;
}
 
template <typename T>
double PairWiseSum(std::vector<T> &v){
  return PairWiseSum(v.data(),v.data() + v.size());
}
 
template <typename T,typename F>
double PairWiseSum(std::vector<T> &v,F value){
  return PairWiseSum(v.data(),v.data() + v.size(),value);
}
 
template <typename T>
class SUMMER{
 
public:
 
  SUMMER():batch(1000),n(0),ntotal(0),current(0.0){};
  SUMMER(size_t batchsize)
  :batch(batchsize),n(0),current(0.0),ntotal(0){};
 
  void operator+=(T x){
    ++n;
    current += x;
    ++ntotal;
    if(n % batch == 0){
      v.push_back(current);
      n=0;
      current = 0;
    }
  }
  
  T operator*(){
    if(v.size() ==0 ) return 0;
    v.push_back(current);
    return Utilities::PairWiseSum(v);
  }
  
  void reset(){
    std::vector<T> dump;
    std::swap(v,dump);
    n = ntotal = 0;
    current = 0;
  }
  
  size_t total_number(){
    return ntotal;
  }
  
private:
  size_t batch;
  size_t n;
  size_t ntotal;
  T current;
  std::vector<T> v;
};
 
/*
This is 1F2(a,b,c:x)
 
 Note convergence restrictions: abs(x) < 1 and c not a negative integer or zero
*/
template <typename T = double>
double hypergeometric( T a, T b, T c, T x )
{
 
  if(abs(x) > 1 ){
    throw std::invalid_argument("out of bounds");
  }
  const double TOLERANCE = 1.0e-10;
   double term = a * b * x / c;
   double value = 1.0 + term;
   int n = 1;
 
   while ( abs( term ) > TOLERANCE )
   {
      a++, b++, c++, n++;
      term *= a * b * x / c / n;
      value += term;
   }
 
   return value;
}
 
class LOGPARAMS{
 
public:
  typedef boost::variant<int,size_t,long,float,double,std::string,bool> MULTITYPE;
  typedef std::map<std::string,MULTITYPE> LINE;
  typedef std::map<std::string,std::vector<double> > VLINE;
  typedef std::map<std::string,std::vector<bool> > BLINE;
 
private:
  time_t to;
  
  LINE lines;
  VLINE vlines;
  BLINE blines;
  std::map<std::string,std::string> input_lines;
  std::string filename;
  std::string blanck_val;
public:
 
  LOGPARAMS(std::string file,std::string blanck_value = "0"):filename(file),blanck_val(blanck_value){
    time(&to);
  };
  
  ~LOGPARAMS(){
    print_file(filename);
  }
  
  void operator()(std::string label,MULTITYPE v){lines[label] = v;}
  void operator()(std::string label,const std::vector<double> &v){vlines[label] = v;}
  void operator()(std::string label,const std::vector<bool> &v){blines[label] = v;}
  void operator()(std::string label){lines[label] = "-----";}
 
  
  MULTITYPE operator[](std::string label){return lines[label];}
  void setlogfile(std::string name){filename = name;}
 
  void setinput(const std::string &input_filename){
    std::ifstream input_file(input_filename);
    if(!input_file){
      throw std::runtime_error("Unable to open parameter file: " + input_filename);
    }
 
    input_lines.clear();
    std::string line;
    const std::string whitespace = " \t\r\n";
    while(std::getline(input_file,line)){
      const std::string::size_type value_end = line.find_last_not_of(whitespace);
      if(value_end == std::string::npos) continue;
 
      const std::string::size_type value_begin = line.find_last_of(whitespace,value_end);
      if(value_begin == std::string::npos) continue;
 
      const std::string::size_type label_end = line.find_last_not_of(whitespace,value_begin);
      if(label_end == std::string::npos) continue;
 
      const std::string label = line.substr(0,label_end+1);
      const std::string value = line.substr(value_begin+1,value_end-value_begin);
      input_lines[label] = value;
    }
  }
 
  bool get(const std::string &label) {
    const std::map<std::string,std::string>::const_iterator entry = input_lines.find(label);
    if(entry == input_lines.end()) return false;
 
    // The text log does not retain the original variant type, so preserve the
    // value exactly as it appeared in the input file.
    lines[label] = entry->second;
    return true;
  }
 
  template<typename T>
  bool get(const std::string &label,T &value) {
    const std::map<std::string,std::string>::const_iterator entry = input_lines.find(label);
    if(entry == input_lines.end()){
      lines[label] = value;
      return false;
    }
 
    std::istringstream stream(entry->second);
    T parsed;
    stream >> parsed;
    if(stream.fail()){
      throw std::runtime_error("Invalid value for input parameter '" + label + "'");
    }
    value = parsed;
    lines[label] = value;
    return true;
  }
  
  void print(){
    for(auto a : lines){
      std::cout << std::left << std::setw(25) << a.first << " " << a.second << std::endl;
    }
    for(auto a : vlines){
      std::cout << std::left << std::setw(25) << a.first << " : ";
      for(double x : a.second ) std::cout << x << " ";
      std::cout << std::endl;
    }
    for(auto a : blines){
      std::cout << std::left << a.first << " : ";
      for(double x : a.second ) std::cout << x << " ";
      std::cout << std::endl;
    }
    time_t now = time(0);
    std::cout << "log has existed for " << difftime(now,to)/60 << " min" << std::endl;
  }
  
  void print_to_file(){
    print_file(filename);
  }
  
  void print_file(std::string filename){
    //printrow(std::cout, it->first, it->second, comment->second);
 
    std::ofstream logfile(filename);
    time_t now = time(0);
    struct tm date = *localtime(&now);
    logfile << date.tm_hour << ":" << date.tm_min << "   " << date.tm_mday << "/" << date.tm_mon << "/" << date.tm_year + 1900 << std::endl;
    logfile << "log has existed for " << difftime(now,to)/60 << " min" << std::endl;
    for(auto a : lines){
      logfile << std::left << std::setw(23) << a.first << " " << a.second << std::endl;
    }
    for(auto a : vlines){
      logfile << std::left << std::setw(23) << a.first << " : ";
      for(double x : a.second ) logfile << x << " ";
      logfile << std::endl;
    }
    for(auto a : blines){
      logfile << std::left << a.first << " : ";
      for(double x : a.second ) logfile << x << " ";
      logfile << std::endl;
    }
  }
 
};
 
 
class LOGDATA{
 
public:
  typedef boost::variant<int,long,size_t,float,double,std::string,bool> MULTITYPE;
  typedef std::map<std::string,MULTITYPE> LINE;
 
private:
  std::vector<LINE> lines;
  std::string filename;
  std::string blank_val;
  std::vector<std::string> header;
  long nbatch;
  int precision;
  long last_line_printed;
  long nlabels;
  std::set<std::string> labels;
  
  std::map<std::string,std::string> label_comments;
  
  long nprints;
  
  void append_file(){
    
    const long default_precision = std::cout.precision();
    std::ofstream logfile;
    logfile.open(filename,std::ios_base::app);
   
    size_t n=labels.size();
    if(n == nlabels && last_line_printed > 0){
      std::cout << std::setprecision(precision);
      for(size_t j=last_line_printed ; j<lines.size() ; ++j ){
        int i=0;
        for(auto &label : labels){
          try{
            if(label == "ID" || label == "id" 
              || label == "galaxy_halo_id" ){
              logfile << std::setprecision(17) << lines[j].at(label);
              std::cout << std::setprecision(precision);
            }else if(label == "RA" || label == "DEC" ){
                logfile << std::setprecision(8) << lines[j].at(label)  ;
                std::cout << std::setprecision(precision);
            }else{
              logfile << std::setprecision(precision) << lines[j].at(label) ;
            }
            if(i<n-1) logfile << ",";
          }catch(std::exception& e){
            logfile << std::setprecision(precision) << blank_val ;
            if(i<n-1) logfile << ",";
          }
          ++i;
        }
        logfile << std::endl;
      }
      std::cout << std::setprecision(default_precision);
      
      last_line_printed = lines.size();
    }else{
      // number of labels has changed print from the beginning
      print_to_file();
    }
  }
  
public:
 
  LOGDATA(std::string file
          ,std::string blanck_value = "0"
          ,long Nbatch=1000
          ):filename(file),blank_val(blanck_value)
  ,nbatch(Nbatch),precision(std::cout.precision())
  ,last_line_printed(0),nlabels(0),nprints(0){
    print_to_file();
  };
  
  ~LOGDATA(){print_to_file();}
  
  std::string output_file(){return filename;}
  
  void set_precision(int p){
    precision = p;
  }
  
  // current number of columns
  size_t ncol(){return labels.size();}
  // names of columns
  std::set<std::string> names = labels;
  
  void print_to_file(){
    
    const int default_precision = std::cout.precision();
    //std::cout << std::setprecision(12);
    
    if(lines.size() == 0 ) return;
    
    std::ofstream logfile(filename);
    
    time_t now = time(0);
    struct tm date = *localtime(&now);
    logfile << "# " << date.tm_hour << ":" << date.tm_min << "   " << date.tm_mday << "/" << date.tm_mon << "/" << date.tm_year + 1900 << std::endl;
    for(std::string &s : header){
      logfile << "# " << s << std::endl;
    }
 
    // print comments on labels
    {
      int i=1;
      for(auto &label : labels){
        auto it = label_comments.find(label);
        if(it == label_comments.end() ){
          logfile << "# " << i << " " << label << " : no comment" << std::endl;
        }else{
          logfile << "# " << i << " " << it->first << " : " << it->second << std::endl;
        }
        ++i;
      }
    }
    size_t i=0;
    nlabels = labels.size();
    for(auto &label : labels){
      logfile << label;
      if(i<nlabels-1) logfile << ",";
      ++i;
    }
    logfile << std::endl;
 
    std::cout << std::setprecision(precision);
    for(LINE &line : lines){
      i=0;
      for(auto &label : labels){
        try{
          if(label == "ID"){
            logfile << std::setprecision(17) << line.at(label)  ;
            std::cout << std::setprecision(precision);
          }else if(label == "RA" || label == "DEC" ){
            logfile << std::setprecision(8) << line.at(label)  ;
            std::cout << std::setprecision(precision);
          }else{
            logfile  << std::setprecision(precision) << line.at(label);
          }
          if(i<nlabels-1) logfile << ",";
        }catch(std::exception& e){
          logfile << std::setprecision(precision) << blank_val;
          if(i<nlabels-1) logfile << ",";
        }
        ++i;
      }
      logfile << std::endl;
    }
    std::cout << std::setprecision(default_precision);
    
    last_line_printed = lines.size();
    ++nprints;
  }
  
  void add(const LINE &line){
    for (auto itr = line.begin(); itr != line.end(); ++itr) labels.insert(itr->first);
    lines.push_back(line);
    if(lines.size() % nbatch == 0) append_file();
  }
  
  void replace(const LINE &line){
    for (auto itr = line.begin(); itr != line.end(); ++itr) labels.insert(itr->first);
    lines.back() = line;
  }
 
  void header_comment(const std::string &s){
    header.push_back(s);
  }
  
  void label_description(const std::string &name
                     ,const std::string &description
                     ){
    label_comments[name] = description;
  }
};
 
//  /// alternative to std::thread that counts threads
//  template<class Function, class ... Args >
//  class MyThread
//  {
//  private:
//    std::thread thread_;
//    
//    static long count;
//    static long unjoined;
//  
//  public:
//    MyThread(Function&& f,Args&&... args)
//    {
//      if(count > 20) throw std::runtime_error("too many threads");
//      thread_ = std::thread(std::forward<Function>(f),std::forward<Args>(args)...);
//      ++count;
//      ++unjoined;
//    }
//    
//    ~MyThread(){
//      --count;
//    }
//    
//    void join(){
//      thread_.join();
//      --unjoined;
//    }
//    
// 
//  };
 
 
template <typename T>
std::valarray<T> AdaptiveSmooth(const std::valarray<T> &map_in,size_t Nx,size_t Ny,T value){
    
    std::valarray<T> map_out(map_in.size());
    long r,area;
    T val;
    for(long i=0;i<Nx;++i){
      for(long j=0;j<Ny;++j){
        r = 0;
        val = map_in[i+j*Nx];
        while(val < value && r < std::min(Nx, Ny) ){
          
          area = 0;
          val = 0;
          long imin,imax,jmin,jmax;
          
          imin = (i-r < 0) ? 0 : i-r;
          imax = (i+r > Nx-1) ? Nx-1 : i+r;
          
          jmin = (j-r < 0) ? 0 : j-r;
          jmax = (j+r > Ny-1) ? Ny-1 : j+r;
          
          
          for(long ii=imin;ii<=imax;++ii){
            for(long jj=jmin;jj<=jmax;++jj){
              if( (ii-i)*(ii-i) + (jj-j)*(jj-j) < r){
                val += map_in[ii+jj*Nx];
                ++area;
              }
            }
          }
          ++r;
        }
        
        map_out[i+j*Nx] = val/area;
      }
    }
    
    return map_out;
  }
}  // Utilities
 
//#endif
 
 
#endif