diff --git a/GPy/__init__.py b/GPy/__init__.py index 94128dfa..b6bf81b7 100644 --- a/GPy/__init__.py +++ b/GPy/__init__.py @@ -4,7 +4,7 @@ import warnings warnings.filterwarnings("ignore", category=DeprecationWarning) import core -from core.parameterization import transformations +from core.parameterization import transformations, priors import models import mappings import inference @@ -15,7 +15,6 @@ import testing from numpy.testing import Tester from nose.tools import nottest import kern -from core import priors import plotting @nottest diff --git a/GPy/core/__init__.py b/GPy/core/__init__.py index 39887284..839529d6 100644 --- a/GPy/core/__init__.py +++ b/GPy/core/__init__.py @@ -2,7 +2,6 @@ # Licensed under the BSD 3-clause license (see LICENSE.txt) from model import * -from parameterization import priors from parameterization.parameterized import * from gp import GP from sparse_gp import SparseGP diff --git a/GPy/core/model.py b/GPy/core/model.py index c944f943..58008a99 100644 --- a/GPy/core/model.py +++ b/GPy/core/model.py @@ -21,8 +21,6 @@ class Model(Parameterized): _allowed_failures = 10 # number of allowed failures def __init__(self, name): super(Model, self).__init__(name) # Parameterized.__init__(self) - self.priors = [] - self._priors = ParameterIndexOperations() self.optimization_runs = [] self.sampling_runs = [] self.preferred_optimizer = 'scg' @@ -67,66 +65,66 @@ class Model(Parameterized): self.priors = state.pop() Parameterized._setstate(self, state) - def set_prior(self, regexp, what): - """ - - Sets priors on the model parameters. - - **Notes** - - Asserts that the prior is suitable for the constraint. If the - wrong constraint is in place, an error is raised. If no - constraint is in place, one is added (warning printed). - - For tied parameters, the prior will only be "counted" once, thus - a prior object is only inserted on the first tied index - - :param regexp: regular expression of parameters on which priors need to be set. - :type param: string, regexp, or integer array - :param what: prior to set on parameter. - :type what: GPy.core.Prior type - - """ - if self.priors is None: - self.priors = [None for i in range(self._get_params().size)] - - which = self.grep_param_names(regexp) - - # check tied situation - tie_partial_matches = [tie for tie in self.tied_indices if (not set(tie).isdisjoint(set(which))) & (not set(tie) == set(which))] - if len(tie_partial_matches): - raise ValueError, "cannot place prior across partial ties" - tie_matches = [tie for tie in self.tied_indices if set(which) == set(tie) ] - if len(tie_matches) > 1: - raise ValueError, "cannot place prior across multiple ties" - elif len(tie_matches) == 1: - which = which[:1] # just place a prior object on the first parameter - - - # check constraints are okay - - if what.domain is _POSITIVE: - constrained_positive_indices = [i for i, t in zip(self.constrained_indices, self.constraints) if t.domain is _POSITIVE] - if len(constrained_positive_indices): - constrained_positive_indices = np.hstack(constrained_positive_indices) - else: - constrained_positive_indices = np.zeros(shape=(0,)) - bad_constraints = np.setdiff1d(self.all_constrained_indices(), constrained_positive_indices) - assert not np.any(which[:, None] == bad_constraints), "constraint and prior incompatible" - unconst = np.setdiff1d(which, constrained_positive_indices) - if len(unconst): - print "Warning: constraining parameters to be positive:" - print '\n'.join([n for i, n in enumerate(self._get_param_names()) if i in unconst]) - print '\n' - self.constrain_positive(unconst) - elif what.domain is _REAL: - assert not np.any(which[:, None] == self.all_constrained_indices()), "constraint and prior incompatible" - else: - raise ValueError, "prior not recognised" - - # store the prior in a local list - for w in which: - self.priors[w] = what +# def set_prior(self, regexp, what): +# """ +# +# Sets priors on the model parameters. +# +# **Notes** +# +# Asserts that the prior is suitable for the constraint. If the +# wrong constraint is in place, an error is raised. If no +# constraint is in place, one is added (warning printed). +# +# For tied parameters, the prior will only be "counted" once, thus +# a prior object is only inserted on the first tied index +# +# :param regexp: regular expression of parameters on which priors need to be set. +# :type param: string, regexp, or integer array +# :param what: prior to set on parameter. +# :type what: GPy.core.Prior type +# +# """ +# if self.priors is None: +# self.priors = [None for i in range(self._get_params().size)] +# +# which = self.grep_param_names(regexp) +# +# # check tied situation +# tie_partial_matches = [tie for tie in self.tied_indices if (not set(tie).isdisjoint(set(which))) & (not set(tie) == set(which))] +# if len(tie_partial_matches): +# raise ValueError, "cannot place prior across partial ties" +# tie_matches = [tie for tie in self.tied_indices if set(which) == set(tie) ] +# if len(tie_matches) > 1: +# raise ValueError, "cannot place prior across multiple ties" +# elif len(tie_matches) == 1: +# which = which[:1] # just place a prior object on the first parameter +# +# +# # check constraints are okay +# +# if what.domain is _POSITIVE: +# constrained_positive_indices = [i for i, t in zip(self.constrained_indices, self.constraints) if t.domain is _POSITIVE] +# if len(constrained_positive_indices): +# constrained_positive_indices = np.hstack(constrained_positive_indices) +# else: +# constrained_positive_indices = np.zeros(shape=(0,)) +# bad_constraints = np.setdiff1d(self.all_constrained_indices(), constrained_positive_indices) +# assert not np.any(which[:, None] == bad_constraints), "constraint and prior incompatible" +# unconst = np.setdiff1d(which, constrained_positive_indices) +# if len(unconst): +# print "Warning: constraining parameters to be positive:" +# print '\n'.join([n for i, n in enumerate(self._get_param_names()) if i in unconst]) +# print '\n' +# self.constrain_positive(unconst) +# elif what.domain is _REAL: +# assert not np.any(which[:, None] == self.all_constrained_indices()), "constraint and prior incompatible" +# else: +# raise ValueError, "prior not recognised" +# +# # store the prior in a local list +# for w in which: +# self.priors[w] = what def get_gradient(self, name, return_names=False): """ @@ -146,22 +144,6 @@ class Model(Parameterized): else: raise AttributeError, "no parameter matches %s" % name - def log_prior(self): - """evaluate the prior""" - if self.priors is not None: - return np.sum([p.lnpdf(x) for p, x in zip(self.priors, self._get_params()) if p is not None]) - else: - return 0. - - def _log_prior_gradients(self): - """evaluate the gradients of the priors""" - if self.priors is None: - return 0. - x = self._get_params() - ret = np.zeros(x.size) - [np.put(ret, i, p.lnpdf_grad(xx)) for i, (p, xx) in enumerate(zip(self.priors, x)) if not p is None] - return ret - def randomize(self): """ Randomize the model. diff --git a/GPy/core/parameterization/index_operations.py b/GPy/core/parameterization/index_operations.py index b816e05f..d5061be7 100644 --- a/GPy/core/parameterization/index_operations.py +++ b/GPy/core/parameterization/index_operations.py @@ -107,7 +107,7 @@ class ParameterIndexOperations(object): def add(self, prop, indices): try: self._properties[prop] = combine_indices(self._properties[prop], indices) - except KeyError: + except KeyError: self._properties[prop] = indices def remove(self, prop, indices): @@ -215,10 +215,10 @@ class ParameterIndexOperationsView(object): def remove(self, prop, indices): removed = self._param_index_ops.remove(prop, indices+self._offset) - if removed.size > 0: - return removed - self._size + 1 if self[prop].size == 0: del self[prop] + if removed.size > 0: + return removed - self._size + 1 return removed diff --git a/GPy/core/parameterization/param.py b/GPy/core/parameterization/param.py index 97fcfd32..a52f2e18 100644 --- a/GPy/core/parameterization/param.py +++ b/GPy/core/parameterization/param.py @@ -10,6 +10,7 @@ from array_core import ObservableArray, ParamList __constraints_name__ = "Constraint" __index_name__ = "Index" __tie_name__ = "Tied to" +__priors_name__ = "Prior" __precision__ = numpy.get_printoptions()['precision'] # numpy printing precision used, sublassing numpy ndarray after all __print_threshold__ = 5 ###### @@ -77,9 +78,10 @@ class Param(ObservableArray, Constrainable, Gradcheckable, Indexable, Parameteri self._name = getattr(obj, 'name', None) self.gradient = getattr(obj, 'gradient', None) self.constraints = getattr(obj, 'constraints', None) + self.priors = getattr(obj, 'priors', None) - def __array_wrap__(self, out_arr, context=None): - return out_arr.view(numpy.ndarray) + #def __array_wrap__(self, out_arr, context=None): + # return out_arr.view(numpy.ndarray) #=========================================================================== # Pickling operations #=========================================================================== @@ -118,154 +120,17 @@ class Param(ObservableArray, Constrainable, Gradcheckable, Indexable, Parameteri #=========================================================================== # get/set parameters #=========================================================================== - def _set_params(self, param, update=True): self.flat = param - self._notify_tied_parameters() + #self._notify_tied_parameters() self._notify_observers() def _get_params(self): return self.flat -# @property -# def name(self): -# """ -# Name of this parameter. -# This can be a callable without parameters. The callable will be called -# every time the name property is accessed. -# """ -# if callable(self.name): -# return self.name() -# return self.name -# @name.setter -# def name(self, new_name): -# from_name = self.name -# self.name = new_name -# self._direct_parent_._name_changed(self, from_name) + def _collect_gradient(self, target): target[:] = self.gradient.flat - #=========================================================================== - # Tying operations -> bugged, TODO - #=========================================================================== - def tie_to(self, param): - """ - :param param: the parameter object to tie this parameter to. - Can be ParamConcatenation (retrieved by regexp search) - Tie this parameter to the given parameter. - Broadcasting is not allowed, but you can tie a whole dimension to - one parameter: self[:,0].tie_to(other), where other is a one-value - parameter. - - Note: For now only one parameter can have ties, so all of a parameter - will be removed, when re-tieing! - """ - #Note: this method will tie to the parameter which is the last in - # the chain of ties. Thus, if you tie to a tied parameter, - # this tie will be created to the parameter the param is tied - # to. - - assert isinstance(param, Param), "Argument {1} not of type {0}".format(Param, param.__class__) - param = numpy.atleast_1d(param) - if param.size != 1: - raise NotImplementedError, "Broadcast tying is not implemented yet" - try: - if self._original_: - self[:] = param - else: # this happens when indexing created a copy of the array - self._direct_parent_._get_original(self)[self._current_slice_] = param - except ValueError: - raise ValueError("Trying to tie {} with shape {} to {} with shape {}".format(self.name, self.shape, param.name, param.shape)) - if param is self: - raise RuntimeError, 'Cyclic tieing is not allowed' -# if len(param._tied_to_) > 0: -# if (self._direct_parent_._get_original(self) is param._direct_parent_._get_original(param) -# and len(set(self._raveled_index())&set(param._tied_to_[0]._raveled_index()))!=0): -# raise RuntimeError, 'Cyclic tieing is not allowed' -# self.tie_to(param._tied_to_[0]) -# return - if not param in self._direct_parent_._get_original(self)._tied_to_: - self._direct_parent_._get_original(self)._tied_to_ += [param] - param._add_tie_listener(self) - self._highest_parent_._set_fixed(self) - cs = self._highest_parent_._constraints_for(param, param._raveled_index()) - for cs in self._highest_parent_._constraints_for(param, param._raveled_index()): - [self.constrain(c, warning=False) for c in cs] -# for t in self._tied_to_me_.keys(): -# if t is not self: -# t.untie(self) -# t.tie_to(param) - - def untie(self, *ties): - """ - remove all ties. - """ - [t._direct_parent_._get_original(t)._remove_tie_listener(self) for t in self._tied_to_] - new_ties = [] - for t in self._direct_parent_._get_original(self)._tied_to_: - for tied in t._tied_to_me_.keys(): - if t._parent_index_ is tied._parent_index_: - new_ties.append(tied) - self._direct_parent_._get_original(self)._tied_to_ = new_ties - self._direct_parent_._get_original(self)._highest_parent_._set_unfixed(self) -# self._direct_parent_._remove_tie(self, *params) - def _notify_tied_parameters(self): - for tied, ind in self._tied_to_me_.iteritems(): - tied._on_tied_parameter_changed(self.base, list(ind)) - def _add_tie_listener(self, tied_to_me): - for t in self._tied_to_me_.keys(): - if tied_to_me._parent_index_ is t._parent_index_: - t_rav_i = t._raveled_index() - tr_rav_i = tied_to_me._raveled_index() - new_index = list(set(t_rav_i) | set(tr_rav_i)) - tmp = t._direct_parent_._get_original(t)[numpy.unravel_index(new_index, t._realshape_)] - self._tied_to_me_[tmp] = self._tied_to_me_[t] | set(self._raveled_index()) - del self._tied_to_me_[t] - return - self._tied_to_me_[tied_to_me] = set(self._raveled_index()) - def _remove_tie_listener(self, to_remove): - for t in self._tied_to_me_.keys(): - if t._parent_index_ == to_remove._parent_index_: - t_rav_i = t._raveled_index() - tr_rav_i = to_remove._raveled_index() - import ipdb;ipdb.set_trace() - new_index = list(set(t_rav_i) - set(tr_rav_i)) - if new_index: - tmp = t._direct_parent_._get_original(t)[numpy.unravel_index(new_index, t._realshape_)] - self._tied_to_me_[tmp] = self._tied_to_me_[t] - del self._tied_to_me_[t] - if len(self._tied_to_me_[tmp]) == 0: - del self._tied_to_me_[tmp] - else: - del self._tied_to_me_[t] - def _on_tied_parameter_changed(self, val, ind): - if not self._updated_: # not fast_array_equal(self, val[ind]): - val = numpy.atleast_1d(val) - self._updated_ = True - if self._original_: - self.__setitem__(slice(None), val[ind], update=False) - else: # this happens when indexing created a copy of the array - self._direct_parent_._get_original(self).__setitem__(self._current_slice_, val[ind], update=False) - self._notify_tied_parameters() - self._updated_ = False - #=========================================================================== - # Prior Operations - #=========================================================================== - def set_prior(self, prior): - """ - :param prior: prior to be set for this parameter - - Set prior for this parameter. - """ - if not hasattr(self._highest_parent_, '_set_prior'): - raise AttributeError("Parent of type {} does not support priors".format(self._highest_parent_.__class__)) - self._highest_parent_._set_prior(self, prior) - def unset_prior(self, *priors): - """ - :param priors: priors to remove from this parameter - - Remove all priors from this parameter - """ - self._highest_parent_._remove_prior(self, *priors) #=========================================================================== # Array operations -> done #=========================================================================== @@ -283,6 +148,7 @@ class Param(ObservableArray, Constrainable, Gradcheckable, Indexable, Parameteri self._notify_tied_parameters() if update: self._highest_parent_.parameters_changed() + #=========================================================================== # Index Operations: #=========================================================================== @@ -328,8 +194,9 @@ class Param(ObservableArray, Constrainable, Gradcheckable, Indexable, Parameteri return numpy.r_[a] return numpy.r_[:b] return itertools.imap(f, itertools.izip_longest(slice_index[:self._realndim_], self._realshape_, fillvalue=slice(self.size))) + #=========================================================================== - # Convienience + # Convenience #=========================================================================== @property def is_fixed(self): @@ -338,11 +205,10 @@ class Param(ObservableArray, Constrainable, Gradcheckable, Indexable, Parameteri view = super(Param, self).round(decimals, out).view(Param) view.__array_finalize__(self) return view - def _has_fixes(self): - return False round.__doc__ = numpy.round.__doc__ def _get_original(self, param): return self + #=========================================================================== # Printing -> done #=========================================================================== @@ -362,6 +228,9 @@ class Param(ObservableArray, Constrainable, Gradcheckable, Indexable, Parameteri def _constraints_str(self): return [' '.join(map(lambda c: str(c[0]) if c[1].size == self._realsize_ else "{" + str(c[0]) + "}", self.constraints.iteritems()))] @property + def _priors_str(self): + return [' '.join(map(lambda c: str(c[0]) if c[1].size == self._realsize_ else "{" + str(c[0]) + "}", self.priors.iteritems()))] + @property def _ties_str(self): return [t._short() for t in self._tied_to_] or [''] @property @@ -385,8 +254,6 @@ class Param(ObservableArray, Constrainable, Gradcheckable, Indexable, Parameteri if len(ind) != 1: ties[i, matches[0][ind_rav_matches]] = numpy.take(tt_rav_index, matches[1], mode='wrap')[ind_rav_matches] else: ties[i, matches[0]] = numpy.take(tt_rav_index, matches[1], mode='wrap') return map(lambda a: sum(a, []), zip(*[[[tie.flatten()] if tx != None else [] for tx in t] for t, tie in zip(ties, self._tied_to_)])) - def _constraints_for(self, rav_index): - return self.constraints.properties_for(rav_index) def _indices(self, slice_index=None): # get a int-array containing all indices in the first axis. if slice_index is None: @@ -404,6 +271,7 @@ class Param(ObservableArray, Constrainable, Gradcheckable, Indexable, Parameteri return numpy.fromiter(itertools.product(*expanded_index), dtype=[('', int)] * self._realndim_, count=reduce(lambda a, b: a * b.size, expanded_index, 1)).view((int, self._realndim_)) def _max_len_names(self, gen, header): + gen = map(lambda x: " ".join(map(str, x)), gen) return reduce(lambda a, b:max(a, len(b)), gen, len(header)) def _max_len_values(self): return reduce(lambda a, b:max(a, len("{x:=.{0}g}".format(__precision__, x=b))), self.flat, len(self.name_hirarchical)) @@ -418,21 +286,26 @@ class Param(ObservableArray, Constrainable, Gradcheckable, Indexable, Parameteri if ind.size > 4: indstr = ','.join(map(str, ind[:2])) + "..." + ','.join(map(str, ind[-2:])) else: indstr = ','.join(map(str, ind)) return name + '[' + indstr + ']' - def __str__(self, constr_matrix=None, indices=None, ties=None, lc=None, lx=None, li=None, lt=None): + def __str__(self, constr_matrix=None, indices=None, prirs=None, ties=None, lc=None, lx=None, li=None, lp=None, lt=None, only_name=False): filter_ = self._current_slice_ vals = self.flat if indices is None: indices = self._indices(filter_) ravi = self._raveled_index(filter_) - if constr_matrix is None: constr_matrix = self._constraints_for(ravi) + if constr_matrix is None: constr_matrix = self.constraints.properties_for(ravi) + if prirs is None: prirs = self.priors.properties_for(ravi) if ties is None: ties = self._ties_for(ravi) ties = [' '.join(map(lambda x: x._short(), t)) for t in ties] if lc is None: lc = self._max_len_names(constr_matrix, __constraints_name__) if lx is None: lx = self._max_len_values() if li is None: li = self._max_len_index(indices) if lt is None: lt = self._max_len_names(ties, __tie_name__) - header = " {i:^{2}s} | \033[1m{x:^{1}s}\033[0;0m | {c:^{0}s} | {t:^{3}s}".format(lc, lx, li, lt, x=self.name_hirarchical, c=__constraints_name__, i=__index_name__, t=__tie_name__) # nice header for printing + if lp is None: lp = self._max_len_names(prirs, __tie_name__) + sep = '-' + header_format = " {i:{5}^{2}s} | \033[1m{x:{5}^{1}s}\033[0;0m | {c:{5}^{0}s} | {p:{5}^{4}s} | {t:{5}^{3}s}" + if only_name: header = header_format.format(lc, lx, li, lt, lp, ' ', x=self.name_hirarchical, c=sep*lc, i=sep*li, t=sep*lt, p=sep*lp) # nice header for printing + else: header = header_format.format(lc, lx, li, lt, lp, ' ', x=self.name_hirarchical, c=__constraints_name__, i=__index_name__, t=__tie_name__, p=__priors_name__) # nice header for printing if not ties: ties = itertools.cycle(['']) - return "\n".join([header] + [" {i!s:^{3}s} | {x: >{1}.{2}g} | {c:^{0}s} | {t:^{4}s} ".format(lc, lx, __precision__, li, lt, x=x, c=" ".join(map(str, c)), t=(t or ''), i=i) for i, x, c, t in itertools.izip(indices, vals, constr_matrix, ties)]) # return all the constraints with right indices + return "\n".join([header] + [" {i!s:^{3}s} | {x: >{1}.{2}g} | {c:^{0}s} | {p:^{5}s} | {t:^{4}s} ".format(lc, lx, __precision__, li, lt, lp, x=x, c=" ".join(map(str, c)), p=" ".join(map(str, p)), t=(t or ''), i=i) for i, x, c, t, p in itertools.izip(indices, vals, constr_matrix, ties, prirs)]) # return all the constraints with right indices # except: return super(Param, self).__str__() class ParamConcatenation(object): @@ -538,53 +411,20 @@ class ParamConcatenation(object): def __str__(self, *args, **kwargs): def f(p): ind = p._raveled_index() - return p._constraints_for(ind), p._ties_for(ind) + return p.constraints.properties_for(ind), p._ties_for(ind), p.priors.properties_for(ind) params = self.params - constr_matrices, ties_matrices = zip(*map(f, params)) + constr_matrices, ties_matrices, prior_matrices = zip(*map(f, params)) indices = [p._indices() for p in params] lc = max([p._max_len_names(cm, __constraints_name__) for p, cm in itertools.izip(params, constr_matrices)]) lx = max([p._max_len_values() for p in params]) li = max([p._max_len_index(i) for p, i in itertools.izip(params, indices)]) lt = max([p._max_len_names(tm, __tie_name__) for p, tm in itertools.izip(params, ties_matrices)]) - strings = [p.__str__(cm, i, tm, lc, lx, li, lt) for p, cm, i, tm in itertools.izip(params,constr_matrices,indices,ties_matrices)] + lp = max([p._max_len_names(pm, __constraints_name__) for p, pm in itertools.izip(params, prior_matrices)]) + strings = [] + start = True + for p, cm, i, tm, pm in itertools.izip(params,constr_matrices,indices,ties_matrices,prior_matrices): + strings.append(p.__str__(constr_matrix=cm, indices=i, prirs=pm, ties=tm, lc=lc, lx=lx, li=li, lp=lp, lt=lt, only_name=(1-start))) + start = False return "\n".join(strings) - return "\n{}\n".format(" -"+"- | -".join(['-'*l for l in [li,lx,lc,lt]])).join(strings) def __repr__(self): - return "\n".join(map(repr,self.params)) - -if __name__ == '__main__': - - - from GPy.core.parameterized import Parameterized - from GPy.core.parameter import Param - - #X = numpy.random.randn(2,3,1,5,2,4,3) - X = numpy.random.randn(3,2) - print "random done" - p = Param("q_mean", X) - p1 = Param("q_variance", numpy.random.rand(*p.shape)) - p2 = Param("Y", numpy.random.randn(p.shape[0], 1)) - - p3 = Param("variance", numpy.random.rand()) - p4 = Param("lengthscale", numpy.random.rand(2)) - - m = Parameterized() - rbf = Parameterized(name='rbf') - - rbf.add_parameter(p3,p4) - m.add_parameter(p,p1,rbf) - - print "setting params" - #print m.q_v[3:5,[1,4,5]] - print "constraining variance" - #m[".*variance"].constrain_positive() - #print "constraining rbf" - #m.rbf_l.constrain_positive() - #m.q_variance[1,[0,5,11,19,2]].tie_to(m.rbf_v) - #m.rbf_v.tie_to(m.rbf_l[0]) - #m.rbf_l[0].tie_to(m.rbf_l[1]) - #m.q_v.tie_to(m.rbf_v) -# m.rbf_l.tie_to(m.rbf_va) - # pt = numpy.array(params._get_params_transformed()) - # ptr = numpy.random.randn(*pt.shape) -# params.X.tie_to(params.rbf_v) + return "\n".join(map(repr,self.params)) \ No newline at end of file diff --git a/GPy/core/parameterization/parameter_core.py b/GPy/core/parameterization/parameter_core.py index 7ad228f7..7e5fea73 100644 --- a/GPy/core/parameterization/parameter_core.py +++ b/GPy/core/parameterization/parameter_core.py @@ -29,7 +29,25 @@ class Parameterizable(object): if add_name: return [adjust_name_for_printing(self.name) + "." + xi for x in self._parameters_ for xi in x.parameter_names(add_name=True)] return [xi for x in self._parameters_ for xi in x.parameter_names(add_name=True)] - + + def _collect_gradient(self, target): + import itertools + [p._collect_gradient(target[s]) for p, s in itertools.izip(self._parameters_, self._param_slices_)] + + def _get_params(self): + import numpy as np + # don't overwrite this anymore! + if not self.size: + return np.empty(shape=(0,), dtype=np.float64) + return np.hstack([x._get_params() for x in self._parameters_ if x.size > 0]) + + def _set_params(self, params, update=True): + # don't overwrite this anymore! + import itertools + [p._set_params(params[s], update=update) for p, s in itertools.izip(self._parameters_, self._param_slices_)] + self.parameters_changed() + + def parameters_changed(self): """ This method gets called when parameters have changed. @@ -130,17 +148,18 @@ class Indexable(object): that is an int array, containing the indexes for the flattened param inside this parameterized logic. """ - raise NotImplementedError, "shouldnt happen, raveld index transformation required from non parameterization object?" - + raise NotImplementedError, "shouldnt happen, raveld index transformation required from non parameterization object?" + class Constrainable(Nameable, Indexable, Parameterizable): def __init__(self, name, default_constraint=None): super(Constrainable,self).__init__(name) self._default_constraint_ = default_constraint from index_operations import ParameterIndexOperations self.constraints = ParameterIndexOperations() + self.priors = ParameterIndexOperations() if self._default_constraint_ is not None: self.constrain(self._default_constraint_) - + #=========================================================================== # Fixing Parameters: #=========================================================================== @@ -186,16 +205,41 @@ class Constrainable(Nameable, Indexable, Parameterizable): self._fixes_[fixed_indices] = FIXED else: self._fixes_ = None + + def _has_fixes(self): + return hasattr(self, "_fixes_") and self._fixes_ is not None + #=========================================================================== + # Prior Operations + #=========================================================================== + def set_prior(self, prior, warning=True, update=True): + repriorized = self.unset_priors() + self._add_to_index_operations(self.priors, repriorized, prior, warning, update) + + def unset_priors(self, *priors): + return self._remove_from_index_operations(self.priors, priors) + + def log_prior(self): + """evaluate the prior""" + import numpy as np + if self.priors.size > 0: + x = self._get_params() + return np.sum([p.lnpdf(x[ind]) for p, ind in self.priors.iteritems()]) + return 0. + + def _log_prior_gradients(self): + """evaluate the gradients of the priors""" + import numpy as np + if self.priors.size > 0: + x = self._get_params() + ret = np.zeros(x.size) + [np.put(ret, ind, p.lnpdf_grad(x[ind])) for p, ind in self.priors.iteritems()] + return ret + return 0. + #=========================================================================== # Constrain operations -> done #=========================================================================== - def _parent_changed(self, parent): - from index_operations import ParameterIndexOperationsView - self.constraints = ParameterIndexOperationsView(parent.constraints, parent._offset_for(self), self.size) - self._fixes_ = None - for p in self._parameters_: - p._parent_changed(parent) def constrain(self, transform, warning=True, update=True): """ @@ -209,11 +253,7 @@ class Constrainable(Nameable, Indexable, Parameterizable): if isinstance(transform, Transformation): self._set_params(transform.initialize(self._get_params()), update=False) reconstrained = self.unconstrain() - self.constraints.add(transform, self._raveled_index()) - if warning and reconstrained.size > 0: - print "WARNING: reconstraining parameters {}".format(self.parameter_names() or self.name) - if update: - self._highest_parent_.parameters_changed() + self._add_to_index_operations(self.constraints, reconstrained, transform, warning, update) def unconstrain(self, *transforms): """ @@ -222,16 +262,7 @@ class Constrainable(Nameable, Indexable, Parameterizable): remove all :py:class:`GPy.core.transformations.Transformation` transformats of this parameter object. """ - if len(transforms) == 0: - transforms = self.constraints.properties() - import numpy as np - removed = np.empty((0,),dtype=int) - for t in transforms: - unconstrained = self.constraints.remove(t, self._raveled_index()) - removed = np.union1d(removed, unconstrained) - if t is __fixed__: - self._highest_parent_._set_unfixed(unconstrained) - return removed + return self._remove_from_index_operations(self.constraints, transforms) def constrain_positive(self, warning=True, update=True): """ @@ -277,3 +308,35 @@ class Constrainable(Nameable, Indexable, Parameterizable): Remove (lower, upper) bounded constrain from this parameter/ """ self.unconstrain(Logistic(lower, upper)) + + def _parent_changed(self, parent): + from index_operations import ParameterIndexOperationsView + self.constraints = ParameterIndexOperationsView(parent.constraints, parent._offset_for(self), self.size) + self.priors = ParameterIndexOperationsView(parent.priors, parent._offset_for(self), self.size) + self._fixes_ = None + for p in self._parameters_: + p._parent_changed(parent) + + def _add_to_index_operations(self, which, reconstrained, transform, warning, update): + if warning and reconstrained.size > 0: + print "WARNING: reconstraining parameters {}".format(self.parameter_names() or self.name) + which.add(transform, self._raveled_index()) + if update: + self._highest_parent_.parameters_changed() + + def _remove_from_index_operations(self, which, transforms): + if len(transforms) == 0: + transforms = which.properties() + import numpy as np + removed = np.empty((0, ), dtype=int) + for t in transforms: + unconstrained = which.remove(t, self._raveled_index()) + removed = np.union1d(removed, unconstrained) + if t is __fixed__: + self._highest_parent_._set_unfixed(unconstrained) + + return removed + + + + diff --git a/GPy/core/parameterization/parameterized.py b/GPy/core/parameterization/parameterized.py index 593c37e9..c8a841c0 100644 --- a/GPy/core/parameterization/parameterized.py +++ b/GPy/core/parameterization/parameterized.py @@ -66,9 +66,6 @@ class Parameterized(Constrainable, Pickleable, Observable, Gradcheckable): self._added_names_ = set() del self._in_init_ - def _has_fixes(self): - return hasattr(self, "_fixes_") and self._fixes_ is not None - def add_parameter(self, param, index=None): """ :param parameters: the parameters to add @@ -88,11 +85,14 @@ class Parameterized(Constrainable, Pickleable, Observable, Gradcheckable): # make sure the size is set if index is None: self.constraints.update(param.constraints, self.size) + self.priors.update(param.priors, self.size) self._parameters_.append(param) else: start = sum(p.size for p in self._parameters_[:index]) self.constraints.shift(start, param.size) + self.priors.shift(start, param.size) self.constraints.update(param.constraints, start) + self.priors.update(param.priors, start) self._parameters_.insert(index, param) self.size += param.size else: @@ -210,6 +210,7 @@ class Parameterized(Constrainable, Pickleable, Observable, Gradcheckable): """ return [ self._fixes_, + self.priors, self.constraints, self._parameters_, self._name, @@ -220,9 +221,10 @@ class Parameterized(Constrainable, Pickleable, Observable, Gradcheckable): self._added_names_ = state.pop() self._name = state.pop() self._parameters_ = state.pop() - self._connect_parameters() self.constraints = state.pop() + self.priors = state.pop() self._fixes_ = state.pop() + self._connect_parameters() self.parameters_changed() #=========================================================================== # Gradient control @@ -248,16 +250,6 @@ class Parameterized(Constrainable, Pickleable, Observable, Gradcheckable): if self._has_fixes(): return n[self._fixes_] return n - def _get_params(self): - # don't overwrite this anymore! - if not self.size: - return np.empty(shape=(0,), dtype=np.float64) - return numpy.hstack([x._get_params() for x in self._parameters_ if x.size > 0]) - - def _set_params(self, params, update=True): - # don't overwrite this anymore! - [p._set_params(params[s], update=update) for p, s in itertools.izip(self._parameters_, self._param_slices_)] - self.parameters_changed() def _get_params_transformed(self): # transformed parameters (apply transformation rules) p = self._get_params() @@ -348,36 +340,6 @@ class Parameterized(Constrainable, Pickleable, Observable, Gradcheckable): return self._direct_parent_.hirarchy_name() + adjust_name_for_printing(self.name) + "." return '' #=========================================================================== - # Constraint Handling: - #=========================================================================== - #=========================================================================== - # def _add_constrain(self, param, transform, warning=True): - # rav_i = self._raveled_index_for(param) - # reconstrained = self._remove_constrain(param, index=rav_i) # remove constraints before - # # if removing constraints before adding new is not wanted, just delete the above line! - # self.constraints.add(transform, rav_i) - # param = self._get_original(param) - # if not (transform == __fixed__): - # param._set_params(transform.initialize(param._get_params()), update=False) - # if warning and any(reconstrained): - # # if you want to print the whole params object, which was reconstrained use: - # # m = str(param[self._backtranslate_index(param, reconstrained)]) - # print "Warning: re-constraining parameters:\n{}".format(param._short()) - # return rav_i - # def _remove_constrain(self, param, *transforms, **kwargs): - # if not transforms: - # transforms = self.constraints.properties() - # removed_indices = numpy.array([]).astype(int) - # if "index" in kwargs: index = kwargs['index'] - # else: index = self._raveled_index_for(param) - # for constr in transforms: - # removed = self.constraints.remove(constr, index) - # if constr is __fixed__: - # self._set_unfixed(removed) - # removed_indices = numpy.union1d(removed_indices, removed) - # return removed_indices - #=========================================================================== - #=========================================================================== # Get/set parameters: #=========================================================================== def grep_param_names(self, regexp): @@ -434,9 +396,6 @@ class Parameterized(Constrainable, Pickleable, Observable, Gradcheckable): return self._direct_parent_.hirarchy_name() + adjust_name_for_printing(self.name) else: return adjust_name_for_printing(self.name) - #parameter_names = property(parameter_names, doc="Names for all parameters handled by this parameterization object -- will add hirarchy name entries for printing") - def _collect_gradient(self, target): - [p._collect_gradient(target[s]) for p, s in itertools.izip(self._parameters_, self._param_slices_)] @property def flattened_parameters(self): return [xi for x in self._parameters_ for xi in x.flattened_parameters] @@ -455,6 +414,9 @@ class Parameterized(Constrainable, Pickleable, Observable, Gradcheckable): def _constraints_str(self): return [cs for p in self._parameters_ for cs in p._constraints_str] @property + def _priors_str(self): + return [cs for p in self._parameters_ for cs in p._priors_str] + @property def _description_str(self): return [xi for x in self._parameters_ for xi in x._description_str] @property @@ -463,20 +425,23 @@ class Parameterized(Constrainable, Pickleable, Observable, Gradcheckable): def __str__(self, header=True): name = adjust_name_for_printing(self.name) + "." - constrs = self._constraints_str; ts = self._ties_str + constrs = self._constraints_str; + ts = self._ties_str + prirs = self._priors_str desc = self._description_str; names = self.parameter_names() nl = max([len(str(x)) for x in names + [name]]) sl = max([len(str(x)) for x in desc + ["Value"]]) cl = max([len(str(x)) if x else 0 for x in constrs + ["Constraint"]]) tl = max([len(str(x)) if x else 0 for x in ts + ["Tied to"]]) - format_spec = " \033[1m{{name:<{0}s}}\033[0;0m | {{desc:^{1}s}} | {{const:^{2}s}} | {{t:^{3}s}}".format(nl, sl, cl, tl) + pl = max([len(str(x)) if x else 0 for x in prirs + ["Prior"]]) + format_spec = " \033[1m{{name:<{0}s}}\033[0;0m | {{desc:^{1}s}} | {{const:^{2}s}} | {{pri:^{3}s}} | {{t:^{4}s}}".format(nl, sl, cl, pl, tl) to_print = [] - for n, d, c, t in itertools.izip(names, desc, constrs, ts): - to_print.append(format_spec.format(name=n, desc=d, const=c, t=t)) + for n, d, c, t, p in itertools.izip(names, desc, constrs, ts, prirs): + to_print.append(format_spec.format(name=n, desc=d, const=c, t=t, pri=p)) # to_print = [format_spec.format(p=p, const=c, t=t) if isinstance(p, Param) else p.__str__(header=False) for p, c, t in itertools.izip(self._parameters_, constrs, ts)] - sep = '-' * (nl + sl + cl + tl + 8 * 2 + 3) + sep = '-' * (nl + sl + cl + + pl + tl + 8 * 2 + 3) if header: - header = " {{0:<{0}s}} | {{1:^{1}s}} | {{2:^{2}s}} | {{3:^{3}s}}".format(nl, sl, cl, tl).format(name, "Value", "Constraint", "Tied to") + header = " {{0:<{0}s}} | {{1:^{1}s}} | {{2:^{2}s}} | {{3:^{3}s}} | {{4:^{4}s}}".format(nl, sl, cl, pl, tl).format(name, "Value", "Constraint", "Prior", "Tied to") # header += '\n' + sep to_print.insert(0, header) return '\n'.format(sep).join(to_print) diff --git a/GPy/core/parameterization/priors.py b/GPy/core/parameterization/priors.py index f1208f18..3155be64 100644 --- a/GPy/core/parameterization/priors.py +++ b/GPy/core/parameterization/priors.py @@ -9,15 +9,16 @@ from domains import _REAL, _POSITIVE import warnings import weakref -class Prior: +class Prior(object): domain = None def pdf(self, x): return np.exp(self.lnpdf(x)) def plot(self): + import sys assert "matplotlib" in sys.modules, "matplotlib package has not been imported." - from ..plotting.matplot_dep import priors_plots + from ...plotting.matplot_dep import priors_plots priors_plots.univariate_plot(self) @@ -150,6 +151,7 @@ class MultivariateGaussian: return np.random.multivariate_normal(self.mu, self.var, n) def plot(self): + import sys assert "matplotlib" in sys.modules, "matplotlib package has not been imported." from ..plotting.matplot_dep import priors_plots priors_plots.multivariate_plot(self)