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.

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