apunkt/GPy
1
0
Fork 0
mirror of https://github.com/SheffieldML/GPy.git synced 2026-06-05 14:55:15 +02:00
Code Issues Projects Releases Packages Wiki Activity Actions

Integrated Laplace and merged Merge remote-tracking branch 'gpy_real/devel' into merge_branch

Browse source 
Conflicts:
	GPy/core/gp.py
	GPy/likelihoods/__init__.py
	GPy/likelihoods/likelihood_functions.py
	GPy/likelihoods/link_functions.py
This commit is contained in:
Alan Saul 2013-10-03 16:52:02 +01:00
commit 8343615098
106 changed files with 5841 additions and 1134 deletions
Download patch file Download diff file Expand all files Collapse all files

38
GPy/util/datasets.py
View file

@ -7,9 +7,21 @@ import urllib as url
import zipfile
import tarfile
import datetime
ipython_notebook = False
if ipython_notebook:
import IPython.core.display
def ipynb_input(varname, prompt=''):
"""Prompt user for input and assign string val to given variable name."""
js_code = ("""
var value = prompt("{prompt}","");
var py_code = "{varname} = '" + value + "'";
IPython.notebook.kernel.execute(py_code);
""").format(prompt=prompt, varname=varname)
return IPython.core.display.Javascript(js_code)
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),
@ -130,14 +142,18 @@ The database was created with funding from NSF EIA-0196217.""",
'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'])
choice = raw_input().lower()
choice = ''
if ipython_notebook:
ipynb_input(choice, prompt='provide your answer here')
else:
choice = raw_input().lower()
if choice in yes:
return True
elif choice in no:
@ -146,6 +162,7 @@ def prompt_user():
sys.stdout.write("Please respond with 'yes', 'y' or 'no', 'n'")
return prompt_user()
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']:
@ -524,11 +541,14 @@ def simulation_BGPLVM():
'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.
"""
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)
num_in = 1
@ -631,11 +651,15 @@ def olympic_marathon_men(data_set='olympic_marathon_men'):
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.
"""
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).
:type num_data: int
:param seed: random seed to be used for data generation.
:type seed: int"""
:type seed: int
"""
np.random.seed(seed=seed)
sqrt2 = np.sqrt(2)
# Rotation matrix

63
GPy/util/erfcx.py Normal file
View file

@ -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
## <http://www.gnu.org/licenses/>.
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)

58
GPy/util/linalg.py
View file

@ -27,31 +27,37 @@ except:
_blas_available = False
def dtrtrs(A, B, lower=0, trans=0, unitdiag=0):
"""Wrapper for lapack dtrtrs function
"""
Wrapper for lapack dtrtrs function
:param A: Matrix A
:param B: Matrix B
:param lower: is matrix lower (true) or upper (false)
:returns:
"""
return lapack.dtrtrs(A, B, lower=lower, trans=trans, unitdiag=unitdiag)
def dpotrs(A, B, lower=0):
"""Wrapper for lapack dpotrs function
"""
Wrapper for lapack dpotrs function
:param A: Matrix A
:param B: Matrix B
:param lower: is matrix lower (true) or upper (false)
:returns:
"""
return lapack.dpotrs(A, B, lower=lower)
def dpotri(A, lower=0):
"""Wrapper for lapack dpotri function
"""
Wrapper for lapack dpotri function
:param A: Matrix A
:param lower: is matrix lower (true) or upper (false)
:returns:
:returns: A inverse
"""
return lapack.dpotri(A, lower=lower)
@ -65,18 +71,20 @@ def pddet(A):
def trace_dot(a, b):
"""
efficiently compute the trace of the matrix product of a and b
Efficiently compute the trace of the matrix product of a and b
"""
return np.sum(a * b)
def mdot(*args):
"""Multiply all the arguments using matrix product rules.
"""
Multiply all the arguments using matrix product rules.
The output is equivalent to multiplying the arguments one by one
from left to right using dot().
Precedence can be controlled by creating tuples of arguments,
for instance mdot(a,((b,c),d)) multiplies a (a*((b*c)*d)).
Note that this means the output of dot(a,b) and mdot(a,b) will differ if
a or b is a pure tuple of numbers.
"""
if len(args) == 1:
return args[0]
@ -123,14 +131,16 @@ def jitchol(A, maxtries=5):
def jitchol_old(A, maxtries=5):
"""
:param A : An almost pd square matrix
:param A: An almost pd square matrix
:rval L: the Cholesky decomposition of A
.. Note:
.. note:
Adds jitter to K, to enforce positive-definiteness
if stuff breaks, please check:
np.allclose(sp.linalg.cholesky(XXT, lower = True), np.triu(sp.linalg.cho_factor(XXT)[0]).T)
"""
try:
return linalg.cholesky(A, lower=True)
@ -150,6 +160,7 @@ def jitchol_old(A, maxtries=5):
def pdinv(A, *args):
"""
:param A: A DxD pd numpy array
:rval Ai: the inverse of A
@ -160,6 +171,7 @@ def pdinv(A, *args):
:rtype Li: np.ndarray
:rval logdet: the log of the determinant of A
:rtype logdet: float64
"""
L = jitchol(A, *args)
logdet = 2.*np.sum(np.log(np.diag(L)))
@ -185,14 +197,13 @@ def chol_inv(L):
def multiple_pdinv(A):
"""
Arguments
---------
:param A: A DxDxN numpy array (each A[:,:,i] is pd)
Returns
-------
invs : the inverses of A
hld: 0.5* the log of the determinants of A
:rval invs: the inverses of A
:rtype invs: np.ndarray
:rval hld: 0.5* the log of the determinants of A
:rtype hld: np.array
"""
N = A.shape[-1]
chols = [jitchol(A[:, :, i]) for i in range(N)]
@ -206,15 +217,13 @@ def PCA(Y, input_dim):
"""
Principal component analysis: maximum likelihood solution by SVD
Arguments
---------
:param Y: NxD np.array of data
:param input_dim: int, dimension of projection
Returns
-------
:rval X: - Nxinput_dim np.array of dimensionality reduced data
W - input_dimxD mapping from X to Y
:rval W: - input_dimxD mapping from X to Y
"""
if not np.allclose(Y.mean(axis=0), 0.0):
print "Y is not zero mean, centering it locally (GPy.util.linalg.PCA)"
@ -281,11 +290,10 @@ def DSYR_blas(A, x, alpha=1.):
Performs a symmetric rank-1 update operation:
A <- A + alpha * np.dot(x,x.T)
Arguments
---------
:param A: Symmetric NxN np.array
:param x: Nx1 np.array
:param alpha: scalar
"""
N = c_int(A.shape[0])
LDA = c_int(A.shape[0])
@ -303,11 +311,10 @@ def DSYR_numpy(A, x, alpha=1.):
Performs a symmetric rank-1 update operation:
A <- A + alpha * np.dot(x,x.T)
Arguments
---------
:param A: Symmetric NxN np.array
:param x: Nx1 np.array
:param alpha: scalar
"""
A += alpha * np.dot(x[:, None], x[None, :])
@ -371,8 +378,9 @@ def cholupdate(L, x):
"""
update the LOWER cholesky factor of a pd matrix IN PLACE
if L is the lower chol. of K, then this function computes L_
where L_ is the lower chol of K + x*x^T
if L is the lower chol. of K, then this function computes L\_
where L\_ is the lower chol of K + x*x^T
"""
support_code = """
#include <math.h>

