Fit of transient photoluminescence (TrPL) with diffusion and recombination model

This example demonstrates how to fit transient photoluminescence (TrPL) using a rate equation model including recombination and diffusion effect using the optimpv package.

The rate equation and model used here are based on the work by Kober-Czerny et al. 2025 see the paper for more details or the GitHub repository. The example data is also taken from the same publication.

[1]:

# Import necessary libraries
import warnings, os, sys, shutil
# remove warnings from the output
os.environ["PYTHONWARNINGS"] = "ignore"
warnings.filterwarnings(action='ignore', category=FutureWarning)
warnings.filterwarnings(action='ignore', category=UserWarning)
import pandas as pd
import matplotlib.pyplot as plt
import numpy as np
from numpy.random import default_rng
import torch, copy, uuid
import pySIMsalabim as sim
from pySIMsalabim.experiments.JV_steady_state import *
import ax, logging
from ax.utils.notebook.plotting import init_notebook_plotting, render
init_notebook_plotting() # for Jupyter notebooks

try:
    from optimpv import *
    from optimpv.RateEqfits.RateEqAgent import RateEqAgent
    from optimpv.RateEqfits.RateEqModel import *
    from optimpv.RateEqfits.Pumps import *
except Exception as e:
    sys.path.append('../') # add the path to the optimpv module
    from optimpv import *
    from optimpv.RateEqfits.RateEqAgent import RateEqAgent
    from optimpv.RateEqfits.RateEqModel import *
    from optimpv.RateEqfits.Pumps import *

