Merge pull request #41 from simontorres/fix_focus_when_data_is_not_good

Fix focus when data is not good

.gitignore

@@ -1,5 +1,6 @@
 docs/_build/
 .idea/*
 dist/*
+build/*
 *__pycache__*
 *.egg-info/*

CHANGES.rst

@@ -1,3 +1,12 @@
+.. _v2.0.2:
+
+2.0.2 08-04-2022
+================
+
+- Added exception handling for when it's not possible to use scipy.optimize.brent
+- Added notes field to report special warnings.
+
+
 .. _v2.0.1:
 
 2.0.1 23-02-2022

docs/conf.py

@@ -170,4 +170,4 @@
 # -- Options for intersphinx extension ---------------------------------------
 
 # Example configuration for intersphinx: refer to the Python standard library.
-intersphinx_mapping = {'https://docs.python.org/': None}
+intersphinx_mapping = {'https://docs.python.org/': None}

goodman_focus/__init__.py

@@ -1,3 +1,3 @@
-from .version import __version__
-from .goodman_focus import GoodmanFocus
-from .goodman_focus import run_goodman_focus
+from .version import __version__  # noqa: F401
+from .goodman_focus import GoodmanFocus  # noqa: F401
+from .goodman_focus import run_goodman_focus  # noqa: F401

goodman_focus/goodman_focus.py

@@ -10,7 +10,7 @@
 
 from astropy.io import fits
 from astropy.stats import sigma_clip
-from astropy.modeling import models, fitting, Model
+from astropy.modeling import models, fitting
 from ccdproc import CCDData
 from ccdproc import ImageFileCollection
 from scipy import optimize
@@ -104,8 +104,6 @@ def get_peaks(ccd, file_name='', plots=False):
 
     clipped_profile = sigma_clip(raw_profile, sigma=1, maxiters=5)
 
-    _mean = np.mean(clipped_profile)
-
     clipped_x_axis = [i for i in range(len(clipped_profile)) if not clipped_profile.mask[i]]
     cleaned_profile = clipped_profile[~clipped_profile.mask]
 
@@ -128,7 +126,6 @@ def get_peaks(ccd, file_name='', plots=False):
 
     peaks = signal.argrelmax(filtered_data, axis=0, order=5)[0]
 
-
     if len(peaks) == 1:
         log.debug(f"Found {len(peaks)} peak in file")
     else:
@@ -283,12 +280,16 @@ def __init__(self,
 
         self.__ccd = None
         self.file_name = None
+        self.__focus = None
+        self.__fwhm = None
+        self.__files = None
         self.__best_focus = None
         self.__best_fwhm = None
         self.__best_image = None
         self.__best_image_focus = None
         self.__best_image_fwhm = None
         self._fwhm = None
+        self.__notes = ''
 
         self.polynomial = models.Polynomial1D(degree=5)
         self.fitter = fitting.LevMarLSQFitter()
@@ -378,47 +379,50 @@ def __call__(self, files=None):
                 self.log.critical('"files" argument must be a list')
                 sys.exit(0)
 
-
         results = []
         for focus_group in self.focus_groups:
             mode_name = self._get_mode_name(focus_group)
-
-            focus_dataframe = self.get_focus_data(group=focus_group)
-
-            self._fit(df=focus_dataframe)
-            self.log.info(f"Best Focus for mode {mode_name} is {self.__best_focus}")
-            results.append({'date': focus_group['DATE'].tolist()[0],
-                            'time': focus_group['DATE-OBS'].tolist()[0],
-                            'mode_name': mode_name,
-                            'focus': round(self.__best_focus, 10),
-                            'fwhm': round(self.__best_fwhm, 10),
-                            'best_image_name': self.__best_image,
-                            'best_image_focus': round(self.__best_image_focus, 10),
-                            'best_image_fwhm': round(self.__best_image_fwhm, 10),
-                            'focus_data': focus_dataframe['focus'].tolist(),
-                            'fwhm_data':  focus_dataframe['fwhm'].tolist()
-                            })
-
-            if self.plot_results:   # pragma: no cover
-
-                fig, ax = plt.subplots()
-
-                focus_list = focus_dataframe['focus'].tolist()
-                fwhm_list = focus_dataframe['fwhm'].tolist()
-                new_x_axis = np.linspace(focus_list[0], focus_list[-1], 1000)
-
-                ax.plot(focus_list, fwhm_list, marker='x', label='Measured FWHM')
-                ax.axvline(self.__best_focus, color='k', label='Best Focus')
-                ax.set_title(f"Best Focus:\n{mode_name} {self.__best_focus:.3f}")
-                ax.set_xlabel("Focus Value")
-                if 'IM_' in mode_name:
-                    ax.set_ylabel("FWHM")
-                else:
-                    ax.set_ylabel("Mean FWHM")
-                ax.plot(new_x_axis,
-                        self.polynomial(new_x_axis), label='Model')
-                ax.legend(loc='best')
-                plt.show()
+            self.notes = ''
+            try:
+                focus_dataframe = self.get_focus_data(group=focus_group)
+
+                self._fit(df=focus_dataframe)
+                self.log.info(f"Best Focus for mode {mode_name} is {self.__best_focus}")
+                results.append({'date': focus_group['DATE'].tolist()[0],
+                                'time': focus_group['DATE-OBS'].tolist()[0],
+                                'mode_name': mode_name,
+                                'notes': self.notes,
+                                'focus': round(self.__best_focus, 10),
+                                'fwhm': round(self.__best_fwhm, 10),
+                                'best_image_name': self.__best_image,
+                                'best_image_focus': round(self.__best_image_focus, 10),
+                                'best_image_fwhm': round(self.__best_image_fwhm, 10),
+                                'focus_data': focus_dataframe['focus'].tolist(),
+                                'fwhm_data': focus_dataframe['fwhm'].tolist()
+                                })
+
+                if self.plot_results:   # pragma: no cover
+
+                    fig, ax = plt.subplots()
+
+                    focus_list = focus_dataframe['focus'].tolist()
+                    fwhm_list = focus_dataframe['fwhm'].tolist()
+                    new_x_axis = np.linspace(focus_list[0], focus_list[-1], 1000)
+
+                    ax.plot(focus_list, fwhm_list, marker='x', label='Measured FWHM')
+                    ax.axvline(self.__best_focus, color='k', label='Best Focus')
+                    ax.set_title(f"Best Focus:\n{mode_name} {self.__best_focus:.3f}")
+                    ax.set_xlabel("Focus Value")
+                    if 'IM_' in mode_name:
+                        ax.set_ylabel("FWHM")
+                    else:
+                        ax.set_ylabel("Mean FWHM")
+                    ax.plot(new_x_axis,
+                            self.polynomial(new_x_axis), label='Model')
+                    ax.legend(loc='best')
+                    plt.show()
+            except ValueError as error:
+                self.log.error(f"Unable to obtain focus due to ValueError: {str(error)}", exc_info=True)
 
         return results
 
@@ -440,23 +444,32 @@ def _fit(self, df):
         Returns:
 
         """
-        focus = df['focus'].tolist()
-        fwhm = df['fwhm'].tolist()
-        files = df['file'].tolist()
-        max_focus = np.max(focus)
-        min_focus = np.min(focus)
-        self.polynomial = self.fitter(self.polynomial, focus, fwhm)
-        self._get_local_minimum(x1=min_focus, x2=max_focus)
-
-        index = np.argmin(np.abs(focus - self.__best_focus))
-
-        self.__best_image = files[index]
-        self.__best_image_focus = focus[index]
-        self.__best_image_fwhm = fwhm[index]
+        self.__focus = df['focus'].tolist()
+        self.__fwhm = df['fwhm'].tolist()
+        self.__files = df['file'].tolist()
+        max_focus = np.max(self.__focus)
+        min_focus = np.min(self.__focus)
+        self.polynomial = self.fitter(self.polynomial, self.__focus, self.__fwhm)
+        try:
+            self._get_local_minimum(x1=min_focus, x2=max_focus)
+            self.notes = f"Focus obtained by using Brent's optimization method."
+        except ValueError as error:
+            self.log.error(f"Error finding local minimum with fitted data: {str(error)}")
+            self.log.warning(f"This method does not guarantee this is the best focus for this setup.")
+            self.notes = f"Warning: This result does not guarantee this is the best focus."
+            lowest_fwhm_index = np.argmin(self.__fwhm)
+            self.__best_focus = self.__focus[lowest_fwhm_index]
+            self.__best_fwhm = self.__fwhm[lowest_fwhm_index]
+
+        index = np.argmin(np.abs(np.array(self.__focus) - self.__best_focus))
+
+        self.__best_image = self.__files[index]
+        self.__best_image_focus = self.__focus[index]
+        self.__best_image_fwhm = self.__fwhm[index]
 
         return self.polynomial
 
-    def _get_local_minimum(self, x1, x2):
+    def _get_local_minimum(self, x1, x2, x_axis_size=2000):
         """Finds best focus
 
         Using scipy.optimize.brent method that uses an inverse parabolic algorithm.
@@ -471,7 +484,7 @@ def _get_local_minimum(self, x1, x2):
             values (i.e. closest to zero).
 
         """
-        x_axis = np.linspace(x1, x2, 2000)
+        x_axis = np.linspace(x1, x2, x_axis_size)
         modeled_data = self.polynomial(x_axis)
         index_of_minimum = np.argmin(modeled_data)
         middle_point = x_axis[index_of_minimum]
@@ -540,7 +553,7 @@ def get_focus_data(self, group):
                                                    self.file_name),
                                       unit='adu')
 
-            peaks, values, x_axis,  profile = get_peaks(
+            peaks, values, x_axis, profile = get_peaks(
                 ccd=self.__ccd,
                 file_name=self.file_name,
                 plots=self.debug)
@@ -591,14 +604,9 @@ def run_goodman_focus(args=None):   # pragma: no cover
                                  debug=args.debug)
 
     results = goodman_focus()
-    print(results)
     log.info("Summary")
     for result in results:
-        print(json.dumps(result, indent=4))
-        for key in result.keys():
-            log.info(f"Mode: {key} Best Focus: {result[key]['focus']} at FWHM: {result[key]['fwhm']}. "
-                     f"Best image: {result[key]['best_image']['file_name']} with focus: "
-                     f"{result[key]['best_image']['focus']} and FWHM: {result[key]['best_image']['fwhm']}")
+        log.info(json.dumps(result, indent=4))
 
 
 if __name__ == '__main__':   # pragma: no cover

goodman_focus/tests/test_goodman_focus.py

@@ -6,13 +6,12 @@
 
 from astropy.io import fits
 from astropy.modeling import models
-from unittest import TestCase, skip
+from unittest import TestCase
 from ccdproc import CCDData
 
 from ..goodman_focus import GoodmanFocus
 from ..goodman_focus import get_args, get_peaks, get_fwhm
 
-import matplotlib.pyplot as plt
 
 logging.disable(logging.CRITICAL)
 
@@ -109,13 +108,9 @@ def test_multiple_peaks(self):
 class GoodmanFocusTests(TestCase):
 
     def setUp(self):
-        arguments = ['--data-path', os.getcwd(),
-                     '--file-pattern', '*.fits',
-                     '--obstype', 'FOCUS',
-                     '--features-model', 'gaussian']
 
         number_of_test_subjects = 21
-        self.file_list = ["file_{}.fits".format(i+1) for i in range(number_of_test_subjects)]
+        self.file_list = ["file_{}.fits".format(i + 1) for i in range(number_of_test_subjects)]
         self.focus_values = list(np.linspace(-2000, 2000, number_of_test_subjects))
         fwhm_model = models.Polynomial1D(degree=5)
         fwhm_model.c0.value = 5
@@ -127,8 +122,6 @@ def setUp(self):
 
         self.list_of_fwhm = fwhm_model(self.focus_values)
 
-
-
         for i in range(number_of_test_subjects):
             now = datetime.datetime.now()
             ccd = CCDData(data=np.ones((100, 1000)),
@@ -153,7 +146,7 @@ def setUp(self):
             gaussian = models.Gaussian1D(
                 mean=500,
                 amplitude=600,
-                stddev=self.list_of_fwhm[i]/2.35482004503)
+                stddev=self.list_of_fwhm[i] / 2.35482004503)
 
             for e in range(100):
                 ccd.data[e] = gaussian(range(1000))
@@ -183,10 +176,8 @@ def test_get_focus_data(self):
                                              self.list_of_fwhm)
 
     def test__fit(self):
-        result = self.goodman_focus._fit(df=self.focus_data_frame)
-
-
-        self.assertIsInstance(result, models.Polynomial1D)
+        self.goodman_focus._fit(df=self.focus_data_frame)
+        self.assertIsInstance(self.goodman_focus.polynomial, models.Polynomial1D)
 
     def test__call__(self):
         self.goodman_focus()
@@ -199,7 +190,7 @@ def test__call__Moffat1D(self):
 
     def test__call__with_list(self):
         self.assertIsNone(self.goodman_focus.fwhm)
-        result = self.goodman_focus(files=self.file_list)
+        self.goodman_focus(files=self.file_list)
         self.assertIsNotNone(self.goodman_focus.fwhm)
 
     def test__call__list_file_no_exist(self):
@@ -224,7 +215,6 @@ def setUp(self):
                      'FILTER2': ['NO FILTER'] * 5,
                      'WAVMODE': ['IMAGING'] * 5}
 
-
     def test_imaging_mode(self):
         df = pandas.DataFrame(self.data)
         expected_name = 'IM__Blue__FILTER-X'
@@ -242,7 +232,6 @@ def test_spectroscopy_mode(self):
         self.assertEqual(mode_name, expected_name)
 
 
-
 class DirectoryAndFilesTest(TestCase):
 
     def setUp(self):
@@ -268,7 +257,7 @@ def setUp(self):
 
             ccd.write(os.path.join(os.getcwd(),
                                    'test_dir_no_focus',
-                                   'file_{}.fits'.format(i+1)))
+                                   'file_{}.fits'.format(i + 1)))
 
     def test_directory_does_not_exists(self):
 

goodman_focus/version.py

@@ -1,2 +1,2 @@
 # This is an automatic generated file please do not edit
-__version__ = '2.0.0'
+__version__ = '2.0.2'

setup.cfg

@@ -17,7 +17,7 @@ package_name = goodman_focus
 description = Finds best focus for Goodman HTS based on a series of images obtained with different focus values
 long_description = Standalone app to get the best focus.
 author = Simon Torres
-author_email = storres@ctio.noao.edu
+author_email = simon.torres@noirlab.edu
 license = BSD-3-Clause
 # url =
 edit_on_github = False
@@ -32,4 +32,4 @@ install_requires =
     sphinx
 
 # version should be PEP440 compatible (http://www.python.org/dev/peps/pep-0440)
-version = 2.0.1
+version = 2.0.2