merged master

GPy/__init__.py

@@ -6,5 +6,6 @@ import kern
 import models
 import inference
 import util
+import examples
 #import examples TODO: discuss!
 from core import priors

GPy/core/model.py

@@ -14,18 +14,18 @@ from ..inference import optimization
 class model(parameterised):
     def __init__(self):
         parameterised.__init__(self)
-        self.priors = [None for i in range(self.get_param().size)]
+        self.priors = [None for i in range(self._get_params().size)]
         self.optimization_runs = []
         self.sampling_runs = []
-        self.set_param(self.get_param())
+        self._set_params(self._get_params())
         self.preferred_optimizer = 'tnc'
-    def get_param(self):
+    def _get_params(self):
         raise NotImplementedError, "this needs to be implemented to utilise the model class"
-    def set_param(self,x):
+    def _set_params(self,x):
         raise NotImplementedError, "this needs to be implemented to utilise the model class"
     def log_likelihood(self):
         raise NotImplementedError, "this needs to be implemented to utilise the model class"
-    def log_likelihood_gradients(self):
+    def _log_likelihood_gradients(self):
         raise NotImplementedError, "this needs to be implemented to utilise the model class"
 
     def set_prior(self,which,what):
@@ -67,7 +67,7 @@ class model(parameterised):
             unconst = np.setdiff1d(which, self.constrained_positive_indices)
             if len(unconst):
                 print "Warning: constraining parameters to be positive:"
-                print '\n'.join([n for i,n in enumerate(self.get_param_names()) if i in unconst])
+                print '\n'.join([n for i,n in enumerate(self._get_param_names()) if i in unconst])
                 print '\n'
                 self.constrain_positive(unconst)
         elif isinstance(what,priors.Gaussian):
@@ -80,48 +80,65 @@ class model(parameterised):
         for w in which:
             self.priors[w] = what
 
-    def get(self,name):
+    def get(self,name, return_names=False):
         """
-        get a model parameter by name
+        Get a model parameter by name. The name is applied as a regular expression and all parameters that match that regular expression are returned. 
         """
         matches = self.grep_param_names(name)
         if len(matches):
-            return self.get_param()[matches]
+            if return_names:
+                return self._get_params()[matches], np.asarray(self._get_param_names())[matches].tolist()
+            else:
+                return self._get_params()[matches]
         else:
             raise AttributeError, "no parameter matches %s"%name
 
     def set(self,name,val):
         """
-        Set a model parameter by name
+        Set model parameter(s) by name. The name is provided as a regular expression. All parameters matching that regular expression are set to ghe given value.
         """
         matches = self.grep_param_names(name)
         if len(matches):
-            x = self.get_param()
+            x = self._get_params()
             x[matches] = val
-            self.set_param(x)
+            self._set_params(x)
+        else:
+            raise AttributeError, "no parameter matches %s"%name
+
+    def get_gradient(self,name, return_names=False):
+        """
+        Get model gradient(s) by name. The name is applied as a regular expression and all parameters that match that regular expression are returned. 
+        """
+        matches = self.grep_param_names(name)
+        if len(matches):
+            if return_names:
+                return self._log_likelihood_gradients()[matches],  np.asarray(self._get_param_names())[matches].tolist()
+            else:
+                return self._log_likelihood_gradients()[matches]
         else:
             raise AttributeError, "no parameter matches %s"%name
 
 
 
+
     def log_prior(self):
         """evaluate the prior"""
-        return np.sum([p.lnpdf(x) for p, x in zip(self.priors,self.get_param()) if p is not None])
+        return np.sum([p.lnpdf(x) for p, x in zip(self.priors,self._get_params()) if p is not None])
 
-    def log_prior_gradients(self):
+    def _log_prior_gradients(self):
         """evaluate the gradients of the priors"""
-        x = self.get_param()
+        x = self._get_params()
         ret = np.zeros(x.size)
         [np.put(ret,i,p.lnpdf_grad(xx)) for i,(p,xx) in enumerate(zip(self.priors,x)) if not p is None]
         return ret
 
-    def extract_gradients(self):
+    def _log_likelihood_gradients_transformed(self):
         """
         Use self.log_likelihood_gradients and self.prior_gradients to get the gradients of the model.
         Adjust the gradient for constraints and ties, return.
         """
-        g = self.log_likelihood_gradients() + self.log_prior_gradients()
-        x = self.get_param()
+        g = self._log_likelihood_gradients() + self._log_prior_gradients()
+        x = self._get_params()
         g[self.constrained_positive_indices] = g[self.constrained_positive_indices]*x[self.constrained_positive_indices]
         g[self.constrained_negative_indices] = g[self.constrained_negative_indices]*x[self.constrained_negative_indices]
         [np.put(g,i,g[i]*(x[i]-l)*(h-x[i])/(h-l)) for i,l,h in zip(self.constrained_bounded_indices, self.constrained_bounded_lowers, self.constrained_bounded_uppers)]
@@ -138,14 +155,14 @@ class model(parameterised):
         Make this draw from the prior if one exists, else draw from N(0,1)
         """
         #first take care of all parameters (from N(0,1))
-        x = self.extract_param()
+        x = self._get_params_transformed()
         x = np.random.randn(x.size)
-        self.expand_param(x)
+        self._set_params_transformed(x)
         #now draw from prior where possible
-        x = self.get_param()
+        x = self._get_params()
         [np.put(x,i,p.rvs(1)) for i,p in enumerate(self.priors) if not p is None]
-        self.set_param(x)
-        self.expand_param(self.extract_param())#makes sure all of the tied parameters get the same init (since there's only one prior object...)
+        self._set_params(x)
+        self._set_params_transformed(self._get_params_transformed())#makes sure all of the tied parameters get the same init (since there's only one prior object...)
 
 
     def optimize_restarts(self, Nrestarts=10, robust=False, verbose=True, **kwargs):
@@ -165,7 +182,7 @@ class model(parameterised):
         :verbose: whether to show informations about the current restart
         """
 
-        initial_parameters = self.extract_param()
+        initial_parameters = self._get_params_transformed()
         for i in range(Nrestarts):
             try:
                 self.randomize()
@@ -181,10 +198,10 @@ class model(parameterised):
                 else:
                     raise e
         if len(self.optimization_runs):
-            i = np.argmax([o.f_opt for o in self.optimization_runs])
-            self.expand_param(self.optimization_runs[i].x_opt)
+            i = np.argmin([o.f_opt for o in self.optimization_runs])
+            self._set_params_transformed(self.optimization_runs[i].x_opt)
         else:
-            self.expand_param(initial_parameters)
+            self._set_params_transformed(initial_parameters)
 
     def ensure_default_constraints(self,warn=False):
         """
@@ -194,7 +211,7 @@ class model(parameterised):
         for s in positive_strings:
             for i in self.grep_param_names(s):
                 if not (i in self.all_constrained_indices()):
-                    name = self.get_param_names()[i]
+                    name = self._get_param_names()[i]
                     self.constrain_positive(name)
                     if warn:
                         print "Warning! constraining %s postive"%name
@@ -214,24 +231,24 @@ class model(parameterised):
             optimizer = self.preferred_optimizer
 
         def f(x):
-            self.expand_param(x)
+            self._set_params_transformed(x)
             return -self.log_likelihood()-self.log_prior()
         def fp(x):
-            self.expand_param(x)
-            return -self.extract_gradients()
+            self._set_params_transformed(x)
+            return -self._log_likelihood_gradients_transformed()
         def f_fp(x):
-            self.expand_param(x)
-            return -self.log_likelihood()-self.log_prior(),-self.extract_gradients()
+            self._set_params_transformed(x)
+            return -self.log_likelihood()-self.log_prior(),-self._log_likelihood_gradients_transformed()
 
         if start == None:
-            start = self.extract_param()
+            start = self._get_params_transformed()
 
         optimizer = optimization.get_optimizer(optimizer)
         opt = optimizer(start, model = self, **kwargs)
         opt.run(f_fp=f_fp, f=f, fp=fp)
         self.optimization_runs.append(opt)
 
-        self.expand_param(opt.x_opt)
+        self._set_params_transformed(opt.x_opt)
 
     def optimize_SGD(self, momentum = 0.1, learning_rate = 0.01, iterations = 20, **kwargs):
         # assert self.Y.shape[1] > 1, "SGD only works with D > 1"
@@ -248,13 +265,13 @@ class model(parameterised):
 
         else:
             print "numerically calculating hessian. please be patient!"
-            x = self.get_param()
+            x = self._get_params()
             def f(x):
-                self.set_param(x)
+                self._set_params(x)
                 return self.log_likelihood()
             h = ndt.Hessian(f)
             A = -h(x)
-            self.set_param(x)
+            self._set_params(x)
         # check for almost zero components on the diagonal which screw up the cholesky
         aa = np.nonzero((np.diag(A)<1e-6) & (np.diag(A)>0.))[0]
         A[aa,aa] = 0.
@@ -268,7 +285,7 @@ class model(parameterised):
             hld = np.sum(np.log(np.diag(jitchol(A)[0])))
         except:
             return np.nan
-        return 0.5*self.get_param().size*np.log(2*np.pi) + self.log_likelihood() - hld
+        return 0.5*self._get_params().size*np.log(2*np.pi) + self.log_likelihood() - hld
 
     def __str__(self):
         s = parameterised.__str__(self).split('\n')
@@ -292,18 +309,18 @@ class model(parameterised):
         If the overall gradient fails, invividual components are tested.
         """
 
-        x = self.extract_param().copy()
+        x = self._get_params_transformed().copy()
 
         #choose a random direction to step in:
         dx = step*np.sign(np.random.uniform(-1,1,x.size))
 
         #evaulate around the point x
-        self.expand_param(x+dx)
-        f1,g1 = self.log_likelihood() + self.log_prior(), self.extract_gradients()
-        self.expand_param(x-dx)
-        f2,g2 = self.log_likelihood() + self.log_prior(), self.extract_gradients()
-        self.expand_param(x)
-        gradient = self.extract_gradients()
+        self._set_params_transformed(x+dx)
+        f1,g1 = self.log_likelihood() + self.log_prior(), self._log_likelihood_gradients_transformed()
+        self._set_params_transformed(x-dx)
+        f2,g2 = self.log_likelihood() + self.log_prior(), self._log_likelihood_gradients_transformed()
+        self._set_params_transformed(x)
+        gradient = self._log_likelihood_gradients_transformed()
 
         numerical_gradient = (f1-f2)/(2*dx)
         ratio = (f1-f2)/(2*np.dot(dx,gradient))
@@ -319,7 +336,7 @@ class model(parameterised):
                 print "Global check failed. Testing individual gradients\n"
 
                 try:
-                    names = self.extract_param_names()
+                    names = self._get_param_names_transformed()
                 except NotImplementedError:
                     names = ['Variable %i'%i for i in range(len(x))]
 
@@ -338,13 +355,13 @@ class model(parameterised):
             for i in range(len(x)):
                 xx = x.copy()
                 xx[i] += step
-                self.expand_param(xx)
-                f1,g1 = self.log_likelihood() + self.log_prior(), self.extract_gradients()[i]
+                self._set_params_transformed(xx)
+                f1,g1 = self.log_likelihood() + self.log_prior(), self._log_likelihood_gradients_transformed()[i]
                 xx[i] -= 2.*step
-                self.expand_param(xx)
-                f2,g2 = self.log_likelihood() + self.log_prior(), self.extract_gradients()[i]
-                self.expand_param(x)
-                gradient = self.extract_gradients()[i]
+                self._set_params_transformed(xx)
+                f2,g2 = self.log_likelihood() + self.log_prior(), self._log_likelihood_gradients_transformed()[i]
+                self._set_params_transformed(x)
+                gradient = self._log_likelihood_gradients_transformed()[i]
 
 
                 numerical_gradient = (f1-f2)/(2*step)

GPy/core/parameterised.py

@@ -66,7 +66,7 @@ class parameterised(object):
         if hasattr(self,'prior'):
             pass
 
-        self.expand_param(self.extract_param())# sets tied parameters to single value
+        self._set_params_transformed(self._get_params_transformed())# sets tied parameters to single value
 
     def untie_everything(self):
         """Unties all parameters by setting tied_indices to an empty list."""
@@ -87,7 +87,7 @@ class parameterised(object):
 
         Returns
         -------
-        the indices of self.get_param_names which match the regular expression.
+        the indices of self._get_param_names which match the regular expression.
 
         Notes
         -----
@@ -96,9 +96,9 @@ class parameterised(object):
 
         if type(expr) is str:
             expr = re.compile(expr)
-            return np.nonzero([expr.search(name) for name in self.get_param_names()])[0]
+            return np.nonzero([expr.search(name) for name in self._get_param_names()])[0]
         elif type(expr) is re._pattern_type:
-            return np.nonzero([expr.search(name) for name in self.get_param_names()])[0]
+            return np.nonzero([expr.search(name) for name in self._get_param_names()])[0]
         else:
             return expr
 
@@ -115,11 +115,11 @@ class parameterised(object):
         assert not np.any(matches[:,None]==self.all_constrained_indices()), "Some indices are already constrained"
         self.constrained_positive_indices = np.hstack((self.constrained_positive_indices, matches))
         #check to ensure constraint is in place
-        x = self.get_param()
+        x = self._get_params()
         for i,xx in enumerate(x):
             if (xx<0) & (i in matches):
                 x[i] = -xx
-        self.set_param(x)
+        self._set_params(x)
 
 
     def unconstrain(self,which):
@@ -163,11 +163,11 @@ class parameterised(object):
         assert not np.any(matches[:,None]==self.all_constrained_indices()), "Some indices are already constrained"
         self.constrained_negative_indices = np.hstack((self.constrained_negative_indices, matches))
         #check to ensure constraint is in place
-        x = self.get_param()
+        x = self._get_params()
         for i,xx in enumerate(x):
             if (xx>0.) and (i in matches):
                 x[i] = -xx
