docs: def loglike(pulsar,efac=1.0,equad=None,jitter=None,Ared=None,gammared=None,marginalize=True,normalize=True,redcomponents=10,usedeleted=True):
      """Returns the Gaussian-process likelihood for 'pulsar'.
  
      The likelihood is evaluated at the current value of the pulsar parameters,
      as given by pulsar[parname].val.
  
      If efac, equad, and/or Ared are set, will compute the likelihood assuming
      the corresponding noise model. EFAC multiplies measurement noise;
      EQUAD adds in quadrature, and is given in us; red-noise is specified with
      the GW-like dimensionless amplitude Ared and exponent gamma, and is
      modeled with 'redcomponents' Fourier components.
  
      If marginalize=True (the default), loglike will marginalize over all the
      parameters in pulsar.fitpars, using an M-matrix formulation.
      """
  
      mask = Mask(pulsar,usedeleted)    
  
      err = 1.0e-6 * mask(pulsar.toaerrs)
      Cdiag = (efac*err)**2
  
      if equad:
          Cdiag = Cdiag + (1e-6*equad)**2 * N.ones(len(err))
  
      if Ared:
          redf, F = _setuprednoise(pulsar,redcomponents)
          F = mask(F)
          phi = Ared**2 * redf**(-gammared)
  
      if jitter:
          # quantize at 1 second; U plays the role of redF
          t, U = _quantize(86400.0 * mask(pulsar.toas()),1.0)
          phi_j = (1e-6*jitter)**2 * N.ones(U.shape[1])
  
          # stack the basis arrays if we're also doing red noise
          phi = N.hstack((phi,phi_j)) if Ared else phi_j
          F   = N.hstack((F,U))       if Ared else U
  
      if Ared or jitter:
          # Lentati formulation for correlated noise
          invphi = N.diag(1/phi)
          Ninv = N.diag(1/Cdiag)
          NinvF = dot(Ninv,F)          # could be accelerated
          X = invphi + dot(F.T,NinvF)  # invphi + FTNinvF
  
          Cinv = Ninv - dot(NinvF,N.linalg.inv(X),NinvF.T)
          logCdet = N.sum(N.log(Cdiag)) + N.sum(N.log(phi)) + N.linalg.slogdet(X)[1] # check
      else:
          # noise is all diagonal
          Cinv = N.diag(1/Cdiag)
          logCdet = N.sum(N.log(Cdiag))
  
      if marginalize:
          M = mask(pulsar.designmatrix())
          res = mask(N.array(pulsar.residuals(updatebats=False),'d'))
          
          CinvM = N.dot(Cinv,M)
          A = dot(M.T,CinvM)
  
          invA = N.linalg.inv(A)
          CinvMres = dot(res,CinvM)
  
          ret = -0.5 * dot(res,Cinv,res) + 0.5 * dot(CinvMres,invA,CinvMres.T) 
  
          if normalize:
              ret = ret - 0.5 * logCdet - 0.5 * N.linalg.slogdet(A)[1] - 0.5 * (M.shape[0] - M.shape[1]) * math.log(2.0*math.pi)
      else:
          res = mask(N.array(pulsar.residuals(),'d'))
  
          ret = -0.5 * dot(res,Cinv,res)
  
          if normalize:
              ret = ret - 0.5 * logCdet - 0.5 * len(res) * math.log(2.0*math.pi)
  
      return ret