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

GPy/core/gp_base.py

@@ -6,8 +6,8 @@ from GPy.core.model import Model
 
 class GPBase(Model):
     """
-    Gaussian Process Model for holding shared behaviour between
-    sprase_GP and GP models
+    Gaussian process base model for holding shared behaviour between
+    sparse_GP and GP models.
     """
 
     def __init__(self, X, likelihood, kernel, normalize_X=False):
@@ -35,8 +35,7 @@ class GPBase(Model):
 
     def getstate(self):
         """
-        Get the current state of the class,
-        here just all the indices, rest can get recomputed
+        Get the current state of the class, here we return everything that is needed to recompute the model.
         """
         return Model.getstate(self) + [self.X,
                 self.num_data,
@@ -63,14 +62,6 @@ class GPBase(Model):
         Plot the GP's view of the world, where the data is normalized and the
         likelihood is Gaussian.
 
-        :param samples: the number of a posteriori samples to plot
-        :param which_data: which if the training data to plot (default all)
-        :type which_data: 'all' or a slice object to slice self.X, self.Y
-        :param plot_limits: The limits of the plot. If 1D [xmin,xmax], if 2D [[xmin,ymin],[xmax,ymax]]. Defaluts to data limits
-        :param which_parts: which of the kernel functions to plot (additively)
-        :type which_parts: 'all', or list of bools
-        :param resolution: the number of intervals to sample the GP on. Defaults to 200 in 1D and 50 (a 50x50 grid) in 2D
-
         Plot the posterior of the GP.
           - In one dimension, the function is plotted with a shaded region identifying two standard deviations.
           - In two dimsensions, a contour-plot shows the mean predicted function
@@ -78,6 +69,22 @@ class GPBase(Model):
 
         Can plot only part of the data and part of the posterior functions
         using which_data and which_functions
+
+        :param samples: the number of a posteriori samples to plot
+        :param plot_limits: The limits of the plot. If 1D [xmin,xmax], if 2D [[xmin,ymin],[xmax,ymax]]. Defaluts to data limits
+        :param which_data: which if the training data to plot (default all)
+        :type which_data: 'all' or a slice object to slice self.X, self.Y
+        :param which_parts: which of the kernel functions to plot (additively)
+        :type which_parts: 'all', or list of bools
+        :param resolution: the number of intervals to sample the GP on. Defaults to 200 in 1D and 50 (a 50x50 grid) in 2D
+        :type resolution: int
+        :param full_cov:
+        :type full_cov: bool
+                :param fignum: figure to plot on.
+        :type fignum: figure number
+        :param ax: axes to plot on.
+        :type ax: axes handle
+
         """
         if which_data == 'all':
             which_data = slice(None)
@@ -118,12 +125,41 @@ class GPBase(Model):
 
     def plot(self, plot_limits=None, which_data='all', which_parts='all', resolution=None, levels=20, samples=0, fignum=None, ax=None, fixed_inputs=[], linecol=Tango.colorsHex['darkBlue'],fillcol=Tango.colorsHex['lightBlue']):
         """
-        TODO: Docstrings!
+        Plot the GP with noise where the likelihood is Gaussian.
 
+        Plot the posterior of the GP.
+          - In one dimension, the function is plotted with a shaded region identifying two standard deviations.
+          - In two dimsensions, a contour-plot shows the mean predicted function
+          - In higher dimensions, we've no implemented this yet !TODO!
+
+        Can plot only part of the data and part of the posterior functions
+        using which_data and which_functions
+
+        :param plot_limits: The limits of the plot. If 1D [xmin,xmax], if 2D [[xmin,ymin],[xmax,ymax]]. Defaluts to data limits
+        :type plot_limits: np.array
+        :param which_data: which if the training data to plot (default all)
+        :type which_data: 'all' or a slice object to slice self.X, self.Y
+        :param which_parts: which of the kernel functions to plot (additively)
+        :type which_parts: 'all', or list of bools
+        :param resolution: the number of intervals to sample the GP on. Defaults to 200 in 1D and 50 (a 50x50 grid) in 2D
+        :type resolution: int
+        :param levels: number of levels to plot in a contour plot.
+        :type levels: int
+        :param samples: the number of a posteriori samples to plot
+        :type samples: int
+        :param fignum: figure to plot on.
+        :type fignum: figure number
+        :param ax: axes to plot on.
+        :type ax: axes handle
+        :param fixed_inputs: a list of tuple [(i,v), (i,v)...], specifying that input index i should be set to value v.
+        :type fixed_inputs: a list of tuples
+        :param linecol: color of line to plot.
+        :type linecol:
+        :param fillcol: color of fill
+        :type fillcol: 
         :param levels: for 2D plotting, the number of contour levels to use
         is ax is None, create a new figure
 
-        fixed_inputs: a list of tuple [(i,v), (i,v)...], specifying that input index i should be set to value v.
         """
         # TODO include samples
         if which_data == 'all':

GPy/core/model.py

@@ -24,15 +24,17 @@ class Model(Parameterized):
         self.preferred_optimizer = 'scg'
         # self._set_params(self._get_params()) has been taken out as it should only be called on leaf nodes
     def log_likelihood(self):
-        raise NotImplementedError, "this needs to be implemented to use the Model class"
+        raise NotImplementedError, "this needs to be implemented to use the model class"
     def _log_likelihood_gradients(self):
-        raise NotImplementedError, "this needs to be implemented to use the Model class"
+        raise NotImplementedError, "this needs to be implemented to use the model class"
 
     def getstate(self):
         """
         Get the current state of the class.
         
         Inherited from Parameterized, so add those parameters to the state
+        :return: list of states from the model.
+
         """
         return Parameterized.getstate(self) + \
             [self.priors, self.optimization_runs,
@@ -41,8 +43,10 @@ class Model(Parameterized):
     def setstate(self, state):
         """
         set state from previous call to getstate
-        
         call Parameterized with the rest of the state
+
+        :param state: the state of the model.
+        :type state: list as returned from getstate.        
         """
         self.preferred_optimizer = state.pop()
         self.sampling_runs = state.pop()
@@ -52,21 +56,22 @@ class Model(Parameterized):
 
     def set_prior(self, regexp, what):
         """
-        Sets priors on the Model parameters.
-
-        Arguments
-        ---------
-        regexp -- string, regexp, or integer array
-        what -- instance of a Prior class
+        Sets priors on the model parameters.
 
         Notes
         -----
-        Asserts that the Prior is suitable for the constraint. If the
+        Asserts that the prior is suitable for the constraint. If the
         wrong constraint is in place, an error is raised.  If no
         constraint is in place, one is added (warning printed).
 
-        For tied parameters, the Prior will only be "counted" once, thus
-        a Prior object is only inserted on the first tied index
+        For tied parameters, the prior will only be "counted" once, thus
+        a prior object is only inserted on the first tied index
+
+        :param regexp: regular expression of parameters on which priors need to be set.
+        :type param: string, regexp, or integer array
+        :param what: prior to set on parameter.
+        :type what: GPy.core.Prior type
+
         """
         if self.priors is None:
             self.priors = [None for i in range(self._get_params().size)]
@@ -76,12 +81,12 @@ class Model(Parameterized):
         # check tied situation
         tie_partial_matches = [tie for tie in self.tied_indices if (not set(tie).isdisjoint(set(which))) & (not set(tie) == set(which))]
         if len(tie_partial_matches):
-            raise ValueError, "cannot place Prior across partial ties"
+            raise ValueError, "cannot place prior across partial ties"
         tie_matches = [tie for tie in self.tied_indices if set(which) == set(tie) ]
         if len(tie_matches) > 1:
-            raise ValueError, "cannot place Prior across multiple ties"
+            raise ValueError, "cannot place prior across multiple ties"
         elif len(tie_matches) == 1:
-            which = which[:1] # just place a Prior object on the first parameter
+            which = which[:1] # just place a prior object on the first parameter
 
 
         # check constraints are okay
@@ -93,7 +98,7 @@ class Model(Parameterized):
             else:
                 constrained_positive_indices = np.zeros(shape=(0,))
             bad_constraints = np.setdiff1d(self.all_constrained_indices(), constrained_positive_indices)
-            assert not np.any(which[:, None] == bad_constraints), "constraint and Prior incompatible"
+            assert not np.any(which[:, None] == bad_constraints), "constraint and prior incompatible"
             unconst = np.setdiff1d(which, constrained_positive_indices)
             if len(unconst):
                 print "Warning: constraining parameters to be positive:"
@@ -101,17 +106,22 @@ class Model(Parameterized):
                 print '\n'
                 self.constrain_positive(unconst)
         elif what.domain is REAL:
-            assert not np.any(which[:, None] == self.all_constrained_indices()), "constraint and Prior incompatible"
+            assert not np.any(which[:, None] == self.all_constrained_indices()), "constraint and prior incompatible"
         else:
-            raise ValueError, "Prior not recognised"
+            raise ValueError, "prior not recognised"
 
-        # store the Prior in a local list
+        # store the prior in a local list
         for w in which:
             self.priors[w] = what
 
     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.
+        Get model gradient(s) by name. The name is applied as a regular expression and all parameters that match that regular expression are returned.
+
+        :param name: the name of parameters required (as a regular expression).
+        :type name: regular expression
+        :param return_names: whether or not to return the names matched (default False)
+        :type return_names: bool
         """
         matches = self.grep_param_names(name)
         if len(matches):
@@ -151,14 +161,14 @@ class Model(Parameterized):
 
     def randomize(self):
         """
-        Randomize the Model.
-        Make this draw from the Prior if one exists, else draw from N(0,1)
+        Randomize the model.
+        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._get_params_transformed()
         x = np.random.randn(x.size)
         self._set_params_transformed(x)
-        # now draw from Prior where possible
+        # now draw from prior where possible
         x = self._get_params()
         if self.priors is not None:
             [np.put(x, i, p.rvs(1)) for i, p in enumerate(self.priors) if not p is None]
@@ -168,21 +178,30 @@ class Model(Parameterized):
 
     def optimize_restarts(self, num_restarts=10, robust=False, verbose=True, parallel=False, num_processes=None, **kwargs):
         """
-        Perform random restarts of the Model, and set the Model to the best
+        Perform random restarts of the model, and set the model to the best
         seen solution.
 
         If the robust flag is set, exceptions raised during optimizations will
-        be handled silently.  If _all_ runs fail, the Model is reset to the
+        be handled silently.  If _all_ runs fail, the model is reset to the
         existing parameter values.
 
         Notes
         -----
+        :param num_restarts: number of restarts to use (default 10)
+        :type num_restarts: int
+        :param robust: whether to handle exceptions silently or not (default False)
+        :type robust: bool
+        :param parallel: whether to run each restart as a separate process. It relies on the multiprocessing module.
+        :type parallel: bool
+        :param num_processes: number of workers in the multiprocessing pool
+        :type numprocesses: int
         **kwargs are passed to the optimizer. They can be:
-        :max_f_eval: maximum number of function evaluations
-        :messages: whether to display during optimisation
-        :verbose: whether to show informations about the current restart
-        :parallel: whether to run each restart as a separate process. It relies on the multiprocessing module.
-        :num_processes: number of workers in the multiprocessing pool
+        :param max_f_eval: maximum number of function evaluations
+        :type max_f_eval: int
+        :param max_iters: maximum number of iterations
+        :type max_iters: int
+        :param messages: whether to display during optimisation
+        :type messages: bool
 
         ..Note: If num_processes is None, the number of workes in the multiprocessing pool is automatically
         set to the number of processors on the current machine.
@@ -231,7 +250,12 @@ class Model(Parameterized):
 
     def ensure_default_constraints(self):
         """       
-        Ensure that any variables which should clearly be positive have been constrained somehow.
+        Ensure that any variables which should clearly be positive
+        have been constrained somehow. The method performs a regular
+        expression search on parameter names looking for the terms
+        'variance', 'lengthscale', 'precision' and 'kappa'. If any of
+        these terms are present in the name the parameter is
+        constrained positive.
         """
         positive_strings = ['variance', 'lengthscale', 'precision', 'kappa']
         # param_names = self._get_param_names()
@@ -246,11 +270,15 @@ class Model(Parameterized):
 
     def objective_function(self, x):
         """
-        The objective function passed to the optimizer. It combines the likelihood and the priors.
+        The objective function passed to the optimizer. It combines
+        the likelihood and the priors.
         
         Failures are handled robustly. The algorithm will try several times to
         return the objective, and will raise the original exception if it
         the objective cannot be computed.
+
+        :param x: the parameters of the model.
+        :parameter type: np.array
         """
         try:
             self._set_params_transformed(x)
@@ -267,8 +295,11 @@ class Model(Parameterized):
         Gets the gradients from the likelihood and the priors.
 
         Failures are handled robustly. The algorithm will try several times to
-        return the objective, and will raise the original exception if it
+        return the gradients, and will raise the original exception if it
         the objective cannot be computed.
+
+        :param x: the parameters of the model.
+        :parameter type: np.array
         """
         try:
             self._set_params_transformed(x)
@@ -282,6 +313,13 @@ class Model(Parameterized):
         return obj_grads
 
     def objective_and_gradients(self, x):
+        """
+        Compute the objective function of the model and the gradient of the model at the point given by x.
+
+        :param x: the point at which gradients are to be computed.
+        :type np.array:
+        """
+
         try:
             self._set_params_transformed(x)
             obj_f = -self.log_likelihood() - self.log_prior()
@@ -297,11 +335,13 @@ class Model(Parameterized):
 
     def optimize(self, optimizer=None, start=None, **kwargs):
         """
-        Optimize the Model using self.log_likelihood and self.log_likelihood_gradient, as well as self.priors.
+        Optimize the model using self.log_likelihood and self.log_likelihood_gradient, as well as self.priors.
         kwargs are passed to the optimizer. They can be:
 
-        :max_f_eval: maximum number of function evaluations
+        :param max_f_eval: maximum number of function evaluations
+        :type max_f_eval: int
         :messages: whether to display during optimisation
+        :type messages: bool
         :param optimzer: which optimizer to use (defaults to self.preferred optimizer)
         :type optimzer: string TODO: valid strings?
         """
@@ -327,14 +367,14 @@ class Model(Parameterized):
         self.optimization_runs.append(sgd)
 
     def Laplace_covariance(self):
-        """return the covariance matric of a Laplace approximatino at the current (stationary) point"""
-        # TODO add in the Prior contributions for MAP estimation
+        """return the covariance matrix of a Laplace approximation at the current (stationary) point."""
+        # TODO add in the prior contributions for MAP estimation
         # TODO fix the hessian for tied, constrained and fixed components
         if hasattr(self, 'log_likelihood_hessian'):
             A = -self.log_likelihood_hessian()
 
         else:
-            print "numerically calculating hessian. please be patient!"
+            print "numerically calculating Hessian. please be patient!"
             x = self._get_params()
             def f(x):
                 self._set_params(x)
@@ -348,8 +388,8 @@ class Model(Parameterized):
         return A
 
     def Laplace_evidence(self):
-        """Returns an estiamte of the Model evidence based on the Laplace approximation.
-        Uses a numerical estimate of the hessian if none is available analytically"""
+        """Returns an estiamte of the model evidence based on the Laplace approximation.
+        Uses a numerical estimate of the Hessian if none is available analytically."""
         A = self.Laplace_covariance()
         try:
             hld = np.sum(np.log(np.diag(jitchol(A)[0])))
@@ -370,27 +410,28 @@ class Model(Parameterized):
         log_prior = self.log_prior()
         obj_funct = '\nLog-likelihood: {0:.3e}'.format(log_like)
         if len(''.join(strs)) != 0:
-            obj_funct += ', Log Prior: {0:.3e}, LL+Prior = {0:.3e}'.format(log_prior, log_like + log_prior)
+            obj_funct += ', Log prior: {0:.3e}, LL+prior = {0:.3e}'.format(log_prior, log_like + log_prior)
         obj_funct += '\n\n'
         s[0] = obj_funct + s[0]
-        s[0] += "|{h:^{col}}".format(h='Prior', col=width)
+        s[0] += "|{h:^{col}}".format(h='prior', col=width)
         s[1] += '-' * (width + 1)
 
         for p in range(2, len(strs) + 2):
-            s[p] += '|{Prior:^{width}}'.format(Prior=strs[p - 2], width=width)
+            s[p] += '|{prior:^{width}}'.format(prior=strs[p - 2], width=width)
 
         return '\n'.join(s)
 
 
     def checkgrad(self, target_param=None, verbose=False, step=1e-6, tolerance=1e-3):
         """
-        Check the gradient of the Model by comparing to a numerical estimate.
-        If the verbose flag is passed, invividual components are tested (and printed)
+        Check the gradient of the ,odel by comparing to a numerical
+        estimate.  If the verbose flag is passed, invividual
+        components are tested (and printed)
 
         :param verbose: If True, print a "full" checking of each parameter
         :type verbose: bool
         :param step: The size of the step around which to linearise the objective
-        :type step: float (defaul 1e-6)
+        :type step: float (default 1e-6)
         :param tolerance: the tolerance allowed (see note)
         :type tolerance: float (default 1e-3)
 
@@ -467,15 +508,15 @@ class Model(Parameterized):
 
     def input_sensitivity(self):
         """
-        return an array describing the sesitivity of the Model to each input
+        return an array describing the sesitivity of the model to each input
 
         NB. Right now, we're basing this on the lengthscales (or
         variances) of the kernel.  TODO: proper sensitivity analysis
-        where we integrate across the Model inputs and evaluate the
-        effect on the variance of the Model output.  """
+        where we integrate across the model inputs and evaluate the
+        effect on the variance of the model output.  """
 
         if not hasattr(self, 'kern'):
-            raise ValueError, "this Model has no kernel"
+            raise ValueError, "this model has no kernel"
 
         k = [p for p in self.kern.parts if p.name in ['rbf', 'linear', 'rbf_inv']]
         if (not len(k) == 1) or (not k[0].ARD):

GPy/core/transformations.py

@@ -4,6 +4,8 @@
 
 import numpy as np
 from GPy.core.domains import POSITIVE, NEGATIVE, BOUNDED
+import sys 
+lim_val = -np.log(sys.float_info.epsilon) 
 
 class transformation(object):
     domain = None
@@ -17,7 +19,7 @@ class transformation(object):
         """ df_dx evaluated at self.f(x)=f"""
         raise NotImplementedError
     def initialize(self, f):
-        """ produce a sensible initial values for f(x)"""
+        """ produce a sensible initial value for f(x)"""
         raise NotImplementedError
     def __str__(self):
         raise NotImplementedError
@@ -25,13 +27,14 @@ class transformation(object):
 class logexp(transformation):
     domain = POSITIVE
     def f(self, x):
-        return np.log(1. + np.exp(x))
+        return np.where(x>lim_val, x, np.log(1. + np.exp(x)))
     def finv(self, f):
-        return np.log(np.exp(f) - 1.)
+        return np.where(f>lim_val, f, np.log(np.exp(f) - 1.))
     def gradfactor(self, f):
-        ef = np.exp(f)
-        return (ef - 1.) / ef
+        return np.where(f>lim_val, 1., 1 - np.exp(-f))
     def initialize(self, f):
+        if np.any(f < 0.):
+            print "Warning: changing parameters to satisfy constraints"
         return np.abs(f)
     def __str__(self):
         return '(+ve)'
@@ -39,18 +42,19 @@ class logexp(transformation):
 class negative_logexp(transformation):
     domain = NEGATIVE
     def f(self, x):
-        return -np.log(1. + np.exp(x))
+        return -logexp.f(x) #np.log(1. + np.exp(x))
     def finv(self, f):
-        return np.log(np.exp(-f) - 1.)
+        return logexp.finv(-f) #np.log(np.exp(-f) - 1.)
     def gradfactor(self, f):
-        ef = np.exp(-f)
-        return -(ef - 1.) / ef
+        return -logexp.gradfactor(-f)
+        #ef = np.exp(-f)
+        #return -(ef - 1.) / ef
     def initialize(self, f):
-        return -np.abs(f)
+        return -logexp.initialize(f) #np.abs(f)
     def __str__(self):
         return '(-ve)'
 
-class logexp_clipped(transformation):
+class logexp_clipped(logexp):
     max_bound = 1e100
     min_bound = 1e-10
     log_max_bound = np.log(max_bound)
@@ -78,9 +82,10 @@ class logexp_clipped(transformation):
         return '(+ve_c)'
 
 class exponent(transformation):
+    # TODO: can't allow this to go to zero, need to set a lower bound. Similar with negative exponent below. See old MATLAB code.
     domain = POSITIVE
     def f(self, x):
-        return np.exp(x)
+        return np.where(x<lim_val, np.where(x>-lim_val, np.exp(x), np.exp(-lim_val)), np.exp(lim_val))
     def finv(self, x):
         return np.log(x)
     def gradfactor(self, f):
@@ -92,18 +97,16 @@ class exponent(transformation):
     def __str__(self):
         return '(+ve)'
 
-class negative_exponent(transformation):
+class negative_exponent(exponent):
     domain = NEGATIVE
     def f(self, x):
-        return -np.exp(x)
-    def finv(self, x):
-        return np.log(-x)
+        return -exponent.f(x)
+    def finv(self, f):
+        return exponent.finv(-f)
     def gradfactor(self, f):
         return f
     def initialize(self, f):
-        if np.any(f > 0.):
-            print "Warning: changing parameters to satisfy constraints"
-        return -np.abs(f)
+        return -exponent.initialize(f) #np.abs(f)
     def __str__(self):
         return '(-ve)'
 

GPy/examples/classification.py

@@ -10,7 +10,7 @@ import numpy as np
 import GPy
 
 default_seed = 10000
-def crescent_data(seed=default_seed): # FIXME
+def crescent_data(seed=default_seed, kernel=None): # FIXME
     """Run a Gaussian process classification on the crescent data. The demonstration calls the basic GP classification model and uses EP to approximate the likelihood.
 
     :param model_type: type of model to fit ['Full', 'FITC', 'DTC'].
@@ -32,33 +32,33 @@ def crescent_data(seed=default_seed): # FIXME
     m.plot()
     return m
 
-def oil(num_inducing=50):
+def oil(num_inducing=50, max_iters=100, kernel=None):
     """
-    Run a Gaussian process classification on the oil data. The demonstration calls the basic GP classification model and uses EP to approximate the likelihood.
+    Run a Gaussian process classification on the three phase oil data. The demonstration calls the basic GP classification model and uses EP to approximate the likelihood.
     """
     data = GPy.util.datasets.oil()
-    X = data['X'][:600,:]
-    X_test = data['X'][600:,:]
-    Y = data['Y'][:600, 0:1]
+    X = data['X']
+    Xtest = data['Xtest']
+    Y = data['Y'][:, 0:1]
+    Ytest = data['Ytest'][:, 0:1]
     Y[Y.flatten()==-1] = 0
-    Y_test = data['Y'][600:, 0:1]
+    Ytest[Ytest.flatten()==-1] = 0
 
     # Create GP model
-    m = GPy.models.SparseGPClassification(X, Y,num_inducing=num_inducing)
+    m = GPy.models.SparseGPClassification(X, Y,kernel=kernel,num_inducing=num_inducing)
 
     # Contrain all parameters to be positive
-    m.constrain_positive('')
     m.tie_params('.*len')
     m['.*len'] = 10.
     m.update_likelihood_approximation()
 
     # Optimize
-    m.optimize()
+    m.optimize(max_iters=max_iters)
     print(m)
 
     #Test
-    probs = m.predict(X_test)[0]
-    GPy.util.classification.conf_matrix(probs,Y_test)
+    probs = m.predict(Xtest)[0]
+    GPy.util.classification.conf_matrix(probs,Ytest)
     return m
 
 def toy_linear_1d_classification(seed=default_seed):
@@ -118,7 +118,7 @@ def sparse_toy_linear_1d_classification(num_inducing=10,seed=default_seed):
 
     return m
 
-def sparse_crescent_data(num_inducing=10, seed=default_seed):
+def sparse_crescent_data(num_inducing=10, seed=default_seed, kernel=None):
     """
     Run a Gaussian process classification with DTC approxiamtion on the crescent data. The demonstration calls the basic GP classification model and uses EP to approximate the likelihood.
 
@@ -133,7 +133,7 @@ def sparse_crescent_data(num_inducing=10, seed=default_seed):
     Y = data['Y']
     Y[Y.flatten()==-1]=0
 
-    m = GPy.models.SparseGPClassification(data['X'], Y,num_inducing=num_inducing)
+    m = GPy.models.SparseGPClassification(data['X'], Y, kernel=kernel, num_inducing=num_inducing)
     m['.*len'] = 10.
     #m.update_likelihood_approximation()
     #m.optimize()

GPy/examples/dimensionality_reduction.py

@@ -352,7 +352,7 @@ def brendan_faces():
 
     return m
 def stick_play(range=None, frame_rate=15):
-    data = GPy.util.datasets.stick()
+    data = GPy.util.datasets.osu_run1()
     # optimize
     if range == None:
         Y = data['Y'].copy()
@@ -364,7 +364,7 @@ def stick_play(range=None, frame_rate=15):
     return Y
 
 def stick():
-    data = GPy.util.datasets.stick()
+    data = GPy.util.datasets.osu_run1()
     # optimize
     m = GPy.models.GPLVM(data['Y'], 2)
     m.optimize(messages=1, max_f_eval=10000)
@@ -378,8 +378,19 @@ def stick():
 
     return m
 
+def robot_wireless():
+    data = GPy.util.datasets.robot_wireless()
+    # optimize
+    m = GPy.models.GPLVM(data['Y'], 2)
+    m.optimize(messages=1, max_f_eval=10000)
+    m._set_params(m._get_params())
+    plt.clf
+    ax = m.plot_latent()
+
+    return m
+
 def stick_bgplvm(model=None):
-    data = GPy.util.datasets.stick()
+    data = GPy.util.datasets.osu_run1()
     Q = 6
     kernel = GPy.kern.rbf(Q, ARD=True) + GPy.kern.bias(Q, np.exp(-2)) + GPy.kern.white(Q, np.exp(-2))
     m = BayesianGPLVM(data['Y'], Q, init="PCA", num_inducing=20, kernel=kernel)

GPy/examples/regression.py

@@ -24,25 +24,45 @@ def toy_rbf_1d(optimizer='tnc', max_nb_eval_optim=100):
     print(m)
     return m
 
-def rogers_girolami_olympics(optim_iters=100):
+def olympic_100m_men(max_iters=100, kernel=None):
     """Run a standard Gaussian process regression on the Rogers and Girolami olympics data."""
-    data = GPy.util.datasets.rogers_girolami_olympics()
+    data = GPy.util.datasets.olympic_100m_men()
 
     # create simple GP Model
-    m = GPy.models.GPRegression(data['X'], data['Y'])
+    m = GPy.models.GPRegression(data['X'], data['Y'], kernel)
 
     # set the lengthscale to be something sensible (defaults to 1)
+    if kernel==None:
         m['rbf_lengthscale'] = 10
 
     # optimize
-    m.optimize(max_f_eval=optim_iters)
+    m.optimize(max_iters=max_iters)
 
     # plot
     m.plot(plot_limits=(1850, 2050))
     print(m)
     return m
 
-def toy_rbf_1d_50(optim_iters=100):
+def olympic_marathon_men(max_iters=100, kernel=None):
+    """Run a standard Gaussian process regression on the Olympic marathon data."""
+    data = GPy.util.datasets.olympic_marathon_men()
+
+    # create simple GP Model
+    m = GPy.models.GPRegression(data['X'], data['Y'], kernel)
+
+    # set the lengthscale to be something sensible (defaults to 1)
+    if kernel==None:
+        m['rbf_lengthscale'] = 10
+
+    # optimize
+    m.optimize(max_iters=max_iters)
+
+    # plot
+    m.plot(plot_limits=(1850, 2050))
+    print(m)
+    return m
+
+def toy_rbf_1d_50(max_iters=100):
     """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()
 
@@ -50,21 +70,21 @@ def toy_rbf_1d_50(optim_iters=100):
     m = GPy.models.GPRegression(data['X'], data['Y'])
 
     # optimize
-    m.optimize(max_f_eval=optim_iters)
+    m.optimize(max_iters=max_iters)
 
     # plot
     m.plot()
     print(m)
     return m
 
-def toy_ARD(optim_iters=1000, kernel_type='linear', N=300, D=4):
+def toy_ARD(max_iters=1000, kernel_type='linear', num_samples=300, D=4):
     # Create an artificial dataset where the values in the targets (Y)
     # only depend in dimensions 1 and 3 of the inputs (X). Run ARD to
     # see if this dependency can be recovered
-    X1 = np.sin(np.sort(np.random.rand(N, 1) * 10, 0))
-    X2 = np.cos(np.sort(np.random.rand(N, 1) * 10, 0))
-    X3 = np.exp(np.sort(np.random.rand(N, 1), 0))
-    X4 = np.log(np.sort(np.random.rand(N, 1), 0))
+    X1 = np.sin(np.sort(np.random.rand(num_samples, 1) * 10, 0))
+    X2 = np.cos(np.sort(np.random.rand(num_samples, 1) * 10, 0))
+    X3 = np.exp(np.sort(np.random.rand(num_samples, 1), 0))
+    X4 = np.log(np.sort(np.random.rand(num_samples, 1), 0))
     X = np.hstack((X1, X2, X3, X4))
 
     Y1 = np.asarray(2 * X[:, 0] + 3).reshape(-1, 1)
@@ -87,20 +107,20 @@ def toy_ARD(optim_iters=1000, kernel_type='linear', N=300, D=4):
     # len_prior = GPy.priors.inverse_gamma(1,18) # 1, 25
     # m.set_prior('.*lengthscale',len_prior)
 
-    m.optimize(optimizer='scg', max_iters=optim_iters, messages=1)
+    m.optimize(optimizer='scg', max_iters=max_iters, messages=1)
 
     m.kern.plot_ARD()
     print(m)
     return m
 
-def toy_ARD_sparse(optim_iters=1000, kernel_type='linear', N=300, D=4):
+def toy_ARD_sparse(max_iters=1000, kernel_type='linear', num_samples=300, D=4):
     # Create an artificial dataset where the values in the targets (Y)
     # only depend in dimensions 1 and 3 of the inputs (X). Run ARD to
     # see if this dependency can be recovered
-    X1 = np.sin(np.sort(np.random.rand(N, 1) * 10, 0))
-    X2 = np.cos(np.sort(np.random.rand(N, 1) * 10, 0))
-    X3 = np.exp(np.sort(np.random.rand(N, 1), 0))
-    X4 = np.log(np.sort(np.random.rand(N, 1), 0))
+    X1 = np.sin(np.sort(np.random.rand(num_samples, 1) * 10, 0))
+    X2 = np.cos(np.sort(np.random.rand(num_samples, 1) * 10, 0))
+    X3 = np.exp(np.sort(np.random.rand(num_samples, 1), 0))
+    X4 = np.log(np.sort(np.random.rand(num_samples, 1), 0))
     X = np.hstack((X1, X2, X3, X4))
 
     Y1 = np.asarray(2 * X[:, 0] + 3)[:, None]
@@ -124,13 +144,13 @@ def toy_ARD_sparse(optim_iters=1000, kernel_type='linear', N=300, D=4):
     # len_prior = GPy.priors.inverse_gamma(1,18) # 1, 25
     # m.set_prior('.*lengthscale',len_prior)
 
-    m.optimize(optimizer='scg', max_iters=optim_iters, messages=1)
+    m.optimize(optimizer='scg', max_iters=max_iters, messages=1)
 
     m.kern.plot_ARD()
     print(m)
     return m
 
-def silhouette(optim_iters=100):
+def silhouette(max_iters=100):
     """Predict the pose of a figure given a silhouette. This is a task from Agarwal and Triggs 2004 ICML paper."""
     data = GPy.util.datasets.silhouette()
 
@@ -138,12 +158,12 @@ def silhouette(optim_iters=100):
     m = GPy.models.GPRegression(data['X'], data['Y'])
 
     # optimize
-    m.optimize(messages=True, max_f_eval=optim_iters)
+    m.optimize(messages=True, max_iters=max_iters)
 
     print(m)
     return m
 
-def coregionalisation_toy2(optim_iters=100):
+def coregionalisation_toy2(max_iters=100):
     """
     A simple demonstration of coregionalisation on two sinusoidal functions.
     """
@@ -161,7 +181,7 @@ def coregionalisation_toy2(optim_iters=100):
     m = GPy.models.GPRegression(X, Y, kernel=k)
     m.constrain_fixed('.*rbf_var', 1.)
     # m.constrain_positive('.*kappa')
-    m.optimize('sim', messages=1, max_f_eval=optim_iters)
+    m.optimize('sim', messages=1, max_iters=max_iters)
 
     pb.figure()
     Xtest1 = np.hstack((np.linspace(0, 9, 100)[:, None], np.zeros((100, 1))))
@@ -174,7 +194,7 @@ def coregionalisation_toy2(optim_iters=100):
     pb.plot(X2[:, 0], Y2[:, 0], 'gx', mew=2)
     return m
 
-def coregionalisation_toy(optim_iters=100):
+def coregionalisation_toy(max_iters=100):
     """
     A simple demonstration of coregionalisation on two sinusoidal functions.
     """
@@ -192,7 +212,7 @@ def coregionalisation_toy(optim_iters=100):
     m = GPy.models.GPRegression(X, Y, kernel=k)
     m.constrain_fixed('.*rbf_var', 1.)
     # m.constrain_positive('kappa')
-    m.optimize(max_f_eval=optim_iters)
+    m.optimize(max_iters=max_iters)
 
     pb.figure()
     Xtest1 = np.hstack((np.linspace(0, 9, 100)[:, None], np.zeros((100, 1))))
@@ -206,7 +226,7 @@ def coregionalisation_toy(optim_iters=100):
     return m
 
 
-def coregionalisation_sparse(optim_iters=100):
+def coregionalisation_sparse(max_iters=100):
     """
     A simple demonstration of coregionalisation on two sinusoidal functions using sparse approximations.
     """
@@ -229,8 +249,8 @@ def coregionalisation_sparse(optim_iters=100):
     m.constrain_fixed('.*rbf_var', 1.)
     m.constrain_fixed('iip')
     m.constrain_bounded('noise_variance', 1e-3, 1e-1)
-#     m.optimize_restarts(5, robust=True, messages=1, max_iters=optim_iters, optimizer='bfgs')
-    m.optimize('bfgs', messages=1, max_iters=optim_iters)
+#     m.optimize_restarts(5, robust=True, messages=1, max_iters=max_iters, optimizer='bfgs')
+    m.optimize('bfgs', messages=1, max_iters=max_iters)
 
     # plotting:
     pb.figure()
@@ -248,7 +268,7 @@ def coregionalisation_sparse(optim_iters=100):
     return m
 
 
-def multiple_optima(gene_number=937, resolution=80, model_restarts=10, seed=10000, optim_iters=300):
+def multiple_optima(gene_number=937, resolution=80, model_restarts=10, seed=10000, max_iters=300):
     """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.
@@ -285,7 +305,7 @@ def multiple_optima(gene_number=937, resolution=80, model_restarts=10, seed=1000
         optim_point_y[0] = np.log10(m['rbf_variance']) - np.log10(m['noise_variance']);
 
         # optimize
-        m.optimize('scg', xtol=1e-6, ftol=1e-6, max_f_eval=optim_iters)
+        m.optimize('scg', xtol=1e-6, ftol=1e-6, max_iters=max_iters)
 
         optim_point_x[1] = m['rbf_lengthscale']
         optim_point_y[1] = np.log10(m['rbf_variance']) - np.log10(m['noise_variance']);
@@ -325,11 +345,32 @@ def _contour_data(data, length_scales, log_SNRs, kernel_call=GPy.kern.rbf):
 
     return np.array(lls)
 
-def sparse_GP_regression_1D(N=400, num_inducing=5, optim_iters=100):
+def robot_wireless(max_iters=100, kernel=None):
+    """Predict the location of a robot given wirelss signal strengthr readings."""
+    data = GPy.util.datasets.robot_wireless()
+
+    # create simple GP Model
+    m = GPy.models.GPRegression(data['Y'], data['X'], kernel=kernel)
+
+    # optimize
+    m.optimize(messages=True, max_iters=max_iters)
+    Xpredict = m.predict(data['Ytest'])[0]
+    pb.plot(data['Xtest'][:, 0], data['Xtest'][:, 1], 'r-')
+    pb.plot(Xpredict[:, 0], Xpredict[:, 1], 'b-')
+    pb.axis('equal')
+    pb.title('WiFi Localization with Gaussian Processes')
+    pb.legend(('True Location', 'Predicted Location'))
+
+    sse = ((data['Xtest'] - Xpredict)**2).sum()
+    print(m)
+    print('Sum of squares error on test data: ' + str(sse))
+    return m
+
+def sparse_GP_regression_1D(num_samples=400, num_inducing=5, max_iters=100):
     """Run a 1D example of a sparse GP regression."""
     # sample inputs and outputs
-    X = np.random.uniform(-3., 3., (N, 1))
-    Y = np.sin(X) + np.random.randn(N, 1) * 0.05
+    X = np.random.uniform(-3., 3., (num_samples, 1))
+    Y = np.sin(X) + np.random.randn(num_samples, 1) * 0.05
     # construct kernel
     rbf = GPy.kern.rbf(1)
     # create simple GP Model
@@ -337,14 +378,14 @@ def sparse_GP_regression_1D(N=400, num_inducing=5, optim_iters=100):
 
 
     m.checkgrad(verbose=1)
-    m.optimize('tnc', messages=1, max_f_eval=optim_iters)
+    m.optimize('tnc', messages=1, max_iters=max_iters)
     m.plot()
     return m
 
-def sparse_GP_regression_2D(N=400, num_inducing=50, optim_iters=100):
+def sparse_GP_regression_2D(num_samples=400, num_inducing=50, max_iters=100):
     """Run a 2D example of a sparse GP regression."""
-    X = np.random.uniform(-3., 3., (N, 2))
-    Y = np.sin(X[:, 0:1]) * np.sin(X[:, 1:2]) + np.random.randn(N, 1) * 0.05
+    X = np.random.uniform(-3., 3., (num_samples, 2))
+    Y = np.sin(X[:, 0:1]) * np.sin(X[:, 1:2]) + np.random.randn(num_samples, 1) * 0.05
 
     # construct kernel
     rbf = GPy.kern.rbf(2)
@@ -358,12 +399,12 @@ def sparse_GP_regression_2D(N=400, num_inducing=50, optim_iters=100):
     m.checkgrad()
 
     # optimize and plot
-    m.optimize('tnc', messages=1, max_f_eval=optim_iters)
+    m.optimize('tnc', messages=1, max_iters=max_iters)
     m.plot()
     print(m)
     return m
 
-def uncertain_inputs_sparse_regression(optim_iters=100):
+def uncertain_inputs_sparse_regression(max_iters=100):
     """Run a 1D example of a sparse GP regression with uncertain inputs."""
     fig, axes = pb.subplots(1, 2, figsize=(12, 5))
 
@@ -378,14 +419,14 @@ def uncertain_inputs_sparse_regression(optim_iters=100):
 
     # create simple GP Model - no input uncertainty on this one
     m = GPy.models.SparseGPRegression(X, Y, kernel=k, Z=Z)
-    m.optimize('scg', messages=1, max_f_eval=optim_iters)
+    m.optimize('scg', messages=1, max_iters=max_iters)
     m.plot(ax=axes[0])
     axes[0].set_title('no input uncertainty')
 
 
     # the same Model with uncertainty
     m = GPy.models.SparseGPRegression(X, Y, kernel=k, Z=Z, X_variance=S)
-    m.optimize('scg', messages=1, max_f_eval=optim_iters)
+    m.optimize('scg', messages=1, max_iters=max_iters)
     m.plot(ax=axes[1])
     axes[1].set_title('with input uncertainty')
     print(m)

GPy/kern/constructors.py

@@ -51,6 +51,44 @@ def linear(input_dim,variances=None,ARD=False):
     part = parts.linear.Linear(input_dim,variances,ARD)
     return kern(input_dim, [part])
 
+def mlp(input_dim,variance=1., weight_variance=None,bias_variance=100.,ARD=False):
+    """
+    Construct an MLP kernel
+
+    :param input_dim: dimensionality of the kernel, obligatory
+    :type input_dim: int
+    :param variance: the variance of the kernel
+    :type variance: float
+    :param weight_scale: the lengthscale of the kernel
+    :type weight_scale: vector of weight variances for input weights in neural network (length 1 if kernel is isotropic)
+    :param bias_variance: the variance of the biases in the neural network.
+    :type bias_variance: float
+    :param ARD: Auto Relevance Determination (allows for ARD version of covariance)
+    :type ARD: Boolean
+    """
+    part = parts.mlp.MLP(input_dim,variance,weight_variance,bias_variance,ARD)
+    return kern(input_dim, [part])
+
+def poly(input_dim,variance=1., weight_variance=None,bias_variance=1.,degree=2, ARD=False):
+    """
+    Construct a polynomial kernel
+
+    :param input_dim: dimensionality of the kernel, obligatory
+    :type input_dim: int
+    :param variance: the variance of the kernel
+    :type variance: float
+    :param weight_scale: the lengthscale of the kernel
+    :type weight_scale: vector of weight variances for input weights.
+    :param bias_variance: the variance of the biases.
+    :type bias_variance: float
+    :param degree: the degree of the polynomial
+    :type degree: int
+    :param ARD: Auto Relevance Determination (allows for ARD version of covariance)
+    :type ARD: Boolean
+    """
+    part = parts.poly.POLY(input_dim,variance,weight_variance,bias_variance,degree,ARD)
+    return kern(input_dim, [part])
+
 def white(input_dim,variance=1.):
     """
      Construct a white kernel.
@@ -253,6 +291,8 @@ def prod(k1,k2,tensor=False):
 
     :param k1, k2: the kernels to multiply
     :type k1, k2: kernpart
+    :param tensor: The kernels are either multiply as functions defined on the same input space (default) or on the product of the input spaces
+    :type tensor: Boolean
     :rtype: kernel object
     """
     part = parts.prod.Prod(k1, k2, tensor)
@@ -260,7 +300,7 @@ def prod(k1,k2,tensor=False):
 
 def symmetric(k):
     """
-    Construct a symmetrical kernel from an existing kernel
+    Construct a symmetric kernel from an existing kernel
     """
     k_ = k.copy()
     k_.parts = [symmetric.Symmetric(p) for p in k.parts]
@@ -291,11 +331,13 @@ def fixed(input_dim, K, variance=1.):
     """
      Construct a Fixed effect kernel.
 
-     Arguments
-     ---------
-     input_dim (int), obligatory
-     K (np.array), obligatory
-     variance (float)
+    :param input_dim: the number of input dimensions
+    :type input_dim: int (input_dim=1 is the only value currently supported)
+    :param K: the variance :math:`\sigma^2`
+    :type K: np.array
+    :param variance: kernel variance
+    :type variance: float
+    :rtype: kern object
     """
     part = parts.fixed.Fixed(input_dim, K, variance)
     return kern(input_dim, [part])
@@ -318,7 +360,7 @@ def independent_outputs(k):
 
 def hierarchical(k):
     """
-    Construct a kernel with independent outputs from an existing kernel
+    TODO THis can't be right! Construct a kernel with independent outputs from an existing kernel
     """
     # for sl in k.input_slices:
     #     assert (sl.start is None) and (sl.stop is None), "cannot adjust input slices! (TODO)"

GPy/kern/kern.py

@@ -14,7 +14,10 @@ class kern(Parameterized):
         """
         This is the main kernel class for GPy. It handles multiple (additive) kernel functions, and keeps track of variaous things like which parameters live where.
 
-        The technical code for kernels is divided into _parts_ (see e.g. rbf.py). This obnject contains a list of parts, which are computed additively. For multiplication, special _prod_ parts are used.
+        The technical code for kernels is divided into _parts_ (see
+        e.g. rbf.py). This object contains a list of parts, which are
+        computed additively. For multiplication, special _prod_ parts
+        are used.
 
         :param input_dim: The dimensionality of the kernel's input space
         :type input_dim: int
@@ -149,7 +152,7 @@ class kern(Parameterized):
             return g
 
     def compute_param_slices(self):
-        """create a set of slices that can index the parameters of each part"""
+        """create a set of slices that can index the parameters of each part."""
         self.param_slices = []
         count = 0
         for p in self.parts:
@@ -200,11 +203,19 @@ class kern(Parameterized):
         """
         return self.prod(other)
 
+    def __pow__(self, other, tensor=False):
+        """
+        Shortcut for tensor `prod`.
+        """
+        return self.prod(other, tensor=True)
+
     def prod(self, other, tensor=False):
         """
-        multiply two kernels (either on the same space, or on the tensor product of the input space)
+        multiply two kernels (either on the same space, or on the tensor product of the input space).
         :param other: the other kernel to be added
         :type other: GPy.kern
+        :param tensor: whether or not to use the tensor space (default is false).
+        :type tensor: bool 
         """
         K1 = self.copy()
         K2 = other.copy()
@@ -273,7 +284,7 @@ class kern(Parameterized):
         [p._set_params(x[s]) for p, s in zip(self.parts, self.param_slices)]
 
     def _get_param_names(self):
-        # this is a bit nasty: we wat to distinguish between parts with the same name by appending a count
+        # this is a bit nasty: we want 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)]
@@ -295,11 +306,13 @@ class kern(Parameterized):
 
     def dK_dtheta(self, dL_dK, X, X2=None):
         """
-        :param dL_dK: An array of dL_dK derivaties, dL_dK
-        :type dL_dK: Np.ndarray (N x num_inducing)
+        Compute the gradient of the covariance function with respect to the parameters.
+        
+        :param dL_dK: An array of gradients of the objective function with respect to the covariance function.
+        :type dL_dK: Np.ndarray (num_samples x num_inducing)
         :param X: Observed data inputs
-        :type X: np.ndarray (N x input_dim)
-        :param X2: Observed dara inputs (optional, defaults to X)
+        :type X: np.ndarray (num_samples x input_dim)
+        :param X2: Observed data inputs (optional, defaults to X)
         :type X2: np.ndarray (num_inducing x input_dim)
         """
         assert X.shape[1] == self.input_dim
@@ -312,6 +325,14 @@ class kern(Parameterized):
         return self._transform_gradients(target)
 
     def dK_dX(self, dL_dK, X, X2=None):
+        """Compute the gradient of the covariance function with respect to X.
+
+        :param dL_dK: An array of gradients of the objective function with respect to the covariance function.
+        :type dL_dK: np.ndarray (num_samples x num_inducing)
+        :param X: Observed data inputs
+        :type X: np.ndarray (num_samples x input_dim)
+        :param X2: Observed data inputs (optional, defaults to X)
+        :type X2: np.ndarray (num_inducing x input_dim)"""
         if X2 is None:
             X2 = X
         target = np.zeros_like(X)
@@ -322,6 +343,7 @@ class kern(Parameterized):
         return target
 
     def Kdiag(self, X, which_parts='all'):
+        """Compute the diagonal of the covariance function for inputs X."""
         if which_parts == 'all':
             which_parts = [True] * self.Nparts
         assert X.shape[1] == self.input_dim
@@ -330,6 +352,7 @@ class kern(Parameterized):
         return target
 
     def dKdiag_dtheta(self, dL_dKdiag, X):
+        """Compute the gradient of the diagonal of the covariance function with respect to the parameters."""
         assert X.shape[1] == self.input_dim
         assert dL_dKdiag.size == X.shape[0]
         target = np.zeros(self.num_params)
@@ -373,16 +396,18 @@ class kern(Parameterized):
         return target
 
     def dpsi1_dmuS(self, dL_dpsi1, Z, mu, S):
-        """return shapes are N,num_inducing,input_dim"""
+        """return shapes are num_samples,num_inducing,input_dim"""
         target_mu, target_S = np.zeros((2, mu.shape[0], mu.shape[1]))
         [p.dpsi1_dmuS(dL_dpsi1, Z[:, i_s], mu[:, i_s], S[:, i_s], target_mu[:, i_s], target_S[:, i_s]) for p, i_s in zip(self.parts, self.input_slices)]
         return target_mu, target_S
 
     def psi2(self, Z, mu, S):
         """
+        Computer the psi2 statistics for the covariance function.
+        
         :param Z: np.ndarray of inducing inputs (num_inducing x input_dim)
-        :param mu, S: np.ndarrays of means and variances (each N x input_dim)
-        :returns psi2: np.ndarray (N,num_inducing,num_inducing)
+        :param mu, S: np.ndarrays of means and variances (each num_samples x input_dim)
+        :returns psi2: np.ndarray (num_samples,num_inducing,num_inducing)
         """
         target = np.zeros((mu.shape[0], Z.shape[0], Z.shape[0]))
         [p.psi2(Z[:, i_s], mu[:, i_s], S[:, i_s], target) for p, i_s in zip(self.parts, self.input_slices)]
@@ -406,6 +431,7 @@ class kern(Parameterized):
         return target
 
     def dpsi2_dtheta(self, dL_dpsi2, Z, mu, S):
+        """Gradient of the psi2 statistics with respect to the parameters."""
         target = np.zeros(self.num_params)
         [p.dpsi2_dtheta(dL_dpsi2, Z[:, i_s], mu[:, i_s], S[:, i_s], target[ps]) for p, i_s, ps in zip(self.parts, self.input_slices, self.param_slices)]
 
@@ -509,3 +535,4 @@ class kern(Parameterized):
             pb.title("k(x1,x2 ; %0.1f,%0.1f)" % (x[0, 0], x[0, 1]))
         else:
             raise NotImplementedError, "Cannot plot a kernel with more than two input dimensions"
+

GPy/kern/parts/__init__.py

@@ -8,9 +8,11 @@ import independent_outputs
 import linear
 import Matern32
 import Matern52
+import mlp
 import periodic_exponential
 import periodic_Matern32
 import periodic_Matern52
+import poly
 import prod_orthogonal
 import prod
 import rational_quadratic

GPy/kern/parts/kernpart.py

@@ -29,7 +29,11 @@ class Kernpart(object):
     def dK_dtheta(self,dL_dK,X,X2,target):
         raise NotImplementedError
     def dKdiag_dtheta(self,dL_dKdiag,X,target):
-        raise NotImplementedError
+        # In the base case compute this by calling dK_dtheta. Need to
+        # override for stationary covariances (for example) to save
+        # time.
+        for i in range(X.shape[0]):
+            self.dK_dtheta(dL_dKdiag[i], X[i, :][None, :], X2=None, target=target)
     def psi0(self,Z,mu,S,target):
         raise NotImplementedError
     def dpsi0_dtheta(self,dL_dpsi0,Z,mu,S,target):
@@ -52,5 +56,21 @@ class Kernpart(object):
         raise NotImplementedError
     def dpsi2_dmuS(self,dL_dpsi2,Z,mu,S,target_mu,target_S):
         raise NotImplementedError
-    def dK_dX(self,X,X2,target):
+    def dK_dX(self, dL_dK, X, X2, target):
         raise NotImplementedError
+
+class Kernpart_inner(Kernpart):
+    def __init__(self,input_dim):
+        """
+        The base class for a kernpart_inner: a positive definite function which forms part of a kernel that is based on the inner product between inputs.
+
+        :param input_dim: the number of input dimensions to the function
+        :type input_dim: int
+
+        Do not instantiate.
+        """
+        Kernpart.__init__(self, input_dim)
+
+        # initialize cache
+        self._Z, self._mu, self._S = np.empty(shape=(3, 1))
+        self._X, self._X2, self._params = np.empty(shape=(3, 1))

GPy/kern/parts/mlp.py

@@ -0,0 +1,164 @@
+# Copyright (c) 2013, GPy authors (see AUTHORS.txt).
+# Licensed under the BSD 3-clause license (see LICENSE.txt)
+
+from kernpart import Kernpart
+import numpy as np
+four_over_tau = 2./np.pi
+
+class MLP(Kernpart):
+    """
+    multi layer perceptron kernel (also known as arc sine kernel or neural network kernel)
+
+    .. math::
+
+       k(x,y) = \sigma^2 \frac{2}{\pi}  \text{asin} \left(\frac{\sigma_w^2 x^\top y+\sigma_b^2}{\sqrt{\sigma_w^2x^\top x + \sigma_b^2 + 1}\sqrt{\sigma_w^2 y^\top y \sigma_b^2 +1}} \right)
+
+    :param input_dim: the number of input dimensions
+    :type input_dim: int 
+    :param variance: the variance :math:`\sigma^2`
+    :type variance: float
+    :param weight_variance: the vector of the variances of the prior over input weights in the neural network :math:`\sigma^2_w`
+    :type weight_variance: array or list of the appropriate size (or float if there is only one weight variance parameter)
+    :param bias_variance: the variance of the prior over bias parameters :math:`\sigma^2_b`
+    :param ARD: Auto Relevance Determination. If equal to "False", the kernel is isotropic (ie. one weight variance parameter \sigma^2_w), otherwise there is one weight variance parameter per dimension.
+    :type ARD: Boolean
+    :rtype: Kernpart object
+
+    """
+
+    def __init__(self, input_dim, variance=1., weight_variance=None, bias_variance=100., ARD=False):
+        ARD = False
+        self.input_dim = input_dim
+        self.ARD = ARD
+        if not ARD:
+            self.num_params=3
+            if weight_variance is not None:
+                weight_variance = np.asarray(weight_variance)
+                assert weight_variance.size == 1, "Only one weight variance needed for non-ARD kernel"
+            else:
+                weight_variance = 100.*np.ones(1)
+        else:
+            self.num_params = self.input_dim + 2
+            if weight_variance is not None:
+                weight_variance = np.asarray(weight_variance)
+                assert weight_variance.size == self.input_dim, "bad number of weight variances"
+            else:
+                weight_variance = np.ones(self.input_dim)
+
+        self.name='mlp'
+        self._set_params(np.hstack((variance, weight_variance.flatten(), bias_variance)))
+
+    def _get_params(self):
+        return np.hstack((self.variance, self.weight_variance.flatten(), self.bias_variance))
+
+    def _set_params(self, x):
+        assert x.size == (self.num_params)
+        self.variance = x[0]
+        self.weight_variance = x[1:-1]
+        self.weight_std = np.sqrt(self.weight_variance)
+        self.bias_variance = x[-1]
+
+    def _get_param_names(self):
+        if self.num_params == 3:
+            return ['variance', 'weight_variance', 'bias_variance']
+        else:
+            return ['variance'] + ['weight_variance_%i' % i for i in range(self.lengthscale.size)] + ['bias_variance']
+
+    def K(self, X, X2, target):
+        """Return covariance between X and X2."""
+        self._K_computations(X, X2)
+        target += self.variance*self._K_dvar
+
+    def Kdiag(self, X, target):
+        """Compute the diagonal of the covariance matrix for X."""
+        self._K_diag_computations(X)
+        target+= self.variance*self._K_diag_dvar
+
+    def dK_dtheta(self, dL_dK, X, X2, target):
+        """Derivative of the covariance with respect to the parameters."""
+        self._K_computations(X, X2)
+        denom3 = self._K_denom*self._K_denom*self._K_denom
+        base = four_over_tau*self.variance/np.sqrt(1-self._K_asin_arg*self._K_asin_arg)
+        base_cov_grad = base*dL_dK
+
+        if X2 is None:
+            vec = np.diag(self._K_inner_prod)
+            target[1] += ((self._K_inner_prod/self._K_denom 
+                           -.5*self._K_numer/denom3
+                           *(np.outer((self.weight_variance*vec+self.bias_variance+1.), vec) 
+                             +np.outer(vec,(self.weight_variance*vec+self.bias_variance+1.))))*base_cov_grad).sum()
+            target[2] += ((1./self._K_denom 
+                           -.5*self._K_numer/denom3 
+                           *((vec[None, :]+vec[:, None])*self.weight_variance
+                           +2.*self.bias_variance + 2.))*base_cov_grad).sum()
+        else:
+            vec1 = (X*X).sum(1)
+            vec2 = (X2*X2).sum(1)
+            target[1] += ((self._K_inner_prod/self._K_denom 
+                           -.5*self._K_numer/denom3
+                           *(np.outer((self.weight_variance*vec1+self.bias_variance+1.), vec2) + np.outer(vec1, self.weight_variance*vec2 + self.bias_variance+1.)))*base_cov_grad).sum()
+            target[2] += ((1./self._K_denom 
+                           -.5*self._K_numer/denom3 
+                           *((vec1[:, None]+vec2[None, :])*self.weight_variance
+                             + 2*self.bias_variance + 2.))*base_cov_grad).sum()
+            
+        target[0] += np.sum(self._K_dvar*dL_dK)
+
+
+    def dK_dX(self, dL_dK, X, X2, target):
+        """Derivative of the covariance matrix with respect to X"""
+        raise NotImplementedError
+        # self._K_computations(X, X2)
+        # gX = np.zeros((X2.shape[0], X.shape[1], X.shape[0]))
+        
+        # for i in range(X.shape[0]):
+        #     gX[:, :, i] = self._dK_dX_point(dL_dK, X, X2, target, i)
+            
+            
+    def _dK_dX_point(self, dL_dK, X, X2, target, i):
+        """Gradient with respect to one point of X"""
+        
+        inner_prod = self._K_inner_prod[i, :].T
+        numer = self._K_numer[i, :].T
+        denom = self._K_denom[i, :].T
+        arg = self._K_asin_arg[i, :].T
+        vec1 = (X[i, :]*X[i, :]).sum()*self.weight_variance + self.bias_variance + 1.
+        vec2 = (X2*X2).sum(1)*self.weight_variance + self.bias_variance + 1.
+        #denom = np.sqrt(np.outer(vec2,vec1))
+        #arg = numer/denom
+        gX = np.zeros(X2.shape)
+        denom3 = denom*denom*denom
+        gX = np.zeros((X2.shape[0], X2.shape[1]))
+        for j in range(X2.shape[1]):
+            gX[:, j] =X2[:, j]/denom - vec2*X[i, j]*numer/denom3
+            gX[:, j] = four_over_tau*self.weight_variance*self.variance*gX[:, j]/np.sqrt(1-arg*arg)
+            target[i, :]
+    
+    
+    def _K_computations(self, X, X2):
+        if self.ARD:
+            pass
+        else:
+            if X2 is None:
+                self._K_inner_prod = np.dot(X,X.T)
+                self._K_numer = self._K_inner_prod*self.weight_variance+self.bias_variance
+                vec = np.diag(self._K_numer) + 1.
+                self._K_denom = np.sqrt(np.outer(vec,vec))
+                self._K_asin_arg = self._K_numer/self._K_denom
+                self._K_dvar = four_over_tau*np.arcsin(self._K_asin_arg)
+            else:
+                self._K_inner_prod = np.dot(X,X2.T)
+                self._K_numer = self._K_inner_prod*self.weight_variance + self.bias_variance
+                vec1 = (X*X).sum(1)*self.weight_variance + self.bias_variance + 1.
+                vec2 = (X2*X2).sum(1)*self.weight_variance + self.bias_variance + 1.
+                self._K_denom = np.sqrt(np.outer(vec1,vec2))
+                self._K_asin_arg = self._K_numer/self._K_denom
+                self._K_dvar = four_over_tau*np.arcsin(self._K_asin_arg)
+
+    def _K_diag_computations(self, X):
+        if self.ARD:
+            pass
+        else:
+            self._K_diag_numer = (X*X).sum(1)*self.weight_variance + self.bias_variance
+            self._K_diag_denom = self._K_diag_numer+1.
+            self._K_diag_dvar = four_over_tau*np.arcsin(self._K_diag_numer/self._K_diag_denom)

GPy/models/gp_regression.py

@@ -37,7 +37,6 @@ class GPRegression(GP):
     def getstate(self):
         return GP.getstate(self)
 
-
     def setstate(self, state):
         return GP.setstate(self, state)
 

GPy/notes.txt

@@ -0,0 +1,31 @@
+Prod.py kernel should take a list of kernels rather than two arguments for kernels.
+transformations.py should have limits on what is fed into exp() particularly for the negative log logistic. 
+
+Load in a model with mlp kernel, plot it, change a parameter, plot it again. It doesn't update the plot.
+
+Tests for kernels which work directly on the kernel implementation (not through GP).
+
+Should stationary covariances have their own kernpart type?
+
+Where do we declare default kernel parameters. In constructors.py or in the definition file for the kernel?
+
+When printing to stdout, can we check that our approach is also working nicely for the ipython notebook? I like the way our optimization ticks over, but at the moment this doesn't seem to work in the ipython notebook, it would be nice if it did.
+
+When we print a model should we also include information such as number of inputs and number of outputs?
+
+Let's not use N for giving the number of data in the model. When it pops up as a help tip it's not as clear as num_samples or num_data. Prefer the second, but oddly I've been using first.
+
+Loving the fact that the * has been overloaded on the kernels (oddly never thought to check this before). Although naming can be a bit confusing. Can we think how to deal with the names in a clearer way when we use a kernel like this one:
+kern = GPy.kern.rbf(30)*(GPy.kern.mlp(30)+GPy.kern.poly(30, degree=5)) + GPy.kern.bias(30). There seems to be some tieing of parameters going on ... should there be? (you can try it as the kernel for the robot wireless model).
+
+Can we comment up some of the list incomprehensions in hierarchical.py??
+
+Need to tidy up classification.py, 
+many examples include help that doesn't apply 
+(it is suggested that you can try different approximation types)
+
+Shall we overload the ** operator to have tensor products? (I've done this now we can see if we like it)
+
+People aren't filling the doc strings in as they go *everyone* needs to get in the habit of this (and modifying them as they edit, or correcting them when there is a problem).
+
+Need some nice way of explaining how to compile documentation and run the unit tests, could this be in a readme or FAQ somewhere? Maybe it's there already somewhere and I've missed it.

GPy/testing/kernel_tests.py

@@ -1,20 +1,118 @@
-# Copyright (c) 2012, GPy authors (see AUTHORS.txt).
+# Copyright (c) 2012, 2013 GPy authors (see AUTHORS.txt).
 # Licensed under the BSD 3-clause license (see LICENSE.txt)
 
 import unittest
 import numpy as np
 import GPy
+from GPy.core.model import Model
+
+class Kern_check_model(Model):
+    """This is a dummy model class used as a base class for checking that the gradients of a given kernel are implemented correctly. It enables checkgradient() to be called independently on a kernel."""
+    def __init__(self, kernel=None, dL_dK=None, X=None, X2=None):
+        num_samples = 20
+        num_samples2 = 10
+        if kernel==None:
+            kernel = GPy.kern.rbf(1)
+        if X==None:
+            X = np.random.randn(num_samples, kernel.input_dim)
+        if X2==None:
+            X2 = np.random.randn(num_samples2, kernel.input_dim)
+        if dL_dK==None:
+            dL_dK = np.ones((X.shape[0], X2.shape[0]))
+        
+        self.kernel=kernel
+        self.X = X
+        self.X2 = X2
+        self.dL_dK = dL_dK
+        #self.constrained_indices=[]
+        #self.constraints=[]
+        Model.__init__(self)
+
+    def is_positive_definite(self):
+        v = np.linalg.eig(self.kernel.K(self.X))[0]
+        if any(v<0):
+            return False
+        else:
+            return True
+        
+    def _get_params(self):
+        return self.kernel._get_params()
+
+    def _get_param_names(self):
+        return self.kernel._get_param_names()
+
+    def _set_params(self, x):
+        self.kernel._set_params(x)
+
+    def log_likelihood(self):
+        return (self.dL_dK*self.kernel.K(self.X, self.X2)).sum()
+
+    def _log_likelihood_gradients(self):
+        raise NotImplementedError, "This needs to be implemented to use the kern_check_model class."
+    
+class Kern_check_dK_dtheta(Kern_check_model):
+    """This class allows gradient checks for the gradient of a kernel with respect to parameters. """
+    def __init__(self, kernel=None, dL_dK=None, X=None, X2=None):
+        Kern_check_model.__init__(self,kernel=kernel,dL_dK=dL_dK, X=X, X2=X2)
+
+    def _log_likelihood_gradients(self):
+        return self.kernel.dK_dtheta(self.dL_dK, self.X, self.X2)
+
+class Kern_check_dKdiag_dtheta(Kern_check_model):
+    """This class allows gradient checks of the gradient of the diagonal of a kernel with respect to the parameters."""
+    def __init__(self, kernel=None, dL_dK=None, X=None):
+        Kern_check_model.__init__(self,kernel=kernel,dL_dK=dL_dK, X=X, X2=None)
+        if dL_dK==None:
+            self.dL_dK = np.ones((self.X.shape[0]))
+        
+    def log_likelihood(self):
+        return (self.dL_dK*self.kernel.Kdiag(self.X)).sum()
+
+    def _log_likelihood_gradients(self):
+        return self.kernel.dKdiag_dtheta(self.dL_dK, self.X)
+
+class Kern_check_dK_dX(Kern_check_model):
+    """This class allows gradient checks for the gradient of a kernel with respect to X. """
+    def __init__(self, kernel=None, dL_dK=None, X=None, X2=None):
+        Kern_check_model.__init__(self,kernel=kernel,dL_dK=dL_dK, X=X, X2=X2)
+
+    def _log_likelihood_gradients(self):
+        return self.kernel.dK_dX(self.dL_dK, self.X, self.X2).flatten()
+
+    def _get_param_names(self):
+        names = []
+        for i in range(self.X.shape[0]):
+            for j in range(self.X.shape[0]):
+                names.append('X_' +str(i) + ','+str(j))
+        return names
+                
+    def _get_params(self):
+        return self.X.flatten()
+
+    def _set_params(self, x):
+        self.X=x.reshape(self.X.shape)
+    
+
 
 class KernelTests(unittest.TestCase):
     def test_kerneltie(self):
         K = GPy.kern.rbf(5, ARD=True)
         K.tie_params('.*[01]')
         K.constrain_fixed('2')
+        
         X = np.random.rand(5,5)
         Y = np.ones((5,1))
         m = GPy.models.GPRegression(X,Y,K)
         self.assertTrue(m.checkgrad())
 
+    def test_rbfkernel(self):
+        verbose = False
+        kern = GPy.kern.rbf(5)
+        self.assertTrue(Kern_check_model(kern).is_positive_definite())
+        self.assertTrue(Kern_check_dK_dtheta(kern).checkgrad(verbose=verbose))
+        self.assertTrue(Kern_check_dKdiag_dtheta(kern).checkgrad(verbose=verbose))
+        self.assertTrue(Kern_check_dK_dX(kern).checkgrad(verbose=verbose))
+
     def test_fixedkernel(self):
         """
         Fixed effect kernel test
@@ -43,7 +141,6 @@ class KernelTests(unittest.TestCase):
 
 
        
-
 if __name__ == "__main__":
     print "Running unit tests, please be (very) patient..."
     unittest.main()

GPy/testing/unit_tests.py

@@ -59,6 +59,16 @@ class GradientTests(unittest.TestCase):
         k = GPy.kern.rbf(2, ARD=True)
         self.check_model_with_white(k, model_type='GPRegression', dimension=2)
 
+    def test_GPRegression_mlp_1d(self):
+        ''' Testing the GP regression with mlp kernel with white kernel on 1d data '''
+        mlp = GPy.kern.mlp(1)
+        self.check_model_with_white(mlp, model_type='GPRegression', dimension=1)
+
+    def test_GPRegression_poly_1d(self):
+        ''' Testing the GP regression with polynomial kernel with white kernel on 1d data '''
+        mlp = GPy.kern.poly(1, degree=5)
+        self.check_model_with_white(mlp, model_type='GPRegression', dimension=1)
+
     def test_GPRegression_matern52_1D(self):
         ''' Testing the GP regression with matern52 kernel on 1d data '''
         matern52 = GPy.kern.Matern52(1)

GPy/util/datasets.py

@@ -5,13 +5,131 @@ import GPy
 import scipy.sparse
 import scipy.io
 import cPickle as pickle
-import urllib2 as url
+import urllib as url
+import zipfile
+import tarfile
+import gzip
+import zlib
 
+import sys, urllib
+def reporthook(a,b,c): 
+    # ',' at the end of the line is important!
+    print "% 3.1f%% of %d bytes\r" % (min(100, float(a * b) / c * 100), c),
+    #you can also use sys.stdout.write
+    #sys.stdout.write("\r% 3.1f%% of %d bytes" 
+    #                 % (min(100, float(a * b) / c * 100), c)
+    sys.stdout.flush()
+     
+# Global variables
 data_path = os.path.join(os.path.dirname(__file__), 'datasets')
 default_seed = 10000
-neil_url = 'http://staffwww.dcs.shef.ac.uk/people/N.Lawrence/'
+overide_manual_authorize=False
+neil_url = 'http://staffwww.dcs.shef.ac.uk/people/N.Lawrence/dataset_mirror/'
+cmu_url = 'http://mocap.cs.cmu.edu/subjects/'
+# Note: there may be a better way of storing data resources. One of the pythonistas will need to take a look.
+data_resources = {'ankur_pose_data' : {'urls' : [neil_url + 'ankur_pose_data/'],
+                                       'files' : [['ankurDataPoseSilhouette.mat']],
+                                       'license' : None,
+                                       'citation' : """3D Human Pose from Silhouettes by Relevance Vector Regression (In CVPR'04). A. Agarwal and B. Triggs.""",
+                                       'details' : """Artificially generated data of silhouettes given poses. Note that the data does not display a left/right ambiguity because across the entire data set one of the arms sticks out more the the other, disambiguating the pose as to which way the individual is facing."""},
+                   
+                  'boston_housing' : {'urls' : ['http://archive.ics.uci.edu/ml/machine-learning-databases/housing/'],
+                                      'files' : [['Index', 'housing.data', 'housing.names']],
+                                      'citation' : """Harrison, D. and Rubinfeld, D.L. 'Hedonic prices and the demand for clean air', J. Environ. Economics & Management, vol.5, 81-102, 1978.""",
+                                      'details' : """The Boston Housing data relates house values in Boston to a range of input variables.""",
+                                      'license' : None,
+                                      'size' : 51276
+                                      },
+                  'brendan_faces' : {'urls' : ['http://www.cs.nyu.edu/~roweis/data/'],
+                                     'files': [['frey_rawface.mat']],
+                                     'citation' : 'Frey, B. J., Colmenarez, A and Huang, T. S. Mixtures of Local Linear Subspaces for Face Recognition. Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition 1998, 32-37, June 1998. Computer Society Press, Los Alamitos, CA.',
+                                     'details' : """A video of Brendan Frey's face popularized as a benchmark for visualization by the Locally Linear Embedding.""",
+                                     'license': None,
+                                     'size' : 1100584},
+                  'cmu_mocap_full' : {'urls' : ['http://mocap.cs.cmu.edu'],
+                                 'files' : [['allasfamc.zip']],
+                                 'citation' : """Please include this in your acknowledgements: The data used in this project was obtained from mocap.cs.cmu.edu.
+The database was created with funding from NSF EIA-0196217.""",
+                                 'details' : """CMU Motion Capture data base. Captured by a Vicon motion capture system consisting of 12 infrared MX-40 cameras, each of which is capable of recording at 120 Hz with images of 4 megapixel resolution. Motions are captured in a working volume of approximately 3m x 8m. The capture subject wears 41 markers and a stylish black garment.""",
+                                 'license' : """From http://mocap.cs.cmu.edu. This data is free for use in research projects. You may include this data in commercially-sold products, but you may not resell this data directly, even in converted form. If you publish results obtained using this data, we would appreciate it if you would send the citation to your published paper to jkh+mocap@cs.cmu.edu, and also would add this text to your acknowledgments section: The data used in this project was obtained from mocap.cs.cmu.edu. The database was created with funding from NSF EIA-0196217.""",
+                                 'size' : None},
+                  'creep_rupture' : {'urls' : ['http://www.msm.cam.ac.uk/map/data/tar/'],
+                                     'files' : [['creeprupt.tar']],
+                                     'citation' : 'Materials Algorithms Project Data Library: MAP_DATA_CREEP_RUPTURE. F. Brun and T. Yoshida.',
+                                     'details' : """Provides 2066 creep rupture test results of steels (mainly of two kinds of steels: 2.25Cr and 9-12 wt% Cr ferritic steels). See http://www.msm.cam.ac.uk/map/data/materials/creeprupt-b.html.""",
+                                     'license' : None,
+                                     'size' : 602797},
+                  'della_gatta' : {'urls' : [neil_url + 'della_gatta/'],
+                                   'files': [['DellaGattadata.mat']],
+                                   'citation' : 'Direct targets of the TRP63 transcription factor revealed by a combination of gene expression profiling and reverse engineering. Giusy Della Gatta, Mukesh Bansal, Alberto Ambesi-Impiombato, Dario Antonini, Caterina Missero, and Diego di Bernardo, Genome Research 2008',
+                                   'details': "The full gene expression data set from della Gatta et al (http://www.ncbi.nlm.nih.gov/pmc/articles/PMC2413161/) processed by RMA.",
+                                   'license':None,
+                                   'size':3729650},
+                  'three_phase_oil_flow': {'urls' : [neil_url + 'three_phase_oil_flow/'],
+                                           'files' : [['DataTrnLbls.txt', 'DataTrn.txt', 'DataTst.txt', 'DataTstLbls.txt', 'DataVdn.txt', 'DataVdnLbls.txt']],
+                                           'citation' : 'Bishop, C. M. and G. D. James (1993). Analysis of multiphase flows using dual-energy gamma densitometry and neural networks. Nuclear Instruments and Methods in Physics Research A327, 580-593',
+                                           'details' : """The three phase oil data used initially for demonstrating the Generative Topographic mapping.""",
+                                           'license' : None,
+                                           'size' : 712796},
+                  'rogers_girolami_data' : {'urls' : ['https://www.dropbox.com/sh/7p6tu1t29idgliq/_XqlH_3nt9/'],
+                                            'files' : [['firstcoursemldata.tar.gz']],
+                                            'suffices' : [['?dl=1']],
+                                            'citation' : 'A First Course in Machine Learning. Simon Rogers and Mark Girolami: Chapman & Hall/CRC, ISBN-13: 978-1439824146',
+                                            'details' : """Data from the textbook 'A First Course in Machine Learning'. Available from http://www.dcs.gla.ac.uk/~srogers/firstcourseml/.""",
+                                            'license' : None,
+                                            'size' : 21949154},
+                  'olympic_marathon_men' : {'urls' : [neil_url + 'olympic_marathon_men/'],
+                                            'files' : [['olympicMarathonTimes.csv']],
+                                            'citation' : None,
+                                            'details' : """Olympic mens' marathon gold medal winning times from 1896 to 2012. Time given in pace (minutes per kilometer). Data is originally downloaded and collated from Wikipedia, we are not responsible for errors in the data""",
+                                            'license': None,
+                                            'size' : 584},
+                  'osu_run1' : {'urls': ['http://accad.osu.edu/research/mocap/data/', neil_url + 'stick/'],
+                                'files': [['run1TXT.ZIP'],['connections.txt']],
+                                'details' : "Motion capture data of a stick man running from the Open Motion Data Project at Ohio State University.",
+                                'citation' : 'The Open Motion Data Project by The Ohio State University Advanced Computing Center for the Arts and Design, http://accad.osu.edu/research/mocap/mocap_data.htm.',
+                                'license' : 'Data is licensed under a Creative Commons Attribution-NonCommercial-ShareAlike 3.0 Unported License (http://creativecommons.org/licenses/by-nc-sa/3.0/).',
+                                'size': 338103},
+                  'osu_accad' : {'urls': ['http://accad.osu.edu/research/mocap/data/', neil_url + 'stick/'],
+                                'files': [['swagger1TXT.ZIP','handspring1TXT.ZIP','quickwalkTXT.ZIP','run1TXT.ZIP','sprintTXT.ZIP','dogwalkTXT.ZIP','camper_04TXT.ZIP','dance_KB3_TXT.ZIP','per20_TXT.ZIP','perTWO07_TXT.ZIP','perTWO13_TXT.ZIP','perTWO14_TXT.ZIP','perTWO15_TXT.ZIP','perTWO16_TXT.ZIP'],['connections.txt']],
+                                'details' : "Motion capture data of different motions from the Open Motion Data Project at Ohio State University.",
+                                'citation' : 'The Open Motion Data Project by The Ohio State University Advanced Computing Center for the Arts and Design, http://accad.osu.edu/research/mocap/mocap_data.htm.',
+                                'license' : 'Data is licensed under a Creative Commons Attribution-NonCommercial-ShareAlike 3.0 Unported License (http://creativecommons.org/licenses/by-nc-sa/3.0/).',
+                                'size': 15922790},
+                  'pumadyn-32nm' : {'urls' : ['ftp://ftp.cs.toronto.edu/pub/neuron/delve/data/tarfiles/pumadyn-family/'],
+                                    'files' : [['pumadyn-32nm.tar.gz']],
+                                    'details' : """Pumadyn non linear 32 input data set with moderate noise. See http://www.cs.utoronto.ca/~delve/data/pumadyn/desc.html for details.""",
+                                    'citation' : """Created by Zoubin Ghahramani using the Matlab Robotics Toolbox of Peter Corke. Corke, P. I. (1996). A Robotics Toolbox for MATLAB. IEEE Robotics and Automation Magazine, 3 (1): 24-32.""",
+                                    'license' : """Data is made available by the Delve system at the University of Toronto""",
+                                    'size' : 5861646},
+                  'robot_wireless' : {'urls' : [neil_url + 'robot_wireless/'],
+                                      'files' : [['uw-floor.txt']],
+                                      'citation' : """WiFi-SLAM using Gaussian Process Latent Variable Models by Brian Ferris, Dieter Fox and Neil Lawrence in IJCAI'07 Proceedings pages 2480-2485. Data used in A Unifying Probabilistic Perspective for Spectral Dimensionality Reduction: Insights and New Models by Neil D. Lawrence, JMLR 13 pg 1609--1638, 2012.""",
+                                      'details' : """Data created by Brian Ferris and Dieter Fox. Consists of WiFi access point strengths taken during a circuit of the Paul Allen building at the University of Washington.""",
+                                      'license' : None,
+                                      'size' : 284390},
+                  'swiss_roll' : {'urls' : ['http://isomap.stanford.edu/'],
+                                  'files' : [['swiss_roll_data.mat']],
+                                  'details' : """Swiss roll data made available by Tenenbaum, de Silva and Langford to demonstrate isomap, available from http://isomap.stanford.edu/datasets.html.""",
+                                  'citation' : 'A Global Geometric Framework for Nonlinear Dimensionality Reduction, J. B. Tenenbaum, V. de Silva and J. C. Langford, Science 290 (5500): 2319-2323, 22 December 2000',
+                                  'license' : None,
+                                  'size' : 800256},
+                  'ripley_prnn_data' : {'urls' : ['http://www.stats.ox.ac.uk/pub/PRNN/'],
+                                        'files' : [['Cushings.dat', 'README', 'crabs.dat', 'fglass.dat', 'fglass.grp', 'pima.te', 'pima.tr', 'pima.tr2', 'synth.te', 'synth.tr', 'viruses.dat', 'virus3.dat']],
+                                        'details' : """Data sets from Brian Ripley's Pattern Recognition and Neural Networks""",
+                                        'citation': """Pattern Recognition and Neural Networks by B.D. Ripley (1996) Cambridge University Press ISBN 0 521 46986 7""",
+                                        'license' : None,
+                                        'size' : 93565},
+                  'isomap_face_data' : {'urls' : [neil_url + 'isomap_face_data/'],
+                                        'files' : [['face_data.mat']],
+                                        'details' : """Face data made available by Tenenbaum, de Silva and Langford to demonstrate isomap, available from http://isomap.stanford.edu/datasets.html.""",
+                                        'citation' : 'A Global Geometric Framework for Nonlinear Dimensionality Reduction, J. B. Tenenbaum, V. de Silva and J. C. Langford, Science 290 (5500): 2319-2323, 22 December 2000',
+                                        'license' : None,
+                                        'size' : 24229368},
+                  }
                   
 def prompt_user():
+    """Ask user for agreeing to data set licenses."""
     # raw_input returns the empty string for "enter"
     yes = set(['yes', 'y'])
     no = set(['no','n'])
@@ -25,45 +143,134 @@ def prompt_user():
         sys.stdout.write("Please respond with 'yes', 'y' or 'no', 'n'")
         return prompt_user()
 
-def download_data(dataset_name=None):
-    """Helper function which contains the resource locations for each data set in one place"""
-
-    # Note: there may be a better way of doing this. One of the pythonistas will need to take a look. Neil
-    data_resources = {'oil': {'urls' : [neil_url + 'oil_data/'],
-                              'files' : [['DataTrnLbls.txt', 'DataTrn.txt']],
-                              'citation' : 'Bishop, C. M. and G. D. James (1993). Analysis of multiphase flows using dual-energy gamma densitometry and neural networks. Nuclear Instruments and Methods in Physics Research A327, 580-593',
-                              'details' : """The three phase oil data used initially for demonstrating the Generative Topographic mapping.""",
-                              'agreement' : None},
-                      'brendan_faces' : {'url' : ['http://www.cs.nyu.edu/~roweis/data/'],
-                                         'files': [['frey_rawface.mat']],
-                                         'citation' : 'Frey, B. J., Colmenarez, A and Huang, T. S. Mixtures of Local Linear Subspaces for Face Recognition. Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition 1998, 32-37, June 1998. Computer Society Press, Los Alamitos, CA.',
-                                         'details' : """A video of Brendan Frey's face popularized as a benchmark for visualization by the Locally Linear Embedding.""",
-                                         'agreement': None}
-                      }
+def data_available(dataset_name=None):
+    """Check if the data set is available on the local machine already."""
+    for file_list in data_resources[dataset_name]['files']:
+        for file in file_list:
+            if not os.path.exists(os.path.join(data_path, dataset_name, file)):
+                return False
+    return True
             
+def download_url(url, store_directory, save_name = None, messages = True, suffix=''):
+    """Download a file from a url and save it to disk."""
+    i = url.rfind('/')
+    file = url[i+1:]
+    print file
+    dir_name = os.path.join(data_path, store_directory)
+    save_name = os.path.join(dir_name, file)
+    print "Downloading ", url, "->", os.path.join(store_directory, file)
+    if not os.path.exists(dir_name):
+        os.makedirs(dir_name)
+    urllib.urlretrieve(url+suffix, save_name, reporthook)
 
+def authorize_download(dataset_name=None):
+    """Check with the user that the are happy with terms and conditions for the data set."""
     print('Acquiring resource: ' + dataset_name)
     # TODO, check resource is in dictionary!
+    print('')
     dr = data_resources[dataset_name]
     print('Details of data: ')
     print(dr['details'])
+    print('')
     if dr['citation']:
         print('Please cite:')
         print(dr['citation'])
-    if dr['agreement']:
-        print('You must also agree to the following:')
-        print(dr['agreement'])
+        print('')
+    if dr['size']:
+        print('After downloading the data will take up ' + str(dr['size']) + ' bytes of space.')
+        print('')
+    print('Data will be stored in ' + os.path.join(data_path, dataset_name) + '.')
+    print('')
+    if overide_manual_authorize:
+        if dr['license']:
+            print('You have agreed to the following license:')
+            print(dr['license'])
+            print('')
+        return True
+    else:
+        if dr['license']:
+            print('You must also agree to the following license:')
+            print(dr['license'])
+            print('')
         print('Do you wish to proceed with the download? [yes/no]')
-    if prompt_user()==False:
+        return prompt_user()
+
+def download_data(dataset_name=None):
+    """Check with the user that the are happy with terms and conditions for the data set, then download it."""
+
+    dr = data_resources[dataset_name]
+    if not authorize_download(dataset_name):
         return False
 
+    if dr.has_key('suffices'):
+        for url, files, suffices in zip(dr['urls'], dr['files'], dr['suffices']):
+            for file, suffix in zip(files, suffices):
+                download_url(os.path.join(url,file), dataset_name, dataset_name, suffix=suffix)
+    else:
         for url, files in zip(dr['urls'], dr['files']):
             for file in files:
-            download_resource(url + file)
+                download_url(os.path.join(url,file), dataset_name, dataset_name)
     return True
                   
+def data_details_return(data, data_set):
+    """Update the data component of the data dictionary with details drawn from the data_resources."""
+    data.update(data_resources[data_set])
+    return data
 
     
+def cmu_urls_files(subj_motions, messages = True):
+    '''
+    Find which resources are missing on the local disk for the requested CMU motion capture motions. 
+    '''
+    
+    subjects_num = subj_motions[0]
+    motions_num = subj_motions[1]
+
+    resource = {'urls' : [], 'files' : []}
+    # Convert numbers to strings
+    subjects = []
+    motions = [list() for _ in range(len(subjects_num))]
+    for i in range(len(subjects_num)):
+        curSubj = str(int(subjects_num[i]))
+        if int(subjects_num[i]) < 10:
+            curSubj = '0' + curSubj
+        subjects.append(curSubj)
+        for j in range(len(motions_num[i])):
+            curMot = str(int(motions_num[i][j]))
+            if int(motions_num[i][j]) < 10:
+                curMot = '0' + curMot
+            motions[i].append(curMot)
+
+    all_skels = []
+    
+    assert len(subjects) == len(motions)
+    
+    all_motions = []
+            
+    for i in range(len(subjects)):
+        skel_dir = os.path.join(data_path, 'cmu_mocap')
+        cur_skel_file = os.path.join(skel_dir, subjects[i] + '.asf')
+        
+        url_required = False
+        file_download = []
+        if not os.path.exists(cur_skel_file):
+            # Current skel file doesn't exist.
+            if not os.path.isdir(skel_dir):
+                os.mkdir(skel_dir)
+            # Add skel file to list.
+            url_required = True
+            file_download.append(subjects[i] + '.asf')
+        for j in range(len(motions[i])):
+            file_name = subjects[i] + '_' + motions[i][j] + '.amc'
+            cur_motion_file = os.path.join(skel_dir, file_name)
+            if not os.path.exists(cur_motion_file):
+                url_required = True
+                file_download.append(subjects[i] + '_' + motions[i][j] + '.amc')
+        if url_required:
+            resource['urls'].append(cmu_url + subjects[i] + '/')
+            resource['files'].append(file_download)
+    return resource
+
 
 # Some general utilities.
 def sample_class(f):
@@ -72,25 +279,25 @@ def sample_class(f):
     c = np.where(c, 1, -1)
     return c
 
-def download_resource(resource, save_name = None, save_file = True, messages = True):
-    if messages:
-        print "Downloading resource: " , resource, " ... ",
-    response = url.urlopen(resource)
-    # TODO: Some error checking...
-    # ...
-    html = response.read()
-    response.close()
-    if save_file:
-        # TODO: Check if already exists...
-        # ...
-        with open(save_name, "w") as text_file:
-            text_file.write("%s"%html)
-            if messages:
-                print "Done!"
-    return html
+def boston_housing(data_set='boston_housing'):
+    if not data_available(data_set):
+        download_data(data_set)
+    all_data = np.genfromtxt(os.path.join(data_path, data_set, 'housing.data'))
+    X = all_data[:, 0:13]
+    Y = all_data[:, 13:14]
+    return data_details_return({'X' : X, 'Y': Y}, data_set)
 
-def della_gatta_TRP63_gene_expression(gene_number=None):
-    mat_data = scipy.io.loadmat(os.path.join(data_path, 'DellaGattadata.mat'))
+def brendan_faces(data_set='brendan_faces'):
+    if not data_available(data_set):
+        download_data(data_set)
+    mat_data = scipy.io.loadmat(os.path.join(data_path, data_set, 'frey_rawface.mat'))
+    Y = mat_data['ff'].T
+    return data_details_return({'Y': Y}, data_set)
+
+def della_gatta_TRP63_gene_expression(data_set='della_gatta', gene_number=None):
+    if not data_available(data_set):
+        download_data(data_set)
+    mat_data = scipy.io.loadmat(os.path.join(data_path, data_set, 'DellaGattadata.mat'))
     X = np.double(mat_data['timepoints'])
     if gene_number == None:
         Y = mat_data['exprs_tp53_RMA']
@@ -98,45 +305,62 @@ def della_gatta_TRP63_gene_expression(gene_number=None):
         Y = mat_data['exprs_tp53_RMA'][:, gene_number]
         if len(Y.shape) == 1:
             Y = Y[:, None]
-    return {'X': X, 'Y': Y, 'info': "The full gene expression data set from della Gatta et al (http://www.ncbi.nlm.nih.gov/pmc/articles/PMC2413161/) processed by RMA."}
+    return data_details_return({'X': X, 'Y': Y, 'gene_number' : gene_number}, data_set)
 
-def simulation_BGPLVM():
-    mat_data = scipy.io.loadmat(os.path.join(data_path, 'BGPLVMSimulation.mat'))
-    Y = np.array(mat_data['Y'], dtype=float)
-    S = np.array(mat_data['initS'], dtype=float)
-    mu = np.array(mat_data['initMu'], dtype=float)
-    return {'Y': Y, 'S': S,
-            'mu' : mu,
-            'info': "Simulated test dataset generated in MATLAB to compare BGPLVM between python and MATLAB"}
 
 
 # The data sets
-def oil():
-    #if download_data('oil'):
-    oil_train_file = os.path.join(data_path, 'oil', 'DataTrn.txt')
-    oil_trainlbls_file = os.path.join(data_path, 'oil', 'DataTrnLbls.txt')
+def oil(data_set='three_phase_oil_flow'):
+    """The three phase oil data from Bishop and James (1993)."""
+    if not data_available(data_set):
+        download_data(data_set)
+    oil_train_file = os.path.join(data_path, data_set, 'DataTrn.txt')
+    oil_trainlbls_file = os.path.join(data_path, data_set, 'DataTrnLbls.txt')
+    oil_test_file = os.path.join(data_path, data_set, 'DataTst.txt')
+    oil_testlbls_file = os.path.join(data_path, data_set, 'DataTstLbls.txt')
+    oil_valid_file = os.path.join(data_path, data_set, 'DataVdn.txt')
+    oil_validlbls_file = os.path.join(data_path, data_set, 'DataVdnLbls.txt')
     fid = open(oil_train_file)
     X = np.fromfile(fid, sep='\t').reshape((-1, 12))
     fid.close()
+    fid = open(oil_test_file)
+    Xtest = np.fromfile(fid, sep='\t').reshape((-1, 12))
+    fid.close()
+    fid = open(oil_valid_file)
+    Xvalid = np.fromfile(fid, sep='\t').reshape((-1, 12))
+    fid.close()
     fid = open(oil_trainlbls_file)
     Y = np.fromfile(fid, sep='\t').reshape((-1, 3)) * 2. - 1.
     fid.close()
-    return {'X': X, 'Y': Y, 'info': "The oil data from Bishop and James (1993)."}
+    fid = open(oil_testlbls_file)
+    Ytest = np.fromfile(fid, sep='\t').reshape((-1, 3)) * 2. - 1.
+    fid.close()
+    fid = open(oil_validlbls_file)
+    Yvalid = np.fromfile(fid, sep='\t').reshape((-1, 3)) * 2. - 1.
+    fid.close()
+    return data_details_return({'X': X, 'Y': Y, 'Xtest': Xtest, 'Ytest': Ytest, 'Xtest' : Xtest, 'Xvalid': Xvalid, 'Yvalid': Yvalid}, data_set)
     #else:
     # throw an error
     
-def oil_100(seed=default_seed):
+def oil_100(seed=default_seed, data_set = 'three_phase_oil_flow'):
     np.random.seed(seed=seed)
     data = oil()
     indices = np.random.permutation(1000)
     indices = indices[0:100]
     X = data['X'][indices, :]
     Y = data['Y'][indices, :]
-    return {'X': X, 'Y': Y, 'info': "Subsample of the oil data extracting 100 values randomly without replacement."}
+    return data_details_return({'X': X, 'Y': Y, 'info': "Subsample of the full oil data extracting 100 values randomly without replacement, here seed was " + str(seed)}, data_set)
 
-def pumadyn(seed=default_seed):
+def pumadyn(seed=default_seed, data_set='pumadyn-32nm'):
+    if not data_available(data_set):
+        download_data(data_set)
+        path = os.path.join(data_path, data_set)
+        tar = tarfile.open(os.path.join(path, 'pumadyn-32nm.tar.gz'))
+        print('Extracting file.')
+        tar.extractall(path=path)
+        tar.close()
     # Data is variance 1, no need to normalize.
-    data = np.loadtxt(os.path.join(data_path, 'pumadyn-32nm/Dataset.data.gz'))
+    data = np.loadtxt(os.path.join(data_path, data_set, 'pumadyn-32nm', 'Dataset.data.gz'))
     indices = np.random.permutation(data.shape[0])
     indicesTrain = indices[0:7168]
     indicesTest = indices[7168:-1]
@@ -146,20 +370,54 @@ def pumadyn(seed=default_seed):
     Y = data[indicesTrain, -1][:, None]
     Xtest = data[indicesTest, 0:-2]
     Ytest = data[indicesTest, -1][:, None]
-    return {'X': X, 'Y': Y, 'Xtest': Xtest, 'Ytest': Ytest, 'info': "The puma robot arm data with 32 inputs. This data is the non linear case with medium noise (pumadyn-32nm). For training 7,168 examples are sampled without replacement."}
+    return data_details_return({'X': X, 'Y': Y, 'Xtest': Xtest, 'Ytest': Ytest, 'seed': seed}, data_set)
 
+def robot_wireless(data_set='robot_wireless'):
+    # WiFi access point strengths on a tour around UW Paul Allen building.
+    if not data_available(data_set):
+        download_data(data_set)
+    file_name = os.path.join(data_path, data_set, 'uw-floor.txt')
+    all_time = np.genfromtxt(file_name, usecols=(0))
+    macaddress = np.genfromtxt(file_name, usecols=(1), dtype='string')
+    x = np.genfromtxt(file_name, usecols=(2))
+    y = np.genfromtxt(file_name, usecols=(3))
+    strength = np.genfromtxt(file_name, usecols=(4))
+    addresses = np.unique(macaddress)
+    times = np.unique(all_time)
+    addresses.sort()
+    times.sort()
+    allY = np.zeros((len(times), len(addresses)))
+    allX = np.zeros((len(times), 2))
+    allY[:]=-92.
+    strengths={}
+    for address, j in zip(addresses, range(len(addresses))):
+        ind = np.nonzero(address==macaddress)
+        temp_strengths=strength[ind]
+        temp_x=x[ind]
+        temp_y=y[ind]
+        temp_times = all_time[ind]
+        for time in temp_times:
+            vals = time==temp_times
+            if any(vals):
+                ind2 = np.nonzero(vals)
+                i = np.nonzero(time==times)
+                allY[i, j] = temp_strengths[ind2]
+                allX[i, 0] = temp_x[ind2]
+                allX[i, 1] = temp_y[ind2]
+    allY = (allY + 85.)/15.
 
-def brendan_faces():
-    mat_data = scipy.io.loadmat(os.path.join(data_path, 'frey_rawface.mat'))
-    Y = mat_data['ff'].T
-    return {'Y': Y, 'info': "Face data made available by Brendan Frey"}
+    X = allX[0:215, :]
+    Y = allY[0:215, :]
 
+    Xtest = allX[215:, :]
+    Ytest = allY[215:, :]
+    return data_details_return({'X': X, 'Y': Y, 'Xtest': Xtest, 'Ytest': Ytest, 'addresses' : addresses, 'times' : times}, data_set)
 
-
-
-def silhouette():
+def silhouette(data_set='ankur_pose_data'):
     # Ankur Agarwal and Bill Trigg's silhoutte data.
-    mat_data = scipy.io.loadmat(os.path.join(data_path, 'mocap', 'ankur', 'ankurDataPoseSilhouette.mat'))
+    if not data_available(data_set):
+        download_data(data_set)
+    mat_data = scipy.io.loadmat(os.path.join(data_path, data_set, 'ankurDataPoseSilhouette.mat'))
     inMean = np.mean(mat_data['Y'])
     inScales = np.sqrt(np.var(mat_data['Y']))
     X = mat_data['Y'] - inMean
@@ -168,22 +426,35 @@ def silhouette():
     Xtest = Xtest / inScales
     Y = mat_data['Z']
     Ytest = mat_data['Z_test']
-    return {'X': X, 'Y': Y, 'Xtest': Xtest, 'Ytest': Ytest, 'info': "Artificial silhouette simulation data developed from Agarwal and Triggs (2004)."}
+    return data_details_return({'X': X, 'Y': Y, 'Xtest': Xtest, 'Ytest': Ytest}, data_set)
 
-def stick():
-    #if download_data('stick'):
-    Y, connect = GPy.util.mocap.load_text_data('run1', data_path)
+def ripley_synth(data_set='ripley_prnn_data'):
+    if not data_available(data_set):
+        download_data(data_set)
+    train = np.genfromtxt(os.path.join(data_path, data_set, 'synth.tr'), skip_header=1)
+    X = train[:, 0:2]
+    y = train[:, 2:3]
+    test = np.genfromtxt(os.path.join(data_path, data_set, 'synth.te'), skip_header=1)
+    Xtest = test[:, 0:2]
+    ytest = test[:, 2:3]
+    return data_details_return({'X': X, 'y': y, 'Xtest': Xtest, 'ytest': ytest, 'info': 'Synthetic data generated by Ripley for a two class classification problem.'}, data_set)
+
+def osu_run1(data_set='osu_run1'):
+    if not data_available(data_set):
+        download_data(data_set)
+    zip = zipfile.ZipFile(os.path.join(data_path, data_set, 'sprintTXT.ZIP'), 'r')
+    path = os.path.join(data_path, data_set)
+    for name in zip.namelist():
+        zip.extract(name, path)
+    Y, connect = GPy.util.mocap.load_text_data('Aug210107', path)
     Y = Y[0:-1:4, :]
-    lbls = 'connect'
-    return {'Y': Y, 'connect' : connect, 'info': "Stick man data from Ohio."}
-    # else:
-    # throw an error.
+    return data_details_return({'Y': Y, 'connect' : connect}, data_set)
 
-def swiss_roll_generated(N=1000, sigma=0.0):
+def swiss_roll_generated(num_samples=1000, sigma=0.0):
     with open(os.path.join(data_path, 'swiss_roll.pickle')) as f:
         data = pickle.load(f)
     Na = data['Y'].shape[0]
-    perm = np.random.permutation(np.r_[:Na])[:N]
+    perm = np.random.permutation(np.r_[:Na])[:num_samples]
     Y = data['Y'][perm, :]
     t = data['t'][perm]
     c = data['colors'][perm, :]
@@ -194,27 +465,49 @@ def swiss_roll_generated(N=1000, sigma=0.0):
     return {'Y':Y, 't':t, 'colors':c}
 
 def swiss_roll_1000():
-    mat_data = scipy.io.loadmat(os.path.join(data_path, 'swiss_roll_data'))
-    Y = mat_data['X_data'][:, 0:1000].transpose()
-    return {'Y': Y, 'info': "Subsample of the swiss roll data extracting only the first 1000 values."}
+    return swiss_roll(num_samples=1000)
 
-def swiss_roll(N=3000):
-    mat_data = scipy.io.loadmat(os.path.join(data_path, 'swiss_roll_data.mat'))
-    Y = mat_data['X_data'][:, 0:N].transpose()
-    return {'Y': Y, 'X': mat_data['X_data'], 'info': "The first 3,000 points from the swiss roll data of Tennenbaum, de Silva and Langford (2001)."}
+def swiss_roll(num_samples=3000, data_set='swiss_roll'):
+    if not data_available(data_set):
+        download_data(data_set)
+    mat_data = scipy.io.loadmat(os.path.join(data_path, data_set, 'swiss_roll_data.mat'))
+    Y = mat_data['X_data'][:, 0:num_samples].transpose()
+    return data_details_return({'Y': Y, 'X': mat_data['X_data'], 'info': "The first " + str(num_samples) + " points from the swiss roll data of Tennenbaum, de Silva and Langford (2001)."}, data_set)
 
-def toy_rbf_1d(seed=default_seed):
+def isomap_faces(num_samples=698, data_set='isomap_face_data'):
+    if not data_available(data_set):
+        download_data(data_set)
+    mat_data = scipy.io.loadmat(os.path.join(data_path, data_set, 'face_data.mat'))
+    Y = mat_data['images'][:, 0:num_samples].transpose()
+    return data_details_return({'Y': Y, 'poses' : mat_data['poses'], 'lights': mat_data['lights'], 'info': "The first " + str(num_samples) + " points from the face data of Tennenbaum, de Silva and Langford (2001)."}, data_set)
+
+def simulation_BGPLVM():
+    mat_data = scipy.io.loadmat(os.path.join(data_path, 'BGPLVMSimulation.mat'))
+    Y = np.array(mat_data['Y'], dtype=float)
+    S = np.array(mat_data['initS'], dtype=float)
+    mu = np.array(mat_data['initMu'], dtype=float)
+    return data_details_return({'S': S, 'Y': Y, 'mu': mu}, data_set)
+    return {'Y': Y, 'S': S,
+            'mu' : mu,
+            'info': "Simulated test dataset generated in MATLAB to compare BGPLVM between python and MATLAB"}
+
+def toy_rbf_1d(seed=default_seed, num_samples=500):
+    """Samples values of a function from an RBF covariance with very small noise for inputs uniformly distributed between -1 and 1.
+    :param seed: seed to use for random sampling.
+    :type seed: int
+    :param num_samples: number of samples to sample in the function (default 500).
+    :type num_samples: int
+    """
     np.random.seed(seed=seed)
-    numIn = 1
-    N = 500
-    X = np.random.uniform(low= -1.0, high=1.0, size=(N, numIn))
+    num_in = 1
+    X = np.random.uniform(low= -1.0, high=1.0, size=(num_samples, num_in))
     X.sort(axis=0)
-    rbf = GPy.kern.rbf(numIn, variance=1., lengthscale=np.array((0.25,)))
-    white = GPy.kern.white(numIn, variance=1e-2)
+    rbf = GPy.kern.rbf(num_in, variance=1., lengthscale=np.array((0.25,)))
+    white = GPy.kern.white(num_in, variance=1e-2)
     kernel = rbf + white
     K = kernel.K(X)
-    y = np.reshape(np.random.multivariate_normal(np.zeros(N), K), (N, 1))
-    return {'X':X, 'Y':y, 'info': "Samples 500 values of a function from an RBF covariance with very small noise for inputs uniformly distributed between -1 and 1."}
+    y = np.reshape(np.random.multivariate_normal(np.zeros(num_samples), K), (num_samples, 1))
+    return {'X':X, 'Y':y, 'info': "Sampled " + str(num_samples) + " values of a function from an RBF covariance with very small noise for inputs uniformly distributed between -1 and 1."}
 
 def toy_rbf_1d_50(seed=default_seed):
     np.random.seed(seed=seed)
@@ -224,7 +517,7 @@ def toy_rbf_1d_50(seed=default_seed):
     indices.sort(axis=0)
     X = data['X'][indices, :]
     Y = data['Y'][indices, :]
-    return {'X': X, 'Y': Y, 'info': "Subsamples the toy_rbf_sample with 50 values randomly taken from the original sample."}
+    return {'X': X, 'Y': Y, 'info': "Subsamples the toy_rbf_sample with 50 values randomly taken from the original sample.", 'seed' : seed}
 
 
 def toy_linear_1d_classification(seed=default_seed):
@@ -232,13 +525,46 @@ def toy_linear_1d_classification(seed=default_seed):
     x1 = np.random.normal(-3, 5, 20)
     x2 = np.random.normal(3, 5, 20)
     X = (np.r_[x1, x2])[:, None]
-    return {'X': X, 'Y':  sample_class(2.*X), 'F': 2.*X}
+    return {'X': X, 'Y':  sample_class(2.*X), 'F': 2.*X, 'seed' : seed}
+
+def olympic_100m_men(data_set='rogers_girolami_data'):
+    if not data_available(data_set):
+        download_data(data_set)
+        path = os.path.join(data_path, data_set)
+        tar_file = os.path.join(path, 'firstcoursemldata.tar.gz')
+        tar = tarfile.open(tar_file)
+        print('Extracting file.')
+        tar.extractall(path=path)
+        tar.close()
+    olympic_data = scipy.io.loadmat(os.path.join(data_path, data_set, 'data', 'olympics.mat'))['male100']
 
-def rogers_girolami_olympics():
-    olympic_data = scipy.io.loadmat(os.path.join(data_path, 'olympics.mat'))['male100']
     X = olympic_data[:, 0][:, None]
     Y = olympic_data[:, 1][:, None]
-    return {'X': X, 'Y': Y, 'info': "Olympic sprint times for 100 m men from 1896 until 2008. Example is from Rogers and Girolami's First Course in Machine Learning."}
+    return data_details_return({'X': X, 'Y': Y, 'info': "Olympic sprint times for 100 m men from 1896 until 2008. Example is from Rogers and Girolami's First Course in Machine Learning."}, data_set)
+
+def olympic_100m_women(data_set='rogers_girolami_data'):
+    if not data_available(data_set):
+        download_data(data_set)
+        path = os.path.join(data_path, data_set)
+        tar_file = os.path.join(path, 'firstcoursemldata.tar.gz')
+        tar = tarfile.open(tar_file)
+        print('Extracting file.')
+        tar.extractall(path=path)
+        tar.close()
+    olympic_data = scipy.io.loadmat(os.path.join(data_path, data_set, 'data', 'olympics.mat'))['female100']
+
+    X = olympic_data[:, 0][:, None]
+    Y = olympic_data[:, 1][:, None]
+    return data_details_return({'X': X, 'Y': Y, 'info': "Olympic sprint times for 100 m women from 1896 until 2008. Example is from Rogers and Girolami's First Course in Machine Learning."}, data_set)
+
+def olympic_marathon_men(data_set='olympic_marathon_men'):
+    if not data_available(data_set):
+        download_data(data_set)
+    olympics = np.genfromtxt(os.path.join(data_path, data_set, 'olympicMarathonTimes.csv'), delimiter=',')
+    X = olympics[:, 0:1]
+    Y = olympics[:, 1:2]
+    return data_details_return({'X': X, 'Y': Y}, data_set)
+
 # def movielens_small(partNo=1,seed=default_seed):
 #     np.random.seed(seed=seed)
 
@@ -272,8 +598,6 @@ def rogers_girolami_olympics():
 #     return {'Y':Y, 'lbls':lbls, 'Ytest':Ytest, 'lblstest':lblstest}
 
 
-
-
 def crescent_data(num_data=200, seed=default_seed):
     """Data set formed from a mixture of four Gaussians. In each class two of the Gaussians are elongated at right angles to each other and offset to form an approximation to the crescent data that is popular in semi-supervised learning as a toy problem.
     :param num_data_part: number of data to be sampled (default is 200).
@@ -302,24 +626,31 @@ def crescent_data(num_data=200, seed=default_seed):
     for i in range(0, 4):
         num_data_part.append(round(((i + 1) * num_data) / 4.))
         num_data_part[i] -= num_data_total
-        # print num_data_part[i]
         part = np.random.normal(size=(num_data_part[i], 2))
         part = np.dot(np.dot(part, scales[i]), R) + means[i]
         Xparts.append(part)
         num_data_total += num_data_part[i]
     X = np.vstack((Xparts[0], Xparts[1], Xparts[2], Xparts[3]))
 
-
     Y = np.vstack((np.ones((num_data_part[0] + num_data_part[1], 1)), -np.ones((num_data_part[2] + num_data_part[3], 1))))
     return {'X':X, 'Y':Y, 'info': "Two separate classes of data formed approximately in the shape of two crescents."}
 
-def creep_data():
-    all_data = np.loadtxt(os.path.join(data_path, 'creep', 'taka'))
+def creep_data(data_set='creep_rupture'):
+    """Brun and Yoshida's metal creep rupture data."""
+    if not data_available(data_set):
+        download_data(data_set)
+        path = os.path.join(data_path, data_set)
+        tar_file = os.path.join(path, 'creeprupt.tar')
+        tar = tarfile.open(tar_file)
+        print('Extracting file.')
+        tar.extractall(path=path)
+        tar.close()
+    all_data = np.loadtxt(os.path.join(data_path, data_set, 'taka'))
     y = all_data[:, 1:2].copy()
     features = [0]
     features.extend(range(2, 31))
     X = all_data[:, features].copy()
-    return {'X': X, 'y' : y}
+    return data_details_return({'X': X, 'y': y}, data_set)
 
 def cmu_mocap_49_balance():
     """Load CMU subject 49's one legged balancing motion that was used by Alvarez, Luengo and Lawrence at AISTATS 2009."""
@@ -341,14 +672,19 @@ def cmu_mocap_35_walk_jog():
     data['info'] = "Walk and jog data from CMU data base subject 35. As used in Tayor, Roweis and Hinton at NIPS 2007, but without their pre-processing (i.e. as used by Lawrence at AISTATS 2007). It consists of " + data['info']
     return data
 
-def cmu_mocap(subject, train_motions, test_motions=[], sample_every=4):
+def cmu_mocap(subject, train_motions, test_motions=[], sample_every=4, data_set='cmu_mocap'):
     """Load a given subject's training and test motions from the CMU motion capture data."""
-
     # Load in subject skeleton.
-    subject_dir = os.path.join(data_path, 'mocap', 'cmu', subject)
+    subject_dir = os.path.join(data_path, data_set)
 
     # Make sure the data is downloaded.
-    mocap.fetch_cmu(([subject], [train_motions]), skel_store_dir=subject_dir,motion_store_dir=subject_dir)
+    all_motions = train_motions + test_motions
+    resource = cmu_urls_files(([subject], [all_motions]))
+    data_resources[data_set] = data_resources['cmu_mocap_full']
+    data_resources[data_set]['files'] = resource['files']
+    data_resources[data_set]['urls'] = resource['urls']
+    if resource['urls']:
+        download_data(data_set)
 
     skel = GPy.util.mocap.acclaim_skeleton(os.path.join(subject_dir, subject + '.asf'))
 
@@ -413,4 +749,4 @@ def cmu_mocap(subject, train_motions, test_motions=[], sample_every=4):
         info += '.'
     if sample_every != 1:
         info += ' Data is sub-sampled to every ' + str(sample_every) + ' frames.'
-    return {'Y': Y, 'lbls' : lbls, 'Ytest': Ytest, 'lblstest' : lblstest, 'info': info, 'skel': skel}
+    return data_details_return({'Y': Y, 'lbls' : lbls, 'Ytest': Ytest, 'lblstest' : lblstest, 'info': info, 'skel': skel}, data_set)

GPy/util/datasets/crabs.dat

@@ -1,201 +0,0 @@
-sp sex index FL RW CL CW BD
-B M  1  8.1  6.7 16.1 19.0  7.0
-B M  2  8.8  7.7 18.1 20.8  7.4
-B M  3  9.2  7.8 19.0 22.4  7.7
-B M  4  9.6  7.9 20.1 23.1  8.2
-B M  5  9.8  8.0 20.3 23.0  8.2
-B M  6 10.8  9.0 23.0 26.5  9.8
-B M  7 11.1  9.9 23.8 27.1  9.8
-B M  8 11.6  9.1 24.5 28.4 10.4
-B M  9 11.8  9.6 24.2 27.8  9.7
-B M 10 11.8 10.5 25.2 29.3 10.3
-B M 11 12.2 10.8 27.3 31.6 10.9
-B M 12 12.3 11.0 26.8 31.5 11.4
-B M 13 12.6 10.0 27.7 31.7 11.4
-B M 14 12.8 10.2 27.2 31.8 10.9
-B M 15 12.8 10.9 27.4 31.5 11.0
-B M 16 12.9 11.0 26.8 30.9 11.4
-B M 17 13.1 10.6 28.2 32.3 11.0
-B M 18 13.1 10.9 28.3 32.4 11.2
-B M 19 13.3 11.1 27.8 32.3 11.3
-B M 20 13.9 11.1 29.2 33.3 12.1
-B M 21 14.3 11.6 31.3 35.5 12.7
-B M 22 14.6 11.3 31.9 36.4 13.7
-B M 23 15.0 10.9 31.4 36.4 13.2
-B M 24 15.0 11.5 32.4 37.0 13.4
-B M 25 15.0 11.9 32.5 37.2 13.6
-B M 26 15.2 12.1 32.3 36.7 13.6
-B M 27 15.4 11.8 33.0 37.5 13.6
-B M 28 15.7 12.6 35.8 40.3 14.5
-B M 29 15.9 12.7 34.0 38.9 14.2
-B M 30 16.1 11.6 33.8 39.0 14.4
-B M 31 16.1 12.8 34.9 40.7 15.7
-B M 32 16.2 13.3 36.0 41.7 15.4
-B M 33 16.3 12.7 35.6 40.9 14.9
-B M 34 16.4 13.0 35.7 41.8 15.2
-B M 35 16.6 13.5 38.1 43.4 14.9
-B M 36 16.8 12.8 36.2 41.8 14.9
-B M 37 16.9 13.2 37.3 42.7 15.6
-B M 38 17.1 12.6 36.4 42.0 15.1
-B M 39 17.1 12.7 36.7 41.9 15.6
-B M 40 17.2 13.5 37.6 43.9 16.1
-B M 41 17.7 13.6 38.7 44.5 16.0
-B M 42 17.9 14.1 39.7 44.6 16.8
-B M 43 18.0 13.7 39.2 44.4 16.2
-B M 44 18.8 15.8 42.1 49.0 17.8
-B M 45 19.3 13.5 41.6 47.4 17.8
-B M 46 19.3 13.8 40.9 46.5 16.8
-B M 47 19.7 15.3 41.9 48.5 17.8
-B M 48 19.8 14.2 43.2 49.7 18.6
-B M 49 19.8 14.3 42.4 48.9 18.3
-B M 50 21.3 15.7 47.1 54.6 20.0
-B F  1  7.2  6.5 14.7 17.1  6.1
-B F  2  9.0  8.5 19.3 22.7  7.7
-B F  3  9.1  8.1 18.5 21.6  7.7
-B F  4  9.1  8.2 19.2 22.2  7.7
-B F  5  9.5  8.2 19.6 22.4  7.8
-B F  6  9.8  8.9 20.4 23.9  8.8
-B F  7 10.1  9.3 20.9 24.4  8.4
-B F  8 10.3  9.5 21.3 24.7  8.9
-B F  9 10.4  9.7 21.7 25.4  8.3
-B F 10 10.8  9.5 22.5 26.3  9.1
-B F 11 11.0  9.8 22.5 25.7  8.2
-B F 12 11.2 10.0 22.8 26.9  9.4
-B F 13 11.5 11.0 24.7 29.2 10.1
-B F 14 11.6 11.0 24.6 28.5 10.4
-B F 15 11.6 11.4 23.7 27.7 10.0
-B F 16 11.7 10.6 24.9 28.5 10.4
-B F 17 11.9 11.4 26.0 30.1 10.9
-B F 18 12.0 10.7 24.6 28.9 10.5
-B F 19 12.0 11.1 25.4 29.2 11.0
-B F 20 12.6 12.2 26.1 31.6 11.2
-B F 21 12.8 11.7 27.1 31.2 11.9
-B F 22 12.8 12.2 26.7 31.1 11.1
-B F 23 12.8 12.2 27.9 31.9 11.5
-B F 24 13.0 11.4 27.3 31.8 11.3
-B F 25 13.1 11.5 27.6 32.6 11.1
-B F 26 13.2 12.2 27.9 32.1 11.5
-B F 27 13.4 11.8 28.4 32.7 11.7
-B F 28 13.7 12.5 28.6 33.8 11.9
-B F 29 13.9 13.0 30.0 34.9 13.1
-B F 30 14.7 12.5 30.1 34.7 12.5
-B F 31 14.9 13.2 30.1 35.6 12.0
-B F 32 15.0 13.8 31.7 36.9 14.0
-B F 33 15.0 14.2 32.8 37.4 14.0
-B F 34 15.1 13.3 31.8 36.3 13.5
-B F 35 15.1 13.5 31.9 37.0 13.8
-B F 36 15.1 13.8 31.7 36.6 13.0
-B F 37 15.2 14.3 33.9 38.5 14.7
-B F 38 15.3 14.2 32.6 38.3 13.8
-B F 39 15.4 13.3 32.4 37.6 13.8
-B F 40 15.5 13.8 33.4 38.7 14.7
-B F 41 15.6 13.9 32.8 37.9 13.4
-B F 42 15.6 14.7 33.9 39.5 14.3
-B F 43 15.7 13.9 33.6 38.5 14.1
-B F 44 15.8 15.0 34.5 40.3 15.3
-B F 45 16.2 15.2 34.5 40.1 13.9
-B F 46 16.4 14.0 34.2 39.8 15.2
-B F 47 16.7 16.1 36.6 41.9 15.4
-B F 48 17.4 16.9 38.2 44.1 16.6
-B F 49 17.5 16.7 38.6 44.5 17.0
-B F 50 19.2 16.5 40.9 47.9 18.1
-O M  1  9.1  6.9 16.7 18.6  7.4
-O M  2 10.2  8.2 20.2 22.2  9.0
-O M  3 10.7  8.6 20.7 22.7  9.2
-O M  4 11.4  9.0 22.7 24.8 10.1
-O M  5 12.5  9.4 23.2 26.0 10.8
-O M  6 12.5  9.4 24.2 27.0 11.2
-O M  7 12.7 10.4 26.0 28.8 12.1
-O M  8 13.2 11.0 27.1 30.4 12.2
-O M  9 13.4 10.1 26.6 29.6 12.0
-O M 10 13.7 11.0 27.5 30.5 12.2
-O M 11 14.0 11.5 29.2 32.2 13.1
-O M 12 14.1 10.4 28.9 31.8 13.5
-O M 13 14.1 10.5 29.1 31.6 13.1
-O M 14 14.1 10.7 28.7 31.9 13.3 
-O M 15 14.2 10.6 28.7 31.7 12.9
-O M 16 14.2 10.7 27.8 30.9 12.7
-O M 17 14.2 11.3 29.2 32.2 13.5
-O M 18 14.6 11.3 29.9 33.5 12.8
-O M 19 14.7 11.1 29.0 32.1 13.1
-O M 20 15.1 11.4 30.2 33.3 14.0
-O M 21 15.1 11.5 30.9 34.0 13.9
-O M 22 15.4 11.1 30.2 33.6 13.5
-O M 23 15.7 12.2 31.7 34.2 14.2
-O M 24 16.2 11.8 32.3 35.3 14.7
-O M 25 16.3 11.6 31.6 34.2 14.5
-O M 26 17.1 12.6 35.0 38.9 15.7
-O M 27 17.4 12.8 36.1 39.5 16.2
-O M 28 17.5 12.0 34.4 37.3 15.3
-O M 29 17.5 12.7 34.6 38.4 16.1
-O M 30 17.8 12.5 36.0 39.8 16.7
-O M 31 17.9 12.9 36.9 40.9 16.5
-O M 32 18.0 13.4 36.7 41.3 17.1
-O M 33 18.2 13.7 38.8 42.7 17.2
-O M 34 18.4 13.4 37.9 42.2 17.7
-O M 35 18.6 13.4 37.8 41.9 17.3
-O M 36 18.6 13.5 36.9 40.2 17.0
-O M 37 18.8 13.4 37.2 41.1 17.5
-O M 38 18.8 13.8 39.2 43.3 17.9
-O M 39 19.4 14.1 39.1 43.2 17.8
-O M 40 19.4 14.4 39.8 44.3 17.9
-O M 41 20.1 13.7 40.6 44.5 18.0
-O M 42 20.6 14.4 42.8 46.5 19.6
-O M 43 21.0 15.0 42.9 47.2 19.4
-O M 44 21.5 15.5 45.5 49.7 20.9
-O M 45 21.6 15.4 45.7 49.7 20.6
-O M 46 21.6 14.8 43.4 48.2 20.1
-O M 47 21.9 15.7 45.4 51.0 21.1
-O M 48 22.1 15.8 44.6 49.6 20.5
-O M 49 23.0 16.8 47.2 52.1 21.5
-O M 50 23.1 15.7 47.6 52.8 21.6
-O F  1 10.7  9.7 21.4 24.0  9.8
-O F  2 11.4  9.2 21.7 24.1  9.7
-O F  3 12.5 10.0 24.1 27.0 10.9
-O F  4 12.6 11.5 25.0 28.1 11.5
-O F  5 12.9 11.2 25.8 29.1 11.9
-O F  6 14.0 11.9 27.0 31.4 12.6
-O F  7 14.0 12.8 28.8 32.4 12.7
-O F  8 14.3 12.2 28.1 31.8 12.5
-O F  9 14.7 13.2 29.6 33.4 12.9
-O F 10 14.9 13.0 30.0 33.7 13.3
-O F 11 15.0 12.3 30.1 33.3 14.0
-O F 12 15.6 13.5 31.2 35.1 14.1
-O F 13 15.6 14.0 31.6 35.3 13.8
-O F 14 15.6 14.1 31.0 34.5 13.8
-O F 15 15.7 13.6 31.0 34.8 13.8
-O F 16 16.1 13.6 31.6 36.0 14.0
-O F 17 16.1 13.7 31.4 36.1 13.9
-O F 18 16.2 14.0 31.6 35.6 13.7
-O F 19 16.7 14.3 32.3 37.0 14.7
-O F 20 17.1 14.5 33.1 37.2 14.6
-O F 21 17.5 14.3 34.5 39.6 15.6
-O F 22 17.5 14.4 34.5 39.0 16.0
-O F 23 17.5 14.7 33.3 37.6 14.6
-O F 24 17.6 14.0 34.0 38.6 15.5
-O F 25 18.0 14.9 34.7 39.5 15.7
-O F 26 18.0 16.3 37.9 43.0 17.2
-O F 27 18.3 15.7 35.1 40.5 16.1
-O F 28 18.4 15.5 35.6 40.0 15.9
-O F 29 18.4 15.7 36.5 41.6 16.4
-O F 30 18.5 14.6 37.0 42.0 16.6
-O F 31 18.6 14.5 34.7 39.4 15.0
-O F 32 18.8 15.2 35.8 40.5 16.6
-O F 33 18.9 16.7 36.3 41.7 15.3
-O F 34 19.1 16.0 37.8 42.3 16.8
-O F 35 19.1 16.3 37.9 42.6 17.2
-O F 36 19.7 16.7 39.9 43.6 18.2
-O F 37 19.9 16.6 39.4 43.9 17.9
-O F 38 19.9 17.9 40.1 46.4 17.9
-O F 39 20.0 16.7 40.4 45.1 17.7
-O F 40 20.1 17.2 39.8 44.1 18.6
-O F 41 20.3 16.0 39.4 44.1 18.0
-O F 42 20.5 17.5 40.0 45.5 19.2
-O F 43 20.6 17.5 41.5 46.2 19.2
-O F 44 20.9 16.5 39.9 44.7 17.5
-O F 45 21.3 18.4 43.8 48.4 20.0
-O F 46 21.4 18.0 41.2 46.2 18.7
-O F 47 21.7 17.1 41.7 47.2 19.6
-O F 48 21.9 17.2 42.6 47.4 19.5
-O F 49 22.5 17.2 43.0 48.7 19.8
-O F 50 23.1 20.2 46.2 52.5 21.1

GPy/util/datasets/pumadyn-32nm/Dataset.spec

@@ -1,44 +0,0 @@
-#
-# Puma forward dynamics -- 32nm = 32 inputs, high nonlinearity, med noise
-#
-#
-Origin: simulated
-
-Usage: assessment
-
-Order: uninformative
-
-Attributes:
-  1  theta1	      u  [-3.1416,3.1416]	# ang position of joint 1 in radians
-  2  theta2	      u  [-3.1416,3.1416]	# ang position of joint 2 in radians
-  3  theta3	      u  [-3.1416,3.1416]	# ang position of joint 3 in radians
-  4  theta4	      u  [-3.1416,3.1416]	# ang position of joint 4 in radians
-  5  theta5	      u  [-3.1416,3.1416]	# ang position of joint 5 in radians
-  6  theta6	      u  [-3.1416,3.1416]	# ang position of joint 6 in radians
-  7  thetad1      u  (-Inf,Inf)	# ang vel of joint 1 in rad/sec
-  8  thetad2      u  (-Inf,Inf)	# ang vel of joint 2 in rad/sec
-  9  thetad3      u  (-Inf,Inf)	# ang vel of joint 3 in rad/sec
-  10 thetad4      u  (-Inf,Inf)	# ang vel of joint 4 in rad/sec
-  11 thetad5      u  (-Inf,Inf)	# ang vel of joint 5 in rad/sec
-  12 thetad6      u  (-Inf,Inf)	# ang vel of joint 6 in rad/sec
-  13 tau1     u  (-Inf,Inf)	# torque on jt 1 
-  14 tau2     u  (-Inf,Inf)	# torque on jt 2
-  15 tau3     u  (-Inf,Inf)	# torque on jt 3
-  16 tau4     u  (-Inf,Inf)	# torque on jt 4
-  17 tau5     u  (-Inf,Inf)	# torque on jt 5
-  18 dm1   u	 [0,Inf)	# proportion change in mass of link 1
-  19 dm2   u	 [0,Inf)	# prop change in mass of link 2
-  20 dm3   u	 [0,Inf)	# prop change in mass of link 3
-  21 dm4   u	 [0,Inf)	# prop change in mass of link 4
-  22 dm5   u	 [0,Inf)	# prop change in mass of link 5
-  23 da1    u	 [0,Inf)	# prop change in length of link 1
-  24 da2    u	 [0,Inf)	# prop change in length of link 2
-  25 da3    u	 [0,Inf)	# prop change in length of link 3
-  26 da4    u	 [0,Inf)	# prop change in length of link 4
-  27 da5    u	 [0,Inf)	# prop change in length of link 5
-  28 db1      u	 [0,Inf)	# prop change in visc friction of link 1
-  29 db2      u	 [0,Inf)	# prop change in visc friction of link 2
-  30 db3      u	 [0,Inf)	# prop change in visc friction of link 3
-  31 db4      u	 [0,Inf)	# prop change in visc friction of link 4
-  32 db5      u	 [0,Inf)	# prop change in visc friction of link 5
-  33 thetadd6     u  (-Inf,Inf)	# ang acceleration of joint 6

GPy/util/datasets/pumadyn-32nm/accel/Prototask.spec

@@ -1,12 +0,0 @@
-#
-# Prototask.spec
-#
-Cases: all
-Origin: simulated 
-Inputs: 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32
-Order: retain
-Targets: 33
-Test-Set-Size: 4096
-Training-Set-Sizes: 64 128 256 512 1024
-Test-Set-Selection: hierarchical
-Maximum-Number-Of-Instances: 8

GPy/util/datasets/pumadyn-32nm/accel/std.prior

@@ -1,33 +0,0 @@
-  1   NLMH   real
-  2   NLMH   real
-  3   NLMH   real
-  4   NLMH   real
-  5   NLMH   real 
-  6   NLMH   real
-  7   NLMH   real
-  8   NLMH   real
-  9   NLMH   real
- 10   NLMH   real
- 11   NLMH   real
- 12   NLMH   real
- 13   NLMH   real
- 14   NLMH   real
- 15   NLMH   real
- 16   NLMH   real
- 17   NLMH   real
- 18   NLMH   real
- 19   NLMH   real
- 20   NLMH   real
- 21   NLMH   real
- 22   NLMH   real
- 23   NLMH   real
- 24   NLMH   real
- 25   NLMH   real
- 26   NLMH   real
- 27   NLMH   real
- 28   NLMH   real
- 29   NLMH   real
- 30   NLMH   real
- 31   NLMH   real
- 32   NLMH   real
- 33   NLMH   real

GPy/util/datasets/synth.tr

@@ -1,251 +0,0 @@
-    xs            ys       yc
-  0.05100797  0.16086164    0
- -0.74807425  0.08904024    0
- -0.77293371  0.26317168    0
-  0.21837360  0.12706142    0
-  0.37268336  0.49656200    0
- -0.62931544  0.63202159    0
- -0.43307167  0.14479166    0
- -0.84151970 -0.19131316    0
-  0.47525648  0.22483671    0
-  0.32082976  0.32721288    0
-  0.32061253  0.33407547    0
- -0.89077472  0.41168783    0
-  0.17850119  0.44691359    0
-  0.31558002  0.38853383    0
-  0.55777224  0.47272748    0
-  0.03191877  0.01222964    0
-  0.25090585  0.30716705    0
-  0.23571547  0.22493837    0
- -0.07236203  0.33376524    0
-  0.50440241  0.08054579    0
- -0.63223351  0.44552458    0
- -0.76784656  0.23614689    0
- -0.70017557  0.21038848    0
- -0.64713491  0.15921366    0
- -0.76739248  0.09259038    0
- -0.51788734  0.03288107    0
-  0.17516644  0.34534871    0
- -0.68031190  0.47612156    0
-  0.01595199  0.32167526    0
- -0.71481078  0.51421443    0
-  0.07837946  0.32284981    0
- -0.80872251  0.47036593    0
- -0.84211234  0.09294232    0
- -0.98591577  0.48309267    0
-  0.29104081  0.34275967    0
-  0.24321541  0.51488295    0
- -0.60104419  0.05060116    0
- -1.24652451  0.45923165    0
- -0.82769016  0.36187460    0
- -0.62117301 -0.10912158    0
- -0.70584105  0.65907662    0
-  0.06718867  0.60574850    0
-  0.30505147  0.47417973    0
-  0.60788138  0.39361588    0
- -0.78937483  0.17591675    0
- -0.53123209  0.42652809    0
-  0.25202071  0.17029707    0
- -0.57880357  0.26553665    0
- -0.83176749  0.54447377    0
- -0.69859164  0.38566851    0
- -0.73642607  0.11857527    0
- -0.93496195  0.11370707    0
-  0.43959309  0.41430638    0
- -0.54690854  0.24956276    0
- -0.08405550  0.36521058    0
-  0.32211458  0.69087105    0
-  0.10764739  0.57946932    0
- -0.71864030  0.25645757    0
- -0.87877752  0.45064757    0
- -0.69846046  0.95053870    0
-  0.39757434  0.11810207    0
- -0.50451354  0.57196376    0
-  0.25023622  0.39783889    0
-  0.61709156  0.10185808    0
-  0.31832860  0.08790562    0
- -0.57453363  0.18624195    0
-  0.09761865  0.55176786    0
-  0.48449339  0.35372973    0
-  0.52400684  0.46616851    0
- -0.78138463 -0.07534713    0
- -0.49704591  0.59948077    0
- -0.96984525  0.46624927    0
-  0.43541407  0.12192386    0
- -0.67942462  0.30753942    0
- -0.62529036  0.07099046    0
- -0.02318116  0.40442601    0
-  0.23200141  0.71066846    0
-  0.09384354  0.46674396    0
-  0.14234301  0.17898711    0
- -0.61686357  0.25507763    0
-  0.23636288  0.51543839    0
-  0.38914177  0.40429568    0
- -0.95178678 -0.03772239    0
-  0.24087822  0.71948890    0
-  0.12446266  0.45178849    0
- -0.60566430  0.26906478    0
- -0.71397188  0.30871780    0
-  0.31008428  0.34675335    0
-  0.18018786  0.46204643    0
- -0.42663885  0.64723225    0
-  0.06143230  0.32491150    0
-  0.07736952  0.32183287    0
-  0.42814970  0.13445957    0
- -0.80250753  0.66878999    0
-  0.40142623  0.42516398    0
-  0.37084776  0.26407123    0
- -0.80774748  0.41485899    0
-  0.50163585  0.23934856    0
-  0.58238323  0.22842741    0
- -0.59136100  0.30230321    0
- -0.87037236  0.26941446    0
- -0.72086765  0.19676678    0
-  0.27778443  0.21792253    0
-  0.33240813  0.27349865    0
- -0.14092068  0.39247351    0
- -0.59759518  0.14790267    0
- -0.85581534  0.14513961    0
- -0.88912232  0.26896001    0
-  0.21345680  0.43611756    0
- -0.53467949  0.57901229    0
-  0.31686848  0.39705856    0
- -0.68121733  0.04209840    0
- -0.97586127  0.45964811    0
-  0.41457183  0.27141230    0
-  0.32751292  0.36780137    0
- -0.93209192  0.09362034    0
-  0.58395341  0.47147282    0
- -0.44437309  0.23010142    0
-  0.29109441  0.19365556    0
- -0.51080722  0.41496003    0
- -0.96597511  0.17931052    0
-  0.18741315  0.29747132    0
-  0.17965417  0.45175449    0
- -0.72689602  0.35728387    0
- -0.54339877  0.41012013    0
- -0.59823393  0.98701425    1
- -0.20194736  0.62101680    1
-  0.47146103  0.48221146    1
- -0.09821987  0.58755577    1
- -0.35657658  0.63709705    1
-  0.63881392  0.42112135    1
-  0.62980614  0.28146085    1
- -0.46223286  0.61661031    1
- -0.07331555  0.55821736    1
- -0.55405533  0.51253129    1
- -0.43761773  0.87811781    1
- -0.22237814  0.88850773    1
-  0.09346162  0.67310494    1
-  0.53174745  0.54372650    1
-  0.40207539  0.51638462    1
-  0.47555171  0.65056336    1
- -0.23383266  0.63642580    1
- -0.31579316  0.75031340    1
- -0.47351720  0.63854125    1
-  0.59239464  0.89256953    1
- -0.22605324  0.79789454    1
- -0.43995011  0.52099256    1
- -0.54645044  0.74577198    1
-  0.46404306  0.51065152    1
- -0.15194296  0.81218439    1
-  0.48536395  0.82018093    1
-  0.34725649  0.70813773    1
-  0.43897015  0.62817158    1
- -0.21415914  0.64363951    1
-  0.57380231  0.63713466    1
-  0.38717361  0.58578395    1
-  0.32038322  0.53529127    1
- -0.20781491  0.65132467    1
- -0.18651283  0.81754816    1
-  0.24752692  0.39081936    1
-  0.66049881  0.89919213    1
- -0.28658801  0.73375946    1
- -0.32588080  0.39865509    1
- -0.25204565  0.67358326    1
-  0.37259022  0.49785904    1
- -0.29096564  1.04372060    1
- -0.30469807  0.86858292    1
- -0.21389978  1.09317811    1
- -0.36830015  0.75639546    1
- -0.46928218  0.88775091    1
-  0.39350146  0.77975197    1
- -0.45639966  0.80523454    1
-  0.51128242  0.76606136    1
-  0.22550468  0.46451215    1
-  0.01462984  0.40190926    1
- -0.19172785  0.80943313    1
-  0.38323479  0.75601744    1
-  0.49791612  0.61334375    1
-  0.35335230  0.77324337    1
- -0.34722575  0.70177856    1
-  0.58380468  0.76357539    1
- -0.13727764  0.71246351    1
-  0.38827268  0.44977123    1
- -0.53172709  0.61934293    1
- -0.11684624  0.87851210    1
-  0.54335864  0.41174865    1
- -0.45399302  0.66512988    1
- -0.21913200  0.83484947    1
-  0.30485742  0.98028760    1
-  0.65676798  0.75766017    1
-  0.61420447  0.75039019    1
- -0.45809964  0.77968606    1
- -0.21617465  0.88626305    1
- -0.26016108  0.81008591    1
-  0.31884531  0.84517725    1
- -0.23727415  0.80178784    1
-  0.58310323  0.77709806    1
-  0.02841337  0.75792620    1
- -0.41840136  0.68041440    1
-  0.67412880  0.60245461    1
- -0.25278281  0.70526103    1
-  0.51609843  0.62092390    1
-  0.20392294  0.91641482    1
- -0.17207124  1.00884096    1
-  0.27274507  0.29346977    1
-  0.07634798  0.56222204    1
- -0.36653499  0.64831007    1
-  0.44290673  0.80087721    1
- -0.19976385  0.54295162    1
- -0.54075738  0.65293033    1
- -0.07060266  1.00296912    1
-  0.50715054  0.35045758    1
- -0.06048611  0.62982713    1
-  0.21532928  0.60260249    1
-  0.46809108  0.87182416    1
- -0.29888511  0.73669866    1
-  0.86129620  0.47289330    1
-  0.70120877  0.74572893    1
- -0.11342797  0.60067099    1
-  0.31234354  0.90756345    1
- -0.12172541  0.84112851    1
-  0.36867857  0.37052586    1
-  0.57311489  0.40949740    1
- -0.25841225  0.67192335    1
-  0.30937186  0.50823318    1
-  0.43319338  0.77016967    1
- -0.30448035  0.57820106    1
-  0.44276338  0.58023403    1
- -0.19442057  0.89876808    1
- -0.06105237  0.74184946    1
-  0.07619347  0.35386246    1
-  0.85826993  0.95819523    1
-  0.37039200  0.72342401    1
-  0.51481515  0.76203996    1
-  0.43127521  0.54259166    1
-  0.42286091  0.65242185    1
-  0.29815001  0.93453682    1
-  0.37128253  0.70089181    1
- -0.51528729  0.76473490    1
-  0.38525783  0.65528189    1
- -0.34825368  0.50529981    1
-  0.68510504  0.78067440    1
- -0.36528923  0.45703265    1
- -0.40903577  0.74230433    1
-  0.43574387  0.44689789    1
-  0.26887846  0.44559230    1
- -0.49254862  1.01443372    1
-  0.07615960  0.63795180    1
-  0.49226224  0.46876241    1
- -0.40249641  0.71301084    1

GPy/util/mocap.py

@@ -692,84 +692,3 @@ skel = acclaim_skeleton()
 
 
 
-    
-def fetch_cmu(subj_motions, base_url = 'http://mocap.cs.cmu.edu:8080/subjects', skel_store_dir = '.', motion_store_dir = '.', store_motions = True, return_motions = True, messages = True):
-    ''' 
-    Download and store the skel. and motions indicated in a tuple (A,B) where A is a list of skeletons and B
-    the corresponding 2-D list of motions, ie B_ij is the j-th motion to download for skeleton A_i
-    The method can optionally store the fetched data and / or return them as arrays.
-    If the data are already stored, they are not fetched but just retrieved.
-
-    e.g.
-    # Download the data, do not return anything
-    GPy.util.mocap.fetch_cmu(subj_motions = ([35],[[1,2,3]]), return_motions = False)
-    # Fetch and return the data in a list. Do not store them anywhere
-    GPy.util.mocap.fetch_cmu(subj_motions = ([35],[[1,2,3]]), return_motions = True, store_motions = False)
-
-    In both cases above, if the data do exist in the given skel_store_dir and motion_store_dir, they are just loaded from there.
-    '''
-    
-    subjectsNum = subj_motions[0]
-    motionsNum = subj_motions[1]
-
-    # Convert numbers to strings
-    subjects = []
-    motions = [list() for _ in range(len(subjectsNum))]
-    for i in range(len(subjectsNum)):
-        curSubj = str(int(subjectsNum[i]))
-        if subjectsNum[i] < 10:
-            curSubj = '0' + curSubj
-        subjects.append(curSubj)
-        for j in range(len(motionsNum[i])):
-            curMot = str(int(motionsNum[i][j]))
-            if motionsNum[i][j] < 10:
-                curMot = '0' + curMot
-            motions[i].append(curMot)
-
-
-    all_skels = []
-    
-    assert len(subjects) == len(motions)
-    
-    if return_motions:
-        all_motions = [list() for _ in range(len(subjects))]
-    else:
-        all_motions = []
-            
-    for i in range(len(subjects)):
-        cur_skel_suffix = '/' + subjects[i] + '/' 
-        cur_skel_dir = skel_store_dir + cur_skel_suffix
-        cur_skel_file = cur_skel_dir + subjects[i] + '.asf'
-        cur_skel_url = base_url + cur_skel_suffix + subjects[i] + '.asf'
-        
-        if os.path.isfile(cur_skel_file):
-            if return_motions:
-                with open(cur_skel_file, 'r') as f:
-                    cur_skel_data = f.read()
-        else:
-            if store_motions:
-                if not os.path.isdir(cur_skel_dir):
-                    os.mkdir(cur_skel_dir)
-                if not os.path.isdir(motion_store_dir + cur_skel_suffix):
-                    os.mkdir(motion_store_dir + cur_skel_suffix)
-            cur_skel_data = dat.download_resource(cur_skel_url, cur_skel_file, store_motions, messages)
-        
-        if return_motions:
-            all_skels.append(cur_skel_data)
-        
-        for j in range(len(motions[i])):
-            cur_motion_url = base_url + cur_skel_suffix + subjects[i] + '_' + motions[i][j] + '.amc'
-            cur_motion_file = motion_store_dir + cur_skel_suffix  + subjects[i] + '_' + motions[i][j] + '.amc'
-            if os.path.isfile(cur_motion_file):
-                with open(cur_motion_file, 'r') as f:
-                    if return_motions:
-                        cur_motion_data = f.read()
-            else:
-                cur_motion_data = dat.download_resource(cur_motion_url, cur_motion_file, store_motions, messages)
-
-            if return_motions:
-                all_motions[i].append(cur_motion_data)
-
-    return all_skels, all_motions
-
-