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-42043: add deltaSkyCorr plot to analysis_tools #174

Merged
merged 1 commit into from
Jan 3, 2024
Merged
Show file tree
Hide file tree
Changes from all 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
8 changes: 8 additions & 0 deletions pipelines/visitQualityExtended.yaml
Original file line number Diff line number Diff line change
Expand Up @@ -5,3 +5,11 @@ tasks:
class: lsst.analysis.tools.tasks.astrometricCatalogMatch.AstrometricCatalogMatchVisitTask
refCatSourceVisit:
class: lsst.analysis.tools.tasks.refCatSourceAnalysis.RefCatSourceAnalysisTask
deltaSkyCorrHist:
class: lsst.analysis.tools.tasks.deltaSkyCorrAnalysis.DeltaSkyCorrHistTask
deltaSkyCorrAnalysis:
class: lsst.analysis.tools.tasks.deltaSkyCorrAnalysis.DeltaSkyCorrAnalysisTask
config:
atools.deltaSkyCorr: DeltaSkyCorrXYPlot
python: |
from lsst.analysis.tools.atools import *
89 changes: 89 additions & 0 deletions python/lsst/analysis/tools/actions/scalar/scalarActions.py
Original file line number Diff line number Diff line change
Expand Up @@ -36,6 +36,8 @@
"FracInRange",
"FracNan",
"SumAction",
"MedianHistAction",
"IqrHistAction",
)

import operator
Expand Down Expand Up @@ -246,3 +248,90 @@ def __call__(self, data: KeyedData, **kwargs) -> Scalar:
mask = self.getMask(**kwargs)
arr = cast(Vector, data[self.vectorKey.format(**kwargs)])[mask]
return cast(Scalar, np.nansum(arr))


class MedianHistAction(ScalarAction):
Copy link
Contributor

Choose a reason for hiding this comment

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

I think you also need to add MedianHistAction and IqrHistAction to the __all__.

Copy link
Contributor

Choose a reason for hiding this comment

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

Done by @aemerywatkins.

"""Calculates the median of the given histogram data."""

histKey = Field[str]("Key of frequency Vector")
midKey = Field[str]("Key of bin midpoints Vector")

def getInputSchema(self) -> KeyedDataSchema:
return (
(self.histKey, Vector),
(self.midKey, Vector),
)

def histMedian(self, hist, bin_mid):
"""Calculates the median of a histogram with binned values

Parameters
----------
hist : `numpy.ndarray`
Frequency array
bin_mid : `numpy.ndarray`
Bin midpoints array

Returns
-------
median : `float`
Median of histogram with binned values
"""
cumulative_sum = np.cumsum(hist)
median_index = np.searchsorted(cumulative_sum, cumulative_sum[-1] / 2)
median = bin_mid[median_index]
return median

def __call__(self, data: KeyedData, **kwargs):
if len(data[self.histKey.format(**kwargs)]) != 0:
hist = cast(Vector, data[self.histKey.format(**kwargs)])
Copy link
Contributor

Choose a reason for hiding this comment

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

You don't really need any casts here, unless you are thinking histMedian should be typed with Vector Scalar etc.

bin_mid = cast(Vector, data[self.midKey.format(**kwargs)])
med = cast(Scalar, float(self.histMedian(hist, bin_mid)))
else:
med = np.NaN
return med


class IqrHistAction(ScalarAction):
Copy link
Contributor

Choose a reason for hiding this comment

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

My preference is for IQRHistAction, but I admit that this style is more in keeping with other aspects of the code base in this file (e.g., sigmaMad). I will leave it up to you as to your preference (also a question for @sr525).

Copy link
Collaborator

Choose a reason for hiding this comment

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

Also elsewhere where we use PsfFlux (for example) rather than PSFFlux. So my vote is for IqrHistAction.

"""Calculates the interquartile range of the given histogram data."""

histKey = Field[str]("Key of frequency Vector")
midKey = Field[str]("Key of bin midpoints Vector")

def getInputSchema(self) -> KeyedDataSchema:
return (
(self.histKey, Vector),
(self.midKey, Vector),
)

def histIqr(self, hist, bin_mid):
"""Calculates the interquartile range of a histogram with binned values

Parameters
----------
hist : `numpy.ndarray`
Frequency array
bin_mid : `numpy.ndarray`
Bin midpoints array

Returns
-------
iqr : `float`
Inter-quartile range of histogram with binned values
"""
cumulative_sum = np.cumsum(hist)
liqr_index = np.searchsorted(cumulative_sum, cumulative_sum[-1] / 4)
uiqr_index = np.searchsorted(cumulative_sum, (3 / 4) * cumulative_sum[-1])
liqr = bin_mid[liqr_index]
uiqr = bin_mid[uiqr_index]
iqr = uiqr - liqr
return iqr

def __call__(self, data: KeyedData, **kwargs):
if len(data[self.histKey.format(**kwargs)]) != 0:
hist = cast(Vector, data[self.histKey.format(**kwargs)])
bin_mid = cast(Vector, data[self.midKey.format(**kwargs)])
iqr = cast(Scalar, float(self.histIqr(hist, bin_mid)))
else:
iqr = np.NaN
return iqr
1 change: 1 addition & 0 deletions python/lsst/analysis/tools/atools/__init__.py
Original file line number Diff line number Diff line change
@@ -1,6 +1,7 @@
from .astrometricRepeatability import *
from .coveragePlots import *
from .deblenderMetric import *
from .deltaSkyCorr import *
from .diaSolarSystemObjectMetrics import *
from .diaSourceMetrics import *
from .diffimMetadataMetrics import *
Expand Down
67 changes: 67 additions & 0 deletions python/lsst/analysis/tools/atools/deltaSkyCorr.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,67 @@
# This file is part of analysis_tools.
#
# 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/>.

__all__ = ("DeltaSkyCorrXYPlot",)

from ..actions.plot.xyPlot import XYPlot
from ..actions.scalar.scalarActions import IqrHistAction, MedianHistAction
from ..actions.vector import ConstantValue, LoadVector
from ..interfaces import AnalysisTool


class DeltaSkyCorrXYPlot(AnalysisTool):
parameterizedBand: bool = False

def setDefaults(self):
super().setDefaults()
self.process.buildActions.x = LoadVector()
self.process.buildActions.y = LoadVector()
self.process.buildActions.x.vectorKey = "bin_mid"
self.process.buildActions.y.vectorKey = "hist"
self.process.buildActions.xerr = ConstantValue(value=0)
self.process.buildActions.yerr = ConstantValue(value=0)

self.process.calculateActions.median = MedianHistAction()
self.process.calculateActions.median.histKey = "hist"
self.process.calculateActions.median.midKey = "bin_mid"
self.process.calculateActions.iqr = IqrHistAction()
self.process.calculateActions.iqr.histKey = "hist"
self.process.calculateActions.iqr.midKey = "bin_mid"

self.produce.plot = XYPlot()
self.produce.plot.xAxisLabel = r"$\Delta$skyCorr Flux (nJy)"
self.produce.plot.yAxisLabel = "Frequency"
self.produce.plot.yScale = "log"
self.produce.plot.xLine = 0
self.produce.plot.strKwargs = {
"fmt": "-",
"color": "b",
}

self.produce.metric.units = {
"median": "nJy",
"iqr": "nJy",
}

self.produce.metric.newNames = {
"median": "delta_skyCorr_median",
"iqr": "delta_skyCorr_iqr",
}