Coverage for python/lsst/analysis/tools/actions/scalar/scalarActions.py: 51%

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/>. 

21 

22from __future__ import annotations 

23 

24__all__ = ( 

25 "MedianAction", 

26 "MeanAction", 

27 "StdevAction", 

28 "ValueAction", 

29 "SigmaMadAction", 

30 "CountAction", 

31 "CountUniqueAction", 

32 "ApproxFloor", 

33 "FracThreshold", 

34 "MaxAction", 

35 "MinAction", 

36 "FracInRange", 

37 "FracNan", 

38 "SumAction", 

39 "WeightedMeanAction", 

40 "MedianHistAction", 

41 "IqrHistAction", 

42 "DivideScalar", 

43 "RmsAction", 

44 "MedianGradientAction", 

45) 

46 

47import logging 

48import operator 

49from math import nan 

50from typing import cast 

51 

52import numpy as np 

53from scipy.optimize import curve_fit 

54from scipy.stats import binned_statistic 

55 

56from lsst.pex.config import ChoiceField, Field 

57from lsst.pex.config.configurableActions import ConfigurableActionField 

58 

59from ...interfaces import KeyedData, KeyedDataSchema, Scalar, ScalarAction, Vector 

60from ...math import nanMax, nanMean, nanMedian, nanMin, nanSigmaMad, nanStd 

61 

62log = logging.getLogger(__name__) 

63 

64 

65def _dataToArray(data): 

66 """Convert input data into a numpy array using the appropriate 

67 protocol. `np.from_dlpack` is used for Tensor-like arrays 

68 where possible. 

69 """ 

70 try: 

71 return np.from_dlpack(data) 

72 except (AttributeError, BufferError): 

73 return np.array(data) 

74 

75 

76class ScalarFromVectorAction(ScalarAction): 

77 """Calculates a statistic from a single vector.""" 

78 

79 vectorKey = Field[str]("Key of Vector to compute statistic from.") 

80 

81 def getInputSchema(self) -> KeyedDataSchema: 

82 return ((self.vectorKey, Vector),) 

83 

84 

85class MedianAction(ScalarFromVectorAction): 

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

87 

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

89 mask = self.getMask(**kwargs) 

90 values = _dataToArray(data[self.vectorKey.format(**kwargs)])[mask] 

91 med = nanMedian(values) if values.size else np.nan 

92 

93 return med 

94 

95 

96class MeanAction(ScalarFromVectorAction): 

97 """Calculates the mean of the given data.""" 

98 

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

100 mask = self.getMask(**kwargs) 

101 values = _dataToArray(data[self.vectorKey.format(**kwargs)])[mask] 

102 mean = nanMean(values) if values.size else np.nan 

103 

104 return mean 

105 

106 

107class WeightedMeanAction(ScalarAction): 

108 """Calculates the weighted mean of a vector using a second vector as 

109 weights.""" 

110 

111 vectorKey = Field[str]("Key of Vector of values to compute the weighted mean of.") 

112 weightsKey = Field[str]("Key of Vector of weights.") 

113 

114 def getInputSchema(self) -> KeyedDataSchema: 

115 return ((self.vectorKey, Vector), (self.weightsKey, Vector)) 

116 

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

118 mask = self.getMask(**kwargs) 

119 values = _dataToArray(data[self.vectorKey.format(**kwargs)])[mask] 

120 weights = _dataToArray(data[self.weightsKey.format(**kwargs)])[mask] 

121 valid = ~np.isnan(values) 

122 total_weight = np.sum(weights[valid]) 

123 if total_weight == 0: 123 ↛ 124line 123 didn't jump to line 124 because the condition on line 123 was never true

124 return cast(Scalar, np.nan) 

125 return cast(Scalar, np.average(values[valid], weights=weights[valid])) 

126 

127 

128class StdevAction(ScalarFromVectorAction): 

129 """Calculates the standard deviation of the given data.""" 

130 

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

132 mask = self.getMask(**kwargs) 

133 return nanStd(_dataToArray(data[self.vectorKey.format(**kwargs)])[mask]) 

134 

135 

136class RmsAction(ScalarFromVectorAction): 

137 """Calculates the root mean square of the given data (without subtracting 

138 the mean as in StdevAction).""" 

139 

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

141 mask = self.getMask(**kwargs) 

142 vector = _dataToArray(data[self.vectorKey.format(**kwargs)])[mask] 

143 vector = vector[~np.isnan(vector)] 

144 

145 return np.sqrt(np.mean(vector**2)) 

146 

147 

148class ValueAction(ScalarFromVectorAction): 

149 """Extracts the first value from a vector.""" 

150 

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

152 return cast(Scalar, float(data[self.vectorKey.format(**kwargs)][0])) 

