Repository: pybamm-team/PyBaMM

Issue:
Make default interpolator linear
I think we should change the default interpolant to be linear. This is normally the default for all other interpolators and is the "safest" option in the sense that you don't get oscillations or introduce maxima. 

Discussion:

Patch:
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 45146f367..9cefbd0e0 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -6,6 +6,7 @@
 
 ## Breaking changes
 
+-   When creating a `pybamm.Interpolant` the default interpolator is now "linear". Passing data directly to `ParameterValues` using the ``[data]`` tag will be still used to create a cubic spline interpolant, as before ([#2258](https://github.com/pybamm-team/PyBaMM/pull/2258))
 -   Events must now be defined in such a way that they are positive at the initial conditions (events will be triggered when they become negative, instead of when they change sign in either direction) ([#2212](https://github.com/pybamm-team/PyBaMM/pull/2212))
 
 # [v22.8](https://github.com/pybamm-team/PyBaMM/tree/v22.8) - 2022-08-31
diff --git a/docs/tutorials/add-parameter-values.rst b/docs/tutorials/add-parameter-values.rst
index 63625044d..168d3c0e0 100644
--- a/docs/tutorials/add-parameter-values.rst
+++ b/docs/tutorials/add-parameter-values.rst
@@ -102,6 +102,9 @@ called (must be in the same folder), with the tag ``[data]``, for example:
 | Example [m2.s-1]    | [data]diffusivity_AuthorYear     | AuthorYear   | some data   |
 +---------------------+----------------------------------+--------------+-------------+
 
+Data passed using the ``[data]`` tag will be used to create a cubic spline interpolant. For more control over the interpolation you can create your own :class:`Interpolant` from the data and pass that to the :class:`ParameterValues` class instead. 
+
+
 Using new parameters
 --------------------
 
diff --git a/pybamm/expression_tree/interpolant.py b/pybamm/expression_tree/interpolant.py
index 78f0cceab..b59dc43e6 100644
--- a/pybamm/expression_tree/interpolant.py
+++ b/pybamm/expression_tree/interpolant.py
@@ -3,6 +3,7 @@
 #
 import numpy as np
 from scipy import interpolate
+import warnings
 
 import pybamm
 
@@ -24,7 +25,8 @@ class Interpolant(pybamm.Function):
         Name of the interpolant. Default is None, in which case the name "interpolating
         function" is given.
     interpolator : str, optional
-        Which interpolator to use ("linear", "pchip", or "cubic spline").
+        Which interpolator to use. Can be "linear", "cubic", or "pchip". Default is
+        "linear".
     extrapolate : bool, optional
         Whether to extrapolate for points that are outside of the parametrisation
         range, or return NaN (following default behaviour from scipy). Default is True.
@@ -38,20 +40,24 @@ class Interpolant(pybamm.Function):
         y,
         children,
         name=None,
-        interpolator=None,
+        interpolator="linear",
         extrapolate=True,
         entries_string=None,
     ):
+        # "cubic spline" has been renamed to "cubic"
+        if interpolator == "cubic spline":
+            interpolator = "cubic"
+            warnings.warn(
+                "The 'cubic spline' interpolator has been renamed to 'cubic'.",
+                DeprecationWarning,
+            )
 
-        # set default dimension value
-        self.dimension = 1
+        # Check interpolator is valid
+        if interpolator not in ["linear", "cubic", "pchip"]:
+            raise ValueError("interpolator '{}' not recognised".format(interpolator))
 
+        # Perform some checks on the data
         if isinstance(x, (tuple, list)) and len(x) == 2:
-            interpolator = interpolator or "linear"
-            if interpolator != "linear":
-                raise ValueError(
-                    "interpolator should be 'linear' if x is two-dimensional"
-                )
             x1, x2 = x
             if y.ndim != 2:
                 raise ValueError("y should be two-dimensional if len(x)=2")
@@ -66,7 +72,6 @@ class Interpolant(pybamm.Function):
                     f"but x2.shape={x2.shape} and y.shape={y.shape}"
                 )
         else:
-            interpolator = interpolator or "cubic spline"
             if isinstance(x, (tuple, list)):
                 x1 = x[0]
             else:
@@ -78,7 +83,7 @@ class Interpolant(pybamm.Function):
                     "len(x1) should equal y=shape[0], "
                     f"but x1.shape={x1.shape} and y.shape={y.shape}"
                 )
-
+        # children should be a list not a symbol
         if isinstance(children, pybamm.Symbol):
             children = [children]
         # Either a single x is provided and there is one child
@@ -92,32 +97,41 @@ class Interpolant(pybamm.Function):
                 "child should have size 1 if y is two-dimensional and len(x)==1"
             )
 
