all parameterization stuff now in seperate module -> GPy.core.parameterization

GPy/core/parameterization/__init__.py

@@ -0,0 +1,5 @@
+# Copyright (c) 2012, GPy authors (see AUTHORS.txt).
+# Licensed under the BSD 3-clause license (see LICENSE.txt)
+
+from param import Param, ObservableArray
+from parameterized import Parameterized

GPy/core/parameterization/array_core.py

@@ -0,0 +1,50 @@
+# Copyright (c) 2012, GPy authors (see AUTHORS.txt).
+# Licensed under the BSD 3-clause license (see LICENSE.txt)
+
+__updated__ = '2013-12-16'
+
+import numpy as np
+from parameter_core import Observable
+
+class ListArray(np.ndarray):
+    """
+    ndarray which can be stored in lists and checked if it is in.
+    WARNING: This overrides the functionality of x==y!!!
+    Use numpy.equal(x,y) for element-wise equality testing.
+    """
+    def __new__(cls, input_array):
+        obj = np.asanyarray(input_array).view(cls)
+        return obj
+    def __eq__(self, other):
+        return other is self
+
+class ObservableArray(ListArray, Observable):
+    """
+    An ndarray which reports changes to its observers.
+    The observers can add themselves with a callable, which
+    will be called every time this array changes. The callable
+    takes exactly one argument, which is this array itself.
+    """
+    def __new__(cls, input_array):
+        obj = super(ObservableArray, cls).__new__(cls, input_array).view(cls)
+        obj._observers_ = {}
+        return obj
+    def __array_finalize__(self, obj):
+        # see InfoArray.__array_finalize__ for comments
+        if obj is None: return
+        self._observers_ = getattr(obj, '_observers_', None)
+    def __setitem__(self, s, val, update=True):
+        if self.ndim:
+            if not np.all(np.equal(self[s], val)):
+                super(ObservableArray, self).__setitem__(s, val)
+                if update:
+                    self._notify_observers()
+        else:
+            if not np.all(np.equal(self, val)):
+                super(ObservableArray, self).__setitem__(Ellipsis, val)
+                if update:
+                    self._notify_observers()
+    def __getslice__(self, start, stop):
+        return self.__getitem__(slice(start, stop))
+    def __setslice__(self, start, stop, val):
+        return self.__setitem__(slice(start, stop), val)  

GPy/core/parameterization/domains.py

@@ -0,0 +1,25 @@
+# Copyright (c) 2012, GPy authors (see AUTHORS.txt).
+# Licensed under the BSD 3-clause license (see LICENSE.txt)
+
+'''
+(Hyper-)Parameter domains defined for :py:mod:`~GPy.core.priors` and :py:mod:`~GPy.kern`.
+These domains specify the legitimate realm of the parameters to live in.
+
+:const:`~GPy.core.domains._REAL` :
+    real domain, all values in the real numbers are allowed
+
+:const:`~GPy.core.domains._POSITIVE`:
+    positive domain, only positive real values are allowed
+    
+:const:`~GPy.core.domains._NEGATIVE`:
+    same as :const:`~GPy.core.domains._POSITIVE`, but only negative values are allowed
+    
+:const:`~GPy.core.domains._BOUNDED`:
+    only values within the bounded range are allowed,
+    the bounds are specified withing the object with the bounded range
+'''
+
+_REAL = 'real'
+_POSITIVE = "positive"
+_NEGATIVE = 'negative'
+_BOUNDED = 'bounded'

GPy/core/parameterization/index_operations.py

