Fit of transient photoluminescence (TrPL) with diffusion and multi trap 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 M. Simmonds.

# 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 copy import deepcopy
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 *

Define the parameters for the simulation

# Define the parameters to be fitted
params = []

Eg = FitParam(name = 'Eg', value = 1.553, bounds = [0.5,2.0], log_scale = False, rescale = True, value_type = 'float', type='fixed', display_name=r'$E_g$', unit='eV', axis_type = 'linear')
params.append(Eg)

L = FitParam(name = 'L', value = 450e-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)

alpha = FitParam(name = 'alpha', value = 64348.30886337494 *1e2, 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)

N_cv = FitParam(name = 'N_cv', value = 2e24, bounds = [1e19,1e25], log_scale = True, rescale = True, value_type = 'float', type='fixed', display_name=r'$N_{cv}$', unit='m$^{-3}$', axis_type = 'log',force_log=True)
params.append(N_cv)

k_direct = FitParam(name = 'k_direct', value = 1.96e-17, bounds = [1e-18,1e-14], 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)

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

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

N_t_bulk_1 = FitParam(name = 'N_t_bulk_1', value = 1.85e23, bounds = [1e19,5e23], log_scale = True, rescale = True, value_type = 'float', type='range', display_name=r'$N_{t,\text{bulk}}$', unit='m$^{-3}$', axis_type = 'log',force_log=True)
params.append(N_t_bulk_1)

C_n_1 = FitParam(name = 'C_n_1', value = 4.24e-15, bounds = [1e-19,1e-12], log_scale = True, rescale = True, value_type = 'float', type='range', display_name=r'$C_{n,1}$', unit='m$^{3}$ s$^{-1}$', axis_type = 'log',force_log=True)
params.append(C_n_1)

C_p_1 = FitParam(name = 'C_p_1', value = 8.85e-19, bounds = [1e-19,1e-12], log_scale = True, rescale = True, value_type = 'float', type='range', display_name=r'$C_{p,1}$', unit='m$^{3}$ s$^{-1}$', axis_type = 'log',force_log=True)
params.append(C_p_1)

E_t_bulk_1 = FitParam(name = 'E_t_bulk_1', value = 0.2, bounds = [0.05,0.3], log_scale = True, rescale = True, value_type = 'float', type='range', display_name=r'$E_{t,\text{bulk}}$', unit='m$^{3}$ s$^{-1}$', axis_type = 'log',force_log=False)
params.append(E_t_bulk_1)

N_t_bulk_2 = FitParam(name = 'N_t_bulk_2', value = 1.36e21, bounds = [1e19,5e23], log_scale = True, rescale = True, value_type = 'float', type='range', display_name=r'$N_{t,\text{bulk}}$', unit='m$^{-3}$', axis_type = 'log',force_log=True)
params.append(N_t_bulk_2)

C_n_2 = FitParam(name = 'C_n_2', value = 1.72e-13, bounds = [1e-19,1e-12], log_scale = True, rescale = True, value_type = 'float', type='range', display_name=r'$C_{n,2}$', unit='m$^{3}$ s$^{-1}$', axis_type = 'log',force_log=True)
params.append(C_n_2)

C_p_2 = FitParam(name = 'C_p_2', value = 1.13e-16, bounds = [1e-19,1e-12], log_scale = True, rescale = True, value_type = 'float', type='range', display_name=r'$C_{p,2}$', unit='m$^{3}$ s$^{-1}$', axis_type = 'log',force_log=True)
params.append(C_p_2)

E_t_bulk_2 = FitParam(name = 'E_t_bulk_2', value = 1.34, bounds = [0.3,Eg.value-0.1], log_scale = True, rescale = True, value_type = 'float', type='range', display_name=r'$E_{t,\text{bulk}}$', unit='m$^{3}$ s$^{-1}$', axis_type = 'log',force_log=False)
params.append(E_t_bulk_2)

I_factor_PL = FitParam(name = 'I_factor_PL', value = 1.275e-22, bounds = [1e-27,1e-20], log_scale = True, rescale = True, value_type = 'float', type='fixed', display_name=r'$I_{\text{PL}}$', unit='-', axis_type = 'log', force_log=True)
params.append(I_factor_PL) # in the following we weill fit the PL with the normalized log transformation so this factor is not useful and can be fixed to any value

N0 = FitParam(name = 'N0', value = 1.39e+20, bounds = [1e20,1.5e20], log_scale = False, rescale = True, value_type = 'float', type='range', display_name=r'$N_0$', unit='m$^{-3}$', axis_type = 'log',force_log=False)
params.append(N0)

# original values
params_orig = 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

