docs: def bin_bblock (widths, counts, p0=0.05):
      """Fundamental Bayesian Blocks algorithm. Arguments:
  
      widths  - Array of consecutive cell widths.
      counts  - Array of numbers of counts in each cell.
      p0=0.05 - Probability of preferring solutions with additional bins.
  
      Returns a Holder with:
  
      blockstarts - Start times of output blocks.
      counts      - Number of events in each output block.
      finalp0     - Final value of p0, after iteration to minimize `nblocks`.
      nblocks     - Number of output blocks.
      ncells      - Number of input cells/bins.
      origp0      - Original value of p0.
      rates       - Event rate associated with each block.
      widths      - Width of each output block.
  
      """
      widths = np.asarray (widths)
      counts = np.asarray (counts)
      ncells = widths.size
      origp0 = p0
  
      if np.any (widths <= 0):
          raise ValueError ('bin widths must be positive')
      if widths.size != counts.size:
          raise ValueError ('widths and counts must have same size')
      if p0 < 0 or p0 >= 1.:
          raise ValueError ('p0 must lie within [0, 1)')
  
      vedges = np.cumsum (np.concatenate (([0], widths))) # size: ncells + 1
      block_remainders = vedges[-1] - vedges # size: nedges = ncells + 1
      ccounts = np.cumsum (np.concatenate (([0], counts)))
      count_remainders = ccounts[-1] - ccounts
  
      prev_blockstarts = None
      best = np.zeros (ncells, dtype=np.float)
      last = np.zeros (ncells, dtype=np.int)
  
      for _ in range (10):
          # Pluggable num-change-points prior-weight expression:
          ncp_prior = 4 - np.log (p0 / (0.0136 * ncells**0.478))
  
          for r in range (ncells):
              tk = block_remainders[:r+1] - block_remainders[r+1]
              nk = count_remainders[:r+1] - count_remainders[r+1]
  
              # Pluggable fitness expression:
              fit_vec = nlogn (nk, tk)
  
              # This incrementally penalizes partitions with more blocks:
              tmp = fit_vec - ncp_prior
              tmp[1:] += best[:r]
  
              imax = np.argmax (tmp)
              last[r] = imax
              best[r] = tmp[imax]
  
          # different semantics than Scargle impl: our blockstarts is similar to
          # their changepoints, but we always finish with blockstarts[0] = 0.
  
          work = np.zeros (ncells, dtype=int)
          workidx = 0
          ind = last[-1]
  
          while True:
              work[workidx] = ind
              workidx += 1
              if ind == 0:
                  break
              ind = last[ind - 1]
  
          blockstarts = work[:workidx][::-1]
  
          if prev_blockstarts is not None:
              if (blockstarts.size == prev_blockstarts.size and
                  (blockstarts == prev_blockstarts).all ()):
                  break # converged
  
          if blockstarts.size == 1:
              break # can't shrink any farther
  
          # Recommended ad-hoc iteration to favor fewer blocks above and beyond
          # the value of p0:
          p0 = 1. - (1. - p0)**(1. / (blockstarts.size - 1))
          prev_blockstarts = blockstarts
  
      assert blockstarts[0] == 0
      nblocks = blockstarts.size
  
      info = Holder ()
      info.ncells = ncells
      info.nblocks = nblocks
      info.origp0 = origp0
      info.finalp0 = p0
      info.blockstarts = blockstarts
      info.counts = np.empty (nblocks, dtype=np.int)
      info.widths = np.empty (nblocks)
  
      for iblk in range (nblocks):
          cellstart = blockstarts[iblk]
          if iblk == nblocks - 1:
              cellend = ncells - 1
          else:
              cellend = blockstarts[iblk+1] - 1
  
          info.widths[iblk] = widths[cellstart:cellend+1].sum ()
          info.counts[iblk] = counts[cellstart:cellend+1].sum ()
  
      info.rates = info.counts / info.widths
      return info