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 "FullRangeAction", 

37 "FracInRange", 

38 "FracNan", 

39 "SumAction", 

40 "WeightedMeanAction", 

41 "MedianHistAction", 

42 "IqrHistAction", 

43 "DivideScalar", 

44 "RmsAction", 

45 "MedianGradientAction", 

46) 

47 

48import logging 

49import operator 

50from math import nan 

51from typing import cast 

52 

53import numpy as np 

54from scipy.optimize import curve_fit 

55from scipy.stats import binned_statistic 

56 

57from lsst.pex.config import ChoiceField, Field 

58from lsst.pex.config.configurableActions import ConfigurableActionField 

59 

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

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

62 

63log = logging.getLogger(__name__) 

64 

65 

66def _dataToArray(data): 

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

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

69 where possible. 

70 """ 

71 try: 

72 return np.from_dlpack(data) 

73 except (AttributeError, BufferError): 

74 return np.array(data) 

75 

76 

77class ScalarFromVectorAction(ScalarAction): 

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

79 

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

81 

82 def getInputSchema(self) -> KeyedDataSchema: 

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

84 

85 

86class MedianAction(ScalarFromVectorAction): 

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

88 

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

90 mask = self.getMask(**kwargs) 

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

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

93 

94 return med 

95 

96 

97class MeanAction(ScalarFromVectorAction): 

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

99 

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

101 mask = self.getMask(**kwargs) 

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

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

104 

105 return mean 

106 

107 

108class WeightedMeanAction(ScalarAction): 

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

110 weights.""" 

111 

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

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

114 

115 def getInputSchema(self) -> KeyedDataSchema: 

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

117 

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

119 mask = self.getMask(**kwargs) 

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

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

122 valid = ~np.isnan(values) 

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

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

125 return cast(Scalar, np.nan) 

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

127 

128 

129class StdevAction(ScalarFromVectorAction): 

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

131 

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

133 mask = self.getMask(**kwargs) 

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

135 

136 

137class RmsAction(ScalarFromVectorAction): 

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

139 the mean as in StdevAction).""" 

140 

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

142 mask = self.getMask(**kwargs) 

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

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

145 

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

147 

148 

149class ValueAction(ScalarFromVectorAction): 

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

151 

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

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

154 

155 

156class SigmaMadAction(ScalarFromVectorAction): 

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

158 

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

160 mask = self.getMask(**kwargs) 

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

162 return nanSigmaMad(values) 

163 

164 

165class CountAction(ScalarAction): 

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

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

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

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

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

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

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

173 use op="le". 

174 """ 

175 

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

177 op = ChoiceField[str]( 

178 doc="Operator name string.", 

179 allowed={ 

180 "lt": "less than threshold", 

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

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

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

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

185 "gt": "greater than threshold", 

186 }, 

187 default="ne", 

188 ) 

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

190 

191 def getInputSchema(self) -> KeyedDataSchema: 

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

193 

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

195 mask = self.getMask(**kwargs) 

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

197 

198 # Count NaNs and non-NaNs 

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

200 if self.op == "eq": 

201 # Count number of NaNs 

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

203 return cast(Scalar, int(result)) 

204 elif self.op == "ne": 

205 # Count number of non-NaNs 

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

207 return cast(Scalar, int(result)) 

208 else: 

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

210 # Count for given threshold ignoring all NaNs 

211 else: 

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

213 result = cast( 

214 Scalar, 

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

216 ) 

217 return result 

218 

219 

220class CountUniqueAction(ScalarFromVectorAction): 

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

222 

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

224 mask = self.getMask(**kwargs) 

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

226 count = np.unique(values).size 

227 return cast(Scalar, count) 

228 

229 

230class ApproxFloor(ScalarFromVectorAction): 

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

232 

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

234 mask = self.getMask(**kwargs) 

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

236 x = values.size // 10 

237 return nanMedian(values[-x:]) 

238 

239 

240class FracThreshold(ScalarFromVectorAction): 

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

242 

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

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

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

246 use op="le". 

