docs: def _calc_covariance(r, pmut, tol=1e-14):
      """Calculate the covariance matrix of the fitted parameters
  
  Parameters:
  r    - n-by-n matrix, the full upper triangle of R
  pmut - n-vector, defines the permutation of R
  tol  - scalar, relative column scale for determining rank
         deficiency. Default 1e-14.
  
  Returns:
  cov  - n-by-n matrix, the covariance matrix C
  
  Given an n-by-n matrix A, the corresponding covariance matrix
  is
  
    C = inverse(A^T A)
  
  This routine is given information relating to the pivoted transposed
  QR factorization of A, which is defined by matrices such that
  
   A P = R Q
  
  where P is a permutation matrix, Q has orthogonal rows, and R is a
  lower triangular matrix with diagonal elements of nonincreasing
  magnitude. In particular we take the full lower triangle of R ('r')
  and a vector describing P ('pmut'). The covariance matrix is then
  
   C = P inverse(R^T R) P^T
  
  If A is nearly rank-deficient, it may be desirable to compute the
  covariance matrix corresponding to the linearly-independent columns of
  A. We use a tolerance, 'tol', to define the numerical rank of A. If j
  is the largest integer such that |R[j,j]| > tol*|R[0,0]|, then we
  compute the covariance matrix for the first j columns of R. For k > j,
  the corresponding covariance entries (pmut[k]) are set to zero.
  """
      # This routine could save an allocation by operating on r in-place,
      # which might be worthwhile for large n, and is what the original
      # Fortran does.
  
      n = r.shape[1]
      assert r.shape[0] >= n
      r = r.copy()
  
      # Form the inverse of R in the full lower triangle of R.
  
      jrank = -1
      abstol = tol * abs(r[0,0])
  
      for i in range(n):
          if abs(r[i,i]) <= abstol:
              break
  
          r[i,i] **= -1
  
          for j in range(i):
              temp = r[i,i] * r[i,j]
              r[i,j] = 0.
              r[i,:j+1] -= temp * r[j,:j+1]
  
          jrank = i
  
      # Form the full lower triangle of the inverse(R^T R) in the full
      # lower triangle of R.
  
      for i in range(jrank + 1):
          for j in range(i):
              r[j,:j+1] += r[i,j] * r[i,:j+1]
          r[i,:i+1] *= r[i,i]
  
      # Form the full upper triangle of the covariance matrix in the
      # strict upper triangle of R and in wa.
  
      wa = np.empty(n)
      wa.fill(r[0,0])
  
      for i in range(n):
          pi = pmut[i]
          sing = i > jrank
  
          for j in range(i + 1):
              if sing:
                  r[i,j] = 0.
  
              pj = pmut[j]
              if pj > pi:
                  r[pi,pj] = r[i,j]
              elif pj < pi:
                  r[pj,pi] = r[i,j]
  
          wa[pi] = r[i,i]
  
      # Symmetrize.
  
      for i in range(n):
          r[i,:i+1] = r[:i+1,i]
          r[i,i] = wa[i]
  
      return r