-        self.set_param(x)
+        self._set_params(x)
 
 
 
@@ -187,11 +187,11 @@ class parameterised(object):
         self.constrained_bounded_uppers.append(upper)
         self.constrained_bounded_lowers.append(lower)
         #check to ensure constraint is in place
-        x = self.get_param()
+        x = self._get_params()
         for i,xx in enumerate(x):
             if ((xx<=lower)|(xx>=upper)) & (i in matches):
                 x[i] = sigmoid(xx)*(upper-lower) + lower
-        self.set_param(x)
+        self._set_params(x)
 
 
     def constrain_fixed(self, which, value = None):
@@ -213,14 +213,14 @@ class parameterised(object):
         if value != None:
             self.constrained_fixed_values.append(value)
         else:
-            self.constrained_fixed_values.append(self.get_param()[self.constrained_fixed_indices[-1]])
+            self.constrained_fixed_values.append(self._get_params()[self.constrained_fixed_indices[-1]])
 
         #self.constrained_fixed_values.append(value)
-        self.expand_param(self.extract_param())
+        self._set_params_transformed(self._get_params_transformed())
 
-    def extract_param(self):
-        """use self.get_param to get the 'true' parameters of the model, which are then tied, constrained and fixed"""
-        x = self.get_param()
+    def _get_params_transformed(self):
+        """use self._get_params to get the 'true' parameters of the model, which are then tied, constrained and fixed"""
+        x = self._get_params()
         x[self.constrained_positive_indices] = np.log(x[self.constrained_positive_indices])
         x[self.constrained_negative_indices] = np.log(-x[self.constrained_negative_indices])
         [np.put(x,i,np.log(np.clip(x[i]-l,1e-10,np.inf)/np.clip(h-x[i],1e-10,np.inf))) for i,l,h in zip(self.constrained_bounded_indices, self.constrained_bounded_lowers, self.constrained_bounded_uppers)]
@@ -232,8 +232,8 @@ class parameterised(object):
             return x
 
 
-    def expand_param(self,x):
-        """ takes the vector x, which is then modified (by untying, reparameterising or inserting fixed values), and then call self.set_param"""
+    def _set_params_transformed(self,x):
+        """ takes the vector x, which is then modified (by untying, reparameterising or inserting fixed values), and then call self._set_params"""
 
         #work out how many places are fixed, and where they are. tricky logic!
         Nfix_places = 0.
@@ -257,14 +257,14 @@ class parameterised(object):
         xx[self.constrained_positive_indices] = np.exp(xx[self.constrained_positive_indices])
         xx[self.constrained_negative_indices] = -np.exp(xx[self.constrained_negative_indices])
         [np.put(xx,i,low+sigmoid(xx[i])*(high-low)) for i,low,high in zip(self.constrained_bounded_indices, self.constrained_bounded_lowers, self.constrained_bounded_uppers)]
-        self.set_param(xx)
+        self._set_params(xx)
 
-    def extract_param_names(self):
+    def _get_param_names_transformed(self):
         """
         Returns the parameter names as propagated after constraining,
-        tying or fixing, i.e. a list of the same length as extract_param()
+        tying or fixing, i.e. a list of the same length as _get_params_transformed()
         """
-        n =  self.get_param_names()
+        n =  self._get_param_names()
 
         #remove/concatenate the tied parameter names
         if len(self.tied_indices):
@@ -294,13 +294,13 @@ class parameterised(object):
         """
         Return a string describing the parameter names and their ties and constraints
         """
-        names = self.get_param_names()
+        names = self._get_param_names()
         N = len(names)
 
         if not N:
             return "This object has no free parameters."
         header = ['Name','Value','Constraints','Ties']
-        values = self.get_param() #map(str,self.get_param())
+        values = self._get_params() #map(str,self._get_params())
         #sort out the constraints
         constraints = ['']*len(names)
         for i in self.constrained_positive_indices:

GPy/examples/GPLVM_demo.py

@@ -1,28 +0,0 @@
-# Copyright (c) 2012, GPy authors (see AUTHORS.txt).
-# Licensed under the BSD 3-clause license (see LICENSE.txt)
-
-
-import numpy as np
-import pylab as pb
-import GPy
-np.random.seed(1)
-print "GPLVM with RBF kernel"
-
-N = 100
-Q = 1
-D = 2
-X = np.random.rand(N, Q)
-k = GPy.kern.rbf(Q, 1.0, 2.0) + GPy.kern.white(Q, 0.00001)
-K = k.K(X)
-Y = np.random.multivariate_normal(np.zeros(N),K,D).T
-
-m = GPy.models.GPLVM(Y, Q)
-m.constrain_positive('(rbf|bias|white)')
-
-pb.figure()
-m.plot()
-pb.title('PCA initialisation')
-pb.figure()
-m.optimize(messages = 1)
-m.plot()
-pb.title('After optimisation')

GPy/examples/GP_regression_demo.py

@@ -1,51 +0,0 @@
-# Copyright (c) 2012, GPy authors (see AUTHORS.txt).
-# Licensed under the BSD 3-clause license (see LICENSE.txt)
-
-
-"""
-Simple Gaussian Processes regression with an RBF kernel
-"""
-import pylab as pb
-import numpy as np
-import GPy
-pb.ion()
-pb.close('all')
-
-
-######################################
-## 1 dimensional example
-
-# sample inputs and outputs
-X = np.random.uniform(-3.,3.,(20,1))
-Y = np.sin(X)+np.random.randn(20,1)*0.05
-
-# create simple GP model
-m = GPy.models.GP_regression(X,Y)
-
-# contrain all parameters to be positive
-m.constrain_positive('')
-
-# optimize and plot
-m.optimize('tnc', max_f_eval = 1000)
-m.plot()
-print(m)
-
-######################################
-## 2 dimensional example
-
-# sample inputs and outputs
-X = np.random.uniform(-3.,3.,(40,2))
-Y = np.sin(X[:,0:1]) * np.sin(X[:,1:2])+np.random.randn(40,1)*0.05
-
-# create simple GP model
-m = GPy.models.GP_regression(X,Y)
-
-# contrain all parameters to be positive
-m.constrain_positive('')
-# optimize and plot
-pb.figure()
-m.optimize('tnc', max_f_eval = 1000)
-m.plot()
-print(m)
-
-

GPy/examples/GP_regression_kern_demo.py

@@ -1,33 +0,0 @@
-# Copyright (c) 2012, GPy authors (see AUTHORS.txt).
-# Licensed under the BSD 3-clause license (see LICENSE.txt)
-
-
-"""
-Simple one-dimensional Gaussian Processes with assorted kernel functions
-"""
-import pylab as pb
-import numpy as np
-import GPy
-
-# sample inputs and outputs
-D = 1
-X = np.random.randn(10,D)*2
-X = np.linspace(-1.5,1.5,5)[:,None]
-X = np.append(X,[[5]],0)
-Y = np.sin(np.pi*X/2) #+np.random.randn(X.shape[0],1)*0.05
-
-models = [GPy.models.GP_regression(X,Y, k) for k in (GPy.kern.rbf(D), GPy.kern.Matern52(D), GPy.kern.Matern32(D), GPy.kern.exponential(D), GPy.kern.linear(D) + GPy.kern.white(D),  GPy.kern.bias(D) + GPy.kern.white(D))]
-
-pb.figure(figsize=(12,8))
-for i,m in enumerate(models):
-    m.constrain_positive('')
-    m.optimize()
-    pb.subplot(3,2,i+1)
-    m.plot()
-    #pb.title(m.kern.parts[0].name)
-
-GPy.util.plot.align_subplots(3,2,(-3,6),(-2.5,2.5))
-
-pb.show()
-
-

GPy/examples/classification.py

@@ -8,8 +8,6 @@ Simple Gaussian Processes classification
 import pylab as pb
 import numpy as np
 import GPy
-pb.ion()
-pb.close('all')
 
 default_seed=10000
 ######################################
@@ -27,7 +25,7 @@ def crescent_data(model_type='Full', inducing=10, seed=default_seed):
     likelihood = GPy.inference.likelihoods.probit(data['Y'])
 
     if model_type=='Full':
-        m = GPy.models.simple_GP_EP(data['X'],likelihood)
+        m = GPy.models.GP_EP(data['X'],likelihood)
     else:
         # create sparse GP EP model
         m = GPy.models.sparse_GP_EP(data['X'],likelihood=likelihood,inducing=inducing,ep_proxy=model_type)
@@ -49,7 +47,7 @@ def oil():
     likelihood = GPy.inference.likelihoods.probit(data['Y'][:, 0:1])
 
     # create simple GP model
-    m = GPy.models.simple_GP_EP(data['X'],likelihood)
+    m = GPy.models.GP_EP(data['X'],likelihood)
 
     # contrain all parameters to be positive
     m.constrain_positive('')

GPy/examples/oil_flow_demo.py

@@ -1,53 +0,0 @@
-# Copyright (c) 2012, GPy authors (see AUTHORS.txt).
-# Licensed under the BSD 3-clause license (see LICENSE.txt)
-
-
-import cPickle as pickle
-import numpy as np
-import pylab as pb
-import GPy
-import pylab as plt
-np.random.seed(1)
-
-def plot_oil(X, theta, labels, label):
-    plt.figure()
-    X = X[:,np.argsort(theta)[:2]]
-    flow_type = (X[labels[:,0]==1])
-    plt.plot(flow_type[:,0], flow_type[:,1], 'rx')
-    flow_type = (X[labels[:,1]==1])
-    plt.plot(flow_type[:,0], flow_type[:,1], 'gx')
-    flow_type = (X[labels[:,2]==1])
-    plt.plot(flow_type[:,0], flow_type[:,1], 'bx')
-    plt.title(label)
-
-data = pickle.load(open('../util/datasets/oil_flow_3classes.pickle', 'r'))
-
-Y = data['DataTrn']
-N, D = Y.shape
-selected = np.random.permutation(N)[:200]
-labels = data['DataTrnLbls'][selected]
-Y = Y[selected]
-N, D = Y.shape
-Y -= Y.mean(axis=0)
-Y /= Y.std(axis=0)
-
-Q = 2
-m1 = GPy.models.sparse_GPLVM(Y, Q, M = 15)
-m1.constrain_positive('(rbf|bias|noise)')
-m1.constrain_bounded('white', 1e-6, 1.0)
-
-plot_oil(m1.X, np.array([1,1]), labels, 'PCA initialization')
-# m.optimize(messages = True)
-m1.optimize('bfgs', messages = True)
-plot_oil(m1.X, np.array([1,1]), labels, 'sparse GPLVM')
-# pb.figure()
-# m.plot()
-# pb.title('PCA initialisation')
-# pb.figure()
-# m.optimize(messages = 1)
-# m.plot()
-# pb.title('After optimisation')
-m = GPy.models.GPLVM(Y, Q)
-m.constrain_positive('(white|rbf|bias|noise)')
-m.optimize()
-plot_oil(m.X, np.array([1,1]), labels, 'GPLVM')

GPy/examples/regression.py

@@ -8,8 +8,6 @@ Gaussian Processes regression examples
 import pylab as pb
 import numpy as np
 import GPy
-pb.ion()
-pb.close('all')
 
 
 def toy_rbf_1d():
@@ -19,10 +17,8 @@ def toy_rbf_1d():
     # create simple GP model
     m = GPy.models.GP_regression(data['X'],data['Y'])
 
-    # contrain all parameters to be positive
-    m.constrain_positive('')
-
     # optimize
+    m.ensure_default_constraints()
     m.optimize()
 
     # plot
@@ -37,10 +33,8 @@ def rogers_girolami_olympics():
     # create simple GP model
     m = GPy.models.GP_regression(data['X'],data['Y'])
 
-    # contrain all parameters to be positive
-    m.constrain_positive('')
-
     # optimize
+    m.ensure_default_constraints()
     m.optimize()
 
     # plot
@@ -48,6 +42,10 @@ def rogers_girolami_olympics():
     print(m)
     return m
 
+def della_gatta_TRP63_gene_expression(number=942):
+    """Run a standard Gaussian process regression on the della Gatta et al TRP63 Gene Expression data set for a given gene number."""
+
+
 def toy_rbf_1d_50():
     """Run a simple demonstration of a standard Gaussian process fitting it to data sampled from an RBF covariance."""
     data = GPy.util.datasets.toy_rbf_1d_50()
@@ -55,10 +53,8 @@ def toy_rbf_1d_50():
     # create simple GP model
     m = GPy.models.GP_regression(data['X'],data['Y'])
 
-    # contrain all parameters to be positive
-    m.constrain_positive('')
-
     # optimize
+    m.ensure_default_constraints()
     m.optimize()
 
     # plot
@@ -73,11 +69,95 @@ def silhouette():
     # create simple GP model
     m = GPy.models.GP_regression(data['X'],data['Y'])
 
-    # contrain all parameters to be positive
-    m.constrain_positive('')
-
     # optimize
+    m.ensure_default_constraints()
     m.optimize()
 
     print(m)
     return m
