[ENH] Improve performance of isotonic regression by tracking region lower/upper bounds

To keep track of which indices were being merged during PAVA, we were
storing a list of indices and appending during pooling. However, due to
the way the pooling is performed, we only ever append contiguous ranges,
and so we were effectively building [1,...,k], [k+1, ..., k+p],
... arrays.

Simply tracking the upper and lower bound allows us to avoid these
operations and perform the update step in O(1) (the array pop is O(N),
however, which is the dominant term for larger problems).

The somewhat crude benchmark in this diff (adapted from one of the
existing unit tests) shows around 20% improvement in speed for problems
of size 10,000, and ~5% improvement for problems of size 100,000.

Unit tests cover the isotonic regression code fairly well, and all pass
before and after the change.