Added "ftol" as optional input parameter for Imfit.doFit() and Imfit.…

…fit().

Added basic unit testing of Imfit.doFit (checking if we get correct chi^2, AIC,
parametr values, errors).

CHANGELOG.md

@@ -12,6 +12,10 @@ Models (including current parameter values) can now be described by a dict-based
 dict-based model descriptions (including current best-fit parameter values) can be
 returned by Imfit instances as well.
 
+Imfit.fit and Imfit.doFit can now take an optional `ftol` parameter (same as the `--ftol`
+parameter for the command-line `imfit` program -- controls fractional tolerance of fit statistics
+as a convergence criterion).
+
 Pre-compiled version for Python versions 3.9 and 3.10 on macOS.
 
 ### Changed

docs/sample_usage.md

@@ -101,11 +101,11 @@ and passing the parent dict to the `ModelObject.dict_to_ModelDescription` functi
         # dict describing the bulge (first define the parameter dict, with initial values
         # and lower & upper limits for each parameter)
         p_bulge = {'PA': [PA_bulge, 0, 180], 'ell_bulge': [ell, 0, 1], 'n': [n, 0.5, 5], 
-                    'I_e': I_e, 0.0, 10*I_e]], 'r_e': [r_e, 0.0, 10*r_e}
+                    'I_e': [I_e, 0.0, 10*I_e], 'r_e': [r_e, 0.0, 10*r_e]}
         bulge_dict = {'name': "Sersic", 'label': "bulge", 'parameters': p_bulge}
         # do the same thing for the disk component
         p_disk = {'PA': [PA_disk, 0, 180], 'ell_disk': [ell, 0, 1], 'I_0': [I_0, 0, 10*I_0],
-                    'h': [h, 0.0, 10*h}
+                    'h': [h, 0.0, 10*h]}
         disk_dict = {'name': "Exponential", 'label': "disk", 'parameters': p_disk}
 
         # make dict for the function set that combines the bulge and disk components

new_notes_and_todo_for_python-wrapper_code.txt

@@ -2,6 +2,14 @@ PLAN:
 
 ** GENERAL/FUTURE TO-DO
 
+[X] Add ftol keyword to fitting.py --
+    Imfit.fit and Imfit.doFit
+
+
+[ ] Modify various methods in Imfit class to check whether _modelObjectWrapper has been
+    instantiated or not.
+
+
     ** Update GitHub access to new methodology
         [X] Visit GitHub web page
         [X]1. Generate Personal Authenticaion Token [PAT] (& copy to pword file)
@@ -64,21 +72,25 @@ PLAN:
                     -- 4. Do new "develop" install and run tests on Linux VM
 
 
-    ** New Update: Make build for Python 3.9
+    ** New Update: Make build for Python 3.10
         howto_new_distribution.txt
 
-        [X] Install Python 3.9
-            [X] Install Python 3.9
+        [X] Install Python 3.10
+            [X] Install Python
             [X] Install necessary libraries via pip
 
-        [X] Add Python 3.9 as an option for GitHub Actions testing workflow
+        [X] Add Python 3.10 as an option for GitHub Actions testing workflow
 
         [X] Upload new version to pip staging area
 
         [ ] Email announcement to Imfit mailing list
 
         [ ] Email announcement to AstroPy mailing list
 
+        [ ] Announce on Facebook Python-Astronomy group
+
+        [ ] Announce on Facebook Astrostatistics group
+
 
     ** Alternate outputs for best-fit parameter values
         https://github.com/perwin/pyimfit/issues/2
@@ -190,6 +202,13 @@ PLAN:
 
 
     * Figure out fix for handling makeimage-style image generation
+        Three levels/stages:
+            0. Instantiate Imfit instance
+                -- optionally supply PSF, subsampling flag, zero point
+            1. Instantiate ModelObjectWrapper *without* image size
+                -- can then be used for e.g. getModelFluxes()
+            2. Specify image size (e.g. via loadData)
+
         -- Currently, we ignore NCOLS/NROWS image-description parameters in a config file
         (config.py) [that is, we load and store the values, but otherwise ignore them]
         -- Remember that model-image size combines with PSF size to set the internal model-image array,

pyimfit/descriptions.py

