Updates and typo to article, adds posterior plot, updates design plot…

…, w/ corresponding code

joss/paper.md

@@ -5,6 +5,7 @@ tags:
   - batteries
   - battery models
   - parameterisation
+  - parameter inference
   - design optimisation
 authors:
   - name: Brady Planden
@@ -32,13 +33,14 @@ affiliations:
    index: 3
  - name:  Mathematics Institute, University of Warwick, Coventry, UK
    index: 4
-date: 30 June 2024
+date: 28 November 2024
 bibliography: paper.bib
+repository: https://github.com/pybop-team/PyBOP
 ---
 
 # Summary
 
-The Python Battery Optimisation and Parameterisation (`PyBOP`) package provides a set of methods for the parameterisation and optimisation of battery models, offering both Bayesian and frequentist approaches with example workflows to assist the user. `PyBOP` has been developed to enable parameter identification of various battery models, including the electrochemical and equivalent circuit models provided by the popular open-source package `PyBaMM` [@Sulzer:2021]. PyBOP
+The Python Battery Optimisation and Parameterisation (`PyBOP`) package provides a set of methods for the parameterisation and optimisation of battery models, offering both Bayesian and frequentist approaches with example workflows to assist the user. `PyBOP` has been developed to enable parameter identification of various battery models, including the electrochemical and equivalent circuit models provided by the popular open-source package `PyBaMM` [@Sulzer:2021].
 
 Similarly, `PyBOP` can be used for parameter design optimisation under user-defined operating conditions across a variety of model structures. `PyBOP` enables battery model parameterisation across a range of methods with diagnostics on the performance and convergence of the identified or optimised parameters. The identified parameters can be used for prediction, on-line control and design optimisation, all of which support improved battery utilisation and development.
 
@@ -84,23 +86,23 @@ Likewise, the current algorithms available for optimisation tasks are presented
 : The currently supported optimisation algorithms classified by candidate solution type, inclusive of gradient information. (*) Scipy minimize has gradient and non-gradient methods. \label{tab:optimisers}
 
 | Gradient-based                                       | Evolutionary Strategies               | (Meta)heuristic      |
-| :--------------------------------------------------- | :------------------------------------ | :------------------- |
+|:-----------------------------------------------------|:--------------------------------------|:---------------------|
 | Adaptive moment estimation with weight decay (AdamW) | Covariance matrix adaptation (CMA-ES) | Particle swarm (PSO) |
 | Improved resilient backpropagation (iRProp-)         | Exponential natural (xNES)            | Nelder-Mead          |
-| Gradient descent                                     | Separable natural (sNES)              |                      |
+| Gradient descent                                     | Separable natural (sNES)              | Cuckoo search        |
 | SciPy minimize (*)                                   | SciPy differential evolution          |                      |
-|                                                      | Cuckoo search                         |                      |
+|                                                      |                                       |                      |
 
  As discussed above, `PyBOP` offers Bayesian inference methods such as Maximum a Posteriori (MAP) presented alongside the point-based methods in \autoref{tab:subclasses}; however, for a full Bayesian framework, Monte Carlo sampling is implemented within `PyBOP`. These methods construct a posterior distribution on the inference parameters which can used for uncertainty and practical identifiability. The individual sampler classes are currently composed within `PyBOP` from the `PINTS` library, with a base sampling class implemented for interoperability and direct integration to the `PyBOP` model, problem, and likelihood classes. The currently supported samplers are presented in \autoref{tab:samplers}.
 
 : PyBOP's supported sampling methods separated based on candidate suggestion method. \label{tab:samplers}
 
