<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 ZARC model. The equivalent circuit consists one ZARC model, *i.e*., a resistor in parallel with a CPE element.<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>R_ct = 50      # Ohm
R_inf = 10.    # Ohm
phi = 0.8
tau_0 = 1.     # sec<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># generate exact
T = tau_0**phi/R_ct
Z_exact = R_inf + 1./(1./R_ct+T*(1j*2.*pi*freq_vec)**phi)

# random
rng = np.random.seed(121295)
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(Z_exact.real, -Z_exact.imag, linewidth=4, color='black', label='exact')
plt.plot(np.real(Z_exp), -Z_exp.imag, 'o', markersize=8, color='red', label='synth exp')
plt.plot(np.real(Z_exp[0:70:20]), -np.imag(Z_exp[0:70:20]), 's', markersize=8, color="black")
plt.plot(np.real(Z_exp[30]), -np.imag(Z_exp[30]), 's', markersize=8, color="black")

plt.annotate(r'$10^{-4}$', xy=(np.real(Z_exp[0]), -np.imag(Z_exp[0])), 
             xytext=(np.real(Z_exp[0])-15, -np.imag(Z_exp[0])), 
             arrowprops=dict(arrowstyle='-',connectionstyle='arc'))
plt.annotate(r'$10^{-1}$', xy=(np.real(Z_exp[20]), -np.imag(Z_exp[20])), 
             xytext=(np.real(Z_exp[20])-5, 10-np.imag(Z_exp[20])), 
             arrowprops=dict(arrowstyle='-',connectionstyle='arc'))
plt.annotate(r'$1$', xy=(np.real(Z_exp[30]), -np.imag(Z_exp[30])), 
             xytext=(np.real(Z_exp[30]), 8-np.imag(Z_exp[30])), 
             arrowprops=dict(arrowstyle='-',connectionstyle='arc'))
plt.annotate(r'$10$', xy=(np.real(Z_exp[40]), -np.imag(Z_exp[40])), 
             xytext=(np.real(Z_exp[40]), 8-np.imag(Z_exp[40])), 
             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, -np.imag(Z_exp[60])), 
             arrowprops=dict(arrowstyle='-',connectionstyle='arc'))
plt.legend(frameon=False, fontsize = 15)
plt.axis('scaled')

plt.xlim(5, 70)
plt.ylim(-2, 32)
plt.xticks(range(5, 70, 10))
plt.yticks(range(0, 40, 10))
plt.xlabel(r'$Z_{\rm re}/\Omega$', fontsize = 20)
plt.ylabel(r'$-Z_{\rm im}/\Omega$', fontsize = 20)
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
8.07576e-01     19.998508     99.999899 
7.67919e-01     19.891836     99.992062 
7.61233e-01     19.635013     99.973366 
7.55739e-01     19.090800     99.933660 
7.50780e-01     17.975445     99.851742 
7.49009e-01     15.531850     99.668509 
7.56260e-01     15.741292     99.674711 
7.63912e-01     15.639933     99.648666 
7.63957e-01     15.673127     99.640903 
7.63984e-01     15.702439     99.613967 
7.64011e-01     15.732748     99.558007 
7.64043e-01     15.770266     99.448787 
7.64085e-01     15.822560     99.236093 
7.64144e-01     15.901025     98.819714 
7.64232e-01     16.024980     98.002129 
7.64360e-01     16.227507     96.393159 
7.64548e-01     16.567962     93.238301 
7.67402e-01     23.043777     25.880260 
7.67433e-01     23.113920     25.140417 
7.68609e-01     22.301051     25.090725 
7.68236e-01     20.485546     25.856902 
7.67165e-01     19.104197     26.477246 
7.66305e-01     18.776957     26.565019 
7.66305e-01     18.8483[...]<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 real parts of impedance for both Bayesian regression and the synthetic experiment<jupyter_code>band = np.sqrt(cov_Z_re)
plt.fill_between(freq_vec, mu_Z_re-3*band, mu_Z_re+3*band,  facecolor='lightgrey')
plt.semilogx(freq_vec, mu_Z_re, 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(5, 65)
plt.xscale('log')
plt.yticks(range(5, 70, 10))
plt.xlabel(r'$f/{\rm Hz}$', fontsize=20)
plt.ylabel(r'$Z_{\rm re}/\Omega$', fontsize=20)
plt.legend(frameon=False, fontsize = 15)
plt.show()<jupyter_output><empty_output><jupyter_text>3.2 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.3) 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(-3, 70)
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.4) plot the imaginary parts of impedance for both Hilbert transform and the synthetic experiment<jupyter_code>mu_Z_H_im_agm = omega_vec*mu_L_0 + mu_Z_H_im
band_agm = np.sqrt((omega_vec**2)*cov_L_0 + cov_Z_H_im + sigma_n_re**2)
plt.fill_between(freq_vec, -mu_Z_H_im_agm-3*band_agm, -mu_Z_H_im_agm+3*band_agm,  facecolor='lightgrey')
plt.semilogx(freq_vec, -mu_Z_H_im_agm, 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'$-\left(\omega L_0 + Z_{\rm H, im}\right)/\Omega$', fontsize=20)
plt.legend(frameon=False, fontsize = 15)
plt.show()<jupyter_output><empty_output><jupyter_text>3.5) 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.6) 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><jupyter_text>3.7) plot the difference between imaginary parts of impedance for Hilbert transform and the synthetic experiment<jupyter_code>difference_im = omega_vec*mu_L_0 + mu_Z_H_im - Z_exp.imag 
band = np.sqrt((omega_vec**2)*cov_L_0 + cov_Z_H_im + sigma_n_re**2)
plt.fill_between(freq_vec, -3*band, 3*band,  facecolor='lightgrey')
plt.plot(freq_vec, difference_im, '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(\omega L_0 + Z_{\rm H, im} - Z_{\rm exp, im}\right)/\Omega$', fontsize=20)
plt.show()<jupyter_output><empty_output><jupyter_text>3.8) plot the density distribution of residuals for the imaginary part<jupyter_code>fig = plt.figure(2)
a = sns.kdeplot(difference_im, shade=True, color='grey')
a = sns.rugplot(difference_im, color='black')
a.set_xlabel(r'$\left(\omega L_0 + Z_{\rm H, im} - Z_{\rm exp, im}\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.16812865497076024, 'fraction_numerical': 0.07972136222910217, 'mean_word_length': 5.3373297002724795, 'pattern_counts': {'":': 0, '<': 76, '<?xml version=': 0, '>': 76, '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': '19', 'dir': 'jupyter-structured-clean-dedup', 'id': '11995', 'license': 'MIT', 'n_tokens_mistral': 5542, 'n_tokens_neox': 4855, 'n_words': 859, 'path': 'tutorials/.ipynb_checkpoints/ex1_ZARC-checkpoint.ipynb', 'repo_name': 'frank-ccc/BHT', 'size': '11628'}