Repository: starsimhub/starsim (language: python)

Issue:
Ports fix from hpvsim.
Ports fix for this hpvsim issue: https://github.com/amath-idm/hpvsim/issues/575

##### TL:DR: 
Depending on the value of dt, the requested results arrays with an annual frequency over a fixed time period (n_years) may have varying lengths.

##### Full description of problem and suggested changes
 PR here : https://github.com/amath-idm/hpvsim/pull/583


Discussion / hints:

Fix patch:
diff --git a/.gitignore b/.gitignore
new file mode 100644
index 00000000..34a00e22
--- /dev/null
+++ b/.gitignore
@@ -0,0 +1,147 @@
+# HPVsim-generated files
+*.obj
+files/
+modules.rst
+
+# Other files
+.idea
+.vscode
+.DS_Store
+
+
+# Byte-compiled / optimized / DLL files
+__pycache__/
+*.py[cod]
+*$py.class
+
+# C extensions
+*.so
+
+# Distribution / packaging
+.Python
+build/
+develop-eggs/
+dist/
+downloads/
+eggs/
+.eggs/
+lib/
+lib64/
+parts/
+sdist/
+var/
+wheels/
+pip-wheel-metadata/
+share/python-wheels/
+*.egg-info/
+.installed.cfg
+*.egg
+MANIFEST
+
+# PyInstaller
+#  Usually these files are written by a python script from a template
+#  before PyInstaller builds the exe, so as to inject date/other infos into it.
+*.manifest
+*.spec
+
+# Installer logs
+pip-log.txt
+pip-delete-this-directory.txt
+
+# Unit test / coverage reports
+htmlcov/
+.tox/
+.nox/
+.coverage
+.coverage.*
+.cache
+nosetests.xml
+coverage.xml
+*.cover
+*.py,cover
+.hypothesis/
+.pytest_cache/
+
+# Translations
+*.mo
+*.pot
+
+# Django stuff:
+*.log
+local_settings.py
+db.sqlite3
+db.sqlite3-journal
+
+# Flask stuff:
+instance/
+.webassets-cache
+*.out
+
+# Scrapy stuff:
+.scrapy
+
+# Sphinx documentation
+docs/_build/
+docs/modules.rst
+docs/hpvsim.*.rst
+
+# PyBuilder
+target/
+
+# Jupyter Notebook
+.ipynb_checkpoints
+
+# IPython
+profile_default/
+ipython_config.py
+
+# pyenv
+.python-version
+
+# pipenv
+#   According to pypa/pipenv#598, it is recommended to include Pipfile.lock in version control.
+#   However, in case of collaboration, if having platform-specific dependencies or dependencies
+#   having no cross-platform support, pipenv may install dependencies that don't work, or not
+#   install all needed dependencies.
+#Pipfile.lock
+
+# PEP 582; used by e.g. github.com/David-OConnor/pyflow
+__pypackages__/
+
+# Celery stuff
+celerybeat-schedule
+celerybeat.pid
+
+# SageMath parsed files
+*.sage.py
+
+# Environments
+.env
+.venv
+env/
+venv/
+ENV/
+env.bak/
+venv.bak/
+
+# Spyder project settings
+.spyderproject
+.spyproject
+
+# Rope project settings
+.ropeproject
+
+# mkdocs documentation
+/site
+
+# mypy
+.mypy_cache/
+.dmypy.json
+dmypy.json
+
+# Pyre type checker
+.pyre/
+
+# Node.js and NPM
+node_modules
+
diff --git a/CHANGELOG.rst b/CHANGELOG.rst
new file mode 100644
index 00000000..fc00a329
--- /dev/null
+++ b/CHANGELOG.rst
@@ -0,0 +1,20 @@
+==========
+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".
+
+.. contents:: **Contents**
+   :local:
+   :depth: 1
+
+
+Version 0.0.2 (2023-06-29)
+--------------------------
+- Adds in basic *sim functionality
+- *GitHub info*: PR `17 <https://github.com/amath-idm/stisim/pull/17>`__
+
+
+Version 0.0.1 (2023-06-22)
+--------------------------
+- Initial version.
diff --git a/CODE_OF_CONDUCT.rst b/CODE_OF_CONDUCT.rst
new file mode 100644
index 00000000..b24839a9
--- /dev/null
+++ b/CODE_OF_CONDUCT.rst
@@ -0,0 +1,57 @@
+===============
+Code of conduct
+===============
+
+Our pledge
+==========
+
+We believe that a diverse, equitable, and inclusive environment is essential for producing the best quality software. In the interest of fostering an open and welcoming environment, we as contributors and maintainers pledge to making participation in STIsim development and the HPVsim community a harassment-free experience for everyone, regardless of age, body size, disability, ethnicity, sex characteristics, gender identity and expression, level of experience, education, socio-economic status, nationality, personal appearance, race, religion, or sexual identity and orientation.
+
+Our standards
+=============
+
+Examples of behavior that contributes to creating a positive environment include:
+
+*   Using welcoming and inclusive language
+*   Being respectful of differing viewpoints and experiences
+*   Gracefully accepting constructive criticism
+*   Focusing on what is best for the community
+*   Showing empathy towards other community members
+
+Examples of unacceptable behavior by participants include:
+
+*   The use of sexualized language or imagery and unwelcome sexual attention or advances
+*   Trolling, insulting/derogatory comments, and personal or political attacks
+*   Public or private harassment
+*   Publishing others' private information, such as a physical or electronic address, without explicit permission
+*   Other conduct which could reasonably be considered inappropriate in a professional setting
+
+Our responsibilities
+====================
+
+STIsim maintainers are responsible for clarifying the standards of acceptable behavior and will take appropriate and fair corrective action in response to any instances of unacceptable behavior.
+
+STIsim maintainers have the right and responsibility to remove, edit, or reject comments, commits, code, wiki edits, issues, and other contributions that are not aligned to this Code of Conduct, or to ban temporarily or permanently any contributor for other behaviors that they deem inappropriate, threatening, offensive, or harmful.
+
+Scope
+=====
+
+This Code of Conduct applies both within project spaces and in public spaces when an individual is representing STIsim or its community. Examples of representing the STIsim project or community include using an official project e-mail address, posting via an official social media account, or acting as an appointed representative at an online or offline event.
+
+Enforcement
+===========
+
+Instances of abusive, harassing, or otherwise unacceptable behavior may be reported by contacting the project team at robyn.stuart@gatesfoundation.org. All complaints will be reviewed and investigated and will result in a response that is deemed necessary and appropriate to the circumstances. The HPVsim team is obligated to maintain confidentiality with regard to the reporter of an incident. Further details of specific enforcement policies may be posted separately.
+
+HPVsim maintainers who do not follow or enforce the Code of Conduct in good faith may face temporary or permanent repercussions as determined by other members of HPVsim's leadership.
+
+Attribution
+===========
+
+This Code of Conduct is adapted from the `Contributor Covenant`_, version 1.4, available at https://www.contributor-covenant.org/version/1/4/code-of-conduct.html.
+
+.. _Contributor Covenant: https://www.contributor-covenant.org
+
+For answers to common questions about this code of conduct, see the FAQ_.
+
+.. _FAQ: https://www.contributor-covenant.org/faq
diff --git a/setup.py b/setup.py
index f0fe6853..2026c153 100644
--- a/setup.py
+++ b/setup.py
@@ -31,6 +31,7 @@
     description="Sexually Transmitted Infections Simulator (STIsim)",
     long_description=long_description,
     long_description_content_type="text/x-rst",
+    keywords=["agent-based model", "simulation"],
     platforms=["OS Independent"],
     classifiers=CLASSIFIERS,
     packages=find_packages(),
@@ -41,5 +42,6 @@
         'pandas>=2.0.0',
         'sciris>=3.0.0',
         'matplotlib',
+        'seaborn',
     ],
 )
diff --git a/stisim/__init__.py b/stisim/__init__.py
index 0464c832..e63db36e 100644
--- a/stisim/__init__.py
+++ b/stisim/__init__.py
@@ -1,7 +1,30 @@
-from .people import *
-from .population import *
-from .modules import *
-from .interventions import *
-from .analyzers import *
-from .utils import *
-from .sim import *
+from .version import __version__, __versiondate__, __license__
+from .settings      import *
+from .misc          import *
+from .parameters    import *
+from .utils         import *
+from .base          import *
+from .people        import *
+from .networks      import *
+from .sim           import *
+from .modules       import *
+# from .hpv           import *
+
+
+# Import the version and print the license
+if settings.options.verbose:
+    print(__license__)
+
+# # Import data and check
+# from . import data
+# if not data.check_downloaded():
+#     try:
+#         data.quick_download(init=True)
+#     except:
+#         import sciris as sc
+#         errormsg = f"Warning: couldn't download data:\n\n{sc.traceback()}\nProceeding anyway..."
+#         print(errormsg)
+
+# Set the root directory for the codebase
+import pathlib
+rootdir = pathlib.Path(__file__).parent
\ No newline at end of file
diff --git a/stisim/analyzers.py b/stisim/analyzers.py
index f6c19ca1..40671d93 100644
--- a/stisim/analyzers.py
+++ b/stisim/analyzers.py
@@ -1,9 +1,7 @@
-import numpy as np
-import sciris as sc
 from .modules import *
 
-class Analyzer():
 
+class Analyzer:
     requires = []
 
     # TODO - what if the analyzer depends on a specific variable? How does the analyzer know which modules are required?
@@ -20,14 +18,15 @@ def finalize(self, sim):
         pass
 
 
-
 class CD4_analyzer(Analyzer):
-
     requires = [HIV]
 
+    def __init__(self):
+        self.cd4 = None
+
     def initialize(self, sim):
         super().initialize(sim)
         self.cd4 = np.zeros((sim.npts, sim.people.n), dtype=int)
 
     def apply(self, sim):
