docs: def computeGeneralExpectations(self, A_in, u_ln, state_list, compute_uncertainty=True,
                                     uncertainty_method=None, warning_cutoff=1.0e-10, return_theta=False):
          """Compute the expectations of multiple observables of phase space functions on multiple states.
  
          Compute the expectations of multiple observables of phase
          space functions.  [A_0(x),A_1(x),...,A_i(x)] along with the
          covariances of their estimates at multiple states.
  
          We calculate all observables at all states that are specified by the choice of state list.
          Generally, it will be called in specific conditions.
  
          Parameters
          ----------
          A_in : np.ndarray, float, shape=(I, N)
              A_in[i,n] = A_i(x_n), the value of phase observable i for configuration n
          u_ln : np.ndarray, float, shape=(L, N)
              u_n[l,n] is the reduced potential of configuration n at state l
              if u_ln = None, we use self.u_kn
  
          state_list : np.ndarray, int, shape (NS,2), where NS is the
                       total number of states of interest.  it will be
                       of the form [[0,0],[1,1],[2,1]] which would
                       indicate we want to output the properties of
                       three observables total: the first property A[0]
                       at the 0th state, the 2nd property A[1] at the
                       1th state, and the 2nd property A[1] at the 2nd
                       state.  This allows us to tailor this to a number of different situations.
  
          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.
  
          Returns
          -------
  
          A_i : np.ndarray, float, shape = (I)
              A_i[i] is the estimate for the expectation of A_state_list[i](x) at the state specified by u_n[state_list[i],:]
  
          d2A_ik : np.ndarray, float, shape = (I, J)
              d2A_ij[i,j] is the COVARIANCE in the estimates of observables A_i and A_j (as determined by the state list)
              (* not the square root *)
  
          General cases this will be used for.
              single observable, multiple states (replacement for computeExpectations)
              multiple observables, single state (replacement for computeMultipleExpectations)
              diagonal cases of multiple states, single states.
  
          Examples
          --------
  
          update this example to be more general
          >>> 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_in = np.array([x_n,x_n**2,x_n**3])
          >>> u_n = u_kn[:2,:]
          >>> state_list = np.array([[0,0],[1,0],[2,0],[2,1]],int)
          >>> [A_i, d2A_ij] = mbar.computeGeneralExpectations(A_in, u_n, state_list)
  
          """
  
          # Retrieve N and K for convenience.
          S = len(state_list) # number of computed expectations examined
          K = self.K
  
          N = self.N  # N is total number of samples
  
          # make observables all positive
          A_list = np.unique(state_list[:,0])
          I = len(A_list) # number of observables used
          A_min = np.zeros([I], dtype=np.float64)
  
          for i in A_list:  # only need to zero the ones we will use.  May be some repetition here.
              A_min[i] = np.min(A_in[i, :]) #find the minimum
              A_in[i, :] = A_in[i,:] - (A_min[i] - 1)  #all values now positive so that we can work in logarithmic scale
  
          # Augment W_nk, N_k, and c_k for q_A(x) for the observables, with one
          # row for the specified state and I rows for the observable at that
          # state.
          # log weight matrix
          sizea = K + 2*S # augmented size
          Log_W_nk = np.zeros([N, sizea], np.float64) # log weight matrix
          N_k = np.zeros([sizea], np.int32)  # counts
          f_k = np.zeros([sizea], np.float64)  # free energies
  
          # Fill in first section 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
          f_k[0:K] = self.f_k
  
          # Compute row of W matrix for the extra states corresponding to u_ln according to the state list
          for s in range(S):
              l = state_list[s,1]
              la = K+s  #l augmented
              Log_W_nk[:, la] = self._computeUnnormalizedLogWeights(u_ln[l,:])
              f_k[la] = -_logsum(Log_W_nk[:, la])
              Log_W_nk[:, la] += f_k[l]
  
          # Compute the remaining rows/columns of W_nk and c_k for the
          # observables.
          for s in range(S):
              sa = K+S+s  # augmented s
              i = state_list[s,0]
              l = state_list[s,1]
              Log_W_nk[:, sa] = np.log(A_in[i, :]) + Log_W_nk[:, K+l]
              f_k[sa] = -_logsum(Log_W_nk[:, sa])
              Log_W_nk[:, sa] += f_k[sa]    # normalize this row
  
          # Compute estimates.
          A_i = np.zeros([S], np.float64)
          for s in range(S):
              A_i[s] = np.exp(-f_k[K + S + s])
  
          if compute_uncertainty or return_theta:
              # Compute augmented asymptotic covariance matrix.
              W_nk = np.exp(Log_W_nk)
              Theta_ij = self._computeAsymptoticCovarianceMatrix(
                  W_nk, N_k, method=uncertainty_method)
  
          if compute_uncertainty:
              # Compute estimates of statistical covariance
              # these variances will be the same whether or not we subtract a different constant from each A_i
              # todo: vectorize
              # compute the covariance component without doing the double loop
              d2A_ij = np.zeros([S, S], np.float64)
              for i in range(S):
                  si = K+S+i
                  li = K+state_list[i,1]
                  for j in range(S):
                      sj = K+S+j
                      lj = K+state_list[j,1]
                      d2A_ij[i, j] = A_i[i] * A_i[j] * (
                          Theta_ij[si, sj] - Theta_ij[si, li] - Theta_ij[lj, sj] + Theta_ij[li, lj])
  
          # Now that covariances are computed, add the constants back to A_i that
          # were required to enforce positivity
          for s in range(S):
              A_i[s] += (A_min[state_list[s,0]] - 1)
  
          # these values may be used outside the routine, so copy back.
          for i in A_list:
              A_in[i, :] = A_in[i,:] + (A_min[i] - 1)
  
          returns = []
          returns.append(A_i)
  
          if compute_uncertainty:
              returns.append(d2A_ij)
  
          if return_theta:
              returns.append(Theta_ij)
  
          # Return expectations and uncertainties.
          return returns