code: def fgrad_y_psi(self, y, return_covar_chain=False):
          
          mpsi = self.psi
  
          w, s, r, d = self.fgrad_y(y, return_precalc=True)
          gradients = np.zeros((y.shape[0], y.shape[1], len(mpsi), 4))
          for i in range(len(mpsi)):
              a,b,c  = mpsi[i]
              gradients[:, :, i, 0] = (b * (1.0/np.cosh(s[i])) ** 2).T
              gradients[:, :, i, 1] = a * (d[i] - 2.0 * s[i] * r[i] * (1.0/np.cosh(s[i])) ** 2).T
              gradients[:, :, i, 2] = (-2.0 * a * (b ** 2) * r[i] * ((1.0 / np.cosh(s[i])) ** 2)).T
          gradients[:, :, 0, 3] = 1.0
  
          if return_covar_chain:
              covar_grad_chain = np.zeros((y.shape[0], y.shape[1], len(mpsi), 4))
              for i in range(len(mpsi)):
                  a,b,c = mpsi[i]
                  covar_grad_chain[:, :, i, 0] = (r[i]).T
                  covar_grad_chain[:, :, i, 1] = (a * (y + c) * ((1.0 / np.cosh(s[i])) ** 2).T)
                  covar_grad_chain[:, :, i, 2] = a * b * ((1.0 / np.cosh(s[i])) ** 2).T
              covar_grad_chain[:, :, 0, 3] = y
              return gradients, covar_grad_chain
  
          return gradients
docs: gradient of f w.r.t to y and psi
  
          :returns: NxIx4 tensor of partial derivatives