smashbox.stats.stats

   1#!/usr/bin/env python3
   2# -*- coding: utf-8 -*-
   3"""
   4Created on Tue Jul 22 10:47:46 2025
   5
   6@author: maxime
   7"""
   8
   9import numpy as np
  10import pandas as pd
  11import warnings
  12import scipy.stats as stats
  13import multiprocessing
  14import matplotlib.pyplot as plt
  15import os
  16from tqdm import tqdm
  17import threading
  18from smashbox.tools import tools
  19from smash.fcore import _mwd_metrics as smash_metrics
  20
  21
  22@tools.autocast_args
  23def _get_mask_colum_nodata(arr, nodata=-99.0, t_axis=0):
  24    """
  25    Test if all data for axis 1 of an array with 2 dimension are nodata
  26    The function return an array with shape equal of the shape of axis 1 of the
  27    orignal array `arr`. This returned array is zeros everywhere but np.nan if all
  28    values of arr are nodata for axis 1.
  29
  30    parameter:
  31    ----------
  32    arr: np.array
  33    nodata: float | int | np.nan
  34    t_axis: array axis for time
  35    """
  36
  37    arr = np.moveaxis(arr, t_axis, 0)
  38    c_axis = 1
  39
  40    res_nodata = None
  41    if len(arr.shape) == 2:
  42        res_nodata = np.zeros(shape=(arr.shape[c_axis]))
  43        for i in range(arr.shape[c_axis]):
  44            if np.all(arr[:, i] <= nodata):
  45                res_nodata[i] = np.nan
  46
  47    # print(res_nodata)
  48
  49    return res_nodata
  50
  51
  52def _test_input_shape(arr):
  53    """
  54    Test the shape of the array `arr`
  55    :param arr: np.ndarray
  56    :type arr: np.ndarray
  57
  58    """
  59
  60    if len(arr.shape) > 2:
  61        raise ValueError(
  62            "Dimension of the input array must not be greater than 2:"
  63            f"given shape is {arr.shape} with {len(arr.shape)} dimension."
  64        )
  65
  66
  67@tools.autocast_args
  68def mse(
  69    obs: np.ndarray | None = None,
  70    sim: np.ndarray | None = None,
  71    nodata: float = -99.0,
  72    t_axis: int = 0,
  73):
  74    """
  75    Compute the misfit criteria mse:
  76        mse = (1.0 / nb_valid_data) * np.sum((obs - sim) ** 2.0)
  77    :param obs: Observed discharges, defaults to None
  78    :type obs: np.ndarray | None, optional
  79    :param sim: simulated discharges, defaults to None
  80    :type sim: np.ndarray | None, optional
  81    :param nodata: No data value, defaults to -99.0
  82    :type nodata: float, optional
  83    :param t_axis: Array axis of the time, defaults to 0
  84    :type t_axis: int, optional
  85
  86    """
  87
  88    t_axis = min(t_axis, len(obs.shape) - 1)
  89    _test_input_shape(obs)
  90    _test_input_shape(sim)
  91
  92    mask_nodata = obs != nodata
  93    nb_valid_data = np.count_nonzero(mask_nodata, axis=t_axis)
  94    res_nodata = _get_mask_colum_nodata(obs, nodata=nodata, t_axis=t_axis)
  95    mse = np.nan
  96
  97    if sum(nb_valid_data) > 0:
  98
  99        if isinstance(obs, np.ndarray) and isinstance(sim, np.ndarray):
 100
 101            if obs.shape == sim.shape:
 102                mse = (1.0 / nb_valid_data) * np.sum(
 103                    (obs - sim) ** 2.0, axis=t_axis, where=mask_nodata
 104                )
 105            else:
 106                raise ValueError(f"Error: len(obs)!=len(sim), {len(obs)}!={len(sim)}")
 107        else:
 108            raise ValueError("Error: obs and sim must be an instance of np.ndarray")
 109
 110    else:
 111        raise ValueError("Error: no valid observation data.")
 112
 113    if res_nodata is not None:
 114        mse = mse + res_nodata
 115
 116    return mse
 117
 118
 119@tools.autocast_args
 120def sm_mse(
 121    obs: np.ndarray | None = None,
 122    sim: np.ndarray | None = None,
 123):
 124    """
 125    Compute the misfit criteria mse:
 126        mse = (1.0 / nb_valid_data) * np.sum((obs - sim) ** 2.0)
 127    :param obs: Observed discharges, defaults to None
 128    :type obs: np.ndarray | None, optional
 129    :param sim: simulated discharges, defaults to None
 130    :type sim: np.ndarray | None, optional
 131    """
 132
 133    mse = smash_metrics.mse(obs, sim)
 134
 135    return mse
 136
 137
 138@tools.autocast_args
 139def rmse(
 140    obs: np.ndarray | None = None,
 141    sim: np.ndarray | None = None,
 142    nodata: float = -99.0,
 143    t_axis: int = 0,
 144):
 145    """
 146    Compute the misfit criteria rmse:
 147        rmse = np.sqrt(res_mse)
 148
 149    :param obs: Observed discharges, defaults to None
 150    :type obs: np.ndarray | None, optional
 151    :param sim: simulated discharges, defaults to None
 152    :type sim: np.ndarray | None, optional
 153    :param nodata: No data value, defaults to -99.0
 154    :type nodata: float, optional
 155    :param t_axis: Array axis of the time, defaults to 0
 156    :type t_axis: int, optional
 157
 158    """
 159    t_axis = min(t_axis, len(obs.shape) - 1)
 160    res_mse = mse(obs, sim, nodata=nodata, t_axis=t_axis)
 161
 162    rmse = np.sqrt(res_mse)
 163
 164    return rmse
 165
 166
 167@tools.autocast_args
 168def sm_rmse(
 169    obs: np.ndarray | None = None,
 170    sim: np.ndarray | None = None,
 171):
 172    """
 173    Compute the misfit criteria rmse:
 174        rmse = np.sqrt(res_mse)
 175
 176    :param obs: Observed discharges, defaults to None
 177    :type obs: np.ndarray | None, optional
 178    :param sim: simulated discharges, defaults to None
 179    :type sim: np.ndarray | None, optional
 180
 181    """
 182
 183    rmse = smash_metrics.rmse(obs, sim)
 184
 185    return rmse
 186
 187
 188@tools.autocast_args
 189def nrmse(
 190    obs: np.ndarray | None = None,
 191    sim: np.ndarray | None = None,
 192    nodata: float = -99.0,
 193    t_axis: int = 0,
 194):
 195    """
 196    Compute the misfit criteria nrmse:
 197        nrmse = res_rmse / mean_obs
 198
 199    :param obs: Observed discharges, defaults to None
 200    :type obs: np.ndarray | None, optional
 201    :param sim: simulated discharges, defaults to None
 202    :type sim: np.ndarray | None, optional
 203    :param nodata: No data value, defaults to -99.0
 204    :type nodata: float, optional
 205    :param t_axis: Array axis of the time, defaults to 0
 206    :type t_axis: int, optional
 207
 208    """
 209
 210    t_axis = min(t_axis, len(obs.shape) - 1)
 211    res_rmse = rmse(obs, sim, nodata=nodata, t_axis=t_axis)
 212    mask_nodata = obs != nodata
 213
 214    with warnings.catch_warnings():
 215        warnings.simplefilter("ignore", category=RuntimeWarning)
 216        mean_obs = np.mean(obs, axis=t_axis, where=mask_nodata)
 217
 218    nrmse = res_rmse / mean_obs
 219
 220    return nrmse
 221
 222
 223@tools.autocast_args
 224def sm_nmse(
 225    obs: np.ndarray | None = None,
 226    sim: np.ndarray | None = None,
 227):
 228    """
 229    Compute the misfit criteria nrmse:
 230        nrmse = res_rmse / mean_obs
 231
 232    :param obs: Observed discharges, defaults to None
 233    :type obs: np.ndarray | None, optional
 234    :param sim: simulated discharges, defaults to None
 235    :type sim: np.ndarray | None, optional
 236
 237    """
 238
 239    mask_nodata = obs < 0.0
 240    mean_obs = np.nanmean(obs, where=mask_nodata)
 241    mse = smash_metrics.mse(obs, sim)
 242
 243    return mse
 244
 245
 246@tools.autocast_args
 247def se(
 248    obs: np.ndarray | None = None,
 249    sim: np.ndarray | None = None,
 250    nodata: float = -99.0,
 251    t_axis: int = 0,
 252):
 253    """
 254    Compute the misfit criteria se:
 255        se = (
 256            np.sum((obs - sim)** 2.0, axis=t_axis, where=mask_nodata)
 257        )
 258
 259    :param obs: Observed discharges, defaults to None
 260    :type obs: np.ndarray | None, optional
 261    :param sim: simulated discharges, defaults to None
 262    :type sim: np.ndarray | None, optional
 263    :param nodata: No data value, defaults to -99.0
 264    :type nodata: float, optional
 265    :param t_axis: Array axis of the time, defaults to 0
 266    :type t_axis: int, optional
 267
 268    """
 269
 270    t_axis = min(t_axis, len(obs.shape) - 1)
 271    _test_input_shape(obs)
 272    _test_input_shape(sim)
 273
 274    mask_nodata = obs != nodata
 275    nb_valid_data = np.count_nonzero(mask_nodata, axis=t_axis)
 276    res_nodata = _get_mask_colum_nodata(obs, nodata=nodata, t_axis=t_axis)
 277    se = np.nan
 278
 279    if sum(nb_valid_data) > 0:
 280
 281        if isinstance(obs, np.ndarray) and isinstance(sim, np.ndarray):
 282
 283            if obs.shape == sim.shape:
 284                se = np.sum((obs - sim) ** 2.0, axis=t_axis, where=mask_nodata)
 285            else:
 286                raise ValueError(f"Error: len(obs)!=len(sim), {len(obs)}!={len(sim)}")
 287        else:
 288            raise ValueError("Error: obs and sim must be an instance of np.ndarray")
 289
 290    else:
 291        raise ValueError("Error: no valid observation data.")
 292
 293    if res_nodata is not None:
 294        se = se + res_nodata
 295
 296    return se
 297
 298
 299@tools.autocast_args
 300def sm_se(
 301    obs: np.ndarray | None = None,
 302    sim: np.ndarray | None = None,
 303):
 304    """
 305    Compute the misfit criteria se:
 306        se = (
 307            np.sum((obs - sim)** 2.0, axis=t_axis, where=mask_nodata)
 308        )
 309
 310    :param obs: Observed discharges, defaults to None
 311    :type obs: np.ndarray | None, optional
 312    :param sim: simulated discharges, defaults to None
 313    :type sim: np.ndarray | None, optional
 314    """
 315
 316    se = smash_metrics.se(obs, sim)
 317
 318    return se
 319
 320
 321@tools.autocast_args
 322def mae(
 323    obs: np.ndarray | None = None,
 324    sim: np.ndarray | None = None,
 325    nodata: float = -99.0,
 326    t_axis: int = 0,
 327):
 328    """
 329    Compute the misfit criteria mae:
 330        mae = 1/n*(np.sum(abs(obs - sim))
 331
 332    :param obs: Observed discharges, defaults to None
 333    :type obs: np.ndarray | None, optional
 334    :param sim: simulated discharges, defaults to None
 335    :type sim: np.ndarray | None, optional
 336    :param nodata: No data value, defaults to -99.0
 337    :type nodata: float, optional
 338    :param t_axis: Array axis of the time, defaults to 0
 339    :type t_axis: int, optional
 340
 341    """
 342
 343    t_axis = min(t_axis, len(obs.shape) - 1)
 344    _test_input_shape(obs)
 345    _test_input_shape(sim)
 346
 347    mask_nodata = obs != nodata
 348    nb_valid_data = np.count_nonzero(mask_nodata, axis=t_axis)
 349    res_nodata = _get_mask_colum_nodata(obs, nodata=nodata, t_axis=t_axis)
 350    mae = np.nan
 351
 352    if sum(nb_valid_data) > 0:
 353
 354        if isinstance(obs, np.ndarray) and isinstance(sim, np.ndarray):
 355
 356            if obs.shape == sim.shape:
 357                mae = (1.0 / nb_valid_data) * (
 358                    np.sum(abs(obs - sim), axis=t_axis, where=mask_nodata)
 359                )
 360            else:
 361                raise ValueError(f"Error: len(obs)!=len(sim), {len(obs)}!={len(sim)}")
 362        else:
 363            raise ValueError("Error: obs and sim must be an instance of np.ndarray")
 364
 365    else:
 366        raise ValueError("Error: no valid observation data.")
 367
 368    if res_nodata is not None:
 369        mae = mae + res_nodata
 370
 371    return mae
 372
 373
 374@tools.autocast_args
 375def sm_mae(
 376    obs: np.ndarray | None = None,
 377    sim: np.ndarray | None = None,
 378):
 379    """
 380    Compute the misfit criteria mae:
 381        mae = np.sqrt(np.sum(abs(obs - sim))
 382
 383    :param obs: Observed discharges, defaults to None
 384    :type obs: np.ndarray | None, optional
 385    :param sim: simulated discharges, defaults to None
 386    :type sim: np.ndarray | None, optional
 387    :param nodata: No data value, defaults to -99.0
 388    :type nodata: float, optional
 389    :param t_axis: Array axis of the time, defaults to 0
 390    :type t_axis: int, optional
 391
 392    """
 393
 394    mae = smash_metrics.mae(obs, sim)
 395
 396    return mae
 397
 398
 399@tools.autocast_args
 400def mape(
 401    obs: np.ndarray | None = None,
 402    sim: np.ndarray | None = None,
 403    nodata: float = -99.0,
 404    t_axis: int = 0,
 405):
 406    """
 407    Compute the misfit criteria mape:
 408        mape = 1/n*(
 409            np.sum(abs((obs - sim) / obs))
 410        )
 411
 412    :param obs: Observed discharges, defaults to None
 413    :type obs: np.ndarray | None, optional
 414    :param sim: simulated discharges, defaults to None
 415    :type sim: np.ndarray | None, optional
 416    :param nodata: No data value, defaults to -99.0
 417    :type nodata: float, optional
 418    :param t_axis: Array axis of the time, defaults to 0
 419    :type t_axis: int, optional
 420
 421    """
 422
 423    t_axis = min(t_axis, len(obs.shape) - 1)
 424    _test_input_shape(obs)
 425    _test_input_shape(sim)
 426
 427    mask_nodata = obs != nodata
 428    nb_valid_data = np.count_nonzero(mask_nodata, axis=t_axis)
 429    res_nodata = _get_mask_colum_nodata(obs, nodata=nodata, t_axis=t_axis)
 430    mape = np.nan
 431
 432    if sum(nb_valid_data) > 0:
 433
 434        if isinstance(obs, np.ndarray) and isinstance(sim, np.ndarray):
 435
 436            if obs.shape == sim.shape:
 437                mape = (1.0 / nb_valid_data) * (
 438                    np.sum(abs((obs - sim) / obs), axis=t_axis, where=mask_nodata)
 439                )
 440            else:
 441                raise ValueError(f"Error: len(obs)!=len(sim), {len(obs)}!={len(sim)}")
 442        else:
 443            raise ValueError("Error: obs and sim must be an instance of np.ndarray")
 444
 445    else:
 446        raise ValueError("Error: no valid observation data.")
 447
 448    if res_nodata is not None:
 449        mape = mape + res_nodata
 450
 451    return mape
 452
 453
 454@tools.autocast_args
 455def sm_mape(
 456    obs: np.ndarray | None = None,
 457    sim: np.ndarray | None = None,
 458):
 459    """
 460    Compute the misfit criteria mape:
 461        mape = np.sqrt(
 462            np.sum(abs((obs - sim) / obs))
 463        )
 464
 465    :param obs: Observed discharges, defaults to None
 466    :type obs: np.ndarray | None, optional
 467    :param sim: simulated discharges, defaults to None
 468    :type sim: np.ndarray | None, optional
 469
 470    """
 471
 472    mape = smash_metrics.mape(obs, sim)
 473
 474    return mape
 475
 476
 477@tools.autocast_args
 478def lgrm(
 479    obs: np.ndarray | None = None,
 480    sim: np.ndarray | None = None,
 481    nodata: float = -99.0,
 482    t_axis: int = 0,
 483):
 484    """
 485    Compute the misfit criteria lgrm:
 486        lgrm = np.sum(
 487            obs * (np.log((obs / sim) ** 2.0)), axis=t_axis, where=mask_nodata
 488        )
 489
 490    :param obs: Observed discharges, defaults to None
 491    :type obs: np.ndarray | None, optional
 492    :param sim: simulated discharges, defaults to None
 493    :type sim: np.ndarray | None, optional
 494    :param nodata: No data value, defaults to -99.0
 495    :type nodata: float, optional
 496    :param t_axis: Array axis of the time, defaults to 0
 497    :type t_axis: int, optional
 498
 499    """
 500
 501    t_axis = min(t_axis, len(obs.shape) - 1)
 502    _test_input_shape(obs)
 503    _test_input_shape(sim)
 504
 505    mask_nodata = obs != nodata
 506    nb_valid_data = np.count_nonzero(mask_nodata, axis=t_axis)
 507    res_nodata = _get_mask_colum_nodata(obs, nodata=nodata, t_axis=t_axis)
 508    lgrm = np.nan
 509
 510    if sum(nb_valid_data) > 0:
 511
 512        if isinstance(obs, np.ndarray) and isinstance(sim, np.ndarray):
 513
 514            if obs.shape == sim.shape:
 515                lgrm = np.sum(
 516                    obs * (np.log(sim / obs) ** 2.0), axis=t_axis, where=mask_nodata
 517                )
 518
 519            else:
 520                raise ValueError(f"Error: len(obs)!=len(sim), {len(obs)}!={len(sim)}")
 521        else:
 522            raise ValueError("Error: obs and sim must be an instance of np.ndarray")
 523
 524    else:
 525        raise ValueError("Error: no valid observation data.")
 526
 527    if res_nodata is not None:
 528        lgrm = lgrm + res_nodata
 529
 530    return lgrm
 531
 532
 533@tools.autocast_args
 534def sm_lgrm(
 535    obs: np.ndarray | None = None,
 536    sim: np.ndarray | None = None,
 537):
 538    """
 539    Compute the misfit criteria lgrm:
 540        lgrm = np.sum(
 541            obs * (np.log((obs / sim) ** 2.0)), axis=t_axis, where=mask_nodata
 542        )
 543
 544    :param obs: Observed discharges, defaults to None
 545    :type obs: np.ndarray | None, optional
 546    :param sim: simulated discharges, defaults to None
 547    :type sim: np.ndarray | None, optional
 548    """
 549
 550    lgrm = smash_metrics.lgrm(obs, sim)
 551
 552    return lgrm
 553
 554
 555@tools.autocast_args
 556def nse(
 557    obs: np.ndarray | None = None,
 558    sim: np.ndarray | None = None,
 559    nodata: float = -99.0,
 560    t_axis: int = 0,
 561):
 562    """
 563    Compute the misfit criteria nse:
 564        numerator = np.sum((obs - sim) ** 2.0, axis=t_axis, where=mask_nodata)
 565
 566        denominator = np.sum(
 567            (obs - mean_obs) ** 2.0, axis=t_axis, where=mask_nodata
 568        )
 569
 570        denominator = np.where(denominator == 0, np.nan, denominator)
 571
 572        nse = 1 - numerator / denominator
 573
 574    :param obs: Observed discharges, defaults to None
 575    :type obs: np.ndarray | None, optional
 576    :param sim: simulated discharges, defaults to None
 577    :type sim: np.ndarray | None, optional
 578    :param nodata: No data value, defaults to -99.0
 579    :type nodata: float, optional
 580    :param t_axis: Array axis of the time, defaults to 0
 581    :type t_axis: int, optional
 582
 583    """
 584
 585    t_axis = min(t_axis, len(obs.shape) - 1)
 586    _test_input_shape(obs)
 587    _test_input_shape(sim)
 588
 589    mask_nodata = obs != nodata
 590    nb_valid_data = np.count_nonzero(mask_nodata, axis=t_axis)
 591
 592    nse = np.nan
 593
 594    if sum(nb_valid_data) > 0:
 595
 596        if isinstance(obs, np.ndarray) and isinstance(sim, np.ndarray):
 597            if obs.shape == sim.shape:
 598
 599                with warnings.catch_warnings():
 600                    warnings.simplefilter("ignore", category=RuntimeWarning)
 601                    mean_obs = np.mean(obs, axis=t_axis, where=mask_nodata, keepdims=True)
 602
 603                numerator = np.sum((obs - sim) ** 2.0, axis=t_axis, where=mask_nodata)
 604
 605                denominator = np.sum(
 606                    (obs - mean_obs) ** 2.0, axis=t_axis, where=mask_nodata
 607                )
 608
 609                denominator = np.where(denominator == 0, np.nan, denominator)
 610                nse = 1 - numerator / denominator
 611            else:
 612                raise ValueError(f"Error: len(obs)!=len(sim), {len(obs)}!={len(sim)}")
 613        else:
 614            raise ValueError("Error: obs and sim must be an instance of np.ndarray")
 615
 616    else:
 617        raise ValueError("Error: no valid observation data.")
 618
 619    return nse
 620
 621
 622@tools.autocast_args
 623def sm_nse(
 624    obs: np.ndarray | None = None,
 625    sim: np.ndarray | None = None,
 626):
 627    """
 628    Compute the misfit criteria nse:
 629        numerator = np.sum((obs - sim) ** 2.0, axis=t_axis, where=mask_nodata)
 630
 631        denominator = np.sum(
 632            (obs - mean_obs) ** 2.0, axis=t_axis, where=mask_nodata
 633        )
 634
 635        denominator = np.where(denominator == 0, np.nan, denominator)
 636
 637        nse = 1 - numerator / denominator
 638
 639    :param obs: Observed discharges, defaults to None
 640    :type obs: np.ndarray | None, optional
 641    :param sim: simulated discharges, defaults to None
 642    :type sim: np.ndarray | None, optional
 643
 644    """
 645
 646    nse = smash_metrics.nse(obs, sim)
 647
 648    return nse
 649
 650
 651@tools.autocast_args
 652def nnse(
 653    obs: np.ndarray | None = None,
 654    sim: np.ndarray | None = None,
 655    nodata: float = -99.0,
 656    t_axis: int = 0,
 657):
 658    """
 659    Compute the misfit criteria nnse:
 660        nnse = 1.0 / (2.0 - res_nse)
 661
 662    :param obs: Observed discharges, defaults to None
 663    :type obs: np.ndarray | None, optional
 664    :param sim: simulated discharges, defaults to None
 665    :type sim: np.ndarray | None, optional
 666    :param nodata: No data value, defaults to -99.0
 667    :type nodata: float, optional
 668    :param t_axis: Array axis of the time, defaults to 0
 669    :type t_axis: int, optional
 670
 671    """
 672
 673    t_axis = min(t_axis, len(obs.shape) - 1)
 674    res_nse = nse(obs, sim, nodata=nodata, t_axis=t_axis)
 675    nnse = 1.0 / (2.0 - res_nse)
 676
 677    return nnse
 678
 679
 680@tools.autocast_args
 681def sm_nnse(
 682    obs: np.ndarray | None = None,
 683    sim: np.ndarray | None = None,
 684):
 685    """
 686    Compute the misfit criteria nnse:
 687        nnse = 1.0 / (2.0 - res_nse)
 688
 689    :param obs: Observed discharges, defaults to None
 690    :type obs: np.ndarray | None, optional
 691    :param sim: simulated discharges, defaults to None
 692    :type sim: np.ndarray | None, optional
 693
 694    """
 695
 696    nnse = smash_metrics.nnse(obs, sim)
 697
 698    return nse
 699
 700
 701@tools.autocast_args
 702def kge(
 703    obs: np.ndarray | None = None,
 704    sim: np.ndarray | None = None,
 705    nodata: float = -99.0,
 706    t_axis: int = 0,
 707):
 708    """
 709    Compute the misfit criteria kge:
 710        std_obs = np.std(
 711            obs, axis=t_axis, where=mask_nodata, keepdims=True, mean=mean_obs
 712        )
 713        std_sim = np.std(
 714            sim, axis=t_axis, where=mask_nodata, keepdims=True, mean=mean_sim
 715        )
 716
 717    beta = mean_sim / mean_obs
 718
 719    alpha = std_sim / std_obs
 720
 721    r2 = (
 722        1
 723        / (sim.shape[t_axis])
 724        * np.sum(
 725            ((sim - mean_sim) / std_sim) * ((obs - mean_obs) / std_obs),
 726            axis=t_axis,
 727            keepdims=True,
 728        )
 729    )
 730
 731    kge = 1 - np.sqrt(
 732        (r2 - 1.0) ** 2.0 + (alpha - 1.0) ** 2.0 + (beta - 1) ** 2.0
 733    )
 734
 735    :param obs: Observed discharges, defaults to None
 736    :type obs: np.ndarray | None, optional
 737    :param sim: simulated discharges, defaults to None
 738    :type sim: np.ndarray | None, optional
 739    :param nodata: No data value, defaults to -99.0
 740    :type nodata: float, optional
 741    :param t_axis: Array axis of the time, defaults to 0
 742    :type t_axis: int, optional
 743
 744    """
 745    t_axis = min(t_axis, len(obs.shape) - 1)
 746    _test_input_shape(obs)
 747    _test_input_shape(sim)
 748
 749    mask_nodata = obs != nodata
 750    nb_valid_data = np.count_nonzero(mask_nodata, axis=t_axis)
 751
 752    kge = np.nan
 753
 754    if sum(nb_valid_data) > 0:
 755
 756        if isinstance(obs, np.ndarray) and isinstance(sim, np.ndarray):
 757            if obs.shape == sim.shape:
 758
 759                with warnings.catch_warnings():
 760                    warnings.simplefilter("ignore", category=RuntimeWarning)
 761                    mean_obs = np.mean(obs, axis=t_axis, where=mask_nodata, keepdims=True)
 762                    mean_sim = np.mean(sim, axis=t_axis, where=mask_nodata, keepdims=True)
 763                    std_obs = np.std(
 764                        obs, axis=t_axis, where=mask_nodata, keepdims=True, mean=mean_obs
 765                    )
 766                    std_sim = np.std(
 767                        sim, axis=t_axis, where=mask_nodata, keepdims=True, mean=mean_sim
 768                    )
 769
 770                mean_obs = np.where(mean_obs == 0, np.nan, mean_obs)
 771                beta = mean_sim / mean_obs
 772
 773                std_obs = np.where(std_obs == 0, np.nan, std_obs)
 774                alpha = std_sim / std_obs
 775
 776                # r2 = np.sum(
 777                #     (sim - mean_sim) * (obs - mean_obs), axis=t_axis, keepdims=True
 778                # ) / ((sim.shape[t_axis]) * (std_obs * std_sim))
 779
 780                r2 = (
 781                    1
 782                    / (sim.shape[t_axis])
 783                    * np.sum(
 784                        ((sim - mean_sim) / std_sim) * ((obs - mean_obs) / std_obs),
 785                        axis=t_axis,
 786                        keepdims=True,
 787                    )
 788                )
 789
 790                kge = 1 - np.sqrt(
 791                    (r2 - 1.0) ** 2.0 + (alpha - 1.0) ** 2.0 + (beta - 1) ** 2.0
 792                )
 793            else:
 794                raise ValueError(f"Error: len(obs)!=len(sim), {len(obs)}!={len(sim)}")
 795        else:
 796            raise ValueError("Error: obs and sim must be an instance of np.ndarray")
 797
 798    else:
 799        raise ValueError("Error: no valid observation data.")
 800
 801    return np.squeeze(kge)
 802
 803
 804@tools.autocast_args
 805def sm_kge(
 806    obs: np.ndarray | None = None,
 807    sim: np.ndarray | None = None,
 808):
 809    """
 810    Compute the misfit criteria kge, see the Smash documentation at https://smash.recover.inrae.fr/math_num_documentation/efficiency_error_metric.html
 811
 812    :param obs: Observed discharges, defaults to None
 813    :type obs: np.ndarray | None, optional
 814    :param sim: simulated discharges, defaults to None
 815    :type sim: np.ndarray | None, optional
 816    """
 817
 818    kge = smash_metrics.kge(obs, sim)
 819
 820    return kge
 821
 822
 823# Fonction de calcul du débit de retour
 824def quantile_gumbel(T: float = 1.0, loc: float = 0.0, scale: float = 0.0):
 825    """
 826    Compute the quantile for a given return period using the Gumbel law.
 827    :param T: The return period, defaults to 1
 828    :type T: float, optional
 829    :param loc: The localisation parameter of the Gumbel law, defaults to 0.0
 830    :type loc: float, optional
 831    :param scale: The scale parameter of the Gumbel law, defaults to 0.0
 832    :type scale: float, optional
 833    :return: The value of the quantile
 834    :rtype: float
 835
 836    """
 837
 838    return stats.gumbel_r.ppf(1 - 1 / T, loc=loc, scale=scale)
 839
 840
 841# Fonction de calcul du débit de retour
 842def quantile_gev(
 843    T: float = 1.0, shape: float = 0.0, loc: float = 0.0, scale: float = 0.0
 844):
 845    """
 846    Compute the quantile for a given return period using the GEV law.
 847    :param T: The return period, defaults to 1
 848    :type T: float, optional
 849    :param shape: The shape parameter of the GEV law, defaults to 0.0
 850    :type shape: float, optional
 851    :param loc: The localisation parameter of the GEV law, defaults to 0.0
 852    :type loc: float, optional
 853    :param scale: The scale parameter of the GEV law, defaults to 0.0
 854    :type scale: float, optional
 855    :return: The value of the quantile
 856    :rtype: float
 857
 858    """
 859    return stats.genextreme.ppf(1 - 1 / T, shape, loc=loc, scale=scale)
 860
 861
 862def genextreme_fit(data: np.ndarray | None = None, estimate_method: str = "MLE"):
 863    """
 864    Return estimates of shape, location, and scale parameters from data. The default
 865     estimation method is Maximum Likelihood Estimation (MLE), but Method of Moments (MM)
 866      is also available.
 867    :param data: Data of maximum values used to fit a GEV, defaults to None
 868    :type data: np.ndarray | None, optional
 869    :param estimate_method: Method to optimize the parameters: MLE for Maximum Likelihood
 870     Estimate, MM for  Method of Moments , defaults to "MLE"
 871    :type estimate_method: str, optional
 872    :return: Estimates for any shape parameters (if applicable), followed by those for
 873     location and scale.
 874    :rtype: tuple of float
 875
 876    """
 877    if data is None:
 878        raise ValueError("input data is None. You must provide valid data.")
 879    res = stats.genextreme.fit(data, method=estimate_method)
 880    return res
 881
 882
 883def gumbel_r_fit(data: np.ndarray | None = None, estimate_method: str = "MLE"):
 884    """
 885    Return estimates of shape, location, and scale parameters from data.
 886     The default estimation method is Maximum Likelihood Estimation (MLE),
 887      but Method of Moments (MM) is also available.
 888    :param data: Data of maximum values used to fit a Gumbel law, defaults to None
 889    :type data: np.ndarray | None, optional
 890    :param estimate_method: Method to optimize the parameters: MLE for Maximum Likelihood
 891     Estimate, MM for  Method of Moments , defaults to "MLE"
 892    :type estimate_method: str, optional
 893    :return: Estimates for any shape parameters (if applicable), followed by those for
 894     location and scale.
 895    :rtype: tuple of float
 896
 897    """
 898    if data is None:
 899        raise ValueError("input data is None. You must provide valid data.")
 900
 901    res = stats.gumbel_r.fit(data, method=estimate_method)
 902    return (0, *res)
 903
 904
 905@tools.autocast_args
 906def time_resample_array(
 907    array: np.ndarray | None,
 908    quantile_duration: int | float = 1,
 909    model_time_step: float = 3600,
 910    quantile_chunk_size: int | None = None,
 911    t_axis: int = 2,
 912):
 913    """
 914    Resample the discharges array for a given time-step.
 915    :param array: the matrix containing the discharge with shape (nbx, nby, nbts)
 916    :type array: np.ndarray | None
 917    :param quantile_duration: The duration of the quantile (hours), defaults to 1
 918    :type quantile_duration: int | float, optional
 919    :param model_time_step: the time-step of the Smash model (seconds), defaults to 3600
 920    :type model_time_step: float, optional
 921    :param quantile_chunk_size: the size of the quantile chunk in days
 922    :type quantile_chunk_size: int | None
 923    :param t_axis: The array axis direction of the time-step, defaults to 2
 924    :type t_axis: int, optional
 925    :return: The resampled array
 926    :rtype: np.ndarray
 927
 928    """
 929    # R workaround as normal number (int) are passed to float to python
 930    # t_axis = int(t_axis)
 931
 932    # resample array to the new duration
 933    if pd.Timedelta(hours=quantile_duration) > pd.Timedelta(seconds=model_time_step):
 934        print(
 935            f"</> Resampling array with time-step `{pd.Timedelta(seconds=model_time_step)}`"
 936            f" to time-step `{pd.Timedelta(hours=quantile_duration)}`"
 937        )
 938        chunk_size = int(
 939            pd.Timedelta(hours=quantile_duration) / pd.Timedelta(seconds=model_time_step)
 940        )
 941
 942        array_trans = np.moveaxis(array, t_axis, 0)  # Axe à la position 0
 943
 944        new_shape = (array_trans.shape[0] // chunk_size, chunk_size) + array_trans.shape[
 945            1:
 946        ]
 947
 948        array_trans_reshaped = array_trans[
 949            0 : chunk_size * (array_trans.shape[0] // chunk_size)
 950        ].reshape(new_shape)
 951
 952        array_trans_reshaped_mean = np.mean(array_trans_reshaped, axis=1)
 953
 954        if quantile_chunk_size is None:
 955            remainder = array_trans.shape[0] % chunk_size
 956        else:
 957            # The final shape must be extended only if it is not a multiple of quantile_chunk_size.
 958            remainder = (
 959                array_trans_reshaped_mean.shape[0]
 960                * pd.Timedelta(hours=quantile_duration)
 961                % pd.Timedelta(days=quantile_chunk_size)
 962            ).total_seconds()
 963
 964        del array_trans
 965
 966        if remainder > 0:
 967            # extend the array with the previous value
 968            array_trans_reshaped_mean = np.insert(
 969                array_trans_reshaped_mean,
 970                obj=-1,
 971                axis=0,
 972                values=array_trans_reshaped_mean[-1, :],
 973            )
 974
 975        del array_trans_reshaped
 976
 977        array = np.moveaxis(array_trans_reshaped_mean, 0, t_axis)
 978
 979        del array_trans_reshaped_mean
 980
 981    return array
 982
 983
 984@tools.autocast_args
 985def compute_maxima(
 986    array: np.ndarray | None,
 987    t_axis: int = 2,
 988    nb_minimum_chunks: int = 4,
 989    chunk_size: int = 365,
 990    quantile_duration: int | float = 1,
 991):
 992    """
 993    Compute the maxima of the discharges for a given chunk_size in the t_axis direction.
 994    :param array: the matrix containing the discharge with shape (nbx, nby, nbts). The
 995    shape can be smaller or higher. But t_axis must be set to target the nbts
 996    (number of time-step) dimension. Ex: if shape=(nbx, nbts), t_axis must be equal to 1.
 997    :type array: np.ndarray | None
 998    :param t_axis: The array axis direction of the time-step, defaults to 2, defaults to 2
 999    :type t_axis: int, optional
1000    :param nb_minimum_chunks: minimal number of chunk required adjust an extrem law and
1001    compute the quantile, defaults to 4
1002    :type nb_minimum_chunks: int, optional
1003    :param chunk_size: The chunk_size (days). It correspond to the 'unit' of the
1004    return period, defaults to 365 (year)
1005    :type chunk_size: int, optional
1006    :param quantile_duration: The duration of the quantile (hour), defaults to 1
1007    :type quantile_duration: int | float, optional
1008    :return: an np.ndarray of spatial maximal values for every chunk
1009    :rtype: np.ndarray
1010
1011    """
1012
1013    if pd.Timedelta(
1014        hours=quantile_duration,
1015    ) > pd.Timedelta(
1016        days=chunk_size,
1017    ):
1018        raise ValueError(
1019            "The chunk_size {chunk_size} (days) must be"
1020            " greater than the quantile duration {quantile_duration} (hours)"
1021        )
1022
1023    # compute the nb of chunk in the input data
1024    nbchunks = int(
1025        array.shape[t_axis]
1026        * pd.Timedelta(hours=quantile_duration)
1027        / pd.Timedelta(days=chunk_size)
1028    )
1029    print(f"</> Data contain {nbchunks} chunks of {pd.Timedelta(days=chunk_size)}")
1030
1031    nb_ts_by_chunk = int(array.shape[t_axis] / nbchunks)
1032
1033    if nbchunks < nb_minimum_chunks:
1034        raise ValueError(
1035            "</> The number of simulated chunks is not enough to compute the quantile:"
1036            f" {nbchunks}<{nb_minimum_chunks}"
1037        )
1038
1039    # nb year limit to compute the quantile
1040    if nbchunks >= nb_minimum_chunks:
1041        print(
1042            f"</> Compute the maxima for every chunk with size"
1043            f" {nb_ts_by_chunk} time-steps"
1044            f" of {quantile_duration} hours."
1045        )
1046        out_shape = []
1047        list_axes = []
1048        for nshape in range(len(array.shape)):
1049            out_shape.append(0)
1050            list_axes.append(nshape)
1051
1052        out_shape[t_axis] = nbchunks
1053        list_axes.remove(t_axis)
1054
1055        for ax in list_axes:
1056            out_shape[ax] = array.shape[ax]
1057
1058        maxima = np.zeros(shape=out_shape)
1059
1060        nb_ts_by_chunk = int(array.shape[t_axis] / nbchunks)
1061        remainder = array.shape[t_axis] % nbchunks
1062        pos = 0
1063
1064        # original_axis = [i for i in range(len(array.shape))]
1065        # other_axis = original_axis.copy()
1066        # other_axis.remove(t_axis)
1067        # destination_axis = [original_axis[t_axis]] + other_axis
1068
1069        # print(original_axis, destination_axis)
1070
1071        # print(maxima.shape)
1072        # print(array.shape)
1073        # move_axis first to extend the possibiltiy to use this function with different array shape
1074        array = np.moveaxis(array, t_axis, 0)
1075        maxima = np.moveaxis(maxima, t_axis, 0)
1076        # print(maxima.shape)
1077        # print(array.shape)
1078        for i in range(nbchunks):
1079            # Alternate one more cell because nb_ts_by_chunk is not a factor of array.shape[t_axis]
1080            if i % 2 > 0 and remainder > 0:
1081                onemore = 1
1082            else:
1083                onemore = 0
1084
1085            # maxima[:, :, i] = np.max(
1086            #     array[:, :, int(i * nb_ts_by_chunk) : int((i + 1) * nb_ts_by_chunk)],
1087            #     axis=t_axis,
1088            # )
1089            # print(array.shape, pos, int(pos + nb_ts_by_chunk + onemore))
1090            # maxima[:, :, i] = np.max(
1091            #     array[:, :, pos : int(pos + nb_ts_by_chunk + onemore)],
1092            #     axis=t_axis,
1093            # )
1094
1095            # Check if value >-99 and check if nb lacuna <20% ? Or remove lowers maxima after all (10%)
1096            # if (array[pos : int(pos + nb_ts_by_chunk + onemore), :])
1097
1098            maxima[i, :] = np.max(
1099                array[pos : int(pos + nb_ts_by_chunk + onemore), :],
1100                axis=0,
1101            )
1102            pos = int(pos + nb_ts_by_chunk + onemore)
1103
1104        array = np.moveaxis(array, 0, t_axis)
1105        maxima = np.moveaxis(maxima, 0, t_axis)
1106
1107    return maxima
1108
1109
1110def quantil_obs(
1111    qobs_directory: str | None = None,
1112    code: np.ndarray | list = [],
1113    model_time_step: float = 3600,
1114    nb_minimum_chunks: int = 4,
1115    chunk_size: int = 365,
1116    quantile_duration: int = 1,
1117):
1118
1119    if qobs_directory is None:
1120        print("</> qobs_directory `{qobs_directory}` is not a valid directory")
1121        return
1122
1123    if isinstance(code, list):
1124        code = np.array(code)
1125
1126    qobs = tools.read_hourly_qobs(qobs_directory, code)
1127
1128    array = time_resample_array(
1129        array=qobs,
1130        quantile_duration=quantile_duration,
1131        model_time_step=model_time_step,
1132        quantile_chunk_size=chunk_size,
1133        t_axis=1,
1134    )
1135
1136    maxima = compute_maxima(
1137        array=array,
1138        t_axis=1,
1139        nb_minimum_chunks=nb_minimum_chunks,
1140        chunk_size=chunk_size,
1141        quantile_duration=quantile_duration,
1142    )
1143
1144    results = empirical_obs_quantile(
1145        maxima=maxima,
1146        nb_minimum_chunks=nb_minimum_chunks,
1147        chunk_size=chunk_size,
1148        quantile_duration=quantile_duration,
1149    )
1150
1151    return results
1152
1153
1154def empirical_obs_quantile(
1155    maxima: np.ndarray | None,
1156    t_axis=1,
1157    nb_minimum_chunks: int = 4,
1158    frac_to_remove: float = 0.1,
1159    quantile_duration: int | float = 1,
1160    chunk_size: int = 365,
1161):
1162    """
1163    :param maxima: Maximal discharge by chunk
1164    :type maxima: np.ndarray | None
1165    :param t_axis: axis of the time series, defaults to 1
1166    :type t_axis: TYPE, optional
1167    :param nb_minimum_chunks: number of minimum chunck required to remove `frac_to_remove` of the maxima distribution, i.e len of the maxima array in the t_axis direction, defaults to 4
1168    :type nb_minimum_chunks: int, optional
1169    :param frac_to_remove: fraction between 0 and 1 to remove lower value of the maxima distribution, incase of incomplete data chunk, defaults to 0.1
1170    :type frac_to_remove: float, default 0.1
1171    :param quantile_duration: duration of the quantile (hour)
1172    :type quantile_duration int | float, default 1
1173    :param chunk_size: size of the chunk used to comute the maxima
1174    :type chunk_size: int, default 365
1175    """
1176
1177    maxima = np.moveaxis(maxima, t_axis, 0)
1178    maxima_sorted = np.sort(maxima, axis=0)
1179    maxima_sorted = np.where(maxima_sorted < 0, np.nan, maxima_sorted)
1180    T_emp = np.zeros(shape=maxima.shape) + np.nan
1181    # n = maxima.shape[0]
1182
1183    for sta in range(maxima_sorted.shape[1]):
1184        n_valid = np.where(maxima_sorted[:, sta] > 0)
1185        n = len(n_valid[0])
1186        # remove 10% of the
1187        if n > nb_minimum_chunks:
1188            n_to_remove = max(1, int(frac_to_remove * n))
1189        else:
1190            n_to_remove = 0
1191
1192        for i in range(n_to_remove):
1193            ind = n_valid[0][i]
1194            maxima_sorted[i] = np.nan
1195
1196        if n > nb_minimum_chunks:
1197            rank = 1
1198            for i in range(n_to_remove, n):
1199                probs = (rank - 0.5) / n
1200                T_emp[n_valid[0][i], sta] = 1 / (1 - probs)
1201                rank = rank + 1
1202
1203    # trim/filter nan cells
1204    index = []
1205    for t in range(maxima_sorted.shape[0]):
1206        if np.any(maxima_sorted[t, :] >= 0):
1207            index.append(t)
1208
1209    maxima_sorted = maxima_sorted[index, :]
1210    T_emp = T_emp[index, :]
1211
1212    maxima_sorted = np.moveaxis(maxima_sorted, 0, t_axis)
1213    T_emp = np.moveaxis(T_emp, 0, t_axis)
1214
1215    results = {
1216        "maxima": maxima_sorted,
1217        "Temp": T_emp,
1218        "chunk_size": chunk_size,
1219        "nb_chunks": maxima.shape[t_axis],
1220        "quantile_duration": quantile_duration,
1221    }
1222    return results
1223
1224
1225@tools.autocast_args
1226def fit_quantile(
1227    maxima: np.ndarray | None,
1228    t_axis: int = 2,
1229    return_periods: list | tuple = [2, 5, 10, 20, 50, 100],
1230    fit: str = "gumbel",
1231    estimate_method: str = "MLE",
1232    quantile_duration: int | float = 1,
1233    chunk_size: int = 365,
1234    ncpu: int | None = None,
1235):
1236    """
1237    Proceed to the qunaitle adjustment.
1238    :param maxima: an np.ndarray of spatial maximal values for every chunk
1239    :type maxima: np.ndarray | None
1240    :param t_axis: The array axis direction of the time-step, defaults to 2,
1241    defaults to 2, defaults to 2
1242    :type t_axis: int, optional
1243    :param return_periods: A list of the return period, defaults to [2, 5, 10, 20, 50, 100]
1244    :type return_periods: list | tuple, optional
1245    :param fit: The extrem law to use (gumbel or gev), defaults to "gumbel"
1246    :type fit: str, optional
1247    :param estimate_method: The methode to use for calibrated the parameters
1248    of the extrem law (MLE or MM), defaults to "MLE"
1249    :type estimate_method: str, optional
1250    :param quantile_duration: The duration of the quantile (hour), defaults to 1
1251    :type quantile_duration: int | float, optional
1252    :param chunk_size: the size of the chunks in days, defaults to 365
1253    :type chunk_size: int, optional
1254    :param ncpu: Number of cpu to use, defaults to None
1255    :type ncpu: int | None, optional
1256    :return: A dictionary containing the results od the quantile computation:
1257     Results include the quantile, the empirical return period,
1258     the maxima for each chunk of chunk_size, the fit parameters of the `fit` extrem law.
1259    :rtype: dict
1260
1261    """
1262
1263    if ncpu is None:
1264        ncpu = int(os.cpu_count() / 2)
1265    else:
1266        ncpu = int(min(ncpu, os.cpu_count() - 1))
1267
1268    if len(maxima.shape) < 2 or len(maxima.shape) > 3:
1269        raise ValueError(
1270            f"Input array `maxima` must have a dimension of 2 or 3. "
1271            f"Actual input have a dimension of {len(maxima.shape)} with shape {maxima.shape}"
1272        )
1273
1274    grid_shape = list(maxima.shape)
1275    del grid_shape[t_axis]
1276
1277    fit_shape = np.zeros(shape=grid_shape)
1278    fit_loc = np.zeros(shape=grid_shape)
1279    fit_scale = np.zeros(shape=grid_shape)
1280    quantile = np.zeros(shape=grid_shape + [len(return_periods)]) * np.nan
1281
1282    # list_axes = [ax for ax in range(len(maxima.shape))]
1283    # list_axes.remove(t_axis)
1284
1285    # fit_shape = np.zeros(shape=(maxima.shape[list_axes[0]], maxima.shape[list_axes[1]]))
1286    # fit_loc = np.zeros(shape=(maxima.shape[list_axes[0]], maxima.shape[list_axes[1]]))
1287    # fit_scale = np.zeros(shape=(maxima.shape[list_axes[0]], maxima.shape[list_axes[1]]))
1288
1289    # quantile = (
1290    #     np.zeros(
1291    #         shape=(
1292    #             maxima.shape[list_axes[0]],
1293    #             maxima.shape[list_axes[1]],
1294    #             len(return_periods),
1295    #         )
1296    #     )
1297    #     * np.nan
1298    # )
1299
1300    # Total length of the grid, linearize the matrix => 1 dim
1301    nd = 1
1302    for d in grid_shape:
1303        nd = nd * d
1304
1305    original_shape = maxima.shape
1306    newshape = (nd, maxima.shape[-1])
1307
1308    maxima = maxima.reshape(newshape)
1309    fit_shape = fit_shape.reshape(nd)
1310    fit_loc = fit_loc.reshape(nd)
1311    fit_scale = fit_scale.reshape(nd)
1312    quantile = quantile.reshape((nd, len(return_periods)))
1313
1314    print(f"</> Fitting {fit} law on data using {estimate_method} method.")
1315
1316    pool = multiprocessing.Pool(ncpu)
1317
1318    # Split by chunk to display a progress bar. chunk have the size of cpu
1319    chunksize = ncpu
1320    nbchunk = int(maxima.shape[0] / chunksize)
1321
1322    for chunk in tqdm(range(nbchunk + 1)):
1323
1324        args = []
1325        index_i = []
1326
1327        for i in range(chunk * chunksize, chunk * chunksize + chunksize):
1328
1329            if i >= nd:
1330                break
1331
1332            if np.all(maxima[i, :] > 0.0):
1333                args.append((maxima[i, :], estimate_method))
1334                index_i.append(i)
1335
1336        if fit == "gumbel":
1337            res = pool.starmap(gumbel_r_fit, args, chunksize=1)
1338
1339        if fit == "gev":
1340            res = pool.starmap(genextreme_fit, args, chunksize=1)
1341
1342        k = 0
1343        for item in res:
1344            i = index_i[k]
1345            fit_shape[i] = item[0]
1346            fit_loc[i] = item[1]
1347            fit_scale[i] = item[2]
1348
1349            for index, T in enumerate(return_periods):
1350                if fit == "gumbel":
1351                    quantile[i, index] = quantile_gumbel(T, item[1], item[2])
1352                if fit == "gev":
1353                    quantile[i, index] = quantile_gev(T, item[0], item[1], item[2])
1354
1355            k = k + 1
1356
1357    maxima = maxima.reshape(original_shape)
1358    fit_shape = fit_shape.reshape(grid_shape)
1359    fit_loc = fit_loc.reshape(grid_shape)
1360    fit_scale = fit_scale.reshape(grid_shape)
1361    quantile = quantile.reshape(grid_shape + [len(return_periods)])
1362
1363    # pool = multiprocessing.Pool(ncpu)
1364
1365    # print(f"</> Fitting {fit} law on data using {estimate_method} method.")
1366
1367    # for i in tqdm(range(maxima.shape[list_axes[0]])):
1368
1369    #     args = []
1370    #     index_j = []
1371
1372    #     for j in range(maxima.shape[list_axes[1]]):
1373    #         if np.all(maxima[i, j, :] > 0.0):
1374    #             args.append((maxima[i, j, :], estimate_method))
1375    #             index_j.append(j)
1376
1377    #     if fit == "gumbel":
1378    #         res = pool.starmap(gumbel_r_fit, args, chunksize=1)
1379
1380    #     if fit == "gev":
1381    #         res = pool.starmap(genextreme_fit, args, chunksize=1)
1382
1383    #     k = 0
1384    #     for item in res:
1385
1386    #         j = index_j[k]
1387
1388    #         fit_shape[i, j] = item[0]
1389    #         fit_loc[i, j] = item[1]
1390    #         fit_scale[i, j] = item[2]
1391
1392    #         for index, T in enumerate(return_periods):
1393    #             if fit == "gumbel":
1394    #                 quantile[i, j, index] = quantile_gumbel(T, item[1], item[2])
1395    #             if fit == "gev":
1396    #                 quantile[i, j, index] = quantile_gev(T, item[0], item[1], item[2])
1397
1398    #         k = k + 1
1399
1400    pool.close()
1401    pool.terminate()
1402
1403    n = maxima.shape[t_axis]
1404    ranks = np.arange(1, n + 1)
1405    # probs = (ranks - 0.44) / (n + 0.12)
1406    probs = (ranks - 0.5) / n
1407    T_emp = 1 / (1 - probs)
1408
1409    results = {
1410        "T": np.array(return_periods),
1411        "Q_th": quantile,
1412        "T_emp": T_emp,
1413        "maxima": maxima,
1414        "nb_chunks": maxima.shape[t_axis],
1415        "fit": fit,
1416        "fit_loc": fit_loc,
1417        "fit_scale": fit_scale,
1418        "fit_shape": fit_shape,
1419        "duration": quantile_duration,
1420        "chunk_size": chunk_size,
1421    }
1422
1423    return results
1424
1425
1426@tools.autocast_args
1427def parametric_bootstrap_uncertainties(
1428    q_results: dict | None,
1429    bootstrap_sample=100,
1430    estimate_method: str = "MLE",
1431    ncpu: int | None = None,
1432):
1433    """
1434    Proceed to the qunaitle adjustment.
1435    :param q_results: A dictionary containing the results of the quantile adjustment
1436    :type q_results: dict | None
1437    :param bootstrap_sample: the number of bootstrap sample, defaults to 100
1438    :type bootstrap_sample: int, optional
1439    :param estimate_method: The methode to use for calibrated the parameters
1440    of the extrem law (MLE or MM), defaults to "MLE"
1441    :type estimate_method: str, optional
1442    :param ncpu: Number of cpu to use, defaults to None
1443    :type ncpu: int | None, optional
1444    :return: A dictionary containing the results od the quantile computation:
1445     Results include the quantile, the empirical return period,
1446     the maxima for each chunk of chunk_size, the fit parameters of the `fit` extrem law.
1447    :rtype: dict
1448
1449    """
1450    # R workaround as normal number (int) are passed to float to python
1451    # ncpu = int(ncpu)
1452
1453    if ncpu is None:
1454        ncpu = int(os.cpu_count() / 2)
1455    else:
1456        ncpu = min(ncpu, os.cpu_count() - 1)
1457
1458    sample_size = q_results["nb_chunks"]
1459    quantile = q_results["Q_th"]
1460    fit_loc = q_results["fit_loc"]
1461    fit_scale = q_results["fit_scale"]
1462    fit_shape = q_results["fit_shape"]
1463
1464    grid_shape = list(quantile.shape)
1465    del grid_shape[-1]  # remove last dimension (return period)
1466
1467    parametric_quantile = (
1468        np.zeros(
1469            shape=(
1470                *grid_shape,
1471                len(q_results["T"]),
1472                bootstrap_sample,
1473            )
1474        )
1475        * np.nan
1476    )
1477
1478    nd = 1
1479    for d in grid_shape:
1480        nd = nd * d
1481
1482    # original_shape = parametric_quantile.shape
1483    newshape = (nd, len(q_results["T"]), bootstrap_sample)
1484
1485    parametric_quantile = parametric_quantile.reshape(newshape)
1486    fit_shape = fit_shape.reshape(nd)
1487    fit_loc = fit_loc.reshape(nd)
1488    fit_scale = fit_scale.reshape(nd)
1489    quantile = quantile.reshape(nd, len(q_results["T"]))
1490
1491    print(f"</> Compute bootstrap uncertainties...")
1492
1493    pool = multiprocessing.Pool(ncpu)
1494
1495    for i in tqdm(range(nd)):
1496
1497        if np.all(quantile[i, :] > 0.0):
1498
1499            args = []
1500            index_k = []
1501            for k in range(bootstrap_sample):
1502
1503                # generate sample
1504                if q_results["fit"] == "gumbel":
1505                    random_q = stats.gumbel_r.rvs(
1506                        loc=fit_loc[i],
1507                        scale=fit_scale[i],
1508                        size=sample_size,
1509                    )
1510                if q_results["fit"] == "gev":
1511                    random_q = stats.genextreme_fit.rvs(
1512                        fit_shape[i],
1513                        loc=fit_loc[i],
1514                        scale=fit_scale[i],
1515                        size=sample_size,
1516                    )
1517
1518                args.append((random_q, q_results["fit"]))
1519                index_k.append(k)
1520
1521            if q_results["fit"] == "gumbel":
1522                res = pool.starmap(gumbel_r_fit, args, chunksize=1)
1523
1524            if q_results["fit"] == "gev":
1525                res = pool.starmap(genextreme_fit, args, chunksize=1)
1526
1527            for k, item in enumerate(res):
1528
1529                ik = index_k[k]
1530
1531                for indexT, T in enumerate(q_results["T"]):
1532                    if q_results["fit"] == "gumbel":
1533                        parametric_quantile[i, indexT, ik] = quantile_gumbel(
1534                            T, item[1], item[2]
1535                        )
1536                    if q_results["fit"] == "gev":
1537                        parametric_quantile[i, indexT, ik] = quantile_gev(
1538                            T, item[0], item[1], item[2]
1539                        )
1540
1541    Umax = np.max(parametric_quantile, axis=2)
1542    Umin = np.min(parametric_quantile, axis=2)
1543
1544    Umax = Umax.reshape(grid_shape + [len(q_results["T"])])
1545    Umin = Umin.reshape(grid_shape + [len(q_results["T"])])
1546
1547    # parametric_quantile = parametric_quantile.reshape(original_shape)
1548    # fit_shape = fit_shape.reshape(grid_shape)
1549    # fit_loc = fit_loc.reshape(grid_shape)
1550    # fit_scale = fit_scale.reshape(grid_shape)
1551    # quantile = quantile.reshape(grid_shape + [len(q_results["T"])])
1552
1553    ##############################################"
1554    # for i in tqdm(range(q_results["Q_th"].shape[0])):
1555    #     for j in range(q_results["Q_th"].shape[1]):
1556
1557    #         if np.all(q_results["Q_th"][i, j] > 0.0):
1558
1559    #             args = []
1560    #             index_k = []
1561    #             for k in range(bootstrap_sample):
1562
1563    #                 # generate sample
1564    #                 if q_results["fit"] == "gumbel":
1565    #                     random_q = stats.gumbel_r.rvs(
1566    #                         loc=q_results["fit_loc"][i, j],
1567    #                         scale=q_results["fit_scale"][i, j],
1568    #                         size=sample_size,
1569    #                     )
1570    #                 if q_results["fit"] == "gev":
1571    #                     random_q = stats.genextreme_fit.rvs(
1572    #                         q_results["fit_shape"][i, j],
1573    #                         loc=q_results["fit_loc"][i, j],
1574    #                         scale=q_results["fit_scale"][i, j],
1575    #                         size=sample_size,
1576    #                     )
1577
1578    #                 args.append((random_q, q_results["fit"]))
1579    #                 index_k.append(k)
1580
1581    #             if q_results["fit"] == "gumbel":
1582    #                 res = pool.starmap(gumbel_r_fit, args, chunksize=1)
1583
1584    #             if q_results["fit"] == "gev":
1585    #                 res = pool.starmap(genextreme_fit, args, chunksize=1)
1586
1587    #             # k = 0
1588    #             for k, item in enumerate(res):
1589
1590    #                 ik = index_k[k]
1591
1592    #                 for indexT, T in enumerate(return_periods):
1593    #                     if q_results["fit"] == "gumbel":
1594    #                         parametric_quantile[i, j, indexT, ik] = quantile_gumbel(
1595    #                             T, item[1], item[2]
1596    #                         )
1597    #                     if q_results["fit"] == "gev":
1598    #                         parametric_quantile[i, j, indexT, ik] = quantile_gev(
1599    #                             T, item[0], item[1], item[2]
1600    #                         )
1601
1602    pool.close()
1603    pool.terminate()
1604
1605    q_results.update({"Umax": Umax, "Umin": Umin})
1606
1607    return q_results
1608
1609
1610@tools.autocast_args
1611def fit_quantile_unparallel(
1612    maxima: np.ndarray | None,
1613    t_axis: int = 2,
1614    return_periods: list | tuple = [2, 5, 10, 20, 50, 100],
1615    fit: str = "gumbel",
1616    estimate_method: str = "MLE",
1617    quantile_duration: int | float = 1.0,
1618    chunk_size: int = 365,
1619):
1620    """
1621    Proceed to the qunaitle adjustment (unparrallel computation).
1622    :param maxima: an np.ndarray of spatial maximal values for every chunk
1623    :type maxima: np.ndarray | None
1624    :param t_axis: The array axis direction of the time-step, defaults to 2,
1625    defaults to 2, defaults to 2
1626    :type t_axis: int, optional
1627    :param return_periods: A list of the return period, defaults to [2, 5, 10, 20, 50, 100]
1628    :type return_periods: list | tuple, optional
1629    :param fit: The extrem law to use (gumbel or gev), defaults to "gumbel"
1630    :type fit: str, optional
1631    :param estimate_method: The methode to use for calibrated the parameters
1632    of the extrem law (MLE or MM), defaults to "MLE"
1633    :type estimate_method: str, optional
1634    :param quantile_duration: The duration of the quantile (hour), defaults to 1
1635    :type quantile_duration: int | float, optional
1636    :param chunk_size: the size of the chunks in days, defaults to 365
1637    :type chunk_size: int, optional
1638    :return: A dictionary containing the results od the quantile computation:
1639     Results include the quantile, the empirical return period,
1640     the maxima for each chunk of chunk_size, the fit parameters of the `fit` extrem law.
1641    :rtype: dict
1642
1643    """
1644    # # R workaround as normal number (int) are passed to float to python
1645    # t_axis = int(t_axis)
1646    # chunk_size = int(chunk_size)
1647
1648    list_axes = [0, 1, 2]
1649    list_axes.remove(t_axis)
1650
1651    fit_shape = np.zeros(shape=(maxima.shape[list_axes[0]], maxima.shape[list_axes[1]]))
1652    fit_loc = np.zeros(shape=(maxima.shape[list_axes[0]], maxima.shape[list_axes[1]]))
1653    fit_scale = np.zeros(shape=(maxima.shape[list_axes[0]], maxima.shape[list_axes[1]]))
1654
1655    quantile = (
1656        np.zeros(
1657            shape=(
1658                maxima.shape[list_axes[0]],
1659                maxima.shape[list_axes[1]],
1660                len(return_periods),
1661            )
1662        )
1663        * np.nan
1664    )
1665
1666    print(f"</> Fitting {fit} law on data using {estimate_method} method.")
1667
1668    for i in range(maxima.shape[list_axes[0]]):
1669
1670        for j in range(maxima.shape[list_axes[1]]):
1671
1672            if np.all(maxima[i, j, :] > 0.0):
1673
1674                if fit == "gumbel":
1675                    res = gumbel_r_fit(maxima[i, j, :], estimate_method)
1676
1677                if fit == "gev":
1678                    res = genextreme_fit(maxima[i, j, :], estimate_method)
1679
1680                fit_shape[i, j] = res[0]
1681                fit_loc[i, j] = res[1]
1682                fit_scale[i, j] = res[2]
1683
1684                for index, T in enumerate(return_periods):
1685                    if fit == "gumbel":
1686                        quantile[i, j, index] = quantile_gumbel(T, res[1], res[2])
1687                    if fit == "gev":
1688                        quantile[i, j, index] = quantile_gev(T, res[0], res[1], res[2])
1689
1690    n = maxima.shape[t_axis]
1691    ranks = np.arange(1, n + 1)
1692    # probs = (ranks - 0.44) / (n + 0.12)
1693    probs = (ranks - 0.5) / n
1694    T_emp = 1 / (1 - probs)
1695
1696    results = {
1697        "T": np.array(return_periods),
1698        "Q_th": quantile,
1699        "T_emp": T_emp,
1700        "maxima": maxima,
1701        "nb_chunks": maxima.shape[t_axis],
1702        "fit": fit,
1703        "fit_loc": fit_loc,
1704        "fit_scale": fit_scale,
1705        "fit_shape": fit_shape,
1706        "duration": quantile_duration,
1707        "chunk_size": chunk_size,
1708    }
1709
1710    return results
1711
1712
1713@tools.autocast_args
1714def parametric_bootstrap_uncertainties_unparallel(
1715    q_results: dict | None,
1716    bootstrap_sample=100,
1717    return_periods: list | tuple = [2, 5, 10, 20, 50, 100],
1718    estimate_method: str = "MLE",
1719):
1720    """
1721    Proceed to the qunaitle adjustment.
1722    :param q_results: A dictionary containing the results of the quantile adjustment
1723    :type q_results: dict | None
1724    :param bootstrap_sample: the number of bootstrap sample, defaults to 100
1725    :type bootstrap_sample: int, optional
1726    :param return_periods: A list of the return period, defaults to [2, 5, 10, 20, 50, 100]
1727    :type return_periods: list | tuple, optional
1728    :param estimate_method: The methode to use for calibrated the parameters
1729    of the extrem law (MLE or MM), defaults to "MLE"
1730    :type estimate_method: str, optional
1731    :return: A dictionary containing the results od the quantile computation:
1732     Results include the quantile, the empirical return period,
1733     the maxima for each chunk of chunk_size, the fit parameters of the `fit` extrem law.
1734    :rtype: dict
1735
1736    """
1737
1738    sample_size = q_results["nb_chunks"]
1739
1740    parametric_quantile = (
1741        np.zeros(
1742            shape=(
1743                q_results["Q_th"].shape[0],
1744                q_results["Q_th"].shape[1],
1745                len(q_results["T"]),
1746                bootstrap_sample,
1747            )
1748        )
1749        * np.nan
1750    )
1751
1752    for i in tqdm(range(q_results["Q_th"].shape[0])):
1753        for j in range(q_results["Q_th"].shape[1]):
1754
1755            if np.all(q_results["Q_th"][i, j] > 0.0):
1756
1757                for k in range(bootstrap_sample):
1758
1759                    # generate sample
1760                    if q_results["fit"] == "gumbel":
1761                        random_q = stats.gumbel_r.rvs(
1762                            loc=q_results["fit_loc"][i, j],
1763                            scale=q_results["fit_scale"][i, j],
1764                            size=sample_size,
1765                        )
1766                    if q_results["fit"] == "gev":
1767                        random_q = stats.genextreme_fit.rvs(
1768                            q_results["shape_loc"][i, j],
1769                            loc=q_results["fit_loc"][i, j],
1770                            scale=q_results["fit_scale"][i, j],
1771                            size=sample_size,
1772                        )
1773
1774                    if q_results["fit"] == "gumbel":
1775                        res = gumbel_r_fit(random_q, estimate_method)
1776
1777                    if q_results["fit"] == "gev":
1778                        res = genextreme_fit(random_q, estimate_method)
1779
1780                    for indexT, T in enumerate(return_periods):
1781                        if q_results["fit"] == "gumbel":
1782                            parametric_quantile[i, j, indexT, k] = quantile_gumbel(
1783                                T, res[1], res[2]
1784                            )
1785                        if q_results["fit"] == "gev":
1786                            parametric_quantile[i, j, indexT, k] = quantile_gev(
1787                                T, res[0], res[1], res[2]
1788                            )
1789
1790    Umax = np.max(parametric_quantile, axis=3)
1791    Umin = np.min(parametric_quantile, axis=3)
1792
1793    q_results.update({"Umax": Umax, "Umin": Umin})
1794
1795    return q_results
1796
1797
1798@tools.autocast_args
1799def spatial_quantiles(
1800    array: np.ndarray | None,
1801    t_axis: int = 2,
1802    return_periods: list | tuple = [2, 5, 10, 20, 50, 100],
1803    fit: str = "gumbel",
1804    nb_minimum_chunks: int = 4,
1805    model_time_step: float = 3600,
1806    estimate_method: str = "MLE",
1807    chunk_size: int = 365,
1808    quantile_duration: int | float = 1,
1809    ncpu: int | None = None,
1810    compute_uncertainties=False,
1811    bootstrap_sample=100,
1812    maxima: np.ndarray | None = None,
1813):
1814    """
1815    Proceed to the quantile adjustment (parrallel computation).
1816    :param array: an np.ndarray of spatial dicharges values for every time-step
1817    :type array: np.ndarray | None
1818    :param t_axis: The array axis direction of the time-step, defaults to 2,
1819    defaults to 2, defaults to 2
1820    :type t_axis: int, optional
1821    :param return_periods: A list of the return period, defaults to [2, 5, 10, 20, 50, 100]
1822    :type return_periods: list | tuple, optional
1823    :param fit: The extrem law to use (gumbel or gev), defaults to "gumbel"
1824    :type fit: str, optional
1825    :param estimate_method: The methode to use for calibrated the parameters
1826    of the extrem law (MLE or MM), defaults to "MLE"
1827    :type estimate_method: str, optional
1828    :param quantile_duration: The duration of the quantile (hour), defaults to 1
1829    :type quantile_duration: int | float, optional
1830    :param chunk_size: the size of the chunks in days, defaults to 365
1831    :type chunk_size: int, optional
1832    :param ncpu: Number of cpu to use, defaults to None
1833    :type ncpu: int | None, optional
1834    :param compute_uncertainties: Compute the uncertainties usung the parametric bootstrap method
1835    :type compute_uncertainties: bool, True
1836    :param bootstrap_sample: the size of bootstrap sample, default is 100
1837    :type bootstrap_sample: int, optional
1838    :return: A dictionary containing the results od the quantile computation:
1839     Results include the quantile, the empirical return period,
1840     the maxima for each chunk of chunk_size, the fit parameters of the `fit` extrem law.
1841    :rtype: dict
1842
1843    """
1844    # # R workaround as normal number (int) are passed to float to python
1845    # t_axis = int(t_axis)
1846    # chunk_size = int(chunk_size)
1847
1848    if ncpu is None:
1849        ncpu = int(os.cpu_count() / 2)
1850    else:
1851        ncpu = int(min(ncpu, os.cpu_count() - 1))
1852
1853    if pd.Timedelta(
1854        hours=quantile_duration,
1855    ) < pd.Timedelta(
1856        seconds=model_time_step,
1857    ):
1858        raise ValueError(
1859            "The quantile duration {quantile_duration} (hours) must be"
1860            "greater or equal than the model time step {model_time_step} (seconds)"
1861        )
1862
1863    if pd.Timedelta(
1864        hours=quantile_duration,
1865    ) > pd.Timedelta(
1866        days=chunk_size,
1867    ):
1868        raise ValueError(
1869            "The chunk_size {chunk_size} (days) must be"
1870            " greater or equal than the quantile duration {quantile_duration} (hours)"
1871        )
1872
1873    if array is None:
1874        raise ValueError("Input array is None, no data to compute quantile.")
1875
1876    if maxima is None:
1877        array = time_resample_array(
1878            array=array,
1879            quantile_duration=quantile_duration,
1880            model_time_step=model_time_step,
1881            quantile_chunk_size=chunk_size,
1882            t_axis=t_axis,
1883        )
1884
1885        maxima = compute_maxima(
1886            array=array,
1887            t_axis=t_axis,
1888            nb_minimum_chunks=nb_minimum_chunks,
1889            chunk_size=chunk_size,
1890            quantile_duration=quantile_duration,
1891        )
1892
1893    results = fit_quantile(
1894        maxima=maxima,
1895        t_axis=t_axis,
1896        return_periods=return_periods,
1897        fit=fit,
1898        estimate_method=estimate_method,
1899        quantile_duration=quantile_duration,
1900        chunk_size=chunk_size,
1901        ncpu=ncpu,
1902    )
1903
1904    if compute_uncertainties:
1905        # calcul des incertitudes méthode bootstrap
1906        results = parametric_bootstrap_uncertainties(
1907            q_results=results,
1908            bootstrap_sample=bootstrap_sample,
1909            estimate_method=estimate_method,
1910            ncpu=ncpu,
1911        )
1912
1913    return results
1914
1915
1916@tools.autocast_args
1917def spatial_quantiles_unparallel(
1918    array: np.ndarray | None = None,
1919    t_axis: int = 2,
1920    return_periods: list | tuple = [2, 5, 10, 20, 50, 100],
1921    fit: str = "gumbel",
1922    nb_minimum_chunks: int = 4,
1923    model_time_step: float = 3600,
1924    estimate_method: str = "MLE",
1925    chunk_size: int = 365,
1926    compute_uncertainties=False,
1927    bootstrap_sample=100,
1928    quantile_duration: int | float = 1,
1929    maxima: np.ndarray | None = None,
1930):
1931    """
1932    Proceed to the quantile adjustment (unparrallel computation).
1933    :param array: an np.ndarray of spatial dicharges values for every time-step
1934    :type array: np.ndarray | None
1935    :param t_axis: The array axis direction of the time-step, defaults to 2,
1936    defaults to 2, defaults to 2
1937    :type t_axis: int, optional
1938    :param return_periods: A list of the return period, defaults to [2, 5, 10, 20, 50, 100]
1939    :type return_periods: list | tuple, optional
1940    :param fit: The extrem law to use (gumbel or gev), defaults to "gumbel"
1941    :type fit: str, optional
1942    :param estimate_method: The methode to use for calibrated the parameters
1943    of the extrem law (MLE or MM), defaults to "MLE"
1944    :type estimate_method: str, optional
1945    :param quantile_duration: The duration of the quantile (hour), defaults to 1
1946    :type quantile_duration: int | float, optional
1947    :param chunk_size: the size of the chunks in days, defaults to 365
1948    :type chunk_size: int, optional
1949    :param compute_uncertainties: Compute the uncertainties usung the parametric bootstrap method
1950    :type compute_uncertainties: bool, True
1951    :param bootstrap_sample: the size of bootstrap sample, default is 100
1952    :type bootstrap_sample: int, optional
1953    :return: A dictionary containing the results od the quantile computation:
1954     Results include the quantile, the empirical return period,
1955     the maxima for each chunk of chunk_size, the fit parameters of the `fit` extrem law.
1956    :rtype: dict
1957
1958    """
1959    # # R workaround as normal number (int) are passed to float to python
1960    # t_axis = int(t_axis)
1961    # chunk_size = int(chunk_size)
1962
1963    if pd.Timedelta(
1964        hours=quantile_duration,
1965    ) < pd.Timedelta(
1966        seconds=model_time_step,
1967    ):
1968        raise ValueError(
1969            "The quantile duration {quantile_duration} (hours) must be"
1970            "greater or equal than the model time step {model_time_step} (seconds)"
1971        )
1972
1973    if pd.Timedelta(
1974        hours=quantile_duration,
1975    ) > pd.Timedelta(
1976        days=chunk_size,
1977    ):
1978        raise ValueError(
1979            "The chunk_size {chunk_size} (days) must be"
1980            " greater or equal than the quantile duration {quantile_duration} (hours)"
1981        )
1982
1983    if array is None:
1984        raise ValueError("Input array is None, no data to compute quantile.")
1985
1986    if maxima is None:
1987        array = time_resample_array(
1988            array=array,
1989            quantile_duration=quantile_duration,
1990            model_time_step=model_time_step,
1991            quantile_chunk_size=chunk_size,
1992            t_axis=t_axis,
1993        )
1994
1995        maxima = compute_maxima(
1996            array=array,
1997            t_axis=t_axis,
1998            nb_minimum_chunks=nb_minimum_chunks,
1999            chunk_size=chunk_size,
2000            quantile_duration=quantile_duration,
2001        )
2002
2003    results = fit_quantile_unparallel(
2004        maxima=maxima,
2005        t_axis=t_axis,
2006        return_periods=return_periods,
2007        fit=fit,
2008        estimate_method=estimate_method,
2009        quantile_duration=quantile_duration,
2010        chunk_size=chunk_size,
2011    )
2012
2013    if compute_uncertainties:
2014        # calcul des incertitudes méthode bootstrap
2015        parametric_bootstrap_uncertainties_unparallel(
2016            q_results=results,
2017            bootstrap_sample=bootstrap_sample,
2018            return_periods=return_periods,
2019            estimate_method=estimate_method,
2020        )
2021
2022    return results