From 627cf9c86c7afbb01531fc4ea5879feca82abfee Mon Sep 17 00:00:00 2001
From: Zhenwen Dai <z.dai@shef.ac.uk>
Date: Thu, 29 May 2014 10:13:16 +0100
Subject: [PATCH 01/14] developing split kernel

---
 GPy/kern/__init__.py                 | 1 +
 GPy/kern/_src/independent_outputs.py | 2 +-
 2 files changed, 2 insertions(+), 1 deletion(-)

diff --git a/GPy/kern/__init__.py b/GPy/kern/__init__.py
index 27e056b4..85c6cf4d 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
 
 # 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..f07db692 100644
--- a/GPy/kern/_src/independent_outputs.py
+++ b/GPy/kern/_src/independent_outputs.py
@@ -22,7 +22,7 @@ def index_to_slices(index):
     """
 
     #contruct the return structure
-    ind = np.asarray(index,dtype=np.int64)
+    ind = np.asarray(index,dtype=np.int)
     ret = [[] for i in range(ind.max()+1)]
 
     #find the switchpoints

From 7a3d13894f14c55f9c7bf6744ba9e7d306f47739 Mon Sep 17 00:00:00 2001
From: James Hensman <james.hensman@gmail.com>
Date: Thu, 29 May 2014 09:55:49 +0100
Subject: [PATCH 02/14] reverting the fixing behaviour.

two reasons: 1) the new behaviour is confusing for new users. Either
something is fixed, or it's not. 2) the fixing didn't work! things that
should have been fixed were passed to the optimizer for optimization.

If we really want to save keystrokes, consider this:

m.foo.fix()
m.foo.unfix()
m.foo.constrain_positive()

is the same as

m.foo.fix()
m.foo.constrain_positive()

but the latter throws a warning.
---
 GPy/core/parameterization/parameter_core.py | 5 ++---
 1 file changed, 2 insertions(+), 3 deletions(-)

diff --git a/GPy/core/parameterization/parameter_core.py b/GPy/core/parameterization/parameter_core.py
index 9cedc46d..dab5c2d0 100644
--- a/GPy/core/parameterization/parameter_core.py
+++ b/GPy/core/parameterization/parameter_core.py
@@ -324,7 +324,6 @@ class Indexable(Nameable, Observable):
         self._default_constraint_ = default_constraint
         from index_operations import ParameterIndexOperations
         self.constraints = ParameterIndexOperations()
-        self._old_constraints = ParameterIndexOperations()
         self.priors = ParameterIndexOperations()
         if self._default_constraint_ is not None:
             self.constrain(self._default_constraint_)
@@ -388,8 +387,8 @@ class Indexable(Nameable, Observable):
         if value is not None:
             self[:] = value
 
-        index = self._raveled_index()
-        # reconstrained = self.unconstrain()
+        #index = self._raveled_index()
+        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)

From 988889247c9aa1f84564c33448a727b1fe165400 Mon Sep 17 00:00:00 2001
From: Zhenwen Dai <z.dai@shef.ac.uk>
Date: Thu, 29 May 2014 10:27:36 +0100
Subject: [PATCH 03/14] developing split kernel

---
 GPy/kern/_src/splitKern.py | 139 +++++++++++++++++++++++++++++++++++++
 1 file changed, 139 insertions(+)
 create mode 100644 GPy/kern/_src/splitKern.py

diff --git a/GPy/kern/_src/splitKern.py b/GPy/kern/_src/splitKern.py
new file mode 100644
index 00000000..eaa43468
--- /dev/null
+++ b/GPy/kern/_src/splitKern.py
@@ -0,0 +1,139 @@
+"""
+A new kernel
+"""
+
+import numpy as np
+from kern import Kern,CombinationKernel
+from .independent_outputs import index_to_slices
+import itertools
+
+class SplitKern(CombinationKernel):
+    """
+    A kernel which can represent several independent functions.  this kernel
+    'switches off' parts of the matrix where the output indexes are different.
+
+    The index of the functions is given by the last column in the input X the
+    rest of the columns of X are passed to the underlying kernel for
+    computation (in blocks).
+
+    :param kernels: either a kernel, or list of kernels to work with. If it is
+    a list of kernels the indices in the index_dim, index the kernels you gave!
+    """
+    def __init__(self, kernel, Xp, index_dim=-1, name='SplitKern'):
+        assert isinstance(index_dim, int), "The index dimension must be an integer!"
+        self.kern = kernel
+        self.kern_cross = SplitKern_cross(kernel,Xp)
+        super(SplitKern, self).__init__(kernels=[self.kern, self.kern_cross], extra_dims=[index_dim], name=name)
+        self.index_dim = index_dim
+
+    def K(self,X ,X2=None):
+        slices = index_to_slices(X[:,self.index_dim])
+        assert len(slices)<=2, 'The Split kernel only support two different indices'
+        if X2 is None:
+            target = np.zeros((X.shape[0], X.shape[0]))
+            # diagonal blocks
+            [[target.__setitem__((s,ss), self.kern.K(X[s,:], X[ss,:])) for s,ss in itertools.product(slices_i, slices_i)] for slices_i in slices]
+            if len(slices)>1:
+                # cross blocks
+                [target.__setitem__((s,ss), self.kern_cross.K(X[s,:], X[ss,:])) for s,ss in itertools.product(slices[0], slices[1])]
+                # cross blocks
+                [target.__setitem__((s,ss), self.kern_cross.K(X[s,:], X[ss,:])) for s,ss in itertools.product(slices[1], slices[0])]
+        else:
+            slices2 = index_to_slices(X2[:,self.index_dim])
+            assert len(slices2)<=2, 'The Split kernel only support two different indices'
+            target = np.zeros((X.shape[0], X2.shape[0]))
+            # diagonal blocks
+            [[target.__setitem__((s,s2), self.kern.K(X[s,:],X2[s2,:])) for s,s2 in itertools.product(slices[i], slices2[i])] for i in xrange(min(len(slices),len(slices)))]
+            if len(slices)>1:
+                [target.__setitem__((s,s2), self.kern.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.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.gradient
+            else:
+                self.kern.update_gradients_full(dL,X,X2)
+                target[:] += self.kern.gradient
+    
+        if X2 is None:
+            [[collate_grads(dL_dK[s,ss], X[s], X[ss]) for s,ss in itertools.product(slices_i, slices_i)] for slices_i in slices]
+            if len(slices)>1:
+                [collate_grads(dL_dK[s,ss], X[s], X[ss], True) for s,ss in itertools.product(slices[0], slices[1])]
+                [collate_grads(dL_dK[s,ss], X[s], X[ss], True) for s,ss in itertools.product(slices[1], slices[0])]
+        else:
+            slices2 = index_to_slices(X2[:,self.index_dim])
+            [[collate_grads(dL_dK[s,s2],X[s],X2[s2]) for s,s2 in itertools.product(slices[i], slices2[i])] for i in xrange(min(len(slices),len(slices)))]
+            if len(slices)>1:
+                [collate_grads(dL_dK[s,ss], X[s], X2[s2], True) for s,s2 in itertools.product(slices[1], slices2[0])]
+            if len(slices2)>1:
+                [collate_grads(dL_dK[s,ss], X[s], X2[s2], True) for s,s2 in itertools.product(slices[0], slices2[1])]
+        self.kern.gradient = target
+
+    def update_gradients_diag(self, dL_dKdiag, X):
+        self.kern.update_gradients_diag(self, dL_dKdiag, X)
+
+class SplitKern_cross(Kern):
+
+    def __init__(self, kernel, Xp, name='SplitKern_cross'):
+        assert isinstance(kernel, Kern)
+        self.kern = kernel
+        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
+        dL_dk2 = np.einsum('ij,i->j',dL_dK,k1[:,0])/k3
+        dL_dk3 = np.einsum('ij,ij->',dL_dK,-np.dot(k1,k2)/(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
+
+    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
+        
+

From 74b62ab3c3fe1758993fefe2caa4aec3c5a4cd2f Mon Sep 17 00:00:00 2001
From: Zhenwen Dai <z.dai@shef.ac.uk>
Date: Thu, 29 May 2014 10:29:18 +0100
Subject: [PATCH 04/14] merge the bug of fixing function

---
 GPy/core/parameterization/parameter_core.py | 3 +--
 1 file changed, 1 insertion(+), 2 deletions(-)

diff --git a/GPy/core/parameterization/parameter_core.py b/GPy/core/parameterization/parameter_core.py
index 16a1d23c..c1194bb0 100644
--- a/GPy/core/parameterization/parameter_core.py
+++ b/GPy/core/parameterization/parameter_core.py
@@ -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)

From f9291fe7da6d0bf8b4bacc154ba024348c84d0be Mon Sep 17 00:00:00 2001
From: Zhenwen Dai <z.dai@shef.ac.uk>
Date: Thu, 29 May 2014 11:21:19 +0100
Subject: [PATCH 05/14] [splitkern] buf fix

---
 GPy/kern/_src/splitKern.py | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/GPy/kern/_src/splitKern.py b/GPy/kern/_src/splitKern.py
index eaa43468..0aed7127 100644
--- a/GPy/kern/_src/splitKern.py
+++ b/GPy/kern/_src/splitKern.py
@@ -106,10 +106,10 @@ class SplitKern_cross(Kern):
         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
-        dL_dk2 = np.einsum('ij,i->j',dL_dK,k1[:,0])/k3
-        dL_dk3 = np.einsum('ij,ij->',dL_dK,-np.dot(k1,k2)/(k3*k3))
-        
+        dL_dk1 = np.einsum('ij,j->i',dL_dK,k2.flat)/k3.flat
+        dL_dk2 = np.einsum('ij,i->j',dL_dK,k1.flat)/k3.flat
+        dL_dk3 = np.einsum('ij,ij->',dL_dK,-np.dot(k1,k2)/(k3.flat*k3.flat))
+                
         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)

From c2568d2e9f81bbed537ef6a2e24076a97a55256c Mon Sep 17 00:00:00 2001
From: Zhenwen Dai <z.dai@shef.ac.uk>
Date: Thu, 29 May 2014 13:21:14 +0100
Subject: [PATCH 06/14] [splitkern] bug fix

---
 GPy/kern/_src/splitKern.py | 28 +++++++++++++++-------------
 1 file changed, 15 insertions(+), 13 deletions(-)

diff --git a/GPy/kern/_src/splitKern.py b/GPy/kern/_src/splitKern.py
index 0aed7127..dc7bef8b 100644
--- a/GPy/kern/_src/splitKern.py
+++ b/GPy/kern/_src/splitKern.py
@@ -6,6 +6,8 @@ import numpy as np
 from kern import Kern,CombinationKernel
 from .independent_outputs import index_to_slices
 import itertools
+from GPy.kern import Linear,RBF
+
 
 class SplitKern(CombinationKernel):
     """
