From 1345e06c572500d9e4dcb289c6de19e93a9bca82 Mon Sep 17 00:00:00 2001 From: ciuccislab Date: Wed, 11 Dec 2019 17:07:38 +0800 Subject: [PATCH] update comments --- GP_DRT.py | 81 ++++++++++++++++++++++++++++++------------------------- 1 file changed, 45 insertions(+), 36 deletions(-) diff --git a/GP_DRT.py b/GP_DRT.py index f7320f8..38b5244 100644 --- a/GP_DRT.py +++ b/GP_DRT.py @@ -4,14 +4,16 @@ @author: Jiapeng Liu, Francesco Ciucci (francesco.ciucci@ust.hk) """ -################################################################################ -# This python file includes all necessary equations for GP-DRT implemented in # -# the paper "Liu, J., & Ciucci, F. (2019). The Gaussian process distribution # -# of relaxation times: A machine learning tool for the analysis and prediction # -# of electrochemical impedance spectroscopy data. Electrochimica Acta, 135316."# -# If you find it is useful, please cite the paper and feel free to use and mod-# -# fy the code. # -################################################################################ + +###################################################################################### +# This python file includes all necessary functions for GP-DRT model implemented in # +# the paper "Liu, J., & Ciucci, F. (2019). The Gaussian process distribution of # +# relaxation times: A machine learning tool for the analysis and prediction of # +# electrochemical impedance spectroscopy data. Electrochimica Acta, 135316." # +# # +# If you find it is useful, please feel free to use it, modify it, and develop based # +# on this code. Please cite the paper if you utlize this code for academic useage. # +###################################################################################### # imports from math import exp @@ -104,7 +106,7 @@ def matrix_L_im_K(xi_n_vec, xi_m_vec, sigma_f, ell): return L_im_K -# assemble the matrix of eq (18b), added the term of $\frac{1}{\sigma_f^2}$ and factor $2\pi$ before $e^{\Delta\xi_{mn}-\chi}$ +# assemble the matrix of eq (18d), added the term of $\frac{1}{\sigma_f^2}$ and factor $2\pi$ before $e^{\Delta\xi_{mn}-\chi}$ # only considering the case that $\xi_n$ and $\xi_m$ are identical, i.e., the matrice are symmetry square def matrix_L2_im_K(xi_n_vec, xi_m_vec, sigma_f, ell): @@ -123,6 +125,24 @@ def matrix_L2_im_K(xi_n_vec, xi_m_vec, sigma_f, ell): return L2_im_K + +# assemble the matrix corresponding to derivative of eq (18d) with respect to $\ell$, similar as above implementation +def der_ell_matrix_L2_im_K(xi_vec, sigma_f, ell): + + N_freqs = xi_vec.size + der_ell_L2_im_K = np.zeros([N_freqs, N_freqs]) + + for n in range(0, N_freqs): + xi = xi_vec[n] + xi_prime = xi_vec[0] + integral, tol = integrate.quad(integrand_der_ell_L2_im, -np.inf, np.inf, epsabs=1E-12, epsrel=1E-12, args=(xi, xi_prime, sigma_f, ell)) + + np.fill_diagonal(der_ell_L2_im_K[n:, :], exp(xi_prime-xi)*(sigma_f**2)/(ell**3)*integral) + np.fill_diagonal(der_ell_L2_im_K[:, n:], exp(xi_prime-xi)*(sigma_f**2)/(ell**3)*integral) + + return der_ell_L2_im_K + + # calculate the negtive marginal log-likelihood (NMLL) of eq (31) def NMLL_fct(theta, Z_exp, xi_vec): # load the initial value for parameters needed to optimize @@ -140,39 +160,27 @@ def NMLL_fct(theta, Z_exp, xi_vec): # solve for alpha alpha = np.linalg.solve(L, Z_exp.imag) alpha = np.linalg.solve(L.T, alpha) # $(\mathcal L^2_{\rm im} \mathbf K + \sigma_n^2 \mathbf I)^{-1} \mathbf Z^{\rm exp}_{\rm im}$ + + # return the final result of $L(\bm \theta)$, i.e., eq (32). Note that $\frac{N}{2} \log 2\pi$ + # is not included as it is a constant. The determinant of $\mathbf K_{\rm im}^{\rm full}$ is + # easily calcualted as product of diagnal element of L return 0.5*np.dot(Z_exp.imag,alpha) + np.sum(np.log(np.diag(L))) -def der_ell_matrix_L2_im_K(xi_vec, sigma_f, ell): - - N_freqs = xi_vec.size - der_ell_L2_im_K = np.zeros([N_freqs, N_freqs]) - - for n in range(0, N_freqs): - xi = xi_vec[n] - xi_prime = xi_vec[0] - integral, tol = integrate.quad(integrand_der_ell_L2_im, -np.inf, np.inf, epsabs=1E-12, epsrel=1E-12, args=(xi, xi_prime, sigma_f, ell)) - - np.fill_diagonal(der_ell_L2_im_K[n:, :], exp(xi_prime-xi)*(sigma_f**2)/(ell**3)*integral) - np.fill_diagonal(der_ell_L2_im_K[:, n:], exp(xi_prime-xi)*(sigma_f**2)/(ell**3)*integral) - - return der_ell_L2_im_K - - - +# gradient of the negative marginal log-likelihhod (NMLL) $L(\bm \theta)$ def grad_NMLL_fct(theta, Z_exp, xi_vec): - + # load the initial value for parameters needed to optimize sigma_n = theta[0] sigma_f = theta[1] ell = theta[2] N_freqs = xi_vec.size - Sigma = (sigma_n**2)*np.eye(N_freqs) - - L2_im_K = matrix_L2_im_K(xi_vec, xi_vec, sigma_f, ell) - K_im_full = L2_im_K + Sigma + Sigma = (sigma_n**2)*np.eye(N_freqs) # $\sigma_n^2 \mathbf I$ - L = np.linalg.cholesky(K_im_full) + L2_im_K = matrix_L2_im_K(xi_vec, xi_vec, sigma_f, ell) # $\mathcal L^2_{\rm im} \mathbf K$ + K_im_full = L2_im_K + Sigma # $\mathbf K_{\rm im}^{\rm full} = \mathcal L^2_{\rm im} \mathbf K + \sigma_n^2 \mathbf I$ + L = np.linalg.cholesky(K_im_full) # Cholesky decomposition to get the inverse of $\mathbf K_{\rm im}^{\rm full}$ + # solve for alpha alpha = np.linalg.solve(L, Z_exp.imag) alpha = np.linalg.solve(L.T, alpha) @@ -181,11 +189,12 @@ def grad_NMLL_fct(theta, Z_exp, xi_vec): inv_L = np.linalg.inv(L) inv_K_im_full = np.dot(inv_L.T, inv_L) - # derivative matrices - der_mat_sigma_n = (2.*sigma_n)*np.eye(N_freqs) - der_mat_sigma_f = (2./sigma_f)*L2_im_K - der_mat_ell = der_ell_matrix_L2_im_K(xi_vec, sigma_f, ell) + # calculate the derivative of matrices with respect to parameters including $sigma_n$, $sigma_f$, and $\ell$ + der_mat_sigma_n = (2.*sigma_n)*np.eye(N_freqs) # derivative of $\sigma_n^2 \mathbf I$ to $sigma_n$ + der_mat_sigma_f = (2./sigma_f)*L2_im_K # note that the derivative is 2*sigma_f*L2_im_K/(sigma_f**2) + der_mat_ell = der_ell_matrix_L2_im_K(xi_vec, sigma_f, ell) # derivative w.r.t $\ell$ + # calculate the derivative of $L(\bm \theta)$ w.r.t $sigma_n$, $sigma_f$, and $\ell$ according to eq (78) d_K_im_full_d_sigma_n = - 0.5*np.dot(alpha.T, np.dot(der_mat_sigma_n, alpha)) \ + 0.5*np.trace(np.dot(inv_K_im_full, der_mat_sigma_n))