Repository: pybamm-team/PyBaMM

Issue:
Replace `interp2d` in `ProcessedVariable`
We currently use `interp2d` when creating processed variables, but this is deprecated. We should switch to `RegularGridInterpolator`.

```
DeprecationWarning: `interp2d` is deprecated!
`interp2d` is deprecated in SciPy 1.10 and will be removed in SciPy 1.12.0.

For legacy code, nearly bug-for-bug compatible replacements are
`RectBivariateSpline` on regular grids, and `bisplrep`/`bisplev` for
scattered 2D data.

In new code, for regular grids use `RegularGridInterpolator` instead.
For scattered data, prefer `LinearNDInterpolator` or
`CloughTocher2DInterpolator`.
```

Discussion:

Patch:
diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml
index 2fab45486..eb2904759 100644
--- a/.pre-commit-config.yaml
+++ b/.pre-commit-config.yaml
@@ -4,12 +4,12 @@ ci:
 
 repos:
   - repo: https://github.com/psf/black
-    rev: 23.3.0
+    rev: 23.1.0
     hooks:
       - id: black
 
   - repo: https://github.com/charliermarsh/ruff-pre-commit
-    rev: "v0.0.260"
+    rev: "v0.0.259"
     hooks:
       - id: ruff
         args: [--ignore=E741, --exclude=__init__.py]
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 27b39d953..482f9cd86 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,17 +1,5 @@
 # [Unreleased](https://github.com/pybamm-team/PyBaMM/)
 