-        if interpolator == "linear":
-            if len(x) == 1:
-                self.dimension = 1
+        # Create interpolating function
+        if len(x) == 1:
+            self.dimension = 1
+            if interpolator == "linear":
                 if extrapolate is False:
-                    interpolating_function = interpolate.interp1d(
-                        x1, y.T, bounds_error=False, fill_value=np.nan
-                    )
+                    fill_value = np.nan
                 elif extrapolate is True:
-                    interpolating_function = interpolate.interp1d(
-                        x1, y.T, bounds_error=False, fill_value="extrapolate"
-                    )
-            elif len(x) == 2:
-                self.dimension = 2
-                interpolating_function = interpolate.interp2d(x1, x2, y)
+                    fill_value = "extrapolate"
+                interpolating_function = interpolate.interp1d(
+                    x1,
+                    y.T,
+                    bounds_error=False,
+                    fill_value=fill_value,
+                )
+            elif interpolator == "cubic":
+                interpolating_function = interpolate.CubicSpline(
+                    x1, y, extrapolate=extrapolate
+                )
+            elif interpolator == "pchip":
+                interpolating_function = interpolate.PchipInterpolator(
+                    x1, y, extrapolate=extrapolate
+                )
+        elif len(x) == 2:
+            self.dimension = 2
+            if interpolator == "pchip":
+                raise ValueError(
+                    "interpolator should be 'linear' or 'cubic' if x is two-dimensional"
+                )
             else:
-                raise ValueError("Invalid dimension of x: {0}".format(len(x)))
-        elif interpolator == "pchip":
-            interpolating_function = interpolate.PchipInterpolator(
-                x1, y, extrapolate=extrapolate
-            )
-        elif interpolator == "cubic spline":
-            interpolating_function = interpolate.CubicSpline(
-                x1, y, extrapolate=extrapolate
-            )
+                interpolating_function = interpolate.interp2d(
+                    x1, x2, y, kind=interpolator
+                )
         else:
-            raise ValueError("interpolator '{}' not recognised".format(interpolator))
+            raise ValueError("Invalid dimension of x: {0}".format(len(x)))
+
         # Set name
         if name is None:
             name = "interpolating_function"
@@ -127,6 +141,7 @@ class Interpolant(pybamm.Function):
         super().__init__(
             interpolating_function, *children, name=name, derivative="derivative"
         )
+
         # Store information as attributes
         self.interpolator = interpolator
         self.extrapolate = extrapolate
diff --git a/pybamm/expression_tree/operations/convert_to_casadi.py b/pybamm/expression_tree/operations/convert_to_casadi.py
index 7405d7f32..1531ac06b 100644
--- a/pybamm/expression_tree/operations/convert_to_casadi.py
+++ b/pybamm/expression_tree/operations/convert_to_casadi.py
@@ -134,12 +134,12 @@ class CasadiConverter(object):
             elif isinstance(symbol, pybamm.Interpolant):
                 if symbol.interpolator == "linear":
                     solver = "linear"
-                elif symbol.interpolator == "cubic spline":
+                elif symbol.interpolator == "cubic":
                     solver = "bspline"
                 elif symbol.interpolator == "pchip":
                     raise NotImplementedError(
                         "The interpolator 'pchip' is not supported by CasAdi. "
-                        "Use 'linear' or 'cubic spline' instead. "
+                        "Use 'linear' or 'cubic' instead. "
                         "Alternatively, set 'model.convert_to_format = 'python'' "
                         "and use a non-CasADi solver. "
                     )
