Merge branch 'devel' of github.com:SheffieldML/GPy into devel

GPy/core/gp.py

@@ -395,7 +395,7 @@ class GP(Model):
         which_data_ycols='all', fixed_inputs=[],
         levels=20, samples=0, fignum=None, ax=None, resolution=None,
         plot_raw=False,
-        linecol=None,fillcol=None, Y_metadata=None, data_symbol='kx'):
+        linecol=None,fillcol=None, Y_metadata=None, data_symbol='kx', predict_kw=None):
         """
         Plot the posterior of the GP.
           - In one dimension, the function is plotted with a shaded region identifying two standard deviations.
@@ -444,7 +444,7 @@ class GP(Model):
                                      which_data_ycols, fixed_inputs,
                                      levels, samples, fignum, ax, resolution,
                                      plot_raw=plot_raw, Y_metadata=Y_metadata,
-                                     data_symbol=data_symbol, **kw)
+                                     data_symbol=data_symbol, predict_kw=predict_kw, **kw)
 
     def input_sensitivity(self, summarize=True):
         """

GPy/core/sparse_gp.py

@@ -132,12 +132,14 @@ class SparseGP(GP):
                 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 = Kxx[:,:,None] - np.tensordot(np.dot(np.atleast_3d(self.posterior.woodbury_inv).T, Kx).T, Kx, [1,0]).swapaxes(1,2)
+                    var = np.empty((Kxx.shape[0],Kxx.shape[1],self.posterior.woodbury_inv.shape[2]))
+                    for i in range(var.shape[1]):
+                        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)
+                    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]):
@@ -147,9 +149,9 @@ class SparseGP(GP):
             if self.mean_function is not None:
                 mu += self.mean_function.f(Xnew)
         else:
-            psi0_star = self.kern.psi0(self.Z, Xnew)
+            psi0_star = kern.psi0(self.Z, Xnew)
-            psi1_star = self.kern.psi1(self.Z, Xnew)
+            psi1_star = kern.psi1(self.Z, Xnew)
-            #psi2_star = self.kern.psi2(self.Z, Xnew) # Only possible if we get NxMxM psi2 out of the code.
+            #psi2_star = kern.psi2(self.Z, Xnew) # Only possible if we get NxMxM psi2 out of the code.
             la = self.posterior.woodbury_vector
             mu = np.dot(psi1_star, la) # TODO: dimensions?
             
@@ -161,7 +163,7 @@ class SparseGP(GP):
                 
             for i in range(Xnew.shape[0]):
                 _mu, _var = Xnew.mean.values[[i]], Xnew.variance.values[[i]]
-                psi2_star = self.kern.psi2(self.Z, NormalPosterior(_mu, _var))
+                psi2_star = kern.psi2(self.Z, NormalPosterior(_mu, _var))
                 tmp = (psi2_star[:, :] - psi1_star[[i]].T.dot(psi1_star[[i]]))
 
                 var_ = mdot(la.T, tmp, la)

GPy/core/verbose_optimization.py

@@ -141,6 +141,13 @@ class VerboseOptimization(object):
 
     def finish(self, opt):
         self.status = opt.status
+        if self.verbose and self.ipython_notebook:
+            if 'conv' in self.status.lower():
+                self.progress.bar_style = 'success'
+            elif self.iteration >= self.maxiters:
+                self.progress.bar_style = 'warning'
+            else:
+                self.progress.bar_style = 'danger'
 
     def __exit__(self, type, value, traceback):
         if self.verbose:

GPy/inference/latent_function_inference/exact_gaussian_inference.py

@@ -52,7 +52,7 @@ class ExactGaussianInference(LatentFunctionInference):
         K = kern.K(X)
 
         Ky = K.copy()
-        diag.add(Ky, likelihood.gaussian_variance(Y_metadata))
+        diag.add(Ky, likelihood.gaussian_variance(Y_metadata)+1e-8)
         Wi, LW, LWi, W_logdet = pdinv(Ky)
 
         alpha, _ = dpotrs(LW, YYT_factor, lower=1)

GPy/kern/_src/basis_funcs.py

@@ -136,3 +136,48 @@ class DomainKernel(LinearSlopeBasisFuncKernel):
     def _phi(self, X):
         phi = np.where((X>self.start)*(X<self.stop), 1, 0)
         return phi#((phi-self.start)/(self.stop-self.start))-.5
