Merge pull request #646 from VChristiaens/master

New features to existing functions and documentation update

docs/source/conf.py

@@ -53,6 +53,8 @@
                           "dollarmath"
                           ]
 
+suppress_warnings = ["mystnb.unknown_mime_type"]
+
 # Add any paths that contain templates here, relative to this directory.
 templates_path = []#['./_templates']
 

docs/source/tutorials/01A_quickstart.ipynb

@@ -105,7 +105,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "```{warning}  \n",
+    "```{important}  \n",
     "Your version of VIP should be >=1.0.3. If not, close the notebook and run for example `pip install vip_hci --upgrade` in your terminal.\n",
     "```"
    ]
@@ -198,7 +198,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "To inspect or retrieve the data, the easiest way is to provide their local path to the relevant functions in VIP - as illustrated below.\n",
+    "To inspect or retrieve the data, the easiest way is to provide their local path to the relevant functions in VIP - as illustrated with the commented code lines below.\n",
     "\n",
     "Here we use the `download_file` utility of `astropy` with the `cache` option - which equivalently stores the local path to the file after downloading the data from an online repository:"
    ]
@@ -260,7 +260,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "```{warning}                                       \n",
+    "```{important}                                       \n",
     "The *convention* in VIP is to always consider you are providing **derotation angles** as input, not **parallactic angles**. Derotation angles are the angles to be applied to the observed images (with positive offsets measured counter-clockwise) in order to align them with north up. These are essentially the *opposite* values of the parallactic angles plus some instrument-specific offsets. The additional offsets can correspond to any internal reflections affecting the orientation of the field of view, or to a correction for the true north orientation which may drift over time and requires dedicated instrumental calibrations. These additional offsets should be checked in the manual of the instrument you are using.\n",
     "```"
    ]
@@ -1351,7 +1351,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "```{warning}                                         \n",
+    "```{important}                                         \n",
     "VIP's convention regarding star centering is:\n",
     "\n",
     "- for **odd** number of pixels along the x and y dimensions: the star centroid should correspond to the central pixel of the image;\n",

docs/source/tutorials/01B_quickstart_with_objects.ipynb

@@ -22,11 +22,11 @@
    "source": [
     "**Table of contents**\n",
     "\n",
-    "* [1.1. Loading ADI data](#1.1.-Loading-ADI-data)\n",
-    "* [1.2. Model PSF subtraction for ADI](#1.2.-Model-PSF-subtraction-for-ADI)\n",
-    "    - [1.2.1. median-ADI](#1.2.1.-median-ADI)\n",
-    "    - [1.2.2. Full-frame PCA](#1.2.2.-Full-frame-PCA)    \n",
-    "* [1.3. Is the point source significant?](#1.3.-Is-it-significant?)"
+    "* [1.1. Loading ADI data](1B.1.-Loading-ADI-data)\n",
+    "* [1.2. Model PSF subtraction for ADI](1B.2.-Model-PSF-subtraction-for-ADI)\n",
+    "    - [1.2.1. median-ADI](1B.2.1.-median-ADI)\n",
+    "    - [1.2.2. Full-frame PCA](1B.2.2.-Full-frame-PCA)    \n",
+    "* [1.3. Is the point source significant?](1B.3.-Is-it-significant?)"
    ]
   },
   {
@@ -108,6 +108,15 @@
     "    raise ValueError(msg)"
    ]
   },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "```{important}\n",
