BenchMark

Tests exp_conv_gau, exp_conv_cauchy, dmp_osc_conv_gau, dmp_osc_conv_cauchy routine

import numpy as np
from scipy.signal import convolve
from TRXASprefitpack import gau_irf, cauchy_irf, exp_conv_gau, exp_conv_cauchy, dmp_osc_conv_gau, dmp_osc_conv_cauchy
import matplotlib.pyplot as plt

def decay(t, k):
    return np.exp(-k*t)*np.heaviside(t, 1)

def dmp_osc(t, k, T, phi):
    return np.exp(-k*t)*np.cos(2*np.pi*t/T+phi)*np.heaviside(t, 1)

exp_conv_gau

Test condition.

1. fwhm = 0.15, tau = 0.05, time range [-1, 1]

2. fwhm = 0.15, tau = 0.15, time range [-1, 1]

3. fwhm = 0.15, tau = 1.5, time range [-10, 10]

fwhm = 0.15
tau_1 = 0.05
tau_2 = 0.15
tau_3 = 1.5

t_1 = np.linspace(-1, 1, 2001)
t_2 = np.linspace(-1, 1, 2001)
t_3 = np.linspace(-10, 10, 20001)

# To test exp_conv_gau routine, calculates convolution numerically
gau_irf_1 = gau_irf(t_1, fwhm)
gau_irf_2 = gau_irf(t_2, fwhm)
gau_irf_3 = gau_irf(t_3, fwhm)

decay_1 = decay(t_1, 1/tau_1)
decay_2 = decay(t_2, 1/tau_2)
decay_3 = decay(t_3, 1/tau_3)

gau_signal_1_ref = 0.001*convolve(gau_irf_1, decay_1, mode='same')
gau_signal_2_ref = 0.001*convolve(gau_irf_2, decay_2, mode='same')
gau_signal_3_ref = 0.001*convolve(gau_irf_3, decay_3, mode='same')

Test -1-

t_1_sample = np.hstack((np.arange(-1, -0.5, 0.1), np.arange(-0.5, 0.5, 0.05), np.arange(0.5, 1.1, 0.1)))
gau_signal_1_tst = exp_conv_gau(t_1_sample, fwhm, 1/tau_1)

plt.plot(t_1, gau_signal_1_ref, label='ref test 1')
plt.plot(t_1_sample, gau_signal_1_tst, 'ro', label='analytic expression')
plt.legend()
plt.show()

%timeit gau_signal_1_ref = 0.001*convolve(gau_irf_1, decay_1, mode='same')

379 µs ± 1.26 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

%timeit gau_signal_1_tst = exp_conv_gau(t_1_sample, fwhm, 1/tau_1)

24 µs ± 106 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)

Analytic implementation matches to numerical one. Moreover it saves much more computation time (about ~18x), since numerial implementation need to compute convolution of 0.001 ps step data to achieve similiar accuracy level of analytic one.

Test -2-

t_2_sample = np.hstack((np.arange(-1, -0.5, 0.1), np.arange(-0.5, 0.5, 0.05), np.arange(0.5, 1.1, 0.1)))
gau_signal_2_tst = exp_conv_gau(t_1_sample, fwhm, 1/tau_2)

plt.plot(t_2, gau_signal_2_ref, label='ref test 2')
plt.plot(t_2_sample, gau_signal_2_tst, 'ro', label='analytic expression')
plt.legend()
plt.show()

Test -3-

t_3_sample = np.hstack((np.arange(-1, 1, 0.1), np.arange(1, 3, 0.2), np.arange(3, 7, 0.4), np.arange(7, 11, 1)))
gau_signal_3_tst = exp_conv_gau(t_3_sample, fwhm, 1/tau_3)

plt.plot(t_3, gau_signal_3_ref, label='ref test 3')
plt.plot(t_3_sample, gau_signal_3_tst, 'ro', label='analytic expression')
plt.legend()
plt.show()

exp_conv_cauchy

Test condition.

1. fwhm = 0.15, tau = 0.05, time range [-1, 1]

2. fwhm = 0.15, tau = 0.15, time range [-1, 1]

3. fwhm = 0.15, tau = 1.5, time range [-10, 10]

fwhm = 0.15
tau_1 = 0.05
tau_2 = 0.15
tau_3 = 1.5

t_1 = np.linspace(-1, 1, 2001)
t_2 = np.linspace(-1, 1, 2001)
t_3 = np.linspace(-10, 10, 20001)

# To test exp_conv_cauchy routine, calculates convolution numerically
cauchy_irf_1 = cauchy_irf(t_1, fwhm)
cauchy_irf_2 = cauchy_irf(t_2, fwhm)
cauchy_irf_3 = cauchy_irf(t_3, fwhm)

