docs: def computeExpectations(self, A_n, output='averages', compute_uncertainty=True, uncertainty_method=None, warning_cutoff=1.0e-10, return_theta=False, useGeneral = False, state_dependent = False):
          """Compute the expectation of an observable of a phase space function.
  
          Compute the expectation of an observable of phase space
          function A(x) at all states where potentials are generated,
          including states for which no samples were drawn.
  
          We assume observables are not function of the state.  u is not
          an observable -- it changes depending on the state. u_k is an
          observable; the energy of state k does not depend on the
          state. To compute the estimators of the energy at all K
          states, use . . .
  
          Parameters
          ----------
          A_n : np.ndarray, float
              A_n (N_max np float64 array) - A_n[n] = A(x_n)
          output : string, optional
              Either output averages, and uncertainties, or output a matrix of differences, with uncertainties.
          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)
          return_theta : bool, optional
              Whether or not to return the theta matrix.  Can be useful for complicated differences.
          useGeneral: bool, whether to use the GeneralExpectations formalism = False,
          state_dependent: bool, whether the expectations are state-dependent.
  
          Returns
          -------
          A : np.ndarray, float
              if output is 'averages'
              A_i  (K np float64 array) -  A_i[i] is the estimate for the expectation of A(x) for state i.
              if output is 'differences'
          dA : np.ndarray, float
              dA_i  (K np float64 array) - dA_i[i] is uncertainty estimate (one standard deviation) for A_i[i]
              or
              dA_ij (K np float64 array) - dA_ij[i,j] is uncertainty estimate (one standard deviation) for the difference in A beteen i and j
  
          Notes
          -----
  
          The reported statistical uncertainty should, in the asymptotic limit,
          reflect one standard deviation for the normal distribution of the estimate.
          The true expectation should fall within the interval [-dA, +dA] centered on the estimate 68% of the time, and within
          the interval [-2 dA, +2 dA] centered on the estimate 95% of the time.
          This will break down in cases where the number of samples is not large enough to reach the asymptotic normal limit.
          This 'breakdown' can be exacerbated by the computation of observables like indicator functions for histograms that are sparsely populated.
  
          References
          ----------
  
          See Section IV of [1].
  
          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)
          >>> A_n = x_n
          >>> (A_ij, dA_ij) = mbar.computeExpectations(A_n)
          >>> A_n = u_kn[0,:]
          >>> (A_ij, dA_ij) = mbar.computeExpectations(A_n, output='differences')
          """
  
          dims = len(np.shape(A_n))
  
          # Retrieve N and K for convenience.
          N = self.N
          K = self.K
  
          if dims == 3:
              print("expecting dim=1 or dim=2")
              return None
  
          if (useGeneral):
  
              state_list = np.zeros([K,2],int)
              if (state_dependent):
                  for k in range(K):
                      state_list[k,0] = k
                      state_list[k,1] = k
                  A_in = A_n
              else:
                  A_in = np.zeros([1,N], dtype=np.float64)
                  if dims == 2:
                      A_n = kn_to_n(A_n, N_k=self.N_k)
                  A_in[0,:] = A_n
  
                  for k in range(K):
                      state_list[k,0] = 0
                      state_list[k,1] = k
  
              general_results = self.computeGeneralExpectations(A_in, self.u_kn, state_list,
                                                                compute_uncertainty=compute_uncertainty,
                                                                uncertainty_method=uncertainty_method,
                                                                warning_cutoff=warning_cutoff,
                                                                return_theta=return_theta)
              returns = []
  
              if output == 'averages':
                  # Return expectations and uncertainties.
                  returns.append(general_results[0])
  
                  if compute_uncertainty:
                      indices = np.eye(K,dtype=bool)
                      returns.append(np.sqrt(general_results[1][indices]))
  
              if output == 'differences':
                  A_im = np.matrix(general_results[0])
                  A_ij = A_im - A_im.transpose()
  
                  returns.append(np.array(A_ij))
                  if compute_uncertainty:
                      return np.sqrt(general_results[1])
  
              if return_theta:
                  returns.append(general_results[2])
  
          else:  # current style
  
              if dims == 2: #convert to 1xN shape
                  A_n = kn_to_n(A_n, N_k=self.N_k)
                  # Convert to np array.
                  A_n = np.array(A_n, np.float64)
  
              # Augment W_nk, N_k, and c_k for q_A(x) for the observable, with one
              # extra row/column for each state (Eq. 13 of [1]).
              # log of weight matrix
              Log_W_nk = np.zeros([N, K * 2], np.float64)
              N_k = np.zeros([K * 2], np.int32)  # counts
              # "free energies" of the new states
              f_k = np.zeros([K], 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
  
              # Make A_n all positive so we can operate logarithmically for
              # robustness
              A_i = np.zeros([K], np.float64)
              A_min = np.min(A_n)
              A_n = A_n - (A_min - 1)
  
              # Compute the remaining rows/columns of W_nk and the rows c_k for the
              # observables.
  
              for l in range(K):
  
                  # this works because all A_n are now positive;
                  Log_W_nk[:, K + l] = np.log(A_n) + self.Log_W_nk[:, l]
                  # we took the min at the beginning.
                  f_k[l] = -_logsum(Log_W_nk[:, K + l])
                  Log_W_nk[:, K + l] += f_k[l]              # normalize the row
                  A_i[l] = np.exp(-f_k[l])
  
              if compute_uncertainty or return_theta:
                  # Compute augmented asymptotic covariance matrix.
                  Theta_ij = self._computeAsymptoticCovarianceMatrix(
                      np.exp(Log_W_nk), N_k, method=uncertainty_method)
  
              returns = []
  
              if output == 'averages':
  
                  if compute_uncertainty:
                      # Compute uncertainties.
                      dA_i = np.zeros([K], np.float64)
                  # just the diagonals
                      for k in range(0, K):
                          dA_i[k] = np.abs(A_i[k]) * np.sqrt(
                              Theta_ij[K + k, K + k] + Theta_ij[k, k] - 2.0 * Theta_ij[k, K + k])
  
                  # add back minima now now that uncertainties are computed.
                      A_i += (A_min - 1)
  
                  # Return expectations and uncertainties.
                  returns.append(np.array(A_i))
  
                  if compute_uncertainty:
                      returns.append(np.array(dA_i))
  
              if output == 'differences':
                  # Return differences of expectations and uncertainties.
  
                  # compute expectation differences
                  A_im = np.matrix(A_i)
                  A_ij = A_im - A_im.transpose()
  
                  returns.append(np.array(A_ij))
  
                  # todo - vectorize the differences!  Faster and less likely to give errors.
  
                  if compute_uncertainty:
                      dA_ij = np.zeros([K, K], dtype=np.float64)
  
                      for i in range(0, K):
                          for j in range(0, K):
                              try:
                                  dA_ij[i, j] = np.sqrt(
                                      + A_i[i] * Theta_ij[i, i] * A_i[i]
                                        - A_i[i] * Theta_ij[i, j] * A_i[j]
                                        - A_i[i] * Theta_ij[i, K + i] * A_i[i]
                                        + A_i[i] * Theta_ij[i, K + j] * A_i[j]
                                        - A_i[j] * Theta_ij[j, i] * A_i[i]
                                        + A_i[j] * Theta_ij[j, j] * A_i[j]
                                        + A_i[j] * Theta_ij[j, K + i] * A_i[i]
                                        - A_i[j] * Theta_ij[j, K + j] * A_i[j]
                                        - A_i[i] * Theta_ij[K + i, i] * A_i[i]
                                        + A_i[i] * Theta_ij[K + i, j] * A_i[j]
                                        + A_i[i] * Theta_ij[K + i, K + i] * A_i[i]
                                        - A_i[i] * Theta_ij[K + i, K + j] * A_i[j]
                                        + A_i[j] * Theta_ij[K + j, i] * A_i[i]
                                        - A_i[j] * Theta_ij[K + j, j] * A_i[j]
                                        - A_i[j] * Theta_ij[K + j, K + i] * A_i[i]
                                        + A_i[j] * Theta_ij[K + j, K + j] * A_i[j]
                                        )
                              except:
                                  dA_ij[i, j] = 0.0
  
                      returns.append(dA_ij)
  
                  if return_theta:
                      returns.append(Theta_ij)
  
          return returns