@@ -0,0 +1,188 @@
+'''
+Created on Oct 2, 2013
+
+@author: maxzwiessele
+'''
+import numpy
+from numpy.lib.function_base import vectorize
+from param import Param
+from collections import defaultdict
+
+class ParamDict(defaultdict):
+    def __init__(self):
+        """
+        Default will be self._default, if not set otherwise
+        """
+        defaultdict.__init__(self, self.default_factory)
+        
+    def __getitem__(self, key):
+        try:
+            return defaultdict.__getitem__(self, key)
+        except KeyError:
+            for a in self.iterkeys():
+                if numpy.all(a==key) and a._parent_index_==key._parent_index_:
+                    return defaultdict.__getitem__(self, a)
+            raise        
+        
+    def __contains__(self, key):
+        if defaultdict.__contains__(self, key):
+            return True
+        for a in self.iterkeys():
+            if numpy.all(a==key) and a._parent_index_==key._parent_index_:
+                return True
+        return False
+
+    def __setitem__(self, key, value):
+        if isinstance(key, Param):
+            for a in self.iterkeys():
+                if numpy.all(a==key) and a._parent_index_==key._parent_index_:
+                    return super(ParamDict, self).__setitem__(a, value)
+        defaultdict.__setitem__(self, key, value)
+
+class SetDict(ParamDict):
+    def default_factory(self):
+        return set()
+
+class IntArrayDict(ParamDict):
+    def default_factory(self):
+        return numpy.int_([])
+
+class ParameterIndexOperations(object):
+    '''
+    Index operations for storing param index _properties
+    This class enables index with slices retrieved from object.__getitem__ calls.
+    Adding an index will add the selected indexes by the slice of an indexarray
+    indexing a shape shaped array to the flattened index array. Remove will
+    remove the selected slice indices from the flattened array.
+    You can give an offset to set an offset for the given indices in the
+    index array, for multi-param handling.
+    '''
+    def __init__(self):
+        self._properties = ParamDict()
+        #self._reverse = collections.defaultdict(list)
+        
+    def __getstate__(self):
+        return self._properties#, self._reverse
+        
+    def __setstate__(self, state):
+        self._properties = state[0]
+        # self._reverse = state[1]
+
+    def iteritems(self):
+        return self._properties.iteritems()
+    
+    def properties(self):
+        return self._properties.keys()
+
+    def iter_properties(self):
+        return self._properties.iterkeys()
+    
+    def clear(self):
+        self._properties.clear()
+    
+    def size(self):
+        return reduce(lambda a,b: a+b.size, self.iterindices(), 0)    
+    
+    def iterindices(self):
+        return self._properties.itervalues()
+    
+    def indices(self):
+        return self._properties.values()
+
+    def properties_for(self, index):
+#         already_seen = dict()
+#         for ni in index:
+#             if ni not in already_seen:
+#                 already_seen[ni] = [prop for prop in self.iter_properties() if ni in self._properties[prop]] 
+#             yield already_seen[ni]
+        return vectorize(lambda i: [prop for prop in self.iter_properties() if i in self._properties[prop]], otypes=[list])(index)
+        
+    def add(self, prop, indices):
+        try:
+            self._properties[prop] = combine_indices(self._properties[prop], indices)
+        except KeyError: 
+            self._properties[prop] = indices
+    
+    def remove(self, prop, indices):
+        if prop in self._properties:
+            diff = remove_indices(self[prop], indices)
+            removed = numpy.intersect1d(self[prop], indices, True)
+            if not index_empty(diff):
+                self._properties[prop] = diff
+            else:
+                del self._properties[prop]
+            #[self._reverse[i].remove(prop) for i in removed if prop in self._reverse[i]] 
+            return removed.astype(int)
+#         else:
+#             for a in self.properties(): 
+#                 if numpy.all(a==prop) and a._parent_index_ == prop._parent_index_:
+#                     ind = create_raveled_indices(indices, shape, offset)
+#                     diff = remove_indices(self[a], ind)
+#                     removed = numpy.intersect1d(self[a], ind, True)
+#                     if not index_empty(diff):
+#                         self._properties[a] = diff
+#                     else:
+#                         del self._properties[a]
+#                     [self._reverse[i].remove(a) for i in removed if a in self._reverse[i]] 
+#                     return removed.astype(int)
+        return numpy.array([]).astype(int)
+    def __getitem__(self, prop):
+        return self._properties[prop]
+       
+# class TieIndexOperations(object):
+#     def __init__(self, params):
+#         self.params = params
+#         self.tied_from = ParameterIndexOperations()
+#         self.tied_to = ParameterIndexOperations()
+#     def add(self, tied_from, tied_to):
+#         rav_from = self.params._raveled_index_for(tied_from)
+#         rav_to = self.params._raveled_index_for(tied_to)
+#         self.tied_from.add(tied_to, rav_from)
+#         self.tied_to.add(tied_to, rav_to)
+#         return rav_from, rav_to
+#     def remove(self, tied_from, tied_to):
+#         rav_from = self.params._raveled_index_for(tied_from)
+#         rav_to = self.params._raveled_index_for(tied_to)
+#         rem_from = self.tied_from.remove(tied_to, rav_from)
+#         rem_to = self.tied_to.remove(tied_to, rav_to)
+#         left_from = self.tied_from._properties.pop(tied_to)
+#         left_to = self.tied_to._properties.pop(tied_to)
+#         self.tied_from[numpy.delete(tied_to, rem_from)] = left_from
+#         self.tied_to[numpy.delete(tied_to, rem_to)] = left_to
+#         return rav_from, rav_to
+#     def from_to_for(self, index):
+#         return self.tied_from.properties_for(index), self.tied_to.properties_for(index)
+#     def iter_from_to_indices(self):
+#         for k, f in self.tied_from.iteritems():
+#             yield f, self.tied_to[k]
+#     def iter_to_indices(self):
+#         return self.tied_to.iterindices()
+#     def iter_from_indices(self):
+#         return self.tied_from.iterindices()
+#     def iter_from_items(self):
+#         for f, i in self.tied_from.iteritems():
+#             yield f, i
+#     def iter_properties(self):
+#         return self.tied_from.iter_properties()
+#     def properties(self):
+#         return self.tied_from.properties()
+#     def from_to_indices(self, param):
+#         return self.tied_from[param], self.tied_to[param]
+#     
+# # def create_raveled_indices(index, shape, offset=0):
+# #     if isinstance(index, (tuple, list)): i = [slice(None)] + list(index)
+# #     else: i = [slice(None), index]
+# #     ind = numpy.array(numpy.ravel_multi_index(numpy.indices(shape)[i], shape)).flat + numpy.int_(offset)
+# #     return ind
+
+def combine_indices(arr1, arr2):
+    return numpy.union1d(arr1, arr2)
+
+def remove_indices(arr, to_remove):
+    return numpy.setdiff1d(arr, to_remove, True)
+
+def index_empty(index):
+    return numpy.size(index) == 0 
+
+
+

GPy/core/parameterization/param.py