@@ -45,9 +47,9 @@ class SplitKern(CombinationKernel):
             # diagonal blocks
             [[target.__setitem__((s,s2), self.kern.K(X[s,:],X2[s2,:])) for s,s2 in itertools.product(slices[i], slices2[i])] for i in xrange(min(len(slices),len(slices)))]
             if len(slices)>1:
-                [target.__setitem__((s,s2), self.kern.K(X[s,:],X2[s2,:])) for s,s2 in itertools.product(slices[1], slices2[0])]
+                [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.K(X[s,:],X2[s2,:])) for s,s2 in itertools.product(slices[0], 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):
@@ -60,7 +62,7 @@ class SplitKern(CombinationKernel):
         def collate_grads(dL, X, X2, cross=False):
             if cross:
                 self.kern_cross.update_gradients_full(dL,X,X2)
-                target[:] += self.kern_cross.gradient
+                target[:] += self.kern_cross.kern.gradient
             else:
                 self.kern.update_gradients_full(dL,X,X2)
                 target[:] += self.kern.gradient
@@ -102,22 +104,22 @@ class SplitKern_cross(Kern):
     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.flat)/k3.flat
-        dL_dk2 = np.einsum('ij,i->j',dL_dK,k1.flat)/k3.flat
-        dL_dk3 = np.einsum('ij,ij->',dL_dK,-np.dot(k1,k2)/(k3.flat*k3.flat))
-                
+        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)
-        grad1 = self.kern.gradient.copy()
-        self.kern.update_gradients_full(dL_dk2[None,:],self.Xp,X)
-        grad2 = self.kern.gradient.copy()
+        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)
-        grad3 = self.kern.gradient.copy()
+        grad += self.kern.gradient.copy()
         