-# Features
-
-- Updated to casadi 3.6, which required some changes to the casadi integrator. ([#2859](https://github.com/pybamm-team/PyBaMM/pull/2859))
-
-## Bug fixes
-
-- Fixed a bug in the discretisation of initial conditions of a scaled variable ([#2856](https://github.com/pybamm-team/PyBaMM/pull/2856))
-
-# Breaking changes
-
-- Made `Jupyter` a development only dependency. Now `Jupyter` would not be a required dependency for users while installing `PyBaMM`. ([#2846](https://github.com/pybamm-team/PyBaMM/pull/2846))
-
 # [v23.3](https://github.com/pybamm-team/PyBaMM/tree/v23.3) - 2023-03-31
 
 ## Features
diff --git a/pybamm/solvers/casadi_solver.py b/pybamm/solvers/casadi_solver.py
index 9383c2d2f..f7a981f19 100644
--- a/pybamm/solvers/casadi_solver.py
+++ b/pybamm/solvers/casadi_solver.py
@@ -516,11 +516,9 @@ class CasadiSolver(pybamm.BaseSolver):
                 if t_eval_shifted_rounded in self.integrators[model]:
                     return self.integrators[model][t_eval_shifted_rounded]
                 else:
-                    method, problem, options, time_args = self.integrator_specs[model]
-                    time_args = [t_eval_shifted[0], t_eval_shifted[1:]]
-                    integrator = casadi.integrator(
-                        "F", method, problem, *time_args, options
-                    )
+                    method, problem, options = self.integrator_specs[model]
+                    options["grid"] = t_eval_shifted
+                    integrator = casadi.integrator("F", method, problem, options)
                     self.integrators[model][t_eval_shifted_rounded] = integrator
                     return integrator
         else:
@@ -544,7 +542,6 @@ class CasadiSolver(pybamm.BaseSolver):
             y_full = casadi.vertcat(y_diff, y_alg)
 
             if use_grid is False:
-                time_args = []
                 # rescale time
                 t_min = casadi.MX.sym("t_min")
                 t_max = casadi.MX.sym("t_max")
@@ -553,7 +550,7 @@ class CasadiSolver(pybamm.BaseSolver):
                 # add time limits as inputs
                 p_with_tlims = casadi.vertcat(p, t_min, t_max)
             else:
-                time_args = [t_eval_shifted[0], t_eval_shifted[1:]]
+                options.update({"grid": t_eval_shifted, "output_t0": True})
                 # rescale time
                 t_min = casadi.MX.sym("t_min")
                 # Set dummy parameters for consistency with rescaled time
@@ -586,8 +583,8 @@ class CasadiSolver(pybamm.BaseSolver):
                         "alg": algebraic(t_scaled, y_full, p),
                     }
                 )
-            integrator = casadi.integrator("F", method, problem, *time_args, options)
-            self.integrator_specs[model] = method, problem, options, time_args
+            integrator = casadi.integrator("F", method, problem, options)
+            self.integrator_specs[model] = method, problem, options
             if use_grid is False:
                 self.integrators[model] = {"no grid": integrator}
             else:
@@ -658,15 +655,13 @@ class CasadiSolver(pybamm.BaseSolver):
             len_alg = len_alg * (num_parameters + 1)
 
         y0_diff = y0[:len_rhs]
-        y0_alg_exact = y0[len_rhs:]
+        y0_alg = y0[len_rhs:]
         if self.perturb_algebraic_initial_conditions and len_alg > 0:
             # Add a tiny perturbation to the algebraic initial conditions
             # For some reason this helps with convergence
             # The actual value of the initial conditions for the algebraic variables
             # doesn't matter
-            y0_alg = y0_alg_exact * (1 + 1e-6 * casadi.DM(np.random.rand(len_alg)))
-        else:
-            y0_alg = y0_alg_exact
+            y0_alg = y0_alg * (1 + 1e-6 * casadi.DM(np.random.rand(len_alg)))
         pybamm.logger.spam("Finished preliminary setup for integrator run")
 
         # Solve
@@ -687,13 +682,7 @@ class CasadiSolver(pybamm.BaseSolver):
                 raise pybamm.SolverError(error.args[0])
             pybamm.logger.debug("Finished casadi integrator")
             integration_time = timer.time()
-            # Manually add initial conditions and concatenate
-            x_sol = casadi.horzcat(y0_diff, casadi_sol["xf"])
-            if len_alg > 0:
-                z_sol = casadi.horzcat(y0_alg_exact, casadi_sol["zf"])
-                y_sol = casadi.vertcat(x_sol, z_sol)
-            else:
-                y_sol = x_sol
+            y_sol = casadi.vertcat(casadi_sol["xf"], casadi_sol["zf"])
             sol = pybamm.Solution(
                 t_eval,
                 y_sol,
@@ -707,7 +696,7 @@ class CasadiSolver(pybamm.BaseSolver):
         else:
             # Repeated calls to the integrator
             x = y0_diff
-            z = y0_alg_exact
+            z = y0_alg
             y_diff = x
             y_alg = z
             for i in range(len(t_eval) - 1):
diff --git a/pybamm/solvers/processed_variable.py b/pybamm/solvers/processed_variable.py
index 610bf92ff..80b2be568 100644
--- a/pybamm/solvers/processed_variable.py
+++ b/pybamm/solvers/processed_variable.py
@@ -219,11 +219,10 @@ class ProcessedVariable(object):
         else:
             # function of space and time. Note that the order of 't' and 'space'
             # is the reverse of what you'd expect
-            self._interpolation_function = interp.interp2d(
-                self.t_pts,
-                pts_for_interp,
-                entries_for_interp,
-                kind="linear",
+            self._interpolation_function = interp.RegularGridInterpolator(
+                (self.t_pts, pts_for_interp),
+                entries_for_interp.T,
+                method="linear",
                 fill_value=np.nan,
                 bounds_error=False,
             )
@@ -595,11 +594,10 @@ class Interpolant2D:
     def __init__(
         self, first_dim_pts_for_interp, second_dim_pts_for_interp, entries_for_interp
     ):
-        self.interpolant = interp.interp2d(
-            second_dim_pts_for_interp,
-            first_dim_pts_for_interp,
-            entries_for_interp[:, :, 0],
-            kind="linear",
+        self.interpolant = interp.RegularGridInterpolator(
+            (second_dim_pts_for_interp, first_dim_pts_for_interp),
+            entries_for_interp.T[:, :, 0],
+            method="linear",
             fill_value=np.nan,
             bounds_error=False,
         )
diff --git a/pybamm/spatial_methods/spatial_method.py b/pybamm/spatial_methods/spatial_method.py
index acb0227bc..a7e5dfc16 100644
--- a/pybamm/spatial_methods/spatial_method.py
+++ b/pybamm/spatial_methods/spatial_method.py
@@ -134,7 +134,7 @@ class SpatialMethod:
             matrix = vstack([identity for _ in range(tertiary_domain_size)])
             out = pybamm.Matrix(matrix) @ symbol
         elif broadcast_type.startswith("full"):
-            out = symbol * pybamm.Vector(np.ones(full_domain_size), domains=domains)
+            out = symbol * pybamm.Vector(np.ones(full_domain_size))
 
         out.domains = domains.copy()
         return out
diff --git a/requirements.txt b/requirements.txt
index 665b77493..b98b6576c 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -6,6 +6,7 @@ autograd >= 1.2
 scikit-fem >= 0.2.0
 casadi >= 3.5.0
 imageio>=2.9.0
+jupyter  # For example notebooks
 pybtex>=0.24.0
 sympy >= 1.8
 bpx
diff --git a/setup.py b/setup.py
index 5d4079d06..f5aa1638a 100644
--- a/setup.py
+++ b/setup.py
@@ -197,7 +197,7 @@ setup(
         "scikit-fem>=0.2.0",
         "casadi>=3.5.0",
         "imageio>=2.9.0",
-        "importlib-metadata",
+        "jupyter",  # For example notebooks
         "pybtex>=0.24.0",
         "sympy>=1.8",
         "bpx",
@@ -218,7 +218,6 @@ setup(
         "dev": [
             "pre-commit",  # For code style checking
             "black",  # For code style auto-formatting
-            "jupyter",  # For example notebooks
         ],
     },
     entry_points={
diff --git a/tox.ini b/tox.ini
index f47e8fa1e..555a7f073 100644
--- a/tox.ini
+++ b/tox.ini
@@ -18,7 +18,6 @@ deps =
      doctests: sphinx-copybutton
      doctests: myst-parser
      !windows-!mac: scikits.odes
-     examples: jupyter   # For example notebooks
 
 commands =
      tests-!windows-!mac: sh -c "pybamm_install_jax"  # install jax, jaxlib for ubuntu


Test patch:
diff --git a/tests/integration/test_models/test_full_battery_models/test_lithium_ion/base_lithium_ion_tests.py b/tests/integration/test_models/test_full_battery_models/test_lithium_ion/base_lithium_ion_tests.py
index cfcd7c0d3..2ea8fa6bb 100644
--- a/tests/integration/test_models/test_full_battery_models/test_lithium_ion/base_lithium_ion_tests.py
+++ b/tests/integration/test_models/test_full_battery_models/test_lithium_ion/base_lithium_ion_tests.py
@@ -21,7 +21,10 @@ class BaseIntegrationTestLithiumIon:
         model = self.model()
         param = pybamm.ParameterValues("Ecker2015")
         modeltest = tests.StandardModelTest(model, parameter_values=param)
-        modeltest.test_sensitivities("Current function [A]", 0.15652)
+        modeltest.test_sensitivities(
+            "Current function [A]",
+            0.15652,
+        )
 
     def test_basic_processing_1plus1D(self):
         options = {"current collector": "potential pair", "dimensionality": 1}
diff --git a/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_dfn.py b/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_dfn.py
index c732234d4..a0bd97186 100644
--- a/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_dfn.py
+++ b/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_dfn.py
@@ -12,10 +12,6 @@ class TestDFN(BaseIntegrationTestLithiumIon, unittest.TestCase):
     def setUp(self):
         self.model = pybamm.lithium_ion.DFN
 
-    def test_sensitivities(self):
-        # skip sensitivities test for now as it is failing with casadi 3.6
-        pass
-
     def test_particle_distribution_in_x(self):
         model = pybamm.lithium_ion.DFN()
         param = model.default_parameter_values
diff --git a/tests/unit/test_solvers/test_base_solver.py b/tests/unit/test_solvers/test_base_solver.py
index e9c8c4f89..9c4e09d72 100644
--- a/tests/unit/test_solvers/test_base_solver.py
+++ b/tests/unit/test_solvers/test_base_solver.py
@@ -228,7 +228,7 @@ class TestBaseSolver(unittest.TestCase):
         # with casadi
         solver = pybamm.BaseSolver(root_method="casadi")
         with self.assertRaisesRegex(
-            pybamm.SolverError, "Could not find acceptable solution: Error in Function"
+            pybamm.SolverError, "Could not find acceptable solution: .../casadi"
         ):
             solver.calculate_consistent_state(Model())
 
diff --git a/tests/unit/test_solvers/test_casadi_algebraic_solver.py b/tests/unit/test_solvers/test_casadi_algebraic_solver.py
index fcba45382..fe4c393f0 100644
--- a/tests/unit/test_solvers/test_casadi_algebraic_solver.py
+++ b/tests/unit/test_solvers/test_casadi_algebraic_solver.py
@@ -69,7 +69,7 @@ class TestCasadiAlgebraicSolver(unittest.TestCase):
 
         solver = pybamm.CasadiAlgebraicSolver()
         with self.assertRaisesRegex(
-            pybamm.SolverError, "Could not find acceptable solution: Error in Function"
+            pybamm.SolverError, "Could not find acceptable solution: .../casadi"
         ):
             solver._integrate(model, np.array([0]), {})
         solver = pybamm.CasadiAlgebraicSolver(extra_options={"error_on_fail": False})
diff --git a/tests/unit/test_spatial_methods/test_finite_volume/test_finite_volume.py b/tests/unit/test_spatial_methods/test_finite_volume/test_finite_volume.py
index 90bbc2905..29d6eedd7 100644
--- a/tests/unit/test_spatial_methods/test_finite_volume/test_finite_volume.py
+++ b/tests/unit/test_spatial_methods/test_finite_volume/test_finite_volume.py
@@ -529,27 +529,6 @@ class TestFiniteVolume(unittest.TestCase):
         self.assertEqual(disc.bcs[var]["right"][0], pybamm.Scalar(0))
         self.assertEqual(disc.bcs[var]["right"][1], "Neumann")
 
-    def test_full_broadcast_domains(self):
-        model = pybamm.BaseModel()
-        var = pybamm.Variable(
-            "var", domain=["negative electrode", "separator"], scale=100
-        )
-        model.rhs = {var: 0}
-        a = pybamm.InputParameter("a")
-        ic = pybamm.concatenation(
-            pybamm.FullBroadcast(a * 100, "negative electrode"),
-            pybamm.FullBroadcast(100, "separator"),
-        )
-        model.initial_conditions = {var: ic}
-
-        mesh = get_mesh_for_testing()
-        spatial_methods = {
-            "negative electrode": pybamm.FiniteVolume(),
-            "separator": pybamm.FiniteVolume(),
-        }
-        disc = pybamm.Discretisation(mesh, spatial_methods)
-        disc.process_model(model)
-
 
 if __name__ == "__main__":
     print("Add -v for more debug output")