Merge branch 'devel' into mrd

GPy/kern/kern.py

@@ -440,7 +440,7 @@ class kern(parameterised):
         # TODO: better looping
         for i1, i2 in itertools.combinations(range(len(self.parts)), 2):
             p1, p2 = self.parts[i1], self.parts[i2]
-            ipsl1, ipsl2 = self.input_slices[i1], self.input_slices[i2]
+#             ipsl1, ipsl2 = self.input_slices[i1], self.input_slices[i2]
             ps1, ps2 = self.param_slices[i1], self.param_slices[i2]
 
             # white doesn;t combine with anything
@@ -455,9 +455,15 @@ class kern(parameterised):
                 p2.dpsi1_dtheta(dL_dpsi2.sum(1) * p1._psi1 * 2., Z, mu, S, target[ps2])
             # linear X bias
             elif p1.name == 'bias' and p2.name == 'linear':
-                p2.dpsi1_dtheta(dL_dpsi2.sum(1)*p1.variance*2., Z, mu, S, target[ps1])
+                p2.dpsi1_dtheta(dL_dpsi2.sum(1) * p1.variance * 2., Z, mu, S, target[ps2])  # [ps1])
+                psi1 = np.zeros((mu.shape[0], Z.shape[0]))
+                p2.psi1(Z, mu, S, psi1)
+                p1.dpsi1_dtheta(dL_dpsi2.sum(1) * psi1 * 2., Z, mu, S, target[ps1])
             elif p2.name == 'bias' and p1.name == 'linear':
                 p1.dpsi1_dtheta(dL_dpsi2.sum(1) * p2.variance * 2., Z, mu, S, target[ps1])
+                psi1 = np.zeros((mu.shape[0], Z.shape[0]))
+                p1.psi1(Z, mu, S, psi1)
+                p2.dpsi1_dtheta(dL_dpsi2.sum(1) * psi1 * 2., Z, mu, S, target[ps2])
             # rbf X linear
             elif p1.name == 'linear' and p2.name == 'rbf':
                 raise NotImplementedError  # TODO
@@ -497,7 +503,7 @@ class kern(parameterised):
                 raise NotImplementedError, "psi2 cannot be computed for this kernel"
 
 
-        return target
+        return target * 2.
 
     def dpsi2_dmuS(self, dL_dpsi2, Z, mu, S, slices1=None, slices2=None):
         """return shapes are N,M,M,Q"""

GPy/kern/linear.py

@@ -114,7 +114,7 @@ class linear(kernpart):
 
     def psi1(self,Z,mu,S,target):
         """the variance, it does nothing"""
-        self.K(mu,Z,target)
+        self._psi1 = self.K(mu, Z, target)
 
     def dpsi1_dtheta(self,dL_dpsi1,Z,mu,S,target):
         """the variance, it does nothing"""

GPy/models/sparse_GP.py

@@ -242,7 +242,7 @@ class sparse_GP(GP):
         dL_dZ = 2.*self.kern.dK_dX(self.dL_dKmm, self.Z)  # factor of two becase of vertical and horizontal 'stripes' in dKmm_dZ
         if self.has_uncertain_inputs:
             dL_dZ += self.kern.dpsi1_dZ(self.dL_dpsi1,self.Z,self.X, self.X_variance)
-            dL_dZ += 2.*self.kern.dpsi2_dZ(self.dL_dpsi2,self.Z,self.X, self.X_variance) # 'stripes'
+            dL_dZ += self.kern.dpsi2_dZ(self.dL_dpsi2, self.Z, self.X, self.X_variance)
         else:
             dL_dZ += self.kern.dK_dX(self.dL_dpsi1,self.Z,self.X)
         return dL_dZ

GPy/testing/bgplvm_tests.py

@@ -60,7 +60,7 @@ class BGPLVMTests(unittest.TestCase):
 
     #@unittest.skip('psi2 cross terms are NotImplemented for this combination')
     def test_linear_bias_kern(self):
-        N, M, Q, D = 10, 3, 2, 4
+        N, M, Q, D = 30, 5, 4, 30
         X = np.random.rand(N, Q)
         k = GPy.kern.linear(Q) +  GPy.kern.bias(Q) + GPy.kern.white(Q, 0.00001)
         K = k.K(X)

GPy/testing/psi_stat_tests.py

@@ -12,18 +12,17 @@ import itertools
 from GPy.core import model
 
 class PsiStatModel(model):
