diaCalculationPlugins.py

Go to the documentation of this file.

# This file is part of ap_association.
#
# Developed for the LSST Data Management System.
# This product includes software developed by the LSST Project
# (https://www.lsst.org).
# See the COPYRIGHT file at the top-level directory of this distribution
# for details of code ownership.
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program.  If not, see <https://www.gnu.org/licenses/>.
 
"""Plugins for use in DiaSource summary statistics.
 
Output columns must be
as defined in the schema of the Apdb both in name and units.
"""
 
import functools
import warnings
 
from astropy.stats import median_absolute_deviation
import numpy as np
import pandas as pd
from scipy.optimize import lsq_linear
 
import lsst.geom as geom
import lsst.pex.config as pexConfig
import lsst.pipe.base as pipeBase
import lsst.sphgeom as sphgeom
from astropy.timeseries import LombScargle
from astropy.timeseries import LombScargleMultiband
import math
import statistics
 
from .diaCalculation import (
    DiaObjectCalculationPluginConfig,
    DiaObjectCalculationPlugin)
from .pluginRegistry import register
 
 
__all__ = ("MeanDiaPositionConfig", "MeanDiaPosition",
           "HTMIndexDiaPosition", "HTMIndexDiaPositionConfig",
           "NumDiaSourcesDiaPlugin", "NumDiaSourcesDiaPluginConfig",
           "SimpleSourceFlagDiaPlugin", "SimpleSourceFlagDiaPluginConfig",
           "WeightedMeanDiaPsfFluxConfig", "WeightedMeanDiaPsfFlux",
           "PercentileDiaPsfFlux", "PercentileDiaPsfFluxConfig",
           "SigmaDiaPsfFlux", "SigmaDiaPsfFluxConfig",
           "Chi2DiaPsfFlux", "Chi2DiaPsfFluxConfig",
           "MadDiaPsfFlux", "MadDiaPsfFluxConfig",
           "SkewDiaPsfFlux", "SkewDiaPsfFluxConfig",
           "MinMaxDiaPsfFlux", "MinMaxDiaPsfFluxConfig",
           "MaxSlopeDiaPsfFlux", "MaxSlopeDiaPsfFluxConfig",
           "ErrMeanDiaPsfFlux", "ErrMeanDiaPsfFluxConfig",
           "LinearFitDiaPsfFlux", "LinearFitDiaPsfFluxConfig",
           "StetsonJDiaPsfFlux", "StetsonJDiaPsfFluxConfig",
           "WeightedMeanDiaTotFlux", "WeightedMeanDiaTotFluxConfig",
           "SigmaDiaTotFlux", "SigmaDiaTotFluxConfig",
           "LombScarglePeriodogram", "LombScarglePeriodogramConfig",
           "LombScarglePeriodogramMulti", "LombScarglePeriodogramMultiConfig",
           "UnphysicalDiaSourceSeparation")
 
 
def catchWarnings(_func=None, *, warns=[]):
    """Decorator for generically catching numpy warnings.
    """
    def decoratorCatchWarnings(func):
        @functools.wraps(func)
        def wrapperCatchWarnings(*args, **kwargs):
            with warnings.catch_warnings():
                for val in warns:
                    warnings.filterwarnings("ignore", val)
                return func(*args, **kwargs)
        return wrapperCatchWarnings
 
    if _func is None:
        return decoratorCatchWarnings
    else:
        return decoratorCatchWarnings(_func)
 
 
def typeSafePandasAssignment(
    target,
    source,
    columns,
    default_dtype=np.float64,
    int_fill_value=0,
    # TODO DM-53254: Remove the force_int_to_float hack.
    force_int_to_float=False,
):
    """
    Assign from a source dataframe to a target dataframe in a type safe way.
 
    Parameters
    ----------
    target : `pd.DataFrame`
        Target pandas dataframe.
    source : `pd.DataFrame` or `pd.Series`
        Grouped source dataframe.
    columns : `list` [`str`]
        List of columns to transfer.
    default_dtype : `np.dtype`, optional
        Default datatype (if not in target).
    int_fill_value : `int`, optional
        Fill value for integer columns to avoid pandas insisting
        that everything should be float-ified as nans.
    force_int_to_float : `bool`, optional
        Force integer columns to float columns? Use this option
        for backwards compatibility for old pandas misfeatures which
        are expected by some downstream processes.
    """
    is_series = isinstance(source, pd.Series)
    for col in columns:
        if is_series:
            source_col = source
        else:
            source_col = source[col]
 
        matched_length = False
        if col in target.columns:
            target_dtype = target[col].dtype
            matched_length = len(target) == len(source)
        else:
            target_dtype = default_dtype
 
        if (matched_length or pd.api.types.is_float_dtype(target_dtype)) and not force_int_to_float:
            # If we have a matched length or float, we can do a
            # straight assignment here.
            target.loc[:, col] = source_col.astype(target_dtype)
        else:
            # With mis-matched integers, we must do this dance to preserve types.
            # Note that this may lose precision with very large numbers.
 
            # Convert to float
            target[col] = target[col].astype(np.float64)
            # Set the column, casting to float.
            target.loc[:, col] = source_col.astype(np.float64)
            if not force_int_to_float:
                # Convert back to integer
                target[col] = target[col].fillna(int_fill_value).astype(target_dtype)
 
 
def compute_optimized_periodogram_grid(x0, oversampling_factor=5, nyquist_factor=100):
    """
    Computes an optimized periodogram frequency grid for a given time series.
 
    Parameters
    ----------
    x0 : `array`
        The input time axis.
    oversampling_factor : `int`, optional
        The oversampling factor for frequency grid.
    nyquist_factor : `int`, optional
        The Nyquist factor for frequency grid.
 
    Returns
    -------
    frequencies : `array`
        The computed optimized periodogram frequency grid.
    """
 
    num_points = len(x0)
    baseline = np.max(x0) - np.min(x0)
 
    # Calculate the frequency resolution based on oversampling factor and baseline
    frequency_resolution = 1. / baseline / oversampling_factor
 
    num_frequencies = int(
        0.5 * oversampling_factor * nyquist_factor * num_points)
    frequencies = frequency_resolution + \
        frequency_resolution * np.arange(num_frequencies)
 
    return frequencies
 
 
class LombScarglePeriodogramConfig(DiaObjectCalculationPluginConfig):
    pass
 
 
