docs: def BARzero(w_F, w_R, DeltaF):
      """A function that when zeroed is equivalent to the solution of
      the Bennett acceptance ratio.
  
      from http://journals.aps.org/prl/pdf/10.1103/PhysRevLett.91.140601
      D_F = M + w_F - Delta F
      D_R = M + w_R - Delta F
  
      we want:
      \sum_N_F (1+exp(D_F))^-1 = \sum N_R N_R <(1+exp(-D_R))^-1>
      ln \sum N_F (1+exp(D_F))^-1>_F = \ln \sum N_R exp((1+exp(-D_R))^(-1)>_R
      ln \sum N_F (1+exp(D_F))^-1>_F - \ln \sum N_R exp((1+exp(-D_R))^(-1)>_R = 0
  
      Parameters
      ----------
      w_F : np.ndarray
          w_F[t] is the forward work value from snapshot t.
          t = 0...(T_F-1)  Length T_F is deduced from vector.
      w_R : np.ndarray
          w_R[t] is the reverse work value from snapshot t.
          t = 0...(T_R-1)  Length T_R is deduced from vector.
      DeltaF : float
          Our current guess
  
      Returns
      -------
      fzero : float
          a variable that is zeroed when DeltaF satisfies BAR.
  
      Examples
      --------
      Compute free energy difference between two specified samples of work values.
  
      >>> from pymbar import testsystems
      >>> [w_F, w_R] = testsystems.gaussian_work_example(mu_F=None, DeltaF=1.0, seed=0)
      >>> DeltaF = BARzero(w_F, w_R, 0.0)
  
      """
  
      np.seterr(over='raise')  # raise exceptions to overflows
      w_F = np.array(w_F, np.float64)
      w_R = np.array(w_R, np.float64)
      DeltaF = float(DeltaF)
  
      # Recommended stable implementation of BAR.
  
      # Determine number of forward and reverse work values provided.
      T_F = float(w_F.size)  # number of forward work values
      T_R = float(w_R.size)  # number of reverse work values
  
      # Compute log ratio of forward and reverse counts.
      M = np.log(T_F / T_R)
  
      # Compute log numerator. We have to watch out for overflows.  We
      # do this by making sure that 1+exp(x) doesn't overflow, choosing
      # to always exponentiate a negative number.
  
      # log f(W) = - log [1 + exp((M + W - DeltaF))]
      #          = - log ( exp[+maxarg] [exp[-maxarg] + exp[(M + W - DeltaF) - maxarg]] )
      #          = - maxarg - log(exp[-maxarg] + exp[(M + W - DeltaF) - maxarg])
      # where maxarg = max((M + W - DeltaF), 0)
  
      exp_arg_F = (M + w_F - DeltaF)
      # use boolean logic to zero out the ones that are less than 0, but not if greater than zero.
      max_arg_F = np.choose(np.less(0.0, exp_arg_F), (0.0, exp_arg_F))
      try:
          log_f_F = - max_arg_F - np.log(np.exp(-max_arg_F) + np.exp(exp_arg_F - max_arg_F))
      except:
          # give up; if there's overflow, return zero
          print("The input data results in overflow in BAR")
          return np.nan
      log_numer = logsumexp(log_f_F)
  
      # Compute log_denominator.
      # log f(R) = - log [1 + exp(-(M + W - DeltaF))]
      #          = - log ( exp[+maxarg] [exp[-maxarg] + exp[(M + W - DeltaF) - maxarg]] )
      #          = - maxarg - log[exp[-maxarg] + (T_F/T_R) exp[(M + W - DeltaF) - maxarg]]
      # where maxarg = max( -(M + W - DeltaF), 0)
  
      exp_arg_R = -(M - w_R - DeltaF)
      # use boolean logic to zero out the ones that are less than 0, but not if greater than zero.
      max_arg_R = np.choose(np.less(0.0, exp_arg_R), (0.0, exp_arg_R))
      try:
          log_f_R = - max_arg_R - np.log(np.exp(-max_arg_R) + np.exp(exp_arg_R - max_arg_R))
      except:
          print("The input data results in overflow in BAR")
          return np.nan
      log_denom = logsumexp(log_f_R)
  
      # This function must be zeroed to find a root
      fzero = log_numer - log_denom
  
      np.seterr(over='warn')  # return options to standard settings so we don't disturb other functionality.
      return fzero