def computeEntropyAndEnthalpy(self, u_kn=None, 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
        ----------
        u_kn : float, NxK array
            The energies of the state that are being used.
        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
        -------
        result_vals : dictionary

        Keys in the result_vals dictionary:

        'Delta_f' : np.ndarray, float, shape=(K, K)
            results['Delta_f'] is the dimensionless free energy difference f_j - f_i
        'dDelta_f' : np.ndarray, float, shape=(K, K)
            uncertainty in results['Delta_f']
        'Delta_u' : np.ndarray, float, shape=(K, K)
            results['Delta_u'] is the reduced potential energy difference u_j - u_i
        'dDelta_u' : np.ndarray, float, shape=(K, K)
            uncertainty in results['Delta_u']
        'Delta_s' : np.ndarray, float, shape=(K, K)
            results['Delta_s'] is the reduced entropy difference S/k between states i and j (s_j - s_i)
        'dDelta_s' : np.ndarray, float, shape=(K, K)
            uncertainty in results['Delta_s']

        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)
        >>> results = mbar.computeEntropyAndEnthalpy()

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

        dims = len(np.shape(u_kn))
        if dims==3:
            u_kn = kln_to_kn(u_kn, N_k=self.N_k)

        if u_kn is None:
            u_kn = self.u_kn

        # Retrieve N and K for convenience.
        [K,N] = np.shape(u_kn)
        A_in = u_kn.copy()
        state_map = np.zeros([2,K],int)
        for k in range(K):
            state_map[0,k] = k
            state_map[1,k] = k

        inner_results = self.computeExpectationsInner(A_in, u_kn, state_map,
                                                      return_theta=True,
                                                      uncertainty_method=uncertainty_method,
                                                      warning_cutoff=warning_cutoff)

        # construct the covariance matrix of exp(ln c_Ua - ln c_a) - ln c_ca

        Theta = np.zeros([3*K,3*K],dtype=np.float64)
        Theta[0:2*K,0:2*K] = inner_results['Theta']
        Theta[2*K:3*K,:] = Theta[K:2*K,:]
        Theta[:,2*K:3*K] = Theta[:,K:2*K]
        diag = np.ones(3*K,dtype=np.float64)
        diag[0:K] = diag[K:2*K] = inner_results['observables']-inner_results['Amin']
        Adiag = np.matrix(np.zeros([3*K,3*K],dtype=np.float64))
        np.fill_diagonal(Adiag,diag)
        Theta = Adiag*Theta*Adiag

        # Compute reduced free energy difference.
        f_k = np.matrix(inner_results['f'])
        Delta_f_ij = np.array(f_k - f_k.transpose())
        # compute uncertainty matrix in free energies:
        covf = Theta[2*K:3*K,2*K:3*K]
        dDelta_f_ij = self._ErrorOfDifferences(covf,warning_cutoff=warning_cutoff)

        # Compute reduced enthalpy difference.
        u_k = np.matrix(inner_results['observables'])
        Delta_u_ij = np.array(u_k - u_k.transpose())
        # compute uncertainty matrix in energies:
        covu = Theta[0:K,0:K]+Theta[K:2*K,K:2*K]-Theta[0:K,K:2*K]-Theta[K:2*K,0:K]
        dDelta_u_ij = self._ErrorOfDifferences(covu,warning_cutoff=warning_cutoff)

        # Compute reduced entropy difference
        s_k = u_k - f_k
        Delta_s_ij = np.array(s_k - s_k.transpose())
        # compute uncertainty matrix in entropies
        #s_i = u_i - f_i
        #cov(s_i) =   cov(u_i - f_i)
        #         =   cov(exp(ln C_a - ln c_a) + ln c_a)
        #         =   cov(exp(ln C_a - ln c_a), exp(ln C_a - ln c_a)) + cov(ln c_a, ln c_a)
        #           + cov(exp(ln C_a - ln c_a), ln c_a) + cov(ln c_a, exp(ln C_a - ln c_a))
        #         = cov(u,u) + cov(f,f)
        #             + A cov(ln C_a - ln c_a, ln c_a) + A cov(ln c_a, ln C_a - ln c_a)
        #         = cov(u,u) + cov(f,f)
        #             + A cov(ln C_a, ln c_a) - A cov(ln c_a, ln c_a) + A cov(ln c_a, ln C_a) - A cov(ln c_a, ln c_a)
        #         = cov(u,u) + cov(f,f) + A cov(ln C_a,ln c_a) + A cov(ln c_a, ln C_a) - 2A cov(ln_ca,ln_ca)
        #
        covs = covu + covf + Theta[0:K,2*K:3*K] + Theta[2*K:3*K,0:K] - Theta[K:2*K,2*K:3*K] - Theta[2*K:3*K,K:2*K]
        # note: not clear that Theta[K:2*K,2*K:3*K] and Theta[K:2*K,2*K:3*K] are symmetric?
        dDelta_s_ij = self._ErrorOfDifferences(covs,warning_cutoff=warning_cutoff)

        result_vals = dict()
        result_vals['Delta_f'] = Delta_f_ij
        result_vals['dDelta_f'] = dDelta_f_ij
        result_vals['Delta_u'] = Delta_u_ij
        result_vals['dDelta_u'] = dDelta_u_ij
        result_vals['Delta_s'] = Delta_s_ij
        result_vals['dDelta_s'] = dDelta_s_ij

        return result_vals
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
        ----------
        u_kn : float, NxK array
            The energies of the state that are being used.
        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
        -------
        result_vals : dictionary

        Keys in the result_vals dictionary:

        'Delta_f' : np.ndarray, float, shape=(K, K)
            results['Delta_f'] is the dimensionless free energy difference f_j - f_i
        'dDelta_f' : np.ndarray, float, shape=(K, K)
            uncertainty in results['Delta_f']
        'Delta_u' : np.ndarray, float, shape=(K, K)
            results['Delta_u'] is the reduced potential energy difference u_j - u_i
        'dDelta_u' : np.ndarray, float, shape=(K, K)
            uncertainty in results['Delta_u']
        'Delta_s' : np.ndarray, float, shape=(K, K)
            results['Delta_s'] is the reduced entropy difference S/k between states i and j (s_j - s_i)
        'dDelta_s' : np.ndarray, float, shape=(K, K)
            uncertainty in results['Delta_s']

        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)
        >>> results = mbar.computeEntropyAndEnthalpy()