@register("ap_lombScarglePeriodogram")
class LombScarglePeriodogram(DiaObjectCalculationPlugin):
    """Compute the single-band period of a DiaObject given a set of DiaSources.
    """
    ConfigClass = LombScarglePeriodogramConfig
 
    plugType = "multi"
    outputCols = ["period", "power"]
    needsFilter = True
 
    @classmethod
    def getExecutionOrder(cls):
        return cls.DEFAULT_CATALOGCALCULATION
 
    @catchWarnings(warns=["All-NaN slice encountered"])
    def calculate(self,
                  diaObjects,
                  diaSources,
                  filterDiaSources,
                  band):
        """Compute the periodogram.
 
        Parameters
        ----------
        diaObjects : `pandas.DataFrame`
            Summary objects to store values in.
        diaSources : `pandas.DataFrame` or `pandas.DataFrameGroupBy`
            Catalog of DiaSources summarized by this DiaObject.
        """
 
        # Check and initialize output columns in diaObjects.
        if (periodCol := f"{band}_period") not in diaObjects.columns:
            diaObjects[periodCol] = np.nan
        if (powerCol := f"{band}_power") not in diaObjects.columns:
            diaObjects[powerCol] = np.nan
 
        def _calculate_period(df, min_detections=5, nterms=1, oversampling_factor=5, nyquist_factor=100):
            """Compute the Lomb-Scargle periodogram given a set of DiaSources.
 
            Parameters
            ----------
            df : `pandas.DataFrame`
                The input DataFrame.
            min_detections : `int`, optional
                The minimum number of detections.
            nterms : `int`, optional
                The number of terms in the Lomb-Scargle model.
            oversampling_factor : `int`, optional
                The oversampling factor for frequency grid.
            nyquist_factor : `int`, optional
                The Nyquist factor for frequency grid.
 
            Returns
            -------
            pd_tab : `pandas.Series`
                The output DataFrame with the Lomb-Scargle parameters.
            """
            tmpDf = df[~np.logical_or(np.isnan(df["psfFlux"]),
                                      np.isnan(df["midpointMjdTai"]))]
 
            if len(tmpDf) < min_detections:
                return pd.Series({periodCol: np.nan, powerCol: np.nan})
 
            time = tmpDf["midpointMjdTai"].to_numpy()
            flux = tmpDf["psfFlux"].to_numpy()
            flux_err = tmpDf["psfFluxErr"].to_numpy()
 
            lsp = LombScargle(time, flux, dy=flux_err, nterms=nterms)
            f_grid = compute_optimized_periodogram_grid(
                time, oversampling_factor=oversampling_factor, nyquist_factor=nyquist_factor)
            period = 1/f_grid
            power = lsp.power(f_grid)
 
            return pd.Series({periodCol: period[np.argmax(power)],
                              powerCol: np.max(power)})
 
        diaObjects.loc[:, [periodCol, powerCol]
                       ] = filterDiaSources.apply(_calculate_period)
 
 
class LombScarglePeriodogramMultiConfig(DiaObjectCalculationPluginConfig):
    pass
 
 
@register("ap_lombScarglePeriodogramMulti")
class LombScarglePeriodogramMulti(DiaObjectCalculationPlugin):
    """Compute the multi-band LombScargle periodogram of a DiaObject given a set of DiaSources.
    """
    ConfigClass = LombScarglePeriodogramMultiConfig
 
    plugType = "multi"
    outputCols = ["multiPeriod", "multiPower",
                  "multiFap", "multiAmp", "multiPhase"]
    needsFilter = True
 
    @classmethod
    def getExecutionOrder(cls):
        return cls.DEFAULT_CATALOGCALCULATION
 
    @staticmethod
    def calculate_baluev_fap(time, n, maxPeriod, zmax):
        """Calculate the False-Alarm probability using the Baluev approximation.
 
        Parameters
        ----------
        time : `array`
            The input time axis.
        n : `int`
            The number of detections.
        maxPeriod : `float`
            The maximum period in the grid.
        zmax : `float`
            The maximum power in the grid.
 
        Returns
        -------
        fap_estimate : `float`
            The False-Alarm probability Baluev approximation.
 
        Notes
        ----------
        .. [1] Baluev, R. V. 2008, MNRAS, 385, 1279
        .. [2] Süveges, M., Guy, L.P., Eyer, L., et al. 2015, MNRAS, 450, 2052
        """
        if n <= 2:
            return np.nan
 
        gam_ratio = math.factorial(int((n - 1)/2)) / math.factorial(int((n - 2)/2))
        fu = 1/maxPeriod
        return gam_ratio * np.sqrt(
            4*np.pi*statistics.variance(time)
        ) * fu * (1-zmax)**((n-4)/2) * np.sqrt(zmax)
 
    @staticmethod
    def generate_lsp_params(lsp_model, fbest, bands):
        """Generate the Lomb-Scargle parameters.
        Parameters
        ----------
        lsp_model : `astropy.timeseries.LombScargleMultiband`
            The Lomb-Scargle model.
        fbest : `float`
            The best period.
        bands : `array`
            The bands of the time series.
 
        Returns
        -------
        Amp : `array`
            The amplitude of the time series.
        Ph : `array`
            The phase of the time series.
 
        Notes
        ----------
        .. [1] VanderPlas, J. T., & Ivezic, Z. 2015, ApJ, 812, 18
        """
        best_params = lsp_model.model_parameters(fbest, units=True)
 
        name_params = [f"theta_base_{i}" for i in range(3)]
        name_params += [f"theta_band_{band}_{i}" for band in np.unique(bands) for i in range(3)]
 
        df_params = pd.DataFrame([best_params], columns=name_params)
 
        unique_bands = np.unique(bands)
 
        amplitude_band = [np.sqrt(df_params[f"theta_band_{band}_1"]**2
                                  + df_params[f"theta_band_{band}_2"]**2)
                          for band in unique_bands]
        phase_bands = [np.arctan2(df_params[f"theta_band_{band}_2"],
                                  df_params[f"theta_band_{band}_1"]) for band in unique_bands]
 
        amp = [a[0] for a in amplitude_band]
        ph = [p[0] for p in phase_bands]
 
        return amp, ph
 
    @catchWarnings(warns=["All-NaN slice encountered"])
    def calculate(self,
                  diaObjects,
                  diaSources,
                  **kwargs):
        """Compute the multi-band LombScargle periodogram of a DiaObject given
        a set of DiaSources.
 
        Parameters
        ----------
        diaObjects : `pandas.DataFrame`
            Summary objects to store values in.
        diaSources : `pandas.DataFrame` or `pandas.DataFrameGroupBy`
            Catalog of DiaSources summarized by this DiaObject.
        **kwargs : `dict`
            Unused kwargs that are always passed to a plugin.
        """
 
        bands_arr = diaSources['band'].unique().values
        unique_bands = np.unique(np.concatenate(bands_arr))
        # Check and initialize output columns in diaObjects.
        if (periodCol := "multiPeriod") not in diaObjects.columns:
            diaObjects[periodCol] = np.nan
        if (powerCol := "multiPower") not in diaObjects.columns:
            diaObjects[powerCol] = np.nan
        if (fapCol := "multiFap") not in diaObjects.columns:
            diaObjects[fapCol] = np.nan
        ampCol = "multiAmp"
        phaseCol = "multiPhase"
        for i in range(len(unique_bands)):
            ampCol_band = f"{unique_bands[i]}_{ampCol}"
            if ampCol_band not in diaObjects.columns:
                diaObjects[ampCol_band] = np.nan
            phaseCol_band = f"{unique_bands[i]}_{phaseCol}"
            if phaseCol_band not in diaObjects.columns:
                diaObjects[phaseCol_band] = np.nan
 
        def _calculate_period_multi(df, all_unique_bands,
                                    min_detections=9, oversampling_factor=5, nyquist_factor=100):
            """Calculate the multi-band Lomb-Scargle periodogram.
 
            Parameters
            ----------
            df : `pandas.DataFrame`
                The input DataFrame.
            all_unique_bands : `list` of `str`
                List of all bands present in the diaSource table that is being worked on.
            min_detections : `int`, optional
                The minimum number of detections, including all bands.
            oversampling_factor : `int`, optional
                The oversampling factor for frequency grid.
            nyquist_factor : `int`, optional
                The Nyquist factor for frequency grid.
 
            Returns
            -------
            pd_tab : `pandas.Series`
                The output DataFrame with the Lomb-Scargle parameters.
            """
            tmpDf = df[~np.logical_or(np.isnan(df["psfFlux"]),
                                      np.isnan(df["midpointMjdTai"]))]
 
            if (len(tmpDf)) < min_detections:
                pd_tab_nodet = pd.Series({periodCol: np.nan,
                                          powerCol: np.nan,
                                          fapCol: np.nan})
                for band in all_unique_bands:
                    pd_tab_nodet[f"{band}_{ampCol}"] = np.nan
                    pd_tab_nodet[f"{band}_{phaseCol}"] = np.nan
 
                return pd_tab_nodet
 
            time = tmpDf["midpointMjdTai"].to_numpy()
            flux = tmpDf["psfFlux"].to_numpy()
            flux_err = tmpDf["psfFluxErr"].to_numpy()
            bands = tmpDf["band"].to_numpy()
 
            lsp = LombScargleMultiband(time, flux, bands, dy=flux_err,
                                       nterms_base=1, nterms_band=1)
 
            f_grid = compute_optimized_periodogram_grid(
                time, oversampling_factor=oversampling_factor, nyquist_factor=nyquist_factor)
            period = 1/f_grid
            power = lsp.power(f_grid)
 
            fap_estimate = self.calculate_baluev_fap(
                time, len(time), period[np.argmax(power)], np.max(power))
 
            params_table_new = self.generate_lsp_params(lsp, f_grid[np.argmax(power)], bands)
 
            pd_tab = pd.Series({periodCol: period[np.argmax(power)],
                                powerCol: np.max(power),
                                fapCol: fap_estimate
                                })
 
            # Initialize the per-band amplitude/phase columns as NaNs
            for band in all_unique_bands:
                pd_tab[f"{band}_{ampCol}"] = np.nan
                pd_tab[f"{band}_{phaseCol}"] = np.nan
 
            # Populate the values of only the bands that have data for this diaSource
            unique_bands = np.unique(bands)
            for i in range(len(unique_bands)):
                pd_tab[f"{unique_bands[i]}_{ampCol}"] = params_table_new[0][i]
                pd_tab[f"{unique_bands[i]}_{phaseCol}"] = params_table_new[1][i]
 
            return pd_tab
 
        columns_list = [periodCol, powerCol, fapCol]
        for i in range(len(unique_bands)):
            columns_list.append(f"{unique_bands[i]}_{ampCol}")
            columns_list.append(f"{unique_bands[i]}_{phaseCol}")
 
        diaObjects.loc[:, columns_list
                       ] = diaSources.apply(_calculate_period_multi, unique_bands)
 
 