cauchy_signal_1_ref = 0.001*convolve(cauchy_irf_1, decay_1, mode='same')
cauchy_signal_2_ref = 0.001*convolve(cauchy_irf_2, decay_2, mode='same')
cauchy_signal_3_ref = 0.001*convolve(cauchy_irf_3, decay_3, mode='same')

Test -1-

t_1_sample = np.hstack((np.arange(-1, -0.5, 0.1), np.arange(-0.5, 0.5, 0.05), np.arange(0.5, 1.1, 0.1)))
cauchy_signal_1_tst = exp_conv_cauchy(t_1_sample, fwhm, 1/tau_1)

plt.plot(t_1, cauchy_signal_1_ref, label='ref test 1')
plt.plot(t_1_sample, cauchy_signal_1_tst, 'ro', label='analytic expression')
plt.legend()
plt.show()

%timeit cauchy_signal_1_ref = 0.001*convolve(cauchy_irf_1, decay_1, mode='same')

338 µs ± 9.16 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

%timeit cauchy_signal_1_tst = exp_conv_cauchy(t_1_sample, fwhm, 1/tau_1)

68.4 µs ± 1.53 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)

Computation time for exp_conv_cauchy routine is about 2-3 times longer than exp_conv_gau routine. Because exp_conv_cauchy routine requires expansive exponential integral function

Test -2-

t_2_sample = np.hstack((np.arange(-1, -0.5, 0.1), np.arange(-0.5, 0.5, 0.05), np.arange(0.5, 1.1, 0.1)))
cauchy_signal_2_tst = exp_conv_cauchy(t_1_sample, fwhm, 1/tau_2)

plt.plot(t_2, cauchy_signal_2_ref, label='ref test 2')
plt.plot(t_2_sample, cauchy_signal_2_tst, 'ro', label='analytic expression')
plt.legend()
plt.show()

Test -3-

t_3_sample = np.hstack((np.arange(-1, 1, 0.1), np.arange(1, 3, 0.2), np.arange(3, 7, 0.4), np.arange(7, 11, 1)))
cauchy_signal_3_tst = exp_conv_cauchy(t_3_sample, fwhm, 1/tau_3)

plt.plot(t_3, cauchy_signal_3_ref, label='ref test 3')
plt.plot(t_3_sample, cauchy_signal_3_tst, 'ro', label='analytic expression')
plt.legend()
plt.show()

dmp_osc_conv_gau

Test condition.

1. fwhm = 0.15, tau = 10, T = 0.1, phase: 0, time range [-50, 50]

2. fwhm = 0.15, tau = 10, T = 0.5, phase: \(\pi/4\), time range [-50, 50]

3. fwhm = 0.15, tau = 10, T = 3, phase: \(\pi/2\), time range [-50, 50]

fwhm = 0.15
tau = 10
T_1 = 0.1
T_2 = 0.5
T_3 = 3
phi_1 = 0
phi_2 = np.pi/4
phi_3 = np.pi/2

t = np.linspace(-50, 50, 1000001)
t_sample = np.hstack((np.arange(-1, -0.5, 0.1), np.arange(-0.5, 0.5, 0.05), np.arange(0.5, 1, 0.1),
                     np.arange(1, 3, 0.2), np.arange(3, 10, 0.4), np.arange(10, 20, 1), np.arange(20, 55, 5)))

# To test exp_conv_gau routine, calculates convolution numerically
gau_irf_0 = gau_irf(t, fwhm)

dmp_osc_1 = dmp_osc(t, 1/tau, T_1, phi_1)
dmp_osc_2 = dmp_osc(t, 1/tau, T_2, phi_2)
dmp_osc_3 = dmp_osc(t, 1/tau, T_3, phi_3)


dmp_osc_signal_1_ref = (t[1]-t[0])*convolve(gau_irf_0, dmp_osc_1, mode='same')
dmp_osc_signal_2_ref = (t[1]-t[0])*convolve(gau_irf_0, dmp_osc_2, mode='same')
dmp_osc_signal_3_ref = (t[1]-t[0])*convolve(gau_irf_0, dmp_osc_3, mode='same')

dmp_osc_signal_1_tst = dmp_osc_conv_gau(t_sample, fwhm, 1/tau, T_1, phi_1)
dmp_osc_signal_2_tst = dmp_osc_conv_gau(t_sample, fwhm, 1/tau, T_2, phi_2)
dmp_osc_signal_3_tst = dmp_osc_conv_gau(t_sample, fwhm, 1/tau, T_3, phi_3)

Test -1-

plt.plot(t, dmp_osc_signal_1_ref, label='ref test 1')
plt.plot(t_sample, dmp_osc_signal_1_tst, 'ro', label='analytic expression')
plt.legend()
plt.xlim(-1, 1)
plt.show()

%timeit dmp_osc_signal_1_ref = (t[1]-t[0])*convolve(gau_irf_0, dmp_osc_1, mode='same')