-        self.kern.gradient = grad1+grad2+grad3
+        self.kern.gradient = grad
 
     def update_gradients_diag(self, dL_dKdiag, X):
         k1 = self.kern.K(X,self.Xp)

From 5f081526513889721614b2e4c4704a58d7f7d3b0 Mon Sep 17 00:00:00 2001
From: Zhenwen Dai <z.dai@shef.ac.uk>
Date: Thu, 29 May 2014 14:03:59 +0100
Subject: [PATCH 07/14] [splitkern] some additional implmentation

---
 GPy/kern/__init__.py       |  2 +-
 GPy/kern/_src/splitKern.py | 69 ++++++++++++++++++++++++++++++++------
 2 files changed, 59 insertions(+), 12 deletions(-)

diff --git a/GPy/kern/__init__.py b/GPy/kern/__init__.py
index 85c6cf4d..2faf57c6 100644
--- a/GPy/kern/__init__.py
+++ b/GPy/kern/__init__.py
@@ -14,7 +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
+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/splitKern.py b/GPy/kern/_src/splitKern.py
index dc7bef8b..6729b56f 100644
--- a/GPy/kern/_src/splitKern.py
+++ b/GPy/kern/_src/splitKern.py
@@ -6,21 +6,66 @@ import numpy as np
 from kern import Kern,CombinationKernel
 from .independent_outputs import index_to_slices
 import itertools
