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

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, 

272 observation and target 

273 

274 Args: 

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

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

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

278 y (np.array): Targets/Observation 

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

280 beta: target miss-coverage of the borns 

281 min_cut (_type_): Bottom extremun percentile. 

282 max_cut (_type_): Bottom upper percentile. 

283 q_var (_type_): Power coefficent 

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

285 

286 Returns: 

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

288 list_ctx_id : list of ctx_values) 

289 """ 

290 if type_res == "res": 

291 residuals = y - pred 

292 

293 if type_res == "w_res": 

294 residuals = UQ_proc.process_UQmeasure_to_residu( 

295 UQ, 

296 type_UQ, 

297 pred, 

298 y, 

299 type_UQ_params=type_UQ_params, 

300 min_cut=min_cut, 

301 max_cut=max_cut, 

302 var_min=var_min, 

303 q_var=q_var, 

304 ) 

305 

306 elif type_res == "cqr": 

307 y_pred_lower = UQ_proc.process_UQmeasure_to_quantile( 

308 UQ, 

309 type_UQ, 

310 pred, 

311 y=None, 

312 type_UQ_params=type_UQ_params, 

313 alpha=beta / 2, 

314 var_min=0, 

315 min_cut=0, 

316 max_cut=1, 

317 q_var=1, 

318 ) 

319 

320 y_pred_upper = UQ_proc.process_UQmeasure_to_quantile( 

321 UQ, 

322 type_UQ, 

323 pred, 

324 y=None, 

325 type_UQ_params=type_UQ_params, 

326 alpha=1 - beta / 2, 

327 var_min=0, 

328 min_cut=0, 

329 max_cut=1, 

330 q_var=1, 

331 ) 

332 

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

334 y_pred_lower - y 

335 ) > (y - y_pred_upper) 

336 

337 elif type_res == "no_calib": 

338 return ([1], None) 

339 

340 else: 

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

342 

343 list_norm_val, list_ctx_id = fit_calibrate( 

344 residuals, 

345 ctx_mask=ctx_mask, 

346 beta=beta, 

347 type_norm="quantiles_local", 

348 empiric_rescale=empiric_rescale, 

349 deg=1, 

350 ) 

351 

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

353 # Add symetric or asymetric parameters mode 

354 # -> to implement in calibrate_residuals 

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

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

357 # empiric_rescale=True, deg=1) 

358 params_ = (list_norm_val, list_ctx_id) 

359 return params_ 

360 

361 

362def compute_born_calibrated( 

363 UQ, 

364 type_UQ, 

365 pred, 

366 y=None, 

367 type_UQ_params=None, 

368 beta=0.1, 

369 type_res="res", 

370 ctx_mask=None, 

371 params_=None, 

372): 

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

374 

375 Args: 

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

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

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

379 y (np.array): Targets/Observation 

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

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

382 type_res : 

383 "no_calib : No calibration 

384 "res" : Calibration based on mean residuals 

385 "w_res" : Calibration based on weigthed mean residuals 

386 "cqr" : Calibration based on quantile residuals 

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

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

389 

390 Raises: 

391 ValueError: _description_ 

392 

393 Returns: 

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

395 """ 

396 if ctx_mask is None: 

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

398 

399 if params_ is None: 

400 params_ = fit_born_calibrated( 

401 UQ, 

402 type_UQ, 

403 pred, 

404 y, 

405 type_UQ_params, 

406 beta=beta, 

407 type_res=type_res, 

408 ctx_mask=ctx_mask, 

409 ) 

410 

411 list_norm_val, list_ctx_id = params_ 

412 

413 if type_res == "res": 

414 UQ = pred * 0 + 1 

415 type_UQ = "var" 

416 

417 y_pred_lower = UQ_proc.process_UQmeasure_to_quantile( 

418 UQ, 

419 type_UQ, 

420 pred, 

421 y=None, 

422 type_UQ_params=type_UQ_params, 

423 alpha=beta / 2, 

424 var_min=0, 

425 min_cut=0, 

426 max_cut=1, 

427 q_var=1, 

428 ) 

429 

430 y_pred_upper = UQ_proc.process_UQmeasure_to_quantile( 

431 UQ, 

432 type_UQ, 

433 pred, 

434 y=None, 

435 type_UQ_params=type_UQ_params, 

436 alpha=1 - beta / 2, 

437 var_min=0, 

438 min_cut=0, 

439 max_cut=1, 

440 q_var=1, 

441 ) 

