Fix #1324 and add more regression tests of maxk, stepk, etc.

(#1324)

galsim/spergel.py

@@ -205,7 +205,7 @@ def _maxk(self):
 
     @property
     def _stepk(self):
-        R = self.calculateFluxRadius(1.0 - self.gsparams.folding_threshold) * self._r0
+        R = self.calculateFluxRadius(1.0 - self.gsparams.folding_threshold)
         # Go to at least 5*hlr
         R = max(R, self.gsparams.stepk_minimum_hlr * self.half_light_radius)
         return math.pi / R

tests/test_airy.py

@@ -222,6 +222,39 @@ def test_airy_flux_scaling():
         obj2.flux, test_flux / 2., decimal=param_decimal,
         err_msg="Flux param inconsistent after __div__ (result).")
 
+@timer
+def test_airy_properties():
+    """Test some basic properties of the Airy profile.
+    """
+    # Regression test based on v2.3.5 version of the code.
+
+    test_loD = 1.9
+    test_obscuration = 0.32
+    test_flux = 17.9
+    psf = galsim.Airy(lam_over_diam=test_loD, flux=test_flux, obscuration=test_obscuration)
+
+    # Check various properties
+    np.testing.assert_equal(psf.centroid, galsim.PositionD(0,0))
+    np.testing.assert_almost_equal(psf.maxk, 3.306939635357677)
+    np.testing.assert_almost_equal(psf.stepk, 0.027742458082373515)
+    np.testing.assert_almost_equal(psf.kValue(0,0), test_flux+0j)
+    np.testing.assert_almost_equal(psf.xValue(0,0), 3.4955744341366577)
+    np.testing.assert_almost_equal(psf.kValue(0,0), (1+0j) * test_flux)
+    np.testing.assert_almost_equal(psf.flux, test_flux)
+    np.testing.assert_almost_equal(psf.xValue(0,0), psf.max_sb)
+
+    # Check that stepk and maxk scale correctly with lam/D
+    psf2 = galsim.Airy(lam_over_diam=5*test_loD, flux=test_flux, obscuration=test_obscuration)
+    np.testing.assert_almost_equal(psf2.maxk, psf.maxk/5)
+    np.testing.assert_almost_equal(psf2.stepk, psf.stepk/5)
+
+    # Check input flux vs output flux
+    for inFlux in np.logspace(-2, 2, 10):
+        psf3 = galsim.Airy(lam_over_diam=test_loD, flux=inFlux, obscuration=test_obscuration)
+        outFlux = psf3.flux
+        np.testing.assert_almost_equal(outFlux, inFlux)
+
+
 @timer
 def test_airy_shoot():
     """Test Airy with photon shooting.  Particularly the flux of the final image.

tests/test_box.py

@@ -241,6 +241,67 @@ def test_tophat():
                 approx_maxsb=True, scale=0.2)
     do_kvalue(conv,im, "Sheared TopHat convolved with pixel in real space")
 
+@timer
+def test_box_properties():
+    """Test some basic properties of the Box profile.
+    """
+    # Regression test based on v2.3.5 version of the code.
+
+    test_flux = 17.9
+    box = galsim.Box(width=0.2, height=0.25, flux=test_flux)
+
+    # Check various properties
+    np.testing.assert_equal(box.centroid, galsim.PositionD(0,0))
+    np.testing.assert_almost_equal(box.maxk, 10000)
+    np.testing.assert_almost_equal(box.stepk, 12.566370614359172)
+    np.testing.assert_almost_equal(box.kValue(0,0), test_flux+0j)
+    np.testing.assert_almost_equal(box.xValue(0,0), 358.0)
+    np.testing.assert_almost_equal(box.kValue(0,0), (1+0j) * test_flux)
+    np.testing.assert_almost_equal(box.flux, test_flux)
+    np.testing.assert_almost_equal(box.xValue(0,0), box.max_sb)
+
+    # Check that stepk and maxk scale correctly with lam/D
+    box2 = galsim.Box(width=5*0.2, height=5*0.25, flux=test_flux)
+    np.testing.assert_almost_equal(box2.maxk, box.maxk/5)
+    np.testing.assert_almost_equal(box2.stepk, box.stepk/5)
+
+    # Check input flux vs output flux
+    for inFlux in np.logspace(-2, 2, 10):
+        box3 = galsim.Box(width=0.2, height=0.25, flux=inFlux)
+        outFlux = box3.flux
+        np.testing.assert_almost_equal(outFlux, inFlux)
+
+
+@timer
+def test_tophat_properties():
+    """Test some basic properties of the TopHat profile.
+    """
+    # Regression test based on v2.3.5 version of the code.
+
+    test_flux = 17.9
+    tophat = galsim.TopHat(radius=0.23, flux=test_flux)
+
+    # Check various properties
+    np.testing.assert_equal(tophat.centroid, galsim.PositionD(0,0))
+    np.testing.assert_almost_equal(tophat.maxk, 593.7252723959091)
+    np.testing.assert_almost_equal(tophat.stepk, 13.659098493868665)
+    np.testing.assert_almost_equal(tophat.kValue(0,0), test_flux+0j)
+    np.testing.assert_almost_equal(tophat.xValue(0,0), 107.70788209243577)
+    np.testing.assert_almost_equal(tophat.kValue(0,0), (1+0j) * test_flux)
+    np.testing.assert_almost_equal(tophat.flux, test_flux)
+    np.testing.assert_almost_equal(tophat.xValue(0,0), tophat.max_sb)
+
+    # Check that stepk and maxk scale correctly with lam/D
+    tophat2 = galsim.TopHat(radius=5*0.23, flux=test_flux)
+    np.testing.assert_almost_equal(tophat2.maxk, tophat.maxk/5)
+    np.testing.assert_almost_equal(tophat2.stepk, tophat.stepk/5)
+
+    # Check input flux vs output flux
+    for inFlux in np.logspace(-2, 2, 10):
+        tophat3 = galsim.TopHat(radius=0.2, flux=inFlux)
+        outFlux = tophat3.flux
+        np.testing.assert_almost_equal(outFlux, inFlux)
+
 
 @timer
 def test_box_shoot():

tests/test_exponential.py

@@ -110,6 +110,11 @@ def test_exponential_properties():
         outFlux = expon.flux
         np.testing.assert_almost_equal(outFlux, inFlux)
 
+    # Check that stepk and maxk scale correctly with radius
+    expon2 = galsim.Exponential(flux=test_flux, scale_radius=5*test_scale)
+    np.testing.assert_almost_equal(expon2.maxk, expon.maxk/5)
+    np.testing.assert_almost_equal(expon2.stepk, expon.stepk/5)
+
 
 @timer
 def test_exponential_radii():

tests/test_gaussian.py

@@ -150,7 +150,7 @@ def test_gaussian_properties():
     np.testing.assert_almost_equal(gauss.xValue(cen), gauss.max_sb)
     # Check input flux vs output flux
     for inFlux in np.logspace(-2, 2, 10):
-        gauss = galsim.Gaussian(flux=inFlux, sigma=2.)
+        gauss = galsim.Gaussian(flux=inFlux, sigma=test_sigma)
         outFlux = gauss.flux
         np.testing.assert_almost_equal(outFlux, inFlux)
 
@@ -169,6 +169,10 @@ def test_gaussian_properties():
     assert_raises(TypeError, gauss.xValue, cen.x, cen.y, invalid=True)
     assert_raises(TypeError, gauss.xValue, pos=cen)
 
+    # Check that stepk and maxk scale correctly with radius
+    gauss2 = galsim.Gaussian(flux=test_flux, sigma=5*test_sigma)
+    np.testing.assert_almost_equal(gauss2.maxk, gauss.maxk/5)
+    np.testing.assert_almost_equal(gauss2.stepk, gauss.stepk/5)
 
 
 @timer

tests/test_inclined.py

@@ -527,6 +527,11 @@ def test_k_limits(run_slow):
             total_flux = np.sum(test_image.array)
             assert (total_flux-contained_flux)/total_flux <= gsparams.folding_threshold
 
+            # Check that stepk and maxk scale correctly with size
+            test_profile2 = get_prof(mode, inc_angle * galsim.radians,
+                                     5*scale_radius, 5*scale_height)
+            np.testing.assert_almost_equal(test_profile2.maxk, test_profile.maxk/5)
+            np.testing.assert_almost_equal(test_profile2.stepk, test_profile.stepk/5)
 
 @timer
 def test_eq_ne():

tests/test_kolmogorov.py

@@ -162,6 +162,11 @@ def test_kolmogorov_properties():
         np.testing.assert_almost_equal(out_flux, test_flux, 3,
                                        err_msg="Flux of Kolmogorov (image array) is incorrect.")
 
+    # Check that stepk and maxk scale correctly with radius
+    psf2 = galsim.Kolmogorov(lam_over_r0=5*lor, flux=test_flux)
+    np.testing.assert_almost_equal(psf2.maxk, psf.maxk/5)
+    np.testing.assert_almost_equal(psf2.stepk, psf.stepk/5)
+
 
 @timer
 def test_kolmogorov_radii():

tests/test_moffat.py

@@ -166,6 +166,19 @@ def test_moffat_properties():
         outFlux = psfFlux.flux
         np.testing.assert_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)
+    np.testing.assert_almost_equal(psf2.maxk, psf.maxk/5)
+    np.testing.assert_almost_equal(psf2.stepk, psf.stepk/5)
+
+    # Repeat without truncation
+    psf = galsim.Moffat(beta=5.0, half_light_radius=1., flux=test_flux)
+    psf2 = galsim.Moffat(beta=5.0, half_light_radius=5., flux=test_flux)
+    np.testing.assert_almost_equal(psf2.maxk, psf.maxk/5)
+    np.testing.assert_almost_equal(psf2.stepk, psf.stepk/5)
+
+
 @timer
 def test_moffat_maxk():
     """Check accuracy of maxk given maxk_threshold