class UnphysicalDiaSourceSeparation(pipeBase.AlgorithmError):
    """Raised if associated DiaSources are unphysically separated.
 
    Parameters
    ----------
    separation : `float`
        Observed separation in arseconds.
    max_allowed_separation : `float`
        Configured maximum separation in arcseconds.
    """
 
    def __init__(self, separation, max_allowed_separation) -> None:
        self._separation = separation
        self._max_allowed_separation = max_allowed_separation
        super().__init__(f"Observed DiaSource separation {separation} exceeds allowed value of "
                         f"{max_allowed_separation}")
 
    @property
    def metadata(self) -> dict:
        return {
            "separation": self._separation,
            "max_allowed_separation": self._max_allowed_separation,
        }
 
 
class MeanDiaPositionConfig(DiaObjectCalculationPluginConfig):
    MaxAllowedDiaSourceSeparation = pexConfig.Field(
        dtype=float,
        default=3.0,
        doc="Max allowed separation of associated DiaSources in arcsec. "
            "Raises if unphysical separation is found. "
    )
 
 
@register("ap_meanPosition")
class MeanDiaPosition(DiaObjectCalculationPlugin):
    """Compute the mean position and position uncertainty of a DiaObject
    from its associated DiaSources.
 
    The reported (ra, dec) is the inverse-variance weighted spherical
    mean of the per-source positions, using only the DiaSources whose
    per-source coordinate errors are finite.  The reported covariance
    is the covariance of that weighted-mean estimator,
 
        C_formal = (sum_i inv(C_i))^-1,
 
    where ``C_i`` is the per-source 2x2 covariance, multiplied by a
    *chi-squared scale factor* that inflates the result when the
    empirical scatter is larger than the per-source errors predict:
 
        chi2  = sum_i r_i^T inv(C_i) r_i      (r_i = tangent-plane
                                               residual from the mean)
        dof   = 2 * (N - 1)
        C_obj = max(1, chi2 / dof) * C_formal
 
    When the data agree with the per-source errors (chi2 ~ dof) the
    scale factor is unity and ``C_obj == C_formal``.  When the residuals
    exceed what the errors predict, the covariance is inflated.
    """
 
    ConfigClass = MeanDiaPositionConfig
 
    plugType = 'multi'
 
    outputCols = ["ra", "dec", "raErr", "decErr", "ra_dec_Cov"]
    needsFilter = False
 
    @classmethod
    def getExecutionOrder(cls):
        return cls.DEFAULT_CATALOGCALCULATION
 
    def calculate(self, diaObjects, diaSources, **kwargs):
        """Compute the mean ra/dec position and its uncertainty for each
        DiaObject.
 
        Parameters
        ----------
        diaObjects : `pandas.DataFrame`
            Summary objects to store values in.
        diaSources : `pandas.DataFrame` or `pandas.DataFrameGroupBy`
            Catalog of DiaSources summarized by this DiaObject.
        **kwargs
            Any additional keyword arguments that may be passed to the plugin.
        """
        for outCol in self.outputCols:
            if outCol not in diaObjects.columns:
                diaObjects[outCol] = np.nan
 
        maxAllowedSep = self.config.MaxAllowedDiaSourceSeparation
 
        def _column(df, name):
            """Read a column as float64, or return all-NaN if absent."""
            if name in df.columns:
                return df[name].to_numpy(dtype=np.float64)
            return np.full(len(df), np.nan)
 
        def _computeMeanPos(df):
            coords = list(geom.SpherePoint(src["ra"], src["dec"], geom.degrees)
                          for idx, src in df.iterrows())
            # Quick check of the unweighted mean to catch unphysical objects
            # and as a fallback.
            unweightedAvg = geom.averageSpherePoint(coords)
            maxSep = max(unweightedAvg.separation(coord).asArcseconds() for coord in coords)
            if maxSep > maxAllowedSep:
                raise UnphysicalDiaSourceSeparation(maxSep, maxAllowedSep)
 
            nSrc = len(coords)
            raErr = _column(df, "raErr")
            decErr = _column(df, "decErr")
            raDecCov = _column(df, "ra_dec_Cov")
            finiteDiag = (np.isfinite(raErr) & (raErr > 0)) & (np.isfinite(decErr) & (decErr > 0))
            nDiag = int(finiteDiag.sum())
 
            # Tangent-plane offsets from the unweighted mean reference,
            # in degrees (arc-length east, north).
            offsets = np.array(
                [[off.asDegrees() for off in unweightedAvg.getTangentPlaneOffset(coord)] for coord in coords]
            ) if nSrc >= 1 else np.zeros((0, 2))
 
            aveCoord = unweightedAvg
            varRa = np.nan
            varDec = np.nan
            raDecCovObj = np.nan
 
            if nDiag >= 1:
                # Weighted-mean path: include only sources with finite
                # per-source coordinate errors.
                idx = np.where(finiteDiag)[0]
                sigmaRaSq = raErr[idx]**2
                sigmaDecSq = decErr[idx]**2
                # Use full 2x2 weights only when every included source
                # has a finite ra_dec_Cov; otherwise use diagonal weights
                # for all of them and emit NaN for the output covariance.
                haveFullCov = bool(np.isfinite(raDecCov[idx]).all())
 
                C = np.zeros((nDiag, 2, 2), dtype=np.float64)
                C[:, 0, 0] = sigmaRaSq
                C[:, 1, 1] = sigmaDecSq
                if haveFullCov:
                    C[:, 0, 1] = raDecCov[idx]
                    C[:, 1, 0] = raDecCov[idx]
                try:
                    W = np.linalg.inv(C)
                except np.linalg.LinAlgError:
                    # Singular per-source covariance: drop off-diagonal
                    # and retry with diagonal weights only.
                    haveFullCov = False
                    C[:, 0, 1] = 0.0
                    C[:, 1, 0] = 0.0
                    W = np.linalg.inv(C)
 
                W_sum = W.sum(axis=0)
                C_formal = np.linalg.inv(W_sum)
 
                # Weighted-mean tangent-plane offset from the unweighted
                # reference: mu_offset = (sum W_i)^-1 (sum W_i r_i).  Then
                # apply that offset to the unweighted mean to get the
                # weighted spherical mean.  SpherePoint.offset takes a
                # bearing measured counter-clockwise from east (bearing 0
                # is due east, bearing pi/2 is due north).
                mu_offset = C_formal @ np.einsum('nij,nj->i', W, offsets[idx])
                east, north = float(mu_offset[0]), float(mu_offset[1])
                sep_deg = float(np.hypot(east, north))
                if sep_deg > 0.0:
                    bearing = geom.Angle(float(np.arctan2(north, east)), geom.radians)
                    aveCoord = unweightedAvg.offset(bearing, sep_deg*geom.degrees)
 
                varRaFormal = float(C_formal[0, 0])
                varDecFormal = float(C_formal[1, 1])
                covFormal = float(C_formal[0, 1]) if haveFullCov else np.nan
 
                if nDiag == 1:
                    # Single included source: pass its own 2x2 covariance through.
                    varRa = varRaFormal
                    varDec = varDecFormal
                    raDecCovObj = covFormal
                else:
                    # Inflate the variance if the reduced chi**2 fit of the
                    # diaSource coordinates is greater than 1, to account for
                    # scatter.
                    residuals = offsets[idx] - mu_offset
                    if haveFullCov:
                        # Vectorized Σₙ rₙᵀ Wₙ rₙ (r = residuals).
                        chi2 = float(np.einsum('ni,nij,nj->', residuals, W, residuals))
                    else:
                        chi2 = float(np.sum(residuals[:, 0]**2/sigmaRaSq + residuals[:, 1]**2/sigmaDecSq))
                    dof = 2*(nDiag - 1)
                    scale = max(1.0, chi2/dof)
                    varRa = varRaFormal*scale
                    varDec = varDecFormal*scale
                    raDecCovObj = (covFormal*scale if np.isfinite(covFormal) else np.nan)
            else:
                # No source has finite per-source coordinate errors.  Fall
                # back to the unweighted mean position plus the empirical
                # SEM (when N >= 2), and warn that the per-source errors
                # were not usable.
                warnings.warn(
                    "No DiaSources with finite coordinate errors; falling "
                    "back to the unweighted mean position with scatter-only "
                    "uncertainty.",
                    stacklevel=2,
                )
                if nSrc >= 2:
                    sampleCov = np.cov(offsets, rowvar=False, ddof=1)
                    semCov = sampleCov/nSrc
                    varRa = float(semCov[0, 0])
                    varDec = float(semCov[1, 1])
                    raDecCovObj = float(semCov[0, 1])
                # else: nSrc == 1 with no error -- outputs stay NaN.
 
            raErrObj = float(np.sqrt(varRa)) if np.isfinite(varRa) else np.nan
            decErrObj = float(np.sqrt(varDec)) if np.isfinite(varDec) else np.nan
 
            return pd.Series({"ra": aveCoord.getRa().asDegrees(),
                              "dec": aveCoord.getDec().asDegrees(),
                              "raErr": raErrObj,
                              "decErr": decErrObj,
                              "ra_dec_Cov": raDecCovObj})
 
        ans = diaSources.apply(_computeMeanPos)
        typeSafePandasAssignment(diaObjects, ans, ["ra", "dec", "raErr", "decErr", "ra_dec_Cov"])
 
 
