migrate to numpy>=2, using explicit int32 dtype

evaluate/problems.py

@@ -57,7 +57,7 @@ def volume_asymgauss(loglike, ndim):
             return np.nan
 
         # compute volume of a n-sphere
-        return nsphere_volume(radius, ndim) * np.product(asym_sigma / asym_sigma_max)
+        return nsphere_volume(radius, ndim) * np.prod(asym_sigma / asym_sigma_max)
 
     gradient_asymgauss = gradient_to_center
 

tests/test_popstepsampling.py

@@ -8,6 +8,7 @@
 from ultranest.popstepsampler import generate_cube_oriented_direction, generate_random_direction, generate_cube_oriented_direction_scaled
 from ultranest.popstepsampler import generate_region_oriented_direction, generate_region_random_direction
 from ultranest.popstepsampler import slice_limit_to_unitcube,slice_limit_to_scale
+from ultranest.popstepsampler import int_dtype
 
 def make_region(ndim, us=None, nlive=400):
     if us is None:
@@ -175,9 +176,9 @@ def test_update_slice_sampler():
     The workers should be split among the 2 unfinished points at the end.
     """
 
-    worker_running = np.array([0,0,0,0,1,1,1,1,2,2,2,2])
+    worker_running = np.array([0,0,0,0,1,1,1,1,2,2,2,2], dtype=int_dtype)
     popsize = 12
-    status = np.zeros(12, dtype=int)
+    status = np.zeros(12, dtype=int_dtype)
     status[3:] = 1
     Lmin = 1.
     shrink = 1.0 
@@ -294,8 +295,8 @@ def test_direction_proposal_values():
     #test_stepsampler_cubegausswalk()
     #test_stepsampler_randomSimSlice()
     #test_direction_proposals()
-    #test_slice_limit()
+    test_slice_limit()
     #test_update_slice_sampler()
-    Test_SimpleSliceSampler(4)
+    #Test_SimpleSliceSampler(4)
 
 

ultranest/integrator.py

@@ -43,7 +43,7 @@
 
 __all__ = ['ReactiveNestedSampler', 'NestedSampler', 'read_file', 'warmstart_from_similar_file']
 
-int_t = int
+int_t = np.int32
 
 
 def _get_cumsum_range(pi, dp):

ultranest/mlfriends.pyx

@@ -22,13 +22,17 @@ cimport cython
 from cython.cimports.libc.math import sqrt
 
 
+ctypedef np.int32_t decl_int_t
+int_dtype = np.int32
+
+
 @cython.boundscheck(False)
 @cython.wraparound(False)
 cdef count_nearby(
     np.ndarray[np.float_t, ndim=2] apts,
     np.ndarray[np.float_t, ndim=2] bpts,
     np.float_t radiussq,
-    np.ndarray[np.int_t, ndim=1] nnearby
+    np.ndarray[decl_int_t, ndim=1] nnearby
 ):
     """Count the number of points in ``apts`` within square radius ``radiussq`` for each point ``b`` in `bpts``.
 
@@ -140,7 +144,7 @@ def find_nearby(
     np.ndarray[np.float_t, ndim=2] apts,
     np.ndarray[np.float_t, ndim=2] bpts,
     np.float_t radiussq,
-    np.ndarray[np.int_t, ndim=1] nnearby
+    np.ndarray[decl_int_t, ndim=1] nnearby
 ):
     """Gets the index of a point in `a` within square radius `radiussq`, for each point `b` in `bpts`.
 
@@ -224,7 +228,7 @@ cdef float compute_maxradiussq(np.ndarray[np.float_t, ndim=2] apts, np.ndarray[n
 @cython.wraparound(False)
 def compute_mean_pair_distance(
     np.ndarray[np.float_t, ndim=2] pts,
-    np.ndarray[np.int_t, ndim=1] clusterids
+    np.ndarray[decl_int_t, ndim=1] clusterids
 ):
     """Compute the average distance between pairs of points.
     Pairs from different clusters are excluded in the computation.
@@ -272,12 +276,12 @@ cdef _update_clusters(
     np.ndarray[np.float_t, ndim=2] upoints,
     np.ndarray[np.float_t, ndim=2] tpoints,
     np.float_t maxradiussq,
-    np.ndarray[np.int_t, ndim=1] clusterids,
+    np.ndarray[decl_int_t, ndim=1] clusterids,
 ):
     """same signature as ``update_clusters()``, see there."""
     assert upoints.shape[0] == tpoints.shape[0], ('different number of points', upoints.shape[0], tpoints.shape[0])
     assert upoints.shape[1] == tpoints.shape[1], ('different dimensionality of points', upoints.shape[1], tpoints.shape[1])
-    clusteridxs = np.zeros(len(tpoints), dtype=int)
+    clusteridxs = np.zeros(len(tpoints), dtype=int_dtype)
     currentclusterid = 1
     i = 0
     # avoid issues when old clusterids are from a longer array
@@ -295,7 +299,7 @@ cdef _update_clusters(
             break
 
         nonmembers = tpoints[nonmembermask,:]
-        idnearby = np.empty(len(nonmembers), dtype=int)
+        idnearby = np.empty(len(nonmembers), dtype=int_dtype)
         members = tpoints[clusteridxs == currentclusterid,:]
         find_nearby(members, nonmembers, maxradiussq, idnearby)
         # print('merging %d into cluster %d of size %d' % (np.count_nonzero(nnearby), currentclusterid, len(members)))
@@ -375,8 +379,9 @@ def update_clusters(
     Clustering is performed on a transformed coordinate space (`tpoints`).
     Returned values are based on upoints.
     """
+    print(clusterids)
     if clusterids is None:
-        clusterids = np.zeros(len(tpoints), dtype=int)
+        clusterids = np.zeros(len(tpoints), dtype=int_dtype)
     return _update_clusters(upoints, tpoints, maxradiussq, clusterids)
 
 
@@ -409,7 +414,7 @@ def make_eigvals_positive(
         raise e
     mask = w < max(1.e-10, 1e-300**(1. / len(a)))
     if np.any(mask):
-        nzprod = np.product(w[~mask])  # product of nonzero eigenvalues
+        nzprod = np.prod(w[~mask])  # product of nonzero eigenvalues
         nzeros = mask.sum()  # number of zero eigenvalues
         w[mask] = (targetprod / nzprod) ** (1. / nzeros)  # adjust zero eigvals
         a = np.dot(np.dot(v, np.diag(w)), np.linalg.inv(v))  # re-form cov
@@ -568,7 +573,7 @@ class ScalingLayer(object):
         """Updates the cluster id assigned to each point."""
         if clusterids is None and self.clusterids is None and npoints is not None:
             # for the beginning, set cluster ids to one for all points
-            clusterids = np.ones(npoints, dtype=int)
+            clusterids = np.ones(npoints, dtype=int_dtype)
         if clusterids is not None:
             # if we have a value, update
             self.clusterids = clusterids
@@ -846,7 +851,7 @@ class LocalAffineLayer(AffineLayer):
         return s
 
 
-def vol_prefactor(np.int_t n):
+def vol_prefactor(int n):
     """Volume constant for an ``n``-dimensional sphere.
 
     for ``n`` even:  $$    (2pi)^(n    /2) / (2 * 4 * ... * n)$$
@@ -1080,7 +1085,7 @@ class MLFriends(object):
         v = self.unormed[idx,:] + v * self.maxradiussq**0.5
 
         # count how many are around
-        nnearby = np.empty(nsamples, dtype=int)
+        nnearby = np.empty(nsamples, dtype=int_dtype)
         count_nearby(self.unormed, v, self.maxradiussq, nnearby)
         vmask = np.random.uniform(high=nnearby) < 1
         w = self.transformLayer.untransform(v[vmask,:])
@@ -1102,7 +1107,7 @@ class MLFriends(object):
         wmask = self.inside_ellipsoid(u)
         # check if inside region in transformed space
         v = self.transformLayer.transform(u[wmask,:])
-        idnearby = np.empty(len(v), dtype=int)
+        idnearby = np.empty(len(v), dtype=int_dtype)
         find_nearby(self.unormed, v, self.maxradiussq, idnearby)
         vmask = idnearby >= 0
         return u[wmask,:][vmask,:]
@@ -1117,7 +1122,7 @@ class MLFriends(object):
         N, ndim = self.u.shape
         # draw from rectangle in transformed space
         v = np.random.uniform(self.bbox_lo - self.maxradiussq, self.bbox_hi + self.maxradiussq, size=(nsamples, ndim))
-        idnearby = np.empty(nsamples, dtype=int)
+        idnearby = np.empty(nsamples, dtype=int_dtype)
         find_nearby(self.unormed, v, self.maxradiussq, idnearby)
         vmask = idnearby >= 0
 
@@ -1149,7 +1154,7 @@ class MLFriends(object):
 
         wmask = np.logical_and(w > 0, w < 1).all(axis=1)
         v = self.transformLayer.transform(w[wmask,:])
-        idnearby = np.empty(len(v), dtype=int)
+        idnearby = np.empty(len(v), dtype=int_dtype)
         find_nearby(self.unormed, v, self.maxradiussq, idnearby)
         vmask = idnearby >= 0
 
@@ -1200,7 +1205,7 @@ class MLFriends(object):
         if mask.any():
             # additionally require points to be near neighbours
             bpts = self.transformLayer.transform(pts[mask,:])
-            idnearby = np.empty(len(bpts), dtype=int)
+            idnearby = np.empty(len(bpts), dtype=int_dtype)
             find_nearby(self.unormed, bpts, self.maxradiussq, idnearby)
             mask[mask] = idnearby >= 0
 

ultranest/popstepsampler.py

@@ -22,7 +22,7 @@
                                  update_vectorised_slice_sampler)
 from ultranest.utils import submasks
 
-int_t = int
+int_dtype = np.int32
 
 
 def unitcube_line_intersection(ray_origin, ray_direction):
@@ -432,7 +432,7 @@ def _setup(self, ndim):
         self.allL = np.zeros((self.popsize, self.nsteps + 1)) + np.nan
         self.currentt = np.zeros(self.popsize) + np.nan
         self.currentv = np.zeros((self.popsize, ndim)) + np.nan
-        self.generation = np.zeros(self.popsize, dtype=int_t) - 1
+        self.generation = np.zeros(self.popsize, dtype=int_dtype) - 1
         self.current_left = np.zeros(self.popsize)
         self.current_right = np.zeros(self.popsize)
         self.searching_left = np.zeros(self.popsize, dtype=bool)
@@ -936,9 +936,9 @@ def __next__(
                 # Slice bounds for each points
                 tleft, tright = self.slice_limit(tleft_unitcube,tright_unitcube)
                 # Index of the workers working concurrently
-                worker_running = np.arange(0, self.popsize, 1, dtype=int_t)
+                worker_running = np.arange(self.popsize, dtype=int_dtype)
                 # Status indicating if a points has already find its next position
-                status = np.zeros(self.popsize, dtype=int_t)  # one for success, zero for running
+                status = np.zeros(self.popsize, dtype=int_dtype)  # one for success, zero for running
 
                 # Loop until each points has found its next position or we reached 100 iterations
                 for _it in range(self.max_it):
@@ -955,6 +955,7 @@ def __next__(
                     nc += self.popsize
 
                     # Updating the pool of points based on the newly sampled points
+                    print(worker_running.dtype, status.dtype)
                     tleft, tright, worker_running, status, allu, allL, allp, n_discarded_it = update_vectorised_slice_sampler(
                         t, tleft, tright, proposed_L, proposed_u, proposed_p, worker_running, status, Lmin, self.shrink_factor,
                         allu, allL, allp, self.popsize)

ultranest/stepfuncs.pyx

@@ -13,6 +13,10 @@ cimport cython
 from cython.parallel import prange
 
 
+ctypedef np.int32_t decl_int_t
+int_dtype = np.int32
+
+
 @cython.boundscheck(False)
 @cython.wraparound(False)
 cdef _within_unit_cube(
@@ -332,7 +336,7 @@ def step_back(Lmin, allL, generation, currentt, log=False):
 
 cdef _fill_directions(
     np.ndarray[np.float_t, ndim=2] v,
-    np.ndarray[np.int_t, ndim=1] indices,
+    np.ndarray[decl_int_t, ndim=1] indices,
     float scale
 ):
     cdef size_t nsamples = v.shape[0]
@@ -361,7 +365,7 @@ def generate_cube_oriented_direction(ui, region, scale=1):
     nsamples, ndim = ui.shape
     v = np.zeros((nsamples, ndim))
     # choose axis
-    j = np.random.randint(ndim, size=nsamples)
+    j = np.random.randint(ndim, size=nsamples, dtype=int_dtype)
     _fill_directions(v, j, scale)
     return v
 
@@ -388,7 +392,7 @@ def generate_cube_oriented_direction_scaled(ui, region, scale=1):
     v = np.zeros((nsamples, ndim))
     scales = region.u.std(axis=0)
     # choose axis
-    j = np.random.randint(ndim, size=nsamples)
+    j = np.random.randint(ndim, size=nsamples, dtype=int_dtype)
     _fill_directions(v, j, scale)
     v *= scales[j].reshape((-1, 1))
     return v
@@ -439,7 +443,7 @@ def generate_region_oriented_direction(ui, region, scale=1):
     """
     nsamples, ndim = ui.shape
     # choose axis in transformed space:
-    j = np.random.randint(ndim, size=nsamples)
+    j = np.random.randint(ndim, size=nsamples, dtype=int_dtype)
     v = region.transformLayer.axes[j] * scale
     return v
 
@@ -490,8 +494,8 @@ def generate_differential_direction(ui, region, scale=1):
     nsamples, ndim = ui.shape
     nlive, ndim = region.u.shape
     # choose pair
-    i = np.random.randint(nlive, size=nsamples)
-    i2 = np.random.randint(nlive - 1, size=nsamples)
+    i = np.random.randint(nlive, size=nsamples, dtype=int_dtype)
+    i2 = np.random.randint(nlive - 1, size=nsamples, dtype=int_dtype)
     i2[i2 >= i] += 1
 
     # compute difference vector
@@ -530,14 +534,22 @@ def generate_mixture_random_direction(ui, region, scale=1):
 
 @cython.boundscheck(False)
 @cython.wraparound(False)
-cpdef tuple update_vectorised_slice_sampler(\
-                          np.ndarray[np.float_t, ndim=1] t, np.ndarray[np.float_t, ndim=1] tleft,\
-                          np.ndarray[np.float_t, ndim=1] tright, np.ndarray[np.float_t, ndim=1] proposed_L,\
-                          np.ndarray[np.float_t, ndim=2] proposed_u, np.ndarray[np.float_t, ndim=2] proposed_p,\
-                          np.ndarray[np.int_t, ndim=1] worker_running, np.ndarray[np.int_t, ndim=1] status,\
-                          np.float_t Likelihood_threshold,np.float_t shrink_factor, np.ndarray[np.float_t, ndim=2] allu,\
-                          np.ndarray[np.float_t, ndim=1] allL, np.ndarray[np.float_t, ndim=2] allp, int popsize):
-
+cpdef tuple update_vectorised_slice_sampler(
+    np.ndarray[np.float_t, ndim=1] t,
+    np.ndarray[np.float_t, ndim=1] tleft,
+    np.ndarray[np.float_t, ndim=1] tright,
+    np.ndarray[np.float_t, ndim=1] proposed_L,
+    np.ndarray[np.float_t, ndim=2] proposed_u,
+    np.ndarray[np.float_t, ndim=2] proposed_p,
+    np.ndarray[decl_int_t, ndim=1] worker_running,
+    np.ndarray[decl_int_t, ndim=1] status,
+    np.float_t Likelihood_threshold,
+    np.float_t shrink_factor,
+    np.ndarray[np.float_t, ndim=2] allu,
+    np.ndarray[np.float_t, ndim=1] allL,
+    np.ndarray[np.float_t, ndim=2] allp,
+    int popsize
+):
     """Update the slice sampler state of each walker in the populations.
 
     Parameters