153 

154 

155class SigmaMadAction(ScalarFromVectorAction): 

156 """Calculates the sigma mad of the given data.""" 

157 

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

159 mask = self.getMask(**kwargs) 

160 values = _dataToArray(data[self.vectorKey.format(**kwargs)])[mask] 

161 return nanSigmaMad(values) 

162 

163 

164class CountAction(ScalarAction): 

165 """Performs count actions, with threshold-based filtering. 

166 The operator is specified as a string, for example, "lt", "le", "ge", 

167 "gt", "ne", and "eq" for the mathematical operations <, <=, >=, >, !=, 

168 and == respectively. To count non-NaN values, only pass the column name 

169 as vector key. To count NaN values, pass threshold = nan (from math.nan). 

170 Optionally to configure from a YAML file, pass "threshold: !!float nan". 

171 To compute the number of elements with values less than a given threshold, 

172 use op="le". 

173 """ 

174 

175 vectorKey = Field[str]("Key of Vector to count") 

176 op = ChoiceField[str]( 

177 doc="Operator name string.", 

178 allowed={ 

179 "lt": "less than threshold", 

180 "le": "less than or equal to threshold", 

181 "ge": "greater than or equal to threshold", 

182 "ne": "not equal to a given value", 

183 "eq": "equal to a given value", 

184 "gt": "greater than threshold", 

185 }, 

186 default="ne", 

187 ) 

188 threshold = Field[float](doc="Threshold to apply.", default=nan) 

189 

190 def getInputSchema(self) -> KeyedDataSchema: 

191 return ((self.vectorKey, Vector),) 

192 

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

194 mask = self.getMask(**kwargs) 

195 arr = _dataToArray(data[self.vectorKey.format(**kwargs)])[mask] 

196 

197 # Count NaNs and non-NaNs 

198 if self.threshold == nan: 198 ↛ 199line 198 didn't jump to line 199 because the condition on line 198 was never true

199 if self.op == "eq": 

200 # Count number of NaNs 

201 result = np.isnan(arr).sum() 

202 return cast(Scalar, int(result)) 

203 elif self.op == "ne": 

204 # Count number of non-NaNs 

205 result = arr.size - np.isnan(arr).sum() 

206 return cast(Scalar, int(result)) 

207 else: 

208 raise ValueError("Invalid operator for counting NaNs.") 

209 # Count for given threshold ignoring all NaNs 

210 else: 

211 result = arr[~np.isnan(arr)] 

212 result = cast( 

213 Scalar, 

214 int(np.sum(getattr(operator, self.op)(result, self.threshold))), 

215 ) 

216 return result 

217 

218 

219class CountUniqueAction(ScalarFromVectorAction): 

220 """Counts the number of unique rows in a given column.""" 

221 

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

223 mask = self.getMask(**kwargs) 

224 values = _dataToArray(data[self.vectorKey.format(**kwargs)])[mask] 

225 count = np.unique(values).size 

226 return cast(Scalar, count) 

227 

228 

229class ApproxFloor(ScalarFromVectorAction): 

230 """Returns the median of the lowest ten values of the sorted input.""" 

231 

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

233 mask = self.getMask(**kwargs) 

234 values = np.sort(_dataToArray(data[self.vectorKey.format(**kwargs)])[mask], axis=None) # type: ignore 

235 x = values.size // 10 

236 return nanMedian(values[-x:]) 

237 

238 

239class FracThreshold(ScalarFromVectorAction): 

240 """Compute the fraction of a distribution above or below a threshold. 

241 

242 The operator is specified as a string, for example, 

243 "lt", "le", "ge", "gt" for the mathematical operations <, <=, >=, >. To 

244 compute the fraction of elements with values less than a given threshold, 

245 use op="le". 

246 """ 

247 

248 op = ChoiceField[str]( 

249 doc="Operator name string.", 

250 allowed={ 

251 "lt": "less than threshold", 

252 "le": "less than or equal to threshold", 

253 "ge": "greater than or equal to threshold", 

254 "gt": "greater than threshold", 

255 }, 

256 ) 

257 threshold = Field[float](doc="Threshold to apply.") 

258 percent = Field[bool](doc="Express result as percentage", default=False) 

259 relative_to_median = Field[bool](doc="Calculate threshold relative to the median?", default=False) 

260 use_absolute_value = Field[bool]( 

261 doc=( 

262 "Calculate threshold after taking absolute value. If relative_to_median" 

263 " is true the absolute value will be applied after the median is subtracted" 

264 ), 

265 default=False, 

266 ) 

267 

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

269 mask = self.getMask(**kwargs) 

270 values = _dataToArray(data[self.vectorKey.format(**kwargs)])[mask] 

271 values = values[np.logical_not(np.isnan(values))] 

