# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python 3
#     language: python
#     name: python3
# ---

# # 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.

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

# ## 1) Define the synthetic impedance experiment $Z_{\rm exp}(\omega)$

# ### 1.1) Define the frequency range

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

# ### 1.2) Define the circuit parameters for the two ZARCs

R_ct = 50      # Ohm
R_inf = 10.    # Ohm
phi = 0.8
tau_0 = 1.     # sec

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

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

# ### 1.4) show the impedance in Nyquist plot

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

# ## 2) Calculate the DRT impedance $Z_{\rm DRT}(\omega)$ and the Hilbert transformed impedance $Z_{\rm H}(\omega)$

# ### 2.1) optimize the hyperparamters

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

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

# +
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]
# -

# #### 2.2.2) Calculate the real part of DRT impedance for both mean and covariance

mu_Z_DRT_re = data_real.get('mu_Z_DRT')
cov_Z_DRT_re = np.diag(data_real.get('Sigma_Z_DRT'))

# #### 2.2.3) Calculate the imaginary part of HT impedance for both mean and covariance

mu_Z_H_im = data_real.get('mu_Z_H')
cov_Z_H_im = np.diag(data_real.get('Sigma_Z_H'))

# #### 2.2.4) Estimate the $\sigma_n$

sigma_n_re = data_real.get('theta')[0]

# ### 2.3) Calculate the imaginary part of the $Z_{\rm DRT}(\omega)$ and the real part of the $Z_{\rm H}(\omega)$

# +
# 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]
# -

# ## 3) Plot the BHT_DRT

# ### 3.1) plot the real parts of impedance for both Bayesian regression and the synthetic experiment

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

# ### 3.2 plot the imaginary parts of impedance for both Bayesian regression and the synthetic experiment

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

# ### 3.3) plot the real parts of impedance for both Hilbert transform and the synthetic experiment

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

# ### 3.4) plot the imaginary parts of impedance for both Hilbert transform and the synthetic experiment

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

# ### 3.5) plot the difference between real parts of impedance for Hilbert transform and the synthetic experiment

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

# ### 3.6) plot the density distribution of residuals for the real part

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

# ### 3.7) plot the difference between imaginary parts of impedance for Hilbert transform and the synthetic experiment

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

# ### 3.8) plot the density distribution of residuals for the imaginary part

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

StarcoderdataJupyterScriptsDedupFiltered