Repository: starsimhub/starsim (language: python)

Issue:
Dist tidying
- [x] Check performance (should be ~30% gains at least)
- [x] Reimplement slots if needed
- [x] Write tests to more carefully check variance reduction
Explore "hash & cache" RNG
Potentially a way to generate random numbers without having to generate extra.
- [x] Replace `pandas` with more direct/faster approach
- [x] Combine single and multi RNG
- [x] Don't remove slots / consider factor of 5 difference
- [x] Simplify/standardize hashing step
- [x] Enable alternatives to SciPy (eg numpy RNG, scikit-learn)
- [x] Simplify `ScipyDistribution`
- [x] Simplify nonstandard distributions with ICDFs

See also https://github.com/amath-idm/starsim/issues/299

New checklist:
- [x] Reimplement `test_random.py`
- [x] Reimplement `test_rngcontainer.py`
- [x] Reimplement `test_distributions.py`
- [x] Get most basic flow working
- [x] Re-enable SciPy distributions
- [x] Define explicit distributions
- [x] Figure out Bernoulli distributions
- [x] Figure out delta distributions
- [x] Refactor code with new API
- [x] Get basic results to match
- [x] Get tests working
- [ ] Check performance

Even newer:
- [x] Figure out difference in indices in demo
- [ ] Figure out slots
- [x] Figure out what lognorm/lognorm_mean are doing
- [x] Figure out rv_discrete replacement
- [x] Figure out SciPyHistogram
- [ ] Figure out why sc.prettyobj breaks Dists
- [ ] Figure out why microsim breaks repr
Dist tidying
- [x] Check performance (should be ~30% gains at least)
- [x] Reimplement slots if needed
- [x] Write tests to more carefully check variance reduction
Refactor rngs ppf
Refactors the RNG refactor to use slots, and to use `ppf` as an alternative solution to the repeated slot handling code.

SUPER ugly code, needs a LOT of tidying, but it seems like it works, maybe?!

In almost all cases, the "standard" solution works fine: take the basic method (e.g. `rng.uniform(low, high, size)`), where `size=n` if not using UIDs, and `size=sim.people.slot[uids].max()+1` if using UIDs.

The only case in which this _doesn't_ work is if there are repeated slots _and_ array (or callable) parameters. Now what happens in this case is that:
1. The number of random numbers drawn are as before: `size = sim.people.slot[uids].max()+1`
2. The _quantile_ for each slot is drawn: `qq = rng.random(sze)`
3. The slots in use are pulled out: `qqs = qq[slots]`
4. These are then mapped to the actual distribution via the inverse CDF, aka PPF, e.g.: `rvs = scipy.stats.uniform(low, high).ppf(qqs)`

Two agents who share the same slot will have exactly the same quantile, but they can still map to different values since the distribution parameter values can be different. Even though calling `ppf` is typically slower than generating new random numbers, since this is only for a corner case, and since `ppf` is only called for the slots that are in use, it should be as fast or faster than the other implementation.

For known issues, see https://github.com/starsimhub/starsim/issues/422


Discussion / hints:

Is there a specific generator you have in mind here, or a more general/efficient approach to using any RNG?



Fix patch:
diff --git a/CHANGELOG.rst b/CHANGELOG.rst
index d866573c..ede4065c 100644
--- a/CHANGELOG.rst
+++ b/CHANGELOG.rst
@@ -7,6 +7,26 @@ What's new
 All notable changes to the codebase are documented in this file. Changes that may result in differences in model output, or are required in order to run an old parameter set with the current version, are flagged with the term "Regression information".
 
 
+Version 0.3.0 (2024-03-30)
+--------------------------
+
+New RNGs & distributions
+~~~~~~~~~~~~~~~~~~~~~~~~
+- Replaces ``ss.SingleRNG()``, ``ss.MultiRNG()``, ``ss.ScipyDistribution()``, and ``ss.ScipyHistogram()`` with a single ``ss.Dist()`` class. The ``starsim.random`` and ``starsim.distributions`` submodules have been removed, and ``starsim.dists`` has been added.
+- The ``ss.Dist`` class uses ``np.random.default_rng()`` rather than ``scipy.stats`` by default, although a ``scipy.stats`` distribution can be supplied as an alternative. This is up to 4x faster (including, critically, for Bernoulli distributions).
+- Also removes ``ss.options.multirng`` (the new version is equivalent to it being always on).
+- Removes duplicate logic for transmission (``make_new_cases()``)
+- Adds new custom distributions such as ``ss.choice()`` and ``ss.delta()``.
+- These distributions can be called directly, e.g. ``dist = ss.weibull(c=2); dist(5)`` will return 5 random variates from a Weibull distribution.
+- Instead of being manually initialized based on the name, the ``Sim`` object is parsed and all distributions will be initialized with a unique identifier based on their place in the object (e.g. ``sim.diseases.sir.pars.dur_inf``), which is used to set their unique seed.
+
+Other changes
+~~~~~~~~~~~~~
+- This PR also fixes bugs with lognormal parameters, and makes it clear whether the parameters are for the *implicit* normal distribution (``ss.lognorm_im()``, the NumPy/SciPy default, equivalent to ``ss.lognorm_mean()`` previously) or the "explicit" lognormal distribution (``ss.lognorm_ex()``, equivalent to ``ss.lognorm()`` previously).
+- Renames ``ss.dx``, ``ss.tx``, ``ss.vx`` to``ss.Dx``, ``ss.Tx``, ``ss.Vx``.
+- Removed ``set_numba_seed()`` as a duplicate of ``set_seed()``.
+- *GitHub info*: PR `392 <https://github.com/amath-idm/stisim/pull/392>`_
+
 Version 0.2.10 (2024-03-18)
 ---------------------------
 - SIR duration of infection now accounts for dt
diff --git a/docs/api/index.rst b/docs/api/index.rst
index f299707d..e050cccd 100644
--- a/docs/api/index.rst
+++ b/docs/api/index.rst
@@ -22,7 +22,6 @@ This page lists Starsim's full API, including all functions and classes.
    starsim.parameters
    starsim.people
    starsim.products
-   starsim.random
    starsim.results
    starsim.run
    starsim.samples
diff --git a/setup.py b/setup.py
index ca628f78..a3866e42 100644
--- a/setup.py
+++ b/setup.py
@@ -19,7 +19,6 @@
     "Programming Language :: Python",
     "Topic :: Software Development :: Libraries :: Python Modules",
     "Development Status :: 5 - Production/Stable",
-    "Programming Language :: Python :: 3.8",
     "Programming Language :: Python :: 3.9",
     "Programming Language :: Python :: 3.10",
     "Programming Language :: Python :: 3.11",
@@ -37,11 +36,11 @@
     classifiers=CLASSIFIERS,
     packages=find_packages(),
     include_package_data=True,
-    install_requires=[
+    install_requires=[ # NB: remember to update __init__.py if these requirements change
         'numpy',
         'scipy',
         'pandas>=2.0.0',
-        'sciris>=3.1.0',
+        'sciris>=3.1.5',
         'matplotlib',
         'numba',
         'networkx',
diff --git a/starsim/__init__.py b/starsim/__init__.py
index 2d48c92b..f86aaa28 100644
--- a/starsim/__init__.py
+++ b/starsim/__init__.py
@@ -4,7 +4,6 @@
 from .utils         import *
 from .distributions import *
 from .states        import *
-from .random        import *
 from .people        import *
 from .modules       import *
 from .network       import *
@@ -20,9 +19,14 @@
 from .run           import *
 from .samples       import *
 
+# Assign the root folder
 import sciris as sc
 root = sc.thispath(__file__).parent
 
 # Import the version and print the license
 if options.verbose:
-    print(__license__)
\ No newline at end of file
+    print(__license__)
+
+# Double-check key requirements -- should match setup.py
+sc.require(['sciris>=3.1.5', 'pandas>=2.0.0', 'scipy', 'numba', 'networkx'], message=f'The following dependencies for Starsim {__version__} were not met: <MISSING>.')
+del sc # Don't keep this in the module
\ No newline at end of file
diff --git a/starsim/demographics.py b/starsim/demographics.py
index 717e12a6..ac71066a 100644
--- a/starsim/demographics.py
+++ b/starsim/demographics.py
@@ -9,6 +9,7 @@
 
 __all__ = ['BaseDemographics', 'Births', 'Deaths', 'Pregnancy']
 
+
 class BaseDemographics(ss.Module):
     # A demographic module typically handles births/deaths/migration and takes
     # place at the start of the timestep, before networks are updated and before
@@ -27,6 +28,7 @@ def update(self, sim):
         # carried out inside this function
         pass
 
+
 class Births(BaseDemographics):
     def __init__(self, pars=None, metadata=None, **kwargs):
         super().__init__(pars, **kwargs)
@@ -237,6 +239,7 @@ def finalize(self, sim):
         self.results['cmr'] = 1/self.pars.units*np.divide(self.results['new'], sim.results['n_alive'], where=sim.results['n_alive']>0)
         return
 
+
 class Pregnancy(BaseDemographics):
 
     def __init__(self, pars=None, par_dists=None, metadata=None, **kwargs):
@@ -275,7 +278,7 @@ def __init__(self, pars=None, par_dists=None, metadata=None, **kwargs):
             data_cols = dict(year='Time', age='AgeGrp', value='ASFR'),
         )
 
-        self.choose_slots = ss.randint(low=0, high=1) # Low and high will be reset upon initialization
+        self.choose_slots = ss.randint() # Low and high will be reset upon initialization
 
         # Process data, which may be provided as a number, dict, dataframe, or series
         # If it's a number it's left as-is; otherwise it's converted to a dataframe
@@ -341,8 +344,9 @@ def standardize_fertility_data(self):
 
     def initialize(self, sim):
         super().initialize(sim)
-        self.choose_slots.kwds['low'] = sim.pars['n_agents']+1
-        self.choose_slots.kwds['high'] = int(sim.pars['slot_scale']*sim.pars['n_agents'])
+        low = sim.pars.n_agents+1 # TODO: or 0?
+        high = int(sim.pars.slot_scale*sim.pars.n_agents)
+        self.choose_slots.set(low=low, high=high)
         return
 
     def init_results(self, sim):
diff --git a/starsim/disease.py b/starsim/disease.py
index 02fe1e09..8ac42a39 100644
--- a/starsim/disease.py
+++ b/starsim/disease.py
@@ -188,10 +188,8 @@ def __init__(self, *args, **kwargs):
             ss.State('ti_infected', int, ss.INT_NAN),
         )
 
-        if ss.options.multirng:
-            # Used only in _make_new_cases_multirng
-            self.rng_acquisition = ss.uniform()
-            self.rng_source = ss.uniform()
+        self.rng_target = ss.random()
+        self.rng_source = ss.random()
 
         return
 
@@ -231,36 +229,43 @@ def init_results(self, sim):
         return
 
     def _check_betas(self, sim):
-        """ Check that there's a network for each beta keys """
-        # Convert keys to lower case
-        self.pars.beta = {k.lower(): v for k, v in self.pars.beta.items()}
+        """ Check that there's a network for each beta key """
+        # Ensure keys are lowercase
+        if isinstance(self.pars.beta, dict): # TODO: check if needed
+            self.pars.beta = {k.lower(): v for k, v in self.pars.beta.items()}
 
+        # Create a mapping between beta and networks, and populate it
         betapars = self.pars.beta
-
         betamap = sc.objdict()
         netkeys = list(sim.networks.keys())
-        for bkey in betapars.keys():
-            orig_bkey = bkey[:]
-            if bkey in netkeys: # TODO: CK: could tidy up logic
-                betamap[bkey] = betapars[orig_bkey]
-            else:
-                if 'net' not in bkey:
-                    bkey += 'net'  # Add 'net' suffix if not already there
-                if bkey in netkeys:
+        if netkeys: # Skip if no networks
+            for bkey in betapars.keys():
+                orig_bkey = bkey[:]
+                if bkey in netkeys: # TODO: CK: could tidy up logic
                     betamap[bkey] = betapars[orig_bkey]
                 else:
-                    errormsg = f'No network for beta parameter "{bkey}"; your beta should match network keys:\n{sc.newlinejoin(netkeys)}'
-                    raise ValueError(errormsg)
+                    if 'net' not in bkey:
+                        bkey += 'net'  # Add 'net' suffix if not already there
+                    if bkey in netkeys:
+                        betamap[bkey] = betapars[orig_bkey]
+                    else:
+                        errormsg = f'No network for beta parameter "{bkey}"; your beta should match network keys:\n{sc.newlinejoin(netkeys)}'
+                        raise ValueError(errormsg)
         return betamap
 
-    def _make_new_cases_singlerng(self, sim):
-        # Not common-random-number-safe, but more efficient for when not using the multirng feature
+    def make_new_cases(self, sim):
+        """
+        Add new cases of module, through transmission, incidence, etc.
+        
+        Common-random-number-safe transmission code works by mapping edges onto
+        slots.
+        """
         new_cases = []
         sources = []
         people = sim.people
         betamap = self._check_betas(sim)
 
-        for nkey, net in sim.networks.items():
+        for i,nkey,net in sim.networks.enumitems():
             if not len(net):
                 break
 
@@ -268,139 +273,46 @@ def _make_new_cases_singlerng(self, sim):
             contacts = net.contacts
             rel_trans = (self.infectious & people.alive) * self.rel_trans
             rel_sus = (self.susceptible & people.alive) * self.rel_sus
-            for a, b, beta in [[contacts.p1, contacts.p2, nbetas[0]],
-                               [contacts.p2, contacts.p1, nbetas[1]]]:
+            p1p2b0 = [contacts.p1, contacts.p2, nbetas[0]]
+            p2p1b1 = [contacts.p2, contacts.p1, nbetas[1]]
+            for src, trg, beta in [p1p2b0, p2p1b1]:
 
                 # Skip networks with no transmission
                 if beta == 0:
                     continue
 
                 # Calculate probability of a->b transmission.
-                beta_per_dt = net.beta_per_dt(disease_beta=beta, dt=people.dt)
-                p_transmit = rel_trans[a] * rel_sus[b] * beta_per_dt
-
-                new_cases_bool = np.random.random(
-                    len(a)) < p_transmit  # As this class is not common-random-number safe anyway, calling np.random is perfectly fine!
-                new_cases.append(b[new_cases_bool])
-                sources.append(a[new_cases_bool])
+                beta_per_dt = net.beta_per_dt(disease_beta=beta, dt=people.dt) # TODO: should this be sim.dt?
+                p_transmit = rel_trans[src] * rel_sus[trg] * beta_per_dt
+
+                slots_s = people.slot[src] # Slots for the possible source
+                slots_t = people.slot[trg] # Slots for the possible target
+                # slots = (slots_s + slots_t) % sim.pars.n_agents # CK: faster, but some pairs would get the same draws
+                # rvs = self.rng_source.rvs(size=np.max(slots)+1)[slots]
+                rvs_s = self.rng_source.rvs(n=np.max(slots_s) + 1)[slots_s]
+                rvs_t = self.rng_target.rvs(n=np.max(slots_t) + 1)[slots_t]
+                
+                # Generate a new random number based on the two other random numbers -- 3x faster than `rvs = np.remainder(rvs_s + rvs_t, 1)`
+                rvs = rvs_s + rvs_t
+                inds = np.where(rvs>1.0)[0]
+                rvs[inds] -= 1
+                
+                new_cases_bool = rvs < p_transmit
+                new_cases.append(trg[new_cases_bool])
+                sources.append(src[new_cases_bool])
+        
+        # Tidy up
         if len(new_cases) and len(sources):
-            return np.concatenate(new_cases), np.concatenate(sources)
-        return np.empty((0,), dtype=int), np.empty((0,), dtype=int)
-
-    def _make_new_cases_multirng(self, sim):
-        """
-        Common-random-number-safe transmission code works by computing the
-        probability of each _node_ acquiring a case rather than checking if each
-        _edge_ transmits.
-        Subsequent step uses a roulette wheel with slotted RNG to determine
-        infection source.
-        """
-        people = sim.people
-        n = len(people.uid)  # TODO: possibly could be shortened to just the people who are alive
-        p_acq_node = np.zeros(n)
-        betamap = self._check_betas(sim)
-
-        avec = []
-        bvec = []
-        pvec = []
-        for nkey, net in sim.networks.items():
-            if not len(net):
-                break
-            nbetas = betamap[nkey]
-            contacts = net.contacts
-            rel_trans = self.rel_trans * (self.infectious & people.alive)
-            rel_sus = self.rel_sus * (self.susceptible & people.alive)
-
-            p1p2 = ['p1', 'p2', nbetas[0]]
-            p2p1 = ['p2', 'p1', nbetas[1]]
-            for source, target, beta in [p1p2, p2p1]:  # Transmission from a --> b
-                if beta == 0:
-                    continue
-
-                a, b, beta_arr = contacts[source], contacts[target], contacts.beta
-                nzi = (rel_trans[a] > 0) & (rel_sus[b] > 0) & (beta_arr > 0)
-                avec.append(a[nzi])
-                bvec.append(b[nzi])
-
-                beta_per_dt = net.beta_per_dt(disease_beta=beta, dt=people.dt, uids=nzi)
-                trans_arr = rel_trans[a[nzi]]
-                sus_arr = rel_sus[b[nzi]]
-                new_pvec = trans_arr * sus_arr * beta_per_dt
-                pvec.append(new_pvec)
-
-        if len(avec):
-            dfp1 = np.concatenate(avec)
-            dfp2 = np.concatenate(bvec)
-            dfp = np.concatenate(pvec)
+            new_cases = np.concatenate(new_cases)
+            sources = np.concatenate(sources)
         else:
-            return np.empty((0,), dtype=int), np.empty((0,), dtype=int)
-
-        if len(dfp) == 0:
-            return np.empty((0,), dtype=int), np.empty((0,), dtype=int)
-
-        p2uniq, p2idx, p2inv, p2cnt = np.unique(dfp2, return_index=True, return_inverse=True, return_counts=True)
-
-        # Pre-draw random numbers
-        slots = people.slot[p2uniq]  # Slots for the possible infectee
-        r = self.rng_acquisition.rvs(size=np.max(slots) + 1)
-        q = self.rng_source.rvs(size=np.max(slots) + 1)
-
-        # Now address nodes with multiple possible infectees
-        degrees = np.unique(p2cnt)
-        new_cases = []
-        sources = []
-        for deg in degrees:
-            if deg == 1:
-                # UIDs that only appear once
-                cnt1 = p2cnt == 1
-                uids = p2uniq[cnt1]
-                idx = p2idx[cnt1]
-                p_acq_node = dfp[idx]
-                cases = r[people.slot[uids]] < p_acq_node
-                if cases.any():
-                    s = dfp1[idx][cases]
-            else:
-                dups = np.argwhere(p2cnt==deg).flatten()
-                uids = p2uniq[dups]
-                inds = [np.argwhere(np.isin(p2inv, d)).flatten() for d in dups]
-                probs = dfp[inds]
-                p_acq_node = 1-np.prod(1-probs, axis=1)
-
-                cases = r[people.slot[uids]] < p_acq_node
-                if cases.any():
-                    # Vectorized roulette wheel
-                    cumsum = probs[cases].cumsum(axis=1)
-                    cumsum /= cumsum[:,-1][:,np.newaxis]
-                    ix = np.argmax(cumsum >= q[people.slot[uids[cases]]][:,np.newaxis], axis=1)
-                    s = np.take_along_axis(dfp1[inds][cases], ix[:,np.newaxis], axis=1).flatten()#dfp1[inds][cases][np.arange(len(cases)),ix]
-
-            if cases.any():
-                new_cases.append(uids[cases])
-                sources.append(s)
-
-        if len(new_cases) == 0:
-            return np.empty((0,), dtype=int), np.empty((0,), dtype=int)
-
-        new_cases = np.concatenate(new_cases)
-        sources = np.concatenate(sources)
-        return new_cases, sources
-
-    def make_new_cases(self, sim):
-        """ Add new cases of module, through transmission, incidence, etc. """
-        if not sim.networks:
-            warnmsg = f'Disease {self.name} does not transmit without a network.'
-            if sim.ti == 0: ss.warn(warnmsg, die=False)
-            return
-
-        if ss.options.multirng:
-            # Determine new cases for multirng
-            new_cases, sources = self._make_new_cases_multirng(sim)
-        else:
-            # Determine new cases for singlerng
-            new_cases, sources = self._make_new_cases_singlerng(sim)
-
+            new_cases = np.empty(0, dtype=int)
+            sources = np.empty(0, dtype=int)
+            
         if len(new_cases):
             self._set_cases(sim, new_cases, sources)
+            
+        return new_cases, sources
 
     def _set_cases(self, sim, target_uids, source_uids=None):
         congenital = sim.people.age[target_uids] <= 0
@@ -414,13 +326,14 @@ def _set_cases(self, sim, target_uids, source_uids=None):
     def set_congenital(self, sim, target_uids, source_uids=None):
         pass
 
-
     def update_results(self, sim):
         super().update_results(sim)
         res = self.results
-        res['prevalence'][sim.ti] = res.n_infected[sim.ti] / np.count_nonzero(sim.people.alive)
-        res['new_infections'][sim.ti] = np.count_nonzero(self.ti_infected == sim.ti)
-        res['cum_infections'][sim.ti] = np.sum(res['new_infections'][:sim.ti])
+        ti = sim.ti
+        res.prevalence[ti] = res.n_infected[ti] / np.count_nonzero(sim.people.alive)
+        res.new_infections[ti] = np.count_nonzero(self.ti_infected == ti)
+        res.cum_infections[ti] = np.sum(res['new_infections'][:ti])
+        return
 
 
 class InfectionLog(nx.MultiDiGraph):
diff --git a/starsim/diseases/cholera.py b/starsim/diseases/cholera.py
index 3993e02b..fe071e26 100644
--- a/starsim/diseases/cholera.py
+++ b/starsim/diseases/cholera.py
@@ -19,10 +19,10 @@ def __init__(self, pars=None, par_dists=None, *args, **kwargs):
 
         pars = ss.omergeleft(pars,
             # Natural history parameters, all specified in days
-            dur_exp2inf = ss.lognorm_mean(mean=2.772, stdev=4.737),  # Calculated from Azman et al. estimates https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3677557/
-            dur_asymp2rec = ss.uniform(loc=1, scale=10),    # From WHO cholera fact sheet, asymptomatic individuals shed bacteria for 1-10 days (https://www.who.int/news-room/fact-sheets/detail/cholera)
-            dur_symp2rec = ss.lognorm_mean(mean=5, stdev=1.8),    # According to Fung most modelling studies assume 5 days of symptoms (https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3926264/), but found a range of 2.9-14 days. Distribution approximately fit to these values
-            dur_symp2dead = ss.lognorm_mean(mean=1, stdev=0.5),   # There does not appear to be differences in timing/duration of mild vs severe disease, but death from severe disease happens rapidly https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5767916/
+            dur_exp2inf = ss.lognorm_ex(mean=2.772, stdev=4.737),  # Calculated from Azman et al. estimates https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3677557/
+            dur_asymp2rec = ss.uniform(low=1, high=10),    # From WHO cholera fact sheet, asymptomatic individuals shed bacteria for 1-10 days (https://www.who.int/news-room/fact-sheets/detail/cholera)
+            dur_symp2rec = ss.lognorm_ex(mean=5, stdev=1.8),    # According to Fung most modelling studies assume 5 days of symptoms (https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3926264/), but found a range of 2.9-14 days. Distribution approximately fit to these values
+            dur_symp2dead = ss.lognorm_ex(mean=1, stdev=0.5),   # There does not appear to be differences in timing/duration of mild vs severe disease, but death from severe disease happens rapidly https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5767916/
             p_death = 0.005,   # Probability of death is typically less than 1% when treated
             p_symp = 0.5,   # Proportion of infected which are symptomatic, mid range of ~25% and 57% estimates from Jaclson et al (https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3795095/) and Nelson et al (https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3842031/), respectively
             asymp_trans = 0.01,    # Reduction in transmission probability for asymptomatic infection, asymptomatic carriers shed 100-1000 times less bacteria than symptomatic carriers (https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3084143/ and https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3842031/). Previous models assume a 10% relative transmissibility (https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4238032/)
@@ -40,10 +40,10 @@ def __init__(self, pars=None, par_dists=None, *args, **kwargs):
         )
 
         par_dists = ss.omergeleft(par_dists,
-            dur_exp2inf    = ss.lognorm,
+            dur_exp2inf    = ss.lognorm_ex,
             dur_asymp2rec  = ss.uniform,
-            dur_symp2rec   = ss.lognorm,
-            dur_symp2dead  = ss.lognorm,
+            dur_symp2rec   = ss.lognorm_ex,
+            dur_symp2dead  = ss.lognorm_ex,
             init_prev      = ss.bernoulli,
             p_death        = ss.bernoulli,
             p_symp         = ss.bernoulli,
@@ -168,16 +168,16 @@ def set_prognoses(self, sim, uids, source_uids=None):
     def make_new_cases(self, sim):
         """ Add indirect transmission """
 
-        p = self.pars
-        r = self.results
+        pars = self.pars
+        res = self.results
 
         # Make new cases via direct transmission
         super().make_new_cases(sim)
 
         # Make new cases via indirect transmission
-        p_transmit = r.env_conc[sim.ti] * p.beta_env
-        p.p_env_transmit.kwds['p'] = p_transmit
-        new_cases = p.p_env_transmit.filter(self.susceptible)
+        p_transmit = res.env_conc[sim.ti] * pars.beta_env
+        pars.p_env_transmit.set(p=p_transmit)
+        new_cases = pars.p_env_transmit.filter(sim.people.uid[self.susceptible]) # TODO: make syntax nicer
         if new_cases.any():
             self.set_prognoses(sim, new_cases, source_uids=None)
         return
diff --git a/starsim/diseases/ebola.py b/starsim/diseases/ebola.py
index 4e68a353..d5376675 100644
--- a/starsim/diseases/ebola.py
+++ b/starsim/diseases/ebola.py
@@ -36,12 +36,12 @@ def __init__(self, pars=None, par_dists=None, *args, **kwargs):
         )
 
         par_dists = ss.omergeleft(par_dists,
-            dur_exp2symp    = ss.lognorm,
-            dur_symp2sev    = ss.lognorm,
-            dur_sev2dead    = ss.lognorm,
-            dur_dead2buried = ss.lognorm,
-            dur_symp2rec    = ss.lognorm,
-            dur_sev2rec     = ss.lognorm,
+            dur_exp2symp    = ss.lognorm_ex,
+            dur_symp2sev    = ss.lognorm_ex,
+            dur_sev2dead    = ss.lognorm_ex,
+            dur_dead2buried = ss.lognorm_ex,
+            dur_symp2rec    = ss.lognorm_ex,
+            dur_sev2rec     = ss.lognorm_ex,
             p_sev           = ss.bernoulli,
             p_death         = ss.bernoulli,
             p_safe_bury     = ss.bernoulli,
diff --git a/starsim/diseases/gonorrhea.py b/starsim/diseases/gonorrhea.py
index c1333cd5..d4b29c68 100644
--- a/starsim/diseases/gonorrhea.py
+++ b/starsim/diseases/gonorrhea.py
@@ -22,14 +22,14 @@ def __init__(self, pars=None, par_dists=None, *args, **kwargs):
 
         # Parameters
         pars = ss.omergeleft(pars,
-            dur_inf_in_days = ss.lognorm(s=0.6, scale=10),  # median of 10 days (IQR 7–15 days) https://sti.bmj.com/content/96/8/556
+            dur_inf_in_days = ss.lognorm_ex(mean=10, stdev=0.6),  # median of 10 days (IQR 7–15 days) https://sti.bmj.com/content/96/8/556
             p_symp = 0.5,  # Share of infections that are symptomatic. Placeholder value
             p_clear = 0.2,  # Share of infections that spontaneously clear: https://sti.bmj.com/content/96/8/556
             init_prev = 0.1,
         )
 
         par_dists = ss.omergeleft(par_dists,
-            dur_inf_in_days = ss.lognorm,
+            dur_inf_in_days = ss.lognorm_ex,
             p_symp          = ss.bernoulli,
             p_clear         = ss.bernoulli,
             init_prev       = ss.bernoulli,
@@ -87,7 +87,7 @@ def set_prognoses(self, sim, target_uids, source_uids=None):
 
         # Set natural clearance
         clear_uids = self.pars.p_clear.filter(target_uids)
-        dur = sim.ti + self.pars['dur_inf_in_days'].rvs(clear_uids)/365/sim.dt # Convert from days to years and then adjust for dt
+        dur = sim.ti + self.pars.dur_inf_in_days.rvs(clear_uids)/365/sim.dt # Convert from days to years and then adjust for dt
         self.ti_clearance[clear_uids] = dur
 
         return
diff --git a/starsim/diseases/hiv.py b/starsim/diseases/hiv.py
index 892e8c14..baad3a43 100644
--- a/starsim/diseases/hiv.py
+++ b/starsim/diseases/hiv.py
@@ -44,7 +44,9 @@ def __init__(self, pars=None, par_dists=None, *args, **kwargs):
     @staticmethod
     def make_death_prob(module, sim, uids):
         p = module.pars
-        return sim.dt * module.death_prob_data / (p.cd4_min - p.cd4_max)**2 *  (module.cd4[uids] - p.cd4_max)**2
+        out = sim.dt * module.death_prob_data / (p.cd4_min - p.cd4_max)**2 *  (module.cd4[uids] - p.cd4_max)**2
+        out = np.array(out)
+        return out
 
     def update_pre(self, sim):
         """ Update CD4 """
diff --git a/starsim/diseases/measles.py b/starsim/diseases/measles.py
index 23d37195..06b377fd 100644
--- a/starsim/diseases/measles.py
+++ b/starsim/diseases/measles.py
@@ -28,8 +28,8 @@ def __init__(self, pars=None, par_dists=None, *args, **kwargs):
         )
 
         par_dists = ss.omergeleft(par_dists,
-            dur_exp   = ss.norm,
-            dur_inf   = ss.norm,
+            dur_exp   = ss.normal,
+            dur_inf   = ss.normal,
             init_prev = ss.bernoulli,
             p_death   = ss.bernoulli,
         )
diff --git a/starsim/diseases/ncd.py b/starsim/diseases/ncd.py
index e05ac976..08d1626c 100644
--- a/starsim/diseases/ncd.py
+++ b/starsim/diseases/ncd.py
@@ -21,7 +21,7 @@ def __init__(self, pars=None):
             initial_risk = ss.bernoulli(p=0.3), # Initial prevalence of risk factors
             #'affection_rate': ss.rate(p=0.1), # Instantaneous rate of acquisition applied to those at risk (units are acquisitions / year)
             dur_risk = ss.expon(scale=10),
-            prognosis = ss.weibull_min(c=2, scale=5), # Time in years between first becoming affected and death
+            prognosis = ss.weibull(c=2, scale=5), # Time in years between first becoming affected and death
         )
 
         super().__init__(ss.omerge(default_pars, pars))
@@ -60,7 +60,7 @@ def update_pre(self, sim):
     def make_new_cases(self, sim):
         new_cases = ss.true(self.ti_affected == sim.ti)
         self.affected[new_cases] = True
-        prog_years = self.pars['prognosis'].rvs(new_cases)
+        prog_years = self.pars.prognosis.rvs(new_cases)
         self.ti_dead[new_cases] = sim.ti + sc.randround(prog_years / sim.dt)
         super().set_prognoses(sim, new_cases)
         return new_cases
diff --git a/starsim/diseases/sir.py b/starsim/diseases/sir.py
index 6f28a433..b682dcb4 100644
--- a/starsim/diseases/sir.py
+++ b/starsim/diseases/sir.py
@@ -2,7 +2,6 @@
 Define example disease modules
 """
 
-import numpy as np
 import pylab as pl
 import starsim as ss
 
@@ -26,7 +25,7 @@ def __init__(self, pars=None, par_dists=None, *args, **kwargs):
         )
 
         par_dists = ss.omergeleft(par_dists,
-            dur_inf   = ss.lognorm,
+            dur_inf   = ss.lognorm_ex,
             init_prev = ss.bernoulli,
             p_death   = ss.bernoulli,
         )
@@ -43,7 +42,6 @@ def __init__(self, pars=None, par_dists=None, *args, **kwargs):
     def update_pre(self, sim):
         # Progress infectious -> recovered
         recovered = ss.true(self.infected & (self.ti_recovered <= sim.ti))
-
         self.infected[recovered] = False
         self.recovered[recovered] = True
 
@@ -94,7 +92,7 @@ def plot(self):
 __all__ += ['sir_vaccine']
 
 
-class sir_vaccine(ss.vx):
+class sir_vaccine(ss.Vx):
     """
     Create a vaccine product that changes susceptible people to recovered (i.e., perfect immunity)
     """
diff --git a/starsim/diseases/syphilis.py b/starsim/diseases/syphilis.py
index cb6d9176..7c1e7a72 100644
--- a/starsim/diseases/syphilis.py
+++ b/starsim/diseases/syphilis.py
@@ -45,11 +45,11 @@ def __init__(self, pars=None):
         # Parameters
         default_pars = dict(
             # Adult syphilis natural history, all specified in years
-            dur_exposed = ss.lognorm_mean(mean=1 / 12, stdev=1 / 36),  # https://pubmed.ncbi.nlm.nih.gov/9101629/
-            dur_primary = ss.lognorm_mean(mean=1.5 / 12, stdev=1 / 36),  # https://pubmed.ncbi.nlm.nih.gov/9101629/
-            dur_secondary = ss.norm(loc=3.6 / 12, scale=1.5 / 12),  # https://pubmed.ncbi.nlm.nih.gov/9101629/
-            dur_latent_temp = ss.lognorm_mean(mean=1, stdev=6 / 12),  # https://pubmed.ncbi.nlm.nih.gov/9101629/
-            dur_latent_long = ss.lognorm_mean(mean=20, stdev=8),  # https://pubmed.ncbi.nlm.nih.gov/9101629/
+            dur_exposed = ss.lognorm_ex(mean=1 / 12, stdev=1 / 36),  # https://pubmed.ncbi.nlm.nih.gov/9101629/
+            dur_primary = ss.lognorm_ex(mean=1.5 / 12, stdev=1 / 36),  # https://pubmed.ncbi.nlm.nih.gov/9101629/
+            dur_secondary = ss.normal(loc=3.6 / 12, scale=1.5 / 12),  # https://pubmed.ncbi.nlm.nih.gov/9101629/
+            dur_latent_temp = ss.lognorm_ex(mean=1, stdev=6 / 12),  # https://pubmed.ncbi.nlm.nih.gov/9101629/
+            dur_latent_long = ss.lognorm_ex(mean=20, stdev=8),  # https://pubmed.ncbi.nlm.nih.gov/9101629/
             p_latent_temp = ss.bernoulli(p=0.25),  # https://pubmed.ncbi.nlm.nih.gov/9101629/
             p_tertiary = ss.bernoulli(p=0.35),  # https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4917057/
 
@@ -61,8 +61,8 @@ def __init__(self, pars=None):
             #   3: Live birth without syphilis-related complications
             # Source: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5973824/)
             birth_outcomes=sc.objdict(
-                active = ss.rv_discrete(values=([0, 1, 2, 3, 4], [0.125, 0.125, 0.20, 0.35, 0.200])),
-                latent = ss.rv_discrete(values=([0, 1, 2, 3, 4], [0.050, 0.075, 0.10, 0.05, 0.725])),
+                active = ss.choice(a=5, p=np.array([0.125, 0.125, 0.20, 0.35, 0.200])), # Probabilities of active by birth outcome
+                latent = ss.choice(a=5, p=np.array([0.050, 0.075, 0.10, 0.05, 0.725])), # Probabilities of latent
             ),
             birth_outcome_keys=['miscarriage', 'nnd', 'stillborn', 'congenital'],
 
@@ -123,10 +123,9 @@ def update_pre(self, sim):
             self.primary[secondary_from_primary] = False
             self.set_secondary_prognoses(sim, ss.true(secondary_from_primary))
 
-        if ss.options.multirng:
-            # Hack to reset the MultiRNGs in set_secondary_prognoses so that they can be called again in this timestep. TODO: Refactor
-            self.pars.p_latent_temp.random_state.step(sim.ti+1)
-            self.pars.dur_secondary.random_state.step(sim.ti+1)
+        # Hack to reset the MultiRNGs in set_secondary_prognoses so that they can be called again in this timestep. TODO: Refactor
+        self.pars.p_latent_temp.jump(sim.ti+1)
+        self.pars.dur_secondary.jump(sim.ti+1)
 
         # Secondary reactivation from latent
         secondary_from_latent = self.latent_temp & (self.ti_secondary <= sim.ti)
@@ -264,7 +263,7 @@ def set_congenital(self, sim, target_uids, source_uids=None):
 
                 # Birth outcomes must be modified to add probability of susceptible birth
                 birth_outcomes = self.pars.birth_outcomes[state]
-                assigned_outcomes = birth_outcomes.rvs(uids)-uids  # WHY??
+                assigned_outcomes = birth_outcomes.rvs(len(uids))-uids  # WHY?? # TODO: TEMP: NOT CRN SAFE
                 time_to_birth = -sim.people.age
 
                 # Schedule events
@@ -280,40 +279,36 @@ def set_congenital(self, sim, target_uids, source_uids=None):
 
 
 # %% Syphilis-related interventions
+
+__all__ += ['syph_screening', 'syph_treatment']
+
 datafiles = sc.objdict()
 for key in ['dx', 'tx', 'vx']:
     datafiles[key] = sc.thispath() / f'../data/products/syph_{key}.csv' # CK: may want to make this more robust if we keep using it
 
-__all__ += ['syph_dx', 'syph_tx', 'syph_screening', 'syph_treatment']
-
 
-def syph_dx(prod_name=None):
+def load_syph_dx():
     """
     Create default diagnostic products
     """
-    dfdx = sc.dataframe.read_csv(datafiles.dx)
+    df = sc.dataframe.read_csv(datafiles.dx)
+    hierarchy = ['positive', 'inadequate', 'negative']
     dxprods = dict(
-        rpr=ss.dx(dfdx[dfdx.name == 'rpr'], hierarchy=['positive', 'inadequate', 'negative']),
-        rst=ss.dx(dfdx[dfdx.name == 'rst'], hierarchy=['positive', 'inadequate', 'negative']),
+        rpr = ss.Dx(df[df.name == 'rpr'], hierarchy=hierarchy),
+        rst = ss.Dx(df[df.name == 'rst'], hierarchy=hierarchy),
     )
-    if prod_name is not None:
-        return dxprods[prod_name]
-    else:
-        return dxprods
+    return dxprods
 
 
-def syph_tx(prod_name=None):
+def load_syph_tx():
     """
     Create default treatment products
     """
-    dftx = sc.dataframe.read_csv(datafiles.tx)  # Read in dataframe with parameters
+    df = sc.dataframe.read_csv(datafiles.tx)  # Read in dataframe with parameters
     txprods = dict()
-    for name in dftx.name.unique():
-        txprods[name] = ss.tx(dftx[dftx.name == name])
-    if prod_name is not None:
-        return txprods[prod_name]
-    else:
-        return txprods
+    for name in df.name.unique():
+        txprods[name] = ss.Tx(df[df.name == name])
+    return txprods
 
 
 class syph_screening(ss.routine_screening):
@@ -324,12 +319,12 @@ def __init__(self, *args, **kwargs):
         return
 
     def _parse_product_str(self, product):
-        try:
-            product = syph_dx(prod_name=product)
-        except:
-            errormsg = f'Could not find product {product} in the standard list.'
+        products = load_syph_dx()
+        if product not in products:
+            errormsg = f'Could not find diagnostic product {product} in the standard list ({sc.strjoin(products.keys())})'
             raise ValueError(errormsg)
-        return product
+        else:
+            return products[product]
 
     def check_eligibility(self, sim):
         """
@@ -359,12 +354,12 @@ def __init__(self, *args, **kwargs):
         return
 
     def _parse_product_str(self, product):
-        try:
-            product = syph_tx(prod_name=product)
-        except:
-            errormsg = f'Could not find product {product} in the standard list.'
+        products = load_syph_tx()
+        if product not in products:
+            errormsg = f'Could not find treatment product {product} in the standard list ({sc.strjoin(products.keys())})'
             raise ValueError(errormsg)
-        return product
+        else:
+            return products[product]
 
     def initialize(self, sim):
         super().initialize(sim)
