Use CBRNG in schedules.

arbor/backends/multicore/rand.cpp

@@ -17,7 +17,7 @@ void generate_random_numbers(
  for (std::size_t n=0; n<num_rv; ++n) {
  for (std::size_t i=0; i<width; ++i) {
  const auto r = cbprng::generator{}({seed, mech_id, n, counter},
- {gid[i], idx[i], 0xdeadf00dull, 0xdeadbeefull});
+  {gid[i], idx[i], 0xdeadf00dull, 0xdeadbeefull});
  const auto [a0, a1] = r123::boxmuller(r[0], r[1]);
  const auto [a2, a3] = r123::boxmuller(r[2], r[3]);
  dst[i + width_padded*(0 + cbprng::cache_size()*n)] = a0;

arbor/include/arbor/schedule.hpp

@@ -5,7 +5,6 @@
 #include <type_traits>
 #include <utility>
 #include <vector>
-#include <random>
 
 #include <arbor/assert.hpp>
 #include <arbor/common_types.hpp>
@@ -17,10 +16,9 @@
 // Time schedules for probe–sampler associations.
 namespace arb {
 
-using engine_type = std::mt19937_64;
-using seed_type = std::remove_cv_t<decltype(engine_type::default_seed)>;
+using seed_type = std::uint64_t;
 
-constexpr static auto default_seed = engine_type::default_seed;
+constexpr static auto default_seed = 0xdeadbeef;
 
 using time_event_span = std::pair<const time_type*, const time_type*>;
 
@@ -111,13 +109,14 @@ schedule ARB_ARBOR_API regular_schedule(const units::quantity& dt);
 schedule ARB_ARBOR_API explicit_schedule(const std::vector<units::quantity>& seq);
 schedule ARB_ARBOR_API explicit_schedule_from_milliseconds(const std::vector<time_type>& seq);
 
+/// Poisson point process schedule.
 schedule ARB_ARBOR_API poisson_schedule(const units::quantity& tstart,
  const units::quantity& rate,
- seed_type seed = default_seed,
+ seed_type seed=default_seed,
  const units::quantity& tstop=terminal_time*units::ms);
 
 schedule ARB_ARBOR_API poisson_schedule(const units::quantity& rate,
- seed_type seed = default_seed,
+ seed_type seed=default_seed,
  const units::quantity& tstop=terminal_time*units::ms);
 
 } // namespace arb

arbor/schedule.cpp

@@ -4,13 +4,15 @@
 #include <arbor/common_types.hpp>
 #include <arbor/schedule.hpp>
 
+#include "util/cbrng.hpp"
+
 namespace arb {
 
 // Schedule at Poisson point process with rate 1/mean_dt,
 // restricted to non-negative times.
 struct poisson_schedule_impl {
  poisson_schedule_impl(time_type tstart, time_type rate_kHz, seed_type seed, time_type tstop):
- tstart_(tstart), rate_(rate_kHz), exp_(rate_kHz), rng_(seed), seed_(seed), next_(tstart), tstop_(tstop) {
+ tstart_(tstart), rate_(rate_kHz), rng_(seed), seed_(seed), next_(tstart), tstop_(tstop) {
  if (!std::isfinite(tstart_)) throw std::domain_error("Poisson schedule: start must be finite and in [ms]");
  if (!std::isfinite(tstop_)) throw std::domain_error("Poisson schedule: stop must be finite and in [ms]");
  if (!std::isfinite(rate_kHz)) throw std::domain_error("Poisson schedule: rate must be finite and in [kHz]");
@@ -20,9 +22,7 @@ struct poisson_schedule_impl {
  }
 
  void reset() {
- rng_ = engine_type{seed_};
- if (discard_ > 0) rng_.discard(discard_);
- exp_ = std::exponential_distribution<time_type>{rate_};
+ rng_ = {seed_, discard_};
  next_ = tstart_;
  step();
  }
@@ -32,23 +32,27 @@ struct poisson_schedule_impl {
  time_event_span events(time_type t0, time_type t1) {
  // if we start after the maximal allowed time, we have nothing to do
  if (t0 >= tstop_) return {};
-
  // restrict by maximal allowed time
  t1 = std::min(t1, tstop_);
-
- times_.clear();
-
+ // advance to start time
  while (next_ < t0) { step(); }
-
+ // record events in [t0, t1)
+ times_.clear();
  while (next_ < t1) {
  times_.push_back(next_);
  step();
  }
-
  return as_time_event_span(times_);
  }
 
- void step() { next_ += exp_(rng_); }
+ // Leverage that for X ~ P(lambda) the interarrival times T = X_i+1 - X_i
+ // are T ~ Exp(lambda), thus an ordered stream of X ~ P(lambda) is the
+ // cumulative sum of T ~ Exp(lambda). We convert U(0, 1) to Exp(lambda) per
+ // inverse transform sampling. See
+ // https://en.wikipedia.org/wiki/Exponential_distribution#Random_variate_generation
+ // Important detail: as rng gives X in [0, 1) and we'd rather not have log(0), we
+ // use 1 - X in (0, 1] (but still uniform).
+ void step() { next_ += -std::log(1 - rng_())/rate_; }
 
  template<typename K>
  void t_serialize(::arb::serializer& ser, const K& k) const {
@@ -71,8 +75,7 @@ struct poisson_schedule_impl {
 
  time_type tstart_;
  time_type rate_;
- std::exponential_distribution<time_type> exp_;
- engine_type rng_;
+ util::uniform_t rng_;
  seed_type seed_;
  time_type next_;
  std::vector<time_type> times_;

arbor/util/cbrng.hpp

@@ -7,6 +7,35 @@
 namespace arb {
 namespace util {
 
+struct uniform_t {
+ typedef r123::Threefry2x64 rng;
+ rng::key_type key_{{0}};
+ rng::ctr_type ctr_{{0,0}};
+ rng g;
+ double cache_ = -1;
+
+ uniform_t(std::uint64_t seed, std::uint64_t discard=0):
+ key_{{seed}}, ctr_{{discard/2, 0}} {
+ // discard was odd, so we prime the cache by generating one pair, discarding the first
+ if (discard % 2) (*this)();
+ }
+
+ double operator()() {
+ // have one value cached, so return that and invalidate
+ if (cache_ > 0) {
+ return std::exchange(cache_, -1);
+ }
+ // nothing in cache, so generate two new numbers caching the second
+ else {
+ auto rand = g(ctr_, key_);
+ ctr_[0] += 1;
+ cache_ = r123::u01<double>(rand[1]);
+ return r123::u01<double>(rand[0]);
+ }
+ }
+};
+
+inline
 std::vector<double> uniform(uint64_t seed, unsigned left, unsigned right) {
  typedef r123::Threefry2x64 cbrng;
  std::vector<double> r;
@@ -38,4 +67,4 @@ std::vector<double> uniform(uint64_t seed, unsigned left, unsigned right) {
 }
 
 }
-}
+}

example/bench/bench.cpp

@@ -21,6 +21,7 @@
 #include <arborenv/default_env.hpp>
 #include <arborenv/gpu_env.hpp>
 
+#include <random>
 #include <sup/ioutil.hpp>
 #include <sup/json_meter.hpp>
 #include <sup/json_params.hpp>

example/brunel/brunel.cpp

@@ -6,6 +6,7 @@
 #include <optional>
 #include <set>
 #include <vector>
+#include <random>
 
 #include <tinyopt/tinyopt.h>
 

example/drybench/drybench.cpp

@@ -1,5 +1,6 @@
 #include <fstream>
 #include <iostream>
+#include <random>
 
 #include <nlohmann/json.hpp>
 

example/lfp/lfp.cpp

@@ -2,6 +2,7 @@
 #include <cassert>
 #include <vector>
 #include <iostream>
+#include <numeric>
 
 #include <arbor/cable_cell.hpp>
 #include <arbor/morph/morphology.hpp>

python/test/unit/test_schedules.py

@@ -6,6 +6,7 @@
 
 import arbor as A
 from arbor import units as U
+import numpy as np
 
 """
 all tests for schedules (regular, explicit, poisson)
@@ -114,29 +115,22 @@ def test_tstart_freq_seed_contor_poisson_schedule(self):
  self.assertEqual(ps.seed, 1000)
 
  def test_events_poisson_schedule(self):
- expected = [17.4107, 502.074, 506.111, 597.116]
- ps = A.poisson_schedule(tstart=0.0 * U.ms, freq=0.01 * U.kHz, seed=0)
- for i in range(len(expected)):
- self.assertAlmostEqual(
- expected[i], ps.events(0.0 * U.ms, 600.0 * U.ms)[i], places=3
- )
- expected = [
- 5030.22,
- 5045.75,
- 5069.84,
- 5091.56,
- 5182.17,
- 5367.3,
- 5566.73,
- 5642.13,
- 5719.85,
- 5796,
- 5808.33,
- ]
- for i in range(len(expected)):
- self.assertAlmostEqual(
- expected[i], ps.events(5000.0 * U.ms, 6000.0 * U.ms)[i], places=2
- )
+ # Sanity check for Poisson Point Process;
+ # lifted and generalized from underlying C++ unit tests (see there).
+ # Test that the number of events in a fixed bucket conforms to
+ # a Poisson distribution.
+ chi2_lb = 888.56352318146696
+ chi2_ub = 1118.9480663231843
+ T = 1001.0
+ for seed in range(10):
+ pss = A.poisson_schedule(tstart=0.0 * U.ms, freq=8.13 * U.kHz, seed=seed)
+ evs = pss.events(0.0 * U.ms, T * U.ms)
+ pdf = np.zeros(shape=int(T))
+ for t in evs:
+ pdf[int(t)] += 1
+ dsp = T * np.mean(pdf) / np.var(pdf)
+ self.assertLess(dsp, chi2_ub)
+ self.assertGreater(dsp, chi2_lb)
 
  def test_exceptions_poisson_schedule(self):
  with self.assertRaises(TypeError):

scripts/run_cpp_examples.sh

@@ -47,8 +47,8 @@ skip_local=(
 
 # Lookup table for expected spike count
 expected_outputs=(
- 972
- 6998
+ 954
+ 6999
  "30"
  ""
  ""

test/ubench/mech_vec.cpp

@@ -6,7 +6,7 @@
 // will need to be reworked in order to compile.
 
 #include <any>
-#include <fstream>
+#include <random>
 
 #include <arbor/cable_cell.hpp>
 #include <arbor/morph/segment_tree.hpp>

test/unit/test_domain_decomposition.cpp

@@ -1,3 +1,5 @@
+#include <numeric>
+
 #include <gtest/gtest.h>
 
 #include <arbor/context.hpp>

test/unit/test_schedule.cpp

@@ -17,6 +17,8 @@ using namespace testing;
 
 using time_range = util::range<const time_type*>;
 
+using engine_type = std::mt19937_64;
+
 // Pull events from n non-contiguous subintervals of [t0, t1) and check for
 // monotonicity and boundedness.
 void run_invariant_checks(schedule S, time_type t0, time_type t1, unsigned n, int seed=0) {
@@ -205,7 +207,7 @@ TEST(schedule, poisson_uniformity) {
 
  constexpr int N = 1001;
  //constexpr double alpha = 0.01;
- constexpr double chi2_lb = 888.56352318146696;
+ constexpr double chi2_lb =  888.56352318146696;
  constexpr double chi2_ub = 1118.9480663231843;
 
  double dispersion = poisson_schedule_dispersion(N, .813);