def Kdiag(self, X):
        """Compute the diagonal of the covariance matrix associated to X."""
        vyt = self.variance_Yt
        vyx = self.variance_Yx

        lyt = 1./(2*self.lengthscale_Yt)
        lyx = 1./(2*self.lengthscale_Yx)

        a = self.a
        b = self.b
        c = self.c

        ## dk^2/dtdt'
        k1 = (2*lyt )*vyt*vyx
        ## dk^2/dx^2
        k2 = ( - 2*lyx )*vyt*vyx
        ## dk^4/dx^2dx'^2
        k3 = ( 4*3*lyx**2 )*vyt*vyx


        Kdiag = np.zeros(X.shape[0])
        slices = index_to_slices(X[:,-1])

        for i, ss1 in enumerate(slices):
            for s1 in ss1:
                if i==0:
                    Kdiag[s1]+= vyt*vyx
                elif i==1:
                    #i=1
                    Kdiag[s1]+= b**2*k1 - 2*a*c*k2 + a**2*k3 + c**2*vyt*vyx
                    #Kdiag[s1]+= Vu*Vy*(k1+k2+k3)
                else:
                    raise ValueError("invalid input/output index")

        return Kdiag