docs: def computeExpectations(self, A_n, u_kn=None, output='averages', state_dependent=False,
                              compute_uncertainty=True, uncertainty_method=None,
                              warning_cutoff=1.0e-10, return_theta=False):
          """Compute the expectation of an observable of a phase space function.
  
          Compute the expectation of an observable of a single phase space
          function A(x) at all states where potentials are generated.
  
          Parameters
          ----------
          A_n : np.ndarray, float
              A_n (N_max np float64 array) - A_n[n] = A(x_n)
  
          u_kn : np.ndarray
              u_kn (energies of state of interest length N)
              default is self.u_kn
  
          output : string, optional
              'averages' outputs expectations of observables and 'differences' outputs
              a matrix of differences in the observables.
  
          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)
  
          state_dependent: bool, whether the expectations are state-dependent.
  
          Returns
          -------
          result_vals : dictionary
  
          Possible keys in the result_vals dictionary:
  
          'mu' : 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'
          'sigma' : 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
              or None, if compute_uncertainty is False.
          'Theta' ((KxK np float64 array): Covariance matrix of log weights
  
          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
          >>> results = mbar.computeExpectations(A_n)
          >>> A_n = u_kn[0,:]
          >>> results = mbar.computeExpectations(A_n, output='differences')
          """
  
          dims = len(np.shape(A_n))
  
          if dims > 2:
              print("Warning: dim=3 for (state_dependent==True) matrices for observables and dim=2 for (state_dependent==False) observables are deprecated; we suggest you convert to NxK form instead of NxKxK form.")
  
          if not state_dependent:
              if dims==2:
                  A_n = kn_to_n(A_n, N_k=self.N_k)
                  if u_kn is not None:
                      if len(np.shape(u_kn)) == 3:
                          u_kn = kln_to_kn(u_kn, N_k=self.N_k)
                      elif len(np.shape(u_kn)) == 2:
                          u_kn = kn_to_n(u_kn, N_k=self.N_k)
          else:
              if dims==3:
                  A_n = kln_to_kn(A_n, N_k=self.N_k)
                  if u_kn is not None:
                      if len(np.shape(u_kn)) == 3:
                          u_kn = kln_to_kn(u_kn, N_k=self.N_k)
                      elif len(np.shape(u_kn)) == 2:
                          u_kn = kn_to_n(u_kn, N_k=self.N_k)
  
          if u_kn is None:
              u_kn = self.u_kn
  
          # Retrieve N and K for convenience.
          N = self.N
          ushape = np.shape(u_kn)
          if len(ushape) == 1:
              K = 1
          else:
              K = np.shape(u_kn)[0] # number of potentials provided.
  
          state_map = np.zeros([2,K],int)
          if state_dependent:
              for k in range(K):
                  # first property at the first state, 2nd property at the 2nd state
                  state_map[0,k] = k
                  state_map[1,k] = k
          else:
              # only one property, evaluate at K different states.
              for k in range(K):
                  state_map[0,k] = k
                  state_map[1,k] = 0
  
          inner_results = self.computeExpectationsInner(A_n,u_kn,state_map,
                                                        return_theta=compute_uncertainty,
                                                        uncertainty_method=uncertainty_method,
                                                        warning_cutoff=warning_cutoff)
  
          result_vals = dict()
          if compute_uncertainty or return_theta:
              # we want the theta matrix for the exponentials of the
              # observables, which means we need to make the
              # transformation.
              Adiag = np.zeros([2*K,2*K],dtype=np.float64)
              diag = np.ones(2*K,dtype=np.float64)
              diag[0:K] = diag[K:2*K] = inner_results['observables']-inner_results['Amin']
              np.fill_diagonal(Adiag,diag)
              Theta = Adiag*inner_results['Theta']*Adiag
              covA_ij = np.array(Theta[0:K,0:K]+Theta[K:2*K,K:2*K]-Theta[0:K,K:2*K]-Theta[K:2*K,0:K])
  
          if output == 'averages':
              result_vals['mu'] = inner_results['observables']
              if compute_uncertainty:
                  result_vals['sigma'] = np.sqrt(covA_ij[0:K,0:K].diagonal())
  
          if output == 'differences':
              A_im = np.matrix(inner_results['observables'])
              A_ij = A_im - A_im.transpose()
  
              result_vals['mu'] = np.array(A_ij)
              if compute_uncertainty:
                  result_vals['sigma'] = self._ErrorOfDifferences(covA_ij,warning_cutoff=warning_cutoff)
  
          if return_theta:
              result_vals['Theta'] = Theta
  
          return result_vals