docs: def detectEquilibration_binary_search(A_t, bs_nodes=10):
      """Automatically detect equilibrated region of a dataset using a heuristic that maximizes number of effectively uncorrelated samples.
  
      Parameters
      ----------
      A_t : np.ndarray 
          timeseries
      
      bs_nodes : int > 4
          number of geometrically distributed binary search nodes
  
      Returns
      -------
      t : int
          start of equilibrated data
      g : float
          statistical inefficiency of equilibrated data
      Neff_max : float
          number of uncorrelated samples
          
      Notes
      -----
      Finds the discard region (t) by a binary search on the range of
      possible lagtime values, with logarithmic spacings.  This will give
      a local maximum.  The global maximum is not guaranteed, but will
      likely be found if the N_eff[t] varies smoothly.
  
      """
      assert bs_nodes > 4, "Number of nodes for binary search must be > 4"
      T = A_t.size
  
      # Special case if timeseries is constant.
      if A_t.std() == 0.0:
          return (0, 1, T)
  
      start = 1
      end = T - 1
      n_grid = min(bs_nodes, T)
      
      while True:
          time_grid = np.unique((10 ** np.linspace(np.log10(start), np.log10(end), n_grid)).round().astype('int')) 
          g_t = np.ones(time_grid.size)
          Neff_t = np.ones(time_grid.size)
          
          for k, t in enumerate(time_grid):
              if t < T-1:
                  g_t[k] = statisticalInefficiency_fft(A_t[t:], memsafe=True)
                  Neff_t[k] = (T - t + 1) / g_t[k]
  
          Neff_max = Neff_t.max()
          k = Neff_t.argmax()
          t = time_grid[k]
          g = g_t[k]
          
          if (end - start < 4):
              break
          
          if k == 0:
              start = time_grid[0]
              end = time_grid[1]
          elif k == time_grid.size - 1:
              start = time_grid[-2]
              end = time_grid[-1]
          else:
              start = time_grid[k - 1]
              end = time_grid[k + 1]
          
      return (t, g, Neff_max)