diff --git a/starsim/distributions.py b/starsim/distributions.py
index e53f871a..83109bd1 100644
--- a/starsim/distributions.py
+++ b/starsim/distributions.py
@@ -1,410 +1,687 @@
-""" 
-Distribution support extending scipy with two key functionalities:
-1. Callable parameters
-2. Ability to use MultiRNG for common random number support
-"""
-
-from copy import deepcopy
-import numpy as np
-import starsim as ss
-from starsim.random import SingleRNG, MultiRNG
-from starsim import options
-from scipy.stats import (bernoulli, expon, lognorm, norm, poisson, randint, rv_discrete, 
-                         uniform, rv_histogram, weibull_min)
-from scipy.stats._discrete_distns import bernoulli_gen, poisson_gen # TODO: can we remove this?
-
-
-__all__ = ['ScipyDistribution', 'ScipyHistogram']
-__all__ += ['bernoulli', 'expon', 'lognorm', 'norm', 'poisson', 'randint', 'rv_discrete', 
-            'uniform', 'weibull_min'] # Add common distributions so they can be imported directly
-
-
-class ScipyDistribution():
-    def __init__(self, gen, rng=None):
-        self._gen = gen
-        class starsim_gen(type(gen.dist)):
-            def __init__(self, *args, **kwargs):
-                super().__init__(*args, **kwargs)
-                self.sim = None
-                return
-
-            def initialize(self, sim, context):
-                self.sim = sim
-                self.context = context
-                return
-
-            def rvs_crn(self, *args, **kwargs):
-                # Bernoulli/binomial and Poisson workarounds
-                # p=0 in binomial (bernoulli) and lam=0 in Poisson do not use random numbers, which throws off CRN
-                # Solution 1: Replace 0 with eps
-                # Solution 2: Replace p (lam) where 0 with any non-zero value (e.g. eps), then revert vals back to 0
-                # --> Solution here is BOTH 1 and 2
-                if options.multirng and isinstance(self, bernoulli_gen) and ('p' in kwargs) and not np.isscalar(kwargs['p']):
-                    inds = np.where(kwargs['p']==0)[0]
-                    kwargs['p'][inds] = np.finfo(float).eps
-
-                    # Actually sample the random values
-                    vals = super().rvs(*args, **kwargs)
-
-                    # Complete the Bernoulli/binomial
-                    vals[inds] = 0
-                elif options.multirng and isinstance(self, poisson_gen) and ('lam' in kwargs) and not np.isscalar(kwargs['lam']):
-                    inds = np.where(kwargs['lam']==0)[0]
-                    kwargs['lam'][inds] = np.finfo(float).eps
-
-                    # Actually sample the random values
-                    vals = super().rvs(*args, **kwargs)
-
-                    # Complete the Poisson fix
-                    vals[inds] = 0
-                else:
-                    # Straightforward call to rvs
-                    vals = super().rvs(*args, **kwargs)
-
-                return vals
-
-            def rvs(self, *args, **kwargs):
-                """
-                Return a specified number of samples from the distribution
-                """
-
-                size = kwargs['size']
-                slots = None
-                repeat_slot_flag = False
-                self.repeat_slot_handling = {}
-                # Work out how many samples to draw. If sampling by UID, this depends on the slots assigned to agents.
-                if np.isscalar(size):
-                    if type(size) not in [int, np.int64, np.int32]: # CK: TODO: need to refactor
-                        raise Exception('Input "size" must be an integer')
-                    if size < 0:
-                        raise Exception('Input "size" cannot be negative')
-                    elif size == 0:
-                        return np.array([], dtype=int) # int dtype allows use as index, e.g. when filtering
-                    else:
-                        n_samples = size
-                elif len(size) == 0:
-                    return np.array([], dtype=int)
-                elif size.dtype == bool:
-                    n_samples = len(size) if options.multirng else size.sum()
-                elif size.dtype in [int, np.int64, np.int32]: # CK: TODO: need to refactor
-                    if not options.multirng:
-                        n_samples = len(size)
-                    else:
-                        try:
-                            slots = self.random_state.slots[size]
-                            max_slot = slots.max()
-                        except AttributeError as e:
-                            if not isinstance(self.random_state, MultiRNG):
-                                raise Exception('With options.multirng and passing agent UIDs to a distribution, the random_state of the distribution must be a MultiRNG.')
-                            else:
-                                if not self.random_state.initialized:
-                                    raise Exception('The MultiRNG instance must be initialized before use.')
-                            raise e
-
-                        if max_slot == ss.INT_NAN:
-                            raise Exception('Attempted to sample from an INT_NAN slot')
-                        n_samples = max_slot + 1
-                else:
-                    raise Exception("Unrecognized input type")
-
-                # Now handle distribution arguments
-                for pname in [p.name for p in self._param_info()]:
-                    if pname in kwargs and callable(kwargs[pname]):
-                        kwargs[pname] = kwargs[pname](self.context, self.sim, size)
-
-                    # Now do slotting if MultiRNG
-                    if options.multirng and (pname in kwargs) and (not np.isscalar(kwargs[pname])) and (len(kwargs[pname]) != n_samples):
-                        # Fill in the blank. The number of UIDs provided is
-                        # hopefully consistent with the length of pars
-                        # provided, but we need to expand out the pars to be
-                        # n_samples in length.
-                        if len(kwargs[pname]) not in [len(size), sum(size)]: # Could handle uid and bool separately? len(size) for uid and sum(size) for bool
-                            raise Exception('When providing an array of parameters, the length of the parameters must match the number of agents for the selected size (uids).')
-                        pars_slots = np.full(n_samples, fill_value=1, dtype=kwargs[pname].dtype) # self.fill_value
-                        if slots is not None:
-                            if len(slots) != len(np.unique(slots)):
-                                # Tricky - repeated slots!
-                                if not repeat_slot_flag:
-                                    repeat_slot_u, repeat_slot_ind, inv, cnt = np.unique(slots, return_index=True, return_inverse=True, return_counts=True)
-                                self.repeat_slot_handling[pname] = kwargs[pname].__array__()  # Save full pars for handling later. Use .__array__() here to provide seamless interoperability with States, UIDArrays, and np.ndarrays
-                                pars_slots[repeat_slot_u] = self.repeat_slot_handling[pname][repeat_slot_ind] # Take first instance of each
-                                repeat_slot_flag = True
-                            else:
-                                pars_slots[slots] = kwargs[pname]
-                        else:
-                            pars_slots[size] = kwargs[pname]
-                        kwargs[pname] = pars_slots
-
-                kwargs['size'] = n_samples
-                
-                # If multirng, make sure the generator is ready to avoid multiple calls without jumping in between
-                if options.multirng and not self.random_state.ready:
-                    raise ss.NotReadyException(self.random_state.name)
-
-                # Get random vals, accounting for inconsistencies in binomial draws
-                vals = self.rvs_crn(*args, **kwargs)
-
-                # Again if multirng, mark the generator as not ready (needs to be jumped)
-                if options.multirng:
-                    self.random_state.ready = False
-
-                if repeat_slot_flag:
-                    # Handle repeated slots
-                    repeat_slot_vals = np.full(len(slots), np.nan)
-                    repeat_slot_vals[repeat_slot_ind] = vals[repeat_slot_u] # Store results
-                    todo_inds = np.where(np.isnan(repeat_slot_vals))[0]
-
-                    if options.verbose > 1 and cnt.max() > 2:
-                        print(f'MultiRNG slots are repeated up to {cnt.max()} times.')
-
-                    #repeat_degree = repeat_slot_cnt.max()
-                    while len(todo_inds):
-                        repeat_slot_u, repeat_slot_ind, inv, cnt = np.unique(slots.values[todo_inds], return_index=True, return_inverse=True, return_counts=True)
-                        cur_inds = todo_inds[repeat_slot_ind] # Absolute positions being filled this pass
-
-                        # Reset RNG, note that ti=0 on initialization and ti+1
-                        # there after, including ti=0. Assuming that repeat
-                        # slots are not encountered during sim initialization.
-                        if self.sim is not None:
-                            self.random_state.step(self.sim.ti+1) 
-                        else:
-                            # Likely from a test? Just reset the random state.
-                            self.random_state.reset()
-
-                        for pname in [p.name for p in self._param_info()]:
-                            if pname in self.repeat_slot_handling:
-                                kwargs_pname = self.repeat_slot_handling[pname][cur_inds]
-                                pars_slots = np.full(n_samples, fill_value=1, dtype=kwargs_pname.dtype) # self.fill_value
-                                pars_slots[repeat_slot_u] = kwargs_pname # Take first instance of each
-                                kwargs[pname] = pars_slots
-
-                        vals = self.rvs_crn(*args, **kwargs)
-                        repeat_slot_vals[cur_inds] = vals[repeat_slot_u]
-                        todo_inds = np.where(np.isnan(repeat_slot_vals))[0]
-
-                    vals = repeat_slot_vals
-                
-                if isinstance(self, bernoulli_gen):
-                    vals = vals.astype(bool)
-
-                # _select:
-                if not options.multirng or repeat_slot_flag:
-                    return vals
-
-                if np.isscalar(size):
-                    return vals
-                elif size.dtype == bool:
-                    return vals[size]
-                else:
-                    return vals[slots] # slots defined above
-
-        rng = self.set_rng(rng, gen)
-        self.gen = starsim_gen(name=gen.dist.name, seed=rng)(**gen.kwds)
-        return
-
-    @staticmethod
-    def set_rng(rng, gen):
-        # Handle random generators
-        ret = gen.random_state # Default
-        if gen.random_state == np.random.mtrand._rand and rng:
-            if options.rng=='multi':
-                # MultiRNG, rng not none, and the current "random_state" is the
-                # numpy global singleton... so let's override
-                if isinstance(rng, str):
-                    ret = MultiRNG(rng) # Crate a new generator with the user-provided string
-                elif isinstance(rng, np.random.Generator):
-                    ret = rng
-                else:
-                    raise Exception(f'The rng must be a string or a np.random.Generator instead of {type(rng)}')
-            elif options.rng=='single':
-                # SingleRNG, rng not none, and the current "random_state" is the
-                # numpy global singleton... so let's override
-                if isinstance(rng, str):
-                    ret = SingleRNG(rng) # Crate a new generator with the user-provided string
-                elif isinstance(rng, np.random.Generator):
-                    ret = rng
-                else:
-                    raise Exception(f'The rng must be a string or a np.random.Generator instead of {type(rng)}')
-        return ret
-
-    def initialize(self, sim, context):
-        # Passing sim and context here allow callables to receive "self" and sim pointers
-        self.gen.dist.initialize(sim, context)
-        if isinstance(self.rng, (SingleRNG, MultiRNG)):
-            self.rng.initialize(sim.rng_container, sim.people.slot)
-        return
-
-    @property
-    def rng(self):
-        return self.dist.random_state
-
-    def __copy__(self):
-        cls = self.__class__
-        result = cls.__new__(cls)
-        result.__dict__.update(self.__dict__)
-        return result
-
-    def __deepcopy__(self, memo):
-        cls = self.__class__
-        result = cls.__new__(cls)
-        memo[id(self)] = result
-        for k, v in self.__dict__.items():
-            setattr(result, k, deepcopy(v, memo))
-
-        if self.gen.random_state == np.random.mtrand._rand:
-            # The gen is using the centralized numpy random number generator
-            # If we copy over the state, we'll get _separate_ random number generators
-            # for each distribution which do not draw from the _centralized_ generator
-            # in the new instance.
-            # Not clear what to do, I suppose keep as centralized?
-            result.gen.random_state = np.random.mtrand._rand
-
-        return result
-
-    def __getstate__(self):
-        dct = self.__dict__.copy()
-        dct.pop('gen')
-        return dct
-
-    def __setstate__(self, state):
-        self.__init__(state['_gen'])
-        return
-    
-    def __getattr__(self, attr):
-        # Returns wrapped generator.(attr) if not a property
-        if attr in ['__await__']:
-            return None
-        try:
-            return self.__getattribute__(attr)
-        except Exception:
-            try:
-                return getattr(self.gen, attr) # .dist?
-            except Exception:
-                errormsg = f'"{attr}" is not a member of this class or the underlying scipy stats class'
-                raise Exception(errormsg)
-
-    def filter(self, size, **kwargs):
-        return size[self.gen.rvs(size, **kwargs)]
-
-
-class ScipyHistogram(rv_histogram):
-    def __init__(self, *args, rng=None, **kwargs):
-        super().__init__(*args, **kwargs)
-        self.random_state = self.set_rng(rng)
-        return
-
-    def initialize(self, sim, context):
-        # Note: context not used here, but maintained for consistency with ScipyDistribution
-        # Passing sim and context here allow callables to receive "self" and sim pointers
-        if isinstance(self.random_state, (SingleRNG, MultiRNG)):
-            self.random_state.initialize(sim.rng_container, sim.people.slot)
-        return
-
-    def set_rng(self, rng):
-        # Handle random generators
-        ret = self.random_state # Default
-        if self.random_state == np.random.mtrand._rand and rng:
-            if options.rng=='multi':
-                # MultiRNG, rng not none, and the current "random_state" is the
-                # numpy global singleton... so let's override
-                if isinstance(rng, str):
-                    ret = MultiRNG(rng) # Crate a new generator with the user-provided string
-                elif isinstance(rng, np.random.Generator):
-                    ret = rng
-                else:
-                    raise Exception(f'The rng must be a string or a np.random.Generator instead of {type(rng)}')
-            elif options.rng=='single':
-                # SingleRNG, rng not none, and the current "random_state" is the
-                # numpy global singleton... so let's override
-                if isinstance(rng, str):
-                    ret = SingleRNG(rng) # Crate a new generator with the user-provided string
-                elif isinstance(rng, np.random.Generator):
-                    ret = rng
-                else:
-                    raise Exception(f'The rng must be a string or a np.random.Generator instead of {type(rng)}')
-        return ret
-
-
-    def rvs(self, *args, **kwargs):
-        """
-        Return a specified number of samples from the distribution
-        """
-
-        size = kwargs['size']
-        slots = None
-        # Work out how many samples to draw. If sampling by UID, this depends on the slots assigned to agents.
-        if np.isscalar(size):
-            if not isinstance(size, int):
-                raise Exception('Input "size" must be an integer')
-            if size < 0:
-                raise Exception('Input "size" cannot be negative')
-            elif size == 0:
-                return np.array([], dtype=int) # int dtype allows use as index, e.g. when filtering
-            else:
-                n_samples = size
-        elif len(size) == 0:
-            return np.array([], dtype=int)
-        elif size.dtype == bool:
-            n_samples = len(size) if options.multirng else size.sum()
-        elif size.dtype in [int, np.int64, np.int32]: # CK: TODO -- need to refactor
-            if not options.multirng:
-                n_samples = len(size)
-            else:
-                try:
-                    slots = self.random_state.slots[size]
-                    max_slot = slots.max()
-                except AttributeError as e:
-                    if not isinstance(self.random_state, MultiRNG):
-                        raise Exception('With options.multirng and passing agent UIDs to a distribution, the random_state of the distribution must be a MultiRNG.')
-                    else:
-                        if not self.random_state.initialized:
-                            raise Exception('The MultiRNG instance must be initialized before use.')
-                    raise e
-
-                if max_slot == ss.INT_NAN:
-                    raise Exception('Attempted to sample from an INT_NAN slot')
-                n_samples = max_slot + 1
-        else:
-            raise Exception("Unrecognized input type")
-
-        kwargs['size'] = n_samples
-        vals = super().rvs(*args, **kwargs)
-        
-        # _select:
-        if not options.multirng:
-            return vals
-
-        if np.isscalar(size):
-            return vals
-        elif size.dtype == bool:
-            return vals[size]
-        else:
-            return vals[slots] # slots defined above
-
-'''
-from scipy.stats import bernoulli
-class rate(ScipyDistribution):
-    """
-    Exponentially distributed, accounts for dt.
-    Assumes the rate is constant over each dt interval.
-    """
-    def __init__(self, p, rng=None):
-        dist = bernoulli(p=p)
-        super().__init__(dist, rng)
-        self.rate = rate
-        self.dt = None
-        return
-
-    def initialize(self, sim, rng):
-        self.dt = sim.dt
-        self.rng = self.set_rng(rng, self.gen)
-        super().initialize(sim)
-        return
-
-    def sample(self, size=None):
-        n_samples, pars = super().sample(size, rate=self.rate)
-        prob = 1 - np.exp(-pars['rate'] * self.dt)
-        vals = self.rng.random(size=n_samples)
-        vals = self._select(vals, size)
-        return vals < prob
-'''
+"""
+Define random-number-safe distributions.
+"""
+
+import numpy as np
+import sciris as sc
+import scipy.stats as sps
+import starsim as ss
+import pylab as pl
+
+__all__ = ['find_dists', 'dist_list', 'Dists', 'Dist']
+
+
+def str2int(string, modulo=10_000_000):
+    """
+    Convert a string to an int
+    
+    Cannot use Python's built-in hash() since it's randomized for strings, but
+    this is almost as fast (and 5x faster than hashlib).
+    """
+    return int.from_bytes(string.encode(), byteorder='big') % modulo
+
+
+def find_dists(obj, verbose=False):
+    """ Find all Dist objects in a parent object """
+    out = sc.objdict()
+    tree = sc.iterobj(obj) # Temporary copy of sc.iterobj(obj) until Sciris 3.1.5 is released
+    if verbose: print(f'Found {len(tree)} objects')
+    for trace,val in tree.items():
+        if isinstance(val, Dist):
+            key = str(trace)
+            out[key] = val
+            if verbose: print(f'  {key} is a dist ({len(out)})')
+    return out
+
+
+class Dists:
+    """ Class for managing a collection of Dist objects """
+
+    def __init__(self, obj=None, base_seed=None, sim=None):
+        self.obj = obj
+        self.dists = None
+        self.base_seed = base_seed
+        self.sim = sim
+        self.initialized = False
+        if self.obj is not None:
+            self.initialize()
+        return
+
+    def initialize(self, obj=None, base_seed=None, sim=None, force=True):
+        """
+        Set the base seed, find and initialize all distributions in an object
+        
+        In practice, the object is usually a Sim, but can be anything.
+        """
+        if base_seed:
+            self.base_seed = base_seed
+        sim = sim if (sim is not None) else self.sim
+        obj = obj if (obj is not None) else self.obj
+        if obj is None:
+            errormsg = 'Must supply a container that contains one or more Dist objects'
+            raise ValueError(errormsg)
+        self.dists = find_dists(obj)
+        for trace,dist in self.dists.items():
+            if not dist.initialized or force:
+                dist.initialize(trace=trace, seed=base_seed, sim=sim)
+        self.check_seeds()
+        self.initialized = True
+        return self
+    
+    def check_seeds(self):
+        """ Check that no two distributions share the same seed """
+        checked = dict()
+        for trace,dist in self.dists.items():
+            seed = dist.seed
+            if seed in checked.keys():
+                raise DistSeedRepeatError(checked[seed], dist)
+            else:
+                checked[seed] = dist
+        return
+
+    def jump(self, to=None, delta=1):
+        """ Advance all RNGs, e.g. to timestep "to", by jumping """
+        out = sc.autolist()
+        for dist in self.dists.values():
+            out += dist.jump(to=to, delta=delta)
+        return out
+
+    def reset(self):
+        """ Reset each RNG """
+        out = sc.autolist()
+        for dist in self.dists.values():
+            out += dist.reset()
+        return out
+
+
+class Dist: # TODO: figure out why subclassing sc.prettyobj breaks isinstance
+    """
+    Base class for tracking one random number generator associated with one distribution,
+    i.e. one decision per timestep.
+    
+    See ss.dist_list for a full list of supported distributions.
+    
+    Args:
+        dist (rv_generic): optional; a scipy.stats distribution (frozen or not) to get the ppf from
+        distname (str): the name for this class of distribution (e.g. "uniform")
+        name (str): the name for this particular distribution (e.g. "age_at_death")
+        seed (int): the user-chosen random seed (e.g. 3)
+        offset (int): the seed offset; will be automatically assigned (based on hashing the name) if None
+        strict (bool): if True, require initialization and invalidate after each call to rvs()
+        auto (bool): whether to auto-reset the state after each draw
+        sim (Sim): usually determined on initialization; the sim to use as input to callable parameters
+        module (Module): usually determined on initialization; the module to use as input to callable parameters
+        kwargs (dict): parameters of the distribution
+        
+    **Examples**::
+        
+        dist = ss.Dist(sps.norm, loc=3)
+        dist.rvs(10) # Return 10 normally distributed random numbers
+    """
+    def __init__(self, dist=None, distname=None, name=None, seed=None, offset=None, strict=False, auto=False, sim=None, module=None, **kwargs): # TODO: switch back to strict=True
+        self.dist = dist # The type of distribution
+        self.distname = distname
+        self.name = name
+        self.pars = sc.dictobj(kwargs) # The user-defined kwargs
+        self.seed = seed # Usually determined once added to the container
+        self.offset = offset
+        self.module = module
+        self.sim = sim
+        self.slots = None # Created on initialization with a sim
+        self.strict = strict
+        self.auto = auto
+        
+        # Auto-generated 
+        self.rvs_func = None # The default function to call in make_rvs() to generate the random numbers
+        self.dynamic_pars = None # Whether or not the distribution has array or callable parameters
+        self._pars = None # Validated and transformed (if necessary) parameters
+        self._n = None # Internal variable to keep track of "n" argument (usually size)
+        self._size = None # Internal variable to keep track of actual number of random variates asked for
+        self._uids = None # Internal variable to track currently-in-use UIDs
+        self._slots = None # Internal variable to track currently-in-use slots
+        
+        # History and random state
+        self.rng = None # The actual RNG generator for generating random numbers
+        self.trace = None # The path of this object within the parent
+        self.ind = 0 # The index of the RNG (usually updated on each timestep)
+        self.called = 0 # The number of times the distribution has been called
+        self.history = [] # Previous states
+        self.ready = True
+        self.initialized = False
+        if not strict: # Otherwise, wait for a sim
+            self.initialize()
+        return
+    
+    def __repr__(self):
+        """ Custom display to show state of object """
+        tracestr = '<no trace>' if self.trace is None else f"{self.trace}"
+        if self.dist is not None:
+            diststr = f'dist={self.dist}, '
+        elif self.distname is not None:
+            diststr = f'dist={self.distname}, '
+        else:
+            diststr = ''
+        string = f'ss.{self.__class__.__name__}({tracestr}, {diststr}pars={dict(self.pars)})'
+        return string
+    
+    def disp(self):
+        """ Return full display of object """
+        return sc.pr(self)
+    
+    def show_state(self):
+        """ Show the state of the object """
+        s = sc.autolist()
+        s += f'  dist = {self.dist}'
+        s += f'  pars = {self.pars}'
+        s += f' trace = {self.trace}'
+        s += f'offset = {self.offset}'
+        s += f'  seed = {self.seed}'
+        s += f'   ind = {self.ind}'
+        s += f'called = {self.called}'
+        s += f' ready = {self.ready}'
+        s += f' state = {self.state}'
+        string = sc.newlinejoin(s)
+        print(string)
+        return
+    
+    def __call__(self, n=1):
+        """ Alias to self.rvs() """
+        return self.rvs(n=n)
+    
+    def set(self, dist=None, **kwargs):
+        """ Set (change) the distribution type, or one or more parameters of the distribution """
+        if dist:
+            self.dist = dist
+            self.process_dist()
+        if kwargs:
+            self.pars.update(kwargs)
+            self.process_pars(call=False)
+        return
+
+    @property
+    def bitgen(self):
+        try:
+            return self.rng.bit_generator
+        except:
+            return None
+    
+    @property
+    def state(self):
+        """ Get the current state """
+        try:
+            return self.bitgen.state
+        except:
+            return None
+    
+    def get_state(self):
+        """ Return a copy of the state """
+        return self.state.copy()
+    
+    def make_history(self, reset=False):
+        """ Store the current state in history """
+        if reset:
+            self.history = []
+        self.history.append(self.get_state()) # Store the initial state
+        return
+
+    def reset(self, state=0):
+        """ Restore state: use 0 for initial, -1 for most recent """
+        if not isinstance(state, dict):
+            state = self.history[state]
+        self.rng.bit_generator.state = state.copy()
+        self.ready = True
+        return self.state
+
+    def jump(self, to=None, delta=1):
+        """ Advance the RNG, e.g. to timestep "to", by jumping """
+        jumps = to if (to is not None) else self.ind + delta
+        self.ind = jumps
+        self.reset() # First reset back to the initial state (used in case of different numbers of calls)
+        if jumps: # Seems to randomize state if jumps=0
+            self.bitgen.state = self.bitgen.jumped(jumps=jumps).state # Now take "jumps" number of jumps
+        return self.state
+    
+    def initialize(self, trace=None, seed=None, module=None, sim=None, slots=None):
+        """ Calculate the starting seed and create the RNG """
+        
+        # Calculate the offset (starting seed)
+        self.process_seed(trace, seed)
+        
+        # Create the actual RNG
+        self.rng = np.random.default_rng(seed=self.seed)
+        self.make_history(reset=True)
+        
+        # Handle the sim, module, and slots
+        self.sim = sim if (sim is not None) else self.sim
+        self.module = module if (module is not None) else self.module
+        if slots is None and self.sim is not None:
+            try:
+                slots = self.sim.people.slot
+            except Exception as E:
+                warnmsg = f'Could not extract slots from sim object, is this an error?\n{E}'
+                ss.warn(warnmsg)
+        if slots is not None:
+            self.slots = slots
+            
+        # Initialize the distribution and finalize
+        self.process_dist()
+        self.process_pars(call=False)
+        self.ready = True
+        self.initialized = True
+        return self
+    
+    def process_seed(self, trace=None, seed=None):
+        """ Obtain the seed offset by hashing the path to this distribution; called automatically """
+        unique_name = trace or self.trace or self.name
+        if unique_name:
+            if not self.name:
+                self.name = unique_name
+            self.trace = unique_name
+            self.offset = str2int(unique_name) # Key step: hash the path to the distribution
+        else:
+            self.offset = self.offset or 0
+        self.seed = self.offset + (seed or self.seed or 0)
+        return
+    
+    def process_dist(self):
+        """ Ensure the distribution works """
+        
+        # Handle a SciPy distribution, if provided
+        if self.dist is not None:
+            
+            # Pull out parameters of an already-frozen distribution
+            if isinstance(self.dist, sps._distn_infrastructure.rv_frozen):
+                if not self.initialized: # Don't do this more than once
+                    self.pars = sc.dictobj(sc.mergedicts(self.pars, self.dist.kwds))
+            
+            # Convert to a frozen distribution
+            if isinstance(self.dist, sps._distn_infrastructure.rv_generic):
+                spars = self.process_pars(call=False)
+                self.dist = self.dist(**spars) 
+                
+            # Override the default random state with the correct one
+            self.dist.random_state = self.rng 
+            
+        # Set the default function for getting the rvs
+        if self.distname is not None and hasattr(self.rng, self.distname): # Don't worry if it doesn't, it's probably being manually overridden
+            self.rvs_func = getattr(self.rng, self.distname) # e.g. self.rng.uniform
+        
+        return
+    
+    def process_size(self, n=1):
+        """ Handle an input of either size or UIDs and calculate size, UIDs, and slots """
+        if np.isscalar(n) or isinstance(n, tuple):  # If passing a non-scalar size, interpret as dimension rather than UIDs iff a tuple
+            uids = None
+            slots = None
+            size = n
+        else:
+            uids = np.asarray(n)
+            if len(uids):
+                slots = self.slots[uids]
+                size = slots.max() + 1
+            else:
+                slots = np.array([])
+                size = 0
+        
+        self._n = n
+        self._size = size
+        self._uids = uids
+        self._slots = slots
+        return size, slots
+    
+    def process_pars(self, call=True):
+        """ Ensure the supplied dist and parameters are valid, and initialize them; called automatically """
+        self._pars = sc.cp(self.pars) # The actual keywords; shallow copy, modified below for special cases
+        if call:
+            self.call_pars()
+        spars = self.sync_pars()
+        return spars
+    
+    def call_pars(self):
+        """ Check if any parameters need to be called to be turned into arrays """
+        
+        # Initialize
+        size, uids = self._size, self._uids
+        if self.dynamic_pars != False: # Allow "False" to prevent ever using dynamic pars (used in ss.choice())
+            self.dynamic_pars = None
+        
+        # Check each parameter
+        for key,val in self._pars.items():
+            
+            # If the parameter is callable, then call it
+            if callable(val): 
+                size_par = uids if uids is not None else size
+                out = val(self.module, self.sim, size_par) # TODO: swap order to sim, module, size?
+                val = np.asarray(out) # Necessary since UIDArrays don't allow slicing # TODO: check if this is correct
+                self._pars[key] = val
+            
+            # If it's iterable and UIDs are provided, then we need to use array-parameter logic
+            if self.dynamic_pars is None and np.iterable(val) and uids is not None:
+                self.dynamic_pars = True
+        return
+    
+    def sync_pars(self):
+        """ Perform any necessary synchronizations or transformations on distribution parameters """
+        self.update_dist_pars()
+        return self._pars
+    
+    def update_dist_pars(self, pars=None):
+        """ Update SciPy distribution parameters """
+        if self.dist is not None:
+            pars = pars if pars is not None else self._pars
+            self.dist.kwds = pars
+        return
+    
+    def rand(self, size):
+        """ Simple way to get simple random numbers """
+        return self.rng.random(size)
+    
+    def make_rvs(self):
+        """ Return default random numbers for scalar parameters """
+        if self.rvs_func is not None:
+            rvs = self.rvs_func(**self._pars, size=self._size)
+        elif self.dist is not None:
+            rvs = self.dist.rvs(self._size)
+        else:
+            errormsg = 'Could not generate random numbers: no valid NumPy function or SciPy distribution found'
+            raise ValueError(errormsg)
+        return rvs
+    
+    def ppf(self, rands):
+        """ Return default random numbers for array parameters """
+        rvs = self.dist.ppf(rands)
+        return rvs
+    
+    def rvs(self, n=1, reset=False):
+        """
+        Get random variables
+        
+        Args:
+            n (int/tuple/arr): if an int or tuple, return this many random variables; if an array, treat as UIDs
+        """
+        # Check for readiness
+        if not self.initialized:
+            raise DistNotInitializedError(self)
+        if not self.ready and self.strict:
+            raise DistNotReadyError(self)
+        
+        # Figure out size, UIDs, and slots
+        size, slots = self.process_size(n)
+        
+        # Check if size is 0, then we can return
+        if size == 0:
+            return np.array([], dtype=int) # int dtype allows use as index, e.g. when filtering
+        
+        # Store the state
+        self.make_history() # Store the pre-call state
+        
+        # Check if any keywords are callable -- parameters shouldn't need to be reprocessed otherwise
+        if True: # self.dynamic_pars: # TODO: fix!!!!
+            self.process_pars()
+        
+        # Actually get the random numbers
+        if self.dynamic_pars:
+            rands = self.rand(size)[slots] # Get random values 
+            rvs = self.ppf(rands) # Convert to actual values via the PPF
+        else:
+            rvs = self.make_rvs() # Or, just get regular values
+            if self._slots is not None:
+                rvs = rvs[self._slots]
+        
+        # Tidy up
+        self.called += 1
+        if reset:
+            self.reset(-1)
+        elif self.auto: # TODO: check
+            self.jump()
+        elif self.strict:
+            self.ready = False
+            
+        return rvs
+
+
+    def plot_hist(self, n=1000, bins=None, fig_kw=None, hist_kw=None):
+        """ Plot the current state of the RNG as a histogram """
+        pl.figure(**sc.mergedicts(fig_kw))
+        rvs = self.rvs(n)
+        self.reset(-1) # As if nothing ever happened
+        pl.hist(rvs, bins=bins, **sc.mergedicts(hist_kw))
+        pl.title(str(self))
+        pl.xlabel('Value')
+        pl.ylabel(f'Count ({n} total)')
+        return rvs
+        
+
+#%% Specific distributions
+
+# Add common distributions so they can be imported directly; assigned to a variable since used in help messages
+dist_list = ['random', 'uniform', 'normal', 'lognorm_ex', 'lognorm_im', 'expon',
+             'poisson', 'weibull', 'delta', 'randint', 'bernoulli', 'choice']
+__all__ += dist_list
+
+
+class random(Dist):
+    """ Random distribution, values on interval (0,1) """
+    def __init__(self, **kwargs):
+        super().__init__(distname='random', **kwargs)
+        return
+    
+    def ppf(self, rands):
+        return rands
+
+
+class uniform(Dist):
+    """ Uniform distribution, values on interval (low, high) """
+    def __init__(self, low=0.0, high=1.0, **kwargs):
+        super().__init__(distname='uniform', low=low, high=high, **kwargs)
+        return
+    
+    def ppf(self, rands):
+        p = self._pars
+        rvs = rands * (p.high - p.low) + p.low
+        return rvs
+
+
+class normal(Dist):
+    """ Normal distribution, with mean=loc and stdev=scale """
+    def __init__(self, loc=0.0, scale=1.0, **kwargs):
+        super().__init__(distname='normal', dist=sps.norm, loc=loc, scale=scale, **kwargs)
+        return
+
+
+class lognorm_im(Dist):
+    """
+    Lognormal distribution, parameterized in terms of the "implicit" (normal)
+    distribution, with mean=loc and stdev=scale (see lognorm_ex for comparison).
+    
+    Note: the "loc" parameter here does *not* correspond to the mean of the resulting
+    random variates!
+    
+    **Example**::
+        
+        ss.lognorm_im(mean=2, sigma=1).rvs(1000).mean() # Should be roughly 10
+    """
+    def __init__(self, mean=0.0, sigma=1.0, **kwargs):
+        super().__init__(distname='lognormal', dist=sps.lognorm, mean=mean, sigma=sigma, **kwargs)
+        return
+    
+    def sync_pars(self):
+        """ Translate between NumPy and SciPy parameters """
+        p = self._pars
+        spars = sc.dictobj()
+        spars.s = p.sigma
+        spars.scale = np.exp(p.mean)
+        spars.loc = 0
+        self.update_dist_pars(spars)
+        return spars
+    
+
+class lognorm_ex(lognorm_im):
+    """
+    Lognormal distribution, parameterized in terms of the "explicit" (lognormal)
+    distribution, with mean=mean and stdev=stdev (see lognorm_im for comparison).
+    
+    **Example**::
+        
+        ss.lognorm_ex(mean=2, stdev=1).rvs(1000).mean() # Should be close to 2
+    """
+    def __init__(self, mean=1.0, stdev=1.0, **kwargs):
+        super().__init__(mean=mean, stdev=stdev, **kwargs)
+        return
+    
+    def convert_ex_to_im(self):
+        """
+        Lognormal distributions can be specified in terms of the mean and standard
+        deviation of the "explicit" lognormal distribution, or the "implicit" normal distribution.
+        This function converts the parameters from the lognormal distribution to the
+        parameters of the underlying (implicit) distribution, which are the form expected by NumPy's
+        and SciPy's lognorm() distributions.
+        """
+        p = self._pars
+        mean = p.pop('mean')
+        stdev = p.pop('stdev')
+        if mean <= 0:
+            errormsg = f'Cannot create a lognorm_ex distribution with mean≤0 (mean={mean}); did you mean to use lognorm_im instead?'
+            raise ValueError(errormsg)
+        std2 = stdev**2
+        mean2 = mean**2
+        sigma_im = np.sqrt(np.log(std2/mean2 + 1)) # Computes stdev for the underlying normal distribution
+        mean_im  = np.log(mean2 / np.sqrt(std2 + mean2)) # Computes the mean of the underlying normal distribution
+        p.mean = mean_im
+        p.sigma = sigma_im
+        return mean_im, sigma_im
+    
+    def sync_pars(self):
+        """ Convert from overlying to underlying parameters, then translate to SciPy """
+        self.convert_ex_to_im()
+        spars = super().sync_pars()
+        return spars
+    
+
+class expon(Dist):
+    """ Exponential distribution """
+    def __init__(self, scale=1.0, **kwargs):
+        super().__init__(distname='exponential', dist=sps.expon, scale=scale, **kwargs)
+        return
+
+
+class poisson(Dist):
+    """ Exponential distribution """
+    def __init__(self, lam=1.0, **kwargs):
+        super().__init__(distname='poisson', dist=sps.poisson, lam=lam, **kwargs)
+        return
+    
+    def sync_pars(self):
+        """ Translate between NumPy and SciPy parameters """
+        spars = dict(mu=self._pars.lam)
+        self.update_dist_pars(spars)
+        return spars
+
+
+class randint(Dist):
+    """ Random integers, values on the interval [low, high-1] (i.e. "high" is excluded) """
+    def __init__(self, low=0, high=2,  **kwargs):
+        super().__init__(distname='integers', dist=sps.randint, low=low, high=high, **kwargs)
+        return
+    
+
+class weibull(Dist):
+    """ Weibull distribution -- NB, uses SciPy rather than NumPy """
+    def __init__(self, c=1, loc=0, scale=1,  **kwargs):
+        super().__init__(distname='weibull', dist=sps.weibull_min, c=c, loc=loc, scale=scale, **kwargs)
+        return
+    
+    def make_rvs(self):
+        """ Use SciPy rather than NumPy to include the scale parameter """
+        rvs = self.dist.rvs(self._size)
+        return rvs
+
+
+class delta(Dist):
+    """ Delta distribution: equivalent to np.full() """
+    def __init__(self, v=0, **kwargs):
+        super().__init__(distname='delta', v=v, **kwargs)
+        return
+    
+    def make_rvs(self):
+        return np.full(self._size, self._pars.v)
+    
+    def ppf(self, rands): # NB: don't actually need to use random numbers here, but not worth the complexity of avoiding this
+        return np.full(rands.shape, self._pars.v)
+
+
+class bernoulli(Dist):
+    """
+    Bernoulli distribution: return True or False with the specified probability (which can be an array)
+    
+    Unlike other distributions, Bernoulli distributions have a filter() method,
+    which returns elements of the array that return True.
+    """
+    def __init__(self, p=5, **kwargs):
+        super().__init__(distname='bernoulli', p=p, **kwargs)
+        return
+    
+    def make_rvs(self):
+        rvs = self.rng.random(self._size) < self._pars.p # 3x faster than using rng.binomial(1, p, size)
+        return rvs
+    
+    def ppf(self, rands):
+        rvs = rands < self._pars.p
+        return rvs
+    
+    def filter(self, uids, both=False):
+        """ Return UIDs that correspond to True, or optionally return both True and False """
+        bools = self.rvs(uids)
+        if both:
+            return uids[bools], uids[~bools]
+        else:
+            return uids[bools]
+
+
+class choice(Dist):
+    """
+    Random choice between discrete options
+    
+    **Examples**::
+        
+        # Simulate 10 die rolls
+        ss.choice(6)(10) + 1 
+        
+        # Choose between specified options each with a specified probability (must sum to 1)
+        ss.choice(a=[30, 70], p=[0.3, 0.7])(10)
+    """
+    def __init__(self, a=2, p=None, **kwargs):
+        super().__init__(distname='choice', a=a, p=p, **kwargs)
+        self.dynamic_pars = False # Set to false since array arguments don't imply dynamic pars here
+        return
+    
+    def ppf(self, rands):
+        """ Shouldn't actually be needed since dynamic pars not supported """
+        pars = self._pars
+        if np.isscalar(pars.a):
+            pars.a = np.arange(pars.a)
+        pcum = np.cumsum(pars.p)
+        inds = np.searchsorted(pcum, rands)
+        rvs = pars.a[inds]
+        return rvs
+
+
+#%% Dist exceptions
+
+class DistNotInitializedError(RuntimeError):
+    """ Raised when Dist object is called when not initialized. """
+    def __init__(self, dist):
+        msg = f'{dist} has not been initialized; please call dist.initialize()'
+        super().__init__(msg)
+
+class DistNotReadyError(RuntimeError):
+    """ Raised when a Dist object is called without being ready. """
+    def __init__(self, dist):
+        msg = f'{dist} is not ready. This is likely caused by calling a distribution multiple times in a single step. Call dist.jump() to reset.'
+        super().__init__(msg)
+        
+class DistSeedRepeatError(RuntimeError):
+    """ Raised when a Dist object shares a seed with another """
+    def __init__(self, dist1, dist2):
+        msg = f'A common seed was found between {dist1} and {dist2}. This is likely caused by incorrect initialization of the parent Dists object.'
+        super().__init__(msg)
\ No newline at end of file
diff --git a/starsim/interventions.py b/starsim/interventions.py
index 6bd92acf..2822d94a 100644
--- a/starsim/interventions.py
+++ b/starsim/interventions.py
@@ -8,16 +8,36 @@
 
 __all__ = ['Analyzer', 'Intervention']
 