442 

443 for n, ctx_val in list_ctx_id: 

444 mask = ctx_mask == ctx_val 

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

446 y_pred_lower[mask] = ( 

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

448 ) + pred[mask] 

449 y_pred_upper[mask] = ( 

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

451 ) + pred[mask] 

452 

453 elif type_res in ["cqr"]: 

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

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

456 else: 

457 raise ValueError( 

458 "compute_UQ_calibration : type_res : " + type_res + " Not covered" 

459 ) 

460 # Hold assymetrical val 

461 

462 return (y_pred_lower, y_pred_upper), params_ 

463 

464 

465# FONCTION PRINCIPALE : CONSTRUCTION DU SCORE. 

466# res = Résidu source 

467# v_mean = biais source 

468# v_std = variance source. 

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

470 

471# Paramètre : 

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

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

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

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

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

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

478# Type_Norm : 

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

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

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

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

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

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

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

486# Dim = dimension de la série. 

487 

488 

489def fit_anom_score( 

490 UQ, 

491 type_UQ, 

492 pred, 

493 y, 

494 type_UQ_params=None, 

495 ctx_mask=None, 

496 beta=0.01, 

497 per_seuil=0.9995, 

498 min_cut=0.01, 

499 max_cut=0.97, 

500 var_min=0, 

501 d=2, 

502 q_var=1, 

503 k_var_e=1, 

504 q_var_e=1, 

505 q_Eratio=3, 

506 type_norm="quantiles_local", 

507 empiric_rescale=False, 

508 global_median_normalization=False, 

509 filt=None, 

510 reduc_filter=None, 

511 roll=0, 

512 debug=False, 

513 **kwargs 

514): 

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

516 

517 Args: 

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

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

520 pred (array): prediction 

521 y (array): Real values 

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

523 ctx_mask : mask for ctx_calib 

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

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

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

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

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

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

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

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

532 by the median residu of all dimension. 

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

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

535 

536 Returns: 

537 params : parameters to calibrate anomaly score 

538 """ 

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

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

541 

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

543 

544 ndUQ_ratio, extremum_var_TOT = None, (None, None) 

545 if type_UQ == "var_A&E": 

546 extremum_var_TOT, ndUQ_ratio = UQ_proc.get_extremum_var_TOT_and_ndUQ_ratio( 

547 UQ, 

548 type_UQ=type_UQ, 

549 min_cut=min_cut, 

550 max_cut=max_cut, 

551 var_min=0, 

552 var_max=None, 

553 factor=2, 

554 q_var=q_var, 

555 q_Eratio=q_Eratio, 

556 mode_multidim=True, 

557 ) 

558 

559 anom_score = UQ_proc.process_UQmeasure_to_residu( 

560 UQ, 

561 type_UQ, 

562 pred, 

563 y=y, 

564 type_UQ_params=type_UQ_params, 

565 d=d, 

566 min_cut=min_cut, 

567 max_cut=max_cut, 

568 q_var=q_var, 

569 var_min=var_min, 

570 k_var_e=k_var_e, 

571 q_var_e=q_var_e, 

572 q_Eratio=q_Eratio, 

573 ndUQ_ratio=ndUQ_ratio, 

574 extremum_var_TOT=extremum_var_TOT, 

575 reduc_filter=reduc_filter, 

576 roll=roll, 

577 debug=debug, 

578 ) 

579 

580 if global_median_normalization: 

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

582 anom_score = anom_score - global_norm 

583 

584 # Dim reduction : 

585 

586 # Seuillage des résidu normalisé 

587 

588 anom_score = apply_conv(anom_score, filt) 

589 anom_score = score_seuil(anom_score, per_seuil, True) 

590 

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

592 calib_params_ = fit_calibrate( 

593 anom_score, 

594 ctx_mask=ctx_mask, 

595 beta=beta, 

596 type_norm=type_norm, 

597 empiric_rescale=empiric_rescale, 

598 d=d, 

599 debug=debug, 

600 ) 

601 

602 return ndUQ_ratio, extremum_var_TOT, calib_params_ 

603 

604 

605def compute_anom_score( 

606 UQ, 

607 type_UQ, 

608 pred, 

609 y, 

610 type_UQ_params=None, 

611 ctx_mask=None, 

612 beta=0.005, 

613 per_seuil=0.9995, 

614 min_cut=0.01, 

615 max_cut=0.97, 

616 var_min=0, 

617 var_max=None, 

618 d=2, 

619 type_norm="quantiles_local", 

620 empiric_rescale=False, 

621 global_median_normalization=False, 

622 filt=None, 

623 reduc_filter=None, 

624 q_var=1, 

625 q_var_e=1, 

626 q_Eratio=3, 

627 k_var_e=1, 

628 with_born=False, 

629 params_=None, 

630 debug=False, 

631 roll=0, 

632 **kwargs 

633): 

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

635 Then apply normalise considering threeshold based on beta miss_covered target 

636 

637 Args: 

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

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

640 pred (array): prediction 

641 y (array): Real values 

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

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

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

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

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

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

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

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

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

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

652 by the median residu of all dimension. 

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

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

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

656 

657 Type_norm : 

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

659 Anomaly distributed homegeneously by dimension. 

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

661 "Chi2": Theorical Chi2 correction 

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

663 heterogeneos anomaly distrubtion by dimension. 

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

665 homogeneous anomaly distribution by dimension. 

666 "None" : No normalisation 

667 

668 Returns: 

669 anom_score: 2D anomaly score matrix 

670 params : parameters provided or computed 

671 """ 