[INFO 08-14 11:02:16] ax.utils.notebook.plotting: Injecting Plotly library into cell. Do not overwrite or delete cell.
[INFO 08-14 11:02:16] ax.utils.notebook.plotting: Please see
    (https://ax.dev/tutorials/visualizations.html#Fix-for-plots-that-are-not-rendering)
    if visualizations are not rendering.

Define the parameters for the simulation

[2]:

# Define the parameters to be fitted
# Note: in general for log spaced value it is better to use the foce_log option when doing the Bayesian inference
params = []

k_direct = FitParam(name = 'k_direct', value = 7.5e-17, bounds = [1e-18,1e-15], log_scale = True, rescale = True, value_type = 'float', type='range', display_name=r'$k_{\text{direct}}$', unit='m$^{3}$ s$^{-1}$', axis_type = 'log',force_log=True)
params.append(k_direct)

k_deep = FitParam(name = 'k_deep', value = 1.6e4, bounds = [1e4,1e6], log_scale = True, rescale = True, value_type = 'float', type='range', display_name=r'$k_{\text{deep}}$', unit='s$^{-1}$', axis_type = 'log', force_log=True)
params.append(k_deep)

k_c = FitParam(name = 'k_c', value = 1.3e6, bounds = [1e4,1e8], log_scale = True, rescale = True, value_type = 'float', type='range', display_name=r'$k_{\text{c}}$', unit='s$^{-1}$', axis_type = 'log', force_log=True)
params.append(k_c)

k_e = FitParam(name = 'k_e', value = 7.5e5, bounds = [1e4,1e8], log_scale = True, rescale = True, value_type = 'float', type='range', display_name=r'$k_{\text{e}}$', unit='m$^{3}$ s$^{-1}$', axis_type = 'log',force_log=True)
params.append(k_e)

mu = FitParam(name = 'mu', value = 4e-1*1e-4, bounds = [1e-6,1e-2], log_scale = True, rescale = True, value_type = 'float', type='range', display_name=r'$\mu$', unit='m$^{2}$ V$^{-1}$ s$^{-1}$', axis_type = 'log',force_log=True)
params.append(mu)

Sfront = FitParam(name = 'S_front', value = 7*1e-2, bounds = [1e-4,1e-1], log_scale = True, rescale = True, value_type = 'float', type='range', display_name=r'$S_{\text{front}}$', unit='m s$^{-1}$', axis_type = 'log',force_log=True)
params.append(Sfront)

L = FitParam(name = 'L', value = 600e-9, bounds = [400e-9,1e-6], log_scale = True, rescale = True, value_type = 'float', type='fixed', display_name=r'$L$', unit='m', axis_type = 'linear',force_log=True)
params.append(L)

Sback = FitParam(name = 'S_back', value = 4*1e-2, bounds = [1e-4,1e-1], log_scale = True, rescale = True, value_type = 'float', type='range', display_name=r'$S_{\text{back}}$', unit='m s$^{-1}$', axis_type = 'log',force_log=True)
params.append(Sback)

I_factor_PL = FitParam(name = 'I_factor_PL', value = 1e-28, bounds = [1e-29,1e-27], log_scale = True, rescale = True, value_type = 'float', type='range', display_name=r'$I_{\text{PL}}$', unit='-', axis_type = 'log', force_log=True)
params.append(I_factor_PL)

N_A = FitParam(name = 'N_A', value = 2.1e21, bounds = [1e20,1e23], log_scale = True, rescale = True, value_type = 'float', type='range', display_name=r'$N_A$', unit='m$^{-3}$', axis_type = 'log',force_log=True)
params.append(N_A)

alpha = FitParam(name = 'alpha', value = 3e8, bounds = [1e6,1e8], log_scale = True, rescale = True, value_type = 'float', type='fixed', display_name=r'$\alpha$', unit='m$^{-1}$', axis_type = 'log',)
params.append(alpha)




# original values
params_orig = copy.deepcopy(params)
num_free_params = 0
dum_dic = {}
for i in range(len(params)):
    if params[i].force_log:
        dum_dic[params[i].name] = np.log10(params[i].value)
    else:
        dum_dic[params[i].name] = params[i].value/params[i].fscale
   # we need this just to run the model to generate some fake data

    if params[i].type != 'fixed'    :
        num_free_params += 1

[3]:

# import Data
parent_dir = os.path.abspath('../') # path to the parent directory if not in Notebooks use os.getcwd()
path2data  = os.path.join(parent_dir,'Data','perovskite_trPL')
data_raw = pd.read_csv(os.path.join(path2data,'Seo_FAPI1_Glass_0.dat'), sep=r'\s+',names=['t',"trPL"], skiprows=1)
# convert time to seconds
data_raw['t'] = data_raw['t'] * 1e-9 # convert to seconds
max_idx = data_raw['trPL'].idxmax()
# remove everything before the maximum
data_raw = data_raw.iloc[max_idx-1:]
# reset the index
data_raw = data_raw.dropna() # remove rows with negative trPL values
data_raw = data_raw[data_raw['trPL'] >= 0]# remove rows with negative trPL values
data_raw = data_raw.reset_index(drop=True)

# interpolate the data to have a logarithmically spaced time axis
t_log = np.logspace(np.log10(2e-9), np.log10(data_raw['t'].max()), num=1000)
t_log = np.insert(t_log, 0, 0)  # add 0 to the time array
# interpolate the trPL values
trPL_log = np.interp(t_log, data_raw['t'] - data_raw['t'].min(), data_raw['trPL'])
data2fit = pd.DataFrame({'t': t_log - t_log.min(), 'trPL': trPL_log})

# plot the data
plt.figure()
plt.plot(data_raw['t'] - data_raw['t'].min(), data_raw['trPL'], 'o', label='Seo_FAPI1_Glass_0')
plt.plot(data2fit['t'], data2fit['trPL'], '*', label='Interpolated trPL')
plt.xscale('log')
plt.yscale('log')
plt.xlabel('Time [s]')
plt.ylabel('trPL [a.u.]')
plt.legend()
plt.show()
../_images/examples_TrPL_diff_4_0.png 
[4]:


time = data2fit['t'].values # time in seconds
X = time
y = data2fit['trPL']

# X = data_raw['t'].values - data_raw['t'].min() # time in seconds
# y = data_raw['trPL'] # dummy data

fpu = 50e3 # Frequency of the pump laser in Hz
Fluence = 4.8e15 # Fluence in m-2
z_array = np.linspace(0, L.value, 100) # z-axis in meters
generation = np.exp(-alpha.value * z_array)
generation_sum = np.trapezoid(generation, z_array)
n_0z = Fluence / generation_sum * generation
N0 = np.mean(n_0z) # mean initial carrier density in m-3
background = 0e28 # Background illumination


# Define the Agent and the target metric/loss function
metric = 'nrmse'
loss = 'linear' # 'nrmse' or 'mse' or 'soft_l1' or 'linear'
pump_args = {'N0': N0, 'fpu': fpu , 'background' : background, }
exp_format = 'trPL' # experiment format
RateEq = RateEqAgent(params, [X], [y], model = DBTD_model, pump_model = initial_carrier_density, pump_args = pump_args, fixed_model_args = {}, metric = metric, loss = loss,minimize=True,exp_format=exp_format,detection_limit=0e-5,  compare_type ='log')
y_test = RateEq.run({},exp_format=exp_format)

plt.figure()
plt.plot(X, y_test)
plt.plot(data2fit['t'], data2fit['trPL'], 'o', label='Seo_FAPI1_Glass_0')
plt.plot(X, y, 'o', label='Seo_FAPI1_Glass_0 original')
# plt.xlim(2e-8, 5e-5)
# plt.ylim(1e-21, 1e1)
plt.xscale('log')
plt.yscale('log')
plt.xlabel('Time [s]')
plt.ylabel(exp_format + ' [a.u.]')
plt.show()
../_images/examples_TrPL_diff_5_0.png

Run the optimization

[5]:

from optimpv.axBOtorch.axBOtorchOptimizer import axBOtorchOptimizer
from botorch.acquisition.logei import qLogNoisyExpectedImprovement
from ax.adapter.transforms.standardize_y import StandardizeY
from ax.adapter.transforms.unit_x import UnitX
from ax.adapter.transforms.remove_fixed import RemoveFixed
from ax.adapter.transforms.log import Log
from ax.generators.torch.botorch_modular.utils import ModelConfig
from ax.generators.torch.botorch_modular.surrogate import SurrogateSpec
from gpytorch.kernels import MaternKernel
from gpytorch.kernels import ScaleKernel
from botorch.models import SingleTaskGP

model_gen_kwargs_list = None
parameter_constraints = None

model_kwargs_list = [{},{"torch_device":torch.device("cuda" if torch.cuda.is_available() else "cpu"),'botorch_acqf_class':qLogNoisyExpectedImprovement,'transforms':[RemoveFixed, Log,UnitX, StandardizeY],'surrogate_spec':SurrogateSpec(model_configs=[ModelConfig(botorch_model_class=SingleTaskGP,covar_module_class=ScaleKernel, covar_module_options={'base_kernel':MaternKernel(nu=2.5, ard_num_dims=num_free_params)})])}]

optimizer = axBOtorchOptimizer(params = params, agents = RateEq, models = ['SOBOL','BOTORCH_MODULAR'],n_batches = [1,100], batch_size = [8,2], ax_client = None,  max_parallelism = 100, model_kwargs_list = model_kwargs_list, model_gen_kwargs_list = model_gen_kwargs_list, name = 'ax_opti',parameter_constraints = parameter_constraints,)

[6]:

# res = optimizer.optimize() # run the optimization with ax
optimizer.optimize_turbo(force_continue=True,kwargs_turbo_state={'failure_tolerance':10}) # run the optimization with turbo

[INFO 08-14 11:02:17] optimpv.axBOtorchOptimizer: Starting optimization with 101 batches and a total of 208 trials
[INFO 08-14 11:02:17] optimpv.axBOtorchOptimizer: Starting Sobol batch 1 with 8 trials
[INFO 08-14 11:02:18] optimpv.axBOtorchOptimizer: Finished Sobol
[INFO 08-14 11:02:20] optimpv.axBOtorchOptimizer: Finished Turbo batch 2 with 2 trials with current best value: 6.03e-02, TR length: 8.00e-01
[INFO 08-14 11:02:20] optimpv.axBOtorchOptimizer: Finished Turbo batch 3 with 2 trials with current best value: 4.06e-02, TR length: 8.00e-01
[INFO 08-14 11:02:21] optimpv.axBOtorchOptimizer: Finished Turbo batch 4 with 2 trials with current best value: 4.06e-02, TR length: 8.00e-01
[INFO 08-14 11:02:21] optimpv.axBOtorchOptimizer: Finished Turbo batch 5 with 2 trials with current best value: 4.06e-02, TR length: 8.00e-01
[INFO 08-14 11:02:22] optimpv.axBOtorchOptimizer: Finished Turbo batch 6 with 2 trials with current best value: 4.06e-02, TR length: 8.00e-01
[INFO 08-14 11:02:22] optimpv.axBOtorchOptimizer: Finished Turbo batch 7 with 2 trials with current best value: 4.06e-02, TR length: 8.00e-01
[INFO 08-14 11:02:23] optimpv.axBOtorchOptimizer: Finished Turbo batch 8 with 2 trials with current best value: 4.06e-02, TR length: 8.00e-01
[INFO 08-14 11:02:23] optimpv.axBOtorchOptimizer: Finished Turbo batch 9 with 2 trials with current best value: 4.06e-02, TR length: 8.00e-01
[INFO 08-14 11:02:24] optimpv.axBOtorchOptimizer: Finished Turbo batch 10 with 2 trials with current best value: 4.06e-02, TR length: 8.00e-01
[INFO 08-14 11:02:24] optimpv.axBOtorchOptimizer: Finished Turbo batch 11 with 2 trials with current best value: 4.06e-02, TR length: 8.00e-01
[INFO 08-14 11:02:25] optimpv.axBOtorchOptimizer: Finished Turbo batch 12 with 2 trials with current best value: 4.06e-02, TR length: 8.00e-01
[INFO 08-14 11:02:25] optimpv.axBOtorchOptimizer: Finished Turbo batch 13 with 2 trials with current best value: 4.06e-02, TR length: 4.00e-01
[INFO 08-14 11:02:26] optimpv.axBOtorchOptimizer: Finished Turbo batch 14 with 2 trials with current best value: 3.65e-02, TR length: 4.00e-01
[INFO 08-14 11:02:26] optimpv.axBOtorchOptimizer: Finished Turbo batch 15 with 2 trials with current best value: 3.41e-02, TR length: 4.00e-01
[INFO 08-14 11:02:27] optimpv.axBOtorchOptimizer: Finished Turbo batch 16 with 2 trials with current best value: 3.34e-02, TR length: 8.00e-01
[INFO 08-14 11:02:27] optimpv.axBOtorchOptimizer: Finished Turbo batch 17 with 2 trials with current best value: 3.34e-02, TR length: 8.00e-01
[INFO 08-14 11:02:28] optimpv.axBOtorchOptimizer: Finished Turbo batch 18 with 2 trials with current best value: 3.34e-02, TR length: 8.00e-01
[INFO 08-14 11:02:28] optimpv.axBOtorchOptimizer: Finished Turbo batch 19 with 2 trials with current best value: 3.34e-02, TR length: 8.00e-01
[INFO 08-14 11:02:29] optimpv.axBOtorchOptimizer: Finished Turbo batch 20 with 2 trials with current best value: 3.34e-02, TR length: 8.00e-01
[INFO 08-14 11:02:29] optimpv.axBOtorchOptimizer: Finished Turbo batch 21 with 2 trials with current best value: 3.34e-02, TR length: 8.00e-01
[INFO 08-14 11:02:30] optimpv.axBOtorchOptimizer: Finished Turbo batch 22 with 2 trials with current best value: 3.34e-02, TR length: 8.00e-01
[INFO 08-14 11:02:30] optimpv.axBOtorchOptimizer: Finished Turbo batch 23 with 2 trials with current best value: 3.34e-02, TR length: 8.00e-01
[INFO 08-14 11:02:31] optimpv.axBOtorchOptimizer: Finished Turbo batch 24 with 2 trials with current best value: 3.34e-02, TR length: 8.00e-01
[INFO 08-14 11:02:31] optimpv.axBOtorchOptimizer: Finished Turbo batch 25 with 2 trials with current best value: 3.34e-02, TR length: 8.00e-01
[INFO 08-14 11:02:32] optimpv.axBOtorchOptimizer: Finished Turbo batch 26 with 2 trials with current best value: 3.34e-02, TR length: 4.00e-01
[INFO 08-14 11:02:32] optimpv.axBOtorchOptimizer: Finished Turbo batch 27 with 2 trials with current best value: 3.07e-02, TR length: 4.00e-01
[INFO 08-14 11:02:33] optimpv.axBOtorchOptimizer: Finished Turbo batch 28 with 2 trials with current best value: 3.07e-02, TR length: 4.00e-01
[INFO 08-14 11:02:33] optimpv.axBOtorchOptimizer: Finished Turbo batch 29 with 2 trials with current best value: 3.07e-02, TR length: 4.00e-01
[INFO 08-14 11:02:34] optimpv.axBOtorchOptimizer: Finished Turbo batch 30 with 2 trials with current best value: 3.07e-02, TR length: 4.00e-01
[INFO 08-14 11:02:34] optimpv.axBOtorchOptimizer: Finished Turbo batch 31 with 2 trials with current best value: 3.07e-02, TR length: 4.00e-01
[INFO 08-14 11:02:34] optimpv.axBOtorchOptimizer: Finished Turbo batch 32 with 2 trials with current best value: 3.07e-02, TR length: 4.00e-01
[INFO 08-14 11:02:35] optimpv.axBOtorchOptimizer: Finished Turbo batch 33 with 2 trials with current best value: 3.07e-02, TR length: 4.00e-01
[INFO 08-14 11:02:35] optimpv.axBOtorchOptimizer: Finished Turbo batch 34 with 2 trials with current best value: 2.76e-02, TR length: 4.00e-01
[INFO 08-14 11:02:36] optimpv.axBOtorchOptimizer: Finished Turbo batch 35 with 2 trials with current best value: 2.76e-02, TR length: 4.00e-01
[INFO 08-14 11:02:36] optimpv.axBOtorchOptimizer: Finished Turbo batch 36 with 2 trials with current best value: 2.76e-02, TR length: 4.00e-01
[INFO 08-14 11:02:37] optimpv.axBOtorchOptimizer: Finished Turbo batch 37 with 2 trials with current best value: 2.76e-02, TR length: 4.00e-01
[INFO 08-14 11:02:37] optimpv.axBOtorchOptimizer: Finished Turbo batch 38 with 2 trials with current best value: 2.76e-02, TR length: 4.00e-01
[INFO 08-14 11:02:38] optimpv.axBOtorchOptimizer: Finished Turbo batch 39 with 2 trials with current best value: 2.52e-02, TR length: 4.00e-01
[INFO 08-14 11:02:38] optimpv.axBOtorchOptimizer: Finished Turbo batch 40 with 2 trials with current best value: 2.51e-02, TR length: 4.00e-01
[INFO 08-14 11:02:39] optimpv.axBOtorchOptimizer: Finished Turbo batch 41 with 2 trials with current best value: 2.51e-02, TR length: 4.00e-01
[INFO 08-14 11:02:39] optimpv.axBOtorchOptimizer: Finished Turbo batch 42 with 2 trials with current best value: 2.51e-02, TR length: 4.00e-01
[INFO 08-14 11:02:40] optimpv.axBOtorchOptimizer: Finished Turbo batch 43 with 2 trials with current best value: 2.51e-02, TR length: 4.00e-01
[INFO 08-14 11:02:40] optimpv.axBOtorchOptimizer: Finished Turbo batch 44 with 2 trials with current best value: 2.51e-02, TR length: 4.00e-01
[INFO 08-14 11:02:41] optimpv.axBOtorchOptimizer: Finished Turbo batch 45 with 2 trials with current best value: 2.51e-02, TR length: 4.00e-01
[INFO 08-14 11:02:42] optimpv.axBOtorchOptimizer: Finished Turbo batch 46 with 2 trials with current best value: 2.51e-02, TR length: 4.00e-01
[INFO 08-14 11:02:42] optimpv.axBOtorchOptimizer: Finished Turbo batch 47 with 2 trials with current best value: 2.51e-02, TR length: 4.00e-01
[INFO 08-14 11:02:43] optimpv.axBOtorchOptimizer: Finished Turbo batch 48 with 2 trials with current best value: 2.51e-02, TR length: 4.00e-01
[INFO 08-14 11:02:43] optimpv.axBOtorchOptimizer: Finished Turbo batch 49 with 2 trials with current best value: 2.51e-02, TR length: 4.00e-01
[INFO 08-14 11:02:44] optimpv.axBOtorchOptimizer: Finished Turbo batch 50 with 2 trials with current best value: 2.51e-02, TR length: 2.00e-01
[INFO 08-14 11:02:44] optimpv.axBOtorchOptimizer: Finished Turbo batch 51 with 2 trials with current best value: 2.51e-02, TR length: 2.00e-01
[INFO 08-14 11:02:45] optimpv.axBOtorchOptimizer: Finished Turbo batch 52 with 2 trials with current best value: 2.47e-02, TR length: 2.00e-01
[INFO 08-14 11:02:45] optimpv.axBOtorchOptimizer: Finished Turbo batch 53 with 2 trials with current best value: 2.35e-02, TR length: 2.00e-01
[INFO 08-14 11:02:46] optimpv.axBOtorchOptimizer: Finished Turbo batch 54 with 2 trials with current best value: 2.35e-02, TR length: 2.00e-01
[INFO 08-14 11:02:47] optimpv.axBOtorchOptimizer: Finished Turbo batch 55 with 2 trials with current best value: 2.35e-02, TR length: 2.00e-01
[INFO 08-14 11:02:47] optimpv.axBOtorchOptimizer: Finished Turbo batch 56 with 2 trials with current best value: 2.35e-02, TR length: 2.00e-01
[INFO 08-14 11:02:48] optimpv.axBOtorchOptimizer: Finished Turbo batch 57 with 2 trials with current best value: 2.35e-02, TR length: 2.00e-01
[INFO 08-14 11:02:48] optimpv.axBOtorchOptimizer: Finished Turbo batch 58 with 2 trials with current best value: 2.35e-02, TR length: 2.00e-01
[INFO 08-14 11:02:48] optimpv.axBOtorchOptimizer: Finished Turbo batch 59 with 2 trials with current best value: 2.35e-02, TR length: 2.00e-01
[INFO 08-14 11:02:49] optimpv.axBOtorchOptimizer: Finished Turbo batch 60 with 2 trials with current best value: 2.35e-02, TR length: 2.00e-01
[INFO 08-14 11:02:49] optimpv.axBOtorchOptimizer: Finished Turbo batch 61 with 2 trials with current best value: 2.35e-02, TR length: 2.00e-01
[INFO 08-14 11:02:50] optimpv.axBOtorchOptimizer: Finished Turbo batch 62 with 2 trials with current best value: 2.35e-02, TR length: 2.00e-01
[INFO 08-14 11:02:50] optimpv.axBOtorchOptimizer: Finished Turbo batch 63 with 2 trials with current best value: 2.35e-02, TR length: 1.00e-01
[INFO 08-14 11:02:51] optimpv.axBOtorchOptimizer: Finished Turbo batch 64 with 2 trials with current best value: 2.35e-02, TR length: 1.00e-01
[INFO 08-14 11:02:51] optimpv.axBOtorchOptimizer: Finished Turbo batch 65 with 2 trials with current best value: 2.32e-02, TR length: 1.00e-01
[INFO 08-14 11:02:52] optimpv.axBOtorchOptimizer: Finished Turbo batch 66 with 2 trials with current best value: 2.32e-02, TR length: 1.00e-01
[INFO 08-14 11:02:53] optimpv.axBOtorchOptimizer: Finished Turbo batch 67 with 2 trials with current best value: 2.32e-02, TR length: 1.00e-01
[INFO 08-14 11:02:53] optimpv.axBOtorchOptimizer: Finished Turbo batch 68 with 2 trials with current best value: 2.32e-02, TR length: 1.00e-01
[INFO 08-14 11:02:54] optimpv.axBOtorchOptimizer: Finished Turbo batch 69 with 2 trials with current best value: 2.32e-02, TR length: 1.00e-01
[INFO 08-14 11:02:54] optimpv.axBOtorchOptimizer: Finished Turbo batch 70 with 2 trials with current best value: 2.32e-02, TR length: 1.00e-01
[INFO 08-14 11:02:55] optimpv.axBOtorchOptimizer: Finished Turbo batch 71 with 2 trials with current best value: 2.32e-02, TR length: 1.00e-01
[INFO 08-14 11:02:56] optimpv.axBOtorchOptimizer: Finished Turbo batch 72 with 2 trials with current best value: 2.30e-02, TR length: 1.00e-01
[INFO 08-14 11:02:56] optimpv.axBOtorchOptimizer: Finished Turbo batch 73 with 2 trials with current best value: 2.30e-02, TR length: 1.00e-01
[INFO 08-14 11:02:57] optimpv.axBOtorchOptimizer: Finished Turbo batch 74 with 2 trials with current best value: 2.30e-02, TR length: 1.00e-01
[INFO 08-14 11:02:57] optimpv.axBOtorchOptimizer: Finished Turbo batch 75 with 2 trials with current best value: 2.30e-02, TR length: 1.00e-01
[INFO 08-14 11:02:58] optimpv.axBOtorchOptimizer: Finished Turbo batch 76 with 2 trials with current best value: 2.30e-02, TR length: 1.00e-01
[INFO 08-14 11:02:59] optimpv.axBOtorchOptimizer: Finished Turbo batch 77 with 2 trials with current best value: 2.30e-02, TR length: 1.00e-01
[INFO 08-14 11:02:59] optimpv.axBOtorchOptimizer: Finished Turbo batch 78 with 2 trials with current best value: 2.30e-02, TR length: 1.00e-01
[INFO 08-14 11:03:00] optimpv.axBOtorchOptimizer: Finished Turbo batch 79 with 2 trials with current best value: 2.30e-02, TR length: 1.00e-01
[INFO 08-14 11:03:00] optimpv.axBOtorchOptimizer: Finished Turbo batch 80 with 2 trials with current best value: 2.30e-02, TR length: 1.00e-01
[INFO 08-14 11:03:01] optimpv.axBOtorchOptimizer: Finished Turbo batch 81 with 2 trials with current best value: 2.30e-02, TR length: 1.00e-01
[INFO 08-14 11:03:02] optimpv.axBOtorchOptimizer: Finished Turbo batch 82 with 2 trials with current best value: 2.30e-02, TR length: 5.00e-02
[INFO 08-14 11:03:02] optimpv.axBOtorchOptimizer: Finished Turbo batch 83 with 2 trials with current best value: 2.30e-02, TR length: 5.00e-02
[INFO 08-14 11:03:03] optimpv.axBOtorchOptimizer: Finished Turbo batch 84 with 2 trials with current best value: 2.28e-02, TR length: 5.00e-02
[INFO 08-14 11:03:04] optimpv.axBOtorchOptimizer: Finished Turbo batch 85 with 2 trials with current best value: 2.28e-02, TR length: 5.00e-02
[INFO 08-14 11:03:04] optimpv.axBOtorchOptimizer: Finished Turbo batch 86 with 2 trials with current best value: 2.28e-02, TR length: 5.00e-02
[INFO 08-14 11:03:05] optimpv.axBOtorchOptimizer: Finished Turbo batch 87 with 2 trials with current best value: 2.28e-02, TR length: 5.00e-02
[INFO 08-14 11:03:05] optimpv.axBOtorchOptimizer: Finished Turbo batch 88 with 2 trials with current best value: 2.28e-02, TR length: 5.00e-02
[INFO 08-14 11:03:06] optimpv.axBOtorchOptimizer: Finished Turbo batch 89 with 2 trials with current best value: 2.28e-02, TR length: 5.00e-02
[INFO 08-14 11:03:07] optimpv.axBOtorchOptimizer: Finished Turbo batch 90 with 2 trials with current best value: 2.28e-02, TR length: 5.00e-02
[INFO 08-14 11:03:07] optimpv.axBOtorchOptimizer: Finished Turbo batch 91 with 2 trials with current best value: 2.28e-02, TR length: 5.00e-02
[INFO 08-14 11:03:08] optimpv.axBOtorchOptimizer: Finished Turbo batch 92 with 2 trials with current best value: 2.28e-02, TR length: 5.00e-02
[INFO 08-14 11:03:08] optimpv.axBOtorchOptimizer: Finished Turbo batch 93 with 2 trials with current best value: 2.28e-02, TR length: 5.00e-02
[INFO 08-14 11:03:09] optimpv.axBOtorchOptimizer: Finished Turbo batch 94 with 2 trials with current best value: 2.28e-02, TR length: 2.50e-02
[INFO 08-14 11:03:10] optimpv.axBOtorchOptimizer: Finished Turbo batch 95 with 2 trials with current best value: 2.28e-02, TR length: 2.50e-02
[INFO 08-14 11:03:11] optimpv.axBOtorchOptimizer: Finished Turbo batch 96 with 2 trials with current best value: 2.28e-02, TR length: 2.50e-02
[INFO 08-14 11:03:11] optimpv.axBOtorchOptimizer: Finished Turbo batch 97 with 2 trials with current best value: 2.28e-02, TR length: 2.50e-02
[INFO 08-14 11:03:12] optimpv.axBOtorchOptimizer: Finished Turbo batch 98 with 2 trials with current best value: 2.28e-02, TR length: 2.50e-02
[INFO 08-14 11:03:13] optimpv.axBOtorchOptimizer: Finished Turbo batch 99 with 2 trials with current best value: 2.28e-02, TR length: 2.50e-02
[INFO 08-14 11:03:13] optimpv.axBOtorchOptimizer: Finished Turbo batch 100 with 2 trials with current best value: 2.28e-02, TR length: 2.50e-02
[INFO 08-14 11:03:14] optimpv.axBOtorchOptimizer: Finished Turbo batch 101 with 2 trials with current best value: 2.28e-02, TR length: 2.50e-02
[INFO 08-14 11:03:15] optimpv.axBOtorchOptimizer: Finished Turbo batch 102 with 2 trials with current best value: 2.28e-02, TR length: 2.50e-02
[INFO 08-14 11:03:15] optimpv.axBOtorchOptimizer: Turbo is terminated as the max number (208) of trials is reached

[7]:

# get the best parameters and update the params list in the optimizer and the agent
ax_client = optimizer.ax_client # get the ax client
optimizer.update_params_with_best_balance() # update the params list in the optimizer with the best parameters
RateEq.params = optimizer.params # update the params list in the agent with the best parameters

# print the best parameters
print('Best parameters:')
for p,po in zip(optimizer.params, params_orig):
    # print(p.name, 'fitted value:', p.value, 'original value:', po.value)
    print(f'{p.name} = {p.value:.2e} {p.unit} (original: {po.value:.2e} {po.unit}), bounds: [{p.bounds[0]:.2e}, {p.bounds[1]:.2e}] {p.unit}')

Best parameters:
k_direct = 1.04e-16 m$^{3}$ s$^{-1}$ (original: 7.50e-17 m$^{3}$ s$^{-1}$), bounds: [1.00e-18, 1.00e-15] m$^{3}$ s$^{-1}$
k_deep = 3.03e+04 s$^{-1}$ (original: 1.60e+04 s$^{-1}$), bounds: [1.00e+04, 1.00e+06] s$^{-1}$
k_c = 1.10e+04 s$^{-1}$ (original: 1.30e+06 s$^{-1}$), bounds: [1.00e+04, 1.00e+08] s$^{-1}$
k_e = 1.47e+07 m$^{3}$ s$^{-1}$ (original: 7.50e+05 m$^{3}$ s$^{-1}$), bounds: [1.00e+04, 1.00e+08] m$^{3}$ s$^{-1}$
mu = 1.52e-04 m$^{2}$ V$^{-1}$ s$^{-1}$ (original: 4.00e-05 m$^{2}$ V$^{-1}$ s$^{-1}$), bounds: [1.00e-06, 1.00e-02] m$^{2}$ V$^{-1}$ s$^{-1}$
S_front = 7.31e-03 m s$^{-1}$ (original: 7.00e-02 m s$^{-1}$), bounds: [1.00e-04, 1.00e-01] m s$^{-1}$
L = 6.00e-07 m (original: 6.00e-07 m), bounds: [4.00e-07, 1.00e-06] m
S_back = 7.67e-03 m s$^{-1}$ (original: 4.00e-02 m s$^{-1}$), bounds: [1.00e-04, 1.00e-01] m s$^{-1}$
I_factor_PL = 1.85e-28 - (original: 1.00e-28 -), bounds: [1.00e-29, 1.00e-27] -
N_A = 1.68e+20 m$^{-3}$ (original: 2.10e+21 m$^{-3}$), bounds: [1.00e+20, 1.00e+23] m$^{-3}$
alpha = 3.00e+08 m$^{-1}$ (original: 3.00e+08 m$^{-1}$), bounds: [1.00e+06, 1.00e+08] m$^{-1}$

[8]:

# Plot optimization results
data = ax_client.summarize()
all_metrics = optimizer.all_metrics
plt.figure()
plt.plot(np.minimum.accumulate(data[all_metrics]), label="Best value seen so far")
plt.yscale("log")
plt.xlabel("Iteration")
plt.ylabel("Target metric: " + metric + " with " + loss + " loss")
plt.legend()
plt.title("Best value seen so far")

print("Best value seen so far is ", min(data[all_metrics[0]]), "at iteration ", int(data[all_metrics[0]].idxmin()))
plt.show()

Best value seen so far is  0.02276154755561851 at iteration  208
../_images/examples_TrPL_diff_10_1.png 
[9]:

# rerun the simulation with the best parameters
yfit = RateEq.run(parameters={},exp_format=exp_format) # run the simulation with the best parameters

# viridis = plt.get_cmap('viridis', len(Gfracs))
plt.figure(figsize=(10,10))
linewidth = 2
plt.plot(X, yfit, label='Best fit', color='black', linewidth=linewidth)
plt.plot(X, y, 'o',label='Fitted data', color='blue', linewidth=linewidth)
# for idx, Gfrac in enumerate(Gfracs[::-1]):
#     plt.plot(X[X[:,1]==Gfrac,0],y[X[:,1]==Gfrac],label='Gfrac = '+str(Gfrac),color=viridis(idx),alpha=0.5,linewidth=linewidth)
#     plt.plot(X[X[:,1]==Gfrac,0],yfit[X[:,1]==Gfrac],label='Gfrac = '+str(Gfrac)+' fit',linestyle='--',color=viridis(idx),linewidth=linewidth)
plt.xscale('log')
plt.yscale('log')
plt.xlabel('Time [s]')
plt.ylabel(exp_format + ' [a.u.]')
plt.legend()
plt.show()
../_images/examples_TrPL_diff_11_0.png 
[10]:

from copy import deepcopy
params_scipy = deepcopy(RateEq.params) # make copies of the parameters
RateEq2 = deepcopy(RateEq) # make a copy of the agent
from optimpv.scipyOpti.scipyOptimizer import ScipyOptimizer
scipyOpti = ScipyOptimizer(params=params_scipy, agents=RateEq2, method='L-BFGS-B', options={'maxiter':200}, name='scipy_opti', parallel_agents=True, max_parallelism=os.cpu_count()-1, verbose_logging=True, )

[11]:

scipyOpti.optimize() # run the optimization with scipy

Starting optimization using L-BFGS-B method
Optimization completed with status: CONVERGENCE: RELATIVE REDUCTION OF F <= FACTR*EPSMCH
Final objective value: 0.022720056750400272

[11]:

  message: CONVERGENCE: RELATIVE REDUCTION OF F <= FACTR*EPSMCH
  success: True
   status: 0
      fun: 0.022720056750400272
        x: [-1.600e+01  4.724e+00  4.041e+00  7.167e+00 -3.846e+00
            -1.923e+00 -1.885e+00 -2.773e+01  2.078e+01]
      nit: 25
      jac: [-1.085e-05  9.302e-07 -8.674e-09 -1.398e-07  1.111e-06
            -1.977e-06 -2.375e-06 -4.292e-06  1.070e-06]
     nfev: 320
     njev: 32
 hess_inv: <9x9 LbfgsInvHessProduct with dtype=float64>

[12]:

# get the best parameters and update the params list in the optimizer and the agent
RateEq2.params = scipyOpti.params # update the params list in the agent with the best parameters.append
# print the best parameters
# print the best parameters
print('Best parameters:')
for p,po in zip(scipyOpti.params, params_orig):
    # print(p.name, 'fitted value:', p.value, 'original value:', po.value)
    print(f'{p.name} = {p.value:.2e} {p.unit} (original: {po.value:.2e} {po.unit}), bounds: [{p.bounds[0]:.2e}, {p.bounds[1]:.2e}] {p.unit}')

Best parameters:
k_direct = 9.90e-17 m$^{3}$ s$^{-1}$ (original: 7.50e-17 m$^{3}$ s$^{-1}$), bounds: [1.00e-18, 1.00e-15] m$^{3}$ s$^{-1}$
k_deep = 5.30e+04 s$^{-1}$ (original: 1.60e+04 s$^{-1}$), bounds: [1.00e+04, 1.00e+06] s$^{-1}$
k_c = 1.10e+04 s$^{-1}$ (original: 1.30e+06 s$^{-1}$), bounds: [1.00e+04, 1.00e+08] s$^{-1}$
k_e = 1.47e+07 m$^{3}$ s$^{-1}$ (original: 7.50e+05 m$^{3}$ s$^{-1}$), bounds: [1.00e+04, 1.00e+08] m$^{3}$ s$^{-1}$
mu = 1.42e-04 m$^{2}$ V$^{-1}$ s$^{-1}$ (original: 4.00e-05 m$^{2}$ V$^{-1}$ s$^{-1}$), bounds: [1.00e-06, 1.00e-02] m$^{2}$ V$^{-1}$ s$^{-1}$
S_front = 1.19e-02 m s$^{-1}$ (original: 7.00e-02 m s$^{-1}$), bounds: [1.00e-04, 1.00e-01] m s$^{-1}$
L = 6.00e-07 m (original: 6.00e-07 m), bounds: [4.00e-07, 1.00e-06] m
S_back = 1.30e-02 m s$^{-1}$ (original: 4.00e-02 m s$^{-1}$), bounds: [1.00e-04, 1.00e-01] m s$^{-1}$
I_factor_PL = 1.85e-28 - (original: 1.00e-28 -), bounds: [1.00e-29, 1.00e-27] -
N_A = 6.05e+20 m$^{-3}$ (original: 2.10e+21 m$^{-3}$), bounds: [1.00e+20, 1.00e+23] m$^{-3}$
alpha = 3.00e+08 m$^{-1}$ (original: 3.00e+08 m$^{-1}$), bounds: [1.00e+06, 1.00e+08] m$^{-1}$

[13]:

# rerun the simulation with the best parameters
yfit = RateEq2.run(parameters={},exp_format=exp_format) # run the simulation with the best parameters

# viridis = plt.get_cmap('viridis', len(Gfracs))
plt.figure(figsize=(10,10))
linewidth = 2
plt.plot(X, yfit, label='Best fit', color='black', linewidth=linewidth)
plt.plot(X, y, 'o',label='Fitted data', color='blue', linewidth=linewidth)

plt.xscale('log')
plt.yscale('log')
plt.xlabel('Time [s]')
plt.ylabel(exp_format + ' [a.u.]')
plt.legend()
plt.show()
../_images/examples_TrPL_diff_15_0.png