-from GPy.kern import Linear,RBF
 
+class DiffGenomeKern(Kern):
+
+    def __init__(self, kernel, idx_p, Xp, index_dim=-1, name='DiffGenomeKern'):
+        self.idx_p = idx_p
+        self.index_dim=index_dim
+        self.kern = SplitKern(kernel,Xp, index_dim=index_dim)
+        super(DiffGenomeKern, self).__init__(input_dim=kernel.input_dim+1, name=name)
+        self.add_parameter(self.kern)
+    
+    def K(self, X, X2=None):
+        assert X2==None
+        K = self.kern.K(X,X2)
+        
+        slices = index_to_slices(X[:,self.index_dim])
+        idx_start = slices[1][0]
+        idx_end = idx_start+self.idx_p
+        K[idx_start:idx_end,:] = K[:self.idx_p,:]
+        K[:,idx_start:idx_end] = K[:,self.idx_p]
+        
+        return K
+    
+    def Kdiag(self,X):
+        Kdiag = self.kern.Kdiag(X)
+
+        slices = index_to_slices(X[:,self.index_dim])
+        idx_start = slices[1][0]
+        idx_end = idx_start+self.idx_p
+        Kdiag[idx_start:idx_end] = Kdiag[:self.idx_p]
+        
+        return Kdiag
+    
+    def update_gradients_full(self,dL_dK,X,X2=None):
+        assert X2==None
+        slices = index_to_slices(X[:,self.index_dim])
+        idx_start = slices[1][0]
+        idx_end = idx_start+self.idx_p
+        
+        self.kern.update_gradients_full(dL_dK, X[:self.idx_p],X)
+        grad_p1 = self.kern.gradient.copy()
+        self.kern.update_gradients_full(dL_dK, X, X[:self.idx_p])
+        grad_p2 = self.kern.gradient.copy()
+        self.kern.update_gradients_full(dL_dK, X[:self.idx_p], X[:self.idx_p])
+        grad_p3 = self.kern.gradient.copy()
+
+        self.kern.update_gradients_full(dL_dK, X[idx_start:idx_end],X)
+        grad_n1 = self.kern.gradient.copy()
+        self.kern.update_gradients_full(dL_dK, X, X[idx_start:idx_end])
+        grad_n2 = self.kern.gradient.copy()
+        self.kern.update_gradients_full(dL_dK, X[idx_start:idx_end], X[idx_start:idx_end])
+        grad_n3 = self.kern.gradient.copy()
+
+        self.kern.update_gradients_full(dL_dK, X)
+        self.kern.gradient += grad_p1+grad_p2+grad_p3-grad_n1-grad_n2-grad_n3
+
+    def update_gradients_diag(self, dL_dKdiag, X):
+        pass
 
 class SplitKern(CombinationKernel):
