Result of Image.calculateFWHM() depends on NumPy version

[welucas2]

Image.calculateFWHM() can return different results depending on the version of NumPy when calling np.argsort() at line

GalSim/galsim/image.py

Line 1587 in 4c04a1e

indx = np.argsort(rsq.ravel())

NumPy 2.0 and more recent has an argument stable to argsort() which guarantees that it will return in order the indices of equal-valued elements of the array. This is by default not enabled, and before 2.0 the option doesn't exist -- in this case, it seems the indices of equal valued elements can be given in any order that seems to be determined by the version of NumPy.

GalSim's FWHM calculation orders the elements by distance from the image center and finds the first element with a surface brightness < half the max. It then interpolates with the brightness of the previous element, which can have been sorted in any order with the other elements the same distance from the center. From tests it does seem that a given NumPy version will always sort the array the same way, but it can and does change when using different versions. Here's a very small reproducer showing different returned values:

import galsim
gal = galsim.Gaussian(flux=1.e5, sigma=5.)
rng = galsim.BaseDeviate(1234)
img = gal.drawImage(method='phot', rng=rng, nx=50, ny=50)
print(img.calculateFWHM())

With NumPy 1.26.4 this returns a FWHM of 13.870431838154461.
With NumPy 2.0.2 it instead returns 13.828419784303065.
With NumPy 2.1.0 it returns 13.870431838154461 once more.

This is likely behind an odd failing test in imSim.