247 """ 

248 

249 op = ChoiceField[str]( 

250 doc="Operator name string.", 

251 allowed={ 

252 "lt": "less than threshold", 

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

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

255 "gt": "greater than threshold", 

256 }, 

257 ) 

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

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

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

261 use_absolute_value = Field[bool]( 

262 doc=( 

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

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

265 ), 

266 default=False, 

267 ) 

268 

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

270 mask = self.getMask(**kwargs) 

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

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

273 n_values = values.size 

274 if n_values == 0: 

275 return np.nan 

276 threshold = self.threshold 

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

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

279 offset = nanMedian(values) 

280 if np.isfinite(offset): 

281 values -= offset 

282 if self.use_absolute_value: 

283 values = np.abs(values) 

284 result = cast( 

285 Scalar, 

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

287 ) 

288 if self.percent: 

289 return 100.0 * result 

290 else: 

291 return result 

292 

293 

294class MaxAction(ScalarFromVectorAction): 

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

296 

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

298 mask = self.getMask(**kwargs) 

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

300 

301 

302class MinAction(ScalarFromVectorAction): 

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

304 

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

306 mask = self.getMask(**kwargs) 

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

308 

309 

310class FullRangeAction(ScalarFromVectorAction): 

311 """Returns the full range (i.e., max-min) of the given data.""" 

312 

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

314 mask = self.getMask(**kwargs) 

315 return nanMax(_dataToArray(data[self.vectorKey.format(**kwargs)][mask])) - nanMin( 

316 _dataToArray(data[self.vectorKey.format(**kwargs)][mask]) 

317 ) 

318 

319 

320class FracInRange(ScalarFromVectorAction): 

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

322 minimum and maximum values, and is not NaN. 

323 """ 

324 

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

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

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

328 

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

330 mask = self.getMask(**kwargs) 

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

332 nvalues = values.size 

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

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

335 result = cast( 

336 Scalar, 

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

338 ) 

339 if self.percent: 

340 return 100.0 * result 

341 else: 

342 return result 

343 

344 

345class FracNan(ScalarFromVectorAction): 

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

347 

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

349 

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

351 mask = self.getMask(**kwargs) 

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

353 nvalues = values.size 

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

355 result = cast( 

356 Scalar, 

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

358 ) 

359 if self.percent: 

360 return 100.0 * result 

361 else: 

362 return result 

363 

364 

365class SumAction(ScalarFromVectorAction): 

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

367 

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

369 mask = self.getMask(**kwargs) 

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

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

372 

373 

374class MedianHistAction(ScalarAction): 

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

376 

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

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

379 

380 def getInputSchema(self) -> KeyedDataSchema: 

381 return ( 

382 (self.histKey, Vector), 

383 (self.midKey, Vector), 

384 ) 

385 

386 def histMedian(self, hist, bin_mid): 

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

388 

389 Parameters 

390 ---------- 

391 hist : `numpy.ndarray` 

392 Frequency array 

393 bin_mid : `numpy.ndarray` 

394 Bin midpoints array 

395 

396 Returns 

397 ------- 

398 median : `float` 

399 Median of histogram with binned values 

400 """ 

401 cumulative_sum = np.cumsum(hist) 

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

403 median = bin_mid[median_index] 

404 return median 

405 

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

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

408 if hist.size != 0: 

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

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

411 else: 

412 med = np.nan 

413 return med 

414 

415 

416class IqrHistAction(ScalarAction): 

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

418 

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

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

421 

422 def getInputSchema(self) -> KeyedDataSchema: 

423 return ( 

424 (self.histKey, Vector), 

425 (self.midKey, Vector), 

426 ) 

427 

428 def histIqr(self, hist, bin_mid): 

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

430 

431 Parameters 

432 ---------- 

433 hist : `numpy.ndarray` 

434 Frequency array 

435 bin_mid : `numpy.ndarray` 

436 Bin midpoints array 

437 

438 Returns 

439 ------- 

440 iqr : `float` 

441 Inter-quartile range of histogram with binned values 

442 """ 

443 cumulative_sum = np.cumsum(hist) 

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

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

446 liqr = bin_mid[liqr_index] 

447 uiqr = bin_mid[uiqr_index] 

448 iqr = uiqr - liqr 

449 return iqr 

450 

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

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

453 if hist.size != 0: 

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

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

456 else: 

457 iqr = np.nan 

458 return iqr 

459 

460 

461class DivideScalar(ScalarAction): 

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

463 

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

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

466 

467 def getInputSchema(self) -> KeyedDataSchema: 

468 yield from self.actionA.getInputSchema() 

469 yield from self.actionB.getInputSchema() 

470 

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

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

473 

474 Parameters 

475 ---------- 

476 data : `KeyedData` 

477 

478 Returns 

479 ------- 

480 result : `Scalar` 

481 The result of dividing A by B. 

482 """ 

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

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

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

486 if scalarA == 0: 

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

488 return np.nan 

489 else: 

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

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

492 return value 

493 else: 

494 return scalarA / scalarB 

495 

496 

497class MedianGradientAction(ScalarAction): 

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

499 

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

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

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

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

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

505 

506 def getInputSchema(self) -> KeyedDataSchema: 

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

508 

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

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

511 

512 Parameters 

513 ---------- 

514 data : `KeyedData` 

515 

516 Returns 

517 ------- 

518 result : `Scalar` 

519 The gradient of the running median 

520 """ 

521 

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

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

524 

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

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

527 

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

529 

530 if np.sum(use) > 2: 

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

532 else: 

533 return np.nan 

534 

535 def func(x, m, b): 

536 return m * x + b 

537 

538 finiteMeds = np.isfinite(meds) 

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

540 

541 return popt[0]