class HTMIndexDiaPositionConfig(DiaObjectCalculationPluginConfig):
 
    htmLevel = pexConfig.Field(
        dtype=int,
        doc="Level of the HTM pixelization.",
        default=20,
    )
 
 
@register("ap_HTMIndex")
class HTMIndexDiaPosition(DiaObjectCalculationPlugin):
    """Compute the mean position of a DiaObject given a set of DiaSources.
 
    Notes
    -----
    This plugin was implemented to satisfy requirements of old APDB interface
    which required ``pixelId`` column in DiaObject with HTM20 index. APDB
    interface had migrated to not need that information, but we keep this
    plugin in case it may be useful for something else.
    """
    ConfigClass = HTMIndexDiaPositionConfig
 
    plugType = 'single'
 
    inputCols = ["ra", "dec"]
    outputCols = ["pixelId"]
    needsFilter = False
 
    def __init__(self, config, name, metadata):
        DiaObjectCalculationPlugin.__init__(self, config, name, metadata)
        self.pixelator = sphgeom.HtmPixelization(self.config.htmLevel)
 
    @classmethod
    def getExecutionOrder(cls):
        return cls.FLUX_MOMENTS_CALCULATED
 
    def calculate(self, diaObjects, diaObjectId, **kwargs):
        """Compute the mean position of a DiaObject given a set of DiaSources
 
        Parameters
        ----------
        diaObjects : `pandas.dataFrame`
            Summary objects to store values in and read ra/dec from.
        diaObjectId : `int`
            Id of the diaObject to update.
        **kwargs
            Any additional keyword arguments that may be passed to the plugin.
        """
        sphPoint = geom.SpherePoint(
            diaObjects.at[diaObjectId, "ra"] * geom.degrees,
            diaObjects.at[diaObjectId, "dec"] * geom.degrees)
        diaObjects.at[diaObjectId, "pixelId"] = self.pixelator.index(
            sphPoint.getVector())
 
 
class NumDiaSourcesDiaPluginConfig(DiaObjectCalculationPluginConfig):
    pass
 
 
@register("ap_nDiaSources")
class NumDiaSourcesDiaPlugin(DiaObjectCalculationPlugin):
    """Compute the total number of DiaSources associated with this DiaObject.
    """
 
    ConfigClass = NumDiaSourcesDiaPluginConfig
    outputCols = ["nDiaSources"]
    plugType = "multi"
    needsFilter = False
 
    @classmethod
    def getExecutionOrder(cls):
        return cls.DEFAULT_CATALOGCALCULATION
 
    def calculate(self, diaObjects, diaSources, **kwargs):
        """Compute the total number of DiaSources associated with this DiaObject.
 
        Parameters
        ----------
        diaObject : `dict`
            Summary object to store values in and read ra/dec from.
        **kwargs
            Any additional keyword arguments that may be passed to the plugin.
        """
        typeSafePandasAssignment(diaObjects, diaSources.diaObjectId.count(), ["nDiaSources"])
 
 
class SimpleSourceFlagDiaPluginConfig(DiaObjectCalculationPluginConfig):
    pass
 
 
