Skip to content

test: changes to run more of moffat tests #4

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
beckermr merged 1 commit into from
Apr 30, 2026
Merged

test: changes to run more of moffat tests #4

Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
216 changes: 122 additions & 94 deletions tests/test_moffat.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -26,48 +26,56 @@
# Directory containing the reference images.
imgdir = os.path.join(path, "SBProfile_comparison_images")

if is_jax_galsim():
USE_TRUNC = 0.0
else:
USE_TRUNC = 1.0


@timer
def test_moffat():
"""Test the generation of a specific Moffat profile against a known result.
"""
fwhm_backwards_compatible = 1.3178976627539716
savedImg = galsim.fits.read(os.path.join(imgdir, "moffat_2_5.fits"))
dx = 0.2
myImg = galsim.ImageF(savedImg.bounds, scale=dx)
myImg.setCenter(0,0)

# Code was formerly:
# moffat = galsim.Moffat(beta=2, truncationFWHM=5, flux=1, half_light_radius=1)
#
# ...but this is no longer quite so simple since we changed the handling of trunc to be in
# physical units. However, the same profile can be constructed using
# fwhm=1.3178976627539716
# as calculated by interval bisection in devutils/external/calculate_moffat_radii.py
fwhm_backwards_compatible = 1.3178976627539716
moffat = galsim.Moffat(beta=2, half_light_radius=1, trunc=5*fwhm_backwards_compatible, flux=1)
moffat.drawImage(myImg, method="sb", use_true_center=False)
np.testing.assert_array_almost_equal(
myImg.array, savedImg.array, 5,
err_msg="Using GSObject Moffat disagrees with expected result")

# Check with default_params
moffat = galsim.Moffat(beta=2, half_light_radius=1, trunc=5*fwhm_backwards_compatible, flux=1,
gsparams=default_params)
moffat.drawImage(myImg, method="sb", use_true_center=False)
np.testing.assert_array_almost_equal(
myImg.array, savedImg.array, 5,
err_msg="Using GSObject Moffat with default_params disagrees with expected result")
moffat = galsim.Moffat(beta=2, half_light_radius=1, trunc=5*fwhm_backwards_compatible, flux=1,
gsparams=galsim.GSParams())
moffat.drawImage(myImg, method="sb", use_true_center=False)
np.testing.assert_array_almost_equal(
myImg.array, savedImg.array, 5,
err_msg="Using GSObject Moffat with GSParams() disagrees with expected result")
if is_jax_galsim():
pass
else:
# Code was formerly:
# moffat = galsim.Moffat(beta=2, truncationFWHM=5, flux=1, half_light_radius=1)
#
# ...but this is no longer quite so simple since we changed the handling of trunc to be in
# physical units. However, the same profile can be constructed using
# fwhm=1.3178976627539716
# as calculated by interval bisection in devutils/external/calculate_moffat_radii.py
moffat = galsim.Moffat(beta=2, half_light_radius=1, trunc=5*fwhm_backwards_compatible, flux=1)
moffat.drawImage(myImg, method="sb", use_true_center=False)
np.testing.assert_array_almost_equal(
myImg.array, savedImg.array, 5,
err_msg="Using GSObject Moffat disagrees with expected result")

# Check with default_params
moffat = galsim.Moffat(beta=2, half_light_radius=1, trunc=5*fwhm_backwards_compatible, flux=1,
gsparams=default_params)
moffat.drawImage(myImg, method="sb", use_true_center=False)
np.testing.assert_array_almost_equal(
myImg.array, savedImg.array, 5,
err_msg="Using GSObject Moffat with default_params disagrees with expected result")
moffat = galsim.Moffat(beta=2, half_light_radius=1, trunc=5*fwhm_backwards_compatible, flux=1,
gsparams=galsim.GSParams())
moffat.drawImage(myImg, method="sb", use_true_center=False)
np.testing.assert_array_almost_equal(
myImg.array, savedImg.array, 5,
err_msg="Using GSObject Moffat with GSParams() disagrees with expected result")

