Circular kernel

RELEASE-NOTES.md

@@ -11,7 +11,7 @@
 
 ### New features
 - `sample_posterior_predictive_w` can now feed on `xarray.Dataset` - e.g. from `InferenceData.posterior`. (see [#4042](https://github.com/pymc-devs/pymc3/pull/4042))
-
+- added `pymc3.gp.cov.Circular` see [#4082](https://github.com/pymc-devs/pymc3/pull/4082), kernel for circular domains (e.g. unit circle)
 
 ## PyMC3 3.9.3 (11 August 2020)
 

pymc3/gp/cov.py

@@ -299,6 +299,73 @@ def full(self, X, Xs=None):
             return tt.alloc(0.0, X.shape[0], Xs.shape[0])
 
 
+class Circular(Covariance):
+    R"""
+    Circular Kernel.
+
+    .. math::
+
+        k_g(x, y) = W_\pi(\operatorname{dist}_{\mathit{geo}}(x, y)),
+
+    with
+
+    .. math::
+
+        W_c = \left(1 + \tau \frac{t}{c}\right)\left(1-\frac{t}{c}\right)
+
+    where :math:`c` is maximum value for :math:`t` and :math:`\tau\ge 4`.
+    The larger :math:`\tau` is the less correlated are neighboring points.
+    See [1]_ for more explanations and use cases.
+
+    Parameters
+    ----------
+    bound : scalar
+        defines the circular interval :math:`[0, \mathit{bound})`
+    tau : scalar
+        :math:`\tau\ge 4` defines correlation strength, the larger,
+        the smaller correlation is. Minimum value is :math:`4`
+
+    References
+    ----------
+    .. [1] Espéran Padonou, O Roustant, "Polar Gaussian Processes for Predicting on Circular Domains"
+    https://hal.archives-ouvertes.fr/hal-01119942v1/document
+    """
+
+    def __init__(self, input_dim, bound, ls=None, ls_inv=None, tau=4, active_dims=None):
+        super().__init__(input_dim, active_dims)
+        if (ls is None and ls_inv is None) or (ls is not None and ls_inv is not None):
+            raise ValueError("Only one of 'ls' or 'ls_inv' must be provided")
+        if len(self.active_dims) != 1:
+            raise ValueError("Only 1 dimension is supported for Circular kernel")
+        elif ls_inv is not None:
+            if isinstance(ls_inv, (list, tuple)):
+                ls = 1.0 / np.asarray(ls_inv)
+            else:
+                ls = 1.0 / ls_inv
+        self.ls = tt.as_tensor_variable(ls)
+        self.c = tt.as_tensor_variable(bound/2) / self.ls
+        self.tau = tau
+
+    def dist(self, X, Xs):
+        X = tt.mul(X, 1.0 / self.ls)
+        if Xs is None:
+            Xs = tt.transpose(X)
+        else:
+            Xs = tt.mul(Xs, 1.0 / self.ls)
+            Xs = tt.transpose(Xs)
+        return tt.abs_((X - Xs + self.c) % (self.c * 2) - self.c)
+
+    def weinland(self, t):
+        return (1 + self.tau * t / self.c) * tt.clip(1 - t / self.c, 0, np.inf) ** self.tau
+
+    def full(self, X, Xs=None):
+        X, Xs = self._slice(X, Xs)
+        return self.weinland(self.dist(X, Xs))
+
+    def diag(self, X):
+        return tt.alloc(1.0, X.shape[0])
+
+
 class Stationary(Covariance):
     r"""
     Base class for stationary kernels/covariance functions.

pymc3/tests/test_gp.py

@@ -1177,3 +1177,17 @@ def test_plot_gp_dist_warn_nan(self):
             )
         plt.close()
         pass
+
+
+class TestCircular:
+    def test_1d(self):
+        X = np.linspace(0, 1, 10)[:, None]
+        with pm.Model():
+            cov = pm.gp.cov.Circular(1, 1, ls=1)
+        K = theano.function([], cov(X))()
+        npt.assert_allclose(K[0, 1], 0.691239, atol=1e-3)
+        K = theano.function([], cov(X, X))()
+        npt.assert_allclose(K[0, 1], 0.691239, atol=1e-3)
+        # check diagonal
+        Kd = theano.function([], cov(X, diag=True))()
+        npt.assert_allclose(np.diag(K), Kd, atol=1e-5)