@@ -0,0 +1,593 @@
+# Copyright (c) 2012, GPy authors (see AUTHORS.txt).
+# Licensed under the BSD 3-clause license (see LICENSE.txt)
+
+import itertools
+import numpy
+from parameter_core import Constrainable, adjust_name_for_printing
+from array_core import ObservableArray
+
+###### printing
+__constraints_name__ = "Constraint"
+__index_name__ = "Index"
+__tie_name__ = "Tied to"
+__precision__ = numpy.get_printoptions()['precision'] # numpy printing precision used, sublassing numpy ndarray after all
+__print_threshold__ = 5
+######      
+
+class Float(numpy.float64, Constrainable):
+    def __init__(self, f, base):
+        super(Float,self).__init__(f)
+        self._base = base
+        
+        
+class Param(ObservableArray, Constrainable):
+    """
+    Parameter object for GPy models.
+
+    :param name:        name of the parameter to be printed
+    :param input_array: array which this parameter handles
+    :param gradient:    callable with one argument, which is the model of this parameter
+    :param args:        additional arguments to gradient
+    :param kwargs:      additional keyword arguments to gradient
+    
+    You can add/remove constraints by calling constrain on the parameter itself, e.g:
+    
+        - self[:,1].constrain_positive()
+        - self[0].tie_to(other)
+        - self.untie()
+        - self[:3,:].unconstrain()
+        - self[1].fix()
+        
+    Fixing parameters will fix them to the value they are right now. If you change
+    the fixed value, it will be fixed to the new value!
+    
+    See :py:class:`GPy.core.parameterized.Parameterized` for more details on constraining etc.
+
+    This ndarray can be stored in lists and checked if it is in.
+
+    >>> import numpy as np
+    >>> x = np.random.normal(size=(10,3))
+    >>> x in [[1], x, [3]]
+    True
+    
+    WARNING: This overrides the functionality of x==y!!!
+    Use numpy.equal(x,y) for element-wise equality testing.
+    """
+    __array_priority__ = 0 # Never give back Param
+    _fixes_ = None
+    def __new__(cls, name, input_array, *args, **kwargs):
+        obj = numpy.atleast_1d(super(Param, cls).__new__(cls, input_array=input_array))
+        obj._current_slice_ = (slice(obj.shape[0]),)
+        obj._realshape_ = obj.shape
+        obj._realsize_ = obj.size
+        obj._realndim_ = obj.ndim
+        obj._updated_ = False
+        from index_operations import SetDict
+        obj._tied_to_me_ = SetDict()
+        obj._tied_to_ = []
+        obj._original_ = True
+        return obj
+
+    def __init__(self, name, input_array):
+        super(Param, self).__init__(name=name)
+        
+    def __array_finalize__(self, obj):
+        # see InfoArray.__array_finalize__ for comments
+        if obj is None: return
+        super(Param, self).__array_finalize__(obj)
+        self._direct_parent_ = getattr(obj, '_direct_parent_', None)
+        self._parent_index_ = getattr(obj, '_parent_index_', None)
+        self._highest_parent_ = getattr(obj, '_highest_parent_', None)
+        self._current_slice_ = getattr(obj, '_current_slice_', None)
+        self._tied_to_me_ = getattr(obj, '_tied_to_me_', None)
+        self._tied_to_ = getattr(obj, '_tied_to_', None)
+        self._realshape_ = getattr(obj, '_realshape_', None)
+        self._realsize_ = getattr(obj, '_realsize_', None)
+        self._realndim_ = getattr(obj, '_realndim_', None)
+        self._updated_ = getattr(obj, '_updated_', None)
+        self._original_ = getattr(obj, '_original_', None)
+        self._name = getattr(obj, 'name', None)
+        
+    def __array_wrap__(self, out_arr, context=None):
+        return out_arr.view(numpy.ndarray)
+    #===========================================================================
+    # Pickling operations
+    #===========================================================================
+    def __reduce_ex__(self):
+        func, args, state = super(Param, self).__reduce__()
+        return func, args, (state, 
+                            (self.name,
+                             self._direct_parent_,
+                             self._parent_index_,
+                             self._highest_parent_,
+                             self._current_slice_,
+                             self._realshape_,
+                             self._realsize_,
+                             self._realndim_,
+                             self._tied_to_me_,
+                             self._tied_to_,
+                             self._updated_,
+                            )
+                            )
+
+    def __setstate__(self, state):
+        super(Param, self).__setstate__(state[0])
+        state = list(state[1])
+        self._updated_ = state.pop()
+        self._tied_to_ = state.pop()
+        self._tied_to_me_ = state.pop()
+        self._realndim_ = state.pop()
+        self._realsize_ = state.pop()
+        self._realshape_ = state.pop()
+        self._current_slice_ = state.pop()
+        self._highest_parent_ = state.pop()
+        self._parent_index_ = state.pop()
+        self._direct_parent_ = state.pop()
+        self.name = state.pop()
+
+    #===========================================================================
+    # get/set parameters
+    #===========================================================================
+
+    def _set_params(self, param, update=True):
+        self.flat = param
+        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)
+    @property
+    def _parameters_(self):
+        return []
+    #===========================================================================
+    # Fixing Parameters:
+    #===========================================================================
+    def constrain_fixed(self, warning=True):
+        """
+        Constrain this paramter to be fixed to the current value it carries.
+        
+        :param warning: print a warning for overwriting constraints.
+        """
+        self._highest_parent_._fix(self,warning)
+    fix = constrain_fixed
+    def unconstrain_fixed(self):
+        """
+        This parameter will no longer be fixed.
+        """
+        self._highest_parent_._unfix(self)
+    unfix = unconstrain_fixed
+    #===========================================================================
+    # 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
+    #===========================================================================
+    def __getitem__(self, s, *args, **kwargs):
+        if not isinstance(s, tuple):
+            s = (s,)
+        if not reduce(lambda a,b: a or numpy.any(b is Ellipsis), s, False) and len(s) <= self.ndim:
+            s += (Ellipsis,)
+        new_arr = super(Param, self).__getitem__(s, *args, **kwargs)
+        try: new_arr._current_slice_ = s; new_arr._original_ = self.base is new_arr.base
+        except AttributeError: pass# returning 0d array or float, double etc
+        return new_arr
+    def __setitem__(self, s, val, update=True):
+        super(Param, self).__setitem__(s, val, update=update)
+        self._notify_tied_parameters()
+        if update:
+            self._highest_parent_.parameters_changed()
+    #===========================================================================
+    # Index Operations:
+    #===========================================================================
+    def _internal_offset(self):
+        internal_offset = 0
+        extended_realshape = numpy.cumprod((1,) + self._realshape_[:0:-1])[::-1]
+        for i, si in enumerate(self._current_slice_[:self._realndim_]):
+            if numpy.all(si == Ellipsis):
+                continue
+            if isinstance(si, slice):
+                a = si.indices(self._realshape_[i])[0] 
+            elif isinstance(si, (list,numpy.ndarray,tuple)):
+                a = si[0]
+            else: a = si
+            if a<0:
+                a = self._realshape_[i]+a
+            internal_offset += a * extended_realshape[i]
+        return internal_offset
+    def _raveled_index(self, slice_index=None):
+        # return an index array on the raveled array, which is formed by the current_slice
+        # of this object
+        extended_realshape = numpy.cumprod((1,) + self._realshape_[:0:-1])[::-1]
+        ind = self._indices(slice_index)
+        if ind.ndim < 2: ind=ind[:,None]
+        return numpy.asarray(numpy.apply_along_axis(lambda x: numpy.sum(extended_realshape*x), 1, ind), dtype=int)
+    def _expand_index(self, slice_index=None):
+        # this calculates the full indexing arrays from the slicing objects given by get_item for _real..._ attributes
+        # it basically translates slices to their respective index arrays and turns negative indices around
+        # it tells you in the second return argument if it has only seen arrays as indices
+        if slice_index is None:
+            slice_index = self._current_slice_
+        def f(a):
+            a, b = a
+            if a not in (slice(None), Ellipsis):
+                if isinstance(a, slice):
+                    start, stop, step = a.indices(b)
+                    return numpy.r_[start:stop:step]
+                elif isinstance(a, (list,numpy.ndarray,tuple)):
+                    a = numpy.asarray(a, dtype=int)
+                    a[a<0] = b + a[a<0]
+                elif a<0:
+                    a = b+a
+                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
+    #===========================================================================
+    @property
+    def is_fixed(self):
+        return self._highest_parent_._is_fixed(self)
+    def round(self, decimals=0, out=None):
+        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
+    #===========================================================================
+    @property
+    def _description_str(self):
+        if self.size <= 1: return ["%f"%self]
+        else: return [str(self.shape)]
+    def _parameter_names(self, add_name):
+        return [self.name]
+    @property
+    def flattened_parameters(self):
+        return [self]
+    @property
+    def parameter_shapes(self):
+        return [self.shape]
+    @property
+    def _constraints_str(self):
+        return [' '.join(map(lambda c: str(c[0]) if c[1].size==self._realsize_ else "{"+str(c[0])+"}", self._highest_parent_._constraints_iter_items(self)))]
+    @property
+    def _ties_str(self):
+        return [t._short() for t in self._tied_to_] or ['']
+    @property
+    def name_hirarchical(self):
+        if self.has_parent():
+            return self._direct_parent_.hirarchy_name()+adjust_name_for_printing(self.name)
+        return adjust_name_for_printing(self.name)
+    def __repr__(self, *args, **kwargs):
+        name = "\033[1m{x:s}\033[0;0m:\n".format(
+                            x=self.name_hirarchical)
+        return name + super(Param, self).__repr__(*args,**kwargs)
+    def _ties_for(self, rav_index):
+        size = sum(p.size for p in self._tied_to_)
+        ties = numpy.empty(shape=(len(self._tied_to_), numpy.size(rav_index)), dtype=Param)
+        for i, tied_to in enumerate(self._tied_to_):
+            for t, ind in tied_to._tied_to_me_.iteritems():
+                if t._parent_index_ == self._parent_index_:
+                    matches = numpy.where(rav_index[:,None] == t._raveled_index()[None, :])
+                    tt_rav_index = tied_to._raveled_index()
+                    ind_rav_matches = numpy.where(tt_rav_index == numpy.array(list(ind)))[0]
+                    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._highest_parent_._constraints_for(self, rav_index)
+    def _indices(self, slice_index=None):
+        # get a int-array containing all indices in the first axis.
+        if slice_index is None:
+            slice_index = self._current_slice_
+        if isinstance(slice_index, (tuple, list)):
+            clean_curr_slice = [s for s in slice_index if numpy.any(s != Ellipsis)]
+            if (all(isinstance(n, (numpy.ndarray, list, tuple)) for n in clean_curr_slice) 
+                and len(set(map(len,clean_curr_slice))) <= 1):
+                return numpy.fromiter(itertools.izip(*clean_curr_slice),
+                    dtype=[('',int)]*self._realndim_,count=len(clean_curr_slice[0])).view((int, self._realndim_))
+        expanded_index = list(self._expand_index(slice_index))
+        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):
+        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))
+    def _max_len_index(self, ind):
+        return reduce(lambda a, b:max(a, len(str(b))), ind, len(__index_name__))
+    def _short(self):
+        # short string to print
+        name = self._direct_parent_.hirarchy_name() + adjust_name_for_printing(self.name)
+        if self._realsize_ < 2:
+            return name
+        ind = self._indices()
+        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):
+        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 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 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
+        #except: return super(Param, self).__str__()
+
+class ParamConcatenation(object):
+    def __init__(self, params):
+        """
+        Parameter concatenation for convienience of printing regular expression matched arrays
+        you can index this concatenation as if it was the flattened concatenation
+        of all the parameters it contains, same for setting parameters (Broadcasting enabled).
+
+        See :py:class:`GPy.core.parameter.Param` for more details on constraining.
+        """
+        #self.params = params
+        self.params = []
+        for p in params:
+            for p in p.flattened_parameters:
+                if p not in self.params:
+                    self.params.append(p)           
+        self._param_sizes = [p.size for p in self.params]
+        startstops = numpy.cumsum([0] + self._param_sizes)
+        self._param_slices_ = [slice(start, stop) for start,stop in zip(startstops, startstops[1:])]
+    #===========================================================================
+    # Get/set items, enable broadcasting
+    #===========================================================================
+    def __getitem__(self, s):
+        ind = numpy.zeros(sum(self._param_sizes), dtype=bool); ind[s] = True; 
+        params = [p._get_params()[ind[ps]] for p,ps in zip(self.params, self._param_slices_) if numpy.any(p._get_params()[ind[ps]])]
+        if len(params)==1: return params[0]
+        return ParamConcatenation(params)
+    def __setitem__(self, s, val, update=True):
+        ind = numpy.zeros(sum(self._param_sizes), dtype=bool); ind[s] = True; 
+        vals = self._vals(); vals[s] = val; del val
+        [numpy.place(p, ind[ps], vals[ps]) and p._notify_tied_parameters() 
+         for p, ps in zip(self.params, self._param_slices_)]
+        if update:
+            self.params[0]._highest_parent_.parameters_changed()
+    def _vals(self):
+        return numpy.hstack([p._get_params() for p in self.params])
+    #===========================================================================
+    # parameter operations:
+    #===========================================================================
+    def constrain(self, constraint, warning=True):
+        [param.constrain(constraint) for param in self.params]
+    constrain.__doc__ = Param.constrain.__doc__
+    def constrain_positive(self, warning=True):
+        [param.constrain_positive(warning) for param in self.params]
+    constrain_positive.__doc__ = Param.constrain_positive.__doc__
+    def constrain_fixed(self, warning=True):
+        [param.constrain_fixed(warning) for param in self.params]
+    constrain_fixed.__doc__ = Param.constrain_fixed.__doc__
+    fix = constrain_fixed
+    def constrain_negative(self, warning=True):
+        [param.constrain_negative(warning) for param in self.params]
+    constrain_negative.__doc__ = Param.constrain_negative.__doc__
+    def constrain_bounded(self, lower, upper, warning=True):
+        [param.constrain_bounded(lower, upper, warning) for param in self.params]
+    constrain_bounded.__doc__ = Param.constrain_bounded.__doc__
+    def unconstrain(self, *constraints):
+        [param.unconstrain(*constraints) for param in self.params]
+    unconstrain.__doc__ = Param.unconstrain.__doc__
+    def unconstrain_negative(self):
+        [param.unconstrain_negative() for param in self.params]
+    unconstrain_negative.__doc__ = Param.unconstrain_negative.__doc__
+    def unconstrain_positive(self):
+        [param.unconstrain_positive() for param in self.params]
+    unconstrain_positive.__doc__ = Param.unconstrain_positive.__doc__
+    def unconstrain_fixed(self):
+        [param.unconstrain_fixed() for param in self.params]
+    unconstrain_fixed.__doc__ = Param.unconstrain_fixed.__doc__
+    unfix = unconstrain_fixed
+    def unconstrain_bounded(self, lower, upper):
+        [param.unconstrain_bounded(lower, upper) for param in self.params]
+    unconstrain_bounded.__doc__ = Param.unconstrain_bounded.__doc__
+    def untie(self, *ties):
+        [param.untie(*ties) for param in self.params]
+    __lt__ = lambda self, val: self._vals()<val
+    __le__ = lambda self, val: self._vals()<=val
+    __eq__ = lambda self, val: self._vals()==val
+    __ne__ = lambda self, val: self._vals()!=val
+    __gt__ = lambda self, val: self._vals()>val
+    __ge__ = lambda self, val: self._vals()>=val
+    def __str__(self, *args, **kwargs):
+        def f(p):
+            ind = p._raveled_index()
+            return p._constraints_for(ind), p._ties_for(ind)
+        params = self.params
+        constr_matrices, ties_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)]
+        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)

