Repository: starsimhub/starsim (language: python)

Issue:
Rethink levels of inheritance
Modularity and reuse is good, right? Except it can get really confusing to figure out what is and isn't defined. Even for SIR, we have `Module → Disease → Infection → SIR`, with each one doing meaningful work. Pars and states and methods get seamlessly merged, which is great, but makes it quite opaque what's going on under the hood (an example was @jamiecohen's attempt to define a `rel_sus` property for HPV, which conflicted with a built-in `Infection` attribute).

Proposal: be more explicit about what is being called and initialized where. Specifically:
- [x] Each class should be self-contained in defining its own states, not relying on a deep stack of `super().__init__()` calls.
- [x] Method calls should be explicit, e.g. each class should define each critical method, even if it just calls a parent method (e.g. `def make_new_cases(self, sim): super().make_new_cases(sim)`)
- [x] Rethink `sim.step()` logic and if anything can be combined/simplified; likewise for disease methods
- [x] Think about if the sim can just be stored in each module, so it doesn't have to get passed around all the time


Discussion / hints:
Another option to consider could be to incorporate more of a Mixin architecture/co-operative multiple inheritance, in that paradigm it is less about a deep tree of inheritance, and more about bringing together features.

Personally I'd quite like to see (or at least try out) class variables making a re-appearance as a means of doing things like declaring which states and parameters are defined by a particular class, that way you can easily see them upfront rather than having to dig through `__init__()` methods
Keen for iteration here as I agree it can be particularly confusing for new (and experienced!) developers alike. I like the mixin approach @RomeshA suggests as having each class own and define its own states will be quite a challenging shift from the current approach, if I understand correctly. A good place to start might be with Networks, where we have SexualNetwork and DynamicNetwork that could be mixed in?

Fix patch:
diff --git a/CHANGELOG.rst b/CHANGELOG.rst
index a527495a..63554a8d 100644
--- a/CHANGELOG.rst
+++ b/CHANGELOG.rst
@@ -7,6 +7,64 @@ What's new
 All notable changes to the codebase are documented in this file. Changes that may result in differences in model output, or are required in order to run an old parameter set with the current version, are flagged with the term "Regression information".
 
 
+Version 2.0.0 (2024-09-30)
+--------------------------
+
+Summary
+~~~~~~~
+This version contains several major breaking changes. These include: module-specific timesteps and time-aware parameters (including a day/year ``unit`` flag for modules, and  ``ss.dur()`` and ``ss.rate()`` classes for parameters), new ``ss.State`` and ``ss.Event`` classes that simplify disease logic; as well as changes to module types and integration (renaming ``update()`` and ``apply()`` methods to ``step()``).
+
+Time-aware parameters
+~~~~~~~~~~~~~~~~~~~~~
+Coming soon!
+
+States and events
+~~~~~~~~~~~~~~~~~
+Coming soon!
+
+Module changes
+~~~~~~~~~~~~~~
+- Functionality has been moved from ``ss.Plugin`` to ``ss.Module``, and the former has been removed.
+- ``ss.Connector`` functionality has been moved to ``ss.Module``. ``ss.Module`` objects can be placed anywhere in the list of modules (e.g., in demographics, networks, diseases, interventions), depending on when you want them to execute. However, ``ss.Connector`` objects are applied after ``Disease.step_state()`` and before
+- Many of the module methods have been renamed; in particular, all modules now have a ``step()`` method, which replaces ``update()`` (for demographics and networks), ``apply()`` (for interventions and analyzers), and ``make_new_cases()`` (for diseases).
+- All modules are treated the same in the integration loop, except for diseases, which have ``step_state()`` and ``step_die()`` methods.
+
+Other changes
+~~~~~~~~~~~~~
+- A new ``ss.multi_random()`` class has been added, that allows random numbers to be generated by two (or more) agents. It largely replaces ``ss.combine_rands()``.
+- Fixed a bug in how results were defined for ``ss.Births`` and ``ss.Deaths``.
+- Fixed a bug in how the ``InfectionLog`` is added to disease modules.
+- ``Infection.validate_beta()`` is now applied on every timestep, so changes to beta during the simulation are now honored.
+- ``sim.get_intervention()`` and ``sim.get_analyzer()`` have been removed; use built-in ``ndict`` operations (e.g., the label) to find the object you're after.
+- ``requires`` has been removed from modules, but ``ss.check_requires()`` is still available if needed.
+- For networks, ``contacts`` has been renamed ``edges`` except in cases where it refers to an *agent's* contacts. For example, ``network.contacts`` has been renamed ``network.edges``, but ``ss.find_contacts()`` remains the same.
+- ``Sim._orig_pars`` has been removed.
+- ``ss.unique()`` has been removed.
+
+Regression information
+~~~~~~~~~~~~~~~~~~~~~~
+- Results from Starsim v2.0 will be stochastically (but not statistically) different from Starsim v1.0.
+- ``ss.Plugin`` has been removed. Use ``ss.Module`` instead.
+- ``init_results()`` is now called by ``init_pre()``, and does not need to be called explicitly.
+- ``default_pars()`` has been renamed ``define_pars()``.
+- ``add_states()`` has been renamed ``define_states()``
+- ``Demographics.update()`` has been renamed ``Demographics.step()``.
+- ``Network.update()`` has been renamed ``Network.step()``.
+- ``Disease.update_pre()`` has been renamed ``Disease.step_pre()``.
+- ``Disease.make_new_cases()`` has been renamed ``Disease.step()``.
+- ``Disease.update_death()`` has been renamed ``Disease.die()``.
+- ``Infection._set_cases()`` has been renamed ``Infection.set_outcomes()``.
+- ``Intervention.apply(sim)`` has been renamed ``Intervention.step()``; ditto for ``Analyzer``.
+- ``Module.step()`` no longer takes ``sim`` as an argument (e.g., replace ``intervention.apply(sim)`` with ``intervention.step()``).
+- All modules now have ``init_results()`` and ``update_results()`` methods.
+- ``Network.contacts`` has been renamed ``Network.edges``.
+- ``sim.get_intervention()`` and ``sim.get_analyzer()`` have been removed; simply call directly instead (e.g. replace ``sim.get_intervention('vaccine')`` with ``sim.interventions['vaccine']``).
+- ``requires`` is no longer an attribute of modules; call the ``ss.check_requires()`` function directly if needed.
+- ``People.resolve_deaths()`` has been renamed ``People.check_deaths()``
+- ``ss.unique()`` has been removed.
+- *GitHub info*: PR TBC
+
+
 Version 1.0.1 (2024-07-22)
 ---------------------------
 - Adds a new distribution, ``ss.rand_raw()``, that samples raw integers from the random number bit generator, for use with calculating transmission. This version is roughly 20-30% faster than the previous implementation.
diff --git a/docs/run_tutorials b/docs/run_tutorials
new file mode 100755
index 00000000..42fd203d
--- /dev/null
+++ b/docs/run_tutorials
@@ -0,0 +1,62 @@
+#!/usr/bin/env python
+"""
+Run all the Jupyter notebook tutorials. Usage:
+  ./run_tutorials
+
+It will execute all notebooks in parallel, print out any errors
+encountered, and print a summary of results.
+
+NB: The execute() function was drafted in part by via ChatGPT.
+"""
+
+import os
+import sciris as sc
+import nbformat
+import nbconvert.preprocessors as nbp
+
+folder = './tutorials'
+timeout = 600
+yay = '✓'
+boo = '😢'
+
+
+def execute(path):
+    """ Executes a single Jupyter notebook and returns success/failure """
+    try:
+        with open(path) as f:
+            print(f'Executing {path}...')
+            nb = nbformat.read(f, as_version=4)
+            ep = nbp.ExecutePreprocessor(timeout=timeout, kernel_name='python3')
+            ep.preprocess(nb, {'metadata': {'path': os.path.dirname(path)}})
+        return f'{yay} {path} executed successfully.'
+    except nbp.CellExecutionError as e:
+        return f'{boo} Execution failed for {path}: {str(e)}'
+    except Exception as e:
+        return f'{boo} Error processing {path}: {str(e)}'
+
+
+def main():
+    """ Executes the notebooks in parallel and prints the results """
+    notebooks = sc.getfilepaths(folder=folder, pattern='*.ipynb')
+    results = sc.parallelize(execute, notebooks)
+    string = sc.newlinejoin(results)
+    
+    sc.heading('Results')
+    print(string)
+    
+    sc.heading('Summary')
+    n_yay = string.count(yay)
+    n_boo = string.count(boo)
+    summary = f'{n_yay} succeeded, {n_boo} failed'
+    if n_boo:
+        for i in range(len(notebooks)):
+            if boo in results[i]:
+                summary += f'\nFailed: {notebooks[i]}'
+    print(summary)
+    
+    return results
+
+
+if __name__ == '__main__':
+    results = main()
+
diff --git a/docs/tutorials/tut_diseases.ipynb b/docs/tutorials/tut_diseases.ipynb
index 2ba4a8e5..770bbb6e 100644
--- a/docs/tutorials/tut_diseases.ipynb
+++ b/docs/tutorials/tut_diseases.ipynb
@@ -60,7 +60,7 @@
    },
    "source": [
     "## Adding multiple diseases\n",
-    "You can add multiple diseases to the same simulation, like so. Here we are making use of a _connector_. A connector is a special thing in Starsim that tells you how two things relate to one another - in this case, how HIV modifies a person's transmissibility and susceptibility to gonorrhea and vice versa "
+    "You can add multiple diseases to the same simulation, like so. Here we are making use of a \"connector\". A connector is a module in Starsim that tells you how two things relate to one another - in this case, how HIV modifies a person's transmissibility and susceptibility to gonorrhea and vice versa. Unlike dieases, networks, interventions, etc., connectors don't have any pre-specified location in the sim. Instead, they can be placed wherever they make the most sense (for example, a connector that mediated how two networks behaved might be placed at the beginning or end of the list of networks; for diseases, it might be placed at the beginning or end of the list of diseases)."
    ]
   },
   {
@@ -79,11 +79,11 @@
     "import pylab as pl\n",
     "sc.options(jupyter=True)\n",
     "\n",
-    "class simple_hiv_ng(ss.Connector):\n",
+    "class simple_hiv_ng(ss.Module):\n",
     "    \"\"\" Simple connector whereby rel_sus to NG doubles if CD4 count is <200\"\"\"\n",
-    "    def __init__(self, pars=None, **kwargs):\n",
-    "        super().__init__(label='HIV-Gonorrhea', requires=[ss.HIV, ss.Gonorrhea])\n",
-    "        self.default_pars(\n",
+    "    def __init__(self, pars=None, label='HIV-Gonorrhea', **kwargs):\n",
+    "        super().__init__()\n",
+    "        self.define_pars(\n",
     "            rel_trans_hiv  = 2,\n",
     "            rel_trans_aids = 5,\n",
     "            rel_sus_hiv    = 2,\n",
@@ -92,7 +92,7 @@
     "        self.update_pars(pars, **kwargs)\n",
     "        return\n",
     "\n",
-    "    def update(self):\n",
+    "    def step(self):\n",
     "        \"\"\" Specify how HIV increases NG rel_sus and rel_trans \"\"\"\n",
     "        people = self.sim.people\n",
     "        people.gonorrhea.rel_sus[people.hiv.cd4 < 500] = self.pars.rel_sus_hiv\n",
@@ -115,7 +115,7 @@
     "\n",
     "# Make the sim, including a connector betweeh HIV and gonorrhea:\n",
     "n_agents = 5_000\n",
-    "sim = ss.Sim(n_agents=n_agents, networks='mf', diseases=[hiv, ng], connectors=simple_hiv_ng())\n",
+    "sim = ss.Sim(n_agents=n_agents, networks='mf', diseases=[simple_hiv_ng(), hiv, ng])\n",
     "sim.run()\n",
     "sim.plot('hiv')\n",
     "sim.plot('gonorrhea');"
diff --git a/starsim/demographics.py b/starsim/demographics.py
index 2c71de9c..80e9a7da 100644
--- a/starsim/demographics.py
+++ b/starsim/demographics.py
@@ -21,24 +21,17 @@ class Demographics(ss.Module):
     """
     def init_pre(self, sim):
         super().init_pre(sim)
-        self.init_results()
         return
 
     def init_results(self):
         pass
 
-    def update(self):
-        pass
-
-    def update_results(self):
-        pass
-
 
 class Births(Demographics):
     """ Create births based on rates, rather than based on pregnancy """
     def __init__(self, pars=None, metadata=None, **kwargs):
         super().__init__()
-        self.default_pars(
+        self.define_pars(
             birth_rate = 30,
             rel_birth = 1,
             units = 1e-3,  # assumes birth rates are per 1000. If using percentages, switch this to 1
@@ -83,11 +76,6 @@ def init_results(self):
         ]
         return
 
-    def update(self):
-        new_uids = self.add_births()
-        self.n_births = len(new_uids)
-        return new_uids
-
     def get_births(self):
         """
         Extract the right birth rates to use and translate it into a number of people to add.
@@ -108,6 +96,11 @@ def get_births(self):
         n_new = int(np.floor(sim.people.alive.count() * scaled_birth_prob))
         return n_new
 
+    def step(self):
+        new_uids = self.add_births()
+        self.n_births = len(new_uids)
+        return new_uids
+
     def add_births(self):
         """ Add n_new births to each state in the sim """
         people = self.sim.people
@@ -123,8 +116,8 @@ def update_results(self):
     def finalize(self):
         super().finalize()
         res = self.sim.results
-        self.results.cumulative = np.cumsum(self.results.new)
-        self.results.cbr = 1/self.pars.units*np.divide(self.results.new/self.sim.dt, res.n_alive, where=res.n_alive>0)
+        self.results.cumulative[:] = np.cumsum(self.results.new)
+        self.results.cbr[:] = 1/self.pars.units*np.divide(self.results.new/self.sim.dt, res.n_alive, where=res.n_alive>0)
         return
 
 
@@ -158,7 +151,7 @@ def __init__(self, pars=None, metadata=None, **kwargs):
                 corresponding column name in data. Similar for "sex_keys". Finally,
         """
         super().__init__()
-        self.default_pars(
+        self.define_pars(
             rel_death = 1,
             death_rate = 20,  # Default = a fixed rate of 2%/year, overwritten if data provided
             units = 1e-3,  # assumes death rates are per 1000. If using percentages, switch this to 1
@@ -237,16 +230,13 @@ def init_results(self):
         ]
         return
 
-    def update(self):
-        self.n_deaths = self.apply_deaths()
-        return
-
-    def apply_deaths(self):
+    def step(self):
         """ Select people to die """
         death_uids = self.pars.death_rate.filter()
         self.sim.people.request_death(death_uids)
-        return len(death_uids)
-
+        self.n_deaths = len(death_uids)
+        return self.n_deaths
+    
     def update_results(self):
         self.results['new'][self.sim.ti] = self.n_deaths
         return
@@ -254,8 +244,8 @@ def update_results(self):
     def finalize(self):
         super().finalize()
         n_alive = self.sim.results.n_alive
-        self.results.cumulative = np.cumsum(self.results.new)
-        self.results.cmr = 1/self.pars.units*np.divide(self.results.new / self.sim.dt, n_alive, where=n_alive>0)
+        self.results.cumulative[:] = np.cumsum(self.results.new)
+        self.results.cmr[:] = 1/self.pars.units*np.divide(self.results.new / self.sim.dt, n_alive, where=n_alive>0)
         return
 
 
@@ -263,7 +253,7 @@ class Pregnancy(Demographics):
     """ Create births via pregnancies """
     def __init__(self, pars=None, metadata=None, **kwargs):
         super().__init__()
-        self.default_pars(
+        self.define_pars(
             dur_pregnancy = 0.75, # Duration for pre-natal transmission
             dur_postpartum = ss.lognorm_ex(0.5, 0.5), # Duration for post-natal transmission (e.g. via breastfeeding)
             fertility_rate = 0, # Can be a number of Pandas DataFrame
@@ -279,7 +269,7 @@ def __init__(self, pars=None, metadata=None, **kwargs):
         self.pars.p_fertility = ss.bernoulli(p=0) # Placeholder, see make_fertility_prob_fn
         
         # Other, e.g. postpartum, on contraception...
-        self.add_states(
+        self.define_states(
             ss.BoolArr('infertile', label='Infertile'),  # Applies to girls and women outside the fertility window
             ss.BoolArr('fecund', default=True, label='Female of childbearing age'),
             ss.BoolArr('pregnant', label='Pregnant'),  # Currently pregnant
@@ -389,7 +379,7 @@ def init_results(self):
         ]
         return
 
-    def update(self):
+    def step(self):
         """ Perform all updates """
         self.update_states()
         conceive_uids = self.make_pregnancies()
diff --git a/starsim/disease.py b/starsim/disease.py
index 72a89939..a827d02a 100644
--- a/starsim/disease.py
+++ b/starsim/disease.py
@@ -19,12 +19,11 @@ class Disease(ss.Module):
 
     def __init__(self, pars=None, *args, **kwargs):
         super().__init__(*args, **kwargs)
-        self.default_pars(
+        self.define_pars(
             log = False,
         )
         self.update_pars(pars, **kwargs)
         self.results = ss.Results(self.name)
-        self.log = InfectionLog()
         return
 
     @property
@@ -43,7 +42,8 @@ def _boolean_states(self):
     def init_pre(self, sim):
         """ Link the disease to the sim, create objects, and initialize results; see Module.init_pre() for details """
         super().init_pre(sim)
-        self.init_results()
+        if self.pars.log:
+            self.log = InfectionLog()
         return
 
     def init_results(self):
@@ -58,13 +58,14 @@ def init_results(self):
             self.results += ss.Result(self.name, f'n_{state.name}', self.sim.npts, dtype=int, scale=True, label=state.label)
         return
 
-    def update_pre(self):
+    def step_state(self):
         """
-        Carry out autonomous updates at the start of the timestep (prior to transmission)
+        Carry out updates at the start of the timestep (prior to transmission);
+        these are typically state changes
         """
         pass
 
-    def update_death(self, uids):
+    def step_die(self, uids):
         """
         Carry out state changes upon death
 
@@ -81,9 +82,9 @@ def update_death(self, uids):
         """
         pass
 
-    def make_new_cases(self):
+    def step(self):
         """
-        Add new cases of the disease
+        Handle the main disease updates, e.g. add new cases
 
         This method is agnostic as to the mechanism by which new cases occur. This
         could be through transmission (parametrized in different ways, which may or
@@ -95,7 +96,7 @@ def make_new_cases(self):
         """
         pass
 
-    def set_prognoses(self, uids, source_uids=None):
+    def set_prognoses(self, uids, sources=None):
         """
         Set prognoses upon infection/acquisition
 
@@ -105,7 +106,7 @@ def set_prognoses(self, uids, source_uids=None):
         the log as part of this operation if logging is enabled (in the
         `Disease` parameters)
 
-        The `from_uids` are relevant for infectious diseases, but would be left
+        The `sources` are relevant for infectious diseases, but would be left
         as `None` for NCDs.
 
         Args:
@@ -113,14 +114,19 @@ def set_prognoses(self, uids, source_uids=None):
             uids (array): UIDs for agents to assign disease progoses to
             from_uids (array): Optionally specify the infecting agent
         """
+        # Track infections
         if self.pars.log:
-            sim = self.sim
-            if source_uids is None:
-                for target in uids:
-                    self.log.append(np.nan, target, sim.year)
-            else:
-                for target, source in zip(uids, source_uids):
-                    self.log.append(source, target, sim.year)
+            self.log_infections(uids, sources)
+        return
+
+    def log_infections(self, uids, sources=None):
+        sim = self.sim
+        if sources is None:
+            for target in uids:
+                self.log.append(np.nan, target, sim.year)
+        else:
+            for target, source in zip(uids, sources):
+                self.log.append(source, target, sim.year)
         return
 
     def update_results(self):
@@ -148,7 +154,7 @@ class Infection(Disease):
 
     def __init__(self, *args, **kwargs):
         super().__init__(*args, **kwargs)
-        self.add_states(
+        self.define_states(
             ss.BoolArr('susceptible', default=True, label='Susceptible'),
             ss.BoolArr('infected', label='Infectious'),
             ss.FloatArr('rel_sus', default=1.0, label='Relative susceptibility'),
@@ -156,38 +162,14 @@ def __init__(self, *args, **kwargs):
             ss.FloatArr('ti_infected', label='Time of infection' ),
         )
 
-        # Define random number generators for make_new_cases
-        self.rng_target = ss.rand_raw(name='target')
-        self.rng_source = ss.rand_raw(name='source')
+        # Define random number generator for determining transmission
+        self.trans_rng = ss.multi_random('source', 'target')
         return
     
     def init_pre(self, sim):
         super().init_pre(sim)
         self.validate_beta()
         return
-    
-    def validate_beta(self):
-        """
-        Perform any parameter validation
-        """
-        networks = self.sim.networks
-        if networks is not None and len(networks) > 0:
-            
-            if 'beta' not in self.pars:
-                errormsg = f'Disease {self.name} is missing beta; pars are: {sc.strjoin(self.pars.keys())}'
-                raise sc.KeyNotFoundError(errormsg)
-
-            # If beta is a scalar, apply this bi-directionally to all networks
-            if sc.isnumber(self.pars.beta):
-                β = self.pars.beta
-                self.pars.beta = sc.objdict({k:[β,β] for k in networks.keys()})
-
-            # If beta is a dict, check all entries are bi-directional
-            elif isinstance(self.pars.beta, dict):
-                for k,β in self.pars.beta.items():
-                    if sc.isnumber(β):
-                        self.pars.beta[k] = [β,β]
-        return
 
     @property
     def infectious(self):
@@ -204,6 +186,7 @@ def init_post(self):
         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.
         """
+        super().init_post()
         if self.pars.init_prev is None:
             return
 
@@ -223,33 +206,52 @@ def init_results(self):
             ss.Result(self.name, 'cum_infections', sim.npts, dtype=int, scale=True, label='Cumulative infections'),
         ]
         return
-
-    def _check_betas(self):
-        """ Check that there's a network for each beta key """
-        # Ensure keys are lowercase
-        if isinstance(self.pars.beta, dict): # TODO: check if needed
-            self.pars.beta = {k.lower(): v for k, v in self.pars.beta.items()}
-
-        # Create a mapping between beta and networks, and populate it
-        betapars = self.pars.beta
-        betamap = sc.objdict()
-        netkeys = list(self.sim.networks.keys())
-        if netkeys: # Skip if no networks
-            for bkey in betapars.keys():
-                orig_bkey = bkey[:]
-                if bkey in netkeys: # TODO: CK: could tidy up logic
-                    betamap[bkey] = betapars[orig_bkey]
+        
+    def validate_beta(self):
+        """ Validate beta and return as a map to match the networks """
+        sim = self.sim
+        
+        if 'beta' not in self.pars:
+            errormsg = f'Disease {self.name} is missing beta; pars are: {sc.strjoin(self.pars.keys())}'
+            raise sc.KeyNotFoundError(errormsg)
+
+        # If beta is a scalar, apply this bi-directionally to all networks
+        β = self.pars.beta
+        if sc.isnumber(β):
+            betamap = {ss.standardize_netkey(k):[β,β] for k in sim.networks.keys()}
+
+        # If beta is a dict, check all entries are bi-directional
+        elif isinstance(β, dict):
+            betamap = dict()
+            for k,thisbeta in β.items():
+                nkey = ss.standardize_netkey(k)
+                if sc.isnumber(thisbeta):
+                    betamap[nkey] = [thisbeta, thisbeta]
                 else:
-                    if 'net' not in bkey:
-                        bkey += 'net'  # Add 'net' suffix if not already there
-                    if bkey in netkeys:
-                        betamap[bkey] = betapars[orig_bkey]
-                    else:
-                        errormsg = f'No network for beta parameter "{bkey}"; your beta should match network keys:\n{sc.newlinejoin(netkeys)}'
-                        raise ValueError(errormsg)
+                    betamap[nkey] = thisbeta
+        
+        # Check that it matches the network
+        netkeys = [ss.standardize_netkey(k) for k in list(sim.networks.keys())]
+        if set(betamap.keys()) != set(netkeys):
+            errormsg = f'Network keys ({netkeys}) and beta keys ({betamap.keys()}) do not match'
+            raise ValueError(errormsg)
+
         return betamap
+    
+    def step(self):
+        """
+        Perform key infection updates, including infection and setting prognoses
+        """
+        # Create new cases
+        new_cases, sources, networks = self.infect() # TODO: store outputs in self or use objdict rather than 3 returns
+        
+        # Set prognoses
+        if len(new_cases):
+            self.set_outcomes(new_cases, sources)
+        
+        return new_cases, sources, networks
 
-    def make_new_cases(self):
+    def infect(self):
         """
         Add new cases of module, through transmission, incidence, etc.
         
@@ -259,68 +261,58 @@ def make_new_cases(self):
         new_cases = []
         sources = []
         networks = []
-        betamap = self._check_betas()
-
+        betamap = self.validate_beta()
+        
+        rel_trans = self.rel_trans.asnew(self.infectious * self.rel_trans)
+        rel_sus   = self.rel_sus.asnew(self.susceptible * self.rel_sus)
+        
         for i, (nkey,net) in enumerate(self.sim.networks.items()):
-            if not len(net):
-                break
-
-            nbetas = betamap[nkey]
-            edges = net.edges
-
-            rel_trans = self.rel_trans.asnew(self.infectious * self.rel_trans)
-            rel_sus   = self.rel_sus.asnew(self.susceptible * self.rel_sus)
-            p1p2b0 = [edges.p1, edges.p2, nbetas[0]]
-            p2p1b1 = [edges.p2, edges.p1, nbetas[1]]
-            for src, trg, beta in [p1p2b0, p2p1b1]:
-
-                # Skip networks with no transmission
-                if beta == 0:
-                    continue
-
-                # Calculate probability of a->b transmission.
-                beta_per_dt = net.beta_per_dt(disease_beta=beta, dt=self.sim.dt)
-                p_transmit = rel_trans[src] * rel_sus[trg] * beta_per_dt
-
-                # Generate a new random number based on the two other random numbers
-                rvs_s = self.rng_source.rvs(src)
-                rvs_t = self.rng_target.rvs(trg)
-                rvs = ss.combine_rands(rvs_s, rvs_t)
-                
-                new_cases_bool = rvs < p_transmit
-                new_cases.append(trg[new_cases_bool])
-                sources.append(src[new_cases_bool])
-                networks.append(np.full(np.count_nonzero(new_cases_bool), dtype=ss_int_, fill_value=i))
+            nk = ss.standardize_netkey(nkey) # TEMP
+            if len(net): # Skip networks with no edges
+                edges = net.edges
+                p1p2b0 = [edges.p1, edges.p2, betamap[nk][0]] # Person 1, person 2, beta 0
+                p2p1b1 = [edges.p2, edges.p1, betamap[nk][1]] # Person 2, person 1, beta 1
+                for src, trg, beta in [p1p2b0, p2p1b1]:
+                    if beta: # Skip networks with no transmission
+    
+                        # Calculate probability of a->b transmission.
+                        beta_per_dt = net.beta_per_dt(disease_beta=beta, dt=self.sim.dt)
+                        p_transmit = rel_trans[src] * rel_sus[trg] * beta_per_dt
+        
+                        # Generate a new random number based on the two other random numbers
+                        randvals = self.trans_rng.rvs(src, trg)
+                        transmitted = p_transmit > randvals
+                        target_uids = trg[transmitted]
+                        source_uids = src[transmitted]
+                        new_cases.append(target_uids)
+                        sources.append(source_uids)
+                        networks.append(np.full(len(target_uids), dtype=ss_int_, fill_value=i))
                 
-        # Tidy up
+        # Finalize
         if len(new_cases) and len(sources):
             new_cases = ss.uids.cat(new_cases)
             new_cases, inds = new_cases.unique(return_index=True)
             sources = ss.uids.cat(sources)[inds]
             networks = np.concatenate(networks)[inds]
         else:
-            new_cases = np.empty(0, dtype=int)
-            sources = np.empty(0, dtype=int)
-            networks = np.empty(0, dtype=int)
-
-        if len(new_cases):
-            self._set_cases(new_cases, sources)
-            
+            new_cases = ss.uids()
+            sources = ss.uids()
+            networks = np.empty(0, dtype=ss_int_)
         return new_cases, sources, networks
 
-    def _set_cases(self, target_uids, source_uids=None):
+    def set_outcomes(self, uids, sources=None):
         sim = self.sim
-        congenital = sim.people.age[target_uids] <= 0
+        congenital = sim.people.age[uids] <= 0
         if np.count_nonzero(congenital):
-            src_c = source_uids[congenital] if source_uids is not None else None
-            self.set_congenital(target_uids[congenital], src_c)
-        src_p = source_uids[~congenital] if source_uids is not None else None
-        self.set_prognoses(target_uids[~congenital], src_p)
+            src_c = sources[congenital] if sources is not None else None
+            self.set_congenital(uids[congenital], src_c)
+        src_p = sources[~congenital] if sources is not None else None
+        self.set_prognoses(uids[~congenital], src_p)
         return
 
-    def set_congenital(self, target_uids, source_uids=None):
+    def set_congenital(self, uids, sources=None):
         pass
-
+    
     def update_results(self):
         super().update_results()
         res = self.results
@@ -359,9 +351,14 @@ class InfectionLog(nx.MultiDiGraph):
 
     A table of outcomes can be returned using `InfectionLog.line_list()`
     """
-
-    # Add entries
-    # Add items to the most recent infection for an agent
+    def add_entries(self, sim, uids, sources=None): # TODO: reconcile with other methods
+        if sources is None:
+            for target in uids:
+                self.log.append(np.nan, target, sim.year)
+        else:
+            for target, source in zip(uids, sources):
+                self.log.append(source, target, sim.year)
+        return
 
     def add_data(self, uids, **kwargs):
         """
@@ -370,8 +367,9 @@ def add_data(self, uids, **kwargs):
         This method can be used to add data to an existing transmission event.
         The most recent transmission event will be used
 
-        :param uid: The UID of the target node (the agent that was infected)
-        :param kwargs: Remaining arguments are stored as edge data
+        Args:
+            uid: The UID of the target node (the agent that was infected)
+            kwargs: Remaining arguments are stored as edge data
         """
         for uid in sc.promotetoarray(uids):
             source, target, key = max(self.in_edges(uid, keys=True),
@@ -384,7 +382,6 @@ def append(self, source, target, t, **kwargs):
         self.add_edge(source, target, key=t, **kwargs)
         return
 
-    @property
     def line_list(self):
         """
         Return a tabular representation of the log
diff --git a/starsim/diseases/cholera.py b/starsim/diseases/cholera.py
index ac8060be..9fead94e 100644
--- a/starsim/diseases/cholera.py
+++ b/starsim/diseases/cholera.py
@@ -17,7 +17,7 @@ class Cholera(ss.Infection):
     def __init__(self, pars=None, *args, **kwargs):
         """ Initialize with parameters """
         super().__init__()
-        self.default_pars(
+        self.define_pars(
             # Initial conditions and beta
             beta = 1.0, # Placeholder value
             init_prev = ss.bernoulli(0.005),
@@ -41,7 +41,7 @@ def __init__(self, pars=None, *args, **kwargs):
         self.update_pars(pars, **kwargs)
 
         # Boolean states
-        self.add_states(
+        self.define_states(
             # Susceptible & infected are added automatically, here we add the rest
             ss.BoolArr('exposed', label='Exposed'),
             ss.BoolArr('symptomatic', label='Symptomatic'),
@@ -77,7 +77,7 @@ def init_results(self):
         ]
         return
 
-    def update_pre(self):
+    def step_state(self):
         """
         Adapted from https://github.com/optimamodel/gavi-outbreaks/blob/main/stisim/gavi/cholera.py
         Original version by Dom Delport
@@ -173,7 +173,7 @@ def make_new_cases(self):
             self.set_prognoses(new_cases, source_uids=None)
         return
 
-    def update_death(self, uids):
+    def step_die(self, uids):
         """ Reset infected/recovered flags for dead agents """
         for state in ['susceptible', 'exposed', 'infected', 'symptomatic', 'recovered']:
             self.statesdict[state][uids] = False
diff --git a/starsim/diseases/ebola.py b/starsim/diseases/ebola.py
index da32cdb3..aaffcade 100644
--- a/starsim/diseases/ebola.py
+++ b/starsim/diseases/ebola.py
@@ -16,7 +16,7 @@ class Ebola(SIR):
     def __init__(self, pars=None, *args, **kwargs):
         """ Initialize with parameters """
         super().__init__()
-        self.default_pars(
+        self.define_pars(
             # Initial conditions and beta
             init_prev       = ss.bernoulli(p=0.005),
             beta            = 1.0, # Placeholder value
@@ -37,7 +37,7 @@ def __init__(self, pars=None, *args, **kwargs):
         self.update_pars(pars=pars, **kwargs)
         
         # Boolean states
-        self.add_states(
+        self.define_states(
             # SIR are added automatically, here we add E
             ss.BoolArr('exposed', label='Exposed'),
             ss.BoolArr('severe', label='Severe'),
@@ -57,7 +57,7 @@ def __init__(self, pars=None, *args, **kwargs):
     def infectious(self):
         return self.infected | self.exposed
 
-    def update_pre(self):
+    def step_state(self):
 
         # Progress exposed -> infected
         ti = self.sim.ti
@@ -131,7 +131,7 @@ def set_prognoses(self, uids, source_uids=None):
         self.rel_trans[unburied_uids] = self.pars['unburied_factor']  # Change for unburied
         return
 
-    def update_death(self, uids):
+    def step_die(self, uids):
         # Reset infected/recovered flags for dead agents
         for state in ['susceptible', 'exposed', 'infected', 'severe', 'recovered']:
             self.statesdict[state][uids] = False
diff --git a/starsim/diseases/gonorrhea.py b/starsim/diseases/gonorrhea.py
index 7b4f2da2..3202aaa2 100644
--- a/starsim/diseases/gonorrhea.py
+++ b/starsim/diseases/gonorrhea.py
@@ -13,7 +13,7 @@ class Gonorrhea(ss.Infection):
     def __init__(self, pars=None, *args, **kwargs):
         # Parameters
         super().__init__()
-        self.default_pars(
+        self.define_pars(
             beta = 1.0, # Placeholder value
             dur_inf_in_days = ss.lognorm_ex(mean=10, stdev=0.6),  # median of 10 days (IQR 7–15 days) https://sti.bmj.com/content/96/8/556
             p_symp    = ss.bernoulli(p=0.5),  # Share of infections that are symptomatic. Placeholder value
@@ -25,7 +25,7 @@ def __init__(self, pars=None, *args, **kwargs):
         # States additional to the default disease states (see base class)
         # Additional states dependent on parameter values, e.g. self.p_symp?
         # These might be useful for connectors to target, e.g. if HIV reduces p_clear
-        self.add_states(
+        self.define_states(
             ss.BoolArr('symptomatic', label='Symptomatic'),
             ss.FloatArr('ti_clearance', label='Time of clearance'),
             ss.FloatArr('p_symp', default=1, label='Probability of symptoms'),
@@ -46,7 +46,7 @@ def update_results(self):
         self.results.new_clearances[ti] = np.count_nonzero(self.ti_clearance == ti)
         return
 
-    def update_pre(self):
+    def step_state(self):
         """ Natural clearance """
         clearances = self.ti_clearance <= self.sim.ti
         self.susceptible[clearances] = True
diff --git a/starsim/diseases/hiv.py b/starsim/diseases/hiv.py
index 706e2da0..b35ce610 100644
--- a/starsim/diseases/hiv.py
+++ b/starsim/diseases/hiv.py
@@ -13,7 +13,7 @@ class HIV(ss.Infection):
 
     def __init__(self, pars=None, *args, **kwargs):
         super().__init__()
-        self.default_pars(
+        self.define_pars(
             beta = 1.0, # Placeholder value
             cd4_min = 100,
             cd4_max = 500,
@@ -27,7 +27,7 @@ def __init__(self, pars=None, *args, **kwargs):
         self.update_pars(pars=pars, **kwargs)
 
         # States
-        self.add_states(
+        self.define_states(
             ss.BoolArr('on_art', label='On ART'),
             ss.FloatArr('ti_art', label='Time of ART initiation'),
             ss.FloatArr('ti_dead', label='Time of death'), # Time of HIV-caused death
@@ -42,7 +42,7 @@ def death_prob_func(module, sim, uids):
         out = np.array(out)
         return out
 
-    def update_pre(self):
+    def step_state(self):
         """ Update CD4 """
         people = self.sim.people
         self.cd4[people.alive & self.infected & self.on_art] += (self.pars.cd4_max - self.cd4[people.alive & self.infected & self.on_art])/self.pars.cd4_rate
@@ -89,7 +89,7 @@ def __init__(self, year, coverage, pars=None, **kwargs):
         self.year = sc.toarray(year)
         self.coverage = sc.toarray(coverage)
         super().__init__()
-        self.default_pars(
+        self.define_pars(
             art_delay = ss.constant(v=1) # Value in years
         )
         self.update_pars(pars=pars, **kwargs)
@@ -104,7 +104,8 @@ def init_pre(self, sim):
         self.initialized = True
         return
 
-    def apply(self, sim):
+    def step(self):
+        sim = self.sim
         if sim.year < self.year[0]:
             return
 
@@ -140,6 +141,7 @@ def init_pre(self, sim):
         self.cd4 = np.zeros((sim.npts, sim.people.n), dtype=int)
         return
 
-    def apply(self, sim):
+    def step(self):
+        sim = self.sim
         self.cd4[sim.t] = sim.people.hiv.cd4
         return
diff --git a/starsim/diseases/measles.py b/starsim/diseases/measles.py
index 7ff86d6e..f6b6d35c 100644
--- a/starsim/diseases/measles.py
+++ b/starsim/diseases/measles.py
@@ -15,7 +15,7 @@ class Measles(SIR):
     def __init__(self, pars=None, *args, **kwargs):
         """ Initialize with parameters """
         super().__init__()
-        self.default_pars(
+        self.define_pars(
             # Initial conditions and beta
             beta = 1.0, # Placeholder value
             init_prev = ss.bernoulli(p=0.005),
@@ -28,7 +28,7 @@ def __init__(self, pars=None, *args, **kwargs):
         self.update_pars(pars=pars, **kwargs)
 
         # SIR are added automatically, here we add E
-        self.add_states(
+        self.define_states(
             ss.BoolArr('exposed', label='Exposed'),
             ss.FloatArr('ti_exposed', label='Time of exposure'),
         )
@@ -39,7 +39,7 @@ def __init__(self, pars=None, *args, **kwargs):
     def infectious(self):
         return self.infected | self.exposed
 
-    def update_pre(self):
+    def step_state(self):
         # Progress exposed -> infected
         ti = self.sim.ti
         infected = (self.exposed & (self.ti_infected <= ti)).uids
@@ -85,7 +85,7 @@ def set_prognoses(self, uids, source_uids=None):
 
         return
 
-    def update_death(self, uids):
+    def step_die(self, uids):
         # Reset infected/recovered flags for dead agents
         for state in ['susceptible', 'exposed', 'infected', 'recovered']:
             self.statesdict[state][uids] = False
diff --git a/starsim/diseases/ncd.py b/starsim/diseases/ncd.py
index 6177f130..9c301359 100644
--- a/starsim/diseases/ncd.py
+++ b/starsim/diseases/ncd.py
@@ -19,14 +19,14 @@ class NCD(ss.Disease):
     """
     def __init__(self, pars=None, **kwargs):
         super().__init__()
-        self.default_pars(
+        self.define_pars(
             initial_risk = ss.bernoulli(p=0.3), # Initial prevalence of risk factors
             dur_risk = ss.expon(scale=10),
             prognosis = ss.weibull(c=2, scale=5), # Time in years between first becoming affected and death
         )
         self.update_pars(pars=pars, **kwargs)
         
-        self.add_states(
+        self.define_states(
             ss.BoolArr('at_risk', label='At risk'),
             ss.BoolArr('affected', label='Affected'),
             ss.FloatArr('ti_affected', label='Time of becoming affected'),
@@ -51,7 +51,7 @@ def init_post(self):
         self.ti_affected[initial_risk] = self.sim.ti + sc.randround(self.pars['dur_risk'].rvs(initial_risk) / self.sim.dt)
         return initial_risk
 
-    def update_pre(self):
+    def step_state(self):
         ti = self.sim.ti
         deaths = (self.ti_dead == ti).uids
         self.sim.people.request_death(deaths)
@@ -59,7 +59,7 @@ def update_pre(self):
         self.results.new_deaths[ti] = len(deaths) # Log deaths attributable to this module
         return
 
-    def make_new_cases(self):
+    def step(self):
         ti = self.sim.ti
         new_cases = (self.ti_affected == ti).uids
         self.affected[new_cases] = True
diff --git a/starsim/diseases/sir.py b/starsim/diseases/sir.py
index 939909cc..9c0e43e8 100644
--- a/starsim/diseases/sir.py
+++ b/starsim/diseases/sir.py
@@ -18,24 +18,67 @@ class SIR(ss.Infection):
     infected/infectious, and recovered. It also includes deaths, and basic
     results.
     """
+    def __future_init__(self, pars=None, **kwargs):
+        """ The plan for implementing events """
+        super().__init__()
+        self.define_pars(
+            beta = 0.1,
+            init_prev = ss.bernoulli(p=0.01),
+            dur_inf = ss.lognorm_ex(mean=6),
+            p_death = ss.bernoulli(p=0.01),
+        )
+        self.update_pars(pars, **kwargs)
+        
+        self.define_states(
+            ss.State('susceptible', default=True),
+            ss.State('infected'),
+            ss.State('recovered'),
+        )
+        
+        # Option 1 for defining events -- using objects
+        self.define_events(
+            ss.Event(self.susceptible, self.infected, func=self.infect,  reskey='infections'),
+            ss.Event(self.infected, self.recovered,   func=self.recover, reskey='recoveries'),
+            ss.Event(self.infected, self.dead,        func=self.step_die,     reskey='deaths'),
+        )
+        
+        # Option 2 for defining events -- using strings
+        self.define_events(
+            ss.Event('susceptible -> infected', self.infect,  reskey='infections'),
+            ss.Event('infected -> recovered',   self.recover, reskey='recoveries'),
+            ss.Event('infected -> dead',        self.step_die,     reskey='deaths'),
+        )
+        
+        self.define_attrs(
+            ss.FloatArr('rel_sus',   default=1.0),
+            ss.FloatArr('rel_trans', default=1.0),
+        )
+        return
+    
     def __init__(self, pars=None, **kwargs):
+        """ The current implementation """
         super().__init__()
-        self.default_pars(
+        self.define_pars(
             beta = 0.1,
             init_prev = ss.bernoulli(p=0.01),
             dur_inf = ss.lognorm_ex(mean=6),
             p_death = ss.bernoulli(p=0.01),
         )
         self.update_pars(pars, **kwargs)
-
-        self.add_states(
+        
+        self.define_states(
+            ss.BoolArr('susceptible', default=True, label='Susceptible'),
+            ss.BoolArr('infected', label='Infectious'),
             ss.BoolArr('recovered', label='Recovered'),
+            ss.FloatArr('ti_infected', label='Time of infection'),
             ss.FloatArr('ti_recovered', label='Time of recovery'),
             ss.FloatArr('ti_dead', label='Time of death'),
+            ss.FloatArr('rel_sus', default=1.0, label='Relative susceptibility'),
+            ss.FloatArr('rel_trans', default=1.0, label='Relative transmission'),
         )
         return
 
-    def update_pre(self):
+    def step_state(self):
         # Progress infectious -> recovered
         sim = self.sim
         recovered = (self.infected & (self.ti_recovered <= sim.ti)).uids
@@ -50,7 +93,6 @@ def update_pre(self):
 
     def set_prognoses(self, uids, source_uids=None):
         """ Set prognoses """
-        super().set_prognoses(uids, source_uids)
         ti = self.sim.ti
         dt = self.sim.dt
         self.susceptible[uids] = False
@@ -71,7 +113,7 @@ def set_prognoses(self, uids, source_uids=None):
         self.ti_recovered[rec_uids] = ti + dur_inf[~will_die] / dt
         return
 
-    def update_death(self, uids):
+    def step_die(self, uids):
         """ Reset infected/recovered flags for dead agents """
         self.susceptible[uids] = False
         self.infected[uids] = False
@@ -102,7 +144,7 @@ class SIS(ss.Infection):
     """
     def __init__(self, pars=None, *args, **kwargs):
         super().__init__()
-        self.default_pars(
+        self.define_pars(
             beta = 0.05,
             init_prev = ss.bernoulli(p=0.01),
             dur_inf = ss.lognorm_ex(mean=10),
@@ -111,13 +153,13 @@ def __init__(self, pars=None, *args, **kwargs):
         )
         self.update_pars(pars=pars, *args, **kwargs)
 
-        self.add_states(
+        self.define_states(
             ss.FloatArr('ti_recovered'),
             ss.FloatArr('immunity', default=0.0),
         )
         return
 
-    def update_pre(self):
+    def step_state(self):
         """ Progress infectious -> recovered """
         recovered = (self.infected & (self.ti_recovered <= self.sim.ti)).uids
         self.infected[recovered] = False
@@ -188,7 +230,7 @@ class sir_vaccine(ss.Vx):
     """
     def __init__(self, pars=None, *args, **kwargs):
         super().__init__()
-        self.default_pars(
+        self.define_pars(
             efficacy = 0.9,
             leaky = True
         )
diff --git a/starsim/diseases/syphilis.py b/starsim/diseases/syphilis.py
index badfb335..fd212c38 100644
--- a/starsim/diseases/syphilis.py
+++ b/starsim/diseases/syphilis.py
@@ -12,10 +12,10 @@
 
 class Syphilis(ss.Infection):
 
-    def __init__(self, pars=None, **kwargs):
+    def __future_init__(self, pars=None, **kwargs):
         # Parameters
         super().__init__()
-        self.default_pars(
+        self.define_pars(
             # Initial conditions
             beta = 1.0, # Placeholder
             init_prev = ss.bernoulli(p=0.03),
@@ -54,7 +54,82 @@ def __init__(self, pars=None, **kwargs):
         )
         self.update_pars(pars, **kwargs)
 
-        self.add_states(
+        self.define_states(
+            # Adult syphilis states
+            ss.State('susceptible', label='Susceptible', default=True),
+            ss.State('exposed', label='Exposed'),  # AKA incubating. Free of symptoms, not transmissible
+            ss.State('primary', label='Primary'),  # Primary chancres
+            ss.State('secondary', label='Secondary'),  # Inclusive of those who may still have primary chancres
+            ss.State('latent_temp', label="Latent temporary"),  # Relapses to secondary (~1y)
+            ss.State('latent_long', label="Latent long"),  # Can progress to tertiary or remain here
+            ss.State('tertiary', label="Tertiary"),  # Includes complications (cardio/neuro/disfigurement)
+            ss.State('immune', label="Immune"),  # After effective treatment people may acquire temp immunity
+            ss.State('ever_exposed', label="Ever exposed", track_time=False),  # Anyone ever exposed - stays true after treatment
+            ss.State('congenital'),  # Congenital syphilis states
+        )
+        
+        self.define_events(
+            ss.Event(src='susceptible', dest=['exposed', 'infected', 'ever_exposed'], func=self.infect),
+            ss.Event('exposed -> primary', func=self.to_primary),
+            ss.Event('primary -> secondary', func=self.to_secondary),
+            ss.Event('latent_temp -> secondary', func=self.to_secondary_latent),
+            ss.Event('secondary -> latent_temp', func=self.to_latent_temp),
+            ss.Event('secondary -> latent_long', func=self.to_latent_long),
+            ss.Event('latent_long -> tertiary', func=self.to_tertiary),
+        )
+    
+        # Timestep of state changes -- not all will be needed in future
+        self.define_attrs(
+            ss.FloatArr('ti_miscarriage', label='Time of miscarriage'),
+            ss.FloatArr('ti_nnd', label='Time of neonatal death'),
+            ss.FloatArr('ti_stillborn', label='Time of stillborn'),
+            ss.FloatArr('ti_congenital', label='Time of congenital syphilis'),
+        )
+        return
+    
+    def __init__(self, pars=None, **kwargs):
+        # Parameters
+        super().__init__()
+        self.define_pars(
+            # Initial conditions
+            beta = 1.0, # Placeholder
+            init_prev = ss.bernoulli(p=0.03),
+            
+            # Adult syphilis natural history, all specified in years
+            dur_exposed = ss.lognorm_ex(mean=1 / 12, stdev=1 / 36),  # https://pubmed.ncbi.nlm.nih.gov/9101629/
+            dur_primary = ss.lognorm_ex(mean=1.5 / 12, stdev=1 / 36),  # https://pubmed.ncbi.nlm.nih.gov/9101629/
+            dur_secondary = ss.normal(loc=3.6 / 12, scale=1.5 / 12),  # https://pubmed.ncbi.nlm.nih.gov/9101629/
+            dur_latent_temp = ss.lognorm_ex(mean=1, stdev=6 / 12),  # https://pubmed.ncbi.nlm.nih.gov/9101629/
+            dur_latent_long = ss.lognorm_ex(mean=20, stdev=8),  # https://pubmed.ncbi.nlm.nih.gov/9101629/
+            p_latent_temp = ss.bernoulli(p=0.25),  # https://pubmed.ncbi.nlm.nih.gov/9101629/
+            p_tertiary = ss.bernoulli(p=0.35),  # https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4917057/
+    
+            # Transmission by stage
+            rel_trans = dict(
+                exposed=1,
+                primary=1,
+                secondary=1,
+                latent_temp=0.075,
+                latent_long=0.075,
+                tertiary=0.05,
+            ),
+    
+            # Congenital syphilis outcomes
+            # Birth outcomes coded as:
+            #   0: Neonatal death
+            #   1: Stillborn
+            #   2: Congenital syphilis
+            #   3: Live birth without syphilis-related complications
+            # Source: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5973824/)
+            birth_outcomes = sc.objdict(
+                active = ss.choice(a=5, p=np.array([0.125, 0.125, 0.20, 0.35, 0.200])), # Probabilities of active by birth outcome
+                latent = ss.choice(a=5, p=np.array([0.050, 0.075, 0.10, 0.05, 0.725])), # Probabilities of latent
+            ),
+            birth_outcome_keys = ['miscarriage', 'nnd', 'stillborn', 'congenital'],
+        )
+        self.update_pars(pars, **kwargs)
+    
+        self.define_states(
             # Adult syphilis states
             ss.BoolArr('exposed', label='Exposed'),  # AKA incubating. Free of symptoms, not transmissible
             ss.BoolArr('primary', label='Primary'),  # Primary chancres
@@ -79,7 +154,7 @@ def __init__(self, pars=None, **kwargs):
             ss.FloatArr('ti_stillborn', label='Time of stillborn'),
             ss.FloatArr('ti_congenital', label='Time of congenital syphilis'),
         )
-
+    
         return
 
     @property
@@ -118,7 +193,7 @@ def init_results(self):
         ]
         return
 
-    def update_pre(self):
+    def step_state(self):
         """ Updates prior to interventions """
 
         # Primary
@@ -195,7 +270,6 @@ def set_prognoses(self, uids, source_uids=None):
         """
         Set initial prognoses for adults newly infected with syphilis
         """
-        super().set_prognoses(uids, source_uids)
         
         ti = self.sim.ti
         dt = self.sim.dt
@@ -334,11 +408,12 @@ def _parse_product_str(self, product):
         else:
             return products[product]
 
-    def check_eligibility(self, sim):
+    def check_eligibility(self):
         """
         Return an array of indices of agents eligible for screening at time t, i.e. sexually active
         females in age range, plus any additional user-defined eligibility
         """
+        sim = self.sim
         if self.eligibility is not None:
             is_eligible = self.eligibility(sim)
         else:
@@ -374,7 +449,8 @@ def init_pre(self, sim):
         self.results += ss.Result('syphilis', 'n_tx', sim.npts, dtype=int, scale=True, label='Number treated')
         return
 
-    def apply(self, sim):
-        treat_inds = super().apply(sim)
+    def step(self):
+        sim = self.sim
+        treat_inds = super().step()
         sim.people.syphilis.infected[treat_inds] = False
         self.results['n_tx'][sim.ti] += len(treat_inds)
diff --git a/starsim/distributions.py b/starsim/distributions.py
index d3512a88..52d74b54 100644
--- a/starsim/distributions.py
+++ b/starsim/distributions.py
@@ -597,6 +597,7 @@ def plot_hist(self, n=1000, bins=None, fig_kw=None, hist_kw=None):
              'poisson', 'weibull', 'constant', 'randint', 'rand_raw', 'bernoulli',
              'choice', 'histogram']
 __all__ += dist_list
+__all__ += ['multi_random'] # Not a dist in the same sense as the others
 
 
 class random(Dist):
@@ -972,7 +973,64 @@ def __init__(self, values=None, bins=None, density=False, data=None, **kwargs):
         super().__init__(dist=dist, distname='histogram', **kwargs)
         self.dynamic_pars = False # Set to false since array arguments don't imply dynamic pars here
         return
+
+
+class multi_random(sc.prettyobj):
+    """
+    A class for holding two or more ss.random() distributions, and generating
+    random numbers linked to each of them. Useful for e.g. pairwise transmission
+    probabilities.
+    
+    See ss.combine_rands() for the manual version; in almost all cases this class
+    should be used instead.
+    
+    Usage:
+        multi = ss.multi_random('source', 'target')
+        rvs = multi.rvs(source_uids, target_uids)
+    """
+    def __init__(self, names, *args, **kwargs):
+        names = sc.mergelists(names, args)
+        self.dists = [ss.rand_raw(name=name, **kwargs) for name in names]
+        return
+    
+    def __len__(self):
+        return len(self.dists)
+        
+    def initialize(self, *args, **kwargs):
+        for dist in self.dists: dist.initialize(*args, **kwargs)
+        return
+    
+    def reset(self, *args, **kwargs):
+        for dist in self.dists: dist.reset(*args, **kwargs)
+        return
+    
+    def jump(self, *args, **kwargs):
+        for dist in self.dists: dist.jump(*args, **kwargs)
+        return
+    
+    def rvs(self, *args):
+        """ Get random variates from each of the underlying distributions and combine them efficiently """
+        # Validation
+        n_args = len(args)
+        n_dists = len(self)
+        if n_args != len(self):
+            errormsg = f'Number of UID lists supplied ({n_args}) does not match number of distributions ({n_dists})'
+            raise ValueError(errormsg)
+        
+        # Generate the random numbers
+        rvs_list = [dist.rvs(arg) for dist,arg in zip(self.dists, args)]
         
+        # Combine using bitwise-or
+        rand_ints = rvs_list[0]
+        for rand_ints2 in rvs_list[1:]:
+            rand_ints = np.bitwise_xor(rand_ints*rand_ints2, rand_ints-rand_ints2)
+            
+        # Normalize
+        rvs = rand_ints / np.iinfo(np.uint64).max
+        return rvs
+
+    
+
 
 #%% Dist exceptions
 
diff --git a/starsim/interventions.py b/starsim/interventions.py
index 4c405002..509b6758 100644
--- a/starsim/interventions.py
+++ b/starsim/interventions.py
@@ -1,60 +1,19 @@
 """
-Define interventions (and analyzers)
+Define interventions
 """
 
 import starsim as ss
 import sciris as sc
 import numpy as np
 
-__all__ = ['Plugin', 'Analyzer', 'Intervention']
+__all__ = ['Intervention']
 
 
-class Plugin(ss.Module):
-    """ Base class for interventions and analyzers """
-    def __init__(self, *args, **kwargs):
-        super().__init__(*args, **kwargs)
-        return
-    
-    def __call__(self, *args, **kwargs):
-        return self.apply(*args, **kwargs)
-    
-    def init_pre(self, sim):
-        return super().init_pre(sim)
-    
-    def apply(self, sim):
-        pass
-
-    def finalize(self):
-        return super().finalize()
-    
-    @classmethod
-    def from_func(cls, func):
-        """ Create an intervention or analyzer from a function """
-        name = func.__name__
-        new = cls(name=name)
-        new.apply = func
-        return new
-
-
-class Analyzer(Plugin):
-    """
-    Base class for analyzers. Analyzers are used to provide more detailed information 
-    about a simulation than is available by default -- for example, pulling states 
-    out of sim.people on a particular timestep before they get updated on the next step.
-    
-    The key method of the analyzer is ``apply()``, which is called with the sim
-    on each timestep.
-    
-    To retrieve a particular analyzer from a sim, use sim.get_analyzer().
-    """
-    pass
-
-
-class Intervention(Plugin):
+class Intervention(ss.Module):
     """
     Base class for interventions.
     
-    The key method of the intervention is ``apply()``, which is called with the sim
+    The key method of the intervention is ``step()``, which is called with the sim
     on each timestep.
     """
 
@@ -79,12 +38,12 @@ def _parse_product(self, product):
     def _parse_product_str(self, product):
         raise NotImplementedError
 
-    def check_eligibility(self, sim):
+    def check_eligibility(self):
         """
         Return an array of indices of agents eligible for screening at time t
         """
         if self.eligibility is not None:
-            is_eligible = self.eligibility(sim)
+            is_eligible = self.eligibility(self.sim)
             if is_eligible is not None and len(is_eligible): # Only worry if non-None/nonzero length
                 if isinstance(is_eligible, ss.BoolArr):
                     is_eligible = is_eligible.uids
@@ -92,7 +51,7 @@ def check_eligibility(self, sim):
                     errormsg = f'Eligibility function must return BoolArr or UIDs, not {type(is_eligible)} {is_eligible}'
                     raise TypeError(errormsg)
         else:
-            is_eligible = sim.people.auids # Everyone
+            is_eligible = self.sim.people.auids # Everyone
         return is_eligible
 
 
@@ -218,20 +177,21 @@ def init_pre(self, sim):
         self.outcomes = {k: np.array([], dtype=int) for k in self.product.hierarchy}
         return
 
-    def deliver(self, sim):
+    def deliver(self):
         """
         Deliver the diagnostics by finding who's eligible, finding who accepts, and applying the product.
         """
+        sim = self.sim
         ti = sc.findinds(self.timepoints, sim.ti)[0]
         prob = self.prob[ti]  # Get the proportion of people who will be tested this timestep
-        eligible_uids = self.check_eligibility(sim)  # Check eligibility
+        eligible_uids = self.check_eligibility()  # Check eligibility
         self.coverage_dist.set(p=prob)
         accept_uids = self.coverage_dist.filter(eligible_uids)
         if len(accept_uids):
-            self.outcomes = self.product.administer(sim, accept_uids)  # Actually administer the diagnostic
+            self.outcomes = self.product.administer(accept_uids)  # Actually administer the diagnostic
         return accept_uids
 
-    def check_eligibility(self, sim):
+    def check_eligibility(self):
         raise NotImplementedError
 
 
@@ -241,19 +201,20 @@ class BaseScreening(BaseTest):
     Args:
         kwargs (dict): passed to BaseTest
     """
-    def check_eligibility(self, sim):
+    def check_eligibility(self):
         """
         Check eligibility
         """
         raise NotImplementedError
 
-    def apply(self, sim, module=None):
+    def step(self):
         """
         Perform screening by finding who's eligible, finding who accepts, and applying the product.
         """
+        sim = self.sim
         accept_uids = ss.uids()
         if sim.ti in self.timepoints:
-            accept_uids = self.deliver(sim)
+            accept_uids = self.deliver()
             self.screened[accept_uids] = True
             self.screens[accept_uids] += 1
             self.ti_screened[accept_uids] = sim.ti
@@ -269,13 +230,13 @@ class BaseTriage(BaseTest):
     Args:
         kwargs (dict): passed to BaseTest
     """
-    def check_eligibility(self, sim):
-        return sc.promotetoarray(self.eligibility(sim))
+    def check_eligibility(self):
+        return sc.promotetoarray(self.eligibility(self.sim))
 
-    def apply(self, sim):
+    def step(self):
         self.outcomes = {k: np.array([], dtype=int) for k in self.product.hierarchy}
         accept_inds = ss.uids()
-        if sim.t in self.timepoints: accept_inds = self.deliver(sim)
+        if self.sim.t in self.timepoints: accept_inds = self.deliver() # TODO: not robust for timestep
         return accept_inds
 
 
@@ -332,7 +293,7 @@ class routine_triage(BaseTriage, RoutineDelivery):
 
     **Example**:
         # Example: Triage positive screens into confirmatory testing
-        screened_pos = lambda sim: sim.get_intervention('screening').outcomes['positive']
+        screened_pos = lambda sim: sim.interventions.screening.outcomes['positive']
         triage = ss.routine_triage(product=my_triage, eligibility=screen_pos, prob=0.9, start_year=2030)
     """
 
@@ -356,7 +317,7 @@ class campaign_triage(BaseTriage, CampaignDelivery):
 
     **Examples**:
         # Example: In 2030, triage all positive screens into confirmatory testing
-        screened_pos = lambda sim: sim.get_intervention('screening').outcomes['positive']
+        screened_pos = lambda sim: sim.interventions.screening.outcomes['positive']
         triage1 = ss.campaign_triage(product=my_triage, eligibility=screen_pos, prob=0.9, years=2030)
     """
 
@@ -399,33 +360,33 @@ def init_pre(self, sim):
         self.outcomes = {k: np.array([], dtype=int) for k in ['unsuccessful', 'successful']} # Store outcomes on each timestep
         return
 
-    def get_accept_inds(self, sim):
+    def get_accept_inds(self):
         """
         Get indices of people who will acccept treatment; these people are then added to a queue or scheduled for receiving treatment
         """
         accept_uids = ss.uids()
-        eligible_uids = self.check_eligibility(sim)  # Apply eligiblity
+        eligible_uids = self.check_eligibility()  # Apply eligiblity
         if len(eligible_uids):
             self.coverage_dist.set(p=self.prob[0])
             accept_uids = self.coverage_dist.filter(eligible_uids)
         return accept_uids
 
-    def get_candidates(self, sim):
+    def get_candidates(self):
         """
         Get candidates for treatment on this timestep. Implemented by derived classes.
         """
         raise NotImplementedError
 
-    def apply(self, sim):
+    def step(self):
         """
         Perform treatment by getting candidates, checking their eligibility, and then treating them.
         """
         # Get indices of who will get treated
-        treat_candidates = self.get_candidates(sim)  # NB, this needs to be implemented by derived classes
-        still_eligible = self.check_eligibility(sim)
+        treat_candidates = self.get_candidates()  # NB, this needs to be implemented by derived classes
+        still_eligible = self.check_eligibility()
         treat_uids = treat_candidates.intersect(still_eligible)
         if len(treat_uids):
-            self.outcomes = self.product.administer(sim, treat_uids)
+            self.outcomes = self.product.administer(treat_uids)
         return treat_uids
 
 
@@ -442,15 +403,15 @@ def __init__(self, max_capacity=None, **kwargs):
         self.max_capacity = max_capacity
         return
 
-    def add_to_queue(self, sim):
+    def add_to_queue(self):
         """
         Add people who are willing to accept treatment to the queue
         """
-        accept_inds = self.get_accept_inds(sim)
+        accept_inds = self.get_accept_inds()
         if len(accept_inds): self.queue += accept_inds.tolist()
         return
 
-    def get_candidates(self, sim):
+    def get_candidates(self):
         """
         Get the indices of people who are candidates for treatment
         """
@@ -462,13 +423,13 @@ def get_candidates(self, sim):
                 treat_candidates = self.queue[:self.max_capacity]
         return ss.uids(treat_candidates) # TODO: Check
 
-    def apply(self, sim):
+    def step(self):
         """
         Apply treatment. On each timestep, this method will add eligible people who are willing to accept treatment to a
         queue, and then will treat as many people in the queue as there is capacity for.
         """
-        self.add_to_queue(sim)
-        treat_inds = BaseTreatment.apply(self, sim) # Apply method from BaseTreatment class
+        self.add_to_queue()
+        treat_inds = BaseTreatment.step(self) # Apply method from BaseTreatment class
         self.queue = [e for e in self.queue if e not in treat_inds] # Recreate the queue, removing people who were treated
         return treat_inds
 
@@ -499,16 +460,17 @@ def __init__(self, product=None, prob=None, label=None, **kwargs):
         self.coverage_dist = ss.bernoulli(p=0)  # Placeholder
         return
 
-    def apply(self, sim):
+    def step(self):
         """
         Deliver the diagnostics by finding who's eligible, finding who accepts, and applying the product.
         """
+        sim = self.sim
         accept_uids = np.array([])
         if sim.ti in self.timepoints:
 
             ti = sc.findinds(self.timepoints, sim.ti)[0]
             prob = self.prob[ti]  # Get the proportion of people who will be tested this timestep
-            is_eligible = self.check_eligibility(sim)  # Check eligibility
+            is_eligible = self.check_eligibility()  # Check eligibility
             self.coverage_dist.set(p=prob)
             accept_uids = self.coverage_dist.filter(is_eligible)
 
diff --git a/starsim/modules.py b/starsim/modules.py
index e3b24234..4aec8122 100644
--- a/starsim/modules.py
+++ b/starsim/modules.py
@@ -1,11 +1,13 @@
 """
-General module class -- base class for diseases, interventions, etc.
+General module class -- base class for diseases, interventions, etc. Also
+defines Analyzers and Connectors.
 """
 
 import sciris as sc
 import starsim as ss
+from functools import partial
 
-__all__ = ['module_map', 'find_modules', 'Module', 'Connector']
+__all__ = ['module_map', 'find_modules', 'Module', 'Analyzer', 'Connector']
 
 
 def module_map(key=None):
@@ -39,13 +41,13 @@ def find_modules(key=None):
             except:
                 pass
     return modules if key is None else modules[key]
-
+    
 
 class Module(sc.quickobj):
 
-    def __init__(self, name=None, label=None, requires=None):
-        self.set_metadata(name, label, requires) # Usually reset as part of self.update_pars()
-        self.pars = ss.Pars() # Usually populated via self.default_pars()
+    def __init__(self, name=None, label=None):
+        self.set_metadata(name, label) # Usually reset as part of self.update_pars()
+        self.pars = ss.Pars() # Usually populated via self.define_pars()
         self.results = ss.Results(self.name)
         self.initialized = False
         self.finalized = False
@@ -55,14 +57,25 @@ def __bool__(self):
         """ Ensure that zero-length modules (e.g. networks) are still truthy """
         return True
     
-    def set_metadata(self, name, label, requires):
+    def __call__(self, *args, **kwargs):
+        """ Allow modules to be called like functions """
+        return self.step(*args, **kwargs)
+
+    def disp(self, output=False):
+        """ Display the full object """
+        out = sc.prepr(self)
+        if not output:
+            print(out)
+        else:
+            return out
+    
+    def set_metadata(self, name, label):
         """ Set metadata for the module """
         self.name = sc.ifelse(name, getattr(self, 'name', self.__class__.__name__.lower())) # Default name is the class name
         self.label = sc.ifelse(label, getattr(self, 'label', self.name))
-        self.requires = sc.mergelists(requires)
         return
     
-    def default_pars(self, inherit=True, **kwargs):
+    def define_pars(self, inherit=True, **kwargs): # TODO: think if inherit should default to true or false
         """ Create or merge Pars objects """
         if inherit: # Merge with existing
             self.pars.update(**kwargs, create=True)
@@ -82,7 +95,7 @@ def update_pars(self, pars, **kwargs):
         self.pars.update(matches)
                 
         # Update module attributes
-        metadata = {key:pars.pop(key, getattr(self, key, None)) for key in ['name', 'label', 'requires']}
+        metadata = {key:pars.pop(key, None) for key in ['name', 'label']}
         self.set_metadata(**metadata)
         
         # Should be no remaining pars
@@ -91,60 +104,49 @@ def update_pars(self, pars, **kwargs):
             raise ValueError(errormsg)
         return
     
-    def check_requires(self, sim):
-        """ Check that the module's requirements (of other modules) are met """
-        errs = sc.autolist()
-        all_names = [m.__class__ for m in sim.modules] + [m.name for m in sim.modules if hasattr(m, 'name')]
-        for req in self.requires:
-            if req not in all_names:
-                errs += req
-        if len(errs):
-            errormsg = f'{self.name} (label={self.label}) requires the following module(s), but the Sim does not contain them.'
-            errormsg += sc.newlinejoin(errs)
-            raise Exception(errormsg)
-        return
-
     def init_pre(self, sim):
         """
         Perform initialization steps
 
         This method is called once, as part of initializing a Sim. Note: after
-        initialization, initialized=False until init_post() is called (which is after
+        initialization, initialized=False until init_vals() is called (which is after
         distributions are initialized).
-        
-        Note: distributions cannot be used here because they aren't initialized 
-        until after init_pre() is called. Use init_post() instead.
         """
-        self.check_requires(sim)
         self.sim = sim # Link back to the sim object
         ss.link_dists(self, sim, skip=ss.Sim) # Link the distributions to sim and module
         sim.pars[self.name] = self.pars
         sim.results[self.name] = self.results
         sim.people.add_module(self) # Connect the states to the people
+        self.init_results()
         return
     
+    def init_results(self):
+        """ Initialize any results required; part of init_pre() """
+        pass
+    
     def init_post(self):
-        """ Initialize the values of the states, including calling distributions; the last step of initialization """
+        """ Initialize the values of the states; the last step of initialization """
         for state in self.states:
             if not state.initialized:
                 state.init_vals()
         self.initialized = True
         return
     
-    def disp(self, output=False):
-        """ Display the full object """
-        out = sc.prepr(self)
-        if not output:
-            print(out)
-        else:
-            return out
+    def step(self):
+        """ Define how the module updates over time """
+        pass
+    
+    def update_results(self):
+        """ Perform any results updates on each timestep """
+        pass
 
     def finalize(self):
+        """ Perform any final operations, such as removing unneeded data """
         self.finalize_results()
         self.finalized = True
         return
 
-    def finalize_results(self):
+    def finalize_results(self): # TODO: this is confusing, needs to be not redefined by the user, or called after a custom finalize_results()
         """ Finalize results """
         # Scale results
         for reskey, res in self.results.items():
@@ -152,7 +154,7 @@ def finalize_results(self):
                 self.results[reskey] = self.results[reskey]*self.sim.pars.pop_scale
         return
     
-    def add_states(self, *args, check=True):
+    def define_states(self, *args, check=True):
         """
         Add states to the module with the same attribute name as the state
         
@@ -209,6 +211,17 @@ def create(cls, name, *args, **kwargs):
         else:
             raise KeyError(f'Module "{name}" did not match any known Starsim modules')
             
+    @classmethod
+    def from_func(cls, func):
+        """ Create an module from a function """
+        def step(mod): # TODO: see if this can be done more simply
+            return mod.func(mod.sim)
+        name = func.__name__
+        mod = cls(name=name)
+        mod.func = func
+        mod.step = partial(step, mod)
+        return mod
+            
     def to_json(self):
         """ Export to a JSON-compatible format """
         out = sc.objdict()
@@ -230,12 +243,25 @@ def plot(self):
         return fig
 
 
+class Analyzer(Module):
+    """
+    Base class for Analyzers. Analyzers are used to provide more detailed information 
+    about a simulation than is available by default -- for example, pulling states 
+    out of sim.people on a particular timestep before they get updated on the next step.
+    
+    The key method of the analyzer is ``step()``, which is called with the sim
+    on each timestep.
+    
+    To retrieve a particular analyzer from a sim, use sim.get_analyzer().
+    """
+    pass
+
+
 class Connector(Module):
     """
-    Define a Connector, which mediates interactions between disease modules
+    Base class for Connectors, which mediate interactions between disease (or other) modules
     
     Because connectors can do anything, they have no specified structure: it is
     up to the user to define how they behave.    
     """
-    def update(self, sim):
-        pass
+    pass
diff --git a/starsim/network.py b/starsim/network.py
index 452bc5ae..8982bcff 100644
--- a/starsim/network.py
+++ b/starsim/network.py
@@ -60,9 +60,9 @@ class Network(ss.Module):
         network2 = ss.Network(**network, index=index, self_conn=self_conn, label=network.label)
     """
 
-    def __init__(self, key_dict=None, prenatal=False, postnatal=False, name=None, label=None, requires=None, **kwargs):
+    def __init__(self, key_dict=None, prenatal=False, postnatal=False, name=None, label=None, **kwargs):
         # Initialize as a module
-        super().__init__(name=name, label=label, requires=requires)
+        super().__init__(name=name, label=label)
 
         # Each relationship is characterized by these default set of keys, plus any user- or network-supplied ones
         default_keys = sc.objdict(
@@ -207,20 +207,20 @@ def pop_inds(self, inds):
         popped_inds = self.get_inds(inds, remove=True)
         return popped_inds
 
-    def append(self, contacts=None, **kwargs):
+    def append(self, edges=None, **kwargs):
         """
-        Append contacts to the current network.
+        Append edges to the current network.
 
         Args:
-            contacts (dict): a dictionary of arrays with keys p1,p2,beta, as returned from network.pop_inds()
+            edges (dict): a dictionary of arrays with keys p1,p2,beta, as returned from network.pop_inds()
         """
-        contacts = sc.mergedicts(contacts, kwargs)
+        edges = sc.mergedicts(edges, kwargs)
         for key in self.meta_keys():
             curr_arr = self.edges[key]
             try:
-                new_arr = contacts[key]
+                new_arr = edges[key]
             except KeyError:
-                errormsg = f'Cannot append contacts since required key "{key}" is missing'
+                errormsg = f'Cannot append edges since required key "{key}" is missing'
                 raise KeyError(errormsg)
             self.edges[key] = np.concatenate([curr_arr, new_arr])  # Resize to make room, preserving dtype
         self.validate_uids()
@@ -244,12 +244,12 @@ def from_df(self, df, keys=None):
             self.edges[key] = df[key].to_numpy()
         return self
 
-    def find_contacts(self, inds, as_array=True):
+    def find_edges(self, inds, as_array=True):
         """
-        Find all contacts of the specified people
+        Find all edges of the specified people
 
-        For some purposes (e.g. contact tracing) it's necessary to find all the contacts
-        associated with a subset of the people in this network. Since contacts are bidirectional
+        For some purposes (e.g. contact tracing) it's necessary to find all the edges
+        associated with a subset of the people in this network. Since edges are bidirectional
         it's necessary to check both p1 and p2 for the target indices. The return type is a Set
         so that there is no duplication of indices (otherwise if the Network has explicit
         symmetric interactions, they could appear multiple times). This is also for performance so
@@ -258,7 +258,7 @@ def find_contacts(self, inds, as_array=True):
         infection, e.g. exposure risk.
 
         Args:
-            inds (array): indices of people whose contacts to return
+            inds (array): indices of people whose edges to return
             as_array (bool): if true, return as sorted array (otherwise, return as unsorted set)
 
         Returns:
@@ -267,7 +267,7 @@ def find_contacts(self, inds, as_array=True):
         Example: If there were a network with
         - p1 = [1,2,3,4]
         - p2 = [2,3,1,4]
-        Then find_contacts([1,3]) would return {1,2,3}
+        Then find_edges([1,3]) would return {1,2,3}
         """
 
         # Check types
@@ -276,8 +276,8 @@ def find_contacts(self, inds, as_array=True):
         if inds.dtype != np.int64:  # pragma: no cover # This is int64 since indices often come from utils.true(), which returns int64
             inds = np.array(inds, dtype=np.int64)
 
-        # Find the contacts
-        contact_inds = ss.find_contacts(self.edges.p1, self.edges.p2, inds)
+        # Find the edges
+        contact_inds = ss.find_edges(self.edges.p1, self.edges.p2, inds)
         if as_array:
             contact_inds = np.fromiter(contact_inds, dtype=ss_int_)
             contact_inds.sort()
@@ -288,10 +288,6 @@ def add_pairs(self):
         """ Define how pairs of people are formed """
         pass
 
-    def update(self):
-        """ Define how pairs/connections evolve (in time) """
-        pass
-
     def remove_uids(self, uids):
         """
         Remove interactions involving specified UIDs
@@ -386,7 +382,7 @@ class StaticNet(Network):
     def __init__(self, graph=None, pars=None, **kwargs):
         super().__init__()
         self.graph = graph
-        self.default_pars(seed=True, p=None, n_contacts=10)
+        self.define_pars(seed=True, p=None, n_contacts=10)
         self.update_pars(pars, **kwargs)
         self.dist = ss.Dist(name='StaticNet')
         return
@@ -412,7 +408,7 @@ def init_post(self):
                 print(f"{str(e)}: networkx {self.graph.name} not supported. Try using ss.NullNet().")
                 raise e
         self.validate_pop(self.n_agents)
-        self.get_contacts()
+        self.get_edges()
         return
 
     def validate_pop(self, n_agents):
@@ -426,15 +422,15 @@ def validate_pop(self, n_agents):
             errmsg = f"The nx generator {self.graph.name} produced a graph with no edges"
             raise ValueError(errmsg)
 
-    def get_contacts(self):
+    def get_edges(self):
         p1s = []
         p2s = []
         for edge in self.graph.edges():
             p1, p2 = edge
             p1s.append(p1)
             p2s.append(p2)
-        contacts = dict(p1=p1s, p2=p2s, beta=np.ones_like(p1s))
-        self.append(contacts)
+        edges = dict(p1=p1s, p2=p2s, beta=np.ones_like(p1s))
+        self.append(edges)
         return
 
 
@@ -444,7 +440,7 @@ class RandomNet(DynamicNetwork):
     def __init__(self, pars=None, key_dict=None, **kwargs):
         """ Initialize """
         super().__init__(key_dict=key_dict)
-        self.default_pars(
+        self.define_pars(
             n_contacts = ss.constant(10),
             dur = 0,
         )
@@ -469,9 +465,9 @@ def get_source(inds, n_contacts):
             count += n
         return source
 
-    def get_contacts(self, inds, n_contacts):
+    def get_edges(self, inds, n_contacts):
         """
-        Efficiently generate contacts
+        Efficiently find edges
 
         Note that because of the shuffling operation, each person is assigned 2N contacts
         (i.e. if a person has 5 contacts, they appear 5 times in the 'source' array and 5
@@ -496,13 +492,13 @@ def get_contacts(self, inds, n_contacts):
         self.dist.jump() # Reset the RNG manually # TODO, think if there's a better way
         return source, target
 
-    def update(self):
+    def step(self):
         self.end_pairs()
         self.add_pairs()
         return
 
     def add_pairs(self):
-        """ Generate contacts """
+        """ Generate edges """
         people = self.sim.people
         born = people.alive & (people.age > 0)
         if isinstance(self.pars.n_contacts, ss.Dist):
@@ -512,7 +508,7 @@ def add_pairs(self):
 
         number_of_contacts = sc.randround(number_of_contacts / 2).astype(ss_int_)  # One-way contacts
 
-        p1, p2 = self.get_contacts(born.uids, number_of_contacts)
+        p1, p2 = self.get_edges(born.uids, number_of_contacts)
         beta = np.ones(len(p1), dtype=ss_float_)
 
         if isinstance(self.pars.dur, ss.Dist):
@@ -543,7 +539,7 @@ class ErdosRenyiNet(DynamicNetwork):
     def __init__(self, pars=None, key_dict=None, **kwargs):
         """ Initialize """
         super().__init__(key_dict=key_dict)
-        self.default_pars(
+        self.define_pars(
             p = 0.1, # Probability of each edge
             dur = 0, # Duration of zero ensures that new random edges are formed on each time step
         )
@@ -574,7 +570,7 @@ def add_pairs(self):
         # Use integers to create random numbers per edge
         i1 = ints[idx1]
         i2 = ints[idx2]
-        r = ss.combine_rands(i1, i2)
+        r = ss.utils.combine_rands(i1, i2) # TODO: use ss.multi_rand()
         edge = r <= self.pars.p
 
         p1 = idx1[edge]
@@ -605,12 +601,12 @@ class DiskNet(Network):
     def __init__(self, pars=None, key_dict=None, **kwargs):
         """ Initialize """
         super().__init__(key_dict=key_dict)
-        self.default_pars(
+        self.define_pars(
             r = 0.1, # Radius
             v = 0.05, # Velocity
         )
         self.update_pars(pars, **kwargs)
-        self.add_states(
+        self.define_states(
             ss.FloatArr('x', default=ss.random(), label='X position'),
             ss.FloatArr('y', default=ss.random(), label='Y position'),
             ss.FloatArr('theta', default=ss.uniform(high=2*np.pi), label='Heading'),
@@ -621,7 +617,7 @@ def init_post(self):
         self.add_pairs()
         return
 
-    def update(self):
+    def step(self):
 
         # Motion step
         vdt = self.pars.v * self.sim.dt
@@ -674,7 +670,7 @@ class NullNet(Network):
 
     Guarantees there's one (1) contact per agent (themselves), and that their connection weight is zero.
 
-    For an empty network (ie, no contacts) use
+    For an empty network (ie, no edges) use
     >> import starsim as ss
     >> import networkx as nx
     >> empty_net_static = ss.StaticNet(nx.empty_graph)
@@ -696,10 +692,10 @@ def init_pre(self, sim):
             if self.n > popsize:
                 errormsg = f'Please ensure the size of the network ({self.n} is less than or equal to the population size ({popsize}).'
                 raise ValueError(errormsg)
-        self.get_contacts()
+        self.get_edges()
         return
 
-    def get_contacts(self):
+    def get_edges(self):
         indices = np.arange(self.n)
         self.append(dict(p1=indices, p2=indices, beta=np.zeros_like(indices)))
         return
@@ -712,7 +708,7 @@ class MFNet(SexualNetwork):
     """
     def __init__(self, pars=None, key_dict=None, **kwargs):
         super().__init__(key_dict=key_dict)
-        self.default_pars(
+        self.define_pars(
             duration = ss.lognorm_ex(mean=15),  # 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 = ss.bernoulli(p=0.9),  # Probability of participating in this network - can vary by individual properties (age, sex, ...) using callable parameter values
             debut = ss.normal(loc=16),  # Age of debut can vary by using callable parameter values
@@ -789,7 +785,7 @@ def add_pairs(self):
 
         return len(p1)
 
-    def update(self):
+    def step(self):
         self.end_pairs()
         self.set_network_states(upper_age=self.sim.dt) # TODO: looks wrong
         self.add_pairs()
@@ -803,7 +799,7 @@ class MSMNet(SexualNetwork):
 
     def __init__(self, pars=None, key_dict=None, **kwargs):
         super().__init__(key_dict=key_dict)
-        self.default_pars(
+        self.define_pars(
             duration = ss.lognorm_ex(mean=15, stdev=15),
             debut = ss.normal(loc=16, scale=2),
             acts = ss.lognorm_ex(mean=80, stdev=20),
@@ -854,10 +850,10 @@ def add_pairs(self):
         
         return len(p1)
 
-    def update(self):
-        self.end_pairs(self.sim)
-        self.set_network_states(self.sim.people, upper_age=self.sim.dt)
-        self.add_pairs(self.sim)
+    def step(self):
+        self.end_pairs()
+        self.set_network_states()
+        self.add_pairs()
         return
 
 
@@ -874,7 +870,7 @@ def __init__(self, pars=None, **kwargs):
             male_shift is the average age that males are older than females in partnerships
         """
         super().__init__()
-        self.default_pars(
+        self.define_pars(
             inherit = True, # The MFNet already comes with pars, we want to keep those
             embedding_func = ss.normal(name='EmbeddingNet', loc=self.embedding_loc, scale=2),
             male_shift = 5,
@@ -905,11 +901,12 @@ def add_pairs(self):
 
         dist_mat = spsp.distance_matrix(loc_m[:, np.newaxis], loc_f[:, np.newaxis])
         ind_m, ind_f = spo.linear_sum_assignment(dist_mat)
+        n_pairs = len(ind_f)
 
         # Finalize pairs
         p1 = available_m[ind_m]
         p2 = available_f[ind_f]
-        beta = np.ones(len(p1)) # TODO: Allow custom beta
+        beta = np.ones(n_pairs) # TODO: Allow custom beta
         dur_vals = self.pars.duration.rvs(p1)
         act_vals = self.pars.acts.rvs(p1)
 
@@ -930,7 +927,7 @@ def __init__(self, key_dict=None, prenatal=True, postnatal=False, **kwargs):
         super().__init__(key_dict=key_dict, prenatal=prenatal, postnatal=postnatal, **kwargs)
         return
 
-    def update(self):
+    def step(self):
         """
         Set beta to 0 for women who complete duration of transmission
         Keep connections for now, might want to consider removing
@@ -941,9 +938,10 @@ def update(self):
 
     def end_pairs(self):
         people = self.sim.people
-        active = (self.edges.end > self.sim.ti) & people.alive[self.edges.p1] & people.alive[self.edges.p2]
+        edges = self.edges
+        active = (edges.end > self.sim.ti) & people.alive[edges.p1] & people.alive[edges.p2]
         for k in self.meta_keys():
-            self.edges[k] = self.edges[k][active]
+            edges[k] = edges[k][active]
         return len(active)
 
     def add_pairs(self, mother_inds=None, unborn_inds=None, dur=None, start=None):
diff --git a/starsim/parameters.py b/starsim/parameters.py
index 87ee237f..193c816b 100644
--- a/starsim/parameters.py
+++ b/starsim/parameters.py
@@ -181,7 +181,7 @@ def __init__(self, **kwargs):
         self.death_rate = None
         self.use_aging  = None # True if demographics, false otherwise
 
-        # Modules: demographics, diseases, connectors, networks, analyzers, and interventions
+        # Modules: demographics, diseases, networks, analyzers, and interventions
         self.people = None
         self.networks      = ss.ndict()
         self.demographics  = ss.ndict()
@@ -422,8 +422,8 @@ def convert_modules(self):
                             mod = expected_cls.from_func(mod)
                     
                     # Do final check
-                    if not isinstance(mod, expected_cls):
-                        errormsg = f'Was expecting {modkey} entry {i} to be class {expected_cls}, but was {type(mod)} instead'
+                    if not isinstance(mod, (expected_cls, ss.Module)): # TEMP: check if this check still works?
+                        errormsg = f'Was expecting {modkey} entry {i} to be class {expected_cls} or Plugin, but was {type(mod)} instead'
                         raise TypeError(errormsg)
                     modlist[i] = mod
         return
diff --git a/starsim/people.py b/starsim/people.py
index 513d56da..5e12d978 100644
--- a/starsim/people.py
+++ b/starsim/people.py
@@ -309,7 +309,7 @@ def request_death(self, uids):
         self.ti_dead[uids] = self.sim.ti
         return
 
-    def resolve_deaths(self):
+    def check_deaths(self):
         """ Carry out any deaths that took place this timestep """
         death_uids = (self.ti_dead <= self.sim.ti).uids
         self.alive[death_uids] = False
@@ -349,6 +349,12 @@ def update_results(self):
         res.cum_deaths[ti] = np.sum(res.new_deaths[:ti]) # TODO: inefficient to compute the cumulative sum on every timestep!
         return
     
+    def finish_step(self):
+        self.update_results()
+        self.remove_dead()
+        self.update_post()
+        return
+    
     def person(self, ind):
         """ Get all the properties for a single person """
         person = Person()
diff --git a/starsim/products.py b/starsim/products.py
index 252df413..0b800139 100644
--- a/starsim/products.py
+++ b/starsim/products.py
@@ -48,7 +48,7 @@ def __init__(self, df, hierarchy=None, *args, **kwargs):
     def default_value(self):
         return len(self.hierarchy) - 1
 
-    def administer(self, sim, uids, return_format='dict'):
+    def administer(self, uids, return_format='dict'):
         """
         Administer a testing product.
         
@@ -63,7 +63,7 @@ def administer(self, sim, uids, return_format='dict'):
 
         for disease in self.diseases:
             for state in self.health_states:
-                this_state = getattr(sim.diseases[disease], state)
+                this_state = getattr(self.sim.diseases[disease], state)
                 true_uids = this_state.uids # Find people for which this state is true
                 these_uids = true_uids.intersect(uids) # Find intersection of people in this state and the supplied UIDs
 
@@ -96,7 +96,7 @@ def __init__(self, df, *args, **kwargs):
         self.efficacy_dist = ss.bernoulli(p=0)
         return
 
-    def administer(self, sim, uids, return_format='dict'):
+    def administer(self, uids, return_format='dict'):
         """
         Loop over treatment states to determine those who are successfully treated and clear infection
         """
@@ -105,7 +105,7 @@ def administer(self, sim, uids, return_format='dict'):
 
         for disease_name in self.diseases:
 
-            disease = sim.diseases[disease_name]
+            disease = self.sim.diseases[disease_name]
 
             for state in self.health_states:
 
@@ -122,7 +122,7 @@ def administer(self, sim, uids, return_format='dict'):
                     self.efficacy_dist.set(p=thisdf.efficacy.values[0])
 
                     # HACK to reset the efficacy_dist as it is called multiple times per timestep. TODO: Refactor
-                    self.efficacy_dist.jump(sim.ti+1)
+                    self.efficacy_dist.jump(self.sim.ti+1)
                     eff_treat_inds = self.efficacy_dist.filter(these_uids)
 
                     post_tx_state_name = thisdf.post_state.values[0]
diff --git a/starsim/sim.py b/starsim/sim.py
index c3f628db..b3026a90 100644
--- a/starsim/sim.py
+++ b/starsim/sim.py
@@ -2,7 +2,6 @@
 Define core Sim classes
 """
 
-# Imports
 import itertools
 import numpy as np
 import sciris as sc
@@ -13,27 +12,36 @@
 
 
 class Sim:
-
+    """
+    The Sim object
+    
+    All Starsim simulations run via the Sim class. It is responsible for initializing
+    and running all modules and generating results.
+    """
     def __init__(self, pars=None, label=None, people=None, demographics=None, diseases=None, networks=None,
                  interventions=None, analyzers=None, connectors=None, copy_inputs=True, **kwargs):
-
         # Make default parameters (using values from parameters.py)
         self.pars = ss.make_pars() # Start with default pars
         args = dict(label=label, people=people, demographics=demographics, diseases=diseases, networks=networks, 
                     interventions=interventions, analyzers=analyzers, connectors=connectors)
         args = {key:val for key,val in args.items() if val is not None} # Remove None inputs
-        self.pars.update(sc.mergedicts(pars, args, kwargs, _copy=copy_inputs))  # Update the parameters
+        input_pars = sc.mergedicts(pars, args, kwargs, _copy=copy_inputs) # TODO: check if copying here is OK
+        self.pars.update(input_pars)  # Update the parameters
         
-        # Set attributes
+        # Set attributes; see also sim.initialize() for more
         self.label = label # Usually overwritten during initialization by the parameters
         self.created = sc.now()  # The datetime the sim was created
+        self.version = ss.__version__ # The Starsim version
+        self.gitinfo = sc.gitinfo(path=__file__, verbose=False)
+        self.dists = ss.Dists(obj=self) # Initialize the random number generator container
         self.initialized = False  # Whether initialization is complete
-        self.complete = False  # Whether a simulation has completed running # TODO: replace with finalized?
+        self.complete = False  # Whether a simulation has completed running
         self.results_ready = False  # Whether results are ready
-        self.dists = ss.Dists(obj=self) # Initialize the random number generator container
         self.results = ss.Results(module='sim')  # For storing results
+        self.elapsed = None # The time required to run
         self.summary = None  # For storing a summary of the results
         self.filename = None # Store the filename, if saved
+
         return
     
     def __getitem__(self, key):
@@ -101,11 +109,11 @@ def initialize(self, **kwargs):
         # Initialize all the modules with the sim
         for mod in self.modules:
             mod.init_pre(self)
-
-        # Initialize products # TODO: think about simplifying
-        for mod in self.interventions:
-            if hasattr(mod, 'product') and isinstance(mod.product, ss.Product):
-                mod.product.init_pre(self)
+                
+        # Initialize products # TODO: think about moving with other modules
+        for intv in self.interventions():
+            if hasattr(intv, 'product'): # TODO: simplify
+                intv.product.init_pre(self)
         
         # Initialize all distributions now that everything else is in place, then set states
         self.dists.initialize(obj=self, base_seed=self.pars.rand_seed, force=True)
@@ -180,87 +188,86 @@ def init_results(self):
     @property
     def modules(self):
         """ Return iterator over all Module instances (stored in standard places) in the Sim """
-        products = [intv.product for intv in self.interventions.values() if
-                    hasattr(intv, 'product') and isinstance(intv.product, ss.Product)]
         return itertools.chain(
-            self.demographics.values(),
-            self.networks.values(),
-            self.diseases.values(),
-            self.connectors.values(),
-            self.interventions.values(),
-            products,
-            self.analyzers.values(),
+            self.demographics(),
+            self.networks(),
+            self.diseases(),
+            self.connectors(),
+            self.interventions(),
+            [intv.product for intv in self.interventions() if hasattr(intv, 'product') and intv.product is not None], # TODO: simplify
+            self.analyzers(),
         )
     
     @property
-    def year(self):
+    def year(self): # TODO: remove when we do the time refactor
         return self.yearvec[self.ti]
 
     def step(self):
-        """ Step through time and update values """
+        """
+        Step through time and update values
+        
+        TODO: assemble every method called during integration in a list, and then call them sequentially        
+        """
 
         # Set the time and if we have reached the end of the simulation, then do nothing
         if self.complete:
-            raise AlreadyRunError('Simulation already complete (call sim.initialize() to re-run)')
+            errormsg = 'Simulation already complete (call sim.initialize() to re-run)'
+            raise AlreadyRunError(errormsg)
 
         # Advance random number generators forward to prepare for any random number calls that may be necessary on this step
         self.dists.jump(to=self.ti+1)  # +1 offset because ti=0 is used on initialization
 
         # Update demographic modules (create new agents from births/immigration, schedule non-disease deaths and emigration)
-        for dem_mod in self.demographics():
-            dem_mod.update()
-
+        for dem in self.demographics():
+            dem.step()
+            
         # Carry out autonomous state changes in the disease modules. This allows autonomous state changes/initializations
         # to be applied to newly created agents
         for disease in self.diseases():
-            disease.update_pre()
+            if isinstance(disease, ss.Disease): # Could be a connector instead -- TODO, rethink this
+                disease.step_state()
 
         # Update connectors -- TBC where this appears in the ordering
         for connector in self.connectors():
-            connector.update()
+            connector.step()
 
-        # Update networks - this takes place here in case autonomous state changes at this timestep affect eligibility for contacts
+        # Update networks - this takes place here in case autonomous state changes at this timestep
+        # affect eligibility for contacts
         for network in self.networks():
-            network.update()
+            network.step()
 
         # Apply interventions - new changes to contacts will be visible and so the final networks can be customized by
         # interventions, by running them at this point
-        for intervention in self.interventions():
-            intervention(self)
-
-        # Carry out transmission/new cases
+        for intv in self.interventions():
+            intv.step()
+        
+        # Carry out autonomous state changes in the disease modules, including transmission (but excluding deaths)
         for disease in self.diseases():
-            disease.make_new_cases()
+            disease.step()
 
         # Execute deaths that took place this timestep (i.e., changing the `alive` state of the agents). This is executed
         # before analyzers have run so that analyzers are able to inspect and record outcomes for agents that died this timestep
-        uids = self.people.resolve_deaths()
+        uids = self.people.check_deaths()
         for disease in self.diseases():
-            disease.update_death(uids)
+            if isinstance(disease, ss.Disease):
+                disease.step_die(uids)
 
         # Update results
         self.people.update_results()
+        for mod in self.modules:
+            mod.update_results()
 
-        for dem_mod in self.demographics():
-            dem_mod.update_results()
-
-        for disease in self.diseases():
-            disease.update_results()
-
-        for analyzer in self.analyzers():
-            analyzer(self)
+        # Apply analyzers
+        for ana in self.analyzers():
+            ana.step()
             
-        # Clean up dead agents
-        self.people.remove_dead()
+        # Clean up dead agents and perform other housekeeping tasks
+        self.people.finish_step()
 
         # Tidy up
         self.ti += 1
-        self.people.ti = self.ti
-        self.people.update_post()
-
         if self.ti == self.npts:
             self.complete = True
-
         return
 
     def run(self, until=None, verbose=None):
@@ -270,14 +277,11 @@ def run(self, until=None, verbose=None):
         T = sc.timer()
         if not self.initialized:
             self.initialize()
-            self._orig_pars = sc.dcp(self.pars)  # Create a copy of the parameters to restore after the run
-
-        if verbose is None:
-            verbose = self.pars.verbose
+        verbose = sc.ifelse(verbose, self.pars.verbose)
 
         # Check for AlreadyRun errors
         errormsg = None
-        if until is None: until = self.npts
+        until = sc.ifelse(until, self.npts)
         if until > self.npts:
             errormsg = f'Requested to run until t={until} but the simulation end is ti={self.npts}'
         if self.ti >= until:  # NB. At the start, self.t is None so this check must occur after initialization
@@ -289,12 +293,8 @@ def run(self, until=None, verbose=None):
 
         # Main simulation loop
         while self.ti < until:
-
-            # Check if we were asked to stop
             elapsed = T.toc(output=True)
-
-            # Print progress
-            if verbose:
+            if verbose: # Print progress
                 simlabel = f'"{self.label}": ' if self.label else ''
                 string = f'  Running {simlabel}{self.yearvec[self.ti]:0.1f} ({self.ti:2.0f}/{self.npts}) ({elapsed:0.2f} s) '
                 if verbose >= 2:
@@ -303,19 +303,19 @@ def run(self, until=None, verbose=None):
                     if not (self.ti % int(1.0 / verbose)):
                         sc.progressbar(self.ti + 1, self.npts, label=string, length=20, newline=True)
 
-            # Actually run the model
+            # Actually run the model -- lots happens on the next line!
             self.step()
 
         # If simulation reached the end, finalize the results
         if self.complete:
-            self.finalize(verbose=verbose)
+            self.ti -= 1  # During the run, this keeps track of the next step; restore this be the final day of the sim
+            self.elapsed = elapsed
+            self.finalize()
             sc.printv(f'Run finished after {elapsed:0.2f} s.\n', 1, verbose)
+        return self # Allows e.g. ss.Sim().run().plot()
 
-        return self
-
-    def finalize(self, verbose=None):
+    def finalize(self):
         """ Compute final results """
-
         if self.results_ready:
             # Because the results are rescaled in-place, finalizing the sim cannot be run more than once or
             # otherwise the scale factor will be applied multiple times
@@ -323,27 +323,26 @@ def finalize(self, verbose=None):
 
         # Scale the results
         for reskey, res in self.results.items():
-            if isinstance(res, ss.Result) and res.scale:
+            if isinstance(res, ss.Result) and res.scale: # NB: disease-specific results are scaled in module.finalize() below
                 self.results[reskey] = self.results[reskey] * self.pars.pop_scale
+        self.results_ready = True # Results are ready to use
 
+        # Finalize each module
         for module in self.modules:
             module.finalize()
 
-        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
+        # Generate the summary and finish up
+        self.summarize() # Create summary
         return
 
     def summarize(self, how='default'):
         """
-        Provide a quick summary of the sim
+        Provide a quick summary of the sim; returns the last entry for count and 
+        cumulative results, and the mean otherwise.
         
         Args:
             how (str): how to summarize: can be 'mean', 'median', 'last', or a mapping of result keys to those
-        
-        Returns the last entry for count and cumulative results, and the mean otherwise
         """
-        
         def get_func(key, how, default='mean'):
             """
             Find the right function by matching the "how" key with the result key
@@ -404,9 +403,9 @@ def shrink(self, skip_attrs=None, in_place=True):
         Returns:
             shrunken (Sim): a Sim object with the listed attributes removed
         """
-        # By default, skip people (~90% of memory), popdict, and _orig_pars (which is just a backup)
+        # By default, skip people
         if skip_attrs is None:
-            skip_attrs = ['people']
+            skip_attrs = ['people'] # TODO: think about skipping all states in all modules
 
         # Create the new object, and copy original dict, skipping the skipped attributes
         if in_place:
@@ -423,127 +422,8 @@ def shrink(self, skip_attrs=None, in_place=True):
         else:
             return shrunken
 
-    def _get_ia(self, which, label=None, partial=False, as_list=False, as_inds=False, die=True, first=False):
-        """ Helper method for get_interventions() and get_analyzers(); see get_interventions() docstring """
-
-        # Handle inputs
-        if which not in ['interventions', 'analyzers']:  # pragma: no cover
-            errormsg = f'This method is only defined for interventions and analyzers, not "{which}"'
-            raise ValueError(errormsg)
-
-        ia_ndict = self.analyzers if which == 'analyzers' else self.interventions  # List of interventions or analyzers
-        n_ia = len(ia_ndict)  # Number of interventions/analyzers
-
-        position = 0 if first else -1  # Choose either the first or last element
-        if label is None:  # Get all interventions if no label is supplied, e.g. sim.get_interventions()
-            label = np.arange(n_ia)
-        if isinstance(label, np.ndarray):  # Allow arrays to be provided
-            label = label.tolist()
-        labels = sc.promotetolist(label)
-
-        # Calculate the matches
-        matches = []
-        match_inds = []
-
-        for label in labels:
-            if sc.isnumber(label):
-                matches.append(ia_ndict[label])
-                label = n_ia + label if label < 0 else label  # Convert to a positive number
-                match_inds.append(label)
-            elif sc.isstring(label) or isinstance(label, type):
-                for ind, ia_key, ia_obj in ia_ndict.enumitems():
-                    if sc.isstring(label) and ia_obj.label == label or (partial and (label in str(ia_obj.label))):
-                        matches.append(ia_obj)
-                        match_inds.append(ind)
-                    elif isinstance(label, type) and isinstance(ia_obj, label):
-                        matches.append(ia_obj)
-                        match_inds.append(ind)
-            else:  # pragma: no cover
-                errormsg = f'Could not interpret label type "{type(label)}": should be str, int, list, or {which} class'
-                raise TypeError(errormsg)
-
-        # Parse the output options
-        if as_inds:
-            output = match_inds
-        elif as_list:  # Used by get_interventions()
-            output = matches
-        else:
-            if len(matches) == 0:  # pragma: no cover
-                if die:
-                    errormsg = f'No {which} matching "{label}" were found'
-                    raise ValueError(errormsg)
-                else:
-                    output = None
-            else:
-                output = matches[
-                    position]  # Return either the first or last match (usually), used by get_intervention()
-
-        return output
-
-    def get_interventions(self, label=None, partial=False, as_inds=False):
-        """
-        Find the matching intervention(s) by label, index, or type. If None, return
-        all interventions. If the label provided is "summary", then print a summary
-        of the interventions (index, label, type).
-
-        Args:
-            label (str, int, Intervention, list): the label, index, or type of intervention to get; if a list, iterate over one of those types
-            partial (bool): if true, return partial matches (e.g. 'beta' will match all beta interventions)
-            as_inds (bool): if true, return matching indices instead of the actual interventions
-        """
-        return self._get_ia('interventions', label=label, partial=partial, as_inds=as_inds, as_list=True)
-
-    def get_intervention(self, label=None, partial=False, first=False, die=True):
-        """
-        Find the matching intervention(s) by label, index, or type.
-        If more than one intervention matches, return the last by default.
-        If no label is provided, return the last intervention in the list.
-
-        Args:
-            label (str, int, Intervention, list): the label, index, or type of intervention to get; if a list, iterate over one of those types
-            partial (bool): if true, return partial matches
-            first (bool): if true, return first matching intervention (otherwise, return last)
-            die (bool): whether to raise an exception if no intervention is found
-        """
-        return self._get_ia('interventions', label=label, partial=partial, first=first, die=die, as_inds=False,
-                            as_list=False)
-
-    def get_analyzers(self, label=None, partial=False, as_inds=False):
-        """
-        Find the matching analyzer(s) by label, index, or type. If None, return
-        all analyzers. If the label provided is "summary", then print a summary
-        of the analyzers (index, label, type).
-
-        Args:
-            label (str, int, Analyzer, list): the label, index, or type of analyzer to get; if a list, iterate over one of those types
-            partial (bool): if true, return partial matches (e.g. 'beta' will match all beta analyzers)
-            as_inds (bool): if true, return matching indices instead of the actual analyzers
-        """
-        return self._get_ia('analyzers', label=label, partial=partial, as_inds=as_inds, as_list=True)
-
-    def get_analyzer(self, label=None, partial=False, first=False, die=True):
-        """
-        Find the matching analyzer(s) by label, index, or type.
-        If more than one analyzer matches, return the last by default.
-        If no label is provided, return the last analyzer in the list.
-
-        Args:
-            label (str, int, Analyzer, list): the label, index, or type of analyzer to get; if a list, iterate over one of those types
-            partial (bool): if true, return partial matches
-            first (bool): if true, return first matching analyzer (otherwise, return last)
-            die (bool): whether to raise an exception if no analyzer is found
-        """
-        return self._get_ia('analyzers', label=label, partial=partial, first=first, die=die, as_inds=False,
-                            as_list=False)
-
     def export_df(self):
-        """
-        Export results as a Pandas dataframe
-
-        :return:
-
-        """
-
+        """ Export results as a Pandas dataframe """
         if not self.results_ready:  # pragma: no cover
             errormsg = 'Please run the sim before exporting the results'
             raise RuntimeError(errormsg)
@@ -559,7 +439,6 @@ def flatten_results(d, prefix=''):
 
         resdict = flatten_results(self.results)
         resdict['t'] = self.yearvec
-
         df = sc.dataframe.from_dict(resdict).set_index('t')
         return df
 
@@ -600,13 +479,11 @@ def save(self, filename=None, keep_people=None, skip_attrs=None, **kwargs):
         else:
             obj = self
         sc.save(filename=filename, obj=obj)
-
         return filename
 
     @staticmethod
     def load(filename, *args, **kwargs):
-        """ Load from disk from a gzipped pickle.  """
-
+        """ Load from disk from a gzipped pickle """
         sim = sc.load(filename, *args, **kwargs)
         if not isinstance(sim, Sim):  # pragma: no cover
             errormsg = f'Cannot load object of {type(sim)} as a Sim object'
@@ -676,27 +553,27 @@ def to_json(self, filename=None, keys=None, tostring=False, indent=2, verbose=Fa
         keys = sc.promotetolist(keys)
 
         # Convert to JSON-compatible format
-        d = {}
+        d = sc.objdict()
         for key in keys:
             if key == 'results':
                 if self.results_ready:
-                    resdict = self.export_results(for_json=True)
-                    d['results'] = resdict
+                    resdict = self.export_results(for_json=True) # TODO: not implemented
+                    d.results = resdict
                 else:
-                    d['results'] = 'Results not available (Sim has not yet been run)'
+                    d.results = 'Results not available (Sim has not yet been run)'
             elif key in ['pars', 'parameters']:
                 pardict = self.export_pars()
-                d['parameters'] = pardict
+                d.parameters = pardict
             elif key == 'summary':
                 if self.results_ready:
-                    d['summary'] = dict(sc.dcp(self.summary))
+                    d.summary = dict(sc.dcp(self.summary))
                 else:
-                    d['summary'] = 'Summary not available (Sim has not yet been run)'
+                    d.summary = 'Summary not available (Sim has not yet been run)'
             elif key == 'short_summary':
                 if self.results_ready:
-                    d['short_summary'] = dict(sc.dcp(self.short_summary))
+                    d.short_summary = dict(sc.dcp(self.short_summary))
                 else:
-                    d['short_summary'] = 'Full summary not available (Sim has not yet been run)'
+                    d.short_summary = 'Full summary not available (Sim has not yet been run)'
             else:  # pragma: no cover
                 try:
                     d[key] = sc.sanitizejson(getattr(self, key))
@@ -772,12 +649,7 @@ def plot(self, key=None, fig=None, style='fancy', fig_kw=None, plot_kw=None):
 
 
 class AlreadyRunError(RuntimeError):
-    """
-    This error is raised if a simulation is run in such a way that no timesteps
-    will be taken. This error is a distinct type so that it can be safely caught
-    and ignored if required, but it is anticipated that most of the time, calling
-    :py:func:`Sim.run` and not taking any timesteps, would be an inadvertent error.
-    """
+    """ Raised if trying to re-run an already-run sim without re-initializing """
     pass
 
 
@@ -888,33 +760,25 @@ def diff_sims(sim1, sim2, skip_key_diffs=False, skip=None, full=False, output=Fa
                 abs_ratio = max(this_ratio, 1.0 / this_ratio)
 
                 # Set the character to use
-                if abs_ratio < small_change:
-                    change_char = '≈'
-                elif new > old:
-                    change_char = '↑'
-                elif new < old:
-                    change_char = '↓'
-                elif new == old:
-                    change_char = '='
+                approx_eq = abs_ratio < small_change
+                if approx_eq:    change_char = '≈'
+                elif new > old:  change_char = '↑'
+                elif new < old:  change_char = '↓'
+                elif new == old: change_char = '='
                 else:
                     errormsg = f'Could not determine relationship between sim1={old} and sim2={new}'
                     raise ValueError(errormsg)
 
                 # Set how many repeats it should have
                 repeats = 1
-                if abs_ratio == 0:
-                    repeats = 0
-                if abs_ratio >= 1.1:
-                    repeats = 2
-                if abs_ratio >= 2:
-                    repeats = 3
-                if abs_ratio >= 10:
-                    repeats = 4
-
+                if abs_ratio == 0:   repeats = 0
+                if abs_ratio >= 1.1: repeats = 2
+                if abs_ratio >= 2:   repeats = 3
+                if abs_ratio >= 10:  repeats = 4
                 this_change = change_char * repeats
             else:  # pragma: no cover
-                this_diff = np.nan
-                this_ratio = np.nan
+                this_diff =  np.nan
+                this_ratio  = np.nan
                 this_change = 'N/A'
 
             diff.append(this_diff)
@@ -945,11 +809,27 @@ def diff_sims(sim1, sim2, skip_key_diffs=False, skip=None, full=False, output=Fa
 
 
 def check_sims_match(*args, full=False):
-    """ Shortcut to using ss.diff_sims() to check if multiple sims match """
-    s1 = args[0]
+    """
+    Shortcut to using ss.diff_sims() to check if multiple sims match
+    
+    Args:
+        args (list): a list of 2 or more sims to compare
+        full (bool): if True, return whether each sim matches the first 
+    
+    **Example**::
+        
+        s1 = ss.Sim(diseases='sir', networks='random')
+        s2 = ss.Sim(pars=dict(diseases='sir', networks='random'))
+        s3 = ss.Sim(diseases=ss.SIR(), networks=ss.RandomNet())
+        assert ss.check_sims_match(s1, s2, s3)
+    """
+    if len(args) < 2:
+        errormsg = 'Must compare at least 2 sims'
+        raise ValueError(errormsg)
+    base = args[0]
     matches = []
-    for s2 in args[1:]:
-        diff = diff_sims(s1, s2, full=False, output=False, die=False)
+    for other in args[1:]:
+        diff = diff_sims(base, other, full=False, output=False, die=False)
         matches.append(not(diff)) # Return the opposite of the diff
     if full:
         return matches
diff --git a/starsim/states.py b/starsim/states.py
index c9c45d17..5a72eee5 100644
--- a/starsim/states.py
+++ b/starsim/states.py
@@ -1,5 +1,15 @@
 """
 Define array-handling classes, including agent states
+
+an event:
+- handles setting states true/false
+- gets uids automatically
+- generates ti_arrays
+- generates results
+
+get working without events/transitions
+define events as scheduling things, then transitions as making them happen
+consider function decorators
 """
 
 import numpy as np
@@ -10,7 +20,7 @@
 ss_int   = ss.dtypes.int
 ss_bool  = ss.dtypes.bool
 
-__all__ = ['check_dtype', 'Arr', 'FloatArr', 'BoolArr', 'IndexArr', 'uids']
+__all__ = ['check_dtype', 'Arr', 'FloatArr', 'BoolArr', 'State', 'IndexArr', 'uids']
 
 
 def check_dtype(dtype, default=None):
@@ -352,6 +362,11 @@ def split(self):
         f_uids = self.false()
         return t_uids, f_uids
 
+
+class State(BoolArr):
+    """ A boolean array being used as a state """
+    pass
+
     
 class IndexArr(Arr):
     """ A special class of Arr used for UIDs and RNG IDs; not to be used as an integer array (for that, use FloatArr) """
@@ -396,15 +411,17 @@ def __new__(cls, arr=None):
         elif isinstance(arr, int): # Convert e.g. ss.uids(0) to ss.uids([0])
             arr = [arr]
         return np.asarray(arr, dtype=ss_int).view(cls) # Handle everything else
-
-    def concat(self, other, **kw): # TODO: why can't they both be called cat()?
-        """ Equivalent to np.concatenate(), but return correct type; see ss.uids.cat() for the class method """
+    
+    def concat(self, other, **kw): # Class and instance methods can't share a name
+        """ Equivalent to np.concatenate(), but return correct type """
         return np.concatenate([self, other], **kw).view(self.__class__)
-
+    
     @classmethod
     def cat(cls, *args, **kw):
-        """ Equivalent to np.concatenate(), but return correct type; see ss.uids.concat() for the instance method """
-        arrs = args[0] if len(args) == 1 else args
+        """ Equivalent to np.concatenate(), but return correct type """
+        if len(args) == 0 or (len(args) == 1 and (args[0] is None or not len(args[0]))):
+            return uids()
+        arrs = args[0] if len(args) == 1 else args # TODO: handle one-array case
         return np.concatenate(arrs, **kw).view(cls)
 
     def remove(self, other, **kw):
@@ -456,4 +473,4 @@ class BooleanOperationError(NotImplementedError):
     """ Raised when a logical operation is performed on a non-logical array """
     def __init__(self, arr):
         msg = f'Logical operations are only valid on Boolean arrays, not {arr.dtype}'
-        super().__init__(msg)
\ No newline at end of file
+        super().__init__(msg)
diff --git a/starsim/utils.py b/starsim/utils.py
index bb7aa20c..7dec2c29 100644
--- a/starsim/utils.py
+++ b/starsim/utils.py
@@ -4,15 +4,15 @@
 
 import warnings
 import numpy as np
-import sciris as sc
-import starsim as ss
 import numba as nb
 import pandas as pd
+import sciris as sc
+import starsim as ss
 
 # %% Helper functions
 
-# What functions are externally visible -- note, this gets populated in each section below
-__all__ = ['ndict', 'warn', 'unique', 'find_contacts', 'combine_rands']
+# What functions are externally visible
+__all__ = ['ndict', 'warn', 'find_contacts', 'set_seed', 'check_requires', 'standardize_netkey', 'standardize_data']
 
 
 class ndict(sc.objdict):
@@ -92,7 +92,7 @@ def _check_key(self, key, overwrite=None):
     def _check_type(self, arg):
         """ Check types """
         if self._type is not None:
-            if not isinstance(arg, self._type):
+            if not isinstance(arg, self._type) and not isinstance(arg, ss.Module): # Module is a valid argument anywhere
                 errormsg = f'The following item does not have the expected type {self._type}:\n{arg}'
                 raise TypeError(errormsg)
         return
@@ -159,18 +159,6 @@ def warn(msg, category=None, verbose=None, die=None):
     return
 
 
-def unique(arr):
-    """
-    Find the unique elements and counts in an array.
-    Equivalent to np.unique(return_counts=True) but ~5x faster, and
-    only works for arrays of positive integers.
-    """
-    counts = np.bincount(arr.ravel())
-    unique = np.flatnonzero(counts)
-    counts = counts[unique]
-    return unique, counts
-
-@nb.njit
 def find_contacts(p1, p2, inds):  # pragma: no cover
     """
     Variation on Network.find_contacts() that avoids sorting.
@@ -189,9 +177,19 @@ def find_contacts(p1, p2, inds):  # pragma: no cover
     return pairing_partners
 
 
-# %% Seed methods
+def check_requires(sim, requires, *args):
+    """ Check that the module's requirements (of other modules) are met """
+    errs = sc.autolist()
+    all_classes = [m.__class__ for m in sim.modules]
+    all_names = [m.name for m in sim.modules]
+    for req in sc.mergelists(requires, *args):
+        if req not in all_classes + all_names:
+            errs += req
+    if len(errs):
+        errormsg = f'The following module(s) are required, but the Sim does not contain them: {sc.strjoin(errs)}'
+        raise AttributeError(errormsg)
+    return
 
-__all__ += ['set_seed']
 
 def set_seed(seed=None):
     '''
@@ -224,8 +222,10 @@ def set_seed_regular(seed):
 
 # %% Data cleaning and processing
 
-__all__ += ['standardize_data']
 
+def standardize_netkey(key):
+    """ Networks can be upper or lowercase, and have a suffix 'net' or not; this function standardizes them """
+    return key.lower().removesuffix('net')
 
 def standardize_data(data=None, metadata=None, min_year=1800, out_of_range=0, default_age=0, default_year=2024):
     """
@@ -319,10 +319,13 @@ def standardize_data(data=None, metadata=None, min_year=1800, out_of_range=0, de
 
 def combine_rands(a, b):
     """
-    Efficient algorithm for combining two arrays of random numbers into one
+    Efficient algorithm for combining two arrays of random integers into an array
+    of floats.
+    
+    See ss.multi_random() for the user-facing version.
     
     Args:
-        a (array): array of random integers between np.iinfo(np.int64).min and np.iinfo(np.int64).max
+        a (array): array of random integers between 0 and np.iinfo(np.uint64).max, as from ss.rand_raw()
         b (array): ditto, same size as a
         
     Returns:
diff --git a/starsim/version.py b/starsim/version.py
index 8a6d0348..7f3908e0 100644
--- a/starsim/version.py
+++ b/starsim/version.py
@@ -4,6 +4,6 @@
 
 __all__ = ['__version__', '__versiondate__', '__license__']
 
-__version__ = '1.0.1'
-__versiondate__ = '2024-07-22'
+__version__ = '1.1.0'
+__versiondate__ = '2024-09-10'
 __license__ = f'Starsim {__version__} ({__versiondate__}) — © 2023-2024 by IDM'


Test patch:
diff --git a/tests/baseline.json b/tests/baseline.json
index 52b867bb..8005c3d2 100644
--- a/tests/baseline.json
+++ b/tests/baseline.json
@@ -1,26 +1,26 @@
 {
   "summary": {
     "yearvec": 2010.1000000000026,
-    "births_new": 65.20588235294117,
-    "births_cumulative": 3261.0,
-    "births_cbr": 29.716475754608183,
-    "deaths_new": 44.72549019607843,
-    "deaths_cumulative": 2227.0686274509803,
-    "deaths_cmr": 20.375999403074815,
-    "sir_n_susceptible": 2397.764705882353,
-    "sir_n_infected": 3586.794117647059,
-    "sir_n_recovered": 4986.558823529412,
-    "sir_prevalence": 0.3357567320197343,
-    "sir_new_infections": 128.8235294117647,
-    "sir_cum_infections": 13140.0,
-    "sis_n_susceptible": 4048.8333333333335,
-    "sis_n_infected": 6922.28431372549,
-    "sis_prevalence": 0.6181125674334483,
-    "sis_new_infections": 207.7941176470588,
-    "sis_cum_infections": 21195.0,
-    "sis_rel_sus": 0.4713089116600518,
-    "n_alive": 10971.117647058823,
-    "new_deaths": 45.96078431372549,
-    "cum_deaths": 4650.0
+    "births_new": 65.24509803921569,
+    "births_cumulative": 3261.872549019608,
+    "births_cbr": 29.0,
+    "deaths_new": 44.6764705882353,
+    "deaths_cumulative": 2225.4607843137255,
+    "deaths_cmr": 19.794117647058822,
+    "sir_n_susceptible": 2416.3627450980393,
+    "sir_n_infected": 3588.0686274509803,
+    "sir_n_recovered": 4975.14705882353,
+    "sir_prevalence": 0.3356119858776437,
+    "sir_new_infections": 128.97058823529412,
+    "sir_cum_infections": 13155.0,
+    "sis_n_susceptible": 4022.2156862745096,
+    "sis_n_infected": 6957.362745098039,
+    "sis_prevalence": 0.6209749044525374,
+    "sis_new_infections": 208.13725490196077,
+    "sis_cum_infections": 21230.0,
+    "sis_rel_sus": 0.4676589142075258,
+    "n_alive": 10979.578431372549,
+    "new_deaths": 45.76470588235294,
+    "cum_deaths": 4630.0
   }
 }
\ No newline at end of file
diff --git a/tests/benchmark.json b/tests/benchmark.json
index 9ec32896..4b148c87 100644
--- a/tests/benchmark.json
+++ b/tests/benchmark.json
@@ -1,12 +1,12 @@
 {
   "time": {
-    "initialize": 0.015,
-    "run": 0.981
+    "initialize": 0.014,
+    "run": 0.907
   },
   "parameters": {
     "n_agents": 10000,
     "n_years": 20,
     "dt": 0.2
   },
-  "cpu_performance": 0.8852149360177684
+  "cpu_performance": 0.8939678660976156
 }
\ No newline at end of file
diff --git a/tests/test_demographics.py b/tests/test_demographics.py
index 2e14c837..f297cbfe 100644
--- a/tests/test_demographics.py
+++ b/tests/test_demographics.py
@@ -88,8 +88,9 @@ def test_nigeria(which='births', dt=1, start=1995, n_years=15, do_plot=False):
             assert np.array_equal(sim.results.pregnancy.pregnancies, sim.results.pregnancy.births)
             print('✓ (births == pregnancies)')
 
-    print("Check final pop size within 5% of data")
-    assert np.isclose(data.n_alive.values[-1], sim.results.n_alive[-1], rtol=0.05)
+    rtol = 0.05
+    print(f'Check final pop size within {rtol*100:n}% of data')
+    assert np.isclose(data.n_alive.values[-1], sim.results.n_alive[-1], rtol=rtol), f'Final population size not within {rtol*100:n}% of data'
     print(f'✓ (simulated/data={sim.results.n_alive[-1] / data.n_alive.values[-1]:.2f})')
 
     # Plots
diff --git a/tests/test_diseases.py b/tests/test_diseases.py
index 4de6f93e..37542a2b 100644
--- a/tests/test_diseases.py
+++ b/tests/test_diseases.py
@@ -123,7 +123,7 @@ def test_ncd():
     sim.run()
 
     assert len(ncd.log.out_edges) == ncd.log.number_of_edges()
-    df = ncd.log.line_list  # Check generation of line-list
+    df = ncd.log.line_list()  # Check generation of line-list
     assert df.source.isna().all()
 
     plt.figure()
diff --git a/tests/test_plugins.py b/tests/test_plugins.py
index 0bc9030f..2729e50b 100644
--- a/tests/test_plugins.py
+++ b/tests/test_plugins.py
@@ -13,8 +13,9 @@
 class hiv_syph(ss.Connector):
     """ Simple connector whereby rel_sus to NG doubles if CD4 count is <200"""
     def __init__(self, pars=None, **kwargs):
-        super().__init__(label='HIV-Syphilis', requires=[ss.HIV, ss.Syphilis])
-        self.default_pars(
+        super().__init__()
+        self.define_pars(
+            label = 'HIV-Syphilis',
             rel_sus_syph_hiv    = 2,   # People with HIV are 2x more likely to acquire syphilis
             rel_sus_syph_aids   = 5,   # People with AIDS are 5x more likely to acquire syphilis
             rel_trans_syph_hiv  = 1.5, # People with HIV are 1.5x more likely to transmit syphilis
@@ -25,7 +26,7 @@ def __init__(self, pars=None, **kwargs):
         self.update_pars(pars, **kwargs)
         return
 
-    def update(self):
+    def step(self):
         """ Specify HIV-syphilis interactions """
         diseases = self.sim.diseases
         syph = diseases.syphilis
@@ -56,7 +57,8 @@ def __init__(self, year=2020, prob=0.8):
         self.year = year
         return
 
-    def apply(self, sim):
+    def step(self):
+        sim = self.sim
         if sim.year > self.year:
             syphilis = sim.diseases.syphilis
 
@@ -154,4 +156,4 @@ def test_connectors(do_plot=False):
     sims = test_connectors(do_plot=do_plot)
     
     T.toc()
-    
\ No newline at end of file
+    
diff --git a/tests/test_randomness.py b/tests/test_randomness.py
index f7255aee..1a9452cb 100644
--- a/tests/test_randomness.py
+++ b/tests/test_randomness.py
@@ -112,27 +112,35 @@ def test_order(n=n):
 class CountInf(ss.Intervention):
     """ Store every infection state in a timepoints x people array """
     def init_pre(self, sim):
+        super().init_pre(sim)
         n_agents = len(sim.people)
         self.arr = np.zeros((sim.npts, n_agents))
         self.n_agents = n_agents
         return
     
-    def apply(self, sim):
-        self.arr[sim.ti, :] = np.array(sim.diseases.sir.infected)[:self.n_agents]
+    def step(self):
+        self.arr[self.sim.ti, :] = np.array(self.sim.diseases.sir.infected)[:self.n_agents]
         return
 
 
 class OneMore(ss.Intervention):
     """ Add one additional agent and infection """
+    def __init__(self, ti_apply=10):
+        super().__init__()
+        self.ti_apply = ti_apply
+        return
+    
     def init_pre(self, sim):
+        super().init_pre(sim)
         one_birth = ss.Pregnancy(name='one_birth', rel_fertility=0) # Ensure no default births
         one_birth.init_pre(sim)
         self.one_birth = one_birth
         return
     
-    def apply(self, sim, ti=10):
+    def step(self):
         """ Create an extra agent """
-        if sim.ti == ti:
+        sim = self.sim
+        if sim.ti == self.ti_apply:
             new_uids = self.one_birth.make_embryos(ss.uids(0)) # Assign 0th agent to be the "mother"
             sim.people.age[new_uids] = -100 # Set to a very low number to never reach debut age
             
@@ -293,7 +301,7 @@ def test_combine_rands(do_plot=False):
     np.random.seed(2)
     a = np.random.randint(np.iinfo(np.uint64).max, size=n, dtype=np.uint64)
     b = np.random.randint(np.iinfo(np.uint64).max, size=n, dtype=np.uint64)
-    c = ss.combine_rands(a, b)
+    c = ss.utils.combine_rands(a, b)
     if do_plot:
         pl.figure()
         for i,k,v in sc.objdict(a=a,b=b,combined=c).enumitems():
diff --git a/tests/test_syphilis.py b/tests/test_syphilis.py
index e0e13564..e36b5ed1 100644
--- a/tests/test_syphilis.py
+++ b/tests/test_syphilis.py
@@ -56,13 +56,11 @@ def __init__(self, *args, **kwargs):
         self.okay = True
         return
 
-    def update_results(self, sim):
-        return self.apply(sim)
-
-    def apply(self, sim):
+    def update_results(self):
         """
         Checks states that should be mutually exlusive and collectively exhaustive
         """
+        sim = self.sim
         sppl = sim.diseases.syphilis
 
         # Infection states: people must be exactly one of these
@@ -138,7 +136,7 @@ def test_syph_intvs(dt=1, n_agents=500, do_plot=False):
         label='syph_screening',
     )
 
-    treat_eligible = lambda sim: ss.uids(sim.get_intervention('syph_screening').outcomes['positive'])
+    treat_eligible = lambda sim: ss.uids(sim.interventions['syph_screening'].outcomes['positive'])
     bpg = ss.syph_treatment(
         prob=0.9,
         product='bpg',