-        self.cd4[sim.t] = sim.people.hiv.cd4
\ No newline at end of file
+        self.cd4[sim.t] = sim.people.hiv.cd4
diff --git a/stisim/base.py b/stisim/base.py
new file mode 100644
index 00000000..a6f05f7f
--- /dev/null
+++ b/stisim/base.py
@@ -0,0 +1,938 @@
+"""
+Base classes for *sim models
+"""
+
+import numpy as np
+import pandas as pd
+import sciris as sc
+import functools
+from . import utils as ssu
+from . import misc as ssm
+from . import settings as sss
+from .version import __version__
+
+# Specify all externally visible classes this file defines
+__all__ = ['ParsObj', 'BaseSim', 'State', 'StochState', 'BasePeople', 'FlexDict', 'Network']
+
+# Default object getter/setter
+obj_set = object.__setattr__
+base_key = 'uid'  # Define the key used by default for getting length, etc.
+
+
+def rsetattr(obj, attr, val):
+    pre, _, post = attr.rpartition('.')
+    return setattr(rgetattr(obj, pre) if pre else obj, post, val)
+
+
+def rgetattr(obj, attr, *args):
+    def _getattr(obj, attr):
+        return getattr(obj, attr, *args)
+
+    return functools.reduce(_getattr, [obj] + attr.split('.'))
+
+
+# %% Define simulation classes
+
+class FlexPretty(sc.prettyobj):
+    """
+    A class that supports multiple different display options: namely obj.brief()
+    for a one-line description and obj.disp() for a full description.
+    """
+
+    def __repr__(self):
+        """ Use brief repr by default """
+        try:
+            string = self._brief()
+        except Exception as E:
+            string = sc.objectid(self)
+            string += f'Warning, something went wrong printing object:\n{str(E)}'
+        return string
+
+    def _disp(self):
+        """ Verbose output -- use Sciris' pretty repr by default """
+        return sc.prepr(self)
+
+    def disp(self, output=False):
+        """ Print or output verbose representation of the object """
+        string = self._disp()
+        if not output:
+            print(string)
+        else:
+            return string
+
+    def _brief(self):
+        """ Brief output -- use a one-line output, a la Python's default """
+        return sc.objectid(self)
+
+    def brief(self, output=False):
+        """ Print or output a brief representation of the object """
+        string = self._brief()
+        if not output:
+            print(string)
+        else:
+            return string
+
+
+class ParsObj(FlexPretty):
+    """
+    A class based around performing operations on a self.pars dict.
+    """
+
+    def __init__(self, pars, *args, **kwargs):
+        super().__init__(*args, **kwargs)
+        self.update_pars(pars, create=True)
+        return
+
+    def __getitem__(self, key):
+        """ Allow sim['par_name'] instead of sim.pars['par_name'] """
+        try:
+            return self.pars[key]
+        except:
+            all_keys = '\n'.join(list(self.pars.keys()))
+            errormsg = f'Key "{key}" not found; available keys:\n{all_keys}'
+            raise sc.KeyNotFoundError(errormsg)
+
+    def __setitem__(self, key, value):
+        """ Ditto """
+        if key in self.pars:
+            self.pars[key] = value
+        else:
+            all_keys = '\n'.join(list(self.pars.keys()))
+            errormsg = f'Key "{key}" not found; available keys:\n{all_keys}'
+            raise sc.KeyNotFoundError(errormsg)
+        return
+
+    def update_pars(self, pars=None, create=False):
+        """
+        Update internal dict with new pars.
+
+        Args:
+            pars (dict): the parameters to update (if None, do nothing)
+            create (bool): if create is False, then raise a KeyNotFoundError if the key does not already exist
+        """
+        if pars is not None:
+            if not isinstance(pars, dict):
+                raise TypeError(f'The pars object must be a dict; you supplied a {type(pars)}')
+            if not hasattr(self, 'pars'):
+                self.pars = pars
+            if not create:
+                available_keys = list(self.pars.keys())
+                mismatches = [key for key in pars.keys() if key not in available_keys]
+                if len(mismatches):
+                    errormsg = f'Key(s) {mismatches} not found; available keys are {available_keys}'
+                    raise sc.KeyNotFoundError(errormsg)
+            self.pars.update(pars)
+        return
+
+
+def set_metadata(obj, **kwargs):
+    """ Set standard metadata for an object """
+    obj.created = kwargs.get('created', sc.now())
+    obj.version = kwargs.get('version', __version__)
+    obj.git_info = kwargs.get('git_info', ssm.git_info())
+    return
+
+
+class BaseSim(ParsObj):
+    """
+    The BaseSim class stores various methods useful for the Sim that are not directly
+    related to simulating.
+    """
+
+    def __init__(self, *args, **kwargs):
+        super().__init__(*args, **kwargs)  # Initialize and set the parameters as attributes
+        self.filename = None
+        self.initialized = None
+        self.results_ready = None
+        return
+
+    def _disp(self):
+        """
+        Print a verbose display of the sim object. Used by repr(). See sim.disp()
+        for the user version. Equivalent to sc.prettyobj().
+        """
+        return sc.prepr(self)
+
+    def update_pars(self, pars=None, create=False, **kwargs):
+        """ Ensure that metaparameters get used properly before being updated """
+
+        # Merge everything together
+        pars = sc.mergedicts(pars, kwargs)
+        if pars:
+            super().update_pars(pars=pars, create=create)
+
+        return
+
+    @property
+    def n(self):
+        """ Count the number of people -- if it fails, assume none """
+        try:  # By default, the length of the people dict
+            return len(self.people)
+        except:  # pragma: no cover # If it's None or missing
+            return 0
+
+    def shrink(self, skip_attrs=None, in_place=True):
+        """
+        "Shrinks" the simulation by removing the people and other memory-intensive
+        attributes (e.g., some interventions and analyzers), and returns a copy of
+        the "shrunken" simulation. Used to reduce the memory required for RAM or
+        for saved files.
+
+        Args:
+            skip_attrs (list): a list of attributes to skip (remove) in order to perform the shrinking; default "people"
+            in_place (bool): whether to perform the shrinking in place (default), or return a shrunken copy instead
+
+        Returns:
+            shrunken (Sim): a Sim object with the listed attributes removed
+        """
+        # By default, skip people (~90% of memory), popdict, and _orig_pars (which is just a backup)
+        if skip_attrs is None:
+            skip_attrs = ['popdict', 'people', '_orig_pars']
+
+        # Create the new object, and copy original dict, skipping the skipped attributes
+        if in_place:
+            shrunken = self
+            for attr in skip_attrs:
+                setattr(self, attr, None)
+        else:
+            shrunken = object.__new__(self.__class__)
+            shrunken.__dict__ = {k: (v if k not in skip_attrs else None) for k, v in self.__dict__.items()}
+
+        # Don't return if in place
+        if in_place:
+            return
+        else:
+            return shrunken
+
+    def save(self, filename=None, keep_people=None, skip_attrs=None, **kwargs):
+        """
+        Save to disk as a gzipped pickle.
+
+        Args:
+            filename (str or None): the name or path of the file to save to; if None, uses stored
+            keep_people (bool or None): whether to keep the people
+            skip_attrs (list): attributes to skip saving
+            kwargs: passed to sc.makefilepath()
+
+        Returns:
+            filename (str): the validated absolute path to the saved file
+
+        **Example**::
+
+            sim.save() # Saves to a .sim file
+        """
+
+        # Set keep_people based on whether we're in the middle of a run
+        if keep_people is None:
+            if self.initialized and not self.results_ready:
+                keep_people = True
+            else:
+                keep_people = False
+
+        # Handle the filename
+        if filename is None:
+            filename = self.simfile
+        filename = sc.makefilepath(filename=filename, **kwargs)
+        self.filename = filename  # Store the actual saved filename
+
+        # Handle the shrinkage and save
+        if skip_attrs or not keep_people:
+            obj = self.shrink(skip_attrs=skip_attrs, in_place=False)
+        else:
+            obj = self
+        ssm.save(filename=filename, obj=obj)
+
+        return filename
+
+    @staticmethod
+    def load(filename, *args, **kwargs):
+        """
+        Load from disk from a gzipped pickle.
+        """
+        sim = ssm.load(filename, *args, **kwargs)
+        if not isinstance(sim, BaseSim):  # pragma: no cover
+            errormsg = f'Cannot load object of {type(sim)} as a Sim object'
+            raise TypeError(errormsg)
+        return sim
+
+    def _get_ia(self, which, label=None, partial=False, as_list=False, as_inds=False, die=True, first=False):
+        """ Helper method for get_interventions() and get_analyzers(); see get_interventions() docstring """
+
+        # Handle inputs
+        if which not in ['interventions', 'analyzers']:  # pragma: no cover
+            errormsg = f'This method is only defined for interventions and analyzers, not "{which}"'
+            raise ValueError(errormsg)
+
+        ia_list = sc.tolist(
+            self.analyzers if which == 'analyzers' else self.interventions)  # List of interventions or analyzers
+        n_ia = len(ia_list)  # Number of interventions/analyzers
+
+        if label == 'summary':  # Print a summary of the interventions
+            df = pd.DataFrame(columns=['ind', 'label', 'type'])
+            for ind, ia_obj in enumerate(ia_list):
+                df = df.append(dict(ind=ind, label=str(ia_obj.label), type=type(ia_obj)), ignore_index=True)
+            print(f'Summary of {which}:')
+            print(df)
+            return
+
+        else:  # Standard usage case
+            position = 0 if first else -1  # Choose either the first or last element
+            if label is None:  # Get all interventions if no label is supplied, e.g. sim.get_interventions()
+                label = np.arange(n_ia)
+            if isinstance(label, np.ndarray):  # Allow arrays to be provided
+                label = label.tolist()
+            labels = sc.promotetolist(label)
+
+            # Calculate the matches
+            matches = []
+            match_inds = []
+            for label in labels:
+                if sc.isnumber(label):
+                    matches.append(ia_list[label])  # This will raise an exception if an invalid index is given
+                    label = n_ia + label if label < 0 else label  # Convert to a positive number
+                    match_inds.append(label)
+                elif sc.isstring(label) or isinstance(label, type):
+                    for ind, ia_obj in enumerate(ia_list):
+                        if sc.isstring(label) and ia_obj.label == label or (partial and (label in str(ia_obj.label))):
+                            matches.append(ia_obj)
+                            match_inds.append(ind)
+                        elif isinstance(label, type) and isinstance(ia_obj, label):
+                            matches.append(ia_obj)
+                            match_inds.append(ind)
+                else:  # pragma: no cover
+                    errormsg = f'Could not interpret label type "{type(label)}": should be str, int, list, or {which} class'
+                    raise TypeError(errormsg)
+
+            # Parse the output options
+            if as_inds:
+                output = match_inds
+            elif as_list:  # Used by get_interventions()
+                output = matches
+            else:
+                if len(matches) == 0:  # pragma: no cover
+                    if die:
+                        errormsg = f'No {which} matching "{label}" were found'
+                        raise ValueError(errormsg)
+                    else:
+                        output = None
+                else:
+                    output = matches[
+                        position]  # Return either the first or last match (usually), used by get_intervention()
+
+            return output
+
+    def get_interventions(self, label=None, partial=False, as_inds=False):
+        """
+        Find the matching intervention(s) by label, index, or type. If None, return
+        all interventions. If the label provided is "summary", then print a summary
+        of the interventions (index, label, type).
+
+        Args:
+            label (str, int, Intervention, list): the label, index, or type of intervention to get; if a list, iterate over one of those types
+            partial (bool): if true, return partial matches (e.g. 'beta' will match all beta interventions)
+            as_inds (bool): if true, return matching indices instead of the actual interventions
+        """
+        return self._get_ia('interventions', label=label, partial=partial, as_inds=as_inds, as_list=True)
+
+    def get_intervention(self, label=None, partial=False, first=False, die=True):
+        """
+        Like get_interventions(), find the matching intervention(s) by label,
+        index, or type. If more than one intervention matches, return the last
+        by default. If no label is provided, return the last intervention in the list.
+
+        Args:
+            label (str, int, Intervention, list): the label, index, or type of intervention to get; if a list, iterate over one of those types
+            partial (bool): if true, return partial matches (e.g. 'beta' will match all beta interventions)
+            first (bool): if true, return first matching intervention (otherwise, return last)
+            die (bool): whether to raise an exception if no intervention is found
+        """
+        return self._get_ia('interventions', label=label, partial=partial, first=first, die=die, as_inds=False,
+                            as_list=False)
+
+    def get_analyzers(self, label=None, partial=False, as_inds=False):
+        """
+        Same as get_interventions(), but for analyzers.
+        """
+        return self._get_ia('analyzers', label=label, partial=partial, as_list=True, as_inds=as_inds)
+
+    def get_analyzer(self, label=None, partial=False, first=False, die=True):
+        """
+        Same as get_intervention(), but for analyzers.
+        """
+        return self._get_ia('analyzers', label=label, partial=partial, first=first, die=die, as_inds=False,
+                            as_list=False)
+
+
+# %% Define people classes
+
+class State(sc.prettyobj):
+    def __init__(self, name, dtype, fill_value=0, shape=None, label=None):
+        """
+        Args:
+            name: name of the result as used in the model
+            dtype: datatype
+            fill_value: default value for this state upon model initialization
+            shape: If not none, set to match a string in `pars` containing the dimensionality
+            label: text used to construct labels for the result for displaying on plots and other outputs
+        """
+        self.name = name
+        self.dtype = dtype
+        self.fill_value = fill_value
+        self.shape = shape
+        self.label = label or name
+        return
+
+    @property
+    def ndim(self):
+        return len(sc.tolist(self.shape)) + 1
+
+    def new(self, n):
+        shape = sc.tolist(self.shape)
+        shape.append(n)
+        return np.full(shape, dtype=self.dtype, fill_value=self.fill_value)
+
+
+class StochState(State):
+    def __init__(self, name, dtype, distdict=None, **kwargs):
+        super().__init__(name, dtype, kwargs)
+        self.distdict = distdict
+        return
+
+    def new(self, n):
+        shape = sc.tolist(self.shape)
+        shape.append(n)
+        return ssu.sample(**self.distdict, size=tuple(shape))
+
+
+base_states = ssu.named_dict(
+    State('uid', sss.default_int),
+    State('age', sss.default_float),
+    State('female', bool, False),
+    State('debut', sss.default_float),
+    State('dead', bool, False),
+    State('ti_dead', sss.default_float, np.nan),  # Time index for death
+    State('scale', sss.default_float, 1.0),
+)
+
+
+class BasePeople(FlexPretty):
+    """
+    A class to handle all the boilerplate for people -- note that as with the
+    BaseSim vs Sim classes, everything interesting happens in the People class,
+    whereas this class exists to handle the less interesting implementation details.
+    """
+
+    def __init__(self, n, states=None, *args, **kwargs):
+        """ Initialize essential attributes """
+
+        super().__init__(*args, **kwargs)
+        self.initialized = False
+        self.version = __version__  # Store version info
+
+        # Initialize states, networks, modules
+        self.states = sc.mergedicts(base_states, states)
+        self.networks = ssu.named_dict()
+        self._modules = sc.autolist()
+
+        # Private variables relating to dynamic allocation
+        self._data = dict()
+        self._n = n  # Number of agents (initial)
+        self._s = self._n  # Underlying array sizes
+        self._inds = None  # No filtering indices
+
+        # Initialize underlying storage and map arrays
+        for state_name, state in self.states.items():
+            self._data[state_name] = state.new(self._n)
+        self._map_arrays()
+        self['uid'][:] = np.arange(self._n)
+
+        # Define lock attribute here, since BasePeople.lock()/unlock() requires it
+        self._lock = False  # Prevent further modification of keys
+
+        return
+
+    def initialize(self, popdict=None):
+        """ Initialize people by setting their attributes """
+        if popdict is None:
+            self['age'][:] = np.random.random(size=self.n) * 100
+            self['female'][:] = np.random.choice([False, True], size=self.n)
+        else:
+            # Use random defaults
+            self['age'][:] = popdict['age']
+            self['female'][:] = popdict['female']
+        self.initialized = True
+        return
+
+    def __len__(self):
+        """ Length of people """
+        try:
+            arr = getattr(self, base_key)
+            return len(arr)
+        except Exception as E:
+            print(f'Warning: could not get length of People (could not get self.{base_key}: {E})')
+            return 0
+
+    @property
+    def n(self):
+        return len(self)
+
+    def _len_arrays(self):
+        """ Length of underlying arrays """
+        return len(self._data[base_key])
+
+    def lock(self):
+        """ Lock the people object to prevent keys from being added """
+        self._lock = True
+        return
+
+    def unlock(self):
+        """ Unlock the people object to allow keys to be added """
+        self._lock = False
+        return
+
+    def _grow(self, n):
+        """
+        Increase the number of agents stored
+
+        Automatically reallocate underlying arrays if required
+        
+        Args:
+            n (int): Number of new agents to add
+        """
+        orig_n = self._n
+        new_total = orig_n + n
+        if new_total > self._s:
+            n_new = max(n, int(self._s / 2))  # Minimum 50% growth
+            for state_name, state in self.states.items():
+                self._data[state_name] = np.concatenate([self._data[state_name], state.new(n_new)],
+                                                        axis=self._data[state_name].ndim - 1)
+            self._s += n_new
+        self._n += n
+        self._map_arrays()
+        new_inds = np.arange(orig_n, self._n)
+        return new_inds
+
+    def _map_arrays(self, keys=None):
+        """
+        Set main simulation attributes to be views of the underlying data
+
+        This method should be called whenever the number of agents required changes
+        (regardless of whether the underlying arrays have been resized)
+        """
+        row_inds = slice(None, self._n)
+
+        # Handle keys
+        if keys is None: keys = self.states.keys()
+        keys = sc.tolist(keys)
+
+        # Map arrays for selected keys
+        for k in keys:
+            arr = self._data[k]
+            if arr.ndim == 1:
+                rsetattr(self, k, arr[row_inds])
+            elif arr.ndim == 2:
+                rsetattr(self, k, arr[:, row_inds])
+            else:
+                errormsg = 'Can only operate on 1D or 2D arrays'
+                raise TypeError(errormsg)
+
+        return
+
+    def __getitem__(self, key):
+        """ Allow people['attr'] instead of getattr(people, 'attr')
+            If the key is an integer, alias `people.person()` to return a `Person` instance
+        """
+        if isinstance(key, int):
+            return self.person(key)
+        else:
+            return self.__getattribute__(key)
+
+    def __setitem__(self, key, value):
+        """ Ditto """
+        if self._lock and key not in self.__dict__:  # pragma: no cover
+            errormsg = f'Key "{key}" is not an attribute of people and the people object is locked; see people.unlock()'
+            raise AttributeError(errormsg)
+        return self.__setattr__(key, value)
+
+    def __iter__(self):
+        """ Iterate over people """
+        for i in range(len(self)):
+            yield self[i]
+
+    def _brief(self):
+        """
+        Return a one-line description of the people -- used internally and by repr();
+        see people.brief() for the user version.
+        """
+        try:
+            string = f'People(n={len(self):0n})'
+        except Exception as E:  # pragma: no cover
+            string = sc.objectid(self)
+            string += f'Warning, sim appears to be malformed:\n{str(E)}'
+        return string
+
+    def set(self, key, value):
+        """
+        Set values. Note that this will raise an exception the shapes don't match,
+        and will automatically cast the value to the existing type
+        """
+        self[key][:] = value[:]
+
+    def get(self, key):
+        """ Convenience method -- key can be string or list of strings """
+        if isinstance(key, str):
+            return self[key]
+        elif isinstance(key, list):
+            arr = np.zeros((len(self), len(key)))
+            for k, ky in enumerate(key):
+                arr[:, k] = self[ky]
+            return arr
+
+    @property
+    def alive(self):
+        """ Alive boolean """
+        return ~self.dead
+
+    @property
+    def f_inds(self):
+        """ Indices of everyone female """
+        return self.true('female')
+
+    @property
+    def m_inds(self):
+        """ Indices of everyone male """
+        return self.false('female')
+
+    @property
+    def active(self):
+        """ Indices of everyone sexually active  """
+        return (self.age >= self.debut) & self.alive
+
+    @property
+    def int_age(self):
+        """ Return ages as an integer """
+        return np.array(self.age, dtype=sss.default_int)
+
+    @property
+    def round_age(self):
+        """ Rounds age up to the next highest integer"""
+        return np.array(np.ceil(self.age))
+
+    @property
+    def alive_inds(self):
+        """ Indices of everyone alive """
+        return self.true('alive')
+
+    @property
+    def n_alive(self):
+        """ Number of people alive """
+        return len(self.alive_inds)
+
+    def true(self, key):
+        """ Return indices matching the condition """
+        return self[key].nonzero()[-1]
+
+    def false(self, key):
+        """ Return indices not matching the condition """
+        return (~self[key]).nonzero()[-1]
+
+    def defined(self, key):
+        """ Return indices of people who are not-nan """
+        return (~np.isnan(self[key])).nonzero()[0]
+
+    def undefined(self, key):
+        """ Return indices of people who are nan """
+        return np.isnan(self[key]).nonzero()[0]
+
+    def count(self, key, weighted=True):
+        """ Count the number of people for a given key """
+        inds = self[key].nonzero()[0]
+        if weighted:
+            out = self.scale[inds].sum()
+        else:
+            out = len(inds)
+        return out
+
+    def count_any(self, key, weighted=True):
+        """ Count the number of people for a given key for a 2D array if any value matches """
+        inds = self[key].sum(axis=0).nonzero()[0]
+        if weighted:
+            out = self.scale[inds].sum()
+        else:
+            out = len(inds)
+        return out
+
+    def keys(self):
+        """ Returns keys for all non-derived properties of the people object """
+        return [state.name for state in self.states]
+
+    def indices(self):
+        """ The indices of each people array """
+        return np.arange(len(self))
+
+    def to_arr(self):
+        """ Return as numpy array """
+        arr = np.empty((len(self), len(self.keys())), dtype=sss.default_float)
+        for k, key in enumerate(self.keys()):
+            if key == 'uid':
+                arr[:, k] = np.arange(len(self))
+            else:
+                arr[:, k] = self[key]
+        return arr
+
+    def to_list(self):
+        """ Return all people as a list """
+        return list(self)
+
+
+class FlexDict(dict):
+    """
+    A dict that allows more flexible element access: in addition to obj['a'],
+    also allow obj[0]. Lightweight implementation of the Sciris odict class.
+    """
+
+    def __getitem__(self, key):
+        """ Lightweight odict -- allow indexing by number, with low performance """
+        try:
+            return super().__getitem__(key)
+        except KeyError as KE:
+            try:  # Assume it's an integer
+                dictkey = self.keys()[key]
+                return self[dictkey]
+            except:
+                raise sc.KeyNotFoundError(KE)  # Raise the original error
+
+    def keys(self):
+        return list(super().keys())
+
+    def values(self):
+        return list(super().values())
+
+    def items(self):
+        return list(super().items())
+
+
+class Network(FlexDict):
+    """
+    A small class holding a single network of contact edges (connections) between people.
+
+    The input is typically arrays including: person 1 of the connection, person 2 of
+    the connection, the weight of the connection, the duration and start/end times of
+    the connection. Connections are undirected; each person is both a source and sink.
+
+    Args:
+        p1 (array): an array of N connections, representing people on one side of the connection
+        p2 (array): an array of people on the other side of the connection
+        beta (array): an array representing relative transmissibility for this network - TODO, do we need this?
+        label (str): the name of the network (optional)
+        kwargs (dict): other keys copied directly into the network
+
+    Note that all arguments (except for label) must be arrays of the same length,
+    although not all have to be supplied at the time of creation (they must all
+    be the same at the time of initialization, though, or else validation will fail).
+
+    **Examples**::
+
+        # Generate an average of 10 contacts for 1000 people
+        n = 10_000
+        n_people = 1000
+        p1 = np.random.randint(n_people, size=n)
+        p2 = np.random.randint(n_people, size=n)
+        beta = np.ones(n)
+        network = ss.Network(p1=p1, p2=p2, beta=beta, label='rand')
+        network = ss.Network(dict(p1=p1, p2=p2, beta=beta), label='rand') # Alternate method
+
+        # Convert one network to another with extra columns
+        index = np.arange(n)
+        self_conn = p1 == p2
+        network2 = ss.Network(**network, index=index, self_conn=self_conn, label=network.label)
+    """
+
+    def __init__(self, *args, key_dict=None, transmission='horizontal', label=None, **kwargs):
+        default_keys = {
+            'p1': sss.default_int,
+            'p2': sss.default_int,
+            'beta': sss.default_float,
+        }
+        self.meta = sc.mergedicts(default_keys, key_dict)
+        self.transmission = transmission  # "vertical" or "horizontal", determines whether transmission is bidirectional
+        self.basekey = 'p1'  # Assign a base key for calculating lengths and performing other operations
+        self.label = label
+        self.initialized = False
+
+        # Handle args
+        kwargs = sc.mergedicts(*args, kwargs)
+
+        # Initialize the keys of the network
+        for key, dtype in self.meta.items():
+            self[key] = np.empty((0,), dtype=dtype)
+
+        # Set data, if provided
+        for key, value in kwargs.items():
+            self[key] = np.array(value, dtype=self.meta.get(key))
+            self.initialized = True
+
+    @property
+    def name(self):
+        # The module name is a lower-case version of its class name
+        return self.__class__.__name__.lower()
+
+    def initialize(self):
+        pass
+
+    def __len__(self):
+        try:
+            return len(self[self.basekey])
+        except:  # pragma: no cover
+            return 0
+
+    def __repr__(self):
+        """ Convert to a dataframe for printing """
+        namestr = self.__class__.__name__
+        labelstr = f'"{self.label}"' if self.label else '<no label>'
+        keys_str = ', '.join(self.keys())
+        output = f'{namestr}({labelstr}, {keys_str})\n'  # e.g. Network("r", f, m, beta)
+        output += self.to_df().__repr__()
+        return output
+
+    def __contains__(self, item):
+        """
+        Check if a person is present in a network
+
+        Args:
+            item: Person index
+
+        Returns: True if person index appears in any interactions
+        """
+        return (item in self['p1']) or (item in self['p2'])
+
+    @property
+    def members(self):
+        """
+        Return sorted array of all members
+        """
+        return np.unique([self['p1'], self['p2']])
+
+    def meta_keys(self):
+        """ Return the keys for the network's meta information """
+        return self.meta.keys()
+
+    def validate(self, force=True):
+        """
+        Check the integrity of the network: right types, right lengths.
+
+        If dtype is incorrect, try to convert automatically; if length is incorrect,
+        do not.
+        """
+        n = len(self[self.basekey])
+        for key, dtype in self.meta.items():
+            if dtype:
+                actual = self[key].dtype
+                expected = dtype
+                if actual != expected:
+                    self[key] = np.array(self[key],
+                                         dtype=expected)  # Probably harmless, so try to convert to correct type
+            actual_n = len(self[key])
+            if n != actual_n:
+                errormsg = f'Expecting length {n} for network key "{key}"; got {actual_n}'  # We can't fix length mismatches
+                raise TypeError(errormsg)
+        return
+
+    def get_inds(self, inds, remove=False):
+        """
+        Get the specified indices from the edgelist and return them as a dict.
+        Args:
+            inds (int, array, slice): the indices to find
+            remove (bool): whether to remove the indices
+        """
+        output = {}
+        for key in self.meta_keys():
+            output[key] = self[key][inds]  # Copy to the output object
+            if remove:
+                self[key] = np.delete(self[key], inds)  # Remove from the original
+        return output
+
+    def pop_inds(self, inds):
+        """
+        "Pop" the specified indices from the edgelist and return them as a dict.
+        Returns arguments in the right format to be used with network.append().
+
+        Args:
+            inds (int, array, slice): the indices to be removed
+        """
+        return self.get_inds(inds, remove=True)
+
+    def append(self, contacts):
+        """
+        Append contacts to the current network.
+
+        Args:
+            contacts (dict): a dictionary of arrays with keys f,m,beta, as returned from network.pop_inds()
+        """
+        for key in self.keys():
+            new_arr = contacts[key]
+            n_curr = len(self[key])  # Current number of contacts
+            n_new = len(new_arr)  # New contacts to add
+            n_total = n_curr + n_new  # New size
+            self[key] = np.resize(self[key], n_total)  # Resize to make room, preserving dtype
+            self[key][n_curr:] = new_arr  # Copy contacts into the network
+        return
+
+    def to_df(self):
+        """ Convert to dataframe """
+        df = pd.DataFrame.from_dict(self)
+        return df
+
+    def from_df(self, df, keys=None):
+        """ Convert from a dataframe """
+        if keys is None:
+            keys = self.meta_keys()
+        for key in keys:
+            self[key] = df[key].to_numpy()
+        return self
+
+    def find_contacts(self, inds, as_array=True):
+        """
+        Find all contacts of the specified people
+
+        For some purposes (e.g. contact tracing) it's necessary to find all the contacts
+        associated with a subset of the people in this network. Since contacts are bidirectional
+        it's necessary to check both P1 and P2 for the target indices. The return type is a Set
+        so that there is no duplication of indices (otherwise if the Network has explicit
+        symmetric interactions, they could appear multiple times). This is also for performance so
+        that the calling code doesn't need to perform its own unique() operation. Note that
+        this cannot be used for cases where multiple connections count differently than a single
+        infection, e.g. exposure risk.
+
+        Args:
+            inds (array): indices of people whose contacts to return
+            as_array (bool): if true, return as sorted array (otherwise, return as unsorted set)
+
+        Returns:
+            contact_inds (array): a set of indices for pairing partners
+
+        Example: If there were a network with
+        - P1 = [1,2,3,4]
+        - P2 = [2,3,1,4]
+        Then find_contacts([1,3]) would return {1,2,3}
+        """
+
+        # Check types
+        if not isinstance(inds, np.ndarray):
+            inds = sc.promotetoarray(inds)
+        if inds.dtype != np.int64:  # pragma: no cover # This is int64 since indices often come from hpv.true(), which returns int64
+            inds = np.array(inds, dtype=np.int64)
+
+        # Find the contacts
+        contact_inds = ssu.find_contacts(self['p1'], self['p2'], inds)
+        if as_array:
+            contact_inds = np.fromiter(contact_inds, dtype=ssd.default_int)
+            contact_inds.sort()  # Sorting ensures that the results are reproducible for a given seed as well as being identical to previous versions of HPVsim
+
+        return contact_inds
+
+    def add_pairs(self):
+        pass
+
+    def update(self):
+        pass
diff --git a/stisim/misc.py b/stisim/misc.py
new file mode 100644
index 00000000..7e0570ae
--- /dev/null
+++ b/stisim/misc.py
@@ -0,0 +1,303 @@
+"""
+Miscellaneous functions that do not belong anywhere else
+"""
+
+import re
+import inspect
+import warnings
+import numpy as np
+import sciris as sc
+import collections as co
+from . import version as ss
+from .settings import options as sso
+
+# %% Loading/saving functions
+
+__all__ = ['load', 'save']
+
+
+def load(*args, **kwargs):
+    """
+    Convenience method for sc.loadobj() which also compares versions
+
+    Args:
+        args (list): passed to sc.loadobj()
+        kwargs (dict): passed to sc.loadobj()
+
+    Returns:
+        Loaded object
+    """
+    obj = sc.loadobj(*args, **kwargs)
+    if hasattr(obj, 'version'):
+        v_curr = ss.__version__
+        v_obj = obj.version
+        cmp = check_version(v_obj, verbose=False)
+        if cmp != 0:
+            print(f'Note: you have STIsim v{v_curr}, but are loading an object from v{v_obj}')
+    return obj
+
+
+def save(*args, **kwargs):
+    """
+    Convenience method for sc.saveobj()
+    Args:
+        args (list): passed to sc.saveobj()
+        kwargs (dict): passed to sc.saveobj()
+
+    Returns:
+        Filename the object is saved to
+    """
+    filepath = sc.saveobj(*args, **kwargs)
+    return filepath
+
+
+# %% Versioning functions
+
+__all__ += ['git_info', 'check_version', 'check_save_version']
+
+
+def git_info(filename=None, check=False, comments=None, old_info=None, die=False, indent=2, verbose=True, frame=2,
+             **kwargs):
+    """
+    Get current git information and optionally write it to disk. Simplest usage
+    is ss.git_info(__file__)
+
+    Args:
+        filename  (str): name of the file to write to or read from
+        check    (bool): whether to compare two git versions
+        comments (dict): additional comments to include in the file
+        old_info (dict): dictionary of information to check against
+        die      (bool): whether to raise an exception if the check fails
+        indent    (int): how many indents to use when writing the file to disk
+        verbose  (bool): detail to print
+        frame     (int): how many frames back to look for caller info
+        kwargs   (dict): passed to sc.loadjson() (if check=True) or sc.savejson() (if check=False)
+    """
+
+    # Handle the case where __file__ is supplied as the argument
+    if isinstance(filename, str) and filename.endswith('.py'):
+        filename = filename.replace('.py', '.gitinfo')
+
+    # Get git info
+    calling_file = sc.makefilepath(sc.getcaller(frame=frame, tostring=False)['filename'])
+    ss_info = {'version': ss.__version__}
+    ss_info.update(sc.gitinfo(__file__, verbose=False))
+    caller_info = sc.gitinfo(calling_file, verbose=False)
+    caller_info['filename'] = calling_file
+    info = {'starsim': ss_info, 'called_by': caller_info}
+    if comments:
+        info['comments'] = comments
+
+    # Just get information and optionally write to disk
+    if not check:
+        if filename is not None:
+            output = sc.savejson(filename, info, indent=indent, **kwargs)
+        else:
+            output = info
+        return output
+
+    # Check if versions match, and optionally raise an error
+    else:
+        if filename is not None:
+            old_info = sc.loadjson(filename, **kwargs)
+        old_ss_info = old_info['stisim'] if 'stisim' in old_info else old_info
+        if ss_info != old_ss_info:  # pragma: no cover
+            string = f'Git information differs: {ss_info} vs. {old_ss_info}'
+            if die:
+                raise ValueError(string)
+            elif verbose:
+                print(string)
+        return
+
+
+def check_version(expected, die=False, verbose=True):
+    """
+    Get current git information and optionally write it to disk. The expected
+    version string may optionally start with '>=' or '<=' (== is implied otherwise),
+    but other operators (e.g. ~=) are not supported. Note that e.g. '>' is interpreted
+    to mean '>='.
+
+    Args:
+        expected (str): expected version information
+        die (bool): whether to raise an exception if the check fails
+        verbose  (bool): detail to print
+    """
+    if expected.startswith('>'):
+        valid = 1
+    elif expected.startswith('<'):
+        valid = -1
+    else:
+        valid = 0  # Assume == is the only valid comparison
+    expected = expected.lstrip('<=>')  # Remove comparator information
+    version = ss.__version__
+    compare = sc.compareversions(version, expected)  # Returns -1, 0, or 1
+    relation = ['older', '', 'newer'][compare + 1]  # Picks the right string
+    if relation:  # Versions mismatch, print warning or raise error
+        string = f'Note: STIsim is {relation} than expected ({version} vs. {expected})'
+        if die and compare != valid:
+            raise ValueError(string)
+        elif verbose:
+            print(string)
+    return compare
+
+
+def check_save_version(expected=None, filename=None, die=False, verbose=True, **kwargs):
+    """
+    A convenience function that bundles check_version with git_info and saves
+    automatically to disk from the calling file. The idea is to put this at the
+    top of an analysis script, and commit the resulting file, to keep track of
+    which version of STIsim was used.
+
+    Args:
+        expected (str): expected version information
+        filename (str): file to save to; if None, guess based on current file name
+        die (bool): whether to raise an exception if the check fails
+        verbose  (bool): detail to print
+        kwargs (dict): passed to git_info(), and thence to sc.savejson()
+    """
+
+    # First, check the version if supplied
+    if expected:
+        check_version(expected, die=die, verbose=verbose)
+
+    # Now, check and save the git info
+    if filename is None:
+        filename = sc.getcaller(tostring=False)['filename']
+    git_info(filename=filename, frame=3, **kwargs)
+
+    return
+
+
+# %% Help and warnings
+
+__all__ += ['help']
+
+
+def help(pattern=None, source=False, ignorecase=True, flags=None, context=False, output=False):
+    """
+    Get help on STIsim in general, or search for a word/expression.
+
+    Args:
+        pattern    (str):  the word, phrase, or regex to search for
+        source     (bool): whether to search source code instead of docstrings for matches
+        ignorecase (bool): whether to ignore case (equivalent to ``flags=re.I``)
+        flags      (list): additional flags to pass to ``re.findall()``
+        context    (bool): whether to show the line(s) of matches
+        output     (bool): whether to return the dictionary of matches
+
+    **Examples**::
+
+        ss.help()
+        ss.help('vaccine')
+        ss.help('contact', ignorecase=False, context=True)
+        ss.help('lognormal', source=True, context=True)
+    """
+    defaultmsg = '''Help is coming.'''
+    # No pattern is provided, print out default help message
+    if pattern is None:
+        print(defaultmsg)
+
+    else:
+        import stisim as ss  # Here to avoid circular import
+
+        # Handle inputs
+        flags = sc.promotetolist(flags)
+        if ignorecase:
+            flags.append(re.I)
+
+        def func_ok(fucname, func):
+            ''' Skip certain functions '''
+            excludes = [
+                fucname.startswith('_'),
+                fucname in ['help', 'options', 'default_float', 'default_int'],
+                inspect.ismodule(func),
+            ]
+            ok = not (any(excludes))
+            return ok
+
+        # Get available functions/classes
+        funcs = [funcname for funcname in dir(ss) if
+                 func_ok(funcname, getattr(ss, funcname))]  # Skip dunder methods and modules
+
+        # Get docstrings or full source code
+        docstrings = dict()
+        for funcname in funcs:
+            f = getattr(ss, funcname)
+            if source:
+                string = inspect.getsource(f)
+            else:
+                string = f.__doc__
+            docstrings[funcname] = string
+
+        # Find matches
+        matches = co.defaultdict(list)
+        linenos = co.defaultdict(list)
+
+        for k, docstring in docstrings.items():
+            for l, line in enumerate(docstring.splitlines()):
+                if re.findall(pattern, line, *flags):
+                    linenos[k].append(str(l))
+                    matches[k].append(line)
+
+        # Assemble output
+        if not len(matches):
+            string = f'No matches for "{pattern}" found among {len(docstrings)} available functions.'
+        else:
+            string = f'Found {len(matches)} matches for "{pattern}" among {len(docstrings)} available functions:\n'
+            maxkeylen = 0
+            for k in matches.keys(): maxkeylen = max(len(k), maxkeylen)
+            for k, match in matches.items():
+                if not context:
+                    keystr = f'  {k:>{maxkeylen}s}'
+                else:
+                    keystr = k
+                matchstr = f'{keystr}: {len(match):>2d} matches'
+                if context:
+                    matchstr = sc.heading(matchstr, output=True)
+                else:
+                    matchstr += '\n'
+                string += matchstr
+                if context:
+                    lineno = linenos[k]
+                    maxlnolen = max([len(l) for l in lineno])
+                    for l, m in zip(lineno, match):
+                        string += sc.colorize(string=f'  {l:>{maxlnolen}s}: ', fg='cyan', output=True)
+                        string += f'{m}\n'
+                    string += '—' * 60 + '\n'
+
+        # Print result and return
+        print(string)
+        if output:
+            return string
+        else:
+            return
+
+
+def warn(msg, category=None, verbose=None, die=None):
+    """ Helper function to handle warnings -- not for the user """
+
+    # Handle inputs
+    warnopt = sso.warnings if not die else 'error'
+    if category is None:
+        category = RuntimeWarning
+    if verbose is None:
+        verbose = sso.verbose
+
+    # Handle the different options
+    if warnopt in ['error', 'errors']:  # Include alias since hard to remember
+        raise category(msg)
+    elif warnopt == 'warn':
+        msg = '\n' + msg
+        warnings.warn(msg, category=category, stacklevel=2)
+    elif warnopt == 'print':
+        if verbose:
+            msg = 'Warning: ' + msg
+            print(msg)
+    elif warnopt == 'ignore':
+        pass
+    else:
+        options = ['error', 'warn', 'print', 'ignore']
+        errormsg = f'Could not understand "{warnopt}": should be one of {options}'
+        raise ValueError(errormsg)
+
+    return
diff --git a/stisim/modules.py b/stisim/modules.py
index 955d0490..240daf62 100644
--- a/stisim/modules.py
+++ b/stisim/modules.py
@@ -1,194 +1,317 @@
-import sciris as sc
-import numpy as np
-from .people import State
-from .results import Result
-
-class Module():
-    # Base module contains states/attributes that all modules have
-    default_pars = {}
-
-    def __init__(self):
-        raise Exception('Modules should not be instantiated but their classes should be used directly')
-
-    @classmethod
-    def initialize(cls, sim):
-        # Add this module to a People instance. This would always involve calling People.add_module
-        # but subsequently modules could have their own logic for initializing the default values
-        # and initializing any outputs that are required
-        sim.people.add_module(cls)
-
-        # Merge parameters
-        if cls.name not in sim.pars:
-            sim.pars[cls.name] = sc.objdict(sc.dcp(cls.default_pars))
-        else:
-            if ~isinstance(sim.pars[cls.name], sc.objdict):
-                sim.pars[cls.name] = sc.objdict(sim.pars[cls.name])
-            for k,v in cls.default_pars.items():
-                if k not in sim.pars[cls.name]:
-                    sim.pars[cls.name][k] = v
-
-        sim.results[cls.name] = sc.objdict()
-
-    @classmethod
-    def update_states_pre(cld, sim):
-        # Carry out any autonomous state changes at the start of the timestep
-        pass
-
-    @classmethod
-    def update_results(cls, sim):
-        pass
-
-    @classmethod
-    def finalize_results(cls, sim):
-        pass
-
-    @classmethod
-    @property
-    def name(cls):
-        # The module name is a lower-case version of its class name
-        return cls.__name__.lower()
-
-
-class HIV(Module):
-    states = [
-        State('susceptible', bool, True),
-        State('infected', bool, False),
-        State('ti_infected', float, 0),
-        State('on_art', bool, False),
-        State("cd4", float, 500),
-    ]
-
-    default_pars = {
-        'cd4_min': 100,
-        'cd4_max': 500,
-        'cd4_rate': 5,
-        'initial': 3,
-    }
-
-    @classmethod
-    def update_states_pre(cls, sim):
-        # Update CD4
-        sim.people.hiv.cd4[~sim.people.dead & sim.people.hiv.infected & sim.people.hiv.on_art] += (sim.pars.hiv.cd4_max - sim.people.hiv.cd4[~sim.people.dead & sim.people.hiv.infected & sim.people.hiv.on_art])/sim.pars.hiv.cd4_rate
-        sim.people.hiv.cd4[~sim.people.dead & sim.people.hiv.infected & ~sim.people.hiv.on_art] += (sim.pars.hiv.cd4_min - sim.people.hiv.cd4[~sim.people.dead & sim.people.hiv.infected & ~sim.people.hiv.on_art])/sim.pars.hiv.cd4_rate
-
-    @classmethod
-    def initialize(cls, sim):
-        super(HIV, cls).initialize(sim)
-
-        # Pick some initially infected agents
-        cls.infect(sim, np.random.choice(sim.people.uid, sim.pars[cls.name]['initial']))
-
-        if 'beta' not in sim.pars[cls.name]:
-            sim.pars[cls.name].beta = sc.objdict({k:1 for k in sim.people.contacts})
-
-        sim.results[cls.name]['n_susceptible'] = Result(cls.name, 'n_susceptible', sim.npts, dtype=int)
-        sim.results[cls.name]['n_infected'] = Result(cls.name, 'n_infected', sim.npts, dtype=int)
-        sim.results[cls.name]['prevalence'] = Result(cls.name, 'prevalence', sim.npts, dtype=float)
-        sim.results[cls.name]['new_infections'] = Result(cls.name, 'n_infected', sim.npts, dtype=int)
-        sim.results[cls.name]['n_art'] = Result(cls.name, 'n_art', sim.npts, dtype=int)
-
-    @classmethod
-    def update_results(cls, sim):
-        sim.results[cls.name]['n_susceptible'][sim.ti] = np.count_nonzero(sim.people.hiv.susceptible)
-        sim.results[cls.name]['n_infected'][sim.ti] = np.count_nonzero(sim.people.hiv.infected)
-        sim.results[cls.name]['prevalence'][sim.ti] = sim.results[cls.name].n_infected[sim.ti] / sim.people.n
-        sim.results[cls.name]['new_infections'] = np.count_nonzero(sim.people[cls.name].ti_infected == sim.ti)
-        sim.results[cls.name]['n_art'] = np.count_nonzero(~sim.people.dead & sim.people[cls.name].on_art)
-
-    @classmethod
-    def transmit(cls, sim):
-        for k, layer in sim.people.contacts.items():
-            if k in sim.pars[cls.name]['beta']:
-                rel_trans = (sim.people[cls.name].infected & ~sim.people.dead).astype(float)
-                rel_sus = (sim.people[cls.name].susceptible & ~sim.people.dead).astype(float)
-                for a,b in [[layer.p1,layer.p2],[layer.p2,layer.p1]]:
-                    # probability of a->b transmission
-                    p_transmit = rel_trans[a]*rel_sus[b]*layer.beta*sim.pars[cls.name]['beta'][k]
-                    cls.infect(sim, b[np.random.random(len(a))<p_transmit])
-
-    @classmethod
-    def infect(cls, sim, uids):
-        sim.people[cls.name].susceptible[uids] = False
-        sim.people[cls.name].infected[uids] = True
-        sim.people[cls.name].ti_infected[uids] = sim.ti
-
-
-class Gonorrhea(Module):
-    states = [
-        State('susceptible', bool, True),
-        State('infected', bool, False),
-        State('ti_infected', float, 0),
-        State('ti_recovered', float, 0),
-        State('ti_dead', float, np.nan), # Death due to gonorrhea
-    ]
-
-    default_pars = {
-        'dur_inf': 3, # not modelling diagnosis or treatment explicitly here
-        'p_death': 0.2,
-        'initial': 3,
-    }
-
-    def __init__(self, pars):
-        super().__init__(pars)
-
-    @classmethod
-    def update_states_pre(cls, sim):
-        # What if something in here should depend on another module?
-        # I guess we could just check for it e.g., 'if HIV in sim.modules' or
-        # 'if 'hiv' in sim.people' or something
-        gonorrhea_deaths = sim.people.gonorrhea.ti_dead <= sim.ti
-        sim.people.dead[gonorrhea_deaths] = True
-        sim.people.ti_dead[gonorrhea_deaths] = sim.ti
-
-    @classmethod
-    def initialize(cls, sim):
-        # Do any steps in this method depend on what other modules are going to be added? We can inspect
-        # them via sim.modules at this point
-        super(Gonorrhea, cls).initialize(sim)
-
-        # Pick some initially infected agents
-        cls.infect(sim, np.random.choice(sim.people.uid, sim.pars[cls.name]['initial']))
-
-        if 'beta' not in sim.pars[cls.name]:
-            sim.pars[cls.name].beta = sc.objdict({k:1 for k in sim.people.contacts})
-
-        # TODO - not a huge fan of having the result key duplicate the Result name
-        sim.results[cls.name]['n_susceptible'] = Result(cls.name, 'n_susceptible', sim.npts, dtype=int)
-        sim.results[cls.name]['n_infected'] = Result(cls.name, 'n_infected', sim.npts, dtype=int)
-        sim.results[cls.name]['prevalence'] = Result(cls.name, 'prevalence', sim.npts, dtype=float)
-        sim.results[cls.name]['new_infections'] = Result(cls.name, 'n_infected', sim.npts, dtype=int)
-
-
-    @classmethod
-    def update_results(cls, sim):
-        sim.results[cls.name]['n_susceptible'][sim.ti] = np.count_nonzero(sim.people.gonorrhea.susceptible)
-        sim.results[cls.name]['n_infected'][sim.ti] = np.count_nonzero(sim.people.gonorrhea.infected)
-        sim.results[cls.name]['prevalence'][sim.ti] = sim.results[cls.name].n_infected[sim.ti] / sim.people.n
-        sim.results[cls.name]['new_infections'] = np.count_nonzero(sim.people[cls.name].ti_infected == sim.ti)
-
-    @classmethod
-    def transmit(cls, sim):
-        for k, layer in sim.people.contacts.items():
-            if k in sim.pars[cls.name]['beta']:
-                rel_trans = (sim.people[cls.name].infected & ~sim.people.dead).astype(float)
-                rel_sus = (sim.people[cls.name].susceptible & ~sim.people.dead).astype(float)
-                for a,b in [[layer.p1,layer.p2],[layer.p2,layer.p1]]:
-                    # probability of a->b transmission
-                    p_transmit = rel_trans[a]*rel_sus[b]*layer.beta*sim.pars[cls.name]['beta'][k]
-                    cls.infect(sim, b[np.random.random(len(a))<p_transmit])
-
-    @classmethod
-    def infect(cls, sim, uids):
-        sim.people[cls.name].susceptible[uids] = False
-        sim.people[cls.name].infected[uids] = True
-        sim.people[cls.name].ti_infected[uids] = sim.ti
-
-        dur = sim.ti+np.random.poisson(sim.pars[cls.name]['dur_inf']/sim.pars.dt, len(uids))
-        dead = np.random.random(len(uids))<sim.pars[cls.name].p_death
-        sim.people[cls.name].ti_recovered[uids[~dead]] = dur[~dead]
-        sim.people[cls.name].ti_dead[uids[dead]] = dur[dead]
-
-        # if HIV in sim.modules:
-        #     # check CD4 count
-        #     # increase susceptibility where relevant
-
+import sciris as sc
+import numpy as np
+from .base import State
+from .results import Result
+from . import utils as ssu
+from . import networks as ssnet
+
+
+class Module(sc.prettyobj):
+    # Base module contains states/attributes that all modules have
+    
+    def __init__(self, pars=None, *args, **kwargs):
+        super().__init__(*args, **kwargs)
+        self.pars = ssu.omerge(pars)
+        self.states = ssu.named_dict(
+            State('rel_sus', float, 1),
+            State('rel_sev', float, 1),
+            State('rel_trans', float, 1),
+        )
+        self.results = sc.objdict()
+        return
+
+    def initialize(self, sim):
+        # Merge parameters
+        sim.pars[self.name] = self.pars
+        sim.results[self.name] = self.results
+
+        # Add this module to a People instance. This would always involve calling People.add_module
+        # but subsequently modules could have their own logic for initializing the default values
+        # and initializing any outputs that are required
+        sim.people.add_module(self)
+
+        # Pick some initially infected agents
+        self.set_prognoses(sim, np.random.choice(sim.people.uid, sim.pars[self.name]['initial']))
+
+        # Validate pars
+        if 'beta' not in sim.pars[self.name]:
+            sim.pars[self.name].beta = sc.objdict({k: [1, 1] for k in sim.people.networks})
+
+        # Initialize results
+        sim.results[self.name]['n_susceptible'] = Result(self.name, 'n_susceptible', sim.npts, dtype=int)
+        sim.results[self.name]['n_infected'] = Result(self.name, 'n_infected', sim.npts, dtype=int)
+        sim.results[self.name]['prevalence'] = Result(self.name, 'prevalence', sim.npts, dtype=float)
+        sim.results[self.name]['new_infections'] = Result(self.name, 'n_infected', sim.npts, dtype=int)
+
+    def update_states(self, sim):
+        # Carry out any autonomous state changes at the start of the timestep
+        pass
+
+    def make_new_cases(self, sim):
+        """ Add new cases of module, through transmission, incidence, etc. """
+        pars = sim.pars[self.name]
+        for k, layer in sim.people.networks.items():
+            if k in pars['beta']:
+                rel_trans = (sim.people[self.name].infected & ~sim.people.dead).astype(float)
+                rel_sus = (sim.people[self.name].susceptible & ~sim.people.dead).astype(float)
+                for a, b, beta in [[layer['p1'], layer['p2'], pars['beta'][k][0]], [layer['p2'], layer['p1'], pars['beta'][k][1]]]:
+                    # probability of a->b transmission
+                    p_transmit = rel_trans[a] * rel_sus[b] * layer['beta'] * beta
+                    new_cases = np.random.random(len(a)) < p_transmit
+                    if new_cases.any():
+                        self.set_prognoses(sim, b[new_cases])
+
+    def set_prognoses(self, sim, uids):
+        pass
+
+    def update_results(self, sim):
+        sim.results[self.name]['n_susceptible'][sim.ti] = np.count_nonzero(sim.people[self.name].susceptible)
+        sim.results[self.name]['n_infected'][sim.ti] = np.count_nonzero(sim.people[self.name].infected)
+        sim.results[self.name]['prevalence'][sim.ti] = sim.results[self.name].n_infected[sim.ti] / sim.people.n
+        sim.results[self.name]['new_infections'] = np.count_nonzero(sim.people[self.name].ti_infected == sim.ti)
+
+    def finalize_results(self, sim):
+        pass
+
+    @property
+    def name(self):
+        # The module name is a lower-case version of its class name
+        return self.__class__.__name__.lower()
+
+
+class HIV(Module):
+    
+    def __init__(self, pars=None):
+        super().__init__(pars)
+        self.states = ssu.omerge(ssu.named_dict(
+            State('susceptible', bool, True),
+            State('infected', bool, False),
+            State('ti_infected', float, 0),
+            State('on_art', bool, False),
+            State('cd4', float, 500),
+        ), self.states)
+    
+        self.pars = ssu.omerge({
+            'cd4_min': 100,
+            'cd4_max': 500,
+            'cd4_rate': 5,
+            'initial': 30,
+            'eff_condoms': 0.7,
+        }, self.pars)
+        return
+
+    def update_states(self, sim):
+        # Update CD4
+        hivppl = sim.people.hiv
+        hivppl.cd4[sim.people.alive & hivppl.infected & hivppl.on_art] += (sim.pars.hiv.cd4_max - hivppl.cd4[sim.people.alive & hivppl.infected & hivppl.on_art])/sim.pars.hiv.cd4_rate
+        hivppl.cd4[sim.people.alive & hivppl.infected & ~hivppl.on_art] += (sim.pars.hiv.cd4_min - hivppl.cd4[sim.people.alive & hivppl.infected & ~hivppl.on_art])/sim.pars.hiv.cd4_rate
+        return
+
+    def initialize(self, sim):
+        Module.initialize(self, sim)
+        sim.results[self.name]['n_art'] = Result('n_art', self.name, sim.npts, dtype=int)
+    
+    def update_results(self, sim):
+        super(HIV, self).update_results(sim)
+        sim.results[self.name]['n_art'] = np.count_nonzero(sim.people.alive & sim.people[self.name].on_art)
+
+    def make_new_cases(self, sim):
+        # eff_condoms = sim.pars[self.name]['eff_condoms'] # TODO figure out how to add this
+        super().make_new_cases(sim)
+    
+    def set_prognoses(self, sim, uids):
+        sim.people[self.name].susceptible[uids] = False
+        sim.people[self.name].infected[uids] = True
+        sim.people[self.name].ti_infected[uids] = sim.ti
+
+
+class Gonorrhea(Module):
+
+    def __init__(self, pars=None):
+        super().__init__(pars)
+        self.states = ssu.omerge(ssu.named_dict(
+            State('susceptible', bool, True),
+            State('infected', bool, False),
+            State('ti_infected', float, 0),
+            State('ti_recovered', float, 0),
+            State('ti_dead', float, np.nan), # Death due to gonorrhea
+        ), self.states)
+
+        self.pars = ssu.omerge({
+            'dur_inf': 3, # not modelling diagnosis or treatment explicitly here
+            'p_death': 0.2,
+            'initial': 3,
+            'eff_condoms': 0.7,
+        }, self.pars)
+        return
+    
+    def update_states(self, sim):
+        # What if something in here should depend on another module?
+        # I guess we could just check for it e.g., 'if HIV in sim.modules' or
+        # 'if 'hiv' in sim.people' or something
+        gonorrhea_deaths = sim.people.gonorrhea.ti_dead <= sim.ti
+        sim.people.alive[gonorrhea_deaths] = False
+        sim.people.ti_dead[gonorrhea_deaths] = sim.ti
+        return
+
+    def initialize(self, sim):
+        # Do any steps in this method depend on what other modules are going to be added? We can inspect
+        # them via sim.modules at this point
+        Module.initialize(self, sim)
+    
+    def update_results(self, sim):
+        super(Gonorrhea, self).update_results(sim)
+    
+    def make_new_cases(self, sim):
+        super(Gonorrhea, self).make_new_cases(sim)
+
+    def set_prognoses(self, sim, uids):
+        sim.people[self.name].susceptible[uids] = False
+        sim.people[self.name].infected[uids] = True
+        sim.people[self.name].ti_infected[uids] = sim.ti
+
+        dur = sim.ti + np.random.poisson(sim.pars[self.name]['dur_inf']/sim.pars.dt, len(uids))
+        dead = np.random.random(len(uids)) < sim.pars[self.name].p_death
+        sim.people[self.name].ti_recovered[uids[~dead]] = dur[~dead]
+        sim.people[self.name].ti_dead[uids[dead]] = dur[dead]
+
+
+class Pregnancy(Module):
+
+    def __init__(self, pars=None):
+        super().__init__(pars)
+
+
+        # Other, e.g. postpartum, on contraception...
+        self.states = ssu.omerge(ssu.named_dict(
+            State('infertile', bool, False),  # Applies to girls and women outside the fertility window
+            State('susceptible', bool, True),  # Applies to girls and women inside the fertility window - needs renaming
+            State('pregnant', bool, False),  # Currently pregnant
+            State('postpartum', bool, False),  # Currently post-partum
+            State('ti_pregnant', float, np.nan),  # Time pregnancy begins
+            State('ti_delivery', float, np.nan),  # Time of delivery
+            State('ti_postpartum', float, np.nan),  # Time postpartum ends
+            State('ti_dead', float, np.nan),  # Maternal mortality
+        ), self.states)
+
+        self.pars = ssu.omerge({
+            'dur_pregnancy': 0.75,  # Make this a distribution?
+            'dur_postpartum': 0.5,  # Make this a distribution?
+            'inci': 0.03,  # Replace this with age-specific rates
+            'p_death': 0.02,  # Probability of maternal death. Question, should this be linked to age and/or duration?
+            'initial': 3,  # Number of women initially pregnant
+        }, self.pars)
+
+        return
+
+    def initialize(self, sim):
+        """
+        Results could include a range of birth outcomes e.g. LGA, stillbirths, etc.
+        Still unclear whether this logic should live in the pregnancy module, the
+        individual disease modules, the connectors, or the sim.
+        """
+        Module.initialize(self, sim)
+        sim.results[self.name]['pregnancies'] = Result('pregnancies', self.name, sim.npts, dtype=int)
+        sim.results[self.name]['births'] = Result('births', self.name, sim.npts, dtype=int)
+        sim['birth_rates'] = None  # This turns off birth rate pars so births only come from this module
+        return
+
+    def update_states(self, sim):
+
+        # Check for new deliveries
+        deliveries = sim.people[self.name].pregnant & (sim.people[self.name].ti_delivery <= sim.ti)
+        sim.people[self.name].pregnant[deliveries] = False
+        sim.people[self.name].postpartum[deliveries] = True
+        sim.people[self.name].susceptible[deliveries] = False
+        sim.people[self.name].ti_delivery[deliveries] = sim.ti
+
+        # Check for new kids emerging from post-partum
+        postpartum = ~sim.people[self.name].pregnant & (sim.people[self.name].ti_postpartum <= sim.ti)
+        sim.people[self.name].postpartum[postpartum] = False
+        sim.people[self.name].susceptible[postpartum] = True
+        sim.people[self.name].ti_postpartum[postpartum] = sim.ti
+
+        delivery_inds = ssu.true(deliveries)
+        if len(delivery_inds):
+            for _, layer in sim.people.networks.items():
+                if isinstance(layer, ssnet.maternal):
+                    # new_birth_inds = layer.find_contacts(delivery_inds)  # Don't think we need this?
+                    new_births = len(delivery_inds) * sim['pop_scale']
+                    # sim.people.demographic_flows['births'] = new_births
+
+        # Maternal deaths
+        maternal_deaths = ssu.true(sim.people[self.name].ti_dead <= sim.ti)
+        if len(maternal_deaths):
+            sim.people.alive[maternal_deaths] = False
+            sim.people.ti_dead[maternal_deaths] = sim.ti
+
+        return
+
+    def make_new_cases(self, sim):
+        """
+        Select people to make pregnancy using incidence data
+        This should use ASFR data from https://population.un.org/wpp/Download/Standard/Fertility/
+        """
+        # Abbreviate key variables
+        cpars = sim.pars[self.name]
+        ppl = sim.people
+        this_inci = cpars['inci']
+
+        # If incidence of pregnancy is non-zero, make some cases
+        # Think about how to deal with age/time-varying fertility
+        if this_inci > 0:
+            demon_conds = ppl.female & ppl.active & ppl[self.name].susceptible
+            inds_to_choose_from = ssu.true(demon_conds)
+            uids = ssu.binomial_filter(this_inci, inds_to_choose_from)
+
+            # Add UIDs for the as-yet-unborn agents so that we can track prognoses and transmission patterns
+            n_unborn_agents = len(uids)
+            if n_unborn_agents > 0:
+                # Grow the arrays and set properties for the unborn agents
+                new_inds = sim.people._grow(n_unborn_agents)
+                sim.people.uid[new_inds] = new_inds
+                sim.people.age[new_inds] = -cpars['dur_pregnancy']
+                sim.people.female[new_inds] = np.random.choice([True, False], size=n_unborn_agents)
+
+                # Add connections to any vertical transmission layers
+                # Placeholder code to be moved / refactored. The maternal network may need to be
+                # handled separately to the sexual networks, TBC how to handle this most elegantly
+                for lkey, layer in sim.people.networks.items():
+                    if layer.transmission == 'vertical':  # What happens if there's more than one vertical layer?
+                        durs = np.full(n_unborn_agents, fill_value=cpars['dur_pregnancy']+cpars['dur_postpartum'])
+                        layer.add_pairs(uids, new_inds, dur=durs)
+
+                # Set prognoses for the pregnancies
+                self.set_prognoses(sim, uids)
+
+        return
+
+    def set_prognoses(self, sim, uids):
+        """
+        Make pregnancies
+        Add miscarriage/termination logic here
+        Also reconciliation with birth rates
+        Q, is this also a good place to check for other conditions and set prognoses for the fetus?
+        """
+        pars = sim.pars[self.name]
+
+        # Change states for the newly pregnant woman
+        sim.people[self.name].susceptible[uids] = False
+        sim.people[self.name].pregnant[uids] = True
+        sim.people[self.name].ti_pregnant[uids] = sim.ti
+
+        # Outcomes for pregnancies
+        dur = np.full(len(uids), sim.ti + pars['dur_pregnancy'] / sim.pars.dt)
+        dead = np.random.random(len(uids)) < sim.pars[self.name].p_death
+        sim.people[self.name].ti_delivery[uids] = dur  # Currently assumes maternal deaths still result in a live baby
+        dur_post_partum = np.full(len(uids), dur + pars['dur_postpartum'] / sim.dt)
+        sim.people[self.name].ti_postpartum[uids] = dur_post_partum
+
+        if len(ssu.true(dead)):
+            sim.people[self.name].ti_dead[uids[dead]] = dur[dead]
+        return
+
+    def update_results(self, sim):
+        mppl = sim.people[self.name]
+        sim.results[self.name]['pregnancies'][sim.ti] = np.count_nonzero(mppl.ti_pregnant == sim.ti)
+        sim.results[self.name]['births'][sim.ti] = np.count_nonzero(mppl.ti_delivery == sim.ti)
diff --git a/stisim/networks.py b/stisim/networks.py
new file mode 100644
index 00000000..9e2ec330
--- /dev/null
+++ b/stisim/networks.py
@@ -0,0 +1,312 @@
+'''
+Networks that connect people within a population
+'''
+
+# %% Imports
+import numpy as np
+import sciris as sc
+from . import utils as ssu
+from . import base as ssb
+from . import settings as sss
+
+
+# Specify all externally visible functions this file defines
+__all__ = ['simple_sexual', 'hpv_network', 'maternal']
+
+
+class simple_sexual(ssb.Network):
+    # Randomly pair males and females with variable relationship durations
+    def __init__(self, mean_dur=5):
+        key_dict = {
+            'p1': sss.default_int,
+            'p2': sss.default_int,
+            'dur': sss.default_float,
+            'beta': sss.default_float,
+        }
+
+        # Call init for the base class, which sets all the keys
+        super().__init__(key_dict=key_dict)
+
+        # Set other parameters
+        self.mean_dur = mean_dur
+
+    def initialize(self, people):
+        self.add_pairs(people, ti=0)
+
+    def add_pairs(self, people, ti=None):
+        # Find unpartnered males and females - could in principle check other contact layers too
+        # by having the People object passed in here
+
+        available_m = np.setdiff1d(people.uid[~people.female], self.members)
+        available_f = np.setdiff1d(people.uid[people.female], self.members)
+
+        if len(available_m) <= len(available_f):
+            p1 = available_m
+            p2 = np.random.choice(available_f, len(p1), replace=False)
+        else:
+            p2 = available_f
+            p1 = np.random.choice(available_m, len(p2), replace=False)
+
+        beta = np.ones_like(p1)
+        dur = np.random.poisson(self.mean_dur, len(p1))
+        self['p1'] = np.concatenate([self['p1'], p1])
+        self['p2'] = np.concatenate([self['p2'], p2])
+        self['beta'] = np.concatenate([self['beta'], beta])
+        self['dur'] = np.concatenate([self['dur'], dur])
+
+    def update(self, people, dt=None):
+        if dt is None: dt = people.dt
+        # First remove any relationships due to end
+        self['dur'] = self['dur'] - dt
+        active = self['dur'] > 0
+        self['p1'] = self['p1'][active]
+        self['p2'] = self['p2'][active]
+        self['beta'] = self['beta'][active]
+        self['dur'] = self['dur'][active]
+
+        # Then add new relationships for unpartnered people
+        self.add_pairs(people)
+
+
+class hpv_network(ssb.Network):
+    def __init__(self, pars=None):
+
+        key_dict = {
+            'p1': sss.default_int,
+            'p2': sss.default_int,
+            'dur': sss.default_float,
+            'acts': sss.default_float,
+            'start': sss.default_float,
+            'beta': sss.default_float,
+        }
+
+        # Call init for the base class, which sets all the keys
+        super().__init__(key_dict=key_dict)
+
+        # Define default parameters
+        self.pars = dict()
+        self.pars['cross_layer'] = 0.05  # Proportion of agents who have concurrent cross-layer relationships
+        self.pars['partners'] = dict(dist='poisson', par1=0.01)  # The number of concurrent sexual partners
+        self.pars['acts'] = dict(dist='neg_binomial', par1=80, par2=40)  # The number of sexual acts per year
+        self.pars['age_act_pars'] = dict(peak=30, retirement=100, debut_ratio=0.5, retirement_ratio=0.1) # Parameters describing changes in coital frequency over agent lifespans
+        self.pars['condoms'] = 0.2  # The proportion of acts in which condoms are used
+        self.pars['dur_pship'] = dict(dist='normal_pos', par1=1, par2=1)  # Duration of partnerships
+        self.pars['participation'] = None  # Incidence of partnership formation by age
+        self.pars['mixing'] = None  # Mixing matrices for storing age differences in partnerships
+
+        self.update_pars(pars)
+        self.get_layer_probs()
+
+    def initialize(self, people):
+        self.add_pairs(people, ti=0)
+
+    def update_pars(self, pars):
+        if pars is not None:
+            for k, v in pars.items():
+                self.pars[k] = v
+        return
+
+    def get_layer_probs(self):
+
+        defaults = {}
+        mixing = np.array([
+            #       0,  5,  10, 15, 20, 25, 30, 35, 40, 45, 50, 55, 60, 65, 70, 75
+            [ 0,    0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0],
+            [ 5,    0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0],
+            [10,    0,  0, .1,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0],
+            [15,    0,  0, .1, .1,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0],
+            [20,    0,  0, .1, .1, .1, .1,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0],
+            [25,    0,  0, .5, .1, .5, .1, .1,  0,  0,  0,  0,  0,  0,  0,  0,  0],
+            [30,    0,  0,  1, .5, .5, .5, .5, .1,  0,  0,  0,  0,  0,  0,  0,  0],
+            [35,    0,  0, .5,  1,  1, .5,  1,  1, .5,  0,  0,  0,  0,  0,  0,  0],
+            [40,    0,  0,  0, .5,  1,  1,  1,  1,  1, .5,  0,  0,  0,  0,  0,  0],
+            [45,    0,  0,  0,  0, .1,  1,  1,  2,  1,  1, .5,  0,  0,  0,  0,  0],
+            [50,    0,  0,  0,  0,  0, .1,  1,  1,  1,  1,  2, .5,  0,  0,  0,  0],
+            [55,    0,  0,  0,  0,  0,  0, .1,  1,  1,  1,  1,  2, .5,  0,  0,  0],
+            [60,    0,  0,  0,  0,  0,  0,  0, .1, .5,  1,  1,  1,  2, .5,  0,  0],
+            [65,    0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  1,  1,  1,  2, .5,  0],
+            [70,    0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  1,  1,  1,  1, .5],
+            [75,    0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  1,  1,  1,  1],
+        ])
+
+        participation = np.array([
+                [ 0,  5,    10,    15,   20,   25,   30,   35,    40,    45,    50,    55,    60,    65,    70,    75],
+                [ 0,  0,  0.10,   0.7,  0.8,  0.6,  0.6,  0.4,   0.1,  0.05, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001], # Share of females of each age newly having casual relationships
+                [ 0,  0,  0.05,   0.7,  0.8,  0.6,  0.6,  0.4,   0.4,   0.3,   0.1,  0.05,  0.01,  0.01, 0.001, 0.001]], # Share of males of each age newly having casual relationships
+            )
+
+        defaults['mixing'] = mixing
+        defaults['participation'] = participation
+
+        for pkey, pval in defaults.items():
+            if self.pars[pkey] is None:
+                self.pars[pkey] = pval
+
+        return
+
+    def add_pairs(self, people, ti=0):
+
+        female = people.female
+        active = people.active
+        f_active = female & active
+        m_active = ~female & active
+
+        # Compute number of partners
+        f_partnered_inds, f_partnered_counts = np.unique(self['p1'], return_counts=True)
+        m_partnered_inds, m_partnered_counts = np.unique(self['p2'], return_counts=True)
+        current_partners = np.zeros((len(people)))
+        current_partners[f_partnered_inds] = f_partnered_counts
+        current_partners[m_partnered_inds] = m_partnered_counts
+        partners = ssu.sample(**self.pars['partners'], size=len(people)) + 1
+        underpartnered = current_partners < partners  # Indices of underpartnered people
+        f_eligible = f_active & underpartnered
+        m_eligible = m_active & underpartnered
+
+        # Bin the agents by age
+        bins = self.pars['participation'][0, :]  # Extract age bins
+
+        # Try randomly select females for pairing
+        f_eligible_inds = ssu.true(f_eligible)  # Inds of all eligible females
+        age_bins_f = np.digitize(people.age[f_eligible_inds], bins=bins) - 1  # Age bins of selected females
+        bin_range_f = np.unique(age_bins_f)  # Range of bins
+        f = []  # Initialize the female partners
+        for ab in bin_range_f:  # Loop over age bins
+            these_f_contacts = ssu.binomial_filter(self.pars['participation'][1][ab], f_eligible_inds[
+                age_bins_f == ab])  # Select females according to their participation rate in this layer
+            f += these_f_contacts.tolist()
+
+        # Select males according to their participation rate in this layer
+        m_eligible_inds = ssu.true(m_eligible)
+        age_bins_m = np.digitize(people.age[m_eligible_inds], bins=bins) - 1
+        bin_range_m = np.unique(age_bins_m)  # Range of bins
+        m = []  # Initialize the male partners
+        for ab in bin_range_m:
+            these_m_contacts = ssu.binomial_filter(self.pars['participation'][2][ab], m_eligible_inds[
+                age_bins_m == ab])  # Select males according to their participation rate in this layer
+            m += these_m_contacts.tolist()
+
+        # Create preference matrix between eligible females and males that combines age and geo mixing
+        age_bins_f = np.digitize(people.age[f], bins=bins) - 1  # Age bins of females that are entering new relationships
+        age_bins_m = np.digitize(people.age[m], bins=bins) - 1  # Age bins of active and participating males
+        age_f, age_m = np.meshgrid(age_bins_f, age_bins_m)
+        pair_probs = self.pars['mixing'][age_m, age_f + 1]
+
+        f_to_remove = pair_probs.max(axis=0) == 0  # list of female inds to remove if no male partners are found for her
+        f = [i for i, flag in zip(f, f_to_remove) if ~flag]  # remove the inds who don't get paired on this timestep
+        selected_males = []
+        if len(f):
+            pair_probs = pair_probs[:, np.invert(f_to_remove)]
+            choices = []
+            fems = np.arange(len(f))
+            f_paired_bools = np.full(len(fems), True, dtype=bool)
+            np.random.shuffle(fems)
+            for fem in fems:
+                m_col = pair_probs[:, fem]
+                if m_col.sum() > 0:
+                    m_col_norm = m_col / m_col.sum()
+                    choice = np.random.choice(len(m_col_norm), 1, replace=False, p=m_col_norm)
+                    choices.append(choice)
+                    pair_probs[choice, :] = 0  # Once male partner is assigned, remove from eligible pool
+                else:
+                    f_paired_bools[fem] = False
+            selected_males = np.array(m)[np.array(choices).flatten()]
+            f = np.array(f)[f_paired_bools]
+
+        p1 = np.array(f)
+        p2 = selected_males
+        n_partnerships = len(p1)
+        dur = ssu.sample(**self.pars['dur_pship'], size=n_partnerships)
+        acts = ssu.sample(**self.pars['acts'], size=n_partnerships)
+        age_p1 = people.age[p1]
+        age_p2 = people.age[p2]
+
+        age_debut_p1 = people.debut[p1]
+        age_debut_p2 = people.debut[p2]
+
+        # For each couple, get the average age they are now and the average age of debut
+        avg_age = np.array([age_p1, age_p2]).mean(axis=0)
+        avg_debut = np.array([age_debut_p1, age_debut_p2]).mean(axis=0)
+
+        # Shorten parameter names
+        dr = self.pars['age_act_pars']['debut_ratio']
+        peak = self.pars['age_act_pars']['peak']
+        rr = self.pars['age_act_pars']['retirement_ratio']
+        retire = self.pars['age_act_pars']['retirement']
+
+        # Get indices of people at different stages
+        below_peak_inds = avg_age <= self.pars['age_act_pars']['peak']
+        above_peak_inds = (avg_age > self.pars['age_act_pars']['peak']) & (avg_age < self.pars['age_act_pars']['retirement'])
+        retired_inds = avg_age > self.pars['age_act_pars']['retirement']
+
+        # Set values by linearly scaling the number of acts for each partnership according to
+        # the age of the couple at the commencement of the relationship
+        below_peak_vals = acts[below_peak_inds] * (dr + (1 - dr) / (peak - avg_debut[below_peak_inds]) * (
+                    avg_age[below_peak_inds] - avg_debut[below_peak_inds]))
+        above_peak_vals = acts[above_peak_inds] * (
+                    rr + (1 - rr) / (peak - retire) * (avg_age[above_peak_inds] - retire))
+        retired_vals = 0
+
+        # Set values and return
+        scaled_acts = np.full(len(acts), np.nan, dtype=sss.default_float)
+        scaled_acts[below_peak_inds] = below_peak_vals
+        scaled_acts[above_peak_inds] = above_peak_vals
+        scaled_acts[retired_inds] = retired_vals
+        start = np.array([ti] * n_partnerships, dtype=sss.default_float)
+        beta = np.ones(n_partnerships)
+
+        new_contacts = dict(
+            p1=p1,
+            p2=p2,
+            dur=dur,
+            acts=acts,
+            start=start,
+            beta=beta
+        )
+        self.append(new_contacts)
+        return
+
+    def update(self, people, ti=None, dt=None):
+        if ti is None: ti = people.ti
+        if dt is None: dt = people.dt
+        # First remove any relationships due to end
+        self['dur'] = self['dur'] - dt
+        active = self['dur'] > 0
+        for key in self.meta.keys():
+            self[key] = self[key][active]
+
+        # Then add new relationships
+        self.add_pairs(people, ti=ti)
+        return
+
+
+class maternal(ssb.Network):
+    def __init__(self, key_dict=None, transmission='vertical'):
+        """
+        Initialized empty and filled with pregnancies throughout the simulation
+        """
+        key_dict = sc.mergedicts({'dur': sss.default_float}, key_dict)
+        super().__init__(key_dict=key_dict, transmission=transmission)
+        return
+
+    def update(self, people, dt=None):
+        if dt is None: dt = people.dt
+        # Set beta to 0 for women who complete post-partum period
+        # Keep connections for now, might want to consider removing
+        self['dur'] = self['dur'] - dt
+        inactive = self['dur'] <= 0
+        self['beta'][inactive] = 0
+
+    def initialize(self, people):
+        """ No pairs added upon initialization """
+        pass
+
+    def add_pairs(self, mother_inds, unborn_inds, dur):
+        """
+        Add connections between pregnant women and their as-yet-unborn babies
+        """
+        beta = np.ones_like(mother_inds)
+        self['p1'] = np.concatenate([self['p1'], mother_inds])
+        self['p2'] = np.concatenate([self['p2'], unborn_inds])
+        self['beta'] = np.concatenate([self['beta'], beta])
+        self['dur'] = np.concatenate([self['dur'], dur])
+        return
diff --git a/stisim/parameters.py b/stisim/parameters.py
new file mode 100644
index 00000000..fc97b1c6
--- /dev/null
+++ b/stisim/parameters.py
@@ -0,0 +1,90 @@
+"""
+Set parameters
+"""
+
+import numpy as np
+import sciris as sc
+from .settings import options as sso  # For setting global options
+# from . import misc as ssm
+# from .data import loaders as ssdata
+
+__all__ = ['make_pars']
+
+
+def make_pars(**kwargs):
+    """
+    Create the parameters for the simulation. Typically, this function is used
+    internally rather than called by the user; e.g. typical use would be to do
+    sim = ss.Sim() and then inspect sim.pars, rather than calling this function
+    directly.
+
+    Args:
+        kwargs        (dict): any additional kwargs are interpreted as parameter names
+    Returns:
+        pars (dict): the parameters of the simulation
+    """
+    pars = sc.objdict()
+
+    # Population parameters
+    pars['n_agents']        = 10e3          # Number of agents
+    pars['total_pop']       = 10e3          # If defined, used for calculating the scale factor
+    pars['pop_scale']       = None          # How much to scale the population
+    pars['location']        = None          # What demographics to use - NOT CURRENTLY FUNCTIONAL
+    pars['birth_rates']     = None          # Birth rates, loaded below
+    pars['death_rates']     = None          # Death rates, loaded below
+    pars['rel_birth']       = 1.0           # Birth rate scale factor
+    pars['rel_death']       = 1.0           # Death rate scale factor
+
+    # Simulation parameters
+    pars['start']           = 1995.         # Start of the simulation
+    pars['end']             = None          # End of the simulation
+    pars['n_years']         = 35            # Number of years to run, if end isn't specified. Note that this includes burn-in
+    pars['burnin']          = 25            # Number of years of burnin. NB, this is doesn't affect the start and end dates of the simulation, but it is possible remove these years from plots
+    pars['dt']              = 1.0           # Timestep (in years)
+    pars['dt_demog']        = 1.0           # Timestep for demographic updates (in years)
+    pars['rand_seed']       = 1             # Random seed, if None, don't reset
+    pars['verbose']         = sso.verbose   # Whether or not to display information during the run -- options are 0 (silent), 0.1 (some; default), 1 (default), 2 (everything)
+    pars['use_migration']   = True          # Whether to estimate migration rates to correct the total population size
+
+    # Events and interventions
+    pars['connectors']      = sc.autolist()
+    pars['interventions']   = sc.autolist() # The interventions present in this simulation; populated by the user
+    pars['analyzers']       = sc.autolist() # The functions present in this simulation; populated by the user
+    pars['timelimit']       = None          # Time limit for the simulation (seconds)
+    pars['stopping_func']   = None          # A function to call to stop the sim partway through
+
+    # Network parameters, generally initialized after the population has been constructed
+    pars['networks']        = sc.autolist()  # Network types and parameters
+    pars['debut']           = dict(f=dict(dist='normal', par1=15.0, par2=2.0),
+                                   m=dict(dist='normal', par1=17.5, par2=2.0))
+
+    # Update with any supplied parameter values and generate things that need to be generated
+    pars.update(kwargs)
+
+    return pars
+
+
+# def get_births_deaths(location, verbose=1, by_sex=True, overall=False, die=True):
+#     """
+#     Get mortality and fertility data by location if provided, or use default
+#
+#     Args:
+#         location (str):  location
+#         verbose (bool):  whether to print progress
+#         by_sex   (bool): whether to get sex-specific death rates (default true)
+#         overall  (bool): whether to get overall values ie not disaggregated by sex (default false)
+#
+#     Returns:
+#         lx (dict): dictionary keyed by sex, storing arrays of lx - the number of people who survive to age x
+#         birth_rates (arr): array of crude birth rates by year
+#     """
+#
+#     if verbose:
+#         print(f'Loading location-specific demographic data for "{location}"')
+#     try:
+#         death_rates = ssdata.get_death_rates(location=location, by_sex=by_sex, overall=overall)
+#         birth_rates = ssdata.get_birth_rates(location=location)
+#         return birth_rates, death_rates
+#     except ValueError as E:
+#         warnmsg = f'Could not load demographic data for requested location "{location}" ({str(E)})'
+#         ssm.warn(warnmsg, die=die)
diff --git a/stisim/people.py b/stisim/people.py
index 54b76d81..d5285923 100644
--- a/stisim/people.py
+++ b/stisim/people.py
@@ -1,103 +1,164 @@
-import sciris as sc
-import numpy as np
-from .results import Result
-
-
-class State():
-    # Simplified version of hpvsim.State()
-    def __init__(self, name, dtype, fill_value=0):
-        self.name = name
-        self.dtype = dtype
-        self.fill_value = fill_value
-
-    def new(self, n):
-        return np.full(n, dtype=self.dtype, fill_value=self.fill_value)
-
-
-class People(sc.prettyobj):
-    # TODO - cater to use cases of
-    #   - Initial contact networks are independent of modules and we want to pre-generate agents and their contact layers
-    #   - Initial contact networks depend on modules, and we need to add the modules before adding subsequent contact layers
-
-    # Define states that every People instance has, regardless of which modules are enabled
-    base_states = [
-        State('uid', int), # TODO: will we support removing agents? It could make indexing much more complicated...
-        State('age', float),
-        State('female', bool, False),
-        State('dead', bool, False),
-        State('ti_dead', float, np.nan), # Time index for death - defaults to natural causes but gets overwritten if they die of something first
-    ]
-
-    def __len__(self):
-        return len(self.uid)
-
-    def __getitem__(self, key):
-        return self.__getattribute__(key)
-
-    def __setitem__(self, key, value):
-        return self.__setattr__(key, value)
-
-    @property
-    def male(self):
-        return ~self.female
-
-    @property
-    def n(self):
-        return len(self)
-
-    @property
-    def indices(self):
-        return self.uid
-
-    def __init__(self, n):
-        # TODO - where is the right place to change how the initial population is generated?
-        for state in self.base_states:
-            self.__setattr__(state.name, state.new(n))
-        self.uid[:] = np.arange(n)
-        self.age[:] = np.random.random(n)
-        self.female[:] = np.random.randint(0,2,n)
-        self.contacts = sc.odict()  # Dict containing Layer instances
-        self._modules = []
-
-    def add_module(self, module):
-        # Initialize all of the states associated with a module
-        # This is implemented as People.add_module rather than
-        # Module.add_to_people(people) or similar because its primary
-        # role is to modify the
-        if hasattr(self,module.name):
-            raise Exception(f'Module {module.name} already added')
-        self.__setattr__(module.name, sc.objdict())
-        for state in module.states:
-            self[module.name][state.name] = state.new(self.n)
-
-    # These methods provide a single function to call that handles
-    # autonomous updates of the people, so it simplifies Sim.run()
-    # However it's a little hacky that People and Module both define
-    # these methods. Would it make sense for the base variables like
-    # UID, age etc. to be in a Module? Or are they so heavily baked into
-    # People that we don't want to separate them out?
-    def update_states_pre(self, sim):
-        self.dead[self.ti_dead <= sim.ti] = True
-
-        for module in sim.modules:
-            module.update_states_pre(sim)
-
-    def initialize(self, sim):
-        sim.results['n_people'] = Result(None, 'n_people', sim.npts, dtype=int)
-        sim.results['new_deaths'] = Result(None, 'new_deaths', sim.npts, dtype=int)
-
-    def update_results(self, sim):
-        sim.results['new_deaths'] = np.count_nonzero(sim.people.ti_dead == sim.ti)
-
-        for module in sim.modules:
-            module.update_results(sim)
-
-    def finalize_results(self, sim):
-        pass
-
-
-def make_people(pars):
-    # A generic function to handle making people via parameters during sim.initialize()
-    # Users could also call this function to pre-generate people using the same routine beforehand
-    # Or otherwise, they can write a custom equivalent of make_people and pass it to the sim directly
-    raise NotImplementedError
\ No newline at end of file
+"""
+Defines the People class and functions associated with making people
+"""
+
+# %% Imports
+import sciris as sc
+import numpy as np
+from . import base as ssb
+from . import misc as ssm
+from . import utils as ssu
+from . import settings as sss
+# from .data import loaders as ssdata
+
+
+__all__ = ['People', 'make_popdict']
+
+
+# %% Main people class
+
+
+class People(ssb.BasePeople):
+    """
+    A class to perform all the operations on the people
+    This class is usually created automatically by the sim. The only required input
+    argument is the population size, but typically the full parameters dictionary
+    will get passed instead since it will be needed before the People object is
+    initialized.
+
+    Note that this class handles the mechanics of updating the actual people, while
+    ``ss.BasePeople`` takes care of housekeeping (saving, loading, exporting, etc.).
+    Please see the BasePeople class for additional methods.
+
+    Args:
+        pars (dict): the sim parameters, e.g. sim.pars -- alternatively, if a number, interpreted as n_agents
+        strict (bool): whether to only create keys that are already in self.meta.person; otherwise, let any key be set
+        pop_trend (dataframe): a dataframe of years and population sizes, if available
+        kwargs (dict): the actual data, e.g. from a popdict, being specified
+
+    **Examples**::
+        ppl = ss.People(2000)
+    """
+
+    # %% Basic methods
+
+    def __init__(self, n, states=None, strict=True, **kwargs):
+        """
+        Initialize
+        """
+        super().__init__(n, states=states)
+        if strict: self.lock()  # If strict is true, stop further keys from being set (does not affect attributes)
+        self.kwargs = kwargs
+        return
+
+    def initialize(self, popdict=None):
+        """ Perform initializations """
+        super().initialize(popdict=popdict)  # Initialize states
+        return
+
+    def add_module(self, module, force=False):
+        # Initialize all the states associated with a module
+        # This is implemented as People.add_module rather than
+        # Module.add_to_people(people) or similar because its primary
+        # role is to modify the People object
+        if hasattr(self, module.name) and not force:
+            raise Exception(f'Module {module.name} already added')
+        self.__setattr__(module.name, sc.objdict())
+
+        for state_name, state in module.states.items():
+            combined_name = module.name + '.' + state_name
+            self._data[combined_name] = state.new(self._n)
+            self._map_arrays(keys=combined_name)
+            self.states[combined_name] = state
+
+        return
+
+    def scale_flows(self, inds):
+        """
+        Return the scaled versions of the flows -- replacement for len(inds)
+        followed by scale factor multiplication
+        """
+        return self.scale[inds].sum()
+
+    def update_states(self, sim):
+        """ Perform all state updates at the current timestep """
+
+        self.age[~self.dead] += sim.dt
+        self.dead[self.ti_dead <= sim.ti] = True
+
+        for module in sim.modules.values():
+            module.update_states(sim)
+
+        # Perform network updates
+        for lkey, layer in self.networks.items():
+            layer.update(self)
+
+        return
+
+
+# %% Helper functions to create popdicts
+
+def set_static(new_n, existing_n=0, pars=None, f_ratio=0.5):
+    """
+    Set static population characteristics that do not change over time.
+    This function can be used when adding new births, in which case the existing popsize can be given as `existing_n`.
+
+    Arguments:
+        new_n (int): Number of new individuals to add to the population.
+        existing_n (int, optional): Number of existing individuals in the population. Default is 0.
+        pars (dict, optional): Dictionary of parameters. Default is None.
+        f_ratio (float, optional): Female ratio in the population. Default is 0.5.
+
+    Returns:
+        uid (ndarray, int): unique identifiers for the individuals.
+        female (ndarray, bool): array indicating whether an individual is female or not.
+        debut (ndarray, bool): array indicating the debut value for each individual.
+    """
+    uid = np.arange(existing_n, existing_n+new_n, dtype=sss.default_int)
+    female = np.random.choice([True, False], size=new_n, p=f_ratio)
+    n_females = len(ssu.true(female))
+    debut = np.full(new_n, np.nan, dtype=sss.default_float)
+    debut[female] = ssu.sample(**pars['debut']['m'], size=n_females)
+    debut[~female] = ssu.sample(**pars['debut']['f'], size=new_n-n_females)
+    return uid, female, debut
+
+
+def make_popdict(n=None, location=None, year=None, verbose=None, f_ratio=0.5, dt_round_age=True, dt=None):
+    """ Create a location-specific population dictionary """
+
+    # Initialize total pop size
+    total_pop = None
+
+    # Load age data for this location if available
+    if verbose:
+        print(f'Loading location-specific data for "{location}"')
+    try:
+        age_data = ssdata.get_age_distribution(location, year=year)
+        total_pop = sum(age_data[:, 2])  # Return the total population
+    except ValueError as E:
+        warnmsg = f'Could not load age data for requested location "{location}" ({str(E)})'
+        ssm.warn(warnmsg, die=True)
+
+    uid, female, debut = set_static(n, f_ratio=f_ratio)
+
+    # Set ages, rounding to nearest timestep if requested
+    age_data_min = age_data[:, 0]
+    age_data_max = age_data[:, 1]
+    age_data_range = age_data_max - age_data_min
+    age_data_prob = age_data[:, 2]
+    age_data_prob /= age_data_prob.sum()  # Ensure it sums to 1
+    age_bins = ssu.n_multinomial(age_data_prob, n)  # Choose age bins
+    if dt_round_age:
+        ages = age_data_min[age_bins] + np.random.randint(
+            age_data_range[age_bins] / dt) * dt  # Uniformly distribute within this age bin
+    else:
+        ages = age_data_min[age_bins] + age_data_range[age_bins] * np.random.random(n)  # Uniform  within this age bin
+
+    # Store output
+    popdict = dict(
+        uid=uid,
+        age=ages,
+        female=female
+    )
+
+    return total_pop, popdict
diff --git a/stisim/population.py b/stisim/population.py
deleted file mode 100644
index 944560ec..00000000
--- a/stisim/population.py
+++ /dev/null
@@ -1,80 +0,0 @@
-import sciris as sc
-import numpy as np
-
-
-class Layer(sc.odict):
-    def __init__(self):
-        super().__init__()
-        self.p1 = np.array([], dtype=int)
-        self.p2 = np.array([], dtype=int)
-        self.beta = np.array([], dtype=float)
-
-    def __repr__(self):
-        return f'<{self.__class__.__name__}, {len(self.members)} members, {len(self.p1)} contacts>'
-
-    @property
-    def members(self):
-        return set(self.p1).union(set(self.p2))
-
-    def update(self):
-        pass
-
-# TODO - modify network to have long-term and casual partnerships with different transmission risk per contact
-class RandomDynamicSexualLayer(Layer):
-    # Randomly pair males and females with variable relationship durations
-    def __init__(self, people, mean_dur=5):
-        super().__init__()
-        self.mean_dur = mean_dur
-        self.dur = np.array([], dtype=float)
-        self.add_partnerships(people)
-
-    def add_partnerships(self, people):
-        # Find unpartnered males and females - could in principle check other contact layers too
-        # by having the People object passed in here
-
-        available_m = np.setdiff1d(people.indices[people.male], self.members)
-        available_f = np.setdiff1d(people.indices[~people.male], self.members)
-
-        if len(available_m) <= len(available_f):
-            p1 = available_m
-            p2 = np.random.choice(available_f, len(p1), replace=False)
-        else:
-            p2 = available_f
-            p1 = np.random.choice(available_m, len(p2), replace=False)
-
-        beta = np.ones_like(p1)
-        dur = np.random.randn(len(p1))*self.mean_dur
-        self.p1 = np.concatenate([self.p1, p1])
-        self.p2 = np.concatenate([self.p2, p2])
-        self.beta = np.concatenate([self.beta, beta])
-        self.dur = np.concatenate([self.dur, dur])
-
-    def update(self, people):
-        # First remove any relationships due to end
-        self.dur = self.dur - people.dt
-        active = self.dur > 0
-        self.p1 = self.p1[active]
-        self.p2 = self.p2[active]
-        self.beta = self.beta[active]
-
-        # Then add new relationships for unpartnered people
-        self.add_partnerships(people)
-
-class StaticLayer(Layer):
-    # Randomly make some partnerships that don't change over time
-    def __init__(self, people, sex='mf'):
-        super().__init__()
-
-        if sex[0]=='m':
-            p1 = people.indices[people.male]
-        else:
-            p1 = people.indices[~people.male]
-
-        if sex[1]=='f':
-            p2 = people.indices[~people.male]
-        else:
-            p2 = people.indices[people.male]
-
-        self.p1 = p1
-        self.p2 = np.random.choice(p2, len(p1), replace=True)
-        self.beta = np.ones_like(p1)
\ No newline at end of file
diff --git a/stisim/results.py b/stisim/results.py
index 671a39d6..779bbaec 100644
--- a/stisim/results.py
+++ b/stisim/results.py
@@ -1,23 +1,33 @@
-import numpy as np
-import sciris as sc
-
-class Result(np.lib.mixins.NDArrayOperatorsMixin):
-    def __init__(self, module, name, n, dtype):
-        self.module = module if module else None
-        self.name = name
-        self.values = np.zeros(n, dtype=dtype)
-
-    def __getitem__(self, *args, **kwargs):  return self.values.__getitem__(*args, **kwargs)
-    def __setitem__(self, *args, **kwargs): return self.values.__setitem__(*args, **kwargs)
-    def __len__(self): return len(self.values)
-    def __repr__(self): return f'Result({self.module},{self.name}): {self.values.__repr__()}'
-
-    # These methods allow automatic use of functions like np.sum, np.exp, etc.
-    # with higher performance in some cases
-    @property
-    def __array_interface__(self): return self.values.__array_interface__
-    def __array__(self): return self.values
-    def __array_ufunc__(self, *args, **kwargs):
-        args = [(x if x is not self else self.values) for x in args]
-        kwargs = {k:v if v is not self else self.values for k,v in kwargs.items()}
-        return self.values.__array_ufunc__(*args, **kwargs)
\ No newline at end of file
+"""
+Results and associated structures. Currently only holds results, could move to base.py
+unless other things get added (e.g. Resultsets, MultiResults, other...)
+"""
+
+import numpy as np
+
+
+class Result(np.lib.mixins.NDArrayOperatorsMixin):
+    def __init__(self, name, module, shape, dtype):
+        self.name = name
+        self.module = module if module else None
+        self.values = np.zeros(shape, dtype=dtype)
+
+    def __getitem__(self, *args, **kwargs):  return self.values.__getitem__(*args, **kwargs)
+
+    def __setitem__(self, *args, **kwargs): return self.values.__setitem__(*args, **kwargs)
+
+    def __len__(self): return len(self.values)
+
+    def __repr__(self): return f'Result({self.module},{self.name}): {self.values.__repr__()}'
+
+    # These methods allow automatic use of functions like np.sum, np.exp, etc.
+    # with higher performance in some cases
+    @property
+    def __array_interface__(self): return self.values.__array_interface__
+
+    def __array__(self): return self.values
+
+    def __array_ufunc__(self, *args, **kwargs):
+        args = [(x if x is not self else self.values) for x in args]
+        kwargs = {k: v if v is not self else self.values for k, v in kwargs.items()}
+        return self.values.__array_ufunc__(*args, **kwargs)
diff --git a/stisim/settings.py b/stisim/settings.py
new file mode 100644
index 00000000..4e2142dc
--- /dev/null
+++ b/stisim/settings.py
@@ -0,0 +1,589 @@
+"""
+Define options for STIsim.
+All options should be set using set() or directly, e.g.::
+
+    ss.options(font_size=18)
+
+To reset default options, use::
+
+    ss.options('default')
+
+Note: "options" is used to refer to the choices available (e.g., DPI), while "settings"
+is used to refer to the choices made (e.g., DPI=150).
+"""
+
+import os
+import pylab as pl
+import sciris as sc
+import numpy as np
+import matplotlib.font_manager as fm
+
+# Specify all externally visible classes this file defines
+__all__ = ['options']
+
+# Define simple plotting options -- similar to Matplotlib default
+rc_simple = {
+    'axes.axisbelow':    True, # So grids show up behind
+    'axes.spines.right': False,
+    'axes.spines.top':   False,
+    'figure.facecolor':  'white',
+    'font.family':       'sans-serif', # Replaced with Mulish in load_fonts() if import succeeds
+    'legend.frameon':    False,
+}
+
+# Define default plotting options -- based on Seaborn
+rc_stisim_extra = {
+    'axes.facecolor': '#f2f2ff',
+    'axes.grid':      True,
+    'grid.color':     'white',
+    'grid.linewidth': 1,
+}
+rc_stisim = sc.mergedicts(rc_simple, rc_stisim_extra)
+
+# %% Define the options class
+
+class Options(sc.objdict):
+    """
+    Set options for STIsim.
+
+    Use ``ss.options.set('defaults')`` to reset all values to default, or ``ss.options.set(dpi='default')``
+    to reset one parameter to default. See ``ss.options.help(detailed=True)`` for
+    more information.
+
+    Options can also be saved and loaded using ``ss.options.save()`` and ``ss.options.load()``.
+    See ``ss.options.context()`` and ``ss.options.with_style()`` to set options
+    temporarily.
+
+    Common options are (see also ``ss.options.help(detailed=True)``):
+
+        - verbose:        default verbosity for simulations to use
+        - style:          the plotting style to use
+        - dpi:            the overall DPI (i.e. size) of the figures
+        - font:           the font family/face used for the plots
+        - fontsize:       the font size used for the plots
+        - interactive:    convenience method to set show, close, and backend
+        - jupyter:        defaults for Jupyter (change backend and figure return)
+        - show:           whether to show figures
+        - close:          whether to close the figures
+        - backend:        which Matplotlib backend to use
+        - warnings:       how to handle warnings (e.g. print, raise as errors, ignore)
+
+    **Examples**::
+
+        ss.options(dpi=150) # Larger size
+        ss.options(style='simple', font='Rosario') # Change to the "simple" STIsim style with a custom font
+        ss.options.set(fontsize=18, show=False, backend='agg', precision=64) # Multiple changes
+        ss.options(interactive=False) # Turn off interactive plots
+        ss.options(jupyter=True) # Defaults for Jupyter
+        ss.options('defaults') # Reset to default options
+
+    """
+
+    def __init__(self):
+        super().__init__()
+        optdesc, options = self.get_orig_options()  # Get the options
+        self.update(options)  # Update this object with them
+        self.setattribute('optdesc', optdesc)  # Set the description as an attribute, not a dict entry
+        self.setattribute('orig_options', sc.dcp(options))  # Copy the default options
+        return
+
+    def __call__(self, *args, **kwargs):
+        """Allow ``ss.options(dpi=150)`` instead of ``ss.options.set(dpi=150)`` """
+        return self.set(*args, **kwargs)
+
+    def to_dict(self):
+        """ Pull out only the settings from the options object """
+        return {k: v for k, v in self.items()}
+
+    def __repr__(self):
+        """ Brief representation """
+        output = sc.objectid(self)
+        output += 'STIsim options (see also ss.options.disp()):\n'
+        output += sc.pp(self.to_dict(), output=True)
+        return output
+
+    def __enter__(self):
+        """ Allow to be used in a with block """
+        return self
+
+    def __exit__(self, *args, **kwargs):
+        """ Allow to be used in a with block """
+        try:
+            reset = {}
+            for k, v in self.on_entry.items():
+                if self[k] != v:  # Only reset settings that have changed
+                    reset[k] = v
+            self.set(**reset)
+            self.delattribute('on_entry')
+        except AttributeError as E:
+            errormsg = 'Please use ss.options.context() if using a with block'
+            raise AttributeError(errormsg) from E
+        return
+
+    def disp(self):
+        """ Detailed representation """
+        output = 'STIsim options (see also ss.options.help()):\n'
+        keylen = 10  # Maximum key length  -- "interactive"
+        for k, v in self.items():
+            keystr = sc.colorize(f'  {k:>{keylen}s}: ', fg='cyan', output=True)
+            reprstr = sc.pp(v, output=True)
+            reprstr = sc.indent(n=keylen + 4, text=reprstr, width=None)
+            output += f'{keystr}{reprstr}'
+        print(output)
+        return
+
+    @staticmethod
+    def get_orig_options():
+        """
+        Set the default options for STIsim -- not to be called by the user, use
+        ``ss.options.set('defaults')`` instead.
+        """
+
+        # Options acts like a class, but is actually an objdict for simplicity
+        optdesc = sc.objdict()  # Help for the 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('STISIM_VERBOSE', 0.1))
+
+        optdesc.style = 'Set the default plotting style -- options are "STIsim" and "simple" plus those in ' \
+                        'pl.style.available; see also options.rc '
+        options.style = os.getenv('STISIM_STYLE', 'stisim')
+
+        optdesc.dpi = 'Set the default DPI -- the larger this is, the larger the figures will be'
+        options.dpi = int(os.getenv('STISIM_DPI', pl.rcParams['figure.dpi']))
+
+        optdesc.font = 'Set the default font family (e.g., sans-serif or Arial)'
+        options.font = os.getenv('STISIM_FONT', pl.rcParams['font.family'])
+
+        optdesc.fontsize = 'Set the default font size'
+        options.fontsize = int(os.getenv('STISIM_FONT_SIZE', pl.rcParams['font.size']))
+
+        optdesc.interactive = 'Convenience method to set figure backend, showing, and closing behavior'
+        options.interactive = os.getenv('STISIM_INTERACTIVE', True)
+
+        optdesc.jupyter = 'Convenience method to set common settings for Jupyter notebooks: set to "retina" or ' \
+                          '"widget" (default) to set backend '
+        options.jupyter = os.getenv('STISIM_JUPYTER', False)
+
+        optdesc.show = 'Set whether or not to show figures (i.e. call pl.show() automatically)'
+        options.show = int(os.getenv('STISIM_SHOW', True))
+
+        optdesc.close = 'Set whether or not to close figures (i.e. call pl.close() automatically)'
+        options.close = int(os.getenv('STISIM_CLOSE', False))
+
+        optdesc.returnfig = 'Set whether or not to return figures from plotting functions'
+        options.returnfig = int(os.getenv('STISIM_RETURNFIG', True))
+
+        optdesc.backend = 'Set the Matplotlib backend (use "agg" for non-interactive)'
+        options.backend = os.getenv('STISIM_BACKEND', pl.get_backend())
+
+        optdesc.rc = 'Matplotlib rc (run control) style parameters used during plotting -- usually set automatically ' \
+                     'by "style" option'
+        options.rc = sc.dcp(rc_stisim)
+
+        optdesc.warnings = 'How warnings are handled: options are "warn" (default), "print", and "error"'
+        options.warnings = str(os.getenv('STISIM_WARNINGS', 'warn'))
+
+        optdesc.sep = 'Set thousands seperator for text output'
+        options.sep = str(os.getenv('STISIM_SEP', ','))
+
+        optdesc.precision = 'Set arithmetic precision -- 32-bit by default for efficiency'
+        options.precision = int(os.getenv('STISIM_PRECISION', 32))
+
+        return optdesc, options
+
+    def set(self, key=None, value=None, use=False, **kwargs):
+        """
+        Actually change the style. See ``ss.options.help()`` for more information.
+
+        Args:
+            key    (str):    the parameter to modify, or 'defaults' to reset everything to default
+            value  (varies): the value to specify; use None or 'default' to reset to default
+            use (bool): whether to use the chosen style
+            kwargs (dict):   if supplied, set multiple key-value pairs
+
+        **Example**::
+            ss.options.set(dpi=50) # Equivalent to ss.options(dpi=50)
+        """
+
+        # Reset to defaults
+        if key in ['default', 'defaults']:
+            kwargs = self.orig_options  # Reset everything to default
+
+        # Handle other keys
+        elif key is not None:
+            kwargs = sc.mergedicts(kwargs, {key: value})
+
+        # Handle Jupyter
+        if 'jupyter' in kwargs.keys() and kwargs['jupyter']:
+            jupyter = kwargs['jupyter']
+            if jupyter:
+                jupyter = 'retina'  # Default option for True
+            try:
+                if not os.environ.get(
+                        'SPHINX_BUILD'):  # Custom check implemented in conf.py to skip this if we're inside Sphinx
+                    if jupyter == 'retina':  # This makes plots much nicer, but isn't available on all systems
+                        import matplotlib_inline
+                        matplotlib_inline.backend_inline.set_matplotlib_formats('retina')
+                    elif jupyter in ['widget', 'interactive']:  # Or use interactive
+                        from IPython import get_ipython
+                        magic = get_ipython().magic
+                        magic('%matplotlib widget')
+            except:
+                pass
+
+        # Handle interactivity
+        if 'interactive' in kwargs.keys():
+            interactive = kwargs['interactive']
+            if interactive in [None, 'default']:
+                interactive = self.orig_options['interactive']
+            if interactive:
+                kwargs['show'] = True
+                kwargs['close'] = False
+                kwargs['backend'] = self.orig_options['backend']
+            else:
+                kwargs['show'] = False
+                kwargs['backend'] = 'agg'
+
+        # Reset options
+        for key, value in kwargs.items():
+
+            # Handle deprecations
+            rename = {'font_size': 'fontsize', 'font_family': 'font'}
+            if key in rename.keys():
+                from . import misc as ssm  # Here to avoid circular import
+                oldkey = key
+                key = rename[oldkey]
+                warnmsg = f'Key "{oldkey}" is deprecated, please use "{key}" instead'
+                ssm.warn(warnmsg, FutureWarning)
+
+            if key not in self:
+                keylist = self.orig_options.keys()
+                keys = '\n'.join(keylist)
+                errormsg = f'Option "{key}" not recognized; options are "defaults" or:\n{keys}\n\nSee help(ss.options.set) for more information.'
+                raise sc.KeyNotFoundError(errormsg)
+            else:
+                if value in [None, 'default']:
+                    value = self.orig_options[key]
+                self[key] = value
+
+                matplotlib_keys = ['fontsize', 'font', 'dpi', 'backend']
+                if key in matplotlib_keys:
+                    self.set_matplotlib_global(key, value)
+
+        if use:
+            self.use_style()
+
+        return
+
+    def set_matplotlib_global(self, key, value):
+        """ Set a global option for Matplotlib -- not for users """
+        if value:  # Don't try to reset any of these to a None value
+            if key == 'fontsize':
+                pl.rcParams['font.size'] = value
+            elif key == 'font':
+                pl.rcParams['font.family'] = value
+            elif key == 'dpi':
+                pl.rcParams['figure.dpi'] = value
+            elif key == 'backend':
+                pl.switch_backend(value)
+            else:
+                raise KeyError(f'Key {key} not found')
+        return
+
+    def context(self, **kwargs):
+        """
+        Alias to set() for non-plotting options, for use in a "with" block.
+
+        Note: for plotting options, use ``ss.options.with_style()``, which is linked
+        to Matplotlib's context manager. If you set plotting options with this,
+        they won't have any effect.
+
+        **Examples**::
+
+            # Silence all output
+            with ss.options.context(verbose=0):
+                ss.Sim().run()
+
+            # Convert warnings to errors
+            with ss.options.context(warnings='error'):
+                ss.Sim(location='not a location').initialize()
+
+            # Use with_style(), not context(), for plotting options
+            with ss.options.with_style(dpi=50):
+                ss.Sim().run().plot()
+
+        """
+
+        # Store current settings
+        on_entry = {k: self[k] for k in kwargs.keys()}
+        self.setattribute('on_entry', on_entry)
+
+        # Make changes
+        self.set(**kwargs)
+        return self
+
+    def get_default(self, key):
+        """ Helper function to get the original default options """
+        return self.orig_options[key]
+
+    def changed(self, key):
+        """ Check if current setting has been changed from default """
+        if key in self.orig_options:
+            return self[key] != self.orig_options[key]
+        else:
+            return None
+
+    def help(self, detailed=False, output=False):
+        """
+        Print information about options.
+
+        Args:
+            detailed (bool): whether to print out full help
+            output (bool): whether to return a list of the options
+
+        **Example**::
+
+            ss.options.help(detailed=True)
+        """
+
+        # If not detailed, just print the docstring for ss.options
+        if not detailed:
+            print(self.__doc__)
+            return
+
+        n = 15  # Size of indent
+        optdict = sc.objdict()
+        for key in self.orig_options.keys():
+            entry = sc.objdict()
+            entry.key = key
+            entry.current = sc.indent(n=n, width=None, text=sc.pp(self[key], output=True)).rstrip()
+            entry.default = sc.indent(n=n, width=None, text=sc.pp(self.orig_options[key], output=True)).rstrip()
+            if not key.startswith('rc'):
+                entry.variable = f'STISIM_{key.upper()}'  # NB, hard-coded above!
+            else:
+                entry.variable = 'No environment variable'
+            entry.desc = sc.indent(n=n, text=self.optdesc[key])
+            optdict[key] = entry
+
+        # Convert to a dataframe for nice printing
+        print('stisim global options ("Environment" = name of corresponding environment variable):')
+        for k, key, entry in optdict.enumitems():
+            sc.heading(f'{k}. {key}', spaces=0, spacesafter=0)
+            changestr = '' if entry.current == entry.default else ' (modified)'
+            print(f'          Key: {key}')
+            print(f'      Current: {entry.current}{changestr}')
+            print(f'      Default: {entry.default}')
+            print(f'  Environment: {entry.variable}')
+            print(f'  Description: {entry.desc}')
+
+        sc.heading('Methods:', spacesafter=0)
+        print('''
+    ss.options(key=value) -- set key to value
+    ss.options[key] -- get or set key
+    ss.options.set() -- set option(s)
+    ss.options.get_default() -- get default setting(s)
+    ss.options.load() -- load settings from file
+    ss.options.save() -- save settings to file
+    ss.options.to_dict() -- convert to dictionary
+    ss.options.style() -- create style context for plotting
+''')
+
+        if output:
+            return optdict
+        else:
+            return
+
+    def load(self, filename, verbose=True, **kwargs):
+        """
+        Load current settings from a JSON file.
+
+        Args:
+            filename (str): file to load
+            kwargs (dict): passed to ``sc.loadjson()``
+        """
+        json = sc.loadjson(filename=filename, **kwargs)
+        current = self.to_dict()
+        new = {k: v for k, v in json.items() if v != current[k]}  # Don't reset keys that haven't changed
+        self.set(**new)
+        if verbose: print(f'Settings loaded from {filename}')
+        return
+
+    def save(self, filename, verbose=True, **kwargs):
+        """
+        Save current settings as a JSON file.
+
+        Args:
+            filename (str): file to save to
+            kwargs (dict): passed to ``sc.savejson()``
+        """
+        json = self.to_dict()
+        output = sc.savejson(filename=filename, obj=json, **kwargs)
+        if verbose: print(f'Settings saved to {filename}')
+        return output
+
+    def _handle_style(self, style=None, reset=False, copy=True):
+        """ Helper function to handle logic for different styles """
+        rc = self.rc  # By default, use current
+        if isinstance(style, dict):  # If an rc-like object is supplied directly
+            rc = sc.dcp(style)
+        elif style is not None:  # Usual use case
+            stylestr = str(style).lower()
+            if stylestr in ['default', 'stisim', 'house']:
+                rc = sc.dcp(rc_stisim)
+            elif stylestr in ['simple', 'stisim_simple', 'plain', 'clean']:
+                rc = sc.dcp(rc_simple)
+            elif style in pl.style.library:
+                rc = sc.dcp(pl.style.library[style])
+            else:
+                errormsg = f'Style "{style}"; not found; options are "stisim" (default), "simple", plus:\n{sc.newlinejoin(pl.style.available)} '
+                raise ValueError(errormsg)
+        if reset:
+            self.rc = rc
+        if copy:
+            rc = sc.dcp(rc)
+        return rc
+
+    def with_style(self, style_args=None, use=False, **kwargs):
+        """
+        Combine all Matplotlib style information, and either apply it directly
+        or create a style context.
+
+        To set globally, use ``ss.options.use_style()``. Otherwise, use ``ss.options.with_style()``
+        as part of a ``with`` block to set the style just for that block (using
+        this function outsde of a with block and with ``use=False`` has no effect, so
+        don't do that!). To set non-style options (e.g. warnings, verbosity) as
+        a context, see ``ss.options.context()``.
+
+        Args:
+            style_args (dict): a dictionary of style arguments
+            use (bool): whether to set as the global style; else, treat as context for use with "with" (default)
+            kwargs (dict): additional style arguments
+
+        Valid style arguments are:
+
+            - ``dpi``:       the figure DPI
+            - ``font``:      font (typeface)
+            - ``fontsize``:  font size
+            - ``grid``:      whether to plot gridlines
+            - ``facecolor``: color of the axes behind the plot
+            - any of the entries in ``pl.rParams``
+
+        **Examples**::
+
+            with ss.options.with_style(dpi=300): # Use default options, but higher DPI
+                pl.plot([1,3,6])
+        """
+        # Handle inputs
+        rc = sc.dcp(self.rc)  # Make a local copy of the currently used settings
+        kwargs = sc.mergedicts(style_args, kwargs)
+
+        # Handle style, overwiting existing
+        style = kwargs.pop('style', None)
+        rc = self._handle_style(style, reset=False)
+
+        def pop_keywords(sourcekeys, rckey):
+            ''' Helper function to handle input arguments '''
+            sourcekeys = sc.tolist(sourcekeys)
+            key = sourcekeys[0]  # Main key
+            value = None
+            changed = self.changed(key)
+            if changed:
+                value = self[key]
+            for k in sourcekeys:
+                kwvalue = kwargs.pop(k, None)
+                if kwvalue is not None:
+                    value = kwvalue
+            if value is not None:
+                rc[rckey] = value
+            return
+
+        # Handle special cases
+        pop_keywords('dpi', rckey='figure.dpi')
+        pop_keywords(['font', 'fontfamily', 'font_family'], rckey='font.family')
+        pop_keywords(['fontsize', 'font_size'], rckey='font.size')
+        pop_keywords('grid', rckey='axes.grid')
+        pop_keywords('facecolor', rckey='axes.facecolor')
+
+        # Handle other keywords
+        for key, value in kwargs.items():
+            if key not in pl.rcParams:
+                errormsg = f'Key "{key}" does not match any value in STIsim options or pl.rcParams'
+                raise sc.KeyNotFoundError(errormsg)
+            elif value is not None:
+                rc[key] = value
+
+        # Tidy up
+        if use:
+            return pl.style.use(sc.dcp(rc))
+        else:
+            return pl.style.context(sc.dcp(rc))
+
+    def use_style(self, **kwargs):
+        """
+        Shortcut to set stisim's current style as the global default.
+
+        **Example**::
+
+            ss.options.use_style() # Set stisim options as default
+            pl.figure()
+            pl.plot([1,3,7])
+
+            pl.style.use('seaborn-whitegrid') # to something else
+            pl.figure()
+            pl.plot([3,1,4])
+        """
+        return self.with_style(use=True, **kwargs)
+
+
+def load_fonts(folder=None, rebuild=False, verbose=False, **kwargs):
+    """
+    Helper function to load custom fonts for plotting -- (usually) not for the user.
+
+    Note: if fonts don't load, try running ``ss.settings.load_fonts(rebuild=True)``,
+    and/or rebooting the system.
+
+    Args:
+        folder (str): the folder to add fonts from
+        rebuild (bool): whether to rebuild the font cache
+        verbose (bool): whether to print out progress/errors
+    """
+
+    if folder is None:
+        folder = str(sc.thisdir(__file__, aspath=True) / 'data' / 'assets')
+    sc.fonts(add=folder, rebuild=rebuild, verbose=verbose, **kwargs)
+
+    # Try to find the font, and if it succeeds, update the styles
+    try:
+        name = 'Mulish'
+        fm.findfont(name, fallback_to_default=False)  # Raise an exception if the font isn't found
+        rc_simple['font.family'] = name  # Need to set both
+        rc_stisim['font.family'] = name
+        if verbose: print(f'Default STIsim font reset to "{name}"')
+    except Exception as E:
+        if verbose: print(f'Could not find font {name}: {str(E)}')
+
+    return
+
+
+# Create the options on module load, and load the fonts
+options = Options()
+
+
+# Default for precision
+# Used in various places throughout the code, generally as:
+#   import stisim.settings as sss
+#   arr = np.full(100, 0, hps.default_float)
+result_float = np.float64  # Always use float64 for results, for simplicity
+if options.precision == 32:
+    default_float = np.float32
+    default_int = np.int32
+elif options.precision == 64:  # pragma: no cover
+    default_float = np.float64
+    default_int = np.int64
+else:
+    raise NotImplementedError(f'Precision must be either 32 bit or 64 bit, not {options.precision}')
+
+
diff --git a/stisim/sim.py b/stisim/sim.py
index 3147988e..5558ec6e 100644
--- a/stisim/sim.py
+++ b/stisim/sim.py
@@ -1,88 +1,468 @@
-import sciris as sc
-import numpy as np
-from .people import make_people
-
-class Sim():
-
-    default_pars = sc.objdict(
-        n=1000,
-        npts=30,
-        dt=1
-    )
-
-    def __init__(self, people=None, modules=None, pars=None, interventions=None, analyzers=None):
-        # TODO - clearer options for time units?
-        #   ti - index self.ti = 0
-        #    t - floating point year e.g., 2020.5
-        #   dt - simulation primary step size (years per step)  # TODO - use-case for model-specific dt?
-        # date - '20200601'
-
-        self.ti = 0
-
-        self.pars = sc.dcp(self.default_pars)
-        if pars is not None:
-            self.pars.update(pars)
-        self.people = people if people is not None else make_people(self.pars['n'])
-        self.modules = sc.promotetolist(modules)
-        self.interventions = sc.promotetolist(interventions)
-        self.analyzers = sc.promotetolist(analyzers)
-        self.results = sc.objdict()
-
-    @property
-    def dt(self):
-        return self.pars['dt']
-
-    @property
-    def t(self):
-        return self.ti*self.dt
-
-    @property
-    def n(self):
-        return len(self.people)
-
-    @property
-    def tvec(self):
-        return np.arange(0,self.npts)*self.dt
-
-    @property
-    def npts(self):
-        return self.pars['npts']
-
-    def initialize(self):
-        self.people.initialize(self)
-
-        for module in self.modules:
-            module.initialize(self)
-
-        for intervention in self.interventions:
-            intervention.initialize(self)
-
-        for analyzer in self.analyzers:
-            analyzer.initialize(self)
-
-
-    def run(self):
-        self.initialize()
-        for i in range(self.pars.npts):
-            self.step()
-        self.people.finalize_results(self)
-
-        for module in self.modules:
-            module.finalize_results(self)
-
-    def step(self):
-
-        self.people.update_states_pre(self)
-
-        for intervention in self.interventions:
-            intervention.apply(self)
-
-        for module in self.modules:
-            module.transmit(self)
-
-        self.people.update_results(self)
-
-        for analyzer in self.analyzers:
-            analyzer.apply(self)
-
-        self.ti += self.pars.dt
+"""
+Define core Sim classes
+"""
+
+# Imports
+import numpy as np
+import sciris as sc
+from . import base as ssb
+from . import misc as ssm
+from . import settings as sss
+from . import utils as ssu
+from . import people as ssppl
+from . import parameters as sspar
+from . import interventions as ssi
+from . import analyzers as ssa
+from .results import Result
+
+
+# Define the model
+class Sim(ssb.BaseSim):
+
+    def __init__(self, pars=None, label=None, people=None, popdict=None, modules=None,
+                 networks=None, version=None, **kwargs):
+
+        # Set attributes
+        self.label = label  # The label/name of the simulation
+        self.created = None  # The datetime the sim was created
+        self.people = people  # People object
+        self.popdict = popdict  # popdict used to create people
+        self.networks = networks  # List of provided networks
+        self.modules = ssu.named_dict(modules)  # List of modules to simulate
+        self.results = sc.objdict()  # For storing results
+        self.summary = None  # For storing a summary of the results
+        self.initialized = False  # Whether initialization is complete
+        self.complete = False  # Whether a simulation has completed running
+        self.results_ready = False  # Whether results are ready
+        self._default_ver = version  # Default version of parameters used
+        self._orig_pars = None  # Store original parameters to optionally restore at the end of the simulation
+
+        # Time indexing
+        self.ti = None  # The time index, e.g. 0, 1, 2
+        self.yearvec = None
+        self.tivec = None
+        self.npts = None
+
+        # Make default parameters (using values from parameters.py)
+        default_pars = sspar.make_pars(version=version)  # Start with default pars
+        super().__init__(default_pars)  # Initialize and set the parameters as attributes
+
+        # Update parameters
+        self.update_pars(pars, **kwargs)  # Update the parameters
+
+        # Initialize other quantities
+        self.interventions = None
+        self.analyzers = None
+
+        return
+
+    @property
+    def dt(self):
+        return self.pars['dt']
+
+    @property
+    def year(self):
+        return self.yearvec[self.ti]
+
+    def initialize(self, popdict=None, reset=False, **kwargs):
+        """
+        Perform all initializations on the sim.
+        """
+        # Validation and initialization
+        self.ti = 0  # The current time index
+        self.validate_pars()  # Ensure parameters have valid values
+        self.validate_dt()
+        self.init_time_vecs()  # Initialize time vecs
+        ssu.set_seed(self['rand_seed'])  # Reset the random seed before the population is created
+
+        # Initialize the core sim components
+        self.init_people(popdict=popdict, reset=reset, **kwargs)  # Create all the people (the heaviest step)
+        self.init_networks()
+        self.init_results()
+        for module in self.modules.values():
+            module.initialize(self)
+        self.init_interventions()
+        self.init_analyzers()
+        self.validate_layer_pars()
+
+        # Reset the random seed to the default run seed, so that if the simulation is run with
+        # reset_seed=False right after initialization, it will still produce the same output
+        ssu.set_seed(self['rand_seed'] + 1)
+
+        # Final steps
+        self.initialized = True
+        self.complete = False
+        self.results_ready = False
+
+        return self
+
+    def layer_keys(self):
+        """
+        Attempt to retrieve the current network names
+        """
+        try:
+            keys = list(self.people['networks'].keys())
+        except:  # pragma: no cover
+            keys = []
+        return keys
+
+    def validate_layer_pars(self):
+        """
+        Check if there is a contact network
+        """
+
+        if self.people is not None:
+            modules = len(self.modules) > 0
+            pop_keys = set(self.people.networks.keys())
+            if modules and not len(pop_keys):
+                warnmsg = f'Warning: your simulation has {len(self.modules)} modules but no contact layers.'
+                ssm.warn(warnmsg, die=False)
+
+        return
+
+    def validate_dt(self):
+        """
+        Check that 1/dt is an integer value, otherwise results and time vectors will have mismatching shapes.
+        init_results explicitly makes this assumption by casting resfrequency = int(1/dt).
+        """
+        dt = self.dt
+        reciprocal = 1.0 / dt  # Compute the reciprocal of dt
+        if not reciprocal.is_integer():  # Check if reciprocal is not a whole (integer) number
+            # Round the reciprocal
+            reciprocal = int(reciprocal)
+            rounded_dt = 1.0 / reciprocal
+            self['dt'] = rounded_dt
+            if self['verbose']:
+                warnmsg = f"Warning: Provided time step dt: {dt} resulted in a non-integer number of steps per year. Rounded to {rounded_dt}."
+                print(warnmsg)
+
+    def validate_pars(self):
+        """
+        Some parameters can take multiple types; this makes them consistent.
+        """
+        # Handle n_agents
+        if self['n_agents'] is not None:
+            self['n_agents'] = int(self['n_agents'])
+        else:
+            if self.people is not None:
+                self['n_agents'] = len(self.people)
+            else:
+                if self.popdict is not None:
+                    self['n_agents'] = len(self.popdict)
+                else:
+                    errormsg = 'Must supply n_agents, a people object, or a popdict'
+                    raise ValueError(errormsg)
+
+        # Handle end and n_years
+        if self['end']:
+            self['n_years'] = int(self['end'] - self['start'])
+            if self['n_years'] <= 0:
+                errormsg = f"Number of years must be >0, but you supplied start={str(self['start'])} and " \
+                           f"end={str(self['end'])}, which gives n_years={self['n_years']}"
+                raise ValueError(errormsg)
+        else:
+            if self['n_years']:
+                self['end'] = self['start'] + self['n_years']
+            else:
+                errormsg = 'You must supply one of n_years and end."'
+                raise ValueError(errormsg)
+
+        # Handle verbose
+        if self['verbose'] == 'brief':
+            self['verbose'] = -1
+        if not sc.isnumber(self['verbose']):  # pragma: no cover
+            errormsg = f'Verbose argument should be either "brief", -1, or a float, not {type(self["verbose"])} "{self["verbose"]}"'
+            raise ValueError(errormsg)
+
+        return
+
+    def init_time_vecs(self):
+        """
+        Construct vectors things that keep track of time
+        """
+        self.yearvec = sc.inclusiverange(start=self['start'], stop=self['end'] + 1 - self['dt'],
+                                         step=self['dt'])  # Includes all the timepoints in the last year
+        self.npts = len(self.yearvec)
+        self.tivec = np.arange(self.npts)
+
+    def init_people(self, popdict=None, reset=False, verbose=None, **kwargs):
+        """
+        Initialize people within the sim
+        Sometimes the people are provided, in which case this just adds a few sim properties to them.
+        Other time people are not provided and this method makes them.
+        Args:
+            popdict         (any):  pre-generated people of various formats.
+            reset           (bool): whether to regenerate the people even if they already exist
+            verbose         (int):  detail to print
+            kwargs          (dict): passed to ss.make_people()
+        """
+
+        # Handle inputs
+        if verbose is None:
+            verbose = self['verbose']
+        if verbose > 0:
+            resetstr = ''
+            if self.people:
+                resetstr = ' (resetting people)' if reset else ' (warning: not resetting sim.people)'
+            print(f'Initializing sim{resetstr} with {self["n_agents"]:0n} agents')
+
+        # If people have not been supplied, make them
+        if self.people is None:
+            self.people = ssppl.People(self['n_agents'], kwargs)  # This just assigns UIDs and length
+
+        # If a popdict has not been supplied, we can make one from location data
+        if popdict is None:
+            if self['location'] is not None:
+                # Check where to get total_pop from
+                if self['total_pop'] is not None:  # If no pop_scale has been provided, try to get it from the location
+                    errormsg = 'You can either define total_pop explicitly or via the location, but not both'
+                    raise ValueError(errormsg)
+                total_pop, popdict = ssppl.make_popdict(n=self['n_agents'], location=self['location'], verbose=self['verbose'])
+
+            else:
+                if self['total_pop'] is not None:  # If no pop_scale has been provided, try to get it from the location
+                    total_pop = self['total_pop']
+                else:
+                    if self['pop_scale'] is not None:
+                        total_pop = self['pop_scale'] * self['n_agents']
+                    else:
+                        total_pop = self['n_agents']
+
+        self['total_pop'] = total_pop
+        if self['pop_scale'] is None:
+            self['pop_scale'] = total_pop / self['n_agents']
+
+        # Finish initialization
+        if not self.people.initialized:
+            self.people.initialize(popdict=popdict)  # Fully initialize the people
+
+        # Add time attributes
+        self.people.ti = self.ti
+        self.people.dt = self.dt
+
+        return self
+
+    def init_networks(self):
+        """ Initialize networks if these have been provided separately from the people """
+
+        # One possible workflow is that users will provide a location and a set of networks but not people.
+        # This means networks will be stored in self['networks'] and we'll need to copy them to the people.
+        if self.people.networks is None or len(self.people.networks) == 0:
+            if self['networks'] is not None:
+                self.people.networks = ssu.named_dict(self['networks'])
+
+        for key, network in self.people.networks.items():
+            if network.label is not None:
+                layer_name = network.label
+            else:
+                layer_name = key
+                network.label = layer_name
+            network.initialize(self.people)
+            self.people.networks[layer_name] = network
+
+        return
+
+    def init_results(self):
+        """
+        Create the main results structure.
+        """
+        # Make results
+        results = ssu.named_dict(
+            Result('births', None, self.npts, sss.default_float),
+            Result('deaths', None, self.npts, sss.default_float),
+            Result('n_alive', None, self.npts, sss.default_int),
+        )
+
+        # Final items
+        self.results = results
+        self.results_ready = False
+
+        return
+
+
+    def init_interventions(self):
+        """ Initialize and validate the interventions """
+
+        # Translate the intervention specs into actual interventions
+        for i, intervention in enumerate(self['interventions']):
+            if isinstance(intervention, type) and issubclass(intervention, ssi.Intervention):
+                intervention = intervention()  # Convert from a class to an instance of a class
+            if isinstance(intervention, ssi.Intervention):
+                intervention.initialize(self)
+                self.interventions += intervention
+            elif callable(intervention):
+                self.interventions += intervention
+            else:
+                errormsg = f'Intervention {intervention} does not seem to be a valid intervention: must be a function or Intervention subclass'
+                raise TypeError(errormsg)
+
+        return
+
+    def init_analyzers(self):
+        """ Initialize the analyzers """
+
+        self.analyzers = sc.autolist()
+
+        # Interpret analyzers
+        for ai, analyzer in enumerate(self['analyzers']):
+            if isinstance(analyzer, type) and issubclass(analyzer, ssa.Analyzer):
+                analyzer = analyzer()  # Convert from a class to an instance of a class
+            if not (isinstance(analyzer, ssa.Analyzer) or callable(analyzer)):
+                errormsg = f'Analyzer {analyzer} does not seem to be a valid analyzer: must be a function or Analyzer subclass'
+                raise TypeError(errormsg)
+            self.analyzers += analyzer  # Add it in
+
+        for analyzer in self.analyzers:
+            if isinstance(analyzer, ssa.Analyzer):
+                analyzer.initialize(self)
+
+        return
+
+    def step(self):
+        """ Step through time and update values """
+
+        # Set the time and if we have reached the end of the simulation, then do nothing
+        if self.complete:
+            raise AlreadyRunError('Simulation already complete (call sim.initialize() to re-run)')
+
+        # Update states, modules, partnerships
+        self.people.update_states(sim=self)  # This runs modules
+        # self.update_connectors()  # TODO: add this when ready
+
+        for module in self.modules.values():
+            module.make_new_cases(self)
+            module.update_results(self)
+
+        # Tidy up
+        self.ti += 1
+        if self.ti == self.npts:
+            self.complete = True
+
+        return
+
+    def run(self, until=None, reset_seed=True, verbose=None):
+        """ Run the model once """
+
+        # Initialization steps
+        T = sc.timer()
+        if not self.initialized:
+            self.initialize()
+            self._orig_pars = sc.dcp(self.pars)  # Create a copy of the parameters to restore after the run
+
+        if verbose is None:
+            verbose = self['verbose']
+
+        if reset_seed:
+            # Reset the RNG. The primary use case (and why it defaults to True) is to ensure that
+            #
+            # >>> sim0.initialize()
+            # >>> sim0.run()
+            # >>> sim1.initialize()
+            # >>> sim1.run()
+            #
+            # produces the same output as
+            #
+            # >>> sim0.initialize()
+            # >>> sim1.initialize()
+            # >>> sim0.run()
+            # >>> sim1.run()
+            #
+            # The seed is offset by 1 to avoid drawing the same random numbers as those used for population generation,
+            # otherwise the first set of random numbers in the model (e.g., deaths) will be correlated with the first
+            # set of random numbers drawn in population generation (e.g., sex)
+            ssu.set_seed(self['rand_seed'] + 1)
+
+        # Check for AlreadyRun errors
+        errormsg = None
+        if until is None: until = self.npts
+        if until > self.npts:
+            errormsg = f'Requested to run until t={until} but the simulation end is ti={self.npts}'
+        if self.ti >= until:  # NB. At the start, self.t is None so this check must occur after initialization
+            errormsg = f'Simulation is currently at t={self.ti}, requested to run until ti={until} which has already been reached'
+        if self.complete:
+            errormsg = 'Simulation is already complete (call sim.initialize() to re-run)'
+        if errormsg:
+            raise AlreadyRunError(errormsg)
+
+        # Main simulation loop
+        while self.ti < until:
+
+            # Check if we were asked to stop
+            elapsed = T.toc(output=True)
+            if self['timelimit'] and elapsed > self['timelimit']:
+                sc.printv(
+                    f"Time limit ({self['timelimit']} s) exceeded; call sim.finalize() to compute results if desired",
+                    1, verbose)
+                return
+            elif self['stopping_func'] and self['stopping_func'](self):
+                sc.printv(
+                    "Stopping function terminated the simulation; call sim.finalize() to compute results if desired", 1,
+                    verbose)
+                return
+
+            # Print progress
+            if verbose:
+                simlabel = f'"{self.label}": ' if self.label else ''
+                string = f'  Running {simlabel}{self.yearvec[self.ti]:0.1f} ({self.ti:2.0f}/{self.npts}) ({elapsed:0.2f} s) '
+                if verbose >= 2:
+                    sc.heading(string)
+                elif verbose > 0:
+                    if not (self.ti % int(1.0 / verbose)):
+                        sc.progressbar(self.ti + 1, self.npts, label=string, length=20, newline=True)
+
+            # Actually run the model
+            self.step()
+
+        # If simulation reached the end, finalize the results
+        if self.complete:
+            self.finalize(verbose=verbose)
+            sc.printv(f'Run finished after {elapsed:0.2f} s.\n', 1, verbose)
+
+        return self
+
+    def finalize(self, verbose=None):
+        """ Compute final results """
+
+        if self.results_ready:
+            # Because the results are rescaled in-place, finalizing the sim cannot be run more than once or
+            # otherwise the scale factor will be applied multiple times
+            raise AlreadyRunError('Simulation has already been finalized')
+
+        # Final settings
+        self.results_ready = True  # Set this first so self.summary() knows to print the results
+        self.ti -= 1  # During the run, this keeps track of the next step; restore this be the final day of the sim
+
+        # Perform calculations on results
+        # self.compute_results(verbose=verbose) # Calculate the rest of the results
+        self.results = sc.objdict(
+            self.results)  # Convert results to a odicts/objdict to allow e.g. sim.results.diagnoses
+
+        return
+
+    def update_connectors(self):
+        """ To do: move this to people? """
+        if len(self.modules) > 1:
+            connectors = self['connectors']
+            if len(connectors) > 0:
+                for connector in connectors:
+                    if callable(connector):
+                        connector(self)
+                    else:
+                        warnmsg = f'Connector must be a callable function'
+                        ssm.warn(warnmsg, die=True)
+            elif self.ti == 0:  # only raise warning on first timestep
+                warnmsg = f'No connectors in sim'
+                ssm.warn(warnmsg, die=False)
+            else:
+                return
+        return
+
+
+class AlreadyRunError(RuntimeError):
+    """
+    This error is raised if a simulation is run in such a way that no timesteps
+    will be taken. This error is a distinct type so that it can be safely caught
+    and ignored if required, but it is anticipated that most of the time, calling
+    :py:func:`Sim.run` and not taking any timesteps, would be an inadvertent error.
+    """
+    pass
diff --git a/stisim/utils.py b/stisim/utils.py
index f0a0dd47..bcdee6d5 100644
--- a/stisim/utils.py
+++ b/stisim/utils.py
@@ -1,732 +1,743 @@
-'''
-Numerical utilities for running hpvsim.
-'''
-
-#%% Housekeeping
-
-import numpy as np # For numerics
-import sciris as sc # For additional utilities
-
-
-# What functions are externally visible -- note, this gets populated in each section below
-__all__ = []
-
-
-#%% The core functions
-
-
-def unique(arr):
-    '''
-    Find the unique elements and counts in an array.
-    Equivalent to np.unique(return_counts=True) but ~5x faster, and
-    only works for arrays of positive integers.
-    '''
-    counts = np.bincount(arr.ravel())
-    unique = np.flatnonzero(counts)
-    counts = counts[unique]
-    return unique, counts
-
-
-def isin(arr, search_inds):
-    ''' Find search_inds in arr. Like np.isin() but faster '''
-    n = len(arr)
-    result = np.full(n, False)
-    set_search_inds = set(search_inds)
-    for i in range(n):
-        if arr[i] in set_search_inds:
-            result[i] = True
-    return result
-
-
-def findinds(arr, vals):
-    ''' Finds indices of vals in arr, accounting for repeats '''
-    return isin(arr,vals).nonzero()[0]
-
-
-def find_contacts(p1, p2, inds): # pragma: no cover
-    """
-    Numba for Layer.find_contacts()
-
-    A set is returned here rather than a sorted array so that custom tracing interventions can efficiently
-    add extra people. For a version with sorting by default, see Layer.find_contacts(). Indices must be
-    an int64 array since this is what's returned by true() etc. functions by default.
-    """
-    pairing_partners = set()
-    inds = set(inds)
-    for i in range(len(p1)):
-        if p1[i] in inds:
-            pairing_partners.add(p2[i])
-        if p2[i] in inds:
-            pairing_partners.add(p1[i])
-    return pairing_partners
-
-
-def logf1(x, k, ttc=25):
-    '''
-    Logistic function passing through (0,0) and (ttc,1).
-    Accepts 1 parameter which determines the growth rate.
-    '''
-    return logf3(x, k, 0, 1, ttc=ttc)
-
-
-def get_asymptotes(k, x_infl, s, ttc=25):
-    '''
-    Get upper asymptotes for logistic functions
-    '''
-    term1 = (1 + np.exp(k*(x_infl-ttc)))**s # Note, this is 1 for most parameter combinations
-    term2 = (1 + np.exp(k*x_infl))**s
-    u_asymp_num = term1*(1-term2)
-    u_asymp_denom = term1 - term2
-    u_asymp = u_asymp_num / u_asymp_denom
-    l_asymp = term1 / (term1 - term2)
-    return l_asymp, u_asymp
-
-
-def logf3(x, k, x_infl, s, ttc=25):
-    '''
-    Logistic function passing through (0,0) and (ttc,1).
-    This version is derived from the 5-parameter version here: https://www.r-bloggers.com/2019/11/five-parameters-logistic-regression/
-    However, since it's constrained to pass through 2 points, there are 3 free parameters remaining.
-    Args:
-         k: growth rate, equivalent to b in https://www.r-bloggers.com/2019/11/five-parameters-logistic-regression/
-         x_infl: a location parameter, equivalent to C in https://www.r-bloggers.com/2019/11/five-parameters-logistic-regression/
-         s: asymmetry parameter, equivalent to s in https://www.r-bloggers.com/2019/11/five-parameters-logistic-regression/
-         ttc (time to cancer): x value for which the curve passess through 1. For x values beyond this, the function returns 1
-    '''
-    l_asymp, u_asymp = get_asymptotes(k, x_infl, s, ttc)
-    return np.minimum(1, l_asymp + (u_asymp-l_asymp)/(1+np.exp(k*(x_infl-x)))**s)
-
-
-def logf2(x, k, x_infl, ttc=25):
-    '''
-    Logistic function constrained to pass through (0,0) and (ttc,1).
-    This version is derived from the 5-parameter version here: https://www.r-bloggers.com/2019/11/five-parameters-logistic-regression/
-    Since it's constrained to pass through 2 points, there are 3 free parameters remaining, and this verison fixes s=1
-    Args:
-         k: growth rate, equivalent to b in https://www.r-bloggers.com/2019/11/five-parameters-logistic-regression/
-         x_infl: point of inflection, equivalent to C in https://www.r-bloggers.com/2019/11/five-parameters-logistic-regression/
-         ttc (time to cancer): x value for which the curve passess through 1. For x values beyond this, the function returns 1
-    '''
-    return logf3(x, k, x_infl, s=1, ttc=ttc)
-
-
-def invlogf3(y, k, x_infl, s, ttc=25):
-    '''
-    Inverse of logf3; see definition there for arguments
-    '''
-    l_asymp, u_asymp = get_asymptotes(k, x_infl, s, ttc)
-    part1 = np.log((u_asymp-l_asymp)/(y-l_asymp))/s
-    part2 = np.log(np.exp(part1)-1)
-    final = 1/k * (k*x_infl - part2)
-    return final
-
-
-def invlogf2(y, k, x_infl, ttc=25):
-    '''
-    Inverse of logf2; see definition there for arguments
-    '''
-    return invlogf3(y, k, x_infl, 1, ttc=ttc)
-
-
-def invlogf1(y, k, ttc=25):
-    '''
-    The inverse of the concave part of a logistic function, with point of inflexion at 0,0
-    and upper asymptote at 1. Accepts 1 parameter which determines the growth rate.
-    '''
-    return invlogf3(y, k, 0, 1, ttc=ttc)
-
-
-def indef_int_logf2(x, k, x_infl, ttc=25):
-    '''
-    Indefinite integral of logf2; see definition there for arguments
-    '''
-    num = np.exp(-x_infl*k)*(np.exp(k*ttc)+np.exp(x_infl*k))*((np.exp(x_infl*k)+1)*np.log(np.exp(k*x)+np.exp(x_infl*k))-k*x)
-    denom = k*(np.exp(k*ttc)-1)
-    return num/denom
-
-
-def intlogf2(upper, k, x_infl, ttc=25):
-    '''
-    Integral of logf2 between 0 and the limit given by upper
-    '''
-    # Find the upper limits not including the part past time to cancer
-    exceeding_ttc_inds = (upper > ttc).nonzero()
-    lims_to_find = np.minimum(ttc, upper)
-
-    # Take the integral
-    val_at_0 = indef_int_logf2(0, k, x_infl, ttc)
-    val_at_lim = indef_int_logf2(lims_to_find, k, x_infl, ttc)
-    integral = val_at_lim - val_at_0
-
-    # Deal with those whose duration of infection exceeds the time to cancer
-    # Note, another option would be to set their transformation probability to 1
-    excess_integral = upper[exceeding_ttc_inds] - ttc
-    integral[exceeding_ttc_inds] += excess_integral
-
-    return integral
-
-
-def indef_int_logf1(x, k, ttc=25):
-    '''
-    Indefinite integral of logf1; see definition there for arguments
-    '''
-    return indef_int_logf2(x, k, 0, ttc=ttc)
-
-
-def intlogf1(upper, k, ttc=25):
-    '''
-    Integral of logf1 between 0 and the limit given by upper
-    '''
-    return intlogf2(x, k, 0, ttc=ttc)
-
-
-def transform_prob(tp, dysp):
-    '''
-    Returns transformation probability given dysplasia
-    Using formula for half an ellipsoid:
-        V = 1/2 * 4/3 * pi * a*b*c
-          = 2 * a*b*c
-          = 2* dysp * (dysp/2)**2, assuming that b = c = 1/2 a
-          = 1/2 * dysp**3
-    '''
-    # return 1-np.power(1-tp, ((dysp*100)**2))
-    return 1-np.power(1-tp, 0.5*((dysp*100)**3))
-
-
-#%% Sampling and seed methods
-
-__all__ += ['sample', 'get_pdf', 'set_seed']
-
-
-def sample(dist=None, par1=None, par2=None, size=None, **kwargs):
-    '''
-    Draw a sample from the distribution specified by the input. The available
-    distributions are:
-
-    - 'uniform'       : uniform distribution from low=par1 to high=par2; mean is equal to (par1+par2)/2
-    - 'normal'        : normal distribution with mean=par1 and std=par2
-    - 'lognormal'     : lognormal distribution with mean=par1 and std=par2 (parameters are for the lognormal distribution, *not* the underlying normal distribution)
-    - 'normal_pos'    : right-sided normal distribution (i.e. only positive values), with mean=par1 and std=par2 *of the underlying normal distribution*
-    - 'normal_int'    : normal distribution with mean=par1 and std=par2, returns only integer values
-    - 'lognormal_int' : lognormal distribution with mean=par1 and std=par2, returns only integer values
-    - 'poisson'       : Poisson distribution with rate=par1 (par2 is not used); mean and variance are equal to par1
-    - 'neg_binomial'  : negative binomial distribution with mean=par1 and k=par2; converges to Poisson with k=∞
-    - 'beta'          : beta distribution with alpha=par1 and beta=par2;
-    - 'gamma'         : gamma distribution with shape=par1 and scale=par2;
-
-    Args:
-        dist (str):   the distribution to sample from
-        par1 (float): the "main" distribution parameter (e.g. mean)
-        par2 (float): the "secondary" distribution parameter (e.g. std)
-        size (int):   the number of samples (default=1)
-        kwargs (dict): passed to individual sampling functions
-
-    Returns:
-        A length N array of samples
-
-    **Examples**::
-
-        hpv.sample() # returns Unif(0,1)
-        hpv.sample(dist='normal', par1=3, par2=0.5) # returns Normal(μ=3, σ=0.5)
-        hpv.sample(dist='lognormal_int', par1=5, par2=3) # returns a lognormally distributed set of values with mean 5 and std 3
-
-    Notes:
-        Lognormal distributions are parameterized with reference to the underlying normal distribution (see:
-        https://docs.scipy.org/doc/numpy-1.14.0/reference/generated/numpy.random.lognormal.html), but this
-        function assumes the user wants to specify the mean and std of the lognormal distribution.
-
-        Negative binomial distributions are parameterized with reference to the mean and dispersion parameter k
-        (see: https://en.wikipedia.org/wiki/Negative_binomial_distribution). The r parameter of the underlying
-        distribution is then calculated from the desired mean and k. For a small mean (~1), a dispersion parameter
-        of ∞ corresponds to the variance and standard deviation being equal to the mean (i.e., Poisson). For a
-        large mean (e.g. >100), a dispersion parameter of 1 corresponds to the standard deviation being equal to
-        the mean.
-    '''
-
-    # Some of these have aliases, but these are the "official" names
-    choices = [
-        'uniform',
-        'normal',
-        'normal_pos',
-        'normal_int',
-        'lognormal',
-        'lognormal_int',
-        'poisson',
-        'neg_binomial',
-        'beta',
-        'gamma',
-    ]
-
-    # Ensure it's an integer
-    if size is not None and not isinstance(size, tuple):
-        size = int(size)
-
-    # Compute distribution parameters and draw samples
-    # NB, if adding a new distribution, also add to choices above
-    if   dist in ['unif', 'uniform']: samples = np.random.uniform(low=par1, high=par2, size=size, **kwargs)
-    elif dist in ['norm', 'normal']:  samples = np.random.normal(loc=par1, scale=par2, size=size, **kwargs)
-    elif dist == 'normal_pos':        samples = np.abs(np.random.normal(loc=par1, scale=par2, size=size, **kwargs))
-    elif dist == 'normal_int':        samples = np.round(np.abs(np.random.normal(loc=par1, scale=par2, size=size, **kwargs)))
-    elif dist == 'poisson':           samples = n_poisson(rate=par1, n=size, **kwargs) # Use Numba version below for speed
-    elif dist == 'neg_binomial':      samples = n_neg_binomial(rate=par1, dispersion=par2, n=size, **kwargs) # Use custom version below
-    elif dist == 'beta':              samples = np.random.beta(a=par1, b=par2, size=size, **kwargs)
-    elif dist == 'gamma':             samples = np.random.gamma(shape=par1, scale=par2, size=size, **kwargs)
-    elif dist in ['lognorm', 'lognormal', 'lognorm_int', 'lognormal_int']:
-        if (sc.isnumber(par1) and par1>0) or (sc.checktype(par1,'arraylike') and (par1>0).all()):
-            mean  = np.log(par1**2 / np.sqrt(par2**2 + par1**2)) # Computes the mean of the underlying normal distribution
-            sigma = np.sqrt(np.log(par2**2/par1**2 + 1)) # Computes sigma for the underlying normal distribution
-            samples = np.random.lognormal(mean=mean, sigma=sigma, size=size, **kwargs)
-        else:
-            samples = np.zeros(size)
-        if '_int' in dist:
-            samples = np.round(samples)
-    elif dist == 'beta_mean': # Calculate a and b using mean (par1) and variance (par2) https://stats.stackexchange.com/questions/12232/calculating-the-parameters-of-a-beta-distribution-using-the-mean-and-variance
-        a       = ((1 - par1)/par2 - 1/par1) * par1**2
-        b       = a * (1 / par1 - 1)
-        samples = np.random.beta(a=a, b=b, size=size, **kwargs)
-    else:
-        errormsg = f'The selected distribution "{dist}" is not implemented; choices are: {sc.newlinejoin(choices)}'
-        raise NotImplementedError(errormsg)
-
-    return samples
-
-
-
-def get_pdf(dist=None, par1=None, par2=None):
-    '''
-    Return a probability density function for the specified distribution. This
-    is used for example by test_num to retrieve the distribution of times from
-    symptom-to-swab for testing. For example, for Washington State, these values
-    are dist='lognormal', par1=10, par2=170.
-    '''
-    import scipy.stats as sps # Import here since slow
-
-    choices = [
-        'none',
-        'uniform',
-        'lognormal',
-    ]
-
-    if dist in ['None', 'none', None]:
-        return None
-    elif dist == 'uniform':
-        pdf = sps.uniform(loc=par1, scale=par2)
-    elif dist == 'lognormal':
-        mean  = np.log(par1**2 / np.sqrt(par2 + par1**2)) # Computes the mean of the underlying normal distribution
-        sigma = np.sqrt(np.log(par2/par1**2 + 1)) # Computes sigma for the underlying normal distribution
-        pdf   = sps.lognorm(sigma, loc=-0.5, scale=np.exp(mean))
-    else:
-        choicestr = '\n'.join(choices)
-        errormsg = f'The selected distribution "{dist}" is not implemented; choices are: {choicestr}'
-        raise NotImplementedError(errormsg)
-
-    return pdf
-
-
-def set_seed(seed=None):
-    '''
-    Reset the random seed. This function also resets Python's built-in random
-    number generated.
-
-    Args:
-        seed (int): the random seed
-    '''
-    # Dies if a float is given
-    if seed is not None:
-        seed = int(seed)
-    np.random.seed(seed) # If None, reinitializes it
-    return
-
-
-#%% Probabilities -- mostly not jitted since performance gain is minimal
-
-__all__ += ['n_binomial', 'binomial_filter', 'binomial_arr', 'n_multinomial',
-            'poisson', 'n_poisson', 'n_neg_binomial', 'choose', 'choose_r', 'choose_w']
-
-def n_binomial(prob, n):
-    '''
-    Perform multiple binomial (Bernolli) trials
-
-    Args:
-        prob (float): probability of each trial succeeding
-        n (int): number of trials (size of array)
-
-    Returns:
-        Boolean array of which trials succeeded
-
-    **Example**::
-
-        outcomes = hpv.n_binomial(0.5, 100) # Perform 100 coin-flips
-    '''
-    return np.random.random(n) < prob
-
-
-def binomial_filter(prob, arr):
-    '''
-    Binomial "filter" -- the same as n_binomial, except return
-    the elements of arr that succeeded.
-
-    Args:
-        prob (float): probability of each trial succeeding
-        arr (array): the array to be filtered
-
-    Returns:
-        Subset of array for which trials succeeded
-
-    **Example**::
-
-        inds = hpv.binomial_filter(0.5, np.arange(20)**2) # Return which values out of the (arbitrary) array passed the coin flip
-    '''
-    return arr[(np.random.random(len(arr)) < prob).nonzero()[0]]
-
-
-def binomial_arr(prob_arr):
-    '''
-    Binomial (Bernoulli) trials each with different probabilities.
-
-    Args:
-        prob_arr (array): array of probabilities
-
-    Returns:
-         Boolean array of which trials on the input array succeeded
-
-    **Example**::
-
-        outcomes = hpv.binomial_arr([0.1, 0.1, 0.2, 0.2, 0.8, 0.8]) # Perform 6 trials with different probabilities
-    '''
-    return np.random.random(prob_arr.shape) < prob_arr
-
-
-def n_multinomial(probs, n): # No speed gain from Numba
-    '''
-    An array of multinomial trials.
-
-    Args:
-        probs (array): probability of each outcome, which usually should sum to 1
-        n (int): number of trials
-
-    Returns:
-        Array of integer outcomes
-
-    **Example**::
-
-        outcomes = hpv.n_multinomial(np.ones(6)/6.0, 50)+1 # Return 50 die-rolls
-    '''
-    return np.searchsorted(np.cumsum(probs), np.random.random(n))
-
-
-def poisson(rate):
-    '''
-    A Poisson trial.
-
-    Args:
-        rate (float): the rate of the Poisson process
-
-    **Example**::
-
-        outcome = hpv.poisson(100) # Single Poisson trial with mean 100
-    '''
-    return np.random.poisson(rate, 1)[0]
-
-
-def n_poisson(rate, n):
-    '''
-    An array of Poisson trials.
-
-    Args:
-        rate (float): the rate of the Poisson process (mean)
-        n (int): number of trials
-
-    **Example**::
-
-        outcomes = hpv.n_poisson(100, 20) # 20 Poisson trials with mean 100
-    '''
-    return np.random.poisson(rate, n)
-
-
-def n_neg_binomial(rate, dispersion, n, step=1): # Numba not used due to incompatible implementation
-    '''
-    An array of negative binomial trials. See hpv.sample() for more explanation.
-
-    Args:
-        rate (float): the rate of the process (mean, same as Poisson)
-        dispersion (float):  dispersion parameter; lower is more dispersion, i.e. 0 = infinite, ∞ = Poisson
-        n (int): number of trials
-        step (float): the step size to use if non-integer outputs are desired
-
-    **Example**::
-
-        outcomes = hpv.n_neg_binomial(100, 1, 50) # 50 negative binomial trials with mean 100 and dispersion roughly equal to mean (large-mean limit)
-        outcomes = hpv.n_neg_binomial(1, 100, 20) # 20 negative binomial trials with mean 1 and dispersion still roughly equal to mean (approximately Poisson)
-    '''
-    nbn_n = dispersion
-    nbn_p = dispersion/(rate/step + dispersion)
-    samples = np.random.negative_binomial(n=nbn_n, p=nbn_p, size=n)*step
-    return samples
-
-
-def choose(max_n, n):
-    '''
-    Choose a subset of items (e.g., people) without replacement.
-
-    Args:
-        max_n (int): the total number of items
-        n (int): the number of items to choose
-
-    **Example**::
-
-        choices = hpv.choose(5, 2) # choose 2 out of 5 people with equal probability (without repeats)
-    '''
-    return np.random.choice(max_n, n, replace=False)
-
-
-def choose_r(max_n, n):
-    '''
-    Choose a subset of items (e.g., people), with replacement.
-
-    Args:
-        max_n (int): the total number of items
-        n (int): the number of items to choose
-
-    **Example**::
-
-        choices = hpv.choose_r(5, 10) # choose 10 out of 5 people with equal probability (with repeats)
-    '''
-    return np.random.choice(max_n, n, replace=True)
-
-
-def choose_w(probs, n, unique=True): # No performance gain from Numba
-    '''
-    Choose n items (e.g. people), each with a probability from the distribution probs.
-
-    Args:
-        probs (array): list of probabilities, should sum to 1
-        n (int): number of samples to choose
-        unique (bool): whether or not to ensure unique indices
-
-    **Example**::
-
-        choices = hpv.choose_w([0.2, 0.5, 0.1, 0.1, 0.1], 2) # choose 2 out of 5 people with nonequal probability.
-    '''
-    probs = np.array(probs)
-    n_choices = len(probs)
-    n_samples = int(n)
-    probs_sum = probs.sum()
-    if probs_sum: # Weight is nonzero, rescale
-        probs = probs/probs_sum
-    else: # Weights are all zero, choose uniformly
-        probs = np.ones(n_choices)/n_choices
-    return np.random.choice(n_choices, n_samples, p=probs, replace=not(unique))
-
-
-
-#%% Simple array operations
-
-__all__ += ['true',   'false',   'defined',   'undefined',
-            'itrue',  'ifalse',  'idefined',  'iundefined',
-            'itruei', 'ifalsei', 'idefinedi', 'iundefinedi',
-            'dtround', 'find_cutoff']
-
-
-def true(arr):
-    '''
-    Returns the indices of the values of the array that are true: just an alias
-    for arr.nonzero()[0].
-
-    Args:
-        arr (array): any array
-
-    **Example**::
-
-        inds = hpv.true(np.array([1,0,0,1,1,0,1])) # Returns array([0, 3, 4, 6])
-    '''
-    return arr.nonzero()[-1]
-
-
-def false(arr):
-    '''
-    Returns the indices of the values of the array that are false.
-
-    Args:
-        arr (array): any array
-
-    **Example**::
-
-        inds = hpv.false(np.array([1,0,0,1,1,0,1]))
-    '''
-    return np.logical_not(arr).nonzero()[-1]
-
-
-def defined(arr):
-    '''
-    Returns the indices of the values of the array that are not-nan.
-
-    Args:
-        arr (array): any array
-
-    **Example**::
-
-        inds = hpv.defined(np.array([1,np.nan,0,np.nan,1,0,1]))
-    '''
-    return (~np.isnan(arr)).nonzero()[-1]
-
-
-def undefined(arr):
-    '''
-    Returns the indices of the values of the array that are not-nan.
-
-    Args:
-        arr (array): any array
-
-    **Example**::
-
-        inds = hpv.defined(np.array([1,np.nan,0,np.nan,1,0,1]))
-    '''
-    return np.isnan(arr).nonzero()[-1]
-
-
-def itrue(arr, inds):
-    '''
-    Returns the indices that are true in the array -- name is short for indices[true]
-
-    Args:
-        arr (array): a Boolean array, used as a filter
-        inds (array): any other array (usually, an array of indices) of the same size
-
-    **Example**::
-
-        inds = hpv.itrue(np.array([True,False,True,True]), inds=np.array([5,22,47,93]))
-    '''
-    return inds[arr]
-
-
-def ifalse(arr, inds):
-    '''
-    Returns the indices that are true in the array -- name is short for indices[false]
-
-    Args:
-        arr (array): a Boolean array, used as a filter
-        inds (array): any other array (usually, an array of indices) of the same size
-
-    **Example**::
-
-        inds = hpv.ifalse(np.array([True,False,True,True]), inds=np.array([5,22,47,93]))
-    '''
-    return inds[np.logical_not(arr)]
-
-
-def idefined(arr, inds):
-    '''
-    Returns the indices that are defined in the array -- name is short for indices[defined]
-
-    Args:
-        arr (array): any array, used as a filter
-        inds (array): any other array (usually, an array of indices) of the same size
-
-    **Example**::
-
-        inds = hpv.idefined(np.array([3,np.nan,np.nan,4]), inds=np.array([5,22,47,93]))
-    '''
-    return inds[~np.isnan(arr)]
-
-
-def iundefined(arr, inds):
-    '''
-    Returns the indices that are undefined in the array -- name is short for indices[undefined]
-
-    Args:
-        arr (array): any array, used as a filter
-        inds (array): any other array (usually, an array of indices) of the same size
-
-    **Example**::
-
-        inds = hpv.iundefined(np.array([3,np.nan,np.nan,4]), inds=np.array([5,22,47,93]))
-    '''
-    return inds[np.isnan(arr)]
-
-
-
-def itruei(arr, inds):
-    '''
-    Returns the indices that are true in the array -- name is short for indices[true[indices]]
-
-    Args:
-        arr (array): a Boolean array, used as a filter
-        inds (array): an array of indices for the original array
-
-    **Example**::
-
-        inds = hpv.itruei(np.array([True,False,True,True,False,False,True,False]), inds=np.array([0,1,3,5]))
-    '''
-    return inds[arr[inds]]
-
-
-def ifalsei(arr, inds):
-    '''
-    Returns the indices that are false in the array -- name is short for indices[false[indices]]
-
-    Args:
-        arr (array): a Boolean array, used as a filter
-        inds (array): an array of indices for the original array
-
-    **Example**::
-
-        inds = hpv.ifalsei(np.array([True,False,True,True,False,False,True,False]), inds=np.array([0,1,3,5]))
-    '''
-    return inds[np.logical_not(arr[inds])]
-
-
-def idefinedi(arr, inds):
-    '''
-    Returns the indices that are defined in the array -- name is short for indices[defined[indices]]
-
-    Args:
-        arr (array): any array, used as a filter
-        inds (array): an array of indices for the original array
-
-    **Example**::
-
-        inds = hpv.idefinedi(np.array([4,np.nan,0,np.nan,np.nan,4,7,4,np.nan]), inds=np.array([0,1,3,5]))
-    '''
-    return inds[~np.isnan(arr[inds])]
-
-
-def iundefinedi(arr, inds):
-    '''
-    Returns the indices that are undefined in the array -- name is short for indices[defined[indices]]
-
-    Args:
-        arr (array): any array, used as a filter
-        inds (array): an array of indices for the original array
-
-    **Example**::
-
-        inds = hpv.iundefinedi(np.array([4,np.nan,0,np.nan,np.nan,4,7,4,np.nan]), inds=np.array([0,1,3,5]))
-    '''
-    return inds[np.isnan(arr[inds])]
-
-
-def dtround(arr, dt, ceil=True):
-    '''
-    Rounds the values in the array to the nearest timestep
-
-    Args:
-        arr (array): any array
-        dt  (float): float, usually representing a timestep in years
-
-    **Example**::
-
-        dtround = hpv.dtround(np.array([0.23,0.61,20.53])) # Returns array([0.2, 0.6, 20.6])
-        dtround = hpv.dtround(np.array([0.23,0.61,20.53]),ceil=True) # Returns array([0.4, 0.8, 20.6])
-    '''
-    if ceil:
-        return np.ceil(arr * (1/dt)) / (1/dt)
-    else:
-        return np.round(arr * (1/dt)) / (1/dt)
-
-
-def find_cutoff(duration_cutoffs, duration):
-    '''
-    Find which duration bin each ind belongs to.
-    '''
-    return np.nonzero(duration_cutoffs <= duration)[0][-1]  # Index of the duration bin to use
+"""
+Numerical utilities
+"""
+
+# %% Housekeeping
+
+import numpy as np  # For numerics
+import sciris as sc  # For additional utilities
+
+# What functions are externally visible -- note, this gets populated in each section below
+__all__ = []
+
+
+# %% Helper functions
+__all__ += ['named_dict', 'omerge']
+
+
+def named_dict(arg=None, *args):
+    """ Create an objdict based on object names """
+    args = sc.mergelists(arg, args, coerce='tuple')  # TODO: Must be a better way, I just forgot
+    return sc.objdict({arg.name: arg for arg in args})  # TODO: use dictobj instead for performance, probably
+
+
+def omerge(*args, **kwargs):
+    """ Merge things into an objdict """
+    return sc.objdict(sc.mergedicts(*args, **kwargs))
+
+
+# %% The core functions
+def unique(arr):
+    """
+    Find the unique elements and counts in an array.
+    Equivalent to np.unique(return_counts=True) but ~5x faster, and
+    only works for arrays of positive integers.
+    """
+    counts = np.bincount(arr.ravel())
+    unique = np.flatnonzero(counts)
+    counts = counts[unique]
+    return unique, counts
+
+
+def find_contacts(p1, p2, inds):  # pragma: no cover
+    """
+    Variation on Network.find_contacts() that avoids sorting.
+
+    A set is returned here rather than a sorted array so that custom tracing interventions can efficiently
+    add extra people. For a version with sorting by default, see Network.find_contacts(). Indices must be
+    an int64 array since this is what's returned by true() etc. functions by default.
+    """
+    pairing_partners = set()
+    inds = set(inds)
+    for i in range(len(p1)):
+        if p1[i] in inds:
+            pairing_partners.add(p2[i])
+        if p2[i] in inds:
+            pairing_partners.add(p1[i])
+    return pairing_partners
+
+
+# %% Sampling and seed methods
+
+__all__ += ['sample', 'get_pdf', 'set_seed']
+
+
+def sample(dist=None, par1=None, par2=None, size=None, **kwargs):
+    """
+    Draw a sample from the distribution specified by the input. The available
+    distributions are:
+
+    - 'uniform'       : uniform from low=par1 to high=par2; mean is equal to (par1+par2)/2
+    - 'choice'        : par1=array of choices, par2=probability of each choice
+    - 'normal'        : normal with mean=par1 and std=par2
+    - 'lognormal'     : lognormal with mean=par1, std=par2 (parameters are for the lognormal, not the underlying normal)
+    - 'normal_pos'    : right-sided normal (i.e. only +ve values), with mean=par1, std=par2 of the underlying normal
+    - 'normal_int'    : normal distribution with mean=par1 and std=par2, returns only integer values
+    - 'lognormal_int' : lognormal distribution with mean=par1 and std=par2, returns only integer values
+    - 'poisson'       : Poisson distribution with rate=par1 (par2 is not used); mean and variance are equal to par1
+    - 'neg_binomial'  : negative binomial distribution with mean=par1 and k=par2; converges to Poisson with k=∞
+    - 'beta'          : beta distribution with alpha=par1 and beta=par2;
+    - 'gamma'         : gamma distribution with shape=par1 and scale=par2;
+
+    Args:
+        dist (str):   the distribution to sample from
+        par1 (float): the "main" distribution parameter (e.g. mean)
+        par2 (float): the "secondary" distribution parameter (e.g. std)
+        size (int):   the number of samples (default=1)
+        kwargs (dict): passed to individual sampling functions
+
+    Returns:
+        A length N array of samples
+
+    **Examples**::
+
+        ss.sample() # returns Unif(0,1)
+        ss.sample(dist='normal', par1=3, par2=0.5) # returns Normal(μ=3, σ=0.5)
+        ss.sample(dist='lognormal_int', par1=5, par2=3) # returns lognormally distributed values with mean 5 and std 3
+
+    Notes:
+        Lognormal distributions are parameterized with reference to the underlying normal distribution (see:
+        https://docs.scipy.org/doc/numpy-1.14.0/reference/generated/numpy.random.lognormal.html), but this
+        function assumes the user wants to specify the mean and std of the lognormal distribution.
+
+        Negative binomial distributions are parameterized with reference to the mean and dispersion parameter k
+        (see: https://en.wikipedia.org/wiki/Negative_binomial_distribution). The r parameter of the underlying
+        distribution is then calculated from the desired mean and k. For a small mean (~1), a dispersion parameter
+        of ∞ corresponds to the variance and standard deviation being equal to the mean (i.e., Poisson). For a
+        large mean (e.g. >100), a dispersion parameter of 1 corresponds to the standard deviation being equal to
+        the mean.
+    """
+
+    # Some of these have aliases, but these are the "official" names
+    choices = [
+        'uniform',
+        'normal',
+        'choice',
+        'normal_pos',
+        'normal_int',
+        'lognormal',
+        'lognormal_int',
+        'poisson',
+        'neg_binomial',
+        'beta',
+        'gamma',
+    ]
+
+    # Ensure it's an integer
+    if size is not None and not isinstance(size, tuple):
+        size = int(size)
+
+    # Compute distribution parameters and draw samples
+    # NB, if adding a new distribution, also add to choices above
+    if dist in ['unif', 'uniform']:
+        samples = np.random.uniform(low=par1, high=par2, size=size)
+    elif dist in ['choice']:
+        samples = np.random.choice(a=par1, p=par2, size=size, **kwargs)
+    elif dist in ['norm', 'normal']:
+        samples = np.random.normal(loc=par1, scale=par2, size=size)
+    elif dist == 'normal_pos':
+        samples = np.abs(np.random.normal(loc=par1, scale=par2, size=size))
+    elif dist == 'normal_int':
+        samples = np.round(np.abs(np.random.normal(loc=par1, scale=par2, size=size)))
+    elif dist == 'poisson':
+        samples = n_poisson(rate=par1, n=size)  # Use Numba version below for speed
+    elif dist == 'neg_binomial':
+        samples = n_neg_binomial(rate=par1, dispersion=par2, n=size, **kwargs)  # Use custom version below
+    elif dist == 'beta':
+        samples = np.random.beta(a=par1, b=par2, size=size)
+    elif dist == 'gamma':
+        samples = np.random.gamma(shape=par1, scale=par2, size=size)
+    elif dist in ['lognorm', 'lognormal', 'lognorm_int', 'lognormal_int']:
+        if (sc.isnumber(par1) and par1 > 0) or (sc.checktype(par1, 'arraylike') and (par1 > 0).all()):
+            mean = np.log(
+                par1 ** 2 / np.sqrt(par2 ** 2 + par1 ** 2))  # Computes the mean of the underlying normal distribution
+            sigma = np.sqrt(np.log(par2 ** 2 / par1 ** 2 + 1))  # Computes sigma for the underlying normal distribution
+            samples = np.random.lognormal(mean=mean, sigma=sigma, size=size)
+        else:
+            samples = np.zeros(size)
+        if '_int' in dist:
+            samples = np.round(samples)
+    # Calculate a and b using mean (par1) and variance (par2)
+    # https://stats.stackexchange.com/questions/12232/calculating-the-parameters-of-a-beta-distribution-using-the-mean-and-variance
+    elif dist == 'beta_mean':
+        a = ((1 - par1) / par2 - 1 / par1) * par1 ** 2
+        b = a * (1 / par1 - 1)
+        samples = np.random.beta(a=a, b=b, size=size)
+    else:
+        errormsg = f'The selected distribution "{dist}" is not implemented; choices are: {sc.newlinejoin(choices)}'
+        raise NotImplementedError(errormsg)
+
+    return samples
+
+
+def get_pdf(dist=None, par1=None, par2=None):
+    """
+    Return a probability density function for the specified distribution. This
+    is used for example by test_num to retrieve the distribution of times from
+    symptom-to-swab for testing. For example, for Washington State, these values
+    are dist='lognormal', par1=10, par2=170.
+    """
+    import scipy.stats as sps  # Import here since slow
+
+    choices = [
+        'none',
+        'uniform',
+        'lognormal',
+    ]
+
+    if dist in ['None', 'none', None]:
+        return None
+    elif dist == 'uniform':
+        pdf = sps.uniform(loc=par1, scale=par2)
+    elif dist == 'lognormal':
+        mean = np.log(par1 ** 2 / np.sqrt(par2 + par1 ** 2))  # Computes the mean of the underlying normal distribution
+        sigma = np.sqrt(np.log(par2 / par1 ** 2 + 1))  # Computes sigma for the underlying normal distribution
+        pdf = sps.lognorm(sigma, loc=-0.5, scale=np.exp(mean))
+    else:
+        choicestr = '\n'.join(choices)
+        errormsg = f'The selected distribution "{dist}" is not implemented; choices are: {choicestr}'
+        raise NotImplementedError(errormsg)
+
+    return pdf
+
+
+def set_seed(seed=None):
+    """
+    Reset the random seed. This function also resets Python's built-in random
+    number generated.
+
+    Args:
+        seed (int): the random seed
+    """
+    # Dies if a float is given
+    if seed is not None:
+        seed = int(seed)
+    np.random.seed(seed)  # If None, reinitializes it
+    return
+
+
+# %% Probabilities -- mostly not jitted since performance gain is minimal
+
+__all__ += ['n_binomial', 'binomial_filter', 'binomial_arr', 'n_multinomial',
+            'poisson', 'n_poisson', 'n_neg_binomial', 'choose', 'choose_r', 'choose_w']
+
+
+def n_binomial(prob, n):
+    """
+    Perform multiple binomial (Bernolli) trials
+
+    Args:
+        prob (float): probability of each trial succeeding
+        n (int): number of trials (size of array)
+
+    Returns:
+        Boolean array of which trials succeeded
+
+    **Example**::
+
+        outcomes = ss.n_binomial(0.5, 100) # Perform 100 coin-flips
+    """
+    return np.random.random(n) < prob
+
+
+def binomial_filter(prob, arr):
+    """
+    Binomial "filter" -- the same as n_binomial, except return
+    the elements of arr that succeeded.
+
+    Args:
+        prob (float): probability of each trial succeeding
+        arr (array): the array to be filtered
+
+    Returns:
+        Subset of array for which trials succeeded
+
+    **Example**::
+
+        inds = ss.binomial_filter(0.5, np.arange(20)**2) # Which values in the (arbitrary) array passed the coin flip
+    """
+    return arr[(np.random.random(len(arr)) < prob).nonzero()[0]]
+
+
+def binomial_arr(prob_arr):
+    """
+    Binomial (Bernoulli) trials each with different probabilities.
+
+    Args:
+        prob_arr (array): array of probabilities
+
+    Returns:
+         Boolean array of which trials on the input array succeeded
+
+    **Example**::
+
+        outcomes = ss.binomial_arr([0.1, 0.1, 0.2, 0.2, 0.8, 0.8]) # Perform 6 trials with different probabilities
+    """
+    return np.random.random(prob_arr.shape) < prob_arr
+
+
+def n_multinomial(probs, n):  # No speed gain from Numba
+    """
+    An array of multinomial trials.
+
+    Args:
+        probs (array): probability of each outcome, which usually should sum to 1
+        n (int): number of trials
+
+    Returns:
+        Array of integer outcomes
+
+    **Example**::
+
+        outcomes = ss.n_multinomial(np.ones(6)/6.0, 50)+1 # Return 50 die-rolls
+    """
+    return np.searchsorted(np.cumsum(probs), np.random.random(n))
+
+
+def poisson(rate):
+    """
+    A Poisson trial.
+
+    Args:
+        rate (float): the rate of the Poisson process
+
+    **Example**::
+
+        outcome = ss.poisson(100) # Single Poisson trial with mean 100
+    """
+    return np.random.poisson(rate, 1)[0]
+
+
+def n_poisson(rate, n):
+    """
+    An array of Poisson trials.
+
+    Args:
+        rate (float): the rate of the Poisson process (mean)
+        n (int): number of trials
+
+    **Example**::
+
+        outcomes = ss.n_poisson(100, 20) # 20 Poisson trials with mean 100
+    """
+    return np.random.poisson(rate, n)
+
+
+def n_neg_binomial(rate, dispersion, n, step=1):  # Numba not used due to incompatible implementation
+    """
+    An array of negative binomial trials. See ss.sample() for more explanation.
+
+    Args:
+        rate (float): the rate of the process (mean, same as Poisson)
+        dispersion (float):  dispersion parameter; lower is more dispersion, i.e. 0 = infinite, ∞ = Poisson
+        n (int): number of trials
+        step (float): the step size to use if non-integer outputs are desired
+
+    **Example**::
+
+        outcomes = ss.n_neg_binomial(100, 1, 50) # 50 negative binomial trials with mean 100 and dispersion roughly equal to mean (large-mean limit)
+        outcomes = ss.n_neg_binomial(1, 100, 20) # 20 negative binomial trials with mean 1 and dispersion still roughly equal to mean (approximately Poisson)
+    """
+    nbn_n = dispersion
+    nbn_p = dispersion / (rate / step + dispersion)
+    samples = np.random.negative_binomial(n=nbn_n, p=nbn_p, size=n) * step
+    return samples
+
+
+def choose(max_n, n):
+    """
+    Choose a subset of items (e.g., people) without replacement.
+
+    Args:
+        max_n (int): the total number of items
+        n (int): the number of items to choose
+
+    **Example**::
+
+        choices = ss.choose(5, 2) # choose 2 out of 5 people with equal probability (without repeats)
+    """
+    return np.random.choice(max_n, n, replace=False)
+
+
+def choose_r(max_n, n):
+    """
+    Choose a subset of items (e.g., people), with replacement.
+
+    Args:
+        max_n (int): the total number of items
+        n (int): the number of items to choose
+
+    **Example**::
+
+        choices = ss.choose_r(5, 10) # choose 10 out of 5 people with equal probability (with repeats)
+    """
+    return np.random.choice(max_n, n, replace=True)
+
+
+def choose_w(probs, n, unique=True):  # No performance gain from Numba
+    """
+    Choose n items (e.g. people), each with a probability from the distribution probs.
+
+    Args:
+        probs (array): list of probabilities, should sum to 1
+        n (int): number of samples to choose
+        unique (bool): whether or not to ensure unique indices
+
+    **Example**::
+
+        choices = ss.choose_w([0.2, 0.5, 0.1, 0.1, 0.1], 2) # choose 2 out of 5 people with nonequal probability.
+    """
+    probs = np.array(probs)
+    n_choices = len(probs)
+    n_samples = int(n)
+    probs_sum = probs.sum()
+    if probs_sum:  # Weight is nonzero, rescale
+        probs = probs / probs_sum
+    else:  # Weights are all zero, choose uniformly
+        probs = np.ones(n_choices) / n_choices
+    return np.random.choice(n_choices, n_samples, p=probs, replace=not (unique))
+
+
+# %% Simple array operations
+
+__all__ += ['true', 'false', 'defined', 'undefined',
+            'itrue', 'ifalse', 'idefined', 'iundefined',
+            'itruei', 'ifalsei', 'idefinedi', 'iundefinedi',
+            'dtround', 'find_cutoff']
+
+
+def true(arr):
+    """
+    Returns the indices of the values of the array that are true: just an alias
+    for arr.nonzero()[0].
+
+    Args:
+        arr (array): any array
+
+    **Example**::
+
+        inds = ss.true(np.array([1,0,0,1,1,0,1])) # Returns array([0, 3, 4, 6])
+    """
+    return arr.nonzero()[-1]
+
+
+def false(arr):
+    """
+    Returns the indices of the values of the array that are false.
+
+    Args:
+        arr (array): any array
+
+    **Example**::
+
+        inds = ss.false(np.array([1,0,0,1,1,0,1]))
+    """
+    return np.logical_not(arr).nonzero()[-1]
+
+
+def defined(arr):
+    """
+    Returns the indices of the values of the array that are not-nan.
+
+    Args:
+        arr (array): any array
+
+    **Example**::
+
+        inds = ss.defined(np.array([1,np.nan,0,np.nan,1,0,1]))
+    """
+    return (~np.isnan(arr)).nonzero()[-1]
+
+
+def undefined(arr):
+    """
+    Returns the indices of the values of the array that are not-nan.
+
+    Args:
+        arr (array): any array
+
+    **Example**::
+
+        inds = ss.defined(np.array([1,np.nan,0,np.nan,1,0,1]))
+    """
+    return np.isnan(arr).nonzero()[-1]
+
+
+def itrue(arr, inds):
+    """
+    Returns the indices that are true in the array -- name is short for indices[true]
+
+    Args:
+        arr (array): a Boolean array, used as a filter
+        inds (array): any other array (usually, an array of indices) of the same size
+
+    **Example**::
+
+        inds = ss.itrue(np.array([True,False,True,True]), inds=np.array([5,22,47,93]))
+    """
+    return inds[arr]
+
+
+def ifalse(arr, inds):
+    """
+    Returns the indices that are true in the array -- name is short for indices[false]
+
+    Args:
+        arr (array): a Boolean array, used as a filter
+        inds (array): any other array (usually, an array of indices) of the same size
+
+    **Example**::
+
+        inds = ss.ifalse(np.array([True,False,True,True]), inds=np.array([5,22,47,93]))
+    """
+    return inds[np.logical_not(arr)]
+
+
+def idefined(arr, inds):
+    """
+    Returns the indices that are defined in the array -- name is short for indices[defined]
+
+    Args:
+        arr (array): any array, used as a filter
+        inds (array): any other array (usually, an array of indices) of the same size
+
+    **Example**::
+
+        inds = ss.idefined(np.array([3,np.nan,np.nan,4]), inds=np.array([5,22,47,93]))
+    """
+    return inds[~np.isnan(arr)]
+
+
+def iundefined(arr, inds):
+    """
+    Returns the indices that are undefined in the array -- name is short for indices[undefined]
+
+    Args:
+        arr (array): any array, used as a filter
+        inds (array): any other array (usually, an array of indices) of the same size
+
+    **Example**::
+
+        inds = ss.iundefined(np.array([3,np.nan,np.nan,4]), inds=np.array([5,22,47,93]))
+    """
+    return inds[np.isnan(arr)]
+
+
+def itruei(arr, inds):
+    """
+    Returns the indices that are true in the array -- name is short for indices[true[indices]]
+
+    Args:
+        arr (array): a Boolean array, used as a filter
+        inds (array): an array of indices for the original array
+
+    **Example**::
+
+        inds = ss.itruei(np.array([True,False,True,True,False,False,True,False]), inds=np.array([0,1,3,5]))
+    """
+    return inds[arr[inds]]
+
+
+def ifalsei(arr, inds):
+    """
+    Returns the indices that are false in the array -- name is short for indices[false[indices]]
+
+    Args:
+        arr (array): a Boolean array, used as a filter
+        inds (array): an array of indices for the original array
+
+    **Example**::
+
+        inds = ss.ifalsei(np.array([True,False,True,True,False,False,True,False]), inds=np.array([0,1,3,5]))
+    """
+    return inds[np.logical_not(arr[inds])]
+
+
+def idefinedi(arr, inds):
+    """
+    Returns the indices that are defined in the array -- name is short for indices[defined[indices]]
+
+    Args:
+        arr (array): any array, used as a filter
+        inds (array): an array of indices for the original array
+
+    **Example**::
+
+        inds = ss.idefinedi(np.array([4,np.nan,0,np.nan,np.nan,4,7,4,np.nan]), inds=np.array([0,1,3,5]))
+    """
+    return inds[~np.isnan(arr[inds])]
+
+
+def iundefinedi(arr, inds):
+    """
+    Returns the indices that are undefined in the array -- name is short for indices[defined[indices]]
+
+    Args:
+        arr (array): any array, used as a filter
+        inds (array): an array of indices for the original array
+
+    **Example**::
+
+        inds = ss.iundefinedi(np.array([4,np.nan,0,np.nan,np.nan,4,7,4,np.nan]), inds=np.array([0,1,3,5]))
+    """
+    return inds[np.isnan(arr[inds])]
+
+
+def dtround(arr, dt, ceil=True):
+    """
+    Rounds the values in the array to the nearest timestep
+
+    Args:
+        arr (array): any array
+        dt  (float): float, usually representing a timestep in years
+        ceil  (bool): whether to always round up
+
+    **Example**::
+
+        dtround = ss.dtround(np.array([0.23,0.61,20.53])) # Returns array([0.2, 0.6, 20.6])
+        dtround = ss.dtround(np.array([0.23,0.61,20.53]),ceil=True) # Returns array([0.4, 0.8, 20.6])
+    """
+    if ceil:
+        return np.ceil(arr * (1 / dt)) / (1 / dt)
+    else:
+        return np.round(arr * (1 / dt)) / (1 / dt)
+
+
+def find_cutoff(duration_cutoffs, duration):
+    """
+    Find which duration bin each ind belongs to.
+    """
+    return np.nonzero(duration_cutoffs <= duration)[0][-1]  # Index of the duration bin to use
+
+
+def logf1(x, k, ttc=25):
+    """
+    Logistic function passing through (0,0) and (ttc,1).
+    Accepts 1 parameter which determines the growth rate.
+    """
+    return logf3(x, k, 0, 1, ttc=ttc)
+
+
+def get_asymptotes(k, x_infl, s, ttc=25):
+    """
+    Get upper asymptotes for logistic functions
+    """
+    term1 = (1 + np.exp(k * (x_infl - ttc))) ** s  # Note, this is 1 for most parameter combinations
+    term2 = (1 + np.exp(k * x_infl)) ** s
+    u_asymp_num = term1 * (1 - term2)
+    u_asymp_denom = term1 - term2
+    u_asymp = u_asymp_num / u_asymp_denom
+    l_asymp = term1 / (term1 - term2)
+    return l_asymp, u_asymp
+
+
+def logf3(x, k, x_infl, s, ttc=25):
+    """
+    Logistic function passing through (0,0) and (ttc,1).
+    This version is derived from the 5-parameter version here: https://www.r-bloggers.com/2019/11/five-parameters-logistic-regression/
+    However, since it's constrained to pass through 2 points, there are 3 free parameters remaining.
+    Args:
+         k: growth rate, equivalent to b in https://www.r-bloggers.com/2019/11/five-parameters-logistic-regression/
+         x_infl: a location parameter, equivalent to C in https://www.r-bloggers.com/2019/11/five-parameters-logistic-regression/
+         s: asymmetry parameter, equivalent to s in https://www.r-bloggers.com/2019/11/five-parameters-logistic-regression/
+         ttc (time to cancer): x value for which the curve passess through 1. For x values beyond this, the function returns 1
+    """
+    l_asymp, u_asymp = get_asymptotes(k, x_infl, s, ttc)
+    return np.minimum(1, l_asymp + (u_asymp - l_asymp) / (1 + np.exp(k * (x_infl - x))) ** s)
+
+
+def logf2(x, k, x_infl, ttc=25):
+    """
+    Logistic function constrained to pass through (0,0) and (ttc,1).
+    This version is derived from the 5-parameter version here: https://www.r-bloggers.com/2019/11/five-parameters-logistic-regression/
+    Since it's constrained to pass through 2 points, there are 3 free parameters remaining, and this verison fixes s=1
+    Args:
+         k: growth rate, equivalent to b in https://www.r-bloggers.com/2019/11/five-parameters-logistic-regression/
+         x_infl: point of inflection, equivalent to C in https://www.r-bloggers.com/2019/11/five-parameters-logistic-regression/
+         ttc (time to cancer): x value for which the curve passess through 1. For x values beyond this, the function returns 1
+    """
+    return logf3(x, k, x_infl, s=1, ttc=ttc)
+
+
+def invlogf3(y, k, x_infl, s, ttc=25):
+    """
+    Inverse of logf3; see definition there for arguments
+    """
+    l_asymp, u_asymp = get_asymptotes(k, x_infl, s, ttc)
+    part1 = np.log((u_asymp - l_asymp) / (y - l_asymp)) / s
+    part2 = np.log(np.exp(part1) - 1)
+    final = 1 / k * (k * x_infl - part2)
+    return final
+
+
+def invlogf2(y, k, x_infl, ttc=25):
+    """
+    Inverse of logf2; see definition there for arguments
+    """
+    return invlogf3(y, k, x_infl, 1, ttc=ttc)
+
+
+def invlogf1(y, k, ttc=25):
+    """
+    The inverse of the concave part of a logistic function, with point of inflexion at 0,0
+    and upper asymptote at 1. Accepts 1 parameter which determines the growth rate.
+    """
+    return invlogf3(y, k, 0, 1, ttc=ttc)
+
+
+def indef_int_logf2(x, k, x_infl, ttc=25):
+    """
+    Indefinite integral of logf2; see definition there for arguments
+    """
+    num = np.exp(-x_infl * k) * (np.exp(k * ttc) + np.exp(x_infl * k)) * (
+                (np.exp(x_infl * k) + 1) * np.log(np.exp(k * x) + np.exp(x_infl * k)) - k * x)
+    denom = k * (np.exp(k * ttc) - 1)
+    return num / denom
+
+
+def intlogf2(upper, k, x_infl, ttc=25):
+    """
+    Integral of logf2 between 0 and the limit given by upper
+    """
+    # Find the upper limits not including the part past time to cancer
+    exceeding_ttc_inds = (upper > ttc).nonzero()
+    lims_to_find = np.minimum(ttc, upper)
+
+    # Take the integral
+    val_at_0 = indef_int_logf2(0, k, x_infl, ttc)
+    val_at_lim = indef_int_logf2(lims_to_find, k, x_infl, ttc)
+    integral = val_at_lim - val_at_0
+
+    # Deal with those whose duration of infection exceeds the time to cancer
+    # Note, another option would be to set their transformation probability to 1
+    excess_integral = upper[exceeding_ttc_inds] - ttc
+    integral[exceeding_ttc_inds] += excess_integral
+
+    return integral
+
+
+def indef_int_logf1(x, k, ttc=25):
+    """
+    Indefinite integral of logf1; see definition there for arguments
+    """
+    return indef_int_logf2(x, k, 0, ttc=ttc)
+
+
+def intlogf1(upper, k, ttc=25):
+    """
+    Integral of logf1 between 0 and the limit given by upper
+    """
+    return intlogf2(upper, k, 0, ttc=ttc)
+
+
+def transform_prob(tp, dysp):
+    """
+    Returns transformation probability given dysplasia
+    Using formula for half an ellipsoid:
+        V = 1/2 * 4/3 * pi * a*b*c
+          = 2 * a*b*c
+          = 2* dysp * (dysp/2)**2, assuming that b = c = 1/2 a
+          = 1/2 * dysp**3
+    """
+    # return 1-np.power(1-tp, ((dysp*100)**2))
+    return 1 - np.power(1 - tp, 0.5 * ((dysp * 100) ** 3))
diff --git a/stisim/version.py b/stisim/version.py
index a5c600d6..060ba3f1 100644
--- a/stisim/version.py
+++ b/stisim/version.py
@@ -4,6 +4,6 @@
 
 __all__ = ['__version__', '__versiondate__', '__license__']
 