+
 class Analyzer(ss.Module):
+    """
+    Base class for analyzers. Analyzers are used to provide more detailed information 
+    about a simulation than is available by default -- for example, pulling states 
+    out of sim.people on a particular timestep before they get updated on the next step.
+    
+    To retrieve a particular analyzer from a sim, use sim.get_analyzer().
+    """
 
     def __call__(self, *args, **kwargs):
         return self.apply(*args, **kwargs)
+    
+    def initialize(self, sim):
+        return super().initialize(sim)
+    
+    def apply(self, sim):
+        pass
 
     def update_results(self, sim):
-        raise NotImplementedError
+        pass
+    
+    def finalize(self, sim):
+        return super().finalize(sim)
 
 
 class Intervention(ss.Module):
+    """
+    Base class for interventions.
+    """
 
     def __init__(self, eligibility=None, *args, **kwargs):
         super().__init__(*args, **kwargs)
@@ -26,6 +46,15 @@ def __init__(self, eligibility=None, *args, **kwargs):
 
     def __call__(self, *args, **kwargs):
         return self.apply(*args, **kwargs)
+    
+    def initialize(self, sim):
+        return super().initialize(sim)
+    
+    def apply(self, sim, *args, **kwargs):
+        raise NotImplementedError
+
+    def finalize(self, sim):
+        return super().finalize(sim)
 
     def _parse_product(self, product):
         """
@@ -53,12 +82,6 @@ def check_eligibility(self, sim):
             is_eligible = sim.people.alive  # Probably not required
         return is_eligible
 
-    def apply(self, sim, *args, **kwargs):
-        raise NotImplementedError
-
-    def finalize(self, sim):
-        super().finalize(sim)
-
 
 # %% Template classes for routine and campaign delivery
 __all__ += ['RoutineDelivery', 'CampaignDelivery']
@@ -188,7 +211,7 @@ def deliver(self, sim):
         ti = sc.findinds(self.timepoints, sim.ti)[0]
         prob = self.prob[ti]  # Get the proportion of people who will be tested this timestep
         is_eligible = self.check_eligibility(sim)  # Check eligibility
-        self.coverage_dist.kwds['p'] = prob
+        self.coverage_dist.set(p=prob)
         accept_uids = self.coverage_dist.filter(ss.true(is_eligible))
         if len(accept_uids):
             self.outcomes = self.product.administer(sim, accept_uids)  # Actually administer the diagnostic
@@ -257,10 +280,12 @@ def __init__(self, product=None, prob=None, eligibility=None,
                  years=None, start_year=None, end_year=None, **kwargs):
         BaseScreening.__init__(self, product=product, eligibility=eligibility, **kwargs)
         RoutineDelivery.__init__(self, prob=prob, start_year=start_year, end_year=end_year, years=years)
+        return
 
     def initialize(self, sim):
         RoutineDelivery.initialize(self, sim)  # Initialize this first, as it ensures that prob is interpolated properly
         BaseScreening.initialize(self, sim)  # Initialize this next
+        return
 
 
 class campaign_screening(BaseScreening, CampaignDelivery):
@@ -278,10 +303,12 @@ def __init__(self, product=None, sex=None, eligibility=None,
                  prob=None, years=None, interpolate=None, **kwargs):
         BaseScreening.__init__(self, product=product, sex=sex, eligibility=eligibility, **kwargs)
         CampaignDelivery.__init__(self, prob=prob, years=years, interpolate=interpolate)
+        return
 
     def initialize(self, sim):
         CampaignDelivery.initialize(self, sim)
         BaseScreening.initialize(self, sim)  # Initialize this next
+        return
 
 
 class routine_triage(BaseTriage, RoutineDelivery):
@@ -300,10 +327,12 @@ def __init__(self, product=None, prob=None, eligibility=None,
         BaseTriage.__init__(self, product=product, eligibility=eligibility, **kwargs)
         RoutineDelivery.__init__(self, prob=prob, start_year=start_year, end_year=end_year, years=years,
                                  annual_prob=annual_prob)
+        return
 
     def initialize(self, sim):
         RoutineDelivery.initialize(self, sim)  # Initialize this first, as it ensures that prob is interpolated properly
         BaseTriage.initialize(self, sim)  # Initialize this next
+        return
 
 
 class campaign_triage(BaseTriage, CampaignDelivery):
@@ -321,10 +350,12 @@ def __init__(self, product=None, sex=None, eligibility=None,
                  prob=None, years=None, interpolate=None, annual_prob=None, **kwargs):
         BaseTriage.__init__(self, product=product, sex=sex, eligibility=eligibility, **kwargs)
         CampaignDelivery.__init__(self, prob=prob, years=years, interpolate=interpolate, annual_prob=annual_prob)
+        return
 
     def initialize(self, sim):
         CampaignDelivery.initialize(self, sim)
         BaseTriage.initialize(self, sim)
+        return
 
 
 #%% Treatment interventions
@@ -361,7 +392,7 @@ def get_accept_inds(self, sim):
         accept_uids = np.array([], dtype=int)
         eligible_uids = self.check_eligibility(sim)  # Apply eligiblity
         if len(eligible_uids):
-            self.coverage_dist.kwds['p'] = self.prob[0]
+            self.coverage_dist.set(p=self.prob[0])
             accept_uids = self.coverage_dist.filter(eligible_uids)
         return accept_uids
 
@@ -403,6 +434,7 @@ def add_to_queue(self, sim):
         """
         accept_inds = self.get_accept_inds(sim)
         if len(accept_inds): self.queue += accept_inds.tolist()
+        return
 
     def get_candidates(self, sim):
         """
@@ -466,7 +498,7 @@ def apply(self, sim):
             ti = sc.findinds(self.timepoints, sim.ti)[0]
             prob = self.prob[ti]  # Get the proportion of people who will be tested this timestep
             is_eligible = self.check_eligibility(sim)  # Check eligibility
-            self.coverage_dist.kwds['p'] = prob
+            self.coverage_dist.set(p=prob)
             accept_uids = self.coverage_dist.filter(ss.true(is_eligible))
 
             if len(accept_uids):
@@ -491,10 +523,12 @@ def __init__(self, product=None, prob=None, eligibility=None,
 
         BaseVaccination.__init__(self, product=product, eligibility=eligibility, **kwargs)
         RoutineDelivery.__init__(self, prob=prob, start_year=start_year, end_year=end_year, years=years)
+        return
 
     def initialize(self, sim):
         RoutineDelivery.initialize(self, sim)  # Initialize this first, as it ensures that prob is interpolated properly
         BaseVaccination.initialize(self, sim)  # Initialize this next
+        return
 
 
 class campaign_vx(BaseVaccination, CampaignDelivery):
@@ -508,9 +542,11 @@ def __init__(self, product=None, prob=None, eligibility=None,
 
         BaseVaccination.__init__(self, product=product, eligibility=eligibility, **kwargs)
         CampaignDelivery.__init__(self, prob=prob, years=years, interpolate=interpolate)
+        return
 
     def initialize(self, sim):
         CampaignDelivery.initialize(self, sim) # Initialize this first, as it ensures that prob is interpolated properly
         BaseVaccination.initialize(self, sim) # Initialize this next
+        return
 
 
diff --git a/starsim/modules.py b/starsim/modules.py
index 113ca03d..6f9b8294 100644
--- a/starsim/modules.py
+++ b/starsim/modules.py
@@ -5,23 +5,30 @@
 import sciris as sc
 import starsim as ss
 from scipy.stats._distn_infrastructure import rv_frozen
-from inspect import signature, _empty
-
 
 __all__ = ['Module']
 
-class Module(sc.prettyobj):
+
+class Module:#(sc.prettyobj): # TODO: replace with sc.qprettyobj
 
     def __init__(self, pars=None, par_dists=None, name=None, label=None, requires=None, **kwargs):
         self.pars = ss.omerge(pars, kwargs)
         self.par_dists = ss.omerge(par_dists)
-        self.name = name if name else self.__class__.__name__.lower() # Default name is the class name
-        self.label = label if label else self.name
+        self.name = name if (name is not None) else self.__class__.__name__.lower() # Default name is the class name
+        self.label = label if (label is not None) else self.name
         self.requires = sc.mergelists(requires)
         self.results = ss.Results(self.name)
         self.initialized = False
         self.finalized = False
         return
+    
+    def disp(self, output=False):
+        """ Display the full object """
+        out = sc.prepr(self)
+        if not output:
+            print(out)
+        else:
+            return out
 
     def check_requires(self, sim):
         """ Check that the module's requirements (of other modules) are met """
@@ -47,49 +54,50 @@ def initialize(self, sim):
         # First, convert any scalar pars to distributions if required
         for key in self.par_dists.keys():
             par = self.pars[key]
-            if isinstance(par, rv_frozen):
+            if isinstance(par, ss.Dist):
                 continue
 
