numpy 2 compatible (#249)

* These changes enable calcos to run under numpy 2, and also lower the differences
between running with numpy 2 and numpy 1.x to more tolerable levels

* Remove npy_promotion_state calls and fix an indent

* Changes for numpy 2:
 - casting changes required for C extensions
 - casting numpy float32 scalars to float64 to ensure that arithmetic is done using float64 precision
 - remove some conditionals in the C extension that were put in for python 2/3 compatibility

* Don't cast arrays that are None

calcos/concurrent.py

@@ -957,7 +957,7 @@ def findFlash(self, output=None, update=True):
         # Get time as a single precision array, for ccos.getstartstop.
         time = cosutil.getColCopy(data=self.events, column="time")
         dq = np.zeros(len(time), dtype=np.int16)
-        nbins = int(math.ceil((time[-1] - time[0]) / self.delta_t))
+        nbins = int(math.ceil((np.float64(time[-1]) - np.float64(time[0])) / self.delta_t))
         istart = np.zeros(nbins, dtype=np.int32)
         istop = np.zeros(nbins, dtype=np.int32)
         src_counts = np.zeros(nbins, dtype=np.int32)
@@ -1109,19 +1109,18 @@ def findLampOn(self, cutoff, t0):
             for k in range(k_on, k_off+1):
                 if self.src_counts[k] > cutoff:
                     lamp_is_on = True
-                    t = t0 + k * self.delta_t
+                    t = np.float64(t0) + k * self.delta_t
                     lamp_on.append(t)
                     break
             if not lamp_is_on:
                 continue
             # search for the actual lamp turn-off time
             for k in range(k_off, k_on-1, -1):
                 if self.src_counts[k] > cutoff:
-                    t = t0 + (k+1) * self.delta_t
+                    t = np.float64(t0) + (k+1) * self.delta_t
                     t = min(t, self.time[-1])
                     lamp_off.append(t)
                     break
-
         if len(lamp_on) != len(lamp_off):
             raise RuntimeError("Internal error:  len(lamp_on) = %d, "
                                "len(lamp_off) = %d" % \
@@ -1168,8 +1167,8 @@ def getIndices(self, nbins, t0, lamp_on_i, lamp_off_i):
             interval (t0, t0+delta_t).
         """
 
-        k_on = (lamp_on_i - t0) / self.delta_t
-        k_off = (lamp_off_i - t0) / self.delta_t
+        k_on = (lamp_on_i - np.float64(t0)) / self.delta_t
+        k_off = (lamp_off_i - np.float64(t0)) / self.delta_t
         k_on = int(math.floor(k_on))
         k_off = int(math.ceil(k_off))
         if k_on >= nbins or k_off < 0:

calcos/cosutil.py

@@ -1715,6 +1715,10 @@ def fuvFlagOutOfBounds(hdr, dq_array, info, switches,
         y_upper = (y1 - 1.) + 0. * x_upper
         x_left  = x0 + 0. * y_left
         x_right = (x1 - 1.) + 0. * y_right
+        y_lower = np.float32(y_lower)
+        y_upper = np.float32(y_upper)
+        x_left = np.float32(x_left)
+        x_right = np.float32(x_right)
         # These are independent variable arrays for interpolation.
         x_lower_uniform = np.arange(nx, dtype=np.float32)
         x_upper_uniform = np.arange(nx, dtype=np.float32)
@@ -1773,7 +1777,10 @@ def fuvFlagOutOfBounds(hdr, dq_array, info, switches,
                                           ny, dy)
         (x_left, x_right)  = applyOffsets(x_left_interp, x_right_interp,
                                           nx, dx, x_offset)
-
+        y_lower = np.float32(y_lower)
+        y_upper = np.float32(y_upper)
+        x_left = np.float32(x_left)
+        x_right = np.float32(x_right)
         ccos.clear_rows(temp, y_lower, y_upper, x_left, x_right)
     elif save_sub is not None:
         (x0, x1, y0, y1) = save_sub
@@ -1980,7 +1987,6 @@ def flagOutsideActiveArea(dq_array, segment, brftab, x_offset,
 
     b_left += x_offset
     b_right += x_offset
-
     (ny, nx) = dq_array.shape
 
     if b_low >= 0:

calcos/dispersion.py

@@ -134,7 +134,7 @@ def evalDisp(self, x):
             Wavelength (or array of wavelengths) at x
         """
 
-        x_prime = x + self.delta
+        x_prime = np.float64(x) + np.float64(self.delta)
 
         sum = self.coeff[self.ncoeff-1]
         for i in range(self.ncoeff-2, -1, -1):
@@ -159,7 +159,7 @@ def evalDerivDisp(self, x):
             Slope at x, in Angstroms per pixel
         """
 
-        x_prime = x + self.delta
+        x_prime = np.float64(x) + np.float64(self.delta)
 
         sum = (self.ncoeff - 1.) * self.coeff[self.ncoeff-1]
         for n in range(self.ncoeff-2, 0, -1):

calcos/timeline.py

@@ -189,7 +189,7 @@ def createTimeline(input, fd, info, reffiles,
         rv_col = tl_data.field("radial_vel")
 
     for i in range(len(tl_time)):
-        mjd = tl_time[i] / SEC_PER_DAY + info["expstart"]
+        mjd = np.float64(tl_time[i]) / SEC_PER_DAY + info["expstart"]
         (rect_hst, vel_hst) = orb.getPos(mjd)
         (r, ra_hst, dec_hst) = rectToSph(rect_hst)
         # Assume that we want geocentric latitude.  The difference from

calcos/timetag.py

@@ -180,7 +180,7 @@ def timetagBasicCalibration(input, inpha, outtag,
     (stim_param, stim_countrate, stim_livetime) = initTempcorr(events,
             input, info, switches, reffiles, headers[1],
             cl_args["stimfile"])
-
+    
     doTempcorr(stim_param, events, info, switches, reffiles, phdr)
 
     doGeocorr(events, info, switches, reffiles, phdr)
@@ -1624,7 +1624,6 @@ def computeThermalParam(time, x, y, dq,
         else:
             rms_s2[0] = math.sqrt(sumstim[8])
             rms_s2[1] = math.sqrt(sumstim[9])
-
     if total_counts1 > 0 and total_counts2 > 0:
         stim_countrate = (total_counts1 + total_counts2) / (2. * exptime)
     elif total_counts1 > 0:
@@ -1771,14 +1770,14 @@ def updateStimSum(sumstim, nevents1, s1, sumsq1, found_s1,
      n2, sum2y, sum2x, sumsq2y, sumsq2x) = sumstim
 
     if found_s1:
-        n1 = n1 + nevents1
+        n1 = n1 + np.float64(nevents1)
         sum1y = sum1y + s1[0] * nevents1
         sum1x = sum1x + s1[1] * nevents1
         sumsq1y = sumsq1y + sumsq1[0]
         sumsq1x = sumsq1x + sumsq1[1]
 
     if found_s2:
-        n2 = n2 + nevents2
+        n2 = n2 + np.float64(nevents2)
         sum2y = sum2y + s2[0] * nevents2
         sum2x = sum2x + s2[1] * nevents2
         sumsq2y = sumsq2y + sumsq2[0]
@@ -2588,8 +2587,8 @@ def  blurDQ(trace_dq, minmax_shift_dict, minmax_doppler, doppler_boundary, widen
                 xshifts.append(int(round(max_shift1 + maxdopp + widen)))
                 xshifts.append(int(round(min_shift1 + mindopp - widen)))
             else:
-                xshifts.append(int(round(max_shift1 + widen)))
-                xshifts.append(int(round(min_shift1 - widen)))
+                xshifts.append(int(round(np.float64(max_shift1) + widen)))
+                xshifts.append(int(round(np.float64(min_shift1) - widen)))
             for xshift in range(min(xshifts), max(xshifts)+1):
                 cosutil.printMsg("Shifting to %d, %d" % (xshift, yshift))
                 shifted_dq = arrayShift(y_shifted_dq[int(lower_y):int(upper_y)], 0, xshift,
@@ -2974,7 +2973,7 @@ def initHelcorr(events, info, hdr):
     # get midpoint of exposure, MJD
     expstart = info["expstart"]
     time = events.field("time")
-    t_mid = expstart + (time[0] + time[len(time)-1]) / 2. / SEC_PER_DAY
+    t_mid = expstart + (np.float64(time[0]) + np.float64(time[len(time)-1])) / 2. / SEC_PER_DAY
 
     # Compute radial velocity and heliocentric correction factor (the latter
     # is actually not used here).
@@ -3995,7 +3994,7 @@ def writeImages(x, y, epsilon, dq,
 
     # Use the Frequentist variance function.
     err_lower, err_upper = cosutil.errFrequentist(C_counts)
-    errC_rate = err_upper / exptime
+    errC_rate = np.float32(err_upper / exptime)
 
     if outcounts is not None:
         C_rate = C_counts / exptime
@@ -4049,9 +4048,14 @@ def makeImage(outimage, phdr, headers, sci_array, err_array, dq_array):
     fd = fits.HDUList(primary_hdu)
     fd[0].header["nextend"] = 3
     cosutil.updateFilename(fd[0].header, outimage)
-
-    makeImageHDU(fd, headers[1], sci_array, name="SCI")
-    makeImageHDU(fd, headers[2], err_array, name="ERR")
+    if sci_array is not None:
+        makeImageHDU(fd, headers[1], np.float32(sci_array), name="SCI")
+    else:
+        makeImageHDU(fd, headers[1], sci_array, name="SCI")
+    if err_array is not None:
+        makeImageHDU(fd, headers[2], np.float32(err_array), name="ERR")
+    else:
+        makeImageHDU(fd, headers[2], err_array, name="ERR")
     makeImageHDU(fd, headers[3], dq_array, name="DQ")
 
     fd.writeto(outimage, output_verify='silentfix')

pyproject.toml

@@ -16,7 +16,6 @@ classifiers = [
 ]
 dependencies = [
     "astropy>=5.0.4",
-    "numpy<2.0",
     "scipy",
     "stsci.tools>=4.0.0",
 ]
@@ -46,7 +45,7 @@ requires = [
     "setuptools>=42.0",
     "setuptools_scm[toml]>=3.4",
     "wheel",
-    "oldest-supported-numpy",
+    "numpy>=2.0.0",
 ]
 build-backend = "setuptools.build_meta"
 

src/ccos.c

@@ -108,6 +108,7 @@ bin2d bins a 2-D image to a smaller 2-D image (block sum).
 2017 Jan 30     Add bilinear_interpolation function for walk correction
 */
 
+# define PY_SSIZE_T_CLEAN
 # include <Python.h>
 
 # include <stdlib.h>
@@ -356,7 +357,6 @@ static PyObject *ccos_binevents(PyObject *self, PyObject *args) {
 	    PyErr_SetString(PyExc_RuntimeError, "can't read arguments");
 	    return NULL;
 	}
-
 	if (PyArray_TYPE(ox) == NPY_INT16) {
 	    x = (PyArrayObject *)PyArray_FROM_OTF(ox, NPY_INT16,
 		NPY_IN_ARRAY);
@@ -4168,7 +4168,6 @@ static PyMethodDef ccos_methods[] = {
 	{NULL, NULL, 0, NULL}
 };
 
-#if defined(NPY_PY3K)
 static struct PyModuleDef moduledef = {
     PyModuleDef_HEAD_INIT,
     "ccos",
@@ -4180,33 +4179,20 @@ static struct PyModuleDef moduledef = {
     NULL,
     NULL
 };
-#endif
 
-#if defined(NPY_PY3K)
 PyObject *PyInit_ccos(void) 
-#else
-PyMODINIT_FUNC initccos(void)
-#endif
+
 {
 	PyObject *mod;		/* the module */
 	PyObject *dict;		/* the module's dictionary */
 
-#if defined(NPY_PY3K)
 	mod = PyModule_Create(&moduledef);
-#else
-	mod = Py_InitModule("ccos", ccos_methods);
-#endif
 	import_array();
 
 	/* set the doc string */
 	dict = PyModule_GetDict(mod);
-#if defined(NPY_PY3K)
+
 	PyDict_SetItemString(dict, "__doc__",
 		PyUnicode_FromString(DocString()));
         return mod;
-#else
-	PyDict_SetItemString(dict, "__doc__",
-		PyString_FromString(DocString()));
-	return;
-#endif
 }