From 8d377809bd745cad84b5e9e8b6d7b965d1704075 Mon Sep 17 00:00:00 2001 From: Mark Farnum Date: Wed, 8 Jan 2025 11:53:12 -0500 Subject: [PATCH 1/3] Fix issue with floating point imprecision & integer overflow --- romanisim/image.py | 13 ++++++++++++- 1 file changed, 12 insertions(+), 1 deletion(-) diff --git a/romanisim/image.py b/romanisim/image.py index 0f84e40..a3d992b 100644 --- a/romanisim/image.py +++ b/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(np.float32(2**31 - 1), 0) 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 From 6c0b0ddb4a9bfdf410727eabe5c1395b90e25fbe Mon Sep 17 00:00:00 2001 From: Mark Farnum Date: Wed, 8 Jan 2025 13:40:08 -0500 Subject: [PATCH 2/3] Add test --- romanisim/tests/test_image.py | 13 +++++++++++++ 1 file changed, 13 insertions(+) diff --git a/romanisim/tests/test_image.py b/romanisim/tests/test_image.py index 63f7adf..db96da8 100644 --- a/romanisim/tests/test_image.py +++ b/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() From 1e2c0ed9cdab439bc15b314c9f7e26ce75d3c9bc Mon Sep 17 00:00:00 2001 From: Mark Farnum Date: Tue, 14 Jan 2025 15:15:17 -0500 Subject: [PATCH 3/3] Fix issue with old version of numpy --- romanisim/image.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/romanisim/image.py b/romanisim/image.py index a45e72f..7f678cf 100644 --- a/romanisim/image.py +++ b/romanisim/image.py @@ -478,7 +478,7 @@ def simulate_counts_generic(image, exptime, objlist=None, psf=None, # 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(np.float32(2**31 - 1), 0) + MAX_SAFE_VALUE = np.nextafter(2**31 - 1, 0, dtype=np.float32) image.array[:, :] = rng_numpy.binomial( np.clip(image.array, 0, MAX_SAFE_VALUE).astype("i4"), flat / maxflat )