-    """
-    A kernel which can represent several independent functions.  this kernel
-    'switches off' parts of the matrix where the output indexes are different.
 
-    The index of the functions is given by the last column in the input X the
-    rest of the columns of X are passed to the underlying kernel for
-    computation (in blocks).
-
-    :param kernels: either a kernel, or list of kernels to work with. If it is
-    a list of kernels the indices in the index_dim, index the kernels you gave!
-    """
     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
@@ -89,6 +134,8 @@ 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)
         

From ea1326c8c8db385ac9249371b696e87253a587f6 Mon Sep 17 00:00:00 2001
From: zhenwen <z.dai@sheffield.ac.uk>
Date: Thu, 29 May 2014 14:54:11 +0100
Subject: [PATCH 08/14] [splitkern] some more changes

---
 GPy/kern/_src/splitKern.py | 14 ++++++++------
 1 file changed, 8 insertions(+), 6 deletions(-)

diff --git a/GPy/kern/_src/splitKern.py b/GPy/kern/_src/splitKern.py
index 6729b56f..54a08914 100644
--- a/GPy/kern/_src/splitKern.py
+++ b/GPy/kern/_src/splitKern.py
@@ -13,7 +13,7 @@ class DiffGenomeKern(Kern):
         self.idx_p = idx_p
         self.index_dim=index_dim
         self.kern = SplitKern(kernel,Xp, index_dim=index_dim)
-        super(DiffGenomeKern, self).__init__(input_dim=kernel.input_dim+1, name=name)
+        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):
@@ -21,10 +21,12 @@ class DiffGenomeKern(Kern):
         K = self.kern.K(X,X2)
         
         slices = index_to_slices(X[:,self.index_dim])
-        idx_start = slices[1][0]
+        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] = K[:,:self.idx_p]
+        K[idx_start:idx_end,idx_start:idx_end] = K_c
         
         return K
     
@@ -32,7 +34,7 @@ class DiffGenomeKern(Kern):
         Kdiag = self.kern.Kdiag(X)
 
         slices = index_to_slices(X[:,self.index_dim])
-        idx_start = slices[1][0]
+        idx_start = slices[1][0].start
         idx_end = idx_start+self.idx_p
         Kdiag[idx_start:idx_end] = Kdiag[:self.idx_p]
         
@@ -41,7 +43,7 @@ class DiffGenomeKern(Kern):
     def update_gradients_full(self,dL_dK,X,X2=None):
         assert X2==None
         slices = index_to_slices(X[:,self.index_dim])
-        idx_start = slices[1][0]
+        idx_start = slices[1][0].start
         idx_end = idx_start+self.idx_p
         
         self.kern.update_gradients_full(dL_dK, X[:self.idx_p],X)
@@ -59,7 +61,7 @@ class DiffGenomeKern(Kern):
         grad_n3 = self.kern.gradient.copy()
 
         self.kern.update_gradients_full(dL_dK, X)
-        self.kern.gradient += grad_p1+grad_p2+grad_p3-grad_n1-grad_n2-grad_n3
+        self.kern.gradient += grad_p1+grad_p2-2*grad_p3-grad_n1-grad_n2+2*grad_n3
 
     def update_gradients_diag(self, dL_dKdiag, X):
         pass

