Merge pull request #80 from pybop-team/72-54-restructure-and-add-pint…

…s-optimisation-methods

Restructure and initial Pints implementation (#54) & (#72)

.gitattributes

@@ -0,0 +1,2 @@
+# Auto detect text files and perform LF normalization
+* text=auto

.github/workflows/scheduled_tests.yaml

@@ -30,7 +30,7 @@ jobs:
       - name: Unit tests with nox
         run: |
           python -m nox -s unit
-  
+
     #M-series Mac Mini
   build-apple-mseries:
     runs-on: [self-hosted, macOS, ARM64]
@@ -64,4 +64,4 @@ jobs:
         run: |
           eval "$(pyenv init -)"
           pyenv activate pybop-${{ matrix.python-version }}
-          pyenv uninstall -f $( python --version )
+          pyenv uninstall -f $( python --version )

.gitignore

@@ -301,3 +301,6 @@ $RECYCLE.BIN/
 *.lnk
 
 # End of https://www.toptal.com/developers/gitignore/api/python,macos,windows,linux,c
+
+# Visual Studio Code settings
+.vscode/*

.pre-commit-config.yaml

@@ -7,7 +7,7 @@ repos:
     rev: "v0.0.292"
     hooks:
       - id: ruff
-        args: [--fix, --ignore=E741, --exclude=__init__.py]
+        args: ["--fix", "--ignore=E501,E402,E741", "--exclude=__init__.py"]
 
   - repo: https://github.com/nbQA-dev/nbQA
     rev: 1.7.0

examples/scripts/grad_descent.py

@@ -0,0 +1,57 @@
+import pybop
+import pints
+import numpy as np
+import matplotlib.pyplot as plt
+
+model = pybop.lithium_ion.SPMe()
+
+inputs = {
+    "Negative electrode active material volume fraction": 0.58,
+    "Positive electrode active material volume fraction": 0.44,
+    "Current function [A]": 1,
+}
+t_eval = np.arange(0, 900, 2)
+model.build(fit_parameters=inputs)
+
+values = model.predict(inputs=inputs, t_eval=t_eval)
+voltage = values["Terminal voltage [V]"].data
+time = values["Time [s]"].data
+
+sigma = 0.001
+CorruptValues = voltage + np.random.normal(0, sigma, len(voltage))
+
+# Show the generated data
+plt.figure()
+plt.xlabel("Time")
+plt.ylabel("Values")
+plt.plot(time, CorruptValues)
+plt.plot(time, voltage)
+plt.show()
+
+
+problem = pints.SingleOutputProblem(model, time, CorruptValues)
+
+# Select a score function
+score = pints.SumOfSquaresError(problem)
+
+x0 = np.array([0.48, 0.55, 1.4])
+opt = pints.OptimisationController(score, x0, method=pints.GradientDescent)
+
+opt.optimiser().set_learning_rate(0.025)
+opt.set_max_unchanged_iterations(50)
+opt.set_max_iterations(200)
+
+x1, f1 = opt.run()
+print("Estimated parameters:")
+print(x1)
+
+# Show the generated data
+simulated_values = problem.evaluate(x1[:3])
+
+plt.figure()
+plt.xlabel("Time")
+plt.ylabel("Values")
+plt.plot(time, CorruptValues)
+plt.fill_between(time, simulated_values - sigma, simulated_values + sigma, alpha=0.2)
+plt.plot(time, simulated_values)
+plt.show()

examples/scripts/mle.py

@@ -4,11 +4,11 @@
 import matplotlib.pyplot as plt
 
 model = pybop.lithium_ion.SPMe()
-model.parameter_set["Current function [A]"] = 2
 
 inputs = {
-    "Negative electrode active material volume fraction": 0.5,
-    "Positive electrode active material volume fraction": 0.6,
+    "Negative electrode active material volume fraction": 0.58,
+    "Positive electrode active material volume fraction": 0.44,
+    "Current function [A]": 1,
 }
 t_eval = np.arange(0, 900, 2)
 model.build(fit_parameters=inputs)
@@ -17,7 +17,7 @@
 voltage = values["Terminal voltage [V]"].data
 time = values["Time [s]"].data
 
-sigma = 0.01
+sigma = 0.001
 CorruptValues = voltage + np.random.normal(0, sigma, len(voltage))
 
 # Show the generated data
@@ -31,18 +31,21 @@
 
 problem = pints.SingleOutputProblem(model, time, CorruptValues)
 log_likelihood = pints.GaussianLogLikelihood(problem)
-boundaries = pints.RectangularBoundaries([0.4, 0.4, 1e-5], [0.6, 0.6, 1e-1])
+boundaries = pints.RectangularBoundaries([0.4, 0.4, 0.7, 1e-5], [0.6, 0.6, 2.1, 1e-1])
 
-x0 = np.array([0.52, 0.47, 1e-3])
-op = pints.OptimisationController(
+x0 = np.array([0.48, 0.55, 1.4, 1e-3])
+opt = pints.OptimisationController(
     log_likelihood, x0, boundaries=boundaries, method=pints.CMAES
 )
-x1, f1 = op.run()
+opt.set_max_unchanged_iterations(50)
+opt.set_max_iterations(200)
+
+x1, f1 = opt.run()
 print("Estimated parameters:")
 print(x1)
 
 # Show the generated data
-simulated_values = problem.evaluate(x1[:2])
+simulated_values = problem.evaluate(x1[:3])
 
 plt.figure()
 plt.xlabel("Time")

examples/scripts/rmse_estimation.py

@@ -1,21 +1,20 @@
 import pybop
 import pandas as pd
-import numpy as np
 
-# Form observations
+# Form dataset
 Measurements = pd.read_csv("examples/scripts/Chen_example.csv", comment="#").to_numpy()
-observations = [
-    pybop.Observed("Time [s]", Measurements[:, 0]),
-    pybop.Observed("Current function [A]", Measurements[:, 1]),
-    pybop.Observed("Voltage [V]", Measurements[:, 2]),
+dataset = [
+    pybop.Dataset("Time [s]", Measurements[:, 0]),
+    pybop.Dataset("Current function [A]", Measurements[:, 1]),
+    pybop.Dataset("Voltage [V]", Measurements[:, 2]),
 ]
 
 # Define model
 # parameter_set = pybop.ParameterSet("pybamm", "Chen2020")
 model = pybop.models.lithium_ion.SPM()
 
 # Fitting parameters
-params = [
+parameters = [
     pybop.Parameter(
         "Negative electrode active material volume fraction",
         prior=pybop.Gaussian(0.75, 0.05),
@@ -28,17 +27,25 @@
     ),
 ]
 
+# Define the cost to optimise
+cost = pybop.RMSE()
+signal = "Voltage [V]"
+
+# Build the optimisation problem
 parameterisation = pybop.Optimisation(
-    model, observations=observations, fit_parameters=params
+    cost=cost,
+    model=model,
+    optimiser=pybop.NLoptOptimize(n_param=len(parameters)),
+    parameters=parameters,
+    dataset=dataset,
+    signal=signal,
 )
 
-# get RMSE estimate using NLOpt
-results, last_optim, num_evals = parameterisation.rmse(
-    signal="Voltage [V]", method="nlopt"
-)
+# Run the optimisation problem
+x, output, final_cost, num_evals = parameterisation.run()
 
 # get MAP estimate, starting at a random initial point in parameter space
-# parameterisation.map(x0=[p.sample() for p in params])
+# parameterisation.map(x0=[p.sample() for p in parameters])
 
 # or sample from posterior
 # parameterisation.sample(1000, n_chains=4, ....)

pybop/__init__.py

@@ -7,7 +7,7 @@
 #
 
 import sys
-import os
+from os import path
 
 #
 # Version info
@@ -20,34 +20,54 @@
 # Float format: a float can be converted to a 17 digit decimal and back without
 # loss of information
 FLOAT_FORMAT = "{: .17e}"
-# Absolute path to the pybop module
-script_path = os.path.dirname(__file__)
+# Absolute path to the pybop repo
+script_path = path.dirname(__file__)
 
 #
-# Model Classes
+# Cost function class
 #
+from .costs.error_costs import RMSE, MLE, PEM, MAP
+
+#
+# Dataset class
+#
+from .datasets.base_dataset import Dataset
+
+#
+# Model classes
+#
+from .models.base_model import BaseModel
 from .models import lithium_ion
-from .models.BaseModel import BaseModel
 
 #
-# Parameterisation class
+# Main optimisation class
+#
+from .optimisation import Optimisation
+
+#
+# Optimiser class
 #
-from .parameters.parameter_set import ParameterSet
-from .parameters.parameter import Parameter
-from .datasets.observations import Observed
+from .optimisers.base_optimiser import BaseOptimiser
+from .optimisers.nlopt_optimize import NLoptOptimize
+from .optimisers.scipy_minimize import SciPyMinimize
+from .optimisers.pints_optimiser import PintsOptimiser, PintsError, PintsBoundaries
 
 #
-# Priors class
+# Parameter classes
 #
+from .parameters.base_parameter import Parameter
+from .parameters.base_parameter_set import ParameterSet
 from .parameters.priors import Gaussian, Uniform, Exponential
 
 #
-# Optimisation class
+# Problem class
 #
-from .optimisation import Optimisation
-from .optimisers import BaseOptimiser
-from .optimisers.NLoptOptimize import NLoptOptimize
-from .optimisers.SciPyMinimize import SciPyMinimize
+from ._problem import SingleOutputProblem
+
+#
+# Plotting class
+#
+from .plotting.quick_plot import QuickPlot
 
 #
 # Remove any imported modules, so we don't expose them as part of pybop

pybop/_problem.py

@@ -0,0 +1,60 @@
+import numpy as np
+
+
+class SingleOutputProblem:
+    """
+    Defines a PyBOP single output problem, follows the PINTS interface.
+    """
+
+    def __init__(self, model, parameters, signal, dataset):
+        self._model = model
+        self.parameters = {o.name: o for o in parameters}
+        self.signal = signal
+        self._dataset = dataset
+
+        if self._model._built_model is None:
+            self._model.build(fit_parameters=self.parameters)
+
+        for i, item in enumerate(self._dataset):
+            if item.name == "Time [s]":
+                self._time_data_available = True
+                self._time_data = self._dataset[i].data
+
+            if item.name == signal:
+                self._ground_truth = self._dataset[i].data
+
+        if self._time_data_available is False:
+            raise ValueError("Dataset must contain time data")
+
+        if np.any(self._time_data < 0):
+            raise ValueError("Times can not be negative.")
+        if np.any(self._time_data[:-1] >= self._time_data[1:]):
+            raise ValueError("Times must be increasing.")
+
+        if len(self._ground_truth) != len(self._time_data):
+            raise ValueError("Time data and signal data must be the same length.")
+
+    def evaluate(self, parameters):
+        """
+        Evaluate the model with the given parameters and return the signal.
+        """
+
+        y = np.asarray(
+            self._model.simulate(inputs=parameters, t_eval=self.model.time_data)[
+                self.signal
+            ].data
+        )
+
+        return y
+
+    def evaluateS1(self, parameters):
+        """
+        Evaluate the model with the given parameters and return the signal and
+        its derivatives.
+        """
+
+        y, dy_dp = self._model.simulateS1(
+            inputs=parameters, t_eval=self.model.time_data, calculate_sensitivities=True
+        )[self.signal]
+
+        return (np.asarray(y), np.asarray(dy_dp))