+
+
+def multiple_optima(gene_number=937,resolution=80, model_restarts=10, seed=10000):
+    """Show an example of a multimodal error surface for Gaussian process regression. Gene 939 has bimodal behaviour where the noisey mode is higher."""
+
+    # Contour over a range of length scales and signal/noise ratios.
+    length_scales = np.linspace(0.1, 60., resolution)
+    log_SNRs = np.linspace(-3., 4., resolution)
+
+    data = GPy.util.datasets.della_gatta_TRP63_gene_expression(gene_number)
+    # Sub sample the data to ensure multiple optima
+    #data['Y'] = data['Y'][0::2, :]
+    #data['X'] = data['X'][0::2, :]
+
+    # 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)
+    pb.contour(length_scales, log_SNRs, np.exp(lls), 20)
+    ax = pb.gca()
+    pb.xlabel('length scale')
+    pb.ylabel('log_10 SNR')
+
+    xlim = ax.get_xlim()
+    ylim = ax.get_ylim()
+    
+    # Now run a few optimizations
+    models = []
+    optim_point_x = np.empty(2)
+    optim_point_y = np.empty(2)
+    np.random.seed(seed=seed)
+    for i in range(0, model_restarts):
+        kern = GPy.kern.rbf(1, variance=np.random.exponential(1.), lengthscale=np.random.exponential(50.)) + GPy.kern.white(1,variance=np.random.exponential(1.))
+        
+        m = GPy.models.GP_regression(data['X'],data['Y'], kernel=kern)
+        optim_point_x[0] = m.get('rbf_lengthscale')
+        optim_point_y[0] = np.log10(m.get('rbf_variance')) - np.log10(m.get('white_variance'));
+        
+        # optimize
+        m.ensure_default_constraints()
+        m.optimize(xtol=1e-6,ftol=1e-6)
+
+        optim_point_x[1] = m.get('rbf_lengthscale')
+        optim_point_y[1] = np.log10(m.get('rbf_variance')) - np.log10(m.get('white_variance'));
+        
+        pb.arrow(optim_point_x[0], optim_point_y[0], optim_point_x[1]-optim_point_x[0], optim_point_y[1]-optim_point_y[0], label=str(i), head_length=1, head_width=0.5, fc='k', ec='k')
+        models.append(m)
+
+    ax.set_xlim(xlim)
+    ax.set_ylim(ylim)
+    return (models, lls)
+
+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.
+    :length_scales: a list of length scales to explore for the contour plot.
+    :log_SNRs: a list of base 10 logarithm signal to noise ratios to explore for the contour plot.
+    :signal_kernel: a kernel to use for the 'signal' portion of the data."""
+
+    lls = []
+    total_var = np.var(data['Y'])
+    for log_SNR in log_SNRs:
+        SNR = 10**log_SNR
+        length_scale_lls = []
+        for length_scale in length_scales:
+            noise_var = 1.
+            signal_var = SNR
+            noise_var = noise_var/(noise_var + signal_var)*total_var
+            signal_var = signal_var/(noise_var + signal_var)*total_var
+
+            signal_kernel = signal_kernel_call(1, variance=signal_var, lengthscale=length_scale)
+            noise_kernel = GPy.kern.white(1, variance=noise_var)
+            kernel = signal_kernel + noise_kernel
+            K = kernel.K(data['X'])
+            total_var = (np.dot(np.dot(data['Y'].T,GPy.util.linalg.pdinv(K)[0]), data['Y'])/data['Y'].shape[0])[0,0]
+            noise_var *= total_var
+            signal_var *= total_var
+            
+            kernel = signal_kernel_call(1, variance=signal_var, lengthscale=length_scale) + GPy.kern.white(1, variance=noise_var)
+
+            model = GPy.models.GP_regression(data['X'], data['Y'], kernel=kernel)
+            model.constrain_positive('')
+            length_scale_lls.append(model.log_likelihood())
+        lls.append(length_scale_lls)
+    return np.array(lls)

GPy/examples/sparse_GPLVM_demo.py

@@ -10,11 +10,11 @@ print "sparse GPLVM with RBF kernel"
 
 N = 100
 M = 4
-Q = 1
+Q = 2
 D = 2
 #generate GPLVM-like data
 X = np.random.rand(N, Q)
-k = GPy.kern.rbf(Q, 1.0, 2.0) + GPy.kern.white(Q, 0.00001)
+k = GPy.kern.rbf(Q,1.,2*np.ones((1,))) + GPy.kern.white(Q, 0.00001)
 K = k.K(X)
 Y = np.random.multivariate_normal(np.zeros(N),K,D).T
 

GPy/examples/unsupervised.py

@@ -0,0 +1,25 @@
+"""
+Usupervised learning with Gaussian Processes.
+"""
+import pylab as pb
+import numpy as np
+import GPy
+
+
+######################################
+## Oil data subsampled to 100 points.
+def oil_100():
+    data = GPy.util.datasets.oil_100()
+
+    # create simple GP model
+    m = GPy.models.GPLVM(data['X'], 2)
+
+
+    # optimize
+    m.ensure_default_constraints()
+    m.optimize()
+
+    # plot
+    print(m)
+    return m
+

GPy/inference/likelihoods.py

@@ -99,6 +99,6 @@ class probit(likelihood):
     def predictive_mean(self,mu,variance):
         return stats.norm.cdf(mu/np.sqrt(1+variance))
 
-    def log_likelihood_gradients():
+    def _log_likelihood_gradients():
         raise NotImplementedError
 

GPy/inference/samplers.py

@@ -17,7 +17,7 @@ class Metropolis_Hastings:
     def __init__(self,model,cov=None):
         """Metropolis Hastings, with tunings according to Gelman et al. """
         self.model = model
-        current = self.model.extract_param()
+        current = self.model._get_params_transformed()
         self.D = current.size
         self.chains = []
         if cov is None:
@@ -32,19 +32,19 @@ class Metropolis_Hastings:
         if start is None:
             self.model.randomize()
         else:
-            self.model.expand_param(start)
+            self.model._set_params_transformed(start)
 
 
 
     def sample(self, Ntotal, Nburn, Nthin, tune=True, tune_throughout=False, tune_interval=400):
-        current = self.model.extract_param()
+        current = self.model._get_params_transformed()
         fcurrent = self.model.log_likelihood() + self.model.log_prior()
         accepted = np.zeros(Ntotal,dtype=np.bool)
         for it in range(Ntotal):
             print "sample %d of %d\r"%(it,Ntotal),
             sys.stdout.flush()
             prop = np.random.multivariate_normal(current, self.cov*self.scale*self.scale)
-            self.model.expand_param(prop)
+            self.model._set_params_transformed(prop)
             fprop = self.model.log_likelihood() + self.model.log_prior()
 
             if fprop>fcurrent:#sample accepted, going 'uphill'
@@ -73,12 +73,12 @@ class Metropolis_Hastings:
 
     def predict(self,function,args):
         """Make a prediction for the function, to which we will pass the additional arguments"""
-        param = self.model.get_param()
+        param = self.model._get_params()
         fs = []
         for p in self.chain:
-            self.model.set_param(p)
+            self.model._set_params(p)
             fs.append(function(*args))
-        self.model.set_param(param)# reset model to starting state
+        self.model._set_params(param)# reset model to starting state
         return fs
 
 

GPy/kern/Brownian.py

@@ -23,16 +23,16 @@ class Brownian(kernpart):
         assert self.D==1, "Brownian motion in 1D only"
         self.Nparam = 1.
         self.name = 'Brownian'
-        self.set_param(np.array([variance]).flatten())
+        self._set_params(np.array([variance]).flatten())
 
-    def get_param(self):
+    def _get_params(self):
         return self.variance
 
-    def set_param(self,x):
+    def _set_params(self,x):
         assert x.shape==(1,)
         self.variance = x
 
-    def get_param_names(self):
+    def _get_param_names(self):
         return ['variance']
 
     def K(self,X,X2,target):

GPy/kern/Matern32.py

@@ -20,43 +20,54 @@ class Matern32(kernpart):
     :type D: int
     :param variance: the variance :math:`\sigma^2`
     :type variance: float
-    :param lengthscale: the lengthscales :math:`\ell_i`
-    :type lengthscale: np.ndarray of size (D,)
+    :param lengthscale: the vector of lengthscale :math:`\ell_i`
+    :type lengthscale: np.ndarray of size (1,) or (D,) depending on ARD
+    :param ARD: Auto Relevance Determination. If equal to "False", the kernel is isotropic (ie. one single lengthscale parameter \ell), otherwise there is one lengthscale parameter per dimension.
+    :type ARD: Boolean
     :rtype: kernel object
 
     """
 
-    def __init__(self,D,variance=1.,lengthscales=None):
+    def __init__(self,D,variance=1.,lengthscale=None,ARD=False):
         self.D = D
-        if lengthscales is not None:
-            assert lengthscales.shape==(self.D,)
+        self.ARD = ARD
+        if ARD == False:
+            self.Nparam = 2
+            self.name = 'Mat32'
+            if lengthscale is not None:
+                assert lengthscale.shape == (1,)
+            else:
+                lengthscale = np.ones(1)
         else:
-            lengthscales = np.ones(self.D)
-        self.Nparam = self.D + 1
-        self.name = 'Mat32'
-        self.set_param(np.hstack((variance,lengthscales)))
+            self.Nparam = self.D + 1
+            self.name = 'Mat32_ARD'
+            if lengthscale is not None:
+                assert lengthscale.shape == (self.D,)
+            else:
+                lengthscale = np.ones(self.D)
+        self._set_params(np.hstack((variance,lengthscale)))
 
-    def get_param(self):
+    def _get_params(self):
         """return the value of the parameters."""
-        return np.hstack((self.variance,self.lengthscales))
+        return np.hstack((self.variance,self.lengthscale))
 
-    def set_param(self,x):
+    def _set_params(self,x):
         """set the value of the parameters."""
-        assert x.size==(self.D+1)
+        assert x.size == self.Nparam
         self.variance = x[0]
-        self.lengthscales = x[1:]
+        self.lengthscale = x[1:]
 
-    def get_param_names(self):
+    def _get_param_names(self):
         """return parameter names."""
-        if self.D==1:
+        if self.Nparam == 2:
             return ['variance','lengthscale']
         else:
-            return ['variance']+['lengthscale_%i'%i for i in range(self.lengthscales.size)]
+            return ['variance']+['lengthscale_%i'%i for i in range(self.lengthscale.size)]
 
     def K(self,X,X2,target):
         """Compute the covariance matrix between X and X2."""
         if X2 is None: X2 = X
-        dist = np.sqrt(np.sum(np.square((X[:,None,:]-X2[None,:,:])/self.lengthscales),-1))
+        dist = np.sqrt(np.sum(np.square((X[:,None,:]-X2[None,:,:])/self.lengthscale),-1))
         np.add(self.variance*(1+np.sqrt(3.)*dist)*np.exp(-np.sqrt(3.)*dist), target,target)
 
     def Kdiag(self,X,target):
@@ -66,26 +77,33 @@ class Matern32(kernpart):
     def dK_dtheta(self,partial,X,X2,target):
         """derivative of the covariance matrix with respect to the parameters."""
         if X2 is None: X2 = X
-        dist = np.sqrt(np.sum(np.square((X[:,None,:]-X2[None,:,:])/self.lengthscales),-1))
+        dist = np.sqrt(np.sum(np.square((X[:,None,:]-X2[None,:,:])/self.lengthscale),-1))
         dvar = (1+np.sqrt(3.)*dist)*np.exp(-np.sqrt(3.)*dist)
         invdist = 1./np.where(dist!=0.,dist,np.inf)
-        dist2M = np.square(X[:,None,:]-X2[None,:,:])/self.lengthscales**3
-        dl = (self.variance* 3 * dist * np.exp(-np.sqrt(3.)*dist))[:,:,np.newaxis] * dist2M*invdist[:,:,np.newaxis]
+        dist2M = np.square(X[:,None,:]-X2[None,:,:])/self.lengthscale**3
+        #dl = (self.variance* 3 * dist * np.exp(-np.sqrt(3.)*dist))[:,:,np.newaxis] * dist2M*invdist[:,:,np.newaxis]
         target[0] += np.sum(dvar*partial)
-        target[1:] += (dl*partial[:,:,None]).sum(0).sum(0)
+        if self.ARD == True:
+            dl = (self.variance* 3 * dist * np.exp(-np.sqrt(3.)*dist))[:,:,np.newaxis] * dist2M*invdist[:,:,np.newaxis]
+            #dl = self.variance*dvar[:,:,None]*dist2M*invdist[:,:,None]
+            target[1:] += (dl*partial[:,:,None]).sum(0).sum(0)
+        else:
+            dl = (self.variance* 3 * dist * np.exp(-np.sqrt(3.)*dist)) * dist2M.sum(-1)*invdist
+            #dl = self.variance*dvar*dist2M.sum(-1)*invdist
+            target[1] += np.sum(dl*partial)
 
     def dKdiag_dtheta(self,partial,X,target):
         """derivative of the diagonal of the covariance matrix with respect to the parameters."""
         target[0] += np.sum(partial)
 
-    def dK_dX(self,X,X2,target):
+    def dK_dX(self,partial,X,X2,target):
         """derivative of the covariance matrix with respect to X."""
         if X2 is None: X2 = X
-        dist = np.sqrt(np.sum(np.square((X[:,None,:]-X2[None,:,:])/self.lengthscales),-1))[:,:,None]
-        ddist_dX = (X[:,None,:]-X2[None,:,:])/self.lengthscales**2/np.where(dist!=0.,dist,np.inf)
-        dK_dX += -  np.transpose(3*self.variance*dist*np.exp(-np.sqrt(3)*dist)*ddist_dX,(1,0,2))
+        dist = np.sqrt(np.sum(np.square((X[:,None,:]-X2[None,:,:])/self.lengthscale),-1))[:,:,None]
+        ddist_dX = (X[:,None,:]-X2[None,:,:])/self.lengthscale**2/np.where(dist!=0.,dist,np.inf)
+        dK_dX = - np.transpose(3*self.variance*dist*np.exp(-np.sqrt(3)*dist)*ddist_dX,(1,0,2))
         target += np.sum(dK_dX*partial.T[:,:,None],0)
-        
+
     def dKdiag_dX(self,X,target):
         pass
 
