<jupyter_start><jupyter_text>Computing the Bayesian Hilbert Transform-DRT In this tutorial example, we will show how the developed BHT-DRT method works using a simple two ZARCs model. The equivalent circuit consists of two ZARCs elements.<jupyter_code># import the libraries
import numpy as np
from math import pi, log10
import matplotlib.pyplot as plt
import seaborn as sns

# core library
import Bayes_HT
import importlib
importlib.reload(Bayes_HT)
# plot standards
plt.rc('font', family='serif', size=15)
plt.rc('text', usetex=True)
plt.rc('xtick', labelsize=15)
plt.rc('ytick', labelsize=15)<jupyter_output><empty_output><jupyter_text>1) Define the synthetic impedance experiment $Z_{\rm exp}(\omega)$  1.1) Define the frequency range<jupyter_code>N_freqs = 81
freq_min = 10**-4   # Hz
freq_max = 10**4    # Hz
freq_vec = np.logspace(log10(freq_min), log10(freq_max), num=N_freqs, endpoint=True)
tau_vec = np.logspace(-log10(freq_max), -log10(freq_min), num=N_freqs, endpoint=True)
omega_vec = 2.*pi*freq_vec<jupyter_output><empty_output><jupyter_text>1.2) Define the circuit parameters for the two ZARCs<jupyter_code># define variables for exact circuit 1
R_inf_1 = 10      # Ohm
R_ct_1 = 50       # Ohm
phi_1 = 0.8
tau_1 = 0.1       # sec

# impedance of circuit 1
T_1 = tau_1**phi_1/R_ct_1
Z_exact_1 = R_inf_1+1./(1./R_ct_1+T_1*(1j*2.*pi*freq_vec)**phi_1)

# define variables for exact circuit 2
R_inf_2 = 10      # Ohm
R_ct_2 = 50       # Ohm
phi_2 = 0.8
tau_2 = 1E1       # sec

# impedance of circuit 2
T_2 = tau_2**phi_2/R_ct_2
Z_exact_2 = R_inf_2+1./(1./R_ct_2+T_2*(1j*2.*pi*freq_vec)**phi_2)<jupyter_output><empty_output><jupyter_text>1.3) Generate exact impedance $Z_{\rm exact}(\omega)$ as well as the stochastic experiment $Z_{\rm exp}(\omega)$, here $Z_{\rm exp}(\omega)=Z_{\rm exact}(\omega)+\sigma_n(\varepsilon_{\rm re}+i\varepsilon_{\rm im})$<jupyter_code>Z_exact = Z_exact_1+Z_exact_2

# random
rng = np.random.seed(1212211)
sigma_n_exp = 0.8       # Ohm
Z_exp = Z_exact + sigma_n_exp*(np.random.normal(0, 1, N_freqs)+1j*np.random.normal(0, 1, N_freqs))<jupyter_output><empty_output><jupyter_text>1.4) show the impedance in Nyquist plot<jupyter_code>fig, ax = plt.subplots()
plt.plot(np.real(Z_exact), -np.imag(Z_exact), linewidth=4, color="black", label="exact")
plt.plot(np.real(Z_exp), -np.imag(Z_exp), "o", markersize=8, color="red", label="synth exp")
plt.plot(np.real(Z_exp[10:70:10]), -np.imag(Z_exp[10:70:10]), 's', markersize=8, color="black")

plt.legend(frameon=False, fontsize = 15)
plt.axis('scaled')

plt.xticks(range(10, 150, 10))
plt.yticks(range(0, 60, 10))
plt.gca().set_aspect('equal', adjustable='box')
plt.xlabel(r'$Z_{\rm re}/\Omega$', fontsize = 20)
plt.ylabel(r'$-Z_{\rm im}/\Omega$', fontsize = 20)
plt.annotate(r'$10^{-3}$', xy=(np.real(Z_exp[10]), -np.imag(Z_exp[10])), 
             xytext=(np.real(Z_exp[10])+10, 5-np.imag(Z_exp[10])), 
             arrowprops=dict(arrowstyle="-", connectionstyle="arc"))
plt.annotate(r'$10^{-2}$', xy=(np.real(Z_exp[20]), -np.imag(Z_exp[20])), 
             xytext=(np.real(Z_exp[20])+10, 5-np.imag(Z_exp[20])), 
             arrowprops=dict(arrowstyle="-", connectionstyle="arc"))
plt.annotate(r'$10^{-1}$', xy=(np.real(Z_exp[30]), -np.imag(Z_exp[30])), 
             xytext=(np.real(Z_exp[30])-5, 10-np.imag(Z_exp[30])), 
             arrowprops=dict(arrowstyle="-", connectionstyle="arc"))
