Coverage for uqmodels/postprocessing/anomaly_processing.py: 42%

1#################################################################### 

2# Ensemble of processing aim to process UQmesures into Intermediate quantity or into Anom-KPI 

3 

4 

5import numpy as np 

6import scipy 

7from sklearn.covariance import EmpiricalCovariance 

8 

9import uqmodels.postprocessing.UQ_processing as UQ_proc 

10from uqmodels.utils import apply_conv, apply_middledim_reduction 

11 

12 

13def score_seuil(s, per_seuil=0.995, local=True): 

14 """Thresholding of extrem values by percentile by dimension if local = True or globally else. 

15 

16 Args: 

17 s (_type_): _description_ 

18 per_seuil (float, optional): _description_. Defaults to 0.995. 

19 local (bool, optional): _description_. Defaults to True. 

20 

21 Returns: 

22 _type_: _description_ 

23 """ 

24 

25 if local: 

26 p_s = np.quantile(np.abs(s), per_seuil, axis=0) 

27 else: 

28 p_s = np.quantile(np.abs(s), per_seuil) 

29 

30 # Seuillage 

31 

32 s = np.sign(s) * np.minimum(np.abs(s), p_s) 

33 

34 return s 

35 

36 

37# Convolution of the signal "s" with a filter "filt" 

38 

39 

40def fit_calibrate( 

41 residuals, 

42 ctx_mask=None, 

43 beta=0.05, 

44 type_norm=None, 

45 empiric_rescale=False, 

46 deg=1, 

47 d=1, 

48 debug=False, 

49): 

50 """Function that estimate calibrate parameters form residu link to some hypothesis 

51 

52 Args: 

53 s (_type_): score 

54 ctx_mask : mask for ctx calib 

55 beta (_type_): Anomaly target ratio 

56 type_norm (_type_, optional): assymption 

57 empiric_rescale (bool, optional): _description_. Defaults to True. 

58 deg (int, optional): _description_. Defaults to 1. 

59 

60 Type_norm : 

61 "Nsigma_local" : Fexible anomalie thresold (N-sgima) by dimension | Hypothesis : 

62 Anomaly distributed homegeneously by dimension. 

63 "Nsigma_global" : Fexible anomalie thresold (N-sgima) global : homogeneous anomaly distribution by dimension. 

64 "Chi2": Theorical Chi2 correction 

65 "quantiles_local" : rigid anomalie thresold (quantiles) by dimension : 

66 heterogeneos anomaly distrubtion by dimension. 

67 "quantiles_global" : (Default) rigid anomalie thresold (quantiles) global : 

68 homogeneous anomaly distribution by dimension. 

69 "None" : No normalisation 

70 

71 Returns: 

72 params : (list_norm_val : calibration coeficient for each ctx_values 

73 list_ctx_id : list of ctx_values) 

74 """ 

75 if ctx_mask is None: 

76 ctx_mask = np.zeros(len(residuals)) 

77 list_ctx_id = np.sort(list(set(ctx_mask))) 

78 

79 if debug: 

80 print( 

81 "fit_calib Start", 

82 (residuals < -1).mean(), 

83 (residuals > 1).mean(), 

84 "type_norm", 

85 type_norm, 

86 ) 

87 

88 list_norm_val = [] 

89 for v_ctx in list_ctx_id: 

90 mask = ctx_mask == v_ctx 

91 # Fexible anomalie thresold (N-sgima) by dimension | Hypothesis : Anomaly 

92 # distributed homegeneously by dimension. 

93 if type_norm == "Nsigma_local": 

94 f_corr = 1 

95 if empiric_rescale: 

96 f_corr = residuals[mask].std(axis=0) 

97 threshold = np.power(scipy.stats.norm.ppf(1 - beta, 0, 1), 2) 

98 

99 # Fexible anomalie thresold (N-sgima) global : heterogenous anomaly distribution by dimension. 

100 elif type_norm == "Nsigma_global": 

101 f_corr = 1 

102 if empiric_rescale: 

103 f_corr = residuals[mask].std() 

104 

105 threshold = np.power(scipy.stats.norm.ppf(1 - beta, 0, 1), 2) 

106 

107 # Theorical Chi2 correction 

108 elif type_norm == "Chi2": 

109 f_corr = 1 

110 threshold = np.power(scipy.stats.chi2.ppf(1 - beta, deg), d / 2) 

111 

112 # Chebyshev's inequality 

113 elif type_norm == "chebyshev_local": 