-    def __init__(self, which, X, X_variance, Z, M, kernel, mu_or_S, dL_=numpy.ones((1, 1))):
+    def __init__(self, which, X, X_variance, Z, M, kernel):
         self.which = which
-        self.dL_ = dL_
         self.X = X
         self.X_variance = X_variance
         self.Z = Z
         self.N, self.Q = X.shape
         self.M, Q = Z.shape
-        self.mu_or_S = mu_or_S
         assert self.Q == Q, "shape missmatch: Z:{!s} X:{!s}".format(Z.shape, X.shape)
         self.kern = kernel
         super(PsiStatModel, self).__init__()
+        self.psi_ = self.kern.__getattribute__(self.which)(self.Z, self.X, self.X_variance)
     def _get_param_names(self):
         Xnames = ["{}_{}_{}".format(what, i, j) for what, i, j in itertools.product(['X', 'X_variance'], range(self.N), range(self.Q))]
         Znames = ["Z_{}_{}".format(i, j) for i, j in itertools.product(range(self.M), range(self.Q))]
@@ -41,13 +40,12 @@ class PsiStatModel(model):
     def log_likelihood(self):
         return self.kern.__getattribute__(self.which)(self.Z, self.X, self.X_variance).sum()
     def _log_likelihood_gradients(self):
-        psi_ = self.kern.__getattribute__(self.which)(self.Z, self.X, self.X_variance)
+        psimu, psiS = self.kern.__getattribute__("d" + self.which + "_dmuS")(numpy.ones_like(self.psi_), self.Z, self.X, self.X_variance)
-        psimu, psiS = self.kern.__getattribute__("d" + self.which + "_dmuS")(numpy.ones_like(psi_), self.Z, self.X, self.X_variance)
         try:
-            psiZ = self.kern.__getattribute__("d" + self.which + "_dZ")(numpy.ones_like(psi_), self.Z, self.X, self.X_variance)
+            psiZ = self.kern.__getattribute__("d" + self.which + "_dZ")(numpy.ones_like(self.psi_), self.Z, self.X, self.X_variance)
         except AttributeError:
             psiZ = numpy.zeros(self.M * self.Q)
-        thetagrad = self.kern.__getattribute__("d" + self.which + "_dtheta")(numpy.ones_like(psi_), self.Z, self.X, self.X_variance).flatten()
+        thetagrad = self.kern.__getattribute__("d" + self.which + "_dtheta")(numpy.ones_like(self.psi_), self.Z, self.X, self.X_variance).flatten()
         return numpy.hstack((psimu.flatten(), psiS.flatten(), psiZ.flatten(), thetagrad))
 
 class Test(unittest.TestCase):
@@ -70,20 +68,55 @@ class Test(unittest.TestCase):
                          M=self.M, kernel=k)
             assert m.checkgrad(), "{} x psi0".format("+".join(map(lambda x: x.name, k.parts)))
 
-    def testPsi1(self):
+#     def testPsi1(self):
-        for k in self.kernels:
+#         for k in self.kernels:
-            m = PsiStatModel('psi0', X=self.X, X_variance=self.X_var, Z=self.Z,
+#             m = PsiStatModel('psi1', X=self.X, X_variance=self.X_var, Z=self.Z,
-                         M=self.M, kernel=k)
+#                      M=self.M, kernel=k)
-            assert m.checkgrad(), "{} x psi1".format("+".join(map(lambda x: x.name, k.parts)))
+#             assert m.checkgrad(), "{} x psi1".format("+".join(map(lambda x: x.name, k.parts)))
 
