#ifndef POISSONIER_H
#define POISSONIER_H

#include <random>       /* for seed_seq */
#include <iostream>
#include <time.h>       /* for clock_gettime */
#include <vector>
#ifdef BOOST_CPD
#  include <boost/math/distributions/poisson.hpp>
#endif

static char const printable[64+1] = "ABCDEFGHIJKLMNOPQRSTUVWXYZ"
                 "abcdefghijklmnopqrstuvwxyz"
                 "0123456789._";

std::string base64(unsigned long int arg) {
  std::string str;
  for (;;) {
    str = printable[arg & 63] + str;
    arg >>= 6;
    if (! arg) break;
  }
  return str;
}

std::string randomstr(unsigned int const nn=15) {
  using namespace std;
// FIXME: ??? Could maybe use std::bitset here:
  random_device rd;
  string str;
  while (str.length() < nn) {
    unsigned int foo(rd());
    int nch = (sizeof(foo) * 8) / 6;
    for (int ii = 0; ii < nch; ii++) {
      str += printable[foo & 63];
      foo >>= 6;
    }
  }
  return str;
}

class poissonier {
public:
  std::mt19937 basegen;
  std::string graine;           // last seed fed to basegen

// constructor
  poissonier () {
    timespec time;
    clock_gettime(CLOCK_REALTIME, &time);
    graine = base64(time.tv_sec);
    graine += base64(time.tv_nsec);
    seed(graine);
  }
  double sample(double const lambda) {
// There is no way to change the intensity after the
// distro is constructed, so construct a new one each time:
    return std::poisson_distribution<>(lambda)(basegen);
  }
  void seed(std::string const seed) {
// Need a temporary variable here, because the seed() function
// modifies it:
    graine = seed;
    std::seed_seq seed1(seed.begin(), seed.end());
    basegen.seed(seed1);
  }
  std::string randomize() {
    graine = randomstr(15);
    seed(graine);
    return graine;
  }
  void check(double const rate) {
    using namespace std;
    cout << "Checking: " << rate
        << "  gives: " << sample(rate)
        << endl;
  }

// Note this is stateless. Static member function.
//
// Log of the probability mass distribution:
  static double lpmd(double const lambda, int const k){
    // https://en.wikipedia.org/wiki/Poisson_distribution#Definition
    if (lambda < 0) throw std::invalid_argument("negative Poisson intensity");
    if (lambda == 0) {
      if (k == 0) return 0;     // 100% probability
      return -INFINITY;         //   0% probability
    }
    return k * log(lambda) - lambda - lgamma(k + 1.0);
  }

// pmd = probability mass distribution
  inline static double pmd(double const lambda, int const k){
    return exp(lpmd(lambda, k));
  }

#ifdef BOOST_CPD
// not implemented
// not needed, and would require boost
// cpd = cumulative probability distribution
  static double cpd(double const lambda, int const k){
    boost::math::poisson_distribution<> pd(lambda);
    return boost::math::cpd(pd, k);
  }
#endif
};

#endif /* POISSONIER_H */