# Use non-unity values.
moffat = galsim.Moffat(beta=3.7, flux=1.7, half_light_radius=2.3, trunc=8.2)
moffat = galsim.Moffat(beta=3.7, flux=1.7, half_light_radius=2.3, trunc=8.2 * USE_TRUNC)
gsp = galsim.GSParams(xvalue_accuracy=1.e-8, kvalue_accuracy=1.e-8)
moffat2 = galsim.Moffat(beta=3.7, flux=1.7, half_light_radius=2.3, trunc=8.2, gsparams=gsp)
moffat2 = galsim.Moffat(beta=3.7, flux=1.7, half_light_radius=2.3, trunc=8.2 * USE_TRUNC, gsparams=gsp)
assert moffat2 != moffat
assert moffat2 == moffat.withGSParams(gsp)
assert moffat2 == moffat.withGSParams(xvalue_accuracy=1.e-8, kvalue_accuracy=1.e-8)
Expand Down Expand Up @@ -106,9 +114,12 @@ def test_moffat():
assert_raises(TypeError, galsim.Moffat, beta=3, scale_radius=3, half_light_radius=1)
assert_raises(TypeError, galsim.Moffat, beta=3)

# beta <= 1.1 needs to be truncated.
assert_raises(ValueError, galsim.Moffat, beta=1.1, scale_radius=3)
assert_raises(ValueError, galsim.Moffat, beta=0.9, scale_radius=3)
if is_jax_galsim():
pass
else:
# beta <= 1.1 needs to be truncated.
assert_raises(ValueError, galsim.Moffat, beta=1.1, scale_radius=3)
assert_raises(ValueError, galsim.Moffat, beta=0.9, scale_radius=3)

# trunc must be > sqrt(2) * hlr
assert_raises(ValueError, galsim.Moffat, beta=3, half_light_radius=1, trunc=1.4)
Expand All @@ -121,61 +132,67 @@ def test_moffat():
def test_moffat_properties():
"""Test some basic properties of the Moffat profile.
"""
# Code was formerly:
# psf = galsim.Moffat(beta=2.0, truncationFWHM=2, flux=test_flux, half_light_radius=1)
#
# ...but this is no longer quite so simple since we changed the handling of trunc to be in
# physical units. However, the same profile can be constructed using
# fwhm=1.4686232496771867,
# as calculated by interval bisection in devutils/external/calculate_moffat_radii.py
test_flux = 1.8
fwhm_backwards_compatible = 1.4686232496771867
psf = galsim.Moffat(beta=2.0, fwhm=fwhm_backwards_compatible,
trunc=2*fwhm_backwards_compatible, flux=test_flux)
# Check that we are centered on (0, 0)
cen = galsim.PositionD(0, 0)
np.testing.assert_equal(psf.centroid, cen)
# Check Fourier properties
if is_jax_galsim():
np.testing.assert_allclose(psf.maxk, 11.634597424960159, atol=0, rtol=0.2)
test_flux = 1.8
fwhm_backwards_compatible = 1.4686232496771867
psf = galsim.Moffat(beta=2.0, half_light_radius=1.,
trunc=2*fwhm_backwards_compatible * USE_TRUNC, flux=test_flux)
else:
np.testing.assert_array_almost_equal(psf.maxk, 11.634597424960159)
np.testing.assert_array_almost_equal(psf.stepk, 0.62831853071795873)
np.testing.assert_array_almost_equal(psf.kValue(cen), test_flux+0j)
np.testing.assert_array_almost_equal(psf.half_light_radius, 1.0)
np.testing.assert_array_almost_equal(psf.fwhm, fwhm_backwards_compatible)
np.testing.assert_array_almost_equal(psf.xValue(cen), 0.50654651638242509)
np.testing.assert_array_almost_equal(psf.kValue(cen), (1+0j) * test_flux)
np.testing.assert_array_almost_equal(psf.flux, test_flux)
np.testing.assert_array_almost_equal(psf.xValue(cen), psf.max_sb)