114 f_corr = residuals[mask].std(axis=0) 

115 threshold = 1 / np.sqrt(beta) 

116 

117 # Chebyshev's inequality 

118 elif type_norm == "chebyshev_global": 

119 f_corr = residuals[mask].std() 

120 threshold = 1 / np.sqrt(beta) 

121 

122 # Cantelli's inequality 

123 elif type_norm == "Cantelli_local": 

124 f_corr = residuals[mask].std(axis=0) 

125 threshold = np.sqrt((1 - beta) / beta) 

126 

127 # Cantelli's inequality 

128 elif type_norm == "Cantelli_global": 

129 f_corr = residuals[mask].std() 

130 threshold = np.sqrt((1 - beta) / beta) 

131 

132 # rigid anomalie thresold (quantiles) by dimension : heterogeneos anomaly distrubtion by dimension. 

133 elif type_norm == "quantiles_local": 

134 f_corr = 1 

135 threshold = np.quantile(np.abs(residuals[mask]), 1 - beta, axis=0) 

136 

137 elif type_norm == "quantiles_global": 

138 f_corr = 1 

139 threshold = np.quantile(np.abs(residuals[mask]), 1 - beta) 

140 

141 else: 

142 f_corr = 1 

143 threshold = 1 

144 

145 list_norm_val.append(f_corr * threshold) 

146 

147 if debug: 

148 print( 

149 "debug fit_calibrate : Val_ctx", 

150 v_ctx, 

151 "target", 

152 beta, 

153 "v_mean", 

154 np.round(np.abs(residuals).mean(), 3), 

155 "emp", 

156 (np.abs(residuals[mask]) >= list_norm_val).mean(), 

157 list_norm_val, 

158 type_norm, 

159 ) 

160 

161 params_ = list_norm_val, list_ctx_id 

162 

163 return params_ 

164 

165 

166def compute_calibrate( 

167 residuals, 

168 ctx_mask=None, 

169 beta=0.05, 

170 type_norm=None, 

171 empiric_rescale=True, 

172 deg=1, 

173 params_=None, 

174 mode="score", 

175 debug=False, 

176): 

177 """Function that apply calibration on residu based on some hypothesis 

178 

179 Args: 

180 res (_type_): residu to calibrate 

181 ctx_mask (None): mask for ctx calib 

182 beta (_type_): Anomaly target ratio 

183 type_norm (_type_, optional): assymption 

184 empiric_rescale (bool, optional): _description_. Defaults to True. 

185 deg (int, optional): Power to apply to the residu. Defaults to 1. 

186 

187 type_norm assymption : 

188 "Nsigma_local" : Fexible anomalie thresold (N-sgima) by dimension | Hypothesis : 

189 Anomaly distributed homegeneously by dimension. 

190 "Nsigma_global" : Fexible anomalie thresold (N-sgima) global : homogeneous anomaly distribution by dimension. 

191 "Chi2": Theorical Chi2 correction 

192 "quantiles_local" : rigid anomalie thresold (quantiles) by dimension : 

193 heterogeneos anomaly distrubtion by dimension. 

194 "quantiles_global" : (Default) rigid anomalie thresold (quantiles) global : 

195 homogeneous anomaly distribution by dimension. 

196 "None" : No normalisation 

197 

198 

199 Returns: 

200 res : calibrated residu 

201 params: Params provided or computed 

202 """ 

203 if ctx_mask is None: 

204 ctx_mask = np.zeros(len(residuals)) 

205 

206 if params_ is None: 

207 params_ = fit_calibrate( 

208 residuals, 

209 ctx_mask=ctx_mask, 

210 beta=beta, 

211 type_norm=type_norm, 

212 empiric_rescale=empiric_rescale, 

213 deg=deg, 

214 debug=debug, 

215 ) 

216 

217 list_norm_val, list_ctx_id = params_ 

218 

219 if mode == "score": 

220 if debug: 

221 print( 

222 "compute_calib Start", 

223 (np.abs(residuals) > 1).mean(), 

224 "type_norm", 

225 type_norm, 

226 "res_mean", 

227 np.abs(residuals).mean(), 

228 list_norm_val, 

229 (np.abs(residuals) > list_norm_val[0]).mean(), 

230 ) 

231 

232 for n, ctx_val in enumerate(list_ctx_id): 

233 mask = ctx_mask == ctx_val 

234 residuals[mask].shape, list_norm_val[n].shape 

235 if mode == "score": 

236 

237 residuals[mask] = residuals[mask] / list_norm_val[n] 

