Repository: starsimhub/starsim (language: python)

Issue:
Starsim tries to attempt an invalid modulo operation
`  File "/Users/devengokhale/Documents/GitHub/starsim/starsim/sim.py", line 292, in start_step
    if not (self.ti % int(1.0 / self.verbose)):
            ~~~~~~~~^~~~~~~~~~~~~~~~~~~~~~~~~
ZeroDivisionError: integer modulo by zero`


Rc2.1.0 time
Adds an `ss.Time()` class that handles time units and conversions.


Discussion / hints:

Fix patch:
diff --git a/.pylintrc b/.pylintrc
index 2f032f13..aca48aca 100644
--- a/.pylintrc
+++ b/.pylintrc
@@ -86,6 +86,7 @@ disable=raw-checker-failed,
         useless-suppression,
         deprecated-pragma,
         use-symbolic-message-instead,
+        C0103, # Starsim: allow lowercase constants
         C0321, # Starsim: multiple statements on a line
         C0325, # Starsim: allow parentheses in if
         C0411, # Starsim: allow different import order
diff --git a/CHANGELOG.rst b/CHANGELOG.rst
index 3edf66af..e43777da 100644
--- a/CHANGELOG.rst
+++ b/CHANGELOG.rst
@@ -7,6 +7,38 @@ What's new
 All notable changes to the codebase are documented in this file. Changes that may result in differences in model output, or are required in order to run an old parameter set with the current version, are flagged with the term "Regression information".
 
 
+Version 2.1.0 (2024-11-07)
+---------------------------
+
+Summary
+~~~~~~~
+- Time in simulations is now handled by an ``ss.Time()`` class, which unifies how time is represented between the ``Sim`` and each module.
+- In addition to networks, there is now a new way of implementing disease transmission via mixing pools.
+
+Time
+~~~~
+- Time handling now performed by the ``ss.Time()`` class. This has inputs similar to before (``start``, ``stop``, ``unit``, ``dt``, with ``dur`` still available as a sim input). However, instead of the previous ``timevec`` and ``abs_tvec`` arrays, there are now multiple ways of representing time (including ``datevec`` and ``yearvec``), regardless of what the inputs were.
+- Dates are now represented in a native format, ``ss.date``, that is based on ``pd.Timestamp``.
+
+Mixing pools
+~~~~~~~~~~~~
+- Adds a new approach to disease transmission called mixing pools. A mixing pool is a "mean field" coupling wherein susceptible agents are exposed to the average infectious agent. The user can create a single mixing pool using the ``ss.MixingPool`` class, or create many pools using ``MixingPools``. Such mixing pools could be used to simulate contact matrices, for example as published by Prem et al.
+- There is a new ``ss.Route`` class, which is the base class for ``ss.Network`` and ``ss.MixingPool``.
+
+Other changes
+~~~~~~~~~~~~~
+- Demographic modules have been updated to fix various bugs around different time units.
+- The method for hashing distribution trace strings into seeds has changed, meaning that results will be stochastically different compared to Starsim v2.0.
+- Fixed a bug with how timepars were updated in parameters.
+- There is a new ``ss.Base`` class, which both ``ss.Sim`` and ``ss.Module`` inherit from.
+- Results now print as a single line rather than the full array. The latter is available as ``result.disp()``.
+- ``sim.to_df()`` now works even if different modules have different numbers of timepoints.
+- The ``timepars`` module has been renamed to ``time``.
+- In demographics modules, ``units`` has been renamed ``rate_units``.
+- There are two new options, ``ss.options.date_sep`` and  ``ss.options.license``. The former sets the date separator (default ``.``, e.g. ``2024.04.0.4``), and the latter sets if the license prints when Starsim is imported.
+- *GitHub info*: PR `724 <https://github.com/starsimhub/starsim/pull/724>`_
+
+
 Version 2.0.0 (2024-10-01)
 --------------------------
 
diff --git a/docs/api/index.rst b/docs/api/index.rst
index 2177c1dd..67630e47 100644
--- a/docs/api/index.rst
+++ b/docs/api/index.rst
@@ -13,6 +13,7 @@ This page lists Starsim's full API, including all functions and classes.
    :nosignatures:
 
    arrays
+   calibration
    demographics
    disease
    distributions
@@ -28,4 +29,5 @@ This page lists Starsim's full API, including all functions and classes.
    samples
    settings
    sim
+   time
    utils
diff --git a/docs/tutorials.rst b/docs/tutorials.rst
index a02506c2..b121e51d 100644
--- a/docs/tutorials.rst
+++ b/docs/tutorials.rst
@@ -9,8 +9,8 @@ User tutorials
     tutorials/tut_intro.ipynb
     tutorials/tut_buildsim.ipynb
     tutorials/tut_demographics.ipynb
-    tutorials/tut_networks.ipynb
     tutorials/tut_diseases.ipynb
+    tutorials/tut_transmission.ipynb
     tutorials/tut_interventions.ipynb
 
 Developer tutorials
diff --git a/docs/tutorials/dev_tut_analyzers.ipynb b/docs/tutorials/dev_tut_analyzers.ipynb
index a764ef39..d5e99753 100644
--- a/docs/tutorials/dev_tut_analyzers.ipynb
+++ b/docs/tutorials/dev_tut_analyzers.ipynb
@@ -88,11 +88,9 @@
     "import starsim as ss\n",
     "import pandas as pd\n",
     "\n",
-    "class HIV_in_pregnancy(ss.Analyzer):\n",
+    "class HIV_preg(ss.Analyzer):\n",
     "    def __init__(self, *args, **kwargs):\n",
     "        super().__init__(*args, **kwargs)\n",
-    "        self.requires = [ss.HIV, ss.Pregnancy]\n",
-    "        self.name = 'hiv_in_pregnancy'\n",
     "        return\n",
     "    \n",
     "    def init_results(self):\n",
@@ -102,7 +100,8 @@
     "        )\n",
     "        return\n",
     "\n",
-    "    def apply(self, sim):\n",
+    "    def step(self):\n",
+    "        sim = self.sim\n",
     "        ti = sim.ti\n",
     "        hiv = sim.diseases.hiv\n",
     "        pregnant = sim.demographics.pregnancy.pregnant\n",
@@ -112,16 +111,16 @@
     "\n",
     "pregnancy = ss.Pregnancy(pars=dict(fertility_rate=pd.read_csv('test_data/nigeria_asfr.csv')))\n",
     "hiv = ss.HIV(beta={'mfnet':[0.5,0.25]})\n",
-    "sim = ss.Sim(diseases=hiv, networks='mfnet', demographics=pregnancy, analyzers=HIV_in_pregnancy())\n",
+    "sim = ss.Sim(diseases=hiv, networks='mfnet', demographics=pregnancy, analyzers=HIV_preg())\n",
     "sim.run()\n",
-    "print(f'Total infections among pregnant women: {sim.results.hiv_in_pregnancy.new_infections_pregnancy.sum()}')\n"
+    "print(f'Total infections among pregnant women: {sim.results.hiv_preg.new_infections_pregnancy.sum()}')\n"
    ]
   },
   {
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "Analyzers are ideal for adding custom results, and because they get added to the sim in the same way as any other result, they also get automatically exported in the same format, e.g. using `sim.export_df()`.\n",
+    "Analyzers are ideal for adding custom results, and because they get added to the sim in the same way as any other result, they also get automatically exported in the same format, e.g. using `sim.to_df()`.\n",
     "\n",
     "Here's a plot of the results from our HIV in pregnancy analyzer:"
    ]
@@ -137,7 +136,7 @@
     "res = sim.results\n",
     "\n",
     "plt.figure()\n",
-    "plt.plot(res.timevec, res.hiv_in_pregnancy.new_infections_pregnancy)\n",
+    "plt.plot(res.timevec, res.hiv_preg.new_infections_pregnancy)\n",
     "plt.title('HIV infections acquired during pregnancy')\n",
     "plt.show();"
    ]
