mirror of
https://github.com/SheffieldML/GPy.git
synced 2026-06-11 15:15:15 +02:00
Merge branch 'devel' of github.com:SheffieldML/GPy into devel
This commit is contained in:
Ricardo
commit
7dd6031a0b
35 changed files with 1352 additions and 445 deletions
|
|
@ -33,13 +33,13 @@ class GP(Model):
|
|||
|
||||
assert X.ndim == 2
|
||||
if isinstance(X, (ObsAr, VariationalPosterior)):
|
||||
self.X = X
|
||||
else: self.X = ObsAr(X)
|
||||
self.X = X.copy()
|
||||
else: self.X = ObsAr(X.copy())
|
||||
|
||||
self.num_data, self.input_dim = self.X.shape
|
||||
|
||||
assert Y.ndim == 2
|
||||
self.Y = ObsAr(Y)
|
||||
self.Y = ObsAr(Y.copy())
|
||||
assert Y.shape[0] == self.num_data
|
||||
_, self.output_dim = self.Y.shape
|
||||
|
||||
|
|
@ -180,40 +180,80 @@ class GP(Model):
|
|||
|
||||
return Ysim
|
||||
|
||||
def plot_f(self, *args, **kwargs):
|
||||
def plot_f(self, plot_limits=None, which_data_rows='all',
|
||||
which_data_ycols='all', fixed_inputs=[],
|
||||
levels=20, samples=0, fignum=None, ax=None, resolution=None,
|
||||
plot_raw=True,
|
||||
linecol=None,fillcol=None, Y_metadata=None, data_symbol='kx'):
|
||||
"""
|
||||
|
||||
Plot the GP's view of the world, where the data is normalized and
|
||||
before applying a likelihood.
|
||||
|
||||
This is a convenience function: arguments are passed to
|
||||
GPy.plotting.matplot_dep.models_plots.plot_f_fit
|
||||
|
||||
Plot the GP's view of the world, where the data is normalized and before applying a likelihood.
|
||||
This is a call to plot with plot_raw=True.
|
||||
Data will not be plotted in this, as the GP's view of the world
|
||||
may live in another space, or units then the data.
|
||||
"""
|
||||
assert "matplotlib" in sys.modules, "matplotlib package has not been imported."
|
||||
from ..plotting.matplot_dep import models_plots
|
||||
return models_plots.plot_fit_f(self,*args,**kwargs)
|
||||
kw = {}
|
||||
if linecol is not None:
|
||||
kw['linecol'] = linecol
|
||||
if fillcol is not None:
|
||||
kw['fillcol'] = fillcol
|
||||
return models_plots.plot_fit(self, plot_limits, which_data_rows,
|
||||
which_data_ycols, fixed_inputs,
|
||||
levels, samples, fignum, ax, resolution,
|
||||
plot_raw=plot_raw, Y_metadata=Y_metadata,
|
||||
data_symbol=data_symbol, **kw)
|
||||
|
||||
def plot(self, *args, **kwargs):
|
||||
def plot(self, plot_limits=None, which_data_rows='all',
|
||||
which_data_ycols='all', fixed_inputs=[],
|
||||
levels=20, samples=0, fignum=None, ax=None, resolution=None,
|
||||
plot_raw=False,
|
||||
linecol=None,fillcol=None, Y_metadata=None, data_symbol='kx'):
|
||||
"""
|
||||
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, use fixed_inputs to plot the GP with some of
|
||||
the inputs fixed.
|
||||
- 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, use fixed_inputs to plot the GP with some of the inputs fixed.
|
||||
|
||||
Can plot only part of the data and part of the posterior functions
|
||||
using which_data_rows which_data_ycols and which_parts
|
||||
|
||||
This is a convenience function: arguments are passed to
|
||||
GPy.plotting.matplot_dep.models_plots.plot_fit
|
||||
using which_data_rowsm which_data_ycols.
|
||||
|
||||
: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_rows: which of the training data to plot (default all)
|
||||
:type which_data_rows: 'all' or a slice object to slice model.X, model.Y
|
||||
:param which_data_ycols: when the data has several columns (independant outputs), only plot these
|
||||
:type which_data_rows: 'all' or a list of integers
|
||||
: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 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
|
||||
:type output: integer (first output is 0)
|
||||
:param linecol: color of line to plot [Tango.colorsHex['darkBlue']]
|
||||
:type linecol:
|
||||
:param fillcol: color of fill [Tango.colorsHex['lightBlue']]
|
||||
:param levels: for 2D plotting, the number of contour levels to use is ax is None, create a new figure
|
||||
"""
|
||||
assert "matplotlib" in sys.modules, "matplotlib package has not been imported."
|
||||
from ..plotting.matplot_dep import models_plots
|
||||
return models_plots.plot_fit(self,*args,**kwargs)
|
||||
kw = {}
|
||||
if linecol is not None:
|
||||
kw['linecol'] = linecol
|
||||
if fillcol is not None:
|
||||
kw['fillcol'] = fillcol
|
||||
return models_plots.plot_fit(self, plot_limits, which_data_rows,
|
||||
which_data_ycols, fixed_inputs,
|
||||
levels, samples, fignum, ax, resolution,
|
||||
plot_raw=plot_raw, Y_metadata=Y_metadata,
|
||||
data_symbol=data_symbol, **kw)
|
||||
|
||||
def input_sensitivity(self):
|
||||
"""
|
||||
|
|
@ -236,5 +276,9 @@ class GP(Model):
|
|||
TODO: valid args
|
||||
"""
|
||||
self.inference_method.on_optimization_start()
|
||||
super(GP, self).optimize(optimizer, start, **kwargs)
|
||||
self.inference_method.on_optimization_end()
|
||||
try:
|
||||
super(GP, self).optimize(optimizer, start, **kwargs)
|
||||
except KeyboardInterrupt:
|
||||
print "KeyboardInterrupt caught, calling on_optimization_end() to round things up"
|
||||
self.inference_method.on_optimization_end()
|
||||
raise
|
||||
|
|
@ -20,7 +20,7 @@ class Model(Parameterized):
|
|||
super(Model, self).__init__(name) # Parameterized.__init__(self)
|
||||
self.optimization_runs = []
|
||||
self.sampling_runs = []
|
||||
self.preferred_optimizer = 'scg'
|
||||
self.preferred_optimizer = 'bfgs'
|
||||
|
||||
def log_likelihood(self):
|
||||
raise NotImplementedError, "this needs to be implemented to use the model class"
|
||||
|
|
@ -61,7 +61,7 @@ class Model(Parameterized):
|
|||
on the current machine.
|
||||
|
||||
"""
|
||||
initial_parameters = self._get_params_transformed()
|
||||
initial_parameters = self.optimizer_array
|
||||
|
||||
if parallel:
|
||||
try:
|
||||
|
|
@ -124,13 +124,15 @@ class Model(Parameterized):
|
|||
|
||||
For probabilistic models this is the negative log_likelihood
|
||||
(including the MAP prior), so we return it here. If your model is not
|
||||
probabilistic, just return your objective here!
|
||||
probabilistic, just return your objective to minimize here!
|
||||
"""
|
||||
return -float(self.log_likelihood()) - self.log_prior()
|
||||
|
||||
def objective_function_gradients(self):
|
||||
"""
|
||||
The gradients for the objective function for the given algorithm.
|
||||
The gradients are w.r.t. the *negative* objective function, as
|
||||
this framework works with *negative* log-likelihoods as a default.
|
||||
|
||||
You can find the gradient for the parameters in self.gradient at all times.
|
||||
This is the place, where gradients get stored for parameters.
|
||||
|
|
@ -141,7 +143,7 @@ class Model(Parameterized):
|
|||
|
||||
For probabilistic models this is the gradient of the negative log_likelihood
|
||||
(including the MAP prior), so we return it here. If your model is not
|
||||
probabilistic, just return your gradient here!
|
||||
probabilistic, just return your *negative* gradient here!
|
||||
"""
|
||||
return -(self._log_likelihood_gradients() + self._log_prior_gradients())
|
||||
|
||||
|
|
@ -157,7 +159,8 @@ class Model(Parameterized):
|
|||
:type x: np.array
|
||||
"""
|
||||
try:
|
||||
self._set_params_transformed(x)
|
||||
# self._set_params_transformed(x)
|
||||
self.optimizer_array = x
|
||||
obj_grads = self._transform_gradients(self.objective_function_gradients())
|
||||
self._fail_count = 0
|
||||
except (LinAlgError, ZeroDivisionError, ValueError):
|
||||
|
|
@ -180,7 +183,7 @@ class Model(Parameterized):
|
|||
:parameter type: np.array
|
||||
"""
|
||||
try:
|
||||
self._set_params_transformed(x)
|
||||
self.optimizer_array = x
|
||||
obj = self.objective_function()
|
||||
self._fail_count = 0
|
||||
except (LinAlgError, ZeroDivisionError, ValueError):
|
||||
|
|
@ -192,7 +195,7 @@ class Model(Parameterized):
|
|||
|
||||
def _objective_grads(self, x):
|
||||
try:
|
||||
self._set_params_transformed(x)
|
||||
self.optimizer_array = x
|
||||
obj_f, obj_grads = self.objective_function(), self._transform_gradients(self.objective_function_gradients())
|
||||
self._fail_count = 0
|
||||
except (LinAlgError, ZeroDivisionError, ValueError):
|
||||
|
|
@ -226,7 +229,7 @@ class Model(Parameterized):
|
|||
optimizer = self.preferred_optimizer
|
||||
|
||||
if start == None:
|
||||
start = self._get_params_transformed()
|
||||
start = self.optimizer_array
|
||||
|
||||
optimizer = optimization.get_optimizer(optimizer)
|
||||
opt = optimizer(start, model=self, **kwargs)
|
||||
|
|
@ -235,7 +238,7 @@ class Model(Parameterized):
|
|||
|
||||
self.optimization_runs.append(opt)
|
||||
|
||||
self._set_params_transformed(opt.x_opt)
|
||||
self.optimizer_array = opt.x_opt
|
||||
|
||||
def optimize_SGD(self, momentum=0.1, learning_rate=0.01, iterations=20, **kwargs):
|
||||
# assert self.Y.shape[1] > 1, "SGD only works with D > 1"
|
||||
|
|
@ -260,7 +263,7 @@ class Model(Parameterized):
|
|||
The gradient is considered correct if the ratio of the analytical
|
||||
and numerical gradients is within <tolerance> of unity.
|
||||
"""
|
||||
x = self._get_params_transformed().copy()
|
||||
x = self.optimizer_array.copy()
|
||||
|
||||
if not verbose:
|
||||
# make sure only to test the selected parameters
|
||||
|
|
@ -270,8 +273,8 @@ class Model(Parameterized):
|
|||
transformed_index = self._raveled_index_for(target_param)
|
||||
if self._has_fixes():
|
||||
indices = np.r_[:self.size]
|
||||
which = (transformed_index[:,None]==indices[self._fixes_][None,:]).nonzero()
|
||||
transformed_index = (indices-(~self._fixes_).cumsum())[transformed_index[which[0]]]
|
||||
which = (transformed_index[:, None] == indices[self._fixes_][None, :]).nonzero()
|
||||
transformed_index = (indices - (~self._fixes_).cumsum())[transformed_index[which[0]]]
|
||||
|
||||
if transformed_index.size == 0:
|
||||
print "No free parameters to check"
|
||||
|
|
@ -290,7 +293,7 @@ class Model(Parameterized):
|
|||
gradient = gradient[transformed_index]
|
||||
|
||||
denominator = (2 * np.dot(dx, gradient))
|
||||
global_ratio = (f1 - f2) / np.where(denominator==0., 1e-32, denominator)
|
||||
global_ratio = (f1 - f2) / np.where(denominator == 0., 1e-32, denominator)
|
||||
global_diff = np.abs(f1 - f2) < tolerance and np.allclose(gradient, 0, atol=tolerance)
|
||||
if global_ratio is np.nan:
|
||||
global_ratio = 0
|
||||
|
|
@ -319,10 +322,10 @@ class Model(Parameterized):
|
|||
param_index = self._raveled_index_for(target_param)
|
||||
if self._has_fixes():
|
||||
indices = np.r_[:self.size]
|
||||
which = (param_index[:,None]==indices[self._fixes_][None,:]).nonzero()
|
||||
which = (param_index[:, None] == indices[self._fixes_][None, :]).nonzero()
|
||||
param_index = param_index[which[0]]
|
||||
transformed_index = (indices-(~self._fixes_).cumsum())[param_index]
|
||||
#print param_index, transformed_index
|
||||
transformed_index = (indices - (~self._fixes_).cumsum())[param_index]
|
||||
# print param_index, transformed_index
|
||||
else:
|
||||
transformed_index = param_index
|
||||
|
||||
|
|
@ -340,7 +343,7 @@ class Model(Parameterized):
|
|||
xx[xind] -= 2.*step
|
||||
f2 = self._objective(xx)
|
||||
numerical_gradient = (f1 - f2) / (2 * step)
|
||||
if np.all(gradient[xind]==0): ratio = (f1-f2) == gradient[xind]
|
||||
if np.all(gradient[xind] == 0): ratio = (f1 - f2) == gradient[xind]
|
||||
else: ratio = (f1 - f2) / (2 * step * gradient[xind])
|
||||
difference = np.abs((f1 - f2) / 2 / step - gradient[xind])
|
||||
|
||||
|
|
@ -358,7 +361,7 @@ class Model(Parameterized):
|
|||
grad_string = "{0:<{c0}}|{1:^{c1}}|{2:^{c2}}|{3:^{c3}}|{4:^{c4}}".format(formatted_name, r, d, g, ng, c0=cols[0] + 9, c1=cols[1], c2=cols[2], c3=cols[3], c4=cols[4])
|
||||
print grad_string
|
||||
|
||||
self._set_params_transformed(x)
|
||||
self.optimizer_array = x
|
||||
return ret
|
||||
|
||||
|
||||
|
|
|
|||
|
|
@ -7,6 +7,20 @@ import numpy
|
|||
from numpy.lib.function_base import vectorize
|
||||
from lists_and_dicts import IntArrayDict
|
||||
|
||||
def extract_properties_to_index(index, props):
|
||||
prop_index = dict()
|
||||
for i, cl in enumerate(props):
|
||||
for c in cl:
|
||||
ind = prop_index.get(c, list())
|
||||
ind.append(index[i])
|
||||
prop_index[c] = ind
|
||||
|
||||
for c, i in prop_index.items():
|
||||
prop_index[c] = numpy.array(i, dtype=int)
|
||||
|
||||
return prop_index
|
||||
|
||||
|
||||
class ParameterIndexOperations(object):
|
||||
'''
|
||||
Index operations for storing param index _properties
|
||||
|
|
@ -66,8 +80,34 @@ class ParameterIndexOperations(object):
|
|||
return self._properties.values()
|
||||
|
||||
def properties_for(self, index):
|
||||
"""
|
||||
Returns a list of properties, such that each entry in the list corresponds
|
||||
to the element of the index given.
|
||||
|
||||
Example:
|
||||
let properties: 'one':[1,2,3,4], 'two':[3,5,6]
|
||||
|
||||
>>> properties_for([2,3,5])
|
||||
[['one'], ['one', 'two'], ['two']]
|
||||
"""
|
||||
return vectorize(lambda i: [prop for prop in self.iterproperties() if i in self[prop]], otypes=[list])(index)
|
||||
|
||||
def properties_to_index_dict(self, index):
|
||||
"""
|
||||
Return a dictionary, containing properties as keys and indices as index
|
||||
Thus, the indices for each constraint, which is contained will be collected as
|
||||
one dictionary
|
||||
|
||||
Example:
|
||||
let properties: 'one':[1,2,3,4], 'two':[3,5,6]
|
||||
|
||||
>>> properties_to_index_dict([2,3,5])
|
||||
{'one':[2,3], 'two':[3,5]}
|
||||
"""
|
||||
props = self.properties_for(index)
|
||||
prop_index = extract_properties_to_index(index, props)
|
||||
return prop_index
|
||||
|
||||
def add(self, prop, indices):
|
||||
self._properties[prop] = combine_indices(self._properties[prop], indices)
|
||||
|
||||
|
|
@ -174,8 +214,32 @@ class ParameterIndexOperationsView(object):
|
|||
|
||||
|
||||
def properties_for(self, index):
|
||||
"""
|
||||
Returns a list of properties, such that each entry in the list corresponds
|
||||
to the element of the index given.
|
||||
|
||||
Example:
|
||||
let properties: 'one':[1,2,3,4], 'two':[3,5,6]
|
||||
|
||||
>>> properties_for([2,3,5])
|
||||
[['one'], ['one', 'two'], ['two']]
|
||||
"""
|
||||
return vectorize(lambda i: [prop for prop in self.iterproperties() if i in self[prop]], otypes=[list])(index)
|
||||
|
||||
def properties_to_index_dict(self, index):
|
||||
"""
|
||||
Return a dictionary, containing properties as keys and indices as index
|
||||
Thus, the indices for each constraint, which is contained will be collected as
|
||||
one dictionary
|
||||
|
||||
Example:
|
||||
let properties: 'one':[1,2,3,4], 'two':[3,5,6]
|
||||
|
||||
>>> properties_to_index_dict([2,3,5])
|
||||
{'one':[2,3], 'two':[3,5]}
|
||||
"""
|
||||
return extract_properties_to_index(index, self.properties_for(index))
|
||||
|
||||
|
||||
def add(self, prop, indices):
|
||||
self._param_index_ops.add(prop, indices+self._offset)
|
||||
|
|
|
|||
|
|
@ -77,8 +77,18 @@ class ObserverList(object):
|
|||
self._poc.insert(ins, (priority, weakref.ref(observer), callble))
|
||||
|
||||
def __str__(self):
|
||||
from . import ObsAr, Param
|
||||
from parameter_core import Parameterizable
|
||||
ret = []
|
||||
curr_p = None
|
||||
|
||||
def frmt(o):
|
||||
if isinstance(o, ObsAr):
|
||||
return 'ObsArr <{}>'.format(hex(id(o)))
|
||||
elif isinstance(o, (Param,Parameterizable)):
|
||||
return '{}'.format(o.hierarchy_name())
|
||||
else:
|
||||
return repr(o)
|
||||
for p, o, c in self:
|
||||
curr = ''
|
||||
if curr_p != p:
|
||||
|
|
@ -87,8 +97,9 @@ class ObserverList(object):
|
|||
else: curr_pre = " "*len(pre)
|
||||
curr_p = p
|
||||
curr += curr_pre
|
||||
ret.append(curr + ", ".join(map(repr, [o,c])))
|
||||
return '\n'.join(ret)
|
||||
|
||||
ret.append(curr + ", ".join([frmt(o), str(c)]))
|
||||
return '\n'.join(ret)
|
||||
|
||||
def flush(self):
|
||||
"""
|
||||
|
|
|
|||
|
|
@ -30,16 +30,22 @@ class ObsAr(np.ndarray, Pickleable, Observable):
|
|||
def __array_wrap__(self, out_arr, context=None):
|
||||
return out_arr.view(np.ndarray)
|
||||
|
||||
def _setup_observers(self):
|
||||
# do not setup anything, as observable arrays do not have default observers
|
||||
pass
|
||||
|
||||
def copy(self):
|
||||
from lists_and_dicts import ObserverList
|
||||
memo = {}
|
||||
memo[id(self)] = self
|
||||
memo[id(self.observers)] = ObserverList()
|
||||
return self.__deepcopy__(memo)
|
||||
|
||||
def __deepcopy__(self, memo):
|
||||
s = self.__new__(self.__class__, input_array=self.view(np.ndarray).copy())
|
||||
memo[id(self)] = s
|
||||
import copy
|
||||
s.__dict__.update(copy.deepcopy(self.__dict__, memo))
|
||||
Pickleable.__setstate__(s, copy.deepcopy(self.__getstate__(), memo))
|
||||
return s
|
||||
|
||||
def __reduce__(self):
|
||||
|
|
|
|||
|
|
@ -4,7 +4,7 @@
|
|||
import itertools
|
||||
import numpy
|
||||
np = numpy
|
||||
from parameter_core import Parameterizable, adjust_name_for_printing
|
||||
from parameter_core import Parameterizable, adjust_name_for_printing, Pickleable
|
||||
from observable_array import ObsAr
|
||||
|
||||
###### printing
|
||||
|
|
@ -173,36 +173,6 @@ class Param(Parameterizable, ObsAr):
|
|||
def _ensure_fixes(self):
|
||||
if not self._has_fixes(): self._fixes_ = numpy.ones(self._realsize_, dtype=bool)
|
||||
|
||||
#===========================================================================
|
||||
# parameterizable
|
||||
#===========================================================================
|
||||
def traverse(self, visit, *args, **kwargs):
|
||||
"""
|
||||
Traverse the hierarchy performing visit(self, *args, **kwargs) at every node passed by.
|
||||
See "visitor pattern" in literature. This is implemented in pre-order fashion.
|
||||
|
||||
This will function will just call visit on self, as Param are leaf nodes.
|
||||
"""
|
||||
self.__visited = True
|
||||
visit(self, *args, **kwargs)
|
||||
self.__visited = False
|
||||
|
||||
def traverse_parents(self, visit, *args, **kwargs):
|
||||
"""
|
||||
Traverse the hierarchy upwards, visiting all parents and their children, except self.
|
||||
See "visitor pattern" in literature. This is implemented in pre-order fashion.
|
||||
|
||||
Example:
|
||||
|
||||
parents = []
|
||||
self.traverse_parents(parents.append)
|
||||
print parents
|
||||
"""
|
||||
if self.has_parent():
|
||||
self.__visited = True
|
||||
self._parent_._traverse_parents(visit, *args, **kwargs)
|
||||
self.__visited = False
|
||||
|
||||
#===========================================================================
|
||||
# Convenience
|
||||
#===========================================================================
|
||||
|
|
@ -217,14 +187,24 @@ class Param(Parameterizable, ObsAr):
|
|||
#===========================================================================
|
||||
# Pickling and copying
|
||||
#===========================================================================
|
||||
def copy(self):
|
||||
return Parameterizable.copy(self, which=self)
|
||||
|
||||
def __deepcopy__(self, memo):
|
||||
s = self.__new__(self.__class__, name=self.name, input_array=self.view(numpy.ndarray).copy())
|
||||
memo[id(self)] = s
|
||||
memo[id(self)] = s
|
||||
import copy
|
||||
s.__dict__.update(copy.deepcopy(self.__dict__, memo))
|
||||
Pickleable.__setstate__(s, copy.deepcopy(self.__getstate__(), memo))
|
||||
return s
|
||||
|
||||
|
||||
def _setup_observers(self):
|
||||
"""
|
||||
Setup the default observers
|
||||
|
||||
1: pass through to parent, if present
|
||||
"""
|
||||
if self.has_parent():
|
||||
self.add_observer(self._parent_, self._parent_._pass_through_notify_observers, -np.inf)
|
||||
|
||||
#===========================================================================
|
||||
# Printing -> done
|
||||
#===========================================================================
|
||||
|
|
|
|||
|
|
@ -16,8 +16,9 @@ Observable Pattern for patameterization
|
|||
from transformations import Logexp, NegativeLogexp, Logistic, __fixed__, FIXED, UNFIXED
|
||||
import numpy as np
|
||||
import re
|
||||
import logging
|
||||
|
||||
__updated__ = '2014-05-15'
|
||||
__updated__ = '2014-05-21'
|
||||
|
||||
class HierarchyError(Exception):
|
||||
"""
|
||||
|
|
@ -49,12 +50,30 @@ class Observable(object):
|
|||
as an observer. Every time the observable changes, it sends a notification with
|
||||
self as only argument to all its observers.
|
||||
"""
|
||||
_updated = True
|
||||
_updates = True
|
||||
def __init__(self, *args, **kwargs):
|
||||
super(Observable, self).__init__()
|
||||
from lists_and_dicts import ObserverList
|
||||
self.observers = ObserverList()
|
||||
|
||||
@property
|
||||
def updates(self):
|
||||
p = getattr(self, '_highest_parent_', None)
|
||||
if p is not None:
|
||||
self._updates = p._updates
|
||||
return self._updates
|
||||
|
||||
@updates.setter
|
||||
def updates(self, ups):
|
||||
assert isinstance(ups, bool), "updates are either on (True) or off (False)"
|
||||
p = getattr(self, '_highest_parent_', None)
|
||||
if p is not None:
|
||||
p._updates = ups
|
||||
else:
|
||||
self._updates = ups
|
||||
if ups:
|
||||
self._trigger_params_changed()
|
||||
|
||||
def add_observer(self, observer, callble, priority=0):
|
||||
"""
|
||||
Add an observer `observer` with the callback `callble`
|
||||
|
|
@ -73,7 +92,7 @@ class Observable(object):
|
|||
for poc in self.observers:
|
||||
_, obs, clble = poc
|
||||
if callble is not None:
|
||||
if (obs == observer) and (callble == clble):
|
||||
if (obs is observer) and (callble == clble):
|
||||
to_remove.append(poc)
|
||||
else:
|
||||
if obs is observer:
|
||||
|
|
@ -91,6 +110,8 @@ class Observable(object):
|
|||
:param min_priority: only notify observers with priority > min_priority
|
||||
if min_priority is None, notify all observers in order
|
||||
"""
|
||||
if not self.updates:
|
||||
return
|
||||
if which is None:
|
||||
which = self
|
||||
if min_priority is None:
|
||||
|
|
@ -157,6 +178,7 @@ class Pickleable(object):
|
|||
"""
|
||||
def __init__(self, *a, **kw):
|
||||
super(Pickleable, self).__init__()
|
||||
|
||||
#===========================================================================
|
||||
# Pickling operations
|
||||
#===========================================================================
|
||||
|
|
@ -177,37 +199,46 @@ class Pickleable(object):
|
|||
#===========================================================================
|
||||
# copy and pickling
|
||||
#===========================================================================
|
||||
def copy(self):
|
||||
def copy(self, memo=None, which=None):
|
||||
"""
|
||||
Returns a (deep) copy of the current parameter handle.
|
||||
|
||||
All connections to parents of the copy will be cut.
|
||||
|
||||
:param dict memo: memo for deepcopy
|
||||
:param Parameterized which: parameterized object which started the copy process [default: self]
|
||||
"""
|
||||
#raise NotImplementedError, "Copy is not yet implemented, TODO: Observable hierarchy"
|
||||
if memo is None:
|
||||
memo = {}
|
||||
import copy
|
||||
memo = {}
|
||||
# the next part makes sure that we do not include parents in any form:
|
||||
parents = []
|
||||
self.traverse_parents(parents.append) # collect parents
|
||||
if which is None:
|
||||
which = self
|
||||
which.traverse_parents(parents.append) # collect parents
|
||||
for p in parents:
|
||||
memo[id(p)] = None # set all parents to be None, so they will not be copied
|
||||
memo[id(self.gradient)] = None # reset the gradient
|
||||
memo[id(self.param_array)] = None # and param_array
|
||||
memo[id(self._fixes_)] = None # fixes have to be reset, as this is now highest parent
|
||||
c = copy.deepcopy(self, memo) # and start the copy
|
||||
c._parent_index_ = None
|
||||
return c
|
||||
if not memo.has_key(id(p)):memo[id(p)] = None # set all parents to be None, so they will not be copied
|
||||
if not memo.has_key(id(self.gradient)):memo[id(self.gradient)] = None # reset the gradient
|
||||
if not memo.has_key(id(self._fixes_)):memo[id(self._fixes_)] = None # fixes have to be reset, as this is now highest parent
|
||||
copy = copy.deepcopy(self, memo) # and start the copy
|
||||
copy._parent_index_ = None
|
||||
copy._trigger_params_changed()
|
||||
return copy
|
||||
|
||||
def __deepcopy__(self, memo):
|
||||
s = self.__new__(self.__class__) # fresh instance
|
||||
memo[id(self)] = s # be sure to break all cycles --> self is already done
|
||||
import copy
|
||||
s.__dict__.update(copy.deepcopy(self.__dict__, memo)) # standard copy
|
||||
s.__setstate__(copy.deepcopy(self.__getstate__(), memo)) # standard copy
|
||||
return s
|
||||
|
||||
def __getstate__(self):
|
||||
ignore_list = ['_param_array_', # parameters get set from bottom to top
|
||||
'_gradient_array_', # as well as gradients
|
||||
'_optimizer_copy_',
|
||||
'logger',
|
||||
'observers',
|
||||
'_fixes_', # and fixes
|
||||
'_Cacher_wrap__cachers', # never pickle cachers
|
||||
]
|
||||
|
|
@ -219,7 +250,11 @@ class Pickleable(object):
|
|||
|
||||
def __setstate__(self, state):
|
||||
self.__dict__.update(state)
|
||||
return self
|
||||
from lists_and_dicts import ObserverList
|
||||
self.observers = ObserverList()
|
||||
self._setup_observers()
|
||||
self._optimizer_copy_transformed = False
|
||||
|
||||
|
||||
class Gradcheckable(Pickleable, Parentable):
|
||||
"""
|
||||
|
|
@ -371,8 +406,9 @@ class Indexable(Nameable, Observable):
|
|||
"""
|
||||
if value is not None:
|
||||
self[:] = value
|
||||
reconstrained = self.unconstrain()
|
||||
index = self._add_to_index_operations(self.constraints, reconstrained, __fixed__, warning)
|
||||
|
||||
index = self.unconstrain()
|
||||
index = self._add_to_index_operations(self.constraints, index, __fixed__, warning)
|
||||
self._highest_parent_._set_fixed(self, index)
|
||||
self.notify_observers(self, None if trigger_parent else -np.inf)
|
||||
return index
|
||||
|
|
@ -595,40 +631,104 @@ class OptimizationHandlable(Indexable):
|
|||
"""
|
||||
This enables optimization handles on an Object as done in GPy 0.4.
|
||||
|
||||
`..._transformed`: make sure the transformations and constraints etc are handled
|
||||
`..._optimizer_copy_transformed`: make sure the transformations and constraints etc are handled
|
||||
"""
|
||||
def __init__(self, name, default_constraint=None, *a, **kw):
|
||||
super(OptimizationHandlable, self).__init__(name, default_constraint=default_constraint, *a, **kw)
|
||||
self._optimizer_copy_ = None
|
||||
self._optimizer_copy_transformed = False
|
||||
|
||||
def _get_params_transformed(self):
|
||||
# transformed parameters (apply un-transformation rules)
|
||||
p = self.param_array.copy()
|
||||
[np.put(p, ind, c.finv(p[ind])) for c, ind in self.constraints.iteritems() if c != __fixed__]
|
||||
if self.has_parent() and self.constraints[__fixed__].size != 0:
|
||||
fixes = np.ones(self.size).astype(bool)
|
||||
fixes[self.constraints[__fixed__]] = FIXED
|
||||
return p[fixes]
|
||||
elif self._has_fixes():
|
||||
return p[self._fixes_]
|
||||
return p
|
||||
|
||||
def _set_params_transformed(self, p):
|
||||
#===========================================================================
|
||||
# Optimizer copy
|
||||
#===========================================================================
|
||||
@property
|
||||
def optimizer_array(self):
|
||||
"""
|
||||
Set parameters p, but make sure they get transformed before setting.
|
||||
This means, the optimizer sees p, whereas the model sees transformed(p),
|
||||
such that, the parameters the model sees are in the right domain.
|
||||
Array for the optimizer to work on.
|
||||
This array always lives in the space for the optimizer.
|
||||
Thus, it is untransformed, going from Transformations.
|
||||
|
||||
Setting this array, will make sure the transformed parameters for this model
|
||||
will be set accordingly. It has to be set with an array, retrieved from
|
||||
this method, as e.g. fixing will resize the array.
|
||||
|
||||
The optimizer should only interfere with this array, such that transofrmations
|
||||
are secured.
|
||||
"""
|
||||
if not(p is self.param_array):
|
||||
if self.__dict__.get('_optimizer_copy_', None) is None or self.size != self._optimizer_copy_.size:
|
||||
self._optimizer_copy_ = np.empty(self.size)
|
||||
|
||||
if not self._optimizer_copy_transformed:
|
||||
self._optimizer_copy_.flat = self.param_array.flat
|
||||
[np.put(self._optimizer_copy_, ind, c.finv(self.param_array[ind])) for c, ind in self.constraints.iteritems() if c != __fixed__]
|
||||
if self.has_parent() and self.constraints[__fixed__].size != 0:
|
||||
fixes = np.ones(self.size).astype(bool)
|
||||
fixes[self.constraints[__fixed__]] = FIXED
|
||||
self.param_array.flat[fixes] = p
|
||||
elif self._has_fixes(): self.param_array.flat[self._fixes_] = p
|
||||
else: self.param_array.flat = p
|
||||
[np.put(self.param_array, ind, c.f(self.param_array.flat[ind]))
|
||||
for c, ind in self.constraints.iteritems() if c != __fixed__]
|
||||
return self._optimizer_copy_[fixes]
|
||||
elif self._has_fixes():
|
||||
return self._optimizer_copy_[self._fixes_]
|
||||
self._optimizer_copy_transformed = True
|
||||
|
||||
return self._optimizer_copy_
|
||||
|
||||
@optimizer_array.setter
|
||||
def optimizer_array(self, p):
|
||||
"""
|
||||
Make sure the optimizer copy does not get touched, thus, we only want to
|
||||
set the values *inside* not the array itself.
|
||||
|
||||
Also we want to update param_array in here.
|
||||
"""
|
||||
f = None
|
||||
if self.has_parent() and self.constraints[__fixed__].size != 0:
|
||||
f = np.ones(self.size).astype(bool)
|
||||
f[self.constraints[__fixed__]] = FIXED
|
||||
elif self._has_fixes():
|
||||
f = self._fixes_
|
||||
if f is None:
|
||||
self.param_array.flat = p
|
||||
[np.put(self.param_array, ind, c.f(self.param_array.flat[ind]))
|
||||
for c, ind in self.constraints.iteritems() if c != __fixed__]
|
||||
else:
|
||||
self.param_array.flat[f] = p
|
||||
[np.put(self.param_array, ind[f[ind]], c.f(self.param_array.flat[ind[f[ind]]]))
|
||||
for c, ind in self.constraints.iteritems() if c != __fixed__]
|
||||
|
||||
self._optimizer_copy_transformed = False
|
||||
self._trigger_params_changed()
|
||||
|
||||
def _get_params_transformed(self):
|
||||
raise DeprecationWarning, "_get|set_params{_optimizer_copy_transformed} is deprecated, use self.optimizer array insetad!"
|
||||
# # transformed parameters (apply un-transformation rules)
|
||||
# p = self.param_array.copy()
|
||||
# [np.put(p, ind, c.finv(p[ind])) for c, ind in self.constraints.iteritems() if c != __fixed__]
|
||||
# if self.has_parent() and self.constraints[__fixed__].size != 0:
|
||||
# fixes = np.ones(self.size).astype(bool)
|
||||
# fixes[self.constraints[__fixed__]] = FIXED
|
||||
# return p[fixes]
|
||||
# elif self._has_fixes():
|
||||
# return p[self._fixes_]
|
||||
# return p
|
||||
#
|
||||
def _set_params_transformed(self, p):
|
||||
raise DeprecationWarning, "_get|set_params{_optimizer_copy_transformed} is deprecated, use self.optimizer array insetad!"
|
||||
|
||||
# """
|
||||
# Set parameters p, but make sure they get transformed before setting.
|
||||
# This means, the optimizer sees p, whereas the model sees transformed(p),
|
||||
# such that, the parameters the model sees are in the right domain.
|
||||
# """
|
||||
# if not(p is self.param_array):
|
||||
# if self.has_parent() and self.constraints[__fixed__].size != 0:
|
||||
# fixes = np.ones(self.size).astype(bool)
|
||||
# fixes[self.constraints[__fixed__]] = FIXED
|
||||
# self.param_array.flat[fixes] = p
|
||||
# elif self._has_fixes(): self.param_array.flat[self._fixes_] = p
|
||||
# else: self.param_array.flat = p
|
||||
# [np.put(self.param_array, ind, c.f(self.param_array.flat[ind]))
|
||||
# for c, ind in self.constraints.iteritems() if c != __fixed__]
|
||||
# self._trigger_params_changed()
|
||||
|
||||
def _trigger_params_changed(self, trigger_parent=True):
|
||||
"""
|
||||
First tell all children to update,
|
||||
|
|
@ -707,7 +807,7 @@ class OptimizationHandlable(Indexable):
|
|||
x = rand_gen(loc=loc, scale=scale, size=self._size_transformed(), *args, **kwargs)
|
||||
# now draw from prior where possible
|
||||
[np.put(x, ind, p.rvs(ind.size)) for p, ind in self.priors.iteritems() if not p is None]
|
||||
self._set_params_transformed(x) # makes sure all of the tied parameters get the same init (since there's only one prior object...)
|
||||
self.optimizer_array = x # makes sure all of the tied parameters get the same init (since there's only one prior object...)
|
||||
|
||||
#===========================================================================
|
||||
# For shared memory arrays. This does nothing in Param, but sets the memory
|
||||
|
|
@ -766,6 +866,7 @@ class Parameterizable(OptimizationHandlable):
|
|||
self.parameters = ArrayList()
|
||||
self._param_array_ = None
|
||||
self._added_names_ = set()
|
||||
self.logger = logging.getLogger(self.__class__.__name__)
|
||||
self.__visited = False # for traversing in reverse order we need to know if we were here already
|
||||
|
||||
@property
|
||||
|
|
@ -876,6 +977,11 @@ class Parameterizable(OptimizationHandlable):
|
|||
self._remove_parameter_name(None, old_name)
|
||||
self._add_parameter_name(param)
|
||||
|
||||
def __setstate__(self, state):
|
||||
super(Parameterizable, self).__setstate__(state)
|
||||
self.logger = logging.getLogger(self.__class__.__name__)
|
||||
return self
|
||||
|
||||
#===========================================================================
|
||||
# notification system
|
||||
#===========================================================================
|
||||
|
|
@ -883,7 +989,16 @@ class Parameterizable(OptimizationHandlable):
|
|||
self.parameters_changed()
|
||||
def _pass_through_notify_observers(self, me, which=None):
|
||||
self.notify_observers(which=which)
|
||||
|
||||
def _setup_observers(self):
|
||||
"""
|
||||
Setup the default observers
|
||||
|
||||
1: parameters_changed_notify
|
||||
2: pass through to parent, if present
|
||||
"""
|
||||
self.add_observer(self, self._parameters_changed_notification, -100)
|
||||
if self.has_parent():
|
||||
self.add_observer(self._parent_, self._parent_._pass_through_notify_observers, -np.inf)
|
||||
#===========================================================================
|
||||
# From being parentable, we have to define the parent_change notification
|
||||
#===========================================================================
|
||||
|
|
|
|||
|
|
@ -272,8 +272,11 @@ class Parameterized(Parameterizable):
|
|||
def __setattr__(self, name, val):
|
||||
# override the default behaviour, if setting a param, so broadcasting can by used
|
||||
if hasattr(self, "parameters"):
|
||||
pnames = self.parameter_names(False, adjust_for_printing=True, recursive=False)
|
||||
if name in pnames: self.parameters[pnames.index(name)][:] = val; return
|
||||
try:
|
||||
pnames = self.parameter_names(False, adjust_for_printing=True, recursive=False)
|
||||
if name in pnames: self.parameters[pnames.index(name)][:] = val; return
|
||||
except AttributeError:
|
||||
pass
|
||||
object.__setattr__(self, name, val);
|
||||
|
||||
#===========================================================================
|
||||
|
|
@ -281,17 +284,24 @@ class Parameterized(Parameterizable):
|
|||
#===========================================================================
|
||||
def __setstate__(self, state):
|
||||
super(Parameterized, self).__setstate__(state)
|
||||
self._connect_parameters()
|
||||
self._connect_fixes()
|
||||
self._notify_parent_change()
|
||||
try:
|
||||
self._connect_parameters()
|
||||
self._connect_fixes()
|
||||
self._notify_parent_change()
|
||||
self.parameters_changed()
|
||||
except Exception as e:
|
||||
print "WARNING: caught exception {!s}, trying to continue".format(e)
|
||||
|
||||
self.parameters_changed()
|
||||
def copy(self):
|
||||
c = super(Parameterized, self).copy()
|
||||
c._connect_parameters()
|
||||
c._connect_fixes()
|
||||
c._notify_parent_change()
|
||||
return c
|
||||
def copy(self, memo=None):
|
||||
if memo is None:
|
||||
memo = {}
|
||||
memo[id(self.optimizer_array)] = None # and param_array
|
||||
memo[id(self.param_array)] = None # and param_array
|
||||
copy = super(Parameterized, self).copy(memo)
|
||||
copy._connect_parameters()
|
||||
copy._connect_fixes()
|
||||
copy._notify_parent_change()
|
||||
return copy
|
||||
|
||||
#===========================================================================
|
||||
# Printing:
|
||||
|
|
|
|||
|
|
@ -66,7 +66,11 @@ class SparseGP(GP):
|
|||
#gradients wrt Z
|
||||
self.Z.gradient = self.kern.gradients_X(dL_dKmm, self.Z)
|
||||
self.Z.gradient += self.kern.gradients_Z_expectations(
|
||||
self.grad_dict['dL_dpsi0'], self.grad_dict['dL_dpsi1'], self.grad_dict['dL_dpsi2'], Z=self.Z, variational_posterior=self.X)
|
||||
self.grad_dict['dL_dpsi0'],
|
||||
self.grad_dict['dL_dpsi1'],
|
||||
self.grad_dict['dL_dpsi2'],
|
||||
Z=self.Z,
|
||||
variational_posterior=self.X)
|
||||
else:
|
||||
#gradients wrt kernel
|
||||
self.kern.update_gradients_diag(self.grad_dict['dL_dKdiag'], self.X)
|
||||
|
|
|
|||
|
|
@ -1,9 +1,14 @@
|
|||
# This is the configuration file for GPy
|
||||
# This is the default configuration file for GPy
|
||||
|
||||
# Do note edit this file.
|
||||
|
||||
# For machine specific changes (i.e. those specific to a given installation) edit GPy/installation.cfg
|
||||
|
||||
# For user specific changes edit $HOME/.gpy_user.cfg
|
||||
[parallel]
|
||||
# Enable openmp support. This speeds up some computations, depending on the number
|
||||
# of cores available. Setting up a compiler with openmp support can be difficult on
|
||||
# some platforms, hence this option.
|
||||
# some platforms, hence by default it is off.
|
||||
openmp=False
|
||||
|
||||
[datasets]
|
||||
|
|
@ -99,7 +99,7 @@ def sparse_gplvm_oil(optimize=True, verbose=0, plot=True, N=100, Q=6, num_induci
|
|||
m.kern.plot_ARD()
|
||||
return m
|
||||
|
||||
def swiss_roll(optimize=True, verbose=1, plot=True, N=1000, num_inducing=15, Q=4, sigma=.2):
|
||||
def swiss_roll(optimize=True, verbose=1, plot=True, N=1000, num_inducing=25, Q=4, sigma=.2):
|
||||
import GPy
|
||||
from GPy.util.datasets import swiss_roll_generated
|
||||
from GPy.models import BayesianGPLVM
|
||||
|
|
@ -144,16 +144,15 @@ def swiss_roll(optimize=True, verbose=1, plot=True, N=1000, num_inducing=15, Q=4
|
|||
m = BayesianGPLVM(Y, Q, X=X, X_variance=S, num_inducing=num_inducing, Z=Z, kernel=kernel)
|
||||
m.data_colors = c
|
||||
m.data_t = t
|
||||
m['noise_variance'] = Y.var() / 100.
|
||||
|
||||
|
||||
if optimize:
|
||||
m.optimize('scg', messages=verbose, max_iters=2e3)
|
||||
m.optimize('bfgs', messages=verbose, max_iters=2e3)
|
||||
|
||||
if plot:
|
||||
fig = plt.figure('fitted')
|
||||
ax = fig.add_subplot(111)
|
||||
s = m.input_sensitivity().argsort()[::-1][:2]
|
||||
ax.scatter(*m.X.T[s], c=c)
|
||||
ax.scatter(*m.X.mean.T[s], c=c)
|
||||
|
||||
return m
|
||||
|
||||
|
|
@ -172,14 +171,14 @@ def bgplvm_oil(optimize=True, verbose=1, plot=True, N=200, Q=7, num_inducing=40,
|
|||
m.data_labels = data['Y'][:N].argmax(axis=1)
|
||||
|
||||
if optimize:
|
||||
m.optimize('scg', messages=verbose, max_iters=max_iters, gtol=.05)
|
||||
m.optimize('bfgs', messages=verbose, max_iters=max_iters, gtol=.05)
|
||||
|
||||
if plot:
|
||||
fig, (latent_axes, sense_axes) = plt.subplots(1, 2)
|
||||
m.plot_latent(ax=latent_axes, labels=m.data_labels)
|
||||
data_show = GPy.plotting.matplot_dep.visualize.vector_show((m.Y[0,:]))
|
||||
lvm_visualizer = GPy.plotting.matplot_dep.visualize.lvm_dimselect(param_to_array(m.X.mean)[0:1,:], # @UnusedVariable
|
||||
m, data_show, latent_axes=latent_axes, sense_axes=sense_axes)
|
||||
m, data_show, latent_axes=latent_axes, sense_axes=sense_axes, labels=m.data_labels)
|
||||
raw_input('Press enter to finish')
|
||||
plt.close(fig)
|
||||
return m
|
||||
|
|
@ -303,9 +302,11 @@ def bgplvm_simulation_missing_data(optimize=True, verbose=1,
|
|||
k = kern.Linear(Q, ARD=True)# + kern.white(Q, _np.exp(-2)) # + kern.bias(Q)
|
||||
|
||||
inan = _np.random.binomial(1, .6, size=Y.shape).astype(bool)
|
||||
m = BayesianGPLVM(Y.copy(), Q, init="random", num_inducing=num_inducing, kernel=k)
|
||||
m.inference_method = VarDTCMissingData()
|
||||
m.Y[inan] = _np.nan
|
||||
Y[inan] = _np.nan
|
||||
|
||||
m = BayesianGPLVM(Y.copy(), Q, init="random", num_inducing=num_inducing,
|
||||
inference_method=VarDTCMissingData(inan=inan), kernel=k)
|
||||
|
||||
m.X.variance[:] = _np.random.uniform(0,.01,m.X.shape)
|
||||
m.likelihood.variance = .01
|
||||
m.parameters_changed()
|
||||
|
|
@ -338,7 +339,40 @@ def mrd_simulation(optimize=True, verbose=True, plot=True, plot_sim=True, **kw):
|
|||
print "Optimizing Model:"
|
||||
m.optimize(messages=verbose, max_iters=8e3, gtol=.1)
|
||||
if plot:
|
||||
m.plot_X_1d("MRD Latent Space 1D")
|
||||
m.X.plot("MRD Latent Space 1D")
|
||||
m.plot_scales("MRD Scales")
|
||||
return m
|
||||
|
||||
def mrd_simulation_missing_data(optimize=True, verbose=True, plot=True, plot_sim=True, **kw):
|
||||
from GPy import kern
|
||||
from GPy.models import MRD
|
||||
from GPy.inference.latent_function_inference.var_dtc import VarDTCMissingData
|
||||
|
||||
D1, D2, D3, N, num_inducing, Q = 60, 20, 36, 60, 6, 5
|
||||
_, _, Ylist = _simulate_sincos(D1, D2, D3, N, num_inducing, Q, plot_sim)
|
||||
|
||||
#Ylist = [Ylist[0]]
|
||||
k = kern.Linear(Q, ARD=True)
|
||||
inanlist = []
|
||||
|
||||
for Y in Ylist:
|
||||
inan = _np.random.binomial(1, .6, size=Y.shape).astype(bool)
|
||||
inanlist.append(inan)
|
||||
Y[inan] = _np.nan
|
||||
|
||||
imlist = []
|
||||
for inan in inanlist:
|
||||
imlist.append(VarDTCMissingData(limit=1, inan=inan))
|
||||
|
||||
m = MRD(Ylist, input_dim=Q, num_inducing=num_inducing,
|
||||
kernel=k, inference_method=imlist,
|
||||
initx="random", initz='permute', **kw)
|
||||
|
||||
if optimize:
|
||||
print "Optimizing Model:"
|
||||
m.optimize('bfgs', messages=verbose, max_iters=8e3, gtol=.1)
|
||||
if plot:
|
||||
m.X.plot("MRD Latent Space 1D")
|
||||
m.plot_scales("MRD Scales")
|
||||
return m
|
||||
|
||||
|
|
@ -351,18 +385,17 @@ def brendan_faces(optimize=True, verbose=True, plot=True):
|
|||
Yn = Y - Y.mean()
|
||||
Yn /= Yn.std()
|
||||
|
||||
m = GPy.models.GPLVM(Yn, Q)
|
||||
m = GPy.models.BayesianGPLVM(Yn, Q, num_inducing=20)
|
||||
|
||||
# optimize
|
||||
m.constrain('rbf|noise|white', GPy.transformations.LogexpClipped())
|
||||
|
||||
if optimize: m.optimize('scg', messages=verbose, max_iters=1000)
|
||||
if optimize: m.optimize('bfgs', messages=verbose, max_iters=1000)
|
||||
|
||||
if plot:
|
||||
ax = m.plot_latent(which_indices=(0, 1))
|
||||
y = m.likelihood.Y[0, :]
|
||||
y = m.Y[0, :]
|
||||
data_show = GPy.plotting.matplot_dep.visualize.image_show(y[None, :], dimensions=(20, 28), transpose=True, order='F', invert=False, scale=False)
|
||||
GPy.plotting.matplot_dep.visualize.lvm(m.X[0, :].copy(), m, data_show, ax)
|
||||
lvm = GPy.plotting.matplot_dep.visualize.lvm(m.X.mean[0, :].copy(), m, data_show, ax)
|
||||
raw_input('Press enter to finish')
|
||||
|
||||
return m
|
||||
|
|
@ -376,13 +409,14 @@ def olivetti_faces(optimize=True, verbose=True, plot=True):
|
|||
Yn = Y - Y.mean()
|
||||
Yn /= Yn.std()
|
||||
|
||||
m = GPy.models.GPLVM(Yn, Q)
|
||||
if optimize: m.optimize('scg', messages=verbose, max_iters=1000)
|
||||
m = GPy.models.BayesianGPLVM(Yn, Q, num_inducing=20)
|
||||
|
||||
if optimize: m.optimize('bfgs', messages=verbose, max_iters=1000)
|
||||
if plot:
|
||||
ax = m.plot_latent(which_indices=(0, 1))
|
||||
y = m.likelihood.Y[0, :]
|
||||
data_show = GPy.plotting.matplot_dep.visualize.image_show(y[None, :], dimensions=(112, 92), transpose=False, invert=False, scale=False)
|
||||
GPy.plotting.matplot_dep.visualize.lvm(m.X[0, :].copy(), m, data_show, ax)
|
||||
lvm = GPy.plotting.matplot_dep.visualize.lvm(m.X.mean[0, :].copy(), m, data_show, ax)
|
||||
raw_input('Press enter to finish')
|
||||
|
||||
return m
|
||||
|
|
@ -414,9 +448,10 @@ def stick(kernel=None, optimize=True, verbose=True, plot=True):
|
|||
ax = m.plot_latent()
|
||||
y = m.Y[0, :]
|
||||
data_show = GPy.plotting.matplot_dep.visualize.stick_show(y[None, :], connect=data['connect'])
|
||||
vis = GPy.plotting.matplot_dep.visualize.lvm(m.X[:1, :].copy(), m, data_show, latent_axes=ax)
|
||||
lvm_visualizer = GPy.plotting.matplot_dep.visualize.lvm(m.X[:1, :].copy(), m, data_show, latent_axes=ax)
|
||||
raw_input('Press enter to finish')
|
||||
|
||||
lvm_visualizer.close()
|
||||
data_show.close()
|
||||
return m
|
||||
|
||||
def bcgplvm_linear_stick(kernel=None, optimize=True, verbose=True, plot=True):
|
||||
|
|
@ -464,9 +499,8 @@ def robot_wireless(optimize=True, verbose=True, plot=True):
|
|||
|
||||
data = GPy.util.datasets.robot_wireless()
|
||||
# optimize
|
||||
m = GPy.models.GPLVM(data['Y'], 2)
|
||||
m = GPy.models.BayesianGPLVM(data['Y'], 4, num_inducing=25)
|
||||
if optimize: m.optimize(messages=verbose, max_f_eval=10000)
|
||||
m._set_params(m._get_params())
|
||||
if plot:
|
||||
m.plot_latent()
|
||||
|
||||
|
|
@ -482,21 +516,26 @@ def stick_bgplvm(model=None, optimize=True, verbose=True, plot=True):
|
|||
Q = 6
|
||||
kernel = GPy.kern.RBF(Q, lengthscale=np.repeat(.5, Q), ARD=True)
|
||||
m = BayesianGPLVM(data['Y'], Q, init="PCA", num_inducing=20, kernel=kernel)
|
||||
|
||||
|
||||
m.data = data
|
||||
m.likelihood.variance = 0.001
|
||||
|
||||
|
||||
# optimize
|
||||
if optimize: m.optimize('bfgs', messages=verbose, max_iters=800, xtol=1e-300, ftol=1e-300)
|
||||
try:
|
||||
if optimize: m.optimize('bfgs', messages=verbose, max_iters=5e3, bfgs_factor=10)
|
||||
except KeyboardInterrupt:
|
||||
print "Keyboard interrupt, continuing to plot and return"
|
||||
|
||||
if plot:
|
||||
plt.clf, (latent_axes, sense_axes) = plt.subplots(1, 2)
|
||||
fig, (latent_axes, sense_axes) = plt.subplots(1, 2)
|
||||
plt.sca(latent_axes)
|
||||
m.plot_latent(ax=latent_axes)
|
||||
y = m.Y[:1, :].copy()
|
||||
data_show = GPy.plotting.matplot_dep.visualize.stick_show(y, connect=data['connect'])
|
||||
GPy.plotting.matplot_dep.visualize.lvm_dimselect(m.X.mean[:1, :].copy(), m, data_show, latent_axes=latent_axes, sense_axes=sense_axes)
|
||||
plt.draw()
|
||||
#raw_input('Press enter to finish')
|
||||
dim_select = GPy.plotting.matplot_dep.visualize.lvm_dimselect(m.X.mean[:1, :].copy(), m, data_show, latent_axes=latent_axes, sense_axes=sense_axes)
|
||||
fig.canvas.draw()
|
||||
fig.canvas.show()
|
||||
raw_input('Press enter to finish')
|
||||
|
||||
return m
|
||||
|
||||
|
|
@ -515,9 +554,10 @@ def cmu_mocap(subject='35', motion=['01'], in_place=True, optimize=True, verbose
|
|||
ax = m.plot_latent()
|
||||
y = m.Y[0, :]
|
||||
data_show = GPy.plotting.matplot_dep.visualize.skeleton_show(y[None, :], data['skel'])
|
||||
lvm_visualizer = GPy.plotting.matplot_dep.visualize.lvm(m.X[0, :].copy(), m, data_show, ax)
|
||||
lvm_visualizer = GPy.plotting.matplot_dep.visualize.lvm(m.X[0].copy(), m, data_show, latent_axes=ax)
|
||||
raw_input('Press enter to finish')
|
||||
lvm_visualizer.close()
|
||||
data_show.close()
|
||||
|
||||
return m
|
||||
|
||||
|
|
|
|||
|
|
@ -38,6 +38,25 @@ class LatentFunctionInference(object):
|
|||
"""
|
||||
pass
|
||||
|
||||
class InferenceMethodList(LatentFunctionInference, list):
|
||||
|
||||
def on_optimization_start(self):
|
||||
for inf in self:
|
||||
inf.on_optimization_start()
|
||||
|
||||
def on_optimization_end(self):
|
||||
for inf in self:
|
||||
inf.on_optimization_end()
|
||||
|
||||
def __getstate__(self):
|
||||
state = []
|
||||
for inf in self:
|
||||
state.append(inf)
|
||||
return state
|
||||
|
||||
def __setstate__(self, state):
|
||||
for inf in state:
|
||||
self.append(inf)
|
||||
|
||||
from exact_gaussian_inference import ExactGaussianInference
|
||||
from laplace import Laplace
|
||||
|
|
|
|||
|
|
@ -95,7 +95,7 @@ class Posterior(object):
|
|||
"""
|
||||
if self._covariance is None:
|
||||
#LiK, _ = dtrtrs(self.woodbury_chol, self._K, lower=1)
|
||||
self._covariance = self._K - (np.tensordot(np.dot(np.atleast_3d(self.woodbury_inv).T, self._K), self._K, [1,0]).T).squeeze()
|
||||
self._covariance = (np.atleast_3d(self._K) - np.tensordot(np.dot(np.atleast_3d(self.woodbury_inv).T, self._K), self._K, [1,0]).T).squeeze()
|
||||
#self._covariance = self._K - self._K.dot(self.woodbury_inv).dot(self._K)
|
||||
return self._covariance
|
||||
|
||||
|
|
|
|||
|
|
@ -202,6 +202,17 @@ class VarDTCMissingData(LatentFunctionInference):
|
|||
def set_limit(self, limit):
|
||||
self._Y.limit = limit
|
||||
|
||||
def __getstate__(self):
|
||||
# has to be overridden, as Cacher objects cannot be pickled.
|
||||
return self._Y.limit, self._inan
|
||||
|
||||
def __setstate__(self, state):
|
||||
# has to be overridden, as Cacher objects cannot be pickled.
|
||||
from ...util.caching import Cacher
|
||||
self.limit = state[0]
|
||||
self._inan = state[1]
|
||||
self._Y = Cacher(self._subarray_computations, self.limit)
|
||||
|
||||
def _subarray_computations(self, Y):
|
||||
if self._inan is None:
|
||||
inan = np.isnan(Y)
|
||||
|
|
@ -272,7 +283,11 @@ class VarDTCMissingData(LatentFunctionInference):
|
|||
else: beta = beta_all
|
||||
|
||||
VVT_factor = (beta*y)
|
||||
VVT_factor_all[v, ind].flat = VVT_factor.flat
|
||||
try:
|
||||
VVT_factor_all[v, ind].flat = VVT_factor.flat
|
||||
except ValueError:
|
||||
mult = np.ravel_multi_index((v.nonzero()[0][:,None],ind[None,:]), VVT_factor_all.shape)
|
||||
VVT_factor_all.flat[mult] = VVT_factor
|
||||
output_dim = y.shape[1]
|
||||
|
||||
psi0 = psi0_all[v]
|
||||
|
|
|
|||
2
GPy/installation.cfg
Normal file
2
GPy/installation.cfg
Normal file
|
|
@ -0,0 +1,2 @@
|
|||
# This is the local configuration file for GPy
|
||||
|
||||
|
|
@ -14,6 +14,7 @@ from _src.ODE_UYC import ODE_UYC
|
|||
from _src.ODE_st import ODE_st
|
||||
from _src.ODE_t import ODE_t
|
||||
from _src.poly import Poly
|
||||
from _src.splitKern import SplitKern,DiffGenomeKern
|
||||
|
||||
# TODO: put this in an init file somewhere
|
||||
#I'm commenting this out because the files were not added. JH. Remember to add the files before commiting
|
||||
|
|
|
|||
|
|
@ -134,7 +134,7 @@ class Add(CombinationKernel):
|
|||
eff_dL_dpsi1 += dL_dpsi2.sum(1) * p2.variance * 2.
|
||||
else:
|
||||
eff_dL_dpsi1 += dL_dpsi2.sum(1) * p2.psi1(Z, variational_posterior) * 2.
|
||||
target += p1.gradients_Z_expectations(eff_dL_dpsi1, dL_dpsi2, Z, variational_posterior)
|
||||
target += p1.gradients_Z_expectations(dL_psi0, eff_dL_dpsi1, dL_dpsi2, Z, variational_posterior)
|
||||
return target
|
||||
|
||||
def gradients_qX_expectations(self, dL_dpsi0, dL_dpsi1, dL_dpsi2, Z, variational_posterior):
|
||||
|
|
|
|||
|
|
@ -22,7 +22,7 @@ def index_to_slices(index):
|
|||
"""
|
||||
|
||||
#contruct the return structure
|
||||
ind = np.asarray(index,dtype=np.int64)
|
||||
ind = np.asarray(index,dtype=np.int)
|
||||
ret = [[] for i in range(ind.max()+1)]
|
||||
|
||||
#find the switchpoints
|
||||
|
|
|
|||
188
GPy/kern/_src/splitKern.py
Normal file
188
GPy/kern/_src/splitKern.py
Normal file
|
|
@ -0,0 +1,188 @@
|
|||
"""
|
||||
A new kernel
|
||||
"""
|
||||
|
||||
import numpy as np
|
||||
from kern import Kern,CombinationKernel
|
||||
from .independent_outputs import index_to_slices
|
||||
import itertools
|
||||
|
||||
class DiffGenomeKern(Kern):
|
||||
|
||||
def __init__(self, kernel, idx_p, Xp, index_dim=-1, name='DiffGenomeKern'):
|
||||
self.idx_p = idx_p
|
||||
self.index_dim=index_dim
|
||||
self.kern = SplitKern(kernel,Xp, index_dim=index_dim)
|
||||
super(DiffGenomeKern, self).__init__(input_dim=kernel.input_dim+1, name=name)
|
||||
self.add_parameter(self.kern)
|
||||
|
||||
def K(self, X, X2=None):
|
||||
assert X2==None
|
||||
K = self.kern.K(X,X2)
|
||||
|
||||
slices = index_to_slices(X[:,self.index_dim])
|
||||
idx_start = slices[1][0]
|
||||
idx_end = idx_start+self.idx_p
|
||||
K[idx_start:idx_end,:] = K[:self.idx_p,:]
|
||||
K[:,idx_start:idx_end] = K[:,self.idx_p]
|
||||
|
||||
return K
|
||||
|
||||
def Kdiag(self,X):
|
||||
Kdiag = self.kern.Kdiag(X)
|
||||
|
||||
slices = index_to_slices(X[:,self.index_dim])
|
||||
idx_start = slices[1][0]
|
||||
idx_end = idx_start+self.idx_p
|
||||
Kdiag[idx_start:idx_end] = Kdiag[:self.idx_p]
|
||||
|
||||
return Kdiag
|
||||
|
||||
def update_gradients_full(self,dL_dK,X,X2=None):
|
||||
assert X2==None
|
||||
slices = index_to_slices(X[:,self.index_dim])
|
||||
idx_start = slices[1][0]
|
||||
idx_end = idx_start+self.idx_p
|
||||
|
||||
self.kern.update_gradients_full(dL_dK, X[:self.idx_p],X)
|
||||
grad_p1 = self.kern.gradient.copy()
|
||||
self.kern.update_gradients_full(dL_dK, X, X[:self.idx_p])
|
||||
grad_p2 = self.kern.gradient.copy()
|
||||
self.kern.update_gradients_full(dL_dK, X[:self.idx_p], X[:self.idx_p])
|
||||
grad_p3 = self.kern.gradient.copy()
|
||||
|
||||
self.kern.update_gradients_full(dL_dK, X[idx_start:idx_end],X)
|
||||
grad_n1 = self.kern.gradient.copy()
|
||||
self.kern.update_gradients_full(dL_dK, X, X[idx_start:idx_end])
|
||||
grad_n2 = self.kern.gradient.copy()
|
||||
self.kern.update_gradients_full(dL_dK, X[idx_start:idx_end], X[idx_start:idx_end])
|
||||
grad_n3 = self.kern.gradient.copy()
|
||||
|
||||
self.kern.update_gradients_full(dL_dK, X)
|
||||
self.kern.gradient += grad_p1+grad_p2+grad_p3-grad_n1-grad_n2-grad_n3
|
||||
|
||||
def update_gradients_diag(self, dL_dKdiag, X):
|
||||
pass
|
||||
|
||||
class SplitKern(CombinationKernel):
|
||||
|
||||
def __init__(self, kernel, Xp, index_dim=-1, name='SplitKern'):
|
||||
assert isinstance(index_dim, int), "The index dimension must be an integer!"
|
||||
self.kern = kernel
|
||||
self.kern_cross = SplitKern_cross(kernel,Xp)
|
||||
super(SplitKern, self).__init__(kernels=[self.kern, self.kern_cross], extra_dims=[index_dim], name=name)
|
||||
self.index_dim = index_dim
|
||||
|
||||
def K(self,X ,X2=None):
|
||||
slices = index_to_slices(X[:,self.index_dim])
|
||||
assert len(slices)<=2, 'The Split kernel only support two different indices'
|
||||
if X2 is None:
|
||||
target = np.zeros((X.shape[0], X.shape[0]))
|
||||
# diagonal blocks
|
||||
[[target.__setitem__((s,ss), self.kern.K(X[s,:], X[ss,:])) for s,ss in itertools.product(slices_i, slices_i)] for slices_i in slices]
|
||||
if len(slices)>1:
|
||||
# cross blocks
|
||||
[target.__setitem__((s,ss), self.kern_cross.K(X[s,:], X[ss,:])) for s,ss in itertools.product(slices[0], slices[1])]
|
||||
# cross blocks
|
||||
[target.__setitem__((s,ss), self.kern_cross.K(X[s,:], X[ss,:])) for s,ss in itertools.product(slices[1], slices[0])]
|
||||
else:
|
||||
slices2 = index_to_slices(X2[:,self.index_dim])
|
||||
assert len(slices2)<=2, 'The Split kernel only support two different indices'
|
||||
target = np.zeros((X.shape[0], X2.shape[0]))
|
||||
# diagonal blocks
|
||||
[[target.__setitem__((s,s2), self.kern.K(X[s,:],X2[s2,:])) for s,s2 in itertools.product(slices[i], slices2[i])] for i in xrange(min(len(slices),len(slices)))]
|
||||
if len(slices)>1:
|
||||
[target.__setitem__((s,s2), self.kern_cross.K(X[s,:],X2[s2,:])) for s,s2 in itertools.product(slices[1], slices2[0])]
|
||||
if len(slices2)>1:
|
||||
[target.__setitem__((s,s2), self.kern_cross.K(X[s,:],X2[s2,:])) for s,s2 in itertools.product(slices[0], slices2[1])]
|
||||
return target
|
||||
|
||||
def Kdiag(self,X):
|
||||
return self.kern.Kdiag(X)
|
||||
|
||||
def update_gradients_full(self,dL_dK,X,X2=None):
|
||||
slices = index_to_slices(X[:,self.index_dim])
|
||||
target = np.zeros(self.kern.size)
|
||||
|
||||
def collate_grads(dL, X, X2, cross=False):
|
||||
if cross:
|
||||
self.kern_cross.update_gradients_full(dL,X,X2)
|
||||
target[:] += self.kern_cross.kern.gradient
|
||||
else:
|
||||
self.kern.update_gradients_full(dL,X,X2)
|
||||
target[:] += self.kern.gradient
|
||||
|
||||
if X2 is None:
|
||||
[[collate_grads(dL_dK[s,ss], X[s], X[ss]) for s,ss in itertools.product(slices_i, slices_i)] for slices_i in slices]
|
||||
if len(slices)>1:
|
||||
[collate_grads(dL_dK[s,ss], X[s], X[ss], True) for s,ss in itertools.product(slices[0], slices[1])]
|
||||
[collate_grads(dL_dK[s,ss], X[s], X[ss], True) for s,ss in itertools.product(slices[1], slices[0])]
|
||||
else:
|
||||
slices2 = index_to_slices(X2[:,self.index_dim])
|
||||
[[collate_grads(dL_dK[s,s2],X[s],X2[s2]) for s,s2 in itertools.product(slices[i], slices2[i])] for i in xrange(min(len(slices),len(slices)))]
|
||||
if len(slices)>1:
|
||||
[collate_grads(dL_dK[s,ss], X[s], X2[s2], True) for s,s2 in itertools.product(slices[1], slices2[0])]
|
||||
if len(slices2)>1:
|
||||
[collate_grads(dL_dK[s,ss], X[s], X2[s2], True) for s,s2 in itertools.product(slices[0], slices2[1])]
|
||||
self.kern.gradient = target
|
||||
|
||||
def update_gradients_diag(self, dL_dKdiag, X):
|
||||
self.kern.update_gradients_diag(self, dL_dKdiag, X)
|
||||
|
||||
class SplitKern_cross(Kern):
|
||||
|
||||
def __init__(self, kernel, Xp, name='SplitKern_cross'):
|
||||
assert isinstance(kernel, Kern)
|
||||
self.kern = kernel
|
||||
if not isinstance(Xp,np.ndarray):
|
||||
Xp = np.array([[Xp]])
|
||||
self.Xp = Xp
|
||||
super(SplitKern_cross, self).__init__(input_dim=kernel.input_dim, active_dims=None, name=name)
|
||||
|
||||
def K(self, X, X2=None):
|
||||
if X2 is None:
|
||||
return np.dot(self.kern.K(X,self.Xp),self.kern.K(self.Xp,X))/self.kern.K(self.Xp,self.Xp)
|
||||
else:
|
||||
return np.dot(self.kern.K(X,self.Xp),self.kern.K(self.Xp,X2))/self.kern.K(self.Xp,self.Xp)
|
||||
|
||||
def Kdiag(self, X):
|
||||
return np.inner(self.kern.K(X,self.Xp),self.kern.K(self.Xp,X).T)/self.kern.K(self.Xp,self.Xp)
|
||||
|
||||
def update_gradients_full(self, dL_dK, X, X2=None):
|
||||
if X2 is None:
|
||||
X2 = X
|
||||
|
||||
k1 = self.kern.K(X,self.Xp)
|
||||
k2 = self.kern.K(self.Xp,X2)
|
||||
k3 = self.kern.K(self.Xp,self.Xp)
|
||||
dL_dk1 = np.einsum('ij,j->i',dL_dK,k2[0])/k3[0,0]
|
||||
dL_dk2 = np.einsum('ij,i->j',dL_dK,k1[:,0])/k3[0,0]
|
||||
dL_dk3 = np.einsum('ij,ij->',dL_dK,-np.dot(k1,k2)/(k3[0,0]*k3[0,0]))
|
||||
|
||||
self.kern.update_gradients_full(dL_dk1[:,None],X,self.Xp)
|
||||
grad = self.kern.gradient.copy()
|
||||
self.kern.update_gradients_full(dL_dk2[None,:],self.Xp,X2)
|
||||
grad += self.kern.gradient.copy()
|
||||
self.kern.update_gradients_full(np.array([[dL_dk3]]),self.Xp,self.Xp)
|
||||
grad += self.kern.gradient.copy()
|
||||
|
||||
self.kern.gradient = grad
|
||||
|
||||
def update_gradients_diag(self, dL_dKdiag, X):
|
||||
k1 = self.kern.K(X,self.Xp)
|
||||
k2 = self.kern.K(self.Xp,X)
|
||||
k3 = self.kern.K(self.Xp,self.Xp)
|
||||
dL_dk1 = dL_dKdiag*k2[0]/k3
|
||||
dL_dk2 = dL_dKdiag*k1[:,0]/k3
|
||||
dL_dk3 = -dL_dKdiag*(k1[:,0]*k2[0]).sum()/(k3*k3)
|
||||
|
||||
self.kern.update_gradients_full(dL_dk1[:,None],X,self.Xp)
|
||||
grad1 = self.kern.gradient.copy()
|
||||
self.kern.update_gradients_full(dL_dk2[None,:],self.Xp,X)
|
||||
grad2 = self.kern.gradient.copy()
|
||||
self.kern.update_gradients_full(np.array([[dL_dk3]]),self.Xp,self.Xp)
|
||||
grad3 = self.kern.gradient.copy()
|
||||
|
||||
self.kern.gradient = grad1+grad2+grad3
|
||||
|
||||
|
||||
|
|
@ -180,7 +180,7 @@ class Stationary(Kern):
|
|||
return np.zeros(X.shape)
|
||||
|
||||
def input_sensitivity(self):
|
||||
return np.ones(self.input_dim)/self.lengthscale
|
||||
return np.ones(self.input_dim)/self.lengthscale**2
|
||||
|
||||
class Exponential(Stationary):
|
||||
def __init__(self, input_dim, variance=1., lengthscale=None, ARD=False, active_dims=None, name='Exponential'):
|
||||
|
|
|
|||
61
GPy/mappings/additive.py
Normal file
61
GPy/mappings/additive.py
Normal file
|
|
@ -0,0 +1,61 @@
|
|||
# Copyright (c) 2013, GPy authors (see AUTHORS.txt).
|
||||
# Licensed under the BSD 3-clause license (see LICENSE.txt)
|
||||
|
||||
import numpy as np
|
||||
from ..core.mapping import Mapping
|
||||
import GPy
|
||||
|
||||
class Additive(Mapping):
|
||||
"""
|
||||
Mapping based on adding two existing mappings together.
|
||||
|
||||
.. math::
|
||||
|
||||
f(\mathbf{x}*) = f_1(\mathbf{x}*) + f_2(\mathbf(x)*)
|
||||
|
||||
:param mapping1: first mapping to add together.
|
||||
:type mapping1: GPy.mappings.Mapping
|
||||
:param mapping2: second mapping to add together.
|
||||
:type mapping2: GPy.mappings.Mapping
|
||||
:param tensor: whether or not to use the tensor product of input spaces
|
||||
:type tensor: bool
|
||||
|
||||
"""
|
||||
|
||||
def __init__(self, mapping1, mapping2, tensor=False):
|
||||
if tensor:
|
||||
input_dim = mapping1.input_dim + mapping2.input_dim
|
||||
else:
|
||||
input_dim = mapping1.input_dim
|
||||
assert(mapping1.input_dim==mapping2.input_dim)
|
||||
assert(mapping1.output_dim==mapping2.output_dim)
|
||||
output_dim = mapping1.output_dim
|
||||
Mapping.__init__(self, input_dim=input_dim, output_dim=output_dim)
|
||||
self.mapping1 = mapping1
|
||||
self.mapping2 = mapping2
|
||||
self.num_params = self.mapping1.num_params + self.mapping2.num_params
|
||||
self.name = self.mapping1.name + '+' + self.mapping2.name
|
||||
def _get_param_names(self):
|
||||
return self.mapping1._get_param_names + self.mapping2._get_param_names
|
||||
|
||||
def _get_params(self):
|
||||
return np.hstack((self.mapping1._get_params() self.mapping2._get_params()))
|
||||
|
||||
def _set_params(self, x):
|
||||
self.mapping1._set_params(x[:self.mapping1.num_params])
|
||||
self.mapping2._set_params(x[self.mapping1.num_params:])
|
||||
|
||||
def randomize(self):
|
||||
self.mapping1._randomize()
|
||||
self.mapping2._randomize()
|
||||
|
||||
def f(self, X):
|
||||
return self.mapping1.f(X) + self.mapping2.f(X)
|
||||
|
||||
def df_dtheta(self, dL_df, X):
|
||||
self._df_dA = (dL_df[:, :, None]*self.kern.K(X, self.X)[:, None, :]).sum(0).T
|
||||
self._df_dbias = (dL_df.sum(0))
|
||||
return np.hstack((self._df_dA.flatten(), self._df_dbias))
|
||||
|
||||
def df_dX(self, dL_df, X):
|
||||
return self.kern.dK_dX((dL_df[:, None, :]*self.A[None, :, :]).sum(2), X, self.X)
|
||||
|
|
@ -10,6 +10,7 @@ from ..util import linalg
|
|||
from ..core.parameterization.variational import NormalPosterior, NormalPrior, VariationalPosterior
|
||||
from ..inference.latent_function_inference.var_dtc_parallel import update_gradients
|
||||
from ..inference.latent_function_inference.var_dtc_gpu import VarDTC_GPU
|
||||
import logging
|
||||
|
||||
class BayesianGPLVM(SparseGP):
|
||||
"""
|
||||
|
|
@ -25,8 +26,10 @@ class BayesianGPLVM(SparseGP):
|
|||
"""
|
||||
def __init__(self, Y, input_dim, X=None, X_variance=None, init='PCA', num_inducing=10,
|
||||
Z=None, kernel=None, inference_method=None, likelihood=None, name='bayesian gplvm', **kwargs):
|
||||
self.logger = logging.getLogger(self.__class__.__name__)
|
||||
if X == None:
|
||||
from ..util.initialization import initialize_latent
|
||||
self.logger.info("initializing latent space X with method {}".format(init))
|
||||
X, fracs = initialize_latent(init, input_dim, Y)
|
||||
else:
|
||||
fracs = np.ones(input_dim)
|
||||
|
|
@ -36,7 +39,6 @@ class BayesianGPLVM(SparseGP):
|
|||
if X_variance is None:
|
||||
X_variance = np.random.uniform(0,.1,X.shape)
|
||||
|
||||
|
||||
if Z is None:
|
||||
Z = np.random.permutation(X.copy())[:num_inducing]
|
||||
assert Z.shape[1] == X.shape[1]
|
||||
|
|
@ -52,11 +54,14 @@ class BayesianGPLVM(SparseGP):
|
|||
X = NormalPosterior(X, X_variance)
|
||||
|
||||
if inference_method is None:
|
||||
if np.any(np.isnan(Y)):
|
||||
inan = np.isnan(Y)
|
||||
if np.any(inan):
|
||||
from ..inference.latent_function_inference.var_dtc import VarDTCMissingData
|
||||
inference_method = VarDTCMissingData()
|
||||
self.logger.debug("creating inference_method with var_dtc missing data")
|
||||
inference_method = VarDTCMissingData(inan=inan)
|
||||
else:
|
||||
from ..inference.latent_function_inference.var_dtc import VarDTC
|
||||
self.logger.debug("creating inference_method var_dtc")
|
||||
inference_method = VarDTC()
|
||||
|
||||
SparseGP.__init__(self, X, Y, Z, kernel, likelihood, inference_method, name, **kwargs)
|
||||
|
|
@ -83,7 +88,7 @@ class BayesianGPLVM(SparseGP):
|
|||
resolution=50, ax=None, marker='o', s=40,
|
||||
fignum=None, plot_inducing=True, legend=True,
|
||||
plot_limits=None,
|
||||
aspect='auto', updates=False, **kwargs):
|
||||
aspect='auto', updates=False, predict_kwargs={}, imshow_kwargs={}):
|
||||
import sys
|
||||
assert "matplotlib" in sys.modules, "matplotlib package has not been imported."
|
||||
from ..plotting.matplot_dep import dim_reduction_plots
|
||||
|
|
@ -91,7 +96,7 @@ class BayesianGPLVM(SparseGP):
|
|||
return dim_reduction_plots.plot_latent(self, labels, which_indices,
|
||||
resolution, ax, marker, s,
|
||||
fignum, plot_inducing, legend,
|
||||
plot_limits, aspect, updates, **kwargs)
|
||||
plot_limits, aspect, updates, predict_kwargs, imshow_kwargs)
|
||||
|
||||
def do_test_latents(self, Y):
|
||||
"""
|
||||
|
|
@ -100,36 +105,41 @@ class BayesianGPLVM(SparseGP):
|
|||
Notes:
|
||||
This will only work with a univariate Gaussian likelihood (for now)
|
||||
"""
|
||||
assert not self.likelihood.is_heteroscedastic
|
||||
N_test = Y.shape[0]
|
||||
input_dim = self.Z.shape[1]
|
||||
|
||||
means = np.zeros((N_test, input_dim))
|
||||
covars = np.zeros((N_test, input_dim))
|
||||
|
||||
dpsi0 = -0.5 * self.input_dim * self.likelihood.precision
|
||||
dpsi2 = self.dL_dpsi2[0][None, :, :] # TODO: this may change if we ignore het. likelihoods
|
||||
V = self.likelihood.precision * Y
|
||||
|
||||
dpsi0 = -0.5 * self.input_dim / self.likelihood.variance
|
||||
dpsi2 = self.grad_dict['dL_dpsi2'][0][None, :, :] # TODO: this may change if we ignore het. likelihoods
|
||||
V = Y/self.likelihood.variance
|
||||
|
||||
#compute CPsi1V
|
||||
if self.Cpsi1V is None:
|
||||
psi1V = np.dot(self.psi1.T, self.likelihood.V)
|
||||
tmp, _ = linalg.dtrtrs(self._Lm, np.asfortranarray(psi1V), lower=1, trans=0)
|
||||
tmp, _ = linalg.dpotrs(self.LB, tmp, lower=1)
|
||||
self.Cpsi1V, _ = linalg.dtrtrs(self._Lm, tmp, lower=1, trans=1)
|
||||
#if self.Cpsi1V is None:
|
||||
# psi1V = np.dot(self.psi1.T, self.likelihood.V)
|
||||
# tmp, _ = linalg.dtrtrs(self._Lm, np.asfortranarray(psi1V), lower=1, trans=0)
|
||||
# tmp, _ = linalg.dpotrs(self.LB, tmp, lower=1)
|
||||
# self.Cpsi1V, _ = linalg.dtrtrs(self._Lm, tmp, lower=1, trans=1)
|
||||
|
||||
dpsi1 = np.dot(self.Cpsi1V, V.T)
|
||||
dpsi1 = np.dot(self.posterior.woodbury_vector, V.T)
|
||||
|
||||
start = np.zeros(self.input_dim * 2)
|
||||
#start = np.zeros(self.input_dim * 2)
|
||||
|
||||
|
||||
from scipy.optimize import minimize
|
||||
|
||||
for n, dpsi1_n in enumerate(dpsi1.T[:, :, None]):
|
||||
args = (self.kern, self.Z, dpsi0, dpsi1_n.T, dpsi2)
|
||||
xopt, fopt, neval, status = SCG(f=latent_cost, gradf=latent_grad, x=start, optargs=args, display=False)
|
||||
|
||||
args = (input_dim, self.kern.copy(), self.Z, dpsi0, dpsi1_n.T, dpsi2)
|
||||
res = minimize(latent_cost_and_grad, jac=True, x0=np.hstack((means[n], covars[n])), args=args, method='BFGS')
|
||||
xopt = res.x
|
||||
mu, log_S = xopt.reshape(2, 1, -1)
|
||||
means[n] = mu[0].copy()
|
||||
covars[n] = np.exp(log_S[0]).copy()
|
||||
|
||||
return means, covars
|
||||
X = NormalPosterior(means, covars)
|
||||
|
||||
return X
|
||||
|
||||
def dmu_dX(self, Xnew):
|
||||
"""
|
||||
|
|
@ -161,57 +171,26 @@ class BayesianGPLVM(SparseGP):
|
|||
return dim_reduction_plots.plot_steepest_gradient_map(self,*args,**kwargs)
|
||||
|
||||
|
||||
def latent_cost_and_grad(mu_S, kern, Z, dL_dpsi0, dL_dpsi1, dL_dpsi2):
|
||||
def latent_cost_and_grad(mu_S, input_dim, kern, Z, dL_dpsi0, dL_dpsi1, dL_dpsi2):
|
||||
"""
|
||||
objective function for fitting the latent variables for test points
|
||||
(negative log-likelihood: should be minimised!)
|
||||
"""
|
||||
mu, log_S = mu_S.reshape(2, 1, -1)
|
||||
mu = mu_S[:input_dim][None]
|
||||
log_S = mu_S[input_dim:][None]
|
||||
S = np.exp(log_S)
|
||||
|
||||
psi0 = kern.psi0(Z, mu, S)
|
||||
psi1 = kern.psi1(Z, mu, S)
|
||||
psi2 = kern.psi2(Z, mu, S)
|
||||
X = NormalPosterior(mu, S)
|
||||
|
||||
lik = dL_dpsi0 * psi0 + np.dot(dL_dpsi1.flatten(), psi1.flatten()) + np.dot(dL_dpsi2.flatten(), psi2.flatten()) - 0.5 * np.sum(np.square(mu) + S) + 0.5 * np.sum(log_S)
|
||||
psi0 = kern.psi0(Z, X)
|
||||
psi1 = kern.psi1(Z, X)
|
||||
psi2 = kern.psi2(Z, X)
|
||||
|
||||
mu0, S0 = kern.dpsi0_dmuS(dL_dpsi0, Z, mu, S)
|
||||
mu1, S1 = kern.dpsi1_dmuS(dL_dpsi1, Z, mu, S)
|
||||
mu2, S2 = kern.dpsi2_dmuS(dL_dpsi2, Z, mu, S)
|
||||
|
||||
dmu = mu0 + mu1 + mu2 - mu
|
||||
lik = dL_dpsi0 * psi0.sum() + np.einsum('ij,kj->...', dL_dpsi1, psi1) + np.einsum('ijk,lkj->...', dL_dpsi2, psi2) - 0.5 * np.sum(np.square(mu) + S) + 0.5 * np.sum(log_S)
|
||||
|
||||
dLdmu, dLdS = kern.gradients_qX_expectations(dL_dpsi0, dL_dpsi1, dL_dpsi2, Z, X)
|
||||
dmu = dLdmu - mu
|
||||
# dS = S0 + S1 + S2 -0.5 + .5/S
|
||||
dlnS = S * (S0 + S1 + S2 - 0.5) + .5
|
||||
dlnS = S * (dLdS - 0.5) + .5
|
||||
|
||||
return -lik, -np.hstack((dmu.flatten(), dlnS.flatten()))
|
||||
|
||||
def latent_cost(mu_S, kern, Z, dL_dpsi0, dL_dpsi1, dL_dpsi2):
|
||||
"""
|
||||
objective function for fitting the latent variables (negative log-likelihood: should be minimised!)
|
||||
This is the same as latent_cost_and_grad but only for the objective
|
||||
"""
|
||||
mu, log_S = mu_S.reshape(2, 1, -1)
|
||||
S = np.exp(log_S)
|
||||
|
||||
psi0 = kern.psi0(Z, mu, S)
|
||||
psi1 = kern.psi1(Z, mu, S)
|
||||
psi2 = kern.psi2(Z, mu, S)
|
||||
|
||||
lik = dL_dpsi0 * psi0 + np.dot(dL_dpsi1.flatten(), psi1.flatten()) + np.dot(dL_dpsi2.flatten(), psi2.flatten()) - 0.5 * np.sum(np.square(mu) + S) + 0.5 * np.sum(log_S)
|
||||
return -float(lik)
|
||||
|
||||
def latent_grad(mu_S, kern, Z, dL_dpsi0, dL_dpsi1, dL_dpsi2):
|
||||
"""
|
||||
This is the same as latent_cost_and_grad but only for the grad
|
||||
"""
|
||||
mu, log_S = mu_S.reshape(2, 1, -1)
|
||||
S = np.exp(log_S)
|
||||
|
||||
mu0, S0 = kern.dpsi0_dmuS(dL_dpsi0, Z, mu, S)
|
||||
mu1, S1 = kern.dpsi1_dmuS(dL_dpsi1, Z, mu, S)
|
||||
mu2, S2 = kern.dpsi2_dmuS(dL_dpsi2, Z, mu, S)
|
||||
|
||||
dmu = mu0 + mu1 + mu2 - mu
|
||||
# dS = S0 + S1 + S2 -0.5 + .5/S
|
||||
dlnS = S * (S0 + S1 + S2 - 0.5) + .5
|
||||
|
||||
return -np.hstack((dmu.flatten(), dlnS.flatten()))
|
||||
|
|
|
|||
|
|
@ -2,23 +2,30 @@
|
|||
# Licensed under the BSD 3-clause license (see LICENSE.txt)
|
||||
|
||||
import numpy as np
|
||||
import itertools
|
||||
import pylab
|
||||
import itertools, logging
|
||||
|
||||
from ..core import Model
|
||||
from ..kern import Kern
|
||||
from ..core.parameterization.variational import NormalPosterior, NormalPrior
|
||||
from ..core.parameterization import Param, Parameterized
|
||||
from ..core.parameterization.observable_array import ObsAr
|
||||
from ..inference.latent_function_inference.var_dtc import VarDTCMissingData, VarDTC
|
||||
from ..inference.latent_function_inference import InferenceMethodList
|
||||
from ..likelihoods import Gaussian
|
||||
from GPy.util.initialization import initialize_latent
|
||||
from ..util.initialization import initialize_latent
|
||||
from ..core.sparse_gp import SparseGP, GP
|
||||
|
||||
class MRD(Model):
|
||||
class MRD(SparseGP):
|
||||
"""
|
||||
!WARNING: This is bleeding edge code and still in development.
|
||||
Functionality may change fundamentally during development!
|
||||
|
||||
Apply MRD to all given datasets Y in Ylist.
|
||||
|
||||
Y_i in [n x p_i]
|
||||
|
||||
If Ylist is a dictionary, the keys of the dictionary are the names, and the
|
||||
values are the different datasets to compare.
|
||||
|
||||
The samples n in the datasets need
|
||||
to match up, whereas the dimensionality p_d can differ.
|
||||
|
||||
|
|
@ -39,40 +46,77 @@ class MRD(Model):
|
|||
:param num_inducing: number of inducing inputs to use
|
||||
:param Z: initial inducing inputs
|
||||
:param kernel: list of kernels or kernel to copy for each output
|
||||
:type kernel: [GPy.kern.kern] | GPy.kern.kern | None (default)
|
||||
:param :class:`~GPy.inference.latent_function_inference inference_method: the inference method to use
|
||||
:param :class:`~GPy.likelihoods.likelihood.Likelihood` likelihood: the likelihood to use
|
||||
:type kernel: [GPy.kernels.kernels] | GPy.kernels.kernels | None (default)
|
||||
:param :class:`~GPy.inference.latent_function_inference inference_method:
|
||||
InferenceMethodList of inferences, or one inference method for all
|
||||
:param :class:`~GPy.likelihoodss.likelihoods.likelihoods` likelihoods: the likelihoods to use
|
||||
:param str name: the name of this model
|
||||
:param [str] Ynames: the names for the datasets given, must be of equal length as Ylist or None
|
||||
"""
|
||||
def __init__(self, Ylist, input_dim, X=None, X_variance=None,
|
||||
initx = 'PCA', initz = 'permute',
|
||||
num_inducing=10, Z=None, kernel=None,
|
||||
inference_method=None, likelihood=None, name='mrd', Ynames=None):
|
||||
super(MRD, self).__init__(name)
|
||||
inference_method=None, likelihoods=None, name='mrd', Ynames=None):
|
||||
super(GP, self).__init__(name)
|
||||
|
||||
self.logger = logging.getLogger(self.__class__.__name__)
|
||||
self.input_dim = input_dim
|
||||
self.num_inducing = num_inducing
|
||||
|
||||
self.Ylist = Ylist
|
||||
if isinstance(Ylist, dict):
|
||||
Ynames, Ylist = zip(*Ylist.items())
|
||||
|
||||
self.logger.debug("creating observable arrays")
|
||||
self.Ylist = [ObsAr(Y) for Y in Ylist]
|
||||
|
||||
if Ynames is None:
|
||||
self.logger.debug("creating Ynames")
|
||||
Ynames = ['Y{}'.format(i) for i in range(len(Ylist))]
|
||||
self.names = Ynames
|
||||
assert len(self.names) == len(self.Ylist), "one name per dataset, or None if Ylist is a dict"
|
||||
|
||||
if inference_method is None:
|
||||
self.inference_method= InferenceMethodList()
|
||||
warned = False
|
||||
for y in Ylist:
|
||||
inan = np.isnan(y)
|
||||
if np.any(inan):
|
||||
if not warned:
|
||||
self.logger.warn("WARNING: NaN values detected, make sure initx method can cope with NaN values or provide starting latent space X")
|
||||
warned = True
|
||||
self.inference_method.append(VarDTCMissingData(limit=1, inan=inan))
|
||||
else:
|
||||
self.inference_method.append(VarDTC(limit=1))
|
||||
self.logger.debug("created inference method <{}>".format(hex(id(self.inference_method[-1]))))
|
||||
else:
|
||||
if not isinstance(inference_method, InferenceMethodList):
|
||||
self.logger.debug("making inference_method an InferenceMethodList")
|
||||
inference_method = InferenceMethodList(inference_method)
|
||||
self.inference_method = inference_method
|
||||
|
||||
|
||||
self._in_init_ = True
|
||||
X, fracs = self._init_X(initx, Ylist)
|
||||
if X is None:
|
||||
X, fracs = self._init_X(initx, Ylist)
|
||||
else:
|
||||
fracs = [X.var(0)]*len(Ylist)
|
||||
self.Z = Param('inducing inputs', self._init_Z(initz, X))
|
||||
self.num_inducing = self.Z.shape[0] # ensure M==N if M>N
|
||||
|
||||
# sort out the kernels
|
||||
self.logger.info("building kernels")
|
||||
if kernel is None:
|
||||
from ..kern import RBF
|
||||
self.kern = [RBF(input_dim, ARD=1, lengthscale=fracs[i], name='rbf'.format(i)) for i in range(len(Ylist))]
|
||||
kernels = [RBF(input_dim, ARD=1, lengthscale=fracs[i]) for i in range(len(Ylist))]
|
||||
elif isinstance(kernel, Kern):
|
||||
self.kern = []
|
||||
kernels = []
|
||||
for i in range(len(Ylist)):
|
||||
k = kernel.copy()
|
||||
self.kern.append(k)
|
||||
kernels.append(k)
|
||||
else:
|
||||
assert len(kernel) == len(Ylist), "need one kernel per output"
|
||||
assert all([isinstance(k, Kern) for k in kernel]), "invalid kernel object detected!"
|
||||
self.kern = kernel
|
||||
kernels = kernel
|
||||
|
||||
if X_variance is None:
|
||||
X_variance = np.random.uniform(0.1, 0.2, X.shape)
|
||||
|
|
@ -80,32 +124,28 @@ class MRD(Model):
|
|||
self.variational_prior = NormalPrior()
|
||||
self.X = NormalPosterior(X, X_variance)
|
||||
|
||||
if likelihood is None:
|
||||
self.likelihood = [Gaussian(name='Gaussian_noise'.format(i)) for i in range(len(Ylist))]
|
||||
else: self.likelihood = likelihood
|
||||
|
||||
if inference_method is None:
|
||||
self.inference_method= []
|
||||
for y in Ylist:
|
||||
if np.any(np.isnan(y)):
|
||||
self.inference_method.append(VarDTCMissingData(limit=1))
|
||||
else:
|
||||
self.inference_method.append(VarDTC(limit=1))
|
||||
else:
|
||||
self.inference_method = inference_method
|
||||
self.inference_method.set_limit(len(Ylist))
|
||||
if likelihoods is None:
|
||||
likelihoods = [Gaussian(name='Gaussian_noise'.format(i)) for i in range(len(Ylist))]
|
||||
else: likelihoods = likelihoods
|
||||
|
||||
self.logger.info("adding X and Z")
|
||||
self.add_parameters(self.X, self.Z)
|
||||
|
||||
if Ynames is None:
|
||||
Ynames = ['Y{}'.format(i) for i in range(len(Ylist))]
|
||||
self.bgplvms = []
|
||||
self.num_data = Ylist[0].shape[0]
|
||||
|
||||
for i, n, k, l in itertools.izip(itertools.count(), Ynames, self.kern, self.likelihood):
|
||||
for i, n, k, l, Y in itertools.izip(itertools.count(), Ynames, kernels, likelihoods, Ylist):
|
||||
assert Y.shape[0] == self.num_data, "All datasets need to share the number of datapoints, and those have to correspond to one another"
|
||||
p = Parameterized(name=n)
|
||||
p.add_parameter(k)
|
||||
p.kern = k
|
||||
p.add_parameter(l)
|
||||
setattr(self, 'Y{}'.format(i), p)
|
||||
p.likelihood = l
|
||||
self.add_parameter(p)
|
||||
self.bgplvms.append(p)
|
||||
|
||||
self.posterior = None
|
||||
self.logger.info("init done")
|
||||
self._in_init_ = False
|
||||
|
||||
def parameters_changed(self):
|
||||
|
|
@ -113,14 +153,15 @@ class MRD(Model):
|
|||
self.posteriors = []
|
||||
self.Z.gradient[:] = 0.
|
||||
self.X.gradient[:] = 0.
|
||||
|
||||
for y, k, l, i in itertools.izip(self.Ylist, self.kern, self.likelihood, self.inference_method):
|
||||
for y, b, i in itertools.izip(self.Ylist, self.bgplvms, self.inference_method):
|
||||
self.logger.info('working on im <{}>'.format(hex(id(i))))
|
||||
k, l = b.kern, b.likelihood
|
||||
posterior, lml, grad_dict = i.inference(k, self.X, self.Z, l, y)
|
||||
|
||||
self.posteriors.append(posterior)
|
||||
self._log_marginal_likelihood += lml
|
||||
|
||||
# likelihood gradients
|
||||
# likelihoods gradients
|
||||
l.update_gradients(grad_dict.pop('dL_dthetaL'))
|
||||
|
||||
#gradients wrt kernel
|
||||
|
|
@ -133,12 +174,19 @@ class MRD(Model):
|
|||
#gradients wrt Z
|
||||
self.Z.gradient += k.gradients_X(dL_dKmm, self.Z)
|
||||
self.Z.gradient += k.gradients_Z_expectations(
|
||||
grad_dict['dL_dpsi1'], grad_dict['dL_dpsi2'], Z=self.Z, variational_posterior=self.X)
|
||||
grad_dict['dL_dpsi0'],
|
||||
grad_dict['dL_dpsi1'],
|
||||
grad_dict['dL_dpsi2'],
|
||||
Z=self.Z, variational_posterior=self.X)
|
||||
|
||||
dL_dmean, dL_dS = k.gradients_qX_expectations(variational_posterior=self.X, Z=self.Z, **grad_dict)
|
||||
self.X.mean.gradient += dL_dmean
|
||||
self.X.variance.gradient += dL_dS
|
||||
|
||||
self.posterior = self.posteriors[0]
|
||||
self.kern = self.bgplvms[0].kern
|
||||
self.likelihood = self.bgplvms[0].likelihood
|
||||
|
||||
# update for the KL divergence
|
||||
self.variational_prior.update_gradients_KL(self.X)
|
||||
self._log_marginal_likelihood -= self.variational_prior.KL_divergence(self.X)
|
||||
|
|
@ -151,7 +199,7 @@ class MRD(Model):
|
|||
Ylist = self.Ylist
|
||||
if init in "PCA_concat":
|
||||
X, fracs = initialize_latent('PCA', self.input_dim, np.hstack(Ylist))
|
||||
fracs = [fracs]*self.input_dim
|
||||
fracs = [fracs]*len(Ylist)
|
||||
elif init in "PCA_single":
|
||||
X = np.zeros((Ylist[0].shape[0], self.input_dim))
|
||||
fracs = []
|
||||
|
|
@ -162,7 +210,7 @@ class MRD(Model):
|
|||
else: # init == 'random':
|
||||
X = np.random.randn(Ylist[0].shape[0], self.input_dim)
|
||||
fracs = X.var(0)
|
||||
fracs = [fracs]*self.input_dim
|
||||
fracs = [fracs]*len(Ylist)
|
||||
X -= X.mean()
|
||||
X /= X.std()
|
||||
return X, fracs
|
||||
|
|
@ -177,10 +225,12 @@ class MRD(Model):
|
|||
return Z
|
||||
|
||||
def _handle_plotting(self, fignum, axes, plotf, sharex=False, sharey=False):
|
||||
import matplotlib.pyplot as plt
|
||||
if axes is None:
|
||||
fig = pylab.figure(num=fignum)
|
||||
fig = plt.figure(num=fignum)
|
||||
sharex_ax = None
|
||||
sharey_ax = None
|
||||
plots = []
|
||||
for i, g in enumerate(self.bgplvms):
|
||||
try:
|
||||
if sharex:
|
||||
|
|
@ -197,26 +247,36 @@ class MRD(Model):
|
|||
ax = axes[i]
|
||||
else:
|
||||
raise ValueError("Need one axes per latent dimension input_dim")
|
||||
plotf(i, g, ax)
|
||||
plots.append(plotf(i, g, ax))
|
||||
if sharey_ax is not None:
|
||||
pylab.setp(ax.get_yticklabels(), visible=False)
|
||||
pylab.draw()
|
||||
plt.setp(ax.get_yticklabels(), visible=False)
|
||||
plt.draw()
|
||||
if axes is None:
|
||||
fig.tight_layout()
|
||||
return fig
|
||||
else:
|
||||
return pylab.gcf()
|
||||
try:
|
||||
fig.tight_layout()
|
||||
except:
|
||||
pass
|
||||
return plots
|
||||
|
||||
def plot_X(self, fignum=None, ax=None):
|
||||
fig = self._handle_plotting(fignum, ax, lambda i, g, ax: ax.imshow(g.X))
|
||||
return fig
|
||||
def predict(self, Xnew, full_cov=False, Y_metadata=None, kern=None, Yindex=0):
|
||||
"""
|
||||
Prediction for data set Yindex[default=0].
|
||||
This predicts the output mean and variance for the dataset given in Ylist[Yindex]
|
||||
"""
|
||||
self.posterior = self.posteriors[Yindex]
|
||||
self.kern = self.bgplvms[0].kern
|
||||
self.likelihood = self.bgplvms[0].likelihood
|
||||
return super(MRD, self).predict(Xnew, full_cov, Y_metadata, kern)
|
||||
|
||||
def plot_predict(self, fignum=None, ax=None, sharex=False, sharey=False, **kwargs):
|
||||
fig = self._handle_plotting(fignum,
|
||||
ax,
|
||||
lambda i, g, ax: ax.imshow(g. predict(g.X)[0], **kwargs),
|
||||
sharex=sharex, sharey=sharey)
|
||||
return fig
|
||||
#===============================================================================
|
||||
# TODO: Predict! Maybe even change to several bgplvms, which share an X?
|
||||
#===============================================================================
|
||||
# def plot_predict(self, fignum=None, ax=None, sharex=False, sharey=False, **kwargs):
|
||||
# fig = self._handle_plotting(fignum,
|
||||
# ax,
|
||||
# lambda i, g, ax: ax.imshow(g.predict(g.X)[0], **kwargs),
|
||||
# sharex=sharex, sharey=sharey)
|
||||
# return fig
|
||||
|
||||
def plot_scales(self, fignum=None, ax=None, titles=None, sharex=False, sharey=True, *args, **kwargs):
|
||||
"""
|
||||
|
|
@ -228,28 +288,55 @@ class MRD(Model):
|
|||
"""
|
||||
if titles is None:
|
||||
titles = [r'${}$'.format(name) for name in self.names]
|
||||
ymax = reduce(max, [np.ceil(max(g.input_sensitivity())) for g in self.bgplvms])
|
||||
ymax = reduce(max, [np.ceil(max(g.kern.input_sensitivity())) for g in self.bgplvms])
|
||||
def plotf(i, g, ax):
|
||||
ax.set_ylim([0,ymax])
|
||||
g.kern.plot_ARD(ax=ax, title=titles[i], *args, **kwargs)
|
||||
return g.kern.plot_ARD(ax=ax, title=titles[i], *args, **kwargs)
|
||||
fig = self._handle_plotting(fignum, ax, plotf, sharex=sharex, sharey=sharey)
|
||||
return fig
|
||||
|
||||
def plot_latent(self, fignum=None, ax=None, *args, **kwargs):
|
||||
fig = self.gref.plot_latent(fignum=fignum, ax=ax, *args, **kwargs) # self._handle_plotting(fignum, ax, lambda i, g, ax: g.plot_latent(ax=ax, *args, **kwargs))
|
||||
return fig
|
||||
def plot_latent(self, labels=None, which_indices=None,
|
||||
resolution=50, ax=None, marker='o', s=40,
|
||||
fignum=None, plot_inducing=True, legend=True,
|
||||
plot_limits=None,
|
||||
aspect='auto', updates=False, predict_kwargs={}, imshow_kwargs={}):
|
||||
"""
|
||||
see plotting.matplot_dep.dim_reduction_plots.plot_latent
|
||||
if predict_kwargs is None, will plot latent spaces for 0th dataset (and kernel), otherwise give
|
||||
predict_kwargs=dict(Yindex='index') for plotting only the latent space of dataset with 'index'.
|
||||
"""
|
||||
import sys
|
||||
assert "matplotlib" in sys.modules, "matplotlib package has not been imported."
|
||||
from matplotlib import pyplot as plt
|
||||
from ..plotting.matplot_dep import dim_reduction_plots
|
||||
if "Yindex" not in predict_kwargs:
|
||||
predict_kwargs['Yindex'] = 0
|
||||
if ax is None:
|
||||
fig = plt.figure(num=fignum)
|
||||
ax = fig.add_subplot(111)
|
||||
else:
|
||||
fig = ax.figure
|
||||
plot = dim_reduction_plots.plot_latent(self, labels, which_indices,
|
||||
resolution, ax, marker, s,
|
||||
fignum, plot_inducing, legend,
|
||||
plot_limits, aspect, updates, predict_kwargs, imshow_kwargs)
|
||||
ax.set_title(self.bgplvms[predict_kwargs['Yindex']].name)
|
||||
try:
|
||||
fig.tight_layout()
|
||||
except:
|
||||
pass
|
||||
|
||||
def _debug_plot(self):
|
||||
self.plot_X_1d()
|
||||
fig = pylab.figure("MRD DEBUG PLOT", figsize=(4 * len(self.bgplvms), 9))
|
||||
fig.clf()
|
||||
axes = [fig.add_subplot(3, len(self.bgplvms), i + 1) for i in range(len(self.bgplvms))]
|
||||
self.plot_X(ax=axes)
|
||||
axes = [fig.add_subplot(3, len(self.bgplvms), i + len(self.bgplvms) + 1) for i in range(len(self.bgplvms))]
|
||||
self.plot_latent(ax=axes)
|
||||
axes = [fig.add_subplot(3, len(self.bgplvms), i + 2 * len(self.bgplvms) + 1) for i in range(len(self.bgplvms))]
|
||||
self.plot_scales(ax=axes)
|
||||
pylab.draw()
|
||||
fig.tight_layout()
|
||||
return plot
|
||||
|
||||
def __getstate__(self):
|
||||
state = super(MRD, self).__getstate__()
|
||||
del state['kern']
|
||||
del state['likelihood']
|
||||
return state
|
||||
|
||||
def __setstate__(self, state):
|
||||
# TODO:
|
||||
super(MRD, self).__setstate__(state)
|
||||
self.kern = self.bgplvms[0].kern
|
||||
self.likelihood = self.bgplvms[0].likelihood
|
||||
self.parameters_changed()
|
||||
|
|
@ -31,7 +31,7 @@ def plot_latent(model, labels=None, which_indices=None,
|
|||
resolution=50, ax=None, marker='o', s=40,
|
||||
fignum=None, plot_inducing=False, legend=True,
|
||||
plot_limits=None,
|
||||
aspect='auto', updates=False, **kwargs):
|
||||
aspect='auto', updates=False, predict_kwargs={}, imshow_kwargs={}):
|
||||
"""
|
||||
:param labels: a np.array of size model.num_data containing labels for the points (can be number, strings, etc)
|
||||
:param resolution: the resolution of the grid on which to evaluate the predictive variance
|
||||
|
|
@ -60,7 +60,7 @@ def plot_latent(model, labels=None, which_indices=None,
|
|||
def plot_function(x):
|
||||
Xtest_full = np.zeros((x.shape[0], model.X.shape[1]))
|
||||
Xtest_full[:, [input_1, input_2]] = x
|
||||
_, var = model.predict(Xtest_full)
|
||||
_, var = model.predict(Xtest_full, **predict_kwargs)
|
||||
var = var[:, :1]
|
||||
return np.log(var)
|
||||
|
||||
|
|
@ -81,7 +81,7 @@ def plot_latent(model, labels=None, which_indices=None,
|
|||
view = ImshowController(ax, plot_function,
|
||||
(xmin, ymin, xmax, ymax),
|
||||
resolution, aspect=aspect, interpolation='bilinear',
|
||||
cmap=pb.cm.binary, **kwargs)
|
||||
cmap=pb.cm.binary, **imshow_kwargs)
|
||||
|
||||
# make sure labels are in order of input:
|
||||
ulabels = []
|
||||
|
|
|
|||
|
|
@ -97,7 +97,7 @@ def plot_fit(model, plot_limits=None, which_data_rows='all',
|
|||
|
||||
for d in which_data_ycols:
|
||||
plots['gpplot'] = gpplot(Xnew, m[:, d], lower[:, d], upper[:, d], ax=ax, edgecol=linecol, fillcol=fillcol)
|
||||
plots['dataplot'] = ax.plot(X[which_data_rows,free_dims], Y[which_data_rows, d], data_symbol, mew=1.5)
|
||||
if not plot_raw: plots['dataplot'] = ax.plot(X[which_data_rows,free_dims], Y[which_data_rows, d], data_symbol, mew=1.5)
|
||||
|
||||
#optionally plot some samples
|
||||
if samples: #NOTE not tested with fixed_inputs
|
||||
|
|
@ -151,7 +151,7 @@ def plot_fit(model, plot_limits=None, which_data_rows='all',
|
|||
for d in which_data_ycols:
|
||||
m_d = m[:,d].reshape(resolution, resolution).T
|
||||
plots['contour'] = ax.contour(x, y, m_d, levels, vmin=m.min(), vmax=m.max(), cmap=pb.cm.jet)
|
||||
plots['dataplot'] = ax.scatter(X[which_data_rows, free_dims[0]], X[which_data_rows, free_dims[1]], 40, Y[which_data_rows, d], cmap=pb.cm.jet, vmin=m.min(), vmax=m.max(), linewidth=0.)
|
||||
if not plot_raw: plots['dataplot'] = ax.scatter(X[which_data_rows, free_dims[0]], X[which_data_rows, free_dims[1]], 40, Y[which_data_rows, d], cmap=pb.cm.jet, vmin=m.min(), vmax=m.max(), linewidth=0.)
|
||||
|
||||
#set the limits of the plot to some sensible values
|
||||
ax.set_xlim(xmin[0], xmax[0])
|
||||
|
|
|
|||
|
|
@ -88,7 +88,6 @@ class vector_show(matplotlib_show):
|
|||
|
||||
|
||||
class lvm(matplotlib_show):
|
||||
|
||||
def __init__(self, vals, model, data_visualize, latent_axes=None, sense_axes=None, latent_index=[0,1]):
|
||||
"""Visualize a latent variable model
|
||||
|
||||
|
|
@ -99,10 +98,11 @@ class lvm(matplotlib_show):
|
|||
"""
|
||||
if vals is None:
|
||||
if isinstance(model.X, VariationalPosterior):
|
||||
vals = param_to_array(model.X.mean)
|
||||
vals = model.X.mean.values
|
||||
else:
|
||||
vals = param_to_array(model.X)
|
||||
|
||||
vals = model.X.values
|
||||
if len(vals.shape)==1:
|
||||
vals = vals[None,:]
|
||||
matplotlib_show.__init__(self, vals, axes=latent_axes)
|
||||
|
||||
if isinstance(latent_axes,mpl.axes.Axes):
|
||||
|
|
@ -133,7 +133,7 @@ class lvm(matplotlib_show):
|
|||
|
||||
def modify(self, vals):
|
||||
"""When latent values are modified update the latent representation and ulso update the output visualization."""
|
||||
self.vals = vals.copy()
|
||||
self.vals = vals.view(np.ndarray).copy()
|
||||
y = self.model.predict(self.vals)[0]
|
||||
self.data_visualize.modify(y)
|
||||
self.latent_handle.set_data(self.vals[0,self.latent_index[0]], self.vals[0,self.latent_index[1]])
|
||||
|
|
@ -146,7 +146,6 @@ class lvm(matplotlib_show):
|
|||
pass
|
||||
|
||||
def on_click(self, event):
|
||||
print 'click!'
|
||||
if event.inaxes!=self.latent_axes: return
|
||||
self.move_on = not self.move_on
|
||||
self.called = True
|
||||
|
|
@ -219,11 +218,11 @@ class lvm_dimselect(lvm):
|
|||
self.labels = labels
|
||||
lvm.__init__(self,vals,model,data_visualize,latent_axes,sense_axes,latent_index)
|
||||
self.show_sensitivities()
|
||||
print "use left and right mouse butons to select dimensions"
|
||||
print self.latent_values
|
||||
print "use left and right mouse buttons to select dimensions"
|
||||
|
||||
|
||||
def on_click(self, event):
|
||||
|
||||
if event.inaxes==self.sense_axes:
|
||||
new_index = max(0,min(int(np.round(event.xdata-0.5)),self.model.input_dim-1))
|
||||
if event.button == 1:
|
||||
|
|
@ -249,6 +248,7 @@ class lvm_dimselect(lvm):
|
|||
|
||||
|
||||
def on_leave(self,event):
|
||||
print type(self.latent_values)
|
||||
latent_values = self.latent_values.copy()
|
||||
y = self.model.predict(latent_values[None,:])[0]
|
||||
self.data_visualize.modify(y)
|
||||
|
|
@ -393,14 +393,13 @@ class mocap_data_show_vpython(vpython_show):
|
|||
def process_values(self):
|
||||
raise NotImplementedError, "this needs to be implemented to use the data_show class"
|
||||
|
||||
|
||||
class mocap_data_show(matplotlib_show):
|
||||
"""Base class for visualizing motion capture data."""
|
||||
|
||||
def __init__(self, vals, axes=None, connect=None):
|
||||
if axes==None:
|
||||
fig = plt.figure()
|
||||
axes = fig.add_subplot(111, projection='3d')
|
||||
axes = fig.add_subplot(111, projection='3d',aspect='equal')
|
||||
matplotlib_show.__init__(self, vals, axes)
|
||||
|
||||
self.connect = connect
|
||||
|
|
@ -445,11 +444,12 @@ class mocap_data_show(matplotlib_show):
|
|||
def process_values(self):
|
||||
raise NotImplementedError, "this needs to be implemented to use the data_show class"
|
||||
|
||||
def initialize_axes(self):
|
||||
def initialize_axes(self, boundary=0.05):
|
||||
"""Set up the axes with the right limits and scaling."""
|
||||
self.x_lim = np.array([self.vals[:, 0].min(), self.vals[:, 0].max()])
|
||||
self.y_lim = np.array([self.vals[:, 1].min(), self.vals[:, 1].max()])
|
||||
self.z_lim = np.array([self.vals[:, 2].min(), self.vals[:, 2].max()])
|
||||
bs = [(self.vals[:, i].max()-self.vals[:, i].min())*boundary for i in xrange(3)]
|
||||
self.x_lim = np.array([self.vals[:, 0].min()-bs[0], self.vals[:, 0].max()+bs[0]])
|
||||
self.y_lim = np.array([self.vals[:, 1].min()-bs[1], self.vals[:, 1].max()+bs[1]])
|
||||
self.z_lim = np.array([self.vals[:, 2].min()-bs[2], self.vals[:, 2].max()+bs[2]])
|
||||
|
||||
def initialize_axes_modify(self):
|
||||
self.points_handle.remove()
|
||||
|
|
@ -472,6 +472,8 @@ class mocap_data_show(matplotlib_show):
|
|||
class stick_show(mocap_data_show):
|
||||
"""Show a three dimensional point cloud as a figure. Connect elements of the figure together using the matrix connect."""
|
||||
def __init__(self, vals, connect=None, axes=None):
|
||||
if len(vals.shape)==1:
|
||||
vals = vals[None,:]
|
||||
mocap_data_show.__init__(self, vals, axes=axes, connect=connect)
|
||||
|
||||
def process_values(self):
|
||||
|
|
|
|||
|
|
@ -94,22 +94,18 @@ class MiscTests(unittest.TestCase):
|
|||
np.testing.assert_equal(m.log_likelihood(), m2.log_likelihood())
|
||||
|
||||
m.kern.lengthscale.randomize()
|
||||
m._trigger_params_changed()
|
||||
m2.kern.lengthscale = m.kern.lengthscale
|
||||
np.testing.assert_equal(m.log_likelihood(), m2.log_likelihood())
|
||||
|
||||
m.kern.lengthscale.randomize()
|
||||
m._trigger_params_changed()
|
||||
m2['.*lengthscale'] = m.kern.lengthscale
|
||||
np.testing.assert_equal(m.log_likelihood(), m2.log_likelihood())
|
||||
|
||||
m.kern.lengthscale.randomize()
|
||||
m._trigger_params_changed()
|
||||
m2['.*lengthscale'] = m.kern['.*lengthscale']
|
||||
np.testing.assert_equal(m.log_likelihood(), m2.log_likelihood())
|
||||
|
||||
m.kern.lengthscale.randomize()
|
||||
m._trigger_params_changed()
|
||||
m2.kern.lengthscale = m.kern['.*lengthscale']
|
||||
np.testing.assert_equal(m.log_likelihood(), m2.log_likelihood())
|
||||
|
||||
|
|
@ -130,6 +126,23 @@ class MiscTests(unittest.TestCase):
|
|||
m2.kern[:] = m.kern[''].values()
|
||||
np.testing.assert_equal(m.log_likelihood(), m2.log_likelihood())
|
||||
|
||||
def test_big_model(self):
|
||||
m = GPy.examples.dimensionality_reduction.mrd_simulation(optimize=0, plot=0, plot_sim=0)
|
||||
m.X.fix()
|
||||
print m
|
||||
m.unfix()
|
||||
m.checkgrad()
|
||||
print m
|
||||
m.fix()
|
||||
print m
|
||||
m.inducing_inputs.unfix()
|
||||
print m
|
||||
m.checkgrad()
|
||||
m.unfix()
|
||||
m.checkgrad()
|
||||
m.checkgrad()
|
||||
print m
|
||||
|
||||
def test_model_set_params(self):
|
||||
m = GPy.models.GPRegression(self.X, self.Y)
|
||||
lengthscale = np.random.uniform()
|
||||
|
|
|
|||
|
|
@ -4,7 +4,8 @@ Created on 13 Mar 2014
|
|||
@author: maxz
|
||||
'''
|
||||
import unittest, itertools
|
||||
import cPickle as pickle
|
||||
#import cPickle as pickle
|
||||
import pickle
|
||||
import numpy as np
|
||||
from GPy.core.parameterization.index_operations import ParameterIndexOperations,\
|
||||
ParameterIndexOperationsView
|
||||
|
|
@ -15,8 +16,7 @@ from GPy.core.parameterization.priors import Gaussian
|
|||
from GPy.kern._src.rbf import RBF
|
||||
from GPy.kern._src.linear import Linear
|
||||
from GPy.kern._src.static import Bias, White
|
||||
from GPy.examples.dimensionality_reduction import mrd_simulation,\
|
||||
bgplvm_simulation
|
||||
from GPy.examples.dimensionality_reduction import mrd_simulation
|
||||
from GPy.examples.regression import toy_rbf_1d_50
|
||||
from GPy.core.parameterization.variational import NormalPosterior
|
||||
from GPy.models.gp_regression import GPRegression
|
||||
|
|
@ -89,6 +89,7 @@ class Test(ListDictTestCase):
|
|||
self.assertIs(pcopy.constraints, pcopy.rbf.lengthscale.constraints._param_index_ops)
|
||||
self.assertIs(pcopy.constraints, pcopy.linear.constraints._param_index_ops)
|
||||
self.assertListEqual(par.param_array.tolist(), pcopy.param_array.tolist())
|
||||
pcopy.gradient = 10 # gradient does not get copied anymore
|
||||
self.assertListEqual(par.gradient_full.tolist(), pcopy.gradient_full.tolist())
|
||||
self.assertSequenceEqual(str(par), str(pcopy))
|
||||
self.assertIsNot(par.param_array, pcopy.param_array)
|
||||
|
|
@ -125,8 +126,8 @@ class Test(ListDictTestCase):
|
|||
def test_modelrecreation(self):
|
||||
par = toy_rbf_1d_50(optimize=0, plot=0)
|
||||
pcopy = GPRegression(par.X.copy(), par.Y.copy(), kernel=par.kern.copy())
|
||||
self.assertListEqual(par.param_array.tolist(), pcopy.param_array.tolist())
|
||||
self.assertListEqual(par.gradient_full.tolist(), pcopy.gradient_full.tolist())
|
||||
np.testing.assert_allclose(par.param_array, pcopy.param_array)
|
||||
np.testing.assert_allclose(par.gradient_full, pcopy.gradient_full)
|
||||
self.assertSequenceEqual(str(par), str(pcopy))
|
||||
self.assertIsNot(par.param_array, pcopy.param_array)
|
||||
self.assertIsNot(par.gradient_full, pcopy.gradient_full)
|
||||
|
|
@ -139,7 +140,7 @@ class Test(ListDictTestCase):
|
|||
par.pickle(f)
|
||||
f.seek(0)
|
||||
pcopy = pickle.load(f)
|
||||
self.assertListEqual(par.param_array.tolist(), pcopy.param_array.tolist())
|
||||
np.testing.assert_allclose(par.param_array, pcopy.param_array)
|
||||
np.testing.assert_allclose(par.gradient_full, pcopy.gradient_full)
|
||||
self.assertSequenceEqual(str(par), str(pcopy))
|
||||
self.assert_(pcopy.checkgrad())
|
||||
|
|
@ -150,6 +151,7 @@ class Test(ListDictTestCase):
|
|||
par = NormalPosterior(X,Xv)
|
||||
par.gradient = 10
|
||||
pcopy = par.copy()
|
||||
pcopy.gradient = 10
|
||||
self.assertListEqual(par.param_array.tolist(), pcopy.param_array.tolist())
|
||||
self.assertListEqual(par.gradient_full.tolist(), pcopy.gradient_full.tolist())
|
||||
self.assertSequenceEqual(str(par), str(pcopy))
|
||||
|
|
@ -174,6 +176,7 @@ class Test(ListDictTestCase):
|
|||
self.assertSequenceEqual(str(par), str(pcopy))
|
||||
self.assertIsNot(par.param_array, pcopy.param_array)
|
||||
self.assertIsNot(par.gradient_full, pcopy.gradient_full)
|
||||
self.assertTrue(par.checkgrad())
|
||||
self.assertTrue(pcopy.checkgrad())
|
||||
self.assert_(np.any(pcopy.gradient!=0.0))
|
||||
with tempfile.TemporaryFile('w+b') as f:
|
||||
|
|
|
|||
|
|
@ -1,9 +1,7 @@
|
|||
from ..core.parameterization.parameter_core import Observable
|
||||
import itertools, collections, weakref
|
||||
import collections, weakref, logging
|
||||
|
||||
class Cacher(object):
|
||||
|
||||
|
||||
def __init__(self, operation, limit=5, ignore_args=(), force_kwargs=()):
|
||||
"""
|
||||
Parameters:
|
||||
|
|
@ -12,13 +10,15 @@ class Cacher(object):
|
|||
:param int limit: depth of cacher
|
||||
:param [int] ignore_args: list of indices, pointing at arguments to ignore in *args of operation(*args). This includes self!
|
||||
:param [str] force_kwargs: list of kwarg names (strings). If a kwarg with that name is given, the cacher will force recompute and wont cache anything.
|
||||
:param int verbose: verbosity level. 0: no print outs, 1: casual print outs, 2: debug level print outs
|
||||
"""
|
||||
self.limit = int(limit)
|
||||
self.ignore_args = ignore_args
|
||||
self.force_kwargs = force_kwargs
|
||||
self.operation=operation
|
||||
self.operation = operation
|
||||
self.order = collections.deque()
|
||||
self.cached_inputs = {} # point from cache_ids to a list of [ind_ids], which where used in cache cache_id
|
||||
self.cached_inputs = {} # point from cache_ids to a list of [ind_ids], which where used in cache cache_id
|
||||
self.logger = logging.getLogger("cache")
|
||||
|
||||
#=======================================================================
|
||||
# point from each ind_id to [ref(obj), cache_ids]
|
||||
|
|
@ -27,81 +27,109 @@ class Cacher(object):
|
|||
self.cached_input_ids = {}
|
||||
#=======================================================================
|
||||
|
||||
self.cached_outputs = {} # point from cache_ids to outputs
|
||||
self.inputs_changed = {} # point from cache_ids to bools
|
||||
self.cached_outputs = {} # point from cache_ids to outputs
|
||||
self.inputs_changed = {} # point from cache_ids to bools
|
||||
|
||||
def id(self, obj):
|
||||
"""returns the self.id of an object, to be used in caching individual self.ids"""
|
||||
return hex(id(obj))
|
||||
|
||||
def combine_inputs(self, args, kw):
|
||||
"Combines the args and kw in a unique way, such that ordering of kwargs does not lead to recompute"
|
||||
self.logger.debug("combining args and kw")
|
||||
return args + tuple(c[1] for c in sorted(kw.items(), key=lambda x: x[0]))
|
||||
|
||||
def prepare_cache_id(self, combined_args_kw, ignore_args):
|
||||
"get the cacheid (conc. string of argument ids in order) ignoring ignore_args"
|
||||
return "".join(str(id(a)) for i,a in enumerate(combined_args_kw) if i not in ignore_args)
|
||||
"get the cacheid (conc. string of argument self.ids in order) ignoring ignore_args"
|
||||
cache_id = "".join(self.id(a) for i, a in enumerate(combined_args_kw) if i not in ignore_args)
|
||||
self.logger.debug("cache_id={} was created".format(cache_id))
|
||||
return cache_id
|
||||
|
||||
def ensure_cache_length(self, cache_id):
|
||||
"Ensures the cache is within its limits and has one place free"
|
||||
self.logger.debug("cache length gets ensured")
|
||||
if len(self.order) == self.limit:
|
||||
self.logger.debug("cache limit of l={} was reached".format(self.limit))
|
||||
# we have reached the limit, so lets release one element
|
||||
cache_id = self.order.popleft()
|
||||
self.logger.debug("cach_id '{}' gets removed".format(cache_id))
|
||||
combined_args_kw = self.cached_inputs[cache_id]
|
||||
for ind in combined_args_kw:
|
||||
ind_id = id(ind)
|
||||
ref, cache_ids = self.cached_input_ids[ind_id]
|
||||
if len(cache_ids) == 1 and ref() is not None:
|
||||
ref().remove_observer(self, self.on_cache_changed)
|
||||
del self.cached_input_ids[ind_id]
|
||||
else:
|
||||
cache_ids.remove(cache_id)
|
||||
self.cached_input_ids[ind_id] = [ref, cache_ids]
|
||||
if ind is not None:
|
||||
ind_id = self.id(ind)
|
||||
tmp = self.cached_input_ids.get(ind_id, None)
|
||||
if tmp is not None:
|
||||
ref, cache_ids = tmp
|
||||
if len(cache_ids) == 1 and ref() is not None:
|
||||
ref().remove_observer(self, self.on_cache_changed)
|
||||
del self.cached_input_ids[ind_id]
|
||||
else:
|
||||
cache_ids.remove(cache_id)
|
||||
self.cached_input_ids[ind_id] = [ref, cache_ids]
|
||||
self.logger.debug("removing caches")
|
||||
del self.cached_outputs[cache_id]
|
||||
del self.inputs_changed[cache_id]
|
||||
del self.cached_inputs[cache_id]
|
||||
|
||||
def add_to_cache(self, cache_id, combined_args_kw, output):
|
||||
def add_to_cache(self, cache_id, inputs, output):
|
||||
"""This adds cache_id to the cache, with inputs and output"""
|
||||
self.inputs_changed[cache_id] = False
|
||||
self.cached_outputs[cache_id] = output
|
||||
self.order.append(cache_id)
|
||||
self.cached_inputs[cache_id] = combined_args_kw
|
||||
for a in combined_args_kw:
|
||||
ind_id = id(a)
|
||||
v = self.cached_input_ids.get(ind_id, [weakref.ref(a), []])
|
||||
v[1].append(cache_id)
|
||||
if len(v[1]) == 1:
|
||||
a.add_observer(self, self.on_cache_changed)
|
||||
self.cached_input_ids[ind_id] = v
|
||||
self.cached_inputs[cache_id] = inputs
|
||||
for a in inputs:
|
||||
if a is not None:
|
||||
ind_id = self.id(a)
|
||||
v = self.cached_input_ids.get(ind_id, [weakref.ref(a), []])
|
||||
self.logger.debug("cache_id '{}' gets stored".format(cache_id))
|
||||
v[1].append(cache_id)
|
||||
if len(v[1]) == 1:
|
||||
self.logger.debug("adding observer to object {}".format(repr(a)))
|
||||
a.add_observer(self, self.on_cache_changed)
|
||||
self.cached_input_ids[ind_id] = v
|
||||
|
||||
def __call__(self, *args, **kw):
|
||||
"""
|
||||
A wrapper function for self.operation,
|
||||
"""
|
||||
|
||||
#=======================================================================
|
||||
# !WARNING CACHE OFFSWITCH!
|
||||
# return self.operation(*args, **kw)
|
||||
#=======================================================================
|
||||
|
||||
# 1: Check whether we have forced recompute arguments:
|
||||
if len(self.force_kwargs) != 0:
|
||||
for k in self.force_kwargs:
|
||||
if k in kw and kw[k] is not None:
|
||||
return self.operation(*args, **kw)
|
||||
|
||||
# 2: prepare_cache_id and get the unique id string for this call
|
||||
|
||||
# 2: prepare_cache_id and get the unique self.id string for this call
|
||||
inputs = self.combine_inputs(args, kw)
|
||||
cache_id = self.prepare_cache_id(inputs, self.ignore_args)
|
||||
|
||||
# 2: if anything is not cachable, we will just return the operation, without caching
|
||||
if reduce(lambda a,b: a or (not isinstance(b, Observable)), inputs, False):
|
||||
if reduce(lambda a, b: a or (not (isinstance(b, Observable) or b is None)), inputs, False):
|
||||
self.logger.info("some inputs are not observable: returning without caching")
|
||||
self.logger.debug(str(map(lambda x: isinstance(x, Observable) or x is None, inputs)))
|
||||
self.logger.debug(str(map(repr, inputs)))
|
||||
return self.operation(*args, **kw)
|
||||
# 3&4: check whether this cache_id has been cached, then has it changed?
|
||||
try:
|
||||
if(self.inputs_changed[cache_id]):
|
||||
# 4: This happens, when elements have changed for this cache id
|
||||
self.logger.debug("{} already seen, but inputs changed. refreshing cacher".format(cache_id))
|
||||
# 4: This happens, when elements have changed for this cache self.id
|
||||
self.inputs_changed[cache_id] = False
|
||||
self.cached_outputs[cache_id] = self.operation(*args, **kw)
|
||||
except KeyError:
|
||||
self.logger.info("{} never seen, creating cache entry".format(cache_id))
|
||||
# 3: This is when we never saw this chache_id:
|
||||
self.ensure_cache_length(cache_id)
|
||||
self.add_to_cache(cache_id, inputs, self.operation(*args, **kw))
|
||||
except:
|
||||
self.logger.error("an error occurred while trying to run caching for {}, resetting".format(cache_id))
|
||||
self.reset()
|
||||
raise
|
||||
# 5: We have seen this cache_id and it is cached:
|
||||
self.logger.info("returning cache {}".format(cache_id))
|
||||
return self.cached_outputs[cache_id]
|
||||
|
||||
def on_cache_changed(self, direct, which=None):
|
||||
|
|
@ -110,10 +138,13 @@ class Cacher(object):
|
|||
|
||||
this function gets 'hooked up' to the inputs when we cache them, and upon their elements being changed we update here.
|
||||
"""
|
||||
for ind_id in [id(direct), id(which)]:
|
||||
_, cache_ids = self.cached_input_ids.get(ind_id, [None, []])
|
||||
for cache_id in cache_ids:
|
||||
self.inputs_changed[cache_id] = True
|
||||
for what in [direct, which]:
|
||||
if what is not None:
|
||||
ind_id = self.id(what)
|
||||
_, cache_ids = self.cached_input_ids.get(ind_id, [None, []])
|
||||
for cache_id in cache_ids:
|
||||
self.logger.info("callback from {} changed inputs from {}".format(ind_id, self.inputs_changed[cache_id]))
|
||||
self.inputs_changed[cache_id] = True
|
||||
|
||||
def reset(self):
|
||||
"""
|
||||
|
|
@ -150,7 +181,7 @@ class Cacher_wrap(object):
|
|||
return partial(self, obj)
|
||||
def __call__(self, *args, **kwargs):
|
||||
obj = args[0]
|
||||
#import ipdb;ipdb.set_trace()
|
||||
# import ipdb;ipdb.set_trace()
|
||||
try:
|
||||
caches = obj.__cachers
|
||||
except AttributeError:
|
||||
|
|
|
|||
|
|
@ -5,18 +5,19 @@ import ConfigParser
|
|||
import os
|
||||
config = ConfigParser.ConfigParser()
|
||||
|
||||
home = os.getenv('HOME') or os.getenv('USERPROFILE')
|
||||
user_file = os.path.join(home,'.gpy_config.cfg')
|
||||
default_file = os.path.abspath(os.path.join(os.path.dirname( __file__ ), '..', 'gpy_config.cfg'))
|
||||
# print user_file, os.path.isfile(user_file)
|
||||
# print default_file, os.path.isfile(default_file)
|
||||
# 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'))
|
||||
|
||||
# 1. check if the user has a ~/.gpy_config.cfg
|
||||
if os.path.isfile(user_file):
|
||||
config.read(user_file)
|
||||
elif os.path.isfile(default_file):
|
||||
# 2. if not, use the default one
|
||||
config.read(default_file)
|
||||
else:
|
||||
#3. panic
|
||||
raise ValueError, "no configuration file found"
|
||||
# These files are optional
|
||||
# This specifies configurations that are typically specific to the machine (it is found alongside the GPy installation).
|
||||
local_file = os.path.abspath(os.path.join(os.path.dirname( __file__ ), '..', 'machine.cfg'))
|
||||
|
||||
# This specifies configurations specific to the user (it is found in the user home directory)
|
||||
home = os.getenv('HOME') or os.getenv('USERPROFILE')
|
||||
user_file = os.path.join(home,'.gpy_user.cfg')
|
||||
|
||||
# Read in the given files.
|
||||
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 + "."
|
||||
|
|
|
|||
|
|
@ -57,6 +57,20 @@
|
|||
"http://www.cs.nyu.edu/~roweis/data/"
|
||||
]
|
||||
},
|
||||
"cifar-10": {
|
||||
"citation": "Learning Multiple Layers of Features from Tiny Images, Alex Krizhevsky, 2009, Tech report available here: http://www.cs.toronto.edu/~kriz/learning-features-2009-TR.pdf",
|
||||
"details": "The CIFAR-10 and CIFAR-100 are labeled subsets of the 80 million tiny images dataset. They were collected by Alex Krizhevsky, Vinod Nair, and Geoffrey Hinton. Details are available on this webpage: http://www.cs.toronto.edu/~kriz/cifar.html. The CIFAR-10 dataset consists of 60000 32x32 colour images in 10 classes, with 6000 images per class. There are 50000 training images and 10000 test images.",
|
||||
"files": [
|
||||
[
|
||||
"cifar-10-python.tar.gz"
|
||||
]
|
||||
],
|
||||
"license": null,
|
||||
"size": 0,
|
||||
"urls": [
|
||||
"http://www.cs.toronto.edu/~kriz/"
|
||||
]
|
||||
},
|
||||
"cmu_mocap_full": {
|
||||
"citation": "Please include this in your acknowledgements: The data used in this project was obtained from mocap.cs.cmu.edu.\\nThe 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.",
|
||||
|
|
@ -123,11 +137,39 @@
|
|||
]
|
||||
],
|
||||
"license": null,
|
||||
"size": 0,
|
||||
"size": 20258,
|
||||
"urls": [
|
||||
"http://staffwww.dcs.shef.ac.uk/people/N.Lawrence/dataset_mirror/drosophila_protein/"
|
||||
]
|
||||
},
|
||||
"spellman_yeast": {
|
||||
"citation": "Paul T. Spellman, Gavin Sherlock, Michael Q. Zhang, Vishwanath R. Iyer, Kirk Anders, Michael B. Eisen, Patrick O. Brown, David Botstein, and Bruce Futcher 'Comprehensive Identification of Cell Cycle-regulated Genes of the Yeast Saccharomyces cerevisiae by Microarray Hybridization.' Molecular Biology of the Cell 9, 3273-3297",
|
||||
"details": "Two colour spotted cDNA array data set of a series of experiments to identify which genes in Yeast are cell cycle regulated.",
|
||||
"files": [
|
||||
[
|
||||
"combined.txt"
|
||||
]
|
||||
],
|
||||
"license": null,
|
||||
"size": 2510955,
|
||||
"urls": [
|
||||
"http://genome-www.stanford.edu/cellcycle/data/rawdata/"
|
||||
]
|
||||
},
|
||||
"lee_yeast_ChIP": {
|
||||
"citation": "Tong Ihn Lee, Nicola J. Rinaldi, Francois Robert, Duncan T. Odom, Ziv Bar-Joseph, Georg K. Gerber, Nancy M. Hannett, Christopher T. Harbison, Craig M. Thompson, Itamar Simon, Julia Zeitlinger, Ezra G. Jennings, Heather L. Murray, D. Benjamin Gordon, Bing Ren, John J. Wyrick, Jean-Bosco Tagne, Thomas L. Volkert, Ernest Fraenkel, David K. Gifford, Richard A. Young 'Transcriptional Regulatory Networks in Saccharomyces cerevisiae' Science 298 (5594) pg 799--804. DOI: 10.1126/science.1075090",
|
||||
"details": "Binding location analysis for 106 regulators in yeast. The data consists of p-values for binding of regulators to genes derived from ChIP-chip experiments.",
|
||||
"files": [
|
||||
[
|
||||
"binding_by_gene.tsv"
|
||||
]
|
||||
],
|
||||
"license": null,
|
||||
"size": 1674161,
|
||||
"urls": [
|
||||
"http://jura.wi.mit.edu/young_public/regulatory_network/"
|
||||
]
|
||||
},
|
||||
"epomeo_gpx": {
|
||||
"citation": "",
|
||||
"details": "Five different GPS traces of the same run up Mount Epomeo in Ischia. The traces are from different sources. endomondo_1 and endomondo_2 are traces from the mobile phone app Endomondo, with a split in the middle. garmin_watch_via_endomondo is the trace from a Garmin watch, with a segment missing about 4 kilometers in. viewranger_phone and viewranger_tablet are traces from a phone and a tablet through the viewranger app. The viewranger_phone data comes from the same mobile phone as the Endomondo data (i.e. there are 3 GPS devices, but one device recorded two traces).",
|
||||
|
|
|
|||
|
|
@ -240,7 +240,7 @@ def cmu_urls_files(subj_motions, messages = True):
|
|||
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)
|
||||
os.makedirs(skel_dir)
|
||||
# Add skel file to list.
|
||||
url_required = True
|
||||
file_download.append(subjects[i] + '.asf')
|
||||
|
|
@ -367,21 +367,57 @@ def sod1_mouse(data_set='sod1_mouse'):
|
|||
if not data_available(data_set):
|
||||
download_data(data_set)
|
||||
from pandas import read_csv
|
||||
dirpath = os.path.join(data_path, data_set)
|
||||
filename = os.path.join(dirpath, 'sod1_C57_129_exprs.csv')
|
||||
dir_path = os.path.join(data_path, data_set)
|
||||
filename = os.path.join(dir_path, 'sod1_C57_129_exprs.csv')
|
||||
Y = read_csv(filename, header=0, index_col=0)
|
||||
num_repeats=4
|
||||
num_time=4
|
||||
num_cond=4
|
||||
X = 1
|
||||
return data_details_return({'X': X, 'Y': Y}, data_set)
|
||||
|
||||
def spellman_yeast(data_set='spellman_yeast'):
|
||||
if not data_available(data_set):
|
||||
download_data(data_set)
|
||||
from pandas import read_csv
|
||||
dir_path = os.path.join(data_path, data_set)
|
||||
filename = os.path.join(dir_path, 'combined.txt')
|
||||
Y = read_csv(filename, header=0, index_col=0, sep='\t')
|
||||
return data_details_return({'Y': Y}, data_set)
|
||||
|
||||
def spellman_yeast_cdc(data_set='spellman_yeast'):
|
||||
if not data_available(data_set):
|
||||
download_data(data_set)
|
||||
from pandas import read_csv
|
||||
dir_path = os.path.join(data_path, data_set)
|
||||
filename = os.path.join(dir_path, 'combined.txt')
|
||||
Y = read_csv(filename, header=0, index_col=0, sep='\t')
|
||||
t = np.asarray([10, 30, 50, 70, 80, 90, 100, 110, 120, 130, 140, 150, 170, 180, 190, 200, 210, 220, 230, 240, 250, 270, 290])
|
||||
times = ['cdc15_'+str(time) for time in t]
|
||||
Y = Y[times].T
|
||||
t = t[:, None]
|
||||
return data_details_return({'Y' : Y, 't': t, 'info': 'Time series of synchronized yeast cells from the CDC-15 experiment of Spellman et al (1998).'}, data_set)
|
||||
|
||||
def lee_yeast_ChIP(data_set='lee_yeast_ChIP'):
|
||||
if not data_available(data_set):
|
||||
download_data(data_set)
|
||||
from pandas import read_csv
|
||||
import zipfile
|
||||
dir_path = os.path.join(data_path, data_set)
|
||||
filename = os.path.join(dir_path, 'binding_by_gene.tsv')
|
||||
X = read_csv(filename, header=1, index_col=0, sep='\t')
|
||||
transcription_factors = [col for col in X.columns if col[:7] != 'Unnamed']
|
||||
annotations = X[['Unnamed: 1', 'Unnamed: 2', 'Unnamed: 3']]
|
||||
X = X[transcription_factors]
|
||||
return data_details_return({'annotations' : annotations, 'X' : X, 'transcription_factors': transcription_factors}, data_set)
|
||||
|
||||
|
||||
def fruitfly_tomancak(data_set='fruitfly_tomancak', gene_number=None):
|
||||
if not data_available(data_set):
|
||||
download_data(data_set)
|
||||
from pandas import read_csv
|
||||
dirpath = os.path.join(data_path, data_set)
|
||||
filename = os.path.join(dirpath, 'tomancak_exprs.csv')
|
||||
dir_path = os.path.join(data_path, data_set)
|
||||
filename = os.path.join(dir_path, 'tomancak_exprs.csv')
|
||||
Y = read_csv(filename, header=0, index_col=0).T
|
||||
num_repeats = 3
|
||||
num_time = 12
|
||||
|
|
@ -395,8 +431,8 @@ def drosophila_protein(data_set='drosophila_protein'):
|
|||
if not data_available(data_set):
|
||||
download_data(data_set)
|
||||
from pandas import read_csv
|
||||
dirpath = os.path.join(data_path, data_set)
|
||||
filename = os.path.join(dirpath, 'becker_et_al.csv')
|
||||
dir_path = os.path.join(data_path, data_set)
|
||||
filename = os.path.join(dir_path, 'becker_et_al.csv')
|
||||
Y = read_csv(filename, header=0)
|
||||
return data_details_return({'Y': Y}, data_set)
|
||||
|
||||
|
|
@ -404,8 +440,8 @@ def drosophila_knirps(data_set='drosophila_protein'):
|
|||
if not data_available(data_set):
|
||||
download_data(data_set)
|
||||
from pandas import read_csv
|
||||
dirpath = os.path.join(data_path, data_set)
|
||||
filename = os.path.join(dirpath, 'becker_et_al.csv')
|
||||
dir_path = os.path.join(data_path, data_set)
|
||||
filename = os.path.join(dir_path, 'becker_et_al.csv')
|
||||
# in the csv file we have facts_kni and ext_kni. We treat facts_kni as protein and ext_kni as mRNA
|
||||
df = read_csv(filename, header=0)
|
||||
t = df['t'][:,None]
|
||||
|
|
@ -426,31 +462,59 @@ def drosophila_knirps(data_set='drosophila_protein'):
|
|||
return data_details_return({'Y': Y, 'X': X}, data_set)
|
||||
|
||||
# This will be for downloading google trends data.
|
||||
def google_trends(query_terms=['big data', 'machine learning', 'data science'], data_set='google_trends'):
|
||||
"""Data downloaded from Google trends for given query terms. Warning, if you use this function multiple times in a row you get blocked due to terms of service violations."""
|
||||
# Inspired by this notebook:
|
||||
# http://nbviewer.ipython.org/github/sahuguet/notebooks/blob/master/GoogleTrends%20meet%20Notebook.ipynb
|
||||
def google_trends(query_terms=['big data', 'machine learning', 'data science'], data_set='google_trends', refresh_data=False):
|
||||
"""Data downloaded from Google trends for given query terms. Warning, if you use this function multiple times in a row you get blocked due to terms of service violations. The function will cache the result of your query, if you wish to refresh an old query set refresh_data to True. The function is inspired by this notebook: http://nbviewer.ipython.org/github/sahuguet/notebooks/blob/master/GoogleTrends%20meet%20Notebook.ipynb"""
|
||||
query_terms.sort()
|
||||
import pandas
|
||||
|
||||
# Create directory name for data
|
||||
dir_path = os.path.join(data_path,'google_trends')
|
||||
if not os.path.isdir(dir_path):
|
||||
os.makedirs(dir_path)
|
||||
dir_name = '-'.join(query_terms)
|
||||
dir_name = dir_name.replace(' ', '_')
|
||||
dir_path = os.path.join(dir_path,dir_name)
|
||||
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."
|
||||
# quote the query terms.
|
||||
quoted_terms = []
|
||||
for term in query_terms:
|
||||
quoted_terms.append(urllib2.quote(term))
|
||||
print "Query terms: ", ', '.join(query_terms)
|
||||
|
||||
# quote the query terms.
|
||||
for i, element in enumerate(query_terms):
|
||||
query_terms[i] = urllib2.quote(element)
|
||||
query = 'http://www.google.com/trends/fetchComponent?q=%s&cid=TIMESERIES_GRAPH_0&export=3' % ",".join(query_terms)
|
||||
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()
|
||||
data = urllib2.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]
|
||||
data = re.sub('new Date\((\d+),(\d+),(\d+)\)', (lambda m: '"%s-%02d-%02d"' % (m.group(1).strip(), 1+int(m.group(2)), int(m.group(3)))), data)
|
||||
timeseries = json.loads(data)
|
||||
columns = [k['label'] for k in timeseries['table']['cols']]
|
||||
rows = map(lambda x: [k['v'] for k in x['c']], timeseries['table']['rows'])
|
||||
df = pandas.DataFrame(rows, columns=columns)
|
||||
if not os.path.isdir(dir_path):
|
||||
os.makedirs(dir_path)
|
||||
|
||||
# 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]
|
||||
data = re.sub('new Date\((\d+),(\d+),(\d+)\)', (lambda m: '"%s-%02d-%02d"' % (m.group(1).strip(), 1+int(m.group(2)), int(m.group(3)))), data)
|
||||
timeseries = json.loads(data)
|
||||
#import pandas as pd
|
||||
columns = [k['label'] for k in timeseries['table']['cols']]
|
||||
rows = map(lambda x: [k['v'] for k in x['c']], timeseries['table']['rows'])
|
||||
terms = len(columns)-1
|
||||
X = np.asarray([(pb.datestr2num(row[0]), i) for i in range(terms) for row in rows ])
|
||||
Y = np.asarray([[row[i+1]] for i in range(terms) for row in rows ])
|
||||
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)
|
||||
|
||||
df = pandas.read_csv(file_name, parse_dates=[0])
|
||||
|
||||
columns = df.columns
|
||||
terms = len(query_terms)
|
||||
import datetime
|
||||
X = np.asarray([(row, i) for i in range(terms) for row in df.index])
|
||||
Y = np.asarray([[df.ix[row][query_terms[i]]] for i in range(terms) for row in df.index ])
|
||||
output_info = columns[1:]
|
||||
return data_details_return({'X': X, 'Y': Y, 'query_terms': output_info, 'info': "Data downloaded from google trends with query terms: " + ', '.join(output_info) + '.'}, data_set)
|
||||
|
||||
return data_details_return({'data frame' : df, 'X': X, 'Y': Y, 'query_terms': output_info, 'info': "Data downloaded from google trends with query terms: " + ', '.join(output_info) + '.'}, data_set)
|
||||
|
||||
# The data sets
|
||||
def oil(data_set='three_phase_oil_flow'):
|
||||
|
|
@ -594,7 +658,7 @@ def ripley_synth(data_set='ripley_prnn_data'):
|
|||
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 mauna_loa(data_set='mauna_loa', num_train=543, refresh_data=False):
|
||||
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'
|
||||
|
|
@ -635,7 +699,7 @@ def osu_run1(data_set='osu_run1', sample_every=4):
|
|||
return data_details_return({'Y': Y, 'connect' : connect}, data_set)
|
||||
|
||||
def swiss_roll_generated(num_samples=1000, sigma=0.0):
|
||||
with open(os.path.join(data_path, 'swiss_roll.pickle')) as f:
|
||||
with open(os.path.join(os.path.dirname(__file__), 'datasets', 'swiss_roll.pickle')) as f:
|
||||
data = pickle.load(f)
|
||||
Na = data['Y'].shape[0]
|
||||
perm = np.random.permutation(np.r_[:Na])[:num_samples]
|
||||
|
|
@ -687,14 +751,20 @@ def hapmap3(data_set='hapmap3'):
|
|||
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"
|
||||
if not data_available(data_set):
|
||||
download_data(data_set)
|
||||
dirpath = os.path.join(data_path,'hapmap3')
|
||||
|
||||
dir_path = os.path.join(data_path,'hapmap3')
|
||||
hapmap_file_name = 'hapmap3_r2_b36_fwd.consensus.qc.poly'
|
||||
preprocessed_data_paths = [os.path.join(dirpath,hapmap_file_name + file_name) for file_name in \
|
||||
unpacked_files = [os.path.join(dir_path, hapmap_file_name+ending) for ending in ['.ped', '.map']]
|
||||
unpacked_files_exist = reduce(lambda a, b:a and b, map(os.path.exists, unpacked_files))
|
||||
|
||||
if not unpacked_files_exist and not data_available(data_set):
|
||||
download_data(data_set)
|
||||
|
||||
preprocessed_data_paths = [os.path.join(dir_path,hapmap_file_name + file_name) for file_name in \
|
||||
['.snps.pickle',
|
||||
'.info.pickle',
|
||||
'.nan.pickle']]
|
||||
|
||||
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]"):
|
||||
|
|
@ -708,8 +778,7 @@ def hapmap3(data_set='hapmap3'):
|
|||
perc="="*int(20.*progress/100.))
|
||||
stdout.write(status); stdout.flush()
|
||||
return status
|
||||
unpacked_files = [os.path.join(dirpath, hapmap_file_name+ending) for ending in ['.ped', '.map']]
|
||||
if not reduce(lambda a,b: a and b, map(os.path.exists, unpacked_files)):
|
||||
if not unpacked_files_exist:
|
||||
status=write_status('unpacking...', 0, '')
|
||||
curr = 0
|
||||
for newfilepath in unpacked_files:
|
||||
|
|
@ -726,6 +795,7 @@ def hapmap3(data_set='hapmap3'):
|
|||
status=write_status('unpacking...', curr+12.*file_processed/(file_size), status)
|
||||
curr += 12
|
||||
status=write_status('unpacking...', curr, status)
|
||||
os.remove(filepath)
|
||||
status=write_status('reading .ped...', 25, status)
|
||||
# Preprocess data:
|
||||
snpstrnp = np.loadtxt(unpacked_files[0], dtype=str)
|
||||
|
|
@ -733,7 +803,7 @@ def hapmap3(data_set='hapmap3'):
|
|||
mapnp = np.loadtxt(unpacked_files[1], dtype=str)
|
||||
status=write_status('reading relationships.txt...', 42, status)
|
||||
# and metainfo:
|
||||
infodf = DataFrame.from_csv(os.path.join(dirpath,'./relationships_w_pops_121708.txt'), header=0, sep='\t')
|
||||
infodf = DataFrame.from_csv(os.path.join(dir_path,'./relationships_w_pops_121708.txt'), header=0, sep='\t')
|
||||
infodf.set_index('IID', inplace=1)
|
||||
status=write_status('filtering nan...', 45, status)
|
||||
snpstr = snpstrnp[:,6:].astype('S1').reshape(snpstrnp.shape[0], -1, 2)
|
||||
|
|
@ -796,14 +866,14 @@ def hapmap3(data_set='hapmap3'):
|
|||
def singlecell(data_set='singlecell'):
|
||||
if not data_available(data_set):
|
||||
download_data(data_set)
|
||||
|
||||
|
||||
from pandas import read_csv
|
||||
dirpath = os.path.join(data_path, data_set)
|
||||
filename = os.path.join(dirpath, 'singlecell.csv')
|
||||
dir_path = os.path.join(data_path, data_set)
|
||||
filename = os.path.join(dir_path, 'singlecell.csv')
|
||||
Y = read_csv(filename, header=0, index_col=0)
|
||||
genes = Y.columns
|
||||
labels = Y.index
|
||||
# data = np.loadtxt(os.path.join(dirpath, 'singlecell.csv'), delimiter=",", dtype=str)
|
||||
# data = np.loadtxt(os.path.join(dir_path, 'singlecell.csv'), delimiter=",", dtype=str)
|
||||
return data_details_return({'Y': Y, 'info' : "qPCR singlecell experiment in Mouse, measuring 48 gene expressions in 1-64 cell states. The labels have been created as in Guo et al. [2010]",
|
||||
'genes': genes, 'labels':labels,
|
||||
}, data_set)
|
||||
|
|
@ -1103,6 +1173,30 @@ def creep_data(data_set='creep_rupture'):
|
|||
X = all_data[:, features].copy()
|
||||
return data_details_return({'X': X, 'y': y}, data_set)
|
||||
|
||||
def cifar10_patches(data_set='cifar-10'):
|
||||
"""The Candian Institute for Advanced Research 10 image data set. Code for loading in this data is taken from this Boris Babenko's blog post, original code available here: http://bbabenko.tumblr.com/post/86756017649/learning-low-level-vision-feautres-in-10-lines-of-code"""
|
||||
dir_path = os.path.join(data_path, data_set)
|
||||
filename = os.path.join(dir_path, 'cifar-10-python.tar.gz')
|
||||
if not data_available(data_set):
|
||||
download_data(data_set)
|
||||
import tarfile
|
||||
# This code is from Boris Babenko's blog post.
|
||||
# http://bbabenko.tumblr.com/post/86756017649/learning-low-level-vision-feautres-in-10-lines-of-code
|
||||
tfile = tarfile.open(filename, 'r:gz')
|
||||
tfile.extractall(dir_path)
|
||||
|
||||
with open(os.path.join(dir_path, 'cifar-10-batches-py','data_batch_1'),'rb') as f:
|
||||
data = pickle.load(f)
|
||||
|
||||
images = data['data'].reshape((-1,3,32,32)).astype('float32')/255
|
||||
images = np.rollaxis(images, 1, 4)
|
||||
patches = np.zeros((0,5,5,3))
|
||||
for x in range(0,32-5,5):
|
||||
for y in range(0,32-5,5):
|
||||
patches = np.concatenate((patches, images[:,x:x+5,y:y+5,:]), axis=0)
|
||||
patches = patches.reshape((patches.shape[0],-1))
|
||||
return data_details_return({'Y': patches, "info" : "32x32 pixel patches extracted from the CIFAR-10 data by Boris Babenko to demonstrate k-means features."}, data_set)
|
||||
|
||||
def cmu_mocap_49_balance(data_set='cmu_mocap'):
|
||||
"""Load CMU subject 49's one legged balancing motion that was used by Alvarez, Luengo and Lawrence at AISTATS 2009."""
|
||||
train_motions = ['18', '19']
|
||||
|
|
|
|||
75
GPy/util/datasets/swiss_roll.pickle
Normal file
75
GPy/util/datasets/swiss_roll.pickle
Normal file
File diff suppressed because one or more lines are too long
|
|
@ -106,12 +106,14 @@ class pca(object):
|
|||
ulabels.append(lab)
|
||||
nlabels = len(ulabels)
|
||||
if colors is None:
|
||||
colors = [cmap(float(i) / nlabels) for i in range(nlabels)]
|
||||
colors = iter([cmap(float(i) / nlabels) for i in range(nlabels)])
|
||||
else:
|
||||
colors = iter(colors)
|
||||
X_ = self.project(X, self.Q)[:,dimensions]
|
||||
kwargs.update(dict(s=s))
|
||||
plots = list()
|
||||
for i, l in enumerate(ulabels):
|
||||
kwargs.update(dict(color=colors[i], marker=marker[i % len(marker)]))
|
||||
kwargs.update(dict(color=colors.next(), marker=marker[i % len(marker)]))
|
||||
plots.append(ax.scatter(*X_[labels == l, :].T, label=str(l), **kwargs))
|
||||
ax.set_xlabel(r"PC$_1$")
|
||||
ax.set_ylabel(r"PC$_2$")
|
||||
|
|
|
|||
|
|
@ -4,9 +4,9 @@
|
|||
.. moduleauthor:: Max Zwiessele <ibinbei@gmail.com>
|
||||
|
||||
'''
|
||||
__updated__ = '2013-12-02'
|
||||
__updated__ = '2014-05-21'
|
||||
|
||||
import numpy as np
|
||||
import numpy as np, logging
|
||||
|
||||
def common_subarrays(X, axis=0):
|
||||
"""
|
||||
|
|
@ -14,11 +14,11 @@ def common_subarrays(X, axis=0):
|
|||
Common subarrays are returned as a dictionary of <subarray, [index]> pairs, where
|
||||
the subarray is a tuple representing the subarray and the index is the index
|
||||
for the subarray in X, where index is the index to the remaining axis.
|
||||
|
||||
|
||||
:param :class:`np.ndarray` X: 2d array to check for common subarrays in
|
||||
:param int axis: axis to apply subarray detection over.
|
||||
When the index is 0, compare rows -- columns, otherwise.
|
||||
|
||||
|
||||
Examples:
|
||||
=========
|
||||
|
||||
|
|
@ -48,7 +48,17 @@ def common_subarrays(X, axis=0):
|
|||
assert X.ndim == 2 and axis in (0,1), "Only implemented for 2D arrays"
|
||||
subarrays = defaultdict(list)
|
||||
cnt = count()
|
||||
np.apply_along_axis(lambda x: iadd(subarrays[tuple(x)], [cnt.next()]), 1-axis, X)
|
||||
logger = logging.getLogger("common_subarrays")
|
||||
def accumulate(x, s, c):
|
||||
logger.debug("creating tuple")
|
||||
t = tuple(x)
|
||||
logger.debug("tuple done")
|
||||
col = c.next()
|
||||
iadd(s[t], [col])
|
||||
logger.debug("added col {}".format(col))
|
||||
return None
|
||||
if axis == 0: [accumulate(x, subarrays, cnt) for x in X]
|
||||
else: [accumulate(x, subarrays, cnt) for x in X.T]
|
||||
return subarrays
|
||||
|
||||
if __name__ == '__main__':
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue