apunkt/GPy
1
0
Fork 0
mirror of https://github.com/SheffieldML/GPy.git synced 2026-05-11 21:12:38 +02:00
Code Issues Projects Releases Packages Wiki Activity Actions

manual merging

Browse source
This commit is contained in:
James Hensman 2015-04-09 15:46:40 +01:00
commit ea787fd376
130 changed files with 982 additions and 787 deletions
Download patch file Download diff file Expand all files Collapse all files

28
GPy/util/__init__.py
View file

@ -2,18 +2,18 @@
# Licensed under the BSD 3-clause license (see LICENSE.txt)
import linalg
import misc
import squashers
import warping_functions
import datasets
import mocap
import decorators
import classification
import subarray_and_sorting
import caching
import diag
import initialization
import multioutput
import linalg_gpu
from . import linalg
from . import misc
from . import squashers
from . import warping_functions
from . import datasets
from . import mocap
from . import decorators
from . import classification
from . import subarray_and_sorting
from . import caching
from . import diag
from . import initialization
from . import multioutput
from . import linalg_gpu

2
GPy/util/block_matrices.py
View file

@ -71,6 +71,6 @@ if __name__=='__main__':
A = np.zeros((5,5))
B = get_blocks(A,[2,3])
B[0,0] += 7
print B
print(B)
assert np.all(unblock(B) == A)

5
GPy/util/caching.py
View file

@ -2,6 +2,7 @@
# Licensed under the BSD 3-clause license (see LICENSE.txt)
from ..core.parameterization.observable import Observable
import collections, weakref
from functools import reduce
class Cacher(object):
def __init__(self, operation, limit=5, ignore_args=(), force_kwargs=()):
@ -148,10 +149,10 @@ class Cacher(object):
return Cacher(self.operation, self.limit, self.ignore_args, self.force_kwargs)
def __getstate__(self, memo=None):
raise NotImplementedError, "Trying to pickle Cacher object with function {}, pickling functions not possible.".format(str(self.operation))
raise NotImplementedError("Trying to pickle Cacher object with function {}, pickling functions not possible.".format(str(self.operation)))
def __setstate__(self, memo=None):
raise NotImplementedError, "Trying to pickle Cacher object with function {}, pickling functions not possible.".format(str(self.operation))
raise NotImplementedError("Trying to pickle Cacher object with function {}, pickling functions not possible.".format(str(self.operation)))
@property
def __name__(self):

57
GPy/util/choleskies.py
View file

@ -2,23 +2,28 @@
# Licensed under the GNU GPL version 3.0
import numpy as np
from scipy import weave
import linalg
from . import linalg
from .config import config
try:
from scipy import weave
except ImportError:
config.set('weave', 'working', 'False')
def safe_root(N):
i = np.sqrt(N)
j = int(i)
if i != j:
raise ValueError, "N is not square!"
raise ValueError("N is not square!")
return j
def flat_to_triang(flat):
def _flat_to_triang_weave(flat):
"""take a matrix N x D and return a M X M x D array where
N = M(M+1)/2
the lower triangluar portion of the d'th slice of the result is filled by the d'th column of flat.
This is the weave implementation
"""
N, D = flat.shape
M = (-1 + safe_root(8*N+1))/2
@ -42,7 +47,24 @@ def flat_to_triang(flat):
weave.inline(code, ['flat', 'ret', 'D', 'M'])
return ret
def triang_to_flat(L):
def _flat_to_triang_pure(flat_mat):
N, D = flat_mat.shape
M = (-1 + safe_root(8*N+1))//2
ret = np.zeros((M, M, D))
count = 0
for m in range(M):
for mm in range(m+1):
for d in range(D):
ret.flat[d + m*D*M + mm*D] = flat_mat.flat[count];
count = count+1
return ret
if config.getboolean('weave', 'working'):
flat_to_triang = _flat_to_triang_weave
else:
flat_to_triang = _flat_to_triang_pure
def _triang_to_flat_weave(L):
M, _, D = L.shape
L = np.ascontiguousarray(L) # should do nothing if L was created by flat_to_triang
@ -66,13 +88,31 @@ def triang_to_flat(L):
weave.inline(code, ['flat', 'L', 'D', 'M'])
return flat
def _triang_to_flat_pure(L):
M, _, D = L.shape
N = M*(M+1)//2
flat = np.empty((N, D))
count = 0;
for m in range(M):
for mm in range(m+1):
for d in range(D):
flat.flat[count] = L.flat[d + m*D*M + mm*D];
count = count +1
return flat
if config.getboolean('weave', 'working'):
triang_to_flat = _triang_to_flat_weave
else:
triang_to_flat = _triang_to_flat_pure
def triang_to_cov(L):
return np.dstack([np.dot(L[:,:,i], L[:,:,i].T) for i in xrange(L.shape[-1])])
return np.dstack([np.dot(L[:,:,i], L[:,:,i].T) for i in range(L.shape[-1])])
def multiple_dpotri_old(Ls):
M, _, D = Ls.shape
Kis = np.rollaxis(Ls, -1).copy()
[dpotri(Kis[i,:,:], overwrite_c=1, lower=1) for i in xrange(D)]
[dpotri(Kis[i,:,:], overwrite_c=1, lower=1) for i in range(D)]
code = """
for(int d=0; d<D; d++)
{
@ -93,9 +133,6 @@ def multiple_dpotri_old(Ls):
def multiple_dpotri(Ls):
return np.dstack([linalg.dpotri(np.asfortranarray(Ls[:,:,i]), lower=1)[0] for i in range(Ls.shape[-1])])
def indexes_to_fix_for_low_rank(rank, size):
"""
work out which indexes of the flatteneed array should be fixed if we want the cholesky to represent a low rank matrix