tests/test_second_kick.py

@@ -156,6 +156,40 @@ def test_limiting_cases():
             rtol=1e-3,
             atol=1e-4)
 
+@timer
+def test_sk_properties():
+    """Test some basic properties of the SecondKick profile.
+    """
+    # Regression test based on v2.3.5 version of the code.
+
+    test_flux = 1.8
+    sk = galsim.SecondKick(lam=500, r0=0.2, diam=8, obscuration=0.6, kcrit=2, flux=test_flux)
+
+    # Check various properties
+    np.testing.assert_equal(sk.centroid, galsim.PositionD(0,0))
+    np.testing.assert_almost_equal(sk.maxk, 487.3878716587337)
+    np.testing.assert_almost_equal(sk.stepk, 1.5080215526716452)
+    np.testing.assert_almost_equal(sk.kValue(0,0), test_flux+0j)
+    np.testing.assert_almost_equal(sk.kValue(0,0), (1+0j) * test_flux)
+    np.testing.assert_almost_equal(sk.flux, test_flux)
+
+    # Check input flux vs output flux
+    for inFlux in np.logspace(-2, 2, 10):
+        skFlux = galsim.SecondKick(lam=500, r0=0.2, diam=8, obscuration=0.6, kcrit=2, flux=inFlux)
+        outFlux = skFlux.flux
+        np.testing.assert_almost_equal(outFlux, inFlux)
+
+    # Check that stepk and maxk scale correctly with r0,L0
+    sk2 = galsim.SecondKick(lam=500, r0=0.2/5, diam=8/5, obscuration=0.6, kcrit=2, flux=test_flux)
+    np.testing.assert_almost_equal(sk2.maxk, sk.maxk/5)
+    np.testing.assert_almost_equal(sk2.stepk, sk.stepk/5)
+
+    # Equivalent if scale lam instead.
+    sk2 = galsim.SecondKick(lam=5*500, r0=0.2, diam=8, obscuration=0.6, kcrit=2, flux=test_flux)
+    np.testing.assert_almost_equal(sk2.maxk, sk.maxk/5)
+    np.testing.assert_almost_equal(sk2.stepk, sk.stepk/5)
+
+
 @timer
 def test_sk_phase_psf(run_slow):
     """Test that analytic second kick profile matches what can be obtained from PhaseScreenPSF with

tests/test_sersic.py

@@ -486,6 +486,90 @@ def test_sersic_1():
         np.testing.assert_almost_equal(sersic.xValue(pos), expon.xValue(pos), decimal=5)
         np.testing.assert_almost_equal(sersic.kValue(pos), expon.kValue(pos), decimal=5)
 
+@timer
+def test_sersic_properties():
+    """Test some basic properties of the Sersic profile.
+    """
+    # Regression test based on v2.3.5 version of the code.
+
+    test_flux = 1.8
+    test_hlr = 2.3
+    test_scale_radius_trunc = 0.04559626861897289
+    test_scale_radius_notrunc = 0.009526646777227634
+    test_n = 3.1
+    prof = galsim.Sersic(n=test_n, half_light_radius=test_hlr, trunc=2*test_hlr, flux=test_flux)
+
+    # Check various properties
+    np.testing.assert_equal(prof.centroid, galsim.PositionD(0,0))
+    np.testing.assert_almost_equal(prof.maxk, 24.663782559047256)
+    np.testing.assert_almost_equal(prof.stepk, 0.6829549246934334)
+    np.testing.assert_almost_equal(prof.kValue(0,0), test_flux+0j)
+    np.testing.assert_almost_equal(prof.half_light_radius, test_hlr)
+    np.testing.assert_almost_equal(prof.xValue(0,0), 1.0238648164347117)
+    np.testing.assert_almost_equal(prof.kValue(0,0), (1+0j) * test_flux)
+    np.testing.assert_almost_equal(prof.flux, test_flux)
+    np.testing.assert_almost_equal(prof.xValue(0,0), prof.max_sb)
+    np.testing.assert_almost_equal(prof.scale_radius, test_scale_radius_trunc)
+    np.testing.assert_almost_equal(prof.half_light_radius, test_hlr)
+
+    # Now create the same profile using scale_radius
+    prof = galsim.Sersic(n=test_n, scale_radius=test_scale_radius_trunc, trunc=2*test_hlr,
+                         flux=test_flux)
+    np.testing.assert_almost_equal(prof.maxk, 24.663782559047256)
+    np.testing.assert_almost_equal(prof.stepk, 0.6829549246934334)
+    np.testing.assert_almost_equal(prof.kValue(0,0), test_flux+0j)
+    np.testing.assert_almost_equal(prof.half_light_radius, test_hlr)
+    np.testing.assert_almost_equal(prof.xValue(0,0), 1.0238648164347117)
+    np.testing.assert_almost_equal(prof.kValue(0,0), (1+0j) * test_flux)
+    np.testing.assert_almost_equal(prof.flux, test_flux)
+    np.testing.assert_almost_equal(prof.xValue(0,0), prof.max_sb)
+    np.testing.assert_almost_equal(prof.scale_radius, test_scale_radius_trunc)
+    np.testing.assert_almost_equal(prof.half_light_radius, test_hlr)
+
+    # Check that stepk and maxk scale correctly with radius
+    prof2 = galsim.Sersic(n=test_n, half_light_radius=5*test_hlr, trunc=5*2*test_hlr)
+    np.testing.assert_almost_equal(prof2.maxk, prof.maxk/5)
+    np.testing.assert_almost_equal(prof2.stepk, prof.stepk/5)
+
+    # Check input flux vs output flux
+    for inFlux in np.logspace(-2, 2, 10):
+        prof3 = galsim.Sersic(n=test_n, half_light_radius=test_hlr, trunc=2*test_hlr, flux=inFlux)
+        outFlux = prof3.flux
+        np.testing.assert_almost_equal(outFlux, inFlux)
+
+    # Repeat without truncation
+    prof = galsim.Sersic(n=test_n, half_light_radius=test_hlr, flux=test_flux)
+
+    # Check various properties
+    np.testing.assert_almost_equal(prof.maxk, 50.35677579928042)
+    np.testing.assert_almost_equal(prof.stepk, 0.08358929429681608)
+    np.testing.assert_almost_equal(prof.kValue(0,0), test_flux+0j)
+    np.testing.assert_almost_equal(prof.half_light_radius, test_hlr)
+    np.testing.assert_almost_equal(prof.xValue(0,0), 6.010654623502727)
+    np.testing.assert_almost_equal(prof.kValue(0,0), (1+0j) * test_flux)
+    np.testing.assert_almost_equal(prof.flux, test_flux)
+    np.testing.assert_almost_equal(prof.xValue(0,0), prof.max_sb)
+    np.testing.assert_almost_equal(prof.scale_radius, test_scale_radius_notrunc)
+    np.testing.assert_almost_equal(prof.half_light_radius, test_hlr)
+
+    # Now create the same profile using scale_radius
+    prof = galsim.Sersic(n=test_n, scale_radius=test_scale_radius_notrunc, flux=test_flux)
+    np.testing.assert_almost_equal(prof.maxk, 50.35677579928042)
+    np.testing.assert_almost_equal(prof.stepk, 0.08358929429681608)
+    np.testing.assert_almost_equal(prof.kValue(0,0), test_flux+0j)
+    np.testing.assert_almost_equal(prof.half_light_radius, test_hlr)
+    np.testing.assert_almost_equal(prof.xValue(0,0), 6.010654623502727)
+    np.testing.assert_almost_equal(prof.kValue(0,0), (1+0j) * test_flux)
+    np.testing.assert_almost_equal(prof.flux, test_flux)
+    np.testing.assert_almost_equal(prof.xValue(0,0), prof.max_sb)
+    np.testing.assert_almost_equal(prof.scale_radius, test_scale_radius_notrunc)
+    np.testing.assert_almost_equal(prof.half_light_radius, test_hlr)
+
+    # Check that stepk and maxk scale correctly with radius
+    prof2 = galsim.Sersic(n=test_n, half_light_radius=5*test_hlr)
+    np.testing.assert_almost_equal(prof2.maxk, prof.maxk/5)
+    np.testing.assert_almost_equal(prof2.stepk, prof.stepk/5)
+
 
 @timer
 def test_sersic_shoot():

tests/test_shapelet.py

@@ -180,6 +180,10 @@ def test_shapelet_properties():
     assert_raises(TypeError, galsim.Shapelet, order=order, bvec=bvec)
     assert_raises(ValueError, galsim.Shapelet, sigma=sigma, order=5, bvec=bvec)
 
+    # Check that stepk and maxk scale correctly with radius
+    shapelet2 = galsim.Shapelet(sigma=5*sigma, order=order, bvec=bvec)
+    np.testing.assert_almost_equal(shapelet2.maxk, shapelet.maxk/5)
+    np.testing.assert_almost_equal(shapelet2.stepk, shapelet.stepk/5)
 
 @timer
 def test_shapelet_fit():