Coverage for python / lsst / analysis / tools / actions / keyedData / calcCompletenessHistogram.py: 22%

1# This file is part of analysis_tools. 

2# 

3# Developed for the LSST Data Management System. 

4# This product includes software developed by the LSST Project 

5# (https://www.lsst.org). 

6# See the COPYRIGHT file at the top-level directory of this distribution 

7# for details of code ownership. 

8# 

9# This program is free software: you can redistribute it and/or modify 

10# it under the terms of the GNU General Public License as published by 

11# the Free Software Foundation, either version 3 of the License, or 

12# (at your option) any later version. 

13# 

14# This program is distributed in the hope that it will be useful, 

15# but WITHOUT ANY WARRANTY; without even the implied warranty of 

16# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 

17# GNU General Public License for more details. 

18# 

19# You should have received a copy of the GNU General Public License 

20# along with this program. If not, see <https://www.gnu.org/licenses/>. 

21from __future__ import annotations 

22 

23__all__ = ("CalcCompletenessHistogramAction", "MagnitudeCompletenessConfig") 

24 

25import copy 

26 

27import numpy as np 

28 

29import lsst.pex.config as pexConfig 

30from lsst.pex.config.configurableActions import ConfigurableActionField 

31 

32from ...interfaces import KeyedData, KeyedDataAction, KeyedDataSchema 

33from ...math import isPercent 

34from ..config import MagnitudeBinConfig 

35from .calcBinnedCompleteness import CalcBinnedCompletenessAction 

36 

37 

38class MagnitudeCompletenessConfig(pexConfig.Config): 

39 """Configuration for measuring magnitudes at given completeness 

40 thresholds.""" 

41 

42 completeness_mag_max = pexConfig.Field[float]( 

43 doc="Brightest magnitude to consider checking if completeness is below a percentile threshold for", 

44 default=18, 

45 ) 

46 completeness_percentiles = pexConfig.ListField[float]( 

47 doc="The percentiles to find the magnitude at.", 

48 default=[90.0, 80.0, 50.0], 

49 itemCheck=isPercent, 

50 ) 

51 

52 

53class CalcCompletenessHistogramAction(KeyedDataAction): 

54 """Action to calculate a histogram of completeness vs magnitude.""" 

55 

56 action = ConfigurableActionField[CalcBinnedCompletenessAction]( 

57 doc="The action to compute completeness/purity", 

58 ) 

59 bins = pexConfig.ConfigField[MagnitudeBinConfig]( 

60 doc="The magnitude bin configuration", 

61 ) 

62 config_metrics = pexConfig.ConfigField[MagnitudeCompletenessConfig](doc="Metric definition configuration") 

63 

64 def __call__(self, data: KeyedData, **kwargs) -> KeyedData: 

65 band = kwargs.get("band") 

66 bins = tuple(x / 1000.0 for x in reversed(self.bins.get_bins())) 

67 bin_width = self.bins.mag_width / 1000.0 

68 n_bins = len(bins) 

69 keys_raw = { 

70 key_formatted: key 

71 for key, key_formatted in self.action.getFormattedOutputKeys(**kwargs).items() 

72 if not self._is_action_key_mask(key) 

73 } 

74 results = {key_formatted: np.zeros(n_bins, dtype=float) for key_formatted in keys_raw.keys()} 

75 action = copy.copy(self.action) 

76 name_completeness = action.name_completeness 

77 if band is not None: 

78 name_completeness = name_completeness.format(band=band) 

79 

80 percentile_completeness = {pc: np.nan for pc in sorted(self.config_metrics.completeness_percentiles)} 

81 percentiles = list(percentile_completeness.keys()) 

82 n_percentiles = len(percentiles) 

83 idx_percentile = 0 

84 # isfinite won't work on None 

85 completeness_last, median_last = np.nan, np.nan 

86 

87 for idx_rev, minimum in enumerate(bins): 

88 maximum = minimum + bin_width 

89 median = (maximum + minimum) / 2.0 

90 action.selector_range_ref.minimum = minimum 

91 action.selector_range_ref.maximum = maximum 

92 action.selector_range_target.minimum = minimum 

93 action.selector_range_target.maximum = maximum 

94 result = action(data, **kwargs) 

95 for key_formatted, array in results.items(): 

96 value = result[key_formatted] 

97 # The implicit float conversion will generate a warning if the 

98 # value is masked, so check for that first 

99 array[n_bins - idx_rev - 1] = np.nan if np.ma.is_masked(value) else value 

100 if median >= self.config_metrics.completeness_mag_max: 

101 completeness = result[name_completeness] * 100 

102 if completeness: 

103 if idx_percentile > 0: 

104 if completeness < percentiles[idx_percentile - 1]: 

105 idx_percentile -= 1 

106 while (idx_percentile < n_percentiles) and (completeness > percentiles[idx_percentile]): 

107 # Crude linear interpolation 

108 # TODO: Replace with a spline or anything better 

109 if np.isfinite(completeness_last) and np.isfinite(median_last): 

110 percentile = percentiles[idx_percentile] 

111 width = completeness - completeness_last 

112 # The abs should be unnecessary but just in case 

113 magnitude = ( 

114 abs(completeness - percentile) / width * median_last 

115 + abs(percentile - completeness_last) / width * median 

116 ) 

117 else: 

118 magnitude = median 

119 percentile_completeness[percentiles[idx_percentile]] = magnitude 

120 idx_percentile += 1 

121 completeness_last, median_last = completeness, median 

122 

123 for percentile, magnitude in percentile_completeness.items(): 

124 name_percentile = self.getPercentileName(percentile) 

125 key = action.name_mag_completeness(name_percentile) 

126 if band is not None: 

127 key = key.format(band=band) 

128 results[key] = magnitude 

129 

130 return results 

131 

132 def _is_action_key_mask(self, key: str): 

133 is_mask = key.startswith(f"{self.action.name_prefix}mask") 

134 return is_mask 

135 

136 def getInputSchema(self) -> KeyedDataSchema: 

137 yield from self.action.getInputSchema() 

138 

139 def getOutputSchema(self) -> KeyedDataSchema: 

140 result = { 

141 (key, typ) for key, typ in self.action.getOutputSchema() if not self._is_action_key_mask(key) 

142 } 

143 return result 

144 

145 def getPercentileName(self, percentile: float) -> str: 

146 return f"{percentile:.2f}_pct".replace(".", "p")