@register("ap_diaObjectFlag")
class SimpleSourceFlagDiaPlugin(DiaObjectCalculationPlugin):
    """Find if any DiaSource is flagged.
 
    Set the DiaObject flag if any DiaSource is flagged.
    """
 
    ConfigClass = NumDiaSourcesDiaPluginConfig
    outputCols = ["flags"]
    plugType = "multi"
    needsFilter = False
 
    @classmethod
    def getExecutionOrder(cls):
        return cls.DEFAULT_CATALOGCALCULATION
 
    def calculate(self, diaObjects, diaSources, **kwargs):
        """Find if any DiaSource is flagged.
 
        Set the DiaObject flag if any DiaSource is flagged.
 
        Parameters
        ----------
        diaObject : `dict`
            Summary object to store values in and read ra/dec from.
        **kwargs
            Any additional keyword arguments that may be passed to the plugin.
        """
        typeSafePandasAssignment(diaObjects, diaSources.flags.any(), ["flags"], default_dtype=np.uint64)
 
 
class WeightedMeanDiaPsfFluxConfig(DiaObjectCalculationPluginConfig):
    pass
 
 
@register("ap_meanFlux")
class WeightedMeanDiaPsfFlux(DiaObjectCalculationPlugin):
    """Compute the weighted mean and mean error on the point source fluxes
    of the DiaSource measured on the difference image.
 
    Additionally store number of usable data points.
    """
 
    ConfigClass = WeightedMeanDiaPsfFluxConfig
    outputCols = ["psfFluxMean", "psfFluxMeanErr", "psfFluxNdata"]
    plugType = "multi"
    needsFilter = True
 
    @classmethod
    def getExecutionOrder(cls):
        return cls.DEFAULT_CATALOGCALCULATION
 
    @catchWarnings(warns=["invalid value encountered",
                          "divide by zero"])
    def calculate(self,
                  diaObjects,
                  diaSources,
                  filterDiaSources,
                  band,
                  **kwargs):
        """Compute the weighted mean and mean error of the point source flux.
 
        Parameters
        ----------
        diaObject : `dict`
            Summary object to store values in.
        diaSources : `pandas.DataFrame`
            DataFrame representing all diaSources associated with this
            diaObject.
        filterDiaSources : `pandas.DataFrame`
            DataFrame representing diaSources associated with this
            diaObject that are observed in the band pass ``band``.
        band : `str`
            Simple, string name of the filter for the flux being calculated.
        **kwargs
            Any additional keyword arguments that may be passed to the plugin.
        """
        meanName = "{}_psfFluxMean".format(band)
        errName = "{}_psfFluxMeanErr".format(band)
        nDataName = "{}_psfFluxNdata".format(band)
        if meanName not in diaObjects.columns:
            diaObjects[meanName] = np.nan
        if errName not in diaObjects.columns:
            diaObjects[errName] = np.nan
        if nDataName not in diaObjects.columns:
            diaObjects[nDataName] = 0
 
        def _weightedMean(df):
            tmpDf = df[~np.logical_or(np.isnan(df["psfFlux"]),
                                      np.isnan(df["psfFluxErr"]))]
            tot_weight = np.nansum(1 / tmpDf["psfFluxErr"] ** 2)
            fluxMean = np.nansum(tmpDf["psfFlux"]
                                 / tmpDf["psfFluxErr"] ** 2)
            fluxMean /= tot_weight
            if tot_weight > 0:
                fluxMeanErr = np.sqrt(1 / tot_weight)
            else:
                fluxMeanErr = np.nan
            nFluxData = len(tmpDf)
 
            return pd.Series({meanName: fluxMean,
                              errName: fluxMeanErr,
                              nDataName: nFluxData},
                             dtype="object")
        df = filterDiaSources.apply(_weightedMean).astype(diaObjects.dtypes[[meanName, errName, nDataName]])
        # TODO DM-53254: Remove the force_int_to_float hack.
        typeSafePandasAssignment(diaObjects, df, [meanName, errName, nDataName], force_int_to_float=True)
 
 
class PercentileDiaPsfFluxConfig(DiaObjectCalculationPluginConfig):
    percentiles = pexConfig.ListField(
        dtype=int,
        default=[5, 25, 50, 75, 95],
        doc="Percentiles to calculate to compute values for. Should be "
            "integer values."
    )
 
 
@register("ap_percentileFlux")
class PercentileDiaPsfFlux(DiaObjectCalculationPlugin):
    """Compute percentiles of diaSource fluxes.
    """
 
    ConfigClass = PercentileDiaPsfFluxConfig
    # Output columns are created upon instantiation of the class.
    outputCols = []
    plugType = "multi"
    needsFilter = True
 
    def __init__(self, config, name, metadata, **kwargs):
        DiaObjectCalculationPlugin.__init__(self,
                                            config,
                                            name,
                                            metadata,
                                            **kwargs)
        self.outputCols = ["psfFluxPercentile{:02d}".format(percent)
                           for percent in self.config.percentiles]
 
    @classmethod
    def getExecutionOrder(cls):
        return cls.DEFAULT_CATALOGCALCULATION
 
    @catchWarnings(warns=["All-NaN slice encountered"])
    def calculate(self,
                  diaObjects,
                  diaSources,
                  filterDiaSources,
                  band,
                  **kwargs):
        """Compute the percentile fluxes of the point source flux.
 
        Parameters
        ----------
        diaObject : `dict`
            Summary object to store values in.
        diaSources : `pandas.DataFrame`
            DataFrame representing all diaSources associated with this
            diaObject.
        filterDiaSources : `pandas.DataFrame`
            DataFrame representing diaSources associated with this
            diaObject that are observed in the band pass ``band``.
        band : `str`
            Simple, string name of the filter for the flux being calculated.
        **kwargs
            Any additional keyword arguments that may be passed to the plugin.
        """
        pTileNames = []
        for tilePercent in self.config.percentiles:
            pTileName = "{}_psfFluxPercentile{:02d}".format(band,
                                                            tilePercent)
            pTileNames.append(pTileName)
            if pTileName not in diaObjects.columns:
                diaObjects[pTileName] = np.nan
 
        def _fluxPercentiles(df):
            pTiles = np.nanpercentile(df["psfFlux"], self.config.percentiles)
            return pd.Series(
                dict((tileName, pTile)
                     for tileName, pTile in zip(pTileNames, pTiles)))
 
        typeSafePandasAssignment(
            diaObjects,
            filterDiaSources.apply(_fluxPercentiles),
            pTileNames,
        )
 
 
class SigmaDiaPsfFluxConfig(DiaObjectCalculationPluginConfig):
    pass
 
 
@register("ap_sigmaFlux")
class SigmaDiaPsfFlux(DiaObjectCalculationPlugin):
    """Compute scatter of diaSource fluxes.
    """
 
    ConfigClass = SigmaDiaPsfFluxConfig
    # Output columns are created upon instantiation of the class.
    outputCols = ["psfFluxSigma"]
    plugType = "multi"
    needsFilter = True
 
    @classmethod
    def getExecutionOrder(cls):
        return cls.DEFAULT_CATALOGCALCULATION
 
    def calculate(self,
                  diaObjects,
                  diaSources,
                  filterDiaSources,
                  band,
                  **kwargs):
        """Compute the sigma fluxes of the point source flux.
 
        Parameters
        ----------
        diaObject : `dict`
            Summary object to store values in.
        diaSources : `pandas.DataFrame`
            DataFrame representing all diaSources associated with this
            diaObject.
        filterDiaSources : `pandas.DataFrame`
            DataFrame representing diaSources associated with this
            diaObject that are observed in the band pass ``band``.
        band : `str`
            Simple, string name of the filter for the flux being calculated.
        **kwargs
            Any additional keyword arguments that may be passed to the plugin.
        """
        # Set "delta degrees of freedom (ddf)" to 1 to calculate the unbiased
        # estimator of scatter (i.e. 'N - 1' instead of 'N').
        column = "{}_psfFluxSigma".format(band)
 
        typeSafePandasAssignment(
            diaObjects,
            filterDiaSources.psfFlux.std(),
            [column],
        )
 
 
