Repository: ModellingWebLab/cellmlmanip

Issue:
Rename/move some Model functions
1. Model. add_number()
This function merely creates a NumberDummy object. It is completely independent of the Model and adds nothing to anything. I think the name could be misleading.

I would prefer: create_number_symbol 

2. Model.check_left_right_units_equal(sympy.Equality)
This functions should be in the one of the unit classes (UnitCalculator); although that would involve passing the model/model units as an argument - so I could let it stay; but it does feel wrong. 

3. Model.find_symbols_and_derivatives(expression)
This creates a set of symbols from the list of expressions passed to it. It does not find symbols/derivatives that are in the model. Again possible misleading.  Also not sure this belongs on Model class (same as add_number) they are helper functions really.

Discussion:

Patch:
diff --git a/cellmlmanip/main.py b/cellmlmanip/main.py
index 42764f1..f7e907b 100644
--- a/cellmlmanip/main.py
+++ b/cellmlmanip/main.py
@@ -3,7 +3,7 @@ from cellmlmanip.parser import Parser
 
 
 def load_model(path, unit_store=None):
-    """Parses a CellML file and returns a :class:`cellmlmanip.Model`.
+    """Parses a CellML file and returns a :class:`cellmlmanip.model.Model`.
 
     :param unit_store: Optional :class:`cellmlmanip.units.UnitStore` instance; if given the model will share the
         underlying registry so that conversions between model units and those from the provided store work.
diff --git a/cellmlmanip/model.py b/cellmlmanip/model.py
index a8e4d91..b76af4e 100644
--- a/cellmlmanip/model.py
+++ b/cellmlmanip/model.py
@@ -1,4 +1,7 @@
-"""Classes to represent a flattened CellML model and metadata about its variables."""
+"""
+The main construct in cellmlmanip is a :class:`cellmlmanip.model.Model`.
+This represents a flattened CellML model and metadata about its variables.
+"""
 import logging
 from enum import Enum
 from io import StringIO
@@ -21,11 +24,33 @@ FLOAT_PRECISION = 17
 
 
 class DataDirectionFlow(Enum):
-    """ Direction of data flow for converting units"""
+    """Direction of data flow for converting units."""
     INPUT = 1
     OUTPUT = 2
 
 
+class VariableType(Enum):
+    """Classification of variables according to their role in the model's mathematics.
+
+    ``UNKNOWN``
+        not yet classified
+    ``STATE``
+        the dependent variable in an ODE
+    ``FREE``
+        the independent variable in an ODE
+    ``PARAMETER``
+        defined directly as a constant number (note that this does not include variables defined by an
+        equation that evaluates as constant)
+    ``COMPUTED``
+        defined by any other equation
+    """
+    UNKNOWN = 0
+    STATE = 1
+    FREE = 2
+    PARAMETER = 3
+    COMPUTED = 4
+
+
 class Model(object):
     """
     A componentless representation of a CellML model, containing a list of equations, units, and RDF metadata about
@@ -34,10 +59,11 @@ class Model(object):
     The main parts of a Model are 1. a list of sympy equation objects; 2. a collection of named units; and 3. an RDF
     graph that stores further meta data about the model.
 
-    Equations are stored as ``Sympy.Eq`` objects, but with the caveat that all variables and numbers must be specified
-    using the ``Sympy.Dummy`` objects returned by :meth:`add_variable()` and :meth:`add_number()`.
+    Equations are stored as :class:`sympy.Eq` objects, but with the caveat that all variables and numbers must be
+    specified using the :class:`sympy.Dummy` objects returned by :meth:`add_variable()` and :meth:`add_number()`.
 
-    Units are handled using the ``units`` property of a model, which is an instance of ``cellmlmanip.units.UnitStore``.
+    Units are handled using the ``units`` property of a model, which is an instance of
+    :class:`cellmlmanip.units.UnitStore`.
 
     RDF meta data can be attached to either the model itself or to model variables, via the CellML ``cmeta:id``
     attribute. Cmeta ids set on any other parts of CellML models are ignored.
@@ -110,17 +136,20 @@ class Model(object):
         Adds an equation to this model.
 
         The left-hand side (LHS) of the equation must be either a variable (as a :class:`VariableDummy`) or a derivative
-        (as a ``sympy.Derivative``).
+        (as a :class:`sympy.Derivative`).
 
         All numbers and variables used in the equation must have been obtained from this model, e.g. via
         :meth:`add_number()`, :meth:`add_variable()`, or :meth:`get_variable_by_ontology_term()`.
 
-        :param equation: A ``sympy.Eq`` object.
+        :param equation: A :class:`sympy.Eq` object.
         :param check_duplicates: whether to check that the equation's LHS is not already defined
         """
         assert isinstance(equation, sympy.Eq), 'The argument `equation` must be a sympy.Eq.'
-        self.equations.append(equation)
         lhs = equation.lhs
+        if lhs.is_Derivative:
+            if len(lhs.args) > 2 or lhs.args[1][1] > 1:
+                raise ValueError('Only first order derivatives wrt a single variable are supported')
+        self.equations.append(equation)
         if lhs.is_Derivative:
             state_var = lhs.free_symbols.pop()
             if check_duplicates:
@@ -138,7 +167,7 @@ class Model(object):
         """
         Assert that a variable doesn't have an existing definition.
 
-        :param var: the VariableDummy, either a state var or normal var
+        :param var: the :class:`VariableDummy`, either a state var or normal var
         :param equation: the new definition being added
         """
         if var in self._ode_definition_map:
@@ -153,7 +182,7 @@ class Model(object):
         Creates and returns a :class:`NumberDummy` to represent a number with units in sympy expressions.
 
         :param number: A number (anything convertible to float).
