docs: def _computeAsymptoticCovarianceMatrix(self, W, N_k, method=None):
          """
          Compute estimate of the asymptotic covariance matrix.
  
          REQUIRED ARGUMENTS
            W (np.array of np.float of dimension [N,K]) - matrix of normalized weights (see Eq. 9 of [1]) - W[n,k] is the weight of snapshot n (n = 1..N) in state k
                                            Note that sum(W(:,k)) = 1 for any k = 1..K, and sum(N_k(:) .* W(n,:)) = 1 for any n.
            N_k (np.array of np.int32 of dimension [K]) - N_k[k] is the number of samples from state K
  
          RETURN VALUES
            Theta (KxK np float64 array) - asymptotic covariance matrix (see Eq. 8 of [1])
  
          OPTIONAL ARGUMENTS
            method (string) - if not None, specified method is used to compute asymptotic covariance method:
                              method must be one of ['generalized-inverse', 'svd', 'svd-ew', 'inverse', 'tan-HGH', 'tan', 'approximate']
                              If None is specified, 'svd-ew' is used.
  
          NOTES
  
          The computational costs of the various 'method' arguments varies:
  
            'generalized-inverse' currently requires computation of the pseudoinverse of an NxN matrix (where N is the total number of samples)
            'svd' computes the generalized inverse using the singular value decomposition -- this should be efficient yet accurate (faster)
            'svd-ev' is the same as 'svd', but uses the eigenvalue decomposition of W'W to bypass the need to perform an SVD (fastest)
            'inverse' only requires standard inversion of a KxK matrix (where K is the number of states), but requires all K states to be different
            'approximate' only requires multiplication of KxN and NxK matrices, but is an approximate underestimate of the uncertainty
            'tan' uses a simplified form that requires two pseudoinversions, but can be unstable
            'tan-HGH' makes weaker assumptions on 'tan' but can occasionally be unstable
  
          REFERENCE
            See Section II and Appendix D of [1].
  
          """
  
          # Set 'svd-ew' as default if uncertainty method specified as None.
          if method == None:
              method = 'svd-ew'
  
          # Get dimensions of weight matrix.
          [N, K] = W.shape
  
          # Check dimensions
          if(K != N_k.size):
              raise ParameterError(
                  'W must be NxK, where N_k is a K-dimensional array.')
          if(np.sum(N_k) != N):
              raise ParameterError('W must be NxK, where N = sum_k N_k.')
  
          # Check to make sure the weight matrix W is properly normalized.
          tolerance = 1.0e-4  # tolerance for checking equality of sums
  
          column_sums = np.sum(W, axis=0)
          badcolumns = (np.abs(column_sums - 1) > tolerance)
          if np.any(badcolumns):
              which_badcolumns = np.arange(K)[badcolumns]
              firstbad = which_badcolumns[0]
              raise ParameterError(
                  'Warning: Should have \sum_n W_nk = 1.  Actual column sum for state %d was %f. %d other columns have similar problems' %
                  (firstbad, column_sums[firstbad], np.sum(badcolumns)))
  
          row_sums = np.sum(W * N_k, axis=1)
          badrows = (np.abs(row_sums - 1) > tolerance)
          if np.any(badrows):
              which_badrows = np.arange(N)[badrows]
              firstbad = which_badrows[0]
              raise ParameterError(
                  'Warning: Should have \sum_k N_k W_nk = 1.  Actual row sum for sample %d was %f. %d other rows have similar problems' %
                  (firstbad, row_sums[firstbad], np.sum(badrows)))
  
          # Compute estimate of asymptotic covariance matrix using specified
          # method.
          if method == 'generalized-inverse':
              # Use generalized inverse (Eq. 8 of [1]) -- most general
              # Theta = W' (I - W N W')^+ W
  
              # Construct matrices
              # Diagonal N_k matrix.
              Ndiag = np.matrix(np.diag(N_k), dtype=np.float64)
              W = np.matrix(W, dtype=np.float64)
              I = np.identity(N, dtype=np.float64)
  
              # Compute covariance
              Theta = W.T * self._pseudoinverse(I - W * Ndiag * W.T) * W
  
          elif method == 'inverse':
              # Use standard inverse method (Eq. D8 of [1]) -- only applicable if all K states are different
              # Theta = [(W'W)^-1 - N + 1 1'/N]^-1
  
              # Construct matrices
              # Diagonal N_k matrix.
              Ndiag = np.matrix(np.diag(N_k), dtype=np.float64)
              W = np.matrix(W, dtype=np.float64)
              # I = np.identity(N, dtype=np.float64)  # Don't unnecessarily allocate O(N^2) memory.  See Issue #106
              # matrix of ones, times 1/N
              O = np.ones([K, K], dtype=np.float64) / float(N)
  
              # Make sure W is nonsingular.
              if (abs(linalg.det(W.T * W)) < tolerance):
                  print("Warning: W'W appears to be singular, yet 'inverse' method of uncertainty estimation requires W contain no duplicate states.")
  
              # Compute covariance
              Theta = ((W.T * W).I - Ndiag + O).I
  
          elif method == 'approximate':
              # Use fast approximate expression from Kong et al. -- this underestimates the true covariance, but may be a good approximation in some cases and requires no matrix inversions
              # Theta = P'P
  
              # Construct matrices
              W = np.matrix(W, dtype=np.float64)
  
              # Compute covariance
              Theta = W.T * W
  
          elif method == 'svd':
              # Use singular value decomposition based approach given in supplementary material to efficiently compute uncertainty
              # See Appendix D.1, Eq. D4 in [1].
  
              # Construct matrices
              Ndiag = np.matrix(np.diag(N_k), dtype=np.float64)
              W = np.matrix(W, dtype=np.float64)
              I = np.identity(K, dtype=np.float64)
  
              # Compute SVD of W
              [U, S, Vt] = linalg.svd(W)
              Sigma = np.matrix(np.diag(S))
              V = np.matrix(Vt).T
  
              # Compute covariance
              Theta = V * Sigma * self._pseudoinverse(
                  I - Sigma * V.T * Ndiag * V * Sigma) * Sigma * V.T
  
          elif method == 'svd-ew':
              # Use singular value decomposition based approach given in supplementary material to efficiently compute uncertainty
              # The eigenvalue decomposition of W'W is used to forego computing the SVD.
              # See Appendix D.1, Eqs. D4 and D5 of [1].
  
              # Construct matrices
              Ndiag = np.matrix(np.diag(N_k), dtype=np.float64)
              W = np.matrix(W, dtype=np.float64)
              I = np.identity(K, dtype=np.float64)
  
              # Compute singular values and right singular vectors of W without using SVD
              # Instead, we compute eigenvalues and eigenvectors of W'W.
              # Note W'W = (U S V')'(U S V') = V S' U' U S V' = V (S'S) V'
              [S2, V] = linalg.eigh(W.T * W)
              # Set any slightly negative eigenvalues to zero.
              S2[np.where(S2 < 0.0)] = 0.0
              # Form matrix of singular values Sigma, and V.
              Sigma = np.matrix(np.diag(np.sqrt(S2)))
              V = np.matrix(V)
  
              # Compute covariance
              Theta = V * Sigma * self._pseudoinverse(
                  I - Sigma * V.T * Ndiag * V * Sigma) * Sigma * V.T
  
          elif method == 'tan-HGH':
              # Use method suggested by Zhiqiang Tan without further simplification.
              # TODO: There may be a problem here -- double-check this.
  
              [N, K] = W.shape
  
              # Estimate O matrix from W'W.
              W = np.matrix(W, dtype=np.float64)
              O = W.T * W
  
              # Assemble the Lambda matrix.
              Lambda = np.matrix(np.diag(N_k), dtype=np.float64)
  
              # Identity matrix.
              I = np.matrix(np.eye(K), dtype=np.float64)
  
              # Compute H and G matrices.
              H = O * Lambda - I
              G = O - O * Lambda * O
  
              # Compute pseudoinverse of H
              Hinv = self._pseudoinverse(H)
  
              # Compute estimate of asymptotic covariance.
              Theta = Hinv * G * Hinv.T
  
          elif method == 'tan':
              # Use method suggested by Zhiqiang Tan.
  
              # Estimate O matrix from W'W.
              W = np.matrix(W, dtype=np.float64)
              O = W.T * W
  
              # Assemble the Lambda matrix.
              Lambda = np.matrix(np.diag(N_k), dtype=np.float64)
  
              # Compute covariance.
              Oinv = self._pseudoinverse(O)
              Theta = self._pseudoinverse(Oinv - Lambda)
  
          else:
              # Raise an exception.
              raise ParameterError('Method ' + method + ' unrecognized.')
  
          return Theta