-            par_dist = self.par_dists[key]
-
-            # If it's a lognormal distribution, initialize assuming the par is the desired mean
-            if par_dist.name == 'lognorm':
-                if sc.isiterable(par):
-                    if isinstance(par, dict):
-                        mu = par['mu']
-                        stdev = par['stdev']
-                    elif isinstance(par, list):
-                        mu = par[0]
-                        stdev = par[1]
-                elif sc.isnumber(par):
-                    mu = par
-                    stdev = 1
-
-                s, scale = ss.lognorm_params(mu, stdev)  # Assume stdev of 1
-                self.pars[key] = self.par_dists[key](s=s, scale=scale)
-
-            # Otherwise, figure out the required arguments and assume the user is trying to set them
+            # Handle arguments
+            args = ()
+            kwargs = {}
+            if isinstance(par, dict):
+                kwargs = par
+            elif isinstance(par, (tuple, list)):
+                args = par
             else:
-                rqrd_args = [x for x, p in signature(par_dist._parse_args).parameters.items() if p.default == _empty]
-                if len(rqrd_args) != 0:
-                    par_dist_arg = rqrd_args[0]
-                else:
-                    par_dist_arg = 'loc'
-                self.pars[key] = self.par_dists[key](**{par_dist_arg: par})
-
-        # Initialize distributions in pars
+                args = [par]
+
+            # Make the distribution
+            par_dist = self.par_dists[key]
+            if isinstance(par_dist, str):
+                try:
+                    par_dist = getattr(ss.distributions, par_dist)
+                except Exception as E:
+                    errormsg = f'"{par_dist}" is not a valid distribution name; valid distributions are: {sc.newlinejoin(ss.dists.dist_list)}'
+                    raise ValueError(errormsg) from E
+            
+            if issubclass(par_dist, ss.Dist): # Main use case
+                par_dist = par_dist(*args, sim=sim, module=self, **kwargs)
+            
+            if not isinstance(par_dist, ss.Dist):
+                print('Warning, should probably be an ss.Dist already')
+                par_dist = ss.Dist(dist=par_dist, *args, **kwargs)
+            
+            self.pars[key] = par_dist
+
+        # Initialize distributions in pars # TODO: refactor
         for key, value in self.pars.items():
             if isinstance(value, rv_frozen):
-                self.pars[key] = ss.ScipyDistribution(value, f'{self.name}_{self.label}_{key}')
-                self.pars[key].initialize(sim, self)
+                self.pars[key] = ss.Dist(value)
 
         for key, value in self.__dict__.items():
             if isinstance(value, rv_frozen):
-                setattr(self, key, ss.ScipyDistribution(value, f'{self.name}_{self.label}_{key}'))
-                getattr(self, key).initialize(sim, self)
-
-            if isinstance(value, ss.MultiRNG):
-                value.initialize(sim.rng_container, sim.people.slot)
+                setattr(self, key, ss.Dist(value))
+        
+        # Initialize everything # TODO: shouldn't be needed, should be able to recurse more
+        for key,val in list(self.pars.items()) + list(self.__dict__.items()):
+            if isinstance(val, ss.Dist):
+                val.initialize(module=self, sim=sim) # Actually a dist
 
         # Connect the states to the sim
         # Will use random numbers, so do after distribution initialization
diff --git a/starsim/network.py b/starsim/network.py
index 60130b74..42bd2a7e 100644
--- a/starsim/network.py
+++ b/starsim/network.py
@@ -483,8 +483,8 @@ def add_pairs(self, people):
         Generate contacts
         """
 
-        if isinstance(self.pars.n_contacts, ss.ScipyDistribution):
-            number_of_contacts = self.pars.n_contacts.rvs(people.alive)  # or people.uid?
+        if isinstance(self.pars.n_contacts, ss.Dist):
+            number_of_contacts = self.pars.n_contacts.rvs(people.uid[people.alive])  # or people.uid?
         else:
             number_of_contacts = np.full(len(people), self.pars.n_contacts)
 
@@ -493,7 +493,7 @@ def add_pairs(self, people):
         p1, p2 = self.get_contacts(people.uid.__array__(), number_of_contacts)
         beta = np.ones(len(p1), dtype=ss_float_)
 
-        if isinstance(self.pars.dur, ss.ScipyDistribution):
+        if isinstance(self.pars.dur, ss.Dist):
             dur = self.pars.dur.rvs(p1)
         else:
             dur = np.full(len(p1), self.pars.dur)
@@ -522,9 +522,9 @@ def __init__(self, pars=None, par_dists=None, key_dict=None, **kwargs):
         )
 
         par_dists = ss.omergeleft(par_dists,
-            duration      = ss.lognorm,
+            duration      = ss.lognorm_ex,
             participation = ss.bernoulli,
-            debut         = ss.norm,
+            debut         = ss.normal,
             acts          = ss.poisson,
         )
 
@@ -575,7 +575,7 @@ def add_pairs(self, people, ti=None):
         beta = np.ones_like(p1)
 
         # Figure out durations and acts
-        if ss.options.multirng and (len(p1) == len(np.unique(p1))):
+        if (len(p1) == len(np.unique(p1))):
             # No duplicates and user has enabled multirng, so use slotting based on p1
             dur_vals = self.pars.duration.rvs(p1)
             act_vals = self.pars.acts.rvs(p1)
@@ -610,10 +610,10 @@ class MSMNet(SexualNetwork, DynamicNetwork):
 
     def __init__(self, pars=None, key_dict=None, **kwargs):
         pars = ss.omergeleft(pars,
-            duration_dist = ss.lognorm_mean(mean=15, stdev=15),
+            duration_dist = ss.lognorm_ex(mean=15, stdev=15),
             participation_dist = ss.bernoulli(p=0.1),  # Probability of participating in this network - can vary by individual properties (age, sex, ...) using callable parameter values
-            debut_dist = ss.norm(loc=16, scale=2),
-            acts = ss.lognorm_mean(mean=80, stdev=20),
+            debut_dist = ss.normal(loc=16, scale=2),
+            acts = ss.lognorm_ex(mean=80, stdev=20),
             rel_part_rates = 1.0,
         )
         DynamicNetwork.__init__(self, key_dict, **kwargs)
@@ -637,7 +637,7 @@ def set_network_states(self, people, upper_age=None):
         # Participation
         self.participant[people.female] = False
         pr = self.pars.rel_part_rates
-        dist = ss.bernoulli.rvs(p=pr, size=len(uids))
+        dist = ss.bernoulli(pr).rvs(size=len(uids)) # TODO: not RNG safe
         self.participant[uids] = dist
 
         # Debut
@@ -652,8 +652,8 @@ def add_pairs(self, people, ti=None):
         p2 = available_m[n_pairs:n_pairs*2]
 
         # Figure out durations
-        if ss.options.multirng and (len(p1) == len(np.unique(p1))):
-            # No duplicates and user has enabled multirng, so use slotting based on p1
+        if (len(p1) == len(np.unique(p1))):
+            # No duplicates, so use slotting based on p1
             dur = self.pars.duration.rvs(p1)
             act_vals = self.pars.acts.rvs(p1)
         else:
@@ -688,16 +688,16 @@ def __init__(self, pars=None, **kwargs):
         
         """
         pars = ss.omerge({
-            'embedding_func': ss.norm(loc=self.embedding_loc, scale=2),
+            'embedding_func': ss.normal(loc=self.embedding_loc, scale=2),
             'male_shift': 5,
         }, pars)
         super().__init__(pars, **kwargs)
         return
 
     @staticmethod
-    def embedding_loc(self, sim, uids):
+    def embedding_loc(module, sim, uids):
         loc = sim.people.age[uids].values
-        loc[sim.people.female[uids]] += self.pars.male_shift  # Shift females so they will be paired with older men
+        loc[sim.people.female[uids]] += module.pars.male_shift  # Shift females so they will be paired with older men
         return loc
 
     def add_pairs(self, people, ti=None):
@@ -1048,7 +1048,7 @@ def set_participation(self, people, upper_age=None):
         remaining_pr = max(mf_pr*len(uids)-len(bi_uids), 0)/len(mf_excl_set)
 
         # Don't love the following new syntax:
-        mf_excl_uids = mf_excl_set[ss.uniform.rvs(size=len(mf_excl_set)) < remaining_pr]
+        mf_excl_uids = mf_excl_set[ss.random().rvs(size=len(mf_excl_set)) < remaining_pr] # TODO: not RNG safe and yes ugly!
 
         mf.participant[bi_uids] = True
         mf.participant[mf_excl_uids] = True
diff --git a/starsim/parameters.py b/starsim/parameters.py
index 3980beb7..be8b6c10 100644
--- a/starsim/parameters.py
+++ b/starsim/parameters.py
@@ -41,7 +41,7 @@ def __init__(self, **kwargs):
         self.dt              = 1.0           # Timestep (in years)
         self.dt_demog        = 1.0           # Timestep for demographic updates (in years)
         self.rand_seed       = 1             # Random seed, if None, don't reset
-        self.slot_scale      = 5             # Random slots will be assigned to newborn agents between min=n_agents and max=slot_scale*n_agents. Choosing a larger value here will reduce the probability of two agents using the same slot (and hence random draws), but increase the number of random numbers that are required.
+        self.slot_scale      = 2             # Random slots will be assigned to newborn agents between min=n_agents and max=slot_scale*n_agents. Choosing a larger value here will reduce the probability of two agents using the same slot (and hence random draws), but increase the number of random numbers that are required.
         self.verbose         = ss.options.verbose # Whether or not to display information during the run -- options are 0 (silent), 0.1 (some; default), 1 (default), 2 (everything)
 
         # Plug-ins: demographics, diseases, connectors, networks, analyzers, and interventions
diff --git a/starsim/people.py b/starsim/people.py
index c57f5b73..dd5eac30 100644
--- a/starsim/people.py
+++ b/starsim/people.py
@@ -152,6 +152,7 @@ def __setstate__(self, state):
         state['_states'] =  {id(v):v for v in state['_states'].values()}
         self.__dict__ = state
 