672 

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

674 

675 # IF not fit_params : 

676 if params_ is None: 

677 params_ = fit_anom_score( 

678 UQ, 

679 type_UQ, 

680 pred, 

681 y, 

682 type_UQ_params=type_UQ_params, 

683 ctx_mask=ctx_mask, 

684 beta=beta, 

685 per_seuil=per_seuil, 

686 min_cut=min_cut, 

687 max_cut=max_cut, 

688 d=d, 

689 q_var=q_var, 

690 k_var_e=k_var_e, 

691 q_var_e=q_var_e, 

692 q_Eratio=q_Eratio, 

693 type_norm=type_norm, 

694 var_min=var_min, 

695 global_median_normalization=global_median_normalization, 

696 filt=filt, 

697 reduc_filter=reduc_filter, 

698 roll=roll, 

699 **kwargs 

700 ) 

701 

702 ndUQ_ratio, extremum_var_TOT, calib_params_ = params_ 

703 

704 # Résidu normalisé 

705 anom_score = UQ_proc.process_UQmeasure_to_residu( 

706 UQ, 

707 type_UQ, 

708 pred, 

709 y=y, 

710 type_UQ_params=type_UQ_params, 

711 d=d, 

712 min_cut=min_cut, 

713 max_cut=max_cut, 

714 q_var=q_var, 

715 var_min=var_min, 

716 k_var_e=k_var_e, 

717 q_var_e=q_var_e, 

718 q_Eratio=q_Eratio, 

719 ndUQ_ratio=ndUQ_ratio, 

720 extremum_var_TOT=extremum_var_TOT, 

721 with_born=with_born, 

722 reduc_filter=reduc_filter, 

723 roll=roll, 

724 debug=debug, 

725 ) 

726 

727 if with_born: 

728 anom_score, born = anom_score 

729 

730 global_norm = 0 

731 if global_median_normalization: 

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

733 anom_score = anom_score - global_norm 

734 

735 if with_born: 

736 if type_UQ == "var_A&E": 

737 sigma, E_penalisation = UQ_proc.process_UQmeasure_to_TOT_and_E_sigma( 

738 UQ, 

739 type_UQ, 

740 pred, 

741 y=y, 

742 type_UQ_params=type_UQ_params, 

743 min_cut=min_cut, 

744 max_cut=max_cut, 

745 var_min=var_min, 

746 var_max=var_max, 

747 q_var=q_var, 

748 q_var_e=q_var_e, 

749 k_var_e=k_var_e, 

750 ndUQ_ratio=ndUQ_ratio, 

751 extremum_var_TOT=extremum_var_TOT, 

752 reduc_filter=reduc_filter, 

753 roll=roll, 

754 ) 

755 else: 

756 sigma = UQ_proc.process_UQmeasure_to_sigma( 

757 UQ, 

758 type_UQ, 

759 pred, 

760 y=y, 

761 type_UQ_params=type_UQ_params, 

762 min_cut=min_cut, 

763 max_cut=max_cut, 

764 q_var=q_var, 

765 var_min=var_min, 

766 reduc_filter=reduc_filter, 

767 ) 