111 ms ± 1.91 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)

%timeit dmp_osc_signal_1_tst = dmp_osc_conv_gau(t_sample, fwhm, 1/tau, T_1, phi_1)

43.8 µs ± 110 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)

Analytic implementation matches to numerical one. Moreover it saves much more computation time (about 1000~3000x), since numerial implementation need to compute convolution of 0.0001 ps step data to achieve similiar accuracy level of analytic one

Test -2-

plt.plot(t, dmp_osc_signal_2_ref, label='ref test 2')
plt.plot(t_sample, dmp_osc_signal_2_tst, 'ro', label='analytic expression')
plt.legend()
plt.xlim(-1, 1)
plt.show()

plt.plot(t, dmp_osc_signal_2_ref, label='ref test 2')
plt.plot(t_sample, dmp_osc_signal_2_tst, 'ro', label='analytic expression')
plt.legend()
plt.xlim(-1, 10)
plt.show()

Test -3-

plt.plot(t, dmp_osc_signal_3_ref, label='ref test 3')
plt.plot(t_sample, dmp_osc_signal_3_tst, 'ro', label='analytic expression')
plt.legend()
plt.xlim(-1, 2)
plt.show()

plt.plot(t, dmp_osc_signal_3_ref, label='ref test 3')
plt.plot(t_sample, dmp_osc_signal_3_tst, 'ro', label='analytic expression')
plt.legend()
plt.xlim(-1, 20)
plt.show()

dmp_osc_conv_cauchy

Test condition.

1. fwhm = 0.15, tau = 10, T = 0.1, phase: 0, time range [-50, 50]

2. fwhm = 0.15, tau = 10, T = 0.5, phase: \(\pi/4\), time range [-50, 50]

3. fwhm = 0.15, tau = 10, T = 3, phase: \(\pi/2\), time range [-50, 50]

# To test exp_conv_cauchy routine, calculates convolution numerically
cauchy_irf_0 = cauchy_irf(t, fwhm)


dmp_osc_cauchy_signal_1_ref = (t[1]-t[0])*convolve(cauchy_irf_0, dmp_osc_1, mode='same')
dmp_osc_cauchy_signal_2_ref = (t[1]-t[0])*convolve(cauchy_irf_0, dmp_osc_2, mode='same')
dmp_osc_cauchy_signal_3_ref = (t[1]-t[0])*convolve(cauchy_irf_0, dmp_osc_3, mode='same')

dmp_osc_cauchy_signal_1_tst = dmp_osc_conv_cauchy(t_sample, fwhm, 1/tau, T_1, phi_1)
dmp_osc_cauchy_signal_2_tst = dmp_osc_conv_cauchy(t_sample, fwhm, 1/tau, T_2, phi_2)
dmp_osc_cauchy_signal_3_tst = dmp_osc_conv_cauchy(t_sample, fwhm, 1/tau, T_3, phi_3)

Test -1-

plt.plot(t, dmp_osc_cauchy_signal_1_ref, label='ref test 1')
plt.plot(t_sample, dmp_osc_cauchy_signal_1_tst, 'ro', label='analytic expression')
plt.legend()
plt.xlim(-1, 1)
plt.show()

%timeit dmp_osc_cauchy_signal_1_ref = (t[1]-t[0])*convolve(cauchy_irf_0, dmp_osc_1, mode='same')

106 ms ± 1.02 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)

%timeit dmp_osc_cauchy_signal_1_tst = dmp_osc_conv_cauchy(t_sample, fwhm, 1/tau, T_1, phi_1)

274 µs ± 1.52 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

dmp_osc_cauchy routine is 5-6 times slower than dmp_osc_gau routine due to the evaluation of expansive exponential integral function

Test -2-

plt.plot(t, dmp_osc_cauchy_signal_2_ref, label='ref test 2')
plt.plot(t_sample, dmp_osc_cauchy_signal_2_tst, 'ro', label='analytic expression')
plt.legend()
plt.xlim(-1, 1)
plt.show()

plt.plot(t, dmp_osc_cauchy_signal_2_ref, label='ref test 2')
plt.plot(t_sample, dmp_osc_cauchy_signal_2_tst, 'ro', label='analytic expression')
plt.legend()
plt.xlim(-1, 10)
plt.show()

Test -3-

plt.plot(t, dmp_osc_cauchy_signal_3_ref, label='ref test 3')
plt.plot(t_sample, dmp_osc_cauchy_signal_3_tst, 'ro', label='analytic expression')
plt.legend()
plt.xlim(-1, 2)
plt.show()

plt.plot(t, dmp_osc_cauchy_signal_3_ref, label='ref test 3')
plt.plot(t_sample, dmp_osc_cauchy_signal_3_tst, 'ro', label='analytic expression')
plt.legend()
plt.xlim(-1, 20)
plt.show()