@@ -104,7 +122,7 @@ class Matern32(kernpart):
         """
         assert self.D == 1
         def L(x,i):
-            return(3./self.lengthscales**2*F[i](x) + 2*np.sqrt(3)/self.lengthscales*F1[i](x) + F2[i](x))
+            return(3./self.lengthscale**2*F[i](x) + 2*np.sqrt(3)/self.lengthscale*F1[i](x) + F2[i](x))
         n = F.shape[0]
         G = np.zeros((n,n))
         for i in range(n):
@@ -114,5 +132,5 @@ class Matern32(kernpart):
         F1lower = np.array([f(lower) for f in F1])[:,None]
         #print "OLD \n", np.dot(F1lower,F1lower.T), "\n \n"
         #return(G)
-        return(self.lengthscales**3/(12.*np.sqrt(3)*self.variance) * G + 1./self.variance*np.dot(Flower,Flower.T) + self.lengthscales**2/(3.*self.variance)*np.dot(F1lower,F1lower.T))
+        return(self.lengthscale**3/(12.*np.sqrt(3)*self.variance) * G + 1./self.variance*np.dot(Flower,Flower.T) + self.lengthscale**2/(3.*self.variance)*np.dot(F1lower,F1lower.T))
 

GPy/kern/Matern52.py

@@ -19,43 +19,53 @@ class Matern52(kernpart):
     :type D: int
     :param variance: the variance :math:`\sigma^2`
     :type variance: float
-    :param lengthscale: the lengthscales :math:`\ell_i`
-    :type lengthscale: np.ndarray of size (D,)
+    :param lengthscale: the vector of lengthscale :math:`\ell_i`
+    :type lengthscale: np.ndarray of size (1,) or (D,) depending on ARD
+    :param ARD: Auto Relevance Determination. If equal to "False", the kernel is isotropic (ie. one single lengthscale parameter \ell), otherwise there is one lengthscale parameter per dimension.
+    :type ARD: Boolean
     :rtype: kernel object
 
     """
-    def __init__(self,D,variance=1.,lengthscales=None):
+    def __init__(self,D,variance=1.,lengthscale=None,ARD=False):
         self.D = D
-        if lengthscales is not None:
-            assert lengthscales.shape==(self.D,)
+        self.ARD = ARD
+        if ARD == False:
+            self.Nparam = 2
+            self.name = 'Mat32'
+            if lengthscale is not None:
+                assert lengthscale.shape == (1,)
+            else:
+                lengthscale = np.ones(1)
         else:
-            lengthscales = np.ones(self.D)
-        self.Nparam = self.D + 1
-        self.name = 'Mat52'
-        self.set_param(np.hstack((variance,lengthscales)))
-        
+            self.Nparam = self.D + 1
+            self.name = 'Mat32_ARD'
+            if lengthscale is not None:
+                assert lengthscale.shape == (self.D,)
+            else:
+                lengthscale = np.ones(self.D)
+        self._set_params(np.hstack((variance,lengthscale)))
 
-    def get_param(self):
+    def _get_params(self):
         """return the value of the parameters."""
-        return np.hstack((self.variance,self.lengthscales))
+        return np.hstack((self.variance,self.lengthscale))
 
-    def set_param(self,x):
+    def _set_params(self,x):
         """set the value of the parameters."""
-        assert x.size==(self.D+1)
+        assert x.size == self.Nparam
         self.variance = x[0]
-        self.lengthscales = x[1:]
+        self.lengthscale = x[1:]
 
-    def get_param_names(self):
+    def _get_param_names(self):
         """return parameter names."""
-        if self.D==1:
+        if self.Nparam == 2:
             return ['variance','lengthscale']
         else:
-            return ['variance']+['lengthscale_%i'%i for i in range(self.lengthscales.size)]
+            return ['variance']+['lengthscale_%i'%i for i in range(self.lengthscale.size)]
 
     def K(self,X,X2,target):
         """Compute the covariance matrix between X and X2."""
         if X2 is None: X2 = X
-        dist = np.sqrt(np.sum(np.square((X[:,None,:]-X2[None,:,:])/self.lengthscales),-1))
+        dist = np.sqrt(np.sum(np.square((X[:,None,:]-X2[None,:,:])/self.lengthscale),-1))
         np.add(self.variance*(1+np.sqrt(5.)*dist+5./3*dist**2)*np.exp(-np.sqrt(5.)*dist), target,target)
 
     def Kdiag(self,X,target):
@@ -65,24 +75,30 @@ class Matern52(kernpart):
     def dK_dtheta(self,partial,X,X2,target):
         """derivative of the covariance matrix with respect to the parameters."""
         if X2 is None: X2 = X
-        dist = np.sqrt(np.sum(np.square((X[:,None,:]-X2[None,:,:])/self.lengthscales),-1))
+        dist = np.sqrt(np.sum(np.square((X[:,None,:]-X2[None,:,:])/self.lengthscale),-1))
         invdist = 1./np.where(dist!=0.,dist,np.inf)
-        dist2M = np.square(X[:,None,:]-X2[None,:,:])/self.lengthscales**3
+        dist2M = np.square(X[:,None,:]-X2[None,:,:])/self.lengthscale**3
         dvar = (1+np.sqrt(5.)*dist+5./3*dist**2)*np.exp(-np.sqrt(5.)*dist)
-        dl = (self.variance * 5./3 * dist * (1 + np.sqrt(5.)*dist ) * np.exp(-np.sqrt(5.)*dist))[:,:,np.newaxis] * dist2M*invdist[:,:,np.newaxis]
         target[0] += np.sum(dvar*partial)
-        target[1:] += (dl*partial[:,:,None]).sum(0).sum(0)
+        if self.ARD:
+            dl = (self.variance * 5./3 * dist * (1 + np.sqrt(5.)*dist ) * np.exp(-np.sqrt(5.)*dist))[:,:,np.newaxis] * dist2M*invdist[:,:,np.newaxis]
+            #dl = (self.variance* 3 * dist * np.exp(-np.sqrt(3.)*dist))[:,:,np.newaxis] * dist2M*invdist[:,:,np.newaxis]
+            target[1:] += (dl*partial[:,:,None]).sum(0).sum(0)
+        else:
+            dl = (self.variance * 5./3 * dist * (1 + np.sqrt(5.)*dist ) * np.exp(-np.sqrt(5.)*dist)) * dist2M.sum(-1)*invdist
+            #dl = (self.variance* 3 * dist * np.exp(-np.sqrt(3.)*dist)) * dist2M.sum(-1)*invdist
+            target[1] += np.sum(dl*partial)
 
     def dKdiag_dtheta(self,X,target):
         """derivative of the diagonal of the covariance matrix with respect to the parameters."""
         target[0] += np.sum(partial)
 
-    def dK_dX(self,X,X2,target):
+    def dK_dX(self,partial,X,X2,target):
         """derivative of the covariance matrix with respect to X."""
         if X2 is None: X2 = X
-        dist = np.sqrt(np.sum(np.square((X[:,None,:]-X2[None,:,:])/self.lengthscales),-1))[:,:,None]
-        ddist_dX = (X[:,None,:]-X2[None,:,:])/self.lengthscales**2/np.where(dist!=0.,dist,np.inf)
-        dK_dX += -  np.transpose(self.variance*5./3*dist*(1+np.sqrt(5)*dist)*np.exp(-np.sqrt(5)*dist)*ddist_dX,(1,0,2))
+        dist = np.sqrt(np.sum(np.square((X[:,None,:]-X2[None,:,:])/self.lengthscale),-1))[:,:,None]
+        ddist_dX = (X[:,None,:]-X2[None,:,:])/self.lengthscale**2/np.where(dist!=0.,dist,np.inf)
+        dK_dX = -  np.transpose(self.variance*5./3*dist*(1+np.sqrt(5)*dist)*np.exp(-np.sqrt(5)*dist)*ddist_dX,(1,0,2))
         target += np.sum(dK_dX*partial.T[:,:,None],0)
 
     def dKdiag_dX(self,X,target):
@@ -97,26 +113,26 @@ class Matern52(kernpart):
         :param F1: vector of derivatives of F
         :type F1: np.array
         :param F2: vector of second derivatives of F
-        :type F2: np.array 
+        :type F2: np.array
         :param F3: vector of third derivatives of F
-        :type F3: np.array    
+        :type F3: np.array
         :param lower,upper: boundaries of the input domain
-        :type lower,upper: floats  
+        :type lower,upper: floats
         """
         assert self.D == 1
         def L(x,i):
-            return(5*np.sqrt(5)/self.lengthscales**3*F[i](x) + 15./self.lengthscales**2*F1[i](x)+ 3*np.sqrt(5)/self.lengthscales*F2[i](x) + F3[i](x))
+            return(5*np.sqrt(5)/self.lengthscale**3*F[i](x) + 15./self.lengthscale**2*F1[i](x)+ 3*np.sqrt(5)/self.lengthscale*F2[i](x) + F3[i](x))
         n = F.shape[0]
         G = np.zeros((n,n))
         for i in range(n):
             for j in range(i,n):
                 G[i,j] = G[j,i] = integrate.quad(lambda x : L(x,i)*L(x,j),lower,upper)[0]
-        G_coef = 3.*self.lengthscales**5/(400*np.sqrt(5))
+        G_coef = 3.*self.lengthscale**5/(400*np.sqrt(5))
         Flower = np.array([f(lower) for f in F])[:,None]
         F1lower = np.array([f(lower) for f in F1])[:,None]
         F2lower = np.array([f(lower) for f in F2])[:,None]
-        orig = 9./8*np.dot(Flower,Flower.T) + 9.*self.lengthscales**4/200*np.dot(F2lower,F2lower.T)
-        orig2 = 3./5*self.lengthscales**2 * ( np.dot(F1lower,F1lower.T) + 1./8*np.dot(Flower,F2lower.T) + 1./8*np.dot(F2lower,Flower.T))
+        orig = 9./8*np.dot(Flower,Flower.T) + 9.*self.lengthscale**4/200*np.dot(F2lower,F2lower.T)
+        orig2 = 3./5*self.lengthscale**2 * ( np.dot(F1lower,F1lower.T) + 1./8*np.dot(Flower,F2lower.T) + 1./8*np.dot(F2lower,Flower.T))
         return(1./self.variance* (G_coef*G + orig + orig2))
 
 

GPy/kern/__init__.py

@@ -2,5 +2,5 @@
 # Licensed under the BSD 3-clause license (see LICENSE.txt)
 
 
-from constructors import rbf, Matern32, Matern52, exponential, linear, white, bias, finite_dimensional, rbf_ARD, spline, Brownian, linear_ARD, rbf_sympy, sympykern
+from constructors import rbf, Matern32, Matern52, exponential, linear, white, bias, finite_dimensional, spline, Brownian, linear_ARD, rbf_sympy, sympykern
 from kern import kern

GPy/kern/bias.py

@@ -17,16 +17,16 @@ class bias(kernpart):
         self.D = D
         self.Nparam = 1
         self.name = 'bias'
-        self.set_param(np.array([variance]).flatten())
+        self._set_params(np.array([variance]).flatten())
 
-    def get_param(self):
+    def _get_params(self):
         return self.variance
 
-    def set_param(self,x):
+    def _set_params(self,x):
         assert x.shape==(1,)
         self.variance = x
 
-    def get_param_names(self):
+    def _get_param_names(self):
         return ['variance']
 
     def K(self,X,X2,target):

GPy/kern/constructors.py

@@ -6,7 +6,6 @@ import numpy as np
 from kern import kern
 
 from rbf import rbf as rbfpart
-from rbf_ARD import rbf_ARD as rbf_ARD_part
 from white import white as whitepart
 from linear import linear as linearpart
 from linear_ARD import linear_ARD as linear_ARD_part
@@ -22,7 +21,7 @@ from Brownian import Brownian as Brownianpart
 #using meta-classes to make the objects construct properly wthout them.
 
 
-def rbf(D,variance=1., lengthscale=1.):
+def rbf(D,variance=1., lengthscale=None,ARD=False):
     """
     Construct an RBF kernel
 
@@ -32,22 +31,10 @@ def rbf(D,variance=1., lengthscale=1.):
     :type variance: float
     :param lengthscale: the lengthscale of the kernel
     :type lengthscale: float
+    :param ARD: Auto Relevance Determination (one lengthscale per dimension)
+    :type ARD: Boolean
     """
-    part = rbfpart(D,variance,lengthscale)
-    return kern(D, [part])
-
-def rbf_ARD(D,variance=1., lengthscales=None):
-    """
-    Construct an RBF kernel with Automatic Relevance Determination (ARD)
-
-    :param D: dimensionality of the kernel, obligatory
-    :type D: int
-    :param variance: the variance of the kernel
-    :type variance: float
-    :param lengthscales: the lengthscales of the kernel
-    :type lengthscales: None|np.ndarray
-    """
-    part = rbf_ARD_part(D,variance,lengthscales)
+    part = rbfpart(D,variance,lengthscale,ARD)
     return kern(D, [part])
 
 def linear(D,lengthscales=None):
@@ -86,43 +73,52 @@ def white(D,variance=1.):
     part = whitepart(D,variance)
     return kern(D, [part])
 
-def exponential(D,variance=1., lengthscales=None):
+def exponential(D,variance=1., lengthscale=None, ARD=False):
     """
-     Construct a exponential kernel.
+    Construct an exponential kernel
 
-     Arguments
-     ---------
-     D (int), obligatory
-     variance (float)
-     lengthscales (np.ndarray)
+    :param D: dimensionality of the kernel, obligatory
+    :type D: int
+    :param variance: the variance of the kernel
+    :type variance: float
+    :param lengthscale: the lengthscale of the kernel
+    :type lengthscale: float
+    :param ARD: Auto Relevance Determination (one lengthscale per dimension)
+    :type ARD: Boolean
     """
-    part = exponentialpart(D,variance, lengthscales)
+    part = exponentialpart(D,variance, lengthscale, ARD)
     return kern(D, [part])
 
-def Matern32(D,variance=1., lengthscales=None):
+def Matern32(D,variance=1., lengthscale=None, ARD=False):
     """
      Construct a Matern 3/2 kernel.
 