@@ -145,7 +144,7 @@
  ],
  "metadata": {
   "kernelspec": {
-   "display_name": "Python 3 (ipykernel)",
+   "display_name": "base",
    "language": "python",
    "name": "python3"
   },
@@ -159,7 +158,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.11.7"
+   "version": "3.12.2"
   },
   "pycharm": {
    "stem_cell": {
diff --git a/docs/tutorials/dev_tut_nonstandard.ipynb b/docs/tutorials/dev_tut_nonstandard.ipynb
index fef9c4f9..74cbff94 100644
--- a/docs/tutorials/dev_tut_nonstandard.ipynb
+++ b/docs/tutorials/dev_tut_nonstandard.ipynb
@@ -105,7 +105,7 @@
     "        res = self.results\n",
     "        kw = dict(lw=2, alpha=0.7)\n",
     "        for rkey in ['S', 'I']:\n",
-    "            plt.plot(self.timevec, res[rkey], label=res[rkey].label, **kw)\n",
+    "            plt.plot(res.timevec, res[rkey], label=res[rkey].label, **kw)\n",
     "        plt.legend(frameon=False)\n",
     "        plt.xlabel('Time')\n",
     "        plt.ylabel('Number of people')\n",
@@ -193,7 +193,7 @@
     "    def init_post(self, bins=10):\n",
     "        \"\"\" Define custom model attributes \"\"\"\n",
     "        super().init_post()\n",
-    "        self.npts = self.sim.npts # Number of timepoints\n",
+    "        self.npts = len(self.sim) # Number of timepoints\n",
     "        self.n_agents = len(sim.people) # Number of agents\n",
     "        self.wealth = np.ones(self.n_agents) # Initial wealth of each agent\n",
     "        self.bins = np.arange(bins+1) # Bins used for plotting\n",
@@ -275,7 +275,7 @@
  ],
  "metadata": {
   "kernelspec": {
-   "display_name": "Python 3 (ipykernel)",
+   "display_name": "base",
    "language": "python",
    "name": "python3"
   },
@@ -289,7 +289,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.11.7"
+   "version": "3.12.2"
   },
   "pycharm": {
    "stem_cell": {
diff --git a/docs/tutorials/tut_demographics.ipynb b/docs/tutorials/tut_demographics.ipynb
index 6f3fbe21..c74e808b 100644
--- a/docs/tutorials/tut_demographics.ipynb
+++ b/docs/tutorials/tut_demographics.ipynb
@@ -11,7 +11,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 2,
+   "execution_count": 1,
    "metadata": {
     "collapsed": false
    },
@@ -40,7 +40,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 3,
+   "execution_count": 2,
    "metadata": {
     "collapsed": false
    },
@@ -67,7 +67,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 4,
+   "execution_count": 3,
    "metadata": {
     "collapsed": false
    },
@@ -107,7 +107,7 @@
     "pregnancy = ss.Pregnancy(pars={'fertility_rate': fertility_rates})\n",
     "\n",
     "death_rates = pd.read_csv('test_data/nigeria_deaths.csv')\n",
-    "death = ss.Deaths(pars={'death_rate': death_rates, 'units': 1})\n",
+    "death = ss.Deaths(pars={'death_rate': death_rates, 'rate_units': 1})\n",
     "\n",
     "demographics = [pregnancy, death]\n",
     "\n",
@@ -127,8 +127,9 @@
     "\n",
     "# Plot the overall population size - simulated vs data\n",
     "fig, ax = plt.subplots(1, 1)\n",
+    "res = sim.results\n",
     "ax.scatter(data.year, data.n_alive, alpha=0.5, label='Data')\n",
-    "ax.plot(sim.timevec, sim.results.n_alive, color='k', label='Model')\n",
+    "ax.plot(res.timevec, res.n_alive, color='k', label='Model')\n",
     "ax.legend()\n",
     "ax.set_title('Population')\n",
     "plt.show();"
@@ -149,7 +150,7 @@
  ],
  "metadata": {
   "kernelspec": {
-   "display_name": "Python 3 (ipykernel)",
+   "display_name": "base",
    "language": "python",
    "name": "python3"
   },
@@ -163,7 +164,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.11.7"
+   "version": "3.12.2"
   },
   "pycharm": {
    "stem_cell": {
diff --git a/docs/tutorials/tut_diseases.ipynb b/docs/tutorials/tut_diseases.ipynb
index 6b9ec07b..d85cbbf2 100644
--- a/docs/tutorials/tut_diseases.ipynb
+++ b/docs/tutorials/tut_diseases.ipynb
@@ -4,7 +4,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "# T5 - Diseases\n",
+    "# T4 - Diseases\n",
     "\n",
     "Diseases are the cornerstone of almost any Starsim analysis. This is where you add details about what disease you are modeling, including transmissibility, natural history, and mortality outcomes.      \n",
     "\n",
diff --git a/docs/tutorials/tut_networks.ipynb b/docs/tutorials/tut_networks.ipynb
deleted file mode 100644
index 6a75c8eb..00000000
--- a/docs/tutorials/tut_networks.ipynb
+++ /dev/null
@@ -1,211 +0,0 @@
-{
- "cells": [
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "# T4 - Networks\n",
-    "\n",
-    "Many infectious diseases transmit when people come into contact with one another. In the basic SIR example that we covered in the first two tutorials, we used a 'random' contact network, which creates random pairings of people. Generally speaking, the structure of network that you use depends on what disease you're modeling. In addition to the 'random' network, Starim contains some basic sexual networks for modeling sexually transmitted infections (STIs) and a 'maternal' network for modeling mother-to-child transmission.  \n",
-    "\n",
-    "In this tutorial, we'll demonstrate a couple of the built-in networks, and give an overview of how you might create your own for your own disease modeling work."
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "## Overview of Starsim's network structure\n",
-    "\n",
-    "Starsim's networks are stored in the `networks.py` file. At the top of this file you'll see the `Network` class, which provides the template for all the other networks. At a high level, all you need to know is that network stores various *parameters* that determine how pairs of people form and dissolve, and a collection of all the *contacts* that exist at a given point in time. Let's look at a simple example of the `random` network:"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": 7,
-   "metadata": {
-    "tags": []
-   },
-   "outputs": [],
-   "source": [
-    "import starsim as ss\n",
-    "\n",
-    "# Make random network\n",
-    "random = ss.RandomNet()"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {
-    "collapsed": false
-   },
-   "source": [
-    "The network we just made stores two very important things: the list of contacts (found in `random.contacts`) and the parameters (`random.pars`).  \n",
-    "\n",
-    "If we just make a network on its own like this, it's not very interesting. Like all components of Starsim, it needs to be initialized within a sim before the database of contacts gets populated:"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "network = ss.RandomNet()\n",
-    "sir = ss.SIR()\n",
-    "sim = ss.Sim(n_agents=2_000, diseases=sir, networks=network)\n",
-    "sim.init()"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {
-    "collapsed": false
-   },
-   "source": [
-    "Now `sim.networks.random.contacts` has all the information about which agents (referenced by their UID) are paired with which other agents.\n",
-    "\n",
-    "## Sexual networks\n",
-    "Starsim contains a few basic sexual networks, including an `mf` network for modeling heterosexual partnerships and an `msm` network for men having sex with men. Modeling STIs is a straightforward extension of the SIR example we looked at previously. Here's an example modeling syphilis:"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": 9,
-   "metadata": {
-    "collapsed": false
-   },
-   "outputs": [],
-   "source": [
-    "# Create the disease\n",
-    "syph = ss.Syphilis(\n",
-    "    beta = {'mf': [0.25, 0.15]},    # See notes below \n",
-    ")\n",
-    "\n",
-    "# Create a heterosexual network\n",
-    "mf = ss.MFNet(\n",
-    "    duration=1/24,  # Using very short relationships so that we get more transmission \n",
-    "    acts=80,        # Distribution of the number of coital acts/year\n",
-    ")\n",
-    "\n",
-    "pars = dict(start=2000, dur=20, dt=1/12)\n",
-    "sim = ss.Sim(pars=pars, diseases=syph, networks=mf)"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {
-    "collapsed": false
-   },
-   "source": [
-    "The main difference here is that we've added the line `syph.pars['beta'] = {'mf': [0.25, 0.15]}`. Because STIs are generally more transmissible from males to females than vice versa, we can specify `beta` as a list, and we can also allow it to vary between networks. For instance, if we also had an MSM network here, we could set \n",
-    "```\n",
-    "syph.pars.update(beta = {'mf': [0.25, 0.15], 'msm': [0.4, 0.4]})\n",
-    "```\n",
-    "The sexual network itself is structurally similar to the `random` network that we looked at previously, and can be customized by changing the `pars` dictionary that's used to construct it.\n",
-    " \n",
-    "## Maternal networks\n",
-    "Maternal networks are useful for capturing mother-to-child transmission. If you want to use them, you also need to model pregnancy. Here's an example showing how you would do this for syphilis:"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": 10,
-   "metadata": {
-    "collapsed": false
-   },
-   "outputs": [],
-   "source": [
-    "syph = ss.Syphilis(\n",
-    "    beta = {'mf': [0.25, 0.15], 'maternal': [0.99, 0]}\n",
-    ")\n",
-    "\n",
-    "# Make demographic modules\n",
-    "pregnancy = ss.Pregnancy(fertility_rate=20)\n",
-    "death = ss.Deaths(death_rate=15)\n",
-    "\n",
-    "# Make maternal network\n",
-    "maternal = ss.MaternalNet()\n",
-    "\n",
-    "# Make sim\n",
-    "sim = ss.Sim(pars=pars, diseases=syph, networks=[mf, maternal], demographics=[pregnancy, death])"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {
-    "collapsed": false
-   },
-   "source": [
-    "## Making your own network\n",
-    "If you want to make your own network, a good place to start is by inheriting from one of the templates in networks.py. For example, if you wanted to make an `mf` network that has more age structure, you could do something like this:"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": 11,
-   "metadata": {
-    "collapsed": false
-   },
-   "outputs": [],
-   "source": [
-    "class age_mf(ss.MFNet):\n",
-    "    def add_pairs(self, people, ti=None):\n",
-    "        # Custom code to select pairs based on age\n",
-    "        return "
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {
-    "collapsed": false
-   },
-   "source": [
-    "This network would operate exactly like the existing `mf` network, but you could control how new pairs were added.\n"
-   ]
-  }
- ],
- "metadata": {
-  "kernelspec": {
-   "display_name": "Python 3 (ipykernel)",
-   "language": "python",
-   "name": "python3"
-  },
-  "language_info": {
-   "codemirror_mode": {
-    "name": "ipython",
-    "version": 3
-   },
-   "file_extension": ".py",
-   "mimetype": "text/x-python",
-   "name": "python",
-   "nbconvert_exporter": "python",
-   "pygments_lexer": "ipython3",
-   "version": "3.11.7"
-  },
-  "pycharm": {
-   "stem_cell": {
-    "cell_type": "raw",
-    "metadata": {
-     "collapsed": false
-    },
-    "source": []
-   }
-  },
-  "toc": {
-   "base_numbering": 1,
-   "nav_menu": {},
-   "number_sections": true,
-   "sideBar": true,
-   "skip_h1_title": false,
-   "title_cell": "Table of Contents",
-   "title_sidebar": "Contents",
-   "toc_cell": false,
-   "toc_position": {},
-   "toc_section_display": true,
-   "toc_window_display": false
-  }
- },
- "nbformat": 4,
- "nbformat_minor": 4
-}
diff --git a/docs/tutorials/tut_transmission.ipynb b/docs/tutorials/tut_transmission.ipynb
new file mode 100644
index 00000000..7d36b65f
--- /dev/null
+++ b/docs/tutorials/tut_transmission.ipynb
@@ -0,0 +1,482 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# T5 - Transmission\n",
+    "\n",
+    "Many infectious diseases transmit when people come into contact with one another. In the basic SIR example that we covered in the first two tutorials, we used a 'random' contact network, which creates random pairings of people. Generally speaking, the structure of network that you use depends on what disease you're modeling.\n",
+    "\n",
+    "Starsim contains two different approaches to disease transmission.\n",
+    "\n",
+    "**Contact network:** Disease transmits on individual person-to-person relationships. These are networks in the sense that each agent (person) can be viewed as a node and each contact is an edge. Networks are undirected in the sense that an edge will allow transmission in both directions, however you can make transmission stronger one way (p1-->p2) than the other (p2-->p1). The framework includes several pre-built contact network algorithms, each of which manages a list of edges that can change on each simulation step. In addition to the 'random' network, the framework include algorithms representing basic sexual networks for modeling sexually transmitted infections (STIs) and a 'maternal' network for modeling mother-to-child transmission.\n",
+    "\n",
+    "**Mixing pools:** Simulate well-mixed transmission from one group of agents to another. Each mixing pool simulates directed transmission from a group of source agents to a group of destination agents. Instead of simulating individual edges like in a contact network, mixing pools first compute the total \"contagion\" that is shed by the source group. This total contagion value is normalized by the number of agents in the source group to arrive at a value that represents a typical infectious source agent. Then, each destination agent is exposed to this typical (average) infectious agent, likely resulting in new infections in the destination group. Multiple mixing pools can be run simultaneously to simulate heterogeneous transmission between various source and destination groups. Mixing pools could be used with age contact matrices like those developed by Prem et al.\n",
+    "\n",
+    "Note that contact networks and mixing pools can be used together in a single simulation.\n",
+    "\n",
+    "In this tutorial, we'll demonstrate a couple of the built-in contact networks, and give an overview of how you might create your own for your own disease modeling work. We'll also demonstrate how to use mixing pools."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Overview of Starsim's contact network structure\n",
+    "\n",
+    "Starsim's networks are stored in the `networks.py` file. At the top of this file you'll see the `Network` class, which provides the template for all the other networks. At a high level, all you need to know is that network stores various *parameters* that determine how pairs of people form and dissolve, and a collection of all the *edges* that exist at a given point in time. Let's look at a simple example of the `random` network:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 12,
+   "metadata": {
+    "tags": []
+   },
+   "outputs": [],
+   "source": [
+    "import starsim as ss\n",
+    "import numpy as np\n",
+    "import matplotlib.pyplot as plt\n",
+    "\n",
+    "# Make random network\n",
+    "random = ss.RandomNet()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {
+    "collapsed": false
+   },
+   "source": [
+    "The network we just made stores two very important things: the list of contacts (found in `random.edges`) and the parameters (`random.pars`).  \n",
+    "\n",
+    "If we just make a network on its own like this, it's not very interesting. Like all components of Starsim, it needs to be initialized within a `sim` before the database of contacts gets populated:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "net = ss.RandomNet()\n",
+    "sir = ss.SIR()\n",
+    "sim = ss.Sim(n_agents=100, diseases=sir, networks=net) # Decreasing n_agents so we can plot the network below\n",
+    "sim.init()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {
+    "collapsed": false
+   },
+   "source": [
+    "We can easily take a look at the list of edges between agents. All networks are stored in `sim.networks`, and individual networks can be found at their class name after conversion to lower case. So for this example, `sim.networks.randomnet.edges` has all the information about which agents (referenced by their UID) are paired with which other agents as well as a per-edge transmission multiplier called `beta` and a per-edge duration, `dur`."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "sim.networks.randomnet"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Networks provide some basic functionality, for example to identify contacts. Additionally, the edge list can be quickly converted to a Pandas data frame for further analysis and manipulation."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "net = sim.networks.randomnet\n",
+    "\n",
+    "# Use the \"find_contacts\" function to find all contacts of agent with uid 0\n",
+    "contacts_of_agent_0 = net.find_contacts([0])\n",
+    "print('The following agents were identified as contacts of agent 0:', contacts_of_agent_0)\n",
+    "\n",
+    "# Convert to a pandas DataFrame and manually find all edges involving agent with uid 0\n",
+    "df = net.to_df()\n",
+    "uid0_edges = df.loc[(df['p1']==0) | (df['p2']==0)]\n",
+    "print('\\nEdges involving UID 0 extracted from the edge list data frame:')\n",
+    "display(uid0_edges)\n",
+    "\n",
+    "import networkx as nx\n",
+    "G = net.to_graph()\n",
+    "nx.draw(G, with_labels=True)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "\n",
+    "## Sexual networks\n",
+    "Starsim contains a few basic sexual networks, including an `mf` network for modeling heterosexual partnerships and an `msm` network for men having sex with men. Modeling STIs is a straightforward extension of the SIR example we looked at previously. Here's an example modeling syphilis:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 16,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [],
+   "source": [
+    "# Create the disease\n",
+    "syph = ss.Syphilis(\n",
+    "    beta = {'mf': [0.25, 0.15]},    # See notes below \n",
+    ")\n",
+    "\n",
+    "# Create a heterosexual network\n",
+    "mf = ss.MFNet(\n",
+    "    duration=1/24,  # Using very short relationships so that we get more transmission \n",
+    "    acts=80,        # Distribution of the number of coital acts/year\n",
+    ")\n",
+    "\n",
+    "pars = dict(start=2000, dur=20, dt=1/12)\n",
+    "sim = ss.Sim(pars=pars, diseases=syph, networks=mf)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {
+    "collapsed": false
+   },
+   "source": [
+    "The main difference here is that we've added the line `syph.pars['beta'] = {'mf': [0.25, 0.15]}`. Because STIs are generally more transmissible from males to females than vice versa, we can specify `beta` as a list, and we can also allow it to vary between networks. For instance, if we also had an MSM network here, we could set \n",
+    "```\n",
+    "syph.pars.update(beta = {'mf': [0.25, 0.15], 'msm': [0.4, 0.4]})\n",
+    "```\n",
+    "The sexual network itself is structurally similar to the `random` network that we looked at previously, and can be customized by changing the `pars` dictionary that's used to construct it.\n",
+    " \n",
+    "## Maternal networks\n",
+    "Maternal networks are useful for capturing mother-to-child transmission. If you want to use them, you also need to model pregnancy. As new agents are born via pregnancy, edges are automatically added between mother and prenatal child in the maternal network. Here's an example showing how you would do this for syphilis:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 17,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [],
+   "source": [
+    "syph = ss.Syphilis(\n",
+    "    beta = {'mf': [0.25, 0.15], 'maternal': [0.99, 0]}\n",
+    ")\n",
+    "\n",
+    "# Make demographic modules\n",
+    "pregnancy = ss.Pregnancy(fertility_rate=20)\n",
+    "death = ss.Deaths(death_rate=15)\n",
+    "\n",
+    "# Make maternal network\n",
+    "maternal = ss.MaternalNet()\n",
+    "\n",
+    "# Make sim\n",
+    "sim = ss.Sim(pars=pars, diseases=syph, networks=[mf, maternal], demographics=[pregnancy, death])"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {
+    "collapsed": false
+   },
+   "source": [
+    "## Making your own contact network\n",
+    "If you want to make your own network, a good place to start is by inheriting from one of the templates in networks.py. For example, if you wanted to make an `mf` network that has more age structure, you could do something like this:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 18,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [],
+   "source": [
+    "class age_mf(ss.MFNet):\n",
+    "    def add_pairs(self, people, ti=None):\n",
+    "        # Custom code to select pairs based on age\n",
+    "        return "
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {
+    "collapsed": false
+   },
+   "source": [
+    "This network would operate exactly like the existing `mf` network, but you could control how new pairs were added.\n"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Mixing Pools"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Let's compare a simple SIR model configured with transmission using a single Mixing Pool to a similar model using a contact network. Mixing pools are not the same as contact networks, so the results are not expected to be identical."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Transmission using a mixing pool\n",
+    "mp = ss.MixingPool(beta=ss.beta(0.1), contacts=ss.poisson(lam=3)) # A single mixing pool\n",
+    "sir = ss.SIR() # Beta doesn't matter, that's just for networks\n",
+    "sim1 = ss.Sim(diseases=sir, networks=mp, verbose=0, label='Mixing Pool')\n",
+    "\n",
+    "# Transmission using a network\n",
+    "net = ss.RandomNet(n_contacts=ss.poisson(lam=3))\n",
+    "sir = ss.SIR(beta=ss.beta(0.1))\n",
+    "sim2 = ss.Sim(diseases=sir, networks=net, verbose=0, label='Contact Network')\n",
+    "\n",
+    "msim = ss.MultiSim([sim1, sim2]).run()\n",
+    "msim.plot();"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Each MixingPool has several configuration parameters, including `src` and `dst`. The value of these parameters can either be an array of uids, like `ss.uids([0,1,2])` or a callable function that takes `sim` as an argument and returns an array of uids (or a BoolArr) like `lambda sim: sim.people.age<10`. One helper class we have created is called `AgeGroup`, which takes `low` and `high` ages as parameters. This class has some performance advantages when using multiple pools via `MixingPools` through caching.\n",
+    "\n",
+    "While mixing pools do not represent individual agent-to-agent connections, each agent does have an \"effective\" number of contacts, typically drawn from a Poisson distribution. Individuals retain their number of contacts within each pool for the duration of the simulation. Effective contacts are used on the acquisition sides. On the transmission side, the \"average\" infection level, called `trans`, against which susceptible agents are challenged is computed as:\n",
+    "```python\n",
+    "trans = np.mean(disease.infectious[self.src_uids] * disease.rel_trans[self.src_uids])\n",
+    "```\n",
+    "And on the acquisition side, each susceptible agent can be thought of as connected to `eff_contacts` average infectious agents, computed as:\n",
+    "```python\n",
+    "acq = self.eff_contacts[self.dst_uids] * disease.susceptible[self.dst_uids] * disease.rel_sus[self.dst_uids]\n",
+    "```\n",
+    "Ultimately, the acquisition probability for each susceptible agent in `dst_uids` is computed as\n",
+    "```python\n",
+    "p = self.pars.beta * trans * acq\n",
+    "```\n",
+    "\n",
+    "You can see the `beta` parameter here, which is another pool parameter and should use `ss.beta()` to automatically adjust for time units. Mixing pools do not currently use the `beta` parameters configured in individual disease modules."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "mp = ss.MixingPool(\n",
+    "    diseases = 'sir', # Use this pool only for SIR, not other diseases (HIV in this example)\n",
+    "    src = lambda sim: sim.people.age < 15, # Infection sources will be agents under age 15\n",
+    "    dst = ss.AgeGroup(low=15, high=None), # Here demonstrating the AgeGroup class, identical to lambda sim: sim.people.age >= 15\n",
+    "    contacts = ss.poisson(lam=2), # Contacts will be poisson distributed with a mean of 2\n",
+    "    beta = ss.beta(0.2), # And we set beta to 0.2, automatically adjusting for time units via the ss.beta wrapper\n",
+    ")\n",
+    "sim = ss.Sim(diseases=['sir', 'hiv'], networks=mp) # Two diseases, only \"sir\" will have transmission\n",
+    "sim.run()\n",
+    "sim.plot();"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "`MixingPools` (plural) allows you to configure multiple individual `MixingPool` pools. This first example shows how to do age mixing. Like `MixingPool`, `MixingPools` also takes `src` and `dst`, however each of these parameters must be a dictionary with key specifying the name of each group and value as in the singlular mixing pool proving the uids by a callable function or list of uids."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "lows = np.arange(0, 80, 5)\n",
+    "groups = ss.ndict([ss.AgeGroup(low=low, high=low+5) for low in lows[:-1]])\n",
+    "groups.append(ss.AgeGroup(low=lows[-1], high=None)) # 80+\n",
+    "n_groups = len(groups)\n",
+    "\n",
+    "# Let's create a random contact matrix. If using contact matrices from Prem et al., put data here.\n",
+    "cm = np.random.random((n_groups, n_groups)) + 3*np.diag(np.random.rand(n_groups))\n",
+    "\n",
+    "print('The first axis (rows) of the contact matrix is the SOURCE group')\n",
+    "print('The second axis (columns) of the contact matrix is the DESTINATION group')\n",
+    "plt.imshow(cm)\n",
+    "\n",
+    "mps = ss.MixingPools(\n",
+    "    contacts = cm,\n",
+    "    beta = ss.beta(0.2),\n",
+    "    src = groups,\n",
+    "    dst = groups,\n",
+    ")\n",
+    "\n",
+    "# Track new infections by age using an Analyzer, you'll learn about these in a subsequent tutorial\n",
+    "class New_Infections_By_Age(ss.Analyzer):\n",
+    "    def __init__(self, bins, **kwargs):\n",
+    "        super().__init__(**kwargs)\n",
+    "        self.bins = bins\n",
+    "\n",
+    "    def init_post(self):\n",
+    "        super().init_post()\n",
+    "        self.new_cases = np.zeros((len(self), len(self.bins)-1))\n",
+    "\n",
+    "    def step(self):\n",
+    "        new_inf = self.sim.diseases.sir.ti_infected == self.ti\n",
+    "        if not new_inf.any(): return\n",
+    "        self.new_cases[self.ti, :] = np.histogram(self.sim.people.age[new_inf], bins=self.bins)[0]\n",
+    "\n",
+    "    def plot(self):\n",
+    "        fig, ax = plt.subplots()\n",
+    "        ax.set_prop_cycle('color', plt.cm.nipy_spectral(np.linspace(0, 1, len(self.bins))))\n",
+    "\n",
+    "        for i, (b1, b2) in enumerate(zip(self.bins[:-1], self.bins[1:])):\n",
+    "            ax.plot(self.timevec, self.new_cases[:,i], label=f'Age {b1}-{b2}')\n",
+    "        ax.legend(); ax.set_xlabel('Year'); ax.set_ylabel('New Infections')\n",
+    "\n",
+    "az = New_Infections_By_Age(np.concatenate([lows, [1000]]))\n",
+    "\n",
+    "sir = ss.SIR()\n",
+    "sim = ss.Sim(diseases=sir, networks=mps, analyzers=az, dur=5, dt=1/4, n_agents=2000, verbose=0)\n",
+    "sim.run()\n",
+    "sim.analyzers[0].plot()\n",
+    "print('^^^ The default age distribution is uniform, so the 75+ group contains many more people than the other age groups, explaining the larger number of new infections.')"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Mixing pools work well with properties other than age. Here's an example that mixes by socio-economic status (SES)."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from enum import IntEnum\n",
+    "class SES(IntEnum):\n",
+    "    LOW = 0\n",
+    "    MID = 1\n",
+    "    HIGH = 2\n",
+    "\n",
+    "# Create 5,000 people adding 50% LOW, 30% MID, and 20% HIGH SES\n",
+    "ses = ss.FloatArr('ses', default=ss.choice(a=[SES.LOW, SES.MID, SES.HIGH], p=[0.5, 0.3, 0.2]))\n",
+    "ppl = ss.People(n_agents=5_000, extra_states=ses)\n",
+    "\n",
+    "# Create the mixing pools, note that the contact matrix need not be symmetric\n",
+    "mps = ss.MixingPools(\n",
+    "    src = {s.name: lambda sim, s=s: ss.uids(sim.people.ses == s) for s in [SES.LOW, SES.MID, SES.HIGH]},\n",
+    "    dst = {s.name: lambda sim, s=s: ss.uids(sim.people.ses == s) for s in [SES.LOW, SES.MID]}, # People in the HIGH group have not route to acquire new infections, for demonstration purposes here\n",
+    "\n",
+    "    # src on rows (1st dimension), dst on cols (2nd dimension)\n",
+    "    contacts = np.array([\n",
+    "        [2.50, 0.00], # LOW->LOW,  LOW->MID\n",
+    "        [0.05, 1.75], # MID->LOW,  MID->MID\n",
+    "        [0.00, 0.15], # HIGH->LOW, HIGH->MID\n",
+    "    ]),\n",
+    "\n",
+    "    beta = ss.beta(0.2),\n",
+    ")\n",
+    "\n",
+    "# Create an analyzer to track results by SES group\n",
+    "class New_Infections_By_SES(ss.Analyzer):\n",
+    "    def init_results(self):\n",
+    "        self.new_cases = np.zeros((len(self), len(SES)))\n",
+    "\n",
+    "    def step(self):\n",
+    "        new_inf = self.sim.diseases.sir.ti_infected == self.ti\n",
+    "        if not new_inf.any():\n",
+    "            return\n",
+    "\n",
+    "        for ses in SES:\n",
+    "            self.new_cases[self.ti, ses] = np.count_nonzero(new_inf & (self.sim.people.ses==ses))\n",
+    "\n",
+    "az = New_Infections_By_SES()\n",
+    "\n",
+    "def seeding(self, sim, uids):\n",
+    "    p = np.zeros(len(uids))\n",
+    "    high_ses = ss.uids(sim.people.ses == SES.HIGH)\n",
+    "    p[high_ses] = 0.1 # Seed infections into 10% of SES HIGH\n",
+    "    return p\n",
+    "sir = ss.SIR(init_prev = ss.bernoulli(p=seeding)) # The probability of being selected as a seed is determined on a per-agent basis by the above 'seeding' function\n",
+    "\n",
+    "sim = ss.Sim(people=ppl, diseases=sir, networks=mps, analyzers=az, dt=1/12, dur=35, verbose=0)\n",
+    "sim.run()\n",
+    "\n",
+    "\n",
+    "fig, ax = plt.subplots()\n",
+    "new_cases = sim.analyzers[0].new_cases\n",
+    "for ses in [SES.LOW, SES.MID, SES.HIGH]:\n",
+    "    ax.plot(sim.results.timevec, new_cases[:,ses], label=ses.name)\n",
+    "ax.legend(); ax.set_xlabel('Year'); ax.set_ylabel('New Infections');"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "In the above figure, infections are seeded into the HIGH group at the beginning of the simulation, explaining the spike in new cases at that time. HIGH SES can only infect MID SES, so the MID (orange) wave follows. LOW SES can only become infected by MID SES in this example, so the blue wave follows. This example intentionally demonstrates the directional nature of mixing pools and includes strong mixing within groups and bottlenecks between groups to make the point."
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "base",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.12.2"
+  },
+  "pycharm": {
+   "stem_cell": {
+    "cell_type": "raw",
+    "metadata": {
+     "collapsed": false
+    },
+    "source": []
+   }
+  },
+  "toc": {
+   "base_numbering": 1,
+   "nav_menu": {},
+   "number_sections": true,
+   "sideBar": true,
+   "skip_h1_title": false,
+   "title_cell": "Table of Contents",
+   "title_sidebar": "Contents",
+   "toc_cell": false,
+   "toc_position": {},
+   "toc_section_display": true,
+   "toc_window_display": false
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 4
+}
diff --git a/pyproject.toml b/pyproject.toml
index 2bf7067d..c5529abc 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -34,12 +34,12 @@ classifiers = [
 ]
 
 dependencies = [
-    "numpy",
-    "scipy",
+    "numpy>=1.24.0",
     "pandas>=2.0.0",
     "sciris>=3.2.0",
-    "matplotlib",
     "numba",
+    "scipy",
+    "matplotlib",
     "networkx",
     "optuna"
 ]
diff --git a/starsim/__init__.py b/starsim/__init__.py
index c8e5ff8a..6d055bcf 100644
--- a/starsim/__init__.py
+++ b/starsim/__init__.py
@@ -1,8 +1,19 @@
+"""
+Import all Starsim modules
+"""
+
+# Start imports: version and settings
 from .version import __version__, __versiondate__, __license__
-from .settings      import *
+from .settings import dtypes, options
+
+# Optionally print the license
+if options.license:
+    print(__license__)
+
+# Finish imports
 from .utils         import *
 from .arrays        import *
-from .timepars      import *
+from .time          import *
 from .parameters    import *
 from .distributions import *
 from .people        import *
@@ -25,12 +36,8 @@
 import sciris as sc
 root = sc.thispath(__file__).parent
 
-# Import the version and print the license
-if options.verbose:
-    print(__license__)
-
-# Double-check key requirements -- should match setup.py
-reqs = ['sciris>=3.2.0', 'pandas>=2.0.0', 'scipy', 'numba', 'networkx']
+# Double-check key requirements -- should match pyproject.toml
+reqs = ['sciris>=3.2.0', 'pandas>=2.0.0']
 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
+del sc, reqs, msg # Don't keep this in the module
diff --git a/starsim/calibration.py b/starsim/calibration.py
index 5e17cb56..3532c114 100644
--- a/starsim/calibration.py
+++ b/starsim/calibration.py
@@ -135,7 +135,7 @@ class Calibration(sc.prettyobj): # pragma: no cover
     """
     def __init__(self, sim, data, calib_pars, n_trials=None, n_workers=None, total_trials=None, reseed=True,
                  weights=None, fit_args=None, sep='.', name=None, db_name=None, keep_db=None, storage=None,
-                 rand_seed=None, sampler=None, label=None, die=False, debug=False, verbose=True):
+                 rand_seed=None, sampler=None, label=None, die=False, debug=False, verbose=True, save_results=False):
 
         # Handle run arguments
         if n_trials  is None: n_trials  = 20
@@ -159,6 +159,7 @@ def __init__(self, sim, data, calib_pars, n_trials=None, n_workers=None, total_t
         self.fit_args   = sc.mergedicts(fit_args)
         self.die        = die
         self.verbose    = verbose
+        self.save_results = save_results
         self.calibrated = False
         self.before_sim = None
         self.after_sim  = None
@@ -272,7 +273,7 @@ def sim_to_df(sim): # TODO: remove this method
         df_res['time'] = np.floor(np.round(df_res.index, 1)).astype(int)
         return df_res
 
-    def run_trial(self, trial, save=False):
+    def run_trial(self, trial):
         """ Define the objective for Optuna """
         if self.calib_pars is not None:
             calib_pars = self.trial_to_sim_pars(self.calib_pars, trial)
@@ -289,7 +290,7 @@ def run_trial(self, trial, save=False):
         sim_results = sc.objdict()
 
         for skey in self.sim_result_list:
-            if 'prevalence' in skey:
+            if 'prevalence' in skey or skey.startswith('n_'):
                 model_output = df_res.groupby(by='time')[skey].mean()
             else:
                 model_output = df_res.groupby(by='time')[skey].sum()
@@ -297,7 +298,7 @@ def run_trial(self, trial, save=False):
 
         sim_results['time'] = model_output.index.values
         # Store results in temporary files
-        if save:
+        if self.save_results:
             filename = self.tmp_filename % trial.number
             sc.save(filename, sim_results)
 
@@ -313,7 +314,7 @@ def compute_fit(self, sim=None, df_res=None):
         if df_res is None:
             df_res = self.sim_to_df(sim)
         for skey in self.sim_result_list:
-            if 'prevalence' in skey:
+            if 'prevalence' in skey or skey.startswith('n_'):
                 model_output = df_res.groupby(by='time')[skey].mean()
             else:
                 model_output = df_res.groupby(by='time')[skey].sum()
@@ -529,7 +530,7 @@ def plot_sims(self, **kwargs):
             kwargs (dict): passed to MultiSim.plot()
         """
         if self.before_sim is None:
-            self.comfirm_fit()
+            self.confirm_fit()
         msim = ss.MultiSim([self.before_sim, self.after_sim])
         fig = msim.plot(**kwargs)
         return fig
diff --git a/starsim/demographics.py b/starsim/demographics.py
index db0357d2..518f2d09 100644
--- a/starsim/demographics.py
+++ b/starsim/demographics.py
@@ -26,9 +26,10 @@ class Births(Demographics):
     def __init__(self, pars=None, metadata=None, **kwargs):
         super().__init__()
         self.define_pars(
-            birth_rate = 30,
+            unit = 'year',
+            birth_rate = ss.peryear(20),
             rel_birth = 1,
-            units = 1e-3,  # assumes birth rates are per 1000. If using percentages, switch this to 1
+            rate_units = 1e-3,  # assumes birth rates are per 1000. If using percentages, switch this to 1
         )
         self.update_pars(pars, **kwargs)
 
@@ -79,15 +80,20 @@ 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.now) # TODO: make robust to timestep
+            year_ind = sc.findnearest(available_years, sim.t.now('year'))
             nearest_year = available_years[year_ind]
             this_birth_rate = p.birth_rate.loc[nearest_year]
         else:
             this_birth_rate = p.birth_rate
 
-        scaled_birth_prob = this_birth_rate * p.units * p.rel_birth * sim.pars.dt
+        if isinstance(this_birth_rate, ss.TimePar):
+            factor = 1.0
+        else:
+            factor = ss.time_ratio(unit1=self.t.unit, dt1=self.t.dt, unit2='year', dt2=1.0)
+
+        scaled_birth_prob = this_birth_rate * p.rate_units * p.rel_birth * factor
         scaled_birth_prob = np.clip(scaled_birth_prob, a_min=0, a_max=1)
-        n_new = int(np.floor(sim.people.alive.count() * scaled_birth_prob))
+        n_new = int(sc.randround(sim.people.alive.count() * scaled_birth_prob))
         return n_new
 
     def step(self):
@@ -108,10 +114,10 @@ def update_results(self):
         self.results.new[self.ti] = self.n_births
 
         # Calculate crude birth rate (CBR)
-        inv_units = 1.0/self.pars.units
-        births_per_year = self.n_births/self.sim.dt_year
-        denom = self.sim.results.n_alive[self.sim.ti]
-        self.results.cbr[self.ti] = inv_units*births_per_year/denom
+        inv_rate_units = 1.0/self.pars.rate_units
+        births_per_year = self.n_births/self.sim.t.dt_year
+        denom = self.sim.people.alive.sum()
+        self.results.cbr[self.ti] = inv_rate_units*births_per_year/denom
         return
 
     def finalize(self):
@@ -143,7 +149,7 @@ def __init__(self, pars=None, metadata=None, **kwargs):
             pars: dict with arguments including:
                 rel_death: constant used to scale all death rates
                 death_rate: float, dict, or pandas dataframe/series containing mortality data
-                units: units for death rates (see in-line comment on par dict below)
+                rate_units: units for death rates (see in-line comment on par dict below)
 
             metadata: data about the data contained within the data input.
                 "data_cols" is is a dictionary mapping standard keys, like "year" to the
@@ -151,9 +157,10 @@ def __init__(self, pars=None, metadata=None, **kwargs):
         """
         super().__init__()
         self.define_pars(
+            unit = 'year',
             rel_death = 1,
-            death_rate = 20,  # Default = a fixed rate of 2%/year, overwritten if data provided
-            units = 1e-3,  # assumes death rates are per 1000. If using percentages, switch this to 1
+            death_rate = ss.peryear(20),  # Default = a fixed rate of 2%/year, overwritten if data provided
+            rate_units = 1e-3,  # assumes death rates are per 1000. If using percentages, switch this to 1
         )
         self.update_pars(pars, **kwargs)
 
@@ -184,12 +191,12 @@ def standardize_death_data(self):
     @staticmethod # Needs to be static since called externally, although it sure looks like a class method!
     def make_death_prob_fn(self, sim, uids):
         """ Take in the module, sim, and uids, and return the probability of death for each UID on this timestep """
-
         drd = self.death_rate_data
-        if sc.isnumber(drd):
+        if sc.isnumber(drd) or isinstance(drd, ss.TimePar):
             death_rate = drd
-        else:
 
+        # Process data
+        else:
             ppl = sim.people
 
             # Performance optimization - the Deaths module checks for deaths for all agents
@@ -197,7 +204,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.now) # TODO: make work with different timesteps
+            year_ind = sc.findnearest(available_years, sim.t.now('year')) # TODO: make work with different timesteps
             nearest_year = available_years[year_ind]
 
             death_rate = np.empty(uids.shape, dtype=ss_float_)
@@ -215,7 +222,11 @@ def make_death_prob_fn(self, sim, uids):
                 death_rate[:] = s.values[binned_ages]
 
         # Scale from rate to probability. Consider an exponential here.
-        death_prob = death_rate * (self.pars.units * self.pars.rel_death * sim.pars.dt)
+        if isinstance(death_rate, ss.TimePar):
+            factor = self.t.dt # TODO: figure out why this isn't 1.0
+        else:
+            factor = ss.time_ratio(unit1=self.t.unit, dt1=self.t.dt, unit2='year', dt2=1.0)
+        death_prob = death_rate * self.pars.rate_units * self.pars.rel_death * factor
         death_prob = np.clip(death_prob, a_min=0, a_max=1)
 
         return death_prob
@@ -242,9 +253,12 @@ def update_results(self):
 
     def finalize(self):
         super().finalize()
-        n_alive = self.sim.results.n_alive
         self.results.cumulative[:] = np.cumsum(self.results.new)
-        self.results.cmr[:] = 1/self.pars.units*np.divide(self.results.new / self.sim.dt_year, n_alive, where=n_alive>0)
+        units = self.pars.rate_units*self.sim.t.dt_year
+        inds = self.match_time_inds()
+        n_alive = self.sim.results.n_alive[inds]
+        deaths = np.divide(self.results.new, n_alive, where=n_alive>0)
+        self.results.cmr[:] = deaths/units
         return
 
 
@@ -253,6 +267,7 @@ class Pregnancy(Demographics):
     def __init__(self, pars=None, metadata=None, **kwargs):
         super().__init__()
         self.define_pars(
+            unit = 'year',
             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
@@ -262,8 +277,10 @@ def __init__(self, pars=None, metadata=None, **kwargs):
             sex_ratio = ss.bernoulli(0.5), # Ratio of babies born female
             min_age = 15, # Minimum age to become pregnant
             max_age = 50, # Maximum age to become pregnant
-            units = 1e-3, # Assumes fertility rates are per 1000. If using percentages, switch this to 1
+            rate_units = 1e-3, # Assumes fertility rates are per 1000. If using percentages, switch this to 1
             burnin = True, # Should we seed pregnancies that would have happened before the start of the simulation?
+            slot_scale = 5, # Random slots will be assigned to newborn agents between min=n_agents and max=slot_scale*n_agents
+            min_slots  = 100, # Minimum number of slots, useful if the population size is very small
         )
         self.update_pars(pars, **kwargs)
 
@@ -303,10 +320,13 @@ def make_fertility_prob_fn(self, sim, uids):
         frd = self.fertility_rate_data
         fertility_rate = np.zeros(len(sim.people.uid.raw), dtype=ss_float_)
 
+        time_factor = ss.time_ratio(unit1=self.t.unit, dt1=self.t.dt, unit2='year', dt2=1.0)
         if sc.isnumber(frd):
             fertility_rate[uids] = self.fertility_rate_data
+            if isinstance(frd, ss.TimePar):
+                time_factor = 1 # Time conversion performed automatically by TimePar
         else:
-            year_ind = sc.findnearest(frd.index, sim.now-self.pars.dur_pregnancy) # TODO: make time-unit-aware
+            year_ind = sc.findnearest(frd.index, self.t.now('year')-self.pars.dur_pregnancy.to('year')) # TODO: make time-unit-aware
             nearest_year = frd.index[year_ind]
 
             # Assign agents to age bins
@@ -333,7 +353,7 @@ def make_fertility_prob_fn(self, sim, uids):
 
         # Scale from rate to probability
         invalid_age = (age < self.pars.min_age) | (age > self.pars.max_age)
-        fertility_prob = fertility_rate * (self.pars.units * self.pars.rel_fertility * sim.pars.dt)
+        fertility_prob = fertility_rate * (self.pars.rate_units * self.pars.rel_fertility) * time_factor
         fertility_prob[(~self.fecund).uids] = 0 # Currently infecund women cannot become pregnant
         fertility_prob[uids[invalid_age]] = 0 # Women too young or old cannot become pregnant
         fertility_prob = np.clip(fertility_prob[uids], a_min=0, a_max=1)
@@ -362,8 +382,8 @@ def init_pre(self, sim):
         self.pars.p_fertility.set(p=self.make_fertility_prob_fn)
 
         low = sim.pars.n_agents + 1
-        high = int(sim.pars.slot_scale*sim.pars.n_agents)
-        high = np.maximum(high, sim.pars.min_slots) # Make sure there are at least min_slots slots to avoid artifacts related to small populations
+        high = int(self.pars.slot_scale*sim.pars.n_agents)
+        high = np.maximum(high, self.pars.min_slots) # Make sure there are at least min_slots slots to avoid artifacts related to small populations
         self.choose_slots = ss.randint(low=low, high=high, sim=sim, module=self)
         return
 
@@ -382,21 +402,20 @@ def init_results(self):
         return
 
     def step(self):
-        if self.sim.ti == 0 and self.pars.burnin:
+        if self.ti == 0 and self.pars.burnin: # TODO: refactor
             dtis = np.arange(np.ceil(-1 * self.pars.dur_pregnancy), 0, 1).astype(int)
             for dti in dtis:
-                self.sim.ti = dti
+                self.t.ti = dti
                 self.do_step()
-            self.sim.ti = 0
+            self.t.ti = 0
         new_uids = self.do_step()
-
         return new_uids
 
     def do_step(self):
         """ Perform all updates """
         self.update_states()
         conceive_uids = self.make_pregnancies()
-        self.n_pregnancies = len(conceive_uids)
+        self.n_pregnancies += len(conceive_uids) # += to handle burn-in
         new_uids = self.make_embryos(conceive_uids)
         return new_uids
 
@@ -454,6 +473,9 @@ def make_pregnancies(self):
         eligible_uids = self.sim.people.female.uids
         conceive_uids = self.pars.p_fertility.filter(eligible_uids)
 
+        if len(conceive_uids) == 0:
+            return ss.uids()
+
         # Validation
         if np.any(self.pregnant[conceive_uids]):
             which_uids = conceive_uids[self.pregnant[conceive_uids]]
@@ -461,8 +483,7 @@ def make_pregnancies(self):
             raise ValueError(errormsg)
 
         # Set prognoses for the pregnancies
-        if len(conceive_uids) > 0:
-            self.set_prognoses(conceive_uids)
+        self.set_prognoses(conceive_uids)
         return conceive_uids
 
     def make_embryos(self, conceive_uids):
@@ -491,8 +512,8 @@ def make_embryos(self, conceive_uids):
                     start = np.full(n_unborn, fill_value=self.ti)
                     layer.add_pairs(conceive_uids, new_uids, dur=durs, start=start)
 
-        if self.sim.ti < 0:
-            people.age[new_uids] += -self.sim.ti * self.dt # Age to ti=0, TODO: Check!
+        if self.ti < 0:
+            people.age[new_uids] += -self.ti * self.sim.t.dt_year # Age to ti=0
 
         return new_uids
 
@@ -511,7 +532,7 @@ def set_prognoses(self, uids):
         self.ti_pregnant[uids] = ti
 
         # Outcomes for pregnancies
-        dur_preg = np.ones(len(uids))*self.pars.dur_pregnancy  # Duration in years
+        dur_preg = np.ones(len(uids))*self.pars.dur_pregnancy
         dur_postpartum = self.pars.dur_postpartum.rvs(uids)
         dead = self.pars.p_maternal_death.rvs(uids)
         self.ti_delivery[uids] = ti + dur_preg # Currently assumes maternal deaths still result in a live baby
@@ -524,7 +545,7 @@ def set_prognoses(self, uids):
 
     def finish_step(self):
         super().finish_step()
-        death_uids = ss.uids(self.sim.people.ti_dead <= self.sim.ti)
+        death_uids = ss.uids(self.sim.people.ti_dead <= self.ti)
         if len(death_uids) == 0:
             return
 
@@ -553,10 +574,18 @@ def update_results(self):
         ti = self.ti
         self.results['pregnancies'][ti] = self.n_pregnancies
         self.results['births'][ti] = self.n_births
+
+        # Reset for the next step
+        self.n_pregnancies = 0
+        self.n_births = 0
+
         return
 
     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_year, n_alive, where=n_alive>0)
+        units = self.pars.rate_units*self.sim.t.dt_year
+        inds = self.match_time_inds()
+        n_alive = self.sim.results.n_alive[inds]
+        births = np.divide(self.results['births'], n_alive, where=n_alive>0)
+        self.results['cbr'][:] = births/units
         return
diff --git a/starsim/disease.py b/starsim/disease.py
index 378aaeba..e340b34f 100644
--- a/starsim/disease.py
+++ b/starsim/disease.py
@@ -152,7 +152,7 @@ class Infection(Disease):
     operate on to capture co-infection
     """
 
-    def __init__(self, *args, **kwargs):
+    def __init__(self, pars=None, *args, **kwargs):
         super().__init__(*args, **kwargs)
         self.define_states(
             ss.State('susceptible', default=True, label='Susceptible'),
@@ -162,6 +162,11 @@ def __init__(self, *args, **kwargs):
             ss.FloatArr('ti_infected', label='Time of infection' ),
         )
 
+        self.define_pars(
+            init_prev = None, # Replace None with a ss.bernoulli to seed infections
+        )
+        self.update_pars(pars, **kwargs)
+
         # Define random number generator for determining transmission
         self.trans_rng = ss.multi_random('source', 'target')
         return
@@ -296,7 +301,7 @@ def infect(self):
 
         for i, (nkey,net) in enumerate(self.sim.networks.items()):
             nk = ss.standardize_netkey(nkey)
-            if len(net): # Skip networks with no edges
+            if isinstance(net, ss.Network) and len(net): # Skip networks with no edges
                 edges = net.edges
                 p1p2b0 = [edges.p1, edges.p2, betamap[nk][0]] # Person 1, person 2, beta 0
                 p2p1b1 = [edges.p2, edges.p1, betamap[nk][1]] # Person 2, person 1, beta 1
diff --git a/starsim/diseases/cholera.py b/starsim/diseases/cholera.py
index b4202871..ccaf2514 100644
--- a/starsim/diseases/cholera.py
+++ b/starsim/diseases/cholera.py
@@ -124,9 +124,9 @@ def calc_environmental_prev(self):
         r.env_conc[ti] = r.env_prev[ti] / (r.env_prev[ti] + p.half_sat_rate)
         return
 
-    def set_prognoses(self, uids, source_uids=None):
+    def set_prognoses(self, uids, sources=None):
         """ Set prognoses for those who get infected """
-        super().set_prognoses(uids, source_uids)
+        super().set_prognoses(uids, sources)
         ti = self.ti
 
         self.susceptible[uids] = False
diff --git a/starsim/diseases/gonorrhea.py b/starsim/diseases/gonorrhea.py
index 6a9dda27..31c2d8a7 100644
--- a/starsim/diseases/gonorrhea.py
+++ b/starsim/diseases/gonorrhea.py
@@ -58,9 +58,9 @@ def step_state(self):
 
         return
 
-    def set_prognoses(self, uids, source_uids=None):
+    def set_prognoses(self, uids, sources=None):
         """ Natural history of gonorrhea for adult infection """
-        super().set_prognoses(uids, source_uids)
+        super().set_prognoses(uids, sources)
         ti = self.ti
 
         # Set infection status
diff --git a/starsim/diseases/ncd.py b/starsim/diseases/ncd.py
index d55fe67e..5f32867b 100644
--- a/starsim/diseases/ncd.py
+++ b/starsim/diseases/ncd.py
@@ -55,7 +55,8 @@ def step_state(self):
         ti = self.ti
         deaths = (self.ti_dead == ti).uids
         self.sim.people.request_death(deaths)
-        self.log.add_data(deaths, died=True)
+        if self.pars.log:
+            self.log.add_data(deaths, died=True)
         self.results.new_deaths[ti] = len(deaths) # Log deaths attributable to this module
         return
 
@@ -64,7 +65,7 @@ def step(self):
         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.dt)
+        self.ti_dead[new_cases] = ti + sc.randround(prog_years / self.t.dt) # TODO: update to allow non-year units
         super().set_prognoses(new_cases)
         return new_cases
 
diff --git a/starsim/diseases/sir.py b/starsim/diseases/sir.py
index 6243d8d7..0c3ebb94 100644
--- a/starsim/diseases/sir.py
+++ b/starsim/diseases/sir.py
@@ -54,9 +54,10 @@ def step_state(self):
             sim.people.request_death(deaths)
         return
 
-    def set_prognoses(self, uids, source_uids=None):
+    def set_prognoses(self, uids, sources=None):
         """ Set prognoses """
-        ti = self.ti
+        super().set_prognoses(uids, sources)
+        ti = self.t.ti
         self.susceptible[uids] = False
         self.infected[uids] = True
         self.ti_infected[uids] = ti
@@ -86,8 +87,9 @@ def plot(self, **kwargs):
         """ Default plot for SIR model """
         fig = plt.figure()
         kw = sc.mergedicts(dict(lw=2, alpha=0.8), kwargs)
+        res = self.results
         for rkey in ['n_susceptible', 'n_infected', 'n_recovered']:
-            plt.plot(self.timevec, self.results[rkey], label=self.results[rkey].label, **kw)
+            plt.plot(res.timevec, res[rkey], label=res[rkey].label, **kw)
         plt.legend(frameon=False)
         plt.xlabel('Time')
         plt.ylabel('Number of people')
@@ -136,9 +138,9 @@ def update_immunity(self):
         self.rel_sus[has_imm] = np.maximum(0, 1 - self.immunity[has_imm])
         return
 
-    def set_prognoses(self, uids, source_uids=None):
+    def set_prognoses(self, uids, sources=None):
         """ Set prognoses """
-        super().set_prognoses(uids, source_uids)
+        super().set_prognoses(uids, sources)
         self.susceptible[uids] = False
         self.infected[uids] = True
         self.ti_infected[uids] = self.ti
@@ -170,8 +172,9 @@ def plot(self, **kwargs):
         """ Default plot for SIS model """
         fig = plt.figure()
         kw = sc.mergedicts(dict(lw=2, alpha=0.8), kwargs)
+        res = self.results
         for rkey in ['n_susceptible', 'n_infected']:
-            plt.plot(self.timevec, self.results[rkey], label=self.results[rkey].label, **kw)
+            plt.plot(res.timevec, res[rkey], label=res[rkey].label, **kw)
         plt.legend(frameon=False)
         plt.xlabel('Time')
         plt.ylabel('Number of people')
diff --git a/starsim/diseases/syphilis.py b/starsim/diseases/syphilis.py
index 77846408..c7b503cb 100644
--- a/starsim/diseases/syphilis.py
+++ b/starsim/diseases/syphilis.py
@@ -4,8 +4,8 @@
 
 import numpy as np
 import sciris as sc
-from sciris import randround as rr # Since used frequently
 import starsim as ss
+rr = sc.randround # Since used frequently
 
 
 __all__ = ['Syphilis']
@@ -248,14 +248,14 @@ def set_latent_long_prognoses(self, uids):
 
         return
 
-    def set_congenital(self, uids, source_uids=None):
+    def set_congenital(self, uids, sources=None):
         """ Natural history of syphilis for congenital infection """
         sim = self.sim
 
         # Determine outcomes
         for state in ['active', 'latent']:
 
-            source_state_inds = getattr(self, state)[source_uids].nonzero()[0]
+            source_state_inds = getattr(self, state)[sources].nonzero()[0]
             state_uids = uids[source_state_inds]
 
             if len(state_uids) > 0:
@@ -266,7 +266,7 @@ 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
+                ratio = ss.time_ratio(unit1='year', dt1=1.0, unit2=self.t.unit, dt2=self.t.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:
diff --git a/starsim/distributions.py b/starsim/distributions.py
index 42bb243f..759fd0e4 100644
--- a/starsim/distributions.py
+++ b/starsim/distributions.py
@@ -11,18 +11,17 @@
 __all__ = ['link_dists', 'make_dist', 'dist_list', 'Dists', 'Dist']
 
 
-def str2int(string, modulo=1_000_000):
+def str2int(string, modulo=1_000_000_000):
     """
-    Convert a string to an int
+    Convert a string to an int to use as a random seed; not for the user
 
     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
+    (sc.sha) is 5x slower than int.from_bytes(string.encode(), byteorder='big'), but should
     only add a couple milliseconds to a typical sim.
     """
-    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
+    integer = sc.sha(string, asint=True) # Hash the string to an integer
+    seed = integer % modulo # Ensure a user-friendly representation user-friendly
+    return seed
 
 
 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
@@ -85,7 +84,10 @@ def init(self, obj=None, base_seed=None, sim=None, force=False):
             raise ValueError(errormsg)
 
         # Do not look for distributions in the people states, since they shadow the "real" states
-        skip = id(sim.people._states) if sim is not None else None
+        skip = dict(
+            ids=id(sim.people._states) if sim is not None else None,
+            keys='module',
+        )
 
         # Find and initialize the distributions
         self.dists = sc.search(obj, type=Dist, skip=skip, flatten=True)
@@ -101,7 +103,7 @@ def init(self, obj=None, base_seed=None, sim=None, force=False):
     def check_seeds(self):
         """ Check that no two distributions share the same seed """
         checked = dict()
-        for trace,dist in self.dists.items():
+        for dist in self.dists.values():
             seed = dist.seed
             if seed in checked.keys():
                 raise DistSeedRepeatError(checked[seed], dist)
@@ -359,7 +361,7 @@ def jump_dt(self, ti=None, force=False):
             ti (int): if specified, jump to this timestep (default: current module timestep plus one)
         """
         if ti is None:
-            ti = self.module.ti + 1
+            ti = self.module.t.ti + 1
         to = self.dt_jump_size*ti
         return self.jump(to=to, force=force)
 
@@ -1017,7 +1019,7 @@ def filter(self, uids=None, both=False):
         """ Return UIDs that correspond to True, or optionally return both True and False """
         if uids is None:
             uids = self.sim.people.auids # All active UIDs
-        elif isinstance(uids, ss.BoolArr):
+        elif isinstance(uids, (ss.BoolArr, ss.IndexArr)):
             uids = uids.uids
 
         bools = self.rvs(uids)
diff --git a/starsim/interventions.py b/starsim/interventions.py
index 0fae1f8b..59815155 100644
--- a/starsim/interventions.py
+++ b/starsim/interventions.py
@@ -21,6 +21,17 @@ def __init__(self, *args, eligibility=None, **kwargs):
         self.eligibility = eligibility
         return
 
+    @property
+    def has_product(self):
+        """ Check if the intervention has a product """
+        return hasattr(self, 'product') and self.product is not None
+
+    def init_pre(self, sim):
+        super().init_pre(sim)
+        if self.has_product:
+            self.product.init_pre(self.sim)
+        return
+
     def _parse_product(self, product):
         """
         Parse the product input
@@ -91,7 +102,8 @@ def init_pre(self, sim):
             self.end_year = self.years[-1]
 
         # More validation
-        if not(any(np.isclose(self.start_year, sim.timevec)) and any(np.isclose(self.end_year, sim.timevec))):
+        yearvec = sim.t.yearvec
+        if not(any(np.isclose(self.start_year, yearvec)) and any(np.isclose(self.end_year, yearvec))):
             errormsg = 'Years must be within simulation start and end dates.'
             raise ValueError(errormsg)
 
@@ -100,11 +112,11 @@ def init_pre(self, sim):
         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.timevec, self.start_year)
-        self.end_point   = sc.findfirst(sim.timevec, self.end_year) + adj_factor
+        self.start_point = sc.findfirst(yearvec, self.start_year)
+        self.end_point   = sc.findfirst(yearvec, 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, dt)
+        self.yearvec     = np.arange(self.start_year, self.end_year + adj_factor, dt) # TODO: integrate with self.t
 
         # Get the probability input into a format compatible with timepoints
         if len(self.years) != len(self.prob):
@@ -464,4 +476,4 @@ class campaign_vx(BaseVaccination, CampaignDelivery):
     Campaign vaccination - an instance of base vaccination combined with campaign delivery.
     See base classes for a description of input arguments.
     """
-    pass
\ No newline at end of file
+    pass
diff --git a/starsim/loop.py b/starsim/loop.py
index 315d7a09..e2753891 100644
--- a/starsim/loop.py
+++ b/starsim/loop.py
@@ -135,12 +135,11 @@ def collect_abs_tvecs(self):
         # Handle the sim and people first
         sim = self.sim
         for key in ['sim', 'people']:
-            self.abs_tvecs[key] = sim.abs_tvec
+            self.abs_tvecs[key] = sim.t.abstvec
 
         # Handle all other modules
         for mod in sim.modules:
-            abs_tvec = ss.make_abs_tvec(mod.timevec, mod.unit, sim.pars.unit)
-            self.abs_tvecs[mod.name] = abs_tvec
+            self.abs_tvecs[mod.name] = mod.t.abstvec
 
         return self.abs_tvecs
 
diff --git a/starsim/modules.py b/starsim/modules.py
index e8039aae..0574fb29 100644
--- a/starsim/modules.py
+++ b/starsim/modules.py
@@ -6,7 +6,9 @@
 import starsim as ss
 from functools import partial
 
-__all__ = ['module_map', 'find_modules', 'Module', 'Analyzer', 'Connector']
+__all__ = ['module_map', 'find_modules', 'Base', 'Module', 'Analyzer', 'Connector']
+
+module_args = ['name', 'label'] # Define allowable module arguments
 
 
 def module_map(key=None):
@@ -44,56 +46,84 @@ def find_modules(key=None, flat=False):
     return modules if key is None else modules[key]
 
 
-class Module(sc.quickobj):
+class Base(sc.quickobj):
+    """
+    The parent class for Sim and Module objects
+    """
+    def __bool__(self):
+        """ Ensure that zero-length modules (e.g. networks) are still truthy """
+        return True
+
+    def __len__(self):
+        """ The length of a module is the number of timepoints; see also len(sim) """
+        try:    return self.t.npts
+        except: return 0
+
+    def disp(self, output=False, **kwargs):
+        """ Display the full object """
+        out = sc.prepr(self, **kwargs)
+        if not output:
+            print(out)
+        else:
+            return out
+
+    @property
+    def ti(self):
+        """ Get the current module timestep """
+        try:    return self.t.ti
+        except: return None
+
+    @property
+    def now(self):
+        """ Shortcut to self.t.now() """
+        try:    return self.t.now()
+        except: return None
+
+    @property
+    def timevec(self):
+        """ Shortcut to self.t.timevec """
+        try:    return self.t.timevec
+        except: return None
+
+
+
+class Module(Base):
     """
     The main base class for all Starsim modules: diseases, networks, interventions, etc.
 
     Args:
         name (str): a short, key-like name for the module (e.g. "randomnet")
         label (str): the full, human-readable name for the module (e.g. "Random network")
-        unit (str): the time unit (e.g. 'day', 'year'); inherits from sim if not supplied
-        dt (float): the timestep (e.g. 1.0, 0.1); inherits from sim if not supplied
+        kwargs (dict): passed to ss.Time() (e.g. start, stop, unit, dt)
     """
 
-    def __init__(self, name=None, label=None, unit=None, dt=None):
+    def __init__(self, name=None, label=None, **kwargs):
         # Handle parameters
         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.t = ss.Time(**kwargs, name=self.name, init=False)
 
         # Properties to be added by init_pre()
         self.sim = None
         self.dists = None # Turned into a Dists object by sim.init_dists() if this module has dists
         self.results = ss.Results(self.name)
 
-        # Time properties, added by init_time_pars()
-        self.unit = None
-        self.dt = None
-        self.timevec = None
-        self.npts = None
-        self.ti = None
-
         # Finish initialization
         self.pre_initialized = False
         self.initialized = False
         self.finalized = False
         return
 
-    def __bool__(self):
-        """ Ensure that zero-length modules (e.g. networks) are still truthy """
-        return True
-
     def __call__(self, *args, **kwargs):
         """ Allow modules to be called like functions """
         return self.step(*args, **kwargs)
 
-    def disp(self, output=False):
-        """ Display the full object """
-        out = sc.prepr(self)
-        if not output:
-            print(out)
-        else:
-            return out
+    def _reconcile(self, key, value=None, default=None):
+        """ Reconcile module attributes, parameters, and input arguments """
+        parval = self.pars.get(key)
+        attrval = getattr(self, key, parval)
+        val = sc.ifelse(value, attrval, default)
+        return val
 
     def set_metadata(self, name=None, label=None):
         """ Set metadata for the module """
@@ -105,14 +135,8 @@ def set_metadata(self, name=None, label=None):
                     raise TypeError(errormsg)
 
         # Set values
-        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')))
+        self.name = self._reconcile('name', name, self.__class__.__name__.lower())
+        self.label = self._reconcile('label', label, self.name)
         return
 
     def define_pars(self, inherit=True, **kwargs): # TODO: think if inherit should default to true or false
@@ -135,14 +159,15 @@ def update_pars(self, pars, **kwargs):
         self.pars.update(matches)
 
         # 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']}
+        metadata = {key:pars.get(key, self.pars.get(key)) for key in module_args}
+        timepars = {key:pars.get(key, self.pars.get(key)) for key in ss.time.time_args}
         self.set_metadata(**metadata)
-        self.set_time_pars(**timepars)
+        self.t.update(**timepars)
 
         # Should be no remaining pars
-        if len(pars):
-            errormsg = f'{len(pars)} unrecognized arguments for {self.name}: {sc.strjoin(pars.keys())}'
+        remaining = set(pars.keys()) - set(module_args + ss.time.time_args)
+        if len(remaining):
+            errormsg = f'{len(pars)} unrecognized arguments for {self.name}: {sc.strjoin(remaining)}'
             raise ValueError(errormsg)
         return
 
@@ -160,14 +185,14 @@ def init_pre(self, sim, force=False):
             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_time() # Initialize the modules' time parameters and link them to the sim
             self.init_results()
             self.pre_initialized = True
         return
 
     def init_results(self):
         """ Initialize any results required; part of init_pre() """
-        self.results.timevec = self.timevec # Store the timevec in the results for plotting
+        self.results.timevec = self.t.timevec # Store the timevec in the results for plotting
         return
 
     def init_post(self):
@@ -178,15 +203,10 @@ def init_post(self):
         self.initialized = True
         return
 
-    def init_time_pars(self, force=False):
+    def init_time(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
+        # Update time and initialize
+        self.t.init(sim=self.sim) # Sets the absolute sim time vector
 
         # Find all time parameters in the module
         timepars = sc.search(self.pars, type=ss.TimePar) # Should it be self or self.pars?
@@ -194,22 +214,19 @@ def init_time_pars(self, force=False):
         # Initialize them with the parent module
         for timepar in timepars.values():
             if force or not timepar.initialized:
-                timepar.init(parent=self, die=False) # In some cases, the values can't be initialized; that's OK here
-
-        # 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
+                timepar.init(parent=self.t, die=False) # In some cases, the values can't be initialized; that's OK here
         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 match_time_inds(self, inds=None):
+         """ Find the nearest matching sim time indices for the current module """
+         if inds is None: inds = Ellipsis
+         self_tvec = self.t.abstvec
+         sim_tvec = self.sim.t.abstvec
+         if len(self_tvec) == len(sim_tvec): # Shortcut to avoid doing matching
+             return inds
+         else:
+             out = sc.findnearest(sim_tvec, [inds])
+             return out
 
     def start_step(self):
         """ Tasks to perform at the beginning of the step """
@@ -219,11 +236,12 @@ def start_step(self):
 
     def step(self):
         """ Define how the module updates over time -- the key part of Starsim!! """
-        pass
+        errormsg = f'Module "{self.name}" does not define a "step" method: use "def step(self): pass" if this is intentional'
+        raise NotImplementedError(errormsg)
 
     def finish_step(self):
         """ Define what should happen at the end of the step; at minimum, increment ti """
-        self.ti += 1
+        self.t.ti += 1
         return
 
     def update_results(self):
@@ -278,7 +296,7 @@ def define_results(self, *args, check=True):
                 result = arg
 
             # Update with module information
-            result.update(module=self.name, shape=self.npts, timevec=self.timevec)
+            result.update(module=self.name, shape=self.t.npts, timevec=self.t.timevec)
 
             # Add the result to the dict of results; does automatic checking
             self.results += result
@@ -346,7 +364,7 @@ def plot(self):
         """ Plot all results in the module """
         with sc.options.with_style('fancy'):
             flat = sc.flattendict(self.results, sep=': ')
-            timevec = self.timevec
+            timevec = self.t.timevec
             fig, axs = sc.getrowscols(len(flat), make=True)
             for ax, (k, v) in zip(axs.flatten(), flat.items()):
                 ax.plot(timevec, v)
diff --git a/starsim/networks.py b/starsim/networks.py
index 80f19b5f..dba7c8ab 100644
--- a/starsim/networks.py
+++ b/starsim/networks.py
@@ -15,11 +15,18 @@
 
 
 # Specify all externally visible functions this file defines; see also more definitions below
-__all__ = ['Network', 'DynamicNetwork', 'SexualNetwork']
+__all__ = ['Route', 'Network', 'DynamicNetwork', 'SexualNetwork']
 
 # %% General network classes
 
-class Network(ss.Module):
+class Route(ss.Module):
+    """
+    A transmission route -- e.g., a network, mixing pool, environmental transmission, etc.
+    """
+    pass
+
+
+class Network(Route):
     """
     A class holding a single network of contact edges (connections) between people
     as well as methods for updating these.
@@ -342,7 +349,7 @@ def step(self):
 
     def end_pairs(self):
         people = self.sim.people
-        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
+        self.edges.dur = self.edges.dur - self.t.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]
@@ -378,7 +385,7 @@ def available(self, people, sex):
 
     def net_beta(self, disease_beta=None, inds=None, disease=None):
         if inds is None: inds = Ellipsis
-        return self.edges.beta[inds] * (1 - (1 - disease_beta) ** (self.edges.acts[inds] * self.dt))
+        return self.edges.beta[inds] * (1 - (1 - disease_beta) ** (self.edges.acts[inds] * self.t.dt))
 
 
 # %% Specific instances of networks
@@ -636,30 +643,29 @@ def init_post(self):
         return
 
     def step(self):
-
         # Motion step
-        vdt = self.pars.v * self.dt
+        vdt = self.pars.v * self.t.dt
         self.x[:] = self.x + vdt * np.cos(self.theta)
         self.y[:] = self.y + vdt * np.sin(self.theta)
 
         # Wall bounce
 
-        # Right edge
+        ## Right edge
         inds = (self.x > 1).uids
         self.x[inds] = 2 - self.x[inds]
         self.theta[inds] = np.pi - self.theta[inds]
 
-        # Left edge
+        ## Left edge
         inds = (self.x < 0).uids
         self.x[inds] =  -self.x[inds]
         self.theta[inds] = np.pi - self.theta[inds]
 
-        # Top edge
+        ## Top edge
         inds = (self.y > 1).uids
         self.y[inds] = 2 - self.y[inds]
         self.theta[inds] = - self.theta[inds]
 
-        # Bottom edge
+        ## Bottom edge
         inds = (self.y < 0).uids
         self.y[inds] = -self.y[inds]
         self.theta[inds] = - self.theta[inds]
@@ -718,6 +724,10 @@ def get_edges(self):
         self.append(dict(p1=indices, p2=indices, beta=np.zeros_like(indices)))
         return
 
+    def step(self):
+        """ Not used for NullNet """
+        pass
+
 
 class MFNet(SexualNetwork):
     """
@@ -805,7 +815,7 @@ def add_pairs(self):
 
     def step(self):
         self.end_pairs()
-        self.set_network_states(upper_age=self.dt) # TODO: check
+        self.set_network_states(upper_age=self.t.dt) # TODO: check
         self.add_pairs()
         return
 
@@ -987,3 +997,265 @@ class PostnatalNet(MaternalNet):
     def __init__(self, key_dict=None, prenatal=False, postnatal=True, **kwargs):
         super().__init__(key_dict=key_dict, prenatal=prenatal, postnatal=postnatal, **kwargs)
         return
+
+
+__all__ += ['AgeGroup', 'MixingPools', 'MixingPool']
+
+class AgeGroup(sc.prettyobj):
+    """ A simple age-based filter that returns uids of agents that match the criteria """
+    def __init__(self, low, high, do_cache=True):
+        self.low = low
+        self.high = high
+
+        self.do_cache = do_cache
+        self.uids = None # Cached
+        self.ti_cache = -1
+
+        self.name = repr(self)
+        return
+
+    def __call__(self, sim):
+        if (not self.do_cache) or (self.ti_cache != sim.ti):
+            in_group = sim.people.age >= self.low
+            if self.high is not None:
+                in_group = in_group & (sim.people.age < self.high)
+            self.uids = ss.uids(in_group)
+            self.ti_cache = sim.ti
+        return self.uids
+
+    def __repr__(self):
+        return f'age({self.low}-{self.high})'
+
+
+class MixingPools(Route):
+    """
+    A container for creating a rectangular array of MixingPool instances
+
+    By default, separates the population into <15 and >15 age groups.
+
+    Args:
+        diseases (str): the diseases that transmit via these mixing pools
+        src (inds): source agents; can be AgeGroup(), ss.uids(), or lambda(sim); None indicates all alive agents
+        dst (inds): destination agents; as above
+        beta (float): overall transmission via these mixing pools
+        contacts (array): the relative connectivity between different mixing pools (can be float or Dist)
+
+    **Example**::
+
+        import starsim as ss
+        mps = ss.MixingPools(
+            diseases = 'sis',
+            beta = 0.1,
+            src = {'0-15': ss.AgeGroup(0, 15), '15+': ss.AgeGroup(15, None)},
+            dst = {'0-15': ss.AgeGroup(0, 15), '15+': ss.AgeGroup(15, None)},
+            contacts = [[2.4, 0.49], [0.91, 0.16]],
+        )
+        sim = ss.Sim(diseases='sis', networks=mps).run()
+        sim.plot()
+    """
+    def __init__(self, pars=None, **kwargs):
+        super().__init__()
+        self.define_pars(
+            diseases = None,
+            src = None,
+            dst = None,
+            beta = ss.beta(0.1),
+            contacts = None,
+        )
+        self.update_pars(pars, **kwargs)
+        self.validate_pars()
+        self.pools = []
+        return
+
+    def __len__(self):
+        try:    return len(self.pools)
+        except: return 0
+
+    def validate_pars(self):
+        """ Check that src and dst have correct types, and contacts is the correct shape """
+        p = self.pars
+
+        # Validate src and dst
+        if not isinstance(p.src, dict):
+            errormsg = f'src must be a provided as a dictionary, not {type(self.pars.src)}'
+            raise TypeError(errormsg)
+        if not isinstance(p.dst, dict):
+            raise TypeError(f'dst must be a provided as a dictionary, not {type(self.pars.src)}')
+        p.src = sc.objdict(p.src)
+        p.dst = sc.objdict(p.dst)
+
+        # Validate the contacts
+        if p.contacts is None:
+            p.contacts = np.ones((len(p.src), len(p.dst)))
+        p.contacts = np.array(p.contacts)
+        actual = p.contacts.shape
+        expected = (len(p.src), len(p.dst))
+        if actual != expected:
+            errormsg = f'The number of source and destination groups must match the number of rows and columns in the mixing matrix, but {actual} != {expected}.'
+            raise ValueError(errormsg)
+
+        return
+
+    def init_pre(self, sim):
+        super().init_pre(sim)
+        p = self.pars
+
+        self.pools = []
+        for i,sk,src in p.src.enumitems():
+            for j,dk,dst in p.dst.enumitems():
+                contacts = p.contacts[i,j]
+                if sc.isnumber(contacts): # If it's a number, convert to a distribution
+                    contacts = ss.poisson(lam=contacts)
+                name = f'pool:{sk}->{dk}'
+                mp = MixingPool(name=name, diseases=p.diseases, beta=p.beta, contacts=contacts, src=src, dst=dst)
+                mp.init_pre(sim) # Initialize the pool
+                self.pools.append(mp)
+        return
+
+    def init_post(self):
+        """ Initialize each mixing pool """
+        for mp in self.pools:
+            mp.init_post()
+        return
+
+    def step(self):
+        """ Step each mixing pool """
+        for mp in self.pools:
+            mp.step()
+        return
+
+    def remove_uids(self, uids):
+        """ Remove UIDs from each mixing pool """
+        for mp in self.pools:
+            mp.remove_uids(uids)
+        return
+
+
+class MixingPool(Route):
+    """
+    Define a single mixing pool; can be used as a drop-in replacement for a network.
+
+    Args:
+        diseases (str): the diseases that transmit via this mixing pool
+        src (inds): source agents; can be AgeGroup(), ss.uids(), or lambda(sim); None indicates all alive agents
+        dst (inds): destination agents; as above
+        beta (float): overall transmission
+        contacts (Dist): the number of effective contacts of the destination agents
+
+    **Example**::
+
+        import starsim as ss
+
+        # Set the parameters
+        mp_pars = dict(
+            src = lambda sim: sim.people.male, # only males are infectious
+            dst = None, # all agents are susceptible
+            beta = ss.beta(0.2),
+            contacts = ss.poisson(lam=4),
+        )
+
+        # Seed 5% of the male population
+        def p_init(self, sim, uids):
+            return 0.05*sim.people.male
+
+        # Create and run the sim
+        sis = ss.SIS(init_prev=p_init)
+        mp = ss.MixingPool(mp_pars)
+        sim = ss.Sim(diseases=sis, networks=mp)
+        sim.run()
+        sim.plot()
+    """
+    def __init__(self, pars=None, **kwargs):
+        super().__init__()
+        self.define_pars(
+            diseases = None,
+            src = None,
+            dst = None, # Same as src
+            beta = ss.beta(0.2),
+            contacts = ss.poisson(1.0),
+        )
+        self.update_pars(pars, **kwargs)
+
+        self.define_states(
+            ss.FloatArr('eff_contacts', default=self.pars.contacts, label='Effective number of contacts')
+        )
+
+        self.pars.diseases = sc.promotetolist(self.pars.diseases)
+        self.diseases = None
+        self.src_uids = None
+        self.dst_uids = None
+
+        self.p_acquire = ss.bernoulli(p=0) # Placeholder value
+
+        return
+
+    def __len__(self):
+        try:
+            return len(self.pars.dst)
+        except:
+            return 0
+
+    def init_post(self):
+        super().init_post()
+
+        if len(self.pars.diseases) == 0:
+            self.diseases = [d for d in self.sim.diseases.values() if isinstance(d, ss.Infection)] # Assume the user wants all communicable diseases
+        else:
+            self.diseases = []
+            for d in self.pars.diseases:
+                if not isinstance(d, str):
+                    raise TypeError(f'Diseases can be specified as ss.Disease objects or strings, not {type(d)}')
+                if d not in self.sim.diseases:
+                    raise KeyError(f'Could not find disease with name {d} in the list of diseases.')
+
+                dis = self.sim.diseases[d]
+                if not isinstance(dis, ss.Infection):
+                    raise TypeError(f'Cannot create a mixing pool for disease {d}. Mixing pools only work for communicable diseases.')
+                self.diseases.append(dis)
+
+            if len(self.diseases) == 0:
+                raise ValueError('You must specify at least one transmissible disease to use mixing pools')
+        return
+
+    def get_uids(self, func_or_array):
+        if func_or_array is None:
+            return self.sim.people.auids
+        elif callable(func_or_array):
+            return func_or_array(self.sim)
+        elif isinstance(func_or_array, ss.uids):
+            return func_or_array
+        raise Exception('src must be either a callable function, e.g. lambda sim: ss.uids(sim.people.age<5), or an array of uids.')
+
+    def remove_uids(self, uids):
+        """ If UIDs are supplied explicitly, remove them if people die """
+        for key in ['src', 'dst']:
+            inds = self.pars[key]
+            if isinstance(inds, ss.uids):
+                self.pars[key] = inds.remove(uids)
+        return
+
+    def step(self):
+        self.src_uids = self.get_uids(self.pars.src)
+        self.dst_uids = self.get_uids(self.pars.dst)
+        beta = self.pars.beta
+        if isinstance(beta, ss.beta):
+            beta = beta.values # Don't use as a time probability
+
+        if beta == 0:
+            return 0
+
+        if len(self.src_uids) == 0 or len(self.dst_uids) == 0:
+            return 0
+
+        n_new_cases = 0
+        for disease in self.diseases:
+            trans = np.mean(disease.infectious[self.src_uids] * disease.rel_trans[self.src_uids])
+            acq = self.eff_contacts[self.dst_uids] * disease.susceptible[self.dst_uids] * disease.rel_sus[self.dst_uids]
+            p = beta*trans*acq
+
+            self.p_acquire.set(p=p)
+            new_cases = self.p_acquire.filter(self.dst_uids)
+            n_new_cases += len(new_cases)
+            disease.set_prognoses(new_cases)
+
+        return n_new_cases
diff --git a/starsim/parameters.py b/starsim/parameters.py
index db378ca4..cdf9ad8b 100644
--- a/starsim/parameters.py
+++ b/starsim/parameters.py
@@ -68,6 +68,8 @@ def update(self, pars=None, create=False, **kwargs):
                     self._update_dist(key, old, new)
                 elif callable(old): # It's a function: update directly
                     self[key] = new
+                elif isinstance(old, dict):
+                    self[key] = new # Take dictionaries directly, without warning the user
                 else: # Everything else; not used currently but could be
                     warnmsg = f'No known mechanism for handling {type(old)} → {type(new)}; using default'
                     ss.warn(warnmsg)
@@ -101,13 +103,11 @@ def _update_timepar(self, key, old, new):
 
         # 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)
+            self[key] = new
+
+        # It's a dataframe, allow the update -- used for demographics
+        elif isinstance(new, (pd.Series, pd.DataFrame)):
+            self[key] = new # TODO: add validation or convert to timepar
 
         # It's a single number, e.g. dur_inf = 6; set parameters
         elif isinstance(new, Number):
@@ -123,6 +123,12 @@ def _update_timepar(self, key, old, new):
                 self[key] = new # TODO: use an actual set here
             else:
                 old.set(**new)
+
+        # Give up
+        else:
+            errormsg = f'Updating timepar from {type(old)} to {type(new)} is not supported'
+            raise TypeError(errormsg)
+
         return
 
     def _update_dist(self, key, old, new):
@@ -144,6 +150,16 @@ def _update_dist(self, key, old, new):
         elif isinstance(new, list):
             old.set(*new)
 
+        # It's a timepar, set it like a number, e.g. dur_inf = ss.dur(6)
+        elif isinstance(new, ss.TimePar):
+            oldpar = old.pars[0]
+            if isinstance(oldpar, ss.TimePar): # If changing from one timepar to another, validate
+                dur_mismatch = isinstance(oldpar, ss.dur) != isinstance(new, ss.dur)
+                if dur_mismatch:
+                    errormsg = f'Cannot change a duration to a non-duration vice versa: old={type(oldpar)}, new={type(new)}'
+                    raise TypeError(errormsg)
+            old.set(new)
+
         # It's a dict, figure out what to do
         elif isinstance(new, dict):
             newtype = new.get('type')
@@ -156,6 +172,12 @@ def _update_dist(self, key, old, new):
                 else:
                     dist = ss.make_dist(new)
                     self[key] = dist
+
+        # Give up
+        else:
+            errormsg = f'Updating dist from {type(old)} to {type(new)} is not supported'
+            raise TypeError(errormsg)
+
         return
 
     def check_key_mismatch(self, pars):
@@ -186,7 +208,27 @@ class SimPars(Pars):
     directly.
 
     Args:
-        kwargs (dict): any additional kwargs are interpreted as parameter names
+        label (str): The name of the simulation
+        n_agents (int/float): The number of agents to run (default 10,000)
+        total_pop (int/float): If provided, scale the agents to this effective population size
+        pop_scale (float): If provided, use this agent-to-population scale factor (total_pop = n_agents*pop_scale)
+        unit (str): The time unit for the simulation (default 'year'; other choices are 'day', 'week', 'month')
+        start (float/str/date): The starting date for the simulation (default 2000); can be a year or date
+        stop (float/str/date): If provided, the ending date for the simulation (if not provided, calculate from "dur")
+        dur (int): How many timesteps to simulate, if "stop" is not provided (default 50)
+        dt (float): The timestep, in units of "unit" (default 1.0)
+        rand_seed (int): The overall random seed for the simulation (used to set module-specific random seeds)
+        birth_rate (float): If provided, include births with this rate (per 1000 people per year)
+        death_rate (float): If provided, include deaths with this rate (per 1000 people per year)
+        use_aging (bool): Specify whether agents age (by default, agents age if and only if births and/or deaths are included)
+        people (People): If provided, use a pre-existing People object rather than creating one (in which case n_agents will be ignored)
+        networks (str/list/Module): The network module(s); can be a string, single module (i.e. Network), or list
+        demographics (str/list/Module): As above
+        diseases (str/list/Module): As above
+        connectors (str/list/Module): As above
+        interventions (str/list/Module): As above
+        analyzers (str/list/Module): As above
+        verbose (float): How much detail to print (1 = every timestep, 0.1 = every 10 timesteps, etc.)
     """
     def __init__(self, **kwargs):
 
@@ -195,19 +237,17 @@ def __init__(self, **kwargs):
         self.verbose = ss.options.verbose # Whether or not to display information during the run -- options are 0 (silent), 0.1 (some; default), 1 (default), 2 (everything)
 
         # Population parameters
-        self.n_agents  = 10e3  # Number of agents
-        self.total_pop = None  # If defined, used for calculating the scale factor
-        self.pop_scale = None  # How much to scale the population
+        self.n_agents  = 10e3 # Number of agents
+        self.total_pop = None # If defined, used for calculating the scale factor
+        self.pop_scale = None # How much to scale the population
 
         # Simulation parameters
-        self.unit       = 'year' # The time unit to use; options are 'year' (default), 'day', and 'none'
-        self.start      = 2000   # Start of the simulation
-        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
-        self.min_slots  = 100  # Minimum number of slots, useful if the population size is very small
+        self.unit      = ''    # The time unit to use; options are 'year' (default), 'day', 'week', 'month', or 'none'
+        self.start     = None  # Start of the simulation (default 2020)
+        self.stop      = None  # End of the simulation
+        self.dur       = None  # Duration of time to run, if stop isn't specified (default 50 steps of self.unit)
+        self.dt        = 1.0   # Timestep (in units of self.unit)
+        self.rand_seed = 1     # Random seed; if None, don't reset
 
         # Demographic parameters
         self.birth_rate = None
@@ -292,12 +332,19 @@ def validate_total_pop(self):
 
     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)
+
+        # Handle the unit
+        if self.unit == '':
+            self.unit = ss.time.default_unit
+        self.unit = ss.time.validate_unit(self.unit)
+
+        # Handle start
+        if self.start is None:
+            self.start = ss.time.default_start[self.unit]
+
+        # Handle stop and dur
         if self.stop is not None:
-            if self.is_default('dur'):
+            if self.dur is None:
                 self.dur = ss.date_diff(self.start, self.stop, self.unit)
             else:
                 errormsg = f'You can supply either stop ({self.stop}) or dur ({self.dur}) but not both, since one is calculated from the other'
@@ -306,11 +353,9 @@ def validate_time(self):
                 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.dur is not None:
-                self.stop = ss.date_add(self.start, self.dur, self.unit)
-            else:
-                errormsg = 'You must supply either "dur" or "stop".'
-                raise ValueError(errormsg)
+            if self.dur is None:
+                self.dur = ss.time.default_dur
+            self.stop = ss.date_add(self.start, self.dur, self.unit)
         return
 
     def validate_modules(self):
@@ -369,8 +414,12 @@ def validate_demographics(self):
     def validate_networks(self):
         """ Validate networks """
         # Don't allow more than one prenatal or postnatal network
-        prenatal_nets  = {k:nw for k,nw in self.networks.items() if nw.prenatal}
-        postnatal_nets = {k:nw for k,nw in self.networks.items() if nw.postnatal}
+        prenatal_nets = []
+        postnatal_nets = []
+        for k,nw in self.networks.items():
+            if isinstance(nw, ss.Network):
+                if nw.prenatal: prenatal_nets.append(k)
+                if nw.postnatal: postnatal_nets.append(k)
         if len(prenatal_nets) > 1:
             errormsg = f'Starsim currently only supports one prenatal network; prenatal networks are: {prenatal_nets.keys()}'
             raise ValueError(errormsg)
diff --git a/starsim/people.py b/starsim/people.py
index 58ac6ff5..37e021d8 100644
--- a/starsim/people.py
+++ b/starsim/people.py
@@ -135,7 +135,6 @@ def link_sim(self, sim):
             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
 
@@ -280,7 +279,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_year
+            self.age[self.alive.uids] += sim.t.dt_year
         return
 
     def request_death(self, uids):
diff --git a/starsim/results.py b/starsim/results.py
index 721fdca0..48dbd959 100644
--- a/starsim/results.py
+++ b/starsim/results.py
@@ -44,10 +44,41 @@ def __init__(self, name=None, label=None, dtype=float, shape=None, scale=True,
         self.init_values()
         return
 
+    def __repr__(self):
+        cls_name = self.__class__.__name__
+        arrstr = super().__repr__().removeprefix(cls_name)
+        out = f'{cls_name}({self.key}):\narray{arrstr}'
+        return out
+
+    def __str__(self):
+        cls_name = self.__class__.__name__
+        try:
+            minval = self.values.min()
+            meanval = self.values.mean()
+            maxval = self.values.max()
+            valstr = f'min={minval:n}, mean={meanval:n}, max={maxval:n}'
+        except:
+            valstr = f'{self.values}'
+        out = f'{cls_name}({self.key}: {valstr})'
+        return out
+
+    def __getitem__(self, key):
+        """ Allow e.g. result['low'] """
+        if isinstance(key, str):
+            return getattr(self, key)
+        else:
+            return super().__getitem__(key)
+
     @property
     def initialized(self):
         return self.values is not None
 
+    @property
+    def has_dates(self):
+        """ Check whether the time vector uses dates (rather than numbers) """
+        try:    return not sc.isnumber(self.timevec[0])
+        except: return False
+
     def init_values(self, values=None, dtype=None, shape=None, force=False):
         """ Handle values """
         if not self.initialized or force:
@@ -96,19 +127,6 @@ def full_label(self):
             full = f'{modlabel}: {reslabel}'
         return full
 
-    def __repr__(self):
-        cls_name = self.__class__.__name__
-        arrstr = super().__repr__().removeprefix(cls_name)
-        out = f'{cls_name}({self.key}):\narray{arrstr}'
-        return out
-
-    def __getitem__(self, key):
-        """ Allow e.g. result['low'] """
-        if isinstance(key, str):
-            return getattr(self, key)
-        else:
-            return super().__getitem__(key)
-
     def to_df(self, sep='_', rename=False):
         """
         Convert to a dataframe with timevec, value, low, and high columns
@@ -151,7 +169,9 @@ def plot(self, fig=None, ax=None, fig_kw=None, plot_kw=None, fill_kw=None, **kwa
         plt.plot(self.timevec, self.values, **plot_kw)
         plt.title(self.full_label)
         plt.xlabel('Time')
-        sc.commaticks()
+        sc.commaticks(ax)
+        if self.has_dates:
+            sc.dateformatter(ax)
         if (self.values.min() >= 0) and (plt.ylim()[0]<0): # Don't allow axis to go negative if results don't
             plt.ylim(bottom=0)
         return fig
@@ -166,8 +186,24 @@ def __init__(self, module, *args, strict=True, **kwargs):
         super().__init__(type=Result, strict=strict, *args, **kwargs)
         return
 
-    def __repr__(self, *args, **kwargs): # TODO: replace with dataframe summary
-        return super().__repr__(*args, **kwargs)
+    def __repr__(self, indent=2, **kwargs): # kwargs are not used, but are needed for disp() to work
+        string = f'Results({self._module})\n'
+        for i,k,v in self.enumitems():
+            if k != 'timevec':
+                entry = f'{v}'
+                if '\n' in entry: # Check if the string is multi-line
+                    lines = entry.splitlines()
+                    entry = f'{i}. {lines[0]}\n'
+                    entry += '\n'.join(' '*indent + f'{i}.' + line for line in lines[1:])
+                    string += entry + '\n'
+                else:
+                    string += f'{i}. {v}\n'
+        string = string.rstrip()
+        return string
+
+    def disp(self, *args, **kwargs):
+        print(super().__repr__(*args, **kwargs))
+        return
 
     def append(self, arg, key=None):
         """ This is activated by adding as well, e.g. results += result """
@@ -196,6 +232,12 @@ def all_results(self):
         """ Iterator over all results, skipping any nested values """
         return iter(res for res in self.values() if isinstance(res, Result))
 
+    @property
+    def equal_len(self):
+        """ Check if all results are equal length """
+        lengths = [len(res) for res in self.flatten().values()]
+        return len(set(lengths)) == 1
+
     def flatten(self, sep='_', only_results=True):
         """ Turn from a nested dictionary into a flat dictionary, keeping only results by default """
         out = sc.flattendict(self, sep=sep)
@@ -203,12 +245,29 @@ def flatten(self, sep='_', only_results=True):
             out = sc.objdict({k:v for k,v in out.items() if isinstance(v, Result)})
         return out
 
-    def to_df(self, sep='_'):
+    def to_df(self, sep='_', descend=False):
         """ Merge all results dataframes into one """
-        dfs = [res.to_df(sep=sep, rename=True) for res in self.all_results]
-        df = dfs[0]
-        for df2 in dfs[1:]:
-            df = df.merge(df2)
+        if not descend:
+            dfs = [res.to_df(sep=sep, rename=True) for res in self.all_results]
+            if len(dfs):
+                df = dfs[0]
+                for df2 in dfs[1:]:
+                    df = df.merge(df2)
+            else:
+                df = None
+        else:
+            if self.equal_len:
+                flat = self.flatten(sep=sep, only_results=True)
+                flat = dict(timevec=self.timevec) | flat # Prepend the timevec
+                df = sc.dataframe.from_dict(flat)
+            else:
+                df = sc.objdict() # For non-equal lengths, actually return an objdict rather than a dataframe
+                df.sim = self.to_df(sep=sep, descend=False)
+                for k,v in self.items():
+                    if isinstance(v, Results):
+                        thisdf = v.to_df(sep=sep, descend=False) # Only allow one level of nesting
+                        if thisdf is not None:
+                            df[k] = thisdf
         return df
 
     def plot(self, style='fancy', fig_kw=None, plot_kw=None):
diff --git a/starsim/run.py b/starsim/run.py
index 064591c1..1b8b608b 100644
--- a/starsim/run.py
+++ b/starsim/run.py
@@ -251,7 +251,7 @@ def reduce(self, quantiles=None, use_mean=False, bounds=None, output=False):
         rflat = reduced_sim.results.flatten()
         rkeys = list(rflat.keys())
         for rkey in rkeys:
-            raw[rkey] = np.zeros((reduced_sim.npts, len(self.sims)))
+            raw[rkey] = np.zeros((len(reduced_sim), len(self.sims)))
             for s, sim in enumerate(self.sims):
                 flat = sim.results.flatten()
                 raw[rkey][:, s] = flat[rkey]
@@ -429,10 +429,6 @@ def single_run(sim, ind=0, reseed=True, keep_people=False, run_args=None, sim_ar
         sim.pars['rand_seed'] += ind  # Reset the seed, otherwise no point of parallel runs
         ss.set_seed() # Note: may not be needed
 
-    if verbose >= 1:
-        verb = 'Running' if do_run else 'Creating'
-        print(f'{verb} a simulation using seed={sim["rand_seed"]}')
-
     # Handle additional arguments
     for key, val in sim_args.items():
         print(f'Processing {key}:{val}')
diff --git a/starsim/settings.py b/starsim/settings.py
index 10524123..3c389375 100644
--- a/starsim/settings.py
+++ b/starsim/settings.py
@@ -4,7 +4,6 @@
 
     ss.options(verbose=False)
 """
-import os
 import numpy as np
 import sciris as sc
 
@@ -64,22 +63,28 @@ def get_orig_options():
         options = sc.objdict()  # The options
 
         optdesc.verbose = 'Set default level of verbosity (i.e. logging detail): e.g., 0.1 is an update every 10 simulated timesteps.'
-        options.verbose = float(os.getenv('STARSIM_VERBOSE', 0.1))
+        options.verbose = sc.parse_env('STARSIM_VERBOSE', 0.1, 'float')
+
+        optdesc.license = 'Whether to print the license on import'
+        options.license = sc.parse_env('STARSIM_LICENSE', False, 'bool')
 
         optdesc.warnings = 'How warnings are handled: options are "warn" (default), "print", and "error"'
-        options.warnings = str(os.getenv('STARSIM_WARNINGS', 'warn'))
+        options.warnings = sc.parse_env('STARSIM_WARNINGS', 'warn', 'str')
 
         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))
+        options.time_eps = sc.parse_env('STARSIM_TIME_EPS', 1e-6, 'float') # If unit = 'year', corresponds to ~30 seconds
 
         optdesc.sep = 'Set thousands seperator for text output'
-        options.sep = str(os.getenv('STARSIM_SEP', ','))
+        options.sep = sc.parse_env('STARSIM_SEP', ',', 'str')
+
+        optdesc.date_sep = 'Set seperator for dates'
+        options.date_sep = sc.parse_env('STARSIM_DATE_SEP', '.', 'str')
 
         optdesc.precision = 'Set arithmetic precision'
-        options.precision = int(os.getenv('STARSIM_PRECISION', 64))
+        options.precision = sc.parse_env('STARSIM_PRECISION', 64, 'int')
 
         optdesc._centralized = 'If True, revert to centralized random number generation (NOT ADVISED).'
-        options._centralized = False
+        options._centralized = sc.parse_env('STARSIM_CENTRALIZED', False, 'bool')
 
         return optdesc, options
 
diff --git a/starsim/sim.py b/starsim/sim.py
index 260ae6c0..ac6cdaf7 100644
--- a/starsim/sim.py
+++ b/starsim/sim.py
@@ -5,12 +5,13 @@
 import numpy as np
 import sciris as sc
 import starsim as ss
+import matplotlib as mpl
 import matplotlib.pyplot as plt
 
 __all__ = ['Sim', 'AlreadyRunError', 'demo', 'diff_sims', 'check_sims_match']
 
 
-class Sim:
+class Sim(ss.Base):
     """
     The Sim object
 
@@ -29,7 +30,7 @@ class Sim:
         connectors (str/Connector/list): as above, for connectors
         copy_inputs (bool): if True, copy modules as they're inserted into the sim (allowing reuse in other sims, but meaning they won't be updated)
         data (df): a dataframe (or dict) of data, with a column "time" plus data of the form "module.result", e.g. "hiv.new_infections" (used for plotting only)
-        kwargs (dict): merged with pars
+        kwargs (dict): merged with pars; see ss.SimPars for all parameter values
 
     **Examples**::
 
@@ -111,16 +112,6 @@ def __repr__(self):
 
         return string
 
-    @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 """
@@ -135,31 +126,20 @@ def modules(self):
         )
 
     def init(self, **kwargs):
-        """ Perform all initializations for the sim; most heavy lifting is done by the parameters """
+        """ Perform all initializations for the sim """
 
-        # Validation and initialization
+        # Validation and initialization -- this is "pre"
         ss.set_seed(self.pars.rand_seed) # Reset the seed before the population is created -- shouldn't matter if only using Dist objects
         self.pars.validate() # Validate parameters
-        self.init_time_attrs() # Initialize time
+        self.init_time() # Initialize time
         self.init_people(**kwargs) # Initialize the people
+        self.init_sim_attrs()
+        self.init_mods_pre()
 
-        # Move initialized modules to the sim
-        keys = ['label', 'demographics', 'networks', 'diseases', 'interventions', 'analyzers', 'connectors']
-        for key in keys:
-            setattr(self, key, self.pars.pop(key))
-
-        # Initialize all the modules with the sim
-        for mod in self.modules:
-            mod.init_pre(self)
-
-        # Initialize products # TODO: think about moving with other modules
-        for intv in self.interventions():
-            if hasattr(intv, 'product'): # TODO: simplify
-                intv.product.init_pre(self)
-
-        # Final initializations
+        # Final initializations -- this is "post"
         self.init_dists() # Initialize distributions
-        self.init_vals() # Initialize the values in all of the states and networks
+        self.init_people_vals() # Initialize the values in all the states and networks
+        self.init_mod_vals() # Initialize the module values
         self.init_results() # Initialize the results
         self.init_data() # Initialize the data
         self.loop.init() # Initialize the integration loop
@@ -170,15 +150,10 @@ def init(self, **kwargs):
         self.initialized = True
         return self
 
-    def init_time_attrs(self):
+    def init_time(self):
         """ Time indexing; derived values live in the sim rather than in the pars """
-        pars = self.pars
-        self.timevec = ss.make_timevec(pars.start, pars.stop, pars.dt, pars.unit)
+        self.t = ss.Time(pars=self.pars, name='sim', sim=True)
         self.results.timevec = self.timevec # Store the timevec in the results for plotting
-        self.npts = len(self.timevec) # The number of points in the sim
-        self.abs_tvec = 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):
@@ -208,6 +183,19 @@ def init_people(self, verbose=None, **kwargs):
         self.people.link_sim(self)
         return self.people
 
+    def init_sim_attrs(self):
+        """ Move initialized modules to the sim """
+        keys = ['label', 'demographics', 'networks', 'diseases', 'interventions', 'analyzers', 'connectors']
+        for key in keys:
+            setattr(self, key, self.pars.pop(key))
+        return
+
+    def init_mods_pre(self):
+        """ Initialize all the modules with the sim """
+        for mod in self.modules:
+            mod.init_pre(self)
+        return
+
     def init_dists(self):
         """ Initialize the distributions """
         # Initialize all distributions now that everything else is in place
@@ -218,13 +206,13 @@ def init_dists(self):
             self.dists.copy_to_module(mod)
         return
 
-    def init_vals(self):
-        """ Initialize the states and other objects with values """
-
-        # Initialize values in people
+    def init_people_vals(self):
+        """ Initialize the People states with actual values """
         self.people.init_vals()
+        return
 
-        # Initialize values in other modules, including networks and time parameters
+    def init_mod_vals(self):
+        """ Initialize values in other modules, including networks and time parameters """
         for mod in self.modules:
             mod.init_post()
         return
@@ -237,7 +225,7 @@ def reset_time_pars(self, force=True):
 
     def init_results(self):
         """ Create initial results that are present in all simulations """
-        kw = dict(shape=self.npts, timevec=self.timevec, dtype=int, scale=True)
+        kw = dict(shape=self.t.npts, timevec=self.t.timevec, dtype=int, scale=True)
         self.results += [
             ss.Result('n_alive',    label='Number alive', **kw),
             ss.Result('new_deaths', label='Deaths', **kw),
@@ -262,20 +250,19 @@ def start_step(self):
         # Print out progress if needed
         self.elapsed = self.timer.toc(output=True)
         if self.verbose: # Print progress
+            t = self.t
             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:
+            string = f'  Running {simlabel}{t.now("str")} ({t.ti:2.0f}/{t.npts}) ({self.elapsed:0.2f} s) '
+            if self.verbose >= 1:
                 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)
+                if not (t.ti % int(1.0 / self.verbose)):
+                    sc.progressbar(t.ti + 1, t.npts, label=string, length=20, newline=True)
         return
 
     def finish_step(self):
         """ Finish the simulation timestep """
-        self.ti += 1
+        self.t.ti += 1
         return
 
     def run_one_step(self, verbose=None):
@@ -288,7 +275,7 @@ def run_one_step(self, verbose=None):
 
         Note: the verbose here is only for the Loop object, not the sim.
         """
-        self.loop.run(self.now, verbose)
+        self.loop.run(self.t.now(), verbose)
         return self
 
     def run(self, until=None, verbose=None):
@@ -315,9 +302,9 @@ def run(self, until=None, verbose=None):
 
         # 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.t.ti -= 1  # During the run, this keeps track of the next step; restore this be the final day of the sim
             for mod in self.modules: # May not be needed, but keeps it consistent with the sim
-                mod.ti -= 1
+                mod.t.ti -= 1
             self.finalize()
             sc.printv(f'Run finished after {self.elapsed:0.2f} s.\n', 1, self.verbose)
         return self # Allows e.g. ss.Sim().run().plot()
@@ -376,14 +363,14 @@ def get_result(res, func):
 
         # Convert "how" from a string to a dict
         if how == 'default':
-            how = {'n_':'mean', 'new_':'mean', 'cum_':'last', '':'mean'}
+            how = {'n_':'mean', 'new_':'mean', 'cum_':'last', 'timevec':'last', '':'mean'}
         elif isinstance(how, str):
             how = {'':how} # Match everything
 
         summary = sc.objdict()
         flat = sc.flattendict(self.results, sep='_')
         for key, res in flat.items():
-            if '_timevec' not in key: # Skip module-specific time vectors
+            if 'timevec' not in key: # Skip module-specific time vectors
                 try:
                     func = get_func(key, how)
                     entry = get_result(res, func)
@@ -393,11 +380,6 @@ def get_result(res, func):
         self.summary = summary
         return summary
 
-    def disp(self):
-        """ Print a full version of the sim """
-        sc.pr(self)
-        return
-
     def shrink(self, skip_attrs=None, in_place=True):
         """
         "Shrinks" the simulation by removing the people and other memory-intensive
@@ -431,14 +413,18 @@ def shrink(self, skip_attrs=None, in_place=True):
         else:
             return shrunken
 
-    def to_df(self, sep='_'):
-        """ Export results as a Pandas dataframe """
+    def check_results_ready(self, errormsg=None):
+        """ Check that results are ready """
+        if errormsg is None:
+            errormsg = 'Please run the sim first'
         if not self.results_ready:  # pragma: no cover
-            errormsg = 'Please run the sim before exporting the results'
             raise RuntimeError(errormsg)
+        return
 
-        flat = self.results.flatten(sep=sep, only_results=False)
-        df = sc.dataframe.from_dict(flat)
+    def to_df(self, sep='_'):
+        """ Export results as a Pandas dataframe """
+        self.check_results_ready('Please run the sim before exporting the results')
+        df = self.results.to_df(sep=sep, descend=True)
         return df
 
     def save(self, filename=None, keep_people=None, skip_attrs=None, **kwargs):
@@ -548,6 +534,8 @@ def plot(self, key=None, fig=None, style='fancy', show_data=True, fig_kw=None, p
             plot_kw (dict): passed to ``plt.plot()``
             scatter_kw (dict): passed to ``plt.scatter()``, for plotting the data
         """
+        self.check_results_ready('Please run the sim before plotting')
+
         # Configuration
         flat = self.results.flatten()
         n_cols = np.ceil(np.sqrt(len(flat))) # Number of columns of axes
@@ -596,8 +584,13 @@ def plot(self, key=None, fig=None, style='fancy', show_data=True, fig_kw=None, p
                 # Plot results
                 ax.plot(res.timevec, res.values, **plot_kw, label=self.label)
                 ax.set_title(res.full_label)
-                ax.set_xlabel('Time')
+                sc.commaticks(ax)
+                if res.has_dates:
+                    locator = mpl.dates.AutoDateLocator(minticks=2, maxticks=5) # Fewer ticks since lots of plots
+                    sc.dateformatter(ax, locator=locator)
 
+        if self.label is not None:
+            fig.suptitle(self.label)
         sc.figlayout(fig=fig)
 
         return fig
diff --git a/starsim/time.py b/starsim/time.py
new file mode 100644
index 00000000..cb8378fa
--- /dev/null
+++ b/starsim/time.py
@@ -0,0 +1,916 @@
+"""
+Functions and classes for handling time
+"""
+import sciris as sc
+import numpy as np
+import pandas as pd
+import datetime as dt
+import dateutil as du
+import starsim as ss
+
+
+#%% Helper functions
+
+# Classes and objects that are externally visible
+__all__ = ['time_units', 'time_ratio', 'date_add', 'date_diff']
+
+# Allowable time arguments
+time_args = ['start', 'stop', 'dt', 'unit']
+
+# Define available time units
+time_units = sc.objdict(
+    day   = 1.0,
+    week  = 7.0,
+    month = 30.4375, # 365.25/12
+    year  = 365.25, # For consistency with months
+)
+
+# Define defaults
+default_dur = 50
+default_unit = 'year'
+default_start_year = 2000
+default_start_date = '2000-01-01'
+default_start = sc.objdict(
+    {k:default_start_date for k in ['day', 'week', 'month']} |
+    {k:default_start_year for k in ['year', 'unitless']}
+)
+
+
+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 is_unitless(unit1) or is_unitless(unit2):
+            errormsg = f'Cannot convert between units if only one is unitless (unit1={unit1}, unit2={unit2})'
+            raise ValueError(errormsg)
+        if unit1 is None or unit2 is None:
+            errormsg = f'Cannot convert between units when only one has been initialized (unit1={unit1}, unit2={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 = (ss.date(stop) - ss.date(start)).days
+        else:
+            raise NotImplementedError
+    return dur
+
+
+def round_tvec(tvec):
+    """ Round time vectors to a certain level of precision, to avoid floating point errors """
+    decimals = int(-np.log10(ss.options.time_eps))
+    tvec = np.round(np.array(tvec), decimals=decimals)
+    return tvec
+
+
+def years_to_dates(yearvec):
+    """ Convert a numeric year vector to a date vector """
+    datevec = np.array([date(sc.datetoyear(y, reverse=True)) for y in yearvec])
+    return datevec
+
+
+def dates_to_years(datevec):
+    """ Convert a date vector to a numeric year vector"""
+    yearvec = round_tvec([sc.datetoyear(d.date()) for d in datevec])
+    return yearvec
+
+
+def dates_to_days(datevec, start_date):
+    """ Convert a date vector into relative days since start date """
+    start_date = date(start_date)
+    dayvec = np.array([(d - start_date).days for d in datevec])
+    return dayvec
+
+
+#%% Time classes
+
+__all__ += ['date', 'Time']
+
+
+class date(pd.Timestamp):
+    """
+    Define a single date; based on ``pd.Timestamp``
+
+    Args:
+        date (int/float/str/datetime): Any type of date input (ints and floats will be interpreted as years)
+        kwargs (dict): passed to pd.Timestamp()
+
+    **Examples**::
+
+        ss.date(2020) # Returns <2020-01-01>
+        ss.date(year=2020) # Returns <2020-01-01>
+        ss.date(year=2024.75) # Returns <2024-10-01>
+        ss.date('2024-04-04') # Returns <2024-04-04>
+        ss.date(year=2024, month=4, day=4) # Returns <2024-04-04>
+    """
+    def __new__(cls, *args, **kwargs):
+        # Check if a year was supplied, and preprocess it
+        single_year_arg = False
+        if len(args) == 1:
+            arg = args[0]
+            if arg is None:
+                return pd.Timestamp(None)
+            elif sc.isnumber(arg):
+                single_year_arg = True
+            elif isinstance(arg, pd.Timedelta):
+                return pd.Timedelta(arg)
+        year_kwarg = len(args) == 0 and len(kwargs) == 1 and 'year' in kwargs
+        if single_year_arg:
+            return cls.from_year(args[0])
+        if year_kwarg:
+            return cls.from_year(kwargs['year'])
+
+        # Otherwise, proceed as normal
+        out = super(date, cls).__new__(cls, *args, **kwargs)
+        out = cls._reset_class(out)
+        return out
+
+    @classmethod
+    def _reset_class(cls, obj):
+        """ Manually reset the class from pd.Timestamp to ss.date """
+        try:
+            obj.__class__ = date
+        except:
+            warnmsg = f'Unable to convert {out} to ss.date(); proceeding with pd.Timestamp'
+            ss.warn(warnmsg)
+        return obj
+
+    def __repr__(self, bracket=True):
+        """ Show the date in brackets, e.g. <2024.04.04> """
+        _ = ss.options.date_sep
+        y = f'{self.year:04d}'
+        m = f'{self.month:02d}'
+        d = f'{self.day:02d}'
+        string = y + _ + m + _ + d
+        if bracket:
+            string = '<' + string + '>'
+        return string
+
+    def __str__(self):
+        """ Like repr, but just the date, e.g. 2024.04.04 """
+        return self.__repr__(bracket=False)
+
+    def disp(self, **kwargs):
+        """ Show the full object """
+        return sc.pr(self, **kwargs)
+
+    @classmethod
+    def from_year(cls, year):
+        """
+        Convert an int or float year to a date.
+
+        **Examples**::
+
+            ss.date.from_year(2020) # Returns <2020-01-01>
+            ss.date.from_year(2024.75) # Returns <2024-10-01>
+        """
+        if isinstance(year, int):
+            return cls(year=year, month=1, day=1)
+        else:
+            dateobj = sc.datetoyear(year, reverse=True)
+            return cls(dateobj)
+
+    def to_year(self):
+        """
+        Convert a date to a floating-point year
+
+        **Examples**::
+
+            ss.date('2020-01-01').to_year() # Returns 2020.0
+            ss.date('2024-10-01').to_year() # Returns 2024.7486
+        """
+        return sc.datetoyear(self.date())
+
+    def to_pandas(self):
+        """ Convert to a standard pd.Timestamp instance """
+        return pd.Timestamp(self.to_numpy()) # Need to convert to NumPy first or it doesn't do anything
+
+    @staticmethod
+    def _convert_other(other):
+        if isinstance(other, du.relativedelta.relativedelta):
+            tu = ss.time.time_units
+            days = other.days + tu.month*other.months + tu.year*other.years
+            int_days = int(round(days))
+            other = dt.timedelta(days=int_days)
+        elif isinstance(other, ss.dur):
+            factor = ss.time.time_units[other.unit]
+            int_days = int(round(factor*other.v))
+            other = dt.timedelta(days=int_days)
+        return other
+
+    def __add__(self, other):
+        other = self._convert_other(other)
+        out = super().__add__(other)
+        out = date(out)
+        return out
+
+    def __sub__(self, other):
+        other = self._convert_other(other)
+        out = super().__sub__(other)
+        out = date(out)
+        return out
+
+    def __radd__(self, other): return self.__add__(other)
+    def __iadd__(self, other): return self.__add__(other)
+    def __rsub__(self, other): return self.__sub__(other) # TODO: check if this should be reversed
+    def __isub__(self, other): return self.__sub__(other)
+
+
+class Time(sc.prettyobj):
+    """
+    Handle time vectors for both simulations and modules.
+
+    Args:
+        start (float/str/date): the start date for the simulation/module
+        stop (float/str/date): the end date for the simulation/module
+        dt (float): the step size, in units of "unit"
+        unit (str): the time unit; choices are "day", "week", "month", "year", or "unitless"
+        pars (dict): if provided, populate parameter values from this dictionary
+        parent (obj): if provided, populate missing parameter values from a 'parent" ``Time`` instance
+        name (str): if provided, name the ``Time`` object
+        init (bool): whether or not to immediately initialize the Time object
+        sim (bool/Sim): if True, initializes as a sim-specific ``Time`` instance; if a Sim instance, initialize the absolute time vector
+
+    The ``Time`` object, after initialization, has the following attributes:
+
+    - ``ti`` (int): the current timestep
+    -  ``dt_year`` (float): the timestep in units of years
+    - ``npts`` (int): the number of timesteps
+    - ``tvec`` (array): time starting at 0, in self units (e.g. ``[0, 0.1, 0.2, ... 10.0]`` if start=0, stop=10, dt=0.1)
+    - ``absvec`` (array): time relative to sim start, in units of sim units (e.g. ``[366, 373, 380, ...]`` if sim-start=2001, start=2002, sim-unit='day', unit='week')
+    - ``yearvec`` (array): time represented as floating-point years (e.g. ``[2000, 2000.1, 2000.2, ... 2010.0]`` if start=2000, stop=2010, dt=0.1)
+    - ``datevec`` (array): time represented as an array of ``ss.date`` objects (e.g. ``[<2000.01.01>, <2000.02.07>, ... <2010.01.01>]`` if start=2000, stop=2010, dt=0.1)
+    - ``timevec`` (array): the "native" time vector, which always matches one of ``tvec``, ``yearvec``, or ``datevec``
+
+    **Examples**::
+
+        t1 = ss.Time(start=2000, stop=2020, dt=1.0, unit='year') # Years, numeric units
+        t2 = ss.Time(start='2021-01-01', stop='2021-04-04', dt=2.0, unit='day') # Days, date units
+    """
+    def __init__(self, start=None, stop=None, dt=None, unit=None, pars=None, parent=None,
+                 name=None, init=True, sim=None):
+        self.name = name
+        self.start = start
+        self.stop = stop
+        self.dt = dt
+        self.unit = unit
+        self.ti = 0 # The time index, e.g. 0, 1, 2
+
+        # Prepare for later initialization
+        self.dt_year = None
+        self.npts    = None
+        self.tvec    = None
+        self.timevec = None
+        self.datevec = None
+        self.yearvec = None
+        self.abstvec = None
+        self.initialized = False
+
+        # Finalize
+        self.update(pars=pars, parent=parent)
+        if init and self.ready:
+            self.init(sim=sim)
+        return
+
+    def __bool__(self):
+        """ Always truthy """
+        return True
+
+    def __len__(self):
+        """ Length is the number of timepoints """
+        return sc.ifelse(self.npts, 0)
+
+    @property
+    def ready(self):
+        """ Check if all parameters are in place to be initialized """
+        return not any([getattr(self, k) is None for k in time_args])
+
+    @property
+    def is_numeric(self):
+        """ Check whether the fundamental simulation unit is numeric (as opposed to date-based) """
+        try:
+            return sc.isnumber(self.start)
+        except:
+            return False
+
+    @property
+    def is_unitless(self):
+        return is_unitless(self.unit)
+
+    def __setstate__(self, state):
+        """ Custom setstate to unpickle ss.date() instances correctly """
+        self.__dict__.update(state)
+        self._convert_timestamps()
+        return
+
+    def _convert_timestamps(self):
+        """ Replace pd.Timestamp instances with ss.date(); required due to pandas limitations with pickling """
+        objs = [self.start, self.stop]
+        objs += sc.tolist(self.datevec, coerce='full')
+        objs += sc.tolist(self.timevec, coerce='full')
+        objs = [obj for obj in objs if type(obj) == pd.Timestamp]
+        for obj in objs:
+            date._reset_class(obj)
+        return
+
+    def update(self, pars=None, parent=None, reset=True, force=None, **kwargs):
+        """ Reconcile different ways of supplying inputs """
+        pars = sc.mergedicts(pars)
+        stale = False
+
+        for key in time_args:
+            current_val = getattr(self, key, None)
+            parent_val = getattr(parent, key, None)
+            kw_val = kwargs.get(key)
+            par_val = pars.get(key)
+
+            # Special handling for dt: don't inherit dt if the units are different
+            if key == 'dt':
+                if isinstance(parent, Time):
+                    if parent.unit != self.unit:
+                        parent_val = 1.0
+
+            if force is False: # Only update missing (None) values
+                val = sc.ifelse(current_val, kw_val, par_val, parent_val)
+            elif force is None: # Prioritize current value
+                val = sc.ifelse(kw_val, par_val, current_val, parent_val)
+            elif force is True: # Prioritize parent value
+                val = sc.ifelse(kw_val, par_val, parent_val, current_val)
+            else:
+                errormsg = f'Invalid value {force} for force: must be False, None, or True'
+                raise ValueError(errormsg)
+
+            if val != current_val:
+                setattr(self, key, val)
+                stale = True
+
+        if stale and reset and self.initialized:
+            self.init()
+        return
+
+    def init(self, sim=None):
+        """ Initialize all vectors """
+        # Initial validation
+        self.unit = validate_unit(self.unit)
+
+        # Copy missing values from sim
+        if isinstance(sim, ss.Sim):
+            self.unit = sc.ifelse(self.unit, sim.t.unit)
+            if self.unit == sim.t.unit: # Units match, use directly
+                sim_dt = sim.t.dt
+                sim_start = sim.t.start
+                sim_stop = sim.t.stop
+            else: # Units don't match, use datevec instead
+                sim_dt = 1.0 # Don't try to reset the dt if the units don't match
+                sim_start = sim.t.datevec[0]
+                sim_stop = sim.t.datevec[-1]
+            self.dt = sc.ifelse(self.dt, sim_dt)
+            self.start = sc.ifelse(self.start, sim_start)
+            self.stop = sc.ifelse(self.stop, sim_stop)
+
+        # Handle start and stop
+        self.start = self.start if self.is_numeric else date(self.start)
+        self.stop  = self.stop  if self.is_numeric else date(self.stop)
+
+        # Convert start and stop to dates
+        date_start = self.start
+        date_stop = self.stop
+        date_unit = 'year' if not has_units(self.unit) else self.unit # Use year by default
+        dt_year = time_ratio(unit1=date_unit, dt1=self.dt, unit2='year', dt2=1.0) # Timestep in units of years
+        offset = 0
+        if self.is_numeric and date_start == 0:
+            date_start = ss.date(ss.time.default_start[date_unit])
+            date_stop = date_start + ss.dur(date_stop, unit=date_unit)
+            offset = date_start.year
+
+        # If numeric, treat that as the ground truth
+        if self.is_numeric:
+            ratio = time_ratio(unit1=date_unit, unit2='year')
+            timevec = round_tvec(sc.inclusiverange(self.start, self.stop, self.dt))
+            yearvec = round_tvec((timevec-timevec[0])*ratio + offset + timevec[0]) # TODO: simplify
+            datevec = years_to_dates(yearvec)
+
+        # If unitless, just use that
+        elif self.is_unitless:
+            timevec = round_tvec(sc.inclusiverange(self.start, self.stop, self.dt))
+            yearvec = round_tvec(timevec)
+            datevec = timevec
+
+        # If the unit is years, handle that
+        elif date_unit == 'year': # For years, the yearvec is the most robust representation
+            start_year = sc.datetoyear(date_start.date())
+            stop_year = sc.datetoyear(date_stop.date())
+            yearvec = round_tvec(sc.inclusiverange(start_year, stop_year, self.dt))
+            datevec = years_to_dates(yearvec)
+            timevec = datevec
+
+        # Otherwise, use dates as the ground truth
+        else:
+            if int(self.dt) == self.dt: # The step is integer-like, use exactly
+                key = date_unit + 's' # e.g. day -> days
+                datelist = sc.daterange(date_start, date_stop, interval={key:int(self.dt)})
+            else: # Convert to the sim unit instead
+                day_delta = time_ratio(unit1=date_unit, dt1=self.dt, unit2='day', dt2=1.0, as_int=True)
+                if day_delta >= 1:
+                    datelist = sc.daterange(date_start, date_stop, interval={'days':day_delta})
+                else:
+                    errormsg = f'Timestep {dt} is too small; must be at least 1 day'
+                    raise ValueError(errormsg)
+            datevec = np.array([ss.date(d) for d in datelist])
+            yearvec = dates_to_years(datevec)
+            timevec = datevec
+
+        # Store things
+        self.dt_year = dt_year
+        self.npts = len(timevec) # The number of points in the sim
+        self.tvec = round_tvec(np.arange(self.npts)*self.dt) # Absolute time array
+        self.timevec = timevec
+        self.datevec = datevec
+        self.yearvec = yearvec
+        if sim == True: # It's the sim itself, the tvec is the absolute time vector
+            self.abstvec = self.tvec
+        elif sim is not None:
+            self.make_abstvec(sim)
+        else:
+            self.abstvec = None # Intentionally set to None, cannot be used in the sim loop until populated
+        self.initialized = True
+        return
+
+    def make_abstvec(self, sim):
+        """ Convert the current time vector into sim units """
+        # Validation
+        if self.is_unitless != sim.t.is_unitless:
+            errormsg = f'Cannot mix units with unitless time: sim.unit={sim.t.unit} {self.name}.unit={self.unit}'
+            raise ValueError(errormsg)
+
+        # Both are unitless or numeric
+        both_unitless = self.is_unitless and sim.t.is_unitless
+        both_numeric = self.is_numeric and sim.t.is_numeric
+        if both_unitless or both_numeric:
+            abstvec = self.tvec.copy() # Start by copying the current time vector
+            ratio = time_ratio(unit1=self.unit, dt1=1.0, unit2=sim.t.unit, dt2=1.0) # tvec has sim units, but not dt
+            if ratio != 1.0:
+                abstvec *= ratio # TODO: CHECK
+            start_diff = self.start - sim.t.start
+            if start_diff != 0.0:
+                abstvec += start_diff
+
+        # The sim uses years; use yearvec
+        elif sim.t.unit == 'year':
+            abstvec = self.yearvec.copy()
+            abstvec -= sim.t.yearvec[0] # Start relative to sim start
+
+        # Otherwise (days, weeks, months), use datevec and convert to days
+        else:
+            dayvec = dates_to_days(self.datevec, start_date=sim.t.datevec[0])
+            ratio = time_ratio(unit1='day', dt1=1.0, unit2=sim.t.unit, dt2=1.0)
+            abstvec = dayvec*ratio # Convert into sim time units
+
+        self.abstvec = round_tvec(abstvec) # Avoid floating point inconsistencies
+        return
+
+    def now(self, key=None):
+        """
+        Get the current simulation time
+
+        Args:
+            which (str): which type of time to get: default (None), "year", "date", "tvec", or "str"
+
+        **Examples**::
+
+            t = ss.Time(start='2021-01-01', stop='2022-02-02', dt=1, unit='week')
+            t.ti = 25
+            t.now() # Returns <2021-06-25>
+            t.now('date') # Returns <2021-06-25>
+            t.now('year') # Returns 2021.479
+            t.now('str') # Returns '2021-06-25'
+        """
+        if key in [None, 'none', 'time', 'str']:
+            vec = self.timevec
+        elif key == 'tvec':
+            vec = self.tvec
+        elif key == 'date':
+            vec = self.datevec
+        elif key == 'year':
+            vec = self.yearvec
+        else:
+            errormsg = f'Invalid key "{key}": must be None, abs, date, or year'
+            raise ValueError(errormsg)
+
+        now = vec[min(self.ti, len(vec)-1)]
+        if key == 'str':
+            now = f'{now:0.1f}' if isinstance(now, float) else str(now)
+        return now
+
+
+#%% TimePar classes
+
+__all__ += ['TimePar', 'dur', 'days', 'years', 'rate', 'perday', 'peryear',
+            'time_prob', 'beta', 'rate_prob']
+
+class TimePar(ss.BaseArr):
+    """
+    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, v=None, *args, **kwargs):
+        """ Allow TimePars to wrap distributions and return the distributions """
+
+        # Special distribution handling
+        if isinstance(v, ss.Dist):
+            if len(args):
+                errormsg = f'When wrapping a distribution with a TimePar, args not allowed ({args}); use kwargs'
+                raise ValueError(errormsg)
+            dist = v
+            dist.pars[0] = cls(dist.pars[0], **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, v, unit=None, parent_unit=None, parent_dt=None, self_dt=1.0):
+        self.v = v
+        self.unit = unit
+        self.parent_unit = parent_unit
+        self.parent_dt = parent_dt
+        self.self_dt = self_dt
+        self.factor = None
+        self.values = None
+        self.initialized = False
+        self.validate_units()
+        return
+
+    def validate_units(self):
+        """ Check that the units entered are valid """
+        try:
+            self.unit = unit_mapping[self.unit]
+        except KeyError:
+            errormsg = f'Invalid unit "{self.unit}"; must be one of: {sc.strjoin(time_units.keys())}'
+            raise ValueError(errormsg)
+        try:
+            self.parent_unit = unit_mapping[self.parent_unit]
+        except KeyError:
+            errormsg = f'Invalid parent unit "{self.parent_unit}"; must be one of: {sc.strjoin(time_units.keys())}'
+            raise ValueError(errormsg)
+        return
+
+    def init(self, parent=None, parent_unit=None, parent_dt=None, update_values=True, die=True):
+        """ 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)
+
+        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) # If unit isn't defined but parent is, set to parent
+        self.parent_unit = sc.ifelse(self.parent_unit, self.unit) # If parent isn't defined but unit is, set to self
+        self.parent_dt = sc.ifelse(self.parent_dt, self.self_dt, 1.0) # If dt isn't defined, assume 1 (self_dt is defined by default)
+
+        # Calculate the actual conversion factor to be used in the calculations
+        self.update_cached(update_values=update_values, die=die)
+        self.initialized = True
+        self.validate_units()
+        return self
+
+    def __repr__(self):
+        name = self.__class__.__name__
+
+        if self.initialized:
+            if self.factor == 1.0:
+                xstr = ''
+            else:
+                xstr = f', values={self.values}'
+        else:
+            xstr = ', initialized=False'
+
+        if (self.parent_unit is not None) and (self.unit != self.parent_unit):
+            parentstr = f', parent={self.parent_unit}'
+        else:
+            parentstr = ''
+
+        default_dt = sc.ifelse(self.self_dt, 1.0) == 1.0
+        if not default_dt:
+            dtstr = f', self_dt={self.self_dt}'
+        else:
+            dtstr = ''
+
+        # Rather than ss.dur(3, unit='day'), dispaly as ss.days(3)
+        prefixstr = 'ss.'
+        key = (name, self.unit)
+        mapping = {
+            ('dur',  'day'):  'days',
+            ('dur',  'year'): 'years',
+            ('rate', 'day'):  'perday',
+            ('rate', 'year'): 'peryear',
+        }
+
+        if key in mapping and default_dt:
+            prefixstr += mapping[key]
+            unitstr = ''
+        else:
+            prefixstr += name
+            unitstr = f', unit={self.unit}'
+
+        suffixstr = unitstr + parentstr + dtstr + xstr
+
+        return f'{prefixstr}({self.v}{suffixstr})'
+
+    @property
+    def isarray(self):
+        """ Check if the value is an array """
+        return isinstance(self.v, np.ndarray)
+
+    def set(self, v=None, unit=None, parent_unit=None, parent_dt=None, self_dt=None, force=False):
+        """ Set the specified parameter values (ignoring None values) and update stored values """
+        if v           is not None: self.v           = v
+        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_dt     is not None: self.self_dt     = self_dt
+        if self.initialized or force: # Don't try to set these unless it's been initialized
+            self.update_cached()
+        self.validate_units()
+        return self
+
+    def update_cached(self, update_values=True, die=True):
+        """ Update the cached factor and values """
+        try:
+            self.update_factor()
+            if update_values:
+                self.update_values()
+        except TypeError as E: # For a known error, skip silently if die=False
+            if die:
+                errormsg = f'Update failed for {self}. Argument v={self.v} should be a number or array; if a function, use update_values=False. Error: {E}'
+                raise TypeError(errormsg) from E
+        except Exception as E: # For other errors, raise a warning
+            if die:
+                raise E
+            else:
+                tb = sc.traceback(E)
+                warnmsg = f'Uncaught error encountered while updating {self}, but die=False. Traceback:\n{tb}'
+                ss.warn(warnmsg)
+
+        return self
+
+    def update_factor(self):
+        """ Set factor used to multiply the value to get the output """
+        self.factor = time_ratio(unit1=self.unit, dt1=self.self_dt, unit2=self.parent_unit, dt2=self.parent_dt)
+        return
+
+    def update_values(self):
+        """ Convert from self.v to self.values based on self.factor -- must be implemented by derived classes """
+        raise NotImplementedError
+
+    def to(self, unit=None, dt=None):
+        """ Create a new timepar based on the current one but with a different unit and/or dt """
+        new = self.asnew()
+        unit = sc.ifelse(unit, self.parent_unit, self.unit)
+        parent_dt = sc.ifelse(dt, 1.0)
+        new.factor = time_ratio(unit1=self.unit, dt1=self.self_dt, unit2=unit, dt2=parent_dt) # Calculate the new factor
+        new.update_values() # Update values
+        new.v = new.values # Reset the base value
+        new.factor = 1.0 # Reset everything else to be 1
+        new.unit = unit
+        new.self_dt = parent_dt
+        new.parent_unit = unit
+        new.parent_dt = parent_dt
+        return new
+
+    def to_parent(self):
+        """ Create a new timepar with the same units as the parent """
+        unit = self.parent_unit
+        dt = self.parent_dt
+        return self.to(unit=unit, dt=dt)
+
+    def to_json(self):
+        """ Export to JSON """
+        attrs = ['v', 'unit', 'parent_unit', 'parent_dt', 'self_dt', 'factor']
+        out = {'classname': self.__class__.__name__}
+        out.update({attr:getattr(self, attr) for attr in attrs})
+        out['values'] = sc.jsonify(self.values)
+        return out
+
+    # Act like a float -- TODO, add type checking
+    def __add__(self, other): return self.values + other
+    def __sub__(self, other): return self.values - other
+    def __mul__(self, other): return self.asnew().set(v=self.v * other)
+    def __pow__(self, other): return self.values ** other
+    def __truediv__(self, other): return self.asnew().set(v=self.v / other)
+
+    # ...from either side
+    def __radd__(self, other): return other + self.values
+    def __rsub__(self, other): return other - self.values
+    def __rmul__(self, other): return self.asnew().set(v= other * self.v)
+    def __rpow__(self, other): return other ** self.values
+    def __rtruediv__(self, other): return other / self.values # TODO: should be a rate?
+
+    # Handle modify-in-place methods
+    def __iadd__(self, other): return self.set(v=self.v + other)
+    def __isub__(self, other): return self.set(v=self.v - other)
+    def __imul__(self, other): return self.set(v=self.v * other)
+    def __itruediv__(self, other): return self.set(v=self.v / other)
+
+    # Other methods
+    def __neg__(self): return self.asnew().set(v=-self.v)
+
+
+class dur(TimePar):
+    """ Any number that acts like a duration """
+    def update_values(self):
+        self.values = self.v*self.factor
+        return self.values
+
+
+def days(v, parent_unit=None, parent_dt=None):
+    """ Shortcut to ss.dur(value, units='day') """
+    return dur(v=v, unit='day', parent_unit=parent_unit, parent_dt=parent_dt)
+
+
+def years(v, parent_unit=None, parent_dt=None):
+    """ Shortcut to ss.dur(value, units='year') """
+    return dur(v=v, unit='year', parent_unit=parent_unit, parent_dt=parent_dt)
+
+
+class rate(TimePar):
+    """ Any number that acts like a rate; can be greater than 1 """
+    def update_values(self):
+        self.values = self.v/self.factor
+        return self.values
+
+
+def perday(v, parent_unit=None, parent_dt=None):
+    """ Shortcut to ss.rate(value, units='day') """
+    return rate(v=v, unit='day', parent_unit=parent_unit, parent_dt=parent_dt)
+
+
+def peryear(v, parent_unit=None, parent_dt=None):
+    """ Shortcut to ss.rate(value, units='year') """
+    return rate(v=v, unit='year', parent_unit=parent_unit, parent_dt=parent_dt)
+
+
+class time_prob(TimePar):
+    """
+    A probability over time (a.k.a. a cumulative hazard rate); must be >=0 and <=1.
+
+    Note: ``ss.time_prob()`` converts one cumulative hazard rate to another with a
+    different time unit. ``ss.rate_prob()`` converts an exponential rate to a cumulative
+    hazard rate.
+    """
+    def update_values(self):
+        v = self.v
+        if self.isarray:
+            self.values = v.copy()
+            inds = np.logical_and(0.0 < v, v < 1.0)
+            if inds.sum():
+                rates = -np.log(1 - v[inds])
+                self.values[inds] = 1 - np.exp(-rates/self.factor)
+            invalid = np.logical_or(v < 0.0, 1.0 < v)
+            if invalid.sum():
+                errormsg = f'Invalid value {self.v} for {self}: must be 0-1. If using in a calculation, use .values instead.'
+                raise ValueError(errormsg)
+        else:
+            if v == 0:
+                self.values = 0
+            elif v == 1:
+                self.values = 1
+            elif 0 <= v <= 1:
+                rate = -np.log(1 - v)
+                self.values = 1 - np.exp(-rate/self.factor)
+            else:
+                errormsg = f'Invalid value {self.v} for {self}: must be 0-1. If using in a calculation, use .values instead.'
+                raise ValueError(errormsg)
+        return self.values
+
+
+class rate_prob(TimePar):
+    """
+    An instantaneous rate converted to a probability; must be >=0.
+
+    Note: ``ss.time_prob()`` converts one cumulative hazard rate to another with a
+    different time unit. ``ss.rate_prob()`` converts an exponential rate to a cumulative
+    hazard rate.
+    """
+    def update_values(self):
+        v = self.v
+        if self.isarray:
+            self.values = v.copy()
+            inds = v > 0.0
+            if inds.sum():
+                self.values[inds] = 1 - np.exp(-v[inds]/self.factor)
+            invalid = v < 0.0
+            if invalid.sum():
+                errormsg = f'Invalid value {self.v} for {self}: must be >=0. If using in a calculation, use .values instead.'
+                raise ValueError(errormsg)
+        else:
+            if v == 0:
+                self.values = 0
+            elif v > 0:
+                self.values = 1 - np.exp(-v/self.factor)
+            else:
+                errormsg = f'Invalid value {self.value} for {self}: must be >=0. If using in a calculation, use .values instead.'
+                raise ValueError(errormsg)
+        return self.values
+
+
+class beta(time_prob):
+    """ A container for beta (i.e. the disease transmission rate) """
+    pass
+
+
+#%% Unit handling
+
+# Map different units onto the time units -- placed at the end to include the functions
+unit_mapping_reverse = {
+    None: [None],
+    'unitless': ['unitless', 'none'],
+    'day': ['d', 'day', 'days', 'perday', days, perday],
+    'year': ['y', 'yr', 'year', 'years', 'peryear', years, peryear],
+    'week': ['w', 'wk', 'week', 'weeks'],
+    'month': ['m', 'mo', 'month', 'months'],
+}
+unit_mapping = {v:k for k,vlist in unit_mapping_reverse.items() for v in vlist}
+
+
+def validate_unit(unit):
+    """ Check that the unit is valid, and convert to a standard type """
+    try:
+        unit = unit_mapping[unit]
+    except KeyError as E:
+        errormsg = f'Invalid unit "{unit}". Valid units are:\n{sc.pp(unit_mapping_reverse, output=True)}'
+        raise sc.KeyNotFoundError(errormsg) from E
+    return unit
+
+
+def is_unitless(unit):
+    """ Check if explicitly unitless (excludes None; use not has_units() for that) """
+    return unit in unit_mapping_reverse['unitless']
+
+
+def has_units(unit):
+    """ Check that the unit is valid and is not unitless or None """
+    unit = validate_unit(unit)
+    if unit is None:
+        return False
+    elif is_unitless(unit):
+        return False
+    else:
+        return True
diff --git a/starsim/timepars.py b/starsim/timepars.py
deleted file mode 100644
index bb730c45..00000000
--- a/starsim/timepars.py
+++ /dev/null
@@ -1,369 +0,0 @@
-"""
-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_abs_tvec',
-           'TimePar', 'dur', 'days', 'years', 'rate', 'perday', 'peryear', '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_abs_tvec(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(ss.BaseArr):
-    """
-    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, v=None, *args, **kwargs):
-        """ Allow TimePars to wrap distributions and return the distributions """
-
-        # Special distribution handling
-        if isinstance(v, ss.Dist):
-            dist = v
-            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, v, unit=None, parent_unit=None, parent_dt=None, self_dt=1.0):
-        self.v = v
-        self.unit = unit
-        self.parent_unit = parent_unit
-        self.parent_dt = parent_dt
-        self.self_dt = self_dt
-        self.factor = None
-        self.values = None
-        self.parent_name = None
-        self.initialized = False
-        return
-
-    def init(self, parent=None, parent_unit=None, parent_dt=None, update_values=True, die=True):
-        """ 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) # If unit isn't defined but parent is, set to parent
-        self.parent_unit = sc.ifelse(self.parent_unit, self.unit) # If parent isn't defined but unit is, set to self
-        self.parent_dt = sc.ifelse(self.parent_dt, self.self_dt, 1.0) # If dt isn't defined, assume 1 (self_dt is defined by default)
-
-        # Calculate the actual conversion factor to be used in the calculations
-        self.update_cached(update_values=update_values, die=die)
-        self.initialized = True
-        return self
-
-    def __repr__(self):
-        name = self.__class__.__name__
-
-        if self.initialized:
-            if self.factor == 1.0:
-                xstr = ''
-            else:
-                xstr = f', values={self.values}'
-        else:
-            xstr = ', initialized=False'
-
-        if (self.parent_unit is not None) and (self.unit != self.parent_unit):
-            parentstr = f', parent={self.parent_unit}'
-        else:
-            parentstr = ''
-
-        return f'ss.{name}({self.v}, unit={self.unit}{parentstr}{xstr})'
-
-    @property
-    def isarray(self):
-        return isinstance(self.v, np.ndarray)
-
-    def set(self, v=None, unit=None, parent_unit=None, parent_dt=None, self_dt=None, force=False):
-        """ Set the specified parameter values (ignoring None values) and update stored values """
-        if v           is not None: self.v           = v
-        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_dt     is not None: self.self_dt     = self_dt
-        if self.initialized or force: # Don't try to set these unless it's been initialized
-            self.update_cached()
-        return self
-
-    def update_cached(self, update_values=True, die=True):
-        """ Update the cached factor and values """
-        try:
-            self.update_factor()
-            if update_values:
-                self.update_values()
-        except TypeError as E: # For a known error, skip silently if die=False
-            if die:
-                errormsg = f'Update failed for {self}. Argument v={self.v} should be a number or array; if a function, use update_values=False. Error: {E}'
-                raise TypeError(errormsg) from E
-        except Exception as E: # For other errors, raise a warning
-            if die:
-                raise E
-            else:
-                tb = sc.traceback(E)
-                warnmsg = f'Uncaught error encountered while updating {self}, but die=False. Traceback:\n{tb}'
-                ss.warn(warnmsg)
-
-        return self
-
-    def update_factor(self):
-        """ Set factor used to multiply the value to get the output """
-        self.factor = time_ratio(unit1=self.unit, dt1=self.self_dt, unit2=self.parent_unit, dt2=self.parent_dt)
-        return
-
-    def update_values(self):
-        """ Convert from self.v to self.values based on self.factor -- must be implemented by derived classes """
-        raise NotImplementedError
-
-    def to(self, unit=None, dt=None):
-        """ Create a new timepar based on the current one but with a different unit and/or dt """
-        new = self.asnew()
-        unit = sc.ifelse(unit, self.parent_unit, self.unit)
-        parent_dt = sc.ifelse(dt, 1.0)
-        new.factor = time_ratio(unit1=self.unit, dt1=self.self_dt, unit2=unit, dt2=parent_dt) # Calculaet the new factor
-        new.update_values() # Update values
-        new.v = new.values # Reset the base value
-        new.factor = 1.0 # Reset everything else to be 1
-        new.parent_unit = unit
-        new.parent_dt = parent_dt
-        return new
-
-    def to_parent(self):
-        """ Create a new timepar with the same units as the parent """
-        unit = self.parent_unit
-        dt = self.parent_dt
-        return self.to(unit=unit, dt=dt)
-
-    def to_json(self):
-        """ Export to JSON """
-        attrs = ['v', 'unit', 'parent_unit', 'parent_dt', 'parent_name', 'self_dt', 'factor']
-        out = {'classname': self.__class__.__name__}
-        out.update({attr:getattr(self, attr) for attr in attrs})
-        out['values'] = sc.jsonify(self.values)
-        return out
-
-    # Act like a float -- TODO, add type checking
-    def __add__(self, other): return self.values + other
-    def __sub__(self, other): return self.values - other
-    def __mul__(self, other): return self.asnew().set(v=self.v * other)
-    def __pow__(self, other): return self.values ** other
-    def __truediv__(self, other): return self.asnew().set(v=self.v / other)
-
-    # ...from either side
-    def __radd__(self, other): return other + self.values
-    def __rsub__(self, other): return other - self.values
-    def __rmul__(self, other): return self.asnew().set(v= other * self.v)
-    def __rpow__(self, other): return other ** self.values
-    def __rtruediv__(self, other): return other / self.values # TODO: should be a rate?
-
-    # Handle modify-in-place methods
-    def __iadd__(self, other): return self.set(v=self.v + other)
-    def __isub__(self, other): return self.set(v=self.v - other)
-    def __imul__(self, other): return self.set(v=self.v * other)
-    def __itruediv__(self, other): return self.set(v=self.v / other)
-
-    # Other methods
-    def __neg__(self): return self.asnew().set(v=-self.v)
-
-
-class dur(TimePar):
-    """ Any number that acts like a duration """
-    def update_values(self):
-        self.values = self.v*self.factor
-        return self.values
-
-
-class days(dur):
-    """ Shortcut to ss.dur(value, units='day') """
-    def __init__(self, v, parent_unit=None, parent_dt=None):
-        super().__init__(v=v, unit='day', parent_unit=parent_unit, parent_dt=parent_dt)
-        return
-
-
-class years(dur):
-    """ Shortcut to ss.dur(value, units='year') """
-    def __init__(self, v, parent_unit=None, parent_dt=None):
-        super().__init__(v=v, 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 update_values(self):
-        self.values = self.v/self.factor
-        return self.values
-
-
-class perday(rate):
-    """ Shortcut to ss.rate(value, units='day') """
-    def __init__(self, v, parent_unit=None, parent_dt=None):
-        super().__init__(v=v, unit='day', parent_unit=parent_unit, parent_dt=parent_dt)
-        return
-
-
-class peryear(rate):
-    """ Shortcut to ss.rate(value, units='year') """
-    def __init__(self, v, parent_unit=None, parent_dt=None):
-        super().__init__(v=v, unit='year', parent_unit=parent_unit, parent_dt=parent_dt)
-        return
-
-
-class time_prob(TimePar):
-    """ A probability over time (a.k.a. a cumulative hazard rate); must be >0 and <1 """
-    def update_values(self):
-        v = self.v
-        if self.isarray:
-            self.values = v.copy()
-            inds = np.logical_and(0.0 < v, v < 1.0)
-            if inds.sum():
-                rates = -np.log(1 - v[inds])
-                self.values[inds] = 1 - np.exp(-rates/self.factor)
-        else:
-            if v == 0:
-                self.values = 0
-            elif v == 1:
-                self.values = 1
-            elif 0 <= v <= 1:
-                rate = -np.log(1 - v)
-                self.values = 1 - np.exp(-rate/self.factor)
-            else:
-                errormsg = f'Invalid value {self.value} for {self}: must be 0-1'
-                raise ValueError(errormsg)
-        return self.values
-
-
-class beta(time_prob):
-    """ A container for beta (i.e. the disease transmission rate) """
-    pass
diff --git a/starsim/utils.py b/starsim/utils.py
index f5f02da9..c6d6411e 100644
--- a/starsim/utils.py
+++ b/starsim/utils.py
@@ -268,10 +268,8 @@ def standardize_data(data=None, metadata=None, min_year=1800, out_of_range=0, de
           they appeared in the metadata (except for year and age appearing first).
         - An ``ss.Dist`` instance - if the ``data`` input is an ``ss.Dist``, that same object will be returned by this function
     """
-
-    if isinstance(data, ss.Dist):
-        return data
-    elif sc.isnumber(data):
+    # It's a format that can be used directly: return immediately
+    if sc.isnumber(data) or isinstance(data, (ss.Dist, ss.TimePar)):
         return data
 
     # Convert series and dataframe inputs into dicts
diff --git a/starsim/version.py b/starsim/version.py
index 12426412..c7d42cbc 100644
--- a/starsim/version.py
+++ b/starsim/version.py
@@ -4,6 +4,6 @@
 
 __all__ = ['__version__', '__versiondate__', '__license__']
 
-__version__ = '2.0.0'
-__versiondate__ = '2024-10-01'
+__version__ = '2.1.0'
+__versiondate__ = '2024-11-07'
 __license__ = f'Starsim {__version__} ({__versiondate__}) — © 2023-2024 by IDM'


Test patch:
diff --git a/tests/baseline.json b/tests/baseline.json
index 8c6fbc9b..0e7720a9 100644
--- a/tests/baseline.json
+++ b/tests/baseline.json
@@ -1,26 +1,26 @@
 {
   "summary": {
-    "timevec": 2010.0,
-    "births_new": 65.02970297029702,
-    "births_cumulative": 3224.019801980198,
-    "births_cbr": 29.723743148428454,
-    "deaths_new": 44.97029702970297,
-    "deaths_cumulative": 2229.108910891089,
-    "deaths_cmr": 20.555949238807358,
-    "sir_n_susceptible": 2520.4356435643563,
-    "sir_n_infected": 3626.366336633663,
-    "sir_n_recovered": 4792.089108910891,
-    "sir_prevalence": 0.33888892397890347,
-    "sir_new_infections": 130.36633663366337,
-    "sir_cum_infections": 13167.0,
-    "sis_n_susceptible": 4740.722772277228,
-    "sis_n_infected": 6244.287128712871,
-    "sis_prevalence": 0.5587172284995694,
-    "sis_new_infections": 183.87128712871288,
-    "sis_cum_infections": 18571.0,
-    "sis_rel_sus": 0.5182920328461298,
-    "n_alive": 10938.89108910891,
-    "new_deaths": 46.118811881188115,
-    "cum_deaths": 4608.0
+    "timevec": 2020.0,
+    "births_new": 46.84158415841584,
+    "births_cumulative": 2267.7425742574255,
+    "births_cbr": 19.93833106141367,
+    "deaths_new": 9.673267326732674,
+    "deaths_cumulative": 468.5742574257426,
+    "deaths_cmr": 4.118825151847284,
+    "sir_n_susceptible": 2440.029702970297,
+    "sir_n_infected": 3630.3960396039606,
+    "sir_n_recovered": 5676.693069306931,
+    "sir_prevalence": 0.32402568798505815,
+    "sir_new_infections": 122.34653465346534,
+    "sir_cum_infections": 12357.0,
+    "sis_n_susceptible": 4784.306930693069,
+    "sis_n_infected": 6973.504950495049,
+    "sis_prevalence": 0.5720906510271361,
+    "sis_new_infections": 193.5742574257426,
+    "sis_cum_infections": 19551.0,
+    "sis_rel_sus": 0.5019197711850157,
+    "n_alive": 11747.118811881188,
+    "new_deaths": 10.693069306930694,
+    "cum_deaths": 1072.0
   }
 }
\ No newline at end of file
diff --git a/tests/benchmark.json b/tests/benchmark.json
index 992ff85b..00f78bb2 100644
--- a/tests/benchmark.json
+++ b/tests/benchmark.json
@@ -1,12 +1,12 @@
 {
   "time": {
-    "initialize": 0.034,
-    "run": 0.913
+    "initialize": 0.055,
+    "run": 1.013
   },
   "parameters": {
     "n_agents": 10000,
     "dur": 20,
     "dt": 0.2
   },
-  "cpu_performance": 0.9555946088102025
+  "cpu_performance": 0.9665005580733697
 }
\ No newline at end of file
diff --git a/tests/benchmark_full.py b/tests/benchmark_full.py
index b3f40f38..1286d8ad 100644
--- a/tests/benchmark_full.py
+++ b/tests/benchmark_full.py
@@ -19,13 +19,13 @@ def make_run_sim():
     """ Make and run a decent-sized default simulation """
     # Make the people
     ppl = ss.People(n_agents=n_agents)
-    
+
     # Make the components
     sir = ss.SIS()
     hiv = ss.HIV()
     hiv.pars['beta'] = {'mf': [0.15, 0.10], 'maternal': [0.2, 0], 'random': [0,0]}
     networks = [ss.RandomNet(), ss.MFNet(), ss.MaternalNet()]
-    
+
     # Make the sim
     sim = ss.Sim(pars=pars, people=ppl, networks=networks, demographics=ss.Pregnancy(), diseases=[sir, hiv])
 
@@ -36,10 +36,10 @@ def make_run_sim():
 
 
 if __name__ == '__main__':
-    
+
     T = sc.timer()
     for r in range(repeats):
         make_run_sim()
         T.tt(f'Trial {r+1}/{repeats}')
-        
+
     sc.heading(f'Average: {T.mean()*1000:0.0f} ± {T.std()/len(T)**0.5*1000:0.0f} ms')
\ No newline at end of file
diff --git a/tests/benchmark_simple.py b/tests/benchmark_simple.py
index cd3f26a1..b1fdee72 100644
--- a/tests/benchmark_simple.py
+++ b/tests/benchmark_simple.py
@@ -10,7 +10,7 @@
 
     hiv = ss.HIV()
     hiv.pars['beta'] = {'simple_sexual': [0.0008, 0.0004], 'maternal': [0.2, 0]}
-    
+
     sim = ss.Sim(start=1950, stop=2050, people=ppl, networks=networks, demographics=ss.Pregnancy(), diseases=hiv)
     sim.init()
     sim.run()
diff --git a/tests/profile_sim.py b/tests/profile_sim.py
index 283182e4..95051fa8 100644
--- a/tests/profile_sim.py
+++ b/tests/profile_sim.py
@@ -1,13 +1,13 @@
 import sciris as sc
 import starsim as ss
 
-kw = dict(n_agents=100e3, start=2000, dur=100, diseases='sis', networks='random')
+kw = dict(n_agents=10e3, start=2000, dur=100, diseases='sis', networks='random', demographics=True)
 sim = ss.Sim(**kw)
 
 prof = ['profile', 'cprofile', 'time'][1]
 
 if prof == 'profile':
-    prf = sc.profile(run=sim.run, follow=[sim.run, ss.disease.Infection.infect, ss.network.RandomNet.step])
+    prf = sc.profile(run=sim.run, follow=[sim.run])
 
 elif prof == 'cprofile':
     cprf = sc.cprofile(sort='selfpct', mintime=1e-3, maxitems=None, maxfunclen=None, maxpathlen=None)
@@ -16,4 +16,6 @@
 
 else:
     with sc.timer() as T:
-        sim.run()
\ No newline at end of file
+        sim.run()
+
+sim.loop.plot_cpu()
\ No newline at end of file
diff --git a/tests/test_demographics.py b/tests/test_demographics.py
index 54e686e5..4e81a3d4 100644
--- a/tests/test_demographics.py
+++ b/tests/test_demographics.py
@@ -33,7 +33,7 @@ def test_nigeria(which='births', dt=1, start=1995, dur=15, do_plot=False):
         demographics += pregnancy
 
     death_rates = pd.read_csv(datadir/'nigeria_deaths.csv')
-    death = ss.Deaths(pars={'death_rate': death_rates, 'units': 1})
+    death = ss.Deaths(pars={'death_rate': death_rates, 'rate_units': 1})
     demographics += death
 
     # Make people
diff --git a/tests/test_dist.py b/tests/test_dist.py
index 820e5fbc..aebb4b4c 100644
--- a/tests/test_dist.py
+++ b/tests/test_dist.py
@@ -11,6 +11,7 @@
 
 n = 1_000_000 # For testing statistical properties and performance of distributions
 m = 5 # For just testing that drawing works
+dpy = 365.25 # Number of days per year
 sc.options(interactive=False)
 
 
@@ -340,7 +341,7 @@ def test_timepar_dists():
 
     # Check special distributions
     print('Testing Poisson distribution ...')
-    lam1 = ss.rate(365, 'year', parent_unit='year').init()
+    lam1 = ss.rate(dpy, 'year', parent_unit='year').init()
     lam2 = ss.rate(1,   'day',  parent_unit='year').init()
     poi1 = ss.poisson(lam=lam1, strict=False).init()
     poi2 = ss.poisson(lam=lam2, strict=False).init()
@@ -377,7 +378,7 @@ def age_year(module, sim, uids):
     def age_day(module, sim, uids):
         """ Extract age as a day """
         out = sim.people.age[uids].copy()
-        out *= 365 # Convert to days manually
+        out *= dpy # Convert to days manually
         return out
 
     scale = 1e-3 # Set to a very small but nonzero scale
@@ -387,7 +388,7 @@ def age_day(module, sim, uids):
     d4 = ss.lognorm_ex(name='callable', mean=mean, std=scale).init(sim=sim)
     draws3 = d3.rvs(uids)
     draws4 = d4.rvs(uids)
-    age_in_days = sim.people.age[uids]*365
+    age_in_days = sim.people.age[uids]*dpy
     draw_diff = np.abs(draws3 - draws4).mean()
     age_diff = np.abs(age_in_days - draws3).mean()
     print(f'Input ages were:\n{sim.people.age[uids]}')
diff --git a/tests/test_interventions.py b/tests/test_interventions.py
index 0c556ff1..8a1abcaa 100644
--- a/tests/test_interventions.py
+++ b/tests/test_interventions.py
@@ -18,7 +18,7 @@ def run_sir_vaccine(efficacy, leaky=True, do_plot=False):
     sim = ss.Sim(
         n_agents = 1000,
         pars = dict(
-          networks = dict(     
+          networks = dict(
                 type = 'random',
                 n_contacts = 4
           ),
@@ -34,14 +34,14 @@ def run_sir_vaccine(efficacy, leaky=True, do_plot=False):
         dt  = 0.01
     )
     sim.init(verbose=False)
-    
+
     # work out who to vaccinate
     in_trial = sim.people.sir.susceptible.uids
     n_vac = round(len(in_trial) * v_frac)
     in_vac = np.random.choice(in_trial, n_vac, replace=False)
     in_pla = np.setdiff1d(in_trial, in_vac)
     uids = ss.uids(in_vac)
-    
+
     # create and apply the vaccination
     vac = ss.sir_vaccine(efficacy=efficacy, leaky=leaky)
     vac.init_pre(sim)
diff --git a/tests/test_mixingpools.py b/tests/test_mixingpools.py
new file mode 100644
index 00000000..43792cb9
--- /dev/null
+++ b/tests/test_mixingpools.py
@@ -0,0 +1,203 @@
+"""
+Test Sim API
+"""
+
+# %% Imports and settings
+import sys
+import numpy as np
+import sciris as sc
+import matplotlib.pyplot as plt
+import starsim as ss
+import pytest
+
+n_agents = 1_000
+do_plot = False
+sc.options(interactive=False) # Assume not running interactively
+
+
+def test_single_defaults(do_plot=do_plot):
+    """ Test a single MixingPool using defaults """
+    test_name = sys._getframe().f_code.co_name
+    sc.heading(f'Testing {test_name}...')
+    mp = ss.MixingPool()
+    sir = ss.SIR()
+    sim = ss.Sim(diseases=sir, networks=mp, label=test_name)
+    sim.run()
+
+    if do_plot: sim.plot()
+
+    assert(sim.results.sir['cum_infections'][-1] > sim.results.sir['cum_infections'][0]) # There were infections
+    return sim
+
+
+def test_single_uids(do_plot=do_plot):
+    """ Test a single MixingPool by UIDS, pre-initialization configuration """
+    test_name = sys._getframe().f_code.co_name
+    sc.heading(f'Testing {test_name}...')
+
+    n_agents = 10_000
+    k = n_agents // 2
+    mp = ss.MixingPool(
+        src = ss.uids(np.arange(k)),
+        dst = ss.uids(np.arange(k,n_agents))
+    )
+    sir = ss.SIR()
+    sim = ss.Sim(n_agents=n_agents, diseases=sir, networks=mp, label=test_name)
+    sim.run()
+
+    if do_plot: sim.plot()
+
+    assert sim.results.sir['cum_infections'][-1] > sim.results.sir['cum_infections'][0], 'There were no new infections in the simulation'
+    return sim
+
+
+def test_single_ncd():
+    """ Test a single MixingPool with a ncd """
+    test_name = sys._getframe().f_code.co_name
+    sc.heading(f'Testing {test_name}...')
+    mp_pars = {
+        'src': ss.AgeGroup(0, 15),
+        'dst': ss.AgeGroup(15, None),
+        'beta': ss.beta(0.15),
+        'contacts': ss.poisson(lam=5),
+        'diseases': 'ncd'
+    }
+    mp = ss.MixingPool(mp_pars)
+
+    ncd = ss.NCD()
+    sim = ss.Sim(diseases=ncd, networks=mp, label=test_name)
+    with pytest.raises(Exception):
+        sim.run()
+    return sim
+
+
+def test_single_missing_disease():
+    """ Test a single MixingPool with a missing disease """
+    test_name = sys._getframe().f_code.co_name
+    sc.heading(f'Testing {test_name}...')
+    mp_pars = {
+        'src': ss.AgeGroup(0, 15),
+        'dst': ss.AgeGroup(15, None),
+        'beta': ss.beta(0.15),
+        'contacts': ss.poisson(lam=5),
+        'diseases': 'hiv'
+    }
+    mp = ss.MixingPool(mp_pars)
+
+    sir = ss.SIR()
+    sim = ss.Sim(diseases=sir, networks=mp, label=test_name)
+    with pytest.raises(Exception):
+        sim.run()
+
+    return sim
+
+
+def test_single_age(do_plot=do_plot):
+    """ Test a single MixingPool by age """
+    # Incidence must decline because 0-15 --> 15+ transmission only
+    test_name = sys._getframe().f_code.co_name
+    sc.heading(f'Testing {test_name}...')
+    mp_pars = {
+        'src': ss.AgeGroup(0, 15),
+        'dst': ss.AgeGroup(15, None),
+        'beta': ss.beta(0.15),
+        'contacts': ss.poisson(lam=5),
+    }
+    mp = ss.MixingPool(mp_pars)
+
+    sir = ss.SIR()
+    sim = ss.Sim(diseases=sir, networks=mp, label=test_name)
+    sim.run()
+
+    if do_plot: sim.plot()
+
+    assert(sim.results.sir['cum_infections'][-1] > sim.results.sir['cum_infections'][0]) # There were infections
+    return sim
+
+
+def test_single_sex(do_plot=do_plot):
+    """ Test a single MixingPool by sex """
+    # Incidence must decline because M --> F transmission only
+    test_name = sys._getframe().f_code.co_name
+    sc.heading(f'Testing {test_name}...')
+    mp_pars = {
+        'src': lambda sim: sim.people.female, # female to male (only) transmission
+        'dst': lambda sim: sim.people.male,
+        'beta': ss.beta(0.2),
+        'contacts': ss.poisson(lam=4),
+    }
+    mp = ss.MixingPool(mp_pars)
+
+    sir = ss.SIR(init_prev=ss.bernoulli(p=lambda self, sim, uids: 0.05*sim.people.female)) # Seed 5% of the female population
+    sim = ss.Sim(diseases=sir, networks=mp, label=test_name)
+    sim.run()
+
+    if do_plot: sim.plot()
+
+    assert(sim.results.sir['cum_infections'][-1] > sim.results.sir['cum_infections'][0]) # There were infections
+    assert(sim.people.male[sim.diseases.sir.ti_infected>0].all()) # All new infections should be in men
+    return sim
+
+
+def test_multi_defaults(do_plot=do_plot):
+    """ Test MixingPools using defaults """
+    test_name = sys._getframe().f_code.co_name
+    sc.heading(f'Testing {test_name}...')
+    mps = ss.MixingPools(src={'all':None}, dst={'all':None}, contacts=[[1]])
+    sir = ss.SIR()
+    sim = ss.Sim(diseases=sir, networks=mps, label=test_name)
+    sim.run()
+
+    if do_plot: sim.plot()
+
+    assert(sim.results.sir['cum_infections'][-1] > sim.results.sir['cum_infections'][0]) # There were infections
+    return sim
+
+
+def test_multi(do_plot=do_plot):
+    """ Test MixingPools """
+    test_name = sys._getframe().f_code.co_name
+    sc.heading(f'Testing {test_name}...')
+
+    groups = {
+        'Female': lambda sim: sim.people.female,
+        'Male': lambda sim: sim.people.male,
+    }
+
+    mps_pars = dict(
+        contacts = np.array([[1.4, 0.5], [1.2, 0.7]]),
+        beta = ss.beta(0.2),
+        src = groups,
+        dst = groups,
+    )
+    mps = ss.MixingPools(mps_pars)
+
+    sir = ss.SIR()
+    sim = ss.Sim(diseases=sir, networks=mps, label=test_name)
+    sim.run()
+
+    if do_plot: sim.plot()
+
+    assert(sim.results.sir['cum_infections'][-1] > sim.results.sir['cum_infections'][0]) # There were infections
+    return sim
+
+
+if __name__ == '__main__':
+    do_plot = True
+    sc.options(interactive=do_plot)
+    T = sc.timer()
+
+    sim0 = test_single_defaults(do_plot)
+    sim1 = test_single_uids(do_plot)
+    sim2 = test_single_ncd()
+    sim3 = test_single_missing_disease()
+    sim4 = test_single_age(do_plot)
+    sim5 = test_single_sex(do_plot)
+
+    sim6 = test_multi_defaults(do_plot)
+    sim7 = test_multi(do_plot)
+
+    T.toc()
+
+    if do_plot:
+        plt.show()
diff --git a/tests/test_other.py b/tests/test_other.py
index be1a5f77..b0a3e6aa 100644
--- a/tests/test_other.py
+++ b/tests/test_other.py
@@ -64,7 +64,7 @@ def test_microsim(do_plot=False):
     return sim
 
 
-def test_ppl_construction():
+def test_ppl_construction(do_plot=False):
     sc.heading('Test making people and providing them to a sim')
 
     def init_debut(module, sim, uids):
@@ -83,9 +83,10 @@ def init_debut(module, sim, uids):
     sim = ss.Sim(pars=sim_pars, diseases=[gon])
     sim.init()
     sim.run()
-    plt.figure()
-    plt.plot(sim.timevec, sim.results.gonorrhea.n_infected)
-    plt.title('Number of gonorrhea infections')
+    if do_plot:
+        plt.figure()
+        plt.plot(sim.timevec, sim.results.gonorrhea.n_infected)
+        plt.title('Number of gonorrhea infections')
 
     return sim
 
@@ -193,7 +194,7 @@ def test_check_reqiures():
     # Run tests
     ppl = test_people()
     sim1 = test_microsim(do_plot)
-    sim2 = test_ppl_construction()
+    sim2 = test_ppl_construction(do_plot)
     sims = test_arrs()
     sims2 = test_deepcopy()
     sims3 = test_deepcopy_until()
diff --git a/tests/test_plugins.py b/tests/test_plugins.py
index d67b052c..6632d502 100644
--- a/tests/test_plugins.py
+++ b/tests/test_plugins.py
@@ -28,12 +28,12 @@ def __init__(self, pars=None, **kwargs):
 
     def step(self):
         """ Specify HIV-syphilis interactions """
-        
+
         diseases = self.sim.diseases
         syph = diseases.syphilis
         hiv = diseases.hiv
         cd4 = self.sim.people.hiv.cd4
-        
+
         # People with HIV are more likely to acquire syphilis
         syph.rel_sus[cd4 < 500] = self.pars.rel_sus_syph_hiv
         syph.rel_sus[cd4 < 200] = self.pars.rel_sus_syph_aids
@@ -90,23 +90,23 @@ def make_args():
 def test_connectors(do_plot=False):
     """ Test connector example """
     sc.heading('Testing connectors')
-    
+
     # Make arguments
     args = make_args()
     sims = sc.objdict() # List of sims
-    
+
     # Make a sim with a connector, and run
     sims.con = ss.Sim(label='With connector', connectors=hiv_syph(), **args)
     sims.con.run()
-    
+
     # Make a sim without a connector, and run
     sims.nocon = ss.Sim(label='Without connector', **args)
     sims.nocon.run()
-    
+
     # Make a sim with a connector and syph treatment, and run
     sims.treat = ss.Sim(label='With treatment', connectors=hiv_syph(), interventions=Penicillin(), **args)
     sims.treat.run()
-    
+
     # Parse results
     results = sc.odict()
     diseases = ['syphilis', 'hiv']
@@ -118,7 +118,7 @@ def test_connectors(do_plot=False):
     # Plot
     if do_plot:
         plt.figure()
-        
+
         plt.subplot(2,1,1)
         x = sims.con.timevec
         for label,res in results.items():
@@ -128,7 +128,7 @@ def test_connectors(do_plot=False):
         plt.ylabel('Count')
         plt.axvline(2020)
         plt.legend()
-        
+
         plt.subplot(2,1,2)
         for label,res in results.items():
             plt.plot(x, res.hiv, label=label)
@@ -137,15 +137,15 @@ def test_connectors(do_plot=False):
         plt.ylabel('Count')
         plt.axvline(2020)
         plt.legend()
-        
+
         sc.figlayout()
         plt.show()
-    
+
     # Check results
     for disease in diseases:
         assert results[0][disease].sum() > results[1][disease].sum(), f'{disease.title()} infections should be higher with connector'
         assert results[0][disease].sum() > results[2][disease].sum(), f'{disease.title()} infections should be lower with treatment'
-   
+
     return sims
 
 
@@ -153,8 +153,8 @@ def test_connectors(do_plot=False):
     do_plot = True
     sc.options(interactive=do_plot)
     T = sc.timer()
-    
+
     sims = test_connectors(do_plot=do_plot)
-    
+
     T.toc()
-    
+
diff --git a/tests/test_randomness.py b/tests/test_randomness.py
index 5d868022..fba2c33c 100644
--- a/tests/test_randomness.py
+++ b/tests/test_randomness.py
@@ -30,7 +30,7 @@ def make_dists(**kwargs):
 def test_seed():
     """ Test assignment of seeds """
     sc.heading('Testing assignment of seeds')
-    
+
     # Create and initialize two distributions
     dists = make_dists()
     dist0, dist1 = dists.dists.values()
@@ -46,7 +46,7 @@ def test_reset(n=n):
     dists = make_dists()
     distlist = dists.dists.values()
 
-    # Reset via the container, but only 
+    # Reset via the container, but only
     before = sc.autolist()
     after = sc.autolist()
     for dist in distlist:
@@ -58,7 +58,7 @@ def test_reset(n=n):
     print(f'Initial sample:\n{before}')
     print(f'After reset sample:\n{after}')
     assert np.array_equal(before, after)
-    
+
     return before, after
 
 
@@ -80,7 +80,7 @@ def test_jump(n=n):
     print(f'Initial sample:\n{before}')
     print(f'After jump sample:\n{after}')
     assert not np.array_equal(before, after)
-    
+
     return before, after
 
 
@@ -93,9 +93,9 @@ def test_order(n=n):
     # Sample d0, d1
     before = d0(n)
     _ = d1(n)
-    
+
     dists.reset()
-    
+
     # Sample d1, d0
     _ = d1(n)
     after = d0(n)
@@ -114,10 +114,10 @@ class CountInf(ss.Intervention):
     def init_pre(self, sim):
         super().init_pre(sim)
         n_agents = len(sim.people)
-        self.arr = np.zeros((sim.npts, n_agents))
+        self.arr = np.zeros((len(sim), n_agents))
         self.n_agents = n_agents
         return
-    
+
     def step(self):
         self.arr[self.sim.ti, :] = np.array(self.sim.diseases.sir.infected)[:self.n_agents]
         return
@@ -129,26 +129,26 @@ def __init__(self, ti_apply=10):
         super().__init__()
         self.ti_apply = ti_apply
         return
-    
+
     def init_pre(self, sim):
         super().init_pre(sim)
         one_birth = ss.Pregnancy(name='one_birth', rel_fertility=0) # Ensure no default births
         one_birth.init_pre(sim)
         self.one_birth = one_birth
         return
-    
+
     def step(self):
         """ Create an extra agent """
         sim = self.sim
         if sim.ti == self.ti_apply:
             new_uids = self.one_birth.make_embryos(ss.uids(0)) # Assign 0th agent to be the "mother"
             sim.people.age[new_uids] = -100 # Set to a very low number to never reach debut age
-            
+
             # Infect that agent and immediately recover
             sir = sim.diseases.sir
             sir.set_prognoses(new_uids)
             sir.ti_recovered[new_uids] = sim.ti + 1 # Reset recovery time to next timestep
-            
+
             # Reset the random states
             p = sir.pars
             for dist in [p.dur_inf, p.p_death]:
@@ -161,22 +161,22 @@ def plot_infs(s1, s2):
     """ Compare infection arrays from two sims """
     a1 = s1.interventions.countinf.arr
     a2 = s2.interventions.countinf.arr
-    
+
     fig = plt.figure()
     plt.subplot(1,3,1)
     plt.pcolormesh(a1.T)
     plt.xlabel('Timestep')
     plt.ylabel('Person')
     plt.title('Baseline')
-    
+
     plt.subplot(1,3,2)
     plt.pcolormesh(a2.T)
     plt.title('OneMore')
-    
+
     plt.subplot(1,3,3)
     plt.pcolormesh(a2.T - a1.T)
     plt.title('Difference')
-    
+
     sc.figlayout()
     return fig
 
@@ -184,9 +184,9 @@ def plot_infs(s1, s2):
 def test_worlds(do_plot=False):
     """ Test that one extra birth leads to one extra infection """
     sc.heading('Testing worlds...')
-    
+
     res = sc.objdict()
-    
+
     pars = dict(
         start = 2000,
         stop = 2100,
@@ -206,20 +206,20 @@ def test_worlds(do_plot=False):
     )
     s1 = ss.Sim(pars=pars, interventions=CountInf())
     s2 = ss.Sim(pars=pars, interventions=[CountInf(), OneMore()])
-    
+
     s1.run()
     s2.run()
-    
+
     sum1 = s1.summarize()
     sum2 = s2.summarize()
     res.sum1 = sum1
     res.sum2 = sum2
-    
+
     if do_plot:
         s1.plot()
         plot_infs(s1, s2)
         plt.show()
-    
+
     l1 = len(s1.people)
     l2 = len(s2.people)
     i1 = sum1.sir_cum_infections
@@ -229,14 +229,14 @@ def test_worlds(do_plot=False):
     assert l2 == l1 + 1, f'Expected one more person in s2 ({l2}) than s1 ({l1})'
     assert i2 == i1 + 1, f'Expected one more infection in s2 ({i2}) than s1 ({i1})'
     assert (a1 != a2).sum() == 0, f'Expected infection arrays to match:\n{a1}\n{a2}'
-        
+
     return res
 
 
 def test_independence(do_plot=False, thresh=0.1):
     """ Test that when variables are created, they are uncorrelated """
     sc.heading('Testing independence...')
-    
+
     # Create the sim and initialize (do not run)
     sim = ss.Sim(
         n_agents = 1000,
@@ -251,7 +251,7 @@ def test_independence(do_plot=False, thresh=0.1):
         ]
     )
     sim.init()
-    
+
     # Assemble measures
     st = sim.people.states
     arrs = sc.objdict()
@@ -265,7 +265,7 @@ def test_independence(do_plot=False, thresh=0.1):
             for uid in network.edges[p]:
                 data[uid] += 1 # Could also use a histogram
         arrs[key] = data
-    
+
     # Compute the correlations
     n = len(arrs)
     stats = np.zeros((n,n))
@@ -273,7 +273,7 @@ def test_independence(do_plot=False, thresh=0.1):
         for j,arr2 in arrs.enumvals():
             if i != j:
                 stats[i,j] = np.corrcoef(arr1, arr2)[0,1]
-                
+
     # Optionally plot
     if do_plot:
         plt.figure()
@@ -286,13 +286,13 @@ def test_independence(do_plot=False, thresh=0.1):
         plt.colorbar()
         sc.figlayout()
         plt.show()
-            
+
     # Test that everything is independent
     max_corr = abs(stats).max()
     assert max_corr < thresh, f'The maximum correlation between variables was {max_corr}, above the threshold {thresh}'
-    
+
     return sim
-    
+
 
 def test_combine_rands(do_plot=False):
     n = int(1e6)
@@ -310,7 +310,7 @@ def test_combine_rands(do_plot=False):
             plt.title(k)
         sc.figlayout()
         plt.show()
-    
+
     mean = c.mean()
     assert np.isclose(mean, target, atol=atol), f'Expected value to be 0.5±{atol}, not {mean}'
     ks = sps.kstest(c, sps.uniform(0,1).cdf)
diff --git a/tests/test_syphilis.py b/tests/test_syphilis.py
index 3971de33..f914cdaa 100644
--- a/tests/test_syphilis.py
+++ b/tests/test_syphilis.py
@@ -22,7 +22,7 @@ def make_syph_sim(dt=1, n_agents=500):
     # Make demographic modules
     fertility_rates = {'fertility_rate': pd.read_csv(datadir/'nigeria_asfr.csv')}
     pregnancy = ss.Pregnancy(pars=fertility_rates)
-    death_rates = {'death_rate': pd.read_csv(datadir/'nigeria_deaths.csv'), 'units': 1}
+    death_rates = {'death_rate': pd.read_csv(datadir/'nigeria_deaths.csv'), 'rate_units': 1}
     death = ss.Deaths(death_rates)
 
     # Make people and networks
@@ -56,7 +56,7 @@ def __init__(self, *args, **kwargs):
         self.okay = True
         return
 
-    def update_results(self):
+    def step(self):
         """
         Checks states that should be mutually exlusive and collectively exhaustive
         """
@@ -79,7 +79,6 @@ def update_results(self):
         checkall = np.array([s1a, s1b, s2a, s2b, s3])
         if not checkall.all():
             self.okay = False
-
         return
 
 
@@ -92,32 +91,33 @@ def test_syph(dt=1, n_agents=500, do_plot=False):
     # Check plots
     burnin = 0
     pi = int(burnin/dt)
-    
+
     if do_plot:
         tvec = sim.timevec[pi:]
+        res = sim.results.syphilis
         fig, ax = plt.subplots(2, 2)
         ax = ax.ravel()
         ax[0].stackplot(
             tvec,
-            # sim.results.syphilis.n_susceptible[pi:],
-            sim.results.syphilis.n_congenital[pi:],
-            sim.results.syphilis.n_exposed[pi:],
-            sim.results.syphilis.n_primary[pi:],
-            sim.results.syphilis.n_secondary[pi:],
-            (sim.results.syphilis.n_latent_temp[pi:]+sim.results.syphilis.n_latent_long[pi:]),
-            sim.results.syphilis.n_tertiary[pi:],
+            # res.n_susceptible[pi:],
+            res.n_congenital[pi:],
+            res.n_exposed[pi:],
+            res.n_primary[pi:],
+            res.n_secondary[pi:],
+            (res.n_latent_temp[pi:]+res.n_latent_long[pi:]),
+            res.n_tertiary[pi:],
         )
         ax[0].legend(['Congenital', 'Exposed', 'Primary', 'Secondary', 'Latent', 'Tertiary'], loc='lower right')
-    
-        ax[1].plot(tvec, sim.results.syphilis.prevalence[pi:])
+
+        ax[1].plot(tvec, res.prevalence[pi:])
         ax[1].set_title('Syphilis prevalence')
-    
+
         ax[2].plot(tvec, sim.results.n_alive[pi:])
         ax[2].set_title('Population')
-    
-        ax[3].plot(tvec, sim.results.syphilis.new_infections[pi:])
+
+        ax[3].plot(tvec, res.new_infections[pi:])
         ax[3].set_title('New infections')
-    
+
         fig.tight_layout()
         plt.show()
 
@@ -159,7 +159,7 @@ def test_syph_intvs(dt=1, n_agents=500, do_plot=False):
         burnin = 10
         syph_b = sim_base.diseases.syphilis
         syph_i = sim_intv.diseases.syphilis
-        pi = int(burnin/syph_b.dt)
+        pi = int(burnin/syph_b.t.dt)
         plt.figure()
         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')
@@ -183,6 +183,6 @@ def test_syph_intvs(dt=1, n_agents=500, do_plot=False):
 
     sim = test_syph(dt=dt, n_agents=n_agents, do_plot=do_plot)
     sim_base, sim_intv = test_syph_intvs(dt=dt, n_agents=n_agents, do_plot=True)
-    
+
     T.toc()
 
diff --git a/tests/test_time.py b/tests/test_time.py
index 0e758e9d..18fd6d3c 100644
--- a/tests/test_time.py
+++ b/tests/test_time.py
@@ -14,57 +14,57 @@
 # %% 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
-    
+
+    assert ss.time_ratio(unit1='year', unit2='day') == 365.25
+    assert ss.time_ratio(unit1='day', unit2='year') == 1/365.25
+
     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
-    
+    assert d5 + d6 == 2*365.25+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.values, rval/7) # A limitation of this approach, not exact!
     assert np.isclose(r6.values, rval/10)
-    
+
     # Test time_prob
     tpval = 0.1
     tp0 = ss.time_prob(tpval).init(parent_dt=1.0)
@@ -75,9 +75,9 @@ def test_classes():
     assert np.isclose(tp2.values, tpval*2, rtol=0.1)
     assert tp1.values > tpval/2
     assert tp2.values < tpval*2
-    
+
     return d3, d4, r3, r4, tp1
-    
+
 
 def test_units(do_plot=False):
     sc.heading('Test behavior of year vs day units')
@@ -100,53 +100,149 @@ def test_units(do_plot=False):
         networks = rnet,
         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)
-    
+    sims.y = ss.Sim(pars, unit='year', label='Year', start=2000, stop=2002, dt=1/365, verbose=0)
+    sims.d = ss.Sim(pars, unit='day', label='Day', start='2000-01-01', stop='2002-01-01', dt=1, verbose=0)
+
     for sim in sims.values():
         sim.run()
         if do_plot:
             sim.plot()
-    
+
     # Check that results match to within stochastic uncertainty
     rtol = 0.05
     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,
+        verbose = 0,
     )
-    
+
     sim = ss.Sim(pars, unit='day', dt=2, start='2000-01-01', stop='2002-01-01')
     sim.run()
-    
+
+    twoyears = 366*2
+    quarters = 2/0.25
+    assert len(sim) == twoyears//2
+    assert len(sim.diseases.sis) == twoyears
+    assert len(sim.demographics.births) == quarters
+
     if do_plot:
         sim.plot()
-        
+
     return sim
 
 
+def test_mixed_timesteps():
+    sc.heading('Test behavior of different commbinations of timesteps')
+
+    siskw = dict(dur_inf=ss.dur(50, 'day'), beta=ss.beta(0.01, 'day'), waning=ss.rate(0.005, 'day'))
+    kw = dict(n_agents=1000, start='2001-01-01', stop='2001-07-01', networks='random', copy_inputs=False, verbose=0)
+
+    print('Year-year')
+    sis1 = ss.SIS(unit='year', dt=1/365, **sc.dcp(siskw))
+    sim1 = ss.Sim(unit='year', dt=1/365, diseases=sis1, label='year-year', **kw)
+
+    print('Day-day')
+    sis2 = ss.SIS(unit='day', dt=1.0, **sc.dcp(siskw))
+    sim2 = ss.Sim(unit='day', dt=1.0, diseases=sis2, label='day-day', **kw)
+
+    print('Day-year')
+    sis3 = ss.SIS(unit='day', dt=1.0, **sc.dcp(siskw))
+    sim3 = ss.Sim(unit='year', dt=1/365, diseases=sis3, label='day-year', **kw)
+
+    print('Year-day')
+    sis4 = ss.SIS(unit='year', dt=1/365, **sc.dcp(siskw))
+    sim4 = ss.Sim(unit='day', dt=1.0, diseases=sis4, label='year-day', **kw)
+
+    msim = ss.parallel(sim1, sim2, sim3, sim4)
+
+    # Check that al results are close
+    threshold = 0.01
+    summary = msim.summarize()
+    for key,res in summary.items():
+        if res.mean:
+            ratio = res.std/res.mean
+            assert ratio < threshold, f'Result {key} exceeds threshold: {ratio:n} > {threshold}'
+            print(f'✓ Result {key} within threshold: {ratio:n} < {threshold}')
+    return msim
+
+
+def test_time_class():
+    sc.heading('Test different instances of ss.Time')
+
+    def sim(start, stop, dt, unit):
+        """ Generate a fake sim """
+        sim = sc.prettyobj()
+        sim.t = ss.Time(start=start, stop=stop, dt=dt, unit=unit, sim=True)
+        return sim
+
+    print('Testing dates vs. numeric')
+    s1 = sim(start=2000, stop=2002, dt=0.1, unit='year')
+    t1 = ss.Time(start='2001-01-01', stop='2001-06-30', dt=2.0, unit='day')
+    t1.init(sim=s1)
+    assert np.array_equal(s1.t.timevec, s1.t.yearvec)
+    assert len(s1.t.timevec) == 21
+    assert t1.npts == sc.daydiff('2001-01-01', '2001-06-30')//2 + 1
+    assert sc.isnumber(s1.t.start)
+    assert isinstance(t1.start, ss.date)
+    assert s1.t.datevec[-1] == ss.date('2002-01-01')
+    assert t1.datevec[-1] == ss.date('2001-06-30')
+    assert t1.abstvec[0] == 1.0
+
+    print('Testing weeks vs. days')
+    s2 = sim(start='2000-06-01', stop='2001-05-01', dt=1, unit='day')
+    t2 = ss.Time(start='2000-06-01', stop='2001-05-01', dt=1, unit='week')
+    t2.init(sim=s2)
+    assert np.array_equal(s2.t.timevec, s2.t.datevec)
+    assert isinstance(s2.t.start, ss.date)
+    assert t2.npts*7 == s2.t.npts + 1
+    assert np.array_equal(t2.abstvec, s2.t.abstvec[::7])
+
+    print('Testing different units and dt')
+    s3 = sim(start='2001-01-01', stop='2003-01-01', dt=0.1, unit='year')
+    t3 = ss.Time(start='2001-01-01', stop='2003-01-01', dt=2.0, unit='day')
+    t3.init(sim=s3)
+    assert np.array_equal(s3.t.timevec, s3.t.datevec)
+    assert s3.t.datevec[-1] == ss.date('2003-01-01')
+    assert s3.t.npts == 21
+    assert np.isclose(s3.t.abstvec.mean(), t3.abstvec.mean(), atol=1e-3)
+
+    print('Testing unitless')
+    s4 = sim(start=0, stop=10, dt=1.0, unit='unitless')
+    t4 = ss.Time(start=2, stop=9, dt=0.1, unit='unitless')
+    t4.init(sim=s4)
+    assert np.array_equal(s4.t.timevec, s4.t.abstvec)
+    assert s4.t.datevec[0] == ss.date(ss.time.default_start_year)
+    assert len(s4.t) == 11
+    assert len(t4) == (9-2)/0.1+1
+
+    return [s1, t1, s2, t2]
+
+
 # %% 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)
-    
+    o5 = test_mixed_timesteps()
+    o6 = test_time_class()
+
     T.toc()