768 

769 E_penalisation = 0 

770 anom_padding = global_norm * (sigma) 

771 

772 else: 

773 anom_padding = 0 

774 E_penalisation = 0 

775 

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

777 if debug: 

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

779 # Apply temp Convolution 

780 

781 anom_score = apply_conv(anom_score, filt) 

782 

783 # Seuillage des résidu normalisé 

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

785 anom_score = score_seuil(anom_score, per_seuil, True) 

786 

787 if debug: 

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

789 if with_born: 

790 print( 

791 "Start", 

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

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

794 ) 

795 

796 # IF not fit_params : 

797 

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

799 

800 anom_score, _ = compute_calibrate( 

801 anom_score, 

802 ctx_mask=ctx_mask, 

803 beta=beta, 

804 type_norm=type_norm, 

805 empiric_rescale=empiric_rescale, 

806 deg=1, 

807 params_=calib_params_, 

808 debug=debug, 

809 ) 

810 

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

812 

813 if debug: 

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

815 

816 if with_born: 

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

818 

819 res_born = np.power(sigma, d) 

820 res_born, _ = compute_calibrate( 

821 res_born, 

822 ctx_mask=ctx_mask, 

823 beta=beta, 

824 type_norm=type_norm, 

825 empiric_rescale=empiric_rescale, 

826 deg=1, 

827 params_=calib_params_, 

828 mode="born", 

829 ) 

830 res_born_bot = pred + np.minimum( 

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

832 ) 

833 

834 res_born_top = pred + np.maximum( 

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

836 ) 

837 

838 born = res_born_bot, res_born_top 

839 

840 if debug: 

841 print( 

842 "Anom_born", 

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

844 "Anom_res", 

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

846 ) 

847 return (anom_score, born), params_ 

848 else: 

849 return anom_score, params_ 

850 

851 

852def fit_score_fusion( 

853 score, 

854 ctx_mask=None, 

855 type_fusion="mahalanobis", 

856 type_norm="quantiles_global", 

857 beta=0.05, 

858 per_seuil=0.995, 

859 filt=[0, 1, 0], 

860 fusion_reduc_filter=None, 

861 fusion_debug=False, 

862 **kwargs 

863): 

864 # middle dim reduction 

865 score = UQ_proc.apply_middledim_reduction(score, fusion_reduc_filter) 

866 

867 # Temporal convolution 

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

869 if type_fusion == "mahalanobis": 

870 # Mahalanobis pour synthese en score global. 

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

872 anom_score = cov_estimator.mahalanobis(score) 

873 # Seuillage du score aggrégées 

874 else: 

875 cov_estimator = None 

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

877 

878 if fusion_debug: 

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

880 

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

882 

883 params_calibrate_ = fit_calibrate( 

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

885 ) 

886 

887 return (params_calibrate_, cov_estimator) 

888 

889 

890def compute_score_fusion( 

891 score, 

892 ctx_mask=None, 

893 type_fusion="mahalanobis", 

894 type_norm="quantiles_global", 

895 beta=0.05, 

896 per_seuil=0.995, 

897 filt=[0, 1, 0], 

898 d=1, 

899 fusion_reduc_filter=None, 

900 params_=None, 

901 fusion_debug=False, 

902 **kwargs 

903): 

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

905 

906 Args: 

907 score (_type_): _description_ 

908 ctx_mask (_type_): _description_ 

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

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

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

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

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

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

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

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

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

918 

919 Returns: 

920 score_agg: 1D anomaly score matrix 

921 """ 

922 if params_ is None: 

923 params_ = fit_score_fusion( 

924 score, 

925 ctx_mask, 

926 type_fusion=type_fusion, 

927 type_norm=type_norm, 

928 beta=beta, 

929 per_seuil=per_seuil, 

930 filt=filt, 

931 ) 

932 

933 params_calibrate_, cov_estimator = params_ 

934 

935 # Middledim reduction 

936 score = UQ_proc.apply_middledim_reduction(score, fusion_reduc_filter) 

937 

938 # Temporal convolution 

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

940 

941 if fusion_debug: 

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

943 

944 if type_fusion == "mahalanobis": 

945 # Mahalanobis pour synthese en score global. 

946 score_agg = cov_estimator.mahalanobis(score) 

947 # Seuillage du score aggrégées 

948 else: # type_fusion == "mean": 

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

950