From 512765235741055f76d172461ba5454a4f7b292f Mon Sep 17 00:00:00 2001
From: Martin Robinson <martinjrobins@gmail.com>
Date: Wed, 4 Mar 2020 10:56:37 +0000
Subject: [PATCH 01/30] #858 update state vector to allow derivative of y or
 ydot

---
 pybamm/expression_tree/state_vector.py | 5 +++++
 1 file changed, 5 insertions(+)

diff --git a/pybamm/expression_tree/state_vector.py b/pybamm/expression_tree/state_vector.py
index d99079b085..fad0f69e92 100644
--- a/pybamm/expression_tree/state_vector.py
+++ b/pybamm/expression_tree/state_vector.py
@@ -118,6 +118,11 @@ def _base_evaluate(self, t=None, y=None, u=None):
             return out
 
     def _jac(self, variable):
+        if isinstance(variable, pybamm.StateVector):
+            return _jac_same_vector(variable)
+        else if isinstance(variable, pybamm.StateVector):
+
+    def _jac_same_vector(self, variable):
         """
         Differentiate a slice of a StateVector of size m with respect to another
         slice of a StateVector of size n. This returns a (sparse) matrix of size

From 4c971a6919a2f9cdeed45d748e2bd05d0c7c5cfe Mon Sep 17 00:00:00 2001
From: Martin Robinson <martinjrobins@gmail.com>
Date: Thu, 5 Mar 2020 17:05:51 +0000
Subject: [PATCH 02/30] #858 add stateVectorDot

---
 pybamm/__init__.py                            |   2 +-
 pybamm/expression_tree/state_vector.py        | 147 +++++++++++++++++-
 .../test_expression_tree/test_state_vector.py |  18 +++
 3 files changed, 158 insertions(+), 9 deletions(-)

diff --git a/pybamm/__init__.py b/pybamm/__init__.py
index da9022e60a..2c3c13bbcd 100644
--- a/pybamm/__init__.py
+++ b/pybamm/__init__.py
@@ -115,7 +115,7 @@ def version(formatted=False):
 )
 from .expression_tree.independent_variable import t
 from .expression_tree.vector import Vector
-from .expression_tree.state_vector import StateVector
+from .expression_tree.state_vector import StateVector, StateVectorDot
 
 from .expression_tree.exceptions import (
     DomainError,
diff --git a/pybamm/expression_tree/state_vector.py b/pybamm/expression_tree/state_vector.py
index fad0f69e92..e8addad924 100644
--- a/pybamm/expression_tree/state_vector.py
+++ b/pybamm/expression_tree/state_vector.py
@@ -7,7 +7,7 @@
 from scipy.sparse import csr_matrix, vstack
 
 
-class StateVector(pybamm.Symbol):
+class StateVectorBase(pybamm.Symbol):
     """
     node in the expression tree that holds a slice to read from an external vector type
 
@@ -32,6 +32,7 @@ class StateVector(pybamm.Symbol):
     def __init__(
         self,
         *y_slices,
+        base_name="y",
         name=None,
         domain=None,
         auxiliary_domains=None,
@@ -42,9 +43,10 @@ def __init__(
                 raise TypeError("all y_slices must be slice objects")
         if name is None:
             if y_slices[0].start is None:
-                name = "y[:{:d}]".format(y_slice.stop)
+                name = base_name + "[:{:d}]".format(y_slice.stop)
             else:
-                name = "y[{:d}:{:d}".format(y_slices[0].start, y_slices[0].stop)
+                name = base_name + \
+                    "[{:d}:{:d}".format(y_slices[0].start, y_slices[0].stop)
             if len(y_slices) > 1:
                 name += ",{:d}:{:d}".format(y_slices[1].start, y_slices[1].stop)
                 if len(y_slices) > 2:
@@ -103,7 +105,7 @@ def set_id(self):
             + tuple(self.domain)
         )
 
-    def _base_evaluate(self, t=None, y=None, u=None):
+    def _base_evaluate(self, t=None, y=None, y_dot=None, u=None):
         """ See :meth:`pybamm.Symbol._base_evaluate()`. """
         if y is None:
             raise TypeError("StateVector cannot evaluate input 'y=None'")
@@ -117,10 +119,27 @@ def _base_evaluate(self, t=None, y=None, u=None):
                 out = out[:, np.newaxis]
             return out
 
-    def _jac(self, variable):
-        if isinstance(variable, pybamm.StateVector):
-            return _jac_same_vector(variable)
-        else if isinstance(variable, pybamm.StateVector):
+    def _jac_diff_vector(self, variable):
+        """
+        Differentiate a slice of a StateVector of size m with respect to another
+        slice of a different StateVector of size n. This returns a (sparse) zero matrix of size
+        m x n
+
+        Parameters
+        ----------
+        variable : :class:`pybamm.Symbol`
+            The variable with respect to which to differentiate
+
+        """
+        if len(variable.y_slices) > 1:
+            raise NotImplementedError(
+                "Jacobian only implemented for a single-slice StateVector"
+            )
+        slices_size = self.y_slices[0].stop - self.y_slices[0].start
+        variable_size = variable.last_point - variable.first_point
+
+        # Return zeros of correct size since no entries match
+        return csr_matrix(slices_size, variable_size)
 
     def _jac_same_vector(self, variable):
         """
@@ -185,3 +204,115 @@ def _evaluate_for_shape(self):
         See :meth:`pybamm.Symbol.evaluate_for_shape()`
         """
         return np.nan * np.ones((self.size, 1))
+
+
+class StateVector(StateVectorBase):
+    """
+    node in the expression tree that holds a slice to read from an external vector type
+
+    Parameters
+    ----------
+
+    y_slice: slice
+        the slice of an external y to read
+    name: str, optional
+        the name of the node
+    domain : iterable of str, optional
+        list of domains the parameter is valid over, defaults to empty list
+    auxiliary_domains : dict of str, optional
+        dictionary of auxiliary domains
+    evaluation_array : list, optional
+        List of boolean arrays representing slices. Default is None, in which case the
+        evaluation_array is computed from y_slices.
+
+    *Extends:* :class:`Array`
+    """
+
+    def __init__(
+        self,
+        *y_slices,
+        name=None,
+        domain=None,
+        auxiliary_domains=None,
+        evaluation_array=None,
+    ):
+        super().__init__(*y_slices,
+                         base_name="y", name=name, domain=domain,
+                         auxiliary_domains=auxiliary_domains,
+                         evaluation_array=evaluation_array)
+
+    def _base_evaluate(self, t=None, y=None, u=None):
+        """ See :meth:`pybamm.Symbol._base_evaluate()`. """
+        if y is None:
+            raise TypeError("StateVector cannot evaluate input 'y=None'")
+        if y.shape[0] < len(self.evaluation_array):
+            raise ValueError(
+                "y is too short, so value with slice is smaller than expected"
+            )
+
+        out = (y[: len(self._evaluation_array)])[self._evaluation_array]
+        if isinstance(out, np.ndarray) and out.ndim == 1:
+            out = out[:, np.newaxis]
+        return out
+
+    def _jac(self, variable):
+        if isinstance(variable, pybamm.StateVector):
+            return self._jac_same_vector(variable)
+        elif isinstance(variable, pybamm.StateVectorDot):
+            return self._jac_diff_vector(variable)
+
+
+class StateVectorDot(StateVectorBase):
+    """
+    node in the expression tree that holds a slice to read from the ydot
+
+    Parameters
+    ----------
+
+    y_slice: slice
+        the slice of an external ydot to read
+    name: str, optional
+        the name of the node
+    domain : iterable of str, optional
+        list of domains the parameter is valid over, defaults to empty list
+    auxiliary_domains : dict of str, optional
+        dictionary of auxiliary domains
+    evaluation_array : list, optional
+        List of boolean arrays representing slices. Default is None, in which case the
+        evaluation_array is computed from y_slices.
+
+    *Extends:* :class:`Array`
+    """
+
+    def __init__(
+        self,
+        *y_slices,
+        name=None,
+        domain=None,
+        auxiliary_domains=None,
+        evaluation_array=None,
+    ):
+        super().__init__(*y_slices,
+                         base_name="y_dot", name=name, domain=domain,
+                         auxiliary_domains=auxiliary_domains,
+                         evaluation_array=evaluation_array)
+
+    def _base_evaluate(self, t=None, y=None, y_dot=None, u=None):
+        """ See :meth:`pybamm.Symbol._base_evaluate()`. """
+        if y_dot is None:
+            raise TypeError("StateVectorDot cannot evaluate input 'y_dot=None'")
+        if y_dot.shape[0] < len(self.evaluation_array):
+            raise ValueError(
+                "y_dot is too short, so value with slice is smaller than expected"
+            )
+
+        out = (y_dot[: len(self._evaluation_array)])[self._evaluation_array]
+        if isinstance(out, np.ndarray) and out.ndim == 1:
+            out = out[:, np.newaxis]
+        return out
+
+    def _jac(self, variable):
+        if isinstance(variable, pybamm.StateVectorDot):
+            return self._jac_same_vector(variable)
+        elif isinstance(variable, pybamm.StateVector):
+            return self._jac_diff_vector(variable)
diff --git a/tests/unit/test_expression_tree/test_state_vector.py b/tests/unit/test_expression_tree/test_state_vector.py
index f6809565d3..3038c4c4fa 100644
--- a/tests/unit/test_expression_tree/test_state_vector.py
+++ b/tests/unit/test_expression_tree/test_state_vector.py
@@ -61,6 +61,24 @@ def test_failure(self):
         with self.assertRaisesRegex(TypeError, "all y_slices must be slice objects"):
             pybamm.StateVector(slice(0, 10), 1)
 
+class TestStateVectorDot(unittest.TestCase):
+    def test_evaluate(self):
+        sv = pybamm.StateVectorDot(slice(0, 10))
+        y_dot = np.linspace(0, 2, 19)
+        np.testing.assert_array_equal(
+            sv.evaluate(y_dot=y_dot), np.linspace(0, 1, 10)[:, np.newaxis]
+        )
+
+        # Try evaluating with a y that is too short
+        y_dot2 = np.ones(5)
+        with self.assertRaisesRegex(
+            ValueError, "y_dot is too short, so value with slice is smaller than expected"
+        ):
+            sv.evaluate(y_dot=y_dot2)
+
+    def test_name(self):
+        sv = pybamm.StateVectorDot(slice(0, 10))
+        self.assertEqual(sv.name, "y_dot[0:10]")
 
 if __name__ == "__main__":
     print("Add -v for more debug output")

From 96c6a417d4b9c7cdfd5115c71ec06fc51ba901ce Mon Sep 17 00:00:00 2001
From: Martin Robinson <martinjrobins@gmail.com>
Date: Fri, 6 Mar 2020 05:29:18 +0000
Subject: [PATCH 03/30] #858 add y_dot arg to all evaluate functions

---
 pybamm/expression_tree/array.py                |  2 +-
 pybamm/expression_tree/independent_variable.py |  2 +-
 pybamm/expression_tree/input_parameter.py      |  2 +-
 pybamm/expression_tree/scalar.py               |  2 +-
 pybamm/expression_tree/state_vector.py         |  2 +-
 pybamm/expression_tree/symbol.py               | 17 +++++++++++------
 pybamm/expression_tree/variable.py             |  2 +-
 7 files changed, 17 insertions(+), 12 deletions(-)

diff --git a/pybamm/expression_tree/array.py b/pybamm/expression_tree/array.py
index 348fbed616..43a8864eb7 100644
--- a/pybamm/expression_tree/array.py
+++ b/pybamm/expression_tree/array.py
@@ -98,6 +98,6 @@ def new_copy(self):
             self.entries_string,
         )
 
-    def _base_evaluate(self, t=None, y=None, u=None):
+    def _base_evaluate(self, t=None, y=None, y_dot=None, u=None):
         """ See :meth:`pybamm.Symbol._base_evaluate()`. """
         return self._entries
diff --git a/pybamm/expression_tree/independent_variable.py b/pybamm/expression_tree/independent_variable.py
index 9dea22f5a7..f04d984a53 100644
--- a/pybamm/expression_tree/independent_variable.py
+++ b/pybamm/expression_tree/independent_variable.py
@@ -51,7 +51,7 @@ def new_copy(self):
         """ See :meth:`pybamm.Symbol.new_copy()`. """
         return Time()
 
-    def _base_evaluate(self, t, y=None, u=None):
+    def _base_evaluate(self, t, y=None, y_dot=None, u=None):
         """ See :meth:`pybamm.Symbol._base_evaluate()`. """
         if t is None:
             raise ValueError("t must be provided")
diff --git a/pybamm/expression_tree/input_parameter.py b/pybamm/expression_tree/input_parameter.py
index 148c062b0c..395b84ebb2 100644
--- a/pybamm/expression_tree/input_parameter.py
+++ b/pybamm/expression_tree/input_parameter.py
@@ -36,7 +36,7 @@ def _jac(self, variable):
         """ See :meth:`pybamm.Symbol._jac()`. """
         return pybamm.Scalar(0)
 
-    def _base_evaluate(self, t=None, y=None, u=None):
+    def _base_evaluate(self, t=None, y=None, y_dot=None, u=None):
         # u should be a dictionary
         # convert 'None' to empty dictionary for more informative error
         if u is None:
diff --git a/pybamm/expression_tree/scalar.py b/pybamm/expression_tree/scalar.py
index b96d618c70..3b3d1225f5 100644
--- a/pybamm/expression_tree/scalar.py
+++ b/pybamm/expression_tree/scalar.py
@@ -51,7 +51,7 @@ def set_id(self):
             (self.__class__, self.name) + tuple(self.domain) + tuple(str(self._value))
         )
 
-    def _base_evaluate(self, t=None, y=None, u=None):
+    def _base_evaluate(self, t=None, y=None, y_dot=None, u=None):
         """ See :meth:`pybamm.Symbol._base_evaluate()`. """
         return self._value
 
diff --git a/pybamm/expression_tree/state_vector.py b/pybamm/expression_tree/state_vector.py
index e8addad924..7968c36320 100644
--- a/pybamm/expression_tree/state_vector.py
+++ b/pybamm/expression_tree/state_vector.py
@@ -241,7 +241,7 @@ def __init__(
                          auxiliary_domains=auxiliary_domains,
                          evaluation_array=evaluation_array)
 
-    def _base_evaluate(self, t=None, y=None, u=None):
+    def _base_evaluate(self, t=None, y=None, y_dot=None, u=None):
         """ See :meth:`pybamm.Symbol._base_evaluate()`. """
         if y is None:
             raise TypeError("StateVector cannot evaluate input 'y=None'")
diff --git a/pybamm/expression_tree/symbol.py b/pybamm/expression_tree/symbol.py
index 3bc390cc26..0c98cb07a6 100644
--- a/pybamm/expression_tree/symbol.py
+++ b/pybamm/expression_tree/symbol.py
@@ -506,7 +506,7 @@ def _jac(self, variable):
         """
         raise NotImplementedError
 
-    def _base_evaluate(self, t=None, y=None, u=None):
+    def _base_evaluate(self, t=None, y=None, y_dot=None, u=None):
         """evaluate expression tree
 
         will raise a ``NotImplementedError`` if this member function has not
@@ -520,7 +520,10 @@ def _base_evaluate(self, t=None, y=None, u=None):
             time at which to evaluate (default None)
 
         y : numpy.array, optional
-            array to evaluate when solving (default None)
+            array with state values to evaluate when solving (default None)
+
+        y_dot : numpy.array, optional
+            array with time derivatives of state values to evaluate when solving (default None)
 
         """
         raise NotImplementedError(
@@ -530,7 +533,7 @@ def _base_evaluate(self, t=None, y=None, u=None):
             )
         )
 
-    def evaluate(self, t=None, y=None, u=None, known_evals=None):
+    def evaluate(self, t=None, y=None, y_dot=None, u=None, known_evals=None):
         """Evaluate expression tree (wrapper to allow using dict of known values).
         If the dict 'known_evals' is provided, the dict is searched for self.id; if
         self.id is in the keys, return that value; otherwise, evaluate using
@@ -541,7 +544,9 @@ def evaluate(self, t=None, y=None, u=None, known_evals=None):
         t : float or numeric type, optional
             time at which to evaluate (default None)
         y : numpy.array, optional
-            array to evaluate when solving (default None)
+            array with state values to evaluate when solving (default None)
+        y_dot : numpy.array, optional
+            array with time derivatives of state values to evaluate when solving (default None)
         u : dict, optional
             dictionary of inputs to use when solving (default None)
         known_evals : dict, optional
@@ -556,10 +561,10 @@ def evaluate(self, t=None, y=None, u=None, known_evals=None):
         """
         if known_evals is not None:
             if self.id not in known_evals:
-                known_evals[self.id] = self._base_evaluate(t, y, u)
+                known_evals[self.id] = self._base_evaluate(t, y, y_dot, u)
             return known_evals[self.id], known_evals
         else:
-            return self._base_evaluate(t, y, u)
+            return self._base_evaluate(t, y, y_dot, u)
 
     def evaluate_for_shape(self):
         """Evaluate expression tree to find its shape. For symbols that cannot be
diff --git a/pybamm/expression_tree/variable.py b/pybamm/expression_tree/variable.py
index 935305c924..bfdf02a200 100644
--- a/pybamm/expression_tree/variable.py
+++ b/pybamm/expression_tree/variable.py
@@ -84,7 +84,7 @@ def _evaluate_for_shape(self):
         """ See :meth:`pybamm.Symbol.evaluate_for_shape_using_domain()` """
         return np.nan * np.ones((self.size, 1))
 
-    def _base_evaluate(self, t=None, y=None, u=None):
+    def _base_evaluate(self, t=None, y=None, y_dot=None, u=None):
         # u should be a dictionary
         # convert 'None' to empty dictionary for more informative error
         if u is None:

From 5489aa3d4112547f68d5b4ac306e8477151fcff9 Mon Sep 17 00:00:00 2001
From: Martin Robinson <martinjrobins@gmail.com>
Date: Fri, 6 Mar 2020 05:37:10 +0000
Subject: [PATCH 04/30] #858 add VariableDot

---
 pybamm/__init__.py                            |  2 +-
 pybamm/expression_tree/variable.py            | 32 ++++++++++++++++++-
 .../test_expression_tree/test_variable.py     | 18 +++++++++++
 3 files changed, 50 insertions(+), 2 deletions(-)

diff --git a/pybamm/__init__.py b/pybamm/__init__.py
index 2c3c13bbcd..47ab5bfe21 100644
--- a/pybamm/__init__.py
+++ b/pybamm/__init__.py
@@ -107,7 +107,7 @@ def version(formatted=False):
     ones_like,
 )
 from .expression_tree.scalar import Scalar
-from .expression_tree.variable import Variable, ExternalVariable
+from .expression_tree.variable import Variable, ExternalVariable, VariableDot
 from .expression_tree.independent_variable import (
     IndependentVariable,
     Time,
diff --git a/pybamm/expression_tree/variable.py b/pybamm/expression_tree/variable.py
index bfdf02a200..db4313171a 100644
--- a/pybamm/expression_tree/variable.py
+++ b/pybamm/expression_tree/variable.py
@@ -48,8 +48,38 @@ def _evaluate_for_shape(self):
         )
 
 
+class VariableDot(Variable):
+    """
+    A node in the expression tree represending the time derviative of a dependent
+    variable
+
+    This node will be discretised by :class:`.Discretisation` and converted
+    to a :class:`pybamm.StateVectorDot` node.
+
+    Parameters
+    ----------
+
+    name : str
+        name of the node
+    domain : iterable of str
+        list of domains that this variable is valid over
+    auxiliary_domains : dict
+        dictionary of auxiliary domains ({'secondary': ..., 'tertiary': ...}). For
+        example, for the single particle model, the particle concentration would be a
+        Variable with domain 'negative particle' and secondary auxiliary domain 'current
+        collector'. For the DFN, the particle concentration would be a Variable with
+        domain 'negative particle', secondary domain 'negative electrode' and tertiary
+        domain 'current collector'
+
+    *Extends:* :class:`Symbol`
+    """
+
+    def __init__(self, name, domain=None, auxiliary_domains=None):
+        super().__init__(name, domain=domain, auxiliary_domains=auxiliary_domains)
+
+
 class ExternalVariable(Variable):
-    """A node in the expression tree represending an external variable variable
+    """A node in the expression tree representing an external variable variable
 
     This node will be discretised by :class:`.Discretisation` and converted
     to a :class:`.Vector` node.
diff --git a/tests/unit/test_expression_tree/test_variable.py b/tests/unit/test_expression_tree/test_variable.py
index 1e4e1b0bee..47ae081ffb 100644
--- a/tests/unit/test_expression_tree/test_variable.py
+++ b/tests/unit/test_expression_tree/test_variable.py
@@ -25,6 +25,24 @@ def test_variable_id(self):
         self.assertNotEqual(a1.id, a3.id)
         self.assertNotEqual(a1.id, a4.id)
 
+class TestVariableDot(unittest.TestCase):
+    def test_variable_init(self):
+        a = pybamm.VariableDot("a'")
+        self.assertEqual(a.name, "a'")
+        self.assertEqual(a.domain, [])
+        a = pybamm.VariableDot("a", domain=["test"])
+        self.assertEqual(a.domain[0], "test")
+        self.assertRaises(TypeError, pybamm.Variable("a", domain="test"))
+
+    def test_variable_id(self):
+        a1 = pybamm.VariableDot("a", domain=["negative electrode"])
+        a2 = pybamm.VariableDot("a", domain=["negative electrode"])
+        self.assertEqual(a1.id, a2.id)
+        a3 = pybamm.VariableDot("b", domain=["negative electrode"])
+        a4 = pybamm.VariableDot("a", domain=["positive electrode"])
+        self.assertNotEqual(a1.id, a3.id)
+        self.assertNotEqual(a1.id, a4.id)
+
 
 class TestExternalVariable(unittest.TestCase):
     def test_external_variable_scalar(self):

From e2b4ba7855b3259fa0e9004d115d5f9ad15a84d6 Mon Sep 17 00:00:00 2001
From: Martin Robinson <martinjrobins@gmail.com>
Date: Fri, 6 Mar 2020 06:43:12 +0000
Subject: [PATCH 05/30] #858 add d_dt unary operator, fixes to jac to support
 this

---
 pybamm/expression_tree/binary_operators.py    | 14 +++++-----
 .../expression_tree/independent_variable.py   |  6 +++--
 pybamm/expression_tree/operations/jacobian.py | 21 ++++++++++++---
 pybamm/expression_tree/symbol.py              | 15 ++++++++---
 pybamm/expression_tree/unary_operators.py     | 26 +++++++++++++++++++
 pybamm/expression_tree/variable.py            | 10 +++++++
 6 files changed, 76 insertions(+), 16 deletions(-)

diff --git a/pybamm/expression_tree/binary_operators.py b/pybamm/expression_tree/binary_operators.py
index 65c459926b..beb322faa8 100644
--- a/pybamm/expression_tree/binary_operators.py
+++ b/pybamm/expression_tree/binary_operators.py
@@ -32,7 +32,7 @@ def is_matrix_zero(expr):
         return False
 
 
-def is_one(expr):
+def is_scalar_one(expr):
     """
     Utility function to test if an expression evaluates to a constant scalar one
     """
@@ -253,7 +253,7 @@ def _binary_simplify(self, left, right):
             return pybamm.Scalar(1)
 
         # anything to the power of one is itself
-        if is_scalar_zero(left):
+        if is_scalar_one(right):
             return left
 
         return self.__class__(left, right)
@@ -425,9 +425,9 @@ def _binary_simplify(self, left, right):
             return zeros_of_shape(shape)
 
         # anything multiplied by a scalar one returns itself
-        if is_one(left):
+        if is_scalar_one(left):
             return right
-        if is_one(right):
+        if is_scalar_one(right):
             return left
 
         return pybamm.simplify_multiplication_division(self.__class__, left, right)
@@ -549,7 +549,7 @@ def _binary_simplify(self, left, right):
                 return pybamm.Array(np.inf * np.ones(left.shape_for_testing))
 
         # anything divided by one is itself
-        if is_one(right):
+        if is_scalar_one(right):
             return left
 
         return pybamm.simplify_multiplication_division(self.__class__, left, right)
@@ -622,9 +622,9 @@ def _binary_simplify(self, left, right):
             return zeros_of_shape(shape)
 
         # anything multiplied by a scalar one returns itself
-        if is_one(left):
+        if is_scalar_one(left):
             return right
-        if is_one(right):
+        if is_scalar_one(right):
             return left
 
         return pybamm.simplify_multiplication_division(self.__class__, left, right)
diff --git a/pybamm/expression_tree/independent_variable.py b/pybamm/expression_tree/independent_variable.py
index f04d984a53..3bc842eebf 100644
--- a/pybamm/expression_tree/independent_variable.py
+++ b/pybamm/expression_tree/independent_variable.py
@@ -35,8 +35,10 @@ def _evaluate_for_shape(self):
 
     def _jac(self, variable):
         """ See :meth:`pybamm.Symbol._jac()`. """
-        return pybamm.Scalar(0)
-
+        if variable.id == self.id:
+            return pybamm.Scalar(1)
+        else:
+            return pybamm.Scalar(0)
 
 class Time(IndependentVariable):
     """A node in the expression tree representing time
diff --git a/pybamm/expression_tree/operations/jacobian.py b/pybamm/expression_tree/operations/jacobian.py
index 42bd3683da..b231c0c784 100644
--- a/pybamm/expression_tree/operations/jacobian.py
+++ b/pybamm/expression_tree/operations/jacobian.py
@@ -5,8 +5,22 @@
 
 
 class Jacobian(object):
-    def __init__(self, known_jacs=None):
+    """
+    Helper class to calculate the jacobian of an expression.
+
+    Parameters
+    ----------
+
+    known_jacs: dict {variable ids -> :class:`pybamm.Symbol`}
+        cached jacobians
+
+    clear_domain: bool
+        wether or not the jacobian clears the domain (default True)
+    """
+
+    def __init__(self, known_jacs=None, clear_domain=True):
         self._known_jacs = known_jacs or {}
+        self._clear_domain = clear_domain
 
     def jac(self, symbol, variable):
         """
@@ -75,6 +89,7 @@ def _jac(self, symbol, variable):
                     )
                 )
 
-        # jacobian removes the domain(s)
-        jac.clear_domains()
+        # jacobian by default removes the domain(s)
+        if self._clear_domain:
+            jac.clear_domains()
         return jac
diff --git a/pybamm/expression_tree/symbol.py b/pybamm/expression_tree/symbol.py
index 0c98cb07a6..f09008d1c0 100644
--- a/pybamm/expression_tree/symbol.py
+++ b/pybamm/expression_tree/symbol.py
@@ -492,12 +492,13 @@ def _diff(self, variable):
         "Default behaviour for differentiation, overriden by Binary and Unary Operators"
         raise NotImplementedError
 
-    def jac(self, variable, known_jacs=None):
+    def jac(self, variable, known_jacs=None, clear_domain=True):
         """
         Differentiate a symbol with respect to a (slice of) a State Vector.
         See :class:`pybamm.Jacobian`.
         """
-        return pybamm.Jacobian(known_jacs).jac(self, variable)
+        jac = pybamm.Jacobian(known_jacs, clear_domain=clear_domain)
+        return jac.jac(self, variable)
 
     def _jac(self, variable):
         """
@@ -606,7 +607,7 @@ def is_constant(self):
     def evaluate_ignoring_errors(self):
         """
         Evaluates the expression. If a node exists in the tree that cannot be evaluated
-        as a scalar or vector (e.g. Parameter, Variable, StateVector, InputParameter),
+        as a scalar or vector (e.g. Time, Parameter, Variable, StateVector, InputParameter),
         then None is returned. Otherwise the result of the evaluation is given
 
         See Also
@@ -615,7 +616,7 @@ def evaluate_ignoring_errors(self):
 
         """
         try:
-            result = self.evaluate(t=0, u="shape test")
+            result = self.evaluate(u="shape test")
         except NotImplementedError:
             # return None if NotImplementedError is raised
             # (there is a e.g. Parameter, Variable, ... in the tree)
@@ -628,6 +629,10 @@ def evaluate_ignoring_errors(self):
             else:
                 raise error
         except ValueError as e:
+            # return None if specific ValueError is raised
+            # (there is a e.g. Time in the tree)
+            if e.args[0] == "t must be provided":
+                return None
             raise pybamm.ShapeError("Cannot find shape (original error: {})".format(e))
         return result
 
@@ -762,3 +767,5 @@ def test_shape(self):
             self.shape_for_testing
         except ValueError as e:
             raise pybamm.ShapeError("Cannot find shape (original error: {})".format(e))
+
+
diff --git a/pybamm/expression_tree/unary_operators.py b/pybamm/expression_tree/unary_operators.py
index 0b4f68e22d..ae3a1af9a7 100644
--- a/pybamm/expression_tree/unary_operators.py
+++ b/pybamm/expression_tree/unary_operators.py
@@ -763,6 +763,32 @@ def __init__(self, child, side):
         super().__init__("boundary flux", child, side)
 
 
+
+def d_dt(expression):
+    """
+    convenience function for taking the time derivative of an expression
+
+    Note that this operator is different to the other unary operators in that it is
+    *not* lazily evaluated, it instead returns the expression tree that is the time
+    derivative of the input
+
+    Parameters
+    ----------
+
+    expression : :class:`Symbol`
+        the time derivative will be performed on this sub-expression
+
+    Returns
+    -------
+
+    :class:`Symbol`
+        the time derivative of ``expression``
+    """
+
+    return expression.jac(pybamm.t, clear_domain=False)
+
+
+
 #
 # Methods to call Gradient, Divergence, Laplacian and Gradient_Squared
 #
diff --git a/pybamm/expression_tree/variable.py b/pybamm/expression_tree/variable.py
index db4313171a..5a081c0778 100644
--- a/pybamm/expression_tree/variable.py
+++ b/pybamm/expression_tree/variable.py
@@ -47,6 +47,16 @@ def _evaluate_for_shape(self):
             self.domain, self.auxiliary_domains
         )
 
+    def _jac(self, variable):
+        if variable == self:
+            return pybamm.Scalar(1)
+        elif variable == pybamm.t:
+            return pybamm.VariableDot(self.name+"'",
+                                      domain=self.domain,
+                                      auxiliary_domains=self.auxiliary_domains)
+        else:
+            return pybamm.Scalar(0)
+
 
 class VariableDot(Variable):
     """

From 130bbb98da3d716236dd2f87ac3dbd4853804878 Mon Sep 17 00:00:00 2001
From: Martin Robinson <martinjrobins@gmail.com>
Date: Fri, 6 Mar 2020 07:04:24 +0000
Subject: [PATCH 06/30] #858 discretisation of VariableDot results in
 StateVectorDot

---
 pybamm/discretisations/discretisation.py             |  7 +++++++
 pybamm/expression_tree/variable.py                   | 12 ++++++++++++
 .../unit/test_discretisations/test_discretisation.py |  9 ++++++++-
 3 files changed, 27 insertions(+), 1 deletion(-)

diff --git a/pybamm/discretisations/discretisation.py b/pybamm/discretisations/discretisation.py
index b74ac08a72..80ef5dc8fe 100644
--- a/pybamm/discretisations/discretisation.py
+++ b/pybamm/discretisations/discretisation.py
@@ -855,6 +855,13 @@ def _process_symbol(self, symbol):
             disc_children = [self.process_symbol(child) for child in symbol.children]
             return symbol._function_new_copy(disc_children)
 
+        elif isinstance(symbol, pybamm.VariableDot):
+            return pybamm.StateVectorDot(
+                *self.y_slices[symbol.get_variable().id],
+                domain=symbol.domain,
+                auxiliary_domains=symbol.auxiliary_domains
+            )
+
         elif isinstance(symbol, pybamm.Variable):
             # Check if variable is a standard variable or an external variable
             if any(symbol.id == var.id for var in self.external_variables.values()):
diff --git a/pybamm/expression_tree/variable.py b/pybamm/expression_tree/variable.py
index 5a081c0778..908a86aa53 100644
--- a/pybamm/expression_tree/variable.py
+++ b/pybamm/expression_tree/variable.py
@@ -87,6 +87,18 @@ class VariableDot(Variable):
     def __init__(self, name, domain=None, auxiliary_domains=None):
         super().__init__(name, domain=domain, auxiliary_domains=auxiliary_domains)
 
+    def get_variable(self):
+        """
+        return a :class:`.Variable` corresponding to this VariableDot
+
+        Note: Variable._jac adds a dash to the name of the corresponding VariableDot, so
+        we remove this here
+
+        """
+        return Variable(self.name[:-1],
+                        domain=self._domain,
+                        auxiliary_domains=self._auxiliary_domains)
+
 
 class ExternalVariable(Variable):
     """A node in the expression tree representing an external variable variable
diff --git a/tests/unit/test_discretisations/test_discretisation.py b/tests/unit/test_discretisations/test_discretisation.py
index d30b6f11f5..17122a3c74 100644
--- a/tests/unit/test_discretisations/test_discretisation.py
+++ b/tests/unit/test_discretisations/test_discretisation.py
@@ -329,6 +329,14 @@ def test_process_symbol_base(self):
         var_disc = disc.process_symbol(var)
         self.assertIsInstance(var_disc, pybamm.StateVector)
         self.assertEqual(var_disc.y_slices[0], disc.y_slices[var.id][0])
+
+        # variable dot
+        var_dot = pybamm.VariableDot("var'")
+        var_vec_dot = pybamm.VariableDot("var vec'", domain=["negative electrode"])
+        var_dot_disc = disc.process_symbol(var_dot)
+        self.assertIsInstance(var_dot_disc, pybamm.StateVectorDot)
+        self.assertEqual(var_dot_disc.y_slices[0], disc.y_slices[var.id][0])
+
         # scalar
         scal = pybamm.Scalar(5)
         scal_disc = disc.process_symbol(scal)
@@ -1086,7 +1094,6 @@ def test_mass_matirx_inverse(self):
             model.mass_matrix_inv.entries.toarray(), mass_inv.toarray()
         )
 
-
 if __name__ == "__main__":
     print("Add -v for more debug output")
     import sys

From 1f03394d217ac0f24f9328ea2cba70113a3ee8d8 Mon Sep 17 00:00:00 2001
From: Martin Robinson <martinjrobins@gmail.com>
Date: Fri, 6 Mar 2020 09:56:55 +0000
Subject: [PATCH 07/30] #858 time derivative of state vector gives state vector
 dot, raise errors if taking time derivative of *Dot classes

---
 pybamm/expression_tree/state_vector.py | 11 ++++++++++-
 pybamm/expression_tree/variable.py     | 22 ++++++++++++++++++++--
 2 files changed, 30 insertions(+), 3 deletions(-)

diff --git a/pybamm/expression_tree/state_vector.py b/pybamm/expression_tree/state_vector.py
index 7968c36320..259f19336c 100644
--- a/pybamm/expression_tree/state_vector.py
+++ b/pybamm/expression_tree/state_vector.py
@@ -256,7 +256,12 @@ def _base_evaluate(self, t=None, y=None, y_dot=None, u=None):
         return out
 
     def _jac(self, variable):
-        if isinstance(variable, pybamm.StateVector):
+        if variable.id == pybamm.t.id:
+            return StateVectorDot(*self._y_slices, name=self.name + "'",
+                                  domain=self.domain,
+                                  auxiliary_domains=self.auxiliary_domains,
+                                  evaluation_array=self.evaluation_array)
+        elif isinstance(variable, pybamm.StateVector):
             return self._jac_same_vector(variable)
         elif isinstance(variable, pybamm.StateVectorDot):
             return self._jac_diff_vector(variable)
@@ -312,6 +317,10 @@ def _base_evaluate(self, t=None, y=None, y_dot=None, u=None):
         return out
 
     def _jac(self, variable):
+        if variable.id == pybamm.t.id:
+            raise pybamm.ModelError(
+                "cannot take second time derivative of a state vector"
+            )
         if isinstance(variable, pybamm.StateVectorDot):
             return self._jac_same_vector(variable)
         elif isinstance(variable, pybamm.StateVector):
diff --git a/pybamm/expression_tree/variable.py b/pybamm/expression_tree/variable.py
index 908a86aa53..7a9b9cd618 100644
--- a/pybamm/expression_tree/variable.py
+++ b/pybamm/expression_tree/variable.py
@@ -48,9 +48,9 @@ def _evaluate_for_shape(self):
         )
 
     def _jac(self, variable):
-        if variable == self:
+        if variable.id == self.id:
             return pybamm.Scalar(1)
-        elif variable == pybamm.t:
+        elif variable.id == pybamm.t.id:
             return pybamm.VariableDot(self.name+"'",
                                       domain=self.domain,
                                       auxiliary_domains=self.auxiliary_domains)
@@ -99,6 +99,14 @@ def get_variable(self):
                         domain=self._domain,
                         auxiliary_domains=self._auxiliary_domains)
 
+    def _jac(self, variable):
+        if variable.id == self.id:
+            return pybamm.Scalar(1)
+        elif variable.id == pybamm.t.id:
+            raise pybamm.ModelError("cannot take second time derivative of a Variable")
+        else:
+            return pybamm.Scalar(0)
+
 
 class ExternalVariable(Variable):
     """A node in the expression tree representing an external variable variable