-    def testPsi2(self):
+    def testPsi2_lin(self):
-        for k in self.kernels:
+        k = self.kernels[0]
-            m = PsiStatModel('psi0', X=self.X, X_variance=self.X_var, Z=self.Z,
+        m = PsiStatModel('psi2', X=self.X, X_variance=self.X_var, Z=self.Z,
+                     M=self.M, kernel=k)
+        assert m.checkgrad(), "{} x psi2".format("+".join(map(lambda x: x.name, k.parts)))
+    def testPsi2_lin_bia(self):
+        k = self.kernels[3]
+        m = PsiStatModel('psi2', X=self.X, X_variance=self.X_var, Z=self.Z,
+                     M=self.M, kernel=k)
+        assert m.checkgrad(), "{} x psi2".format("+".join(map(lambda x: x.name, k.parts)))
+    def testPsi2_rbf(self):
+        k = self.kernels[1]
+        m = PsiStatModel('psi2', X=self.X, X_variance=self.X_var, Z=self.Z,
+                     M=self.M, kernel=k)
+        assert m.checkgrad(), "{} x psi2".format("+".join(map(lambda x: x.name, k.parts)))
+    def testPsi2_rbf_bia(self):
+        k = self.kernels[-1]
+        m = PsiStatModel('psi2', X=self.X, X_variance=self.X_var, Z=self.Z,
+                     M=self.M, kernel=k)
+        assert m.checkgrad(), "{} x psi2".format("+".join(map(lambda x: x.name, k.parts)))
+    def testPsi2_bia(self):
+        k = self.kernels[2]
+        m = PsiStatModel('psi2', X=self.X, X_variance=self.X_var, Z=self.Z,
                      M=self.M, kernel=k)
         assert m.checkgrad(), "{} x psi2".format("+".join(map(lambda x: x.name, k.parts)))
 
 
 if __name__ == "__main__":
+    import sys
+    interactive = 'i' in sys.argv
+    if interactive:
+        N, M, Q, D = 30, 5, 4, 30
+        X = numpy.random.rand(N, Q)
+        k = GPy.kern.linear(Q) + GPy.kern.bias(Q) + GPy.kern.white(Q, 0.00001)
+        K = k.K(X)
+        Y = numpy.random.multivariate_normal(numpy.zeros(N), K, D).T
+        Y -= Y.mean(axis=0)
+        k = GPy.kern.linear(Q) + GPy.kern.bias(Q) + GPy.kern.white(Q, 0.00001)
+        m = GPy.models.Bayesian_GPLVM(Y, Q, kernel=k, M=M)
+        m.ensure_default_constraints()
+        m.randomize()
+#         self.assertTrue(m.checkgrad())
+
         Q = 5
         N = 50
         M = 10
@@ -92,11 +125,28 @@ if __name__ == "__main__":
         X_var = .5 * numpy.ones_like(X) + .4 * numpy.clip(numpy.random.randn(*X.shape), 0, 1)
         Z = numpy.random.permutation(X)[:M]
         Y = X.dot(numpy.random.randn(Q, D))
-    kernel = GPy.kern.linear(Q)  # GPy.kern.bias(Q)  # GPy.kern.linear(Q) + GPy.kern.rbf(Q)
+        kernel = GPy.kern.bias(Q)
-    m0 = PsiStatModel('psi0', X=X, X_variance=X_var, Z=Z,
-                     M=M, kernel=kernel, mu_or_S=0, dL_=numpy.ones((1)))
-    m1 = PsiStatModel('psi0', X=X, X_variance=X_var, Z=Z,
-                     M=M, kernel=kernel, mu_or_S=0, dL_=numpy.ones((1)))
-    m2 = PsiStatModel('psi2', X=X, X_variance=X_var, Z=Z,
-                     M=M, kernel=kernel, mu_or_S=0, dL_=numpy.ones((1, 1, 1)))
 
+        kernels = [GPy.kern.linear(Q), GPy.kern.rbf(Q), GPy.kern.bias(Q),
+               GPy.kern.linear(Q) + GPy.kern.bias(Q),
+               GPy.kern.rbf(Q) + GPy.kern.bias(Q)]
+
+#         for k in kernels:
+#             m = PsiStatModel('psi1', X=X, X_variance=X_var, Z=Z,
+#                      M=M, kernel=k)
+#             assert m.checkgrad(), "{} x psi1".format("+".join(map(lambda x: x.name, k.parts)))
+#
+#         m0 = PsiStatModel('psi0', X=X, X_variance=X_var, Z=Z,
+#                          M=M, kernel=GPy.kern.linear(Q))
+#         m1 = PsiStatModel('psi1', X=X, X_variance=X_var, Z=Z,
+#                          M=M, kernel=kernel)
+#         m1 = PsiStatModel('psi1', X=X, X_variance=X_var, Z=Z,
+#                          M=M, kernel=kernel)
+        m2 = PsiStatModel('psi2', X=X, X_variance=X_var, Z=Z,
+                         M=M, kernel=GPy.kern.rbf(Q))
+        m3 = PsiStatModel('psi2', X=X, X_variance=X_var, Z=Z,
+                         M=M, kernel=GPy.kern.linear(Q) + GPy.kern.bias(Q))
+        m4 = PsiStatModel('psi2', X=X, X_variance=X_var, Z=Z,
+                         M=M, kernel=GPy.kern.rbf(Q) + GPy.kern.bias(Q))
+    else:
+        unittest.main()