From 02cbb14255fcfe2095d4e1fb41c12ec60bb4ba39 Mon Sep 17 00:00:00 2001 From: NicolaCourtier <45851982+NicolaCourtier@users.noreply.github.com> Date: Wed, 5 Jun 2024 16:25:27 +0100 Subject: [PATCH] Add parameters class (#322) * Add Parameters class * Update error type * Update parameters in tests * Update parameters in examples * Update parameters in notebooks * Update parameter name * Add parameters len and tests * Delay unpacking of parameter properties * Replace sample_initial_conditions * Update test_thevenin_parameterisation.py * Test without flaky * Move x0 def to parameters * Add parameters.update * Define n_parameters for model * Make fit_keys a property * Add parameters as_dict * Update exp_UKF.py * Remove fit_keys * Update comment * Fixes for notebooks * Check stricter of relative/absolute tolerance * Change Thevenin parameters * Update ground truth * Pass parameters to get_data * Update max values to help GradientDescent * Reset Thevenin test but with new C1 * Equalise integration test settings * Update default SciPy DE tol * Relax atol a little * Move sigma0, bounds to BaseOptimiser, update likelihoods * Remove unused bounds_for_scipy * Update from list to args * Fix test_model_experiment_changes * style: pre-commit fixes * Add ValueError for None sigma * Rename add_parameter as add * Rename remove_parameter to remove * Apply suggestions from code review Co-authored-by: Brady Planden <55357039+BradyPlanden@users.noreply.github.com> * Update error match * Update tolerances and threshold * Rename update_bounds to set_bounds * Update all_samples to list * Update set_bounds to get_bounds * Move check on sigma * Update description --------- Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com> Co-authored-by: Brady Planden <55357039+BradyPlanden@users.noreply.github.com> --- .../1-single-pulse-circuit-model.ipynb | 6 +- .../equivalent_circuit_identification.ipynb | 8 +- .../multi_model_identification.ipynb | 4 +- .../multi_optimiser_identification.ipynb | 4 +- .../notebooks/optimiser_calibration.ipynb | 4 +- examples/notebooks/optimiser_interface.ipynb | 12 +- .../notebooks/pouch_cell_identification.ipynb | 4 +- examples/notebooks/spm_electrode_design.ipynb | 4 +- examples/scripts/BPX_spm.py | 12 +- examples/scripts/ecm_CMAES.py | 4 +- examples/scripts/exp_UKF.py | 26 +- examples/scripts/gitt.py | 14 +- examples/scripts/spm_CMAES.py | 12 +- examples/scripts/spm_IRPropMin.py | 4 +- examples/scripts/spm_MAP.py | 4 +- examples/scripts/spm_MLE.py | 4 +- examples/scripts/spm_NelderMead.py | 4 +- examples/scripts/spm_SNES.py | 4 +- examples/scripts/spm_UKF.py | 4 +- examples/scripts/spm_XNES.py | 4 +- examples/scripts/spm_adam.py | 4 +- examples/scripts/spm_descent.py | 4 +- examples/scripts/spm_pso.py | 4 +- examples/scripts/spm_scipymin.py | 4 +- examples/scripts/spme_max_energy.py | 4 +- examples/standalone/cost.py | 31 +-- examples/standalone/problem.py | 3 +- pybop/__init__.py | 2 +- pybop/_dataset.py | 4 +- pybop/costs/_likelihoods.py | 73 +++--- pybop/costs/base_cost.py | 30 +-- pybop/costs/design_costs.py | 14 +- pybop/costs/fitting_costs.py | 2 +- pybop/models/base_model.py | 29 ++- pybop/observers/observer.py | 17 +- pybop/observers/unscented_kalman.py | 5 +- pybop/optimisers/base_optimiser.py | 13 +- pybop/optimisers/base_pints_optimiser.py | 83 ++++-- pybop/optimisers/scipy_optimisers.py | 13 +- pybop/parameters/parameter.py | 241 ++++++++++++++++++ pybop/plotting/plot2d.py | 36 +-- pybop/plotting/plot_parameters.py | 7 +- pybop/problems/base_problem.py | 97 ++----- pybop/problems/design_problem.py | 17 +- pybop/problems/fitting_problem.py | 14 +- pyproject.toml | 1 - .../test_model_experiment_changes.py | 39 +-- .../integration/test_optimisation_options.py | 23 +- .../integration/test_spm_parameterisations.py | 65 ++--- .../test_thevenin_parameterisation.py | 37 ++- tests/unit/test_cost.py | 83 +++--- tests/unit/test_dataset.py | 8 +- tests/unit/test_likelihoods.py | 65 +++-- tests/unit/test_models.py | 19 +- tests/unit/test_observer_unscented_kalman.py | 18 +- tests/unit/test_observers.py | 22 +- tests/unit/test_optimisation.py | 42 +-- tests/unit/test_parameter_sets.py | 4 +- tests/unit/test_parameters.py | 90 ++++++- tests/unit/test_plots.py | 4 +- tests/unit/test_problem.py | 37 +-- tests/unit/test_standalone.py | 8 +- 62 files changed, 834 insertions(+), 624 deletions(-) diff --git a/examples/notebooks/LG_M50_ECM/1-single-pulse-circuit-model.ipynb b/examples/notebooks/LG_M50_ECM/1-single-pulse-circuit-model.ipynb index 3181d952..365eb6e1 100644 --- a/examples/notebooks/LG_M50_ECM/1-single-pulse-circuit-model.ipynb +++ b/examples/notebooks/LG_M50_ECM/1-single-pulse-circuit-model.ipynb @@ -1467,7 +1467,7 @@ "id": "bf63b4f9-de38-4e70-9472-1de4973a0954", "metadata": {}, "source": [ - "In this example, we are going to try to fit all five parameters at once. To do this, we define a `pybop.parameter` for each fitting parameter as," + "In this example, we are going to try to fit all five parameters at once. To do this, we define a `pybop.Parameter` for each fitting parameter and compile them in pybop.Parameters," ] }, { @@ -1484,7 +1484,7 @@ }, "outputs": [], "source": [ - "parameters = [\n", + "parameters = pybop.Parameters(\n", " pybop.Parameter(\n", " \"R0 [Ohm]\",\n", " prior=pybop.Gaussian(0.005, 0.0001),\n", @@ -1510,7 +1510,7 @@ " prior=pybop.Gaussian(3000, 2500),\n", " bounds=[0.5, 1e4],\n", " ),\n", - "]" + ")" ] }, { diff --git a/examples/notebooks/equivalent_circuit_identification.ipynb b/examples/notebooks/equivalent_circuit_identification.ipynb index 95511e84..8a13a199 100644 --- a/examples/notebooks/equivalent_circuit_identification.ipynb +++ b/examples/notebooks/equivalent_circuit_identification.ipynb @@ -239,7 +239,7 @@ "id": "bf63b4f9-de38-4e70-9472-1de4973a0954", "metadata": {}, "source": [ - "In this example, we are going to try to fit all five parameters at once. This isn't recommend for real-life application as identifiablity is challenging to guarantee with this large a parameter space. To do this, we define the `pybop.parameters` as," + "In this example, we are going to try to fit all five parameters at once. This isn't recommend for real-life application as identifiablity is challenging to guarantee with this large a parameter space. To do this, we define the `pybop.Parameters` as," ] }, { @@ -256,7 +256,7 @@ }, "outputs": [], "source": [ - "parameters = [\n", + "parameters = pybop.Parameters(\n", " pybop.Parameter(\n", " \"R0 [Ohm]\",\n", " prior=pybop.Gaussian(0.0002, 0.0001),\n", @@ -278,11 +278,11 @@ " bounds=[2500, 5e4],\n", " ),\n", " pybop.Parameter(\n", - " \"C1 [F]\",\n", + " \"C2 [F]\",\n", " prior=pybop.Gaussian(10000, 2500),\n", " bounds=[2500, 5e4],\n", " ),\n", - "]" + ")" ] }, { diff --git a/examples/notebooks/multi_model_identification.ipynb b/examples/notebooks/multi_model_identification.ipynb index f2ab1822..699b2eda 100644 --- a/examples/notebooks/multi_model_identification.ipynb +++ b/examples/notebooks/multi_model_identification.ipynb @@ -3729,7 +3729,7 @@ }, "outputs": [], "source": [ - "parameters = [\n", + "parameters = pybop.Parameters(\n", " pybop.Parameter(\n", " \"Positive electrode thickness [m]\",\n", " prior=pybop.Gaussian(7.56e-05, 0.05e-05),\n", @@ -3742,7 +3742,7 @@ " bounds=[65e-06, 90e-06],\n", " true_value=parameter_set[\"Negative electrode thickness [m]\"],\n", " ),\n", - "]" + ")" ] }, { diff --git a/examples/notebooks/multi_optimiser_identification.ipynb b/examples/notebooks/multi_optimiser_identification.ipynb index 5c38e963..8bb6f353 100644 --- a/examples/notebooks/multi_optimiser_identification.ipynb +++ b/examples/notebooks/multi_optimiser_identification.ipynb @@ -259,7 +259,7 @@ }, "outputs": [], "source": [ - "parameters = [\n", + "parameters = pybop.Parameters(\n", " pybop.Parameter(\n", " \"Negative electrode active material volume fraction\",\n", " prior=pybop.Gaussian(0.6, 0.02),\n", @@ -270,7 +270,7 @@ " prior=pybop.Gaussian(0.48, 0.02),\n", " bounds=[0.4, 0.7],\n", " ),\n", - "]" + ")" ] }, { diff --git a/examples/notebooks/optimiser_calibration.ipynb b/examples/notebooks/optimiser_calibration.ipynb index f94ecec6..3199fadb 100644 --- a/examples/notebooks/optimiser_calibration.ipynb +++ b/examples/notebooks/optimiser_calibration.ipynb @@ -232,7 +232,7 @@ }, "outputs": [], "source": [ - "parameters = [\n", + "parameters = pybop.Parameters(\n", " pybop.Parameter(\n", " \"Negative electrode active material volume fraction\",\n", " prior=pybop.Gaussian(0.7, 0.025),\n", @@ -243,7 +243,7 @@ " prior=pybop.Gaussian(0.6, 0.025),\n", " bounds=[0.5, 0.8],\n", " ),\n", - "]" + ")" ] }, { diff --git a/examples/notebooks/optimiser_interface.ipynb b/examples/notebooks/optimiser_interface.ipynb index 1200e399..efe4b71f 100644 --- a/examples/notebooks/optimiser_interface.ipynb +++ b/examples/notebooks/optimiser_interface.ipynb @@ -100,13 +100,11 @@ ")\n", "\n", "# Define the parameters\n", - "parameters = [\n", - " pybop.Parameter(\n", - " \"R0 [Ohm]\",\n", - " prior=pybop.Gaussian(0.0002, 0.0001),\n", - " bounds=[1e-4, 1e-2],\n", - " )\n", - "]\n", + "parameters = pybop.Parameter(\n", + " \"R0 [Ohm]\",\n", + " prior=pybop.Gaussian(0.0002, 0.0001),\n", + " bounds=[1e-4, 1e-2],\n", + ")\n", "\n", "# Generate synthetic data\n", "t_eval = np.arange(0, 900, 2)\n", diff --git a/examples/notebooks/pouch_cell_identification.ipynb b/examples/notebooks/pouch_cell_identification.ipynb index 0cf64bf5..c24300ea 100644 --- a/examples/notebooks/pouch_cell_identification.ipynb +++ b/examples/notebooks/pouch_cell_identification.ipynb @@ -323,7 +323,7 @@ }, "outputs": [], "source": [ - "parameters = [\n", + "parameters = pybop.Parameters(\n", " pybop.Parameter(\n", " \"Negative electrode active material volume fraction\",\n", " prior=pybop.Gaussian(0.7, 0.05),\n", @@ -334,7 +334,7 @@ " prior=pybop.Gaussian(0.58, 0.05),\n", " bounds=[0.5, 0.8],\n", " ),\n", - "]" + ")" ] }, { diff --git a/examples/notebooks/spm_electrode_design.ipynb b/examples/notebooks/spm_electrode_design.ipynb index fd537352..950cee32 100644 --- a/examples/notebooks/spm_electrode_design.ipynb +++ b/examples/notebooks/spm_electrode_design.ipynb @@ -143,7 +143,7 @@ }, "outputs": [], "source": [ - "parameters = [\n", + "parameters = pybop.Parameters(\n", " pybop.Parameter(\n", " \"Positive electrode thickness [m]\",\n", " prior=pybop.Gaussian(7.56e-05, 0.05e-05),\n", @@ -154,7 +154,7 @@ " prior=pybop.Gaussian(5.22e-06, 0.05e-06),\n", " bounds=[2e-06, 9e-06],\n", " ),\n", - "]" + ")" ] }, { diff --git a/examples/scripts/BPX_spm.py b/examples/scripts/BPX_spm.py index 0d108935..6fdb7649 100644 --- a/examples/scripts/BPX_spm.py +++ b/examples/scripts/BPX_spm.py @@ -10,7 +10,7 @@ model = pybop.lithium_ion.SPM(parameter_set=parameter_set) # Fitting parameters -parameters = [ +parameters = pybop.Parameters( pybop.Parameter( "Negative particle radius [m]", prior=pybop.Gaussian(6e-06, 0.1e-6), @@ -23,7 +23,7 @@ bounds=[1e-7, 9e-7], true_value=parameter_set["Positive particle radius [m]"], ), -] +) # Generate data sigma = 0.001 @@ -47,13 +47,7 @@ # Run the optimisation x, final_cost = optim.run() -print( - "True parameters:", - [ - parameters[0].true_value, - parameters[1].true_value, - ], -) +print("True parameters:", parameters.true_value()) print("Estimated parameters:", x) # Plot the timeseries output diff --git a/examples/scripts/ecm_CMAES.py b/examples/scripts/ecm_CMAES.py index 5d241a6b..fc711cab 100644 --- a/examples/scripts/ecm_CMAES.py +++ b/examples/scripts/ecm_CMAES.py @@ -45,7 +45,7 @@ ) # Fitting parameters -parameters = [ +parameters = pybop.Parameters( pybop.Parameter( "R0 [Ohm]", prior=pybop.Gaussian(0.0002, 0.0001), @@ -56,7 +56,7 @@ prior=pybop.Gaussian(0.0001, 0.0001), bounds=[1e-5, 1e-2], ), -] +) sigma = 0.001 t_eval = np.arange(0, 900, 3) diff --git a/examples/scripts/exp_UKF.py b/examples/scripts/exp_UKF.py index 5a61436b..d469c781 100644 --- a/examples/scripts/exp_UKF.py +++ b/examples/scripts/exp_UKF.py @@ -7,42 +7,41 @@ # Parameter set and model definition parameter_set = pybamm.ParameterValues({"k": "[input]", "y0": "[input]"}) model = ExponentialDecay(parameter_set=parameter_set, n_states=1) -x0 = np.array([0.1, 1.0]) # Fitting parameters -parameters = [ +parameters = pybop.Parameters( pybop.Parameter( "k", prior=pybop.Gaussian(0.1, 0.05), bounds=[0, 1], + true_value=0.1, ), pybop.Parameter( "y0", prior=pybop.Gaussian(1, 0.05), bounds=[0, 3], + true_value=1.0, ), -] - -# Verification: save fixed inputs for testing -inputs = dict() -for i, param in enumerate(parameters): - inputs[param.name] = x0[i] +) # Make a prediction with measurement noise sigma = 1e-2 t_eval = np.linspace(0, 20, 10) -values = model.predict(t_eval=t_eval, inputs=inputs) +model.parameters = parameters +values = model.predict(t_eval=t_eval, inputs=parameters.true_value()) values = values["2y"].data corrupt_values = values + np.random.normal(0, sigma, len(t_eval)) # Verification step: compute the analytical solution for 2y -expected_values = 2 * inputs["y0"] * np.exp(-inputs["k"] * t_eval) +expected_values = ( + 2 * parameters["y0"].true_value * np.exp(-parameters["k"].true_value * t_eval) +) # Verification step: make another prediction using the Observer class model.build(parameters=parameters) -simulator = pybop.Observer(parameters, model, signal=["2y"], x0=x0) +simulator = pybop.Observer(parameters, model, signal=["2y"]) simulator._time_data = t_eval -measurements = simulator.evaluate(x0) +measurements = simulator.evaluate(parameters.true_value()) # Verification step: Compare by plotting go = pybop.PlotlyManager().go @@ -82,11 +81,10 @@ measurement_noise, dataset, signal=signal, - x0=x0, ) # Verification step: Find the maximum likelihood estimate given the true parameters -estimation = observer.evaluate(x0) +estimation = observer.evaluate(parameters.true_value()) # Verification step: Add the estimate to the plot line4 = go.Scatter( diff --git a/examples/scripts/gitt.py b/examples/scripts/gitt.py index b81f7433..52517fdb 100644 --- a/examples/scripts/gitt.py +++ b/examples/scripts/gitt.py @@ -35,14 +35,12 @@ # Define the cost to optimise model = pybop.lithium_ion.WeppnerHuggins(parameter_set=parameter_set) -parameters = [ - pybop.Parameter( - "Positive electrode diffusivity [m2.s-1]", - prior=pybop.Gaussian(5e-14, 1e-13), - bounds=[1e-16, 1e-11], - true_value=parameter_set["Positive electrode diffusivity [m2.s-1]"], - ), -] +parameters = pybop.Parameter( + "Positive electrode diffusivity [m2.s-1]", + prior=pybop.Gaussian(5e-14, 1e-13), + bounds=[1e-16, 1e-11], + true_value=parameter_set["Positive electrode diffusivity [m2.s-1]"], +) problem = pybop.FittingProblem( model, diff --git a/examples/scripts/spm_CMAES.py b/examples/scripts/spm_CMAES.py index b60bc019..1fc051cc 100644 --- a/examples/scripts/spm_CMAES.py +++ b/examples/scripts/spm_CMAES.py @@ -7,7 +7,7 @@ model = pybop.lithium_ion.SPM(parameter_set=parameter_set) # Fitting parameters -parameters = [ +parameters = pybop.Parameters( pybop.Parameter( "Negative particle radius [m]", prior=pybop.Gaussian(6e-06, 0.1e-6), @@ -20,7 +20,7 @@ bounds=[1e-6, 9e-6], true_value=parameter_set["Positive particle radius [m]"], ), -] +) # Generate data sigma = 0.001 @@ -46,13 +46,7 @@ # Run the optimisation x, final_cost = optim.run() -print( - "True parameters:", - [ - parameters[0].true_value, - parameters[1].true_value, - ], -) +print("True parameters:", parameters.true_value()) print("Estimated parameters:", x) # Plot the time series diff --git a/examples/scripts/spm_IRPropMin.py b/examples/scripts/spm_IRPropMin.py index 2aedff90..3b38668c 100644 --- a/examples/scripts/spm_IRPropMin.py +++ b/examples/scripts/spm_IRPropMin.py @@ -7,7 +7,7 @@ model = pybop.lithium_ion.SPM(parameter_set=parameter_set) # Fitting parameters -parameters = [ +parameters = pybop.Parameters( pybop.Parameter( "Negative electrode active material volume fraction", prior=pybop.Gaussian(0.6, 0.05), @@ -16,7 +16,7 @@ "Positive electrode active material volume fraction", prior=pybop.Gaussian(0.48, 0.05), ), -] +) # Generate data sigma = 0.001 diff --git a/examples/scripts/spm_MAP.py b/examples/scripts/spm_MAP.py index e09ce231..191f93d8 100644 --- a/examples/scripts/spm_MAP.py +++ b/examples/scripts/spm_MAP.py @@ -7,7 +7,7 @@ model = pybop.lithium_ion.SPM(parameter_set=parameter_set) # Fitting parameters -parameters = [ +parameters = pybop.Parameters( pybop.Parameter( "Negative electrode active material volume fraction", prior=pybop.Gaussian(0.6, 0.05), @@ -18,7 +18,7 @@ prior=pybop.Gaussian(0.48, 0.05), bounds=[0.4, 0.7], ), -] +) # Set initial parameter values parameter_set.update( diff --git a/examples/scripts/spm_MLE.py b/examples/scripts/spm_MLE.py index 9a3636de..7e1b3c93 100644 --- a/examples/scripts/spm_MLE.py +++ b/examples/scripts/spm_MLE.py @@ -7,7 +7,7 @@ model = pybop.lithium_ion.SPM(parameter_set=parameter_set) # Fitting parameters -parameters = [ +parameters = pybop.Parameters( pybop.Parameter( "Negative electrode active material volume fraction", prior=pybop.Gaussian(0.6, 0.05), @@ -18,7 +18,7 @@ prior=pybop.Gaussian(0.48, 0.05), bounds=[0.4, 0.7], ), -] +) # Set initial parameter values parameter_set.update( diff --git a/examples/scripts/spm_NelderMead.py b/examples/scripts/spm_NelderMead.py index 90a95bff..82639632 100644 --- a/examples/scripts/spm_NelderMead.py +++ b/examples/scripts/spm_NelderMead.py @@ -7,7 +7,7 @@ model = pybop.lithium_ion.SPM(parameter_set=parameter_set) # Fitting parameters -parameters = [ +parameters = pybop.Parameters( pybop.Parameter( "Negative electrode active material volume fraction", prior=pybop.Gaussian(0.68, 0.05), @@ -16,7 +16,7 @@ "Positive electrode active material volume fraction", prior=pybop.Gaussian(0.58, 0.05), ), -] +) # Generate data init_soc = 0.5 diff --git a/examples/scripts/spm_SNES.py b/examples/scripts/spm_SNES.py index a304dd04..d2afcc85 100644 --- a/examples/scripts/spm_SNES.py +++ b/examples/scripts/spm_SNES.py @@ -7,7 +7,7 @@ model = pybop.lithium_ion.SPM(parameter_set=parameter_set) # Fitting parameters -parameters = [ +parameters = pybop.Parameters( pybop.Parameter( "Negative electrode active material volume fraction", prior=pybop.Gaussian(0.6, 0.05), @@ -18,7 +18,7 @@ prior=pybop.Gaussian(0.48, 0.05), bounds=[0.4, 0.7], ), -] +) sigma = 0.001 t_eval = np.arange(0, 900, 3) diff --git a/examples/scripts/spm_UKF.py b/examples/scripts/spm_UKF.py index 0814a22c..e9972bd0 100644 --- a/examples/scripts/spm_UKF.py +++ b/examples/scripts/spm_UKF.py @@ -7,7 +7,7 @@ model = pybop.lithium_ion.SPM(parameter_set=parameter_set) # Fitting parameters -parameters = [ +parameters = pybop.Parameters( pybop.Parameter( "Negative electrode active material volume fraction", prior=pybop.Gaussian(0.6, 0.05), @@ -18,7 +18,7 @@ prior=pybop.Gaussian(0.48, 0.05), bounds=[0.4, 0.7], ), -] +) # Make a prediction with measurement noise sigma = 0.001 diff --git a/examples/scripts/spm_XNES.py b/examples/scripts/spm_XNES.py index bcca73de..59b6eca8 100644 --- a/examples/scripts/spm_XNES.py +++ b/examples/scripts/spm_XNES.py @@ -7,7 +7,7 @@ model = pybop.lithium_ion.SPM(parameter_set=parameter_set) # Fitting parameters -parameters = [ +parameters = pybop.Parameters( pybop.Parameter( "Negative electrode active material volume fraction", prior=pybop.Gaussian(0.68, 0.05), @@ -18,7 +18,7 @@ prior=pybop.Gaussian(0.58, 0.05), bounds=[0.4, 0.7], ), -] +) sigma = 0.001 t_eval = np.arange(0, 900, 3) diff --git a/examples/scripts/spm_adam.py b/examples/scripts/spm_adam.py index 7523384a..82f884e2 100644 --- a/examples/scripts/spm_adam.py +++ b/examples/scripts/spm_adam.py @@ -7,7 +7,7 @@ model = pybop.lithium_ion.SPM(parameter_set=parameter_set) # Fitting parameters -parameters = [ +parameters = pybop.Parameters( pybop.Parameter( "Negative electrode active material volume fraction", prior=pybop.Gaussian(0.68, 0.05), @@ -16,7 +16,7 @@ "Positive electrode active material volume fraction", prior=pybop.Gaussian(0.58, 0.05), ), -] +) # Generate data init_soc = 0.5 diff --git a/examples/scripts/spm_descent.py b/examples/scripts/spm_descent.py index 9bacfdc8..df57a7ca 100644 --- a/examples/scripts/spm_descent.py +++ b/examples/scripts/spm_descent.py @@ -7,7 +7,7 @@ model = pybop.lithium_ion.SPM(parameter_set=parameter_set) # Fitting parameters -parameters = [ +parameters = pybop.Parameters( pybop.Parameter( "Negative electrode active material volume fraction", prior=pybop.Gaussian(0.68, 0.05), @@ -16,7 +16,7 @@ "Positive electrode active material volume fraction", prior=pybop.Gaussian(0.58, 0.05), ), -] +) # Generate data sigma = 0.001 diff --git a/examples/scripts/spm_pso.py b/examples/scripts/spm_pso.py index 4b99bd12..acb3e1c6 100644 --- a/examples/scripts/spm_pso.py +++ b/examples/scripts/spm_pso.py @@ -7,7 +7,7 @@ model = pybop.lithium_ion.SPM(parameter_set=parameter_set) # Fitting parameters -parameters = [ +parameters = pybop.Parameters( pybop.Parameter( "Negative electrode active material volume fraction", prior=pybop.Gaussian(0.6, 0.05), @@ -18,7 +18,7 @@ prior=pybop.Gaussian(0.48, 0.05), bounds=[0.4, 0.7], ), -] +) sigma = 0.001 t_eval = np.arange(0, 900, 3) diff --git a/examples/scripts/spm_scipymin.py b/examples/scripts/spm_scipymin.py index 12a15234..8c7b80c5 100644 --- a/examples/scripts/spm_scipymin.py +++ b/examples/scripts/spm_scipymin.py @@ -19,7 +19,7 @@ ) # Fitting parameters -parameters = [ +parameters = pybop.Parameters( pybop.Parameter( "Negative electrode active material volume fraction", prior=pybop.Gaussian(0.6, 0.05), @@ -30,7 +30,7 @@ prior=pybop.Gaussian(0.48, 0.05), bounds=[0.4, 0.7], ), -] +) # Define the cost to optimise signal = ["Voltage [V]"] diff --git a/examples/scripts/spme_max_energy.py b/examples/scripts/spme_max_energy.py index 8590f68c..800a535c 100644 --- a/examples/scripts/spme_max_energy.py +++ b/examples/scripts/spme_max_energy.py @@ -20,7 +20,7 @@ model = pybop.lithium_ion.SPMe(parameter_set=parameter_set) # Fitting parameters -parameters = [ +parameters = pybop.Parameters( pybop.Parameter( "Positive electrode thickness [m]", prior=pybop.Gaussian(7.56e-05, 0.1e-05), @@ -31,7 +31,7 @@ prior=pybop.Gaussian(5.22e-06, 0.1e-06), bounds=[2e-06, 9e-06], ), -] +) # Define test protocol experiment = pybop.Experiment( diff --git a/examples/standalone/cost.py b/examples/standalone/cost.py index c3763ff4..806bc0ea 100644 --- a/examples/standalone/cost.py +++ b/examples/standalone/cost.py @@ -1,5 +1,3 @@ -import numpy as np - import pybop @@ -17,13 +15,9 @@ class StandaloneCost(pybop.BaseCost): A dummy problem instance used to initialize the superclass. This is not used in the current class but is accepted for compatibility with the BaseCost interface. - x0 : array-like - The initial guess for the optimization problem, set to [4.2]. - _n_parameters : int - The number of parameters in the model, which is 1 in this case. - bounds : dict - A dictionary containing the lower and upper bounds for the parameter, - set to [-1] and [10], respectively. + parameters : pybop.Parameters + A pybop.Parameters object storing a dictionary of parameters and their + properties, for example their initial value and bounds. Methods ------- @@ -33,20 +27,21 @@ class StandaloneCost(pybop.BaseCost): def __init__(self, problem=None): """ - Initialize the StandaloneCost class with optional problem instance. + Initialise the StandaloneCost class with optional problem instance. - The problem parameter is not utilized in this subclass. The initial guess, - number of parameters, and bounds are predefined for the standalone cost function. + The problem object is not utilised in this subclass. The parameters, including + their initial value and bounds, are defined within this standalone cost object. """ super().__init__(problem) - self.x0 = np.array([4.2]) - self._n_parameters = len(self.x0) - - self.bounds = dict( - lower=[-1], - upper=[10], + self.parameters = pybop.Parameters( + pybop.Parameter( + "x", + initial_value=4.2, + bounds=[-1, 10], + ), ) + self.x0 = self.parameters.initial_value() def _evaluate(self, x, grad=None): """ diff --git a/examples/standalone/problem.py b/examples/standalone/problem.py index 6fdcf97f..d6d1f4b0 100644 --- a/examples/standalone/problem.py +++ b/examples/standalone/problem.py @@ -17,10 +17,9 @@ def __init__( signal=None, additional_variables=None, init_soc=None, - x0=None, ): super().__init__( - parameters, model, check_model, signal, additional_variables, init_soc, x0 + parameters, model, check_model, signal, additional_variables, init_soc ) self._dataset = dataset.data diff --git a/pybop/__init__.py b/pybop/__init__.py index ecd42019..14c316df 100644 --- a/pybop/__init__.py +++ b/pybop/__init__.py @@ -117,7 +117,7 @@ # # Parameter classes # -from .parameters.parameter import Parameter +from .parameters.parameter import Parameter, Parameters from .parameters.parameter_set import ParameterSet from .parameters.priors import BasePrior, Gaussian, Uniform, Exponential diff --git a/pybop/_dataset.py b/pybop/_dataset.py index 7d18e6fe..310cc744 100644 --- a/pybop/_dataset.py +++ b/pybop/_dataset.py @@ -23,7 +23,7 @@ def __init__(self, data_dictionary): if isinstance(data_dictionary, pybamm.solvers.solution.Solution): data_dictionary = data_dictionary.get_data_dict() if not isinstance(data_dictionary, dict): - raise ValueError("The input to pybop.Dataset must be a dictionary.") + raise TypeError("The input to pybop.Dataset must be a dictionary.") self.data = data_dictionary self.names = self.data.keys() @@ -64,7 +64,7 @@ def __getitem__(self, key): Returns ------- list or np.ndarray - The data series corresonding to the key. + The data series corresponding to the key. Raises ------ diff --git a/pybop/costs/_likelihoods.py b/pybop/costs/_likelihoods.py index 91374cc0..cd5e4a9c 100644 --- a/pybop/costs/_likelihoods.py +++ b/pybop/costs/_likelihoods.py @@ -8,14 +8,43 @@ class BaseLikelihood(BaseCost): Base class for likelihoods """ - def __init__(self, problem, sigma=None): - super(BaseLikelihood, self).__init__(problem, sigma) + def __init__(self, problem): + super(BaseLikelihood, self).__init__(problem) self.n_time_data = problem.n_time_data + +class GaussianLogLikelihoodKnownSigma(BaseLikelihood): + """ + This class represents a Gaussian Log Likelihood with a known sigma, + which assumes that the data follows a Gaussian distribution and computes + the log-likelihood of observed data under this assumption. + + Parameters + ---------- + sigma : scalar or array + Initial standard deviation around ``x0``. Either a scalar value (one + standard deviation for all coordinates) or an array with one entry + per dimension. Not all methods will use this information. + """ + + def __init__(self, problem, sigma): + super(GaussianLogLikelihoodKnownSigma, self).__init__(problem) + self.sigma = None + self.set_sigma(sigma) + self._offset = -0.5 * self.n_time_data * np.log(2 * np.pi / self.sigma) + self._multip = -1 / (2.0 * self.sigma**2) + self.sigma2 = self.sigma**-2 + self._dl = np.ones(self.n_parameters) + def set_sigma(self, sigma): """ Setter for sigma parameter """ + if sigma is None: + raise ValueError( + "The GaussianLogLikelihoodKnownSigma cost requires sigma to be " + + "either a scalar value or an array with one entry per dimension." + ) if not isinstance(sigma, np.ndarray): sigma = np.array(sigma) @@ -24,41 +53,15 @@ def set_sigma(self, sigma): raise ValueError("Sigma must contain only numeric values") if np.any(sigma <= 0): - raise ValueError("Sigma must not be negative") + raise ValueError("Sigma must be positive") else: - self.sigma0 = sigma + self.sigma = sigma def get_sigma(self): """ Getter for sigma parameter """ - return self.sigma0 - - def get_n_parameters(self): - """ - Returns the number of parameters - """ - return self._n_parameters - - -class GaussianLogLikelihoodKnownSigma(BaseLikelihood): - """ - This class represents a Gaussian Log Likelihood with a known sigma, - which assumes that the data follows a Gaussian distribution and computes - the log-likelihood of observed data under this assumption. - - Attributes: - _logpi (float): Precomputed offset value for the log-likelihood function. - """ - - def __init__(self, problem, sigma): - super(GaussianLogLikelihoodKnownSigma, self).__init__(problem, sigma) - if sigma is not None: - self.set_sigma(sigma) - self._offset = -0.5 * self.n_time_data * np.log(2 * np.pi / self.sigma0) - self._multip = -1 / (2.0 * self.sigma0**2) - self.sigma2 = self.sigma0**-2 - self._dl = np.ones(self._n_parameters) + return self.sigma def _evaluate(self, x, grad=None): """ @@ -111,14 +114,16 @@ class GaussianLogLikelihood(BaseLikelihood): data follows a Gaussian distribution and computes the log-likelihood of observed data under this assumption. - Attributes: - _logpi (float): Precomputed offset value for the log-likelihood function. + Attributes + ---------- + _logpi : float + Precomputed offset value for the log-likelihood function. """ def __init__(self, problem): super(GaussianLogLikelihood, self).__init__(problem) self._logpi = -0.5 * self.n_time_data * np.log(2 * np.pi) - self._dl = np.ones(self._n_parameters + self.n_outputs) + self._dl = np.ones(self.n_parameters + self.n_outputs) def _evaluate(self, x, grad=None): """ diff --git a/pybop/costs/base_cost.py b/pybop/costs/base_cost.py index 025cb2e4..04d0a393 100644 --- a/pybop/costs/base_cost.py +++ b/pybop/costs/base_cost.py @@ -1,5 +1,3 @@ -import numpy as np - from pybop import BaseProblem @@ -21,36 +19,24 @@ class BaseCost: An array containing the target data to fit. x0 : array-like The initial guess for the model parameters. - bounds : tuple - The bounds for the model parameters. - sigma0 : scalar or array - Initial standard deviation around ``x0``. Either a scalar value (one - standard deviation for all coordinates) or an array with one entry - per dimension. Not all methods will use this information. - _n_parameters : int - The number of parameters in the model. n_outputs : int The number of outputs in the model. """ - def __init__(self, problem=None, sigma=None): + def __init__(self, problem=None): + self.parameters = None self.problem = problem self.x0 = None - self.bounds = None - self.sigma0 = sigma if isinstance(self.problem, BaseProblem): - self._target = problem._target - self.parameters = problem.parameters - self.x0 = problem.x0 - self.bounds = problem.bounds - self.n_outputs = problem.n_outputs - self.signal = problem.signal - self._n_parameters = problem.n_parameters - self.sigma0 = sigma or problem.sigma0 or np.zeros(self._n_parameters) + self._target = self.problem._target + self.parameters = self.problem.parameters + self.x0 = self.problem.x0 + self.n_outputs = self.problem.n_outputs + self.signal = self.problem.signal @property def n_parameters(self): - return self._n_parameters + return len(self.parameters) def __call__(self, x, grad=None): """ diff --git a/pybop/costs/design_costs.py b/pybop/costs/design_costs.py index f6364cdc..60064c65 100644 --- a/pybop/costs/design_costs.py +++ b/pybop/costs/design_costs.py @@ -44,20 +44,20 @@ def __init__(self, problem, update_capacity=False): warnings.warn(nominal_capacity_warning, UserWarning) self.update_capacity = update_capacity self.parameter_set = problem.model.parameter_set - self.update_simulation_data(problem.x0) + self.update_simulation_data(self.x0) - def update_simulation_data(self, initial_conditions): + def update_simulation_data(self, x0): """ - Updates the simulation data based on the initial conditions. + Updates the simulation data based on the initial parameter values. Parameters ---------- - initial_conditions : array - The initial conditions for the simulation. + x0 : array + The initial parameter values for the simulation. """ if self.update_capacity: - self.problem.model.approximate_capacity(self.problem.x0) - solution = self.problem.evaluate(initial_conditions) + self.problem.model.approximate_capacity(x0) + solution = self.problem.evaluate(x0) if "Time [s]" not in solution: raise ValueError("The solution does not contain time data.") diff --git a/pybop/costs/fitting_costs.py b/pybop/costs/fitting_costs.py index b7b26659..03b4488a 100644 --- a/pybop/costs/fitting_costs.py +++ b/pybop/costs/fitting_costs.py @@ -253,7 +253,7 @@ def _evaluate(self, x, grad=None): float The observer cost (negative of the log likelihood). """ - inputs = {key: x[i] for i, key in enumerate(self._observer._model.fit_keys)} + inputs = self._observer.parameters.as_dict(x) log_likelihood = self._observer.log_likelihood( self._target, self._observer.time_data(), inputs ) diff --git a/pybop/models/base_model.py b/pybop/models/base_model.py index 729c826a..4b40386a 100644 --- a/pybop/models/base_model.py +++ b/pybop/models/base_model.py @@ -70,10 +70,13 @@ def __init__(self, name="Base Model", parameter_set=None): self.additional_variables = [] self.matched_parameters = {} self.non_matched_parameters = {} - self.fit_keys = [] self.param_check_counter = 0 self.allow_infeasible_solutions = True + @property + def n_parameters(self): + return len(self.parameters) + def build( self, dataset=None, @@ -103,9 +106,6 @@ def build( self.parameters = parameters if self.parameters is not None: self.classify_and_update_parameters(self.parameters) - self.fit_keys = [param.name for param in self.parameters] - else: - self.fit_keys = [] if init_soc is not None: self.set_init_soc(init_soc) @@ -173,7 +173,10 @@ def set_params(self, rebuild=False): self._parameter_set[key] = "[input]" if self.dataset is not None and (not self.matched_parameters or not rebuild): - if "Current function [A]" not in self.fit_keys: + if ( + self.parameters is None + or "Current function [A]" not in self.parameters.keys() + ): self._parameter_set["Current function [A]"] = pybamm.Interpolant( self.dataset["Time [s]"], self.dataset["Current function [A]"], @@ -219,8 +222,8 @@ def rebuild( The initial state of charge to be used in simulations. """ self.dataset = dataset - self.parameters = parameters if parameters is not None: + self.parameters = parameters self.classify_and_update_parameters(parameters) if init_soc is not None: @@ -250,7 +253,7 @@ def classify_and_update_parameters(self, parameters): parameters : pybop.ParameterSet """ - parameter_dictionary = {param.name: param.value for param in parameters} + parameter_dictionary = parameters.as_dict() matched_parameters = { param: parameter_dictionary[param] for param in parameter_dictionary @@ -280,7 +283,7 @@ def reinit( raise ValueError("Model must be built before calling reinit") if not isinstance(inputs, dict): - inputs = {key: inputs[i] for i, key in enumerate(self.fit_keys)} + inputs = self.parameters.as_dict(inputs) self._solver.set_up(self._built_model, inputs=inputs) @@ -349,7 +352,7 @@ def simulate(self, inputs, t_eval) -> np.ndarray[np.float64]: else: if not isinstance(inputs, dict): - inputs = {key: inputs[i] for i, key in enumerate(self.fit_keys)} + inputs = self.parameters.as_dict(inputs) if self.check_params( inputs=inputs, @@ -405,7 +408,7 @@ def simulateS1(self, inputs, t_eval): ) if not isinstance(inputs, dict): - inputs = {key: inputs[i] for i, key in enumerate(self.fit_keys)} + inputs = self.parameters.as_dict(inputs) if self.check_params( inputs=inputs, @@ -433,7 +436,7 @@ def simulateS1(self, inputs, t_eval): dy[:, i, :] = np.stack( [ sol[signal].sensitivities[key].toarray()[:, 0] - for key in self.fit_keys + for key in self.parameters.keys() ], axis=-1, ) @@ -498,7 +501,7 @@ def predict( parameter_set = parameter_set or self._unprocessed_parameter_set if inputs is not None: if not isinstance(inputs, dict): - inputs = {key: inputs[i] for i, key in enumerate(self.fit_keys)} + inputs = self.parameters.as_dict(inputs) parameter_set.update(inputs) if self.check_params( @@ -555,7 +558,7 @@ def check_params( + f" or None, but received a list with type: {type(inputs)}" ) else: - inputs = {key: inputs[i] for i, key in enumerate(self.fit_keys)} + inputs = self.parameters.as_dict(inputs) return self._check_params( inputs=inputs, allow_infeasible_solutions=allow_infeasible_solutions diff --git a/pybop/observers/observer.py b/pybop/observers/observer.py index 38c7ef86..162d03de 100644 --- a/pybop/observers/observer.py +++ b/pybop/observers/observer.py @@ -1,10 +1,10 @@ -from typing import List, Optional +from typing import Optional import numpy as np from pybop import BaseProblem from pybop.models.base_model import BaseModel, Inputs, TimeSeriesState -from pybop.parameters.parameter import Parameter +from pybop.parameters.parameter import Parameters class Observer(BaseProblem): @@ -16,8 +16,8 @@ class Observer(BaseProblem): Parameters ---------- - parameters : list - List of parameters for the problem. + parameters : pybop.Parameter or pybop.Parameters + An object or list of the parameters for the problem. model : BaseModel The model to observe. check_model : bool, optional @@ -28,8 +28,6 @@ class Observer(BaseProblem): Additional variables to observe and store in the solution (default: []). init_soc : float, optional Initial state of charge (default: None). - x0 : np.ndarray, optional - Initial parameter values (default: None). """ # define a subtype for covariance matrices for use by derived classes @@ -37,16 +35,15 @@ class Observer(BaseProblem): def __init__( self, - parameters: List[Parameter], + parameters: Parameters, model: BaseModel, check_model=True, signal=["Voltage [V]"], additional_variables=[], init_soc=None, - x0=None, ) -> None: super().__init__( - parameters, model, check_model, signal, additional_variables, init_soc, x0 + parameters, model, check_model, signal, additional_variables, init_soc ) if model._built_model is None: raise ValueError("Only built models can be used in Observers") @@ -160,7 +157,7 @@ def evaluate(self, x): The model output y(t) simulated with inputs x. """ inputs = dict() - if isinstance(x[0], Parameter): + if isinstance(x, Parameters): for param in x: inputs[param.name] = param.value else: # x is an array of parameter values diff --git a/pybop/observers/unscented_kalman.py b/pybop/observers/unscented_kalman.py index 562bf658..b7ea0f35 100644 --- a/pybop/observers/unscented_kalman.py +++ b/pybop/observers/unscented_kalman.py @@ -35,8 +35,6 @@ class UnscentedKalmanFilterObserver(Observer): The signal to observe. init_soc : float, optional Initial state of charge (default: None). - x0 : np.ndarray, optional - Initial parameter values (default: None). """ Covariance = np.ndarray @@ -53,10 +51,9 @@ def __init__( signal=["Voltage [V]"], additional_variables=[], init_soc=None, - x0=None, ) -> None: super().__init__( - parameters, model, check_model, signal, additional_variables, init_soc, x0 + parameters, model, check_model, signal, additional_variables, init_soc ) if dataset is not None: self._dataset = dataset.data diff --git a/pybop/optimisers/base_optimiser.py b/pybop/optimisers/base_optimiser.py index 1a600d7e..dfe60d36 100644 --- a/pybop/optimisers/base_optimiser.py +++ b/pybop/optimisers/base_optimiser.py @@ -28,7 +28,9 @@ class BaseOptimiser: bounds : dict Dictionary containing the parameter bounds with keys 'lower' and 'upper'. sigma0 : float or sequence - Initial step size or standard deviation for the optimiser. + Initial step size or standard deviation around ``x0``. Either a scalar value (one + standard deviation for all coordinates) or an array with one entry per dimension. + Not all methods will use this information. verbose : bool, optional If True, the optimisation progress is printed (default: False). minimising : bool, optional @@ -62,11 +64,16 @@ def __init__( if isinstance(cost, BaseCost): self.cost = cost self.x0 = cost.x0 - self.bounds = cost.bounds - self.sigma0 = cost.sigma0 self.set_allow_infeasible_solutions() if isinstance(cost, (BaseLikelihood, DesignCost)): self.minimising = False + + # Set default bounds (for all or no parameters) + self.bounds = cost.parameters.get_bounds() + + # Set default initial standard deviation (for all or no parameters) + self.sigma0 = cost.parameters.get_sigma0() or self.sigma0 + else: try: cost_test = cost(optimiser_kwargs.get("x0", [])) diff --git a/pybop/optimisers/base_pints_optimiser.py b/pybop/optimisers/base_pints_optimiser.py index 543d32b6..2c9d17f0 100644 --- a/pybop/optimisers/base_pints_optimiser.py +++ b/pybop/optimisers/base_pints_optimiser.py @@ -31,8 +31,9 @@ def __init__(self, cost, pints_optimiser, **optimiser_kwargs): self._callback = None self._max_iterations = None self._min_iterations = 2 - self._unchanged_threshold = 1e-5 # smallest significant f change self._unchanged_max_iterations = 15 + self._absolute_tolerance = 1e-5 + self._relative_tolerance = 1e-2 self._max_evaluations = None self._threshold = None self._evaluations = None @@ -79,17 +80,20 @@ def _set_up_optimiser(self): elif key == "min_iterations": self.set_min_iterations(self.unset_options.pop(key)) elif key == "max_unchanged_iterations": - if "threshold" in self.unset_options.keys(): - self.set_max_unchanged_iterations( - self.unset_options.pop(key), - self.unset_options.pop("threshold"), + max_unchanged_kwargs = {"iterations": self.unset_options.pop(key)} + if "absolute_tolerance" in self.unset_options.keys(): + max_unchanged_kwargs["absolute_tolerance"] = self.unset_options.pop( + "absolute_tolerance" ) - else: - self.set_max_unchanged_iterations(self.unset_options.pop(key)) - elif key == "threshold": - pass # only used with unchanged_max_iterations + if "relative_tolerance" in self.unset_options.keys(): + max_unchanged_kwargs["relative_tolerance"] = self.unset_options.pop( + "relative_tolerance" + ) + self.set_max_unchanged_iterations(**max_unchanged_kwargs) elif key == "max_evaluations": self.set_max_evaluations(self.unset_options.pop(key)) + elif key == "threshold": + self.set_threshold(self.unset_options.pop(key)) def _sanitise_inputs(self): """ @@ -108,12 +112,8 @@ def _sanitise_inputs(self): self.unset_options.pop("options") # Check for duplicate keywords - expected_keys = [ - "max_iterations", - "popsize", - "threshold", - ] - alternative_keys = ["maxiter", "population_size", "tol"] + expected_keys = ["max_iterations", "popsize"] + alternative_keys = ["maxiter", "population_size"] for exp_key, alt_key in zip(expected_keys, alternative_keys): if alt_key in self.unset_options.keys(): if exp_key in self.unset_options.keys(): @@ -238,9 +238,11 @@ def f(x, grad=None): fg = self.pints_optimiser.f_guessed() fg_user = (fb, fg) if self.minimising else (-fb, -fg) - # Check for significant changes + # Check for significant changes against the absolute and relative tolerance f_new = fg if self._use_f_guessed else fb - if np.abs(f_new - f_sig) >= self._unchanged_threshold: + if np.abs(f_new - f_sig) >= np.maximum( + self._absolute_tolerance, self._relative_tolerance * np.abs(f_sig) + ): unchanged_iterations = 0 f_sig = f_new else: @@ -429,7 +431,9 @@ def set_min_iterations(self, iterations=2): raise ValueError("Minimum number of iterations cannot be negative.") self._min_iterations = iterations - def set_max_unchanged_iterations(self, iterations=15, threshold=1e-5): + def set_max_unchanged_iterations( + self, iterations=15, absolute_tolerance=1e-5, relative_tolerance=1e-2 + ): """ Set the maximum number of iterations without significant change as a stopping criterion. Credit: PINTS @@ -439,21 +443,29 @@ def set_max_unchanged_iterations(self, iterations=15, threshold=1e-5): iterations : int, optional The maximum number of unchanged iterations to run (default: 15). Set to `None` to remove this stopping criterion. - threshold : float, optional - The minimum significant change in the objective function value that resets the - unchanged iteration counter (default: 1e-5). + absolute_tolerance : float, optional + The minimum significant change (absolute tolerance) in the objective function value + that resets the unchanged iteration counter (default: 1e-5). + relative_tolerance : float, optional + The minimum significant proportional change (relative tolerance) in the objective function + value that resets the unchanged iteration counter (default: 1e-2). """ if iterations is not None: iterations = int(iterations) if iterations < 0: raise ValueError("Maximum number of iterations cannot be negative.") - threshold = float(threshold) - if threshold < 0: - raise ValueError("Minimum significant change cannot be negative.") + absolute_tolerance = float(absolute_tolerance) + if absolute_tolerance < 0: + raise ValueError("Absolute tolerance cannot be negative.") + + relative_tolerance = float(relative_tolerance) + if relative_tolerance < 0: + raise ValueError("Relative tolerance cannot be negative.") self._unchanged_max_iterations = iterations - self._unchanged_threshold = threshold + self._absolute_tolerance = absolute_tolerance + self._relative_tolerance = relative_tolerance def set_max_evaluations(self, evaluations=None): """ @@ -472,6 +484,27 @@ def set_max_evaluations(self, evaluations=None): raise ValueError("Maximum number of evaluations cannot be negative.") self._max_evaluations = evaluations + def set_threshold(self, threshold=None): + """ + Adds a stopping criterion, allowing the routine to halt once the + objective function goes below a set ``threshold``. + + This criterion is disabled by default, but can be enabled by calling + this method with a valid ``threshold``. To disable it, use + ``set_threshold(None)``. + Credit: PINTS + + Parameters + ---------- + threshold : float, optional + The threshold below which the objective function value is considered optimal + (default: None). + """ + if threshold is None: + self._threshold = None + else: + self._threshold = float(threshold) + class Result: """ diff --git a/pybop/optimisers/scipy_optimisers.py b/pybop/optimisers/scipy_optimisers.py index 6c249809..4d9a0e90 100644 --- a/pybop/optimisers/scipy_optimisers.py +++ b/pybop/optimisers/scipy_optimisers.py @@ -39,12 +39,8 @@ def _sanitise_inputs(self): self.unset_options.pop("options") # Check for duplicate keywords - expected_keys = ["maxiter", "popsize", "tol"] - alternative_keys = [ - "max_iterations", - "population_size", - "threshold", - ] + expected_keys = ["maxiter", "popsize"] + alternative_keys = ["max_iterations", "population_size"] for exp_key, alt_key in zip(expected_keys, alternative_keys): if alt_key in self.unset_options.keys(): if exp_key in self.unset_options.keys(): @@ -165,7 +161,7 @@ def callback(x): self._cost0 = np.abs(self.cost(self.x0)) if np.isinf(self._cost0): for i in range(1, self.num_resamples): - x0 = self.cost.problem.sample_initial_conditions(seed=i) + x0 = self.cost.parameters.rvs(1) self._cost0 = np.abs(self.cost(x0)) if not np.isinf(self._cost0): break @@ -259,9 +255,10 @@ def _set_up_optimiser(self): ): raise ValueError("Bounds must be specified for differential_evolution.") - # Apply default maxiter + # Apply default maxiter and tolerance self._options = dict() self._options["maxiter"] = self.default_max_iterations + self._options["tol"] = 1e-5 # Apply additional options and remove them from the options dictionary key_list = list(self.unset_options.keys()) diff --git a/pybop/parameters/parameter.py b/pybop/parameters/parameter.py index 52b700bb..d8117f8f 100644 --- a/pybop/parameters/parameter.py +++ b/pybop/parameters/parameter.py @@ -1,3 +1,6 @@ +from collections import OrderedDict +from typing import Dict, List + import numpy as np @@ -82,6 +85,7 @@ def update(self, value=None, initial_value=None): if value is not None: self.value = value elif initial_value is not None: + self.initial_value = initial_value self.value = initial_value else: raise ValueError("No value provided to update parameter") @@ -143,3 +147,240 @@ def set_bounds(self, bounds=None): self.upper_bound = bounds[1] self.bounds = bounds + + +class Parameters: + """ + Represents a set of uncertain parameters within the PyBOP framework. + + This class encapsulates the definition of a parameter, including its name, prior + distribution, initial value, bounds, and a margin to ensure the parameter stays + within feasible limits during optimisation or sampling. + + Parameters + ---------- + parameter_list : pybop.Parameter or Dict + """ + + def __init__(self, *args): + self.param = OrderedDict() + for param in args: + self.add(param) + + def __getitem__(self, key: str): + """ + Return the parameter dictionary corresponding to a particular key. + + Parameters + ---------- + key : str + The name of a parameter. + + Returns + ------- + pybop.Parameter + The Parameter object. + + Raises + ------ + ValueError + The key must be the name of one of the parameters. + """ + if key not in self.param.keys(): + raise ValueError(f"The key {key} is not the name of a parameter.") + + return self.param[key] + + def __len__(self) -> int: + return len(self.param) + + def keys(self) -> List: + """ + A list of parameter names + """ + return list(self.param.keys()) + + def __iter__(self): + self.index = 0 + return self + + def __next__(self): + parameter_names = self.keys() + if self.index == len(parameter_names): + raise StopIteration + name = parameter_names[self.index] + self.index = self.index + 1 + return self.param[name] + + def add(self, parameter): + """ + Construct the parameter class with a name, initial value, prior, and bounds. + """ + if isinstance(parameter, Parameter): + if parameter.name in self.param.keys(): + raise ValueError( + f"There is already a parameter with the name {parameter.name} " + + "in the Parameters object. Please remove the duplicate entry." + ) + self.param[parameter.name] = parameter + elif isinstance(parameter, dict): + if "name" not in parameter.keys(): + raise Exception("Parameter requires a name.") + name = parameter["name"] + if name in self.param.keys(): + raise ValueError( + f"There is already a parameter with the name {name} " + + "in the Parameters object. Please remove the duplicate entry." + ) + self.param[name] = Parameter(**parameter) + else: + raise TypeError("Each parameter input must be a Parameter or a dictionary.") + + def remove(self, parameter_name): + """ + Remove the `Parameter` object from the `Parameters` dictionary. + """ + if not isinstance(parameter_name, str): + raise TypeError("The input parameter_name is not a string.") + if parameter_name not in self.param.keys(): + raise ValueError("This parameter does not exist in the Parameters object.") + + # Remove the parameter + self.param.pop(parameter_name) + + def get_bounds(self) -> Dict: + """ + Get bounds, for either all or no parameters. + """ + all_unbounded = True # assumption + bounds = {"lower": [], "upper": []} + + for param in self.param.values(): + if param.bounds is not None: + bounds["lower"].append(param.bounds[0]) + bounds["upper"].append(param.bounds[1]) + all_unbounded = False + else: + bounds["lower"].append(-np.inf) + bounds["upper"].append(np.inf) + if all_unbounded: + bounds = None + + return bounds + + def update(self, values): + """ + Set value of each parameter. + """ + for i, param in enumerate(self.param.values()): + param.update(value=values[i]) + + def rvs(self, n_samples: int) -> List: + """ + Draw random samples from each parameter's prior distribution. + + The samples are constrained to be within the parameter's bounds, excluding + a predefined margin at the boundaries. + + Parameters + ---------- + n_samples : int + The number of samples to draw. + + Returns + ------- + array-like + An array of samples drawn from the prior distribution within each parameter's bounds. + """ + all_samples = [] + + for param in self.param.values(): + samples = param.rvs(n_samples) + + # Constrain samples to be within bounds + if param.bounds is not None: + offset = param.margin * (param.upper_bound - param.lower_bound) + samples = np.clip( + samples, param.lower_bound + offset, param.upper_bound - offset + ) + + all_samples.append(samples) + + return all_samples + + def get_sigma0(self) -> List: + """ + Get the standard deviation, for either all or no parameters. + """ + all_have_sigma = True # assumption + sigma0 = [] + + for param in self.param.values(): + if hasattr(param.prior, "sigma"): + sigma0.append(param.prior.sigma) + else: + all_have_sigma = False + if not all_have_sigma: + sigma0 = None + + return sigma0 + + def initial_value(self) -> List: + """ + Return the initial value of each parameter. + """ + initial_values = [] + + for param in self.param.values(): + if param.initial_value is None: + initial_value = param.rvs(1) + param.update(initial_value=initial_value[0]) + initial_values.append(param.initial_value) + + return initial_values + + def current_value(self) -> List: + """ + Return the current value of each parameter. + """ + current_values = [] + + for param in self.param.values(): + current_values.append(param.value) + + return current_values + + def true_value(self) -> List: + """ + Return the true value of each parameter. + """ + true_values = [] + + for param in self.param.values(): + true_values.append(param.true_value) + + return true_values + + def get_bounds_for_plotly(self): + """ + Retrieve parameter bounds in the format expected by Plotly. + + Returns + ------- + bounds : numpy.ndarray + An array of shape (n_parameters, 2) containing the bounds for each parameter. + """ + bounds = np.empty((len(self), 2)) + + for i, param in enumerate(self.param.values()): + if param.bounds is not None: + bounds[i] = param.bounds + else: + raise ValueError("All parameters require bounds for plotting.") + + return bounds + + def as_dict(self, values=None) -> Dict: + if values is None: + values = self.current_value() + return {key: values[i] for i, key in enumerate(self.param.keys())} diff --git a/pybop/plotting/plot2d.py b/pybop/plotting/plot2d.py index 95727961..eb7cc387 100644 --- a/pybop/plotting/plot2d.py +++ b/pybop/plotting/plot2d.py @@ -23,7 +23,8 @@ def plot2d( gradient : bool, optional If True, the gradient is shown (default: False). bounds : numpy.ndarray, optional - A 2x2 array specifying the [min, max] bounds for each parameter. If None, uses `get_param_bounds`. + A 2x2 array specifying the [min, max] bounds for each parameter. If None, uses + `cost.parameters.get_bounds_for_plotly`. steps : int, optional The number of intervals to divide the parameter space into along each dimension (default is 10). show : bool, optional @@ -55,9 +56,7 @@ def plot2d( # Set up parameter bounds if bounds is None: - bounds = get_param_bounds(cost) - else: - bounds = bounds + bounds = cost.parameters.get_bounds_for_plotly() # Generate grid x = np.linspace(bounds[0, 0], bounds[0, 1], steps) @@ -102,8 +101,9 @@ def plot2d( yaxis=dict(range=bounds[1]), ) if hasattr(cost, "parameters"): - layout_options["xaxis_title"] = cost.parameters[0].name - layout_options["yaxis_title"] = cost.parameters[1].name + name = cost.parameters.keys() + layout_options["xaxis_title"] = name[0] + layout_options["yaxis_title"] = name[1] layout = go.Layout(layout_options) # Create contour plot and update the layout @@ -180,27 +180,3 @@ def plot2d( return fig, grad_figs return fig - - -def get_param_bounds(cost): - """ - Retrieve parameter bounds from a cost object's associated problem parameters. - - Parameters - ---------- - cost : object - The cost object with an associated 'problem' attribute containing 'parameters'. - - Returns - ------- - bounds : numpy.ndarray - An array of shape (n_parameters, 2) containing the bounds for each parameter. - """ - bounds = np.empty((len(cost.parameters), 2)) - for i, param in enumerate(cost.parameters): - if param.bounds is not None: - bounds[i] = param.bounds - else: - raise ValueError("plot2d could not find bounds required for plotting") - - return bounds diff --git a/pybop/plotting/plot_parameters.py b/pybop/plotting/plot_parameters.py index cbc1718f..1784569d 100644 --- a/pybop/plotting/plot_parameters.py +++ b/pybop/plotting/plot_parameters.py @@ -47,10 +47,9 @@ def plot_parameters(optim, show=True, **layout_kwargs): # Create lists of axis titles and trace names axis_titles = [] - trace_names = [] - for param in parameters: - axis_titles.append(("Function Call", param.name)) - trace_names.append(param.name) + trace_names = parameters.keys() + for name in trace_names: + axis_titles.append(("Function Call", name)) # Set subplot layout options layout_options = dict( diff --git a/pybop/problems/base_problem.py b/pybop/problems/base_problem.py index 530bad10..0e2d4ac7 100644 --- a/pybop/problems/base_problem.py +++ b/pybop/problems/base_problem.py @@ -1,5 +1,3 @@ -import numpy as np - import pybop @@ -9,8 +7,8 @@ class BaseProblem: Parameters ---------- - parameters : list - List of parameters for the problem. + parameters : pybop.Parameter or pybop.Parameters + An object or list of the parameters for the problem. model : object, optional The model to be used for the problem (default: None). check_model : bool, optional @@ -21,8 +19,6 @@ class BaseProblem: Additional variables to observe and store in the solution (default: []). init_soc : float, optional Initial state of charge (default: None). - x0 : np.ndarray, optional - Initial parameter values (default: None). """ def __init__( @@ -33,8 +29,24 @@ def __init__( signal=["Voltage [V]"], additional_variables=[], init_soc=None, - x0=None, ): + # Check if parameters is a list of pybop.Parameter objects + if isinstance(parameters, list): + if all(isinstance(param, pybop.Parameter) for param in parameters): + parameters = pybop.Parameters(*parameters) + else: + raise TypeError( + "All elements in the list must be pybop.Parameter objects." + ) + # Check if parameters is a single pybop.Parameter object + elif isinstance(parameters, pybop.Parameter): + parameters = pybop.Parameters(parameters) + # Check if parameters is already a pybop.Parameters object + elif not isinstance(parameters, pybop.Parameters): + raise TypeError( + "The input parameters must be a pybop Parameter, a list of pybop.Parameter objects, or a pybop Parameters object." + ) + self.parameters = parameters self._model = model self.check_model = check_model @@ -44,8 +56,6 @@ def __init__( raise ValueError("Signal should be either a string or list of strings.") self.signal = signal self.init_soc = init_soc - self.x0 = x0 - self.n_parameters = len(self.parameters) self.n_outputs = len(self.signal) self._time_data = None self._target = None @@ -55,71 +65,12 @@ def __init__( else: self.additional_variables = [] - # Set bounds (for all or no parameters) - all_unbounded = True # assumption - self.bounds = {"lower": [], "upper": []} - for param in self.parameters: - if param.bounds is not None: - self.bounds["lower"].append(param.bounds[0]) - self.bounds["upper"].append(param.bounds[1]) - all_unbounded = False - else: - self.bounds["lower"].append(-np.inf) - self.bounds["upper"].append(np.inf) - if all_unbounded: - self.bounds = None - - # Set initial standard deviation (for all or no parameters) - all_have_sigma = True # assumption - self.sigma0 = [] - for param in self.parameters: - if hasattr(param.prior, "sigma"): - self.sigma0.append(param.prior.sigma) - else: - all_have_sigma = False - if not all_have_sigma: - self.sigma0 = None - - # Set initial conditions - if self.x0 is None: - self.x0 = self.sample_initial_conditions() - elif len(self.x0) != self.n_parameters: - raise ValueError("x0 dimensions do not match number of parameters") - - # Add the initial values to the parameter definitions - for i, param in enumerate(self.parameters): - param.update(initial_value=self.x0[i]) - - def sample_initial_conditions( - self, num_samples: int = 1, seed: int = None - ) -> np.ndarray: - """ - Sample initial conditions for the problem. - - Parameters - ---------- - num_samples : int, optional - The number of initial condition samples to generate (default is 1). - seed : int, optional - The seed value for the random number generator (default is None). + # Set initial values + self.x0 = self.parameters.initial_value() - Returns - ------- - np.ndarray - An array of initial conditions for the problem. Each row corresponds to a sample. - """ - # Create a local random generator - rng = np.random.default_rng(seed) - - # Initialize - x0 = np.zeros((num_samples, self.n_parameters)) - - # Vectorized sampling from the parameters - for j, param in enumerate(self.parameters): - x0[:, j] = param.rvs(num_samples, random_state=rng) - - # Return a flattened array if only one sample, otherwise return the 2D array - return x0[0] if num_samples == 1 else x0 + @property + def n_parameters(self): + return len(self.parameters) def evaluate(self, x): """ diff --git a/pybop/problems/design_problem.py b/pybop/problems/design_problem.py index 29f58ba9..3217ca95 100644 --- a/pybop/problems/design_problem.py +++ b/pybop/problems/design_problem.py @@ -13,8 +13,8 @@ class DesignProblem(BaseProblem): ---------- model : object The model to apply the design to. - parameters : list - List of parameters for the problem. + parameters : pybop.Parameter or pybop.Parameters + An object or list of the parameters for the problem. experiment : object The experimental setup to apply the model to. check_model : bool, optional @@ -25,8 +25,6 @@ class DesignProblem(BaseProblem): Additional variables to observe and store in the solution (default additions are: ["Time [s]", "Current [A]"]). init_soc : float, optional Initial state of charge (default: None). - x0 : np.ndarray, optional - Initial parameter values (default: None). """ def __init__( @@ -38,14 +36,18 @@ def __init__( signal=["Voltage [V]"], additional_variables=[], init_soc=None, - x0=None, ): # Add time and current and remove duplicates additional_variables.extend(["Time [s]", "Current [A]"]) additional_variables = list(set(additional_variables)) super().__init__( - parameters, model, check_model, signal, additional_variables, init_soc, x0 + parameters, + model, + check_model, + signal, + additional_variables, + init_soc, ) self.experiment = experiment @@ -53,8 +55,6 @@ def __init__( if experiment is not None: # Leave the build until later to apply the experiment self._model.parameters = self.parameters - if self.parameters is not None: - self._model.fit_keys = [param.name for param in self.parameters] elif self._model._built_model is None: self._model.build( @@ -84,7 +84,6 @@ def evaluate(self, x): y : np.ndarray The model output y(t) simulated with inputs x. """ - sol = self._model.predict( inputs=x, experiment=self.experiment, diff --git a/pybop/problems/fitting_problem.py b/pybop/problems/fitting_problem.py index 6bfa3814..75f32bac 100644 --- a/pybop/problems/fitting_problem.py +++ b/pybop/problems/fitting_problem.py @@ -13,8 +13,8 @@ class FittingProblem(BaseProblem): ---------- model : object The model to fit. - parameters : list - List of parameters for the problem. + parameters : pybop.Parameter or pybop.Parameters + An object or list of the parameters for the problem. dataset : Dataset Dataset object containing the data to fit the model to. signal : str, optional @@ -23,8 +23,6 @@ class FittingProblem(BaseProblem): Additional variables to observe and store in the solution (default additions are: ["Time [s]"]). init_soc : float, optional Initial state of charge (default: None). - x0 : np.ndarray, optional - Initial parameter values (default: None). """ def __init__( @@ -36,14 +34,13 @@ def __init__( signal=["Voltage [V]"], additional_variables=[], init_soc=None, - x0=None, ): # Add time and remove duplicates additional_variables.extend(["Time [s]"]) additional_variables = list(set(additional_variables)) super().__init__( - parameters, model, check_model, signal, additional_variables, init_soc, x0 + parameters, model, check_model, signal, additional_variables, init_soc ) self._dataset = dataset.data self.x = self.x0 @@ -60,7 +57,6 @@ def __init__( if model is not None: self._model.signal = self.signal self._model.additional_variables = self.additional_variables - self._model.n_parameters = self.n_parameters self._model.n_outputs = self.n_outputs self._model.n_time_data = self.n_time_data @@ -93,9 +89,7 @@ def evaluate(self, x): The model output y(t) simulated with inputs x. """ if np.any(x != self.x) and self._model.matched_parameters: - for i, param in enumerate(self.parameters): - param.update(value=x[i]) - + self.parameters.update(values=x) self._model.rebuild(parameters=self.parameters) self.x = x diff --git a/pyproject.toml b/pyproject.toml index bf6080e0..14ec41be 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -58,7 +58,6 @@ dev = [ "pytest-cov", "pytest-mock", "pytest-xdist", - "flaky", "ruff", ] scifem = [ diff --git a/tests/integration/test_model_experiment_changes.py b/tests/integration/test_model_experiment_changes.py index b3c82277..6902f873 100644 --- a/tests/integration/test_model_experiment_changes.py +++ b/tests/integration/test_model_experiment_changes.py @@ -11,29 +11,34 @@ class TestModelAndExperimentChanges: @pytest.fixture( params=[ - pybop.Parameter( # geometric parameter - "Negative particle radius [m]", - prior=pybop.Gaussian(6e-06, 0.1e-6), - bounds=[1e-6, 9e-6], - true_value=5.86e-6, + pybop.Parameters( + pybop.Parameter( # geometric parameter + "Negative particle radius [m]", + prior=pybop.Gaussian(6e-06, 0.1e-6), + bounds=[1e-6, 9e-6], + true_value=5.86e-6, + initial_value=5.86e-6, + ), ), - pybop.Parameter( # non-geometric parameter - "Positive electrode diffusivity [m2.s-1]", - prior=pybop.Gaussian(3.43e-15, 1e-15), - bounds=[1e-15, 5e-15], - true_value=4e-15, + pybop.Parameters( + pybop.Parameter( # non-geometric parameter + "Positive electrode diffusivity [m2.s-1]", + prior=pybop.Gaussian(3.43e-15, 1e-15), + bounds=[1e-15, 5e-15], + true_value=4e-15, + initial_value=4e-15, + ), ), ] ) - def parameter(self, request): + def parameters(self, request): return request.param @pytest.mark.integration - def test_changing_experiment(self, parameter): + def test_changing_experiment(self, parameters): # Change the experiment and check that the results are different. parameter_set = pybop.ParameterSet.pybamm("Chen2020") - parameters = [parameter] init_soc = 0.5 model = pybop.lithium_ion.SPM(parameter_set=parameter_set) @@ -43,7 +48,7 @@ def test_changing_experiment(self, parameter): experiment = pybop.Experiment(["Charge at 1C until 4.1 V (2 seconds period)"]) solution_2 = model.predict( - init_soc=init_soc, experiment=experiment, inputs=[parameter.true_value] + init_soc=init_soc, experiment=experiment, inputs=parameters.true_value() ) cost_2 = self.final_cost(solution_2, model, parameters, init_soc) @@ -58,11 +63,10 @@ def test_changing_experiment(self, parameter): np.testing.assert_allclose(cost_2, 0, atol=1e-5) @pytest.mark.integration - def test_changing_model(self, parameter): + def test_changing_model(self, parameters): # Change the model and check that the results are different. parameter_set = pybop.ParameterSet.pybamm("Chen2020") - parameters = [parameter] init_soc = 0.5 experiment = pybop.Experiment(["Charge at 1C until 4.1 V (2 seconds period)"]) @@ -86,7 +90,6 @@ def test_changing_model(self, parameter): def final_cost(self, solution, model, parameters, init_soc): # Compute the cost corresponding to a particular solution - x0 = np.array([parameters[0].true_value]) dataset = pybop.Dataset( { "Time [s]": solution["Time [s]"].data, @@ -96,7 +99,7 @@ def final_cost(self, solution, model, parameters, init_soc): ) signal = ["Voltage [V]"] problem = pybop.FittingProblem( - model, parameters, dataset, signal=signal, x0=x0, init_soc=init_soc + model, parameters, dataset, signal=signal, init_soc=init_soc ) cost = pybop.RootMeanSquaredError(problem) optim = pybop.PSO(cost) diff --git a/tests/integration/test_optimisation_options.py b/tests/integration/test_optimisation_options.py index 1505a37d..dcd94276 100644 --- a/tests/integration/test_optimisation_options.py +++ b/tests/integration/test_optimisation_options.py @@ -24,7 +24,7 @@ def model(self): @pytest.fixture def parameters(self): - return [ + return pybop.Parameters( pybop.Parameter( "Negative electrode active material volume fraction", prior=pybop.Gaussian(0.55, 0.05), @@ -35,7 +35,7 @@ def parameters(self): prior=pybop.Gaussian(0.55, 0.05), # no bounds ), - ] + ) @pytest.fixture( params=[ @@ -54,7 +54,7 @@ def noise(self, sigma, values): def spm_costs(self, model, parameters, cost_class): # Form dataset init_soc = 0.5 - solution = self.getdata(model, self.ground_truth, init_soc) + solution = self.get_data(model, parameters, self.ground_truth, init_soc) dataset = pybop.Dataset( { "Time [s]": solution["Time [s]"].data, @@ -85,9 +85,9 @@ def test_optimisation_f_guessed(self, f_guessed, spm_costs): optim = pybop.XNES( cost=spm_costs, sigma0=0.05, - max_iterations=125, + max_iterations=250, max_unchanged_iterations=35, - threshold=1e-5, + absolute_tolerance=1e-5, use_f_guessed=f_guessed, ) @@ -104,15 +104,10 @@ def test_optimisation_f_guessed(self, f_guessed, spm_costs): assert initial_cost > final_cost else: assert initial_cost < final_cost - np.testing.assert_allclose(x, self.ground_truth, atol=2.5e-2) + np.testing.assert_allclose(x, self.ground_truth, atol=1.5e-2) - def getdata(self, model, x, init_soc): - model.parameter_set.update( - { - "Negative electrode active material volume fraction": x[0], - "Positive electrode active material volume fraction": x[1], - } - ) + def get_data(self, model, parameters, x, init_soc): + model.parameters = parameters experiment = pybop.Experiment( [ ( @@ -122,5 +117,5 @@ def getdata(self, model, x, init_soc): ] * 2 ) - sim = model.predict(init_soc=init_soc, experiment=experiment) + sim = model.predict(init_soc=init_soc, experiment=experiment, inputs=x) return sim diff --git a/tests/integration/test_spm_parameterisations.py b/tests/integration/test_spm_parameterisations.py index 470bfe0d..6e5108d7 100644 --- a/tests/integration/test_spm_parameterisations.py +++ b/tests/integration/test_spm_parameterisations.py @@ -1,7 +1,5 @@ import numpy as np import pytest -from flaky import flaky -from pybamm import __version__ as pybamm_version import pybop @@ -24,7 +22,7 @@ def model(self): @pytest.fixture def parameters(self): - return [ + return pybop.Parameters( pybop.Parameter( "Negative electrode active material volume fraction", prior=pybop.Uniform(0.4, 0.7), @@ -35,7 +33,7 @@ def parameters(self): prior=pybop.Uniform(0.4, 0.7), # no bounds ), - ] + ) @pytest.fixture(params=[0.4, 0.7]) def init_soc(self, request): @@ -58,7 +56,7 @@ def noise(self, sigma, values): @pytest.fixture def spm_costs(self, model, parameters, cost_class, init_soc): # Form dataset - solution = self.getdata(model, self.ground_truth, init_soc) + solution = self.get_data(model, parameters, self.ground_truth, init_soc) dataset = pybop.Dataset( { "Time [s]": solution["Time [s]"].data, @@ -91,7 +89,6 @@ def spm_costs(self, model, parameters, cost_class, init_soc): pybop.XNES, ], ) - @flaky(max_runs=3, min_passes=1) @pytest.mark.integration def test_spm_optimisers(self, optimiser, spm_costs): x0 = spm_costs.x0 @@ -99,25 +96,22 @@ def test_spm_optimisers(self, optimiser, spm_costs): if optimiser in [ pybop.SciPyDifferentialEvolution, ]: - spm_costs.problem.parameters[1].set_bounds( - [0.375, 0.725] - ) # Large range to ensure IC within bounds - bounds = {"lower": [], "upper": []} - for param in spm_costs.problem.parameters: - bounds["lower"].append(param.bounds[0]) - bounds["upper"].append(param.bounds[1]) + spm_costs.problem.parameters[ + "Positive electrode active material volume fraction" + ].set_bounds([0.375, 0.725]) # Large range to ensure IC within bounds + bounds = spm_costs.problem.parameters.get_bounds() spm_costs.problem.bounds = bounds spm_costs.bounds = bounds # Test each optimiser if optimiser in [pybop.PSO]: optim = pybop.Optimisation( - cost=spm_costs, optimiser=optimiser, sigma0=0.05, max_iterations=125 + cost=spm_costs, optimiser=optimiser, sigma0=0.05, max_iterations=250 ) else: - optim = optimiser(cost=spm_costs, sigma0=0.05, max_iterations=125) + optim = optimiser(cost=spm_costs, sigma0=0.05, max_iterations=250) if issubclass(optimiser, pybop.BasePintsOptimiser): - optim.set_max_unchanged_iterations(iterations=35, threshold=1e-5) + optim.set_max_unchanged_iterations(iterations=35, absolute_tolerance=1e-5) initial_cost = optim.cost(x0) x, final_cost = optim.run() @@ -128,16 +122,13 @@ def test_spm_optimisers(self, optimiser, spm_costs): assert initial_cost > final_cost else: assert initial_cost < final_cost - if pybamm_version <= "23.9": - np.testing.assert_allclose(x, self.ground_truth, atol=2.5e-2) - else: - np.testing.assert_allclose(x, self.ground_truth, atol=1.75e-2) + np.testing.assert_allclose(x, self.ground_truth, atol=1.5e-2) @pytest.fixture def spm_two_signal_cost(self, parameters, model, cost_class): # Form dataset init_soc = 0.5 - solution = self.getdata(model, self.ground_truth, init_soc) + solution = self.get_data(model, parameters, self.ground_truth, init_soc) dataset = pybop.Dataset( { "Time [s]": solution["Time [s]"].data, @@ -177,22 +168,21 @@ def test_multiple_signals(self, multi_optimiser, spm_two_signal_cost): x0 = spm_two_signal_cost.x0 # Some optimisers require a complete set of bounds if multi_optimiser in [pybop.SciPyDifferentialEvolution]: - spm_two_signal_cost.problem.parameters[1].set_bounds( - [0.375, 0.725] - ) # Large range to ensure IC within bounds - bounds = {"lower": [], "upper": []} - for param in spm_two_signal_cost.problem.parameters: - bounds["lower"].append(param.bounds[0]) - bounds["upper"].append(param.bounds[1]) + spm_two_signal_cost.problem.parameters[ + "Positive electrode active material volume fraction" + ].set_bounds([0.375, 0.725]) # Large range to ensure IC within bounds + bounds = spm_two_signal_cost.problem.parameters.get_bounds() spm_two_signal_cost.problem.bounds = bounds spm_two_signal_cost.bounds = bounds # Test each optimiser optim = multi_optimiser( - cost=spm_two_signal_cost, sigma0=0.03, max_iterations=125 + cost=spm_two_signal_cost, + sigma0=0.03, + max_iterations=250, ) if issubclass(multi_optimiser, pybop.BasePintsOptimiser): - optim.set_max_unchanged_iterations(iterations=35, threshold=5e-4) + optim.set_max_unchanged_iterations(iterations=35, absolute_tolerance=1e-5) initial_cost = optim.cost(spm_two_signal_cost.x0) x, final_cost = optim.run() @@ -203,7 +193,7 @@ def test_multiple_signals(self, multi_optimiser, spm_two_signal_cost): assert initial_cost > final_cost else: assert initial_cost < final_cost - np.testing.assert_allclose(x, self.ground_truth, atol=2.5e-2) + np.testing.assert_allclose(x, self.ground_truth, atol=1.5e-2) @pytest.mark.parametrize("init_soc", [0.4, 0.6]) @pytest.mark.integration @@ -214,7 +204,7 @@ def test_model_misparameterisation(self, parameters, model, init_soc): second_model = pybop.lithium_ion.SPMe(parameter_set=second_parameter_set) # Form dataset - solution = self.getdata(second_model, self.ground_truth, init_soc) + solution = self.get_data(second_model, parameters, self.ground_truth, init_soc) dataset = pybop.Dataset( { "Time [s]": solution["Time [s]"].data, @@ -244,13 +234,8 @@ def test_model_misparameterisation(self, parameters, model, init_soc): with np.testing.assert_raises(AssertionError): np.testing.assert_allclose(x, self.ground_truth, atol=2e-2) - def getdata(self, model, x, init_soc): - model.parameter_set.update( - { - "Negative electrode active material volume fraction": x[0], - "Positive electrode active material volume fraction": x[1], - } - ) + def get_data(self, model, parameters, x, init_soc): + model.parameters = parameters experiment = pybop.Experiment( [ ( @@ -260,5 +245,5 @@ def getdata(self, model, x, init_soc): ] * 2 ) - sim = model.predict(init_soc=init_soc, experiment=experiment) + sim = model.predict(init_soc=init_soc, experiment=experiment, inputs=x) return sim diff --git a/tests/integration/test_thevenin_parameterisation.py b/tests/integration/test_thevenin_parameterisation.py index 5ac6e84e..57bb0689 100644 --- a/tests/integration/test_thevenin_parameterisation.py +++ b/tests/integration/test_thevenin_parameterisation.py @@ -1,19 +1,16 @@ import numpy as np import pytest -from flaky import flaky import pybop class TestTheveninParameterisation: """ - A class to test the flaky optimisers on a simple model. + A class to test a subset of optimisers on a simple model. """ @pytest.fixture(autouse=True) def setup(self): - # Set random seed for reproducibility - np.random.seed(0) self.ground_truth = np.array([0.05, 0.05]) + np.random.normal( loc=0.0, scale=0.01, size=2 ) @@ -24,22 +21,23 @@ def model(self): json_path="examples/scripts/parameters/initial_ecm_parameters.json" ) parameter_set.import_parameters() + parameter_set.params.update({"C1 [F]": 1000}) return pybop.empirical.Thevenin(parameter_set=parameter_set) @pytest.fixture def parameters(self): - return [ + return pybop.Parameters( pybop.Parameter( "R0 [Ohm]", prior=pybop.Gaussian(0.05, 0.01), - bounds=[0.01, 0.1], + bounds=[0, 0.1], ), pybop.Parameter( "R1 [Ohm]", prior=pybop.Gaussian(0.05, 0.01), - bounds=[0.01, 0.1], + bounds=[0, 0.1], ), - ] + ) @pytest.fixture(params=[pybop.RootMeanSquaredError, pybop.SumSquaredError]) def cost_class(self, request): @@ -48,7 +46,7 @@ def cost_class(self, request): @pytest.fixture def cost(self, model, parameters, cost_class): # Form dataset - solution = self.getdata(model, self.ground_truth) + solution = self.get_data(model, parameters, self.ground_truth) dataset = pybop.Dataset( { "Time [s]": solution["Time [s]"].data, @@ -65,7 +63,6 @@ def cost(self, model, parameters, cost_class): "optimiser", [pybop.SciPyMinimize, pybop.GradientDescent, pybop.PSO], ) - @flaky(max_runs=8, min_passes=1) # These can be very flaky @pytest.mark.integration def test_optimisers_on_simple_model(self, optimiser, cost): x0 = cost.x0 @@ -82,7 +79,7 @@ def test_optimisers_on_simple_model(self, optimiser, cost): max_iterations=250, ) if isinstance(optimiser, pybop.BasePintsOptimiser): - optim.set_max_unchanged_iterations(iterations=55, threshold=1e-5) + optim.set_max_unchanged_iterations(iterations=35, absolute_tolerance=1e-5) initial_cost = optim.cost(x0) x, final_cost = optim.run() @@ -93,19 +90,17 @@ def test_optimisers_on_simple_model(self, optimiser, cost): assert initial_cost > final_cost else: assert initial_cost < final_cost - np.testing.assert_allclose(x, self.ground_truth, atol=1e-2) + np.testing.assert_allclose(x, self.ground_truth, atol=1.5e-2) - def getdata(self, model, x): - model.parameter_set.update( - { - "R0 [Ohm]": x[0], - "R1 [Ohm]": x[1], - } - ) + def get_data(self, model, parameters, x): + model.parameters = parameters experiment = pybop.Experiment( [ - ("Discharge at 0.5C for 2 minutes (1 second period)",), + ( + "Discharge at 0.5C for 2 minutes (4 seconds period)", + "Charge at 0.5C for 2 minutes (4 seconds period)", + ), ] ) - sim = model.predict(experiment=experiment) + sim = model.predict(experiment=experiment, inputs=x) return sim diff --git a/tests/unit/test_cost.py b/tests/unit/test_cost.py index f68df92b..3c0d8151 100644 --- a/tests/unit/test_cost.py +++ b/tests/unit/test_cost.py @@ -14,14 +14,17 @@ def model(self): return pybop.lithium_ion.SPM() @pytest.fixture - def parameters(self): - return [ - pybop.Parameter( - "Negative electrode active material volume fraction", - prior=pybop.Gaussian(0.5, 0.01), - bounds=[0.375, 0.625], - ), - ] + def ground_truth(self): + return 0.52 + + @pytest.fixture + def parameters(self, ground_truth): + return pybop.Parameter( + "Negative electrode active material volume fraction", + prior=pybop.Gaussian(0.5, 0.01), + bounds=[0.375, 0.625], + initial_value=ground_truth, + ) @pytest.fixture def experiment(self): @@ -32,15 +35,11 @@ def experiment(self): ) @pytest.fixture - def x0(self): - return np.array([0.52]) - - @pytest.fixture - def dataset(self, model, experiment, x0): + def dataset(self, model, experiment, ground_truth): model.parameter_set = model.pybamm_model.default_parameter_values model.parameter_set.update( { - "Negative electrode active material volume fraction": x0[0], + "Negative electrode active material volume fraction": ground_truth, } ) solution = model.predict(experiment=experiment) @@ -57,11 +56,11 @@ def signal(self): return "Voltage [V]" @pytest.fixture(params=[2.5, 3.777]) - def problem(self, model, parameters, dataset, signal, x0, request): + def problem(self, model, parameters, dataset, signal, request): cut_off = request.param model.parameter_set.update({"Lower voltage cut-off [V]": cut_off}) problem = pybop.FittingProblem( - model, parameters, dataset, signal=signal, x0=x0, init_soc=1.0 + model, parameters, dataset, signal=signal, init_soc=1.0 ) problem.dataset = dataset # add this to pass the pybop dataset to cost return problem @@ -81,7 +80,7 @@ def cost(self, problem, request): elif cls in [pybop.MAP]: return cls(problem, pybop.GaussianLogLikelihoodKnownSigma) elif cls in [pybop.ObserverCost]: - inputs = {p.name: problem.x0[i] for i, p in enumerate(problem.parameters)} + inputs = problem.parameters.initial_value() state = problem._model.reinit(inputs) n = len(state) sigma_diag = [0.0] * n @@ -114,12 +113,6 @@ def test_base(self, problem): with pytest.raises(NotImplementedError): base_cost.evaluateS1([0.5]) - @pytest.mark.unit - def test_design_base(self, problem): - design_cost = pybop.DesignCost(problem) - with pytest.raises(NotImplementedError): - design_cost([0.5]) - @pytest.mark.unit def test_MAP(self, problem): # Incorrect likelihood @@ -190,10 +183,14 @@ def test_costs(self, cost): @pytest.mark.parametrize( "cost_class", - [pybop.GravimetricEnergyDensity, pybop.VolumetricEnergyDensity], + [ + pybop.DesignCost, + pybop.GravimetricEnergyDensity, + pybop.VolumetricEnergyDensity, + ], ) @pytest.mark.unit - def test_energy_density_costs( + def test_design_costs( self, cost_class, model, @@ -209,21 +206,27 @@ def test_energy_density_costs( # Construct Cost cost = cost_class(problem) - # Test type of returned value - assert np.isscalar(cost([0.5])) - assert cost([0.4]) >= 0 # Should be a viable design - assert cost([0.8]) == -np.inf # Should exceed active material + porosity < 1 - assert cost([1.4]) == -np.inf # Definitely not viable - assert cost([-0.1]) == -np.inf # Should not be a viable design + if cost_class in [pybop.DesignCost]: + with pytest.raises(NotImplementedError): + cost([0.5]) + else: + # Test type of returned value + assert np.isscalar(cost([0.5])) + assert cost([0.4]) >= 0 # Should be a viable design + assert ( + cost([0.8]) == -np.inf + ) # Should exceed active material + porosity < 1 + assert cost([1.4]) == -np.inf # Definitely not viable + assert cost([-0.1]) == -np.inf # Should not be a viable design - # Test infeasible locations - cost.problem._model.allow_infeasible_solutions = False - assert cost([1.1]) == -np.inf + # Test infeasible locations + cost.problem._model.allow_infeasible_solutions = False + assert cost([1.1]) == -np.inf - # Test exception for non-numeric inputs - with pytest.raises(ValueError): - cost(["StringInputShouldNotWork"]) + # Test exception for non-numeric inputs + with pytest.raises(ValueError): + cost(["StringInputShouldNotWork"]) - # Compute after updating nominal capacity - cost = cost_class(problem, update_capacity=True) - cost([0.4]) + # Compute after updating nominal capacity + cost = cost_class(problem, update_capacity=True) + cost([0.4]) diff --git a/tests/unit/test_dataset.py b/tests/unit/test_dataset.py index fd58e8bc..9c4eac13 100644 --- a/tests/unit/test_dataset.py +++ b/tests/unit/test_dataset.py @@ -32,9 +32,13 @@ def test_dataset(self): assert np.all(dataset["Time [s]"] == solution["Time [s]"].data) # Test exception for non-dictionary inputs - with pytest.raises(ValueError): + with pytest.raises( + TypeError, match="The input to pybop.Dataset must be a dictionary." + ): pybop.Dataset(["StringInputShouldNotWork"]) - with pytest.raises(ValueError): + with pytest.raises( + TypeError, match="The input to pybop.Dataset must be a dictionary." + ): pybop.Dataset(solution["Time [s]"].data) # Test conversion of pybamm solution into dictionary diff --git a/tests/unit/test_likelihoods.py b/tests/unit/test_likelihoods.py index a590808c..41ee3667 100644 --- a/tests/unit/test_likelihoods.py +++ b/tests/unit/test_likelihoods.py @@ -14,14 +14,17 @@ def model(self): return pybop.lithium_ion.SPM() @pytest.fixture - def parameters(self): - return [ - pybop.Parameter( - "Negative electrode active material volume fraction", - prior=pybop.Gaussian(0.5, 0.01), - bounds=[0.375, 0.625], - ), - ] + def ground_truth(self): + return 0.52 + + @pytest.fixture + def parameters(self, ground_truth): + return pybop.Parameter( + "Negative electrode active material volume fraction", + prior=pybop.Gaussian(0.5, 0.01), + bounds=[0.375, 0.625], + initial_value=ground_truth, + ) @pytest.fixture def experiment(self): @@ -32,15 +35,11 @@ def experiment(self): ) @pytest.fixture - def x0(self): - return np.array([0.52]) - - @pytest.fixture - def dataset(self, model, experiment, x0): + def dataset(self, model, experiment, ground_truth): model.parameter_set = model.pybamm_model.default_parameter_values model.parameter_set.update( { - "Negative electrode active material volume fraction": x0[0], + "Negative electrode active material volume fraction": ground_truth, } ) solution = model.predict(experiment=experiment) @@ -53,17 +52,17 @@ def dataset(self, model, experiment, x0): ) @pytest.fixture - def one_signal_problem(self, model, parameters, dataset, x0): + def one_signal_problem(self, model, parameters, dataset): signal = ["Voltage [V]"] return pybop.FittingProblem( - model, parameters, dataset, signal=signal, x0=x0, init_soc=1.0 + model, parameters, dataset, signal=signal, init_soc=1.0 ) @pytest.fixture - def two_signal_problem(self, model, parameters, dataset, x0): + def two_signal_problem(self, model, parameters, dataset): signal = ["Time [s]", "Voltage [V]"] return pybop.FittingProblem( - model, parameters, dataset, signal=signal, x0=x0, init_soc=1.0 + model, parameters, dataset, signal=signal, init_soc=1.0 ) @pytest.mark.parametrize( @@ -73,17 +72,13 @@ def two_signal_problem(self, model, parameters, dataset, x0): @pytest.mark.unit def test_base_likelihood_init(self, problem_name, n_outputs, request): problem = request.getfixturevalue(problem_name) - likelihood = pybop.BaseLikelihood(problem, sigma=np.array([0.2])) + likelihood = pybop.BaseLikelihood(problem) assert likelihood.problem == problem assert likelihood.n_outputs == n_outputs assert likelihood.n_time_data == problem.n_time_data - assert np.array_equal(likelihood.get_sigma(), np.array([0.2])) assert likelihood.x0 == problem.x0 - assert likelihood.bounds == problem.bounds - assert likelihood._n_parameters == 1 + assert likelihood.n_parameters == 1 assert np.array_equal(likelihood._target, problem._target) - with pytest.raises(ValueError): - likelihood.set_sigma("Test") @pytest.mark.unit def test_base_likelihood_call_raises_not_implemented_error( @@ -94,24 +89,22 @@ def test_base_likelihood_call_raises_not_implemented_error( likelihood(np.array([0.5, 0.5])) @pytest.mark.unit - def test_base_likelihood_set_get_sigma(self, one_signal_problem): - likelihood = pybop.BaseLikelihood(one_signal_problem) + def test_set_get_sigma(self, one_signal_problem): + likelihood = pybop.GaussianLogLikelihoodKnownSigma(one_signal_problem, 0.1) likelihood.set_sigma(np.array([0.3])) assert np.array_equal(likelihood.get_sigma(), np.array([0.3])) - @pytest.mark.unit - def test_base_likelihood_set_sigma_raises_value_error_for_negative_sigma( - self, one_signal_problem - ): - likelihood = pybop.BaseLikelihood(one_signal_problem) + with pytest.raises( + ValueError, + match="The GaussianLogLikelihoodKnownSigma cost requires sigma to be " + + "either a scalar value or an array with one entry per dimension.", + ): + pybop.GaussianLogLikelihoodKnownSigma(one_signal_problem, sigma=None) + + likelihood = pybop.GaussianLogLikelihoodKnownSigma(one_signal_problem, 0.1) with pytest.raises(ValueError): likelihood.set_sigma(np.array([-0.2])) - @pytest.mark.unit - def test_base_likelihood_get_n_parameters(self, one_signal_problem): - likelihood = pybop.BaseLikelihood(one_signal_problem) - assert likelihood.get_n_parameters() == 1 - @pytest.mark.unit def test_base_likelihood_n_parameters_property(self, one_signal_problem): likelihood = pybop.BaseLikelihood(one_signal_problem) diff --git a/tests/unit/test_models.py b/tests/unit/test_models.py index e3c9df9a..d8686690 100644 --- a/tests/unit/test_models.py +++ b/tests/unit/test_models.py @@ -184,7 +184,7 @@ def test_parameter_set_definition(self): @pytest.mark.unit def test_rebuild_geometric_parameters(self): parameter_set = pybop.ParameterSet.pybamm("Chen2020") - parameters = [ + parameters = pybop.Parameters( pybop.Parameter( "Positive particle radius [m]", prior=pybop.Gaussian(4.8e-06, 0.05e-06), @@ -197,7 +197,7 @@ def test_rebuild_geometric_parameters(self): bounds=[30e-06, 50e-06], initial_value=48e-06, ), - ] + ) model = pybop.lithium_ion.SPM(parameter_set=parameter_set) model.build(parameters=parameters) @@ -209,8 +209,8 @@ def test_rebuild_geometric_parameters(self): out_init = initial_built_model.predict(t_eval=t_eval) # Test that the model can be rebuilt with different geometric parameters - parameters[0].update(5e-06) - parameters[1].update(45e-06) + parameters["Positive particle radius [m]"].update(5e-06) + parameters["Negative electrode thickness [m]"].update(45e-06) model.rebuild(parameters=parameters) rebuilt_model = model assert rebuilt_model._built_model is not None @@ -253,11 +253,14 @@ def test_reinit(self): state = model.reinit(inputs={}) np.testing.assert_array_almost_equal(state.as_ndarray(), np.array([[y0]])) - state = model.reinit(inputs=[]) + model.parameters = pybop.Parameters(pybop.Parameter("y0")) + state = model.reinit(inputs=[1]) np.testing.assert_array_almost_equal(state.as_ndarray(), np.array([[y0]])) model = ExponentialDecay(pybamm.ParameterValues({"k": k, "y0": y0})) - with pytest.raises(ValueError): + with pytest.raises( + ValueError, match="Model must be built before calling reinit" + ): model.reinit(inputs={}) @pytest.mark.unit @@ -302,7 +305,7 @@ def test_check_params(self): @pytest.mark.unit def test_non_converged_solution(self): model = pybop.lithium_ion.DFN() - parameters = [ + parameters = pybop.Parameters( pybop.Parameter( "Negative electrode active material volume fraction", prior=pybop.Gaussian(0.2, 0.01), @@ -311,7 +314,7 @@ def test_non_converged_solution(self): "Positive electrode active material volume fraction", prior=pybop.Gaussian(0.2, 0.01), ), - ] + ) dataset = pybop.Dataset( { "Time [s]": np.linspace(0, 100, 100), diff --git a/tests/unit/test_observer_unscented_kalman.py b/tests/unit/test_observer_unscented_kalman.py index b886f61b..2a947e71 100644 --- a/tests/unit/test_observer_unscented_kalman.py +++ b/tests/unit/test_observer_unscented_kalman.py @@ -25,26 +25,24 @@ def model(self, request): @pytest.fixture def parameters(self): - return [ + return pybop.Parameters( pybop.Parameter( "k", prior=pybop.Gaussian(0.1, 0.05), bounds=[0, 1], + initial_value=0.1, ), pybop.Parameter( "y0", prior=pybop.Gaussian(1, 0.05), bounds=[0, 3], + initial_value=1.0, ), - ] - - @pytest.fixture - def x0(self): - return np.array([0.1, 1.0]) + ) @pytest.fixture - def dataset(self, model: pybop.BaseModel, parameters, x0): - observer = pybop.Observer(parameters, model, signal=["2y"], x0=x0) + def dataset(self, model: pybop.BaseModel, parameters): + observer = pybop.Observer(parameters, model, signal=["2y"]) measurements = [] t_eval = np.linspace(0, 20, 10) for t in t_eval: @@ -57,7 +55,7 @@ def dataset(self, model: pybop.BaseModel, parameters, x0): return {"Time [s]": t_eval, "y": measurements} @pytest.fixture - def observer(self, model: pybop.BaseModel, parameters, x0): + def observer(self, model: pybop.BaseModel, parameters): n = model.n_states sigma0 = np.diag([self.measure_noise] * n) process = np.diag([1e-6] * n) @@ -69,7 +67,7 @@ def observer(self, model: pybop.BaseModel, parameters, x0): process[1, 1] = 0 measure = np.diag([1e-4]) observer = pybop.UnscentedKalmanFilterObserver( - parameters, model, sigma0, process, measure, signal=["2y"], x0=x0 + parameters, model, sigma0, process, measure, signal=["2y"] ) return observer diff --git a/tests/unit/test_observers.py b/tests/unit/test_observers.py index 9e77682d..46987bae 100644 --- a/tests/unit/test_observers.py +++ b/tests/unit/test_observers.py @@ -22,29 +22,29 @@ def model(self, request): @pytest.fixture def parameters(self): - return [ + return pybop.Parameters( pybop.Parameter( "k", prior=pybop.Gaussian(0.1, 0.05), bounds=[0, 1], + initial_value=0.1, ), pybop.Parameter( "y0", prior=pybop.Gaussian(1, 0.05), bounds=[0, 3], + initial_value=1.0, ), - ] - - @pytest.fixture - def x0(self): - return np.array([0.1, 1.0]) + ) @pytest.mark.unit - def test_observer(self, model, parameters, x0): + def test_observer(self, model, parameters): n = model.n_states - observer = pybop.Observer(parameters, model, signal=["2y"], x0=x0) + observer = pybop.Observer(parameters, model, signal=["2y"]) t_eval = np.linspace(0, 1, 100) - expected = x0[1] * np.exp(-x0[0] * t_eval) + expected = parameters["y0"].initial_value * np.exp( + -parameters["k"].initial_value * t_eval + ) for y, t in zip(expected, t_eval): observer.observe(t) np.testing.assert_array_almost_equal( @@ -72,7 +72,7 @@ def test_observer(self, model, parameters, x0): # Test evaluate with different inputs observer._time_data = t_eval - observer.evaluate(x0) + observer.evaluate(parameters.initial_value()) observer.evaluate(parameters) # Test evaluate with dataset @@ -83,7 +83,7 @@ def test_observer(self, model, parameters, x0): } ) observer._target = {"2y": expected} - observer.evaluate(x0) + observer.evaluate(parameters.initial_value()) @pytest.mark.unit def test_unbuilt_model(self, parameters): diff --git a/tests/unit/test_optimisation.py b/tests/unit/test_optimisation.py index 2f407885..c8494e09 100644 --- a/tests/unit/test_optimisation.py +++ b/tests/unit/test_optimisation.py @@ -23,17 +23,15 @@ def dataset(self): @pytest.fixture def one_parameter(self): - return [ - pybop.Parameter( - "Negative electrode active material volume fraction", - prior=pybop.Gaussian(0.6, 0.2), - bounds=[0.58, 0.62], - ) - ] + return pybop.Parameter( + "Negative electrode active material volume fraction", + prior=pybop.Gaussian(0.6, 0.2), + bounds=[0.58, 0.62], + ) @pytest.fixture def two_parameters(self): - return [ + return pybop.Parameters( pybop.Parameter( "Negative electrode active material volume fraction", prior=pybop.Gaussian(0.6, 0.2), @@ -44,7 +42,7 @@ def two_parameters(self): prior=pybop.Gaussian(0.55, 0.05), bounds=[0.53, 0.57], ), - ] + ) @pytest.fixture def model(self): @@ -122,6 +120,7 @@ def test_optimiser_classes(self, two_param_cost, optimiser, expected_name): @pytest.mark.unit def test_optimiser_kwargs(self, cost, optimiser): optim = optimiser(cost=cost, maxiter=1) + cost_bounds = cost.parameters.get_bounds() # Check maximum iterations optim.run() @@ -129,10 +128,10 @@ def test_optimiser_kwargs(self, cost, optimiser): if optimiser in [pybop.GradientDescent, pybop.Adam, pybop.NelderMead]: # Ignored bounds - optim = optimiser(cost=cost, bounds=cost.bounds) + optim = optimiser(cost=cost, bounds=cost_bounds) assert optim.bounds is None elif optimiser in [pybop.PSO]: - assert optim.bounds == cost.bounds + assert optim.bounds == cost_bounds # Cannot accept infinite bounds bounds = {"upper": [np.inf], "lower": [0.57]} with pytest.raises( @@ -142,7 +141,7 @@ def test_optimiser_kwargs(self, cost, optimiser): optim = optimiser(cost=cost, bounds=bounds) else: # Check and update bounds - assert optim.bounds == cost.bounds + assert optim.bounds == cost_bounds bounds = {"upper": [0.63], "lower": [0.57]} optim = optimiser(cost=cost, bounds=bounds) assert optim.bounds == bounds @@ -154,8 +153,10 @@ def test_optimiser_kwargs(self, cost, optimiser): parallel=False, min_iterations=3, max_unchanged_iterations=5, - threshold=1e-2, + absolute_tolerance=1e-2, + relative_tolerance=1e-4, max_evaluations=20, + threshold=1e-4, ) with pytest.warns( UserWarning, @@ -289,8 +290,7 @@ def test_prior_sampling(self, cost): # Tests prior sampling for i in range(50): optim = pybop.Optimisation(cost=cost) - - assert optim.x0 <= 0.62 and optim.x0 >= 0.58 + assert optim.x0[0] < 0.62 and optim.x0[0] > 0.58 @pytest.mark.unit @pytest.mark.parametrize( @@ -311,7 +311,7 @@ def test_scipy_prior_resampling( ) # Define the problem and cost - problem = pybop.FittingProblem(model, [parameter], dataset) + problem = pybop.FittingProblem(model, parameter, dataset) cost = pybop.SumSquaredError(problem) # Create the optimisation class with infeasible solutions disabled @@ -352,19 +352,21 @@ def test_halting(self, cost): # Test invalid values with pytest.raises(ValueError): optim.set_max_iterations(-1) - with pytest.raises(ValueError): - optim.set_max_evaluations(-1) with pytest.raises(ValueError): optim.set_min_iterations(-1) with pytest.raises(ValueError): optim.set_max_unchanged_iterations(-1) with pytest.raises(ValueError): - optim.set_max_unchanged_iterations(1, threshold=-1) + optim.set_max_unchanged_iterations(1, absolute_tolerance=-1) + with pytest.raises(ValueError): + optim.set_max_unchanged_iterations(1, relative_tolerance=-1) + with pytest.raises(ValueError): + optim.set_max_evaluations(-1) optim = pybop.Optimisation(cost=cost) # Trigger threshold - optim._threshold = np.inf + optim.set_threshold(np.inf) optim.run() optim.set_max_unchanged_iterations() diff --git a/tests/unit/test_parameter_sets.py b/tests/unit/test_parameter_sets.py index f24e7bea..a71a28d2 100644 --- a/tests/unit/test_parameter_sets.py +++ b/tests/unit/test_parameter_sets.py @@ -74,7 +74,7 @@ def test_ecm_parameter_sets(self): assert params() == params.params # Test exporting a json file - parameters = [ + parameters = pybop.Parameters( pybop.Parameter( "R0 [Ohm]", prior=pybop.Gaussian(0.0002, 0.0001), @@ -87,7 +87,7 @@ def test_ecm_parameter_sets(self): bounds=[1e-5, 1e-2], initial_value=0.0002, ), - ] + ) params.export_parameters( "examples/scripts/parameters/fit_ecm_parameters.json", fit_params=parameters ) diff --git a/tests/unit/test_parameters.py b/tests/unit/test_parameters.py index 19649cb9..195fbdef 100644 --- a/tests/unit/test_parameters.py +++ b/tests/unit/test_parameters.py @@ -3,7 +3,7 @@ import pybop -class TestParameters: +class TestParameter: """ A class to test the parameter classes. """ @@ -65,13 +65,95 @@ def test_no_bounds(self): @pytest.mark.unit def test_invalid_inputs(self, parameter): # Test error with invalid value - with pytest.raises(ValueError): + with pytest.raises(ValueError, match="Margin must be between 0 and 1"): parameter.set_margin(margin=-1) # Test error with no parameter value - with pytest.raises(ValueError): + with pytest.raises(ValueError, match="No value provided to update parameter"): parameter.update() # Test error with opposite bounds - with pytest.raises(ValueError): + with pytest.raises( + ValueError, match="Lower bound must be less than upper bound" + ): pybop.Parameter("Name", bounds=[0.7, 0.3]) + + +class TestParameters: + """ + A class to test the parameter classes. + """ + + @pytest.fixture + def parameter(self): + return pybop.Parameter( + "Negative electrode active material volume fraction", + prior=pybop.Gaussian(0.6, 0.02), + bounds=[0.375, 0.7], + initial_value=0.6, + ) + + @pytest.mark.unit + def test_parameters_construction(self, parameter): + params = pybop.Parameters(parameter) + assert parameter.name in params.param.keys() + assert parameter in params.param.values() + + # Test parameter addition via Parameter class + params = pybop.Parameters() # empty + params.add(parameter) + assert parameter.name in params.param.keys() + assert parameter in params.param.values() + + with pytest.raises( + ValueError, + match="There is already a parameter with the name " + + "Negative electrode active material volume fraction" + + " in the Parameters object. Please remove the duplicate entry.", + ): + params.add(parameter) + + params.remove(parameter_name=parameter.name) + + # Test parameter addition via dict + params.add( + dict( + name="Negative electrode active material volume fraction", + initial_value=0.6, + ) + ) + with pytest.raises( + ValueError, + match="There is already a parameter with the name " + + "Negative electrode active material volume fraction" + + " in the Parameters object. Please remove the duplicate entry.", + ): + params.add( + dict( + name="Negative electrode active material volume fraction", + initial_value=0.6, + ) + ) + + params.remove(parameter_name=parameter.name) + with pytest.raises( + ValueError, match="This parameter does not exist in the Parameters object." + ): + params.remove(parameter_name=parameter.name) + + with pytest.raises( + TypeError, match="Each parameter input must be a Parameter or a dictionary." + ): + params.add(parameter="Invalid string") + with pytest.raises( + TypeError, match="The input parameter_name is not a string." + ): + params.remove(parameter_name=parameter) + + @pytest.mark.unit + def test_get_sigma(self, parameter): + params = pybop.Parameters(parameter) + assert params.get_sigma0() == [0.02] + + parameter.prior = None + assert params.get_sigma0() is None diff --git a/tests/unit/test_plots.py b/tests/unit/test_plots.py index f5998c0a..e36b8ba8 100644 --- a/tests/unit/test_plots.py +++ b/tests/unit/test_plots.py @@ -28,7 +28,7 @@ def model(self): @pytest.fixture def parameters(self): - return [ + return pybop.Parameters( pybop.Parameter( "Negative electrode active material volume fraction", prior=pybop.Gaussian(0.68, 0.05), @@ -39,7 +39,7 @@ def parameters(self): prior=pybop.Gaussian(0.58, 0.05), bounds=[0.4, 0.7], ), - ] + ) @pytest.fixture def dataset(self, model): diff --git a/tests/unit/test_problem.py b/tests/unit/test_problem.py index 4fe2e769..9af00164 100644 --- a/tests/unit/test_problem.py +++ b/tests/unit/test_problem.py @@ -16,7 +16,7 @@ def model(self): @pytest.fixture def parameters(self): - return [ + return pybop.Parameters( pybop.Parameter( "Negative particle radius [m]", prior=pybop.Gaussian(2e-05, 0.1e-5), @@ -27,7 +27,7 @@ def parameters(self): prior=pybop.Gaussian(0.5e-05, 0.1e-5), bounds=[1e-6, 5e-5], ), - ] + ) @pytest.fixture def experiment(self): @@ -68,10 +68,6 @@ def signal(self): @pytest.mark.unit def test_base_problem(self, parameters, model, dataset): - # Test incorrect number of initial parameter values - with pytest.raises(ValueError): - pybop.BaseProblem(parameters, model=model, x0=np.array([])) - # Construct Problem problem = pybop.BaseProblem(parameters, model=model) @@ -85,12 +81,6 @@ def test_base_problem(self, parameters, model, dataset): with pytest.raises(ValueError): pybop.BaseProblem(parameters, model=model, signal=[1e-5, 1e-5]) - # Test without bounds - for param in parameters: - param.bounds = None - problem = pybop.BaseProblem(parameters, model=model) - assert problem.bounds is None - # Incorrect set target with pytest.raises(ValueError, match="Dataset must be a pybop Dataset object."): problem.set_target("This is not a dataset") @@ -100,14 +90,18 @@ def test_base_problem(self, parameters, model, dataset): with pytest.raises(ValueError, match="Signal must be defined to set target."): problem.set_target(dataset) + # Different types of parameters + parameter_list = list(parameters.param.values()) + problem = pybop.BaseProblem(parameters=parameter_list) + problem = pybop.BaseProblem(parameters=parameter_list[0]) + with pytest.raises( + TypeError, + match="The input parameters must be a pybop Parameter, a list of pybop.Parameter objects, or a pybop Parameters object.", + ): + problem = pybop.BaseProblem(parameters="Invalid string") + @pytest.mark.unit def test_fitting_problem(self, parameters, dataset, model, signal): - # Test incorrect number of initial parameter values - with pytest.raises(ValueError): - pybop.FittingProblem( - model, parameters, dataset, signal=signal, x0=np.array([]) - ) - # Construct Problem problem = pybop.FittingProblem(model, parameters, dataset, signal=signal) @@ -163,10 +157,6 @@ def test_fitting_problem(self, parameters, dataset, model, signal): @pytest.mark.unit def test_design_problem(self, parameters, experiment, model): - # Test incorrect number of initial parameter values - with pytest.raises(ValueError): - pybop.DesignProblem(model, parameters, experiment, x0=np.array([])) - # Construct Problem problem = pybop.DesignProblem(model, parameters, experiment) @@ -184,6 +174,7 @@ def test_problem_construct_with_model_predict( self, parameters, model, dataset, signal ): # Construct model and predict + model.parameters = parameters out = model.predict(inputs=[1e-5, 1e-5], t_eval=np.linspace(0, 10, 100)) problem = pybop.FittingProblem( @@ -198,5 +189,5 @@ def test_problem_construct_with_model_predict( assert_allclose( out["Voltage [V]"].data, problem_output["Voltage [V]"], - atol=1e-5, + atol=1e-6, ) diff --git a/tests/unit/test_standalone.py b/tests/unit/test_standalone.py index d524555a..a164312f 100644 --- a/tests/unit/test_standalone.py +++ b/tests/unit/test_standalone.py @@ -42,7 +42,7 @@ def test_optimisation_on_standalone_cost(self): @pytest.mark.unit def test_standalone_problem(self): # Define parameters to estimate - parameters = [ + parameters = pybop.Parameters( pybop.Parameter( "Gradient", prior=pybop.Gaussian(4.2, 0.02), @@ -53,7 +53,7 @@ def test_standalone_problem(self): prior=pybop.Gaussian(3.3, 0.02), bounds=[-1, 10], ), - ] + ) # Define target data t_eval = np.linspace(0, 1, 100) @@ -84,10 +84,6 @@ def test_standalone_problem(self): np.testing.assert_allclose(x[0], 934.006734006734, atol=1e-2) np.testing.assert_allclose(x[1], [-334.006734, 0.0], atol=1e-2) - # Test incorrect number of initial parameter values - with pytest.raises(ValueError): - StandaloneProblem(parameters, dataset, signal=signal, x0=np.array([])) - # Test problem construction errors for bad_dataset in [ pybop.Dataset({"Time [s]": np.array([0])}),