-     Arguments
-     ---------
-     D (int), obligatory
-     variance (float)
-     lengthscales (np.ndarray)
+    :param D: dimensionality of the kernel, obligatory
+    :type D: int
+    :param variance: the variance of the kernel
+    :type variance: float
+    :param lengthscale: the lengthscale of the kernel
+    :type lengthscale: float
+    :param ARD: Auto Relevance Determination (one lengthscale per dimension)
+    :type ARD: Boolean
     """
-    part = Matern32part(D,variance, lengthscales)
+    part = Matern32part(D,variance, lengthscale, ARD)
     return kern(D, [part])
 
-def Matern52(D,variance=1., lengthscales=None):
+def Matern52(D,variance=1., lengthscale=None, ARD=False):
     """
      Construct a Matern 5/2 kernel.
 
-     Arguments
-     ---------
-     D (int), obligatory
-     variance (float)
-     lengthscales (np.ndarray)
+    :param D: dimensionality of the kernel, obligatory
+    :type D: int
+    :param variance: the variance of the kernel
+    :type variance: float
+    :param lengthscale: the lengthscale of the kernel
+    :type lengthscale: float
+    :param ARD: Auto Relevance Determination (one lengthscale per dimension)
+    :type ARD: Boolean
     """
-    part = Matern52part(D,variance, lengthscales)
+    part = Matern52part(D,variance, lengthscale, ARD)
     return kern(D, [part])
 
 def bias(D,variance=1.):

GPy/kern/exponential.py

@@ -19,42 +19,53 @@ class exponential(kernpart):
     :type D: int
     :param variance: the variance :math:`\sigma^2`
     :type variance: float
-    :param lengthscale: the lengthscales :math:`\ell_i`
-    :type lengthscale: np.ndarray of size (D,)
+    :param lengthscale: the vector of lengthscale :math:`\ell_i`
+    :type lengthscale: np.ndarray of size (1,) or (D,) depending on ARD
+    :param ARD: Auto Relevance Determination. If equal to "False", the kernel is isotropic (ie. one single lengthscale parameter \ell), otherwise there is one lengthscale parameter per dimension.
+    :type ARD: Boolean
     :rtype: kernel object
 
     """
-    def __init__(self,D,variance=1.,lengthscales=None):
+    def __init__(self,D,variance=1.,lengthscale=None,ARD=False):
         self.D = D
-        if lengthscales is not None:
-            assert lengthscales.shape==(self.D,)
+        self.ARD = ARD
+        if ARD == False:
+            self.Nparam = 2
+            self.name = 'exp'
+            if lengthscale is not None:
+                assert lengthscale.shape == (1,)
+            else:
+                lengthscale = np.ones(1)
         else:
-            lengthscales = np.ones(self.D)
-        self.Nparam = self.D + 1
-        self.name = 'exp'
-        self.set_param(np.hstack((variance,lengthscales)))
+            self.Nparam = self.D + 1
+            self.name = 'exp_ARD'
+            if lengthscale is not None:
+                assert lengthscale.shape == (self.D,)
+            else:
+                lengthscale = np.ones(self.D)
+        self._set_params(np.hstack((variance,lengthscale)))
 
-    def get_param(self):
+    def _get_params(self):
         """return the value of the parameters."""
-        return np.hstack((self.variance,self.lengthscales))
+        return np.hstack((self.variance,self.lengthscale))
 
-    def set_param(self,x):
+    def _set_params(self,x):
         """set the value of the parameters."""
-        assert x.size==(self.D+1)
+        assert x.size == self.Nparam
         self.variance = x[0]
-        self.lengthscales = x[1:]
+        self.lengthscale = x[1:]
 
-    def get_param_names(self):
+    def _get_param_names(self):
         """return parameter names."""
-        if self.D==1:
+        if self.Nparam == 2:
             return ['variance','lengthscale']
         else:
-            return ['variance']+['lengthscale_%i'%i for i in range(self.lengthscales.size)]
+            return ['variance']+['lengthscale_%i'%i for i in range(self.lengthscale.size)]
 
     def K(self,X,X2,target):
         """Compute the covariance matrix between X and X2."""
         if X2 is None: X2 = X
-        dist = np.sqrt(np.sum(np.square((X[:,None,:]-X2[None,:,:])/self.lengthscales),-1))
+        dist = np.sqrt(np.sum(np.square((X[:,None,:]-X2[None,:,:])/self.lengthscale),-1))
         np.add(self.variance*np.exp(-dist), target,target)
 
     def Kdiag(self,X,target):
@@ -64,24 +75,28 @@ class exponential(kernpart):
     def dK_dtheta(self,partial,X,X2,target):
         """derivative of the covariance matrix with respect to the parameters."""
         if X2 is None: X2 = X
-        dist = np.sqrt(np.sum(np.square((X[:,None,:]-X2[None,:,:])/self.lengthscales),-1))
+        dist = np.sqrt(np.sum(np.square((X[:,None,:]-X2[None,:,:])/self.lengthscale),-1))
         invdist = 1./np.where(dist!=0.,dist,np.inf)
-        dist2M = np.square(X[:,None,:]-X2[None,:,:])/self.lengthscales**3
+        dist2M = np.square(X[:,None,:]-X2[None,:,:])/self.lengthscale**3
         dvar = np.exp(-dist)
-        dl = self.variance*dvar[:,:,None]*dist2M*invdist[:,:,None]
         target[0] += np.sum(dvar*partial)
-        target[1:] += (dl*partial[:,:,None]).sum(0).sum(0)
+        if self.ARD == True:
+            dl = self.variance*dvar[:,:,None]*dist2M*invdist[:,:,None]
+            target[1:] += (dl*partial[:,:,None]).sum(0).sum(0)
+        else:
+            dl = self.variance*dvar*dist2M.sum(-1)*invdist
+            target[1] += np.sum(dl*partial)
 
     def dKdiag_dtheta(self,partial,X,target): 
         """derivative of the diagonal of the covariance matrix with respect to the parameters."""
         #NB: derivative of diagonal elements wrt lengthscale is 0
         target[0] += np.sum(partial)
 
-    def dK_dX(self,X,X2,target):
+    def dK_dX(self,partial,X,X2,target):
         """derivative of the covariance matrix with respect to X."""
         if X2 is None: X2 = X
-        dist = np.sqrt(np.sum(np.square((X[:,None,:]-X2[None,:,:])/self.lengthscales),-1))[:,:,None]
-        ddist_dX = (X[:,None,:]-X2[None,:,:])/self.lengthscales**2/np.where(dist!=0.,dist,np.inf)
+        dist = np.sqrt(np.sum(np.square((X[:,None,:]-X2[None,:,:])/self.lengthscale),-1))[:,:,None]
+        ddist_dX = (X[:,None,:]-X2[None,:,:])/self.lengthscale**2/np.where(dist!=0.,dist,np.inf)
         dK_dX = - np.transpose(self.variance*np.exp(-dist)*ddist_dX,(1,0,2))
         target += np.sum(dK_dX*partial.T[:,:,None],0)
 
@@ -101,14 +116,14 @@ class exponential(kernpart):
         """
         assert self.D == 1
         def L(x,i):
-            return(1./self.lengthscales*F[i](x) + F1[i](x))
+            return(1./self.lengthscale*F[i](x) + F1[i](x))
         n = F.shape[0]
         G = np.zeros((n,n))
         for i in range(n):
             for j in range(i,n):
                 G[i,j] = G[j,i] = integrate.quad(lambda x : L(x,i)*L(x,j),lower,upper)[0]
         Flower = np.array([f(lower) for f in F])[:,None]
-        return(self.lengthscales/2./self.variance * G + 1./self.variance * np.dot(Flower,Flower.T))
+        return(self.lengthscale/2./self.variance * G + 1./self.variance * np.dot(Flower,Flower.T))
 
 
 

GPy/kern/finite_dimensional.py

@@ -27,15 +27,15 @@ class finite_dimensional(kernpart):
             weights = np.ones(self.n)
         self.Nparam = self.n + 1
         self.name = 'finite_dim'
-        self.set_param(np.hstack((variance,weights)))
+        self._set_params(np.hstack((variance,weights)))
 
-    def get_param(self):
+    def _get_params(self):
         return np.hstack((self.variance,self.weights))
-    def set_param(self,x):
+    def _set_params(self,x):
         assert x.size == (self.Nparam)
         self.variance = x[0]
         self.weights = x[1:]
-    def get_param_names(self):
+    def _get_param_names(self):
         if self.n==1:
             return ['variance','weight']
         else:

GPy/kern/kern.py

@@ -133,20 +133,20 @@ class kern(parameterised):
         newkern.tied_indices = self.tied_indices + [self.Nparam + x for x in other.tied_indices]
         return newkern
 
-    def get_param(self):
-        return np.hstack([p.get_param() for p in self.parts])
+    def _get_params(self):
+        return np.hstack([p._get_params() for p in self.parts])
 
-    def set_param(self,x):
-        [p.set_param(x[s]) for p, s in zip(self.parts, self.param_slices)]
+    def _set_params(self,x):
+        [p._set_params(x[s]) for p, s in zip(self.parts, self.param_slices)]
 
-    def get_param_names(self):
+    def _get_param_names(self):
         #this is a bit nasty: we wat to distinguish between parts with the same name by appending a count
         part_names = np.array([k.name for k in self.parts],dtype=np.str)
         counts = [np.sum(part_names==ni) for i, ni in enumerate(part_names)]
         cum_counts = [np.sum(part_names[i:]==ni) for i, ni in enumerate(part_names)]
         names = [name+'_'+str(cum_count) if count>1 else name for name,count,cum_count in zip(part_names,counts,cum_counts)]
 
-        return sum([[name+'_'+n for n in k.get_param_names()] for name,k in zip(names,self.parts)],[])
+        return sum([[name+'_'+n for n in k._get_param_names()] for name,k in zip(names,self.parts)],[])
 
     def K(self,X,X2=None,slices1=None,slices2=None):
         assert X.shape[1]==self.D

GPy/kern/kernpart.py

@@ -16,11 +16,11 @@ class kernpart(object):
         self.Nparam = 1
         self.name = 'unnamed'
 
-    def get_param(self):
+    def _get_params(self):
         raise NotImplementedError
-    def set_param(self,x):
+    def _set_params(self,x):
         raise NotImplementedError
-    def get_param_names(self):
+    def _get_param_names(self):
         raise NotImplementedError
     def K(self,X,X2,target):
         raise NotImplementedError

GPy/kern/linear.py

@@ -20,16 +20,16 @@ class linear(kernpart):
             variance = 1.0
         self.Nparam = 1
         self.name = 'linear'
-        self.set_param(variance)
+        self._set_params(variance)
         self._Xcache, self._X2cache = np.empty(shape=(2,))
 
-    def get_param(self):
+    def _get_params(self):
         return self.variance
 
-    def set_param(self,x):
+    def _set_params(self,x):
         self.variance = x
 
-    def get_param_names(self):
+    def _get_param_names(self):
         return ['variance']
 
     def K(self,X,X2,target):

GPy/kern/linear_ARD.py

@@ -23,16 +23,16 @@ class linear_ARD(kernpart):
             variances = np.ones(self.D)
         self.Nparam = int(self.D)
         self.name = 'linear'
-        self.set_param(variances)
+        self._set_params(variances)
 
-    def get_param(self):
+    def _get_params(self):
         return self.variances
 
-    def set_param(self,x):
+    def _set_params(self,x):
         assert x.size==(self.Nparam)
         self.variances = x
 
-    def get_param_names(self):
+    def _get_param_names(self):
         if self.D==1:
             return ['variance']
         else:

GPy/kern/rbf.py

@@ -8,46 +8,67 @@ import hashlib
 
 class rbf(kernpart):
     """
-    Radial Basis Function kernel, aka squared-exponential, exponentiated quadratic or Gaussian kernel.
+    Radial Basis Function kernel, aka squared-exponential, exponentiated quadratic or Gaussian kernel:
 
     .. math::
 
-       k(r) = \sigma^2 \exp(- \frac{r^2}{2\ell}) \qquad \qquad \\text{ where  } r = \sqrt{\frac{\sum_{i=1}^d (x_i-x^\prime_i)^2}{\ell^2}}
+       k(r) = \sigma^2 \exp(- \frac{1}{2}r^2) \qquad \qquad \\text{ where  } r^2 = \sum_{i=1}^d \frac{ (x_i-x^\prime_i)^2}{\ell_i^2}}
 
-    where \ell is the lengthscale, \alpha the smoothness, \sigma^2 the variance and d the dimensionality of the input.
+    where \ell_i is the lengthscale, \sigma^2 the variance and d the dimensionality of the input.
 
     :param D: the number of input dimensions
     :type D: int
     :param variance: the variance of the kernel
     :type variance: float
-    :param lengthscale: the lengthscale of the kernel
-    :type lengthscale: float
+    :param lengthscale: the vector of lengthscale of the kernel
+    :type lengthscale: np.ndarray od size (1,) or (D,) depending on ARD
+    :param ARD: Auto Relevance Determination. If equal to "False", the kernel is isotropic (ie. one single lengthscale parameter \ell), otherwise there is one lengthscale parameter per dimension.
+    :type ARD: Boolean
+    :rtype: kernel object
 