class Chi2DiaPsfFluxConfig(DiaObjectCalculationPluginConfig):
    pass
 
 
@register("ap_chi2Flux")
class Chi2DiaPsfFlux(DiaObjectCalculationPlugin):
    """Compute chi2 of diaSource fluxes.
    """
 
    ConfigClass = Chi2DiaPsfFluxConfig
 
    # Required input Cols
    inputCols = ["psfFluxMean"]
    # Output columns are created upon instantiation of the class.
    outputCols = ["psfFluxChi2"]
    plugType = "multi"
    needsFilter = True
 
    @classmethod
    def getExecutionOrder(cls):
        return cls.FLUX_MOMENTS_CALCULATED
 
    @catchWarnings(warns=["All-NaN slice encountered"])
    def calculate(self,
                  diaObjects,
                  diaSources,
                  filterDiaSources,
                  band,
                  **kwargs):
        """Compute the chi2 of the point source fluxes.
 
        Parameters
        ----------
        diaObject : `dict`
            Summary object to store values in.
        diaSources : `pandas.DataFrame`
            DataFrame representing all diaSources associated with this
            diaObject.
        filterDiaSources : `pandas.DataFrame`
            DataFrame representing diaSources associated with this
            diaObject that are observed in the band pass ``band``.
        band : `str`
            Simple, string name of the filter for the flux being calculated.
        **kwargs
            Any additional keyword arguments that may be passed to the plugin.
        """
        meanName = "{}_psfFluxMean".format(band)
        column = "{}_psfFluxChi2".format(band)
 
        def _chi2(df):
            delta = (df["psfFlux"]
                     - diaObjects.at[df.diaObjectId.iat[0], meanName])
            return np.nansum((delta / df["psfFluxErr"]) ** 2)
 
        typeSafePandasAssignment(
            diaObjects,
            filterDiaSources.apply(_chi2),
            [column],
        )
 
 
class MadDiaPsfFluxConfig(DiaObjectCalculationPluginConfig):
    pass
 
 
@register("ap_madFlux")
class MadDiaPsfFlux(DiaObjectCalculationPlugin):
    """Compute median absolute deviation of diaSource fluxes.
    """
 
    ConfigClass = MadDiaPsfFluxConfig
 
    # Required input Cols
    # Output columns are created upon instantiation of the class.
    outputCols = ["psfFluxMAD"]
    plugType = "multi"
    needsFilter = True
 
    @classmethod
    def getExecutionOrder(cls):
        return cls.DEFAULT_CATALOGCALCULATION
 
    @catchWarnings(warns=["All-NaN slice encountered"])
    def calculate(self,
                  diaObjects,
                  diaSources,
                  filterDiaSources,
                  band,
                  **kwargs):
        """Compute the median absolute deviation of the point source fluxes.
 
        Parameters
        ----------
        diaObject : `dict`
            Summary object to store values in.
        diaSources : `pandas.DataFrame`
            DataFrame representing all diaSources associated with this
            diaObject.
        filterDiaSources : `pandas.DataFrame`
            DataFrame representing diaSources associated with this
            diaObject that are observed in the band pass ``band``.
        band : `str`
            Simple, string name of the filter for the flux being calculated.
        **kwargs
            Any additional keyword arguments that may be passed to the plugin.
        """
        column = "{}_psfFluxMAD".format(band)
 
        typeSafePandasAssignment(
            diaObjects,
            filterDiaSources.psfFlux.apply(median_absolute_deviation, ignore_nan=True),
            [column],
        )
 
 
class SkewDiaPsfFluxConfig(DiaObjectCalculationPluginConfig):
    pass
 
 
@register("ap_skewFlux")
class SkewDiaPsfFlux(DiaObjectCalculationPlugin):
    """Compute the skew of diaSource fluxes.
    """
 
    ConfigClass = SkewDiaPsfFluxConfig
 
    # Required input Cols
    # Output columns are created upon instantiation of the class.
    outputCols = ["psfFluxSkew"]
    plugType = "multi"
    needsFilter = True
 
    @classmethod
    def getExecutionOrder(cls):
        return cls.DEFAULT_CATALOGCALCULATION
 
    def calculate(self,
                  diaObjects,
                  diaSources,
                  filterDiaSources,
                  band,
                  **kwargs):
        """Compute the skew of the point source fluxes.
 
        Parameters
        ----------
        diaObject : `dict`
            Summary object to store values in.
        diaSources : `pandas.DataFrame`
            DataFrame representing all diaSources associated with this
            diaObject.
        filterDiaSources : `pandas.DataFrame`
            DataFrame representing diaSources associated with this
            diaObject that are observed in the band pass ``band``.
        band : `str`
            Simple, string name of the filter for the flux being calculated.
        **kwargs
            Any additional keyword arguments that may be passed to the plugin.
        """
        column = "{}_psfFluxSkew".format(band)
 
        typeSafePandasAssignment(
            diaObjects,
            filterDiaSources.psfFlux.skew(),
            [column],
        )
 
 
class MinMaxDiaPsfFluxConfig(DiaObjectCalculationPluginConfig):
    pass
 
 
@register("ap_minMaxFlux")
class MinMaxDiaPsfFlux(DiaObjectCalculationPlugin):
    """Compute min/max of diaSource fluxes.
    """
 
    ConfigClass = MinMaxDiaPsfFluxConfig
 
    # Required input Cols
    # Output columns are created upon instantiation of the class.
    outputCols = ["psfFluxMin", "psfFluxMax"]
    plugType = "multi"
    needsFilter = True
 
    @classmethod
    def getExecutionOrder(cls):
        return cls.DEFAULT_CATALOGCALCULATION
 
    def calculate(self,
                  diaObjects,
                  diaSources,
                  filterDiaSources,
                  band,
                  **kwargs):
        """Compute min/max of the point source fluxes.
 
        Parameters
        ----------
        diaObject : `dict`
            Summary object to store values in.
        diaSources : `pandas.DataFrame`
            DataFrame representing all diaSources associated with this
            diaObject.
        filterDiaSources : `pandas.DataFrame`
            DataFrame representing diaSources associated with this
            diaObject that are observed in the band pass ``band``.
        band : `str`
            Simple, string name of the filter for the flux being calculated.
        **kwargs
            Any additional keyword arguments that may be passed to the plugin.
        """
        minName = "{}_psfFluxMin".format(band)
        if minName not in diaObjects.columns:
            diaObjects[minName] = np.nan
        maxName = "{}_psfFluxMax".format(band)
        if maxName not in diaObjects.columns:
            diaObjects[maxName] = np.nan
 
        typeSafePandasAssignment(
            diaObjects,
            filterDiaSources.psfFlux.min(),
            [minName],
        )
        typeSafePandasAssignment(
            diaObjects,
            filterDiaSources.psfFlux.max(),
            [maxName],
        )
 
 
class MaxSlopeDiaPsfFluxConfig(DiaObjectCalculationPluginConfig):
    pass
 
 