110
GPy/util/ln_diff_erfs.py Normal file
View file

@ -0,0 +1,110 @@
# Copyright (c) 2013, GPy authors (see AUTHORS.txt).
# Licensed under the BSD 3-clause license (see LICENSE.txt)
#Only works for scipy 0.12+
try:
from scipy.special import erfcx, erf
except ImportError:
from scipy.special import erf
from erfcx import erfcx
import numpy as np
def ln_diff_erfs(x1, x2, return_sign=False):
"""Function for stably computing the log of difference of two erfs in a numerically stable manner.
:param x1 : argument of the positive erf
:type x1: ndarray
:param x2 : argument of the negative erf
: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 and originally derived from code by Neil Lawrence.
"""
x1 = np.require(x1).real
x2 = np.require(x2).real
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 x1.size == 1:
x1 = np.tile(x1, x2.shape)
if x2.size == 1:
x2 = np.tile(x2, x1.shape)
sign = np.sign(x1 - x2)
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
with np.errstate(divide='ignore'):
# switch off log of zero warnings.
# 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 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:
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

10
GPy/util/misc.py
View file

@ -17,12 +17,9 @@ def linear_grid(D, n = 100, min_max = (-100, 100)):
"""
Creates a D-dimensional grid of n linearly spaced points
Parameters:
D: dimension of the grid
n: number of points
min_max: (min, max) list
:param D: dimension of the grid
:param n: number of points
:param min_max: (min, max) list
"""
@ -39,6 +36,7 @@ def kmm_init(X, m = 10):
:param X: data
:param m: number of inducing points
"""
# compute the distances

39
GPy/util/mocap.py
View file

@ -92,13 +92,15 @@ class tree:
def swap_vertices(self, i, j):
"""Swap two vertices in the tree structure array.
"""
Swap two vertices in the tree structure array.
swap_vertex swaps the location of two vertices in a tree structure array.
ARG tree : the tree for which two vertices are to be swapped.
ARG i : the index of the first vertex to be swapped.
ARG j : the index of the second vertex to be swapped.
RETURN tree : the tree structure with the two vertex locations
swapped.
:param tree: the tree for which two vertices are to be swapped.
:param i: the index of the first vertex to be swapped.
:param j: the index of the second vertex to be swapped.
:rval tree: the tree structure with the two vertex locations swapped.
"""
store_vertex_i = self.vertices[i]
store_vertex_j = self.vertices[j]
@ -117,12 +119,17 @@ class tree:
def rotation_matrix(xangle, yangle, zangle, order='zxy', degrees=False):
"""Compute the rotation matrix for an angle in each direction.
"""
Compute the rotation matrix for an angle in each direction.
This is a helper function for computing the rotation matrix for a given set of angles in a given order.
ARG xangle : rotation for x-axis.
ARG yangle : rotation for y-axis.
ARG zangle : rotation for z-axis.
ARG order : the order for the rotations."""
:param xangle: rotation for x-axis.
:param yangle: rotation for y-axis.
:param zangle: rotation for z-axis.
:param order: the order for the rotations.
"""
if degrees:
xangle = math.radians(xangle)
yangle = math.radians(yangle)
@ -301,10 +308,12 @@ class acclaim_skeleton(skeleton):
def load_skel(self, file_name):
"""Loads an ASF file into a skeleton structure.
loads skeleton structure from an acclaim skeleton file.
ARG file_name : the file name to load in.
RETURN skel : the skeleton for the file."""
"""
Loads an ASF file into a skeleton structure.
:param file_name: The file name to load in.
"""
fid = open(file_name, 'r')
self.read_skel(fid)