-    .. Note: for rbf with different lengthscale on each dimension, see rbf_ARD
     """
 
-    def __init__(self,D,variance=1.,lengthscale=1.):
+    def __init__(self,D,variance=1.,lengthscale=None,ARD=False):
         self.D = D
-        self.Nparam = 2
-        self.name = 'rbf'
-        self.set_param(np.hstack((variance,lengthscale)))
+        self.ARD = ARD
+        if ARD == False:
+            self.Nparam = 2
+            self.name = 'rbf'
+            if lengthscale is not None:
+                assert lengthscale.shape == (1,)
+            else:
+                lengthscale = np.ones(1)     
+        else:
+            self.Nparam = self.D + 1
+            self.name = 'rbf_ARD'
+            if lengthscale is not None:
+                assert lengthscale.shape == (self.D,)
+            else:
+                lengthscale = np.ones(self.D)
+
+        self._set_params(np.hstack((variance,lengthscale)))        
 
         #initialize cache
         self._Z, self._mu, self._S = np.empty(shape=(3,1))
         self._X, self._X2, self._params = np.empty(shape=(3,1))
 
-    def get_param(self):
+    def _get_params(self):
         return np.hstack((self.variance,self.lengthscale))
 
-    def set_param(self,x):
-        self.variance, self.lengthscale = x
+    def _set_params(self,x):
+        assert x.size==(self.Nparam)
+        self.variance = x[0]
+        self.lengthscale = x[1:]
         self.lengthscale2 = np.square(self.lengthscale)
         #reset cached results
         self._X, self._X2, self._params = np.empty(shape=(3,1))
         self._Z, self._mu, self._S = np.empty(shape=(3,1)) # cached versions of Z,mu,S
 
-    def get_param_names(self):
-        return ['variance','lengthscale']
+    def _get_param_names(self):
+        if self.Nparam == 2:
+            return ['variance','lengthscale']
+        else:
+            return ['variance']+['lengthscale_%i'%i for i in range(self.lengthscale.size)]        
 
     def K(self,X,X2,target):
         if X2 is None:
@@ -61,7 +82,12 @@ class rbf(kernpart):
     def dK_dtheta(self,partial,X,X2,target):
         self._K_computations(X,X2)
         target[0] += np.sum(self._K_dvar*partial)
-        target[1] += np.sum(self._K_dvar*self.variance*self._K_dist2/self.lengthscale*partial)
+        if self.ARD == True:
+            dl = self._K_dvar[:,:,None]*self.variance*self._K_dist2/self.lengthscale
+            target[1:] += (dl*partial[:,:,None]).sum(0).sum(0)
+        else:
+            target[1] += np.sum(self._K_dvar*self.variance*(self._K_dist2.sum(-1))/self.lengthscale*partial)
+        #np.sum(self._K_dvar*self.variance*self._K_dist2/self.lengthscale*partial)
 
     def dKdiag_dtheta(self,partial,X,target):
         #NB: derivative of diagonal elements wrt lengthscale is 0
@@ -81,15 +107,12 @@ class rbf(kernpart):
             self._X = X
             self._X2 = X2
             if X2 is None: X2 = X
-            XXT = np.dot(X,X2.T)
-            if X is X2:
-                self._K_dist2 = (-2.*XXT + np.diag(XXT)[:,np.newaxis] + np.diag(XXT)[np.newaxis,:])/self.lengthscale2
-            else:
-                self._K_dist2 = (-2.*XXT + np.sum(np.square(X),1)[:,None] + np.sum(np.square(X2),1)[None,:])/self.lengthscale2
-            # TODO Remove comments if this is fine.
-            # Commented out by Neil as doesn't seem to be used elsewhere.
-            #self._K_exponent = -0.5*self._K_dist2
-            self._K_dvar = np.exp(-0.5*self._K_dist2)
+            self._K_dist = X[:,None,:]-X2[None,:,:] # this can be computationally heavy
+            self._params = np.empty(shape=(1,0))    #ensure the next section gets called
+        if not np.all(self._params == self._get_params()):
+            self._params == self._get_params()
+            self._K_dist2 = np.square(self._K_dist/self.lengthscale)
+            self._K_dvar = np.exp(-0.5*self._K_dist2.sum(-1))
 
     def psi0(self,Z,mu,S,target):
         target += self.variance
@@ -132,7 +155,7 @@ class rbf(kernpart):
         d_length = self._psi2[:,:,:,None]*(0.5*self._psi2_Zdist_sq*self._psi2_denom + 2.*self._psi2_mudist_sq + 2.*S[:,None,None,:]/self.lengthscale2)/(self.lengthscale*self._psi2_denom)
         d_length = d_length.sum(0)
         target[0] += np.sum(partial*d_var)
-        target[1] += np.sum(d_length*partial)
+        target[1:] += (d_length*partial[:,:,None]).sum(0).sum(0)
 
     def dpsi2_dZ(self,partial,Z,mu,S,target):
         """Returns shape N,M,M,Q"""
@@ -175,4 +198,3 @@ class rbf(kernpart):
             self._psi2 = np.square(self.variance)*np.exp(self._psi2_exponent) # N,M,M
 
             self._Z, self._mu, self._S = Z, mu,S
-

GPy/kern/rbf_ARD.py

@@ -1,251 +0,0 @@
-# Copyright (c) 2012, GPy authors (see AUTHORS.txt).
-# Licensed under the BSD 3-clause license (see LICENSE.txt)
-
-
-from kernpart import kernpart
-import numpy as np
-import hashlib
-
-class rbf_ARD(kernpart):
-    def __init__(self,D,variance=1.,lengthscales=None):
-        """
-        Arguments
-        ----------
-        D: int - the number of input dimensions
-        variance: float
-        lengthscales : np.ndarray of shape (D,)
-        """
-        self.D = D
-        if lengthscales is not None:
-            assert lengthscales.shape==(self.D,)
-        else:
-            lengthscales = np.ones(self.D)
-        self.Nparam = self.D + 1
-        self.name = 'rbf_ARD'
-        self.set_param(np.hstack((variance,lengthscales)))
-
-        #initialize cache
-        self._Z, self._mu, self._S = np.empty(shape=(3,1))
-        self._X, self._X2, self._params = np.empty(shape=(3,1))
-
-    def get_param(self):
-        return np.hstack((self.variance,self.lengthscales))
-
-    def set_param(self,x):
-        assert x.size==(self.D+1)
-        self.variance = x[0]
-        self.lengthscales = x[1:]
-        self.lengthscales2 = np.square(self.lengthscales)
-        #reset cached results
-        self._Z, self._mu, self._S = np.empty(shape=(3,1)) # cached versions of Z,mu,S
-
-    def get_param_names(self):
-        if self.D==1:
-            return ['variance','lengthscale']
-        else:
-            return ['variance']+['lengthscale_%i'%i for i in range(self.lengthscales.size)]
-
-    def K(self,X,X2,target):
-        self._K_computations(X,X2)
-        np.add(self.variance*self._K_dvar, target,target)
-
-    def Kdiag(self,X,target):
-        np.add(target,self.variance,target)
-
-    def dK_dtheta(self,partial,X,X2,target):
-        self._K_computations(X,X2)
-        dl = self._K_dvar[:,:,None]*self.variance*self._K_dist2/self.lengthscales
-        target[0] += np.sum(self._K_dvar*partial)
-        target[1:] += (dl*partial[:,:,None]).sum(0).sum(0)
-
-    def dKdiag_dtheta(self,X,target):
-        target[0] += np.sum(partial)
-
-    def dK_dX(self,partial,X,X2,target):
-        self._K_computations(X,X2)
-        dZ = self.variance*self._K_dvar[:,:,None]*self._K_dist/self.lengthscales2
-        dK_dX = -dZ.transpose(1,0,2)
-        target += np.sum(dK_dX*partial.T[:,:,None],0)
-
-    def dKdiag_dX(self,partial,X,target):
-        pass
-
-    def psi0(self,Z,mu,S,target):
-        target += self.variance
-
-    def dpsi0_dtheta(self,partial,Z,mu,S,target):
-        target[0] += 1.
-
-    def dpsi0_dmuS(self,Z,mu,S,target_mu,target_S):
-        pass
-
-    def psi1(self,Z,mu,S,target):
-        self._psi_computations(Z,mu,S)
-        np.add(target, self._psi1,target)
-
-    def dpsi1_dtheta(self,partial,Z,mu,S,target):
-        self._psi_computations(Z,mu,S)
-        denom_deriv = S[:,None,:]/(self.lengthscales**3+self.lengthscales*S[:,None,:])
-        d_length = self._psi1[:,:,None]*(self.lengthscales*np.square(self._psi1_dist/(self.lengthscales2+S[:,None,:])) + denom_deriv)
-        target[0] += np.sum(partial*self._psi1/self.variance)
-        target[1:] += (d_length*partial[:,:,None]).sum(0).sum(0)
-
-    def dpsi1_dZ(self,partial,Z,mu,S,target):
-        self._psi_computations(Z,mu,S)
-        np.add(target,-self._psi1[:,:,None]*self._psi1_dist/self.lengthscales2/self._psi1_denom,target)
-        target += np.sum(partial[:,:,None]*-self._psi1[:,:,None]*self._psi1_dist/self.lengthscales2/self._psi1_denom,0)
-
-    def dpsi1_dmuS(self,partial,Z,mu,S,target_mu,target_S):
-        """return shapes are N,M,Q"""
-        self._psi_computations(Z,mu,S)
-        tmp = self._psi1[:,:,None]/self.lengthscales2/self._psi1_denom
-        target_mu += np.sum(partial*tmp*self._psi1_dist,1)
-        target_S += np.sum(partial*0.5*tmp*(self._psi1_dist_sq-1),1)
-
-    def psi2(self,Z,mu,S,target):
-        self._psi_computations(Z,mu,S)
-        target += self._psi2.sum(0) #TODO: psi2 should be NxMxM (for het. noise)
-
-    def dpsi2_dtheta(self,Z,mu,S,target):
-        """Shape N,M,M,Ntheta"""
-        self._psi_computations(Z,mu,S)
-        d_var = np.sum(2.*self._psi2/self.variance,0)
-        d_length = self._psi2[:,:,:,None]*(0.5*self._psi2_Zdist_sq*self._psi2_denom + 2.*self._psi2_mudist_sq + 2.*S[:,None,None,:]/self.lengthscales2)/(self.lengthscales*self._psi2_denom)
-        d_length = d_length.sum(0)
-        target[0] += np.sum(partial*d_var)
-        target[1:] += (d_length*partial[:,:,None]).sum(0).sum(0)
-
-    def dpsi2_dZ(self,Z,mu,S,target):
-        """Returns shape N,M,M,Q"""
-        self._psi_computations(Z,mu,S)
-        dZ = self._psi2[:,:,:,None]/self.lengthscales2*(-0.5*self._psi2_Zdist + self._psi2_mudist/self._psi2_denom)
-        target += np.sum(partial[None,:,:,None]*dZ,0).sum(1)
-
-    def dpsi2_dmuS(self,Z,mu,S,target_mu,target_S):
-        """Think N,M,M,Q """
-        self._psi_computations(Z,mu,S)
-        tmp = self._psi2[:,:,:,None]/self.lengthscales2/self._psi2_denom
-        target_mu += (partial*-tmp*2.*self._psi2_mudist).sum(1).sum(1)
-        target_S += (partial*tmp*(2.*self._psi2_mudist_sq-1)).sum(1).sum(1)
-
-    def _K_computations(self,X,X2):
-        if not (np.all(X==self._X) and np.all(X2==self._X2)):
-            self._X = X
-            self._X2 = X2
-            if X2 is None: X2 = X
-            self._K_dist = X[:,None,:]-X2[None,:,:] # this can be computationally heavy
-            self._params = np.empty(shape=(1,0))#ensure the next section gets called
-        if not np.all(self._params == self.get_param()):
-            self._params == self.get_param()
-            self._K_dist2 = np.square(self._K_dist/self.lengthscales)
-            self._K_exponent = -0.5*self._K_dist2.sum(-1)
-            self._K_dvar = np.exp(-0.5*self._K_dist2.sum(-1))
-
-    def _psi_computations(self,Z,mu,S):
-        #here are the "statistics" for psi1 and psi2
-        if not np.all(Z==self._Z):
-            #Z has changed, compute Z specific stuff
-            self._psi2_Zhat = 0.5*(Z[:,None,:] +Z[None,:,:]) # M,M,Q
-            self._psi2_Zdist = Z[:,None,:]-Z[None,:,:] # M,M,Q
-            self._psi2_Zdist_sq = np.square(self._psi2_Zdist)/self.lengthscales2 # M,M,Q
-            self._Z = Z
-
-        if not (np.all(Z==self._Z) and np.all(mu==self._mu) and np.all(S==self._S)):
-            #something's changed. recompute EVERYTHING
-
-            #psi1
-            self._psi1_denom = S[:,None,:]/self.lengthscales2 + 1.
-            self._psi1_dist = Z[None,:,:]-mu[:,None,:]
-            self._psi1_dist_sq = np.square(self._psi1_dist)/self.lengthscales2/self._psi1_denom
-            self._psi1_exponent = -0.5*np.sum(self._psi1_dist_sq+np.log(self._psi1_denom),-1)
-            self._psi1 = self.variance*np.exp(self._psi1_exponent)
-
-            #psi2
-            self._psi2_denom = 2.*S[:,None,None,:]/self.lengthscales2+1. # N,M,M,Q
-            self._psi2_mudist = mu[:,None,None,:]-self._psi2_Zhat #N,M,M,Q
-            self._psi2_mudist_sq = np.square(self._psi2_mudist)/(self.lengthscales2*self._psi2_denom)
-            self._psi2_exponent = np.sum(-self._psi2_Zdist_sq/4. -self._psi2_mudist_sq -0.5*np.log(self._psi2_denom),-1) #N,M,M
-            self._psi2 = np.square(self.variance)*np.exp(self._psi2_exponent) # N,M,M
-
-            self._Z, self._mu, self._S = Z, mu,S
-
-
-if __name__=='__main__':
-    #run some simple tests on the kernel (TODO:move these to unititest)
-    #TODO: these are broken in this new structure!
-    N = 10
-    M = 5
-    Q = 3
-
-    Z = np.random.randn(M,Q)
-    mu = np.random.randn(N,Q)
-    S = np.random.rand(N,Q)
-
-    var = 2.5
-    lengthscales = np.ones(Q)*0.7
-
-    k = rbf(Q,var,lengthscales)
-
-    from checkgrad import checkgrad
-
-    def k_theta_test(param,k):
-        k.set_param(param)
-        K = k.K(Z)
-        dK_dtheta = k.dK_dtheta(Z)
-        f = np.sum(K)
-        df = dK_dtheta.sum(0).sum(0)
-        return f,np.array(df)
-    print "dk_dtheta_test"
-    checkgrad(k_theta_test,np.random.randn(1+Q),args=(k,))
-
-
-    def psi1_mu_test(mu,k):
-        mu = mu.reshape(N,Q)
-        f = np.sum(k.psi1(Z,mu,S))
-        df = k.dpsi1_dmuS(Z,mu,S)[0].sum(1)
-        return f,df.flatten()
-    print "psi1_mu_test"
-    checkgrad(psi1_mu_test,np.random.randn(N*Q),args=(k,))
-
-    def psi1_S_test(S,k):
-        S = S.reshape(N,Q)
-        f = np.sum(k.psi1(Z,mu,S))
-        df = k.dpsi1_dmuS(Z,mu,S)[1].sum(1)
-        return f,df.flatten()
-    print "psi1_S_test"
-    checkgrad(psi1_S_test,np.random.rand(N*Q),args=(k,))
-
-    def psi1_theta_test(theta,k):
-        k.set_param(theta)
-        f = np.sum(k.psi1(Z,mu,S))
-        df = np.array([np.sum(grad) for grad in k.dpsi1_dtheta(Z,mu,S)])
-        return f,df
-    print "psi1_theta_test"
-    checkgrad(psi1_theta_test,np.random.rand(1+Q),args=(k,))
-
-
-    def psi2_mu_test(mu,k):
-        mu = mu.reshape(N,Q)
-        f = np.sum(k.psi2(Z,mu,S))
-        df = k.dpsi2_dmuS(Z,mu,S)[0].sum(1).sum(1)
-        return f,df.flatten()
-    print "psi2_mu_test"
-    checkgrad(psi2_mu_test,np.random.randn(N*Q),args=(k,))
-
-    def psi2_S_test(S,k):
-        S = S.reshape(N,Q)
-        f = np.sum(k.psi2(Z,mu,S))
-        df = k.dpsi2_dmuS(Z,mu,S)[1].sum(1).sum(1)
-        return f,df.flatten()
-    print "psi2_S_test"
-    checkgrad(psi2_S_test,np.random.rand(N*Q),args=(k,))
-
-    def psi2_theta_test(theta,k):
-        k.set_param(theta)
-        f = np.sum(k.psi2(Z,mu,S))
-        df = np.array([np.sum(grad) for grad in k.dpsi2_dtheta(Z,mu,S)])
-        return f,df
-    print "psi2_theta_test"
-    checkgrad(psi2_theta_test,np.random.rand(1+Q),args=(k,))
-
-

GPy/kern/spline.py

@@ -25,15 +25,15 @@ class spline(kernpart):
         assert self.D==1
         self.Nparam = 1
         self.name = 'spline'
-        self.set_param(np.squeeze(variance))
+        self._set_params(np.squeeze(variance))
 
-    def get_param(self):
+    def _get_params(self):
         return self.variance
 
-    def set_param(self,x):
+    def _set_params(self,x):
         self.variance = x
 
-    def get_param_names(self):
+    def _get_param_names(self):
         return ['variance']
 
     def K(self,X,X2,target):

GPy/kern/sympykern.py

@@ -44,7 +44,7 @@ class spkern(kernpart):
         if param is None:
             param = np.ones(self.Nparam)
         assert param.size==self.Nparam
-        self.set_param(param)
+        self._set_params(param)
 
         #Differentiate!
         self._sp_dk_dtheta = [sp.diff(k,theta).simplify() for theta in self._sp_theta]
@@ -247,12 +247,12 @@ class spkern(kernpart):
         Z = X
         weave.inline(self._dKdiag_dX_code,arg_names=['target','X','Z','param','partial'],**self.weave_kwargs)
 
-    def set_param(self,param):
+    def _set_params(self,param):
         #print param.flags['C_CONTIGUOUS']
         self._param = param.copy()
 
-    def get_param(self):
+    def _get_params(self):
         return self._param
 
-    def get_param_names(self):
+    def _get_param_names(self):
         return [x.name for x in self._sp_theta]

GPy/kern/white.py

@@ -17,16 +17,16 @@ class white(kernpart):
         self.D = D
         self.Nparam = 1
         self.name = 'white'
-        self.set_param(np.array([variance]).flatten())
+        self._set_params(np.array([variance]).flatten())
 
-    def get_param(self):
+    def _get_params(self):
         return self.variance
 
-    def set_param(self,x):
+    def _set_params(self,x):
         assert x.shape==(1,)
         self.variance = x
 
-    def get_param_names(self):
+    def _get_param_names(self):
         return ['variance']
 
     def K(self,X,X2,target):

GPy/models/GPLVM.py

@@ -33,18 +33,18 @@ class GPLVM(GP_regression):
         else:
             return np.random.randn(Y.shape[0], Q)
 
-    def get_param_names(self):
+    def _get_param_names(self):
         return (sum([['X_%i_%i'%(n,q) for n in range(self.N)] for q in range(self.Q)],[])
-                + self.kern.extract_param_names())
+                + self.kern._get_param_names_transformed())
 
-    def get_param(self):
-        return np.hstack((self.X.flatten(), self.kern.extract_param()))
+    def _get_params(self):
+        return np.hstack((self.X.flatten(), self.kern._get_params_transformed()))
 
-    def set_param(self,x):
+    def _set_params(self,x):
         self.X = x[:self.X.size].reshape(self.N,self.Q).copy()
-        GP_regression.set_param(self, x[self.X.size:])
+        GP_regression._set_params(self, x[self.X.size:])
 
-    def log_likelihood_gradients(self):
+    def _log_likelihood_gradients(self):
         dL_dK = self.dL_dK()
 
         dL_dtheta = self.kern.dK_dtheta(dL_dK,self.X)

GPy/models/GP_EP.py

@@ -41,14 +41,14 @@ class GP_EP(model):
         self.K = self.kernel.K(self.X)
         model.__init__(self)
 
-    def set_param(self,p):
-        self.kernel.expand_param(p)
+    def _set_params(self,p):
+        self.kernel._set_params_transformed(p)
 
-    def get_param(self):
-        return self.kernel.extract_param()
+    def _get_params(self):
+        return self.kernel._get_params_transformed()
 
-    def get_param_names(self):
-        return self.kernel.extract_param_names()
+    def _get_param_names(self):
+        return self.kernel._get_param_names_transformed()
 
     def approximate_likelihood(self):
         self.ep_approx = Full(self.K,self.likelihood,epsilon=self.epsilon_ep,powerep=[self.eta,self.delta])
@@ -78,7 +78,7 @@ class GP_EP(model):
         L3 = sum(np.log(self.ep_approx.Z_hat))
         return L1 + L2A + L2B + L3
 
-    def log_likelihood_gradients(self):
+    def _log_likelihood_gradients(self):
         dK_dp = self.kernel.dK_dtheta(self.X)
         self.dK_dp = dK_dp
         aux1,info_1 = linalg.flapack.dtrtrs(self.L,np.dot(self.Sroot_tilde_K,self.ep_approx.v_tilde),lower=1)
@@ -138,7 +138,7 @@ class GP_EP(model):
         """
         self.epsilon_em = epsilon
         log_likelihood_change = self.epsilon_em + 1.
