fixed typos and inconsistencies

fixed typos
polished description
fixed negative sign issue (from previous version incorrectly updated)

tutorials/GP_DRT.py

@@ -1,4 +1,3 @@
-# -*- coding: utf-8 -*-
 """
 Created on Tue Dec 10 2019
 
@@ -28,7 +27,7 @@ def kernel(xi, xi_prime, sigma_f, ell):
     return (sigma_f**2)*exp(-0.5/(ell**2)*((xi-xi_prime)**2))
 
 
-# the function to be integrated in eq (65) of the main text, omitting the constant part and minus sign!!!
+# the function to be integrated in eq (65) of the main text.
 # $\frac{\displaystyle e^{\Delta\xi_{mn}-\chi}}{1+\left(\displaystyle e^{\Delta\xi_{mn}-\chi}\right)^2} \frac{k(\chi)}{\sigma_f^2}$
 def integrand_L_im(x, delta_xi, sigma_f, ell):
     kernel_part = 0.0
@@ -41,19 +40,25 @@ def integrand_L_im(x, delta_xi, sigma_f, ell):
     return kernel_part*sqr_exp
 
 
-# the function to be integrated in eq (76) of the main text, omitting the constant part
+# the function to be integrated in eq (76) of the main text.
 # $\frac{1}{2} \left(\chi+\Delta\xi_{mn}\right){\rm csch}\left(\chi+\Delta\xi_{mn}\right) \frac{k(\chi)}{\sigma_f^2}$
 def integrand_L2_im(x, xi, xi_prime, sigma_f, ell):
-    f = exp(xi)
-    f_prime = exp(xi_prime)
-    numerator = exp(x-0.5/(ell**2)*(x**2))*(x+xi_prime-xi)
-    if x<0:
-        numerator = exp(x-0.5/(ell**2)*(x**2))*(x+xi_prime-xi)
-        denominator = (-1+((f_prime/f)**2)*exp(2*x))
+    delta_xi = xi_prime-xi
+
+    y = x + delta_xi
+    eps = 1E-8
+    if y<-eps:
+        factor_1 = exp(-0.5/(ell**2)*(x**2))
+        factor_2 = y*exp(y)/(exp(2*y)-1.)
+    elif y>eps:
+        factor_1 = exp(-0.5/(ell**2)*(x**2))
+        factor_2 = y*exp(-y)/(1.-exp(-2*y))
     else:
-        numerator = exp(-x-0.5/(ell**2)*(x**2))*(x+xi_prime-xi)
-        denominator = (-exp(-2*x)+((f_prime/f)**2))
-    return numerator/denominator
+        factor_1 = exp(-0.5/(ell**2)*(x**2))
+        factor_2 = 0.5
+
+    out_val = factor_1*factor_2
+    return out_val
 
 
 # the derivative of the integrand in eq (76) with respect to \ell 
@@ -70,7 +75,6 @@ def integrand_der_ell_L2_im(x, xi, xi_prime, sigma_f, ell):
         denominator = (-exp(-2*x)+((f_prime/f)**2))
     return numerator/denominator
 
-
 # assemble the covariance matrix K as shown in eq (18a), which calculates the kernel distance between $\xi_n$ and $\xi_m$
 def matrix_K(xi_n_vec, xi_m_vec, sigma_f, ell):
     N_n_freqs = xi_n_vec.size
@@ -90,25 +94,28 @@ def matrix_L_im_K(xi_n_vec, xi_m_vec, sigma_f, ell):
         xi_vec = xi_n_vec
         N_freqs = xi_vec.size
         L_im_K = np.zeros([N_freqs, N_freqs])
-
-        for n in range(0, N_freqs):
+
+        delta_xi = log(2*pi)
+        integral, tol = integrate.quad(integrand_L_im, -np.inf, np.inf, epsabs=1E-12, epsrel=1E-12, args=(delta_xi, sigma_f, ell))
+        np.fill_diagonal(L_im_K, -(sigma_f**2)*(integral))
+        for n in range(1, N_freqs):
             delta_xi = xi_vec[n]-xi_vec[0] + log(2*pi)
             integral, tol = integrate.quad(integrand_L_im, -np.inf, np.inf, epsabs=1E-12, epsrel=1E-12, args=(delta_xi, sigma_f, ell))
-            np.fill_diagonal(L_im_K[n:, :], (sigma_f**2)*(integral))
+            np.fill_diagonal(L_im_K[:, n:], -(sigma_f**2)*(integral))
 
             delta_xi = xi_vec[0]-xi_vec[n] + log(2*pi)
             integral, tol = integrate.quad(integrand_L_im, -np.inf, np.inf, epsabs=1E-12, epsrel=1E-12, args=(delta_xi, sigma_f, ell))
-            np.fill_diagonal(L_im_K[:, n:], (sigma_f**2)*(integral))
+            np.fill_diagonal(L_im_K[n:, :], -(sigma_f**2)*(integral))
     else:
         N_n_freqs = xi_n_vec.size
         N_m_freqs = xi_m_vec.size
         L_im_K = np.zeros([N_n_freqs, N_m_freqs])
 
         for n in range(0, N_n_freqs):
             for m in range(0, N_m_freqs):
-                delta_xi = xi_n_vec[n]-xi_m_vec[m] + log(2*pi)
+                delta_xi = xi_m_vec[m] - xi_n_vec[n] + log(2*pi)
                 integral, tol = integrate.quad(integrand_L_im, -np.inf, np.inf, epsabs=1E-12, epsrel=1E-12, args=(delta_xi, sigma_f, ell))
-                L_im_K[n,m] =  (sigma_f**2)*(integral);
+                L_im_K[n,m] = -(sigma_f**2)*(integral)
     return L_im_K
 
 
@@ -124,9 +131,11 @@ def matrix_L2_im_K(xi_n_vec, xi_m_vec, sigma_f, ell):
             xi = xi_vec[n]
             xi_prime = xi_vec[0]
             integral, tol = integrate.quad(integrand_L2_im, -np.inf, np.inf, epsabs=1E-12, epsrel=1E-12, args=(xi, xi_prime, sigma_f, ell))
-
-            np.fill_diagonal(L2_im_K[n:, :], exp(xi_prime-xi)*(sigma_f**2)*integral)
-            np.fill_diagonal(L2_im_K[:, n:], exp(xi_prime-xi)*(sigma_f**2)*integral)
+            if n == 0:
+                np.fill_diagonal(L2_im_K, (sigma_f**2)*integral)
+            else:
+                np.fill_diagonal(L2_im_K[n:, :], (sigma_f**2)*integral)
+                np.fill_diagonal(L2_im_K[:, n:], (sigma_f**2)*integral)
     else:
         N_n_freqs = xi_n_vec.size
         N_m_freqs = xi_m_vec.size
@@ -137,23 +146,15 @@ def matrix_L2_im_K(xi_n_vec, xi_m_vec, sigma_f, ell):
             for m in range(0, N_m_freqs):
                 xi_prime = xi_m_vec[m]
                 integral, tol = integrate.quad(integrand_L2_im, -np.inf, np.inf, epsabs=1E-12, epsrel=1E-12, args=(xi, xi_prime, sigma_f, ell))
-                L2_im_K[n,m] = exp(xi_prime-xi)*(sigma_f**2)*integral
+                L2_im_K[n,m] = (sigma_f**2)*integral
     return L2_im_K
 
+def compute_h_L(xi):
 
-# assemble the matrix corresponding to the derivative of eq (18d) with respect to $\ell$, similar to the 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])
+    h_L = 2.*pi*np.exp(xi)
+    h_L = h_L.reshape(xi.size, 1)
 
-    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
+    return h_L
 
 
 # calculate the negative marginal log-likelihood (NMLL) of eq (31)
@@ -179,13 +180,47 @@ def NMLL_fct(theta, Z_exp, xi_vec):
     # easily calculated as the product of the diagonal element of L
     return 0.5*np.dot(Z_exp.imag,alpha) + np.sum(np.log(np.diag(L)))
 
+def NMLL_fct_L(theta, Z_exp, xi_vec):
+
+    sigma_n = theta[0] 
+    sigma_f = theta[1]  
+    ell = theta[2]
+    sigma_L = theta[3]
+
+    N_freqs = xi_vec.size
+
+    L2_im_K = matrix_L2_im_K(xi_vec, xi_vec, sigma_f, ell)
+    Sigma = (sigma_n**2)*np.eye(N_freqs)
+    h_L = compute_h_L(xi_vec)
+
+    K_im_full = L2_im_K + Sigma + (sigma_L**2)*np.outer(h_L, h_L) 
+    L = np.linalg.cholesky(K_im_full)
+
+    # solve for alpha
+    alpha = np.linalg.solve(L, Z_exp.imag)
+    alpha = np.linalg.solve(L.T, alpha)
+    return 0.5*np.dot(Z_exp.imag,alpha) + np.sum(np.log(np.diag(L)))
+
+# assemble the matrix corresponding to the derivative of eq (18d) with respect to $\ell$, similar to the 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
 
 # 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]
+    ell     = theta[2]
 
     N_freqs = xi_vec.size
     Sigma = (sigma_n**2)*np.eye(N_freqs)                    # $\sigma_n^2 \mathbf I$
@@ -219,3 +254,49 @@ def grad_NMLL_fct(theta, Z_exp, xi_vec):
 
     grad = np.array([d_K_im_full_d_sigma_n, d_K_im_full_d_sigma_f, d_K_im_full_d_ell])
     return grad
+
+def grad_NMLL_fct_L(theta, Z_exp, xi_vec):
+
+    sigma_n = theta[0] 
+    sigma_f = theta[1]  
+    ell     = theta[2]
+    sigma_L = theta[3]
+
+    N_freqs = xi_vec.size
+    L2_im_K = matrix_L2_im_K(xi_vec, xi_vec, sigma_f, ell)
+    Sigma = (sigma_n**2)*np.eye(N_freqs)
+    h_L = compute_h_L(xi_vec)
+
+    K_im_full = L2_im_K + Sigma + (sigma_L**2)*np.outer(h_L, h_L)
+
+    L = np.linalg.cholesky(K_im_full)
+
+    # solve for alpha
+    alpha = np.linalg.solve(L, Z_exp.imag)
+    alpha = np.linalg.solve(L.T, alpha)
+
+    # compute inverse of K_im_full
+    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)
+    der_mat_sigma_L = (2.*sigma_L)*np.outer(h_L, h_L)
+
+    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))    
+
+    d_K_im_full_d_sigma_f = - 0.5*np.dot(alpha.T, np.dot(der_mat_sigma_f, alpha)) \
+        + 0.5*np.trace(np.dot(inv_K_im_full, der_mat_sigma_f))
+
+    d_K_im_full_d_ell     = - 0.5*np.dot(alpha.T, np.dot(der_mat_ell, alpha)) \
+        + 0.5*np.trace(np.dot(inv_K_im_full, der_mat_ell))
+
+    d_K_im_full_d_sigma_L = - 0.5*np.dot(alpha.T, np.dot(der_mat_sigma_L, alpha)) \
+        + 0.5*np.trace(np.dot(inv_K_im_full, der_mat_sigma_L))   
+
+    grad = np.array([d_K_im_full_d_sigma_n, d_K_im_full_d_sigma_f, d_K_im_full_d_ell, d_K_im_full_d_sigma_L])
+
+    return grad