+    "Your version of VIP should be >=1.4.2 to run this notebook. If not, close the notebook and run for example `pip install vip_hci --upgrade` in your terminal. Also make sure the Python version used is 3.10 or higher. \n",
+    "```"
+   ]
+  },
   {
    "cell_type": "markdown",
    "metadata": {},
@@ -119,15 +128,40 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
+    "(1B.1.-Loading-ADI-data)=\n",
     "## 1.1. Loading ADI data"
    ]
   },
   {
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "In order to demonstrate the capabilities of ``VIP``, you can find in the 'dataset' folder of the `VIP_extras` repository: a toy Angular Differential Imaging (ADI) datacube of images and a point spread function (PSF) associated to the NACO instrument of the Very Large Telescope (VLT). This is an L'-band VLT/NACO dataset of beta Pictoris published in [Absil et al. (2013)](https://ui.adsabs.harvard.edu/abs/2013A%26A...559L..12A/abstract), obtained using the Annular Groove Phase Mask (AGPM) Vortex coronagraph. The sequence has been heavily sub-sampled temporarily to make it smaller. The frames were also cropped to the central 101x101 area. In case you want to plug-in your cube just change the path of the following cells.\n",
+    "In order to demonstrate the capabilities of ``VIP``, you can find in the 'dataset' folder of the `VIP_extras` repository: a toy Angular Differential Imaging (ADI) datacube of images and a point spread function (PSF) associated to the NACO instrument of the Very Large Telescope (VLT). This is an L'-band VLT/NACO dataset of beta Pictoris published in [Absil et al. (2013)](https://ui.adsabs.harvard.edu/abs/2013A%26A...559L..12A/abstract), obtained using the Annular Groove Phase Mask (AGPM) Vortex coronagraph. The sequence has been heavily sub-sampled temporarily to make it smaller. The frames were also cropped to the central 101x101 area. In case you want to plug-in your cube just change the path of the following cells."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from astropy.utils.data import download_file\n",
+    "\n",
+    "url_d = \"https://github.com/vortex-exoplanet/VIP_extras/raw/master/datasets\"\n",
+    "f1 = download_file(\"{}/naco_betapic_cube_cen.fits\".format(url_d), cache=True)\n",
+    "f2 = download_file(\"{}/naco_betapic_psf.fits\".format(url_d), cache=True)\n",
+    "f3 = download_file(\"{}/naco_betapic_derot_angles.fits\".format(url_d), cache=True)\n",
     "\n",
+    "# alternatively, for local files simply provide their full or relative path. E.g.:\n",
+    "#f1 = '../datasets/naco_betapic_cube_cen.fits'\n",
+    "#f2 = '../datasets/naco_betapic_psf.fits'\n",
+    "#f3 = '../datasets/naco_betapic_derot_angles.fits'"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
     "Let's first inspect the fits file:"
    ]
   },