238 

239 if mode == "born": 

240 residuals[mask] = residuals[mask] * list_norm_val[n] 

241 

242 if mode == "score": 

243 if debug: 

244 print( 

245 "compute_calib out", 

246 (np.abs(residuals) > 1).mean(), 

247 "type_norm", 

248 type_norm, 

249 "res_mean", 

250 np.abs(residuals).mean(), 

251 ) 

252 

253 return residuals, params_ 

254 

255 

256def fit_born_calibrated( 

257 UQ, 

258 type_UQ, 

259 pred, 

260 y, 

261 type_UQ_params=None, 

262 beta=0.1, 

263 type_res="res", 

264 ctx_mask=None, 

265 var_min=0, 

266 min_cut=0, 

267 max_cut=0, 

268 q_var=1, 

269 empiric_rescale=True, 

270): 

271 """!!!!Depreciated !!! Estimate calibration parameters in order to calibrate born from UQ measure, pred, observation and target 

272 

273 Args: 

274 UQ (np.array or list): UQmeasure obtain from UQEstimator 

275 type_UQ (_type_): Type UQ that the nature of UQmeasure 

276 pred (np.array): prediction provide by a predictor or an UQEstimator 

277 y (np.array): Targets/Observation 

278 type_UQ_params : additional parameters link to type paradigm (ex : alpha for quantile) 

279 beta: target miss-coverage of the borns 

280 min_cut (_type_): Bottom extremun percentile. 

281 max_cut (_type_): Bottom upper percentile. 

282 q_var (_type_): Power coefficent 

283 sigma_min (float, optional): Minimum values of UQ considered. 

284 

285 Returns: 

286 params = (list_norm_val, : list of calibration term for each ctx_values, 

287 list_ctx_id : list of ctx_values) 

288 """ 

289 if type_res == "res": 

290 residuals = y - pred 

291 

292 if type_res == "w_res": 

293 residuals = UQ_proc.process_UQmeasure_to_residu( 

294 UQ, 

295 type_UQ, 

296 pred, 

297 y, 

298 type_UQ_params=type_UQ_params, 

299 min_cut=min_cut, 

300 max_cut=max_cut, 

301 var_min=var_min, 

302 q_var=q_var, 

303 ) 

304 

305 elif type_res == "cqr": 

306 y_pred_lower = UQ_proc.process_UQmeasure_to_quantile( 

307 UQ, 

308 type_UQ, 

309 pred, 

310 y=None, 

311 type_UQ_params=type_UQ_params, 

312 alpha=beta / 2, 

313 var_min=0, 

314 min_cut=0, 

315 max_cut=1, 

316 q_var=1, 

317 ) 

318 

319 y_pred_upper = UQ_proc.process_UQmeasure_to_quantile( 

320 UQ, 

321 type_UQ, 

322 pred, 

323 y=None, 

324 type_UQ_params=type_UQ_params, 

325 alpha=1 - beta / 2, 

326 var_min=0, 

327 min_cut=0, 

328 max_cut=1, 

329 q_var=1, 

330 ) 

331 

332 residuals = np.maximum(y_pred_lower - y, y - y_pred_upper) * -( 

333 y_pred_lower - y 

334 ) > (y - y_pred_upper) 

335 

336 elif type_res == "no_calib": 

337 return ([1], None) 

338 

339 else: 

340 raise ValueError("type_res" + type_res + " not covered") 

341 

342 list_norm_val, list_ctx_id = fit_calibrate( 

343 residuals, 

344 ctx_mask=ctx_mask, 

345 beta=beta, 

346 type_norm="quantiles_local", 

347 empiric_rescale=empiric_rescale, 

348 deg=1, 

349 ) 

350 

351 # To do asymetrical : calibration distinct for negative and positive residuals 

352 # Add symetric or asymetric parameters mode 

353 # -> to implement in calibrate_residuals 

354 # list_norm_val_neg&pos, list_ctx_id = calibrate_residuals(residuals, 

355 # ctx_mask=None, alpha=alpha, type_norm="quantiles_local", 

356 # empiric_rescale=True, deg=1) 

357 params_ = (list_norm_val, list_ctx_id) 

358 return params_ 

359 

360 

361def compute_born_calibrated( 

362 UQ, 

363 type_UQ, 

364 pred, 

365 y=None, 

366 type_UQ_params=None, 

367 beta=0.1, 

368 type_res="res", 

369 ctx_mask=None, 

370 params_=None, 

371): 

