docs: 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