Fix floating point issues (#194)

* Fix issue with floating point imprecision & integer overflow

* Add test

* Fix issue with old version of numpy

---------

Co-authored-by: Eddie Schlafly <[email protected]>

romanisim/image.py

@@ -469,8 +469,19 @@ def simulate_counts_generic(image, exptime, objlist=None, psf=None,
         image.quantize()
         rng_numpy_seed = rng.raw()
         rng_numpy = np.random.default_rng(rng_numpy_seed)
+
+        # NOTE: Here, we convert the array of float32 values to int32 (i4), so
+        # we need to "clip" them to be in the valid range first. However, just
+        # clipping with the maximum int32 value (2^31 - 1) isn't sufficient.
+        # This is because `np.float32(2**31 - 1)` is the same value as
+        # `np.float32(2**31)` due to floating point imprecision, and on some
+        # architectures, converting that value to int32 rolls over to a negative
+        # number. To resolve, we use `np.nextafter` to get the previous floating
+        # point number, which is roughly 2^31 - 128.
+        MAX_SAFE_VALUE = np.nextafter(2**31 - 1, 0, dtype=np.float32)
         image.array[:, :] = rng_numpy.binomial(
-            np.clip(image.array, 0, 2**31 - 1).astype('i4'), flat / maxflat)
+            np.clip(image.array, 0, MAX_SAFE_VALUE).astype("i4"), flat / maxflat
+        )
 
     if dark is not None:
         workim = image * 0

romanisim/tests/test_image.py

@@ -411,6 +411,19 @@ def test_simulate_counts_generic():
     assert np.sum(objinfo['counts'] > 0) == 0
     # these sources should be out of bounds
 
+    ### Assert that out-of-range values do not cause integer overflow:
+    try:
+        # Have numpy raise on warnings:
+        np.seterr(all='raise')
+        # Test simulating with a value that is out of range for an int32
+        im9 = im.copy()
+        im9.array[0][0] = np.float32(2**40)
+        image.simulate_counts_generic(im9, exptime, sky=sky, flat=0.5, zpflux=zpflux)
+    finally:
+        # revert numpy warning settings to default
+        np.seterr(all='print')
+
+
 
 def test_simulate_counts():
     imdict = set_up_image_rendering_things()