docs: def computePMF_states(self, u_n, bin_n, nbins):
          """Compute the free energy of occupying a number of bins.
  
          This implementation defines each bin as a separate thermodynamic state.
  
          Parameters
          ----------
          u_n : np.ndarray, float, shape=(K, N)
              u_n[n] is the reduced potential energy of snapshot n for which the PMF is to be computed.
          bin_n : np.ndarray, int, shape=(N)
              bin_n[n] is the bin index of snapshot n.  bin_n can assume a value in range(0,nbins)
          nbins : int
              The number of bins
          fmax : float, optional
              The maximum value of the free energy, used for an empty bin (default: 1000)
  
          Returns
          -------
          f_i : np.ndarray, float, shape=(K)
              f_i[i] is the dimensionless free energy of state i, relative to the state of lowest free energy
          d2f_ij : np.ndarray, float, shape=(K)
              d2f_ij[i,j] is the uncertainty in the difference of (f_i - f_j)
  
          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.
  
          WARNING!
          This method is EXPERIMENTAL and should be used at your own risk.
  
          """
  
          # 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
  
          # Unroll to n-indices
          log_w_n = self._computeUnnormalizedLogWeights(u_n)
  
          # Compute the free energies for these states.
          f_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[self.indices] == i)[0]
  
              # 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] = - _logsum(log_w_n[indices])
  
          # Shift so that f_i.min() = 0
          f_i_min = f_i.min()
          f_i -= f_i.min()
  
          if self.verbose:
              print("bins f_i = ")
              print(f_i)
  
          # Compute uncertainties by forming matrix of W_nk.
          if self.verbose:
              print("Forming W_nk matrix...")
          N_k = np.zeros([self.K + nbins], np.int32)
          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[self.indices] == i)[0]
  
              if self.verbose:
                  print("bin %5d count = %10d" % (i, len(indices)))
  
              # Compute normalized weights for this state.
              W_nk[indices, K + i] = np.exp(
                  log_w_n[indices] + f_i[i] + f_i_min)
  
          # Compute asymptotic covariance matrix using specified method.
          Theta_ij = self._computeAsymptoticCovarianceMatrix(W_nk, N_k)
  
          # Compute uncertainties with respect to difference in free energy from
          # this state j.
          diag = Theta_ij.diagonal()
          dii = diag[0, K:K + nbins]
          d2f_ij = dii + dii.transpose() - 2 * Theta_ij[K:K + nbins, K:K + nbins]
  
          # Return dimensionless free energy and uncertainty.
          return (f_i, d2f_ij)