Merge branch 'master' of github.com:SheffieldML/GPy

GPy/examples/regression.py

@@ -194,7 +194,7 @@ def multiple_optima(gene_number=937,resolution=80, model_restarts=10, seed=10000
     # Remove the mean (no bias kernel to ensure signal/noise is in RBF/white)
     data['Y'] = data['Y'] - np.mean(data['Y'])
 
-    lls = GPy.examples.regression.contour_data(data, length_scales, log_SNRs, GPy.kern.rbf)
+    lls = GPy.examples.regression._contour_data(data, length_scales, log_SNRs, GPy.kern.rbf)
     pb.contour(length_scales, log_SNRs, np.exp(lls), 20)
     ax = pb.gca()
     pb.xlabel('length scale')
@@ -229,7 +229,7 @@ def multiple_optima(gene_number=937,resolution=80, model_restarts=10, seed=10000
     ax.set_ylim(ylim)
     return (models, lls)
 
-def contour_data(data, length_scales, log_SNRs, signal_kernel_call=GPy.kern.rbf):
+def _contour_data(data, length_scales, log_SNRs, signal_kernel_call=GPy.kern.rbf):
     """Evaluate the GP objective function for a given data set for a range of signal to noise ratios and a range of lengthscales.
 
     :data_set: A data set from the utils.datasets director.

GPy/examples/tutorials.py

@@ -6,14 +6,14 @@
 Code of Tutorials
 """
 
+import pylab as pb
+pb.ion()
+import numpy as np
+import GPy
+
 def tuto_GP_regression():
     """The detailed explanations of the commands used in this file can be found in the tutorial section"""
 
-    import pylab as pb
-    pb.ion()
-    import numpy as np
-    import GPy
-
     X = np.random.uniform(-3.,3.,(20,1))
     Y = np.sin(X) + np.random.randn(20,1)*0.05
 
@@ -39,11 +39,6 @@ def tuto_GP_regression():
     #  2-dimensional example  #
     ###########################
 
-    import pylab as pb
-    pb.ion()
-    import numpy as np
-    import GPy
-
     # sample inputs and outputs
     X = np.random.uniform(-3.,3.,(50,2))
     Y = np.sin(X[:,0:1]) * np.sin(X[:,1:2])+np.random.randn(50,1)*0.05
@@ -67,9 +62,6 @@ def tuto_GP_regression():
 
 def tuto_kernel_overview():
     """The detailed explanations of the commands used in this file can be found in the tutorial section"""
-    import pylab as pb
-    import numpy as np
-    import GPy
     pb.ion()
 
     ker1 = GPy.kern.rbf(1)  # Equivalent to ker1 = GPy.kern.rbf(D=1, variance=1., lengthscale=1.)

GPy/kern/rational_quadratic.py

@@ -11,7 +11,7 @@ class rational_quadratic(kernpart):
 
     .. math::
 
-       k(r) = \sigma^2 \left(1 + \frac{r^2}{2 \ell^2})^{- \alpha} \ \ \ \ \  \\text{ where  } r^2 = (x-y)^2
+       k(r) = \sigma^2 \\bigg( 1 + \\frac{r^2}{2 \ell^2} \\bigg)^{- \\alpha} \ \ \ \ \  \\text{ where  } r^2 = (x-y)^2
 
     :param D: the number of input dimensions
     :type D: int (D=1 is the only value currently supported)
@@ -19,6 +19,8 @@ class rational_quadratic(kernpart):
     :type variance: float
     :param lengthscale: the lengthscale :math:`\ell`
     :type lengthscale: float
+    :param power: the power :math:`\\alpha`
+    :type power: float
     :rtype: kernpart object
 
     """
@@ -76,4 +78,3 @@ class rational_quadratic(kernpart):
 
     def dKdiag_dX(self,dL_dKdiag,X,target):
         pass
-    

GPy/kern/rbf.py

@@ -12,7 +12,7 @@ class rbf(kernpart):
 
     .. math::
 
-       k(r) = \sigma^2 \exp(- \frac{1}{2}r^2) \ \ \ \ \  \\text{ where  } r^2 = \sum_{i=1}^d \frac{ (x_i-x^\prime_i)^2}{\ell_i^2}}
+       k(r) = \sigma^2 \exp \\bigg(- \\frac{1}{2} r^2 \\bigg) \ \ \ \ \  \\text{ where  } r^2 = \sum_{i=1}^d \\frac{ (x_i-x^\prime_i)^2}{\ell_i^2}
 
     where \ell_i is the lengthscale, \sigma^2 the variance and d the dimensionality of the input.
 

GPy/models/Bayesian_GPLVM.py

@@ -83,3 +83,7 @@ class Bayesian_GPLVM(sparse_GP, GPLVM):
 
     def _log_likelihood_gradients(self):
         return np.hstack((self.dL_dmuS().flatten(), sparse_GP._log_likelihood_gradients(self)))
+
+    def plot_latent(self, *args, **kwargs):
+        input_1, input_2 = GPLVM.plot_latent(self, *args, **kwargs)
+        pb.plot(self.Z[:, input_1], self.Z[:, input_2], '^w')

GPy/models/GPLVM.py

@@ -81,13 +81,16 @@ class GPLVM(GP):
                     raise ValueError, "cannot Atomatically determine which dimensions to plot, please pass 'which_indices'"
                 k = k[0]
                 if k.name=='rbf':
-                    input_1, input_2 = np.argsort(k.lengthscales)[:2]
+                    input_1, input_2 = np.argsort(k.lengthscale)[:2]
                 elif k.name=='linear':
                     input_1, input_2 = np.argsort(k.variances)[::-1][:2]
 
         #first, plot the output variance as a function of the latent space
         Xtest, xx,yy,xmin,xmax = util.plot.x_frame2D(self.X[:,[input_1, input_2]],resolution=resolution)
-        mu, var, low, up = self.predict(Xtest)
+	Xtest_full = np.zeros((Xtest.shape[0], self.X.shape[1]))
+	Xtest_full[:, :2] = Xtest        
+	mu, var, low, up = self.predict(Xtest_full)
+	var = var[:, :2]
         pb.imshow(var.reshape(resolution,resolution).T[::-1,:],extent=[xmin[0],xmax[0],xmin[1],xmax[1]],cmap=pb.cm.binary,interpolation='bilinear')
 
 
@@ -117,6 +120,4 @@ class GPLVM(GP):
         pb.xlim(xmin[0],xmax[0])
         pb.ylim(xmin[1],xmax[1])
 
-
+        return input_1, input_2
-
-        

GPy/models/sparse_GPLVM.py

@@ -55,3 +55,7 @@ class sparse_GPLVM(sparse_GP_regression, GPLVM):
         #passing Z without a small amout of jitter will induce the white kernel where we don;t want it!
         mu, var, upper, lower = sparse_GP_regression.predict(self, self.Z+np.random.randn(*self.Z.shape)*0.0001)
         pb.plot(mu[:, 0] , mu[:, 1], 'ko')
+
+    def plot_latent(self, *args, **kwargs):
+        input_1, input_2 = GPLVM.plot_latent(*args, **kwargs)
+        pb.plot(m.Z[:, input_1], m.Z[:, input_2], '^w')

GPy/testing/examples_tests.py

@@ -4,22 +4,73 @@
 import unittest
 import numpy as np
 import GPy
+import inspect
+import pkgutil
+import os
+import random
+
 
 class ExamplesTests(unittest.TestCase):
-    def test_check_model_returned(self):
+    def _checkgrad(self, model):
-        pass
+        self.assertTrue(model.checkgrad())
 
-    def test_model_checkgrads(self):
+    def _model_instance(self, model):
-        pass
+        self.assertTrue(isinstance(model, GPy.models))
 
-    def test_all_examples(self):
+"""
-        pass
+def model_instance_generator(model):
-        #Load models
+    def check_model_returned(self):
+        self._model_instance(model)
+    return check_model_returned
 
-        #Loop through models
+def checkgrads_generator(model):
-        #for model in models:
+    def model_checkgrads(self):
-            #self.assertTrue(m.checkgrad())
+        self._checkgrad(model)
+    return model_checkgrads
+"""
 
+def model_checkgrads(model):
+    model.randomize()
+    assert model.checkgrad()
+
+
+def model_instance(model):
+    assert isinstance(model, GPy.core.model)
+
+
+def test_models():
+    examples_path = os.path.dirname(GPy.examples.__file__)
+    #Load modules
+    for loader, module_name, is_pkg in pkgutil.iter_modules([examples_path]):
+        #Load examples
+        module_examples = loader.find_module(module_name).load_module(module_name)
+        print "MODULE", module_examples
+        print "Before"
+        print inspect.getmembers(module_examples, predicate=inspect.isfunction)
+        functions = [ func for func in inspect.getmembers(module_examples, predicate=inspect.isfunction) if func[0].startswith('_') is False ][::-1]
+        print "After"
+        print functions
+        for example in functions:
+            print "Testing example: ", example[0]
+            #Generate model
+            model = example[1]()
+            print model
+
+            #Create tests for instance check
+            """
+            test = model_instance_generator(model)
+            test.__name__ = 'test_instance_%s' % example[0]
+            setattr(ExamplesTests, test.__name__, test)
+
+            #Create tests for checkgrads check
+            test = checkgrads_generator(model)
+            test.__name__ = 'test_checkgrads_%s' % example[0]
+            setattr(ExamplesTests, test.__name__, test)
+            """
+            model_checkgrads.description = 'test_checkgrads_%s' % example[0]
+            yield model_checkgrads, model
+            model_instance.description = 'test_instance_%s' % example[0]
+            yield model_instance, model
 
 if __name__ == "__main__":
     print "Running unit tests, please be (very) patient..."

doc/GPy.testing.rst

@@ -0,0 +1,59 @@
+testing Package
+===============
+
+:mod:`bgplvm_tests` Module
+--------------------------
+
+.. automodule:: GPy.testing.bgplvm_tests
+    :members:
+    :undoc-members:
+    :show-inheritance:
+
+:mod:`examples_tests` Module
+----------------------------
+
+.. automodule:: GPy.testing.examples_tests
+    :members:
+    :undoc-members:
+    :show-inheritance:
+
+:mod:`gplvm_tests` Module
+-------------------------
+
+.. automodule:: GPy.testing.gplvm_tests
+    :members:
+    :undoc-members:
+    :show-inheritance:
+
+:mod:`kernel_tests` Module
+--------------------------
+
+.. automodule:: GPy.testing.kernel_tests
+    :members:
+    :undoc-members:
+    :show-inheritance:
+
+:mod:`prior_tests` Module
+-------------------------
+
+.. automodule:: GPy.testing.prior_tests
+    :members:
+    :undoc-members:
+    :show-inheritance:
+
+:mod:`sparse_gplvm_tests` Module
+--------------------------------
+
+.. automodule:: GPy.testing.sparse_gplvm_tests
+    :members:
+    :undoc-members:
+    :show-inheritance:
+
+:mod:`unit_tests` Module
+------------------------
+
+.. automodule:: GPy.testing.unit_tests
+    :members:
+    :undoc-members:
+    :show-inheritance:
+

doc/index.rst

@@ -10,8 +10,7 @@ For a quick start, you can have a look at one of the tutorials:
 * `Basic Gaussian process regression <tuto_GP_regression.html>`_  
 * `Interacting with models <tuto_interacting_with_models.html>`_
 * `A kernel overview <tuto_kernel_overview.html>`_ 
-* Advanced GP regression (Forthcoming)
+* `Writing new kernels <tuto_creating_new_kernels.html>`_
-* Writing kernels (Forthcoming)
 
 You may also be interested by some examples in the GPy/examples folder.
 

doc/tuto_creating_new_kernels.rst

@@ -0,0 +1,183 @@
+********************
+Creating new kernels
+********************
+
+We will see in this tutorial how to create new kernels in GPy. We will also give details on how to implement each function of the kernel and illustrate with a running example: the rational quadratic kernel. 
+
+Structure of a kernel in GPy
+============================
+
+In GPy a kernel object is made of a list of kernpart objects, which correspond to symetric positive definite functions. More precisely, the kernel should be understood as the sum of the kernparts. In order to implement a new covariance, the following steps must be followed
+
+    1. implement the new covariance as a kernpart object
+    2. update the constructors that allow to use the kernpart as a kern object
+    3. update the __init__.py file
+
+Theses three steps are detailed below.
+
+Implementing a kernpart object
+==============================
+
+We advise the reader to start with copy-pasting an existing kernel and to modify the new file. We will now give a description of the various functions that can be found in a kernpart object.
+
+**Header**
+
+The header is similar to all kernels: ::
+
+    from kernpart import kernpart
+    import numpy as np
+
+    class rational_quadratic(kernpart):
+
+**__init__(self,D, param1, param2, ...)**
+
+The implementation of this function in mandatory.
+
+For all kernparts the first parameter ``D`` corresponds to the dimension of the input space, and the following parameters stand for the parameterization of the kernel.
+
+The following attributes are compulsory: ``self.D`` (the dimension, integer), ``self.name`` (name of the kernel, string), ``self.Nparam`` (number of parameters, integer). ::
+
+    def __init__(self,D,variance=1.,lengthscale=1.,power=1.):
+        assert D == 1, "For this kernel we assume D=1"
+        self.D = D
+        self.Nparam = 3
+        self.name = 'rat_quad'
+        self.variance = variance
+        self.lengthscale = lengthscale
+        self.power = power
+
+**_get_params(self)**
+
+The implementation of this function in mandatory.
+
+This function returns a one dimensional array of length ``self.Nparam`` containing the value of the parameters. ::
+
+    def _get_params(self):
+        return np.hstack((self.variance,self.lengthscale,self.power))
+
+**_set_params(self,x)**
+
+The implementation of this function in mandatory.
+
+The input is a one dimensional array of length ``self.Nparam`` containing the value of the parameters. The function has no output but it updates the values of the attribute associated to the parameters (such as ``self.variance``, ``self.lengthscale``, ...). ::
+
+    def _set_params(self,x):
+        self.variance = x[0]
+        self.lengthscale = x[1]
+        self.power = x[2]
+
+**_get_param_names(self)**
+
+The implementation of this function in mandatory.
+
+It returns a list of strings of length ``self.Nparam`` corresponding to the parameter names. ::
+
+    def _get_param_names(self):
+        return ['variance','lengthscale','power']
+
+**K(self,X,X2,target)**
+
+The implementation of this function in mandatory.
+
+This function is used to compute the covariance matrix associated with the inputs X, X2 (np.arrays with arbitrary number of line (say :math:`n_1`, :math:`n_2`) and ``self.D`` columns). This function does not returns anything but it adds the :math:`n_1 \times n_2` covariance matrix to the kernpart to the object ``target`` (a :math:`n_1 \times n_2` np.array). This trick allows to compute the covariance matrix of a kernel containing many kernparts with a limited memory use. ::
+
+    def K(self,X,X2,target):
+        if X2 is None: X2 = X
+        dist2 = np.square((X-X2.T)/self.lengthscale)
+        target += self.variance*(1 + dist2/2.)**(-self.power)
+
+**Kdiag(self,X,target)**
+
+The implementation of this function in mandatory.
+
+This function is similar to ``K`` but it computes only the values of the kernel on the diagonal. Thus, ``target`` is a 1-dimensional np.array of length :math:`n_1`. ::
+
+    def Kdiag(self,X,target):
+        target += self.variance    
+
+
+**dK_dtheta(self,dL_dK,X,X2,target)**
+
+This function is required for the optimization of the parameters.
+
+Computes the derivative of the likelihood. As previously, the values are added to the object target which is a 1-dimensional np.array of length ``self.Nparam``. For example, if the kernel is parameterized by :math:`\sigma^2,\ \theta`, then :math:`\frac{dL}{d\sigma^2} = \frac{dL}{d K} \frac{dK}{d\sigma^2}` is added to the first element of target and :math:`\frac{dL}{d\theta} = \frac{dL}{d K} \frac{dK}{d\theta}` to the second. ::
+
+    def dK_dtheta(self,dL_dK,X,X2,target):
+        if X2 is None: X2 = X
+        dist2 = np.square((X-X2.T)/self.lengthscale)
+
+        dvar = (1 + dist2/2.)**(-self.power)
+        dl = self.power * self.variance * dist2 * self.lengthscale**(-3) * (1 + dist2/2./self.power)**(-self.power-1)
+        dp = - self.variance * np.log(1 + dist2/2.) * (1 + dist2/2.)**(-self.power)
+
+        target[0] += np.sum(dvar*dL_dK)
+        target[1] += np.sum(dl*dL_dK)
+        target[2] += np.sum(dp*dL_dK)
+
+
+**dKdiag_dtheta(self,dL_dKdiag,X,target)**
+
+This function is required for BGPLVM, sparse models and uncertain inputs.
+
+As previously, target is an ``self.Nparam`` array and :math:`\frac{dL}{d Kdiag} \frac{dKdiag}{dparam}` is added to each element. ::
+
+    def dKdiag_dtheta(self,dL_dKdiag,X,target):
+        target[0] += np.sum(dL_dKdiag)
+        # here self.lengthscale and self.power have no influence on Kdiag so target[1:] are unchanged
+
+**dK_dX(self,dL_dK,X,X2,target)**
+
+This function is required for GPLVM, BGPLVM, sparse models and uncertain inputs.
+
+Computes the derivative of the likelihood with respect to the inputs ``X`` (a :math:`n \times D` np.array). The result is added to target which is a :math:`n \times D` np.array. ::
+
+    def dK_dX(self,dL_dK,X,X2,target):
+        """derivative of the covariance matrix with respect to X."""
+        if X2 is None: X2 = X
+        dist2 = np.square((X-X2.T)/self.lengthscale)
+
+        dX = -self.variance*self.power * (X-X2.T)/self.lengthscale**2 *  (1 + dist2/2./self.power)**(-self.power-1)
+        target += np.sum(dL_dK*dX)
+
+**dKdiag_dX(self,dL_dKdiag,X,target)**
+
+This function is required for BGPLVM, sparse models and uncertain inputs. As for ``dKdiag_dtheta``, :math:`\frac{dL}{d Kdiag} \frac{dKdiag}{dX}` is added to each element of target. ::
+
+    def dKdiag_dX(self,dL_dKdiag,X,target):
+        pass
+
+**Psi statistics**
+
+The psi statistics and their derivatives are required for BGPLVM and GPS with uncertain inputs.
+
+The expressions of the psi statistics are:
+
+TODO
+
+For the rational quadratic we have:
+
+TODO
+
+Update the constructor 
+======================
+
+Once the required functions have been implemented as a kernpart object, the file GPy/kern/constructors.py has to be updated to allow to build a kernel based on the kernpart object.
+
+The following line should be added in the preamble of the file::
+
+    from rational_quadratic import rational_quadratic as rational_quadratic_part
+
+as well as the following block ::
+
+    def rational_quadratic(D,variance=1., lengthscale=1., power=1.):
+        part = rational_quadraticpart(D,variance, lengthscale, power)
+        return kern(D, [part])
+
+
+Update initialization
+=====================
+
+The last step is to update the list of kernels imported from constructor in GPy/kern/__init__.py.
+
+
+

doc/tuto_kernel_overview.rst

@@ -133,7 +133,7 @@ Various constrains can be applied to the parameters of a kernel
     * ``constrain_fixed`` to fix the value of a parameter (the value will not be modified during optimisation)
     * ``constrain_positive`` to make sure the parameter is greater than 0.
     * ``constrain_bounded`` to impose the parameter to be in a given range.
-    * ``tie_param`` to impose the value of two (or more) parameters to be equal.
+    * ``tie_params`` to impose the value of two (or more) parameters to be equal.
 
 When calling one of these functions, the parameters to constrain can either by specified by a regular expression that matches its name or by a number that corresponds to the rank of the parameter. Here is an example ::
 
@@ -146,7 +146,7 @@ When calling one of these functions, the parameters to constrain can either by s
 
     k.constrain_positive('var')
     k.constrain_fixed(np.array([1]),1.75)
-    k.tie_param('len')
+    k.tie_params('len')
     k.unconstrain('white')
     k.constrain_bounded('white',lower=1e-5,upper=.5)
     print k