docs: def computePMF(self, u_n, bin_n, nbins, uncertainties='from-lowest', pmf_reference=None):
          """
          Compute the free energy of occupying a number of bins.
  
          This implementation computes the expectation of an indicator-function observable for each bin.
  
          Parameters
          ----------
          u_n : np.ndarray, float, shape=(N)
              u_n[n] is the reduced potential energy of snapshot n of state k for which the PMF is to be computed.
          bin_n : np.ndarray, float, shape=(N)
              bin_n[n] is the bin index of snapshot n of state k.  bin_n can assume a value in range(0,nbins)
          nbins : int
              The number of bins
          uncertainties : string, optional
              Method for reporting uncertainties (default: 'from-lowest')
  
              * 'from-lowest' - the uncertainties in the free energy difference with lowest point on PMF are reported
              * 'from-specified' - same as from lowest, but from a user specified point
              * 'from-normalization' - the normalization \sum_i p_i = 1 is used to determine uncertainties spread out through the PMF
              * 'all-differences' - the nbins x nbins matrix df_ij of uncertainties in free energy differences is returned instead of df_i
  
          pmf_reference : int, optional
              the reference state that is zeroed when uncertainty = 'from-specified'
  
          Returns
          -------
          result_vals : dictionary
  
          Possible keys in the result_vals dictionary:
  
          'f_i' : np.ndarray, float, shape=(K)
              result_vals['f_i'][i] is the dimensionless free energy of state i, relative to the state of lowest free energy
          'df_i' : np.ndarray, float, shape=(K)
              result_vals['df_i'][i] is the uncertainty in the difference of f_i with respect to the state of lowest free energy
  
          Notes
          -----
          * All bins must have some samples in them from at least one of the states -- this will not work if bin_n.sum(0) == 0. Empty bins should be removed before calling computePMF().
          * This method works by computing the free energy of localizing the system to each bin for the given potential by aggregating the log weights for the given potential.
          * To estimate uncertainties, the NxK weight matrix W_nk is augmented to be Nx(K+nbins) in order to accomodate the normalized weights of states where
          * the potential is given by u_kn within each bin and infinite potential outside the bin.  The uncertainties with respect to the bin of lowest free energy are then computed in the standard way.
  
          Examples
          --------
  
          >>> # Generate some test data
          >>> from pymbar import testsystems
          >>> (x_n, u_kn, N_k, s_n) = testsystems.HarmonicOscillatorsTestCase().sample(mode='u_kn')
          >>> # Initialize MBAR on data.
          >>> mbar = MBAR(u_kn, N_k)
          >>> # Select the potential we want to compute the PMF for (here, condition 0).
          >>> u_n = u_kn[0, :]
          >>> # Sort into nbins equally-populated bins
          >>> nbins = 10 # number of equally-populated bins to use
          >>> import numpy as np
          >>> N_tot = N_k.sum()
          >>> x_n_sorted = np.sort(x_n) # unroll to n-indices
          >>> bins = np.append(x_n_sorted[0::int(N_tot/nbins)], x_n_sorted.max()+0.1)
          >>> bin_widths = bins[1:] - bins[0:-1]
          >>> bin_n = np.zeros(x_n.shape, np.int64)
          >>> bin_n = np.digitize(x_n, bins) - 1
          >>> # Compute PMF for these unequally-sized bins.
          >>> results = mbar.computePMF(u_n, bin_n, nbins)
          >>> # If we want to correct for unequally-spaced bins to get a PMF on uniform measure
          >>> f_i_corrected = results['f_i'] - np.log(bin_widths)
  
          """
  
          # Verify that no PMF bins are empty -- we can't deal with empty bins,
          # because the free energy is infinite.
          for i in range(nbins):
              if np.sum(bin_n == i) == 0:
                  raise ParameterError(
                      "At least one bin in provided bin_n argument has no samples.  All bins must have samples for free energies to be finite.  Adjust bin sizes or eliminate empty bins to ensure at least one sample per bin.")
          K = self.K
  
          if len(np.shape(u_n)) == 2:
              u_n = kn_to_n(u_n, N_k = self.N_k)
  
          if len(np.shape(bin_n)) == 2:
              bin_n = kn_to_n(bin_n, N_k = self.N_k)
  
          # Compute unnormalized log weights for the given reduced potential
          # u_n.
          log_w_n = self._computeUnnormalizedLogWeights(u_n)
  
          # Compute the free energies for these states.
          f_i = np.zeros([nbins], np.float64)
          df_i = np.zeros([nbins], np.float64)
          for i in range(nbins):
              # Get linear n-indices of samples that fall in this bin.
              indices = np.where(bin_n == i)
  
              # Sanity check.
              if (len(indices) == 0):
                  raise DataError("WARNING: bin %d has no samples -- all bins must have at least one sample." % i)
  
              # Compute dimensionless free energy of occupying state i.
              f_i[i] = - logsumexp(log_w_n[indices])
  
          # Compute uncertainties by forming matrix of W_nk.
          N_k = np.zeros([self.K + nbins], np.int64)
          N_k[0:K] = self.N_k
          W_nk = np.zeros([self.N, self.K + nbins], np.float64)
          W_nk[:, 0:K] = np.exp(self.Log_W_nk)
          for i in range(nbins):
              # Get indices of samples that fall in this bin.
              indices = np.where(bin_n == i)
  
              # Compute normalized weights for this state.
              W_nk[indices, K + i] = np.exp(log_w_n[indices] + f_i[i])
  
          # Compute asymptotic covariance matrix using specified method.
          Theta_ij = self._computeAsymptoticCovarianceMatrix(W_nk, N_k)
  
          # create dictionary to return results
          result_vals = dict()
  
          if (uncertainties == 'from-lowest') or (uncertainties == 'from-specified'):
              # Report uncertainties in free energy difference from a given point
              # on PMF.
  
              if (uncertainties == 'from-lowest'):
                  # Determine bin index with lowest free energy.
                  j = f_i.argmin()
              elif (uncertainties == 'from-specified'):
                  if pmf_reference == None:
                      raise ParameterError(
                          "no reference state specified for PMF using uncertainties = from-specified")
                  else:
                      j = pmf_reference
              # Compute uncertainties with respect to difference in free energy
              # from this state j.
              for i in range(nbins):
                  df_i[i] = math.sqrt(
                      Theta_ij[K + i, K + i] + Theta_ij[K + j, K + j] - 2.0 * Theta_ij[K + i, K + j])
  
              # Shift free energies so that state j has zero free energy.
              f_i -= f_i[j]
  
              # Return dimensionless free energy and uncertainty.
              result_vals['f_i'] = f_i
              result_vals['df_i'] = df_i
  
              return result_vals
  
          elif (uncertainties == 'all-differences'):
              # Report uncertainties in all free energy differences.
  
              diag = Theta_ij.diagonal()
              dii = diag[K, K + nbins]
              d2f_ij = dii + dii.transpose() - 2 * Theta_ij[K:K + nbins, K:K + nbins]
  
              # unsquare uncertainties
              df_ij = np.sqrt(d2f_ij)
  
              # Return dimensionless free energy and uncertainty.
              result_vals['f_i'] = f_i
              result_vals['df_ij'] = df_ij
  
              return result_vals
  
          elif (uncertainties == 'from-normalization'):
              # Determine uncertainties from normalization that \sum_i p_i = 1.
  
              # Compute bin probabilities p_i
              p_i = np.exp(-f_i - logsumexp(-f_i))
  
              # todo -- eliminate triple loop over nbins!
              # Compute uncertainties in bin probabilities.
              d2p_i = np.zeros([nbins], np.float64)
              for k in range(nbins):
                  for i in range(nbins):
                      for j in range(nbins):
                          delta_ik = 1.0 * (i == k)
                          delta_jk = 1.0 * (j == k)
                          d2p_i[k] += p_i[k] * (p_i[i] - delta_ik) * p_i[
                              k] * (p_i[j] - delta_jk) * Theta_ij[K + i, K + j]
  
              # Transform from d2p_i to df_i
              d2f_i = d2p_i / p_i ** 2
              df_i = np.sqrt(d2f_i)
  
              # return free energy and uncertainty
              # Return dimensionless free energy and uncertainty.
              result_vals['f_i'] = f_i
              result_vals['df_i'] = df_i
  
          else:
              raise ParameterError("Uncertainty method '%s' not recognized." % uncertainties)