basically functional, but still a work in progress

1 files changed, 108 insertions, 0 deletions

diff --git a/poissonier.h b/poissonier.h
new file mode 100644
index 0000000..95da15e
--- /dev/null
+++ b/poissonier.h
@@ -0,0 +1,108 @@
+#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 */