@@ -161,3 +169,13 @@ def _base_evaluate(self, t=None, y=None, y_dot=None, u=None):
         # raise more informative error if can't find name in dict
         except KeyError:
             raise KeyError("External variable '{}' not found".format(self.name))
+
+    def _jac(self, variable):
+        if variable.id == self.id:
+            return pybamm.Scalar(1)
+        elif variable.id == pybamm.t.id:
+            raise pybamm.ModelError("cannot take time derivative of an external variable")
+        else:
+            return pybamm.Scalar(0)
+
+

From 4ce9dbdfebc36b9f5b35c212d58e6cd0b347e940 Mon Sep 17 00:00:00 2001
From: Martin Robinson <martinjrobins@gmail.com>
Date: Fri, 6 Mar 2020 10:16:37 +0000
Subject: [PATCH 08/30] #858 add check for time derivs in model.rhs, fix some
 evaluate arg bugs

---
 pybamm/discretisations/discretisation.py       | 11 +++++++++++
 pybamm/expression_tree/binary_operators.py     | 10 +++++-----
 pybamm/expression_tree/concatenations.py       |  6 +++---
 pybamm/expression_tree/functions.py            |  6 +++---
 pybamm/expression_tree/unary_operators.py      |  6 +++---
 .../test_discretisation.py                     | 18 ++++++++++++++++++
 6 files changed, 43 insertions(+), 14 deletions(-)

diff --git a/pybamm/discretisations/discretisation.py b/pybamm/discretisations/discretisation.py
index 80ef5dc8fe..efbfe46772 100644
--- a/pybamm/discretisations/discretisation.py
+++ b/pybamm/discretisations/discretisation.py
@@ -514,6 +514,7 @@ def process_rhs_and_algebraic(self, model):
             equations) and processed_concatenated_algebraic
 
         """
+
         # Discretise right-hand sides, passing domain from variable
         processed_rhs = self.process_dict(model.rhs)
 
@@ -996,6 +997,16 @@ def check_model(self, model):
         self.check_initial_conditions(model)
         self.check_initial_conditions_rhs(model)
         self.check_variables(model)
+        self.check_for_time_derivatives_in_rhs(model)
+
+    def check_for_time_derivatives_in_rhs(self, model):
+        # Check that no variable time derivatives exist in the rhs equations
+        for eq in model.rhs.values():
+            for node in eq.pre_order():
+                if isinstance(node, pybamm.VariableDot):
+                    raise pybamm.ModelError(
+                        "time derivative of variable found ({}) in rhs equation"
+                    )
 
     def check_initial_conditions(self, model):
         """Check initial conditions are a numpy array"""
diff --git a/pybamm/expression_tree/binary_operators.py b/pybamm/expression_tree/binary_operators.py
index beb322faa8..90ea3c41c0 100644
--- a/pybamm/expression_tree/binary_operators.py
+++ b/pybamm/expression_tree/binary_operators.py
@@ -162,21 +162,21 @@ def _binary_new_copy(self, left, right):
         "Default behaviour for new_copy"
         return self.__class__(left, right)
 
-    def evaluate(self, t=None, y=None, u=None, known_evals=None):
+    def evaluate(self, t=None, y=None, y_dot=None, u=None, known_evals=None):
         """ See :meth:`pybamm.Symbol.evaluate()`. """
         if known_evals is not None:
             id = self.id
             try:
                 return known_evals[id], known_evals
             except KeyError:
-                left, known_evals = self.left.evaluate(t, y, u, known_evals)
-                right, known_evals = self.right.evaluate(t, y, u, known_evals)
+                left, known_evals = self.left.evaluate(t, y, y_dot, u, known_evals)
+                right, known_evals = self.right.evaluate(t, y, y_dot, u, known_evals)
                 value = self._binary_evaluate(left, right)
                 known_evals[id] = value
                 return value, known_evals
         else:
-            left = self.left.evaluate(t, y, u)
-            right = self.right.evaluate(t, y, u)
+            left = self.left.evaluate(t, y, y_dot, u)
+            right = self.right.evaluate(t, y, y_dot, u)
             return self._binary_evaluate(left, right)
 
     def _evaluate_for_shape(self):
diff --git a/pybamm/expression_tree/concatenations.py b/pybamm/expression_tree/concatenations.py
index e1225bf03a..4a74de84cf 100644
--- a/pybamm/expression_tree/concatenations.py
+++ b/pybamm/expression_tree/concatenations.py
@@ -54,7 +54,7 @@ def _concatenation_evaluate(self, children_eval):
         else:
             return self.concatenation_function(children_eval)
 
-    def evaluate(self, t=None, y=None, u=None, known_evals=None):
+    def evaluate(self, t=None, y=None, y_dot=None, u=None, known_evals=None):
         """ See :meth:`pybamm.Symbol.evaluate()`. """
         children = self.cached_children
         if known_evals is not None:
@@ -62,14 +62,14 @@ def evaluate(self, t=None, y=None, u=None, known_evals=None):
                 children_eval = [None] * len(children)
                 for idx, child in enumerate(children):
                     children_eval[idx], known_evals = child.evaluate(
-                        t, y, u, known_evals
+                        t, y, y_dot, u, known_evals
                     )
                 known_evals[self.id] = self._concatenation_evaluate(children_eval)
             return known_evals[self.id], known_evals
         else:
             children_eval = [None] * len(children)
             for idx, child in enumerate(children):
-                children_eval[idx] = child.evaluate(t, y, u)
+                children_eval[idx] = child.evaluate(t, y, y_dot, u)
             return self._concatenation_evaluate(children_eval)
 
     def new_copy(self):
diff --git a/pybamm/expression_tree/functions.py b/pybamm/expression_tree/functions.py
index 2601034b98..c21b473ef4 100644
--- a/pybamm/expression_tree/functions.py
+++ b/pybamm/expression_tree/functions.py
@@ -152,19 +152,19 @@ def _function_jac(self, children_jacs):
 
         return jacobian
 
-    def evaluate(self, t=None, y=None, u=None, known_evals=None):
+    def evaluate(self, t=None, y=None, y_dot=None, u=None, known_evals=None):
         """ See :meth:`pybamm.Symbol.evaluate()`. """
         if known_evals is not None:
             if self.id not in known_evals:
                 evaluated_children = [None] * len(self.children)
                 for i, child in enumerate(self.children):
                     evaluated_children[i], known_evals = child.evaluate(
-                        t, y, u, known_evals=known_evals
+                        t, y, y_dot, u, known_evals=known_evals
                     )
                 known_evals[self.id] = self._function_evaluate(evaluated_children)
             return known_evals[self.id], known_evals
         else:
-            evaluated_children = [child.evaluate(t, y, u) for child in self.children]
+            evaluated_children = [child.evaluate(t, y, y_dot, u) for child in self.children]
             return self._function_evaluate(evaluated_children)
 
     def _evaluate_for_shape(self):
diff --git a/pybamm/expression_tree/unary_operators.py b/pybamm/expression_tree/unary_operators.py
index ae3a1af9a7..a4c7f8574b 100644
--- a/pybamm/expression_tree/unary_operators.py
+++ b/pybamm/expression_tree/unary_operators.py
@@ -63,15 +63,15 @@ def _unary_evaluate(self, child):
         """Perform unary operation on a child. """
         raise NotImplementedError
 
-    def evaluate(self, t=None, y=None, u=None, known_evals=None):
+    def evaluate(self, t=None, y=None, y_dot=None, u=None, known_evals=None):
         """ See :meth:`pybamm.Symbol.evaluate()`. """
         if known_evals is not None:
             if self.id not in known_evals:
-                child, known_evals = self.child.evaluate(t, y, u, known_evals)
+                child, known_evals = self.child.evaluate(t, y, y_dot, u, known_evals)
                 known_evals[self.id] = self._unary_evaluate(child)
             return known_evals[self.id], known_evals
         else:
-            child = self.child.evaluate(t, y, u)
+            child = self.child.evaluate(t, y, y_dot, u)
             return self._unary_evaluate(child)
 
     def _evaluate_for_shape(self):
diff --git a/tests/unit/test_discretisations/test_discretisation.py b/tests/unit/test_discretisations/test_discretisation.py
index 17122a3c74..28df731ec9 100644
--- a/tests/unit/test_discretisations/test_discretisation.py
+++ b/tests/unit/test_discretisations/test_discretisation.py
@@ -701,6 +701,24 @@ def test_process_model_ode(self):
         with self.assertRaises(pybamm.ModelError):
             disc.process_model(model)
 
+        # test that ill possed model with time derivatives of variables in rhs raises an
+        # error
+        model = pybamm.BaseModel()
+        model.rhs = {c: pybamm.div(N) + pybamm.d_dt(c), T: pybamm.div(q), S: pybamm.div(p)}
+        model.initial_conditions = {
+            c: pybamm.Scalar(2),
+            T: pybamm.Scalar(5),
+            S: pybamm.Scalar(8),
+        }
+        model.boundary_conditions = {
+            c: {"left": (0, "Neumann"), "right": (0, "Neumann")},
+            T: {"left": (0, "Neumann"), "right": (0, "Neumann")},
+            S: {"left": (0, "Neumann"), "right": (0, "Neumann")},
+        }
+        model.variables = {"ST": S * T}
+        with self.assertRaises(pybamm.ModelError):
+            disc.process_model(model)
+
     def test_process_model_dae(self):
         # one rhs equation and one algebraic
         whole_cell = ["negative electrode", "separator", "positive electrode"]

From d330e67f8bad1d5c517970cd50b01f8ee8c3d7e3 Mon Sep 17 00:00:00 2001
From: Martin Robinson <martinjrobins@gmail.com>
Date: Fri, 6 Mar 2020 18:08:44 +0000
Subject: [PATCH 09/30] #858 added make_semi_explicit to convert implicit dae
 equations to algebraic

---
 pybamm/__init__.py                            |  1 +
 pybamm/discretisations/discretisation.py      | 12 +++++--
 .../operations/convert_to_casadi.py           | 23 +++++++-----
 pybamm/expression_tree/symbol.py              |  6 ++--
 pybamm/expression_tree/variable.py            | 31 ++++++++++++++--
 pybamm/models/base_model.py                   |  9 +++++
 pybamm/processed_variable.py                  | 22 ++++++------
 pybamm/solvers/base_solver.py                 | 11 +++---
 pybamm/solvers/casadi_solver.py               |  8 +++--
 .../test_discretisation.py                    | 19 +++++++++-
 tests/unit/test_solvers/test_casadi_solver.py | 35 +++++++++++++++++++
 11 files changed, 140 insertions(+), 37 deletions(-)

diff --git a/pybamm/__init__.py b/pybamm/__init__.py
index 47ab5bfe21..f2d3d14bc6 100644
--- a/pybamm/__init__.py
+++ b/pybamm/__init__.py
@@ -153,6 +153,7 @@ def version(formatted=False):
 from .models import standard_variables
 from .models.event import Event
 from .models.event import EventType
+from .models.make_semi_explicit import make_semi_explicit
 
 # Battery models
 from .models.full_battery_models.base_battery_model import BaseBatteryModel
diff --git a/pybamm/discretisations/discretisation.py b/pybamm/discretisations/discretisation.py
index efbfe46772..01fc1d6e0b 100644
--- a/pybamm/discretisations/discretisation.py
+++ b/pybamm/discretisations/discretisation.py
@@ -997,9 +997,9 @@ def check_model(self, model):
         self.check_initial_conditions(model)
         self.check_initial_conditions_rhs(model)
         self.check_variables(model)
-        self.check_for_time_derivatives_in_rhs(model)
+        self.check_for_time_derivatives(model)
 
-    def check_for_time_derivatives_in_rhs(self, model):
+    def check_for_time_derivatives(self, model):
         # Check that no variable time derivatives exist in the rhs equations
         for eq in model.rhs.values():
             for node in eq.pre_order():
@@ -1008,6 +1008,14 @@ def check_for_time_derivatives_in_rhs(self, model):
                         "time derivative of variable found ({}) in rhs equation"
                     )
 
+        # Check that no variable time derivatives exist in the algebraic equations
+        for eq in model.algebraic.values():
+            for node in eq.pre_order():
+                if isinstance(node, pybamm.VariableDot):
+                    raise pybamm.ModelError(
+                        "time derivative of variable found ({}) in algebraic equation"
+                    )
+
     def check_initial_conditions(self, model):
         """Check initial conditions are a numpy array"""
         # Individual
diff --git a/pybamm/expression_tree/operations/convert_to_casadi.py b/pybamm/expression_tree/operations/convert_to_casadi.py
index 5e89ab75cf..cb82217129 100644
--- a/pybamm/expression_tree/operations/convert_to_casadi.py
+++ b/pybamm/expression_tree/operations/convert_to_casadi.py
@@ -13,7 +13,7 @@ def __init__(self, casadi_symbols=None):
 
         pybamm.citations.register("Andersson2019")
 
-    def convert(self, symbol, t=None, y=None, u=None):
+    def convert(self, symbol, t=None, y=None, y_dot=None, u=None):
         """
         This function recurses down the tree, converting the PyBaMM expression tree to
         a CasADi expression tree
@@ -39,12 +39,12 @@ def convert(self, symbol, t=None, y=None, u=None):
         except KeyError:
             # Change u to empty dictionary if it's None
             u = u or {}
-            casadi_symbol = self._convert(symbol, t, y, u)
+            casadi_symbol = self._convert(symbol, t, y, y_dot, u)
             self._casadi_symbols[symbol.id] = casadi_symbol
 
             return casadi_symbol
 
-    def _convert(self, symbol, t=None, y=None, u=None):
+    def _convert(self, symbol, t=None, y=None, y_dot=None, u=None):
         """ See :meth:`CasadiConverter.convert()`. """
         if isinstance(
             symbol,
@@ -56,30 +56,35 @@ def _convert(self, symbol, t=None, y=None, u=None):
                 pybamm.ExternalVariable,
             ),
         ):
-            return casadi.MX(symbol.evaluate(t, y, u))
+            return casadi.MX(symbol.evaluate(t, y, y_dot, u))
 
         elif isinstance(symbol, pybamm.StateVector):
             if y is None:
                 raise ValueError("Must provide a 'y' for converting state vectors")
             return casadi.vertcat(*[y[y_slice] for y_slice in symbol.y_slices])
 
+        elif isinstance(symbol, pybamm.StateVectorDot):
+            if y_dot is None:
+                raise ValueError("Must provide a 'y_dot' for converting state vectors")
+            return casadi.vertcat(*[y_dot[y_slice] for y_slice in symbol.y_slices])
+
         elif isinstance(symbol, pybamm.BinaryOperator):
             left, right = symbol.children
             # process children
-            converted_left = self.convert(left, t, y, u)
-            converted_right = self.convert(right, t, y, u)
+            converted_left = self.convert(left, t, y, y_dot, u)
+            converted_right = self.convert(right, t, y, y_dot, u)
             # _binary_evaluate defined in derived classes for specific rules
             return symbol._binary_evaluate(converted_left, converted_right)
 
         elif isinstance(symbol, pybamm.UnaryOperator):
-            converted_child = self.convert(symbol.child, t, y, u)
+            converted_child = self.convert(symbol.child, t, y, y_dot, u)
             if isinstance(symbol, pybamm.AbsoluteValue):
                 return casadi.fabs(converted_child)
             return symbol._unary_evaluate(converted_child)
 
         elif isinstance(symbol, pybamm.Function):
             converted_children = [
-                self.convert(child, t, y, u) for child in symbol.children
+                self.convert(child, t, y, y_dot, u) for child in symbol.children
             ]
             # Special functions
             if symbol.function == np.min:
@@ -110,7 +115,7 @@ def _convert(self, symbol, t=None, y=None, u=None):
                 return symbol._function_evaluate(converted_children)
         elif isinstance(symbol, pybamm.Concatenation):
             converted_children = [
-                self.convert(child, t, y, u) for child in symbol.children
+                self.convert(child, t, y, y_dot, u) for child in symbol.children
             ]
             if isinstance(symbol, (pybamm.NumpyConcatenation, pybamm.SparseStack)):
                 return casadi.vertcat(*converted_children)
diff --git a/pybamm/expression_tree/symbol.py b/pybamm/expression_tree/symbol.py
index f09008d1c0..912b102b6d 100644
--- a/pybamm/expression_tree/symbol.py
+++ b/pybamm/expression_tree/symbol.py
@@ -677,12 +677,12 @@ def simplify(self, simplified_symbols=None):
         """ Simplify the expression tree. See :class:`pybamm.Simplification`. """
         return pybamm.Simplification(simplified_symbols).simplify(self)
 
-    def to_casadi(self, t=None, y=None, u=None, casadi_symbols=None):
+    def to_casadi(self, t=None, y=None, ydot=None, u=None, casadi_symbols=None):
         """
         Convert the expression tree to a CasADi expression tree.
         See :class:`pybamm.CasadiConverter`.
         """
-        return pybamm.CasadiConverter(casadi_symbols).convert(self, t, y, u)
+        return pybamm.CasadiConverter(casadi_symbols).convert(self, t, y, ydot, u)
 
     def new_copy(self):
         """
@@ -712,7 +712,7 @@ def shape(self):
         # Try with some large y, to avoid having to use pre_order (slow)
         try:
             y = np.linspace(0.1, 0.9, int(1e4))
-            evaluated_self = self.evaluate(0, y, u="shape test")
+            evaluated_self = self.evaluate(0, y, y, u="shape test")
         # If that fails, fall back to calculating how big y should really be
         except ValueError:
             state_vectors_in_node = [
diff --git a/pybamm/expression_tree/variable.py b/pybamm/expression_tree/variable.py
index 7a9b9cd618..8650c3e0c3 100644
--- a/pybamm/expression_tree/variable.py
+++ b/pybamm/expression_tree/variable.py
@@ -6,7 +6,7 @@
 import numpy as np
 
 
-class Variable(pybamm.Symbol):
+class VariableBase(pybamm.Symbol):
     """A node in the expression tree represending a dependent variable
 
     This node will be discretised by :class:`.Discretisation` and converted
@@ -47,6 +47,32 @@ def _evaluate_for_shape(self):
             self.domain, self.auxiliary_domains
         )
 
+class Variable(VariableBase):
+    """A node in the expression tree represending a dependent variable
+
+    This node will be discretised by :class:`.Discretisation` and converted
+    to a :class:`pybamm.StateVector` node.
+
+    Parameters
+    ----------
+
+    name : str
+        name of the node
+    domain : iterable of str
+        list of domains that this variable is valid over
+    auxiliary_domains : dict
+        dictionary of auxiliary domains ({'secondary': ..., 'tertiary': ...}). For
+        example, for the single particle model, the particle concentration would be a
+        Variable with domain 'negative particle' and secondary auxiliary domain 'current
+        collector'. For the DFN, the particle concentration would be a Variable with
+        domain 'negative particle', secondary domain 'negative electrode' and tertiary
+        domain 'current collector'
+
+    *Extends:* :class:`Symbol`
+    """
+    def __init__(self, name, domain=None, auxiliary_domains=None):
+        super().__init__(name, domain=domain, auxiliary_domains=auxiliary_domains)
+
     def _jac(self, variable):
         if variable.id == self.id:
             return pybamm.Scalar(1)
@@ -57,8 +83,7 @@ def _jac(self, variable):
         else:
             return pybamm.Scalar(0)
 
-
-class VariableDot(Variable):
+class VariableDot(VariableBase):
     """
     A node in the expression tree represending the time derviative of a dependent
     variable
diff --git a/pybamm/models/base_model.py b/pybamm/models/base_model.py
index f8ae8c5c17..158a896fc4 100644
--- a/pybamm/models/base_model.py
+++ b/pybamm/models/base_model.py
@@ -410,6 +410,9 @@ def check_well_determined(self, post_discretisation):
             vars_in_eqns.update(
                 [x.id for x in eqn.pre_order() if isinstance(x, pybamm.Variable)]
             )
+            vars_in_eqns.update(
+                [x.get_variable().id for x in eqn.pre_order() if isinstance(x, pybamm.VariableDot)]
+            )
         for var, eqn in self.algebraic.items():
             vars_in_algebraic_keys.update(
                 [x.id for x in var.pre_order() if isinstance(x, pybamm.Variable)]
@@ -417,11 +420,17 @@ def check_well_determined(self, post_discretisation):
             vars_in_eqns.update(
                 [x.id for x in eqn.pre_order() if isinstance(x, pybamm.Variable)]
             )
+            vars_in_eqns.update(
+                [x.get_variable().id for x in eqn.pre_order() if isinstance(x, pybamm.VariableDot)]
+            )
         for var, side_eqn in self.boundary_conditions.items():
             for side, (eqn, typ) in side_eqn.items():
                 vars_in_eqns.update(
                     [x.id for x in eqn.pre_order() if isinstance(x, pybamm.Variable)]
                 )
+                vars_in_eqns.update(
+                    [x.get_variable().id for x in eqn.pre_order() if isinstance(x, pybamm.VariableDot)]
+                )
         # If any keys are repeated between rhs and algebraic then the model is
         # overdetermined
         if not set(vars_in_rhs_keys).isdisjoint(vars_in_algebraic_keys):
diff --git a/pybamm/processed_variable.py b/pybamm/processed_variable.py
index 23d74a461c..19a3df714a 100644
--- a/pybamm/processed_variable.py
+++ b/pybamm/processed_variable.py
@@ -41,14 +41,14 @@ def __init__(self, base_variable, solution, known_evals=None):
             self.base_eval, self.known_evals[solution.t[0]] = base_variable.evaluate(
                 solution.t[0],
                 solution.y[:, 0],
-                {name: inp[0] for name, inp in solution.inputs.items()},
+                u={name: inp[0] for name, inp in solution.inputs.items()},
                 known_evals=self.known_evals[solution.t[0]],
             )
         else:
             self.base_eval = base_variable.evaluate(
                 solution.t[0],
                 solution.y[:, 0],
-                {name: inp[0] for name, inp in solution.inputs.items()},
+                u={name: inp[0] for name, inp in solution.inputs.items()},
             )
 
         # handle 2D (in space) finite element variables differently
@@ -96,10 +96,10 @@ def initialise_1D(self):
             inputs = {name: inp[idx] for name, inp in self.inputs.items()}
             if self.known_evals:
                 entries[idx], self.known_evals[t] = self.base_variable.evaluate(
-                    t, u, inputs, known_evals=self.known_evals[t]
+                    t, u, u=inputs, known_evals=self.known_evals[t]
                 )
             else:
-                entries[idx] = self.base_variable.evaluate(t, u, inputs)
+                entries[idx] = self.base_variable.evaluate(t, u, u=inputs)
 
         # No discretisation provided, or variable has no domain (function of t only)
         self._interpolation_function = interp.interp1d(
@@ -120,12 +120,12 @@ def initialise_2D(self):
             inputs = {name: inp[idx] for name, inp in self.inputs.items()}
             if self.known_evals:
                 eval_and_known_evals = self.base_variable.evaluate(
-                    t, u, inputs, known_evals=self.known_evals[t]
+                    t, u, u=inputs, known_evals=self.known_evals[t]
                 )
                 entries[:, idx] = eval_and_known_evals[0][:, 0]
                 self.known_evals[t] = eval_and_known_evals[1]
             else:
-                entries[:, idx] = self.base_variable.evaluate(t, u, inputs)[:, 0]
+                entries[:, idx] = self.base_variable.evaluate(t, u, u=inputs)[:, 0]
 
         # Process the discretisation to get x values
         nodes = self.mesh[0].nodes
@@ -226,7 +226,7 @@ def initialise_3D(self):
             inputs = {name: inp[idx] for name, inp in self.inputs.items()}
             if self.known_evals:
                 eval_and_known_evals = self.base_variable.evaluate(
-                    t, u, inputs, known_evals=self.known_evals[t]
+                    t, u, u=inputs, known_evals=self.known_evals[t]
                 )
                 entries[:, :, idx] = np.reshape(
                     eval_and_known_evals[0],
@@ -236,7 +236,7 @@ def initialise_3D(self):
                 self.known_evals[t] = eval_and_known_evals[1]
             else:
                 entries[:, :, idx] = np.reshape(
-                    self.base_variable.evaluate(t, u, inputs),
+                    self.base_variable.evaluate(t, u, u=inputs),
                     [first_dim_size, second_dim_size],
                     order="F",
                 )
@@ -265,7 +265,7 @@ def initialise_2Dspace_scikit_fem(self):
         inputs = {name: inp[0] for name, inp in self.inputs.items()}
 
         entries = np.reshape(
-            self.base_variable.evaluate(0, self.u_sol, inputs), [len_y, len_z]
+            self.base_variable.evaluate(0, self.u_sol, u=inputs), [len_y, len_z]
         )
 
         # assign attributes for reference
@@ -296,13 +296,13 @@ def initialise_3D_scikit_fem(self):
 
             if self.known_evals:
                 eval_and_known_evals = self.base_variable.evaluate(
-                    t, u, inputs, known_evals=self.known_evals[t]
+                    t, u, u=inputs, known_evals=self.known_evals[t]
                 )
                 entries[:, :, idx] = np.reshape(eval_and_known_evals[0], [len_y, len_z])
                 self.known_evals[t] = eval_and_known_evals[1]
             else:
                 entries[:, :, idx] = np.reshape(
-                    self.base_variable.evaluate(t, u, inputs), [len_y, len_z]
+                    self.base_variable.evaluate(t, u, u=inputs), [len_y, len_z]
                 )
 
         # assign attributes for reference
diff --git a/pybamm/solvers/base_solver.py b/pybamm/solvers/base_solver.py
index 648edf6469..a29172f147 100644
--- a/pybamm/solvers/base_solver.py
+++ b/pybamm/solvers/base_solver.py
@@ -148,10 +148,10 @@ def set_up(self, model, inputs=None):
             # Convert model attributes to casadi
             t_casadi = casadi.MX.sym("t")
             y_diff = casadi.MX.sym(
-                "y_diff", len(model.concatenated_rhs.evaluate(0, y0, inputs))
+                "y_diff", len(model.concatenated_rhs.evaluate(0, y0, u=inputs))
             )
             y_alg = casadi.MX.sym(
-                "y_alg", len(model.concatenated_algebraic.evaluate(0, y0, inputs))
+                "y_alg", len(model.concatenated_algebraic.evaluate(0, y0, u=inputs))
             )
             y_casadi = casadi.vertcat(y_diff, y_alg)
             u_casadi = {}
@@ -292,7 +292,7 @@ def report(string):
             model.residuals_eval = residuals_eval
             model.jacobian_eval = jacobian_eval
             y0_guess = y0.flatten()
-            model.y0 = self.calculate_consistent_state(model, 0, y0_guess, inputs)
+            model.y0 = self.calculate_consistent_state(model, 0, y0_guess,inputs)
         else:
             # can use DAE solver to solve ODE model
             model.residuals_eval = Residuals(rhs, "residuals", model)
@@ -812,11 +812,12 @@ def __call__(self, t, y):
 
     def function(self, t, y):
         if self.form == "casadi":
+            states_eval = self._function(t, y, self.inputs)
             if self.name in ["RHS", "algebraic", "residuals", "event"]:
-                return self._function(t, y, self.inputs).full()
+                return states_eval.full()
             else:
                 # keep jacobians sparse
-                return self._function(t, y, self.inputs)
+                return states_eval
         else:
             return self._function(t, y, self.inputs, known_evals={})[0]
 
diff --git a/pybamm/solvers/casadi_solver.py b/pybamm/solvers/casadi_solver.py
index 1029a8560e..c1415007fb 100644
--- a/pybamm/solvers/casadi_solver.py
+++ b/pybamm/solvers/casadi_solver.py
@@ -173,6 +173,7 @@ def get_integrator(self, model, t_eval, inputs):
         # Only set up problem once
         if model not in self.problems:
             y0 = model.y0
+            ydot0 = model.ydot0
             rhs = model.casadi_rhs
             algebraic = model.casadi_algebraic
             u_stacked = casadi.vertcat(*[x for x in inputs.values()])
@@ -190,21 +191,22 @@ def get_integrator(self, model, t_eval, inputs):
             u = casadi.MX.sym("u", u_stacked.shape[0])
             y_diff = casadi.MX.sym("y_diff", rhs(t_eval[0], y0, u).shape[0])
             problem = {"t": t, "x": y_diff, "p": u}
-            if algebraic(t_eval[0], y0, u).is_empty():
+            if algebraic(t_eval[0], y0, ydot0, u).is_empty():
                 method = "cvodes"
                 problem.update({"ode": rhs(t, y_diff, u)})
             else:
                 options["calc_ic"] = True
                 method = "idas"
-                y_alg = casadi.MX.sym("y_alg", algebraic(t_eval[0], y0, u).shape[0])
+                y_alg = casadi.MX.sym("y_alg", algebraic(t_eval[0], y0, ydot0, u).shape[0])
                 y_full = casadi.vertcat(y_diff, y_alg)
                 problem.update(
                     {
                         "z": y_alg,
                         "ode": rhs(t, y_full, u),
-                        "alg": algebraic(t, y_full, u),
+                        "alg": algebraic(t, y_full, ydot0, u),
                     }
                 )
+            print('using method' ,method)
             self.problems[model] = problem
             self.options[model] = options
             self.methods[model] = method
diff --git a/tests/unit/test_discretisations/test_discretisation.py b/tests/unit/test_discretisations/test_discretisation.py
index 28df731ec9..0e9bbe9ba2 100644
--- a/tests/unit/test_discretisations/test_discretisation.py
+++ b/tests/unit/test_discretisations/test_discretisation.py
@@ -701,7 +701,7 @@ def test_process_model_ode(self):
         with self.assertRaises(pybamm.ModelError):
             disc.process_model(model)
 
-        # test that ill possed model with time derivatives of variables in rhs raises an
+        # test that any time derivatives of variables in rhs raises an
         # error
         model = pybamm.BaseModel()
         model.rhs = {c: pybamm.div(N) + pybamm.d_dt(c), T: pybamm.div(q), S: pybamm.div(p)}
@@ -719,6 +719,7 @@ def test_process_model_ode(self):
         with self.assertRaises(pybamm.ModelError):
             disc.process_model(model)
 
+
     def test_process_model_dae(self):
         # one rhs equation and one algebraic
         whole_cell = ["negative electrode", "separator", "positive electrode"]
@@ -816,6 +817,22 @@ def test_process_model_dae(self):
         jacobian = expr.evaluate(0, y0, known_evals=known_evals)[0]
         np.testing.assert_array_equal(jacobian_actual, jacobian.toarray())
 
+        # check that any time derivatives of variables in algebraic raises an
+        # error
+        model = pybamm.BaseModel()
+        model.rhs = {c: pybamm.div(N)}
+        model.algebraic = {d: d - 2 * pybamm.d_dt(c)}
+        model.initial_conditions = {d: pybamm.Scalar(6), c: pybamm.Scalar(3)}
+        model.boundary_conditions = {
+            c: {"left": (0, "Neumann"), "right": (0, "Neumann")}
+        }
+        model.variables = {"c": c, "N": N, "d": d}
+
+        with self.assertRaises(pybamm.ModelError):
+            disc.process_model(model)
+
+
+
     def test_process_model_concatenation(self):
         # concatenation of variables as the key
         cn = pybamm.Variable("c", domain=["negative electrode"])
diff --git a/tests/unit/test_solvers/test_casadi_solver.py b/tests/unit/test_solvers/test_casadi_solver.py
index 7319768462..27094449c4 100644
--- a/tests/unit/test_solvers/test_casadi_solver.py
+++ b/tests/unit/test_solvers/test_casadi_solver.py
@@ -309,6 +309,40 @@ def test_model_solver_with_non_identity_mass(self):
         np.testing.assert_allclose(solution.y[0], np.exp(0.1 * solution.t))
         np.testing.assert_allclose(solution.y[-1], 2 * np.exp(0.1 * solution.t))
 
+    def test_model_solver_with_dvdt(self):
+        model = pybamm.BaseModel()
+        var1 = pybamm.Variable("var1", domain="negative electrode")
+        var2 = pybamm.Variable("var2", domain="negative electrode")
+        model.rhs = {var1: -2 * var1 * pybamm.t}
+        model.algebraic = {var2: var2 - pybamm.d_dt(var1)}
+        model.initial_conditions = {var1: 1, var2: 0}
+        print('before semi explicit')
+        for key,value in model.rhs.items():
+            print('{}: {}'.format(key, value))
+        for key,value in model.algebraic.items():
+            print('{}: {}'.format(key, value))
+        pybamm.make_semi_explicit(model)
+        print('after semi explicit')
+        for key,value in model.rhs.items():
+            print('{}: {}'.format(key, value))
+        for key,value in model.algebraic.items():
+            print('{}: {}'.format(key, value))
+        disc = get_discretisation_for_testing()
+        disc.process_model(model)
+
+        # Solve
+        solver = pybamm.CasadiSolver(rtol=1e-8, atol=1e-8)
+        t_eval = np.linspace(0, 1, 100)
+        solution = solver.solve(model, t_eval)
+        np.testing.assert_array_equal(solution.t, t_eval)
+        import matplotlib.pyplot as plt
+        plt.plot(solution.y[0])
+        plt.plot(np.exp(-solution.t**2))
+        plt.show()
+        np.testing.assert_allclose(solution.y[0], np.exp(-solution.t ** 2), rtol=1e-06)
+        np.testing.assert_allclose(solution.y[-1],
+                                -2 * solution.t * np.exp(-solution.t**2))
+
 
 if __name__ == "__main__":
     print("Add -v for more debug output")
@@ -317,4 +351,5 @@ def test_model_solver_with_non_identity_mass(self):
     if "-v" in sys.argv:
         debug = True
     pybamm.settings.debug_mode = True
+    #pybamm.set_logging_level("DEBUG")
     unittest.main()

From eb1b3810e4f6367d1e5a7340f54bf46b7ccd47e2 Mon Sep 17 00:00:00 2001
From: Martin Robinson <martinjrobins@gmail.com>
Date: Fri, 6 Mar 2020 18:28:53 +0000
Subject: [PATCH 10/30] #858 make_semi_explicit updates initial conditions

---
 pybamm/solvers/casadi_solver.py               |  8 +++-----
 tests/unit/test_solvers/test_casadi_solver.py | 18 ++----------------
 2 files changed, 5 insertions(+), 21 deletions(-)

diff --git a/pybamm/solvers/casadi_solver.py b/pybamm/solvers/casadi_solver.py
index c1415007fb..1029a8560e 100644
--- a/pybamm/solvers/casadi_solver.py
+++ b/pybamm/solvers/casadi_solver.py
@@ -173,7 +173,6 @@ def get_integrator(self, model, t_eval, inputs):
         # Only set up problem once
         if model not in self.problems:
             y0 = model.y0
-            ydot0 = model.ydot0
             rhs = model.casadi_rhs
             algebraic = model.casadi_algebraic
             u_stacked = casadi.vertcat(*[x for x in inputs.values()])
@@ -191,22 +190,21 @@ def get_integrator(self, model, t_eval, inputs):
             u = casadi.MX.sym("u", u_stacked.shape[0])
             y_diff = casadi.MX.sym("y_diff", rhs(t_eval[0], y0, u).shape[0])
             problem = {"t": t, "x": y_diff, "p": u}
-            if algebraic(t_eval[0], y0, ydot0, u).is_empty():
+            if algebraic(t_eval[0], y0, u).is_empty():
                 method = "cvodes"
                 problem.update({"ode": rhs(t, y_diff, u)})
             else:
                 options["calc_ic"] = True
                 method = "idas"
-                y_alg = casadi.MX.sym("y_alg", algebraic(t_eval[0], y0, ydot0, u).shape[0])
+                y_alg = casadi.MX.sym("y_alg", algebraic(t_eval[0], y0, u).shape[0])
                 y_full = casadi.vertcat(y_diff, y_alg)
                 problem.update(
                     {
                         "z": y_alg,
                         "ode": rhs(t, y_full, u),
-                        "alg": algebraic(t, y_full, ydot0, u),
+                        "alg": algebraic(t, y_full, u),
                     }
                 )
-            print('using method' ,method)
             self.problems[model] = problem
             self.options[model] = options
             self.methods[model] = method
diff --git a/tests/unit/test_solvers/test_casadi_solver.py b/tests/unit/test_solvers/test_casadi_solver.py
index 27094449c4..39781ed088 100644
--- a/tests/unit/test_solvers/test_casadi_solver.py
+++ b/tests/unit/test_solvers/test_casadi_solver.py
@@ -316,17 +316,7 @@ def test_model_solver_with_dvdt(self):
         model.rhs = {var1: -2 * var1 * pybamm.t}
         model.algebraic = {var2: var2 - pybamm.d_dt(var1)}
         model.initial_conditions = {var1: 1, var2: 0}
-        print('before semi explicit')
-        for key,value in model.rhs.items():
-            print('{}: {}'.format(key, value))
-        for key,value in model.algebraic.items():
-            print('{}: {}'.format(key, value))
         pybamm.make_semi_explicit(model)
-        print('after semi explicit')
-        for key,value in model.rhs.items():
-            print('{}: {}'.format(key, value))
-        for key,value in model.algebraic.items():
-            print('{}: {}'.format(key, value))
         disc = get_discretisation_for_testing()
         disc.process_model(model)
 
@@ -335,13 +325,9 @@ def test_model_solver_with_dvdt(self):
         t_eval = np.linspace(0, 1, 100)
         solution = solver.solve(model, t_eval)
         np.testing.assert_array_equal(solution.t, t_eval)
-        import matplotlib.pyplot as plt
-        plt.plot(solution.y[0])
-        plt.plot(np.exp(-solution.t**2))
-        plt.show()
-        np.testing.assert_allclose(solution.y[0], np.exp(-solution.t ** 2), rtol=1e-06)
+        np.testing.assert_allclose(solution.y[0], np.exp(-solution.t**2), rtol=1e-06)
         np.testing.assert_allclose(solution.y[-1],
-                                -2 * solution.t * np.exp(-solution.t**2))
+                                -2 * solution.t * np.exp(-solution.t**2), rtol=1e-06)
 
 
 if __name__ == "__main__":

From 5a0da61c883ff80b482aada8972dc0b02325dcd8 Mon Sep 17 00:00:00 2001
From: Martin Robinson <martinjrobins@gmail.com>
Date: Sat, 7 Mar 2020 05:32:22 +0000
Subject: [PATCH 11/30] #858 some fixes, move time deriv model checks to
 base_model

---
 pybamm/discretisations/discretisation.py      | 18 ---------
 .../operations/convert_to_casadi.py           |  4 +-
 pybamm/models/base_model.py                   | 39 +++++++++++++++++--
 pybamm/solvers/base_solver.py                 |  7 ++--
 tests/unit/test_solvers/test_casadi_solver.py |  2 +-
 5 files changed, 43 insertions(+), 27 deletions(-)

diff --git a/pybamm/discretisations/discretisation.py b/pybamm/discretisations/discretisation.py
index 01fc1d6e0b..847f5b662f 100644
--- a/pybamm/discretisations/discretisation.py
+++ b/pybamm/discretisations/discretisation.py
@@ -997,24 +997,6 @@ def check_model(self, model):
         self.check_initial_conditions(model)
         self.check_initial_conditions_rhs(model)
         self.check_variables(model)
-        self.check_for_time_derivatives(model)
-
-    def check_for_time_derivatives(self, model):
-        # Check that no variable time derivatives exist in the rhs equations
-        for eq in model.rhs.values():
-            for node in eq.pre_order():
-                if isinstance(node, pybamm.VariableDot):
-                    raise pybamm.ModelError(
-                        "time derivative of variable found ({}) in rhs equation"
-                    )
-
-        # Check that no variable time derivatives exist in the algebraic equations
-        for eq in model.algebraic.values():
-            for node in eq.pre_order():
-                if isinstance(node, pybamm.VariableDot):
-                    raise pybamm.ModelError(
-                        "time derivative of variable found ({}) in algebraic equation"
-                    )
 
     def check_initial_conditions(self, model):
         """Check initial conditions are a numpy array"""
diff --git a/pybamm/expression_tree/operations/convert_to_casadi.py b/pybamm/expression_tree/operations/convert_to_casadi.py
index cb82217129..70ed35b694 100644
--- a/pybamm/expression_tree/operations/convert_to_casadi.py
+++ b/pybamm/expression_tree/operations/convert_to_casadi.py
@@ -13,7 +13,7 @@ def __init__(self, casadi_symbols=None):
 
         pybamm.citations.register("Andersson2019")
 
-    def convert(self, symbol, t=None, y=None, y_dot=None, u=None):
+    def convert(self, symbol, t, y, y_dot, u):
         """
         This function recurses down the tree, converting the PyBaMM expression tree to
         a CasADi expression tree
