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.

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