# 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','FAPI_trPL')
filenames = ['2M11-FAPI27_3600s_10kHz_ND3_1280cps_0.272112uW.dat','1M3-FAPI27_3600s_10kHz_ND2_500cps_0.0737403uW.dat','0M4-FAPI27_3600s_10kHz_ND1_940cps_0.048888uW.dat',] #'3M8-FAPI27_3600s_10kHz_ND5_570cps_1.2156099999999999uW.dat',
fig, axes = plt.subplots(nrows=1, ncols=len(filenames), figsize=(12, 4))
power = []
pmax = 1.2156099999999999
pmax = max([float(file.split('_')[-1].replace('uW.dat', '')) for file in filenames])
for idx, file in enumerate(filenames):
    data_raw = pd.read_csv(os.path.join(path2data,file), sep=r'\s+',names=['idx','t',"trPL"], skiprows=1)
    data_raw = data_raw[['t', 'trPL']]
    # convert time to seconds
    data_raw['t'] = data_raw['t'] * 1e-12 # convert to seconds
    max_idx = data_raw['trPL'].idxmax()

    # remove everything before the maximum
    data_raw = data_raw.iloc[max_idx:]
    # 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['t'] = data_raw['t'] - data_raw['t'][max_idx]
    data_raw = data_raw[data_raw['t'] >= 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(data_raw['t'][1]), 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['trPL'])

    power.append(float(file.split('_')[-1].replace('uW.dat', '')))
    if idx == 0:
        data2fit = {'t': t_log , 'trPL': trPL_log, 'G_frac': power[-1] / pmax * np.ones_like(t_log)}
    else:
        data2fit['t'] = np.concatenate((data2fit['t'], t_log ))
        data2fit['trPL'] = np.concatenate((data2fit['trPL'], trPL_log))
        data2fit['G_frac'] = np.concatenate((data2fit['G_frac'], power[-1] / pmax * np.ones_like(t_log)))

    # plot the data
    ax = axes[idx // len(axes[0]), idx % len(axes[0])] if len(axes.shape) > 1 else axes[idx]
    ax.plot(data_raw['t'] - data_raw['t'].min(), data_raw['trPL'], 'o', label='raw')
    ax.plot(t_log, trPL_log, '*', label='Interpolated trPL')
    ax.set_xscale('log')
    ax.set_yscale('log')
    ax.set_xlabel('Time [s]')
    ax.set_ylabel('trPL [a.u.]')
    ax.set_title(f'P = {power[-1]:.2e} $\\mu$W')
    ax.legend()
plt.tight_layout()
plt.show()
data2fit= pd.DataFrame(data2fit)

# Plot the data to be fitted and the initial guess
time = data2fit['t'].values # time in seconds
X = np.asarray(data2fit[['t', 'G_frac']])
y = np.asarray(data2fit['trPL'])
fpu = 10e3 # Frequency of the pump laser in Hz
N0 = 1.39e+20
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 = { 'fpu': fpu , 'background' : background, 'N0': N0,}
exp_format = 'trPL' # experiment format
model = partial(DBTD_multi_trap, method='LSODA', dimensionless=False, timeout=30, timeout_solve=30)
RateEq = RateEqAgent(params, [X], [y], model = 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 =['normalize','log'],do_G_frac_transform=True,parallel=False)

y_test = RateEq.run({},exp_format=exp_format) # get the initial prediction with the original parameters


from optimpv.general.general import *
y_transformed, y_pred_transformed = transform_data(y,y_test, transforms=['normalize','log'], do_G_frac_transform=True, X=X)


plt.figure(figsize=(10,10))

viridis = plt.get_cmap('viridis', len(np.unique(X[:, 1])))
plt.figure()
for i, G in enumerate(np.unique(X[:, 1])):
    plt.plot(X[:, 0][X[:, 1] == G], 10**y_transformed[X[:, 1] == G],'o',color=viridis(i),alpha=0.1)
    plt.plot(X[:, 0][X[:, 1] == G], 10**y_pred_transformed[X[:, 1] == G], color=viridis(i), linestyle='-', label=f'{G:.2f}')

plt.xscale('log')
plt.yscale('log')
plt.xlabel('Time [s]')
plt.ylabel(exp_format + ' [a.u.]')
plt.legend(title='Intensity fraction:')
plt.show()

Run the optimization

from optimpv.axBOtorch.axBOtorchOptimizer import axBOtorchOptimizer
from optimpv.axBOtorch.axUtils import get_VMLC_default_model_kwargs_list

# Here we add some constraints to the parameters to help the optimizer
parameter_constraints = [f' -N_t_bulk_1 - 0.5 * C_n_1 - 0.5 * C_p_1  - N_t_bulk_2 - 0.5 * C_n_2 - 0.5 * C_p_2<= -5']

optimizer = axBOtorchOptimizer(params = params, agents = RateEq, models = ['SOBOL','BOTORCH_MODULAR'],n_batches = [1,150], batch_size = [10,4], max_parallelism = 100, model_kwargs_list = get_VMLC_default_model_kwargs_list(num_free_params), parameter_constraints = parameter_constraints,parallel_agents= True)

%%time
optimizer.optimize_turbo(force_continue=False,kwargs_turbo_state={'failure_tolerance':4}) # run the optimization with turbo

# 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}')

# 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()

# 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(np.unique(X[:, 1])))
plt.figure()
for i, G in enumerate(np.unique(X[:, 1])):
    plt.plot(X[:, 0][X[:, 1] == G], y[X[:, 1] == G]/max(y[X[:, 1] == G]),'o',color=viridis(i),alpha=0.05)
    plt.plot(X[:, 0][X[:, 1] == G], yfit[X[:, 1] == G]/max(yfit[X[:, 1] == G]), color=viridis(i), linestyle='-', label=f'G = {G}')

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