diff --git a/CHANGELOG.md b/CHANGELOG.md
index b7ae3616..be08914b 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -2,6 +2,7 @@
## Features
+- [#393](https://github.com/pybop-team/PyBOP/pull/383) - Adds Minkowski and SumofPower cost classes, with an example and corresponding tests.
- [#403](https://github.com/pybop-team/PyBOP/pull/403/) - Adds lychee link checking action.
## Bug Fixes
diff --git a/examples/notebooks/comparing_cost_functions.ipynb b/examples/notebooks/comparing_cost_functions.ipynb
new file mode 100644
index 00000000..54faa801
--- /dev/null
+++ b/examples/notebooks/comparing_cost_functions.ipynb
@@ -0,0 +1,609 @@
+{
+ "cells": [
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Investigating different cost functions\n",
+ "\n",
+ "In this notebook, we take a look at the different fitting cost functions offered in PyBOP. Cost functions for fitting problems conventionally describe the distance between two points (the target and the prediction) which is to be minimised via PyBOP's optimisation algorithms. \n",
+ "\n",
+ "First, we install and import the required packages below."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 1,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "Requirement already satisfied: pip in /Users/engs2510/.pyenv/versions/3.12.2/envs/pybop-3.12/lib/python3.12/site-packages (24.1.2)\n",
+ "Requirement already satisfied: ipywidgets in /Users/engs2510/.pyenv/versions/3.12.2/envs/pybop-3.12/lib/python3.12/site-packages (8.1.3)\n",
+ "Requirement already satisfied: comm>=0.1.3 in /Users/engs2510/.pyenv/versions/3.12.2/envs/pybop-3.12/lib/python3.12/site-packages (from ipywidgets) (0.2.2)\n",
+ "Requirement already satisfied: ipython>=6.1.0 in /Users/engs2510/.pyenv/versions/3.12.2/envs/pybop-3.12/lib/python3.12/site-packages (from ipywidgets) (8.23.0)\n",
+ "Requirement already satisfied: traitlets>=4.3.1 in /Users/engs2510/.pyenv/versions/3.12.2/envs/pybop-3.12/lib/python3.12/site-packages (from ipywidgets) (5.14.2)\n",
+ "Requirement already satisfied: widgetsnbextension~=4.0.11 in /Users/engs2510/.pyenv/versions/3.12.2/envs/pybop-3.12/lib/python3.12/site-packages (from ipywidgets) (4.0.11)\n",
+ "Requirement already satisfied: jupyterlab-widgets~=3.0.11 in /Users/engs2510/.pyenv/versions/3.12.2/envs/pybop-3.12/lib/python3.12/site-packages (from ipywidgets) (3.0.11)\n",
+ "Requirement already satisfied: decorator in /Users/engs2510/.pyenv/versions/3.12.2/envs/pybop-3.12/lib/python3.12/site-packages (from ipython>=6.1.0->ipywidgets) (5.1.1)\n",
+ "Requirement already satisfied: jedi>=0.16 in /Users/engs2510/.pyenv/versions/3.12.2/envs/pybop-3.12/lib/python3.12/site-packages (from ipython>=6.1.0->ipywidgets) (0.19.1)\n",
+ "Requirement already satisfied: matplotlib-inline in /Users/engs2510/.pyenv/versions/3.12.2/envs/pybop-3.12/lib/python3.12/site-packages (from ipython>=6.1.0->ipywidgets) (0.1.6)\n",
+ "Requirement already satisfied: prompt-toolkit<3.1.0,>=3.0.41 in /Users/engs2510/.pyenv/versions/3.12.2/envs/pybop-3.12/lib/python3.12/site-packages (from ipython>=6.1.0->ipywidgets) (3.0.43)\n",
+ "Requirement already satisfied: pygments>=2.4.0 in /Users/engs2510/.pyenv/versions/3.12.2/envs/pybop-3.12/lib/python3.12/site-packages (from ipython>=6.1.0->ipywidgets) (2.17.2)\n",
+ "Requirement already satisfied: stack-data in /Users/engs2510/.pyenv/versions/3.12.2/envs/pybop-3.12/lib/python3.12/site-packages (from ipython>=6.1.0->ipywidgets) (0.6.3)\n",
+ "Requirement already satisfied: pexpect>4.3 in /Users/engs2510/.pyenv/versions/3.12.2/envs/pybop-3.12/lib/python3.12/site-packages (from ipython>=6.1.0->ipywidgets) (4.9.0)\n",
+ "Requirement already satisfied: parso<0.9.0,>=0.8.3 in /Users/engs2510/.pyenv/versions/3.12.2/envs/pybop-3.12/lib/python3.12/site-packages (from jedi>=0.16->ipython>=6.1.0->ipywidgets) (0.8.4)\n",
+ "Requirement already satisfied: ptyprocess>=0.5 in /Users/engs2510/.pyenv/versions/3.12.2/envs/pybop-3.12/lib/python3.12/site-packages (from pexpect>4.3->ipython>=6.1.0->ipywidgets) (0.7.0)\n",
+ "Requirement already satisfied: wcwidth in /Users/engs2510/.pyenv/versions/3.12.2/envs/pybop-3.12/lib/python3.12/site-packages (from prompt-toolkit<3.1.0,>=3.0.41->ipython>=6.1.0->ipywidgets) (0.2.13)\n",
+ "Requirement already satisfied: executing>=1.2.0 in /Users/engs2510/.pyenv/versions/3.12.2/envs/pybop-3.12/lib/python3.12/site-packages (from stack-data->ipython>=6.1.0->ipywidgets) (2.0.1)\n",
+ "Requirement already satisfied: asttokens>=2.1.0 in /Users/engs2510/.pyenv/versions/3.12.2/envs/pybop-3.12/lib/python3.12/site-packages (from stack-data->ipython>=6.1.0->ipywidgets) (2.4.1)\n",
+ "Requirement already satisfied: pure-eval in /Users/engs2510/.pyenv/versions/3.12.2/envs/pybop-3.12/lib/python3.12/site-packages (from stack-data->ipython>=6.1.0->ipywidgets) (0.2.2)\n",
+ "Requirement already satisfied: six>=1.12.0 in /Users/engs2510/.pyenv/versions/3.12.2/envs/pybop-3.12/lib/python3.12/site-packages (from asttokens>=2.1.0->stack-data->ipython>=6.1.0->ipywidgets) (1.16.0)\n",
+ "Note: you may need to restart the kernel to use updated packages.\n",
+ "Note: you may need to restart the kernel to use updated packages.\n"
+ ]
+ }
+ ],
+ "source": [
+ "%pip install --upgrade pip ipywidgets\n",
+ "%pip install pybop -q\n",
+ "\n",
+ "import numpy as np\n",
+ "import plotly.graph_objects as go\n",
+ "import plotly.io as pio\n",
+ "\n",
+ "pio.renderers.default = \"svg\"\n",
+ "import pybop"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "For this notebook, we need to construct parameters, a model and a problem class before we can compare differing cost functions. We start with two parameters, but this is an arbitrary selection and can be expanded given the model and data in question."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 2,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "parameters = pybop.Parameters(\n",
+ " pybop.Parameter(\n",
+ " \"Positive electrode thickness [m]\",\n",
+ " prior=pybop.Gaussian(7.56e-05, 0.5e-05),\n",
+ " bounds=[65e-06, 10e-05],\n",
+ " ),\n",
+ " pybop.Parameter(\n",
+ " \"Positive particle radius [m]\",\n",
+ " prior=pybop.Gaussian(5.22e-06, 0.5e-06),\n",
+ " bounds=[2e-06, 9e-06],\n",
+ " ),\n",
+ ")"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Next, we will construct the Single Particle Model (SPM) with the Chen2020 parameter set, but like the above, this is an arbitrary selection and can be replaced with any PyBOP model."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 3,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "parameter_set = pybop.ParameterSet.pybamm(\"Chen2020\")\n",
+ "model = pybop.lithium_ion.SPM(parameter_set=parameter_set)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Next, as we will need reference data to compare our model predictions to (via the cost function), we will create synthetic data from the model constructed above. "
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 4,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "t_eval = np.arange(0, 900, 2)\n",
+ "values = model.predict(t_eval=t_eval)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "We can then construct the PyBOP dataset class with the synthetic data as,"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 5,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "dataset = pybop.Dataset(\n",
+ " {\n",
+ " \"Time [s]\": t_eval,\n",
+ " \"Current function [A]\": values[\"Current [A]\"].data,\n",
+ " \"Voltage [V]\": values[\"Voltage [V]\"].data,\n",
+ " }\n",
+ ")"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Now, we can put this all together and construct the problem class. In this situation, we are going to compare differing fitting cost functions, so we construct the `FittingProblem`."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 6,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "problem = pybop.FittingProblem(model, parameters, dataset)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "### Sum of Squared Errors and Root Mean Squared Error\n",
+ "\n",
+ "First, let's start with two commonly-used cost functions: the sum of squared errors (SSE) and the root mean squared error (RMSE). Constructing these classes is very concise in PyBOP, and only requires the problem class."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 7,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "cost_SSE = pybop.SumSquaredError(problem)\n",
+ "cost_RMSE = pybop.RootMeanSquaredError(problem)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Now, we can investigate how these functions differ when fitting the parameters. To acquire the cost value for each of these, we can simply use the call method of the constructed class, such as:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 8,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "1.2690291451182834e-09"
+ ]
+ },
+ "execution_count": 8,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "cost_SSE([7.56e-05, 5.22e-06])"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Alternatively, we can use the `Parameters` class for this,"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 9,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "[7.56e-05 5.22e-06]\n"
+ ]
+ },
+ {
+ "data": {
+ "text/plain": [
+ "1.2690291451182834e-09"
+ ]
+ },
+ "execution_count": 9,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "print(parameters.current_value())\n",
+ "cost_SSE(parameters.current_value())"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "If we want to generate a random sample of candidate solutions from the parameter class prior, we can also do that as:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 10,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "[7.60957550e-05 5.48691392e-06]\n"
+ ]
+ },
+ {
+ "data": {
+ "text/plain": [
+ "0.014466013735507651"
+ ]
+ },
+ "execution_count": 10,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "sample = parameters.rvs()\n",
+ "print(sample)\n",
+ "cost_SSE(sample)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "#### Comparing RMSE and SSE\n",
+ "\n",
+ "Now, let's vary one of the parameters, and keep a fixed value for the other, to create a scatter plot comparing the cost values for the RMSE and SSE functions."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 11,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "image/svg+xml": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "x_range = np.linspace(4.72e-06, 5.72e-06, 75)\n",
+ "y_SSE = []\n",
+ "y_RMSE = []\n",
+ "for i in x_range:\n",
+ " y_SSE.append(cost_SSE([7.56e-05, i]))\n",
+ " y_RMSE.append(cost_RMSE([7.56e-05, i]))\n",
+ "\n",
+ "fig = go.Figure()\n",
+ "fig.add_trace(go.Scatter(x=x_range, y=y_SSE, mode=\"lines\", name=\"SSE\"))\n",
+ "fig.add_trace(go.Scatter(x=x_range, y=y_RMSE, mode=\"lines\", name=\"RMSE\"))\n",
+ "fig.show()"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "In this situation, it's clear that the curvature of the SSE cost is greater than that of the RMSE. This can improve the rate of convergence for certain optimisation algorithms. However, with incorrect hyperparameter values, larger gradients can also result in the algorithm not converging due to sampling locations outside of the \"cost valley\", e.g. infeasible parameter values."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "### Minkowski distance\n",
+ "\n",
+ "Next, let's investigate the Minkowski distance. The Minkowski cost takes a general form, which allows for hyperparameter calibration on the cost function itself, given by\n",
+ "\n",
+ "$\\mathcal{L_p} = \\displaystyle \\Big(\\sum_i |\\hat{y_i}-y_i|^p\\Big)^{1/p}$\n",
+ "\n",
+ "where $p ≥ 0$ is the order of the Minkowski distance.\n",
+ "\n",
+ "For $p = 1$, it is the Manhattan distance. \n",
+ "For $p = 2$, it is the Euclidean distance. \n",
+ "For $p ≥ 1$, the Minkowski distance is a metric, but for $04.8μ5μ5.2μ5.4μ5.6μ00.050.10.150.20.25RMSE*Nsqrt(SSE)Minkowski"
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "y_minkowski = []\n",
+ "for i in x_range:\n",
+ " y_minkowski.append(cost_minkowski([7.56e-05, i]))\n",
+ "\n",
+ "fig = go.Figure()\n",
+ "fig.add_trace(\n",
+ " go.Scatter(\n",
+ " x=x_range,\n",
+ " y=np.asarray(y_RMSE) * np.sqrt(len(t_eval)),\n",
+ " mode=\"lines\",\n",
+ " name=\"RMSE*N\",\n",
+ " )\n",
+ ")\n",
+ "fig.add_trace(\n",
+ " go.Scatter(\n",
+ " x=x_range,\n",
+ " y=np.sqrt(y_SSE),\n",
+ " mode=\"lines\",\n",
+ " line=dict(dash=\"dash\"),\n",
+ " name=\"sqrt(SSE)\",\n",
+ " )\n",
+ ")\n",
+ "fig.add_trace(\n",
+ " go.Scatter(\n",
+ " x=x_range, y=y_minkowski, mode=\"lines\", line=dict(dash=\"dot\"), name=\"Minkowski\"\n",
+ " )\n",
+ ")\n",
+ "fig.show()"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "As expected, these lines lie on top of one another. Now, let's take a look at how the Minkowski cost changes for different orders, `p`."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 14,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "p_orders = np.append(0.75, np.linspace(1, 3, 5))\n",
+ "y_minkowski = tuple(\n",
+ " [pybop.Minkowski(problem, p=j)([7.56e-05, i]) for i in x_range] for j in p_orders\n",
+ ")"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 15,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "image/svg+xml": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "fig = go.Figure()\n",
+ "for k, _ in enumerate(p_orders):\n",
+ " fig.add_trace(\n",
+ " go.Scatter(x=x_range, y=y_minkowski[k], mode=\"lines\", name=f\"Minkowski {_}\")\n",
+ " )\n",
+ "fig.update_yaxes(range=[0, np.max(y_minkowski[2])])\n",
+ "fig.show()"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "As seen above, the Minkowski cost allows for a range of different cost functions to be created. This provides users with another hyperparameter to calibrate for optimisation algorithm convergence. This addition does expand the global search space, and should be carefully considered before deciding upon.\n",
+ "\n",
+ "### Sum of Power\n",
+ "Next, we introduce a similar cost function, the `SumofPower` implementation. This cost function is the $p$-th power of the Minkowski distance of order $p$. It provides a generalised formulation for the Sum of Squared Errors (SSE) cost function, and is given by,\n",
+ "\n",
+ "$\\mathcal{L_p} = \\displaystyle \\sum_i |\\hat{y_i}-y_i|^p$\n",
+ "\n",
+ "where $p ≥ 0$ is the power order. A few special cases include,\n",
+ "\n",
+ "$p = 1$: Sum of Absolute Differences\n",
+ "$p = 2$: Sum of Squared Differences\n",
+ "\n",
+ "Next we repeat the above examples with the addition of the `SumofPower` class."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 16,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "cost_sumofpower = pybop.SumofPower(problem, p=2)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 17,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "image/svg+xml": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "y_sumofpower = []\n",
+ "for i in x_range:\n",
+ " y_sumofpower.append(cost_sumofpower([7.56e-05, i]))\n",
+ "\n",
+ "fig = go.Figure()\n",
+ "fig.add_trace(\n",
+ " go.Scatter(\n",
+ " x=x_range,\n",
+ " y=np.asarray(y_RMSE) * np.sqrt(len(t_eval)),\n",
+ " mode=\"lines\",\n",
+ " name=\"RMSE*N\",\n",
+ " )\n",
+ ")\n",
+ "fig.add_trace(\n",
+ " go.Scatter(\n",
+ " x=x_range,\n",
+ " y=y_SSE,\n",
+ " mode=\"lines\",\n",
+ " line=dict(dash=\"dash\"),\n",
+ " name=\"SSE\",\n",
+ " )\n",
+ ")\n",
+ "fig.add_trace(\n",
+ " go.Scatter(\n",
+ " x=x_range,\n",
+ " y=y_sumofpower,\n",
+ " mode=\"lines\",\n",
+ " line=dict(dash=\"dot\"),\n",
+ " name=\"Sum of Power\",\n",
+ " )\n",
+ ")\n",
+ "fig.show()"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "As expected, the `SumofPower` with order `p=2` equates to the `SSE` implementation. Next, we compare the `Minkowski` to the `SumofPower`,"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 18,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "p_orders = np.append(0.75, np.linspace(1, 2, 2))\n",
+ "\n",
+ "y_minkowski = tuple(\n",
+ " [pybop.Minkowski(problem, p=j)([7.56e-05, i]) for i in x_range] for j in p_orders\n",
+ ")\n",
+ "\n",
+ "y_sumofpower = tuple(\n",
+ " [pybop.SumofPower(problem, p=j)([7.56e-05, i]) for i in x_range] for j in p_orders\n",
+ ")"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 19,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "image/svg+xml": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "fig = go.Figure()\n",
+ "for k, _ in enumerate(p_orders):\n",
+ " fig.add_trace(\n",
+ " go.Scatter(x=x_range, y=y_minkowski[k], mode=\"lines\", name=f\"Minkowski {_}\")\n",
+ " )\n",
+ " fig.add_trace(\n",
+ " go.Scatter(\n",
+ " x=x_range,\n",
+ " y=y_sumofpower[k],\n",
+ " mode=\"lines\",\n",
+ " line=dict(dash=\"dash\"),\n",
+ " name=f\"Sum of Power {_}\",\n",
+ " )\n",
+ " )\n",
+ "fig.update_yaxes(range=[0, 2.5])\n",
+ "fig.show()"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "The figure demonstrates the distinct behaviour of the `Minkowski` distance and the `SumofPower` function. One notable difference is the effect of the `1/p` exponent in the `Minkowski` distance, which has a linearising impact on the response. This linearisation can enhance the robustness of certain optimisation algorithms, potentially making them less sensitive to outliers or extreme values. However, this increased robustness may come at the cost of a slower convergence rate, as the linearised response might require more iterations to reach the optimal solution. In contrast, the `SumofPower` function does not exhibit this linearising effect, which can lead to faster convergence in some cases but may be more susceptible to the influence of outliers or extreme values.\n",
+ "\n",
+ "In this notebook, we've shown the different fitting cost functions offered in PyBOP. Selection between these functions can affect the optimisation result in the case that the optimiser hyperparameter values are not properly calibrated. "
+ ]
+ }
+ ],
+ "metadata": {
+ "kernelspec": {
+ "display_name": "pybop-3.12",
+ "language": "python",
+ "name": "python3"
+ },
+ "language_info": {
+ "codemirror_mode": {
+ "name": "ipython",
+ "version": 3
+ },
+ "file_extension": ".py",
+ "mimetype": "text/x-python",
+ "name": "python",
+ "nbconvert_exporter": "python",
+ "pygments_lexer": "ipython3",
+ "version": "3.12.2"
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}
diff --git a/examples/scripts/spm_AdamW.py b/examples/scripts/spm_AdamW.py
index 796849be..46c96350 100644
--- a/examples/scripts/spm_AdamW.py
+++ b/examples/scripts/spm_AdamW.py
@@ -53,7 +53,7 @@ def noise(sigma):
problem = pybop.FittingProblem(
model, parameters, dataset, signal=signal, init_soc=init_soc
)
-cost = pybop.RootMeanSquaredError(problem)
+cost = pybop.Minkowski(problem, p=2)
optim = pybop.AdamW(
cost,
verbose=True,
diff --git a/examples/scripts/spm_MAP.py b/examples/scripts/spm_MAP.py
index f613eccb..40f7a279 100644
--- a/examples/scripts/spm_MAP.py
+++ b/examples/scripts/spm_MAP.py
@@ -2,12 +2,16 @@
import pybop
+# Set variables
+sigma = 0.002
+init_soc = 0.7
+
# Construct and update initial parameter values
parameter_set = pybop.ParameterSet.pybamm("Chen2020")
parameter_set.update(
{
- "Negative electrode active material volume fraction": 0.63,
- "Positive electrode active material volume fraction": 0.51,
+ "Negative electrode active material volume fraction": 0.43,
+ "Positive electrode active material volume fraction": 0.56,
}
)
@@ -18,38 +22,51 @@
parameters = pybop.Parameters(
pybop.Parameter(
"Negative electrode active material volume fraction",
- prior=pybop.Gaussian(0.6, 0.05),
- bounds=[0.5, 0.8],
+ prior=pybop.Uniform(0.3, 0.8),
+ bounds=[0.3, 0.8],
+ initial_value=0.653,
true_value=parameter_set["Negative electrode active material volume fraction"],
),
pybop.Parameter(
"Positive electrode active material volume fraction",
- prior=pybop.Gaussian(0.48, 0.05),
+ prior=pybop.Uniform(0.3, 0.8),
bounds=[0.4, 0.7],
+ initial_value=0.657,
true_value=parameter_set["Positive electrode active material volume fraction"],
),
)
-# Generate data
-sigma = 0.005
-t_eval = np.arange(0, 900, 3)
-values = model.predict(t_eval=t_eval)
-corrupt_values = values["Voltage [V]"].data + np.random.normal(0, sigma, len(t_eval))
+# Generate data and corrupt it with noise
+experiment = pybop.Experiment(
+ [
+ (
+ "Discharge at 0.5C for 3 minutes (4 second period)",
+ "Charge at 0.5C for 3 minutes (4 second period)",
+ ),
+ ]
+)
+values = model.predict(
+ init_soc=init_soc, experiment=experiment, inputs=parameters.true_value()
+)
+corrupt_values = values["Voltage [V]"].data + np.random.normal(
+ 0, sigma, len(values["Voltage [V]"].data)
+)
# Form dataset
dataset = pybop.Dataset(
{
- "Time [s]": t_eval,
+ "Time [s]": values["Time [s]"].data,
"Current function [A]": values["Current [A]"].data,
"Voltage [V]": corrupt_values,
}
)
# Generate problem, cost function, and optimisation class
-problem = pybop.FittingProblem(model, parameters, dataset)
+problem = pybop.FittingProblem(model, parameters, dataset, init_soc=init_soc)
cost = pybop.MAP(problem, pybop.GaussianLogLikelihoodKnownSigma, sigma0=sigma)
optim = pybop.AdamW(
cost,
+ sigma0=0.05,
max_unchanged_iterations=20,
min_iterations=20,
max_iterations=100,
@@ -61,7 +78,7 @@
print("Estimated parameters:", x)
# Plot the timeseries output
-pybop.quick_plot(problem, problem_inputs=x[0:2], title="Optimised Comparison")
+pybop.quick_plot(problem, problem_inputs=x, title="Optimised Comparison")
# Plot convergence
pybop.plot_convergence(optim)
@@ -73,5 +90,5 @@
pybop.plot2d(cost, steps=15)
# Plot the cost landscape with optimisation path
-bounds = np.asarray([[0.55, 0.77], [0.48, 0.68]])
+bounds = np.asarray([[0.35, 0.7], [0.45, 0.625]])
pybop.plot2d(optim, bounds=bounds, steps=15)
diff --git a/examples/scripts/spm_SNES.py b/examples/scripts/spm_SNES.py
index 93046d63..7cddfad9 100644
--- a/examples/scripts/spm_SNES.py
+++ b/examples/scripts/spm_SNES.py
@@ -35,7 +35,7 @@
# Generate problem, cost function, and optimisation class
problem = pybop.FittingProblem(model, parameters, dataset)
-cost = pybop.SumSquaredError(problem)
+cost = pybop.SumofPower(problem, p=2)
optim = pybop.SNES(cost, max_iterations=100)
x, final_cost = optim.run()
diff --git a/pybop/__init__.py b/pybop/__init__.py
index 92c869f4..c9aabb8d 100644
--- a/pybop/__init__.py
+++ b/pybop/__init__.py
@@ -85,6 +85,8 @@
from .costs.fitting_costs import (
RootMeanSquaredError,
SumSquaredError,
+ Minkowski,
+ SumofPower,
ObserverCost,
)
from .costs.design_costs import (
diff --git a/pybop/costs/_likelihoods.py b/pybop/costs/_likelihoods.py
index c6e5916d..e65f02fc 100644
--- a/pybop/costs/_likelihoods.py
+++ b/pybop/costs/_likelihoods.py
@@ -38,19 +38,17 @@ def __init__(self, problem: BaseProblem, sigma0: Union[list[float], float]):
self.sigma2 = sigma0**2.0
self._offset = -0.5 * self.n_time_data * np.log(2 * np.pi * self.sigma2)
self._multip = -1 / (2.0 * self.sigma2)
- self._dl = np.ones(self.n_parameters)
def _evaluate(self, inputs: Inputs, grad: Union[None, np.ndarray] = None) -> float:
"""
Evaluates the Gaussian log-likelihood for the given parameters with known sigma.
"""
y = self.problem.evaluate(inputs)
- if any(
- len(y.get(key, [])) != len(self._target.get(key, [])) for key in self.signal
- ):
- return -np.inf # prediction length doesn't match target
- e = np.sum(
+ if not self.verify_prediction(y):
+ return -np.inf
+
+ e = np.asarray(
[
np.sum(
self._offset
@@ -60,7 +58,7 @@ def _evaluate(self, inputs: Inputs, grad: Union[None, np.ndarray] = None) -> flo
]
)
- return e if self.n_outputs != 1 else e.item()
+ return e.item() if self.n_outputs == 1 else np.sum(e)
def _evaluateS1(self, inputs: Inputs) -> tuple[float, np.ndarray]:
"""
@@ -68,10 +66,8 @@ def _evaluateS1(self, inputs: Inputs) -> tuple[float, np.ndarray]:
"""
y, dy = self.problem.evaluateS1(inputs)
- if any(
- len(y.get(key, [])) != len(self._target.get(key, [])) for key in self.signal
- ):
- return -np.inf, -self._dl
+ if not self.verify_prediction(y):
+ return -np.inf, -self._de * np.ones(self.n_parameters)
likelihood = self._evaluate(inputs)
@@ -125,7 +121,6 @@ def __init__(
self.sigma = Parameters()
self._add_sigma_parameters(sigma0)
self.parameters.join(self.sigma)
- self._dl = np.ones(self.n_parameters)
def _add_sigma_parameters(self, sigma0):
sigma0 = [sigma0] if not isinstance(sigma0, list) else sigma0
@@ -195,12 +190,10 @@ def _evaluate(self, inputs: Inputs, grad: Union[None, np.ndarray] = None) -> flo
return -np.inf
y = self.problem.evaluate(self.problem.parameters.as_dict())
- if any(
- len(y.get(key, [])) != len(self._target.get(key, [])) for key in self.signal
- ):
- return -np.inf # prediction length doesn't match target
+ if not self.verify_prediction(y):
+ return -np.inf
- e = np.sum(
+ e = np.asarray(
[
np.sum(
self._logpi
@@ -212,7 +205,7 @@ def _evaluate(self, inputs: Inputs, grad: Union[None, np.ndarray] = None) -> flo
]
)
- return e if self.n_outputs != 1 else e.item()
+ return e.item() if self.n_outputs == 1 else np.sum(e)
def _evaluateS1(self, inputs: Inputs) -> tuple[float, np.ndarray]:
"""
@@ -232,13 +225,11 @@ def _evaluateS1(self, inputs: Inputs) -> tuple[float, np.ndarray]:
sigma = self.sigma.current_value()
if np.any(sigma <= 0):
- return -np.inf, -self._dl
+ return -np.inf, -self._de * np.ones(self.n_parameters)
y, dy = self.problem.evaluateS1(self.problem.parameters.as_dict())
- if any(
- len(y.get(key, [])) != len(self._target.get(key, [])) for key in self.signal
- ):
- return -np.inf, -self._dl
+ if not self.verify_prediction(y):
+ return -np.inf, -self._de * np.ones(self.n_parameters)
likelihood = self._evaluate(inputs)
@@ -302,11 +293,14 @@ def _evaluate(self, inputs: Inputs, grad=None) -> float:
float
The maximum a posteriori cost.
"""
- log_likelihood = self.likelihood._evaluate(inputs)
log_prior = sum(
self.parameters[key].prior.logpdf(value) for key, value in inputs.items()
)
+ if not np.isfinite(log_prior).any():
+ return -np.inf
+
+ log_likelihood = self.likelihood._evaluate(inputs)
posterior = log_likelihood + log_prior
return posterior
@@ -331,10 +325,13 @@ def _evaluateS1(self, inputs: Inputs) -> tuple[float, np.ndarray]:
ValueError
If an error occurs during the calculation of the cost or gradient.
"""
- log_likelihood, dl = self.likelihood._evaluateS1(inputs)
log_prior = sum(
self.parameters[key].prior.logpdf(value) for key, value in inputs.items()
)
+ if not np.isfinite(log_prior).any():
+ return -np.inf, -self._de * np.ones(self.n_parameters)
+
+ log_likelihood, dl = self.likelihood._evaluateS1(inputs)
# Compute a finite difference approximation of the gradient of the log prior
delta = self.parameters.initial_value() * self.gradient_step
diff --git a/pybop/costs/base_cost.py b/pybop/costs/base_cost.py
index c7f7a7bb..d40bbb99 100644
--- a/pybop/costs/base_cost.py
+++ b/pybop/costs/base_cost.py
@@ -1,3 +1,5 @@
+from typing import Union
+
from pybop import BaseProblem
from pybop.parameters.parameter import Inputs, Parameters
@@ -25,6 +27,7 @@ class BaseCost:
def __init__(self, problem=None):
self.parameters = Parameters()
self.problem = problem
+ self.set_fail_gradient()
if isinstance(self.problem, BaseProblem):
self._target = self.problem._target
self.parameters.join(self.problem.parameters)
@@ -35,20 +38,20 @@ def __init__(self, problem=None):
def n_parameters(self):
return len(self.parameters)
- def __call__(self, x, grad=None):
+ def __call__(self, inputs: Union[Inputs, list], grad=None):
"""
Call the evaluate function for a given set of parameters.
"""
- return self.evaluate(x, grad)
+ return self.evaluate(inputs, grad)
- def evaluate(self, x, grad=None):
+ def evaluate(self, inputs: Union[Inputs, list], grad=None):
"""
Call the evaluate function for a given set of parameters.
Parameters
----------
- x : array-like
- The parameters for which to evaluate the cost.
+ inputs : Inputs or array-like
+ The parameters for which to compute the cost and gradient.
grad : array-like, optional
An array to store the gradient of the cost function with respect
to the parameters.
@@ -63,7 +66,7 @@ def evaluate(self, x, grad=None):
ValueError
If an error occurs during the calculation of the cost.
"""
- inputs = self.parameters.verify(x)
+ inputs = self.parameters.verify(inputs)
try:
return self._evaluate(inputs, grad)
@@ -100,27 +103,27 @@ def _evaluate(self, inputs: Inputs, grad=None):
"""
raise NotImplementedError
- def evaluateS1(self, x):
+ def evaluateS1(self, inputs: Union[Inputs, list]):
"""
Call _evaluateS1 for a given set of parameters.
Parameters
----------
- x : array-like
+ inputs : Inputs or array-like
The parameters for which to compute the cost and gradient.
Returns
-------
tuple
A tuple containing the cost and the gradient. The cost is a float,
- and the gradient is an array-like of the same length as `x`.
+ and the gradient is an array-like of the same length as `inputs`.
Raises
------
ValueError
If an error occurs during the calculation of the cost or gradient.
"""
- inputs = self.parameters.verify(x)
+ inputs = self.parameters.verify(inputs)
try:
return self._evaluateS1(inputs)
@@ -144,7 +147,7 @@ def _evaluateS1(self, inputs: Inputs):
-------
tuple
A tuple containing the cost and the gradient. The cost is a float,
- and the gradient is an array-like of the same length as `x`.
+ and the gradient is an array-like of the same length as `inputs`.
Raises
------
@@ -152,3 +155,40 @@ def _evaluateS1(self, inputs: Inputs):
If the method has not been implemented by the subclass.
"""
raise NotImplementedError
+
+ def set_fail_gradient(self, de: float = 1.0):
+ """
+ Set the fail gradient to a specified value.
+
+ The fail gradient is used if an error occurs during the calculation
+ of the gradient. This method allows updating the default gradient value.
+
+ Parameters
+ ----------
+ de : float
+ The new fail gradient value to be used.
+ """
+ if not isinstance(de, float):
+ de = float(de)
+ self._de = de
+
+ def verify_prediction(self, y):
+ """
+ Verify that the prediction matches the target data.
+
+ Parameters
+ ----------
+ y : dict
+ The model predictions.
+
+ Returns
+ -------
+ bool
+ True if the prediction matches the target data, otherwise False.
+ """
+ if any(
+ len(y.get(key, [])) != len(self._target.get(key, [])) for key in self.signal
+ ):
+ return False
+
+ return True
diff --git a/pybop/costs/fitting_costs.py b/pybop/costs/fitting_costs.py
index a0cba766..ac17f3ea 100644
--- a/pybop/costs/fitting_costs.py
+++ b/pybop/costs/fitting_costs.py
@@ -20,9 +20,6 @@ class RootMeanSquaredError(BaseCost):
def __init__(self, problem):
super().__init__(problem)
- # Default fail gradient
- self._de = 1.0
-
def _evaluate(self, inputs: Inputs, grad=None):
"""
Calculate the root mean square error for a given set of parameters.
@@ -43,9 +40,8 @@ def _evaluate(self, inputs: Inputs, grad=None):
"""
prediction = self.problem.evaluate(inputs)
- for key in self.signal:
- if len(prediction.get(key, [])) != len(self._target.get(key, [])):
- return np.float64(np.inf) # prediction doesn't match target
+ if not self.verify_prediction(prediction):
+ return np.inf
e = np.asarray(
[
@@ -54,10 +50,7 @@ def _evaluate(self, inputs: Inputs, grad=None):
]
)
- if self.n_outputs == 1:
- return e.item()
- else:
- return np.sum(e)
+ return e.item() if self.n_outputs == 1 else np.sum(e)
def _evaluateS1(self, inputs: Inputs):
"""
@@ -72,7 +65,7 @@ def _evaluateS1(self, inputs: Inputs):
-------
tuple
A tuple containing the cost and the gradient. The cost is a float,
- and the gradient is an array-like of the same length as `x`.
+ and the gradient is an array-like of the same length as `inputs`.
Raises
------
@@ -80,12 +73,8 @@ def _evaluateS1(self, inputs: Inputs):
If an error occurs during the calculation of the cost or gradient.
"""
y, dy = self.problem.evaluateS1(inputs)
-
- for key in self.signal:
- if len(y.get(key, [])) != len(self._target.get(key, [])):
- e = np.float64(np.inf)
- de = self._de * np.ones(self.n_parameters)
- return e, de
+ if not self.verify_prediction(y):
+ return np.inf, self._de * np.ones(self.n_parameters)
r = np.asarray([y[signal] - self._target[signal] for signal in self.signal])
e = np.sqrt(np.mean(r**2, axis=1))
@@ -96,21 +85,6 @@ def _evaluateS1(self, inputs: Inputs):
else:
return np.sum(e), np.sum(de, axis=1)
- def set_fail_gradient(self, de):
- """
- Set the fail gradient to a specified value.
-
- The fail gradient is used if an error occurs during the calculation
- of the gradient. This method allows updating the default gradient value.
-
- Parameters
- ----------
- de : float
- The new fail gradient value to be used.
- """
- de = float(de)
- self._de = de
-
class SumSquaredError(BaseCost):
"""
@@ -133,9 +107,6 @@ class SumSquaredError(BaseCost):
def __init__(self, problem):
super().__init__(problem)
- # Default fail gradient
- self._de = 1.0
-
def _evaluate(self, inputs: Inputs, grad=None):
"""
Calculate the sum of squared errors for a given set of parameters.
@@ -151,13 +122,12 @@ def _evaluate(self, inputs: Inputs, grad=None):
Returns
-------
float
- The sum of squared errors.
+ The Sum of Squared Error.
"""
prediction = self.problem.evaluate(inputs)
- for key in self.signal:
- if len(prediction.get(key, [])) != len(self._target.get(key, [])):
- return np.float64(np.inf) # prediction doesn't match target
+ if not self.verify_prediction(prediction):
+ return np.inf
e = np.asarray(
[
@@ -165,10 +135,8 @@ def _evaluate(self, inputs: Inputs, grad=None):
for signal in self.signal
]
)
- if self.n_outputs == 1:
- return e.item()
- else:
- return np.sum(e)
+
+ return e.item() if self.n_outputs == 1 else np.sum(e)
def _evaluateS1(self, inputs: Inputs):
"""
@@ -183,7 +151,7 @@ def _evaluateS1(self, inputs: Inputs):
-------
tuple
A tuple containing the cost and the gradient. The cost is a float,
- and the gradient is an array-like of the same length as `x`.
+ and the gradient is an array-like of the same length as `inputs`.
Raises
------
@@ -191,11 +159,8 @@ def _evaluateS1(self, inputs: Inputs):
If an error occurs during the calculation of the cost or gradient.
"""
y, dy = self.problem.evaluateS1(inputs)
- for key in self.signal:
- if len(y.get(key, [])) != len(self._target.get(key, [])):
- e = np.float64(np.inf)
- de = self._de * np.ones(self.n_parameters)
- return e, de
+ if not self.verify_prediction(y):
+ return np.inf, self._de * np.ones(self.n_parameters)
r = np.asarray([y[signal] - self._target[signal] for signal in self.signal])
e = np.sum(np.sum(r**2, axis=0), axis=0)
@@ -203,20 +168,211 @@ def _evaluateS1(self, inputs: Inputs):
return e, de
- def set_fail_gradient(self, de):
+
+class Minkowski(BaseCost):
+ """
+ The Minkowski distance is a generalisation of several distance metrics,
+ including the Euclidean and Manhattan distances. It is defined as:
+
+ .. math::
+ L_p(x, y) = ( \\sum_i |x_i - y_i|^p )^(1/p)
+
+ where p > 0 is the order of the Minkowski distance. For p ≥ 1, the
+ Minkowski distance is a metric. For 0 < p < 1, it is not a metric, as it
+ does not satisfy the triangle inequality, although a metric can be
+ obtained by removing the (1/p) exponent.
+
+ Special cases:
+
+ * p = 1: Manhattan distance
+ * p = 2: Euclidean distance
+ * p → ∞: Chebyshev distance (not implemented as yet)
+
+ This class implements the Minkowski distance as a cost function for
+ optimisation problems, allowing for flexible distance-based optimisation
+ across various problem domains.
+
+ Attributes
+ ----------
+ p : float, optional
+ The order of the Minkowski distance.
+ """
+
+ def __init__(self, problem, p: float = 2.0):
+ super().__init__(problem)
+ if p < 0:
+ raise ValueError(
+ "The order of the Minkowski distance must be greater than 0."
+ )
+ elif not np.isfinite(p):
+ raise ValueError(
+ "For p = infinity, an implementation of the Chebyshev distance is required."
+ )
+ self.p = float(p)
+
+ def _evaluate(self, inputs: Inputs, grad=None):
"""
- Set the fail gradient to a specified value.
+ Calculate the Minkowski cost for a given set of parameters.
- The fail gradient is used if an error occurs during the calculation
- of the gradient. This method allows updating the default gradient value.
+ Parameters
+ ----------
+ inputs : Inputs
+ The parameters for which to compute the cost and gradient.
+
+ Returns
+ -------
+ float
+ The Minkowski cost.
+ """
+ prediction = self.problem.evaluate(inputs)
+ if not self.verify_prediction(prediction):
+ return np.inf
+
+ e = np.asarray(
+ [
+ np.sum(np.abs(prediction[signal] - self._target[signal]) ** self.p)
+ ** (1 / self.p)
+ for signal in self.signal
+ ]
+ )
+
+ return e.item() if self.n_outputs == 1 else np.sum(e)
+
+ def _evaluateS1(self, inputs):
+ """
+ Compute the cost and its gradient with respect to the parameters.
Parameters
----------
- de : float
- The new fail gradient value to be used.
+ inputs : Inputs
+ The parameters for which to compute the cost and gradient.
+
+ Returns
+ -------
+ tuple
+ A tuple containing the cost and the gradient. The cost is a float,
+ and the gradient is an array-like of the same length as `inputs`.
+
+ Raises
+ ------
+ ValueError
+ If an error occurs during the calculation of the cost or gradient.
"""
- de = float(de)
- self._de = de
+ y, dy = self.problem.evaluateS1(inputs)
+ if not self.verify_prediction(y):
+ return np.inf, self._de * np.ones(self.n_parameters)
+
+ r = np.asarray([y[signal] - self._target[signal] for signal in self.signal])
+ e = np.asarray(
+ [
+ np.sum(np.abs(y[signal] - self._target[signal]) ** self.p)
+ ** (1 / self.p)
+ for signal in self.signal
+ ]
+ )
+ de = np.sum(
+ np.sum(r ** (self.p - 1) * dy.T, axis=2)
+ / (e ** (self.p - 1) + np.finfo(float).eps),
+ axis=1,
+ )
+
+ return np.sum(e), de
+
+
+class SumofPower(BaseCost):
+ """
+ The Sum of Power [1] is a generalised cost function based on the p-th power
+ of absolute differences between two vectors. It is defined as:
+
+ .. math::
+ C_p(x, y) = \\sum_i |x_i - y_i|^p
+
+ where p ≥ 0 is the power order.
+
+ This class implements the Sum of Power as a cost function for
+ optimisation problems, allowing for flexible power-based optimisation
+ across various problem domains.
+
+ Special cases:
+
+ * p = 1: Sum of Absolute Differences
+ * p = 2: Sum of Squared Differences
+ * p → ∞: Maximum Absolute Difference
+
+ Note that this is not normalised, unlike distance metrics. To get a
+ distance metric, you would need to take the p-th root of the result.
+
+ [1]: https://mathworld.wolfram.com/PowerSum.html
+
+ Attributes:
+ p : float, optional
+ The power order for Sum of Power.
+ """
+
+ def __init__(self, problem, p: float = 2.0):
+ super().__init__(problem)
+ if p < 0:
+ raise ValueError("The order of 'p' must be greater than 0.")
+ elif not np.isfinite(p):
+ raise ValueError("p = np.inf is not yet supported.")
+ self.p = float(p)
+
+ def _evaluate(self, inputs: Inputs, grad=None):
+ """
+ Calculate the Sum of Power cost for a given set of parameters.
+
+ Parameters
+ ----------
+ inputs : Inputs
+ The parameters for which to compute the cost and gradient.
+
+ Returns
+ -------
+ float
+ The Sum of Power cost.
+ """
+ prediction = self.problem.evaluate(inputs)
+ if not self.verify_prediction(prediction):
+ return np.inf
+
+ e = np.asarray(
+ [
+ np.sum(np.abs(prediction[signal] - self._target[signal]) ** self.p)
+ for signal in self.signal
+ ]
+ )
+
+ return e.item() if self.n_outputs == 1 else np.sum(e)
+
+ def _evaluateS1(self, inputs):
+ """
+ Compute the cost and its gradient with respect to the parameters.
+
+ Parameters
+ ----------
+ inputs : Inputs
+ The parameters for which to compute the cost and gradient.
+
+ Returns
+ -------
+ tuple
+ A tuple containing the cost and the gradient. The cost is a float,
+ and the gradient is an array-like of the same length as `inputs`.
+
+ Raises
+ ------
+ ValueError
+ If an error occurs during the calculation of the cost or gradient.
+ """
+ y, dy = self.problem.evaluateS1(inputs)
+ if not self.verify_prediction(y):
+ return np.inf, self._de * np.ones(self.n_parameters)
+
+ r = np.asarray([y[signal] - self._target[signal] for signal in self.signal])
+ e = np.sum(np.sum(np.abs(r) ** self.p))
+ de = self.p * np.sum(np.sum(r ** (self.p - 1) * dy.T, axis=2), axis=1)
+
+ return e, de
class ObserverCost(BaseCost):
@@ -269,7 +425,7 @@ def evaluateS1(self, inputs: Inputs):
-------
tuple
A tuple containing the cost and the gradient. The cost is a float,
- and the gradient is an array-like of the same length as `x`.
+ and the gradient is an array-like of the same length as `inputs`.
Raises
------
diff --git a/pybop/optimisers/scipy_optimisers.py b/pybop/optimisers/scipy_optimisers.py
index 30499cdb..7528fbe0 100644
--- a/pybop/optimisers/scipy_optimisers.py
+++ b/pybop/optimisers/scipy_optimisers.py
@@ -141,6 +141,22 @@ def _set_up_optimiser(self):
# Nest this option within an options dictionary for SciPy minimize
self._options["options"]["maxiter"] = self.unset_options.pop(key)
+ def cost_wrapper(self, x):
+ """
+ Scale the cost function, preserving the sign convention, and eliminate nan values
+ """
+ self.log["x"].append([x])
+
+ if not self._options["jac"]:
+ cost = self.cost(x) / self._cost0
+ if np.isinf(cost):
+ self.inf_count += 1
+ cost = 1 + 0.9**self.inf_count # for fake finite gradient
+ return cost if self.minimising else -cost
+
+ L, dl = self.cost.evaluateS1(x)
+ return (L, dl) if self.minimising else (-L, -dl)
+
def _run_optimiser(self):
"""
Executes the optimisation process using SciPy's minimize function.
@@ -150,6 +166,7 @@ def _run_optimiser(self):
result : scipy.optimize.OptimizeResult
The result of the optimisation including the optimised parameter values and cost.
"""
+ self.inf_count = 0
# Add callback storing history of parameter values
def callback(intermediate_result: OptimizeResult):
@@ -162,7 +179,7 @@ def callback(intermediate_result: OptimizeResult):
self._cost0 = np.abs(self.cost(self.x0))
if np.isinf(self._cost0):
for _i in range(1, self.num_resamples):
- self.x0 = self.parameters.rvs(1)[0]
+ self.x0 = self.parameters.rvs()
self._cost0 = np.abs(self.cost(self.x0))
if not np.isinf(self._cost0):
break
@@ -171,27 +188,8 @@ def callback(intermediate_result: OptimizeResult):
"The initial parameter values return an infinite cost."
)
- # Scale the cost function, preserving the sign convention, and eliminate nan values
- self.inf_count = 0
-
- if not self._options["jac"]:
-
- def cost_wrapper(x):
- self.log["x"].append([x])
- cost = self.cost(x) / self._cost0
- if np.isinf(cost):
- self.inf_count += 1
- cost = 1 + 0.9**self.inf_count # for fake finite gradient
- return cost if self.minimising else -cost
- elif self._options["jac"] is True:
-
- def cost_wrapper(x):
- self.log["x"].append([x])
- L, dl = self.cost.evaluateS1(x)
- return L, dl if self.minimising else -L, -dl
-
return minimize(
- cost_wrapper,
+ self.cost_wrapper,
self.x0,
bounds=self._scipy_bounds,
callback=callback,
diff --git a/pybop/parameters/parameter.py b/pybop/parameters/parameter.py
index 093abaed..2614027c 100644
--- a/pybop/parameters/parameter.py
+++ b/pybop/parameters/parameter.py
@@ -51,7 +51,7 @@ def __init__(
self.set_bounds(bounds)
self.margin = 1e-4
- def rvs(self, n_samples, random_state=None):
+ def rvs(self, n_samples: int = 1, random_state=None):
"""
Draw random samples from the parameter's prior distribution.
@@ -61,7 +61,7 @@ def rvs(self, n_samples, random_state=None):
Parameters
----------
n_samples : int
- The number of samples to draw.
+ The number of samples to draw (default: 1).
Returns
-------
@@ -322,7 +322,7 @@ def update(self, initial_values=None, values=None, bounds=None):
else:
param.set_bounds(bounds=bounds[i])
- def rvs(self, n_samples: int) -> list:
+ def rvs(self, n_samples: int = 1) -> np.ndarray:
"""
Draw random samples from each parameter's prior distribution.
@@ -332,7 +332,7 @@ def rvs(self, n_samples: int) -> list:
Parameters
----------
n_samples : int
- The number of samples to draw.
+ The number of samples to draw (default: 1).
Returns
-------
@@ -353,7 +353,7 @@ def rvs(self, n_samples: int) -> list:
all_samples.append(samples)
- return all_samples
+ return np.concatenate(all_samples)
def get_sigma0(self) -> list:
"""
@@ -417,7 +417,7 @@ def get_bounds_for_plotly(self):
bounds : numpy.ndarray
An array of shape (n_parameters, 2) containing the bounds for each parameter.
"""
- bounds = np.empty((len(self), 2))
+ bounds = np.zeros((len(self), 2))
for i, param in enumerate(self.param.values()):
if param.applied_prior_bounds:
@@ -427,7 +427,7 @@ def get_bounds_for_plotly(self):
UserWarning,
stacklevel=2,
)
- elif param.bounds is not None:
+ if param.bounds is not None:
bounds[i] = param.bounds
else:
raise ValueError("All parameters require bounds for plotting.")
diff --git a/pybop/plotting/plot2d.py b/pybop/plotting/plot2d.py
index 781b697b..d5c85574 100644
--- a/pybop/plotting/plot2d.py
+++ b/pybop/plotting/plot2d.py
@@ -1,5 +1,6 @@
import sys
import warnings
+from typing import Union
import numpy as np
from scipy.interpolate import griddata
@@ -10,7 +11,7 @@
def plot2d(
cost_or_optim,
gradient: bool = False,
- bounds: np.ndarray = None,
+ bounds: Union[np.ndarray, None] = None,
steps: int = 10,
show: bool = True,
use_optim_log: bool = False,
diff --git a/tests/integration/test_spm_parameterisations.py b/tests/integration/test_spm_parameterisations.py
index d7acaf7d..5e9d6b00 100644
--- a/tests/integration/test_spm_parameterisations.py
+++ b/tests/integration/test_spm_parameterisations.py
@@ -46,6 +46,8 @@ def init_soc(self, request):
pybop.GaussianLogLikelihood,
pybop.RootMeanSquaredError,
pybop.SumSquaredError,
+ pybop.SumofPower,
+ pybop.Minkowski,
pybop.MAP,
]
)
@@ -78,6 +80,8 @@ def spm_costs(self, model, parameters, cost_class, init_soc):
return cost_class(
problem, pybop.GaussianLogLikelihoodKnownSigma, sigma0=self.sigma0
)
+ elif cost_class in [pybop.SumofPower, pybop.Minkowski]:
+ return cost_class(problem, p=2)
else:
return cost_class(problem)
@@ -101,6 +105,7 @@ def test_spm_optimisers(self, optimiser, spm_costs):
"cost": spm_costs,
"max_iterations": 250,
"absolute_tolerance": 1e-6,
+ "max_unchanged_iterations": 55,
}
# Add sigma0 to ground truth for GaussianLogLikelihood
@@ -111,21 +116,18 @@ def test_spm_optimisers(self, optimiser, spm_costs):
if isinstance(spm_costs, pybop.MAP):
for i in spm_costs.parameters.keys():
spm_costs.parameters[i].prior = pybop.Uniform(
- 0.4, 2.0
+ 0.2, 2.0
) # Increase range to avoid prior == np.inf
# Set sigma0 and create optimiser
sigma0 = 0.05 if isinstance(spm_costs, pybop.MAP) else None
optim = optimiser(sigma0=sigma0, **common_args)
- # Set max unchanged iterations for BasePintsOptimisers
- if issubclass(optimiser, pybop.BasePintsOptimiser):
- optim.set_max_unchanged_iterations(iterations=55)
-
# AdamW will use lowest sigma0 for learning rate, so allow more iterations
if issubclass(optimiser, (pybop.AdamW, pybop.IRPropMin)) and isinstance(
spm_costs, pybop.GaussianLogLikelihood
):
- optim = optimiser(max_unchanged_iterations=75, **common_args)
+ common_args["max_unchanged_iterations"] = 75
+ optim = optimiser(**common_args)
initial_cost = optim.cost(x0)
x, final_cost = optim.run()
diff --git a/tests/unit/test_cost.py b/tests/unit/test_cost.py
index 80e44bea..09372376 100644
--- a/tests/unit/test_cost.py
+++ b/tests/unit/test_cost.py
@@ -74,6 +74,8 @@ def problem(self, model, parameters, dataset, signal, request):
params=[
pybop.RootMeanSquaredError,
pybop.SumSquaredError,
+ pybop.Minkowski,
+ pybop.SumofPower,
pybop.ObserverCost,
pybop.MAP,
]
@@ -84,6 +86,8 @@ def cost(self, problem, request):
return cls(problem)
elif cls in [pybop.MAP]:
return cls(problem, pybop.GaussianLogLikelihoodKnownSigma)
+ elif cls in [pybop.Minkowski, pybop.SumofPower]:
+ return cls(problem, p=2)
elif cls in [pybop.ObserverCost]:
inputs = problem.parameters.initial_value()
state = problem._model.reinit(inputs)
@@ -153,6 +157,20 @@ def test_MAP(self, problem):
with pytest.raises(ValueError, match="must be a subclass of BaseLikelihood"):
pybop.MAP(problem, self.InvalidLikelihood, sigma0=0.1)
+ # Non finite prior
+ parameter = pybop.Parameter(
+ "Negative electrode active material volume fraction",
+ prior=pybop.Uniform(0.55, 0.6),
+ )
+ problem_non_finite = pybop.FittingProblem(
+ problem.model, parameter, problem.dataset, signal=problem.signal
+ )
+ likelihood = pybop.MAP(
+ problem_non_finite, pybop.GaussianLogLikelihoodKnownSigma, sigma0=0.01
+ )
+ assert not np.isfinite(likelihood([0.7]))
+ assert not np.isfinite(likelihood.evaluateS1([0.7])[0])
+
@pytest.mark.unit
def test_costs(self, cost):
if isinstance(cost, pybop.BaseLikelihood):
@@ -181,7 +199,7 @@ def test_costs(self, cost):
# Test option setting
cost.set_fail_gradient(1)
- if isinstance(cost, pybop.SumSquaredError):
+ if not isinstance(cost, (pybop.ObserverCost, pybop.MAP)):
e, de = cost.evaluateS1([0.5])
assert np.isscalar(e)
@@ -210,8 +228,27 @@ def test_costs(self, cost):
with pytest.raises(TypeError, match="Inputs must be a dictionary or numeric."):
cost(["StringInputShouldNotWork"])
- # Test treatment of simulations that terminated early
- # by variation of the cut-off voltage.
+ @pytest.mark.unit
+ def test_minkowski(self, problem):
+ # Incorrect order
+ with pytest.raises(ValueError, match="The order of the Minkowski distance"):
+ pybop.Minkowski(problem, p=-1)
+ with pytest.raises(
+ ValueError,
+ match="For p = infinity, an implementation of the Chebyshev distance is required.",
+ ):
+ pybop.Minkowski(problem, p=np.inf)
+
+ @pytest.mark.unit
+ def test_SumofPower(self, problem):
+ # Incorrect order
+ with pytest.raises(
+ ValueError, match="The order of 'p' must be greater than 0."
+ ):
+ pybop.SumofPower(problem, p=-1)
+
+ with pytest.raises(ValueError, match="p = np.inf is not yet supported."):
+ pybop.SumofPower(problem, p=np.inf)
@pytest.mark.parametrize(
"cost_class",
diff --git a/tests/unit/test_optimisation.py b/tests/unit/test_optimisation.py
index c61d2fae..bb6d08fa 100644
--- a/tests/unit/test_optimisation.py
+++ b/tests/unit/test_optimisation.py
@@ -352,6 +352,7 @@ class RandomClass:
[
(0.85, 0.2, False),
(0.85, 0.001, True),
+ (1.0, 0.5, False),
],
)
def test_scipy_prior_resampling(
@@ -385,6 +386,10 @@ def test_scipy_prior_resampling(
else:
opt.run()
+ # Test cost_wrapper inf return
+ cost = opt.cost_wrapper(np.array([0.9]))
+ assert cost in [1.729, 1.81, 1.9]
+
@pytest.mark.unit
def test_halting(self, cost):
# Test max evalutions