From 0509b6f9d04bf4bb6d73f7e3120d0d3d39d4eaa0 Mon Sep 17 00:00:00 2001 From: Neil Lawrence Date: Tue, 24 Sep 2013 06:15:09 +0200 Subject: [PATCH] Added missing files. --- GPy/FAQ.txt | 8 ++ GPy/coding_style_guide.txt | 10 ++ GPy/kern/parts/eq_ode1.py | 217 +++++++++++++++-------------------- GPy/testing/bcgplvm_tests.py | 50 ++++++++ GPy/util/erfcx.py | 63 ++++++++++ GPy/util/ln_diff_erfs.py | 109 ++++++++++++------ 6 files changed, 297 insertions(+), 160 deletions(-) create mode 100644 GPy/FAQ.txt create mode 100644 GPy/coding_style_guide.txt create mode 100644 GPy/testing/bcgplvm_tests.py create mode 100644 GPy/util/erfcx.py diff --git a/GPy/FAQ.txt b/GPy/FAQ.txt new file mode 100644 index 00000000..66ba4834 --- /dev/null +++ b/GPy/FAQ.txt @@ -0,0 +1,8 @@ +Frequently Asked Questions +-------------------------- + +Unit tests are run through Travis-Ci. They can be run locally through entering the GPy route diretory and writing + +nosetests testing/ + +Documentation is handled by Sphinx. To build the documentation: diff --git a/GPy/coding_style_guide.txt b/GPy/coding_style_guide.txt new file mode 100644 index 00000000..0cc732e4 --- /dev/null +++ b/GPy/coding_style_guide.txt @@ -0,0 +1,10 @@ +In this text document we will describe coding conventions to be used in GPy to keep things consistent. + +All arrays containing data are two dimensional. The first dimension is the number of data, the second dimension is number of features. This keeps things consistent with the idea of a design matrix. + +Input matrices are either X or t, output matrices are Y. + +Input dimensionality is input_dim, output dimensionality is output_dim, number of data is num_data. + +Data sets are preprocessed in the datasets.py file. This file also records where the data set was obtained from in the dictionary stored in the file. Long term we should move this dictionary to sqlite or similar. + diff --git a/GPy/kern/parts/eq_ode1.py b/GPy/kern/parts/eq_ode1.py index 1e5345fe..6fd5fb91 100644 --- a/GPy/kern/parts/eq_ode1.py +++ b/GPy/kern/parts/eq_ode1.py @@ -39,11 +39,12 @@ class Eq_ode1(Kernpart): self.name = 'eq_ode1' self.output_dim = output_dim self.lengthscale = lengthscale - self.num_params = self.output_dim*(1. + self.rank) + 1 + self.num_params = self.output_dim*self.rank + 1 + (self.output_dim - 1) if kappa is not None: self.num_params+=self.output_dim if delay is not None: - self.num_params+=self.output_dim + assert delay.shape==(self.output_dim-1,) + self.num_params+=self.output_dim-1 self.rank = rank if W is None: self.W = 0.5*np.random.randn(self.output_dim,self.rank)/np.sqrt(self.rank) @@ -51,18 +52,17 @@ class Eq_ode1(Kernpart): assert W.shape==(self.output_dim,self.rank) self.W = W if decay is None: - self.decay = np.ones(self.output_dim) + self.decay = np.ones(self.output_dim-1) if kappa is not None: assert kappa.shape==(self.output_dim,) self.kappa = kappa - if delay is not None: - assert delay.shape==(self.output_dim,) self.delay = delay self.is_normalized = True self.is_stationary = False + self.gaussian_initial = False self._set_params(self._get_params()) - + def _get_params(self): param_list = [self.W.flatten()] if self.kappa is not None: @@ -84,11 +84,11 @@ class Eq_ode1(Kernpart): self.kappa = x[start:end] self.B += np.diag(self.kappa) start=end - end+=self.output_dim + end+=self.output_dim-1 self.decay = x[start:end] start=end if self.delay is not None: - end+=self.output_dim + end+=self.output_dim-1 self.delay = x[start:end] start=end end+=1 @@ -100,9 +100,9 @@ class Eq_ode1(Kernpart): param_names = sum([['W%i_%i'%(i,j) for j in range(self.rank)] for i in range(self.output_dim)],[]) if self.kappa is not None: param_names += ['kappa_%i'%i for i in range(self.output_dim)] - param_names += ['decay_%i'%i for i in range(self.output_dim)] + param_names += ['decay_%i'%i for i in range(1,self.output_dim)] if self.delay is not None: - param_names += ['delay_%i'%i for i in range(self.output_dim)] + param_names += ['delay_%i'%i for i in 1+range(1,self.output_dim)] param_names+= ['lengthscale'] return param_names @@ -112,7 +112,7 @@ class Eq_ode1(Kernpart): raise ValueError('Input matrix for ode1 covariance should have at most two columns, one containing times, the other output indices') self._K_computations(X, X2) - target += self._scales*self._dK_dvar + target += self._scale*self._K_dvar if self.gaussian_initial: # Add covariance associated with initial condition. @@ -140,9 +140,8 @@ class Eq_ode1(Kernpart): # TODO: some fast checking here to see if this needs recomputing? self._t = X[:, 0] - if X.shape[1]==1: - # No index passed, assume single output of ode model. - self._index = np.ones_like(X[:, 1],dtype=np.int) + if not X.shape[1] == 2: + raise ValueError('Input matrix for ode1 covariance should have two columns, one containing times, the other output indices') self._index = np.asarray(X[:, 1],dtype=np.int) # Sort indices so that outputs are in blocks for computational # convenience. @@ -156,15 +155,12 @@ class Eq_ode1(Kernpart): self._index2 = None self._rorder2 = self._rorder else: - if X2.shape[1] > 2: - raise ValueError('Input matrix for ode1 covariance should have at most two columns, one containing times, the other output indices') + if not X2.shape[1] == 2: + raise ValueError('Input matrix for ode1 covariance should have two columns, one containing times, the other output indices') self._t2 = X2[:, 0] - if X.shape[1]==1: - # No index passed, assume single output of ode model. - self._index2 = np.ones_like(X2[:, 1],dtype=np.int) self._index2 = np.asarray(X2[:, 1],dtype=np.int) self._order2 = self._index2.argsort() - slef._index2 = self._index2[self._order2] + self._index2 = self._index2[self._order2] self._t2 = self._t2[self._order2] self._rorder2 = self._order2.argsort() # rorder2 is for reversing order @@ -180,25 +176,65 @@ class Eq_ode1(Kernpart): else: self._K_compute_ode_eq(transpose=True) self._K_compute_ode() - + + if X2 is None: + self._K_dvar = np.zeros((self._t.shape[0], self._t.shape[0])) + else: + self._K_dvar = np.zeros((self._t.shape[0], self._t2.shape[0])) + # Reorder values of blocks for placing back into _K_dvar. - self._K_dvar[self._rorder, :] = np.vstack((np.hstack((self._K_eq, self._Keq_ode)), - np.hstack((self._K_ode_eq, self.K_ode))))[:, self._rorder2] + self._K_dvar = np.vstack((np.hstack((self._K_eq, self._K_eq_ode)), + np.hstack((self._K_ode_eq, self._K_ode)))) + self._K_dvar = self._K_dvar[self._rorder, :] + self._K_dvar = self._K_dvar[:, self._rorder2] + + + if X2 is None: + # Matrix giving scales of each output + self._scale = np.zeros((self._t.size, self._t.size)) + code=""" + for(int i=0;i0] index_ode = self._index2[self._index2>0]-1 - if t_ode.shape[0]==0: - self._K_eq_ode = np.zeros((0, 0)) - return else: - self._K_eq_ode = np.zeros((0, 0)) - return - - t_eq = self._t[self._index==0] - if t_eq.shape[0]==0: - self._K_eq_ode = np.zeros((0, 0)) - return + t_eq = self._t2[self._index2==0] + t_ode = self._t[self._index>0] + index_ode = self._index[self._index>0]-1 else: + t_eq = self._t[self._index==0] t_ode = self._t[self._index>0] index_ode = self._index[self._index>0]-1 - if t_ode.shape[0]==0: - self._K_ode_eq = np.zeros((0, 0)) - return - if self._t2 is not None: - t_eq = self._t2[self._index2==0] - if t_eq.shape[0]==0: - self._K_ode_eq = np.zeros((0, 0)) - return + + if t_ode.size==0 or t_eq.size==0: + if transpose: + self._K_eq_ode = np.zeros((t_eq.shape[0], t_ode.shape[0])) else: - self._K_ode_eq = np.zeros((0, 0)) - return - # Matrix giving scales of each output - # self._scale = np.zeros((t_ode.shape[0], t_eq.shape[0])) - # code=""" - # for(int i=0;i0] index_ode = self._index[self._index>0]-1 - if t_ode.shape[0]==0: - self._K_ode = np.zeros((0, 0)) - return - if self._t2 is None: + if t_ode.size==0: + self._K_ode = np.zeros((0, 0)) + return t2_ode = t_ode index2_ode = index_ode else: t2_ode = self._t2[self._index2>0] - index2_ode = self._index2[self._index2>0]-1 - if t2_eq.shape[0]==0: - self._K_ode = np.zeros((0, 0)) + if t_ode.size==0 or t2_ode.size==0: + self._K_ode = np.zeros((t_ode.size, t2_ode.size)) return - - if self._index2 is None: - # Matrix giving scales of each output - self._scale = np.zeros((t_ode.shape[0], t_ode.shape[0])) - code=""" - for(int i=0;i0]-1 # When index is identical if self.is_stationary: @@ -360,10 +327,10 @@ class Eq_ode1(Kernpart): h2 = self._compute_H_stat(t2_ode, index2_ode, t_ode, index_ode) else: h2 = self._compute_H(t2_ode, index2_ode, t_ode, index_ode) - self._K_ode += 0.5 * (h + h2.T) + self._K_ode = 0.5 * (h + h2.T) if not self.is_normalized: - self._K_ode *= np.sqrt(np.pi)*sigma + self._K_ode *= np.sqrt(np.pi)*self.sigma def _compute_H(self, t, index, t2, index2, update_derivatives=False): """Helper function for computing part of the ode1 covariance function. @@ -441,7 +408,7 @@ class Eq_ode1(Kernpart): -Decay[:, None]*t_mat-Decay2[None, :]*t2_mat+ln_part_2 -np.log(Decay[:, None] + Decay2[None, :])) - + return h # if update_derivatives: # sigma2 = self.sigma*self.sigma # # Update ith decay gradient diff --git a/GPy/testing/bcgplvm_tests.py b/GPy/testing/bcgplvm_tests.py new file mode 100644 index 00000000..94282a0b --- /dev/null +++ b/GPy/testing/bcgplvm_tests.py @@ -0,0 +1,50 @@ +# Copyright (c) 2013, GPy authors (see AUTHORS.txt) +# Licensed under the BSD 3-clause license (see LICENSE.txt) + +import unittest +import numpy as np +import GPy + +class BCGPLVMTests(unittest.TestCase): + def test_kernel_backconstraint(self): + num_data, num_inducing, input_dim, output_dim = 10, 3, 2, 4 + X = np.random.rand(num_data, input_dim) + k = GPy.kern.rbf(input_dim) + GPy.kern.white(input_dim, 0.00001) + K = k.K(X) + Y = np.random.multivariate_normal(np.zeros(num_data),K,output_dim).T + k = GPy.kern.mlp(input_dim) + GPy.kern.bias(input_dim) + bk = GPy.kern.rbf(output_dim) + mapping = GPy.mappings.Kernel(output_dim=input_dim, X=Y, kernel=bk) + m = GPy.models.BCGPLVM(Y, input_dim, kernel = k, mapping=mapping) + m.randomize() + self.assertTrue(m.checkgrad()) + + def test_linear_backconstraint(self): + num_data, num_inducing, input_dim, output_dim = 10, 3, 2, 4 + X = np.random.rand(num_data, input_dim) + k = GPy.kern.rbf(input_dim) + GPy.kern.white(input_dim, 0.00001) + K = k.K(X) + Y = np.random.multivariate_normal(np.zeros(num_data),K,output_dim).T + k = GPy.kern.mlp(input_dim) + GPy.kern.bias(input_dim) + bk = GPy.kern.rbf(output_dim) + mapping = GPy.mappings.Linear(output_dim=input_dim, input_dim=output_dim) + m = GPy.models.BCGPLVM(Y, input_dim, kernel = k, mapping=mapping) + m.randomize() + self.assertTrue(m.checkgrad()) + + def test_mlp_backconstraint(self): + num_data, num_inducing, input_dim, output_dim = 10, 3, 2, 4 + X = np.random.rand(num_data, input_dim) + k = GPy.kern.rbf(input_dim) + GPy.kern.white(input_dim, 0.00001) + K = k.K(X) + Y = np.random.multivariate_normal(np.zeros(num_data),K,output_dim).T + k = GPy.kern.mlp(input_dim) + GPy.kern.bias(input_dim) + bk = GPy.kern.rbf(output_dim) + mapping = GPy.mappings.MLP(output_dim=input_dim, input_dim=output_dim, hidden_dim=[5, 4, 7]) + m = GPy.models.BCGPLVM(Y, input_dim, kernel = k, mapping=mapping) + m.randomize() + self.assertTrue(m.checkgrad()) + +if __name__ == "__main__": + print "Running unit tests, please be (very) patient..." + unittest.main() diff --git a/GPy/util/erfcx.py b/GPy/util/erfcx.py new file mode 100644 index 00000000..f42e49f3 --- /dev/null +++ b/GPy/util/erfcx.py @@ -0,0 +1,63 @@ +## Copyright (C) 2010 Soren Hauberg +## +## Copyright James Hensman 2011 +## +## This program is free software; you can redistribute it and/or modify it +## under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 3 of the License, or (at +## your option) any later version. +## +## This program is distributed in the hope that it will be useful, but +## WITHOUT ANY WARRANTY; without even the implied warranty of +## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +## General Public License for more details. +## +## You should have received a copy of the GNU General Public License +## along with this program; see the file COPYING. If not, see +## . + +import numpy as np + +def erfcx (arg): + arg = np.atleast_1d(arg) + assert(np.all(np.isreal(arg)),"erfcx: input must be real") + + ## Get precision dependent thresholds -- or not :p + xneg = -26.628; + xmax = 2.53e+307; + + ## Allocate output + result = np.zeros (arg.shape) + + ## Find values where erfcx can be evaluated + idx_neg = (arg < xneg); + idx_max = (arg > xmax); + idx = ~(idx_neg | idx_max); + + arg = arg [idx]; + + ## Perform the actual computation + t = 3.97886080735226 / (np.abs (arg) + 3.97886080735226); + u = t - 0.5; + y = (((((((((u * 0.00127109764952614092 + 1.19314022838340944e-4) * u \ + - 0.003963850973605135) * u - 8.70779635317295828e-4) * u + \ + 0.00773672528313526668) * u + 0.00383335126264887303) * u - \ + 0.0127223813782122755) * u - 0.0133823644533460069) * u + \ + 0.0161315329733252248) * u + 0.0390976845588484035) * u + \ + 0.00249367200053503304; + y = ((((((((((((y * u - 0.0838864557023001992) * u - \ + 0.119463959964325415) * u + 0.0166207924969367356) * u + \ + 0.357524274449531043) * u + 0.805276408752910567) * u + \ + 1.18902982909273333) * u + 1.37040217682338167) * u + \ + 1.31314653831023098) * u + 1.07925515155856677) * u + \ + 0.774368199119538609) * u + 0.490165080585318424) * u + \ + 0.275374741597376782) * t; + + y [arg < 0] = 2 * np.exp (arg [arg < 0]**2) - y [arg < 0]; + + ## Put the results back into something with the same size is the original input + result [idx] = y; + result [idx_neg] = np.inf; + ## result (idx_max) = 0; # not needed as we initialise with zeros + return(result) + diff --git a/GPy/util/ln_diff_erfs.py b/GPy/util/ln_diff_erfs.py index 13b442f7..bb9cfe03 100644 --- a/GPy/util/ln_diff_erfs.py +++ b/GPy/util/ln_diff_erfs.py @@ -18,54 +18,93 @@ def ln_diff_erfs(x1, x2, return_sign=False): :type x2: ndarray :return: tuple containing (log(abs(erf(x1) - erf(x2))), sign(erf(x1) - erf(x2))) - Based on MATLAB code that was written by Antti Honkela and modified by David Luengo, originally derived from code by Neil Lawrence. + Based on MATLAB code that was written by Antti Honkela and modified by David Luengo and originally derived from code by Neil Lawrence. """ x1 = np.require(x1).real x2 = np.require(x2).real - - v = np.zeros(np.max((x1.size, x2.size))) + if x1.size==1: + x1 = np.reshape(x1, (1, 1)) + if x2.size==1: + x2 = np.reshape(x2, (1, 1)) + + if x1.shape==x2.shape: + v = np.zeros_like(x1) + else: + if x1.size==1: + v = np.zeros(x2.shape) + elif x2.size==1: + v = np.zeros(x1.shape) + else: + raise ValueError, "This function does not broadcast unless provided with a scalar." - # if numel(x1) == 1: - # x1 = x1 * ones(size(x2)) + if x1.size == 1: + x1 = np.tile(x1, x2.shape) - # if numel(x2) == 1: - # x2 = x2 * ones(size(x1)) + if x2.size == 1: + x2 = np.tile(x2, x1.shape) sign = np.sign(x1 - x2) - I = sign == -1 - swap = x1[I] - x1[I] = x2[I] - x2[I] = swap + if x1.size == 1: + if sign== -1: + swap = x1 + x1 = x2 + x2 = swap + else: + I = sign == -1 + swap = x1[I] + x1[I] = x2[I] + x2[I] = swap - # TODO: switch off log of zero warnings. - # Case 1: arguments of different sign, no problems with loss of accuracy - I1 = np.logical_or(np.logical_and(x1>0, x2<0), np.logical_and(x2>0, x1<0)) # I1=(x1*x2)<0 - v[I1] = np.log( erf(x1[I1]) - erf(x2[I1]) ) + with np.errstate(divide='ignore'): + # switch off log of zero warnings. - # Case 2: x1 = x2 so we have log of zero. - I2 = (x1 == x2) - v[I2] = -np.inf + # Case 0: arguments of different sign, no problems with loss of accuracy + I0 = np.logical_or(np.logical_and(x1>0, x2<0), np.logical_and(x2>0, x1<0)) # I1=(x1*x2)<0 - # Case 3: Both arguments are non-negative - I3 = np.logical_and(x1 > 0, np.logical_and(np.logical_not(I1), - np.logical_not(I2))) - v[I3] = np.log(erfcx(x2[I3]) - -erfcx(x1[I3])*np.exp(x2[I3]**2 - -x1[I3]**2)) - x2[I3]**2 - - # Case 4: Both arguments are non-positive - I4 = np.logical_and(np.logical_and(np.logical_not(I1), - np.logical_not(I2)), - np.logical_not(I3)) - v[I4] = np.log(erfcx(-x1[I4]) - -erfcx(-x2[I4])*np.exp(x1[I4]**2 - -x2[I4]**2))-x1[I4]**2 - + # Case 1: x1 = x2 so we have log of zero. + I1 = (x1 == x2) + + # Case 2: Both arguments are non-negative + I2 = np.logical_and(x1 > 0, np.logical_and(np.logical_not(I0), + np.logical_not(I1))) + # Case 3: Both arguments are non-positive + I3 = np.logical_and(np.logical_and(np.logical_not(I0), + np.logical_not(I1)), + np.logical_not(I2)) + _x2 = x2.flatten() + _x1 = x1.flatten() + for group, flags in zip((0, 1, 2, 3), (I0, I1, I2, I3)): + + if np.any(flags): + if not x1.size==1: + _x1 = x1[flags] + if not x2.size==1: + _x2 = x2[flags] + if group==0: + v[flags] = np.log( erf(_x1) - erf(_x2) ) + elif group==1: + v[flags] = -np.inf + elif group==2: + v[flags] = np.log(erfcx(_x2) + -erfcx(_x1)*np.exp(_x2**2 + -_x1**2)) - _x2**2 + elif group==3: + v[flags] = np.log(erfcx(-_x1) + -erfcx(-_x2)*np.exp(_x1**2 + -_x2**2))-_x1**2 + # TODO: switch back on log of zero warnings. if return_sign: return v, sign else: - # Need to add in a complex part because argument is negative. - v[I] += np.pi*1j + if v.size==1: + if sign==-1: + v = v.view('complex64') + v += np.pi*1j + else: + # Need to add in a complex part because argument is negative. + v = v.view('complex64') + v[I] += np.pi*1j + return v