plt.annotate(r'$1$', xy=(np.real(Z_exp[40]), -np.imag(Z_exp[40])), 
             xytext=(np.real(Z_exp[40])-5, 10-np.imag(Z_exp[40])), 
             arrowprops=dict(arrowstyle="-", connectionstyle="arc"))
plt.annotate(r'$10$', xy=(np.real(Z_exp[50]), -np.imag(Z_exp[50])), 
             xytext=(np.real(Z_exp[50])-5, 10-np.imag(Z_exp[50])), 
             arrowprops=dict(arrowstyle="-", connectionstyle="arc"))
plt.annotate(r'$10^2$', xy=(np.real(Z_exp[60]), -np.imag(Z_exp[60])), 
             xytext=(np.real(Z_exp[60])-5, 10-np.imag(Z_exp[60])), 
             arrowprops=dict(arrowstyle="-", connectionstyle="arc"))
plt.show()<jupyter_output><empty_output><jupyter_text>2) Calculate the DRT impedance $Z_{\rm DRT}(\omega)$ and the Hilbert transformed impedance $Z_{\rm H}(\omega)$  2.1) optimize the hyperparamters<jupyter_code># set the intial parameters
sigma_n = 1
sigma_beta = 20
sigma_lambda = 100

theta_0 = np.array([sigma_n, sigma_beta, sigma_lambda])
data_real, data_imag, scores = Bayes_HT.HT_est(theta_0, Z_exp, freq_vec, tau_vec)<jupyter_output>sigma_n; sigma_beta; sigma_lambda
7.94346e-01     20.001648     99.999931 
7.47692e-01     20.116695     99.995787 
7.46375e-01     20.357757     99.986730 
7.46979e-01     20.862391     99.965553 
7.47733e-01     22.004108     99.911132 
7.48413e-01     23.376315     99.833287 
7.48424e-01     23.829517     99.794591 
7.48436e-01     23.944283     99.771499 
7.48440e-01     23.993170     99.743990 
7.48446e-01     24.044109     99.691365 
7.48456e-01     24.107922     99.591287 
7.48473e-01     24.198055     99.398419 
7.48503e-01     24.334707     99.022957 
7.48559e-01     24.551895     98.287721 
7.48666e-01     24.908624     96.842662 
7.48873e-01     25.509549     94.006608 
7.49280e-01     26.547787     88.486577 
7.53209e-01     35.522491     37.606194 
7.53424e-01     35.642164     36.872877 
7.61023e-01     35.017626     37.075223 
7.66783e-01     33.715271     37.433105 
7.71521e-01     31.111735     38.091642 
7.69752e-01     26.372727     38.943711 
7.63337e-01     28.5062[...]<jupyter_text>2.2) Calculate the real part of the $Z_{\rm DRT}(\omega)$ and the imaginary part of the $Z_{\rm H}(\omega)$  2.2.1) Bayesian regression to obtain the real part of impedance for both mean and covariance<jupyter_code>mu_Z_re = data_real.get('mu_Z')
cov_Z_re = np.diag(data_real.get('Sigma_Z'))

# the mean and covariance of $R_\infty$
mu_R_inf = data_real.get('mu_gamma')[0]
cov_R_inf = np.diag(data_real.get('Sigma_gamma'))[0]<jupyter_output><empty_output><jupyter_text>2.2.2) Calculate the real part of DRT impedance for both mean and covariance<jupyter_code>mu_Z_DRT_re = data_real.get('mu_Z_DRT')
cov_Z_DRT_re = np.diag(data_real.get('Sigma_Z_DRT'))<jupyter_output><empty_output><jupyter_text>2.2.3) Calculate the imaginary part of HT impedance for both mean and covariance<jupyter_code>mu_Z_H_im = data_real.get('mu_Z_H')
cov_Z_H_im = np.diag(data_real.get('Sigma_Z_H'))<jupyter_output><empty_output><jupyter_text>2.2.4) Estimate the $\sigma_n$<jupyter_code>sigma_n_re = data_real.get('theta')[0]<jupyter_output><empty_output><jupyter_text>2.3) Calculate the imaginary part of the $Z_{\rm DRT}(\omega)$ and the real part of the $Z_{\rm H}(\omega)$<jupyter_code># 2.3.1 Bayesian regression
mu_Z_im = data_imag.get('mu_Z')
cov_Z_im = np.diag(data_imag.get('Sigma_Z'))

# the mean and covariance of the inductance $L_0$
mu_L_0 = data_imag.get('mu_gamma')[0]
cov_L_0 = np.diag(data_imag.get('Sigma_gamma'))[0]