@register("ap_maxSlopeFlux")
class MaxSlopeDiaPsfFlux(DiaObjectCalculationPlugin):
    """Compute the maximum ratio time ordered deltaFlux / deltaTime.
    """
 
    ConfigClass = MinMaxDiaPsfFluxConfig
 
    # Required input Cols
    # Output columns are created upon instantiation of the class.
    outputCols = ["psfFluxMaxSlope"]
    plugType = "multi"
    needsFilter = True
 
    @classmethod
    def getExecutionOrder(cls):
        return cls.DEFAULT_CATALOGCALCULATION
 
    def calculate(self,
                  diaObjects,
                  diaSources,
                  filterDiaSources,
                  band,
                  **kwargs):
        """Compute the maximum ratio time ordered deltaFlux / deltaTime.
 
        Parameters
        ----------
        diaObject : `dict`
            Summary object to store values in.
        diaSources : `pandas.DataFrame`
            DataFrame representing all diaSources associated with this
            diaObject.
        filterDiaSources : `pandas.DataFrame`
            DataFrame representing diaSources associated with this
            diaObject that are observed in the band pass ``band``.
        band : `str`
            Simple, string name of the filter for the flux being calculated.
        **kwargs
            Any additional keyword arguments that may be passed to the plugin.
        """
 
        def _maxSlope(df):
            tmpDf = df[~np.logical_or(np.isnan(df["psfFlux"]),
                                      np.isnan(df["midpointMjdTai"]))]
            if len(tmpDf) < 2:
                return np.nan
            times = tmpDf["midpointMjdTai"].to_numpy()
            timeArgs = times.argsort()
            times = times[timeArgs]
            fluxes = tmpDf["psfFlux"].to_numpy()[timeArgs]
            return (np.diff(fluxes) / np.diff(times)).max()
 
        column = "{}_psfFluxMaxSlope".format(band)
 
        typeSafePandasAssignment(
            diaObjects,
            filterDiaSources.apply(_maxSlope),
            [column],
        )
 
 
class ErrMeanDiaPsfFluxConfig(DiaObjectCalculationPluginConfig):
    pass
 
 
@register("ap_meanErrFlux")
class ErrMeanDiaPsfFlux(DiaObjectCalculationPlugin):
    """Compute the mean of the dia source errors.
    """
 
    ConfigClass = ErrMeanDiaPsfFluxConfig
 
    # Required input Cols
    # Output columns are created upon instantiation of the class.
    outputCols = ["psfFluxErrMean"]
    plugType = "multi"
    needsFilter = True
 
    @classmethod
    def getExecutionOrder(cls):
        return cls.DEFAULT_CATALOGCALCULATION
 
    def calculate(self,
                  diaObjects,
                  diaSources,
                  filterDiaSources,
                  band,
                  **kwargs):
        """Compute the mean of the dia source errors.
 
        Parameters
        ----------
        diaObject : `dict`
            Summary object to store values in.
        diaSources : `pandas.DataFrame`
            DataFrame representing all diaSources associated with this
            diaObject.
        filterDiaSources : `pandas.DataFrame`
            DataFrame representing diaSources associated with this
            diaObject that are observed in the band pass ``band``.
        band : `str`
            Simple, string name of the filter for the flux being calculated.
        **kwargs
            Any additional keyword arguments that may be passed to the plugin.
        """
        column = "{}_psfFluxErrMean".format(band)
 
        typeSafePandasAssignment(
            diaObjects,
            filterDiaSources.psfFluxErr.mean(),
            [column],
            # Note that the schemas expect this to be single-precision.
            default_dtype=np.float32,
        )
 
 
class LinearFitDiaPsfFluxConfig(DiaObjectCalculationPluginConfig):
    pass
 
 
@register("ap_linearFit")
class LinearFitDiaPsfFlux(DiaObjectCalculationPlugin):
    """Compute fit a linear model to flux vs time.
    """
 
    ConfigClass = LinearFitDiaPsfFluxConfig
 
    # Required input Cols
    # Output columns are created upon instantiation of the class.
    outputCols = ["psfFluxLinearSlope", "psfFluxLinearIntercept"]
    plugType = "multi"
    needsFilter = True
 
    @classmethod
    def getExecutionOrder(cls):
        return cls.DEFAULT_CATALOGCALCULATION
 
    def calculate(self,
                  diaObjects,
                  diaSources,
                  filterDiaSources,
                  band,
                  **kwargs):
        """Compute fit a linear model to flux vs time.
 
        Parameters
        ----------
        diaObject : `dict`
            Summary object to store values in.
        diaSources : `pandas.DataFrame`
            DataFrame representing all diaSources associated with this
            diaObject.
        filterDiaSources : `pandas.DataFrame`
            DataFrame representing diaSources associated with this
            diaObject that are observed in the band pass ``band``.
        band : `str`
            Simple, string name of the filter for the flux being calculated.
        **kwargs
            Any additional keyword arguments that may be passed to the plugin.
        """
 
        mName = "{}_psfFluxLinearSlope".format(band)
        if mName not in diaObjects.columns:
            diaObjects[mName] = np.nan
        bName = "{}_psfFluxLinearIntercept".format(band)
        if bName not in diaObjects.columns:
            diaObjects[bName] = np.nan
        dtype = diaObjects[mName].dtype
 
        def _linearFit(df):
            tmpDf = df[~np.logical_or(
                np.isnan(df["psfFlux"]),
                np.logical_or(np.isnan(df["psfFluxErr"]),
                              np.isnan(df["midpointMjdTai"])))]
            if len(tmpDf) < 2:
                return pd.Series({mName: np.nan, bName: np.nan})
            fluxes = tmpDf["psfFlux"].to_numpy()
            errors = tmpDf["psfFluxErr"].to_numpy()
            times = tmpDf["midpointMjdTai"].to_numpy()
            A = np.array([times / errors, 1 / errors]).transpose()
            m, b = lsq_linear(A, fluxes / errors).x
            return pd.Series({mName: m, bName: b}, dtype=dtype)
 
        typeSafePandasAssignment(
            diaObjects,
            filterDiaSources.apply(_linearFit),
            [mName, bName],
        )
 
 
class StetsonJDiaPsfFluxConfig(DiaObjectCalculationPluginConfig):
    pass
 
 