372 """!!!!Depreciated !!! Compute_UQ_calibration from UQ measure, pred, observation and target 

373 

374 Args: 

375 UQ (np.array or list): UQmeasure obtain from UQEstimator 

376 type_UQ (str): Type UQ that the nature of UQmeasure 

377 pred (np.array): prediction provide by a predictor or an UQEstimator 

378 y (np.array): Targets/Observation 

379 type_UQ_params (_type_, optional): additional parameters link to type paradigm (ex : alpha for quantile) 

380 beta (float, optional): target miss-coverage of the borns 

381 type_res : 

382 "no_calib : No calibration 

383 "res" : Calibration based on mean residuals 

384 "w_res" : Calibration based on weigthed mean residuals 

385 "cqr" : Calibration based on quantile residuals 

386 ctx_mask (_type_, optional): _description_. Defaults to None. 

387 params (_type_, optional): _description_. Defaults to None. 

388 

389 Raises: 

390 ValueError: _description_ 

391 

392 Returns: 

393 (y_pred_lower, y_pred_upper), params : (upper & lower calibrated bound) & compute params 

394 """ 

395 if ctx_mask is None: 

396 ctx_mask = np.zeros(len(y)) 

397 

398 if params_ is None: 

399 params_ = fit_born_calibrated( 

400 UQ, 

401 type_UQ, 

402 pred, 

403 y, 

404 type_UQ_params, 

405 beta=beta, 

406 type_res=type_res, 

407 ctx_mask=ctx_mask, 

408 ) 

409 

410 list_norm_val, list_ctx_id = params_ 

411 

412 if type_res == "res": 

413 UQ = pred * 0 + 1 

414 type_UQ = "var" 

415 

416 y_pred_lower = UQ_proc.process_UQmeasure_to_quantile( 

417 UQ, 

418 type_UQ, 

419 pred, 

420 y=None, 

421 type_UQ_params=type_UQ_params, 

422 alpha=beta / 2, 

423 var_min=0, 

424 min_cut=0, 

425 max_cut=1, 

426 q_var=1, 

427 ) 

428 

429 y_pred_upper = UQ_proc.process_UQmeasure_to_quantile( 

430 UQ, 

431 type_UQ, 

432 pred, 

433 y=None, 

434 type_UQ_params=type_UQ_params, 

435 alpha=1 - beta / 2, 

436 var_min=0, 

437 min_cut=0, 

438 max_cut=1, 

439 q_var=1, 

440 ) 

441 

442 for n, ctx_val in list_ctx_id: 

443 mask = ctx_mask == ctx_val 

444 if type_res in ["res", "w_res", "no_calib"]: 

445 y_pred_lower[mask] = ( 

446 (y_pred_lower[mask] - pred[mask]) * list_norm_val[n] 

447 ) + pred[mask] 

448 y_pred_upper[mask] = ( 

449 (y_pred_upper[mask] - pred[mask]) * list_norm_val[n] 

450 ) + pred[mask] 

451 

452 elif type_res in ["cqr"]: 

453 y_pred_lower[mask] = y_pred_lower[mask] - list_norm_val[n] 

454 y_pred_upper[mask] = y_pred_upper[mask] + list_norm_val[n] 

455 else: 

456 raise ValueError( 

457 "compute_UQ_calibration : type_res : " + type_res + " Not covered" 

458 ) 

459 # Hold assymetrical val 

460 

461 return (y_pred_lower, y_pred_upper), params_ 

462 

463 

464# FONCTION PRINCIPALE : CONSTRUCTION DU SCORE. 

465# res = Résidu source 

466# v_mean = biais source 

467# v_std = variance source. 

468# Train = mask lié à l'échantillon de test (Pour pouvoir adapater le seuil entre l'échantillon train et test) 

469 

470# Paramètre : 

471# Per_nom : ratio de donnée normale (local) à priori (par défault = 95%) 

472# beta : Target anomaly (par défault = 0.05) 

473# Per_seuil = seuillage percentage maximum sur le score d'anomalie.f 

474# min_cut = Valeur pour la coupe basse de la variance et du biais.100 

475# max_cut = Valeur pour la coupe haute de la variance et du biais. 

476# d = Mise en forme puissance du score final. 

477# Type_Norm : 

478# - True = Constante de normalisation selon une hypothèse de distribution d'anomalie homgène par dimension. 

479# - False = Constante de normalisation selon une hypothèse de distribution d'anomalie hétérogène par dimension. 

480# filt = filtre de convolution a appliquées. 

