diff --git a/.travis.yml b/.travis.yml
index 1f796285..e9ae5831 100644
--- a/.travis.yml
+++ b/.travis.yml
@@ -24,3 +24,4 @@ install:
 # command to run tests, e.g. python setup.py test
 script:
   - nosetests GPy/testing
+  #- yes | nosetests GPy/testing
diff --git a/GPy/__init__.py b/GPy/__init__.py
index 5c4838da..2ccd109d 100644
--- a/GPy/__init__.py
+++ b/GPy/__init__.py
@@ -4,17 +4,17 @@ import warnings
 warnings.filterwarnings("ignore", category=DeprecationWarning)
 import os
 
-import core
-import models
-import mappings
-import inference
 import util
-import examples
+import core
+import kern
+import mappings
 import likelihoods
+import inference
+import models
+import examples
 import testing
 from numpy.testing import Tester
 from nose.tools import nottest
-import kern
 from core import priors
 
 @nottest
diff --git a/GPy/core/gp_base.py b/GPy/core/gp_base.py
index 2577e06c..7d58c82c 100644
--- a/GPy/core/gp_base.py
+++ b/GPy/core/gp_base.py
@@ -218,8 +218,8 @@ class GPBase(Model):
                 Y = self.likelihood.data
             for d in which_data_ycols:
                 m_d = m[:,d].reshape(resolution, resolution).T
-                ax.contour(x, y, m_d, levels, vmin=m.min(), vmax=m.max(), cmap=pb.cm.jet)
-                ax.scatter(self.X[which_data_rows, free_dims[0]], self.X[which_data_rows, free_dims[1]], 40, Y[which_data_rows, d], cmap=pb.cm.jet, vmin=m.min(), vmax=m.max(), linewidth=0.)
+                contour = ax.contour(x, y, m_d, levels, vmin=m.min(), vmax=m.max(), cmap=pb.cm.jet)
+                scatter = ax.scatter(self.X[which_data_rows, free_dims[0]], self.X[which_data_rows, free_dims[1]], 40, Y[which_data_rows, d], cmap=pb.cm.jet, vmin=m.min(), vmax=m.max(), linewidth=0.)
 
             #set the limits of the plot to some sensible values
             ax.set_xlim(xmin[0], xmax[0])
@@ -227,7 +227,7 @@ class GPBase(Model):
 
             if samples:
                 warnings.warn("Samples are rather difficult to plot for 2D inputs...")
-
+            return contour, scatter
         else:
             raise NotImplementedError, "Cannot define a frame with more than two input dimensions"
 
diff --git a/GPy/core/sparse_gp.py b/GPy/core/sparse_gp.py
index 5e381110..43af97aa 100644
--- a/GPy/core/sparse_gp.py
+++ b/GPy/core/sparse_gp.py
@@ -381,7 +381,7 @@ class SparseGP(GPBase):
             which_data_ycols='all', which_parts='all', fixed_inputs=[],
             plot_raw=False,
             levels=20, samples=0, fignum=None, ax=None, resolution=None):
-        """ 
+        """
         Plot the posterior of the sparse GP.
           - In one dimension, the function is plotted with a shaded region identifying two standard deviations.
           - In two dimsensions, a contour-plot shows the mean predicted function
@@ -417,6 +417,11 @@ class SparseGP(GPBase):
         :param levels: for 2D plotting, the number of contour levels to use is ax is None, create a new figure
         """
         #deal work out which ax to plot on
+        #Need these because we use which_data_rows in this function not just base
+        if which_data_rows == 'all':
+            which_data_rows = slice(None)
+        if which_data_ycols == 'all':
+            which_data_ycols = np.arange(self.output_dim)
         if ax is None:
             fig = pb.figure(num=fignum)
             ax = fig.add_subplot(111)
diff --git a/GPy/core/svigp.py b/GPy/core/svigp.py
index fdd95aa8..d47d46f1 100644
--- a/GPy/core/svigp.py
+++ b/GPy/core/svigp.py
@@ -37,7 +37,6 @@ class SVIGP(GPBase):
         self.Y = self.likelihood.Y.copy()
         self.Z = Z
         self.num_inducing = Z.shape[0]
-
         self.batchcounter = 0
         self.epochs = 0
         self.iterations = 0
@@ -78,6 +77,8 @@ class SVIGP(GPBase):
         self._param_steplength_trace = []
         self._vb_steplength_trace = []
 
+        self.ensure_default_constraints()
+
     def getstate(self):
         steplength_params = [self.hbar_t, self.tau_t, self.gbar_t, self.gbar_t1, self.gbar_t2, self.hbar_tp, self.tau_tp, self.gbar_tp, self.adapt_param_steplength, self.adapt_vb_steplength, self.vb_steplength, self.param_steplength]
         return GPBase.getstate(self) + \
@@ -303,12 +304,12 @@ class SVIGP(GPBase):
 
         #Iterate!
         for i in range(iterations):
-
+            
             #store the current configuration for plotting later
             self._param_trace.append(self._get_params())
             self._ll_trace.append(self.log_likelihood() + self.log_prior())
 
-            #load a batch
+            #load a batch and do the appropriate computations (kernel matrices, etc)
             self.load_batch()
 
             #compute the (stochastic) gradient
@@ -318,7 +319,8 @@ class SVIGP(GPBase):
 
             #compute the steps in all parameters
             vb_step = self.vb_steplength*natgrads[0]
-            if (self.epochs>=1):#only move the parameters after the first epoch
+            #only move the parameters after the first epoch and only if the steplength is nonzero
+            if (self.epochs>=1) and (self.param_steplength > 0):
                 param_step = self.momentum*param_step + self.param_steplength*grads
             else:
                 param_step = 0.
@@ -340,6 +342,8 @@ class SVIGP(GPBase):
 
             if self.epochs > 10:
                 self._adapt_steplength()
+            self._vb_steplength_trace.append(self.vb_steplength)
+            self._param_steplength_trace.append(self.param_steplength)
 
             self.iterations += 1
 
@@ -348,17 +352,20 @@ class SVIGP(GPBase):
         if self.adapt_vb_steplength:
             # self._adaptive_vb_steplength()
             self._adaptive_vb_steplength_KL()
-        self._vb_steplength_trace.append(self.vb_steplength)
-        assert self.vb_steplength > 0
+        #self._vb_steplength_trace.append(self.vb_steplength)
+        assert self.vb_steplength >= 0
 
         if self.adapt_param_steplength:
             self._adaptive_param_steplength()
             # self._adaptive_param_steplength_log()
             # self._adaptive_param_steplength_from_vb()
-        self._param_steplength_trace.append(self.param_steplength)
+        #self._param_steplength_trace.append(self.param_steplength)
 
     def _adaptive_param_steplength(self):
-        decr_factor = 0.02
+        if hasattr(self, 'adapt_param_steplength_decr'):
+            decr_factor = self.adapt_param_steplength_decr
+        else:
+            decr_factor = 0.02
         g_tp = self._transform_gradients(self._log_likelihood_gradients())
         self.gbar_tp = (1-1/self.tau_tp)*self.gbar_tp + 1/self.tau_tp * g_tp
         self.hbar_tp = (1-1/self.tau_tp)*self.hbar_tp + 1/self.tau_tp * np.dot(g_tp.T, g_tp)
diff --git a/GPy/core/variational.py b/GPy/core/variational.py
deleted file mode 100644
index 74287dcf..00000000
--- a/GPy/core/variational.py
+++ /dev/null
@@ -1,19 +0,0 @@
-'''
-Created on 6 Nov 2013
-
-@author: maxz
-'''
-from parameterized import Parameterized
-from parameter import Param
-
-class Normal(Parameterized):
-    '''
-    Normal distribution for variational approximations.
-    
-    holds the means and variances for a factorizing multivariate normal distribution
-    '''
-    def __init__(self, name, means, variances):
-        Parameterized.__init__(self, name=name)
-        self.means = Param("mean", means)
-        self.variances = Param('variance', variances)
-        self.add_parameters(self.means, self.variances)
\ No newline at end of file
diff --git a/GPy/examples/dimensionality_reduction.py b/GPy/examples/dimensionality_reduction.py
index 94bb4955..46fc6797 100644
--- a/GPy/examples/dimensionality_reduction.py
+++ b/GPy/examples/dimensionality_reduction.py
@@ -37,20 +37,32 @@ def bgplvm_test_model(seed=default_seed, optimize=False, verbose=1, plot=False):
     # k = GPy.kern.rbf(input_dim, .5, _np.ones(input_dim) * 2., ARD=True) + GPy.kern.linear(input_dim, _np.ones(input_dim) * .2, ARD=True)
 
     m = GPy.models.BayesianGPLVM(lik, input_dim, kernel=k, num_inducing=num_inducing)
+    #===========================================================================
+    # randomly obstruct data with percentage p
+    p = .8
+    Y_obstruct = Y.copy()
+    Y_obstruct[_np.random.uniform(size=(Y.shape)) < p] = _np.nan
+    #===========================================================================
+    m2 = GPy.models.BayesianGPLVMWithMissingData(Y_obstruct, input_dim, kernel=k, num_inducing=num_inducing)
     m.lengthscales = lengthscales
 
     if plot:
         import matplotlib.pyplot as pb
         m.plot()
         pb.title('PCA initialisation')
+        m2.plot()
+        pb.title('PCA initialisation')
 
     if optimize:
         m.optimize('scg', messages=verbose)
+        m2.optimize('scg', messages=verbose)
         if plot:
             m.plot()
             pb.title('After optimisation')
+            m2.plot()
+            pb.title('After optimisation')
 
-    return m
+    return m, m2
 
 def gplvm_oil_100(optimize=True, verbose=1, plot=True):
     import GPy
@@ -217,7 +229,7 @@ def _simulate_sincos(D1, D2, D3, N, num_inducing, Q, plot_sim=False):
         ax.legend()
         for i, Y in enumerate(Ylist):
             ax = fig.add_subplot(2, len(Ylist), len(Ylist) + 1 + i)
-            ax.imshow(Y, aspect='auto', cmap=cm.gray)
+            ax.imshow(Y, aspect='auto', cmap=cm.gray) # @UndefinedVariable
             ax.set_title("Y{}".format(i + 1))
         pylab.draw()
         pylab.tight_layout()
@@ -451,12 +463,9 @@ def cmu_mocap(subject='35', motion=['01'], in_place=True, optimize=True, verbose
     if in_place:
         # Make figure move in place.
         data['Y'][:, 0:3] = 0.0
-
     m = GPy.models.GPLVM(data['Y'], 2, normalize_Y=True)
 
-    if optimize:
-        m.optimize(messages=verbose, max_f_eval=10000)
-
+    if optimize: m.optimize(messages=verbose, max_f_eval=10000)
     if plot:
         ax = m.plot_latent()
         y = m.likelihood.Y[0, :]
diff --git a/GPy/examples/regression.py b/GPy/examples/regression.py
index 9b910005..65a50f0e 100644
--- a/GPy/examples/regression.py
+++ b/GPy/examples/regression.py
@@ -273,27 +273,6 @@ def toy_rbf_1d_50(optimize=True, plot=True):
 
     return m
 
-
-def toy_poisson_rbf_1d(optimize=True, plot=True):
-    """Run a simple demonstration of a standard Gaussian process fitting it to data sampled from an RBF covariance."""
-    x_len = 400
-    X = np.linspace(0, 10, x_len)[:, None]
-    f_true = np.random.multivariate_normal(np.zeros(x_len), GPy.kern.rbf(1).K(X))
-    Y = np.array([np.random.poisson(np.exp(f)) for f in f_true]).reshape(x_len,1)
-
-    noise_model = GPy.likelihoods.poisson()
-    likelihood = GPy.likelihoods.EP(Y,noise_model)
-
-    # create simple GP Model
-    m = GPy.models.GPRegression(X, Y, likelihood=likelihood)
-
-    if optimize:
-        m.optimize('bfgs')
-    if plot:
-        m.plot()
-
-    return m
-
 def toy_poisson_rbf_1d_laplace(optimize=True, plot=True):
     """Run a simple demonstration of a standard Gaussian process fitting it to data sampled from an RBF covariance."""
     optimizer='scg'
diff --git a/GPy/kern/parts/ODE_UY.py b/GPy/kern/parts/ODE_UY.py
index bb736cc5..f251d8e2 100644
--- a/GPy/kern/parts/ODE_UY.py
+++ b/GPy/kern/parts/ODE_UY.py
@@ -7,7 +7,7 @@ import numpy as np
 
 def index_to_slices(index):
     """
-    take a numpy array of integers (index) and return a  nested list of slices such that the slices describe the start, stop points for each integer in the index. 
+    take a numpy array of integers (index) and return a  nested list of slices such that the slices describe the start, stop points for each integer in the index.
 
     e.g.
     >>> index = np.asarray([0,0,0,1,1,1,2,2,2])
@@ -38,7 +38,7 @@ class ODE_UY(Kernpart):
     :param input_dim: the number of input dimension, has to be equal to one
     :type input_dim: int
     :param input_lengthU: the number of input U length
-    :type input_dim: int   
+    :type input_dim: int
     :param varianceU: variance of the driving GP
     :type varianceU: float
     :param lengthscaleU: lengthscale of the driving GP  (sqrt(3)/lengthscaleU)
@@ -96,7 +96,7 @@ class ODE_UY(Kernpart):
     def K(self, X, X2, target):
         """Compute the covariance matrix between X and X2."""
         # model :   a * dy/dt + b * y = U
-        #lu=sqrt(3)/theta1  ly=1/theta2  theta2= a/b :thetay   sigma2=1/(2ab) :sigmay   
+        #lu=sqrt(3)/theta1  ly=1/theta2  theta2= a/b :thetay   sigma2=1/(2ab) :sigmay
 
         X,slices = X[:,:-1],index_to_slices(X[:,-1])
         if X2 is None:
@@ -110,24 +110,31 @@ class ODE_UY(Kernpart):
         ly=1/self.lengthscaleY
         lu=np.sqrt(3)/self.lengthscaleU
         #iu=self.input_lengthU  #dimention of U
-        
+
         Vu=self.varianceU
         Vy=self.varianceY
 
+        # kernel for kuu  matern3/2
         kuu = lambda dist:Vu * (1 + lu* np.abs(dist)) * np.exp(-lu * np.abs(dist))
 
+        # kernel for kyy
         k1 = lambda dist:np.exp(-ly*np.abs(dist))*(2*lu+ly)/(lu+ly)**2
-        k2 = lambda dist:(np.exp(-lu*dist)*(ly-2*lu+lu*ly*dist-lu**2*dist) + np.exp(-ly*dist)*(2*lu-ly) ) / (ly-lu)**2 
+        k2 = lambda dist:(np.exp(-lu*dist)*(ly-2*lu+lu*ly*dist-lu**2*dist) + np.exp(-ly*dist)*(2*lu-ly) ) / (ly-lu)**2
         k3 = lambda dist:np.exp(-lu*dist) * ( (1+lu*dist)/(lu+ly) + (lu)/(lu+ly)**2 )
         kyy = lambda dist:Vu*Vy*(k1(dist) + k2(dist) + k3(dist))
 