@@ -131,7 +131,8 @@ def setValue( self, value: float, limits: Optional[Sequence[float]]=None, fixed:
             elif value > upper_limit:
                 upper_limit = value
             if lower_limit >= upper_limit:
-                raise ValueError("Lower limit must be < upper limit.")
+                err_msg = "Lower limit ({0}) must be < upper limit ({1}).".format(lower_limit, upper_limit)
+                raise ValueError(err_msg)
             self._limits = (float(lower_limit), float(upper_limit))
 
         self._value = float(value)

pyimfit/fitting.py

@@ -530,7 +530,7 @@ def loadData(self, image, error=None, mask=None, **kwargs ):
         self._dataSet = True
 
 
-    def doFit( self, solver='LM', verbose=None ):
+    def doFit( self, solver='LM', ftol=1e-8, verbose=None ):
         """
         Fit the model to previously supplied data image.
 
@@ -542,6 +542,9 @@ def doFit( self, solver='LM', verbose=None ):
                 * ``'NM'`` : Nelder-Mead Simplex.
                 * ``'DE'`` : Differential Evolution.
 
+        ftol : float, optional
+            fractional tolerance in fit statistic for determining fit convergence
+
         verbose : int or None, optional
             set this to an integer to specify a feedback level for the fit (this overrides
             the Imfit object's internal verbosity setting)
@@ -567,15 +570,15 @@ def doFit( self, solver='LM', verbose=None ):
             verboseLevel = verbose
         else:
             verboseLevel = self._verboseLevel
-        self._modelObjectWrapper.fit(verbose=verboseLevel, mode=solver)
+        self._modelObjectWrapper.fit(verbose=verboseLevel, mode=solver, ftol=ftol)
         self._lastSolverUsed = solver
         if not self.fitError:
             self._fitDone = True
             self._fitStatComputed = True
         return self.getFitResult()
 
 
-    def fit( self, image, error=None, mask=None, solver='LM', verbose=None, **kwargs ):
+    def fit( self, image, error=None, mask=None, solver='LM', ftol=1e-8, verbose=None, **kwargs ):
         """
         Supply data image (and optionally mask and/or error images) and image info, then
         fit the model to the data.
@@ -603,6 +606,9 @@ def fit( self, image, error=None, mask=None, solver='LM', verbose=None, **kwargs
                 * ``'NM'`` : Nelder-Mead Simplex.
                 * ``'DE'`` : Differential Evolution.
 
+        ftol : float, optional
+            fractional tolerance in fit statistic for determining fit convergence
+
         verbose : int or None, optional
             set this to an integer to specify a feedback level for the fit (this overrides
             the Imfit object's internal verbosity setting)
@@ -614,7 +620,7 @@ def fit( self, image, error=None, mask=None, solver='LM', verbose=None, **kwargs
         result : FitResult object
         """
         self.loadData(image, error, mask, **kwargs)
-        return self.doFit(solver=solver, verbose=verbose)
+        return self.doFit(solver=solver, ftol=ftol, verbose=verbose)
 
 
     def getFitResult( self ):

pyimfit/pyimfit_lib.pyx

@@ -1,6 +1,6 @@
 # Cython implementation file for wrapping Imfit code, by PE; based on code
 # by Andre Luis de Amorim.
-# Copyright André Luiz de Amorim, 2013; Peter Erwin, 2018-2019.
+# Copyright André Luiz de Amorim, 2013; Peter Erwin, 2018-2022.
 
 # Note that we are using "typed memoryviews" to translate numpy arrays into
 # C-style double * arrays; this apparently the preferred (newer, more flexible)

pyimfit/tests/test_fitting.py

@@ -23,6 +23,9 @@
 testDataDir = "../data/"
 imageFile = testDataDir + "ic3478rss_256.fits"
 configFile = testDataDir + "config_exponential_ic3478_256.dat"
+true_bestfit_params_ic3478 = np.array([128.8540,129.1028, 19.7266,0.23152,316.313,20.522])
+true_errs_ic3478 = np.array([0.0239,0.0293, 0.21721,0.0015524,0.61962,0.034674])
+
 imageFile2 = testDataDir + "n3073rss_small.fits"
 maskFile2 = testDataDir + "n3073rss_small_mask.fits"
 configFile2 = testDataDir + "config_n3073.dat"
@@ -109,6 +112,33 @@ def test_Imfit_optionsDict_updates( self ):
         imfit_fitter2._updateModelDescription(keywords_new)
         assert imfit_fitter2._modelDescr.optionsDict == optionsDict_correct3
 
+    def test_Imfit_simple_fit( self ):
+        # Fitting Exponential to 256x256-pixel SDSS r-band image of IC 3478 (no PSF convolution)
+        imfit_fitter = Imfit(self.modelDesc)
+        imfit_fitter.loadData(image_ic3478, gain=4.725, read_noise=4.3, original_sky=130.14)
+        # fit with defautl LM solver
+        imfit_fitter.doFit()
+        assert imfit_fitter.fitTerminated == False
+        assert imfit_fitter.fitConverged == True
+        assert imfit_fitter.fitStatistic == approx(136470.399329, rel=1e-10)
+        assert imfit_fitter.AIC == approx(136482.400611, rel=1e-10)
+        bestfit_params = imfit_fitter.getRawParameters()
+        assert_allclose(bestfit_params, true_bestfit_params_ic3478, rtol=1e-5)
+        bestfit_errs = imfit_fitter.getParameterErrors()
+        assert_allclose(bestfit_errs, true_errs_ic3478, rtol=1e-2)
+
+    def test_Imfit_simple_fit_altftol( self ):
+        # Fitting Exponential to 256x256-pixel SDSS r-band image of IC 3478 (no PSF convolution)
+        # this time we explicitly use ftol=1e-6 (instead of default 1e-8)
+        imfit_fitter = Imfit(self.modelDesc)
+        imfit_fitter.loadData(image_ic3478, gain=4.725, read_noise=4.3, original_sky=130.14)
+        # fit with defautl LM solver
+        imfit_fitter.doFit(ftol=1e-6)
+        assert imfit_fitter.fitTerminated == False
+        assert imfit_fitter.fitConverged == True
+        assert imfit_fitter.fitStatistic == approx(136470.402986, rel=1e-10)
+        assert imfit_fitter.AIC == approx(136482.404268, rel=1e-10)
+
     def test_Imfit_get_fluxes( self ):
         # Fitting Exponential to 256x256-pixel SDSS r-band image of IC 3478 (no PSF convolution)
         imfit_fitter = Imfit(self.modelDesc)