Repository: starsimhub/starsim (language: python)

Issue:
Per-individual distribution parameters with rng slotting
With scipy now providing distribution support, at least on rng-merge-scipy, rng slotting has run into a roadblock for the case in which a user has provided unique parameters for each uid.

What we did before was to sample a large number of random numbers to subsequently select the random numbers we wanted (using slots). Two agents with the same slot would have received the same random number draw. These random numbers were then used to sample from a standardized form of the distribution, like a normal with mu=1 and sigma=1. Finally, the per-agent parameters were applied to transform the resulting draws.

Now with scipy providing the frozen distributions, the user calls `rvs` or `filter` and provides uids and either an array of parameters or a callable that returns an array of parameters. (This for the case of per-agent distribution parameters.) To account for slotting, the current approach expands the parameters and fills BY SLOT with the user-provided values:
```python
pars_slots = np.full(n_samples, fill_value=self.fill_value, dtype=kwargs[pname].dtype)
pars_slots[slots] = kwargs[pname] # <-- PROBLEM
kwargs[pname] = pars_slots
```
where `kwargs[pname]` is the array of parameters.

Finally, we ask scipy to do the heavy lifting: `vals = super().rvs(*args, **kwargs)`.

The PROBLEM occurs in the unlikely event that two agents are assigned the same slot. With the approach above, these two agents _must share parameter values too!_ This is clearly a big bad bug.

Possible solutions include:
1. Do not allow slots to be reused. However, doing so would knock off common random numbers.
2. When possible, draw from standard distributions and apply parameters later. However, not possible for all distributions.
3. Scipy `rvs` is calling into the MultiRNG generator, perhaps there's a way to catch it there?
4. Make several calls to `.rvs` to evaluate repeat slots, resetting the MultiRNG each time. The number of calls would have to equal the max slot multiplicity.
Use numpy normal in streams.py
np.random.normal(size=size, loc=mu, scale=std)
Rng merge states ck
Replaces #162 with a merge keeping history but generally resolving in favor of `rng-merge-scipy`.
Record cause of death
Considering we have demographics as well as multiple disease modules, it will be useful to record which module (or demographics) was responsible for the death of each agent during the course of a simulation.
Consider making check_requires more robust
Current code is based on the class type:
```python
def check_requires(self, sim):
        for req in self.requires:
            if req not in [m.__class__ for m in sim.modules]:
                raise Exception(f'{self.name} (label={self.label}) requires module {req} but the Sim did not contain a module of this type.')
        return
```

Consider connecting instead by name or label, e.g. in case a user wants to connect two module instances of the same class.
All rates need to be multiplied by dt to get probs
Looking specifically at transmission within `modules.py`, and not seeing how dt enters into the probability of transmission - it should! All rates must be multiplied by `dt` to get per-timestep probabilities.

Technically it's not correct to set `prob = rate * dt`, unless either the `rate` or `dt` are small. The correct conversion from continuous to discrete time involves an exponential, which is more costly. The higher-order terms are insignificant if the step is small.

I have seen some placed in the code that already account for `dt`, as in some of the demographics. But I suspect there are other (important) places where `dt` has been left off.
Per-individual distribution parameters with rng slotting
With scipy now providing distribution support, at least on rng-merge-scipy, rng slotting has run into a roadblock for the case in which a user has provided unique parameters for each uid.

What we did before was to sample a large number of random numbers to subsequently select the random numbers we wanted (using slots). Two agents with the same slot would have received the same random number draw. These random numbers were then used to sample from a standardized form of the distribution, like a normal with mu=1 and sigma=1. Finally, the per-agent parameters were applied to transform the resulting draws.

Now with scipy providing the frozen distributions, the user calls `rvs` or `filter` and provides uids and either an array of parameters or a callable that returns an array of parameters. (This for the case of per-agent distribution parameters.) To account for slotting, the current approach expands the parameters and fills BY SLOT with the user-provided values:
```python
pars_slots = np.full(n_samples, fill_value=self.fill_value, dtype=kwargs[pname].dtype)
pars_slots[slots] = kwargs[pname] # <-- PROBLEM
kwargs[pname] = pars_slots
```
where `kwargs[pname]` is the array of parameters.

Finally, we ask scipy to do the heavy lifting: `vals = super().rvs(*args, **kwargs)`.

The PROBLEM occurs in the unlikely event that two agents are assigned the same slot. With the approach above, these two agents _must share parameter values too!_ This is clearly a big bad bug.

Possible solutions include:
1. Do not allow slots to be reused. However, doing so would knock off common random numbers.
2. When possible, draw from standard distributions and apply parameters later. However, not possible for all distributions.
3. Scipy `rvs` is calling into the MultiRNG generator, perhaps there's a way to catch it there?
4. Make several calls to `.rvs` to evaluate repeat slots, resetting the MultiRNG each time. The number of calls would have to equal the max slot multiplicity.
Per-individual distribution parameters with rng slotting
With scipy now providing distribution support, at least on rng-merge-scipy, rng slotting has run into a roadblock for the case in which a user has provided unique parameters for each uid.

What we did before was to sample a large number of random numbers to subsequently select the random numbers we wanted (using slots). Two agents with the same slot would have received the same random number draw. These random numbers were then used to sample from a standardized form of the distribution, like a normal with mu=1 and sigma=1. Finally, the per-agent parameters were applied to transform the resulting draws.

Now with scipy providing the frozen distributions, the user calls `rvs` or `filter` and provides uids and either an array of parameters or a callable that returns an array of parameters. (This for the case of per-agent distribution parameters.) To account for slotting, the current approach expands the parameters and fills BY SLOT with the user-provided values:
```python
pars_slots = np.full(n_samples, fill_value=self.fill_value, dtype=kwargs[pname].dtype)
pars_slots[slots] = kwargs[pname] # <-- PROBLEM
kwargs[pname] = pars_slots
```
where `kwargs[pname]` is the array of parameters.

Finally, we ask scipy to do the heavy lifting: `vals = super().rvs(*args, **kwargs)`.

The PROBLEM occurs in the unlikely event that two agents are assigned the same slot. With the approach above, these two agents _must share parameter values too!_ This is clearly a big bad bug.

Possible solutions include:
1. Do not allow slots to be reused. However, doing so would knock off common random numbers.
2. When possible, draw from standard distributions and apply parameters later. However, not possible for all distributions.
3. Scipy `rvs` is calling into the MultiRNG generator, perhaps there's a way to catch it there?
4. Make several calls to `.rvs` to evaluate repeat slots, resetting the MultiRNG each time. The number of calls would have to equal the max slot multiplicity.
Update msm and fix mf_msm connector
mf_msm connector code is out of date an no longer functional. Also needs a test, not just a devtest.


Discussion / hints:

Fix patch:
diff --git a/.gitignore b/.gitignore
index 34a00e22..063dbe1d 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1,6 +1,9 @@
 # HPVsim-generated files
 *.obj
 files/