@@ -44,7 +44,7 @@ def convert(self, symbol, t=None, y=None, y_dot=None, u=None):
 
             return casadi_symbol
 
-    def _convert(self, symbol, t=None, y=None, y_dot=None, u=None):
+    def _convert(self, symbol, t, y, y_dot, u):
         """ See :meth:`CasadiConverter.convert()`. """
         if isinstance(
             symbol,
diff --git a/pybamm/models/base_model.py b/pybamm/models/base_model.py
index 158a896fc4..6f90a0d266 100644
--- a/pybamm/models/base_model.py
+++ b/pybamm/models/base_model.py
@@ -385,12 +385,42 @@ def check_well_posedness(self, post_discretisation=False):
         self.check_algebraic_equations(post_discretisation)
         self.check_ics_bcs()
         self.check_default_variables_dictionaries()
+        self.check_for_time_derivatives()
         # Can't check variables after discretising, since Variable objects get replaced
         # by StateVector objects
         # Checking variables is slow, so only do it in debug mode
         if pybamm.settings.debug_mode is True and post_discretisation is False:
             self.check_variables()
 
+    def check_for_time_derivatives(self):
+        # Check that no variable time derivatives exist in the rhs equations
+        for key, eq in self.rhs.items():
+            for node in eq.pre_order():
+                if isinstance(node, pybamm.VariableDot):
+                    raise pybamm.ModelError(
+                        "time derivative of variable found ({}) in rhs equation {}"
+                        .format(node, key)
+                    )
+                if isinstance(node, pybamm.StateVectorDot):
+                    raise pybamm.ModelError(
+                        "time derivative of state vector found ({}) in rhs equation {}"
+                        .format(node, key)
+                    )
+
+        # Check that no variable time derivatives exist in the algebraic equations
+        for key, eq in self.algebraic.items():
+            for node in eq.pre_order():
+                if isinstance(node, pybamm.VariableDot):
+                    raise pybamm.ModelError(
+                        "time derivative of variable found ({}) in algebraic equation {}"
+                        .format(node, key)
+                    )
+                if isinstance(node, pybamm.StateVectorDot):
+                    raise pybamm.ModelError(
+                        "time derivative of state vector found ({}) in algebraic"
+                        "equation {}".format(node, key)
+                    )
+
     def check_well_determined(self, post_discretisation):
         """ Check that the model is not under- or over-determined. """
         # Equations (differential and algebraic)
@@ -411,7 +441,8 @@ def check_well_determined(self, post_discretisation):
                 [x.id for x in eqn.pre_order() if isinstance(x, pybamm.Variable)]
             )
             vars_in_eqns.update(
-                [x.get_variable().id for x in eqn.pre_order() if isinstance(x, pybamm.VariableDot)]
+                [x.get_variable().id for x in eqn.pre_order()
+                 if isinstance(x, pybamm.VariableDot)]
             )
         for var, eqn in self.algebraic.items():
             vars_in_algebraic_keys.update(
@@ -421,7 +452,8 @@ def check_well_determined(self, post_discretisation):
                 [x.id for x in eqn.pre_order() if isinstance(x, pybamm.Variable)]
             )
             vars_in_eqns.update(
-                [x.get_variable().id for x in eqn.pre_order() if isinstance(x, pybamm.VariableDot)]
+                [x.get_variable().id for x in eqn.pre_order()
+                 if isinstance(x, pybamm.VariableDot)]
             )
         for var, side_eqn in self.boundary_conditions.items():
             for side, (eqn, typ) in side_eqn.items():
@@ -429,7 +461,8 @@ def check_well_determined(self, post_discretisation):
                     [x.id for x in eqn.pre_order() if isinstance(x, pybamm.Variable)]
                 )
                 vars_in_eqns.update(
-                    [x.get_variable().id for x in eqn.pre_order() if isinstance(x, pybamm.VariableDot)]
+                    [x.get_variable().id for x in eqn.pre_order()
+                     if isinstance(x, pybamm.VariableDot)]
                 )
         # If any keys are repeated between rhs and algebraic then the model is
         # overdetermined
diff --git a/pybamm/solvers/base_solver.py b/pybamm/solvers/base_solver.py
index a29172f147..33f94d65d0 100644
--- a/pybamm/solvers/base_solver.py
+++ b/pybamm/solvers/base_solver.py
@@ -117,7 +117,7 @@ def set_up(self, model, inputs=None):
 
         """
         inputs = inputs or {}
-        y0 = model.concatenated_initial_conditions.evaluate(0, None, inputs)
+        y0 = model.concatenated_initial_conditions.evaluate(0, None, u=inputs)
 
         # Set model timescale
         model.timescale_eval = model.timescale.evaluate(u=inputs)
@@ -194,7 +194,8 @@ def report(string):
             else:
                 # Process with CasADi
                 report(f"Converting {name} to CasADi")
-                func = func.to_casadi(t_casadi, y_casadi, u_casadi)
+                print(u_casadi)
+                func = func.to_casadi(t_casadi, y_casadi, u=u_casadi)
                 if use_jacobian:
                     report(f"Calculating jacobian for {name} using CasADi")
                     jac_casadi = casadi.jacobian(func, y_casadi)
@@ -767,7 +768,7 @@ def get_termination_reason(self, solution, events):
                         event.expression.evaluate(
                             solution.t_event,
                             solution.y_event,
-                            {k: v[-1] for k, v in solution.inputs.items()},
+                            u={k: v[-1] for k, v in solution.inputs.items()},
                         )
                     )
             termination_event = min(final_event_values, key=final_event_values.get)
diff --git a/tests/unit/test_solvers/test_casadi_solver.py b/tests/unit/test_solvers/test_casadi_solver.py
index 39781ed088..e63e88ce59 100644
--- a/tests/unit/test_solvers/test_casadi_solver.py
+++ b/tests/unit/test_solvers/test_casadi_solver.py
@@ -337,5 +337,5 @@ def test_model_solver_with_dvdt(self):
     if "-v" in sys.argv:
         debug = True
     pybamm.settings.debug_mode = True
-    #pybamm.set_logging_level("DEBUG")
+    pybamm.set_logging_level("DEBUG")
     unittest.main()

From 9ba63f684bd2a8c09e904cf51e071885ea848cb0 Mon Sep 17 00:00:00 2001
From: Martin Robinson <martinjrobins@gmail.com>
Date: Sat, 7 Mar 2020 06:03:35 +0000
Subject: [PATCH 12/30] #858 misc fixes

---
 pybamm/discretisations/discretisation.py          | 4 ++++
 pybamm/expression_tree/symbol.py                  | 3 ++-
 pybamm/solvers/base_solver.py                     | 1 -
 tests/unit/test_parameters/test_parameters_cli.py | 8 ++++++++
 tests/unit/test_solvers/test_casadi_solver.py     | 1 -
 5 files changed, 14 insertions(+), 3 deletions(-)

diff --git a/pybamm/discretisations/discretisation.py b/pybamm/discretisations/discretisation.py
index 847f5b662f..bcfe9949da 100644
--- a/pybamm/discretisations/discretisation.py
+++ b/pybamm/discretisations/discretisation.py
@@ -1030,6 +1030,10 @@ def check_initial_conditions_rhs(self, model):
         y0 = model.concatenated_initial_conditions
         # Individual
         for var in model.rhs.keys():
+            print('rhs')
+            print(model.rhs[var])
+            print('init')
+            print(model.initial_conditions[var])
             assert (
                 model.rhs[var].shape == model.initial_conditions[var].shape
             ), pybamm.ModelError(
diff --git a/pybamm/expression_tree/symbol.py b/pybamm/expression_tree/symbol.py
index 912b102b6d..06b9693557 100644
--- a/pybamm/expression_tree/symbol.py
+++ b/pybamm/expression_tree/symbol.py
@@ -713,6 +713,7 @@ def shape(self):
         try:
             y = np.linspace(0.1, 0.9, int(1e4))
             evaluated_self = self.evaluate(0, y, y, u="shape test")
+            print('evaluated self is ',evaluated_self)
         # If that fails, fall back to calculating how big y should really be
         except ValueError:
             state_vectors_in_node = [
@@ -726,7 +727,7 @@ def shape(self):
                 )
                 # Pick a y that won't cause RuntimeWarnings
                 y = np.linspace(0.1, 0.9, min_y_size)
-            evaluated_self = self.evaluate(0, y)
+            evaluated_self = self.evaluate(0, y, y)
 
         # Return shape of evaluated object
         if isinstance(evaluated_self, numbers.Number):
diff --git a/pybamm/solvers/base_solver.py b/pybamm/solvers/base_solver.py
index 33f94d65d0..ce8c115ea6 100644
--- a/pybamm/solvers/base_solver.py
+++ b/pybamm/solvers/base_solver.py
@@ -194,7 +194,6 @@ def report(string):
             else:
                 # Process with CasADi
                 report(f"Converting {name} to CasADi")
-                print(u_casadi)
                 func = func.to_casadi(t_casadi, y_casadi, u=u_casadi)
                 if use_jacobian:
                     report(f"Calculating jacobian for {name} using CasADi")
diff --git a/tests/unit/test_parameters/test_parameters_cli.py b/tests/unit/test_parameters/test_parameters_cli.py
index 469a100bd0..85584b7c84 100644
--- a/tests/unit/test_parameters/test_parameters_cli.py
+++ b/tests/unit/test_parameters/test_parameters_cli.py
@@ -119,3 +119,11 @@ def test_list_params(self):
         # ./input/parameters/lithium-ion/cathodes/tmp_dir
         # but must not intefere with existing input dir if it exists
         # in the current dir...
+
+if __name__ == "__main__":
+    print("Add -v for more debug output")
+    import sys
+
+    if "-v" in sys.argv:
+        debug = True
+    unittest.main()
diff --git a/tests/unit/test_solvers/test_casadi_solver.py b/tests/unit/test_solvers/test_casadi_solver.py
index e63e88ce59..71632874f2 100644
--- a/tests/unit/test_solvers/test_casadi_solver.py
+++ b/tests/unit/test_solvers/test_casadi_solver.py
@@ -337,5 +337,4 @@ def test_model_solver_with_dvdt(self):
     if "-v" in sys.argv:
         debug = True
     pybamm.settings.debug_mode = True
-    pybamm.set_logging_level("DEBUG")
     unittest.main()

From c8eeb81ebf64f70c956471d27eb5a1d8254c22d4 Mon Sep 17 00:00:00 2001
From: Martin Robinson <martinjrobins@gmail.com>
Date: Sat, 7 Mar 2020 07:01:36 +0000
Subject: [PATCH 13/30] #858 style fixes, remove d_dt helper function, fix bugs

---
 pybamm/discretisations/discretisation.py      |  4 ---
 pybamm/expression_tree/binary_operators.py    |  6 ++--
 pybamm/expression_tree/functions.py           |  3 +-
 .../expression_tree/independent_variable.py   |  1 +
 pybamm/expression_tree/state_vector.py        | 28 +++++++++++++------
 pybamm/expression_tree/symbol.py              | 16 ++++++-----
 pybamm/expression_tree/unary_operators.py     | 26 -----------------
 pybamm/expression_tree/variable.py            | 16 ++++++-----
 pybamm/models/base_model.py                   |  4 +--
 pybamm/solvers/base_solver.py                 |  2 +-
 .../test_discretisation.py                    |  4 +--
 tests/unit/test_solvers/test_casadi_solver.py |  2 +-
 12 files changed, 50 insertions(+), 62 deletions(-)

diff --git a/pybamm/discretisations/discretisation.py b/pybamm/discretisations/discretisation.py
index bcfe9949da..847f5b662f 100644
--- a/pybamm/discretisations/discretisation.py
+++ b/pybamm/discretisations/discretisation.py
@@ -1030,10 +1030,6 @@ def check_initial_conditions_rhs(self, model):
         y0 = model.concatenated_initial_conditions
         # Individual
         for var in model.rhs.keys():
-            print('rhs')
-            print(model.rhs[var])
-            print('init')
-            print(model.initial_conditions[var])
             assert (
                 model.rhs[var].shape == model.initial_conditions[var].shape
             ), pybamm.ModelError(
diff --git a/pybamm/expression_tree/binary_operators.py b/pybamm/expression_tree/binary_operators.py
index 90ea3c41c0..eec915f832 100644
--- a/pybamm/expression_tree/binary_operators.py
+++ b/pybamm/expression_tree/binary_operators.py
@@ -13,7 +13,7 @@ def is_scalar_zero(expr):
     Utility function to test if an expression evaluates to a constant scalar zero
     """
     if expr.is_constant():
-        result = expr.evaluate_ignoring_errors()
+        result = expr.evaluate_ignoring_errors(t=None)
         return isinstance(result, numbers.Number) and result == 0
     else:
         return False
@@ -24,7 +24,7 @@ def is_matrix_zero(expr):
     Utility function to test if an expression evaluates to a constant matrix zero
     """
     if expr.is_constant():
-        result = expr.evaluate_ignoring_errors()
+        result = expr.evaluate_ignoring_errors(t=None)
         return (issparse(result) and result.count_nonzero() == 0) or (
             isinstance(result, np.ndarray) and np.all(result == 0)
         )
@@ -37,7 +37,7 @@ def is_scalar_one(expr):
     Utility function to test if an expression evaluates to a constant scalar one
     """
     if expr.is_constant():
-        result = expr.evaluate_ignoring_errors()
+        result = expr.evaluate_ignoring_errors(t=None)
         return isinstance(result, numbers.Number) and result == 1
     else:
         return False
diff --git a/pybamm/expression_tree/functions.py b/pybamm/expression_tree/functions.py
index c21b473ef4..0e96a753e7 100644
--- a/pybamm/expression_tree/functions.py
+++ b/pybamm/expression_tree/functions.py
@@ -164,7 +164,8 @@ def evaluate(self, t=None, y=None, y_dot=None, u=None, known_evals=None):
                 known_evals[self.id] = self._function_evaluate(evaluated_children)
             return known_evals[self.id], known_evals
         else:
-            evaluated_children = [child.evaluate(t, y, y_dot, u) for child in self.children]
+            evaluated_children = [child.evaluate(t, y, y_dot, u)
+                                  for child in self.children]
             return self._function_evaluate(evaluated_children)
 
     def _evaluate_for_shape(self):
diff --git a/pybamm/expression_tree/independent_variable.py b/pybamm/expression_tree/independent_variable.py
index 3bc842eebf..74a78cff4f 100644
--- a/pybamm/expression_tree/independent_variable.py
+++ b/pybamm/expression_tree/independent_variable.py
@@ -40,6 +40,7 @@ def _jac(self, variable):
         else:
             return pybamm.Scalar(0)
 
+
 class Time(IndependentVariable):
     """A node in the expression tree representing time
 
diff --git a/pybamm/expression_tree/state_vector.py b/pybamm/expression_tree/state_vector.py
index 259f19336c..8e40b23400 100644
--- a/pybamm/expression_tree/state_vector.py
+++ b/pybamm/expression_tree/state_vector.py
@@ -121,9 +121,9 @@ def _base_evaluate(self, t=None, y=None, y_dot=None, u=None):
 
     def _jac_diff_vector(self, variable):
         """
-        Differentiate a slice of a StateVector of size m with respect to another
-        slice of a different StateVector of size n. This returns a (sparse) zero matrix of size
-        m x n
+        Differentiate a slice of a StateVector of size m with respect to another slice
+        of a different StateVector of size n. This returns a (sparse) zero matrix of
+        size m x n
 
         Parameters
         ----------
@@ -255,13 +255,19 @@ def _base_evaluate(self, t=None, y=None, y_dot=None, u=None):
             out = out[:, np.newaxis]
         return out
 
-    def _jac(self, variable):
-        if variable.id == pybamm.t.id:
+    def diff(self, variable):
+        if variable.id == self.id:
+            return pybamm.Scalar(1)
+        elif variable.id == pybamm.t.id:
             return StateVectorDot(*self._y_slices, name=self.name + "'",
                                   domain=self.domain,
                                   auxiliary_domains=self.auxiliary_domains,
                                   evaluation_array=self.evaluation_array)
-        elif isinstance(variable, pybamm.StateVector):
+        else:
+            return pybamm.Scalar(0)
+
+    def _jac(self, variable):
+        if isinstance(variable, pybamm.StateVector):
             return self._jac_same_vector(variable)
         elif isinstance(variable, pybamm.StateVectorDot):
             return self._jac_diff_vector(variable)
@@ -316,11 +322,17 @@ def _base_evaluate(self, t=None, y=None, y_dot=None, u=None):
             out = out[:, np.newaxis]
         return out
 
-    def _jac(self, variable):
-        if variable.id == pybamm.t.id:
+    def diff(self, variable):
+        if variable.id == self.id:
+            return pybamm.Scalar(1)
+        elif variable.id == pybamm.t.id:
             raise pybamm.ModelError(
                 "cannot take second time derivative of a state vector"
             )
+        else:
+            return pybamm.Scalar(0)
+
+    def _jac(self, variable):
         if isinstance(variable, pybamm.StateVectorDot):
             return self._jac_same_vector(variable)
         elif isinstance(variable, pybamm.StateVector):
diff --git a/pybamm/expression_tree/symbol.py b/pybamm/expression_tree/symbol.py
index 06b9693557..79496fa9c5 100644
--- a/pybamm/expression_tree/symbol.py
+++ b/pybamm/expression_tree/symbol.py
@@ -524,7 +524,8 @@ def _base_evaluate(self, t=None, y=None, y_dot=None, u=None):
             array with state values to evaluate when solving (default None)
 
         y_dot : numpy.array, optional
-            array with time derivatives of state values to evaluate when solving (default None)
+            array with time derivatives of state values to evaluate when solving
+            (default None)
 
         """
         raise NotImplementedError(
@@ -547,7 +548,8 @@ def evaluate(self, t=None, y=None, y_dot=None, u=None, known_evals=None):
         y : numpy.array, optional
             array with state values to evaluate when solving (default None)
         y_dot : numpy.array, optional
-            array with time derivatives of state values to evaluate when solving (default None)
+            array with time derivatives of state values to evaluate when solving
+            (default None)
         u : dict, optional
             dictionary of inputs to use when solving (default None)
         known_evals : dict, optional
@@ -604,11 +606,12 @@ def is_constant(self):
         # do the search, return true if no relevent nodes are found
         return not any((isinstance(n, search_types)) for n in self.pre_order())
 
-    def evaluate_ignoring_errors(self):
+    def evaluate_ignoring_errors(self, t=0):
         """
         Evaluates the expression. If a node exists in the tree that cannot be evaluated
-        as a scalar or vector (e.g. Time, Parameter, Variable, StateVector, InputParameter),
-        then None is returned. Otherwise the result of the evaluation is given
+        as a scalar or vector (e.g. Time, Parameter, Variable, StateVector,
+        InputParameter), then None is returned. Otherwise the result of the evaluation
+        is given
 
         See Also
         --------
@@ -616,7 +619,7 @@ def evaluate_ignoring_errors(self):
 
         """
         try:
-            result = self.evaluate(u="shape test")
+            result = self.evaluate(t=t, u="shape test")
         except NotImplementedError:
             # return None if NotImplementedError is raised
             # (there is a e.g. Parameter, Variable, ... in the tree)
@@ -713,7 +716,6 @@ def shape(self):
         try:
             y = np.linspace(0.1, 0.9, int(1e4))
             evaluated_self = self.evaluate(0, y, y, u="shape test")
-            print('evaluated self is ',evaluated_self)
         # If that fails, fall back to calculating how big y should really be
         except ValueError:
             state_vectors_in_node = [
diff --git a/pybamm/expression_tree/unary_operators.py b/pybamm/expression_tree/unary_operators.py
index a4c7f8574b..dab197fede 100644
--- a/pybamm/expression_tree/unary_operators.py
+++ b/pybamm/expression_tree/unary_operators.py
@@ -763,32 +763,6 @@ def __init__(self, child, side):
         super().__init__("boundary flux", child, side)
 
 
-
-def d_dt(expression):
-    """
-    convenience function for taking the time derivative of an expression
-
-    Note that this operator is different to the other unary operators in that it is
-    *not* lazily evaluated, it instead returns the expression tree that is the time
-    derivative of the input
-
-    Parameters
-    ----------
-
-    expression : :class:`Symbol`
-        the time derivative will be performed on this sub-expression
-
-    Returns
-    -------
-
-    :class:`Symbol`
-        the time derivative of ``expression``
-    """
-
-    return expression.jac(pybamm.t, clear_domain=False)
-
-
-
 #
 # Methods to call Gradient, Divergence, Laplacian and Gradient_Squared
 #
diff --git a/pybamm/expression_tree/variable.py b/pybamm/expression_tree/variable.py
index 8650c3e0c3..6568b7c65e 100644
--- a/pybamm/expression_tree/variable.py
+++ b/pybamm/expression_tree/variable.py
@@ -47,6 +47,7 @@ def _evaluate_for_shape(self):
             self.domain, self.auxiliary_domains
         )
 
+
 class Variable(VariableBase):
     """A node in the expression tree represending a dependent variable
 
@@ -70,19 +71,21 @@ class Variable(VariableBase):
 
     *Extends:* :class:`Symbol`
     """
+
     def __init__(self, name, domain=None, auxiliary_domains=None):
         super().__init__(name, domain=domain, auxiliary_domains=auxiliary_domains)
 
-    def _jac(self, variable):
+    def diff(self, variable):
         if variable.id == self.id:
             return pybamm.Scalar(1)
         elif variable.id == pybamm.t.id:
-            return pybamm.VariableDot(self.name+"'",
+            return pybamm.VariableDot(self.name + "'",
                                       domain=self.domain,
                                       auxiliary_domains=self.auxiliary_domains)
         else:
             return pybamm.Scalar(0)
 
+
 class VariableDot(VariableBase):
     """
     A node in the expression tree represending the time derviative of a dependent
@@ -124,7 +127,7 @@ def get_variable(self):
                         domain=self._domain,
                         auxiliary_domains=self._auxiliary_domains)
 
-    def _jac(self, variable):
+    def diff(self, variable):
         if variable.id == self.id:
             return pybamm.Scalar(1)
         elif variable.id == pybamm.t.id:
@@ -195,12 +198,11 @@ def _base_evaluate(self, t=None, y=None, y_dot=None, u=None):
         except KeyError:
             raise KeyError("External variable '{}' not found".format(self.name))
 
-    def _jac(self, variable):
+    def diff(self, variable):
         if variable.id == self.id:
             return pybamm.Scalar(1)
         elif variable.id == pybamm.t.id:
-            raise pybamm.ModelError("cannot take time derivative of an external variable")
+            raise pybamm.ModelError(
+                "cannot take time derivative of an external variable")
         else:
             return pybamm.Scalar(0)