35
GPy/util/multioutput.py Normal file
View file

@ -0,0 +1,35 @@
import numpy as np
import warnings
from .. import kern
def build_lcm(input_dim, num_outputs, CK = [], NC = [], W_columns=1,W=None,kappa=None):
#TODO build_icm or build_lcm
"""
Builds a kernel for a linear coregionalization model
:input_dim: Input dimensionality
:num_outputs: Number of outputs
:param CK: List of coregionalized kernels (i.e., this will be multiplied by a coregionalize kernel).
:param K: List of kernels that will be added up together with CK, but won't be multiplied by a coregionalize kernel
:param W_columns: number tuples of the corregionalization parameters 'coregion_W'
:type W_columns: integer
"""
for k in CK:
if k.input_dim <> input_dim:
k.input_dim = input_dim
warnings.warn("kernel's input dimension overwritten to fit input_dim parameter.")
for k in NC:
if k.input_dim <> input_dim + 1:
k.input_dim = input_dim + 1
warnings.warn("kernel's input dimension overwritten to fit input_dim parameter.")
kernel = CK[0].prod(kern.coregionalize(num_outputs,W_columns,W,kappa),tensor=True)
for k in CK[1:]:
k_coreg = kern.coregionalize(num_outputs,W_columns,W,kappa)
kernel += k.prod(k_coreg,tensor=True)
for k in NC:
kernel += k
return kernel

2
GPy/util/plot_latent.py
View file

@ -15,7 +15,7 @@ def most_significant_input_dimensions(model, which_indices):
try:
input_1, input_2 = np.argsort(model.input_sensitivity())[::-1][:2]
except:
raise ValueError, "cannot Atomatically determine which dimensions to plot, please pass 'which_indices'"
raise ValueError, "cannot automatically determine which dimensions to plot, please pass 'which_indices'"
else:
input_1, input_2 = which_indices
return input_1, input_2

32
GPy/util/symbolic.py Normal file
View file

