Merge branch 'develop' of github.com:pybamm-team/PyBaMM into issue-93…

…3-remove-reactions-dict

CHANGELOG.md

@@ -3,6 +3,8 @@
 ## Features
 
 -   Added functionality to solver to automatically discretise a 0D model ([#947](https://github.com/pybamm-team/PyBaMM/pull/947))
+-   Added sensitivity to `CasadiAlgebraicSolver` ([#940](https://github.com/pybamm-team/PyBaMM/pull/940))
+-   Added `ProcessedSymbolicVariable` class, which can handle symbolic variables (i.e. variables for which the inputs are symbolic) ([#940](https://github.com/pybamm-team/PyBaMM/pull/940))
 -   Made `QuickPlot` compatible with Google Colab ([#935](https://github.com/pybamm-team/PyBaMM/pull/935))
 -   Added `BasicFull` model for lead-acid ([#932](https://github.com/pybamm-team/PyBaMM/pull/932))
 

docs/index.rst

@@ -33,7 +33,6 @@ Contents
     source/experiments/index
     source/simulation
     source/quick_plot
-    source/processed_variable
     source/util
     source/citations
     source/parameters_cli

docs/source/solvers/index.rst

@@ -10,3 +10,5 @@ Solvers
   casadi_solver
   algebraic_solvers
   solution
+  processed_variable
+

docs/source/processed_variable.rst → docs/source/solvers/processed_variable.rst

@@ -3,3 +3,6 @@ Post-Process Variables
 
 .. autoclass:: pybamm.ProcessedVariable
   :members:
+
+.. autoclass:: pybamm.ProcessedSymbolicVariable
+  :members:

docs/source/solvers/solution.rst

@@ -1,5 +1,5 @@
-Solution
-========
+Solutions
+=========
 
 .. autoclass:: pybamm._BaseSolution
   :members:

pybamm/__init__.py

@@ -201,6 +201,8 @@ def version(formatted=False):
 # Solver classes
 #
 from .solvers.solution import Solution, _BaseSolution
+from .solvers.processed_variable import ProcessedVariable
+from .solvers.processed_symbolic_variable import ProcessedSymbolicVariable
 from .solvers.base_solver import BaseSolver
 from .solvers.dummy_solver import DummySolver
 from .solvers.algebraic_solver import AlgebraicSolver
@@ -220,7 +222,6 @@ def version(formatted=False):
 #
 # other
 #
-from .processed_variable import ProcessedVariable
 from .quick_plot import QuickPlot, dynamic_plot, ax_min, ax_max
 
 from .simulation import Simulation, load_sim, is_notebook

pybamm/discretisations/discretisation.py

@@ -950,6 +950,14 @@ def _process_symbol(self, symbol):
 
             return new_symbol
 
+        elif isinstance(symbol, pybamm.InputParameter):
+            # Return a new copy of the input parameter, but set the expected size
+            # according to the domain of the input parameter
+            expected_size = self._get_variable_size(symbol)
+            new_input_parameter = symbol.new_copy()
+            new_input_parameter.set_expected_size(expected_size)
+            return new_input_parameter
+
         else:
             # Backup option: return new copy of the object
             try:

pybamm/expression_tree/input_parameter.py

@@ -1,6 +1,7 @@
 #
 # Parameter classes
 #
+import numbers
 import numpy as np
 import pybamm
 
@@ -15,22 +16,32 @@ class InputParameter(pybamm.Symbol):
     ----------
     name : str
         name of the node
-
+    domain : iterable of str, or str
+        list of domains over which the node is valid (empty list indicates the symbol
+        is valid over all domains)
     """
 
-    def __init__(self, name):
-        super().__init__(name)
+    def __init__(self, name, domain=None):
+        # Expected shape defaults to 1
+        self._expected_size = 1
+        super().__init__(name, domain=domain)
 
     def new_copy(self):
         """ See :meth:`pybamm.Symbol.new_copy()`. """
-        return InputParameter(self.name)
+        new_input_parameter = InputParameter(self.name, self.domain)
+        new_input_parameter._expected_size = self._expected_size
+        return new_input_parameter
+
+    def set_expected_size(self, size):
+        "Specify the size that the input parameter should be"
+        self._expected_size = size
 
     def _evaluate_for_shape(self):
         """
         Returns the scalar 'NaN' to represent the shape of a parameter.
         See :meth:`pybamm.Symbol.evaluate_for_shape()`
         """
-        return np.nan
+        return np.nan * np.ones_like(self._expected_size)
 
     def _jac(self, variable):
         """ See :meth:`pybamm.Symbol._jac()`. """
@@ -44,10 +55,26 @@ def _base_evaluate(self, t=None, y=None, y_dot=None, inputs=None):
         if not isinstance(inputs, dict):
             # if the special input "shape test" is passed, just return 1
             if inputs == "shape test":
-                return 1
+                return np.ones_like(self._expected_size)
             raise TypeError("inputs should be a dictionary")
         try:
-            return inputs[self.name]
+            input_eval = inputs[self.name]
         # raise more informative error if can't find name in dict
         except KeyError:
             raise KeyError("Input parameter '{}' not found".format(self.name))
+
+        if isinstance(input_eval, numbers.Number):
+            input_shape = 1
+        else:
+            input_shape = input_eval.shape[0]
+        if input_shape == self._expected_size:
+            return input_eval
+        else:
+            raise ValueError(
+                "Input parameter '{}' was given an object of size '{}'".format(
+                    self.name, input_shape
+                )
+                + " but was expecting an object of size '{}'.".format(
+                    self._expected_size
+                )
+            )

pybamm/expression_tree/operations/jacobian.py

@@ -76,7 +76,9 @@ def _jac(self, symbol, variable):
             jac = symbol._function_jac(children_jacs)
 
         elif isinstance(symbol, pybamm.Concatenation):
-            children_jacs = [child.jac(variable) for child in symbol.cached_children]
+            children_jacs = [
+                self.jac(child, variable) for child in symbol.cached_children
+            ]
             jac = symbol._concatenation_jac(children_jacs)
 
         else:

pybamm/expression_tree/operations/unpack_symbols.py

@@ -22,7 +22,7 @@ def __init__(self, classes_to_find, unpacked_symbols=None):
 
     def unpack_list_of_symbols(self, list_of_symbols):
         """
-        Unpack a list of symbols. See :meth:`EquationUnpacker.unpack()`
+        Unpack a list of symbols. See :meth:`SymbolUnpacker.unpack()`
 
         Parameters
         ----------
@@ -65,7 +65,7 @@ def unpack_symbol(self, symbol):
             return unpacked
 
     def _unpack(self, symbol):
-        """ See :meth:`EquationUnpacker.unpack()`. """
+        """ See :meth:`SymbolUnpacker.unpack()`. """
 
         children = symbol.children
 
@@ -86,4 +86,3 @@ def _unpack(self, symbol):
                 child_vars = self.unpack_symbol(child)
                 found_vars.update(child_vars)
             return found_vars
-

pybamm/expression_tree/symbol.py

@@ -660,7 +660,9 @@ def evaluates_to_number(self):
         """
         result = self.evaluate_ignoring_errors()
 
-        if isinstance(result, numbers.Number):
+        if isinstance(result, numbers.Number) or (
+            isinstance(result, np.ndarray) and result.shape == ()
+        ):
             return True
         else:
             return False

pybamm/models/base_model.py

@@ -115,6 +115,7 @@ def __init__(self, name="Unnamed model"):
         self._jacobian = None
         self._jacobian_algebraic = None
         self.external_variables = []
+        self._input_parameters = None
 
         # Default behaviour is to use the jacobian and simplify
         self.use_jacobian = True
@@ -320,6 +321,25 @@ def timescale(self, value):
         "Set the timescale"
         self._timescale = value
 
+    @property
+    def input_parameters(self):
+        "Returns all the input parameters in the model"
+        if self._input_parameters is None:
+            self._input_parameters = self._find_input_parameters()
+        return self._input_parameters
+
+    def _find_input_parameters(self):
+        "Find all the input parameters in the model"
+        unpacker = pybamm.SymbolUnpacker(pybamm.InputParameter)
+        all_input_parameters = unpacker.unpack_list_of_symbols(
+            list(self.rhs.values())
+            + list(self.algebraic.values())
+            + list(self.initial_conditions.values())
+            + list(self.variables.values())
+            + [event.expression for event in self.events]
+        )
+        return list(all_input_parameters.values())
+
     def __getitem__(self, key):
         return self.rhs[key]
 

pybamm/solvers/base_solver.py

@@ -183,15 +183,30 @@ def set_up(self, model, inputs=None):
             # set up Jacobian object, for re-use of dict
             jacobian = pybamm.Jacobian()
         else:
+            # Create placeholder inputs for evaluating rhs and algebraic sizes
+            placeholder_inputs = {}
+            for k, v in inputs.items():
+                if isinstance(v, casadi.MX):
+                    placeholder_inputs[k] = np.zeros(v.shape[0])
+                else:
+                    placeholder_inputs[k] = v
             # Convert model attributes to casadi
             t_casadi = casadi.MX.sym("t")
             y_diff = casadi.MX.sym(
                 "y_diff",
-                len(model.concatenated_rhs.evaluate(0, model.y0, inputs=inputs)),
+                len(
+                    model.concatenated_rhs.evaluate(
+                        0, model.y0, inputs=placeholder_inputs
+                    )
+                ),
             )
             y_alg = casadi.MX.sym(
                 "y_alg",
-                len(model.concatenated_algebraic.evaluate(0, model.y0, inputs=inputs)),
+                len(
+                    model.concatenated_algebraic.evaluate(
+                        0, model.y0, inputs=placeholder_inputs
+                    )
+                ),
             )
             y_casadi = casadi.vertcat(y_diff, y_alg)
             p_casadi = {}
@@ -575,7 +590,8 @@ def solve(self, model, t_eval=None, external_variables=None, inputs=None):
         if (np.diff(t_eval) < 0).any():
             raise pybamm.SolverError("t_eval must increase monotonically")
 
-        # Non-dimensionalise t_eval
+        # Set up external variables and inputs
+        ext_and_inputs = self._set_up_ext_and_inputs(model, external_variables, inputs)
 
         # Make sure t_eval is monotonic
         if (np.diff(t_eval) < 0).any():
@@ -584,11 +600,6 @@ def solve(self, model, t_eval=None, external_variables=None, inputs=None):
         # Set up
         timer = pybamm.Timer()
 
-        # Set up external variables and inputs
-        external_variables = external_variables or {}
-        inputs = inputs or {}
-        ext_and_inputs = {**external_variables, **inputs}
-
         # Set up (if not done already)
         if model not in self.models_set_up:
             self.set_up(model, ext_and_inputs)
@@ -699,13 +710,6 @@ def solve(self, model, t_eval=None, external_variables=None, inputs=None):
         solution.set_up_time = set_up_time
         solution.solve_time = timer.time()
 
-        # Add model and inputs to solution
-        solution.model = model
-        solution.inputs = inputs
-
-        # Identify the event that caused termination
-        termination = self.get_termination_reason(solution, model.events)
-
         # Add model and inputs to solution
         solution.model = model
         solution.inputs = ext_and_inputs
@@ -880,6 +884,27 @@ def get_termination_reason(self, solution, events):
             solution.termination = "event: {}".format(termination_event)
             return "the termination event '{}' occurred".format(termination_event)
 
+    def _set_up_ext_and_inputs(self, model, external_variables, inputs):
+        "Set up external variables and input parameters"
+        inputs = inputs or {}
+
+        # Go through all input parameters that can be found in the model
+        # If any of them are *not* provided by "inputs", a symbolic input parameter is
+        # created, with appropriate size
+        for input_param in model.input_parameters:
+            name = input_param.name
+            if name not in inputs:
+                # Only allow symbolic inputs for CasadiAlgebraicSolver
+                if not isinstance(self, pybamm.CasadiAlgebraicSolver):
+                    raise pybamm.SolverError(
+                        "Only CasadiAlgebraicSolver can have symbolic inputs"
+                    )
+                inputs[name] = casadi.MX.sym(name, input_param._expected_size)
+
+        external_variables = external_variables or {}
+        ext_and_inputs = {**external_variables, **inputs}
+        return ext_and_inputs
+
 
 class SolverCallable:
     "A class that will be called by the solver when integrating"

pybamm/solvers/casadi_algebraic_solver.py

@@ -46,13 +46,20 @@ def _integrate(self, model, t_eval, inputs=None):
         t_eval : :class:`numpy.array`, size (k,)
             The times at which to compute the solution
         inputs : dict, optional
-            Any input parameters to pass to the model when solving
+            Any input parameters to pass to the model when solving. If any input
+            parameters that are present in the model are missing from "inputs", then
+            the solution will consist of `ProcessedSymbolicVariable` objects, which must
+            be provided with inputs to obtain their value.
         """
         y0 = model.y0
 
-        y = np.empty((len(y0), len(t_eval)))
+        y = None
 
         # Set up
+        # Record whether there are any symbolic inputs
+        has_symbolic_inputs = any(isinstance(v, casadi.MX) for v in inputs.values())
+
+        # Create casadi objects for the root-finder
         inputs = casadi.vertcat(*[x for x in inputs.values()])
         t_sym = casadi.MX.sym("t")
         y_sym = casadi.MX.sym("y_alg", y0.shape[0])
@@ -71,17 +78,23 @@ def _integrate(self, model, t_eval, inputs=None):
         for idx, t in enumerate(t_eval):
             # Evaluate algebraic with new t and previous y0, if it's already close
             # enough then keep it
-            if np.all(abs(model.algebraic_eval(t, y0, inputs)) < self.tol):
+            # We can't do this if there are symbolic inputs
+            if has_symbolic_inputs is False and np.all(
+                abs(model.casadi_algebraic(t, y0, inputs).full()) < self.tol
+            ):
                 pybamm.logger.debug(
                     "Keeping same solution at t={}".format(t * model.timescale_eval)
                 )
-                y[:, idx] = y0
-            # Otherwise calculate new y0
+                if y is None:
+                    y = y0
+                else:
+                    y = casadi.horzcat(y, y0)
+            # Otherwise calculate new y_sol
             else:
                 t_inputs = casadi.vertcat(t, inputs)
                 # Solve
                 try:
-                    y_sol = roots(y0, t_inputs).full().flatten()
+                    y_sol = roots(y0, t_inputs)
                     success = True
                     message = None
                     # Check final output
@@ -91,11 +104,18 @@ def _integrate(self, model, t_eval, inputs=None):
                     message = err.args[0]
                     fun = None
 
-                if success and np.all(casadi.fabs(fun) < self.tol):
+                # If there are no symbolic inputs, check the function is below the tol
+                # Skip this check if there are symbolic inputs
+                if success and (
+                    has_symbolic_inputs is True or np.all(casadi.fabs(fun) < self.tol)
+                ):
                     # update initial guess for the next iteration
                     y0 = y_sol
                     # update solution array
-                    y[:, idx] = y_sol
+                    if y is None:
+                        y = y_sol
+                    else:
+                        y = casadi.horzcat(y, y_sol)
                 elif not success:
                     raise pybamm.SolverError(
                         "Could not find acceptable solution: {}".format(message)