Compatibility of the MCMCsamplingModule with the Hessian minimization (…

…#409)

* Testing MCMC with Hessian

* Print output

* Hessian variance

* Updated documentation

* Updated print output

* Updated documentation

* Indentation fix

docs/coding.rst

@@ -5,6 +5,8 @@ Coding a New Module
 
 .. _constructor:
 
+There are three different types of pipeline modules: :class:`~pynpoint.core.processing.ReadingModule`, :class:`~pynpoint.core.processing.WritingModule`, and :class:`~pynpoint.core.processing.ProcessingModule`. The concept is similar for these three modules so here we will explain only how to code a processing module.
+
 Class Constructor
 -----------------
 

docs/python.rst

@@ -26,8 +26,6 @@ The modular architecture of PynPoint allows for easy implementation of new pipel
 
    <a href="https://docs.python.org/3/library/abc.html" target="_blank">Abstract classes</a>
 
-There are three different types of pipeline modules: :class:`~pynpoint.core.processing.ReadingModule`, :class:`~pynpoint.core.processing.WritingModule`, and :class:`~pynpoint.core.processing.ProcessingModule`. The concept is similar for the three types of modules so here we will explain only how to code a processing module.
-
 .. _conventions:
 
 Conventions
@@ -52,4 +50,24 @@ Before we start writing a new PynPoint module, please take notice of the followi
 
    <a href="https://pypi.org/project/pycodestyle/" target="_blank">pycodestyle</a>
 
-Now we are ready to code!
+Unit tests
+----------
+
+PynPoint is a robust pipeline package with 95% of the code covered by |unittest|. Testing of the package is done by running ``make test`` in the cloned repository. This requires the installation of:
+
+   * |pytest|
+   * |pytest-cov|
+
+The unit tests ensure that the output from existing functionalities will not change whenever new code. With these things in mind, we are now ready to code!
+
+.. |unittest| raw:: html
+
+   <a href="https://docs.python.org/3/library/unittest.html" target="_blank">unit tests</a>
+
+.. |pytest| raw:: html
+
+   <a href="https://docs.pytest.org/en/latest/getting-started.html" target="_blank">pytest</a>
+
+.. |pytest-cov| raw:: html
+
+   <a href="https://pytest-cov.readthedocs.io/en/latest/readme.html" target="_blank">pytest-cov</a>

pull_request_template.md

@@ -3,7 +3,9 @@ Thank you for your contribution to the PynPoint repo! Before submitting this PR,
 - [ ] To read the documentation page on the [Python guidelines](https://pynpoint.readthedocs.io/en/latest/python.html).
 - [ ] That your branch of the PR is up to date with the master branch of the upstream PynPoint repo.
 - [ ] To run both `pycodestyle` and `pylint` on the code that has been added and/or changed.
-- [ ] That the documentation is successfully build after running `make docs` in your local repo folder.
-- [ ] That all unit tests are finishing after running `make test` in your local repo folder.
+- [ ] That the documentation is successfully build after running `make docs` in your local repo folder. This requires the installation of `sphinx` and `sphinx_rtd_theme`.
+- [ ] That all unit tests are finishing after running `make test` in your local repo folder. This requires the installation of `pytest` and `pytest-cov`.
 - [ ] To add unit tests in case there are new pipeline modules and/or functionalities added.
-- [ ] Make sure that only text files have been added and changed in the commits of the PR. Binary files will clutter up the repo because even after removing such files they will remain in the repo history.
+- [ ] That only text files have been added and changed in the commits of the PR. Binary files will clutter up the repo because even after removing such files they will remain in the repo history.
+- [ ] To add and/or update the docstrings.
+- [ ] To add and/or update the typehints and typechecks.

pynpoint/processing/fluxposition.py

@@ -19,7 +19,7 @@
 from photutils.aperture import Aperture
 
 from pynpoint.core.processing import ProcessingModule
-from pynpoint.util.analysis import fake_planet, merit_function, false_alarm, gaussian_noise
+from pynpoint.util.analysis import fake_planet, merit_function, false_alarm, pixel_variance
 from pynpoint.util.image import create_mask, polar_to_cartesian, cartesian_to_polar, \
                                 center_subpixel, rotate_coordinates
 from pynpoint.util.mcmc import lnprob
@@ -158,8 +158,8 @@ def run(self) -> None:
 
 class SimplexMinimizationModule(ProcessingModule):
     """
-    Pipeline module to measure the flux and position of a planet by injecting negative fake planets
-    and minimizing a figure of merit.
+    Pipeline module to retrieve the contrast and position of a planet by injecting negative
+    artificial planets and using a downhill simplex method.
     """
 
     __author__ = 'Tomas Stolker'
@@ -213,11 +213,11 @@ def __init__(self,
             Additional scaling factor of the planet flux (e.g., to correct for a neutral density
             filter). Should be negative in order to inject negative fake planets.
         merit : str
-            Figure of merit for the minimization. Can be 'hessian', to minimize the sum of the
-            absolute values of the determinant of the Hessian matrix, or 'poisson', to minimize the
-            sum of the absolute pixel values, assuming a Poisson distribution for the noise
-            (Wertz et al. 2017), or 'gaussian', to minimize the ratio of the squared pixel values
-            and the variance of the pixels within an annulus but excluding the aperture area.
+            Figure of merit for the minimization ('hessian', 'gaussian', or 'poisson'). Either the
+            determinant of the Hessian matrix is minimized ('hessian') or the flux of each pixel
+            ('gaussian' or 'poisson'). For the latter case, the estimate noise is assumed to follow
+            a Poisson (see Wertz et al. 2017) or Gaussian distribution (see Wertz et al. 2017 and
+            Stolker et al. 2020).
         aperture : float
             Aperture radius (arcsec) at the position specified at *position*.
         sigma : float
@@ -309,8 +309,8 @@ def __init__(self,
     def run(self) -> None:
         """
         Run method of the module. The position and contrast of a planet is measured by injecting
-        negative copies of the PSF template and applying a simplex method (Nelder-Mead) for
-        minimization of a figure of merit at the planet location.
+        negative copies of the PSF template and applying a downhill simplex method (Nelder-Mead)
+        for minimization of a figure of merit at the planet location.
 
         Returns
         -------
@@ -351,7 +351,7 @@ def _objective(arg: np.ndarray,
                        count: int,
                        n_components: int,
                        sklearn_pca: Optional[PCA],
-                       noise: Optional[float]) -> float:
+                       var_noise: Optional[float]) -> float:
 
             pos_y = arg[0]
             pos_x = arg[1]
@@ -406,7 +406,7 @@ def _objective(arg: np.ndarray,
                                         merit=self.m_merit,
                                         aperture=aperture,
                                         sigma=self.m_sigma,
-                                        noise=noise)
+                                        var_noise=var_noise)
 
             position = rotate_coordinates(center, (pos_y, pos_x), -self.m_extra_rot)
 
@@ -464,21 +464,23 @@ def _objective(arg: np.ndarray,
 
                 sklearn_pca.components_ = q_ortho.T
 
-            if self.m_merit in ('poisson', 'hessian'):
-                noise = None
+            if self.m_merit == 'poisson':
+                var_noise = None
 
-            elif self.m_merit == 'gaussian':
-                noise = gaussian_noise(images=images,
-                                       parang=parang,
-                                       cent_size=self.m_cent_size,
-                                       edge_size=self.m_edge_size,
-                                       pca_number=n_components,
-                                       residuals=self.m_residuals,
-                                       aperture=aperture)
+            elif self.m_merit in ['gaussian', 'hessian']:
+                var_noise = pixel_variance(var_type=self.m_merit,
+                                           images=images,
+                                           parang=parang,
+                                           cent_size=self.m_cent_size,
+                                           edge_size=self.m_edge_size,
+                                           pca_number=n_components,
+                                           residuals=self.m_residuals,
+                                           aperture=aperture,
+                                           sigma=self.m_sigma)
 
             minimize(fun=_objective,
                      x0=np.array([pos_init[0], pos_init[1], self.m_magnitude]),
-                     args=(i, n_components, sklearn_pca, noise),
+                     args=(i, n_components, sklearn_pca, var_noise),
                      method='Nelder-Mead',
                      tol=None,
                      options={'xatol': self.m_tolerance, 'fatol': float('inf')})
@@ -620,7 +622,7 @@ def _snr_optimize(arg: np.ndarray) -> float:
                                            size=self.m_aperture,
                                            ignore=self.m_ignore)
 
-            return -snr
+            return -1.*snr
 
         pixscale = self.m_image_in_port.get_attribute('PIXSCALE')
         self.m_aperture /= pixscale
@@ -741,12 +743,11 @@ def __init__(self,
         extra_rot : float
             Additional rotation angle of the images (deg).
         merit : str
-            Figure of merit that is used for the likelihood function ('gaussian' or 'poisson').
-            Pixels are assumed to be independent measurements which are expected to be equal to
-            zero in case the best-fit negative PSF template is injected. With 'gaussian', the
-            variance is estimated from the pixel values within an annulus at the separation of
-            the aperture (but excluding the pixels within the aperture). With 'poisson', a Poisson
-            distribution is assumed for the variance of each pixel value (see Wertz et al. 2017).
+            Figure of merit for the minimization ('hessian', 'gaussian', or 'poisson'). Either the
+            determinant of the Hessian matrix is minimized ('hessian') or the flux of each pixel
+            ('gaussian' or 'poisson'). For the latter case, the estimate noise is assumed to follow
+            a Poisson (see Wertz et al. 2017) or Gaussian distribution (see Wertz et al. 2017 and
+            Stolker et al. 2020).
         residuals : str
             Method used for combining the residuals ('mean', 'median', 'weighted', or 'clipped').
 
@@ -854,17 +855,28 @@ def run(self) -> None:
         elif isinstance(self.m_aperture, tuple):
             aperture = (self.m_aperture[1], self.m_aperture[0], self.m_aperture[2]/pixscale)
 
+        print(f'Aperture position [x, y]: [{aperture[1]}, {aperture[0]}]')
+        print(f'Aperture radius (pix): {aperture[2]:.2f}')
+        print(f'Figure of merit: {self.m_merit}')
+
         if self.m_merit == 'poisson':
-            noise = None
+            var_noise = None
+
+        elif self.m_merit in ['gaussian', 'hessian']:
+            var_noise = pixel_variance(var_type=self.m_merit,
+                                       images=images,
+                                       parang=parang,
+                                       cent_size=self.m_mask[0],
+                                       edge_size=self.m_mask[1],
+                                       pca_number=self.m_pca_number,
+                                       residuals=self.m_residuals,
+                                       aperture=aperture,
+                                       sigma=0.)
 
-        elif self.m_merit == 'gaussian':
-            noise = gaussian_noise(images=images,
-                                   parang=parang,
-                                   cent_size=self.m_mask[0],
-                                   edge_size=self.m_mask[1],
-                                   pca_number=self.m_pca_number,
-                                   residuals=self.m_residuals,
-                                   aperture=aperture)
+            if self.m_merit == 'gaussian':
+                print(f'Gaussian standard deviation (counts): {np.sqrt(var_noise):.2e}')
+            if self.m_merit == 'hessian':
+                print(f'Hessian standard deviation: {np.sqrt(var_noise):.2e}')
 
         initial = np.zeros((self.m_nwalkers, ndim))
 
@@ -892,7 +904,7 @@ def run(self) -> None:
                                                    indices,
                                                    self.m_merit,
                                                    self.m_residuals,
-                                                   noise]))
+                                                   var_noise]))
 
             sampler.run_mcmc(initial, self.m_nsteps, progress=True)
 
@@ -1089,12 +1101,11 @@ def __init__(self,
             Additional scaling factor of the planet flux (e.g., to correct for a neutral density
             filter). Should be a positive value.
         merit : str
-            Figure of merit that is used for the likelihood function ('gaussian' or 'poisson').
-            Pixels are assumed to be independent measurements which are expected to be equal to
-            zero in case the best-fit negative PSF template is injected. With 'gaussian', the
-            variance is estimated from the pixel values within an annulus at the separation of
-            the aperture (but excluding the pixels within the aperture). With 'poisson', a Poisson
-            distribution is assumed for the variance of each pixel value (see Wertz et al. 2017).
+            Figure of merit for the minimization ('hessian', 'gaussian', or 'poisson'). Either the
+            determinant of the Hessian matrix is minimized ('hessian') or the flux of each pixel
+            ('gaussian' or 'poisson'). For the latter case, the estimate noise is assumed to follow
+            a Poisson (see Wertz et al. 2017) or Gaussian distribution (see Wertz et al. 2017 and
+            Stolker et al. 2020).
         aperture : float
             Aperture radius (arcsec) that is used for measuring the figure of merit.
         tolerance : float

pynpoint/processing/resizing.py

@@ -83,7 +83,12 @@ def run(self) -> None:
 
         # Convert size parameter from arcseconds to pixels
         pixscale = self.m_image_in_port.get_attribute('PIXSCALE')
+        print(f'New image size (arcsec) = {self.m_size}')
         self.m_size = int(math.ceil(self.m_size / pixscale))
+        print(f'New image size (pixels) = {self.m_size}')
+
+        if self.m_center is not None:
+            print(f'New image center (x, y) = {self.m_center}')
 
         # Crop images chunk by chunk
         start_time = time.time()
@@ -100,7 +105,7 @@ def run(self) -> None:
             self.m_image_out_port.append(images, data_dim=3)
 
         # Save history and copy attributes
-        history = f'image size [pix] = {self.m_size}'
+        history = f'image size (pix) = {self.m_size}'
         self.m_image_out_port.add_history('CropImagesModule', history)
         self.m_image_out_port.copy_attributes(self.m_image_in_port)
         self.m_image_out_port.close_port()

pynpoint/util/analysis.py

@@ -258,7 +258,7 @@ def merit_function(residuals: np.ndarray,
                    merit: str,
                    aperture: Tuple[int, int, float],
                    sigma: float,
-                   noise: Optional[float]) -> float:
+                   var_noise: Optional[float]) -> float:
     """
     Function to calculate the figure of merit at a given position in the image residuals.
 
@@ -273,13 +273,13 @@ def merit_function(residuals: np.ndarray,
     sigma : float
         Standard deviation (pix) of the Gaussian kernel which is used to smooth the residuals
         before the chi-square is calculated.
-    noise : float, None
-        Variance of the noise which is required when `merit` is set to 'gaussian'.
+    var_noise : float, None
+        Variance of the noise which is required when `merit` is set to 'gaussian' or 'hessian'.
 
     Returns
     -------
     float
-        Chi-square ('poisson' and 'gaussian') or sum of the absolute values ('hessian').
+        Chi-square value.
     """
 
     rr_grid = pixel_distance(im_shape=residuals.shape,
@@ -289,8 +289,6 @@ def merit_function(residuals: np.ndarray,
 
     if merit == 'hessian':
 
-        # This is not the chi-square but simply the sum of the absolute values
-
         hessian_rr, hessian_rc, hessian_cc = hessian_matrix(image=residuals,
                                                             sigma=sigma,
                                                             mode='constant',
@@ -299,7 +297,7 @@ def merit_function(residuals: np.ndarray,
 
         hes_det = (hessian_rr*hessian_cc) - (hessian_rc*hessian_rc)
 
-        chi_square = np.sum(np.abs(hes_det[indices]))
+        chi_square = np.sum(hes_det[indices]**2)/var_noise
 
     elif merit == 'poisson':
 
@@ -310,7 +308,7 @@ def merit_function(residuals: np.ndarray,
 
     elif merit == 'gaussian':
 
-        chi_square = np.sum(residuals[indices]**2)/noise
+        chi_square = np.sum(residuals[indices]**2)/var_noise
 
     else:
 
@@ -322,20 +320,24 @@ def merit_function(residuals: np.ndarray,
 
 
 @typechecked
-def gaussian_noise(images: np.ndarray,
+def pixel_variance(var_type: str,
+                   images: np.ndarray,
                    parang: np.ndarray,
                    cent_size: Optional[float],
                    edge_size: Optional[float],
                    pca_number: int,
                    residuals: str,
-                   aperture: Tuple[int, int, float]) -> float:
+                   aperture: Tuple[int, int, float],
+                   sigma: float) -> float:
     """
     Function to calculate the variance of the noise. After the PSF subtraction, images are rotated
     in opposite direction of the regular derotation, therefore dispersing any companion or disk
     signal. The noise is measured within an annulus.
 
     Parameters
     ----------
+    var_type : str
+        Variance type ('gaussian' or 'hessian').
     images : numpy.ndarray
         Input images (3D).
     parang : numpy.ndarray
@@ -351,11 +353,14 @@ def gaussian_noise(images: np.ndarray,
         Method for combining the residuals ('mean', 'median', 'weighted', or 'clipped').
     aperture : tuple(int, int, float)
         Aperture position (y, x) and radius (pix).
+    sigma : float, None
+        Standard deviation (pix) of the Gaussian kernel which is used to smooth the images.
 
     Returns
     -------
     float
-        Variance of the pixel values.
+        Variance of the pixel values. Either the variance of the pixel values ('gaussian') or
+        the variance of the determinant of the Hessian ('hessian').
     """
 
     mask = create_mask(images.shape[-2:], (cent_size, edge_size))
@@ -366,6 +371,18 @@ def gaussian_noise(images: np.ndarray,
 
     sep_ang = cartesian_to_polar(center_subpixel(res_noise), aperture[0], aperture[1])
 
-    selected = select_annulus(res_noise[0, ], sep_ang[0]-aperture[2], sep_ang[0]+aperture[2])
+    if var_type == 'gaussian':
+        selected = select_annulus(res_noise[0, ], sep_ang[0]-aperture[2], sep_ang[0]+aperture[2])
+
+    elif var_type == 'hessian':
+        hessian_rr, hessian_rc, hessian_cc = hessian_matrix(image=res_noise[0, ],
+                                                            sigma=sigma,
+                                                            mode='constant',
+                                                            cval=0.,
+                                                            order='rc')
+
+        hes_det = (hessian_rr*hessian_cc) - (hessian_rc*hessian_rc)
+
+        selected = select_annulus(hes_det, sep_ang[0]-aperture[2], sep_ang[0]+aperture[2])
 
     return float(np.var(selected))