-        :param units: A string unit name or a ``Unit`` object.
+        :param units: A string unit name or a :class:`~cellmlmanip.units.UnitStore.Unit` object.
 
         :return: A :class:`NumberDummy` object.
         """
@@ -170,7 +199,7 @@ class Model(object):
         Adds a variable to the model and returns a :class:`VariableDummy` to represent it in sympy expressions.
 
         :param name: A string name.
-        :param units: A string unit name or a ``Unit`` object.
+        :param units: A string unit name or a :class:`~cellmlmanip.units.UnitStore.Unit` object.
         :param initial_value: An optional initial value.
         :param public_interface: An optional public interface specifier (only required when parsing CellML).
         :param private_interface: An optional private interface specifier (only required when parsing CellML).
@@ -211,36 +240,25 @@ class Model(object):
 
         return var
 
-    def transform_constants(self):
-        """
-        Called by CellML parser to standardise handling of 'constants'.
-
-        Once this has been called, the only variables with an initial_value attribute will be state variables,
-        and the initial value will do what it implies - hold the value the state variable should take at t=0.
-
-        Non state variables with an initial value are actually just constants. For consistent processing later on we add
-        equations defining them, and remove the initial_value attribute.
-        """
-        for var in self._name_to_variable.values():
-            if var in self._ode_definition_map:
-                assert var.initial_value is not None, 'State variable {} has no initial_value set'.format(var)
-            elif var.initial_value is not None:
-                value = self.add_number(var.initial_value, var.units)
-                self.add_equation(sympy.Eq(var, value))
-                var.initial_value = None
+    def connect_variables(self, source_name, target_name):
+        """Combine two variables that represent the same entity within the model.
 
-    def connect_variables(self, source_name: str, target_name: str):
-        """Tells this model that the variable indicated by ``target_name`` should get its value from the variable
-        indicated by ``source_name``.
+        This method is used to implement CellML connections between variables in different components. It tells the
+        model that the variable indicated by ``target_name`` should get its value from the variable indicated by
+        ``source_name``. In general therefore this method should only be called by the :class:`Parser`.
 
-        If the variables' units are not the same, an equation will be added to the target component.
+        If the units of both variables match, then any equations referencing the target variable are changed to
+        reference the source directly, and metadata annotations are moved over. If the units differ, an equation is
+        added defining target in terms of source with a conversion factor.
 
         This method will only work if the source variable has been assigned a value, either through an equation or by
         connecting it (directly or indirectly) to a variable with an equation. If this is not yet the case, ``False`` is
         returned. If successful the method returns ``True``.
 
-        :param source_name: The source variable name
-        :param target_name: The target variable name
+        :param str source_name: the source variable name
+        :param str target_name: the target variable name
+        :raises ValueError: if a logically impossible connection is attempted
+        :raises DimensionalityError: if the units are incompatible
         """
         logger.debug('connect_variables(%s ⟶ %s)', source_name, target_name)
 
@@ -263,11 +281,12 @@ class Model(object):
             raise ValueError('Multiple definitions for {} ({} and {})'.format(
                 target, source, self._var_definition_map[target]))
 
-        # If source/target variable is in the same unit
-        if source.units == target.units:
+        # Check whether we need a unit conversion
+        cf = self.units.get_conversion_factor(from_unit=source.units, to_unit=target.units)
+        if cf == 1:
             # Direct substitution is possible
             target.assigned_to = source.assigned_to
-            # everywhere the target variable is used, replace with source variable,
+            # Everywhere the target variable is used, replace with source variable,
             # updating the definition maps accordingly
             for index, equation in enumerate(self.equations):
                 self.equations[index] = new_eq = equation.xreplace({target: source.assigned_to})
@@ -289,13 +308,10 @@ class Model(object):
 
         # Otherwise, this connection requires a conversion
         else:
-            # Get the scaling factor required to convert source units to target units
-            factor = self.units.convert(1 * source.units, target.units).magnitude
-
             # Dummy to represent this factor in equations, having units for conversion
-            factor_dummy = self.add_number(factor, target.units / source.units)
+            factor_dummy = self.add_number(cf, target.units / source.units)
 
-            # Add an equations making the connection with the required conversion
+            # Add an equation making the connection with the required conversion
             self.add_equation(sympy.Eq(target, source.assigned_to * factor_dummy))
 
             logger.info('Connection req. unit conversion: %s', self.equations[-1])
@@ -310,22 +326,10 @@ class Model(object):
 
         return True
 
-    def add_rdf(self, rdf: str):
-        """ Takes an RDF string and stores it in the model's RDF graph. """
+    def add_rdf(self, rdf):
+        """Takes an RDF string and stores it in the model's RDF graph."""
         self.rdf.parse(StringIO(rdf))
 
-    def check_left_right_units_equal(self, equality):
-        """
-        Checks whether the LHS and RHS in a ``sympy.Eq`` have the same units.
-        :param equality: A ``sympy.Eq``.
-        """
-        lhs_units = self.units.evaluate_units(equality.lhs)
-        rhs_units = self.units.evaluate_units(equality.rhs)
-        assert self.units.is_equivalent(rhs_units, lhs_units), 'Units %s %s != %s %s' % (
-            self.units.format(lhs_units), self.units.format(lhs_units, True),
-            self.units.format(rhs_units), self.units.format(rhs_units, True),
-        )
-
     def get_equations_for(self, variables, recurse=True, strip_units=True):
         """Get all equations for a given collection of variables.
 
@@ -333,8 +337,8 @@ class Model(object):
 
         :param variables: The variables to get the equations for (as :class:`VariableDummy` objects).
         :param recurse: Indicates whether to recurse the equation graph, or to return only the top level equations.
-        :param strip_units: If ``True``, all ``sympy.Dummy`` objects representing number with units will be replaced
-            with ordinary sympy number objects.
+        :param strip_units: If ``True``, all :class:`NumberDummy` objects representing number with units will be
+            replaced with ordinary sympy number objects.
         """
         # Get graph
         if strip_units:
@@ -372,19 +376,23 @@ class Model(object):
         return eqs
 
     def get_derivatives(self):
-        """Returns a list of ``sympy.Derivative`` objects found as LHS in the given model graph.
-        The list is ordered by appearance in the cellml document.
+        """Returns a list of :class:`sympy.Derivative` objects found as LHS in the given model graph.
+
+        The list is ordered by appearance in the CellML document.
         """
         derivatives = [v for v in self.graph if isinstance(v, sympy.Derivative)]
         return sorted(derivatives, key=lambda deriv: deriv.args[0].order_added)
 
     def get_derived_quantities(self):
         """Returns a list of derived quantities found in the given model graph.