-
-
diff --git a/pybamm/models/base_model.py b/pybamm/models/base_model.py
index 6f90a0d266..a7285fada7 100644
--- a/pybamm/models/base_model.py
+++ b/pybamm/models/base_model.py
@@ -412,8 +412,8 @@ def check_for_time_derivatives(self):
             for node in eq.pre_order():
                 if isinstance(node, pybamm.VariableDot):
                     raise pybamm.ModelError(
-                        "time derivative of variable found ({}) in algebraic equation {}"
-                        .format(node, key)
+                        "time derivative of variable found ({}) in algebraic"
+                        "equation {}".format(node, key)
                     )
                 if isinstance(node, pybamm.StateVectorDot):
                     raise pybamm.ModelError(
diff --git a/pybamm/solvers/base_solver.py b/pybamm/solvers/base_solver.py
index ce8c115ea6..efc82253f3 100644
--- a/pybamm/solvers/base_solver.py
+++ b/pybamm/solvers/base_solver.py
@@ -292,7 +292,7 @@ def report(string):
             model.residuals_eval = residuals_eval
             model.jacobian_eval = jacobian_eval
             y0_guess = y0.flatten()
-            model.y0 = self.calculate_consistent_state(model, 0, y0_guess,inputs)
+            model.y0 = self.calculate_consistent_state(model, 0, y0_guess, inputs)
         else:
             # can use DAE solver to solve ODE model
             model.residuals_eval = Residuals(rhs, "residuals", model)
diff --git a/tests/unit/test_discretisations/test_discretisation.py b/tests/unit/test_discretisations/test_discretisation.py
index 0e9bbe9ba2..1e695d7347 100644
--- a/tests/unit/test_discretisations/test_discretisation.py
+++ b/tests/unit/test_discretisations/test_discretisation.py
@@ -704,7 +704,7 @@ def test_process_model_ode(self):
         # test that any time derivatives of variables in rhs raises an
         # error
         model = pybamm.BaseModel()
-        model.rhs = {c: pybamm.div(N) + pybamm.d_dt(c), T: pybamm.div(q), S: pybamm.div(p)}
+        model.rhs = {c: pybamm.div(N) + c.diff(pybamm.t), T: pybamm.div(q), S: pybamm.div(p)}
         model.initial_conditions = {
             c: pybamm.Scalar(2),
             T: pybamm.Scalar(5),
@@ -821,7 +821,7 @@ def test_process_model_dae(self):
         # error
         model = pybamm.BaseModel()
         model.rhs = {c: pybamm.div(N)}
-        model.algebraic = {d: d - 2 * pybamm.d_dt(c)}
+        model.algebraic = {d: d - 2 * c.diff(pybamm.t)}
         model.initial_conditions = {d: pybamm.Scalar(6), c: pybamm.Scalar(3)}
         model.boundary_conditions = {
             c: {"left": (0, "Neumann"), "right": (0, "Neumann")}
diff --git a/tests/unit/test_solvers/test_casadi_solver.py b/tests/unit/test_solvers/test_casadi_solver.py
index 71632874f2..653650430c 100644
--- a/tests/unit/test_solvers/test_casadi_solver.py
+++ b/tests/unit/test_solvers/test_casadi_solver.py
@@ -314,7 +314,7 @@ def test_model_solver_with_dvdt(self):
         var1 = pybamm.Variable("var1", domain="negative electrode")
         var2 = pybamm.Variable("var2", domain="negative electrode")
         model.rhs = {var1: -2 * var1 * pybamm.t}
-        model.algebraic = {var2: var2 - pybamm.d_dt(var1)}
+        model.algebraic = {var2: var2 - var1.diff(pybamm.t)}
         model.initial_conditions = {var1: 1, var2: 0}
         pybamm.make_semi_explicit(model)
         disc = get_discretisation_for_testing()

From 06ccfe87bcfcafd92e02091655c2c074983f5239 Mon Sep 17 00:00:00 2001
From: Martin Robinson <martinjrobins@gmail.com>
Date: Sat, 7 Mar 2020 08:27:39 +0000
Subject: [PATCH 14/30] #858 fixes for diff

---
 pybamm/__init__.py               | 3 ++-
 pybamm/expression_tree/symbol.py | 8 ++++++--
 2 files changed, 8 insertions(+), 3 deletions(-)

diff --git a/pybamm/__init__.py b/pybamm/__init__.py
index f2d3d14bc6..5e9790389c 100644
--- a/pybamm/__init__.py
+++ b/pybamm/__init__.py
@@ -108,6 +108,7 @@ def version(formatted=False):
 )
 from .expression_tree.scalar import Scalar
 from .expression_tree.variable import Variable, ExternalVariable, VariableDot
+from .expression_tree.variable import VariableBase
 from .expression_tree.independent_variable import (
     IndependentVariable,
     Time,
@@ -115,7 +116,7 @@ def version(formatted=False):
 )
 from .expression_tree.independent_variable import t
 from .expression_tree.vector import Vector
-from .expression_tree.state_vector import StateVector, StateVectorDot
+from .expression_tree.state_vector import StateVectorBase, StateVector, StateVectorDot
 
 from .expression_tree.exceptions import (
     DomainError,
diff --git a/pybamm/expression_tree/symbol.py b/pybamm/expression_tree/symbol.py
index 79496fa9c5..eb863f0d88 100644
--- a/pybamm/expression_tree/symbol.py
+++ b/pybamm/expression_tree/symbol.py
@@ -485,6 +485,12 @@ def diff(self, variable):
             return pybamm.Scalar(1)
         elif any(variable.id == x.id for x in self.pre_order()):
             return self._diff(variable)
+        elif variable.id == pybamm.t.id and \
+            any(
+                isinstance(x, (pybamm.VariableBase, pybamm.StateVectorBase))
+                for x in self.pre_order()
+            ):
+            return self._diff(variable)
         else:
             return pybamm.Scalar(0)
 
@@ -770,5 +776,3 @@ def test_shape(self):
             self.shape_for_testing
         except ValueError as e:
             raise pybamm.ShapeError("Cannot find shape (original error: {})".format(e))
-
-

From b3416169dc7abe9d9cc189ca9ba1332e2fbdd92f Mon Sep 17 00:00:00 2001
From: Martin Robinson <martinjrobins@gmail.com>
Date: Sat, 14 Mar 2020 06:14:26 +0000
Subject: [PATCH 15/30] #858 add vectordot and state_vectordot to api docs

---
 docs/source/expression_tree/state_vector.rst | 3 +++
 docs/source/expression_tree/variable.rst     | 4 ++++
 pybamm/expression_tree/symbol.py             | 8 ++++----
 3 files changed, 11 insertions(+), 4 deletions(-)

diff --git a/docs/source/expression_tree/state_vector.rst b/docs/source/expression_tree/state_vector.rst
index 4a0792b3b7..4370f6e980 100644
--- a/docs/source/expression_tree/state_vector.rst
+++ b/docs/source/expression_tree/state_vector.rst
@@ -3,3 +3,6 @@ State Vector
 
 .. autoclass:: pybamm.StateVector
   :members:
+
+.. autoclass:: pybamm.StateVectorDot
+  :members:
diff --git a/docs/source/expression_tree/variable.rst b/docs/source/expression_tree/variable.rst
index 0e610f5f28..510e366cf1 100644
--- a/docs/source/expression_tree/variable.rst
+++ b/docs/source/expression_tree/variable.rst
@@ -4,5 +4,9 @@ Variable
 .. autoclass:: pybamm.Variable
   :members:
 
+.. autoclass:: pybamm.VariableDot
+  :members:
+
 .. autoclass:: pybamm.ExternalVariable
   :members:
+
diff --git a/pybamm/expression_tree/symbol.py b/pybamm/expression_tree/symbol.py
index eb863f0d88..5c9bdb2ac2 100644
--- a/pybamm/expression_tree/symbol.py
+++ b/pybamm/expression_tree/symbol.py
@@ -486,10 +486,10 @@ def diff(self, variable):
         elif any(variable.id == x.id for x in self.pre_order()):
             return self._diff(variable)
         elif variable.id == pybamm.t.id and \
-            any(
-                isinstance(x, (pybamm.VariableBase, pybamm.StateVectorBase))
-                for x in self.pre_order()
-            ):
+                any(
+                    isinstance(x, (pybamm.VariableBase, pybamm.StateVectorBase))
+                    for x in self.pre_order()
+                ):
             return self._diff(variable)
         else:
             return pybamm.Scalar(0)

From 6ca684d3b6fc541e041128bc65280cadaa2dbb89 Mon Sep 17 00:00:00 2001
From: Martin Robinson <martinjrobins@gmail.com>
Date: Sat, 14 Mar 2020 06:18:32 +0000
Subject: [PATCH 16/30] #858 remove make_semi_explicit

---
 pybamm/__init__.py                           |  1 -
 tests/unit/test_expression_tree/test_d_dt.py | 69 ++++++++++++++++++++
 2 files changed, 69 insertions(+), 1 deletion(-)
 create mode 100644 tests/unit/test_expression_tree/test_d_dt.py

diff --git a/pybamm/__init__.py b/pybamm/__init__.py
index 5e9790389c..bbb9413b5c 100644
--- a/pybamm/__init__.py
+++ b/pybamm/__init__.py
@@ -154,7 +154,6 @@ def version(formatted=False):
 from .models import standard_variables
 from .models.event import Event
 from .models.event import EventType
-from .models.make_semi_explicit import make_semi_explicit
 
 # Battery models
 from .models.full_battery_models.base_battery_model import BaseBatteryModel
diff --git a/tests/unit/test_expression_tree/test_d_dt.py b/tests/unit/test_expression_tree/test_d_dt.py
new file mode 100644
index 0000000000..475cb733d8
--- /dev/null
+++ b/tests/unit/test_expression_tree/test_d_dt.py
@@ -0,0 +1,69 @@
+#
+# Tests for the Scalar class
+#
+import pybamm
+
+import unittest
+import numpy as np
+
+
+class TestDDT(unittest.TestCase):
+    def test_time_derivative(self):
+        a = pybamm.Scalar(5).diff(pybamm.t)
+        self.assertIsInstance(a, pybamm.Scalar)
+        self.assertEqual(a.value, 0)
+
+        a = pybamm.t.diff(pybamm.t)
+        self.assertIsInstance(a, pybamm.Scalar)
+        self.assertEqual(a.value, 1)
+
+        a = (pybamm.t**2).diff(pybamm.t)
+        self.assertEqual(a.id, (2 * pybamm.t ** 1 * 1).id)
+        self.assertEqual(a.simplify().id, (2 * pybamm.t).id)
+        self.assertEqual(a.evaluate(t=1), 2)
+
+        a =(2 + pybamm.t**2).diff(pybamm.t)
+        self.assertEqual(a.simplify().id, (2*pybamm.t).id)
+        self.assertEqual(a.evaluate(t=1), 2)
+
+    def test_time_derivative_of_variable(self):
+
+        a = (pybamm.Variable('a')).diff(pybamm.t)
+        self.assertIsInstance(a, pybamm.VariableDot)
+        self.assertEqual(a.name, "a'")
+
+        p = pybamm.Parameter('p')
+        a = (1 + p*pybamm.Variable('a')).diff(pybamm.t).simplify()
+        self.assertIsInstance(a, pybamm.Multiplication)
+        self.assertEqual(a.children[0].name, 'p')
+        self.assertEqual(a.children[1].name, "a'")
+
+        with self.assertRaises(pybamm.ModelError):
+            a = (pybamm.Variable('a')).diff(pybamm.t).diff(pybamm.t)
+
+        with self.assertRaises(pybamm.ModelError):
+            a = pybamm.ExternalVariable("a", 1).diff(pybamm.t)
+
+    def test_time_derivative_of_state_vector(self):
+
+        sv = pybamm.StateVector(slice(0, 10))
+        y_dot = np.linspace(0, 2, 19)
+
+        a = sv.diff(pybamm.t)
+        self.assertIsInstance(a, pybamm.StateVectorDot)
+        self.assertEqual(a.name[-1], "'")
+        np.testing.assert_array_equal(
+            a.evaluate(y_dot=y_dot), np.linspace(0, 1, 10)[:, np.newaxis]
+        )
+
+        with self.assertRaises(pybamm.ModelError):
+            a = (sv).diff(pybamm.t).diff(pybamm.t)
+
+if __name__ == "__main__":
+    print("Add -v for more debug output")
+    import sys
+
+    if "-v" in sys.argv:
+        debug = True
+    pybamm.settings.debug_mode = True
+    unittest.main()

From f3747facbe2f60ee1e7b767ea5a73cdc9ee37596 Mon Sep 17 00:00:00 2001
From: Martin Robinson <martinjrobins@gmail.com>
Date: Sat, 14 Mar 2020 06:32:25 +0000
Subject: [PATCH 17/30] #858 fix static analysis errors

---
 pybamm/expression_tree/independent_variable.py         | 2 +-
 tests/unit/test_discretisations/test_discretisation.py | 1 -
 2 files changed, 1 insertion(+), 2 deletions(-)

diff --git a/pybamm/expression_tree/independent_variable.py b/pybamm/expression_tree/independent_variable.py
index 74a78cff4f..008ddf6113 100644
--- a/pybamm/expression_tree/independent_variable.py
+++ b/pybamm/expression_tree/independent_variable.py
@@ -54,7 +54,7 @@ def new_copy(self):
         """ See :meth:`pybamm.Symbol.new_copy()`. """
         return Time()
 
-    def _base_evaluate(self, t, y=None, y_dot=None, u=None):
+    def _base_evaluate(self, t=None, y=None, y_dot=None, u=None):
         """ See :meth:`pybamm.Symbol._base_evaluate()`. """
         if t is None:
             raise ValueError("t must be provided")
diff --git a/tests/unit/test_discretisations/test_discretisation.py b/tests/unit/test_discretisations/test_discretisation.py
index 1e695d7347..1720711794 100644
--- a/tests/unit/test_discretisations/test_discretisation.py
+++ b/tests/unit/test_discretisations/test_discretisation.py
@@ -332,7 +332,6 @@ def test_process_symbol_base(self):
 
         # variable dot
         var_dot = pybamm.VariableDot("var'")
-        var_vec_dot = pybamm.VariableDot("var vec'", domain=["negative electrode"])
         var_dot_disc = disc.process_symbol(var_dot)
         self.assertIsInstance(var_dot_disc, pybamm.StateVectorDot)
         self.assertEqual(var_dot_disc.y_slices[0], disc.y_slices[var.id][0])

From 6da78cb208f726534510fc2270446f0be3140aab Mon Sep 17 00:00:00 2001
From: Martin Robinson <martinjrobins@gmail.com>
Date: Sun, 15 Mar 2020 06:27:18 +0000
Subject: [PATCH 18/30] #858 fix test failures

---
 .../test_solvers/test_external_variables.py          |  2 +-
 .../test_operations/test_convert_to_casadi.py        | 12 +++++++-----
 2 files changed, 8 insertions(+), 6 deletions(-)

diff --git a/tests/integration/test_solvers/test_external_variables.py b/tests/integration/test_solvers/test_external_variables.py
index cd2a43abb9..f8f796284f 100644
--- a/tests/integration/test_solvers/test_external_variables.py
+++ b/tests/integration/test_solvers/test_external_variables.py
@@ -50,7 +50,7 @@ def test_external_variables_SPMe(self):
         t = sim.solution.t[-1]
         y = sim.solution.y[:, -1]
         u = external_variables
-        sim.built_model.variables[var].evaluate(t, y, u)
+        sim.built_model.variables[var].evaluate(t, y, u=u)
         sim.solution[var](t)
 
 
diff --git a/tests/unit/test_expression_tree/test_operations/test_convert_to_casadi.py b/tests/unit/test_expression_tree/test_operations/test_convert_to_casadi.py
index 017cae7f7a..9dbd48492e 100644
--- a/tests/unit/test_expression_tree/test_operations/test_convert_to_casadi.py
+++ b/tests/unit/test_expression_tree/test_operations/test_convert_to_casadi.py
@@ -154,6 +154,7 @@ def myfunction(x, y):
     def test_convert_input_parameter(self):
         casadi_t = casadi.MX.sym("t")
         casadi_y = casadi.MX.sym("y", 10)
+        casadi_ydot = casadi.MX.sym("ydot", 10)
         casadi_us = {
             "Input 1": casadi.MX.sym("Input 1"),
             "Input 2": casadi.MX.sym("Input 2"),
@@ -165,18 +166,19 @@ def test_convert_input_parameter(self):
 
         # Input only
         self.assert_casadi_equal(
-            pybamm_u1.to_casadi(casadi_t, casadi_y, casadi_us), casadi_us["Input 1"]
+            pybamm_u1.to_casadi(casadi_t, casadi_y, casadi_ydot, casadi_us),
+            casadi_us["Input 1"]
         )
 
         # More complex
         expr = pybamm_u1 + pybamm_y
         self.assert_casadi_equal(
-            expr.to_casadi(casadi_t, casadi_y, casadi_us),
+            expr.to_casadi(casadi_t, casadi_y, casadi_ydot, casadi_us),
             casadi_us["Input 1"] + casadi_y,
         )
         expr = pybamm_u2 * pybamm_y
         self.assert_casadi_equal(
-            expr.to_casadi(casadi_t, casadi_y, casadi_us),
+            expr.to_casadi(casadi_t, casadi_y, casadi_ydot, casadi_us),
             casadi_us["Input 2"] * casadi_y,
         )
 
@@ -194,13 +196,13 @@ def test_convert_external_variable(self):
 
         # External only
         self.assert_casadi_equal(
-            pybamm_u1.to_casadi(casadi_t, casadi_y, casadi_us), casadi_us["External 1"]
+            pybamm_u1.to_casadi(casadi_t, casadi_y, u=casadi_us), casadi_us["External 1"]
         )
 
         # More complex
         expr = pybamm_u2 + pybamm_y
         self.assert_casadi_equal(
-            expr.to_casadi(casadi_t, casadi_y, casadi_us),
+            expr.to_casadi(casadi_t, casadi_y, u=casadi_us),
             casadi_us["External 2"] + casadi_y,
         )
 

From 40feb08eb26f55760ecee50e174f782e6ac9cd65 Mon Sep 17 00:00:00 2001
From: Martin Robinson <martinjrobins@gmail.com>
Date: Sun, 15 Mar 2020 06:30:38 +0000
Subject: [PATCH 19/30] #858 not testing solver anymore

---
 tests/unit/test_solvers/test_casadi_solver.py | 20 -------------------
 1 file changed, 20 deletions(-)

diff --git a/tests/unit/test_solvers/test_casadi_solver.py b/tests/unit/test_solvers/test_casadi_solver.py
index 653650430c..7319768462 100644
--- a/tests/unit/test_solvers/test_casadi_solver.py
+++ b/tests/unit/test_solvers/test_casadi_solver.py
@@ -309,26 +309,6 @@ def test_model_solver_with_non_identity_mass(self):
         np.testing.assert_allclose(solution.y[0], np.exp(0.1 * solution.t))
         np.testing.assert_allclose(solution.y[-1], 2 * np.exp(0.1 * solution.t))
 
-    def test_model_solver_with_dvdt(self):
-        model = pybamm.BaseModel()
-        var1 = pybamm.Variable("var1", domain="negative electrode")
-        var2 = pybamm.Variable("var2", domain="negative electrode")
-        model.rhs = {var1: -2 * var1 * pybamm.t}
-        model.algebraic = {var2: var2 - var1.diff(pybamm.t)}
-        model.initial_conditions = {var1: 1, var2: 0}
-        pybamm.make_semi_explicit(model)
-        disc = get_discretisation_for_testing()
-        disc.process_model(model)
-
-        # Solve
-        solver = pybamm.CasadiSolver(rtol=1e-8, atol=1e-8)
-        t_eval = np.linspace(0, 1, 100)
-        solution = solver.solve(model, t_eval)
-        np.testing.assert_array_equal(solution.t, t_eval)
-        np.testing.assert_allclose(solution.y[0], np.exp(-solution.t**2), rtol=1e-06)
-        np.testing.assert_allclose(solution.y[-1],
-                                -2 * solution.t * np.exp(-solution.t**2), rtol=1e-06)
-
 
 if __name__ == "__main__":
     print("Add -v for more debug output")

From 3cdb6dc8041fa3addfef36e2f1aea1e05a350dd8 Mon Sep 17 00:00:00 2001
From: Martin Robinson <martinjrobins@gmail.com>
Date: Sun, 15 Mar 2020 06:52:40 +0000
Subject: [PATCH 20/30] #858 fix simplify bug, fix expression tree notebook

---
 .../expression_tree/expression-tree.ipynb     |  17 ++++++++++++-----
 .../expression_tree/expression_tree2.png      | Bin 2563 -> 40853 bytes
 pybamm/expression_tree/binary_operators.py    |   4 ++++
 pybamm/expression_tree/symbol.py              |   2 ++
 .../test_operations/test_simplify.py          |  13 +++++++++++++
 5 files changed, 31 insertions(+), 5 deletions(-)

diff --git a/examples/notebooks/expression_tree/expression-tree.ipynb b/examples/notebooks/expression_tree/expression-tree.ipynb
index 9bcfeb3194..bf28c8fe2e 100644
--- a/examples/notebooks/expression_tree/expression-tree.ipynb
+++ b/examples/notebooks/expression_tree/expression-tree.ipynb
@@ -64,7 +64,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "We can also calculate the expression tree representing the gradient of the equation with respect to $t$ (which is of course simply the scalar value 1),"
+    "We can also calculate the expression tree representing the gradient of the equation with respect to $t$,"
    ]
   },
   {
@@ -84,7 +84,7 @@
     "![](expression_tree2.png)\n",
     "\n",
     "\n",
-    "...and evaluate this expression, which will again give 1."
+    "...and evaluate this expression,"
    ]
   },
   {
@@ -95,7 +95,7 @@
     {
      "data": {
       "text/plain": [
-       "1.0"
+       "array([[-11.]])"
       ]
      },
      "execution_count": 4,
@@ -104,7 +104,7 @@
     }
    ],
    "source": [
-    "diff_wrt_equation.evaluate(1, np.array([2]))"
+    "diff_wrt_equation.evaluate(t=1, y=np.array([2]), y_dot=np.array([2]))"
    ]
   },
   {
@@ -202,6 +202,13 @@
     "\n",
     "After the third stage, our expression tree is now able to be evaluated by one of the solver classes. Note that we have used a single equation above to illustrate the different types of expression trees in PyBaMM, but any given models will consist of many RHS or algebraic equations, along with boundary conditions. See [here](https://github.com/pybamm-team/PyBaMM/blob/master/examples/notebooks/add-model.ipynb) for more details of PyBaMM models."
    ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": []
   }
  ],
  "metadata": {
@@ -220,7 +227,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.7.3"
+   "version": "3.6.7"
   }
  },
  "nbformat": 4,