481# Q paramètre q lié à l'hypthèse de distrubution des anomalies par contexte : 

482# q < 1 : Anomalie lié au contexte de forte variance. 

483# q = 1 : Anomalie distribué de manière homogène par contexte 

484# q > 1 : Anomalie lié au contexte de faible variance. 

485# Dim = dimension de la série. 

486 

487 

488def fit_anom_score( 

489 UQ, 

490 type_UQ, 

491 pred, 

492 y, 

493 type_UQ_params=None, 

494 ctx_mask=None, 

495 beta=0.01, 

496 per_seuil=0.9995, 

497 min_cut=0.01, 

498 max_cut=0.97, 

499 var_min=0, 

500 d=2, 

501 q_var=1, 

502 k_var_e=1, 

503 q_var_e=1, 

504 q_Eratio=3, 

505 type_norm="quantiles_local", 

506 empiric_rescale=False, 

507 global_median_normalization=False, 

508 filt=None, 

509 reduc_filter=None, 

510 roll=0, 

511 debug=False, 

512 **kwargs 

513): 

514 """Fit parameters link to anomaly score calibration using Prediction and UQ measure. 

515 

516 Args: 

517 UQ (UQ-measure): UQ measure from UQ estimateurs 

518 type_UQ (str, optional): _description_. Defaults to "sigma". 

519 pred (array): prediction 

520 y (array): Real values 

521 type_UQ_params (str, optional): additional information about UQ. Defaults to None. 

522 ctx_mask : mask for ctx_calib 

523 beta (float, optional): Anom target Defaults to 0.05. 

524 per_seuil (float, optional): Threshold Defaults to 0.995. 

525 min_cut (float, optional): extremum cut Defaults to 0.005. 

526 max_cut (float, optional): extremum cut Defaults to 0.995. 

527 var_min (float, optional): Threshold of minimum var, Default to 0 

528 d (int, optional): residual power Defaults to 2. 

529 type_norm (str, optional): Type of normalisation see Defaults to "quantiles_local". 

530 global_median_normalization (bool,optional): Normalisation of residu by each step 

531 by the median residu of all dimension. 

532 filt (list, optional): _description_. Defaults to [0, 1, 0]. 

533 q_var (int, optional): _description_. Defaults to 1. 

534 

535 Returns: 

536 params : parameters to calibrate anomaly score 

537 """ 

538 # print('debug', UQ.shape, type_UQ, pred.shape, 

539 # y.shape, reduc_filter, filt, roll) 

540 

541 y, pred, UQ = UQ_proc.check_y_vs_pred_and_UQ_shape(y, pred, UQ) 

542 

543 ndUQ_ratio, extremum_var_TOT = None, (None, None) 

544 if type_UQ == "var_A&E": 

545 extremum_var_TOT, ndUQ_ratio = UQ_proc.get_extremum_var_TOT_and_ndUQ_ratio( 

546 UQ, 

547 type_UQ=type_UQ, 

548 min_cut=min_cut, 

549 max_cut=max_cut, 

550 var_min=0, 

551 var_max=None, 

552 factor=2, 

553 q_var=q_var, 

554 q_Eratio=q_Eratio, 

555 mode_multidim=True, 

556 ) 

557 

558 anom_score = UQ_proc.process_UQmeasure_to_residu( 

559 UQ, 

560 type_UQ, 

561 pred, 

562 y=y, 

563 type_UQ_params=type_UQ_params, 

564 d=d, 

565 min_cut=min_cut, 

566 max_cut=max_cut, 

567 q_var=q_var, 

568 var_min=var_min, 

569 k_var_e=k_var_e, 

570 q_var_e=q_var_e, 

571 q_Eratio=q_Eratio, 

572 ndUQ_ratio=ndUQ_ratio, 

573 extremum_var_TOT=extremum_var_TOT, 

574 reduc_filter=reduc_filter, 

575 roll=roll, 

576 debug=debug, 

577 ) 

578 

579 if global_median_normalization: 

580 global_norm = np.quantile(anom_score, 0.5, axis=1)[:, None] 

581 anom_score = anom_score - global_norm 

582 

583 # Dim reduction : 

584 

585 # Seuillage des résidu normalisé 

586 

587 anom_score = apply_conv(anom_score, filt) 

588 anom_score = score_seuil(anom_score, per_seuil, True) 

589 

590 # Normalisation lié au seuil d'anomalie (en fonction de type_norm) 

