diff --git a/GPy/kern/src/independent_outputs.py b/GPy/kern/src/independent_outputs.py
index 5fd22b15..de636c53 100644
--- a/GPy/kern/src/independent_outputs.py
+++ b/GPy/kern/src/independent_outputs.py
@@ -56,7 +56,10 @@ class IndependentOutputs(CombinationKernel):
             self.single_kern = False
             self.kern = kernels
         super(IndependentOutputs, self).__init__(kernels=kernels, extra_dims=[index_dim], name=name)
-        self.index_dim = index_dim
+        # The combination kernel ALLWAYS puts the extra dimension last.
+        # Thus, the index dimension of this kernel is always the last dimension
+        # after slicing. This is why the index_dim is just the last column:
+        self.index_dim = -1
 
     def K(self,X ,X2=None):
         slices = index_to_slices(X[:,self.index_dim])
@@ -103,13 +106,10 @@ class IndependentOutputs(CombinationKernel):
         target = np.zeros(X.shape)
         kerns = itertools.repeat(self.kern) if self.single_kern else self.kern
         if X2 is None:
-            # TODO: make use of index_to_slices
-            # FIXME: Broken as X is already sliced out
-            # print("Warning, gradients_X may not be working, I believe X has already been sliced out by the slicer!")
             values = np.unique(X[:,self.index_dim])
             slices = [X[:,self.index_dim]==i for i in values]
-            [target.__setitem__(s, kern.gradients_X(dL_dK[s,s],X[s],None))
-              for kern, s in zip(kerns, slices)]
+            for kern, s in zip(kerns, slices):
+                target[s] += kern.gradients_X(dL_dK[s, :][:, s],X[s], None)
             #slices = index_to_slices(X[:,self.index_dim])
             #[[np.add(target[s], kern.gradients_X(dL_dK[s,s], X[s]), out=target[s])
             #  for s in slices_i] for kern, slices_i in zip(kerns, slices)]
@@ -121,8 +121,8 @@ class IndependentOutputs(CombinationKernel):
             values = np.unique(X[:,self.index_dim])
             slices = [X[:,self.index_dim]==i for i in values]
             slices2 = [X2[:,self.index_dim]==i for i in values]
-            [target.__setitem__(s, kern.gradients_X(dL_dK[s, :][:, s2],X[s],X2[s2]))
-              for kern, s, s2 in zip(kerns, slices, slices2)]
+            for kern, s, s2 in zip(kerns, slices, slices2):
+                target[s] += kern.gradients_X(dL_dK[s, :][:, s2],X[s],X2[s2])
             # TODO: make work with index_to_slices
             #slices = index_to_slices(X[:,self.index_dim])
             #slices2 = index_to_slices(X2[:,self.index_dim])
@@ -133,7 +133,9 @@ class IndependentOutputs(CombinationKernel):
         slices = index_to_slices(X[:,self.index_dim])
         kerns = itertools.repeat(self.kern) if self.single_kern else self.kern
         target = np.zeros(X.shape)
-        [[target.__setitem__(s, kern.gradients_X_diag(dL_dKdiag[s],X[s])) for s in slices_i] for kern, slices_i in zip(kerns, slices)]
+        for kern, slices_i in zip(kerns, slices):
+            for s in slices_i:
+                target[s] += kern.gradients_X_diag(dL_dKdiag[s],X[s])
         return target
 
     def update_gradients_diag(self, dL_dKdiag, X):
diff --git a/GPy/testing/kernel_tests.py b/GPy/testing/kernel_tests.py
index 093ed3e5..eafde0ec 100644
--- a/GPy/testing/kernel_tests.py
+++ b/GPy/testing/kernel_tests.py
@@ -402,7 +402,7 @@ class KernelTestsNonContinuous(unittest.TestCase):
         np.testing.assert_array_equal(kern.active_dims, [-1,0,1,2])
         np.testing.assert_array_equal(kern._all_dims_active, [0,1,2,-1])
 
-    @skip('Gradients for independend outputs with different X do not work correctly')
+    #@skip('Gradients for independend outputs with different X do not work correctly')
     def testIndependendGradients(self):
         k = GPy.kern.RBF(self.D, active_dims=range(self.D))
         kern = GPy.kern.IndependentOutputs(k, -1, 'ind_single')
@@ -418,7 +418,7 @@ class KernelTestsNonContinuous(unittest.TestCase):
         np.testing.assert_array_equal(kern.active_dims, [-1,0,2])
         np.testing.assert_array_equal(kern._all_dims_active, [0,1,2,-1])
 
-    @skip('Gradients for independend outputs with different X do not work correctly')
+    #@skip('Gradients for independend outputs with different X do not work correctly')
     def test_Hierarchical_gradients(self):
         k = [GPy.kern.RBF(2, active_dims=[0,2], name='rbf1'), GPy.kern.RBF(2, active_dims=[0,2], name='rbf2')]
         kern = GPy.kern.IndependentOutputs(k, -1, name='ind_split')
@@ -431,6 +431,10 @@ class KernelTestsNonContinuous(unittest.TestCase):
         X2 = self.X2[self.X2[:,-1]!=2]
         self.assertTrue(check_kernel_gradient_functions(kern, X=X, X2=X2, verbose=verbose, fixed_X_dims=-1))
 
+    def test_Coregionalize(self):
+        kern = GPy.kern.Coregionalize(1, output_dim=3, active_dims=[-1])
+        self.assertTrue(check_kernel_gradient_functions(kern, X=self.X, X2=self.X2, verbose=verbose, fixed_X_dims=-1))
+
 @unittest.skipIf(not config.getboolean('cython', 'working'),"Cython modules have not been built on this machine")
 class Coregionalize_cython_test(unittest.TestCase):
     """