Repository: pybamm-team/PyBaMM

Issue:
DAE Multiple inputs with fewer sensitivites
- PyBaMM version - 21.12:
- Python version - 3.8.10:

**Describe the bug**
An error arises within the algebraic solver if the number of sensitivities is different to the number of inputs. 

Error is: 
``Slice (start=0, stop=3783, step=1) out of bounds with supplied length of 2722``

I believe as a result of line 86 of casadi_algebraic_solver.py:
                 ```len_rhs = model.len_rhs * (inputs.shape[0] + 1)```
which calls inputs.shape[0]. Since there are more inputs than sensitivities, this results in trying to extract too many elements of rhs because less sensitivity equations have been added.

**To Reproduce**
```
import pybamm
import numpy as np

chemistry = pybamm.parameter_sets.Marquis2019
parameter_values = pybamm.ParameterValues(chemistry=chemistry)
parameter_values["Negative electrode diffusivity [m2.s-1]"] = "[input]"
parameter_values["Negative electrode Bruggeman coefficient (electrolyte)"] = "[input]"

model = pybamm.lithium_ion.DFN()
sim = pybamm.Simulation(model, parameter_values=parameter_values)

inputs = {
    "Negative electrode diffusivity [m2.s-1]": 1e-13,
    "Negative electrode Bruggeman coefficient (electrolyte)": 1.5,
}

sensitivities = ["Negative electrode diffusivity [m2.s-1]"]

t_eval = np.linspace(0, 30, 1000)
sim.solve(t_eval=t_eval, inputs=inputs, calculate_sensitivities=sensitivities)
```


**Expected behaviour**
Solving should not through an error when the sensitivities are fewer than the number of inputs. The length of the extracted rhs in casadi_algebraic_solver.py should depend upon the number of sensitivities and not the number of inputs. 



Discussion:

Patch:
diff --git a/CHANGELOG.md b/CHANGELOG.md
index b75419b52..6950254ed 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,3 +1,4 @@
+
 # [Unreleased](https://github.com/pybamm-team/PyBaMM/)
 
 ## Features
@@ -7,6 +8,7 @@
 
 ## Bug fixes
 