-| Hamiltonian-based | Adaptive                   | Slice Sampling       | Evolutionary           | Other                        |
-|:------------------|:---------------------------|:---------------------|:-----------------------|:-----------------------------|
-| Monomial Gamma    | Delayed Rejection Adaptive | Slice Doubling       | Differential Evolution | Metropolis Random Walk       |
-| No-U-Turn         | Haario Bardenet            | Slice Rank Shrinking |                        | Emcee Hammer                 |
-| Hamiltonian       | Rao Blackwell              | Slice Stepout        |                        | Metropolis Adjusted Langevin |
-| Relativistic      | Haario                     |                      |                        |                              |
+| Hamiltonian-based | Adaptive                   | Slice Sampling | Evolutionary           | Other                        |
+|:------------------|:---------------------------|:---------------|:-----------------------|:-----------------------------|
+| Monomial Gamma    | Delayed Rejection Adaptive | Doubling       | Differential Evolution | Metropolis Random Walk       |
+| No-U-Turn         | Haario Bardenet            | Rank Shrinking |                        | Emcee Hammer                 |
+| Hamiltonian       | Haario                     | Stepout        |                        | Metropolis Adjusted Langevin |
+| Relativistic      | Rao Blackwell              |                |                        |                              |
 
 # Background
 
@@ -134,7 +136,7 @@ Battery model parameterisation is challenging due to the high number of paramete
 A generic data fitting optimisation problem may be formulated as:
 \begin{equation}
 \min_{\mathbf{\theta}} ~ \mathcal{L}_{(\mathbf{y}_i)}(\mathbf{\theta}) ~~~
-\textrm{subject to equations (\ref{dynamics})\textrm{-}(\autoref{initial_conditions})}
+\textrm{subject to equations (\ref{dynamics})\textrm{-}(\ref{initial_conditions})}
 \label{eqn:parameterisation}
 \end{equation}
 
@@ -146,13 +148,13 @@ We next demonstrate the fitting of synthetic data of which the system parameters
 
  As gradient information is available for this problem, the choice of distance-based cost function and optimiser is not constrained. Due to the vastly different magnitudes of the two parameters, we apply two of the parameter transformations offered by `PyBOP`, namely the log transformation for the negative particle diffusivity and the scaled transformation (with a coefficient of 100) for the contact resistance. This application transforms the optimisers search space, enabling a shared step-size between the parameters; however, in general is not required. As a demonstration of `PyBOP`'s parameterisation capabilities, \autoref{fig:convergence-min-max} (left) shows the rate of convergence for each of the distance-minimising cost functions, while \autoref{fig:convergence-min-max} (right) displays analogous results for maximising a likelihood. Here, the optimisation is performed with SciPy Minimize using the gradient-based L-BFGS-B method.
 
-![Convergence in the likelihood functions obtained using various likelihood functions and SciPy minimize. \label{fig:convergence-min-max}](figures/joss/converge.png){ width=100% }
+![Convergence in the likelihood functions obtained using various likelihood functions and L-BFGS-B algorithm. \label{fig:convergence-min-max}](figures/joss/converge.png){ width=100% }
 
 Furthermore, we can also compare the performance of the various optimisation algorithms divided by category: gradient-based in \autoref{fig:optimiser-inference} (left), evolution strategies in \autoref{fig:optimiser-inference} (middle) and (meta)heuristics in \autoref{fig:optimiser-inference} (right) for a sum squared error cost function. Note that optimiser performance depends on the cost landscape, prior information, and corresponding hyperparameters for each specific problem.
 
-![Convergence in the parameter values obtained for various (meta)heuristics.  \label{fig:optimiser-inference}](figures/joss/optimisers_parameters.png){ width=100% }
+![Convergence in the parameter values obtained for the various optimisation algorithms provided by `PyBOP`.  \label{fig:optimiser-inference}](figures/joss/optimisers_parameters.png){ width=100% }
 
-![Counter plot for all optimisers.  \label{fig:optimiser-inference}](figures/joss/contour_total.png){ width=100% }
+![Cost landscape contour plot with corresponding optimisation traces. The top row represents the gradient-based optimisers, the middle is the evolution-based, and the bottom is the (meta)heuristics. The order left to right aligns with the entries in \autoref{tab:optimisers}.  \label{fig:optimiser-inference}](figures/joss/contour_total.png){ width=100% }
 
 This parameterisation task can also be approached from the Bayesian perspective, which we will present below using PyBOP's sampler methods. The optimisation equation presented in equation \autoref{eqn:parameterisation} does not represent the Bayesian parameter identification task, and as such we introduce Bayes Theorem as,
 