GPy/core/parameterization/parameter_core.py

@@ -0,0 +1,155 @@
+# Copyright (c) 2012, GPy authors (see AUTHORS.txt).
+# Licensed under the BSD 3-clause license (see LICENSE.txt)
+
+from transformations import Logexp, NegativeLogexp, Logistic
+
+__updated__ = '2013-12-16'
+
+def adjust_name_for_printing(name):
+    if name is not None:
+        return name.replace(" ", "_").replace(".", "_").replace("-","").replace("+","").replace("!","").replace("*","").replace("/","")
+    return ''
+
+class Observable(object):
+    _observers_ = {}
+    def add_observer(self, observer, callble):
+        self._observers_[observer] = callble
+        callble(self)
+    def remove_observer(self, observer):
+        del self._observers_[observer]
+    def _notify_observers(self):
+        [callble(self) for callble in self._observers_.itervalues()]
+        
+class Pickleable(object):
+    def getstate(self):
+        """
+        Returns the state of this class in a memento pattern.
+        The state must be a list-like structure of all the fields
+        this class needs to run.
+
+        See python doc "pickling" (`__getstate__` and `__setstate__`) for details.
+        """
+        raise NotImplementedError, "To be able to use pickling you need to implement this method"
+    def setstate(self, state):
+        """
+        Set the state (memento pattern) of this class to the given state.
+        Usually this is just the counterpart to getstate, such that
+        an object is a copy of another when calling
+    
+            copy = <classname>.__new__(*args,**kw).setstate(<to_be_copied>.getstate())
+            
+        See python doc "pickling" (`__getstate__` and `__setstate__`) for details.
+        """
+        raise NotImplementedError, "To be able to use pickling you need to implement this method"
+
+#===============================================================================
+# Foundation framework for parameterized and param objects:
+#===============================================================================
+
+class Parentable(object):
+    def __init__(self, direct_parent=None, highest_parent=None, parent_index=None):
+        super(Parentable,self).__init__()        
+        self._direct_parent_ = direct_parent
+        self._parent_index_ = parent_index
+        self._highest_parent_ = highest_parent
+        
+    def has_parent(self):
+        return self._direct_parent_ is not None
+    
+class Nameable(Parentable):
+    _name = None
+    def __init__(self, name, direct_parent=None, highest_parent=None, parent_index=None):
+        super(Nameable,self).__init__(direct_parent, highest_parent, parent_index)
+        self._name = name or self.__class__.__name__
+
+    @property
+    def name(self):
+        return self._name
+    @name.setter
+    def name(self, name):
+        from_name = self.name
+        self._name = name        
+        if self.has_parent():
+            self._direct_parent_._name_changed(self, from_name)
+            
+class Constrainable(Nameable):
+    def __init__(self, name):
+        super(Constrainable,self).__init__(name)
+    #===========================================================================
+    # Constrain operations -> done
+    #===========================================================================
+    def constrain(self, transform, warning=True, update=True):
+        """
+        :param transform: the :py:class:`GPy.core.transformations.Transformation`
+                          to constrain the this parameter to.
+        :param warning: print a warning if re-constraining parameters.
+        
+        Constrain the parameter to the given
+        :py:class:`GPy.core.transformations.Transformation`.
+        """
+        if self.has_parent():
+            self._highest_parent_._add_constrain(self, transform, warning)
+            if update:
+                self._highest_parent_.parameters_changed()
+        else:
+            for p in self._parameters_:
+                self._add_constrain(p, transform, warning)
+            if update:
+                self.parameters_changed()
+    def constrain_positive(self, warning=True):
+        """
+        :param warning: print a warning if re-constraining parameters.
+        
+        Constrain this parameter to the default positive constraint.
+        """
+        self.constrain(Logexp(), warning)
+
+    def constrain_negative(self, warning=True):
+        """
+        :param warning: print a warning if re-constraining parameters.
+        
+        Constrain this parameter to the default negative constraint.
+        """
+        self.constrain(NegativeLogexp(), warning)
+
+    def constrain_bounded(self, lower, upper, warning=True):
+        """
+        :param lower, upper: the limits to bound this parameter to
+        :param warning: print a warning if re-constraining parameters.
+        
+        Constrain this parameter to lie within the given range.
+        """
+        self.constrain(Logistic(lower, upper), warning)
+
+    def unconstrain(self, *transforms):
+        """
+        :param transforms: The transformations to unconstrain from.
+        
+        remove all :py:class:`GPy.core.transformations.Transformation` 
+        transformats of this parameter object.
+        """
+        if self.has_parent():
+            self._highest_parent_._remove_constrain(self, *transforms)
+        else:
+            for p in self._parameters_:
+                self._remove_constrain(p, *transforms)
+
+    def unconstrain_positive(self):
+        """
+        Remove positive constraint of this parameter. 
+        """
+        self.unconstrain(Logexp())
+
+    def unconstrain_negative(self):
+        """
+        Remove negative constraint of this parameter. 
+        """
+        self.unconstrain(NegativeLogexp())
+
+    def unconstrain_bounded(self, lower, upper):
+        """
+        :param lower, upper: the limits to unbound this parameter from
+        
+        Remove (lower, upper) bounded constrain from this parameter/
+        """
+        self.unconstrain(Logistic(lower, upper))