diff --git a/examples/notebooks/expression_tree/expression_tree2.png b/examples/notebooks/expression_tree/expression_tree2.png
index 949690cfc08cd224e5d082d65963e3bd4ae8ee9e..09f6e2900e7ecb4557528a9889e7479d6eec742e 100644
GIT binary patch
literal 40853
zcmce;2|Sfu`!>9#kTNEfF+(y{D)W#+L?WT2GAE=#2_aIMLy|~?GLOkrgeX%q2nl7D
zLWD9;`#V<M_y7IB-#dNpH~pUH`8{p>+SjnIwa)W6kK;I3;C?+VCPpqs3WdU?y+`8!
zg+lX|LRsR$K!?Aao`3WS|3_!FS4)FJCI5L?k{m&ytfy#e>@sq@H`IRd=*J6-^rOvX
zD-Q)M5qEp&c7=zNFO&01O5rgpvoM+ShP9vCPO0P?x^CjPFTC@(`FUUV#e`IrchZJZ
zyA>*4uX_8d=K#m*tvkalZFiKa&NFHA#Qa+Dda?bt`wqT2x38By&5K=PTt`OUdTMRF
z#=2fmFz&gqe2v~t4Go(W@=}2T0Y7QzMcCNtRX)67-uc;5_O8MH{g<y?S#q@0i`O=H
z-G&Vtr0(dBj*ebVPnT$GZ?8`=5RyLmY2}U`J930pv8DRaHqt~!Mhd8?@K%TMYuVb0
z=-rd|$vsim^Q-TT)ZKSe^XJ>z+CJnQ-@A0#G66}+m5&}h>TgI~p{AzxuyKH6o+%*U
zV{&i0&2Vg`)ytPx>9eh?nT16}1aZE;zP{11u~Z)CUAuPKeacuDBX@>bRaJF*db+Hp
zW;qj=^kuWG>UVeL49v_zk32iw{W<^S>$<w?%HTEOQBnT<-t%l5_3o9`L~OK_=BiK7
zUNt`0v|{VltpOn+yVDP5+8o-cG5(NUUTWAxM~Cs#<72^@N1h3)s`3wi&Q~AlC}HE}
z4e9UicjvbeiCC@`##gzxFqeI%bH~t^XDpSKmHnUdCHz*dyRydIlMVNf)ojhqmU@l0
zvE$b5+g43!#wO35#kOtR_M!0niIwZNExk&|q~_*UkaqaGrqJY2>(xz1o?86=_O8Hd
zL6L6d+RIM8Ro#IbAMBS|!B$V(Q<k?TF)`8Pg{Q)f`O)&)m*Y{T^UD@yCn9*6%7(iO
zFZ?*@IXjR%{j+F<m5(pfzA?$N|HC8C7ZCxDH1tWCxh9g+Bi_{DfyU(cxVR;~)nU~M
z+7Z$VjGxN=S$a#RI?A7(`tqT~!&TBWkD<NDHHtSlAZ2rLL)a!I&uJ<5@qVA>t3=js
z-pqXX@Zo5M^GQkm8z{~L4a+dZm5<c@ySlr5_*C32+6tU-_s@syG&J&Y$Gg-R-rm_x
zui`#_xu{4flt<x;_reeJ82($OtJ&CmR6HiTe%Gm1A2Pq>tHpP5X=-X}jx`%wQvhva
z702q;-8CCyEAs3cM_<zOuiv_r4Ra*Nx?S$Ro`1NImY|Rjt)-=9Kxk-H>Egte*6e_Q
zv~J%Oj(g0^_#U_qigB!2BeY>d<>LoaC)*0P;@SKz{{Ge-BG<9=;6YB#%|{tdoH)^1
z9}qw}Yw`K|_7m^di-^#5bacp{?V_0G*$F5eC}|vUj#j?7?eo)9!E8dBUA2)CCeKdu
zj=Oj43*58gh0Z6sqs4A);o;$SIKRvAXJ5?_uA0z<f`kNN<(aP$Yj=3c-`;ej>e&O<
za7WL`j8E;w7gy=4OmG}MdUW)iT&H}+lckr<O581bt3pgl++{jEMpRmz0s^Mv_%rWr
zJI;9S{Q2zj-}R(Vx0$%`7^^IdH<}cj6|jAOFR-FQ0|&l-Id&_<@6?y4OtkdO6;&bJ
zTWwymwmK!OlXYC%(9m!rJiJ?qTb9k9Icq%LSTmkYKtQ0sKA!$l#$h3)P@~<ucQY_D
zzQSL&o_NngPPcVTK{wMh$7cDj-@mtdOggq2B>bM7tQsF^%s$(#_Vw%6*P-ef%t5q<
zS0y=De0z7-ACpsH{ra6QJp0xvxk$#yo(jSedlRbm_aK*%j*zytqs6~!C352J{DOm5
z^z`&RacB}ZGBWyerU0DjQ<izy0~HTJd3h6uP&`Ig(a4L4m~K5MC#M3>nQfCl+a;73
z*%a0@>oV!+=rmN&OO}o|?2VaG-bb-(e=hU+*_p7fUu{|q0$$TNSvb6HN<ZW@Idph@
zsMXMgW1ss_cJ<k=GMd({f4|n%(__BIil@HGyE{VBWhh{1$nnp2V`-nb&QA;)Q!Zy^
zWerxpudTgOZ<jHD{^CUiO4s*T@5jF`ZSyb8%G$!XVnrfvYyi(1`lPB`k&!eM1U_Tc
z$wqzEkXWzT)X&4Em&U^1o;>E}?%snJ(h#bz5gMFqBUAO|3wPOY`j?d9uQd^@YDSB7
zilQUM`zC!K;f&qAy+SgyYziXGx*^#Qf=bTy1$`C%?A%woS}Zv|FhE{PCQ03&`Rw`g
z<(Tv~TJ;H7u0qPnM_pD5p4-V2&DA?N7!(k|kdTyA8ycse#*mPhIJQYnLwFf|65Yh`
zmo;CWp1RuVcsr(i_FmG#q0UnN=Euj>V|QG#(Gn1(ukcCTMF#t9_Y!t?cAWiI-<zm-
zkEx#wD^^7C`Udc@;w<^6zHk~C7(~mRc`zgz<Fq(`(eZQMa_g$=;p7Mo)t#pClOvsB
z{rz@BpYvBN%uPk{YCSI~kVm}kNsjfplAfMEB(L$+4M9~#QdXAp`}_OlEUNPqIXOAZ
zSM9xf)opBq4`!PBVwI7p5gJHq_~zHIU!$YtEY@FNhnqaH@EMXnfCw75H)^xpm#3Rv
zEKDBhcUEi-$F$urs-VlXf`Osz=7zmb?CJ$QXZoUYdQ;w8Wm|kqF%T6Kdx=;0-cnoR
zz?f+cmz?MHSsrCK)_BbjH7vaPXJ^b?!;OA_KAB}%9m*RJ7Umxt8%yrCtHAQ}<-+sd
ztq?=kO53sinTw%4*I3;aCYr|vl4GsfpI_*2N{<@aqoKhQ{gmV$i9?Sow$~|NO@83c
z_@g;1idTCp7M_q+XvHhv71@~k%2RFUbA;H~Xu8Vs*3n(1IobTUvU71#?)uuTOU{%`
zaEgnI+f=6=QczF`2n>8VIs9b|fns!Zwnl;&Pqb_WE^+?*`^uc-Z^(mc?B3ly)S8=P
z-4wj<uB?Ed;F8k$Un>q68rEl-t72#Z(r}esTS~^+8yXu84j%Ly9X<Q@p28~RPejqE
zj10+@Yo#dVeoW<0EX!{^co6#G!~VZfv@FlyzTyQncXve+lt}zWpf$<04H;kYv@;z~
z*_HXML(R@k^!<G${}PX>o>J<<7OaKiT0(+pZ1hP?+t;p)R=gmt_dq$Yn%|oTGmL?W
zsgk-lrP}uV!d{nRF5JVPP>ji$eXdXEZH&wc<kH>9uYb<5Dr?+Y>-{~2t1LUFDdwJI
z%am1A#(ILcSiOzj+Mw@r(8$OFvv_*EQJ>%fY5RtyS*44LSaS3%EG((5055D?LQaG+
zGBQ4~s$}$>Y(4S1wpQ)zSt*j5zd!I`W@BSJ>*8`ySzTLuIo@nAQnfYiMO0i|4NF7H
zWyoYDk38$<&6|TxtQSnv;)$MoYn8oz!-i#d<<3+gs^5r+xa!BmEx2tPI}!i^F0Wp_
z`cQoFH@#wLjm<0M*zh}d0`u~ug*R;IejP6K#6+^cHFqN&3*U8s86hdDCthA#1qCOZ
zEEI3PeXBdvR=7&yV8%`Z1Gde_o-d^U-|WJ4>>e2jyL_1<ed;p@6SwR&AT4EO<u0uE
zn>X2`bwzI7y7jQ2VC}_=7ez%ycb-1IS?|8$viN$260<PbQ=cDdM~FnrI|;_c#kpxK
z3FfR>UVJ-)cKM1G<w^VQ-iVD2L->n+XKQ0)(^V0;O8)$}mALe^c&)IwhYu-0GNN0y
zREKcOSx9qr4-727`tmXz6Ia=rH!C)o=3Oo=RTUHzWaZ(xjyZGf+O@79S*2B@-Q{dN
zJa4hmTzt8cXt4~<KP2tr*1adsMxCx!<xucrT@t4Ji$U?ipiOhssF0w4dDzrBVZfcM
zz<ASRy{nb{^peYw@1HpJYT~{`U3dfq^PS3_Ssfgs4z_Ize9epx9RGgc!v~7P>A&B+
zKU?m+^xqF1XUp@>Z+p5|;`AL28!aqkrFCBtS3lDg%-Pr`H5~nRxFP6>2wUs&)6tW<
z-(2h~H$?;tEtAt#O7#8b8ybA0v*MrM(Bn5ZZvS((UYt#6)2Y7~{p3`hx81aDb7sMD
z_Hd?NKl?Lz1;<RC-UtgGsbXpjwYyl9*K=E2LpF#zNu~bR-wQde#)T~ZJ@+o0+jrAH
zPsfSV(LT-ldj{@urbfP&hJSv<B&W-LBKV&h5@?e;b$`=84}KfZYWcMIpC7$sYP{Y0
z@0IXK+{an|`4J6%6shv>mGGlolmA>P1UDW$WdHYGO{<w2t6F;hURL@#cT#3{_1{6U
zUx(9K)%`tPgXrmKQ<;;WB##!kWY(UV#Pyr<kMk+HtV}<Yb^F%!03Kc>j@5v0?XI1>
zckZMl>Boj5+I2^nd%e7}Y=zXtUpCSXANROu>}KLs<Un|t?z+lCY2&*E6riD@vGm~M
zqhboYYzpEJZZ<tc1bv9ePob2Uv8-9M#`(u5CK4_5^g`QNbSW(bPC|Yx{GrG{$T(_Q
z%>>?Nn&zwF`yn^v9XFkW!8tiO;?!1~ZQVdFd-m<C1gIbF_UAvSWPI|>w(Z-;nhu%o
zV(<ZUV-VPHU_e7b#o%9Q6@KqV01pBh-SOkcUshHo)}8|BM@WzL{Iwh)ks#uS4<GJ4
zpmyBG*47HCs=uvpyX@&U27#rzB8HJ8Tkq)NL{YbJmOExHl(fxCN=2WtjE#&Q;ucqM
zN|pgQKSA*zucVZ|Mo&PMI=_8!)L->sMn*=h=HUp0=sd^P1C`Di8f{WC0|@p~Tp{^h
z?G{q&m6dq_yz_%Cms3y%wY<2bboAJ<AM^L`-Mg1|FGoc1zzW>;{T-K{f67|7k`7Ch
z;O7@FUWg|hFw@mt33SFJC-=<MKk21K%6Zz2f(ICJ{IP<wkE*^f(a}D484f}oFR!lF
z>N6U2n;)wtB@SSXT2|KLb)QzpZ)L)QVo2xTP|b--NTj<InPq8B{cLCB=V#fnWy|uK
zkB#{!58-uul9%2mvI>D|gV;JYg>Aq%sRJ_#h)@JAbd+$TT)S~Y=MsqvBy*vPT)Waw
z++4}^=e-0Yjjyk94G|?JZpHhLf7QKz|Gxjrvn@~0bTYfTx(3#KY)&=YD7<kaZZo;|
zgQTcv%5UR54Go>HvX9Y9omW$ve~eYDF0x^qeDB`9`_<+gDK$53V)F9xid1hzJm0cK
z|J*jgbJCvtwkq51z9Ug|*DjjJ=B3QJt=Ip?ST{GfvZkihKY#wT`<TKOBXja#?+VwP
z++5U--9tkpm|LSTlzA~-8NQT`PEbs2d5QbPD*$0s%;AxdzNkA6&b*lIj}1tX4PDsd
z8CzRhOL8sBFQM(*SCLAzHTQ(m%&+6o@4k5~Dx5xj+O9Qc!<8#nSUEZU5M6n&rqmdU
zM_;Wl$+q%sFRuHwYZ)#iBeNQi?DhNinx3A@$g_t>7F<z`0puI&1+QYGDW0BH$u^Oc
zcN<;rGSuROk|7~EIlj1<2L*a0OWywrrE~YUZ+-;mT$&x<<>)ABn(wfoXi4hNA=Mo_
z_yz_BSk**7c}Ixs`;dNUOGQP+6K}yByqFC2|7lU}y*@J~r7XoDUcB>UA!ZzExOha3
zwcBhM%sj_-OP}c|HVaEVo<G5l;&Fa{zBY`Sw&iafo;lP}KesVXTwcmMJui;~mHexQ
z#Jxz)J^(4!uWtwt5cjIQe6`9`*0|c#o$>Vns6H*;+}_kud~rM8c_<An9jRX3dIC3c
zyc2v=GD-dN_*PfI!YYGI_Vs%h1T?g@1EW>FUyM#ot-))uYI$NQC?K%QB|`^j#tJwp
z&b-u1Ex9oe7Ys(eW(cknR5C5GplT#IjWMfAy|c|V$*~C>otRi{W@cvfInO@s$&>vq
z8H&7$&U-4Vsu=I<ZREQ&ZZNhm)%hmWxwA!;t(BpgFWR4(mmT%^8oH~yiWg^p#CShk
z%tQ?k+)?7OV&le*V}t2gL2+?{t>Ip4CiOLBE4I=nRa&v}Um6dWTDwwETwDjSsTSvx
z<P1ignNMyX(T`J1Oe{yJef8W%LGG_^w7<cU#CcN}Us*^6D4hS+^6vcpLywQL)hYkJ
z{^7@j6|h)N@5i8pRpBA$GJNjby(^)nu8z2V`PhyTn{1Pcw!-tixw*M(<xbxNEtDf<
z{BkkywO^&u+ej}0<D#S2x_Nl?mif@y)yMe`?U~R^ZnXUMwT57M&FgD@04Xg1r(+88
z_;5c*zO#;$7rTuKqQtK+aIz*<MF^M7O<rv$ZEbD0Uo{)I0>w)nd8#{@*wIqjwiuXV
zCM;hUz17-ph(qW3#or>c;{&vmmJ*L=<wwsJUHI`SARwSjX{4Aoxbg}{ZEwxX#Xbzx
z#J79x-o?$g#0+oTy&HUA$u$loElP($p!+Hy_O0i?=^}*b)Rg@_y!x?z&FTM1z5Qr#
zVE?4|-3)R+tt0t1A(}#<vp&RYv1A<1)pX%W4GIck<>0thAFpMNQAy*8SfJ`?QZ<Z5
z*oP)>4WWJ^unghgO{n1gXV2tF&h%ZeMsVXstw`$}s(B5f$YAWT)6ty4%T9gHd->s!
z`c~JUM=L`k@Uvq@E|w^)UcY&x?%*(}vQV;G;-Cfo(BG0R1RVS!S)X5-I{!PUZ|UhX
zo|tiEv5VsyE1&V?&pzAWfWZiQQl~G`ntw8ivL4|ibvhgT!!im0ou_^$b#ZPaxU8(K
z_D3k^>eYc|Wjn{fBps=VR$~xYoT>9>1pS79xLPGNX5kBa{_IA<-4$Q@#9AkQG{>K)
zd%(cVd~n``hK6Qpz64+u<dH~KDDRP(o-i;b-I0gPUVb@K8uIFuTB6kw3Rr@(h=xP$
zMThEpQ(C=`G4se@9(ggR=isN5+N8=RR=%0u_P~LgBC1`k?_y+9Qc`wyzPOa$)`16X
z1H<wrYRjSVX-sg<xyj-1d-pQ>o$?x}d$*Ka^IDj#2{q8pTnyZ^&P7ft%d8{_d}JJM
zFel!~`qD$k)|}&!UowvjqNd=e8*Q-;>&Ao)TJYR7m9<P~N?)TkyU6EGL<FtCQiNw4
ztyB@Pf@Q;p_I1>_PuxYWuH5cC)WX)<+WKGrCSA?E#C<iT>&FiX1?>x1`jq{{&*MI+
za))SEZl9i=6`wtt`gSolCnwF>#LSGDGClJ>Hr+G!t>t59gbtn}?r5M73I!>zMmuc<
zGjkOfN2i6^3)3?*#vNuxMoTGn9mRz`>rf+M{xVg1y?DC_I>LS)<x$gV>Tqz_wF=h*
zV(l2a>YsWwruwd?{x@!4^YdFWJ@GN43Au`ff;bs;LR3=HXnaZkwmH>{zl5)`MpAQ`
zYn9GKC`i)uVO93Ew3p>QU@ZQ_PdvNx?pv+*;_q1Mb-Y7kk0UQZ<8<}(TtNvh`rtdL
zvY-c>n@p*u@&7&prrhPbce%j;#?_8xB(ERgxqdL$_MIA1g!mT!Z^OeO&pN%hf%10S
z+lx^iCL}DyGl-#7IXK*ut0NBzk_7L$8HcGeCvW|X4Y+0PO52Fcj7o0G^Uibn@aX7(
zz*QTUP>`I*vb-14YmcatYN@YpS;=B^a!{vRFH`Zvrz6up%&9@(M!Qj-RDDb}+<Ix&
zx$-xoz>V9tllqVB_&C)0DW#qa)#j`XO89ObGhN-qMfTqP`^T&6uV0#*eB?UP@%r`a
zohByHR@IhN8N+mM!2{PWd6{Y-J<7BBXu;)c*Os;#NFSg1)z{xqQq;rpX{glJVFL5#
z$s3>6r>BmfaJTFzz6e^ZYUcO%Z2Lx$>+7~_j#+0@2R{zs$t<IiZD=g0iwlI?e0;3v
zL#eknh?3yoV09}iL9f|?jRnu2i+d?gnn!*{b}HHX<crU8N_KX(yz5UMujy|3g2KX@
zP@#RMrd)`igm*CRcy3j>cqRhtr%7UtFB+@RZMxfk>$$IcBu$=x!V}ZkySF}Z@0vVN
zJ)K^2r}E$Y5(`{-GmL-}dFe%ewfOq=l)bklU&ieY9>u`bhQ=!*mKwxsJ|Zv5d$Dj3
zPicM{v8;6Hr)f^9%zsU@q+jApRPh+jZ+?Puy`sAMP@zj^%-_h(x18QW_HN~hEye4(
zr0v#&1P49YgIpE8?RXGsbj~AB9m9G3k$F@CTmyW>uCZ#y8?5BtaTH?^gwS3S8Xx`Z
z*SW1$6%5@yJ%S+YC6x3wkVLq-<6;#C*5*v11I|;`_3|oNHM;i|S7SY4XrHD_)PvSs
z4@Rn{rl#fj1qnK49#N$Oe}4S)-m%jD@6Q0Ir4V9~YgPtxZc;Nc;&}f+MMCLdNlTvn
zDwM#ZZ?=@K2NwIo!l!cKvQOVWqDyKb@n(G786=?qkUH39jf)oxhfJn6NJy-prCV7B
z#zo4pMGBbgkk`+ogV{hw+By1+Z)~KwD(@+rzm}YwY}fgsgyd_IPp#V{5Z8j`Iz0^w
zii(17++fXi9Wg`3mm?2a;bb;$7*$-Xyw%~O=%eNzyZ7(EVN!T55x_ixH(j%^u+VDo
zW9qBMMnf<OMsM%R>8}mZY|Z9lAi2NH=VEeWvVJg<(hAg8gpGt$BZ`{ncr`nQ+VRWN
zec*|9qC{EY?d?4_R?Q#GEqD5L=)o)i6nm56W9D2VUaL7cEGmOmPY>nR-2m$C>arl}
zMxDwODtJMoi^os=uX&amzIexUtyB3n{I@Etb^#~Kp6R#_L1g{Ljr4$18XyS?AiylW
z=Ii?iWIl^CPejc%+&n?j_V)IBbaY;UBelpj`H*`;e<h#t+N<;|niduo3-ea8S8#T}
zvlfedckbM2d3J`MFqYtszkdHtOY!vd{E%&}H9a>cgn14!Z8y+d_uwES*d$D|ILvW?
zOh{g$n>TAhUO+x2h`}&f@1>nIDLC;jOejrEk3Vcnz7i=1k$eH!mm%q&9cb7<nfT&;
zSN7Bg2u@%`!WUN6-<Ox-*RZq{fE;wTr(!u?J6+V~W56A2_Cspx>9M*M6$Oft8#UkP
z`1rM#7CWI9A)fhiZZa+F?BvJLXn@iJ5MyI!e^sYb9K3G3t;OQ}45|XZR>#Lvi=~UR
zC$Us`NU;o+X~;Y(Anz(|gk$xV()lZ(h%9<40>R8IK6<?syyn+%tq&-uGR#ZO!mU}u
z_vMU)O$p9zJJX?TBFTB3U)3zz8l0={U0Ei)e+^AdKdWj1!7t@Z5sp{5{;~pTG;+E*
zO#kpRDI#rB+P^>kWk>4IoGlj~Ha!2AWGJm<ijdw{^_OJ$S;iER_a7fB@Xg+TIGII$
zIo*cMcm5$Tdi=Cw68<4*e7HsIBRT(4GD~r0riQ;n$?I+>T&DgX)}+AlW|E!LpG$uI
zfDcdq{m_CBw*Rm*U!UNE%fBBA@Ij>I)t@_iEr64qdh+*@;nui*L+ihn!8K<b|57Gj
z-S8Wur+<l`bE0_Y%}4$|v;ybpXyK4eXa0~yJalrpH>+<y`$I~FUgZ|iw)}^h@y7ss
zed7OzWeK<L^3zVsP5%3#z1L5B@ZS#$_>h+Kmy8KNM!xy?!zJ=f_V<5&gZr$FmO06Z
z#Bs1s)}v!{=HTaTVZ2IyP)GLaXlk$o(MAj&1Jl@XY3kV0rAsXUmJ3{lxx8lj)*u;B
zA~x!4f*rTfGWrOB$%ZPFkg#RLyA;lLK^@kx=&Oyq*~nJE&?9y@%~dOuXBTEIh2mbq
z2fDQXQzmyPpGwfZdpt$_(G)`DkqKB|U%y9J_Zm2Y-Mg2eQuM|5tz!ZGQ@uEcHw48<
zs^hjhXMca&2N|R1Q|6IrQ2*)o)AHr8TBlLqtiQJNzybDatbltI3f647RrVgeJsMM(
z&4b6F=Rs3TYqQ;U2W;gSM8k)>gE^z6ESF%R3U4sdw6zt{4CNtpnKc#^t6G`OHscvS
z#E$uSPwnM_{k!^))32Y4xPSlHk7iKtAb@<9t>Em!<C92%%-^&k`3ADwI`8Rn{{06J
z9y}hOp^dAWK6OmBQWQ*MX_HFrxUb^DsTrUBDT^C?i@(41k58E{TP11Mi;6OQd3MG>
zJ^keU(ujx%!$`wV8a7Ju<QM-Vdd5qIhT~r^z<*2C4EK~JxVv`0^5r@sw{|5FoDs`d
z`1$!G4O4p&o~vJ6n(G2(Tk>LdJD+Xt9omT3nobr6N}fYe1N)HH{<!%^J#=2=;qo_c
z_Nr~0a321=IxsL$hIt>x1hGt;Op_}ZTH}3nbpsRb8d_S`Ab0FOJoLF{H0XgqaU&)s
z2*B^9g_P0EXQ(8?iGp$TFMYHp#Q*omnuL8k(MDN6Jz0{V+4EDUwH(`@uSRiwGopqK
zl@bVaf^wPn0c;?@)zo0h^vp@6-fVT$1XF|8dlo4jEH<ccywsQF<+P6<KR)Z~s)ZT+
z0aX{5c>!K46G(BWy(A-mO~@Rp+!87*CB=$j5(3Ve&1S_6;Cl$!1~s*Is!5B$?rJaL
z^ZNl5_UP+ZkM&j)Iv;FJ7dSYgOun`KtwQC2fCYYz9b^jRGIi)7p}b0T0!z)@x|d~`
zWWPGyabo7-5?b1D3|2~tsLhEJLf}R=>OaVsS}bxNxN_sh4e{j4z!b01o*&F*R@q&E
zC@tC6^z5rvg@Ymjs)47$XvYV-Wy>BN&b0-~R{-bH)!R#pvR~Tn16gx&XF32a6i&B4
z_r;?Lpdwkb#Y_V^%4>d1>z=&R5?t<Dz1=lHZz76A&7}uWv};O>9AZp(o=cscb^{1e
zXLw{n-4`Rh3q=v@s#WESfrI)>=#zN&YHPm)F^~;w0Tx|6YRvxjBDuKuc(@>rjj$wu
z;9V=@VASWaX=?b3KT;N%N$P$pUj}b^5dfyO`N!7$jHio2f`X`Q&1c5UCn`a~LUfL7
zS8}InOkDq886)DVIn|65R@+H)WuK_yLxuK7fI(fk)uN1Yom-AsUwQ0W++mx-=0)!>
zk?cGGQp+~%ym1V&AKUNmb!!$E7hk=5w{Io4tgf)I@as_H)Z*f!jxA4ijsIxYM>ea6
z4<p;MoX#@a#6~N9%~0&?pOxn=KFjB}O5T<{Ea-XZl5bR0RJ81=7!*ghc+>R=8G`{N
zG>;blB_mJj!W38(z_71U@=~W-!RKmVrTD(IV1**}{2PCY_YdeIaVQyH)zoMg&km-?
zrKId}N$3RHQQ%<KJwn(MoJn{k-PM4gAa=w%az0Q;2eZsJ%`YtMMjY9#qZ2f=$F*`X
zkcS;kjzX8=nEuYvQqVhum~r{d^D$E55K$h5#@$Pzx1R_dlFzrIYA^@WHs1dF?OPCN
zL^f7di>NxK!;+kYG5~cc31qQRS0rL+k1@D8f6h9Ig_eN->o$n{BAYX-7Xt%E%JXQ*
zi~~61k~@7JGvpDff*c{S&7eZ{cgzSAX8!n_+u9``x#yz3sD^9)2UBBI1Tig6KXwP@
zGUzR(x67V^sSSaYUHtpK5)sVWi=N%$`BPPJoeSk3Npvzy(%l7_%jIWVAnyK6XnY#v
zWv&%#HdjESL;)zef%eZeZjb&ioDu(psgd?;zS$%WG7y-0yuY6I<fn{qUad6nX&{R&
z@fR|j;?*|P(wp_F;dC2tX3b59bC%A`%tWhrJRK5kxai?wjTB1E8RyQOGg1(lqefm|
z49p|`hX0vULc&XO()@nR!YcdKpx%eUg#)O&^+g_n3cnIayrQzw5{4mUy#3=QL#2u`
za~m7PVT3R(^-^)490~-+ZVEeQV&Z0QZm3QwPxBYop~Z^3^Gup=2M6=#V7Ml3{-P;e
zwrrVQ^J88l)<eT5p)Ql80fi;ieRDeg7LTHKkDlJaarwW6YJxFr1;h|>4=D{KsnJ3U
zD04SV^*<VHW759JXW#q!1YnpUL=p@<S?-UOR3=-jVV=2cr5d}?yl(&_aD7eS$7p1F
z1_7{M!hO@=ytNRkK={!R6$haLwkw8}D=V7Ky`e_Pi4|rLTgI(h9F<->=CbZKU~JaG
zYH{>m`eai%3U&}pTTwa@vl&>PDv0FWKYwoC6cGvYMqFKAHuX8Ej_}(s$NceVTpp*d
z-fCS%luJA~0R><M624{s8XvEKyvZ>)GgNvS>`uaQ;-wfKX6l-3GPWsnpAes&_~8#*
zk7eyGv9${4gkQ}4RI<PL)5ydmrGMFnR_}JO5PRN0Te8kbUXpdR@Cvw!G_QHg9Obc!
zRc9|;s3fKP6H9eePC}?LJ3W7gjO%9>fIo>#-uv(G3oIRU5mUNArbFdQ@ft!99UtlB
zS7ShWj<+oLqg)Q#G0R|G9a@=Nw?krK(=vY+eocV+|1`|BhSm^&K!o^w{Akn^HVhgD
z5rt#<^5so>MVCOA5*fI-xY(%f_=amcO-*Cs&dgX)^9|#@1vm93N7{WrRv3eX&jRWS
zmGW&~t?d|9e=wy&e?3+~sK1+aL7Yhk#U6~?{ogf=1gV<jQe(i-g&%|!XFtSrw$-bB
zQKi*{#;yBa7h8?tk+OZewlb`8G2JlNTQC#&G30Zpw<^*}_Q~eWcxrPOo&oTGH$dfp
z{8^4v9>gEEa46A*s=;-u%qP|IKXDubQ5q3q5Ocs>ndCdH#-w`m*{PL!DhD%vEwk>X
zH*emQ!)0<E1)&ub$DH0h2BxM#khmWr$xlyyIrFG-z(`5<PrRra-K}cM$;x_p!`|B#
z^>MqWf4$jqyxQPCWNBy-`wtwrQm0BiIzfMOa?3N^cz2|i^y%r*ZX&J#1o<Ks#h>YX
zkt3A)Ex*~U93X%dG8i~fv2}kq4(Z)DzTJ6eJN9bDmH>$8M0a={8ZVC+pQL+75b&}R
zjshD{Ng80?#$fb#WhWifUGBf)(4j-wr`r`<oj455%r@f@1%PuakdTj8r&b^tm=*s9
z!D0Kl<J`bcoMx-rs0FcS0>SM*e0Uwkb$C!mnwUvoZ9I2w+tbtSK|ZwfwV_6VNIu2q
zJ?q1moX(wl1@b=Iv2`2q{XtWb_n6EWGMTW#GAFb-DHSoTNWOGQaA;8tx9<JF;VKy<
zw|i?S;*~Q9!*~)?QPOe@e|feBdNHR+gtkkD$2rixW8dzcVFc5OidOKC<ECMx3cqAD
zAGc<KBd7we169NtATBGEA#oWQx-J=9eMoVshRKZj`ug>G_J&(gu_8zuuXbOLd@cZc
zUj!mgFx)SMAVRHQ0l&j;Tic`|6Ee`44uOz<2tb72V#|)jX=_5Gw$=|e!NySuup-l0
zir51E=tld5eu5OUt_L}w9)xo@INh5!Z&m<lYM{7ZO89a~Qda@LJ^n0Kj>Sx{6X9hD
zvEr<F9|C3;grlo)R~*`KRT8CBgEJMzm{MQ8<W(^D7;jh+em~%QzRhqMJSlA4+`+`|
zB_Ke9%)0FCT!SrS8zNXeb!@$Y&?@cN!ip6TaQ(UlPDLa`pf1%9V?=X}-4Oz73G843
zNX4kre88&IJy0RDgdkxbAD@ZQFQ1t=O{sb>dg{hVF9#j_3QYXM<x;V5^{&e1Q^!t~
z=ZO6O!k5@(dNc2f>+}U+S$JT1I09n8|H+twEDx#uA7hFZIARoPq{{hkj43F}NNxnb
zEhsE}EjDsm$?g$3xkr_ceLA1BNk=-!q%M8$*;Q_G;p!oa3etf;wU?-7R?wx&NY-p6
zb`^@0eS<Ki>w|C_OmrdLPlPv|wPSCGqk)P1+sTf-DTKCB^Q_5DE_vxdhiV^tT6szB
zeLOK5GQDLV<(c>1^ihvUQADopeIk4ZhotqCeH3Nh+pq+;VUFm?>c7E$xxYyLDWsfo
zOKsDNN7^vriSgV<5{miS56$~Z`AtU}e-R>>Gj%-?KJx03b||cu+nD#}?8YH*1#alU
zpR4d^KUza}{Fxqq-h)5mXb*Dq6<WjNJhyGa<fUua&RV4Pzqee+&tGS$?c%E*F@#6a
z@luk6G$x|NZcx(uK^UWxx`cV}Ibog{n=E;$)m4_-p#eA&=7%ftk~QH@b~j79%ZwQO
zH?aAOw60;+U3S)leIW=(ETgwo!E>gC;@>p*H|{Hq7sn_vzjc3QsqO5FGrp0R`a$Er
zp`KjnNTturWqHX!FIDD!c$jFG%jGxw7jgc~?y_S2NqNb=d-sk3QyeNBMqr{`hWXLB
z5SE9BSRHA{&o2#aQN`KW@3X4CX0|*z<YcmBr<2oG6h>!(l#9GZ&ey5<gBnY@QR^gs
zj*23H*qo(cRbk%OkhXMV>vG4C66U>>(O<t_0Ws5-EdJcc_rQDc;%asUPOYEbpu}Kw
z<m<Tjo%ihnl_M@Zaa>Loi_DYCl9{BggG#ShTHW5x3l?|3Sh%&*S4<%7AJaVFRL)(4
z8AD%nL7LCKUjSysPtzALIZ+Y3YHvT*XA}Cee(wAEn1thPQX^7WaYtdc_&glCO_GKf
z4Dgh`!NGe&%(EbX2e1hl>GjFCP)kn5m|vP$ec-?WBRx%K{Tr-m0Z>SmBQLP6T`T%N
zR@MCC96*c)6foas&z_-Hjf&o{GB!50WM>HV<j0h7pvN1$^^dcuUl12IQUXv$8W&!I
z=>4j-_3-%cmqQ1(o!gI;Lmh$9kOK-c+!wWv1_H_tFbKT0&4shsDWCr2)tBoRC!dwp
zgr1mcp=Nn0NoSTHNN!wCVV)S-pj}Nkd0$^&JY`8Ot)Yd*Lr|fdu$uJ1_m(4+c6xq#
z>=nE-pRc-{JC`$L$LzB}b@<s<NJ^+?Gmm|1SJHKV{c1EbV(;jfG354DzGr5KE@ivz
z+cJ87PfA3M(Z#!<K2gs<v3<uyK!TN`0{qROc<aQU!Gg4M&E}(*^^zlQ)jB<)rd5UU
zhf(LhtD=6ue<Pp$a8@J^itg@i_Qvq5U2{@r+6sK2x*S?}yKH#PA(AeEtun`acTPSA
zO^U+UxGxqC807W78!Vr;J)xRk`u$E0wtf42yk1n?h>|yNUFQiYXL$?hta+AwF(`-;
zFOvpwm{~y@ffg}|(&e&da?<u<rZx4Mf8WDhSD|sFGwtni*}^UBScbf^Awfv1_d~MH
zo7!5^1VL<6??+%b2OQP_(*<+UqXUXFVBsb7FXW%;e88e3etvwhbaqf;CHztlfG?=%
z!q_tp^|cBf9AjNA=GWQ!H7Ejj-~!!yzcar12h!aY<a-pJ*@G||t0(K-r%+(<yUkri
zqu~Ch-iH+kLf>7NeTb&Soc*=0FEu)pCp_u#xb^==%g5BE3=9I+OHoy$GN?j(skNiE
z5?Oz2q8P8Ab0$yYk`{~LDk8WJ#RRl@*OCeSrfg~&pZB~IYV$5IAHW67u!OT96~iTN
zG-n2<Kg<Xg;JT(RDa-JCzn0Zj*l+|^R%k0cs$M+EP4Z4XOQFdGLm%m!P>gYfY^{#@
zP1vI+Po$6qh-a8#`SQXa!9svpeh_yLa!12T<ktVlpTv3HXek=_{--9C9PV(=y4nm;
z7YT2I=~aaafK1I(XSr6=K`CmRP(9tM+Mc}Bm<|q$|9)3WaXf0R|EFw?#((k6K%IJ8
zo6h9(g}R~j#U&3^^RKM3cP2N!gclf6@>*3dWs!X`*Wg7JKNjr3mq*3OD5E=M;!|n7
zC3QjqdJ_zm<=)gebIRp<yQOD6EcB4-;SMo~9BKaC%hH;$boZ8{H}2j0LkOUkW-}!K
zRZk#{dWSPi^L6mtALqs}Jv7|ipGPaXJ_eXvjSR}6cwY}jynQ4{V)pN15Fo%*qLbn0
z*RR6lG@M(^BtbV0szm6T4tE+n@jl-7^=klG@dX713<xO{1^_(+6BE&E7F45=L8^`_
z6c($oo=QQu+bOK@V{nd?Og)hIns>XzKlr%z&6`R{#1f0YqDlh*xmT2!Z2_7VUb1A#
zQc69Hro>1L+-F>u#Brn_E+n%FNe4^nAegqBFe5wkH9da-ReVxkC#BTq@PPw@-|Aww
zfM6lsu(h&|(R=pnNy*};mXYgVI5i;OEThz?8B2^o@MQe!<#`ghI5+WQBf)WPg}K0S
zv~!07g~nBpU0{+G!{3R^XfAnnvWvFSQ8pO$1cSM`d3Rsm2BKtN)X@kgRyd;30@_L7
zc(?I-t;(rR?=4E7F6HJg`m*+bbsjPo^cYcv?E|FSo6v(29)q1Pm0QY`o%?SPMH#$L
z&qT9%^E_?k8QvY9CgVTbkH{-11fq+^7+o_pq2onQybD`Oj}FBexp>pU=@x(*at{J6
zO0o@zDL<5=9sUPyps3)5H<gRpH#=3OvW^vr%ubDL0W(Dad;7hsrKnxqzTG)NAo%%j
zZ^L<$4-T1cp4_r#ZeyTVXW83#?^dHlM{xc6vdwixNB;X%)8G3NE}HVWI6H1%kJ!h?
z$F~`&EXzi80H{_FKJrS80UA_%hvd~5Cbxh$3SNkyYx`36wytgs;QDGjO-i})))H^3
z7XrJ0nft(of#gdqs&j{M?B2Z#hKX+9Q5U^U`C+<~5TvklQz%8QBQDV^n4g_)=R&FM
z_Wiyq5ekc3&xk4Dmr*_oZxSOHQ-jXCwo4)B8g~y5Z-SMeDY`2r*LmOo;(#zINcyIS
zVi(?2EcV@OS2Z@65^Dv`X9^;bsQCFy>>J3PBm5Yv@>VRKv8JBO_AazRocp*#3I*wI
zM8FCt-4AWo?}NoM#;zeD5K}YFf$&RKe-9qlmr1*1@g5>AQF_rqBSUNd365DXxOW|r
z{aU<@{|tH5=LN$M2fk?A+{1YeO**=80M^HEfh9`$*J~jn^^zc91xa>leCqZV^Uavi
zBBU_ciV}RRa*erYjA6hRYULNqZ{tM?{mjKvfhv=AbKdXVs9AdN7E=$(u)Iadv7FL_
zmX#ph85sWAK~Qi6ui3)6N<=qh-tGZ}%1<*>9p2PS03#+Bca(zIga<MRv8FTVc>e5n
zRjL?Z83z@3cJ|rs$&acd)coFsKVFFZw2Ky1m6d^TjflH^1UKPOx6g#S!5%;+(|e1Y
z+#yU>G$9R%Kl4WO7n~pJAc2g-iw}aLAix2)Cg#qQIG45e7_WHyNPEj(YB9ArWa<wC
z(InoHnk09ZSno;m+|WBo!u%+?o{s#2`5cCL*@cvuzgPnFK#aJczcSBxV3H64O49tr
zQC#Q0G7#PRSt|QKrqCoj24bhE&&L)C?Uwv`z}utzKQ$s{UqwaPE7Nfr&4#FTO{oti
zLL1P*A^wbh=~4;>WWhoG<tL~9i`z(MHY+4HmY<Ak@pMlx7W-O7XHg2wsw|+((2YV%
z&}!59_8Pv~?+c|acy6$5si&LhZ2xU=Iv}&h!9g4?h3FPL{A22m#$8-elHPW@WLE3o
z!Q1C(sOD%d2)}tVuB9%jsk8~^d(tfl^G;c1<uVNJZU6%)7h?!B<>1eCHcNO6;~~Fb
zc7T>MY<FmrqA6=`K9tXjUeJ-^-^>RyOng8Wn3mRoXrrKoGp+6PLT<gC(dMtz-xCwN
zEiE5eshVq2OJyz$8o_=5Mt&viTk62<*_Y;AFxU6Ku{j>L@OzOe(Lw!@<V-bs=1^yO
z8(mhH(_BNRIy_kLnkvDi|CoP)@{Upl5g@(oKkQH2Xhu*C211z@o&0@#+&tpH$NA+{
zFJ7clSC3FH7(Vw-JCg5XuhMS45VhYP0$sH_i^??sv(5<*b#--`d3T4}!Jf)s%4Jx6
zW?81WfpCM&D5k@UC(nTy1s%cy-)}j*gHkZjcuuvuHvDR<=^$@J{6FrWw$WYuEAc|{
zE2DxZB-q_EFAn!jD#d=3z56bb;8f~5I2&JqYD@Q|{wE6+_u6M32EUxcV-}hdQc@V8
zax4YKM`$unn6anlI;o*8Ehc}=iOw562JvbJ|9v%SOmCTa70ceeV*(U47%%I*G{uX4
zQFW+fl4iwgv4E@K;vok7p_Xh~3aN}p0}Yl4<utYVJF(&=?v09rhn7N^FBCK!pn}pT
z2X}4Y+2Oes0T=eR{MjsRR9#q3ECNe=dVA}^b5hE{I$i_5>S;X@8@YJX<<Gi^!Mc!3
zoJ)Rg>`@<v4i@-9h<FHEPW;>lnEi?U1~Iw<51;0_FS)TBR6d0Q<AvMBi@UO@GT*q*
z6#F3E_rlX9L8N1_Tb8i15T=}?#(kr}QWOKDb>54H+14rb-lxEjjrPZfS;6|#q&MdY
z_A`>|{n}~^W|p~92lZsKCe?Y3%(3?t#D5E6;+4&|)bk5Vxruy00gtkBcUbx2$*RRZ
z#AC3TO=;gG7C^-Vg028U*ni&TWCwx=(NEz*K@0-fP7+|$-ps{R7a8WaFnq&v%LveA
zs1#;<U)7rr#Lz2vQ#n*Q=dwaccI`Wm>z$^4s^~j}+NBh;{wZn*pNc67oe)4$&%xKV
znxwop;607vKcl^lK>*Fb>xtv`@ZmeAKGFm5C6WRz0}_8zZw~bf$@S)&4=)(XTwr;N
zoI#S;B1on5)qBoE0?+`%TtU=T$XrG$vU4Y>IV3MlA}^^(`LdEuHWF0HielSih`MjB
zvfC5l69%Vr^QO;2uJ>VF!L&a%tzfQ<Lyh?C7zsD&F?BkT*1opDsh4D@mQwH1k31wj
z!Y|vPLp2)cpi63X<lIs`dWA^GEx9Mc8XhH}57`QoI=O%7dAnVBCbz8nli@Xy;W^Y8
z;_F4B8EwXNuY7*t$E}qcjA5NxyVdI2StqCZj}b%vA`Yg4DIc+lVQRx&<%>BV2Aq3H
z1q;Jcs>J~)t-;lN9=e!9fy9jGlxF!SxwfckPtL)eg%?<f<t~h-W&sFjhMkWRS_+)v
zM<*va34YnW-O%NpB1Y&Bgc2?+B^9pf;&eRm%=s1-6%i@q)Fe|rAmXaJ9dj9!^?yJq
z<Sh8C_5Z$%NqzP<Y8gk*<ZnhmCAJ3-9>gJlSxpQ!$)11t^?&5%A~OkKr<Tvr!W6FZ
zK!f<u(i2MmdI2<-1M|XrT7j4%<8dJrZ+a`jDS@77zZFM`f8uwY?|)#w1FE!*gW&e9
zAY-8d-widHygG*-&*FMJ4lwkqKw_|>6nfR!X-?t_X<E&RZ%S91o5E0lG(ABT?^yuC
z(`QWY^s&DH8`D;V3j(%jb}pr-4|eUmbO`csx=U?BL4gUhd*ka|GCH^p1HA0Blfx3S
zjxFl}sEupZ?=>@f0K0c4W&+Jz)a5Df!o1Jki2jv>P3glgyp=e!R;e*0?pR*=|9+nx
z(<(L6w@zC3u;TWcn$|9E=*DAhn&`|tB1gDT^f0_c!L$6n-XZ<c4XL8L>V0O8{C`e-
zc76Yz{3$DWIjS5!($}5!VwMrp$uTmk$(4$wrv@2+u2id-`sTv3#fAC#4*+YGv-)~^
znt(}TkTB4KyY8mxL;8n5%jPrx7nJG$te5Y<qD+@UcFFPL9e4p6fjGuZ^6c0aXPZkq
zCg%7Ps0|GVki8KQv#E0xbEQXJC&8F4Cmoikv0nl;as2QGGw+7(zijX|XlzwW?7TG?
z*amONMj{C#lZog?6H#l=SnmQ7;7bzVp_=tixL1a0B_y1AF++!WNjh<2r0wJSvq^o6
zj3kH#CgQ*0P~mJ!Cp=#NOJcbpZUZP|!LZGzpL>gb`bWrI*mB}fpPPn`4#w*2*(W{a
z?38cs?+3#HZu#NiE_mb(=iQ9FA=n><|Bm<~@=z6^W#mdkgd#P{&(H5wbu~Rqp=mr}
zd*BiwHjyatYRJ2<u&)CRDmu~~DtEeVe_z3P7f#cH2bY393Z#VvZ0Fb(WB+iA;(g{Y
zmC3dABCW-+a_}fRGk_sV<IUZ5m2P_WM;7U$CIGK%yY*8a(8<dM2xYP!iW^pD#4iWn
zjy-^|sa{9pX2;PI(hH?Fjaips*O<gaoRJVu$gT~B6D`i<Nj{Q&5>TmHA@h*t?VFDd
zWtkN?53nP^3H42DOb`_Sy_-9cj*L`Ypx58td@Kl+CqMdYeFv+PJN;=J=`xQNvJ>Xx
zMC>I$tQDO%^;jY3ff~(|Y|Vo(+w`C`gV>zOE5w3XN6*4{>-N>3q1dm1Y=D6oEKa%>
zU;)s62PN+cd?`*AH{a!$NTTV^c)@w-%vfYDd0fAH$`{$;XM6>^Hu;W(FZ_h(aA580
z*96!lvgp1z_emDgO#m>nRoP`4(m$MLoGyKS6AsyPPw@g{)E4vGyuSUed-e4E;L*ry
zrfCbzy4SRWS(h6YEk8&)I~jbSp`0B2bs|i4Q65&_g5x*$VM<i>hN&`>-5(4nqmj9$
z)jp|CAk;{kHt^rr@Uzaea~p1U!W*>$V23B#CIb9xLTu>5rWfZ97P+zQ5Ao;?k$VEv
zO4@dLm0WJ|CTn`~ML43zniC>b!V9&hvd~CDT*#Yw%^c4=*&Nz@uZiJwG#BdIZ3L2-
zyzo@8va+(NPBq4hAiH8<U11x8pNApSA@dP<(Q4pOap=7R3w@F%i11Y~u9X8C5zBU~
z)8I#RRiOj;h5vO&9sC4|l8$uK|Hn`0?t{XDa#Rftb_^SW)P|lvgt5@V%Xhn@oFj(0
z9iz<au@O&ngtYMi>ADT3CXc=r3bPn!NX+SN%f(U}yQrJ94%4U-o2AH1U5SmsaD9I1
zlPW*}8(the6gQH?4bSiw-o<Ny*7PVYPR)fQu1H5)n}$L*VmNoxqa2Sxi^l*0C5`UL
zn&m|CKXc{`9L!Pp9av^skgt(=osVxuM8I=M2RmD<!OwMB3KIPIB~hoqYE!E_?pN<N
z9--2D-#^^qO4rIpvwD|?0f&K;-0sBPZ;vuc<Y!2}^Zl7WqW`<&$kUUzCG3ATOK)Ow
zHQ-pbd;RW@ExTHG-LO1i+fk}otecu5x&6Ft*Kbc{&vCnGh1h1#u6t86_kJsU??RP<
zCETT?mwwURo$!Z4$|CWr@WPyX&JM-5yglrdeI<1ncG9>yQcC3`$lyJXpgV`hL+N-4
zhm$7(B~v1nMZCM$w?-mbIr#eZ3bgi&qUa`8OV4Wis3dxw+&iAvWr?F-UHz&ar|xAO
zyb1^F=;~HrAxmNZ4avhf*GA-rWu5!?Lrg%!U*w4B%HJBK284S9Fn2n%!PfDcq^ilw
zui9kK6I6wLL(f9X>gvK&_Q%w_KX+gJW#ITN<6HKX8roE&AXH7!ckfywcd!O3R&Bae
zX|)ynSp}m=eE!#S3Sf@g1Q<oiY~MR@NknsrMAWcyagj!I9u;?vZ{NOoY8^<09U=Qp
z&h<~=F1UU|6rFILwr7$(fAGA?%a2b;D8bE5qh!1=%TpUcc<OcQ)}jCPyAor$@;3nJ
zCwX}qaF|Vp#%&|__2f$GF1RWc;6%Pgdhxo-d_1cSg3(i{F3c`AlJH9EbHc!U;?57*
zvX^kYt#d{x*CrS2=rNDJlM3`=j?VwC`^~F$N>^7G-8FCCTW>Sf=H_3UDUPFO@rr0L
zXAld^L7zOac)8?7ZwhY_n^JrT*$-j={tuQb$m6DQOGG=M&h*73x2_0aBe8iVG@cEK
z#X#?|#oAL&-LICdUAy)ny5l^x3?%TpmrPz7#n;yApslv7zMc&c|LT6xuSz%nRJG>{
z$$kf(T9#`Q(e}2IOL|Z8gC!pi;+EX;z#<_bA#MYSQTtJDqtSv%L`1|>t6}YtBS)M^
zIu^3ab6_k2o~F}Jjeh?%LVw3g)L+B~*>{q+&u`hLm$+Yg-0w~e@EMyQ!*zn?1MY8(
zSG(7bnKuaLT`SN}ZfLdW8X8)urKN=?+q)xnx*N}eZ6oc2WDgi&_T^|yyK(2va$NUy
z-?<mP_Dk<7q`UKEn>nr$TvF}yZqbqKRRmv!*TTdhtUm4Xn(6wMmH_mTR#a4U&CmQM
zy(Z!k5(4ajWllU?+}r{jj3nU~qf+;*KCV~*rO6W8IqVnZpzvKUkPp3jbhH1XM*?7L
zrO>KD%vBicc|YS^4AONz6)v(xz!NaIFM9iP40(GWI4YBlaN-!jC@YW8A{&vfiQ(>k
zk_q0|BQ$QI!{Mrlq~Foahfb&1Y`gR85(b}{nc1__Mv)EA9#CoZvALP{@~;s#w^6q&
z2#93E47Z7<LrR`69b**c;5LXu;6Zg(Wj&IbmUbmP+`zwRhjEk1NS=?8_--_!o1yPC
z<H$3;nOX2%L|_LL>ACOP(Y^JN?cpOwtkE>a%COYV(b3bccU%+~WfHMgl%S{gPhR9*
za=EbZV)C8B#k@PYsIN&ID@+*|b&}Ns?DNrP5WZhjn6mUTp?YSac-B-^`N0vQcCPT=
zNB~!JOUoWK>k6|6-dbpAD+PV|i#Kv$DZIY=Jk=^=BR6{GKbLsOb#9f>*&Njfjq`O!
z#{;MzfdT!9<~FBJ?W6RgtG!R*(I=&v_m)q0Vd0a0xVEdS8GWB*XZz&lZpYnN+S$E6
zujT~-iU}Rx-aQG7K07bZ&dxIbcIMRdsI00I#AX$(v4NFV#K;W^so|{nqITi|l&flz
zoW0ngqYT?#tdDvl5cS5=!h-ZIM90MZ;x#;S-^$<S$N)4HxXVAQ`+FeBKo2-eUA%;N
z_A5Hn92y~9bKR>}#6?D~mbiUgQ+fD9>CW|KpvtAAH9c?|LWb|zySE%~4YMz$?d2&u
zJK{)rbzW__iqTf$R7deTcX#)&k`ByGm>r073UEzXxRIHlxvN*zbLH$UzHeE#FZ<cE
zGcKG(I6LYMv@Tbvi!T}<I%EyGAtg1H=5qCB8SHu1C3O3}#j@Ph2;zS>7Lk^go;!a>
z*AE-wkiE2^HC^*-w5NwHiZs`S?JMY#X9io=C$_q=6_d@x)O5#q%pPH~hlpb3!o9l@
z*l-Npk_e%Bt;-D<eU}hauc=ALGT=-r5uJ(MiF{6UY`xkm|MAnO9*7$r-y};y9k7o8
zG4NMdW&e^1-?J1*h0k#7{de(EQfzGS@O-vEZ7|~s_~VByiB7i_jD1Kn3u6z|Xem6O
zOcyc0c~I&I*;7it`0%bq-+Bx(6q`!C2Bis^eF!<55G`*lR66k}3BR%X*coOQpdkpq
zdOVt!l|_2OuuB)2@~rjhN;zQ;j*etAh9Zd%tDHMFxl1G_B?&^1YF+nGFq}bFR(2Yl
zvhX^E0n(1$*k5FF<cO+sNSdyQ^|FY1QQSAiguz!0@u?e->zfk$8AKl9Tmh(i1@^G;
zYUe%Jx?T2tTlM5hs~Tweuy>;z!IdAHo2yZylm2<=MZ<yK*Y4L7?1ODtaI>Nuj4fZj
z{N@cD(BR?WSszjJaKGpi$O3vEh0UD^pWUdATb`ccmXwr?;l-v@R^(m9uq4wu&noH`
z9t{Q5B3lWir7<2z<tV%QMOj8xb_oSXBan?xOjuf3EeX`Xe;d*wle4qgCMPEmW@tyd
zp6|tv#SD|kW+%20haU*a$*n<IpOTqbiK^Pd($WVh0XNt8d)vEr?*v#GENpG*PM$m|
zFE9Tjgu5{#Y2ui8q&U?pEj>MQ)mV(k>{>3a9y}Il60)s(5Q0Y-yM6dbWK@(jX{EyM
z6ze${BfscUq@<*FVpgX=$PtcuV_;--4b>7Vjh&mL)Wsv!L?T|ke$9XtD<v%rN^05A
z@UZ3a<B6eJiF(O-FDyVUt^wMEH9!md`JmSbVLuRy^f=AQ71q|)WCu6QV8Uwrij9g|
zgC<gbu6>IzzyWa^hw>q}1@n7N)1Z^;^wiW3JJQ0uW{p2e0c^ZLN)jwjLp|#=VsLFh
z*CRGd1T8QrRYdeYc``le>vdT&L_qaUKR#{hQ;dbR9WNcPV8?}5wT)Oauin3BhUXvB
zmXflvGTDNT(5skU;)|TmH=t}GTO(n6s+}6^1sRquAuR&fFeI#`^yyPhWO%Y!9+rU*
z?g^eD?1)oftbb*zl+@*jh*iWWixDWt23$~`;rVFzBy{1ku~d3m+F7J^vVj{KA3*F8
z4+({AGJ=h^3^`?y7h67mepOXP#~!%U*49=^R+et`;ab^a5{J>E+x2pHA?6O`%G#2T
zsJlrQA?h}QcHX@+-ti91_!?Y}Xn|z5Bg5?xX7@ze8T#;})9BElUTEF&mu96Aex2qf
zPo3qV2b+cl`FPmG!E(I>k%P5B!_Cc&7<i>u-o(U&LV*OkQzIoK1KSWRCk<!_Eii!f
zG2TIcgP<BQ?nClM)Utx#qP9~Z{w}sWVTA<~8*vfsL?O7btc-$vEXXiJV}k{JB`$!c
zBmh?ml1oQ;lKrKPWM)w)gM)+UPPhzcOABor-ggZHgL?Jd1@`_$WIqNe3ksA|(sPN{
zz=w_R5{yMTNCPD)DF+7~UfXOV?SI6Wq4uZ@JBm4t_Z!?(a%INbO3BLd$HtI3A0~!c
zS7B|PJ$v?YX67dBc(NWV#=(Ixe9yI*7%t@P_{Wc#;iXsi^4dWfJ-)x^CH>v-zZG2g
z!O6wNWnpc76^PUl;S$_~4|bUAY;U~Qc|}_T9$Q9IRH2gXK`^*cV*|Hg#bp{=FfR1)
zWy1Kf29VaFt`TViw2*l4g+^~~oYCGsaE+N84H`Wvt}S0>VPlgpsr)$&_d(ChQ?`h<
z>%os=$O#r85y<X&s1va1-i_+hS;l~qkIzTr!3wZ%0n|*gzsb?UbIXvO+<q2*mjX?v
zduQ7+EYslyhwCv-O-<bpUF6-zMNyw$Me171&25c37fscRi<pBar}R>fZuWNQme{mu
zx0>3L$LLl^)Fh8Xwx+>4Cml`jA$Fs#7&6k5`Pp9NkCC=U`M+Ke+hoxK9__}^;Z^^6
zPe|UV2M=8oh^VodMl`<l!_C=Q8VIf)-NoeAvPveoI`x$lE8o0pwzjge1LOU5Xh@fW
zDIJ0dsf8Svh`GK)P^JFP_7irf;_^Ia&T%QbNiFwPLtKeNm%-;E*Dc7Zzj+O;!OGhq
zxe^>=l^u;#EG8&;B`7#p=0RdXK{QkVQF;033ahQ!!t7+z!Nla`lgcZ8I|Cj&e|=3h
zTHyr$7QS=+Lu`OGiyIz*8O+Gsyv|i45pS#k%%c`exypB<qi^gNO@N%x0ErXJ<NAo3
z!^TLR0e^n4^-vL(?0ow4={Q=#J|F<$+Am%_faIo$Q4$K&P!kLfAjwYg;;)3N-1$1~
zkkxNN)j%0oH8VXe8vJ!`q*NcX_=C~XNYPE35<$5TptVPmGvdyjb*QJqVE&vAO_IgR
z$KE8$NWx$H`}cq~s7>azyNcaXO|h|K;q&J=@f>idHsAy^p-lra3uY*xbSNOAg3;!Y
zxm`aaH8V32@24?6>xE8IQc^-zmOg}zMEq}zcVAz^0kO=(xqd@KL$AxqE+^^T*Fjqk
zgdx!v@oe>1du9M08v$YidgTWjlUWu}7A5Q#m0qX5uxz`K(A95m@2Jz!(V-_*6Lz=}
zsmjEZlz<)7`SX1pCj0lZ@~V2>8!-_#Ja9k_+?xpACLq!pG$b5J=Rbxzki6=TXsBgp
z^p!^Mi1?;WeM$$DBVuB72;Bo<n3J9T2=3b2w{O2H9W2_Dc;GX3WSW_sWv6h-pIyIz
z%_BUk9pnm+>S0g3jYx=aAjE@_G5Mj@0FOEcJG&ZBbv-1F0Zi@~2h&T-wr|=Lg1j$7
z#<RP77v6VLTAKI;nQ94NwOz^ju}Ro>iIwteYAP9|%_wjBI}XNWr&@D-Gcq%4>g&S+
zHxrwDe;I7Ma~*r%>_S06W+5VKOvkF>Te!NIq~vuCK{gck4KTW=K7JgE=Xb?y^_OjV
z{ygH|jT<|$3WYg??Q3psFf|wDkBpAqBcyqKz0b~_@W;v`MW+{OH?4V#gb(g?9k#m(
zPfp$lnkE93ka7F<PeEZ}$=G)AfqZB@LUIG9%1r1-aU={#Z9w!IdDy=%EhFQr(t*3}
zrQWgFj&u*^qKN*3x0-hy4jehMm+UaLj*II+7@rCmf8mIl>Ci+Ma50J4){TwgzBuz;
zOiAhP$evAc`$hZl{x@cu=)xB53}^BSq5L~)j~>E6bMyXv-V4K@BhliP^8TCA5veIu
z4`e)!AK%@Hsv}waR#;GAAlk-M*@G+d;B&a&()a=qOhHlc0Cw0D6A|$nkrzLYRC&(H
zNdRqUK_MZ1N_zh6yuADP`1owBtkggsE|Za#9z_L;P^Jz<ge<IK?lqHmkpEqQ>j)nV
zilJEOmtZ?Ey`L`wJ<1LyON>H4y2SP$+MXsBeIWH#hv)BoOLjskvGn*w6~e4T3pSat
zz^}J%-hAbnuu;*upA|#W07Jq6ayqDj!knCKG3af;@~g*tsdDAHzcK2KJxW($4n{G!
ze_=B$cNI2kYHn^sf5REn{qflSFBBkVI&?f5^2x8sNfF@GeGjG@*J5s&doRdf{>NaK
zpZ}w|Gmp!0ecSzGv9KbfSVD%(nUW->45bVWQW+vLHVKhZp)5lg3Plk`QV~fs5K4wJ
zSDL64Aq|qsP&DlCwXF5~yzhQLd+&eu`uv6!J@@lG_kCUGb)Ls@oX4r*yD=`ssyKPu
z30R_O-S=W2@^yjvuF4IBy5R!`cy^YSwuM9An(itnznn(>hzYAy{*DQ|rVn%|Ne-@2
z*SK)z%wIU2%pvN>u_(KYpKY=m=EoEb%7fI@X&3g*!=i!0@bJNd2cJEex9&nja4S*N
zy#M>P%Tf>UT-Yp99Xp;TT{7QYd5E=T3pjVru`d?7t)-g2&k>i-94_~qT)ENr%s50{
zGJACuHSo$PF>x{iK*#?xkQ}+34ZSrzOZXJ2F2B=GW!$(o^+D6Z>3*eN3_~Jlw_?SC
z`RX1mlf|6`j}LvMkd*M}6SRVUJ9f0Wh53}3scBf$RE1*dh*`X&d(EwM7dgWDptfz>
za=WxkOjOt5X|c-2=A3z$!~sB7pOh{K?;C19KI(jX8=KQ$0tn;##Y9TTfgk4Q2a;``
zKYMmpV&APWuY*~^m>|y0RV<vZJa$>KL0VdxDjJY3+d3Kx7nrHjLl0>;zLPChDPqgu
z9&K}H*B=_2A9_cTZ_s?!*LBO5@-GHfhKFZ|BYDtgCEhH0IqT06a#BGE3&mICfOLKM
z8Ro<vBjsmLUYp0WI7M(KQ;1sO7(MnFht8^%D|N^KKIUAtDzu0C=6;i!`viVlMo2;u
zl>O-8A3^ix&YfzR<P-|+P_HoK0`DR9#~-1zmy|b2cT|mzjXexRF3z6bDVlloa_oJp
zu0eB@USX~p6e9`a*?o7jqDLO1)`!r0ke(iSw|Li7e-w&(eC^&Tt*Bu~IrdE-YAcfk
zAjZp!^(0Bc8=9q^yUk~@m-(Vavp~(x(wu64@#&8rZW~=(A~~7uE*d&)t|`_MQ5mcv
z!Nr8P5eu6rmDS&)DXVyoyxV8}>Z+<ah>mm1muC&s*^`OnST!#s0d;CGlrFUjIwR?-
z@Z*$2N`0VaxJET*TD-F@j9qWr8bks^KrW+`K8!yl`l3{nm2GN5D+2x$$EV7~iIMap
zCPb#HwVnFg4o~r&NL?c6B2iTQL4&d_^u5R(#8at$c{Oqd(x}H_P$<uF(=edY*j{dp
z6Lg&0S38T#R&Cl0cs*^Rxw*Ls0jYw_S?T9es-eN%I@+|hG%GDG%U+V?v{!PuQ$wr3
z7~6RoO6<}WufH{e^=Gzt@rBxi_~pS1!oj5u_Ek{GZE@<!Nkpsq7YeMXcm>OO%a$4Q
zW3|HSgqv0^Y%Il1n>IzAJ7;6O)1kd51xY_Q2O`S5%F3w8YL-iv8o`f=-Q08Z_=yt(
zC5Iv-KggC~!W*^3wX*8(%>h91JZ*i}O=`utfG!g@Y87Jcz%$mbt}YMK-u5Gt`}fqP
zbM!njGhN80IJ)GQ7Ub{Qy!7a&;#+e_8SVP74pDdfOOnY=$tm!?w1^V%WkEs0hjnJ7
zM;?>f{TD@E#i`aN>&1QS67^dj@QG)Q>X{8|^DXnIaOJ`M*<^o>+Lwz*)3dYdLytz1
zy(y)9UqL3>X710hXN0On3m0Z$C=!<EY*EU?bpjtuUO9XlRmv(`+m|D(J5PDY_4Fvi
zZXOu3H_n=Sh@h53cfsQ;Cj^ps(C4)L8B9b=)qU}H&&r!|QZDY@xNPbvjP!2vD2qPn
z&%f|Ubit8k<1aGh%MZ~`Woo*UWoR<G+RfFKj%}Ur$6xB}cY1ksATf)0{$5*Y(j<A&
zSgh%M(W20H@5JIQZ6P|}XzqCqb}3{^6ne>FHa_svD>+E0j$dm^Nm}w;&4aMbrcRu@
zWXZ(@HJo@(_o%v-Pyt%bQ>3jqyKHsbIjazeCyu~!d_zw3v3o||ynOB2&c1IyW4oj?
zCa|pp3!1juwryKmlH9!al*5nbMn*<mZ+)U%k_ejwiD&^Aelm#2&HWm&x7^JUH9O#|
z3!g`Hikl}FrR?|d*?s-`By?>&y)kvQwL8wAA3<=q=OPz9*V57l?GDY$M&XZh@AJ^f
z<o}@RA)SjKl>&Q7q!I!Qi-4FaG*}2l3(GPR&sjx8Ecj!=$bYo}5oc&Qd#!%ysq_SE
z5<Fp;&_;MS?3p{fokXas2w}Pl7rv04a^lpft`ea4t>v4M151Scxn!!+(NfM(@@2DG
zYc^M3CpqrQma|*r>`ze%;Wu@FO;*e;QyHn<qMIIJS&2#9bZn5QyMWp+JbzzVsu43A
z^jg%*g-IdUVku6iyD0a{R?B>%J`#+?#`-e>H_u<b+>P7qu<o5Lz9qZj?ekZl`w}fo
zo*mpnNjS?)r&t4&NDlU>jc8qOV<V9O_@2IUW#07_b|c;0-ObF+L&0izD?8h;ix)fR
z<>lEe`gWBRW`xHu$=qejx=I?zfBL)iKQ%H*S8*DhZ|ke={lQ2iQ-Ryel0QC6epXme
zp!9P=E-puHJOjTr@zv=Ncdc~i6h|e0CRF#})r5p~`5MJ%Lid;{nHTC6J$fW13_8gV
zkGdbSes)qA)Smh=3C2q@b;gYQIx!q^ZEx%RzI5?|lPTjTO!!?A|EWiQPR?x(aEHep
zL9i`b+}xB}C-}5)X|c$zIWm^Gel9Aibi&SrbC)i)m&9NE9=t7h;)DsidBew7IO@em
z+Y0%6FgoL4%Ly#@(>65+vz_YJ-xzISri_re?unMyH}R-LBO|>~FAhLfhK6cuS^t~(
z#d)T5$yOaPe7Hmc_hP2CJXc7^`lhF*rRh01fU_t~924XW2Gi|u-LmYonwl6HnXxZ}
zif*_ZOq98naCqimNlWk8a~9T1KR`?^2fGPaJyGFO*-EWtnDMQwst;=wd#ZMc$AbCu
zBPQ=E={Z{In6tApRgC?nK^ZF?9DLXl!eK*<Ao!ftXZ_^*@uEZkTHC;;@>+Uu-@d;%
z|JOV@7P>XscEt|0$CR?ViaQpLn!mUYS#h+ccG4yCfgVp<Z3|9Mu{dcKK0c<N^T=*e
zj&#%J;Q`Sfp4!*&kp7xP`@d+ArBr$<t@RkdJ1aBOlY|Z-;gtI3fPetuYqNvXU8!pr
zX!TlvVuDWW#-WxFt^yeh-cjJxFul+GOBCuD|6GU{d?W;B>|Nidhh~QZYl&$Oq*XoO
zW0wOpAe@Ve_IY;g*@t)UUSlN!2;ylGPy$*7IW%Czh*_f_=Wpg|;oWRo|6#Feg*8CT
zxeFI`D6`20GRx9J((>|picVe_K8vQ~h|i`e9)1m%Ag~MEPh-iyp9E_zXzWM$ArsMZ
z9`OKgjmua)NxH6She1oAFj#zK{d>YLgQ0gh>pvrkn-g6vu=t_=12LDN9x}6EoFr^n
z4?TWQGs!#6ee>ZbdTF_-iUAz>qMQ#GU*Z%|kOgV?Q!_h1Uq5Agv+yXwhI`IvhjJyw
z=}`PAI<qeAj+=>j?=55s4!3mfjHkBu_wq`Y?XTx!T@<Wind6#&a45{g(8Duz$z$yt
zyu8nWJA9$!yKqavvYvG??{C5{)Ex2joE5PP9Wv&9JW#cji|*dDCu;ReL)5V=5mv0o
zP}+WvFWXTvWZ1C8%`&5JWOw~+mGKFp_&!guU}^BFQcC*UP<bI5H!f7nYiH>9aFv~<
zw#*yUA9qnXO;CBI!5k4jO<DNZcymij1u;;A*9X@Ow-}<TDn}G*fbhw=()29h*(uN{
zXu3yLb93LrSn5NdqRbD^b9Q=4Y1+}lqxD;`GoDfaTkqb!)j{$+_od;Tx22`$Ff0&7
zDVjO=O3S9&lwloj;*Z!>>-toMppYe(VfahcjiLuVA#oDO1cz08e60(PYoeZno9k_B
z{3i|+JHCARa_;nLPd=PAN>Y(HU-=q+^-?_<AvWO`W78|8`)i*>krjD3!6}tu?m=ee
zA*j-vko$3;T#p%MI2GP>-Q?ce97V%W@WV<@*vPF#W`uW=*etY|XS19Ue_C$bu)(Tp
zn_-IQ&1SP(;~3_^=^0K?s4PjDf|qR6yR<LCn}gT81$-75$@8Tw*fus9K53dq$Zq}7
zY=TzBTaP1+)YS#s5JeveShel7><j15hmxK8CnmbMd+d$Im;K(H&^$f9<secLr4b|E
zXzt6S;8P^EZ>3cU=GUJC8fl9+e*fWF&cGN30E{=a9P}1rf$GN%!wVvjoOSH6uV25a
zx_z;^l9(9z<jIq7?aPmlwd~%t>rj0B1ZHXY2X;oyph!thbbb@Q&0)i*qDf#BaxyZh
z9JAk%%*duY#<6*l!^)44r-%yyn3=3x8KzhG?CH}XK-T1!JIKk0Qh+A8)EHq{c2db`
zNJ6qzcNtNaQ|;*Ud8}t#JxXq8%Q-7*TTYeI36GKfJ^*YGC8m8BZ8G{pw{B_OYN~s@
zRX<A}$b6uDuU?tx1D?6a1q<*J!-2slhmt&6wa3bf-2mtz+H#f_eDGf~CNglT+g|XU
zb6CRm*uQph>R?6U=3ol*EAjE+xU`t>s`w1&ZcHjP(A3lv6Lq-~$8$Gto^<>E$&gMU
zo;&rg5Af(f&L{x`VKe<%LEjpFS0vgLh6T7pmpH~1kW@ane}6ZM@DXFj9?Ty_Ywchn
zPPCiVjiFCYO*z%MbLVmx-YAPzg`I|G?R@izqwnF17srJ2pr}KUXgn+|JV6c`U)fZ{
z{W%XN{H~_P7}t2J@*5nzVrgN{=ff-!xb#inGtFeBgMHuH;ELFL%9C@Q1PVu0O@b4@
zLf;Gq6DH(gPR<bw)Ea!loXK%%ATprjL(F{k(Tv0xtL|FIW+ZN@T|vIN)b)Ljg!C4;
zrD+r-q6Y1)vdUX;ol&L*dzL0CJw8<=VqPeJVS{-6#tnanEkH)(0+9kL66HJnJHd5_
zwr?V*IT9F7@Y_#W+UR>WCo^+u%*{HfimED?^7>4%4p*ji20SpmKlqe>e2SIYB%01%
z78d%4g$)uo+3I||$}50d*;G*={aN5x@owvHUcP(?fwXb=y=Qeo!bDXLU+j-{KFMvZ
zLZpu1Er6i>!o!C|<GD5@Ft|L5#IO>JYf-Dhtu*>(kD^!+i^;LCK3k7ZcH3p}(MQjo
zML;WpWSQoV7@OO};7o@Pdl0zIq1plk)dMCne_y&;52oN2u?eWf2BuMr6agCYHd7Lb
zJa>_N?908y0Tq7W?H1UV@E+8`Yw|;)7Mh#)7;Tw!1lXqW=TB2kgObCl)I`G5N1r#*
zJJ`?Aqor|J$>ER+20RRgNCKxQSJt->^fagMyXtB~p^pRI8GXzUfw=AG*QpfJh9}iy
zvzMG+LiNmfI#gCxwvyQi@a*O9mVReRJYWTYSTiEYuSVU$86+5l0(-yl^p^cu;;0da
zz&{h%{SOIlMd$Atw;MHPOaxV~Uu^8SN$TpEsK^@t+R!J7y{4k7de;2HeA_0A4dxG@
zJOR1|I-@r&B*~>N*HrlgqO<3CF8v|#={w19JCIsrvDBOeme)7#7LX~G*{b#HjVYbu
zpB-E!Dl@<U({&Qr?%jEh6>u4Hefy%vL)3GVQ_TS1v)L%x6j|f|6JNlNtKJ#C^@p_Z
z7h^O)T?Rrc{Xkw4U%6zpJv>4tYQkwGw&Fa&@jT~vp9G9iOt)-N2j}vGjYHr#U-oHK
zA}bEM$S`E`alj=C!Fl`&tev5(;RQn!jLl?(ISl>yZ7K#Z5d?9UYd7ppYHtgIqyae)
z$OIfz@*u4&611c4C-{_=a^ruC7veJ?`SE)|DG};rQSL<?l`XcIu5|P{H_`@Tg+~Ue
zk_k0^;VOu#%cr7H)()wf0E8KgbRutEI!{QzQjfy6z%EA*&A4wDVWq}V-UTtCA3CLu
z5-u^9PfU2IxEErG8(XyHjR0T|uq(gAhh=VkE3mh))HNHzQ%4_iM6&~%It-dNd8G-B
z{@ZA-5bH@IU^i~^BbQNS9_vby4f}8X+-8x5MfY?^xo+L23V)TnlF#p;Y2vUBmy(ti
z5CtDq0L&nDRdWZGr4aJN#GytOitAh!ctNhK_T~FTfO3>_kER6zn2Hj3>&MGB307(*
z!{x9b*M}e*g1qjxsXcr4`~tY1&1$5qw6_=+raEkx7kP#bz@<r<D!?{<XPF>s8c+JI
zu_r^#7L`A-6S!I2Y4p%jvH9u!m;oX>M2tw3Cx>PZR%9at(V!F=yegs-CyiMD3>0|)
z-*WwK6L?iNm)znp&tUWB3A7Z1w3M{|^pNY-OMu4W$_knRNJ5Lwf10pl;li21IG5mS
z*l6$^z{_;&)|+YV)u5)n06Wlc*^6&3lkO<a)@0+V;CJ5M1uYwta}l~u`(<F+jS*+y
zI70J8vh%sVUI!F|=+lONO?6F8o#Ap(am_p)&O$+Z3U=ARKo@IT$}99h$~Kt;kbz6l
zvVDi$)!*IE-FooR*)sC-d%<=}L2xL8rh|bB!V>K1?>DKd`|JF9+Z%->7O5;?$J){v
zC+VG_-Aj4gxM|p2B8>S<DA>@}m7f?_nMb<MV|hd3Czw>|F{$fIx&Tm%nKlSnyx2Xz
z$mH*l+GDQ@g&XVia6rJGDXmSizc^pvuz_j?qtAVUH-$74X)+{&d0xUL9=)LQy7Zgu
z4NL~~-2ui>bYda15`Nb-UG^cPthUo9^Bz4~nrJtbI4S&pnIeG<K<LWF@Ey^eh>2)a
z+5Y%fBcg{N^8liDq|t137GghL#G9PzT_u^>**nR!$N-=<+Y%=@M|#Q2r#^q)U4nzD
zgr&?u&LRJsBq}OOU?~5{81N>gHR-{_hpsj+#gsK#$7s6CBpTqH7wc0yJ|{PKAU!s;
zmdxN0k!k}#Nc+BJjEoe)25rU&4{O#nN*4ZC;z=KOHyOC*BPoHT(s++|Ah92y$bA>B
zVZsNBc6i*xBpkXSDU-sapcM3(K_Cc9dM)VW9)Ja);fbL<3l_^=D0CpC%U7@N1OlN$
z6kWuTsv8DFcbX~vF5&x}NlNMs7%L|4K*7_dE>s_T_OzFHA_^k^iX8PEE~enJNuGtc
z#pI@fT_v1}uC^_i2gP@0(_()koM90->MUE<Ur$dD4T%?pb<1`rPDj!&8e<0{*Z}0P
zv9a+)qeESyj7zbY34ui;NbO{gDleO5saxUnx{IZ)rKP|*3{EhOa}=ms_y=)PzTcCG
z?qKug&6HKp<{cz(^mm}vi8PHF6nMYX6n|6fAPaR__SVyM#n=kBua&zA&x&+pi!Xw|
zEve&9pnW+xIR<`o=018&4+UHJU2>a!-+04s=3sH#P%$7~evLsks6`l$QFJ~X2{YU=
z^>+QlH94?XT|jn)zJo*u`SB6WBzk8bNyjW1qfZfplnH6xa)Ge$X_4TCUU^XPOWHt=
zw7_>s^u>!?S^|>zzs_PO69v$Ti06R7KnP8EM9Inhz+ch0?b^Hd4{`znLql;v#BT51
zKH<Sa<xH!M?awUDS0j!ca4y8uNn~hrl=P<>Hks3h^LW?6gI$HI6!1uVT&nVc$|eFS
z#(KRqW*(n(-#PoZ^B~10`r3pickX;lF;J55Vev5@A}6)GpkR5gJ+@s&9y7S1puR!E
z4ZKZ>n_HZ>ug6GyKSndsYi_l}xXyQ*6vXe?*g3u2<)b3M=S`=-f86-*kEpH;PT0Er
zw;yZAC9S=qcf9jm8-%Y0Cq7U3?>|oZ_s8!1?S)#c|MsKp-SR7c^jWRaIc&s5E0yRo
zqr)B#U3=uA@_~rUe?5%br{*)}#ATDAca?)~OAXClQK0tki+|%<`)!dr{(A8`UOf6t
z^uI4oLa`g?t#<IA7gwUutv4vp{^!N#5DGsTUvcN(7mwk^#^EFX^J3n4^3gH>zH=q-
zT(6t-?>nF44m^>qDE;@vcD(rM^W3b}vTBc3_0r2-n7S$`{PDlnQsRS5*2~uqG<-4K
ze|azS6~Ftt8tC`3R9hZquxGj1?j8nZ0~J>Pb2m+T%vV0{<#hJnYxUq-M>GoleJ3g2
zX|{snKfi1{cD}N_M8V}hKgyJ}=~bPnnp7kbJ(HX@NT66d!el|Fckf3oefLz%CGdTQ
zRMq}sK?)%O`+5_akBFK4D=c<x)Eya_T7|!|P7-~KlDG9CSKB%D^J__4nNWaI;q%$T
z>Z&_A_ZITPJyhjbe^@p)_mKLXRVjq*XpRnk*ws1g1OWr^YkyKwl6CU;-G45<4a=qq
zBxDq;b0t{J805_EjI>=AOnvwC(DKX=WBPynzUesfKK^8euvp?3(A*rY+0r@eDz&S!
zW2vsUZie;-vN$K{!Sl4{D{CIdZxneP6z2+(O`*!gaPqIr%*-yXpXxxlm*RgeB@k0G
zG4UG@o;=yLaFia^eMxOEy~w8wa@%;l>3q#NRn>biq7#nHjDIodiP2E3us2>`GO1QL
zjT`E}cVYMoK+*K!a&_`14%A^-XdgauWEXJ<*^%PB?H)Xbh;8UIK*0}ihu3(F>A9a<
zF2#u>J<L>E?as$lU{n2rhD6OBK2`Uc`}Z+ItR_AcNmWMCoSY$iaav7xu)&1}op1D)
zcfTqx=$?M2DhECZd<kn*;EteX*Oe;=3H0xcS=3#rHN&vO&Hh+rxSRZ5R*ZlGE>b;d
zuq}_FtJpCX(bW(;z0rjnAV@k*Ss@dzGwx63#|-d0Aq<utHcx>&iM=ZITi}%XZT~pb
zs~(y;UfenA0`2;3ke(<!v!*KUsklbnBE46)?nTv#*j{*DuFX0-z~Me01fJX{(sBmb
zuCNG!znQo5yy{!%>Dw4g{>b^>ElXDab?3Wr2hfl5Df(%SDnyQY_3FwK=oa@*wx9Qy
zICL4lN;*92Caq27$4d?xHZ4E)?F4MkO0U(`)z$RKKQL9yMO`ZU5P$B4qt<=17QSDx
z45*VH_WULS3ZW2tt^O%ih;5zw`7bvcfEzi~)T7ov_uogyZtr|6d&osAXGIdYgUzWM
z`w9ryY2S&Ss*7%J8AocFxpvuQnW1;jrQt_*M6+Ds&4y98rPf43<Lw3RdEdIJuNQgY
z8#KGG$Mu*}>{c}zT(mLt%uBw^&F@3*Dp!wZfxiEi1}ihHM&ajhad(CI{A_<j!N3Nk
z&>in5_3&Ul?&S2)8PIL``V~x0NOKvrr=lB3gb!UKZ6qv^r|Wz5$a>XyE7-unVKgX!
zm_PVS9tnBdUt+axU6&nGB|LaEqjytNLtkpzsy>zU=-BQBnr3l$_-yaspvtG&(ANvv
zc-ol2DYbnglzyCAt`SvGU?BvYFeRKy=*c?Bgd$$Ow(X}Hxtk+Amc)}>1DWf9!!!Ah
zQ{}T2A8m(NM1Y0^`H2&btV=lN!XfvxoYLq*!HVo|2Pva)f~Vuj;FPuk3RlFBcX9GR
z8?pZM36TggHe%XP|K;p(#tF`|ZH1j!i4^Q|e7t$~M5|t34?kQbv}xkfU=RP~5Ts8c
zhah1%5jgtrOp)GEOF*&Q#um!8V}4ea$?lRfuZ#$_SJBQ1w<07idL%m3IR7MS9Izta
z6sZM}1B|}73Mdd;xKpQ2f{^&#yCWICcq3?kvH&rlkXo6Mg)>~GBanKt!-akpANPLX
zzLUR8Qd){zX-H#3U$Xd_0%U-qZAEj2O+1z8zc$A>5j$2FOhBDHeQ@{GmoMef4+=gA
zCU&{Cb!6Ys1utIQM!YuZWjcu}?u^2g_R6{W3l@Z3xzb&vbr|N8)*-k!vRj>60R@ai
zLUv1TnwpW(RW$E03c+f{iaQGuIk&u$k~I2qzW@AQla5}B&F+b`SmbA-m+8<z9Rz|=
z&(*prXzW51BqY2ziwm_M(dBc|rVP;-F5)IKUh{&81UHM;#TFK4_-;xbLDN^S9)b?W
z8wr%lj~44VV?#q}5QjbgxH^+Sd0~lnW#Lye<fOTr82Pee4Nb%OEQ8RX{wf)!LZG%Z
z)f-cFi55Un;t+iC4UBm>$J}T~FCV6vhG7fk4@Su@mg>?a0z+G+F=OVt7fS6e4Q(9@
zUcsO67em2B0l<j*)5+iLGVTNhd}|Hc7TzT|wpG4*o~5NBr%5>J5`q)O(77}Da9>EM
zvWe#B5)uMxyt_j6oJIkVNr`!mys(`o#kP>`$;rxw?d>>hQEHAS%MiSyAK>3cC#SPQ
zr1kXad(D7(ERyJ@w}PNHW)0|;Ny=5hDY!L#o_t<w=lVb{W<tF#;3@!FArC;0H+;&J
zQ#Cfx8xa^;Kt}^;%SgbSJBaH5@#f?X9(k<N?0<;JPEdK3#D>fq>YoK^avpX-$OuS%
z&FI5*xu~|XIYZfF5;{oc0GcgYBDzj}`o4M&2iVT7DmqaBLg0}$yx_>;!=aUHXuatG
z?8j6Pp^g#Xid|v*b<dCgDfj9<$>6-w#a8*Uu^dkW8D8<*)YQ~e3W40^tszM29Ez?T
z25drsD5_#s>v_A{i&jw<r2C*t+%p;aCA~=p&~B|EJ&&XVgrxiojY3uTY`nfsxjXrg
zpWk$W$ojr5CsR0$%%M(U)=DANA$Gj0sF($rFW0Z1(Y*6<xAw38-Bc*CX8QX2<}e-4
zFnCv<H`wbT*dKy~34)Ij#}IEW+buUQ?>U?noW3a<hYJ36tA42=hbJ-Echu<7Lghw4
zdJ%G;E+-|r`xQB@gP2uKP(FxY3xZCMbv>9yg&*`YHu%=};=$@1IeF3mb@r-t>x|ew
zdFM@KD$SNGF%Yv;+nffmnMXWb-1Dwfgc4Ife1Mvsf(8S+mpJC=@%-rX5Nab<3YhLa
zLb10%gayb382OLPM*~k(v2KvVHAARJH>;04R=v4Ab0EmP(0o8WJb{|f1-R(aU;BxR
zocLV2%>M!#9zJ>UNp2k%&*gFJ2rLbaCTP3eQ~gD$YY(DiMrd*x3=7H4TLDD^(Fs(u
zSmkodhs0L-79sUPj3lsbG@K=lPI>?)f-<1<h>plJ;r%rquJg&)BUF~>$MhHKaAa@_
zzfJq#HZO>J37i2?RrK%GEpF86ROPY#=c%kNh>}DeASNQd<te|`^ffyC_`R3|?8p~m
znonv;N$aBi%M;V=FVL^xBVc^6(wgBn#SVq`&qOjXp;+zj%e^UqUz*uB-GsghR9f?O
zblMsoQ{Iu}R56`IC&Nusb8ul#A}P5@Gf?4H{Wq_E>3^{Azmit4s=4(~?H{M!EDI^!
znk13*?>$R*ni*y<@_qVbb6*h{1TQc*pCNES+*~Hx&mSkQC)@=Wj!Fah5`-qm8?6qA
z&nG<ZSJmoeo75`5{o#8t9Wj(*k|VOWtJ_TdKwNmy-fl&<C%Sa}i_`km<XZ5_U=?*)
zs~~kbc}6IV1+R$q;Vs*ne&c!jwG**C74B@5=3=6)u*jvc_QL345_-q1K;4Phv*$c(
z8@E5=$Fe}1D?B(MSw+Cw9;=+L!p<tpVN>4@<Qy1i-G+%5I12nsaar}V?1W7~{K9&d
zQ|>bExDts>C>R1P5*`kN@c)asQuEED^}9y{EaYKy3FF6KEi1l52UOlaDx5(MXq&oA
zAY|5N{ag@jF#|=98Nwc2U2;WqKC_`S%&iR+6$OiOtNB}@-Fkih=9oFExviaB0awZy
zS@tH}K>}=Bn`vMD>Qme*%@4e+YvmA=ZQHaE5dNSMyP#oo&;k>aX#z9E5JjtE^vAwy
z?2)QnxqSIFx(B_L4heBqfJ%M~M{WNY*o>UV>%f5nIota`hcw^EF3)wdmeKyPY=O-g
z?$8W<eOU={wlut1KT$ub;L~J#%?}k7=SM-xvFjr_?Bi6LukQBpf|Ypdzs|_x%BH9H
z*GDqrMs%g)^|Dpp^~wK9j1ygr&{RxiCOa}&+^3+Zb!)b6>-Z*30yFx@9dzH-Df}++
zC{$p&C%MbHt=GwtRd@s&ATQT=Jfj&4{_ajuQC~_I-ihq1vxJ%06`}dT!NF}_@AZ&G
zK?9JqWcNH}@6gHNlSi_x_+D(w>ER!$C00oNMVCdoV*_geSgwr(?ca{%WNqesw%p4|
zaJsRH+3r9}Z{ent9yU=dOb6p1^7zM%?ae<@=^<g5rp>nx7p)SQ{MrTZFRo6p5sjn7
zGd%pgPR!y(i`q!&iO4BPo;!ROY!FabPTn51?^8-tvv{z=Qxeqd=op~f|C^kq@Oibb
znbLiB*xW^n&ig&7)Of^RWl{ijNWsU98FseuJ2^9Uk^z*A3(8f!u>b8fuchmMv;u0y
zTLbpeO=MbF*Z8rO^PTZo@BXggb|Jag(_PExRddGY3*=5vu7E$LjaUBs-$gTnGEi;R
zpFd?^;XG{c;C6(vM=LHqg^Uv|tsgU9Tyh(k^S2D;P))C=tnqddp8w=8h^;0`ZlgM_
zmOtR*V`6R|p76f@Wh-)%ohaB$TPZ001J@enylU;760|+B@|v{h%Ii65p(mpCk@`N6
ztgdBmO7hmaIE60!t*i5IXXn{oowAZ<IacaXkB$X@54m6FH8V)HdW1O3tL2eaMdry2
zIMtj2Yu7Eyw3Yd?JN|tSI>FR^bsQaIEdp9m;!QU;?uz&RH1aE^a5X$pjFc`^FdsQ;
z)UN1ghhOhAxckHLT_4-d%YNu85j~t#jI}Rbt%oo?jT#2r{hXFDBRVZ6z9`6REzlEp
zat>EzT#GP#cw8499X$@;)wM;}WP~SD@s_tf*4rZ{6b7fPkBM8feS3>+SNE%by_6!7
zp#I2yu{?%CL-MpWblZ=LNE#+hjOtKP-J8-<+m&mW95+@vD$EudXin&V;(ml39YI7U
zyIb|%5kCoQVsT3|I4fFT1t#+VPhXffzIxSsyk<R*j>+Lf20Q{9RzBWR;B1sF8pjVE
z`g~7wP3){G(cOCW>PX`u#P5N=ljcP`!2fq)z5FggSoIcd-94Iw5n^?cD66Rag-lcT
zOv?7%@a{Wx6;s8eOQ7Ae5fLlkWrSy=#ns<L8xFgQGW1-6n(E8eBt2Pf(0Z7v2KtQY
zqU&ez-TeN46Q~5NbrJn=Yu0>w5k4TJ;F)a!2Z+gn1zu>VFTFf<)%{_DALo3C`=?1b
z5e1-Rc#k244kQ9n;PZLT3tJ}|oAM^Ygpv@IPVf;H7w5J%PM9!(F{q*^>_o6tx?^84
zx)2oaXkxRf7@y0E3;ne5D+Avdy`}5AtzwYCUZUmbzU4Y$^l<QAqCJQRkxDn-rEduf
zTc;TLi@ihcilv^C_Ro*w<K+&D)-_~YTPoK&Y3FGjH;DMAhMB{LhR@4sd7RQ{sFWhV
ztb-6M)1yMG)*XZcVjPGNiqpLjpLqXAt7Fhdk3AmwG0TP_HQ?&z59@HY_%AY$Fy!6C
zlK+&nIZX@vDf&4K0hyH^W}XBA5&a)RL96aNXmfz*UO*icMljW`(7vMkv8>8LNM!tW
z?rcluRjaVu1swv}gGl0g+xvfQUDi;}dkdzNN*UI7-S;bke9z52Dn#j=#Y!_H-EmAc
zCY|C#g?v++uUyE05qe>tyvTmGmr=s2MYc!T^<S<RA0Lm;>kNvoP@oeUuH$IyhMo)+
z106HmMvgq@5RaW<_m1tH%AUw4hkY)i_eLn{6&0tyb-4N~)e1j*mPv0tgf!^y=7cbR
z>J%oA`677#nWK}{y2<aUpM!vVF^d^EW{gpF{v^8-D73`LDRKZYS!t~N!q)miVy~#n
z1O-F+hKZc#!ck?Ke$P;e+!iJiaCYQa@$Uk5MlMWU>jgm}<c)w%bmq)04owN8Me06P
zq=*6(cXx}5P6qb9*xFd_p{FTaM#R(*E2{x0$WI5jKXq1Ru5GF)+OQ`<;7*TU#E||g
z1CPe(X9tf&05Fa%Y2;X#K#!y-mSO+dr0Ow!?DJZizD)^JK#}7mn!9ig<_<2;FV((F
z1AsCX?}B=wSHs>mw%u*A%ys1Xf}0{UR$ld>8@L#?fdfm6J;<mRDPdwl1Kd`9s%~Wo
zc?FG?9Uub5Wzb>I{W+s0DT0a3Aenc`Ts}Y3ta%L)yE8V{v8v~pvQ1IJdT_s#tEY&i
zzq)Kpu2UzA+EcLhq+TVC{qBl9ovgZ|wWX$2r0+;@hfkVxyk_Q2)8TS)n?4w!4;s*+
zLkAQkW{dy4zvf|--+xsI<P3~(N7k4LRSK0H#LScV%YHuKqW!%0?+?Reg%j-z6T(hs
zFLjpd-`|*2bBJh@p=Fb!Dc;(Vib{+`_<u0_76(M#Ll#<aWpyX;s6#LfLV*Ag(dZj4
zLlJNW_V2f;Fo+RUSy3xux$a}4fi0PSJ5M}<bm=A9s&5xd$b8?w5MEG-TWV%Tr-{!B
zeEWrF_`@(oQ0vdBWX#N`#?2;MYyO_n8je6;RwC88^H~-|)RjNxE9X{o0U>E2QJ*=p
z<NbU0J}Xr0g?Tc@npaW6@TCD@8l{@^Zp5}&AM0ReXIBp2KMb7ax2co~kr?F81@|b)
zDxiOLFKv;tA;H8J8a!CnvI-s}sJp!UEJotC^A!4>|G{@ZIttl6HP?N$#SZRE(-8V_
zmSRf!4-`Nv(rJLa^Y4n|DrwLvN50J*Dis<4hDqL)mQJUfItzIhb=7vh@`Taa+96~#
zI-ov62m;F)M}z|h$X&Pxt*;m-BbaHbU3%#<+3A|^mIvvJua!cpsvw?)!hzMRCwcJd
zRj`nma1iB|yVVYsmX>DlmM-UZC^fM^#J0)@tSte3KG_=gG8y7*KRKCpFk}{`5W(us
zAA12X3w%x3KQNLT4R%*F(goY{yWmgq3zJu1n$SH;JRX%an*o=xmAJ(*LU<1d%)aE~
z>#W!<T`$DMqyz7Y{#InuqR~?2cxP@gWc5xMQ0}QTWHn+?9mT`v;<nMa4XP!`YC)CF
zo;@8Vlt9W^CZdZlHsHahu`^mc++F&}<+@&WY)TmrbLDcD<+z@zdfwei{+gR+6xK~D
zEueR{s^5UDf`x{rVQ&|oj5V}Yag6YZjZh0UoL+cfXd79>f}sks=e!4N7-mnkJ)@`F
zOJ?x!{rfzpZh5&)dsUIP@6bz08=T+#c=~zARqL8fvE!di_@HfFbC)XpQ`XxX)8@Q&
z-FfNK7-86`-SE1tu!d-9uJI5ZSO!C76N@%>LO?o1NogNyzm->C;&F}*5Mzcuzwhpd
z33@<4&yUWUkw(G7w)G?}1OnZ1()goOr$Ov_(SJiP6$m!>Q>|g__bxx0W{k-%E`CPj
z@1JiBjL13U<l@q{8zMR~j3upq&cD#(#8W-giT2Yx!H-asRDAq+_QHiu4IA6+@b&E?
z-j0sE>(`Idt7x)#vHQap{L5;V3{<c(0Peb!W@ldI)2N($_RbeQJ|xBBCJJZje*-ok
zi%3;XZ6ha*Aa?J*IUOGUd&9=P>Tv9IiC(sd#vLSxhfJP#Ox3nK68AL-Yx&kYTUAO<
zXDoa7i1Qbn?WW-0^Kx^q&#oyoUa;U#2L}hqrPiaw44h7uv<?$gV)T4biauR)cg<B9
zhMdvi>%DhxZNtk8Ugwl&Tmm!cDPfhUs_*sgEE%0VMp3b}OAo+0?>=UtUv=Hw<@>z7
zuT9&{F<9FE|Np<}wZ=&zvkZN_F*VnaLaOr)Q(avfhZO4j^q%45^u=?v?BN}ucjO%J
zqL_jCy<TdM_#&_=<kYDQRJQFE1H3yA@W+zAx5D&$h4G1rEVuH?${z*IK*_O%)p!j^
z&sC9;AK>mJZF*0dY#abqi_?AUdu@96R?-coH0SWN*}F%SMi6HJ!jkL0(kc;iGNFe|
zN#uw{G84ak`C=7sVe8;<tF($ig=J-HkOOv5KhFI(e7Nu0Nc2ARg^x6P+`e7A&$z)v
zjhfMZT#twc6Bc6U(<5@VyPKPn^~DDRxY%L)*BjKYHI|(i`_m}m>6$Q%8vy)v`t6-L
zCHu(52k``kx1nDBl2%ntP*y%9WYAn|lDL)-)7`jH_|@ITWg-Yn`=Q+r?>Nl`t2O{X
zZft6BlyKtfJ2*tg?$AA(Vmo{OJ2xif2xtTZWZ}&%v009~t0ZCti80YyT3P~&zTLKK
zF9m(w_MSbL{ds*KAfHuYAx`$4S(abEeS7}u)l+0JZA4<h@)zmM)t65YAjQ61OE4;W
z{Z6rTR??mkIw9E=XQQI_ai*GBT0SL1qH`bu^}&r-M*8~cn|<h6-}iLQO7;WD%C+sn
zJO*TjBjK*cUhgkqT1p1;8!Z+uW=_(3`fHg+46$frZ$A=zTUTHI#y$sLNpxO8AbGi#
z{lPNlMll4}lRzRNtnhWz`*soM-lk3Coo^`Ic=Zfi9zTchq@*C&3baa&>V5GDVP?}z
z<`%dVyW;XUg5<hS(?jnE`E!>)%EXN=O*A|qUudLO;`cEW+Pi*xy>MuTi@XQ?jH5wi
zeqT#VONjhKbaj$5iMnbMqDek0Rud?;o>5%PVK=9gbvtYyapy>C6^9DA#!hM@7|Oe&
z6N`km3K1R||MGYXAK#8s`PAgTXt?SRdNl1M7<TK^XAT>kgk=DcL`~z2_c`U&mn%Mf
zDkB09C1W7%0~u>f?YjLoI71_2FREqprMrSvQoaaZCCJ2q<Hj9Ey&sE|ntu-X4pT8Z
zue^GL&HDAjXbCG}Rt9DVzRbR=#ff%fzNJf}?yRMIvny87QE-!HN-^M$vd$0s?lPF@
zs~sa%KARmYErfEW+y4D=>yXe37xwoUIbZxO9*--XR1ezO*bLy}B+&T4F5=Lpk#v~4
z`ipWKMV;MJrd#uG)FeMfcq)W~!gi39VkH7Tz4XcjxRErKO-w=l63D4h54=+Us%CqB
z`9ejV9RX>wCD}PS_Ym3X!>V(}BKWAI3(hSwaENEF-FqBf%3!Bq#DZ(K)wmaVv^J|w
zoY+Ufa8+iud_=<$WJmxhV1k3WdC(1EgXD!!RKK6=uHxdQ7&l;L6OeUSxBbZA<@;)c
zdpt~wHD3x;qI2VwcGrG>x^e#8xd9wcFpbn;E1U)&o;GLqh@}~CHi(fF@S3;5dR^I#
zoOXKYb-gFPa*`GKCGzhT4Bp0_sf%IXy<0-qYcXX`8m?HUxB||;WRxaMsh)n_?4ka9
zX*^PaGZv!aIdeLBG`;P=?c2x+ip4W|r|Q6zrXjSMi06U~IlWwc=~$YpIqp8)`o5Hd
zXA7QvKE8uhY;0!FznMpm{-c=fUq=dM*U|SAVCLU-^=N2oto2n6`x(7-8F^Kwz>8bh
zA{y+-JI34VWM4hl7>wfmK9zzxZOn|>;@4_7TU?uZQz!(*5DQen77Ul!?Ck7!=+JEr
zCz#MgV$0T7{Veaw#fbKNGJ4pk{-3X?9ybu`F3?Ygf~AUQ$F>)M98$0%=Ryae{p4=a
ztQw<HjAX)ySMpbwYSiPDoc*<cP}1;tSqn?4t#DEB)UR15{l^~zd4yzPBI%$?o~)Me
zVD8!R2DKWPElflp&Nk}W5X)#OY~G+LpTs**&E7=uE&#^Zg*LXf0|yS=fk5Kw#`ljd
z_%6b&fhf0k|Nggxwgp9^9|k7z+vEIcY}Y@erQgFQe{Fc7sN>ypOt%`bg;5nOWD?49
z`90D<OYY4fI{?;X+-yI`ZQmJVa2a@6EC(pCOB%E(d@=uCkuhE<kWZqR^6fvRczt>|
zc`3G%Id#Lh6Hm`ZM)pzo-BW<_C_a<xgKL&t(ZXDSyrl*6?LJ4IY(w#!wpPv5%xsmD
zoXn|+J907ephS>Tk(~GIGak~5;$Wr9@fQ#Z+ye<}Ep&_bAVE^<Zea@|x%P2P_I_(0
z(XcOp^ckd~5kz#kz4E%&5aK8p0P@h4|GBogx#<3xtRcMG$(!3)!F14m)K+;@nE}*8
zW7fg#UyjxYIDE(0z*9e4dA@C42bHb0w|@(X5SNhPrg|gN;DlPWU)2@Gjn(3-1_cE%
zg!M>zx+P!($0!_8Qo{QWAFi%?=ohnX8ARUg)pMR$Hr={$qmQTir|VvMdGnK7n`_)v
z64D9&8}T#=Tj|}P+9P7i(Qn_J(Vb4#{)CcU_!41^?h-7#MZQOk?1+B1P3H6G8y|kd
zNv(}gVB{z4x4u;TV69qqUfvp~&cE*8ig=5k)yp67?FMqx@@eLA=xp3rXFg}$^Ru&V
zw??cPi?39>*VMBg122Z;<>csMw30PWcWaQ8>CHz38_U)2Z0TfrM-$6?>`mTX{}A#P
zj4-*<9is$9Gt1GV2Cf}Ts{?Pqq~>c0mN$<-=!1{bx5s9Q59XfXDf%5dR-QCIeu_Y*
zK7Cq;p$)Si@0fp2sfais26v22Zcs#)aXWBg17HXf8;!4h3VDV2py&SG+=9nx<sUz`
zJ1{zjYR>=2kyS5hXS8oS+Gp6()jZCo3S*B_9>*n36D=XtCG7Rv^14l8$<^)6ZWrf|
zobvN!2JNe(NTw&lXpFK<%856|X#)){m5*gzF+me&iZ)ZxMnJ@^;gbDQWE-l81bYzP
zZU=_x`eNLS_JhP5hR4SrUBC5hA27$mh}m7YHnDoqFr`6B&7L{)!RnY31`#oqcP^zM
zH#S|g=;hcsZ?%4YygIq@oxkv^Q9Gffo&2C%Xn5_t8=T8G`c-VM9oZ-33u>qkj+iB}
zqsIOKX`8&g=t2cIB+!wps)6$%ol4FuKyl(s>tm&xxqaR^KX4eV!NGds<YK`!@eAeU
z<)>p}4x;QJbBnoFbH4_3<+>jNsB0nfpZ6c59sQ2`B)ab+vny)Hi>9d0`<FHSMe3j>
zYx*v&-1S&#urT1IvFRS&FeQD)TqN5OGitdY6K^yb<MDa!?3(#(90ix6R9WN6p;dv3
z6ODF{C_WG11HC~yO3Kk#Fv9Hm5^L*WfN!K&b*vrBx^_Kfu8$#+DjvaB6s=&`hzMoS
z-8ae_ei<9o`_i_9mR8^seARwp-Aq8MCG3vtzy*TAOfRUu5_V*K*Dd%oHx!n%q&XzM
zS|biA9gYqMTtIKuru7koDLqCypRB(25RMoeJVc7iz$>Q~B^cjZa8CJtl`xDUP5~S+
zI`-S@sw64kDL)4GK<l7?T|p?&)9)_t#HN-d{^F-Th0?0bSM$?rQlIe}e)&Qi1I)nz
zya^jjJfm&J=a{v7gtjP>$?wKlKkG{vS#)L9tQ4)6_x^lQ?Yy0%ZanvsKOQl1<R&f4
z8M>2=dBqS+P!=q8-A#Y>PIR8Yhrr}jiAOiyG+DoXtYF+=0txw789KnRl}6WvkMB^G
zbLW-s4->pF8}u?RCjHz$j)T~7MU;B;hs@%MQHhCx;J302n?n;temoBH{f*a@rqHY`
z-97DkL^O;d*QIA~4tfkUSOlpdEZlmE2jPeu_SRLL&)}>>z%3~8RWKYEqJM;0NT9}R
zggiJDdUX4$OwzFBpzo-F-#Te5!X^Nev}|$7j{tk0h0ES-P#0aYFL8n6SSP;LXX<nn
z$K|?PU+0It-P?JFMQ`Kwr$7tq87^PPV9w7_bPTnAO3rtkmYG?&1ecVQlzT7W%uG*d
z`8rnklZgKsg-m`(ZQV!+cmyW)N7t?sHO@SkEAr85k&og5gCIc|UXeb#yxROH8e))m
z6_!;~(}C+Md^EbQ&&I_aA=Bg((dXN8qiJUr<it?pkhL?_WUe2OGTr0L9j29?0ZZ<}
zd?-3v`szyq_z|i!4gk?Zjw3C3dq1s&SL-6l$64IM9|OO25~|bD(LN+FoT?8WJOKIF
ziASB=mlukqG86?tZtm`u_svlW&@GK=z;XKut^N1!Q1OcriYnyAS~cPrU?M4IuV0n8
z=h#S3g)SY)(Z#ScDv-Y`f?$MYFc+2+2XdpwVQzzZdZvNFU{Eg#1-l(iGBuHR<gPz>
zvw=!9<3ANl0FAqIq+q!GpFK08ez$RO_<g|S$$hEW$f$+cvt2yr2%m%#N@ZP<dj4^H
zC#SgOrneP0>ijClGwYS<8!DoLm637WDP^jwdr$hVr=^}3*hTHeQ!V|!-*o*gvb>_N
QgnwquFw(s=-RiIZ1%zqWh5!Hn

literal 2563
zcmWkw2{=^iAD*%kLo*~Q{h8ddU%NzHGqyA#I~iRXbQuwmoiN$5jTlQzBf}kA*_SbP
zS;G*rWS!y4UYKnE<DBQ5?|IJi`+nzpfA8;ozc1DlX~4@Z!VLm}c#RC<=D_$6cr-cL
zfj5U>jT$hp|9#H@4mv%%^C_e>5QwYJ2(D`pkh7i>X7X@cv@4HC)>9x3n$jD$fD<!M
zaVxwWGmqmv%}h14saZ{+z88Q%h>Dd7QErL@6R)Q(CMuah{2oU<>9`!Ga0vl_BJtE>
z2oAqTHbk_To(m%EPl|n?bAjGmb|i!zA9ZVOYB4mxIr}yHO&k%hZFY^Dd!p`a<q)ra
zT#6dH_w~~Pwx)!><;WQfb_UZ;BR!b+sh@kMJ~T6v*wUiorYy4JdpqH(rVvEZB)rv0
z*u94$ETE*RscB|u87-=Y$!};7QBhGTB9q$}yOV&qcQy>jV{I;q|C5vB1F~T!J-SM4
zE0dFxLAw0VoQR+(>uu!>7+2WwuS>{xcVa(%QnawNY+dzjb5dwhIGY6gNXNwmCN3`S
zq#!)UU|{FwoQf=@e0_aeR(%(ky9_sXcU?U_mTOvptl0xG_O%#7VIh01*Y}v4ahzaq
zchDvU3q(msNU%WuKw1RGflMz*86psDm6es3FJCqS<tRSk9xt(NZ^wcwo0}c?wi_B6
zY%DE1BO}kFU0f7ku&BuETe<$T^JR}354Ps@CcA@U#U}@R)vrufx`H7P2v@igx_@qO
zed^_67jHdQ1WKayJR~8#XBzyjs9~;HN(pAj*tXq0InaxpUs9sYKifLoq>_f+3arS>
z-C|+=O0I*S$Brw!Z|BNLNiA@Q1?%YYt2=jH5E2wLONG_P2sh*9GU9{wH=fwo2>%My
zdj9<R<n*)*R=u;%XPH|`QSm0{c}73JoqBp={|19~rKd|<s|W3Q{5$p(1FT^*`v6tu
zNa@}G9+npx8tOm|)v~v@$I^2B(Z;DAUd_yyr`1`-#l<CMW!K)k7C@xJu4(v0=@V1Z
zB=0YD#$TeWKV3oBR8=WxY2^ZuOeWKZo+I#~xHwiU<fj0DZg+Q=vii3h+OE_pYXGf}
zK=h=$>I<gAKIZ3#pPZcZ4-TUJ{r62V`_uj3zS(7&NCMYR>njdD>;Lp{IGibo_NTOT
z{G&&Y-p$hScwS!K;KTJ=4cEb(c}e#?LQ_*tL_~xGC9s-Q?D0%#YHA9P!*v`S1gNU0
z@S;9|y%##T0HpNg<(C(=>+oyi<*sgSH)VzTDmt0shkNV9ygaclxteM=D=p}glfz9Y
zYXtS=H`!w$IHjaSCd;I-)!g9}eLEM4#4ASky_HQv<t7b}jI?cSJ)<uzU4VY2;v4-p
z3y~z^*wb`WE_i1#>8?3>j1RXFH{JXN2jfcn<AJ`Fl@K>KH#F$n#L!Uq$cTB4c8G{-
zHgAU9qn(xJR~+DlpR2UjiHXsGq$1gPwkBQkZWTg6a5K1x$pvlZ5@d<NXyoxe)vm0n
z!hr$~_dGv@4$D{Am9n+~K$8~ye@wKltgKvOlB=uJhlhuc=CNWpoKb3jU!OA^&WbSe
z9sd!EFSo0X>*?v~l?vi|Fi>SzYC#^GuKvfmd<_+Ltir;=!d&Wse1N6GWQBOxbe-JY
ziyP_UwT)?)<Q!VCVxdKI@3-F&4U*oI$vp1v?vyfln+nj?t5?;Z&4s^Te2}Hg1@trm
zfry^k(?)B*++p?e^D{&utt>LtewqLR|J86!P>>B&7kb>VAVMF^jql5lZ}ePqMu&uG
z8H&rvA-x~#eIO9df%rv4zO4m*dAoUhxGy-)FDl9t-?>mTu#S=hMMX#7lIj)x8K?#N
zIoY+&w1px%Re+C#`ua&vd$BX0IR^r$vbA*N7ZU@60)GE=FJCJyEd_;t{%qK(!pg=*
zpipMcLWAC#CbbMlO?t3zZTf7^4-*YQ;W!+$Q$<-#ZQtH#q`Uv!^33RHL5URjWqq47
zb7xst^M%Y{uI9Kj{l{b_Cnsmi$cKA^<FwJyQ7`Tv&S-$CDmz)g0}AUD=Dt@gl*ynO
zgSmji7ytmnPfN4psYe%qervCE{Q?5<VOb9MHw`ctOdhG}BGA2pv9Yn3so$d_?#d$l
zIT{e4HQE??VP!=XfaO7jgwBIxgm{{d_v{NuB#B$MxIp~E!kiRGaQc5JJSupbtt(L!
z%PXqtw9pY-_ng7!<Kt8DV%-i=eI8_FVPz#QFJJyWTziR5A8}FS^#d<elS)bwKmoh-
z`~CwPROHvOF*APFcenkFZT+;nc9#1zg7$u#WF3aFsd>%bDI}4kiM<>jCp|`2eUDa2
zdB1s9zxVGlzCd)sQ2(87+SL2pAQFkpBj)bOU+l9rHcbJ44xrE2Shi;9IW!030hu2$
z8F?Jsq7x?MvEHE}Q_OUonWLX}RQtRWP&uH-)F<9<&NZe*M@OG2@MtGKm^*|uqp!2S
z-$W!)UsY9=PiVtdpsnzJcG0lg0Z<w~6ak!)n|r_1s&8k<=Rup65^dSXd-Bth3w(U^
z$x2r&X^0~k7HePQL7-Ca-ZM77LmAk$lQ4ra0R-m6;fjiiA<cB{aN6n;z~p6qeuW7Z
z2K7QI_aKmFmp%Y2W2(k;&C9vR^`nhUj;6oduRx00jip(^#%t=&!kxPk<{LMg`vCC=
z(mIQXb#<8p0zo-@U{+`|5ET=$_ElE$#tkzP&Aa=}RhE#<t7FrR{z>sdvf`jJ``R~v
za%NwsG-bUogW~h@zU1d$N9Nv6H+XYp8=oPOBzBix<2hACTg^;Po-=y<)#>wd-%fO;
zcW0cSMX^OfeEcFo539dkD+v@cC+7+~8=Lph!S*yIsMet#HLGV?U0v+}kiv5Qyqkh>
zoja}2%sW5_inw++G$;O&d?8SeM12YB;f7{SZEc@M!qFe2hMzuuOwsh;SX$~yLs|S2
z{A<j*<llb-pSb6jmL>uS7UKk^fn2~=LK_*ev9N%XXvH1L8$7j3f2G14ngRr<$9tm2
z#>S>7R8n#>e?dXP^&2<XD$xUNfVE*M<O$$_0Jyt_wKeh6r;BHvQvJ{Pb)$=^X=y2V
zd@R5++`ynu>Cb&AHT5Ay3y50;+Ucgeya58y8UINCRcrXfgpF{bJ{*B~3}EHt;)>8G
zzVse!VocY!&+FLRi+ZSv{&=BTVkNyXUC&|$C72g0cz#pv&sGBiyO-hp`v^F^gPb)>
zR@tQjB44BrWET||x6U;B4^K?I#)8G#Ftx<eQeGqpvcA5~j}Xw+)iouLaq{uSIF@(Y
zdr{bk!<=b7QGnxS4xnK$m=19`VboXuYyWtH+NqX4ayIjp1aM3e5DZ^P5!%i&WrAPl
WS4A>nB<}&II*<_p2`~H0Df0gu+5yb~

diff --git a/pybamm/expression_tree/binary_operators.py b/pybamm/expression_tree/binary_operators.py
index eec915f832..520ebbdee5 100644
--- a/pybamm/expression_tree/binary_operators.py
+++ b/pybamm/expression_tree/binary_operators.py
@@ -252,6 +252,10 @@ def _binary_simplify(self, left, right):
         if is_scalar_zero(right):
             return pybamm.Scalar(1)
 
+        # zero to the power of anything is zero
+        if is_scalar_zero(left):
+            return pybamm.Scalar(0)
+
         # anything to the power of one is itself
         if is_scalar_one(right):
             return left
diff --git a/pybamm/expression_tree/symbol.py b/pybamm/expression_tree/symbol.py
index 5c9bdb2ac2..2a069ba786 100644
--- a/pybamm/expression_tree/symbol.py
+++ b/pybamm/expression_tree/symbol.py
@@ -635,6 +635,8 @@ def evaluate_ignoring_errors(self, t=0):
             # (there is a e.g. StateVector in the tree)
             if error.args[0] == "StateVector cannot evaluate input 'y=None'":
                 return None
+            elif error.args[0] == "StateVectorDot cannot evaluate input 'y_dot=None'":
+                return None
             else:
                 raise error
         except ValueError as e:
diff --git a/tests/unit/test_expression_tree/test_operations/test_simplify.py b/tests/unit/test_expression_tree/test_operations/test_simplify.py
index 129aef5e98..2cd37f7537 100644
--- a/tests/unit/test_expression_tree/test_operations/test_simplify.py
+++ b/tests/unit/test_expression_tree/test_operations/test_simplify.py
@@ -16,6 +16,7 @@ def test_symbol_simplify(self):
         d = pybamm.Scalar(-1)
         e = pybamm.Scalar(2)
         g = pybamm.Variable("g")
+        gdot = pybamm.VariableDot("g'")
 
         # negate
         self.assertIsInstance((-a).simplify(), pybamm.Scalar)
@@ -175,6 +176,18 @@ def myfunction(x, y):
         self.assertIsInstance(expr.children[1], pybamm.Negate)
         self.assertIsInstance(expr.children[1].children[0], pybamm.Parameter)
 
+        expr = (e * g * b).simplify()
+        self.assertIsInstance(expr, pybamm.Multiplication)
+        self.assertIsInstance(expr.children[0], pybamm.Scalar)
+        self.assertEqual(expr.children[0].evaluate(), 2.0)
+        self.assertIsInstance(expr.children[1], pybamm.Variable)
+
+        expr = (e * gdot * b).simplify()
+        self.assertIsInstance(expr, pybamm.Multiplication)
+        self.assertIsInstance(expr.children[0], pybamm.Scalar)
+        self.assertEqual(expr.children[0].evaluate(), 2.0)
+        self.assertIsInstance(expr.children[1], pybamm.VariableDot)
+
         expr = (e + (g - c)).simplify()
         self.assertIsInstance(expr, pybamm.Addition)
         self.assertIsInstance(expr.children[0], pybamm.Scalar)

From 00990f9c7133120a2065609e39f6dc2993be6594 Mon Sep 17 00:00:00 2001
From: Martin Robinson <martinjrobins@gmail.com>
Date: Sun, 15 Mar 2020 07:23:36 +0000
Subject: [PATCH 21/30] #858 improve coverage

---
 .../expression_tree/independent_variable.py   |  5 +--
 pybamm/expression_tree/state_vector.py        | 16 +-------
 .../test_operations/test_convert_to_casadi.py |  5 +++
 .../test_operations/test_jac.py               | 41 +++++++++++++++++++
 .../test_expression_tree/test_variable.py     | 25 +++++++++++
 5 files changed, 73 insertions(+), 19 deletions(-)

diff --git a/pybamm/expression_tree/independent_variable.py b/pybamm/expression_tree/independent_variable.py
index 008ddf6113..76789dfc08 100644
--- a/pybamm/expression_tree/independent_variable.py
+++ b/pybamm/expression_tree/independent_variable.py
@@ -35,10 +35,7 @@ def _evaluate_for_shape(self):
 
     def _jac(self, variable):
         """ See :meth:`pybamm.Symbol._jac()`. """
-        if variable.id == self.id:
-            return pybamm.Scalar(1)
-        else:
-            return pybamm.Scalar(0)
+        return pybamm.Scalar(0)
 
 
 class Time(IndependentVariable):
diff --git a/pybamm/expression_tree/state_vector.py b/pybamm/expression_tree/state_vector.py
index 8e40b23400..8abae5414f 100644
--- a/pybamm/expression_tree/state_vector.py
+++ b/pybamm/expression_tree/state_vector.py
@@ -105,20 +105,6 @@ def set_id(self):
             + tuple(self.domain)
         )
 
-    def _base_evaluate(self, t=None, y=None, y_dot=None, u=None):
-        """ See :meth:`pybamm.Symbol._base_evaluate()`. """
-        if y is None:
-            raise TypeError("StateVector cannot evaluate input 'y=None'")
-        if y.shape[0] < len(self.evaluation_array):
-            raise ValueError(
-                "y is too short, so value with slice is smaller than expected"
-            )
-        else:
-            out = (y[: len(self._evaluation_array)])[self._evaluation_array]
-            if isinstance(out, np.ndarray) and out.ndim == 1:
-                out = out[:, np.newaxis]
-            return out
-
     def _jac_diff_vector(self, variable):
         """
         Differentiate a slice of a StateVector of size m with respect to another slice
@@ -139,7 +125,7 @@ def _jac_diff_vector(self, variable):
         variable_size = variable.last_point - variable.first_point
 
         # Return zeros of correct size since no entries match
-        return csr_matrix(slices_size, variable_size)
+        return pybamm.Matrix(csr_matrix((slices_size, variable_size)))
 
     def _jac_same_vector(self, variable):
         """
diff --git a/tests/unit/test_expression_tree/test_operations/test_convert_to_casadi.py b/tests/unit/test_expression_tree/test_operations/test_convert_to_casadi.py
index 9dbd48492e..5746aaa9ec 100644
--- a/tests/unit/test_expression_tree/test_operations/test_convert_to_casadi.py
+++ b/tests/unit/test_expression_tree/test_operations/test_convert_to_casadi.py
@@ -212,6 +212,11 @@ def test_errors(self):
             ValueError, "Must provide a 'y' for converting state vectors"
         ):
             y.to_casadi()
+        y_dot = pybamm.StateVectorDot(slice(0, 10))
+        with self.assertRaisesRegex(
+            ValueError, "Must provide a 'y_dot' for converting state vectors"
+        ):
+            y_dot.to_casadi()
         var = pybamm.Variable("var")
         with self.assertRaisesRegex(TypeError, "Cannot convert symbol of type"):
             var.to_casadi()
diff --git a/tests/unit/test_expression_tree/test_operations/test_jac.py b/tests/unit/test_expression_tree/test_operations/test_jac.py
index 270b030c44..a45a32ca5d 100644
--- a/tests/unit/test_expression_tree/test_operations/test_jac.py
+++ b/tests/unit/test_expression_tree/test_operations/test_jac.py
@@ -109,6 +109,47 @@ def test_nonlinear(self):
         with self.assertRaises(pybamm.UndefinedOperationError):
             func.jac(y)
 
+    def test_linear_ydot(self):
+        y = pybamm.StateVector(slice(0, 4))
+        y_dot = pybamm.StateVectorDot(slice(0, 4))
+        u = pybamm.StateVector(slice(0, 2))
+        v = pybamm.StateVector(slice(2, 4))
+        u_dot = pybamm.StateVectorDot(slice(0, 2))
+        v_dot = pybamm.StateVectorDot(slice(2, 4))
+
+        y0 = np.ones(4)
+        y_dot0 = np.ones(4)
+
+        func = u_dot
+        jacobian = np.array([[1, 0, 0, 0], [0, 1, 0, 0]])
+        dfunc_dy = func.jac(y_dot).evaluate(y=y0, y_dot=y_dot0)
+        np.testing.assert_array_equal(jacobian, dfunc_dy.toarray())
+
+        func = -v_dot
+        jacobian = np.array([[0, 0, -1, 0], [0, 0, 0, -1]])
+        dfunc_dy = func.jac(y_dot).evaluate(y=y0, y_dot=y_dot0)
+        np.testing.assert_array_equal(jacobian, dfunc_dy.toarray())
+
+        func = u_dot
+        jacobian = np.array([[0, 0, 0, 0], [0, 0, 0, 0]])
+        dfunc_dy = func.jac(y).evaluate(y=y0, y_dot=y_dot0)
+        np.testing.assert_array_equal(jacobian, dfunc_dy.toarray())
+
+        func = -v_dot
+        jacobian = np.array([[0, 0, 0, 0], [0, 0, 0, 0]])
+        dfunc_dy = func.jac(y).evaluate(y=y0, y_dot=y_dot0)
+        np.testing.assert_array_equal(jacobian, dfunc_dy.toarray())
+
+        func = u
+        jacobian = np.array([[0, 0, 0, 0], [0, 0, 0, 0]])
+        dfunc_dy = func.jac(y_dot).evaluate(y=y0, y_dot=y_dot0)
+        np.testing.assert_array_equal(jacobian, dfunc_dy.toarray())
+
+        func = -v
+        jacobian = np.array([[0, 0, 0, 0], [0, 0, 0, 0]])
+        dfunc_dy = func.jac(y_dot).evaluate(y=y0, y_dot=y_dot0)
+        np.testing.assert_array_equal(jacobian, dfunc_dy.toarray())
+
     def test_functions(self):
         y = pybamm.StateVector(slice(0, 4))
         u = pybamm.StateVector(slice(0, 2))
diff --git a/tests/unit/test_expression_tree/test_variable.py b/tests/unit/test_expression_tree/test_variable.py
index 47ae081ffb..33f83f02b2 100644
--- a/tests/unit/test_expression_tree/test_variable.py
+++ b/tests/unit/test_expression_tree/test_variable.py
@@ -16,6 +16,14 @@ def test_variable_init(self):
         self.assertEqual(a.domain[0], "test")
         self.assertRaises(TypeError, pybamm.Variable("a", domain="test"))
 
+    def test_variable_diff(self):
+        a = pybamm.Variable("a")
+        b = pybamm.Variable("b")
+        self.assertIsInstance(a.diff(a), pybamm.Scalar)
+        self.assertEqual(a.diff(a).evaluate(), 1)
+        self.assertIsInstance(a.diff(b), pybamm.Scalar)
+        self.assertEqual(a.diff(b).evaluate(), 0)
+
     def test_variable_id(self):
         a1 = pybamm.Variable("a", domain=["negative electrode"])
         a2 = pybamm.Variable("a", domain=["negative electrode"])
@@ -25,6 +33,7 @@ def test_variable_id(self):
         self.assertNotEqual(a1.id, a3.id)
         self.assertNotEqual(a1.id, a4.id)
 
+
 class TestVariableDot(unittest.TestCase):
     def test_variable_init(self):
         a = pybamm.VariableDot("a'")
@@ -43,6 +52,14 @@ def test_variable_id(self):
         self.assertNotEqual(a1.id, a3.id)
         self.assertNotEqual(a1.id, a4.id)
 
+    def test_variable_diff(self):
+        a = pybamm.VariableDot("a")
+        b = pybamm.Variable("b")
+        self.assertIsInstance(a.diff(a), pybamm.Scalar)
+        self.assertEqual(a.diff(a).evaluate(), 1)
+        self.assertIsInstance(a.diff(b), pybamm.Scalar)
+        self.assertEqual(a.diff(b).evaluate(), 0)
+
 
 class TestExternalVariable(unittest.TestCase):
     def test_external_variable_scalar(self):
@@ -68,6 +85,14 @@ def test_external_variable_vector(self):
         with self.assertRaisesRegex(ValueError, "External variable"):
             a.evaluate(u={"a": np.ones((5, 1))})
 
+    def test_external_variable_diff(self):
+        a = pybamm.ExternalVariable("a", 10)
+        b = pybamm.Variable("b")
+        self.assertIsInstance(a.diff(a), pybamm.Scalar)
+        self.assertEqual(a.diff(a).evaluate(), 1)
+        self.assertIsInstance(a.diff(b), pybamm.Scalar)
+        self.assertEqual(a.diff(b).evaluate(), 0)
+
 
 if __name__ == "__main__":
     print("Add -v for more debug output")

From 445ae5445e445397415acb006fce794c50f1d596 Mon Sep 17 00:00:00 2001
From: Martin Robinson <martinjrobins@gmail.com>
Date: Sun, 15 Mar 2020 09:12:58 +0000
Subject: [PATCH 22/30] #858 fix external variables test

---
 pybamm/solvers/base_solver.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/pybamm/solvers/base_solver.py b/pybamm/solvers/base_solver.py
index b6a6a4ff99..f797909580 100644
--- a/pybamm/solvers/base_solver.py
+++ b/pybamm/solvers/base_solver.py
@@ -130,7 +130,7 @@ def set_up(self, model, inputs=None):
             )
 
         inputs = inputs or {}
-        y0 = model.concatenated_initial_conditions.evaluate(0, None, inputs)
+        y0 = model.concatenated_initial_conditions.evaluate(0, None, u=inputs)
 
         # Set model timescale
         model.timescale_eval = model.timescale.evaluate(u=inputs)

From 7ac78480e04be8f310c11b46e217fd1e33955af7 Mon Sep 17 00:00:00 2001
From: Martin Robinson <martinjrobins@gmail.com>
Date: Sun, 15 Mar 2020 10:36:26 +0000
Subject: [PATCH 23/30] #858 fix coverage and flake8

---
 pybamm/expression_tree/state_vector.py              | 12 ++++--------
 .../test_discretisations/test_discretisation.py     |  6 +++---
 tests/unit/test_expression_tree/test_d_dt.py        |  7 ++++---
 .../test_operations/test_convert_to_casadi.py       |  3 ++-
 .../unit/test_expression_tree/test_state_vector.py  | 13 ++++++++++++-
 tests/unit/test_parameters/test_parameters_cli.py   |  1 +
 6 files changed, 26 insertions(+), 16 deletions(-)

diff --git a/pybamm/expression_tree/state_vector.py b/pybamm/expression_tree/state_vector.py
index 8abae5414f..01b2bb6eca 100644
--- a/pybamm/expression_tree/state_vector.py
+++ b/pybamm/expression_tree/state_vector.py
@@ -242,15 +242,13 @@ def _base_evaluate(self, t=None, y=None, y_dot=None, u=None):
         return out
 
     def diff(self, variable):
-        if variable.id == self.id:
-            return pybamm.Scalar(1)
-        elif variable.id == pybamm.t.id:
+        if variable.id == pybamm.t.id:
             return StateVectorDot(*self._y_slices, name=self.name + "'",
                                   domain=self.domain,
                                   auxiliary_domains=self.auxiliary_domains,
                                   evaluation_array=self.evaluation_array)
         else:
-            return pybamm.Scalar(0)
+            raise NotImplementedError
 
     def _jac(self, variable):
         if isinstance(variable, pybamm.StateVector):
@@ -309,14 +307,12 @@ def _base_evaluate(self, t=None, y=None, y_dot=None, u=None):
         return out
 
     def diff(self, variable):
-        if variable.id == self.id:
-            return pybamm.Scalar(1)
-        elif variable.id == pybamm.t.id:
+        if variable.id == pybamm.t.id:
             raise pybamm.ModelError(
                 "cannot take second time derivative of a state vector"
             )
         else:
-            return pybamm.Scalar(0)
+            raise NotImplementedError
 
     def _jac(self, variable):
         if isinstance(variable, pybamm.StateVectorDot):
diff --git a/tests/unit/test_discretisations/test_discretisation.py b/tests/unit/test_discretisations/test_discretisation.py
index 787aad46cd..cb0f015a99 100644
--- a/tests/unit/test_discretisations/test_discretisation.py
+++ b/tests/unit/test_discretisations/test_discretisation.py
@@ -703,7 +703,8 @@ def test_process_model_ode(self):
         # test that any time derivatives of variables in rhs raises an
         # error
         model = pybamm.BaseModel()
-        model.rhs = {c: pybamm.div(N) + c.diff(pybamm.t), T: pybamm.div(q), S: pybamm.div(p)}
+        model.rhs = {c: pybamm.div(N) + c.diff(pybamm.t),
+                     T: pybamm.div(q), S: pybamm.div(p)}
         model.initial_conditions = {
             c: pybamm.Scalar(2),
             T: pybamm.Scalar(5),
@@ -846,8 +847,6 @@ def test_process_model_dae(self):
         with self.assertRaises(pybamm.ModelError):
             disc.process_model(model)
 
-
-
     def test_process_model_concatenation(self):
         # concatenation of variables as the key
         cn = pybamm.Variable("c", domain=["negative electrode"])
@@ -1144,6 +1143,7 @@ def test_mass_matirx_inverse(self):
             model.mass_matrix_inv.entries.toarray(), mass_inv.toarray()
         )
 
+
 if __name__ == "__main__":
     print("Add -v for more debug output")
     import sys
diff --git a/tests/unit/test_expression_tree/test_d_dt.py b/tests/unit/test_expression_tree/test_d_dt.py
index 475cb733d8..f9eb5a2cd9 100644
--- a/tests/unit/test_expression_tree/test_d_dt.py
+++ b/tests/unit/test_expression_tree/test_d_dt.py
@@ -22,8 +22,8 @@ def test_time_derivative(self):
         self.assertEqual(a.simplify().id, (2 * pybamm.t).id)
         self.assertEqual(a.evaluate(t=1), 2)
 
-        a =(2 + pybamm.t**2).diff(pybamm.t)
-        self.assertEqual(a.simplify().id, (2*pybamm.t).id)
+        a = (2 + pybamm.t**2).diff(pybamm.t)
+        self.assertEqual(a.simplify().id, (2 * pybamm.t).id)
         self.assertEqual(a.evaluate(t=1), 2)
 
     def test_time_derivative_of_variable(self):
@@ -33,7 +33,7 @@ def test_time_derivative_of_variable(self):
         self.assertEqual(a.name, "a'")
 
         p = pybamm.Parameter('p')
-        a = (1 + p*pybamm.Variable('a')).diff(pybamm.t).simplify()
+        a = (1 + p * pybamm.Variable('a')).diff(pybamm.t).simplify()
         self.assertIsInstance(a, pybamm.Multiplication)
         self.assertEqual(a.children[0].name, 'p')
         self.assertEqual(a.children[1].name, "a'")
@@ -59,6 +59,7 @@ def test_time_derivative_of_state_vector(self):
         with self.assertRaises(pybamm.ModelError):
             a = (sv).diff(pybamm.t).diff(pybamm.t)
 
+
 if __name__ == "__main__":
     print("Add -v for more debug output")
     import sys
diff --git a/tests/unit/test_expression_tree/test_operations/test_convert_to_casadi.py b/tests/unit/test_expression_tree/test_operations/test_convert_to_casadi.py
index ef8c923bdb..0ebdb31ebf 100644
--- a/tests/unit/test_expression_tree/test_operations/test_convert_to_casadi.py
+++ b/tests/unit/test_expression_tree/test_operations/test_convert_to_casadi.py
@@ -201,7 +201,8 @@ def test_convert_external_variable(self):
 
         # External only
         self.assert_casadi_equal(
-            pybamm_u1.to_casadi(casadi_t, casadi_y, u=casadi_us), casadi_us["External 1"]
+            pybamm_u1.to_casadi(casadi_t, casadi_y, u=casadi_us),
+            casadi_us["External 1"]
         )
 
         # More complex
diff --git a/tests/unit/test_expression_tree/test_state_vector.py b/tests/unit/test_expression_tree/test_state_vector.py
index 3038c4c4fa..c16a038ddb 100644
--- a/tests/unit/test_expression_tree/test_state_vector.py
+++ b/tests/unit/test_expression_tree/test_state_vector.py
@@ -35,6 +35,14 @@ def test_evaluate_list(self):
         y = np.linspace(0, 3, 31)
         np.testing.assert_array_almost_equal(sv.evaluate(y=y), y[:, np.newaxis])
 
+    def test_diff(self):
+        a = pybamm.StateVector(slice(0, 10))
+        with self.assertRaises(NotImplementedError):
+            a.diff(a)
+        b = pybamm.StateVectorDot(slice(0, 10))
+        with self.assertRaises(NotImplementedError):
+            a.diff(b)
+
     def test_name(self):
         sv = pybamm.StateVector(slice(0, 10))
         self.assertEqual(sv.name, "y[0:10]")
@@ -61,6 +69,7 @@ def test_failure(self):
         with self.assertRaisesRegex(TypeError, "all y_slices must be slice objects"):
             pybamm.StateVector(slice(0, 10), 1)
 
+
 class TestStateVectorDot(unittest.TestCase):
     def test_evaluate(self):
         sv = pybamm.StateVectorDot(slice(0, 10))
@@ -72,7 +81,8 @@ def test_evaluate(self):
         # Try evaluating with a y that is too short
         y_dot2 = np.ones(5)
         with self.assertRaisesRegex(
-            ValueError, "y_dot is too short, so value with slice is smaller than expected"
+            ValueError,
+            "y_dot is too short, so value with slice is smaller than expected"
         ):
             sv.evaluate(y_dot=y_dot2)
 
@@ -80,6 +90,7 @@ def test_name(self):
         sv = pybamm.StateVectorDot(slice(0, 10))
         self.assertEqual(sv.name, "y_dot[0:10]")
 
+
 if __name__ == "__main__":
     print("Add -v for more debug output")
     import sys
diff --git a/tests/unit/test_parameters/test_parameters_cli.py b/tests/unit/test_parameters/test_parameters_cli.py
index 1ee50bb4b3..05bb15cf75 100644
--- a/tests/unit/test_parameters/test_parameters_cli.py
+++ b/tests/unit/test_parameters/test_parameters_cli.py
@@ -120,6 +120,7 @@ def test_list_params(self):
         # but must not intefere with existing input dir if it exists
         # in the current dir...
 
+
 if __name__ == "__main__":
     print("Add -v for more debug output")
     import sys

From 345be6a6a8fba1612fae454f47b06cf337a09211 Mon Sep 17 00:00:00 2001
From: Martin Robinson <martinjrobins@gmail.com>
Date: Sun, 15 Mar 2020 12:17:31 +0000
Subject: [PATCH 24/30] #858 fix test errors

---
 pybamm/expression_tree/state_vector.py               | 10 +++++++---
 tests/unit/test_expression_tree/test_state_vector.py |  8 --------
 2 files changed, 7 insertions(+), 11 deletions(-)

diff --git a/pybamm/expression_tree/state_vector.py b/pybamm/expression_tree/state_vector.py
index 01b2bb6eca..b19ee9108d 100644
--- a/pybamm/expression_tree/state_vector.py
+++ b/pybamm/expression_tree/state_vector.py
@@ -242,13 +242,15 @@ def _base_evaluate(self, t=None, y=None, y_dot=None, u=None):
         return out
 
     def diff(self, variable):
+        if variable.id == self.id:
+            return pybamm.Scalar(1)
         if variable.id == pybamm.t.id:
             return StateVectorDot(*self._y_slices, name=self.name + "'",
                                   domain=self.domain,
                                   auxiliary_domains=self.auxiliary_domains,
                                   evaluation_array=self.evaluation_array)
         else:
-            raise NotImplementedError
+            return pybamm.Scalar(0)
 
     def _jac(self, variable):
         if isinstance(variable, pybamm.StateVector):
@@ -307,12 +309,14 @@ def _base_evaluate(self, t=None, y=None, y_dot=None, u=None):
         return out
 
     def diff(self, variable):
-        if variable.id == pybamm.t.id:
+        if variable.id == self.id:
+            return pybamm.Scalar(1)
+        elif variable.id == pybamm.t.id:
             raise pybamm.ModelError(
                 "cannot take second time derivative of a state vector"
             )
         else:
-            raise NotImplementedError
+            return pybamm.Scalar(0)
 
     def _jac(self, variable):
         if isinstance(variable, pybamm.StateVectorDot):
diff --git a/tests/unit/test_expression_tree/test_state_vector.py b/tests/unit/test_expression_tree/test_state_vector.py
index c16a038ddb..38f3cd2c39 100644
--- a/tests/unit/test_expression_tree/test_state_vector.py
+++ b/tests/unit/test_expression_tree/test_state_vector.py
@@ -35,14 +35,6 @@ def test_evaluate_list(self):
         y = np.linspace(0, 3, 31)
         np.testing.assert_array_almost_equal(sv.evaluate(y=y), y[:, np.newaxis])
 
-    def test_diff(self):
-        a = pybamm.StateVector(slice(0, 10))
-        with self.assertRaises(NotImplementedError):
-            a.diff(a)
-        b = pybamm.StateVectorDot(slice(0, 10))
-        with self.assertRaises(NotImplementedError):
-            a.diff(b)
-
     def test_name(self):
         sv = pybamm.StateVector(slice(0, 10))
         self.assertEqual(sv.name, "y[0:10]")

From 0c5270efb5b41c1a64353b51dbb6ceefa1bdde61 Mon Sep 17 00:00:00 2001
From: Martin Robinson <martinjrobins@gmail.com>
Date: Mon, 16 Mar 2020 10:16:13 +0000
Subject: [PATCH 25/30] #858 improve coverage

---
 .../unit/test_expression_tree/test_operations/test_jac.py | 8 ++++++++
 tests/unit/test_expression_tree/test_symbolic_diff.py     | 6 ++++++
 2 files changed, 14 insertions(+)

diff --git a/tests/unit/test_expression_tree/test_operations/test_jac.py b/tests/unit/test_expression_tree/test_operations/test_jac.py
index 7c34949629..7e24c4b66a 100644
--- a/tests/unit/test_expression_tree/test_operations/test_jac.py
+++ b/tests/unit/test_expression_tree/test_operations/test_jac.py
@@ -105,6 +105,14 @@ def test_nonlinear(self):
         dfunc_dy = func.jac(y).evaluate(y=y0)
         np.testing.assert_array_equal(jacobian, dfunc_dy.toarray())
 
+    def test_multislice_raises(self):
+        y1 = pybamm.StateVector(slice(0, 4), slice(7, 8))
+        y2 = pybamm.StateVector(slice(4, 7))
+        with self.assertRaises(NotImplementedError):
+            y1.jac(y1)
+        with self.assertRaises(NotImplementedError):
+            y2.jac(y1)
+
     def test_linear_ydot(self):
         y = pybamm.StateVector(slice(0, 4))
         y_dot = pybamm.StateVectorDot(slice(0, 4))
diff --git a/tests/unit/test_expression_tree/test_symbolic_diff.py b/tests/unit/test_expression_tree/test_symbolic_diff.py
index 68370ddfd3..25c6b76749 100644
--- a/tests/unit/test_expression_tree/test_symbolic_diff.py
+++ b/tests/unit/test_expression_tree/test_symbolic_diff.py
@@ -60,6 +60,12 @@ def test_diff_zero(self):
         self.assertEqual(func.diff(b).id, pybamm.Scalar(0).id)
         self.assertNotEqual(func.diff(a).id, pybamm.Scalar(0).id)
 
+    def test_diff_state_vector_dot(self):
+        a = pybamm.StateVectorDot(slice(0, 1))
+        b = pybamm.StateVector(slice(1, 2))
+        self.assertEqual(a.diff(a).id, pybamm.Scalar(1).id)
+        self.assertEqual(a.diff(b).id, pybamm.Scalar(0).id)
+
     def test_diff_heaviside(self):
         a = pybamm.Scalar(1)
         b = pybamm.StateVector(slice(0, 1))

From fa5373ef505a59f2a24cd56addb2644985a1df86 Mon Sep 17 00:00:00 2001
From: Martin Robinson <martinjrobins@gmail.com>
Date: Wed, 18 Mar 2020 06:43:44 +0000
Subject: [PATCH 26/30] #858 add coverage

---
 tests/unit/test_expression_tree/test_operations/test_jac.py | 3 +++
 1 file changed, 3 insertions(+)

diff --git a/tests/unit/test_expression_tree/test_operations/test_jac.py b/tests/unit/test_expression_tree/test_operations/test_jac.py
index 7e24c4b66a..56895d26b4 100644
--- a/tests/unit/test_expression_tree/test_operations/test_jac.py
+++ b/tests/unit/test_expression_tree/test_operations/test_jac.py
@@ -107,11 +107,14 @@ def test_nonlinear(self):
 
     def test_multislice_raises(self):
         y1 = pybamm.StateVector(slice(0, 4), slice(7, 8))
+        y_dot1 = pybamm.StateVectorDot(slice(0, 4), slice(7, 8))
         y2 = pybamm.StateVector(slice(4, 7))
         with self.assertRaises(NotImplementedError):
             y1.jac(y1)
         with self.assertRaises(NotImplementedError):
             y2.jac(y1)
+        with self.assertRaises(NotImplementedError):
+            y_dot1.jac(y1)
 
     def test_linear_ydot(self):
         y = pybamm.StateVector(slice(0, 4))

From d21f0b604e357d0efd8b8104d19b80ae37fc990f Mon Sep 17 00:00:00 2001
From: Martin Robinson <martinjrobins@gmail.com>
Date: Wed, 18 Mar 2020 06:54:26 +0000
Subject: [PATCH 27/30] #858 improve coverage, document exception for
 InputParameter in evaluate_ignoring_errors

---
 pybamm/expression_tree/symbol.py               | 7 ++++---
 tests/unit/test_expression_tree/test_symbol.py | 7 +++++++
 2 files changed, 11 insertions(+), 3 deletions(-)

diff --git a/pybamm/expression_tree/symbol.py b/pybamm/expression_tree/symbol.py
index d3f2dfba2f..c59c69060f 100644
--- a/pybamm/expression_tree/symbol.py
+++ b/pybamm/expression_tree/symbol.py
@@ -615,9 +615,10 @@ def is_constant(self):
     def evaluate_ignoring_errors(self, t=0):
         """
         Evaluates the expression. If a node exists in the tree that cannot be evaluated
-        as a scalar or vector (e.g. Time, Parameter, Variable, StateVector,
-        InputParameter), then None is returned. Otherwise the result of the evaluation
-        is given
+        as a scalar or vector (e.g. Time, Parameter, Variable, StateVector), then None
+        is returned. If there is an InputParameter in the tree then a 1 is returned.
+        Otherwise the result of the evaluation is given.
+
 
         See Also
         --------
diff --git a/tests/unit/test_expression_tree/test_symbol.py b/tests/unit/test_expression_tree/test_symbol.py
index dbff68dd3a..a8462f2559 100644
--- a/tests/unit/test_expression_tree/test_symbol.py
+++ b/tests/unit/test_expression_tree/test_symbol.py
@@ -168,6 +168,13 @@ def test_symbol_evaluation(self):
         with self.assertRaises(NotImplementedError):
             a.evaluate()
 
+    def test_evaluate_ignoring_errors(self):
+        self.assertIsNone(pybamm.t.evaluate_ignoring_errors(t=None))
+        self.assertEqual(pybamm.t.evaluate_ignoring_errors(t=0), 0)
+        self.assertIsNone(pybamm.Parameter("a").evaluate_ignoring_errors())
+        self.assertIsNone(pybamm.StateVector(slice(0,1)).evaluate_ignoring_errors())
+        self.assertEqual(pybamm.InputParameter("a").evaluate_ignoring_errors(), 1)
+
     def test_symbol_is_constant(self):
         a = pybamm.Variable("a")
         self.assertFalse(a.is_constant())

From d8808a6f9722a605a6ae19b072c8ad1060e08c20 Mon Sep 17 00:00:00 2001
From: Martin Robinson <martinjrobins@gmail.com>
Date: Wed, 18 Mar 2020 06:57:59 +0000
Subject: [PATCH 28/30] #858 improve coverage

---
 .../test_operations/test_convert_to_casadi.py             | 8 ++++++++
 1 file changed, 8 insertions(+)

diff --git a/tests/unit/test_expression_tree/test_operations/test_convert_to_casadi.py b/tests/unit/test_expression_tree/test_operations/test_convert_to_casadi.py
index 0ebdb31ebf..fa9992ad24 100644
--- a/tests/unit/test_expression_tree/test_operations/test_convert_to_casadi.py
+++ b/tests/unit/test_expression_tree/test_operations/test_convert_to_casadi.py
@@ -66,9 +66,11 @@ def test_convert_array_symbols(self):
 
         casadi_t = casadi.MX.sym("t")
         casadi_y = casadi.MX.sym("y", 10)
+        casadi_y_dot = casadi.MX.sym("y_dot", 10)
 
         pybamm_t = pybamm.Time()
         pybamm_y = pybamm.StateVector(slice(0, 10))
+        pybamm_y_dot = pybamm.StateVectorDot(slice(0, 10))
 
         # Time
         self.assertEqual(pybamm_t.to_casadi(casadi_t, casadi_y), casadi_t)
@@ -76,6 +78,12 @@ def test_convert_array_symbols(self):
         # State Vector
         self.assert_casadi_equal(pybamm_y.to_casadi(casadi_t, casadi_y), casadi_y)
 
+        # State Vector Dot
+        self.assert_casadi_equal(
+            pybamm_y_dot.to_casadi(casadi_t, casadi_y, casadi_y_dot),
+            casadi_y_dot
+        )
+
     def test_special_functions(self):
         a = pybamm.Array(np.array([1, 2, 3, 4, 5]))
         self.assert_casadi_equal(pybamm.max(a).to_casadi(), casadi.MX(5), evalf=True)

From bc328835491234a523a9af1f12cfa4a74f3b39d7 Mon Sep 17 00:00:00 2001
From: Martin Robinson <martinjrobins@gmail.com>
Date: Wed, 18 Mar 2020 07:10:17 +0000
Subject: [PATCH 29/30] #858 improve coverage

---
 pybamm/models/base_model.py               |  2 +-
 tests/unit/test_models/test_base_model.py | 55 +++++++++++++++++++++++
 2 files changed, 56 insertions(+), 1 deletion(-)

diff --git a/pybamm/models/base_model.py b/pybamm/models/base_model.py
index a7285fada7..641ef347fa 100644
--- a/pybamm/models/base_model.py
+++ b/pybamm/models/base_model.py
@@ -381,11 +381,11 @@ def check_well_posedness(self, post_discretisation=False):
         post_discretisation : boolean
             A flag indicating tests to be skipped after discretisation
         """
+        self.check_for_time_derivatives()
         self.check_well_determined(post_discretisation)
         self.check_algebraic_equations(post_discretisation)
         self.check_ics_bcs()
         self.check_default_variables_dictionaries()
-        self.check_for_time_derivatives()
         # Can't check variables after discretising, since Variable objects get replaced
         # by StateVector objects
         # Checking variables is slow, so only do it in debug mode
diff --git a/tests/unit/test_models/test_base_model.py b/tests/unit/test_models/test_base_model.py
index b9c5e6c703..2f2b37d05a 100644
--- a/tests/unit/test_models/test_base_model.py
+++ b/tests/unit/test_models/test_base_model.py
@@ -248,6 +248,61 @@ def test_check_well_posedness_variables(self):
         ):
             model.check_well_posedness(post_discretisation=True)
 
+        # model must be in semi-explicit form
+        model = pybamm.BaseModel()
+        model.rhs = {c: d.diff(pybamm.t), d: -1}
+        model.initial_conditions = {c: 1, d: 1}
+        with self.assertRaisesRegex(
+            pybamm.ModelError,
+            "time derivative of variable found",
+        ):
+            model.check_well_posedness()
+
+        # model must be in semi-explicit form
+        model = pybamm.BaseModel()
+        model.algebraic = {
+            c: 2 * d - c,
+            d: c * d.diff(pybamm.t) - d,
+        }
+        model.initial_conditions = {c: 1, d: 1}
+        with self.assertRaisesRegex(
+            pybamm.ModelError,
+            "time derivative of variable found",
+        ):
+            model.check_well_posedness()
+
+        # model must be in semi-explicit form
+        model = pybamm.BaseModel()
+        model.rhs = {c: d.diff(pybamm.t), d: -1}
+        model.initial_conditions = {c: 1, d: 1}
+        with self.assertRaisesRegex(
+            pybamm.ModelError,
+            "time derivative of variable found",
+        ):
+            model.check_well_posedness()
+
+        # model must be in semi-explicit form
+        model = pybamm.BaseModel()
+        model.algebraic = {
+            d: 5 * pybamm.StateVector(slice(0, 15)) - 1,
+            c: 5 * pybamm.StateVectorDot(slice(0, 15)) - 1
+        }
+        with self.assertRaisesRegex(
+            pybamm.ModelError,
+            "time derivative of state vector found",
+        ):
+            model.check_well_posedness(post_discretisation=True)
+
+        # model must be in semi-explicit form
+        model = pybamm.BaseModel()
+        model.rhs = {c: 5 * pybamm.StateVectorDot(slice(0, 15)) - 1}
+        model.initial_conditions = {c: 1}
+        with self.assertRaisesRegex(
+            pybamm.ModelError,
+            "time derivative of state vector found",
+        ):
+            model.check_well_posedness(post_discretisation=True)
+
     def test_check_well_posedness_initial_boundary_conditions(self):
         # Well-posed model - Dirichlet
         whole_cell = ["negative electrode", "separator", "positive electrode"]

From 088e62589cba7e0612fa8632b56b257b191b5624 Mon Sep 17 00:00:00 2001
From: Martin Robinson <martinjrobins@gmail.com>
Date: Wed, 18 Mar 2020 11:48:21 +0000
Subject: [PATCH 30/30] #858 flake8

---
 tests/unit/test_expression_tree/test_symbol.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/tests/unit/test_expression_tree/test_symbol.py b/tests/unit/test_expression_tree/test_symbol.py
index a8462f2559..7faa4992c4 100644
--- a/tests/unit/test_expression_tree/test_symbol.py
+++ b/tests/unit/test_expression_tree/test_symbol.py
@@ -172,7 +172,7 @@ def test_evaluate_ignoring_errors(self):
         self.assertIsNone(pybamm.t.evaluate_ignoring_errors(t=None))
         self.assertEqual(pybamm.t.evaluate_ignoring_errors(t=0), 0)
         self.assertIsNone(pybamm.Parameter("a").evaluate_ignoring_errors())
-        self.assertIsNone(pybamm.StateVector(slice(0,1)).evaluate_ignoring_errors())
+        self.assertIsNone(pybamm.StateVector(slice(0, 1)).evaluate_ignoring_errors())
         self.assertEqual(pybamm.InputParameter("a").evaluate_ignoring_errors(), 1)
 
     def test_symbol_is_constant(self):