code: def Kdiag(self, 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
docs: Compute the diagonal of the covariance matrix associated to X.