@register("ap_stetsonJ")
class StetsonJDiaPsfFlux(DiaObjectCalculationPlugin):
    """Compute the StetsonJ statistic on the DIA point source fluxes.
    """
 
    ConfigClass = LinearFitDiaPsfFluxConfig
 
    # Required input Cols
    inputCols = ["psfFluxMean"]
    # Output columns are created upon instantiation of the class.
    outputCols = ["psfFluxStetsonJ"]
    plugType = "multi"
    needsFilter = True
 
    @classmethod
    def getExecutionOrder(cls):
        return cls.FLUX_MOMENTS_CALCULATED
 
    def calculate(self,
                  diaObjects,
                  diaSources,
                  filterDiaSources,
                  band,
                  **kwargs):
        """Compute the StetsonJ statistic on the DIA point source fluxes.
 
        Parameters
        ----------
        diaObject : `dict`
            Summary object to store values in.
        diaSources : `pandas.DataFrame`
            DataFrame representing all diaSources associated with this
            diaObject.
        filterDiaSources : `pandas.DataFrame`
            DataFrame representing diaSources associated with this
            diaObject that are observed in the band pass ``band``.
        band : `str`
            Simple, string name of the filter for the flux being calculated.
        **kwargs
            Any additional keyword arguments that may be passed to the plugin.
        """
        meanName = "{}_psfFluxMean".format(band)
 
        def _stetsonJ(df):
            tmpDf = df[~np.logical_or(np.isnan(df["psfFlux"]),
                                      np.isnan(df["psfFluxErr"]))]
            if len(tmpDf) < 2:
                return np.nan
            fluxes = tmpDf["psfFlux"].to_numpy()
            errors = tmpDf["psfFluxErr"].to_numpy()
 
            return self._stetson_J(
                fluxes,
                errors,
                diaObjects.at[tmpDf.diaObjectId.iat[0], meanName])
 
        column = "{}_psfFluxStetsonJ".format(band)
        typeSafePandasAssignment(
            diaObjects,
            filterDiaSources.apply(_stetsonJ),
            [column],
        )
 
    def _stetson_J(self, fluxes, errors, mean=None):
        """Compute the single band stetsonJ statistic.
 
        Parameters
        ----------
        fluxes : `numpy.ndarray` (N,)
            Calibrated lightcurve flux values.
        errors : `numpy.ndarray` (N,)
            Errors on the calibrated lightcurve fluxes.
        mean : `float`
            Starting mean from previous plugin.
 
        Returns
        -------
        stetsonJ : `float`
            stetsonJ statistic for the input fluxes and errors.
 
        References
        ----------
        .. [1] Stetson, P. B., "On the Automatic Determination of Light-Curve
           Parameters for Cepheid Variables", PASP, 108, 851S, 1996
        """
        n_points = len(fluxes)
        flux_mean = self._stetson_mean(fluxes, errors, mean)
        delta_val = (
            np.sqrt(n_points / (n_points - 1)) * (fluxes - flux_mean) / errors)
        p_k = delta_val ** 2 - 1
 
        return np.mean(np.sign(p_k) * np.sqrt(np.fabs(p_k)))
 
    def _stetson_mean(self,
                      values,
                      errors,
                      mean=None,
                      alpha=2.,
                      beta=2.,
                      n_iter=20,
                      tol=1e-6):
        """Compute the stetson mean of the fluxes which down-weights outliers.
 
        Weighted biased on an error weighted difference scaled by a constant
        (1/``a``) and raised to the power beta. Higher betas more harshly
        penalize outliers and ``a`` sets the number of sigma where a weighted
        difference of 1 occurs.
 
        Parameters
        ----------
        values : `numpy.dnarray`, (N,)
            Input values to compute the mean of.
        errors : `numpy.ndarray`, (N,)
            Errors on the input values.
        mean : `float`
            Starting mean value or None.
        alpha : `float`
            Scalar down-weighting of the fractional difference. lower->more
            clipping. (Default value is 2.)
        beta : `float`
            Power law slope of the used to down-weight outliers. higher->more
            clipping. (Default value is 2.)
        n_iter : `int`
            Number of iterations of clipping.
        tol : `float`
            Fractional and absolute tolerance goal on the change in the mean
            before exiting early. (Default value is 1e-6)
 
        Returns
        -------
        mean : `float`
            Weighted stetson mean result.
 
        References
        ----------
        .. [1] Stetson, P. B., "On the Automatic Determination of Light-Curve
           Parameters for Cepheid Variables", PASP, 108, 851S, 1996
        """
        n_points = len(values)
        n_factor = np.sqrt(n_points / (n_points - 1))
        inv_var = 1 / errors ** 2
 
        if mean is None:
            mean = np.average(values, weights=inv_var)
        for iter_idx in range(n_iter):
            chi = np.fabs(n_factor * (values - mean) / errors)
            tmp_mean = np.average(
                values,
                weights=inv_var / (1 + (chi / alpha) ** beta))
            diff = np.fabs(tmp_mean - mean)
            mean = tmp_mean
            if diff / mean < tol and diff < tol:
                break
        return mean
 
 
class WeightedMeanDiaTotFluxConfig(DiaObjectCalculationPluginConfig):
    pass
 
 
@register("ap_meanTotFlux")
class WeightedMeanDiaTotFlux(DiaObjectCalculationPlugin):
    """Compute the weighted mean and mean error on the point source fluxes
    forced photometered at the DiaSource location in the calibrated image.
 
    Additionally store number of usable data points.
    """
 
    ConfigClass = WeightedMeanDiaPsfFluxConfig
    outputCols = ["scienceFluxMean", "scienceFluxMeanErr"]
    plugType = "multi"
    needsFilter = True
 
    @classmethod
    def getExecutionOrder(cls):
        return cls.DEFAULT_CATALOGCALCULATION
 
    @catchWarnings(warns=["invalid value encountered",
                          "divide by zero"])
    def calculate(self,
                  diaObjects,
                  diaSources,
                  filterDiaSources,
                  band,
                  **kwargs):
        """Compute the weighted mean and mean error of the point source flux.
 
        Parameters
        ----------
        diaObject : `dict`
            Summary object to store values in.
        diaSources : `pandas.DataFrame`
            DataFrame representing all diaSources associated with this
            diaObject.
        filterDiaSources : `pandas.DataFrame`
            DataFrame representing diaSources associated with this
            diaObject that are observed in the band pass ``band``.
        band : `str`
            Simple, string name of the filter for the flux being calculated.
        **kwargs
            Any additional keyword arguments that may be passed to the plugin.
        """
        totMeanName = "{}_scienceFluxMean".format(band)
        if totMeanName not in diaObjects.columns:
            diaObjects[totMeanName] = np.nan
        totErrName = "{}_scienceFluxMeanErr".format(band)
        if totErrName not in diaObjects.columns:
            diaObjects[totErrName] = np.nan
 
        def _meanFlux(df):
            tmpDf = df[~np.logical_or(np.isnan(df["scienceFlux"]),
                                      np.isnan(df["scienceFluxErr"]))]
            tot_weight = np.nansum(1 / tmpDf["scienceFluxErr"] ** 2)
            fluxMean = np.nansum(tmpDf["scienceFlux"]
                                 / tmpDf["scienceFluxErr"] ** 2)
            fluxMean /= tot_weight
            fluxMeanErr = np.sqrt(1 / tot_weight)
 
            return pd.Series({totMeanName: fluxMean,
                              totErrName: fluxMeanErr})
 
        df = filterDiaSources.apply(_meanFlux).astype(diaObjects.dtypes[[totMeanName, totErrName]])
        typeSafePandasAssignment(diaObjects, df, [totMeanName, totErrName])
 
 
class SigmaDiaTotFluxConfig(DiaObjectCalculationPluginConfig):
    pass
 
 
@register("ap_sigmaTotFlux")
class SigmaDiaTotFlux(DiaObjectCalculationPlugin):
    """Compute scatter of diaSource fluxes.
    """
 
    ConfigClass = SigmaDiaPsfFluxConfig
    # Output columns are created upon instantiation of the class.
    outputCols = ["scienceFluxSigma"]
    plugType = "multi"
    needsFilter = True
 
    @classmethod
    def getExecutionOrder(cls):
        return cls.DEFAULT_CATALOGCALCULATION
 
    def calculate(self,
                  diaObjects,
                  diaSources,
                  filterDiaSources,
                  band,
                  **kwargs):
        """Compute the sigma fluxes of the point source flux measured on the
        calibrated image.
 
        Parameters
        ----------
        diaObject : `dict`
            Summary object to store values in.
        diaSources : `pandas.DataFrame`
            DataFrame representing all diaSources associated with this
            diaObject.
        filterDiaSources : `pandas.DataFrame`
            DataFrame representing diaSources associated with this
            diaObject that are observed in the band pass ``band``.
        band : `str`
            Simple, string name of the filter for the flux being calculated.
        **kwargs
            Any additional keyword arguments that may be passed to the plugin.
        """
        # Set "delta degrees of freedom (ddf)" to 1 to calculate the unbiased
        # estimator of scatter (i.e. 'N - 1' instead of 'N').
        column = "{}_scienceFluxSigma".format(band)
        typeSafePandasAssignment(diaObjects, filterDiaSources.scienceFlux.std(), [column])