-        self.parameters_path = [self.kernel.get_param()]
+        self.parameters_path = [self.kernel._get_params()]
         self.approximate_likelihood()
         self.site_approximations_path = [[self.ep_approx.tau_tilde,self.ep_approx.v_tilde]]
         self.log_likelihood_path = [self.log_likelihood()]
@@ -150,11 +150,11 @@ class GP_EP(model):
             log_likelihood_change = log_likelihood_new - self.log_likelihood_path[-1]
             if log_likelihood_change < 0:
                 print 'log_likelihood decrement'
-                self.kernel.expand_param(self.parameters_path[-1])
-                self.kernM.expand_param(self.parameters_path[-1])
+                self.kernel._set_params_transformed(self.parameters_path[-1])
+                self.kernM._set_params_transformed(self.parameters_path[-1])
             else:
                 self.approximate_likelihood()
                 self.log_likelihood_path.append(self.log_likelihood())
-                self.parameters_path.append(self.kernel.get_param())
+                self.parameters_path.append(self.kernel._get_params())
                 self.site_approximations_path.append([self.ep_approx.tau_tilde,self.ep_approx.v_tilde])
             iteration += 1

GPy/models/GP_regression.py

@@ -63,47 +63,47 @@ class GP_regression(model):
             self._Ystd = np.ones((1,self.Y.shape[1]))
 
         if self.D > self.N:
-            # then it's more efficient to store Youter
-            self.Youter = np.dot(self.Y, self.Y.T)
+            # then it's more efficient to store YYT
+            self.YYT = np.dot(self.Y, self.Y.T)
         else:
-            self.Youter = None
+            self.YYT = None
 
         model.__init__(self)
 
-    def set_param(self,p):
-        self.kern.expand_param(p)
+    def _set_params(self,p):
+        self.kern._set_params_transformed(p)
         self.K = self.kern.K(self.X,slices1=self.Xslices)
         self.Ki, self.L, self.Li, self.K_logdet = pdinv(self.K)
 
-    def get_param(self):
-        return self.kern.extract_param()
+    def _get_params(self):
+        return self.kern._get_params_transformed()
 
-    def get_param_names(self):
-        return self.kern.extract_param_names()
+    def _get_param_names(self):
+        return self.kern._get_param_names_transformed()
 
     def _model_fit_term(self):
         """
-        Computes the model fit using Youter if it's available
+        Computes the model fit using YYT if it's available
         """
-        if self.Youter is None:
+        if self.YYT is None:
             return -0.5*np.sum(np.square(np.dot(self.Li,self.Y)))
         else:
-            return -0.5*np.sum(np.multiply(self.Ki, self.Youter))
+            return -0.5*np.sum(np.multiply(self.Ki, self.YYT))
 
     def log_likelihood(self):
         complexity_term = -0.5*self.N*self.D*np.log(2.*np.pi) - 0.5*self.D*self.K_logdet
         return complexity_term + self._model_fit_term()
 
     def dL_dK(self):
-        if self.Youter is None:
+        if self.YYT is None:
             alpha = np.dot(self.Ki,self.Y)
             dL_dK = 0.5*(np.dot(alpha,alpha.T)-self.D*self.Ki)
         else:
-            dL_dK = 0.5*(mdot(self.Ki, self.Youter, self.Ki) - self.D*self.Ki)
+            dL_dK = 0.5*(mdot(self.Ki, self.YYT, self.Ki) - self.D*self.Ki)
 
         return dL_dK
 
-    def log_likelihood_gradients(self):
+    def _log_likelihood_gradients(self):
         return self.kern.dK_dtheta(partial=self.dL_dK(),X=self.X)
 
     def predict(self,Xnew, slices=None, full_cov=False):

GPy/models/generalized_FITC.py

@@ -42,15 +42,15 @@ class generalized_FITC(model):
         self.jitter = 1e-12
         model.__init__(self)
 
-    def set_param(self,p):
-        self.kernel.expand_param(p[0:-self.Z.size])
+    def _set_params(self,p):
+        self.kernel._set_params_transformed(p[0:-self.Z.size])
         self.Z = p[-self.Z.size:].reshape(self.M,self.D)
 
-    def get_param(self):
-        return np.hstack([self.kernel.extract_param(),self.Z.flatten()])
+    def _get_params(self):
+        return np.hstack([self.kernel._get_params_transformed(),self.Z.flatten()])
 
-    def get_param_names(self):
-        return self.kernel.extract_param_names()+['iip_%i'%i for i in range(self.Z.size)]
+    def _get_param_names(self):
+        return self.kernel._get_param_names_transformed()+['iip_%i'%i for i in range(self.Z.size)]
 
     def approximate_likelihood(self):
         self.Kmm = self.kernel.K(self.Z)
@@ -91,15 +91,15 @@ class generalized_FITC(model):
 
     def log_likelihood(self):
         self.posterior_param()
-        self.Youter = np.dot(self.mu_tilde,self.mu_tilde.T)
+        self.YYT = np.dot(self.mu_tilde,self.mu_tilde.T)
         A = -self.hld
-        B = -.5*np.sum(self.Qi*self.Youter)
+        B = -.5*np.sum(self.Qi*self.YYT)
         C = sum(np.log(self.ep_approx.Z_hat))
         D = .5*np.sum(np.log(1./self.ep_approx.tau_tilde + 1./self.ep_approx.tau_))
         E = .5*np.sum((self.ep_approx.v_/self.ep_approx.tau_ - self.mu_tilde.flatten())**2/(1./self.ep_approx.tau_ + 1./self.ep_approx.tau_tilde))
         return  A + B + C + D + E
 
-    def log_likelihood_gradients(self):
+    def _log_likelihood_gradients(self):
         dKmm_dtheta = self.kernel.dK_dtheta(self.Z)
         dKnn_dtheta = self.kernel.dK_dtheta(self.X)
         dKmn_dtheta = self.kernel.dK_dtheta(self.Z,self.X)