10
GPy/util/classification.py
View file

@ -25,9 +25,9 @@ def conf_matrix(p,labels,names=['1','0'],threshold=.5,show=True):
true_0 = labels.size - true_1 - false_0 - false_1
error = (false_1 + false_0)/np.float(labels.size)
if show:
print 100. - error * 100,'% instances correctly classified'
print '%-10s| %-10s| %-10s| ' % ('',names[0],names[1])
print '----------|------------|------------|'
print '%-10s| %-10s| %-10s| ' % (names[0],true_1,false_0)
print '%-10s| %-10s| %-10s| ' % (names[1],false_1,true_0)
print(100. - error * 100,'% instances correctly classified')
print('%-10s| %-10s| %-10s| ' % ('',names[0],names[1]))
print('----------|------------|------------|')
print('%-10s| %-10s| %-10s| ' % (names[0],true_1,false_0))
print('%-10s| %-10s| %-10s| ' % (names[1],false_1,true_0))
return error,true_1, false_1, true_0, false_0

15
GPy/util/config.py
View file

@ -1,9 +1,18 @@
#
# This loads the configuration
#
import ConfigParser
import os
config = ConfigParser.ConfigParser()
try:
#Attempt Python 2 ConfigParser setup
import ConfigParser
config = ConfigParser.ConfigParser()
except ImportError:
#Attempt Python 3 ConfigParser setup
import configparser
config = configparser.ConfigParser()
# This is the default configuration file that always needs to be present.
default_file = os.path.abspath(os.path.join(os.path.dirname( __file__ ), '..', 'defaults.cfg'))
@ -20,4 +29,4 @@ user_file = os.path.join(home,'.gpy_user.cfg')
config.readfp(open(default_file))
config.read([local_file, user_file])
if not config:
raise ValueError, "No configuration file found at either " + user_file + " or " + local_file + " or " + default_file + "."
raise ValueError("No configuration file found at either " + user_file + " or " + local_file + " or " + default_file + ".")

92
GPy/util/datasets.py
View file

