Merge branch 'devel' of https://github.com/SheffieldML/GPy into devel

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()
+        index = self.unconstrain()
-        #reconstrained = 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
         """

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

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

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
+        
+

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)