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