@ -1,17 +1,17 @@
from __future__ import print_function
import csv
import os
import copy
import numpy as np
import GPy
import scipy.io
import cPickle as pickle
import zipfile
import tarfile
import datetime
import json
import re
from config import *
import sys
from .config import *
ipython_available=True
try:
@ -19,8 +19,20 @@ try:
except ImportError:
ipython_available=False
try:
#In Python 2, cPickle is faster. It does not exist in Python 3 but the underlying code is always used
#if available
import cPickle as pickle
except ImportError:
import pickle
import sys, urllib2
#A Python2/3 import handler - urllib2 changed its name in Py3 and was also reorganised
try:
from urllib2 import urlopen
from urllib2 import URLError
except ImportError:
from urllib.request import urlopen
from urllib.error import URLError
def reporthook(a,b,c):
# ',' at the end of the line is important!
@ -75,7 +87,7 @@ def prompt_user(prompt):
elif choice in no:
return False
else:
print("Your response was a " + choice)
print(("Your response was a " + choice))
print("Please respond with 'yes', 'y' or 'no', 'n'")
#return prompt_user()
@ -99,7 +111,7 @@ 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
print(file)
dir_name = os.path.join(data_path, store_directory)
if save_name is None: save_name = os.path.join(dir_name, file)
@ -107,12 +119,12 @@ def download_url(url, store_directory, save_name=None, messages=True, suffix='')
if suffix is None: suffix=''
print "Downloading ", url, "->", save_name
print("Downloading ", url, "->", save_name)
if not os.path.exists(dir_name):
os.makedirs(dir_name)
try:
response = urllib2.urlopen(url+suffix)
except urllib2.URLError, e:
response = urlopen(url+suffix)
except URLError as e:
if not hasattr(e, "code"):
raise
response = e
@ -150,7 +162,7 @@ def download_url(url, store_directory, save_name=None, messages=True, suffix='')
sys.stdout.write(status)
sys.stdout.flush()
sys.stdout.write(" "*(len(status)) + "\r")
print status
print(status)
# if we wanted to get more sophisticated maybe we should check the response code here again even for successes.
#with open(save_name, 'wb') as f:
# f.write(response.read())
@ -159,32 +171,32 @@ def download_url(url, store_directory, save_name=None, messages=True, suffix='')
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)
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((dr['details']))
print('')
if dr['citation']:
print('Please cite:')
print(dr['citation'])
print((dr['citation']))
print('')
if dr['size']:
print('After downloading the data will take up ' + str(dr['size']) + ' bytes of space.')
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(('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((dr['license']))
print('')
return True
else:
if dr['license']:
print('You must also agree to the following license:')
print(dr['license'])
print((dr['license']))
print('')
return prompt_user('Do you wish to proceed with the download? [yes/no]')
@ -495,18 +507,18 @@ def google_trends(query_terms=['big data', 'machine learning', 'data science'],
file = 'data.csv'
file_name = os.path.join(dir_path,file)
if not os.path.exists(file_name) or refresh_data:
print "Accessing Google trends to acquire the data. Note that repeated accesses will result in a block due to a google terms of service violation. Failure at this point may be due to such blocks."
print("Accessing Google trends to acquire the data. Note that repeated accesses will result in a block due to a google terms of service violation. Failure at this point may be due to such blocks.")
# quote the query terms.
quoted_terms = []
for term in query_terms:
quoted_terms.append(urllib2.quote(term))
print "Query terms: ", ', '.join(query_terms)
print("Query terms: ", ', '.join(query_terms))
print "Fetching query:"
print("Fetching query:")
query = 'http://www.google.com/trends/fetchComponent?q=%s&cid=TIMESERIES_GRAPH_0&export=3' % ",".join(quoted_terms)
data = urllib2.urlopen(query).read()
print "Done."
data = urlopen(query).read()
print("Done.")
# In the notebook they did some data cleaning: remove Javascript header+footer, and translate new Date(....,..,..) into YYYY-MM-DD.
header = """// Data table response\ngoogle.visualization.Query.setResponse("""
data = data[len(header):-2]
@ -520,8 +532,8 @@ def google_trends(query_terms=['big data', 'machine learning', 'data science'],
df.to_csv(file_name)
else:
print "Reading cached data for google trends. To refresh the cache set 'refresh_data=True' when calling this function."
print "Query terms: ", ', '.join(query_terms)
print("Reading cached data for google trends. To refresh the cache set 'refresh_data=True' when calling this function.")
print("Query terms: ", ', '.join(query_terms))
df = pandas.read_csv(file_name, parse_dates=[0])
@ -679,11 +691,11 @@ def ripley_synth(data_set='ripley_prnn_data'):
def global_average_temperature(data_set='global_temperature', num_train=1000, refresh_data=False):
path = os.path.join(data_path, data_set)
if data_available(data_set) and not refresh_data:
print 'Using cached version of the data set, to use latest version set refresh_data to True'
print('Using cached version of the data set, to use latest version set refresh_data to True')
else:
download_data(data_set)
data = np.loadtxt(os.path.join(data_path, data_set, 'GLBTS.long.data'))
print 'Most recent data observation from month ', data[-1, 1], ' in year ', data[-1, 0]
print('Most recent data observation from month ', data[-1, 1], ' in year ', data[-1, 0])
allX = data[data[:, 3]!=-99.99, 2:3]
allY = data[data[:, 3]!=-99.99, 3:4]
X = allX[:num_train, 0:1]
@ -695,11 +707,11 @@ def global_average_temperature(data_set='global_temperature', num_train=1000, re
def mauna_loa(data_set='mauna_loa', num_train=545, refresh_data=False):
path = os.path.join(data_path, data_set)
if data_available(data_set) and not refresh_data:
print 'Using cached version of the data set, to use latest version set refresh_data to True'
print('Using cached version of the data set, to use latest version set refresh_data to True')
else:
download_data(data_set)
data = np.loadtxt(os.path.join(data_path, data_set, 'co2_mm_mlo.txt'))
print 'Most recent data observation from month ', data[-1, 1], ' in year ', data[-1, 0]
print('Most recent data observation from month ', data[-1, 1], ' in year ', data[-1, 0])
allX = data[data[:, 3]!=-99.99, 2:3]
allY = data[data[:, 3]!=-99.99, 3:4]
X = allX[:num_train, 0:1]
@ -784,7 +796,7 @@ def hapmap3(data_set='hapmap3'):
from sys import stdout
import bz2
except ImportError as i:
raise i, "Need pandas for hapmap dataset, make sure to install pandas (http://pandas.pydata.org/) before loading the hapmap dataset"
raise i("Need pandas for hapmap dataset, make sure to install pandas (http://pandas.pydata.org/) before loading the hapmap dataset")
dir_path = os.path.join(data_path,'hapmap3')
hapmap_file_name = 'hapmap3_r2_b36_fwd.consensus.qc.poly'
@ -802,10 +814,10 @@ def hapmap3(data_set='hapmap3'):
if not reduce(lambda a,b: a and b, map(os.path.exists, preprocessed_data_paths)):
if not overide_manual_authorize and not prompt_user("Preprocessing requires ~25GB "
"of memory and can take a (very) long time, continue? [Y/n]"):
print "Preprocessing required for further usage."
print("Preprocessing required for further usage.")
return
status = "Preprocessing data, please be patient..."
print status
print(status)
def write_status(message, progress, status):
stdout.write(" "*len(status)); stdout.write("\r"); stdout.flush()
status = r"[{perc: <{ll}}] {message: <13s}".format(message=message, ll=20,
@ -873,13 +885,13 @@ def hapmap3(data_set='hapmap3'):
inandf = DataFrame(index=metadf.index, data=inan, columns=mapnp[:,1])
inandf.to_pickle(preprocessed_data_paths[2])
status=write_status('done :)', 100, status)
print ''
print('')
else:
print "loading snps..."
print("loading snps...")
snpsdf = read_pickle(preprocessed_data_paths[0])
print "loading metainfo..."
print("loading metainfo...")
metadf = read_pickle(preprocessed_data_paths[1])
print "loading nan entries..."
print("loading nan entries...")
inandf = read_pickle(preprocessed_data_paths[2])
snps = snpsdf.values
populations = metadf.population.values.astype('S3')
@ -1001,7 +1013,7 @@ def singlecell_rna_seq_deng(dataset='singlecell_deng'):
# Extract the tar file
filename = os.path.join(dir_path, 'GSE45719_Raw.tar')
with tarfile.open(filename, 'r') as files:
print "Extracting Archive {}...".format(files.name)
print("Extracting Archive {}...".format(files.name))
data = None
gene_info = None
message = ''
@ -1010,9 +1022,9 @@ def singlecell_rna_seq_deng(dataset='singlecell_deng'):
for i, file_info in enumerate(members):
f = files.extractfile(file_info)
inner = read_csv(f, sep='\t', header=0, compression='gzip', index_col=0)
print ' '*(len(message)+1) + '\r',
print(' '*(len(message)+1) + '\r', end=' ')
message = "{: >7.2%}: Extracting: {}".format(float(i+1)/overall, file_info.name[:20]+"...txt.gz")
print message,
print(message, end=' ')
if data is None:
data = inner.RPKM.to_frame()
data.columns = [file_info.name[:-18]]
@ -1035,8 +1047,8 @@ def singlecell_rna_seq_deng(dataset='singlecell_deng'):
sys.stdout.write(' '*len(message) + '\r')
sys.stdout.flush()
print
print "Read Archive {}".format(files.name)
print()
print("Read Archive {}".format(files.name))
return data_details_return({'Y': data,
'series_info': info,

8
GPy/util/debug.py
View file

@ -13,7 +13,7 @@ def checkFinite(arr, name=None):
if np.any(np.logical_not(np.isfinite(arr))):
idx = np.where(np.logical_not(np.isfinite(arr)))[0]
print name+' at indices '+str(idx)+' have not finite values: '+str(arr[idx])+'!'
print(name+' at indices '+str(idx)+' have not finite values: '+str(arr[idx])+'!')
return False
return True
@ -23,13 +23,13 @@ def checkFullRank(m, tol=1e-10, name=None, force_check=False):
assert len(m.shape)==2 and m.shape[0]==m.shape[1], 'The input of checkFullRank has to be a square matrix!'
if not force_check and m.shape[0]>=10000:
print 'The size of '+name+'is too big to check (>=10000)!'
print('The size of '+name+'is too big to check (>=10000)!')
return True
s = np.real(np.linalg.eigvals(m))
if s.min()/s.max()<tol:
print name+' is close to singlar!'
print 'The eigen values of '+name+' is '+str(s)
print(name+' is close to singlar!')
print('The eigen values of '+name+' is '+str(s))
return False
return True

2
GPy/util/gpu_init.py
View file

@ -23,7 +23,7 @@ try:
import pycuda.driver
pycuda.driver.init()
if gpuid>=pycuda.driver.Device.count():
print '['+MPI.Get_processor_name()+'] more processes than the GPU numbers!'
print('['+MPI.Get_processor_name()+'] more processes than the GPU numbers!')
#MPI.COMM_WORLD.Abort()
raise
gpu_device = pycuda.driver.Device(gpuid)

93
GPy/util/linalg.py
View file

@ -6,16 +6,22 @@
# http://homepages.inf.ed.ac.uk/imurray2/code/tdot/tdot.py
import numpy as np
from scipy import linalg, weave
from scipy import linalg
import types
import ctypes
from ctypes import byref, c_char, c_int, c_double # TODO
import scipy
import warnings
import os
from config import config
from .config import config
import logging
try:
from scipy import weave
except ImportError:
config.set('weave', 'working', 'False')
_scipyversion = np.float64((scipy.__version__).split('.')[:2])
_fix_dpotri_scipy_bug = True
if np.all(_scipyversion >= np.array([0, 14])):
@ -34,7 +40,7 @@ if config.getboolean('anaconda', 'installed') and config.getboolean('anaconda',
dsyrk = mkl_rt.dsyrk
dsyr = mkl_rt.dsyr
_blas_available = True
print 'anaconda installed and mkl is loaded'
print('anaconda installed and mkl is loaded')
except:
_blas_available = False
else:
@ -64,7 +70,7 @@ def force_F_ordered(A):
"""
if A.flags['F_CONTIGUOUS']:
return A
print "why are your arrays not F order?"
print("why are your arrays not F order?")
return np.asfortranarray(A)
# def jitchol(A, maxtries=5):
@ -91,19 +97,19 @@ def jitchol(A, maxtries=5):
else:
diagA = np.diag(A)
if np.any(diagA <= 0.):
raise linalg.LinAlgError, "not pd: non-positive diagonal elements"
raise linalg.LinAlgError("not pd: non-positive diagonal elements")
jitter = diagA.mean() * 1e-6
num_tries = 1
while num_tries <= maxtries and np.isfinite(jitter):
try:
print jitter
print(jitter)
L = linalg.cholesky(A + np.eye(A.shape[0]) * jitter, lower=True)
return L
except:
jitter *= 10
finally:
num_tries += 1
raise linalg.LinAlgError, "not positive definite, even with jitter."
raise linalg.LinAlgError("not positive definite, even with jitter.")
import traceback
try: raise
except:
@ -213,12 +219,12 @@ def mdot(*args):
def _mdot_r(a, b):
"""Recursive helper for mdot"""
if type(a) == types.TupleType:
if type(a) == tuple:
if len(a) > 1:
a = mdot(*a)
else:
a = a[0]
if type(b) == types.TupleType:
if type(b) == tuple:
if len(b) > 1:
b = mdot(*b)
else:
@ -293,7 +299,7 @@ def pca(Y, input_dim):
"""
if not np.allclose(Y.mean(axis=0), 0.0):
print "Y is not zero mean, centering it locally (GPy.util.linalg.pca)"
print("Y is not zero mean, centering it locally (GPy.util.linalg.pca)")
# Y -= Y.mean(axis=0)
@ -352,16 +358,16 @@ def tdot_blas(mat, out=None):
# of C order. However, I tried that and had errors with large matrices:
# http://homepages.inf.ed.ac.uk/imurray2/code/tdot/tdot_broken.py
mat = np.asfortranarray(mat)
TRANS = c_char('n')
TRANS = c_char('n'.encode('ascii'))
N = c_int(mat.shape[0])
K = c_int(mat.shape[1])
LDA = c_int(mat.shape[0])
UPLO = c_char('l')
UPLO = c_char('l'.encode('ascii'))
ALPHA = c_double(1.0)
A = mat.ctypes.data_as(ctypes.c_void_p)
BETA = c_double(0.0)
C = out.ctypes.data_as(ctypes.c_void_p)
LDC = c_int(np.max(out.strides) / 8)
LDC = c_int(np.max(out.strides) // 8)
dsyrk(byref(UPLO), byref(TRANS), byref(N), byref(K),
byref(ALPHA), A, byref(LDA), byref(BETA), C, byref(LDC))
@ -388,7 +394,7 @@ def DSYR_blas(A, x, alpha=1.):
"""
N = c_int(A.shape[0])
LDA = c_int(A.shape[0])
UPLO = c_char('l')
UPLO = c_char('l'.encode('ascii'))
ALPHA = c_double(alpha)
A_ = A.ctypes.data_as(ctypes.c_void_p)
x_ = x.ctypes.data_as(ctypes.c_void_p)
@ -428,7 +434,7 @@ def symmetrify(A, upper=False):
try:
symmetrify_weave(A, upper)
except:
print "\n Weave compilation failed. Falling back to (slower) numpy implementation\n"
print("\n Weave compilation failed. Falling back to (slower) numpy implementation\n")
config.set('weave', 'working', 'False')
symmetrify_numpy(A, upper)
else:
@ -494,34 +500,35 @@ def symmetrify_numpy(A, upper=False):
else:
A[triu] = A.T[triu]
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
"""
support_code = """
#include <math.h>
"""
code = """
double r,c,s;
int j,i;
for(j=0; j<N; j++){
r = sqrt(L(j,j)*L(j,j) + x(j)*x(j));
c = r / L(j,j);
s = x(j) / L(j,j);
L(j,j) = r;
for (i=j+1; i<N; i++){
L(i,j) = (L(i,j) + s*x(i))/c;
x(i) = c*x(i) - s*L(i,j);
}
}
"""
x = x.copy()
N = x.size
weave.inline(code, support_code=support_code, arg_names=['N', 'L', 'x'], type_converters=weave.converters.blitz)
#This function appears to be unused. It's use of weave makes it problematic
#Commenting out for now
#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
# """
# support_code = """
# #include <math.h>
# """
# code = """
# double r,c,s;
# int j,i;
# for(j=0; j<N; j++){
# r = sqrt(L(j,j)*L(j,j) + x(j)*x(j));
# c = r / L(j,j);
# s = x(j) / L(j,j);
# L(j,j) = r;
# for (i=j+1; i<N; i++){
# L(i,j) = (L(i,j) + s*x(i))/c;
# x(i) = c*x(i) - s*L(i,j);
# }
# }
# """
# x = x.copy()
# N = x.size
# weave.inline(code, support_code=support_code, arg_names=['N', 'L', 'x'], type_converters=weave.converters.blitz)
def backsub_both_sides(L, X, transpose='left'):
""" Return L^-T * X * L^-1, assumuing X is symmetrical and L is lower cholesky"""

4
GPy/util/ln_diff_erfs.py
View file

@ -6,7 +6,7 @@ try:
from scipy.special import erfcx, erf
except ImportError:
from scipy.special import erf
from erfcx import erfcx
from .erfcx import erfcx
import numpy as np
@ -35,7 +35,7 @@ def ln_diff_erfs(x1, x2, return_sign=False):
elif x2.size==1:
v = np.zeros(x1.shape)
else:
raise ValueError, "This function does not broadcast unless provided with a scalar."
raise ValueError("This function does not broadcast unless provided with a scalar.")
if x1.size == 1:
x1 = np.tile(x1, x2.shape)

4
GPy/util/misc.py
View file

@ -2,7 +2,7 @@
# Licensed under the BSD 3-clause license (see LICENSE.txt)
import numpy as np
from config import *
from .config import *
_lim_val = np.finfo(np.float64).max
@ -24,7 +24,7 @@ def chain_1(df_dg, dg_dx):
if np.all(dg_dx==1.):
return df_dg
if len(df_dg) > 1 and len(df_dg.shape)>1 and df_dg.shape[-1] > 1:
import ipdb; ipdb.set_trace() # XXX BREAKPOINT
#import ipdb; ipdb.set_trace() # XXX BREAKPOINT
raise NotImplementedError('Not implemented for matricies yet')
return df_dg * dg_dx

3
GPy/util/mocap.py
View file

@ -2,7 +2,6 @@ import os
import numpy as np
import math
from GPy.util import datasets as dat
import urllib2
class vertex:
def __init__(self, name, id, parents=[], children=[], meta = {}):
@ -174,7 +173,7 @@ class skeleton(tree):
return connection
def to_xyz(self, channels):
raise NotImplementedError, "this needs to be implemented to use the skeleton class"
raise NotImplementedError("this needs to be implemented to use the skeleton class")
def finalize(self):

2
GPy/util/multioutput.py
View file

@ -51,7 +51,7 @@ def ICM(input_dim, num_outputs, kernel, W_rank=1,W=None,kappa=None,name='ICM'):
:param W_rank: number tuples of the corregionalization parameters 'W'
:type W_rank: integer
"""
if kernel.input_dim <> input_dim:
if kernel.input_dim != input_dim:
kernel.input_dim = input_dim
warnings.warn("kernel's input dimension overwritten to fit input_dim parameter.")

4
GPy/util/parallel.py
View file

@ -27,7 +27,7 @@ def divide_data(datanum, rank, size):
residue = (datanum)%size
datanum_list = np.empty((size),dtype=np.int32)
for i in xrange(size):
for i in range(size):
if i<residue:
datanum_list[i] = int(datanum/size)+1
else:
@ -38,4 +38,4 @@ def divide_data(datanum, rank, size):
else:
size = datanum/size
offset = size*rank+residue
return offset, offset+size, datanum_list
return offset, offset+size, datanum_list

5
GPy/util/pca.py
View file

@ -13,6 +13,7 @@ except:
from numpy.linalg.linalg import LinAlgError
from operator import setitem
import itertools
from functools import reduce
class PCA(object):
"""
@ -47,7 +48,7 @@ class PCA(object):
X_ = numpy.ma.masked_array(X, inan)
self.mu = X_.mean(0).base
self.sigma = X_.std(0).base
reduce(lambda y,x: setitem(x[0], x[1], x[2]), itertools.izip(X.T, inan.T, self.mu), None)
reduce(lambda y,x: setitem(x[0], x[1], x[2]), zip(X.T, inan.T, self.mu), None)
X = X - self.mu
X = X / numpy.where(self.sigma == 0, 1e-30, self.sigma)
return X
@ -138,4 +139,4 @@ class PCA(object):
pylab.tight_layout()
except:
pass
return plots
return plots

1
GPy/util/univariate_Gaussian.py
View file

@ -3,7 +3,6 @@
# Licensed under the BSD 3-clause license (see LICENSE.txt)
import numpy as np
from scipy import weave
from scipy.special import ndtr as std_norm_cdf
#define a standard normal pdf

2
GPy/util/warping_functions.py
View file

@ -220,7 +220,7 @@ class TanhWarpingFunction_d(WarpingFunction):
y -= update
it += 1
if it == max_iterations:
print "WARNING!!! Maximum number of iterations reached in f_inv "
print("WARNING!!! Maximum number of iterations reached in f_inv ")
return y