docs: def computePerturbedFreeEnergies(self, u_ln, compute_uncertainty=True, uncertainty_method=None, warning_cutoff=1.0e-10, return_theta=False):
          """Compute the free energies for a new set of states.
  
          Here, we desire the free energy differences among a set of new states, as well as the uncertainty estimates in these differences.
  
          Parameters
          ----------
          u_ln : np.ndarray, float, shape=(L, Nmax)
              u_ln[l,n] is the reduced potential energy of uncorrelated
              configuration n evaluated at new state k.  Can be completely indepednent of the original number of states.
          compute_uncertainty : bool, optional
              If False, the uncertainties will not be computed (default: True)
          uncertainty_method : string, optional
              Choice of method used to compute asymptotic covariance method, or None to use default
              See help for computeAsymptoticCovarianceMatrix() for more information on various methods. (default: None)
          warning_cutoff : float, optional
              Warn if squared-uncertainty is negative and larger in magnitude than this number (default: 1.0e-10)
  
          Returns
          -------
          Deltaf_ij : np.ndarray, float, shape=(L, L)
              Deltaf_ij[i,j] = f_j - f_i, the dimensionless free energy difference between new states i and j
          dDeltaf_ij : np.ndarray, float, shape=(L, L)
              dDeltaf_ij[i,j] is the estimated statistical uncertainty in Deltaf_ij[i,j]
  
          Examples
          --------
          >>> from pymbar import testsystems
          >>> (x_n, u_kn, N_k, s_n) = testsystems.HarmonicOscillatorsTestCase().sample(mode='u_kn')
          >>> mbar = MBAR(u_kn, N_k)
          >>> [Deltaf_ij, dDeltaf_ij] = mbar.computePerturbedFreeEnergies(u_kn)
          """
  
          # Convert to np matrix.
          u_ln = np.array(u_ln, dtype=np.float64)
  
          # Get the dimensions of the matrix of reduced potential energies, and convert if necessary
          if len(np.shape(u_ln)) == 3:
              u_ln = kln_to_kn(u_ln, N_k=self.N_k)
  
          [L, N] = u_ln.shape
  
          # Check dimensions.
          if (N < self.N):
              raise DataError("There seems to be too few samples in u_kn. You must evaluate at the new potential with all of the samples used originally.")
  
          # Retrieve N and K for convenience.
          N = self.N
          K = self.K
  
          # Augment W_nk, N_k, and c_k for the new states.
          W_nk = np.zeros([N, K + L], dtype=np.float64)  # weight matrix
          N_k = np.zeros([K + L], dtype=np.int32)  # counts
          f_k = np.zeros([K + L], dtype=np.float64)  # free energies
  
          # Fill in first half of matrix with existing q_k(x) from states.
          W_nk[:, 0:K] = np.exp(self.Log_W_nk)
          N_k[0:K] = self.N_k
          f_k[0:K] = self.f_k
  
          # Compute normalized weights.
          for l in range(0, L):
              # Compute unnormalized log weights.
              log_w_n = self._computeUnnormalizedLogWeights(u_ln[l, :])
              # Compute free energies
              f_k[K + l] = - _logsum(log_w_n)
              # Store normalized weights.  Keep in exponential not log form
              # because we will not store W_nk
              W_nk[:, K + l] = np.exp(log_w_n + f_k[K + l])
  
          if (compute_uncertainty or return_theta):
              # Compute augmented asymptotic covariance matrix.
              Theta_ij = self._computeAsymptoticCovarianceMatrix(
                  W_nk, N_k, method = uncertainty_method)
  
          # Compute matrix of free energy differences between states and
          # associated uncertainties.
          # makes matrix operations easier to recast
          f_k = np.matrix(f_k[K:K + L])
  
          Deltaf_ij = f_k - f_k.transpose()
  
          returns = []
          returns.append(Deltaf_ij)
  
          if (compute_uncertainty):
              diag = Theta_ij.diagonal()
              dii = diag[0, K:K + L]
              d2DeltaF = dii + dii.transpose() - 2 * Theta_ij[K:K + L, K:K + L]
  
              # check for any numbers below zero.
              if (np.any(d2DeltaF < 0.0)):
                  if(np.any(d2DeltaF) < warning_cutoff):
                      print("A squared uncertainty is negative.  d2DeltaF = %e" % d2DeltaF[(np.any(d2DeltaF) < warning_cutoff)])
                  else:
                      d2DeltaF[(np.any(d2DeltaF) < warning_cutoff)] = 0.0
  
              # take the square root of entries of the matrix
              dDeltaf_ij = np.sqrt(d2DeltaF)
  
              returns.append(dDeltaf_ij)
  
          if (return_theta):
              returns.append(Theta_ij)
  
          # Return matrix of free energy differences and uncertainties.
          return returns