+tests/figs/*
+tests/devtests/figs/*
+figs/*
 modules.rst
 
 # Other files
diff --git a/stisim/__init__.py b/stisim/__init__.py
index f709e345..15440353 100644
--- a/stisim/__init__.py
+++ b/stisim/__init__.py
@@ -4,6 +4,7 @@
 from .utils         import *
 from .distributions import *
 from .states        import *
+from .random        import *
 from .people        import *
 from .modules       import *
 from .networks      import *
diff --git a/stisim/analyzers.py b/stisim/analyzers.py
index d86d1b9e..c24a57d3 100644
--- a/stisim/analyzers.py
+++ b/stisim/analyzers.py
@@ -14,5 +14,4 @@ def __call__(self, *args, **kwargs):
         return self.apply(*args, **kwargs)
 
     def update_results(self, sim):
-        raise NotImplementedError
-
+        raise NotImplementedError
\ No newline at end of file
diff --git a/stisim/connectors.py b/stisim/connectors.py
index 17eefdce..daad46cd 100644
--- a/stisim/connectors.py
+++ b/stisim/connectors.py
@@ -7,9 +7,7 @@
 
 class Connector(ss.Module):
     def __init__(self, pars=None, diseases=None, *args, **kwargs):
-        self.pars = ss.omerge(pars)
-        self.results = ss.ndict()
-        self.diseases = ss.ndict(diseases)
+        super().__init__(pars=pars, requires=diseases, *args, **kwargs)
         return
 
 
@@ -19,22 +17,13 @@ class simple_hiv_ng(Connector):
     """ Simple connector whereby rel_sus to NG doubles if CD4 count is <200"""
 
     def __init__(self, pars=None):
-        super().__init__(pars=pars)
+        super().__init__(pars=pars, label='HIV-Gonorrhea', diseases=[ss.HIV, ss.Gonorrhea])
         self.pars = ss.omerge({
             'rel_trans_hiv': 2,
             'rel_trans_aids': 5,
             'rel_sus_hiv': 2,
             'rel_sus_aids': 5,
-        }, self.pars)  # Could add pars here
-        self.diseases = ['hiv', 'gonorrhea']
-        return
-
-    def initialize(self, sim):
-        # Make sure the sim has the modules that this connector deals with
-        # TODO replace this placeholder code with something robust.
-        if not set(self.diseases).issubset(sim.diseases.keys()):
-            errormsg = f'Missing required diseases {set(self.diseases).difference(sim.diseases.keys())}'
-            raise ValueError(errormsg)
+        }, self.pars)
         return
 
     def update(self, sim):
@@ -46,4 +35,4 @@ def update(self, sim):
         sim.people.gonorrhea.rel_trans[sim.people.hiv.cd4 < 500] = self.pars.rel_trans_hiv
         sim.people.gonorrhea.rel_trans[sim.people.hiv.cd4 < 200] = self.pars.rel_trans_aids
 
-        return
+        return
\ No newline at end of file
diff --git a/stisim/demographics.py b/stisim/demographics.py
index 2d86ffbd..4bf115be 100644
--- a/stisim/demographics.py
+++ b/stisim/demographics.py
@@ -6,6 +6,7 @@
 import stisim as ss
 import sciris as sc
 import pandas as pd
+import scipy.stats as sps
 
 __all__ = ['DemographicModule', 'births', 'background_deaths', 'Pregnancy']
 
@@ -68,9 +69,9 @@ def init_results(self, sim):
         return
 
     def update(self, sim):
-        n_new = self.add_births(sim)
-        self.update_results(n_new, sim)
-        return
+        new_uids = self.add_births(sim)
+        self.update_results(len(new_uids), sim)
+        return new_uids
 
     def get_birth_rate(self, sim):
         """
@@ -87,8 +88,8 @@ def get_birth_rate(self, sim):
     def add_births(self, sim):
         # Add n_new births to each state in the sim
         n_new = self.get_birth_rate(sim)
-        sim.people.grow(n_new)
-        return n_new
+        new_uids = sim.people.grow(n_new)
+        return new_uids
 
     def update_results(self, n_new, sim):
         self.results['new'][sim.ti] = n_new
@@ -96,7 +97,7 @@ def update_results(self, n_new, sim):
     def finalize(self, sim):
         super().finalize(sim)
         self.results['cumulative'] = np.cumsum(self.results['new'])
-        self.results['cbr'] = self.results['new'] / sim.results['pop_size']
+        self.results['cbr'] = np.divide(self.results['new'], sim.results['n_alive'], where=sim.results['n_alive']>0)
 
 
 class background_deaths(DemographicModule):
@@ -104,75 +105,26 @@ def __init__(self, pars=None):
         super().__init__(pars)
 
         self.pars = ss.omerge({
-            'death_rates': 0,
-            'rel_death': 1,
-            'data_cols': {'year': 'Time', 'sex': 'Sex', 'age': 'AgeGrpStart', 'value': 'mx'},
-            'sex_keys': {'f': 'Female', 'm': 'Male'},
-            'units_per_100': 1e-3  # assumes birth rates are per 1000. If using percentages, switch this to 1
+            #'rel_death': 1,
+            #'units_per_100': 1e-3,  # assumes birth rates are per 1000. If using percentages, switch this to 1
+            'death_prob_func': self.death_prob,
         }, self.pars)
-
-        # Validate death rate inputs
-        self.set_death_rates(pars['death_rates'])
-
-        # Set death probs
-        self.death_probs = ss.State('death_probs', float, 0)
-
-    def set_death_rates(self, death_rates):
-        """Standardize/validate death rates"""
-
-        if sc.checktype(death_rates, pd.DataFrame):
-            if not set(self.pars.data_cols.values()).issubset(death_rates.columns):
-                errormsg = 'Please ensure the columns of the death rate data match the values in pars.data_cols.'
-                raise ValueError(errormsg)
-            df = death_rates
-
-        elif sc.checktype(death_rates, pd.Series):
-            if (death_rates.index < 120).all():  # Assume index is age bins
-                df = pd.DataFrame({
-                    self.pars.data_cols['year']: 2000,
-                    self.pars.data_cols['age']: death_rates.index.values,
-                    self.pars.data_cols['value']: death_rates.values,
-                })
-            elif (death_rates.index > 1900).all():  # Assume index year
-                df = pd.DataFrame({
-                    self.pars.data_cols['year']: death_rates.index.values,
-                    self.pars.data_cols['age']: 0,
-                    self.pars.data_cols['value']: death_rates.values,
-
-                })
-            else:
-                errormsg = 'Could not understand index of death rate series: should be age or year.'
-                raise ValueError(errormsg)
-
-            df = pd.concat([df, df])
-            df[self.pars.data_cols['sex']] = np.repeat(list(self.pars.sex_keys.values()), len(death_rates))
-
-        elif sc.checktype(death_rates, dict):
-            if not set(self.pars.data_cols.values()).issubset(death_rates.keys()):
-                errormsg = 'Please ensure the keys of the death rate data dict match the values in pars.data_cols.'
-                raise ValueError(errormsg)
-            df = pd.DataFrame(death_rates)
-
-        elif sc.isnumber(death_rates):
-            df = pd.DataFrame({
-                self.pars.data_cols['year']: [2000, 2000],
-                self.pars.data_cols['age']: [0, 0],
-                self.pars.data_cols['sex']: self.pars.sex_keys.values(),
-                self.pars.data_cols['value']: [death_rates, death_rates],
-            })
-
-        else:
-            errormsg = f'Death rate data type {type(death_rates)} not understood.'
-            raise ValueError(errormsg)
-
-        self.pars.death_rates = df
-
+        self.death_prob_dist = sps.bernoulli(p=self.pars['death_prob_func'])
         return
 
+    @staticmethod
+    def death_prob(self, sim, uids):
+        #p = self.pars
+        #year = p.death_rates[p.data_cols['year']]
+        #val = p.death_rates[p.data_cols['cbr']]
+        #this_death_rate = np.interp(sim.year, year, val) * p.units_per_100 * p.rel_death * sim.pars.dt_demog
+        #n_new = int(np.floor(np.count_nonzero(sim.people.alive) * this_death_rate))
+        return 0.01
+
     def init_results(self, sim):
         self.results += ss.Result(self.name, 'new', sim.npts, dtype=int)
         self.results += ss.Result(self.name, 'cumulative', sim.npts, dtype=int)
-        self.results += ss.Result(self.name, 'mortality_rate', sim.npts, dtype=int)
+        self.results += ss.Result(self.name, 'cmr', sim.npts, dtype=int)
         return
 
     def update(self, sim):
@@ -182,31 +134,8 @@ def update(self, sim):
 
     def apply_deaths(self, sim):
         """ Select people to die """
-
-        p = self.pars
-        year_label = p.data_cols['year']
-        age_label = p.data_cols['age']
-        sex_label = p.data_cols['sex']
-        val_label = p.data_cols['value']
-        sex_keys = p.sex_keys
-
-        available_years = p.death_rates[year_label].unique()
-        year_ind = sc.findnearest(available_years, sim.year)
-        nearest_year = available_years[year_ind]
-
-        df = p.death_rates.loc[p.death_rates[year_label] == nearest_year]
-        age_bins = df[age_label].unique()
-        age_inds = np.digitize(sim.people.age, age_bins) - 1
-
-        f_arr = df[val_label].loc[df[sex_label] == sex_keys['f']].values
-        m_arr = df[val_label].loc[df[sex_label] == sex_keys['m']].values
-        self.death_probs[sim.people.female] = f_arr[age_inds[sim.people.female]]
-        self.death_probs[sim.people.male] = m_arr[age_inds[sim.people.male]]
-        self.death_probs *= p.rel_death  # Adjust overall death probabilities
-
-        # Get indices of people who die of other causes
-        death_uids = ss.true(ss.binomial_arr(self.death_probs))
-        death_uids = ss.true(sim.people.alive[death_uids])
+        alive_uids = ss.true(sim.people.alive)
+        death_uids = self.death_prob_dist.filter(alive_uids)
         sim.people.request_death(death_uids)
 
         return len(death_uids)
@@ -216,7 +145,7 @@ def update_results(self, n_deaths, sim):
 
     def finalize(self, sim):
         self.results['cumulative'] = np.cumsum(self.results['new'])
-        self.results['mortality_rate'] = self.results['new'] / sim.results['pop_size']
+        self.results['cmr'] = np.divide(self.results['new'], sim.results['n_alive'], where=sim.results['n_alive']>0)
 
 
 class Pregnancy(DemographicModule):
@@ -237,15 +166,19 @@ def __init__(self, pars=None):
         self.pars = ss.omerge({
             'dur_pregnancy': 0.75,  # Make this a distribution?
             'dur_postpartum': 0.5,  # Make this a distribution?
-            'inci': 0.03,  # Replace this with age-specific rates
-            'p_death': 0,  # Probability of maternal death. Question, should this be linked to age and/or duration?
-            'initial': 3,  # Number of women initially pregnant
+            'pregnancy_prob_per_dt': sps.bernoulli(p=0.3),  # Probabilty of acquiring pregnancy on each time step. Can replace with callable parameters for age-specific rates, etc. NOTE: You will manually need to adjust for the simulation timestep dt!
+            'p_death': sps.bernoulli(p=0.15),  # Probability of maternal death.
+            'p_female': sps.bernoulli(p=0.5),
+            'init_prev': 0.3,  # Number of women initially pregnant # TODO: Default value
         }, self.pars)
 
+        self.choose_slots = sps.randint(low=0, high=1) # Low and high will be reset upon initialization
         return
 
     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'])
         sim.pars['birth_rates'] = None  # This turns off birth rate pars so births only come from this module
         return
 
@@ -302,32 +235,38 @@ def make_pregnancies(self, sim):
 
         # If incidence of pregnancy is non-zero, make some cases
         # Think about how to deal with age/time-varying fertility
-        if self.pars.inci > 0:
-            denom_conds = ppl.female & ppl.active & self.susceptible
-            inds_to_choose_from = ss.true(denom_conds)
-            uids = ss.binomial_filter(self.pars.inci, inds_to_choose_from)
-
-            # Add UIDs for the as-yet-unborn agents so that we can track prognoses and transmission patterns
-            n_unborn_agents = len(uids)
-            if n_unborn_agents > 0:
-                # Grow the arrays and set properties for the unborn agents
-                new_uids = sim.people.grow(n_unborn_agents)
-                sim.people.age[new_uids] = -self.pars.dur_pregnancy
-
-                # Add connections to any vertical transmission layers
-                # Placeholder code to be moved / refactored. The maternal network may need to be
-                # handled separately to the sexual networks, TBC how to handle this most elegantly
-                for lkey, layer in sim.people.networks.items():
-                    if layer.vertical:  # What happens if there's more than one vertical layer?
-                        durs = np.full(n_unborn_agents, fill_value=self.pars.dur_pregnancy + self.pars.dur_postpartum)
-                        layer.add_pairs(uids, new_uids, dur=durs)
-
-                # Set prognoses for the pregnancies
-                self.set_prognoses(sim, uids)
+        denom_conds = ppl.female & ppl.active & self.susceptible
+        inds_to_choose_from = ss.true(denom_conds)
+        uids = self.pars['pregnancy_prob_per_dt'].filter(inds_to_choose_from)
+
+        # Add UIDs for the as-yet-unborn agents so that we can track prognoses and transmission patterns
+        n_unborn_agents = len(uids)
+        if n_unborn_agents > 0:
+
+            # Choose slots for the unborn agents
+            new_slots = self.choose_slots.rvs(uids)
+
+            # Grow the arrays and set properties for the unborn agents
+            new_uids = sim.people.grow(len(new_slots))
+
+            sim.people.age[new_uids] = -self.pars.dur_pregnancy
+            sim.people.slot[new_uids] = new_slots # Before sampling female_dist
+            sim.people.female[new_uids] = self.pars['p_female'].rvs(uids)
+
+            # Add connections to any vertical transmission layers
+            # Placeholder code to be moved / refactored. The maternal network may need to be
+            # handled separately to the sexual networks, TBC how to handle this most elegantly
+            for lkey, layer in sim.people.networks.items():
+                if layer.vertical:  # What happens if there's more than one vertical layer?
+                    durs = np.full(n_unborn_agents, fill_value=self.pars.dur_pregnancy + self.pars.dur_postpartum)
+                    layer.add_pairs(uids, new_uids, dur=durs)
+
+            # Set prognoses for the pregnancies
+            self.set_prognoses(sim, uids) # Could set from_uids to network partners?
 
         return
 
-    def set_prognoses(self, sim, uids):
+    def set_prognoses(self, sim, uids, from_uids=None):
         """
         Make pregnancies
         Add miscarriage/termination logic here
@@ -342,7 +281,7 @@ def set_prognoses(self, sim, uids):
 
         # Outcomes for pregnancies
         dur = np.full(len(uids), sim.ti + self.pars.dur_pregnancy / sim.dt)
-        dead = np.random.random(len(uids)) < self.pars.p_death
+        dead = self.pars.p_death.rvs(uids)
         self.ti_delivery[uids] = dur  # Currently assumes maternal deaths still result in a live baby
         dur_post_partum = np.full(len(uids), dur + self.pars.dur_postpartum / sim.dt)
         self.ti_postpartum[uids] = dur_post_partum
diff --git a/stisim/diseases/__init__.py b/stisim/diseases/__init__.py
index 673a487a..0b0629ba 100644
--- a/stisim/diseases/__init__.py
+++ b/stisim/diseases/__init__.py
@@ -1,4 +1,4 @@
 from .disease       import *
 from .hiv           import *
 from .gonorrhea     import *
-from .examples      import *
+from .examples      import *
\ No newline at end of file
diff --git a/stisim/diseases/disease.py b/stisim/diseases/disease.py
index 929e485e..49c16135 100644
--- a/stisim/diseases/disease.py
+++ b/stisim/diseases/disease.py
@@ -5,6 +5,7 @@
 import numpy as np
 import sciris as sc
 import stisim as ss
+import scipy.stats as sps
 import networkx as nx
 from operator import itemgetter
 import pandas as pd
@@ -103,11 +104,11 @@ def __init__(self, *args, **kwargs):
     def _boolean_states(self):
         """
         Iterator over states with boolean type
-        
+
         For diseases, these states typically represent attributes like 'susceptible',
-        'infectious', 'diagnosed' etc. These variables are typically useful to  
-        
-        :return: 
+        'infectious', 'diagnosed' etc. These variables are typically useful to
+
+        :return:
         """
         for state in self.states:
             if state.dtype == bool:
@@ -248,7 +249,6 @@ def update_results(self, sim):
         return
 
 
-
 class STI(Disease):
     """
     Base class for STIs used in STIsim
@@ -260,12 +260,13 @@ class STI(Disease):
 
     def __init__(self, *args, **kwargs):
         super().__init__(*args, **kwargs)
-        self.rel_sus = ss.State('rel_sus', float, 1)
-        self.rel_sev = ss.State('rel_sev', float, 1)
-        self.rel_trans = ss.State('rel_trans', float, 1)
+        self.rel_sus     = ss.State('rel_sus', float, 1)
+        self.rel_sev     = ss.State('rel_sev', float, 1)
+        self.rel_trans   = ss.State('rel_trans', float, 1)
         self.susceptible = ss.State('susceptible', bool, True)
-        self.infected = ss.State('infected', bool, False)
+        self.infected    = ss.State('infected', bool, False)
         self.ti_infected = ss.State('ti_infected', int, ss.INT_NAN)
+
         return
 
     def validate_pars(self, sim):
@@ -284,9 +285,13 @@ def set_initial_states(self, sim):
         i.e., creating their dynamic array, linking them to a People instance. That should have already
         taken place by the time this method is called.
         """
-        n_init_cases = int(self.pars['init_prev'] * len(sim.people))
-        initial_cases = np.random.choice(sim.people.uid, n_init_cases, replace=False)
-        self.set_prognoses(sim, initial_cases)
+        if self.pars['seed_infections'] is None:
+            return
+
+        alive_uids = ss.true(sim.people.alive) # Maybe just sim.people.uid?
+        initial_cases = self.pars['seed_infections'].filter(alive_uids)
+
+        self.set_prognoses(sim, initial_cases, from_uids=None)  # TODO: sentinel value to indicate seeds?
         return
 
     def init_results(self, sim):
@@ -307,22 +312,132 @@ def update_pre(self, sim):
         """
         pass
 
-    def make_new_cases(self, sim):
-        """ Add new cases of module, through transmission, incidence, etc. """
-        pars = self.pars
+    def _make_new_cases_singlerng(self, sim):
+        # Not common-random-number-safe, but more efficient for when not using the multirng feature
+        new_cases = []
+        sources = []
         for k, layer in sim.people.networks.items():
-            if k in pars['beta']:
+            if k in self.pars['beta']:
                 contacts = layer.contacts
-                rel_trans = (self.infected & sim.people.alive).astype(float) * self.rel_trans
-                rel_sus = (self.susceptible & sim.people.alive).astype(float) * self.rel_sus
-                for a, b, beta in [[contacts.p1, contacts.p2, pars.beta[k][0]],
-                                   [contacts.p2, contacts.p1, pars.beta[k][1]]]:
+                rel_trans = (self.infected & sim.people.alive) * self.rel_trans
+                rel_sus = (self.susceptible & sim.people.alive) * self.rel_sus
+                for a, b, beta in [[contacts.p1, contacts.p2, self.pars.beta[k][0]],
+                                   [contacts.p2, contacts.p1, self.pars.beta[k][1]]]:
+                    if beta == 0:
+                        continue
                     # probability of a->b transmission
-                    p_transmit = rel_trans[a] * rel_sus[b] * contacts.beta * beta
-                    new_cases = np.random.random(len(a)) < p_transmit
-                    if np.any(new_cases):
-                        self.set_prognoses(sim, b[new_cases], a[new_cases])
+                    p_transmit = rel_trans[a] * rel_sus[b] * contacts.beta * beta # TODO: Needs 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])
+        return np.concatenate(new_cases), np.concatenate(sources)
+
+    def _choose_new_cases_multirng(self, people):
+        '''
+        Common-random-number-safe transmission code works by computing the
+        probability of each _node_ acquiring a case rather than checking if each
+        _edge_ transmits.
+        '''
+        n = len(people.uid) # TODO: possibly could be shortened to just the people who are alive
+        p_acq_node = np.zeros(n)
+
+        for lkey, layer in people.networks.items():
+            if lkey in self.pars['beta']:
+                contacts = layer.contacts
+                rel_trans = self.rel_trans * (self.infected & people.alive)
+                rel_sus = self.rel_sus * (self.susceptible & people.alive)
+
+                a_to_b = [contacts.p1, contacts.p2, self.pars.beta[lkey][0]]
+                b_to_a = [contacts.p2, contacts.p1, self.pars.beta[lkey][1]]
+                for a, b, beta in [a_to_b, b_to_a]: # Transmission from a --> b
+                    if beta == 0:
+                        continue
+                    
+                    # Accumulate acquisition probability for each person (node)
+                    node_from_edge = np.ones((n, len(a)))
+                    bi = people._uid_map[b] # Indices of b (rather than uid)
+                    node_from_edge[bi, np.arange(len(a))] = 1 - rel_trans[a] * rel_sus[b] * contacts.beta * beta # TODO: Needs DT
+                    p_acq_node = 1 - (1-p_acq_node) * node_from_edge.prod(axis=1) # (1-p1)*(1-p2)*...
+
+
+        # Slotted draw, need to find a long-term place for this logic
+        slots = people.slot[people.uid]
+        #p = np.full(np.max(slots)+1, 0, dtype=p_acq_node.dtype)
+        #p[slots] = p_acq_node
+        #new_cases_bool = sps.bernoulli.rvs(p=p).astype(bool)[slots]
+        new_cases_bool = sps.uniform.rvs(size=np.max(slots)+1)[slots] < p_acq_node
+        new_cases = people.uid[new_cases_bool]
+        return new_cases
+
+    def _determine_case_source_multirng(self, people, new_cases):
+        '''
+        Given the uids of new cases, determine which agents are the source of each case
+        '''
+        sources = np.zeros_like(new_cases)
+
+        # Slotted draw, need to find a long-term place for this logic
+        slots = people.slot[new_cases]
+        r = sps.uniform.rvs(size=np.max(slots)+1)[slots]
+
+        for i, uid in enumerate(new_cases):
+            p_acqs = []
+            possible_sources = []
+
+            for lkey, layer in people.networks.items():
+                if lkey in self.pars['beta']:
+                    contacts = layer.contacts
+                    a_to_b = [contacts.p1, contacts.p2, self.pars.beta[lkey][0]]
+                    b_to_a = [contacts.p2, contacts.p1, self.pars.beta[lkey][1]]
+                    for a, b, beta in [a_to_b, b_to_a]: # Transmission from a --> b
+                        if beta == 0:
+                            continue
+                    
+                        inds = np.where(b == uid)[0]
+                        if len(inds) == 0:
+                            continue
+                        neighbors = a[inds]
+
+                        rel_trans = self.rel_trans[neighbors] * (self.infected[neighbors] & people.alive[neighbors])
+                        rel_sus = self.rel_sus[uid] * (self.susceptible[uid] & people.alive[uid])
+                        beta_combined = contacts.beta[inds] * beta
+
+                        # Compute acquisition probabilities from neighbors --> uid
+                        # TODO: Could remove zeros
+                        p_acqs.append((rel_trans * rel_sus * beta_combined).__array__()) # Needs DT
+                        possible_sources.append(neighbors.__array__())
+
+            # Concatenate across layers and directions (p1->p2 vs p2->p1)
+            p_acqs = np.concatenate(p_acqs)
+            possible_sources = np.concatenate(possible_sources)
+
+            if len(possible_sources) == 1: # Easy if only one possible source
+                sources[i] = possible_sources[0]
+            else:
+                # Roulette selection using slotted draw r associated with this new case
+                cumsum = p_acqs / p_acqs.sum()
+                source_idx = np.argmax(cumsum >= r[i])
+                sources[i] = possible_sources[source_idx]
+        return sources
 
+    def make_new_cases(self, sim):
+        """ Add new cases of module, through transmission, incidence, etc. """
+        if not ss.options.multirng:
+            # Determine new cases for singlerng
+            new_cases, sources = self._make_new_cases_singlerng(sim)
+        else:
+            # Determine new cases for multirng
+            new_cases = self._choose_new_cases_multirng(sim.people)
+            if len(new_cases):
+                # Now determine whom infected each case
+                sources = self._determine_case_source_multirng(sim.people, new_cases)
+
+        if len(new_cases):
+            self.set_prognoses(sim, new_cases, sources)
+        
+        return len(new_cases) # number of new cases made
+
+    def set_prognoses(self, sim, uids, from_uids):
+        pass
 
     def set_congenital(self, sim, uids):
         # Need to figure out whether we would have a methods like this here or make it
@@ -333,3 +448,4 @@ def update_results(self, sim):
         super().update_results(sim)
         self.results['prevalence'][sim.ti] = self.results.n_infected[sim.ti] / np.count_nonzero(sim.people.alive)
         self.results['new_infections'][sim.ti] = np.count_nonzero(self.ti_infected == sim.ti)
+        return
diff --git a/stisim/diseases/examples.py b/stisim/diseases/examples.py
index 2ffcd3e9..1ea1ac0d 100644
--- a/stisim/diseases/examples.py
+++ b/stisim/diseases/examples.py
@@ -6,6 +6,7 @@
 import stisim as ss
 import sciris as sc
 from .disease import Disease
+import scipy.stats as sps
 
 class SIR(Disease):
     """
@@ -13,31 +14,32 @@ class SIR(Disease):
 
     This class implements a basic SIR model with states for susceptible,
     infected/infectious, and recovered. It also includes deaths, and basic
-    results
-    """
+    results.
 
-    default_pars = {
-        'dur_inf': ss.normal_pos(mean=5, std=1),
-        'p_death': 0.2,
-        'initial': 3,
-        'beta': None,
-    }
+    Note that this class is not fully compatible with common random numbers.
+    """
 
     def __init__(self, pars=None, *args, **kwargs):
+        default_pars = {
+            #'dur_inf': sps.weibull_min(c=lambda self, sim, uids: sim.people.age[uids], scale=10),#, seed='Duration of SIR Infection'),
+            #'dur_inf': sps.norm(loc=lambda self, sim, uids: sim.people.age[uids], scale=2),
+            'dur_inf': sps.lognorm(s=1, loc=10),
+            'seed_infections': sps.bernoulli(p=0.1),
+            'death_given_infection': sps.bernoulli(p=0.2),
+            'beta': None,
+        }
 
-        super().__init__(pars=ss.omerge(self.default_pars, pars), *args, **kwargs)
+        super().__init__(pars=ss.omerge(default_pars, pars), *args, **kwargs)
 
         self.susceptible = ss.State('susceptible', bool, True)
         self.infected = ss.State('infected', bool, False)
         self.recovered = ss.State('recovered', bool, False)
-
         self.t_infected = ss.State('t_infected', float, np.nan)
         self.t_recovered = ss.State('t_recovered', float, np.nan)
         self.t_dead = ss.State('t_dead', float, np.nan)
 
         return
 
-
     def init_results(self, sim):
         """
         Initialize results
@@ -45,7 +47,7 @@ def init_results(self, sim):
         super().init_results(sim)
         self.results += ss.Result(self.name, 'prevalence', sim.npts, dtype=float)
         self.results += ss.Result(self.name, 'new_infections', sim.npts, dtype=int)
-
+        return
 
     def update_pre(self, sim):
         # Progress infectious -> recovered
@@ -55,8 +57,9 @@ def update_pre(self, sim):
 
         # Trigger deaths
         deaths = ss.true(self.t_dead <= sim.year)
-        sim.people.request_death(deaths)
-
+        if len(deaths):
+            sim.people.request_death(deaths)
+        return len(deaths)
 
     def update_death(self, sim, uids):
         # Reset infected/recovered flags for dead agents
@@ -80,7 +83,8 @@ def set_initial_states(self, sim):
         i.e., creating their dynamic array, linking them to a People instance. That should have already
         taken place by the time this method is called.
         """
-        initial_cases = np.random.choice(sim.people.uid, self.pars['initial'], replace=False)
+        alive_uids = ss.true(sim.people.alive)
+        initial_cases = self.pars['seed_infections'].filter(alive_uids)
         self.infect(sim, initial_cases, None)
         return
 
@@ -93,17 +97,17 @@ def infect(self, sim, uids, from_uids):
         self.t_infected[uids] = sim.year
 
         # Calculate and schedule future outcomes for recovery/death
-        dur_inf = self.pars['dur_inf'].sample(len(uids))
-        dead = np.random.random(len(uids)) < self.pars.p_death
-        self.t_recovered[uids[~dead]] = sim.year + dur_inf[~dead]
-        self.t_dead[uids[dead]] = sim.year + dur_inf[dead]
+        dur_inf = self.pars['dur_inf'].rvs(uids)
+        will_die = self.pars['death_given_infection'].rvs(uids)
+        self.t_recovered[uids[~will_die]] = sim.year + dur_inf[~will_die]
+        self.t_dead[uids[will_die]] = sim.year + dur_inf[will_die]
 
-        # Update result count of new infections - important to use += because infect() may be called multiple times
-        # per timestep
+        # Update result count of new infections - important to use += because
+        # infect() may be called multiple times per timestep
         self.results['new_infections'][sim.ti] += len(uids)
         return
 
-    def make_new_cases(self, sim):
+    def make_new_cases(self, sim): # TODO: Use function from STI
         for k, layer in sim.people.networks.items():
             if k in self.pars['beta']:
                 rel_trans = (self.infected & sim.people.alive).astype(float)
@@ -112,7 +116,7 @@ def make_new_cases(self, sim):
                                    [layer.contacts['p2'], layer.contacts['p1'], self.pars['beta'][k]]]:
                     # probability of a->b transmission
                     p_transmit = rel_trans[a] * rel_sus[b] * layer.contacts['beta'] * beta * sim.dt
-                    new_cases = np.random.random(len(a)) < p_transmit
+                    new_cases = np.random.random(len(a)) < p_transmit # As this class is not common-random-number safe anyway, calling np.random is perfectly fine!
                     if new_cases.any():
                         self.infect(sim, b[new_cases], a[new_cases])
         return
@@ -122,6 +126,7 @@ def update_results(self, sim):
         self.results['prevalence'][sim.ti] = self.results.n_infected[sim.ti] / np.count_nonzero(sim.people.alive)
         return
 
+
 class NCD(Disease):
     """
     Example non-communicable disease
@@ -130,17 +135,19 @@ class NCD(Disease):
     (e.g., hypertension, diabetes), a state for having the condition, and associated
     mortality.
     """
-    default_pars = {
-        'risk_prev': 0.3, # Initial prevalence of risk factors
-        'p_affected_given_risk': 0.1, # 10% chance per year of acquiring
-        'p_death_given_risk': 0.05,  # 5% chance per year of death
-    }
-
     def __init__(self, pars=None):
-        super().__init__(ss.omerge(self.default_pars, pars))
-        self.at_risk = ss.State('at_risk', bool, False)
-        self.affected = ss.State('affected', bool, False)
-        self.ti_dead = ss.State('ti_dead', float, np.nan)
+        default_pars = {
+            'initial_risk': sps.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': sps.expon(scale=10),
+            'prognosis': sps.weibull_min(c=2, scale=5), # Time in years between first becoming affected and death
+        }
+
+        super().__init__(ss.omerge(default_pars, pars))
+        self.at_risk      = ss.State('at_risk', bool, False)
+        self.affected     = ss.State('affected', bool, False)
+        self.ti_affected  = ss.State('ti_affected', int, ss.INT_NAN)
+        self.ti_dead      = ss.State('ti_dead', int, ss.INT_NAN)
         return
 
     @property
@@ -154,23 +161,28 @@ def set_initial_states(self, sim):
         i.e., creating their dynamic array, linking them to a People instance. That should have already
         taken place by the time this method is called.
         """
-        initial_cases = np.random.choice(sim.people.uid, int(len(sim.people)*self.pars['risk_prev']), replace=False)
-        self.at_risk[initial_cases] = True
-        return initial_cases
+        alive_uids = ss.true(sim.people.alive)
+        initial_risk = self.pars['initial_risk'].filter(alive_uids)
+        self.at_risk[initial_risk] = True
+        self.ti_affected[initial_risk] = sim.ti + sc.randround(self.pars['dur_risk'].rvs(initial_risk) / sim.dt)
+
+        return initial_risk
 
     def update_pre(self, sim):
-        deaths = ss.binomial_filter(sim.dt*self.pars['p_death_given_risk'], ss.true(self.affected))
+        deaths = ss.true(self.ti_dead == sim.ti)
         sim.people.request_death(deaths)
         self.log.add_data(deaths, died=True)
         self.results.new_deaths[sim.ti] = len(deaths) # Log deaths attributable to this module
         return
 
     def make_new_cases(self, sim):
-        new_cases = ss.binomial_filter(self.pars['p_affected_given_risk'], ss.true(self.at_risk))
+        #atrisk_uids = ss.true(self.at_risk)
+        #new_cases = self.pars['affection_rate'].filter(atrisk_uids)
+        new_cases = ss.true(self.ti_affected == sim.ti)
         self.affected[new_cases] = True
-
+        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
 
     def init_results(self, sim):
@@ -187,4 +199,5 @@ def update_results(self, sim):
         super().update_results(sim)
         self.results['n_not_at_risk'][sim.ti] = np.count_nonzero(self.not_at_risk & sim.people.alive)
         self.results['prevalence'][sim.ti] = np.count_nonzero(self.affected & sim.people.alive)/np.count_nonzero(sim.people.alive)
-        return
\ No newline at end of file
+        self.results['new_deaths'][sim.ti] = np.count_nonzero(self.ti_dead == sim.ti)
+        return
diff --git a/stisim/diseases/gonorrhea.py b/stisim/diseases/gonorrhea.py
index 5b19359f..ad971c2e 100644
--- a/stisim/diseases/gonorrhea.py
+++ b/stisim/diseases/gonorrhea.py
@@ -5,6 +5,7 @@
 import numpy as np
 import stisim as ss
 from .disease import STI
+import scipy.stats as sps
 
 
 __all__ = ['Gonorrhea']
@@ -16,16 +17,16 @@ def __init__(self, pars=None):
         super().__init__(pars)
 
         # States additional to the default disease states (see base class)
-        self.symptomatic = ss.State('symptomatic', float, False)
+        self.symptomatic = ss.State('symptomatic', bool, False)
         self.ti_clearance = ss.State('ti_clearance', int, ss.INT_NAN)
         self.p_symp = ss.State('p_symp', float, 1)
 
         # Parameters
         self.pars = ss.omerge({
-            'dur_inf': 10/365,  # Median for those who spontaneously clear: 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.03,
+            'dur_inf_in_days': sps.lognorm(s=0.6, scale=10),  # median of 10 days (IQR 7–15 days) https://sti.bmj.com/content/96/8/556
+            'p_symp': sps.bernoulli(p=0.5),  # Share of infections that are symptomatic. Placeholder value
+            'p_clear': sps.bernoulli(p=0.2),  # Share of infections that spontaneously clear: https://sti.bmj.com/content/96/8/556
+            'seed_infections': sps.bernoulli(p=0.1),
         }, self.pars)
 
         # Additional states dependent on parameter values, e.g. self.p_symp?
@@ -44,11 +45,11 @@ def init_results(self, sim):
 
     def update_results(self, sim):
         super(Gonorrhea, self).update_results(sim)
-        self.results['n_sympotmatic'][sim.ti] = np.count_nonzero(self.symptomatic)
+        self.results['n_symptomatic'][sim.ti] = np.count_nonzero(self.symptomatic)
         self.results['new_clearances'][sim.ti] = np.count_nonzero(self.ti_clearance == sim.ti)
         return
 
-    def update_states_pre(self, sim):
+    def update_pre(self, sim):
         # What if something in here should depend on another module?
         # I guess we could just check for it e.g., 'if HIV in sim.modules' or
         # 'if 'hiv' in sim.people' or something
@@ -81,18 +82,15 @@ def set_prognoses(self, sim, uids, from_uids=None):
         self.ti_infected[uids] = sim.ti
 
         # Set infection status
-        n_symptomatic = int(self.pars.p_symp * len(uids))
-        symptomatic_uids = np.random.choice(uids, n_symptomatic, replace=False)
-        self.symptomatic[symptomatic_uids] = True
+        symp_uids = self.pars.p_symp.filter(uids)
+        self.symptomatic[symp_uids] = True
 
         # Set natural clearance
-        n_clear = int(self.pars.p_clear * len(uids))
-        clear_uids = np.random.choice(uids, n_clear, replace=False)
-        dur = sim.ti + np.random.poisson(self.pars['dur_inf']/sim.pars.dt, len(clear_uids))
+        clear_uids = self.pars.p_clear.filter(uids)
+        dur = sim.ti + self.pars['dur_inf_in_days'].rvs(clear_uids)/365/sim.pars.dt # Convert from days to years and then adjust for dt
         self.ti_clearance[clear_uids] = dur
 
         return
 
     def set_congenital(self, sim, uids):
         pass
-
diff --git a/stisim/diseases/hiv.py b/stisim/diseases/hiv.py
index 9811eb79..0b45edc8 100644
--- a/stisim/diseases/hiv.py
+++ b/stisim/diseases/hiv.py
@@ -6,6 +6,7 @@
 import sciris as sc
 import stisim as ss
 from .disease import STI
+import scipy.stats as sps
 
 __all__ = ['HIV', 'ART', 'CD4_analyzer']
 
@@ -15,34 +16,54 @@ class HIV(STI):
     def __init__(self, pars=None):
         super().__init__(pars)
 
-        # States additional to the default disease states (see base class)
-        self.on_art = ss.State('on_art', bool, False)
-        self.cd4 = ss.State('cd4', float, 500)
+        self.on_art      = ss.State('on_art', bool, False)
+        self.ti_art      = ss.State('ti_art', int, ss.INT_NAN)
+        self.cd4         = ss.State('cd4', float, 500)
+        self.ti_dead     = ss.State('ti_dead', int, ss.INT_NAN) # Time of HIV-cause death
+
+        self.death_prob_per_dt = sps.bernoulli(p=self.death_prob)
 
         self.pars = ss.omerge({
             'cd4_min': 100,
             'cd4_max': 500,
             'cd4_rate': 5,
-            'init_prev': 0.05,
+            'seed_infections': sps.bernoulli(p=0.05),
             'eff_condoms': 0.7,
+            'art_efficacy': 0.96,
         }, self.pars)
 
         return
 
-    def update_states_pre(self, sim):
+    @staticmethod
+    def death_prob(self, sim, uids):
+        return 0.05 / (self.pars.cd4_min - self.pars.cd4_max)**2 *  (self.cd4[uids] - self.pars.cd4_max)**2
+
+    def update_pre(self, sim):
         """ Update CD4 """
         self.cd4[sim.people.alive & self.infected & self.on_art] += (self.pars.cd4_max - self.cd4[sim.people.alive & self.infected & self.on_art])/self.pars.cd4_rate
         self.cd4[sim.people.alive & self.infected & ~self.on_art] += (self.pars.cd4_min - self.cd4[sim.people.alive & self.infected & ~self.on_art])/self.pars.cd4_rate
+
+        self.rel_trans[sim.people.alive & self.infected & self.on_art] = 1 - self.pars['art_efficacy']
+
+        can_die = ss.true(sim.people.alive & sim.people.hiv.infected)
+        hiv_deaths = self.death_prob_per_dt.filter(can_die)
+        
+        sim.people.request_death(hiv_deaths)
+        self.ti_dead[hiv_deaths] = sim.ti
         return
 
     def init_results(self, sim):
         """
         Initialize results
         """
-        return super().init_results(sim)
+        super().init_results(sim)
+        self.results += ss.Result(self.name, 'new_deaths', sim.npts, dtype=int)
+        return
 
     def update_results(self, sim):
-        return super(HIV, self).update_results(sim)
+        super(HIV, self).update_results(sim)
+        self.results['new_deaths'][sim.ti] = np.count_nonzero(self.ti_dead == sim.ti)
+        return 
 
     def make_new_cases(self, sim):
         # eff_condoms = sim.pars[self.name]['eff_condoms'] # TODO figure out how to add this
@@ -65,41 +86,40 @@ def set_congenital(self, sim, uids):
 
 class ART(ss.Intervention):
 
-    def __init__(self, t: np.array, capacity: np.array):
+    def __init__(self, t: np.array, coverage: np.array, **kwargs):
         self.requires = HIV
         self.t = sc.promotetoarray(t)
-        self.capacity = sc.promotetoarray(capacity)
+        self.coverage = sc.promotetoarray(coverage)
+
+        super().__init__(**kwargs)
+
+        self.prob_art_at_infection = sps.bernoulli(p=lambda self, sim, uids: np.interp(sim.year, self.t, self.coverage))
         return
 
     def initialize(self, sim):
-        sim.hiv.results += ss.Result(self.name, 'n_art', sim.npts, dtype=int)
+        super().initialize(sim)
+        sim.results.hiv += ss.Result(self.name, 'n_art', sim.npts, dtype=int)
+        self.initialized = True
         return
 
     def apply(self, sim):
-        if sim.t < self.t[0]:
+        if sim.ti < self.t[0]:
             return
 
-        capacity = self.capacity[np.where(self.t <= sim.t)[0][-1]]
-        on_art = sim.people.alive & sim.people.hiv.on_art
-
-        n_change = capacity - np.count_nonzero(on_art)
-        if n_change > 0:
-            # Add more ART
-            eligible = sim.people.alive & sim.people.hiv.infected & ~sim.people.hiv.on_art
-            n_eligible = np.count_nonzero(eligible)
-            if n_eligible:
-                inds = np.random.choice(ss.true(eligible), min(n_eligible, n_change), replace=False)
-                sim.people.hiv.on_art[inds] = True
-        elif n_change < 0:
-            # Take some people off ART
-            eligible = sim.people.alive & sim.people.hiv.infected & sim.people.hiv.on_art
-            inds = np.random.choice(ss.true(eligible), min(n_change), replace=False)
-            sim.people.hiv.on_art[inds] = False
+        ti_delay = 1 # 1 time step delay TODO
+        recently_infected = ss.true((sim.people.hiv.ti_infected == sim.ti-ti_delay) & sim.people.alive)
+
+        n_added = 0
+        if len(recently_infected) > 0:
+            inds = self.prob_art_at_infection.filter(recently_infected)
+            sim.people.hiv.on_art[inds] = True
+            sim.people.hiv.ti_art[inds] = sim.ti
+            n_added = len(inds)
 
         # Add result
         sim.results.hiv.n_art = np.count_nonzero(sim.people.alive & sim.people.hiv.on_art)
 
-        return
+        return n_added
 
 
 #%% Analyzers
diff --git a/stisim/distributions.py b/stisim/distributions.py
index e5fd0a1b..235be0f6 100644
--- a/stisim/distributions.py
+++ b/stisim/distributions.py
@@ -1,268 +1,248 @@
 """ 
-stisim distributions
-
-Example usage
-
->>> dist = stisim.normal(1,1) # Make a distribution
->>> dist()  # Draw a sample
->>> dist(10) # Draw several samples
->>> dist.sample(10) # Same as above
->>> stisim.State('foo', float, fill_value=dist)  # Use distribution as the fill value for a state
->>> disease.pars['immunity'] = dist  # Store the distribution as a parameter
->>> disease.pars['immunity'].sample(5)  # Draw some samples from the parameter
->>> stisim.poisson(rate=1).sample(n=10)  # Sample from a temporary distribution
+Distribution support extending scipy with two key functionalities:
+1. Callable parameters
+2. Ability to use MultiRNG for common random number support
 """
 
 import numpy as np
-import sciris as sc
-from .utils import get_subclasses
-
-__all__ = [
-    'Distribution', 'uniform', 'choice', 'normal', 'normal_pos', 'normal_int', 'lognormal', 'lognormal_int',
-    'poisson', 'neg_binomial', 'beta', 'gamma', 'data_dist',
-]
-
-
-
-
-class Distribution():
-
-    def __init__(self):
+from stisim.utils import INT_NAN
+from stisim.random import SingleRNG, MultiRNG
+from stisim import options, int_
+from copy import deepcopy
+
+__all__ = ['ScipyDistribution']
+
+from scipy.stats._discrete_distns import bernoulli_gen
+
+
+class ScipyDistribution():
+    def __init__(self, gen, rng=None):
+        self.gen = None
+        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(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 not isinstance(size, (int, np.int64, 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, int_]:
+                    if not options.multirng:
+                        n_samples = len(size)
+                    else:
+                        v = size.__array__() # TODO - check if this works without calling __array__()?
+                        try:
+                            slots = self.random_state.slots[v].__array__()
+                            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 == 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__().copy() # Save full pars for handling later
+                                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
+                vals = super().rvs(*args, **kwargs)
+                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[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 = super().rvs(*args, **kwargs) # Draw again for slot repeat
+                        #assert np.allclose(slots[cur_inds], repeat_slot_u) # TEMP: Check alignment
+                        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
+            
+
+
+        self.rng = self.set_rng(rng, gen)
+        self.gen = starsim_gen(name=gen.dist.name, seed=self.rng)(**gen.kwds)
         return
 
-    def mean(self):
-        """
-        Return the mean value
-        """
-        raise NotImplementedError
-
-    def __call__(self, n=1, **kwargs):
-        return self.sample(n, **kwargs)
-
-    @property
-    def name(self):
-        return self.__class__.__name__
-
-    def sample(self, n=1, **kwargs):
-        """
-        Return a specified number of samples from the distribution
-        """
-        raise NotImplementedError
-
-    @classmethod
-    def create(cls, name, *args, **kwargs):
-        """
-        Create a distribution instance by name
-
-        :param name: A string with the name of a distribution class e.g., 'normal'
-        :param args:
-        :param kwargs:
-        :return:
-        """
-        for subcls in get_subclasses(cls):
-            if subcls.__name__ == name:
-                return subcls(*args, **kwargs)
-        else:
-            raise KeyError(f'Distribution "{name}" did not match any known distributions')
-
-class delta(Distribution):
-    """
-    Delta function at specified value
-    """
-
-    def __init__(self, value):
-        self.value = value
-
-    def sample(self, n=1):
-        return np.full(n, fill_value=self.value)
-
-class data_dist(Distribution):
-    """ Sample from data """
-
-    def __init__(self, vals, bins):
-        self.vals = vals
-        self.bins = bins
-
-    def mean(self):
+    @staticmethod
+    def set_rng(rng, gen):
+        # Handle random generators
+        ret = gen.random_state # Default
+        if options.multirng and rng and (gen.random_state == np.random.mtrand._rand):
+            # 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)}')
+        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
 
-    def sample(self, n=1):
-        """ Sample using CDF """
-        bin_midpoints = self.bins[:-1] + np.diff(self.bins) / 2
-        cdf = np.cumsum(self.vals)
-        cdf = cdf / cdf[-1]
-        values = np.random.rand(n)
-        value_bins = np.searchsorted(cdf, values)
-        return bin_midpoints[value_bins]
-
-
-class uniform(Distribution):
-    """
-    Uniform distribution
-    """
-
-    def __init__(self, low, high):
-        self.low = low
-        self.high = high
-
-    def mean(self):
-        return (self.low + self.high) / 2
-
-    def sample(self, n=1):
-        return np.random.uniform(low=self.low, high=self.high, size=n)
-
-
-class choice(Distribution):
-    """
-    Choose from samples, optionally with specific probabilities
-    """
-
-    def __init__(self, choices, probabilities=None, replace=True):
-        self.choices = choices
-        self.probabilities = probabilities
-        self.replace = replace
-
-    def sample(self, n, replace=True):
-        return np.random.choice(a=self.choices, p=self.probabilities, replace=self.replace, size=n)
-
-
-class normal(Distribution):
-    """
-    Normal distribution
-    """
-
-    def __init__(self, mean, std):
-        self.mean = mean
-        self.std = std
-
-    def sample(self, n=1):
-        return np.random.normal(loc=self.mean, scale=self.std, size=n)
-
-
-class normal_pos(normal):
-    """
-    right-sided normal (i.e. only +ve values), with mean=par1, std=par2 of the underlying normal
-
-    WARNING - this function came from hpvsim but confirm that the implementation is correct?
-    """
-
-    def sample(self, n=1):
-        return np.abs(super().sample(n))
-
-
-class normal_int(normal):
-    """
-    Normal distribution returning only integer values
-    """
-
-    def sample(self, n=1):
-        return np.round(super().sample(n))
-
-
-class lognormal(Distribution):
-    """
-    lognormal with mean=par1, std=par2 (parameters are for the lognormal, not the underlying normal)
-
-    Lognormal distributions are parameterized with reference to the underlying normal distribution (see:
-    https://docs.scipy.org/doc/numpy-1.14.0/reference/generated/numpy.random.lognormal.html), but this
-    function assumes the user wants to specify the mean and std of the lognormal distribution.
-    """
-
-    def __init__(self, mean, std):
-        self.mean = mean
-        self.std = std
-        self.underlying_mean = np.log(mean ** 2 / np.sqrt(std ** 2 + mean ** 2))  # Computes the mean of the underlying normal distribution
-        self.underlying_std = np.sqrt(np.log(std ** 2 / mean ** 2 + 1))  # Computes sigma for the underlying normal distribution
-
-    def sample(self, n=1):
-
-        if (sc.isnumber(self.mean) and self.mean > 0) or (sc.checktype(self.mean, 'arraylike') and (self.mean > 0).all()):
-            return np.random.lognormal(mean=self.underlying_mean, sigma=self.underlying_std, size=n)
-        else:
-            return np.zeros(n)
-
-
-class lognormal_int(lognormal):
-    """
-    Lognormal returning only integer values
-    """
-
-    def sample(self, n=1):
-        return np.round(super().sample(n))
-
-
-class poisson(Distribution):
-    """
-    Poisson distribution
-    """
-
-    def __init__(self, rate):
+    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))
+        return result
+    
+    def __getattr__(self, attr):
+        # Returns wrapped generator.(attr) if not a property
+        if attr == '__getstate__':
+            # Must be from pickle, return a callable function that returns None
+            return lambda: None
+        #elif attr == '__deepcopy__':
+        #    return self
+        elif attr in ['__setstate__', '__await__']:
+            # Must be from pickle, async programming, copy
+            return None
+        try:
+            return self.__getattribute__(attr)
+        except Exception:
+            try:
+                return getattr(self.gen, attr) # .dist?
+            except Exception as e:
+                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)]
+
+'''
+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 mean(self):
-        return self.rate
-
-    def sample(self, n=1):
-        return np.random.poisson(self.rate, n)
-
-
-class neg_binomial(Distribution):
-    """
-    Negative binomial distribution
-
-    Negative binomial distributions are parameterized with reference to the mean and dispersion parameter k
-    (see: https://en.wikipedia.org/wiki/Negative_binomial_distribution). The r parameter of the underlying
-    distribution is then calculated from the desired mean and k. For a small mean (~1), a dispersion parameter
-    of ∞ corresponds to the variance and standard deviation being equal to the mean (i.e., Poisson). For a
-    large mean (e.g. >100), a dispersion parameter of 1 corresponds to the standard deviation being equal to
-    the mean.
-    """
-
-    def __init__(self, mean, dispersion, step=1):
-        """
-        mean (float): the rate of the process (same as Poisson)
-        dispersion (float):  dispersion parameter; lower is more dispersion, i.e. 0 = infinite, ∞ = Poisson
-        n (int): number of trials
-        step (float): the step size to use if non-integer outputs are desired
-        """
-        self.mean = mean
-        self.dispersion = dispersion
-        self.step = step
-
-    def sample(self, n=1):
-        nbn_n = self.dispersion
-        nbn_p = self.dispersion / (self.mean / self.step + self.dispersion)
-        return np.random.negative_binomial(n=nbn_n, p=nbn_p, size=n) * self.step
-
-
-class beta(Distribution):
-    """
-    Beta distribution
-    """
-
-    def __init__(self, alpha, beta):
-        self.alpha = alpha
-        self.beta = beta
-
-    def mean(self):
-        return self.alpha / (self.alpha + self.beta)
-
-    def sample(self, n=1):
-        return np.random.beta(a=self.alpha, b=self.beta, size=n)
-
-
-class gamma(Distribution):
-    """
-    Gamma distribution
-    """
-
-    def __init__(self, shape, scale):
-        self.shape = shape
-        self.scale = scale
-
-    def mean(self):
-        return self.shape * self.scale
-
-    def sample(self, n=1):
-        return np.random.gamma(shape=self.shape, scale=self.scale, size=n)
+    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
+'''
diff --git a/stisim/interventions.py b/stisim/interventions.py
index 16dbd298..26d0e86e 100644
--- a/stisim/interventions.py
+++ b/stisim/interventions.py
@@ -9,10 +9,8 @@
 
 
 class Intervention(ss.Module):
-
-
     def __call__(self, *args, **kwargs):
         return self.apply(*args, **kwargs)
 
     def apply(self, sim, *args, **kwargs):
-        raise NotImplementedError
+        raise NotImplementedError
\ No newline at end of file
diff --git a/stisim/modules.py b/stisim/modules.py
index a43e6eac..17a2d0a7 100644
--- a/stisim/modules.py
+++ b/stisim/modules.py
@@ -2,16 +2,16 @@
 Disease modules
 '''
 
-import numpy as np
 import sciris as sc
 import stisim as ss
+from scipy.stats._distn_infrastructure import rv_frozen
 
 __all__ = ['Module']
 
 
 class Module(sc.prettyobj):
 
-    def __init__(self, pars=None, name = None, label=None, requires=None, *args, **kwargs):
+    def __init__(self, pars=None, name=None, label=None, requires=None, *args, **kwargs):
         self.pars = ss.omerge(pars)
         self.name = name if name else self.__class__.__name__.lower() # Default name is the class name
         self.label = label if label else ''
@@ -19,12 +19,13 @@ def __init__(self, pars=None, name = None, label=None, requires=None, *args, **k
         self.results = ss.ndict(type=ss.Result)
         self.initialized = False
         self.finalized = False
+
         return
 
     def check_requires(self, sim):
         for req in self.requires:
-            if req not in sim.modules:
-                raise Exception(f'{self.name} (label={self.label}) requires module {req} but the Sim did not contain this module')
+            if req not in [m.__class__ for m in sim.modules]:
+                raise Exception(f'{self.name} (label={self.label}) requires module {req} but the Sim did not contain a module of this type.')
         return
 
     def initialize(self, sim):
@@ -38,9 +39,29 @@ def initialize(self, sim):
         """
         self.check_requires(sim)
 
+        # Connect the random number generators to the sim. The RNGs for this module should be initialized
+        # first as some of the module's State instances may require them to generate initial values
+        for rng in self.rngs:
+            if not rng.initialized:
+                rng.initialize(sim.rng_container, sim.people.slot)
+
+        # Initialize distributions in pars
+        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)
+            #elif isinstance(value, ss.rate):
+            #    self.pars[key].initialize(sim, f'{self.name}_{self.label}_{key}')
+
+        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)
+
         # Connect the states to the sim
+        # Will use random numbers, so do after distribution initialization
         for state in self.states:
-            state.initialize(sim.people)
+            state.initialize(sim)
 
         self.initialized = True
         return
@@ -64,3 +85,21 @@ def states(self):
         """
         return [x for x in self.__dict__.values() if isinstance(x, ss.State)]
 
+    @property
+    def rngs(self):
+        """
+        Return a flat collection of all random number generators, as with states above
+
+        :return:
+        """
+        return [x for x in self.__dict__.values() if isinstance(x, (ss.MultiRNG, ss.SingleRNG))]
+
+    @property
+    def scipy_dbns(self):
+        """
+        Return a flat collection of all ScipyDistributions
+
+        :return:
+        """
+        return [x for x in self.__dict__.values() if isinstance(x, ss.ScipyDistribution)] \
+             + [x for x in self.pars.values()     if isinstance(x, ss.ScipyDistribution)]
diff --git a/stisim/networks/__init__.py b/stisim/networks/__init__.py
index eebb18e6..ee4ed0e5 100644
--- a/stisim/networks/__init__.py
+++ b/stisim/networks/__init__.py
@@ -1,2 +1,2 @@
 from .networks import *
-from .random import *
\ No newline at end of file
+from .randnet import *
\ No newline at end of file
diff --git a/stisim/networks/networks.py b/stisim/networks/networks.py
index cd0c48ef..216de2b5 100644
--- a/stisim/networks/networks.py
+++ b/stisim/networks/networks.py
@@ -6,10 +6,16 @@
 import numpy as np
 import sciris as sc
 import stisim as ss
+import scipy.optimize as spo
+import scipy.stats as sps
+import scipy.spatial as spsp
 import pandas as pd
+from scipy.stats._distn_infrastructure import rv_frozen
+import scipy.stats as sps
+
 
 # Specify all externally visible functions this file defines
-__all__ = ['Networks', 'Network', 'NetworkConnector', 'DynamicNetwork', 'SexualNetwork', 'mf', 'msm', 'mf_msm', 'hpv_network', 'maternal', 'static']
+__all__ = ['Networks', 'Network', 'NetworkConnector', 'SexualNetwork', 'mf', 'msm', 'mf_msm', 'hpv_network', 'maternal', 'embedding', 'static']
 
 
 class Network(ss.Module):
@@ -276,7 +282,10 @@ def __init__(self, pars=None, key_dict=None):
     def end_pairs(self, people):
         dt = people.dt
         self.contacts.dur = self.contacts.dur - dt
-        active = self.contacts.dur > 0
+
+        # Non-alive agents are removed
+        active = (self.contacts.dur > 0) & people.alive[self.contacts.p1] & people.alive[self.contacts.p2]
+
         self.contacts.p1 = self.contacts.p1[active]
         self.contacts.p2 = self.contacts.p2[active]
         self.contacts.beta = self.contacts.beta[active]
@@ -294,12 +303,14 @@ def initialize(self, sim):
             nw.initialize(sim)
         for cn in self._connectors.values():
             cn.initialize(sim)
+        return
 
     def update(self, people):
         for nw in self.values():
             nw.update(people)
         for cn in self._connectors.values():
             cn.update(people)
+        return
 
 
 class SexualNetwork(Network):
@@ -308,7 +319,8 @@ def __init__(self, pars=None):
         super().__init__(pars)
 
     def active(self, people):
-        return self.participant & (people.age > self.debut)
+        # Exclude people who are not alive
+        return self.participant & (people.age > self.debut) & people.alive
 
     def available(self, people, sex):
         # Currently assumes unpartnered people are available
@@ -316,7 +328,7 @@ def available(self, people, sex):
         # This property could also be overwritten by a NetworkConnector
         # which could incorporate information about membership in other
         # contact networks
-        return np.setdiff1d(people.uid[people[sex] & self.active(people)], self.members)
+        return np.setdiff1d(ss.true(people[sex] & self.active(people)), self.members) # ss.true instead of people.uid[]?
 
 
 class mf(SexualNetwork, DynamicNetwork):
@@ -325,138 +337,85 @@ class mf(SexualNetwork, DynamicNetwork):
     relationship durations.
     """
 
-    def __init__(self, pars=None):
-        DynamicNetwork.__init__(self)
-        SexualNetwork.__init__(self, pars)
-        self.pars = ss.omerge({
-            'dur': ss.lognormal(15, 15),  # Can vary by age, year, and individual pair
-            'part_rates': 0.9,  # Participation rates - can vary by sex and year
+    def __init__(self, pars=None, key_dict=None):
+        desired_mean = 15
+        desired_std = 15
+        mu = np.log(desired_mean**2 / np.sqrt(desired_mean**2 + desired_std**2))
+        sigma = np.sqrt(np.log(1 + desired_std**2 / desired_mean**2))
+        pars = ss.omerge({
+            'duration_dist': sps.lognorm(s=sigma, scale=np.exp(mu)), # Can vary by age, year, and individual pair. Set scale=exp(mu) and s=sigma where mu,sigma are of the underlying normal distribution.
+            'participation_dist': sps.bernoulli(p=0.9),  # Probability of participating in this network - can vary by individual properties (age, sex, ...) using callable parameter values
+            'debut_dist': sps.norm(loc=16, scale=2),  # Age of debut can vary by using callable parameter values
             'rel_part_rates': 1.0,
-            'debut': 16,  # Age of debut can vary by sex, year, and individual
-            'rel_debut': 1.0,
-        }, self.pars)
-
-        # Set metadata for select parameters - this tells us which dimensions are allowed to vary
-        self.meta_pars = sc.objdict({
-            'part_rates': {'year': 2000, 'sex': ['m', 'f'], 'part_rates': None},
-            'debut': {'year': 2000, 'sex': ['m', 'f'], 'debut': None, 'std': 0, 'dist': 'lognormal'},
-            'dur': {'year': 2000, 'age': 0, 'dur': [None], 'std': 0, 'dist': 'lognormal'},
-        })
-        self.validate_pars()
-        return
-
-    def validate_pars(self):
-        """ Validate parameters and expand assumptions """
-
-        for parname, defaults in self.meta_pars.items():
-            cols = defaults.keys()
-            par = self.pars[parname]
-
-            if sc.checktype(par, pd.DataFrame):
-                if not set(cols).issubset(par.columns):
-                    errormsg = f'Please ensure the columns of the data for parameter {parname} are {cols}.'
-                    raise ValueError(errormsg)
-                df = par
-
-            elif sc.checktype(par, dict):
-                if not set(cols).issubset(par.keys()):
-                    errormsg = f'Please ensure the keys of the data dict for parameters {parname} are {cols}.'
-                    raise ValueError(errormsg)
-                df = pd.DataFrame(par)
-
-            elif sc.checktype(par, ss.Distribution):
-                df = pd.DataFrame(defaults)
-                df[parname] = par.mean
-                df['std'] = par.std
-                df['dist'] = par.name
-
-            elif sc.isnumber(par):
-                dd = {}
-                for col in cols:
-                    dd[col] = defaults[col]
-                    if col == parname:
-                        dd[col] = par
-                df = pd.DataFrame(dd)
-
-            else:
-                errormsg = f'Data type {type(par)} not understood for parameter {parname}.'
-                raise ValueError(errormsg)
-
-            self.pars[parname] = df
-
-        # Check values
-        if (self.pars.part_rates['part_rates'] < 0).any() or (self.pars.part_rates['part_rates'] > 1).any():
-            raise ValueError('Invalid participation rate for network.')
+        }, pars)
 
+        DynamicNetwork.__init__(self, key_dict)
+        SexualNetwork.__init__(self, pars)
         return
 
     def initialize(self, sim):
         super().initialize(sim)
+        
         self.set_network_states(sim.people)
         self.add_pairs(sim.people, ti=0)
+        return
 
     def set_network_states(self, people, upper_age=None):
         """ Set network states including age of entry into network and participation rates """
-        self.set_participation(people, upper_age=upper_age)
-        self.set_debut(people, upper_age=upper_age)
+        self.set_debut(people, upper_age) # Debut before participation!
+        self.set_participation(people, upper_age)
+        return
 
     def set_participation(self, people, upper_age=None):
-        # Set people who will participate in the network at some point
-        year = people.year
-        if upper_age is None: uids = people.uid
-        else: uids = people.uid[(people.age < upper_age)]
+        if upper_age is None:
+            uids = people.uid
+        else:
+            uids = people.uid[(people.age < upper_age)]
 
-        for sk in ['f', 'm']:
-            uids = ss.true(people[sk][uids])
-            df = self.pars.part_rates.loc[self.pars.part_rates.sex == sk]
-            pr = np.interp(year, df['year'], df['part_rates']) * self.pars.rel_part_rates
-            dist = ss.choice([True, False], probabilities=[pr, 1-pr])(len(uids))
-            self.participant[uids] = dist
+        # Set people who will participate in the network at some point
+        self.participant[uids] = self.pars.participation_dist.rvs(uids)
+        return
 
     def set_debut(self, people, upper_age=None):
         # Set debut age
-        if upper_age is None: uids = people.uid
-        else: uids = people.uid[(people.age < upper_age)]
-
-        for sk in ['f', 'm']:
-            uids = ss.true(people[sk][uids])
-            df = self.pars.debut.loc[self.pars.debut.sex == sk]
-            all_years = self.pars.debut.year.values
-            year_ind = sc.findnearest(all_years, people.year)
-            nearest_year = all_years[year_ind]
-
-            mean = np.interp(people.year, df['year'], df['debut'])
-            std = np.interp(people.year, df['year'], df['std'])
-            dist = df.loc[df.year == nearest_year].dist.iloc[0]
-            debut_vals = ss.Distribution.create(dist, mean, std)(len(uids)) * self.pars.rel_debut
-            self.debut[uids] = debut_vals
+        if upper_age is None:
+            uids = people.uid
+        else:
+            uids = people.uid[(people.age < upper_age)]
+
+        self.debut[uids] = self.pars.debut_dist.rvs(uids)
+        uids_to_update = uids[np.isnan(people.debut[uids])]
+        people.debut[uids_to_update] = self.debut[uids_to_update]
+        return
 
     def add_pairs(self, people, ti=None):
         available_m = self.available(people, 'male')
         available_f = self.available(people, 'female')
+
+        # random.choice is not common-random-number safe, and therefore we do
+        # not try to Stream-ify the following draws at this time.
         if len(available_m) <= len(available_f):
             p1 = available_m
-            p2 = np.random.choice(available_f, len(p1), replace=False)
+            p2 = np.random.choice(a=available_f, size=len(p1), replace=False)
         else:
             p2 = available_f
-            p1 = np.random.choice(available_m, len(p2), replace=False)
+            p1 = np.random.choice(a=available_m, size=len(p2), replace=False)
 
         beta = np.ones_like(p1)
 
         # Figure out durations
-        all_years = self.pars.dur.year.values
-        year_ind = sc.findnearest(all_years, people.year)
-        nearest_year = all_years[year_ind]
-        mean = np.interp(people.year, self.pars.dur['year'], self.pars.dur['dur'])
-        std = np.interp(people.year, self.pars.dur['year'], self.pars.dur['std'])
-        dur_dist = self.pars.dur.loc[self.pars.dur.year == nearest_year].dist.iloc[0]
-        dur_vals = ss.Distribution.create(dur_dist, mean, std)(len(p1))
+        # print('DJK TODO')
+        if ss.options.multirng and (len(p1) == len(np.unique(p1))):
+            # No duplicates and user has enabled multirng, so use slotting based on p1
+            dur_vals = self.pars.duration_dist.rvs(p1)
+        else:
+            dur_vals = self.pars.duration_dist.rvs(len(p1)) # Just use len(p1) to say how many draws are needed
 
         self.contacts.p1 = np.concatenate([self.contacts.p1, p1])
         self.contacts.p2 = np.concatenate([self.contacts.p2, p2])
         self.contacts.beta = np.concatenate([self.contacts.beta, beta])
         self.contacts.dur = np.concatenate([self.contacts.dur, dur_vals])
-        return
+        return len(p1)
 
     def update(self, people, dt=None):
         self.end_pairs(people)
@@ -471,15 +430,28 @@ class msm(SexualNetwork, DynamicNetwork):
     """
 
     def __init__(self, pars=None):
-        DynamicNetwork.__init__(self)
-        SexualNetwork.__init__(self, pars)
-        self.pars = ss.omerge({
-            'dur': ss.lognormal(5, 3),
+        default_pars = {
             'part_rates': 0.1,  # Participation rates - can vary by sex and year
             'rel_part_rates': 1.0,
-            'debut': ss.lognormal(18, 2),
-            'rel_debut': 1.0,
-        }, self.pars)
+        }
+
+        desired_mean = 5
+        desired_std = 3
+        mu = np.log(desired_mean**2 / np.sqrt(desired_mean**2 + desired_std**2))
+        sigma = np.sqrt(np.log(1 + desired_std**2 / desired_mean**2))
+        default_pars['duration_dist'] = sps.lognorm(s=sigma, scale=np.exp(mu)) # Can vary by age, year, and individual pair. Set scale=exp(mu) and s=sigma where mu,sigma are of the underlying normal distribution.
+
+        desired_mean = 18
+        desired_std = 2
+        mu = np.log(desired_mean**2 / np.sqrt(desired_mean**2 + desired_std**2))
+        sigma = np.sqrt(np.log(1 + desired_std**2 / desired_mean**2))
+        default_pars['debut_dist'] = sps.lognorm(s=sigma, scale=np.exp(mu))
+
+        pars = ss.omerge(default_pars, pars)
+        DynamicNetwork.__init__(self)
+        SexualNetwork.__init__(self, pars)
+
+        return
 
     def initialize(self, sim):
         # Add more here in line with MF network, e.g. age of debut
@@ -488,6 +460,7 @@ def initialize(self, sim):
         super().initialize(sim)
         self.set_network_states(sim.people)
         self.add_pairs(sim.people, ti=0)
+        return
 
     def set_network_states(self, people, upper_age=None):
         """ Set network states including age of entry into network and participation rates """
@@ -497,11 +470,14 @@ def set_network_states(self, people, upper_age=None):
         # Participation
         self.participant[people.female] = False
         pr = self.pars.part_rates
-        dist = ss.choice([True, False], probabilities=[pr, 1-pr])(len(uids))
+        dist = sps.bernoulli.rvs(p=pr, size=len(uids))
         self.participant[uids] = dist
 
         # Debut
-        self.debut[uids] = self.pars.debut(len(uids))
+        self.debut[uids] = self.pars.debut_dist.rvs(len(uids)) # Just pass len(uids) as this network is not crn safe anyway
+        uids_to_update = uids[np.isnan(people.debut[uids])]
+        people.debut[uids_to_update] = self.debut[uids_to_update]
+        return
 
     def add_pairs(self, people, ti=None):
         # Pair all unpartnered MSM
@@ -509,15 +485,85 @@ def add_pairs(self, people, ti=None):
         n_pairs = int(len(available_m)/2)
         p1 = available_m[:n_pairs]
         p2 = available_m[n_pairs:n_pairs*2]
+
+        # Figure out durations
+        # print('DJK TODO')
+        if ss.options.multirng and (len(p1) == len(np.unique(p1))):
+            # No duplicates and user has enabled multirng, so use slotting based on p1
+            dur = self.pars['duration_dist'].rvs(p1)
+        else:
+            dur = self.pars['duration_dist'].rvs(len(p1)) # Just use len(p1) to say how many draws are needed
+
         self.contacts.p1 = np.concatenate([self.contacts.p1, p1])
         self.contacts.p2 = np.concatenate([self.contacts.p2, p2])
         self.contacts.beta = np.concatenate([self.contacts.beta, np.ones_like(p1)])
-        self.contacts.dur = np.concatenate([self.contacts.dur, self.pars.dur.sample(len(p1))])
+        self.contacts.dur = np.concatenate([self.contacts.dur, dur])
+        return len(p1)
 
     def update(self, people, dt=None):
         self.end_pairs(people)
         self.set_network_states(people, upper_age=people.dt)
         self.add_pairs(people)
+        return
+
+
+class embedding(mf):
+    """
+    Heterosexual age-assortative network based on a one-dimensional embedding. Could be made more generic.
+    """
+
+    def __init__(self, pars=None, **kwargs):
+        """
+        Create a sexual network from a 1D embedding based on age
+
+        male_shift is the average age that males are older than females in partnerships
+        std is the standard deviation of noise added to the age of each individual when seeking a pair. Larger values will created more diversity in age gaps.
+        
+        """
+        pars = ss.omerge({
+            'embedding_func': sps.norm(loc=self.embedding_loc, scale=2),
+        }, pars)
+        super().__init__(pars, **kwargs)
+        return
+
+    @staticmethod
+    def embedding_loc(self, sim, uids):
+        loc = sim.people.age[uids].values
+        loc[sim.people.female[uids]] += 5 # Shift females so they will be paired with older men
+        return loc
+
+    def add_pairs(self, people, ti=None):
+        available_m = self.available(people, 'male')
+        available_f = self.available(people, 'female')
+
+        if not len(available_m) or not len(available_f):
+            if ss.options.verbose > 1:
+                print('No pairs to add')
+            return 0
+
+        available = np.concatenate((available_m, available_f))
+        loc = self.pars['embedding_func'].rvs(available)
+        loc_f = loc[people.female[available]]
+        loc_m = loc[~people.female[available]]
+
+        dist_mat = spsp.distance_matrix(loc_m[:, np.newaxis], loc_f[:, np.newaxis])
+
+        ind_m, ind_f = spo.linear_sum_assignment(dist_mat)
+        # loc_f[ind_f[0]] is close to loc_m[ind_m[0]]
+
+        n_pairs = len(ind_f)
+
+        beta = np.ones(n_pairs)
+
+        # Figure out durations
+        p1 = available_m[ind_m]
+        dur_vals = self.pars['duration_dist'].rvs(p1)
+
+        self.contacts.p1 = np.concatenate([self.contacts.p1, p1])
+        self.contacts.p2 = np.concatenate([self.contacts.p2, available_f[ind_f]])
+        self.contacts.beta = np.concatenate([self.contacts.beta, beta])
+        self.contacts.dur = np.concatenate([self.contacts.dur, dur_vals])
+        return len(beta)
 
 
 class NetworkConnector(ss.Module):
@@ -525,19 +571,10 @@ class NetworkConnector(ss.Module):
     Template for a connector between networks.
     """
     def __init__(self, *args, networks=None, pars=None, **kwargs):
-        super().__init__(pars, *args, **kwargs)
-        self.networks = networks
-        self.pars = ss.omerge(pars)
-
-    def initialize(self, sim):
-        # Check that the requested networks are in the sim
-        avail = set(sim.people.networks.keys())
-        if not set(self.networks).issubset(avail):
-            errormsg = f'Connection between {self.networks} has been requested, but available networks are {avail}.'
-            raise ValueError(errormsg)
+        super().__init__(pars, requires=networks, *args, **kwargs)
         return
 
-    def set_participation(self, people):
+    def set_participation(self, people, upper_age=None):
         pass
 
     def update(self, people):
@@ -546,12 +583,15 @@ def update(self, people):
 
 class mf_msm(NetworkConnector):
     """ Combines the MF and MSM networks """
-    def __init__(self, networks=None, pars=None):
-        super().__init__(networks=networks, pars=pars)
-        self.pars = ss.omerge({
+    def __init__(self, pars=None):
+        networks = [ss.mf, ss.msm]
+        pars = ss.omerge({
             'prop_bi': 0.5,  # Could vary over time -- but not by age or sex or individual
-        }, self.pars)
-        if self.networks is None: self.networks = ['mf', 'msm']
+        }, pars)
+        super().__init__(networks=networks, pars=pars)
+
+        self.bi_dist = sps.bernoulli(p=self.pars.prop_bi)
+        return
 
     def initialize(self, sim):
         super().initialize(sim)
@@ -574,28 +614,32 @@ def set_participation(self, people, upper_age=None):
         # Male participation rate uses info about cross-network participation.
         # First, we determine who's participating in the MSM network
         pr = msm.pars.part_rates
-        dist = ss.choice([True, False], probabilities=[pr, 1-pr])(len(uids))
+        dist = sps.bernoulli.rvs(p=pr, size=len(uids))
         msm.participant[uids] = dist
 
         # Now we take the MSM participants and determine which are also in the MF network
         msm_uids = ss.true(msm.participant[uids])  # Males in the MSM network
-        bi_uids = ss.binomial_filter(self.pars.prop_bi, msm_uids)  # Males in both MSM and MF networks
+        bi_uids = self.bi_dist.filter(msm_uids)  # Males in both MSM and MF networks
         mf_excl_set = np.setdiff1d(uids, msm_uids)  # Set of males who aren't in the MSM network
 
         # What remaining share to we need?
         mf_df = mf.pars.part_rates.loc[mf.pars.part_rates.sex == 'm']  # Male participation in the MF network
         mf_pr = np.interp(people.year, mf_df['year'], mf_df['part_rates']) * mf.pars.rel_part_rates
         remaining_pr = max(mf_pr*len(uids)-len(bi_uids), 0)/len(mf_excl_set)
-        mf_excl_uids = ss.binomial_filter(remaining_pr, mf_excl_set)  # Males in MF network only
+
+        # Don't love the following new syntax:
+        mf_excl_uids = mf_excl_set[sps.uniform.rvs(size=len(mf_excl_set)) < remaining_pr]
+
         mf.participant[bi_uids] = True
         mf.participant[mf_excl_uids] = True
+        return
 
     def update(self, people):
         self.set_participation(people, upper_age=people.dt)
         return
 
 
-class hpv_network(SexualNetwork, DynamicNetwork):
+class hpv_network(mf):
     def __init__(self, pars=None):
 
         key_dict = {
@@ -603,27 +647,56 @@ def __init__(self, pars=None):
             'start': ss.float_,
         }
 
-        DynamicNetwork.__init__(self, key_dict=key_dict)
-        SexualNetwork.__init__(self, pars=pars)
-
         # Define default parameters
-        self.pars = dict()
-        self.pars['cross_layer'] = 0.05  # Proportion of agents who have concurrent cross-layer relationships
-        self.pars['partners'] = ss.poisson(rate=0.01)  # The number of concurrent sexual partners
-        self.pars['acts'] = ss.neg_binomial(mean=80, dispersion=40)  # The number of sexual acts per year
-        self.pars['age_act_pars'] = dict(peak=30, retirement=100, debut_ratio=0.5,
+        default_pars = dict()
+        default_pars['cross_layer']   = 0.05  # Proportion of agents who have concurrent cross-layer relationships
+        default_pars['partner_dist']  = sps.poisson(mu=0.01)  # The number of concurrent sexual partners
+
+        # TODO: Wrap so user can provide mean and dispersion directly - see #168
+        mu = 80 # Mean
+        alpha = 40 # Dispersion
+        sigma2 = mu + alpha * mu**2
+        n = mu**2 / (sigma2 - mu)
+        p = mu / sigma2
+        default_pars['act_dist']      = sps.nbinom(n=n, p=p)  # The number of sexual acts per year
+
+        default_pars['age_act_pars']  = dict(peak=30, retirement=100, debut_ratio=0.5,
                                          retirement_ratio=0.1)  # Parameters describing changes in coital frequency over agent lifespans
-        self.pars['condoms'] = 0.2  # The proportion of acts in which condoms are used
-        self.pars['dur_pship'] = ss.normal_pos(mean=1, std=1)  # Duration of partnerships
-        self.pars['participation'] = None  # Incidence of partnership formation by age
-        self.pars['mixing'] = None  # Mixing matrices for storing age differences in partnerships
+        default_pars['condoms']       = 0.2  # The proportion of acts in which condoms are used
+
+        low = 0
+        loc = 1
+        scale = 1
+        a = (low - loc) / scale
+        default_pars['duration_dist'] = sps.truncnorm(a=a, b=np.inf, loc=loc, scale=scale)  # Duration of partnerships
+
+        #default_pars['participation'] = None  # Incidence of partnership formation by age
+        default_pars['mixing']        = None  # Mixing matrices for storing age differences in partnerships
+
+        self.agebins = [0, 5, 10, 15, 20, 25, 30, 35, 40, 45, 50, 55, 60, 65, 70, 75]
+        # Share of females of each age newly having casual relationships
+        self.f_participation = [0, 0, 0.10, 0.7, 0.8, 0.6, 0.6, 0.4, 0.1, 0.05, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001]
+        # Share of males of each age newly having casual relationships
+        self.m_participation = [0, 0, 0.05, 0.7, 0.8, 0.6, 0.6, 0.4, 0.4, 0.3, 0.1, 0.05, 0.01, 0.01, 0.001, 0.001]
+
+        default_pars['participation_dist'] = sps.bernoulli(p=self.participation)
 
-        self.update_pars(pars)
+        pars = ss.omerge(default_pars, pars)
+        super().__init__(pars, key_dict)
+ 
         self.get_layer_probs()
+        self.validate_pars()
+
         return
 
-    def initialize(self, people):
-        return self.add_pairs(people, ti=0)
+    def initialize(self, sim):
+        super().initialize(sim)
+        return self.add_pairs(sim.people, ti=0)
+
+    def validate_pars(self):
+        for par in ['partner_dist', 'act_dist', 'duration_dist', 'participation_dist']:
+            if not isinstance(self.pars[par], rv_frozen):
+                raise Exception(f'Network parameter {par} must be an instance of a scipy frozen distribution')
 
     def update_pars(self, pars):
         if pars is not None:
@@ -631,6 +704,14 @@ def update_pars(self, pars):
                 self.pars[k] = v
         return
 
+    @staticmethod
+    def participation(self, sim, uids):
+        p = np.ones_like(uids, dtype=ss.float_)
+        fem = sim.people.female[uids]
+        p[fem] = np.interp(sim.people.age[uids[fem]], self.agebins, self.f_participation)
+        p[~fem] = np.interp(sim.people.age[uids[~fem]], self.agebins, self.m_participation)
+        return p
+
     def get_layer_probs(self):
 
         defaults = {}
@@ -654,16 +735,7 @@ def get_layer_probs(self):
             [75, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1],
         ])
 
-        participation = np.array([
-            [0, 5, 10, 15, 20, 25, 30, 35, 40, 45, 50, 55, 60, 65, 70, 75],
-            [0, 0, 0.10, 0.7, 0.8, 0.6, 0.6, 0.4, 0.1, 0.05, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001],
-            # Share of females of each age newly having casual relationships
-            [0, 0, 0.05, 0.7, 0.8, 0.6, 0.6, 0.4, 0.4, 0.3, 0.1, 0.05, 0.01, 0.01, 0.001, 0.001]],
-            # Share of males of each age newly having casual relationships
-        )
-
         defaults['mixing'] = mixing
-        defaults['participation'] = participation
 
         for pkey, pval in defaults.items():
             if self.pars[pkey] is None:
@@ -671,12 +743,11 @@ def get_layer_probs(self):
 
         return
 
-    def add_pairs(self, people, ti=0):
-
-        female = people.female
-        active = people.active
-        f_active = female & active
-        m_active = ~female & active
+    def set_participation(self, people, upper_age=None):
+        if upper_age is None:
+            uids = people.uid
+        else:
+            uids = people.uid[(people.age < upper_age)]
 
         # Compute number of partners
         f_partnered_inds, f_partnered_counts = np.unique(self.contacts.p1, return_counts=True)
@@ -684,67 +755,53 @@ def add_pairs(self, people, ti=0):
         current_partners = np.zeros((len(people)))
         current_partners[f_partnered_inds] = f_partnered_counts
         current_partners[m_partnered_inds] = m_partnered_counts
-        partners = self.pars['partners'].sample(len(people)) + 1
+        partners = self.pars.partner_dist.rvs(len(people)) + 1
         underpartnered = current_partners < partners  # Indices of underpartnered people
-        f_eligible = f_active & underpartnered
-        m_eligible = m_active & underpartnered
-
-        # Bin the agents by age
-        bins = self.pars['participation'][0, :]  # Extract age bins
-
-        # Try randomly select females for pairing
-        f_eligible_inds = ss.true(f_eligible)  # Inds of all eligible females
-        age_bins_f = np.digitize(people.age[f_eligible_inds], bins=bins) - 1  # Age bins of selected females
-        bin_range_f = np.unique(age_bins_f)  # Range of bins
-        f = []  # Initialize the female partners
-        for ab in bin_range_f:  # Loop over age bins
-            these_f_contacts = ss.binomial_filter(self.pars['participation'][1][ab], f_eligible_inds[
-                age_bins_f == ab])  # Select females according to their participation rate in this layer
-            f += these_f_contacts.tolist()
-
-        # Select males according to their participation rate in this layer
-        m_eligible_inds = ss.true(m_eligible)
-        age_bins_m = np.digitize(people.age[m_eligible_inds], bins=bins) - 1
-        bin_range_m = np.unique(age_bins_m)  # Range of bins
-        m = []  # Initialize the male partners
-        for ab in bin_range_m:
-            these_m_contacts = ss.binomial_filter(self.pars['participation'][2][ab], m_eligible_inds[
-                age_bins_m == ab])  # Select males according to their participation rate in this layer
-            m += these_m_contacts.tolist()
+
+        # Set people who will participate in the network at some point
+        can_participate = ss.true(people.active * underpartnered)
+        self.participant[uids] = self.pars.participation_dist.rvs(can_participate)
+        return
+
+    def add_pairs(self, people, ti=0):
+        participating = ss.true(self.participant) # Will be the same people each time, with participation decided once per person
+        f = participating[people.female[participating]]
+        m = participating[~people.female[participating]]
 
         # Create preference matrix between eligible females and males that combines age and geo mixing
         age_bins_f = np.digitize(people.age[f],
-                                 bins=bins) - 1  # Age bins of females that are entering new relationships
-        age_bins_m = np.digitize(people.age[m], bins=bins) - 1  # Age bins of active and participating males
+                                 bins=self.agebins) - 1  # Age bins of females that are entering new relationships
+        age_bins_m = np.digitize(people.age[m], bins=self.agebins) - 1  # Age bins of active and participating males
         age_f, age_m = np.meshgrid(age_bins_f, age_bins_m)
         pair_probs = self.pars['mixing'][age_m, age_f + 1]
 
         f_to_remove = pair_probs.max(axis=0) == 0  # list of female inds to remove if no male partners are found for her
-        f = [i for i, flag in zip(f, f_to_remove) if ~flag]  # remove the inds who don't get paired on this timestep
+        #f = [i for i, flag in zip(f, f_to_remove) if ~flag]  # remove the inds who don't get paired on this timestep
+        f = f[~f_to_remove]
         selected_males = []
         if len(f):
             pair_probs = pair_probs[:, np.invert(f_to_remove)]
             choices = []
             fems = np.arange(len(f))
             f_paired_bools = np.full(len(fems), True, dtype=bool)
-            np.random.shuffle(fems)
+            np.random.shuffle(fems) # TODO: Stream-ify?
             for fem in fems:
                 m_col = pair_probs[:, fem]
                 if m_col.sum() > 0:
                     m_col_norm = m_col / m_col.sum()
-                    choice = np.random.choice(len(m_col_norm), 1, replace=False, p=m_col_norm)
+                    choice = np.random.choice(len(m_col_norm), p=m_col_norm) # TODO: Stream-ify?
                     choices.append(choice)
                     pair_probs[choice, :] = 0  # Once male partner is assigned, remove from eligible pool
                 else:
                     f_paired_bools[fem] = False
-            selected_males = np.array(m)[np.array(choices).flatten()]
-            f = np.array(f)[f_paired_bools]
+            selected_males = m[np.array(choices).flatten()]
+            f = f[f_paired_bools]
 
-        p1 = np.array(f)
+        p1 = f
         p2 = selected_males
         n_partnerships = len(p1)
-        dur = self.pars['dur_pship'].sample(n_partnerships)
-        acts = self.pars['acts'].sample(n_partnerships)
+        dur = self.pars.duration_dist.rvs(n_partnerships)
+        acts = self.pars.act_dist.rvs(n_partnerships)
         age_p1 = people.age[p1]
         age_p2 = people.age[p2]
 
@@ -792,7 +849,7 @@ def add_pairs(self, people, ti=0):
             beta=beta
         )
         self.append(new_contacts)
-        return
+        return len(new_contacts)
 
     def update(self, people, ti=None, dt=None):
         if ti is None: ti = people.ti
@@ -809,12 +866,12 @@ def update(self, people, ti=None, dt=None):
 
 
 class maternal(Network):
-    def __init__(self, key_dict=None, vertical=True):
+    def __init__(self, key_dict=None, vertical=True, **kwargs):
         """
         Initialized empty and filled with pregnancies throughout the simulation
         """
         key_dict = sc.mergedicts({'dur': ss.float_}, key_dict)
-        super().__init__(key_dict=key_dict, vertical=vertical)
+        super().__init__(key_dict=key_dict, vertical=vertical, **kwargs)
         return
 
     def update(self, people, dt=None):
@@ -839,7 +896,7 @@ def add_pairs(self, mother_inds, unborn_inds, dur):
         self.contacts.p2 = np.concatenate([self.contacts.p2, unborn_inds])
         self.contacts.beta = np.concatenate([self.contacts.beta, beta])
         self.contacts.dur = np.concatenate([self.contacts.dur, dur])
-        return
+        return len(mother_inds)
 
 class static(Network):
     """
@@ -892,4 +949,3 @@ def get_contacts(self):
         self.contacts.beta = np.concatenate([self.contacts.beta, np.ones_like(p1s)])
         return
 
-
diff --git a/stisim/networks/random.py b/stisim/networks/randnet.py
similarity index 72%
rename from stisim/networks/random.py
rename to stisim/networks/randnet.py
index 89132f9d..67cb07ed 100644
--- a/stisim/networks/random.py
+++ b/stisim/networks/randnet.py
@@ -4,22 +4,30 @@
 import numpy as np
 import numba as nb
 from .networks import Network
+from scipy.stats._distn_infrastructure import rv_frozen
+from typing import Union
 
 __all__ = ['RandomNetwork']
 
 class RandomNetwork(Network):
 
-    def __init__(self, *, n_contacts: ss.Distribution, dynamic=True, **kwargs):
+    def __init__(self, *, n_contacts : Union[int, rv_frozen], dynamic=True, **kwargs):
         """
-        :param n_contacts: A distribution of contacts e.g., ss.delta(5), ss.neg_binomial(5,2)
+        :param n_contacts: A SciPy distribution from which the number of
+        contacts are sampled e.g., sps.nbinom(n=5, p=2). Note distribution
+        parameter values can be callable for greater functionality.
         :param dynamic: If True, regenerate contacts each timestep
         """
         super().__init__(**kwargs)
-        self.n_contacts = n_contacts
+        if isinstance(n_contacts, rv_frozen):
+            self.n_contacts = ss.ScipyDistribution(n_contacts, f'{self.__class__.__name__}_{self.name}_{self.label}')
+
         self.dynamic = dynamic
 
     def initialize(self, sim):
         super().initialize(sim)
+        if isinstance(self.n_contacts, ss.ScipyDistribution):
+            self.n_contacts.initialize(sim, self)
         self.update(sim.people, force=True)
 
     @staticmethod
@@ -45,8 +53,6 @@ def get_contacts(inds, number_of_contacts):
             of random contacts.
 
         Returns: Two arrays, for source and target
-
-
         """
 
         total_number_of_half_edges = np.sum(number_of_contacts)
@@ -71,7 +77,11 @@ def update(self, people: ss.People, force: bool = True) -> None:
         if not self.dynamic and not force:
             return
 
-        number_of_contacts = self.n_contacts.sample(len(people))
+        if isinstance(self.n_contacts, ss.ScipyDistribution):
+            number_of_contacts = self.n_contacts.rvs(people.alive) # or people.uid?
+        else:
+            number_of_contacts = np.full(len(people), self.n_contacts)
+
         number_of_contacts = np.round(number_of_contacts / 2).astype(ss.int_)  # One-way contacts
         self.contacts.p1, self.contacts.p2 = self.get_contacts(people.uid.__array__(), number_of_contacts)
         self.contacts.beta = np.ones(len(self.contacts.p1), dtype=ss.float_)
\ No newline at end of file
diff --git a/stisim/parameters.py b/stisim/parameters.py
index 6b7f37bd..1cd66ed2 100644
--- a/stisim/parameters.py
+++ b/stisim/parameters.py
@@ -46,6 +46,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.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)
 
         # Events and interventions
@@ -55,12 +56,13 @@ def __init__(self, **kwargs):
 
         # Network parameters, generally initialized after the population has been constructed
         self.networks        = sc.autolist()  # Network types and parameters
-        self.debut           = dict(f=ss.normal(mean=15.0, std=2.0),
-                                    m=ss.normal(mean=17.5, std=2.0))
 
         # Update with any supplied parameter values and generate things that need to be generated
         self.update(kwargs)
 
+        if self.slot_scale < 1:
+            raise Exception('The value of the "slot_scale" parameter must be a number >= 1.0')
+
         return
 
     def update_pars(self, pars=None, create=False, **kwargs):
@@ -86,4 +88,4 @@ def update_pars(self, pars=None, create=False, **kwargs):
 
 
 def make_pars(**kwargs):
-    return Parameters(**kwargs)
+    return Parameters(**kwargs)
\ No newline at end of file
diff --git a/stisim/people.py b/stisim/people.py
index 2b2b792f..593549f7 100644
--- a/stisim/people.py
+++ b/stisim/people.py
@@ -7,6 +7,7 @@
 import pandas as pd
 import sciris as sc
 import stisim as ss
+from scipy.stats import bernoulli, uniform
 
 __all__ = ['BasePeople', 'People']
 
@@ -24,36 +25,58 @@ def __init__(self, n):
         self._uid_map = ss.DynamicView(int, fill_value=ss.INT_NAN)  # This variable tracks all UIDs ever created
         self.uid = ss.DynamicView(int, fill_value=ss.INT_NAN)  # This variable tracks all UIDs currently in use
 
-        self._uid_map.initialize(n)
+        self._uid_map.grow(n)
         self._uid_map[:] = np.arange(0, n)
-        self.uid.initialize(n)
+        self.uid.grow(n)
         self.uid[:] = np.arange(0, n)
 
+        # A slot is a special state managed internally by BasePeople
+        # This is because it needs to be updated separately from any other states, as other states
+        # might have fill_values that depend on the slot
+        self.slot = ss.State('slot', int, ss.INT_NAN)
+
         self.ti = None  # Track simulation time index
         self.dt = np.nan  # Track simulation time step
 
-        # We internally store states in a dict keyed by the memory ID of the state, so that we can have colliding names
+        # User-facing collection of states
+        self.states = ss.ndict(type=ss.State)
+
+        # We also internally store states in a dict keyed by the memory ID of the state, so that we can have colliding names
         # e.g., across modules, but we will never change the size of a State multiple times in the same iteration over
         # _states. This is a hidden variable because it is internally used to synchronize the size of all States contained
         # within the sim, regardless of where they are. In contrast, `People.states` offers a more user-friendly way to access
         # a selection of the states e.g., module states could be added in there, while intervention states might not
         self._states = {}
 
+    @property
+    def rngs(self):
+        return [x for x in self.__dict__.values() if isinstance(x, (ss.MultiRNG, ss.SingleRNG))]
+
     def __len__(self):
         """ Length of people """
         return len(self.uid)
 
-    def add_state(self, state):
+    def add_state(self, state, die=True):
         if id(state) not in self._states:
             self._states[id(state)] = state
+            self.states.append(state)  # Expose these states with their original names
+            setattr(self, state.name, state)
+        elif die:
+            errormsg = f'Cannot add state {state} since already added'
+            raise ValueError(errormsg)
         return
 
-    def grow(self, n):
+    def grow(self, n, new_slots=None):
         """
         Increase the number of agents
 
         :param n: Integer number of agents to add
+        :param new_slots: Optionally specify the slots to assign for the new agents. Otherwise, it will default to the new UIDs
         """
+
+        if n == 0:
+            return np.array([], dtype=ss.int_)
+
         start_uid = len(self._uid_map)
         start_idx = len(self.uid)
 
@@ -66,8 +89,12 @@ def grow(self, n):
         self.uid.grow(n)
         self.uid[new_inds] = new_uids
 
+        # We need to grow the slots as well
+        self.slot.grow(new_uids)
+        self.slot[new_uids] = new_slots if new_slots is not None else new_uids
+
         for state in self._states.values():
-            state.grow(n)
+            state.grow(new_uids)
 
         return new_uids
 
@@ -77,13 +104,14 @@ def remove(self, uids_to_remove):
 
         :param uids_to_remove: An int/list/array containing the UID(s) to remove
         """
+
         # Calculate the *indices* to keep
         keep_uids = self.uid[~np.in1d(self.uid, uids_to_remove)]  # Calculate UIDs to keep
         keep_inds = self._uid_map[keep_uids]  # Calculate indices to keep
 
         # Trim the UIDs and states
         self.uid._trim(keep_inds)
-        for state in self._states.values():
+        for state in self._states.values(): # includes self.slot
             state._trim(keep_inds)
 
         # Update the UID map
@@ -138,9 +166,7 @@ class People(BasePeople):
         ppl = ss.People(2000)
     """
 
-    # %% Basic methods
-
-    def __init__(self, n, age_data=None, extra_states=None, networks=None):
+    def __init__(self, n, age_data=None, extra_states=None, networks=None, rand_seed=0):
         """ Initialize """
 
         super().__init__(n)
@@ -148,44 +174,57 @@ def __init__(self, n, age_data=None, extra_states=None, networks=None):
         self.initialized = False
         self.version = ss.__version__  # Store version info
 
+        # Handle states
         states = [
-            ss.State('age', float, 0),
-            ss.State('female', bool, ss.choice([True, False])),
+            ss.State('age', float, np.nan), # NaN until conceived
+            ss.State('female', bool, bernoulli(p=0.5)),
             ss.State('debut', float),
-            ss.State('alive', bool, True),
             ss.State('ti_dead', int, ss.INT_NAN),  # Time index for death
+            ss.State('alive', bool, True),  # Time index for death
             ss.State('scale', float, 1.0),
         ]
         states.extend(sc.promotetolist(extra_states))
-
-        self.states = ss.ndict()
-        self._initialize_states(states)
+        for state in states:
+            self.add_state(state)
+        self._initialize_states(sim=None) # No sim yet, but initialize what we can
+        
         self.networks = ss.Networks(networks)
 
         # Set initial age distribution - likely move this somewhere else later
-        age_data_dist = self.validate_age_data(age_data)
-        self.age[:] = age_data_dist.sample(len(self))
+        self.age_data_dist = self.get_age_dist(age_data)
 
         return
 
     @staticmethod
-    def validate_age_data(age_data):
-        """ Validate age data """
-        if age_data is None: return ss.uniform(0, 100)
+    def get_age_dist(age_data):
+        """ Return an age distribution based on provided data """
+        if age_data is None:
+            return uniform(loc=0, scale=100) # low and width
         if sc.checktype(age_data, pd.DataFrame):
             return ss.data_dist(vals=age_data['value'].values, bins=age_data['age'].values)
 
-    def initialize(self):
-        """ Initialization - TBC what needs to go here """
-        self.initialized = True
+    def _initialize_states(self, sim=None):
+        for state in self.states.values():
+            state.initialize(sim=sim, people=self)  # Connect the state to this people instance
         return
 
-    def _initialize_states(self, states):
-        for state in states:
-            self.add_state(state)  # Register the state internally for dynamic growth
-            self.states.append(state)  # Expose these states with their original names
-            state.initialize(self)  # Connect the state to this people instance
-            setattr(self, state.name, state)
+    def initialize(self, sim):
+        """ Initialization """
+    
+        # For People initialization, first initialize slots, then initialize RNGs, then initialize remaining states
+        # This is because some states may depend on RNGs being initialized to generate initial values
+        self.slot.initialize(sim)
+        self.slot[:] = self.uid
+    
+        # Initialize all RNGs (noting that includes those that are declared in child classes)
+        for rng in self.rngs:
+            rng.initialize(sim.rng_container, self.slot)
+            
+        # Define age (CK: why is age handled differently than sex?)
+        self._initialize_states(sim=sim) # Now initialize with the sim
+        self.age[:] = self.age_data_dist.rvs(len(self))
+        
+        self.initialized = True
         return
 
     def add_module(self, module, force=False):
@@ -227,9 +266,8 @@ def remove_dead(self, sim):
         Remove dead agents
         """
         uids_to_remove = ss.true(self.dead)
-        self.remove(uids_to_remove)
-        for network in self.networks.values():
-            network.remove_uids(uids_to_remove)
+        if len(uids_to_remove):
+            self.remove(uids_to_remove)
         return
 
     def update_post(self, sim):
diff --git a/stisim/random.py b/stisim/random.py
new file mode 100644
index 00000000..930395d8
--- /dev/null
+++ b/stisim/random.py
@@ -0,0 +1,260 @@
+import hashlib
+import numpy as np
+import stisim as ss
+
+__all__ = ['RNGContainer', 'MultiRNG', 'SingleRNG', 'RNG']
+
+
+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 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.multirng:
+        rng = MultiRNG(*args, **kwargs)
+    else:
+        rng = SingleRNG(*args, **kwargs)
+    return rng
+
+
+class MultiRNG(np.random.Generator):
+    """
+    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 stisim 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), an np.arange 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"
+        """
+
+        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.seed = None # Will be determined once added to the RNG Container
+        self.initialized = False
+        self.ready = True
+        return
+
+    def initialize(self, container, slots):
+        if self.initialized:
+            return
+
+        if container is not None:
+            self.seed = container.add(self) # base_seed + seed_offset
+        else:
+            # Enable use of MultiRNG without a container
+            self.seed = self.seed_offset
+
+        if isinstance(slots, int):
+            # Handle edge case in which the user wants n sequential slots, as used in testing.
+            self.slots = np.arange(slots)
+        else:
+            self.slots = slots # E.g. sim.people.slots (instead of using uid as the slots directly)
+
+        if 'bit_generator' not in self.kwargs:
+            self.kwargs['bit_generator'] = np.random.PCG64DXSM(seed=self.seed)
+        super().__init__(**self.kwargs)
+
+        self._init_state = self.bit_generator.state # Store the initial state
+
+        self.initialized = True
+        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
+
+
+class SingleRNG(np.random.Generator):
+    """
+    Dummy class that mirrors MultiRNG but provides the numpy singleton rng instead
+    """
+    
+    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:
+            container.add(self) # base_seed + seed_offset
+
+        if 'bit_generator' not in self.kwargs:
+            self.kwargs['bit_generator'] = np.random.mtrand._rand._bit_generator
+        super().__init__(**self.kwargs)
+
+        self.initialized = True
+        return
+
+    def reset(self):
+        pass
+
+    def step(self, ti):
+        pass
diff --git a/stisim/results.py b/stisim/results.py
index b81173a5..3716781f 100644
--- a/stisim/results.py
+++ b/stisim/results.py
@@ -25,7 +25,8 @@ def __repr__(self):
         return out
 
     def __array_finalize__(self, obj):
-        if obj is None: return
+        if obj is None:
+            return
         self.name = getattr(obj, 'name', None)
         self.module = getattr(obj, 'module', None)
         return
diff --git a/stisim/settings.py b/stisim/settings.py
index f7d221d1..dd60282f 100644
--- a/stisim/settings.py
+++ b/stisim/settings.py
@@ -56,6 +56,12 @@ def get_orig_options():
         optdesc = sc.objdict()  # Help for the options
         options = sc.objdict()  # The options
 
+        optdesc.multirng = 'Set True to use multiple random number generators\
+            via the MultiRNG class, which is compatible with common random\
+            numbers. The default value is False, which instead uses a single\
+            centralized random number generator.'
+        options.multirng = False
+
         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('STISIM_VERBOSE', 0.1))
 
diff --git a/stisim/sim.py b/stisim/sim.py
index c2288fea..64ed911b 100644
--- a/stisim/sim.py
+++ b/stisim/sim.py
@@ -7,9 +7,15 @@
 import sciris as sc
 import stisim as ss
 import itertools
+import numba as nb
 
 __all__ = ['Sim', 'AlreadyRunError', 'diff_sims']
 
+@nb.njit
+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)
 
 class Sim(sc.prettyobj):
 
@@ -43,6 +49,9 @@ def __init__(self, pars=None, label=None, people=None, demographics=None, diseas
         self.interventions = ss.ndict(type=ss.Intervention)
         self.analyzers = ss.ndict(type=ss.Analyzer)
 
+        # Initialize the random number generator container
+        self.rng_container = ss.RNGContainer()
+
         return
 
     @property
@@ -75,8 +84,10 @@ def initialize(self, popdict=None, 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.init_people(popdict=popdict, reset=reset, **kwargs)  # Create all the people (the heaviest step)
         self.init_networks()
         self.init_demographics()
@@ -90,7 +101,7 @@ def initialize(self, popdict=None, reset=False, **kwargs):
 
         # Reset the random seed to the default run seed, so that if the simulation is run with
         # reset_seed=False right after initialization, it will still produce the same output
-        ss.set_seed(self.pars['rand_seed'] + 1)
+        ss.set_seed(self.pars['rand_seed'] + 1) # Hopefully not used now that we can use multiple random number generators
 
         # Final steps
         self.initialized = True
@@ -121,17 +132,15 @@ def validate_pars(self):
         Some parameters can take multiple types; this makes them consistent.
         """
         # Handle n_agents
-        if self.pars['n_agents'] is not None:
+        if self.people is not None:
+            self.pars['n_agents'] = len(self.people)
+        #elif self.popdict is not None: # STIsim does not currenlty support self.popdict
+            #self.pars['n_agents'] = len(self.popdict)
+        elif self.pars['n_agents'] is not None:
             self.pars['n_agents'] = int(self.pars['n_agents'])
         else:
-            if self.people is not None:
-                self.pars['n_agents'] = len(self.people)
-            else:
-                if self.popdict is not None:
-                    self.pars['n_agents'] = len(self.popdict)
-                else:
-                    errormsg = 'Must supply n_agents, a people object, or a popdict'
-                    raise ValueError(errormsg)
+            errormsg = 'Must supply n_agents, a people object, or a popdict'
+            raise ValueError(errormsg)
 
         # Handle end and n_years
         if self.pars['end']:
@@ -183,27 +192,24 @@ def init_people(self, popdict=None, reset=False, verbose=None, **kwargs):
             verbose = self.pars['verbose']
         if verbose > 0:
             resetstr = ''
-            if self.people:
-                resetstr = ' (resetting people)' if reset else ' (warning: not resetting sim.people)'
+            if self.people and reset:
+                resetstr = ' (resetting people)'
             print(f'Initializing sim{resetstr} with {self.pars["n_agents"]:0n} agents')
 
         # If people have not been supplied, make them
-        if self.people is None:
+        if self.people is None or reset:
             self.people = ss.People(n=self.pars['n_agents'], **kwargs)  # This just assigns UIDs and length
 
         # If a popdict has not been supplied, we can make one from location data
         if popdict is None:
             if self.pars['location'] is not None:
                 # Check where to get total_pop from
-                if self.pars[
-                    'total_pop'] is not None:  # If no pop_scale has been provided, try to get it from the location
+                if self.pars['total_pop'] is not None:  # If no pop_scale has been provided, try to get it from the location
                     errormsg = 'You can either define total_pop explicitly or via the location, but not both'
                     raise ValueError(errormsg)
                 total_pop, popdict = ss.make_popdict(n=self.pars['n_agents'], location=self.pars['location'], verbose=self.pars['verbose'])
-
             else:
-                if self.pars[
-                    'total_pop'] is not None:  # If no pop_scale has been provided, try to get it from the location
+                if self.pars['total_pop'] is not None:  # If no pop_scale has been provided, try to get it from the location
                     total_pop = self.pars['total_pop']
                 else:
                     if self.pars['pop_scale'] is not None:
@@ -217,7 +223,7 @@ def init_people(self, popdict=None, reset=False, verbose=None, **kwargs):
 
         # Any other initialization
         if not self.people.initialized:
-            self.people.initialize()
+            self.people.initialize(self)
 
         # Set time attributes
         self.people.ti = self.ti
@@ -293,6 +299,9 @@ def init_interventions(self):
                 errormsg = f'Intervention {intervention} does not seem to be a valid intervention: must be a function or Intervention subclass'
                 raise TypeError(errormsg)
 
+            for rng in intervention.rngs:
+                rng.initialize(self.rng_container, self.people.slot)
+
         return
 
     def init_analyzers(self):
@@ -331,6 +340,9 @@ def step(self):
         if self.complete:
             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
+
         # Clean up dead agents, if removing agents is enabled
         if self.pars.remove_dead:
             self.people.remove_dead(self)
@@ -449,7 +461,7 @@ def finalize(self, verbose=None):
 
         for module in self.modules:
             module.finalize(self)
-            
+
         self.summarize()
         self.results_ready = True  # Set this first so self.summary() knows to print the results
         self.ti -= 1  # During the run, this keeps track of the next step; restore this be the final day of the sim
@@ -459,7 +471,7 @@ def summarize(self):
         summary = sc.objdict()
         flat = sc.flattendict(self.results, sep='_')
         for k,v in flat.items():
-            summary[k] = np.array(v).mean() # To fix later
+            summary[k] = v.mean()
         self.summary = summary
         return summary
         
@@ -779,4 +791,4 @@ def diff_sims(sim1, sim2, skip_key_diffs=False, skip=None, full=False, output=Fa
     else:
         if not output:
             print('Sims match')
-    return
\ No newline at end of file
+    return
diff --git a/stisim/states.py b/stisim/states.py
index af8ea5f8..75100089 100644
--- a/stisim/states.py
+++ b/stisim/states.py
@@ -3,7 +3,10 @@
 import sciris as sc
 import numba as nb
 from stisim.utils import INT_NAN
+from stisim.distributions import ScipyDistribution
+from stisim.utils import warn
 from numpy.lib.mixins import NDArrayOperatorsMixin  # Inherit from this to automatically gain operators like +, -, ==, <, etc.
+from scipy.stats._distn_infrastructure import rv_frozen
 
 __all__ = ['State', 'DynamicView']
 
@@ -46,7 +49,11 @@ def dtype(self):
     @nb.njit
     def _get_vals_uids(vals, key, uid_map):
         """
-        Extract valued from a collection of UIDs
+        Extract values from a collection of UIDs
+
+        This function is used to retrieve values based on UID. As indexing a FusedArray returns a new FusedArray,
+        this method also populates the new UID map for use in the subsequently created FusedArray, avoiding the
+        need to re-compute it separately.
 
         :param vals: A 1D np.ndarray containing the values
         :param key: A 1D np.ndnarray of integers containing the UIDs to query
@@ -76,10 +83,15 @@ def _set_vals_uids_multiple(vals, key, uid_map, value):
         :param value: A 1D np.ndarray the same length as ``key`` containing values to insert
         :return:
         """
+
         for i in range(len(key)):
+            if key[i] >= len(uid_map):
+                raise IndexError('UID not present in array (requested UID is larger than the maximum UID in use)')
             idx = uid_map[key[i]]
             if idx == INT_NAN:
                 raise IndexError('UID not present in array')
+            elif idx >= len(vals):
+                raise Exception(f'Attempted to write to a non-existant index - this can happen if attempting to write to new entries that have not yet been allocated via grow()')
             vals[idx] = value[i]
 
     @staticmethod
@@ -95,9 +107,13 @@ def _set_vals_uids_single(vals, key, uid_map, value):
         :return:
         """
         for i in range(len(key)):
+            if key[i] >= len(uid_map):
+                raise IndexError('UID not present in array (requested UID is larger than the maximum UID in use)')
             idx = uid_map[key[i]]
             if idx == INT_NAN:
                 raise IndexError('UID not present in array')
+            elif idx >= len(vals):
+                raise Exception('Attempted to write to a non-existant index - this can happen if attempting to write to new entries that have not yet been allocated via grow()')
             vals[idx] = value
 
     def __getitem__(self, key):
@@ -105,7 +121,7 @@ def __getitem__(self, key):
             if isinstance(key, (int, np.integer)):
                 # Handle getting a single item by UID
                 return self.values[self._uid_map[key]]
-            elif isinstance(key, (np.ndarray, FusedArray)):
+            elif isinstance(key, (np.ndarray, FusedArray, DynamicView)):
                 if key.dtype.kind == 'b':
                     # Handle accessing items with a logical array. Surprisingly, it seems faster to use nonzero() to convert
                     # it to indices first. Also, the pure Python implementation is difficult to improve upon using numba
@@ -116,7 +132,7 @@ def __getitem__(self, key):
                     values = self.values[mapped_key]
                 else:
                     # Access items by an array of integers. We do get a decent performance boost from using numba here
-                    values, uids, new_uid_map = self._get_vals_uids(self.values, key, self._uid_map.__array__())
+                    values, uids, new_uid_map = self._get_vals_uids(self.values, key.__array__(), self._uid_map.__array__())
             elif isinstance(key, slice):
                 if key.start is None and key.stop is None and key.step is None:
                     return sc.dcp(self)
@@ -235,20 +251,20 @@ def __array_wrap__(self, out_arr, context=None):
 
 
 class DynamicView(NDArrayOperatorsMixin):
-    def __init__(self, dtype, fill_value=0):
+    def __init__(self, dtype, fill_value=None):
         """
         Args:
             name: name of the result as used in the model
             dtype: datatype
-            fill_value: default value for this state upon model initialization. If callable, it is called with a single argument for the number of samples
+            fill_value: default value for this state upon model initialization. If not provided, it will use the default value for the dtype
             shape: If not none, set to match a string in `pars` containing the dimensionality
             label: text used to construct labels for the result for displaying on plots and other outputs
         """
-        self.dtype = dtype
-        self.fill_value = fill_value
+        self.fill_value = fill_value if fill_value is not None else dtype()
         self.n = 0  # Number of agents currently in use
-        self._data = None  # The underlying memory array (length at least equal to n)
+        self._data = np.empty(0, dtype=dtype)  # The underlying memory array (length at least equal to n)
         self._view = None  # The view corresponding to what is actually accessible (length equal to n)
+        self._map_arrays()
         return
 
     @property
@@ -256,6 +272,15 @@ def _s(self):
         # 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
+        # DynamicView(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 DynamicView'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 DynamicView._dtype() and use
+        return self._data.dtype
+
     def __len__(self):
         # Return the number of active elements
         return self.n
@@ -264,26 +289,11 @@ def __repr__(self):
         # Print out the numpy view directly
         return self._view.__repr__()
 
-    def _new_items(self, n):
-        # Create new arrays of the correct dtype and fill them based on the (optionally callable) fill value
-        if callable(self.fill_value):
-            new = np.empty(n, dtype=self.dtype)
-            new[:] = self.fill_value(n)
-        else:
-            new = np.full(n, dtype=self.dtype, fill_value=self.fill_value)
-        return new
-
-    def initialize(self, n):
-        self._data = self._new_items(n)
-        self.n = n
-        self._map_arrays()
-
     def grow(self, n):
+        # If the total number of agents exceeds the array size, extend the underlying arrays
         if self.n + n > self._s:
-            # If the total number of agents exceeds the array size, extend the storage array
             n_new = max(n, int(self._s / 2))  # Minimum 50% growth
-            self._data = np.concatenate([self._data, self._new_items(n_new)], axis=0)
-
+            self._data = np.concatenate([self._data, np.full(n_new, dtype=self.dtype, fill_value=self.fill_value)], axis=0)
         self.n += n  # Increase the count of the number of agents by `n` (the requested number of new agents)
         self._map_arrays()
 
@@ -292,7 +302,7 @@ def _trim(self, inds):
         # Note that these are indices, not UIDs!
         n = len(inds)
         self._data[:n] = self._data[inds]
-        self._data[n:self.n] = self.fill_value(self.n-n) if callable(self.fill_value) else self.fill_value
+        self._data[n:self.n] = self.fill_value
         self.n = n
         self._map_arrays()
 
@@ -326,9 +336,23 @@ def __array_ufunc__(self, *args, **kwargs):
 
 class State(FusedArray):
 
-    def __init__(self, name, dtype, fill_value=0, label=None):
+    def __init__(self, name, dtype, fill_value=None, label=None):
+        """
+
+        :param name: A string name for the state
+        :param dtype: The dtype to use for this instance
+        :param fill_value: 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
+        :param label:
+        """
+
         super().__init__(values=None, uid=None, uid_map=None)  # Call the FusedArray constructor
-        self._data = DynamicView(dtype=dtype, fill_value=fill_value)
+
+        self.fill_value = fill_value
+
+        self._data = DynamicView(dtype=dtype)
         self.name = name
         self.label = label or name
         self.values = self._data._view
@@ -340,22 +364,57 @@ def __repr__(self):
         else:
             return FusedArray.__repr__(self)
 
-    def initialize(self, people):
+    def _new_vals(self, uids):
+        if isinstance(self.fill_value, ScipyDistribution):
+            new_vals = self.fill_value.rvs(uids)
+        elif callable(self.fill_value):
+            new_vals = self.fill_value(len(uids))
+        else:
+            new_vals = self.fill_value
+        return new_vals
+
+    def initialize(self, sim=None, people=None):
         if self._initialized:
             return
 
-        people.add_state(self)
-        self._uid_map = people._uid_map
-        self.uid = people.uid
-        self._data.initialize(len(self.uid))
-        self.values = self._data._view
-        self._initialized = True
+        if sim is not None and people is None:
+            people = sim.people
 
-    def grow(self, n):
+        sim_still_needed = False
+        if isinstance(self.fill_value, rv_frozen):
+            if sim is not None:
+                self.fill_value = ScipyDistribution(self.fill_value, f'{self.__class__.__name__}_{self.label}')
+                self.fill_value.initialize(sim, self)
+            else:
+                sim_still_needed = True
+
+        people.add_state(self, die=False) # CK: should not be needed
+        if not sim_still_needed:
+            self._uid_map = people._uid_map
+            self.uid = people.uid
+            self._data.grow(len(self.uid))
+            self._data[:len(self.uid)] = self._new_vals(self.uid)
+            self.values = self._data._view
+            self._initialized = True
+        return
+
+    def grow(self, uids):
+        """
+        Add state for new agents
+
+        This method is normally only called via `People.grow()`.
+
+        :param uids: Numpy array of UIDs for the new agents being added This array should have length n
+        """
+
+        n = len(uids)
         self._data.grow(n)
         self.values = self._data._view
+        self._data[-n:] = self._new_vals(uids)
+        return
 
     def _trim(self, inds):
+        # Trim arrays to remove agents - should only be called via `People.remove()`
         self._data._trim(inds)
         self.values = self._data._view
-
+        return
diff --git a/stisim/utils.py b/stisim/utils.py
index 04081f4d..25036d70 100644
--- a/stisim/utils.py
+++ b/stisim/utils.py
@@ -34,9 +34,9 @@ class ndict(sc.objdict):
 
     **Examples**::
 
-        networks = ss.ndict(ss.simple_sexual(), ss.maternal())
-        networks = ss.ndict([ss.simple_sexual(), ss.maternal()])
-        networks = ss.ndict({'simple_sexual':ss.simple_sexual(), 'maternal':ss.maternal()})
+        networks = ss.ndict(ss.mf(), ss.maternal())
+        networks = ss.ndict([ss.mf(), ss.maternal()])
+        networks = ss.ndict({'mf':ss.mf(), 'maternal':ss.maternal()})
 
     """
 
@@ -216,100 +216,6 @@ def set_seed_regular(seed):
     return
 
 
-# %% Probabilities -- mostly not jitted since performance gain is minimal
-
-__all__ += ['binomial_arr', 'binomial_filter', 'n_poisson', 'n_neg_binomial']
-
-
-def binomial_arr(prob_arr):
-    '''
-    Binomial (Bernoulli) trials each with different probabilities.
-
-    Args:
-        prob_arr (array): array of probabilities
-
-    Returns:
-         Boolean array of which trials on the input array succeeded
-
-    **Example**::
-
-        outcomes = ss.binomial_arr([0.1, 0.1, 0.2, 0.2, 0.8, 0.8]) # Perform 6 trials with different probabilities
-    '''
-    return np.random.random(prob_arr.shape) < prob_arr
-
-
-def binomial_filter(prob, arr):
-    """
-    Binomial "filter" -- the same as n_binomial, except return
-    the elements of arr that succeeded.
-
-    Args:
-        prob (float): probability of each trial succeeding
-        arr (array): the array to be filtered
-
-    Returns:
-        Subset of array for which trials succeeded
-
-    **Example**::
-
-        inds = ss.binomial_filter(0.5, np.arange(20)**2) # Which values in the (arbitrary) array passed the coin flip
-    """
-    return arr[(np.random.random(len(arr)) < prob).nonzero()[0]]
-
-
-def binomial_arr(prob_arr):
-    """
-    Binomial (Bernoulli) trials each with different probabilities.
-
-    Args:
-        prob_arr (array): array of probabilities
-
-    Returns:
-         Boolean array of which trials on the input array succeeded
-
-    **Example**::
-
-        outcomes = ss.binomial_arr([0.1, 0.1, 0.2, 0.2, 0.8, 0.8]) # Perform 6 trials with different probabilities
-    """
-    return np.random.random(prob_arr.shape) < prob_arr
-
-
-def n_poisson(rate, n):
-    """
-    An array of Poisson trials.
-
-    Args:
-        rate (float): the rate of the Poisson process (mean)
-        n (int): number of trials
-
-    **Example**::
-
-        outcomes = ss.n_poisson(100, 20) # 20 Poisson trials with mean 100
-    """
-    return np.random.poisson(rate, n)
-
-
-def n_neg_binomial(rate, dispersion, n, step=1):  # Numba not used due to incompatible implementation
-    """
-    An array of negative binomial trials. See ss.sample() for more explanation.
-
-    Args:
-        rate (float): the rate of the process (mean, same as Poisson)
-        dispersion (float):  dispersion parameter; lower is more dispersion, i.e. 0 = infinite, ∞ = Poisson
-        n (int): number of trials
-        step (float): the step size to use if non-integer outputs are desired
-
-    **Example**::
-
-        outcomes = ss.n_neg_binomial(100, 1, 50) # 50 negative binomial trials with mean 100 and dispersion roughly equal to mean (large-mean limit)
-        outcomes = ss.n_neg_binomial(1, 100, 20) # 20 negative binomial trials with mean 1 and dispersion still roughly equal to mean (approximately Poisson)
-    """
-    nbn_n = dispersion
-    nbn_p = dispersion / (rate / step + dispersion)
-    samples = np.random.negative_binomial(n=nbn_n, p=nbn_p, size=n) * step
-    return samples
-
-
 # %% Simple array operations
 
 __all__ += ['true', 'false', 'defined', 'undefined']
@@ -396,4 +302,4 @@ def undefined(arr):
 
         inds = ss.defined(np.array([1,np.nan,0,np.nan,1,0,1]))
     """
-    return np.isnan(arr).nonzero()[-1]
+    return np.isnan(arr).nonzero()[-1]
\ No newline at end of file
diff --git a/stisim/version.py b/stisim/version.py
index 0683d0a0..04511fc5 100644
--- a/stisim/version.py
+++ b/stisim/version.py
@@ -4,6 +4,6 @@
 
 __all__ = ['__version__', '__versiondate__', '__license__']
 
-__version__ = '0.0.8'
-__versiondate__ = '2023-10-04'
+__version__ = '0.1.0'
+__versiondate__ = '2023-12-10'
 __license__ = f'STIsim {__version__} ({__versiondate__}) — © 2023 by IDM'
diff --git a/vignettes/v09_custom_networks.py b/vignettes/v09_custom_networks.py
index 88f2a309..8166efd2 100644
--- a/vignettes/v09_custom_networks.py
+++ b/vignettes/v09_custom_networks.py
@@ -58,9 +58,9 @@ def update(self, sim):
 
 # Case 2: adding a custom network
 ppl = ss.People(100)
-nwk = ss.simple_sexual()
+nwk = ss.mf()
 syph = ss.syphilis()
-sim['syphilis']['beta']['simple_sexual'] *= 2  # Increase beta for syph transmission
+sim['syphilis']['beta']['mf'] *= 2  # Increase beta for syph transmission
 sim = ss.Sim(people=ppl, networks=nwk, modules=syph)
 
 


Test patch:
diff --git a/tests/baseline.json b/tests/baseline.json
index 1986c88b..c97384d8 100644
--- a/tests/baseline.json
+++ b/tests/baseline.json
@@ -1,13 +1,14 @@
 {
   "summary": {
-    "n_alive": 20824.4,
-    "new_deaths": 0.0,
-    "pregnancy_pregnancies": 257.8,
-    "pregnancy_births": 244.22857142857143,
-    "hiv_n_susceptible": 20002.742857142857,
-    "hiv_n_infected": 821.6571428571428,
+    "n_alive": 24216.114285714284,
+    "new_deaths": 61.20952380952381,
+    "pregnancy_pregnancies": 320.6095238095238,
+    "pregnancy_births": 303.6952380952381,
+    "hiv_n_susceptible": 23795.180952380953,
+    "hiv_n_infected": 420.93333333333334,
     "hiv_n_on_art": 0.0,
-    "hiv_prevalence": 0.0410060553015017,
-    "hiv_new_infections": 12.0
+    "hiv_prevalence": 0.022111498205925655,
+    "hiv_new_infections": 18.057142857142857,
+    "hiv_new_deaths": 61.20952380952381
   }
 }
\ No newline at end of file
diff --git a/tests/benchmark.json b/tests/benchmark.json
index 27ac21c4..9687bac1 100644
--- a/tests/benchmark.json
+++ b/tests/benchmark.json
@@ -1,13 +1,13 @@
 {
   "time": {
     "people": 0.018,
-    "initialize": 0.038,
-    "run": 0.347
+    "initialize": 0.043,
+    "run": 0.355
   },
   "parameters": {
     "n_agents": 10000,
     "n_years": 20,
     "dt": 0.2
   },
-  "cpu_performance": 0.9276295186623464
+  "cpu_performance": 0.9446053779432078
 }
\ No newline at end of file
diff --git a/tests/conftest.py b/tests/conftest.py
new file mode 100644
index 00000000..3af303b5
--- /dev/null
+++ b/tests/conftest.py
@@ -0,0 +1,15 @@
+# Fixtures defined in this file will be used by all tests
+
+import stisim as ss
+import pytest
+
+# Run all tests with both single and multi RNG streams
+@pytest.fixture(params=['single','multi'], autouse=True)
+def set_rng_type(request):
+    original_rng_type = ss.options.multirng
+    if request.param == 'single':
+        ss.options.multirng = False
+    else:
+        ss.options.multirng = True
+    yield
+    ss.options.multirng = original_rng_type
\ No newline at end of file
diff --git a/tests/devtests/devtest_art_impact_viz.py b/tests/devtests/devtest_art_impact_viz.py
new file mode 100644
index 00000000..84a16b2f
--- /dev/null
+++ b/tests/devtests/devtest_art_impact_viz.py
@@ -0,0 +1,375 @@
+"""
+Compare two HIV simulations, one baseline and the other with ART
+"""
+
+# %% Imports and settings
+import stisim as ss
+import matplotlib.pyplot as plt
+import pandas as pd
+import numpy as np
+import networkx as nx
+import sys
+import os
+import argparse
+import sciris as sc
+
+default_n_agents = 25
+# Three choices for network here, note that only the first two are common-random-number safe
+network = ['stable_monogamy', 'embedding', 'hpv_network'][1]
+
+do_plot_graph = True
+# Several choices for how to layout the graph when plotting
+kind = ['radial', 'bipartite', 'spring', 'multipartite'][1]
+
+do_plot_longitudinal = True
+do_plot_timeseries = True
+
+ss.options(multirng = True) # Can set multirng to False for comparison
+
+figdir = os.path.join(os.getcwd(), 'figs', network)
+sc.path(figdir).mkdir(parents=True, exist_ok=True)
+
+
+class stable_monogamy(ss.SexualNetwork):
+    """
+    Very simple network for debugging in which edges are:
+    1-2, 3-4, 5-6, ...
+    """
+    def __init__(self, **kwargs):
+        # Call init for the base class, which sets all the keys
+        super().__init__(**kwargs)
+        return
+
+    def initialize(self, sim):
+        n = len(sim.people._uid_map)
+        n_edges = n//2
+        self.contacts.p1 = np.arange(0, 2*n_edges, 2) # EVEN
+        self.contacts.p2 = np.arange(1, 2*n_edges, 2) # ODD
+        self.contacts.beta = np.ones(n_edges)
+        return
+
+
+class Graph():
+    def __init__(self, nodes, edges):
+        self.graph = nx.from_pandas_edgelist(df=edges, source='p1', target='p2', edge_attr=True)
+        self.graph.add_nodes_from(nodes.index)
+        nx.set_node_attributes(self.graph, nodes.transpose().to_dict())
+        return
+
+    def draw_nodes(self, filter, pos, ax, **kwargs):
+        inds = [i for i,n in self.graph.nodes.data() if filter(n)]
+        nc = ['red' if nd['hiv'] else 'lightgray' for i, nd in self.graph.nodes.data() if i in inds]
+        ec = ['green' if nd['on_art'] else 'black' for i, nd in self.graph.nodes.data() if i in inds]
+        if inds:
+            nx.draw_networkx_nodes(self.graph, nodelist=inds, pos=pos, ax=ax, node_color=nc, edgecolors=ec, **kwargs)
+        return
+
+    def plot(self, pos, edge_labels=False, ax=None):
+        kwargs = dict(node_shape='x', node_size=250, linewidths=2, ax=ax, pos=pos)
+        self.draw_nodes(lambda n: n['dead'], **kwargs)
+
+        kwargs['node_shape'] = 'o'
+        self.draw_nodes(lambda n: not n['dead'] and n['female'], **kwargs)
+        
+        kwargs['node_shape'] = 's'
+        self.draw_nodes(lambda n: not n['dead'] and not n['female'], **kwargs)
+
+        nx.draw_networkx_edges(self.graph, pos=pos, ax=ax)
+        nx.draw_networkx_labels(self.graph, labels={i:int(a['cd4']) for i,a in self.graph.nodes.data()}, font_size=8, pos=pos, ax=ax)
+        if edge_labels:
+            nx.draw_networkx_edge_labels(self.graph, edge_labels={(i,j): int(a['dur']) if 'dur' in a else np.nan for i,j,a in self.graph.edges.data()}, font_size=8, pos=pos, ax=ax)
+        return
+
+
+class GraphAnalyzer(ss.Analyzer):
+    ''' Simple analyzer to assess if common random numbers are working '''
+
+    def __init__(self, **kwargs):
+        super().__init__(**kwargs) # Initialize the Analyzer object
+
+        self.graphs = {}
+        return
+
+    def initialize(self, sim):
+        self.initialized = True
+        self.update_results(sim, init=True)
+        return
+
+    def update_results(self, sim, init=False):
+        nodes = pd.DataFrame({
+            'age': sim.people.age.values,
+            'female': sim.people.female.values,
+            'dead': sim.people.dead.values,
+            'hiv': sim.people.hiv.infected.values,
+            'on_art': sim.people.hiv.on_art.values,
+            'cd4': sim.people.hiv.cd4.values,
+        })
+
+        edges = sim.people.networks[network].to_df()
+
+        idx = sim.ti if not init else -1
+        self.graphs[idx] = Graph(nodes, edges)
+        return
+
+    def finalize(self, sim):
+        super().finalize(sim)
+        return
+
+
+def run_sim(n=25, rand_seed=0, intervention=False, analyze=False, lbl=None):
+    ppl = ss.People(n)
+
+    try:
+        net_class = eval(network)
+    except Exception as e:
+        net_class = getattr(ss, network)
+    ppl.networks = ss.ndict(net_class(), ss.maternal())
+
+    hiv_pars = {
+        'beta': {network: [0.3, 0.25], 'maternal': [0.2, 0]},
+        'init_prev': 0.25,
+        'art_efficacy': 0.96,
+    }
+    hiv = ss.HIV(hiv_pars)
+
+    art = ss.hiv.ART(0, 0.6) # 60% coverage from day 0
+
+    pars = {
+        'start': 1980,
+        'end': 2020,
+        'remove_dead': False, # So we can see who dies, sim results should not change with True
+        'interventions': [art] if intervention else [],
+        'rand_seed': rand_seed,
+        'analyzers': [GraphAnalyzer()] if analyze else [],
+    }
+    sim = ss.Sim(people=ppl, diseases=[hiv], demographics=[ss.Pregnancy(), ss.background_deaths()], pars=pars, label=lbl)
+    sim.initialize()
+
+    # Infect every other person, useful for exploration in conjunction with init_prev=0
+    #sim.diseases['hiv'].set_prognoses(sim, np.arange(0,n,2), from_uids=None)
+
+    sim.run()
+
+    return sim
+
+
+def run_scenario(n=10, rand_seed=0, analyze=True):
+    sims = sc.parallelize(run_sim,
+                          kwargs={'n':n, 'analyze': analyze, 'rand_seed': rand_seed},
+                          iterkwargs=[{'intervention':False, 'lbl':'Baseline'}, {'intervention':True, 'lbl':'Intervention'}], die=True)
+
+    for i, sim in enumerate(sims):
+        sim.save(os.path.join(figdir, f'sim{i}.obj'))
+
+    return sims
+
+
+def getpos(ti, g1, g2, guess=None, kind='bipartite'):
+
+    n1 = dict(g1[ti].graph.nodes.data())
+    n2 = dict(g2[ti].graph.nodes.data())
+    nodes = sc.mergedicts(n2, n1)
+    n = len(nodes)
+
+    if kind == 'radial':
+        pos = {i:(np.cos(2*np.pi*i/n), np.sin(2*np.pi*i/n)) for i in range(n)}
+        if guess:
+            if len(guess) < n:
+                pos = {i:(np.cos(2*np.pi*i/n), np.sin(2*np.pi*i/n)) for i in range(n)}
+
+    elif kind == 'spring':
+        pos = nx.spring_layout(g1[ti].graph, k=None, pos=guess, fixed=None, iterations=50, threshold=0.0001, weight=None, scale=1, center=None, dim=2, seed=None)
+        if guess:
+            pos = sc.mergedicts(pos, guess)
+
+    elif kind == 'multipartite':
+        pos = nx.multipartite_layout(g1[ti].graph, subset_key='female', align='vertical', scale=10, center=None)
+        if guess:
+            pos = sc.mergedicts(pos, guess)
+
+        if guess:
+            for i in guess.keys():
+                pos[i] = (pos[i][0], guess[i][1]) # Keep new x but carry over y
+
+    elif kind == 'bipartite':
+        pos = {i:(nd['age'], 2*nd['female']-1 + np.random.uniform(-0.3, 0.3)) for i, nd in nodes.items()}
+
+        if guess:
+            for i in guess.keys():
+                pos[i] = (pos[i][0], guess[i][1]) # Keep new x but carry over y
+
+    return pos
+
+
+def plot_graph(sim1, sim2):
+    g1 = sim1.analyzers[0].graphs
+    g2 = sim2.analyzers[0].graphs
+
+    n = len(g1[-1].graph)
+    el = n <= 25 # Draw edge labels
+
+    fig, axv = plt.subplots(1, 2, figsize=(10,5))
+    global ti
+    timax = sim1.tivec[-1]
+
+    global pos
+    pos = {}
+    pos[-1] = getpos(0, g1, g2, kind=kind)
+    for ti in range(timax+1):
+        pos[ti] = getpos(ti, g1, g2, guess=pos[ti-1], kind=kind)
+
+    ti = -1 # Initial state is -1, representing the state before the first step
+
+    def on_press(event):
+        print('press', event.key)
+        sys.stdout.flush()
+        global ti, pos
+        if event.key == 'right':
+            ti = min(ti+1, timax)
+        elif event.key == 'left':
+            ti = max(ti-1, -1)
+
+        # Clear
+        axv[0].clear()
+        axv[1].clear()
+
+        g1[ti].plot(pos[ti], edge_labels=el, ax=axv[0])
+        g2[ti].plot(pos[ti], edge_labels=el, ax=axv[1])
+        fig.suptitle(f'Time is {ti} (use the arrow keys to change)')
+        axv[0].set_title(sim1.label)
+        axv[1].set_title(sim2.label)
+        fig.canvas.draw()
+
+    fig.canvas.mpl_connect('key_press_event', on_press)
+
+    g1[ti].plot(pos[ti], edge_labels=el, ax=axv[0])
+    g2[ti].plot(pos[ti], edge_labels=el, ax=axv[1])
+    fig.suptitle(f'Time is {ti} (use the arrow keys to change)')
+    axv[0].set_title(sim1.label)
+    axv[1].set_title(sim2.label)
+    return fig
+
+
+def plot_ts():
+    # Plot timeseries summary
+    fig, axv = plt.subplots(2,1, sharex=True)
+    axv[0].plot(sim1.tivec, sim1.results.hiv.n_infected, label=sim1.label)
+    axv[0].plot(sim2.tivec, sim2.results.hiv.n_infected, ls=':', label=sim2.label)
+    axv[0].set_title('HIV number of infections')
+
+    axv[1].plot(sim1.tivec, sim1.results.hiv.new_deaths, label=sim1.label)
+    axv[1].plot(sim2.tivec, sim2.results.hiv.new_deaths, ls=':', label=sim2.label)
+    axv[1].set_title('HIV Deaths')
+
+    plt.legend()
+    return fig
+
+
+def analyze_people(sim):
+    p = sim.people
+    ever_alive = ss.false(np.isnan(p.age))
+    years_lived = np.full(len(p), sim.ti+1) # Actually +1 dt here, I think
+    years_lived[p.dead] = p.ti_dead[p.dead]
+    years_lived = years_lived[ever_alive] # Trim, could be more efficient
+    age_initial = p.age[ever_alive].values - years_lived
+    age_initial = age_initial.astype(np.float32) # For better hash comparability, there are small differences at float64
+
+    df = pd.DataFrame({
+        'id': [hash((p.slot[i], age_initial[i], p.female[i])) for i in ever_alive], # if slicing, don't need ._view,
+        'age_initial': age_initial,
+        'years_lived': years_lived,
+        'ti_infected': p.hiv.ti_infected[ever_alive].values,
+        'ti_art': p.hiv.ti_art[ever_alive].values,
+        'ti_dead': p.ti_dead[ever_alive].values,
+
+        # Useful for debugging, but not needed for plotting
+        'slot': p.slot[ever_alive].values,
+        'female': p.female[ever_alive].values,
+    })
+    df.replace(to_replace=ss.INT_NAN, value=np.nan, inplace=True)
+    df['age_infected'] = df['age_initial'] + df['ti_infected']
+    df['age_art']      = df['age_initial'] + df['ti_art']
+    df['age_dead']     = df['age_initial'] + df['ti_dead']
+    return df
+
+
+def plot_longitudinal(sim1, sim2):
+
+    df1 = analyze_people(sim1)
+    df1['sim'] = 'Baseline'
+    df2 = analyze_people(sim2)
+    df2['sim'] = 'With ART'
+
+    df = pd.concat([df1, df2]).set_index('id')
+
+    df['ypos'] = pd.factorize(df.index.values)[0]
+    N = df['sim'].nunique()
+    height = 0.5/N
+
+    fig, ax = plt.subplots(figsize=(10,6))
+
+    # For the legend:
+    plt.barh(y=0, left=0, width=1e-6, color='k', height=height, label='Alive')
+    plt.barh(y=0, left=0, width=1e-6, color='m', height=height, label='Infected before birth')
+    plt.barh(y=0, left=0, width=1e-6, color='r', height=height, label='Infected')
+    plt.barh(y=0, left=0, width=1e-6, color='g', height=height, label='ART')
+    plt.scatter(y=0, x=0, color='c', marker='|', label='Death')
+
+    for n, (lbl, data) in enumerate(df.groupby('sim')):
+        yp = data['ypos'] + n/(N+1) # Leave space
+
+        ti_initial = np.maximum(-data['age_initial'], 0)
+        ti_final = data['ti_dead'].fillna(40)
+        plt.barh(y=yp, left=ti_initial, width=ti_final - ti_initial, color='k', height=height)
+
+        # Infected before birth
+        vertical = data['age_infected']<0
+        plt.barh(y=yp[vertical], left=data.loc[vertical]['ti_infected'], width=ti_final[vertical]-data.loc[vertical]['ti_infected'], color='m', height=height)
+
+        # Infected
+        infected = ~data['ti_infected'].isna()
+        ai = data.loc[infected]['age_infected'].values # Adjust for vertical transmission
+        ai[~(ai<0)] = 0
+        plt.barh(y=yp[infected], left=data.loc[infected]['ti_infected']-ai, width=ti_final[infected]-data.loc[infected]['ti_infected']+ai, color='r', height=height)
+
+        # ART
+        art = ~data['ti_art'].isna()
+        plt.barh(y=yp[art], left=data.loc[art]['ti_art'], width=ti_final[art]-data.loc[art]['ti_art'], color='g', height=height)
+
+        # Dead
+        dead = ~data['ti_dead'].isna()
+        plt.scatter(y=yp[dead], x=data.loc[dead]['ti_dead'], color='c', marker='|')
+
+    ax.set_xlabel('Age (years)')
+    ax.set_ylabel('UID')
+    ax.legend()
+
+    return fig
+
+
+if __name__ == '__main__':
+    parser = argparse.ArgumentParser()
+    parser.add_argument('-p', '--plot', help='Plot from a cached CSV file', type=str)
+    parser.add_argument('-n', help='Number of agents', type=int, default=default_n_agents)
+    parser.add_argument('-s', help='Rand seed', type=int, default=0)
+    args = parser.parse_args()
+
+    if args.plot:
+        print('Reading files', args.plot)
+        sim1 = sc.load(os.path.join(args.plot, 'sim1.obj'))
+        sim2 = sc.load(os.path.join(args.plot, 'sim2.obj'))
+    else:
+        print('Running scenarios')
+        [sim1, sim2] = run_scenario(n=args.n, rand_seed=args.s)
+
+    if do_plot_longitudinal:
+        plot_longitudinal(sim1, sim2)
+
+    if do_plot_graph:
+        plot_graph(sim1, sim2)
+
+    if do_plot_timeseries:
+        plot_ts()
+
+    plt.show()
+    print('Done')
\ No newline at end of file
diff --git a/tests/devtests/devtest_births.py b/tests/devtests/devtest_births.py
index 98e801b9..ac03b11e 100644
--- a/tests/devtests/devtest_births.py
+++ b/tests/devtests/devtest_births.py
@@ -17,15 +17,15 @@
 realistic_birth = {'birth_rates': pd.read_csv(ss.root/'tests/test_data/nigeria_births.csv')}
 realistic_death = {'death_rates': pd.read_csv(ss.root/'tests/test_data/nigeria_deaths.csv')}
 
-series_death = {'death_rates': pd.Series(
+asmr = {'death_rates': pd.Series(
             index=[0, 10, 20, 25, 30, 35, 40, 45, 50, 55, 60, 65, 70, 75, 80, 85, 90, 95],
             data=[0.0046355, 0.000776, 0.0014232, 0.0016693, 0.0021449, 0.0028822, 0.0039143, 0.0053676, 0.0082756, 0.01, 0.02, 0.03, 0.04, 0.06, 0.11, 0.15, 0.21, 0.30],
         )}
 
 births = ss.births(realistic_birth)
-deaths = ss.background_deaths(series_death)
-gon = ss.Gonorrhea({'p_death': 0.5, 'initial': 1000})
-sim = ss.Sim(people=ppl, demographics=[births, deaths], diseases=gon, networks=ss.simple_sexual(), n_years=100)
+deaths = ss.background_deaths(asmr)
+gon = ss.Gonorrhea({'p_death': 0.5, 'init_prev': 0.1})
+sim = ss.Sim(people=ppl, demographics=[births, deaths], diseases=gon, networks=ss.mf(), n_years=100)
 sim.initialize()
 sim.run()
 
@@ -37,5 +37,6 @@
 ax[0].set_title('Births and deaths')
 ax[1].set_title('Population size')
 ax[0].legend()
+fig.tight_layout()
 
 plt.show()
\ No newline at end of file
diff --git a/tests/devtests/devtest_connectors.py b/tests/devtests/devtest_connectors.py
index 992dba7d..f557542f 100644
--- a/tests/devtests/devtest_connectors.py
+++ b/tests/devtests/devtest_connectors.py
@@ -7,28 +7,30 @@
 import matplotlib.pyplot as plt
 
 ng = ss.Gonorrhea()
-ng.pars['beta'] = {'simple_sexual': [0.05, 0.025]}
+ng.pars['beta'] = {'mf': [0.05, 0.025]}
 ng.pars['init_prev'] = 0.025
 
 ppl1 = ss.People(10000)
-ppl1.networks = ss.ndict(ss.simple_sexual())
+ppl1.networks = ss.ndict(ss.mf())
 sim_nohiv = ss.Sim(people=ppl1, diseases=ng)
 sim_nohiv.run()
 
 hiv = ss.HIV()
-hiv.pars['beta'] = {'simple_sexual': [0.0008, 0.0004]}
+hiv.pars['beta'] = {'mf': [0.0008, 0.0004]}
 hiv.pars['init_prev'] = 0.05
 ng = ss.Gonorrhea()
-ng.pars['beta'] = {'simple_sexual': [0.05, 0.025]}
+ng.pars['beta'] = {'mf': [0.05, 0.025]}
 ng.pars['init_prev'] = 0.025
 ppl2 = ss.People(10000)
-ppl2.networks = ss.ndict(ss.simple_sexual())
+ppl2.networks = ss.ndict(ss.mf())
 sim_hiv = ss.Sim(people=ppl2, diseases=[hiv, ng], connectors=ss.simple_hiv_ng())
 sim_hiv.run()
 
 plt.figure()
-plt.plot(sim_hiv.tivec, sim_hiv.results.gonorrhea.n_infected, label='With HIV')
-plt.plot(sim_nohiv.tivec, sim_nohiv.results.gonorrhea.n_infected, label='Without HIV')
+plt.plot(sim_hiv.yearvec, sim_hiv.results.gonorrhea.n_infected, label='With HIV')
+plt.plot(sim_nohiv.yearvec, sim_nohiv.results.gonorrhea.n_infected, label='Without HIV')
 plt.title('Gonorrhea infections')
+plt.xlabel('Year')
+plt.ylabel('Count')
 plt.legend()
-plt.show()
+plt.show()
\ No newline at end of file
diff --git a/tests/devtests/devtest_dists.py b/tests/devtests/devtest_dists.py
index 489dba04..600e160a 100644
--- a/tests/devtests/devtest_dists.py
+++ b/tests/devtests/devtest_dists.py
@@ -12,5 +12,5 @@
 plt.subplot(121)
 plt.hist(data, 50)
 plt.subplot(122)
-plt.hist(ss.data_dist(vals=hist, bins=bins).sample(10_000), 50)
-plt.show()
+plt.hist(ss.data_dist(vals=hist, bins=bins)._sample(10_000), 50)
+plt.show()
\ No newline at end of file
diff --git a/tests/devtests/devtest_networks.py b/tests/devtests/devtest_networks.py
index 4c779412..536623fb 100644
--- a/tests/devtests/devtest_networks.py
+++ b/tests/devtests/devtest_networks.py
@@ -32,6 +32,8 @@
 sim.run()
 
 plt.figure()
-plt.plot(sim.tivec, sim.results.hiv.n_infected)
+plt.plot(sim.yearvec, sim.results.hiv.n_infected)
+plt.xlabel('Year')
+plt.ylabel('Count')
 plt.title('HIV number of infections')
-plt.show()
+plt.show()
\ No newline at end of file
diff --git a/tests/devtests/devtest_orientation.py b/tests/devtests/devtest_orientation.py
index a7d1194c..80fd9782 100644
--- a/tests/devtests/devtest_orientation.py
+++ b/tests/devtests/devtest_orientation.py
@@ -2,7 +2,6 @@
 Experimenting with a sexual orentation class
 """
 
-import sciris
 import stisim as ss
 import numpy as np
 
@@ -100,6 +99,3 @@ def get_preferences(self, potentials):
 # Pairs:
 # [p1, p2], [p3, p4]
 # [p1, p4], [na, na]
-
-
-
diff --git a/tests/devtests/devtest_remove_people.py b/tests/devtests/devtest_remove_people.py
index 87cf5ccb..37aa8757 100644
--- a/tests/devtests/devtest_remove_people.py
+++ b/tests/devtests/devtest_remove_people.py
@@ -33,10 +33,9 @@ def run_test(remove_dead=True, do_plot=False, rand_seed=0, verbose=False):
 
     births = ss.births(realistic_birth)
     deaths = ss.background_deaths(series_death)
-    gon = ss.Gonorrhea({'p_death': 0.5, 'initial': 1000})
+    gon = ss.Gonorrhea({'p_death': 0.5, 'init_prev': 0.1})
 
-    sim = ss.Sim(people=ppl, demographics=[births,deaths], diseases=[gon], networks=ss.simple_sexual(), remove_dead=remove_dead, n_years=100, analyzers=[agent_analyzer()], rand_seed=rand_seed, verbose=verbose)
-    # sim = ss.Sim(people=ppl, modules=[births, deaths], networks=ss.simple_sexual(), remove_dead=remove_dead, n_years=100, analyzers=[agent_analyzer()], rand_seed=rand_seed, verbose=False)
+    sim = ss.Sim(people=ppl, demographics=[births,deaths], diseases=[gon], networks=ss.mf(), remove_dead=remove_dead, n_years=100, analyzers=[agent_analyzer()], rand_seed=rand_seed, verbose=verbose)
 
     sim.initialize()
     sim.run()
@@ -97,3 +96,4 @@ def compare_totals(rand_seed=0):
     ax[3].set_title('Number of agents')
 
     fig.tight_layout()
+    plt.show()
diff --git a/tests/devtests/devtest_rng.py b/tests/devtests/devtest_rng.py
new file mode 100644
index 00000000..379c1bcb
--- /dev/null
+++ b/tests/devtests/devtest_rng.py
@@ -0,0 +1,153 @@
+"""
+Exploring singlerng vs multirng for a model HIV system sweeping coverage of ART.
+"""
+
+# %% Imports and settings
+import os
+import stisim as ss
+import sciris as sc
+import pandas as pd
+import seaborn as sns
+import numpy as np
+
+# Suppress warning from seaborn
+import warnings
+warnings.filterwarnings("ignore", "is_categorical_dtype")
+warnings.filterwarnings("ignore", "use_inf_as_na")
+
+n = 100 # Agents
+n_rand_seeds = 25
+intv_cov_levels = [0.01, 0.10, 0.25, 0.73] + [0] # Must include 0 as that's the baseline
+
+figdir = os.path.join(os.getcwd(), 'figs', 'ART')
+sc.path(figdir).mkdir(parents=True, exist_ok=True)
+
+def run_sim(n, idx, intv_cov, rand_seed, multirng):
+
+    print(f'Starting sim {idx} with rand_seed={rand_seed} and intv_cov={intv_cov}, multirng={multirng}')
+
+    ppl = ss.People(n)
+
+    ppl.networks = ss.ndict( ss.embedding(), ss.maternal())
+
+    hiv_pars = {
+        'beta': {'embedding': [0.2, 0.15], 'maternal': [0.3, 0]},
+        'init_prev': np.maximum(10/n, 0.01),
+        'art_efficacy': 0.96,
+    }
+    hiv = ss.HIV(hiv_pars)
+
+    pregnancy = ss.Pregnancy()
+    deaths = ss.background_deaths()
+
+    pars = {
+        'start': 1980,
+        'end': 2070,
+        'rand_seed': rand_seed,
+        'verbose': 0,
+        'remove_dead': True,
+        'slot_scale': 10 # Increase slot scale to reduce repeated slots
+    }
+
+    if intv_cov > 0:
+        pars['interventions'] = [ ss.hiv.ART(t=[0, 10, 20], coverage=[0, intv_cov/3, intv_cov]) ]
+
+    sim = ss.Sim(people=ppl, diseases=[hiv], demographics=[pregnancy, deaths], pars=pars, label=f'Sim with {n} agents and intv_cov={intv_cov}')
+    sim.initialize()
+    sim.run()
+
+    df = pd.DataFrame( {
+        'ti': sim.tivec,
+        #'hiv.n_infected': sim.results.hiv.n_infected, # Optional, but mostly redundant with prevalence
+        'hiv.prevalence': sim.results.hiv.prevalence,
+        'hiv.cum_deaths': sim.results.hiv.new_deaths.cumsum(),
+        'pregnancy.cum_births': sim.results.pregnancy.births.cumsum(),
+    })
+    df['intv_cov'] = intv_cov
+    df['rand_seed'] = rand_seed
+    df['multirng'] = multirng
+
+    print(f'Finishing sim {idx} with rand_seed={rand_seed} and intv_cov={intv_cov}, multirng={multirng}')
+
+    return df
+
+def run_scenarios():
+    results = []
+    times = {}
+    for multirng in [True, False]:
+        ss.options(multirng=multirng)
+        cfgs = []
+        for rs in range(n_rand_seeds):
+            for intv_cov in intv_cov_levels:
+                cfgs.append({'intv_cov':intv_cov, 'rand_seed':rs, 'multirng':multirng, 'idx':len(cfgs)})
+        T = sc.tic()
+        results += sc.parallelize(run_sim, kwargs={'n': n}, iterkwargs=cfgs, die=True, serial=False)
+        times[f'multirng={multirng}'] = sc.toc(T, output=True)
+
+    print('Timings:', times)
+
+    df = pd.concat(results)
+    df.to_csv(os.path.join(figdir, 'result.csv'))
+    return df
+
+def plot_scenarios(df):
+    d = pd.melt(df, id_vars=['ti', 'rand_seed', 'intv_cov', 'multirng'], var_name='channel', value_name='Value')
+    d['baseline'] = d['intv_cov']==0
+    bl = d.loc[d['baseline']]
+    scn = d.loc[~d['baseline']]
+    bl = bl.set_index(['ti', 'channel', 'rand_seed', 'intv_cov', 'multirng'])[['Value']].reset_index('intv_cov')
+    scn = scn.set_index(['ti', 'channel', 'rand_seed', 'intv_cov', 'multirng'])[['Value']].reset_index('intv_cov')
+    mrg = scn.merge(bl, on=['ti', 'channel', 'rand_seed', 'multirng'], suffixes=('', '_ref'))
+    mrg['Value - Reference'] = mrg['Value'] - mrg['Value_ref']
+    mrg = mrg.sort_index()
+
+    fkw = {'sharey': False, 'sharex': 'col', 'margin_titles': True}
+
+    ## TIMESERIES
+    g = sns.relplot(kind='line', data=d, x='ti', y='Value', hue='intv_cov', col='channel', row='multirng',
+        height=5, aspect=1.2, palette='Set1', errorbar='sd', lw=2, facet_kws=fkw)
+    g.set_titles(col_template='{col_name}', row_template='multirng: {row_name}')
+    g.set_xlabels(r'$t_i$')
+    g.fig.savefig(os.path.join(figdir, 'timeseries.png'), bbox_inches='tight', dpi=300)
+
+    ## DIFF TIMESERIES
+    for ms, mrg_by_ms in mrg.groupby('multirng'):
+        g = sns.relplot(kind='line', data=mrg_by_ms, x='ti', y='Value - Reference', hue='intv_cov', col='channel', row='intv_cov',
+            height=3, aspect=1.0, palette='Set1', estimator=None, units='rand_seed', lw=0.5, facet_kws=fkw) #errorbar='sd', lw=2, 
+        g.set_titles(col_template='{col_name}', row_template='Coverage: {row_name}')
+        g.fig.suptitle('MultiRNG' if ms else 'SingleRNG')
+        g.fig.subplots_adjust(top=0.88)
+        g.set_xlabels(r'Timestep $t_i$')
+        g.fig.savefig(os.path.join(figdir, 'diff_multi.png' if ms else 'diff_single.png'), bbox_inches='tight', dpi=300)
+
+    ## FINAL TIME
+    tf = df['ti'].max()
+    mtf = mrg.loc[tf]
+    g = sns.displot(data=mtf.reset_index(), kind='kde', fill=True, rug=True, cut=0, hue='intv_cov', x='Value - Reference', col='channel', row='multirng',
+        height=5, aspect=1.2, facet_kws=fkw, palette='Set1')
+    g.set_titles(col_template='{col_name}', row_template='multirng: {row_name}')
+    g.set_xlabels(f'Value - Reference at $t_i={{{tf}}}$')
+    g.fig.savefig(os.path.join(figdir, 'final.png'), bbox_inches='tight', dpi=300)
+
+    print('Figures saved to:', os.path.join(os.getcwd(), figdir))
+
+    return
+
+
+if __name__ == '__main__':
+    import argparse
+
+    parser = argparse.ArgumentParser()
+    parser.add_argument('-p', '--plot', help='plot from a cached CSV file', type=str)
+    args = parser.parse_args()
+
+    if args.plot:
+        print('Reading CSV file', args.plot)
+        df = pd.read_csv(args.plot, index_col=0)
+    else:
+        print('Running scenarios')
+        df = run_scenarios()
+
+    plot_scenarios(df)
+
+    print('Done')
\ No newline at end of file
diff --git a/tests/devtests/devtest_sensitivity.py b/tests/devtests/devtest_sensitivity.py
new file mode 100644
index 00000000..973631d5
--- /dev/null
+++ b/tests/devtests/devtest_sensitivity.py
@@ -0,0 +1,158 @@
+"""
+A simple beta sweep with multisim on and off using HIV as an example pathogen
+"""
+
+# %% Imports and settings
+import os
+import stisim as ss
+import sciris as sc
+import pandas as pd
+import seaborn as sns
+
+# Suppress warning from seaborn
+import warnings
+warnings.filterwarnings("ignore", "is_categorical_dtype")
+warnings.filterwarnings("ignore", "use_inf_as_na")
+
+n = 1_000 # Agents
+n_rand_seeds = 25
+xf_levels = [0.5, 0.8, 1.26, 2.0] + [1] # Must include 1 as that's the baseline | roughly np.logspace(np.log2(0.5), np.log2(20), 4, base=2)
+
+figdir = os.path.join(os.getcwd(), 'figs', 'Sweep')
+sc.path(figdir).mkdir(parents=True, exist_ok=True)
+
+def run_sim(n, xf, rand_seed, multirng):
+    ppl = ss.People(n)
+
+    rel_pars = {
+        'male_shift': 5,
+        'std': 3,
+        'pars': {
+            'dur': ss.lognormal(15*xf, 15*xf),  # Can vary by age, year, and individual pair
+            'part_rates': 0.9,  # Participation rates - can vary by sex and year
+            'rel_part_rates': 1.0,
+            'debut': 16,  # Age of debut can vary by sex, year, and individual
+        }} 
+    ppl.networks = ss.ndict(ss.embedding(**rel_pars), ss.maternal())
+
+    hiv_pars = {
+        #'beta': {'embedding': [xf * 0.30, xf * 0.25], 'maternal': [xf * 0.2, 0]},
+        'beta': {'embedding': [0.30, 0.25], 'maternal': [0.2, 0]},
+        'init_prev': 0.05
+    }
+    hiv = ss.HIV(hiv_pars)
+
+    pregnancy = ss.Pregnancy()
+
+    pars = {
+        'start': 1980,
+        'end': 2020,
+        'rand_seed': rand_seed,
+    }
+    sim = ss.Sim(people=ppl, diseases=[hiv], demographics=[pregnancy], pars=pars, label=f'Sim with {n} agents and xf={xf}')
+    sim.initialize()
+    sim.run()
+
+    df = pd.DataFrame( {
+        'ti': sim.tivec,
+        #'hiv.n_infected': sim.results.hiv.n_infected,
+        'hiv.prevalence': sim.results.hiv.prevalence,
+        'hiv.cum_deaths': sim.results.hiv.new_deaths.cumsum(),
+        'pregnancy.cum_births': sim.results.pregnancy.births.cumsum(),
+    })
+    df['xf'] = xf
+    df['rand_seed'] = rand_seed
+    df['multirng'] = multirng
+
+    return df
+
+def run_scenarios(figdir):
+    results = []
+    times = {}
+    for multirng in [True, False]:
+        ss.options(multirng=multirng)
+        cfgs = []
+        for rs in range(n_rand_seeds):
+            for xf in xf_levels:
+                cfgs.append({'xf':xf, 'rand_seed':rs, 'multirng':multirng})
+        T = sc.tic()
+        results += sc.parallelize(run_sim, kwargs={'n': n}, iterkwargs=cfgs, die=True)
+        times[f'multirng={multirng}'] = sc.toc(T, output=True)
+
+    print('Timings:', times)
+
+    df = pd.concat(results)
+    df.to_csv(os.path.join(figdir, 'result.csv'))
+    return df
+
+
+def plot_scenarios(df, figdir):
+    d = pd.melt(df, id_vars=['ti', 'rand_seed', 'xf', 'multirng'], var_name='channel', value_name='Value')
+    d['baseline'] = d['xf']==1
+    bl = d.loc[d['baseline']]
+    scn = d.loc[~d['baseline']]
+    bl = bl.set_index(['ti', 'channel', 'rand_seed', 'xf', 'multirng'])[['Value']].reset_index('xf')
+    scn = scn.set_index(['ti', 'channel', 'rand_seed', 'xf', 'multirng'])[['Value']].reset_index('xf')
+    mrg = scn.merge(bl, on=['ti', 'channel', 'rand_seed', 'multirng'], suffixes=('', '_ref'))
+    mrg['Value - Reference'] = mrg['Value'] - mrg['Value_ref']
+    mrg = mrg.sort_index()
+
+    fkw = {'sharey': False, 'sharex': 'col', 'margin_titles': True}
+
+    ## TIMESERIES
+    g = sns.relplot(kind='line', data=d, x='ti', y='Value', hue='xf', col='channel', row='multirng',
+        height=5, aspect=1.2, palette='Set1', errorbar='sd', lw=2, facet_kws=fkw)
+    g.set_titles(col_template='{col_name}', row_template='multirng: {row_name}')
+    g.set_xlabels(r'$t_i$')
+    g.fig.savefig(os.path.join(figdir, 'timeseries.png'), bbox_inches='tight', dpi=300)
+
+    ## DIFF TIMESERIES
+    for ms, mrg_by_ms in mrg.groupby('multirng'):
+        g = sns.relplot(kind='line', data=mrg_by_ms, x='ti', y='Value - Reference', hue='xf', col='channel', row='xf',
+            height=3, aspect=1.0, palette='Set1', estimator=None, units='rand_seed', lw=0.5, facet_kws=fkw) #errorbar='sd', lw=2, 
+        g.set_titles(col_template='{col_name}', row_template='Beta: {row_name}')
+        g.fig.suptitle('MultiRNG' if ms else 'SingleRNG')
+        g.fig.subplots_adjust(top=0.88)
+        g.set_xlabels(r'Value - Reference at $t_i$')
+        g.fig.savefig(os.path.join(figdir, 'diff_multi.png' if ms else 'diff_single.png'), bbox_inches='tight', dpi=300)
+
+    ## FINAL TIME
+    tf = df['ti'].max()
+    mtf = mrg.loc[tf]
+    g = sns.displot(data=mtf.reset_index(), kind='kde', fill=True, rug=True, cut=0, hue='xf', x='Value - Reference', col='channel', row='multirng',
+        height=5, aspect=1.2, facet_kws=fkw, palette='Set1')
+    g.set_titles(col_template='{col_name}', row_template='multirng: {row_name}')
+    g.set_xlabels(f'Value - Reference at $t_i={{{tf}}}$')
+    g.fig.savefig(os.path.join(figdir, 'final.png'), bbox_inches='tight', dpi=300)
+
+    ## FINAL TIME function of beta
+    dtf = d.set_index(['ti', 'rand_seed']).sort_index().loc[tf]
+    g = sns.relplot(kind='line', data=dtf.reset_index(), x='xf', y='Value', col='channel', row='multirng',
+        height=5, aspect=1.2, facet_kws=fkw, estimator=None, units='rand_seed', lw=0.25)
+    g.set_titles(col_template='{col_name}', row_template='multirng: {row_name}')
+    g.set_ylabels(f'Value at $t_i={{{tf}}}$')
+    g.fig.savefig(os.path.join(figdir, 'final_beta.png'), bbox_inches='tight', dpi=300)
+
+    print('Figures saved to:', os.path.join(os.getcwd(), figdir))
+
+    return
+
+
+if __name__ == '__main__':
+    import argparse
+
+    parser = argparse.ArgumentParser()
+    parser.add_argument('-p', '--plot', help='Folder containing cached results.csv', type=str)
+    args = parser.parse_args()
+
+    if args.plot:
+        print('Reading CSV file', args.plot)
+        df = pd.read_csv(os.path.join(args.plot, 'result.csv'), index_col=0)
+        figdir = args.plot
+    else:
+        print('Running scenarios')
+        df = run_scenarios(figdir)
+
+    plot_scenarios(df, figdir)
+
+    print('Done')
\ No newline at end of file
diff --git a/tests/devtests/devtest_simple.py b/tests/devtests/devtest_simple.py
new file mode 100644
index 00000000..95ca23e4
--- /dev/null
+++ b/tests/devtests/devtest_simple.py
@@ -0,0 +1,40 @@
+"""
+Run simplest tests
+"""
+
+# %% Imports and settings
+import stisim as ss
+import matplotlib.pyplot as plt
+
+
+ppl = ss.People(9000)
+ppl.networks = ss.ndict(ss.mf(), ss.maternal())
+
+hiv = ss.HIV()
+hiv.pars['beta'] = {'mf': [0.08, 0.04], 'maternal': [0.2, 0]}
+
+ng = ss.Gonorrhea()
+
+sim = ss.Sim(people=ppl, demographics=[ss.Pregnancy()], diseases=[hiv, ng])
+sim.initialize()
+sim.run()
+
+# Plotting
+fig, axv = plt.subplots(3,1, figsize=(10,8))
+for ax, module in zip(axv[:2], ['HIV', 'Gonorrhea']):
+    results = sim.results[module.lower()]
+    for ch in results.keys():
+        if ch[:2] == 'n_':
+            ax.plot(sim.yearvec, results[ch], label=ch)
+    ax.legend()
+    ax.set_title(module)
+
+results = sim.results.pregnancy
+for ch in results.keys():
+    axv[2].plot(sim.yearvec, sim.results.pregnancy[ch], label=ch)
+axv[2].set_title('Pregnancy')
+axv[2].legend()
+fig.supylabel('Count')
+fig.supxlabel('Year')
+fig.tight_layout()
+plt.show()
\ No newline at end of file
diff --git a/tests/removing_people.py b/tests/removing_people.py
index f1fa4a42..8dc4e61d 100644
--- a/tests/removing_people.py
+++ b/tests/removing_people.py
@@ -35,9 +35,9 @@ def __len__(self):
 
     def initialize(self, n):
 
-        self._uid_map.initialize(n)
+        self._uid_map.grow(n)
         self._uid_map[:] = np.arange(0, n)
-        self.uid.initialize(n)
+        self.uid.grow(n)
         self.uid[:] = np.arange(0, n)
 
         for state in self.states:
@@ -136,8 +136,8 @@ def assign(lbl, uids, vals, n=2_000):
                     s.values[uid_map[uids]] = vals
 
     lookup('SINGLE ITEM LOOKUP', single_item_uid)
-    lookup('MULTIPLE ITEM LOOKUP', multiple_items_uid)
-    lookup('MULTIPLE ITEM LOOKUP (FEW)', multiple_items_few_uid)
+    lookup('MULTIPLE ITEM LOOKUP (state lookup also includes construction of a FusedArray)', multiple_items_uid)
+    lookup('MULTIPLE ITEM LOOKUP (FEW) (state lookup also includes construction of a FusedArray)', multiple_items_few_uid)
     lookup('BOOLEAN ARRAY LOOKUP', boolean)
 
     assign('SINGLE ITEM ASSIGNMENT', single_item_uid, 1)
diff --git a/tests/simple.py b/tests/simple.py
deleted file mode 100644
index 2c5fa2de..00000000
--- a/tests/simple.py
+++ /dev/null
@@ -1,22 +0,0 @@
-"""
-Run simplest tests
-"""
-
-# %% Imports and settings
-import stisim as ss
-import matplotlib.pyplot as plt
-
-
-ppl = ss.People(10000)
-ppl.networks = ss.ndict(ss.simple_sexual(), ss.maternal())
-
-hiv = ss.HIV()
-hiv.pars['beta'] = {'simple_sexual': [0.0008, 0.0004], 'maternal': [0.2, 0]}
-
-sim = ss.Sim(people=ppl, demographics=[ss.Pregnancy()], diseases=[hiv, ss.Gonorrhea()])
-sim.initialize()
-sim.run()
-
-plt.figure()
-plt.plot(sim.tivec, sim.results.hiv.n_infected)
-plt.title('HIV number of infections')
\ No newline at end of file
diff --git a/tests/simple_benchmark.py b/tests/simple_benchmark.py
index 88be5929..479c6e9e 100644
--- a/tests/simple_benchmark.py
+++ b/tests/simple_benchmark.py
@@ -1,20 +1,20 @@
 import sciris as sc
 import pylab as pl
-import stisim as sti
+import stisim as ss
 
 do_plot = False
 
 with sc.timer():
-    ppl = sti.People(int(1e3))
-    ppl.networks = sti.ndict(sti.simple_sexual(), sti.maternal())
-    
-    hiv = sti.HIV()
+    ppl = ss.People(int(1e3))
+    ppl.networks = ss.ndict(ss.mf(), ss.maternal())
+
+    hiv = ss.HIV()
     hiv.pars['beta'] = {'simple_sexual': [0.0008, 0.0004], 'maternal': [0.2, 0]}
     
-    sim = sti.Sim(start=1950, end=2050, people=ppl, modules=[hiv, sti.Pregnancy()])
+    sim = ss.Sim(start=1950, end=2050, people=ppl, demographics=ss.Pregnancy(), diseases=hiv)
     sim.initialize()
     sim.run()
 
 if do_plot:
     pl.plot(sim.tivec, sim.results.hiv.n_infected)
-    pl.show()
\ No newline at end of file
+    pl.show()
diff --git a/tests/test_base.py b/tests/test_base.py
index 8f44b67d..8b7b329e 100644
--- a/tests/test_base.py
+++ b/tests/test_base.py
@@ -7,6 +7,7 @@
 import numpy as np
 import stisim as ss
 import matplotlib.pyplot as plt
+import scipy.stats as sps
 
 
 # %% Define the tests
@@ -41,15 +42,15 @@ def test_networks():
     nw1 = ss.Network(p1=p1, p2=p2, beta=beta, label='rand')
     nw2 = ss.Network(dict(p1=p1, p2=p2, beta=beta), label='rand')  # Alternate method
 
-    # Make people, then make a dynamic sexual layer and update it
-    ppl = ss.People(100)  # BasePeople
-    ppl.initialize()  # This seems to be necessary, although not completely clear why...
+    sim = ss.Sim()
+    sim.initialize()
+    
     nw3 = ss.hpv_network()
-    nw3.initialize(ppl)
-    nw3.update(ppl, ti=1, dt=1)  # Update by providing a timestep & current time index
+    nw3.initialize(sim)
+    sim.people.networks.update(sim.people)  # Update by providing a timestep & current time index
 
     nw4 = ss.maternal()
-    nw4.initialize(ppl)
+    nw4.initialize(sim)
     nw4.add_pairs(mother_inds=[1, 2, 3], unborn_inds=[100, 101, 102], dur=[1, 1, 1])
 
     return nw1, nw2, nw3, nw4
@@ -66,7 +67,7 @@ def test_microsim():
     # Set beta. The first entry represents transmission risk from infected p1 -> susceptible p2
     # Need to be careful to get the ordering right. The set-up here assumes that in the simple
     # sexual  network, p1 is male and p2 is female. In the maternal network, p1=mothers, p2=babies.
-    hiv.pars['beta'] = {'mf': [0.0008, 0.0004], 'maternal': [0.2, 0]}
+    hiv.pars['beta'] = {'mf': [0.15, 0.10], 'maternal': [0.2, 0]}
 
     sim = ss.Sim(people=ppl, demographics=ss.Pregnancy(), diseases=hiv)
     sim.initialize()
@@ -78,9 +79,19 @@ def test_microsim():
 
     return sim
 
+
 def test_ppl_construction():
 
-    sim_pars = {'networks': [ss.mf()], 'n_agents': 100}
+    def init_debut(self, sim, uids):
+        #loc = 16
+        loc = np.full(len(uids), 16)
+        loc[sim.people.female[uids]] = 21
+        return loc
+
+    mf_pars = {
+        'debut_dist': sps.norm(loc=init_debut, scale=2),  # Age of debut can vary by using callable parameter values
+    }
+    sim_pars = {'networks': [ss.mf(mf_pars)], 'n_agents': 100}
     gon_pars = {'beta': {'mf': [0.08, 0.04]}, 'p_death': 0.2}
     gon = ss.Gonorrhea(pars=gon_pars)
 
@@ -100,6 +111,8 @@ def test_ppl_construction():
     # Start timing
     T = sc.tic()
 
+    ss.options.multirng = True
+
     # Run tests
     ppl = test_people()
     nw1, nw2, nw3, nw4 = test_networks()
@@ -107,4 +120,5 @@ def test_ppl_construction():
     sim2 = test_ppl_construction()
 
     sc.toc(T)
-    print('Done.')
+    plt.show()
+    print('Done.')
\ No newline at end of file
diff --git a/tests/test_baselines.py b/tests/test_baselines.py
index e86c2084..6bdf74ad 100644
--- a/tests/test_baselines.py
+++ b/tests/test_baselines.py
@@ -78,6 +78,10 @@ def save_baseline():
 
 def test_baseline():
     ''' Compare the current default sim against the saved baseline '''
+    
+    # Do not run with multi-RNG
+    if ss.options.multirng:
+        return
 
     # Load existing baseline
     baseline = sc.loadjson(baseline_filename)
@@ -95,6 +99,10 @@ 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:
diff --git a/tests/test_distributions.py b/tests/test_distributions.py
new file mode 100644
index 00000000..7c1fe33e
--- /dev/null
+++ b/tests/test_distributions.py
@@ -0,0 +1,192 @@
+"""
+Test Distributions from distributions.py
+"""
+
+# %% Imports and settings
+import numpy as np
+import sciris as sc
+import stisim as ss
+from stisim.random import RNG
+import scipy.stats as sps
+from stisim.distributions import ScipyDistribution
+import pytest
+import matplotlib.pyplot as plt
+
+@pytest.fixture(params=[5, 50])
+def n(request):
+    yield request.param
+
+
+# %% Define the tests
+def test_basic():
+    dist = sps.norm(loc=1, scale=1) # Make a distribution
+    d = ScipyDistribution(dist)
+
+    sample = d.rvs(1)  # Draw a sample
+    samples = d.rvs(10) # Draw several samples
+
+    if ss.options.multirng:
+        print('Attempting to sample for specific UIDs, but have not provided a MultiRNG so an Exception should be raised.')
+        with pytest.raises(Exception):
+            d.rvs(np.array([1,3,8])) # Draw samples for specific uids
+    else:
+        d.rvs(np.array([1,3,8])) # Draw three samples
+
+    mu = 5 # mean (mu) of the underlying normal
+    sigma = 1 # stdev (sigma) of the underlying normal
+    s = sigma 
+    loc = 0
+    scale = np.exp(mu)
+    fig, ax = plt.subplots(1, 1)
+    x = np.linspace(sps.lognorm.ppf(0.01, s=s, loc=loc, scale=scale), sps.lognorm.ppf(0.99, s=s, loc=loc, scale=scale), 100)
+    ax.plot(x, sps.lognorm.pdf(x, s=s, loc=loc, scale=scale), 'r-', lw=5, alpha=0.6, label='lognorm pdf')
+
+    mean, var, skew, kurt = sps.lognorm.stats(s, loc=loc, scale=scale, moments='mvsk')
+    print('mean', mean, 'var', var, 'skew', skew, 'kurt', kurt)
+    print('calc mean', np.exp(mu + sigma**2/2), 'calc var', (np.exp(sigma**2)-1)*np.exp(2*mu+sigma**2)) # Check against math
+
+    temp_poisson_samples = sps.poisson(mu=2).rvs(10)  # Sample from a temporary distribution
+    return sample, samples, temp_poisson_samples
+
+def test_uniform_scalar(n):
+    """ Create a uniform distribution """
+    sc.heading('test_uniform: 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])
+    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: Testing uniform with scalar parameters')
+
+    dist = sps.uniform(loc=1, scale=4)
+    d = ScipyDistribution(dist, 'Uniform') # String here!
+    if ss.options.multirng:
+        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])
+    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: Testing uniform with callable parameters')
+
+    sim = ss.Sim().initialize()
+
+    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)
+
+    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
+
+    d = ScipyDistribution(dist)
+    draws = d.rvs(uids)
+    print(f'Uniform sample for uids {uids} returned {draws}')
+
+    assert len(draws) == len(uids)
+    return draws
+
+
+def test_repeat_slot():
+    """ Test behavior of repeated slots """
+    sc.heading('test_repeat_slot: Test behavior of repeated slots')
+
+    rng = ss.MultiRNG('Uniform')
+    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
+
+    d = ScipyDistribution(dist)
+    draws = d.rvs(uids)
+    print(f'Uniform sample for uids {uids} returned {draws}')
+
+    assert len(draws) == len(uids)
+
+    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])
+    return draws
+
+
+
+# %% Run as a script
+if __name__ == '__main__':
+    # Start timing
+    ###T = sc.tic()
+
+    n = 5
+    nTrials = 3
+
+    for multirng in [True, False]:
+        ss.options(multirng=multirng)
+        sc.heading('Testing with multirng set to', multirng)
+
+        times = []
+        for trial in range(nTrials):
+            T = sc.tic()
+
+            # Run tests - some will only pass if multirng is True
+            test_basic()
+            test_uniform_scalar(n)
+            test_uniform_scalar_str(n)
+            test_uniform_callable(n)
+            test_uniform_array(n)
+
+            times.append(sc.toc(T, doprint=False, output=True))
+
+    print(times)
+    ###sc.toc(T)
+
+    plt.show()
+    print('Done.')
\ No newline at end of file
diff --git a/tests/test_examples.py b/tests/test_examples.py
index 188b1524..e660d8ac 100644
--- a/tests/test_examples.py
+++ b/tests/test_examples.py
@@ -5,30 +5,47 @@
 # %% Imports and settings
 import stisim as ss
 import matplotlib.pyplot as plt
+import scipy.stats as sps
 import numpy as np
 
 def test_sir():
     ppl = ss.People(10000)
-    ppl.networks = ss.ndict(ss.RandomNetwork(n_contacts=ss.poisson(5)))
-    sir = ss.SIR()
+    ppl.networks = ss.ndict(ss.RandomNetwork(n_contacts=sps.poisson(mu=4)))
+
+    sir_pars = {
+        'dur_inf': sps.norm(loc=10), # Override the default distribution
+    }
+    sir = ss.SIR(sir_pars)
+
+    # You can also change the parameters of the default lognormal distribution directly!
+    #sir.pars['dur_inf'].kwds['loc'] = 5 
+    
+    # Or why not put a lambda here for fun!
+    sir.pars['dur_inf'].kwds['loc'] = lambda self, sim, uids: sim.people.age[uids]/10
+
+    sir.pars['beta'] = {'randomnetwork': 0.1}
     sim = ss.Sim(people=ppl, diseases=sir)
     sim.run()
 
-    assert len(sir.log.out_edges(np.nan)) == sir.pars.initial # Log should match initial infections
+    # CK: parameters changed
+    # assert len(sir.log.out_edges(np.nan)) == sir.pars.initial # Log should match initial infections
     df = sir.log.line_list # Check generation of line-list
-    assert df.source.isna().sum() == sir.pars.initial # Check seed infections in line list
+    # assert df.source.isna().sum() == sir.pars.initial # Check seed infections in line list
 
     plt.figure()
     plt.stackplot(
-        sim.tivec,
+        sim.yearvec,
         sir.results.n_susceptible,
         sir.results.n_infected,
         sir.results.n_recovered,
         sim.results.new_deaths.cumsum(),
     )
     plt.legend(['Susceptible', 'Infected', 'Recovered', 'Dead'])
+    plt.xlabel('Year')
+    plt.title('SIR')
+    return
 
-
+#@pytest.mark.skip(reason="Haven't converted yet")
 def test_ncd():
     ppl = ss.People(10000)
     ppl.networks = None
@@ -42,17 +59,20 @@ def test_ncd():
 
     plt.figure()
     plt.stackplot(
-        sim.tivec,
+        sim.yearvec,
         ncd.results.n_not_at_risk,
         ncd.results.n_at_risk-ncd.results.n_affected,
         ncd.results.n_affected,
         sim.results.new_deaths.cumsum(),
     )
     plt.legend(['Not at risk','At risk','Affected', 'Dead'])
+    plt.xlabel('Year')
+    plt.title('NCD')
+    return
 
 
 if __name__ == '__main__':
+    ss.options(multirng=False)
     sim1 = test_sir()
     sim2 = test_ncd()
-
-
+    plt.show()
diff --git a/tests/test_random.py b/tests/test_random.py
new file mode 100644
index 00000000..c0021d4e
--- /dev/null
+++ b/tests/test_random.py
@@ -0,0 +1,121 @@
+"""
+Test the RNG object from random.py
+"""
+
+# %% Imports and settings
+import numpy as np
+import sciris as sc
+import stisim as ss
+from stisim.random import RNG
+import pytest
+
+@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)
+    rng.initialize(rng_container, slots=slots)
+    return rng
+
+
+# %% Define the tests
+def test_random(rng, n=5):
+    """ Simple random draw """
+    sc.heading('test_random: Testing simple random draw from a RNG object')
+    draws = rng.random(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.multirng:
+        pytest.skip('Skipping multirng mode')
+    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):
+    """ Sample, reset, sample """
+    sc.heading('test_reset: Testing sample, reset, sample')
+
+    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 """
+    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}')
+
+    assert not np.array_equal(draws1, draws2)
+    return draws1, draws2
+
+
+def test_seed(n=5):
+    """ Changing seeds """
+    sc.heading('test_seed: Testing sample with seeds 0 and 1')
+
+    # NOTE: SingleRNG is using the numpy random singleton, and thus ignores base_seed here
+
+    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}')
+
+    assert not np.array_equal(draws0, draws1)
+    return draws0, draws1
+
+
+# %% Run as a script
+if __name__ == '__main__':
+    # Start timing
+    T = sc.tic()
+
+    n=5
+
+    for multirng in [True, False]:
+        ss.options(multirng=multirng)
+        sc.heading('Testing with multirng set to', multirng)
+
+        # Run tests - some will only pass if multirng is True
+        test_random(n)
+        test_reset(n)
+        test_step(n)
+        test_seed(n)
+
+    sc.toc(T)
+    print('Done.')
diff --git a/tests/test_rngcontainer.py b/tests/test_rngcontainer.py
new file mode 100644
index 00000000..8bf21707
--- /dev/null
+++ b/tests/test_rngcontainer.py
@@ -0,0 +1,194 @@
+"""
+Test the RNGContainer object from random.py
+"""
+
+# %% Imports and settings
+import numpy as np
+import sciris as sc
+import stisim as ss
+from stisim.random import NotInitializedException, SeedRepeatException, RepeatNameException
+from stisim.random import SingleRNG, MultiRNG
+import scipy.stats as sps
+import pytest
+
+
+# %% Define the tests
+
+@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)
+    return dists
+
+
+def test_rng(rng_container, rngs, n=5):
+    """ Simple sample from rng """
+    sc.heading('test_container: Testing RNGContainer object')
+
+    rng0 = rngs[0]
+    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_container: 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
+
+
+def test_seed(rng_container, rngs, n=5):
+    """ 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)
+
+    print(f'Random generators rng0 and rng1 were assigned seeds {rng0.seed} and {rng1.seed}, respectively')
+
+    assert rng1.seed != rng0.seed
+    return rng0, rng1
+
+
+def test_reset(rng_container, dists, n=5):
+    """ Sample, reset, sample """
+    sc.heading('test_reset: Testing sample, reset, sample')
+
+    d0 = dists[0]
+    if ss.options.multirng:
+        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.multirng:
+        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 """
+    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):
+    """ 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.multirng:
+        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()
+
+    _ = d1.rvs(uids)
+    s_after = d0.rvs(uids)
+
+    print(f'When sampling rng0 before rng1:', s_before)
+    print('Reset')
+    print(f'When sampling rng0 after rng1:', s_after)
+
+    if ss.options.multirng:
+        assert np.array_equal(s_before, s_after)
+    else:
+        assert not np.array_equal(s_before, s_after)
+
+    return s_before, s_after
+
+
+def test_repeatname(rng_container, n=5):
+    """ Test two random number generators with the same name """
+    sc.heading('test_repeatname: Testing if two random number generators with the same name are allowed')
+
+    rng0 = ss.RNG('test')
+    rng0.initialize(rng_container, slots=n)
+
+    rng1 = ss.RNG('test')
+    with pytest.raises(RepeatNameException):
+        rng1.initialize(rng_container, slots=n)
+    return rng0, rng1