+
+class LogisticBasisFuncKernel(BasisFuncKernel):
+    def __init__(self, input_dim, centers, variance=1., slope=1., active_dims=None, ARD=False, ARD_slope=True, name='logistic'):
+        self.centers = np.atleast_2d(centers)
+        self.ARD_slope = ARD_slope
+        if self.ARD_slope:
+            self.slope = Param('slope', slope * np.ones(self.centers.size), Logexp())
+        else:
+            self.slope = Param('slope', slope, Logexp())
+        super(LogisticBasisFuncKernel, self).__init__(input_dim, variance, active_dims, ARD, name)
+        self.link_parameter(self.slope)
+    
+    @Cache_this(limit=3, ignore_args=())
+    def _phi(self, X):
+        import scipy as sp
+        phi = 1/(1+np.exp(-((X-self.centers)*self.slope)))
+        return np.where(np.isnan(phi), 0, phi)#((phi-self.start)/(self.stop-self.start))-.5
+
+    def parameters_changed(self):
+        BasisFuncKernel.parameters_changed(self)
+    
+    def update_gradients_full(self, dL_dK, X, X2=None):
+        super(LogisticBasisFuncKernel, self).update_gradients_full(dL_dK, X, X2)
+        if X2 is None or X is X2:
+            phi1 = self.phi(X)
+            if phi1.ndim != 2:
+                phi1 = phi1[:, None]
+            dphi1_dl = (phi1**2) * (np.exp(-((X-self.centers)*self.slope)) * (X-self.centers))
+            if self.ARD_slope:
+                self.slope.gradient = self.variance * 2 * np.einsum('ij,iq,jq->q', dL_dK, phi1, dphi1_dl)
+            else:
+                self.slope.gradient = self.variance * 2 * (dL_dK * phi1.dot(dphi1_dl.T)).sum()
+        else:
+            phi1 = self.phi(X)
+            phi2 = self.phi(X2)
+            if phi1.ndim != 2:
+                phi1 = phi1[:, None]
+                phi2 = phi2[:, None]
+            dphi1_dl = (phi1**2) * (np.exp(-((X-self.centers)*self.slope)) * (X-self.centers))
+            dphi2_dl = (phi2**2) * (np.exp(-((X2-self.centers)*self.slope)) * (X2-self.centers))
+            if self.ARD_slope:
+                self.slope.gradient = (self.variance * np.einsum('ij,iq,jq->q', dL_dK, phi1, dphi2_dl) + np.einsum('ij,iq,jq->q', dL_dK, phi2, dphi1_dl))
+            else:
+                self.slope.gradient = self.variance * (dL_dK * phi1.dot(dphi2_dl.T)).sum() + (dL_dK * phi2.dot(dphi1_dl.T)).sum()
+        self.slope.gradient = np.where(np.isnan(self.slope.gradient), 0, self.slope.gradient)

GPy/plotting/matplot_dep/models_plots.py

@@ -17,7 +17,7 @@ def plot_fit(model, plot_limits=None, which_data_rows='all',
         levels=20, samples=0, fignum=None, ax=None, resolution=None,
         plot_raw=False,
         linecol=Tango.colorsHex['darkBlue'],fillcol=Tango.colorsHex['lightBlue'], Y_metadata=None, data_symbol='kx',
-        apply_link=False, samples_f=0, plot_uncertain_inputs=True):
+        apply_link=False, samples_f=0, plot_uncertain_inputs=True, predict_kw=None):
     """
     Plot the posterior of the GP.
       - In one dimension, the function is plotted with a shaded region identifying two standard deviations.
@@ -76,6 +76,9 @@ def plot_fit(model, plot_limits=None, which_data_rows='all',
 
     if hasattr(model, 'Z'): Z = model.Z
 
+    if predict_kw is None:
+        predict_kw = {}
+    
     #work out what the inputs are for plotting (1D or 2D)
     fixed_dims = np.array([i for i,v in fixed_inputs])
     free_dims = np.setdiff1d(np.arange(model.input_dim),fixed_dims)
@@ -92,7 +95,7 @@ def plot_fit(model, plot_limits=None, which_data_rows='all',
 
         #make a prediction on the frame and plot it
         if plot_raw:
-            m, v = model._raw_predict(Xgrid)
+            m, v = model._raw_predict(Xgrid, **predict_kw)
             if apply_link:
                 lower = model.likelihood.gp_link.transf(m - 2*np.sqrt(v))
                 upper = model.likelihood.gp_link.transf(m + 2*np.sqrt(v))
@@ -106,7 +109,7 @@ def plot_fit(model, plot_limits=None, which_data_rows='all',
                 meta = {'output_index': Xgrid[:,-1:].astype(np.int)}
             else:
                 meta = None
-            m, v = model.predict(Xgrid, full_cov=False, Y_metadata=meta)
+            m, v = model.predict(Xgrid, full_cov=False, Y_metadata=meta, **predict_kw)
             lower, upper = model.predict_quantiles(Xgrid, Y_metadata=meta)
 
 
@@ -178,13 +181,13 @@ def plot_fit(model, plot_limits=None, which_data_rows='all',
 
         #predict on the frame and plot
         if plot_raw:
-            m, _ = model._raw_predict(Xgrid)
+            m, _ = model._raw_predict(Xgrid, **predict_kw)
         else:
             if isinstance(model,GPCoregionalizedRegression) or isinstance(model,SparseGPCoregionalizedRegression):
                 meta = {'output_index': Xgrid[:,-1:].astype(np.int)}
             else:
                 meta = None
-            m, v = model.predict(Xgrid, full_cov=False, Y_metadata=meta)
+            m, v = model.predict(Xgrid, full_cov=False, Y_metadata=meta, **predict_kw)
         for d in which_data_ycols:
             m_d = m[:,d].reshape(resolution, resolution).T
             plots['contour'] = ax.contour(x, y, m_d, levels, vmin=m.min(), vmax=m.max(), cmap=pb.cm.jet)