fixes to cov.py, fix to how tests handle active_dims and lengthscales

Author: bwengals

pymc3/gp/cov.py

@@ -2,6 +2,7 @@
 import theano.tensor as tt
 import numpy as np
 from functools import reduce
+from operator import mul, add
 
 __all__ = ['ExpQuad',
            'RatQuad',
@@ -91,9 +92,8 @@ def __array_wrap__(self, result):
 
 class Combination(Covariance):
     def __init__(self, factor_list):
-        input_dim = np.max([factor.input_dim for factor in
-                            filter(lambda x: isinstance(x, Covariance),
-                                   factor_list)])
+        input_dim = max([factor.input_dim for factor in factor_list
+                             if isinstance(factor, Covariance)])
         super(Combination, self).__init__(input_dim=input_dim)
         self.factor_list = []
         for factor in factor_list:
@@ -103,45 +103,36 @@ def __init__(self, factor_list):
                 self.factor_list.append(factor)
 
     def merge_factors(self, X, Xs=None, diag=False):
-        # this function makes sure diag=True is handled properly
         factor_list = []
         for factor in self.factor_list:
-
-            # if factor is a Covariance
+            # make sure diag=True is handled properly
             if isinstance(factor, Covariance):
                 factor_list.append(factor(X, Xs, diag))
-                continue
-
-            # if factor is a numpy array
-            if isinstance(factor, np.ndarray):
-                if np.ndim(factor) == 2:
-                    if diag:
-                        factor_list.append(np.diag(factor))
-                        continue
-
-            # if factor is a theano variable with ndim attribute
-            if isinstance(factor, (tt.TensorConstant,
+            elif isinstance(factor, np.ndarray):
+                if np.ndim(factor) == 2 and diag:
+                    factor_list.append(np.diag(factor))
+                else:
+                    factor_list.append(factor)
+            elif isinstance(factor, (tt.TensorConstant,
                                      tt.TensorVariable,
                                      tt.sharedvar.TensorSharedVariable)):
-                if factor.ndim == 2:
-                    if diag:
-                        factor_list.append(tt.diag(factor))
-                        continue
-
-            # othewise
-            factor_list.append(factor)
-
+                if factor.ndim == 2 and diag:
+                    factor_list.append(tt.diag(factor))
+                else:
+                    factor_list.append(factor)
+            else:
+                factor_list.append(factor)
         return factor_list
 
 
 class Add(Combination):
     def __call__(self, X, Xs=None, diag=False):
-        return reduce((lambda x, y: x + y), self.merge_factors(X, Xs, diag))
+        return reduce(add, self.merge_factors(X, Xs, diag))
 
 
 class Prod(Combination):
     def __call__(self, X, Xs=None, diag=False):
-        return reduce((lambda x, y: x * y), self.merge_factors(X, Xs, diag))
+        return reduce(mul, self.merge_factors(X, Xs, diag))
 
 
 class Constant(Covariance):
@@ -205,7 +196,7 @@ def __init__(self, input_dim, ls=None, ls_inv=None, active_dims=None):
         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")
         elif ls_inv is not None:
-            if isinstance(ls_inv, (np.ndarray, list, tuple)):
+            if isinstance(ls_inv, (list, tuple)):
                 ls = 1.0 / np.asarray(ls_inv)
             else:
                 ls = 1.0 / ls_inv
@@ -229,7 +220,8 @@ def euclidean_dist(self, X, Xs):
         return tt.sqrt(r2 + 1e-12)
 
     def diag(self, X):
-        return tt.ones(tt.stack([X.shape[0], ]))
+        return tt.alloc(1.0, X.shape[0])
+        #return tt.ones(tt.stack([X.shape[0], ]))
 
     def full(self, X, Xs=None):
         raise NotImplementedError
@@ -274,7 +266,7 @@ class RatQuad(Stationary):
        k(x, x') = \left(1 + \frac{(x - x')^2}{2\alpha\ell^2} \right)^{-\alpha}
     """
 
-    def __init__(self, input_dim, alpha, ls, ls_inv, active_dims=None):
+    def __init__(self, input_dim, alpha, ls=None, ls_inv=None, active_dims=None):
         super(RatQuad, self).__init__(input_dim, ls, ls_inv, active_dims)
         self.alpha = alpha
 
@@ -337,7 +329,7 @@ class Cosine(Stationary):
     The Cosine kernel.
 
     .. math::
-       k(x, x') = \mathrm{cos}\left( \frac{||x - x'||}{ \ell^2} \right)
+       k(x, x') = \mathrm{cos}\left( \pi \frac{||x - x'||}{ \ell^2} \right)
     """
 
     def full(self, X, Xs=None):

pymc3/tests/test_gp.py

@@ -7,6 +7,7 @@
 import numpy.testing as npt
 import pytest
 
+
 class TestZeroMean(object):
     def test_value(self):
         X = np.linspace(0, 1, 10)[:, None]
@@ -125,6 +126,17 @@ def test_rightadd_matrix(self):
         Kd = theano.function([], cov(X, diag=True))()
         npt.assert_allclose(np.diag(K), Kd, atol=1e-5)
 
+    def test_leftadd_matrixt(self):
+        X = np.linspace(0, 1, 10)[:, None]
+        M = 2 * tt.ones((10, 10))
+        with Model() as model:
+            cov = M + gp.cov.ExpQuad(1, 0.1)
+        K = theano.function([], cov(X))()
+        npt.assert_allclose(K[0, 1], 2.53940, atol=1e-3)
+        # check diagonal
+        Kd = theano.function([], cov(X, diag=True))()
+        npt.assert_allclose(np.diag(K), Kd, atol=1e-5)
+
     def test_leftprod_matrix(self):
         X = np.linspace(0, 1, 3)[:, None]
         M = np.array([[1, 2, 3], [2, 1, 2], [3, 2, 1]])
@@ -222,7 +234,7 @@ def test_slice1(self):
     def test_slice2(self):
         X = np.linspace(0, 1, 30).reshape(10, 3)
         with Model() as model:
-            cov = gp.cov.ExpQuad(3, [0.1, 0.1], active_dims=[False, True, True])
+            cov = gp.cov.ExpQuad(3, ls=[0.1, 0.1], active_dims=[1,2])
         K = theano.function([], cov(X))()
         npt.assert_allclose(K[0, 1], 0.34295549, atol=1e-3)
         # check diagonal
@@ -232,7 +244,7 @@ def test_slice2(self):
     def test_slice3(self):
         X = np.linspace(0, 1, 30).reshape(10, 3)
         with Model() as model:
-            cov = gp.cov.ExpQuad(3, np.array([0.1, 0.1]), active_dims=[False, True, True])
+            cov = gp.cov.ExpQuad(3, ls=np.array([0.1, 0.1]), active_dims=[1,2])
         K = theano.function([], cov(X))()
         npt.assert_allclose(K[0, 1], 0.34295549, atol=1e-3)
         # check diagonal
@@ -242,7 +254,7 @@ def test_slice3(self):
     def test_diffslice(self):
         X = np.linspace(0, 1, 30).reshape(10, 3)
         with Model() as model:
-            cov = gp.cov.ExpQuad(3, 0.1, [1, 0, 0]) + gp.cov.ExpQuad(3, [0.1, 0.2, 0.3])
+            cov = gp.cov.ExpQuad(3, ls=0.1, active_dims=[1, 0, 0]) + gp.cov.ExpQuad(3, ls=[0.1, 0.2, 0.3])
         K = theano.function([], cov(X))()
         npt.assert_allclose(K[0, 1], 0.683572, atol=1e-3)
         # check diagonal
@@ -303,7 +315,7 @@ class TestRatQuad(object):
     def test_1d(self):
         X = np.linspace(0, 1, 10)[:, None]
         with Model() as model:
-            cov = gp.cov.RatQuad(1, 0.1, 0.5)
+            cov = gp.cov.RatQuad(1, ls=0.1, alpha=0.5)
         K = theano.function([], cov(X))()
         npt.assert_allclose(K[0, 1], 0.66896, atol=1e-3)
         K = theano.function([], cov(X, X))()
@@ -442,7 +454,7 @@ def test_raises(self):
         with pytest.raises(NotImplementedError):
             gp.cov.Gibbs(2, lambda x: x)
         with pytest.raises(NotImplementedError):
-            gp.cov.Gibbs(3, lambda x: x, active_dims=[True, True, False])
+            gp.cov.Gibbs(3, lambda x: x, active_dims=[0,1])
 
 
 class TestHandleArgs(object):