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 lsst.pex.config as pexConfig 

28import numpy as np 

29from lsst.pex.config.configurableActions import ConfigurableActionField 

30 

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

32from ...math import isPercent 

33from ..config import MagnitudeBinConfig 

34from .calcBinnedCompleteness import CalcBinnedCompletenessAction 

35 

36 

37class MagnitudeCompletenessConfig(pexConfig.Config): 

38 """Configuration for measuring magnitudes at given completeness 

39 thresholds.""" 

40 

41 completeness_mag_max = pexConfig.Field[float]( 

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

43 default=18, 

44 ) 

45 completeness_percentiles = pexConfig.ListField[float]( 

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

47 default=[90.0, 80.0, 50.0], 

48 itemCheck=isPercent, 

49 ) 

50 

51 

52class CalcCompletenessHistogramAction(KeyedDataAction): 

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

54 

55 action = ConfigurableActionField[CalcBinnedCompletenessAction]( 

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

57 ) 

58 bins = pexConfig.ConfigField[MagnitudeBinConfig]( 

59 doc="The magnitude bin configuration", 

60 ) 

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

62 

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

64 band = kwargs.get("band") 

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

66 bin_width = self.bins.mag_width / 1000.0 

67 n_bins = len(bins) 

68 keys_raw = { 

69 key_formatted: key 

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

71 if not self._is_action_key_mask(key) 

72 } 

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

74 action = copy.copy(self.action) 

75 name_completeness = action.name_completeness 

76 if band is not None: 

77 name_completeness = name_completeness.format(band=band) 

78 

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

80 percentiles = list(percentile_completeness.keys()) 

81 n_percentiles = len(percentiles) 

82 idx_percentile = 0 

83 # isfinite won't work on None 

84 completeness_last, median_last = np.nan, np.nan 

85 

86 for idx_rev, minimum in enumerate(bins): 

87 maximum = minimum + bin_width 

88 median = (maximum + minimum) / 2.0 

89 action.selector_range_ref.minimum = minimum 

90 action.selector_range_ref.maximum = maximum 

91 action.selector_range_target.minimum = minimum 

92 action.selector_range_target.maximum = maximum 

93 result = action(data, **kwargs) 

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

95 value = result[key_formatted] 

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

97 # value is masked, so check for that first 

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

99 if median >= self.config_metrics.completeness_mag_max: 

100 completeness = result[name_completeness] * 100 

101 if completeness: 

102 if idx_percentile > 0: 

103 if completeness < percentiles[idx_percentile - 1]: 

104 idx_percentile -= 1 

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

106 # Crude linear interpolation 

107 # TODO: Replace with a spline or anything better 

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

109 percentile = percentiles[idx_percentile] 

110 width = completeness - completeness_last 

111 # The abs should be unnecessary but just in case 

112 magnitude = ( 

113 abs(completeness - percentile) / width * median_last 

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

115 ) 

116 else: 

117 magnitude = median 

118 percentile_completeness[percentiles[idx_percentile]] = magnitude 

119 idx_percentile += 1 

120 completeness_last, median_last = completeness, median 

121 

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

123 name_percentile = self.getPercentileName(percentile) 

124 key = action.name_mag_completeness(name_percentile) 

125 if band is not None: 

126 key = key.format(band=band) 

127 results[key] = magnitude 

128 

129 return results 

130 

131 def _is_action_key_mask(self, key: str): 

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

133 return is_mask 

134 

135 def getInputSchema(self) -> KeyedDataSchema: 

136 yield from self.action.getInputSchema() 

137 

138 def getOutputSchema(self) -> KeyedDataSchema: 

139 result = { 

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

141 } 

142 return result 

143 

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

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