@@ -164,7 +166,7 @@ P(\theta|D) = \frac{P(D|\theta)P(\theta)}{P(D)}
 where, $P(\theta|D)$ is the posterior and represents the probability density function of the parameter. $P(D|\theta)$ is the likelihood function and assesses the parameter values alongside a noise model. $P(\theta)$ encapsulates the prior knowledge on the parameters, and finally $P(D)$ is the model evidence and acts as a normalising constant such that the final posterior is a correctly scaled density function.
 Our aim in parameter inference is to identify the parameter values with the highest probability, which can be presented from point-based metric or as the posterior distribution, which provides additional information on the uncertainty of the identified parameters. To acquire this posterior distribution, we provide Monte-Carlo sampling methods. These methods sample from the posterior through a variety of methods, including gradient-based such as No-U-Turn [@NUTS:2011] and Hamiltonian [@Hamiltonian:2011] as well as heuristic methods such as Differential Evolution [@DiffEvolution:2006], and finally conventional methods based on random sampling with rejection criteria [@metropolis:1953]. PyBOP offers a sampling class which provides an interface for these samplers, which are supported from the Probabilistic Inference of Noise Time-Series (PINTS) package. \autoref{fig:posteriors} below presents the sampled posterior for the synthetic workflow described above, using an adaptive covariance based sampler, Haario Bardenet [@Haario:2001].
 
-![[Placeholder]: Posterior distributions for inferred parameters.  \label{fig:posteriors}](figures/joss/optimisers_parameters.png){ width=100% }
+![Posterior distributions for model parameters alongside identified noise on the observations. Shaded area denotes confidence bounds for each parameter.  \label{fig:posteriors}](figures/joss/posteriors.png){ width=100% }
 
 ## Design optimisation
 
@@ -179,24 +181,22 @@ Design optimisation can be written in the form of a constrained optimisation pro
 in which $\mathcal{L} : \mathbf{\theta} \mapsto [0,\infty)$ is a cost function that quantifies the desirability
 of the design and $\Omega$ is the set of allowable parameter values.
 
-As an example, let us consider the problem of maximising gravimetric energy density subject to constraints on two of the geometric electrode parameters [@Couto:2023]. For this example, we use `PyBaMM`'s implementation of the single particle model with electrolyte (SPMe) to investigate the impact of the positive electrode thickness and active material volume fraction on the target cost. In this problem, we estimate the 1C rate from the theoretical capacity for each iteration of the design.
+As an example, let us consider the problem of maximising gravimetric energy density subject to constraints on two of the geometric electrode parameters [@Couto:2023]. For this example, we use `PyBaMM`'s implementation of the single particle model with electrolyte (SPMe) to investigate the impact of the positive electrode thickness and active material volume fraction on the target cost. As the active material volume fraction is linked to the electrode porosity, the porosity is defined with a driven reference from the volume fraction. In this problem, we estimate the 1C rate from the theoretical capacity for each iteration of the design. For this example, we employ the particle swarm optimisation (PSO) algorithm.
 
-![The voltage profiles for a 1C discharge predicted using the initial and optimised parameter values. \label{fig:design_gravimetric}](figures/joss/design.png){ width=100% }
+![The gravimetric landscape alongside the corresponding initial and optimised voltage profiles for a 1C discharge. \label{fig:design_gravimetric}](figures/joss/design.png){ width=100% }
 