+
 class People(BasePeople):
     """
     A class to perform all the operations on the people
@@ -207,13 +208,14 @@ def __init__(self, n_agents, age_data=None, extra_states=None, rand_seed=0):
     def get_age_dist(age_data):
         """ Return an age distribution based on provided data """
         if age_data is None:
-            dist = ss.uniform(loc=0, scale=100)  # loc and width
-            return ss.ScipyDistribution(dist, 'Age distribution')
+            dist = ss.uniform(low=0, high=100, name='Age distribution')
+            return dist
 
         if sc.checktype(age_data, pd.DataFrame):
-            bb = np.append(age_data['age'].values, age_data['age'].values[-1] + 1)
-            vv = age_data['value'].values
-            return ss.ScipyHistogram((vv, bb), density=False, rng='Age distribution')
+            age_bins = age_data['age'].values
+            age_props = age_data['value'].values
+            age_props /= age_props.sum()
+            return ss.choice(a=age_bins, p=age_props)
 
 
     def initialize(self, sim):
@@ -229,7 +231,7 @@ def initialize(self, sim):
         self.slot.initialize(sim)
         self.slot[:] = self.uid
 
-        self.age_data_dist.initialize(sim, self)
+        self.age_data_dist.initialize(module=self, sim=sim)
 
         # Initialize states
         # Age is handled separately because the default value for new agents is NaN until they are concieved/born whereas
@@ -239,7 +241,7 @@ def initialize(self, sim):
             state.initialize(sim)
 
         # Assign initial ages based on the current age distribution
-        self.age[:] = self.age_data_dist.rvs(size=self.uid)
+        self.age[:] = self.age_data_dist.rvs(self.uid)
         return
 
     def add_module(self, module, force=False):
diff --git a/starsim/products.py b/starsim/products.py
index 292b602c..a8d29607 100644
--- a/starsim/products.py
+++ b/starsim/products.py
@@ -7,7 +7,7 @@
 import numpy as np
 
 
-__all__ = ['Product', 'dx', 'tx', 'vx']
+__all__ = ['Product', 'Dx', 'Tx', 'Vx']
 
 
 class Product(ss.Module):
@@ -18,7 +18,7 @@ def administer(self, people, inds):
         raise NotImplementedError
 
 
-class dx(Product):
+class Dx(Product):
     """
     Generic class for diagnostics 
     """
@@ -36,7 +36,8 @@ def __init__(self, df, hierarchy=None, *args, **kwargs):
 
         # Create placehold for multinomial sampling
         n_results = len(self.hierarchy)
-        self.result_dist = ss.rv_discrete(values=(np.arange(n_results), 1/n_results*np.ones(n_results)))
+        self.result_dist = ss.choice(a=n_results)
+        return
 
     @property
     def default_value(self):
@@ -67,7 +68,7 @@ def administer(self, sim, uids, return_format='dict'):
                 self.result_dist.pk = probs  # Overwrite distribution probabilities
 
                 # Sort people into one of the possible result states and then update their overall results
-                this_result = self.result_dist.rvs(these_uids)-these_uids
+                this_result = self.result_dist.rvs(these_uids)-these_uids # TODO: check!
                 row_inds = results.uids.isin(these_uids)
                 results.loc[row_inds, 'result'] = np.minimum(this_result, results.loc[row_inds, 'result'])
 
@@ -79,7 +80,7 @@ def administer(self, sim, uids, return_format='dict'):
         return output
 
 
-class tx(Product):
+class Tx(Product):
     """
     Treatment products change fundamental properties about People, including their prognoses and infectiousness.
     """
@@ -113,11 +114,10 @@ def administer(self, sim, uids, return_format='dict'):
                     thisdf = self.df[df_filter]  # apply filter to get the results for this state & genotype
 
                     # Determine whether treatment is successful
-                    self.efficacy_dist.kwds['p'] = thisdf.efficacy.values[0]
+                    self.efficacy_dist.set(p=thisdf.efficacy.values[0])
 
-                    if ss.options.multirng:
                     # HACK to reset the efficacy_dist as it is called multiple times per timestep. TODO: Refactor
-                        self.efficacy_dist.random_state.step(sim.ti+1)
+                    self.efficacy_dist.jump(sim.ti+1)
                     eff_treat_inds = self.efficacy_dist.filter(these_uids)
 
                     post_tx_state_name = thisdf.post_state.values[0]
@@ -138,7 +138,7 @@ def administer(self, sim, uids, return_format='dict'):
         return output
 
 
-class vx(Product):
+class Vx(Product):
     """ Vaccine product """
     def __init__(self, diseases=None, pars=None, par_dists=None, *args, **kwargs):
         pars = ss.omerge({}, pars)
diff --git a/starsim/random.py b/starsim/random.py
deleted file mode 100644
index 12c48770..00000000
--- a/starsim/random.py
+++ /dev/null
@@ -1,279 +0,0 @@
-import hashlib
-import numpy as np
-import starsim as ss
-from starsim.states import UIDArray
-from copy import deepcopy
-
-__all__ = ['RNGContainer', 'MultiRNG', 'SingleRNG', 'RNG', 'NotReadyException']
-
-
-class RNGContainer:
-    """
-    Class for managing a collection MultiRNG random number generators
-    """
-
-    def __init__(self):
-        self._rngs = ss.ndict()
-        self.used_seed_offsets = set()
-        self.initialized = False
-        return
-
-    def initialize(self, base_seed):
-        self.base_seed = base_seed
-        self.initialized = True
-        return
-    
-    def add(self, rng, check_repeats=True):
-        """
-        Add a random number generator
-        
-        Can request an offset, will check for overlap
-        Otherwise, return value will be used as the seed offset for this rng
-        """
-        if not self.initialized:
-            raise NotInitializedException()
-
-        if rng.name in self._rngs:
-            raise RepeatNameException(rng.name)
-
-        if check_repeats:
-            if rng.seed_offset in self.used_seed_offsets:
-                raise SeedRepeatException(rng.name, rng.seed_offset)
-            self.used_seed_offsets.add(rng.seed_offset)
-
-        self._rngs.append(rng)
-
-        return self.base_seed + rng.seed_offset # Add in the base seed
-
-    def step(self, ti):
-        for rng in self._rngs.dict_values():
-            rng.step(ti)
-        return
-
-    def reset(self):
-        for rng in self._rngs.dict_values():
-            rng.reset()
-        return
-
-
-class NotInitializedException(Exception):
-    "Raised when a random number generator or a RNGContainer object is called when not initialized."
-    def __init__(self, obj_name=None):
-        if obj_name is None: 
-            msg = f'An object is being used without proper initialization.'
-        else:
-            msg = f'The object named {obj_name} is being used prior to initialization.'
-        super().__init__(msg)
-        return
-
-
-class SeedRepeatException(Exception):
-    "Raised when two random number generators have the same seed."
-    def __init__(self, rng_name, seed_offset):
-        msg = f'Requested seed offset {seed_offset} for the random number generator named {rng_name} has already been used.'
-        super().__init__(msg)
-        return
-
-
-class NotReadyException(Exception):
-    "Raised when a random generator is called without being ready."
-    def __init__(self, rng_name):
-        msg = f'The random generator named "{rng_name}" was not ready when called. This error is likely caused by calling a distribution or underlying MultiRNG generator two or more times in a single step.'
-        super().__init__(msg)
-        return
-
-
-class RepeatNameException(Exception):
-    "Raised when adding a random number generator to a RNGContainer when the rng name has already been used."
-    def __init__(self, rng_name):
-        msg = f'A random number generator with name {rng_name} has already been added.'
-        super().__init__(msg)
-        return
-
-
-def RNG(*args, **kwargs):
-    """
-    Class to choose a random number generator class.
-    
-    Returns:
-    Numpy RandomState singleton or a MultiRNG: Instance of a random number generator
-    """
-    if ss.options.rng == 'multi':
-        return MultiRNG(*args, **kwargs)
-    elif ss.options.rng == 'single':
-        return SingleRNG(*args, **kwargs)
-    
-    # Centralized
-    return np.random.mtrand._rand
-
-class SingleRNG(np.random.Generator):
-    """
-    Enable a single random number generator stream per distribution without the complexity of MultiRNG.
-    Used with options.rng = 'single'
-    """
-    
-    def __init__(self, name, seed_offset=None, **kwargs):
-        """
-        Create a random number generator
-
-        seed_offset will be automatically assigned (based on hashing the name) if None
-        
-        name: a name for this random number generator, like "coin_flip"
-        """
-
-        self.name = name
-        self.kwargs = kwargs
-
-        if seed_offset is None:
-            # Obtain the seed offset by hashing the class name. Don't use python's hash because it is randomized.
-            self.seed_offset = int(hashlib.sha256(self.name.encode('utf-8')).hexdigest(), 16) % 10**8
-        else:
-            # Use user-provided seed_offset (unlikely)
-            self.seed_offset = seed_offset
-
-        self.initialized = False
-        return
-
-    def initialize(self, container, slots=None):
-        if self.initialized:
-            return
-
-        if container is not None:
-            self.seed = container.add(self) # base_seed + seed_offset
-        else:
-            # Enable use of SingleRNG without a container
-            self.seed = self.seed_offset
-
-        if 'bit_generator' not in self.kwargs:
-            self.kwargs['bit_generator'] = np.random.PCG64DXSM(seed=self.seed)
-        super().__init__(**self.kwargs)
-
-        self.initialized = True
-        return
-
-    def reset(self):
-        pass
-
-    def step(self, ti):
-        pass
-
-    def __deepcopy__(self, memo):
-        cls = self.__class__
-        result = cls.__new__(cls)
-        memo[id(self)] = result
-        for k, v in self.__dict__.items():
-            setattr(result, k, deepcopy(v, memo))
-
-        #'bit_generator' kwarg has changed
-        super(SingleRNG, result).__init__(**result.kwargs)
-
-        # and so has the _init_state
-        result._init_state = result.bit_generator.state # Store the initial state
-
-        self.seed = None # Will be determined once added to the RNG Container
-        result.initialized = True
-        result.ready = True
-
-        return result
-
-class MultiRNG(SingleRNG):
-    """
-    Class for tracking one random number generators associated with one decision per timestep.
-
-    The main use case is to sample random numbers from various distributions
-    that are specific to each agent (per decision and timestep) so as to enable
-    variance reduction between simulations through the use of common random
-    numbers. For example, the user might create a random number generator called
-    rng and ultimately ask for randomly distributed random numbers for agents
-    with UIDs 1 and 4:
-
-        >>> import starsim as ss
-        >>> import numpy as np
-        >>> rng = ss.MultiRNG('Test') # The hashed name determines the seed offset.
-        >>> rng.initialize(container=None, slots=5) # In practice, slots will be sim.people.slots. When scalar (for testing), a state-like UIDArray will be used.
-        >>> uids = np.array([1,4])
-        >>> rng.random(uids)
-        array([0.88110549, 0.86915719])
-
-    In theory, what this is doing is drawing 5 random numbers and returning the
-    draws at positions 1 and 4.
-
-    In practice, using UIDs as "slots" (the indices into the larger draw) falls
-    apart when new agents are born.  The issue is that one simulation might have
-    more births than another, so an agent born in one simulation may not
-    get the same UID as that same agent in a comparison simulation.
-    
-    The solution applied here is for each agent to have a property called "slot"
-    that is precisely the index used when selecting from an array of random
-    numbers.  When new agents are born, the mother uses her UID to sample a
-    random integer for the newborn that is used as the "slot".  With this
-    approach, newborns will be identical between two different simulations,
-    unless an intervention mechanistically drove a change.
-
-    The slot-based approach is not without challenges.
-    - Two newborn agents may received the same "slot," and thus will receive the
-    same random draws.
-    - The chance of overlapping slots can be reduced by
-    allowing mothers to choose from a larger number of possible slots (say up
-    to one-million). However, as slots are used as indices, the number of
-    random variables drawn for each query must number the maximum slot. So if
-    one agent has a slot of 1M, then 1M random numbers will be drawn,
-    consuming more time than would be necessary if the maximum slot was
-    smaller.
-    - The maximum slot is now determined by a new configure parameter named
-    "slot_scale". A value of 5 will mean that new agents will be assigned
-    slots between 1*N and 5*N, where N is sim.pars['n_agents'].
-    """
-    
-    def __init__(self, name, seed_offset=None, **kwargs):
-        """
-        Create a random number generator
-
-        seed_offset will be automatically assigned (based on hashing the name) if None
-        
-        name: a name for this random number generator, like "coin_flip"
-        """
-
-        super().__init__(name, seed_offset, **kwargs)
-
-        if seed_offset is None:
-            # Obtain the seed offset by hashing the class name. Don't use python's hash because it is randomized.
-            self.seed_offset = int(hashlib.sha256(self.name.encode('utf-8')).hexdigest(), 16) % 10**8
-        else:
-            # Use user-provided seed_offset (unlikely)
-            self.seed_offset = seed_offset
-
-        self.seed = None # Will be determined once added to the RNG Container
-        self.ready = True
-        return
-
-    def initialize(self, container, slots):
-        super().initialize(container, slots)
-
-        if isinstance(slots, int):
-            # Handle edge case in which the user wants n sequential slots, as used in testing.
-            self.slots = UIDArray(np.arange(slots),np.arange(slots))
-        else:
-            self.slots = slots # E.g. sim.people.slots (instead of using uid as the slots directly)
-
-        self._init_state = self.bit_generator.state # Store the initial state
-
-        self.ready = True
-        return
-
-    def reset(self):
-        """ Restore initial state """
-        self.bit_generator.state = self._init_state
-        self.ready = True
-        return
-
-    def step(self, ti):
-        """ Advance to time ti step by jumping """
-        
-        # First reset back to the initial state
-        self.reset()
-
-        # Now take ti jumps
-        # jumped returns a new bit_generator, use directly instead of setting state?
-        self.bit_generator.state = self.bit_generator.jumped(jumps=ti).state
-        return
diff --git a/starsim/regression/pars_v0.3.0.json b/starsim/regression/pars_v0.3.0.json
new file mode 100644
index 00000000..294a098c
--- /dev/null
+++ b/starsim/regression/pars_v0.3.0.json
@@ -0,0 +1,24 @@
+{
+  "n_agents": 10000.0,
+  "total_pop": null,
+  "pop_scale": null,
+  "remove_dead": 10,
+  "location": null,
+  "birth_rate": null,
+  "death_rate": null,
+  "start": 2000,
+  "end": null,
+  "n_years": 20,
+  "burnin": 0,
+  "dt": 0.2,
+  "dt_demog": 1.0,
+  "rand_seed": 2,
+  "slot_scale": 2,
+  "verbose": 0,
+  "demographics": {},
+  "diseases": {},
+  "networks": {},
+  "connectors": {},
+  "interventions": [],
+  "analyzers": {}
+}
\ No newline at end of file
diff --git a/starsim/results.py b/starsim/results.py
index 8447402f..67bb956c 100644
--- a/starsim/results.py
+++ b/starsim/results.py
@@ -20,7 +20,7 @@ def __new__(cls, module=None, name=None, shape=None, dtype=None, scale=None):
         return arr
     
     def __repr__(self):
-        modulestr = f'{self.module}.' if self.module else ''
+        modulestr = f'{self.module}.' if (self.module is not None) else ''
         cls_name = self.__class__.__name__
         arrstr = super().__repr__().removeprefix(cls_name)
         out = f'{cls_name}({modulestr}{self.name}):\narray{arrstr}'
diff --git a/starsim/settings.py b/starsim/settings.py
index 46e59c9b..f224779c 100644
--- a/starsim/settings.py
+++ b/starsim/settings.py
@@ -65,9 +65,6 @@ def get_orig_options():
         optdesc = sc.objdict()  # Help for the options
         options = sc.objdict()  # The options
 
-        optdesc.rng = 'Set how random numbers are handled in Starsim with three options: 1) "centralized" uses the centralized numpy random number generator for all distributions. 2) "single" uses a separate (SingleRNG) random number generator for each distribution. 3) "multi" uses a separate (MultiRNG) random number generator for each distribution. In comparing two simulations, the "single" option may be slightly better than the "centralized" option without any real disadvantages. Only "multi" can achieve full common random number (CRN) coherence, but is more computationally expensive and only produces CRN results when using CRN-safe code. The default value is "single".'
-        options.rng = ['centralized', 'single', 'multi'][1]
-
         optdesc.verbose = 'Set default level of verbosity (i.e. logging detail): e.g., 0.1 is an update every 10 simulated timesteps.'
         options.verbose = float(os.getenv('STARSIM_VERBOSE', 0.1))
 
@@ -214,9 +211,6 @@ def set_precision(self):
             raise ValueError(errormsg)
         return
 
-    @property
-    def multirng(self):
-        return self['rng'] == 'multi'
 
 # Create the options on module load
 options = Options()
diff --git a/starsim/sim.py b/starsim/sim.py
index f4c2e66f..89ed806b 100644
--- a/starsim/sim.py
+++ b/starsim/sim.py
@@ -7,21 +7,11 @@
 import sciris as sc
 import starsim as ss
 import itertools
-import numba as nb
 import pylab as pl
 
-
 __all__ = ['Sim', 'AlreadyRunError', 'demo', 'diff_sims']
 
 
-@nb.njit(cache=True)
-def set_numba_seed(value):
-    # Needed to ensure reproducibility when using random calls in numba, e.g. RandomNetwork
-    # Note, these random numbers are not currently common-random-number safe
-    np.random.seed(value)
-    return
-
-
 class Sim(sc.prettyobj):
 
     def __init__(self, pars=None, label=None, people=None, demographics=None, diseases=None, networks=None,
@@ -59,7 +49,7 @@ def __init__(self, pars=None, label=None, people=None, demographics=None, diseas
         self.analyzers = ss.ndict(analyzers, type=ss.Analyzer)
 
         # Initialize the random number generator container
-        self.rng_container = ss.RNGContainer()
+        self.dists = ss.Dists(obj=self)
 
         return
 
@@ -102,10 +92,9 @@ def initialize(self, reset=False, **kwargs):
         self.validate_dt()
         self.init_time_vecs()  # Initialize time vecs
         ss.set_seed(self.pars.rand_seed)  # Reset the random seed before the population is created
-        set_numba_seed(self.pars.rand_seed)
 
         # Initialize the core sim components
-        self.rng_container.initialize(self.pars.rand_seed + 2)  # +2 ensures that seeds from the above population initialization and the +1-offset below are not reused within the rng_container
+        self.dists.initialize(obj=self, base_seed=self.pars.rand_seed + 2)  # +2 ensures that seeds from the above population initialization and the +1-offset below are not reused within the rngs
         self.init_people(reset=reset, **kwargs)  # Create all the people (the heaviest step)
 
         # Initialize plug-ins
@@ -117,6 +106,7 @@ def initialize(self, reset=False, **kwargs):
         self.init_analyzers()
 
         # Perform post-initialization validation
+        self.dists.initialize(obj=self, base_seed=self.pars.rand_seed + 2)  # TEMP # TODO Should not be here!!!
         self.validate_post_init()
 
         # Reset the random seed to the default run seed, so that if the simulation is run with
@@ -470,7 +460,7 @@ def step(self):
             raise AlreadyRunError('Simulation already complete (call sim.initialize() to re-run)')
 
         # Advance random number generators forward to prepare for any random number calls that may be necessary on this step
-        self.rng_container.step(self.ti + 1)  # +1 offset because ti=0 is used on initialization
+        self.dists.jump(to=self.ti+1)  # +1 offset because ti=0 is used on initialization
 
         # Clean up dead agents, if removing agents is enabled
         if self.pars.remove_dead and (self.ti % self.pars.remove_dead == 0):
@@ -1003,7 +993,7 @@ def diff_sims(sim1, sim2, skip_key_diffs=False, skip=None, full=False, output=Fa
         sim2 (sim/dict): ditto
         skip_key_diffs (bool): whether to skip keys that don't match between sims
         skip (list): a list of values to skip
-        full (bool): whether to use the full summary (else, brief)
+        full (bool): whether to print out all values (not just those that differ)
         output (bool): whether to return the output as a string (otherwise print)
         die (bool): whether to raise an exception if the sims don't match
         require_run (bool): require that the simulations have been run
@@ -1040,16 +1030,19 @@ def diff_sims(sim1, sim2, skip_key_diffs=False, skip=None, full=False, output=Fa
     # Compare values
     valmatchmsg = ''
     mismatches = {}
+    n_mismatch = 0
     skip = sc.tolist(skip)
     for key in sim2.keys():  # To ensure order
         if key in sim1_keys and key not in skip:  # If a key is missing, don't count it as a mismatch
             sim1_val = sim1[key] if key in sim1 else 'not present'
             sim2_val = sim2[key] if key in sim2 else 'not present'
-            if not np.isclose(sim1_val, sim2_val, equal_nan=True):
+            mm = not np.isclose(sim1_val, sim2_val, equal_nan=True)
+            n_mismatch += mm
+            if mm or full:
                 mismatches[key] = {'sim1': sim1_val, 'sim2': sim2_val}
 
     if len(mismatches):
-        valmatchmsg = '\nThe following values differ between the two simulations:\n'
+        valmatchmsg = '\nThe following values differ between the two simulations:\n' if not full else ''
         df = sc.dataframe.from_dict(mismatches).transpose()
         diff = []
         ratio = []
@@ -1071,12 +1064,16 @@ def diff_sims(sim1, sim2, skip_key_diffs=False, skip=None, full=False, output=Fa
                     change_char = '↑'
                 elif new < old:
                     change_char = '↓'
+                elif new == old:
+                    change_char = '='
                 else:
                     errormsg = f'Could not determine relationship between sim1={old} and sim2={new}'
                     raise ValueError(errormsg)
 
                 # Set how many repeats it should have
                 repeats = 1
+                if abs_ratio == 0:
+                    repeats = 0
                 if abs_ratio >= 1.1:
                     repeats = 2
                 if abs_ratio >= 2:
@@ -1104,10 +1101,10 @@ def diff_sims(sim1, sim2, skip_key_diffs=False, skip=None, full=False, output=Fa
     # Raise an error if mismatches were found
     mismatchmsg = keymatchmsg + valmatchmsg
     if mismatchmsg:  # pragma: no cover
-        if die:
+        if die and n_mismatch: # To catch full=True case
             raise ValueError(mismatchmsg)
         elif output:
-            return mismatchmsg
+            return df
         else:
             print(mismatchmsg)
     else:
diff --git a/starsim/states.py b/starsim/states.py
index 4d8569ab..26733fc5 100644
--- a/starsim/states.py
+++ b/starsim/states.py
@@ -9,10 +9,9 @@
 from starsim.settings import INT_NAN
 from starsim.settings import dtypes as sdt
 from numpy.lib.mixins import NDArrayOperatorsMixin  # Inherit from this to automatically gain operators like +, -, ==, <, etc.
-from scipy.stats._distn_infrastructure import rv_frozen
 
 
-__all__ = ['check_dtype', 'State', 'ArrayView']
+__all__ = ['check_dtype', 'UIDArray', 'State', 'ArrayView']
 
 
 def check_dtype(dtype, default=None):
@@ -189,10 +188,12 @@ def __getitem__(self, key):
                 raise e
 
     def __setitem__(self, key, value):
-        # nb. the use of .__array__() calls is to access the array interface and thereby treat both np.ndarray and ArrayView instances
-        # in the same way without needing an additional type check. This is also why the UIDArray.dtype property is defined. Noting
-        # that for a State, the uid_map is a dynamic view attached to the People, but after an indexing operation, it will be a bare
-        # UIDArray that has an ordinary numpy array as the uid_map
+        """
+        nb. the use of .__array__() calls is to access the array interface and thereby treat both np.ndarray and ArrayView instances
+        in the same way without needing an additional type check. This is also why the UIDArray.dtype property is defined. Noting
+        that for a State, the uid_map is a dynamic view attached to the People, but after an indexing operation, it will be a bare
+        UIDArray that has an ordinary numpy array as the uid_map
+        """
         try:
             if isinstance(key, (int, np.integer)):
                 return self.values.__setitem__(self._uid_map[key], value)
@@ -318,56 +319,58 @@ def __init__(self, dtype, default=None, coerce=True):
 
     @property
     def _s(self):
-        # Return the size of the underlying array (maximum number of agents that can be stored without reallocation)
+        """ Return the size of the underlying array (maximum number of agents that can be stored without reallocation) """
         return len(self._data)
 
     @property
     def dtype(self):
-        # The specified dtype and the underlying array dtype can be different. For instance, the user might pass in
-        # ArrayView(dtype=int) but the underlying array's dtype will be np.dtype('int32'). This distinction is important
-        # because the numpy dtype has attributes like 'kind' that the input dtype may not have. We need the ArrayView's
-        # dtype to match that of the underlying array so that it can be more seamlessly exchanged with direct numpy arrays
-        # Therefore, we retain the original dtype in ArrayView._dtype() and use
+        """
+        The specified dtype and the underlying array dtype can be different. For instance, the user might pass in
+        ArrayView(dtype=int) but the underlying array's dtype will be np.dtype('int32'). This distinction is important
+        because the numpy dtype has attributes like 'kind' that the input dtype may not have. We need the ArrayView's
+        dtype to match that of the underlying array so that it can be more seamlessly exchanged with direct numpy arrays
+        Therefore, we retain the original dtype in ArrayView._dtype() and use
+        """
         return self._data.dtype
 
     def __len__(self):
-        # Return the number of active elements
+        """ Return the number of active elements """
         return self.n
 
     def __repr__(self):
-        # Print out the numpy view directly
+        """ Print out the numpy view directly """
         return self._view.__repr__()
 
     def grow(self, n):
-        # If the total number of agents exceeds the array size, extend the underlying arrays
+        """ If the total number of agents exceeds the array size, extend the underlying arrays """
         if self.n + n > self._s:
             n_new = max(n, int(self._s / 2))  # Minimum 50% growth
             self._data = np.concatenate([self._data, np.full(n_new, dtype=self.dtype, fill_value=self.default)], axis=0)
         self.n += n  # Increase the count of the number of agents by `n` (the requested number of new agents)
         self._map_arrays()
+        return
 
     def _trim(self, inds):
-        # Keep only specified indices
+        """ Keep only specified indices """
         # Note that these are indices, not UIDs!
         n = len(inds)
         self._data[:n] = self._data[inds]
         self._data[n:self.n] = self.default
         self.n = n
         self._map_arrays()
-
+        return
 
     def __getstate__(self):
-        # When pickling, skip storing the `._view` attribute, which should be re-linked after unpickling
+        """ When pickling, skip storing the `._view` attribute, which should be re-linked after unpickling """
         return {k:getattr(self, k) for k in self.__slots__ if k != '_view'}
 
     def __setstate__(self, state):
-        # Re-map arrays after unpickling so that `.view` is a reference to the correct array in-memory
+        """ Re-map arrays after unpickling so that `.view` is a reference to the correct array in-memory """
         for k,v in state.items():
             setattr(self, k, v)
         self._map_arrays()
         return
 
-
     def _map_arrays(self):
         """
         Set main simulation attributes to be views of the underlying data
@@ -382,6 +385,13 @@ def __getitem__(self, key):
 
     def __setitem__(self, key, value):
         self._view.__setitem__(key, value)
+        
+    def __getattr__(self, attr):
+        """ Make it behave like a regular array mostly -- enables things like sum(), mean(), etc. """
+        if attr in ['__deepcopy__', '__getstate__', '__setstate__']:
+            return self.__getattribute__(attr)
+        else:
+            return getattr(self._view, attr)
 
     @property
     def __array_interface__(self):
@@ -410,7 +420,7 @@ def __init__(self, name, dtype=None, default=None, label=None, coerce=True):
             default (any): Specify default value for new agents. This can be
             - A scalar with the same dtype (or castable to the same dtype) as the State
             - A callable, with a single argument for the number of values to produce
-            - An ``ss.ScipyDistribution`` instance
+            - A ``ss.Dist`` instance
             label (str): The human-readable name for the state
             coerce (bool): Whether to ensure the the data is one of the supported data types
         """
@@ -438,7 +448,7 @@ def __repr__(self):
             return UIDArray.__repr__(self)
 
     def _new_vals(self, uids):
-        if isinstance(self.default, ss.ScipyDistribution):
+        if isinstance(self.default, ss.Dist):
             new_vals = self.default.rvs(uids)
         elif callable(self.default):
             new_vals = self.default(len(uids))
@@ -454,7 +464,7 @@ def initialize(self, sim):
         specifically, `People.initialize()` and `Module.initialize()`. Initialization of a State object
         involves two processes:
 
-        - Converting any distribution objects to a ScipyDistribution instance and linking it to RNGs stored in a `Sim`
+        - Converting any distribution objects to a Dist instance and linking it to RNGs stored in a `Sim`
         - Establishing a bidirectional connection with a `People` object for the purpose of UID indexing and resizing
 
         Since State objects can be stored in `People` or in a `Module` and the collection of all states in a `Sim` should
@@ -471,11 +481,10 @@ def initialize(self, sim):
 
         people = sim.people
         assert people.initialized, 'People must be initialized before initializing states'
-
+        
         # Connect any distributions in the default to RNGs in the Sim
-        if isinstance(self.default, rv_frozen):
-            self.default = ss.ScipyDistribution(self.default, f'{self.__class__.__name__}_{self.label}')
-            self.default.initialize(sim, self)
+        if isinstance(self.default, ss.Dist):
+            self.default.initialize(module=self, sim=sim)
 
         # Establish connection with the People object
         people.register_state(self)
diff --git a/starsim/utils.py b/starsim/utils.py
index 226b9554..faa4b977 100644
--- a/starsim/utils.py
+++ b/starsim/utils.py
@@ -77,7 +77,7 @@ def _check_key(self, key, overwrite=None):
         if overwrite is None: overwrite = self._overwrite
         if key in self:
             if not overwrite:
-                typestr = ' "{self._type"}' if self._type else ''
+                typestr = f' "{self._type}"' if self._type else ''
                 errormsg = f'Cannot add object "{key}" since already present in ndict{typestr} with keys:\n{sc.newlinejoin(self.keys())}'
                 raise ValueError(errormsg)
             else:
@@ -217,7 +217,6 @@ def all_subclasses(cls):
 
 __all__ += ['set_seed']
 
-
 def set_seed(seed=None):
     '''
     Reset the random seed -- complicated because of Numba, which requires special
@@ -247,29 +246,6 @@ def set_seed_regular(seed):
     return
 
 
-# %% Helper functions related to distributions
-__all__ += ['lognorm_params', 'lognorm_mean']
-
-
-def lognorm_params(mean, stdev):
-    """
-    Returns the shape and scale parameters for scipy's parameterization of the
-    lognormal distribution which will give the specified mean and stdev
-    """
-    s = np.sqrt(np.log(stdev ** 2 / mean ** 2 + 1))
-    mu = np.log(mean ** 2 / np.sqrt(stdev ** 2 + mean ** 2))
-    scale = np.exp(mu)
-    return s, scale
-
-
-def lognorm_mean(mean, stdev):
-    """
-    Wrapper for scipy lognorm but using mean and stdev
-    """
-    s, scale = lognorm_params(mean, stdev)
-    return ss.lognorm(s=s, scale=scale)
-
-
 # %% Simple array operations
 
 __all__ += ['true', 'false', 'defined', 'undefined']
diff --git a/starsim/version.py b/starsim/version.py
index fa8ae612..c2ed42c2 100644
--- a/starsim/version.py
+++ b/starsim/version.py
@@ -4,6 +4,6 @@
 
 __all__ = ['__version__', '__versiondate__', '__license__']
 
-__version__ = '0.2.10'
-__versiondate__ = '2024-03-18'
+__version__ = '0.3.0'
+__versiondate__ = '2024-03-30'
 __license__ = f'Starsim {__version__} ({__versiondate__}) — © 2023-2024 by IDM'


Test patch:
diff --git a/tests/baseline.json b/tests/baseline.json
index 08dcfaa6..598ba2de 100644
--- a/tests/baseline.json
+++ b/tests/baseline.json
@@ -1,17 +1,17 @@
 {
   "summary": {
-    "n_alive": 9673.058823529413,
-    "new_deaths": 6.990196078431373,
-    "cum_deaths": 707.0,
+    "n_alive": 9674.078431372549,
+    "new_deaths": 6.647058823529412,
+    "cum_deaths": 669.0,
     "pregnancy_pregnancies": 0.0,
     "pregnancy_births": 0.0,
     "pregnancy_cbr": 0.0,
     "hiv_n_on_art": 0.0,
-    "hiv_n_susceptible": 8893.872549019608,
-    "hiv_n_infected": 779.1862745098039,
-    "hiv_prevalence": 0.0806248275622018,
-    "hiv_new_infections": 15.049019607843137,
-    "hiv_cum_infections": 1532.0,
-    "hiv_new_deaths": 6.990196078431373
+    "hiv_n_susceptible": 8900.813725490196,
+    "hiv_n_infected": 773.2647058823529,
+    "hiv_prevalence": 0.07997036887697175,
+    "hiv_new_infections": 14.558823529411764,
+    "hiv_cum_infections": 1482.0,
+    "hiv_new_deaths": 6.647058823529412
   }
 }
\ No newline at end of file
diff --git a/tests/benchmark.json b/tests/benchmark.json
index b8078b49..eb4dc6c7 100644
--- a/tests/benchmark.json
+++ b/tests/benchmark.json
@@ -1,13 +1,13 @@
 {
   "time": {
-    "people": 0.003,
-    "initialize": 0.014,
-    "run": 0.155
+    "people": 0.001,
+    "initialize": 0.028,
+    "run": 0.276
   },
   "parameters": {
     "n_agents": 10000,
     "n_years": 20,
     "dt": 0.2
   },
-  "cpu_performance": 0.7936101956793515
+  "cpu_performance": 0.9093911904728096
 }
\ No newline at end of file
diff --git a/tests/conftest.py b/tests/conftest.py
deleted file mode 100644
index 161b190e..00000000
--- a/tests/conftest.py
+++ /dev/null
@@ -1,12 +0,0 @@
-# Fixtures defined in this file will be used by all tests
-
-import starsim as ss
-import pytest
-
-# Run all tests with both single and multi RNG streams
-@pytest.fixture(params=['centralized', 'single', 'multi'], autouse=True)
-def set_rng_type(request):
-    original_rng_type = ss.options.rng
-    ss.options.rng = request.param
-    yield
-    ss.options.rng = original_rng_type
\ No newline at end of file
diff --git a/tests/test_base.py b/tests/test_base.py
index a8eafc36..588cfd52 100644
--- a/tests/test_base.py
+++ b/tests/test_base.py
@@ -7,9 +7,7 @@
 import numpy as np
 import starsim as ss
 import matplotlib.pyplot as plt
-import scipy.stats as sps
 
-do_plot = True
 sc.options(interactive=False) # Assume not running interactively
 
 # %% Define the tests
@@ -34,6 +32,7 @@ def test_people():
 
 
 def test_networks():
+    sc.heading('Testing networks')
 
     # Make completely abstract layers
     n_edges = 10_000
@@ -60,7 +59,7 @@ def test_networks():
     return nw1, nw2, nw3
 
 
-def test_microsim():
+def test_microsim(do_plot=False):
     sc.heading('Test making people and providing them to a sim')
 
     # Make HIV module
@@ -79,23 +78,24 @@ def test_microsim():
     sim.initialize()
     sim.run()
 
-    plt.figure()
-    plt.plot(sim.tivec, sim.results.hiv.n_infected)
-    plt.title('HIV number of infections')
+    if do_plot:
+        plt.figure()
+        plt.plot(sim.tivec, sim.results.hiv.n_infected)
+        plt.title('HIV number of infections')
 
     return sim
 
 
 def test_ppl_construction():
 
-    def init_debut(self, sim, uids):
+    def init_debut(module, sim, uids):
         # Test setting the mean debut age by sex, 16 for men and 21 for women.
         loc = np.full(len(uids), 16)
         loc[sim.people.female[uids]] = 21
         return loc
 
     mf_pars = {
-        'debut': sps.norm(loc=init_debut, scale=2),  # Age of debut can vary by using callable parameter values
+        'debut': ss.normal(loc=init_debut, scale=2),  # Age of debut can vary by using callable parameter values
     }
     sim_pars = {'networks': [ss.MFNet(mf_pars)], 'n_agents': 100}
     gon_pars = {'beta': {'mf': [0.08, 0.04]}, 'p_death': 0.2}
@@ -114,6 +114,7 @@ def init_debut(self, sim, uids):
 
 # %% Run as a script
 if __name__ == '__main__':
+    do_plot = True
     sc.options(interactive=do_plot)
 
     # Start timing
@@ -122,7 +123,7 @@ def init_debut(self, sim, uids):
     # Run tests
     ppl = test_people()
     nw1, nw2, nw3 = test_networks()
-    sim1 = test_microsim()
+    sim1 = test_microsim(do_plot)
     sim2 = test_ppl_construction()
 
     sc.toc(T)
diff --git a/tests/test_baselines.py b/tests/test_baselines.py
index 95a9c088..3ad89e3e 100644
--- a/tests/test_baselines.py
+++ b/tests/test_baselines.py
@@ -6,7 +6,6 @@
 import numpy as np
 import sciris as sc
 import starsim as ss
-import pytest
 
 do_plot = True
 do_save = False
@@ -80,10 +79,6 @@ def save_baseline():
 def test_baseline():
     ''' Compare the current default sim against the saved baseline '''
     
-    # Do not run with multi-RNG
-    if ss.options.rng != 'single':
-        pytest.skip('This test is only for the single rng.')
-
     # Load existing baseline
     baseline = sc.loadjson(baseline_filename)
     old = baseline['summary']
@@ -93,7 +88,7 @@ def test_baseline():
     new.run()
 
     # Compute the comparison
-    ss.diff_sims(old, new, full=True, die=True)
+    ss.diff_sims(old, new, die=True)
 
     return new
 
@@ -101,10 +96,6 @@ def test_baseline():
 def test_benchmark(do_save=do_save, repeats=1, verbose=True):
     ''' Compare benchmark performance '''
     
-    # Do not run with multi-RNG
-    if ss.options.multirng:
-        return
-
     if verbose: print('Running benchmark...')
     try:
         previous = sc.loadjson(benchmark_filename)
diff --git a/tests/test_dcp.py b/tests/test_dcp.py
index c0c1962d..a3bbbf37 100644
--- a/tests/test_dcp.py
+++ b/tests/test_dcp.py
@@ -8,6 +8,8 @@
 n_agents = 250
 
 
+
+
 def test_dcp():
     s1 = ss.Sim(pars=dict(diseases='sir', networks='embedding'), n_agents=n_agents)
     s1.initialize()
@@ -18,7 +20,8 @@ def test_dcp():
     s2.run()
     s1.plot()
     s2.plot()
-
+    
+    ss.diff_sims(s1, s2, full=True)
     assert np.allclose(s1.summary[:], s2.summary[:], rtol=0, atol=0, equal_nan=True)
 
     return s1, s2
diff --git a/tests/test_demographics.py b/tests/test_demographics.py
index 4589f849..8af28891 100644
--- a/tests/test_demographics.py
+++ b/tests/test_demographics.py
@@ -141,7 +141,8 @@ def test_module_adding():
     deaths = ss.Deaths(pars={'death_rate': 10})
     demographics = [births, deaths]
     with pytest.raises(Exception): # CK: should be ValueError, but that fails for now, and this is OK
-        sim = ss.Sim(n_agents=1e3, demographics=demographics, birth_rate=10, death_rate=10).run()
+        ss.Sim(n_agents=1e3, demographics=demographics, birth_rate=10, death_rate=10).run()
+    return demographics
 
 
 if __name__ == '__main__':
diff --git a/tests/test_diseases.py b/tests/test_diseases.py
index 63aee654..dc8bf92c 100644
--- a/tests/test_diseases.py
+++ b/tests/test_diseases.py
@@ -4,10 +4,8 @@
 
 # %% Imports and settings
 import starsim as ss
-import matplotlib.pyplot as plt
-import scipy.stats as sps
 import sciris as sc
-
+import matplotlib.pyplot as plt
 
 test_run = True
 n_agents = [10_000, 2_000][test_run]
@@ -18,12 +16,12 @@
 def test_sir():
     ppl = ss.People(n_agents)
     network_pars = {
-        'n_contacts': sps.poisson(mu=4), # Contacts Poisson distributed with a mean of 4
+        'n_contacts': ss.poisson(4), # Contacts Poisson distributed with a mean of 4
     }
     networks = ss.RandomNet(pars=network_pars)
 
     sir_pars = {
-        'dur_inf': sps.norm(loc=10),  # Override the default distribution
+        'dur_inf': ss.normal(loc=10),  # Override the default distribution
     }
     sir = ss.SIR(sir_pars)
 
@@ -31,10 +29,10 @@ def test_sir():
     sir.pars.beta = {'random': 0.1}
 
     # You can also change the parameters of the default lognormal distribution directly
-    # sir.pars['dur_inf'].kwds['loc'] = 5
+    sir.pars.dur_inf.set(loc=5)
 
     # Or use a function, here a lambda that makes the mean for each agent equal to their age divided by 10
-    sir.pars.dur_inf.kwds['loc'] = lambda self, sim, uids: sim.people.age[uids] / 10
+    sir.pars.dur_inf.set(loc = lambda self, sim, uids: sim.people.age[uids] / 10)
 
     sim = ss.Sim(people=ppl, diseases=sir, networks=networks)
     sim.run()
@@ -103,7 +101,7 @@ def test_multidisease():
 
     sir1.pars.beta = {'randomnet': 0.1}
     sir2.pars.beta = {'randomnet': 0.2}
-    networks = ss.RandomNet(pars=dict(n_contacts=sps.poisson(mu=4)))
+    networks = ss.RandomNet(pars=dict(n_contacts=ss.poisson(4)))
 
     sim = ss.Sim(people=ppl, diseases=[sir1, sir2], networks=networks)
     sim.run()
@@ -116,4 +114,3 @@ def test_multidisease():
     sim2 = test_ncd()
     sims = test_gavi()
     sim = test_multidisease()
-    plt.show()
diff --git a/tests/test_distributions.py b/tests/test_distributions.py
index 0bb55377..d6fe2a49 100644
--- a/tests/test_distributions.py
+++ b/tests/test_distributions.py
@@ -1,189 +1,151 @@
 """