GPy/core/parameterization/priors.py

@@ -0,0 +1,268 @@
+# Copyright (c) 2012, GPy authors (see AUTHORS.txt).
+# Licensed under the BSD 3-clause license (see LICENSE.txt)
+
+
+import numpy as np
+import pylab as pb
+from scipy.special import gammaln, digamma
+from ...util.linalg import pdinv
+from domains import _REAL, _POSITIVE
+import warnings
+import weakref
+
+class Prior:
+    domain = None
+    
+    def pdf(self, x):
+        return np.exp(self.lnpdf(x))
+
+    def plot(self):
+        rvs = self.rvs(1000)
+        pb.hist(rvs, 100, normed=True)
+        xmin, xmax = pb.xlim()
+        xx = np.linspace(xmin, xmax, 1000)
+        pb.plot(xx, self.pdf(xx), 'r', linewidth=2)
+
+
+class Gaussian(Prior):
+    """
+    Implementation of the univariate Gaussian probability function, coupled with random variables.
+
+    :param mu: mean
+    :param sigma: standard deviation
+
+    .. Note:: Bishop 2006 notation is used throughout the code
+
+    """
+    domain = _REAL
+    _instances = []
+    def __new__(cls, mu, sigma): # Singleton:
+        if cls._instances:
+            cls._instances[:] = [instance for instance in cls._instances if instance()]
+            for instance in cls._instances:
+                if instance().mu == mu and instance().sigma == sigma:
+                    return instance()
+        o = super(Prior, cls).__new__(cls, mu, sigma)
+        cls._instances.append(weakref.ref(o))
+        return cls._instances[-1]()
+    def __init__(self, mu, sigma):
+        self.mu = float(mu)
+        self.sigma = float(sigma)
+        self.sigma2 = np.square(self.sigma)
+        self.constant = -0.5 * np.log(2 * np.pi * self.sigma2)
+
+    def __str__(self):
+        return "N(" + str(np.round(self.mu)) + ', ' + str(np.round(self.sigma2)) + ')'
+
+    def lnpdf(self, x):
+        return self.constant - 0.5 * np.square(x - self.mu) / self.sigma2
+
+    def lnpdf_grad(self, x):
+        return -(x - self.mu) / self.sigma2
+
+    def rvs(self, n):
+        return np.random.randn(n) * self.sigma + self.mu
+
+
+class LogGaussian(Prior):
+    """
+    Implementation of the univariate *log*-Gaussian probability function, coupled with random variables.
+
+    :param mu: mean
+    :param sigma: standard deviation
+
+    .. Note:: Bishop 2006 notation is used throughout the code
+
+    """
+    domain = _POSITIVE
+    _instances = []
+    def __new__(cls, mu, sigma): # Singleton:
+        if cls._instances:
+            cls._instances[:] = [instance for instance in cls._instances if instance()]
+            for instance in cls._instances:
+                if instance().mu == mu and instance().sigma == sigma:
+                    return instance()
+        o = super(Prior, cls).__new__(cls, mu, sigma)
+        cls._instances.append(weakref.ref(o))
+        return cls._instances[-1]()
+    def __init__(self, mu, sigma):
+        self.mu = float(mu)
+        self.sigma = float(sigma)
+        self.sigma2 = np.square(self.sigma)
+        self.constant = -0.5 * np.log(2 * np.pi * self.sigma2)
+
+    def __str__(self):
+        return "lnN(" + str(np.round(self.mu)) + ', ' + str(np.round(self.sigma2)) + ')'
+
+    def lnpdf(self, x):
+        return self.constant - 0.5 * np.square(np.log(x) - self.mu) / self.sigma2 - np.log(x)
+
+    def lnpdf_grad(self, x):
+        return -((np.log(x) - self.mu) / self.sigma2 + 1.) / x
+
+    def rvs(self, n):
+        return np.exp(np.random.randn(n) * self.sigma + self.mu)
+
+
+class MultivariateGaussian:
+    """
+    Implementation of the multivariate Gaussian probability function, coupled with random variables.
+
+    :param mu: mean (N-dimensional array)
+    :param var: covariance matrix (NxN)
+
+    .. Note:: Bishop 2006 notation is used throughout the code
+
+    """
+    domain = _REAL
+    _instances = []
+    def __new__(cls, mu, var): # Singleton:
+        if cls._instances:
+            cls._instances[:] = [instance for instance in cls._instances if instance()]
+            for instance in cls._instances:
+                if np.all(instance().mu == mu) and np.all(instance().var == var):
+                    return instance()
+        o = super(Prior, cls).__new__(cls, mu, var)
+        cls._instances.append(weakref.ref(o))
+        return cls._instances[-1]()
+    def __init__(self, mu, var):
+        self.mu = np.array(mu).flatten()
+        self.var = np.array(var)
+        assert len(self.var.shape) == 2
+        assert self.var.shape[0] == self.var.shape[1]
+        assert self.var.shape[0] == self.mu.size
+        self.input_dim = self.mu.size
+        self.inv, self.hld = pdinv(self.var)
+        self.constant = -0.5 * self.input_dim * np.log(2 * np.pi) - self.hld
+
+    def summary(self):
+        raise NotImplementedError
+
+    def pdf(self, x):
+        return np.exp(self.lnpdf(x))
+
+    def lnpdf(self, x):
+        d = x - self.mu
+        return self.constant - 0.5 * np.sum(d * np.dot(d, self.inv), 1)
+
+    def lnpdf_grad(self, x):
+        d = x - self.mu
+        return -np.dot(self.inv, d)
+
+    def rvs(self, n):
+        return np.random.multivariate_normal(self.mu, self.var, n)
+
+    def plot(self):
+        if self.input_dim == 2:
+            rvs = self.rvs(200)
+            pb.plot(rvs[:, 0], rvs[:, 1], 'kx', mew=1.5)
+            xmin, xmax = pb.xlim()
+            ymin, ymax = pb.ylim()
+            xx, yy = np.mgrid[xmin:xmax:100j, ymin:ymax:100j]
+            xflat = np.vstack((xx.flatten(), yy.flatten())).T
+            zz = self.pdf(xflat).reshape(100, 100)
+            pb.contour(xx, yy, zz, linewidths=2)
+
+
+def gamma_from_EV(E, V):
+    warnings.warn("use Gamma.from_EV to create Gamma Prior", FutureWarning)
+    return Gamma.from_EV(E, V)
+
+
+class Gamma(Prior):
+    """
+    Implementation of the Gamma probability function, coupled with random variables.
+
+    :param a: shape parameter
+    :param b: rate parameter (warning: it's the *inverse* of the scale)
+
+    .. Note:: Bishop 2006 notation is used throughout the code
+
+    """
+    domain = _POSITIVE
+    _instances = []
+    def __new__(cls, a, b): # Singleton:
+        if cls._instances:
+            cls._instances[:] = [instance for instance in cls._instances if instance()]
+            for instance in cls._instances:
+                if instance().a == a and instance().b == b:
+                    return instance()
+        o = super(Prior, cls).__new__(cls, a, b)
+        cls._instances.append(weakref.ref(o))
+        return cls._instances[-1]()
+    def __init__(self, a, b):
+        self.a = float(a)
+        self.b = float(b)
+        self.constant = -gammaln(self.a) + a * np.log(b)
+
+    def __str__(self):
+        return "Ga(" + str(np.round(self.a)) + ', ' + str(np.round(self.b)) + ')'
+
+    def summary(self):
+        ret = {"E[x]": self.a / self.b, \
+            "E[ln x]": digamma(self.a) - np.log(self.b), \
+            "var[x]": self.a / self.b / self.b, \
+            "Entropy": gammaln(self.a) - (self.a - 1.) * digamma(self.a) - np.log(self.b) + self.a}
+        if self.a > 1:
+            ret['Mode'] = (self.a - 1.) / self.b
+        else:
+            ret['mode'] = np.nan
+        return ret
+
+    def lnpdf(self, x):
+        return self.constant + (self.a - 1) * np.log(x) - self.b * x
+
+    def lnpdf_grad(self, x):
+        return (self.a - 1.) / x - self.b
+
+    def rvs(self, n):
+        return np.random.gamma(scale=1. / self.b, shape=self.a, size=n)
+    @staticmethod
+    def from_EV(E, V):
+        """
+        Creates an instance of a Gamma Prior  by specifying the Expected value(s)
+        and Variance(s) of the distribution.
+    
+        :param E: expected value
+        :param V: variance
+        """
+        a = np.square(E) / V
+        b = E / V
+        return Gamma(a, b)
+
+class inverse_gamma(Prior):
+    """
+    Implementation of the inverse-Gamma probability function, coupled with random variables.
+
+    :param a: shape parameter
+    :param b: rate parameter (warning: it's the *inverse* of the scale)
+
+    .. Note:: Bishop 2006 notation is used throughout the code
+
+    """
+    domain = _POSITIVE
+    def __new__(cls, a, b): # Singleton:
+        if cls._instances:
+            cls._instances[:] = [instance for instance in cls._instances if instance()]
+            for instance in cls._instances:
+                if instance().a == a and instance().b == b:
+                    return instance()
+        o = super(Prior, cls).__new__(cls, a, b)
+        cls._instances.append(weakref.ref(o))
+        return cls._instances[-1]()
+    def __init__(self, a, b):
+        self.a = float(a)
+        self.b = float(b)
+        self.constant = -gammaln(self.a) + a * np.log(b)
+
+    def __str__(self):
+        return "iGa(" + str(np.round(self.a)) + ', ' + str(np.round(self.b)) + ')'
+
+    def lnpdf(self, x):
+        return self.constant - (self.a + 1) * np.log(x) - self.b / x
+
+    def lnpdf_grad(self, x):
+        return -(self.a + 1.) / x + self.b / x ** 2
+
+    def rvs(self, n):
+        return 1. / np.random.gamma(scale=1. / self.b, shape=self.a, size=n)

