From 838d2e3706ed41e1f4de45142929102390516157 Mon Sep 17 00:00:00 2001
From: justinbois <bois@caltech.edu>
Date: Sat, 6 Mar 2021 15:00:05 -0800
Subject: [PATCH] Fixed bug in periodic kernel and added linear and polynomial
 kernels

---
 bebi103/__init__.py    |  2 +-
 bebi103/gp.py          | 85 +++++++++++++++++++++++++++++++++++++++---
 codemeta.json          |  2 +-
 doc/user_guide/api.rst |  2 +
 4 files changed, 83 insertions(+), 8 deletions(-)

diff --git a/bebi103/__init__.py b/bebi103/__init__.py
index 42c0895..082e1c3 100644
--- a/bebi103/__init__.py
+++ b/bebi103/__init__.py
@@ -47,4 +47,4 @@
 
 __author__ = """Justin Bois"""
 __email__ = "bois@caltech.edu"
-__version__ = "0.1.4"
+__version__ = "0.1.5"
diff --git a/bebi103/gp.py b/bebi103/gp.py
index 260422d..7af433c 100644
--- a/bebi103/gp.py
+++ b/bebi103/gp.py
@@ -199,7 +199,7 @@ def _vector_to_array(x):
     return x
 
 
-def cov_from_kernel(X1, X2, kernel, **kernel_params):
+def cov_from_kernel(X1, X2, kernel, delta=1e-8, **kernel_params):
     """Return covariance matrix for specified kernel."""
     X1 = _vector_to_array(X1)
     X2 = _vector_to_array(X2) if X2 is not None else None
@@ -220,7 +220,7 @@ def cov_from_kernel(X1, X2, kernel, **kernel_params):
             for j in range(m):
                 K[i, j] = kernel(X1[i, :], X2[j, :], **kernel_params)
 
-    return K
+    return K + np.diag(np.ones(len(K))) * delta
 
 
 @numba.njit
@@ -237,6 +237,77 @@ def _se_covariance_matrix_sym(X, alpha, rho):
     return K
 
 
+def linear_kernel(x1, x2, sigma_b=0.0, sigma=1.0):
+    """Linear kernel. k = sigma_b ** 2 + sigma ** 2 * np.dot(x1, x2)
+
+    Parameters
+    ----------
+    x1 : float or array
+        Point in the space of covariates.
+    x2 : float or array
+        Point in the space of covariates.
+    sigma_b : float, default 0.0
+        Offset of linear kernel.
+    sigma : float or array, default 1.0
+        scale of linear kernel. If an array, must be of the same length
+        as x1 and x2.
+
+    Returns
+    -------
+    output : float
+        Value of returned by kernel evaluated at x1, x2.
+
+    """
+    # Make sure input is vectorial
+    try:
+        len_x1 = len(x1)
+    except:
+        x1 = np.array([x1])
+
+    try:
+        len_x2 = len(x2)
+    except:
+        x2 = np.array([x2])
+
+    return sigma_b**2 + sigma**2 * np.dot(x1, x2)
+
+
+def polynomial_kernel(x1, x2, sigma_b=1.0, sigma_p=1.0, d=1):
+    """Polynomial kernel: (sigma_0^2 + sigma_p^2 x1 . x2)^d
+
+    Parameters
+    ----------
+    x1 : float or array
+        Point in the space of covariates.
+    x2 : float or array
+        Point in the space of covariates.
+    sigma_b : float, default 1.0
+        Offset of polynomial kernel.
+    sigma_p : float, default 1.0
+        Scale of polynomial kernel.
+    d : int, default 1
+        Degree of polynomial kernel
+
+    Returns
+    -------
+    output : float
+        Value of returned by kernel evaluated at x1, x2.
+
+    """
+    # Make sure input is vectorial
+    try:
+        len_x1 = len(x1)
+    except:
+        x1 = np.array([x1])
+
+    try:
+        len_x2 = len(x2)
+    except:
+        x2 = np.array([x2])
+
+    return (sigma_b**2 + sigma_p**2 * np.dot(x1, x2))**d
+
+
 def se_kernel(x1, x2, alpha, rho):
     """Squared exponential kernel.
 
@@ -398,7 +469,7 @@ def matern_kernel(x1, x2, alpha, rho, nu):
     return _matern_kernel(x1, x2, alpha, rho, nu)
 
 
-def periodic_kernel(x1, x2, alpha, rho):
+def periodic_kernel(x1, x2, alpha, rho, T):
     """Periodic kernel.
 
     Parameters
@@ -411,6 +482,8 @@ def periodic_kernel(x1, x2, alpha, rho):
         Marginalized standard deviation of the kernel.
     rho : float
         Length scale of the kernel.
+    T : float
+        Period of kernel.
 
     Returns
     -------
@@ -429,7 +502,7 @@ def periodic_kernel(x1, x2, alpha, rho):
     except:
         x2 = np.array([x2])
 
-    return _periodic_kernel(x1, x2, alpha, rho)
+    return _periodic_kernel(x1, x2, alpha, rho, T)
 
 
 @numba.njit
@@ -514,7 +587,7 @@ def _matern_5_kernel(x1, x2, alpha, rho):
 def _periodic_kernel(x1, x2, alpha, rho, T):
     """Perdioic kernel."""
     x_diff = x1 - x2
-    exp_arg = 2.0 / rho ** 2 * sin(np.pi / T * np.sqrt(np.dot(x_diff, x_diff)))
+    exp_arg = 2.0 / rho ** 2 * np.sin(np.pi / T * np.sqrt(np.dot(x_diff, x_diff))) ** 2
 
     return alpha ** 2 * np.exp(-exp_arg)
 
@@ -710,7 +783,7 @@ def _periodic_covariance_matrix(X1, X2, alpha, rho, T):
 
 
 @numba.njit
-def _periodic_covariance_matrix_sym(X, alpha, rho):
+def _periodic_covariance_matrix_sym(X, alpha, rho, T):
     """Return covariance matrix for squared exponential kernel."""
     n = X.shape[0]
     K = np.empty((n, n))
diff --git a/codemeta.json b/codemeta.json
index 846302d..16488fa 100644
--- a/codemeta.json
+++ b/codemeta.json
@@ -6,7 +6,7 @@
   "codeRepository": "https://github.com/justinbois/bebi103",
   "issueTracker": "https://github.com/justinbois/bebi103/issues",
   "license": "https://spdx.org/licenses/BSD-3-Clause",
-  "version": "0.1.4",
+  "version": "0.1.5",
   "author": [
     {
       "@type": "Person",
diff --git a/doc/user_guide/api.rst b/doc/user_guide/api.rst
index 3d48b5a..7fdea1f 100644
--- a/doc/user_guide/api.rst
+++ b/doc/user_guide/api.rst
@@ -88,6 +88,8 @@ Gaussian process utilities
    :toctree: generated/gp
    :nosignatures:
 
+   linear_kernel
+   polynomial_kernel
    se_kernel
    d1_se_kernel
    d2_se_kernel