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

1# This file is part of analysis_tools.

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.

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

22from __future__ import annotations

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)

47import logging

48import operator

49from math import nan

50from typing import cast

52import numpy as np

53from scipy.optimize import curve_fit

54from scipy.stats import binned_statistic

56from lsst.pex.config import ChoiceField, Field

57from lsst.pex.config.configurableActions import ConfigurableActionField

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

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

62log = logging.getLogger(__name__)

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)

76class ScalarFromVectorAction(ScalarAction):

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

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

81 def getInputSchema(self) -> KeyedDataSchema:

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

85class MedianAction(ScalarFromVectorAction):

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

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

93 return med

96class MeanAction(ScalarFromVectorAction):

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

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]