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, copy
# 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
from copy import deepcopy
# 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 12-05 10:27:53] ax.utils.notebook.plotting: Injecting Plotly library into cell. Do not overwrite or delete cell.
[INFO 12-05 10:27:53] 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]:

# Plot the data to be fitted and the initial guess
tim = data2fit['t'].values # time in seconds
X = tim
y = data2fit['trPL']

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,  transforms = ['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 optimpv.axBOtorch.axUtils import get_VMLC_default_model_kwargs_list

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

[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 12-05 10:27:54] optimpv.axBOtorchOptimizer: Starting optimization with 101 batches and a total of 208 trials
[INFO 12-05 10:27:54] optimpv.axBOtorchOptimizer: Starting Sobol batch 1 with 8 trials
[INFO 12-05 10:27:55] optimpv.axBOtorchOptimizer: Finished Sobol with best value of 0.103077
[INFO 12-05 10:27:57] optimpv.axBOtorchOptimizer: Finished Turbo batch 2 with 2 trials with current best value: 4.22e-02, TR length: 8.00e-01
[INFO 12-05 10:27:57] optimpv.axBOtorchOptimizer: Finished Turbo batch 3 with 2 trials with current best value: 4.22e-02, TR length: 8.00e-01
[INFO 12-05 10:27:58] optimpv.axBOtorchOptimizer: Finished Turbo batch 4 with 2 trials with current best value: 4.22e-02, TR length: 8.00e-01
[INFO 12-05 10:27:58] optimpv.axBOtorchOptimizer: Finished Turbo batch 5 with 2 trials with current best value: 4.22e-02, TR length: 8.00e-01
[INFO 12-05 10:27:59] optimpv.axBOtorchOptimizer: Finished Turbo batch 6 with 2 trials with current best value: 4.22e-02, TR length: 8.00e-01
[INFO 12-05 10:27:59] optimpv.axBOtorchOptimizer: Finished Turbo batch 7 with 2 trials with current best value: 4.22e-02, TR length: 8.00e-01
[INFO 12-05 10:28:00] optimpv.axBOtorchOptimizer: Finished Turbo batch 8 with 2 trials with current best value: 4.22e-02, TR length: 8.00e-01
[INFO 12-05 10:28:00] optimpv.axBOtorchOptimizer: Finished Turbo batch 9 with 2 trials with current best value: 4.22e-02, TR length: 8.00e-01
[INFO 12-05 10:28:01] optimpv.axBOtorchOptimizer: Finished Turbo batch 10 with 2 trials with current best value: 4.22e-02, TR length: 8.00e-01
[INFO 12-05 10:28:01] optimpv.axBOtorchOptimizer: Finished Turbo batch 11 with 2 trials with current best value: 4.22e-02, TR length: 8.00e-01
[INFO 12-05 10:28:01] optimpv.axBOtorchOptimizer: Finished Turbo batch 12 with 2 trials with current best value: 4.22e-02, TR length: 5.00e-01
[INFO 12-05 10:28:02] optimpv.axBOtorchOptimizer: Finished Turbo batch 13 with 2 trials with current best value: 4.22e-02, TR length: 5.00e-01
[INFO 12-05 10:28:02] optimpv.axBOtorchOptimizer: Finished Turbo batch 14 with 2 trials with current best value: 3.94e-02, TR length: 5.00e-01
[INFO 12-05 10:28:03] optimpv.axBOtorchOptimizer: Finished Turbo batch 15 with 2 trials with current best value: 3.19e-02, TR length: 5.00e-01
[INFO 12-05 10:28:03] optimpv.axBOtorchOptimizer: Finished Turbo batch 16 with 2 trials with current best value: 3.19e-02, TR length: 5.00e-01
[INFO 12-05 10:28:04] optimpv.axBOtorchOptimizer: Finished Turbo batch 17 with 2 trials with current best value: 3.19e-02, TR length: 5.00e-01
[INFO 12-05 10:28:04] optimpv.axBOtorchOptimizer: Finished Turbo batch 18 with 2 trials with current best value: 3.19e-02, TR length: 5.00e-01
[INFO 12-05 10:28:05] optimpv.axBOtorchOptimizer: Finished Turbo batch 19 with 2 trials with current best value: 3.19e-02, TR length: 5.00e-01
[INFO 12-05 10:28:05] optimpv.axBOtorchOptimizer: Finished Turbo batch 20 with 2 trials with current best value: 2.58e-02, TR length: 5.00e-01
[INFO 12-05 10:28:06] optimpv.axBOtorchOptimizer: Finished Turbo batch 21 with 2 trials with current best value: 2.58e-02, TR length: 5.00e-01
[INFO 12-05 10:28:06] optimpv.axBOtorchOptimizer: Finished Turbo batch 22 with 2 trials with current best value: 2.58e-02, TR length: 5.00e-01
[INFO 12-05 10:28:07] optimpv.axBOtorchOptimizer: Finished Turbo batch 23 with 2 trials with current best value: 2.58e-02, TR length: 5.00e-01
[INFO 12-05 10:28:07] optimpv.axBOtorchOptimizer: Finished Turbo batch 24 with 2 trials with current best value: 2.58e-02, TR length: 5.00e-01
[INFO 12-05 10:28:08] optimpv.axBOtorchOptimizer: Finished Turbo batch 25 with 2 trials with current best value: 2.58e-02, TR length: 5.00e-01
[INFO 12-05 10:28:08] optimpv.axBOtorchOptimizer: Finished Turbo batch 26 with 2 trials with current best value: 2.58e-02, TR length: 5.00e-01
[INFO 12-05 10:28:09] optimpv.axBOtorchOptimizer: Finished Turbo batch 27 with 2 trials with current best value: 2.58e-02, TR length: 5.00e-01
[INFO 12-05 10:28:09] optimpv.axBOtorchOptimizer: Finished Turbo batch 28 with 2 trials with current best value: 2.58e-02, TR length: 5.00e-01
[INFO 12-05 10:28:10] optimpv.axBOtorchOptimizer: Finished Turbo batch 29 with 2 trials with current best value: 2.58e-02, TR length: 5.00e-01
[INFO 12-05 10:28:10] optimpv.axBOtorchOptimizer: Finished Turbo batch 30 with 2 trials with current best value: 2.58e-02, TR length: 3.12e-01
[INFO 12-05 10:28:11] optimpv.axBOtorchOptimizer: Finished Turbo batch 31 with 2 trials with current best value: 2.58e-02, TR length: 3.12e-01
[INFO 12-05 10:28:12] optimpv.axBOtorchOptimizer: Finished Turbo batch 32 with 2 trials with current best value: 2.58e-02, TR length: 3.12e-01
[INFO 12-05 10:28:12] optimpv.axBOtorchOptimizer: Finished Turbo batch 33 with 2 trials with current best value: 2.58e-02, TR length: 3.12e-01
[INFO 12-05 10:28:13] optimpv.axBOtorchOptimizer: Finished Turbo batch 34 with 2 trials with current best value: 2.34e-02, TR length: 3.12e-01
[INFO 12-05 10:28:13] optimpv.axBOtorchOptimizer: Finished Turbo batch 35 with 2 trials with current best value: 2.34e-02, TR length: 3.12e-01
[INFO 12-05 10:28:14] optimpv.axBOtorchOptimizer: Finished Turbo batch 36 with 2 trials with current best value: 2.34e-02, TR length: 3.12e-01
[INFO 12-05 10:28:14] optimpv.axBOtorchOptimizer: Finished Turbo batch 37 with 2 trials with current best value: 2.33e-02, TR length: 3.12e-01
[INFO 12-05 10:28:15] optimpv.axBOtorchOptimizer: Finished Turbo batch 38 with 2 trials with current best value: 2.33e-02, TR length: 3.12e-01
[INFO 12-05 10:28:15] optimpv.axBOtorchOptimizer: Finished Turbo batch 39 with 2 trials with current best value: 2.33e-02, TR length: 3.12e-01
[INFO 12-05 10:28:16] optimpv.axBOtorchOptimizer: Finished Turbo batch 40 with 2 trials with current best value: 2.33e-02, TR length: 3.12e-01
[INFO 12-05 10:28:16] optimpv.axBOtorchOptimizer: Finished Turbo batch 41 with 2 trials with current best value: 2.33e-02, TR length: 3.12e-01
[INFO 12-05 10:28:17] optimpv.axBOtorchOptimizer: Finished Turbo batch 42 with 2 trials with current best value: 2.33e-02, TR length: 3.12e-01
[INFO 12-05 10:28:17] optimpv.axBOtorchOptimizer: Finished Turbo batch 43 with 2 trials with current best value: 2.33e-02, TR length: 3.12e-01
[INFO 12-05 10:28:18] optimpv.axBOtorchOptimizer: Finished Turbo batch 44 with 2 trials with current best value: 2.33e-02, TR length: 3.12e-01
[INFO 12-05 10:28:18] optimpv.axBOtorchOptimizer: Finished Turbo batch 45 with 2 trials with current best value: 2.33e-02, TR length: 3.12e-01
[INFO 12-05 10:28:19] optimpv.axBOtorchOptimizer: Finished Turbo batch 46 with 2 trials with current best value: 2.33e-02, TR length: 3.12e-01
[INFO 12-05 10:28:19] optimpv.axBOtorchOptimizer: Finished Turbo batch 47 with 2 trials with current best value: 2.33e-02, TR length: 1.95e-01
[INFO 12-05 10:28:20] optimpv.axBOtorchOptimizer: Finished Turbo batch 48 with 2 trials with current best value: 2.27e-02, TR length: 1.95e-01
[INFO 12-05 10:28:20] optimpv.axBOtorchOptimizer: Finished Turbo batch 49 with 2 trials with current best value: 2.26e-02, TR length: 1.95e-01
[INFO 12-05 10:28:21] optimpv.axBOtorchOptimizer: Finished Turbo batch 50 with 2 trials with current best value: 2.26e-02, TR length: 1.95e-01
[INFO 12-05 10:28:21] optimpv.axBOtorchOptimizer: Finished Turbo batch 51 with 2 trials with current best value: 2.26e-02, TR length: 1.95e-01
[INFO 12-05 10:28:22] optimpv.axBOtorchOptimizer: Finished Turbo batch 52 with 2 trials with current best value: 2.26e-02, TR length: 1.95e-01
[INFO 12-05 10:28:22] optimpv.axBOtorchOptimizer: Finished Turbo batch 53 with 2 trials with current best value: 2.26e-02, TR length: 1.95e-01
[INFO 12-05 10:28:23] optimpv.axBOtorchOptimizer: Finished Turbo batch 54 with 2 trials with current best value: 2.26e-02, TR length: 1.95e-01
[INFO 12-05 10:28:23] optimpv.axBOtorchOptimizer: Finished Turbo batch 55 with 2 trials with current best value: 2.26e-02, TR length: 1.95e-01
[INFO 12-05 10:28:24] optimpv.axBOtorchOptimizer: Finished Turbo batch 56 with 2 trials with current best value: 2.26e-02, TR length: 1.95e-01
[INFO 12-05 10:28:25] optimpv.axBOtorchOptimizer: Finished Turbo batch 57 with 2 trials with current best value: 2.26e-02, TR length: 1.95e-01
[INFO 12-05 10:28:25] optimpv.axBOtorchOptimizer: Finished Turbo batch 58 with 2 trials with current best value: 2.26e-02, TR length: 1.95e-01
[INFO 12-05 10:28:26] optimpv.axBOtorchOptimizer: Finished Turbo batch 59 with 2 trials with current best value: 2.26e-02, TR length: 1.22e-01
[INFO 12-05 10:28:26] optimpv.axBOtorchOptimizer: Finished Turbo batch 60 with 2 trials with current best value: 2.25e-02, TR length: 1.22e-01
[INFO 12-05 10:28:27] optimpv.axBOtorchOptimizer: Finished Turbo batch 61 with 2 trials with current best value: 2.24e-02, TR length: 1.22e-01
[INFO 12-05 10:28:27] optimpv.axBOtorchOptimizer: Finished Turbo batch 62 with 2 trials with current best value: 2.24e-02, TR length: 1.22e-01
[INFO 12-05 10:28:28] optimpv.axBOtorchOptimizer: Finished Turbo batch 63 with 2 trials with current best value: 2.24e-02, TR length: 1.22e-01
[INFO 12-05 10:28:28] optimpv.axBOtorchOptimizer: Finished Turbo batch 64 with 2 trials with current best value: 2.24e-02, TR length: 1.22e-01
[INFO 12-05 10:28:29] optimpv.axBOtorchOptimizer: Finished Turbo batch 65 with 2 trials with current best value: 2.24e-02, TR length: 1.22e-01
[INFO 12-05 10:28:29] optimpv.axBOtorchOptimizer: Finished Turbo batch 66 with 2 trials with current best value: 2.24e-02, TR length: 1.22e-01
[INFO 12-05 10:28:30] optimpv.axBOtorchOptimizer: Finished Turbo batch 67 with 2 trials with current best value: 2.20e-02, TR length: 1.22e-01
[INFO 12-05 10:28:30] optimpv.axBOtorchOptimizer: Finished Turbo batch 68 with 2 trials with current best value: 2.20e-02, TR length: 1.22e-01
[INFO 12-05 10:28:31] optimpv.axBOtorchOptimizer: Finished Turbo batch 69 with 2 trials with current best value: 2.20e-02, TR length: 1.22e-01
[INFO 12-05 10:28:31] optimpv.axBOtorchOptimizer: Finished Turbo batch 70 with 2 trials with current best value: 2.20e-02, TR length: 1.22e-01
[INFO 12-05 10:28:32] optimpv.axBOtorchOptimizer: Finished Turbo batch 71 with 2 trials with current best value: 2.19e-02, TR length: 1.22e-01
[INFO 12-05 10:28:32] optimpv.axBOtorchOptimizer: Finished Turbo batch 72 with 2 trials with current best value: 2.19e-02, TR length: 1.22e-01
[INFO 12-05 10:28:33] optimpv.axBOtorchOptimizer: Finished Turbo batch 73 with 2 trials with current best value: 2.19e-02, TR length: 1.22e-01
[INFO 12-05 10:28:34] optimpv.axBOtorchOptimizer: Finished Turbo batch 74 with 2 trials with current best value: 2.19e-02, TR length: 1.22e-01
[INFO 12-05 10:28:34] optimpv.axBOtorchOptimizer: Finished Turbo batch 75 with 2 trials with current best value: 2.18e-02, TR length: 1.22e-01
[INFO 12-05 10:28:35] optimpv.axBOtorchOptimizer: Finished Turbo batch 76 with 2 trials with current best value: 2.18e-02, TR length: 1.22e-01
[INFO 12-05 10:28:35] optimpv.axBOtorchOptimizer: Finished Turbo batch 77 with 2 trials with current best value: 2.18e-02, TR length: 1.22e-01
[INFO 12-05 10:28:36] optimpv.axBOtorchOptimizer: Finished Turbo batch 78 with 2 trials with current best value: 2.18e-02, TR length: 1.22e-01
[INFO 12-05 10:28:36] optimpv.axBOtorchOptimizer: Finished Turbo batch 79 with 2 trials with current best value: 2.18e-02, TR length: 1.22e-01
[INFO 12-05 10:28:37] optimpv.axBOtorchOptimizer: Finished Turbo batch 80 with 2 trials with current best value: 2.18e-02, TR length: 1.22e-01
[INFO 12-05 10:28:38] optimpv.axBOtorchOptimizer: Finished Turbo batch 81 with 2 trials with current best value: 2.18e-02, TR length: 1.22e-01
[INFO 12-05 10:28:38] optimpv.axBOtorchOptimizer: Finished Turbo batch 82 with 2 trials with current best value: 2.18e-02, TR length: 1.22e-01
[INFO 12-05 10:28:39] optimpv.axBOtorchOptimizer: Finished Turbo batch 83 with 2 trials with current best value: 2.18e-02, TR length: 1.22e-01
[INFO 12-05 10:28:39] optimpv.axBOtorchOptimizer: Finished Turbo batch 84 with 2 trials with current best value: 2.18e-02, TR length: 1.22e-01
[INFO 12-05 10:28:40] optimpv.axBOtorchOptimizer: Finished Turbo batch 85 with 2 trials with current best value: 2.18e-02, TR length: 7.63e-02
[INFO 12-05 10:28:40] optimpv.axBOtorchOptimizer: Finished Turbo batch 86 with 2 trials with current best value: 2.18e-02, TR length: 7.63e-02
[INFO 12-05 10:28:41] optimpv.axBOtorchOptimizer: Finished Turbo batch 87 with 2 trials with current best value: 2.18e-02, TR length: 7.63e-02
[INFO 12-05 10:28:42] optimpv.axBOtorchOptimizer: Finished Turbo batch 88 with 2 trials with current best value: 2.18e-02, TR length: 7.63e-02
[INFO 12-05 10:28:42] optimpv.axBOtorchOptimizer: Finished Turbo batch 89 with 2 trials with current best value: 2.18e-02, TR length: 7.63e-02
[INFO 12-05 10:28:43] optimpv.axBOtorchOptimizer: Finished Turbo batch 90 with 2 trials with current best value: 2.18e-02, TR length: 7.63e-02
[INFO 12-05 10:28:43] optimpv.axBOtorchOptimizer: Finished Turbo batch 91 with 2 trials with current best value: 2.17e-02, TR length: 7.63e-02
[INFO 12-05 10:28:44] optimpv.axBOtorchOptimizer: Finished Turbo batch 92 with 2 trials with current best value: 2.17e-02, TR length: 7.63e-02
[INFO 12-05 10:28:45] optimpv.axBOtorchOptimizer: Finished Turbo batch 93 with 2 trials with current best value: 2.17e-02, TR length: 7.63e-02
[INFO 12-05 10:28:45] optimpv.axBOtorchOptimizer: Finished Turbo batch 94 with 2 trials with current best value: 2.17e-02, TR length: 7.63e-02
[INFO 12-05 10:28:46] optimpv.axBOtorchOptimizer: Finished Turbo batch 95 with 2 trials with current best value: 2.17e-02, TR length: 7.63e-02
[INFO 12-05 10:28:46] optimpv.axBOtorchOptimizer: Finished Turbo batch 96 with 2 trials with current best value: 2.17e-02, TR length: 7.63e-02
[INFO 12-05 10:28:47] optimpv.axBOtorchOptimizer: Finished Turbo batch 97 with 2 trials with current best value: 2.17e-02, TR length: 7.63e-02
[INFO 12-05 10:28:48] optimpv.axBOtorchOptimizer: Finished Turbo batch 98 with 2 trials with current best value: 2.17e-02, TR length: 7.63e-02
[INFO 12-05 10:28:48] optimpv.axBOtorchOptimizer: Finished Turbo batch 99 with 2 trials with current best value: 2.17e-02, TR length: 7.63e-02
[INFO 12-05 10:28:49] optimpv.axBOtorchOptimizer: Finished Turbo batch 100 with 2 trials with current best value: 2.17e-02, TR length: 7.63e-02
[INFO 12-05 10:28:50] optimpv.axBOtorchOptimizer: Finished Turbo batch 101 with 2 trials with current best value: 2.17e-02, TR length: 4.77e-02
[INFO 12-05 10:28:50] optimpv.axBOtorchOptimizer: Finished Turbo batch 102 with 2 trials with current best value: 2.17e-02, TR length: 4.77e-02
[INFO 12-05 10:28:51] optimpv.axBOtorchOptimizer: Turbo is terminated.

[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(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.37e-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 = 1.11e+04 s$^{-1}$ (original: 1.60e+04 s$^{-1}$), bounds: [1.00e+04, 1.00e+06] s$^{-1}$
k_c = 1.30e+04 s$^{-1}$ (original: 1.30e+06 s$^{-1}$), bounds: [1.00e+04, 1.00e+08] s$^{-1}$
k_e = 1.04e+04 m$^{3}$ s$^{-1}$ (original: 7.50e+05 m$^{3}$ s$^{-1}$), bounds: [1.00e+04, 1.00e+08] m$^{3}$ s$^{-1}$
mu = 7.10e-05 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 = 6.55e-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 = 3.19e-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.32e-28 - (original: 1.00e-28 -), bounds: [1.00e-29, 1.00e-27] -
N_A = 6.26e+21 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.021723093489163972 at iteration  186
../_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

# Plot the fitted data
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_11_0.png 
[10]:

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='dogbox', options={'max_nfev':200}, name='scipy_opti', parallel_agents=True, max_parallelism=os.cpu_count()-1, verbose_logging=True, )

[11]:

scipyOpti.optimize_least_squares() # run the optimization with scipy

Starting optimization using dogbox method
Optimization completed with status: The maximum number of function evaluations is exceeded.
Final objective value: [0.02171124]

[11]:

     message: The maximum number of function evaluations is exceeded.
     success: False
      status: 0
         fun: [ 2.171e-02]
           x: [-1.603e+01  4.044e+00  4.108e+00  4.031e+00 -4.139e+00
               -2.185e+00 -1.497e+00 -2.788e+01  2.179e+01]
        cost: 0.00023568902102490268
         jac: [[ 8.073e-06 -8.725e-07  2.790e-05 -1.309e-05 -2.357e-05
                 3.002e-07 -1.020e-05  2.789e-05  2.380e-05]]
        grad: [ 1.753e-07 -1.894e-08  6.057e-07 -2.841e-07 -5.118e-07
                6.518e-09 -2.215e-07  6.055e-07  5.168e-07]
  optimality: 6.057454955461765e-07
 active_mask: [0 0 0 0 0 0 0 0 0]
        nfev: 200
        njev: 191

[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('Best parameters:')
for p,po in zip(scipyOpti.params, params_orig):
    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.44e-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 = 1.11e+04 s$^{-1}$ (original: 1.60e+04 s$^{-1}$), bounds: [1.00e+04, 1.00e+06] s$^{-1}$
k_c = 1.28e+04 s$^{-1}$ (original: 1.30e+06 s$^{-1}$), bounds: [1.00e+04, 1.00e+08] s$^{-1}$
k_e = 1.07e+04 m$^{3}$ s$^{-1}$ (original: 7.50e+05 m$^{3}$ s$^{-1}$), bounds: [1.00e+04, 1.00e+08] m$^{3}$ s$^{-1}$
mu = 7.26e-05 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 = 6.53e-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 = 3.18e-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.32e-28 - (original: 1.00e-28 -), bounds: [1.00e-29, 1.00e-27] -
N_A = 6.19e+21 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