+-   Fix bug where sensitivity calculation failed if len of `calculate_sensitivities` was less than `inputs` ([#1897](https://github.com/pybamm-team/PyBaMM/pull/1897))
 -   Fixed a bug in the eSOH variable calculation when OCV is given as data ([#1975](https://github.com/pybamm-team/PyBaMM/pull/1975))
 -   Fixed a bug where isothermal models did not compute any heat source terms ([#1958](https://github.com/pybamm-team/PyBaMM/pull/1958))
 
@@ -50,6 +52,7 @@
 -   Parameters can now be imported from any given path in `Windows` ([#1900](https://github.com/pybamm-team/PyBaMM/pull/1900))
 -   Fixed initial conditions for the EC SEI model ([#1895](https://github.com/pybamm-team/PyBaMM/pull/1895))
 -   Fixed issue in extraction of sensitivites ([#1894](https://github.com/pybamm-team/PyBaMM/pull/1894))
+>>>>>>> develop
 
 # [v21.12](https://github.com/pybamm-team/PyBaMM/tree/v21.11) - 2021-12-29
 
diff --git a/pybamm/solvers/base_solver.py b/pybamm/solvers/base_solver.py
index 6bc05e129..a8f9d9cb7 100644
--- a/pybamm/solvers/base_solver.py
+++ b/pybamm/solvers/base_solver.py
@@ -248,7 +248,32 @@ class BaseSolver(object):
                 y_and_S = y_casadi
 
         # if we will change the equations to include the explicit sensitivity
-        # equations, then we also need to update the mass matrix and bounds
+        # equations, then we also need to update the mass matrix and bounds.
+        # First, we reset the mass matrix and bounds back to their original form
+        # if they have been extended
+        if model.bounds[0].shape[0] > model.len_rhs_and_alg:
+            model.bounds = (
+                model.bounds[0][: model.len_rhs_and_alg],
+                model.bounds[1][: model.len_rhs_and_alg],
+            )
+        if (
+            model.mass_matrix is not None
+            and model.mass_matrix.shape[0] > model.len_rhs_and_alg
+        ):
+            if model.mass_matrix_inv is not None:
+                model.mass_matrix_inv = pybamm.Matrix(
+                    model.mass_matrix_inv.entries[
+                        : model.len_rhs, : model.len_rhs
+                    ]
+                )
+            model.mass_matrix = pybamm.Matrix(
+                model.mass_matrix.entries[
+                    : model.len_rhs_and_alg, : model.len_rhs_and_alg
+                ]
+            )
+
+        # now we can extend them by the number of sensitivity parameters
+        # if needed
         if calculate_sensitivities_explicit:
             if model.len_rhs != 0:
                 n_inputs = model.len_rhs_sens // model.len_rhs
@@ -276,28 +301,6 @@ class BaseSolver(object):
                         [model.mass_matrix.entries] * (n_inputs + 1), format="csr"
                     )
                 )
-        else:
-            # take care if calculate_sensitivites used then not used
-            if model.bounds[0].shape[0] > model.len_rhs_and_alg:
-                model.bounds = (
-                    model.bounds[0][: model.len_rhs_and_alg],
-                    model.bounds[1][: model.len_rhs_and_alg],
-                )
-            if (
-                model.mass_matrix is not None
-                and model.mass_matrix.shape[0] > model.len_rhs_and_alg
-            ):
-                if model.mass_matrix_inv is not None:
-                    model.mass_matrix_inv = pybamm.Matrix(
-                        model.mass_matrix_inv.entries[
-                            : model.len_rhs, : model.len_rhs
-                        ]
-                    )
-                model.mass_matrix = pybamm.Matrix(
-                    model.mass_matrix.entries[
-                        : model.len_rhs_and_alg, : model.len_rhs_and_alg
-                    ]
-                )
 
         def process(symbol, name, use_jacobian=None):
             def report(string):
diff --git a/pybamm/solvers/casadi_algebraic_solver.py b/pybamm/solvers/casadi_algebraic_solver.py
index 5125c7c4e..052df9e80 100644
--- a/pybamm/solvers/casadi_algebraic_solver.py
+++ b/pybamm/solvers/casadi_algebraic_solver.py
@@ -83,7 +83,7 @@ class CasadiAlgebraicSolver(pybamm.BaseSolver):
             if model.len_rhs_and_alg == y0.shape[0]:
                 len_rhs = model.len_rhs
             else:
-                len_rhs = model.len_rhs * (inputs.shape[0] + 1)
+                len_rhs = model.len_rhs + model.len_rhs_sens
             y0_diff = y0[:len_rhs]
             y0_alg = y0[len_rhs:]
 


Test patch:
diff --git a/tests/unit/test_solvers/test_casadi_solver.py b/tests/unit/test_solvers/test_casadi_solver.py
index 1b4444511..55d2f070e 100644
--- a/tests/unit/test_solvers/test_casadi_solver.py
+++ b/tests/unit/test_solvers/test_casadi_solver.py
@@ -865,6 +865,64 @@ class TestCasadiSolverODEsWithForwardSensitivityEquations(unittest.TestCase):
             solution["var squared"].data, np.exp(0.1 * solution.t) ** 2
         )
 
+    def test_solve_sensitivity_subset(self):
+        # Create model
+        model = pybamm.BaseModel()
+        var = pybamm.Variable("var")
+        p = pybamm.InputParameter("p")
+        q = pybamm.InputParameter("q")
+        r = pybamm.InputParameter("r")
+        model.rhs = {var: p * q}
+        model.initial_conditions = {var: r}
+
+        # only calculate the sensitivities of a subset of parameters
+        solver = pybamm.CasadiSolver(rtol=1e-10, atol=1e-10)
+        t_eval = np.linspace(0, 1, 80)
+        solution = solver.solve(
+            model,
+            t_eval,
+            inputs={"q": 2, "r": -1, "p": 0.1},
+            calculate_sensitivities=["q", "p"],
+        )
+        np.testing.assert_allclose(solution.y[0], -1 + 0.2 * solution.t)
+        np.testing.assert_allclose(
+            solution.sensitivities["p"],
+            (2 * solution.t)[:, np.newaxis],
+        )
+        np.testing.assert_allclose(
+            solution.sensitivities["q"],
+            (0.1 * solution.t)[:, np.newaxis],
+        )
+        self.assertTrue("r" not in solution.sensitivities)
+        np.testing.assert_allclose(
+            solution.sensitivities["all"],
+            np.hstack(
+                [
+                    solution.sensitivities["p"],
+                    solution.sensitivities["q"],
+                ]
+            ),
+        )
+
+        solution = solver.solve(
+            model,
+            t_eval,
+            inputs={"q": 2, "r": -1, "p": 0.1},
+            calculate_sensitivities=["r"],
+        )
+        np.testing.assert_allclose(solution.y[0], -1 + 0.2 * solution.t)
+        self.assertTrue("p" not in solution.sensitivities)
+        self.assertTrue("q" not in solution.sensitivities)
+        np.testing.assert_allclose(solution.sensitivities["r"], 1)
+        np.testing.assert_allclose(
+            solution.sensitivities["all"],
+            np.hstack(
+                [
+                    solution.sensitivities["r"],
+                ]
+            ),
+        )
+
 
 class TestCasadiSolverDAEsWithForwardSensitivityEquations(unittest.TestCase):
     def test_solve_sensitivity_scalar_var_scalar_input(self):
@@ -939,6 +997,49 @@ class TestCasadiSolverDAEsWithForwardSensitivityEquations(unittest.TestCase):
             atol=1e-7,
         )
 
+    def test_solve_sensitivity_subset(self):
+        # Create model
+        model = pybamm.BaseModel()
+        var = pybamm.Variable("var")
+        var2 = pybamm.Variable("var2")
+        p = pybamm.InputParameter("p")
+        q = pybamm.InputParameter("q")
+        r = pybamm.InputParameter("r")
+        model.rhs = {var: p * q}
+        model.algebraic = {var2: 2 * var - var2}
+        model.initial_conditions = {var: r, var2: 2 * r}
+
+        # only calculate the sensitivities of a subset of parameters
+        solver = pybamm.CasadiSolver(rtol=1e-10, atol=1e-10)
+        t_eval = np.linspace(0, 1, 80)
+        solution = solver.solve(
+            model,
+            t_eval,
+            inputs={"p": 0.1, "q": 2, "r": -1, "s": 0.5},
+            calculate_sensitivities=["p", "q"],
+        )
+        np.testing.assert_allclose(solution.y[0], -1 + 0.2 * solution.t)
+        np.testing.assert_allclose(solution.y[-1], 2 * (-1 + 0.2 * solution.t))
+        np.testing.assert_allclose(
+            solution.sensitivities["p"][::2],
+            (2 * solution.t)[:, np.newaxis],
+        )
+        np.testing.assert_allclose(
+            solution.sensitivities["q"][::2],
+            (0.1 * solution.t)[:, np.newaxis],
+        )
+        self.assertTrue("r" not in solution.sensitivities)
+        self.assertTrue("s" not in solution.sensitivities)
+        np.testing.assert_allclose(
+            solution.sensitivities["all"],
+            np.hstack(
+                [
+                    solution.sensitivities["p"],
+                    solution.sensitivities["q"],
+                ]
+            ),
+        )
+
 
 if __name__ == "__main__":
     print("Add -v for more debug output")