<SYSTEM_TASK:>
Solve the following problem using Python, implementing the functions described below, one line at a time
<END_TASK>
<USER_TASK:>
Description:
 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