Repository: starsimhub/starsim (language: python)

Issue:
Make jump robust
Makes the changes to `jump` discussed on the call -- which already caught two bugs in the code (both of which were due to ancient holdovers from the `MultiRNG` days).

If these changes look good, I'll backport `jump_dt()` and the updated `jump()` to Starsim v1.0.


Discussion / hints:

Fix patch:
diff --git a/CHANGELOG.rst b/CHANGELOG.rst
index eed7b2a9..d98225ba 100644
--- a/CHANGELOG.rst
+++ b/CHANGELOG.rst
@@ -12,21 +12,20 @@ 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()``).
+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), and changes to module types and integration (e.g. renaming ``update()`` and ``apply()`` methods to ``step()``).
 
-Time-aware parameters
-~~~~~~~~~~~~~~~~~~~~~
-Coming soon!
-
-States and events
-~~~~~~~~~~~~~~~~~
-Coming soon!
+Time-aware parameters and modules
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+- Added ``ss.dur()``, ``ss.rate()``, and ``ss.time_prob()`` classes, for automatic handling of time units in simulations. There are also convenience classes ``ss.days()``, ``ss.years()``, and ``ss.beta()`` for special cases of these.
+- Durations and rates, along with modules and the sim itself, now have a ``unit`` parameter which can be ``'day'``, ``'week'``, ``'month'``, or ``'year'`` (default). Modules now also have their own timestep ``dt``. Different units and timesteps can be mixed and matched.
+- There is a new ``Loop`` class which handles the integration loop. You can view the integration plan via ``sim.loop.to_df()`` or ``sim.loop.plot()``.
+- There are more advanced debugging tools. You can run a single sim timestep with ``sim.run_one_step()`` (which in turn calls multiple functions), and you can run a single function from the integration loop with ``sim.loop.run_one_step()``.
 
 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).
+- ``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 ``Network.step()``.
+- 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). For both the sim and modules, ``initialize()`` has been renamed ``init()``.
 - All modules are treated the same in the integration loop, except for diseases, which have ``step_state()`` and ``step_die()`` methods.
 
 Other changes
@@ -38,21 +37,25 @@ Other changes
 - ``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. Call it manually from ``init_pre()`` if desired, e.g. ``ss.check_requires(self.sim, self.requires)``.
 - 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.