GPy/core/parameterization/transformations.py

@@ -0,0 +1,170 @@
+# Copyright (c) 2012, GPy authors (see AUTHORS.txt).
+# Licensed under the BSD 3-clause license (see LICENSE.txt)
+
+
+import numpy as np
+from domains import _POSITIVE,_NEGATIVE, _BOUNDED
+import sys 
+import weakref
+_lim_val = -np.log(sys.float_info.epsilon) 
+
+class Transformation(object):
+    domain = None
+    _instance = None
+    def __new__(cls, *args, **kwargs):
+        if not cls._instance or cls._instance.__class__ is not cls:
+            cls._instance = super(Transformation, cls).__new__(cls, *args, **kwargs)
+        return cls._instance
+    def f(self, x):
+        raise NotImplementedError
+    def finv(self, x):
+        raise NotImplementedError
+    def gradfactor(self, f):
+        """ df_dx evaluated at self.f(x)=f"""
+        raise NotImplementedError
+    def initialize(self, f):
+        """ produce a sensible initial value for f(x)"""
+        raise NotImplementedError
+    def __str__(self):
+        raise NotImplementedError
+
+class Logexp(Transformation):
+    domain = _POSITIVE
+    def f(self, x):
+        return np.where(x>_lim_val, x, np.log(1. + np.exp(x)))
+    def finv(self, f):
+        return np.where(f>_lim_val, f, np.log(np.exp(f) - 1.))
+    def gradfactor(self, f):
+        return np.where(f>_lim_val, 1., 1 - np.exp(-f))
+    def initialize(self, f):
+        if np.any(f < 0.):
+            print "Warning: changing parameters to satisfy constraints"
+        return np.abs(f)
+    def __str__(self):
+        return '+ve'
+
+class NegativeLogexp(Transformation):
+    domain = _NEGATIVE
+    logexp = Logexp()
+    def f(self, x):
+        return -self.logexp.f(x)  # np.log(1. + np.exp(x))
+    def finv(self, f):
+        return self.logexp.finv(-f)  # np.log(np.exp(-f) - 1.)
+    def gradfactor(self, f):
+        return -self.logexp.gradfactor(-f)
+    def initialize(self, f):
+        return -self.logexp.initialize(f)  # np.abs(f)
+    def __str__(self):
+        return '-ve'
+
+class LogexpClipped(Logexp):
+    max_bound = 1e100
+    min_bound = 1e-10
+    log_max_bound = np.log(max_bound)
+    log_min_bound = np.log(min_bound)
+    domain = _POSITIVE
+    _instances = []
+    def __new__(cls, lower=1e-6, *args, **kwargs):
+        if cls._instances:
+            cls._instances[:] = [instance for instance in cls._instances if instance()]
+            for instance in cls._instances:
+                if instance().lower == lower:
+                    return instance()
+        o = super(Transformation, cls).__new__(cls, lower, *args, **kwargs)
+        cls._instances.append(weakref.ref(o))
+        return cls._instances[-1]()
+    def __init__(self, lower=1e-6):
+        self.lower = lower
+    def f(self, x):
+        exp = np.exp(np.clip(x, self.log_min_bound, self.log_max_bound))
+        f = np.log(1. + exp)
+#         if np.isnan(f).any():
+#             import ipdb;ipdb.set_trace()
+        return np.clip(f, self.min_bound, self.max_bound)
+    def finv(self, f):
+        return np.log(np.exp(f - 1.))
+    def gradfactor(self, f):
+        ef = np.exp(f) # np.clip(f, self.min_bound, self.max_bound))
+        gf = (ef - 1.) / ef
+        return gf # np.where(f < self.lower, 0, gf)
+    def initialize(self, f):
+        if np.any(f < 0.):
+            print "Warning: changing parameters to satisfy constraints"
+        return np.abs(f)
+    def __str__(self):
+        return '+ve_c'
+
+    
+class Exponent(Transformation):
+    # TODO: can't allow this to go to zero, need to set a lower bound. Similar with negative Exponent below. See old MATLAB code.
+    domain = _POSITIVE
+    def f(self, x):
+        return np.where(x<_lim_val, np.where(x>-_lim_val, np.exp(x), np.exp(-_lim_val)), np.exp(_lim_val))
+    def finv(self, x):
+        return np.log(x)
+    def gradfactor(self, f):
+        return f
+    def initialize(self, f):
+        if np.any(f < 0.):
+            print "Warning: changing parameters to satisfy constraints"
+        return np.abs(f)
+    def __str__(self):
+        return '+ve'
+
+class NegativeExponent(Exponent):
+    domain = _NEGATIVE
+    def f(self, x):
+        return -Exponent.f(x)
+    def finv(self, f):
+        return Exponent.finv(-f)
+    def gradfactor(self, f):
+        return f
+    def initialize(self, f):
+        return -Exponent.initialize(f) #np.abs(f)
+    def __str__(self):
+        return '-ve'
+
+class Square(Transformation):
+    domain = _POSITIVE
+    def f(self, x):
+        return x ** 2
+    def finv(self, x):
+        return np.sqrt(x)
+    def gradfactor(self, f):
+        return 2 * np.sqrt(f)
+    def initialize(self, f):
+        return np.abs(f)
+    def __str__(self):
+        return '+sq'
+
+class Logistic(Transformation):
+    domain = _BOUNDED
+    _instances = []
+    def __new__(cls, lower=1e-6, upper=1e-6, *args, **kwargs):
+        if cls._instances:
+            cls._instances[:] = [instance for instance in cls._instances if instance()]
+            for instance in cls._instances:
+                if instance().lower == lower and instance().upper == upper:
+                    return instance()
+        o = super(Transformation, cls).__new__(cls, lower, upper, *args, **kwargs)
+        cls._instances.append(weakref.ref(o))
+        return cls._instances[-1]()
+    def __init__(self, lower, upper):
+        assert lower < upper
+        self.lower, self.upper = float(lower), float(upper)
+        self.difference = self.upper - self.lower
+    def f(self, x):
+        return self.lower + self.difference / (1. + np.exp(-x))
+    def finv(self, f):
+        return np.log(np.clip(f - self.lower, 1e-10, np.inf) / np.clip(self.upper - f, 1e-10, np.inf))
+    def gradfactor(self, f):
+        return (f - self.lower) * (self.upper - f) / self.difference
+    def initialize(self, f):
+        if np.any(np.logical_or(f < self.lower, f > self.upper)):
+            print "Warning: changing parameters to satisfy constraints"
+        return np.where(np.logical_or(f < self.lower, f > self.upper), self.f(f * 0.), f)
+    def __str__(self):
+        return '{},{}'.format(self.lower, self.upper)
+
+ 
+