diff --git a/pybamm/parameters/parameter_values.py b/pybamm/parameters/parameter_values.py
index 6252b80c9..ac7073d4c 100644
--- a/pybamm/parameters/parameter_values.py
+++ b/pybamm/parameters/parameter_values.py
@@ -617,8 +617,14 @@ class ParameterValues:
                     else:
                         input_data = data
 
+                    # For parameters provided as data we use a cubic interpolant
+                    # Note: the cubic interpolant can be differentiated
                     function = pybamm.Interpolant(
-                        input_data[0], input_data[-1], new_children, name=name
+                        input_data[0],
+                        input_data[-1],
+                        new_children,
+                        interpolator="cubic",
+                        name=name,
                     )
                     # Define event to catch extrapolation. In these events the sign is
                     # important: it should be positive inside of the range and negative


Test patch:
diff --git a/tests/unit/test_expression_tree/test_interpolant.py b/tests/unit/test_expression_tree/test_interpolant.py
index 1310714db..327149923 100644
--- a/tests/unit/test_expression_tree/test_interpolant.py
+++ b/tests/unit/test_expression_tree/test_interpolant.py
@@ -35,9 +35,9 @@ class TestInterpolant(unittest.TestCase):
             pybamm.Interpolant(
                 (np.ones(12), np.ones(10)), np.ones((10, 12)), pybamm.Symbol("a")
             )