-Test Distributions from distributions.py
+Test SciPy distributions
 """
 
 # %% Imports and settings
 import numpy as np
 import sciris as sc
 import starsim as ss
-from starsim.random import RNG
 import scipy.stats as sps
-from starsim.distributions import ScipyDistribution
-from starsim.states import UIDArray
 import pytest
-import matplotlib.pyplot as plt
 
-
-do_plot = True
-sc.options(interactive=False) # Assume not running interactively
-
-
-@pytest.fixture(params=[5, 50])
-def n(request):
-    yield request.param
+n = 5
 
 
 # %% Define the tests
-def test_basic(n):
-    dist = sps.norm(loc=1, scale=1) # Make a distribution
-    if ss.options.multirng:
-        rng = ss.MultiRNG('Uniform')
-        rng.initialize(container=None, slots=n)
-        dist.random_state = rng
-    d = ScipyDistribution(dist)
 
+def test_basic():
+    """ Basic scipy.stats test """
+    sc.heading('Test basic scipy.stats usage')
+    spsdist = sps.norm(loc=1, scale=1) # Make a distribution
+    d = ss.Dist(dist=spsdist).initialize(slots=np.arange(n)) # Convert it to Starsim
     sample = d.rvs(1)  # Draw a sample
 
-    # Reset before the next call
-    if ss.options.multirng: d.random_state.reset()
-    samples = d.rvs(10) # Draw several samples
-
-    if ss.options.multirng: d.random_state.reset()
-    samples_uid = d.rvs(size=np.array([1,3,4])) # Draw three samples
-
-    return sample, samples, samples_uid
-
-
-def test_uniform_scalar(n):
-    """ Create a uniform distribution """
-    sc.heading('test_uniform_scalar: Testing uniform with scalar parameters')
-
-    rng = ss.MultiRNG('Uniform')
-    rng.initialize(container=None, slots=n)
-    dist = sps.uniform(loc=1, scale=4)
-    dist.random_state = rng
-    d = ScipyDistribution(dist)
-
-    uids = np.array([1,3])
+    # Draw some samples
+    d.reset()
+    m = 10
+    samples = d.rvs(m) # Draw several samples
+    
+    # Draw UID samples
+    d.reset()
+    uids = np.array([0,3,4])
+    samples_uid = d.rvs(uids)
+    
+    # Print and test
+    for s in [sample, samples, samples_uid]:
+        print(s)
+    assert sample == samples[0] == samples_uid[0], 'Samples should match after reset'
+    assert len(samples) == m, 'Incorrect number of samples'
+    assert len(samples_uid) == len(uids), 'Incorrect number of samples'
+    
+    return d
+
+
+def test_scalar(n=n):
+    """ Test a basic scalar distribution """
+    sc.heading('Testing basic uniform distribution with scalar parameters')
+
+    loc = 1
+    scale = 4
+    uids = np.array([1,3,5,9])
+    
+    spsdist = sps.uniform(loc=loc, scale=scale)
+    d = ss.Dist(spsdist).initialize(slots=np.arange(uids.max()+1))
+    
     draws = d.rvs(uids)
     print(f'Uniform sample for uids {uids} returned {draws}')
 
-    assert len(draws) == len(uids)
-    return draws
-
-
-def test_uniform_scalar_str(n):
-    """ Create a uniform distribution """
-    sc.heading('test_uniform_scalar_str: Testing uniform with scalar parameters')
-
-    dist = sps.uniform(loc=1, scale=4)
-    d = ScipyDistribution(dist, 'Uniform') # String here!
-    if ss.options.rng in ['single', 'multi']:
-        d.rng.initialize(container=None, slots=n) # Only really needed for testing as initializing the distribution will do something similar.
-
-    uids = np.array([1,3])
+    assert len(draws) == len(uids), 'Incorrect number of draws'
+    assert (draws.min() > loc) and (draws.max() < loc + scale), 'Values are out of bounds'
+    
+    return d
+
+
+def test_callable(n=n):
+    """ Test callable parameters """
+    sc.heading('Testing a uniform distribution with callable parameters')
+    
+    # Define a fake people object
+    np.random.seed(1) # Since not random number safe here!
+    sim = sc.prettyobj()
+    sim.n = 10
+    sim.people = sc.prettyobj()
+    sim.people.uid = np.arange(sim.n)
+    sim.people.slot = np.arange(sim.n)
+    sim.people.age = np.random.uniform(0, 90, size=sim.n)
+
+    # Define a parameter as a lambda function
+    loc = lambda module, sim, uids: sim.people.age[uids]
+    scale = 1
+    d = ss.normal(loc=loc).initialize(sim=sim)
+
+    uids = np.array([1, 3, 7, 9])
     draws = d.rvs(uids)
-    print(f'Uniform sample for uids {uids} returned {draws}')
-
-    assert len(draws) == len(uids)
-    return draws
-
-
-def test_uniform_callable(n):
-    """ Create a uniform distribution """
-    sc.heading('test_uniform_callable: Testing uniform with callable parameters')
+    print(f'Input ages were: {sim.people.age[uids]}')
+    print(f'Output samples were: {draws}')
 
-    sim = ss.Sim().initialize()
+    meandiff = np.abs(sim.people.age[uids] - draws).mean()
+    assert meandiff < scale*3
+    return d
 
-    loc = lambda self, sim, uids: sim.people.age[uids] # Low
-    scale = 1 # Width, could also be a lambda
-    dist = sps.uniform(loc=loc, scale=scale)
 
-    d = ScipyDistribution(dist, 'Uniform')
-    d.initialize(sim, context=None)
-
-    uids = np.array([1,3])
-    draws = d.rvs(uids)
-    print(sim.people.age[uids])
-    print(f'Uniform sample for uids {uids} returned {draws}')
-
-    assert len(draws) == len(uids)
-    return draws
-
-
-def test_uniform_array(n):
-    """ Create a uniform distribution """
-    sc.heading('test_uniform: Testing uniform with a array parameters')
-
-    rng = ss.MultiRNG('Uniform')
-    rng.initialize(container=None, slots=n)
+def test_array(n=n):
+    """ Test array parameters """
+    sc.heading('Testing uniform with a array parameters')
 
     uids = np.array([1, 3])
-    loc = np.array([1, 100]) # Low
-    scale = np.array([2, 25]) # Width
-
-    dist = sps.uniform(loc=loc, scale=scale)
-    dist.random_state = rng
+    low  = np.array([1, 100]) # Low
+    high = np.array([3, 125]) # High
 
-    d = ScipyDistribution(dist)
+    d = ss.uniform(low=low, high=high).initialize(slots=np.arange(uids.max()+1))
     draws = d.rvs(uids)
     print(f'Uniform sample for uids {uids} returned {draws}')
 
     assert len(draws) == len(uids)
+    for i in range(len(uids)):
+        assert low[i] < draws[i] < low[i] + high[i], 'Invalid value'
     return draws
 
 
+@pytest.mark.skip
 def test_repeat_slot():
     """ Test behavior of repeated slots """
-    sc.heading('test_repeat_slot: Test behavior of repeated slots')
+    sc.heading('Test behavior of repeated slots')
 
-    rng = ss.MultiRNG('Uniform')
-    slots = UIDArray(values=np.array([4,2,3,2,2,3]), uid=np.arange(6))
+    # Initialize parameters
+    slots = np.array([4,2,3,2,2,3])
     n = len(slots)
-    rng.initialize(container=None, slots=slots)
-
     uids = np.arange(n)
-    loc = np.arange(n) # Low
-    scale = 1 # Width
-
-    dist = sps.uniform(loc=loc, scale=scale)
-    dist.random_state = rng
+    low = np.arange(n)
+    high = low + 1
 
-    d = ScipyDistribution(dist)
+    # Draw values
+    d = ss.uniform(low=low, high=high).initialize()
     draws = d.rvs(uids)
-    print(f'Uniform sample for uids {uids} returned {draws}')
-
-    assert len(draws) == len(uids)
+    
+    # Print and test
+    print(f'Uniform sample for slots {slots} returned {draws}')
+    assert len(draws) == len(slots)
 
     unique_slots = np.unique(slots)
     for s in unique_slots:
         inds = np.where(slots==s)[0]
         frac, integ = np.modf(draws[inds])
-        assert np.allclose(integ, loc[inds]) # Integral part should match the loc
-        if ss.options.multirng:
-            # Same random numbers, so should be same fractional part
-            assert np.allclose(frac, frac[0])
+        assert np.allclose(integ, low[inds]), 'Integral part should match the low parameter'
+        assert np.allclose(frac, frac[0]), 'Same random numbers, so should be same fractional part'
     return draws
 
 
 
 # %% Run as a script
 if __name__ == '__main__':
-    sc.options(interactive=do_plot)
-
-    n = 5
-    nTrials = 3
-
-    times = {}
-    for multirng in [True, False]:
-        ss.options(multirng=multirng)
-        key = f'MultiRNG: {multirng}'
-        times[key] = []
-        sc.heading('Testing with multirng set to', multirng)
-
-        for trial in range(nTrials):
-            T = sc.tic()
-
-            # Run tests - some will only pass if multirng is True
-            test_basic(n)
-            test_uniform_scalar(n)
-            test_uniform_scalar_str(n)
-            test_uniform_callable(n)
-            test_uniform_array(n)
-
-            times[key].append(sc.toc(T, doprint=False, output=True))
+    
+    T = sc.timer()
 
-    print('Times:', times)
+    o1 = test_basic()
+    o2 = test_scalar(n)
+    o3 = test_callable(n)
+    o4 = test_array(n)
+    # o5 = test_repeat_slot() # TODO: confirm behavior and reimplement
 
-    plt.show()
-    print('Done.')
+    T.toc()
diff --git a/tests/test_dists.py b/tests/test_dists.py
new file mode 100644
index 00000000..32e044c4
--- /dev/null
+++ b/tests/test_dists.py
@@ -0,0 +1,165 @@
+"""
+Test Dist and Dists
+"""
+
+import numpy as np
+import sciris as sc
+import scipy.stats as sps
+import starsim as ss
+import pylab as pl
+import pytest
+
+n = 1_000_000 # For testing statistical properties and performance of distributions
+m = 5 # For just testing that drawing works
+sc.options(interactive=False)
+
+
+def plot_rvs(rvs, times=None, nrows=None):
+    fig = pl.figure(figsize=(12,12))
+    nrows, ncols = sc.getrowscols(len(rvs), nrows=nrows)
+    for i,name,r in rvs.enumitems():
+        pl.subplot(nrows, ncols, i+1)
+        pl.hist(r)
+        title = times[name] if times else name
+        pl.title(title)
+        sc.commaticks()
+    sc.figlayout()
+    return fig
+
+
+# %% Define the tests
+def test_dist(m=m):
+    """ Test the Dist class """
+    sc.heading('Testing the basic Dist call')
+    dist = ss.Dist(distname='random', name='test')
+    dist.initialize()
+    rvs = dist(m)
+    print(rvs)
+    assert 0 < rvs.min() < 1, 'Values should be between 0 and 1'
+    return rvs
+
+
+def test_custom_dists(n=n, do_plot=False):
+    """ Test all custom dists """
+    sc.heading('Testing all custom distributions')
+    
+    o = sc.objdict()
+    dists = sc.objdict()
+    rvs = sc.objdict()
+    times = sc.objdict()
+    for name in ss.dist_list:
+        func = getattr(ss, name)
+        dist = func(name='test')
+        dist.initialize()
+        dists[name] = dist
+        sc.tic()
+        rvs[name] = dist.rvs(n)
+        times[name] = sc.toc(name, unit='ms', output='message')
+        print(f'{name:10s}: mean = {rvs[name].mean():n}')
+    
+    if do_plot:
+        plot_rvs(rvs, times=times, nrows=5)
+    
+    o.dists = dists
+    o.rvs = rvs
+    o.times = times
+    return o
+        
+
+def test_dists(n=n, do_plot=False):
+    """ Test the Dists container """
+    sc.heading('Testing Dists container')
+    testvals = np.zeros((2,2))
+
+    # Create the objects twice
+    for i in range(2):
+        
+        # Create a complex object containing various distributions
+        obj = sc.prettyobj()
+        obj.a = sc.objdict()
+        obj.a.mylist = [ss.random(), ss.Dist(distname='uniform', low=2, high=3)]
+        obj.b = dict(d3=ss.weibull(c=2), d4=ss.delta(v=0.3))
+        dists = ss.Dists(obj)
+        
+        # Call each distribution twice
+        for j in range(2):
+            rvs = sc.objdict()
+            for key,dist in dists.dists.items():
+                rvs[str(dist)] = dist(n)
+                dist.jump() # Reset
+            testvals[i,j] = rvs[0][283] # Pick one of the random values and store it
+    
+    # Check that results are as expected
+    print(testvals)
+    assert np.all(testvals[0,:] == testvals[1,:]), 'Newly initialized objects should match'
+    assert np.all(testvals[:,0] != testvals[:,1]), 'After jumping, values should be different'
+            
+    if do_plot:
+        plot_rvs(rvs)
+    
+    o = sc.objdict()
+    o.dists = dists
+    o.rvs = rvs
+    return dists
+
+
+def test_scipy(m=m):
+    """ Test that SciPy distributions also work """
+    sc.heading('Testing SciPy distributions')
+    
+    # Make SciPy distributions in two different ways
+    dist1 = ss.Dist(dist=sps.expon, name='scipy', scale=2).initialize() # Version 1: callable
+    dist2 = ss.Dist(dist=sps.expon(scale=2), name='scipy').initialize() # Version 2: frozen
+    rvs1 = dist1(m)
+    rvs2 = dist2(m)
+    
+    # Check that they match
+    print(rvs1)
+    assert np.array_equal(rvs1, rvs2), 'Arrays should match'
+    
+    return dist1, dist2
+
+@pytest.mark.skip
+def test_exceptions(m=m):
+    """ Check that exceptions are being appropriately raised """
+    sc.heading('Testing exceptions and strict')
+    
+    # Create a strict distribution
+    dist = ss.random()
+    with pytest.raises(ss.dists.DistNotInitializedError):
+        dist(m) # Check that we can't call an uninitialized
+    
+    # Initialize and check we can't call repeatedly
+    dist.initialize()
+    rvs = dist(m)
+    with pytest.raises(ss.dists.DistNotReadyError):
+        dist(m) # Check that we can't call an already-used distribution
+    
+    # Check that we can with a non-strict Dist
+    dist2 = ss.random(strict=False)
+    rvs2 = sc.autolist()
+    for i in range(2):
+        rvs2 += dist2(m) # We should be able to call multiple times with no problem
+    
+    print(rvs)
+    print(rvs2)
+    assert np.array_equal(rvs, rvs2[0]), 'Separate dists should match'
+    assert not np.array_equal(rvs2[0], rvs2[1]), 'Multiple calls to the same dist should not match'
+    
+    return dist, dist2
+    
+
+# %% Run as a script
+if __name__ == '__main__':
+    do_plot = True
+    sc.options(interactive=do_plot)
+    
+    T = sc.timer()
+    
+    o1 = test_dist()
+    o2 = test_custom_dists(do_plot=do_plot)
+    o3 = test_dists(do_plot=do_plot)
+    o4 = test_scipy()
+    # o5 = test_exceptions() # TODO: re-enable once strict=True
+    
+    T.toc()
diff --git a/tests/test_random.py b/tests/test_random.py
index 85e2b988..d81b2f76 100644
--- a/tests/test_random.py
+++ b/tests/test_random.py
@@ -6,122 +6,108 @@
 import numpy as np
 import sciris as sc
 import starsim as ss
-from starsim.random import RNG
-import pytest
 
+n = 5 # Default number of samples
 
-@pytest.fixture
-def rng(request, slots=5, base_seed=1, name='Test', **kwargs):
-    return make_rng(slots, base_seed, name, **kwargs)
-
-
-def make_rng(slots=5, base_seed=1, name='Test', **kwargs):
-    rng_container = ss.RNGContainer()
-    rng_container.initialize(base_seed=base_seed)
-    rng = ss.RNG(name, **kwargs)
-    if ss.options.rng in ['single', 'multi']:
-        rng.initialize(rng_container, slots=slots)
-    return rng
+def make_dist(seed=1, name='test', **kwargs):
+    """ Make a default distribution for testing """
+    dist = ss.Dist(distname='random', name=name, seed=seed, **kwargs)
+    return dist
 
 
 # %% Define the tests
-def test_random(rng, n=5):
+def test_random(n=n):
     """ Simple random draw """
-    sc.heading('test_random: Testing simple random draw from a RNG object')
-    draws = rng.random(n)
+    sc.heading('Testing simple random draw from a RNG object')
+    
+    dist = make_dist()
+    draws = dist(n)
     print(f'Sampled {n} random draws', draws)
     assert len(draws) == n
     return draws
 
 
-def test_SingleRNG_using_np_singleton(rng, n=5):
-    """ Make sure the SingleRNG is using the numpy.random singleton """
-    if ss.options.rng != 'centralized':
-        pytest.skip('This test is only for the centralized rng.')
-    sc.heading('test_singleton_rng: Testing RNG object')
-
-    np.random.seed(0)
-    draws_np = np.random.rand(n)
-    np.random.seed(0)
-    draws_SingleRNG = rng.random(n)
-    assert np.array_equal(draws_np, draws_SingleRNG)
-    return draws_np, draws_SingleRNG
-
-
-def test_reset(rng, n=5):
+def test_reset(n=n):
     """ Sample, reset, sample """
-    if ss.options.rng not in ['single', 'multi']:
-        pytest.skip('Reset is only for SingleRNG and MultiRNG.')
-    sc.heading('test_reset: Testing sample, reset, sample')
+    sc.heading('Testing sample, reset, sample')
+    
+    # Make dist and draw twice
+    dist = make_dist()
+    draws1 = dist(n)
+    dist.reset()
+    draws2 = dist(n)
+    
+    # Print and test results
+    print(f'Random sample of size {n} returned:\n{draws1}')
+    print(f'After reset, random sample of size {n} returned:\n{draws2}')
+    assert np.array_equal(draws1, draws2)
 
-    draws1 = rng.random(n)
-    print(f'Random sample of size {n} returned {draws1}')
-
-    print('Reset')
-    rng.reset()
-
-    draws2 = rng.random(n)
-    print(f'After reset, random sample of size {n} returned {draws2}')
-
-    if isinstance(rng, ss.MultiRNG):
-        assert np.array_equal(draws1, draws2)
-    else:
-        assert not np.array_equal(draws1, draws2)
     return draws1, draws2
 
 
-def test_step(rng, n=5):
-    """ Sample, step, sample """
-    if ss.options.rng not in ['single', 'multi']:
-        pytest.skip('Step only for SingleRNG and MultiRNG.')
-    sc.heading('test_step: Testing sample, step, sample')
-
-    draws1 = rng.random(n)
-    print(f'Random sample of size {n} returned {draws1}')
-
-    print('Reset')
-    rng.step(1)
-
-    draws2 = rng.random(n)
-    print(f'After reset, random sample of size {n} returned {draws2}')
-
+def test_jump(n=n):
+    """ Sample, jump, sample """
+    sc.heading('Testing sample, jump, sample')
+    
+    dist = make_dist()
+    draws1 = dist(n)
+    dist.jump()
+    draws2 = dist(n)
+    
+    print(f'Random sample of size {n} returned:\n{draws1}')
+    print(f'After jump, random sample of size {n} returned:\n{draws2}')
     assert not np.array_equal(draws1, draws2)
+    
     return draws1, draws2
 
 
-def test_seed(n=5):
+def test_seed(n=n):
     """ Changing seeds """
-    sc.heading('test_seed: Testing sample with seeds 0 and 1')
-
-    rng0 = make_rng(n, base_seed=0)
-    draws0 = rng0.uniform(size=n)
-    print(f'Random sample of size {n} for rng0 with base_seed 0 returned {draws0}')
-
-    rng1 = make_rng(n, base_seed=1)
-    draws1 = rng1.uniform(size=n)
-    print(f'Random sample of size {n} for rng1 with base_seed 1 returned {draws1}')
-
+    sc.heading('Testing sample with seeds 0 and 1')
+
+    dist0 = make_dist(seed=0)
+    dist1 = make_dist(seed=1)
+    draws0 = dist0(n)
+    draws1 = dist1(n)
+    
+    print(f'Random sample of size {n} for dist0 with seed=0 returned:\n{draws0}')
+    print(f'Random sample of size {n} for dist1 with seed=1 returned:\n{draws1}')
     assert not np.array_equal(draws0, draws1)
+    
     return draws0, draws1
 
 
+def test_rvs(n=n):
+    """ Simple sample from distribution by UID """
+    sc.heading('Testing UID sample')
+    
+    slots = np.arange(n+1)
+    dist = make_dist()
+    dist.initialize(slots=slots)
+    uids = np.arange(0, n, 2) # every other to make it interesting
+    draws = dist.rvs(uids)
+    print(f'Created seed and sampled: {draws}')
+    assert len(draws) == len(uids)
+    
+    # Draws without UIDs should match the first element only
+    dist2 = make_dist()
+    dist2.initialize(slots=slots)
+    draws2 = dist2.rvs(len(uids))
+    assert draws[0] == draws2[0]
+    assert not np.array_equal(draws, draws2)
+    
+    return draws
+
+
 # %% Run as a script
 if __name__ == '__main__':
-    # Start timing
-    T = sc.tic()
-
-    n=5
-    rng = make_rng()
 
-    for multirng in [True, False]:
-        ss.options(multirng=multirng)
-        sc.heading('Testing with multirng set to', multirng)
+    T = sc.timer()
 
-        # Run tests - some will only pass if multirng is True
-        test_random(rng, n)
-        test_reset(rng, n)
-        test_step(rng, n)
-        test_seed(n)
+    o1 = test_random(n)
+    o2 = test_reset(n)
+    o3 = test_jump(n)
+    o4 = test_seed(n)
+    o5 = test_rvs(n)
 
-    sc.toc(T)
-    print('Done.')
+    T.toc()
diff --git a/tests/test_rngcontainer.py b/tests/test_rngcontainer.py
index 20ba712c..eb7d64d6 100644
--- a/tests/test_rngcontainer.py
+++ b/tests/test_rngcontainer.py
@@ -1,199 +1,117 @@
 """