591 calib_params_ = fit_calibrate( 

592 anom_score, 

593 ctx_mask=ctx_mask, 

594 beta=beta, 

595 type_norm=type_norm, 

596 empiric_rescale=empiric_rescale, 

597 d=d, 

598 debug=debug, 

599 ) 

600 

601 return ndUQ_ratio, extremum_var_TOT, calib_params_ 

602 

603 

604def compute_anom_score( 

605 UQ, 

606 type_UQ, 

607 pred, 

608 y, 

609 type_UQ_params=None, 

610 ctx_mask=None, 

611 beta=0.005, 

612 per_seuil=0.9995, 

613 min_cut=0.01, 

614 max_cut=0.97, 

615 var_min=0, 

616 var_max=None, 

617 d=2, 

618 type_norm="quantiles_local", 

619 empiric_rescale=False, 

620 global_median_normalization=False, 

621 filt=None, 

622 reduc_filter=None, 

623 q_var=1, 

624 q_var_e=1, 

625 q_Eratio=3, 

626 k_var_e=1, 

627 with_born=False, 

628 params_=None, 

629 debug=False, 

630 roll=0, 

631 **kwargs 

632): 

633 """Compute contextual deviation score from observation, Prediction and UQ measure. 

634 Then apply normalise considering threeshold based on beta miss_covered target 

635 

636 Args: 

637 UQ (UQ-measure): UQ measure from UQ estimateurs 

638 type_UQ (str, optional): _description_. Defaults to "sigma". 

639 pred (array): prediction 

640 y (array): Real values 

641 type_UQ_params (str, optional): _description_. Defaults to "sigma". 

642 beta (float, optional): Anom target Defaults to 0.05. 

643 per_seuil (float, optional): Threshold Defaults to 0.995. 

644 min_cut (float, optional): extremum cut Defaults to 0.005. 

645 max_cut (float, optional): extremum cut Defaults to 0.995. 

646 var_min (float, optional): Threshold of minimum var, Default to 0 

647 d (int, optional): residual power Defaults to 2. 

648 type_norm (str, optional): Type of normalisation see norm Defaults to "quantiles_local". 

649 empiric_rescale (bool,optional): Force empiric rescale based on train percentile 

650 global_median_normalization (bool,optional): Normalisation of residu by each step 

651 by the median residu of all dimension. 

652 filt (list, optional): fitler applied to the score . Defaults to [0, 1, 0]. 

653 q_var (int, optional): power coefficent apply to the score. Defaults to 1. 

654 params : params provide by fit function. Defaults to None imply internal call to fit 

655 

656 Type_norm : 

657 "Nsigma_local" : Fexible anomalie thresold (N-sgima) by dimension | Hypothesis : 

658 Anomaly distributed homegeneously by dimension. 

659 "Nsigma_global" : Fexible anomalie thresold (N-sgima) global : homogeneous anomaly distribution by dimension. 

660 "Chi2": Theorical Chi2 correction 

661 "quantiles_local" : rigid anomalie thresold (quantiles) by dimension : 

662 heterogeneos anomaly distrubtion by dimension. 

663 "quantiles_global" : (Default) rigid anomalie thresold (quantiles) global : 

664 homogeneous anomaly distribution by dimension. 

665 "None" : No normalisation 

666 

667 Returns: 

668 anom_score: 2D anomaly score matrix 

669 params : parameters provided or computed 

670 """ 

671 

672 y, pred, UQ = UQ_proc.check_y_vs_pred_and_UQ_shape(y, pred, UQ) 

673 

674 # IF not fit_params : 

675 if params_ is None: 

676 params_ = fit_anom_score( 

677 UQ, 

678 type_UQ, 

679 pred, 

680 y, 

681 type_UQ_params=type_UQ_params, 

682 ctx_mask=ctx_mask, 

683 beta=beta, 

684 per_seuil=per_seuil, 

685 min_cut=min_cut, 

686 max_cut=max_cut, 

687 d=d, 

688 q_var=q_var, 

689 k_var_e=k_var_e, 

690 q_var_e=q_var_e, 

691 q_Eratio=q_Eratio, 

692 type_norm=type_norm, 

693 var_min=var_min, 

694 global_median_normalization=global_median_normalization, 

695 filt=filt, 

696 reduc_filter=reduc_filter, 

697 roll=roll, 

698 **kwargs 

699 ) 

700 

701 ndUQ_ratio, extremum_var_TOT, calib_params_ = params_ 

702 

703 # Résidu normalisé 