-        with self.assertRaisesRegex(ValueError, "interpolator should be 'linear'"):
+        with self.assertWarns(DeprecationWarning):
             pybamm.Interpolant(
-                (np.ones(10), np.ones(12)),
+                (np.ones(12), np.ones(10)),
                 np.ones((10, 12)),
                 (pybamm.Symbol("a"), pybamm.Symbol("b")),
                 interpolator="cubic spline",
@@ -47,23 +47,23 @@ class TestInterpolant(unittest.TestCase):
         x = np.linspace(0, 1, 200)
         y = pybamm.StateVector(slice(0, 2))
         # linear
-        for interpolator in ["linear", "pchip", "cubic spline"]:
+        for interpolator in ["linear", "cubic", "pchip"]:
             interp = pybamm.Interpolant(x, 2 * x, y, interpolator=interpolator)
             np.testing.assert_array_almost_equal(
                 interp.evaluate(y=np.array([0.397, 1.5]))[:, 0], np.array([0.794, 3])
             )
         # square
         y = pybamm.StateVector(slice(0, 1))
-        for interpolator in ["linear", "pchip", "cubic spline"]:
-            interp = pybamm.Interpolant(x, x ** 2, y, interpolator=interpolator)
+        for interpolator in ["linear", "cubic", "pchip"]:
+            interp = pybamm.Interpolant(x, x**2, y, interpolator=interpolator)
             np.testing.assert_array_almost_equal(
-                interp.evaluate(y=np.array([0.397]))[:, 0], np.array([0.397 ** 2])
+                interp.evaluate(y=np.array([0.397]))[:, 0], np.array([0.397**2])
             )
 
         # with extrapolation set to False
-        for interpolator in ["linear", "pchip", "cubic spline"]:
+        for interpolator in ["linear", "cubic", "pchip"]:
             interp = pybamm.Interpolant(
-                x, x ** 2, y, interpolator=interpolator, extrapolate=False
+                x, x**2, y, interpolator=interpolator, extrapolate=False
             )
             np.testing.assert_array_equal(
                 interp.evaluate(y=np.array([2]))[:, 0], np.array([np.nan])
@@ -74,7 +74,7 @@ class TestInterpolant(unittest.TestCase):
         y = np.tile(2 * x, (10, 1)).T
         var = pybamm.StateVector(slice(0, 1))
         # linear
-        for interpolator in ["linear", "pchip", "cubic spline"]:
+        for interpolator in ["linear", "cubic", "pchip"]:
             interp = pybamm.Interpolant(x, y, var, interpolator=interpolator)
             np.testing.assert_array_almost_equal(
                 interp.evaluate(y=np.array([0.397])), 0.794 * np.ones((10, 1))
@@ -83,7 +83,7 @@ class TestInterpolant(unittest.TestCase):
     def test_interpolation_2_x_2d_y(self):
         x = (np.arange(-5.01, 5.01, 0.05), np.arange(-5.01, 5.01, 0.01))
         xx, yy = np.meshgrid(x[0], x[1])
-        z = np.sin(xx ** 2 + yy ** 2)
+        z = np.sin(xx**2 + yy**2)
         var1 = pybamm.StateVector(slice(0, 1))
         var2 = pybamm.StateVector(slice(1, 2))
         # linear
@@ -91,6 +91,11 @@ class TestInterpolant(unittest.TestCase):
         np.testing.assert_array_almost_equal(
             interp.evaluate(y=np.array([0, 0])), 0, decimal=3
         )
+        # cubic
+        interp = pybamm.Interpolant(x, z, (var1, var2), interpolator="cubic")
+        np.testing.assert_array_almost_equal(
+            interp.evaluate(y=np.array([0, 0])), 0, decimal=3
+        )
 
     def test_name(self):
         a = pybamm.Symbol("a")
@@ -105,7 +110,7 @@ class TestInterpolant(unittest.TestCase):
         y = pybamm.StateVector(slice(0, 2))
         # linear (derivative should be 2)
         # linear interpolator cannot be differentiated
-        for interpolator in ["pchip", "cubic spline"]:
+        for interpolator in ["cubic", "pchip"]:
             interp_diff = pybamm.Interpolant(
                 x, 2 * x, y, interpolator=interpolator
             ).diff(y)
@@ -113,9 +118,9 @@ class TestInterpolant(unittest.TestCase):
                 interp_diff.evaluate(y=np.array([0.397, 1.5]))[:, 0], np.array([2, 2])
             )
         # square (derivative should be 2*x)
-        for interpolator in ["pchip", "cubic spline"]:
+        for interpolator in ["cubic", "pchip"]:
             interp_diff = pybamm.Interpolant(
-                x, x ** 2, y, interpolator=interpolator
+                x, x**2, y, interpolator=interpolator
             ).diff(y)
             np.testing.assert_array_almost_equal(
                 interp_diff.evaluate(y=np.array([0.397, 0.806]))[:, 0],
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 a644e0c90..0286c5efd 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
@@ -47,7 +47,7 @@ class TestCasadiConverter(unittest.TestCase):
 
         # function
         def square_plus_one(x):
-            return x ** 2 + 1
+            return x**2 + 1
 
         f = pybamm.Function(square_plus_one, b)
         self.assertEqual(f.to_casadi(), 2)
@@ -160,14 +160,14 @@ class TestCasadiConverter(unittest.TestCase):
         casadi_y = casadi.MX.sym("y", 2)
         # linear
         y_test = np.array([0.4, 0.6])
-        for interpolator in ["linear", "cubic spline"]:
+        for interpolator in ["linear", "cubic"]:
             interp = pybamm.Interpolant(x, 2 * x, y, interpolator=interpolator)
             interp_casadi = interp.to_casadi(y=casadi_y)
             f = casadi.Function("f", [casadi_y], [interp_casadi])
             np.testing.assert_array_almost_equal(interp.evaluate(y=y_test), f(y_test))
         # square
         y = pybamm.StateVector(slice(0, 1))
-        interp = pybamm.Interpolant(x, x ** 2, y, interpolator="cubic spline")
+        interp = pybamm.Interpolant(x, x**2, y, interpolator="cubic")
         interp_casadi = interp.to_casadi(y=casadi_y)
         f = casadi.Function("f", [casadi_y], [interp_casadi])
         np.testing.assert_array_almost_equal(interp.evaluate(y=y_test), f(y_test))
@@ -177,7 +177,7 @@ class TestCasadiConverter(unittest.TestCase):
         casadi_y = casadi.MX.sym("y", 1)
         data = np.tile(2 * x, (10, 1)).T
         y_test = np.array([0.4])
-        for interpolator in ["linear", "cubic spline"]:
+        for interpolator in ["linear", "cubic"]:
             interp = pybamm.Interpolant(x, data, y, interpolator=interpolator)
             interp_casadi = interp.to_casadi(y=casadi_y)
             f = casadi.Function("f", [casadi_y], [interp_casadi])
@@ -224,7 +224,7 @@ class TestCasadiConverter(unittest.TestCase):
             np.testing.assert_array_almost_equal(interp.evaluate(y=y_test), f(y_test))
         # square
         y = (pybamm.StateVector(slice(0, 1)), pybamm.StateVector(slice(0, 1)))
-        Y = (x ** 2).sum(axis=1).reshape(*[len(el) for el in x_])
+        Y = (x**2).sum(axis=1).reshape(*[len(el) for el in x_])
         interp = pybamm.Interpolant(x_, Y, y, interpolator="linear")
         interp_casadi = interp.to_casadi(y=casadi_y)
         f = casadi.Function("f", [casadi_y], [interp_casadi])
@@ -285,7 +285,7 @@ class TestCasadiConverter(unittest.TestCase):
         b = pybamm.Scalar(1)
 
         def myfunction(x, y):
-            return x + y ** 3
+            return x + y**3
 
         f = pybamm.Function(myfunction, a, b).diff(a)
         self.assert_casadi_equal(f.to_casadi(), casadi.MX(1), evalf=True)
diff --git a/tests/unit/test_parameters/test_parameter_values.py b/tests/unit/test_parameters/test_parameter_values.py
index 8aa8b5642..3ab5f613a 100644
--- a/tests/unit/test_parameters/test_parameter_values.py
+++ b/tests/unit/test_parameters/test_parameter_values.py
@@ -420,7 +420,7 @@ class TestParameterValues(unittest.TestCase):
 
     def test_process_inline_function_parameters(self):
         def D(c):
-            return c ** 2
+            return c**2
 
         parameter_values = pybamm.ParameterValues({"Diffusivity": D})
 
@@ -534,7 +534,7 @@ class TestParameterValues(unittest.TestCase):
 
         processed_func = parameter_values.process_symbol(func)
         self.assertIsInstance(processed_func, pybamm.Interpolant)
-        self.assertEqual(processed_func.evaluate(), 14.82)
+        self.assertAlmostEqual(processed_func.evaluate()[0][0], 14.82)
 
         # process differentiated function parameter
         # diff_func = func.diff(a)
@@ -600,7 +600,7 @@ class TestParameterValues(unittest.TestCase):
         processed_func = parameter_values.process_symbol(func)
         processed_interp = parameter_values.process_symbol(interp)
         np.testing.assert_array_almost_equal(
-            processed_func.evaluate(), processed_interp.evaluate(), decimal=4
+            processed_func.evaluate(), processed_interp.evaluate(), decimal=3
         )
 
         # process differentiated function parameter
@@ -807,7 +807,7 @@ class TestParameterValues(unittest.TestCase):
         var_av = pybamm.size_average(var)
 
         def dist(R):
-            return R ** 2
+            return R**2
 
         param = pybamm.ParameterValues(
             {
@@ -836,7 +836,7 @@ class TestParameterValues(unittest.TestCase):
         var2 = pybamm.Variable("var2")
         par1 = pybamm.Parameter("par1")
         par2 = pybamm.Parameter("par2")
-        expression = (3 * (par1 ** var2)) / ((var1 - par2) + var2)
+        expression = (3 * (par1**var2)) / ((var1 - par2) + var2)
 
         param = pybamm.ParameterValues({"par1": 1, "par2": 2})
         exp_param = param.process_symbol(expression)