# Now create the same profile using the half_light_radius:
psf = galsim.Moffat(beta=2.0, half_light_radius=1.,
trunc=2*fwhm_backwards_compatible, flux=test_flux)
np.testing.assert_equal(psf.centroid, cen)
if is_jax_galsim():
np.testing.assert_allclose(psf.maxk, 11.634597424960159, atol=0, rtol=0.2)
else:
np.testing.assert_array_almost_equal(psf.maxk, 11.634597424960159)
np.testing.assert_array_almost_equal(psf.stepk, 0.62831853071795862)
np.testing.assert_array_almost_equal(psf.kValue(cen), test_flux+0j)
np.testing.assert_array_almost_equal(psf.half_light_radius, 1.0)
np.testing.assert_array_almost_equal(psf.fwhm, fwhm_backwards_compatible)
np.testing.assert_array_almost_equal(psf.xValue(cen), 0.50654651638242509)
np.testing.assert_array_almost_equal(psf.kValue(cen), (1+0j) * test_flux)
np.testing.assert_array_almost_equal(psf.flux, test_flux)
np.testing.assert_array_almost_equal(psf.xValue(cen), psf.max_sb)
# Code was formerly:
# psf = galsim.Moffat(beta=2.0, truncationFWHM=2, flux=test_flux, half_light_radius=1)
#
# ...but this is no longer quite so simple since we changed the handling of trunc to be in
# physical units. However, the same profile can be constructed using
# fwhm=1.4686232496771867,
# as calculated by interval bisection in devutils/external/calculate_moffat_radii.py
test_flux = 1.8
fwhm_backwards_compatible = 1.4686232496771867
psf = galsim.Moffat(beta=2.0, fwhm=fwhm_backwards_compatible,
trunc=2*fwhm_backwards_compatible, flux=test_flux)
# Check that we are centered on (0, 0)
cen = galsim.PositionD(0, 0)
np.testing.assert_equal(psf.centroid, cen)
# Check Fourier properties
if is_jax_galsim():
np.testing.assert_allclose(psf.maxk, 11.634597424960159, atol=0, rtol=0.2)
else:
np.testing.assert_array_almost_equal(psf.maxk, 11.634597424960159)
np.testing.assert_array_almost_equal(psf.stepk, 0.62831853071795873)
np.testing.assert_array_almost_equal(psf.kValue(cen), test_flux+0j)
np.testing.assert_array_almost_equal(psf.half_light_radius, 1.0)
np.testing.assert_array_almost_equal(psf.fwhm, fwhm_backwards_compatible)
np.testing.assert_array_almost_equal(psf.xValue(cen), 0.50654651638242509)
np.testing.assert_array_almost_equal(psf.kValue(cen), (1+0j) * test_flux)
np.testing.assert_array_almost_equal(psf.flux, test_flux)
np.testing.assert_array_almost_equal(psf.xValue(cen), psf.max_sb)

# Now create the same profile using the half_light_radius:
psf = galsim.Moffat(beta=2.0, half_light_radius=1.,
trunc=2*fwhm_backwards_compatible, flux=test_flux)
np.testing.assert_equal(psf.centroid, cen)
if is_jax_galsim():
np.testing.assert_allclose(psf.maxk, 11.634597424960159, atol=0, rtol=0.2)
else:
np.testing.assert_array_almost_equal(psf.maxk, 11.634597424960159)
np.testing.assert_array_almost_equal(psf.stepk, 0.62831853071795862)
np.testing.assert_array_almost_equal(psf.kValue(cen), test_flux+0j)
np.testing.assert_array_almost_equal(psf.half_light_radius, 1.0)
np.testing.assert_array_almost_equal(psf.fwhm, fwhm_backwards_compatible)
np.testing.assert_array_almost_equal(psf.xValue(cen), 0.50654651638242509)
np.testing.assert_array_almost_equal(psf.kValue(cen), (1+0j) * test_flux)
np.testing.assert_array_almost_equal(psf.flux, test_flux)
np.testing.assert_array_almost_equal(psf.xValue(cen), psf.max_sb)

# Check input flux vs output flux
for inFlux in np.logspace(-2, 2, 10):
psfFlux = galsim.Moffat(2.0, fwhm=fwhm_backwards_compatible,
trunc=2*fwhm_backwards_compatible, flux=inFlux)
trunc=2*fwhm_backwards_compatible * USE_TRUNC, flux=inFlux)
outFlux = psfFlux.flux
np.testing.assert_array_almost_equal(outFlux, inFlux)

