Fitting with time delay scan (model: sum of exponential decay and damped oscillation)¶
Objective¶
Fitting with sum of exponential decay model and damped oscillation model
Save and Load fitting result
Calculates species associated coefficent from fitting result
Evaluates F-test based confidence interval
In this example, we only deal with gaussian irf
# import needed module
import numpy as np
import matplotlib.pyplot as plt
import TRXASprefitpack
from TRXASprefitpack import solve_seq_model, rate_eq_conv, dmp_osc_conv_gau
plt.rcParams["figure.figsize"] = (12,9)
Version information¶
print(TRXASprefitpack.__version__)
0.6.0
Detecting oscillation feature¶
# Generates fake experiment data
# Model: 1 -> 2 -> 3 -> GS
# lifetime tau1: 500 fs, tau2: 10 ps, tau3: 1000 ps
# oscillation: tau_osc: 1 ps, period_osc: 300 fs, phase: pi/4
# fwhm paramter of gaussian IRF: 100 fs
tau_1 = 0.5
tau_2 = 10
tau_3 = 1000
fwhm = 0.100
tau_osc = 1
period_osc = 0.3
phase = np.pi/4
# initial condition
y0 = np.array([1, 0, 0, 0])
# set time range (mixed step)
t_seq1 = np.arange(-2, -1, 0.2)
t_seq2 = np.arange(-1, 1, 0.02)
t_seq3 = np.arange(1, 5, 0.2)
t_seq4 = np.arange(5, 10, 1)
t_seq5 = np.arange(10, 100, 10)
t_seq6 = np.arange(100, 1000, 100)
t_seq7 = np.linspace(1000, 2000, 2)
t_seq = np.hstack((t_seq1, t_seq2, t_seq3, t_seq4, t_seq5, t_seq6, t_seq7))
eigval_seq, V_seq, c_seq = solve_seq_model(np.array([tau_1, tau_2, tau_3]), y0)
# Now generates measured transient signal
# Last element is ground state
abs_1 = [1, 1, 1, 0]; abs_1_osc = 0.05
abs_2 = [0.5, 0.8, 0.2, 0]; abs_2_osc = 0.001
abs_3 = [-0.5, 0.7, 0.9, 0]; abs_3_osc = -0.002
abs_4 = [0.6, 0.3, -1, 0]; abs_4_osc = 0.0018
t0 = np.random.normal(0, fwhm, 4) # perturb time zero of each scan
# generate measured data
y_obs_1 = rate_eq_conv(t_seq-t0[0], fwhm, abs_1, eigval_seq, V_seq, c_seq, irf='g')+\
abs_1_osc*dmp_osc_conv_gau(t_seq-t0[0], fwhm, 1/tau_osc, period_osc, phase)
y_obs_2 = rate_eq_conv(t_seq-t0[1], fwhm, abs_2, eigval_seq, V_seq, c_seq, irf='g')+\
abs_2_osc*dmp_osc_conv_gau(t_seq-t0[1], fwhm, 1/tau_osc, period_osc, phase)
y_obs_3 = rate_eq_conv(t_seq-t0[2], fwhm, abs_3, eigval_seq, V_seq, c_seq, irf='g')+\
abs_3_osc*dmp_osc_conv_gau(t_seq-t0[2], fwhm, 1/tau_osc, period_osc, phase)
y_obs_4 = rate_eq_conv(t_seq-t0[3], fwhm, abs_4, eigval_seq, V_seq, c_seq, irf='g')+\
abs_4_osc*dmp_osc_conv_gau(t_seq-t0[3], fwhm, 1/tau_osc, period_osc, phase)
# generate random noise with (S/N = 200)
# Define noise level (S/N=200) w.r.t peak
eps_obs_1 = np.max(np.abs(y_obs_1))/200*np.ones_like(y_obs_1)
eps_obs_2 = np.max(np.abs(y_obs_2))/200*np.ones_like(y_obs_2)
eps_obs_3 = np.max(np.abs(y_obs_3))/200*np.ones_like(y_obs_3)
eps_obs_4 = np.max(np.abs(y_obs_4))/200*np.ones_like(y_obs_4)
# generate random noise
noise_1 = np.random.normal(0, eps_obs_1, t_seq.size)
noise_2 = np.random.normal(0, eps_obs_2, t_seq.size)
noise_3 = np.random.normal(0, eps_obs_3, t_seq.size)
noise_4 = np.random.normal(0, eps_obs_4, t_seq.size)
# generate measured intensity
i_obs_1 = y_obs_1 + noise_1
i_obs_2 = y_obs_2 + noise_2
i_obs_3 = y_obs_3 + noise_3
i_obs_4 = y_obs_4 + noise_4
# print real values
print('-'*24)
print(f'fwhm: {fwhm}')
print(f'tau_1: {tau_1}')
print(f'tau_2: {tau_2}')
print(f'tau_3: {tau_3}')
print(f'tau_osc: {tau_osc}')
print(f'period_osc: {period_osc}')
print(f'phase_osc: {phase}')
for i in range(4):
print(f't_0_{i+1}: {t0[i]}')
print('-'*24)
print('Excited Species contribution')
print(f'scan 1: {abs_1[0]} \t {abs_1[1]} \t {abs_1[2]}')
print(f'scan 2: {abs_2[0]} \t {abs_2[1]} \t {abs_2[2]}')
print(f'scan 3: {abs_3[0]} \t {abs_3[1]} \t {abs_3[2]}')
print(f'scan 4: {abs_4[0]} \t {abs_4[1]} \t {abs_4[2]}')
param_exact = [fwhm, t0[0], t0[1], t0[2], t0[3], tau_1, tau_2, tau_3, tau_osc, period_osc, phase]
data1 = np.vstack((t_seq, i_obs_1, eps_obs_1)).T
data2 = np.vstack((t_seq, i_obs_2, eps_obs_2)).T
data3 = np.vstack((t_seq, i_obs_3, eps_obs_3)).T
data4 = np.vstack((t_seq, i_obs_4, eps_obs_4)).T
np.savetxt('./fit_tscan/osc/example_osc_1.txt', data1)
np.savetxt('./fit_tscan/osc/example_osc_2.txt', data2)
np.savetxt('./fit_tscan/osc/example_osc_3.txt', data3)
np.savetxt('./fit_tscan/osc/example_osc_4.txt', data4)
------------------------
fwhm: 0.1
tau_1: 0.5
tau_2: 10
tau_3: 1000
tau_osc: 1
period_osc: 0.3
phase_osc: 0.7853981633974483
t_0_1: 0.05053956334459484
t_0_2: -0.07163546673750752
t_0_3: -0.004390946806632458
t_0_4: 0.13670418243246096
------------------------
Excited Species contribution
scan 1: 1 1 1
scan 2: 0.5 0.8 0.2
scan 3: -0.5 0.7 0.9
scan 4: 0.6 0.3 -1
# plot model experimental data
plt.errorbar(t_seq, i_obs_1, eps_obs_1, label='1')
plt.errorbar(t_seq, i_obs_2, eps_obs_2, label='2')
plt.errorbar(t_seq, i_obs_3, eps_obs_3, label='3')
plt.errorbar(t_seq, i_obs_4, eps_obs_4, label='4')
plt.legend()
plt.show()
plt.errorbar(t_seq, i_obs_1, eps_obs_1, label='1')
plt.errorbar(t_seq, i_obs_2, eps_obs_2, label='2')
plt.errorbar(t_seq, i_obs_3, eps_obs_3, label='3')
plt.errorbar(t_seq, i_obs_4, eps_obs_4, label='4')
plt.legend()
plt.xlim(-10*fwhm, 20*fwhm)
plt.show()
We can show oscillation feature at scan 1. First try fitting without oscillation.
# import needed module for fitting
from TRXASprefitpack import fit_transient_exp
# time, intensity, eps should be sequence of numpy.ndarray
t = [t_seq]
intensity = [np.vstack((i_obs_1, i_obs_2, i_obs_3, i_obs_4)).T]
eps = [np.vstack((eps_obs_1, eps_obs_2, eps_obs_3, eps_obs_4)).T]
# set initial guess
irf = 'g' # shape of irf function
fwhm_init = 0.15
t0_init = np.array([0, 0, 0, 0])
# test with one decay module
tau_init = np.array([0.2, 20, 1500])
fit_result_decay = fit_transient_exp(irf, fwhm_init, t0_init, tau_init, False, do_glb=True, t=t, intensity=intensity, eps=eps)
# print fitting result
print(fit_result_decay)
[Model information]
model : decay
irf: g
fwhm: 0.1008
eta: 0.0000
base: False
[Optimization Method]
global: basinhopping
leastsq: trf
[Optimization Status]
nfev: 3456
status: 0
global_opt msg: requested number of basinhopping iterations completed successfully
leastsq_opt msg: `ftol` termination condition is satisfied.
[Optimization Results]
Total Data points: 600
Number of effective parameters: 20
Degree of Freedom: 580
Chi squared: 1060.8466
Reduced chi squared: 1.829
AIC (Akaike Information Criterion statistic): 381.9358
BIC (Bayesian Information Criterion statistic): 469.8743
[Parameters]
fwhm_G: 0.10077973 +/- 0.00091351 ( 0.91%)
t_0_1_1: 0.05003520 +/- 0.00042857 ( 0.86%)
t_0_1_2: -0.07119986 +/- 0.00060856 ( 0.85%)
t_0_1_3: -0.00538671 +/- 0.00076649 ( 14.23%)
t_0_1_4: 0.13664093 +/- 0.00066225 ( 0.48%)
tau_1: 0.50197077 +/- 0.00270331 ( 0.54%)
tau_2: 10.00304080 +/- 0.05645921 ( 0.56%)
tau_3: 1003.82449967 +/- 4.76296335 ( 0.47%)
[Parameter Bound]
fwhm_G: 0.075 <= 0.10077973 <= 0.3
t_0_1_1: -0.3 <= 0.05003520 <= 0.3
t_0_1_2: -0.3 <= -0.07119986 <= 0.3
t_0_1_3: -0.3 <= -0.00538671 <= 0.3
t_0_1_4: -0.3 <= 0.13664093 <= 0.3
tau_1: 0.075 <= 0.50197077 <= 1.2
tau_2: 4.8 <= 10.00304080 <= 76.8
tau_3: 307.2 <= 1003.82449967 <= 4915.2
[Component Contribution]
DataSet dataset_1:
#tscan tscan_1 tscan_2 tscan_3 tscan_4
decay 1 -1.29% -28.41% -51.13% 8.89%
decay 2 -0.71% 54.18% -9.62% 52.56%
decay 3 97.99% 17.41% 39.25% -38.55%
[Parameter Correlation]
Parameter Correlations > 0.1 are reported.
(t_0_1_1, fwhm_G) = 0.106
(t_0_1_2, fwhm_G) = 0.103
(tau_1, t_0_1_2) = 0.118
(tau_1, t_0_1_3) = -0.347
(tau_2, t_0_1_2) = 0.101
(tau_2, t_0_1_4) = 0.145
(tau_3, tau_2) = -0.156
# plot fitting result and experimental data
color_lst = ['red', 'blue', 'green', 'black']
for i in range(4):
plt.errorbar(t[0], intensity[0][:, i], eps[0][:, i], label=f'expt {i+1}', color=color_lst[i])
plt.errorbar(t[0], fit_result_decay['fit'][0][:, i], label=f'fit {i+1}', color=color_lst[i])
plt.legend()
plt.show()
# plot with shorter time range
for i in range(4):
plt.errorbar(t[0], intensity[0][:, i], eps[0][:, i], label=f'expt {i+1}', color=color_lst[i])
plt.errorbar(t[0], fit_result_decay['fit'][0][:, i], label=f'fit {i+1}', color=color_lst[i])
plt.legend()
plt.xlim(-10*fwhm_init, 20*fwhm_init)
plt.show()
There may exists oscillation in experimental scan 1. To show oscillation feature plot residual (expt-fit)
# To show oscillation feature plot residual
for i in range(4):
plt.errorbar(t[0], fit_result_decay['res'][0][:, i], eps[0][:, i], label=f'res {i+1}', color=color_lst[i])
plt.legend()
plt.xlim(-10*fwhm_init, 20*fwhm_init)
plt.show()
Only residual for experimental scan 1 shows clear oscillation feature, Now add oscillation feature.
from TRXASprefitpack import fit_transient_both
tau_osc_init = np.array([1.5])
period_osc_init = np.array([0.5])
fit_result_decay_osc = fit_transient_both(irf, fwhm_init, t0_init, tau_init,
tau_osc_init, period_osc_init,
False, do_glb=True, t=t, intensity=intensity, eps=eps)
# print fitting result
print(fit_result_decay_osc)
[Model information]
model : both
irf: g
fwhm: 0.0998
eta: 0.0000
base: False
[Optimization Method]
global: basinhopping
leastsq: trf
[Optimization Status]
nfev: 13494
status: 0
global_opt msg: requested number of basinhopping iterations completed successfully
leastsq_opt msg: `ftol` termination condition is satisfied.
[Optimization Results]
Total Data points: 600
Number of effective parameters: 30
Degree of Freedom: 570
Chi squared: 638.1717
Reduced chi squared: 1.1196
AIC (Akaike Information Criterion statistic): 97.0066
BIC (Bayesian Information Criterion statistic): 228.9145
[Parameters]
fwhm_G: 0.09984398 +/- 0.00075909 ( 0.76%)
t_0_1_1: 0.05061665 +/- 0.00039406 ( 0.78%)
t_0_1_2: -0.07132916 +/- 0.00051803 ( 0.73%)
t_0_1_3: -0.00474142 +/- 0.00064774 ( 13.66%)
t_0_1_4: 0.13670938 +/- 0.00056220 ( 0.41%)
tau_1: 0.50172568 +/- 0.00215798 ( 0.43%)
tau_2: 10.00631278 +/- 0.04434765 ( 0.44%)
tau_3: 1002.83212275 +/- 3.72071882 ( 0.37%)
tau_osc_1: 1.28461622 +/- 0.25058768 ( 19.51%)
period_osc_1: 0.29797454 +/- 0.00220269 ( 0.74%)
[Parameter Bound]
fwhm_G: 0.075 <= 0.09984398 <= 0.3
t_0_1_1: -0.3 <= 0.05061665 <= 0.3
t_0_1_2: -0.3 <= -0.07132916 <= 0.3
t_0_1_3: -0.3 <= -0.00474142 <= 0.3
t_0_1_4: -0.3 <= 0.13670938 <= 0.3
tau_1: 0.075 <= 0.50172568 <= 1.2
tau_2: 4.8 <= 10.00631278 <= 76.8
tau_3: 307.2 <= 1002.83212275 <= 4915.2
tau_osc_1: 0.3 <= 1.28461622 <= 4.8
period_osc_1: 0.075 <= 0.29797454 <= 1.2
[Phase Factor]
DataSet dataset_1:
#tscan tscan_1 tscan_2 tscan_3 tscan_4
dmp_osc 1 0.2120 π 0.4132 π 0.8410 π 0.2091 π
[Component Contribution]
DataSet dataset_1:
#tscan tscan_1 tscan_2 tscan_3 tscan_4
decay 1 0.01% -28.40% -51.00% 8.90%
decay 2 -1.21% 54.11% -9.61% 52.46%
decay 3 94.45% 17.39% 39.18% -38.48%
dmp_osc 1 4.34% 0.11% 0.22% 0.16%
[Parameter Correlation]
Parameter Correlations > 0.1 are reported.
(t_0_1_1, fwhm_G) = 0.241
(t_0_1_2, fwhm_G) = 0.154
(t_0_1_4, fwhm_G) = 0.116
(tau_1, t_0_1_2) = 0.111
(tau_1, t_0_1_3) = -0.348
(tau_2, t_0_1_2) = 0.103
(tau_2, t_0_1_4) = 0.133
(tau_3, tau_2) = -0.156
(period_osc_1, fwhm_G) = -0.251
(period_osc_1, t_0_1_1) = -0.451
# plot residual and oscilation fit
for i in range(1):
plt.errorbar(t[0], intensity[0][:, i]-fit_result_decay_osc['fit_decay'][0][:, i], eps[0][:, i], label=f'res {i+1}', color='black')
plt.errorbar(t[0], fit_result_decay_osc['fit_osc'][0][:, i], label=f'osc {i+1}', color='red')
plt.legend()
plt.xlim(-10*fwhm_init, 20*fwhm_init)
plt.show()
print()
# Compare fitting value and exact value
for i in range(len(fit_result_decay_osc['x'])):
print(f"{fit_result_decay_osc['param_name'][i]}: {fit_result_decay_osc['x'][i]} (fit) \t {param_exact[i]} (exact)")
fwhm_G: 0.09984398023073539 (fit) 0.1 (exact)
t_0_1_1: 0.050616647691810734 (fit) 0.05053956334459484 (exact)
t_0_1_2: -0.07132916197304762 (fit) -0.07163546673750752 (exact)
t_0_1_3: -0.004741422504451592 (fit) -0.004390946806632458 (exact)
t_0_1_4: 0.1367093818502816 (fit) 0.13670418243246096 (exact)
tau_1: 0.5017256765256979 (fit) 0.5 (exact)
tau_2: 10.006312781399181 (fit) 10 (exact)
tau_3: 1002.8321227541483 (fit) 1000 (exact)
tau_osc_1: 1.284616223329954 (fit) 1 (exact)
period_osc_1: 0.29797453736969626 (fit) 0.3 (exact)
# save fitting result to file
from TRXASprefitpack import save_TransientResult, load_TransientResult
save_TransientResult(fit_result_decay_osc, 'example_decay_osc') # save fitting result to example_decay_2.h5
loaded_result = load_TransientResult('example_decay_osc') # load fitting result from example_decay_2.h5
Now deduce species associated difference coefficient from sequential decay model
y0 = np.array([1, 0, 0, 0]) # initial cond
eigval, V, c = solve_seq_model(loaded_result['x'][5:-2], y0)
# compute scaled V matrix
V_scale = np.einsum('j,ij->ij', c, V)
diff_abs_fit = np.linalg.solve(V_scale[:-1, :-1].T, loaded_result['c'][0][:-1,:])
# slice last column and row corresponding to ground state
# exclude oscillation factor
# compare with exact result
print('-'*24)
print('[Species Associated Difference Coefficent]')
print('scan # \t ex 1 (fit) \t ex 1 (exact) \t ex 2 (fit) \t ex 2 (exact) \t ex 3 (exact)')
print(f'1 \t {diff_abs_fit[0,0]} \t {abs_1[0]} \t {diff_abs_fit[1,0]} \t {abs_1[1]} \t {diff_abs_fit[2,0]} \t {abs_1[2]}')
print(f'2 \t {diff_abs_fit[0,1]} \t {abs_2[0]} \t {diff_abs_fit[1,1]} \t {abs_2[1]} \t {diff_abs_fit[2,1]} \t {abs_2[2]}')
print(f'3 \t {diff_abs_fit[0,2]} \t {abs_3[0]} \t {diff_abs_fit[1,2]} \t {abs_3[1]} \t {diff_abs_fit[2,2]} \t {abs_3[2]}')
print(f'4 \t {diff_abs_fit[0,3]} \t {abs_4[0]} \t {diff_abs_fit[1,3]} \t {abs_4[1]} \t {diff_abs_fit[2,3]} \t {abs_4[2]}')
------------------------
[Species Associated Difference Coefficent]
scan # ex 1 (fit) ex 1 (exact) ex 2 (fit) ex 2 (exact) ex 3 (exact)
1 0.999814728016317 1 0.9998813520856183 1 1.0020621741915658 1
2 0.5013994541791791 0.5 0.8002106071258839 0.8 0.20017964124729995 0.2
3 -0.498563790521078 -0.5 0.6989863593941319 0.7 0.9022594839616441 0.9
4 0.6007902340276545 0.6 0.29855248807181206 0.3 -1.0000321282616873 -1
It also matches well, as expected.
Now calculates F-test based confidence interval.
from TRXASprefitpack import confidence_interval
ci_result = confidence_interval(loaded_result, 0.05) # set significant level: 0.05 -> 95% confidence level
print(ci_result) # report confidence interval
[Report for Confidence Interval]
Method: f
Significance level: 5.000000e-02
[Confidence interval]
0.09984398 - 0.00147569 <= b'fwhm_G' <= 0.09984398 + 0.00147812
0.05061665 - 0.00076623 <= b't_0_1_1' <= 0.05061665 + 0.00076005
-0.07132916 - 0.00102483 <= b't_0_1_2' <= -0.07132916 + 0.00099725
-0.00474142 - 0.00126688 <= b't_0_1_3' <= -0.00474142 + 0.00127079
0.13670938 - 0.00110796 <= b't_0_1_4' <= 0.13670938 + 0.001104
0.50172568 - 0.00418296 <= b'tau_1' <= 0.50172568 + 0.00422175
10.00631278 - 0.08543495 <= b'tau_2' <= 10.00631278 + 0.08629124
1002.83212275 - 7.26121914 <= b'tau_3' <= 1002.83212275 + 7.33143627
1.28461622 - 0.40931793 <= b'tau_osc_1' <= 1.28461622 + 0.75237185
0.29797454 - 0.00397536 <= b'period_osc_1' <= 0.29797454 + 0.00455939
# compare with ase
from scipy.stats import norm
factor = norm.ppf(1-0.05/2)
print('[Confidence interval (from ASE)]')
for i in range(loaded_result['param_name'].size):
print(f"{loaded_result['x'][i]: .8f} - {factor*loaded_result['x_eps'][i] :.8f}",
f"<= {loaded_result['param_name'][i]} <=", f"{loaded_result['x'][i] :.8f} + {factor*loaded_result['x_eps'][i]: .8f}")
[Confidence interval (from ASE)]
0.09984398 - 0.00148779 <= b'fwhm_G' <= 0.09984398 + 0.00148779
0.05061665 - 0.00077233 <= b't_0_1_1' <= 0.05061665 + 0.00077233
-0.07132916 - 0.00101533 <= b't_0_1_2' <= -0.07132916 + 0.00101533
-0.00474142 - 0.00126955 <= b't_0_1_3' <= -0.00474142 + 0.00126955
0.13670938 - 0.00110190 <= b't_0_1_4' <= 0.13670938 + 0.00110190
0.50172568 - 0.00422955 <= b'tau_1' <= 0.50172568 + 0.00422955
10.00631278 - 0.08691980 <= b'tau_2' <= 10.00631278 + 0.08691980
1002.83212275 - 7.29247489 <= b'tau_3' <= 1002.83212275 + 7.29247489
1.28461622 - 0.49114282 <= b'tau_osc_1' <= 1.28461622 + 0.49114282
0.29797454 - 0.00431719 <= b'period_osc_1' <= 0.29797454 + 0.00431719
However, as you can see, in many case, ASE does not much different from more sophisticated f-test based error estimation.