code: def dchisq(psr,formbats=False,renormalize=True):
      
      
      if formbats:
          psr.formbats()
  
      res, err = psr.residuals(removemean=False)[psr.deleted == 0], psr.toaerrs[psr.deleted == 0]
  
      res -= numpy.sum(res/err**2) / numpy.sum(1/err**2)
      
      # bats already updated by residuals(); skip constant-phase column
      M = psr.designmatrix(updatebats=False,fixunits=True,fixsigns=True)[psr.deleted==0,1:]
  
      # renormalize design-matrix columns
      if renormalize:
          norm = numpy.sqrt(numpy.sum(M**2,axis=0))
          M /= norm
      else:
          norm = 1.0
  
      # compute chisq derivative, de-renormalize
      dr = -2 * numpy.dot(M.T,res / (1e-12 * err**2)) * norm
      
      return dr
docs: Return gradient of total chisq for the current timing solution,
      after removing noise-averaged mean residual, and ignoring deleted points.