def computeEntropyAndEnthalpy(self, uncertainty_method=None, verbose=False, warning_cutoff=1.0e-10):
        """Decompose free energy differences into enthalpy and entropy differences.

        Compute the decomposition of the free energy difference between
        states 1 and N into reduced free energy differences, reduced potential
        (enthalpy) differences, and reduced entropy (S/k) differences.

        Parameters
        ----------
        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
        -------
        Delta_f_ij : np.ndarray, float, shape=(K, K)
            Delta_f_ij[i,j] is the dimensionless free energy difference f_j - f_i
        dDelta_f_ij : np.ndarray, float, shape=(K, K)
            uncertainty in Delta_f_ij
        Delta_u_ij : np.ndarray, float, shape=(K, K)
            Delta_u_ij[i,j] is the reduced potential energy difference u_j - u_i
        dDelta_u_ij : np.ndarray, float, shape=(K, K)
            uncertainty in Delta_f_ij
        Delta_s_ij : np.ndarray, float, shape=(K, K)
            Delta_s_ij[i,j] is the reduced entropy difference S/k between states i and j (s_j - s_i)
        dDelta_s_ij : np.ndarray, float, shape=(K, K)
            uncertainty in Delta_s_ij

        Notes
        -----
        This method is EXPERIMENTAL and should be used at your own risk.

        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)
        >>> [Delta_f_ij, dDelta_f_ij, Delta_u_ij, dDelta_u_ij, Delta_s_ij, dDelta_s_ij] = mbar.computeEntropyAndEnthalpy()

        """

        if verbose:
            print("Computing average energy and entropy by MBAR.")

        # Retrieve N and K for convenience.
        N = self.N
        K = self.K

        # Augment W_nk, N_k, and c_k for q_A(x) for the potential energies,
        # with one extra row/column for each state.
        # weight matrix
        Log_W_nk = np.zeros([N, K * 2], dtype=np.float64)
        N_k = np.zeros([K * 2], dtype=np.int32)  # counts
        # "free energies" of average states
        f_k = np.zeros(K, dtype=np.float64)

        # Fill in first half of matrix with existing q_k(x) from states.
        Log_W_nk[:, 0:K] = self.Log_W_nk
        N_k[0:K] = self.N_k

        # Compute the remaining rows/columns of W_nk and c_k for the potential
        # energy observable.

        u_min = self.u_kn.min()
        u_i = np.zeros([K], dtype=np.float64)
        for l in range(0, K):
            u_kn = self.u_kn[l, :] - (u_min-1)  # all positive now!  Subtracting off arbitrary constants doesn't affect results
            # since they are all differences.
            # Compute unnormalized weights.
            # A(x_n) exp[f_{k} - q_{k}(x_n)] / \sum_{k'=1}^K N_{k'} exp[f_{k'} - q_{k'}(x_n)]
            # harden for over/underflow with logarithms

            Log_W_nk[:, K + l] = np.log(u_kn) + self.Log_W_nk[:, l]

            f_k[l] = -_logsum(Log_W_nk[:, K + l])
            Log_W_nk[:, K + l] += f_k[l]              # normalize the row
            u_i[l] = np.exp(-f_k[l])

            # print "MBAR u_i[%d]: %10.5f,%10.5f" % (l,u_i[l]+u_min, u_i[l])

        # Compute augmented asymptotic covariance matrix.
        W_nk = np.exp(Log_W_nk)
        Theta_ij = self._computeAsymptoticCovarianceMatrix(
            W_nk, N_k, method=uncertainty_method)

        # Compute estimators and uncertainties.
        dDelta_f_ij = np.zeros([K, K], dtype=np.float64)
        dDelta_u_ij = np.zeros([K, K], dtype=np.float64)
        dDelta_s_ij = np.zeros([K, K], dtype=np.float64)

        # Compute reduced free energy difference.
        f_k = np.matrix(self.f_k)
        Delta_f_ij = f_k - f_k.transpose()

        # Compute reduced enthalpy difference.
        u_k = np.matrix(u_i)
        Delta_u_ij = u_k - u_k.transpose()

        # Compute reduced entropy difference
        s_k = u_k - f_k
        Delta_s_ij = s_k - s_k.transpose()

        # compute uncertainty matrix in free energies:
        # d2DeltaF = Theta_ij[i,i] + Theta_ij[j,j] - 2.0 * Theta_ij[i,j]

        diag = Theta_ij.diagonal()
        dii = diag[0:K, 0:K]
        d2DeltaF = dii + dii.transpose() - 2 * Theta_ij[0:K, 0:K]

        # check for any numbers below zero.
        if (np.any(d2DeltaF < 0.0)):
            if(np.any(d2DeltaF) < warning_cutoff):
                # Hmm.  Will this print correctly?
                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 the entries of matrix
        dDelta_f_ij = np.sqrt(d2DeltaF)
        # TODO -- vectorize this calculation for entropy and enthalpy!

        for i in range(0, K):
            for j in range(0, K):
                try:
                    dDelta_u_ij[i, j] = math.sqrt(
                        + u_i[i] * Theta_ij[i, i] * u_i[i] - u_i[i] * Theta_ij[i, j] * u_i[j] - u_i[
                            i] * Theta_ij[i, K + i] * u_i[i] + u_i[i] * Theta_ij[i, K + j] * u_i[j]
                        - u_i[j] * Theta_ij[j, i] * u_i[i] + u_i[j] * Theta_ij[j, j] * u_i[j] + u_i[
                            j] * Theta_ij[j, K + i] * u_i[i] - u_i[j] * Theta_ij[j, K + j] * u_i[j]
                        - u_i[i] * Theta_ij[K + i, i] * u_i[i] + u_i[i] * Theta_ij[K + i, j] * u_i[
                            j] + u_i[i] * Theta_ij[K + i, K + i] * u_i[i] - u_i[i] * Theta_ij[K + i, K + j] * u_i[j]
                        + u_i[j] * Theta_ij[K + j, i] * u_i[i] - u_i[j] * Theta_ij[K + j, j] * u_i[
                            j] - u_i[j] * Theta_ij[K + j, K + i] * u_i[i] + u_i[j] * Theta_ij[K + j, K + j] * u_i[j]
                    )
                except:
                    dDelta_u_ij[i, j] = 0.0

                # Compute reduced entropy difference.
                try:
                    dDelta_s_ij[i, j] = math.sqrt(
                        + (u_i[i] - 1) * Theta_ij[i, i] * (u_i[i] - 1) + (u_i[i] - 1) * Theta_ij[i, j] * (-u_i[j] + 1) + (
                            u_i[i] - 1) * Theta_ij[i, K + i] * (-u_i[i]) + (u_i[i] - 1) * Theta_ij[i, K + j] * u_i[j]
                        + (-u_i[j] + 1) * Theta_ij[j, i] * (u_i[i] - 1) + (-u_i[j] + 1) * Theta_ij[j, j] * (-u_i[j] + 1) +
                        (-u_i[j] + 1) * Theta_ij[j, K + i] * (-u_i[i]) +
                        (-u_i[j] + 1) * Theta_ij[j, K + j] * u_i[j]
                        + (-u_i[i]) * Theta_ij[K + i, i] * (u_i[i] - 1) + (-u_i[i]) * Theta_ij[K + i, j] * (-u_i[j] + 1) +
                        (-u_i[i]) * Theta_ij[K + i, K + i] * (-u_i[i]) +
                        (-u_i[i]) * Theta_ij[K + i, K + j] * u_i[j]
                        + u_i[j] * Theta_ij[K + j, i] * (u_i[i] - 1) + u_i[j] * Theta_ij[K + j, j] * (-u_i[j] + 1) + u_i[
                            j] * Theta_ij[K + j, K + i] * (-u_i[i]) + u_i[j] * Theta_ij[K + j, K + j] * u_i[j]
                    )
                except:
                    dDelta_s_ij[i, j] = 0.0

        # Return expectations and uncertainties.
        return (Delta_f_ij, dDelta_f_ij, Delta_u_ij, dDelta_u_ij, Delta_s_ij, dDelta_s_ij)