# Check that stepk and maxk scale correctly with radius
psf2 = galsim.Moffat(beta=2.0, half_light_radius=5.,
trunc=5*2*fwhm_backwards_compatible, flux=test_flux)
trunc=5*2*fwhm_backwards_compatible * USE_TRUNC, flux=test_flux)
np.testing.assert_almost_equal(psf2.maxk, psf.maxk/5)
np.testing.assert_almost_equal(psf2.stepk, psf.stepk/5)

Expand Down Expand Up @@ -205,19 +222,24 @@ def test_moffat_maxk():
galsim.Moffat(beta=1.22, scale_radius=23, flux=23),
galsim.Moffat(beta=3.6, scale_radius=2, flux=23),
galsim.Moffat(beta=12.9, scale_radius=5, flux=23),
galsim.Moffat(beta=1.22, scale_radius=7, flux=23, trunc=30),
galsim.Moffat(beta=3.6, scale_radius=9, flux=23, trunc=50),
galsim.Moffat(beta=12.9, scale_radius=11, flux=23, trunc=1000),
galsim.Moffat(beta=1.22, scale_radius=7, flux=23, trunc=30 * USE_TRUNC),
galsim.Moffat(beta=3.6, scale_radius=9, flux=23, trunc=50 * USE_TRUNC),
galsim.Moffat(beta=12.9, scale_radius=11, flux=23, trunc=1000 * USE_TRUNC),
]
threshs = [1.e-3, 1.e-4, 0.03]
print('beta \t trunc \t thresh \t kValue(maxk) \t maxk')
print('\nbeta \t trunc \t thresh \t kValue(maxk) \t maxk')
for psf in psfs:
for thresh in threshs:
psf = psf.withGSParams(maxk_threshold=thresh)
rtol = 1.e-7 if psf.trunc == 0 else 3.e-3
if is_jax_galsim():
atol = 5e-6
rtol = 0
else:
atol = 0
fk = psf.kValue(psf.maxk,0).real/psf.flux
print(f'{psf.beta} \t {int(psf.trunc)} \t {thresh:.1e} \t {fk:.3e} \t {psf.maxk:.3e}')
np.testing.assert_allclose(abs(psf.kValue(psf.maxk,0).real)/psf.flux, thresh, rtol=rtol)
np.testing.assert_allclose(abs(psf.kValue(psf.maxk,0).real)/psf.flux, thresh, rtol=rtol, atol=atol)


@timer
Expand Down Expand Up @@ -326,7 +348,7 @@ def test_moffat_radii():

