docs: def detectEquilibration(A_t, fast=True, nskip=1):
      """Automatically detect equilibrated region of a dataset using a heuristic that maximizes number of effectively uncorrelated samples.
  
      Parameters
      ----------
      A_t : np.ndarray 
          timeseries
      nskip : int, optional, default=1
          number of samples to sparsify data by in order to speed equilibration detection
  
      Returns
      -------
      t : int
          start of equilibrated data
      g : float
          statistical inefficiency of equilibrated data
      Neff_max : float
          number of uncorrelated samples   
  
      ToDo
      ----
      Consider implementing a binary search for Neff_max.
  
      Notes
      -----
      If your input consists of some period of equilibration followed by
      a constant sequence, this function treats the trailing constant sequence
      as having Neff = 1.  
  
      Examples
      --------
  
      Determine start of equilibrated data for a correlated timeseries.
  
      >>> from pymbar import testsystems
      >>> A_t = testsystems.correlated_timeseries_example(N=1000, tau=5.0) # generate a test correlated timeseries
      >>> [t, g, Neff_max] = detectEquilibration(A_t) # compute indices of uncorrelated timeseries
  
      Determine start of equilibrated data for a correlated timeseries with a shift.
  
      >>> from pymbar import testsystems
      >>> A_t = testsystems.correlated_timeseries_example(N=1000, tau=5.0) + 2.0 # generate a test correlated timeseries
      >>> B_t = testsystems.correlated_timeseries_example(N=10000, tau=5.0) # generate a test correlated timeseries
      >>> C_t = np.concatenate([A_t, B_t])
      >>> [t, g, Neff_max] = detectEquilibration(C_t, nskip=50) # compute indices of uncorrelated timeseries
  
      """
      T = A_t.size
  
      # Special case if timeseries is constant.
      if A_t.std() == 0.0:
          return (0, 1, 1)  # Changed from Neff=N to Neff=1 after issue #122
  
      g_t = np.ones([T - 1], np.float32)
      Neff_t = np.ones([T - 1], np.float32)
      for t in range(0, T - 1, nskip):
          try:
              g_t[t] = statisticalInefficiency(A_t[t:T], fast=fast)
          except ParameterError:  # Fix for issue https://github.com/choderalab/pymbar/issues/122
              g_t[t] = (T - t + 1)
          Neff_t[t] = (T - t + 1) / g_t[t]
      Neff_max = Neff_t.max()
      t = Neff_t.argmax()
      g = g_t[t]
  
      return (t, g, Neff_max)