GPy/core/parameterization/variational.py

@@ -0,0 +1,66 @@
+'''
+Created on 6 Nov 2013
+
+@author: maxz
+'''
+import numpy as np
+from parameterized import Parameterized
+from param import Param
+from ...util.misc import param_to_array
+
+class Normal(Parameterized):
+    '''
+    Normal distribution for variational approximations.
+    
+    holds the means and variances for a factorizing multivariate normal distribution
+    '''
+    def __init__(self, means, variances, name='latent space'):
+        Parameterized.__init__(self, name=name)
+        self.means = Param("mean", means)
+        self.variances = Param('variance', variances)
+        self.add_parameters(self.means, self.variances)
+
+    def plot(self, fignum=None, ax=None, colors=None):
+        """
+        Plot latent space X in 1D:
+
+            - if fig is given, create input_dim subplots in fig and plot in these
+            - if ax is given plot input_dim 1D latent space plots of X into each `axis`
+            - if neither fig nor ax is given create a figure with fignum and plot in there
+
+        colors:
+            colors of different latent space dimensions input_dim
+
+        """
+        import pylab
+        if ax is None:
+            fig = pylab.figure(num=fignum, figsize=(8, min(12, (2 * self.means.shape[1]))))
+        if colors is None:
+            colors = pylab.gca()._get_lines.color_cycle
+            pylab.clf()
+        else:
+            colors = iter(colors)
+        plots = []
+        means, variances = param_to_array(self.means, self.variances)
+        x = np.arange(means.shape[0])
+        for i in range(means.shape[1]):
+            if ax is None:
+                a = fig.add_subplot(means.shape[1], 1, i + 1)
+            elif isinstance(ax, (tuple, list)):
+                a = ax[i]
+            else:
+                raise ValueError("Need one ax per latent dimnesion input_dim")
+            a.plot(means, c='k', alpha=.3)
+            plots.extend(a.plot(x, means.T[i], c=colors.next(), label=r"$\mathbf{{X_{{{}}}}}$".format(i)))
+            a.fill_between(x,
+                            means.T[i] - 2 * np.sqrt(variances.T[i]),
+                            means.T[i] + 2 * np.sqrt(variances.T[i]),
+                            facecolor=plots[-1].get_color(),
+                            alpha=.3)
+            a.legend(borderaxespad=0.)
+            a.set_xlim(x.min(), x.max())
+            if i < means.shape[1] - 1:
+                a.set_xticklabels('')
+        pylab.draw()
+        fig.tight_layout(h_pad=.01) # , rect=(0, 0, 1, .95))
+        return fig