+
+        # cross covariance function
         kyu3 = lambda dist:np.exp(-lu*dist)/(lu+ly)*(1+lu*(dist+1/(lu+ly)))
+
+        # cross covariance kyu
         kyup = lambda dist:Vu*Vy*(k1(dist)+k2(dist))    #t>0 kyu
         kyun = lambda dist:Vu*Vy*(kyu3(dist))       #t<0 kyu
 
+        # cross covariance kuy
         kuyp = lambda dist:Vu*Vy*(kyu3(dist))       #t>0 kuy
         kuyn = lambda dist:Vu*Vy*(k1(dist)+k2(dist))      #t<0 kuy
-        
+
         for i, s1 in enumerate(slices):
             for j, s2 in enumerate(slices2):
                 for ss1 in s1:
@@ -135,12 +142,13 @@ class ODE_UY(Kernpart):
                         if i==0 and j==0:
                             target[ss1,ss2] = kuu(np.abs(rdist[ss1,ss2]))
                         elif i==0 and j==1:
-                            target[ss1,ss2] = np.where(  rdist[ss1,ss2]>0 , kuyp(np.abs(rdist[ss1,ss2])), kuyn(np.abs(rdist[s1[0],s2[0]]) )   )
+                            #target[ss1,ss2] = np.where(  rdist[ss1,ss2]>0 , kuyp(np.abs(rdist[ss1,ss2])), kuyn(np.abs(rdist[s1[0],s2[0]]) )   )
+                            target[ss1,ss2] = np.where(  rdist[ss1,ss2]>0 , kuyp(np.abs(rdist[ss1,ss2])), kuyn(np.abs(rdist[ss1,ss2]) )   )
                         elif i==1 and j==1:
                             target[ss1,ss2] = kyy(np.abs(rdist[ss1,ss2]))
                         else:
-                            target[ss1,ss2] = np.where(  rdist[ss1,ss2]>0 , kyup(np.abs(rdist[ss1,ss2])), kyun(np.abs(rdist[s1[0],s2[0]]) )   )
-
+                            #target[ss1,ss2] = np.where(  rdist[ss1,ss2]>0 , kyup(np.abs(rdist[ss1,ss2])), kyun(np.abs(rdist[s1[0],s2[0]]) )   )
+                            target[ss1,ss2] = np.where(  rdist[ss1,ss2]>0 , kyup(np.abs(rdist[ss1,ss2])), kyun(np.abs(rdist[ss1,ss2]) )   )
 
         #KUU = kuu(np.abs(rdist[:iu,:iu]))
 
@@ -160,22 +168,22 @@ class ODE_UY(Kernpart):
         lu=np.sqrt(3)/self.lengthscaleU
         #ly=self.lengthscaleY
         #lu=self.lengthscaleU
-        
+
         k1 = (2*lu+ly)/(lu+ly)**2
-        k2 = (ly-2*lu + 2*lu-ly ) / (ly-lu)**2 
-        k3 = 1/(lu+ly) + (lu)/(lu+ly)**2 
+        k2 = (ly-2*lu + 2*lu-ly ) / (ly-lu)**2
+        k3 = 1/(lu+ly) + (lu)/(lu+ly)**2
 
         slices = index_to_slices(X[:,-1])
 
         for i, ss1 in enumerate(slices):
             for s1 in ss1:
                 if i==0:
-                    target[s1]+= self.varianceU 
+                    target[s1]+= self.varianceU
                 elif i==1:
                     target[s1]+= self.varianceU*self.varianceY*(k1+k2+k3)
                 else:
                     raise ValueError, "invalid input/output index"
-        
+
         #target[slices[0][0]]+= self.varianceU   #matern32 diag
         #target[slices[1][0]]+= self.varianceU*self.varianceY*(k1+k2+k3)  #  diag
 
@@ -186,20 +194,30 @@ class ODE_UY(Kernpart):
 
     def dK_dtheta(self, dL_dK, X, X2, target):
         """derivative of the covariance matrix with respect to the parameters."""
-        if X2 is None: X2 = X
-        dist = np.abs(X - X2.T)
 
         X,slices = X[:,:-1],index_to_slices(X[:,-1])
         if X2 is None:
             X2,slices2 = X,slices
         else:
             X2,slices2 = X2[:,:-1],index_to_slices(X2[:,-1])
-
-
+        #rdist = X[:,0][:,None] - X2[:,0][:,None].T
+        rdist = X - X2.T
         ly=1/self.lengthscaleY
         lu=np.sqrt(3)/self.lengthscaleU
-        #ly=self.lengthscaleY
-        #lu=self.lengthscaleU
+
+        rd=rdist.shape[0]
+        dktheta1 = np.zeros([rd,rd])
+        dktheta2 = np.zeros([rd,rd])
+        dkUdvar = np.zeros([rd,rd])
+        dkYdvar = np.zeros([rd,rd])
+
+        # dk dtheta for UU
+        UUdtheta1 = lambda dist: np.exp(-lu* dist)*dist + (-dist)*np.exp(-lu* dist)*(1+lu*dist)
+        UUdtheta2 = lambda dist: 0
+        #UUdvar = lambda dist: (1 + lu*dist)*np.exp(-lu*dist)
+        UUdvar = lambda dist: (1 + lu* np.abs(dist)) * np.exp(-lu * np.abs(dist))
+
+        # dk dtheta for YY
 
         dk1theta1 = lambda dist: np.exp(-ly*dist)*2*(-lu)/(lu+ly)**3
         #c=np.sqrt(3)
@@ -207,16 +225,16 @@ class ODE_UY(Kernpart):
         #t2=1/ly
         #dk1theta1=np.exp(-dist*ly)*t2*( (2*c*t2+2*t1)/(c*t2+t1)**2 -2*(2*c*t2*t1+t1**2)/(c*t2+t1)**3   )
 