Decompose free energy differences into enthalpy and entropy differences.

        Compute the decomposition of the free energy difference between
        states 1 and N into reduced free energy differences, reduced potential
        (enthalpy) differences, and reduced entropy (S/k) differences.

        Parameters
        ----------
        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
        -------
        Delta_f_ij : np.ndarray, float, shape=(K, K)
            Delta_f_ij[i,j] is the dimensionless free energy difference f_j - f_i
        dDelta_f_ij : np.ndarray, float, shape=(K, K)
            uncertainty in Delta_f_ij
        Delta_u_ij : np.ndarray, float, shape=(K, K)
            Delta_u_ij[i,j] is the reduced potential energy difference u_j - u_i
        dDelta_u_ij : np.ndarray, float, shape=(K, K)
            uncertainty in Delta_f_ij
        Delta_s_ij : np.ndarray, float, shape=(K, K)
            Delta_s_ij[i,j] is the reduced entropy difference S/k between states i and j (s_j - s_i)
        dDelta_s_ij : np.ndarray, float, shape=(K, K)
            uncertainty in Delta_s_ij

        Notes
        -----
        This method is EXPERIMENTAL and should be used at your own risk.

        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)
        >>> [Delta_f_ij, dDelta_f_ij, Delta_u_ij, dDelta_u_ij, Delta_s_ij, dDelta_s_ij] = mbar.computeEntropyAndEnthalpy()