@@ -151,7 +185,7 @@
    "source": [
     "from vip_hci.fits import info_fits\n",
     "\n",
-    "info_fits('../datasets/naco_betapic_cube_cen.fits')"
+    "info_fits(f1)"
    ]
   },
   {
@@ -160,20 +194,23 @@
    "source": [
     "The fits file contains an ADI datacube of 61 images 101x101 px across.\n",
     "\n",
-    "The parallactic angle vector for this observation is saved in the `naco_betapic_pa.fits` file. \n",
-    "\n",
-    "Now we load all the data in memory with the `Dataset` object:"
+    "The parallactic angle vector for this observation is saved in the `naco_betapic_pa.fits` file. "
    ]
   },
   {
-   "cell_type": "code",
-   "execution_count": 4,
+   "cell_type": "markdown",
    "metadata": {},
-   "outputs": [],
    "source": [
-    "psfname = '../datasets/naco_betapic_psf.fits'\n",
-    "cubename = '../datasets/naco_betapic_cube_cen.fits'\n",
-    "angname = '../datasets/naco_betapic_pa.fits'"
+    "```{important}                                       \n",
+    "The *convention* in VIP is to always consider you are providing **derotation angles** as input, not **parallactic angles**. Derotation angles are the angles to be applied to the observed images (with positive offsets measured counter-clockwise) in order to align them with north up. These are essentially the *opposite* values of the parallactic angles plus some instrument-specific offsets. The additional offsets can correspond to any internal reflections affecting the orientation of the field of view, or to a correction for the true north orientation which may drift over time and requires dedicated instrumental calibrations. These additional offsets should be checked in the manual of the instrument you are using.\n",
+    "```"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Now we load all the data in memory with the `Dataset` object:"
    ]
   },
   {
@@ -193,7 +230,7 @@
    ],
    "source": [
     "from vip_hci.objects import Dataset\n",
-    "betapic = Dataset(cube=cubename, psf=psfname, angles=angname)"
+    "betapic = Dataset(cube=f1, psf=f2, angles=f3)"
    ]
   },
   {
@@ -1156,12 +1193,9 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "---\n",
-    "**NOTE**\n",
-    "\n",
+    "```{note}\n",
     "The `normalize_psf` function performs internally a fine centering of the PSF. The input PSF should nevertheless already be centered within a few pixel accuracy for the fine centering to work.\n",
-    "\n",
-    "---"
+    "```"
    ]
   },
   {
@@ -1223,7 +1257,8 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "## 1.2. Model PSF subtraction for ADI"
+    "(1B.2.-PSF-subtraction)=\n",
+    "## 1.2. PSF subtraction"
    ]
   },
   {
@@ -1237,15 +1272,19 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "**IMPORTANT** - VIP's convention regarding centering is:\n",
-    "- for **odd** number of pixels along the x and y directions: the star is placed on the central pixel;\n",
-    "- for **even** number of pixels: the star is placed on coordinates (dim/2,dim/2) in 0-based indexing."
+    "```{important}                                 \n",
+    "VIP's convention regarding star centering is:\n",
+    "\n",
+    "- for **odd** number of pixels along the x and y dimensions: the star centroid should correspond to the central pixel of the image;\n",
+    "- for **even** number of pixels: the star centroid should be located on coordinates (dim/2,dim/2) in 0-based indexing (i.e. on the top-right pixel of the 4 central pixels).\n",
+    "```"
    ]
   },
   {
    "cell_type": "markdown",
    "metadata": {},
    "source": [
+    "(1B.2.1.-median-ADI)=\n",
     "### 1.2.1. median-ADI"
    ]
   },
@@ -1365,6 +1404,16 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
+    "```{note}\n",
+    "As per usual convention, north is up and east is left in all images shown in the tutorial. Note that for the images of some instruments, te x-axis may be flipped (i.e. east would be right instead of left). This is for example the case for Gemini/NICI or VLT/ERIS-NIX. Double-check your instrument manual.\n",
+    "```"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "(1B.2.2.-Full-frame-PCA)=\n",
     "### 1.2.2. Full-frame PCA"
    ]
   },
@@ -1470,6 +1519,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
+    "(1B.3.-Is-it-significant?)=\n",
     "## 1.3. Is it significant?"
    ]
   },
@@ -1575,6 +1625,13 @@
    "toc_position": {},
    "toc_section_display": true,
    "toc_window_display": true
+  },
+  "widgets": {
+   "application/vnd.jupyter.widget-state+json": {
+    "state": {},
+    "version_major": 2,
+    "version_minor": 0
+   }
   }
  },
  "nbformat": 4,

vip_hci/fm/fakecomp.py

@@ -793,7 +793,7 @@ def psf_norm_2d(psf, fwhm, threshold, mask_core, full_output, verbose):
             restemp = psf_norm_2d(array[fr], fwhm[fr], threshold, mask_core,
                                   True, False)
             array_out.append(restemp[0])
-            fwhm_flux[fr] = restemp[1]
+            fwhm_flux[fr] = restemp[1][0]
 
         array_out = np.array(array_out)
         if verbose:

vip_hci/fm/negfc_fmerit.py