@@ -214,7 +214,7 @@ class generalized_FITC(model):
         """
         self.epsilon_em = epsilon
         log_likelihood_change = self.epsilon_em + 1.
-        self.parameters_path = [self.kernel.get_param()]
+        self.parameters_path = [self.kernel._get_params()]
         self.approximate_likelihood()
         self.site_approximations_path = [[self.ep_approx.tau_tilde,self.ep_approx.v_tilde]]
         self.inducing_inputs_path = [self.Z]
@@ -227,7 +227,7 @@ class generalized_FITC(model):
             log_likelihood_change = log_likelihood_new - self.log_likelihood_path[-1]
             if log_likelihood_change < 0:
                 print 'log_likelihood decrement'
-                self.kernel.expand_param(self.parameters_path[-1])
+                self.kernel._set_params_transformed(self.parameters_path[-1])
                 self.kernM = self.kernel.copy()
                 slef.kernM.expand_X(self.iducing_inputs_path[-1])
                 self.__init__(self.kernel,self.likelihood,kernM=self.kernM,powerep=[self.eta,self.delta],epsilon_ep = self.epsilon_ep, epsilon_em = self.epsilon_em)
@@ -235,7 +235,7 @@ class generalized_FITC(model):
             else:
                 self.approximate_likelihood()
                 self.log_likelihood_path.append(self.log_likelihood())
-                self.parameters_path.append(self.kernel.get_param())
+                self.parameters_path.append(self.kernel._get_params())
                 self.site_approximations_path.append([self.ep_approx.tau_tilde,self.ep_approx.v_tilde])
                 self.inducing_inputs_path.append(self.Z)
             iteration += 1

GPy/models/sparse_GPLVM.py

@@ -27,16 +27,16 @@ class sparse_GPLVM(sparse_GP_regression, GPLVM):
         X = self.initialise_latent(init, Q, Y)
         sparse_GP_regression.__init__(self, X, Y, **kwargs)
 
-    def get_param_names(self):
+    def _get_param_names(self):
         return (sum([['X_%i_%i'%(n,q) for n in range(self.N)] for q in range(self.Q)],[])
-                + sparse_GP_regression.get_param_names(self))
+                + sparse_GP_regression._get_param_names(self))
 
-    def get_param(self):
-        return np.hstack((self.X.flatten(), sparse_GP_regression.get_param(self)))
+    def _get_params(self):
+        return np.hstack((self.X.flatten(), sparse_GP_regression._get_params(self)))
 
-    def set_param(self,x):
+    def _set_params(self,x):
         self.X = x[:self.X.size].reshape(self.N,self.Q).copy()
-        sparse_GP_regression.set_param(self, x[self.X.size:])
+        sparse_GP_regression._set_params(self, x[self.X.size:])
 
     def log_likelihood(self):
         return sparse_GP_regression.log_likelihood(self)
@@ -49,8 +49,8 @@ class sparse_GPLVM(sparse_GP_regression, GPLVM):
 
         return dL_dX
 
-    def log_likelihood_gradients(self):
-        return np.hstack((self.dL_dX().flatten(), sparse_GP_regression.log_likelihood_gradients(self)))
+    def _log_likelihood_gradients(self):
+        return np.hstack((self.dL_dX().flatten(), sparse_GP_regression._log_likelihood_gradients(self)))
 
     def plot(self):
         GPLVM.plot(self)

GPy/models/sparse_GP_regression.py

@@ -59,10 +59,10 @@ class sparse_GP_regression(GP_regression):
         if self.has_uncertain_inputs:
             self.X_uncertainty /= np.square(self._Xstd)
 
-    def set_param(self, p):
+    def _set_params(self, p):
         self.Z = p[:self.M*self.Q].reshape(self.M, self.Q)
         self.beta = p[self.M*self.Q]
-        self.kern.set_param(p[self.Z.size + 1:])
+        self.kern._set_params(p[self.Z.size + 1:])
         self.beta2 = self.beta**2
         self._compute_kernel_matrices()
         self._computations()
@@ -106,11 +106,11 @@ class sparse_GP_regression(GP_regression):
         self.dL_dKmm += -0.5 * self.D * (- self.LBL_inv - 2.*self.beta*mdot(self.LBL_inv, self.psi2, self.Kmmi) + self.Kmmi) # dC
         self.dL_dKmm +=  np.dot(np.dot(self.G,self.beta*self.psi2) - np.dot(self.LBL_inv, self.psi1VVpsi1), self.Kmmi) + 0.5*self.G # dE
 
-    def get_param(self):
-        return np.hstack([self.Z.flatten(),self.beta,self.kern.extract_param()])
+    def _get_params(self):
+        return np.hstack([self.Z.flatten(),self.beta,self.kern._get_params_transformed()])
 
-    def get_param_names(self):
-        return sum([['iip_%i_%i'%(i,j) for i in range(self.Z.shape[0])] for j in range(self.Z.shape[1])],[]) + ['noise_precision']+self.kern.extract_param_names()
+    def _get_param_names(self):
+        return sum([['iip_%i_%i'%(i,j) for i in range(self.Z.shape[0])] for j in range(self.Z.shape[1])],[]) + ['noise_precision']+self.kern._get_param_names_transformed()
 
     def log_likelihood(self):
         """
@@ -168,7 +168,7 @@ class sparse_GP_regression(GP_regression):
             dL_dZ += self.kern.dK_dX(dL_dpsi1,self.Z,self.X)
         return dL_dZ
 
-    def log_likelihood_gradients(self):
+    def _log_likelihood_gradients(self):
         return np.hstack([self.dL_dZ().flatten(), self.dL_dbeta(), self.dL_dtheta()])
 
     def _raw_predict(self, Xnew, slices, full_cov=False):

GPy/models/warped_GP.py

@@ -13,7 +13,7 @@ from GP_regression import GP_regression
 
 class warpedGP(GP_regression):
     """
-    TODO: fucking docstrings!
+    TODO: fecking docstrings!
 
     @nfusi: I'#ve hacked a little on this, but no guarantees. J.
     """
@@ -30,17 +30,17 @@ class warpedGP(GP_regression):
         self.transform_data()
         GP_regression.__init__(self, X, self.Y, **kwargs)
 
-    def set_param(self, x):
+    def _set_params(self, x):
         self.warping_params = x[:self.warping_function.num_parameters]
         self.transform_data()
-        GP_regression.set_param(self, x[self.warping_function.num_parameters:].copy())
+        GP_regression._set_params(self, x[self.warping_function.num_parameters:].copy())
 
-    def get_param(self):
-        return np.hstack((self.warping_params.flatten().copy(), GP_regression.get_param(self).copy()))
+    def _get_params(self):
+        return np.hstack((self.warping_params.flatten().copy(), GP_regression._get_params(self).copy()))
 
-    def get_param_names(self):
-        warping_names = self.warping_function.get_param_names()
-        param_names = GP_regression.get_param_names(self)
+    def _get_param_names(self):
+        warping_names = self.warping_function._get_param_names()
+        param_names = GP_regression._get_param_names(self)
         return warping_names + param_names
 
     def transform_data(self):
@@ -48,9 +48,9 @@ class warpedGP(GP_regression):
 
         # this supports the 'smart' behaviour in GP_regression
         if self.D > self.N:
-            self.Youter = np.dot(self.Y, self.Y.T)
+            self.YYT = np.dot(self.Y, self.Y.T)
         else:
-            self.Youter = None
+            self.YYT = None
 
         return self.Y
 
@@ -59,8 +59,8 @@ class warpedGP(GP_regression):
         jacobian = self.warping_function.fgrad_y(self.Z, self.warping_params)
         return ll + np.log(jacobian).sum()
 
-    def log_likelihood_gradients(self):
-        ll_grads = GP_regression.log_likelihood_gradients(self)
+    def _log_likelihood_gradients(self):
+        ll_grads = GP_regression._log_likelihood_gradients(self)
         alpha = np.dot(self.Ki, self.Y.flatten())
         warping_grads = self.warping_function_gradients(alpha)
 
@@ -83,7 +83,7 @@ class warpedGP(GP_regression):
     def predict(self, X, in_unwarped_space = False, **kwargs):
         mu, var = GP_regression.predict(self, X, **kwargs)
 
-        # The plot() function calls set_param() before calling predict()
+        # The plot() function calls _set_params() before calling predict()
         # this is causing the observations to be plotted in the transformed
         # space (where Y lives), making the plot looks very wrong
         # if the predictions are made in the untransformed space

GPy/testing/unit_tests.py

@@ -42,51 +42,66 @@ class GradientTests(unittest.TestCase):
         # contrain all parameters to be positive
         self.assertTrue(m.checkgrad())
 
-    def test_gp_regression_rbf_white_kern_1d(self):
+    def test_gp_regression_rbf_1d(self):
         ''' Testing the GP regression with rbf kernel with white kernel on 1d data '''
         rbf = GPy.kern.rbf(1)
         self.check_model_with_white(rbf, model_type='GP_regression', dimension=1)
 
-    def test_GP_regression_rbf_ARD_white_kern_2D(self):
-        ''' Testing the GP regression with rbf and white kernel on 2d data '''
-        k = GPy.kern.rbf_ARD(2)
-        self.check_model_with_white(k, model_type='GP_regression', dimension=2)
-
-    def test_GP_regression_rbf_white_kern_2D(self):
+    def test_GP_regression_rbf_2D(self):
         ''' Testing the GP regression with rbf and white kernel on 2d data '''
         rbf = GPy.kern.rbf(2)
         self.check_model_with_white(rbf, model_type='GP_regression', dimension=2)
 
-    def test_GP_regression_matern52_kern_1D(self):
+    def test_GP_regression_rbf_ARD_2D(self):
+        ''' Testing the GP regression with rbf and white kernel on 2d data '''
+        k = GPy.kern.rbf(2,ARD=True)
+        self.check_model_with_white(k, model_type='GP_regression', dimension=2)
+
+    def test_GP_regression_matern52_1D(self):
         ''' Testing the GP regression with matern52 kernel on 1d data '''
         matern52 = GPy.kern.Matern52(1)
         self.check_model_with_white(matern52, model_type='GP_regression', dimension=1)
 
-    def test_GP_regression_matern52_kern_2D(self):
+    def test_GP_regression_matern52_2D(self):
         ''' Testing the GP regression with matern52 kernel on 2d data '''
         matern52 = GPy.kern.Matern52(2)
         self.check_model_with_white(matern52, model_type='GP_regression', dimension=2)
 
-    def test_GP_regression_matern32_kern_1D(self):
+    def test_GP_regression_matern52_ARD_2D(self):
+        ''' Testing the GP regression with matern52 kernel on 2d data '''
+        matern52 = GPy.kern.Matern52(2,ARD=True)
+        self.check_model_with_white(matern52, model_type='GP_regression', dimension=2)
+
+    def test_GP_regression_matern32_1D(self):
         ''' Testing the GP regression with matern32 kernel on 1d data '''
         matern32 = GPy.kern.Matern32(1)
         self.check_model_with_white(matern32, model_type='GP_regression', dimension=1)
 
-    def test_GP_regression_matern32_kern_2D(self):
+    def test_GP_regression_matern32_2D(self):
         ''' Testing the GP regression with matern32 kernel on 2d data '''
         matern32 = GPy.kern.Matern32(2)
         self.check_model_with_white(matern32, model_type='GP_regression', dimension=2)
 
-    def test_GP_regression_exponential_kern_1D(self):
+    def test_GP_regression_matern32_ARD_2D(self):
+        ''' Testing the GP regression with matern32 kernel on 2d data '''
+        matern32 = GPy.kern.Matern32(2,ARD=True)
+        self.check_model_with_white(matern32, model_type='GP_regression', dimension=2)
+
+    def test_GP_regression_exponential_1D(self):
         ''' Testing the GP regression with exponential kernel on 1d data '''
         exponential = GPy.kern.exponential(1)
         self.check_model_with_white(exponential, model_type='GP_regression', dimension=1)
 
-    def test_GP_regression_exponential_kern_2D(self):
+    def test_GP_regression_exponential_2D(self):
         ''' Testing the GP regression with exponential kernel on 2d data '''
         exponential = GPy.kern.exponential(2)
         self.check_model_with_white(exponential, model_type='GP_regression', dimension=2)
 
+    def test_GP_regression_exponential_ARD_2D(self):
+        ''' Testing the GP regression with exponential kernel on 2d data '''
+        exponential = GPy.kern.exponential(2,ARD=True)
+        self.check_model_with_white(exponential, model_type='GP_regression', dimension=2)
+
     def test_GP_regression_bias_kern_1D(self):
         ''' Testing the GP regression with bias kernel on 1d data '''
         bias = GPy.kern.bias(1)
@@ -121,7 +136,7 @@ class GradientTests(unittest.TestCase):
         """ Testing GPLVM with rbf + bias and white kernel """
         N, Q, D = 50, 1, 2
         X = np.random.rand(N, Q)
-        k = GPy.kern.rbf(Q, 0.5, 0.9) + GPy.kern.bias(Q, 0.1) + GPy.kern.white(Q, 0.05)
+        k = GPy.kern.rbf(Q, 0.5, 0.9*np.ones((1,))) + GPy.kern.bias(Q, 0.1) + GPy.kern.white(Q, 0.05)
         K = k.K(X)
         Y = np.random.multivariate_normal(np.zeros(N),K,D).T
         m = GPy.models.GPLVM(Y, Q, kernel = k)

GPy/util/datasets.py

@@ -90,7 +90,7 @@ def toy_rbf_1d(seed=default_seed):
     N = 500
     X = np.random.uniform(low=-1.0, high=1.0, size=(N, numIn))
     X.sort(axis=0)
-    rbf = GPy.kern.rbf(numIn, variance=1., lengthscale=0.25)
+    rbf = GPy.kern.rbf(numIn, variance=1., lengthscale=np.array((0.25,)))
     white = GPy.kern.white(numIn, variance=1e-2)
     kernel = rbf + white
     K = kernel.K(X)

GPy/util/warping_functions.py

@@ -33,7 +33,7 @@ class WarpingFunction(object):
         """inverse function transformation"""
         raise NotImplementedError
 
-    def get_param_names(self):
+    def _get_param_names(self):
         raise NotImplementedError
 
     def plot(self, psi, xmin, xmax):
@@ -151,7 +151,7 @@ class TanhWarpingFunction(WarpingFunction):
 
         return gradients
 
-    def get_param_names(self):
+    def _get_param_names(self):
         variables = ['a', 'b', 'c']
         names = sum([['warp_tanh_%s_t%i' % (variables[n],q) for n in range(3)] for q in range(self.n_terms)],[])
         return names