docs: def bs_tt_bblock (times, tstarts, tstops, p0=0.05, nbootstrap=512):
      """Bayesian Blocks for time-tagged events with bootstrapping uncertainty
      assessment. THE UNCERTAINTIES ARE NOT VERY GOOD! Arguments:
  
      tstarts        - Array of input bin start times.
      tstops         - Array of input bin stop times.
      times          - Array of event arrival times.
      p0=0.05        - Probability of preferring solutions with additional bins.
      nbootstrap=512 - Number of bootstrap runs to perform.
  
      Returns a Holder with:
  
      blockstarts - Start times of output blocks.
      bsrates     - Mean event rate in each bin from bootstrap analysis.
      bsrstds     - ~Uncertainty: stddev of event rate in each bin from bootstrap analysis.
      counts      - Number of events in each output block.
      finalp0     - Final value of p0, after iteration to minimize `nblocks`.
      ledges      - Times of left edges of output blocks.
      midpoints   - Times of midpoints of output blocks.
      nblocks     - Number of output blocks.
      ncells      - Number of input cells/bins.
      origp0      - Original value of p0.
      rates       - Event rate associated with each block.
      redges      - Times of right edges of output blocks.
      widths      - Width of each output block.
  
      """
      times = np.asarray (times)
      tstarts = np.asarray (tstarts)
      tstops = np.asarray (tstops)
  
      nevents = times.size
      if nevents < 1:
          raise ValueError ('must be given at least 1 event')
  
      info = tt_bblock (tstarts, tstops, times, p0)
  
      # Now bootstrap resample to assess uncertainties on the bin heights. This
      # is the approach recommended by Scargle+.
  
      bsrsums = np.zeros (info.nblocks)
      bsrsumsqs = np.zeros (info.nblocks)
  
      for _ in range (nbootstrap):
          bstimes = times[np.random.randint (0, times.size, times.size)]
          bstimes.sort ()
          bsinfo = tt_bblock (tstarts, tstops, bstimes, p0)
          blocknums = np.minimum (np.searchsorted (bsinfo.redges, info.midpoints),
                                  bsinfo.nblocks - 1)
          samprates = bsinfo.rates[blocknums]
          bsrsums += samprates
          bsrsumsqs += samprates**2
  
      bsrmeans = bsrsums / nbootstrap
      mask = bsrsumsqs / nbootstrap <= bsrmeans**2
      bsrstds = np.sqrt (np.where (mask, 0, bsrsumsqs / nbootstrap - bsrmeans**2))
      info.bsrates = bsrmeans
      info.bsrstds = bsrstds
      return info