@@ -305,7 +305,7 @@ def chisquare(
                 msg += "calculation. Consider larger input images."
                 raise ValueError(msg)
 
-            subim = frame_crop(frpca, crop_sz, cenxy=(xround, yround),
+            subim = frame_crop(frpca, crop_sz, xy=(xround, yround),
                                force=True, verbose=False)
             H = hessian(subim)
             dets = np.zeros([ndet, ndet])

vip_hci/invprob/fmmf.py

@@ -460,15 +460,15 @@ def _snr_contrast_esti(
                 psfm = frame_crop(
                     psfm_temp[j],
                     crop,
-                    cenxy=[int(psfm_temp.shape[-1] / 2),
+                    xy=[int(psfm_temp.shape[-1] / 2),
                            int(psfm_temp.shape[-1] / 2)],
                     verbose=False,
                 )
 
             num.append(
                 np.multiply(
                     frame_crop(
-                        mcube[j], crop, cenxy=[
+                        mcube[j], crop, xy=[
                             poscentx, poscenty], verbose=False
                     ),
                     psfm,

vip_hci/metrics/contrcurve.py

@@ -703,12 +703,12 @@ def throughput(
                 raise TypeError(msg)
             if psf_template.ndim != 3:
                 raise TypeError("Template PSF is not a frame, 3d array")
-            if "scale_list" not in algo_dict:
-                raise ValueError("Vector of wavelength not found")
-            else:
+            if "scale_list" in algo_dict:
+            #     raise ValueError("Vector of wavelength not found")
+            # else:
                 if algo_dict["scale_list"].shape[0] != array.shape[0]:
                     raise TypeError("Input wavelength vector has wrong length")
-                if isinstance(fwhm, float) or isinstance(fwhm, int):
+                if np.isscalar(fwhm):
                     maxfcsep = int((array.shape[2] / 2.0) / fwhm) - 1
                 else:
                     maxfcsep = int((array.shape[2] / 2.0) / np.amin(fwhm)) - 1

vip_hci/metrics/snr_source.py

@@ -452,7 +452,8 @@ def snr(array, source_xy, fwhm, full_output=False, array2=None, use2alone=False,
         return snr_vale
 
 
-def significance(snr, rad, fwhm, student_to_gauss=True, verbose=True):
+def significance(snr, rad, fwhm, n_ap=None, student_to_gauss=True,
+                 verbose=True):
     """Convert a S/N ratio (measured as in [MAW14]_) into a Gaussian\
     significance (n-sigma) with equivalent false alarm probability for a\
     point-source detection measured at a given separation, or the opposite.
@@ -467,6 +468,11 @@ def significance(snr, rad, fwhm, student_to_gauss=True, verbose=True):
         to each snr measurement.
     fwhm : float
         Full Width Half Maximum of the PSF.
+    n_ap : int, optional
+        Number of independent samples available. If provided, will take
+        precedence over using input rad and FWHM. This can be useful when the
+        limited azimuthal coverage is available (e.g. when using APP or 4QPM
+        coronagraphs).
     student_to_gauss : bool, optional
         Whether the conversion is from Student SNR to Gaussian significance
         (True), or the opposite (False).
@@ -477,9 +483,12 @@ def significance(snr, rad, fwhm, student_to_gauss=True, verbose=True):
         Equivalent Gaussian significance [student_to_gauss=True] or equivalent
         Student S/N ratio [student_to_gauss=False].
     """
+    if n_ap is None:
+        n_ap = (rad/fwhm)*2*np.pi-2
+
     if student_to_gauss:
-        sig = norm.ppf(t.cdf(snr, (rad/fwhm)*2*np.pi-2))
-        if t.cdf(snr, (rad/fwhm)*2*np.pi-2) == 1.0:
+        sig = norm.ppf(t.cdf(snr, n_ap))
+        if t.cdf(snr, n_ap) == 1.0:
             print("Warning high S/N! cdf>0.9999999999999999 is rounded to 1")
             msg = "Returning 8.2 sigma, but quote significance > 8.2 sigma."
             print(msg)
@@ -490,7 +499,7 @@ def significance(snr, rad, fwhm, student_to_gauss=True, verbose=True):
             msg += r"Gaussian false alarm probability."
             print(msg.format(rad, rad/fwhm, snr, sig))
     else:
-        sig = t.ppf(norm.cdf(snr), (rad/fwhm)*2*np.pi-2)
+        sig = t.ppf(norm.cdf(snr), n_ap)
         if verbose:
             msg = r"At a separation of {:.1f} px ({:.1f} FWHM), a {:.1f}-sigma "
             msg += r"detection in terms of Gaussian false alarm probability "