bispect

Bispectrum extraction routines for aperture-masking interferometry.

The module samples matched Fourier peaks, estimates calibrated visibility and bispectrum statistics, and assembles AMICAL interferometric observables.

extract_bs(cube, filename, maskname, filtname=None, targetname=None, instrum=None, bs_multi_tri=False, peakmethod='gauss', hole_diam=0.8, cutoff=0.0001, fw_splodge=0.7, naive_err=False, n_wl=3, n_blocks=0, theta_detector=0, scaling_uv=1, i_wl=None, unbias_v2=True, compute_cp_cov=True, expert_plot=False, save_to=None, verbose=False, display=True)

Extract calibrated aperture-masking interferometric observables from a data cube.

Parameters:

Name	Type	Description	Default
cube	(ndarray, shape(n_frames, ny, nx))	Cleaned image cube ready for non-redundant-mask extraction.	required
filename	str or path - like	FITS file from which observation headers are read.	required
maskname	str	Aperture-mask identifier.	required
filtname	str	Filter-name fallback when the FITS header lacks FILTER.	None
targetname	str	Target-name fallback when the FITS header lacks OBJECT.	None
instrum	str	Instrument-name fallback when the FITS header lacks INSTRUME.	None
bs_multi_tri	bool	Whether to compute bispectra from multiple pixel closure triangles.	False
peakmethod	(fft, square, unique, gauss)	Method used to sample Fourier splodges.	"fft"
hole_diam	float	Aperture-hole diameter in metres.	0.8
cutoff	float	Minimum matched-filter weight retained for a peak pixel.	1e-4
fw_splodge	float	Relative Fourier-splodge size and Gaussian-method FWHM scale.	0.7
naive_err	bool	Whether to use frame-wise standard deviations instead of covariance-based uncertainties.	False
n_wl	int	Number of wavelengths used to sample filter bandwidth.	3
n_blocks	int	Number of frame blocks used for uncertainty estimation; 0 uses all frames.	0
theta_detector	float	Mask rotation relative to the detector in degrees.	0
scaling_uv	float	Multiplicative scale applied to aperture-mask coordinates.	1
i_wl	int or sequence of int	SPHERE-IFS spectral-channel selector.	None
unbias_v2	bool	Whether to subtract the correlated-noise bias from squared visibilities.	True
compute_cp_cov	bool	Whether to compute the closure-phase covariance matrix.	True
expert_plot	bool	Whether to show expert-level diagnostic plots.	False
save_to	str or path - like	Directory in which diagnostic figures are saved.	None
verbose	bool	Whether to print processing progress.	False
display	bool	Whether to display diagnostic figures.	True

Returns:

Type	Description
object or None	AMI observables including squared visibilities, closure phases in degrees, uncertainties, baseline coordinates in metres, wavelength in metres, and mask, metadata, and covariance information. Returns None if the mask cannot be loaded or the matched filter cannot be constructed.

Raises:

Type	Description
ValueError	If SPHERE-IFS data are processed without i_wl.

Source code in src/amical/mf_pipeline/bispect.py

