Repository: dfm/emcee

Issue:
Incompatible with numpy 2.0.0rc1
Numpy 2 is around the corner with having its first release candidate: https://pypi.org/project/numpy/2.0.0rc1

`emcee` at least uses `VisibleDeprecationWarning` from numpy that will be removed with numpy 2: https://numpy.org/devdocs/release/2.0.0-notes.html

> Warnings and exceptions present in [numpy.exceptions](https://numpy.org/devdocs/reference/routines.exceptions.html#module-numpy.exceptions) (e.g, [ComplexWarning](https://numpy.org/devdocs/reference/generated/numpy.exceptions.ComplexWarning.html#numpy.exceptions.ComplexWarning), [VisibleDeprecationWarning](https://numpy.org/devdocs/reference/generated/numpy.exceptions.VisibleDeprecationWarning.html#numpy.exceptions.VisibleDeprecationWarning)) are no longer exposed in the main namespace.

So with the release of numpy 2, emcee will run into some troubles.

Discussion:
MuellerSeb: Here it is:
```python
import emcee as mc
import numpy as np
import numpy.random as rand
from emcee.state import State

sample_around = 1.0
nwalkers = 50
burn_in = 20
oversampling_factor = 10


def ln_pdf(x):
    return np.log(np.sqrt(np.pi) * np.exp(-((x / 2.0) ** 2)))


init_guess = (
    rand.RandomState(123).rand(nwalkers).reshape((nwalkers, 1)) * sample_around
)
# initialize the sampler
sampler = mc.EnsembleSampler(nwalkers, 1, ln_pdf)
# burn in phase with saving of last position
initial_state = State(init_guess, copy=True)
initial_state.random_state = rand.RandomState(456).get_state()
burn_in_state = sampler.run_mcmc(initial_state=initial_state, nsteps=burn_in)
```
MuellerSeb: In backend.py in L227, i added some print statements to look at the shapes:
```python
        if state.blobs is not None:
            print(self.blobs.shape)
            print(state.blobs.shape)
            self.blobs[self.iteration, :] = state.blobs
```
With np 1.x, I get
```
(20, 50, 0)
(50, 0)
```
With np 2 I get:
```
(20, 50)
(50, 0)
```
And `self.blobs` looks like this:
```
[[b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b''
  b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b''
  b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'']
 [b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b''
  b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b''
  b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'']
 [b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b''
  b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b''
  b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'']
 [b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b''
  b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b''
  b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'']
 [b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b''
  b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b''
  b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'']
 [b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b''
  b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b''
  b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'']
 [b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b''
  b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b''
  b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'']
 [b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b''
  b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b''
  b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'']
 [b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b''
  b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b''
  b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'']
 [b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b''
  b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b''
  b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'']
 [b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b''
  b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b''
  b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'']
 [b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b''
  b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b''
  b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'']
 [b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b''
  b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b''
  b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'']
 [b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b''
  b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b''
  b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'']
 [b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b''
  b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b''
  b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'']
 [b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b''
  b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b''
  b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'']
 [b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b''
  b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b''
  b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'']
 [b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b''
  b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b''
  b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'']
 [b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b''
  b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b''
  b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'']
 [b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b''
  b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b''
  b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'' b'']]
  ```
So `self.blobs` is not correctly generated with numpy 2.
MuellerSeb: And this seems to come from L179 in backend.py:
```python
        if blobs is not None:
            dt = np.dtype((blobs.dtype, blobs.shape[1:]))
            a = np.empty((i, self.nwalkers), dtype=dt)
```
Here, `a` has shape (20, 50, 0) for np 1.x and (20, 50) with numpy 2.
MuellerSeb: So I guess the special dtype is the root of this error: `np.dtype((blobs.dtype, blobs.shape[1:]))`
MuellerSeb: And that is maybe related to these changes in numpy:
- https://numpy.org/devdocs/numpy_2_0_migration_guide.html#changes-to-numpy-data-type-promotion
- https://numpy.org/neps/nep-0050-scalar-promotion.html#nep50
MuellerSeb: Here is what ChatGPT says about the problem:

The Problem

In NumPy 1.x, creating an array with a dtype shaped like (datatype, (0,)) and then attempting to concatenate it along a new axis could lead to an array that still retains the original datatype (float64 in your case), but with an additional dimension of size 0. This behavior effectively results in an array with a shape that includes a 0-sized dimension, which might not be what was intended but works due to how NumPy 1.x handles such cases.

In NumPy 2, the behavior seems to have changed to more strictly interpret the dtype, leading to the creation of an array where the dtype is V (void) or similar, representing a block of memory of the specified size but without a direct numerical datatype associated with it. This change reflects a stricter interpretation of the dtype specification, particularly around structured arrays and zero-sized dimensions.

The Solution

To address this issue and ensure compatibility with both NumPy 1 and 2, you should revise the approach to creating and manipulating these arrays. If the goal is to have self.blobs be a structured array with dimensions that match blobs but starting empty, you need to reconsider how to define dt and how to initialize a.

If blobs.shape[1:] results in a (0,) shape and you still wish to use this structure for consistency, consider explicitly handling this case or rethinking the necessity of the zero-sized dimension. It might be more straightforward to avoid creating a dtype that includes a zero-sized shape unless it's essential for your application.

For applications that don't strictly require the structured array with a zero-sized dimension, you can simplify the dtype definition. If blobs is expected to have additional dimensions that are meaningful, use those directly. If the zero-sized dimension is not essential, you might avoid including it in your dtype definition.

Here's a conceptual fix, assuming you want to preserve the non-zero dimensions of blobs:

```python

if blobs is not None:
    # If blobs.shape[1:] is (0,), this needs handling or reconsideration
    if blobs.shape[1:]:  # This checks if there are meaningful dimensions beyond the first
        dt = np.dtype((blobs.dtype, blobs.shape[1:]))
    else:
        # Adjust the dtype definition to suit your needs when blobs.shape[1:] is (0,)
        dt = blobs.dtype  # Or another appropriate adjustment

    a = np.empty((i, self.nwalkers), dtype=dt)

    if self.blobs is None:
        self.blobs = a
    else:
        self.blobs = np.concatenate((self.blobs, a), axis=0)
```
This adjustment ensures you're consciously deciding on the dtype structure based on the presence of meaningful dimensions in blobs. If handling zero-sized dimensions is crucial for your application, you may need to include conditional logic to manage these cases explicitly, aligning with the behaviors of NumPy 2.

Given the migration to NumPy 2 and its stricter behaviors, it's also beneficial to consult the NumPy 2.0 migration guide for specific guidance on handling dtype and structured array changes.

MuellerSeb: Maybe it is sufficient to always use:

```python
a = np.empty((i, self.nwalkers) + blobs.shape[1:], dtype=blobs.dtype)
```
To not depend on this custom dtype.
MuellerSeb: See: https://github.com/andyfaff/emcee/pull/1
dfm: Thanks for figuring this all out so quickly, @andyfaff! Much appreciated!

Patch:
diff --git a/src/emcee/ensemble.py b/src/emcee/ensemble.py
index 3212099..cfef4a5 100644
--- a/src/emcee/ensemble.py
+++ b/src/emcee/ensemble.py
@@ -489,10 +489,19 @@ class EnsembleSampler(object):
             results = list(map_func(self.log_prob_fn, p))
 
         try:
-            log_prob = np.array([float(l[0]) for l in results])
-            blob = [l[1:] for l in results]
+            # perhaps log_prob_fn returns blobs?
+
+            # deal with the blobs first
+            # if l does not have a len attribute (i.e. not a sequence, no blob)
+            # then a TypeError is raised. However, no error will be raised if
+            # l is a length-1 array, np.array([1.234]). In that case blob
+            # will become an empty list.
+            blob = [l[1:] for l in results if len(l) > 1]
+            if not len(blob):
+                raise IndexError
+            log_prob = np.array([_scalar(l[0]) for l in results])
         except (IndexError, TypeError):
-            log_prob = np.array([float(l) for l in results])
+            log_prob = np.array([_scalar(l) for l in results])
             blob = None
         else:
             # Get the blobs dtype
@@ -502,7 +511,7 @@ class EnsembleSampler(object):
                 try:
                     with warnings.catch_warnings(record=True):
                         warnings.simplefilter(
-                            "error", np.VisibleDeprecationWarning
+                            "error", np.exceptions.VisibleDeprecationWarning
                         )
                         try:
                             dt = np.atleast_1d(blob[0]).dtype
@@ -682,3 +691,16 @@ def ndarray_to_list_of_dicts(
       list of dictionaries of parameters
     """
     return [{key: xi[val] for key, val in key_map.items()} for xi in x]
+
+
+def _scalar(fx):
+    # Make sure a value is a true scalar
+    # 1.0, np.float64(1.0), np.array([1.0]), np.array(1.0)
+    if not np.isscalar(fx):
+        try:
+            fx = np.asarray(fx).item()
+        except (TypeError, ValueError) as e:
+            raise ValueError("log_prob_fn should return scalar") from e
+        return float(fx)
+    else:
+        return float(fx)
diff --git a/tox.ini b/tox.ini
index d759fa6..c566174 100644
--- a/tox.ini
+++ b/tox.ini
@@ -1,12 +1,12 @@
 [tox]
-envlist = py{37,38,39,310}{,-extras},lint
+envlist = py{39,310,311,312}{,-extras},lint
 
 [gh-actions]
 python =
-    3.7: py37
-    3.8: py38
-    3.9: py39-extras
+    3.9: py39
     3.10: py310
+    3.11: py311-extras
+    3.12: py312
 
 [testenv]
 deps = coverage[toml]


Test patch:
diff --git a/.github/workflows/tests.yml b/.github/workflows/tests.yml
index 3a3dc22..063f0c0 100644
--- a/.github/workflows/tests.yml
+++ b/.github/workflows/tests.yml
@@ -16,7 +16,7 @@ jobs:
     runs-on: ${{ matrix.os }}
     strategy:
       matrix:
-        python-version: ["3.7", "3.8", "3.9", "3.10"]
+        python-version: ["3.9", "3.10", "3.11", "3.12"]
         os: ["ubuntu-latest"]
         include:
           - python-version: "3.9"
@@ -49,6 +49,30 @@ jobs:
           COVERALLS_PARALLEL: true
           COVERALLS_FLAG_NAME: ${{ matrix.python-version }}-${{ matrix.os }}
 
+  leading_edge:
+    runs-on: ${{ matrix.os }}
+    strategy:
+      matrix:
+        python-version: ["3.12"]
+        os: ["ubuntu-latest"]
+
+    steps:
+      - name: Checkout
+        uses: actions/checkout@v4
+        with:
+          fetch-depth: 0
+      - name: Setup Python
+        uses: actions/setup-python@v5
+        with:
+          python-version: ${{ matrix.python-version }}
+      - name: Install dependencies
+        run: |
+          python -m pip install -U pip
+          python -m pip install pip install pytest "numpy>=2.0.0rc1"
+          python -m pip install -e.
+      - name: Run tests
+        run: pytest
+
   coverage:
     needs: tests
     runs-on: ubuntu-latest
diff --git a/src/emcee/tests/unit/test_ensemble.py b/src/emcee/tests/unit/test_ensemble.py
index f0568b2..e7981c8 100644
--- a/src/emcee/tests/unit/test_ensemble.py
+++ b/src/emcee/tests/unit/test_ensemble.py
@@ -183,3 +183,37 @@ class TestNamedParameters(TestCase):
         assert results.coords.shape == (n_walkers, len(self.names))
         chain = sampler.chain
         assert chain.shape == (n_walkers, n_steps, len(self.names))
+
+
+class TestLnProbFn(TestCase):
+    # checks that the log_prob_fn can deal with a variety of 'scalar-likes'
+    def lnpdf(self, x):
+        v = np.log(np.sqrt(np.pi) * np.exp(-((x / 2.0) ** 2)))
+        v = float(v[0])
+        assert np.isscalar(v)
+        return v
+
+    def lnpdf_arr1(self, x):
+        v = self.lnpdf(x)
+        return np.array([v])
+
+    def lnpdf_float64(self, x):
+        v = self.lnpdf(x)
+        return np.float64(v)
+
+    def lnpdf_arr0D(self, x):
+        v = self.lnpdf(x)
+        return np.array(v)
+
+    def test_deal_with_scalar_likes(self):
+        rng = np.random.default_rng()
+        fns = [
+            self.lnpdf,
+            self.lnpdf_arr1,
+            self.lnpdf_float64,
+            self.lnpdf_arr0D,
+        ]
+        for fn in fns:
+            init = rng.random((50, 1))
+            sampler = EnsembleSampler(50, 1, fn)
+            _ = sampler.run_mcmc(initial_state=init, nsteps=20)