+- Distributions have a new ``jump_dt`` method that jumps by much more than a single state update.
+- ``ss.diff_sims()`` can now handle ``MultiSim` objects.
 - ``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.
+- All duration and rate parameters should now be wrapped with ``ss.dur()`` and ``ss.rate()``. Events that represent probabilities over time (i.e. hazard rates) can also be wrapped with ``ss.time_prob()``, although this is similar to ``ss.rate()`` unless the value is relatively large.
 - ``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()``
+- ``initialize()`` has been renamed ``init()``.
 - ``Demographics.update()`` has been renamed ``Demographics.step()``.
 - ``Network.update()`` has been renamed ``Network.step()``.
-- ``Disease.update_pre()`` has been renamed ``Disease.step_pre()``.
+- ``Disease.update_pre()`` has been renamed ``Disease.step_state()``.
 - ``Disease.make_new_cases()`` has been renamed ``Disease.step()``.
-- ``Disease.update_death()`` has been renamed ``Disease.die()``.
+- ``Disease.update_death()`` has been renamed ``Disease.step_die()`` (which is now called by ``People.step_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()``).
diff --git a/docs/tutorials/dev_tut_abm.ipynb b/docs/tutorials/dev_tut_abm.ipynb
index 6f3cb382..4442029a 100644
--- a/docs/tutorials/dev_tut_abm.ipynb
+++ b/docs/tutorials/dev_tut_abm.ipynb
@@ -46,11 +46,11 @@
     "        self.wealth_dist = np.zeros((self.npts, len(self.bins)-1)) # Wealth distribution over time\n",
     "        return\n",
     "        \n",
-    "    def apply(self, sim):\n",
+    "    def step(self):\n",
     "        \"\"\" Transfer wealth between agents -- core model logic \"\"\"\n",
     "        self.wealth_hist() # Store the wealth at this time point\n",
     "        givers = self.wealth > 0 # People need wealth to be givers\n",
-    "        receivers = np.random.choice(sim.people.uid, size=givers.sum()) # Anyone can be a receiver\n",
+    "        receivers = np.random.choice(self.sim.people.uid, size=givers.sum()) # Anyone can be a receiver\n",
     "        self.wealth[givers] -= 1 # Givers are unique, so can use vectorized version\n",
     "        for receive in receivers: # Vectorized version is: np.add.at(sim.people.wealth.raw, receivers, 1)\n",
     "            self.wealth[receive] += 1\n",
@@ -88,8 +88,8 @@
     "pars = dict(\n",
     "    n_agents = 100, # Number of agents\n",
     "    start = 0,\n",
-    "    end = 100,\n",
-    "    interventions = wealth,\n",
+    "    stop = 100,\n",
+    "    demographics = wealth,\n",
     ")\n",
     "\n",
     "# Run the model\n",
diff --git a/docs/tutorials/dev_tut_analyzers.ipynb b/docs/tutorials/dev_tut_analyzers.ipynb
index facfa12c..048cce7b 100644
--- a/docs/tutorials/dev_tut_analyzers.ipynb
+++ b/docs/tutorials/dev_tut_analyzers.ipynb
@@ -77,7 +77,7 @@
     "res = sim.results\n",
     "\n",
     "plt.figure()\n",
-    "plt.plot(res.yearvec, res.hiv_in_pregnancy.new_infections_pregnancy)\n",
+    "plt.plot(res.timevec, res.hiv_in_pregnancy.new_infections_pregnancy)\n",
     "plt.title('HIV infections acquired during pregnancy')\n",
     "plt.show();"
    ]
@@ -99,7 +99,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.10.11"
+   "version": "3.11.7"
   },
   "pycharm": {
    "stem_cell": {
diff --git a/docs/tutorials/tut_buildsim.ipynb b/docs/tutorials/tut_buildsim.ipynb
index 0a374c47..bf6c0a2c 100644
--- a/docs/tutorials/tut_buildsim.ipynb
+++ b/docs/tutorials/tut_buildsim.ipynb
@@ -96,7 +96,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": null,
+   "execution_count": 9,
    "metadata": {
     "collapsed": false
    },
@@ -130,13 +130,13 @@
   },
   {
    "cell_type": "code",
-   "execution_count": null,
+   "execution_count": 10,
    "metadata": {
     "collapsed": false
    },
    "outputs": [],
    "source": [
-    "sim = ss.Sim(start=2020, n_years=1, dt=1/12)  # Simulate from 2020 for 1 year with a monthly timestep"
+    "sim = ss.Sim(start=2020, dur=1, dt=1/12)  # Simulate from 2020 for 1 year with a monthly timestep"
    ]
   },
   {
diff --git a/docs/tutorials/tut_demographics.ipynb b/docs/tutorials/tut_demographics.ipynb
index 93c3981c..6f3fbe21 100644
--- a/docs/tutorials/tut_demographics.ipynb
+++ b/docs/tutorials/tut_demographics.ipynb
@@ -11,7 +11,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": null,
+   "execution_count": 2,
    "metadata": {
     "collapsed": false
    },
@@ -40,7 +40,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": null,
+   "execution_count": 3,
    "metadata": {
     "collapsed": false
    },
@@ -67,7 +67,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": null,
+   "execution_count": 4,
    "metadata": {
     "collapsed": false
    },
@@ -100,6 +100,7 @@
    "source": [
     "import starsim as ss\n",
     "import pandas as pd\n",
+    "import matplotlib.pyplot as plt\n",
     "\n",
     "# Read in age-specific fertility rates\n",
     "fertility_rates = pd.read_csv('test_data/nigeria_asfr.csv')\n",
@@ -125,10 +126,9 @@
     "data = nigeria_popsize[(nigeria_popsize.year >= 1995) & (nigeria_popsize.year <= 2030)]\n",
     "\n",
     "# Plot the overall population size - simulated vs data\n",
-    "import matplotlib.pyplot as plt\n",
     "fig, ax = plt.subplots(1, 1)\n",
     "ax.scatter(data.year, data.n_alive, alpha=0.5, label='Data')\n",
-    "ax.plot(sim.yearvec, sim.results.n_alive, color='k', label='Model')\n",
+    "ax.plot(sim.timevec, sim.results.n_alive, color='k', label='Model')\n",
     "ax.legend()\n",
     "ax.set_title('Population')\n",
     "plt.show();"
diff --git a/docs/tutorials/tut_interventions.ipynb b/docs/tutorials/tut_interventions.ipynb
index 053294c9..aa038198 100644
--- a/docs/tutorials/tut_interventions.ipynb
+++ b/docs/tutorials/tut_interventions.ipynb
@@ -35,7 +35,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": null,
+   "execution_count": 2,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -114,8 +114,8 @@
     "res_intv = sim_intv.results\n",
     "\n",
     "plt.figure()\n",
-    "plt.plot(res_base.yearvec, res_base.sir.prevalence, label='Baseline')\n",
-    "plt.plot(res_intv.yearvec, res_intv.sir.prevalence, label='Vax')\n",
+    "plt.plot(res_base.timevec, res_base.sir.prevalence, label='Baseline')\n",
+    "plt.plot(res_intv.timevec, res_intv.sir.prevalence, label='Vax')\n",
     "plt.axvline(x=2015, color='k', ls='--')\n",
     "plt.title('Prevalence')\n",
     "plt.legend()\n",
diff --git a/docs/tutorials/tut_networks.ipynb b/docs/tutorials/tut_networks.ipynb
index 939134cb..6a75c8eb 100644
--- a/docs/tutorials/tut_networks.ipynb
+++ b/docs/tutorials/tut_networks.ipynb
@@ -22,7 +22,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": null,
+   "execution_count": 7,
    "metadata": {
     "tags": []
    },
@@ -54,7 +54,7 @@
     "network = ss.RandomNet()\n",
     "sir = ss.SIR()\n",
     "sim = ss.Sim(n_agents=2_000, diseases=sir, networks=network)\n",
-    "sim.initialize()"
+    "sim.init()"
    ]
   },
   {
@@ -71,7 +71,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": null,
+   "execution_count": 9,
    "metadata": {
     "collapsed": false
    },
@@ -88,7 +88,7 @@
     "    acts=80,        # Distribution of the number of coital acts/year\n",
     ")\n",
     "\n",
-    "pars = dict(dt=1/12, start=2000, n_years=20)\n",
+    "pars = dict(start=2000, dur=20, dt=1/12)\n",
     "sim = ss.Sim(pars=pars, diseases=syph, networks=mf)"
    ]
   },
@@ -110,7 +110,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": null,
+   "execution_count": 10,
    "metadata": {
     "collapsed": false
    },
@@ -143,7 +143,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": null,
+   "execution_count": 11,
    "metadata": {
     "collapsed": false
    },
diff --git a/setup.py b/setup.py
index 4d4ae2f9..95046d9d 100644
--- a/setup.py
+++ b/setup.py
@@ -41,7 +41,7 @@
         'numpy',
         'scipy',
         'pandas>=2.0.0',
-        'sciris>=3.1.6',
+        'sciris>=3.2.0',
         'matplotlib',
         'numba',
         'networkx',
diff --git a/starsim/__init__.py b/starsim/__init__.py
index 19c32818..f6e19630 100644
--- a/starsim/__init__.py
+++ b/starsim/__init__.py
@@ -1,6 +1,7 @@
 from .version import __version__, __versiondate__, __license__
 from .settings      import *
 from .utils         import *
+from .timepars      import *
 from .parameters    import *
 from .distributions import *
 from .states        import *
@@ -14,6 +15,7 @@
 from .demographics  import *
 from .disease       import *
 from .diseases      import *
+from .loop          import *
 from .sim           import *
 from .run           import *
 from .calibration   import *
@@ -28,7 +30,7 @@
     print(__license__)
 
 # Double-check key requirements -- should match setup.py
-reqs = ['sciris>=3.1.6', 'pandas>=2.0.0', 'scipy', 'numba', 'networkx']
+reqs = ['sciris>=3.2.0', 'pandas>=2.0.0', 'scipy', 'numba', 'networkx']
 msg = f'The following dependencies for Starsim {__version__} were not met: <MISSING>.'
 sc.require(reqs, message=msg)
 del sc, reqs, msg # Don't keep this in the module
\ No newline at end of file
diff --git a/starsim/calibration.py b/starsim/calibration.py
index 677ce47b..84a3d678 100644
--- a/starsim/calibration.py
+++ b/starsim/calibration.py
@@ -177,7 +177,7 @@ def __init__(self, sim, data, calib_pars=None, weights=None, fit_args=None, par_
 
         # Initialize sim
         if not self.sim.initialized:
-            self.sim.initialize()
+            self.sim.init()
 
         # Figure out which sim results to get
         self.sim_result_list = self.target_data.columns.values.tolist()
diff --git a/starsim/demographics.py b/starsim/demographics.py
index 80e9a7da..a9fb51db 100644
--- a/starsim/demographics.py
+++ b/starsim/demographics.py
@@ -56,7 +56,7 @@ def init_pre(self, sim):
         if isinstance(self.pars.birth_rate, pd.DataFrame):
             br_year = self.pars.birth_rate[self.metadata.data_cols['year']]
             br_val = self.pars.birth_rate[self.metadata.data_cols['cbr']]
-            all_birth_rates = np.interp(sim.yearvec, br_year, br_val)
+            all_birth_rates = np.interp(self.timevec, br_year, br_val) # This assumes a year timestep -- probably ok?
             self.pars.birth_rate = all_birth_rates
         return
 
@@ -68,11 +68,10 @@ def standardize_birth_data(self):
         return birth_rate
 
     def init_results(self):
-        npts = self.sim.npts
         self.results += [
-            ss.Result(self.name, 'new', npts, dtype=int, scale=True, label='New births'),
-            ss.Result(self.name, 'cumulative', npts, dtype=int, scale=True, label='Cumulative births'),
-            ss.Result(self.name, 'cbr', npts, dtype=int, scale=False, label='Crude birth rate'),
+            ss.Result(self.name, 'new',        self.npts, dtype=int, scale=True,  label='New births'),
+            ss.Result(self.name, 'cumulative', self.npts, dtype=int, scale=True,  label='Cumulative births'),
+            ss.Result(self.name, 'cbr',        self.npts, dtype=int, scale=False, label='Crude birth rate'),
         ]
         return
 
@@ -85,7 +84,7 @@ def get_births(self):
 
         if isinstance(p.birth_rate, (pd.Series, pd.DataFrame)):
             available_years = p.birth_rate.index
-            year_ind = sc.findnearest(available_years, sim.year)
+            year_ind = sc.findnearest(available_years, sim.now) # TODO: make robust to timestep
             nearest_year = available_years[year_ind]
             this_birth_rate = p.birth_rate.loc[nearest_year]
         else:
@@ -110,14 +109,14 @@ def add_births(self):
         return new_uids
 
     def update_results(self):
-        self.results['new'][self.sim.ti] = self.n_births
+        self.results['new'][self.ti] = self.n_births
         return
 
     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.cbr[:] = 1/self.pars.units*np.divide(self.results.new/self.sim.dt_year, res.n_alive, where=res.n_alive>0)
         return
 
 
@@ -198,7 +197,7 @@ def make_death_prob_fn(self, sim, uids):
             assert len(uids) == len(ppl.auids)
 
             available_years = drd.index.get_level_values('year')
-            year_ind = sc.findnearest(available_years, sim.year)
+            year_ind = sc.findnearest(available_years, sim.now) # TODO: make work with different timesteps
             nearest_year = available_years[year_ind]
 
             death_rate = np.empty(uids.shape, dtype=ss_float_)
@@ -222,11 +221,10 @@ def make_death_prob_fn(self, sim, uids):
         return death_prob
 
     def init_results(self):
-        npts = self.sim.npts
         self.results += [
-            ss.Result(self.name, 'new', npts, dtype=int, scale=True, label='Deaths'),
-            ss.Result(self.name, 'cumulative', npts, dtype=int, scale=True, label='Cumulative deaths'),
-            ss.Result(self.name, 'cmr', npts, dtype=int, scale=False, label='Crude mortality rate'),
+            ss.Result(self.name, 'new',        self.npts, dtype=int, scale=True, label='Deaths'),
+            ss.Result(self.name, 'cumulative', self.npts, dtype=int, scale=True, label='Cumulative deaths'),
+            ss.Result(self.name, 'cmr',        self.npts, dtype=int, scale=False, label='Crude mortality rate'),
         ]
         return
 
@@ -238,14 +236,14 @@ def step(self):
         return self.n_deaths
     
     def update_results(self):
-        self.results['new'][self.sim.ti] = self.n_deaths
+        self.results['new'][self.ti] = self.n_deaths
         return
 
     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.cmr[:] = 1/self.pars.units*np.divide(self.results.new / self.sim.dt_year, n_alive, where=n_alive>0)
         return
 
 
@@ -254,8 +252,8 @@ class Pregnancy(Demographics):
     def __init__(self, pars=None, metadata=None, **kwargs):
         super().__init__()
         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)
+            dur_pregnancy = ss.years(0.75), # Duration for pre-natal transmission
+            dur_postpartum = ss.lognorm_ex(mean=ss.years(0.5), std=ss.years(0.5)), # Duration for post-natal transmission (e.g. via breastfeeding)
             fertility_rate = 0, # Can be a number of Pandas DataFrame
             rel_fertility = 1,
             maternal_death_prob = ss.bernoulli(0),
@@ -286,7 +284,7 @@ def __init__(self, pars=None, metadata=None, **kwargs):
             sc.objdict(data_cols=dict(year='Time', age='AgeGrp', value='ASFR')),
             metadata,
         )
-        self.choose_slots = None # Distribution for choosing slots; set in self.initialize()
+        self.choose_slots = None # Distribution for choosing slots; set in self.init()
 
         # For results tracking
         self.n_pregnancies = 0
@@ -305,7 +303,7 @@ def make_fertility_prob_fn(self, sim, uids):
         if sc.isnumber(frd):
             fertility_rate[uids] = self.fertility_rate_data
         else:
-            year_ind = sc.findnearest(frd.index, sim.year-self.pars.dur_pregnancy)
+            year_ind = sc.findnearest(frd.index, sim.now-self.pars.dur_pregnancy) # TODO: make time-unit-aware
             nearest_year = frd.index[year_ind]
 
             # Assign agents to age bins
@@ -371,11 +369,10 @@ def init_results(self):
         Still unclear whether this logic should live in the pregnancy module, the
         individual disease modules, the connectors, or the sim.
         """
-        npts = self.sim.npts
         self.results += [
-            ss.Result(self.name, 'pregnancies', npts, dtype=int, scale=True, label='New pregnancies'),
-            ss.Result(self.name, 'births', npts, dtype=int, scale=True, label='New births'),
-            ss.Result(self.name, 'cbr', npts, dtype=int, scale=False, label='Crude birth rate'),
+            ss.Result(self.name, 'pregnancies', self.npts, dtype=int, scale=True,  label='New pregnancies'),
+            ss.Result(self.name, 'births',      self.npts, dtype=int, scale=True,  label='New births'),
+            ss.Result(self.name, 'cbr',         self.npts, dtype=int, scale=False, label='Crude birth rate'),
         ]
         return
 
@@ -390,7 +387,7 @@ def step(self):
     def update_states(self):
         """ Update states """
         # Check for new deliveries
-        ti = self.sim.ti
+        ti = self.ti
         deliveries = self.pregnant & (self.ti_delivery <= ti)
         self.n_births = np.count_nonzero(deliveries)
         self.pregnant[deliveries] = False
@@ -406,7 +403,7 @@ def update_states(self):
                 prenatalnet = [nw for nw in self.sim.networks.values() if nw.prenatal][0]
 
                 # Find the prenatal connections that are ending
-                prenatal_ending = prenatalnet.edges.end<=self.sim.ti
+                prenatal_ending = prenatalnet.edges.end <= ti
                 new_mother_uids = prenatalnet.edges.p1[prenatal_ending]
                 new_infant_uids = prenatalnet.edges.p2[prenatal_ending]
 
@@ -417,7 +414,7 @@ def update_states(self):
 
                 # Create durations and start dates, and add connections
                 durs = self.dur_postpartum[new_mother_uids]
-                start = np.full(self.n_births, fill_value=self.sim.ti)
+                start = np.full(self.n_births, fill_value=ti)
 
                 # # Remove pairs from prenatal network and add to postnatal
                 prenatalnet.end_pairs()
@@ -443,7 +440,7 @@ def make_pregnancies(self):
         # Validation
         if np.any(self.pregnant[conceive_uids]):
             which_uids = conceive_uids[self.pregnant[conceive_uids]]
-            errormsg = f'New conceptions registered in {len(which_uids)} pregnant agent(s) at timestep {self.sim.ti}.'
+            errormsg = f'New conceptions registered in {len(which_uids)} pregnant agent(s) at timestep {self.ti}.'
             raise ValueError(errormsg)
 
         # Set prognoses for the pregnancies
@@ -472,7 +469,7 @@ def make_embryos(self, conceive_uids):
             for lkey, layer in self.sim.networks.items():
                 if layer.prenatal:
                     durs = np.full(n_unborn, fill_value=self.pars.dur_pregnancy)
-                    start = np.full(n_unborn, fill_value=self.sim.ti)
+                    start = np.full(n_unborn, fill_value=self.ti)
                     layer.add_pairs(conceive_uids, new_uids, dur=durs, start=start)
 
         return new_uids
@@ -486,18 +483,18 @@ def set_prognoses(self, uids):
         """
 
         # Change states for the newly pregnant woman
-        ti = self.sim.ti
-        dt = self.sim.dt
+        ti = self.ti
+        dt = self.dt # TODO: CHECK
         self.fecund[uids] = False
         self.pregnant[uids] = True
         self.ti_pregnant[uids] = ti
 
         # Outcomes for pregnancies
-        dur_preg = np.full(len(uids), self.pars.dur_pregnancy)  # Duration in years
+        dur_preg = np.ones(len(uids))*self.pars.dur_pregnancy  # Duration in years
         dur_postpartum = self.pars.dur_postpartum.rvs(uids)
         dead = self.pars.maternal_death_prob.rvs(uids)
-        self.ti_delivery[uids] = ti + dur_preg/dt # Currently assumes maternal deaths still result in a live baby
-        self.ti_postpartum[uids] = self.ti_delivery[uids] + dur_postpartum/dt
+        self.ti_delivery[uids] = ti + dur_preg # Currently assumes maternal deaths still result in a live baby
+        self.ti_postpartum[uids] = self.ti_delivery[uids] + dur_postpartum
         self.dur_postpartum[uids] = dur_postpartum
 
         if np.any(dead): # NB: 100x faster than np.sum(), 10x faster than np.count_nonzero()
@@ -505,7 +502,7 @@ def set_prognoses(self, uids):
         return
 
     def update_results(self):
-        ti = self.sim.ti
+        ti = self.ti
         self.results['pregnancies'][ti] = self.n_pregnancies
         self.results['births'][ti] = self.n_births
         return
@@ -513,5 +510,5 @@ def update_results(self):
     def finalize(self):
         super().finalize()
         n_alive = self.sim.results.n_alive
-        self.results['cbr'] = 1/self.pars.units * np.divide(self.results['births'] / self.sim.dt, n_alive, where=n_alive>0)
+        self.results['cbr'] = 1/self.pars.units * np.divide(self.results['births'] / self.sim.dt_year, n_alive, where=n_alive>0)
         return
diff --git a/starsim/disease.py b/starsim/disease.py
index a827d02a..9a593070 100644
--- a/starsim/disease.py
+++ b/starsim/disease.py
@@ -55,7 +55,7 @@ def init_results(self):
         Result for 'n_susceptible'
         """
         for state in self._boolean_states:
-            self.results += ss.Result(self.name, f'n_{state.name}', self.sim.npts, dtype=int, scale=True, label=state.label)
+            self.results += ss.Result(self.name, f'n_{state.name}', self.npts, dtype=int, scale=True, label=state.label)
         return
 
     def step_state(self):
@@ -76,6 +76,10 @@ def step_die(self, uids):
         allows an intervention to avert a death scheduled on the same timestep, without having
         to undo any state changes that have already been applied (because they only run via this
         function if the death actually occurs).
+        
+        Unlike other methods during the integration loop, this method is not called directly
+        by the sim; instead, it is called by people.step_die(), which reconciles the UIDs of
+        the agents who will die.
 
         Depending on the module and the results it produces, it may or may not be necessary
         to implement this.
@@ -120,13 +124,7 @@ def set_prognoses(self, uids, sources=None):
         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)
+        self.log.add_entries(uids, sources, self.now)
         return
 
     def update_results(self):
@@ -137,9 +135,8 @@ def update_results(self):
         This allows result updates at this point to capture outcomes dependent on multiple
         modules, where relevant.
         """
-        sim = self.sim
         for state in self._boolean_states:
-            self.results[f'n_{state.name}'][sim.ti] = np.count_nonzero(state & sim.people.alive)
+            self.results[f'n_{state.name}'][self.ti] = np.count_nonzero(state & self.sim.people.alive)
         return
 
 
@@ -168,7 +165,7 @@ def __init__(self, *args, **kwargs):
     
     def init_pre(self, sim):
         super().init_pre(sim)
-        self.validate_beta()
+        self.validate_beta(run_checks=True)
         return
 
     @property
@@ -189,6 +186,7 @@ def init_post(self):
         super().init_post()
         if self.pars.init_prev is None:
             return
+        
 
         initial_cases = self.pars.init_prev.filter()
         self.set_prognoses(initial_cases)  # TODO: sentinel value to indicate seeds?
@@ -199,25 +197,41 @@ def init_results(self):
         Initialize results
         """
         super().init_results()
-        sim = self.sim
         self.results += [
-            ss.Result(self.name, 'prevalence',     sim.npts, dtype=float, scale=False, label='Prevalence'),
-            ss.Result(self.name, 'new_infections', sim.npts, dtype=int, scale=True, label='New infections'),
-            ss.Result(self.name, 'cum_infections', sim.npts, dtype=int, scale=True, label='Cumulative infections'),
+            ss.Result(self.name, 'prevalence',     self.npts, dtype=float, scale=False, label='Prevalence'),
+            ss.Result(self.name, 'new_infections', self.npts, dtype=int, scale=True, label='New infections'),
+            ss.Result(self.name, 'cum_infections', self.npts, dtype=int, scale=True, label='Cumulative infections'),
         ]
         return
         
-    def validate_beta(self):
+    def validate_beta(self, run_checks=False):
         """ Validate beta and return as a map to match the networks """
         sim = self.sim
+        β = self.pars.beta
+        
+        def scalar_beta(β):
+            return isinstance(β, ss.TimePar) or sc.isnumber(β)
         
-        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 run_checks:
+            scalar_warn = f'Beta is defined as a number ({β}); convert it to a rate to handle timestep conversions'
+
+            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 sc.isnumber(β):
+                ss.warn(scalar_warn)
+            elif isinstance(β, dict):
+                for netbeta in β.values():
+                    if sc.isnumber(netbeta):
+                        ss.warn(scalar_warn)
+                    elif isinstance(netbeta, (list, tuple)):
+                        for thisbeta in netbeta:
+                            if sc.isnumber(netbeta):
+                                ss.warn(scalar_warn)
 
         # If beta is a scalar, apply this bi-directionally to all networks
-        β = self.pars.beta
-        if sc.isnumber(β):
+        if scalar_beta(β):
             betamap = {ss.standardize_netkey(k):[β,β] for k in sim.networks.keys()}
 
         # If beta is a dict, check all entries are bi-directional
@@ -225,11 +239,15 @@ def validate_beta(self):
             betamap = dict()
             for k,thisbeta in β.items():
                 nkey = ss.standardize_netkey(k)
-                if sc.isnumber(thisbeta):
+                if scalar_beta(thisbeta):
                     betamap[nkey] = [thisbeta, thisbeta]
                 else:
                     betamap[nkey] = thisbeta
         
+        else:
+            errormsg = f'Invalid type {type(β)} for beta'
+            raise TypeError(errormsg)
+        
         # Check that it matches the network
         netkeys = [ss.standardize_netkey(k) for k in list(sim.networks.keys())]
         if set(betamap.keys()) != set(netkeys):
@@ -275,8 +293,8 @@ def infect(self):
                 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)
+                        # Calculate probability of a->b transmission
+                        beta_per_dt = net.net_beta(disease_beta=beta) # TODO: potentially refactor
                         p_transmit = rel_trans[src] * rel_sus[trg] * beta_per_dt
         
                         # Generate a new random number based on the two other random numbers
@@ -298,6 +316,7 @@ def infect(self):
             new_cases = ss.uids()
             sources = ss.uids()
             networks = np.empty(0, dtype=ss_int_)
+
         return new_cases, sources, networks
 
     def set_outcomes(self, uids, sources=None):
@@ -316,10 +335,10 @@ def set_congenital(self, uids, sources=None):
     def update_results(self):
         super().update_results()
         res = self.results
-        ti = self.sim.ti
+        ti = self.ti
         res.prevalence[ti] = res.n_infected[ti] / np.count_nonzero(self.sim.people.alive)
         res.new_infections[ti] = np.count_nonzero(self.ti_infected == ti)
-        res.cum_infections[ti] = np.sum(res['new_infections'][:ti+1])
+        res.cum_infections[ti] = np.sum(res['new_infections'][:ti+1]) # TODO: can compute at end
         return
 
 
@@ -351,13 +370,13 @@ class InfectionLog(nx.MultiDiGraph):
 
     A table of outcomes can be returned using `InfectionLog.line_list()`
     """
-    def add_entries(self, sim, uids, sources=None): # TODO: reconcile with other methods
+    def add_entries(self, uids, sources=None, time=np.nan):
         if sources is None:
             for target in uids:
-                self.log.append(np.nan, target, sim.year)
+                self.append(np.nan, target, time)
         else:
             for target, source in zip(uids, sources):
-                self.log.append(source, target, sim.year)
+                self.append(source, target, time)
         return
 
     def add_data(self, uids, **kwargs):
diff --git a/starsim/diseases/cholera.py b/starsim/diseases/cholera.py
index 9fead94e..d65cd9a1 100644
--- a/starsim/diseases/cholera.py
+++ b/starsim/diseases/cholera.py
@@ -19,23 +19,23 @@ def __init__(self, pars=None, *args, **kwargs):
         super().__init__()
         self.define_pars(
             # Initial conditions and beta
-            beta = 1.0, # Placeholder value
+            beta = ss.beta(1.0), # Placeholder value
             init_prev = ss.bernoulli(0.005),
             
             # Natural history parameters, all specified in days
-            dur_exp2inf   = ss.lognorm_ex(mean=2.772, stdev=4.737), # Calculated from Azman et al. estimates https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3677557/
-            dur_asymp2rec = ss.uniform(low=1, high=10), # From WHO cholera fact sheet, asymptomatic individuals shed bacteria for 1-10 days (https://www.who.int/news-room/fact-sheets/detail/cholera)
-            dur_symp2rec  = ss.lognorm_ex(mean=5, stdev=1.8), # According to Fung most modelling studies assume 5 days of symptoms (https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3926264/), but found a range of 2.9-14 days. Distribution approximately fit to these values
-            dur_symp2dead = ss.lognorm_ex(mean=1, stdev=0.5), # There does not appear to be differences in timing/duration of mild vs severe disease, but death from severe disease happens rapidly https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5767916/
+            dur_exp2inf   = ss.days(ss.lognorm_ex(mean=2.772, std=4.737)), # Calculated from Azman et al. estimates https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3677557/
+            dur_asymp2rec = ss.days(ss.uniform(low=1, high=10)), # From WHO cholera fact sheet, asymptomatic individuals shed bacteria for 1-10 days (https://www.who.int/news-room/fact-sheets/detail/cholera)
+            dur_symp2rec  = ss.days(ss.lognorm_ex(mean=5, std=1.8)), # According to Fung most modelling studies assume 5 days of symptoms (https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3926264/), but found a range of 2.9-14 days. Distribution approximately fit to these values
+            dur_symp2dead = ss.days(ss.lognorm_ex(mean=1, std=0.5)), # There does not appear to be differences in timing/duration of mild vs severe disease, but death from severe disease happens rapidly https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5767916/
             p_death       = ss.bernoulli(p=0.005), # Probability of death is typically less than 1% when treated
             p_symp        = ss.bernoulli(p=0.5), # Proportion of infected which are symptomatic, mid range of ~25% and 57% estimates from Jaclson et al (https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3795095/) and Nelson et al (https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3842031/), respectively
             asymp_trans   = 0.01, # Reduction in transmission probability for asymptomatic infection, asymptomatic carriers shed 100-1000 times less bacteria than symptomatic carriers (https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3084143/ and https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3842031/). Previous models assume a 10% relative transmissibility (https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4238032/)
 
             # Environmental parameters
-            beta_env = 0.5 / 3, # Scaling factor for transmission from environment,
+            beta_env = ss.beta(0.5 / 3), # Scaling factor for transmission from environment,
             half_sat_rate = 1_000_000, # Infectious dose in water sufficient to produce infection in 50% of  exposed, from Mukandavire et al. (https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3102413/)
-            shedding_rate = 10, # Rate at which infectious people shed bacteria to the environment (per day), from Mukandavire et al. (https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3102413/)
-            decay_rate = 0.033, # Rate at which bacteria in the environment dies (per day), from Chao et al. and Mukandavire et al. citing https://pubmed.ncbi.nlm.nih.gov/8882180/
+            shedding_rate = ss.rate(10, unit='day'), # Rate at which infectious people shed bacteria to the environment (per day), from Mukandavire et al. (https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3102413/)
+            decay_rate = ss.rate(0.033, unit='day'), # Rate at which bacteria in the environment dies (per day), from Chao et al. and Mukandavire et al. citing https://pubmed.ncbi.nlm.nih.gov/8882180/
             p_env_transmit = ss.bernoulli(p=0), # Probability of environmental transmission - filled out later
         )
         self.update_pars(pars, **kwargs)
@@ -64,16 +64,13 @@ def asymptomatic(self):
         return self.infected & ~self.symptomatic
 
     def init_results(self):
-        """
-        Initialize results
-        """
+        """ Initialize results """
         super().init_results()
-        npts = self.sim.npts
         self.results += [
-            ss.Result(self.name, 'new_deaths', npts, dtype=int, label='Deaths'),
-            ss.Result(self.name, 'cum_deaths', npts, dtype=int, label='Cumulative deaths'),
-            ss.Result(self.name, 'env_prev', npts, dtype=float, label='Environmental prevalence'),
-            ss.Result(self.name, 'env_conc', npts, dtype=float, label='Environmental concentration'),
+            ss.Result(self.name, 'new_deaths', self.npts, dtype=int,   label='Deaths'),
+            ss.Result(self.name, 'cum_deaths', self.npts, dtype=int,   label='Cumulative deaths'),
+            ss.Result(self.name, 'env_prev',   self.npts, dtype=float, label='Environmental prevalence'),
+            ss.Result(self.name, 'env_conc',   self.npts, dtype=float, label='Environmental concentration'),
         ]
         return
 
@@ -83,7 +80,7 @@ def step_state(self):
         Original version by Dom Delport
         """
         # Progress exposed -> infected
-        ti = self.sim.ti
+        ti = self.ti
         infected = (self.exposed & (self.ti_infected <= ti)).uids
         self.infected[infected] = True
 
@@ -114,7 +111,7 @@ def calc_environmental_prev(self):
         """
         p = self.pars
         r = self.results
-        ti = self.sim.ti
+        ti = self.ti
 
         n_symptomatic = self.symptomatic.sum()
         n_asymptomatic = self.asymptomatic.sum()
@@ -130,8 +127,7 @@ def calc_environmental_prev(self):
     def set_prognoses(self, uids, source_uids=None):
         """ Set prognoses for those who get infected """
         super().set_prognoses(uids, source_uids)
-        ti = self.sim.ti
-        dt = self.sim.dt
+        ti = self.ti
 
         self.susceptible[uids] = False
         self.exposed[uids] = True
@@ -140,7 +136,7 @@ def set_prognoses(self, uids, source_uids=None):
         p = self.pars
 
         # Determine when exposed become infected
-        self.ti_infected[uids] = ti + p.dur_exp2inf.rvs(uids) / dt
+        self.ti_infected[uids] = ti + p.dur_exp2inf.rvs(uids)
 
         # Determine who becomes symptomatic and when
         symp_uids = p.p_symp.filter(uids)
@@ -148,13 +144,13 @@ def set_prognoses(self, uids, source_uids=None):
 
         # Determine who dies and when
         dead_uids = p.p_death.filter(symp_uids)
-        self.ti_dead[dead_uids] = self.ti_symptomatic[dead_uids] + p.dur_symp2dead.rvs(dead_uids) / dt
+        self.ti_dead[dead_uids] = self.ti_symptomatic[dead_uids] + p.dur_symp2dead.rvs(dead_uids)
         symp_rev_uids = np.setdiff1d(symp_uids, dead_uids)
         asymp_uids = np.setdiff1d(uids, symp_uids)
 
         # Determine when agents recover
-        self.ti_recovered[symp_rev_uids] = self.ti_exposed[symp_rev_uids] + p.dur_symp2rec.rvs(symp_rev_uids) / dt
-        self.ti_recovered[asymp_uids] = self.ti_exposed[asymp_uids] + p.dur_asymp2rec.rvs(asymp_uids) / dt
+        self.ti_recovered[symp_rev_uids] = self.ti_exposed[symp_rev_uids] + p.dur_symp2rec.rvs(symp_rev_uids)
+        self.ti_recovered[asymp_uids] = self.ti_exposed[asymp_uids] + p.dur_asymp2rec.rvs(asymp_uids)
 
         return
 
@@ -166,7 +162,7 @@ def make_new_cases(self):
         # Make new cases via indirect transmission
         pars = self.pars
         res = self.results
-        p_transmit = res.env_conc[self.sim.ti] * pars.beta_env
+        p_transmit = res.env_conc[self.ti] * pars.beta_env
         pars.p_env_transmit.set(p=p_transmit)
         new_cases = pars.p_env_transmit.filter(self.sim.people.uid[self.susceptible]) # TODO: make syntax nicer
         if new_cases.any():
@@ -182,7 +178,7 @@ def step_die(self, uids):
     def update_results(self):
         super().update_results()
         res = self.results
-        ti = self.sim.ti
+        ti = self.ti
         res.new_deaths[ti] = np.count_nonzero(self.ti_dead == ti)
         res.cum_deaths[ti] = np.sum(res.new_deaths[:ti+1])
         return
diff --git a/starsim/diseases/ebola.py b/starsim/diseases/ebola.py
index aaffcade..a46dba2f 100644
--- a/starsim/diseases/ebola.py
+++ b/starsim/diseases/ebola.py
@@ -19,17 +19,17 @@ def __init__(self, pars=None, *args, **kwargs):
         self.define_pars(
             # Initial conditions and beta
             init_prev       = ss.bernoulli(p=0.005),
-            beta            = 1.0, # Placeholder value
+            beta            = ss.beta(1.0), # Placeholder value
             sev_factor      = 2.2,
             unburied_factor = 2.1,
             
             # Natural history parameters, all specified in days
-            dur_exp2symp    = ss.lognorm_ex(mean=12.7), # Add source
-            dur_symp2sev    = ss.lognorm_ex(mean=6), # Add source
-            dur_sev2dead    = ss.lognorm_ex(mean=1.5), # Add source
-            dur_dead2buried = ss.lognorm_ex(mean=2), # Add source
-            dur_symp2rec    = ss.lognorm_ex(mean=10), # Add source
-            dur_sev2rec     = ss.lognorm_ex(mean=10.4), # Add source
+            dur_exp2symp    = ss.lognorm_ex(mean=ss.days(12.7)), # Add source
+            dur_symp2sev    = ss.lognorm_ex(mean=ss.days(6)), # Add source
+            dur_sev2dead    = ss.lognorm_ex(mean=ss.days(1.5)), # Add source
+            dur_dead2buried = ss.lognorm_ex(mean=ss.days(2)), # Add source
+            dur_symp2rec    = ss.lognorm_ex(mean=ss.days(10)), # Add source
+            dur_sev2rec     = ss.lognorm_ex(mean=ss.days(10.4)), # Add source
             p_sev           = ss.bernoulli(p=0.7), # Add source
             p_death         = ss.bernoulli(p=0.55), # Add source
             p_safe_bury     = ss.bernoulli(p=0.25), # Probability of a safe burial - should be linked to diagnoses
@@ -60,7 +60,7 @@ def infectious(self):
     def step_state(self):
 
         # Progress exposed -> infected
-        ti = self.sim.ti
+        ti = self.ti
         infected = (self.exposed & (self.ti_infected <= ti)).uids
         self.exposed[infected] = False
         self.infected[infected] = True
@@ -95,8 +95,7 @@ def set_prognoses(self, uids, source_uids=None):
         # Do not call set_prognoses on the parent
         #super().set_prognoses(sim, uids, source_uids)
 
-        ti = self.sim.ti
-        dt = self.sim.dt
+        ti = self.ti
         self.susceptible[uids] = False
         self.exposed[uids] = True
         self.ti_exposed[uids] = ti
@@ -104,25 +103,25 @@ def set_prognoses(self, uids, source_uids=None):
         p = self.pars
 
         # Determine when exposed become infected
-        self.ti_infected[uids] = ti + p.dur_exp2symp.rvs(uids) / dt
+        self.ti_infected[uids] = ti + p.dur_exp2symp.rvs(uids)
 
         # Determine who progresses to sever and when
         sev_uids = p.p_sev.filter(uids)
-        self.ti_severe[sev_uids] = self.ti_infected[sev_uids] + p.dur_symp2sev.rvs(sev_uids) / dt
+        self.ti_severe[sev_uids] = self.ti_infected[sev_uids] + p.dur_symp2sev.rvs(sev_uids)
 
         # Determine who dies and who recovers and when
         dead_uids = p.p_death.filter(sev_uids)
-        self.ti_dead[dead_uids] = self.ti_severe[dead_uids] + p.dur_sev2dead.rvs(dead_uids) / dt
+        self.ti_dead[dead_uids] = self.ti_severe[dead_uids] + p.dur_sev2dead.rvs(dead_uids)
         rec_sev_uids = np.setdiff1d(sev_uids, dead_uids)
-        self.ti_recovered[rec_sev_uids] = self.ti_severe[rec_sev_uids] + p.dur_sev2rec.rvs(rec_sev_uids) / dt
+        self.ti_recovered[rec_sev_uids] = self.ti_severe[rec_sev_uids] + p.dur_sev2rec.rvs(rec_sev_uids)
         rec_symp_uids = np.setdiff1d(uids, sev_uids)
-        self.ti_recovered[rec_symp_uids] = self.ti_infected[rec_symp_uids] + p.dur_symp2rec.rvs(rec_symp_uids) / dt
+        self.ti_recovered[rec_symp_uids] = self.ti_infected[rec_symp_uids] + p.dur_symp2rec.rvs(rec_symp_uids)
 
         # Determine time of burial - either immediate (safe burials) or after a delay (unsafe)
         safe_buried = p.p_safe_bury.filter(dead_uids)
         unsafe_buried = np.setdiff1d(dead_uids, safe_buried)
         self.ti_buried[safe_buried] = self.ti_dead[safe_buried]
-        self.ti_buried[unsafe_buried] = self.ti_dead[unsafe_buried] + p.dur_dead2buried.rvs(unsafe_buried) / dt
+        self.ti_buried[unsafe_buried] = self.ti_dead[unsafe_buried] + p.dur_dead2buried.rvs(unsafe_buried)
 
         # Change rel_trans values
         self.rel_trans[self.infectious] = 1
diff --git a/starsim/diseases/gonorrhea.py b/starsim/diseases/gonorrhea.py
index 3202aaa2..ab2a6347 100644
--- a/starsim/diseases/gonorrhea.py
+++ b/starsim/diseases/gonorrhea.py
@@ -14,8 +14,8 @@ def __init__(self, pars=None, *args, **kwargs):
         # Parameters
         super().__init__()
         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
+            beta = ss.beta(1.0), # Placeholder value
+            dur_inf_in_days = ss.lognorm_ex(mean=ss.days(10), std=ss.days(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
             p_clear   = ss.bernoulli(p=0.2),  # Share of infections that spontaneously clear: https://sti.bmj.com/content/96/8/556
             init_prev = ss.bernoulli(p=0.1),
@@ -36,30 +36,30 @@ def __init__(self, pars=None, *args, **kwargs):
     def init_results(self):
         """ Initialize results """
         super().init_results()
-        self.results += ss.Result(self.name, 'new_clearances', self.sim.npts, dtype=int)
+        self.results += ss.Result(self.name, 'new_clearances', self.npts, dtype=int)
         return
 
     def update_results(self):
         super().update_results()
-        ti = self.sim.ti
+        ti = self.ti
         self.results.n_symptomatic[ti] = self.symptomatic.count()
         self.results.new_clearances[ti] = np.count_nonzero(self.ti_clearance == ti)
         return
 
     def step_state(self):
         """ Natural clearance """
-        clearances = self.ti_clearance <= self.sim.ti
+        clearances = self.ti_clearance <= self.ti
         self.susceptible[clearances] = True
         self.infected[clearances] = False
         self.symptomatic[clearances] = False
-        self.ti_clearance[clearances] = self.sim.ti
+        self.ti_clearance[clearances] = self.ti
 
         return
 
     def set_prognoses(self, uids, source_uids=None):
         """ Natural history of gonorrhea for adult infection """
         super().set_prognoses(uids, source_uids)
-        ti = self.sim.ti
+        ti = self.ti
 
         # Set infection status
         self.susceptible[uids] = False
@@ -72,6 +72,6 @@ def set_prognoses(self, uids, source_uids=None):
 
         # Set natural clearance
         clear_uids = self.pars.p_clear.filter(uids)
-        dur = ti + self.pars.dur_inf_in_days.rvs(clear_uids)/365/self.sim.dt # Convert from days to years and then adjust for dt
+        dur = ti + self.pars.dur_inf_in_days.rvs(clear_uids) # Convert from days to years and then adjust for dt automatically # TODO: check!
         self.ti_clearance[clear_uids] = dur
         return
\ No newline at end of file
diff --git a/starsim/diseases/hiv.py b/starsim/diseases/hiv.py
index b35ce610..565e0dfd 100644
--- a/starsim/diseases/hiv.py
+++ b/starsim/diseases/hiv.py
@@ -14,7 +14,8 @@ class HIV(ss.Infection):
     def __init__(self, pars=None, *args, **kwargs):
         super().__init__()
         self.define_pars(
-            beta = 1.0, # Placeholder value
+            unit = 'year',
+            beta = ss.beta(1.0), # Placeholder value
             cd4_min = 100,
             cd4_max = 500,
             cd4_rate = 5,
@@ -22,7 +23,7 @@ def __init__(self, pars=None, *args, **kwargs):
             art_efficacy = 0.96,
             init_prev = ss.bernoulli(p=0.05),
             death_dist = ss.bernoulli(p=self.death_prob_func), # Uses p_death by default, modulated by CD4
-            p_death = 0.05,
+            p_death = ss.rate(0.05), # NB: this is death per unit time, not death per infection
         )
         self.update_pars(pars=pars, **kwargs)
 
@@ -38,7 +39,7 @@ def __init__(self, pars=None, *args, **kwargs):
     @staticmethod
     def death_prob_func(module, sim, uids):
         p = module.pars
-        out = sim.dt * p.p_death / (p.cd4_min - p.cd4_max)**2 *  (module.cd4[uids] - p.cd4_max)**2
+        out = p.p_death / (p.cd4_min - p.cd4_max)**2 *  (module.cd4[uids] - p.cd4_max)**2
         out = np.array(out)
         return out
 
@@ -54,18 +55,18 @@ def step_state(self):
         hiv_deaths = self.pars.death_dist.filter(can_die)
         
         people.request_death(hiv_deaths)
-        self.ti_dead[hiv_deaths] = self.sim.ti
+        self.ti_dead[hiv_deaths] = self.ti
         return
 
     def init_results(self):
         """ Initialize results """
         super().init_results()
-        self.results += ss.Result(self.name, 'new_deaths', self.sim.npts, dtype=int, label='Deaths')
+        self.results += ss.Result(self.name, 'new_deaths', self.npts, dtype=int, label='Deaths')
         return
 
     def update_results(self):
         super().update_results()
-        ti = self.sim.ti
+        ti = self.ti
         self.results['new_deaths'][ti] = np.count_nonzero(self.ti_dead == ti)
         return 
 
@@ -73,7 +74,7 @@ def set_prognoses(self, uids, source_uids=None):
         super().set_prognoses(uids, source_uids)
         self.susceptible[uids] = False
         self.infected[uids] = True
-        self.ti_infected[uids] = self.sim.ti
+        self.ti_infected[uids] = self.ti
         return
 
     def set_congenital(self, uids, source_uids):
@@ -90,7 +91,7 @@ def __init__(self, year, coverage, pars=None, **kwargs):
         self.coverage = sc.toarray(coverage)
         super().__init__()
         self.define_pars(
-            art_delay = ss.constant(v=1) # Value in years
+            art_delay = ss.constant(v=ss.years(1.0)) # Value in years
         )
         self.update_pars(pars=pars, **kwargs)
 
@@ -100,7 +101,7 @@ def __init__(self, year, coverage, pars=None, **kwargs):
 
     def init_pre(self, sim):
         super().init_pre(sim)
-        self.results += ss.Result(self.name, 'n_art', sim.npts, dtype=int)
+        self.results += ss.Result(self.name, 'n_art', self.npts, dtype=int)
         self.initialized = True
         return
 
@@ -111,7 +112,7 @@ def step(self):
 
         hiv = sim.people.hiv
         infected = hiv.infected.uids
-        ti_delay = np.round(self.pars.art_delay.rvs(infected)/sim.dt).astype(int)
+        ti_delay = np.round(self.pars.art_delay.rvs(infected)).astype(int)
         recently_infected = infected[hiv.ti_infected[infected] == sim.ti-ti_delay]
 
         n_added = 0
@@ -138,7 +139,7 @@ def __init__(self):
 
     def init_pre(self, sim):
         super().init_pre(sim)
-        self.cd4 = np.zeros((sim.npts, sim.people.n), dtype=int)
+        self.cd4 = np.zeros((self.npts, sim.people.n), dtype=int)
         return
 
     def step(self):
diff --git a/starsim/diseases/measles.py b/starsim/diseases/measles.py
index f6b6d35c..099d4746 100644
--- a/starsim/diseases/measles.py
+++ b/starsim/diseases/measles.py
@@ -21,8 +21,8 @@ def __init__(self, pars=None, *args, **kwargs):
             init_prev = ss.bernoulli(p=0.005),
             
             # Natural history parameters, all specified in days
-            dur_exp = ss.normal(loc=8),        # (days) - source: US CDC
-            dur_inf = ss.normal(loc=11),       # (days) - source: US CDC
+            dur_exp = ss.normal(loc=ss.days(8)),        # (days) - source: US CDC
+            dur_inf = ss.normal(loc=ss.days(11)),       # (days) - source: US CDC
             p_death = ss.bernoulli(p=0.005), # Probability of death
         )
         self.update_pars(pars=pars, **kwargs)
@@ -41,7 +41,7 @@ def infectious(self):
 
     def step_state(self):
         # Progress exposed -> infected
-        ti = self.sim.ti
+        ti = self.ti
         infected = (self.exposed & (self.ti_infected <= ti)).uids
         self.exposed[infected] = False
         self.infected[infected] = True
@@ -60,8 +60,7 @@ def step_state(self):
     def set_prognoses(self, uids, source_uids=None):
         """ Set prognoses for those who get infected """
         super().set_prognoses(uids, source_uids)
-        ti = self.sim.ti
-        dt = self.sim.dt
+        ti = self.ti
 
         self.susceptible[uids] = False
         self.exposed[uids] = True
@@ -70,7 +69,7 @@ def set_prognoses(self, uids, source_uids=None):
         p = self.pars
 
         # Determine when exposed become infected
-        self.ti_infected[uids] = ti + p.dur_exp.rvs(uids) / dt
+        self.ti_infected[uids] = ti + p.dur_exp.rvs(uids)
 
         # Sample duration of infection, being careful to only sample from the
         # distribution once per timestep.
@@ -80,8 +79,8 @@ def set_prognoses(self, uids, source_uids=None):
         will_die = p.p_death.rvs(uids)
         dead_uids = uids[will_die]
         rec_uids = uids[~will_die]
-        self.ti_dead[dead_uids] = self.ti_infected[dead_uids] + dur_inf[will_die] / dt
-        self.ti_recovered[rec_uids] = self.ti_infected[rec_uids] + dur_inf[~will_die] / dt
+        self.ti_dead[dead_uids] = self.ti_infected[dead_uids] + dur_inf[will_die]
+        self.ti_recovered[rec_uids] = self.ti_infected[rec_uids] + dur_inf[~will_die]
 
         return
 
diff --git a/starsim/diseases/ncd.py b/starsim/diseases/ncd.py
index 9c301359..8bb4a9a0 100644
--- a/starsim/diseases/ncd.py
+++ b/starsim/diseases/ncd.py
@@ -21,8 +21,8 @@ def __init__(self, pars=None, **kwargs):
         super().__init__()
         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
+            dur_risk = ss.expon(scale=ss.dur(10)),
+            prognosis = ss.weibull(c=ss.years(2), scale=5), # Time in years between first becoming affected and death
         )
         self.update_pars(pars=pars, **kwargs)
         
@@ -48,11 +48,11 @@ def init_post(self):
         super().init_post()
         initial_risk = self.pars['initial_risk'].filter()
         self.at_risk[initial_risk] = True
-        self.ti_affected[initial_risk] = self.sim.ti + sc.randround(self.pars['dur_risk'].rvs(initial_risk) / self.sim.dt)
+        self.ti_affected[initial_risk] = self.ti + sc.randround(self.pars['dur_risk'].rvs(initial_risk))
         return initial_risk
 
     def step_state(self):
-        ti = self.sim.ti
+        ti = self.ti
         deaths = (self.ti_dead == ti).uids
         self.sim.people.request_death(deaths)
         self.log.add_data(deaths, died=True)
@@ -60,11 +60,11 @@ def step_state(self):
         return
 
     def step(self):
-        ti = self.sim.ti
+        ti = self.ti
         new_cases = (self.ti_affected == ti).uids
         self.affected[new_cases] = True
         prog_years = self.pars.prognosis.rvs(new_cases)
-        self.ti_dead[new_cases] = ti + sc.randround(prog_years / self.sim.dt)
+        self.ti_dead[new_cases] = ti + sc.randround(prog_years / self.dt)
         super().set_prognoses(new_cases)
         return new_cases
 
@@ -73,17 +73,16 @@ def init_results(self):
         Initialize results
         """
         super().init_results()
-        npts = self.sim.npts
         self.results += [
-            ss.Result(self.name, 'n_not_at_risk', npts, dtype=int, label='Not at risk'),
-            ss.Result(self.name, 'prevalence', npts, dtype=float, label='Prevalence'),
-            ss.Result(self.name, 'new_deaths', npts, dtype=int, label='Deaths'),
+            ss.Result(self.name, 'n_not_at_risk', self.npts, dtype=int,   label='Not at risk'),
+            ss.Result(self.name, 'prevalence',    self.npts, dtype=float, label='Prevalence'),
+            ss.Result(self.name, 'new_deaths',    self.npts, dtype=int,   label='Deaths'),
         ]
         return
 
     def update_results(self):
         super().update_results()
-        ti = self.sim.ti
+        ti = self.ti
         self.results.n_not_at_risk[ti] = np.count_nonzero(self.not_at_risk)
         self.results.prevalence[ti]    = np.count_nonzero(self.affected)/len(self.sim.people)
         self.results.new_deaths[ti]    = np.count_nonzero(self.ti_dead == ti)
diff --git a/starsim/diseases/sir.py b/starsim/diseases/sir.py
index 9c0e43e8..af00c42c 100644
--- a/starsim/diseases/sir.py
+++ b/starsim/diseases/sir.py
@@ -18,50 +18,50 @@ 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)
+    # def __future_init__(self, pars=None, **kwargs):
+    #     """ The plan for implementing events """
+    #     super().__init__()
+    #     self.define_pars(
+    #         beta = ss.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'),
-        )
+    #     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 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'),
-        )
+    #     # 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
+    #     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.define_pars(
-            beta = 0.1,
+            beta = ss.beta(0.1),
             init_prev = ss.bernoulli(p=0.01),
-            dur_inf = ss.lognorm_ex(mean=6),
+            dur_inf = ss.lognorm_ex(mean=ss.dur(6)),
             p_death = ss.bernoulli(p=0.01),
         )
         self.update_pars(pars, **kwargs)
@@ -93,8 +93,7 @@ def step_state(self):
 
     def set_prognoses(self, uids, source_uids=None):
         """ Set prognoses """
-        ti = self.sim.ti
-        dt = self.sim.dt
+        ti = self.ti
         self.susceptible[uids] = False
         self.infected[uids] = True
         self.ti_infected[uids] = ti
@@ -109,8 +108,8 @@ def set_prognoses(self, uids, source_uids=None):
         will_die = p.p_death.rvs(uids)
         dead_uids = uids[will_die]
         rec_uids = uids[~will_die]
-        self.ti_dead[dead_uids] = ti + dur_inf[will_die] / dt # Consider rand round, but not CRN safe
-        self.ti_recovered[rec_uids] = ti + dur_inf[~will_die] / dt
+        self.ti_dead[dead_uids] = ti + dur_inf[will_die] # Consider rand round, but not CRN safe
+        self.ti_recovered[rec_uids] = ti + dur_inf[~will_die]
         return
 
     def step_die(self, uids):
@@ -145,10 +144,10 @@ class SIS(ss.Infection):
     def __init__(self, pars=None, *args, **kwargs):
         super().__init__()
         self.define_pars(
-            beta = 0.05,
+            beta = ss.beta(0.05),
             init_prev = ss.bernoulli(p=0.01),
-            dur_inf = ss.lognorm_ex(mean=10),
-            waning = 0.05,
+            dur_inf = ss.lognorm_ex(mean=ss.dur(10)),
+            waning = ss.rate(0.05),
             imm_boost = 1.0,
         )
         self.update_pars(pars=pars, *args, **kwargs)
@@ -161,7 +160,7 @@ def __init__(self, pars=None, *args, **kwargs):
 
     def step_state(self):
         """ Progress infectious -> recovered """
-        recovered = (self.infected & (self.ti_recovered <= self.sim.ti)).uids
+        recovered = (self.infected & (self.ti_recovered <= self.ti)).uids
         self.infected[recovered] = False
         self.susceptible[recovered] = True
         self.update_immunity()
@@ -169,7 +168,7 @@ def step_state(self):
     
     def update_immunity(self):
         has_imm = (self.immunity > 0).uids
-        self.immunity[has_imm] = (self.immunity[has_imm])*(1 - self.pars.waning*self.sim.dt)
+        self.immunity[has_imm] = (self.immunity[has_imm])*(1 - self.pars.waning)
         self.rel_sus[has_imm] = np.maximum(0, 1 - self.immunity[has_imm])
         return
 
@@ -178,27 +177,27 @@ def set_prognoses(self, uids, source_uids=None):
         super().set_prognoses(uids, source_uids)
         self.susceptible[uids] = False
         self.infected[uids] = True
-        self.ti_infected[uids] = self.sim.ti
+        self.ti_infected[uids] = self.ti
         self.immunity[uids] += self.pars.imm_boost
 
         # Sample duration of infection
         dur_inf = self.pars.dur_inf.rvs(uids)
 
         # Determine when people recover
-        self.ti_recovered[uids] = self.sim.ti + dur_inf / self.sim.dt
+        self.ti_recovered[uids] = self.ti + dur_inf
 
         return
     
     def init_results(self):
         """ Initialize results """
         super().init_results()
-        self.results += ss.Result(self.name, 'rel_sus', self.sim.npts, dtype=float, label='Relative susceptibility')
+        self.results += ss.Result(self.name, 'rel_sus', self.npts, dtype=float, label='Relative susceptibility')
         return
 
     def update_results(self):
         """ Store the population immunity (susceptibility) """
         super().update_results()
-        self.results['rel_sus'][self.sim.ti] = self.rel_sus.mean()
+        self.results['rel_sus'][self.ti] = self.rel_sus.mean()
         return 
 
     def plot(self):
diff --git a/starsim/diseases/syphilis.py b/starsim/diseases/syphilis.py
index fd212c38..155e0e5b 100644
--- a/starsim/diseases/syphilis.py
+++ b/starsim/diseases/syphilis.py
@@ -17,15 +17,15 @@ def __future_init__(self, pars=None, **kwargs):
         super().__init__()
         self.define_pars(
             # Initial conditions
-            beta = 1.0, # Placeholder
+            beta = ss.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_exposed = ss.lognorm_ex(mean=1 / 12, std=1 / 36),  # https://pubmed.ncbi.nlm.nih.gov/9101629/
+            dur_primary = ss.lognorm_ex(mean=1.5 / 12, std=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/
+            dur_latent_temp = ss.lognorm_ex(mean=1, std=6 / 12),  # https://pubmed.ncbi.nlm.nih.gov/9101629/
+            dur_latent_long = ss.lognorm_ex(mean=20, std=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/
 
@@ -96,11 +96,11 @@ def __init__(self, pars=None, **kwargs):
             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_exposed = ss.lognorm_ex(mean=1 / 12, std=1 / 36),  # https://pubmed.ncbi.nlm.nih.gov/9101629/
+            dur_primary = ss.lognorm_ex(mean=1.5 / 12, std=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/
+            dur_latent_temp = ss.lognorm_ex(mean=1, std=6 / 12),  # https://pubmed.ncbi.nlm.nih.gov/9101629/
+            dur_latent_long = ss.lognorm_ex(mean=20, std=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/
     
@@ -185,11 +185,10 @@ def infectious(self):
     def init_results(self):
         """ Initialize results """
         super().init_results()
-        npts = self.sim.npts
         self.results += [
-            ss.Result(self.name, 'new_nnds',       npts, dtype=int, scale=True, label='New neonatal deaths'),
-            ss.Result(self.name, 'new_stillborns', npts, dtype=int, scale=True, label='New stillborns'),
-            ss.Result(self.name, 'new_congenital', npts, dtype=int, scale=True, label='New congenital syphilis'),
+            ss.Result(self.name, 'new_nnds',       self.npts, dtype=int, scale=True, label='New neonatal deaths'),
+            ss.Result(self.name, 'new_stillborns', self.npts, dtype=int, scale=True, label='New stillborns'),
+            ss.Result(self.name, 'new_congenital', self.npts, dtype=int, scale=True, label='New congenital syphilis'),
         ]
         return
 
@@ -197,7 +196,7 @@ def step_state(self):
         """ Updates prior to interventions """
 
         # Primary
-        ti = self.sim.ti
+        ti = self.ti
         primary = self.exposed & (self.ti_primary <= ti)
         self.primary[primary] = True
         self.exposed[primary] = False
@@ -211,10 +210,6 @@ def step_state(self):
             self.set_secondary_prognoses(secondary_from_primary)
             self.rel_trans[secondary_from_primary] = self.pars.rel_trans['secondary']
 
-        # Hack to reset the MultiRNGs in set_secondary_prognoses so that they can be called again in this timestep. TODO: Refactor
-        self.pars.p_latent_temp.jump(ti+1)
-        self.pars.dur_secondary.jump(ti+1)
-
         # Secondary reactivation from latent
         secondary_from_latent = (self.latent_temp & (self.ti_secondary <= ti)).uids
         if len(secondary_from_latent) > 0:
@@ -260,7 +255,7 @@ def step_state(self):
 
     def update_results(self):
         super().update_results()
-        ti = self.sim.ti
+        ti = self.ti
         self.results.new_nnds[ti]       = np.count_nonzero(self.ti_nnd == ti)
         self.results.new_stillborns[ti] = np.count_nonzero(self.ti_stillborn == ti)
         self.results.new_congenital[ti] = np.count_nonzero(self.ti_congenital == ti)
@@ -271,8 +266,7 @@ def set_prognoses(self, uids, source_uids=None):
         Set initial prognoses for adults newly infected with syphilis
         """
         
-        ti = self.sim.ti
-        dt = self.sim.dt
+        ti = self.ti
 
         self.susceptible[uids] = False
         self.ever_exposed[uids] = True
@@ -284,17 +278,16 @@ def set_prognoses(self, uids, source_uids=None):
         # Set future dates and probabilities
         # Exposed to primary
         dur_exposed = self.pars.dur_exposed.rvs(uids)
-        self.ti_primary[uids] = ti + rr(dur_exposed / dt)
+        self.ti_primary[uids] = ti + rr(dur_exposed)
 
         # Primary to secondary
         dur_primary = self.pars.dur_primary.rvs(uids)
-        self.ti_secondary[uids] = self.ti_primary[uids] + rr(dur_primary / dt)
+        self.ti_secondary[uids] = self.ti_primary[uids] + rr(dur_primary)
         return
 
     def set_secondary_prognoses(self, uids):
         """ Set prognoses for people who have just progressed to secondary infection """
 
-        dt = self.sim.dt
         dur_secondary = self.pars.dur_secondary.rvs(uids)
 
         # Secondary to latent_temp or latent_long
@@ -303,31 +296,30 @@ def set_secondary_prognoses(self, uids):
         latent_long_uids = uids[~latent_temp]
 
         dur_secondary_temp = dur_secondary[latent_temp]
-        self.ti_latent_temp[latent_temp_uids] = self.ti_secondary[latent_temp_uids] + rr(dur_secondary_temp / dt)
+        self.ti_latent_temp[latent_temp_uids] = self.ti_secondary[latent_temp_uids] + rr(dur_secondary_temp)
 
         dur_secondary_long = dur_secondary[~latent_temp]
-        self.ti_latent_long[latent_long_uids] = self.ti_secondary[latent_long_uids] + rr(dur_secondary_long / dt)
+        self.ti_latent_long[latent_long_uids] = self.ti_secondary[latent_long_uids] + rr(dur_secondary_long)
 
         return
 
     def set_latent_temp_prognoses(self, uids):
         # Primary to secondary
         dur_latent_temp = self.pars.dur_latent_temp.rvs(uids)
-        self.ti_secondary[uids] = self.ti_latent_temp[uids] + rr(dur_latent_temp / self.sim.dt)
+        self.ti_secondary[uids] = self.ti_latent_temp[uids] + rr(dur_latent_temp)
         return
 
     def set_latent_long_prognoses(self, uids):
-        dt = self.sim.dt
         dur_latent = self.pars.dur_latent_long.rvs(uids)
 
         # Primary to secondary
         dur_latent_long = dur_latent
-        self.ti_secondary[uids] = self.ti_latent_temp[uids] + rr(dur_latent_long / dt)
+        self.ti_secondary[uids] = self.ti_latent_temp[uids] + rr(dur_latent_long)
 
         # Latent_long to tertiary
         tertiary = self.pars.p_tertiary.rvs(uids)
         tertiary_uids = uids[tertiary]
-        self.ti_tertiary[tertiary_uids] = self.ti_latent_long[tertiary_uids] + rr(dur_latent_long[tertiary] / dt)
+        self.ti_tertiary[tertiary_uids] = self.ti_latent_long[tertiary_uids] + rr(dur_latent_long[tertiary])
 
         return
 
@@ -349,12 +341,13 @@ def set_congenital(self, uids, source_uids=None):
                 time_to_birth = -sim.people.age.raw # TODO: make nicer
 
                 # Schedule events
+                ratio = ss.time_ratio(unit1='year', dt1=1.0, unit2=self.unit, dt2=self.dt) # TODO: think about simplifying
                 for oi, outcome in enumerate(self.pars.birth_outcome_keys):
                     o_uids = state_uids[assigned_outcomes == oi]
                     if len(o_uids) > 0:
                         ti_outcome = f'ti_{outcome}'
                         vals = getattr(self, ti_outcome)
-                        vals[o_uids] = sim.ti + rr(time_to_birth[o_uids] / sim.dt)
+                        vals[o_uids] = sim.ti + rr(time_to_birth[o_uids] * ratio)
                         setattr(self, ti_outcome, vals)
 
         return
@@ -423,8 +416,8 @@ def check_eligibility(self):
     def init_pre(self, sim):
         super().init_pre(sim)
         self.results += [
-            ss.Result('syphilis', 'n_screened', sim.npts, dtype=int, scale=True, label='Number screened'),
-            ss.Result('syphilis', 'n_dx', sim.npts, dtype=int, scale=True, label='Number diagnosed'),
+            ss.Result('syphilis', 'n_screened', self.npts, dtype=int, scale=True, label='Number screened'),
+            ss.Result('syphilis', 'n_dx',       self.npts, dtype=int, scale=True, label='Number diagnosed'),
         ]
         return
 
@@ -446,7 +439,7 @@ def _parse_product_str(self, product):
 
     def init_pre(self, sim):
         super().init_pre(sim)
-        self.results += ss.Result('syphilis', 'n_tx', sim.npts, dtype=int, scale=True, label='Number treated')
+        self.results += ss.Result('syphilis', 'n_tx', self.npts, dtype=int, scale=True, label='Number treated')
         return
 
     def step(self):
diff --git a/starsim/distributions.py b/starsim/distributions.py
index e360de0a..a413d14b 100644
--- a/starsim/distributions.py
+++ b/starsim/distributions.py
@@ -8,42 +8,33 @@
 import starsim as ss
 import matplotlib.pyplot as pl
 
-__all__ = ['find_dists', 'link_dists', 'make_dist', 'dist_list', 'Dists', 'Dist']
+__all__ = ['link_dists', 'make_dist', 'dist_list', 'Dists', 'Dist']
 
 
 def str2int(string, modulo=1_000_000):
     """
     Convert a string to an int
     
-    Cannot use Python's built-in hash() since it's randomized for strings, but
-    this is almost as fast (and 5x faster than hashlib).
+    Cannot use Python's built-in hash() since it's randomized for strings. Hashlib
+    is 5x slower than int.from_bytes(string.encode(), byteorder='big'), but should
+    only add a couple milliseconds to a typical sim.
     """
-    return int.from_bytes(string.encode(), byteorder='big') % modulo
-
-
-def find_dists(obj, verbose=False, **kwargs):
-    """ Find all Dist objects in a parent object """
-    out = sc.objdict()
-    tree = sc.iterobj(obj, depthfirst=False, flatten=True, **kwargs)
-    if verbose: print(f'Found {len(tree)} objects')
-    for trace,val in tree.items():
-        if isinstance(val, Dist):
-            out[trace] = val
-            if verbose: print(f'  {trace} is a dist ({len(out)})')
+    integer = int.from_bytes(string.encode(), byteorder='big')
+    # integer = sc.sha(string, asint=True) # Hash the string to an integer
+    out = integer % modulo # Don't need all of it, this is more user-friendly
     return out
 
 
-def link_dists(obj, sim, module=None, overwrite=False, init=False, **kwargs):
-    """ Link distributions to the sim and the module; used in module.initialize() and people.initialize() """
+def link_dists(obj, sim, module=None, overwrite=False, init=False, **kwargs): # TODO: actually link the distributions to the modules! Currently this only does the opposite, but should have mod.dists as well
+    """ Link distributions to the sim and the module; used in module.init() and people.init() """
     if module is None and isinstance(obj, ss.Module):
         module = obj
-    dists = ss.find_dists(obj, **kwargs) # Important that this comes first, before the sim is linked to the dist!
-    for key,val in dists.items():
-        if isinstance(val, ss.Dist):
-            val.link_sim(sim, overwrite=overwrite)
-            val.link_module(module, overwrite=overwrite)
-            if init: # Usually this is false since usually these are initialized centrally by the sim
-                val.initialize()
+    dists = sc.search(obj, type=Dist, **kwargs)
+    for dist in dists.values():
+        dist.link_sim(sim, overwrite=overwrite)
+        dist.link_module(module, overwrite=overwrite)
+        if init: # Usually this is false since usually these are initialized centrally by the sim
+            dist.init()
     return
 
 
@@ -75,7 +66,7 @@ def __init__(self, obj=None, *args, base_seed=None, sim=None):
         self.initialized = False
         return
 
-    def initialize(self, obj=None, base_seed=None, sim=None, force=False):
+    def init(self, obj=None, base_seed=None, sim=None, force=False):
         """
         Set the base seed, find and initialize all distributions in an object
         
@@ -97,10 +88,10 @@ def initialize(self, obj=None, base_seed=None, sim=None, force=False):
         skip = id(sim.people._states) if sim is not None else None
         
         # Find and initialize the distributions
-        self.dists = find_dists(obj, skip=skip)
+        self.dists = sc.search(obj, type=Dist, skip=skip, flatten=True)
         for trace,dist in self.dists.items():
             if not dist.initialized or force:
-                dist.initialize(trace=trace, seed=base_seed, sim=sim, force=force)
+                dist.init(trace=trace, seed=base_seed, sim=sim, force=force)
         
         # Confirm the seeds are unique
         self.check_seeds()
@@ -118,8 +109,8 @@ def check_seeds(self):
                 checked[seed] = dist
         return
 
-    def jump(self, to=None, delta=1):
-        """ Advance all RNGs, e.g. to timestep "to", by jumping """
+    def jump(self, to=None, delta=1, force=False):
+        """ Advance all RNGs, e.g. to call "to", by jumping """
         out = sc.autolist()
 
         # Do not jump if centralized
@@ -127,7 +118,24 @@ def jump(self, to=None, delta=1):
             return out
 
         for dist in self.dists.values():
-            out += dist.jump(to=to, delta=delta)
+            out += dist.jump(to=to, delta=delta, force=force)
+        return out
+    
+    def jump_dt(self, ti=None, force=False): # Could this be simplified with jump(), or nice to have the parallel with Dist?
+        """
+        Advance all RNGs to the next timestep
+        
+        Args:
+            ti (int): if specified, jump to this timestep (default: current sim timestep)
+        """
+        out = sc.autolist()
+        
+        # Do not jump if centralized
+        if ss.options._centralized:
+            return out
+        
+        for dist in self.dists.values():
+            out += dist.jump_dt(ti=ti, force=force)
         return out
 
     def reset(self):
@@ -179,7 +187,7 @@ def __init__(self, dist=None, distname=None, name=None, seed=None, offset=None,
         self.dist = dist # The type of distribution
         self.distname = distname
         self.name = name
-        self.pars = sc.dictobj(kwargs) # The user-defined kwargs
+        self.pars = sc.objdict(kwargs) # The user-defined kwargs
         self.seed = seed # Usually determined once added to the container
         self.offset = offset
         self.module = module
@@ -199,6 +207,7 @@ def __init__(self, dist=None, distname=None, name=None, seed=None, offset=None,
         self._slots = None # Internal variable to track currently-in-use slots
         
         # History and random state
+        self.dt_jump_size = 1000 # How much to advance the RNG for each timestep (must be larger than the number times dist is called per timestep)
         self.rng = None # The actual RNG generator for generating random numbers
         self.trace = None # The path of this object within the parent
         self.ind = 0 # The index of the RNG (usually updated on each timestep)
@@ -207,7 +216,7 @@ def __init__(self, dist=None, distname=None, name=None, seed=None, offset=None,
         self.ready = True
         self.initialized = False
         if not strict: # Otherwise, wait for a sim
-            self.initialize()
+            self.init()
         return
     
     def __repr__(self):
@@ -304,7 +313,7 @@ def reset(self, state=0):
         
         **Example**::
             
-            dist = ss.random(seed=5).initialize()
+            dist = ss.random(seed=5).init()
             r1 = dist(5)
             r2 = dist(5)
             dist.reset(-1)
@@ -321,21 +330,41 @@ def reset(self, state=0):
         self.ready = True
         return self.state
 
-    def jump(self, to=None, delta=1):
+    def jump(self, to=None, delta=1, force=False):
         """ Advance the RNG, e.g. to timestep "to", by jumping """
         
-        # Do not jump if centralized
+        # Do not jump if centralized # TODO: remove
         if ss.options._centralized:
             return self.state
-
+        
+        # Validation
         jumps = to if (to is not None) else self.ind + delta
+        if self.ind >= jumps and not force:
+            errormsg = f'You tried to jump the distribution "{self}" to state {jumps}, but the ' \
+                       f'RNG state is already at state {self.ind}, meaning you will draw the same ' \
+                        'random numbers twice. If you are sure you want to do this, set force=True.'
+            raise DistSeedRepeatError(msg=errormsg)
+        
+        # Do the jumping
         self.ind = jumps
         self.reset() # First reset back to the initial state (used in case of different numbers of calls)
         if jumps: # Seems to randomize state if jumps=0
             self.bitgen.state = self.bitgen.jumped(jumps=jumps).state # Now take "jumps" number of jumps
         return self.state
     
-    def initialize(self, trace=None, seed=None, module=None, sim=None, slots=None, force=False):
+    def jump_dt(self, ti=None, force=False):
+        """
+        Automatically jump on the next value of dt
+        
+        Args:
+            ti (int): if specified, jump to this timestep (default: current sim timestep plus one)
+        """
+        if ti is None:
+            ti = self.sim.ti + 1
+        to = self.dt_jump_size*ti
+        return self.jump(to=to, force=force)
+    
+    def init(self, trace=None, seed=None, module=None, sim=None, slots=None, force=False):
         """ Calculate the starting seed and create the RNG """
         
         if self.initialized is True and not force: # Don't warn if we have a partially initialized distribution
@@ -380,7 +409,7 @@ def initialize(self, trace=None, seed=None, module=None, sim=None, slots=None, f
         return self
     
     def link_sim(self, sim=None, overwrite=False):
-        """ Shortcut for linking the sim, only overwriting an existing one if overwrite=True """
+        """ Shortcut for linking the sim, only overwriting an existing one if overwrite=True; not for the user """
         if (not self.sim or overwrite) and sim is not None:
             self.sim = sim
         return
@@ -392,7 +421,7 @@ def link_module(self, module=None, overwrite=False):
         return
     
     def process_seed(self, trace=None, seed=None):
-        """ Obtain the seed offset by hashing the path to this distribution; called automatically """
+        """ Obtain the seed offset by hashing the path to this distribution; not for the user """
         unique_name = trace or self.trace or self.name
         if unique_name:
             if not self.name:
@@ -405,7 +434,7 @@ def process_seed(self, trace=None, seed=None):
         return
     
     def process_dist(self):
-        """ Ensure the distribution works """
+        """ Ensure the distribution works; not for the user """
         
         # Handle a SciPy distribution, if provided
         if self.dist is not None:
@@ -430,7 +459,7 @@ def process_dist(self):
         return
     
     def process_size(self, n=1):
-        """ Handle an input of either size or UIDs and calculate size, UIDs, and slots """
+        """ Handle an input of either size or UIDs and calculate size, UIDs, and slots; not for the user """
         if np.isscalar(n) or isinstance(n, tuple):  # If passing a non-scalar size, interpret as dimension rather than UIDs iff a tuple
             uids = None
             slots = None
@@ -454,7 +483,8 @@ def process_size(self, n=1):
         return size, slots
     
     def process_pars(self, call=True):
-        """ Ensure the supplied dist and parameters are valid, and initialize them; called automatically """
+        """ Ensure the supplied dist and parameters are valid, and initialize them; not for the user """
+        self._time_factor = None # Time rescalings need to be done after distributions are calculated; store the correction factor here 
         self._pars = sc.cp(self.pars) # The actual keywords; shallow copy, modified below for special cases
         if call:
             self.call_pars() # Convert from function to values if needed
@@ -462,7 +492,7 @@ def process_pars(self, call=True):
         return spars
     
     def call_pars(self):
-        """ Check if any parameters need to be called to be turned into arrays """
+        """ Check if any parameters need to be called to be turned into arrays; not for the user """
         
         # Initialize
         size, uids = self._size, self._uids
@@ -472,11 +502,21 @@ def call_pars(self):
         # Check each parameter
         for key,val in self._pars.items():
             
+            # If it's a time parameter, transform it to a float now
+            if isinstance(val, ss.TimePar):
+                self._pars[key] = val.value # Use the raw value
+                factor = val.x/val.value # Calculate the ratio here; NB, this is usually val.factor, but not for ss.time_prob
+                if self._time_factor is None:
+                    self._time_factor = factor
+                elif factor != self._time_factor:
+                    errormsg = f'Cannot have time parameters in the same distribution with inconsistent unit/dt values: {self._pars}'
+                    raise ValueError(errormsg)                    
+            
             # If the parameter is callable, then call it
-            if callable(val) and not isinstance(val, type): # Types can appear as callable
+            elif callable(val) and not isinstance(val, type): # Types can appear as callable
                 size_par = uids if uids is not None else size
                 out = val(self.module, self.sim, size_par) # TODO: swap order to sim, module, size?
-                val = np.asarray(out) # Necessary since UIDArrays don't allow slicing # TODO: check if this is correct
+                val = np.asarray(out) # Necessary since FloatArrs don't allow slicing # TODO: check if this is correct
                 self._pars[key] = val
             
             # If it's iterable and UIDs are provided, then we need to use array-parameter logic
@@ -485,12 +525,12 @@ def call_pars(self):
         return
     
     def sync_pars(self):
-        """ Perform any necessary synchronizations or transformations on distribution parameters """
+        """ Perform any necessary synchronizations or transformations on distribution parameters; not for the user """
         self.update_dist_pars()
         return self._pars
     
     def update_dist_pars(self, pars=None):
-        """ Update SciPy distribution parameters """
+        """ Update SciPy distribution parameters; not for the user """
         if self.dist is not None:
             pars = pars if pars is not None else self._pars
             self.dist.kwds = pars
@@ -501,7 +541,7 @@ def rand(self, size):
         return self.rng.random(size)
     
     def make_rvs(self):
-        """ Return default random numbers for scalar parameters """
+        """ Return default random numbers for scalar parameters; not for the user """
         if self.rvs_func is not None:
             rvs_func = getattr(self.rng, self.rvs_func) # Can't store this because then it references the wrong RNG after copy
             rvs = rvs_func(**self._pars, size=self._size)
@@ -513,13 +553,13 @@ def make_rvs(self):
         return rvs
     
     def ppf(self, rands):
-        """ Return default random numbers for array parameters """
+        """ Return default random numbers for array parameters; not for the user """
         rvs = self.dist.ppf(rands)
         return rvs
     
     def rvs(self, n=1, reset=False):
         """
-        Get random variables
+        Get random variables -- use this!
         
         Args:
             n (int/tuple/arr): if an int or tuple, return this many random variables; if an array, treat as UIDs
@@ -555,6 +595,10 @@ def rvs(self, n=1, reset=False):
             rvs = self.make_rvs() # Or, just get regular values
             if self._slots is not None:
                 rvs = rvs[self._slots]
+                
+        # Scale by time if needed
+        if self._time_factor is not None:
+            rvs = rvs*self._time_factor # TODO: Can't use rvs *= in case it's not a float; this is not ideal though since if not a float (e.g. for Poisson) it's probably a bug
         
         # Tidy up
         self.called += 1
@@ -645,7 +689,7 @@ def __init__(self, loc=0.0, scale=1.0, **kwargs):
 class lognorm_im(Dist):
     """
     Lognormal distribution, parameterized in terms of the "implicit" (normal)
-    distribution, with mean=loc and stdev=scale (see lognorm_ex for comparison).
+    distribution, with mean=loc and std=scale (see lognorm_ex for comparison).
     
     Note: the "loc" parameter here does *not* correspond to the mean of the resulting
     random variates!
@@ -678,20 +722,20 @@ def sync_pars(self, call=True):
 class lognorm_ex(Dist):
     """
     Lognormal distribution, parameterized in terms of the "explicit" (lognormal)
-    distribution, with mean=mean and stdev=stdev for this distribution (see lognorm_im for comparison).
+    distribution, with mean=mean and std=std for this distribution (see lognorm_im for comparison).
     Note that a mean ≤ 0.0 is impossible, since this is the parameter of the distribution
     after the log transform.
     
     Args:
         mean (float): the mean of this distribution (not the underlying distribution) (default 1.0)
-        stdev (float): the standard deviation of this distribution (not the underlying distribution) (default 1.0)
+        std (float): the standard deviation of this distribution (not the underlying distribution) (default 1.0)
     
     **Example**::
         
-        ss.lognorm_ex(mean=2, stdev=1, strict=False).rvs(1000).mean() # Should be close to 2
+        ss.lognorm_ex(mean=2, std=1, strict=False).rvs(1000).mean() # Should be close to 2
     """
-    def __init__(self, mean=1.0, stdev=1.0, **kwargs):
-        super().__init__(distname='lognormal', dist=sps.lognorm, mean=mean, stdev=stdev, **kwargs)
+    def __init__(self, mean=1.0, std=1.0, **kwargs):
+        super().__init__(distname='lognormal', dist=sps.lognorm, mean=mean, std=std, **kwargs)
         return
     
     def convert_ex_to_im(self):
@@ -705,13 +749,13 @@ def convert_ex_to_im(self):
         self.call_pars() # Since can't work with functions
         p = self._pars
         mean = p.pop('mean')
-        stdev = p.pop('stdev')
+        std = p.pop('std')
         if np.isscalar(mean) and mean <= 0:
             errormsg = f'Cannot create a lognorm_ex distribution with mean≤0 (mean={mean}); did you mean to use lognorm_im instead?'
             raise ValueError(errormsg)
-        std2 = stdev**2
+        std2 = std**2
         mean2 = mean**2
-        sigma_im = np.sqrt(np.log(std2/mean2 + 1)) # Computes stdev for the underlying normal distribution
+        sigma_im = np.sqrt(np.log(std2/mean2 + 1)) # Computes std for the underlying normal distribution
         mean_im  = np.log(mean2 / np.sqrt(std2 + mean2)) # Computes the mean of the underlying normal distribution
         p.mean = mean_im
         p.sigma = sigma_im
@@ -737,7 +781,7 @@ def __init__(self, scale=1.0, **kwargs):
         return
 
 
-class poisson(Dist):
+class poisson(Dist): # TODO: does not currently scale correctly with dt
     """
     Poisson distribution
     
@@ -790,6 +834,7 @@ def ppf(self, rands):
         rvs = rvs.astype(self.pars['dtype'])
         return rvs
 
+
 class rand_raw(Dist):
     """
     Directly sample raw integers (uint64) from the random number generator.
@@ -996,15 +1041,18 @@ def __init__(self, names, *args, **kwargs):
     def __len__(self):
         return len(self.dists)
         
-    def initialize(self, *args, **kwargs):
-        for dist in self.dists: dist.initialize(*args, **kwargs)
+    def init(self, *args, **kwargs):
+        """ Not usually needed since each dist will handle this automatically; for completeness only """
+        for dist in self.dists: dist.init(*args, **kwargs)
         return
     
     def reset(self, *args, **kwargs):
+        """ Not usually needed since each dist will handle this automatically; for completeness only """
         for dist in self.dists: dist.reset(*args, **kwargs)
         return
     
     def jump(self, *args, **kwargs):
+        """ Not usually needed since each dist will handle this automatically; for completeness only """
         for dist in self.dists: dist.jump(*args, **kwargs)
         return
     
@@ -1036,23 +1084,26 @@ def rvs(self, *args):
 
 class DistNotInitializedError(RuntimeError):
     """ Raised when Dist object is called when not initialized. """
-    def __init__(self, dist):
-        msg = f'{dist} has not been initialized; please set strict=False when creating the distribution, or call dist.initialize()'
+    def __init__(self, dist=None, msg=None):
+        if msg is None:
+            msg = f'{dist} has not been initialized; please set strict=False when creating the distribution, or call dist.init()'
         super().__init__(msg)
         return
 
 
 class DistNotReadyError(RuntimeError):
     """ Raised when a Dist object is called without being ready. """
-    def __init__(self, dist):
-        msg = f'{dist} is not ready. This is likely caused by calling a distribution multiple times in a single step. Call dist.jump() to reset.'
+    def __init__(self, dist=None, msg=None):
+        if msg is None:
+            msg = f'{dist} is not ready. This is likely caused by calling a distribution multiple times in a single step. Call dist.jump() to reset.'
         super().__init__(msg)
         return
 
 
 class DistSeedRepeatError(RuntimeError):
     """ Raised when a Dist object shares a seed with another """
-    def __init__(self, dist1, dist2):
-        msg = f'A common seed was found between {dist1} and {dist2}. This is likely caused by incorrect initialization of the parent Dists object.'
+    def __init__(self, dist1=None, dist2=None, msg=None):
+        if msg is None:
+            msg = f'A common seed was found between {dist1} and {dist2}. This is likely caused by incorrect initialization of the parent Dists object.'
         super().__init__(msg)
         return
diff --git a/starsim/interventions.py b/starsim/interventions.py
index 509b6758..11101aff 100644
--- a/starsim/interventions.py
+++ b/starsim/interventions.py
@@ -83,26 +83,27 @@ def init_pre(self, sim):
 
         # If start_year and end_year are not provided, figure them out from the provided years or the sim
         if self.years is None:
-            if self.start_year is None: self.start_year = sim.pars['start']
-            if self.end_year is None:   self.end_year = sim.pars['end']
+            if self.start_year is None: self.start_year = sim.pars.start
+            if self.end_year is None:   self.end_year = sim.pars.stop
         else:
             self.start_year = self.years[0]
             self.end_year = self.years[-1]
 
         # More validation
-        if not(any(np.isclose(self.start_year, sim.yearvec)) and any(np.isclose(self.end_year, sim.yearvec))):
+        if not(any(np.isclose(self.start_year, sim.timevec)) and any(np.isclose(self.end_year, sim.timevec))):
             errormsg = 'Years must be within simulation start and end dates.'
             raise ValueError(errormsg)
 
         # Adjustment to get the right end point
-        adj_factor = int(1 / sim.dt) - 1 if sim.dt < 1 else 1
+        dt = sim.pars.dt # TODO: need to eventually replace with own timestep, but not initialized yet since super().init_pre() hasn't been called
+        adj_factor = int(1/dt) - 1 if dt < 1 else 1
 
         # Determine the timepoints at which the intervention will be applied
-        self.start_point = sc.findfirst(sim.yearvec, self.start_year)
-        self.end_point   = sc.findfirst(sim.yearvec, self.end_year) + adj_factor
+        self.start_point = sc.findfirst(sim.timevec, self.start_year)
+        self.end_point   = sc.findfirst(sim.timevec, self.end_year) + adj_factor
         self.years       = sc.inclusiverange(self.start_year, self.end_year)
         self.timepoints  = sc.inclusiverange(self.start_point, self.end_point)
-        self.yearvec     = np.arange(self.start_year, self.end_year + adj_factor, sim.dt)
+        self.yearvec     = np.arange(self.start_year, self.end_year + adj_factor, dt)
 
         # Get the probability input into a format compatible with timepoints
         if len(self.years) != len(self.prob):
@@ -115,7 +116,7 @@ def init_pre(self, sim):
             self.prob = sc.smoothinterp(self.yearvec, self.years, self.prob, smoothness=0)
 
         # Lastly, adjust the probability by the sim's timestep, if it's an annual probability
-        if self.annual_prob: self.prob = 1 - (1 - self.prob) ** sim.dt
+        if self.annual_prob: self.prob = 1 - (1 - self.prob) ** dt
 
         return
 
@@ -134,7 +135,7 @@ def __init__(self, years, interpolate=None, prob=None, *args, **kwargs):
 
     def init_pre(self, sim):
         # Decide whether to apply the intervention at every timepoint throughout the year, or just once.
-        self.timepoints = sc.findnearest(sim.yearvec, self.years)
+        self.timepoints = sc.findnearest(sim.timevec, self.years)
 
         if len(self.prob) == 1:
             self.prob = np.array([self.prob[0]] * len(self.timepoints))
@@ -213,7 +214,7 @@ def step(self):
         """
         sim = self.sim
         accept_uids = ss.uids()
-        if sim.ti in self.timepoints:
+        if sim.ti in self.timepoints: # TODO: change to self.ti
             accept_uids = self.deliver()
             self.screened[accept_uids] = True
             self.screens[accept_uids] += 1
diff --git a/starsim/loop.py b/starsim/loop.py
new file mode 100644
index 00000000..a4e1fa1a
--- /dev/null
+++ b/starsim/loop.py
@@ -0,0 +1,257 @@
+"""
+Parent class for the integration loop.
+"""
+
+import numpy as np
+import pandas as pd
+import sciris as sc
+import starsim as ss
+import matplotlib.pyplot as plt
+
+# What classes are externally visible
+__all__ = ['Loop']
+
+
+
+#%% Loop class
+
+class Loop:
+    """ Base class for integration loop """
+    def __init__(self, sim): # TODO: consider eps=1e-6 and round times to this value
+        self.sim = sim
+        self.funcs = None
+        self.timearrays = None
+        self.plan = None
+        self.index = 0 # The next function to execute
+        self.initialized = False
+        return
+    
+    def init(self):
+        """ Parse the sim into the integration plan """
+        self.collect_funcs()
+        self.collect_timearrays()
+        self.make_plan()
+        self.initialized = True
+        return
+    
+    def __len__(self):
+        if self.initialized:
+            return len(self.plan)
+        else:
+            return 0 # Or None?
+    
+    def __iadd__(self, func):
+        """ Allow functions to be added to the function list """
+        parent = func.__self__
+        func_name = func.__name__
+        
+        # Get the name if it's defined, the class otherwise; these must match timearrays
+        module = parent.name if isinstance(parent, ss.Module) else parent.__class__.__name__.lower()
+        
+        # Create the row and append it to the function list
+        row = dict(func_order=len(self.funcs), module=module, func_name=func_name, func=func)
+        self.funcs.append(row)
+        return self
+    
+    def __repr__(self):
+        if self.initialized:
+            arrs = list({len(arr) for arr in self.timearrays.values()})
+            if len(arrs) == 1: arrs = arrs[0] # If all are the same, just use that
+            string = f'Loop(n={len(self)}, funcs={len(self.funcs)}, npts={arrs}, index={self.index})'
+        else:
+            string = 'Loop(initialized=False)'
+        return string
+    
+    def disp(self):
+        return sc.pr(self)
+
+    def collect_funcs(self):
+        """ Collect all the callable functions (methods) that comprise the step """
+        
+        # Run the simulation step first (updates the distributions)
+        self.funcs = [] # Reset, just in case
+        sim = self.sim
+        self += sim.start_step # Note special __iadd__() method above, which appends these to the funcs list
+        
+        # Update demographic modules (create new agents from births/immigration, schedule non-disease deaths and emigration)
+        for dem in sim.demographics():
+            self += 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 sim.diseases():
+            if isinstance(disease, ss.Disease): # Could be a connector instead -- TODO, rethink this
+                self += disease.step_state
+
+        # Update connectors
+        for connector in sim.connectors():
+            self += connector.step
+
+        # Update networks - this takes place here in case autonomous state changes at this timestep
+        # affect eligibility for contacts
+        for network in sim.networks():
+            self += 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 intv in sim.interventions():
+            self += intv.step
+        
+        # Carry out autonomous state changes in the disease modules, including transmission (but excluding deaths)
+        for disease in sim.diseases():
+            self += disease.step
+
+        # Update people who died -- calls disease.step_die() internally
+        self += sim.people.step_die
+
+        # Update results
+        self += sim.people.update_results
+        for mod in sim.modules:
+            self += mod.update_results
+
+        # Apply analyzers
+        for ana in sim.analyzers():
+            self += ana.step
+            
+        # Clean up dead agents, increment the time index, and perform other housekeeping tasks
+        for mod in sim.modules:
+            self += mod.finish_step
+        self += sim.people.finish_step
+        self += sim.finish_step
+        
+        return self.funcs
+    
+    def collect_timearrays(self):
+        """ Collect numerical time arrays for each module """
+        self.timearrays = sc.objdict()
+        
+        # Handle the sim and people first
+        sim = self.sim
+        for key in ['sim', 'people']:
+            self.timearrays[key] = sim.timearray
+        
+        # Handle all other modules
+        for mod in sim.modules:
+            timearray = ss.make_timearray(mod.timevec, mod.unit, sim.pars.unit)
+            self.timearrays[mod.name] = timearray
+            
+        return self.timearrays
+    
+    def make_plan(self):
+        """ Combine the module ordering and the time vectors into the integration plan """
+        
+        # Assemble the list of dicts
+        raw = []
+        for func_row in self.funcs:
+            for time in self.timearrays[func_row['module']]:
+                row = func_row.copy()
+                row['time'] = time # Add time column
+                raw.append(row)
+        
+        # Turn it into a dataframe
+        self.plan = sc.dataframe(raw)
+        
+        # Sort it by step_order, a combination of time and function order
+        self.plan['step_order'] = self.plan.time + ss.options.time_eps*self.plan.func_order
+        self.plan['func_label'] = self.plan.module + '.' + self.plan.func_name
+        col_order = ['time', 'func_order', 'step_order', 'func', 'func_label', 'module', 'func_name'] # Func in the middle to hide it
+        self.plan = self.plan.sort_values('step_order').reset_index(drop=True)[col_order]
+        return
+    
+    def run_one_step(self):
+        """
+        Take a single step, i.e. call a single function; only used for debugging purposes.
+        
+        Compare sim.run_one_step(), which runs a full timestep (which involves multiple function calls).        
+        """
+        f = self.plan.func[self.index] # Get the next function
+        f() # Call it
+        self.index += 1 # Increment the time
+        return
+    
+    def run(self, until=None, verbose=None):
+        """ Actually run the integration loop; usually called by sim.run() """
+        # Loop over every function in the integration loop, e.g. disease.step()
+        for f in self.plan.func[self.index:]:
+            if verbose:
+                row = self.plan[self.index]
+                print(f'Running t={row.time:n}, {row.func_label}()')
+            
+            f() # Execute the function -- this is where all of Starsim happens!!
+
+            # Tidy up
+            self.index += 1 # Increment the count
+            if until and self.sim.now > until: # Terminate if asked to
+                break
+        return
+    
+    def manual_reset(self): # TODO: do we need this? I feel if we don't have it, people will be tempted to manually set loop.index = 0.
+        """
+        Reset the loop to run again. Note, does not reset sim quantities so should 
+        only be used for debugging.
+        """
+        self.index = 0
+        self.sim.complete = False
+        return
+    
+    def to_df(self):
+        """ Return a user-friendly version of the plan, omitting object columns """
+        cols = ['time', 'func_order', 'module', 'func_name', 'func_label']
+        df = self.plan[cols]
+        return df
+    
+    def plot(self, simplify=False, fig_kw=None, plot_kw=None, scatter_kw=None):
+        """
+        Plot a diagram of all the events
+        
+        Args:
+            simplify (bool): if True, skip update_results and finish_step events, which are automatically applied
+            fig_kw (dict): passed to ``plt.figure()``
+            plot_kw (dict): passed to ``plt.plot()``
+            scatter_kw (dict): passed to ``plt.scatter()``
+        """
+        
+        # Assemble data
+        df = self.to_df()
+        if simplify:
+            filter_out = ['update_results', 'finish_step']
+            df = df[~df.func_name.isin(filter_out)]
+        yticks = df.func_order.unique()
+        ylabels = df.func_label.unique()
+        x = df.time
+        y = df.func_order
+        
+        # Convert module names to integers for plotting colors
+        mod_int, _ = pd.factorize(df.module) 
+        colors = sc.gridcolors(np.unique(mod_int), asarray=True)
+        
+        # Do the plotting
+        plot_kw = sc.mergedicts(dict(lw=2, alpha=0.2, c='k'), plot_kw)
+        scatter_kw = sc.mergedicts(dict(s=200, alpha=0.6), scatter_kw)
+        fig = plt.figure(**sc.mergedicts(fig_kw))
+        plt.plot(x, y, **plot_kw)
+        plt.scatter(x, y, c=colors[mod_int], **scatter_kw)
+        plt.yticks(yticks, ylabels)
+        plt.title(f'Integration plan ({len(df)} events)')
+        plt.xlabel(f'Time since simulation start (in {self.sim.pars.unit}s)')
+        plt.grid(True)
+        sc.figlayout()
+        sc.boxoff()
+        return fig
+    
+    def __deepcopy__(self, memo):
+        """ A dataframe that has functions in it doesn't copy well; convert to a dict first """
+        cls = self.__class__
+        new = cls.__new__(cls)
+        memo[id(self)] = new
+        for k, v in self.__dict__.items():
+            if k == 'plan':
+                origdict = v.to_dict() # Convert to a dictionary
+                newdict = sc.dcp(origdict, memo=memo) # Copy the dict
+                newdf = sc.dataframe(newdict)
+                setattr(new, k, newdf)
+            else:
+                setattr(new, k, sc.dcp(v, memo=memo)) # Regular deepcopy
+        return new
+        
+    
\ No newline at end of file
diff --git a/starsim/modules.py b/starsim/modules.py
index 4aec8122..000a34ca 100644
--- a/starsim/modules.py
+++ b/starsim/modules.py
@@ -45,10 +45,12 @@ def find_modules(key=None):
 
 class Module(sc.quickobj):
 
-    def __init__(self, name=None, label=None):
-        self.set_metadata(name, label) # Usually reset as part of self.update_pars()
+    def __init__(self, name=None, label=None, unit=None, dt=None):
         self.pars = ss.Pars() # Usually populated via self.define_pars()
+        self.set_metadata(name, label) # Usually reset as part of self.update_pars()
+        self.set_time_pars(unit, dt)
         self.results = ss.Results(self.name)
+        self.pre_initialized = False
         self.initialized = False
         self.finalized = False
         return
@@ -69,10 +71,16 @@ def disp(self, output=False):
         else:
             return out
     
-    def set_metadata(self, name, label):
+    def set_metadata(self, name=None, label=None):
         """ 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.name  = sc.ifelse(name,  getattr(self, 'name',  self.pars.get('name', self.__class__.__name__.lower()))) # Default name is the class name
+        self.label = sc.ifelse(label, getattr(self, 'label', self.pars.get('label', self.name)))
+        return
+    
+    def set_time_pars(self, unit=None, dt=None):
+        """ Set time units for the module """
+        self.unit  = sc.ifelse(unit,  getattr(self, 'unit', self.pars.get('unit')))
+        self.dt    = sc.ifelse(dt,    getattr(self, 'dt',   self.pars.get('dt')))
         return
     
     def define_pars(self, inherit=True, **kwargs): # TODO: think if inherit should default to true or false
@@ -96,7 +104,9 @@ def update_pars(self, pars, **kwargs):
                 
         # Update module attributes
         metadata = {key:pars.pop(key, None) for key in ['name', 'label']}
+        timepars = {key:pars.pop(key, None) for key in ['unit', 'dt']}
         self.set_metadata(**metadata)
+        self.set_time_pars(**timepars)
         
         # Should be no remaining pars
         if len(pars):
@@ -104,7 +114,7 @@ def update_pars(self, pars, **kwargs):
             raise ValueError(errormsg)
         return
     
-    def init_pre(self, sim):
+    def init_pre(self, sim, force=False):
         """
         Perform initialization steps
 
@@ -112,12 +122,15 @@ def init_pre(self, sim):
         initialization, initialized=False until init_vals() is called (which is after
         distributions are initialized).
         """
-        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()
+        if force or not self.pre_initialized:
+            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_time_pars() # Initialize the modules' time parameters and link them to the sim
+            self.init_results()
+            self.pre_initialized = True
         return
     
     def init_results(self):
@@ -132,10 +145,48 @@ def init_post(self):
         self.initialized = True
         return
     
+    def init_time_pars(self, force=False):
+        """ Initialize all time parameters by ensuring all parameters are initialized; part of init_post() """
+        pars = self.sim.pars
+        
+        # Find all modules and set the timestep
+        if force or self.unit is None:
+            self.unit = pars.unit
+        if force or self.dt is None:
+            self.dt = pars.dt
+        
+        # Find all time parameters in the module
+        timepars = sc.search(self.pars, type=ss.TimePar) # Should it be self or self.pars?
+        
+        # Initialize them with the parent module
+        for timepar in timepars.values():
+            if force or not timepar.initialized:
+                timepar.init(parent=self)
+        
+        # Create the module-specific time vector
+        self.timevec = ss.make_timevec(pars.start, pars.stop, self.dt, self.unit)
+        self.npts = len(self.timevec)
+        self.ti = 0 # Track the current timestep, which may or may not match the sim's
+        return
+    
+    @property
+    def now(self):
+        """ Return the current time, i.e. the time vector at the current timestep """
+        try:
+            return self.timevec[self.ti]
+        except Exception as E:
+            ss.warn(f'Encountered exception when getting the current time in {self.name}: {E}')
+            return None
+    
     def step(self):
         """ Define how the module updates over time """
         pass
     
+    def finish_step(self):
+        """ Define what should happen at the end of the step; at minimum, increment ti """
+        self.ti += 1
+        return
+    
     def update_results(self):
         """ Perform any results updates on each timestep """
         pass
@@ -220,6 +271,8 @@ def step(mod): # TODO: see if this can be done more simply
         mod = cls(name=name)
         mod.func = func
         mod.step = partial(step, mod)
+        mod.step.__name__ = name # Manually add these in as for a regular class method
+        mod.step.__self__ = mod 
         return mod
             
     def to_json(self):
@@ -232,12 +285,13 @@ def to_json(self):
         return out
 
     def plot(self):
+        """ Plot all results in the module """
         with sc.options.with_style('fancy'):
             flat = sc.flattendict(self.results, sep=': ')
-            yearvec = self.sim.yearvec
+            timevec = self.timevec
             fig, axs = sc.getrowscols(len(flat), make=True)
             for ax, (k, v) in zip(axs.flatten(), flat.items()):
-                ax.plot(yearvec, v)
+                ax.plot(timevec, v)
                 ax.set_title(k)
                 ax.set_xlabel('Year')
         return fig
diff --git a/starsim/network.py b/starsim/network.py
index 8982bcff..ae5cc0d5 100644
--- a/starsim/network.py
+++ b/starsim/network.py
@@ -114,11 +114,14 @@ def __len__(self):
 
     def __repr__(self, **kwargs):
         """ Convert to a dataframe for printing """
-        namestr = self.name
-        labelstr = f'"{self.label}"' if self.label else '<no label>'
-        keys_str = ', '.join(self.edges.keys())
-        output = f'{namestr}({labelstr}, {keys_str})\n'  # e.g. Network("r", p1, p2, beta)
-        output += self.to_df().__repr__()
+        try:
+            namestr = self.name
+            labelstr = f'"{self.label}"' if self.label else '<no label>'
+            keys_str = ', '.join(self.edges.keys())
+            output = f'{namestr}({labelstr}, {keys_str})\n'  # e.g. Network("r", p1, p2, beta)
+            output += self.to_df().__repr__()
+        except:
+            output = sc.prepr(self, vals=False)
         return output
 
     def __contains__(self, item):
@@ -300,9 +303,9 @@ def remove_uids(self, uids):
 
         return
 
-    def beta_per_dt(self, disease_beta=None, dt=None, uids=None):
+    def net_beta(self, disease_beta=None, dt=None, uids=None):
         if uids is None: uids = Ellipsis
-        return self.edges.beta[uids] * disease_beta * dt
+        return self.edges.beta[uids] * disease_beta # Beta should already include dt if desired
 
 
 class DynamicNetwork(Network):
@@ -311,10 +314,15 @@ def __init__(self, key_dict=None, **kwargs):
         key_dict = sc.mergedicts({'dur': ss_float_}, key_dict)
         super().__init__(key_dict=key_dict, **kwargs)
         return
+    
+    def step(self):
+        self.end_pairs()
+        self.add_pairs()
+        return
 
     def end_pairs(self):
         people = self.sim.people
-        self.edges.dur = self.edges.dur - self.sim.dt
+        self.edges.dur = self.edges.dur - self.dt # TODO: think about whether this is right # Update: it is, if duration is *NOT* a ss.dur! Otherwise it should be -1, in timestep units
 
         # Non-alive agents are removed
         active = (self.edges.dur > 0) & people.alive[self.edges.p1] & people.alive[self.edges.p2]
@@ -348,9 +356,9 @@ def available(self, people, sex):
         available[self.edges.p2] = False
         return available.uids
 
-    def beta_per_dt(self, disease_beta=None, dt=None, uids=None):
+    def net_beta(self, disease_beta=None, uids=None):
         if uids is None: uids = Ellipsis
-        return self.edges.beta[uids] * (1 - (1 - disease_beta) ** (self.edges.acts[uids] * dt))
+        return self.edges.beta[uids] * (1 - (1 - disease_beta) ** (self.edges.acts[uids] * self.dt))
 
 
 # %% Specific instances of networks
@@ -442,7 +450,7 @@ def __init__(self, pars=None, key_dict=None, **kwargs):
         super().__init__(key_dict=key_dict)
         self.define_pars(
             n_contacts = ss.constant(10),
-            dur = 0,
+            dur = ss.dur(0),
         )
         self.update_pars(pars, **kwargs)
         self.dist = ss.Dist(distname='RandomNet') # Default RNG
@@ -492,11 +500,6 @@ def get_edges(self, inds, n_contacts):
         self.dist.jump() # Reset the RNG manually # TODO, think if there's a better way
         return source, target
 
-    def step(self):
-        self.end_pairs()
-        self.add_pairs()
-        return
-
     def add_pairs(self):
         """ Generate edges """
         people = self.sim.people
@@ -514,7 +517,7 @@ def add_pairs(self):
         if isinstance(self.pars.dur, ss.Dist):
             dur = self.pars.dur.rvs(p1)
         else:
-            dur = np.full(len(p1), self.pars.dur)
+            dur = np.ones(len(p1))*self.pars.dur # Other option would be np.full(len(p1), self.pars.dur.x), but this is harder to read
         
         self.append(p1=p1, p2=p2, beta=beta, dur=dur)
         return
@@ -551,11 +554,6 @@ def init_post(self):
         self.add_pairs()
         return
 
-    def update(self):
-        self.end_pairs()
-        self.add_pairs()
-        return
-
     def add_pairs(self):
         """ Generate contacts """
         people = self.sim.people
@@ -620,7 +618,7 @@ def init_post(self):
     def step(self):
 
         # Motion step
-        vdt = self.pars.v * self.sim.dt
+        vdt = self.pars.v * self.dt
         self.x[:] = self.x + vdt * np.cos(self.theta)
         self.y[:] = self.y + vdt * np.sin(self.theta)
 
@@ -712,7 +710,7 @@ def __init__(self, pars=None, key_dict=None, **kwargs):
             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
-            acts = ss.poisson(lam=80),
+            acts = ss.poisson(lam=80), # TODO: make this work with ss.rate, which it currently does not due to network initialization limitations
             rel_part_rates = 1.0,
         )
         self.update_pars(pars=pars, **kwargs)
@@ -787,7 +785,7 @@ def add_pairs(self):
 
     def step(self):
         self.end_pairs()
-        self.set_network_states(upper_age=self.sim.dt) # TODO: looks wrong
+        self.set_network_states(upper_age=self.dt) # TODO: check
         self.add_pairs()
         return
 
@@ -800,9 +798,9 @@ class MSMNet(SexualNetwork):
     def __init__(self, pars=None, key_dict=None, **kwargs):
         super().__init__(key_dict=key_dict)
         self.define_pars(
-            duration = ss.lognorm_ex(mean=15, stdev=15),
+            duration = ss.lognorm_ex(mean=15, std=15),
             debut = ss.normal(loc=16, scale=2),
-            acts = ss.lognorm_ex(mean=80, stdev=20),
+            acts = ss.lognorm_ex(mean=80, std=20),
             participation = ss.bernoulli(p=0.1),
         )
         self.update_pars(pars, **kwargs)
@@ -931,15 +929,18 @@ def step(self):
         """
         Set beta to 0 for women who complete duration of transmission
         Keep connections for now, might want to consider removing
+        
+        NB: add_pairs() and end_pairs() are NOT called here; this is done separately
+        in ss.Pregnancy.update_states().
         """
-        inactive = self.edges.end <= self.sim.ti
+        inactive = self.edges.end <= self.ti
         self.edges.beta[inactive] = 0
         return
 
     def end_pairs(self):
         people = self.sim.people
         edges = self.edges
-        active = (edges.end > self.sim.ti) & people.alive[edges.p1] & people.alive[edges.p2]
+        active = (edges.end > self.ti) & people.alive[edges.p1] & people.alive[edges.p2]
         for k in self.meta_keys():
             edges[k] = edges[k][active]
         return len(active)
@@ -950,10 +951,10 @@ def add_pairs(self, mother_inds=None, unborn_inds=None, dur=None, start=None):
             return 0
         else:
             if start is None:
-                start = np.full_like(dur, fill_value=self.sim.ti)
+                start = np.full_like(dur, fill_value=self.ti)
             n = len(mother_inds)
             beta = np.ones(n)
-            end = start + sc.promotetoarray(dur) / self.sim.dt
+            end = start + sc.promotetoarray(dur)
             self.append(p1=mother_inds, p2=unborn_inds, beta=beta, dur=dur, start=start, end=end)
             return n
 
diff --git a/starsim/parameters.py b/starsim/parameters.py
index 193c816b..c02cefc7 100644
--- a/starsim/parameters.py
+++ b/starsim/parameters.py
@@ -56,12 +56,15 @@ def update(self, pars=None, create=False, **kwargs):
                 old = self[key] # Get the existing object we're about to update
                 if isinstance(old, atomic_classes): # It's a number, string, etc: update directly
                     self[key] = new
+
                 elif isinstance(old, Pars): # It's a Pars object: update recursively
                     old.update(new, create=create)
                 elif isinstance(old, ss.ndict): # Update module containers
                     self._update_ndict(key, old, new)
                 elif isinstance(old, ss.Module):  # Update modules
                     self._update_module(key, old, new)
+                elif isinstance(old, ss.TimePar):
+                    self._update_timepar(key, old, new)
                 elif isinstance(old, ss.Dist): # Update a distribution
                     self._update_dist(key, old, new)
                 elif callable(old): # It's a function: update directly
@@ -71,7 +74,7 @@ def update(self, pars=None, create=False, **kwargs):
                     ss.warn(warnmsg)
                     self[key] = new
         return self
-    
+        
     def _update_ndict(self, key, old, new):
         """ Update an ndict object in the parameters, e.g. sim.pars.diseases """
         if not len(old): # It's empty, just overwrite
@@ -94,8 +97,38 @@ def _update_module(self, key, old, new):
             raise TypeError(errormsg)
         return
     
+    def _update_timepar(self, key, old, new):
+        """ Update a time parameter (duration or rate) """
+        
+        # It's a TimePar, e.g. dur_inf = ss.dur(6); use directly
+        if isinstance(new, ss.TimePar): 
+            new_cls = new.__class__
+            old_cls = old.__class__
+            if new_cls == old_cls:
+                self[key] = new
+            else:
+                errormsg = f'When updating a time parameter, the new class ({new_cls}) should match the original class ({old_cls})'
+                raise TypeError(errormsg)
+        
+        # It's a single number, e.g. dur_inf = 6; set parameters
+        elif isinstance(new, Number):
+            old.set(new)
+        
+        # It's a list of numbers, e.g. dur_inf = [6, 2]; set parameters
+        elif isinstance(new, list):
+            old.set(*new)
+        
+        # It's a dict, figure out what to do
+        elif isinstance(new, dict):
+            if isinstance(old, ss.beta):
+                self[key] = new # TODO: use an actual set here
+            else:
+                old.set(**new)
+        return
+    
     def _update_dist(self, key, old, new):
         """ Update a Dist object in the parameters, e.g. sim.pars.diseases.sir.dur_inf """
+        
         # It's a Dist, e.g. dur_inf = ss.normal(6,2); use directly
         if isinstance(new, ss.Dist): 
             if isinstance(old, ss.bernoulli) and not isinstance(new, ss.bernoulli):
@@ -108,7 +141,7 @@ def _update_dist(self, key, old, new):
         elif isinstance(new, Number):
             old.set(new)
         
-        # It's a list number, e.g. dur_inf = [6, 2]; set parameters
+        # It's a list of numbers, e.g. dur_inf = [6, 2]; set parameters
         elif isinstance(new, list):
             old.set(*new)
         
@@ -168,12 +201,12 @@ def __init__(self, **kwargs):
         self.pop_scale = None  # How much to scale the population
 
         # Simulation parameters
-        self.unit       = 'year' # The time unit to use (NOT YET IMPLEMENTED)
+        self.unit       = 'year' # The time unit to use; options are 'year' (default), 'day', and 'none'
         self.start      = 2000   # Start of the simulation
-        self.end        = None   # End of the simulation
-        self.n_years    = 50     # Number of years to run, if end isn't specified. Note that this includes burn-in
-        self.dt         = 1.0    # Timestep
-        self.rand_seed  = 1      # Random seed, if None, don't reset
+        self.stop       = None   # End of the simulation
+        self.dur        = 50     # Duration of time to run, if stop isn't specified
+        self.dt         = 1.0    # Timestep (in units of self.unit)
+        self.rand_seed  = 1      # Random seed; if None, don't reset
         self.slot_scale = 5      # Random slots will be assigned to newborn agents between min=n_agents and max=slot_scale*n_agents
 
         # Demographic parameters
@@ -213,8 +246,7 @@ def validate_sim_pars(self):
         self.validate_verbose()
         self.validate_agents()
         self.validate_total_pop()
-        self.validate_start_end()
-        self.validate_dt()
+        self.validate_time()
         return
     
     def validate_verbose(self):
@@ -258,43 +290,28 @@ def validate_total_pop(self):
             self.pop_scale = total_pop / self.n_agents
         return
     
-    def validate_start_end(self):
-        """ Ensure at least one of n_years and end is defined, but not both """
-        if self.end is not None:
-            if self.is_default('n_years'):
-                self.n_years = self.end - self.start
+    def validate_time(self):
+        """ Ensure at least one of dur and stop is defined, but not both """
+        if isinstance(self.start, str):
+            self.start = sc.date(self.start)
+        if isinstance(self.stop, str):
+            self.stop = sc.date(self.stop)
+        if self.stop is not None:
+            if self.is_default('dur'):
+                self.dur = ss.date_diff(self.start, self.stop, self.unit)
             else:
-                errormsg = f'You can supply either end ({self.end}) or n_years ({self.n_years}) but not both, since one is calculated from the other'
+                errormsg = f'You can supply either stop ({self.stop}) or dur ({self.dur}) but not both, since one is calculated from the other'
                 raise ValueError(errormsg)
-            if self.n_years <= 0:
-                errormsg = f"Number of years must be >0, but you supplied start={str(self.start)} and " \
-                           f"end={str(self.end)}, which gives n_years={self.n_years}"
+            if self.dur <= 0:
+                errormsg = f"Duration must be >0, but you supplied start={str(self.start)} and stop={str(self.stop)}, which gives dur={self.dur}"
                 raise ValueError(errormsg)
         else:
-            if self.n_years is not None:
-                self.end = self.start + self.n_years
+            if self.dur is not None:
+                self.stop = ss.date_add(self.start, self.dur, self.unit)
             else:
-                errormsg = 'You must supply one of n_years and end."'
+                errormsg = 'You must supply either "dur" or "stop".'
                 raise ValueError(errormsg)
         return
-
-    def validate_dt(self):
-        """
-        Check that 1/dt is an integer value, otherwise results and time vectors will have mismatching shapes.
-        init_results explicitly makes this assumption by casting resfrequency = int(1/dt).
-        """
-        if self.unit == 'year': # TODO: implement properly
-            dt = self.dt
-            reciprocal = 1.0 / dt  # Compute the reciprocal of dt
-            if not reciprocal.is_integer():  # Check if reciprocal is not a whole (integer) number
-                # Round the reciprocal
-                reciprocal = int(reciprocal)
-                rounded_dt = 1.0 / reciprocal
-                self.dt = rounded_dt
-                if self.verbose:
-                    warnmsg = f'Warning: Provided time step dt={dt} resulted in a non-integer number of steps per year. Rounded to {rounded_dt}.'
-                    ss.warn(warnmsg)
-        return
     
     def validate_modules(self):
         """ Validate modules passed in pars"""
diff --git a/starsim/people.py b/starsim/people.py
index 06246b46..d97cf737 100644
--- a/starsim/people.py
+++ b/starsim/people.py
@@ -132,9 +132,10 @@ def get_age_dist(age_data):
     def link_sim(self, sim):
         """ Initialization """
         if self.initialized:
-            errormsg = 'Cannot re-initialize a People object directly; use sim.initialize(reset=True)'
+            errormsg = 'Cannot re-initialize a People object directly; use sim.init(reset=True)'
             raise RuntimeError(errormsg)
         self.sim = sim # Store the sim
+        self.timevec = sim.timevec # Also store the timevec in the People object (not used, but for completeness)
         ss.link_dists(obj=self.states, sim=sim, module=self, skip=[ss.Sim, ss.Module])
         return
     
@@ -190,7 +191,7 @@ def _link_state(self, state, die=True):
 
         All states should be registered by this function for the purpose of connecting them to the
         People's UIDs and to have them be automatically resized when the number of agents changes.
-        This operation is normally triggered as part of initializing the state (via `State.initialize()`)
+        This operation is normally triggered as part of initializing the state (via `State.init()`)
         """
         if id(state) not in self._states:
             self._states[id(state)] = state
@@ -277,7 +278,7 @@ def update_post(self):
         """ Final updates at the very end of the timestep """
         sim = self.sim
         if sim.pars.use_aging:
-            self.age[self.alive.uids] += sim.dt
+            self.age[self.alive.uids] += sim.dt_year
         return
 
     def request_death(self, uids):
@@ -308,10 +309,17 @@ def request_death(self, uids):
         self.ti_dead[uids] = self.sim.ti
         return
 
-    def check_deaths(self):
+    def step_die(self):
         """ Carry out any deaths that took place this timestep """
         death_uids = (self.ti_dead <= self.sim.ti).uids
         self.alive[death_uids] = False
+        
+        # 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
+        for disease in self.sim.diseases():
+            if isinstance(disease, ss.Disease):
+                disease.step_die(death_uids)
+                
         return death_uids
     
     def remove_dead(self):
@@ -349,7 +357,7 @@ def update_results(self):
         return
     
     def finish_step(self):
-        self.update_results()
+        # self.update_results() # This is called separately
         self.remove_dead()
         self.update_post()
         return
diff --git a/starsim/products.py b/starsim/products.py
index 0b800139..bffd3bd1 100644
--- a/starsim/products.py
+++ b/starsim/products.py
@@ -104,26 +104,20 @@ def administer(self, uids, return_format='dict'):
         tx_successful = []  # Initialize list of successfully treated individuals
 
         for disease_name in self.diseases:
-
             disease = self.sim.diseases[disease_name]
 
             for state in self.health_states:
-
                 pre_tx_state = getattr(disease, state)
                 true_uids = pre_tx_state.uids # People in this state
                 these_uids = true_uids.intersect(uids)
 
                 if len(these_uids):
-
                     df_filter = (self.df.state == state) & (self.df.disease == disease_name)  # Filter by state
                     thisdf = self.df[df_filter]  # apply filter to get the results for this state & genotype
 
                     # Determine whether treatment is successful
                     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(self.sim.ti+1)
-                    eff_treat_inds = self.efficacy_dist.filter(these_uids)
+                    eff_treat_inds = self.efficacy_dist.filter(these_uids) # TODO: think if there's a way of not calling this inside a loop like this
 
                     post_tx_state_name = thisdf.post_state.values[0]
                     post_tx_state = getattr(disease, post_tx_state_name)
diff --git a/starsim/run.py b/starsim/run.py
index 5305f877..655d1c66 100644
--- a/starsim/run.py
+++ b/starsim/run.py
@@ -64,13 +64,13 @@ def init_sims(self, **kwargs):
 
         return
 
-    def run(self, reduce=False, **kwargs):
+    def run(self, n_runs=4, **kwargs):
         """
         Run the sims
 
         Args:
-            reduce  (bool): whether to reduce after running (see reduce())
-            kwargs  (dict): passed to multi_run(); use run_args to pass arguments to sim.run()
+            n_runs (int): whether to reduce after running (see reduce())
+            kwargs (dict): passed to multi_run(); use run_args to pass arguments to sim.run()
 
         Returns:
             None (modifies MultiSim object in place)
@@ -89,11 +89,7 @@ def run(self, reduce=False, **kwargs):
 
         # Run
         kwargs = sc.mergedicts(self.run_args, kwargs)
-        self.sims = multi_run(sims, **kwargs)
-
-        # Reduce
-        if reduce:
-            self.reduce()
+        self.sims = multi_run(sims, n_runs=n_runs, **kwargs)
 
         return self
 
@@ -172,7 +168,7 @@ def reduce(self, quantiles=None, use_mean=False, bounds=None, output=False):
 
         # Calculate the statistics
         raw = {}
-        exclude = ['yearvec'] # Exclude things that we don't want to include
+        exclude = ['timevec'] # Exclude things that we don't want to include
 
         rflat = reduced_sim.results.flatten()
         rkeys = list(rflat.keys())
@@ -229,6 +225,39 @@ def median(self, quantiles=None, **kwargs):
         """
         return self.reduce(use_mean=False, quantiles=quantiles, **kwargs)
     
+    def summarize(self, method='mean', quantiles=None, how='default'):
+        """
+        Summarize the simulations statistically.
+        
+        Args:
+            method (str): one of 'mean' (default: [mean, 2*std]), 'median' ([median, min, max]), or 'all' (all results)
+            quantiles (dict): if method='median', use these quantiles
+            how (str): passed to sim.summarize()
+        """
+        
+        # Compute the summaries
+        summaries = []
+        for sim in self.sims:
+            summaries.append(sim.summarize(how=how))
+            
+        summary = sc.dcp(summaries[0]) # Use the first one as a template
+        for k in summary.keys():
+            arr = np.array([s[k] for s in summaries])
+            if method == 'all':
+                summary[k] = arr
+            elif method == 'mean':
+                summary[k] = sc.objdict({'mean':arr.mean(), 'std':arr.std(), 'sem':sc.sem(arr)})
+            elif method == 'median':
+                if quantiles is None:
+                    quantiles = sc.objdict({'median':0.5, 'min':0, 'max':1, 'q25':0.25, 'q75':0.75})
+                elif isinstance(quantiles, list):
+                    quantiles = {q:q for q in quantiles}
+                summary[k] = {q:v for q,v in zip(quantiles, np.quantile(arr, quantiles))}
+        
+        self.summary = summary # Could reconcile with reduce()'s summary
+                                  
+        return summary
+    
     def plot(self, key=None, fig=None, fig_kw=None, plot_kw=None, fill_kw=None):
         """ 
         Plot all results in the MultiSim object.
@@ -255,7 +284,7 @@ def plot(self, key=None, fig=None, fig_kw=None, plot_kw=None, fill_kw=None):
         # Has been reduced, plot with uncertainty bounds
         else:
             flat = self.results
-            yearvec = flat.pop('yearvec')
+            timevec = flat.pop('timevec')
             n_cols = np.ceil(np.sqrt(len(flat))) # TODO: remove duplication with sim.plot()
             default_figsize = np.array([8, 6])
             figsize_factor = np.clip((n_cols-3)/6+1, 1, 1.5) # Scale the default figure size based on the number of rows and columns
@@ -274,8 +303,8 @@ def plot(self, key=None, fig=None, fig_kw=None, plot_kw=None, fill_kw=None):
                     
                 # Do the plotting
                 for ax, (key, res) in zip(axs.flatten(), flat.items()):
-                    ax.fill_between(yearvec, res.low, res.high, **fill_kw)
-                    ax.plot(yearvec, res, **plot_kw)
+                    ax.fill_between(timevec, res.low, res.high, **fill_kw)
+                    ax.plot(timevec, res, **plot_kw)
                     ax.set_title(getattr(res, 'label', key)) 
                     ax.set_xlabel('Year')
                 
@@ -354,11 +383,12 @@ def single_run(sim, ind=0, reseed=True, keep_people=False, run_args=None, sim_ar
 def multi_run(sim, n_runs=4, reseed=None, iterpars=None, keep_people=None, run_args=None, sim_args=None,
               par_args=None, do_run=True, parallel=True, n_cpus=None, verbose=None, **kwargs):
     """
-    For running multiple runs in parallel. If the first argument is a list of sims,
-    exactly these will be run and most other arguments will be ignored.
+    For running multiple sims in parallel. If the first argument is a list of sims
+    rather than a single sim, exactly these will be run and most other arguments 
+    will be ignored.
 
     Args:
-        sim         (Sim)   : the sim instance to be run, or a list of sims.
+        sim         (Sim/list): the sim instance to be run, or a list of sims.
         n_runs      (int)   : the number of parallel runs
         reseed      (bool)  : whether or not to generate a fresh seed for each run (default: true for single, false for list of sims)
         iterpars    (dict)  : any other parameters to iterate over the runs; see sc.parallelize() for syntax
diff --git a/starsim/settings.py b/starsim/settings.py
index 1bccd52e..3c7ac87d 100644
--- a/starsim/settings.py
+++ b/starsim/settings.py
@@ -68,11 +68,14 @@ def get_orig_options():
 
         optdesc.warnings = 'How warnings are handled: options are "warn" (default), "print", and "error"'
         options.warnings = str(os.getenv('STARSIM_WARNINGS', 'warn'))
+        
+        optdesc.time_eps = 'Set size of smallest possible time unit (in units of sim time, e.g. "year" or "day")'
+        options.time_eps = float(os.getenv('STARSIM_TIME_EPS', 1e-6))
 
         optdesc.sep = 'Set thousands seperator for text output'
         options.sep = str(os.getenv('STARSIM_SEP', ','))
 
-        optdesc.precision = 'Set arithmetic precision -- 32-bit by default for efficiency'
+        optdesc.precision = 'Set arithmetic precision'
         options.precision = int(os.getenv('STARSIM_PRECISION', 64))
 
         optdesc._centralized = 'If True, revert to centralized random number generation (NOT ADVISED).'
@@ -174,7 +177,7 @@ def context(self, **kwargs):
 
             # Convert warnings to errors
             with ss.options.context(warnings='error'):
-                ss.Sim(location='not a location').initialize()
+                ss.Sim(location='not a location').init()
 
             # Use with_style(), not context(), for plotting options
             with ss.options.with_style(dpi=50):
diff --git a/starsim/sim.py b/starsim/sim.py
index b3026a90..f7e21e5e 100644
--- a/starsim/sim.py
+++ b/starsim/sim.py
@@ -28,16 +28,17 @@ def __init__(self, pars=None, label=None, people=None, demographics=None, diseas
         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; see also sim.initialize() for more
+        # Set attributes; see also sim.init() 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.loop = ss.Loop(self) # Initialize the integration loop
+        self.results = ss.Results(module='sim')  # For storing results
         self.initialized = False  # Whether initialization is complete
         self.complete = False  # Whether a simulation has completed running
         self.results_ready = False  # Whether results are ready
-        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
@@ -56,7 +57,10 @@ def __repr__(self):
         """ Show a quick version of the sim """
         # Try a more custom repr first
         try:
+            labelstr = f'{self.label}; ' if self.label else ''
             n = int(self.pars.n_agents)
+            timestr = f'{self.pars.start}—{self.pars.stop}'
+            
             moddict = {}
             for modkey in ss.module_map().keys():
                 if hasattr(self, modkey):
@@ -65,8 +69,12 @@ def __repr__(self):
                     thismodtype = self.pars[modkey]
                 else:
                     thismodtype = {}
-                if sc.isiterable(thismodtype) and len(thismodtype):
+                    
+                if isinstance(thismodtype, dict) and len(thismodtype):
                     moddict[modkey] = sc.strjoin(thismodtype.keys())
+                elif isinstance(thismodtype, str):
+                    moddict[modkey] = thismodtype
+                
             if len(moddict):
                 modulestr = ''
                 for k,mstr in moddict.items():
@@ -75,7 +83,8 @@ def __repr__(self):
                 modulestr = ''
             if not self.initialized:
                 modulestr += '; not initialized'
-            string = f'Sim(n={n:n}{modulestr})'
+                
+            string = f'Sim({labelstr}n={n:n}; {timestr}{modulestr})'
         
         # Or just use default
         except Exception as E:
@@ -84,22 +93,37 @@ def __repr__(self):
             
         return string
     
-    def initialize(self, **kwargs):
+    @property
+    def now(self):
+        """ Return the current time, i.e. the time vector at the current timestep """
+        try:
+            ti = min(self.ti, len(self.timevec)-1) # During integration, ti can go one past the end of the time vector
+            return self.timevec[ti]
+        except Exception as E:
+            ss.warn(f'Encountered exception in sim when getting the current time: {E}')
+            return None
+        
+    @property
+    def modules(self):
+        """ Return iterator over all Module instances (stored in standard places) in the Sim """
+        return itertools.chain(
+            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(),
+        )
+    
+    def init(self, **kwargs):
         """ Perform all initializations for the sim; most heavy lifting is done by the parameters """
+        
         # Validation and initialization
         ss.set_seed(self.pars.rand_seed) # Reset the seed before the population is created -- shouldn't matter if only using Dist objects
-        
-        # Validate parameters
-        self.pars.validate()
-
-        # Initialize time
-        self.init_time_attrs()
-        
-        # Initialize the people
-        self.init_people(**kwargs)  # Create all the people
-        
-        # # Initialize the modules within the parameters
-        # self.pars.validate_modules(self)
+        self.pars.validate() # Validate parameters
+        self.init_time_attrs() # Initialize time
+        self.init_people(**kwargs) # Initialize the people
         
         # Move initialized modules to the sim
         keys = ['label', 'demographics', 'networks', 'diseases', 'interventions', 'analyzers', 'connectors']
@@ -115,29 +139,29 @@ def initialize(self, **kwargs):
             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)
-        
-        # Initialize the values in all of the states and networks
-        self.init_vals()
+        # Final initializations
+        self.dists.init(obj=self, base_seed=self.pars.rand_seed, force=True) # Initialize all distributions now that everything else is in place
+        self.init_vals() # Initialize the values in all of the states and networks
+        self.init_results() # Initialize the results
+        self.loop.init() # Initialize the integration loop
+        self.timer = sc.timer(start=False) # Store a timer for keeping track of how long the run takes
+        self.verbose = self.pars.verbose # Store a run-specific value of verbose
         
-        # Initialize the results
-        self.init_results()
-
         # It's initialized
         self.initialized = True
         return self
     
     def init_time_attrs(self):
         """ Time indexing; derived values live in the sim rather than in the pars """
-        self.dt = self.pars.dt # Shortcut to dt since used a lot
-        self.yearvec = np.arange(start=self.pars.start, stop=self.pars.end + self.pars.dt, step=self.pars.dt) # The time points of the sim
-        self.results.yearvec = self.yearvec # Store the yearvec in the results for plotting
-        self.npts = len(self.yearvec) # The number of points in the sim
-        self.tivec = np.arange(self.npts) # The vector of time indices
+        pars = self.pars
+        self.timevec = ss.make_timevec(pars.start, pars.stop, pars.dt, pars.unit)
+        self.results.timevec = self.timevec # Store the timevec in the results for plotting # TODO: instead, store a timevec with each result
+        self.npts = len(self.timevec) # The number of points in the sim
+        self.timearray = np.arange(self.npts)*pars.dt # Absolute time array
         self.ti = 0  # The time index, e.g. 0, 1, 2
+        self.dt_year = ss.time_ratio(pars.unit, pars.dt, 'year', 1.0) # Figure out what dt is in years; used for demographics # TODO: handle None
         return
-
+    
     def init_people(self, verbose=None, **kwargs):
         """
         Initialize people within the sim
@@ -171,11 +195,17 @@ def init_vals(self):
         # Initialize values in people
         self.people.init_vals()
         
-        # Initialize values in other modules, including networks
+        # Initialize values in other modules, including networks and time parameters
         for mod in self.modules:
             mod.init_post()
         return
     
+    def reset_time_pars(self, force=True):
+        """ Reset the time parameters in the modules; used for imposing the sim's timestep on the modules """
+        for mod in self.modules:
+            mod.init_time_pars(force=force)
+        return
+    
     def init_results(self):
         """ Create initial results that are present in all simulations """
         self.results += [ # TODO: refactor with self.add_results()
@@ -185,133 +215,78 @@ def init_results(self):
         ]
         return
 
-    @property
-    def modules(self):
-        """ Return iterator over all Module instances (stored in standard places) in the Sim """
-        return itertools.chain(
-            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): # TODO: remove when we do the time refactor
-        return self.yearvec[self.ti]
-
-    def step(self):
-        """
-        Step through time and update values
-        
-        TODO: assemble every method called during integration in a list, and then call them sequentially        
-        """
+    def start_step(self):
+        """ Step through time and update values """
 
         # Set the time and if we have reached the end of the simulation, then do nothing
         if self.complete:
-            errormsg = 'Simulation already complete (call sim.initialize() to re-run)'
+            errormsg = 'Simulation already complete (call sim.init() 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 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():
-            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.step()
-
-        # Update networks - this takes place here in case autonomous state changes at this timestep
-        # affect eligibility for contacts
-        for network in self.networks():
-            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 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.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.check_deaths()
-        for disease in self.diseases():
-            if isinstance(disease, ss.Disease):
-                disease.step_die(uids)
-
-        # Update results
-        self.people.update_results()
-        for mod in self.modules:
-            mod.update_results()
+        # Print out progress if needed
+        self.elapsed = self.timer.toc(output=True)
+        if self.verbose: # Print progress
+            simlabel = f'"{self.label}": ' if self.label else ''
+            timepoint = self.timevec[self.ti]
+            timelabel = f'{timepoint:0.1f}' if isinstance(timepoint, float) else str(timepoint) # TODO: fix
+            string = f'  Running {simlabel}{timelabel} ({self.ti:2.0f}/{self.npts}) ({self.elapsed:0.2f} s) '
+            if self.verbose >= 2:
+                sc.heading(string)
+            elif self.verbose > 0:
+                if not (self.ti % int(1.0 / self.verbose)):
+                    sc.progressbar(self.ti + 1, self.npts, label=string, length=20, newline=True)
 
-        # Apply analyzers
-        for ana in self.analyzers():
-            ana.step()
-            
-        # Clean up dead agents and perform other housekeeping tasks
-        self.people.finish_step()
+        # Advance random number generators forward to prepare for any random number calls that may be necessary on this step
+        self.dists.jump_dt() # TODO: each module should do this, but should be ok as-is with auto-jump
+        return
 
-        # Tidy up
+    def finish_step(self):
+        """ Finish the simulation timestep """
         self.ti += 1
-        if self.ti == self.npts:
-            self.complete = True
         return
+    
+    def run_one_step(self, verbose=None):
+        """
+        Run a single sim step; only used for debugging purposes.
+        
+        Note: sim.run_one_step() runs a single simulation timestep, which involves
+        multiple function calls. In contrast, loop.run_one_step() runs a single
+        function call.
+        
+        Note: the verbose here is only for the Loop object, not the sim.        
+        """
+        self.loop.run(self.now, verbose)
+        return self
 
     def run(self, until=None, verbose=None):
         """ Run the model once """
 
         # Initialization steps
-        T = sc.timer()
-        if not self.initialized:
-            self.initialize()
-        verbose = sc.ifelse(verbose, self.pars.verbose)
+        if not self.initialized: self.init()
+        self.verbose = sc.ifelse(verbose, self.pars.verbose)
+        self.timer.start()
 
         # Check for AlreadyRun errors
         errormsg = None
-        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
-            errormsg = f'Simulation is currently at t={self.ti}, requested to run until ti={until} which has already been reached'
         if self.complete:
-            errormsg = 'Simulation is already complete (call sim.initialize() to re-run)'
+            errormsg = 'Simulation is already complete (call sim.init() to re-run)'
         if errormsg:
             raise AlreadyRunError(errormsg)
 
-        # Main simulation loop
-        while self.ti < until:
-            elapsed = T.toc(output=True)
-            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:
-                    sc.heading(string)
-                elif verbose > 0:
-                    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 -- lots happens on the next line!
-            self.step()
+        # Main simulation loop -- just one line!!!
+        self.loop.run(until)
+        
+        # Check if the simulation is complete
+        if self.loop.index == len(self.loop.plan):
+            self.complete = True
 
         # If simulation reached the end, finalize the results
         if self.complete:
             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
+            for mod in self.modules: # May not be needed, but keeps it consistent with the sim
+                mod.ti -= 1
             self.finalize()
-            sc.printv(f'Run finished after {elapsed:0.2f} s.\n', 1, verbose)
+            sc.printv(f'Run finished after {self.elapsed:0.2f} s.\n', 1, self.verbose)
         return self # Allows e.g. ss.Sim().run().plot()
 
     def finalize(self):
@@ -428,7 +403,7 @@ def export_df(self):
             errormsg = 'Please run the sim before exporting the results'
             raise RuntimeError(errormsg)
 
-        def flatten_results(d, prefix=''):
+        def flatten_results(d, prefix=''): # TODO: won't work with different timesteps
             flat = {}
             for key, val in d.items():
                 if isinstance(val, dict):
@@ -438,7 +413,7 @@ def flatten_results(d, prefix=''):
             return flat
 
         resdict = flatten_results(self.results)
-        resdict['t'] = self.yearvec
+        resdict['t'] = self.timevec
         df = sc.dataframe.from_dict(resdict).set_index('t')
         return df
 
@@ -614,7 +589,7 @@ def plot(self, key=None, fig=None, style='fancy', fig_kw=None, plot_kw=None):
         # Do the plotting
         with sc.options.with_style(style):
             
-            yearvec = flat.pop('yearvec')
+            timevec = flat.pop('timevec')
             if key is not None:
                 flat = {k:v for k,v in flat.items() if k.startswith(key)}
             
@@ -630,7 +605,7 @@ def plot(self, key=None, fig=None, style='fancy', fig_kw=None, plot_kw=None):
                 
             # Do the plotting
             for ax, (key, res) in zip(axs, flat.items()):
-                ax.plot(yearvec, res, **plot_kw, label=self.label)
+                ax.plot(timevec, res, **plot_kw, label=self.label)
                 title = getattr(res, 'label', key)
                 if res.module != 'sim':
                     try:
@@ -691,8 +666,8 @@ def diff_sims(sim1, sim2, skip_key_diffs=False, skip=None, full=False, output=Fa
     values which differ.
 
     Args:
-        sim1 (sim/dict): either a simulation object or the sim.summary dictionary
-        sim2 (sim/dict): ditto
+        sim1 (Sim/MultiSim/dict): either a simulation/MultiSim object or the sim.summary dictionary
+        sim2 (im/dict): ditto
         skip_key_diffs (bool): whether to skip keys that don't match between sims
         skip (list): a list of values to skip
         full (bool): whether to print out all values (not just those that differ)
@@ -706,15 +681,25 @@ def diff_sims(sim1, sim2, skip_key_diffs=False, skip=None, full=False, output=Fa
         s2 = ss.Sim(rand_seed=2).run()
         ss.diff_sims(s1, s2)
     '''
-
-    if isinstance(sim1, Sim):
+    
+    # Convert to dict
+    if isinstance(sim1, (Sim, ss.MultiSim)):
         sim1 = sim1.summarize()
-    if isinstance(sim2, Sim):
+    if isinstance(sim2, (Sim, ss.MultiSim)):
         sim2 = sim2.summarize()
     for sim in [sim1, sim2]:
         if not isinstance(sim, dict):  # pragma: no cover
             errormsg = f'Cannot compare object of type {type(sim)}, must be a sim or a sim.summary dict'
             raise TypeError(errormsg)
+    
+    
+    # Check if it's a multisim
+    sim1 = sc.objdict(sim1)
+    sim2 = sc.objdict(sim2)
+    multi = isinstance(sim1[0], dict) # If it's a dict, then it's a multisim
+    if multi:
+        for sim in [sim1, sim2]:
+            assert 'mean' in sim[0], f"Can only compare multisims with summarize(method='mean'), but you have keys {sim[0].keys()}"
 
     # Compare keys
     keymatchmsg = ''
@@ -725,7 +710,7 @@ def diff_sims(sim1, sim2, skip_key_diffs=False, skip=None, full=False, output=Fa
         missing = list(sim1_keys - sim2_keys)
         extra = list(sim2_keys - sim1_keys)
         if missing:
-            keymatchmsg += f'  Missing sim1 keys: {missing}\ns'
+            keymatchmsg += f'  Missing sim1 keys: {missing}\n'
         if extra:
             keymatchmsg += f'  Extra sim2 keys: {extra}\n'
 
@@ -736,29 +721,49 @@ def diff_sims(sim1, sim2, skip_key_diffs=False, skip=None, full=False, output=Fa
     skip = sc.tolist(skip)
     for key in sim2.keys():  # To ensure order
         if key in sim1_keys and key not in skip:  # If a key is missing, don't count it as a mismatch
-            sim1_val = sim1[key] if key in sim1 else 'not present'
-            sim2_val = sim2[key] if key in sim2 else 'not present'
-            mm = not np.isclose(sim1_val, sim2_val, equal_nan=True)
+            d = sc.objdict()
+            if multi:
+                d.sim1 = sim1[key]['mean']
+                d.sim2 = sim2[key]['mean']
+                d.sim1_sem = sim1[key]['sem']
+                d.sim2_sem = sim2[key]['sem']
+            else:
+                d.sim1 = sim1[key]
+                d.sim2 = sim2[key]
+            mm = not np.isclose(d.sim1, d.sim2, equal_nan=True)
             n_mismatch += mm
             if mm or full:
-                mismatches[key] = {'sim1': sim1_val, 'sim2': sim2_val}
+                mismatches[key] = d
 
+    df = sc.dataframe() # Preallocate in case there were no mismatches
     if len(mismatches):
         valmatchmsg = '\nThe following values differ between the two simulations:\n' if not full else ''
         df = sc.dataframe.from_dict(mismatches).transpose()
         diff = []
         ratio = []
         change = []
-        small_change = 1e-3  # Define a small change, e.g. a rounding error
-        for mdict in mismatches.values():
-            old = mdict['sim1']
-            new = mdict['sim2']
-            numeric = sc.isnumber(sim1_val) and sc.isnumber(sim2_val)
+        zscore = []
+        small_change = 1 + 1e-3  # Define a small change, e.g. a rounding error
+        for d in mismatches.values():
+            old = d.sim1
+            new = d.sim2
+            if multi:
+                old_sem = d.sim1_sem
+                new_sem = d.sim1_sem
+                sem = old_sem + new_sem
+                small_change = 1.96 # 95% CI, roughly speaking
+                
+            numeric = sc.isnumber(old) and sc.isnumber(new) # Should all be numeric, but just in case
             if numeric and old > 0:
                 this_diff = new - old
                 this_ratio = new / old
-                abs_ratio = max(this_ratio, 1.0 / this_ratio)
-
+                if multi:
+                    abs_ratio = abs(this_diff)/sem
+                    this_zscore = abs_ratio
+                else:
+                    abs_ratio = max(this_ratio, 1.0/this_ratio)
+                    this_zscore = None
+                
                 # Set the character to use
                 approx_eq = abs_ratio < small_change
                 if approx_eq:    change_char = '≈'
@@ -780,18 +785,26 @@ def diff_sims(sim1, sim2, skip_key_diffs=False, skip=None, full=False, output=Fa
                 this_diff =  np.nan
                 this_ratio  = np.nan
                 this_change = 'N/A'
+                this_zscore = np.nan
 
             diff.append(this_diff)
             ratio.append(this_ratio)
             change.append(this_change)
+            zscore.append(this_zscore)
 
         df['diff'] = diff
         df['ratio'] = ratio
-        for col in ['sim1', 'sim2', 'diff', 'ratio']:
+        numeric_cols = ['sim1', 'sim2', 'diff', 'ratio']
+        if multi:
+            numeric_cols += ['sim1_sem', 'sim2_sem']
+        for col in numeric_cols:
             df[col] = df[col].round(decimals=3)
         df['change'] = change
+        if multi:
+            df['zscore'] = zscore
+            df['statsig'] = df.zscore > small_change
         valmatchmsg += str(df)
-
+        
     # Raise an error if mismatches were found
     mismatchmsg = keymatchmsg + valmatchmsg
     if mismatchmsg:  # pragma: no cover
diff --git a/starsim/states.py b/starsim/states.py
index 5a72eee5..ffbaa063 100644
--- a/starsim/states.py
+++ b/starsim/states.py
@@ -170,25 +170,23 @@ def __or__(self, other):  raise BooleanOperationError(self)
     def __xor__(self, other): raise BooleanOperationError(self)
     def __invert__(self):     raise BooleanOperationError(self)
 
-    def __array_ufunc__(self, *args, **kwargs):
-        if args[1] != '__call__':
+    def __array_ufunc__(self, ufunc, method, *args, out=None, **kwargs):
+        args = [(arg.values if arg is self else arg) for arg in args]
+        if method != '__call__':
             # This is a catch-all for ufuncs that are not being applied with '__call__' (e.g., operations returning a scalar like 'np.sum()' use reduce instead)
-            args = [(x if x is not self else self.values) for x in args]
-            kwargs = {k: v if v is not self else self.values for k, v in kwargs.items()}
-            return self.values.__array_ufunc__(*args, **kwargs)
+            kwargs = {k:v.values if v is self else v for k,v in kwargs.items()}
+            return self.values.__array_ufunc__(ufunc, method, *args, out=out, **kwargs)
         else:
-            args = [(x if x is not self else self.values) for x in args] # Convert any operands that are Arr instances to their value arrays
-            if 'out' in kwargs and kwargs['out'][0] is self:
+            if out is self:
                 # In-place operations like += applied to the entire Arr instance
                 # use this branch. Therefore, we perform our computation on a new
                 # array with the same size as self.values, and then write it back
                 # to the appropriate entries in `self.raw` via `self[:]`
-                del kwargs['out']
-                self[:] = args[0](*args[2:], **kwargs)
+                self[:] = ufunc(*args, **kwargs)
                 return self
             else:
                 # Otherwise, just run the ufunc
-                return args[0](*args[2:], **kwargs)
+                return ufunc(*args, **kwargs)
 
     @property
     def auids(self):
diff --git a/starsim/timepars.py b/starsim/timepars.py
new file mode 100644
index 00000000..fd964bf9
--- /dev/null
+++ b/starsim/timepars.py
@@ -0,0 +1,297 @@
+"""
+Functions and classes for handling time
+"""
+
+import numpy as np
+import sciris as sc
+import starsim as ss
+
+# Classes that are externally visible
+__all__ = ['time_units', 'time_ratio', 'date_add', 'date_diff', 'make_timevec', 'make_timearray',
+           'TimePar', 'dur', 'days', 'years', 'rate', 'time_prob', 'beta']
+
+    
+#%% Helper functions
+
+# Define available time units
+time_units = sc.dictobj(
+    day = 1,
+    week = 7,
+    month = 30.4375, # 365.25/12 -- more accurate and nicer fraction
+    year = 365, # For simplicity with days
+)
+
+def time_ratio(unit1='day', dt1=1.0, unit2='day', dt2=1.0, as_int=False):
+    """
+    Calculate the relationship between two sets of time factors
+    
+    Args:
+        unit1 (str): units for the numerator
+        dt1 (float): timestep for the numerator
+        unit2 (str): units for the denominator
+        dt2 (float): timestep for the denominator
+        as_int (bool): round and convert to an integer
+    
+    **Example**::
+        ss.time_ratio(unit1='week', dt1=2, unit2='day', dt2=1, as_int=True) # Returns 14
+    """
+    if dt1 == dt2:
+        dt_ratio = 1.0
+    else:
+        if dt1 is None or dt2 is None:
+            errormsg = f'Cannot convert between dt when one is None ({dt1}, {dt2})'
+            raise ValueError(errormsg)
+        dt_ratio = dt1/dt2
+    
+    if unit1 == unit2:
+        unit_ratio = 1.0
+    else:
+        if unit1 is None or unit2 is None:
+            errormsg = f'Cannot convert between units when one is None ({unit1}, {unit2})'
+            raise ValueError(errormsg)
+        u1 = time_units[unit1]
+        u2 = time_units[unit2]
+        unit_ratio = u1/u2
+        
+    factor = dt_ratio * unit_ratio
+    if as_int:
+        factor = int(round(factor))
+    return factor
+
+
+def date_add(start, dur, unit):
+    """ Add two dates (or integers) together """
+    if sc.isnumber(start):
+        stop = start + dur
+    else:
+        if unit in time_units:
+            ndays = int(round(time_units[unit]*dur))
+            stop = sc.datedelta(start, days=ndays)
+        else:
+            errormsg = f'Unknown unit {unit}, choices are: {sc.strjoin(time_units.keys())}'
+            raise ValueError(errormsg)
+    return stop
+        
+
+def date_diff(start, stop, unit):
+    """ Find the difference between two dates (or integers) """
+    if sc.isnumber(start) and sc.isnumber(stop):
+        dur = stop - start
+    else:
+        if unit == 'year':
+            dur = sc.datetoyear(stop) - sc.datetoyear(start) # TODO: allow non-integer amounts
+        elif unit == 'day':
+            dur = (stop - start).days
+        else:
+            raise NotImplementedError
+    return dur
+
+
+def make_timevec(start, stop, dt, unit):
+    """ Parse start, stop, and dt into an appropriate time vector """
+    if sc.isnumber(start):
+        try:
+            timevec = sc.inclusiverange(start=start, stop=stop, step=dt) # The time points of the sim
+        except Exception as E:
+            errormsg = f'Incompatible set of time inputs: start={start}, stop={stop}, dt={dt}. You can use dates or numbers but not both.'
+            raise ValueError(errormsg) from E
+    else:
+        day_delta = time_ratio(unit1=unit, dt1=dt, unit2='day', dt2=1.0, as_int=True)
+        if day_delta > 0:
+            timevec = sc.daterange(start, stop, interval={'days':day_delta})
+        else:
+            errormsg = f'Timestep {dt} is too small; must be at least 1 day'
+            raise ValueError(errormsg)
+    return timevec
+
+
+def make_timearray(tv, unit, sim_unit):
+    """ Convert a module time vector into a numerical time array with the same units as the sim """
+    
+    # It's an array of days or years: easy
+    if sc.isarray(tv):
+        ratio = time_ratio(unit1=unit, unit2=sim_unit)
+        abstv = tv*ratio # Get the units right
+        abstv -= abstv[0] # Start at 0
+    
+    # It's a date: convert to fractional years and then subtract the 
+    else:
+        yearvec = [sc.datetoyear(d) for d in tv]
+        absyearvec = np.array(yearvec) - yearvec[0] # Subtract start date
+        abstv = absyearvec*time_ratio(unit1='year', unit2=sim_unit)
+    
+    # Round to the value of epsilon; alternative to np.round(abstv/eps)*eps, which has floating point error
+    decimals = int(-np.log10(ss.options.time_eps))
+    abstv = np.round(abstv, decimals=decimals)
+        
+    return abstv
+
+
+#%% Time classes
+
+class TimePar:
+    """
+    Base class for time-aware parameters, durations and rates
+    
+    NB, because the factor needs to be recalculated, do not set values directly.
+    """
+    def __new__(cls, value=None, *args, **kwargs):
+        """ Allow TimePars to wrap distributions and return the distributions """
+        
+        # Special distribution handling
+        if isinstance(value, ss.Dist):
+            dist = value
+            dist.pars[0] = cls(dist.pars[0], *args, **kwargs) # Convert the first parameter to a TimePar (the same scale is applied to all parameters)
+            return dist
+        
+        # Otherwise, do the usual initialization
+        else:
+            return super().__new__(cls)
+    
+    def __init__(self, value, unit=None, parent_unit=None, parent_dt=None):
+        self.value = value
+        self.unit = unit
+        self.parent_unit = parent_unit
+        self.parent_dt = parent_dt
+        self.factor = None
+        self.x = np.nan
+        self.parent_name = None
+        self.initialized = False
+        return
+    
+    def init(self, parent=None, parent_unit=None, parent_dt=None):
+        """ Link to the sim and/or module units """
+        if parent is None:
+            parent = sc.dictobj(unit=parent_unit, dt=parent_dt)
+        else:
+            if parent_dt is not None:
+                errormsg = f'Cannot override parent {parent} by setting parent_dt; set in parent object instead'
+                raise ValueError(errormsg)
+            self.parent_name = parent.__class__.__name__ # TODO: or parent.name for Starsim objects?
+        
+        if parent.unit is not None:
+            self.parent_unit = parent.unit
+            
+        if parent.dt is not None:
+           self.parent_dt = parent.dt
+        
+        # Set defaults if not yet set
+        self.unit = sc.ifelse(self.unit, self.parent_unit)
+        self.parent_dt = sc.ifelse(self.parent_dt, 1.0)
+        
+        # Calculate the actual conversion factor to be used in the calculations
+        self.set_factor()
+        self.set_x()
+        self.initialized = True
+        return self
+        
+    def __repr__(self):
+        name = self.__class__.__name__
+        if self.initialized:
+            if self.factor == 1.0:
+                xstr = ''
+            else:
+                xstr = f', x={self.x}'
+        else:
+            xstr = ', initialized=False'
+        return f'ss.{name}({self.value}, unit={self.unit}{xstr})'
+    
+    def disp(self):
+        return sc.pr(self)
+    
+    def set(self, value=None, unit=None, parent_unit=None, parent_dt=None):
+        """ Reset the parameter values """
+        if value       is not None: self.value       = value
+        if unit        is not None: self.unit        = unit
+        if parent_unit is not None: self.parent_unit = parent_unit
+        if parent_dt   is not None: self.parent_dt   = parent_dt
+        if self.initialized: # Don't try to set these unless it's been initialized
+            self.set_factor()
+            self.set_x()
+        return self
+    
+    def set_factor(self):
+        """ Set factor used to multiply the value to get the output """
+        self.factor = time_ratio(unit1=self.unit, dt1=1.0, unit2=self.parent_unit, dt2=self.parent_dt)
+        return
+    
+    # Act like a float
+    def __add__(self, other): return self.x + other
+    def __sub__(self, other): return self.x - other
+    def __mul__(self, other): return self.x * other
+    def __pow__(self, other): return self.x ** other
+    def __truediv__(self, other): return self.x / other
+    
+    # ...from either side
+    def __radd__(self, other): return other + self.x
+    def __rsub__(self, other): return other - self.x
+    def __rmul__(self, other): return other * self.x
+    def __rpow__(self, other): return other ** self.x
+    def __rtruediv__(self, other): return other / self.x
+    
+    # Handle modify-in-place methods
+    def __iadd__(self, other): self.value += other; return self
+    def __isub__(self, other): self.value -= other; return self
+    def __imul__(self, other): self.value *= other; return self
+    def __itruediv__(self, other): self.value /= other; return self
+    
+    # Other methods
+    def __neg__(self): return -self.x
+    
+    # Unfortunately, floats don't define the above methods, so we can't *just* use this
+    def __getattr__(self, attr):
+        """ Make it behave like a regular float mostly """
+        if attr in ['__deepcopy__', '__getstate__', '__setstate__']:
+            return self.__getattribute__(attr)
+        else:
+            return getattr(self.x, attr)
+        
+    # Similar to Arr -- required for doing efficient array operations
+    def __array_ufunc__(self, ufunc, method, *args, **kwargs):
+        args = [(arg.x if arg is self else arg) for arg in args] # Use self.x for the operation
+        if method == '__call__':
+            return ufunc(*args, **kwargs)
+        else:
+            return self.x.__array_ufunc__(ufunc, method, *args, **kwargs) # Probably not needed
+
+
+class dur(TimePar):
+    """ Any number that acts like a duration """
+    def set_x(self):
+        self.x = self.value*self.factor
+        return
+
+
+class days(dur):
+    """ Shortcut to ss.dur(value, units='day') """
+    def __init__(self, value, parent_unit=None, parent_dt=None):
+        super().__init__(value=value, unit='day', parent_unit=parent_unit, parent_dt=parent_dt)
+        return
+
+
+class years(dur):
+    """ Shortcut to ss.dur(value, units='year') """
+    def __init__(self, value, parent_unit=None, parent_dt=None):
+        super().__init__(value=value, unit='year', parent_unit=parent_unit, parent_dt=parent_dt)
+        return
+
+
+class rate(TimePar): # TODO: should all rates just be time_prob?
+    """ Any number that acts like a rate; can be greater than 1 """
+    def set_x(self):
+        self.x = self.value/self.factor
+        return
+
+
+class time_prob(TimePar):
+    """ A probability over time (a.k.a. a cumulative hazard rate); must be >0 and <1 """
+    def set_x(self):
+        rate = -np.log(1 - self.value)
+        self.x = 1 - np.exp(-rate/self.factor)
+        return
+        
+    
+class beta(time_prob):
+    """ A container for beta (i.e. the disease transmission rate) """
+    pass
+    
\ No newline at end of file
diff --git a/starsim/version.py b/starsim/version.py
index 7f3908e0..c9fc55b3 100644
--- a/starsim/version.py
+++ b/starsim/version.py
@@ -4,6 +4,6 @@
 
 __all__ = ['__version__', '__versiondate__', '__license__']
 
-__version__ = '1.1.0'
-__versiondate__ = '2024-09-10'
+__version__ = '1.2.0'
+__versiondate__ = '2024-09-24'
 __license__ = f'Starsim {__version__} ({__versiondate__}) — © 2023-2024 by IDM'


Test patch:
diff --git a/tests/baseline.json b/tests/baseline.json
index 8005c3d2..a63ed09f 100644
--- a/tests/baseline.json
+++ b/tests/baseline.json
@@ -1,26 +1,26 @@
 {
   "summary": {
-    "yearvec": 2010.1000000000026,
-    "births_new": 65.24509803921569,
-    "births_cumulative": 3261.872549019608,
+    "timevec": 2010.0,
+    "births_new": 65.10891089108911,
+    "births_cumulative": 3224.9504950495048,
     "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
+    "deaths_new": 44.51485148514851,
+    "deaths_cumulative": 2212.0495049504952,
+    "deaths_cmr": 19.801980198019802,
+    "sir_n_susceptible": 2531.2871287128714,
+    "sir_n_infected": 3618.5049504950493,
+    "sir_n_recovered": 4801.425742574257,
+    "sir_prevalence": 0.33809576643794426,
+    "sir_new_infections": 129.88118811881188,
+    "sir_cum_infections": 13118.0,
+    "sis_n_susceptible": 4807.6435643564355,
+    "sis_n_infected": 6143.574257425743,
+    "sis_prevalence": 0.5480454559430055,
+    "sis_new_infections": 181.21782178217822,
+    "sis_cum_infections": 18303.0,
+    "sis_rel_sus": 0.5275405966379092,
+    "n_alive": 10951.217821782178,
+    "new_deaths": 45.742574257425744,
+    "cum_deaths": 4571.0
   }
 }
\ No newline at end of file
diff --git a/tests/baseline_multisim.json b/tests/baseline_multisim.json
new file mode 100644
index 00000000..dfca24f2
--- /dev/null
+++ b/tests/baseline_multisim.json
@@ -0,0 +1,112 @@
+{
+  "yearvec": {
+    "mean": 2010.1000000000026,
+    "std": 0.0,
+    "sem": 0.0
+  },
+  "births_new": {
+    "mean": 65.42156862745097,
+    "std": 0.1636403441859318,
+    "sem": 0.05174762047214366
+  },
+  "births_cumulative": {
+    "mean": 3268.022549019608,
+    "std": 7.404566545623425,
+    "sem": 2.3415295370455103
+  },
+  "births_cbr": {
+    "mean": 29.722512541939885,
+    "std": 0.006276951230389925,
+    "sem": 0.001984946264982848
+  },
+  "deaths_new": {
+    "mean": 44.380392156862754,
+    "std": 0.31770757109406206,
+    "sem": 0.10046795545371096
+  },
+  "deaths_cumulative": {
+    "mean": 2206.050980392157,
+    "std": 20.506215616069515,
+    "sem": 6.484634753729258
+  },
+  "deaths_cmr": {
+    "mean": 20.155330390777614,
+    "std": 0.14640525779385213,
+    "sem": 0.04629740760526911
+  },
+  "sir_n_susceptible": {
+    "mean": 2440.7539215686274,
+    "std": 28.746250762671917,
+    "sem": 9.090362660039563
+  },
+  "sir_n_infected": {
+    "mean": 3592.6529411764704,
+    "std": 15.832047055140144,
+    "sem": 5.006532871720426
+  },
+  "sir_n_recovered": {
+    "mean": 4971.801960784314,
+    "std": 31.176074975909064,
+    "sem": 9.858740542805167
+  },
+  "sir_prevalence": {
+    "mean": 0.33524628593554295,
+    "std": 0.0006949779189371683,
+    "sem": 0.0002197713147365318
+  },
+  "sir_new_infections": {
+    "mean": 128.78137254901958,
+    "std": 0.6174299338251982,
+    "sem": 0.19524848864546648
+  },
+  "sir_cum_infections": {
+    "mean": 13135.7,
+    "std": 62.97785325017042,
+    "sem": 19.915345841837645
+  },
+  "sis_n_susceptible": {
+    "mean": 4016.530392156862,
+    "std": 103.15173055519776,
+    "sem": 32.61944131424099
+  },
+  "sis_n_infected": {
+    "mean": 6988.678431372548,
+    "std": 112.41894818744964,
+    "sem": 35.54999284327985
+  },
+  "sis_prevalence": {
+    "mean": 0.622135221279589,
+    "std": 0.00994225067508044,
+    "sem": 0.003144015720160086
+  },
+  "sis_new_infections": {
+    "mean": 208.82450980392156,
+    "std": 2.531950580258016,
+    "sem": 0.8006730756600289
+  },
+  "sis_cum_infections": {
+    "mean": 21300.1,
+    "std": 258.2589591863175,
+    "sem": 81.6686537173229
+  },
+  "sis_rel_sus": {
+    "mean": 0.46661987224910195,
+    "std": 0.00995264046990363,
+    "sem": 0.003147301261766397
+  },
+  "n_alive": {
+    "mean": 11005.20882352941,
+    "std": 27.11539023756586,
+    "sem": 8.574639279500229
+  },
+  "new_deaths": {
+    "mean": 45.478431372549025,
+    "std": 0.3597865572558359,
+    "sem": 0.11377449924390214
+  },
+  "cum_deaths": {
+    "mean": 4589.7,
+    "std": 36.32093060481793,
+    "sem": 11.485686744814172
+  }
+}
\ No newline at end of file
diff --git a/tests/benchmark.json b/tests/benchmark.json
index 4b148c87..d28967fa 100644
--- a/tests/benchmark.json
+++ b/tests/benchmark.json
@@ -1,12 +1,12 @@
 {
   "time": {
-    "initialize": 0.014,
-    "run": 0.907
+    "initialize": 0.03,
+    "run": 1.036
   },
   "parameters": {
     "n_agents": 10000,
-    "n_years": 20,
+    "dur": 20,
     "dt": 0.2
   },
-  "cpu_performance": 0.8939678660976156
+  "cpu_performance": 0.9791453304037594
 }
\ No newline at end of file
diff --git a/tests/benchmark_full.py b/tests/benchmark_full.py
index 9cc54afa..d811aa18 100644
--- a/tests/benchmark_full.py
+++ b/tests/benchmark_full.py
@@ -10,9 +10,9 @@
 repeats = 10
 n_agents = 100_000
 pars = sc.objdict(
-    n_years       = 100,        # Number of years to simulate
-    dt            = 0.5,        # Timestep
-    verbose       = 0,          # Don't print details of the run
+    dur     = 100, # Number of years to simulate
+    dt      = 0.5, # Timestep
+    verbose = 0,   # Don't print details of the run
 )
 
 def make_run_sim():
diff --git a/tests/benchmark_simple.py b/tests/benchmark_simple.py
index 54e2a445..f3fcac6b 100644
--- a/tests/benchmark_simple.py
+++ b/tests/benchmark_simple.py
@@ -11,8 +11,8 @@
     hiv = ss.HIV()
     hiv.pars['beta'] = {'simple_sexual': [0.0008, 0.0004], 'maternal': [0.2, 0]}
     
-    sim = ss.Sim(start=1950, end=2050, people=ppl, networks=networks, demographics=ss.Pregnancy(), diseases=hiv)
-    sim.initialize()
+    sim = ss.Sim(start=1950, stop=2050, people=ppl, networks=networks, demographics=ss.Pregnancy(), diseases=hiv)
+    sim.init()
     sim.run()
 
 if do_plot:
diff --git a/tests/multisim_baseline.py b/tests/multisim_baseline.py
new file mode 100755
index 00000000..ede2f48d
--- /dev/null
+++ b/tests/multisim_baseline.py
@@ -0,0 +1,22 @@
+#!/usr/bin/env python3
+
+"""
+Run this script to save or check the multisim baseline.
+
+You can supply command line arguments for "save" and "n_runs" as integers, e.g.
+    ./multisim_baseline.py 0 20
+"""
+
+import sys
+import test_baselines
+
+save = False
+n_runs = 10
+
+if len(sys.argv) > 1:
+    save = int(sys.argv[1])
+    assert save in [0,1] # Avoid confusing with n_runs
+if len(sys.argv) > 2:
+    n_runs = int(sys.argv[2])
+
+test_baselines.multisim_baseline(save, n_runs)
diff --git a/tests/quick_benchmark b/tests/quick_benchmark
new file mode 100755
index 00000000..a41c8991
--- /dev/null
+++ b/tests/quick_benchmark
@@ -0,0 +1,8 @@
+#!/usr/bin/env python3
+"""
+Quickly check the performance benchmark
+"""
+
+import test_baselines
+
+test_baselines.test_benchmark(do_save=False, repeats=5)
diff --git a/tests/test_baselines.py b/tests/test_baselines.py
index 45f8714c..5b41d991 100644
--- a/tests/test_baselines.py
+++ b/tests/test_baselines.py
@@ -8,29 +8,30 @@
 import starsim as ss
 
 baseline_filename  = sc.thisdir(__file__, 'baseline.json')
+multisim_filename  = sc.thisdir(__file__, 'baseline_multisim.json')
 benchmark_filename = sc.thisdir(__file__, 'benchmark.json')
 parameters_filename = sc.thisdir(ss.__file__, 'regression', f'pars_v{ss.__version__}.json')
 sc.options(interactive=False) # Assume not running interactively
 
 # Define the parameters
 pars = sc.objdict(
-    n_agents   = 10e3, # Number of agents
-    start      = 2000, # Starting year
-    n_years    = 20,   # Number of years to simulate
-    dt         = 0.2,  # Timestep
-    verbose    = 0,    # Don't print details of the run
-    rand_seed  = 2,    # Set a non-default seed
+    n_agents  = 10e3, # Number of agents
+    start     = 2000, # Starting year
+    dur       = 20,   # Number of years to simulate
+    dt        = 0.2,  # Timestep
+    verbose   = 0,    # Don't print details of the run
+    rand_seed = 2,    # Set a non-default seed
 )
 
 
-def make_sim(run=False):
+def make_sim(run=False, **kwargs):
     """
     Define a default simulation for testing the baseline. If run directly (not 
     via pytest), also plot the sim by default.
     """
     diseases = ['sir', 'sis']
     networks = ['random', 'mf', 'maternal']
-    sim = ss.Sim(pars=pars, networks=networks, diseases=diseases, demographics=True)
+    sim = ss.Sim(pars=pars | kwargs, networks=networks, diseases=diseases, demographics=True)
     
     # Optionally run and plot
     if run:
@@ -61,6 +62,31 @@ def save_baseline():
     return
 
 
+def multisim_baseline(save=False, n_runs=10, **kwargs):
+    """
+    Check or update the multisim baseline results; not part of the integration tests
+    """
+    word = 'Saving' if save else 'Checking'
+    sc.heading(f'{word} MultiSim baseline...')
+
+    # Make and run sims
+    kwargs.setdefault('verbose', ss.options.verbose/n_runs*2)
+    sim = make_sim(**kwargs)
+    msim = ss.MultiSim(base_sim=sim)
+    msim.run(n_runs)
+    summary = msim.summarize()
+
+    # Export results
+    if save:
+        sc.savejson(multisim_filename, summary)
+        print('Done.')
+    else:
+        baseline = sc.loadjson(multisim_filename)
+        ss.diff_sims(baseline, summary, die=False)
+
+    return
+
+
 def test_baseline():
     """ Compare the current default sim against the saved baseline """
     
@@ -121,7 +147,7 @@ def normalize_performance():
         # Time initialization
         t0 = sc.tic()
         sim = make_sim()
-        sim.initialize()
+        sim.init()
         t_init = sc.toc(t0, output=True)
 
         # Time running
@@ -147,7 +173,7 @@ def normalize_performance():
                 },
             'parameters': {
                 'n_agents': sim.pars.n_agents,
-                'n_years':  sim.pars.n_years,
+                'dur':      sim.pars.dur,
                 'dt':       sim.pars.dt,
                 },
             'cpu_performance': ratio,
diff --git a/tests/test_calibration.py b/tests/test_calibration.py
index 31999001..0d127d70 100644
--- a/tests/test_calibration.py
+++ b/tests/test_calibration.py
@@ -28,7 +28,7 @@ def make_sim():
         n_agents = n_agents,
         total_pop = 9980999,
         start = 1990,
-        n_years = 40,
+        dur = 40,
         diseases = [hiv],
         networks = [random, maternal],
         demographics = [pregnancy, death],
diff --git a/tests/test_demographics.py b/tests/test_demographics.py
index f297cbfe..0b9e4db4 100644
--- a/tests/test_demographics.py
+++ b/tests/test_demographics.py
@@ -15,7 +15,7 @@
 datadir = ss.root / 'tests/test_data'
 
 
-def test_nigeria(which='births', dt=1, start=1995, n_years=15, do_plot=False):
+def test_nigeria(which='births', dt=1, start=1995, dur=15, do_plot=False):
     """
     Make a Nigeria sim with demographic modules
     Switch between which='births' or 'pregnancy' to determine which demographic module to use
@@ -49,13 +49,13 @@ def test_nigeria(which='births', dt=1, start=1995, n_years=15, do_plot=False):
         dt=dt,
         total_pop=nga_pop_1995,
         start=start,
-        n_years=n_years,
+        dur=dur,
         people=ppl,
         demographics=demographics,
     )
 
     if do_plot:
-        sim.initialize()
+        sim.init()
         # Plot histograms of the age distributions - simulated vs data
         bins = np.arange(0, 101, 1)
         init_scale = nga_pop_1995 / n_agents
@@ -66,15 +66,15 @@ def test_nigeria(which='births', dt=1, start=1995, n_years=15, do_plot=False):
 
     sim.run()
 
-    end = start + n_years
+    stop = start + dur
     nigeria_popsize = pd.read_csv(datadir/'nigeria_popsize.csv')
-    data = nigeria_popsize[(nigeria_popsize.year >= start) & (nigeria_popsize.year <= end)]
+    data = nigeria_popsize[(nigeria_popsize.year >= start) & (nigeria_popsize.year <= stop)]
 
     nigeria_cbr = pd.read_csv(datadir/'nigeria_births.csv')
-    cbr_data = nigeria_cbr[(nigeria_cbr.Year >= start) & (nigeria_cbr.Year <= end)]
+    cbr_data = nigeria_cbr[(nigeria_cbr.Year >= start) & (nigeria_cbr.Year <= stop)]
 
     nigeria_cmr = pd.read_csv(datadir/'nigeria_cmr.csv')
-    cmr_data = nigeria_cmr[(nigeria_cmr.Year >= start) & (nigeria_cmr.Year <= end)]
+    cmr_data = nigeria_cmr[(nigeria_cmr.Year >= start) & (nigeria_cmr.Year <= stop)]
 
     # Tests
     if which == 'pregnancy':
@@ -95,28 +95,29 @@ def test_nigeria(which='births', dt=1, start=1995, n_years=15, do_plot=False):
 
     # Plots
     if do_plot:
+        tvec = sim.timevec
         fig, ax = plt.subplots(2, 2)
         ax = ax.ravel()
         ax[0].scatter(data.year, data.n_alive, alpha=0.5)
-        ax[0].plot(sim.yearvec, sim.results.n_alive, color='k')
+        ax[0].plot(tvec, sim.results.n_alive, color='k')
         ax[0].set_title('Population')
 
-        ax[1].plot(sim.yearvec, 1000 * sim.results.deaths.cmr / dt, label='Simulated CMR')
+        ax[1].plot(tvec, 1000 * sim.results.deaths.cmr, label='Simulated CMR')
         ax[1].scatter(cmr_data.Year, cmr_data.CMR, label='Data CMR')
         ax[1].set_title('CMR')
         ax[1].legend()
 
         if which == 'births':
-            ax[2].plot(sim.yearvec, sim.results.births.cbr / dt, label='Simulated CBR')
+            ax[2].plot(tvec, sim.results.births.cbr, label='Simulated CBR')
         elif which == 'pregnancy':
-            ax[2].plot(sim.yearvec, sim.results.pregnancy.cbr / dt, label='Simulated CBR')
+            ax[2].plot(tvec, sim.results.pregnancy.cbr, label='Simulated CBR')
         ax[2].scatter(cbr_data.Year, cbr_data.CBR, label='Data CBR')
         ax[2].set_title('CBR')
         ax[2].legend()
 
         if which == 'pregnancy':
-            ax[3].plot(sim.yearvec, sim.results.pregnancy.pregnancies / dt, label='Pregnancies')
-            ax[3].plot(sim.yearvec, sim.results.pregnancy.births / dt, label='Births')
+            ax[3].plot(tvec, sim.results.pregnancy.pregnancies, label='Pregnancies')
+            ax[3].plot(tvec, sim.results.pregnancy.births, label='Births')
             ax[3].set_title('Pregnancies and births')
             ax[3].legend()
 
@@ -128,7 +129,7 @@ def test_nigeria(which='births', dt=1, start=1995, n_years=15, do_plot=False):
 def test_constant_pop(do_plot=False):
     """ Test pars for constant pop size """
     sc.heading('Testing constant population size')
-    sim = ss.Sim(n_agents=10e3, birth_rate=10, death_rate=10/1010*1000, n_years=200, rand_seed=1).run()
+    sim = ss.Sim(n_agents=10e3, birth_rate=10, death_rate=10/1010*1000, dur=200, rand_seed=1).run()
     print("Check final pop size within 5% of starting pop")
     assert np.isclose(sim.results.n_alive[0], sim.results.n_alive[-1], rtol=0.05)
     print(f'✓ (final pop / starting pop={sim.results.n_alive[-1] / sim.results.n_alive[0]:.2f})')
@@ -156,7 +157,7 @@ def test_aging():
     n_agents = int(1e3)
     
     # With no demograhpics, people shouldn't age
-    s1 = ss.Sim(n_agents=n_agents).initialize()
+    s1 = ss.Sim(n_agents=n_agents).init()
     orig_ages = s1.people.age.raw.copy()
     orig_age = orig_ages.mean()
     s1.run()
diff --git a/tests/test_diseases.py b/tests/test_diseases.py
index 37542a2b..124f635d 100644
--- a/tests/test_diseases.py
+++ b/tests/test_diseases.py
@@ -26,7 +26,7 @@ def test_sir():
     sir = ss.SIR(sir_pars)
 
     # Change pars after creating the SIR instance
-    sir.pars.beta = {'random': 0.1}
+    sir.pars.beta = {'random': ss.rate(0.1)}
 
     # You can also change the parameters of the default lognormal distribution directly
     sir.pars.dur_inf.set(loc=5)
@@ -40,7 +40,7 @@ def test_sir():
     plt.figure()
     res = sim.results
     plt.stackplot(
-        sim.yearvec,
+        sim.timevec,
         res.sir.n_susceptible,
         res.sir.n_infected,
         res.sir.n_recovered,
@@ -128,7 +128,7 @@ def test_ncd():
 
     plt.figure()
     plt.stackplot(
-        sim.yearvec,
+        sim.timevec,
         ncd.results.n_not_at_risk,
         ncd.results.n_at_risk - ncd.results.n_affected,
         ncd.results.n_affected,
diff --git a/tests/test_dist.py b/tests/test_dist.py
index eaedd018..97c19bed 100644
--- a/tests/test_dist.py
+++ b/tests/test_dist.py
@@ -29,7 +29,7 @@ def plot_rvs(rvs, times=None, nrows=None):
 
 def make_sim():
     """ Make a tiny sim for initializing the distributions """
-    sim = ss.Sim(n_agents=100).initialize() # Need an empty sim to initialize properly
+    sim = ss.Sim(n_agents=100).init() # Need an empty sim to initialize properly
     return sim
 
 
@@ -38,7 +38,7 @@ def test_dist(m=m):
     """ Test the Dist class """
     sc.heading('Testing the basic Dist call')
     dist = ss.Dist(distname='random', name='test', strict=False)
-    dist.initialize()
+    dist.init()
     rvs = dist(m)
     print(rvs)
     assert 0 < rvs.min() < 1, 'Values should be between 0 and 1'
@@ -56,7 +56,7 @@ def test_custom_dists(n=n, do_plot=False):
     for name in ss.dist_list:
         func = getattr(ss, name)
         dist = func(name='test', strict=False)
-        dist.initialize()
+        dist.init()
         dists[name] = dist
         sc.tic()
         rvs[name] = dist.rvs(n)
@@ -85,7 +85,7 @@ def test_dists(n=n, do_plot=False):
         obj.a = sc.objdict()
         obj.a.mylist = [ss.random(), ss.Dist(distname='uniform', low=2, high=3)]
         obj.b = dict(d3=ss.weibull(c=2), d4=ss.constant(v=0.3))
-        dists = ss.Dists(obj).initialize(sim=make_sim())
+        dists = ss.Dists(obj).init(sim=make_sim())
         
         # Call each distribution twice
         for j in range(2):
@@ -114,8 +114,8 @@ def test_scipy(m=m):
     sc.heading('Testing SciPy distributions')
     
     # Make SciPy distributions in two different ways
-    dist1 = ss.Dist(dist=sps.expon, name='scipy', scale=2, strict=False).initialize() # Version 1: callable
-    dist2 = ss.Dist(dist=sps.expon(scale=2), name='scipy', strict=False).initialize() # Version 2: frozen
+    dist1 = ss.Dist(dist=sps.expon, name='scipy', scale=2, strict=False).init() # Version 1: callable
+    dist2 = ss.Dist(dist=sps.expon(scale=2), name='scipy', strict=False).init() # Version 2: frozen
     rvs1 = dist1(m)
     rvs2 = dist2(m)
     
@@ -136,14 +136,14 @@ def test_exceptions(m=m):
         dist(m) # Check that we can't call an uninitialized
     
     # Initialize and check we can't call repeatedly
-    dist.initialize(trace='test', sim=make_sim())
+    dist.init(trace='test', sim=make_sim())
     rvs = dist(m)
     with pytest.raises(ss.distributions.DistNotReadyError):
         dist(m) # Check that we can't call an already-used distribution
     
     # Check that we can with a non-strict Dist
     dist2 = ss.random(strict=False)
-    dist2.initialize(trace='test')
+    dist2.init(trace='test')
     rvs2 = sc.autolist()
     for i in range(2):
         rvs2 += dist2(m) # We should be able to call multiple times with no problem
@@ -161,7 +161,7 @@ def test_reset(m=m):
     sc.heading('Testing reset')
     
     # Create and draw two sets of random numbers
-    dist = ss.random(seed=533, strict=False).initialize()
+    dist = ss.random(seed=533, strict=False).init()
     r1 = dist.rvs(m)
     r2 = dist.rvs(m)
     assert all(r1 != r2)
@@ -200,8 +200,8 @@ def custom_loc(module, sim, uids):
         return out
     
     scale = 1
-    d1 = ss.normal(name='callable', loc=custom_loc, scale=scale).initialize(sim=sim)
-    d2 = ss.lognorm_ex(name='callable', mean=custom_loc, stdev=scale).initialize(sim=sim)
+    d1 = ss.normal(name='callable', loc=custom_loc, scale=scale).init(sim=sim)
+    d2 = ss.lognorm_ex(name='callable', mean=custom_loc, std=scale).init(sim=sim)
 
     uids = np.array([1, 3, 7, 9])
     draws1 = d1.rvs(uids)
@@ -223,7 +223,7 @@ def test_array(n=n):
     low  = np.array([1, 100]) # Low
     high = np.array([3, 125]) # High
 
-    d = ss.uniform(low=low, high=high, strict=False).initialize(slots=np.arange(uids.max()+1))
+    d = ss.uniform(low=low, high=high, strict=False).init(slots=np.arange(uids.max()+1))
     draws = d.rvs(uids)
     print(f'Uniform sample for uids {uids} returned {draws}')
 
@@ -245,7 +245,7 @@ def test_repeat_slot():
     high = low + 1
 
     # Draw values
-    d = ss.uniform(low=low, high=high, strict=False).initialize(slots=slots)
+    d = ss.uniform(low=low, high=high, strict=False).init(slots=slots)
     draws = d.rvs(uids)
     
     # Print and test
diff --git a/tests/test_interventions.py b/tests/test_interventions.py
index 1f7f05e7..f2f7d95e 100644
--- a/tests/test_interventions.py
+++ b/tests/test_interventions.py
@@ -30,10 +30,10 @@ def run_sir_vaccine(efficacy, leaky=True):
                 beta      = 6,       
           )
         ),
-        n_years = 10,
-        dt      = 0.01
+        dur = 10,
+        dt  = 0.01
     )
-    sim.initialize(verbose=False)
+    sim.init(verbose=False)
     
     # work out who to vaccinate
     in_trial = sim.people.sir.susceptible.uids
@@ -62,7 +62,7 @@ def run_sir_vaccine(efficacy, leaky=True):
     # run the simulation until sufficient cases
     old_cases = []
     for idx in range(1000):
-        sim.step()
+        sim.run_one_step()
         susc = sim.people.sir.susceptible.uids
         cases = np.setdiff1d(in_trial, susc)
         if len(cases) > total_cases:
@@ -97,7 +97,7 @@ def test_sir_vaccine_all_or_nothing():
 if __name__ == '__main__':
     T = sc.timer()
 
-    sir_vaccine_leaky   = test_sir_vaccine_leaky(leaky=True)
-    sir_vaccine_a_or_n  = test_sir_vaccine_all_or_nothing(leaky=False)
+    leaky  = test_sir_vaccine_leaky()
+    a_or_n = test_sir_vaccine_all_or_nothing()
 
     T.toc()
diff --git a/tests/test_networks.py b/tests/test_networks.py
index 1ceb0eae..186006df 100644
--- a/tests/test_networks.py
+++ b/tests/test_networks.py
@@ -28,7 +28,7 @@ def test_manual():
 
     # Create a maternal network
     sim = ss.Sim(n_agents=n_agents)
-    sim.initialize()
+    sim.init()
     nw2 = ss.MaternalNet()
     nw2.init_pre(sim)
     nw2.add_pairs(mother_inds=[1, 2, 3], unborn_inds=[100, 101, 102], dur=[1, 1, 1])
@@ -43,15 +43,15 @@ def test_random():
 
     # Manual creation
     nw1 = ss.RandomNet()
-    ss.Sim(n_agents=small, networks=nw1, copy_inputs=False).initialize() # This initializes the network
+    ss.Sim(n_agents=small, networks=nw1, copy_inputs=False).init() # This initializes the network
     
     # Automatic creation as part of sim
-    s2 = ss.Sim(n_agents=small, networks='random').initialize()
+    s2 = ss.Sim(n_agents=small, networks='random').init()
     nw2 = s2.networks[0]
     
     # Increase the number of contacts
     nwdict = dict(type='random', n_contacts=20)
-    s3 = ss.Sim(n_agents=small, networks=nwdict).initialize()
+    s3 = ss.Sim(n_agents=small, networks=nwdict).init()
     nw3 = s3.networks[0]
     
     # Checks
@@ -99,17 +99,17 @@ def test_ER(n, p, nw, alpha=0.01):
     # Manual creation
     p = 0.1
     nw1 = ss.ErdosRenyiNet(p=p)
-    ss.Sim(n_agents=small, networks=nw1, copy_inputs=False).initialize() # This initializes the network
+    ss.Sim(n_agents=small, networks=nw1, copy_inputs=False).init() # This initializes the network
     test_ER(small, p, nw1)
 
     # Automatic creation as part of sim
-    s2 = ss.Sim(n_agents=small, networks='erdosrenyi').initialize()
+    s2 = ss.Sim(n_agents=small, networks='erdosrenyi').init()
     nw2 = s2.networks[0]
     
     # Larger example with higher p
     p=0.6
     nwdict = dict(type='erdosrenyi', p=p)
-    s3 = ss.Sim(n_agents=medium, networks=nwdict).initialize()
+    s3 = ss.Sim(n_agents=medium, networks=nwdict).init()
     nw3 = s3.networks[0]
     test_ER(medium, p, nw3)
     
@@ -126,7 +126,7 @@ def test_disk():
 
     # Visualize the path of agents
     nw1 = ss.DiskNet()
-    s1 = ss.Sim(n_agents=5, n_years=50, networks=nw1, copy_inputs=False).initialize() # This initializes the network
+    s1 = ss.Sim(n_agents=5, dur=50, networks=nw1, copy_inputs=False).init() # This initializes the network
 
     if sc.options.interactive:
         # Visualize motion:
@@ -139,7 +139,7 @@ def test_disk():
         cmap = mpl.colormaps['plasma']
         colors = cmap(np.linspace(0, 1, s1.pars.n_agents))
         ax.scatter(nw1.x, nw1.y, s=50, c=colors)
-        for i in range(s1.pars.n_years):
+        for i in range(s1.pars.dur):
             ax.plot([0,1,1,0,0], [0,0,1,1,0], 'k-', lw=1)
             ax.quiver(nw1.x, nw1.y, vdt * np.cos(nw1.theta), vdt * np.sin(nw1.theta), color=colors)
             ax.set_aspect('equal', adjustable='box') #ax.set_xlim([0,1]); ax.set_ylim([0,1])
@@ -147,7 +147,7 @@ def test_disk():
 
     # Simulate SIR on a DiskNet
     nw2 = ss.DiskNet(r=0.15, v=0.05)
-    s2 = ss.Sim(n_agents=small, networks=nw2, diseases='sir').initialize() # This initializes the network
+    s2 = ss.Sim(n_agents=small, networks=nw2, diseases='sir').init() # This initializes the network
     s2.run()
 
     if sc.options.interactive:
@@ -165,11 +165,11 @@ def test_static():
     n = 100
     nc = p*n
     nd1 = dict(type='static', p=p)
-    nw1 = ss.Sim(n_agents=n, networks=nd1).initialize().networks[0]
+    nw1 = ss.Sim(n_agents=n, networks=nd1).init().networks[0]
     
     # Create with n_contacts
     nd2 = dict(type='static', n_contacts=nc)
-    nw2 = ss.Sim(n_agents=n, networks=nd2).initialize().networks[0]
+    nw2 = ss.Sim(n_agents=n, networks=nd2).init().networks[0]
     
     # Check
     assert len(nw1) == len(nw2), 'Networks should be the same length'
diff --git a/tests/test_other.py b/tests/test_other.py
index 4ab4309d..78590e9f 100644
--- a/tests/test_other.py
+++ b/tests/test_other.py
@@ -51,14 +51,15 @@ def test_microsim(do_plot=False):
         people=ss.People(small),
         networks=[ss.MFNet(), ss.MaternalNet()],
         demographics=ss.Pregnancy(),
-        diseases=hiv
+        diseases=hiv,
+        copy_inputs = False, # So we can reuse hiv
     )
-    sim.initialize()
+    sim.init()
     sim.run()
 
     if do_plot:
         pl.figure()
-        pl.plot(sim.tivec, sim.results.hiv.n_infected)
+        pl.plot(hiv.timevec, hiv.results.n_infected)
         pl.title('HIV number of infections')
 
     return sim
@@ -81,10 +82,10 @@ def init_debut(module, sim, uids):
     gon = ss.Gonorrhea(pars=gon_pars)
 
     sim = ss.Sim(pars=sim_pars, diseases=[gon])
-    sim.initialize()
+    sim.init()
     sim.run()
     pl.figure()
-    pl.plot(sim.tivec, sim.results.gonorrhea.n_infected)
+    pl.plot(sim.timevec, sim.results.gonorrhea.n_infected)
     pl.title('Number of gonorrhea infections')
 
     return sim
@@ -118,7 +119,7 @@ def test_arrs():
 
 def test_deepcopy():
     s1 = ss.Sim(pars=dict(diseases='sir', networks='embedding'), n_agents=small)
-    s1.initialize()
+    s1.init()
     
     s2 = sc.dcp(s1)
 
@@ -135,7 +136,7 @@ def test_deepcopy():
 
 def test_deepcopy_until():
     s1 = ss.Sim(pars=dict(diseases='sir', networks='embedding'), n_agents=small)
-    s1.initialize()
+    s1.init()
 
     s1.run(until=5)
 
diff --git a/tests/test_plugins.py b/tests/test_plugins.py
index 2729e50b..93d34e27 100644
--- a/tests/test_plugins.py
+++ b/tests/test_plugins.py
@@ -28,6 +28,7 @@ def __init__(self, pars=None, **kwargs):
 
     def step(self):
         """ Specify HIV-syphilis interactions """
+        
         diseases = self.sim.diseases
         syph = diseases.syphilis
         hiv = diseases.hiv
@@ -59,7 +60,7 @@ def __init__(self, year=2020, prob=0.8):
 
     def step(self):
         sim = self.sim
-        if sim.year > self.year:
+        if sim.now > self.year:
             syphilis = sim.diseases.syphilis
 
             # Define who is eligible for treatment
@@ -79,8 +80,8 @@ def step(self):
 def make_args():
     """ Make people, HIV, syphilis, and network """
     pars = dict(n_agents=2000, verbose=0)
-    mf = ss.MFNet(duration=ss.lognorm_ex(mean=5, stdev=0.5))
-    hiv = ss.HIV(beta={'mf': [0.0008, 0.0004]}, init_prev=0.2)
+    mf = ss.MFNet(duration=ss.lognorm_ex(mean=5, std=0.5)) # TODO: think about whether these should be ss.dur(); currently they are not since stored in natural units with -self.dt
+    hiv = ss.HIV(beta={'mf': [0.0008, 0.0004]}, init_prev=0.2) # TODO: beta should wrap the other way
     syph = ss.Syphilis(beta={'mf': [0.1, 0.05]}, init_prev=0.05)
     args = dict(pars=pars, networks=mf, diseases=[hiv, syph])
     return args
@@ -119,7 +120,7 @@ def test_connectors(do_plot=False):
         pl.figure()
         
         pl.subplot(2,1,1)
-        x = sims.con.yearvec
+        x = sims.con.timevec
         for label,res in results.items():
             pl.plot(x, res.syphilis, label=label)
         pl.title('Syphilis infections')
diff --git a/tests/test_randomness.py b/tests/test_randomness.py
index 1a9452cb..82169521 100644
--- a/tests/test_randomness.py
+++ b/tests/test_randomness.py
@@ -18,10 +18,10 @@ def make_dist(name='test', **kwargs):
 
 def make_dists(**kwargs):
     """ Make a Dists object with two distributions in it """
-    sim = ss.Sim(n_agents=100).initialize() # Need an empty sim to initialize properly
+    sim = ss.Sim(n_agents=100).init() # Need an empty sim to initialize properly
     distlist = [make_dist(), make_dist()]
     dists = ss.Dists(distlist)
-    dists.initialize(sim=sim)
+    dists.init(sim=sim)
     return dists
 
 
@@ -152,7 +152,7 @@ def step(self):
             # Reset the random states
             p = sir.pars
             for dist in [p.dur_inf, p.p_death]:
-                dist.jump(sim.ti+1)
+                dist.jump_dt(force=True) # Already has correct dt value, but we need to force going back in time
 
         return
 
@@ -189,7 +189,7 @@ def test_worlds(do_plot=False):
     
     pars = dict(
         start = 2000,
-        end = 2100,
+        stop = 2100,
         n_agents = 200,
         verbose = 0.05,
         diseases = dict(
@@ -250,7 +250,7 @@ def test_independence(do_plot=False, thresh=0.1):
             dict(type='mf', debut=ss.constant(0), participation=0.5), # To avoid age correlations
         ]
     )
-    sim.initialize()
+    sim.init()
     
     # Assemble measures
     st = sim.people.states
diff --git a/tests/test_samples.py b/tests/test_samples.py
index 105662a4..ab502043 100644
--- a/tests/test_samples.py
+++ b/tests/test_samples.py
@@ -14,7 +14,7 @@ def get_outputs(p_death):
         ppl = ss.People(1000)
         network = ss.RandomNet(n_contacts=ss.poisson(5))
         sir = ss.SIR(pars={'p_death':p_death})
-        sim = ss.Sim(people=ppl, networks=network, diseases=sir, rand_seed=0, n_years=5)
+        sim = ss.Sim(people=ppl, networks=network, diseases=sir, rand_seed=0, dur=5)
         sim.run(verbose=0)
         df = sim.export_df()
         summary = {}
diff --git a/tests/test_sim.py b/tests/test_sim.py
index b45e127c..e510ca64 100644
--- a/tests/test_sim.py
+++ b/tests/test_sim.py
@@ -80,7 +80,7 @@ def test_api():
     # Check different ways of setting a distribution
     kw = dict(n_agents=n_agents, networks='random')
     d1 = ss.lognorm_ex(10) # Create a distribution with an argument
-    d2 = ss.lognorm_ex(mean=10, stdev=2) # Create a distribution with kwargs
+    d2 = ss.lognorm_ex(mean=10, std=2) # Create a distribution with kwargs
     d3 = ss.normal(loc=10) # Create a different type of distribution
     
     # Check specifying dist with a scalar
@@ -90,7 +90,7 @@ def test_api():
     
     # Check specifying dist with a list and dict
     s6 = ss.Sim(diseases=dict(type='sir', dur_inf=[10,2]), **kw).run() # Supply values as a list
-    s7 = ss.Sim(diseases=dict(type='sir', dur_inf=dict(mean=10, stdev=2)), **kw).run() # Supply values as a dict
+    s7 = ss.Sim(diseases=dict(type='sir', dur_inf=dict(mean=10, std=2)), **kw).run() # Supply values as a dict
     s8 = ss.Sim(diseases=dict(type='sir', dur_inf=d2), **kw).run() # Supply as a distribution
     ss.check_sims_match(s6, s7, s8), 'Sims should match'
     
@@ -101,7 +101,7 @@ def test_api():
     
     # Check that Bernoulli distributions can't be changed
     with pytest.raises(TypeError):
-        ss.Sim(diseases=dict(type='sir', init_prev=dict(type='normal', loc=10)), **kw).initialize()
+        ss.Sim(diseases=dict(type='sir', init_prev=dict(type='normal', loc=10)), **kw).init()
     
     return s1
 
@@ -116,11 +116,11 @@ def jump_age(sim):
             sim.people.age[:] = sim.people.age[:] + 1000
     
     # Specify parameters as a dictionary
-    p = dict(
+    pars1 = dict(
         n_agents = 1000,
         label = 'v1',
         verbose = 'brief',
-        end = 2020,
+        stop = 2020,
         networks = [
             ss.RandomNet(name='random1', n_contacts=6),
             dict(type='random', name='random2', n_contacts=4)
@@ -137,7 +137,7 @@ def jump_age(sim):
     )
     
     # Test with explicit initialization
-    pars = ss.SimPars(n_agents=1000, label='v1', verbose='brief', end=2020)
+    pars2 = ss.SimPars(n_agents=1000, label='v1', verbose='brief', stop=2020)
     
     net1 = ss.RandomNet(name='random1', n_contacts=6)
     net2 = ss.RandomNet(name='random2', n_contacts=4)
@@ -155,8 +155,8 @@ def jump_age(sim):
     interventions = ss.ndict(int1)
 
     # Assemble
-    s1 = ss.Sim(p)
-    s2 = ss.Sim(pars=pars, networks=networks, diseases=diseases, demographics=demographics, interventions=interventions)
+    s1 = ss.Sim(pars1)
+    s2 = ss.Sim(pars=pars2, networks=networks, diseases=diseases, demographics=demographics, interventions=interventions)
     
     # Run
     s1.run()
diff --git a/tests/test_syphilis.py b/tests/test_syphilis.py
index e36b5ed1..a566e239 100644
--- a/tests/test_syphilis.py
+++ b/tests/test_syphilis.py
@@ -37,7 +37,7 @@ def make_syph_sim(dt=1, n_agents=500):
         dt=dt,
         total_pop=93963392,
         start=1990,
-        n_years=40,
+        dur=40,
         people=ppl,
         diseases=syph,
         networks=ss.ndict(mf, maternal),
@@ -91,13 +91,14 @@ def test_syph(dt=1, n_agents=500, do_plot=False):
 
     # Check plots
     burnin = 0
-    pi = int(burnin/sim.dt)
+    pi = int(burnin/dt)
     
     if do_plot:
+        tvec = sim.timevec[pi:]
         fig, ax = plt.subplots(2, 2)
         ax = ax.ravel()
         ax[0].stackplot(
-            sim.yearvec[pi:],
+            tvec,
             # sim.results.syphilis.n_susceptible[pi:],
             sim.results.syphilis.n_congenital[pi:],
             sim.results.syphilis.n_exposed[pi:],
@@ -108,13 +109,13 @@ def test_syph(dt=1, n_agents=500, do_plot=False):
         )
         ax[0].legend(['Congenital', 'Exposed', 'Primary', 'Secondary', 'Latent', 'Tertiary'], loc='lower right')
     
-        ax[1].plot(sim.yearvec[pi:], sim.results.syphilis.prevalence[pi:])
+        ax[1].plot(tvec, sim.results.syphilis.prevalence[pi:])
         ax[1].set_title('Syphilis prevalence')
     
-        ax[2].plot(sim.yearvec[pi:], sim.results.n_alive[pi:])
+        ax[2].plot(tvec, sim.results.n_alive[pi:])
         ax[2].set_title('Population')
     
-        ax[3].plot(sim.yearvec[pi:], sim.results.syphilis.new_infections[pi:])
+        ax[3].plot(tvec, sim.results.syphilis.new_infections[pi:])
         ax[3].set_title('New infections')
     
         fig.tight_layout()
@@ -156,10 +157,12 @@ def test_syph_intvs(dt=1, n_agents=500, do_plot=False):
         sim_base.run()
 
         burnin = 10
-        pi = int(burnin/sim_base.dt)
+        syph_b = sim_base.diseases.syphilis
+        syph_i = sim_intv.diseases.syphilis
+        pi = int(burnin/syph_b.dt)
         plt.figure()
-        plt.plot(sim_base.yearvec[pi:], sim_base.results.syphilis.prevalence[pi:], label='Baseline')
-        plt.plot(sim_intv.yearvec[pi:], sim_intv.results.syphilis.prevalence[pi:], label='S&T')
+        plt.plot(syph_b.timevec[pi:], syph_b.results.prevalence[pi:], label='Baseline')
+        plt.plot(syph_i.timevec[pi:], syph_i.results.prevalence[pi:], label='S&T')
         plt.axvline(x=2020, color='k', ls='--')
         plt.title('Syphilis prevalence')
         plt.legend()
diff --git a/tests/test_time.py b/tests/test_time.py
new file mode 100644
index 00000000..5c945db6
--- /dev/null
+++ b/tests/test_time.py
@@ -0,0 +1,139 @@
+"""
+Test different time units and timesteps
+"""
+
+import numpy as np
+import sciris as sc
+import starsim as ss
+
+small = 100
+medium = 1000
+sc.options(interactive=False)
+
+
+# %% Define the tests
+def test_ratio():
+    sc.heading('Test behavior of time_ratio()')
+    
+    assert ss.time_ratio() == 1.0
+    assert ss.time_ratio(dt1=1.0, dt2=0.1) == 10.0
+    assert ss.time_ratio(dt1=0.5, dt2=5) == 0.1
+    assert ss.time_ratio(dt1=0.5, dt2=5) == 0.1
+    
+    assert ss.time_ratio(unit1='year', unit2='day') == 365
+    assert ss.time_ratio(unit1='day', unit2='year') == 1/365
+    
+    return
+
+
+def test_classes():
+    sc.heading('Test behavior of dur() and rate()')
+    
+    # Test duration dt
+    d1 = ss.dur(2)
+    d2 = ss.dur(3)
+    d3 = ss.dur(2, parent_dt=0.1)
+    d4 = ss.dur(3, parent_dt=0.2)
+    for d in [d1,d2,d3,d4]: d.init()
+    
+    assert d1 + d2 == 2+3
+    assert d3 + d4 == 2/0.1 + 3/0.2
+    assert d3 * 2 == 2/0.1*2
+    assert d3 / 2 == 2/0.1/2
+    
+    # Test rate dt
+    r1 = ss.rate(2)
+    r2 = ss.rate(3)
+    r3 = ss.rate(2, parent_dt=0.1)
+    r4 = ss.rate(3, parent_dt=0.2)
+    for r in [r1,r2,r3,r4]: r.init()
+    
+    assert r1 + r2 == 2 + 3
+    assert r3 + r4 == 2*0.1 + 3*0.2
+    assert r3 * 2 == 2*0.1*2
+    assert r3 / 2 == 2*0.1/2
+    
+    # Test duration units
+    d5 = ss.dur(2, unit='year').init(parent_unit='day')
+    d6 = ss.dur(3, unit='day').init(parent_unit='day')
+    assert d5 + d6 == 2*365+3
+    
+    # Test rate units
+    rval = 0.7
+    r5 = ss.rate(rval, unit='week').init(parent_unit='day')
+    r6 = ss.rate(rval).init(parent_dt=0.1)
+    assert np.isclose(r5.x, rval/7) # A limitation of this approach, not exact!
+    assert np.isclose(r6.x, rval/10)
+    
+    # Test time_prob
+    tpval = 0.1
+    tp0 = ss.time_prob(tpval).init(parent_dt=1.0)
+    tp1 = ss.time_prob(tpval).init(parent_dt=0.5)
+    tp2 = ss.time_prob(tpval).init(parent_dt=2)
+    assert np.isclose(tp0.x, tpval)
+    assert np.isclose(tp1.x, tpval/2, rtol=0.1)
+    assert np.isclose(tp2.x, tpval*2, rtol=0.1)
+    assert tp1.x > tpval/2
+    assert tp2.x < tpval*2
+    
+    return d3, d4, r3, r4, tp1
+    
+
+def test_units(do_plot=False):
+    sc.heading('Test behavior of year vs day units')
+    
+    pars = dict(
+        diseases = dict(type='sis', init_prev=0.1),
+        networks = 'random',
+        n_agents = small,
+    )
+    
+    sims = sc.objdict()
+    sims.y = ss.Sim(pars, unit='year', label='Year', start=2000, stop=2002, dt=1/365)
+    sims.d = ss.Sim(pars, unit='day', label='Day', start='2000-01-01', stop='2002-01-01', dt=1)
+    
+    for sim in sims.values():
+        sim.run()
+        if do_plot:
+            sim.plot()
+    
+    # Uncomment this test once it might potentially pass lol
+    rtol = 0.01
+    vals = [sim.summary.sis_cum_infections for sim in [sims.y, sims.d]]
+    # assert np.isclose(*vals, rtol=rtol), f'Values for cum_infections do not match ({vals})'
+        
+    return sims
+
+
+def test_multi_timestep(do_plot=False):
+    sc.heading('Test behavior of different modules having different timesteps')
+    
+    pars = dict(
+        diseases = ss.SIS(unit='day', dt=1.0, init_prev=0.1, beta=ss.beta(0.01)),
+        # demographics = ss.Births(unit='year', dt=0.25),
+        networks = ss.RandomNet(unit='week'),
+        n_agents = small,
+    )
+    
+    sim = ss.Sim(pars, unit='day', dt=2, start='2000-01-01', stop='2002-01-01')
+    sim.run()
+    
+    if do_plot:
+        sim.plot()
+        
+    return sim
+
+
+# %% Run as a script
+if __name__ == '__main__':
+    do_plot = True
+    sc.options(interactive=do_plot)
+    
+    T = sc.timer()
+    
+    o1 = test_ratio()
+    o2 = test_classes()
+    o3 = test_units(do_plot)
+    o4 = test_multi_timestep(do_plot)
+    
+    T.toc()
diff --git a/tests/update_baseline.py b/tests/update_baseline.py
index 67af0b86..39b85147 100755
--- a/tests/update_baseline.py
+++ b/tests/update_baseline.py
@@ -1,10 +1,9 @@
 #!/usr/bin/env python3
-
-'''
+"""
 Run this script to regenerate the baseline.
-'''
+"""
 
-import test_baselines as tb
+import test_baselines
 
-tb.save_baseline()
-tb.test_benchmark(do_save=True, repeats=5)
+test_baselines.save_baseline()
+test_baselines.test_benchmark(do_save=True, repeats=5)