Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

DM-22299: Speed up specific diaCalculation plugins using fast pandas functionality #65

Merged
merged 3 commits into from
Dec 3, 2019
Merged
Show file tree
Hide file tree
Changes from 2 commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
43 changes: 10 additions & 33 deletions python/lsst/ap/association/diaCalculationPlugins.py
Original file line number Diff line number Diff line change
Expand Up @@ -29,7 +29,6 @@
import numpy as np
import pandas as pd
from scipy.optimize import lsq_linear
from scipy.stats import skew

import lsst.geom as geom
from lsst.meas.algorithms.indexerRegistry import IndexerRegistry
Expand Down Expand Up @@ -184,7 +183,7 @@ def calculate(self, diaObjects, diaSources, **kwargs):
diaObject : `dict`
Summary object to store values in and read ra/decl from.
"""
diaObjects.loc[:, "nDiaSources"] = diaSources.apply(len)
diaObjects.loc[:, "nDiaSources"] = diaSources.diaObjectId.count()


class SimpleSourceFlagDiaPluginConfig(DiaObjectCalculationPluginConfig):
Expand Down Expand Up @@ -216,13 +215,7 @@ def calculate(self, diaObjects, diaSources, **kwargs):
diaObject : `dict`
Summary object to store values in and read ra/decl from.
"""

def _flagDiaObject(df):
if np.any(df["flags"] > 0):
return 1
return 0

diaObjects.loc[:, "flags"] = diaSources.apply(_flagDiaObject)
diaObjects.loc[:, "flags"] = diaSources.flags.any()


class WeightedMeanDiaPsFluxConfig(DiaObjectCalculationPluginConfig):
Expand Down Expand Up @@ -404,12 +397,8 @@ def calculate(self,
"""
# Set "delta degrees of freedom (ddf)" to 1 to calculate the unbiased
# estimator of scatter (i.e. 'N - 1' instead of 'N').

def _sigma(df):
return np.nanstd(df["psFlux"], ddof=1)

diaObjects.loc[:, "{}PSFluxSigma".format(filterName)] = \
filterDiaSources.apply(_sigma)
filterDiaSources.psFlux.agg(np.nanstd)


class Chi2DiaPsFluxConfig(DiaObjectCalculationPluginConfig):
Expand Down Expand Up @@ -506,12 +495,9 @@ def calculate(self,
filterName : `str`
Simple, string name of the filter for the flux being calculated.
"""

def _mad(df):
return median_absolute_deviation(df["psFlux"], ignore_nan=True)

diaObjects.loc[:, "{}PSFluxMAD".format(filterName)] = \
filterDiaSources.apply(_mad)
filterDiaSources.psFlux.apply(median_absolute_deviation,
ignore_nan=True)


class SkewDiaPsFluxConfig(DiaObjectCalculationPluginConfig):
Expand Down Expand Up @@ -556,7 +542,7 @@ def calculate(self,
Simple, string name of the filter for the flux being calculated.
"""
diaObjects.loc[:, "{}PSFluxSkew".format(filterName)] = \
filterDiaSources.psFlux.apply(skew, nan_policy='omit')
filterDiaSources.psFlux.skew()


class MinMaxDiaPsFluxConfig(DiaObjectCalculationPluginConfig):
Expand Down Expand Up @@ -607,11 +593,8 @@ def calculate(self,
if maxName not in diaObjects.columns:
diaObjects[maxName] = np.nan

def _minMax(df):
return pd.Series({minName: df["psFlux"].min(),
maxName: df["psFlux"].max()})

diaObjects.loc[:, [minName, maxName]] = filterDiaSources.apply(_minMax)
diaObjects.loc[:, minName] = filterDiaSources.psFlux.min()
diaObjects.loc[:, maxName] = filterDiaSources.psFlux.max()


class MaxSlopeDiaPsFluxConfig(DiaObjectCalculationPluginConfig):
Expand Down Expand Up @@ -712,11 +695,8 @@ def calculate(self,
filterName : `str`
Simple, string name of the filter for the flux being calculated.
"""
def _meanErr(df):
return np.nanmean(df["psFluxErr"])

diaObjects.loc[:, "{}PSFluxErrMean".format(filterName)] = \
filterDiaSources.apply(_meanErr)
filterDiaSources.psFluxErr.agg(np.nanmean)


class LinearFitDiaPsFluxConfig(DiaObjectCalculationPluginConfig):
Expand Down Expand Up @@ -1044,8 +1024,5 @@ def calculate(self,
"""
# Set "delta degrees of freedom (ddf)" to 1 to calculate the unbiased
# estimator of scatter (i.e. 'N - 1' instead of 'N').
def _sigma(df):
return np.nanstd(df["totFlux"], ddof=1)

diaObjects.loc[:, "{}TOTFluxSigma".format(filterName)] = \
filterDiaSources.apply(_sigma)
filterDiaSources.totFlux.agg(np.nanstd)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

check to see if std is nan safe and statistically unbiased / pandas function

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Checked and the pandas functions are nan-safe and unbiased. I have thus changed to using these functions directly rather than agg. Checked for any performance difference using timeit and found the run time to be nearly identical with slightly less variance per loop for the pandas builtin.

19 changes: 12 additions & 7 deletions tests/test_dia_calculation_plugins.py
Original file line number Diff line number Diff line change
Expand Up @@ -543,7 +543,9 @@ def testCalculate(self):
"ap_skewFlux",
None)
run_multi_plugin(diaObjects, diaSources, "u", plug)
self.assertAlmostEqual(diaObjects.loc[objId, "uPSFluxSkew"], skew(fluxes))
self.assertAlmostEqual(
diaObjects.loc[objId, "uPSFluxSkew"],
skew(fluxes, bias=False, nan_policy="omit"))

# Test expected skew value with a nan set.
fluxes[4] = np.nan
Expand All @@ -555,8 +557,11 @@ def testCalculate(self):
"psFlux": fluxes,
"psFluxErr": np.ones(n_sources)})
run_multi_plugin(diaObjects, diaSources, "r", plug)
self.assertAlmostEqual(diaObjects.at[objId, "rPSFluxSkew"],
skew(fluxes, nan_policy="omit"))
# Skew returns a named tuple when called on an array
# with nan values.
self.assertAlmostEqual(
diaObjects.at[objId, "rPSFluxSkew"],
skew(fluxes, bias=False, nan_policy="omit").data)


class TestMinMaxDiaPsFlux(unittest.TestCase):
Expand Down Expand Up @@ -674,8 +679,8 @@ def testCalculate(self):
"ap_errMeanFlux",
None)
run_multi_plugin(diaObjects, diaSources, "u", plug)
self.assertEqual(diaObjects.at[objId, "uPSFluxErrMean"],
np.nanmean(errors))
self.assertAlmostEqual(diaObjects.at[objId, "uPSFluxErrMean"],
np.nanmean(errors))

# Test mean of the errors with input nan value.
errors[4] = np.nan
Expand All @@ -687,8 +692,8 @@ def testCalculate(self):
"psFlux": fluxes,
"psFluxErr": errors})
run_multi_plugin(diaObjects, diaSources, "r", plug)
self.assertEqual(diaObjects.at[objId, "rPSFluxErrMean"],
np.nanmean(errors))
self.assertAlmostEqual(diaObjects.at[objId, "rPSFluxErrMean"],
np.nanmean(errors))


class TestLinearFitDiaPsFlux(unittest.TestCase):
Expand Down