-        A derived quantity is any variable that is not a state variable or parameter/constant.
+
+        A derived quantity is any variable that is not a state variable, free variable, or parameter/constant.
         """
-        derived_quantities = [v for v, node in self.graph.nodes.items()
-                              if not isinstance(v, sympy.Derivative)
-                              and node.get('variable_type', '') not in ('state', 'free', 'parameter')]
+        derived_quantities = [
+            v for v, node in self.graph.nodes.items()
+            if not isinstance(v, sympy.Derivative)
+            and node.get('variable_type', VariableType.UNKNOWN) not in (
+                VariableType.FREE, VariableType.STATE, VariableType.PARAMETER)]
         return sorted(derived_quantities, key=lambda var: var.order_added)
 
     def get_display_name(self, var, ontology=None):
@@ -394,6 +402,7 @@ class Model(object):
         specified), then cmeta:id if present, or the variable's name attribute if not.
 
         Dollar symbols in the name are replaced by a double underscore.
+
         :param var: the variable for which to get the display name.
         :param ontology: the base URL of an ontology if only annotations within that ontology should be considered
 
@@ -412,23 +421,25 @@ class Model(object):
 
     def get_free_variable(self):
         """Returns the free variable in this model (if any)."""
-        for v, node in self.graph.nodes.items():
-            if node.get('variable_type', '') == 'free':
-                return v
+        for ode in self._ode_definition_map.values():
+            free_variable = ode.lhs.variables[0]
+            return free_variable
 
         raise ValueError('No free variable set in model.')  # pragma: no cover
 
     def get_rdf_annotations(self, subject=None, predicate=None, object_=None):
-        """Searches the RDF graph and returns 'triples matching the given parameters'
+        """Searches the RDF graph and returns triples matching the given parameters.
 
         :param subject: the subject of the triples returned
         :param predicate: the predicate of the triples returned
         :param object_: the object of the triples returned
 
-        ``subject`` ``predicate`` and ``object_`` are optional, if None then any triple matches
-        if all are none, all triples are returned
-        ``subject`` ``predicate`` and ``object_`` can be anything valid as input to create_rdf_node
-        typically an (NS, local) pair, a string or None"""
+        Each of ``subject``, ``predicate`` and ``object_`` are optional; if ``None`` then any triple matches.
+        If all are ``None``, then all triples are returned.
+
+        The arguments can be anything valid as input to :meth:`cellmlmanip.rdf.create_rdf_node`,
+        typically a (namespace URI, local name) pair, a string or ``None``.
+        """
         subject = create_rdf_node(subject)
         predicate = create_rdf_node(predicate)
         object_ = create_rdf_node(object_)
@@ -440,7 +451,8 @@ class Model(object):
         :param subject: the object of the triple returned
         :param predicate: the object of the triple returned
 
-        Note: expects exactly one triple to match and the result to be a literal. It's string value is  returned."""
+        Note: expects exactly one triple to match and the result to be a literal. Its string value is returned.
+        """
         triples = list(self.get_rdf_annotations(subject, predicate))
         assert len(triples) == 1
         assert isinstance(triples[0][2], rdflib.Literal)
@@ -451,7 +463,7 @@ class Model(object):
         """
         Searches the model and returns the variable with the given cmeta id.
 
-        To get variables from e.g. an oxmeta ontology term, use :meth:`get_variable_by_ontology_term()`.
+        To get variables from e.g. an oxmeta ontology term, use :meth:`get_variable_by_ontology_term`.
 
         :param cmeta_id: Either a string id or :class:`rdflib.URIRef` instance.
         :returns: A :class:`VariableDummy` object
@@ -475,21 +487,21 @@ class Model(object):
             raise KeyError('No variable with cmeta id "%s" found.' % str(cmeta_id))
 
     def get_variable_by_name(self, name):
-        """ Returns the variable with the given ``name``. """
+        """Returns the variable with the given ``name``."""
         return self._name_to_variable[name]
 
     def get_variable_by_ontology_term(self, term):
         """Searches the RDF graph for a variable annotated with the given ``term`` and returns it.
 
         Specifically, this method searches for a unique variable annotated with
-        predicate ``http://biomodels.net/biology-qualifiers/is`` and the object
+        predicate http://biomodels.net/biology-qualifiers/is and the object
         specified by ``term``.
 
         Will raise a ``KeyError`` if no variable with the given annotation is
         found, and a ``ValueError`` if more than one variable with the given
         annotation is found.
 
-        :param term: anything suitable as an input to :meth:`create_rdf_node`; typically either an RDF
+        :param term: anything suitable as an input to :meth:`cellmlmanip.rdf.create_rdf_node`; typically either an RDF
             node already, or a tuple ``(namespace_uri, local_name)``.
         """
         variables = self.get_variables_by_rdf(('http://biomodels.net/biology-qualifiers/', 'is'), term)
@@ -501,10 +513,12 @@ class Model(object):
             raise ValueError('Multiple variables annotated with {}'.format(term))
 
     def get_variables_by_rdf(self, predicate, object_=None):
-        """Searches the RDF graph for variables annotated with the given predicate and object (e.g. ``is oxmeta:time``)
-        and returns the associated variables sorted in document order.
+        """Find variables annotated with the given predicate and object (e.g. ``is oxmeta:time``) in our RDF graph.
+
+        Both ``predicate`` and ``object_`` (if given) must be suitable as an input to
+        :meth:`cellmlmanip.rdf.create_rdf_node`; typically either ``(namespace, local_name)`` tuples or string literals.
 
-        Both ``predicate`` and ``object_`` (if given) must be ``(namespace, local_name)`` tuples or string literals.
+        :return: the associated variables sorted in document order
         """
         predicate = create_rdf_node(predicate)
         object_ = create_rdf_node(object_)
@@ -515,13 +529,13 @@ class Model(object):
 
     def get_ontology_terms_by_variable(self, variable, namespace_uri=None):
         """