-Test the RNGContainer object from random.py
+Test the RNGs object from random.py
 """
 
 # %% Imports and settings
 import numpy as np
 import sciris as sc
 import starsim as ss
-from starsim.random import NotInitializedException, SeedRepeatException, RepeatNameException
-from starsim.random import SingleRNG, MultiRNG
-import scipy.stats as sps
-import pytest
 
+n = 5 # Default number of samples
 
-# %% Define the tests
+def make_dist(name='test', **kwargs):
+    """ Make a default Dist for testing """
+    dist = ss.random(name=name, **kwargs).initialize()
+    return dist
 
-@pytest.fixture
-def rng_container():
-    rng_container = ss.RNGContainer()
-    rng_container.initialize(base_seed=10)
-    return rng_container
-
-@pytest.fixture
-def rngs():
-    return [ss.RNG('rng0'), ss.RNG('rng1')]
-
-@pytest.fixture
-def dists():
-    names = ['rng1', 'rng2']
-    dists = []
-    for rng_name in names:
-        d = ss.ScipyDistribution(sps.uniform(), rng=rng_name)
-        dists.append(d)
+def make_dists(**kwargs):
+    """ Make a Dists object with two distributions in it """
+    distlist = [make_dist(), make_dist()]
+    dists = ss.Dists(distlist)
+    dists.initialize()
     return dists
 
 
-def test_rng(rng_container, rngs, n=5):
-    """ Simple sample from rng """
-    sc.heading('test_rng: Testing RNGContainer object')
-
-    rng0 = rngs[0]
-    if ss.options.rng in ['single', 'multi']:
-        rng0.initialize(rng_container, slots=n)
-
-    draws = rng0.random(n)
-    print(f'Created seed and sampled: {draws}')
-
-    assert len(draws) == n
-    return draws
-
-
-def test_dists(rng_container, dists, n=5):
-    """ Simple sample from distribution """
-    sc.heading('test_dists: Testing RNGContainer object')
-
-    d0 = dists[0]
-    if isinstance(d0.rng, (SingleRNG, MultiRNG)):
-        d0.rng.initialize(rng_container, slots=n)
-
-    uids = np.arange(0,n,2) # every other to make it interesting
-    draws = d0.rvs(uids)
-    print(f'Created seed and sampled: {draws}')
-
-    assert len(draws) == len(uids)
-    return draws
-
+# %% Define the tests
 
-def test_seed(rng_container, rngs, n=5):
+def test_seed():
     """ Test assignment of seeds """
-    if not ss.options.multirng:
-        pytest.skip('Seed assignment is only for multirng.')
-
-    sc.heading('test_seed: Testing assignment of seeds')
-
-    rng0, rng1 = rngs
-    rng0.initialize(rng_container, slots=n)
-    rng1.initialize(rng_container, slots=n)
+    sc.heading('Testing assignment of seeds')
+    
+    # Create and initialize two distributions
+    dists = make_dists()
+    dist0, dist1 = dists.dists.values()
 
-    print(f'Random generators rng0 and rng1 were assigned seeds {rng0.seed} and {rng1.seed}, respectively')
+    print(f'Dists dist0 and dist1 were assigned seeds {dist0.seed} and {dist1.seed}, respectively')
+    assert dist0.seed != dist1.seed
+    return dist0, dist1
 
-    assert rng1.seed != rng0.seed
-    return rng0, rng1
 
-
-def test_reset(rng_container, dists, n=5):
+def test_reset(n=n):
     """ Sample, reset, sample """
-    sc.heading('test_reset: Testing sample, reset, sample')
-
-    d0 = dists[0]
-    if ss.options.rng in ['single', 'multi']:
-        d0.rng.initialize(rng_container, slots=n)
-
-    uids = np.arange(0,n,2) # every other to make it interesting
-    s_before = d0.rvs(uids)
-    rng_container.reset() # Return to step 0
-    s_after = d0.rvs(uids)
-
-    print(f'Initial sample', s_before)
-    print('Reset')
-    print(f'After reset sample', s_after)
-
-    if ss.options.multirng:
-        assert np.array_equal(s_after, s_before)
-    else:
-        assert not np.array_equal(s_after, s_before)
-    return s_before, s_after
-
-
-def test_step(rng_container, dists, n=5):
-    """ Sample, step, sample """
-    sc.heading('test_step: Testing sample, step, sample')
-
-    d0 = dists[0]
-    if ss.options.rng in ['single', 'multi']:
-        d0.rng.initialize(rng_container, slots=n)
-
-    uids = np.arange(0,n,2) # every other to make it interesting
-    s_before = d0.rvs(uids)
-    rng_container.step(10) # 10 steps
-    s_after = d0.rvs(uids)
-
-    print(f'Initial sample', s_before)
-    print('Step')
-    print(f'Subsequent sample', s_after)
-
-    assert not np.array_equal(s_after, s_before)
-    return s_before, s_after
-
-
-def test_initialize(rngs, n=5):
-    """ Sample without initializing, should raise exception """
-    if ss.options.rng not in ['single', 'multi']:
-        pytest.skip('Initialize only for SingleRNG and MultiRNG.')
-    sc.heading('test_initialize: Testing without initializing, should raise exception.')
-    rng_container = ss.RNGContainer()
-    #rng_container.initialize(base_seed=3) # Do not initialize
-
-    with pytest.raises(NotInitializedException):
-        rngs[0].initialize(rng_container, slots=n)
-    return rngs[0]
-
-
-def test_seedrepeat(rng_container, n=5):
-    """ Two random number generators with the same seed, should raise exception """
-    sc.heading('test_seedrepeat: Testing two random number generators with the same seed, should raise exception.')
-
-    rng0 = ss.MultiRNG('rng0', seed_offset=0)
-    rng0.initialize(rng_container, slots=n)
-
-    with pytest.raises(SeedRepeatException):
-        rng1 = ss.MultiRNG('rng1', seed_offset=0)
-        rng1.initialize(rng_container, slots=n)
-    return rng0, rng1
-
-
-def test_samplingorder(rng_container, dists, n=5):
+    sc.heading('Testing sample, reset, sample')
+    dists = make_dists()
+    distlist = dists.dists.values()
+
+    # Reset via the container, but only 
+    before = sc.autolist()
+    after = sc.autolist()
+    for dist in distlist:
+        before += list(dist(n))
+    dists.reset() # Return to step 0
+    for dist in distlist:
+        after += list(dist(n))
+
+    print(f'Initial sample:\n{before}')
+    print(f'After reset sample:\n{after}')
+    assert np.array_equal(before, after)
+    
+    return before, after
+
+
+def test_jump(n=n):
+    """ Sample, jump, sample """
+    sc.heading('Testing sample, jump, sample')
+    dists = make_dists()
+    distlist = dists.dists.values()
+
+    # Jump via the contianer
+    before = sc.autolist()
+    after = sc.autolist()
+    for dist in distlist:
+        before += list(dist(n))
+    dists.jump(to=10) # Jump to 10th step
+    for dist in distlist:
+        after += list(dist(n))
+
+    print(f'Initial sample:\n{before}')
+    print(f'After jump sample:\n{after}')
+    assert not np.array_equal(before, after)
+    
+    return before, after
+
+
+def test_order(n=n):
     """ Ensure sampling from one RNG doesn't affect another """
-    sc.heading('test_samplingorder: Testing from multiple random number generators to test if sampling order matters')
-
-    uids = np.arange(0,n,2) # every other to make it interesting
-
-    d0, d1 = dists
-    if ss.options.rng in ['single', 'multi']:
-        d0.rng.initialize(rng_container, slots=n)
-        d1.rng.initialize(rng_container, slots=n)
-
-    s_before = d0.rvs(uids)
-    _ = d1.rvs(uids)
-
-    rng_container.reset()
+    sc.heading('Testing from multiple random number generators to test if sampling order matters')
+    dists = make_dists()
+    d0, d1 = dists.dists.values()
 
-    _ = d1.rvs(uids)
-    s_after = d0.rvs(uids)
+    # Sample d0, d1
+    before = d0(n)
+    _ = d1(n)
+    
+    dists.reset()
+    
+    # Sample d1, d0
+    _ = d1(n)
+    after = d0(n)
 
-    print(f'When sampling rng0 before rng1:', s_before)
-    print('Reset')
-    print(f'When sampling rng0 after rng1:', s_after)
+    print(f'When sampling rng0 before rng1: {before}')
+    print(f'When sampling rng0 after rng1: {after}')
+    assert np.array_equal(before, after)
 
-    if ss.options.multirng:
-        assert np.array_equal(s_before, s_after)
-    else:
-        assert not np.array_equal(s_before, s_after)
+    return before, after
 
-    return s_before, s_after
 
+# %% Run as a script
+if __name__ == '__main__':
 
-def test_repeatname(rng_container, n=5):
-    """ Test two random number generators with the same name """
-    if ss.options.rng not in ['single', 'multi']:
-        pytest.skip('Repeate name is only for SingleRNG and MultiRNG.')
-    sc.heading('test_repeatname: Testing if two random number generators with the same name are allowed')
+    T = sc.timer()
 
-    rng0 = ss.RNG('test')
-    rng0.initialize(rng_container, slots=n)
+    o1 = test_seed()
+    o2 = test_reset(n)
+    o3 = test_jump(n)
+    o4 = test_order(n)
 
-    rng1 = ss.RNG('test')
-    with pytest.raises(RepeatNameException):
-        rng1.initialize(rng_container, slots=n)
-    return rng0, rng1
+    T.toc()
diff --git a/tests/test_samples.py b/tests/test_samples.py
index 86ea8f6b..105662a4 100644
--- a/tests/test_samples.py
+++ b/tests/test_samples.py
@@ -1,5 +1,8 @@
+"""
+Test the Samples class
+"""
+
 import starsim as ss
-import pytest
 
 testdir = ss.root/'tests'
 tempdir = testdir/'temp'
@@ -9,9 +12,9 @@ def get_outputs(p_death):
     outputs = []
     for i in range(3):
         ppl = ss.People(1000)
-        ppl.networks = ss.ndict(ss.RandomNet(n_contacts=ss.poisson(5)))
+        network = ss.RandomNet(n_contacts=ss.poisson(5))
         sir = ss.SIR(pars={'p_death':p_death})
-        sim = ss.Sim(people=ppl, diseases=sir, rand_seed=0, n_years=5)
+        sim = ss.Sim(people=ppl, networks=network, diseases=sir, rand_seed=0, n_years=5)
         sim.run(verbose=0)
         df = sim.export_df()
         summary = {}
@@ -23,9 +26,6 @@ def get_outputs(p_death):
     return outputs
 
 def test_samples():
-    if ss.options.multirng:
-        # Not necessary to run with both single and multi RNG
-        pytest.skip('Skipping multirng mode as not necessary to test with both single and multi RNGs')
     outputs = get_outputs(0.2)
     resultsdir = tempdir/'samples_results'
     resultsdir.mkdir(exist_ok=True)
@@ -33,8 +33,6 @@ def test_samples():
     return s
 
 def test_dataset():
-    if ss.options.multirng:
-        pytest.skip('Skipping multirng mode as not necessary to test with both single and multi RNGs')
     resultsdir = tempdir / 'dataset_results'
     resultsdir.mkdir(exist_ok=True)
     for p_death in [0.25, 0.5]:
diff --git a/tests/test_syphilis.py b/tests/test_syphilis.py
index 4c9eb4ea..eae90a10 100644
--- a/tests/test_syphilis.py
+++ b/tests/test_syphilis.py
@@ -9,7 +9,8 @@
 import matplotlib.pyplot as plt
 
 quick_run = True
-ss.options['multirng']=False
+
+
 def make_syph_sim(dt=1, n_agents=500):
     """ Make a sim with syphilis - used by several subsequent tests """
     syph = ss.Syphilis()
@@ -29,8 +30,8 @@ def make_syph_sim(dt=1, n_agents=500):
     # Marital
     mf = ss.MFNet(
         pars = dict(
-            duration = ss.lognorm_mean(mean=1/24, stdev=0.5),
-            acts = ss.lognorm_mean(mean=80, stdev=30),
+            duration = ss.lognorm_ex(mean=1/24, stdev=0.5),
+            acts = ss.lognorm_ex(mean=80, stdev=30),
         )
     )
     maternal = ss.MaternalNet()
diff --git a/tests/test_worlds.py b/tests/test_worlds.py
new file mode 100644
index 00000000..4ab92dd0
--- /dev/null
+++ b/tests/test_worlds.py
@@ -0,0 +1,121 @@
+""" Test minimally different worlds """
+
+import sciris as sc
+import numpy as np
+import starsim as ss
+import pylab as pl
+
+
+class CountInf(ss.Intervention):
+    """ Store every infection state in a timepoints x people array """
+    def initialize(self, sim):
+        n_agents = len(sim.people)
+        self.arr = np.zeros((sim.npts, n_agents))
+        self.n_agents = n_agents
+        return
+    
+    def apply(self, sim):
+        self.arr[sim.ti, :] = np.array(sim.diseases.sir.infected)[:self.n_agents]
+        return
+
+
+class OneMore(ss.Intervention):
+    """ Add one additional agent and infection """
+    def apply(self, sim):
+        if sim.ti == 0:
+            # Create an extra agent
+            preg = ss.Pregnancy(rel_fertility=0) # Ensure no default births
+            preg.initialize(sim)
+            new_uids = np.array([len(sim.people)]) # Hack since make_embryos doesn't return UIDs
+            preg.make_embryos(sim, np.array([0])) # Assign 0th agent to be the "mother"
+            assert len(new_uids) == 1
+            sim.people.age[new_uids] = -100 # Set to a very low number to never reach debut age
+            
+            # Infect that agent
+            sir = sim.diseases.sir
+            sir.set_prognoses(sim, new_uids)
+            sir.ti_recovered[new_uids] = sim.ti + 1 # Reset recovery time to next timestep
+            
+            # Reset the random states
+            p = sir.pars
+            for dist in [p.dur_inf, p.p_death]:
+                dist.jump(sim.ti+1)
+
+        return
+
+
+def plot_infs(s1, s2):
+    """ Compare infection arrays from two sims """
+    a1 = s1.interventions.countinf.arr
+    a2 = s2.interventions.countinf.arr
+    
+    fig = pl.figure()
+    pl.subplot(1,3,1)
+    pl.pcolormesh(a1.T)
+    pl.xlabel('Timestep')
+    pl.ylabel('Person')
+    pl.title('Baseline')
+    
+    pl.subplot(1,3,2)
+    pl.pcolormesh(a2.T)
+    pl.title('OneMore')
+    
+    pl.subplot(1,3,3)
+    pl.pcolormesh(a2.T - a1.T)
+    pl.title('Difference')
+    
+    sc.figlayout()
+    return fig
+
+
+def test_worlds(do_plot=False):
+    
+    res = sc.objdict()
+    
+    pars = dict(
+        start = 2000,
+        end = 2100,
+        n_agents = 200,
+        verbose = 0.05,
+        diseases = dict(
+            type = 'sir',
+            init_prev = 0.1,
+            beta = 1.0,
+            dur_inf = 20,
+            p_death = 0, # Here since analyzer can't handle variable numbers of people
+        ),
+        networks = dict(
+            type = 'embedding',
+            duration = 5, # Must be shorter than dur_inf for SIR transmission to occur
+        ),
+    )
+    s1 = ss.Sim(pars=pars, interventions=CountInf())
+    s2 = ss.Sim(pars=pars, interventions=[CountInf(), OneMore()])
+    
+    s1.run()
+    s2.run()
+    
+    sum1 = s1.summarize()
+    sum2 = s2.summarize()
+    res.sum1 = sum1
+    res.sum2 = sum2
+    
+    if do_plot:
+        s1.plot()
+        plot_infs(s1, s2)
+        pl.show()
+    
+    assert len(s2.people) == len(s1.people) + 1
+    assert sum2.sir_cum_infections == sum1.sir_cum_infections + 1
+    assert (s1.interventions.countinf.arr != s2.interventions.countinf.arr).sum() == 0
+        
+    return res
+    
+
+if __name__ == '__main__':
+    T = sc.timer()
+    do_plot = True
+    
+    res = test_worlds(do_plot=do_plot)
+    
+    T.toc()
\ No newline at end of file
diff --git a/tests/update_baseline.py b/tests/update_baseline.py
old mode 100644
new mode 100755