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