-        Returns all ontology terms linked to the variable ``variable`` via the
-        ``http://biomodels.net/biology-qualifiers/is`` predicate.
+        Returns all ontology terms linked to the given ``variable`` via the
+        http://biomodels.net/biology-qualifiers/is predicate.
 
         :param variable: The variable to search for (as a :class:`VariableDummy` object).
         :param namespace_uri: An optional namespace URI. If given, only terms within the given namespace will be
             returned.
-        :returns: A list of term names.
+        :returns: A list of term local names.
         """
         ontology_terms = []
         if variable.rdf_identity:
@@ -536,7 +550,7 @@ class Model(object):
     def get_definition(self, variable):
         """Get the equation (if any) defining the given variable.
 
-        :param variable: The variable to look up (as a class:`VariableDummy`. If this appears as the LHS of a straight
+        :param variable: The variable to look up (as a :class:`VariableDummy`). If this appears as the LHS of a straight
             assignment, or the state variable in an ODE, the corresponding equation will be returned.
         :returns: A Sympy equation, or ``None`` if the variable is not defined by an equation.
         """
@@ -546,13 +560,12 @@ class Model(object):
         return defn
 
     def get_value(self, variable):
-        """ Returns the evaluated value of the given variable's RHS. """
+        """Returns the evaluated value of the given variable's RHS."""
         return float(self.graph.nodes[variable]['equation'].rhs)
 
     @property
     def graph(self):
-        """ A ``networkx.DiGraph`` containing the model equations. """
-        # TODO: Set the parameters of the model (parameters rather than use initial values)
+        """A :class:`networkx.DiGraph` containing the model equations."""
 
         # Return cached graph
         if self._graph is not None:
@@ -575,15 +588,15 @@ class Model(object):
             if lhs.is_Derivative:
                 # Get the state symbol and update the variable information
                 state_symbol = lhs.free_symbols.pop()
-                state_symbol.type = 'state'
+                state_symbol.type = VariableType.STATE
 
                 # Get the free symbol and update the variable information
                 free_symbol = lhs.variables[0]
-                free_symbol.type = 'free'
+                free_symbol.type = VariableType.FREE
             elif isinstance(equation.rhs, NumberDummy):
-                lhs.type = 'parameter'
+                lhs.type = VariableType.PARAMETER
             else:
-                lhs.type = 'computed'
+                lhs.type = VariableType.COMPUTED
 
         # Sanity check: none of the lhs have the same hash
         assert len(graph.nodes) == equation_count
@@ -601,7 +614,7 @@ class Model(object):
                 if rhs in graph.nodes:
                     # If the symbol maps to a node in the graph just add the dependency edge
                     graph.add_edge(rhs, lhs)
-                elif rhs.type in ['state', 'free']:
+                elif rhs.type in [VariableType.STATE, VariableType.FREE]:
                     # If the variable is a state or free variable of a derivative
                     graph.add_node(rhs, equation=None, variable_type=rhs.type)
                     graph.add_edge(rhs, lhs)
@@ -630,8 +643,8 @@ class Model(object):
     @property
     def graph_with_sympy_numbers(self):
         """
-        A ``networkx.DiGraph`` containing the model equations, but with numbers represented as sympy ``Number`` objects
-        instead of dummies.
+        A :class:`networkx.DiGraph` containing the model equations,
+        but with numbers represented as :class:`sympy.Number` objects instead of :class:`NumberDummy`.
         """
         if self._graph_with_sympy_numbers is not None:
             return self._graph_with_sympy_numbers
@@ -684,18 +697,19 @@ class Model(object):
         return cmeta_id in self._cmeta_id_to_variable
 
     def has_ontology_annotation(self, variable, namespace_uri=None):
-        """Checks that there is at least one result for
-        :meth:`Model.get_ontology_terms_by_variable(variable, namespace_uri)`.
+        """
+        Checks that there is at least one result for
+        :meth:`Model.get_ontology_terms_by_variable` with the given arguments.
         """
         return len(self.get_ontology_terms_by_variable(variable, namespace_uri)) != 0
 
     def _invalidate_cache(self):
-        """ Removes cached graphs: should be called after manipulating variables or equations. """
+        """Removes cached graphs: should be called after manipulating variables or equations."""
         self._graph = None
         self._graph_with_sympy_numbers = None
 
     def is_state(self, variable):
-        """ Checks if the given ``variable`` is a state variable (i.e. if it's defined by an ODE). """
+        """Checks if the given ``variable`` is a state variable (i.e. if it's defined by an ODE)."""
         return variable in self._ode_definition_map
 
     def is_constant(self, variable):
@@ -708,13 +722,14 @@ class Model(object):
         return defn is not None and len(defn.rhs.atoms(VariableDummy)) == 0
 
     def find_variables_and_derivatives(self, expressions):
-        """ Returns a set containing all variables and derivatives referenced in a list of expressions.
+        """Returns a set containing all variables and derivatives referenced in a list of expressions.
 
         Note that we can't just use ``.atoms(VariableDummy, sympy.Derivative)`` for this, because it
         will return the state and free variables from inside derivatives, which is not what we want.
 
         :param expressions: an iterable of expressions to get variables for.
-        :return: a set of variables and derivatives, as :class:`VariableDummy` and ``sympy.Derivative`` respectively.
+        :return: a set of variables and derivatives, as :class:`VariableDummy` and :class:`sympy.Derivative`
+            objects respectively.
         """
         variables = set()
         for expr in expressions:
@@ -789,7 +804,7 @@ class Model(object):
         self._cmeta_id_to_variable[target._cmeta_id] = target
 
     def variables(self):
-        """ Returns an iterator over this model's variables. """
+        """Returns an iterator over this model's variables."""
         return self._name_to_variable.values()
 
     def convert_variable(self, original_variable, units, direction):
@@ -859,11 +874,12 @@ class Model(object):
 
         :param original_variable: the :class:`VariableDummy` object representing the variable in the model to be
                                   converted
-        :param units: a Pint unit object representing the units to convert variable to (note if variable is already
-                      in these units, model remains unchanged and the original variable is returned
-        :param direction: either DataDirectionFlow.INPUT; the variable to be changed is an input and all affected
-                          equations will be adjusted
-                          or DataDirectionFlow.OUTPUT; the variable to be changed is an output, equations
+        :param units: a :class:`~cellmlmanip.units.UnitStore.Unit` object representing the units to convert variable to
+                      (note if variable is already in these units, model remains unchanged and the original variable is
+                      returned)
+        :param direction: either DataDirectionFlow.INPUT: the variable to be changed is an input and all affected
+                          equations will be adjusted;
+                          or DataDirectionFlow.OUTPUT: the variable to be changed is an output, equations
                           are unaffected apart from converting the actual output
         :return: new variable with desired units, or original unchanged if conversion was not necessary
         :raises DimensionalityError: if the unit conversion is impossible
@@ -1057,7 +1073,7 @@ class Model(object):
 
     def get_unique_name(self, name):
         """
-        Creates and returns a unique name, not used in the model.
+        Creates and returns a unique variable name, not used in the model.
 
         :param str name: Suggested unique name.
         :return str: Guaranteed unique name.
@@ -1152,7 +1168,6 @@ class VariableDummy(sympy.Dummy):
         self._set_cmeta_id(cmeta_id)
 
         # This variable's type
-        # TODO: Define allowed types via enum
         self.type = None
 
     def __str__(self):
diff --git a/cellmlmanip/parser.py b/cellmlmanip/parser.py
index 3699aeb..f6df255 100644
--- a/cellmlmanip/parser.py
+++ b/cellmlmanip/parser.py
@@ -1,12 +1,14 @@
-"""This module contains the CellML parser and related classes. It reads a CellML model and stores
-model information in the cellmlmanip.Model class. MathML equations are translated to Sympy. RDF is
-handled by RDFLib.
+"""
+This :mod:`cellmlmanip.parser` module contains the CellML parser and related classes. It reads a CellML model and
+stores model information in the :class:`cellmlmanip.model.Model` class. MathML equations are translated to Sympy. RDF
+is handled by RDFLib.
 """
 import itertools
 import os
 from collections import OrderedDict, deque
 from enum import Enum
 
+import sympy
 from lxml import etree
 
 from cellmlmanip.model import SYMPY_SYMBOL_DELIMITER, Model
@@ -131,7 +133,7 @@ class Parser(object):
         self._add_connection(model_xml)
 
         # Canonicalise representation
-        self.model.transform_constants()
+        self.transform_constants()
 
         return self.model
 
@@ -485,6 +487,24 @@ class Parser(object):
         raise ValueError('Cannot determine the source & target for connection (%s, %s) - (%s, %s)' %
                          (comp_1, var_1, comp_2, var_2))
 
+    def transform_constants(self):
+        """
+        Standardise handling of 'constants'.
+
+        Once this has been called, the only variables with an initial_value attribute will be state variables,
+        and the initial value will do what it implies - hold the value the state variable should take at t=0.
+
+        Non state variables with an initial value are treated as constants. For consistent processing later on we add
+        equations defining them, and remove the initial_value attribute.
+        """
+        for var in self.model.variables():
+            if var in self.model.get_state_variables():
+                assert var.initial_value is not None, 'State variable {} has no initial_value set'.format(var)
+            elif var.initial_value is not None:
+                value = self.model.add_number(var.initial_value, var.units)
+                self.model.add_equation(sympy.Eq(var, value))
+                var.initial_value = None
+
     def _validate(self, parser, tree):
         """
         Validates the given lxml ``tree`` against the CellML 1.0 RELAX NG schema.
diff --git a/cellmlmanip/rdf.py b/cellmlmanip/rdf.py
index 7fbfab5..979b775 100644
--- a/cellmlmanip/rdf.py
+++ b/cellmlmanip/rdf.py
@@ -8,10 +8,11 @@ def create_rdf_node(node_content):
     :param node_content: the content for the node.
 
     The ``node_content``, if given, must either be a :class:`rdflib.term.Node` instance,
-    a tuple ``(namespace, local_name)``, or a string, in which case it is interpreted
+    a tuple ``(namespace_uri, local_name)``, or a string, in which case it is interpreted
     as either a URI ref if it starts with # otherwise a literal RDF node.
 
-    ``node_content`` may also be ``None`` to return `None`, allowing easy handling of wildcard options to queries."""
+    ``node_content`` may also be ``None`` to return ``None``, allowing easy handling of wildcard options to queries.
+    """
     if node_content is None or isinstance(node_content, rdflib.term.Node):
         return node_content
     elif isinstance(node_content, tuple):
diff --git a/cellmlmanip/transpiler.py b/cellmlmanip/transpiler.py
index b2b1f79..64165a6 100644
--- a/cellmlmanip/transpiler.py
+++ b/cellmlmanip/transpiler.py
@@ -1,4 +1,5 @@
-"""Parses Content MathML and returns equivalent SymPy expressions
+"""
+The :mod:`cellmlmanip.transpiler` module parses Content MathML and returns equivalent SymPy expressions.
 
 Content Markup specification: https://www.w3.org/TR/MathML2/chapter4.html
 """
diff --git a/cellmlmanip/units.py b/cellmlmanip/units.py
index 6c37439..4bdd090 100644
--- a/cellmlmanip/units.py
+++ b/cellmlmanip/units.py
@@ -1,4 +1,7 @@
-"""Unit handling for CellML models, using the Pint unit library."""
+"""
+The :mod:`cellmlmanip.units` module provides unit handling for CellML models, using the
+`Pint unit library <https://pint.readthedocs.io/>`_.
+"""
 import logging
 import math
 import numbers
@@ -180,13 +183,13 @@ class UnitStore(object):
     Creates and stores units, and has functions for unit conversion.
 
     Within a UnitStore ``store``, units are represented as ``store.Unit`` objects, while numbers with units are
-    represented as ``store.Quantity`` objects. Both classes are inherited from the ``pint`` package.
+    represented as ``store.Quantity`` objects. Both classes are inherited from the :mod:`pint` package.
 
     Each UnitStore has its own unique namespace for storing units.
 
-    By default, each ``UnitStore`` maintains an internal ``pint.UnitRegistry`` and units from different stores cannot be
-    compared or interact with each other. To allow comparison and conversion between unit stores, an existing
-    ``UnitStore`` can be passed in at construction time, so that the underlying registry will be shared. (This will
+    By default, each UnitStore maintains an internal :class:`pint.UnitRegistry` and units from different stores cannot
+    be compared or interact with each other. To allow comparison and conversion between unit stores, an existing
+    UnitStore can be passed in at construction time, so that the underlying registry will be shared. (This will
     allow units from different stores to interact, but the unit stores will maintain independent namespaces).
 
     For example::
@@ -207,6 +210,14 @@ class UnitStore(object):
     stores, allowing e.g. units from ``c`` to be converted to units from ``a``.
 
     :param store: An existing :class:`UnitStore` to share a unit registry with.
+
+    .. class:: Unit
+
+        A :class:`pint.Unit` class tied to this unit store.
+
+    .. class:: Quantity
+
+        A :class:`pint.Quantity` class tied to this unit store.
     """
     _next_id = 0
 
@@ -302,9 +313,9 @@ class UnitStore(object):
     def convert(self, quantity, unit):
         """Converts the given ``quantity`` to be in the given ``unit``.
 
-        :param quantity: a ``Quantity`` to convert.
-        :param unit: a ``Unit`` to convert it to.
-        :returns: the converted ``Quantity``.
+        :param quantity: a :class:`.Quantity` to convert.
+        :param unit: a :class:`.Unit` to convert it to.
+        :returns: the converted :class:`.Quantity`.
         """
         assert isinstance(quantity, self._registry.Quantity)
         assert isinstance(unit, self._registry.Unit)
@@ -314,7 +325,7 @@ class UnitStore(object):
         """
         Returns a string representation of a unit.
 
-        :param unit: The ``Unit`` object to format.
+        :param unit: The :class:`.Unit` object to format.
         :param base_units: If this optional argument is set to ``True`` the string will show the base units expansion of
             the given unit.
         """
@@ -327,10 +338,10 @@ class UnitStore(object):
             return _STORE_PREFIX.sub('', str(unit))
 
     def get_unit(self, name):
-        """Retrieves the ``Unit`` object with the given name.
+        """Retrieves the :class:`.Unit` object with the given name.
 
         :param unit_name: string name of the unit
-        :returns: A ``Unit`` object.
+        :returns: A :class:`Unit` object.
         :raises KeyError: If the unit is not defined in this unit store.
         """
         if name not in self._known_units:
@@ -345,8 +356,8 @@ class UnitStore(object):
     def is_equivalent(self, unit1, unit2):
         """Tests whether two units are equivalent.
 
-        :param unit1: The first ``Unit`` object to compare.
-        :param unit2: The second ``Unit`` object to compare.
+        :param unit1: The first :class:`Unit` object to compare.
+        :param unit2: The second :class:`Unit` object to compare.
         :returns: ``True`` if both units are equivalent, ``False`` otherwise.
         """
         assert isinstance(unit1, self.Unit)
@@ -359,10 +370,10 @@ class UnitStore(object):
 
     def evaluate_units(self, expr):
         """
-        Evaluates and returns the ``Unit`` a Sympy expression is in.
+        Evaluates and returns the :class:`Unit` a Sympy expression is in.
 
         :param expr: the Sympy expression whose units to evaluate.
-        :returns: The calculated ``Unit`` object.
+        :returns: The calculated :class:`Unit` object.
         :raises UnitError: if there are unit errors when evaluating the expression's units.
         """
         found = self._calculator.traverse(expr).units
@@ -371,12 +382,12 @@ class UnitStore(object):
 
     def evaluate_units_and_fix(self, expr):
         """
-        Evaluates and returns the ``Unit`` a Sympy expression is in; but will also attempt to fix inconsistencies in
-        ``expr`` and return an updated expression if needed.
+        Evaluates and returns the :class:`Unit` a Sympy expression is in; but will also attempt to fix inconsistencies
+        in ``expr`` and return an updated expression if needed.
 
         :param expr: the Sympy expression whose units to evaluate.
-        :returns: A tuple ``(units, new_expr)`` where ``units`` is the calculated ``Unit`` object, and ``new_expr`` is
-            either ``expr`` or a copy made internally consistent.
+        :returns: A tuple ``(units, new_expr)`` where ``units`` is the calculated :class:`Unit` object, and ``new_expr``
+            is either ``expr`` or a copy made internally consistent.
         :raises UnitError: if there are unfixable unit errors when evaluating the expression's units.
         """
         try:
@@ -389,8 +400,8 @@ class UnitStore(object):
     def get_conversion_factor(self, from_unit, to_unit):
         """Returns the magnitude multiplier required to convert a unit to the specified unit.
 
-        :param from_unit: the ``Unit`` to be converted
-        :param to_unit: ``Unit`` object into which the units should be converted
+        :param from_unit: the :class:`Unit` to be converted
+        :param to_unit: :class:`Unit` object into which the units should be converted
         :returns: the magnitude of the resulting conversion factor
         """
         assert isinstance(from_unit, self._registry.Unit), 'from_unit must be a unit, not ' + str(from_unit)
@@ -424,20 +435,21 @@ class UnitStore(object):
         desired by the LHS, if  the ``Eq`` expression is passed in as ``expr`` and ``to_units`` is given as ``None``.
 
         The conversion strategy for each (sub-)expression depends on the operator:
-        - for relational operators, all operands are converted to the units of the first operand
-        - for Mul the operands can be in any units, and we convert the result if needed
-        - for Add all operands are converted to the desired units (or the units of the first operand if no desired units
+
+        * for relational operators, all operands are converted to the units of the first operand
+        * for Mul the operands can be in any units, and we convert the result if needed
+        * for Add all operands are converted to the desired units (or the units of the first operand if no desired units
           are given)
-        - for Pow the exponent must be dimensionless while the operand can be in any units, and we convert the result if
+        * for Pow the exponent must be dimensionless while the operand can be in any units, and we convert the result if
           needed
-        - for trig functions, exp, log, etc. the operands have to have dimensionless units
-        - for piecewise, the conditions must be dimensionless, and the pieces are set to the desired units (or the units
+        * for trig functions, exp, log, etc. the operands have to have dimensionless units
+        * for piecewise, the conditions must be dimensionless, and the pieces are set to the desired units (or the units
           of the first piece)
-        - for derivatives, numbers, variables, etc. we just convert to the desired units
+        * for derivatives, numbers, variables, etc. we just convert to the desired units
 
         :param expr: the Sympy expression to convert
-        :param to_units: the desired units of the expression, or ``None`` if we don't care or for converting an
-            assignment expression.
+        :param to_units: the desired units of the expression as a :class:`Unit` object, or ``None`` if we don't care or
+            for converting an assignment expression.
         :returns: a Sympy expression in the desired units; the input ``expr`` if no conversion was needed.
         :raises UnitError: if conversion is not possible, using a suitable subclass depending on the exact reason
         """
diff --git a/docs/source/code_generation.rst b/docs/source/code_generation.rst
index 34dfbfc..5cc76cc 100644
--- a/docs/source/code_generation.rst
+++ b/docs/source/code_generation.rst
@@ -6,4 +6,4 @@ Sympy expressions can be converted to code strings using a :class:`cellmlmanip.p
 
 .. currentmodule:: cellmlmanip.printer
 
-.. autoclass :: Printer
+.. automodule:: cellmlmanip.printer
diff --git a/docs/source/conf.py b/docs/source/conf.py
index 6607a9e..36daa0f 100644
--- a/docs/source/conf.py
+++ b/docs/source/conf.py
@@ -20,14 +20,14 @@
 # import sys
 # sys.path.insert(0, os.path.abspath('.'))
 # import cellmlmanip
-import sphinx
+# import sphinx
 
 
 # -- General configuration ------------------------------------------------
 
 # If your documentation needs a minimal Sphinx version, state it here.
 #
-# needs_sphinx = '1.0'
+needs_sphinx = '2.0'
 
 # Add any Sphinx extension module names here, as strings. They can be
 # extensions coming with Sphinx (named 'sphinx.ext.*') or your custom
@@ -35,23 +35,29 @@ import sphinx
 extensions = [
     'sphinx.ext.autodoc',
     'sphinx.ext.doctest',
+    'sphinx.ext.intersphinx',
     # 'sphinx.ext.mathjax',
     # 'sphinx.ext.napoleon',
     'sphinx.ext.viewcode',
 ]
 
 # Autodoc defaults
-if int(sphinx.__version__.split('.')[1]) < 8:
-    autodoc_default_flags = [
-        'members',
-        'inherited-members',
-        # 'show-inheritance',
-    ]
-else:
-    autodoc_default_options = {
-        'members': None,
-        'inherited-members': None,
-    }
+autodoc_default_options = {
+    'members': True,
+    'member-order': 'bysource',
+    'inherited-members': False,
+    'show-inheritance': True,
+}
+
+# Where to find documentation for projects we use
+intersphinx_mapping = {
+    'networkx': ('https://networkx.github.io/documentation/stable/', None),
+    'rdflib': ('https://rdflib.readthedocs.io/en/stable/', None),
+    'sympy': ('https://docs.sympy.org/latest/', None),
+    'pint': ('https://pint.readthedocs.io/en/stable/', None),
+    'python': ('https://docs.python.org/3', None),
+}
+intersphinx_timeout = 10
 
 # Add any paths that contain templates here, relative to this directory.
 templates_path = ['_templates']
@@ -67,7 +73,7 @@ master_doc = 'index'
 
 # General information about the project.
 project = u'Cellmlmanip'
-copyright = u'2019, Cellmlmanip authors'
+copyright = u'2019-, Cellmlmanip authors'
 author = u'Cellmlmanip authors'
 
 # The version info for the project you're documenting, acts as replacement for
@@ -75,7 +81,7 @@ author = u'Cellmlmanip authors'
 # built documents.
 #
 # The short X.Y version.
-version = '0.0.1'
+version = '0.1.0'
 # The full version, including alpha/beta/rc tags.
 release = version
 
diff --git a/docs/source/model_api.rst b/docs/source/model_api.rst
index bc72f93..5127016 100644
--- a/docs/source/model_api.rst
+++ b/docs/source/model_api.rst
@@ -2,6 +2,7 @@
 The Model API
 **************
 
-.. currentmodule:: cellmlmanip.model
 
-.. autoclass :: Model
+.. automodule:: cellmlmanip.model
+
+.. autofunction:: cellmlmanip.rdf.create_rdf_node
diff --git a/docs/source/reading_cellml.rst b/docs/source/reading_cellml.rst
index 81668da..4f4819c 100644
--- a/docs/source/reading_cellml.rst
+++ b/docs/source/reading_cellml.rst
@@ -7,3 +7,10 @@ Reading CellML
 To start using cellmlmanip, load a model using
 
 .. autofunction :: load_model
+
+This uses the :class:`cellmlmanip.parser.Parser` class.
+Most users won't need to worry about this, but its documentation is included below for completeness.
+
+.. automodule:: cellmlmanip.parser
+
+.. automodule:: cellmlmanip.transpiler
diff --git a/docs/source/units.rst b/docs/source/units.rst
index 3583623..78d419c 100644
--- a/docs/source/units.rst
+++ b/docs/source/units.rst
@@ -3,23 +3,8 @@ Units
 *****
 
 Units are stored in a :class:`cellmlmanip.units.UnitStore`.
-When a model is parsed, a `UnitStore` is created automatically.
+When a model is parsed, a ``UnitStore`` is created automatically.
 
 .. currentmodule:: cellmlmanip.units
 
-.. autoclass :: UnitStore
-
-.. autoclass :: UnitError
-
-.. autoclass :: UnitsCannotBeCalculatedError
-
-.. autoclass :: UnexpectedMathUnitsError
-
-.. autoclass :: InputArgumentsInvalidUnitsError
-
-.. autoclass :: InputArgumentsMustBeDimensionlessError
-
-.. autoclass :: InputArgumentMustBeNumberError
-
-.. autoclass :: BooleanUnitsError
-
+.. automodule:: cellmlmanip.units
diff --git a/setup.py b/setup.py
index 918b045..6923e3c 100644
--- a/setup.py
+++ b/setup.py
@@ -9,7 +9,7 @@ with open('LICENSE') as f:
 
 setup(
     name='cellmlmanip',
-    version='0.0.1',
+    version='0.1.0',
     description='CellML loading and model equation manipulation',
     long_description=readme,
     author='',
@@ -28,7 +28,7 @@ setup(
     ],
     extras_require={
         'docs': [
-            'sphinx>=1.5, !=1.7.3',
+            'sphinx>=2.0',
         ],
         'test': [
             'codecov',


Test patch:
diff --git a/tests/shared.py b/tests/shared.py
index f100e2d..1b717d1 100644
--- a/tests/shared.py
+++ b/tests/shared.py
@@ -14,3 +14,17 @@ def load_model(name, unit_store=None):
     return cellmlmanip.load_model(
         os.path.join(os.path.dirname(__file__), 'cellml_files', name),
         unit_store=unit_store)
+
+
+def check_left_right_units_equal(unit_store, equality):
+    """
+    Checks whether the LHS and RHS in a ``sympy.Eq`` have the same units.
+    :param unit_store: a :class:`UnitStore`.
+    :param equality: A ``sympy.Eq``.
+    """
+    lhs_units = unit_store.evaluate_units(equality.lhs)
+    rhs_units = unit_store.evaluate_units(equality.rhs)
+    assert unit_store.is_equivalent(rhs_units, lhs_units), 'Units %s %s != %s %s' % (
+        unit_store.format(lhs_units), unit_store.format(lhs_units, True),
+        unit_store.format(rhs_units), unit_store.format(rhs_units, True),
+    )
diff --git a/tests/test_model_functions.py b/tests/test_model_functions.py
index 799d12a..388b6a8 100644
--- a/tests/test_model_functions.py
+++ b/tests/test_model_functions.py
@@ -434,6 +434,12 @@ class TestModelFunctions():
         model.add_equation(sp.Eq(sv1, 2.0), check_duplicates=False)
         model.add_equation(sp.Eq(sv1_def.lhs, 1.0), check_duplicates=False)
 
+        # Only first-order derivatives allowed
+        with pytest.raises(ValueError, match='Only first order derivatives'):
+            model.add_equation(sp.Eq(sp.Derivative(symbol, symbol1, symbol2), 1.0))
+        with pytest.raises(ValueError, match='Only first order derivatives'):
+            model.add_equation(sp.Eq(sp.Derivative(symbol, symbol1, symbol1), 1.0))
+
     def test_remove_equation(self, local_hh_model):
         """ Tests the Model.remove_equation method. """
 
@@ -700,20 +706,3 @@ class TestModelFunctions():
         c = aslanidi_model.get_variable_by_ontology_term((shared.OXMETA, 'membrane_capacitance'))
         assert not aslanidi_model.is_state(c)
         assert aslanidi_model.is_constant(c)
-
-    def test_transform_constants(self):
-        """ Tests Model.transform_constants(). """
-
-        # Parse model, which should convert an initial value to an equation
-        model = shared.load_model('initial_value_constant.cellml')
-        v = model.get_variable_by_name('A$a')
-
-        # Initial value should be None, if transform_constants worked
-        assert v.initial_value is None
-
-        # And its RHS value should be 1
-        assert model.get_value(v) == 1
-
-        # And the equations should have matching units
-        for eq in model.equations:
-            model.check_left_right_units_equal(eq)
diff --git a/tests/test_parser.py b/tests/test_parser.py
index ee1d22d..e536288 100644
--- a/tests/test_parser.py
+++ b/tests/test_parser.py
@@ -4,7 +4,9 @@ import os
 import pytest
 import sympy
 
-from cellmlmanip import load_model, parser
+from cellmlmanip import parser
+
+from .shared import check_left_right_units_equal, load_model
 
 
 logger = logging.getLogger(__name__)
@@ -250,7 +252,7 @@ class TestParser(object):
         p = parser.Parser(cellml)
         model = p.parse()
         for e in model.equations:
-            model.check_left_right_units_equal(e)
+            check_left_right_units_equal(model.units, e)
 
     def test_units_with_multiplier(self):
         """Tests parsing a unit with a multiplier."""
@@ -274,10 +276,7 @@ class TestParser(object):
 
     def test_multiple_math_elements(self):
         """ Tests parser correctly handles a component with more than one math element. """
-        example_cellml = os.path.join(
-            os.path.dirname(__file__), "cellml_files", "3.4.2.1.component_with_maths.cellml"
-        )
-        model = load_model(example_cellml)
+        model = load_model('3.4.2.1.component_with_maths.cellml')
         assert len(list(model.equations)) == 2
 
     def test_bad_unit_definitions(self):
@@ -310,6 +309,23 @@ class TestParser(object):
         with pytest.raises(ValueError, match='Reactions are not supported'):
             p.parse()
 
+    def test_transform_constants(self):
+        """ Tests Parser.transform_constants(). """
+
+        # Parse model, which should convert an initial value to an equation
+        model = load_model('initial_value_constant.cellml')
+        v = model.get_variable_by_name('A$a')
+
+        # Initial value should be None, if transform_constants worked
+        assert v.initial_value is None
+
+        # And its RHS value should be 1
+        assert model.get_value(v) == 1
+
+        # And the equations should have matching units
+        for eq in model.equations:
+            check_left_right_units_equal(model.units, eq)
+
     def test_validation(self):
         """ Tests that RNG validation can pick stuff up. """