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
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)