Parameter estimation¶
Parameter estimation is given tool when working with mathematical modelling. This page will give some examples of how one can combine the SUND package with common optimizers to perform parameter estimations.
Objective function¶
All optimization tools needs an objective function to evaluate how the parameter values affect the residuals for a given problem. Typically, one formulate the objective function to be minimized. A general example is given below, where we have simulation objects (in the dictionary sims) corresponding to the experiments in the data dictionary. For every experiment, the simulation result is compared towards the data and the residuals is evaluated using a least square method, and normalized using the standard error of the mean.
def fcost(params, sims, data):
cost=0
for experiment in data.keys():
try:
sims[experiment].simulate(time_vector = data[experiment]["time"],
parameter_values = params,
reset = True)
y_sim = sims[experiment].feature_data[:,0]
y = data[experiment]["mean"]
SEM = data[experiment]["SEM"]
cost += np.sum(np.square(((y_sim - y) / SEM)))
except Exception as e:
# print(str(e)) # print exception, not generally useful
cost+=1e29 # in-case of exception add a "high cost"
return cost
Note that a try-except is included, the purpose is to catch faulty simulations (when bad parameters are given to the objective function), use a high cost to discourage the that theses parameters are tested again, and continue the parameter estimation. Without this except clause, the optimization would fail due to CVODE errors which is thrown when the integration fails during simulation.
Scipy¶
In python, Scipy is a well established tool for various numerical operations. By formulating an objective function (like fcost above), we can evaluate the simulations and parameters compared to data. Scipy offers a simple approach to perform optimizations through scipy.optimize. It includes different solvers for various problems - we will offer an example how parameter estimation can be setup up.
First, a method needs to chosen, like differential_evolution. Thereafter, one can call the the optimizer, and provide various arguments - see the documentation, using differential_evolution().
from scipy.optimize import Bounds, differential_evolution
res = differential_evolution(func=fcost, bounds=bounds, args=args, x0=x0,
callback=callback, disp=True) # This is the optimization
func argument and the bounds (upper and lower) of the parameters is defined in the bounds argument - Scipy offers a Bounds structure that can set up the bounds. Thereafter, there is multiple optional, but potentially useful, arguments that can be provided:
- args (tuple) - packages variables that is needed for the objective function
- x0 (array) - parameter values that are used as the initial guess
- callback (function) - allows us to operations between iterations, such as saving the currently best parameters
- disp (bool) - enables/disable progress printouts
Full example¶
Putting this all together, a python implementation is given below. For this to run, one needs the following packages installed.
uv add numpy scipy
pip install numpy scipy
The example is given below and the python file can be downloaded. Additionally you can download the model file and the data file used in this example.
#%% Import packages
import json
import numpy as np
from scipy.optimize import Bounds, differential_evolution
import sund
#%% Load and print the data
with open("extras/data_M1.json", "r") as f:
data = json.load(f)
#%% Install and load the model
sund.install_model('extras/M1.txt')
m1 = sund.load_model("M1")
#%% Create a activity objects
activity1 = sund.Activity(time_unit='m')
activity1.add_output(sund.PIECEWISE_CONSTANT, "A_in", t = data['dose1']['input']['A_in']['t'], f=data['dose1']['input']['A_in']['f'])
activity2 = sund.Activity(time_unit='m')
activity2.add_output(sund.PIECEWISE_CONSTANT, "A_in", t = data['dose2']['input']['A_in']['t'], f=data['dose2']['input']['A_in']['f'])
#%% Create simulation objects
m1_sims = {}
m1_sims['dose1'] = sund.Simulation(models = m1, activities = activity1, time_unit = 'm')
m1_sims['dose2'] = sund.Simulation(models = m1, activities = activity2, time_unit = 'm')
#%% Define a objective function
def fcost(params, sims, data):
cost=0
for experiment in data.keys():
try:
sims[experiment].simulate(time_vector = data[experiment]["time"],
parameter_values = params,
reset = True)
y_sim = sims[experiment].feature_data[:,0]
y = data[experiment]["mean"]
SEM = data[experiment]["SEM"]
cost += np.sum(np.square(((y_sim - y) / SEM)))
except Exception as e:
# print(str(e)) # print exception, not generally useful
cost+=1e30 # in-case of exception add a "high cost"
return cost
# Define a function to use parameters in log-space to calculate the cost
def fcost_log(params_log, sim, data):
params = np.exp(params_log.copy())
return fcost(params, sim, data)
# %% Prepare input for the optimization
# Define a initial parameter set
params0_M1 = [1, 0.0001, 1000000] # = [k1, k2, kfeed]
x0_log_M1 = np.log(params0_M1)
# parameter bounds
bounds_M1_log = Bounds([np.log(1e-6)]*len(x0_log_M1), [np.log(1e6)]*len(x0_log_M1)) # Create a bounds object
# callback function to save the current parameter values and cost
def callback_M1(x,convergence):
with open("./M1-temp.json",'w') as file:
out = {"x": list(x)}
json.dump(out,file)
def callback_M1_log(x,convergence):
callback_M1(np.exp(x.copy()),convergence)
# arguments for the objective function
args_M1 = (m1_sims, data)
# %% call the optimization function
res = differential_evolution(func=fcost_log, bounds=bounds_M1_log, args=args_M1, x0=x0_log_M1, callback=callback_M1_log, disp=True) # This is the optimization
# %% Print and extract the results
res['x'] = np.exp(res['x']) # convert back from log scale
print(f"\nOptimized parameter values: {res['x']}")
print(f"Optimized cost: {res['fun']}")
# save the results to a json file
file_name = f"M1 ({res['fun']:.3f}).json"
with open(file_name,'w') as file:
res_out = {}
res_out['f'] = float(res['fun']) # We save the file as a .json file, and for that we have to convert the cost value that is currently stored as a ndarray into a traditional float
res_out['x'] = list(res['x']) # We save the file as a .json file, and for that we have to convert the parameter values that is currently stored as a ndarray into a traditional python list
json.dump(res_out, file)