Repository: astropy/astropy

Issue:
Bayesian blocks returns a single edge
I stumbled on an edge case which is causing problems in my app. First, here's an example of expected behavior:

```python
import numpy as np
from astropy.stats import bayesian_blocks

np.random.seed(0)
values = np.concatenate([np.random.rand(100), 0.5 * np.random.rand(100)])
bins = bayesian_blocks(values)

assert values.min() == bins[0]
assert values.max() == bins[-1]
```
These assertions pass successfully. However, they fail for this edge case:
```python
values = np.array([1,1,1,1,1,1,1,1,2])
bins = bayesian_blocks(values)
bins    # produces array([1.])
```
The "edge" in this case turns out to be the difference between 7 and 8 ones. Deleting just one 1 gives
```python
values = np.array([1,1,1,1,1,1,1,2])
bins = bayesian_blocks(values)
bins    # produces array([1., 2.])
```

EDIT: Snippet formatting

Discussion:
codecov[bot]: # [Codecov](https://codecov.io/gh/astropy/astropy/pull/8560?src=pr&el=h1) Report
> Merging [#8560](https://codecov.io/gh/astropy/astropy/pull/8560?src=pr&el=desc) into [master](https://codecov.io/gh/astropy/astropy/commit/8824789d857638fc320d35b53d9c17f4fa9bb889?src=pr&el=desc) will **decrease** coverage by `<.01%`.
> The diff coverage is `66.66%`.

[![Impacted file tree graph](https://codecov.io/gh/astropy/astropy/pull/8560/graphs/tree.svg?width=650&token=eEhbJxwU0j&height=150&src=pr)](https://codecov.io/gh/astropy/astropy/pull/8560?src=pr&el=tree)

```diff
@@            Coverage Diff             @@
##           master    #8560      +/-   ##
==========================================
- Coverage   86.81%   86.81%   -0.01%     
==========================================
  Files         386      386              
  Lines       58162    58164       +2     
  Branches     1060     1060              
==========================================
  Hits        50493    50493              
- Misses       7054     7056       +2     
  Partials      615      615
```


| [Impacted Files](https://codecov.io/gh/astropy/astropy/pull/8560?src=pr&el=tree) | Coverage Δ | |
|---|---|---|
| [astropy/stats/bayesian\_blocks.py](https://codecov.io/gh/astropy/astropy/pull/8560/diff?src=pr&el=tree#diff-YXN0cm9weS9zdGF0cy9iYXllc2lhbl9ibG9ja3MucHk=) | `94.7% <66.66%> (-0.61%)` | :arrow_down: |
| [astropy/samp/hub.py](https://codecov.io/gh/astropy/astropy/pull/8560/diff?src=pr&el=tree#diff-YXN0cm9weS9zYW1wL2h1Yi5weQ==) | `77.38% <0%> (-0.14%)` | :arrow_down: |

------

[Continue to review full report at Codecov](https://codecov.io/gh/astropy/astropy/pull/8560?src=pr&el=continue).
> **Legend** - [Click here to learn more](https://docs.codecov.io/docs/codecov-delta)
> `Δ = absolute <relative> (impact)`, `ø = not affected`, `? = missing data`
> Powered by [Codecov](https://codecov.io/gh/astropy/astropy/pull/8560?src=pr&el=footer). Last update [8824789...67cc8bb](https://codecov.io/gh/astropy/astropy/pull/8560?src=pr&el=lastupdated). Read the [comment docs](https://docs.codecov.io/docs/pull-request-comments).

codecov[bot]: # [Codecov](https://codecov.io/gh/astropy/astropy/pull/8560?src=pr&el=h1) Report
> Merging [#8560](https://codecov.io/gh/astropy/astropy/pull/8560?src=pr&el=desc) into [master](https://codecov.io/gh/astropy/astropy/commit/8824789d857638fc320d35b53d9c17f4fa9bb889?src=pr&el=desc) will **decrease** coverage by `<.01%`.
> The diff coverage is `66.66%`.

[![Impacted file tree graph](https://codecov.io/gh/astropy/astropy/pull/8560/graphs/tree.svg?width=650&token=eEhbJxwU0j&height=150&src=pr)](https://codecov.io/gh/astropy/astropy/pull/8560?src=pr&el=tree)

```diff
@@            Coverage Diff             @@
##           master    #8560      +/-   ##
==========================================
- Coverage   86.81%   86.81%   -0.01%     
==========================================
  Files         386      386              
  Lines       58162    58164       +2     
  Branches     1060     1060              
==========================================
  Hits        50493    50493              
- Misses       7054     7056       +2     
  Partials      615      615
```


| [Impacted Files](https://codecov.io/gh/astropy/astropy/pull/8560?src=pr&el=tree) | Coverage Δ | |
|---|---|---|
| [astropy/stats/bayesian\_blocks.py](https://codecov.io/gh/astropy/astropy/pull/8560/diff?src=pr&el=tree#diff-YXN0cm9weS9zdGF0cy9iYXllc2lhbl9ibG9ja3MucHk=) | `94.7% <66.66%> (-0.61%)` | :arrow_down: |
| [astropy/samp/hub.py](https://codecov.io/gh/astropy/astropy/pull/8560/diff?src=pr&el=tree#diff-YXN0cm9weS9zYW1wL2h1Yi5weQ==) | `77.38% <0%> (-0.14%)` | :arrow_down: |

------

[Continue to review full report at Codecov](https://codecov.io/gh/astropy/astropy/pull/8560?src=pr&el=continue).
> **Legend** - [Click here to learn more](https://docs.codecov.io/docs/codecov-delta)
> `Δ = absolute <relative> (impact)`, `ø = not affected`, `? = missing data`
> Powered by [Codecov](https://codecov.io/gh/astropy/astropy/pull/8560?src=pr&el=footer). Last update [8824789...67cc8bb](https://codecov.io/gh/astropy/astropy/pull/8560?src=pr&el=lastupdated). Read the [comment docs](https://docs.codecov.io/docs/pull-request-comments).

astrofrog: @jakevdp - as the original author of the bayesian blocks code, could you review this?

@zkurtz - even if the fix has to be changed, it would be good to add a test of a case where this issue happened to make sure that it is fixed and remains so in future.
jakevdp: It's unclear to me which changes I'm meant to be reviewing. Could you clean up the branch?
bsipocz: @jakevdp - I've rebase to keep only the relevant changes here.

@zkurtz - the usual suggestion is to use feature branches for new development otherwise one can end up in a messy situation with pulling in new changes from upstream master. Here is a detailed guideline about the howto-s http://docs.astropy.org/en/latest/development/workflow/development_workflow.html
jakevdp: The change looks plausible, but I have to admit I can't reason very well about the implications. While I did author this code, it was mostly just a transcription of the algorithm laid out in the Scargle paper.
astropy-bot[bot]:  Hi humans :wave: - this pull request hasn't had any new commits for approximately 5 months. **I plan to close this in a month if the pull request doesn't have any new commits by then.**

In lieu of a stalled pull request, please consider closing this and open an issue instead if a reminder is needed to revisit in the future. Maintainers may also choose to add `keep-open` label to keep this PR open but it is discouraged unless absolutely necessary.

If this PR still needs to be reviewed, as an author, you can rebase it to reset the clock.

*If you believe I commented on this pull request incorrectly, please report this [here](https://github.com/OpenAstronomy/baldrick/issues).*

Patch:
diff --git a/CHANGES.rst b/CHANGES.rst
index 7035d4a298..f1a92dbe57 100644
--- a/CHANGES.rst
+++ b/CHANGES.rst
@@ -476,6 +476,8 @@ astropy.samp
 astropy.stats
 ^^^^^^^^^^^^^
 
+- Fixed a bug where ``bayesian_blocks`` returned a single edge. [#8560]
+
 astropy.table
 ^^^^^^^^^^^^^
 
diff --git a/astropy/stats/bayesian_blocks.py b/astropy/stats/bayesian_blocks.py
index 90b16bf15f..d4f5406552 100644
--- a/astropy/stats/bayesian_blocks.py
+++ b/astropy/stats/bayesian_blocks.py
@@ -37,7 +37,6 @@ References
 .. [1] http://adsabs.harvard.edu/abs/2012arXiv1207.5578S
 .. [2] http://astroml.org/ https://github.com//astroML/astroML/
 """
-
 import warnings
 
 import numpy as np
@@ -385,12 +384,14 @@ class FitnessFunc:
         change_points = np.zeros(N, dtype=int)
         i_cp = N
         ind = N
-        while True:
+        while i_cp > 0:
             i_cp -= 1
             change_points[i_cp] = ind
             if ind == 0:
                 break
             ind = last[ind - 1]
+        if i_cp == 0:
+            change_points[i_cp] = 0
         change_points = change_points[i_cp:]
 
         return edges[change_points]


Test patch:
diff --git a/astropy/stats/tests/test_bayesian_blocks.py b/astropy/stats/tests/test_bayesian_blocks.py
index b2593eb9aa..cd69000aca 100644
--- a/astropy/stats/tests/test_bayesian_blocks.py
+++ b/astropy/stats/tests/test_bayesian_blocks.py
@@ -162,3 +162,13 @@ def test_fitness_function_results():
     edges = bayesian_blocks(t, x_obs, sigma, fitness='measures',
                             gamma=gamma_sel)
     assert_allclose(edges, expected)
+
+
+def test_zero_change_points(rseed=0):
+    ''' Ensure that edges contains both endpoints when there are no change points '''
+    np.random.seed(rseed)
+    # Using the failed edge case from https://github.com/astropy/astropy/issues/8558
+    values = np.array([1, 1, 1, 1, 1, 1, 1, 1, 2])
+    bins = bayesian_blocks(values)
+    assert values.min() == bins[0]
+    assert values.max() == bins[-1]