testcases, optimized dists excluding noncentral t, new rng

RunDEMC/de.py

@@ -16,6 +16,7 @@ def njit(func):
         return func
 
 import numpy as np
+from numpy.random import default_rng
 import sys
 import random
 import time

RunDEMC/demc.py

@@ -10,6 +10,7 @@
 
 # global imports
 import numpy as np
+from numpy.random import default_rng
 import sys
 import random
 import time
@@ -106,10 +107,16 @@ def __init__(self, name, params, like_fun,
                  pop_recarray=True,
                  pop_parallel=False,
                  n_jobs=-1, backend=None,
-                 use_dask=False):
+                 use_dask=False,
+                 rng=None):
         """
         DEMC
         """
+        self._rng = rng
+
+        if self._rng is None:
+            self._rng = default_rng()
+
         # set the vars
         self._name = name
         self._params = params  # can't be None
@@ -288,7 +295,7 @@ def _get_split_ind(self):
 
         if self._rand_split:
             # randomize the indices
-            all_ind = np.random.permutation(num_particles)
+            all_ind = self._rng.permutation(num_particles)
         else:
             all_ind = np.arange(num_particles)
 
@@ -300,7 +307,7 @@ def _get_split_ind(self):
 
     def _migrate(self):
         # pick which items will migrate
-        num_to_migrate = np.random.random_integers(2, self._num_chains)
+        num_to_migrate = self._rng.integers(2, self._num_chains)
         to_migrate = random.sample(
             list(range(self._num_chains)), num_to_migrate)
 
@@ -328,7 +335,7 @@ def _migrate(self):
             mh_prob = np.exp(log_diff)
             if np.isnan(mh_prob):
                 mh_prob = 0.0
-            keep = (mh_prob - np.random.rand()) > 0.0
+            keep = (mh_prob - self._rng.random()) > 0.0
             if keep:
                 keepers.append({'ind': j,
                                 'particle': self._particles[-1][i],
@@ -362,7 +369,7 @@ def sample(self, num_iter, burnin=False, migration_prob=0.0):
         else:
             progress = range(num_iter)
         for i in progress:
-            if np.random.rand() < migration_prob:
+            if self._rng.random() < migration_prob:
                 # migrate, which is deterministic and done in place
                 # if self._verbose:
                 #     sys.stdout.write('x ')
@@ -401,7 +408,7 @@ def sample(self, num_iter, burnin=False, migration_prob=0.0):
             evo_args = self._get_evo_args(burnin=burnin)
 
             # get new state
-            evo_res = _evolve(**evo_args)
+            evo_res = _evolve(self._rng, **evo_args)
 
             # save the state
             self._save_evolution(**evo_res)
@@ -743,7 +750,7 @@ def _calc_log_prior(*props, priors, cur_split=None):
     return log_priors
 
 
-def _evolve(prop_gen=None, parts_ind=None, split_ind=None,
+def _evolve(rng, prop_gen=None, parts_ind=None, split_ind=None,
             particles=None, param_names=None,
             transforms=None, priors=None, fixed=None,
             log_likes=None, weights=None, use_priors=True, recalc_likes=False,
@@ -812,15 +819,15 @@ def _evolve(prop_gen=None, parts_ind=None, split_ind=None,
             prev_weights = prev_log_likes
 
         # calc acceptance in log space
-        keep = np.log(np.random.rand(len(log_diff))) < log_diff
+        keep = np.log(rng.random(len(log_diff))) < log_diff
 
         # handle much greater than one
         #log_diff[log_diff > 0.0] = 0.0
 
         # now exp so we can get the other probs
         # mh_prob = np.exp(log_diff)
         # mh_prob[np.isnan(mh_prob)] = 0.0
-        # keep = (mh_prob - np.random.rand(len(mh_prob))) > 0.0
+        # keep = (mh_prob - rng.random(len(mh_prob))) > 0.0
 
         # modify the relevant particles
         # use mask the mask approach
@@ -891,7 +898,7 @@ def rvs(self, size, cur_split):
             chains = np.arange(self._num_chains)
         else:
             # pick randomly
-            chains = np.random.randint(0, self._num_chains, size[0])
+            chains = self._rng.integers(0, self._num_chains, size[0])
 
         # generate the random vars using the likelihood func
         # pop = self._particles[-1][chains]
@@ -910,7 +917,8 @@ def __init__(self, name, dist, params,
                  proposal_gen=None,
                  burnin_proposal_gen=None,
                  use_priors=True,
-                 verbose=False):
+                 verbose=False,
+                 rng=None):
 
         # handle the dist
         self._dist = dist
@@ -925,7 +933,8 @@ def __init__(self, name, dist, params,
                                          initial_zeros_ok=True,
                                          use_priors=use_priors,
                                          verbose=verbose,
-                                         pop_recarray=False)
+                                         pop_recarray=False,
+                                         rng=rng)
 
         # We need to recalc likes of prev iteration in case of hyperprior
         self._recalc_likes = True
@@ -967,7 +976,7 @@ def pdf(self, vals, cur_split):
                 cur_split = np.arange(self._num_chains)
             else:
                 # pick randomly
-                cur_split = np.random.randint(0, self._num_chains, len(vals))
+                cur_split = self._rng.integers(0, self._num_chains, len(vals))
 
         # generate the pdf using the likelihood func
         pop = _apply_param_transform(self._particles[-1][cur_split],
@@ -997,7 +1006,7 @@ def logpdf(self, vals, cur_split):
                 cur_split = np.arange(self._num_chains)
             else:
                 # pick randomly
-                cur_split = np.random.randint(0, self._num_chains, len(vals))
+                cur_split = self._rng.integers(0, self._num_chains, len(vals))
 
         # generate the pdf using the likelihood func
         pop = _apply_param_transform(self._particles[-1][cur_split],
@@ -1027,7 +1036,7 @@ def rvs(self, size, cur_split):
             chains = np.arange(num_chains)[cur_split]
         else:
             # pick randomly
-            chains = np.random.randint(0, num_chains, size[0])
+            chains = self._rng.integers(0, num_chains, size[0])
 
         # generate the random vars using the likelihood func
         # pop = self._particles[-1][chains]