-__version__ = '0.0.1'
-__versiondate__ = '2023-06-14'
+__version__ = '0.0.2'
+__versiondate__ = '2023-06-29'
 __license__ = f'STIsim {__version__} ({__versiondate__}) — © 2023 by IDM'


Test patch:
diff --git a/tests/test_base.py b/tests/test_base.py
new file mode 100644
index 00000000..e0a87f8c
--- /dev/null
+++ b/tests/test_base.py
@@ -0,0 +1,133 @@
+"""
+Test objects from base.py
+"""
+
+# %% Imports and settings
+import sciris as sc
+import numpy as np
+import stisim as ss
+import pytest
+from stisim import utils as ssu
+import matplotlib.pyplot as plt
+
+
+# %% Define the tests
+
+def test_parsobj():
+    sc.heading('Testing parameters object')
+
+    pars1 = {'a': 1, 'b': 2}
+    parsobj = ss.ParsObj(pars1)
+
+    # Once created, you cannot directly add new keys to a parsobj, and a nonexistent key works like a dict
+    with pytest.raises(KeyError): parsobj['c'] = 3
+    with pytest.raises(KeyError): parsobj['c']
+
+    # Only a dict is allowed
+    with pytest.raises(TypeError):
+        pars2 = ['a', 'b']
+        ss.ParsObj(pars2)
+
+    return parsobj
+
+
+def test_people():
+    sc.heading('Testing base people object')
+
+    # Base people contains only the states defined in base.base_states
+    ppl = ss.BasePeople(100)  # BasePeople
+    del ppl
+
+    # Possible to initialize people with extra states, e.g. a geolocation
+    extra_states = ss.named_dict(
+        ss.StochState('geolocation', int, distdict=dict(dist='choice', par1=[1, 2, 3])),
+    )
+    ppl = ss.People(100, states=extra_states)
+
+    # Possible to add a module to people outside a sim (not typical workflow)
+    ppl.add_module(ss.HIV())
+
+    return ppl
+
+
+def test_networks():
+
+    # Make completely abstract layers
+    n_edges = 10_000
+    n_people = 1000
+    p1 = np.random.randint(n_people, size=n_edges)
+    p2 = np.random.randint(n_people, size=n_edges)
+    beta = np.ones(n_edges)
+    nw1 = ss.Network(p1=p1, p2=p2, beta=beta, label='rand')
+    nw2 = ss.Network(dict(p1=p1, p2=p2, beta=beta), label='rand')  # Alternate method
+
+    # Make people, then make a dynamic sexual layer and update it
+    ppl = ss.People(100)  # BasePeople
+    ppl.initialize()  # This seems to be necessary, although not completely clear why...
+    nw3 = ss.hpv_network()
+    nw3.initialize(ppl)
+    nw3.update(ppl, ti=1, dt=1)  # Update by providing a timestep & current time index
+
+    nw4 = ss.maternal()
+    nw4.initialize(ppl)
+    nw4.add_pairs(mother_inds=[1, 2, 3], unborn_inds=[100, 101, 102], dur=[1, 1, 1])
+
+    return nw1, nw2, nw3, nw4
+
+
+def test_microsim():
+    sc.heading('Test making people and providing them to a sim')
+
+    ppl = ss.People(100)
+    ppl.networks = ssu.named_dict(ss.simple_sexual(), ss.maternal())
+
+    # Make HIV module
+    hiv = ss.HIV()
+    # Set beta. The first entry represents transmission risk from infected p1 -> susceptible p2
+    # Need to be careful to get the ordering right. The set-up here assumes that in the simple
+    # sexual  network, p1 is male and p2 is female. In the maternal network, p1=mothers, p2=babies.
+    hiv.pars['beta'] = {'simple_sexual': [0.0008, 0.0004], 'maternal': [0.2, 0]}
+
+    sim = ss.Sim(people=ppl, modules=[hiv, ss.Pregnancy()])
+    sim.initialize()
+    sim.run()
+
+    plt.figure()
+    plt.plot(sim.tivec, sim.results.hiv.n_infected)
+    plt.title('HIV number of infections')
+    plt.show()
+
+    return sim
+
+def test_ppl_construction():
+
+    sim_pars = {'networks': [ss.simple_sexual()], 'n_agents': 100}
+    gon_pars = {'beta': {'simple_sexual': [0.08, 0.04]}}
+    gon = ss.Gonorrhea(pars=gon_pars)
+
+    sim = ss.Sim(pars=sim_pars, modules=[gon])
+    sim.initialize()
+    sim.run()
+    plt.figure()
+    plt.plot(sim.tivec, sim.results.gonorrhea.n_infected)
+    plt.title('Number of gonorrhea infections')
+    plt.show()
+
+    return sim
+
+
+
+# %% Run as a script
+if __name__ == '__main__':
+    # Start timing
+    T = sc.tic()
+
+    # Run tests
+    parsobj = test_parsobj()
+    ppl = test_people()
+    nw1, nw2, nw3, nw4 = test_networks()
+    sim1 = test_microsim()
+    sim2 = test_ppl_construction()
+
+    sc.toc(T)
+    print('Done.')
diff --git a/tests/test_prototype.py b/tests/test_prototype.py
deleted file mode 100644
index cffa5b9c..00000000
--- a/tests/test_prototype.py
+++ /dev/null
@@ -1,87 +0,0 @@
-# TOY MODEL
-# This toy model has HIV and gonorrhea with coinfection. The diseases interact
-# such that gonorrhea susceptibility depends on the HIV CD4 value, and the HIV
-# module includes ART as a state which decreases the impact of HIV coinfection by
-# raising the CD4 value.
-
-# Some challenges
-# - How would you make it so that people infected with gonorrhea reduce sexual activity
-#   until their infection is cleared?
-
-import numpy as np
-import sciris as sc
-from collections import defaultdict
-import matplotlib.pyplot as plt
-import stisim as ss
-
-def make_people(n):
-    people = ss.People(n)
-    people.contacts['random'] = ss.RandomDynamicSexualLayer(people)
-    people.contacts['msm'] = ss.StaticLayer(people, 'mm')
-
-    # relationship types, homosexual/heterosexual, births?
-    # do we model pregnancies? track parents? households?
-    return people
-
-
-
-#######
-
-# TODO - should the module be stateful or stateless?
-
-#### GONORRHEA SIMULATION
-
-# people = make_people(100)
-# pars = defaultdict(dict)
-# pars['gonorrhea']['beta'] = {'random':0.3,'msm':0.5}
-# sim = ss.Sim(people, [ss.Gonorrhea], pars=pars)
-# sim.run()
-# plt.figure()
-# plt.plot(sim.tvec, sim.results.gonorrhea.n_infected)
-
-#### HIV SIMULATION
-
-people = make_people(100)
-pars = defaultdict(dict)
-pars['hiv']['beta'] = {'random':0.3,'msm':0.5}
-sim = ss.Sim(people, [ss.HIV], pars=pars, analyzers=ss.CD4_analyzer())
-sim.run()
-plt.figure()
-plt.plot(sim.tvec, sim.results.hiv.n_infected)
-plt.title('HIV number of infections')
-plt.figure()
-plt.plot(sim.tvec,sim.analyzers[0].cd4)
-plt.title('CD4 counts')
-
-people = make_people(100)
-sim = ss.Sim(people, [ss.HIV], pars=pars, interventions=ss.ART([10,20],[20,40]),analyzers=ss.CD4_analyzer())
-sim.run()
-plt.figure()
-plt.plot(sim.tvec, sim.analyzers[0].cd4)
-plt.title('CD4 counts (ART)')
-
-print('Done.')
-#
-#
-# # # Custom module by user
-# class Gonorrhea_DR(cs.Gonorrhea):
-#     default_pars = sc.dcp(cs.Gonorrhea.default_pars)
-#     default_pars['p_death'] = 0.3
-
-    # We need to make it so that infection with Gonorrhea and Gonorrhea_DR is mutually exclusive
-
-
-#
-# How should coinfection with these two modules be prevented?
-#
-#
-# people = make_people(100)
-# pars = defaultdict(dict)
-# pars['gonorrhea']['beta'] = {'random':0.1,'msm':0.15}
-# pars['gonorrhea_dr']['beta'] = {'random':0.05,'msm':0.1}
-# sim = Sim(people, [Gonorrhea, Gonorrhea_DR], pars=pars)
-# sim.run()
-# plt.figure()
-# plt.plot(sim.tvec, sim.results.gonorrhea.prevalence)
-# plt.plot(sim.tvec, sim.results.gonorrhea_dr.prevalence, color='r')
-# plt.show()
\ No newline at end of file