704 anom_score = UQ_proc.process_UQmeasure_to_residu( 

705 UQ, 

706 type_UQ, 

707 pred, 

708 y=y, 

709 type_UQ_params=type_UQ_params, 

710 d=d, 

711 min_cut=min_cut, 

712 max_cut=max_cut, 

713 q_var=q_var, 

714 var_min=var_min, 

715 k_var_e=k_var_e, 

716 q_var_e=q_var_e, 

717 q_Eratio=q_Eratio, 

718 ndUQ_ratio=ndUQ_ratio, 

719 extremum_var_TOT=extremum_var_TOT, 

720 with_born=with_born, 

721 reduc_filter=reduc_filter, 

722 roll=roll, 

723 debug=debug, 

724 ) 

725 

726 if with_born: 

727 anom_score, born = anom_score 

728 

729 global_norm = 0 

730 if global_median_normalization: 

731 global_norm = np.quantile(anom_score, 0.5, axis=-1)[:, None] 

732 anom_score = anom_score - global_norm 

733 

734 if with_born: 

735 if type_UQ == "var_A&E": 

736 sigma, E_penalisation = UQ_proc.process_UQmeasure_to_TOT_and_E_sigma( 

737 UQ, 

738 type_UQ, 

739 pred, 

740 y=y, 

741 type_UQ_params=type_UQ_params, 

742 min_cut=min_cut, 

743 max_cut=max_cut, 

744 var_min=var_min, 

745 var_max=var_max, 

746 q_var=q_var, 

747 q_var_e=q_var_e, 

748 k_var_e=k_var_e, 

749 ndUQ_ratio=ndUQ_ratio, 

750 extremum_var_TOT=extremum_var_TOT, 

751 reduc_filter=reduc_filter, 

752 roll=roll, 

753 ) 

754 else: 

755 sigma = UQ_proc.process_UQmeasure_to_sigma( 

756 UQ, 

757 type_UQ, 

758 pred, 

759 y=y, 

760 type_UQ_params=type_UQ_params, 

761 min_cut=min_cut, 

762 max_cut=max_cut, 

763 q_var=q_var, 

764 var_min=var_min, 

765 reduc_filter=reduc_filter, 

766 ) 

767 

768 E_penalisation = 0 

769 anom_padding = global_norm * (sigma) 

770 

771 else: 

772 anom_padding = 0 

773 E_penalisation = 0 

774 

775 # Reduc middle dim for 3D object (ex : multi-step prediction) 

776 if debug: 

777 print("res_val post res", np.abs(anom_score).mean()) 

778 # Apply temp Convolution 

779 

780 anom_score = apply_conv(anom_score, filt) 

781 

782 # Seuillage des résidu normalisé 

783 # -> Warning : difference between anom & born. 

784 anom_score = score_seuil(anom_score, per_seuil, True) 

785 

786 if debug: 

787 print("Start", (anom_score < -1).mean(), (anom_score > 1).mean(), beta) 

788 if with_born: 

789 print( 

790 "Start", 

791 ((born[0] + anom_padding) > y).mean(), 

792 ((born[1] + anom_padding) < y).mean(), 

793 ) 

794 

795 # IF not fit_params : 

796 

797 # Normalisation lié au seuil d'anomalie (en fonction de type_norm) 

798 

799 anom_score, _ = compute_calibrate( 

800 anom_score, 

801 ctx_mask=ctx_mask, 

802 beta=beta, 

803 type_norm=type_norm, 

804 empiric_rescale=empiric_rescale, 

805 deg=1, 

806 params_=calib_params_, 

807 debug=debug, 

808 ) 

809 

810 anom_score = np.sign(anom_score) * np.power(np.abs(anom_score), d) 

811 

812 if debug: 

813 print("End", (anom_score < -1).mean(), (anom_score > 1).mean()) 

814 

815 if with_born: 

816 pred = apply_middledim_reduction(pred, reduc_filter=reduc_filter, roll=roll) 

817 

818 res_born = np.power(sigma, d) 

819 res_born, _ = compute_calibrate( 

820 res_born, 

821 ctx_mask=ctx_mask, 

822 beta=beta, 

823 type_norm=type_norm, 

824 empiric_rescale=empiric_rescale, 

825 deg=1, 

826 params_=calib_params_, 

827 mode="born", 

828 ) 

829 res_born_bot = pred + np.minimum( 

830 0, -np.power(res_born, 1 / d) + anom_padding + E_penalisation 

831 ) 

832 

833 res_born_top = pred + np.maximum( 

834 0, np.power(res_born, 1 / d) + anom_padding - E_penalisation 

835 ) 

