diff --git a/GPy/core/parameterization/parameter_core.py b/GPy/core/parameterization/parameter_core.py index 16a1d23c..0357eb39 100644 --- a/GPy/core/parameterization/parameter_core.py +++ b/GPy/core/parameterization/parameter_core.py @@ -76,14 +76,14 @@ class Observable(object): def add_observer(self, observer, callble, priority=0): """ - Add an observer `observer` with the callback `callble` + Add an observer `observer` with the callback `callble` and priority `priority` to this observers list. """ self.observers.add(priority, observer, callble) def remove_observer(self, observer, callble=None): """ - Either (if callble is None) remove all callables, + Either (if callble is None) remove all callables, which were added alongside observer, or remove callable `callble` which was added alongside the observer `observer`. @@ -201,12 +201,12 @@ class Pickleable(object): #=========================================================================== def copy(self, memo=None, which=None): """ - Returns a (deep) copy of the current parameter handle. + Returns a (deep) copy of the current parameter handle. All connections to parents of the copy will be cut. - + :param dict memo: memo for deepcopy - :param Parameterized which: parameterized object which started the copy process [default: self] + :param Parameterized which: parameterized object which started the copy process [default: self] """ #raise NotImplementedError, "Copy is not yet implemented, TODO: Observable hierarchy" if memo is None: @@ -247,7 +247,7 @@ class Pickleable(object): if k not in ignore_list: dc[k] = v return dc - + def __setstate__(self, state): self.__dict__.update(state) from lists_and_dicts import ObserverList @@ -407,8 +407,7 @@ class Indexable(Nameable, Observable): if value is not None: self[:] = value - index = self._raveled_index() - #reconstrained = self.unconstrain() + index = self.unconstrain() index = self._add_to_index_operations(self.constraints, index, __fixed__, warning) self._highest_parent_._set_fixed(self, index) self.notify_observers(self, None if trigger_parent else -np.inf) @@ -641,24 +640,24 @@ class OptimizationHandlable(Indexable): #=========================================================================== # Optimizer copy - #=========================================================================== + #=========================================================================== @property def optimizer_array(self): """ Array for the optimizer to work on. This array always lives in the space for the optimizer. Thus, it is untransformed, going from Transformations. - + Setting this array, will make sure the transformed parameters for this model will be set accordingly. It has to be set with an array, retrieved from - this method, as e.g. fixing will resize the array. - + this method, as e.g. fixing will resize the array. + The optimizer should only interfere with this array, such that transofrmations are secured. """ if self.__dict__.get('_optimizer_copy_', None) is None or self.size != self._optimizer_copy_.size: self._optimizer_copy_ = np.empty(self.size) - + if not self._optimizer_copy_transformed: self._optimizer_copy_.flat = self.param_array.flat [np.put(self._optimizer_copy_, ind, c.finv(self.param_array[ind])) for c, ind in self.constraints.iteritems() if c != __fixed__] @@ -669,32 +668,32 @@ class OptimizationHandlable(Indexable): elif self._has_fixes(): return self._optimizer_copy_[self._fixes_] self._optimizer_copy_transformed = True - + return self._optimizer_copy_ - + @optimizer_array.setter def optimizer_array(self, p): """ - Make sure the optimizer copy does not get touched, thus, we only want to + Make sure the optimizer copy does not get touched, thus, we only want to set the values *inside* not the array itself. - + Also we want to update param_array in here. """ f = None if self.has_parent() and self.constraints[__fixed__].size != 0: f = np.ones(self.size).astype(bool) f[self.constraints[__fixed__]] = FIXED - elif self._has_fixes(): + elif self._has_fixes(): f = self._fixes_ if f is None: self.param_array.flat = p - [np.put(self.param_array, ind, c.f(self.param_array.flat[ind])) + [np.put(self.param_array, ind, c.f(self.param_array.flat[ind])) for c, ind in self.constraints.iteritems() if c != __fixed__] else: self.param_array.flat[f] = p - [np.put(self.param_array, ind[f[ind]], c.f(self.param_array.flat[ind[f[ind]]])) + [np.put(self.param_array, ind[f[ind]], c.f(self.param_array.flat[ind[f[ind]]])) for c, ind in self.constraints.iteritems() if c != __fixed__] - + self._optimizer_copy_transformed = False self._trigger_params_changed() @@ -710,13 +709,13 @@ class OptimizationHandlable(Indexable): # elif self._has_fixes(): # return p[self._fixes_] # return p -# +# def _set_params_transformed(self, p): raise DeprecationWarning, "_get|set_params{_optimizer_copy_transformed} is deprecated, use self.optimizer array insetad!" # """ # Set parameters p, but make sure they get transformed before setting. -# This means, the optimizer sees p, whereas the model sees transformed(p), +# This means, the optimizer sees p, whereas the model sees transformed(p), # such that, the parameters the model sees are in the right domain. # """ # if not(p is self.param_array): @@ -726,7 +725,7 @@ class OptimizationHandlable(Indexable): # self.param_array.flat[fixes] = p # elif self._has_fixes(): self.param_array.flat[self._fixes_] = p # else: self.param_array.flat = p -# [np.put(self.param_array, ind, c.f(self.param_array.flat[ind])) +# [np.put(self.param_array, ind, c.f(self.param_array.flat[ind])) # for c, ind in self.constraints.iteritems() if c != __fixed__] # self._trigger_params_changed() @@ -886,7 +885,7 @@ class Parameterizable(OptimizationHandlable): def traverse(self, visit, *args, **kwargs): """ - Traverse the hierarchy performing visit(self, *args, **kwargs) + Traverse the hierarchy performing visit(self, *args, **kwargs) at every node passed by downwards. This function includes self! See "visitor pattern" in literature. This is implemented in pre-order fashion. @@ -993,7 +992,7 @@ class Parameterizable(OptimizationHandlable): def _setup_observers(self): """ Setup the default observers - + 1: parameters_changed_notify 2: pass through to parent, if present """ diff --git a/GPy/kern/__init__.py b/GPy/kern/__init__.py index 27e056b4..2faf57c6 100644 --- a/GPy/kern/__init__.py +++ b/GPy/kern/__init__.py @@ -14,6 +14,7 @@ from _src.ODE_UYC import ODE_UYC from _src.ODE_st import ODE_st from _src.ODE_t import ODE_t from _src.poly import Poly +from _src.splitKern import SplitKern,DiffGenomeKern # TODO: put this in an init file somewhere #I'm commenting this out because the files were not added. JH. Remember to add the files before commiting diff --git a/GPy/kern/_src/independent_outputs.py b/GPy/kern/_src/independent_outputs.py index 64314197..21958267 100644 --- a/GPy/kern/_src/independent_outputs.py +++ b/GPy/kern/_src/independent_outputs.py @@ -20,9 +20,11 @@ def index_to_slices(index): returns >>> [[slice(0,2,None),slice(4,5,None)],[slice(2,4,None),slice(8,10,None)],[slice(5,8,None)]] """ + if len(index)==0: + return[] #contruct the return structure - ind = np.asarray(index,dtype=np.int64) + ind = np.asarray(index,dtype=np.int) ret = [[] for i in range(ind.max()+1)] #find the switchpoints diff --git a/GPy/kern/_src/splitKern.py b/GPy/kern/_src/splitKern.py new file mode 100644 index 00000000..27e4f76b --- /dev/null +++ b/GPy/kern/_src/splitKern.py @@ -0,0 +1,204 @@ +""" +A new kernel +""" + +import numpy as np +from kern import Kern,CombinationKernel +from .independent_outputs import index_to_slices +import itertools + +class DiffGenomeKern(Kern): + + def __init__(self, kernel, idx_p, Xp, index_dim=-1, name='DiffGenomeKern'): + self.idx_p = idx_p + self.index_dim=index_dim + self.kern = SplitKern(kernel,Xp, index_dim=index_dim) + super(DiffGenomeKern, self).__init__(input_dim=kernel.input_dim+1, active_dims=None, name=name) + self.add_parameter(self.kern) + + def K(self, X, X2=None): + assert X2==None + K = self.kern.K(X,X2) + + if self.idx_p<=0 or self.idx_p>X.shape[0]/2: + return K + + slices = index_to_slices(X[:,self.index_dim]) + idx_start = slices[1][0].start + idx_end = idx_start+self.idx_p + K_c = K[idx_start:idx_end,idx_start:idx_end].copy() + K[idx_start:idx_end,:] = K[:self.idx_p,:] + K[:,idx_start:idx_end] = K[:,:self.idx_p] + K[idx_start:idx_end,idx_start:idx_end] = K_c + + return K + + def Kdiag(self,X): + Kdiag = self.kern.Kdiag(X) + + if self.idx_p<=0 or self.idx_p>X.shape[0]/2: + return Kdiag + + slices = index_to_slices(X[:,self.index_dim]) + idx_start = slices[1][0].start + idx_end = idx_start+self.idx_p + Kdiag[idx_start:idx_end] = Kdiag[:self.idx_p] + + return Kdiag + + def update_gradients_full(self,dL_dK,X,X2=None): + assert X2==None + if self.idx_p<=0 or self.idx_p>X.shape[0]/2: + self.kern.update_gradients_full(dL_dK, X) + return + + slices = index_to_slices(X[:,self.index_dim]) + idx_start = slices[1][0].start + idx_end = idx_start+self.idx_p + + self.kern.update_gradients_full(dL_dK[idx_start:idx_end,:], X[:self.idx_p],X) + grad_p1 = self.kern.gradient.copy() + self.kern.update_gradients_full(dL_dK[:,idx_start:idx_end], X, X[:self.idx_p]) + grad_p2 = self.kern.gradient.copy() + self.kern.update_gradients_full(dL_dK[idx_start:idx_end,idx_start:idx_end], X[:self.idx_p],X[idx_start:idx_end]) + grad_p3 = self.kern.gradient.copy() + self.kern.update_gradients_full(dL_dK[idx_start:idx_end,idx_start:idx_end], X[idx_start:idx_end], X[:self.idx_p]) + grad_p4 = self.kern.gradient.copy() + + self.kern.update_gradients_full(dL_dK[idx_start:idx_end,:], X[idx_start:idx_end],X) + grad_n1 = self.kern.gradient.copy() + self.kern.update_gradients_full(dL_dK[:,idx_start:idx_end], X, X[idx_start:idx_end]) + grad_n2 = self.kern.gradient.copy() + self.kern.update_gradients_full(dL_dK[idx_start:idx_end,idx_start:idx_end], X[idx_start:idx_end], X[idx_start:idx_end]) + grad_n3 = self.kern.gradient.copy() + + self.kern.update_gradients_full(dL_dK, X) + self.kern.gradient += grad_p1+grad_p2-grad_p3-grad_p4-grad_n1-grad_n2+2*grad_n3 + + def update_gradients_diag(self, dL_dKdiag, X): + pass + +class SplitKern(CombinationKernel): + + def __init__(self, kernel, Xp, index_dim=-1, name='SplitKern'): + assert isinstance(index_dim, int), "The index dimension must be an integer!" + self.kern = kernel + self.kern_cross = SplitKern_cross(kernel,Xp) + super(SplitKern, self).__init__(kernels=[self.kern, self.kern_cross], extra_dims=[index_dim], name=name) + self.index_dim = index_dim + + def K(self,X ,X2=None): + slices = index_to_slices(X[:,self.index_dim]) + assert len(slices)<=2, 'The Split kernel only support two different indices' + if X2 is None: + target = np.zeros((X.shape[0], X.shape[0])) + # diagonal blocks + [[target.__setitem__((s,ss), self.kern.K(X[s,:], X[ss,:])) for s,ss in itertools.product(slices_i, slices_i)] for slices_i in slices] + if len(slices)>1: + # cross blocks + [target.__setitem__((s,ss), self.kern_cross.K(X[s,:], X[ss,:])) for s,ss in itertools.product(slices[0], slices[1])] + # cross blocks + [target.__setitem__((s,ss), self.kern_cross.K(X[s,:], X[ss,:])) for s,ss in itertools.product(slices[1], slices[0])] + else: + slices2 = index_to_slices(X2[:,self.index_dim]) + assert len(slices2)<=2, 'The Split kernel only support two different indices' + target = np.zeros((X.shape[0], X2.shape[0])) + # diagonal blocks + [[target.__setitem__((s,s2), self.kern.K(X[s,:],X2[s2,:])) for s,s2 in itertools.product(slices[i], slices2[i])] for i in xrange(min(len(slices),len(slices2)))] + if len(slices)>1: + [target.__setitem__((s,s2), self.kern_cross.K(X[s,:],X2[s2,:])) for s,s2 in itertools.product(slices[1], slices2[0])] + if len(slices2)>1: + [target.__setitem__((s,s2), self.kern_cross.K(X[s,:],X2[s2,:])) for s,s2 in itertools.product(slices[0], slices2[1])] + return target + + def Kdiag(self,X): + return self.kern.Kdiag(X) + + def update_gradients_full(self,dL_dK,X,X2=None): + slices = index_to_slices(X[:,self.index_dim]) + target = np.zeros(self.kern.size) + + def collate_grads(dL, X, X2, cross=False): + if cross: + self.kern_cross.update_gradients_full(dL,X,X2) + target[:] += self.kern_cross.kern.gradient + else: + self.kern.update_gradients_full(dL,X,X2) + target[:] += self.kern.gradient + + if X2 is None: + assert dL_dK.shape==(X.shape[0],X.shape[0]) + [[collate_grads(dL_dK[s,ss], X[s], X[ss]) for s,ss in itertools.product(slices_i, slices_i)] for slices_i in slices] + if len(slices)>1: + [collate_grads(dL_dK[s,ss], X[s], X[ss], True) for s,ss in itertools.product(slices[0], slices[1])] + [collate_grads(dL_dK[s,ss], X[s], X[ss], True) for s,ss in itertools.product(slices[1], slices[0])] + else: + assert dL_dK.shape==(X.shape[0],X2.shape[0]) + slices2 = index_to_slices(X2[:,self.index_dim]) + [[collate_grads(dL_dK[s,s2],X[s],X2[s2]) for s,s2 in itertools.product(slices[i], slices2[i])] for i in xrange(min(len(slices),len(slices2)))] + if len(slices)>1: + [collate_grads(dL_dK[s,s2], X[s], X2[s2], True) for s,s2 in itertools.product(slices[1], slices2[0])] + if len(slices2)>1: + [collate_grads(dL_dK[s,s2], X[s], X2[s2], True) for s,s2 in itertools.product(slices[0], slices2[1])] + self.kern.gradient = target + + def update_gradients_diag(self, dL_dKdiag, X): + self.kern.update_gradients_diag(self, dL_dKdiag, X) + +class SplitKern_cross(Kern): + + def __init__(self, kernel, Xp, name='SplitKern_cross'): + assert isinstance(kernel, Kern) + self.kern = kernel + if not isinstance(Xp,np.ndarray): + Xp = np.array([[Xp]]) + self.Xp = Xp + super(SplitKern_cross, self).__init__(input_dim=kernel.input_dim, active_dims=None, name=name) + + def K(self, X, X2=None): + if X2 is None: + return np.dot(self.kern.K(X,self.Xp),self.kern.K(self.Xp,X))/self.kern.K(self.Xp,self.Xp) + else: + return np.dot(self.kern.K(X,self.Xp),self.kern.K(self.Xp,X2))/self.kern.K(self.Xp,self.Xp) + + def Kdiag(self, X): + return np.inner(self.kern.K(X,self.Xp),self.kern.K(self.Xp,X).T)/self.kern.K(self.Xp,self.Xp) + + def update_gradients_full(self, dL_dK, X, X2=None): + if X2 is None: + X2 = X + + k1 = self.kern.K(X,self.Xp) + k2 = self.kern.K(self.Xp,X2) + k3 = self.kern.K(self.Xp,self.Xp) + dL_dk1 = np.einsum('ij,j->i',dL_dK,k2[0])/k3[0,0] + dL_dk2 = np.einsum('ij,i->j',dL_dK,k1[:,0])/k3[0,0] + dL_dk3 = np.einsum('ij,ij->',dL_dK,-np.dot(k1,k2)/(k3[0,0]*k3[0,0])) + + self.kern.update_gradients_full(dL_dk1[:,None],X,self.Xp) + grad = self.kern.gradient.copy() + self.kern.update_gradients_full(dL_dk2[None,:],self.Xp,X2) + grad += self.kern.gradient.copy() + self.kern.update_gradients_full(np.array([[dL_dk3]]),self.Xp,self.Xp) + grad += self.kern.gradient.copy() + + self.kern.gradient = grad + + def update_gradients_diag(self, dL_dKdiag, X): + k1 = self.kern.K(X,self.Xp) + k2 = self.kern.K(self.Xp,X) + k3 = self.kern.K(self.Xp,self.Xp) + dL_dk1 = dL_dKdiag*k2[0]/k3 + dL_dk2 = dL_dKdiag*k1[:,0]/k3 + dL_dk3 = -dL_dKdiag*(k1[:,0]*k2[0]).sum()/(k3*k3) + + self.kern.update_gradients_full(dL_dk1[:,None],X,self.Xp) + grad1 = self.kern.gradient.copy() + self.kern.update_gradients_full(dL_dk2[None,:],self.Xp,X) + grad2 = self.kern.gradient.copy() + self.kern.update_gradients_full(np.array([[dL_dk3]]),self.Xp,self.Xp) + grad3 = self.kern.gradient.copy() + + self.kern.gradient = grad1+grad2+grad3 + + diff --git a/GPy/util/linalg.py b/GPy/util/linalg.py index 661a2b47..b6b4a204 100644 --- a/GPy/util/linalg.py +++ b/GPy/util/linalg.py @@ -30,6 +30,7 @@ if config.getboolean('anaconda', 'installed') and config.getboolean('anaconda', dsyrk = mkl_rt.dsyrk dsyr = mkl_rt.dsyr _blas_available = True + print 'anaconda installed and mkl is loaded' except: _blas_available = False else: @@ -150,7 +151,7 @@ def dpotri(A, lower=1): assert lower==1, "scipy linalg behaviour is very weird. please use lower, fortran ordered arrays" A = force_F_ordered(A) - R, info = lapack.dpotri(A, lower=0) + R, info = lapack.dpotri(A, lower=0) #needs to be zero here, seems to be a scipy bug symmetrify(R) return R, info @@ -217,7 +218,7 @@ def pdinv(A, *args): L = jitchol(A, *args) logdet = 2.*np.sum(np.log(np.diag(L))) Li = dtrtri(L) - Ai, _ = lapack.dpotri(L) + Ai, _ = dpotri(L, lower=1) # Ai = np.tril(Ai) + np.tril(Ai,-1).T symmetrify(Ai)