[Merge] merge devel

.coveragerc

@@ -2,7 +2,7 @@
 [run]
 branch = True
 source = GPy
-omit = ./GPy/testing/*.py, travis_tests.py, setup.py, ./GPy/__version__.py
+omit = ./GPy/testing/*.py, travis_tests.py, setup.py, ./GPy/__version__.py, ./GPy/plotting/*
 
 [report]
 # Regexes for lines to exclude from consideration

.travis.yml

@@ -30,6 +30,8 @@ install:
 - source install_retry.sh
 - pip install codecov
 - pip install pypandoc
+- pip install git+git://github.com/BRML/climin.git
+- pip install autograd
 - python setup.py develop
 
 script:

GPy/core/__init__.py

@@ -43,4 +43,4 @@ def randomize(self, rand_gen=None, *args, **kwargs):
     
 Model.randomize = randomize
 Param.randomize = randomize
-Parameterized.randomize = randomize
+Parameterized.randomize = randomize

GPy/core/gp.py

@@ -212,36 +212,9 @@ class GP(Model):
                         = N(f*| K_{x*x}(K_{xx} + \Sigma)^{-1}Y, K_{x*x*} - K_{xx*}(K_{xx} + \Sigma)^{-1}K_{xx*}
             \Sigma := \texttt{Likelihood.variance / Approximate likelihood covariance}
         """
-        if kern is None:
-            kern = self.kern
-
-        Kx = kern.K(self._predictive_variable, Xnew)
-        mu = np.dot(Kx.T, self.posterior.woodbury_vector)
-        if len(mu.shape)==1:
-            mu = mu.reshape(-1,1)
-        if full_cov:
-            Kxx = kern.K(Xnew)
-            if self.posterior.woodbury_inv.ndim == 2:
-                var = Kxx - np.dot(Kx.T, np.dot(self.posterior.woodbury_inv, Kx))
-            elif self.posterior.woodbury_inv.ndim == 3: # Missing data
-                var = np.empty((Kxx.shape[0],Kxx.shape[1],self.posterior.woodbury_inv.shape[2]))
-                from ..util.linalg import mdot
-                for i in range(var.shape[2]):
-                    var[:, :, i] = (Kxx - mdot(Kx.T, self.posterior.woodbury_inv[:, :, i], Kx))
-            var = var
-        else:
-            Kxx = kern.Kdiag(Xnew)
-            if self.posterior.woodbury_inv.ndim == 2:
-                var = (Kxx - np.sum(np.dot(self.posterior.woodbury_inv.T, Kx) * Kx, 0))[:,None]
-            elif self.posterior.woodbury_inv.ndim == 3: # Missing data
-                var = np.empty((Kxx.shape[0],self.posterior.woodbury_inv.shape[2]))
-                for i in range(var.shape[1]):
-                    var[:, i] = (Kxx - (np.sum(np.dot(self.posterior.woodbury_inv[:, :, i].T, Kx) * Kx, 0)))
-            var = var
-        #add in the mean function
+        mu, var = self.posterior._raw_predict(kern=self.kern if kern is None else kern, Xnew=Xnew, pred_var=self._predictive_variable, full_cov=full_cov)
         if self.mean_function is not None:
             mu += self.mean_function.f(Xnew)
-
         return mu, var
 
     def predict(self, Xnew, full_cov=False, Y_metadata=None, kern=None, likelihood=None):
@@ -395,9 +368,9 @@ class GP(Model):
             var_jac = compute_cov_inner(self.posterior.woodbury_inv)
         return mean_jac, var_jac
 
-    def predict_wishard_embedding(self, Xnew, kern=None, mean=True, covariance=True):
+    def predict_wishart_embedding(self, Xnew, kern=None, mean=True, covariance=True):
         """
-        Predict the wishard embedding G of the GP. This is the density of the
+        Predict the wishart embedding G of the GP. This is the density of the
         input of the GP defined by the probabilistic function mapping f.
         G = J_mean.T*J_mean + output_dim*J_cov.
 
@@ -425,6 +398,10 @@ class GP(Model):
             G += Sigma
         return G
 
+    def predict_wishard_embedding(self, Xnew, kern=None, mean=True, covariance=True):
+        warnings.warn("Wrong naming, use predict_wishart_embedding instead. Will be removed in future versions!", DeprecationWarning)
+        return self.predict_wishart_embedding(Xnew, kern, mean, covariance)
+
     def predict_magnification(self, Xnew, kern=None, mean=True, covariance=True):
         """
         Predict the magnification factor as

GPy/core/model.py

@@ -45,4 +45,4 @@ class Model(ParamzModel, Priorizable):
         (including the MAP prior), so we return it here. If your model is not
         probabilistic, just return your *negative* gradient here!
         """
-        return -(self._log_likelihood_gradients() + self._log_prior_gradients())
+        return -(self._log_likelihood_gradients() + self._log_prior_gradients())

GPy/core/parameterization/__init__.py

@@ -3,7 +3,7 @@
 
 from .param import Param
 from .parameterized import Parameterized
-from paramz import transformations
+from . import transformations
 
 from paramz.core import lists_and_dicts, index_operations, observable_array, observable
-from paramz import ties_and_remappings, ObsAr
+from paramz import ties_and_remappings, ObsAr

GPy/core/parameterization/param.py

@@ -7,4 +7,4 @@ from paramz.transformations import __fixed__
 import logging, numpy as np
 
 class Param(Param, Priorizable):
-    pass
+    pass

GPy/core/parameterization/parameterized.py

@@ -49,4 +49,4 @@ class Parameterized(Parameterized, Priorizable):
         If you want to operate on all parameters use m[''] to wildcard select all paramters
         and concatenate them. Printing m[''] will result in printing of all parameters in detail.
     """     
-    pass
+    pass

GPy/core/parameterization/transformations.py

@@ -1,4 +1,5 @@
 # Copyright (c) 2014, Max Zwiessele, James Hensman
 # Licensed under the BSD 3-clause license (see LICENSE.txt)
 
-from paramz.transformations import *
+from paramz.transformations import *
+from paramz.transformations import __fixed__

GPy/core/sparse_gp.py

@@ -44,7 +44,7 @@ class SparseGP(GP):
         #pick a sensible inference method
         if inference_method is None:
             if isinstance(likelihood, likelihoods.Gaussian):
-                inference_method = var_dtc.VarDTC(limit=1)
+                inference_method = var_dtc.VarDTC(limit=3)
             else:
                 #inference_method = ??
                 raise NotImplementedError("what to do what to do?")
@@ -113,85 +113,3 @@ class SparseGP(GP):
             self.Z.gradient += self.kern.gradients_X(self.grad_dict['dL_dKnm'].T, self.Z, self.X)
         self._Zgrad = self.Z.gradient.copy()
 
-
-    def _raw_predict(self, Xnew, full_cov=False, kern=None):
-        """
-        Make a prediction for the latent function values.
-
-        For certain inputs we give back a full_cov of shape NxN,
-        if there is missing data, each dimension has its own full_cov of shape NxNxD, and if full_cov is of,
-        we take only the diagonal elements across N.
-
-        For uncertain inputs, the SparseGP bound produces cannot predict the full covariance matrix full_cov for now.
-        The implementation of that will follow. However, for each dimension the
-        covariance changes, so if full_cov is False (standard), we return the variance
-        for each dimension [NxD].
-        """
-
-        if kern is None: kern = self.kern
-
-        if not isinstance(Xnew, VariationalPosterior):
-            # Kx = kern.K(self._predictive_variable, Xnew)
-            # mu = np.dot(Kx.T, self.posterior.woodbury_vector)
-            # if full_cov:
-            #     Kxx = kern.K(Xnew)
-            #     if self.posterior.woodbury_inv.ndim == 2:
-            #         var = Kxx - np.dot(Kx.T, np.dot(self.posterior.woodbury_inv, Kx))
-            #     elif self.posterior.woodbury_inv.ndim == 3:
-            #         var = np.empty((Kxx.shape[0],Kxx.shape[1],self.posterior.woodbury_inv.shape[2]))
-            #         for i in range(var.shape[2]):
-            #             var[:, :, i] = (Kxx - mdot(Kx.T, self.posterior.woodbury_inv[:, :, i], Kx))
-            #     var = var
-            # else:
-            #     Kxx = kern.Kdiag(Xnew)
-            #     if self.posterior.woodbury_inv.ndim == 2:
-            #         var = (Kxx - np.sum(np.dot(self.posterior.woodbury_inv.T, Kx) * Kx, 0))[:,None]
-            #     elif self.posterior.woodbury_inv.ndim == 3:
-            #         var = np.empty((Kxx.shape[0],self.posterior.woodbury_inv.shape[2]))
-            #         for i in range(var.shape[1]):
-            #             var[:, i] = (Kxx - (np.sum(np.dot(self.posterior.woodbury_inv[:, :, i].T, Kx) * Kx, 0)))
-            #     var = var
-            # #add in the mean function
-            # if self.mean_function is not None:
-            #     mu += self.mean_function.f(Xnew)
-            mu, var = super(SparseGP, self)._raw_predict(Xnew, full_cov, kern)
-        else:
-            psi0_star = kern.psi0(self._predictive_variable, Xnew)
-            psi1_star = kern.psi1(self._predictive_variable, Xnew)
-            psi2_star = kern.psi2n(self._predictive_variable, Xnew) 
-            la = self.posterior.woodbury_vector
-            mu = np.dot(psi1_star, la) # TODO: dimensions?
-            N,M,D = psi0_star.shape[0],psi1_star.shape[1], la.shape[1]
-
-            if full_cov:
-                raise NotImplementedError("Full covariance for Sparse GP predicted with uncertain inputs not implemented yet.")
-                var = np.zeros((Xnew.shape[0], la.shape[1], la.shape[1]))
-                di = np.diag_indices(la.shape[1])
-            else:
-                tmp = psi2_star - psi1_star[:,:,None]*psi1_star[:,None,:]
-                var = (tmp.reshape(-1,M).dot(la).reshape(N,M,D)*la[None,:,:]).sum(1) + psi0_star[:,None] 
-                if self.posterior.woodbury_inv.ndim==2:
-                    var += -psi2_star.reshape(N,-1).dot(self.posterior.woodbury_inv.flat)[:,None]
-                else:
-                    var += -psi2_star.reshape(N,-1).dot(self.posterior.woodbury_inv.reshape(-1,D))
-            assert np.all(var>=-1e-5), "The predicted variance goes negative!: "+str(var)
-            var = np.clip(var,1e-15,np.inf)
-
-#             for i in range(Xnew.shape[0]):
-#                 _mu, _var = Xnew.mean.values[[i]], Xnew.variance.values[[i]]
-#                 psi2_star = kern.psi2(self._predictive_variable, NormalPosterior(_mu, _var))
-#                 tmp = (psi2_star[:, :] - psi1_star[[i]].T.dot(psi1_star[[i]]))
-# 
-#                 var_ = mdot(la.T, tmp, la)
-#                 p0 = psi0_star[i]
-#                 t = np.atleast_3d(self.posterior.woodbury_inv)
-#                 t2 = np.trace(t.T.dot(psi2_star), axis1=1, axis2=2)
-# 
-#                 if full_cov:
-#                     var_[di] += p0
-#                     var_[di] += -t2
-#                     var[i] = var_
-#                 else:
-#                     var[i] = np.diag(var_)+p0-t2
-
-        return mu, var

GPy/examples/dimensionality_reduction.py

@@ -459,7 +459,7 @@ def mrd_simulation(optimize=True, verbose=True, plot=True, plot_sim=True, **kw):
     D1, D2, D3, N, num_inducing, Q = 60, 20, 36, 60, 6, 5
     _, _, Ylist = _simulate_sincos(D1, D2, D3, N, num_inducing, plot_sim)
 
-    k = kern.Linear(Q) + kern.White(Q, variance=1e-4)
+    k = kern.Linear(Q, ARD=True) + kern.White(Q, variance=1e-4)
     m = MRD(Ylist, input_dim=Q, num_inducing=num_inducing, kernel=k, initx="PCA_concat", initz='permute', **kw)
 
     m['.*noise'] = [Y.var() / 40. for Y in Ylist]
@@ -479,7 +479,7 @@ def mrd_simulation_missing_data(optimize=True, verbose=True, plot=True, plot_sim
     D1, D2, D3, N, num_inducing, Q = 60, 20, 36, 60, 6, 5
     _, _, Ylist = _simulate_matern(D1, D2, D3, N, num_inducing, plot_sim)
 
-    k = kern.Linear(Q) + kern.White(Q, variance=1e-4)
+    k = kern.Linear(Q, ARD=True) + kern.White(Q, variance=1e-4)
     inanlist = []
 
     for Y in Ylist:

GPy/inference/latent_function_inference/exact_gaussian_inference.py

@@ -1,14 +1,13 @@
 # Copyright (c) 2012-2014, GPy authors (see AUTHORS.txt).
 # Licensed under the BSD 3-clause license (see LICENSE.txt)
 
-from .posterior import Posterior
+from .posterior import PosteriorExact as Posterior
 from ...util.linalg import pdinv, dpotrs, tdot
 from ...util import diag
 import numpy as np
 from . import LatentFunctionInference
 log_2_pi = np.log(2*np.pi)
 
-
 class ExactGaussianInference(LatentFunctionInference):
     """
     An object for inference when the likelihood is Gaussian.

GPy/inference/latent_function_inference/posterior.py

@@ -2,7 +2,8 @@
 # Licensed under the BSD 3-clause license (see LICENSE.txt)
 
 import numpy as np
-from ...util.linalg import pdinv, dpotrs, dpotri, symmetrify, jitchol
+from ...util.linalg import pdinv, dpotrs, dpotri, symmetrify, jitchol, dtrtrs, tdot
+from GPy.core.parameterization.variational import VariationalPosterior
 
 class Posterior(object):
     """
@@ -187,3 +188,85 @@ class Posterior(object):
         if self._K_chol is None:
             self._K_chol = jitchol(self._K)
         return self._K_chol
+
+    def _raw_predict(self, kern, Xnew, pred_var, full_cov=False):
+        woodbury_vector = self.woodbury_vector
+        woodbury_inv = self.woodbury_inv
+        
+        if not isinstance(Xnew, VariationalPosterior):
+            Kx = kern.K(pred_var, Xnew)
+            mu = np.dot(Kx.T, woodbury_vector)
+            if len(mu.shape)==1:
+                mu = mu.reshape(-1,1)
+            if full_cov:
+                Kxx = kern.K(Xnew)
+                if woodbury_inv.ndim == 2:
+                    var = Kxx - np.dot(Kx.T, np.dot(woodbury_inv, Kx))
+                elif woodbury_inv.ndim == 3: # Missing data
+                    var = np.empty((Kxx.shape[0],Kxx.shape[1],woodbury_inv.shape[2]))
+                    from ...util.linalg import mdot
+                    for i in range(var.shape[2]):
+                        var[:, :, i] = (Kxx - mdot(Kx.T, woodbury_inv[:, :, i], Kx))
+                var = var
+            else:
+                Kxx = kern.Kdiag(Xnew)
+                if woodbury_inv.ndim == 2:
+                    var = (Kxx - np.sum(np.dot(woodbury_inv.T, Kx) * Kx, 0))[:,None]
+                elif woodbury_inv.ndim == 3: # Missing data
+                    var = np.empty((Kxx.shape[0],woodbury_inv.shape[2]))
+                    for i in range(var.shape[1]):
+                        var[:, i] = (Kxx - (np.sum(np.dot(woodbury_inv[:, :, i].T, Kx) * Kx, 0)))
+                var = var
+        else:
+            psi0_star = kern.psi0(pred_var, Xnew)
+            psi1_star = kern.psi1(pred_var, Xnew)
+            psi2_star = kern.psi2n(pred_var, Xnew) 
+            la = woodbury_vector
+            mu = np.dot(psi1_star, la) # TODO: dimensions?
+            N,M,D = psi0_star.shape[0],psi1_star.shape[1], la.shape[1]
+
+            if full_cov:
+                raise NotImplementedError("Full covariance for Sparse GP predicted with uncertain inputs not implemented yet.")
+                var = np.zeros((Xnew.shape[0], la.shape[1], la.shape[1]))
+                di = np.diag_indices(la.shape[1])
+            else:
+                tmp = psi2_star - psi1_star[:,:,None]*psi1_star[:,None,:]
+                var = (tmp.reshape(-1,M).dot(la).reshape(N,M,D)*la[None,:,:]).sum(1) + psi0_star[:,None] 
+                if woodbury_inv.ndim==2:
+                    var += -psi2_star.reshape(N,-1).dot(woodbury_inv.flat)[:,None]
+                else:
+                    var += -psi2_star.reshape(N,-1).dot(woodbury_inv.reshape(-1,D))
+        var = np.clip(var,1e-15,np.inf)
+        return mu, var
+    
+class PosteriorExact(Posterior):
+     
+    def _raw_predict(self, kern, Xnew, pred_var, full_cov=False):
+        
+        Kx = kern.K(pred_var, Xnew)
+        mu = np.dot(Kx.T, self.woodbury_vector)
+        if len(mu.shape)==1:
+            mu = mu.reshape(-1,1)
+        if full_cov:
+            Kxx = kern.K(Xnew)
+            if self._woodbury_chol.ndim == 2:
+                tmp = dtrtrs(self._woodbury_chol, Kx)[0]
+                var = Kxx - tdot(tmp.T)
+            elif self._woodbury_chol.ndim == 3: # Missing data
+                var = np.empty((Kxx.shape[0],Kxx.shape[1],self._woodbury_chol.shape[2]))
+                for i in range(var.shape[2]):
+                    tmp = dtrtrs(self._woodbury_chol[:,:,i], Kx)[0]
+                    var[:, :, i] = (Kxx - tdot(tmp.T))
+            var = var
+        else:
+            Kxx = kern.Kdiag(Xnew)
+            if self._woodbury_chol.ndim == 2:
+                tmp = dtrtrs(self._woodbury_chol, Kx)[0]
+                var = (Kxx - np.square(tmp).sum(0))[:,None]
+            elif self._woodbury_chol.ndim == 3: # Missing data
+                var = np.empty((Kxx.shape[0],self._woodbury_chol.shape[2]))
+                for i in range(var.shape[1]):
+                    tmp = dtrtrs(self._woodbury_chol[:,:,i], Kx)[0]
+                    var[:, i] = (Kxx - np.square(tmp).sum(0))
+            var = var
+        return mu, var

GPy/inference/latent_function_inference/var_dtc.py

@@ -22,7 +22,7 @@ class VarDTC(LatentFunctionInference):
 
     """
     const_jitter = 1e-8
-    def __init__(self, limit=1):
+    def __init__(self, limit=3):
         from paramz.caching import Cacher
         self.limit = limit
         self.get_trYYT = Cacher(self._get_trYYT, limit)

GPy/inference/latent_function_inference/var_dtc_parallel.py

@@ -21,7 +21,7 @@ class VarDTC_minibatch(LatentFunctionInference):
 
     """
     const_jitter = 1e-8
-    def __init__(self, batchsize=None, limit=1, mpi_comm=None):
+    def __init__(self, batchsize=None, limit=3, mpi_comm=None):
 
         self.batchsize = batchsize
         self.mpi_comm = mpi_comm

GPy/inference/mcmc/samplers.py

@@ -37,16 +37,14 @@ class Metropolis_Hastings(object):
 
     def sample(self, Ntotal=10000, Nburn=1000, Nthin=10, tune=True, tune_throughout=False, tune_interval=400):
         current = self.model.optimizer_array
-        fcurrent = self.model.log_likelihood() + self.model.log_prior() + \
-                   self.model._log_det_jacobian()
+        fcurrent = self.model.log_likelihood() + self.model.log_prior() 
         accepted = np.zeros(Ntotal,dtype=np.bool)
         for it in range(Ntotal):
-            print("sample %d of %d\r"%(it,Ntotal),end="\t")
+            print("sample %d of %d\r"%(it+1,Ntotal),end="")
             sys.stdout.flush()
             prop = np.random.multivariate_normal(current, self.cov*self.scale*self.scale)
             self.model.optimizer_array = prop
-            fprop = self.model.log_likelihood() + self.model.log_prior() + \
-                    self.model._log_det_jacobian()
+            fprop = self.model.log_likelihood() + self.model.log_prior() 
 
             if fprop>fcurrent:#sample accepted, going 'uphill'
                 accepted[it] = True

GPy/inference/optimization/__init__.py

@@ -1,5 +1,8 @@
-from paramz.optimization import stochastics, Optimizer
+from paramz.optimization import Optimizer
+from . import stochastics
+
 from paramz.optimization import *
 import sys
+
 sys.modules['GPy.inference.optimization.stochastics'] = stochastics
-sys.modules['GPy.inference.optimization.Optimizer'] = Optimizer
+sys.modules['GPy.inference.optimization.Optimizer'] = Optimizer

GPy/inference/optimization/stochastics.py

@@ -0,0 +1,119 @@
+#===============================================================================
+# Copyright (c) 2015, Max Zwiessele
+# All rights reserved.
+# 
+# Redistribution and use in source and binary forms, with or without
+# modification, are permitted provided that the following conditions are met:
+# 
+# * Redistributions of source code must retain the above copyright notice, this
+#   list of conditions and the following disclaimer.
+# 
+# * Redistributions in binary form must reproduce the above copyright notice,
+#   this list of conditions and the following disclaimer in the documentation
+#   and/or other materials provided with the distribution.
+# 
+# * Neither the name of paramax nor the names of its
+#   contributors may be used to endorse or promote products derived from
+#   this software without specific prior written permission.
+# 
+# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
+# DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
+# FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
+# DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
+# SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
+# CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
+# OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+# OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+#===============================================================================
+
+class StochasticStorage(object):
+    '''
+    This is a container for holding the stochastic parameters,
+    such as subset indices or step length and so on.
+
+    self.d has to be a list of lists:
+    [dimension indices, nan indices for those dimensions]
+    so that the minibatches can be used as efficiently as possible.
+    '''
+    def __init__(self, model):
+        """
+        Initialize this stochastic container using the given model
+        """
+
+    def do_stochastics(self):
+        """
+        Update the internal state to the next batch of the stochastic
+        descent algorithm.
+        """
+        pass
+
+    def reset(self):
+        """
+        Reset the state of this stochastics generator.
+        """
+
+class SparseGPMissing(StochasticStorage):
+    def __init__(self, model, batchsize=1):
+        """
+        Here we want to loop over all dimensions everytime.
+        Thus, we can just make sure the loop goes over self.d every
+        time. We will try to get batches which look the same together
+        which speeds up calculations significantly.
+        """
+        import numpy as np
+        self.Y = model.Y_normalized
+        bdict = {}
+        #For N > 1000 array2string default crops
+        opt = np.get_printoptions()
+        np.set_printoptions(threshold=np.inf)
+        for d in range(self.Y.shape[1]):
+            inan = np.isnan(self.Y)[:, d]
+            arr_str = np.array2string(inan, np.inf, 0, True, '', formatter={'bool':lambda x: '1' if x else '0'})
+            try:
+                bdict[arr_str][0].append(d)
+            except:
+                bdict[arr_str] = [[d], ~inan]
+        np.set_printoptions(**opt)
+        self.d = bdict.values()
+
+class SparseGPStochastics(StochasticStorage):
+    """
+    For the sparse gp we need to store the dimension we are in,
+    and the indices corresponding to those
+    """
+    def __init__(self, model, batchsize=1, missing_data=True):
+        self.batchsize = batchsize
+        self.output_dim = model.Y.shape[1]
+        self.Y = model.Y_normalized
+        self.missing_data = missing_data
+        self.reset()
+        self.do_stochastics()
+
+    def do_stochastics(self):
+        import numpy as np
+        if self.batchsize == 1:
+            self.current_dim = (self.current_dim+1)%self.output_dim
+            self.d = [[[self.current_dim], np.isnan(self.Y[:, self.current_dim]) if self.missing_data else None]]
+        else:
+            self.d = np.random.choice(self.output_dim, size=self.batchsize, replace=False)
+            bdict = {}
+            if self.missing_data:
+                opt = np.get_printoptions()
+                np.set_printoptions(threshold=np.inf)
+                for d in self.d:
+                    inan = np.isnan(self.Y[:, d])
+                    arr_str = np.array2string(inan,np.inf, 0,True, '',formatter={'bool':lambda x: '1' if x else '0'})
+                    try:
+                        bdict[arr_str][0].append(d)
+                    except:
+                        bdict[arr_str] = [[d], ~inan]
+                np.set_printoptions(**opt)
+                self.d = bdict.values()
+            else:
+                self.d = [[self.d, None]]
+
+    def reset(self):
+        self.current_dim = -1
+        self.d = None

GPy/installation.cfg

@@ -15,4 +15,4 @@
 
 
 # [plotting]
-# library = matplotlib # plotly
+# library = matplotlib # plotly, none

GPy/kern/__init__.py

@@ -10,7 +10,7 @@ from .src.add import Add
 from .src.prod import Prod
 from .src.rbf import RBF
 from .src.linear import Linear, LinearFull
-from .src.static import Bias, White, Fixed
+from .src.static import Bias, White, Fixed, WhiteHeteroscedastic
 from .src.brownian import Brownian
 from .src.stationary import Exponential, OU, Matern32, Matern52, ExpQuad, RatQuad, Cosine
 from .src.mlp import MLP
@@ -28,4 +28,4 @@ from .src.trunclinear import TruncLinear,TruncLinear_inf
 from .src.splitKern import SplitKern,DEtime
 from .src.splitKern import DEtime as DiffGenomeKern
 from .src.spline import Spline
-from .src.basis_funcs import LogisticBasisFuncKernel, LinearSlopeBasisFuncKernel, BasisFuncKernel, ChangePointBasisFuncKernel, DomainKernel
+from .src.basis_funcs import LogisticBasisFuncKernel, LinearSlopeBasisFuncKernel, BasisFuncKernel, ChangePointBasisFuncKernel, DomainKernel

GPy/kern/src/ODE_t.py

@@ -162,4 +162,4 @@ class ODE_t(Kern):
 
                 self.lengthscale_Yt.gradient =  np.sum(dkYdlent*(-0.5*self.lengthscale_Yt**(-2)) * dL_dK)
 
-                self.ubias.gradient = np.sum(dkdubias * dL_dK) 
+                self.ubias.gradient = np.sum(dkdubias * dL_dK) 

GPy/kern/src/__init__.py

@@ -1 +1 @@
-from . import psi_comp
+from . import psi_comp

GPy/kern/src/add.py

@@ -19,8 +19,8 @@ class Add(CombinationKernel):
             if isinstance(kern, Add):
                 del subkerns[i]
                 for part in kern.parts[::-1]:
-                    kern.unlink_parameter(part)
-                    subkerns.insert(i, part)
+                    #kern.unlink_parameter(part)
+                    subkerns.insert(i, part.copy())
         super(Add, self).__init__(subkerns, name)
         self._exact_psicomp = self._check_exact_psicomp()
 
@@ -37,7 +37,7 @@ class Add(CombinationKernel):
         else:
             return False
 
-    @Cache_this(limit=2, force_kwargs=['which_parts'])
+    @Cache_this(limit=3, force_kwargs=['which_parts'])
     def K(self, X, X2=None, which_parts=None):
         """
         Add all kernels together.
@@ -51,7 +51,7 @@ class Add(CombinationKernel):
             which_parts = [which_parts]
         return reduce(np.add, (p.K(X, X2) for p in which_parts))
 
-    @Cache_this(limit=2, force_kwargs=['which_parts'])
+    @Cache_this(limit=3, force_kwargs=['which_parts'])
     def Kdiag(self, X, which_parts=None):
         if which_parts is None:
             which_parts = self.parts
@@ -98,17 +98,17 @@ class Add(CombinationKernel):
         [target.__iadd__(p.gradients_XX_diag(dL_dKdiag, X)) for p in self.parts]
         return target
 
-    @Cache_this(limit=1, force_kwargs=['which_parts'])
+    @Cache_this(limit=3, force_kwargs=['which_parts'])
     def psi0(self, Z, variational_posterior):
         if not self._exact_psicomp: return Kern.psi0(self,Z,variational_posterior)
         return reduce(np.add, (p.psi0(Z, variational_posterior) for p in self.parts))
 
-    @Cache_this(limit=1, force_kwargs=['which_parts'])
+    @Cache_this(limit=3, force_kwargs=['which_parts'])
     def psi1(self, Z, variational_posterior):
         if not self._exact_psicomp: return Kern.psi1(self,Z,variational_posterior)
         return reduce(np.add, (p.psi1(Z, variational_posterior) for p in self.parts))
 
-    @Cache_this(limit=1, force_kwargs=['which_parts'])
+    @Cache_this(limit=3, force_kwargs=['which_parts'])
     def psi2(self, Z, variational_posterior):
         if not self._exact_psicomp: return Kern.psi2(self,Z,variational_posterior)
         psi2 = reduce(np.add, (p.psi2(Z, variational_posterior) for p in self.parts))
@@ -144,7 +144,7 @@ class Add(CombinationKernel):
                 raise NotImplementedError("psi2 cannot be computed for this kernel")
         return psi2
 
-    @Cache_this(limit=1, force_kwargs=['which_parts'])
+    @Cache_this(limit=3, force_kwargs=['which_parts'])
     def psi2n(self, Z, variational_posterior):
         if not self._exact_psicomp: return Kern.psi2n(self, Z, variational_posterior)
         psi2 = reduce(np.add, (p.psi2n(Z, variational_posterior) for p in self.parts))
@@ -241,16 +241,20 @@ class Add(CombinationKernel):
             [np.add(target_grads[i],grads[i],target_grads[i]) for i in range(len(grads))]
         return target_grads
 
-    def add(self, other):
-        if isinstance(other, Add):
-            other_params = other.parameters[:]
-            for p in other_params:
-                other.unlink_parameter(p)
-            self.link_parameters(*other_params)
-        else:
-            self.link_parameter(other)
-        self.input_dim, self._all_dims_active = self.get_input_dim_active_dims(self.parts)
-        return self
+    #def add(self, other):
+    #    parts = self.parts
+    #    if 0:#isinstance(other, Add):
+    #        #other_params = other.parameters[:]
+    #        for p in other.parts[:]:
+    #            other.unlink_parameter(p)
+    #        parts.extend(other.parts)
+    #        #self.link_parameters(*other_params)
+    #        
+    #    else:
+    #        #self.link_parameter(other)
+    #        parts.append(other)
+    #    #self.input_dim, self._all_dims_active = self.get_input_dim_active_dims(parts)
+    #    return Add([p for p in parts], self.name)
 
     def input_sensitivity(self, summarize=True):
         if summarize:

GPy/kern/src/eq_ode2.py

@@ -64,7 +64,7 @@ class EQ_ODE2(Kern):
         self.W = Param('W', W)
         self.link_parameters(self.lengthscale, self.C, self.B, self.W)
 
-    @Cache_this(limit=2)
+    @Cache_this(limit=3)
     def K(self, X, X2=None):
         #This way is not working, indexes are lost after using k._slice_X
         #index = np.asarray(X, dtype=np.int)

GPy/kern/src/kern.py

@@ -48,11 +48,12 @@ class Kern(Parameterized):
 
         if active_dims is None:
             active_dims = np.arange(input_dim)
-
-        self.active_dims = active_dims
-        self._all_dims_active = np.atleast_1d(active_dims).astype(int)
-
-        assert self._all_dims_active.size == self.input_dim, "input_dim={} does not match len(active_dim)={}, _all_dims_active={}".format(self.input_dim, self._all_dims_active.size, self._all_dims_active)
+        
+        self.active_dims = np.asarray(active_dims, np.int_)
+        
+        self._all_dims_active = np.atleast_1d(self.active_dims).astype(int)
+        
+        assert self.active_dims.size == self.input_dim, "input_dim={} does not match len(active_dim)={}".format(self.input_dim, self._all_dims_active.size)
 
         self._sliced_X = 0
         self.useGPU = self._support_GPU and useGPU
@@ -68,9 +69,12 @@ class Kern(Parameterized):
     def _effective_input_dim(self):
         return np.size(self._all_dims_active)
 
-    @Cache_this(limit=20)
+    @Cache_this(limit=3)
     def _slice_X(self, X):
-        return X[:, self._all_dims_active]
+        try:
+            return X[:, self._all_dims_active].astype('float')
+        except:
+            return X[:, self._all_dims_active]
 
     def K(self, X, X2):
         """
@@ -319,10 +323,20 @@ class CombinationKernel(Kern):
         :param array-like extra_dims: if needed extra dimensions for the combination kernel to work on
         """
         assert all([isinstance(k, Kern) for k in kernels])
-        extra_dims = np.array(extra_dims, dtype=int)
-        input_dim, active_dims = self.get_input_dim_active_dims(kernels, extra_dims)
+        extra_dims = np.asarray(extra_dims, dtype=int)
+        
+        active_dims = reduce(np.union1d, (np.r_[x.active_dims] for x in kernels), np.array([], dtype=int))
+        
+        input_dim = active_dims.size
+        if extra_dims is not None:
+            input_dim += extra_dims.size
+
         # initialize the kernel with the full input_dim
         super(CombinationKernel, self).__init__(input_dim, active_dims, name)
+
+        effective_input_dim = reduce(max, (k._all_dims_active.max() for k in kernels)) + 1
+        self._all_dims_active = np.array(np.concatenate((np.arange(effective_input_dim), extra_dims if extra_dims is not None else [])), dtype=int)
+        
         self.extra_dims = extra_dims
         self.link_parameters(*kernels)
 
@@ -330,16 +344,8 @@ class CombinationKernel(Kern):
     def parts(self):
         return self.parameters
 
-    def get_input_dim_active_dims(self, kernels, extra_dims = None):
-        self.active_dims = reduce(np.union1d, (np.r_[x.active_dims] for x in kernels), np.array([], dtype=int))
-        #_all_dims_active = np.array(np.concatenate((_all_dims_active, extra_dims if extra_dims is not None else [])), dtype=int)
-        input_dim = reduce(max, (k._all_dims_active.max() for k in kernels)) + 1
-
-        if extra_dims is not None:
-            input_dim += extra_dims.size
-
-        _all_dims_active = np.arange(input_dim)
-        return input_dim, _all_dims_active
+    def _set_all_dims_ative(self):
+        self._all_dims_active = np.atleast_1d(self.active_dims).astype(int)        
 
     def input_sensitivity(self, summarize=True):
         """

GPy/kern/src/linear.py

@@ -51,7 +51,7 @@ class Linear(Kern):
         self.link_parameter(self.variances)
         self.psicomp = PSICOMP_Linear()
 
-    @Cache_this(limit=2)
+    @Cache_this(limit=3)
     def K(self, X, X2=None):
         if self.ARD:
             if X2 is None:
@@ -62,7 +62,7 @@ class Linear(Kern):
         else:
             return self._dot_product(X, X2) * self.variances
 
-    @Cache_this(limit=1, ignore_args=(0,))
+    @Cache_this(limit=3, ignore_args=(0,))
     def _dot_product(self, X, X2=None):
         if X2 is None:
             return tdot(X)

GPy/kern/src/mlp.py

@@ -45,7 +45,7 @@ class MLP(Kern):
         self.link_parameters(self.variance, self.weight_variance, self.bias_variance)
 
 
-    @Cache_this(limit=20, ignore_args=())
+    @Cache_this(limit=3, ignore_args=())
     def K(self, X, X2=None):
         if X2 is None:
             X_denom = np.sqrt(self._comp_prod(X)+1.)
@@ -57,7 +57,7 @@ class MLP(Kern):
         XTX = self._comp_prod(X,X2)/X_denom[:,None]/X2_denom[None,:]
         return self.variance*four_over_tau*np.arcsin(XTX)
 
-    @Cache_this(limit=20, ignore_args=())
+    @Cache_this(limit=3, ignore_args=())
     def Kdiag(self, X):
         """Compute the diagonal of the covariance matrix for X."""
         X_prod = self._comp_prod(X)
@@ -88,14 +88,14 @@ class MLP(Kern):
         """Gradient of diagonal of covariance with respect to X"""
         return self._comp_grads_diag(dL_dKdiag, X)[3]
 
-    @Cache_this(limit=50, ignore_args=())
+    @Cache_this(limit=3, ignore_args=())
     def _comp_prod(self, X, X2=None):
         if X2 is None:
             return (np.square(X)*self.weight_variance).sum(axis=1)+self.bias_variance
         else:
             return (X*self.weight_variance).dot(X2.T)+self.bias_variance
     
-    @Cache_this(limit=20, ignore_args=(1,))
+    @Cache_this(limit=3, ignore_args=(1,))
     def _comp_grads(self, dL_dK, X, X2=None):
         var,w,b = self.variance, self.weight_variance, self.bias_variance
         K = self.K(X, X2)
@@ -130,7 +130,7 @@ class MLP(Kern):
             dX2 = common.T.dot(X)*w-((common*XTX).sum(axis=0)/(X2_prod+1.))[:,None]*X2*w
         return dvar, dw, db, dX, dX2
     
-    @Cache_this(limit=20, ignore_args=(1,))
+    @Cache_this(limit=3, ignore_args=(1,))
     def _comp_grads_diag(self, dL_dKdiag, X):
         var,w,b = self.variance, self.weight_variance, self.bias_variance
         K = self.Kdiag(X)

GPy/kern/src/poly.py

@@ -27,7 +27,7 @@ class Poly(Kern):
         _, _, B = self._AB(X, X2)
         return B * self.variance
 
-    @Cache_this(limit=2)
+    @Cache_this(limit=3)
     def _AB(self, X, X2=None):
         if X2 is None:
             dot_prod = np.dot(X, X.T)

GPy/kern/src/prod.py

@@ -39,7 +39,7 @@ class Prod(CombinationKernel):
                     kernels.insert(i, part)
         super(Prod, self).__init__(kernels, name)
 
-    @Cache_this(limit=2, force_kwargs=['which_parts'])
+    @Cache_this(limit=3, force_kwargs=['which_parts'])
     def K(self, X, X2=None, which_parts=None):
         if which_parts is None:
             which_parts = self.parts
@@ -48,7 +48,7 @@ class Prod(CombinationKernel):
             which_parts = [which_parts]
         return reduce(np.multiply, (p.K(X, X2) for p in which_parts))
 
-    @Cache_this(limit=2, force_kwargs=['which_parts'])
+    @Cache_this(limit=3, force_kwargs=['which_parts'])
     def Kdiag(self, X, which_parts=None):
         if which_parts is None:
             which_parts = self.parts

GPy/kern/src/psi_comp/__init__.py

@@ -21,7 +21,7 @@ from .gaussherm import PSICOMP_GH
 from . import rbf_psi_comp, linear_psi_comp, ssrbf_psi_comp, sslinear_psi_comp
 
 class PSICOMP_RBF(PSICOMP):
-    @Cache_this(limit=10, ignore_args=(0,))
+    @Cache_this(limit=3, ignore_args=(0,))
     def psicomputations(self, kern, Z, variational_posterior, return_psi2_n=False):
         variance, lengthscale = kern.variance, kern.lengthscale
         if isinstance(variational_posterior, variational.NormalPosterior):
@@ -31,7 +31,7 @@ class PSICOMP_RBF(PSICOMP):
         else:
             raise ValueError("unknown distriubtion received for psi-statistics")
 
-    @Cache_this(limit=10, ignore_args=(0,2,3,4))
+    @Cache_this(limit=3, ignore_args=(0,2,3,4))
     def psiDerivativecomputations(self, kern, dL_dpsi0, dL_dpsi1, dL_dpsi2, Z, variational_posterior):
         variance, lengthscale = kern.variance, kern.lengthscale
         if isinstance(variational_posterior, variational.NormalPosterior):
@@ -43,7 +43,7 @@ class PSICOMP_RBF(PSICOMP):
 
 class PSICOMP_Linear(PSICOMP):
 
-    @Cache_this(limit=10, ignore_args=(0,))
+    @Cache_this(limit=3, ignore_args=(0,))
     def psicomputations(self, kern, Z, variational_posterior, return_psi2_n=False):
         variances = kern.variances
         if isinstance(variational_posterior, variational.NormalPosterior):
@@ -53,7 +53,7 @@ class PSICOMP_Linear(PSICOMP):
         else:
             raise ValueError("unknown distriubtion received for psi-statistics")
 
-    @Cache_this(limit=10, ignore_args=(0,2,3,4))
+    @Cache_this(limit=3, ignore_args=(0,2,3,4))
     def psiDerivativecomputations(self, kern, dL_dpsi0, dL_dpsi1, dL_dpsi2, Z, variational_posterior):
         variances = kern.variances
         if isinstance(variational_posterior, variational.NormalPosterior):

GPy/kern/src/psi_comp/gaussherm.py

@@ -27,7 +27,7 @@ class PSICOMP_GH(PSICOMP):
     def _setup_observers(self):
         pass
     
-    @Cache_this(limit=10, ignore_args=(0,))
+    @Cache_this(limit=3, ignore_args=(0,))
     def comp_K(self, Z, qX):
         if self.Xs is None or self.Xs.shape != qX.mean.shape:
             from paramz import ObsAr
@@ -38,7 +38,7 @@ class PSICOMP_GH(PSICOMP):
             self.Xs[i] = self.locs[i]*S_sq+mu
         return self.Xs
     
-    @Cache_this(limit=10, ignore_args=(0,))
+    @Cache_this(limit=3, ignore_args=(0,))
     def psicomputations(self, kern, Z, qX, return_psi2_n=False):
         mu, S = qX.mean.values, qX.variance.values
         N,M,Q = mu.shape[0],Z.shape[0],mu.shape[1]
@@ -62,7 +62,7 @@ class PSICOMP_GH(PSICOMP):
                 psi2 += self.weights[i]* tdot(Kfu.T)
         return psi0, psi1, psi2
     
-    @Cache_this(limit=10, ignore_args=(0, 2,3,4))
+    @Cache_this(limit=3, ignore_args=(0, 2,3,4))
     def psiDerivativecomputations(self, kern, dL_dpsi0, dL_dpsi1, dL_dpsi2, Z, qX):
         mu, S = qX.mean.values, qX.variance.values
         if self.cache_K: Xs = self.comp_K(Z, qX)

GPy/kern/src/psi_comp/rbf_psi_comp.py

@@ -132,5 +132,5 @@ def _psi2compDer(dL_dpsi2, variance, lengthscale, Z, mu, S):
 
     return _dL_dvar, _dL_dl, _dL_dZ, _dL_dmu, _dL_dS
 
-_psi1computations = Cacher(__psi1computations, limit=5)
-_psi2computations = Cacher(__psi2computations, limit=5)
+_psi1computations = Cacher(__psi1computations, limit=3)
+_psi2computations = Cacher(__psi2computations, limit=3)

GPy/kern/src/psi_comp/rbf_psi_gpucomp.py

@@ -5,7 +5,6 @@ The module for psi-statistics for RBF kernel
 import numpy as np
 from paramz.caching import Cache_this
 from . import PSICOMP_RBF
-from ....util import gpu_init
 
 gpu_code = """
     // define THREADNUM
@@ -238,8 +237,7 @@ class PSICOMP_RBF_GPU(PSICOMP_RBF):
         self.fall_back = PSICOMP_RBF()
         
         from pycuda.compiler import SourceModule
-        from ....util.gpu_init import initGPU
-        initGPU()
+        import GPy.util.gpu_init
         
         self.GPU_direct = GPU_direct
         self.gpuCache = None
@@ -326,7 +324,7 @@ class PSICOMP_RBF_GPU(PSICOMP_RBF):
         except:
             return self.fall_back.psicomputations(kern, Z, variational_posterior, return_psi2_n)
 
-    @Cache_this(limit=10, ignore_args=(0,))
+    @Cache_this(limit=3, ignore_args=(0,))
     def _psicomputations(self, kern, Z, variational_posterior, return_psi2_n=False):
         """
         Z - MxQ
@@ -371,7 +369,7 @@ class PSICOMP_RBF_GPU(PSICOMP_RBF):
         except:
             return self.fall_back.psiDerivativecomputations(kern, dL_dpsi0, dL_dpsi1, dL_dpsi2, Z, variational_posterior)
 
-    @Cache_this(limit=10, ignore_args=(0,2,3,4))
+    @Cache_this(limit=3, ignore_args=(0,2,3,4))
     def _psiDerivativecomputations(self, kern, dL_dpsi0, dL_dpsi1, dL_dpsi2, Z, variational_posterior):
         # resolve the requirement of dL_dpsi2 to be symmetric
         if len(dL_dpsi2.shape)==2: dL_dpsi2 = (dL_dpsi2+dL_dpsi2.T)/2

GPy/kern/src/psi_comp/ssrbf_psi_comp.py

@@ -88,7 +88,7 @@ try:
         return psi0,psi1,psi2,psi2n
 
     from GPy.util.caching import Cacher
-    psicomputations = Cacher(_psicomputations, limit=1)
+    psicomputations = Cacher(_psicomputations, limit=3)
 
     def psiDerivativecomputations(dL_dpsi0, dL_dpsi1, dL_dpsi2, variance, lengthscale, Z, variational_posterior):
         ARD = (len(lengthscale)!=1)

GPy/kern/src/psi_comp/ssrbf_psi_gpucomp.py

@@ -7,7 +7,6 @@ import numpy as np
 from paramz.caching import Cache_this
 from . import PSICOMP_RBF
 
-
 gpu_code = """
     // define THREADNUM
 
@@ -287,8 +286,7 @@ class PSICOMP_SSRBF_GPU(PSICOMP_RBF):
     def __init__(self, threadnum=128, blocknum=15, GPU_direct=False):
         
         from pycuda.compiler import SourceModule
-        from ....util.gpu_init import initGPU
-        initGPU()
+        import GPy.util.gpu_init
         
         self.GPU_direct = GPU_direct
         self.gpuCache = None
@@ -375,7 +373,7 @@ class PSICOMP_SSRBF_GPU(PSICOMP_RBF):
     def get_dimensions(self, Z, variational_posterior):
         return variational_posterior.mean.shape[0], Z.shape[0], Z.shape[1]
 
-    @Cache_this(limit=1, ignore_args=(0,))
+    @Cache_this(limit=3, ignore_args=(0,))
     def psicomputations(self, kern, Z, variational_posterior, return_psi2_n=False):
         """
         Z - MxQ
@@ -409,7 +407,7 @@ class PSICOMP_SSRBF_GPU(PSICOMP_RBF):
         else:
             return psi0, psi1_gpu.get(), psi2_gpu.get()
 
-    @Cache_this(limit=1, ignore_args=(0,2,3,4))
+    @Cache_this(limit=3, ignore_args=(0,2,3,4))
     def psiDerivativecomputations(self, kern, dL_dpsi0, dL_dpsi1, dL_dpsi2, Z, variational_posterior):
         variance, lengthscale = kern.variance, kern.lengthscale
         from ....util.linalg_gpu import sum_axis

GPy/kern/src/standard_periodic.py

@@ -1,6 +1,5 @@
 # -*- coding: utf-8 -*-
-
-# Copyright (c) 2014, GPy authors (see AUTHORS.txt).
+# Copyright (c) 2015, Alex Grigorevskiy
 # Licensed under the BSD 3-clause license (see LICENSE.txt)
 """
 The standard periodic kernel which mentioned in:
@@ -9,7 +8,7 @@ The standard periodic kernel which mentioned in:
 The MIT Press, 2005.
 
 
-[2] Introduction to Gaussian processes. D. J. C. MacKay. In C. M. Bishop, editor, 
+[2] Introduction to Gaussian processes. D. J. C. MacKay. In C. M. Bishop, editor,
 Neural Networks and Machine Learning, pages 133-165. Springer, 1998.
 """
 
@@ -25,56 +24,56 @@ class StdPeriodic(Kern):
 
     .. math::
 
-       k(x,y) = \theta_1 \exp \left[  - \frac{1}{2} {}\sum_{i=1}^{input\_dim}  
-       \left( \frac{\sin(\frac{\pi}{\lambda_i} (x_i - y_i) )}{l_i} \right)^2 \right] }
+       k(x,y) = \theta_1 \exp \left[  - \frac{1}{2} \sum_{i=1}^{input\_dim}
+       \left( \frac{\sin(\frac{\pi}{T_i} (x_i - y_i) )}{l_i} \right)^2 \right] }
 
     :param input_dim: the number of input dimensions
     :type input_dim: int
     :param variance: the variance :math:`\theta_1` in the formula above
     :type variance: float
-    :param wavelength: the vector of wavelengths :math:`\lambda_i`. If None then 1.0 is assumed.
-    :type wavelength: array or list of the appropriate size (or float if there is only one wavelength parameter)
+    :param period: the vector of periods :math:`\T_i`. If None then 1.0 is assumed.
+    :type period: array or list of the appropriate size (or float if there is only one period parameter)
     :param lengthscale: the vector of lengthscale :math:`\l_i`. If None then 1.0 is assumed.
     :type lengthscale: array or list of the appropriate size (or float if there is only one lengthscale parameter)
-    :param ARD1: Auto Relevance Determination with respect to wavelength. 
-        If equal to "False" one single wavelength parameter :math:`\lambda_i` for 
-        each dimension is assumed, otherwise there is one lengthscale 
+    :param ARD1: Auto Relevance Determination with respect to period.
+        If equal to "False" one single period parameter :math:`\T_i` for
+        each dimension is assumed, otherwise there is one lengthscale
         parameter per dimension.
     :type ARD1: Boolean
-    :param ARD2: Auto Relevance Determination with respect to lengthscale. 
-        If equal to "False" one single wavelength parameter :math:`l_i` for 
-        each dimension is assumed, otherwise there is one lengthscale 
+    :param ARD2: Auto Relevance Determination with respect to lengthscale.
+        If equal to "False" one single lengthscale parameter :math:`l_i` for
+        each dimension is assumed, otherwise there is one lengthscale
         parameter per dimension.
     :type ARD2: Boolean
     :param active_dims: indices of dimensions which are used in the computation of the kernel
-    :type wavelength: array or list of the appropriate size
+    :type active_dims: array or list of the appropriate size
     :param name: Name of the kernel for output
     :type String
     :param useGPU: whether of not use GPU
     :type Boolean
     """
-    
-    def __init__(self, input_dim, variance=1., wavelength=None, lengthscale=None, ARD1=False, ARD2=False, active_dims=None, name='std_periodic',useGPU=False):
+
+    def __init__(self, input_dim, variance=1., period=None, lengthscale=None, ARD1=False, ARD2=False, active_dims=None, name='std_periodic',useGPU=False):
         super(StdPeriodic, self).__init__(input_dim, active_dims, name, useGPU=useGPU)
         self.input_dim = input_dim
-        self.ARD1 = ARD1 # correspond to wavelengths        
+        self.ARD1 = ARD1 # correspond to periods
         self.ARD2 = ARD2 # correspond to lengthscales
-        
+
         self.name = name
-        
+
         if self.ARD1 == False:
-            if wavelength is not None:
-                wavelength = np.asarray(wavelength)
-                assert wavelength.size == 1, "Only one wavelength needed for non-ARD kernel"
+            if period is not None:
+                period = np.asarray(period)
+                assert period.size == 1, "Only one period needed for non-ARD kernel"
             else:
-                wavelength = np.ones(1)
+                period = np.ones(1.0)
         else:
-            if wavelength is not None:
-                wavelength = np.asarray(wavelength)
-                assert wavelength.size == input_dim, "bad number of wavelengths"
+            if period is not None:
+                period = np.asarray(period)
+                assert period.size == input_dim, "bad number of periods"
             else:
-                wavelength = np.ones(input_dim)
-        
+                period = np.ones(input_dim)
+
         if self.ARD2 == False:
             if lengthscale is not None:
                 lengthscale = np.asarray(lengthscale)
@@ -87,33 +86,33 @@ class StdPeriodic(Kern):
                 assert lengthscale.size == input_dim, "bad number of lengthscales"
             else:
                 lengthscale = np.ones(input_dim)
-        
+
         self.variance = Param('variance', variance, Logexp())
         assert self.variance.size==1, "Variance size must be one"
-        self.wavelengths =  Param('wavelengths', wavelength, Logexp())
-        self.lengthscales =  Param('lengthscales', lengthscale, Logexp())
-        
-        self.link_parameters(self.variance,  self.wavelengths, self.lengthscales)
+        self.period =  Param('period', period, Logexp())
+        self.lengthscale =  Param('lengthscale', lengthscale, Logexp())
+
+        self.link_parameters(self.variance,  self.period, self.lengthscale)
 
     def parameters_changed(self):
         """
-        This functions deals as a callback for each optimization iteration. 
+        This functions deals as a callback for each optimization iteration.
         If one optimization step was successfull and the parameters
-        this callback function will be called to be able to update any 
+        this callback function will be called to be able to update any
         precomputations for the kernel.
         """
-        
+
         pass
-        
-        
+
+
     def K(self, X, X2=None):
         """Compute the covariance matrix between X and X2."""
-        if X2 is None: 
+        if X2 is None:
             X2 = X
-            
-        base = np.pi * (X[:, None, :] - X2[None, :, :]) / self.wavelengths
-        exp_dist = np.exp( -0.5* np.sum( np.square(  np.sin( base ) / self.lengthscales ), axis = -1 ) ) 
-            
+
+        base = np.pi * (X[:, None, :] - X2[None, :, :]) / self.period
+        exp_dist = np.exp( -0.5* np.sum( np.square(  np.sin( base ) / self.lengthscale ), axis = -1 ) )
+
         return self.variance * exp_dist
 
 
@@ -125,42 +124,47 @@ class StdPeriodic(Kern):
 
     def update_gradients_full(self, dL_dK, X, X2=None):
         """derivative of the covariance matrix with respect to the parameters."""
-        if X2 is None: 
+        if X2 is None:
             X2 = X
-        
-        base = np.pi * (X[:, None, :] - X2[None, :, :]) / self.wavelengths
-        
-        sin_base = np.sin( base )         
-        exp_dist = np.exp( -0.5* np.sum( np.square(  sin_base / self.lengthscales ), axis = -1 ) ) 
-        
-        dwl = self.variance * (1.0/np.square(self.lengthscales)) * sin_base*np.cos(base) * (base / self.wavelengths)
-        
-        dl = self.variance * np.square( sin_base) / np.power( self.lengthscales, 3) 
-        
-        self.variance.gradient = np.sum(exp_dist * dL_dK)    
-        #target[0] += np.sum( exp_dist * dL_dK)        
-        
-        if self.ARD1: # different wavelengths
-            self.wavelengths.gradient = (dwl * exp_dist[:,:,None] * dL_dK[:, :, None]).sum(0).sum(0)
-        else:  # same wavelengths
-            self.wavelengths.gradient = np.sum(dwl.sum(-1) * exp_dist * dL_dK)
-            
+
+        base = np.pi * (X[:, None, :] - X2[None, :, :]) / self.period
+
+        sin_base = np.sin( base )
+        exp_dist = np.exp( -0.5* np.sum( np.square(  sin_base / self.lengthscale ), axis = -1 ) )
+
+        dwl = self.variance * (1.0/np.square(self.lengthscale)) * sin_base*np.cos(base) * (base / self.period)
+
+        dl = self.variance * np.square( sin_base) / np.power( self.lengthscale, 3)
+
+        self.variance.gradient = np.sum(exp_dist * dL_dK)
+        #target[0] += np.sum( exp_dist * dL_dK)
+
+        if self.ARD1: # different periods
+            self.period.gradient = (dwl * exp_dist[:,:,None] * dL_dK[:, :, None]).sum(0).sum(0)
+        else:  # same period
+            self.period.gradient = np.sum(dwl.sum(-1) * exp_dist * dL_dK)
+
         if self.ARD2: # different lengthscales
-            self.lengthscales.gradient = (dl * exp_dist[:,:,None] * dL_dK[:, :, None]).sum(0).sum(0)
+            self.lengthscale.gradient = (dl * exp_dist[:,:,None] * dL_dK[:, :, None]).sum(0).sum(0)
         else: # same lengthscales
-            self.lengthscales.gradient = np.sum(dl.sum(-1) * exp_dist * dL_dK)
-        
+            self.lengthscale.gradient = np.sum(dl.sum(-1) * exp_dist * dL_dK)
+
     def update_gradients_diag(self, dL_dKdiag, X):
         """derivative of the diagonal of the covariance matrix with respect to the parameters."""
         self.variance.gradient = np.sum(dL_dKdiag)
-        self.wavelengths.gradient = 0
-        self.lengthscales.gradient = 0
+        self.period.gradient = 0
+        self.lengthscale.gradient = 0
 
-#    def gradients_X(self, dL_dK, X, X2=None):
-#        """derivative of the covariance matrix with respect to X."""
-#    
-#        raise NotImplemented("Periodic kernel: dK_dX not implemented")
-#
-#    def gradients_X_diag(self, dL_dKdiag, X):
-#        
-#        raise NotImplemented("Periodic kernel: dKdiag_dX not implemented")
+    def gradients_X(self, dL_dK, X, X2=None):
+        K = self.K(X, X2)
+        if X2 is None:
+            dL_dK = dL_dK+dL_dK.T
+            X2 = X
+        dX = -np.pi*((dL_dK*K)[:,:,None]*np.sin(2*np.pi/self.period*(X[:,None,:] - X2[None,:,:]))/(2.*np.square(self.lengthscale)*self.period)).sum(1)
+        return dX
+    
+    def gradients_X_diag(self, dL_dKdiag, X):
+        return np.zeros(X.shape)
+    
+    def input_sensitivity(self, summarize=True):
+        return self.variance*np.ones(self.input_dim)/self.lengthscale**2

GPy/kern/src/static.py

@@ -81,6 +81,52 @@ class White(Static):
     def update_gradients_expectations(self, dL_dpsi0, dL_dpsi1, dL_dpsi2, Z, variational_posterior):
         self.variance.gradient = dL_dpsi0.sum()
 
+class WhiteHeteroscedastic(Static):
+    def __init__(self, input_dim, num_data, variance=1., active_dims=None, name='white_hetero'):
+        """
+        A heteroscedastic White kernel (nugget/noise).
+        It defines one variance (nugget) per input sample. 
+        
+        Prediction excludes any noise learnt by this Kernel, so be careful using this kernel.
+        
+        You can plot the errors learnt by this kernel by something similar as:
+        plt.errorbar(m.X, m.Y, yerr=2*np.sqrt(m.kern.white.variance))
+        """
+        super(Static, self).__init__(input_dim, active_dims, name)
+        self.variance = Param('variance', np.ones(num_data) * variance, Logexp())
+        self.link_parameters(self.variance)
+
+    def Kdiag(self, X):
+        if X.shape[0] == self.variance.shape[0]:
+            # If the input has the same number of samples as 
+            # the number of variances, we return the variances
+            return self.variance
+        return 0.
+
+    def K(self, X, X2=None):
+        if X2 is None and X.shape[0] == self.variance.shape[0]:
+            return np.eye(X.shape[0]) * self.variance
+        else:
+            return 0.
+
+    def psi2(self, Z, variational_posterior):
+        return np.zeros((Z.shape[0], Z.shape[0]), dtype=np.float64)
+
+    def psi2n(self, Z, variational_posterior):
+        return np.zeros((1, Z.shape[0], Z.shape[0]), dtype=np.float64)
+
+    def update_gradients_full(self, dL_dK, X, X2=None):
+        if X2 is None:
+            self.variance.gradient = np.diagonal(dL_dK)
+        else:
+            self.variance.gradient = 0.
+
+    def update_gradients_diag(self, dL_dKdiag, X):
+        self.variance.gradient = dL_dKdiag
+
+    def update_gradients_expectations(self, dL_dpsi0, dL_dpsi1, dL_dpsi2, Z, variational_posterior):
+        self.variance.gradient = dL_dpsi0
+
 class Bias(Static):
     def __init__(self, input_dim, variance=1., active_dims=None, name='bias'):
         super(Bias, self).__init__(input_dim, variance, active_dims, name)

GPy/kern/src/stationary.py

@@ -81,11 +81,11 @@ class Stationary(Kern):
     def dK_dr(self, r):
         raise NotImplementedError("implement derivative of the covariance function wrt r to use this class")
 
-    @Cache_this(limit=20, ignore_args=())
+    @Cache_this(limit=3, ignore_args=())
     def dK2_drdr(self, r):
         raise NotImplementedError("implement second derivative of covariance wrt r to use this method")
 
-    @Cache_this(limit=5, ignore_args=())
+    @Cache_this(limit=3, ignore_args=())
     def K(self, X, X2=None):
         """
         Kernel function applied on inputs X and X2.
@@ -99,6 +99,9 @@ class Stationary(Kern):
 
     @Cache_this(limit=3, ignore_args=())
     def dK_dr_via_X(self, X, X2):
+        """
+        compute the derivative of K wrt X going through X
+        """
         #a convenience function, so we can cache dK_dr
         return self.dK_dr(self._scaled_dist(X, X2))
 
@@ -486,18 +489,21 @@ class RatQuad(Stationary):
         self.link_parameters(self.power)
 
     def K_of_r(self, r):
-        r2 = np.power(r, 2.)
-        return self.variance*np.power(1. + r2/2., -self.power)
+        r2 = np.square(r)
+#         return self.variance*np.power(1. + r2/2., -self.power)
+        return self.variance*np.exp(-self.power*np.log1p(r2/2.))
 
     def dK_dr(self, r):
-        r2 = np.power(r, 2.)
-        return -self.variance*self.power*r*np.power(1. + r2/2., - self.power - 1.)
+        r2 = np.square(r)
+#         return -self.variance*self.power*r*np.power(1. + r2/2., - self.power - 1.)
+        return-self.variance*self.power*r*np.exp(-(self.power+1)*np.log1p(r2/2.))
 
     def update_gradients_full(self, dL_dK, X, X2=None):
         super(RatQuad, self).update_gradients_full(dL_dK, X, X2)
         r = self._scaled_dist(X, X2)
-        r2 = np.power(r, 2.)
-        dK_dpow = -self.variance * np.power(2., self.power) * np.power(r2 + 2., -self.power) * np.log(0.5*(r2+2.))
+        r2 = np.square(r)
+#        dK_dpow = -self.variance * np.power(2., self.power) * np.power(r2 + 2., -self.power) * np.log(0.5*(r2+2.))
+        dK_dpow = -self.variance * np.exp(self.power*(np.log(2.)-np.log1p(r2+1)))*np.log1p(r2/2.)
         grad = np.sum(dL_dK*dK_dpow)
         self.power.gradient = grad
 

GPy/kern/src/trunclinear.py

@@ -54,12 +54,12 @@ class TruncLinear(Kern):
         self.add_parameter(self.variances)
         self.add_parameter(self.delta)
 
-    @Cache_this(limit=2)
+    @Cache_this(limit=3)
     def K(self, X, X2=None):
         XX = self.variances*self._product(X, X2)
         return XX.sum(axis=-1)
 
-    @Cache_this(limit=2)
+    @Cache_this(limit=3)
     def _product(self, X, X2=None):
         if X2 is None:
             X2 = X
@@ -149,12 +149,12 @@ class TruncLinear_inf(Kern):
         self.add_parameter(self.variances)
 
 
-#     @Cache_this(limit=2)
+#     @Cache_this(limit=3)
     def K(self, X, X2=None):
         tmp = self._product(X, X2)
         return (self.variances*tmp).sum(axis=-1)
 
-#     @Cache_this(limit=2)
+#     @Cache_this(limit=3)
     def _product(self, X, X2=None):
         if X2 is None:
             X2 = X

GPy/models/bayesian_gplvm.py

@@ -61,7 +61,7 @@ class BayesianGPLVM(SparseGP_MPI):
             else:
                 from ..inference.latent_function_inference.var_dtc import VarDTC
                 self.logger.debug("creating inference_method var_dtc")
-                inference_method = VarDTC(limit=1 if not missing_data else Y.shape[1])
+                inference_method = VarDTC(limit=3 if not missing_data else Y.shape[1])
         if isinstance(inference_method,VarDTC_minibatch):
             inference_method.mpi_comm = mpi_comm
 

GPy/models/bayesian_gplvm_minibatch.py

@@ -40,12 +40,13 @@ class BayesianGPLVMMiniBatch(SparseGPMiniBatch):
             Z = np.random.permutation(X.copy())[:num_inducing]
         assert Z.shape[1] == X.shape[1]
 
-        if X_variance == False:
+        if X_variance is False:
             self.logger.info('no variance on X, activating sparse GPLVM')
             X = Param("latent space", X)
-        elif X_variance is None:
-            self.logger.info("initializing latent space variance ~ uniform(0,.1)")
-            X_variance = np.random.uniform(0,.1,X.shape)
+        else:
+            if X_variance is None:
+                self.logger.info("initializing latent space variance ~ uniform(0,.1)")
+                X_variance = np.random.uniform(0,.1,X.shape)
             self.variational_prior = NormalPrior()
             X = NormalPosterior(X, X_variance)
 
@@ -61,7 +62,7 @@ class BayesianGPLVMMiniBatch(SparseGPMiniBatch):
         if inference_method is None:
             from ..inference.latent_function_inference.var_dtc import VarDTC
             self.logger.debug("creating inference_method var_dtc")
-            inference_method = VarDTC(limit=1 if not missing_data else Y.shape[1])
+            inference_method = VarDTC(limit=3 if not missing_data else Y.shape[1])
 
         super(BayesianGPLVMMiniBatch,self).__init__(X, Y, Z, kernel, likelihood=likelihood,
                                            name=name, inference_method=inference_method,
@@ -71,13 +72,13 @@ class BayesianGPLVMMiniBatch(SparseGPMiniBatch):
         self.X = X
         self.link_parameter(self.X, 0)
 
-    def set_X_gradients(self, X, X_grad):
-        """Set the gradients of the posterior distribution of X in its specific form."""
-        X.mean.gradient, X.variance.gradient = X_grad
+    #def set_X_gradients(self, X, X_grad):
+    #    """Set the gradients of the posterior distribution of X in its specific form."""
+    #    X.mean.gradient, X.variance.gradient = X_grad
 
-    def get_X_gradients(self, X):
-        """Get the gradients of the posterior distribution of X in its specific form."""
-        return X.mean.gradient, X.variance.gradient
+    #def get_X_gradients(self, X):
+    #    """Get the gradients of the posterior distribution of X in its specific form."""
+    #    return X.mean.gradient, X.variance.gradient
 
     def _outer_values_update(self, full_values):
         """
@@ -106,7 +107,7 @@ class BayesianGPLVMMiniBatch(SparseGPMiniBatch):
         super(BayesianGPLVMMiniBatch,self).parameters_changed()
 
         kl_fctr = self.kl_factr
-        if kl_fctr > 0:
+        if kl_fctr > 0 and self.has_uncertain_inputs():
             Xgrad = self.X.gradient.copy()
             self.X.gradient[:] = 0
             self.variational_prior.update_gradients_KL(self.X)
@@ -122,8 +123,8 @@ class BayesianGPLVMMiniBatch(SparseGPMiniBatch):
 
             if self.missing_data or not self.stochastics:
                 self._log_marginal_likelihood -= kl_fctr*self.variational_prior.KL_divergence(self.X)
-            elif self.stochastics:
+            else: #self.stochastics is given:
                 d = self.output_dim
                 self._log_marginal_likelihood -= kl_fctr*self.variational_prior.KL_divergence(self.X)*self.stochastics.batchsize/d
 
-        self._Xgrad = self.X.gradient.copy()
+        self._Xgrad = self.X.gradient.copy()

GPy/models/gplvm.py

@@ -41,4 +41,4 @@ class GPLVM(GP):
 
     def parameters_changed(self):
         super(GPLVM, self).parameters_changed()
-        self.X.gradient = self.kern.gradients_X(self.grad_dict['dL_dK'], self.X, None)
+        self.X.gradient = self.kern.gradients_X(self.grad_dict['dL_dK'], self.X, None)

GPy/models/mrd.py

@@ -5,14 +5,14 @@ import numpy as np
 import itertools, logging
 
 from ..kern import Kern
-from GPy.core.parameterization.variational import NormalPrior
+from ..core.parameterization.variational import NormalPrior
 from ..core.parameterization import Param
 from paramz import ObsAr
 from ..inference.latent_function_inference.var_dtc import VarDTC
 from ..inference.latent_function_inference import InferenceMethodList
 from ..likelihoods import Gaussian
 from ..util.initialization import initialize_latent
-from GPy.models.bayesian_gplvm_minibatch import BayesianGPLVMMiniBatch
+from ..models.bayesian_gplvm_minibatch import BayesianGPLVMMiniBatch
 
 class MRD(BayesianGPLVMMiniBatch):
     """
@@ -215,40 +215,6 @@ class MRD(BayesianGPLVMMiniBatch):
             Z = np.random.randn(self.num_inducing, self.input_dim) * X.var()
         return Z
 
-    def _handle_plotting(self, fignum, axes, plotf, sharex=False, sharey=False):
-        import matplotlib.pyplot as plt
-        if axes is None:
-            fig = plt.figure(num=fignum)
-        sharex_ax = None
-        sharey_ax = None
-        plots = []
-        for i, g in enumerate(self.bgplvms):
-            try:
-                if sharex:
-                    sharex_ax = ax # @UndefinedVariable
-                    sharex = False # dont set twice
-                if sharey:
-                    sharey_ax = ax # @UndefinedVariable
-                    sharey = False # dont set twice
-            except:
-                pass
-            if axes is None:
-                ax = fig.add_subplot(1, len(self.bgplvms), i + 1, sharex=sharex_ax, sharey=sharey_ax)
-            elif isinstance(axes, (tuple, list, np.ndarray)):
-                ax = axes[i]
-            else:
-                raise ValueError("Need one axes per latent dimension input_dim")
-            plots.append(plotf(i, g, ax))
-            if sharey_ax is not None:
-                plt.setp(ax.get_yticklabels(), visible=False)
-        plt.draw()
-        if axes is None:
-            try:
-                fig.tight_layout()
-            except:
-                pass
-        return plots
-
     def predict(self, Xnew, full_cov=False, Y_metadata=None, kern=None, Yindex=0):
         """
         Prediction for data set Yindex[default=0].
@@ -270,59 +236,50 @@ class MRD(BayesianGPLVMMiniBatch):
     #                                     sharex=sharex, sharey=sharey)
     #         return fig
 
-    def plot_scales(self, fignum=None, ax=None, titles=None, sharex=False, sharey=True, *args, **kwargs):
+    def plot_scales(self, titles=None, fig_kwargs={}, **kwargs):
         """
-
-        TODO: Explain other parameters
+        Plot input sensitivity for all datasets, to see which input dimensions are
+        significant for which dataset.
 
         :param titles: titles for axes of datasets
 
+        kwargs go into plot_ARD for each kernel.
         """
+        from ..plotting import plotting_library as pl
+
         if titles is None:
             titles = [r'${}$'.format(name) for name in self.names]
-        ymax = reduce(max, [np.ceil(max(g.kern.input_sensitivity())) for g in self.bgplvms])
-        def plotf(i, g, ax):
-            #ax.set_ylim([0,ymax])
-            return g.kern.plot_ARD(ax=ax, title=titles[i], *args, **kwargs)
-        fig = self._handle_plotting(fignum, ax, plotf, sharex=sharex, sharey=sharey)
-        return fig
+
+        M = len(self.bgplvms)
+        fig = pl().figure(rows=1, cols=M, **fig_kwargs)
+        for c in range(M):
+            canvas = self.bgplvms[c].kern.plot_ARD(title=titles[c], figure=fig, col=c+1, **kwargs)
+        return canvas
 
     def plot_latent(self, labels=None, which_indices=None,
-                resolution=50, ax=None, marker='o', s=40,
-                fignum=None, plot_inducing=True, legend=True,
+                resolution=60, legend=True,
                 plot_limits=None,
-                aspect='auto', updates=False, predict_kwargs={}, imshow_kwargs={}):
+                updates=False,
+                kern=None, marker='<>^vsd',
+                num_samples=1000, projection='2d',
+                predict_kwargs={},
+                scatter_kwargs=None, **imshow_kwargs):
         """
         see plotting.matplot_dep.dim_reduction_plots.plot_latent
         if predict_kwargs is None, will plot latent spaces for 0th dataset (and kernel), otherwise give
         predict_kwargs=dict(Yindex='index') for plotting only the latent space of dataset with 'index'.
         """
-        import sys
-        assert "matplotlib" in sys.modules, "matplotlib package has not been imported."
-        from matplotlib import pyplot as plt
-        from ..plotting.matplot_dep import dim_reduction_plots
+        from ..plotting.gpy_plot.latent_plots import plot_latent
+
         if "Yindex" not in predict_kwargs:
             predict_kwargs['Yindex'] = 0
 
         Yindex = predict_kwargs['Yindex']
-        if ax is None:
-            fig = plt.figure(num=fignum)
-            ax = fig.add_subplot(111)
-        else:
-            fig = ax.figure
+
         self.kern = self.bgplvms[Yindex].kern
         self.likelihood = self.bgplvms[Yindex].likelihood
-        plot = dim_reduction_plots.plot_latent(self, labels, which_indices,
-                                        resolution, ax, marker, s,
-                                        fignum, plot_inducing, legend,
-                                        plot_limits, aspect, updates, predict_kwargs, imshow_kwargs)
-        ax.set_title(self.bgplvms[Yindex].name)
-        try:
-            fig.tight_layout()
-        except:
-            pass
 
-        return plot
+        return plot_latent(self, labels, which_indices, resolution, legend, plot_limits, updates, kern, marker, num_samples, projection, scatter_kwargs)
 
     def __getstate__(self):
         state = super(MRD, self).__getstate__()

GPy/models/sparse_gp_minibatch.py

@@ -41,11 +41,12 @@ class SparseGPMiniBatch(SparseGP):
     def __init__(self, X, Y, Z, kernel, likelihood, inference_method=None,
                  name='sparse gp', Y_metadata=None, normalizer=False,
                  missing_data=False, stochastic=False, batchsize=1):
+        self._update_stochastics = False
 
         # pick a sensible inference method
         if inference_method is None:
             if isinstance(likelihood, likelihoods.Gaussian):
-                inference_method = var_dtc.VarDTC(limit=1 if not missing_data else Y.shape[1])
+                inference_method = var_dtc.VarDTC(limit=3 if not missing_data else Y.shape[1])
             else:
                 #inference_method = ??
                 raise NotImplementedError("what to do what to do?")
@@ -73,7 +74,14 @@ class SparseGPMiniBatch(SparseGP):
         logger.info("Adding Z as parameter")
         self.link_parameter(self.Z, index=0)
         self.posterior = None
-
+        
+    def optimize(self, optimizer=None, start=None, **kwargs):
+        try:
+            self._update_stochastics = True
+            SparseGP.optimize(self, optimizer=optimizer, start=start, **kwargs)
+        finally:
+            self._update_stochastics = False
+            
     def has_uncertain_inputs(self):
         return isinstance(self.X, VariationalPosterior)
 
@@ -226,16 +234,16 @@ class SparseGPMiniBatch(SparseGP):
             woodbury_inv = self.posterior._woodbury_inv
             woodbury_vector = self.posterior._woodbury_vector
 
-        if not self.stochastics:
-            m_f = lambda i: "Inference with missing_data: {: >7.2%}".format(float(i+1)/self.output_dim)
-            message = m_f(-1)
-            print(message, end=' ')
+        #if not self.stochastics:
+        #    m_f = lambda i: "Inference with missing_data: {: >7.2%}".format(float(i+1)/self.output_dim)
+        #    message = m_f(-1)
+        #    print(message, end=' ')
 
         for d, ninan in self.stochastics.d:
-            if not self.stochastics:
-                print(' '*(len(message)) + '\r', end=' ')
-                message = m_f(d)
-                print(message, end=' ')
+            #if not self.stochastics:
+            #    print(' '*(len(message)) + '\r', end=' ')
+            #    message = m_f(d)
+            #    print(message, end=' ')
 
             psi0ni = self.psi0[ninan]
             psi1ni = self.psi1[ninan]
@@ -262,8 +270,8 @@ class SparseGPMiniBatch(SparseGP):
             woodbury_vector[:, d] = posterior.woodbury_vector
             self._log_marginal_likelihood += log_marginal_likelihood
 
-        if not self.stochastics:
-            print('')
+        #if not self.stochastics:
+        #    print('')
 
         if self.posterior is None:
             self.posterior = Posterior(woodbury_inv=woodbury_inv, woodbury_vector=woodbury_vector,
@@ -314,6 +322,8 @@ class SparseGPMiniBatch(SparseGP):
         if self.missing_data:
             self._outer_loop_for_missing_data()
         elif self.stochastics:
+            if self._update_stochastics:
+                self.stochastics.do_stochastics()
             self._outer_loop_without_missing_data()
         else:
             self.posterior, self._log_marginal_likelihood, self.grad_dict = self._inner_parameters_changed(self.kern, self.X, self.Z, self.likelihood, self.Y_normalized, self.Y_metadata)

GPy/models/sparse_gp_regression.py

@@ -62,4 +62,4 @@ class SparseGPRegression(SparseGP_MPI):
         if isinstance(self.inference_method,VarDTC_minibatch):
             update_gradients_sparsegp(self, mpi_comm=self.mpi_comm)
         else:
-            super(SparseGPRegression, self).parameters_changed()
+            super(SparseGPRegression, self).parameters_changed()

GPy/models/sparse_gplvm.py

@@ -4,6 +4,7 @@
 
 import sys
 from .sparse_gp_regression import SparseGPRegression
+from ..core import Param
 
 class SparseGPLVM(SparseGPRegression):
     """
@@ -21,7 +22,9 @@ class SparseGPLVM(SparseGPRegression):
         if X is None:
             from ..util.initialization import initialize_latent
             X, fracs = initialize_latent(init, input_dim, Y)
+        X = Param('latent space', X)
         SparseGPRegression.__init__(self, X, Y, kernel=kernel, num_inducing=num_inducing)
+        self.link_parameter(self.X, 0)
 
     def parameters_changed(self):
         super(SparseGPLVM, self).parameters_changed()

GPy/models/ss_gplvm.py

@@ -94,6 +94,86 @@ class IBPPrior(VariationalPrior):
         variational_posterior.tau.gradient[:,0] = -((tau[:,0]-gamma-ad)*polygamma(1,tau[:,0])+common)
         variational_posterior.tau.gradient[:,1] = -((tau[:,1]+gamma-2)*polygamma(1,tau[:,1])+common)
 
+class SLVMPosterior(SpikeAndSlabPosterior):
+    '''
+    The SpikeAndSlab distribution for variational approximations.
+    '''
+    def __init__(self, means, variances, binary_prob, tau=None,  name='latent space'):
+        """
+        binary_prob : the probability of the distribution on the slab part.
+        """
+        from paramz.transformations import Logexp
+        super(SLVMPosterior, self).__init__(means, variances, binary_prob, group_spike=False, name=name)
+        self.tau = Param("tau_", np.ones((self.gamma.shape[1],2)), Logexp())
+        self.link_parameter(self.tau)
+
+    def set_gradients(self, grad):
+        self.mean.gradient, self.variance.gradient, self.gamma.gradient, self.tau.gradient = grad
+
+    def __getitem__(self, s):
+        if isinstance(s, (int, slice, tuple, list, np.ndarray)):
+            import copy
+            n = self.__new__(self.__class__, self.name)
+            dc = self.__dict__.copy()
+            dc['mean'] = self.mean[s]
+            dc['variance'] = self.variance[s]
+            dc['binary_prob'] = self.binary_prob[s]
+            dc['tau'] = self.tau
+            dc['parameters'] = copy.copy(self.parameters)
+            n.__dict__.update(dc)
+            n.parameters[dc['mean']._parent_index_] = dc['mean']
+            n.parameters[dc['variance']._parent_index_] = dc['variance']
+            n.parameters[dc['binary_prob']._parent_index_] = dc['binary_prob']
+            n.parameters[dc['tau']._parent_index_] = dc['tau']
+            n._gradient_array_ = None
+            oversize = self.size - self.mean.size - self.variance.size - self.gamma.size - self.tau.size
+            n.size = n.mean.size + n.variance.size + n.gamma.size+ n.tau.size + oversize
+            n.ndim = n.mean.ndim
+            n.shape = n.mean.shape
+            n.num_data = n.mean.shape[0]
+            n.input_dim = n.mean.shape[1] if n.ndim != 1 else 1
+            return n
+        else:
+            return super(IBPPosterior, self).__getitem__(s)
+
+class SLVMPrior(VariationalPrior):
+    def __init__(self, input_dim, alpha =1., beta=1., Z=None, name='SLVMPrior', **kw):
+        super(SLVMPrior, self).__init__(name=name, **kw)
+        self.input_dim = input_dim
+        self.variance = 1.
+        self.alpha = alpha
+        self.beta = beta
+        self.Z = Z
+        if Z is not None:
+            assert np.all(np.unique(Z)==np.array([0,1]))
+
+    def KL_divergence(self, variational_posterior):
+        mu, S, gamma, tau = variational_posterior.mean.values, variational_posterior.variance.values, variational_posterior.gamma.values, variational_posterior.tau.values
+        
+        var_mean = np.square(mu)/self.variance
+        var_S = (S/self.variance - np.log(S))
+        part1 = (gamma* (np.log(self.variance)-1. +var_mean + var_S)).sum()/2.
+        
+        from scipy.special import betaln,digamma            
+        part2 = (gamma*np.log(gamma)).sum() + ((1.-gamma)*np.log(1.-gamma)).sum() + betaln(self.alpha,self.beta)*self.input_dim \
+                -betaln(tau[:,0], tau[:,1]).sum() + ((tau[:,0]-(gamma*self.Z).sum(0)-self.alpha)*digamma(tau[:,0])).sum() + \
+                ((tau[:,1]-((1-gamma)*self.Z).sum(0)-self.beta)*digamma(tau[:,1])).sum() + ((self.Z.sum(0)+self.alpha+self.beta-tau[:,0]-tau[:,1])*digamma(tau.sum(axis=1))).sum()        
+
+        return part1+part2
+
+    def update_gradients_KL(self, variational_posterior):
+        mu, S, gamma, tau = variational_posterior.mean.values, variational_posterior.variance.values, variational_posterior.gamma.values, variational_posterior.tau.values
+
+        variational_posterior.mean.gradient -= gamma*mu/self.variance
+        variational_posterior.variance.gradient -= (1./self.variance - 1./S) * gamma /2.
+        from scipy.special import digamma,polygamma
+        dgamma = np.log(gamma/(1.-gamma))+ (digamma(tau[:,1])-digamma(tau[:,0]))*self.Z
+        variational_posterior.binary_prob.gradient -= dgamma+((np.square(mu)+S)/self.variance-np.log(S)+np.log(self.variance)-1.)/2.
+        common = (self.Z.sum(0)+self.alpha+self.beta-tau[:,0]-tau[:,1])*polygamma(1,tau.sum(axis=1))
+        variational_posterior.tau.gradient[:,0] = -((tau[:,0]-(gamma*self.Z).sum(0)-self.alpha)*polygamma(1,tau[:,0])+common)
+        variational_posterior.tau.gradient[:,1] = -((tau[:,1]-((1-gamma)*self.Z).sum(0)-self.beta)*polygamma(1,tau[:,1])+common)
+
+
 class SSGPLVM(SparseGP_MPI):
     """
     Spike-and-Slab Gaussian Process Latent Variable Model
@@ -107,7 +187,7 @@ class SSGPLVM(SparseGP_MPI):
 
     """
     def __init__(self, Y, input_dim, X=None, X_variance=None, Gamma=None, init='PCA', num_inducing=10,
-                 Z=None, kernel=None, inference_method=None, likelihood=None, name='Spike_and_Slab GPLVM', group_spike=False, IBP=False, alpha=2., tau=None, mpi_comm=None, pi=None, learnPi=False,normalizer=False, sharedX=False, variational_prior=None,**kwargs):
+                 Z=None, kernel=None, inference_method=None, likelihood=None, name='Spike_and_Slab GPLVM', group_spike=False, IBP=False,SLVM=False, alpha=2., beta=2., connM=None, tau=None, mpi_comm=None, pi=None, learnPi=False,normalizer=False, sharedX=False, variational_prior=None,**kwargs):
 
         self.group_spike = group_spike
         self.init = init
@@ -152,6 +232,9 @@ class SSGPLVM(SparseGP_MPI):
         if IBP:
             self.variational_prior = IBPPrior(input_dim=input_dim, alpha=alpha) if variational_prior is None else variational_prior
             X = IBPPosterior(X, X_variance, gamma, tau=tau,sharedX=sharedX)
+        elif SLVM:
+            self.variational_prior = SLVMPrior(input_dim=input_dim, alpha=alpha, beta=beta, Z=connM) if variational_prior is None else variational_prior
+            X = SLVMPosterior(X, X_variance, gamma, tau=tau)
         else:
             self.variational_prior = SpikeAndSlabPrior(pi=pi,learnPi=learnPi, group_spike=group_spike)  if variational_prior is None else variational_prior
             X = SpikeAndSlabPosterior(X, X_variance, gamma, group_spike=group_spike,sharedX=sharedX)

GPy/plotting/Tango.py

@@ -104,4 +104,4 @@ cdict_Alu = {'red' :((0./5,colorsRGB['Aluminium1'][0]/256.,colorsRGB['Aluminium1
                      (2./5,colorsRGB['Aluminium3'][2]/256.,colorsRGB['Aluminium3'][2]/256.),
                      (3./5,colorsRGB['Aluminium4'][2]/256.,colorsRGB['Aluminium4'][2]/256.),
                      (4./5,colorsRGB['Aluminium5'][2]/256.,colorsRGB['Aluminium5'][2]/256.),
-                     (5./5,colorsRGB['Aluminium6'][2]/256.,colorsRGB['Aluminium6'][2]/256.))}
+                     (5./5,colorsRGB['Aluminium6'][2]/256.,colorsRGB['Aluminium6'][2]/256.))}

GPy/plotting/__init__.py

@@ -25,18 +25,66 @@ def change_plotting_library(lib):
             current_lib[0] = PlotlyPlots()
         if lib == 'none':
             current_lib[0] = None
+        inject_plotting()
         #===========================================================================
     except (ImportError, NameError):
         config.set('plotting', 'library', 'none')
+        raise
         import warnings
         warnings.warn(ImportWarning("You spevified {} in your configuration, but is not available. Install newest version of {} for plotting".format(lib, lib)))
 
-from ..util.config import config, NoOptionError
-try:
-    lib = config.get('plotting', 'library')
-    change_plotting_library(lib)
-except NoOptionError:
-    print("No plotting library was specified in config file. \n{}".format(error_suggestion))
+def inject_plotting():
+    if current_lib[0] is not None:
+        # Inject the plots into classes here:
+
+        # Already converted to new style:
+        from . import gpy_plot
+
+        from ..core import GP
+        GP.plot_data = gpy_plot.data_plots.plot_data
+        GP.plot_data_error = gpy_plot.data_plots.plot_data_error
+        GP.plot_errorbars_trainset = gpy_plot.data_plots.plot_errorbars_trainset
+        GP.plot_mean = gpy_plot.gp_plots.plot_mean
+        GP.plot_confidence = gpy_plot.gp_plots.plot_confidence
+        GP.plot_density = gpy_plot.gp_plots.plot_density
+        GP.plot_samples = gpy_plot.gp_plots.plot_samples
+        GP.plot = gpy_plot.gp_plots.plot
+        GP.plot_f = gpy_plot.gp_plots.plot_f
+        GP.plot_magnification = gpy_plot.latent_plots.plot_magnification
+
+        from ..core import SparseGP
+        SparseGP.plot_inducing = gpy_plot.data_plots.plot_inducing
+
+        from ..models import GPLVM, BayesianGPLVM, bayesian_gplvm_minibatch, SSGPLVM, SSMRD
+        GPLVM.plot_latent = gpy_plot.latent_plots.plot_latent
+        GPLVM.plot_scatter = gpy_plot.latent_plots.plot_latent_scatter
+        GPLVM.plot_inducing = gpy_plot.latent_plots.plot_latent_inducing
+        GPLVM.plot_steepest_gradient_map = gpy_plot.latent_plots.plot_steepest_gradient_map
+        BayesianGPLVM.plot_latent = gpy_plot.latent_plots.plot_latent
+        BayesianGPLVM.plot_scatter = gpy_plot.latent_plots.plot_latent_scatter
+        BayesianGPLVM.plot_inducing = gpy_plot.latent_plots.plot_latent_inducing
+        BayesianGPLVM.plot_steepest_gradient_map = gpy_plot.latent_plots.plot_steepest_gradient_map
+        bayesian_gplvm_minibatch.BayesianGPLVMMiniBatch.plot_latent = gpy_plot.latent_plots.plot_latent
+        bayesian_gplvm_minibatch.BayesianGPLVMMiniBatch.plot_scatter = gpy_plot.latent_plots.plot_latent_scatter
+        bayesian_gplvm_minibatch.BayesianGPLVMMiniBatch.plot_inducing = gpy_plot.latent_plots.plot_latent_inducing
+        bayesian_gplvm_minibatch.BayesianGPLVMMiniBatch.plot_steepest_gradient_map = gpy_plot.latent_plots.plot_steepest_gradient_map
+        SSGPLVM.plot_latent = gpy_plot.latent_plots.plot_latent
+        SSGPLVM.plot_scatter = gpy_plot.latent_plots.plot_latent_scatter
+        SSGPLVM.plot_inducing = gpy_plot.latent_plots.plot_latent_inducing
+        SSGPLVM.plot_steepest_gradient_map = gpy_plot.latent_plots.plot_steepest_gradient_map
+
+        from ..kern import Kern
+        Kern.plot_covariance = gpy_plot.kernel_plots.plot_covariance
+        def deprecate_plot(self, *args, **kwargs):
+            import warnings
+            warnings.warn(DeprecationWarning('Kern.plot is being deprecated and will not be available in the 1.0 release. Use Kern.plot_covariance instead'))
+            return self.plot_covariance(*args, **kwargs)
+        Kern.plot = deprecate_plot
+        Kern.plot_ARD = gpy_plot.kernel_plots.plot_ARD
+
+        from ..inference.optimization import Optimizer
+        Optimizer.plot = gpy_plot.inference_plots.plot_optimizer
+        # Variational plot!
 
 def plotting_library():
     if current_lib[0] is None:
@@ -53,54 +101,10 @@ def show(figure, **kwargs):
     """
     return plotting_library().show_canvas(figure, **kwargs)
 
-if config.get('plotting', 'library') is not 'none':
-    # Inject the plots into classes here:
 
-    # Already converted to new style:
-    from . import gpy_plot
-
-    from ..core import GP
-    GP.plot_data = gpy_plot.data_plots.plot_data
-    GP.plot_data_error = gpy_plot.data_plots.plot_data_error
-    GP.plot_errorbars_trainset = gpy_plot.data_plots.plot_errorbars_trainset
-    GP.plot_mean = gpy_plot.gp_plots.plot_mean
-    GP.plot_confidence = gpy_plot.gp_plots.plot_confidence
-    GP.plot_density = gpy_plot.gp_plots.plot_density
-    GP.plot_samples = gpy_plot.gp_plots.plot_samples
-    GP.plot = gpy_plot.gp_plots.plot
-    GP.plot_f = gpy_plot.gp_plots.plot_f
-    GP.plot_magnification = gpy_plot.latent_plots.plot_magnification
-
-    from ..core import SparseGP
-    SparseGP.plot_inducing = gpy_plot.data_plots.plot_inducing
-
-    from ..models import GPLVM, BayesianGPLVM, bayesian_gplvm_minibatch, SSGPLVM, SSMRD
-    GPLVM.plot_latent = gpy_plot.latent_plots.plot_latent
-    GPLVM.plot_scatter = gpy_plot.latent_plots.plot_latent_scatter
-    GPLVM.plot_inducing = gpy_plot.latent_plots.plot_latent_inducing
-    GPLVM.plot_steepest_gradient_map = gpy_plot.latent_plots.plot_steepest_gradient_map
-    BayesianGPLVM.plot_latent = gpy_plot.latent_plots.plot_latent
-    BayesianGPLVM.plot_scatter = gpy_plot.latent_plots.plot_latent_scatter
-    BayesianGPLVM.plot_inducing = gpy_plot.latent_plots.plot_latent_inducing
-    BayesianGPLVM.plot_steepest_gradient_map = gpy_plot.latent_plots.plot_steepest_gradient_map
-    bayesian_gplvm_minibatch.BayesianGPLVMMiniBatch.plot_latent = gpy_plot.latent_plots.plot_latent
-    bayesian_gplvm_minibatch.BayesianGPLVMMiniBatch.plot_scatter = gpy_plot.latent_plots.plot_latent_scatter
-    bayesian_gplvm_minibatch.BayesianGPLVMMiniBatch.plot_inducing = gpy_plot.latent_plots.plot_latent_inducing
-    bayesian_gplvm_minibatch.BayesianGPLVMMiniBatch.plot_steepest_gradient_map = gpy_plot.latent_plots.plot_steepest_gradient_map
-    SSGPLVM.plot_latent = gpy_plot.latent_plots.plot_latent
-    SSGPLVM.plot_scatter = gpy_plot.latent_plots.plot_latent_scatter
-    SSGPLVM.plot_inducing = gpy_plot.latent_plots.plot_latent_inducing
-    SSGPLVM.plot_steepest_gradient_map = gpy_plot.latent_plots.plot_steepest_gradient_map
-
-    from ..kern import Kern
-    Kern.plot_covariance = gpy_plot.kernel_plots.plot_covariance
-    def deprecate_plot(self, *args, **kwargs):
-        import warnings
-        warnings.warn(DeprecationWarning('Kern.plot is being deprecated and will not be available in the 1.0 release. Use Kern.plot_covariance instead'))
-        return self.plot_covariance(*args, **kwargs)
-    Kern.plot = deprecate_plot
-    Kern.plot_ARD = gpy_plot.kernel_plots.plot_ARD
-
-    from ..inference.optimization import Optimizer
-    Optimizer.plot = gpy_plot.inference_plots.plot_optimizer
-    # Variational plot!
+from ..util.config import config, NoOptionError
+try:
+    lib = config.get('plotting', 'library')
+    change_plotting_library(lib)
+except NoOptionError:
+    print("No plotting library was specified in config file. \n{}".format(error_suggestion))

GPy/plotting/gpy_plot/gp_plots.py

@@ -420,4 +420,4 @@ def _plot(self, canvas, plots, helper_data, helper_prediction, levels, plot_indu
 
         if helper_prediction[2] is not None:
             plots.update(_plot_samples(self, canvas, helper_data, helper_prediction, projection, "Samples"))
-        return plots
+        return plots

GPy/plotting/gpy_plot/kernel_plots.py

@@ -33,7 +33,7 @@ from .. import Tango
 from .plot_util import update_not_existing_kwargs, helper_for_plot_data
 from ...kern.src.kern import Kern, CombinationKernel
 
-def plot_ARD(kernel, filtering=None, legend=False, **kwargs):
+def plot_ARD(kernel, filtering=None, legend=False, canvas=None, **kwargs):
     """
     If an ARD kernel is present, plot a bar representation using matplotlib
 
@@ -62,7 +62,11 @@ def plot_ARD(kernel, filtering=None, legend=False, **kwargs):
 
     bars = []
     kwargs = update_not_existing_kwargs(kwargs, pl().defaults.ard)
-    canvas, kwargs = pl().new_canvas(xlim=(-.5, kernel._effective_input_dim-.5), xlabel='input dimension', ylabel='sensitivity', **kwargs)
+
+
+    if canvas is None:
+        canvas, kwargs = pl().new_canvas(xlim=(-.5, kernel._effective_input_dim-.5), xlabel='input dimension', ylabel='sensitivity', **kwargs)
+
     for i in range(ard_params.shape[0]):
         if parts[i].name in filtering:
             c = Tango.nextMedium()
@@ -96,7 +100,7 @@ def plot_covariance(kernel, x=None, label=None,
     """
     X = np.ones((2, kernel._effective_input_dim)) * [[-3], [3]]
     _, free_dims, Xgrid, xx, yy, _, _, resolution = helper_for_plot_data(kernel, X, plot_limits, visible_dims, None, resolution)
-    
+
     from numbers import Number
     if x is None:
         from ...kern.src.stationary import Stationary
@@ -104,7 +108,7 @@ def plot_covariance(kernel, x=None, label=None,
     elif isinstance(x, Number):
         x = np.ones((1, kernel._effective_input_dim))*x
     K = kernel.K(Xgrid, x)
-    
+
     if projection == '3d':
         xlabel = 'X[:,0]'
         ylabel = 'X[:,1]'
@@ -136,4 +140,4 @@ def plot_covariance(kernel, x=None, label=None,
         return pl().add_to_canvas(canvas, plots)
 
     else:
-        raise NotImplementedError("Cannot plot a kernel with more than two input dimensions")
+        raise NotImplementedError("Cannot plot a kernel with more than two input dimensions")

GPy/plotting/gpy_plot/latent_plots.py

@@ -147,6 +147,7 @@ def _plot_magnification(self, canvas, which_indices, Xgrid,
     def plot_function(x):
         Xtest_full = np.zeros((x.shape[0], Xgrid.shape[1]))
         Xtest_full[:, which_indices] = x
+
         mf = self.predict_magnification(Xtest_full, kern=kern, mean=mean, covariance=covariance)
         return mf.reshape(resolution, resolution).T
     imshow_kwargs = update_not_existing_kwargs(imshow_kwargs, pl().defaults.magnification)
@@ -215,7 +216,12 @@ def _plot_latent(self, canvas, which_indices, Xgrid,
     def plot_function(x):
         Xtest_full = np.zeros((x.shape[0], Xgrid.shape[1]))
         Xtest_full[:, which_indices] = x
-        mf = np.log(self.predict(Xtest_full, kern=kern)[1])
+        mf = self.predict(Xtest_full, kern=kern)[1]
+        if mf.shape[1]==self.output_dim:
+            mf = mf.sum(-1)
+        else:
+            mf *= self.output_dim
+        mf = np.log(mf)
         return mf.reshape(resolution, resolution).T
 
     imshow_kwargs = update_not_existing_kwargs(imshow_kwargs, pl().defaults.latent)

GPy/plotting/matplot_dep/__init__.py

@@ -18,4 +18,4 @@
 
 
 from .util import align_subplot_array, align_subplots, fewerXticks, removeRightTicks, removeUpperTicks
-from . import controllers, base_plots
+from . import controllers, base_plots

GPy/plotting/matplot_dep/controllers/__init__.py

@@ -1 +1 @@
-from .imshow_controller import ImshowController, ImAnnotateController
+from .imshow_controller import ImshowController, ImAnnotateController

GPy/plotting/matplot_dep/controllers/imshow_controller.py

@@ -72,4 +72,4 @@ class ImAnnotateController(ImshowController):
             text.set_x(x+xoffset)
             text.set_y(y+yoffset)
             text.set_text("{}".format(X[1][j, i]))
-        return view
+        return view

GPy/plotting/matplot_dep/defaults.py

@@ -72,4 +72,4 @@ latent = dict(aspect='auto', cmap='Greys', interpolation='bicubic')
 gradient = dict(aspect='auto', cmap='RdBu', interpolation='nearest', alpha=.7)
 magnification = dict(aspect='auto', cmap='Greys', interpolation='bicubic')
 latent_scatter = dict(s=40, linewidth=.2, edgecolor='k', alpha=.9)
-annotation = dict(fontdict=dict(family='sans-serif', weight='light', fontsize=9), zorder=.3, alpha=.7)
+annotation = dict(fontdict=dict(family='sans-serif', weight='light', fontsize=9), zorder=.3, alpha=.7)

GPy/plotting/matplot_dep/plot_definitions.py

@@ -42,10 +42,11 @@ class MatplotlibPlots(AbstractPlottingLibrary):
         super(MatplotlibPlots, self).__init__()
         self._defaults = defaults.__dict__
 
-    def figure(self, rows=1, cols=1, **kwargs):
-        fig = plt.figure(**kwargs)
+    def figure(self, rows=1, cols=1, gridspec_kwargs={}, tight_layout=True, **kwargs):
+        fig = plt.figure(tight_layout=tight_layout, **kwargs)
         fig.rows = rows
         fig.cols = cols
+        fig.gridspec = plt.GridSpec(rows, cols, **gridspec_kwargs)
         return fig
 
     def new_canvas(self, figure=None, row=1, col=1, projection='2d', xlabel=None, ylabel=None, zlabel=None, title=None, xlim=None, ylim=None, zlim=None, **kwargs):
@@ -56,7 +57,9 @@ class MatplotlibPlots(AbstractPlottingLibrary):
         if 'ax' in kwargs:
             ax = kwargs.pop('ax')
         else:
-            if 'num' in kwargs and 'figsize' in kwargs:
+            if figure is not None:
+                fig = figure
+            elif 'num' in kwargs and 'figsize' in kwargs:
                 fig = self.figure(num=kwargs.pop('num'), figsize=kwargs.pop('figsize'))
             elif 'num' in kwargs:
                 fig = self.figure(num=kwargs.pop('num'))
@@ -66,7 +69,7 @@ class MatplotlibPlots(AbstractPlottingLibrary):
                 fig = self.figure()
 
             #if hasattr(fig, 'rows') and hasattr(fig, 'cols'):
-            ax = fig.add_subplot(fig.rows, fig.cols, (col,row), projection=projection)
+            ax = fig.add_subplot(fig.gridspec[row-1, col-1], projection=projection)
 
         if xlim is not None: ax.set_xlim(xlim)
         if ylim is not None: ax.set_ylim(ylim)
@@ -79,7 +82,7 @@ class MatplotlibPlots(AbstractPlottingLibrary):
         return ax, kwargs
 
     def add_to_canvas(self, ax, plots, legend=False, title=None, **kwargs):
-        ax.autoscale_view()
+        #ax.autoscale_view()
         fontdict=dict(family='sans-serif', weight='light', size=9)
         if legend is True:
             ax.legend(*ax.get_legend_handles_labels())
@@ -89,9 +92,7 @@ class MatplotlibPlots(AbstractPlottingLibrary):
         if title is not None: ax.figure.suptitle(title)
         return ax
 
-    def show_canvas(self, ax, tight_layout=False, **kwargs):
-        if tight_layout:
-            ax.figure.tight_layout()
+    def show_canvas(self, ax):
         ax.figure.canvas.draw()
         return ax.figure
 

GPy/plotting/matplot_dep/ssgplvm.py

@@ -13,16 +13,16 @@ class SSGPLVM_plot(object):
         self.model = model
         self.imgsize= imgsize
         assert model.Y.shape[1] == imgsize[0]*imgsize[1]
-        
+
     def plot_inducing(self):
         fig1 = pylab.figure()
         mean = self.model.posterior.mean
         arr = mean.reshape(*(mean.shape[0],self.imgsize[1],self.imgsize[0]))
         plot_2D_images(fig1, arr)
         fig1.gca().set_title('The mean of inducing points')
-        
+
         fig2 = pylab.figure()
         covar = self.model.posterior.covariance
         plot_2D_images(fig2, covar)
         fig2.gca().set_title('The variance of inducing points')
-        
+

GPy/plotting/matplot_dep/util.py

@@ -116,4 +116,4 @@ def align_subplot_array(axes,xlim=None, ylim=None):
         if i<(M*(N-1)):
             ax.set_xticks([])
         else:
-            removeUpperTicks(ax)
+            removeUpperTicks(ax)

GPy/plotting/plotly_dep/defaults.py

@@ -73,4 +73,4 @@ latent = dict(colorscale='Greys', reversescale=True, zsmooth='best')
 gradient = dict(colorscale='RdBu', opacity=.7)
 magnification = dict(colorscale='Greys', zsmooth='best', reversescale=True)
 latent_scatter = dict(marker_kwargs=dict(size='5', opacity=.7))
-# annotation = dict(fontdict=dict(family='sans-serif', weight='light', fontsize=9), zorder=.3, alpha=.7)
+# annotation = dict(fontdict=dict(family='sans-serif', weight='light', fontsize=9), zorder=.3, alpha=.7)

GPy/plotting/plotly_dep/plot_definitions.py

@@ -254,7 +254,7 @@ class PlotlyPlots(AbstractPlottingLibrary):
                         font=dict(color='white' if np.abs(var) > 0.8 else 'black', size=10),
                         opacity=.5,
                         showarrow=False,
-                        hoverinfo='x'))
+                        ))
         return imshow, annotations
 
     def annotation_heatmap_interact(self, ax, plot_function, extent, label=None, resolution=15, imshow_kwargs=None, **annotation_kwargs):

GPy/testing/__init__.py

@@ -1,14 +1,9 @@
-# Copyright (c) 2014, Max Zwiessele
+# Copyright (c) 2014, Max Zwiessele, GPy Authors
 # Licensed under the BSD 3-clause license (see LICENSE.txt)
-"""
-
-MaxZ
-
-"""
 import unittest
 import sys
 
 def deepTest(reason):
     if reason:
         return lambda x:x
-    return unittest.skip("Not deep scanning, enable deepscan by adding 'deep' argument")
+    return unittest.skip("Not deep scanning, enable deepscan by adding 'deep' argument to unittest call")