-        dk2theta1 = lambda dist: 1*( 
-            np.exp(-lu*dist)*dist*(-ly+2*lu-lu*ly*dist+dist*lu**2)*(ly-lu)**(-2) + np.exp(-lu*dist)*(-2+ly*dist-2*dist*lu)*(ly-lu)**(-2) 
-            +np.exp(-dist*lu)*(ly-2*lu+ly*lu*dist-dist*lu**2)*2*(ly-lu)**(-3) 
+        dk2theta1 = lambda dist: 1*(
+            np.exp(-lu*dist)*dist*(-ly+2*lu-lu*ly*dist+dist*lu**2)*(ly-lu)**(-2) + np.exp(-lu*dist)*(-2+ly*dist-2*dist*lu)*(ly-lu)**(-2)
+            +np.exp(-dist*lu)*(ly-2*lu+ly*lu*dist-dist*lu**2)*2*(ly-lu)**(-3)
             +np.exp(-dist*ly)*2*(ly-lu)**(-2)
             +np.exp(-dist*ly)*2*(2*lu-ly)*(ly-lu)**(-3)
             )
-      
+
         dk3theta1 = lambda dist: np.exp(-dist*lu)*(lu+ly)**(-2)*((2*lu+ly+dist*lu**2+lu*ly*dist)*(-dist-2/(lu+ly))+2+2*lu*dist+ly*dist)
 
-        dktheta1 = lambda dist: self.varianceU*self.varianceY*(dk1theta1+dk2theta1+dk3theta1)
+        #dktheta1 = lambda dist: self.varianceU*self.varianceY*(dk1theta1+dk2theta1+dk3theta1)
 
 
 
@@ -230,14 +248,35 @@ class ODE_UY(Kernpart):
 
         dk3theta2 = lambda dist: np.exp(-dist*lu) * (-3*lu-ly-dist*lu**2-lu*ly*dist)/(lu+ly)**3
 
-        dktheta2 = lambda dist: self.varianceU*self.varianceY*(dk1theta2 + dk2theta2 +dk3theta2)
-
-
+        #dktheta2 = lambda dist: self.varianceU*self.varianceY*(dk1theta2 + dk2theta2 +dk3theta2)
 
+        # kyy kernel
+        #k1 = lambda dist: np.exp(-ly*dist)*(2*lu+ly)/(lu+ly)**2
+        #k2 = lambda dist: (np.exp(-lu*dist)*(ly-2*lu+lu*ly*dist-lu**2*dist) + np.exp(-ly*dist)*(2*lu-ly) ) / (ly-lu)**2
+        #k3 = lambda dist: np.exp(-lu*dist) * ( (1+lu*dist)/(lu+ly) + (lu)/(lu+ly)**2 )
         k1 = lambda dist: np.exp(-ly*dist)*(2*lu+ly)/(lu+ly)**2
-        k2 = lambda dist: (np.exp(-lu*dist)*(ly-2*lu+lu*ly*dist-lu**2*dist) + np.exp(-ly*dist)*(2*lu-ly) ) / (ly-lu)**2 
+        k2 = lambda dist: (np.exp(-lu*dist)*(ly-2*lu+lu*ly*dist-lu**2*dist) + np.exp(-ly*dist)*(2*lu-ly) ) / (ly-lu)**2
         k3 = lambda dist: np.exp(-lu*dist) * ( (1+lu*dist)/(lu+ly) + (lu)/(lu+ly)**2 )
-        dkdvar = k1+k2+k3
+        #dkdvar = k1+k2+k3
+
+        #cross covariance kernel
+        kyu3 = lambda dist:np.exp(-lu*dist)/(lu+ly)*(1+lu*(dist+1/(lu+ly)))
+
+        # dk dtheta for UY
+        dkcrtheta2 = lambda dist: np.exp(-lu*dist) * ( (-1)*(lu+ly)**(-2)*(1+lu*dist+lu*(lu+ly)**(-1)) + (lu+ly)**(-1)*(-lu)*(lu+ly)**(-2) )
+        dkcrtheta1 = lambda dist: np.exp(-lu*dist)*(lu+ly)**(-1)* ( (-dist)*(1+dist*lu+lu*(lu+ly)**(-1)) - (lu+ly)**(-1)*(1+dist*lu+lu*(lu+ly)**(-1)) +dist+(lu+ly)**(-1)-lu*(lu+ly)**(-2) )
+        #dkuyp dtheta
+        #dkuyp dtheta1 = self.varianceU*self.varianceY* (dk1theta1() + dk2theta1())
+        #dkuyp dtheta2 = self.varianceU*self.varianceY* (dk1theta2() + dk2theta2())
+        #dkuyp dVar = k1() + k2()
+
+
+        #dkyup dtheta
+        #dkyun dtheta1 = self.varianceU*self.varianceY* (dk1theta1() + dk2theta1())
+        #dkyun dtheta2 = self.varianceU*self.varianceY* (dk1theta2() + dk2theta2())
+        #dkyup dVar = k1() + k2()        #
+
+
 
 
         for i, s1 in enumerate(slices):
@@ -246,20 +285,35 @@ class ODE_UY(Kernpart):
                     for ss2 in s2:
                         if i==0 and j==0:
                             #target[ss1,ss2] = kuu(np.abs(rdist[ss1,ss2]))
+                            dktheta1[ss1,ss2] = self.varianceU*self.varianceY*UUdtheta1(np.abs(rdist[ss1,ss2]))
+                            dktheta2[ss1,ss2] = 0
+                            dkUdvar[ss1,ss2] = UUdvar(np.abs(rdist[ss1,ss2]))
+                            dkYdvar[ss1,ss2] = 0
                         elif i==0 and j==1:
                             #target[ss1,ss2] = np.where(  rdist[ss1,ss2]>0 , kuyp(np.abs(rdist[ss1,ss2])), kuyn(np.abs(rdist[s1[0],s2[0]]) )   )
+                            #dktheta1[ss1,ss2] =
+                            #dktheta2[ss1,ss2] =
+                            #dkdvar[ss1,ss2] =           np.where(  rdist[ss1,ss2]>0 , kuyp(np.abs(rdist[ss1,ss2])), kuyn(np.abs(rdist[s1[0],s2[0]]) )   )
+                            dktheta1[ss1,ss2] = np.where(  rdist[ss1,ss2]>0 , dkcrtheta1(np.abs(rdist[ss1,ss2])) ,self.varianceU*self.varianceY*(dk1theta1(np.abs(rdist[ss1,ss2]))+dk2theta1(np.abs(rdist[ss1,ss2])))    )
+                            dktheta2[ss1,ss2] = np.where(  rdist[ss1,ss2]>0 , dkcrtheta2(np.abs(rdist[ss1,ss2])) ,self.varianceU*self.varianceY*(dk1theta2(np.abs(rdist[ss1,ss2]))+dk2theta2(np.abs(rdist[ss1,ss2])))    )
+                            dkUdvar[ss1,ss2] = np.where(  rdist[ss1,ss2]>0 ,  kyu3(np.abs(rdist[ss1,ss2]))  ,k1(np.abs(rdist[ss1,ss2]))+k2(np.abs(rdist[ss1,ss2]))  )
+                            dkYdvar[ss1,ss2] = dkUdvar[ss1,ss2]
                         elif i==1 and j==1:
                             #target[ss1,ss2] = kyy(np.abs(rdist[ss1,ss2]))
+                            dktheta1[ss1,ss2] = self.varianceU*self.varianceY*(dk1theta1(np.abs(rdist[ss1,ss2]))+dk2theta1(np.abs(rdist[ss1,ss2]))+dk3theta1(np.abs(rdist[ss1,ss2])))
+                            dktheta2[ss1,ss2] = self.varianceU*self.varianceY*(dk1theta2(np.abs(rdist[ss1,ss2])) + dk2theta2(np.abs(rdist[ss1,ss2])) +dk3theta2(np.abs(rdist[ss1,ss2])))
+                            dkUdvar[ss1,ss2] = (k1(np.abs(rdist[ss1,ss2]))+k2(np.abs(rdist[ss1,ss2]))+k3(np.abs(rdist[ss1,ss2])) )
+                            dkYdvar[ss1,ss2] = dkUdvar[ss1,ss2]
                         else:
                             #target[ss1,ss2] = np.where(  rdist[ss1,ss2]>0 , kyup(np.abs(rdist[ss1,ss2])), kyun(np.abs(rdist[s1[0],s2[0]]) )   )
+                            dktheta1[ss1,ss2] = np.where(  rdist[ss1,ss2]>0 ,self.varianceU*self.varianceY*(dk1theta1(np.abs(rdist[ss1,ss2]))+dk2theta1(np.abs(rdist[ss1,ss2]))) , dkcrtheta1(np.abs(rdist[ss1,ss2])) )
+                            dktheta2[ss1,ss2] = np.where(  rdist[ss1,ss2]>0 ,self.varianceU*self.varianceY*(dk1theta2(np.abs(rdist[ss1,ss2]))+dk2theta2(np.abs(rdist[ss1,ss2]))) , dkcrtheta2(np.abs(rdist[ss1,ss2])) )
+                            dkUdvar[ss1,ss2] = np.where(  rdist[ss1,ss2]>0 , k1(np.abs(rdist[ss1,ss2]))+k2(np.abs(rdist[ss1,ss2])), kyu3(np.abs(rdist[ss1,ss2])) )
+                            dkYdvar[ss1,ss2] = dkUdvar[ss1,ss2]
 
 
-
-
-
-
-        target[0] += np.sum(self.varianceY*dkdvar * dL_dK)
-        target[1] += np.sum(self.varianceU*dkdvar * dL_dK)
+        target[0] += np.sum(self.varianceY*dkUdvar * dL_dK)
+        target[1] += np.sum(self.varianceU*dkYdvar * dL_dK)
         target[2] += np.sum(dktheta1*(-np.sqrt(3)*self.lengthscaleU**(-2)) * dL_dK)
         target[3] += np.sum(dktheta2*(-self.lengthscaleY**(-2)) * dL_dK)
 
diff --git a/GPy/kern/parts/sympykern.py b/GPy/kern/parts/sympykern.py
index bcd52fe2..a839437b 100644
--- a/GPy/kern/parts/sympykern.py
+++ b/GPy/kern/parts/sympykern.py
@@ -6,7 +6,7 @@ from scipy import weave
 import re
 import os
 import sys
-current_dir = os.path.dirname(os.path.abspath(os.path.dirname(__file__)))
+current_dir = os.path.dirname(os.path.abspath(__file__))
 import tempfile
 import pdb
 import ast
@@ -107,9 +107,9 @@ class spkern(Kernpart):
             
         self.weave_kwargs = {
             'support_code':self._function_code,
-            'include_dirs':[tempfile.gettempdir(), os.path.join(current_dir,'parts/')],
+            'include_dirs':[tempfile.gettempdir(), current_dir],
             'headers':['"sympy_helpers.h"'],
-            'sources':[os.path.join(current_dir,"parts/sympy_helpers.cpp")],
+            'sources':[os.path.join(current_dir,"sympy_helpers.cpp")],
             'extra_compile_args':extra_compile_args,
             'extra_link_args':[],
             'verbose':True}
diff --git a/GPy/likelihoods/laplace.py b/GPy/likelihoods/laplace.py
index 0def0c8b..76bfc629 100644
--- a/GPy/likelihoods/laplace.py
+++ b/GPy/likelihoods/laplace.py
@@ -250,8 +250,11 @@ class Laplace(likelihood):
         self.W = -self.noise_model.d2logpdf_df2(self.f_hat, self.data, extra_data=self.extra_data)
 
         if not self.noise_model.log_concave:
-            #print "Under 1e-10: {}".format(np.sum(self.W < 1e-6))
-            self.W[self.W < 1e-6] = 1e-6  # FIXME-HACK: This is a hack since GPy can't handle negative variances which can occur
+            i = self.W < 1e-6
+            if np.any(i):
+                warnings.warn('truncating non log-concave likelihood curvature')
+                # FIXME-HACK: This is a hack since GPy can't handle negative variances which can occur
+                self.W[i] = 1e-6
 
         self.W12BiW12, self.ln_B_det = self._compute_B_statistics(self.K, self.W, np.eye(self.N))
 
@@ -270,14 +273,14 @@ class Laplace(likelihood):
         :type W: Vector of diagonal values of hessian (1xN)
         :param a: Matrix to calculate W12BiW12a
         :type a: Matrix NxN
-        :returns: (W12BiW12, ln_B_det)
+        :returns: (W12BiW12a, ln_B_det)
         """
         if not self.noise_model.log_concave:
             #print "Under 1e-10: {}".format(np.sum(W < 1e-6))
-            W[W < 1e-6] = 1e-6  # FIXME-HACK: This is a hack since GPy can't handle negative variances which can occur
-                                # If the likelihood is non-log-concave. We wan't to say that there is a negative variance
-                                # To cause the posterior to become less certain than the prior and likelihood,
-                                # This is a property only held by non-log-concave likelihoods
+            W[W < 1e-10] = 1e-10  # FIXME-HACK: This is a hack since GPy can't handle negative variances which can occur
+                                  # If the likelihood is non-log-concave. We wan't to say that there is a negative variance
+                                  # To cause the posterior to become less certain than the prior and likelihood,
+                                  # This is a property only held by non-log-concave likelihoods
 
 
         #W is diagonal so its sqrt is just the sqrt of the diagonal elements
diff --git a/GPy/likelihoods/noise_models/noise_distributions.py b/GPy/likelihoods/noise_models/noise_distributions.py
index a67d8792..b65b7750 100644
--- a/GPy/likelihoods/noise_models/noise_distributions.py
+++ b/GPy/likelihoods/noise_models/noise_distributions.py
@@ -153,9 +153,11 @@ class NoiseDistribution(object):
         :param sigma: standard deviation of posterior
 
         """
+        #import ipdb; ipdb.set_trace()
         def int_mean(f,m,v):
             return self._mean(f)*np.exp(-(0.5/v)*np.square(f - m))
-        scaled_mean = [quad(int_mean, -np.inf, np.inf,args=(mj,s2j))[0] for mj,s2j in zip(mu,variance)]
+        #scaled_mean = [quad(int_mean, -np.inf, np.inf,args=(mj,s2j))[0] for mj,s2j in zip(mu,variance)]
+        scaled_mean = [quad(int_mean, mj-6*np.sqrt(s2j), mj+6*np.sqrt(s2j), args=(mj,s2j))[0] for mj,s2j in zip(mu,variance)]
         mean = np.array(scaled_mean)[:,None] / np.sqrt(2*np.pi*(variance))
 
         return mean
@@ -172,16 +174,16 @@ class NoiseDistribution(object):
         :predictive_mean: output's predictive mean, if None _predictive_mean function will be called.
 
         """
-        #sigma2 = sigma**2
         normalizer = np.sqrt(2*np.pi*variance)
 
         # E( V(Y_star|f_star) )
         def int_var(f,m,v):
             return self._variance(f)*np.exp(-(0.5/v)*np.square(f - m))
-        scaled_exp_variance = [quad(int_var, -np.inf, np.inf,args=(mj,s2j))[0] for mj,s2j in zip(mu,variance)]
+        #Most of the weight is within 6 stds and this avoids some negative infinity and infinity problems of taking f^2
+        scaled_exp_variance = [quad(int_var, mj-6*np.sqrt(s2j), mj+6*np.sqrt(s2j), args=(mj,s2j))[0] for mj,s2j in zip(mu,variance)]
         exp_var = np.array(scaled_exp_variance)[:,None] / normalizer
 
-        #V( E(Y_star|f_star) ) =  E( E(Y_star|f_star)**2 ) - E( E(Y_star|f_star) )**2
+        #V( E(Y_star|f_star) ) = E( E(Y_star|f_star)**2 ) - E( E(Y_star|f_star) )**2
 
         #E( E(Y_star|f_star) )**2
         if predictive_mean is None:
@@ -189,9 +191,9 @@ class NoiseDistribution(object):
         predictive_mean_sq = predictive_mean**2
 
         #E( E(Y_star|f_star)**2 )
-        def int_pred_mean_sq(f,m,v,predictive_mean_sq):
+        def int_pred_mean_sq(f,m,v):
             return self._mean(f)**2*np.exp(-(0.5/v)*np.square(f - m))
-        scaled_exp_exp2 = [quad(int_pred_mean_sq, -np.inf, np.inf,args=(mj,s2j,pm2j))[0] for mj,s2j,pm2j in zip(mu,variance,predictive_mean_sq)]
+        scaled_exp_exp2 = [quad(int_pred_mean_sq, mj-6*np.sqrt(s2j), mj+6*np.sqrt(s2j), args=(mj,s2j))[0] for mj,s2j in zip(mu,variance)]
         exp_exp2 = np.array(scaled_exp_exp2)[:,None] / normalizer
 
         var_exp = exp_exp2 - predictive_mean_sq
@@ -408,17 +410,16 @@ class NoiseDistribution(object):
             axis=-1
 
             #Calculate mean, variance and precentiles from samples
-            print "WARNING: Using sampling to calculate mean, variance and predictive quantiles."
+            warnings.warn("Using sampling to calculate mean, variance and predictive quantiles.")
             pred_mean = np.mean(samples, axis=axis)[:,None]
             pred_var = np.var(samples, axis=axis)[:,None]
             q1 = np.percentile(samples, 2.5, axis=axis)[:,None]
             q3 = np.percentile(samples, 97.5, axis=axis)[:,None]
 
         else:
-
             pred_mean = self.predictive_mean(mu, var)
             pred_var = self.predictive_variance(mu, var, pred_mean)
-            print "WARNING: Predictive quantiles are only computed when sampling."
+            warnings.warn("Predictive quantiles are only computed when sampling.")
             q1 = np.repeat(np.nan,pred_mean.size)[:,None]
             q3 = q1.copy()
 
diff --git a/GPy/models.py b/GPy/models.py
index 3b2683ea..0aea59a0 100644
--- a/GPy/models.py
+++ b/GPy/models.py
@@ -1,18 +1,20 @@
 '''
-GPy Models
-==========
+.. module:: GPy.models
 
 Implementations for common models used in GP regression and classification.
 The different models can be viewed in :mod:`GPy.models_modules`, which holds
 detailed explanations for the different models.
 
-:warning: This module is a convienince module for endusers to use. For developers 
-see :mod:`GPy.models_modules`, which holds the implementions for each model. 
+.. note::
+    This module is a convienince module for endusers to use. For developers 
+    see :mod:`GPy.models_modules`, which holds the implementions for each model.: 
+
+.. moduleauthor:: Max Zwiessele <ibinbei@gmail.com>
 '''
 
 __updated__ = '2013-11-28'
 
-from models_modules.bayesian_gplvm import BayesianGPLVM
+from models_modules.bayesian_gplvm import BayesianGPLVM, BayesianGPLVMWithMissingData
 from models_modules.gp_regression import GPRegression
 from models_modules.gp_classification import GPClassification#; _gp_classification = gp_classification ; del gp_classification 
 from models_modules.sparse_gp_regression import SparseGPRegression#; _sparse_gp_regression = sparse_gp_regression ; del sparse_gp_regression 
diff --git a/GPy/models_modules/bayesian_gplvm.py b/GPy/models_modules/bayesian_gplvm.py
index 90e54111..57e50955 100644
--- a/GPy/models_modules/bayesian_gplvm.py
+++ b/GPy/models_modules/bayesian_gplvm.py
@@ -2,17 +2,18 @@
 # Licensed under the BSD 3-clause license (see LICENSE.txt)
 
 import numpy as np
+import itertools
+from matplotlib import pyplot
+
 from ..core.sparse_gp import SparseGP
 from ..likelihoods import Gaussian
 from .. import kern
-import itertools
-from matplotlib.colors import colorConverter
-from GPy.inference.optimization import SCG
-from GPy.util import plot_latent, linalg
-from .gplvm import GPLVM
-from GPy.util.plot_latent import most_significant_input_dimensions
-from matplotlib import pyplot
-from GPy.core.model import Model
+from ..inference.optimization import SCG
+from ..util import plot_latent, linalg
+from .gplvm import GPLVM, initialise_latent
+from ..util.plot_latent import most_significant_input_dimensions
+from ..core.model import Model
+from ..util.subarray_and_sorting import common_subarrays
 
 class BayesianGPLVM(SparseGP, GPLVM):
     """
@@ -34,7 +35,7 @@ class BayesianGPLVM(SparseGP, GPLVM):
             likelihood = likelihood_or_Y
 
         if X == None:
-            X = self.initialise_latent(init, input_dim, likelihood.Y)
+            X = initialise_latent(init, input_dim, likelihood.Y)
         self.init = init
 
         if X_variance is None:
@@ -308,14 +309,36 @@ class BayesianGPLVMWithMissingData(Model):
     :type init: 'PCA' | 'random'
     """
     def __init__(self, likelihood_or_Y, input_dim, X=None, X_variance=None, init='PCA', num_inducing=10,
-                 Z=None, kernel=None, missing=np.nan, **kwargs):
+                 Z=None, kernel=None, **kwargs):
+        #=======================================================================
+        # Filter Y, such that same missing data is at same positions. 
+        # If full rows are missing, delete them entirely!
         if type(likelihood_or_Y) is np.ndarray:
-            likelihood = Gaussian(likelihood_or_Y)
+            Y = likelihood_or_Y
+            likelihood = Gaussian
+            params = 1.
+            normalize=None
         else:
-            likelihood = likelihood_or_Y
+            Y = likelihood_or_Y.Y
+            likelihood = likelihood_or_Y.__class__
+            params = likelihood_or_Y._get_params()
+            if isinstance(likelihood_or_Y, Gaussian):
+                normalize = True
+                scale = likelihood_or_Y._scale
+                offset = likelihood_or_Y._offset
+        # Get common subrows
+        filter_ = np.isnan(Y)
+        self.subarray_indices = common_subarrays(filter_,axis=1)
+        likelihoods = [likelihood(Y[~np.array(v,dtype=bool),:][:,ind]) for v,ind in self.subarray_indices.iteritems()]
+        for l in likelihoods:
+            l._set_params(params)
+            if normalize: # get normalization in common
+                l._scale = scale
+                l._offset = offset
+        #=======================================================================
 
         if X == None:
-            X = self.initialise_latent(init, input_dim, likelihood.Y)
+            X = initialise_latent(init, input_dim, Y[:,np.any(np.isnan(Y),1)])
         self.init = init
 
         if X_variance is None:
@@ -328,13 +351,52 @@ class BayesianGPLVMWithMissingData(Model):
         if kernel is None:
             kernel = kern.rbf(input_dim) # + kern.white(input_dim)
 
-        SparseGP.__init__(self, X, likelihood, kernel, Z=Z, X_variance=X_variance, **kwargs)
+        self.submodels = [BayesianGPLVM(l, input_dim, X, X_variance, init, num_inducing, Z, kernel) for l in likelihoods]
+        self.gref = self.submodels[0] 
+        #:type self.gref: BayesianGPLVM 
         self.ensure_default_constraints()
 
+    def log_likelihood(self):
+        ll = -self.gref.KL_divergence()
+        for g in self.submodels:
+            ll += SparseGP.log_likelihood(g)
+        return ll
+
+    def _log_likelihood_gradients(self):
+        dLdmu, dLdS = reduce(lambda a, b: [a[0] + b[0], a[1] + b[1]], (g.dL_dmuS() for g in self.bgplvms))
+        dKLmu, dKLdS = self.gref.dKL_dmuS()
+        dLdmu -= dKLmu
+        dLdS -= dKLdS
+        dLdmuS = np.hstack((dLdmu.flatten(), dLdS.flatten())).flatten()
+        dldzt1 = reduce(lambda a, b: a + b, (SparseGP._log_likelihood_gradients(g)[:self.gref.num_inducing*self.gref.input_dim] for g in self.submodels))
+
+        return np.hstack((dLdmuS,
+                             dldzt1,
+                np.hstack([np.hstack([g.dL_dtheta(),
+                                            g.likelihood._gradients(\
+                                                partial=g.partial_for_likelihood)]) \
+                              for g in self.submodels])))
+
+    def getstate(self):
+        return Model.getstate(self)+[self.submodels,self.subarray_indices]
+
+    def setstate(self, state):
+        self.subarray_indices = state.pop()
+        self.submodels = state.pop()
+        self.gref = self.submodels[0]
+        Model.setstate(self, state)
+        self._set_params(self._get_params())
+
     def _get_param_names(self):
         X_names = sum([['X_%i_%i' % (n, q) for q in range(self.input_dim)] for n in range(self.num_data)], [])
         S_names = sum([['X_variance_%i_%i' % (n, q) for q in range(self.input_dim)] for n in range(self.num_data)], [])
-        return (X_names + S_names + SparseGP._get_param_names(self))
+        return (X_names + S_names + SparseGP._get_param_names(self.gref))
+
+    def _get_params(self):
+        return self.gref._get_params()
+    def _set_params(self, x):
+        [g._set_params(x) for g in self.submodels]
+    
 
     pass
 
diff --git a/GPy/models_modules/gplvm.py b/GPy/models_modules/gplvm.py
index f27f861c..a5fe4284 100644
--- a/GPy/models_modules/gplvm.py
+++ b/GPy/models_modules/gplvm.py
@@ -9,7 +9,14 @@ from ..core import priors
 from ..core import GP
 from ..likelihoods import Gaussian
 from .. import util
+from ..util.linalg import pca
 
+def initialise_latent(init, input_dim, Y):
+    Xr = np.random.randn(Y.shape[0], input_dim)
+    if init.lower() == 'pca':
+        PC = pca(Y, input_dim)[0]
+        Xr[:PC.shape[0], :PC.shape[1]] = PC
+    return Xr
 
 class GPLVM(GP):
     """
@@ -20,12 +27,12 @@ class GPLVM(GP):
     :param input_dim: latent dimensionality
     :type input_dim: int
     :param init: initialisation method for the latent space
-    :type init: 'PCA'|'random'
+    :type init: 'pca'|'random'
 
     """
     def __init__(self, Y, input_dim, init='PCA', X=None, kernel=None, normalize_Y=False):
         if X is None:
-            X = self.initialise_latent(init, input_dim, Y)
+            X = initialise_latent(init, input_dim, Y)
         if kernel is None:
             kernel = kern.rbf(input_dim, ARD=input_dim > 1) + kern.bias(input_dim, np.exp(-2))
         likelihood = Gaussian(Y, normalize=normalize_Y, variance=np.exp(-2.))
@@ -33,14 +40,6 @@ class GPLVM(GP):
         self.set_prior('.*X', priors.Gaussian(0, 1))
         self.ensure_default_constraints()
 
-    def initialise_latent(self, init, input_dim, Y):
-        Xr = np.random.randn(Y.shape[0], input_dim)
-        if init == 'PCA':
-            from ..util.linalg import PCA
-            PC = PCA(Y, input_dim)[0]
-            Xr[:PC.shape[0], :PC.shape[1]] = PC
-        return Xr
-
     def _get_param_names(self):
         return sum([['X_%i_%i' % (n, q) for q in range(self.input_dim)] for n in range(self.num_data)], []) + GP._get_param_names(self)
 
@@ -61,7 +60,7 @@ class GPLVM(GP):
         for i in range(self.output_dim):
             target[:,:,i] = self.kern.dK_dX(np.dot(self.Ki,self.likelihood.Y[:,i])[None, :],X,self.X)
         return target
-   
+
     def magnification(self,X):
         target=np.zeros(X.shape[0])
         J = np.zeros((X.shape[0],X.shape[1],self.output_dim))
diff --git a/GPy/models_modules/mrd.py b/GPy/models_modules/mrd.py
index b9c99a64..862b697b 100644
--- a/GPy/models_modules/mrd.py
+++ b/GPy/models_modules/mrd.py
@@ -5,7 +5,7 @@ Created on 10 Apr 2013
 '''
 from GPy.core import Model
 from GPy.core import SparseGP
-from GPy.util.linalg import PCA
+from GPy.util.linalg import pca
 import numpy
 import itertools
 import pylab
@@ -26,8 +26,8 @@ class MRD(Model):
     :type input_dim: int
     :param initx: initialisation method for the latent space :
 
-        * 'concat' - PCA on concatenation of all datasets
-        * 'single' - Concatenation of PCA on datasets, respectively
+        * 'concat' - pca on concatenation of all datasets
+        * 'single' - Concatenation of pca on datasets, respectively
         * 'random' - Random draw from a normal
 
     :type initx: ['concat'|'single'|'random']
@@ -248,11 +248,11 @@ class MRD(Model):
                 Ylist.append(likelihood_or_Y.Y)
         del likelihood_list
         if init in "PCA_concat":
-            X = PCA(numpy.hstack(Ylist), self.input_dim)[0]
+            X = pca(numpy.hstack(Ylist), self.input_dim)[0]
         elif init in "PCA_single":
             X = numpy.zeros((Ylist[0].shape[0], self.input_dim))
             for qs, Y in itertools.izip(numpy.array_split(numpy.arange(self.input_dim), len(Ylist)), Ylist):
-                X[:, qs] = PCA(Y, len(qs))[0]
+                X[:, qs] = pca(Y, len(qs))[0]
         else: # init == 'random':
             X = numpy.random.randn(Ylist[0].shape[0], self.input_dim)
         self.X = X
diff --git a/GPy/models_modules/sparse_gplvm.py b/GPy/models_modules/sparse_gplvm.py
index 4e401ee3..04d3415a 100644
--- a/GPy/models_modules/sparse_gplvm.py
+++ b/GPy/models_modules/sparse_gplvm.py
@@ -6,10 +6,7 @@ import numpy as np
 import pylab as pb
 import sys, pdb
 from sparse_gp_regression import SparseGPRegression
-from gplvm import GPLVM
-# from .. import kern
-# from ..core import model
-# from ..util.linalg import pdinv, PCA
+from gplvm import GPLVM, initialise_latent
 
 class SparseGPLVM(SparseGPRegression, GPLVM):
     """
@@ -24,7 +21,7 @@ class SparseGPLVM(SparseGPRegression, GPLVM):
 
     """
     def __init__(self, Y, input_dim, kernel=None, init='PCA', num_inducing=10):
-        X = self.initialise_latent(init, input_dim, Y)
+        X = initialise_latent(init, input_dim, Y)
         SparseGPRegression.__init__(self, X, Y, kernel=kernel, num_inducing=num_inducing)
         self.ensure_default_constraints()
 
diff --git a/GPy/testing/examples_tests.py b/GPy/testing/examples_tests.py
index 9998590a..be26fff6 100644
--- a/GPy/testing/examples_tests.py
+++ b/GPy/testing/examples_tests.py
@@ -19,25 +19,12 @@ class ExamplesTests(unittest.TestCase):
     def _model_instance(self, Model):
         self.assertTrue(isinstance(Model, GPy.models))
 
-"""
-def model_instance_generator(model):
-    def check_model_returned(self):
-        self._model_instance(model)
-    return check_model_returned
-
-def checkgrads_generator(model):
-    def model_checkgrads(self):
-        self._checkgrad(model)
-    return model_checkgrads
-"""
-
 def model_checkgrads(model):
     model.randomize()
-    #assert model.checkgrad()
-    return model.checkgrad()
+    #NOTE: Step as 1e-4, this should be acceptable for more peaky models
+    return model.checkgrad(step=1e-4)
 
 def model_instance(model):
-    #assert isinstance(model, GPy.core.model)
     return isinstance(model, GPy.core.model.Model)
 
 def flatten_nested(lst):
@@ -49,7 +36,7 @@ def flatten_nested(lst):
             result.append(element)
     return result
 
-#@nottest
+@nottest
 def test_models():
     optimize=False
     plot=True
@@ -66,9 +53,11 @@ def test_models():
         print "After"
         print functions
         for example in functions:
-            #if example[0] in ['oil', 'silhouette', 'GPLVM_oil_100', 'brendan_faces']:
-                #print "SKIPPING"
-                #continue
+            if example[0] in ['epomeo_gpx']:
+                #These are the edge cases that we might want to handle specially
+                if example[0] == 'epomeo_gpx' and not GPy.util.datasets.gpxpy_available:
+                    print "Skipping as gpxpy is not available to parse GPS"
+                    continue
 
             print "Testing example: ", example[0]
             # Generate model
diff --git a/GPy/testing/likelihoods_tests.py b/GPy/testing/likelihood_tests.py
similarity index 99%
rename from GPy/testing/likelihoods_tests.py
rename to GPy/testing/likelihood_tests.py
index 58c9a64b..d14c9a41 100644
--- a/GPy/testing/likelihoods_tests.py
+++ b/GPy/testing/likelihood_tests.py
@@ -593,7 +593,6 @@ class LaplaceTests(unittest.TestCase):
         grad.checkgrad(verbose=1)
         self.assertTrue(grad.checkgrad())
 
-    #@unittest.skip('Not working yet, needs to be checked')
     def test_laplace_log_likelihood(self):
         debug = False
         real_std = 0.1
diff --git a/GPy/util/data_resources.json b/GPy/util/data_resources.json
index d86d9088..c999b796 100644
--- a/GPy/util/data_resources.json
+++ b/GPy/util/data_resources.json
@@ -29,7 +29,8 @@
       "urls":[
          "http://staffwww.dcs.shef.ac.uk/people/N.Lawrence/dataset_mirror/ankur_pose_data/"
       ],
-      "details":"Artificially generated data of silhouettes given poses. Note that the data does not display a left/right ambiguity because across the entire data set one of the arms sticks out more the the other, disambiguating the pose as to which way the individual is facing."
+      "details":"Artificially generated data of silhouettes given poses. Note that the data does not display a left/right ambiguity because across the entire data set one of the arms sticks out more the the other, disambiguating the pose as to which way the individual is facing.",
+      "size":1
    },
    "osu_accad":{
       "files":[
@@ -316,4 +317,4 @@
       ],
       "size":2031872
    }
-}
\ No newline at end of file
+}
diff --git a/GPy/util/datasets.py b/GPy/util/datasets.py
index 7fd1b6c5..83983832 100644
--- a/GPy/util/datasets.py
+++ b/GPy/util/datasets.py
@@ -26,13 +26,16 @@ def reporthook(a,b,c):
 # Global variables
 data_path = os.path.join(os.path.dirname(__file__), 'datasets')
 default_seed = 10000
-overide_manual_authorize=False
+overide_manual_authorize=True
 neil_url = 'http://staffwww.dcs.shef.ac.uk/people/N.Lawrence/dataset_mirror/'
 
 # Read data resources from json file.
-path = os.path.join(os.path.dirname(__file__), 'data_resources.json')
-json_data=open(path).read()
-data_resources = json.loads(json_data)
+# Don't do this when ReadTheDocs is scanning as it breaks things
+on_rtd = os.environ.get('READTHEDOCS', None) == 'True' #Checks if RTD is scanning
+if not (on_rtd):
+    path = os.path.join(os.path.dirname(__file__), 'data_resources.json')
+    json_data=open(path).read()
+    data_resources = json.loads(json_data)
 
 
 def prompt_user(prompt):
@@ -94,7 +97,7 @@ def download_url(url, store_directory, save_name = None, messages = True, suffix
     # if we wanted to get more sophisticated maybe we should check the response code here again even for successes.
     with open(save_name, 'wb') as f:
         f.write(response.read())
-    
+
     #urllib.urlretrieve(url+suffix, save_name, reporthook)
 
 def authorize_download(dataset_name=None):
@@ -232,7 +235,7 @@ if gpxpy_available:
             gpx_file.close()
         return data_details_return({'X' : X, 'info' : 'Data is an array containing time in seconds, latitude, longitude and elevation in that order.'}, data_set)
 
-del gpxpy_available
+#del gpxpy_available
 
 
 
diff --git a/GPy/util/diag.py b/GPy/util/diag.py
new file mode 100644
index 00000000..3d6b4dc9
--- /dev/null
+++ b/GPy/util/diag.py
@@ -0,0 +1,114 @@
+'''
+.. module:: GPy.util.diag
+
+.. moduleauthor:: Max Zwiessele <ibinbei@gmail.com>
+
+'''
+__updated__ = '2013-12-03'
+
+import numpy as np
+
+def view(A, offset=0):
+    """
+    Get a view on the diagonal elements of a 2D array.
+    
+    This is actually a view (!) on the diagonal of the array, so you can 
+    in-place adjust the view.
+    
+    :param :class:`ndarray` A: 2 dimensional numpy array
+    :param int offset: view offset to give back (negative entries allowed)
+    :rtype: :class:`ndarray` view of diag(A)
+    
+    >>> import numpy as np
+    >>> X = np.arange(9).reshape(3,3)
+    >>> view(X)
+    array([0, 4, 8])
+    >>> d = view(X)
+    >>> d += 2
+    >>> view(X)
+    array([ 2,  6, 10])
+    >>> view(X, offset=-1)
+    array([3, 7])
+    >>> subtract(X, 3, offset=-1)
+    array([[ 2,  1,  2],
+           [ 0,  6,  5],
+           [ 6,  4, 10]])
+    """
+    from numpy.lib.stride_tricks import as_strided
+    assert A.ndim == 2, "only implemented for 2 dimensions"
+    assert A.shape[0] == A.shape[1], "attempting to get the view of non-square matrix?!" 
+    if offset > 0:
+        return as_strided(A[0, offset:], shape=(A.shape[0] - offset, ), strides=((A.shape[0]+1)*A.itemsize, ))
+    elif offset < 0:
+        return as_strided(A[-offset:, 0], shape=(A.shape[0] + offset, ), strides=((A.shape[0]+1)*A.itemsize, ))
+    else:
+        return as_strided(A, shape=(A.shape[0], ), strides=((A.shape[0]+1)*A.itemsize, ))
+
+def _diag_ufunc(A,b,offset,func):
+    dA = view(A, offset); func(dA,b,dA)
+    return A
+
+def times(A, b, offset=0):
+    """
+    Times the view of A with b in place (!).
+    Returns modified A 
+    Broadcasting is allowed, thus b can be scalar.
+    
+    if offset is not zero, make sure b is of right shape!
+    
+    :param ndarray A: 2 dimensional array
+    :param ndarray-like b: either one dimensional or scalar
+    :param int offset: same as in view.
+    :rtype: view of A, which is adjusted inplace
+    """
+    return _diag_ufunc(A, b, offset, np.multiply)
+multiply = times
+
+def divide(A, b, offset=0):
+    """
+    Divide the view of A by b in place (!).
+    Returns modified A 
+    Broadcasting is allowed, thus b can be scalar.
+    
+    if offset is not zero, make sure b is of right shape!
+    
+    :param ndarray A: 2 dimensional array
+    :param ndarray-like b: either one dimensional or scalar
+    :param int offset: same as in view.
+    :rtype: view of A, which is adjusted inplace
+    """
+    return _diag_ufunc(A, b, offset, np.divide)
+
+def add(A, b, offset=0):
+    """
+    Add b to the view of A in place (!).
+    Returns modified A.
+    Broadcasting is allowed, thus b can be scalar.
+    
+    if offset is not zero, make sure b is of right shape!
+    
+    :param ndarray A: 2 dimensional array
+    :param ndarray-like b: either one dimensional or scalar
+    :param int offset: same as in view.
+    :rtype: view of A, which is adjusted inplace
+    """
+    return _diag_ufunc(A, b, offset, np.add)
+
+def subtract(A, b, offset=0):
+    """
+    Subtract b from the view of A in place (!).
+    Returns modified A.
+    Broadcasting is allowed, thus b can be scalar.
+    
+    if offset is not zero, make sure b is of right shape!
+    
+    :param ndarray A: 2 dimensional array
+    :param ndarray-like b: either one dimensional or scalar
+    :param int offset: same as in view.
+    :rtype: view of A, which is adjusted inplace
+    """
+    return _diag_ufunc(A, b, offset, np.subtract)
+        
+if __name__ == '__main__':
+    import doctest
+    doctest.testmod()
\ No newline at end of file
diff --git a/GPy/util/linalg.py b/GPy/util/linalg.py
index e3e421f6..842178e2 100644
--- a/GPy/util/linalg.py
+++ b/GPy/util/linalg.py
@@ -217,7 +217,7 @@ def multiple_pdinv(A):
     return np.dstack(invs), np.array(halflogdets)
 
 
-def PCA(Y, input_dim):
+def pca(Y, input_dim):
     """
     Principal component analysis: maximum likelihood solution by SVD
 
@@ -230,7 +230,7 @@ def PCA(Y, input_dim):
 
     """
     if not np.allclose(Y.mean(axis=0), 0.0):
-        print "Y is not zero mean, centering it locally (GPy.util.linalg.PCA)"
+        print "Y is not zero mean, centering it locally (GPy.util.linalg.pca)"
 
         # Y -= Y.mean(axis=0)
 
@@ -241,6 +241,124 @@ def PCA(Y, input_dim):
     W *= v;
     return X, W.T
 
+def ppca(Y, Q, iterations=100):
+    """
+    EM implementation for probabilistic pca.
+
+    :param array-like Y: Observed Data
+    :param int Q: Dimensionality for reduced array
+    :param int iterations: number of iterations for EM
+    """
+    from numpy.ma import dot as madot
+    N, D = Y.shape
+    # Initialise W randomly
+    W = np.random.randn(D, Q) * 1e-3
+    Y = np.ma.masked_invalid(Y, copy=0)
+    mu = Y.mean(0)
+    Ycentered = Y - mu
+    try:
+        for _ in range(iterations):
+            exp_x = np.asarray_chkfinite(np.linalg.solve(W.T.dot(W), madot(W.T, Ycentered.T))).T
+            W = np.asarray_chkfinite(np.linalg.solve(exp_x.T.dot(exp_x), madot(exp_x.T, Ycentered))).T
+    except np.linalg.linalg.LinAlgError:
+        #"converged"
+        pass
+    return np.asarray_chkfinite(exp_x), np.asarray_chkfinite(W)
+
+def ppca_missing_data_at_random(Y, Q, iters=100):
+    """
+    EM implementation of Probabilistic pca for when there is missing data.
+    
+    Taken from <SheffieldML, https://github.com/SheffieldML>
+
+    .. math:
+        \\mathbf{Y} = \mathbf{XW} + \\epsilon \\text{, where}
+        \\epsilon = \\mathcal{N}(0, \\sigma^2 \mathbf{I})
+        
+    :returns: X, W, sigma^2 
+    """
+    from numpy.ma import dot as madot
+    import diag
+    from GPy.util.subarray_and_sorting import common_subarrays
+    import time
+    debug = 1
+    # Initialise W randomly
+    N, D = Y.shape
+    W = np.random.randn(Q, D) * 1e-3
+    Y = np.ma.masked_invalid(Y, copy=1)
+    nu = 1.
+    #num_obs_i = 1./Y.count()
+    Ycentered = Y - Y.mean(0)
+    
+    X = np.zeros((N,Q))
+    cs = common_subarrays(Y.mask)
+    cr = common_subarrays(Y.mask, 1)
+    Sigma = np.zeros((N, Q, Q))
+    Sigma2 = np.zeros((N, Q, Q))
+    mu = np.zeros(D)
+    if debug:
+        import matplotlib.pyplot as pylab
+        fig = pylab.figure("FIT MISSING DATA"); 
+        ax = fig.gca()
+        ax.cla()
+        lines = pylab.plot(np.zeros((N,Q)).dot(W))
+    W2 = np.zeros((Q,D))
+
+    for i in range(iters):
+#         Sigma = np.linalg.solve(diag.add(madot(W,W.T), nu), diag.times(np.eye(Q),nu))
+#         exp_x = madot(madot(Ycentered, W.T),Sigma)/nu
+#         Ycentered = (Y - exp_x.dot(W).mean(0))
+#         #import ipdb;ipdb.set_trace()
+#         #Ycentered = mu
+#         W = np.linalg.solve(madot(exp_x.T,exp_x) + Sigma, madot(exp_x.T, Ycentered))
+#         nu = (((Ycentered - madot(exp_x, W))**2).sum(0) + madot(W.T,madot(Sigma,W)).sum(0)).sum()/N
+        for csi, (mask, index) in enumerate(cs.iteritems()):
+            mask = ~np.array(mask)
+            Sigma2[index, :, :] = nu * np.linalg.inv(diag.add(W2[:,mask].dot(W2[:,mask].T), nu))
+            #X[index,:] = madot((Sigma[csi]/nu),madot(W,Ycentered[index].T))[:,0]
+        X2 = ((Sigma2/nu) * (madot(Ycentered,W2.T).base)[:,:,None]).sum(-1)
+        mu2 = (Y - X.dot(W)).mean(0)
+        for n in range(N):
+            Sigma[n] = nu * np.linalg.inv(diag.add(W[:,~Y.mask[n]].dot(W[:,~Y.mask[n]].T), nu))
+            X[n, :] = (Sigma[n]/nu).dot(W[:,~Y.mask[n]].dot(Ycentered[n,~Y.mask[n]].T))
+        for d in range(D):
+            mu[d] = (Y[~Y.mask[:,d], d] - X[~Y.mask[:,d]].dot(W[:, d])).mean()
+        Ycentered = (Y - mu)
+        nu3 = 0.
+        for cri, (mask, index) in enumerate(cr.iteritems()):
+            mask = ~np.array(mask)
+            W2[:,index] = np.linalg.solve(X[mask].T.dot(X[mask]) + Sigma[mask].sum(0), madot(X[mask].T, Ycentered[mask,index]))[:,None]
+            W2[:,index] = np.linalg.solve(X.T.dot(X) + Sigma.sum(0), madot(X.T, Ycentered[:,index]))
+            #nu += (((Ycentered[mask,index] - X[mask].dot(W[:,index]))**2).sum(0) + W[:,index].T.dot(Sigma[mask].sum(0).dot(W[:,index])).sum(0)).sum()
+            nu3 += (((Ycentered[index] - X.dot(W[:,index]))**2).sum(0) + W[:,index].T.dot(Sigma.sum(0).dot(W[:,index])).sum(0)).sum()
+        nu3 /= N
+        nu = 0.
+        nu2 = 0.
+        W = np.zeros((Q,D))
+        for j in range(D):
+            W[:,j] = np.linalg.solve(X[~Y.mask[:,j]].T.dot(X[~Y.mask[:,j]]) + Sigma[~Y.mask[:,j]].sum(0), madot(X[~Y.mask[:,j]].T, Ycentered[~Y.mask[:,j],j]))
+            nu2f = np.tensordot(W[:,j].T, Sigma[~Y.mask[:,j],:,:], [0,1]).dot(W[:,j])
+            nu2s = W[:,j].T.dot(Sigma[~Y.mask[:,j],:,:].sum(0).dot(W[:,j]))
+            nu2 += (((Ycentered[~Y.mask[:,j],j] - X[~Y.mask[:,j],:].dot(W[:,j]))**2) + nu2f).sum()
+            for i in range(N):
+                if not Y.mask[i,j]:
+                    nu += ((Ycentered[i,j] - X[i,:].dot(W[:,j]))**2) + W[:,j].T.dot(Sigma[i,:,:].dot(W[:,j]))
+        nu /= N
+        nu2 /= N
+        nu4 = (((Ycentered - X.dot(W))**2).sum(0) + W.T.dot(Sigma.sum(0).dot(W)).sum(0)).sum()/N
+        import ipdb;ipdb.set_trace()
+        if debug:
+            #print Sigma[0]
+            print "nu:", nu, "sum(X):", X.sum()
+            pred_y = X.dot(W)
+            for x, l in zip(pred_y.T, lines):
+                l.set_ydata(x)
+            ax.autoscale_view()
+            ax.set_ylim(pred_y.min(), pred_y.max())
+            fig.canvas.draw()
+            time.sleep(.3)
+    return np.asarray_chkfinite(X), np.asarray_chkfinite(W), nu
+
 
 def tdot_numpy(mat, out=None):
     return np.dot(mat, mat.T, out)
diff --git a/GPy/util/plot_latent.py b/GPy/util/plot_latent.py
index 62442650..207a7974 100644
--- a/GPy/util/plot_latent.py
+++ b/GPy/util/plot_latent.py
@@ -20,8 +20,8 @@ def most_significant_input_dimensions(model, which_indices):
         input_1, input_2 = which_indices
     return input_1, input_2
 
-def plot_latent(model, labels=None, which_indices=None, 
-                resolution=50, ax=None, marker='o', s=40, 
+def plot_latent(model, labels=None, which_indices=None,
+                resolution=50, ax=None, marker='o', s=40,
                 fignum=None, plot_inducing=False, legend=True,
                 aspect='auto', updates=False):
     """
@@ -48,10 +48,10 @@ def plot_latent(model, labels=None, which_indices=None,
         var = var[:, :1]
         return np.log(var)
     view = ImshowController(ax, plot_function,
-                            tuple(model.X.min(0)[:, [input_1, input_2]]) + tuple(model.X.max(0)[:, [input_1, input_2]]),
+                            tuple(model.X[:, [input_1, input_2]].min(0)) + tuple(model.X[:, [input_1, input_2]].max(0)),
                             resolution, aspect=aspect, interpolation='bilinear',
                             cmap=pb.cm.binary)
-    
+
 #     ax.imshow(var.reshape(resolution, resolution).T,
 #               extent=[xmin[0], xmax[0], xmin[1], xmax[1]], cmap=pb.cm.binary, interpolation='bilinear', origin='lower')
 
@@ -100,8 +100,8 @@ def plot_latent(model, labels=None, which_indices=None,
         raw_input('Enter to continue')
     return ax
 
-def plot_magnification(model, labels=None, which_indices=None, 
-                resolution=60, ax=None, marker='o', s=40, 
+def plot_magnification(model, labels=None, which_indices=None,
+                resolution=60, ax=None, marker='o', s=40,
                 fignum=None, plot_inducing=False, legend=True,
                 aspect='auto', updates=False):
     """
diff --git a/GPy/util/subarray_and_sorting.py b/GPy/util/subarray_and_sorting.py
new file mode 100644
index 00000000..49385771
--- /dev/null
+++ b/GPy/util/subarray_and_sorting.py
@@ -0,0 +1,56 @@
+'''
+.. module:: GPy.util.subarray_and_sorting
+
+.. moduleauthor:: Max Zwiessele <ibinbei@gmail.com>
+
+'''
+__updated__ = '2013-12-02'
+
+import numpy as np
+
+def common_subarrays(X, axis=0):
+    """
+    Find common subarrays of 2 dimensional X, where axis is the axis to apply the search over.
+    Common subarrays are returned as a dictionary of <subarray, [index]> pairs, where
+    the subarray is a tuple representing the subarray and the index is the index
+    for the subarray in X, where index is the index to the remaining axis.
+    
+    :param :class:`np.ndarray` X: 2d array to check for common subarrays in
+    :param int axis: axis to apply subarray detection over. 
+        When the index is 0, compare rows, columns, otherwise.   
+    
+    Examples:
+    =========
+
+    In a 2d array:    
+    >>> import numpy as np
+    >>> X = np.zeros((3,6), dtype=bool)
+    >>> X[[1,1,1],[0,4,5]] = 1; X[1:,[2,3]] = 1
+    >>> X
+    array([[False, False, False, False, False, False],
+           [ True, False,  True,  True,  True,  True],
+           [False, False,  True,  True, False, False]], dtype=bool)
+    >>> d = common_subarrays(X,axis=1)
+    >>> len(d)
+    3
+    >>> X[:, d[tuple(X[:,0])]]
+    array([[False, False, False],
+           [ True,  True,  True],
+           [False, False, False]], dtype=bool)
+    >>> d[tuple(X[:,4])] == d[tuple(X[:,0])] == [0, 4, 5]
+    True
+    >>> d[tuple(X[:,1])]
+    [1]
+    """
+    from collections import defaultdict
+    from itertools import count
+    from operator import iadd
+    assert X.ndim == 2 and axis in (0,1), "Only implemented for 2D arrays"
+    subarrays = defaultdict(list)
+    cnt = count()
+    np.apply_along_axis(lambda x: iadd(subarrays[tuple(x)], [cnt.next()]), 1-axis, X)
+    return subarrays
+
+if __name__ == '__main__':
+    import doctest
+    doctest.testmod()
\ No newline at end of file
diff --git a/doc/GPy.core.rst b/doc/GPy.core.rst
index d7f18192..c4f1849d 100644
--- a/doc/GPy.core.rst
+++ b/doc/GPy.core.rst
@@ -1,107 +1,102 @@
-core Package
-============
+GPy.core package
+================
 
-:mod:`core` Package
--------------------
+Submodules
+----------
 
-.. automodule:: GPy.core
-    :members:
-    :undoc-members:
-    :show-inheritance:
-
-:mod:`domains` Module
----------------------
+GPy.core.domains module
+-----------------------
 
 .. automodule:: GPy.core.domains
     :members:
     :undoc-members:
     :show-inheritance:
 
-:mod:`fitc` Module
-------------------
+GPy.core.fitc module
+--------------------
 
 .. automodule:: GPy.core.fitc
     :members:
     :undoc-members:
     :show-inheritance:
 
-:mod:`gp` Module
-----------------
+GPy.core.gp module
+------------------
 
 .. automodule:: GPy.core.gp
     :members:
     :undoc-members:
     :show-inheritance:
 
-:mod:`gp_base` Module
----------------------
+GPy.core.gp_base module
+-----------------------
 
 .. automodule:: GPy.core.gp_base
     :members:
     :undoc-members:
     :show-inheritance:
 
-:mod:`mapping` Module
----------------------
+GPy.core.mapping module
+-----------------------
 
 .. automodule:: GPy.core.mapping
     :members:
     :undoc-members:
     :show-inheritance:
 
-:mod:`model` Module
--------------------
+GPy.core.model module
+---------------------
 
 .. automodule:: GPy.core.model
     :members:
     :undoc-members:
     :show-inheritance:
 
-:mod:`parameterized` Module
----------------------------
+GPy.core.parameterized module
+-----------------------------
 
 .. automodule:: GPy.core.parameterized
     :members:
     :undoc-members:
     :show-inheritance:
 
-:mod:`priors` Module
---------------------
+GPy.core.priors module
+----------------------
 
 .. automodule:: GPy.core.priors
     :members:
     :undoc-members:
     :show-inheritance:
 
-:mod:`sparse_gp` Module
------------------------
+GPy.core.sparse_gp module
+-------------------------
 
 .. automodule:: GPy.core.sparse_gp
     :members:
     :undoc-members:
     :show-inheritance:
 
-:mod:`svigp` Module
--------------------
+GPy.core.svigp module
+---------------------
 
 .. automodule:: GPy.core.svigp
     :members:
     :undoc-members:
     :show-inheritance:
 
-:mod:`transformations` Module
------------------------------
+GPy.core.transformations module
+-------------------------------
 
 .. automodule:: GPy.core.transformations
     :members:
     :undoc-members:
     :show-inheritance:
 
-:mod:`variational` Module
--------------------------
 
-.. automodule:: GPy.core.variational
+Module contents
+---------------
+
+.. automodule:: GPy.core
     :members:
     :undoc-members:
     :show-inheritance:
-
diff --git a/doc/GPy.examples.rst b/doc/GPy.examples.rst
index 176ae396..bde015dd 100644
--- a/doc/GPy.examples.rst
+++ b/doc/GPy.examples.rst
@@ -1,59 +1,62 @@
-examples Package
-================
+GPy.examples package
+====================
 
-:mod:`examples` Package
------------------------
+Submodules
+----------
 
-.. automodule:: GPy.examples
-    :members:
-    :undoc-members:
-    :show-inheritance:
-
-:mod:`classification` Module
-----------------------------
+GPy.examples.classification module
+----------------------------------
 
 .. automodule:: GPy.examples.classification
     :members:
     :undoc-members:
     :show-inheritance:
 
-:mod:`dimensionality_reduction` Module
---------------------------------------
+GPy.examples.dimensionality_reduction module
+--------------------------------------------
 
 .. automodule:: GPy.examples.dimensionality_reduction
     :members:
     :undoc-members:
     :show-inheritance:
 
-:mod:`laplace_approximations` Module
-------------------------------------
+GPy.examples.non_gaussian module
+--------------------------------
 
-.. automodule:: GPy.examples.laplace_approximations
+.. automodule:: GPy.examples.non_gaussian
     :members:
     :undoc-members:
     :show-inheritance:
 
-:mod:`regression` Module
-------------------------
+GPy.examples.regression module
+------------------------------
 
 .. automodule:: GPy.examples.regression
     :members:
     :undoc-members:
     :show-inheritance:
 
-:mod:`stochastic` Module
-------------------------
+GPy.examples.stochastic module
+------------------------------
 
 .. automodule:: GPy.examples.stochastic
     :members:
     :undoc-members:
     :show-inheritance:
 
-:mod:`tutorials` Module
------------------------
+GPy.examples.tutorials module
+-----------------------------
 
 .. automodule:: GPy.examples.tutorials
     :members:
     :undoc-members:
     :show-inheritance:
 
+
+Module contents
+---------------
+
+.. automodule:: GPy.examples
+    :members:
+    :undoc-members:
+    :show-inheritance:
diff --git a/doc/GPy.inference.rst b/doc/GPy.inference.rst
index 6a1bef4a..28f42994 100644
--- a/doc/GPy.inference.rst
+++ b/doc/GPy.inference.rst
@@ -1,51 +1,62 @@
-inference Package
-=================
+GPy.inference package
+=====================
 
-:mod:`conjugate_gradient_descent` Module
-----------------------------------------
+Submodules
+----------
+
+GPy.inference.conjugate_gradient_descent module
+-----------------------------------------------
 
 .. automodule:: GPy.inference.conjugate_gradient_descent
     :members:
     :undoc-members:
     :show-inheritance:
 
-:mod:`gradient_descent_update_rules` Module
--------------------------------------------
+GPy.inference.gradient_descent_update_rules module
+--------------------------------------------------
 
 .. automodule:: GPy.inference.gradient_descent_update_rules
     :members:
     :undoc-members:
     :show-inheritance:
 
-:mod:`optimization` Module
---------------------------
+GPy.inference.optimization module
+---------------------------------
 
 .. automodule:: GPy.inference.optimization
     :members:
     :undoc-members:
     :show-inheritance:
 
-:mod:`samplers` Module
-----------------------
+GPy.inference.samplers module
+-----------------------------
 
 .. automodule:: GPy.inference.samplers
     :members:
     :undoc-members:
     :show-inheritance:
 
-:mod:`scg` Module
------------------
+GPy.inference.scg module
+------------------------
 
 .. automodule:: GPy.inference.scg
     :members:
     :undoc-members:
     :show-inheritance:
 
-:mod:`sgd` Module
------------------
+GPy.inference.sgd module
+------------------------
 
 .. automodule:: GPy.inference.sgd
     :members:
     :undoc-members:
     :show-inheritance:
 
+
+Module contents
+---------------
+
+.. automodule:: GPy.inference
+    :members:
+    :undoc-members:
+    :show-inheritance:
diff --git a/doc/GPy.kern.parts.rst b/doc/GPy.kern.parts.rst
index 45d3e235..59c48d96 100644
--- a/doc/GPy.kern.parts.rst
+++ b/doc/GPy.kern.parts.rst
@@ -1,275 +1,278 @@
-parts Package
-=============
+GPy.kern.parts package
+======================
 
-:mod:`parts` Package
---------------------
+Submodules
+----------
 
-.. automodule:: GPy.kern.parts
-    :members:
-    :undoc-members:
-    :show-inheritance:
-
-:mod:`Brownian` Module
-----------------------
+GPy.kern.parts.Brownian module
+------------------------------
 
 .. automodule:: GPy.kern.parts.Brownian
     :members:
     :undoc-members:
     :show-inheritance:
 
-:mod:`Matern32` Module
-----------------------
+GPy.kern.parts.Matern32 module
+------------------------------
 
 .. automodule:: GPy.kern.parts.Matern32
     :members:
     :undoc-members:
     :show-inheritance:
 
-:mod:`Matern52` Module
-----------------------
+GPy.kern.parts.Matern52 module
+------------------------------
 
 .. automodule:: GPy.kern.parts.Matern52
     :members:
     :undoc-members:
     :show-inheritance:
 
-:mod:`ODE_1` Module
--------------------
+GPy.kern.parts.ODE_1 module
+---------------------------
 
 .. automodule:: GPy.kern.parts.ODE_1
     :members:
     :undoc-members:
     :show-inheritance:
 
-:mod:`ODE_UY` Module
---------------------
+GPy.kern.parts.ODE_UY module
+----------------------------
 
 .. automodule:: GPy.kern.parts.ODE_UY
     :members:
     :undoc-members:
     :show-inheritance:
 
-:mod:`bias` Module
-------------------
+GPy.kern.parts.bias module
+--------------------------
 
 .. automodule:: GPy.kern.parts.bias
     :members:
     :undoc-members:
     :show-inheritance:
 
-:mod:`coregionalize` Module
----------------------------
+GPy.kern.parts.coregionalize module
+-----------------------------------
 
 .. automodule:: GPy.kern.parts.coregionalize
     :members:
     :undoc-members:
     :show-inheritance:
 
-:mod:`eq_ode1` Module
----------------------
+GPy.kern.parts.eq_ode1 module
+-----------------------------
 
 .. automodule:: GPy.kern.parts.eq_ode1
     :members:
     :undoc-members:
     :show-inheritance:
 
-:mod:`exponential` Module
--------------------------
+GPy.kern.parts.exponential module
+---------------------------------
 
 .. automodule:: GPy.kern.parts.exponential
     :members:
     :undoc-members:
     :show-inheritance:
 
-:mod:`finite_dimensional` Module
---------------------------------
+GPy.kern.parts.finite_dimensional module
+----------------------------------------
 
 .. automodule:: GPy.kern.parts.finite_dimensional
     :members:
     :undoc-members:
     :show-inheritance:
 
-:mod:`fixed` Module
--------------------
+GPy.kern.parts.fixed module
+---------------------------
 
 .. automodule:: GPy.kern.parts.fixed
     :members:
     :undoc-members:
     :show-inheritance:
 
-:mod:`gibbs` Module
--------------------
+GPy.kern.parts.gibbs module
+---------------------------
 
 .. automodule:: GPy.kern.parts.gibbs
     :members:
     :undoc-members:
     :show-inheritance:
 
-:mod:`hetero` Module
---------------------
+GPy.kern.parts.hetero module
+----------------------------
 
 .. automodule:: GPy.kern.parts.hetero
     :members:
     :undoc-members:
     :show-inheritance:
 
-:mod:`hierarchical` Module
---------------------------
+GPy.kern.parts.hierarchical module
+----------------------------------
 
 .. automodule:: GPy.kern.parts.hierarchical
     :members:
     :undoc-members:
     :show-inheritance:
 
-:mod:`independent_outputs` Module
----------------------------------
+GPy.kern.parts.independent_outputs module
+-----------------------------------------
 
 .. automodule:: GPy.kern.parts.independent_outputs
     :members:
     :undoc-members:
     :show-inheritance:
 
-:mod:`kernpart` Module
-----------------------
+GPy.kern.parts.kernpart module
+------------------------------
 
 .. automodule:: GPy.kern.parts.kernpart
     :members:
     :undoc-members:
     :show-inheritance:
 
-:mod:`linear` Module
---------------------
+GPy.kern.parts.linear module
+----------------------------
 
 .. automodule:: GPy.kern.parts.linear
     :members:
     :undoc-members:
     :show-inheritance:
 
-:mod:`mlp` Module
------------------
+GPy.kern.parts.mlp module
+-------------------------
 
 .. automodule:: GPy.kern.parts.mlp
     :members:
     :undoc-members:
     :show-inheritance:
 
-:mod:`periodic_Matern32` Module
--------------------------------
+GPy.kern.parts.periodic_Matern32 module
+---------------------------------------
 
 .. automodule:: GPy.kern.parts.periodic_Matern32
     :members:
     :undoc-members:
     :show-inheritance:
 
-:mod:`periodic_Matern52` Module
--------------------------------
+GPy.kern.parts.periodic_Matern52 module
+---------------------------------------
 
 .. automodule:: GPy.kern.parts.periodic_Matern52
     :members:
     :undoc-members:
     :show-inheritance:
 
-:mod:`periodic_exponential` Module
-----------------------------------
+GPy.kern.parts.periodic_exponential module
+------------------------------------------
 
 .. automodule:: GPy.kern.parts.periodic_exponential
     :members:
     :undoc-members:
     :show-inheritance:
 
-:mod:`poly` Module
-------------------
+GPy.kern.parts.poly module
+--------------------------
 
 .. automodule:: GPy.kern.parts.poly
     :members:
     :undoc-members:
     :show-inheritance:
 
-:mod:`prod` Module
-------------------
+GPy.kern.parts.prod module
+--------------------------
 
 .. automodule:: GPy.kern.parts.prod
     :members:
     :undoc-members:
     :show-inheritance:
 
-:mod:`prod_orthogonal` Module
------------------------------
+GPy.kern.parts.prod_orthogonal module
+-------------------------------------
 
 .. automodule:: GPy.kern.parts.prod_orthogonal
     :members:
     :undoc-members:
     :show-inheritance:
 
-:mod:`rational_quadratic` Module
---------------------------------
+GPy.kern.parts.rational_quadratic module
+----------------------------------------
 
 .. automodule:: GPy.kern.parts.rational_quadratic
     :members:
     :undoc-members:
     :show-inheritance:
 
-:mod:`rbf` Module
------------------
+GPy.kern.parts.rbf module
+-------------------------
 
 .. automodule:: GPy.kern.parts.rbf
     :members:
     :undoc-members:
     :show-inheritance:
 
-:mod:`rbf_inv` Module
----------------------
+GPy.kern.parts.rbf_inv module
+-----------------------------
 
 .. automodule:: GPy.kern.parts.rbf_inv
     :members:
     :undoc-members:
     :show-inheritance:
 
-:mod:`rbfcos` Module
---------------------
+GPy.kern.parts.rbfcos module
+----------------------------
 
 .. automodule:: GPy.kern.parts.rbfcos
     :members:
     :undoc-members:
     :show-inheritance:
 
-:mod:`spline` Module
---------------------
+GPy.kern.parts.spline module
+----------------------------
 
 .. automodule:: GPy.kern.parts.spline
     :members:
     :undoc-members:
     :show-inheritance:
 
-:mod:`symmetric` Module
------------------------
+GPy.kern.parts.symmetric module
+-------------------------------
 
 .. automodule:: GPy.kern.parts.symmetric
     :members:
     :undoc-members:
     :show-inheritance:
 
-:mod:`sympy_helpers` Module
----------------------------
+GPy.kern.parts.sympy_helpers module
+-----------------------------------
 
 .. automodule:: GPy.kern.parts.sympy_helpers
     :members:
     :undoc-members:
     :show-inheritance:
 
-:mod:`sympykern` Module
------------------------
+GPy.kern.parts.sympykern module
+-------------------------------
 
 .. automodule:: GPy.kern.parts.sympykern
     :members:
     :undoc-members:
     :show-inheritance:
 
-:mod:`white` Module
--------------------
+GPy.kern.parts.white module
+---------------------------
 
 .. automodule:: GPy.kern.parts.white
     :members:
     :undoc-members:
     :show-inheritance:
 
+
+Module contents
+---------------
+
+.. automodule:: GPy.kern.parts
+    :members:
+    :undoc-members:
+    :show-inheritance:
diff --git a/doc/GPy.kern.rst b/doc/GPy.kern.rst
index 35d9ec00..b4b9d9aa 100644
--- a/doc/GPy.kern.rst
+++ b/doc/GPy.kern.rst
@@ -1,29 +1,5 @@
-kern Package
-============
-
-:mod:`kern` Package
--------------------
-
-.. automodule:: GPy.kern
-    :members:
-    :undoc-members:
-    :show-inheritance:
-
-:mod:`constructors` Module
---------------------------
-
-.. automodule:: GPy.kern.constructors
-    :members:
-    :undoc-members:
-    :show-inheritance:
-
-:mod:`kern` Module
-------------------
-
-.. automodule:: GPy.kern.kern
-    :members:
-    :undoc-members:
-    :show-inheritance:
+GPy.kern package
+================
 
 Subpackages
 -----------
@@ -32,3 +8,30 @@ Subpackages
 
     GPy.kern.parts
 
+Submodules
+----------
+
+GPy.kern.constructors module
+----------------------------
+
+.. automodule:: GPy.kern.constructors
+    :members:
+    :undoc-members:
+    :show-inheritance:
+
+GPy.kern.kern module
+--------------------
+
+.. automodule:: GPy.kern.kern
+    :members:
+    :undoc-members:
+    :show-inheritance:
+
+
+Module contents
+---------------
+
+.. automodule:: GPy.kern
+    :members:
+    :undoc-members:
+    :show-inheritance:
diff --git a/doc/GPy.likelihoods.noise_models.rst b/doc/GPy.likelihoods.noise_models.rst
index 19e5e9fe..6fec5aff 100644
--- a/doc/GPy.likelihoods.noise_models.rst
+++ b/doc/GPy.likelihoods.noise_models.rst
@@ -1,75 +1,78 @@
-noise_models Package
-====================
+GPy.likelihoods.noise_models package
+====================================
 
-:mod:`noise_models` Package
----------------------------
+Submodules
+----------
 
-.. automodule:: GPy.likelihoods.noise_models
-    :members:
-    :undoc-members:
-    :show-inheritance:
-
-:mod:`bernoulli_noise` Module
------------------------------
+GPy.likelihoods.noise_models.bernoulli_noise module
+---------------------------------------------------
 
 .. automodule:: GPy.likelihoods.noise_models.bernoulli_noise
     :members:
     :undoc-members:
     :show-inheritance:
 
-:mod:`exponential_noise` Module
--------------------------------
+GPy.likelihoods.noise_models.exponential_noise module
+-----------------------------------------------------
 
 .. automodule:: GPy.likelihoods.noise_models.exponential_noise
     :members:
     :undoc-members:
     :show-inheritance:
 
-:mod:`gamma_noise` Module
--------------------------
+GPy.likelihoods.noise_models.gamma_noise module
+-----------------------------------------------
 
 .. automodule:: GPy.likelihoods.noise_models.gamma_noise
     :members:
     :undoc-members:
     :show-inheritance:
 
-:mod:`gaussian_noise` Module
-----------------------------
+GPy.likelihoods.noise_models.gaussian_noise module
+--------------------------------------------------
 
 .. automodule:: GPy.likelihoods.noise_models.gaussian_noise
     :members:
     :undoc-members:
     :show-inheritance:
 
-:mod:`gp_transformations` Module
---------------------------------
+GPy.likelihoods.noise_models.gp_transformations module
+------------------------------------------------------
 
 .. automodule:: GPy.likelihoods.noise_models.gp_transformations
     :members:
     :undoc-members:
     :show-inheritance:
 
-:mod:`noise_distributions` Module
----------------------------------
+GPy.likelihoods.noise_models.noise_distributions module
+-------------------------------------------------------
 
 .. automodule:: GPy.likelihoods.noise_models.noise_distributions
     :members:
     :undoc-members:
     :show-inheritance:
 
-:mod:`poisson_noise` Module
----------------------------
+GPy.likelihoods.noise_models.poisson_noise module
+-------------------------------------------------
 
 .. automodule:: GPy.likelihoods.noise_models.poisson_noise
     :members:
     :undoc-members:
     :show-inheritance:
 
-:mod:`student_t_noise` Module
------------------------------
+GPy.likelihoods.noise_models.student_t_noise module
+---------------------------------------------------
 
 .. automodule:: GPy.likelihoods.noise_models.student_t_noise
     :members:
     :undoc-members:
     :show-inheritance:
 
+
+Module contents
+---------------
+
+.. automodule:: GPy.likelihoods.noise_models
+    :members:
+    :undoc-members:
+    :show-inheritance:
diff --git a/doc/GPy.likelihoods.rst b/doc/GPy.likelihoods.rst
index 5dcabbd1..34d98739 100644
--- a/doc/GPy.likelihoods.rst
+++ b/doc/GPy.likelihoods.rst
@@ -1,69 +1,5 @@
-likelihoods Package
-===================
-
-:mod:`likelihoods` Package
---------------------------
-
-.. automodule:: GPy.likelihoods
-    :members:
-    :undoc-members:
-    :show-inheritance:
-
-:mod:`ep` Module
-----------------
-
-.. automodule:: GPy.likelihoods.ep
-    :members:
-    :undoc-members:
-    :show-inheritance:
-
-:mod:`ep_mixed_noise` Module
-----------------------------
-
-.. automodule:: GPy.likelihoods.ep_mixed_noise
-    :members:
-    :undoc-members:
-    :show-inheritance:
-
-:mod:`gaussian` Module
-----------------------
-
-.. automodule:: GPy.likelihoods.gaussian
-    :members:
-    :undoc-members:
-    :show-inheritance:
-
-:mod:`gaussian_mixed_noise` Module
-----------------------------------
-
-.. automodule:: GPy.likelihoods.gaussian_mixed_noise
-    :members:
-    :undoc-members:
-    :show-inheritance:
-
-:mod:`laplace` Module
----------------------
-
-.. automodule:: GPy.likelihoods.laplace
-    :members:
-    :undoc-members:
-    :show-inheritance:
-
-:mod:`likelihood` Module
-------------------------
-
-.. automodule:: GPy.likelihoods.likelihood
-    :members:
-    :undoc-members:
-    :show-inheritance:
-
-:mod:`noise_model_constructors` Module
---------------------------------------
-
-.. automodule:: GPy.likelihoods.noise_model_constructors
-    :members:
-    :undoc-members:
-    :show-inheritance:
+GPy.likelihoods package
+=======================
 
 Subpackages
 -----------
@@ -72,3 +8,70 @@ Subpackages
 
     GPy.likelihoods.noise_models
 
+Submodules
+----------
+
+GPy.likelihoods.ep module
+-------------------------
+
+.. automodule:: GPy.likelihoods.ep
+    :members:
+    :undoc-members:
+    :show-inheritance:
+
+GPy.likelihoods.ep_mixed_noise module
+-------------------------------------
+
+.. automodule:: GPy.likelihoods.ep_mixed_noise
+    :members:
+    :undoc-members:
+    :show-inheritance:
+
+GPy.likelihoods.gaussian module
+-------------------------------
+
+.. automodule:: GPy.likelihoods.gaussian
+    :members:
+    :undoc-members:
+    :show-inheritance:
+
+GPy.likelihoods.gaussian_mixed_noise module
+-------------------------------------------
+
+.. automodule:: GPy.likelihoods.gaussian_mixed_noise
+    :members:
+    :undoc-members:
+    :show-inheritance:
+
+GPy.likelihoods.laplace module
+------------------------------
+
+.. automodule:: GPy.likelihoods.laplace
+    :members:
+    :undoc-members:
+    :show-inheritance:
+
+GPy.likelihoods.likelihood module
+---------------------------------
+
+.. automodule:: GPy.likelihoods.likelihood
+    :members:
+    :undoc-members:
+    :show-inheritance:
+
+GPy.likelihoods.noise_model_constructors module
+-----------------------------------------------
+
+.. automodule:: GPy.likelihoods.noise_model_constructors
+    :members:
+    :undoc-members:
+    :show-inheritance:
+
+
+Module contents
+---------------
+
+.. automodule:: GPy.likelihoods
+    :members:
+    :undoc-members:
+    :show-inheritance:
diff --git a/doc/GPy.mappings.rst b/doc/GPy.mappings.rst
index b7444808..c48cb06e 100644
--- a/doc/GPy.mappings.rst
+++ b/doc/GPy.mappings.rst
@@ -1,35 +1,38 @@
-mappings Package
-================
+GPy.mappings package
+====================
 
-:mod:`mappings` Package
------------------------
+Submodules
+----------
 
-.. automodule:: GPy.mappings
-    :members:
-    :undoc-members:
-    :show-inheritance:
-
-:mod:`kernel` Module
---------------------
+GPy.mappings.kernel module
+--------------------------
 
 .. automodule:: GPy.mappings.kernel
     :members:
     :undoc-members:
     :show-inheritance:
 
-:mod:`linear` Module
---------------------
+GPy.mappings.linear module
+--------------------------
 
 .. automodule:: GPy.mappings.linear
     :members:
     :undoc-members:
     :show-inheritance:
 
-:mod:`mlp` Module
------------------
+GPy.mappings.mlp module
+-----------------------
 
 .. automodule:: GPy.mappings.mlp
     :members:
     :undoc-members:
     :show-inheritance:
 
+
+Module contents
+---------------
+
+.. automodule:: GPy.mappings
+    :members:
+    :undoc-members:
+    :show-inheritance:
diff --git a/doc/GPy.models_modules.rst b/doc/GPy.models_modules.rst
index 4169ec3a..c16941b1 100644
--- a/doc/GPy.models_modules.rst
+++ b/doc/GPy.models_modules.rst
@@ -1,131 +1,134 @@
-models_modules Package
-======================
+GPy.models_modules package
+==========================
 
-:mod:`models_modules` Package
------------------------------
+Submodules
+----------
 
-.. automodule:: GPy.models_modules
-    :members:
-    :undoc-members:
-    :show-inheritance:
-
-:mod:`bayesian_gplvm` Module
-----------------------------
+GPy.models_modules.bayesian_gplvm module
+----------------------------------------
 
 .. automodule:: GPy.models_modules.bayesian_gplvm
     :members:
     :undoc-members:
     :show-inheritance:
 
-:mod:`bcgplvm` Module
----------------------
+GPy.models_modules.bcgplvm module
+---------------------------------
 
 .. automodule:: GPy.models_modules.bcgplvm
     :members:
     :undoc-members:
     :show-inheritance:
 
-:mod:`fitc_classification` Module
----------------------------------
+GPy.models_modules.fitc_classification module
+---------------------------------------------
 
 .. automodule:: GPy.models_modules.fitc_classification
     :members:
     :undoc-members:
     :show-inheritance:
 
-:mod:`gp_classification` Module
--------------------------------
+GPy.models_modules.gp_classification module
+-------------------------------------------
 
 .. automodule:: GPy.models_modules.gp_classification
     :members:
     :undoc-members:
     :show-inheritance:
 
-:mod:`gp_multioutput_regression` Module
----------------------------------------
+GPy.models_modules.gp_multioutput_regression module
+---------------------------------------------------
 
 .. automodule:: GPy.models_modules.gp_multioutput_regression
     :members:
     :undoc-members:
     :show-inheritance:
 
-:mod:`gp_regression` Module
----------------------------
+GPy.models_modules.gp_regression module
+---------------------------------------
 
 .. automodule:: GPy.models_modules.gp_regression
     :members:
     :undoc-members:
     :show-inheritance:
 
-:mod:`gplvm` Module
--------------------
+GPy.models_modules.gplvm module
+-------------------------------
 
 .. automodule:: GPy.models_modules.gplvm
     :members:
     :undoc-members:
     :show-inheritance:
 
-:mod:`gradient_checker` Module
-------------------------------
+GPy.models_modules.gradient_checker module
+------------------------------------------
 
 .. automodule:: GPy.models_modules.gradient_checker
     :members:
     :undoc-members:
     :show-inheritance:
 
-:mod:`mrd` Module
------------------
+GPy.models_modules.mrd module
+-----------------------------
 
 .. automodule:: GPy.models_modules.mrd
     :members:
     :undoc-members:
     :show-inheritance:
 
-:mod:`sparse_gp_classification` Module
---------------------------------------
+GPy.models_modules.sparse_gp_classification module
+--------------------------------------------------
 
 .. automodule:: GPy.models_modules.sparse_gp_classification
     :members:
     :undoc-members:
     :show-inheritance:
 
-:mod:`sparse_gp_multioutput_regression` Module
-----------------------------------------------
+GPy.models_modules.sparse_gp_multioutput_regression module
+----------------------------------------------------------
 
 .. automodule:: GPy.models_modules.sparse_gp_multioutput_regression
     :members:
     :undoc-members:
     :show-inheritance:
 
-:mod:`sparse_gp_regression` Module
-----------------------------------
+GPy.models_modules.sparse_gp_regression module
+----------------------------------------------
 
 .. automodule:: GPy.models_modules.sparse_gp_regression
     :members:
     :undoc-members:
     :show-inheritance:
 
-:mod:`sparse_gplvm` Module
---------------------------
+GPy.models_modules.sparse_gplvm module
+--------------------------------------
 
 .. automodule:: GPy.models_modules.sparse_gplvm
     :members:
     :undoc-members:
     :show-inheritance:
 
-:mod:`svigp_regression` Module
-------------------------------
+GPy.models_modules.svigp_regression module
+------------------------------------------
 
 .. automodule:: GPy.models_modules.svigp_regression
     :members:
     :undoc-members:
     :show-inheritance:
 
-:mod:`warped_gp` Module
------------------------
+GPy.models_modules.warped_gp module
+-----------------------------------
 
 .. automodule:: GPy.models_modules.warped_gp
     :members:
     :undoc-members:
     :show-inheritance:
 
+
+Module contents
+---------------
+
+.. automodule:: GPy.models_modules
+    :members:
+    :undoc-members:
+    :show-inheritance:
diff --git a/doc/GPy.rst b/doc/GPy.rst
index 31ec3562..cd1afd29 100644
--- a/doc/GPy.rst
+++ b/doc/GPy.rst
@@ -1,22 +1,6 @@
-GPy Package
+GPy package
 ===========
 
-:mod:`GPy` Package
-------------------
-
-.. automodule:: GPy.__init__
-    :members:
-    :undoc-members:
-    :show-inheritance:
-
-:mod:`models` Module
---------------------
-
-.. automodule:: GPy.models
-    :members:
-    :undoc-members:
-    :show-inheritance:
-
 Subpackages
 -----------
 
@@ -32,3 +16,22 @@ Subpackages
     GPy.testing
     GPy.util
 
+Submodules
+----------
+
+GPy.models module
+-----------------
+
+.. automodule:: GPy.models
+    :members:
+    :undoc-members:
+    :show-inheritance:
+
+
+Module contents
+---------------
+
+.. automodule:: GPy
+    :members:
+    :undoc-members:
+    :show-inheritance:
diff --git a/doc/GPy.testing.rst b/doc/GPy.testing.rst
index 15b0cc79..fcf9dc30 100644
--- a/doc/GPy.testing.rst
+++ b/doc/GPy.testing.rst
@@ -1,131 +1,134 @@
-testing Package
-===============
+GPy.testing package
+===================
 
-:mod:`testing` Package
-----------------------
+Submodules
+----------
 
-.. automodule:: GPy.testing
-    :members:
-    :undoc-members:
-    :show-inheritance:
-
-:mod:`bcgplvm_tests` Module
----------------------------
+GPy.testing.bcgplvm_tests module
+--------------------------------
 
 .. automodule:: GPy.testing.bcgplvm_tests
     :members:
     :undoc-members:
     :show-inheritance:
 
-:mod:`bgplvm_tests` Module
---------------------------
+GPy.testing.bgplvm_tests module
+-------------------------------
 
 .. automodule:: GPy.testing.bgplvm_tests
     :members:
     :undoc-members:
     :show-inheritance:
 
-:mod:`cgd_tests` Module
------------------------
+GPy.testing.cgd_tests module
+----------------------------
 
 .. automodule:: GPy.testing.cgd_tests
     :members:
     :undoc-members:
     :show-inheritance:
 
-:mod:`examples_tests` Module
-----------------------------
+GPy.testing.examples_tests module
+---------------------------------
 
 .. automodule:: GPy.testing.examples_tests
     :members:
     :undoc-members:
     :show-inheritance:
 
-:mod:`gp_transformation_tests` Module
--------------------------------------
+GPy.testing.gp_transformation_tests module
+------------------------------------------
 
 .. automodule:: GPy.testing.gp_transformation_tests
     :members:
     :undoc-members:
     :show-inheritance:
 
-:mod:`gplvm_tests` Module
--------------------------
+GPy.testing.gplvm_tests module
+------------------------------
 
 .. automodule:: GPy.testing.gplvm_tests
     :members:
     :undoc-members:
     :show-inheritance:
 
-:mod:`kernel_tests` Module
---------------------------
+GPy.testing.kernel_tests module
+-------------------------------
 
 .. automodule:: GPy.testing.kernel_tests
     :members:
     :undoc-members:
     :show-inheritance:
 
-:mod:`likelihoods_tests` Module
--------------------------------
+GPy.testing.likelihood_tests module
+-----------------------------------
 
-.. automodule:: GPy.testing.likelihoods_tests
+.. automodule:: GPy.testing.likelihood_tests
     :members:
     :undoc-members:
     :show-inheritance:
 
-:mod:`mapping_tests` Module
----------------------------
+GPy.testing.mapping_tests module
+--------------------------------
 
 .. automodule:: GPy.testing.mapping_tests
     :members:
     :undoc-members:
     :show-inheritance:
 
-:mod:`mrd_tests` Module
------------------------
+GPy.testing.mrd_tests module
+----------------------------
 
 .. automodule:: GPy.testing.mrd_tests
     :members:
     :undoc-members:
     :show-inheritance:
 
-:mod:`prior_tests` Module
--------------------------
+GPy.testing.prior_tests module
+------------------------------
 
 .. automodule:: GPy.testing.prior_tests
     :members:
     :undoc-members:
     :show-inheritance:
 
-:mod:`psi_stat_expectation_tests` Module
-----------------------------------------
+GPy.testing.psi_stat_expectation_tests module
+---------------------------------------------
 
 .. automodule:: GPy.testing.psi_stat_expectation_tests
     :members:
     :undoc-members:
     :show-inheritance:
 
-:mod:`psi_stat_gradient_tests` Module
--------------------------------------
+GPy.testing.psi_stat_gradient_tests module
+------------------------------------------
 
 .. automodule:: GPy.testing.psi_stat_gradient_tests
     :members:
     :undoc-members:
     :show-inheritance:
 
-:mod:`sparse_gplvm_tests` Module
---------------------------------
+GPy.testing.sparse_gplvm_tests module
+-------------------------------------
 
 .. automodule:: GPy.testing.sparse_gplvm_tests
     :members:
     :undoc-members:
     :show-inheritance:
 
-:mod:`unit_tests` Module
-------------------------
+GPy.testing.unit_tests module
+-----------------------------
 
 .. automodule:: GPy.testing.unit_tests
     :members:
     :undoc-members:
     :show-inheritance:
 
+
+Module contents
+---------------
+
+.. automodule:: GPy.testing
+    :members:
+    :undoc-members:
+    :show-inheritance:
diff --git a/doc/GPy.util.latent_space_visualizations.controllers.rst b/doc/GPy.util.latent_space_visualizations.controllers.rst
index e78ade7b..a88c1f5c 100644
--- a/doc/GPy.util.latent_space_visualizations.controllers.rst
+++ b/doc/GPy.util.latent_space_visualizations.controllers.rst
@@ -1,27 +1,30 @@
-controllers Package
-===================
+GPy.util.latent_space_visualizations.controllers package
+========================================================
 
-:mod:`controllers` Package
---------------------------
+Submodules
+----------
 
-.. automodule:: GPy.util.latent_space_visualizations.controllers
-    :members:
-    :undoc-members:
-    :show-inheritance:
-
-:mod:`axis_event_controller` Module
------------------------------------
+GPy.util.latent_space_visualizations.controllers.axis_event_controller module
+-----------------------------------------------------------------------------
 
 .. automodule:: GPy.util.latent_space_visualizations.controllers.axis_event_controller
     :members:
     :undoc-members:
     :show-inheritance:
 
-:mod:`imshow_controller` Module
--------------------------------
+GPy.util.latent_space_visualizations.controllers.imshow_controller module
+-------------------------------------------------------------------------
 
 .. automodule:: GPy.util.latent_space_visualizations.controllers.imshow_controller
     :members:
     :undoc-members:
     :show-inheritance:
 
+
+Module contents
+---------------
+
+.. automodule:: GPy.util.latent_space_visualizations.controllers
+    :members:
+    :undoc-members:
+    :show-inheritance:
diff --git a/doc/GPy.util.latent_space_visualizations.rst b/doc/GPy.util.latent_space_visualizations.rst
index 4b440f61..d8cbd843 100644
--- a/doc/GPy.util.latent_space_visualizations.rst
+++ b/doc/GPy.util.latent_space_visualizations.rst
@@ -1,13 +1,5 @@
-latent_space_visualizations Package
-===================================
-
-:mod:`latent_space_visualizations` Package
-------------------------------------------
-
-.. automodule:: GPy.util.latent_space_visualizations
-    :members:
-    :undoc-members:
-    :show-inheritance:
+GPy.util.latent_space_visualizations package
+============================================
 
 Subpackages
 -----------
@@ -15,5 +7,11 @@ Subpackages
 .. toctree::
 
     GPy.util.latent_space_visualizations.controllers
-    GPy.util.latent_space_visualizations.views
 
+Module contents
+---------------
+
+.. automodule:: GPy.util.latent_space_visualizations
+    :members:
+    :undoc-members:
+    :show-inheritance:
diff --git a/doc/GPy.util.rst b/doc/GPy.util.rst
index 2e20c006..1c35a7ba 100644
--- a/doc/GPy.util.rst
+++ b/doc/GPy.util.rst
@@ -1,181 +1,5 @@
-util Package
-============
-
-:mod:`util` Package
--------------------
-
-.. automodule:: GPy.util
-    :members:
-    :undoc-members:
-    :show-inheritance:
-
-:mod:`Tango` Module
--------------------
-
-.. automodule:: GPy.util.Tango
-    :members:
-    :undoc-members:
-    :show-inheritance:
-
-:mod:`block_matrices` Module
-----------------------------
-
-.. automodule:: GPy.util.block_matrices
-    :members:
-    :undoc-members:
-    :show-inheritance:
-
-:mod:`classification` Module
-----------------------------
-
-.. automodule:: GPy.util.classification
-    :members:
-    :undoc-members:
-    :show-inheritance:
-
-:mod:`config` Module
---------------------
-
-.. automodule:: GPy.util.config
-    :members:
-    :undoc-members:
-    :show-inheritance:
-
-:mod:`datasets` Module
-----------------------
-
-.. automodule:: GPy.util.datasets
-    :members:
-    :undoc-members:
-    :show-inheritance:
-
-:mod:`decorators` Module
-------------------------
-
-.. automodule:: GPy.util.decorators
-    :members:
-    :undoc-members:
-    :show-inheritance:
-
-:mod:`erfcx` Module
--------------------
-
-.. automodule:: GPy.util.erfcx
-    :members:
-    :undoc-members:
-    :show-inheritance:
-
-:mod:`linalg` Module
---------------------
-
-.. automodule:: GPy.util.linalg
-    :members:
-    :undoc-members:
-    :show-inheritance:
-
-:mod:`ln_diff_erfs` Module
---------------------------
-
-.. automodule:: GPy.util.ln_diff_erfs
-    :members:
-    :undoc-members:
-    :show-inheritance:
-
-:mod:`misc` Module
-------------------
-
-.. automodule:: GPy.util.misc
-    :members:
-    :undoc-members:
-    :show-inheritance:
-
-:mod:`mocap` Module
--------------------
-
-.. automodule:: GPy.util.mocap
-    :members:
-    :undoc-members:
-    :show-inheritance:
-
-:mod:`multioutput` Module
--------------------------
-
-.. automodule:: GPy.util.multioutput
-    :members:
-    :undoc-members:
-    :show-inheritance:
-
-:mod:`netpbmfile` Module
-------------------------
-
-.. automodule:: GPy.util.netpbmfile
-    :members:
-    :undoc-members:
-    :show-inheritance:
-
-:mod:`pca` Module
------------------
-
-.. automodule:: GPy.util.pca
-    :members:
-    :undoc-members:
-    :show-inheritance:
-
-:mod:`plot` Module
-------------------
-
-.. automodule:: GPy.util.plot
-    :members:
-    :undoc-members:
-    :show-inheritance:
-
-:mod:`plot_latent` Module
--------------------------
-
-.. automodule:: GPy.util.plot_latent
-    :members:
-    :undoc-members:
-    :show-inheritance:
-
-:mod:`squashers` Module
------------------------
-
-.. automodule:: GPy.util.squashers
-    :members:
-    :undoc-members:
-    :show-inheritance:
-
-:mod:`symbolic` Module
-----------------------
-
-.. automodule:: GPy.util.symbolic
-    :members:
-    :undoc-members:
-    :show-inheritance:
-
-:mod:`univariate_Gaussian` Module
----------------------------------
-
-.. automodule:: GPy.util.univariate_Gaussian
-    :members:
-    :undoc-members:
-    :show-inheritance:
-
-:mod:`visualize` Module
------------------------
-
-.. automodule:: GPy.util.visualize
-    :members:
-    :undoc-members:
-    :show-inheritance:
-
-:mod:`warping_functions` Module
--------------------------------
-
-.. automodule:: GPy.util.warping_functions
-    :members:
-    :undoc-members:
-    :show-inheritance:
+GPy.util package
+================
 
 Subpackages
 -----------
@@ -184,3 +8,190 @@ Subpackages
 
     GPy.util.latent_space_visualizations
 
+Submodules
+----------
+
+GPy.util.Tango module
+---------------------
+
+.. automodule:: GPy.util.Tango
+    :members:
+    :undoc-members:
+    :show-inheritance:
+
+GPy.util.block_matrices module
+------------------------------
+
+.. automodule:: GPy.util.block_matrices
+    :members:
+    :undoc-members:
+    :show-inheritance:
+
+GPy.util.classification module
+------------------------------
+
+.. automodule:: GPy.util.classification
+    :members:
+    :undoc-members:
+    :show-inheritance:
+
+GPy.util.config module
+----------------------
+
+.. automodule:: GPy.util.config
+    :members:
+    :undoc-members:
+    :show-inheritance:
+
+GPy.util.datasets module
+------------------------
+
+.. automodule:: GPy.util.datasets
+    :members:
+    :undoc-members:
+    :show-inheritance:
+
+GPy.util.decorators module
+--------------------------
+
+.. automodule:: GPy.util.decorators
+    :members:
+    :undoc-members:
+    :show-inheritance:
+
+GPy.util.diag module
+--------------------
+
+.. automodule:: GPy.util.diag
+    :members:
+    :undoc-members:
+    :show-inheritance:
+
+GPy.util.erfcx module
+---------------------
+
+.. automodule:: GPy.util.erfcx
+    :members:
+    :undoc-members:
+    :show-inheritance:
+
+GPy.util.linalg module
+----------------------
+
+.. automodule:: GPy.util.linalg
+    :members:
+    :undoc-members:
+    :show-inheritance:
+
+GPy.util.ln_diff_erfs module
+----------------------------
+
+.. automodule:: GPy.util.ln_diff_erfs
+    :members:
+    :undoc-members:
+    :show-inheritance:
+
+GPy.util.misc module
+--------------------
+
+.. automodule:: GPy.util.misc
+    :members:
+    :undoc-members:
+    :show-inheritance:
+
+GPy.util.mocap module
+---------------------
+
+.. automodule:: GPy.util.mocap
+    :members:
+    :undoc-members:
+    :show-inheritance:
+
+GPy.util.multioutput module
+---------------------------
+
+.. automodule:: GPy.util.multioutput
+    :members:
+    :undoc-members:
+    :show-inheritance:
+
+GPy.util.netpbmfile module
+--------------------------
+
+.. automodule:: GPy.util.netpbmfile
+    :members:
+    :undoc-members:
+    :show-inheritance:
+
+GPy.util.plot module
+--------------------
+
+.. automodule:: GPy.util.plot
+    :members:
+    :undoc-members:
+    :show-inheritance:
+
+GPy.util.plot_latent module
+---------------------------
+
+.. automodule:: GPy.util.plot_latent
+    :members:
+    :undoc-members:
+    :show-inheritance:
+
+GPy.util.squashers module
+-------------------------
+
+.. automodule:: GPy.util.squashers
+    :members:
+    :undoc-members:
+    :show-inheritance:
+
+GPy.util.subarray_and_sorting module
+------------------------------------
+
+.. automodule:: GPy.util.subarray_and_sorting
+    :members:
+    :undoc-members:
+    :show-inheritance:
+
+GPy.util.symbolic module
+------------------------
+
+.. automodule:: GPy.util.symbolic
+    :members:
+    :undoc-members:
+    :show-inheritance:
+
+GPy.util.univariate_Gaussian module
+-----------------------------------
+
+.. automodule:: GPy.util.univariate_Gaussian
+    :members:
+    :undoc-members:
+    :show-inheritance:
+
+GPy.util.visualize module
+-------------------------
+
+.. automodule:: GPy.util.visualize
+    :members:
+    :undoc-members:
+    :show-inheritance:
+
+GPy.util.warping_functions module
+---------------------------------
+
+.. automodule:: GPy.util.warping_functions
+    :members:
+    :undoc-members:
+    :show-inheritance:
+
+
+Module contents
+---------------
+
+.. automodule:: GPy.util
+    :members:
+    :undoc-members:
+    :show-inheritance:
diff --git a/setup.py b/setup.py
index 3b493022..80f3648a 100644
--- a/setup.py
+++ b/setup.py
@@ -20,7 +20,7 @@ setup(name = 'GPy',
       url = "http://sheffieldml.github.com/GPy/",
       packages = ['GPy', 'GPy.core', 'GPy.kern', 'GPy.util', 'GPy.models_modules', 'GPy.inference', 'GPy.examples', 'GPy.likelihoods', 'GPy.testing', 'GPy.util.latent_space_visualizations', 'GPy.util.latent_space_visualizations.controllers', 'GPy.likelihoods.noise_models', 'GPy.kern.parts', 'GPy.mappings'],
       package_dir={'GPy': 'GPy'},
-      package_data = {'GPy': ['GPy/examples', 'gpy_config.cfg']},
+      package_data = {'GPy': ['GPy/examples', 'gpy_config.cfg', 'util/data_resources.json']},
       py_modules = ['GPy.__init__'],
       long_description=read('README.md'),
       install_requires=['numpy>=1.6', 'scipy>=0.9','matplotlib>=1.1', 'nose'],