-\autoref{fig:design_gravimetric} (left) shows the optimiser's search for the maximum gravimetric energy density within the parameter space. For this example, we employ the particle swarm optimisation (PSO) algorithm. The predicted improvement in the discharge profile between the initial and optimised parameter values is shown in \autoref{fig:design_gravimetric} (right).
+\autoref{fig:design_gravimetric} (left) shows the optimiser's search across the gravimetric energy density parameter space. The predicted improvement in the discharge profile between the initial and optimised parameter values (right) for their respective applied 1C current.
 
 # Acknowledgements
 
 We gratefully acknowledge all [contributors](https://github.com/pybop-team/PyBOP?tab=readme-ov-file#contributors-) to this package. This work was supported by the Faraday Institution Multiscale Modelling (MSM) project (ref. FIRG059), UKRI's Horizon Europe Guarantee (ref. 10038031), and EU IntelLiGent project (ref. 101069765).
 
 [//]: # (# Discussion Points)
 
+[//]: # (EIS numerical identification)
+
 [//]: # (- Performance discussion (multiprocessing / JAX))
 
 [//]: # (- Feasibility checks on identified parameters)
 
-[//]: # (- Add posterior figure)
-
-[//]: # (- Change optimiser comparison to cost convergence alongside optimiser landscape plots (surface))
-
 # References

joss/param_plots.py

@@ -1,8 +1,10 @@
 # A script to generate parameterisation plots for the JOSS paper.
 
+import matplotlib.pyplot as plt
 import numpy as np
 import plotly
 import pybamm
+from matplotlib.ticker import ScalarFormatter
 
 import pybop
 from pybop.plot import PlotlyManager
@@ -12,14 +14,14 @@
 
 # Choose which plots to show and save
 create_plot = {}
-create_plot["simulation"] = True
-create_plot["landscape"] = True
-create_plot["minimising"] = True
-create_plot["maximising"] = True
-create_plot["gradient"] = True  # takes longest
+create_plot["simulation"] = False
+create_plot["landscape"] = False
+create_plot["minimising"] = False
+create_plot["maximising"] = False
+create_plot["gradient"] = True
 create_plot["evolution"] = True
 create_plot["heuristic"] = True
-create_plot["posteriors"] = True
+create_plot["posteriors"] = False
 
 
 # Parameter set and model definition
@@ -286,15 +288,15 @@
 # Categorise the optimisers
 gradient_optimiser_classes = [
     pybop.AdamW,
-    pybop.GradientDescent,
     pybop.IRPropMin,
+    pybop.GradientDescent,
     pybop.SciPyMinimize,  # with L-BFGS-B and jac=True
 ]
 evolution_optimiser_classes = [
-    pybop.SciPyDifferentialEvolution,  # most iterations
-    pybop.XNES,
     pybop.CMAES,
-    pybop.SNES,  # least iterations
+    pybop.XNES,
+    pybop.SNES,
+    pybop.SciPyDifferentialEvolution,
 ]
 heuristic_optimiser_classes = [
     pybop.PSO,
@@ -317,7 +319,7 @@
         elif optimiser is pybop.GradientDescent:
             kwargs["sigma0"] = [1.2, 0.3]
         elif optimiser is pybop.AdamW:
-            kwargs["sigma0"] = [0.25, 0.08]
+            kwargs["sigma0"] = 0.2
         # elif optimiser is pybop.IRPropMin:
         #     kwargs["sigma0"] = [4e-3,2e-2]
 
@@ -332,6 +334,11 @@
             **kwargs,
         )
 
+        if optimiser is pybop.AdamW:
+            optim.optimiser.b1 = 0.8
+            optim.optimiser.b2 = 0.9
+            optim.optimiser.lam = 0.01
+
         # Run optimisation
         results = optim.run()
         print("True parameter values:", parameters.true_value())
@@ -490,3 +497,111 @@
     summary.plot_trace()
     summary.plot_chains()
     summary.plot_posterior()
+
+    # Enable LaTeX for text rendering
+    plt.rcParams.update(
+        {
+            "text.usetex": True,
+            "font.family": "serif",  # Use a LaTeX-compatible font
+            "font.serif": ["Computer Modern"],  # Default LaTeX font
+        }
+    )
+
+    # Sample from the distributions
+    neg_particle_samples = summary.all_samples[:, 0]
+    contact_samples = summary.all_samples[:, 1]
+    sigma_samples = summary.all_samples[:, 2]
+
+    # Create a grid for subplots
+    fig = plt.figure(figsize=(13, 6))
+    gs = fig.add_gridspec(1, 3, height_ratios=[1])
+
+    # Function to format axis
+    def format_axis(ax):
+        formatter = ScalarFormatter(useMathText=True)
+        formatter.set_powerlimits((-1, 1))  # Set limits to use scientific notation
+        formatter.set_scientific(True)  # Ensure scientific notation
+        ax.xaxis.set_major_formatter(formatter)
+        ax.yaxis.set_major_formatter(formatter)
+        ax.tick_params(axis="both", which="major", labelsize=14)
+
+    # Top subplot for Neg Particle
+    ax1 = fig.add_subplot(gs[0, 0])
+    ax1.hist(
+        neg_particle_samples,
+        bins=50,
+        density=False,
+        alpha=0.6,
+        color="tab:red",
+        label="$R_n$",
+    )
+    ax1.set_title("Negative Particle Radius", fontsize=22)
+    ax1.set_xlabel(r"m", fontsize=18)
+    ax1.set_ylabel("Density", fontsize=18)
+    ax1.set_xlim(
+        parameter_set["Negative particle diffusivity [m2.s-1]"] * 0.95,
+        parameter_set["Negative particle diffusivity [m2.s-1]"] * 1.07,
+    )
+    # ax1.set_ylim(0, 1e4)
+    ax1.tick_params(axis="both", which="major", labelsize=18)
+    ax1.axvspan(
+        summary.get_summary_statistics()[("ci_lower")][0],
+        summary.get_summary_statistics()[("ci_upper")][0],
+        alpha=0.1,
+        color="tab:red",
+    )
+    format_axis(ax1)
+
+    # Top right subplot for Contact Resistance
+    ax2 = fig.add_subplot(gs[0, 1])
+    ax2.hist(
+        contact_samples,
+        bins=50,
+        density=False,
+        alpha=0.6,
+        color="tab:blue",
+        label="$R_c$",
+    )
+    ax2.set_title("Contact Resistance", fontsize=22)
+    ax2.set_xlabel(r"$\Omega$", fontsize=18)
+    ax2.set_xlim(
+        parameter_set["Contact resistance [Ohm]"] * 0.95,
+        parameter_set["Contact resistance [Ohm]"] * 1.05,
+    )
+    # ax2.set_ylim(0, 1e3)
+    ax2.tick_params(axis="both", which="major", labelsize=18)
+    ax2.axvspan(
+        summary.get_summary_statistics()[("ci_lower")][1],
+        summary.get_summary_statistics()[("ci_upper")][1],
+        alpha=0.1,
+        color="tab:blue",
+    )
+    format_axis(ax2)
+
+    # Bottom right subplot for sigma
+    ax3 = fig.add_subplot(gs[0, 2])
+    ax3.hist(
+        sigma_samples,
+        bins=50,
+        density=False,
+        alpha=0.6,
+        color="tab:purple",
+        label=r"$\sigma$",
+    )
+    ax3.set_title(r"Observation Noise, $\sigma$", fontsize=22)
+    ax3.set_xlabel("V", fontsize=18)
+    ax3.set_xlim(4e-3, 5.5e-3)
+    # ax3.set_ylim(0, 10 * 1.1)
+    ax3.tick_params(axis="both", which="major", labelsize=18)
+    ax3.axvspan(
+        summary.get_summary_statistics()[("ci_lower")][2],
+        summary.get_summary_statistics()[("ci_upper")][2],
+        alpha=0.1,
+        color="tab:purple",
+    )
+    format_axis(ax3)
+
+    # Adjust layout
+    plt.tight_layout()
+    plt.savefig("joss/figures/joss/posteriors.png")
+    plt.show()