Repository: myokit/myokit (language: python)

Issue:
Split `_aux` into two modules
## `myokit._aux`

- Methods that are loaded into the main `myokit` module
- For example model loading and saving

## `myokit.util`

- Methods that are used by myokit but not really myokit related
- Should contain the methods that are now hidden in aux, for example the float comparison ones
- Imported by default
- Included in the docs

If any classes/methods (e.g. myokit.Benchmarker) used to be public, but now move into utils, they should have a deprecated alias



Discussion / hints:

Fix patch:
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 976afd2a5..8c6bb7184 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -13,6 +13,10 @@ This page lists the main changes made to Myokit in each release.
   - [#682](https://github.com/MichaelClerx/myokit/pull/682) OpenCL code now recognises error code -1001.
   - [#683](https://github.com/MichaelClerx/myokit/pull/683) Now testing OpenCL code in CI.
   - [#704](https://github.com/MichaelClerx/myokit/pull/704) Added first example notebook to the github repository.
+  - [#728](https://github.com/MichaelClerx/myokit/pull/728) Added a context manager `myokit.tools.capture` that can redirect and capture output from Python `stdout` and `stderr` via a Python-only method or file-descriptor redirection and should be thread-safe.
+  - [#735](https://github.com/MichaelClerx/myokit/pull/735) The new module `myokit.float` contains several functions related to floating point numbers that were previously in the main `myokit` namespace or not part of the public API.
+  - [#735](https://github.com/MichaelClerx/myokit/pull/735) The new module `myokit.tools` contains several stand-alone tools that are used by Myokit and were previously in the main `myokit` namespace or not part of the public API.
+  - [#735](https://github.com/MichaelClerx/myokit/pull/735) The method `pid_hash`, which is used by compiled code to generate (hopefully) unique ids is now part of the public API.
 - Changed
   - [#581](https://github.com/MichaelClerx/myokit/pull/581) Powers are now rendered without spaces in mmt code, e.g. `x^2` instead of `x ^ 2`.
   - [#595](https://github.com/MichaelClerx/myokit/pull/595) The `Simulation` class now uses CVODES instead of CVODE as backend, which may require changes to your installation.
@@ -22,15 +26,19 @@ This page lists the main changes made to Myokit in each release.
   - [#595](https://github.com/MichaelClerx/myokit/pull/595) Fix to `Simulation` for Sundials 4.1.0.
   - [#689](https://github.com/MichaelClerx/myokit/pull/689) In Python 2, an `ImportError` is now raised if `myokit.ini` contains the sequence " ;" in any of its value (as this cannot be processed by Python 2's `ConfigParser`).
   - [#728](https://github.com/MichaelClerx/myokit/pull/728) The `LegacySimulation` class no longer suppresses CVODE warnings but passes them on to the `warnings` module (backported from new `Simulation` class).
-  - [#728](https://github.com/MichaelClerx/myokit/pull/728) The classes `PyCapture` and `SubCapture` have been replaced by a single context manager ``capture`` that can do Python-only or low-level redirection and should be thread-safe.
 - Deprecated
-  - [#595](https://github.com/MichaelClerx/myokit/pull/595) The classes `myokit.PSimulation` and `myokit.ICSimulation` are deprecated. Sensitivities can now be obtained much quicker with the single cell `Simulation` class.
+  - [#595](https://github.com/MichaelClerx/myokit/pull/595) The classes `myokit.PSimulation` and `myokit.ICSimulation` have been deprecated in favour of the new `Simulation` class.
+  - [#735](https://github.com/MichaelClerx/myokit/pull/735) The class `myokit.Benchmarker` is depreacted in favour of `myokit.tools.Benchmarker`.
+  - [#735](https://github.com/MichaelClerx/myokit/pull/735) The method `myokit.format_float_dict` is depreacted and will be removed in future versions.
+  - [#735](https://github.com/MichaelClerx/myokit/pull/735) The method `myokit.format_path` is depreacted in favour of `myokit.tools.format_path`.
+  - [#735](https://github.com/MichaelClerx/myokit/pull/735)The method `myokit.strfloat` is depreacted in favour of `myokit.float.str`.
 - Removed
   - [#683](https://github.com/MichaelClerx/myokit/pull/683) No longer testing on Python 2.7.6 on linux (still testing latest 2.7). No longer testing any Python 2.7 on Windows.
   - [#730](https://github.com/MichaelClerx/myokit/pull/730) Removed the method `myokit.pack_snapshot`.
+  - [#728](https://github.com/MichaelClerx/myokit/pull/728) The classes `myokit.PyCapture` and `myokit.SubCapture` have been replaced by a single context manager `myokit.tools.capture`.
 - Fixed
   - [#684](https://github.com/MichaelClerx/myokit/pull/684) Fixed OpenCL loading issue on OS/X (with special thanks to Martin Aguilar and David Augustin).
-  - [#686](https://github.com/MichaelClerx/myokit/pull/686) Fixed a (windows only) bug in `myokit.format_path()`.
+  - [#686](https://github.com/MichaelClerx/myokit/pull/686) Fixed a (windows only) bug in `myokit.tools.format_path()`.
   - [#689](https://github.com/MichaelClerx/myokit/pull/689) Path lists read from `myokit.ini` are now filtered for empty entries and closing semicolons.
 
 ## [1.32.0] - 2021-01-19
diff --git a/docs/source/api_core/Model.rst b/docs/source/api_core/Model.rst
index 89d3a90df..a09882664 100644
--- a/docs/source/api_core/Model.rst
+++ b/docs/source/api_core/Model.rst
@@ -23,6 +23,10 @@ Once created and validated, a model can be asked to
 which can then be used by an exporter to generate runnable simulation or
 analysis code.
 
+A :class:`ModelComparison` class is available that can compare models
+syntactically, while the :meth:`step` method can compare model output (either
+to other models' output or to stored reference files from external
+implementations).
 
 ``myokit.Model``
 ================
@@ -40,3 +44,9 @@ Name validation
 
 .. autofunction:: check_name
 
+Model comparison
+================
+
+.. autoclass:: ModelComparison
+
+.. autofunction:: step
diff --git a/docs/source/api_core/exceptions.rst b/docs/source/api_core/exceptions.rst
index 9d818deac..109766124 100644
--- a/docs/source/api_core/exceptions.rst
+++ b/docs/source/api_core/exceptions.rst
@@ -7,7 +7,7 @@ Exceptions
 .. currentmodule:: myokit
 
 Myokit tries to raise errors in a sensible manner. The following classes are
-used.
+used:
 
 Base classes
 ============
diff --git a/docs/source/api_core/float.rst b/docs/source/api_core/float.rst
new file mode 100644
index 000000000..1d80dd5fa
--- /dev/null
+++ b/docs/source/api_core/float.rst
@@ -0,0 +1,35 @@
+.. _api/float:
+
+**********************
+Floating point numbers
+**********************
+
+Myokit includes a module ``myokit.float`` that contains methods to work with
+floating point numbers.
+These are used for example when checking if a pacing event has started or ended
+*to within machine precision*, or when checking if two units are equivalent.
+This module is imported automatically when ``import myokit`` is run.
+
+.. currentmodule:: myokit.float
+
+Comparison
+==========
+
+.. autofunction:: eq
+
+.. autofunction:: geq
+
+.. autofunction:: close
+
+Rounding to integer
+===================
+
+.. autofunction:: cround
+
+.. autofunction:: round
+
+Conversion to string
+====================
+
+.. autofunction:: str
+
diff --git a/docs/source/api_core/index.rst b/docs/source/api_core/index.rst
index 9ec18cbbb..4eb343b34 100644
--- a/docs/source/api_core/index.rst
+++ b/docs/source/api_core/index.rst
@@ -45,6 +45,7 @@ section of the documentation.
     equations
     units
     io
-    misc
+    float
+    tools
     settings
     exceptions
diff --git a/docs/source/api_core/io.rst b/docs/source/api_core/io.rst
index 052dd700e..aa9ac5afe 100644
--- a/docs/source/api_core/io.rst
+++ b/docs/source/api_core/io.rst
@@ -1,8 +1,8 @@
 .. _api/io:
 
-*****************************
-Reading and writing mmt files
-*****************************
+*******************
+Reading and writing
+*******************
 
 .. currentmodule:: myokit
 
@@ -15,6 +15,7 @@ Myokit models can be loaded form disk using :func:`load()` and stored using
 ``(model, protocol, script)`` tuple. Support for partial loading and saving
 is provided through functions like :func:`load_model()` and :func:`split()`.
 
+Similar functions are available to load save model states.
 
 Basic commands
 ==============
@@ -25,6 +26,8 @@ Basic commands
 
 .. autofunction:: parse
 
+.. autofunction:: run
+
 Using parts of mmt files
 ========================
 
@@ -42,6 +45,16 @@ Using parts of mmt files
 
 .. autofunction:: save_script
 
+Default mmt file parts
+======================
+When generating new ``mmt`` files, default protocols and scripts can be
+inserted:
+
+.. autofunction:: default_protocol
+
+.. autofunction:: default_script
+
+
 Parsing parts mmt files
 ========================
 
@@ -68,6 +81,18 @@ should be of the type :class:`ParseError`.
 
 .. autofunction:: format_parse_error
 
+Loading and saving states
+=========================
+
+.. autofunction:: load_state
+
+.. autofunction:: load_state_bin
+
+.. autofunction:: save_state
+
+.. autofunction:: save_state_bin
+
+.. autofunction:: parse_state
 
 A note about the format
 =======================
diff --git a/docs/source/api_core/misc.rst b/docs/source/api_core/misc.rst
deleted file mode 100644
index bca2a0221..000000000
--- a/docs/source/api_core/misc.rst
+++ /dev/null
@@ -1,65 +0,0 @@
-.. _api/misc:
-
-*************
-Miscellaneous
-*************
-
-.. currentmodule:: myokit
-
-Model comparison
-================
-
-.. autoclass:: ModelComparison
-
-.. _api/modelcomparison:
-
-Benchmarking
-============
-
-.. autoclass:: Benchmarker
-
-.. _api/logging:
-
-State i/o
-=========
-
-.. autofunction:: load_state
-
-.. autofunction:: load_state_bin
-
-.. autofunction:: parse_state
-
-.. autofunction:: save_state
-
-.. autofunction:: save_state_bin
-
-String functions
-================
-
-.. autofunction:: date
-
-.. autofunction:: format_float_dict
-
-.. autofunction:: format_path
-
-.. autofunction:: strfloat
-
-.. autofunction:: time
-
-.. autofunction:: version
-
-Other functions
-===============
-
-.. autofunction:: default_protocol
-
-.. autofunction:: default_script
-
-.. autofunction:: numpy_writer
-
-.. autofunction:: python_writer
-
-.. autofunction:: run
-
-.. autofunction:: strip_expression_units
-
diff --git a/docs/source/api_core/settings.rst b/docs/source/api_core/settings.rst
index 5e0b58a0c..bbb38c9c5 100644
--- a/docs/source/api_core/settings.rst
+++ b/docs/source/api_core/settings.rst
@@ -4,11 +4,16 @@
 Settings
 ********
 
-A configuration file ``~/.config/myokit/myokit.ini`` (where ~ denotes the user
-home directory) can be used to configure Myokit. Specifically, this file can be
-used to indicate the location of Sundials and OpenCL shared libraries and
-header files.
+.. currentmodule:: myokit
+
+myokit.ini
+----------
 
+Myokit can be configured by editing the configuration file
+``~/.config/myokit/myokit.ini``, where ``~`` denotes the user home directory.
+This file is generated automatically when Myokit is installed, and can be
+modified to indicate e.g. the version and location of the Sundials (CVODES)
+library or the location of OpenCL libraries and header files.
 
 System information
 ------------------
@@ -16,6 +21,17 @@ System information
 The :meth:`system()` method can be used to check if the system meets the
 requirements to use Myokit and run simulations.
 
-.. currentmodule:: myokit
-
 .. autofunction:: system
+
+The :meth:`version()` method provides information about the current Myokit
+version.
+
+.. autofunction:: version
+
+Myokit uses standard formats for date and time. The current date and time can
+be obtained using :meth:`date()` and :meth:`time()`.
+
+.. autofunction:: date
+
+.. autofunction:: time
+
diff --git a/docs/source/api_core/tools.rst b/docs/source/api_core/tools.rst
new file mode 100644
index 000000000..efb26504b
--- /dev/null
+++ b/docs/source/api_core/tools.rst
@@ -0,0 +1,51 @@
+.. _api/tools:
+
+*****
+Tools
+*****
+
+Myokit contains a module ``myokit.tools`` with functions and classes that are
+used throughout Myokit, but are not particularly Myokit-specific.
+This module is imported automatically when ``import myokit`` is run.
+
+.. currentmodule:: myokit.tools
+
+Benchmarking
+============
+
+.. autoclass:: Benchmarker
+
+Capturing printed output
+========================
+
+.. autoclass:: capture
+
+File system
+===========
+
+.. autofunction:: format_path
+
+.. autofunction:: rmtree
+
+String comparison
+=================
+
+.. autofunction:: lvsd
+
+.. autofunction:: natural_sort_key
+
+Deprecated functions
+====================
+Some of the functions in ``tools`` have deprecated aliases, that will be
+removed in future Myokit releases.
+
+.. currentmodule:: myokit
+
+.. autoclass:: Benchmarker
+
+.. autofunction:: format_float_dict
+
+.. autofunction:: format_path
+
+.. autofunction:: strfloat
+
diff --git a/docs/source/api_core/units.rst b/docs/source/api_core/units.rst
index c452ccf50..52f817497 100644
--- a/docs/source/api_core/units.rst
+++ b/docs/source/api_core/units.rst
@@ -11,7 +11,7 @@ Model variables can declare their units using the syntax explained in
 :class:`Unit` objects.
 
 A large number of units can be accessed through ``myokit.units``, this module
-is imported as part of the main ``myokit`` package::
+is imported automatically when ``import myokit`` is run::
 
     >>> import myokit
     >>> print myokit.units.m
@@ -182,3 +182,8 @@ Unit arithmetic
 
 .. autoclass:: Quantity
 
+For people who don't like units
+===============================
+
+.. autofunction:: strip_expression_units
+
diff --git a/docs/source/api_formats/index.rst b/docs/source/api_formats/index.rst
index 816f92aa8..e0cc17565 100644
--- a/docs/source/api_formats/index.rst
+++ b/docs/source/api_formats/index.rst
@@ -124,3 +124,15 @@ obtain e.g. an :class:`Exporter` via
 
 .. autofunction:: register_external_ewriter
 
+Default expression writers
+==========================
+
+Finally, Myokit contains two default expression writers, which are used
+internally to write expressions in Python format with or without NumPy support.
+
+.. currentmodule:: myokit
+
+.. autofunction:: python_writer
+
+.. autofunction:: numpy_writer
+
diff --git a/docs/source/api_index/index.rst b/docs/source/api_index/index.rst
index fb1ed7e96..5b857db55 100644
--- a/docs/source/api_index/index.rst
+++ b/docs/source/api_index/index.rst
@@ -10,9 +10,11 @@ sorted by module
 ..  toctree::
 
     myokit
+    myokit.float
     myokit.formats
     myokit.gui
     myokit.lib
     myokit.pacing
     myokit.pype
+    myokit.tools
 
diff --git a/docs/source/api_index/myokit.float.rst b/docs/source/api_index/myokit.float.rst
new file mode 100644
index 000000000..08c1fae4e
--- /dev/null
+++ b/docs/source/api_index/myokit.float.rst
@@ -0,0 +1,13 @@
+.. _api/index/myokit/float:
+
+============
+myokit.float
+============
+
+- :meth:`myokit.float.close`
+- :meth:`myokit.float.cround`
+- :meth:`myokit.float.eq`
+- :meth:`myokit.float.geq`
+- :meth:`myokit.float.round`
+- :meth:`myokit.float.str`
+
diff --git a/docs/source/api_index/myokit.rst b/docs/source/api_index/myokit.rst
index bb9a62bfc..4a6966169 100644
--- a/docs/source/api_index/myokit.rst
+++ b/docs/source/api_index/myokit.rst
@@ -9,7 +9,6 @@ myokit
 - :class:`myokit.ASin`
 - :class:`myokit.ATan`
 - :class:`myokit.Benchmarker`
-- :class:`myokit.capture`
 - :class:`myokit.Ceil`
 - :meth:`myokit.check_name`
 - :class:`myokit.CModel`
@@ -112,6 +111,7 @@ myokit
 - :meth:`myokit.parse_state`
 - :meth:`myokit.parse_unit`
 - :class:`myokit.PartialDerivative`
+- :meth:`myokit.pid_hash`
 - :class:`myokit.Piecewise`
 - :class:`myokit.Plus`
 - :class:`myokit.Power`
diff --git a/docs/source/api_index/myokit.tools.rst b/docs/source/api_index/myokit.tools.rst
new file mode 100644
index 000000000..00d2c867f
--- /dev/null
+++ b/docs/source/api_index/myokit.tools.rst
@@ -0,0 +1,13 @@
+.. _api/index/myokit/tools:
+
+============
+myokit.tools
+============
+
+- :class:`myokit.tools.Benchmarker`
+- :class:`myokit.tools.capture`
+- :meth:`myokit.tools.format_path`
+- :meth:`myokit.tools.lvsd`
+- :meth:`myokit.tools.natural_sort_key`
+- :meth:`myokit.tools.rmtree`
+
diff --git a/docs/source/api_simulations/Protocol.rst b/docs/source/api_simulations/Protocol.rst
index dee14fda8..e00f700e1 100644
--- a/docs/source/api_simulations/Protocol.rst
+++ b/docs/source/api_simulations/Protocol.rst
@@ -19,7 +19,7 @@ Protocol factory
 
 The ``myokit.pacing`` module provides a factory methods to facilitate the
 creation of :class:`Protocol <myokit.Protocol>` objects directly from python.
-This module is imported as part of the main ``myokit`` package.
+This module is imported automatically when ``import myokit`` is run.
 
 Periodic pacing
 ---------------
diff --git a/docs/source/api_simulations/RhsBenchmarker.rst b/docs/source/api_simulations/RhsBenchmarker.rst
index 87d20497d..93300b6ec 100644
--- a/docs/source/api_simulations/RhsBenchmarker.rst
+++ b/docs/source/api_simulations/RhsBenchmarker.rst
@@ -1,12 +1,10 @@
 .. _api/simulations/myokit.RhsBenchmarker:
 
-***************************
-Right-hand-side evaluations
-***************************
+****************************************
+Benchmarking right-hand-side evaluations
+****************************************
 
 .. currentmodule:: myokit
 
-.. autofunction:: step
-
 .. autoclass:: RhsBenchmarker
 
diff --git a/docs/source/api_simulations/backend.rst b/docs/source/api_simulations/backend.rst
index 3b6695ba2..dfd269b41 100644
--- a/docs/source/api_simulations/backend.rst
+++ b/docs/source/api_simulations/backend.rst
@@ -19,12 +19,10 @@ Back-end classes
     :private-members:
     :inherited-members:
 
-.. autoclass:: capture
-    :private-members:
-    :inherited-members:
-
 .. autoclass:: Compiler
 
+.. autofunction:: pid_hash
+
 
 Templating engine
 -----------------
diff --git a/myokit/__init__.py b/myokit/__init__.py
index e8d0bf89e..b6b11fd07 100644
--- a/myokit/__init__.py
+++ b/myokit/__init__.py
@@ -374,6 +374,10 @@
 del(ex, name, clas, _globals, inspect)  # Prevent public visibility
 del(_err)
 
+# Tools
+from . import float  # noqa
+from . import tools  # noqa
+
 # Model api
 from ._model_api import ( # noqa
     check_name,
@@ -448,6 +452,7 @@
     Protocol,
     ProtocolEvent,
 )
+from . import pacing  # noqa
 
 # Parser functions
 from ._parsing import (  # noqa
@@ -456,7 +461,6 @@
     parse,
     parse_expression_string as parse_expression,
     parse_model,
-    #parse_number_string as parse_number,
     parse_protocol,
     parse_state,
     parse_unit_string as parse_unit,
@@ -464,69 +468,43 @@
     strip_expression_units,
 )
 
-# Auxillary functions
-from ._aux import (  # noqa
-    # Version info
-    version,
-
-    # Global date and time formats
-    date,
-    time,
-
-    # Default mmt parts
-    default_protocol,
-    default_script,
-
+# Load/save functions
+from ._io import (  # noqa
     # Reading, writing
     load,
     load_model,
     load_protocol,
     load_script,
+    load_state,
+    load_state_bin,
     save,
     save_model,
     save_protocol,
     save_script,
-    load_state,
     save_state,
-    load_state_bin,
     save_state_bin,
+)
 
-    # Running scripts
-    run,
-
-    # Test step
-    step,
-
-    # Output capturing
-    capture,
-
-    # Sorting
-    _natural_sort_key,
+# Common units
+from . import units  # noqa, also loads all common unit names
 
-    # Dyanmic generation of Python/NumPy expressions
-    python_writer,
-    numpy_writer,
-
-    # Model comparison
+# Auxillary functions
+from ._aux import (  # noqa
+    default_protocol,
+    default_script,
     ModelComparison,
-
-    # Benchmarking
+    numpy_writer,
+    python_writer,
+    run,
+    step,
+    version,
+    # Deprecated and/or moved to myokit.tools
     Benchmarker,
-
-    # Floating point
-    _close,
-    _cround,
-    _feq,
-    _fgeq,
-    _fround,
+    date,
     format_float_dict,
-    strfloat,
-
-    # Misc
     format_path,
-    _lvsd,
-    _pid_hash,
-    _rmtree,
+    strfloat,
+    time,
 )
 
 # System information
@@ -555,11 +533,11 @@
     DataBlock2d,
 )
 
-
 # Simulations
 from ._sim import (  # noqa
     CModule,
     CppModule,
+    pid_hash,
 )
 from ._sim.compiler import (  # noqa
     Compiler,
@@ -584,15 +562,6 @@
 from ._sim.openclsim import SimulationOpenCL                    # noqa
 from ._sim.fiber_tissue import FiberTissueSimulation            # noqa
 
-# Import whole modules
-# This allows these modules to be used after myokit was imported, without
-# importing the modules specifically (like os and os.path).
-# All modules imported here must report so in their docs
-from . import ( # noqa
-    pacing,
-    units,  # Also loads all common unit names
-)
-
 
 #
 # Globally shared progress reporter
diff --git a/myokit/__main__.py b/myokit/__main__.py
index 9048497f1..e6daf1730 100644
--- a/myokit/__main__.py
+++ b/myokit/__main__.py
@@ -338,7 +338,7 @@ def mmt_export(exporter, source, target):
 
     # Parse input file
     try:
-        print('Reading model from ' + myokit.format_path(source))
+        print('Reading model from ' + myokit.tools.format_path(source))
         model, protocol, script = myokit.load(source)
     except myokit.ParseError as ex:
         print(myokit.format_parse_error(ex, source))
@@ -609,7 +609,7 @@ def install_windows():
         print('Done')
 
     finally:
-        myokit._rmtree(tdir)
+        myokit.tools.rmtree(tdir)
 
 
 def add_icon_parser(subparsers):
@@ -979,7 +979,7 @@ def reset(force=False):
     if remove:
         print('Removing')
         print('  ' + myokit.DIR_USER)
-        myokit._rmtree(myokit.DIR_USER)
+        myokit.tools.rmtree(myokit.DIR_USER)
         print('Done')
     else:
         print('Aborting.')
@@ -1021,7 +1021,7 @@ def run(source, debug, debugfile):
     # Read mmt file
     try:
         print('Reading model from ' + source)
-        b = myokit.Benchmarker()
+        b = myokit.tools.Benchmarker()
         (model, protocol, script) = myokit.load(source)
         print('File loaded in ' + str(b.time()) + ' seconds')
         if model is None:
@@ -1163,7 +1163,7 @@ def step(source, ref, ini, raw):
     try:
         if raw:
             derivs = model.eval_state_derivatives(state=ini)
-            print('\n'.join([myokit.strfloat(x) for x in derivs]))
+            print('\n'.join([myokit.float.str(x) for x in derivs]))
         else:
             print(myokit.step(model, initial=ini, reference=ref))
     except myokit.NumericalError as ee:
@@ -1791,7 +1791,7 @@ def test_examples_single(root, path):
     import subprocess
     import sys
 
-    b = myokit.Benchmarker()
+    b = myokit.tools.Benchmarker()
     print('Running ' + path + ' ... ', end='')
     sys.stdout.flush()
 
diff --git a/myokit/_aux.py b/myokit/_aux.py
index ff96a5215..e5d459938 100644
--- a/myokit/_aux.py
+++ b/myokit/_aux.py
@@ -9,437 +9,37 @@
 from __future__ import absolute_import, division
 from __future__ import print_function, unicode_literals
 
-import array
-import os
-import re
-import shutil
-import stat
-import sys
-import tempfile
-import threading
-import timeit
-
-# StringIO in Python 2 and 3
-try:
-    from cStringIO import StringIO
-except ImportError:  # pragma: no python 2 cover
-    from io import StringIO
-
-# String types in Python 2 and 3
-try:
-    basestring
-except NameError:   # pragma: no python 2 cover
-    basestring = str
-
 import myokit
 
+import sys
+
 # Globally shared numpy expression writer
 _numpywriter_ = None
 
 # Globally shared python expression writer
 _pywriter_ = None
 
-# Natural sort regex
-_natural_sort_regex = re.compile('([0-9]+)')
+
+class Benchmarker(myokit.tools.Benchmarker):
+    """Deprecated alias of :class:`myokit.tools.Benchmarker`."""
+
+    def __init__(self):
+        # Deprecated since 2021-03-24
+        import warnings
+        warnings.warn(
+            'The class `myokit.Benchmarker` is deprecated.'
+            ' Please use `myokit.tools.Benchmarker` instead.')
+        super(Benchmarker, self).__init__()
 
 
 def date():
     """
-    Returns the current date and time, formatted as specified in
-    ``myokit.settings``.
+    Returns the current date and time, in the format used throughout Myokit.
     """
     import time
     return time.strftime(myokit.DATE_FORMAT)
 
 
-def time():
-    """
-    Returns the current time, formatted as specified in ``myokit.settings``.
-    """
-    import time as t
-    return t.strftime(myokit.TIME_FORMAT)
-
-
-def _natural_sort_key(s):
-    """
-    Natural sort key, from: http://stackoverflow.com/questions/4836710/
-
-    Usage example::
-
-        values.sort(key=lambda x: myokit._natural_sort_key(x))
-
-    """
-    return [
-        int(text) if text.isdigit() else text.lower()
-        for text in _natural_sort_regex.split(s)]
-
-
-class Benchmarker(object):
-    """
-    Allows benchmarking using the with statement.
-
-    Example::
-
-        m,p,x = myokit.load('example')
-        s = myokit.Simulation(m, p)
-        b = myokit.Benchmarker()
-        s.run()
-        print(b.time())
-        b.reset()
-        s.run()
-        print(b.time())
-
-    """
-    def __init__(self, output=None):
-        self._start = timeit.default_timer()
-
-    def format(self, time):
-        """
-        Formats a (non-integer) number of seconds, returns a string like
-        "5 weeks, 3 days, 1 hour, 4 minutes, 9 seconds", or "0.0019 seconds".
-        """
-        if time < 60:
-            return '1 second' if time == 1 else str(time) + ' seconds'
-        output = []
-        time = int(round(time))
-        units = [
-            (604800, 'week'),
-            (86400, 'day'),
-            (3600, 'hour'),
-            (60, 'minute'),
-        ]
-        for k, name in units:
-            f = time // k
-            if f > 0 or output:
-                output.append(str(f) + ' ' + (name if f == 1 else name + 's'))
-            time -= f * k
-        output.append('1 second' if time == 1 else str(time) + ' seconds')
-        return ', '.join(output)
-
-    def reset(self):
-        """
-        Resets this timer's start time.
-        """
-        self._start = timeit.default_timer()
-
-    def time(self):
-        """
-        Returns the time since benchmarking started.
-        """
-        return timeit.default_timer() - self._start
-
-
-class capture(object):
-    """
-    Context manager that temporarily redirects the current standard output and
-    error streams, and captures anything that's written to them.
-
-    Example::
-
-        with myokit.capture() as a:
-            print('This will be captured')
-
-    Within a single thread, captures can be nested, for example::
-
-        with myokit.capture() as a:
-            print('This will be captured by a')
-            with myokit.capture() as b:
-                print('This will be captured by b, not a')
-            print('This will be captured by a again')
-
-    Capturing is thread-safe: a lock is used to ensure only a single thread is
-    capturing at any time. For example, if we have a function::
-
-        def f(i):
-            with myokit.capture() as a:
-                print(i)
-                ...
-
-            return a.text()
-
-    and this is called from several threads, the ``capture`` acts as a lock (a
-    ``threading.RLock``) so that one thread will need to finish executing the
-    code within the ``with`` statement before a second thread can start
-    capturing.
-
-    In multiprocessing, no locks are used, and no memory or streams are shared
-    so that this should also be safe.
-
-    By default, this method works by simply redirecting ``sys.stdout`` and
-    ``sys.stderr``. This captures any output written by the Python interpreter,
-    but does not catch output from C/C++ extensions or subrocesses. To also
-    catch that output start the capture with the optional argument
-    ``fd=True``, which enables a file descriptor duplication method of
-    redirection.
-    """
-    # Note: It seems we need to capture both streams to make the file
-    # descriptor method work, and we want both anyway throughout Myokit, so
-    # there are no options to choose stdout or stderr here.
-
-    # Lock to stop other threads from capturing while this thread is capturing.
-    _rlock = threading.RLock()
-
-    def __init__(self, fd=False):
-
-        # Are we already capturing? This is needed in case someone enters the
-        # same context twice.
-        self._active_count = 0
-
-        # Captured text
-        self._txt_out = None
-        self._txt_err = None
-
-        # Capturing mode
-        self._fd = bool(fd)
-
-        # Shared by both modes
-        self._org_out = None     # Original stdout object
-        self._org_err = None     # Original stderr object
-
-        # Python stream redirects only
-        self._tmp_out = None    # Temporary stdout StringIO object
-        self._tmp_err = None    # Temporary stderr StringIO object
-
-        # File descriptor redirects only
-        self._out_fd = None      # File descriptor used for output
-        self._err_fd = None      # File descriptor used for errors
-        self._org_out_fd = None  # Back-up of original stdout file descriptor
-        self._org_err_fd = None  # Back-up of original stderr file descriptor
-        self._file_out = None    # Temporary file to write output to
-        self._file_err = None    # Temporary file to write errors to
-
-    def __enter__(self):
-        """Called when the context is entered."""
-        # Avoid entering the same context object twice
-        self._active_count += 1
-        if self._active_count == 1:
-
-            # Wait until other threads have stopped capturing
-            capture._rlock.acquire()
-
-            # Set up redirection
-            self._start()
-
-        # Return context
-        return self
-
-    def __exit__(self, exc_type, exc_value, traceback):
-        """Called when exiting the context."""
-        self._active_count -= 1
-        if self._active_count == 0:
-            self._stop()
-            capture._rlock.release()
-
-    def _start(self):
-        """Starts capturing output to stdout and stderr."""
-
-        # Clear any cached text
-        self._txt_out = None
-        self._txt_err = None
-
-        # Save the current output / error streams
-        self._org_out = sys.stdout
-        self._org_err = sys.stderr
-
-        # Redirect
-        if not self._fd:
-            # Create temporary output and error streams
-            self._tmp_out = StringIO()
-            self._tmp_err = StringIO()
-
-            # Redirect, attempting to flush first
-            try:
-                sys.stdout.flush()
-            except AttributeError:  # pragma: no cover
-                pass
-            finally:
-                sys.stdout = self._tmp_out
-
-            try:
-                sys.stderr.flush()
-            except AttributeError:  # pragma: no cover
-                pass
-            finally:
-                sys.stderr = self._tmp_err
-
-        else:
-            # File-descriptor redirection
-
-            # Get file descriptors used for output and errors.
-            #
-            # These represent open files, and are the low-level equivalent of
-            # sys.stdout and sys.stderr: all subprocesses etc. use these
-            # descriptors to write output to.
-            #
-            # Note that we get these file descriptors from __stdout__ and
-            # __stderr__, to get the original streams that subprocesses etc.
-            # will use, even if someone has redirected the Python-level
-            # streams.
-            #
-            # On https://docs.python.org/3/library/sys.html#module-sys, it says
-            # that stdout/err as well as __stdout__ can be None (e.g. in spyder
-            # on windows).  In other cases (pythonw.exe) they can be set but
-            # return a negative file descriptor (indicating it's invalid).
-            # So here we check if __stdout__ is None and if so set a negative
-            # fileno so that we can catch both cases at once in the rest of the
-            # code.
-            #
-            if sys.__stdout__ is not None:
-                self._out_fd = sys.__stdout__.fileno()
-            else:   # pragma: no cover
-                self._out_fd = -1
-            if sys.__stderr__ is not None:
-                self._err_fd = sys.__stderr__.fileno()
-            else:   # pragma: no cover
-                self._err_fd = -1
-
-            # We start by making a copy of these file descriptors using dup(),
-            # meaning that two file numbers will point to the same file. In a
-            # later step we use dup2() to write a different value into the
-            # descriptor, so that all output gets redirected. When exiting, we
-            # will use dup2() again to copy the original value back in.
-            #
-            # We check if the fd >= 0 first, because the previous code could
-            # get a negative (invalid) value either when __stdxxx__ is None or
-            # when a negative value is returned.
-            #
-            # On windows, the order is important: First dup both stdout and
-            # stderr, then dup2 the new descriptors in. This prevents a weird
-            # infinite recursion on windows ipython / python shell.
-            #
-            if self._out_fd >= 0:
-                self._org_out_fd = os.dup(self._out_fd)
-            if self._err_fd >= 0:
-                self._org_err_fd = os.dup(self._err_fd)
-
-            # Create temporary files to redirect the output to. Make sure they
-            # aren't opened in binary mode, and specify + for reading and
-            # writing.
-            self._file_out = tempfile.TemporaryFile(mode='w+')
-            self._file_err = tempfile.TemporaryFile(mode='w+')
-
-            # Apply Python-level redirection of sys.stdxxx in addition to the
-            # lower-level redirect of sys.__stdxxx___.
-            try:
-                sys.stdout.flush()
-            except AttributeError:  # pragma: no cover
-                pass
-            finally:
-                sys.stdout = self._file_out
-            try:
-                sys.stderr.flush()
-            except AttributeError:  # pragma: no cover
-                pass
-            finally:
-                sys.stderr = self._file_err
-
-            # If possible, overwrite the stdout and stderr file descriptors
-            # with the file descriptor for our files. C/C++ code etc. will
-            # still look in this location for a descriptor, so that output is
-            # redirected.
-            if self._out_fd >= 0:
-                sys.__stdout__.flush()
-                os.dup2(self._file_out.fileno(), self._out_fd)
-            if self._err_fd >= 0:
-                sys.__stderr__.flush()
-                os.dup2(self._file_err.fileno(), self._err_fd)
-
-    def _stop(self):
-        """Stops capturing output."""
-
-        # Flush any remaining output
-        sys.stdout.flush()
-        sys.stderr.flush()
-
-        # Undo redirect
-        if not self._fd:
-            # Direct back at original
-            sys.stdout = self._org_out
-            sys.stderr = self._org_err
-
-            # Store text
-            self._txt_out = self._tmp_out.getvalue()
-            self._txt_err = self._tmp_err.getvalue()
-
-            # Tidy
-            self._tmp_out = self._tmp_err = None
-
-        else:
-            # Undo redirection
-            sys.stdout = self._org_out
-            if self._org_out_fd is not None:
-                # Copy our backup of the original fd into out_fd
-                os.dup2(self._org_out_fd, self._out_fd)
-                # And close our backup
-                os.close(self._org_out_fd)
-
-            sys.stderr = self._org_err
-            if self._org_err_fd is not None:
-                os.dup2(self._org_err_fd, self._err_fd)
-                os.close(self._org_err_fd)
-
-            # Close temporary files and store capture output
-            try:
-                self._file_out.seek(0)
-                self._txt_out = self._file_out.read()
-                self._file_out.close()
-            except ValueError:  # pragma: no cover
-                # In rare cases, I've seen a ValueError, "underlying buffer has
-                # been detached".
-                pass
-            try:
-                self._file_err.seek(0)
-                self._txt_err = self._file_err.read()
-                self._file_err.close()
-            except ValueError:  # pragma: no cover
-                pass
-
-            # Tidy
-            self._out_fd = self._err_fd = None
-            self._org_out_fd = self._org_err_fd = None
-            self._file_out = self._file_err = None
-
-        # Unset reference to original streams (for proper decreffing etc.)
-        self._stdout = self._stderr = None
-
-    def err(self):
-        """
-        Returns the text captured from stderr, or an empty string if nothing
-        was captured or capturing is still active.
-        """
-        if self._txt_err is None:
-            return ''
-        text = self._txt_err
-
-        # In Python 2, the text needs to be decoded from ascii
-        if sys.hexversion < 0x03000000:  # pragma: no python 3 cover
-            text = text.decode('ascii', 'ignore')
-
-        return text
-
-    def out(self):
-        """
-        Returns the text captured from stdout, or an empty string if nothing
-        was captured or capturing is still active.
-        """
-        if self._txt_out is None:
-            return ''
-        text = self._txt_out
-
-        # In Python 2, the text needs to be decoded from ascii
-        if sys.hexversion < 0x03000000:  # pragma: no python 3 cover
-            text = text.decode('ascii', 'ignore')
-
-        return text
-
-    def text(self):
-        """
-        Returns the combined text captured from output and error text, if any
-        (output first, then error text).
-        """
-        return self.out() + self.err()
-
-
 def default_protocol(model=None):
     """
     Returns a default protocol to use when no embedded one is available.
@@ -519,212 +119,35 @@ def default_script(model=None):
     ))
 
 
-def _examplify(filename):
-    """
-    If ``filename`` is equal to "example" and there isn't a file with that
-    name, this function returns the file specified by myokit.EXAMPLE. In all
-    other cases, the original filename is returned.
-    """
-    if filename == 'example' and not os.path.exists(filename):
-        return myokit.EXAMPLE
-    else:
-        return os.path.expanduser(filename)
-
-
 def format_float_dict(d):
     """
-    Takes a dictionary of ``string : float`` mappings and returns a formatted
-    string.
+    Takes a dictionary of ``string : float`` mappings and returns a single
+    multi-line string, where each line is of the form ``key = value``, where
+    key is padded with spaces to make the equals signs align.
+
+    This method is deprecated and will be removed in future versions of Myokit.
     """
+    # Deprecated since 2021-03-24
+    import warnings
+    warnings.warn(
+        'The method `myokit.format_float_dict` is deprecated.'
+        ' It will be removed in future versions of Myokit.')
+
     keys = [str(k) for k in d.keys()]
     keys.sort()
     n = max([len(k) for k in keys])
-    return '\n'.join(
-        [k + ' ' * (n - len(k)) + ' = ' + strfloat(d[k]) for k in keys])
+    return '\n'.join([
+        k + ' ' * (n - len(k)) + ' = ' + myokit.float.str(d[k]) for k in keys])
 
 
 def format_path(path, root='.'):
-    """
-    Formats a path for use in user messages. If the given path is a
-    subdirectory of the current directory this part is chopped off.
-
-    Alternatively, a ``root`` directory may be given explicitly: any
-    subdirectory of this path will be formatted relative to ``root``.
-
-    This function differs from os.path.relpath() in the way it handles paths
-    *outside* the root: In these cases relpath returns a relative path such as
-    '../../' while this function returns an absolute path.
-    """
-    if path == '':
-        path = '.'
-    try:
-        path = os.path.relpath(path, root)
-    except ValueError:  # pragma: no cover
-        # This can happen on windows, if `path` is on a different drive than
-        # root (so that no relative path from one to the other can be made).
-        return path
-
-    if '..' in path:
-        path = os.path.abspath(os.path.join(root, path))
-    return path
-
-
-def load(filename):
-    """
-    Reads an ``mmt`` file and returns a tuple ``(model, protocol, embedded
-    script)``.
-
-    If the file specified by ``filename`` doesn't contain one of these parts
-    the corresponding entry in the tuple will be ``None``.
-    """
-    f = open(_examplify(filename), 'r')
-    try:
-        return myokit.parse(f)
-    finally:
-        f.close()
-
-
-def load_model(filename):
-    """
-    Loads the model section from an ``mmt`` file.
-
-    Raises a :class:`SectionNotFoundError` if no model section is found.
-    """
-    filename = _examplify(filename)
-    with open(filename, 'r') as f:
-        section = myokit.split(f)[0]
-        if not section.strip():
-            raise myokit.SectionNotFoundError('Model section not found.')
-        return myokit.parse(section.splitlines())[0]
-
-
-def load_protocol(filename):
-    """
-    Loads the protocol section from an ``mmt`` file.
-
-    Raises a :class:`SectionNotFoundError` if no protocol section is found.
-    """
-    filename = _examplify(filename)
-    with open(filename, 'r') as f:
-        section = myokit.split(f)[1]
-        if not section.strip():
-            raise myokit.SectionNotFoundError('Protocol section not found.')
-        return myokit.parse(section.splitlines())[1]
-
-
-def load_script(filename):
-    """
-    Loads the script section from an ``mmt`` file.
-
-    Raises a :class:`SectionNotFoundError` if no script section is found.
-    """
-    filename = _examplify(filename)
-    with open(filename, 'r') as f:
-        section = myokit.split(f)[2]
-        if not section.strip():
-            raise myokit.SectionNotFoundError('Script section not found.')
-        return myokit.parse(section.splitlines())[2]
-
-
-def load_state(filename, model=None):
-    """
-    Loads an initial state from a file in one of the formats specified by
-    :func:`myokit.parse_state()`.
-
-    If a :class:`Model` is provided the state will be run through
-    :meth:`Model.map_to_state()` and returned as a list of floating point
-    numbers.
-    """
-    filename = os.path.expanduser(filename)
-    with open(filename, 'r') as f:
-        s = myokit.parse_state(f)
-        if model:
-            s = model.map_to_state(s)
-        return s
-
-
-def load_state_bin(filename):
-    """
-    Loads an initial state from a file in the binary format used by myokit.
-    See :meth:`save_state_bin` for details.
-    """
-    filename = os.path.expanduser(filename)
-
-    # Load compression modules
-    import zipfile
-    try:
-        import zlib
-        del(zlib)
-    except ImportError:
-        raise Exception(
-            'This method requires the `zlib` module to be installed.')
-
-    # Open file
-    with zipfile.ZipFile(filename, 'r') as f:
-        info = f.infolist()
-
-        if len(info) != 1:  # pragma: no cover
-            raise Exception('Invalid state file format [10].')
-
-        # Split into parts, get data type and array size
-        info = info[0]
-        parts = info.filename.split('_')
-
-        if len(parts) != 3:     # pragma: no cover
-            raise Exception('Invalid state file format [20].')
-
-        if parts[0] != 'state':     # pragma: no cover
-            raise Exception('Invalid state file format [30].')
-
-        code = parts[1]
-        if code not in ['d', 'f']:  # pragma: no cover
-            raise Exception('Invalid state file format [40].')
-        # Convert code to str for Python 2.7.10 (see #225)
-        code = str(code)
-
-        size = int(parts[2])
-        if size < 0:    # pragma: no cover
-            raise Exception('Invalid state file format [50].')
-
-        # Create array, read bytes into array
-        ar = array.array(code)
-        try:
-            ar.frombytes(f.read(info))
-        except AttributeError:  # pragma: no python 3 cover
-            ar.fromstring(f.read(info))
-
-        # Always store as little endian
-        if sys.byteorder == 'big':  # pragma: no cover
-            ar.byteswap()
-
-    return list(ar)
-
-
-def _lvsd(s1, s2):
-    """
-    Calculates a Levenshtein distance, as found on wikibooks
-
-    :param s1: The first string to compare
-    :param s2: The second string to compare
-    :returns: int The distance between s1 and s2
-    """
-    if len(s1) < len(s2):
-        return _lvsd(s2, s1)
-
-    if not s1:
-        return len(s2)
-
-    previous_row = range(len(s2) + 1)
-    for i, c1 in enumerate(s1):
-        current_row = [i + 1]
-        for j, c2 in enumerate(s2):
-            insertions = previous_row[j + 1] + 1
-            deletions = current_row[j] + 1
-            substitutions = previous_row[j] + (c1 != c2)
-            current_row.append(min(insertions, deletions, substitutions))
-        previous_row = current_row
-
-    return previous_row[-1]
+    """Deprecated alias of :class:`myokit.tools.format_path`."""
+    # Deprecated since 2021-03-24
+    import warnings
+    warnings.warn(
+        'The method `myokit.format_path` is deprecated.'
+        ' Please use `myokit.tools.format_path` instead.')
+    return myokit.tools.format_path(path, root)
 
 
 class ModelComparison(object):
@@ -1105,156 +528,6 @@ def get_protocol():
     gc.collect()
 
 
-def save(filename=None, model=None, protocol=None, script=None):
-    """
-    Saves a model, protocol, and embedded script to an ``mmt`` file.
-
-    The ``model`` argument can be given as plain text or a
-    :class:`myokit.Model` object. Similarly, ``protocol`` can be either a
-    :class:`myokit.Protocol` or its textual represenation.
-
-    If no filename is given the ``mmt`` code is returned as a string.
-    """
-    if filename:
-        filename = os.path.expanduser(filename)
-        f = open(filename, 'w')
-    else:
-        f = StringIO()
-    out = None
-    try:
-        if model is not None:
-            if isinstance(model, myokit.Model):
-                model = model.code()
-            else:
-                model = model.strip()
-                if model != '' and model[:9] != '[[model]]':
-                    f.write('[[model]]\n')
-            model = model.strip()
-            if model:
-                f.write(model)
-                f.write('\n\n')
-
-        if protocol is not None:
-            if isinstance(protocol, myokit.Protocol):
-                protocol = protocol.code()
-            else:
-                protocol = protocol.strip()
-                if protocol != '' and protocol[:12] != '[[protocol]]':
-                    f.write('[[protocol]]\n')
-            protocol = protocol.strip()
-            if protocol:
-                f.write(protocol)
-                f.write('\n\n')
-
-        if script is not None:
-            script = script.strip()
-            if script != '' and script[:10] != '[[script]]':
-                f.write('[[script]]\n')
-            if script:
-                f.write(script)
-                f.write('\n\n')
-    finally:
-        if filename:
-            f.close()
-        else:
-            out = f.getvalue()
-    return out
-
-
-def save_model(filename, model):
-    """
-    Saves a model to a file
-    """
-    return save(filename, model)
-
-
-def save_protocol(filename, protocol):
-    """
-    Saves a protocol to a file
-    """
-    return save(filename, protocol=protocol)
-
-
-def save_script(filename, script):
-    """
-    Saves an embedded script to a file
-    """
-    return save(filename, script=script)
-
-
-def save_state(filename, state, model=None):
-    """
-    Stores the given state in the file at ``filename``.
-
-    If no ``model`` is specified ``state`` should be given as a list of
-    floating point numbers and will be stored by simply placing each number on
-    a new line.
-
-    If a :class:`Model <myokit.Model>` is provided the state can be in any
-    format accepted by :meth:`Model.map_to_state() <myokit.Model.map_to_state>`
-    and will be stored in the format returned by
-    :meth:`Model.format_state() <myokit.Model.format_state>`.
-    """
-    # Check filename
-    filename = os.path.expanduser(filename)
-
-    # Format
-    if model is not None:
-        state = model.map_to_state(state)
-        state = model.format_state(state)
-    else:
-        state = [myokit.strfloat(s) for s in state]
-        state = '\n'.join(state)
-
-    # Store
-    with open(filename, 'w') as f:
-        f.write(state)
-
-
-def save_state_bin(filename, state, precision=myokit.DOUBLE_PRECISION):
-    """
-    Stores the given state (or any list of floating point numbers) in the file
-    at ``filename``, using a binary format.
-
-    The used format is a zip file, containing a single entry: ``state_x_y``,
-    where ``x`` is the used data type (``d`` or ``f``) and ``y`` is the number
-    of entries. All entries are stored little-endian.
-    """
-    # Check filename
-    filename = os.path.expanduser(filename)
-
-    # Load compression modules
-    import zipfile
-    try:
-        import zlib
-        del(zlib)
-    except ImportError:
-        raise Exception(
-            'This method requires the `zlib` module to be installed.')
-
-    # Data type
-    # Convert code to str for Python 2.7.10 (see #225)
-    code = str('d' if precision == myokit.DOUBLE_PRECISION else 'f')
-
-    # Create array, ensure it's little-endian
-    ar = array.array(code, state)
-    if sys.byteorder == 'big':  # pragma: no cover
-        ar.byteswap()
-
-    # Store precision and data type in internal filename
-    name = 'state_' + code + '_' + str(len(state))
-    info = zipfile.ZipInfo(name)
-    info.compress_type = zipfile.ZIP_DEFLATED
-
-    # Write to compressed file
-    try:
-        ar = ar.tobytes()
-    except AttributeError:  # pragma: no python 3 cover
-        ar = ar.tostring()
-    with zipfile.ZipFile(filename, 'w') as f:
-        f.writestr(info, ar)
-
-
 def step(model, initial=None, reference=None, ignore_errors=False):
     """
     Evaluates the state derivatives in a model and compares the results with a
@@ -1330,7 +603,7 @@ def step(model, initial=None, reference=None, ignore_errors=False):
             if xx[0] != yy[0]:
 
                 # Ignore if zero
-                if (_feq(x, 0) and _feq(y, 0)):
+                if (myokit.float.eq(x, 0) and myokit.float.eq(y, 0)):
                     log.append(line)
                     log.append('')
                 else:
@@ -1397,35 +670,21 @@ def step(model, initial=None, reference=None, ignore_errors=False):
 
 
 def strfloat(number, full=False, precision=myokit.DOUBLE_PRECISION):
-    """
-    Turns the given number into a string.
-    """
-    # Force full precision output
-    if full:
-        if precision == myokit.SINGLE_PRECISION:
-            return myokit.SFSINGLE.format(float(number))
-        else:
-            return myokit.SFDOUBLE.format(float(number))
-
-    # Pass through strings
-    if isinstance(number, str):
-        return number
-
-    # Handle myokit.Numbers
-    if isinstance(number, myokit.Number):
-        number = number.eval()
+    """Deprecated alias of :class:`myokit.float.str`."""
+    # Deprecated since 2021-03-24
+    import warnings
+    warnings.warn(
+        'The method `myokit.strfloat` is deprecated.'
+        ' Please use `myokit.float.str` instead.')
+    return myokit.float.str(number, full, precision)
 
-    # For most numbers, allow python to format the float
-    s = str(number)
-    if len(s) < 10:
-        return s
 
-    # But if the number is given with lots of decimals, use the representation
-    # with enough digits to prevent loss of information
-    if precision == myokit.SINGLE_PRECISION:
-        return myokit.SFSINGLE.format(float(number))
-    else:
-        return myokit.SFDOUBLE.format(float(number))
+def time():
+    """
+    Returns the current time, in the format used throughout Myokit.
+    """
+    import time as t
+    return t.strftime(myokit.TIME_FORMAT)
 
 
 def version(raw=False):
@@ -1450,82 +709,3 @@ def version(raw=False):
         t2 += '|  |' + '_' * 5
         return '\n' + t1 + '\n' + t2
 
-
-def _feq(a, b):
-    """
-    Checks if floating point numbers ``a`` and ``b`` are equal, or so close to
-    each other that the difference could be a single rounding error.
-    """
-    # Note the initial == check handles infinity
-    return a == b or abs(a - b) < max(abs(a), abs(b)) * sys.float_info.epsilon
-
-
-def _fgeq(a, b):
-    """
-    Checks if ``a >= b``, but using :meth:`myokit._feq` instead of ``=``.
-    """
-    return a >= b or abs(a - b) < max(abs(a), abs(b)) * sys.float_info.epsilon
-
-
-def _fround(x):
-    """
-    Checks if a float ``x`` is within 1 rounding error of an integer, and if
-    so, converts it to that integer.
-    """
-    ix = int(round(x))
-    return ix if _feq(x, ix) else x
-
-
-def _close(a, b, reltol=1e-9, abstol=1e-9):
-    """
-    Test whether two numbers are close enough to be considered equal.
-
-    Differs from :meth:`myokit._feq` in that it tries to answer the question
-    "are two number resulting from various calculations close enough to be
-    considered equal". Whereas `_feq` aims to deal with numbers that are
-    numerically indistinguishable but still have a slightly different floating
-    point representation.
-    """
-    # Note the initial == check handles infinity
-    return a == b or abs(a - b) < max(reltol * max(abs(a), abs(b)), abstol)
-
-
-def _cround(x, reltol=1e-9, abstol=1e-9):
-    """
-    Checks if a float ``x`` is close to an integer with :meth:`close()`, and if
-    so, converts it to that integer.
-    """
-    ix = int(round(x))
-    return ix if _close(x, ix, reltol, abstol) else x
-
-
-def _rmtree(path):
-    """
-    Version of ``shutil.rmtree`` that handles access denied errors (when the
-    user is lacking write permissions). This seems to happen on Windows some
-    times.
-
-    The solution here is based on answers given on stackoverflow:
-    https://stackoverflow.com/questions/2656322
-    """
-    def onerror(function, path, excinfo):   # pragma: no cover
-        if not os.access(path, os.W_OK):
-            # Give user write permissions (remove read-only flag)
-            os.chmod(path, stat.S_IWUSR)
-            function(path)
-        else:
-            raise
-
-    shutil.rmtree(path, ignore_errors=False, onerror=onerror)
-
-
-def _pid_hash():
-    """
-    Returns an integer hash that depends on the current time as well as the
-    process id, so that it's likely to return a different number when called
-    twice.
-    """
-    x = os.getpid() * timeit.default_timer()
-    x = abs(hash(str(x - int(x))))
-    return x
-
diff --git a/myokit/_datablock.py b/myokit/_datablock.py
index 7066e84c0..e49dd9bb8 100644
--- a/myokit/_datablock.py
+++ b/myokit/_datablock.py
@@ -1532,7 +1532,7 @@ def save_frame_csv(
         text = [delimx.join('"' + str(x) + '"' for x in [xname, yname, zname])]
         for y, row in enumerate(data):
             for x, z in enumerate(row):
-                text.append(delimx.join([str(x), str(y), myokit.strfloat(z)]))
+                text.append(delimx.join([str(x), str(y), myokit.float.str(z)]))
         text = delimy.join(text)
         with open(filename, 'w') as f:
             f.write(text)
@@ -1555,7 +1555,7 @@ def save_frame_grid(self, filename, name, frame, delimx=' ', delimy='\n'):
         data = data[frame]
         text = []
         for row in data:
-            text.append(delimx.join([myokit.strfloat(x) for x in row]))
+            text.append(delimx.join([myokit.float.str(x) for x in row]))
         text = delimy.join(text)
         with open(filename, 'w') as f:
             f.write(text)
diff --git a/myokit/_expressions.py b/myokit/_expressions.py
index 6bec32d84..a66652c25 100644
--- a/myokit/_expressions.py
+++ b/myokit/_expressions.py
@@ -811,7 +811,7 @@ def __init__(self, value, unit=None):
                 raise ValueError(
                     'Unit in myokit.Number should be a myokit.Unit or None.')
         # Create nice string representation
-        self._str = myokit.strfloat(self._value)
+        self._str = myokit.float.str(self._value)
         if self._str[-2:] == '.0':
             # Turn 5.0 into 5
             self._str = self._str[:-2]
diff --git a/myokit/_io.py b/myokit/_io.py
new file mode 100644
index 000000000..b11b5b511
--- /dev/null
+++ b/myokit/_io.py
@@ -0,0 +1,315 @@
+#
+# Myokit auxillary functions: This module can be used to gather any
+# functions that are important enough to warrant inclusion in the main
+# myokit module but don't belong to any specific hidden module.
+#
+# This file is part of Myokit.
+# See http://myokit.org for copyright, sharing, and licensing details.
+#
+from __future__ import absolute_import, division
+from __future__ import print_function, unicode_literals
+
+import array
+import os
+import sys
+
+# StringIO in Python 2 and 3
+try:
+    from cStringIO import StringIO
+except ImportError:  # pragma: no python 2 cover
+    from io import StringIO
+
+import myokit
+
+
+def _examplify(filename):
+    """
+    If ``filename`` is equal to "example" and there isn't a file with that
+    name, this function returns the file specified by myokit.EXAMPLE. In all
+    other cases, the original filename is returned.
+    """
+    if filename == 'example' and not os.path.exists(filename):
+        return myokit.EXAMPLE
+    else:
+        return os.path.expanduser(filename)
+
+
+def load(filename):
+    """
+    Reads an ``mmt`` file and returns a tuple ``(model, protocol, embedded
+    script)``.
+
+    If the file specified by ``filename`` doesn't contain one of these parts
+    the corresponding entry in the tuple will be ``None``.
+    """
+    f = open(_examplify(filename), 'r')
+    try:
+        return myokit.parse(f)
+    finally:
+        f.close()
+
+
+def load_model(filename):
+    """
+    Loads the model section from an ``mmt`` file.
+
+    Raises a :class:`SectionNotFoundError` if no model section is found.
+    """
+    filename = _examplify(filename)
+    with open(filename, 'r') as f:
+        section = myokit.split(f)[0]
+        if not section.strip():
+            raise myokit.SectionNotFoundError('Model section not found.')
+        return myokit.parse(section.splitlines())[0]
+
+
+def load_protocol(filename):
+    """
+    Loads the protocol section from an ``mmt`` file.
+
+    Raises a :class:`SectionNotFoundError` if no protocol section is found.
+    """
+    filename = _examplify(filename)
+    with open(filename, 'r') as f:
+        section = myokit.split(f)[1]
+        if not section.strip():
+            raise myokit.SectionNotFoundError('Protocol section not found.')
+        return myokit.parse(section.splitlines())[1]
+
+
+def load_script(filename):
+    """
+    Loads the script section from an ``mmt`` file.
+
+    Raises a :class:`SectionNotFoundError` if no script section is found.
+    """
+    filename = _examplify(filename)
+    with open(filename, 'r') as f:
+        section = myokit.split(f)[2]
+        if not section.strip():
+            raise myokit.SectionNotFoundError('Script section not found.')
+        return myokit.parse(section.splitlines())[2]
+
+
+def load_state(filename, model=None):
+    """
+    Loads an initial state from a file in one of the formats specified by
+    :func:`myokit.parse_state()`.
+
+    If a :class:`Model` is provided the state will be run through
+    :meth:`Model.map_to_state()` and returned as a list of floating point
+    numbers.
+    """
+    filename = os.path.expanduser(filename)
+    with open(filename, 'r') as f:
+        s = myokit.parse_state(f)
+        if model:
+            s = model.map_to_state(s)
+        return s
+
+
+def load_state_bin(filename):
+    """
+    Loads an initial state from a file in the binary format used by myokit.
+    See :meth:`save_state_bin` for details.
+    """
+    filename = os.path.expanduser(filename)
+
+    # Load compression modules
+    import zipfile
+    try:
+        import zlib
+        del(zlib)
+    except ImportError:
+        raise Exception(
+            'This method requires the `zlib` module to be installed.')
+
+    # Open file
+    with zipfile.ZipFile(filename, 'r') as f:
+        info = f.infolist()
+
+        if len(info) != 1:  # pragma: no cover
+            raise Exception('Invalid state file format [10].')
+
+        # Split into parts, get data type and array size
+        info = info[0]
+        parts = info.filename.split('_')
+
+        if len(parts) != 3:     # pragma: no cover
+            raise Exception('Invalid state file format [20].')
+
+        if parts[0] != 'state':     # pragma: no cover
+            raise Exception('Invalid state file format [30].')
+
+        code = parts[1]
+        if code not in ['d', 'f']:  # pragma: no cover
+            raise Exception('Invalid state file format [40].')
+        # Convert code to str for Python 2.7.10 (see #225)
+        code = str(code)
+
+        size = int(parts[2])
+        if size < 0:    # pragma: no cover
+            raise Exception('Invalid state file format [50].')
+
+        # Create array, read bytes into array
+        ar = array.array(code)
+        try:
+            ar.frombytes(f.read(info))
+        except AttributeError:  # pragma: no python 3 cover
+            ar.fromstring(f.read(info))
+
+        # Always store as little endian
+        if sys.byteorder == 'big':  # pragma: no cover
+            ar.byteswap()
+
+    return list(ar)
+
+
+def save(filename=None, model=None, protocol=None, script=None):
+    """
+    Saves a model, protocol, and embedded script to an ``mmt`` file.
+
+    The ``model`` argument can be given as plain text or a
+    :class:`myokit.Model` object. Similarly, ``protocol`` can be either a
+    :class:`myokit.Protocol` or its textual represenation.
+
+    If no filename is given the ``mmt`` code is returned as a string.
+    """
+    if filename:
+        filename = os.path.expanduser(filename)
+        f = open(filename, 'w')
+    else:
+        f = StringIO()
+    out = None
+    try:
+        if model is not None:
+            if isinstance(model, myokit.Model):
+                model = model.code()
+            else:
+                model = model.strip()
+                if model != '' and model[:9] != '[[model]]':
+                    f.write('[[model]]\n')
+            model = model.strip()
+            if model:
+                f.write(model)
+                f.write('\n\n')
+
+        if protocol is not None:
+            if isinstance(protocol, myokit.Protocol):
+                protocol = protocol.code()
+            else:
+                protocol = protocol.strip()
+                if protocol != '' and protocol[:12] != '[[protocol]]':
+                    f.write('[[protocol]]\n')
+            protocol = protocol.strip()
+            if protocol:
+                f.write(protocol)
+                f.write('\n\n')
+
+        if script is not None:
+            script = script.strip()
+            if script != '' and script[:10] != '[[script]]':
+                f.write('[[script]]\n')
+            if script:
+                f.write(script)
+                f.write('\n\n')
+    finally:
+        if filename:
+            f.close()
+        else:
+            out = f.getvalue()
+    return out
+
+
+def save_model(filename, model):
+    """
+    Saves a model to a file
+    """
+    return save(filename, model)
+
+
+def save_protocol(filename, protocol):
+    """
+    Saves a protocol to a file
+    """
+    return save(filename, protocol=protocol)
+
+
+def save_script(filename, script):
+    """
+    Saves an embedded script to a file
+    """
+    return save(filename, script=script)
+
+
+def save_state(filename, state, model=None):
+    """
+    Stores the given state in the file at ``filename``.
+
+    If no ``model`` is specified ``state`` should be given as a list of
+    floating point numbers and will be stored by simply placing each number on
+    a new line.
+
+    If a :class:`Model <myokit.Model>` is provided the state can be in any
+    format accepted by :meth:`Model.map_to_state() <myokit.Model.map_to_state>`
+    and will be stored in the format returned by
+    :meth:`Model.format_state() <myokit.Model.format_state>`.
+    """
+    # Check filename
+    filename = os.path.expanduser(filename)
+
+    # Format
+    if model is not None:
+        state = model.map_to_state(state)
+        state = model.format_state(state)
+    else:
+        state = '\n'.join([myokit.float.str(s) for s in state])
+
+    # Store
+    with open(filename, 'w') as f:
+        f.write(state)
+
+
+def save_state_bin(filename, state, precision=myokit.DOUBLE_PRECISION):
+    """
+    Stores the given state (or any list of floating point numbers) in the file
+    at ``filename``, using a binary format.
+
+    The used format is a zip file, containing a single entry: ``state_x_y``,
+    where ``x`` is the used data type (``d`` or ``f``) and ``y`` is the number
+    of entries. All entries are stored little-endian.
+    """
+    # Check filename
+    filename = os.path.expanduser(filename)
+
+    # Load compression modules
+    import zipfile
+    try:
+        import zlib
+        del(zlib)
+    except ImportError:
+        raise Exception(
+            'This method requires the `zlib` module to be installed.')
+
+    # Data type
+    # Convert code to str for Python 2.7.10 (see #225)
+    code = str('d' if precision == myokit.DOUBLE_PRECISION else 'f')
+
+    # Create array, ensure it's little-endian
+    ar = array.array(code, state)
+    if sys.byteorder == 'big':  # pragma: no cover
+        ar.byteswap()
+
+    # Store precision and data type in internal filename
+    name = 'state_' + code + '_' + str(len(state))
+    info = zipfile.ZipInfo(name)
+    info.compress_type = zipfile.ZIP_DEFLATED
+
+    # Write to compressed file
+    try:
+        ar = ar.tobytes()
+    except AttributeError:  # pragma: no python 3 cover
+        ar = ar.tostring()
+    with zipfile.ZipFile(filename, 'w') as f:
+        f.writestr(info, ar)
+
diff --git a/myokit/_model_api.py b/myokit/_model_api.py
index 3595121c8..e7207ec7e 100644
--- a/myokit/_model_api.py
+++ b/myokit/_model_api.py
@@ -1491,8 +1491,8 @@ def format_state(self, state=None, state2=None,
             An optional second state, to be shown next to ``state`` for
             comparison.
         ``precision=myokit.DOUBLE_PRECISION``
-            An optional precision argument to pass into :meth:`myokit.strfloat`
-            when formatting the state values.
+            An optional precision argument to pass into
+            :meth:`myokit.float.str` when formatting the state values.
 
         """
         n = len(self._state)
@@ -1515,13 +1515,13 @@ def format_state(self, state=None, state2=None,
         for k, var in enumerate(self.states()):
             out.append(
                 var.qname() + ' ' * (n - len(var.qname()))
-                + ' = ' + myokit.strfloat(state[k], precision=precision))
+                + ' = ' + myokit.float.str(state[k], precision=precision))
         if state2 is not None:
             n = max([len(x) for x in out])
             for k, var in enumerate(self.states()):
                 out[k] += (
                     ' ' * (4 + n - len(out[k]))
-                    + myokit.strfloat(state2[k], precision=precision))
+                    + myokit.float.str(state2[k], precision=precision))
 
         return '\n'.join(out)
 
@@ -1542,7 +1542,7 @@ def format_state_derivatives(self, state=None, derivatives=None,
             will be calculed from ``state`` using :meth:`eval_derivatives()`.
         ``precision=myokit.DOUBLE_PRECISION``
             An optional precision argument to use when evaluating the state
-            derivatives, and to pass into :meth:`myokit.strfloat` when
+            derivatives, and to pass into :meth:`myokit.float.str` when
             formatting the state values and derivatives.
 
         """
@@ -1565,8 +1565,8 @@ def format_state_derivatives(self, state=None, derivatives=None,
         out = []
         n = max([len(x.qname()) for x in self.states()])
         for i, var in enumerate(self.states()):
-            s = myokit.strfloat(state[i], precision=precision)
-            d = myokit.strfloat(derivatives[i], precision=precision)
+            s = myokit.float.str(state[i], precision=precision)
+            d = myokit.float.str(derivatives[i], precision=precision)
             out.append(
                 var.qname() + ' ' * (n - len(var.qname())) + ' = ' + s
                 + ' ' * (24 - len(s)) + '   dot = ' + d)
@@ -2930,10 +2930,10 @@ def suggest_variable(self, name):
         for v in self.variables(deep=True):
             n1 = v.name()
             n2 = v.qname()
-            d = min(myokit._lvsd(name, n1),
-                    myokit._lvsd(qname, n2),
-                    myokit._lvsd(name_low, n1.lower()),
-                    myokit._lvsd(qname_low, n2.lower()))
+            d = min(myokit.tools.lvsd(name, n1),
+                    myokit.tools.lvsd(qname, n2),
+                    myokit.tools.lvsd(name_low, n1.lower()),
+                    myokit.tools.lvsd(qname_low, n2.lower()))
             if d < mn:
                 mn = d
                 sg = v
diff --git a/myokit/_progress.py b/myokit/_progress.py
index c087cae87..7c72ebb79 100644
--- a/myokit/_progress.py
+++ b/myokit/_progress.py
@@ -98,7 +98,7 @@ class ProgressPrinter(ProgressReporter):
     """
     def __init__(self, digits=1):
         super(ProgressPrinter, self).__init__()
-        self._b = myokit.Benchmarker()
+        self._b = myokit.tools.Benchmarker()
         self._f = None
         self._d = int(digits)
 
diff --git a/myokit/_protocol.py b/myokit/_protocol.py
index 54c63ad48..f86d8165a 100644
--- a/myokit/_protocol.py
+++ b/myokit/_protocol.py
@@ -94,8 +94,9 @@ def characteristic_time(self):
         Returns the characteristic time associated with this protocol.
 
         The characteristic time is defined as the maximum characteristic time
-        of all events in the protocol. For a sequence of events, this is simply
-        the protocol duration.
+        of all events in the protocol
+        (see :meth:`ProtocolEvent.characteristic_time()`). For a sequence of
+        events, this is simply the protocol duration.
 
         Examples:
 
@@ -283,7 +284,7 @@ def is_sequence_exception(self):
 
             # Calculated position indistinguishable from user-specified next
             # even start? Then jump there instead
-            if e and myokit._feq(t, e._start):
+            if e and myokit.float.eq(t, e._start):
                 t = e._start
 
         return True
@@ -326,7 +327,7 @@ def is_unbroken_sequence_exception(self):
 
             # Calculated position indistinguishable from user-specified next
             # even start? Then jump there instead
-            if myokit._feq(t, e._start):
+            if myokit.float.eq(t, e._start):
                 t = e._start
 
             # Check for periodic events
@@ -794,16 +795,16 @@ def advance(self, new_time):
         self._time = new_time
 
         # Advance pacing system
-        while myokit._fgeq(new_time, self._tnext):
+        while myokit.float.geq(new_time, self._tnext):
 
             # Active event finished
-            if self._fire and myokit._fgeq(self._tnext, self._tdown):
+            if self._fire and myokit.float.geq(self._tnext, self._tdown):
                 self._fire = None
                 self._pace = 0
 
             # New event starting
             e = self._protocol._head
-            if e and myokit._fgeq(new_time, e._start):
+            if e and myokit.float.geq(new_time, e._start):
                 self._protocol.pop()
                 self._fire = e
                 self._tdown = e._start + e._duration
@@ -820,7 +821,7 @@ def advance(self, new_time):
                 # If so, then set tdown (which is always calculated) to the
                 # next event start (which may be user-specified).
                 x = self._protocol._head
-                if x and myokit._feq(self._tdown, x._start):
+                if x and myokit.float.eq(self._tdown, x._start):
                     self._tdown = x._start
 
             # Next stopping time
diff --git a/myokit/_sim/__init__.py b/myokit/_sim/__init__.py
index 5b6e779a1..a377611d4 100644
--- a/myokit/_sim/__init__.py
+++ b/myokit/_sim/__init__.py
@@ -13,6 +13,8 @@
 import platform
 import sys
 import tempfile
+import threading
+import timeit
 import traceback
 
 
@@ -187,7 +189,7 @@ def _compile(self, name, tpl, tpl_vars, libs, libd=None, incd=None,
 
             # Compile in build directory, catch output
             error, trace = None, None
-            with myokit.capture(fd=True) as s:
+            with myokit.tools.capture(fd=True) as s:
                 try:
                     os.chdir(d_build)
                     setup(
@@ -219,7 +221,7 @@ def _compile(self, name, tpl, tpl_vars, libs, libd=None, incd=None,
 
             # Delete cached module
             try:
-                myokit._rmtree(d_cache)
+                myokit.tools.rmtree(d_cache)
             except Exception:   # pragma: no cover
                 pass
 
@@ -279,3 +281,15 @@ class CppModule(CModule):
     def _source_file(self):
         return 'source.cpp'
 
+
+def pid_hash():
+    """
+    Returns a positive integer hash that depends on the current time as well as
+    the process and thread id, so that it's likely to return a different number
+    when called twice.
+    """
+    pid = 1 + os.getpid()                       # Range 0 to 99999
+    tid = threading.current_thread().ident      # Non-zero integer
+    x = pid * tid * timeit.default_timer()
+    return abs(hash(str(x - int(x))))
+
diff --git a/myokit/_sim/cable.py b/myokit/_sim/cable.py
index c26539173..8af6f6393 100644
--- a/myokit/_sim/cable.py
+++ b/myokit/_sim/cable.py
@@ -182,7 +182,7 @@ def __init__(self, model, protocol=None, ncells=50, rl=False):
         # Unique simulation id
         Simulation1d._index += 1
         module_name = 'myokit_sim1d_' + str(Simulation1d._index)
-        module_name += '_' + str(myokit._pid_hash())
+        module_name += '_' + str(myokit.pid_hash())
 
         # Arguments
         args = {
@@ -196,8 +196,8 @@ def __init__(self, model, protocol=None, ncells=50, rl=False):
 
         # Debug
         if myokit.DEBUG:
-            print(self._code(fname, args,
-                             line_numbers=myokit.DEBUG_LINE_NUMBERS))
+            print(self._code(
+                fname, args, line_numbers=myokit.DEBUG_LINE_NUMBERS))
             import sys
             sys.exit(1)
 
diff --git a/myokit/_sim/cmodel.h b/myokit/_sim/cmodel.h
index a1ed274cb..602c2714e 100644
--- a/myokit/_sim/cmodel.h
+++ b/myokit/_sim/cmodel.h
@@ -1228,7 +1228,7 @@ for eq in parameters.values():
     /* States */
 <?
 for var in model.states():
-    print(tab + v(var) + ' = ' + myokit.strfloat(var.state_value()) + ';')
+    print(tab + v(var) + ' = ' + myokit.float.str(var.state_value()) + ';')
 ?>
     /*
      * Caching.
diff --git a/myokit/_sim/compiler.py b/myokit/_sim/compiler.py
index 5e2caebd8..31960ce03 100644
--- a/myokit/_sim/compiler.py
+++ b/myokit/_sim/compiler.py
@@ -44,7 +44,7 @@ def __init__(self):
 
         # Create back-end
         mname = 'myokit_compiler_info_' + str(Compiler._index)
-        mname += '_' + str(myokit._pid_hash())
+        mname += '_' + str(myokit.pid_hash())
         fname = os.path.join(myokit.DIR_CFUNC, SOURCE_FILE)
         args = {'module_name': mname}
         try:
diff --git a/myokit/_sim/cvodesim.c b/myokit/_sim/cvodesim.c
index 7d9ea8d07..fe6dfec4a 100644
--- a/myokit/_sim/cvodesim.c
+++ b/myokit/_sim/cvodesim.c
@@ -212,7 +212,7 @@ static long engine_steps = 0;
 <?
 for var in model.variables(state=False, deep=True):
     if var.is_literal():
-        print('static realtype ' + v(var) + ' = ' + myokit.strfloat(var.rhs().eval()) + ';')
+        print('static realtype ' + v(var) + ' = ' + myokit.float.str(var.rhs().eval()) + ';')
     else:
         print('static realtype ' + v(var) + ';')
 ?>
diff --git a/myokit/_sim/cvodesim.py b/myokit/_sim/cvodesim.py
index 23a5eb84c..48abaec33 100644
--- a/myokit/_sim/cvodesim.py
+++ b/myokit/_sim/cvodesim.py
@@ -116,7 +116,7 @@ def __init__(self, model, protocol=None, apd_var=None):
         # Unique simulation id
         Simulation._index += 1
         module_name = 'myokit_sim_' + str(Simulation._index)
-        module_name += '_' + str(myokit._pid_hash())
+        module_name += '_' + str(myokit.pid_hash())
 
         # Arguments
         args = {
@@ -438,10 +438,10 @@ def _run(
                 txt.extend(['  ' + x for x in
                             self._model.format_state(rstate).splitlines()])
                 txt.append('Inputs for binding: ')
-                txt.append('  time        = ' + myokit.strfloat(rbound[0]))
-                txt.append('  pace        = ' + myokit.strfloat(rbound[1]))
-                txt.append('  realtime    = ' + myokit.strfloat(rbound[2]))
-                txt.append('  evaluations = ' + myokit.strfloat(rbound[3]))
+                txt.append('  time        = ' + myokit.float.str(rbound[0]))
+                txt.append('  pace        = ' + myokit.float.str(rbound[1]))
+                txt.append('  realtime    = ' + myokit.float.str(rbound[2]))
+                txt.append('  evaluations = ' + myokit.float.str(rbound[3]))
                 txt.append(str(e))
                 try:
                     self._model.eval_state_derivatives(rstate)
diff --git a/myokit/_sim/cvodessim.py b/myokit/_sim/cvodessim.py
index 41079e931..b9f166f43 100644
--- a/myokit/_sim/cvodessim.py
+++ b/myokit/_sim/cvodessim.py
@@ -187,7 +187,7 @@ def _create_simulation(self, cmodel_code):
         # Unique simulation id
         Simulation._index += 1
         module_name = 'myokit_sim_' + str(Simulation._index)
-        module_name += '_' + str(myokit._pid_hash())
+        module_name += '_' + str(myokit.pid_hash())
 
         # Arguments
         args = {
@@ -666,10 +666,10 @@ def _run(self, duration, log, log_interval, log_times, sensitivities,
                 txt.extend(['  ' + x for x
                             in self._model.format_state(state).splitlines()])
                 txt.append('Inputs for binding:')
-                txt.append('  time        = ' + myokit.strfloat(bound[0]))
-                txt.append('  pace        = ' + myokit.strfloat(bound[1]))
-                txt.append('  realtime    = ' + myokit.strfloat(bound[2]))
-                txt.append('  evaluations = ' + myokit.strfloat(bound[3]))
+                txt.append('  time        = ' + myokit.float.str(bound[0]))
+                txt.append('  pace        = ' + myokit.float.str(bound[1]))
+                txt.append('  realtime    = ' + myokit.float.str(bound[2]))
+                txt.append('  evaluations = ' + myokit.float.str(bound[3]))
                 txt.append(str(e))
 
                 # Check if state derivatives can be evaluated in Python, if
diff --git a/myokit/_sim/fiber_tissue.py b/myokit/_sim/fiber_tissue.py
index 1dc53b3de..58220787d 100644
--- a/myokit/_sim/fiber_tissue.py
+++ b/myokit/_sim/fiber_tissue.py
@@ -341,7 +341,7 @@ def _create_backend(self):
         # Unique simulation id
         FiberTissueSimulation._index += 1
         mname = 'myokit_sim_fiber_tissue_' + str(FiberTissueSimulation._index)
-        mname += '_' + str(myokit._pid_hash())
+        mname += '_' + str(myokit.pid_hash())
 
         # Arguments
         args = {
diff --git a/myokit/_sim/icsim.py b/myokit/_sim/icsim.py
index 769865b5e..a69eb63d8 100644
--- a/myokit/_sim/icsim.py
+++ b/myokit/_sim/icsim.py
@@ -117,7 +117,7 @@ def __init__(self, model, protocol=None):
         # Unique simulation id
         ICSimulation._index += 1
         module_name = 'myokit_ICSimulation_' + str(ICSimulation._index)
-        module_name += '_' + str(myokit._pid_hash())
+        module_name += '_' + str(myokit.pid_hash())
 
         # Arguments
         args = {
diff --git a/myokit/_sim/jacobian.py b/myokit/_sim/jacobian.py
index 0c1fb7b43..a4a9ffd23 100644
--- a/myokit/_sim/jacobian.py
+++ b/myokit/_sim/jacobian.py
@@ -59,7 +59,7 @@ def __init__(self, model):
         # Extension module id
         JacobianTracer._index += 1
         module_name = 'myokit_JacobianTracer_' + str(JacobianTracer._index)
-        module_name += '_' + str(myokit._pid_hash())
+        module_name += '_' + str(myokit.pid_hash())
 
         # Template arguments
         args = {
@@ -238,7 +238,7 @@ def __init__(self, model):
         JacobianCalculator._index += 1
         module_name = 'myokit_JacobianCalculator_'
         module_name += str(JacobianCalculator._index)
-        module_name += '_' + str(myokit._pid_hash())
+        module_name += '_' + str(myokit.pid_hash())
 
         # Template arguments
         args = {
diff --git a/myokit/_sim/opencl.py b/myokit/_sim/opencl.py
index 0620c3cd4..a08074c2c 100644
--- a/myokit/_sim/opencl.py
+++ b/myokit/_sim/opencl.py
@@ -45,7 +45,7 @@ def __init__(self):
         # Create back-end and cache it
         OpenCL._index += 1
         mname = 'myokit_opencl_info_' + str(OpenCL._index)
-        mname += '_' + str(myokit._pid_hash())
+        mname += '_' + str(myokit.pid_hash())
         fname = os.path.join(myokit.DIR_CFUNC, SOURCE_FILE)
         args = {'module_name': mname}
 
diff --git a/myokit/_sim/openclsim.py b/myokit/_sim/openclsim.py
index 7f23d2949..673f1c011 100644
--- a/myokit/_sim/openclsim.py
+++ b/myokit/_sim/openclsim.py
@@ -286,7 +286,7 @@ def __init__(
         # Create back-end
         SimulationOpenCL._index += 1
         mname = 'myokit_sim_opencl_' + str(SimulationOpenCL._index)
-        mname += '_' + str(myokit._pid_hash())
+        mname += '_' + str(myokit.pid_hash())
         fname = os.path.join(myokit.DIR_CFUNC, SOURCE_FILE)
         args = {
             'module_name': mname,
@@ -1049,9 +1049,9 @@ def _run(self, duration, log, log_interval, report_nan, progress, msg):
                 icell_str = '(' + ','.join([str(x) for x in icell]) + ')'
                 txt.append(
                     'Encountered numerical error at t='
-                    + myokit.strfloat(time, precision=self._precision)
+                    + myokit.float.str(time, precision=self._precision)
                     + ' in cell ' + icell_str + ' when ' + var + '='
-                    + myokit.strfloat(value, precision=self._precision) + '.')
+                    + myokit.float.str(value, precision=self._precision) + '.')
 
                 # Check if this cell was paced.
                 is_paced = self.is_paced(*icell)
@@ -1084,16 +1084,16 @@ def _run(self, duration, log, log_interval, report_nan, progress, msg):
 
                 # Show bound variables
                 txt.append('Simulation variables during:')
-                txt.append('  Time: ' + myokit.strfloat(
+                txt.append('  Time: ' + myokit.float.str(
                     bounds[0][vtime], precision=self._precision))
                 if vpace:
                     txt.append(
-                        '  Pacing variable: ' + myokit.strfloat(
+                        '  Pacing variable: ' + myokit.float.str(
                             bounds[0][vpace], precision=self._precision)
                         + is_paced_str)
                 if vdiff:
                     txt.append(
-                        '  Diffusion current: ' + myokit.strfloat(
+                        '  Diffusion current: ' + myokit.float.str(
                             bounds[0][vdiff], precision=self._precision))
                 if neighbours:
                     txt.append('  Connected cells: ' + neighbours_str)
@@ -1131,16 +1131,16 @@ def _run(self, duration, log, log_interval, report_nan, progress, msg):
 
                     # Show bound variables
                     txt.append('Simulation variables before:')
-                    txt.append('  Time: ' + myokit.strfloat(
+                    txt.append('  Time: ' + myokit.float.str(
                         bounds[1][vtime], precision=self._precision))
                     if vpace:
                         txt.append(
-                            '  Pacing variable: ' + myokit.strfloat(
+                            '  Pacing variable: ' + myokit.float.str(
                                 bounds[1][vpace], precision=self._precision)
                             + is_paced_str)
                     if vdiff:
                         txt.append(
-                            '  Diffusion current: ' + myokit.strfloat(
+                            '  Diffusion current: ' + myokit.float.str(
                                 bounds[1][vdiff], precision=self._precision))
                     if neighbours:
                         txt.append('  Connected cells: ' + neighbours_str)
@@ -1148,16 +1148,16 @@ def _run(self, duration, log, log_interval, report_nan, progress, msg):
                     # Show all variables with non-finite values
                     txt.append(
                         'Logged all variables for points: ' + ', '.join(
-                            myokit.strfloat(t, precision=self._precision)
+                            myokit.float.str(t, precision=self._precision)
                             for t in d.time()))
                     txt.append(
-                        'Non-finite valued variables at t=' + myokit.strfloat(
+                        'Non-finite valued variables at t=' + myokit.float.str(
                             d.time()[-2], precision=self._precision))
                     for key, values in d.items():
                         value = values[-2]
                         if not np.isfinite(value):
                             txt.append(
-                                '  ' + str(key) + ' = ' + myokit.strfloat(
+                                '  ' + str(key) + ' = ' + myokit.float.str(
                                     value, precision=self._precision))
 
                     # Show any error in evaluating derivatives
diff --git a/myokit/_sim/psim.py b/myokit/_sim/psim.py
index eff93dbc6..25ea7326f 100644
--- a/myokit/_sim/psim.py
+++ b/myokit/_sim/psim.py
@@ -164,7 +164,7 @@ def __init__(
         # Unique simulation id
         PSimulation._index += 1
         module_name = 'myokit_PSimulation_' + str(PSimulation._index)
-        module_name += '_' + str(myokit._pid_hash())
+        module_name += '_' + str(myokit.pid_hash())
 
         # Arguments
         args = {
diff --git a/myokit/_sim/rhs.c b/myokit/_sim/rhs.c
index 048a91a0a..5e8211107 100644
--- a/myokit/_sim/rhs.c
+++ b/myokit/_sim/rhs.c
@@ -67,7 +67,7 @@ static double engine_time = 0;
 <?
 for var in model.variables(deep=True):
     if var.is_literal():
-        print('static double ' + v(var) + ' = ' + myokit.strfloat(var.rhs().eval()) + ';')
+        print('static double ' + v(var) + ' = ' + myokit.float.str(var.rhs().eval()) + ';')
     else:
         print('static double ' + v(var) + ';')
         if var.is_state():
diff --git a/myokit/_sim/rhs.py b/myokit/_sim/rhs.py
index 95a0c44ec..cd3d103ad 100644
--- a/myokit/_sim/rhs.py
+++ b/myokit/_sim/rhs.py
@@ -56,7 +56,7 @@ def __init__(self, model, variables=None, exclude_selected=False):
         # Extension module id
         RhsBenchmarker._index += 1
         module_name = 'myokit_RhsBenchmarker_' + str(RhsBenchmarker._index)
-        module_name += '_' + str(myokit._pid_hash())
+        module_name += '_' + str(myokit.pid_hash())
 
         # Distutils arguments
         args = {
diff --git a/myokit/_sim/sundials.py b/myokit/_sim/sundials.py
index b50a503b2..ae1d1b668 100644
--- a/myokit/_sim/sundials.py
+++ b/myokit/_sim/sundials.py
@@ -45,7 +45,7 @@ def __init__(self):
 
         # Create Sundials back-end
         mname = 'myokit_sundials_info_' + str(Sundials._index)
-        mname += '_' + str(myokit._pid_hash())
+        mname += '_' + str(myokit.pid_hash())
         fname = os.path.join(myokit.DIR_CFUNC, SOURCE_FILE)
         args = {'module_name': mname}
         try:
diff --git a/myokit/_unit.py b/myokit/_unit.py
index 56edcad01..268ec9b66 100644
--- a/myokit/_unit.py
+++ b/myokit/_unit.py
@@ -186,7 +186,7 @@ def close(unit1, unit2, reltol=1e-9, abstol=1e-9):
         # Compare multipliers in log10 space. This means we can still easily
         # distinguish between e.g. pF and nF, but we lose the ability to tell
         # the difference between e.g. [F (1)] and [F (1.00000000000001)]
-        return myokit._close(unit1._m, unit2._m, reltol, abstol)
+        return myokit.float.close(unit1._m, unit2._m, reltol, abstol)
 
     @staticmethod
     def close_exponent(unit1, unit2, reltol=1e-9, abstol=1e-9):
@@ -202,7 +202,7 @@ def close_exponent(unit1, unit2, reltol=1e-9, abstol=1e-9):
 
         """
         for i, a in enumerate(unit1._x):
-            if not myokit._close(a, unit2._x[i], reltol, abstol):
+            if not myokit.float.close(a, unit2._x[i], reltol, abstol):
                 return False
         return True
 
@@ -319,7 +319,7 @@ def conversion_factor(unit1, unit2, helpers=None):
             raise myokit.IncompatibleUnitError(msg + '.')
 
         # Create Quantity and return
-        fw = myokit._fround(fw)
+        fw = myokit.float.round(fw)
         return Quantity(fw, unit2 / unit1)
 
     @staticmethod
@@ -346,9 +346,9 @@ def __eq__(self, other):
             return False
 
         for i, x in enumerate(self._x):
-            if not myokit._feq(x, other._x[i]):
+            if not myokit.float.eq(x, other._x[i]):
                 return False
-        return myokit._feq(self._m, other._m)
+        return myokit.float.eq(self._m, other._m)
 
     def exponents(self):
         """
@@ -481,7 +481,7 @@ def __pow__(self, f):
         # Evaluates ``self ^ other``
 
         f = float(f)
-        e = [myokit._fround(f * x) for x in self._x]
+        e = [myokit.float.round(f * x) for x in self._x]
         return Unit(e, self._m * f)
 
     def __rdiv__(self, other):  # pragma: no cover    rtruediv used instead
@@ -568,10 +568,13 @@ def __rtruediv__(self, other):
     def _str(self, use_close_for_rounding):
         """
         String representation without preferred representations, that will
-        round to nearby ints if ``use_close_for_rounding=True``.
+        round to nearby integers if ``use_close_for_rounding=True``.
         """
         # Rounding
-        rnd = myokit._cround if use_close_for_rounding else myokit._fround
+        if use_close_for_rounding:
+            rnd = myokit.float.cround
+        else:
+            rnd = myokit.float.round
 
         # SI unit names
         si = ['g', 'm', 's', 'A', 'K', 'cd', 'mol']
@@ -610,7 +613,7 @@ def _str(self, use_close_for_rounding):
 
     def __str__(self):
 
-        # Strategy 1: Try simple look-up (using _feq float comparison)
+        # Strategy 1: Try simple look-up (using float.eq comparison)
         try:
             return '[' + Unit._preferred_representations[self] + ']'
         except KeyError:
@@ -636,10 +639,10 @@ def __str__(self):
 
             # Add multiplier part
             if rep is not None:
-                m = myokit._cround(self._m - m)
+                m = myokit.float.cround(self._m - m)
                 m = 10**m
                 if m >= 1:
-                    m = myokit._cround(m)
+                    m = myokit.float.cround(m)
                 if m < 1e6:
                     m = str(m)
                 else:
diff --git a/myokit/float.py b/myokit/float.py
new file mode 100644
index 000000000..0b3708427
--- /dev/null
+++ b/myokit/float.py
@@ -0,0 +1,107 @@
+#
+# Functions related to floating point numbers. Imported by default.
+#
+# This file is part of Myokit.
+# See http://myokit.org for copyright, sharing, and licensing details.
+#
+from __future__ import absolute_import, division
+from __future__ import print_function, unicode_literals
+
+import sys
+
+import myokit
+
+# Store ref to built-in `round` and `str` functions
+_round = round
+_str = str
+
+
+def close(a, b, reltol=1e-9, abstol=1e-9):
+    """
+    Test whether two numbers are close enough to be considered equal.
+
+    Differs from :meth:`myokit.float.eq` in that it tries to answer the
+    question "are two numbers resulting from various calculations close enough
+    to be considered equal". Whereas :meth`myokit.float.eq` aims to deal with
+    umbers that are numerically indistinguishable but still have a slightly
+    different floating point representation.
+    """
+    # Note the initial == check handles infinity
+    return a == b or abs(a - b) < max(reltol * max(abs(a), abs(b)), abstol)
+
+
+def cround(x, reltol=1e-9, abstol=1e-9):
+    """
+    Checks if a float ``x`` is close to an integer with :meth:`close()`, and if
+    so, returns that integer.
+    """
+    ix = int(_round(x))
+    return ix if close(x, ix, reltol, abstol) else x
+
+
+def eq(a, b):
+    """
+    Checks if floating point numbers ``a`` and ``b`` are equal, or so close to
+    each other that the difference could be a single rounding error.
+    """
+    # Note the initial == check handles infinity
+    return a == b or abs(a - b) < max(abs(a), abs(b)) * sys.float_info.epsilon
+
+
+def geq(a, b):
+    """
+    Checks if ``a >= b``, but using :meth:`myokit.float.eq` instead of ``==``.
+    """
+    return a >= b or abs(a - b) < max(abs(a), abs(b)) * sys.float_info.epsilon
+
+
+def round(x):
+    """
+    Checks if a float ``x`` is within a single rounding error of an integer,
+    and if so, returns that integer.
+    """
+    ix = int(_round(x))
+    return ix if eq(x, ix) else x
+
+
+def str(number, full=False, precision=myokit.DOUBLE_PRECISION):
+    """
+    Converts the given number to a string, using a full-precision
+    representation if needed.
+
+    Arguments:
+
+    ``full``
+        Set to ``True`` to force full precision output.
+    ``precision``
+        Set to ``myokit.SINGLE_PRECISION`` to assume single-precision numbers
+        when formatting with full precision.
+
+    """
+    # Force full precision output
+    if full:
+        if precision == myokit.SINGLE_PRECISION:
+            return myokit.SFSINGLE.format(float(number))
+        else:
+            return myokit.SFDOUBLE.format(float(number))
+
+    # Pass through strings
+    if isinstance(number, _str):
+        return number
+
+    # Handle myokit.Numbers
+    if isinstance(number, myokit.Number):
+        number = number.eval()
+
+    # For most numbers, allow python to format the float
+    s = _str(number)
+    if len(s) < 10:
+        return s
+
+    # But if the number is given with lots of decimals, use the representation
+    # with enough digits to prevent loss of information
+    if precision == myokit.SINGLE_PRECISION:
+        return myokit.SFSINGLE.format(float(number))
+    else:
+        return myokit.SFDOUBLE.format(float(number))
+
diff --git a/myokit/formats/__init__.py b/myokit/formats/__init__.py
index f1b4fabb3..49861e58f 100644
--- a/myokit/formats/__init__.py
+++ b/myokit/formats/__init__.py
@@ -422,7 +422,7 @@ def runnable(self, path, model, protocol=None, *args):
         """
         # Get and test path
         path = os.path.abspath(os.path.expanduser(path))
-        path = myokit.format_path(path)
+        path = myokit.tools.format_path(path)
         self._test_writable_dir(path)
 
         # Clone the model, allowing changes to be made during export
@@ -456,7 +456,7 @@ def runnable(self, path, model, protocol=None, *args):
                 if os.path.exists(file_dir):
                     if not os.path.isdir(file_dir):
                         msg = 'Failed to create directory at: '
-                        msg += myokit.format_path(file_dir)
+                        msg += myokit.tools.format_path(file_dir)
                         msg += ' A file or link with that name already exists.'
                         raise myokit.ExportError(msg)
                 else:   # pragma: no cover
@@ -464,7 +464,7 @@ def runnable(self, path, model, protocol=None, *args):
                         os.makedirs(file_dir)
                     except IOError as e:
                         msg = 'Failed to create directory at: '
-                        msg += myokit.format_path(file_dir)
+                        msg += myokit.tools.format_path(file_dir)
                         msg += ' IOError:' + str(e)
                         raise myokit.ExportError(msg)
 
@@ -472,19 +472,20 @@ def runnable(self, path, model, protocol=None, *args):
             out_name = os.path.join(path, out_name)
             if os.path.exists(out_name):
                 if os.path.isdir(out_name):
-                    msg = 'Directory exists at ' + myokit.format_path(out_name)
+                    msg = 'Directory exists at ' \
+                          + myokit.tools.format_path(out_name)
                     raise myokit.ExportError(msg)
 
             # Check template file
             tpl_name = os.path.join(tpl_dir, tpl_name)
             if not os.path.exists(tpl_name):    # pragma: no cover
                 # Cover pragma: If this happens it's a bug in the exporter
-                msg = 'File not found: ' + myokit.format_path(tpl_name)
+                msg = 'File not found: ' + myokit.tools.format_path(tpl_name)
                 raise myokit.ExportError(msg)
             if not os.path.isfile(tpl_name):    # pragma: no cover
                 # Cover pragma: If this happens it's a bug in the exporter
                 msg = 'Directory found, expecting file at '
-                msg += myokit.format_path(tpl_name)
+                msg += myokit.tools.format_path(tpl_name)
                 raise myokit.ExportError(msg)
 
             # Render
@@ -496,7 +497,7 @@ def runnable(self, path, model, protocol=None, *args):
                 except Exception as e:      # pragma: no cover
                     warnings.warn(
                         'An error ocurred while processing the template at '
-                        + myokit.format_path(tpl_name))
+                        + myokit.tools.format_path(tpl_name))
                     warnings.warn(traceback.format_exc())
                     if isinstance(e, myokit.pype.PypeError):
                         # Pype error? Then add any error details
diff --git a/myokit/formats/ansic/template/euler.c b/myokit/formats/ansic/template/euler.c
index 765536e69..2b915efad 100644
--- a/myokit/formats/ansic/template/euler.c
+++ b/myokit/formats/ansic/template/euler.c
@@ -220,7 +220,7 @@ int main()
     Real state[N_STATE];
 <?
 for var in model.states():
-    print(tab + v(var) + ' = ' + myokit.strfloat(var.state_value()) + ';')
+    print(tab + v(var) + ' = ' + myokit.float.str(var.state_value()) + ';')
 
 ?>
 
diff --git a/myokit/formats/axon/_atf.py b/myokit/formats/axon/_atf.py
index 76b72a399..d786afb39 100644
--- a/myokit/formats/axon/_atf.py
+++ b/myokit/formats/axon/_atf.py
@@ -317,6 +317,6 @@ def save_atf(log, filename, fields=None):
         # Write data
         for i in range(nd):
             f.write((
-                delim.join([myokit.strfloat(next(d)) for d in data]) + eol
+                delim.join([myokit.float.str(next(d)) for d in data]) + eol
             ).encode(_ENC))
 
diff --git a/myokit/formats/cellml/v1/_api.py b/myokit/formats/cellml/v1/_api.py
index 2ac9f6260..b189fd2cf 100644
--- a/myokit/formats/cellml/v1/_api.py
+++ b/myokit/formats/cellml/v1/_api.py
@@ -96,11 +96,11 @@ def create_unit_name(unit):
 
         # Use e-notation if multiple of 10
         multiplier10 = unit.multiplier_log_10()
-        if myokit._feq(multiplier10, int(multiplier10)):
+        if myokit.float.eq(multiplier10, int(multiplier10)):
             multiplier = '1e' + str(int(multiplier10))
 
         # Format as integer
-        elif myokit._feq(multiplier, int(multiplier)):
+        elif myokit.float.eq(multiplier, int(multiplier)):
             multiplier = str(int(multiplier))
 
         # Format as float
diff --git a/myokit/formats/cellml/v1/_writer.py b/myokit/formats/cellml/v1/_writer.py
index 6ca1fed80..7caab379c 100644
--- a/myokit/formats/cellml/v1/_writer.py
+++ b/myokit/formats/cellml/v1/_writer.py
@@ -335,10 +335,10 @@ def _units(self, parent, units):
         # Add multiplier or prefix to first row
         multiplier = myokit_unit.multiplier()
         if multiplier != 1:
-            if myokit._feq(multiplier, int(multiplier)):
+            if myokit.float.eq(multiplier, int(multiplier)):
                 rows[0].attrib['multiplier'] = str(int(multiplier))
             else:
-                rows[0].attrib['multiplier'] = myokit.strfloat(
+                rows[0].attrib['multiplier'] = myokit.float.str(
                     multiplier).strip()
 
     def _variable(self, parent, variable):
@@ -362,7 +362,7 @@ def _variable(self, parent, variable):
 
         # Add initial value
         if variable.initial_value() is not None:
-            element.attrib['initial_value'] = myokit.strfloat(
+            element.attrib['initial_value'] = myokit.float.str(
                 variable.initial_value()).strip()
 
         # Add cmeta id
diff --git a/myokit/formats/cellml/v2/_api.py b/myokit/formats/cellml/v2/_api.py
index 9b43e8e3c..4ec82841e 100644
--- a/myokit/formats/cellml/v2/_api.py
+++ b/myokit/formats/cellml/v2/_api.py
@@ -122,11 +122,11 @@ def create_unit_name(unit):
 
         # Use e-notation if multiple of 10
         multiplier10 = unit.multiplier_log_10()
-        if myokit._feq(multiplier10, int(multiplier10)):
+        if myokit.float.eq(multiplier10, int(multiplier10)):
             multiplier = '1e' + str(int(multiplier10))
 
         # Format as integer
-        elif myokit._feq(multiplier, int(multiplier)):
+        elif myokit.float.eq(multiplier, int(multiplier)):
             multiplier = str(int(multiplier))
 
         # Format as float
diff --git a/myokit/formats/cellml/v2/_writer.py b/myokit/formats/cellml/v2/_writer.py
index b4c39a191..954dd973b 100644
--- a/myokit/formats/cellml/v2/_writer.py
+++ b/myokit/formats/cellml/v2/_writer.py
@@ -257,10 +257,10 @@ def _units(self, parent, units):
         # Add multiplier or prefix to first row
         multiplier = myokit_unit.multiplier()
         if multiplier != 1:
-            if myokit._feq(multiplier, int(multiplier)):
+            if myokit.float.eq(multiplier, int(multiplier)):
                 rows[0].attrib['multiplier'] = str(int(multiplier))
             else:
-                rows[0].attrib['multiplier'] = myokit.strfloat(
+                rows[0].attrib['multiplier'] = myokit.float.str(
                     multiplier).strip()
 
     def _variable(self, parent, variable):
@@ -282,7 +282,7 @@ def _variable(self, parent, variable):
 
         # Add initial value
         if variable is variable.initial_value_variable():
-            element.attrib['initial_value'] = myokit.strfloat(
+            element.attrib['initial_value'] = myokit.float.str(
                 variable.initial_value()).strip()
 
     def write(self, model):
diff --git a/myokit/formats/cuda/_ewriter.py b/myokit/formats/cuda/_ewriter.py
index 85243d1ee..943c1dce6 100644
--- a/myokit/formats/cuda/_ewriter.py
+++ b/myokit/formats/cuda/_ewriter.py
@@ -25,7 +25,7 @@ def __init__(self, precision=myokit.SINGLE_PRECISION):
     #def _ex_derivative(self, e):
 
     def _ex_number(self, e):
-        return myokit.strfloat(e) + 'f' if self._sp else myokit.strfloat(e)
+        return myokit.float.str(e) + 'f' if self._sp else myokit.float.str(e)
 
     #def _ex_prefix_plus(self, e):
     #def _ex_prefix_minus(self, e):
diff --git a/myokit/formats/cuda/template/kernel.cu b/myokit/formats/cuda/template/kernel.cu
index ff6e9adc5..fc004983e 100644
--- a/myokit/formats/cuda/template/kernel.cu
+++ b/myokit/formats/cuda/template/kernel.cu
@@ -60,7 +60,7 @@ model.reserve_unique_names(
     #'time',
     #'pace',
     'I_diff',
-    'dt',    
+    'dt',
     'parameters',
     'state',
     )
@@ -109,7 +109,7 @@ for comp, clist in comp_in.items():
     for bound in bound_variables:
         lhs = bound.lhs()
         if lhs in clist:
-            continue        
+            continue
         for var in comp.variables(deep=True):
             if var.rhs().depends_on(lhs):
                 clist.append(lhs)
@@ -234,7 +234,7 @@ for comp, ilist in comp_in.items():
 ////////////////////////////////////////////////////////////////////////////////
 //! Compute an Euler step of the model.
 ////////////////////////////////////////////////////////////////////////////////
-__device__ int iterate_euler_cu(const Real dt, Real *state, Real I_diff, 
+__device__ int iterate_euler_cu(const Real dt, Real *state, Real I_diff,
                                     Real *parameters)
 {
 <?
@@ -273,7 +273,7 @@ if last_component:
             pre += 'Real '
         if var not in bound_variables:
             print(pre + w.eq(eq) + ';')
-            
+
 ?>
     /* Perform update */
 <?
@@ -291,16 +291,16 @@ for var in model.states():
 
 /* Set the standard initial conditions */
 int get_default_initial_state(Real *state)
-{  
+{
     if (state == 0) return(-1);
-  
+
 <?
 for var in model.states():
     if 'desc' in var.meta:
         print(tab + '// ' + var.meta['desc'])
-    print(tab + v(var) + ' = ' + myokit.strfloat(var.state_value()) + ';')
+    print(tab + v(var) + ' = ' + myokit.float.str(var.state_value()) + ';')
 ?>
-  
+
     return(0);
 }
 
diff --git a/myokit/formats/easyml/_ewriter.py b/myokit/formats/easyml/_ewriter.py
index 1661b1c5b..9689936d3 100644
--- a/myokit/formats/easyml/_ewriter.py
+++ b/myokit/formats/easyml/_ewriter.py
@@ -26,7 +26,7 @@ def __init__(self):
     #def _ex_derivative(self, e):
 
     def _ex_number(self, e):
-        return myokit.strfloat(e)
+        return myokit.float.str(e)
 
     #def _ex_prefix_plus(self, e):
     #def _ex_prefix_minus(self, e):
diff --git a/myokit/formats/easyml/_exporter.py b/myokit/formats/easyml/_exporter.py
index 73df87e34..eb1f493cc 100644
--- a/myokit/formats/easyml/_exporter.py
+++ b/myokit/formats/easyml/_exporter.py
@@ -372,7 +372,7 @@ def lhs(e):
 
             # Write initial conditions
             for v in model.states():
-                f.write(lhs(v) + '_init = ' + myokit.strfloat(v.state_value())
+                f.write(lhs(v) + '_init = ' + myokit.float.str(v.state_value())
                         + eos)
             f.write(eol)
 
diff --git a/myokit/formats/latex/_ewriter.py b/myokit/formats/latex/_ewriter.py
index 49be2c58f..489c92377 100644
--- a/myokit/formats/latex/_ewriter.py
+++ b/myokit/formats/latex/_ewriter.py
@@ -118,7 +118,7 @@ def _ex_derivative(self, e, b):
         b.append(self._flhs(e))
 
     def _ex_number(self, e, b):
-        b.append(myokit.strfloat(e))
+        b.append(myokit.float.str(e))
 
     def _ex_prefix_plus(self, e, b):
         self._ex(e[0], b)
diff --git a/myokit/formats/mathml/_ewriter.py b/myokit/formats/mathml/_ewriter.py
index 06f4f85ac..cf701e4e6 100644
--- a/myokit/formats/mathml/_ewriter.py
+++ b/myokit/formats/mathml/_ewriter.py
@@ -41,7 +41,7 @@ def flhs(lhs):
         self._flhs = flhs
 
         # Default number conversion function
-        self._fnum = lambda x: myokit.strfloat(x.eval())
+        self._fnum = lambda x: myokit.float.str(x.eval())
 
         # Default time variable
         self._tvar = myokit.Name('time')
diff --git a/myokit/formats/opencl/_ewriter.py b/myokit/formats/opencl/_ewriter.py
index 9076dd9ee..391f4c322 100644
--- a/myokit/formats/opencl/_ewriter.py
+++ b/myokit/formats/opencl/_ewriter.py
@@ -48,7 +48,7 @@ def _ex_infix_logical(self, e, op):
     #def _ex_derivative(self, e):
 
     def _ex_number(self, e):
-        return myokit.strfloat(e) + 'f' if self._sp else myokit.strfloat(e)
+        return myokit.float.str(e) + 'f' if self._sp else myokit.float.str(e)
 
     #def _ex_prefix_plus(self, e):
     #def _ex_prefix_minus(self, e):
diff --git a/myokit/formats/opencl/template/cable.c b/myokit/formats/opencl/template/cable.c
index 2d6301ab4..17ea546bf 100644
--- a/myokit/formats/opencl/template/cable.c
+++ b/myokit/formats/opencl/template/cable.c
@@ -75,7 +75,7 @@ if precision == myokit.SINGLE_PRECISION:
 ?>
 /*
  * Prints an opencl error
- */ 
+ */
 void ocl_print_error(const int e, const char* s)
 {
     switch(e) {
@@ -120,7 +120,7 @@ cl_device_id get_device_id()
 	int i;
 	cl_int err;
 	/* Get preferred type */
-	for (i=0; i<platforms_n; i++) {		
+	for (i=0; i<platforms_n; i++) {
 	    err = clGetDeviceIDs(platforms[i], CL_DEVICE_TYPE_GPU, 100, devices, &devices_n);
 	    if(err == CL_SUCCESS) {
 	        return devices[0];
@@ -135,7 +135,7 @@ cl_device_id get_device_id()
 	        return devices[0];
         } else if (err != CL_DEVICE_NOT_FOUND) {
             ocl_check("clGetDeviceIDs", err);
-        }        
+        }
 	}
 	printf("OpenCL Error: Unable to detect an opencl device.\n");
 	abort();
@@ -144,7 +144,7 @@ cl_device_id get_device_id()
 /*
  * Rounds up to the nearest multiple of ws_size.
  */
-static int round_total_size(const int ws_size, const int total_size) 
+static int round_total_size(const int ws_size, const int total_size)
 {
     int size = (total_size / ws_size) * ws_size;
     if(size < total_size) size += ws_size;
@@ -162,7 +162,7 @@ char* read_kernel(const char* fname)
         fprintf(stderr, "Error: Unable to read kernel file.\n");
         abort();
     }
-    
+
     /* Get length */
     fseek(fp, 0L, SEEK_END);
     long n_chars = ftell(fp);
@@ -176,12 +176,12 @@ char* read_kernel(const char* fname)
         fprintf(stderr, "Error: can't allocate enough memory to read file.");
         abort();
     }
-    
+
     /* Read bytes */
     fread(text, n_chars, 1, fp);
-    
+
     /* Close file and return */
-    fclose(fp);    
+    fclose(fp);
     return text;
 }
 
@@ -194,10 +194,10 @@ void set_initial_values(int n_cells, Real* s)
     for(i=0; i<n_cells; i++) {
 <?
 for v in model.states():
-    print(tab*2 + '*s = ' + myokit.strfloat(v.state_value()) + '; s++; // ' + str(v.qname()))
+    print(tab*2 + '*s = ' + myokit.float.str(v.state_value()) + '; s++; // ' + str(v.qname()))
 ?>
     }
-}   
+}
 
 /*
  * Runs a simulation
@@ -209,7 +209,7 @@ for v in model.states():
  * g            The cell-to-cell conductance (for example 1.5)
  */
 int run(const int n_cells, const Real time_start, const Real time_end, const Real dt, const Real g)
-{   
+{
     /* Steps per log action */
     const int steps_per_log = (int)(1.0 / dt);
 
@@ -217,10 +217,10 @@ int run(const int n_cells, const Real time_start, const Real time_end, const Rea
     Real time = time_start;
 
     /* Work group size (in this case, one-dimensional) */
-    size_t local_work_size  = 64; 
+    size_t local_work_size  = 64;
     /* Total number of work items rounded up to a multiple of the local size */
     size_t global_work_size = round_total_size(local_work_size, n_cells);
-    
+
     /* Objects needing destruction */
     cl_context context = NULL;
     cl_command_queue command_queue = NULL;
@@ -231,7 +231,7 @@ int run(const int n_cells, const Real time_start, const Real time_end, const Rea
     cl_mem mbuf_idiff = NULL;
     Real *state = NULL;
     Real *idiff = NULL;
-    
+
     /* Get device id */
     cl_device_id device_id = get_device_id();
     #ifdef __MYOKIT_DEBUG
@@ -255,12 +255,12 @@ int run(const int n_cells, const Real time_start, const Real time_end, const Rea
     int dsize_state = n_cells * n_state * sizeof(Real);
     state = (Real*)malloc(dsize_state);
     set_initial_values(n_cells, state);
-    
+
     /* Create diffusion current vector and copy to device */
     int dsize_idiff = n_cells * sizeof(Real);
     idiff = (Real*)malloc(dsize_idiff);
     for(i=0; i<n_cells; i++) idiff[i] = 0.0;
-        
+
     /* Create memory buffers on the device */
     mbuf_state = clCreateBuffer(context, CL_MEM_READ_WRITE, dsize_state, NULL, &err);
     ocl_check("clCreateBuffer mbuf_state", err);
@@ -268,7 +268,7 @@ int run(const int n_cells, const Real time_start, const Real time_end, const Rea
     ocl_check("clCreateBuffer mbuf_idiff", err);
     #ifdef __MYOKIT_DEBUG
     printf("Created buffers.\n");
-    #endif  
+    #endif
 
     /* Copy data into buffers */
     ocl(clEnqueueWriteBuffer(command_queue, mbuf_state, CL_TRUE, 0, dsize_state, state, 0, NULL, NULL));
@@ -276,7 +276,7 @@ int run(const int n_cells, const Real time_start, const Real time_end, const Rea
     #ifdef __MYOKIT_DEBUG
     printf("Set initial state.\n");
     #endif
-    
+
     /* Load and compile the kernel program(s) */
     char* source = read_kernel("kernel.cl");
     program = clCreateProgramWithSource(context, 1, (const char**)&source, NULL, &err);
@@ -285,7 +285,7 @@ int run(const int n_cells, const Real time_start, const Real time_end, const Rea
     if (err == CL_BUILD_PROGRAM_FAILURE) {
         fprintf(stderr, "OpenCL Error: Kernel failed to compile.\n");
         fprintf(stderr, "----------------------------------------");
-        fprintf(stderr, "---------------------------------------\n");        
+        fprintf(stderr, "---------------------------------------\n");
         /* Extract build log */
         size_t blog_size;
         clGetProgramBuildInfo(program, device_id, CL_PROGRAM_BUILD_LOG, 0, NULL, &blog_size);
@@ -300,7 +300,7 @@ int run(const int n_cells, const Real time_start, const Real time_end, const Rea
     #ifdef __MYOKIT_DEBUG
     printf("Program created and built.\n");
     #endif
-    
+
     /* Create the kernels */
     kernel_cell = clCreateKernel(program, "cell_step", &err);
     ocl_check("clCreateKernel", err);
@@ -316,14 +316,14 @@ int run(const int n_cells, const Real time_start, const Real time_end, const Rea
     ocl(clSetKernelArg(kernel_cell, 2, sizeof(dt), &dt));
     ocl(clSetKernelArg(kernel_cell, 3, sizeof(mbuf_state), &mbuf_state));
     ocl(clSetKernelArg(kernel_cell, 4, sizeof(mbuf_idiff), &mbuf_idiff));
-    ocl(clSetKernelArg(kernel_diff, 0, sizeof(n_cells), &n_cells));    
+    ocl(clSetKernelArg(kernel_diff, 0, sizeof(n_cells), &n_cells));
     ocl(clSetKernelArg(kernel_diff, 1, sizeof(g), &g));
     ocl(clSetKernelArg(kernel_diff, 2, sizeof(mbuf_state), &mbuf_state));
     ocl(clSetKernelArg(kernel_diff, 3, sizeof(mbuf_idiff), &mbuf_idiff));
     #ifdef __MYOKIT_DEBUG
     printf("Arguments passed into kernels.\n");
     #endif
-    
+
     /* Add log header */
     printf("engine.time");
     for(i=0; i<n_cells; i++) {
@@ -333,7 +333,7 @@ for v in model.states():
 ?>
     }
     printf("\n");
-    
+
     /* Start simulation */
     int steps = 0;
     int steps_till_log = 1;
@@ -342,7 +342,7 @@ for v in model.states():
         ocl(clSetKernelArg(kernel_cell, 1, sizeof(time), &time));
         ocl(clEnqueueNDRangeKernel(command_queue, kernel_diff, 1, NULL, &global_work_size, &local_work_size, 0, NULL, NULL));
         ocl(clEnqueueNDRangeKernel(command_queue, kernel_cell, 1, NULL, &global_work_size, &local_work_size, 0, NULL, NULL));
-    
+
         /* Show results for first cell */
         if (--steps_till_log == 0) {
             ocl(clEnqueueReadBuffer(command_queue, mbuf_state, CL_TRUE, 0, dsize_state, state, 0, NULL, NULL));
@@ -353,16 +353,16 @@ for v in model.states():
             printf("\n");
             steps_till_log = steps_per_log;
         }
-        
+
         steps++;
         time = time_start + steps * dt;
     }
-        
+
 done:
     /* Tidy up */
     #ifdef __MYOKIT_DEBUG
     printf("Tidying up.\n");
-    #endif  
+    #endif
 
     ocl(clFlush(command_queue));
     ocl(clFinish(command_queue));
@@ -370,15 +370,15 @@ for v in model.states():
     ocl(clReleaseMemObject(mbuf_idiff));
     ocl(clReleaseKernel(kernel_cell));
     ocl(clReleaseKernel(kernel_diff));
-    ocl(clReleaseProgram(program)); 
-    ocl(clReleaseCommandQueue(command_queue));    
+    ocl(clReleaseProgram(program));
+    ocl(clReleaseCommandQueue(command_queue));
     ocl(clReleaseContext(context));
     free(state);
     free(idiff);
-    
+
     #ifdef __MYOKIT_DEBUG
     printf("Done.\n");
-    #endif 
+    #endif
 }
 
 int main()
@@ -390,6 +390,6 @@ int main()
         0.005,      // dt
         1.5         // g
         );
-        
+
     return 0;
 }
diff --git a/myokit/formats/opencl/template/kernel.cl b/myokit/formats/opencl/template/kernel.cl
index ee3a631db..fe01b49f5 100644
--- a/myokit/formats/opencl/template/kernel.cl
+++ b/myokit/formats/opencl/template/kernel.cl
@@ -54,7 +54,7 @@ for comp, clist in comp_in.items():
     for bound in bound_variables:
         lhs = bound.lhs()
         if lhs in clist:
-            continue        
+            continue
         for var in comp.variables(deep=True):
             if var.rhs().depends_on(lhs):
                 clist.append(lhs)
@@ -69,7 +69,7 @@ def set_pointers(names=None):
     """
     Tells the expression writer to write the given variable names (given as
     LhsExpression objects) as pointers.
-    
+
     Calling set_pointers a second time clears the first list. Calling with
     ``names=None`` unsets all pointers.
     """
@@ -84,13 +84,13 @@ def v(var):
     """
     if isinstance(var, myokit.Derivative):
         # Explicitly asked for derivative
-        pre = '*' if var in ptrs else ''        
+        pre = '*' if var in ptrs else ''
         return pre + 'D_' + var.var().uname()
     if isinstance(var, myokit.Name):
         var = var.var()
     if var in bound_variables:
         return bound_variables[var]
-    pre = '*' if myokit.Name(var) in ptrs else ''        
+    pre = '*' if myokit.Name(var) in ptrs else ''
     return pre + var.uname()
 w.set_lhs_function(v)
 
@@ -190,10 +190,10 @@ __kernel void cell_step(const int n_cells, const Real time, const Real dt, __glo
     const int cid = get_global_id(0);
     const int offset = cid * n_state;
     if(cid >= n_cells) return;
-    
+
     /* Diffusion */
     Real idiff = idiff_vec[cid];
-    
+
     /* Pacing */
     Real pace;
     if(cid < 4) {
@@ -201,7 +201,7 @@ __kernel void cell_step(const int n_cells, const Real time, const Real dt, __glo
     } else {
         pace = 0;
     }
-    
+
 <?
 print(tab + '/* Evaluate derivatives */')
 for comp in comp_order:
@@ -238,7 +238,7 @@ if last_component:
             pre += 'Real '
         if var not in bound_variables:
             print(pre + w.eq(eq) + ';')
-            
+
 ?>
     /* Perform update */
 <?
diff --git a/myokit/formats/opencl/template/minilog.py b/myokit/formats/opencl/template/minilog.py
index fab9a28ee..25eadf915 100644
--- a/myokit/formats/opencl/template/minilog.py
+++ b/myokit/formats/opencl/template/minilog.py
@@ -237,14 +237,14 @@ def save_csv(self, filename, pad=None):
                 for i in range(0, n):
                     line = []
                     for d in data:
-                        line.append(myokit.strfloat(next(d)))
+                        line.append(myokit.float.str(next(d)))
                     f.write(delim.join(line) + eol)
             else:
                 for i in range(0, n):
                     line = []
                     for d in data:
                         try:
-                            line.append(myokit.strfloat(next(d)))
+                            line.append(myokit.float.str(next(d)))
                         except StopIteration:
                             line.append('0')
                     f.write(delim.join(line) + eol)
diff --git a/myokit/formats/python/_ewriter.py b/myokit/formats/python/_ewriter.py
index 3eba0e719..42c936e7c 100644
--- a/myokit/formats/python/_ewriter.py
+++ b/myokit/formats/python/_ewriter.py
@@ -73,7 +73,7 @@ def _ex_partial_derivative(self, e):
         return self._flhs(e)
 
     def _ex_number(self, e):
-        return myokit.strfloat(e)
+        return myokit.float.str(e)
 
     def _ex_prefix_plus(self, e):
         return self.ex(e[0])
diff --git a/myokit/lib/hh.py b/myokit/lib/hh.py
index ce187a2bc..a621816e8 100644
--- a/myokit/lib/hh.py
+++ b/myokit/lib/hh.py
@@ -410,7 +410,8 @@ def from_component(
             parameters = [
                 x for x in component.variables(const=True) if x.is_literal()]
             # Sort by qname, using natural sort
-            parameters.sort(key=lambda x: myokit._natural_sort_key(x.qname()))
+            parameters.sort(
+                key=lambda x: myokit.tools.natural_sort_key(x.qname()))
 
         # Get current
         if current is None:
diff --git a/myokit/lib/markov.py b/myokit/lib/markov.py
index 4b3bbf8df..e09cf5557 100644
--- a/myokit/lib/markov.py
+++ b/myokit/lib/markov.py
@@ -481,7 +481,8 @@ def from_component(
             parameters = [
                 x for x in component.variables(const=True) if x.is_literal()]
             # Sort by qname, using natural sort
-            parameters.sort(key=lambda x: myokit._natural_sort_key(x.qname()))
+            parameters.sort(
+                key=lambda x: myokit.tools.natural_sort_key(x.qname()))
 
         # Get current
         if current is None:
@@ -1848,7 +1849,7 @@ def find_markov_models(model):
 
         # Get sorted list of states for output
         states = list(states)
-        states.sort(key=lambda x: myokit._natural_sort_key(x.qname()))
+        states.sort(key=lambda x: myokit.tools.natural_sort_key(x.qname()))
 
         # Check all members of `group` are a linear combination of states
         try:
diff --git a/myokit/pacing.py b/myokit/pacing.py
index 2035ac436..25fac11d9 100644
--- a/myokit/pacing.py
+++ b/myokit/pacing.py
@@ -16,7 +16,7 @@
 
 def blocktrain(period, duration, offset=0, level=1.0, limit=0):
     """
-    Creates an train of block pulses.
+    Creates a train of block pulses.
 
     Each pulse lasts ``duration`` time units and a pulse is initiated every
     ``period`` time units.
diff --git a/myokit/tools.py b/myokit/tools.py
new file mode 100644
index 000000000..53620e17f
--- /dev/null
+++ b/myokit/tools.py
@@ -0,0 +1,480 @@
+#
+# Myokit utility functions. This module gathers functions that are used in
+# Myokit, but are not particularly Myokit-dependent, i.e. they could easily be
+# stand-alone. Imported by default.
+#
+# This file is part of Myokit.
+# See http://myokit.org for copyright, sharing, and licensing details.
+#
+from __future__ import absolute_import, division
+from __future__ import print_function, unicode_literals
+
+import os
+import re
+import shutil
+import stat
+import sys
+import tempfile
+import threading
+import timeit
+
+# StringIO in Python 2 and 3
+try:
+    from cStringIO import StringIO
+except ImportError:  # pragma: no python 2 cover
+    from io import StringIO
+
+
+# Natural sort regex
+_natural_sort_regex = re.compile('([0-9]+)')
+
+
+class Benchmarker(object):
+    """
+    Allows benchmarking using the with statement.
+
+    Example::
+
+        m,p,x = myokit.load('example')
+        s = myokit.Simulation(m, p)
+        b = myokit.tools.Benchmarker()
+        s.run()
+        print(b.time())
+        b.reset()
+        s.run()
+        print(b.time())
+
+    """
+    def __init__(self):
+        self._start = timeit.default_timer()
+
+    def format(self, time):
+        """
+        Formats a (non-integer) number of seconds, returns a string like
+        "5 weeks, 3 days, 1 hour, 4 minutes, 9 seconds", or "0.0019 seconds".
+        """
+        if time < 60:
+            return '1 second' if time == 1 else str(time) + ' seconds'
+        output = []
+        time = int(round(time))
+        units = [
+            (604800, 'week'),
+            (86400, 'day'),
+            (3600, 'hour'),
+            (60, 'minute'),
+        ]
+        for k, name in units:
+            f = time // k
+            if f > 0 or output:
+                output.append(str(f) + ' ' + (name if f == 1 else name + 's'))
+            time -= f * k
+        output.append('1 second' if time == 1 else str(time) + ' seconds')
+        return ', '.join(output)
+
+    def reset(self):
+        """
+        Resets this timer's start time.
+        """
+        self._start = timeit.default_timer()
+
+    def time(self):
+        """
+        Returns the time since benchmarking started.
+        """
+        return timeit.default_timer() - self._start
+
+
+class capture(object):
+    """
+    Context manager that temporarily redirects the current standard output and
+    error streams, and captures anything that's written to them.
+
+    Example::
+
+        with myokit.tools.capture() as a:
+            print('This will be captured')
+
+    Within a single thread, captures can be nested, for example::
+
+        with myokit.tools.capture() as a:
+            print('This will be captured by a')
+            with myokit.tools.capture() as b:
+                print('This will be captured by b, not a')
+            print('This will be captured by a again')
+
+    Capturing is thread-safe: a lock is used to ensure only a single thread is
+    capturing at any time. For example, if we have a function::
+
+        def f(i):
+            with myokit.tools.capture() as a:
+                print(i)
+                ...
+
+            return a.text()
+
+    and this is called from several threads, the ``capture`` acts as a lock (a
+    ``threading.RLock``) so that one thread will need to finish executing the
+    code within the ``with`` statement before a second thread can start
+    capturing.
+
+    In multiprocessing, no locks are used, and no memory or streams are shared
+    so that this should also be safe.
+
+    By default, this method works by simply redirecting ``sys.stdout`` and
+    ``sys.stderr``. This captures any output written by the Python interpreter,
+    but does not catch output from C/C++ extensions or subrocesses. To also
+    catch that output start the capture with the optional argument
+    ``fd=True``, which enables a file descriptor duplication method of
+    redirection.
+    """
+    # Note: It seems we need to capture both streams to make the file
+    # descriptor method work, and we want both anyway throughout Myokit, so
+    # there are no options to choose stdout or stderr here.
+
+    # Lock to stop other threads from capturing while this thread is capturing.
+    _rlock = threading.RLock()
+
+    def __init__(self, fd=False):
+
+        # Are we already capturing? This is needed in case someone enters the
+        # same context twice.
+        self._active_count = 0
+
+        # Captured text
+        self._txt_out = None
+        self._txt_err = None
+
+        # Capturing mode
+        self._fd = bool(fd)
+
+        # Shared by both modes
+        self._org_out = None     # Original stdout object
+        self._org_err = None     # Original stderr object
+
+        # Python stream redirects only
+        self._tmp_out = None    # Temporary stdout StringIO object
+        self._tmp_err = None    # Temporary stderr StringIO object
+
+        # File descriptor redirects only
+        self._out_fd = None      # File descriptor used for output
+        self._err_fd = None      # File descriptor used for errors
+        self._org_out_fd = None  # Back-up of original stdout file descriptor
+        self._org_err_fd = None  # Back-up of original stderr file descriptor
+        self._file_out = None    # Temporary file to write output to
+        self._file_err = None    # Temporary file to write errors to
+
+    def __enter__(self):
+        """Called when the context is entered."""
+        # Avoid entering the same context object twice
+        self._active_count += 1
+        if self._active_count == 1:
+
+            # Wait until other threads have stopped capturing
+            capture._rlock.acquire()
+
+            # Set up redirection
+            self._start()
+
+        # Return context
+        return self
+
+    def __exit__(self, exc_type, exc_value, traceback):
+        """Called when exiting the context."""
+        self._active_count -= 1
+        if self._active_count == 0:
+            self._stop()
+            capture._rlock.release()
+
+    def _start(self):
+        """Starts capturing output to stdout and stderr."""
+
+        # Clear any cached text
+        self._txt_out = None
+        self._txt_err = None
+
+        # Save the current output / error streams
+        self._org_out = sys.stdout
+        self._org_err = sys.stderr
+
+        # Redirect
+        if not self._fd:
+            # Create temporary output and error streams
+            self._tmp_out = StringIO()
+            self._tmp_err = StringIO()
+
+            # Redirect, attempting to flush first
+            try:
+                sys.stdout.flush()
+            except AttributeError:  # pragma: no cover
+                pass
+            finally:
+                sys.stdout = self._tmp_out
+
+            try:
+                sys.stderr.flush()
+            except AttributeError:  # pragma: no cover
+                pass
+            finally:
+                sys.stderr = self._tmp_err
+
+        else:
+            # File-descriptor redirection
+
+            # Get file descriptors used for output and errors.
+            #
+            # These represent open files, and are the low-level equivalent of
+            # sys.stdout and sys.stderr: all subprocesses etc. use these
+            # descriptors to write output to.
+            #
+            # Note that we get these file descriptors from __stdout__ and
+            # __stderr__, to get the original streams that subprocesses etc.
+            # will use, even if someone has redirected the Python-level
+            # streams.
+            #
+            # On https://docs.python.org/3/library/sys.html#module-sys, it says
+            # that stdout/err as well as __stdout__ can be None (e.g. in spyder
+            # on windows).  In other cases (pythonw.exe) they can be set but
+            # return a negative file descriptor (indicating it's invalid).
+            # So here we check if __stdout__ is None and if so set a negative
+            # fileno so that we can catch both cases at once in the rest of the
+            # code.
+            #
+            if sys.__stdout__ is not None:
+                self._out_fd = sys.__stdout__.fileno()
+            else:   # pragma: no cover
+                self._out_fd = -1
+            if sys.__stderr__ is not None:
+                self._err_fd = sys.__stderr__.fileno()
+            else:   # pragma: no cover
+                self._err_fd = -1
+
+            # We start by making a copy of these file descriptors using dup(),
+            # meaning that two file numbers will point to the same file. In a
+            # later step we use dup2() to write a different value into the
+            # descriptor, so that all output gets redirected. When exiting, we
+            # will use dup2() again to copy the original value back in.
+            #
+            # We check if the fd >= 0 first, because the previous code could
+            # get a negative (invalid) value either when __stdxxx__ is None or
+            # when a negative value is returned.
+            #
+            # On windows, the order is important: First dup both stdout and
+            # stderr, then dup2 the new descriptors in. This prevents a weird
+            # infinite recursion on windows ipython / python shell.
+            #
+            if self._out_fd >= 0:
+                self._org_out_fd = os.dup(self._out_fd)
+            if self._err_fd >= 0:
+                self._org_err_fd = os.dup(self._err_fd)
+
+            # Create temporary files to redirect the output to. Make sure they
+            # aren't opened in binary mode, and specify + for reading and
+            # writing.
+            self._file_out = tempfile.TemporaryFile(mode='w+')
+            self._file_err = tempfile.TemporaryFile(mode='w+')
+
+            # Apply Python-level redirection of sys.stdxxx in addition to the
+            # lower-level redirect of sys.__stdxxx___.
+            try:
+                sys.stdout.flush()
+            except AttributeError:  # pragma: no cover
+                pass
+            finally:
+                sys.stdout = self._file_out
+            try:
+                sys.stderr.flush()
+            except AttributeError:  # pragma: no cover
+                pass
+            finally:
+                sys.stderr = self._file_err
+
+            # If possible, overwrite the stdout and stderr file descriptors
+            # with the file descriptor for our files. C/C++ code etc. will
+            # still look in this location for a descriptor, so that output is
+            # redirected.
+            if self._out_fd >= 0:
+                sys.__stdout__.flush()
+                os.dup2(self._file_out.fileno(), self._out_fd)
+            if self._err_fd >= 0:
+                sys.__stderr__.flush()
+                os.dup2(self._file_err.fileno(), self._err_fd)
+
+    def _stop(self):
+        """Stops capturing output."""
+
+        # Flush any remaining output
+        sys.stdout.flush()
+        sys.stderr.flush()
+
+        # Undo redirect
+        if not self._fd:
+            # Direct back at original
+            sys.stdout = self._org_out
+            sys.stderr = self._org_err
+
+            # Store text
+            self._txt_out = self._tmp_out.getvalue()
+            self._txt_err = self._tmp_err.getvalue()
+
+            # Tidy
+            self._tmp_out = self._tmp_err = None
+
+        else:
+            # Undo redirection
+            sys.stdout = self._org_out
+            if self._org_out_fd is not None:
+                # Copy our backup of the original fd into out_fd
+                os.dup2(self._org_out_fd, self._out_fd)
+                # And close our backup
+                os.close(self._org_out_fd)
+
+            sys.stderr = self._org_err
+            if self._org_err_fd is not None:
+                os.dup2(self._org_err_fd, self._err_fd)
+                os.close(self._org_err_fd)
+
+            # Close temporary files and store capture output
+            try:
+                self._file_out.seek(0)
+                self._txt_out = self._file_out.read()
+                self._file_out.close()
+            except ValueError:  # pragma: no cover
+                # In rare cases, I've seen a ValueError, "underlying buffer has
+                # been detached".
+                pass
+            try:
+                self._file_err.seek(0)
+                self._txt_err = self._file_err.read()
+                self._file_err.close()
+            except ValueError:  # pragma: no cover
+                pass
+
+            # Tidy
+            self._out_fd = self._err_fd = None
+            self._org_out_fd = self._org_err_fd = None
+            self._file_out = self._file_err = None
+
+        # Unset reference to original streams (for proper decreffing etc.)
+        self._stdout = self._stderr = None
+
+    def err(self):
+        """
+        Returns the text captured from stderr, or an empty string if nothing
+        was captured or capturing is still active.
+        """
+        if self._txt_err is None:
+            return ''
+        text = self._txt_err
+
+        # In Python 2, the text needs to be decoded from ascii
+        if sys.hexversion < 0x03000000:  # pragma: no python 3 cover
+            text = text.decode('ascii', 'ignore')
+
+        return text
+
+    def out(self):
+        """
+        Returns the text captured from stdout, or an empty string if nothing
+        was captured or capturing is still active.
+        """
+        if self._txt_out is None:
+            return ''
+        text = self._txt_out
+
+        # In Python 2, the text needs to be decoded from ascii
+        if sys.hexversion < 0x03000000:  # pragma: no python 3 cover
+            text = text.decode('ascii', 'ignore')
+
+        return text
+
+    def text(self):
+        """
+        Returns the combined text captured from output and error text, if any
+        (output first, then error text).
+        """
+        return self.out() + self.err()
+
+
+def format_path(path, root='.'):
+    """
+    Formats a path for use in user messages. If the given path is a
+    subdirectory of the current directory this part is chopped off.
+
+    Alternatively, a ``root`` directory may be given explicitly: any
+    subdirectory of this path will be formatted relative to ``root``.
+
+    This function differs from os.path.relpath() in the way it handles paths
+    *outside* the root: In these cases relpath returns a relative path such as
+    '../../' while this function returns an absolute path.
+    """
+    if path == '':
+        path = '.'
+    try:
+        path = os.path.relpath(path, root)
+    except ValueError:  # pragma: no cover
+        # This can happen on windows, if `path` is on a different drive than
+        # root (so that no relative path from one to the other can be made).
+        return path
+
+    if '..' in path:
+        path = os.path.abspath(os.path.join(root, path))
+    return path
+
+
+def lvsd(s, t):
+    """Returns the Levenshtein distance between two strings ``s`` and ``t``."""
+    # This implementation is adapted from wikipedia:
+    # en.wikipedia.org/wiki/Levenshtein_distance#Iterative_with_full_matrix
+    if s == t:
+        return 0
+    if len(s) == 0:
+        return len(t)
+    if len(t) == 0:
+        return len(s)
+
+    n = len(t) + 1
+    v0 = list(range(n))
+    v1 = [None] * n
+    for i, si in enumerate(s):
+        v1[0] = i + 1
+        for j, tj in enumerate(t):
+            v1[j + 1] = min(
+                v0[j + 1] + 1, v1[j] + 1, v0[j] if si == tj else v0[j] + 1)
+        v0, v1 = v1, v0
+    return v0[n - 1]
+
+
+def natural_sort_key(s):
+    """
+    Function to use as ``key`` in a sort, to get natural sorting of strings
+    (e.g. "2" before "10").
+
+    Usage examples::
+
+        names.sort(key=myokit.tools.natural_sort_key)
+
+        variables.sort(key=lambda v: myokit.tools.natural_sort_key(v.qname()))
+
+    """
+    # Code adapted from: http://stackoverflow.com/questions/4836710/
+    return [
+        int(text) if text.isdigit() else text.lower()
+        for text in _natural_sort_regex.split(s)]
+
+
+def rmtree(path):
+    """
+    Version of ``shutil.rmtree`` that handles Windows "access denied" errors
+    (when the user is lacking write permissions, but is allowed to set them).
+    """
+    # From https://stackoverflow.com/questions/2656322
+    def onerror(function, path, excinfo):   # pragma: no cover
+        if not os.access(path, os.W_OK):
+            # Give user write permissions (remove read-only flag)
+            os.chmod(path, stat.S_IWUSR)
+            function(path)
+        else:
+            raise
+
+    shutil.rmtree(path, ignore_errors=False, onerror=onerror)
+


Test patch:
diff --git a/myokit/tests/data/cv1d.mmt b/myokit/tests/data/cv1d.mmt
index 71eea6b55..d9b59f7f6 100644
--- a/myokit/tests/data/cv1d.mmt
+++ b/myokit/tests/data/cv1d.mmt
@@ -5,7 +5,7 @@ Implementation of the Luo-Rudy model I for the ventricular myocyte.
 Based on an updated formulation downloaded from http://rudylab.wustl.edu
 
 Original copyright notice:
-- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
+- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 Copyright (C) 2006 by Leonid Livshitz and Yoram Rudy
 Email rudy@wustl.edu
 
@@ -66,8 +66,8 @@ ENa = cell.RTF * log(cell.Na_o / cell.Na_i)
 a = 1 - 1 / (1 + exp(-(V + 40) / 0.24))
     desc: Used instead of if statement.
 dot(m) =  alpha * (1 - m) - beta * m
-    alpha = if(V == -47.13, -3.2, 
-            0.32 * (V + 47.13) / (1 - exp(-0.1 * (V + 47.13))))    
+    alpha = if(V == -47.13, -3.2,
+            0.32 * (V + 47.13) / (1 - exp(-0.1 * (V + 47.13))))
     beta  = 0.08 * exp(-V / 11)
     desc: """
     Activation parameter of the fast sodium channel, based on Ebihara-Johnson
@@ -185,5 +185,5 @@ d = s.run(500, log=['engine.time', 'membrane.V', 'engine.pace'])
 b = d.block1d()
 b.save('cv1d.zip')
 
-print(myokit.strfloat(b.cv('membrane.V')))
+print(myokit.float.str(b.cv('membrane.V')))
 
diff --git a/myokit/tests/data/io/block2d.mmt b/myokit/tests/data/io/block2d.mmt
index 63870414c..96287ac67 100644
--- a/myokit/tests/data/io/block2d.mmt
+++ b/myokit/tests/data/io/block2d.mmt
@@ -5,7 +5,7 @@ Implementation of the Luo-Rudy model I for the ventricular myocyte.
 Based on an updated formulation downloaded from http://rudylab.wustl.edu
 
 Original copyright notice:
-- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
+- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 Copyright (C) 2006 by Leonid Livshitz and Yoram Rudy
 Email rudy@wustl.edu
 
@@ -66,8 +66,8 @@ ENa = cell.RTF * log(cell.Na_o / cell.Na_i)
 a = 1 - 1 / (1 + exp(-(V + 40) / 0.24))
     desc: Used instead of if statement.
 dot(m) =  alpha * (1 - m) - beta * m
-    alpha = if(V == -47.13, -3.2, 
-            0.32 * (V + 47.13) / (1 - exp(-0.1 * (V + 47.13))))    
+    alpha = if(V == -47.13, -3.2,
+            0.32 * (V + 47.13) / (1 - exp(-0.1 * (V + 47.13))))
     beta  = 0.08 * exp(-V / 11)
     desc: """
     Activation parameter of the fast sodium channel, based on Ebihara-Johnson
@@ -185,5 +185,3 @@ d = s.run(500, log=['engine.time', 'membrane.V', 'engine.pace'])
 b = d.block2d()
 b.save('block2d.zip')
 
-#print(myokit.strfloat(b.cv('membrane.V')))
-
diff --git a/myokit/tests/shared.py b/myokit/tests/shared.py
index 38e187eb4..5451df45f 100644
--- a/myokit/tests/shared.py
+++ b/myokit/tests/shared.py
@@ -80,7 +80,7 @@ def path(self, path=None):
 
     def __exit__(self, exc_type, exc_value, traceback):
         try:
-            myokit._rmtree(self._dir)
+            myokit.tools.rmtree(self._dir)
         finally:
             self._dir = None
 
diff --git a/myokit/tests/test_aux.py b/myokit/tests/test_aux.py
index b65eecbe0..0aaae52d8 100755
--- a/myokit/tests/test_aux.py
+++ b/myokit/tests/test_aux.py
@@ -9,7 +9,6 @@
 from __future__ import print_function, unicode_literals
 
 import os
-import sys
 import unittest
 
 # StringIO in Python 2 and 3
@@ -26,7 +25,7 @@
 
 import myokit
 
-from shared import DIR_DATA
+from shared import DIR_DATA, WarningCollector
 
 
 class AuxTest(unittest.TestCase):
@@ -35,41 +34,14 @@ class AuxTest(unittest.TestCase):
     """
 
     def test_benchmarker(self):
-        # Test the benchmarker.
-
-        b = myokit.Benchmarker()
-        x = [0] * 1000
-        t0 = b.time()
-        self.assertTrue(t0 >= 0)
-        x = [0] * 1000
-        t1 = b.time()
-        self.assertTrue(t1 >= t0)
-        x = [0] * 1000
-        t2 = b.time()
-        self.assertTrue(t2 >= t1)
-        for i in range(1000):
-            x = [0] * 1000
-        t3 = b.time()
-        self.assertTrue(t3 >= t2)
-        b.reset()
-        t4 = b.time()
-        self.assertTrue(t4 < t3)
-
-        self.assertEqual(b.format(1), '1 second')
-        self.assertEqual(b.format(61), '1 minute, 1 second')
-        self.assertEqual(b.format(60), '1 minute, 0 seconds')
-        self.assertEqual(b.format(180), '3 minutes, 0 seconds')
-        self.assertEqual(b.format(3600), '1 hour, 0 minutes, 0 seconds')
-        self.assertEqual(b.format(3661), '1 hour, 1 minute, 1 second')
-        self.assertEqual(
-            b.format(3600 * 24), '1 day, 0 hours, 0 minutes, 0 seconds')
-        self.assertEqual(
-            b.format(3600 * 24 * 7),
-            '1 week, 0 days, 0 hours, 0 minutes, 0 seconds')
+        # Deprecated alias of myokit.tools.Benchmarker
+        with WarningCollector() as c:
+            b = myokit.Benchmarker()
+        self.assertIn('`myokit.Benchmarker` is deprecated', c.text())
+        self.assertIsInstance(b, myokit.tools.Benchmarker)
 
     def test_date(self):
         # Test date formatting method.
-
         import time
         for i in range(3):
             a = time.strftime(myokit.DATE_FORMAT)
@@ -149,130 +121,26 @@ def test_default_script(self):
 
         # TODO: Run with tiny model?
 
-    def test_examplify(self):
-        # Test examplify.
-
-        self.assertEqual(myokit._aux._examplify('test.txt'), 'test.txt')
-        self.assertEqual(myokit._aux._examplify('example'), myokit.EXAMPLE)
-
-    def test_float_functions(self):
-        # Test the floating point comparison methods
-
-        # Test that test is going to work
-        x = 49
-        y = 1 / (1 / x)
-        self.assertNotEqual(x, y)
-
-        # Test if feq allows 1 floating point error
-        # Test if
-        self.assertTrue(myokit._feq(x, y))
-        self.assertTrue(myokit._fgeq(x, y))
-        self.assertTrue(myokit._fgeq(y, x))
-        x += x * sys.float_info.epsilon
-        self.assertTrue(myokit._feq(x, y))
-        self.assertTrue(myokit._fgeq(x, y))
-        self.assertTrue(myokit._fgeq(y, x))
-        x += x * sys.float_info.epsilon
-        self.assertFalse(myokit._feq(x, y))
-        self.assertTrue(myokit._fgeq(x, y))
-        self.assertFalse(myokit._fgeq(y, x))
-
-        # Test rounding
-        self.assertNotEqual(49, y)
-        self.assertEqual(49, myokit._fround(y))
-        self.assertNotEqual(49, myokit._fround(x))
-        self.assertEqual(0.5, myokit._fround(0.5))
-        self.assertIsInstance(myokit._fround(y), int)
-
-        # Try with negative numbers
-        self.assertNotEqual(-49, -y)
-        self.assertEqual(-49, myokit._fround(-y))
-        self.assertNotEqual(-49, myokit._fround(-x))
-        self.assertEqual(-0.5, myokit._fround(-0.5))
-
-        # Test that _close allows bigger errors
-        x = 49
-        y = x * (1 + 1e-11)
-        self.assertNotEqual(x, y)
-        self.assertFalse(myokit._feq(x, y))
-        self.assertTrue(myokit._close(x, y))
-
-        # And that close thinks everything small is equal
-        x = 1e-16
-        y = 1e-12
-        self.assertNotEqual(x, y)
-        self.assertFalse(myokit._feq(x, y))
-        self.assertTrue(myokit._close(x, y))
-
-        # Test rounding based on closeness
-        x = 49
-        y = x * (1 + 1e-11)
-        self.assertNotEqual(x, y)
-        self.assertEqual(x, myokit._cround(y))
-        self.assertIsInstance(myokit._cround(y), int)
-        self.assertNotEqual(x, myokit._cround(49.001))
-
     def test_format_float_dict(self):
-        # Test myokit.format_float_dict.
-
+        # Test myokit.format_float_dict, which is deprecated
         d = {'one': 1, 'Definitely two': 2, 'Three-ish': 3.1234567}
-        x = myokit.format_float_dict(d).splitlines()
+        with WarningCollector() as c:
+            x = myokit.format_float_dict(d).splitlines()
+        self.assertIn('`myokit.format_float_dict` is deprecated', c.text())
         self.assertEqual(len(x), 3)
         self.assertEqual(x[0], 'Definitely two = 2')
         self.assertEqual(x[1], 'Three-ish      = 3.1234567')
         self.assertEqual(x[2], 'one            = 1')
 
     def test_format_path(self):
-        # Test format_path().
-
-        # Normal use
-        self.assertEqual(
-            myokit.format_path(os.path.join('a', 'b', 'c')),
-            os.path.join('a', 'b', 'c'))
-
-        # No trailing slash
-        self.assertEqual(
-            myokit.format_path('a'), 'a')
-        self.assertEqual(
-            myokit.format_path('a/b/'), os.path.join('a', 'b'))
-
-        # Use with custom root
+        # Deprecated alias of myokit.tools.format_path
         root = os.path.join(os.path.abspath('.'), 'a')
-        self.assertEqual(
-            myokit.format_path(os.path.join(root, 'b', 'c'), root),
-            os.path.join('b', 'c'))
-
-        # Empty path
-        self.assertEqual(
-            myokit.format_path(''), '.')
-        self.assertEqual(
-            myokit.format_path('.'), '.')
-
-        # Filesystem root
-        self.assertEqual(
-            myokit.format_path('/'), os.path.abspath('/'))
-        self.assertEqual(
-            myokit.format_path('/', root='/'), '.')
-
-        # Path outside of root
-        self.assertEqual(
-            myokit.format_path(
-                os.path.abspath('test'),
-                os.path.abspath('test/tost')),
-            os.path.abspath('test'))
-
-    def test_levenshtein_distance(self):
-        # Test the levenshtein distance method.
-
-        self.assertEqual(myokit._lvsd('kitten', 'sitting'), 3)
-        self.assertEqual(myokit._lvsd('sitting', 'kitten'), 3)
-        self.assertEqual(myokit._lvsd('saturday', 'sunday'), 3)
-        self.assertEqual(myokit._lvsd('sunday', 'saturday'), 3)
-        self.assertEqual(myokit._lvsd('michael', 'jennifer'), 7)
-        self.assertEqual(myokit._lvsd('jennifer', 'michael'), 7)
-        self.assertEqual(myokit._lvsd('jennifer', ''), 8)
-        self.assertEqual(myokit._lvsd('', 'jennifer'), 8)
-        self.assertEqual(myokit._lvsd('', ''), 0)
+        a = os.path.join(root, 'b', 'c')
+        b = os.path.join(os.path.abspath('.'), 'a')
+        with WarningCollector() as c:
+            x = myokit.format_path(a, b)
+        self.assertIn('`myokit.format_path` is deprecated', c.text())
+        self.assertEqual(x, myokit.tools.format_path(a, b))
 
     def test_model_comparison(self):
         # Test the model comparison class.
@@ -282,7 +150,7 @@ def test_model_comparison(self):
         m1 = myokit.load_model(m1)
         m2 = myokit.load_model(m2)
 
-        with myokit.capture() as capture:
+        with myokit.tools.capture() as capture:
             c = myokit.ModelComparison(m1, m2, live=True)
 
         differences = [
@@ -359,17 +227,6 @@ def test_model_comparison(self):
         self.assertEqual(len(c), 2)
         self.assertEqual(len(c), len(myokit.ModelComparison(m3, m1)))
 
-    def test_natural_sort_key(self):
-        # Test natural sort key method.
-
-        a = ['a12', 'a3', 'a11', 'a2', 'a10', 'a1']
-        b = ['a1', 'a2', 'a3', 'a10', 'a11', 'a12']
-        self.assertNotEqual(a, b)
-        a.sort()
-        self.assertNotEqual(a, b)
-        a.sort(key=lambda x: myokit._natural_sort_key(x))
-        self.assertEqual(a, b)
-
     def test_numpy_writer(self):
         # Test NumPy expression writer obtaining method.
 
@@ -421,16 +278,16 @@ def test_run(self):
             's = myokit.Simulation(m, p)',
             's.run(200)',
         ])
-        with myokit.capture():
+        with myokit.tools.capture():
             myokit.run(m, p, x)
-        with myokit.capture():
+        with myokit.tools.capture():
             myokit.run(m, p, '[[script]]\n' + x)
         self.assertRaises(ZeroDivisionError, myokit.run, m, p, 'print(1 / 0)')
 
         # Test with stringio
         x = "print('Hi there')"
         s = StringIO()
-        with myokit.capture() as c:
+        with myokit.tools.capture() as c:
             myokit.run(m, p, x, stderr=s, stdout=s)
         self.assertEqual(c.text(), '')
         self.assertEqual(s.getvalue(), 'Hi there\n')
@@ -655,66 +512,18 @@ def test_step(self):
             self.assertEqual(line, x[i])
         self.assertEqual(len(x), len(y))
 
-    def test_myokit_strfloat(self):
-        # Test float to string conversion.
-
-        # String should be passed through
-        # Note: convert to str() to test in python 2 and 3.
-        self.assertEqual(myokit.strfloat(str('123')), '123')
-
-        # Simple numbers
-        self.assertEqual(myokit.strfloat(0), '0')
-        self.assertEqual(myokit.strfloat(0.0000), '0.0')
-        self.assertEqual(myokit.strfloat(1.234), '1.234')
-        self.assertEqual(
-            myokit.strfloat(0.12432656245e12), ' 1.24326562450000000e+11')
-        self.assertEqual(myokit.strfloat(-0), '0')
-        self.assertEqual(myokit.strfloat(-0.0000), '-0.0')
-        self.assertEqual(myokit.strfloat(-1.234), '-1.234')
-        self.assertEqual(
-            myokit.strfloat(-0.12432656245e12), '-1.24326562450000000e+11')
-
-        # Strings are not converted
-        x = '1.234'
-        self.assertEqual(x, myokit.strfloat(x))
-
-        # Myokit Numbers are converted
-        x = myokit.Number(1.23)
-        self.assertEqual(myokit.strfloat(x), '1.23')
-
-        # Single and double precision
-        self.assertEqual(
-            myokit.strfloat(-1.234, precision=myokit.SINGLE_PRECISION),
-            '-1.234')
-        self.assertEqual(
-            myokit.strfloat(
-                -0.124326562458734682153498731245756e12,
-                precision=myokit.SINGLE_PRECISION),
-            '-1.243265625e+11')
-        self.assertEqual(
-            myokit.strfloat(-1.234, precision=myokit.DOUBLE_PRECISION),
-            '-1.234')
-        self.assertEqual(
-            myokit.strfloat(
-                -0.124326562458734682153498731245756e12,
-                precision=myokit.DOUBLE_PRECISION),
-            '-1.24326562458734680e+11')
-
-        # Full precision override
-        self.assertEqual(
-            myokit.strfloat(1.23, True), ' 1.22999999999999998e+00')
-        self.assertEqual(
-            myokit.strfloat(1.23, True, myokit.DOUBLE_PRECISION),
-            ' 1.22999999999999998e+00')
-        self.assertEqual(
-            myokit.strfloat(1.23, True, myokit.SINGLE_PRECISION),
-            ' 1.230000000e+00')
+    def test_strfloat(self):
+        # Deprecated alias of myokit.float.str
+        args = ['-1.234', True, myokit.SINGLE_PRECISION]
+        with WarningCollector() as c:
+            x = myokit.strfloat(*args)
+        self.assertIn('`myokit.strfloat` is deprecated', c.text())
+        self.assertEqual(x, myokit.float.str(*args))
 
     def test_time(self):
         # Test time formatting method.
-
         import time
-        for i in range(3):
+        for i in range(6):
             a = time.strftime(myokit.TIME_FORMAT)
             b = myokit.time()
             if a == b:
diff --git a/myokit/tests/test_float.py b/myokit/tests/test_float.py
new file mode 100755
index 000000000..a21dbb9bd
--- /dev/null
+++ b/myokit/tests/test_float.py
@@ -0,0 +1,147 @@
+#!/usr/bin/env python3
+#
+# Tests the methods and classes in myokit.floats
+#
+# This file is part of Myokit.
+# See http://myokit.org for copyright, sharing, and licensing details.
+#
+from __future__ import absolute_import, division
+from __future__ import print_function, unicode_literals
+
+import sys
+import unittest
+
+import myokit
+
+
+class FloatTest(unittest.TestCase):
+    """Tests the methods in ``myokit.float``."""
+
+    def test_close(self):
+        # Test that close allows bigger errors than eq
+        x = 49
+        y = x * (1 + 1e-11)
+        self.assertNotEqual(x, y)
+        self.assertFalse(myokit.float.eq(x, y))
+        self.assertTrue(myokit.float.close(x, y))
+
+        # And that close thinks everything small is equal
+        x = 1e-16
+        y = 1e-12
+        self.assertNotEqual(x, y)
+        self.assertFalse(myokit.float.eq(x, y))
+        self.assertTrue(myokit.float.close(x, y))
+
+    def test_cround(self):
+        # Test rounding based on closeness
+        x = 49
+        y = x * (1 + 1e-11)
+        self.assertNotEqual(x, y)
+        self.assertEqual(x, myokit.float.cround(y))
+        self.assertIsInstance(myokit.float.cround(y), int)
+        self.assertNotEqual(x, myokit.float.cround(49.001))
+
+    def test_eq_and_geq(self):
+        # Test the floating point comparison methods
+
+        # Check that test is going to work
+        x = 49
+        y = 1 / (1 / x)
+        self.assertNotEqual(x, y)
+
+        # Test if eq allows 1 floating point error
+        self.assertTrue(myokit.float.eq(x, y))
+        self.assertTrue(myokit.float.geq(x, y))
+        self.assertTrue(myokit.float.geq(y, x))
+        x += x * sys.float_info.epsilon
+        self.assertTrue(myokit.float.eq(x, y))
+        self.assertTrue(myokit.float.geq(x, y))
+        self.assertTrue(myokit.float.geq(y, x))
+        x += x * sys.float_info.epsilon
+        self.assertFalse(myokit.float.eq(x, y))
+        self.assertTrue(myokit.float.geq(x, y))
+        self.assertFalse(myokit.float.geq(y, x))
+
+    def test_round(self):
+        # Test rounding, ignoring single rounding errors
+
+        # Check that test is going to work
+        x = 49
+        y = 1 / (1 / x)
+        self.assertNotEqual(x, y)
+        x += 2 * x * sys.float_info.epsilon
+        self.assertFalse(myokit.float.eq(x, y))
+        self.assertTrue(myokit.float.geq(x, y))
+        self.assertFalse(myokit.float.geq(y, x))
+
+        # Test rounding
+        self.assertNotEqual(49, y)
+        self.assertEqual(49, myokit.float.round(y))
+        self.assertNotEqual(49, myokit.float.round(x))
+        self.assertEqual(0.5, myokit.float.round(0.5))
+        self.assertIsInstance(myokit.float.round(y), int)
+
+        # Try with negative numbers
+        self.assertNotEqual(-49, -y)
+        self.assertEqual(-49, myokit.float.round(-y))
+        self.assertNotEqual(-49, myokit.float.round(-x))
+        self.assertEqual(-0.5, myokit.float.round(-0.5))
+
+    def test_str(self):
+        # Test float to string conversion.
+
+        # String should be passed through
+        # Note: convert to str() to test in python 2 and 3.
+        self.assertEqual(myokit.float.str(str('123')), '123')
+
+        # Simple numbers
+        self.assertEqual(myokit.float.str(0), '0')
+        self.assertEqual(myokit.float.str(0.0000), '0.0')
+        self.assertEqual(myokit.float.str(1.234), '1.234')
+        self.assertEqual(
+            myokit.float.str(0.12432656245e12), ' 1.24326562450000000e+11')
+        self.assertEqual(myokit.float.str(-0), '0')
+        self.assertEqual(myokit.float.str(-0.0000), '-0.0')
+        self.assertEqual(myokit.float.str(-1.234), '-1.234')
+        self.assertEqual(
+            myokit.float.str(-0.12432656245e12), '-1.24326562450000000e+11')
+
+        # Strings are not converted
+        x = '1.234'
+        self.assertEqual(x, myokit.float.str(x))
+
+        # Myokit Numbers are converted
+        x = myokit.Number(1.23)
+        self.assertEqual(myokit.float.str(x), '1.23')
+
+        # Single and double precision
+        self.assertEqual(
+            myokit.float.str(-1.234, precision=myokit.SINGLE_PRECISION),
+            '-1.234')
+        self.assertEqual(
+            myokit.float.str(
+                -0.124326562458734682153498731245756e12,
+                precision=myokit.SINGLE_PRECISION),
+            '-1.243265625e+11')
+        self.assertEqual(
+            myokit.float.str(-1.234, precision=myokit.DOUBLE_PRECISION),
+            '-1.234')
+        self.assertEqual(
+            myokit.float.str(
+                -0.124326562458734682153498731245756e12,
+                precision=myokit.DOUBLE_PRECISION),
+            '-1.24326562458734680e+11')
+
+        # Full precision override
+        self.assertEqual(
+            myokit.float.str(1.23, True), ' 1.22999999999999998e+00')
+        self.assertEqual(
+            myokit.float.str(1.23, True, myokit.DOUBLE_PRECISION),
+            ' 1.22999999999999998e+00')
+        self.assertEqual(
+            myokit.float.str(1.23, True, myokit.SINGLE_PRECISION),
+            ' 1.230000000e+00')
+
+
+if __name__ == '__main__':
+    unittest.main()
diff --git a/myokit/tests/test_formats_exporters.py b/myokit/tests/test_formats_exporters.py
index 01da03d8b..549c6d42a 100755
--- a/myokit/tests/test_formats_exporters.py
+++ b/myokit/tests/test_formats_exporters.py
@@ -265,7 +265,7 @@ def test_stan_exporter(self):
         for v in m.get('ina').variables(const=True):
             if v.name()[:1] == 'p':
                 parameters.append(v)
-        parameters.sort(key=lambda v: myokit._natural_sort_key(v.name()))
+        parameters.sort(key=lambda v: myokit.tools.natural_sort_key(v.name()))
 
         # Set output
         output = 'ina.INa'
diff --git a/myokit/tests/test_loadsave.py b/myokit/tests/test_io.py
similarity index 96%
rename from myokit/tests/test_loadsave.py
rename to myokit/tests/test_io.py
index 705f83c5a..1732557e5 100755
--- a/myokit/tests/test_loadsave.py
+++ b/myokit/tests/test_io.py
@@ -1,6 +1,7 @@
 #!/usr/bin/env python3
 #
-# Tests the mmt file load/save functionality
+# Tests the loading/saveing mmt files and states.
+# See also test_parsing.py
 #
 # This file is part of Myokit.
 # See http://myokit.org for copyright, sharing, and licensing details.
@@ -23,11 +24,15 @@
     basestring = str
 
 
-class LoadSaveTest(unittest.TestCase):
-    """
-    Tests various parts of load/save behavior for model/protocol/script.
-    """
-    def test_multiline_string_indent(self):
+class LoadSaveMmtTest(unittest.TestCase):
+    """Tests load/save functions for mmt files."""
+
+    def test_examplify(self):
+        # Test _examplify.
+        self.assertEqual(myokit._io._examplify('test.txt'), 'test.txt')
+        self.assertEqual(myokit._io._examplify('example'), myokit.EXAMPLE)
+
+    def test_load_multiline_string_indent(self):
         # Test what happens when you load save a string that gets auto-indented
 
         # Create model with multi-line meta-data property
@@ -81,80 +86,39 @@ def test_multiline_string_indent(self):
             d2 = m2.meta['desc']
             self.assertEqual(d2, dr)
 
-    def test_save_model(self):
-        # Test if the correct parts are saved/loaded from disk using the
-        # ``save_model()`` method.
+    def test_load_partial(self):
+        # Test loading of partial files.
+        mpath = os.path.join(DIR_DATA, 'beeler-1977-model.mmt')
+        ppath = os.path.join(DIR_DATA, 'beeler-1977-protocol.mmt')
+        spath = os.path.join(DIR_DATA, 'beeler-1977-script.mmt')
 
-        ipath = os.path.join(DIR_DATA, 'lr-1991.mmt')
-        # Test example loading
-        m = myokit.load_model('example')
-        self.assertIsInstance(m, myokit.Model)
-        # Test file loading
-        m = myokit.load_model(ipath)
+        m, p, x = myokit.load(mpath)
         self.assertIsInstance(m, myokit.Model)
-        with TemporaryDirectory() as d:
-            opath = d.path('loadsave.mmt')
-            myokit.save_model(opath, m)
-            # Test no other parts were written
-            with open(opath, 'r') as f:
-                text = f.read()
-            self.assertTrue('[[model]]' in text)
-            self.assertFalse('[[protocol]]' in text)
-            self.assertFalse('[[script]]' in text)
-            # Test reloading
-            mm = myokit.load_model(opath)
-            self.assertIsInstance(mm, myokit.Model)
-            self.assertEqual(mm.code(), m.code())
-
-    def test_save_protocol(self):
-        # Test if the correct parts are saved/loaded from disk using the
-        # ``save_protocol()`` method.
+        self.assertIsNone(p)
+        self.assertIsNone(x)
 
-        ipath = os.path.join(DIR_DATA, 'lr-1991.mmt')
-        # Test example loading
-        p = myokit.load_protocol('example')
-        self.assertIsInstance(p, myokit.Protocol)
-        # Test file loading
-        p = myokit.load_protocol(ipath)
+        m, p, x = myokit.load(ppath)
+        self.assertIsNone(m)
         self.assertIsInstance(p, myokit.Protocol)
-        with TemporaryDirectory() as d:
-            opath = d.path('test.mmt')
-            myokit.save_protocol(opath, p)
-            # Test no other parts were written
-            with open(opath, 'r') as f:
-                text = f.read()
-            self.assertFalse('[[model]]' in text)
-            self.assertTrue('[[protocol]]' in text)
-            self.assertFalse('[[script]]' in text)
-            # Test reloading
-            pp = myokit.load_protocol(opath)
-            self.assertIsInstance(pp, myokit.Protocol)
-            self.assertEqual(pp.code(), p.code())
-
-    def test_save_script(self):
-        # Test if the correct parts are saved/loaded from disk using the
-        # ``save_script()`` method.
+        self.assertIsNone(x)
 
-        ipath = os.path.join(DIR_DATA, 'lr-1991.mmt')
-        # Test example loading
-        x = myokit.load_script('example')
-        self.assertTrue(isinstance(x, basestring))
-        # Test file loading
-        x = myokit.load_script(ipath)
+        m, p, x = myokit.load(spath)
+        self.assertIsNone(m)
+        self.assertIsNone(p)
         self.assertTrue(isinstance(x, basestring))
-        with TemporaryDirectory() as d:
-            opath = d.path('test.mmt')
-            myokit.save_script(opath, x)
-            # Test no other parts were written
-            with open(opath, 'r') as f:
-                text = f.read()
-            self.assertFalse('[[model]]' in text)
-            self.assertFalse('[[protocol]]' in text)
-            self.assertTrue('[[script]]' in text)
-            # Test reloading
-            xx = myokit.load_script(opath)
-            self.assertTrue(isinstance(xx, basestring))
-            self.assertEqual(x, xx)
+
+        self.assertRaises(
+            myokit.SectionNotFoundError, myokit.load_model, ppath)
+        self.assertRaises(
+            myokit.SectionNotFoundError, myokit.load_model, spath)
+        self.assertRaises(
+            myokit.SectionNotFoundError, myokit.load_protocol, mpath)
+        self.assertRaises(
+            myokit.SectionNotFoundError, myokit.load_protocol, spath)
+        self.assertRaises(
+            myokit.SectionNotFoundError, myokit.load_script, mpath)
+        self.assertRaises(
+            myokit.SectionNotFoundError, myokit.load_script, ppath)
 
     def test_save(self):
         # Test if the correct parts are saved/loaded from disk using the
@@ -291,39 +255,84 @@ def test_save(self):
                 text = f.read()
             self.assertEqual(text, myokit.save(model=m, protocol=p, script=x))
 
-    def test_load_partial(self):
-        # Test loading of partial files.
-        mpath = os.path.join(DIR_DATA, 'beeler-1977-model.mmt')
-        ppath = os.path.join(DIR_DATA, 'beeler-1977-protocol.mmt')
-        spath = os.path.join(DIR_DATA, 'beeler-1977-script.mmt')
+    def test_save_model(self):
+        # Test if the correct parts are saved/loaded from disk using the
+        # ``save_model()`` method.
 
-        m, p, x = myokit.load(mpath)
+        ipath = os.path.join(DIR_DATA, 'lr-1991.mmt')
+        # Test example loading
+        m = myokit.load_model('example')
         self.assertIsInstance(m, myokit.Model)
-        self.assertIsNone(p)
-        self.assertIsNone(x)
+        # Test file loading
+        m = myokit.load_model(ipath)
+        self.assertIsInstance(m, myokit.Model)
+        with TemporaryDirectory() as d:
+            opath = d.path('loadsave.mmt')
+            myokit.save_model(opath, m)
+            # Test no other parts were written
+            with open(opath, 'r') as f:
+                text = f.read()
+            self.assertTrue('[[model]]' in text)
+            self.assertFalse('[[protocol]]' in text)
+            self.assertFalse('[[script]]' in text)
+            # Test reloading
+            mm = myokit.load_model(opath)
+            self.assertIsInstance(mm, myokit.Model)
+            self.assertEqual(mm.code(), m.code())
 
-        m, p, x = myokit.load(ppath)
-        self.assertIsNone(m)
+    def test_save_protocol(self):
+        # Test if the correct parts are saved/loaded from disk using the
+        # ``save_protocol()`` method.
+
+        ipath = os.path.join(DIR_DATA, 'lr-1991.mmt')
+        # Test example loading
+        p = myokit.load_protocol('example')
         self.assertIsInstance(p, myokit.Protocol)
-        self.assertIsNone(x)
+        # Test file loading
+        p = myokit.load_protocol(ipath)
+        self.assertIsInstance(p, myokit.Protocol)
+        with TemporaryDirectory() as d:
+            opath = d.path('test.mmt')
+            myokit.save_protocol(opath, p)
+            # Test no other parts were written
+            with open(opath, 'r') as f:
+                text = f.read()
+            self.assertFalse('[[model]]' in text)
+            self.assertTrue('[[protocol]]' in text)
+            self.assertFalse('[[script]]' in text)
+            # Test reloading
+            pp = myokit.load_protocol(opath)
+            self.assertIsInstance(pp, myokit.Protocol)
+            self.assertEqual(pp.code(), p.code())
 
-        m, p, x = myokit.load(spath)
-        self.assertIsNone(m)
-        self.assertIsNone(p)
+    def test_save_script(self):
+        # Test if the correct parts are saved/loaded from disk using the
+        # ``save_script()`` method.
+
+        ipath = os.path.join(DIR_DATA, 'lr-1991.mmt')
+        # Test example loading
+        x = myokit.load_script('example')
+        self.assertTrue(isinstance(x, basestring))
+        # Test file loading
+        x = myokit.load_script(ipath)
         self.assertTrue(isinstance(x, basestring))
+        with TemporaryDirectory() as d:
+            opath = d.path('test.mmt')
+            myokit.save_script(opath, x)
+            # Test no other parts were written
+            with open(opath, 'r') as f:
+                text = f.read()
+            self.assertFalse('[[model]]' in text)
+            self.assertFalse('[[protocol]]' in text)
+            self.assertTrue('[[script]]' in text)
+            # Test reloading
+            xx = myokit.load_script(opath)
+            self.assertTrue(isinstance(xx, basestring))
+            self.assertEqual(x, xx)
 
-        self.assertRaises(
-            myokit.SectionNotFoundError, myokit.load_model, ppath)
-        self.assertRaises(
-            myokit.SectionNotFoundError, myokit.load_model, spath)
-        self.assertRaises(
-            myokit.SectionNotFoundError, myokit.load_protocol, mpath)
-        self.assertRaises(
-            myokit.SectionNotFoundError, myokit.load_protocol, spath)
-        self.assertRaises(
-            myokit.SectionNotFoundError, myokit.load_script, mpath)
-        self.assertRaises(
-            myokit.SectionNotFoundError, myokit.load_script, ppath)
+
+class LoadSaveStateTest(unittest.TestCase):
+    """Tests loading and saving states."""
 
     def test_load_save_state(self):
         # Test loading/saving state.
diff --git a/myokit/tests/test_pacing_system_c.py b/myokit/tests/test_pacing_system_c.py
index 5bbe9cbfa..3906a3806 100755
--- a/myokit/tests/test_pacing_system_c.py
+++ b/myokit/tests/test_pacing_system_c.py
@@ -134,7 +134,7 @@ def test_fixed_form_pacing_ansic(self):
             del(t3, t4)
             # Get the pacing value at the points, measure how long it takes
             pacing = AnsicFixedFormPacing(list(t), list(v))
-            b = myokit.Benchmarker()
+            b = myokit.tools.Benchmarker()
             v2 = [pacing.pace(x) for x in t2]
             print(b.time())
             # Plot the new points
diff --git a/myokit/tests/test_parsing.py b/myokit/tests/test_parsing.py
index 34752bc0d..40ca0e83d 100755
--- a/myokit/tests/test_parsing.py
+++ b/myokit/tests/test_parsing.py
@@ -1,6 +1,7 @@
 #!/usr/bin/env python3
 #
 # Tests the parsing module.
+# See also test_io.py
 #
 # This file is part of Myokit.
 # See http://myokit.org for copyright, sharing, and licensing details.
diff --git a/myokit/tests/test_progress_reporters.py b/myokit/tests/test_progress_reporters.py
index f3ad778d8..55098ee66 100755
--- a/myokit/tests/test_progress_reporters.py
+++ b/myokit/tests/test_progress_reporters.py
@@ -23,7 +23,7 @@ def test_progress_printer(self):
         # Test basic functionality.
 
         # Zero call
-        with myokit.capture() as c:
+        with myokit.tools.capture() as c:
             p = myokit.ProgressPrinter()
             self.assertTrue(p.update(0))
         pattern1 = re.compile(
@@ -34,7 +34,7 @@ def test_progress_printer(self):
             self.assertTrue(pattern1.match(line))
 
         # Normal call
-        with myokit.capture() as c:
+        with myokit.tools.capture() as c:
             p = myokit.ProgressPrinter()
             self.assertTrue(p.update(0.49))
         pattern2 = re.compile(
@@ -46,7 +46,7 @@ def test_progress_printer(self):
             self.assertTrue(pattern2.match(line))
 
         # Call that will take minutes
-        with myokit.capture() as c:
+        with myokit.tools.capture() as c:
             p = myokit.ProgressPrinter()
             time.sleep(0.1)
             self.assertTrue(p.update(1e-3))
@@ -56,7 +56,7 @@ def test_progress_printer(self):
         # will print to stdout (which won't have been redirected yet).
 
         # Status every ten percent
-        with myokit.capture() as c:
+        with myokit.tools.capture() as c:
             p = myokit.ProgressPrinter(digits=-1)
             self.assertTrue(p.update(0))
             self.assertTrue(p.update(0.08))
diff --git a/myokit/tests/test_pype.py b/myokit/tests/test_pype.py
index ea36d4a07..43b2dbf4d 100755
--- a/myokit/tests/test_pype.py
+++ b/myokit/tests/test_pype.py
@@ -34,7 +34,7 @@ def test_process_errors(self):
         self.assertRaisesRegex(
             ValueError, 'dict', e.process, 'file.txt', [])
 
-        with myokit.capture():
+        with myokit.tools.capture():
 
             # Test not-a-file
             self.assertRaises(IOError, e.process, 'file.txt', {})
diff --git a/myokit/tests/test_simulation_1d.py b/myokit/tests/test_simulation_1d.py
index 2b7552a57..c57def906 100755
--- a/myokit/tests/test_simulation_1d.py
+++ b/myokit/tests/test_simulation_1d.py
@@ -110,7 +110,7 @@ def test_with_progress_reporter(self):
         # Test using a progress reporter
         s = myokit.Simulation1d(m, p, ncells=5)
         s.set_step_size(0.05)
-        with myokit.capture() as c:
+        with myokit.tools.capture() as c:
             s.run(110, progress=myokit.ProgressPrinter())
         c = c.text().splitlines()
         self.assertTrue(len(c) > 0)
diff --git a/myokit/tests/test_simulation_cvode.py b/myokit/tests/test_simulation_cvode.py
index da874cd9a..c17a3f3bb 100755
--- a/myokit/tests/test_simulation_cvode.py
+++ b/myokit/tests/test_simulation_cvode.py
@@ -15,11 +15,10 @@
 import re
 import unittest
 import sys
-import warnings
 
 import myokit
 
-from shared import DIR_DATA, CancellingReporter
+from shared import DIR_DATA, CancellingReporter, WarningCollector
 
 # Unit testing in Python 2 and 3
 try:
@@ -214,7 +213,7 @@ def test_progress_reporter(self):
 
         # Test if it works
         sim = myokit.LegacySimulation(self.model, self.protocol)
-        with myokit.capture() as c:
+        with myokit.tools.capture() as c:
             sim.run(110, progress=myokit.ProgressPrinter())
         c = c.text().splitlines()
         self.assertEqual(len(c), 2)
@@ -423,7 +422,7 @@ def test_simulation_error_2(self):
         v = m.get('membrane.V')
         v.set_rhs(myokit.Multiply(v.rhs(), myokit.Number(1e18)))
         s = myokit.LegacySimulation(m, self.protocol)
-        with warnings.catch_warnings(record=True):
+        with WarningCollector():
             self.assertRaisesRegex(
                 myokit.SimulationError, 'numerical error', s.run, 5000)
 
diff --git a/myokit/tests/test_simulation_cvodes.py b/myokit/tests/test_simulation_cvodes.py
index 4b17189c9..5ca54561f 100755
--- a/myokit/tests/test_simulation_cvodes.py
+++ b/myokit/tests/test_simulation_cvodes.py
@@ -15,11 +15,15 @@
 import re
 import unittest
 import sys
-import warnings
 
 import myokit
 
-from shared import DIR_DATA, CancellingReporter, test_case_pk_model
+from shared import (
+    CancellingReporter,
+    DIR_DATA,
+    test_case_pk_model,
+    WarningCollector,
+)
 
 # Unit testing in Python 2 and 3
 try:
@@ -213,7 +217,7 @@ def test_progress_reporter(self):
 
         # Test if it works
         sim = myokit.Simulation(self.model, self.protocol)
-        with myokit.capture() as c:
+        with myokit.tools.capture() as c:
             sim.run(110, progress=myokit.ProgressPrinter())
         c = c.text().splitlines()
         self.assertEqual(len(c), 2)
@@ -574,7 +578,7 @@ def test_simulation_error_2(self):
         v = m.get('membrane.V')
         v.set_rhs(myokit.Multiply(v.rhs(), myokit.Number(1e18)))
         s = myokit.Simulation(m, self.protocol)
-        with warnings.catch_warnings(record=True):
+        with WarningCollector():
             self.assertRaisesRegex(
                 myokit.SimulationError, 'numerical error', s.run, 5000)
 
diff --git a/myokit/tests/test_simulation_ic.py b/myokit/tests/test_simulation_ic.py
index eba80b71a..41e3aceea 100755
--- a/myokit/tests/test_simulation_ic.py
+++ b/myokit/tests/test_simulation_ic.py
@@ -82,7 +82,7 @@ def test_progress_reporter(self):
         # Test using a progress reporter
         with WarningCollector() as c:
             s = myokit.ICSimulation(m, p)
-        with myokit.capture() as c:
+        with myokit.tools.capture() as c:
             s.run(110, progress=myokit.ProgressPrinter())
         c = c.text().splitlines()
         self.assertTrue(len(c) > 0)
diff --git a/myokit/tests/test_simulation_opencl.py b/myokit/tests/test_simulation_opencl.py
index 1aaa8ff80..2c501d823 100755
--- a/myokit/tests/test_simulation_opencl.py
+++ b/myokit/tests/test_simulation_opencl.py
@@ -787,7 +787,7 @@ def test_basic(self):
             'ica.ICa',
         ]
         # Run simulation
-        with myokit.capture():
+        with myokit.tools.capture():
             logf, logt = s.run(run, logf=logf, logt=logt, log_interval=0.01)
 
         self.assertEqual(len(logf), 1 + 2 * nfx * nfy)
@@ -836,7 +836,7 @@ def test_against_cvode(self):
         # Run simulation
         tmax = 100
         dlog = 0.1
-        with myokit.capture():
+        with myokit.tools.capture():
             d1, logt = s1.run(tmax, logf=logvars, logt=logt, log_interval=dlog)
         del(logt)
         d1 = d1.npview()
diff --git a/myokit/tests/test_simulation_p.py b/myokit/tests/test_simulation_p.py
index f9e74f37f..744921134 100755
--- a/myokit/tests/test_simulation_p.py
+++ b/myokit/tests/test_simulation_p.py
@@ -203,7 +203,7 @@ def test_progress_reporter(self):
         with WarningCollector() as c:
             s = myokit.PSimulation(
                 m, p, variables=['membrane.V'], parameters=['ina.gNa'])
-        with myokit.capture() as c:
+        with myokit.tools.capture() as c:
             s.run(2, progress=myokit.ProgressPrinter())
         c = c.text().splitlines()
         self.assertEqual(len(c), 2)
diff --git a/myokit/tests/test_system_info.py b/myokit/tests/test_system_info.py
index f791dc91a..cdd1250da 100755
--- a/myokit/tests/test_system_info.py
+++ b/myokit/tests/test_system_info.py
@@ -30,7 +30,7 @@ def test_system_info(self):
         matplotlib.use('template')
 
         self.assertIsInstance(myokit.system(), basestring)
-        with myokit.capture():
+        with myokit.tools.capture():
             myokit.system(live_printing=True)
 
 
diff --git a/myokit/tests/test_aux_capture.py b/myokit/tests/test_tools.py
similarity index 55%
rename from myokit/tests/test_aux_capture.py
rename to myokit/tests/test_tools.py
index b2372dbf6..c108aaa58 100755
--- a/myokit/tests/test_aux_capture.py
+++ b/myokit/tests/test_tools.py
@@ -1,6 +1,6 @@
 #!/usr/bin/env python3
 #
-# Tests the capture methods in myokit._aux
+# Tests the methods and classes in myokit.tools
 #
 # This file is part of Myokit.
 # See http://myokit.org for copyright, sharing, and licensing details.
@@ -8,6 +8,7 @@
 from __future__ import absolute_import, division
 from __future__ import print_function, unicode_literals
 
+import os
 import sys
 import threading
 import time
@@ -22,18 +23,53 @@
 import myokit
 
 
-class AuxCaptureTest(unittest.TestCase):
-    """
-    Test the capture methods in myokit.aux
-    """
+class BenchMarkerTest(unittest.TestCase):
+    """Tests the ``Benchmarker``."""
+
+    def test_benchmarker(self):
+        # Test the benchmarker.
+
+        b = myokit.tools.Benchmarker()
+        x = [0] * 1000
+        t0 = b.time()
+        self.assertTrue(t0 >= 0)
+        x = [0] * 1000
+        t1 = b.time()
+        self.assertTrue(t1 >= t0)
+        x = [0] * 1000
+        t2 = b.time()
+        self.assertTrue(t2 >= t1)
+        for i in range(1000):
+            x = [0] * 1000
+        t3 = b.time()
+        self.assertTrue(t3 >= t2)
+        b.reset()
+        t4 = b.time()
+        self.assertTrue(t4 < t3)
+
+        self.assertEqual(b.format(1), '1 second')
+        self.assertEqual(b.format(61), '1 minute, 1 second')
+        self.assertEqual(b.format(60), '1 minute, 0 seconds')
+        self.assertEqual(b.format(180), '3 minutes, 0 seconds')
+        self.assertEqual(b.format(3600), '1 hour, 0 minutes, 0 seconds')
+        self.assertEqual(b.format(3661), '1 hour, 1 minute, 1 second')
+        self.assertEqual(
+            b.format(3600 * 24), '1 day, 0 hours, 0 minutes, 0 seconds')
+        self.assertEqual(
+            b.format(3600 * 24 * 7),
+            '1 week, 0 days, 0 hours, 0 minutes, 0 seconds')
+
+
+class CaptureTest(unittest.TestCase):
+    """Test the ``capture`` context manager."""
 
     def test_capture_nested(self):
         """Tests capturing in a nested pattern."""
-        r = myokit.capture(False)
-        q = myokit.capture(True)
+        r = myokit.tools.capture(False)
+        q = myokit.tools.capture(True)
         self.assertEqual(r.out(), '')
         self.assertEqual(r.err(), '')
-        with myokit.capture(False) as p:
+        with myokit.tools.capture(False) as p:
             self.assertEqual(p.out(), '')
             self.assertEqual(p.err(), '')
             print('1', end='')
@@ -44,7 +80,7 @@ def test_capture_nested(self):
                 with r:
                     print('3', end='')
                     print('c', end='', file=sys.stderr)
-                    with myokit.capture(True) as s:
+                    with myokit.tools.capture(True) as s:
                         print('4', end='')
                         print('d', end='', file=sys.stderr)
                     print('5', end='')
@@ -73,8 +109,8 @@ def test_capture_nested(self):
 
     def test_capture_repeated_use(self):
         """Tests capturing in a nested pattern with repeated enters/exits."""
-        x = myokit.capture()
-        y = myokit.capture(True)
+        x = myokit.tools.capture()
+        y = myokit.tools.capture(True)
         self.assertEqual(x.out(), '')
         with x:
             print('1', end='')
@@ -128,7 +164,7 @@ def test_capture_with_threads(self):
 
         # Function to call inside threads
         def f(i):
-            with myokit.capture() as c:
+            with myokit.tools.capture() as c:
                 print(str(i) + 'a ', end='')
                 time.sleep(times[i])
                 print(str(i) + 'b ', end='')
@@ -154,5 +190,64 @@ def f(i):
         self.assertEqual(captured[2], '2a 2b 2e')
 
 
+class ToolsTest(unittest.TestCase):
+    """Tests various tools in myokit.tools"""
+
+    def test_format_path(self):
+        # Test format_path().
+        fp = myokit.tools.format_path
+
+        # Normal use
+        self.assertEqual(
+            fp(os.path.join('a', 'b', 'c')), os.path.join('a', 'b', 'c'))
+
+        # No trailing slash
+        self.assertEqual(fp('a'), 'a')
+        self.assertEqual(fp('a/b/'), os.path.join('a', 'b'))
+
+        # Use with custom root
+        root = os.path.join(os.path.abspath('.'), 'a')
+        self.assertEqual(
+            fp(os.path.join(root, 'b', 'c'), root), os.path.join('b', 'c'))
+
+        # Empty path
+        self.assertEqual(fp(''), '.')
+        self.assertEqual(fp('.'), '.')
+
+        # Filesystem root
+        self.assertEqual(fp('/'), os.path.abspath('/'))
+        self.assertEqual(fp('/', root='/'), '.')
+
+        # Path outside of root
+        self.assertEqual(
+            fp(os.path.abspath('test'), os.path.abspath('test/tost')),
+            os.path.abspath('test'))
+
+    def test_levenshtein_distance(self):
+        # Test the levenshtein distance method.
+        self.assertEqual(myokit.tools.lvsd('kitten', 'sitting'), 3)
+        self.assertEqual(myokit.tools.lvsd('sitting', 'kitten'), 3)
+        self.assertEqual(myokit.tools.lvsd('saturday', 'sunday'), 3)
+        self.assertEqual(myokit.tools.lvsd('sunday', 'saturday'), 3)
+        self.assertEqual(myokit.tools.lvsd('michael', 'jennifer'), 7)
+        self.assertEqual(myokit.tools.lvsd('jennifer', 'michael'), 7)
+        self.assertEqual(myokit.tools.lvsd('jennifer', ''), 8)
+        self.assertEqual(myokit.tools.lvsd('', 'jennifer'), 8)
+        self.assertEqual(myokit.tools.lvsd('', ''), 0)
+
+    def test_natural_sort_key(self):
+        # Test natural sort key method.
+
+        a = ['a12', 'a3', 'a11', 'a2', 'a10', 'a1']
+        b = ['a1', 'a2', 'a3', 'a10', 'a11', 'a12']
+        self.assertNotEqual(a, b)
+        a.sort()
+        self.assertNotEqual(a, b)
+        a.sort(key=myokit.tools.natural_sort_key)
+        self.assertEqual(a, b)
+
+    # TODO: rmtree
+
+
 if __name__ == '__main__':
     unittest.main()