836 

837 born = res_born_bot, res_born_top 

838 

839 if debug: 

840 print( 

841 "Anom_born", 

842 ((y < born[0]) | (y > born[1])).mean(axis=0), 

843 "Anom_res", 

844 (np.abs(anom_score) > 1).mean(axis=0), 

845 ) 

846 return (anom_score, born), params_ 

847 else: 

848 return anom_score, params_ 

849 

850 

851def fit_score_fusion( 

852 score, 

853 ctx_mask=None, 

854 type_fusion="mahalanobis", 

855 type_norm="quantiles_global", 

856 beta=0.05, 

857 per_seuil=0.995, 

858 filt=[0, 1, 0], 

859 fusion_reduc_filter=None, 

860 fusion_debug=False, 

861 **kwargs 

862): 

863 # middle dim reduction 

864 score = UQ_proc.apply_middledim_reduction(score, fusion_reduc_filter) 

865 

866 # Temporal convolution 

867 score = score_seuil(apply_conv(score, filt), per_seuil, True) 

868 if type_fusion == "mahalanobis": 

869 # Mahalanobis pour synthese en score global. 

870 cov_estimator = EmpiricalCovariance(assume_centered=True).fit(score) 

871 anom_score = cov_estimator.mahalanobis(score) 

872 # Seuillage du score aggrégées 

873 else: 

874 cov_estimator = None 

875 anom_score = np.abs(score).mean(axis=1) 

876 

877 if fusion_debug: 

878 print("Start_fit", (anom_score < -1).mean(), (anom_score > 1).mean()) 

879 

880 anom_score = score_seuil(anom_score, per_seuil, True).reshape((-1, 1)) 

881 

882 params_calibrate_ = fit_calibrate( 

883 anom_score, ctx_mask, beta, type_norm, empiric_rescale=True, debug=fusion_debug 

884 ) 

885 

886 return (params_calibrate_, cov_estimator) 

887 

888 

889def compute_score_fusion( 

890 score, 

891 ctx_mask=None, 

892 type_fusion="mahalanobis", 

893 type_norm="quantiles_global", 

894 beta=0.05, 

895 per_seuil=0.995, 

896 filt=[0, 1, 0], 

897 d=1, 

898 fusion_reduc_filter=None, 

899 params_=None, 

900 fusion_debug=False, 

901 **kwargs 

902): 

903 """Compute and aggregate and 2D anomaly score matrix into an aggregated 1D 

904 

905 Args: 

906 score (_type_): _description_ 

907 ctx_mask (_type_): _description_ 

908 type_fusion (str, optional): mahalanobis or mean. Defaults to "mahalanobis". 

909 type_norm (str, optional): _description_. Defaults to "quantiles_global". 

910 beta (float, optional): anom ratio target. Defaults to 0.95. 

911 per_seuil (float, optional): _description_. Defaults to 0.995. 

912 filt (list, optional): _description_. Defaults to [0, 1, 0]. 

913 d (int, optional): _description_. Defaults to 2. 

914 list_norm_val (_type_, optional): _description_. Defaults to None. 

915 list_ctx_val (_type_, optional): _description_. Defaults to None. 

916 params_: params provided by fit function. Defaults to None imply interal call to fit 

917 

918 Returns: 

919 score_agg: 1D anomaly score matrix 

920 """ 

921 if params_ is None: 

922 params_ = fit_score_fusion( 

923 score, 

924 ctx_mask, 

925 type_fusion=type_fusion, 

926 type_norm=type_norm, 

927 beta=beta, 

928 per_seuil=per_seuil, 

929 filt=filt, 

930 ) 

931 

932 params_calibrate_, cov_estimator = params_ 

933 

934 # Middledim reduction 

935 score = UQ_proc.apply_middledim_reduction(score, fusion_reduc_filter) 

936 

937 # Temporal convolution 

938 score = score_seuil(apply_conv(score, filt), per_seuil, True) 

939 

940 if fusion_debug: 

941 print("Start_compute", (score < -1).mean(), (score > 1).mean()) 

942 

943 if type_fusion == "mahalanobis": 

944 # Mahalanobis pour synthese en score global. 

945 score_agg = cov_estimator.mahalanobis(score) 

946 # Seuillage du score aggrégées 

947 else: # type_fusion == "mean": 

948 score_agg = np.abs(score).mean(axis=1) 

949 

950 score_agg = score_seuil(score_agg, per_seuil, True).reshape((-1, 1))