From 4c538efb64192470449b6c7247d0b445f4404e5a Mon Sep 17 00:00:00 2001
From: Zhenwen Dai <z.dai@shef.ac.uk>
Date: Thu, 29 May 2014 17:10:14 +0100
Subject: [PATCH 09/14] DiffGenomeKern bug fix

---
 GPy/kern/_src/splitKern.py | 28 +++++++++++++++++-----------
 1 file changed, 17 insertions(+), 11 deletions(-)

diff --git a/GPy/kern/_src/splitKern.py b/GPy/kern/_src/splitKern.py
index 54a08914..624bae58 100644
--- a/GPy/kern/_src/splitKern.py
+++ b/GPy/kern/_src/splitKern.py
@@ -6,6 +6,7 @@ import numpy as np
 from kern import Kern,CombinationKernel
 from .independent_outputs import index_to_slices
 import itertools
+from .rbf import RBF
 
 class DiffGenomeKern(Kern):
 
@@ -13,6 +14,7 @@ class DiffGenomeKern(Kern):
         self.idx_p = idx_p
         self.index_dim=index_dim
         self.kern = SplitKern(kernel,Xp, index_dim=index_dim)
+#         self.kern = RBF(1)
         super(DiffGenomeKern, self).__init__(input_dim=kernel.input_dim+1, active_dims=None, name=name)
         self.add_parameter(self.kern)
     
@@ -46,22 +48,24 @@ class DiffGenomeKern(Kern):
         idx_start = slices[1][0].start
         idx_end = idx_start+self.idx_p
         
-        self.kern.update_gradients_full(dL_dK, X[:self.idx_p],X)
+        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, X, X[:self.idx_p])
+        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, X[:self.idx_p], X[:self.idx_p])
+        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, X[idx_start:idx_end],X)
+        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, X, X[idx_start:idx_end])
+        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, X[idx_start:idx_end], X[idx_start:idx_end])
+        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-2*grad_p3-grad_n1-grad_n2+2*grad_n3
+        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
@@ -92,7 +96,7 @@ class SplitKern(CombinationKernel):
             assert len(slices2)<=2, 'The Split kernel only support two different indices'
             target = np.zeros((X.shape[0], X2.shape[0]))
             # diagonal blocks
-            [[target.__setitem__((s,s2), self.kern.K(X[s,:],X2[s2,:])) for s,s2 in itertools.product(slices[i], slices2[i])] for i in xrange(min(len(slices),len(slices)))]
+            [[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:
@@ -115,17 +119,19 @@ class SplitKern(CombinationKernel):
                 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(slices)))]
+            [[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,ss], X[s], X2[s2], True) for s,s2 in itertools.product(slices[1], slices2[0])]
+                [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,ss], X[s], X2[s2], True) for s,s2 in itertools.product(slices[0], 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):

From 47ba2542c2e9c65d5cc797051dd1ed755422d0c0 Mon Sep 17 00:00:00 2001
From: Zhenwen Dai <z.dai@shef.ac.uk>
Date: Thu, 29 May 2014 17:11:26 +0100
Subject: [PATCH 10/14] minor changes

---
 GPy/kern/_src/splitKern.py | 2 --
 1 file changed, 2 deletions(-)

diff --git a/GPy/kern/_src/splitKern.py b/GPy/kern/_src/splitKern.py
index 624bae58..dfaf5108 100644
--- a/GPy/kern/_src/splitKern.py
+++ b/GPy/kern/_src/splitKern.py
@@ -6,7 +6,6 @@ import numpy as np
 from kern import Kern,CombinationKernel
 from .independent_outputs import index_to_slices
 import itertools
-from .rbf import RBF
 
 class DiffGenomeKern(Kern):
 
@@ -14,7 +13,6 @@ class DiffGenomeKern(Kern):
         self.idx_p = idx_p
         self.index_dim=index_dim
         self.kern = SplitKern(kernel,Xp, index_dim=index_dim)
