From 9e2027effd4efbb46adc402bcc2058eaaf8605b7 Mon Sep 17 00:00:00 2001 From: derek-slack Date: Wed, 5 Aug 2026 15:53:56 -0400 Subject: [PATCH 1/9] Add FoKL GPs for parameter estimation --- .../SPMe_FoKLGPy_Fitting.py | 75 +++++ pybop/__init__.py | 1 + pybop/parameters/gp_parameter.py | 300 ++++++++++++++++++ pyproject.toml | 5 +- 4 files changed, 380 insertions(+), 1 deletion(-) create mode 100644 examples/scripts/battery_parameterisation/SPMe_FoKLGPy_Fitting.py create mode 100644 pybop/parameters/gp_parameter.py diff --git a/examples/scripts/battery_parameterisation/SPMe_FoKLGPy_Fitting.py b/examples/scripts/battery_parameterisation/SPMe_FoKLGPy_Fitting.py new file mode 100644 index 000000000..19fb6f5ba --- /dev/null +++ b/examples/scripts/battery_parameterisation/SPMe_FoKLGPy_Fitting.py @@ -0,0 +1,75 @@ +import numpy as np +import pybamm + +import pybop +import matplotlib.pyplot as plt +import pandas as pd +from pathlib import Path + +""" +In this example, we demonstrate the use of FoKL decomposed Gaussian Processes for estimation of +some parameters. +""" + + +# Define model and parameter values +model = pybamm.lithium_ion.DFN() + + +parameter_values = pybamm.ParameterValues("Chen2020") + +# Generate a synthetic dataset +sim = pybamm.Simulation(model, parameter_values=parameter_values) +t_eval = np.linspace(0, 3400, 240) +solution = sim.solve(t_eval=t_eval) +original = solution["Positive electrode exchange current density [A.m-2]"].entries + +sigma = 5e-3 +dataset = pybop.Dataset( + { + "Time [s]": t_eval, + "Current [A]": solution["Current [A]"](t_eval), + "Voltage [V]": pybop.add_noise(solution["Voltage [V]"](t_eval), sigma), + "Bulk open-circuit voltage [V]": solution["Bulk open-circuit voltage [V]"]( + t_eval + ), + } +) + +# Create GP terms +GP_options = {'Number of terms':5,'Constant mean':8,'div_arg':[[1,2],[0,2]], 'exp':False} +GP_param_neg = pybop.FoKLGP("Positive electrode exchange-current density [A.m-2]", parameter_values=parameter_values, options=GP_options, twoway=True) +new_parameters = GP_param_neg.get_parameter_values() + +simulator = pybop.pybamm.Simulator(model, new_parameters, protocol=dataset) +target = ["Voltage [V]", "Bulk open-circuit voltage [V]"] +cost = pybop.RootMeanSquaredError(dataset, target=target) +problem = pybop.Problem(simulator, cost) + +# Set up the optimiser +options = pybop.PintsOptions( + max_iterations=300, + max_unchanged_iterations=50, + verbose=True +) +optim = pybop.IRPropPlus(problem, options=options) + +# Run the optimisation +result = optim.run() +print(result) + +# Plot the timeseries output +pybop.plot.problem(problem, inputs=result.best_inputs, title="Optimised Comparison") + +# +# sim = pybamm.Simulation(model, parameter_values=new_parameters) +# +# solution = sim.solve(t_eval=t_eval, inputs=result.best_inputs) +# solution.all_inputs = [result.best_inputs] +# new = solution["Positive electrode exchange current density [A.m-2]"].entries +# +# print(np.sum((new-original)**2)/np.size(original)) + +# Plot the optimisation result +result.plot_parameters() + diff --git a/pybop/__init__.py b/pybop/__init__.py index 11436f11c..eb0091ea2 100644 --- a/pybop/__init__.py +++ b/pybop/__init__.py @@ -51,6 +51,7 @@ from .parameters.parameter import Parameter, Parameters from .parameters.distributions import Distribution, Exponential, Gaussian, JointDistribution, LogNormal, LogUniform, Unbounded, Uniform from .parameters.multivariate_distributions import MultivariateNonparametric, MultivariateUniform, MultivariateGaussian, MultivariateLogNormal, MarginalDistribution +from .parameters.gp_parameter import FoKLGP # # Model classes diff --git a/pybop/parameters/gp_parameter.py b/pybop/parameters/gp_parameter.py new file mode 100644 index 000000000..2f779471f --- /dev/null +++ b/pybop/parameters/gp_parameter.py @@ -0,0 +1,300 @@ +import numpy as np +import pybamm + +import FoKL +from FoKL.getKernels import sp500, bernoulli +import itertools + +import pybop +from pybop.parameters import parameter + +class FoKLGP: + """ + Creates parameter functions as decomposed GPs + options: [dict] + """ + def __init__(self,name, options=None, parameter_values=None, kernel='Bernoulli Polynomial', twoway=True): + GP_dict_list = self._process_options(name, options, parameter_values) + self.evaluate_func = self._evaluate_parameter(kernel) + self.damtx = self._create_interaction_matrix(GP_dict_list['Number of terms'], GP_dict_list['Number of inputs'], twoway) + self.add_to_params(GP_dict_list,self.damtx) + def __call__(self, *args): + return self.pybamm_function(*args) + + def _process_options(self, name, options, parameter_values): + """ + + """ + default_options = {'arg_inds':None, 'exp':True, 'Number of inputs':1, 'div_arg':None, 'div_const':None, + 'Number of terms':1, 'Constant standard deviation':0.5,'Bi mean':0, 'Bi standard deviation':0.5} + default_options.update({'Name':name}) + if options is not None: + default_options.update(options) + if 'Constant mean' not in default_options: + # If no beta constant term distribution described grab from parameters, if this is a function then user supplied + try: + if default_options['exp']: + B0_mean = np.log(parameter_values[name]) + + else: + B0_mean = parameter_values[name] + default_options.update({'Constant mean': B0_mean}) + except: + raise ValueError(f'Default parameter value for {name} is not a constant, please supply an estimate') + num_inputs = 0 + if default_options['arg_inds'] is not None: + num_inputs += len(default_options['arg_inds']) + if default_options['div_arg'] is not None: + num_inputs += len(default_options['div_arg']) + + default_options['Number of inputs']=num_inputs + self.parameter_values = parameter_values + return default_options + + def _create_interaction_matrix(self, number_of_terms, number_of_inputs, twoway, damtx=[]): + def perms(x): + """Python equivalent of MATLAB perms.""" + a = np.array(np.vstack(list(itertools.permutations(x)))[::-1]) + return a + + if twoway: + sett = 2 + else: + sett = 1 + + + for ind in range(1, number_of_terms+1): + indvec = np.zeros((number_of_inputs)) + summ = ind + while summ: + for j in range(0, sett): + indvec[j] = indvec[j] + 1 + summ = summ - 1 + if summ == 0: + break + + vecs = np.unique(perms(indvec), axis=0) + + if np.size(damtx) == 0: + damtx = vecs + else: + damtx = np.append(damtx, vecs, axis=0) + + print(damtx) + return damtx.astype(int) + + def _set_kernel(self, kernel = 'Bernoulli Polynomial'): + if kernel == 'Cubic Splines': + self.phis = sp500() + elif kernel == 'Bernoulli Polynomial': + self.phis = bernoulli() + + def _evaluate_parameter(self, kernel): + """ + The symbolic evaluation of the decomposed GP. + + kernel: Kernel function for evaluation, must be `Cubic Splines` or `Bernoulli Polynomial` + """ + self._set_kernel(kernel) + + if kernel == 'Cubic Splines': + def evaluate_pybamm( + betas, + mtx, + inputs, + coeff=None): + + + num_basis_terms = len(mtx) + num_inputs = len(mtx[0]) + X_sol = [] + + mtx = np.array(mtx) + phind = [] + for i in range(num_inputs): + phind_temp = inputs[i] * 499 + sett = (phind_temp == 0) + phind_temp = phind_temp + sett + r = 1 / 499 # interval of when basis function changes (i.e., when next cubic function defines spline) + phind.append(phind_temp - 1) + + A = [1, 2, 3] + + X_sc = [(1 - inputs[0]) ** a for a in A] + + lspace = [] + for i in range(num_inputs): + lspace.append(np.linspace(0, 499, 499)) + lspace = np.array(lspace) + + for j in range(num_basis_terms): + phi = 1 + for k in range(num_inputs): + num = mtx[j][k] + + if num > 0: + nid = int(num - 1) + if coeff is None: + coeff = [] + for jj in range(4): + phispace = self.phis[nid][jj].reshape(1, -1) + phi_interp = pybamm.Interpolant(lspace[0], phispace[0], + phind[k]) + coeff.append(phi_interp) + + phi *= coeff[0] + coeff[1] * X_sc[0] + coeff[2] * X_sc[1] + coeff[3] * X_sc[2] + X_sol.append(phi) + + X_sol_ones = betas[0] + mean = X_sol_ones + for i in range(len(X_sol)): + X_sol_betas = X_sol[i] * betas[i + 1] + mean += X_sol_betas + + return mean + elif kernel == 'Bernoulli Polynomial': + def evaluate_pybamm(betas, + mtx, + inputs, + coeff=None): + """ + Pybamm Function evaluation + betas: indexed from beta list that relates to + """ + n = 1 + num_basis_terms = len(mtx) + num_inputs = len(mtx[0]) + X_sol = [] + + mtx = np.array(mtx) + + def bernoulli_func(phis, num, x): + if num > 0: + coeff = phis[num - 1] + result = coeff[0] + sum(coeff[k] * (x ** k) for k in range(1, len(coeff))) + else: + result = 1.0 + return result + + for j in range(num_basis_terms): + phi = 1 + for k in range(num_inputs): + num = mtx[j][k] + phi *= bernoulli_func(self.phis, num, inputs[k][0]) + X_sol.append(phi) + + X_sol_ones = betas[0] + mean = X_sol_ones + for i in range(len(X_sol)): + X_sol_betas = X_sol[i] * betas[i + 1] + mean += X_sol_betas + return mean + else: + raise NotImplementedError("Kernel must be either `Cubic Splines` or `Bernoulli Polynomial`") + + return evaluate_pybamm + + def add_function(self, name, mtx, arg_inds, betas_function, exp=False, div_arg=None, div_const=None, + new_children=None): + """ + Creates Parameter function specified as a GP object + inputs: + name: str, Name of parameter to be estimated + arg_inds: list of int, Index of inputs to parameter function from PyBaMM + list_index: int, Index of function hyperparameters, betas and mtx, supplied in initialization + exp: Bool, if function should be exponential + div_arg: list of lists of int, optional, Index of arguments to be used as div + ex: div_arg = [[1,4],[3,2]] + inputs to GP would be GP((1/2),(3,2)) + """ + + beta_func = betas_function + + if div_arg: + if exp: + def pybamm_function(*args): + xs = [] + for x in div_arg: + xs.append([args[x[0]] / args[x[1]]]) + for x in arg_inds: + xs.append([args[x]]) + + res = np.exp(self.evaluate_func(beta_func, mtx, xs)) + return res + else: + def pybamm_function(*args): + xs = [] + for x in div_arg: + xs.append([args[x[0]] / args[x[1]]]) + for x in arg_inds: + xs.append([args[x]]) + + res = self.evaluate_func(beta_func, mtx, xs) + return res + else: + if exp: + def pybamm_function(*args): + xs = [] + for x in arg_inds: + if div_const: + xs.append([args[x] / div_const[0]]) + else: + xs.append([args[x]]) + + res = np.exp(self.evaluate_func(beta_func, mtx, xs)) + return res + else: + def pybamm_function(*args): + xs = [] + for x in arg_inds: + xs.append([args[x]]) + + res = self.evaluate_func(beta_func, mtx, xs) + return res + if type(self.parameter_values[name]) is not float: + function_args = self.parameter_values[name].__code__.co_varnames + function_args_mod = [] + if div_arg is not None: + for x in div_arg: + function_args_mod.append(function_args[x[0]] + str('/') + function_args[x[1]]) + for x in arg_inds: + function_args_mod.append(function_args[x]) + print(f"GP function created for {name} \n inputs are {function_args_mod}") + self.pybamm_function = pybamm_function + return pybamm_function + + + def create_beta_inputs(self,len_mtx, GP): + """ + Generates Input Variables + """ + betas_symbolic = [] + beta_parameters = {} + for i in range(len_mtx): + key_str = GP['Name'] + ' Beta ' + str(i) + betas_symbolic.append(pybamm.InputParameter(key_str)) + if i == 0: + beta_parameters[key_str] = pybop.Parameter( + distribution=pybop.Gaussian(GP['Constant mean'], GP['Constant standard deviation']), + ) + + else: + beta_parameters[key_str] = pybop.Parameter( + distribution=pybop.Gaussian(GP['Bi mean'], GP['Bi standard deviation']), + ) + + self.beta_parameters = beta_parameters + return betas_symbolic, beta_parameters + + def add_to_params(self, GP, damtxs): + + betas_function, beta_parameters = self.create_beta_inputs(len(damtxs) + 1, GP) + self.parameter_values.update(beta_parameters) + pybamm_function = self.add_function(GP['Name'], damtxs, GP['arg_inds'], betas_function, exp=GP['exp'], div_arg=GP['div_arg'], + div_const=GP['div_const']) + self.parameter_values.update({GP['Name']:pybamm_function}) + + def get_parameter_values(self): + return self.parameter_values + + diff --git a/pyproject.toml b/pyproject.toml index acd729723..a5515df04 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -51,7 +51,10 @@ ep-bolfi = [ pyprobe = [ "PyProBE-Data>=2.5.0;python_version >= '3.11' and python_version < '3.13'" ] -all = ["pybop[plot,salib,scifem,bpx,pyprobe,ep-bolfi]"] +fokl = [ + "FoKL>=1.2.0" +] +all = ["pybop[plot,salib,scifem,bpx,pyprobe,ep-bolfi,fokl]"] [dependency-groups] docs = [ From 96d54b34cca977ec9338722bb85d7a1e76c68184 Mon Sep 17 00:00:00 2001 From: derek-slack Date: Wed, 5 Aug 2026 15:53:56 -0400 Subject: [PATCH 2/9] Add FoKL GPs for parameter estimation --- .../SPMe_FoKLGPy_Fitting.py | 75 +++++ pybop/__init__.py | 1 + pybop/parameters/gp_parameter.py | 300 ++++++++++++++++++ pyproject.toml | 19 +- 4 files changed, 389 insertions(+), 6 deletions(-) create mode 100644 examples/scripts/battery_parameterisation/SPMe_FoKLGPy_Fitting.py create mode 100644 pybop/parameters/gp_parameter.py diff --git a/examples/scripts/battery_parameterisation/SPMe_FoKLGPy_Fitting.py b/examples/scripts/battery_parameterisation/SPMe_FoKLGPy_Fitting.py new file mode 100644 index 000000000..19fb6f5ba --- /dev/null +++ b/examples/scripts/battery_parameterisation/SPMe_FoKLGPy_Fitting.py @@ -0,0 +1,75 @@ +import numpy as np +import pybamm + +import pybop +import matplotlib.pyplot as plt +import pandas as pd +from pathlib import Path + +""" +In this example, we demonstrate the use of FoKL decomposed Gaussian Processes for estimation of +some parameters. +""" + + +# Define model and parameter values +model = pybamm.lithium_ion.DFN() + + +parameter_values = pybamm.ParameterValues("Chen2020") + +# Generate a synthetic dataset +sim = pybamm.Simulation(model, parameter_values=parameter_values) +t_eval = np.linspace(0, 3400, 240) +solution = sim.solve(t_eval=t_eval) +original = solution["Positive electrode exchange current density [A.m-2]"].entries + +sigma = 5e-3 +dataset = pybop.Dataset( + { + "Time [s]": t_eval, + "Current [A]": solution["Current [A]"](t_eval), + "Voltage [V]": pybop.add_noise(solution["Voltage [V]"](t_eval), sigma), + "Bulk open-circuit voltage [V]": solution["Bulk open-circuit voltage [V]"]( + t_eval + ), + } +) + +# Create GP terms +GP_options = {'Number of terms':5,'Constant mean':8,'div_arg':[[1,2],[0,2]], 'exp':False} +GP_param_neg = pybop.FoKLGP("Positive electrode exchange-current density [A.m-2]", parameter_values=parameter_values, options=GP_options, twoway=True) +new_parameters = GP_param_neg.get_parameter_values() + +simulator = pybop.pybamm.Simulator(model, new_parameters, protocol=dataset) +target = ["Voltage [V]", "Bulk open-circuit voltage [V]"] +cost = pybop.RootMeanSquaredError(dataset, target=target) +problem = pybop.Problem(simulator, cost) + +# Set up the optimiser +options = pybop.PintsOptions( + max_iterations=300, + max_unchanged_iterations=50, + verbose=True +) +optim = pybop.IRPropPlus(problem, options=options) + +# Run the optimisation +result = optim.run() +print(result) + +# Plot the timeseries output +pybop.plot.problem(problem, inputs=result.best_inputs, title="Optimised Comparison") + +# +# sim = pybamm.Simulation(model, parameter_values=new_parameters) +# +# solution = sim.solve(t_eval=t_eval, inputs=result.best_inputs) +# solution.all_inputs = [result.best_inputs] +# new = solution["Positive electrode exchange current density [A.m-2]"].entries +# +# print(np.sum((new-original)**2)/np.size(original)) + +# Plot the optimisation result +result.plot_parameters() + diff --git a/pybop/__init__.py b/pybop/__init__.py index 228d1e7f2..719231a9f 100644 --- a/pybop/__init__.py +++ b/pybop/__init__.py @@ -51,6 +51,7 @@ from .parameters.parameter import Parameter, Parameters from .parameters.distributions import Distribution, Exponential, Gaussian, JointDistribution, LogNormal, LogUniform, Unbounded, Uniform from .parameters.multivariate_distributions import MultivariateNonparametric, MultivariateUniform, MultivariateGaussian, MultivariateLogNormal, MarginalDistribution +from .parameters.gp_parameter import FoKLGP # # Model classes diff --git a/pybop/parameters/gp_parameter.py b/pybop/parameters/gp_parameter.py new file mode 100644 index 000000000..2f779471f --- /dev/null +++ b/pybop/parameters/gp_parameter.py @@ -0,0 +1,300 @@ +import numpy as np +import pybamm + +import FoKL +from FoKL.getKernels import sp500, bernoulli +import itertools + +import pybop +from pybop.parameters import parameter + +class FoKLGP: + """ + Creates parameter functions as decomposed GPs + options: [dict] + """ + def __init__(self,name, options=None, parameter_values=None, kernel='Bernoulli Polynomial', twoway=True): + GP_dict_list = self._process_options(name, options, parameter_values) + self.evaluate_func = self._evaluate_parameter(kernel) + self.damtx = self._create_interaction_matrix(GP_dict_list['Number of terms'], GP_dict_list['Number of inputs'], twoway) + self.add_to_params(GP_dict_list,self.damtx) + def __call__(self, *args): + return self.pybamm_function(*args) + + def _process_options(self, name, options, parameter_values): + """ + + """ + default_options = {'arg_inds':None, 'exp':True, 'Number of inputs':1, 'div_arg':None, 'div_const':None, + 'Number of terms':1, 'Constant standard deviation':0.5,'Bi mean':0, 'Bi standard deviation':0.5} + default_options.update({'Name':name}) + if options is not None: + default_options.update(options) + if 'Constant mean' not in default_options: + # If no beta constant term distribution described grab from parameters, if this is a function then user supplied + try: + if default_options['exp']: + B0_mean = np.log(parameter_values[name]) + + else: + B0_mean = parameter_values[name] + default_options.update({'Constant mean': B0_mean}) + except: + raise ValueError(f'Default parameter value for {name} is not a constant, please supply an estimate') + num_inputs = 0 + if default_options['arg_inds'] is not None: + num_inputs += len(default_options['arg_inds']) + if default_options['div_arg'] is not None: + num_inputs += len(default_options['div_arg']) + + default_options['Number of inputs']=num_inputs + self.parameter_values = parameter_values + return default_options + + def _create_interaction_matrix(self, number_of_terms, number_of_inputs, twoway, damtx=[]): + def perms(x): + """Python equivalent of MATLAB perms.""" + a = np.array(np.vstack(list(itertools.permutations(x)))[::-1]) + return a + + if twoway: + sett = 2 + else: + sett = 1 + + + for ind in range(1, number_of_terms+1): + indvec = np.zeros((number_of_inputs)) + summ = ind + while summ: + for j in range(0, sett): + indvec[j] = indvec[j] + 1 + summ = summ - 1 + if summ == 0: + break + + vecs = np.unique(perms(indvec), axis=0) + + if np.size(damtx) == 0: + damtx = vecs + else: + damtx = np.append(damtx, vecs, axis=0) + + print(damtx) + return damtx.astype(int) + + def _set_kernel(self, kernel = 'Bernoulli Polynomial'): + if kernel == 'Cubic Splines': + self.phis = sp500() + elif kernel == 'Bernoulli Polynomial': + self.phis = bernoulli() + + def _evaluate_parameter(self, kernel): + """ + The symbolic evaluation of the decomposed GP. + + kernel: Kernel function for evaluation, must be `Cubic Splines` or `Bernoulli Polynomial` + """ + self._set_kernel(kernel) + + if kernel == 'Cubic Splines': + def evaluate_pybamm( + betas, + mtx, + inputs, + coeff=None): + + + num_basis_terms = len(mtx) + num_inputs = len(mtx[0]) + X_sol = [] + + mtx = np.array(mtx) + phind = [] + for i in range(num_inputs): + phind_temp = inputs[i] * 499 + sett = (phind_temp == 0) + phind_temp = phind_temp + sett + r = 1 / 499 # interval of when basis function changes (i.e., when next cubic function defines spline) + phind.append(phind_temp - 1) + + A = [1, 2, 3] + + X_sc = [(1 - inputs[0]) ** a for a in A] + + lspace = [] + for i in range(num_inputs): + lspace.append(np.linspace(0, 499, 499)) + lspace = np.array(lspace) + + for j in range(num_basis_terms): + phi = 1 + for k in range(num_inputs): + num = mtx[j][k] + + if num > 0: + nid = int(num - 1) + if coeff is None: + coeff = [] + for jj in range(4): + phispace = self.phis[nid][jj].reshape(1, -1) + phi_interp = pybamm.Interpolant(lspace[0], phispace[0], + phind[k]) + coeff.append(phi_interp) + + phi *= coeff[0] + coeff[1] * X_sc[0] + coeff[2] * X_sc[1] + coeff[3] * X_sc[2] + X_sol.append(phi) + + X_sol_ones = betas[0] + mean = X_sol_ones + for i in range(len(X_sol)): + X_sol_betas = X_sol[i] * betas[i + 1] + mean += X_sol_betas + + return mean + elif kernel == 'Bernoulli Polynomial': + def evaluate_pybamm(betas, + mtx, + inputs, + coeff=None): + """ + Pybamm Function evaluation + betas: indexed from beta list that relates to + """ + n = 1 + num_basis_terms = len(mtx) + num_inputs = len(mtx[0]) + X_sol = [] + + mtx = np.array(mtx) + + def bernoulli_func(phis, num, x): + if num > 0: + coeff = phis[num - 1] + result = coeff[0] + sum(coeff[k] * (x ** k) for k in range(1, len(coeff))) + else: + result = 1.0 + return result + + for j in range(num_basis_terms): + phi = 1 + for k in range(num_inputs): + num = mtx[j][k] + phi *= bernoulli_func(self.phis, num, inputs[k][0]) + X_sol.append(phi) + + X_sol_ones = betas[0] + mean = X_sol_ones + for i in range(len(X_sol)): + X_sol_betas = X_sol[i] * betas[i + 1] + mean += X_sol_betas + return mean + else: + raise NotImplementedError("Kernel must be either `Cubic Splines` or `Bernoulli Polynomial`") + + return evaluate_pybamm + + def add_function(self, name, mtx, arg_inds, betas_function, exp=False, div_arg=None, div_const=None, + new_children=None): + """ + Creates Parameter function specified as a GP object + inputs: + name: str, Name of parameter to be estimated + arg_inds: list of int, Index of inputs to parameter function from PyBaMM + list_index: int, Index of function hyperparameters, betas and mtx, supplied in initialization + exp: Bool, if function should be exponential + div_arg: list of lists of int, optional, Index of arguments to be used as div + ex: div_arg = [[1,4],[3,2]] + inputs to GP would be GP((1/2),(3,2)) + """ + + beta_func = betas_function + + if div_arg: + if exp: + def pybamm_function(*args): + xs = [] + for x in div_arg: + xs.append([args[x[0]] / args[x[1]]]) + for x in arg_inds: + xs.append([args[x]]) + + res = np.exp(self.evaluate_func(beta_func, mtx, xs)) + return res + else: + def pybamm_function(*args): + xs = [] + for x in div_arg: + xs.append([args[x[0]] / args[x[1]]]) + for x in arg_inds: + xs.append([args[x]]) + + res = self.evaluate_func(beta_func, mtx, xs) + return res + else: + if exp: + def pybamm_function(*args): + xs = [] + for x in arg_inds: + if div_const: + xs.append([args[x] / div_const[0]]) + else: + xs.append([args[x]]) + + res = np.exp(self.evaluate_func(beta_func, mtx, xs)) + return res + else: + def pybamm_function(*args): + xs = [] + for x in arg_inds: + xs.append([args[x]]) + + res = self.evaluate_func(beta_func, mtx, xs) + return res + if type(self.parameter_values[name]) is not float: + function_args = self.parameter_values[name].__code__.co_varnames + function_args_mod = [] + if div_arg is not None: + for x in div_arg: + function_args_mod.append(function_args[x[0]] + str('/') + function_args[x[1]]) + for x in arg_inds: + function_args_mod.append(function_args[x]) + print(f"GP function created for {name} \n inputs are {function_args_mod}") + self.pybamm_function = pybamm_function + return pybamm_function + + + def create_beta_inputs(self,len_mtx, GP): + """ + Generates Input Variables + """ + betas_symbolic = [] + beta_parameters = {} + for i in range(len_mtx): + key_str = GP['Name'] + ' Beta ' + str(i) + betas_symbolic.append(pybamm.InputParameter(key_str)) + if i == 0: + beta_parameters[key_str] = pybop.Parameter( + distribution=pybop.Gaussian(GP['Constant mean'], GP['Constant standard deviation']), + ) + + else: + beta_parameters[key_str] = pybop.Parameter( + distribution=pybop.Gaussian(GP['Bi mean'], GP['Bi standard deviation']), + ) + + self.beta_parameters = beta_parameters + return betas_symbolic, beta_parameters + + def add_to_params(self, GP, damtxs): + + betas_function, beta_parameters = self.create_beta_inputs(len(damtxs) + 1, GP) + self.parameter_values.update(beta_parameters) + pybamm_function = self.add_function(GP['Name'], damtxs, GP['arg_inds'], betas_function, exp=GP['exp'], div_arg=GP['div_arg'], + div_const=GP['div_const']) + self.parameter_values.update({GP['Name']:pybamm_function}) + + def get_parameter_values(self): + return self.parameter_values + + diff --git a/pyproject.toml b/pyproject.toml index 9e7540668..a5515df04 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -29,16 +29,14 @@ classifiers = [ # versions in the tests in nightly_dependency_tests.yml as appropriate requires-python = ">=3.10, <3.15" dependencies = [ - "pybamm[plot]>=26.3.0", + "pybamm>=26.3.0", "numpy>=1.26", "scipy>=1.12", "pints>=0.6.0", - "polars>=1.18.0", ] [project.optional-dependencies] -plotly = ["plotly>=6"] -scienceplots = ["SciencePlots>=2.2"] +plot = ["plotly>=6"] salib = ["SALib>=1.5"] scifem = [ "scikit-fem>=8.1.0" # scikit-fem is a dependency for the multi-dimensional pybamm models @@ -51,9 +49,12 @@ ep-bolfi = [ "ep-bolfi>=3.0.2; python_version < '3.13'" ] pyprobe = [ - "PyProBE-Data>=2.6.0;python_version >= '3.11' and python_version < '3.13'" + "PyProBE-Data>=2.5.0;python_version >= '3.11' and python_version < '3.13'" ] -all = ["pybop[plotly,scienceplots,salib,scifem,bpx,pyprobe,ep-bolfi]"] +fokl = [ + "FoKL>=1.2.0" +] +all = ["pybop[plot,salib,scifem,bpx,pyprobe,ep-bolfi,fokl]"] [dependency-groups] docs = [ @@ -82,10 +83,13 @@ dev = [ exclude = [ "papers/", "tests/", + "benchmarks/", "docs/", "assets/", "scripts/", + "multiprocessing_bench.py", "noxfile.py", + "asv.conf.json", "conftest.py", "CONTRIBUTING.md", ] @@ -96,11 +100,14 @@ exclude = [ exclude = [ "papers/", "tests/", + "benchmarks/", "docs/", "assets/", "scripts/", "examples/", + "multiprocessing_bench.py", "noxfile.py", + "asv.conf.json", "conftest.py", "CONTRIBUTING.md", "CHANGELOG.md", From b5c19491829ef4f93cb63fcc5067db3a264985c1 Mon Sep 17 00:00:00 2001 From: derek-slack Date: Thu, 13 Aug 2026 17:18:35 -0400 Subject: [PATCH 3/9] Documentation --- .../FoKLGPy_Fitting.py | 159 +++++++++++++ .../SPMe_FoKLGPy_Fitting.py | 75 ------- pybop/parameters/gp_parameter.py | 211 +++++++++++++----- 3 files changed, 312 insertions(+), 133 deletions(-) create mode 100644 examples/scripts/battery_parameterisation/FoKLGPy_Fitting.py delete mode 100644 examples/scripts/battery_parameterisation/SPMe_FoKLGPy_Fitting.py diff --git a/examples/scripts/battery_parameterisation/FoKLGPy_Fitting.py b/examples/scripts/battery_parameterisation/FoKLGPy_Fitting.py new file mode 100644 index 000000000..a6bec82a0 --- /dev/null +++ b/examples/scripts/battery_parameterisation/FoKLGPy_Fitting.py @@ -0,0 +1,159 @@ +import numpy as np +import pybamm + +import pybop +import matplotlib.pyplot as plt + +""" +In this example, we demonstrate the use of FoKL decomposed Gaussian Processes for estimation of +some parameters. +""" + + +# Define model and parameter values +model = pybamm.lithium_ion.DFN() + +parameter_values = pybamm.ParameterValues("Chen2020") + +# 1. Define a dynamic pulse experiment to excite spatial concentration gradients + +experiment = pybamm.Experiment( + [ + + "Discharge at 1C for 10 minutes", + "Rest for 5 minutes", + "Discharge at 2C for 3 minutes or until 3.0 V", + "Rest for 10 minutes", + "Charge at 1C for 5 minutes or until 4.2 V", + "Rest for 10 minutes", + + ] * 3 +) + +# 2. Generate a synthetic dataset using the experiment +sim = pybamm.Simulation(model, experiment=experiment, parameter_values=parameter_values) + +solution = sim.solve() + +t_eval = solution["Time [s]"].entries +original = solution["Positive electrode exchange current density [A.m-2]"].entries +sigma = 5e-3 +dataset = pybop.Dataset( + { + "Time [s]": t_eval, + "Current [A]": solution["Current [A]"].entries, + "Voltage [V]": pybop.add_noise(solution["Voltage [V]"].entries, sigma), + "Bulk open-circuit voltage [V]": solution["Bulk open-circuit voltage [V]"].entries, + } +) + +# Create GP terms +counter = 0 +tolerance = 1 +num_of_terms = 7 +bic_i = 1e10 + +GP_options = {'Number of terms':num_of_terms,'arg_inds':[0],'Normalization min-max':{'0':(400,2500)}, + 'Constant mean':3.4,'div_arg':[[1,2]], 'exp':True} + +GP_param_neg = pybop.FoKLGP("Positive electrode exchange-current density [A.m-2]", parameter_values=parameter_values.copy(), options=GP_options, twoway=True) +new_parameters = GP_param_neg.get_parameter_values() + + +simulator = pybop.pybamm.Simulator(model, new_parameters, protocol=dataset) +target = ["Voltage [V]", "Bulk open-circuit voltage [V]"] +cost = pybop.GaussianLogLikelihoodKnownSigma(dataset,sigma=sigma, target=target) +problem = pybop.Problem(simulator, cost) + +# Set up the optimiser +options = pybop.PintsOptions( + max_iterations=100, + max_unchanged_iterations=30, + verbose=True +) +optim = pybop.IRPropPlus(problem, options=options) + + +result = optim.run() +L = result.best_cost +k = len(result.x) +n = len(dataset.data['Time [s]']) +bic = -2*L + k*np.log(n) +print(bic) + + +# Plot the timeseries output +pybop.plot.problem(problem, inputs=result.best_inputs, title="Optimised Comparison") + + +# 1. Re-initialize the simulation using the EXACT SAME experiment used for training +sim_final = pybamm.Simulation( + model, + experiment=experiment, # CRITICAL: Must match the training protocol + parameter_values=new_parameters +) + +# 2. Run the final solve (let the experiment handle the time steps natively) +solution_final = sim_final.solve(inputs=result.best_inputs) + +solution_final.all_inputs = [result.best_inputs] * len(solution_final.all_ys) + +new = solution_final["Positive electrode exchange current density [A.m-2]"].entries +t_max = min(solution["Time [s]"].entries[-1], solution_final["Time [s]"].entries[-1]) + +t_eval_valid = t_eval[t_eval <= t_max] + +j0_original_var = solution["Positive electrode exchange current density [A.m-2]"] +j0_new_var = solution_final["Positive electrode exchange current density [A.m-2]"] + +original_aligned = j0_original_var(t=t_eval_valid) +new_aligned = j0_new_var(t=t_eval_valid) + + +idx_top = 0 +idx_mid = new_aligned.shape[0] // 2 +idx_bot = -1 + +gp_top = new_aligned[idx_top, :] +gp_mid = new_aligned[idx_mid, :] +gp_bot = new_aligned[idx_bot, :] + +ref_top = original_aligned[idx_top, :] +ref_mid = original_aligned[idx_mid, :] +ref_bot = original_aligned[idx_bot, :] + +import plotly.graph_objects as go + +fig = go.Figure() + +# 1. Plot the Analytical Baseline (Chen2020) as dashed lines +fig.add_trace(go.Scatter(x=t_eval_valid, y=ref_top, name='Chen2020 - Near Separator', + line=dict(color='blue', dash='dash'), opacity=0.7)) +fig.add_trace(go.Scatter(x=t_eval_valid, y=ref_mid, name='Chen2020 - Middle', + line=dict(color='orange', dash='dash'), opacity=0.7)) +fig.add_trace(go.Scatter(x=t_eval_valid, y=ref_bot, name='Chen2020 - Near Collector', + line=dict(color='green', dash='dash'), opacity=0.7)) + +# 2. Plot the GP Estimate as solid lines +fig.add_trace(go.Scatter(x=t_eval_valid, y=gp_top, name='GP - Near Separator', + line=dict(color='blue'))) +fig.add_trace(go.Scatter(x=t_eval_valid, y=gp_mid, name='GP - Middle', + line=dict(color='orange'))) +fig.add_trace(go.Scatter(x=t_eval_valid, y=gp_bot, name='GP - Near Collector', + line=dict(color='green'))) + +# 3. Formatting for readability +fig.update_layout( + title='Positive Electrode J0 Fitting', + xaxis_title='Time [s]', + yaxis_title='J0 [A.m-2]', + template='plotly_white', + # Place legend outside the plot + legend=dict(x=1.02, y=0.5) +) + +# 4. Use a dotted grid +fig.update_xaxes(showgrid=True, griddash='dot') +fig.update_yaxes(showgrid=True, griddash='dot') + +fig.show() \ No newline at end of file diff --git a/examples/scripts/battery_parameterisation/SPMe_FoKLGPy_Fitting.py b/examples/scripts/battery_parameterisation/SPMe_FoKLGPy_Fitting.py deleted file mode 100644 index 19fb6f5ba..000000000 --- a/examples/scripts/battery_parameterisation/SPMe_FoKLGPy_Fitting.py +++ /dev/null @@ -1,75 +0,0 @@ -import numpy as np -import pybamm - -import pybop -import matplotlib.pyplot as plt -import pandas as pd -from pathlib import Path - -""" -In this example, we demonstrate the use of FoKL decomposed Gaussian Processes for estimation of -some parameters. -""" - - -# Define model and parameter values -model = pybamm.lithium_ion.DFN() - - -parameter_values = pybamm.ParameterValues("Chen2020") - -# Generate a synthetic dataset -sim = pybamm.Simulation(model, parameter_values=parameter_values) -t_eval = np.linspace(0, 3400, 240) -solution = sim.solve(t_eval=t_eval) -original = solution["Positive electrode exchange current density [A.m-2]"].entries - -sigma = 5e-3 -dataset = pybop.Dataset( - { - "Time [s]": t_eval, - "Current [A]": solution["Current [A]"](t_eval), - "Voltage [V]": pybop.add_noise(solution["Voltage [V]"](t_eval), sigma), - "Bulk open-circuit voltage [V]": solution["Bulk open-circuit voltage [V]"]( - t_eval - ), - } -) - -# Create GP terms -GP_options = {'Number of terms':5,'Constant mean':8,'div_arg':[[1,2],[0,2]], 'exp':False} -GP_param_neg = pybop.FoKLGP("Positive electrode exchange-current density [A.m-2]", parameter_values=parameter_values, options=GP_options, twoway=True) -new_parameters = GP_param_neg.get_parameter_values() - -simulator = pybop.pybamm.Simulator(model, new_parameters, protocol=dataset) -target = ["Voltage [V]", "Bulk open-circuit voltage [V]"] -cost = pybop.RootMeanSquaredError(dataset, target=target) -problem = pybop.Problem(simulator, cost) - -# Set up the optimiser -options = pybop.PintsOptions( - max_iterations=300, - max_unchanged_iterations=50, - verbose=True -) -optim = pybop.IRPropPlus(problem, options=options) - -# Run the optimisation -result = optim.run() -print(result) - -# Plot the timeseries output -pybop.plot.problem(problem, inputs=result.best_inputs, title="Optimised Comparison") - -# -# sim = pybamm.Simulation(model, parameter_values=new_parameters) -# -# solution = sim.solve(t_eval=t_eval, inputs=result.best_inputs) -# solution.all_inputs = [result.best_inputs] -# new = solution["Positive electrode exchange current density [A.m-2]"].entries -# -# print(np.sum((new-original)**2)/np.size(original)) - -# Plot the optimisation result -result.plot_parameters() - diff --git a/pybop/parameters/gp_parameter.py b/pybop/parameters/gp_parameter.py index 2f779471f..8ffbd1bf1 100644 --- a/pybop/parameters/gp_parameter.py +++ b/pybop/parameters/gp_parameter.py @@ -1,32 +1,81 @@ -import numpy as np -import pybamm +from typing import Any -import FoKL -from FoKL.getKernels import sp500, bernoulli +import numpy as np import itertools +import pybamm +from pybop.parameters.parameter import Parameter +from pybop.parameters.distributions import Gaussian +from pybop.pybamm.parameter_utils import ParameterValues -import pybop -from pybop.parameters import parameter +try: + import FoKL + from FoKL.getKernels import sp500, bernoulli + FOKL_AVAILABLE = True +except ImportError: + FOKL_AVAILABLE = False class FoKLGP: """ Creates parameter functions as decomposed GPs - options: [dict] """ def __init__(self,name, options=None, parameter_values=None, kernel='Bernoulli Polynomial', twoway=True): + if not FOKL_AVAILABLE: + raise ModuleNotFoundError( + "The `FoKL` package is required to use FoKLGP objects. " + "Please install it using: pip install FoKL" + ) GP_dict_list = self._process_options(name, options, parameter_values) self.evaluate_func = self._evaluate_parameter(kernel) self.damtx = self._create_interaction_matrix(GP_dict_list['Number of terms'], GP_dict_list['Number of inputs'], twoway) self.add_to_params(GP_dict_list,self.damtx) + def __call__(self, *args): return self.pybamm_function(*args) - def _process_options(self, name, options, parameter_values): + def _process_options(self, name: str, options: dict[Any], parameter_values: ParameterValues): """ + Process configuration options for a FoKLGP object and apply defaults. + + Parameters + ---------- + name : str + Name of the parameter being processed. + options : dict or None + User-supplied configuration options to override defaults. Expected + dictionary keys include: + * 'arg_inds' (list[int]): Indices of PyBAMM function arguments. + * 'div_arg' (list[list[int]]): Indices for division terms, e.g., + [[3, 2]] computes input 3 divided by input 2. + * 'div_const' (int): Normalization term to scale inputs between 0 and 1. + * 'inv_arg' (list[int]): Indices to invert, e.g., [3] returns 1 / input 3. + * 'Number of terms' (int): Model order depth (e.g., 3 creates 7 terms). + * 'Constant mean' (float): Mean for the first Beta parameter (B0). + Inferred from `parameter_values` if not provided. + * 'Constant standard deviation' (float): Standard deviation for B0. + * 'Bi mean' (float): Mean for subsequent Beta parameters (Bi). + * 'Bi standard deviation' (float): Standard deviation for Bi. + * 'exp' (bool): If True, applies log-transformation to B0_mean. + parameter_values : dict or Mapping + Dictionary containing base parameter values, used to calculate + 'Constant mean' if it is missing from options. + + Returns + ------- + dict + The finalized options dictionary containing both defaults and + user-defined overrides. + + Raises + ------ + ValueError + If 'Constant mean' is missing and the value in `parameter_values` + cannot be resolved to a constant. + """ - default_options = {'arg_inds':None, 'exp':True, 'Number of inputs':1, 'div_arg':None, 'div_const':None, - 'Number of terms':1, 'Constant standard deviation':0.5,'Bi mean':0, 'Bi standard deviation':0.5} + default_options = {'arg_inds':None, 'exp':True, 'Number of inputs':1, 'div_arg':None, 'div_const':None, 'inv_arg':None, + 'Number of terms':1, 'Constant standard deviation':0.2,'Bi mean':0, 'Bi standard deviation':0.2, + 'Normalization min-max':{}} default_options.update({'Name':name}) if options is not None: default_options.update(options) @@ -41,45 +90,68 @@ def _process_options(self, name, options, parameter_values): default_options.update({'Constant mean': B0_mean}) except: raise ValueError(f'Default parameter value for {name} is not a constant, please supply an estimate') + num_inputs = 0 if default_options['arg_inds'] is not None: num_inputs += len(default_options['arg_inds']) if default_options['div_arg'] is not None: num_inputs += len(default_options['div_arg']) + if default_options['inv_arg'] is not None: + num_inputs += len(default_options['inv_arg']) + + for arg in default_options['arg_inds']: + if str(arg) not in default_options['Normalization min-max']: + default_options['Normalization min-max'] = (0,1) default_options['Number of inputs']=num_inputs self.parameter_values = parameter_values return default_options def _create_interaction_matrix(self, number_of_terms, number_of_inputs, twoway, damtx=[]): + """ + Creates interaction matrix, defines terms of model expansion + """ def perms(x): """Python equivalent of MATLAB perms.""" a = np.array(np.vstack(list(itertools.permutations(x)))[::-1]) return a + def sum_to_n(n, size, limit=None): + """Produce all lists of `size` positive integers in decreasing order + that add up to `n`.""" + if size == 1: + yield [n] + return + if limit is None: + limit = n + start = (n + size - 1) // size + stop = min(limit, n - size + 1) + 1 + for i in range(start, stop): + for tail in sum_to_n(n - i, size - 1, i): + yield [i] + tail + if twoway: sett = 2 else: sett = 1 - + principle = np.zeros((number_of_inputs,)) for ind in range(1, number_of_terms+1): - indvec = np.zeros((number_of_inputs)) - summ = ind - while summ: - for j in range(0, sett): - indvec[j] = indvec[j] + 1 - summ = summ - 1 - if summ == 0: - break - - vecs = np.unique(perms(indvec), axis=0) - - if np.size(damtx) == 0: - damtx = vecs - else: - damtx = np.append(damtx, vecs, axis=0) + indvecs = [i for i in sum_to_n(ind, size=min(number_of_inputs, sett))] + principle[0] = ind + indvecs.append(list(principle)) + for indvec in indvecs: + new_term = perms(indvec) + if np.size(damtx) == 0: + damtx = new_term + else: + if all(new_term[0] == new_term[1]): + damtx = np.vstack([damtx, new_term[0]]) + else: + damtx = np.vstack([damtx, new_term]) + indvec[0]+=1 + damtx = np.array(damtx) print(damtx) return damtx.astype(int) @@ -94,6 +166,11 @@ def _evaluate_parameter(self, kernel): The symbolic evaluation of the decomposed GP. kernel: Kernel function for evaluation, must be `Cubic Splines` or `Bernoulli Polynomial` + + Returns: + --------- + evaluate_pybamm : function + Symbolic GP function for basis function """ self._set_kernel(kernel) @@ -115,7 +192,6 @@ def evaluate_pybamm( phind_temp = inputs[i] * 499 sett = (phind_temp == 0) phind_temp = phind_temp + sett - r = 1 / 499 # interval of when basis function changes (i.e., when next cubic function defines spline) phind.append(phind_temp - 1) A = [1, 2, 3] @@ -158,8 +234,8 @@ def evaluate_pybamm(betas, inputs, coeff=None): """ - Pybamm Function evaluation - betas: indexed from beta list that relates to + Pybamm Function evaluation, creates symbolic + betas: indexed from beta list that scales basis function expansion """ n = 1 num_basis_terms = len(mtx) @@ -194,18 +270,12 @@ def bernoulli_func(phis, num, x): return evaluate_pybamm - def add_function(self, name, mtx, arg_inds, betas_function, exp=False, div_arg=None, div_const=None, - new_children=None): + def add_function(self, name, mtx, arg_inds, + betas_function, norm_bounds, exp=False, + div_arg=None, div_const=None, + ): """ Creates Parameter function specified as a GP object - inputs: - name: str, Name of parameter to be estimated - arg_inds: list of int, Index of inputs to parameter function from PyBaMM - list_index: int, Index of function hyperparameters, betas and mtx, supplied in initialization - exp: Bool, if function should be exponential - div_arg: list of lists of int, optional, Index of arguments to be used as div - ex: div_arg = [[1,4],[3,2]] - inputs to GP would be GP((1/2),(3,2)) """ beta_func = betas_function @@ -217,7 +287,8 @@ def pybamm_function(*args): for x in div_arg: xs.append([args[x[0]] / args[x[1]]]) for x in arg_inds: - xs.append([args[x]]) + xs.append( + [(args[x] - norm_bounds[str(x)][0]) / (norm_bounds[str(x)][1] - norm_bounds[str(x)][0])]) res = np.exp(self.evaluate_func(beta_func, mtx, xs)) return res @@ -227,8 +298,8 @@ def pybamm_function(*args): for x in div_arg: xs.append([args[x[0]] / args[x[1]]]) for x in arg_inds: - xs.append([args[x]]) - + xs.append( + [(args[x] - norm_bounds[str(x)][0]) / (norm_bounds[str(x)][1] - norm_bounds[str(x)][0])]) res = self.evaluate_func(beta_func, mtx, xs) return res else: @@ -236,29 +307,30 @@ def pybamm_function(*args): def pybamm_function(*args): xs = [] for x in arg_inds: - if div_const: - xs.append([args[x] / div_const[0]]) - else: - xs.append([args[x]]) - + xs.append( + [(args[x] - norm_bounds[str(x)][0]) / (norm_bounds[str(x)][1] - norm_bounds[str(x)][0])]) res = np.exp(self.evaluate_func(beta_func, mtx, xs)) return res else: def pybamm_function(*args): xs = [] for x in arg_inds: - xs.append([args[x]]) - + xs.append( + [(args[x] - norm_bounds[str(x)][0]) / (norm_bounds[str(x)][1] - norm_bounds[str(x)][0])]) res = self.evaluate_func(beta_func, mtx, xs) return res + if type(self.parameter_values[name]) is not float: function_args = self.parameter_values[name].__code__.co_varnames function_args_mod = [] if div_arg is not None: for x in div_arg: function_args_mod.append(function_args[x[0]] + str('/') + function_args[x[1]]) - for x in arg_inds: - function_args_mod.append(function_args[x]) + if arg_inds is not None: + for x in arg_inds: + function_args_mod.append(str('(') + function_args[x] + + str(f' - {norm_bounds[str(x)][0]}) / ')+ + str(f'({norm_bounds[str(x)][1] - norm_bounds[str(x)][0]})')) print(f"GP function created for {name} \n inputs are {function_args_mod}") self.pybamm_function = pybamm_function return pybamm_function @@ -266,7 +338,22 @@ def pybamm_function(*args): def create_beta_inputs(self,len_mtx, GP): """ - Generates Input Variables + Generates Input Variables as PyBOP Gaussian Parameters + + Attributes + ------------ + len_mtx : int + Number of basis function expansions + GP : dict + GP dictionary structure + + Returns + -------- + betas_symbolic : list[pybamm.InputParameter] + List of PyBAMM InputParameters added + beta_parameters: dict + Dictionary containing Beta Parameters as Gaussian distributions + """ betas_symbolic = [] beta_parameters = {} @@ -274,27 +361,35 @@ def create_beta_inputs(self,len_mtx, GP): key_str = GP['Name'] + ' Beta ' + str(i) betas_symbolic.append(pybamm.InputParameter(key_str)) if i == 0: - beta_parameters[key_str] = pybop.Parameter( - distribution=pybop.Gaussian(GP['Constant mean'], GP['Constant standard deviation']), + beta_parameters[key_str] = Parameter( + distribution=Gaussian(GP['Constant mean'], GP['Constant standard deviation']), ) else: - beta_parameters[key_str] = pybop.Parameter( - distribution=pybop.Gaussian(GP['Bi mean'], GP['Bi standard deviation']), + beta_parameters[key_str] = Parameter( + distribution=Gaussian(GP['Bi mean'], GP['Bi standard deviation']), ) self.beta_parameters = beta_parameters return betas_symbolic, beta_parameters def add_to_params(self, GP, damtxs): + """ + Updates parameter dictionary with function, input terms + """ betas_function, beta_parameters = self.create_beta_inputs(len(damtxs) + 1, GP) self.parameter_values.update(beta_parameters) - pybamm_function = self.add_function(GP['Name'], damtxs, GP['arg_inds'], betas_function, exp=GP['exp'], div_arg=GP['div_arg'], - div_const=GP['div_const']) + pybamm_function = self.add_function(GP['Name'], damtxs, GP['arg_inds'], betas_function, + GP['Normalization min-max'], exp=GP['exp'], div_arg=GP['div_arg']) self.parameter_values.update({GP['Name']:pybamm_function}) def get_parameter_values(self): + """ + returns parameter values + """ return self.parameter_values + + From 7beb47d1902ab5c76ffa5190d634ef1bce539a20 Mon Sep 17 00:00:00 2001 From: derek-slack Date: Tue, 18 Aug 2026 17:51:00 -0400 Subject: [PATCH 4/9] Update validation script, add documentation --- .../FoKLGPy_Fitting.py | 263 +++++++++++------- pybop/parameters/gp_parameter.py | 39 +-- 2 files changed, 191 insertions(+), 111 deletions(-) diff --git a/examples/scripts/battery_parameterisation/FoKLGPy_Fitting.py b/examples/scripts/battery_parameterisation/FoKLGPy_Fitting.py index a6bec82a0..2631b3586 100644 --- a/examples/scripts/battery_parameterisation/FoKLGPy_Fitting.py +++ b/examples/scripts/battery_parameterisation/FoKLGPy_Fitting.py @@ -1,8 +1,7 @@ import numpy as np import pybamm - import pybop -import matplotlib.pyplot as plt +import plotly.graph_objects as go """ In this example, we demonstrate the use of FoKL decomposed Gaussian Processes for estimation of @@ -17,143 +16,219 @@ # 1. Define a dynamic pulse experiment to excite spatial concentration gradients -experiment = pybamm.Experiment( - [ +experiment = pybamm.Experiment([ + "Discharge at 2C until 2.5 V", +]) - "Discharge at 1C for 10 minutes", - "Rest for 5 minutes", - "Discharge at 2C for 3 minutes or until 3.0 V", - "Rest for 10 minutes", - "Charge at 1C for 5 minutes or until 4.2 V", - "Rest for 10 minutes", - ] * 3 -) +experiment_validate = pybamm.Experiment([ + "Discharge at 3C until 2.7 V", +]) # 2. Generate a synthetic dataset using the experiment sim = pybamm.Simulation(model, experiment=experiment, parameter_values=parameter_values) +sim_validate = pybamm.Simulation(model, experiment=experiment_validate, parameter_values=parameter_values) solution = sim.solve() +solution_validate = sim_validate.solve() +t_eval_validate = solution_validate["Time [s]"].entries t_eval = solution["Time [s]"].entries -original = solution["Positive electrode exchange current density [A.m-2]"].entries + sigma = 5e-3 + +# Construct a PyBOP dataset for training and validation +V_train_data = pybop.add_noise(solution["Voltage [V]"].entries, sigma) dataset = pybop.Dataset( { "Time [s]": t_eval, "Current [A]": solution["Current [A]"].entries, - "Voltage [V]": pybop.add_noise(solution["Voltage [V]"].entries, sigma), + "Voltage [V]": V_train_data, "Bulk open-circuit voltage [V]": solution["Bulk open-circuit voltage [V]"].entries, } ) +V_data = pybop.add_noise(solution_validate["Voltage [V]"].entries, sigma) +dataset_validate = pybop.Dataset( + { + "Time [s]": t_eval_validate, + "Current [A]": solution_validate["Current [A]"].entries, + "Voltage [V]": V_data, + "Bulk open-circuit voltage [V]": solution_validate["Bulk open-circuit voltage [V]"].entries, + } +) + + # Create GP terms -counter = 0 -tolerance = 1 -num_of_terms = 7 -bic_i = 1e10 -GP_options = {'Number of terms':num_of_terms,'arg_inds':[0],'Normalization min-max':{'0':(400,2500)}, - 'Constant mean':3.4,'div_arg':[[1,2]], 'exp':True} +num_of_terms = 2 + +GP_options = {'Number of terms':num_of_terms,'arg_inds':[0],'Normalization min-max':{'0':(-1,4500)}, + 'Constant mean':1.8e-10, 'Constant standard deviation':1e-11,'Bi mean':0,'Bi standard deviation':5e-11,'exp':False} -GP_param_neg = pybop.FoKLGP("Positive electrode exchange-current density [A.m-2]", parameter_values=parameter_values.copy(), options=GP_options, twoway=True) +GP_param_neg = pybop.FoKLGP("Electrolyte diffusivity [m2.s-1]", parameter_values=parameter_values.copy(), options=GP_options, twoway=True) new_parameters = GP_param_neg.get_parameter_values() +# define constant parameter estimation +new_parameters_constant = parameter_values.copy() +new_parameters_constant.update({"Electrolyte diffusivity [m2.s-1]": + pybop.Parameter( + pybop.Gaussian(1.79e-10,2e-11, + truncated_at=[1e-14, 1e-7]) + ) + +}) simulator = pybop.pybamm.Simulator(model, new_parameters, protocol=dataset) target = ["Voltage [V]", "Bulk open-circuit voltage [V]"] cost = pybop.GaussianLogLikelihoodKnownSigma(dataset,sigma=sigma, target=target) problem = pybop.Problem(simulator, cost) +simulator_constant = pybop.pybamm.Simulator(model, new_parameters_constant, protocol=dataset) +problem_constant = pybop.Problem(simulator_constant, cost) + + # Set up the optimiser options = pybop.PintsOptions( - max_iterations=100, - max_unchanged_iterations=30, - verbose=True + max_iterations=1000, + max_unchanged_iterations=100, + verbose=True, ) -optim = pybop.IRPropPlus(problem, options=options) +optim_GP = pybop.XNES(problem, options=options) +result = optim_GP.run() +optim_constant = pybop.XNES(problem_constant, options=options) -result = optim.run() -L = result.best_cost -k = len(result.x) -n = len(dataset.data['Time [s]']) -bic = -2*L + k*np.log(n) -print(bic) +result_constant = optim_constant.run() -# Plot the timeseries output -pybop.plot.problem(problem, inputs=result.best_inputs, title="Optimised Comparison") +new_parameters_constant.update(result_constant.best_inputs) -# 1. Re-initialize the simulation using the EXACT SAME experiment used for training -sim_final = pybamm.Simulation( +sim_final_train_constant = pybamm.Simulation( model, - experiment=experiment, # CRITICAL: Must match the training protocol - parameter_values=new_parameters + experiment=experiment, + parameter_values=new_parameters_constant ) -# 2. Run the final solve (let the experiment handle the time steps natively) -solution_final = sim_final.solve(inputs=result.best_inputs) - -solution_final.all_inputs = [result.best_inputs] * len(solution_final.all_ys) - -new = solution_final["Positive electrode exchange current density [A.m-2]"].entries -t_max = min(solution["Time [s]"].entries[-1], solution_final["Time [s]"].entries[-1]) - -t_eval_valid = t_eval[t_eval <= t_max] - -j0_original_var = solution["Positive electrode exchange current density [A.m-2]"] -j0_new_var = solution_final["Positive electrode exchange current density [A.m-2]"] - -original_aligned = j0_original_var(t=t_eval_valid) -new_aligned = j0_new_var(t=t_eval_valid) - - -idx_top = 0 -idx_mid = new_aligned.shape[0] // 2 -idx_bot = -1 - -gp_top = new_aligned[idx_top, :] -gp_mid = new_aligned[idx_mid, :] -gp_bot = new_aligned[idx_bot, :] +sim_final_train = pybamm.Simulation( + model, + experiment=experiment, + parameter_values=new_parameters +) -ref_top = original_aligned[idx_top, :] -ref_mid = original_aligned[idx_mid, :] -ref_bot = original_aligned[idx_bot, :] +sim_final_test_constant = pybamm.Simulation( + model, + experiment=experiment_validate, + parameter_values=new_parameters_constant +) -import plotly.graph_objects as go +sim_final_test = pybamm.Simulation( + model, + experiment=experiment_validate, + parameter_values=new_parameters +) -fig = go.Figure() - -# 1. Plot the Analytical Baseline (Chen2020) as dashed lines -fig.add_trace(go.Scatter(x=t_eval_valid, y=ref_top, name='Chen2020 - Near Separator', - line=dict(color='blue', dash='dash'), opacity=0.7)) -fig.add_trace(go.Scatter(x=t_eval_valid, y=ref_mid, name='Chen2020 - Middle', - line=dict(color='orange', dash='dash'), opacity=0.7)) -fig.add_trace(go.Scatter(x=t_eval_valid, y=ref_bot, name='Chen2020 - Near Collector', - line=dict(color='green', dash='dash'), opacity=0.7)) - -# 2. Plot the GP Estimate as solid lines -fig.add_trace(go.Scatter(x=t_eval_valid, y=gp_top, name='GP - Near Separator', - line=dict(color='blue'))) -fig.add_trace(go.Scatter(x=t_eval_valid, y=gp_mid, name='GP - Middle', - line=dict(color='orange'))) -fig.add_trace(go.Scatter(x=t_eval_valid, y=gp_bot, name='GP - Near Collector', - line=dict(color='green'))) - -# 3. Formatting for readability -fig.update_layout( - title='Positive Electrode J0 Fitting', +# Run the final solve +solution_final_train = sim_final_train.solve(inputs=result.best_inputs) +solution_final_train_constant = sim_final_train_constant.solve() + +solution_final_test = sim_final_test.solve(inputs=result.best_inputs) +solution_final_test_constant = sim_final_test_constant.solve() + +V_test = solution_final_test['Voltage [V]'].entries +V_train = solution_final_train['Voltage [V]'].entries +t_eval_validate_GP = solution_final_test['Time [s]'].entries +t_eval_validate_GP_train = solution_final_train['Time [s]'].entries + +V_test_constant = solution_final_test_constant['Voltage [V]'].entries +V_train_constant = solution_final_train_constant['Voltage [V]'].entries +t_eval_validate_constant = solution_final_test_constant['Time [s]'].entries +t_eval_validate_constant_train = solution_final_train_constant['Time [s]'].entries + +c_e_train = solution['Electrolyte concentration [mol.m-3]'].entries +c_e_test = solution_validate['Electrolyte concentration [mol.m-3]'].entries +c_e_train_GP = solution_final_train['Electrolyte concentration [mol.m-3]'].entries +c_e_test_GP = solution_final_test['Electrolyte concentration [mol.m-3]'].entries + + +for i in range(4): + fig = go.Figure() + fig.add_trace(go.Scatter( + x=t_eval, y=c_e_train[i * 19, :], + mode='lines', name='Training' + )) + fig.add_trace(go.Scatter( + x=t_eval_validate, y=c_e_test[i * 19, :], + mode='lines', name='Testing' + )) + fig.add_trace(go.Scatter( + x=t_eval_validate_GP_train, y=c_e_train_GP[i * 19, :], + mode='lines', name='GP train' + )) + fig.add_trace(go.Scatter( + x=t_eval_validate_GP, y=c_e_test_GP[i * 19, :], + mode='lines', name='GP test' + )) + fig.update_layout( + xaxis_title='Time [s]', + yaxis_title='Concentration [mol.m-3]', + title=f'Concentration profile (index {i * 19})', + legend=dict(x=0.01, y=0.99) + ) + fig.show() + +# --- Validation experiment (3C discharge) --- +fig_val = go.Figure() +fig_val.add_trace(go.Scatter( + x=t_eval_validate, y=V_data, + mode='markers', name='Data', + marker=dict(size=4) +)) +fig_val.add_trace(go.Scatter( + x=t_eval_validate_GP, y=V_test, + mode='lines', name='FoKL GP', + line=dict(color='green') +)) +fig_val.add_trace(go.Scatter( + x=t_eval_validate_constant, y=V_test_constant, + mode='lines', name='PyBOP Constant', + line=dict(color='red') +)) +fig_val.update_layout( + xaxis_title='Time [s]', + yaxis_title='Voltage [V]', + title='Validation experiment (3C discharge)', + legend=dict(x=0.01, y=0.01) +) +fig_val.show() + +# --- Training experiment (2C discharge) --- +fig_train = go.Figure() +fig_train.add_trace(go.Scatter( + x=t_eval, y=V_train_data, + mode='markers', name='Data', + marker=dict(size=4) +)) +fig_train.add_trace(go.Scatter( + x=t_eval_validate_GP_train, y=V_train, + mode='lines', name='FoKL GP', + line=dict(color='green') +)) +fig_train.add_trace(go.Scatter( + x=t_eval_validate_constant_train, y=V_train_constant, + mode='lines', name='PyBOP Constant', + line=dict(color='red') +)) +fig_train.update_layout( xaxis_title='Time [s]', - yaxis_title='J0 [A.m-2]', - template='plotly_white', - # Place legend outside the plot - legend=dict(x=1.02, y=0.5) + yaxis_title='Voltage [V]', + title='Training experiment (2C discharge)', + legend=dict(x=0.01, y=0.01) ) +fig_train.show() -# 4. Use a dotted grid -fig.update_xaxes(showgrid=True, griddash='dot') -fig.update_yaxes(showgrid=True, griddash='dot') +# Plot the optimisation result +result.plot_convergence() +result.plot_parameters() -fig.show() \ No newline at end of file diff --git a/pybop/parameters/gp_parameter.py b/pybop/parameters/gp_parameter.py index 8ffbd1bf1..f2e7fdafe 100644 --- a/pybop/parameters/gp_parameter.py +++ b/pybop/parameters/gp_parameter.py @@ -56,6 +56,9 @@ def _process_options(self, name: str, options: dict[Any], parameter_values: Para * 'Bi mean' (float): Mean for subsequent Beta parameters (Bi). * 'Bi standard deviation' (float): Standard deviation for Bi. * 'exp' (bool): If True, applies log-transformation to B0_mean. + * 'Normalization min-max' (dict[str(int):tuple]) : Normalization minimum and maximum for `arg_ind` terms + (e.g, {'0':(10,20)} results in argument 0 being + normalized between 10 - 20. parameter_values : dict or Mapping Dictionary containing base parameter values, used to calculate 'Constant mean' if it is missing from options. @@ -134,23 +137,25 @@ def sum_to_n(n, size, limit=None): sett = 2 else: sett = 1 - - principle = np.zeros((number_of_inputs,)) - for ind in range(1, number_of_terms+1): - indvecs = [i for i in sum_to_n(ind, size=min(number_of_inputs, sett))] - principle[0] = ind - indvecs.append(list(principle)) - for indvec in indvecs: - new_term = perms(indvec) - if np.size(damtx) == 0: - damtx = new_term - else: - if all(new_term[0] == new_term[1]): - damtx = np.vstack([damtx, new_term[0]]) + if number_of_inputs == 1: + damtx = np.linspace(1,number_of_terms, number_of_terms).astype(int).reshape(-1,1) + else: + principle = np.zeros((number_of_inputs,)) + for ind in range(1, number_of_terms+1): + indvecs = [i for i in sum_to_n(ind, size=min(number_of_inputs, sett))] + principle[0] = ind + indvecs.append(list(principle)) + for indvec in indvecs: + new_term = perms(indvec) + if np.size(damtx) == 0: + damtx = new_term else: - damtx = np.vstack([damtx, new_term]) + if all(new_term[0] == new_term[1]): + damtx = np.vstack([damtx, new_term[0]]) + else: + damtx = np.vstack([damtx, new_term]) - indvec[0]+=1 + indvec[0]+=1 damtx = np.array(damtx) print(damtx) return damtx.astype(int) @@ -362,12 +367,12 @@ def create_beta_inputs(self,len_mtx, GP): betas_symbolic.append(pybamm.InputParameter(key_str)) if i == 0: beta_parameters[key_str] = Parameter( - distribution=Gaussian(GP['Constant mean'], GP['Constant standard deviation']), + Gaussian(GP['Constant mean'], GP['Constant standard deviation']), ) else: beta_parameters[key_str] = Parameter( - distribution=Gaussian(GP['Bi mean'], GP['Bi standard deviation']), + Gaussian(GP['Bi mean'], GP['Bi standard deviation']), ) self.beta_parameters = beta_parameters From 06f0656687a220ce71cd123ed49df38f0373a4af Mon Sep 17 00:00:00 2001 From: derek-slack Date: Tue, 25 Aug 2026 10:40:32 -0400 Subject: [PATCH 5/9] Preparation for PR --- .../FoKLGPy_Fitting.py | 257 ++++++++------- pybop/parameters/gp_parameter.py | 302 +++++++++++------- 2 files changed, 329 insertions(+), 230 deletions(-) diff --git a/examples/scripts/battery_parameterisation/FoKLGPy_Fitting.py b/examples/scripts/battery_parameterisation/FoKLGPy_Fitting.py index 2631b3586..246012436 100644 --- a/examples/scripts/battery_parameterisation/FoKLGPy_Fitting.py +++ b/examples/scripts/battery_parameterisation/FoKLGPy_Fitting.py @@ -1,10 +1,10 @@ -import numpy as np +import plotly.graph_objects as go import pybamm + import pybop -import plotly.graph_objects as go """ -In this example, we demonstrate the use of FoKL decomposed Gaussian Processes for estimation of +In this example, we demonstrate the use of FoKL decomposed Gaussian Processes for estimation of some parameters. """ @@ -16,18 +16,26 @@ # 1. Define a dynamic pulse experiment to excite spatial concentration gradients -experiment = pybamm.Experiment([ - "Discharge at 2C until 2.5 V", -]) +experiment = pybamm.Experiment( + [ + "Discharge at 2C until 2.5 V", + ] +) -experiment_validate = pybamm.Experiment([ - "Discharge at 3C until 2.7 V", -]) +experiment_validate = pybamm.Experiment( + [ + "Discharge at 3C until 2.7 V", + "Charge at 1C until 4.0 V", + "Discharge at 1C until 2.7 V", + ] +) # 2. Generate a synthetic dataset using the experiment sim = pybamm.Simulation(model, experiment=experiment, parameter_values=parameter_values) -sim_validate = pybamm.Simulation(model, experiment=experiment_validate, parameter_values=parameter_values) +sim_validate = pybamm.Simulation( + model, experiment=experiment_validate, parameter_values=parameter_values +) solution = sim.solve() solution_validate = sim_validate.solve() @@ -44,7 +52,9 @@ "Time [s]": t_eval, "Current [A]": solution["Current [A]"].entries, "Voltage [V]": V_train_data, - "Bulk open-circuit voltage [V]": solution["Bulk open-circuit voltage [V]"].entries, + "Bulk open-circuit voltage [V]": solution[ + "Bulk open-circuit voltage [V]" + ].entries, } ) @@ -54,7 +64,9 @@ "Time [s]": t_eval_validate, "Current [A]": solution_validate["Current [A]"].entries, "Voltage [V]": V_data, - "Bulk open-circuit voltage [V]": solution_validate["Bulk open-circuit voltage [V]"].entries, + "Bulk open-circuit voltage [V]": solution_validate[ + "Bulk open-circuit voltage [V]" + ].entries, } ) @@ -63,28 +75,45 @@ num_of_terms = 2 -GP_options = {'Number of terms':num_of_terms,'arg_inds':[0],'Normalization min-max':{'0':(-1,4500)}, - 'Constant mean':1.8e-10, 'Constant standard deviation':1e-11,'Bi mean':0,'Bi standard deviation':5e-11,'exp':False} - -GP_param_neg = pybop.FoKLGP("Electrolyte diffusivity [m2.s-1]", parameter_values=parameter_values.copy(), options=GP_options, twoway=True) +GP_options = { + "Number of terms": num_of_terms, + "arg_inds": [0], # Argument [0] corresponds to concentration in the electrolyte + "Normalization min-max": { + "0": (-1, 4500) + }, # Normalizing such that inputs are between 0-1 is necessary + "Constant mean": 1.8e-10, # Beta 0 mean + "Constant standard deviation": 1e-11, + "Bi standard deviation": 5e-11, + "exp": False, + "Verbose": True, +} + +GP_param_neg = pybop.FoKLGP( + "Electrolyte diffusivity [m2.s-1]", + parameter_values=parameter_values.copy(), + options=GP_options, + twoway=True, +) new_parameters = GP_param_neg.get_parameter_values() # define constant parameter estimation new_parameters_constant = parameter_values.copy() -new_parameters_constant.update({"Electrolyte diffusivity [m2.s-1]": - pybop.Parameter( - pybop.Gaussian(1.79e-10,2e-11, - truncated_at=[1e-14, 1e-7]) - ) - -}) +new_parameters_constant.update( + { + "Electrolyte diffusivity [m2.s-1]": pybop.Parameter( + pybop.Gaussian(1.79e-10, 2e-11, truncated_at=[1e-14, 1e-7]) + ) + } +) simulator = pybop.pybamm.Simulator(model, new_parameters, protocol=dataset) target = ["Voltage [V]", "Bulk open-circuit voltage [V]"] -cost = pybop.GaussianLogLikelihoodKnownSigma(dataset,sigma=sigma, target=target) +cost = pybop.GaussianLogLikelihoodKnownSigma(dataset, sigma=sigma, target=target) problem = pybop.Problem(simulator, cost) -simulator_constant = pybop.pybamm.Simulator(model, new_parameters_constant, protocol=dataset) +simulator_constant = pybop.pybamm.Simulator( + model, new_parameters_constant, protocol=dataset +) problem_constant = pybop.Problem(simulator_constant, cost) @@ -106,27 +135,19 @@ new_parameters_constant.update(result_constant.best_inputs) sim_final_train_constant = pybamm.Simulation( - model, - experiment=experiment, - parameter_values=new_parameters_constant + model, experiment=experiment, parameter_values=new_parameters_constant ) sim_final_train = pybamm.Simulation( - model, - experiment=experiment, - parameter_values=new_parameters + model, experiment=experiment, parameter_values=new_parameters ) sim_final_test_constant = pybamm.Simulation( - model, - experiment=experiment_validate, - parameter_values=new_parameters_constant + model, experiment=experiment_validate, parameter_values=new_parameters_constant ) sim_final_test = pybamm.Simulation( - model, - experiment=experiment_validate, - parameter_values=new_parameters + model, experiment=experiment_validate, parameter_values=new_parameters ) # Run the final solve @@ -136,99 +157,119 @@ solution_final_test = sim_final_test.solve(inputs=result.best_inputs) solution_final_test_constant = sim_final_test_constant.solve() -V_test = solution_final_test['Voltage [V]'].entries -V_train = solution_final_train['Voltage [V]'].entries -t_eval_validate_GP = solution_final_test['Time [s]'].entries -t_eval_validate_GP_train = solution_final_train['Time [s]'].entries +V_test = solution_final_test["Voltage [V]"].entries +V_train = solution_final_train["Voltage [V]"].entries +t_eval_validate_GP = solution_final_test["Time [s]"].entries +t_eval_validate_GP_train = solution_final_train["Time [s]"].entries -V_test_constant = solution_final_test_constant['Voltage [V]'].entries -V_train_constant = solution_final_train_constant['Voltage [V]'].entries -t_eval_validate_constant = solution_final_test_constant['Time [s]'].entries -t_eval_validate_constant_train = solution_final_train_constant['Time [s]'].entries +V_test_constant = solution_final_test_constant["Voltage [V]"].entries +V_train_constant = solution_final_train_constant["Voltage [V]"].entries +t_eval_validate_constant = solution_final_test_constant["Time [s]"].entries +t_eval_validate_constant_train = solution_final_train_constant["Time [s]"].entries -c_e_train = solution['Electrolyte concentration [mol.m-3]'].entries -c_e_test = solution_validate['Electrolyte concentration [mol.m-3]'].entries -c_e_train_GP = solution_final_train['Electrolyte concentration [mol.m-3]'].entries -c_e_test_GP = solution_final_test['Electrolyte concentration [mol.m-3]'].entries +c_e_train = solution["Electrolyte concentration [mol.m-3]"].entries +c_e_test = solution_validate["Electrolyte concentration [mol.m-3]"].entries +c_e_train_GP = solution_final_train["Electrolyte concentration [mol.m-3]"].entries +c_e_test_GP = solution_final_test["Electrolyte concentration [mol.m-3]"].entries for i in range(4): fig = go.Figure() - fig.add_trace(go.Scatter( - x=t_eval, y=c_e_train[i * 19, :], - mode='lines', name='Training' - )) - fig.add_trace(go.Scatter( - x=t_eval_validate, y=c_e_test[i * 19, :], - mode='lines', name='Testing' - )) - fig.add_trace(go.Scatter( - x=t_eval_validate_GP_train, y=c_e_train_GP[i * 19, :], - mode='lines', name='GP train' - )) - fig.add_trace(go.Scatter( - x=t_eval_validate_GP, y=c_e_test_GP[i * 19, :], - mode='lines', name='GP test' - )) + fig.add_trace( + go.Scatter(x=t_eval, y=c_e_train[i * 19, :], mode="lines", name="Training") + ) + fig.add_trace( + go.Scatter( + x=t_eval_validate, y=c_e_test[i * 19, :], mode="lines", name="Testing" + ) + ) + fig.add_trace( + go.Scatter( + x=t_eval_validate_GP_train, + y=c_e_train_GP[i * 19, :], + mode="lines", + name="GP train", + ) + ) + fig.add_trace( + go.Scatter( + x=t_eval_validate_GP, y=c_e_test_GP[i * 19, :], mode="lines", name="GP test" + ) + ) fig.update_layout( - xaxis_title='Time [s]', - yaxis_title='Concentration [mol.m-3]', - title=f'Concentration profile (index {i * 19})', - legend=dict(x=0.01, y=0.99) + xaxis_title="Time [s]", + yaxis_title="Concentration [mol.m-3]", + title=f"Concentration profile (index {i * 19})", + legend=dict(x=0.01, y=0.99), ) fig.show() # --- Validation experiment (3C discharge) --- fig_val = go.Figure() -fig_val.add_trace(go.Scatter( - x=t_eval_validate, y=V_data, - mode='markers', name='Data', - marker=dict(size=4) -)) -fig_val.add_trace(go.Scatter( - x=t_eval_validate_GP, y=V_test, - mode='lines', name='FoKL GP', - line=dict(color='green') -)) -fig_val.add_trace(go.Scatter( - x=t_eval_validate_constant, y=V_test_constant, - mode='lines', name='PyBOP Constant', - line=dict(color='red') -)) +fig_val.add_trace( + go.Scatter( + x=t_eval_validate, y=V_data, mode="markers", name="Data", marker=dict(size=4) + ) +) +fig_val.add_trace( + go.Scatter( + x=t_eval_validate_GP, + y=V_test, + mode="lines", + name="FoKL GP", + line=dict(color="green"), + ) +) +fig_val.add_trace( + go.Scatter( + x=t_eval_validate_constant, + y=V_test_constant, + mode="lines", + name="PyBOP Constant", + line=dict(color="red"), + ) +) fig_val.update_layout( - xaxis_title='Time [s]', - yaxis_title='Voltage [V]', - title='Validation experiment (3C discharge)', - legend=dict(x=0.01, y=0.01) + xaxis_title="Time [s]", + yaxis_title="Voltage [V]", + title="Validation experiment (3C discharge)", + legend=dict(x=0.01, y=0.01), ) fig_val.show() # --- Training experiment (2C discharge) --- fig_train = go.Figure() -fig_train.add_trace(go.Scatter( - x=t_eval, y=V_train_data, - mode='markers', name='Data', - marker=dict(size=4) -)) -fig_train.add_trace(go.Scatter( - x=t_eval_validate_GP_train, y=V_train, - mode='lines', name='FoKL GP', - line=dict(color='green') -)) -fig_train.add_trace(go.Scatter( - x=t_eval_validate_constant_train, y=V_train_constant, - mode='lines', name='PyBOP Constant', - line=dict(color='red') -)) +fig_train.add_trace( + go.Scatter( + x=t_eval, y=V_train_data, mode="markers", name="Data", marker=dict(size=4) + ) +) +fig_train.add_trace( + go.Scatter( + x=t_eval_validate_GP_train, + y=V_train, + mode="lines", + name="FoKL GP", + line=dict(color="green"), + ) +) +fig_train.add_trace( + go.Scatter( + x=t_eval_validate_constant_train, + y=V_train_constant, + mode="lines", + name="PyBOP Constant", + line=dict(color="red"), + ) +) fig_train.update_layout( - xaxis_title='Time [s]', - yaxis_title='Voltage [V]', - title='Training experiment (2C discharge)', - legend=dict(x=0.01, y=0.01) + xaxis_title="Time [s]", + yaxis_title="Voltage [V]", + title="Training experiment (2C discharge)", + legend=dict(x=0.01, y=0.01), ) fig_train.show() # Plot the optimisation result result.plot_convergence() result.plot_parameters() - diff --git a/pybop/parameters/gp_parameter.py b/pybop/parameters/gp_parameter.py index f2e7fdafe..455cbe6ab 100644 --- a/pybop/parameters/gp_parameter.py +++ b/pybop/parameters/gp_parameter.py @@ -1,24 +1,34 @@ +import itertools from typing import Any import numpy as np -import itertools import pybamm -from pybop.parameters.parameter import Parameter + from pybop.parameters.distributions import Gaussian +from pybop.parameters.parameter import Parameter from pybop.pybamm.parameter_utils import ParameterValues try: - import FoKL - from FoKL.getKernels import sp500, bernoulli + from FoKL.getKernels import bernoulli, sp500 + FOKL_AVAILABLE = True except ImportError: FOKL_AVAILABLE = False + class FoKLGP: """ Creates parameter functions as decomposed GPs """ - def __init__(self,name, options=None, parameter_values=None, kernel='Bernoulli Polynomial', twoway=True): + + def __init__( + self, + name, + options=None, + parameter_values=None, + kernel="Bernoulli Polynomial", + twoway=True, + ): if not FOKL_AVAILABLE: raise ModuleNotFoundError( "The `FoKL` package is required to use FoKLGP objects. " @@ -26,13 +36,17 @@ def __init__(self,name, options=None, parameter_values=None, kernel='Bernoulli P ) GP_dict_list = self._process_options(name, options, parameter_values) self.evaluate_func = self._evaluate_parameter(kernel) - self.damtx = self._create_interaction_matrix(GP_dict_list['Number of terms'], GP_dict_list['Number of inputs'], twoway) - self.add_to_params(GP_dict_list,self.damtx) + self.damtx = self._create_interaction_matrix( + GP_dict_list["Number of terms"], GP_dict_list["Number of inputs"], twoway + ) + self.add_to_params(GP_dict_list, self.damtx) def __call__(self, *args): return self.pybamm_function(*args) - def _process_options(self, name: str, options: dict[Any], parameter_values: ParameterValues): + def _process_options( + self, name: str, options: dict[Any], parameter_values: ParameterValues + ): """ Process configuration options for a FoKLGP object and apply defaults. @@ -59,6 +73,7 @@ def _process_options(self, name: str, options: dict[Any], parameter_values: Para * 'Normalization min-max' (dict[str(int):tuple]) : Normalization minimum and maximum for `arg_ind` terms (e.g, {'0':(10,20)} results in argument 0 being normalized between 10 - 20. + * 'Verbose' (bool): Show debugging print statements parameter_values : dict or Mapping Dictionary containing base parameter values, used to calculate 'Constant mean' if it is missing from options. @@ -76,44 +91,61 @@ def _process_options(self, name: str, options: dict[Any], parameter_values: Para cannot be resolved to a constant. """ - default_options = {'arg_inds':None, 'exp':True, 'Number of inputs':1, 'div_arg':None, 'div_const':None, 'inv_arg':None, - 'Number of terms':1, 'Constant standard deviation':0.2,'Bi mean':0, 'Bi standard deviation':0.2, - 'Normalization min-max':{}} - default_options.update({'Name':name}) + default_options = { + "arg_inds": None, + "exp": True, + "Number of inputs": 1, + "div_arg": None, + "div_const": None, + "inv_arg": None, + "Number of terms": 1, + "Constant standard deviation": 0.2, + "Bi mean": 0, + "Bi standard deviation": 0.2, + "Normalization min-max": {}, + "Verbose": False, + } + default_options.update({"Name": name}) if options is not None: default_options.update(options) - if 'Constant mean' not in default_options: - # If no beta constant term distribution described grab from parameters, if this is a function then user supplied + if "Constant mean" not in default_options: + # If no beta constant term distribution described grab from parameters, + # if this is a function then user supplied try: - if default_options['exp']: + if default_options["exp"]: B0_mean = np.log(parameter_values[name]) else: B0_mean = parameter_values[name] - default_options.update({'Constant mean': B0_mean}) - except: - raise ValueError(f'Default parameter value for {name} is not a constant, please supply an estimate') + default_options.update({"Constant mean": B0_mean}) + except (TypeError, KeyError) as err: + raise ValueError( + f"Default parameter value for {name} is not a constant, please supply an estimate" + ) from err num_inputs = 0 - if default_options['arg_inds'] is not None: - num_inputs += len(default_options['arg_inds']) - if default_options['div_arg'] is not None: - num_inputs += len(default_options['div_arg']) - if default_options['inv_arg'] is not None: - num_inputs += len(default_options['inv_arg']) - - for arg in default_options['arg_inds']: - if str(arg) not in default_options['Normalization min-max']: - default_options['Normalization min-max'] = (0,1) - - default_options['Number of inputs']=num_inputs + if default_options["arg_inds"] is not None: + num_inputs += len(default_options["arg_inds"]) + if default_options["div_arg"] is not None: + num_inputs += len(default_options["div_arg"]) + if default_options["inv_arg"] is not None: + num_inputs += len(default_options["inv_arg"]) + + for arg in default_options["arg_inds"] or []: + if str(arg) not in default_options["Normalization min-max"]: + default_options["Normalization min-max"] = {str(arg): (0, 1)} + self.verbose = default_options["Verbose"] + default_options["Number of inputs"] = num_inputs self.parameter_values = parameter_values return default_options - def _create_interaction_matrix(self, number_of_terms, number_of_inputs, twoway, damtx=[]): + def _create_interaction_matrix( + self, number_of_terms, number_of_inputs, twoway, damtx=None + ): """ Creates interaction matrix, defines terms of model expansion """ + def perms(x): """Python equivalent of MATLAB perms.""" a = np.array(np.vstack(list(itertools.permutations(x)))[::-1]) @@ -138,10 +170,16 @@ def sum_to_n(n, size, limit=None): else: sett = 1 if number_of_inputs == 1: - damtx = np.linspace(1,number_of_terms, number_of_terms).astype(int).reshape(-1,1) + damtx = ( + np.linspace(1, number_of_terms, number_of_terms) + .astype(int) + .reshape(-1, 1) + ) else: principle = np.zeros((number_of_inputs,)) - for ind in range(1, number_of_terms+1): + if damtx is None: + damtx = [] + for ind in range(1, number_of_terms + 1): indvecs = [i for i in sum_to_n(ind, size=min(number_of_inputs, sett))] principle[0] = ind indvecs.append(list(principle)) @@ -155,15 +193,16 @@ def sum_to_n(n, size, limit=None): else: damtx = np.vstack([damtx, new_term]) - indvec[0]+=1 + indvec[0] += 1 damtx = np.array(damtx) - print(damtx) + if self.verbose: + print(damtx) return damtx.astype(int) - def _set_kernel(self, kernel = 'Bernoulli Polynomial'): - if kernel == 'Cubic Splines': + def _set_kernel(self, kernel="Bernoulli Polynomial"): + if kernel == "Cubic Splines": self.phis = sp500() - elif kernel == 'Bernoulli Polynomial': + elif kernel == "Bernoulli Polynomial": self.phis = bernoulli() def _evaluate_parameter(self, kernel): @@ -179,13 +218,9 @@ def _evaluate_parameter(self, kernel): """ self._set_kernel(kernel) - if kernel == 'Cubic Splines': - def evaluate_pybamm( - betas, - mtx, - inputs, - coeff=None): + if kernel == "Cubic Splines": + def evaluate_pybamm(betas, mtx, inputs, coeff=None): num_basis_terms = len(mtx) num_inputs = len(mtx[0]) @@ -195,7 +230,7 @@ def evaluate_pybamm( phind = [] for i in range(num_inputs): phind_temp = inputs[i] * 499 - sett = (phind_temp == 0) + sett = phind_temp == 0 phind_temp = phind_temp + sett phind.append(phind_temp - 1) @@ -204,7 +239,7 @@ def evaluate_pybamm( X_sc = [(1 - inputs[0]) ** a for a in A] lspace = [] - for i in range(num_inputs): + for _i in range(num_inputs): lspace.append(np.linspace(0, 499, 499)) lspace = np.array(lspace) @@ -219,11 +254,17 @@ def evaluate_pybamm( coeff = [] for jj in range(4): phispace = self.phis[nid][jj].reshape(1, -1) - phi_interp = pybamm.Interpolant(lspace[0], phispace[0], - phind[k]) + phi_interp = pybamm.Interpolant( + lspace[0], phispace[0], phind[k] + ) coeff.append(phi_interp) - phi *= coeff[0] + coeff[1] * X_sc[0] + coeff[2] * X_sc[1] + coeff[3] * X_sc[2] + phi *= ( + coeff[0] + + coeff[1] * X_sc[0] + + coeff[2] * X_sc[1] + + coeff[3] * X_sc[2] + ) X_sol.append(phi) X_sol_ones = betas[0] @@ -233,16 +274,14 @@ def evaluate_pybamm( mean += X_sol_betas return mean - elif kernel == 'Bernoulli Polynomial': - def evaluate_pybamm(betas, - mtx, - inputs, - coeff=None): + elif kernel == "Bernoulli Polynomial": + + def evaluate_pybamm(betas, mtx, inputs, coeff=None): """ Pybamm Function evaluation, creates symbolic betas: indexed from beta list that scales basis function expansion """ - n = 1 + num_basis_terms = len(mtx) num_inputs = len(mtx[0]) X_sol = [] @@ -252,7 +291,9 @@ def evaluate_pybamm(betas, def bernoulli_func(phis, num, x): if num > 0: coeff = phis[num - 1] - result = coeff[0] + sum(coeff[k] * (x ** k) for k in range(1, len(coeff))) + result = coeff[0] + sum( + coeff[k] * (x**k) for k in range(1, len(coeff)) + ) else: result = 1.0 return result @@ -271,77 +312,91 @@ def bernoulli_func(phis, num, x): mean += X_sol_betas return mean else: - raise NotImplementedError("Kernel must be either `Cubic Splines` or `Bernoulli Polynomial`") + raise NotImplementedError( + "Kernel must be either `Cubic Splines` or `Bernoulli Polynomial`" + ) return evaluate_pybamm - def add_function(self, name, mtx, arg_inds, - betas_function, norm_bounds, exp=False, - div_arg=None, div_const=None, - ): + def add_function( + self, + name, + mtx, + arg_inds, + betas_function, + norm_bounds, + exp=False, + div_arg=None, + ): """ Creates Parameter function specified as a GP object """ beta_func = betas_function - if div_arg: - if exp: - def pybamm_function(*args): - xs = [] - for x in div_arg: - xs.append([args[x[0]] / args[x[1]]]) - for x in arg_inds: - xs.append( - [(args[x] - norm_bounds[str(x)][0]) / (norm_bounds[str(x)][1] - norm_bounds[str(x)][0])]) - - res = np.exp(self.evaluate_func(beta_func, mtx, xs)) - return res - else: - def pybamm_function(*args): - xs = [] - for x in div_arg: - xs.append([args[x[0]] / args[x[1]]]) - for x in arg_inds: - xs.append( - [(args[x] - norm_bounds[str(x)][0]) / (norm_bounds[str(x)][1] - norm_bounds[str(x)][0])]) - res = self.evaluate_func(beta_func, mtx, xs) - return res + if arg_inds is None: + arg_inds = [] + if div_arg is None: + div_arg = [] + + if exp: + + def pybamm_function(*args): + xs = [] + for x in div_arg: + xs.append([args[x[0]] / args[x[1]]]) + + for x in arg_inds: + xs.append( + [ + (args[x] - norm_bounds[str(x)][0]) + / (norm_bounds[str(x)][1] - norm_bounds[str(x)][0]) + ] + ) + + res = np.exp(self.evaluate_func(beta_func, mtx, xs)) + return res else: - if exp: - def pybamm_function(*args): - xs = [] - for x in arg_inds: - xs.append( - [(args[x] - norm_bounds[str(x)][0]) / (norm_bounds[str(x)][1] - norm_bounds[str(x)][0])]) - res = np.exp(self.evaluate_func(beta_func, mtx, xs)) - return res - else: - def pybamm_function(*args): - xs = [] - for x in arg_inds: - xs.append( - [(args[x] - norm_bounds[str(x)][0]) / (norm_bounds[str(x)][1] - norm_bounds[str(x)][0])]) - res = self.evaluate_func(beta_func, mtx, xs) - return res - - if type(self.parameter_values[name]) is not float: - function_args = self.parameter_values[name].__code__.co_varnames - function_args_mod = [] - if div_arg is not None: + + def pybamm_function(*args): + xs = [] + for x in div_arg: + xs.append([args[x[0]] / args[x[1]]]) + for x in arg_inds: + xs.append( + [ + (args[x] - norm_bounds[str(x)][0]) + / (norm_bounds[str(x)][1] - norm_bounds[str(x)][0]) + ] + ) + res = self.evaluate_func(beta_func, mtx, xs) + return res + + if self.verbose: + # Try to pull function string arguments + if type(self.parameter_values[name]) is not float: + function_args = self.parameter_values[name].__code__.co_varnames + function_args_mod = [] + for x in div_arg: - function_args_mod.append(function_args[x[0]] + str('/') + function_args[x[1]]) - if arg_inds is not None: + function_args_mod.append( + function_args[x[0]] + "/" + function_args[x[1]] + ) + for x in arg_inds: - function_args_mod.append(str('(') + function_args[x] + - str(f' - {norm_bounds[str(x)][0]}) / ')+ - str(f'({norm_bounds[str(x)][1] - norm_bounds[str(x)][0]})')) - print(f"GP function created for {name} \n inputs are {function_args_mod}") + function_args_mod.append( + "(" + + function_args[x] + + str(f" - {norm_bounds[str(x)][0]}) / ") + + str(f"({norm_bounds[str(x)][1] - norm_bounds[str(x)][0]})") + ) + print( + f"GP function created for {name} \n inputs are {function_args_mod}" + ) self.pybamm_function = pybamm_function return pybamm_function - - def create_beta_inputs(self,len_mtx, GP): + def create_beta_inputs(self, len_mtx, GP): """ Generates Input Variables as PyBOP Gaussian Parameters @@ -363,16 +418,16 @@ def create_beta_inputs(self,len_mtx, GP): betas_symbolic = [] beta_parameters = {} for i in range(len_mtx): - key_str = GP['Name'] + ' Beta ' + str(i) + key_str = GP["Name"] + " Beta " + str(i) betas_symbolic.append(pybamm.InputParameter(key_str)) if i == 0: beta_parameters[key_str] = Parameter( - Gaussian(GP['Constant mean'], GP['Constant standard deviation']), + Gaussian(GP["Constant mean"], GP["Constant standard deviation"]), ) else: beta_parameters[key_str] = Parameter( - Gaussian(GP['Bi mean'], GP['Bi standard deviation']), + Gaussian(GP["Bi mean"], GP["Bi standard deviation"]), ) self.beta_parameters = beta_parameters @@ -385,16 +440,19 @@ def add_to_params(self, GP, damtxs): betas_function, beta_parameters = self.create_beta_inputs(len(damtxs) + 1, GP) self.parameter_values.update(beta_parameters) - pybamm_function = self.add_function(GP['Name'], damtxs, GP['arg_inds'], betas_function, - GP['Normalization min-max'], exp=GP['exp'], div_arg=GP['div_arg']) - self.parameter_values.update({GP['Name']:pybamm_function}) + pybamm_function = self.add_function( + GP["Name"], + damtxs, + GP["arg_inds"], + betas_function, + GP["Normalization min-max"], + exp=GP["exp"], + div_arg=GP["div_arg"], + ) + self.parameter_values.update({GP["Name"]: pybamm_function}) def get_parameter_values(self): """ returns parameter values """ return self.parameter_values - - - - From b16dc26a8cabed7e912964cff246f546cd187328 Mon Sep 17 00:00:00 2001 From: derek-slack Date: Fri, 28 Aug 2026 10:34:39 -0400 Subject: [PATCH 6/9] Added string calls --- .../FoKLGPy_Fitting.py | 48 +---- pybop/parameters/gp_parameter.py | 186 ++++++++++++------ 2 files changed, 131 insertions(+), 103 deletions(-) diff --git a/examples/scripts/battery_parameterisation/FoKLGPy_Fitting.py b/examples/scripts/battery_parameterisation/FoKLGPy_Fitting.py index 246012436..e37e7f177 100644 --- a/examples/scripts/battery_parameterisation/FoKLGPy_Fitting.py +++ b/examples/scripts/battery_parameterisation/FoKLGPy_Fitting.py @@ -26,8 +26,6 @@ experiment_validate = pybamm.Experiment( [ "Discharge at 3C until 2.7 V", - "Charge at 1C until 4.0 V", - "Discharge at 1C until 2.7 V", ] ) @@ -73,13 +71,15 @@ # Create GP terms -num_of_terms = 2 +num_of_terms = 3 GP_options = { "Number of terms": num_of_terms, - "arg_inds": [0], # Argument [0] corresponds to concentration in the electrolyte + "Arguments": [ + "Electrolyte concentration [mol.m-3]" + ], # Argument [0] corresponds to concentration in the electrolyte "Normalization min-max": { - "0": (-1, 4500) + "Electrolyte concentration [mol.m-3]": (-1, 5000) }, # Normalizing such that inputs are between 0-1 is necessary "Constant mean": 1.8e-10, # Beta 0 mean "Constant standard deviation": 1e-11, @@ -93,6 +93,7 @@ parameter_values=parameter_values.copy(), options=GP_options, twoway=True, + model=model, ) new_parameters = GP_param_neg.get_parameter_values() @@ -167,43 +168,6 @@ t_eval_validate_constant = solution_final_test_constant["Time [s]"].entries t_eval_validate_constant_train = solution_final_train_constant["Time [s]"].entries -c_e_train = solution["Electrolyte concentration [mol.m-3]"].entries -c_e_test = solution_validate["Electrolyte concentration [mol.m-3]"].entries -c_e_train_GP = solution_final_train["Electrolyte concentration [mol.m-3]"].entries -c_e_test_GP = solution_final_test["Electrolyte concentration [mol.m-3]"].entries - - -for i in range(4): - fig = go.Figure() - fig.add_trace( - go.Scatter(x=t_eval, y=c_e_train[i * 19, :], mode="lines", name="Training") - ) - fig.add_trace( - go.Scatter( - x=t_eval_validate, y=c_e_test[i * 19, :], mode="lines", name="Testing" - ) - ) - fig.add_trace( - go.Scatter( - x=t_eval_validate_GP_train, - y=c_e_train_GP[i * 19, :], - mode="lines", - name="GP train", - ) - ) - fig.add_trace( - go.Scatter( - x=t_eval_validate_GP, y=c_e_test_GP[i * 19, :], mode="lines", name="GP test" - ) - ) - fig.update_layout( - xaxis_title="Time [s]", - yaxis_title="Concentration [mol.m-3]", - title=f"Concentration profile (index {i * 19})", - legend=dict(x=0.01, y=0.99), - ) - fig.show() - # --- Validation experiment (3C discharge) --- fig_val = go.Figure() fig_val.add_trace( diff --git a/pybop/parameters/gp_parameter.py b/pybop/parameters/gp_parameter.py index 455cbe6ab..562f9d23c 100644 --- a/pybop/parameters/gp_parameter.py +++ b/pybop/parameters/gp_parameter.py @@ -28,12 +28,25 @@ def __init__( parameter_values=None, kernel="Bernoulli Polynomial", twoway=True, + model=None, ): + """ + Parameters + ---------- + name : str + Name of the parameter being processed + options : dict + options dictionary for GP creation + kernel : str + sets GP kernel, only "Bernoulli Polynomials + + """ if not FOKL_AVAILABLE: raise ModuleNotFoundError( "The `FoKL` package is required to use FoKLGP objects. " "Please install it using: pip install FoKL" ) + self.model = model GP_dict_list = self._process_options(name, options, parameter_values) self.evaluate_func = self._evaluate_parameter(kernel) self.damtx = self._create_interaction_matrix( @@ -92,12 +105,10 @@ def _process_options( """ default_options = { - "arg_inds": None, + "Arguments": None, "exp": True, "Number of inputs": 1, - "div_arg": None, - "div_const": None, - "inv_arg": None, + "Division Arguments": None, "Number of terms": 1, "Constant standard deviation": 0.2, "Bi mean": 0, @@ -124,21 +135,41 @@ def _process_options( ) from err num_inputs = 0 - if default_options["arg_inds"] is not None: - num_inputs += len(default_options["arg_inds"]) - if default_options["div_arg"] is not None: - num_inputs += len(default_options["div_arg"]) - if default_options["inv_arg"] is not None: - num_inputs += len(default_options["inv_arg"]) - - for arg in default_options["arg_inds"] or []: - if str(arg) not in default_options["Normalization min-max"]: - default_options["Normalization min-max"] = {str(arg): (0, 1)} + + default_options["Input names"] = self._check_arguments(default_options) + + if default_options["Arguments"] is not None: + num_inputs += len(default_options["Arguments"]) + if default_options["Division Arguments"] is not None: + num_inputs += len(default_options["Division Arguments"]) + + for arg in default_options["Arguments"] or []: + if arg not in default_options["Normalization min-max"]: + default_options["Normalization min-max"] = {arg: (0, 1)} self.verbose = default_options["Verbose"] default_options["Number of inputs"] = num_inputs self.parameter_values = parameter_values return default_options + def _check_arguments(self, GP_options): + input_names = self._get_function_parameter_input_names( + self.model, GP_options["Name"] + ) + if GP_options["Arguments"] is not None: + for n in GP_options["Arguments"]: + if n not in input_names: + raise ValueError( + f"Input argument {n} not found. Possible inputs are {input_names}" + ) + if GP_options["Division Arguments"] is not None: + for div_arg in GP_options["Division Arguments"]: + for n in div_arg: + if n not in input_names: + raise ValueError( + f"Input argument {n} not found. Possible inputs are {input_names}" + ) + return input_names + def _create_interaction_matrix( self, number_of_terms, number_of_inputs, twoway, damtx=None ): @@ -214,7 +245,7 @@ def _evaluate_parameter(self, kernel): Returns: --------- evaluate_pybamm : function - Symbolic GP function for basis function + Symbolic GP function for defined kernel """ self._set_kernel(kernel) @@ -228,20 +259,16 @@ def evaluate_pybamm(betas, mtx, inputs, coeff=None): mtx = np.array(mtx) phind = [] + X_sc = [] + A = [1, 2, 3] for i in range(num_inputs): - phind_temp = inputs[i] * 499 - sett = phind_temp == 0 + phind_temp = inputs[i][0] * 499 + sett = pybamm.EqualHeaviside(0, phind_temp) phind_temp = phind_temp + sett phind.append(phind_temp - 1) + X_sc.append([(1 - inputs[0][i]) ** a for a in A]) - A = [1, 2, 3] - - X_sc = [(1 - inputs[0]) ** a for a in A] - - lspace = [] - for _i in range(num_inputs): - lspace.append(np.linspace(0, 499, 499)) - lspace = np.array(lspace) + lspace = np.linspace(0, 499, 499) for j in range(num_basis_terms): phi = 1 @@ -255,17 +282,18 @@ def evaluate_pybamm(betas, mtx, inputs, coeff=None): for jj in range(4): phispace = self.phis[nid][jj].reshape(1, -1) phi_interp = pybamm.Interpolant( - lspace[0], phispace[0], phind[k] + lspace, phispace[0], phind[k] ) coeff.append(phi_interp) phi *= ( coeff[0] - + coeff[1] * X_sc[0] - + coeff[2] * X_sc[1] - + coeff[3] * X_sc[2] + + coeff[1] * X_sc[k][0] + + coeff[2] * X_sc[k][1] + + coeff[3] * X_sc[k][2] ) - X_sol.append(phi) + coeff = None + X_sol.append(phi) X_sol_ones = betas[0] mean = X_sol_ones @@ -277,10 +305,6 @@ def evaluate_pybamm(betas, mtx, inputs, coeff=None): elif kernel == "Bernoulli Polynomial": def evaluate_pybamm(betas, mtx, inputs, coeff=None): - """ - Pybamm Function evaluation, creates symbolic - betas: indexed from beta list that scales basis function expansion - """ num_basis_terms = len(mtx) num_inputs = len(mtx[0]) @@ -318,15 +342,54 @@ def bernoulli_func(phis, num, x): return evaluate_pybamm + @staticmethod + def _get_function_parameter_input_names(model, name): + """ + Return the ordered, PyBaMM-standard descriptive input names for the + FunctionParameter called `name`, e.g. + ["Electrolyte concentration [mol.m-3]", "Temperature [K]"]. + + `model` must be the un-built/un-discretised pybamm.BaseModel — once + parameters are processed, FunctionParameter nodes are replaced by + plain Function nodes and this info is gone. + """ + info = model.get_parameter_info() + for var_symbol, _ in info.values(): + if ( + isinstance(var_symbol, pybamm.FunctionParameter) + and var_symbol.name == name + ): + return list(var_symbol.input_names) + raise ValueError( + f"No FunctionParameter named '{name}' found in the model. " + "Call model.print_parameter_info() to see the available names." + ) + + def _unpack_str_inputs(self, arguments, input_names): + arg_inds = [] + for arg in arguments: + pos = input_names.index(arg) + arg_inds.append(pos) + return arg_inds + + def _unpack_div_str_inputs(self, div_args_str, input_names): + div_arg = [] + for term in div_args_str: + num = term[0] + dom = term[1] + div_arg.append([num, dom]) + return div_arg + def add_function( self, name, mtx, - arg_inds, betas_function, - norm_bounds, + input_names, + arguments=None, + division_arguments=None, + norm_bounds=None, exp=False, - div_arg=None, ): """ Creates Parameter function specified as a GP object @@ -334,10 +397,14 @@ def add_function( beta_func = betas_function - if arg_inds is None: + if arguments is None: arg_inds = [] - if div_arg is None: + else: + arg_inds = self._unpack_str_inputs(arguments, input_names) + if division_arguments is None: div_arg = [] + else: + div_arg = self._unpack_div_str_inputs(division_arguments, input_names) if exp: @@ -365,34 +432,30 @@ def pybamm_function(*args): for x in arg_inds: xs.append( [ - (args[x] - norm_bounds[str(x)][0]) - / (norm_bounds[str(x)][1] - norm_bounds[str(x)][0]) + (args[x] - norm_bounds[input_names[x]][0]) + / ( + norm_bounds[input_names[x]][1] + - norm_bounds[input_names[x]][0] + ) ] ) res = self.evaluate_func(beta_func, mtx, xs) return res if self.verbose: - # Try to pull function string arguments - if type(self.parameter_values[name]) is not float: - function_args = self.parameter_values[name].__code__.co_varnames - function_args_mod = [] - - for x in div_arg: - function_args_mod.append( - function_args[x[0]] + "/" + function_args[x[1]] - ) - - for x in arg_inds: + function_args_mod = [] + if division_arguments is not None: + for x in division_arguments: + function_args_mod.append(x[0] + "/" + x[1]) + if arguments is not None: + for x in arguments: function_args_mod.append( "(" - + function_args[x] - + str(f" - {norm_bounds[str(x)][0]}) / ") - + str(f"({norm_bounds[str(x)][1] - norm_bounds[str(x)][0]})") + + x + + str(f" - {norm_bounds[x][0]}) / ") + + str(f"({norm_bounds[x][1] - norm_bounds[x][0]})") ) - print( - f"GP function created for {name} \n inputs are {function_args_mod}" - ) + print(f"GP function created for {name} \n inputs are {function_args_mod}") self.pybamm_function = pybamm_function return pybamm_function @@ -443,11 +506,12 @@ def add_to_params(self, GP, damtxs): pybamm_function = self.add_function( GP["Name"], damtxs, - GP["arg_inds"], betas_function, - GP["Normalization min-max"], + GP["Input names"], + arguments=GP["Arguments"], + division_arguments=GP["Division Arguments"], + norm_bounds=GP["Normalization min-max"], exp=GP["exp"], - div_arg=GP["div_arg"], ) self.parameter_values.update({GP["Name"]: pybamm_function}) From 837355c1db72c4578bec29d180ba91b5b13ef4e9 Mon Sep 17 00:00:00 2001 From: derek-slack Date: Fri, 28 Aug 2026 10:35:51 -0400 Subject: [PATCH 7/9] Reformatting --- pybop/parameters/gp_parameter.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/pybop/parameters/gp_parameter.py b/pybop/parameters/gp_parameter.py index 562f9d23c..190e3fcd7 100644 --- a/pybop/parameters/gp_parameter.py +++ b/pybop/parameters/gp_parameter.py @@ -365,14 +365,16 @@ def _get_function_parameter_input_names(model, name): "Call model.print_parameter_info() to see the available names." ) - def _unpack_str_inputs(self, arguments, input_names): + @staticmethod + def _unpack_str_inputs(arguments, input_names): arg_inds = [] for arg in arguments: pos = input_names.index(arg) arg_inds.append(pos) return arg_inds - def _unpack_div_str_inputs(self, div_args_str, input_names): + @staticmethod + def _unpack_div_str_inputs(div_args_str, input_names): div_arg = [] for term in div_args_str: num = term[0] From 4a4de643fe233ab4aa2b5cff0beb9036b73cbd9a Mon Sep 17 00:00:00 2001 From: derek-slack Date: Fri, 28 Aug 2026 10:43:48 -0400 Subject: [PATCH 8/9] Updating pyproject.toml --- pyproject.toml | 14 +++++--------- 1 file changed, 5 insertions(+), 9 deletions(-) diff --git a/pyproject.toml b/pyproject.toml index a5515df04..6bb2233e3 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -29,14 +29,16 @@ classifiers = [ # versions in the tests in nightly_dependency_tests.yml as appropriate requires-python = ">=3.10, <3.15" dependencies = [ - "pybamm>=26.3.0", + "pybamm[plot]>=26.3.0", "numpy>=1.26", "scipy>=1.12", "pints>=0.6.0", + "polars>=1.18.0", ] [project.optional-dependencies] -plot = ["plotly>=6"] +plotly = ["plotly>=6"] +scienceplots = ["SciencePlots>=2.2"] salib = ["SALib>=1.5"] scifem = [ "scikit-fem>=8.1.0" # scikit-fem is a dependency for the multi-dimensional pybamm models @@ -49,7 +51,7 @@ ep-bolfi = [ "ep-bolfi>=3.0.2; python_version < '3.13'" ] pyprobe = [ - "PyProBE-Data>=2.5.0;python_version >= '3.11' and python_version < '3.13'" + "PyProBE-Data>=2.6.0;python_version >= '3.11' and python_version < '3.13'" ] fokl = [ "FoKL>=1.2.0" @@ -83,13 +85,10 @@ dev = [ exclude = [ "papers/", "tests/", - "benchmarks/", "docs/", "assets/", "scripts/", - "multiprocessing_bench.py", "noxfile.py", - "asv.conf.json", "conftest.py", "CONTRIBUTING.md", ] @@ -100,14 +99,11 @@ exclude = [ exclude = [ "papers/", "tests/", - "benchmarks/", "docs/", "assets/", "scripts/", "examples/", - "multiprocessing_bench.py", "noxfile.py", - "asv.conf.json", "conftest.py", "CONTRIBUTING.md", "CHANGELOG.md", From b8b9163420cd94f7b6c1e750d7494b26d3158454 Mon Sep 17 00:00:00 2001 From: derek-slack <99365545+derek-slack@users.noreply.github.com> Date: Thu, 3 Sep 2026 11:55:00 -0400 Subject: [PATCH 9/9] Update pyproject.toml Co-authored-by: NicolaCourtier <45851982+NicolaCourtier@users.noreply.github.com> --- pyproject.toml | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/pyproject.toml b/pyproject.toml index 6bb2233e3..e7f730cdb 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -53,10 +53,8 @@ ep-bolfi = [ pyprobe = [ "PyProBE-Data>=2.6.0;python_version >= '3.11' and python_version < '3.13'" ] -fokl = [ - "FoKL>=1.2.0" -] -all = ["pybop[plot,salib,scifem,bpx,pyprobe,ep-bolfi,fokl]"] +fokl = ["FoKL>=1.2.0"] +all = ["pybop[plotly,scienceplots,salib,scifem,bpx,pyprobe,ep-bolfi,fokl]"] [dependency-groups] docs = [