def extract_bs(
    cube,
    filename,
    maskname,
    filtname=None,
    targetname=None,
    instrum=None,
    bs_multi_tri=False,
    peakmethod="gauss",
    hole_diam=0.8,
    cutoff=1e-4,
    fw_splodge=0.7,
    naive_err=False,
    n_wl=3,
    n_blocks=0,
    theta_detector=0,
    scaling_uv=1,
    i_wl=None,
    unbias_v2=True,
    compute_cp_cov=True,
    expert_plot=False,
    save_to=None,
    verbose=False,
    display=True,
):
    """Extract calibrated aperture-masking interferometric observables from a data cube.

    Parameters
    ----------
    cube : numpy.ndarray, shape (n_frames, ny, nx)
        Cleaned image cube ready for non-redundant-mask extraction.
    filename : str or path-like
        FITS file from which observation headers are read.
    maskname : str
        Aperture-mask identifier.
    filtname : str, optional
        Filter-name fallback when the FITS header lacks ``FILTER``.
    targetname : str, optional
        Target-name fallback when the FITS header lacks ``OBJECT``.
    instrum : str, optional
        Instrument-name fallback when the FITS header lacks ``INSTRUME``.
    bs_multi_tri : bool, default=False
        Whether to compute bispectra from multiple pixel closure triangles.
    peakmethod : {"fft", "square", "unique", "gauss"}, default="gauss"
        Method used to sample Fourier splodges.
    hole_diam : float, default=0.8
        Aperture-hole diameter in metres.
    cutoff : float, default=1e-4
        Minimum matched-filter weight retained for a peak pixel.
    fw_splodge : float, default=0.7
        Relative Fourier-splodge size and Gaussian-method FWHM scale.
    naive_err : bool, default=False
        Whether to use frame-wise standard deviations instead of covariance-based
        uncertainties.
    n_wl : int, default=3
        Number of wavelengths used to sample filter bandwidth.
    n_blocks : int, default=0
        Number of frame blocks used for uncertainty estimation; 0 uses all frames.
    theta_detector : float, default=0
        Mask rotation relative to the detector in degrees.
    scaling_uv : float, default=1
        Multiplicative scale applied to aperture-mask coordinates.
    i_wl : int or sequence of int, optional
        SPHERE-IFS spectral-channel selector.
    unbias_v2 : bool, default=True
        Whether to subtract the correlated-noise bias from squared visibilities.
    compute_cp_cov : bool, default=True
        Whether to compute the closure-phase covariance matrix.
    expert_plot : bool, default=False
        Whether to show expert-level diagnostic plots.
    save_to : str or path-like, optional
        Directory in which diagnostic figures are saved.
    verbose : bool, default=False
        Whether to print processing progress.
    display : bool, default=True
        Whether to display diagnostic figures.

    Returns
    -------
    object or None
        AMI observables including squared visibilities, closure phases in degrees,
        uncertainties, baseline coordinates in metres, wavelength in metres, and
        mask, metadata, and covariance information. Returns ``None`` if the mask
        cannot be loaded or the matched filter cannot be constructed.

    Raises
    ------
    ValueError
        If SPHERE-IFS data are processed without ``i_wl``."""
    from astropy.io import fits

    if verbose:
        rprint("[cyan]\n-- Starting extraction of observables --")
    start_time = time.time()

    if save_to is not None and not display:
        warnings.warn(
            "Diagnostic figures are only saved when display=True: "
            "`save_to` is ignored.",
            stacklevel=2,
        )
        save_to = None

    if save_to is not None:
        if not os.path.exists(save_to):
            os.mkdir(save_to)

    with fits.open(filename) as hdu:
        hdr = hdu[0].header
        try:
            sci_hdr = hdu["SCI"].header
        except KeyError:
            sci_hdr = None

    infos = _check_input_infos(
        hdr, targetname=targetname, filtname=filtname, instrum=instrum, verbose=False
    )

    if infos.instrument == "SPHERE-IFS":
        if i_wl is None:
            raise ValueError(
                "Your file seems to be obtained with an IFU instrument: spectral "
                "channel index `i_wl` must be specified."
            )

    if "INSTRUME" not in hdr.keys():
        hdr["INSTRUME"] = infos["instrument"]
    # 1. Open the data cube and perform a series of roll (both axis) to avoid
    # grid artefact (negative fft values).
    # ------------------------------------------------------------------------
    ft_arr, n_ps, npix = _construct_ft_arr(cube)

    # Number of aperture in the mask
    try:
        n_holes = len(get_mask(infos.instrument, maskname))
    except TypeError:
        return None

    # 2. Determine the number of different baselines (bl), bispectrums (bs) or
    # covariance matrices (cov) and associates each holes as couple for bl or
    # triplet for bs (or cp) using compute_index_mask function (see ami_function.py).
    # ------------------------------------------------------------------------
    index_mask = compute_index_mask(n_holes)

    n_baselines = index_mask.n_baselines

    closing_tri = _format_closing_triangle(index_mask)

    # 3. Compute the match filter mf
    # ------------------------------------------------------------------------
    mf = make_mf(
        maskname,
        infos.instrument,
        infos.filtname,
        npix,
        peakmethod=peakmethod,
        fw_splodge=fw_splodge,
        n_wl=n_wl,
        cutoff=cutoff,
        hole_diam=hole_diam,
        scaling=scaling_uv,
        theta_detector=theta_detector,
        i_wl=i_wl,
        display=display,
        save_to=save_to,
        filename=filename,
        index_mask=index_mask,
    )

    ifig = 2
    if save_to is not None:
        import matplotlib.pyplot as plt

        figname = os.path.join(save_to, Path(filename).stem)
        plt.savefig(f"{figname}_{ifig}.pdf")

    ifig += 1

    if mf is None:
        return None

    # We store the principal results in the new dictionnary to be save at the end
    obs_result = {"u": mf.u, "v": mf.v, "wl": mf.wl, "e_wl": mf.e_wl}

    # 4. Compute indices for the multiple triangle technique (tri_pix function)
    # -------------------------------------------------------------------------
    l_B = np.sqrt(mf.u**2 + mf.v**2)  # Length of different bl [m]
    minbl = np.min(l_B)

    if n_holes >= 15:
        sampledisk_r = minbl / 2 / mf.wl * mf.pixelSize * npix * 0.9
    else:
        sampledisk_r = minbl / 2 / mf.wl * mf.pixelSize * npix * fw_splodge

    if bs_multi_tri:
        closing_tri_pix = tri_pix(npix, sampledisk_r, display=display, verbose=verbose)
    else:
        closing_tri_pix = None

    # 5. Display the power spectrum of the first frame to check the computed
    # positions of the peaks.
    # ------------------------------------------------------------------------
    if display:
        _show_complex_ps(ft_arr)
        if save_to is not None:
            plt.savefig(f"{figname}_{ifig}.pdf")
        ifig += 1

        _show_peak_position(ft_arr, n_baselines, mf, maskname, peakmethod)
        if save_to is not None:
            plt.savefig(f"{figname}_{ifig}.pdf")
        ifig += 1

    if verbose:
        print(f"\nFilename: {filename}")
        print(f"# of frames = {n_ps}")

    n_blocks = _set_good_nblocks(n_blocks, n_ps)

    # 6. Extract the complex quantities from the fft_arr (complex vis, bispectrum,
    # phase, etc.)
    # ------------------------------------------------------------------------
    fringe_peak = give_peak_info2d(mf, n_baselines, npix, npix)

    if verbose:
        print("\nCalculating V^2 and BS...")

    complex_bs = _compute_complex_bs(
        ft_arr,
        index_mask,
        fringe_peak,
        mf,
        dark_ps=None,
        closing_tri_pix=closing_tri_pix,
        bs_multi_tri=bs_multi_tri,
        verbose=verbose,
    )

    cvis_arr = complex_bs["vis_arr"]["complex"]
    v2_arr = complex_bs["vis_arr"]["squared"]
    bs_arr = complex_bs["bs_arr"]
    fluxes = complex_bs["fluxes"]

    # 7. Compute correlated noise and bias at the peak position
    # ---------------------------------------------------------
    bias, dark_bias, autocor_noise = _compute_corr_noise(
        complex_bs, ft_arr, fringe_peak
    )

    v2_arr_unbiased, bias_arr = _unbias_v2_arr(
        v2_arr, npix, fringe_peak, bias, dark_bias, autocor_noise, unbias=unbias_v2
    )

    # 8. Turn Arrays into means and covariance matrices
    # -------------------------------------------------
    v2_quantities = _compute_v2_quantities(v2_arr_unbiased, bias_arr, n_blocks)

    bs_quantities = _compute_bs_quantities(
        bs_arr, v2_quantities["v2"], fluxes, index_mask, n_blocks
    )

    bs_v2_cov = _compute_bs_v2_cov(
        bs_arr, v2_arr_unbiased, v2_quantities["v2"], bs_quantities["bs"], index_mask
    )

    if compute_cp_cov:
        cp_cov = _compute_cp_cov(
            bs_arr, bs_quantities["bs"], index_mask, disable=np.invert(verbose)
        )
    else:
        cp_cov = None

    # 9. Now normalize all extracted observables
    vis2_norm, obs_norm = _normalize_all_obs(
        bs_quantities,
        v2_quantities,
        cvis_arr,
        cp_cov,
        bs_v2_cov,
        fluxes,
        index_mask,
        infos,
        expert_plot=display,
    )
    if save_to is not None:
        plt.savefig(f"{figname}_{ifig}.pdf")
    ifig += 1

    obs_result["vis2"] = vis2_norm

    # 10. Now we compute the cp quantities and store them with the other observables
    obs_result = _compute_cp(obs_result, obs_norm, infos, expert_plot=display)

    if save_to is not None:
        plt.savefig(f"{figname}_{ifig}.pdf")
    ifig += 1

    if display:
        _show_norm_matrices(obs_norm, expert_plot=expert_plot)
        if save_to is not None:
            plt.savefig(f"{figname}_{ifig}.pdf")
        ifig += 1

    t3_coord, bl_cp = _compute_t3_coord(mf, index_mask)
    bl_v2 = np.sqrt(mf.u**2 + mf.v**2)
    obs_result["bl"] = bl_v2
    obs_result["bl_cp"] = bl_cp

    # 11. Now we compute the uncertainties using the covariance matrix (for v2)
    # and the variance matrix for the cp.
    obs_result = _compute_uncertainties(obs_result, obs_norm, naive_err=naive_err)

    # 12. Compute scaling error due to phase error (piston) between holes.
    fitmat = _compute_phs_piston(complex_bs, index_mask, display=expert_plot)
    phs_v2corr = _compute_phs_error(complex_bs, fitmat, index_mask, npix)
    obs_norm["phs_v2corr"] = phs_v2corr

    # 13. Compute the absolute oriention (North-up, East-left)
    # ------------------------------------------------------------------------
    pa = compute_pa(hdr, n_ps, display=display, verbose=verbose, sci_hdr=sci_hdr)

    # Compile informations in the storage infos class
    infos = _add_infos_header(infos, hdr, mf, pa, filename, maskname, npix)

    mask = {
        "bl2h_ix": index_mask.bl2h_ix,
        "bs2bl_ix": index_mask.bs2bl_ix,
        "closing_tri": closing_tri,
        "xycoord": mf.xy_coords,
        "n_holes": index_mask.n_holes,
        "n_baselines": index_mask.n_baselines,
        "t3_coord": t3_coord,
    }

    # Finally we store the computed matrices (cov, var, arr, etc,), the informations
    # and the mask parameters to the final output.
    obs_result["mask"] = mask
    obs_result["infos"] = infos
    obs_result["matrix"] = obs_norm

    t = time.time() - start_time
    m = t // 60

    if save_to is not None:
        produce_result_pdf(save_to, Path(filename).stem)

    if verbose:
        rprint(f"[magenta]\nDone (exec time: {m} min {t - m * 60:2.1f} s).")
    return dict2class(obs_result)