272 n_values = values.size 

273 if n_values == 0: 

274 return np.nan 

275 threshold = self.threshold 

276 # If relative_to_median is set, shift the threshold to be median+thresh 

277 if self.relative_to_median and values.size > 0: 

278 offset = nanMedian(values) 

279 if np.isfinite(offset): 

280 values -= offset 

281 if self.use_absolute_value: 

282 values = np.abs(values) 

283 result = cast( 

284 Scalar, 

285 float(np.sum(getattr(operator, self.op)(values, threshold)) / n_values), # type: ignore 

286 ) 

287 if self.percent: 

288 return 100.0 * result 

289 else: 

290 return result 

291 

292 

293class MaxAction(ScalarFromVectorAction): 

294 """Returns the maximum of the given data.""" 

295 

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

297 mask = self.getMask(**kwargs) 

298 return nanMax(_dataToArray(data[self.vectorKey.format(**kwargs)])[mask]) 

299 

300 

301class MinAction(ScalarFromVectorAction): 

302 """Returns the minimum of the given data.""" 

303 

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

305 mask = self.getMask(**kwargs) 

306 return nanMin(_dataToArray(data[self.vectorKey.format(**kwargs)])[mask]) 

307 

308 

309class FracInRange(ScalarFromVectorAction): 

310 """Compute the fraction of a distribution that is between specified 

311 minimum and maximum values, and is not NaN. 

312 """ 

313 

314 maximum = Field[float](doc="The maximum value", default=np.nextafter(np.inf, 0.0)) 

315 minimum = Field[float](doc="The minimum value", default=np.nextafter(-np.inf, 0.0)) 

316 percent = Field[bool](doc="Express result as percentage", default=False) 

317 

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

319 mask = self.getMask(**kwargs) 

320 values = _dataToArray(data[self.vectorKey.format(**kwargs)])[mask] 

321 nvalues = values.size 

322 values = values[np.logical_not(np.isnan(values))] 

323 sel_range = (values >= self.minimum) & (values < self.maximum) 

324 result = cast( 

325 Scalar, 

326 float(values[sel_range].size / nvalues), # type: ignore 

327 ) 

328 if self.percent: 

329 return 100.0 * result 

330 else: 

331 return result 

332 

333 

334class FracNan(ScalarFromVectorAction): 

335 """Compute the fraction of vector entries that are NaN.""" 

336 

337 percent = Field[bool](doc="Express result as percentage", default=False) 

338 

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

340 mask = self.getMask(**kwargs) 

341 values = _dataToArray(data[self.vectorKey.format(**kwargs)])[mask] 

342 nvalues = values.size 

343 values = values[np.isnan(values)] 

344 result = cast( 

345 Scalar, 

346 float(values.size / nvalues), # type: ignore 

347 ) 

348 if self.percent: 

349 return 100.0 * result 

350 else: 

351 return result 

352 

353 

354class SumAction(ScalarFromVectorAction): 

355 """Returns the sum of all values in the column.""" 

356 

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

358 mask = self.getMask(**kwargs) 

359 arr = _dataToArray(data[self.vectorKey.format(**kwargs)])[mask] 

360 return cast(Scalar, np.nansum(arr)) 

361 

362 

363class MedianHistAction(ScalarAction): 

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

365 

366 histKey = Field[str]("Key of frequency Vector") 

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

368 

369 def getInputSchema(self) -> KeyedDataSchema: 

370 return ( 

371 (self.histKey, Vector), 

372 (self.midKey, Vector), 

373 ) 

374 

375 def histMedian(self, hist, bin_mid): 

376 """Calculates the median of a histogram with binned values 

377 

378 Parameters 

379 ---------- 

380 hist : `numpy.ndarray` 

381 Frequency array 

382 bin_mid : `numpy.ndarray` 

383 Bin midpoints array 

384 

385 Returns 

386 ------- 

387 median : `float` 

388 Median of histogram with binned values 

389 """ 

390 cumulative_sum = np.cumsum(hist) 

391 median_index = np.searchsorted(cumulative_sum, cumulative_sum[-1] / 2) 

392 median = bin_mid[median_index] 

393 return median 

394 

395 def __call__(self, data: KeyedData, **kwargs): 

396 hist = _dataToArray(data[self.histKey.format(**kwargs)]) 

397 if hist.size != 0: 

398 bin_mid = _dataToArray(data[self.midKey.format(**kwargs)]) 

399 med = cast(Scalar, float(self.histMedian(hist, bin_mid))) 

400 else: 

401 med = np.nan 

402 return med 

403 

404 

405class IqrHistAction(ScalarAction): 

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

407 

408 histKey = Field[str]("Key of frequency Vector") 

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

410 

