From 16cf08085804bec88848af3d31aea19b91a70887 Mon Sep 17 00:00:00 2001 From: Mike Jarvis Date: Fri, 3 Oct 2025 18:00:49 -0400 Subject: [PATCH 1/2] Update calculateFWHM to be insensitive to np.argsort algorithmic details (#1336) --- galsim/image.py | 9 +++++++-- tests/test_calc.py | 25 +++++++++++++++++++++++++ 2 files changed, 32 insertions(+), 2 deletions(-) diff --git a/galsim/image.py b/galsim/image.py index 747c97a46b..49c01871fe 100644 --- a/galsim/image.py +++ b/galsim/image.py @@ -1590,13 +1590,18 @@ def calculateFWHM(self, center=None, Imax=0.): # Find the first value with I < 0.5 * Imax k = np.argmax(data < 0.5 * Imax) - Ik = data[k] / Imax + # If there are duplicate rsq values, argmax won't be deterministic across numpy versions. + # To achieve consistent results, we take the average of all data values with the same rsq. + k2 = np.searchsorted(rsqf, rsqf[k], side='right') + Ik = np.mean(data[k:k2]) / Imax # Interpolate (linearly) between this and the previous value. if k == 0: rsqhm = rsqf[0] * (0.5 / Ik) else: - Ikm1 = data[k-1] / Imax + # Similarly, use the mean data for all rsq equal to rsqf[k-1]. + k1 = np.searchsorted(rsqf, rsqf[k-1], side='left') + Ikm1 = np.mean(data[k1:k]) / Imax rsqhm = (rsqf[k-1] * (Ik-0.5) + rsqf[k] * (0.5-Ikm1)) / (Ik-Ikm1) # This has all been done in pixels. So normalize according to the pixel scale. diff --git a/tests/test_calc.py b/tests/test_calc.py index a80e82c71f..655e8b683f 100644 --- a/tests/test_calc.py +++ b/tests/test_calc.py @@ -369,5 +369,30 @@ def test_fwhm(): err_msg="non-square image.calculateFWHM is not accurate.") +def test_fwhm_deterministic_fwhm(): + # This test was proposed by Wiliam Lucas in issue #1336. + # Numpy changed some details about how argsort works, so different numpy versions had + # given different answers to the following calculation. + + gal = galsim.Gaussian(flux=1.e5, sigma=5.) + rng = galsim.BaseDeviate(1234) + img = gal.drawImage(method='phot', rng=rng, nx=50, ny=50) + fwhm = img.calculateFWHM() + # This value used to be different on systems with different numpy versions. + # Note: this only affects situations where the center is a pixel center or pixel corner + # (or some other very special locations) such that the rsq values include repeated values. + # Then argsort could return a different data index for what pixel was the first one past + # the half max value. For symmetric profiles, even this probalby wasn't a problem, since all + # data values at that radius have the same value. But for photon shooting, there is photon + # noise, which means the values are not identical, leading to indeterminacy in the FWHM value. + print(fwhm) + + # The code was changed in GalSim version 2.7 to be insensitive to the details of what argsort + # does with repeated values. Specifically, we now take the mean of all data values with the + # same rsq value. + # Thus, this value should be consistent across any algorithmic changes numpy might make to + # how argsort works in detail. + assert np.isclose(fwhm, 13.715780373024351) + if __name__ == "__main__": runtests(__file__) From d23bcd9adc43d00c9d303a58fcd5683f1d112939 Mon Sep 17 00:00:00 2001 From: Mike Jarvis Date: Sun, 7 Sep 2025 16:48:21 -0400 Subject: [PATCH 2/2] Fix for numpy 2.3 --- galsim/image.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/galsim/image.py b/galsim/image.py index 49c01871fe..a6d20b95bc 100644 --- a/galsim/image.py +++ b/galsim/image.py @@ -1794,7 +1794,10 @@ def __pow__(self, other): def __ipow__(self, other): if not isinstance(other, int) and not isinstance(other, float): raise TypeError("Can only raise an image to a float or int power!") - self.array[:,:] **= other + if not self.isinteger or isinstance(other, int): + self.array[:,:] **= other + else: + self.array[:,:] = self._safe_cast(self.array ** other) return self def __neg__(self):