-#         self.kern = RBF(1)
         super(DiffGenomeKern, self).__init__(input_dim=kernel.input_dim+1, active_dims=None, name=name)
         self.add_parameter(self.kern)
     

From ed74a817326065673bc799ee10901cb4d211df8c Mon Sep 17 00:00:00 2001
From: James Hensman <james.hensman@gmail.com>
Date: Fri, 30 May 2014 09:22:43 +0100
Subject: [PATCH 11/14] editied whitespace

---
 GPy/core/parameterization/parameter_core.py | 48 ++++++++++-----------
 1 file changed, 24 insertions(+), 24 deletions(-)

diff --git a/GPy/core/parameterization/parameter_core.py b/GPy/core/parameterization/parameter_core.py
index c1194bb0..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
@@ -640,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__]
@@ -668,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()
 
@@ -709,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):
@@ -725,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()
 
@@ -885,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.
@@ -992,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
         """

From 6e30f3c6c37f6e00a3d03597d8abb945af28f49e Mon Sep 17 00:00:00 2001
From: Zhenwen Dai <z.dai@shef.ac.uk>
Date: Fri, 30 May 2014 18:56:23 +0100
Subject: [PATCH 12/14] [splitkern] support idx_p==0

---
 GPy/kern/_src/independent_outputs.py |  2 ++
 GPy/kern/_src/splitKern.py           | 10 ++++++++++
 2 files changed, 12 insertions(+)

diff --git a/GPy/kern/_src/independent_outputs.py b/GPy/kern/_src/independent_outputs.py
index f07db692..21958267 100644
--- a/GPy/kern/_src/independent_outputs.py
+++ b/GPy/kern/_src/independent_outputs.py
@@ -20,6 +20,8 @@ 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.int)
diff --git a/GPy/kern/_src/splitKern.py b/GPy/kern/_src/splitKern.py
index dfaf5108..27e4f76b 100644
--- a/GPy/kern/_src/splitKern.py
+++ b/GPy/kern/_src/splitKern.py
@@ -20,6 +20,9 @@ class DiffGenomeKern(Kern):
         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
@@ -33,6 +36,9 @@ class DiffGenomeKern(Kern):
     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
@@ -42,6 +48,10 @@ class DiffGenomeKern(Kern):
     
     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

From 38ed60385aa0c9a2c1b54150e6efcf0fd5350620 Mon Sep 17 00:00:00 2001
From: Max Zwiessele <t-mazwie@microsoft.com>
Date: Tue, 3 Jun 2014 15:42:32 -0700
Subject: [PATCH 13/14] linalg had lowers missing for windows libraries to work
 correctly

---
 GPy/util/linalg.py | 5 +++--
 1 file changed, 3 insertions(+), 2 deletions(-)

diff --git a/GPy/util/linalg.py b/GPy/util/linalg.py
index 661a2b47..6a38a7c3 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=1)
     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, _ = lapack.dpotri(L, lower=1)
     # Ai = np.tril(Ai) + np.tril(Ai,-1).T
     symmetrify(Ai)
 

From d20222a90bac8181f39fac50f851d048e91c5bac Mon Sep 17 00:00:00 2001
From: James Hensman <james.hensman@gmail.com>
Date: Wed, 4 Jun 2014 15:43:03 +0100
Subject: [PATCH 14/14] reverting Maxs linalg changes

---
 GPy/util/linalg.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/GPy/util/linalg.py b/GPy/util/linalg.py
index 6a38a7c3..b6b4a204 100644
--- a/GPy/util/linalg.py
+++ b/GPy/util/linalg.py
@@ -151,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=1)
+    R, info = lapack.dpotri(A, lower=0) #needs to be zero here, seems to be a scipy bug
     symmetrify(R)
     return R, info
 
@@ -218,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, lower=1)
+    Ai, _ = dpotri(L, lower=1)
     # Ai = np.tril(Ai) + np.tril(Ai,-1).T
     symmetrify(Ai)