# 2.3.2 DRT part
mu_Z_DRT_im = data_imag.get('mu_Z_DRT')
cov_Z_DRT_im = np.diag(data_imag.get('Sigma_Z_DRT'))

# 2.3.3 HT prediction
mu_Z_H_re = data_imag.get('mu_Z_H')
cov_Z_H_re = np.diag(data_imag.get('Sigma_Z_H'))

# 2.3.4 estimated sigma_n
sigma_n_im = data_imag.get('theta')[0]<jupyter_output><empty_output><jupyter_text>3) Plot the BHT_DRT  3.1) plot the imaginary parts of impedance for both Bayesian regression and the synthetic experiment<jupyter_code>band = np.sqrt(cov_Z_im)
plt.fill_between(freq_vec, -mu_Z_im-3*band, -mu_Z_im+3*band,  facecolor='lightgrey')
plt.semilogx(freq_vec, -mu_Z_im, linewidth=4, color='black', label='mean')
plt.semilogx(freq_vec, -Z_exp.imag, 'o', markersize=8, color='red', label='synth exp')
plt.xlim(1E-4, 1E4)
plt.ylim(-3, 30)
plt.xscale('log')
plt.xlabel(r'$f/{\rm Hz}$', fontsize=20)
plt.ylabel(r'$-Z_{\rm im}/\Omega$', fontsize=20)
plt.legend(frameon=False, fontsize = 15) 
plt.show()<jupyter_output><empty_output><jupyter_text>3.2) plot the real parts of impedance for both Hilbert transform and the synthetic experiment<jupyter_code>mu_Z_H_re_agm = mu_R_inf + mu_Z_H_re 
band_agm = np.sqrt(cov_R_inf + cov_Z_H_re + sigma_n_im**2)
plt.fill_between(freq_vec, mu_Z_H_re_agm-3*band_agm, mu_Z_H_re_agm+3*band_agm,  facecolor='lightgrey')
plt.semilogx(freq_vec, mu_Z_H_re_agm, linewidth=4, color='black', label='mean')
plt.semilogx(freq_vec, Z_exp.real, 'o', markersize=8, color='red', label='synth exp')
plt.xlim(1E-4, 1E4)
plt.ylim(7, 130)
plt.xscale('log')
plt.xlabel(r'$f/{\rm Hz}$', fontsize=20)
plt.ylabel(r'$\left(R_\infty + Z_{\rm H, re}\right)/\Omega$', fontsize=20)
plt.legend(frameon=False, fontsize = 15)
plt.show()<jupyter_output><empty_output><jupyter_text>3.3) plot the difference between real parts of impedance for Hilbert transform and the synthetic experiment<jupyter_code>difference_re = mu_R_inf + mu_Z_H_re - Z_exp.real 
band = np.sqrt(cov_R_inf + cov_Z_H_re + sigma_n_im**2)
plt.fill_between(freq_vec, -3*band, 3*band,  facecolor='lightgrey')
plt.plot(freq_vec, difference_re, 'o', markersize=8, color='red')
plt.xlim(1E-4, 1E4)
plt.ylim(-10, 10)
plt.xscale('log')
plt.xlabel(r'$f/{\rm Hz}$', fontsize=20)
plt.ylabel(r'$\left(R_\infty + Z_{\rm H, re} - Z_{\rm exp, re}\right)/\Omega$', fontsize=20)
plt.show()<jupyter_output><empty_output><jupyter_text>3.4) plot the density distribution of residuals for the real part<jupyter_code>fig = plt.figure(1)
a = sns.kdeplot(difference_re, shade=True, color='grey')
a = sns.rugplot(difference_re, color='black')
a.set_xlabel(r'$\left(R_\infty + Z_{\rm H, re} - Z_{\rm exp, re}\right)/\Omega$',fontsize=20)
a.set_ylabel(r'pdf',fontsize=20)
a.tick_params(labelsize=15)
plt.xlim(-5, 5)
plt.ylim(0, 0.5)
plt.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.16470832900072788, 'fraction_numerical': 0.09410419049599668, 'mean_word_length': 4.868212324588163, 'pattern_counts': {'":': 0, '<': 60, '<?xml version=': 0, '>': 60, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '15', 'dir': 'jupyter-structured-clean-dedup', 'id': '461508', 'license': 'MIT', 'n_tokens_mistral': 4630, 'n_tokens_neox': 4014, 'n_words': 706, 'path': 'tutorials/ex2_two-ZARCs.ipynb', 'repo_name': 'frank-ccc/BHT', 'size': '9617'}