# Test constructor using half-light-radius:
test_gal = galsim.Moffat(flux=1., beta=test_beta, half_light_radius=test_hlr,
trunc=2*test_hlr)
trunc=2*test_hlr * USE_TRUNC)
hlr_sum = radial_integrate(test_gal, 0., test_hlr)
print('hlr_sum = ',hlr_sum)
np.testing.assert_almost_equal(
Expand Down Expand Up @@ -356,12 +378,15 @@ def test_moffat_radii():
np.testing.assert_equal(
test_gal.scale_radius, got_scale,
err_msg="Moffat scale_radius returned wrong value")
np.testing.assert_equal(
test_gal.trunc, 2*test_hlr,
err_msg="Moffat trunc returned wrong value")
if is_jax_galsim():
pass
else:
np.testing.assert_equal(
test_gal.trunc, 2*test_hlr,
err_msg="Moffat trunc returned wrong value")

# Test constructor using scale radius:
test_gal = galsim.Moffat(flux=1., beta=test_beta, trunc=2*test_scale,
test_gal = galsim.Moffat(flux=1., beta=test_beta, trunc=2*test_scale * USE_TRUNC,
scale_radius=test_scale)
center = test_gal.xValue(galsim.PositionD(0,0))
ratio = test_gal.xValue(galsim.PositionD(test_scale,0)) / center
Expand Down Expand Up @@ -389,7 +414,7 @@ def test_moffat_radii():
err_msg="Error in FWHM for truncated Moffat initialized with scale radius")

# Test constructor using FWHM:
test_gal = galsim.Moffat(flux=1, beta=test_beta, trunc=2.*test_fwhm,
test_gal = galsim.Moffat(flux=1, beta=test_beta, trunc=2.*test_fwhm * USE_TRUNC,
fwhm = test_fwhm)
center = test_gal.xValue(galsim.PositionD(0,0))
ratio = test_gal.xValue(galsim.PositionD(test_fwhm/2.,0)) / center
Expand Down Expand Up @@ -436,19 +461,19 @@ def test_moffat_flux_scaling():
for test_trunc in [ 0., 8.5 ]:

# init with scale_radius only (should be ok given last tests)
obj = galsim.Moffat(scale_radius=test_scale, beta=test_beta, trunc=test_trunc,
obj = galsim.Moffat(scale_radius=test_scale, beta=test_beta, trunc=test_trunc * USE_TRUNC,
flux=test_flux)
obj *= 2.
np.testing.assert_almost_equal(
obj.flux, test_flux * 2., decimal=param_decimal,
err_msg="Flux param inconsistent after __imul__.")
obj = galsim.Moffat(scale_radius=test_scale, beta=test_beta, trunc=test_trunc,
obj = galsim.Moffat(scale_radius=test_scale, beta=test_beta, trunc=test_trunc * USE_TRUNC,
flux=test_flux)
obj /= 2.
np.testing.assert_almost_equal(
obj.flux, test_flux / 2., decimal=param_decimal,
err_msg="Flux param inconsistent after __idiv__.")
obj = galsim.Moffat(scale_radius=test_scale, beta=test_beta, trunc=test_trunc,
obj = galsim.Moffat(scale_radius=test_scale, beta=test_beta, trunc=test_trunc * USE_TRUNC,
flux=test_flux)
obj2 = obj * 2.
# First test that original obj is unharmed...
Expand All @@ -459,7 +484,7 @@ def test_moffat_flux_scaling():
np.testing.assert_almost_equal(
obj2.flux, test_flux * 2., decimal=param_decimal,
err_msg="Flux param inconsistent after __rmul__ (result).")
obj = galsim.Moffat(scale_radius=test_scale, beta=test_beta, trunc=test_trunc,
obj = galsim.Moffat(scale_radius=test_scale, beta=test_beta, trunc=test_trunc * USE_TRUNC,
flux=test_flux)
obj2 = 2. * obj
# First test that original obj is unharmed...
Expand All @@ -470,7 +495,7 @@ def test_moffat_flux_scaling():
np.testing.assert_almost_equal(
obj2.flux, test_flux * 2., decimal=param_decimal,
err_msg="Flux param inconsistent after __mul__ (result).")
obj = galsim.Moffat(scale_radius=test_scale, beta=test_beta, trunc=test_trunc,
obj = galsim.Moffat(scale_radius=test_scale, beta=test_beta, trunc=test_trunc * USE_TRUNC,
flux=test_flux)
obj2 = obj / 2.
# First test that original obj is unharmed...
Expand Down Expand Up @@ -522,7 +547,7 @@ def test_moffat_shoot():


@timer
def test_ne():
def test_moffat_ne():
"""Test base.py GSObjects for not-equals."""
# Define some universal gsps
gsp = galsim.GSParams(maxk_threshold=1.1e-3, folding_threshold=5.1e-3)
Expand All @@ -534,9 +559,12 @@ def test_ne():
galsim.Moffat(beta=3.0, scale_radius=1.1),
galsim.Moffat(beta=3.0, half_light_radius=1.2),
galsim.Moffat(beta=3.0, fwhm=1.3),
galsim.Moffat(beta=3.0, scale_radius=1.1, trunc=2.0),
galsim.Moffat(beta=3.0, scale_radius=1.0, flux=1.1),
galsim.Moffat(beta=3.0, scale_radius=1.0, gsparams=gsp)]
if is_jax_galsim():
pass
else:
gals += [galsim.Moffat(beta=3.0, scale_radius=1.1, trunc=2.0)]
check_all_diff(gals)


Expand Down
Loading