@ -0,0 +1,32 @@
from sympy import Function, S, oo, I, cos, sin
class sinc_grad(Function):
nargs = 1
def fdiff(self, argindex=1):
return ((2-x*x)*sin(self.args[0]) - 2*x*cos(x))/(x*x*x)
@classmethod
def eval(cls, x):
if x is S.Zero:
return S.Zero
else:
return (x*cos(x) - sin(x))/(x*x)
class sinc(Function):
nargs = 1
def fdiff(self, argindex=1):
return sinc_grad(self.args[0])
@classmethod
def eval(cls, x):
if x is S.Zero:
return S.One
else:
return sin(x)/x
def _eval_is_real(self):
return self.args[0].is_real

11
GPy/util/univariate_Gaussian.py
View file

@ -32,4 +32,15 @@ def std_norm_cdf(x):
x = float(x)
return weave.inline(code,arg_names=['x'],support_code=support_code)
def inv_std_norm_cdf(x):
"""
Inverse cumulative standard Gaussian distribution
Based on Winitzki, S. (2008)
"""
z = 2*x -1
ln1z2 = np.log(1-z**2)
a = 8*(np.pi -3)/(3*np.pi*(4-np.pi))
b = 2/(np.pi * a) + ln1z2/2
inv_erf = np.sign(z) * np.sqrt( np.sqrt(b**2 - ln1z2/a) - b )
return np.sqrt(2) * inv_erf

13
GPy/util/visualize.py
View file

@ -502,11 +502,14 @@ def data_play(Y, visualizer, frame_rate=30):
This example loads in the CMU mocap database (http://mocap.cs.cmu.edu) subject number 35 motion number 01. It then plays it using the mocap_show visualize object.
data = GPy.util.datasets.cmu_mocap(subject='35', train_motions=['01'])
Y = data['Y']
Y[:, 0:3] = 0. # Make figure walk in place
visualize = GPy.util.visualize.skeleton_show(Y[0, :], data['skel'])
GPy.util.visualize.data_play(Y, visualize)
.. code-block:: python
data = GPy.util.datasets.cmu_mocap(subject='35', train_motions=['01'])
Y = data['Y']
Y[:, 0:3] = 0. # Make figure walk in place
visualize = GPy.util.visualize.skeleton_show(Y[0, :], data['skel'])
GPy.util.visualize.data_play(Y, visualize)
"""

27
GPy/util/warping_functions.py
View file

@ -53,9 +53,11 @@ class TanhWarpingFunction(WarpingFunction):
self.num_parameters = 3 * self.n_terms
def f(self,y,psi):
"""transform y with f using parameter vector psi
"""
transform y with f using parameter vector psi
psi = [[a,b,c]]
f = \sum_{terms} a * tanh(b*(y+c))
::math::`f = \\sum_{terms} a * tanh(b*(y+c))`
"""
#1. check that number of params is consistent
@ -77,8 +79,7 @@ class TanhWarpingFunction(WarpingFunction):
"""
calculate the numerical inverse of f
== input ==
iterations: number of N.R. iterations
:param iterations: number of N.R. iterations
"""
@ -165,9 +166,11 @@ class TanhWarpingFunction_d(WarpingFunction):
self.num_parameters = 3 * self.n_terms + 1
def f(self,y,psi):
"""transform y with f using parameter vector psi
"""
Transform y with f using parameter vector psi
psi = [[a,b,c]]
f = \sum_{terms} a * tanh(b*(y+c))
:math:`f = \\sum_{terms} a * tanh(b*(y+c))`
"""
#1. check that number of params is consistent
@ -189,8 +192,7 @@ class TanhWarpingFunction_d(WarpingFunction):
"""
calculate the numerical inverse of f
== input ==
iterations: number of N.R. iterations
:param max_iterations: maximum number of N.R. iterations
"""
@ -214,12 +216,13 @@ class TanhWarpingFunction_d(WarpingFunction):
def fgrad_y(self, y, psi, return_precalc = False):
"""
gradient of f w.r.t to y ([N x 1])
returns: Nx1 vector of derivatives, unless return_precalc is true,
then it also returns the precomputed stuff
:returns: Nx1 vector of derivatives, unless return_precalc is true, then it also returns the precomputed stuff
"""
mpsi = psi.copy()
mpsi = psi.coSpy()
d = psi[-1]
mpsi = mpsi[:self.num_parameters-1].reshape(self.n_terms, 3)
@ -242,7 +245,7 @@ class TanhWarpingFunction_d(WarpingFunction):
"""
gradient of f w.r.t to y and psi
returns: NxIx4 tensor of partial derivatives
:returns: NxIx4 tensor of partial derivatives
"""