411 def getInputSchema(self) -> KeyedDataSchema: 

412 return ( 

413 (self.histKey, Vector), 

414 (self.midKey, Vector), 

415 ) 

416 

417 def histIqr(self, hist, bin_mid): 

418 """Calculates the interquartile range of a histogram with binned values 

419 

420 Parameters 

421 ---------- 

422 hist : `numpy.ndarray` 

423 Frequency array 

424 bin_mid : `numpy.ndarray` 

425 Bin midpoints array 

426 

427 Returns 

428 ------- 

429 iqr : `float` 

430 Inter-quartile range of histogram with binned values 

431 """ 

432 cumulative_sum = np.cumsum(hist) 

433 liqr_index = np.searchsorted(cumulative_sum, cumulative_sum[-1] / 4) 

434 uiqr_index = np.searchsorted(cumulative_sum, (3 / 4) * cumulative_sum[-1]) 

435 liqr = bin_mid[liqr_index] 

436 uiqr = bin_mid[uiqr_index] 

437 iqr = uiqr - liqr 

438 return iqr 

439 

440 def __call__(self, data: KeyedData, **kwargs): 

441 hist = _dataToArray(data[self.histKey.format(**kwargs)]) 

442 if hist.size != 0: 

443 bin_mid = _dataToArray(data[self.midKey.format(**kwargs)]) 

444 iqr = cast(Scalar, float(self.histIqr(hist, bin_mid))) 

445 else: 

446 iqr = np.nan 

447 return iqr 

448 

449 

450class DivideScalar(ScalarAction): 

451 """Calculate (A/B) for scalars.""" 

452 

453 actionA = ConfigurableActionField[ScalarAction](doc="Action which supplies scalar A") 

454 actionB = ConfigurableActionField[ScalarAction](doc="Action which supplies scalar B") 

455 

456 def getInputSchema(self) -> KeyedDataSchema: 

457 yield from self.actionA.getInputSchema() 

458 yield from self.actionB.getInputSchema() 

459 

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

461 """Return the result of A/B. 

462 

463 Parameters 

464 ---------- 

465 data : `KeyedData` 

466 

467 Returns 

468 ------- 

469 result : `Scalar` 

470 The result of dividing A by B. 

471 """ 

472 scalarA = self.actionA(data, **kwargs) 

473 scalarB = self.actionB(data, **kwargs) 

474 if scalarB == 0: 474 ↛ 475line 474 didn't jump to line 475 because the condition on line 474 was never true

475 if scalarA == 0: 

476 log.warning("Both numerator and denominator are zero! Returning NaN.") 

477 return np.nan 

478 else: 

479 value = np.sign(scalarA) * np.inf 

480 log.warning("Non-zero scalar divided by zero! Returning %f.", value) 

481 return value 

482 else: 

483 return scalarA / scalarB 

484 

485 

486class MedianGradientAction(ScalarAction): 

487 """Calculate the gradient of a running median""" 

488 

489 lowerBinLimit = Field[float](doc="Percentile of data to start the bining at", default=5.0) 

490 upperBinLimit = Field[float](doc="Percentile of data to end the binning at", default=95.0) 

491 nBins = Field[int](doc="Number of bins to use for running median", default=50) 

492 xsVectorKey = Field[str]("Key of Vector that gives the x location of the points.") 

493 ysVectorKey = Field[str]("Key of the Vector to compute the statistic from.") 

494 

495 def getInputSchema(self) -> KeyedDataSchema: 

496 return ((self.xsVectorKey, Vector), (self.ysVectorKey, Vector)) 

497 

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

499 """Return the gradient of the running median. 

500 

501 Parameters 

502 ---------- 

503 data : `KeyedData` 

504 

505 Returns 

506 ------- 

507 result : `Scalar` 

508 The gradient of the running median 

509 """ 

510 

511 xs = _dataToArray(data[self.xsVectorKey.format(**kwargs)]) 

512 ys = _dataToArray(data[self.ysVectorKey.format(**kwargs)]) 

513 

514 lowerLim = np.nanpercentile(xs, self.lowerBinLimit) 

515 upperLim = np.nanpercentile(xs, self.upperBinLimit) 

516 

517 use = (xs > lowerLim) & (xs < upperLim) & np.isfinite(xs) & np.isfinite(ys) 

518 

519 if np.sum(use) > 2: 

520 meds, binEdges, _ = binned_statistic(xs[use], ys[use], statistic="median", bins=self.nBins) 

521 else: 

522 return np.nan 

523 

524 def func(x, m, b): 

525 return m * x + b 

526 

527 finiteMeds = np.isfinite(meds) 

528 popt, _ = curve_fit(func, binEdges[:-1][finiteMeds], meds[finiteMeds]) 

529 

530 return popt[0]