Merge pull request #197 from robelgeda/tests

General Tests [v0.5.0]

petrofit/modeling/models.py

@@ -261,7 +261,7 @@ def oversample(self, oversample):
             else:
                 raise TypeError("oversample should be a single int factor or a tuple (x, y, size, factor).")
             if grid_factor <= 0:
-                raise TypeError("oversample should be a positive int factor.")
+                raise ValueError("oversample should be a positive int factor.")
             if grid_factor % psf_oversample != 0:
                 raise ValueError('oversample should be equal to or an integer multiple of psf_oversample')
         self._oversample = oversample
@@ -280,7 +280,7 @@ def psf_oversample(self, psf_oversample):
             if not isinstance(psf_oversample, int) or psf_oversample <= 0:
                 raise TypeError("psf_oversample should be a single positive int factor.")
             if grid_factor % psf_oversample != 0:
-                raise ValueError('oversample should be equal to or an integer multiple of psf_oversample.'
+                raise ValueError('oversample should be equal to or an integer multiple of psf_oversample. ' 
                                  'Set PSFConvolvedModel2D.oversample value first.')
         self._psf_oversample = psf_oversample
 
@@ -406,14 +406,14 @@ def evaluate(self, x, y, *params, **kwargs):
                 # Sample the sub-model onto the sub-grid
                 sub_model_oversampled_image = self._model.evaluate(x_sub_grid, y_sub_grid, *sub_model_params)
 
-                # Block reduce the window to the main image resolution
+                # Block reduce the window to the psf resolution
                 sub_model_image = block_reduce(sub_model_oversampled_image, sub_grid_to_psf_factor) / sub_grid_to_psf_factor ** 2
 
-                # Compute window indices in main image frame
+                # Compute window indices in main image frame at data resolution first
                 i_sub_min = int(np.round(sub_grid_origin[0]))
                 j_sub_min = int(np.round(sub_grid_origin[1]))
-                i_sub_max = i_sub_min + sub_grid_size
-                j_sub_max = j_sub_min + sub_grid_size
+                i_sub_max = (i_sub_min + sub_grid_size)
+                j_sub_max = (j_sub_min + sub_grid_size)
 
                 # Clip window indices
                 if i_sub_min < 0:
@@ -425,6 +425,12 @@ def evaluate(self, x, y, *params, **kwargs):
                 if j_sub_max > j.max():
                     j_sub_max = j.max() + 1
 
+                # Convert to PSF resolution
+                i_sub_min *= psf_factor
+                j_sub_min *= psf_factor
+                i_sub_max *= psf_factor
+                j_sub_max *= psf_factor
+
                 # Add oversampled window to image
                 model_image[
                     j_sub_min:j_sub_max,

petrofit/tests/conftest.py

@@ -23,7 +23,7 @@ def output_dir():
 @pytest.fixture
 def galfit_images(data_dir):
     """fixture for galfit test images used for testing"""
-    path = os.path.join(data_dir, 'galfit_*fits.gz')
+    path = os.path.join(data_dir, 'galfit_sersic*fits.gz')
     fb = glob(path)
 
     galfit_image_dict = {}
@@ -35,17 +35,25 @@ def galfit_images(data_dir):
     return galfit_image_dict
 
 
+@pytest.fixture
+def galfit_psf_images(data_dir):
+    """fixture for galfit test images used for testing"""
+    path = os.path.join(data_dir, 'galfit_n4_psf*.gz')
+    fb = glob(path)
 
+    galfit_psf_image_dict = {}
+    for f in fb:
+        name = int(os.path.basename(f).replace('galfit_n4_psf', '').replace('.fits.gz', ''))
+        image = fits.getdata(f)
+        galfit_psf_image_dict[name] = image
 
+    return galfit_psf_image_dict
 
-# @pytest.fixture
-# def sersic_2d_image():
-#     """fixture for Sersic 2D image """
-#     path = "sersic_2d_image.fits.gz"
-#     sersic_2d_path =  os.path.join(os.path.dirname(__file__), path)
-#
-#     if not os.path.isfile(sersic_2d_path):
-#         make_image()
-#
-#     return CCDData.read(sersic_2d_path)
 
+@pytest.fixture
+def psf_image(data_dir):
+    """fixture for PSF test images used for testing"""
+    path = os.path.join(data_dir, 'f105w_psf.fits.gz')
+    PSF = fits.getdata(path)
+    PSF /= PSF.sum()
+    return PSF

petrofit/tests/test_fitting.py

@@ -1,10 +1,13 @@
 import numpy as np
+import pytest
 
 from astropy.convolution import convolve
 from astropy.modeling import models
 
 import petrofit as pf
 
+from matplotlib import pyplot as plt
+
 
 def test_psf_convolved_image_model():
     """Test fitting and PSF convolution"""
@@ -73,3 +76,159 @@ def test_psf_convolved_image_model():
     # Check if fit is close to actual
     error_arr = abs(fitted_model_image - psf_sersic_image) / psf_sersic_image
     assert np.max(error_arr) < 0.01  # max error less than 1% error
+
+
+def test_make_grid():
+    # Test default arguments
+    x, y = pf.make_grid(3)
+    expected_x = np.array([[0., 1., 2.],
+                           [0., 1., 2.],
+                           [0., 1., 2.]])
+    expected_y = np.array([[0., 0., 0.],
+                           [1., 1., 1.],
+                           [2., 2., 2.]])
+    assert np.allclose(x, expected_x)
+    assert np.allclose(y, expected_y)
+
+    # Test with origin and factor
+    x, y = pf.make_grid(3, origin=(4, 5))
+    assert np.allclose(x, expected_x + 4)
+    assert np.allclose(y, expected_y + 5)
+
+    x, y = pf.make_grid(2, origin=(0.25, 0.25), factor=2)
+    expected_x, expected_y = [
+        np.array([
+            [0., 0.5, 1., 1.5],
+            [0., 0.5, 1., 1.5],
+            [0., 0.5, 1., 1.5],
+            [0., 0.5, 1., 1.5]
+        ]),
+         np.array([
+             [0., 0., 0., 0.],
+             [0.5, 0.5, 0.5, 0.5],
+             [1., 1., 1., 1.],
+             [1.5, 1.5, 1.5, 1.5]
+         ])]
+    assert np.allclose(x, expected_x)
+    assert np.allclose(y, expected_y)
+
+    # Test factor assertion
+    with pytest.raises(AssertionError):
+        pf.make_grid(3, factor=1.5)
+
+    # Test origin size
+    with pytest.raises(IndexError):
+        pf.make_grid(3, origin=(1,))
+
+
+def test_psf_convolved_model_2d_init():
+    with pytest.raises(AssertionError):
+        pf.PSFConvolvedModel2D(models.Sersic1D())
+
+    with pytest.raises(AssertionError):
+        pf.PSFConvolvedModel2D(models.Sersic1D())
+
+    I_e = 1
+    n = 1
+    base_model = models.Sersic2D(amplitude=I_e, r_eff=5, n=n, x_0=15., y_0=15.)
+    psf = pf.model_to_image(models.Gaussian2D(x_mean=15, y_mean=15, x_stddev=5, y_stddev=5), size=30)
+    psf /= psf.sum()
+
+    pf.PSFConvolvedModel2D(base_model)
+    pf.PSFConvolvedModel2D(base_model, psf=None, oversample=None, psf_oversample=None)
+    pf.PSFConvolvedModel2D(base_model, psf=None, oversample=5, psf_oversample=None)
+
+    pf.PSFConvolvedModel2D(base_model, psf=psf, oversample=None, psf_oversample=None)
+    pf.PSFConvolvedModel2D(base_model, psf=psf, oversample=5, psf_oversample=None)
+    pf.PSFConvolvedModel2D(base_model, psf=psf, oversample=10, psf_oversample=5)
+    model = pf.PSFConvolvedModel2D(base_model, psf=psf, oversample=5, psf_oversample=5)
+
+    assert model.oversample == 5
+    assert model.psf_oversample == 5
+
+    with pytest.raises(ValueError):
+        # psf_oversample provided but PSF is None
+        pf.PSFConvolvedModel2D(base_model, psf=None, oversample=5, psf_oversample=5)
+
+    with pytest.raises(ValueError):
+        # oversample should be equal to or an integer multiple of psf_oversample
+        pf.PSFConvolvedModel2D(base_model, psf=psf, oversample=1, psf_oversample=5)
+        pf.PSFConvolvedModel2D(base_model, psf=psf, oversample=9, psf_oversample=5)
+
+    psf_model = pf.PSFConvolvedModel2D(base_model)
+    psf_model.psf = psf
+    psf_model.oversample = 4
+    psf_model.psf_oversample = 4
+
+    psf_model = pf.PSFConvolvedModel2D(base_model)
+    psf_model.oversample = 4
+    with pytest.raises(ValueError):
+        psf_model.psf_oversample = 4
+
+    psf_model = pf.PSFConvolvedModel2D(base_model)
+    psf_model.psf = psf
+    psf_model.oversample = None
+    with pytest.raises(ValueError):
+        # oversample should be equal to or an integer multiple of psf_oversample
+        psf_model.psf_oversample = 4
+
+    # Additional Tests
+    with pytest.raises(TypeError):
+        # Test for float oversample value
+        pf.PSFConvolvedModel2D(base_model, psf=psf, oversample=5.5, psf_oversample=5)
+
+    with pytest.raises(ValueError):
+        # Test for negative oversample value
+        pf.PSFConvolvedModel2D(base_model, psf=psf, oversample=-5, psf_oversample=5)
+
+    psf_model = pf.PSFConvolvedModel2D(base_model)
+    psf_model.psf = psf
+    psf_model.oversample = 4
+    psf_model.psf_oversample = 4
+
+    psf_model = pf.PSFConvolvedModel2D(base_model)
+    psf_model.psf = psf
+    psf_model.oversample = 8
+    psf_model.psf_oversample = 4
+    assert psf_model._get_psf_factor() == 4
+    assert psf_model._get_oversample_factor() == 8
+    psf_model.oversample = 16
+    assert psf_model._get_oversample_factor() == 16
+    psf_model.oversample = (15, 15, 100, 24)
+    assert psf_model._get_oversample_factor() == 24
+    psf_model.psf_oversample = None
+    assert psf_model._get_psf_factor() == 1
+    psf_model.oversample = 3
+    assert psf_model._get_oversample_factor() == 3
+    psf_model.oversample = None
+    assert psf_model._get_psf_factor() == psf_model._get_oversample_factor() == 1
+
+
+def test_psf_sampling(galfit_psf_images, psf_image):
+    for psf_index, image in galfit_psf_images.items():
+        base_model = models.Sersic2D(
+            amplitude=0.020625826413226116,
+            r_eff=30,
+            n=4,
+            x_0=99.,
+            y_0=99.,
+            ellip=0,
+            theta=0.0,
+            bounds=pf.get_default_sersic_bounds(),
+            fixed={'x_0': False, 'y_0': False, 'n': True, 'r_eff': True, 'ellip': True, 'theta': True}
+        )
+
+        PSF = psf_image if psf_index > 0 else None
+        psf_oversample = psf_index if psf_index > 0 else None
+        model = pf.PSFConvolvedModel2D(base_model, psf=PSF,
+                                       oversample=('x_0', 'y_0', 100, psf_oversample * 5 if PSF is not None else 50),
+                                       psf_oversample=psf_oversample if PSF is not None else None
+                                       )
+        model.fixed.update({'psf_pa': True})
+
+        fitted_model, fit = pf.fit_model(image, model, acc=1e-12, )
+
+        *_, residual = pf.plot_fit(fitted_model, image, vmax=0.039)
+        assert abs(residual).max() < 0.05
+        pf.print_model_params(fitted_model)
+        plt.show()