smashbox.stats.mystats

   1from smashbox.stats import stats
   2from smashbox.tools import tools
   3import numpy as np
   4import multiprocessing
   5import os
   6from functools import partial
   7from tqdm import tqdm
   8import pandas as pd
   9from smash.fcore import _mwd_metrics as smash_metrics
  10
  11
  12class mystats:
  13    """
  14    The class mystats includes functions and children classes to compute basics statistics
  15    with the results of the smash model.
  16    """
  17
  18    def __init__(self, parent_class):
  19        self._parent_class = parent_class
  20        """_parent_class attribute stores the parent_class src.model() to be able to
  21        access to the result of the smash simulation"""
  22
  23        self.misfit_stats = misfit_stats(self)
  24        """Attribute misfit_stats stores the class src.mystats.misfit_stats(). Its goal
  25         is to compute misfit criteria between simulated and observed discharges"""
  26        self.quantile_stats = spatial_quantile()
  27        """Attribute quantile_stats owns the class src.mystats.spatial_quantile().
  28         Its goal is to compute the discharges quantiles for various return period."""
  29        self.spatial_stats = spatial_stats(self)
  30        """Atrribute spatial_stats owns the class src.mystats.spatial_stats(). Its goal
  31        is to provide basic statistics on the discharges field over the time
  32        (mean, median, q20, q80, maximum, minimum)"""
  33        self.outlets_stats = outlets_stats(self)
  34        """Atrribute outlets_stats owns the class src.mystats.outlets_stats().
  35        Its goal is to provide basic statistics on the discharges at every outlets over the
  36        time (mean, median, q20, q80, maximum, minimum)"""
  37
  38    def fmisfit_stats(
  39        self, nodata=-99.0, column=[], ret=False, use_smash_metrics=True
  40    ):
  41        """
  42        Compute the misfit for every outlets between the simulated and the
  43         observed discharges
  44        :param nodata: No data values, defaults to -99.0
  45        :type nodata: TYPE, optional
  46        :param column: column on which to compute the statistics (gauge), defaults to []
  47        :type column: TYPE, optional
  48        :param ret: return the result, defaults to False
  49        :type ret: TYPE, optional
  50        :return: object with attributes with different statistics.
  51        :rtype: class src.mystats.misfit.results()
  52
  53        """
  54
  55        if not use_smash_metrics:
  56            self.misfit_stats.se(nodata=nodata, column=column)
  57            self.misfit_stats.mse(nodata=nodata, column=column)
  58            self.misfit_stats.rmse(nodata=nodata, column=column)
  59            self.misfit_stats.nrmse(nodata=nodata, column=column)
  60            self.misfit_stats.mae(nodata=nodata, column=column)
  61            self.misfit_stats.mape(nodata=nodata, column=column)
  62            self.misfit_stats.lgrm(nodata=nodata, column=column)
  63            self.misfit_stats.nse(nodata=nodata, column=column)
  64            self.misfit_stats.nnse(nodata=nodata, column=column)
  65            self.misfit_stats.kge(nodata=nodata, column=column)
  66
  67            self.misfit_stats.pearson(nodata=nodata, column=column)
  68        else:
  69            self.misfit_stats.sm_se(column=column)
  70            self.misfit_stats.sm_mse(column=column)
  71            self.misfit_stats.sm_rmse(column=column)
  72            self.misfit_stats.sm_nrmse(column=column)
  73            self.misfit_stats.sm_mae(column=column)
  74            self.misfit_stats.sm_mape(column=column)
  75            self.misfit_stats.sm_lgrm(column=column)
  76            self.misfit_stats.sm_nse(column=column)
  77            self.misfit_stats.sm_nnse(column=column)
  78            self.misfit_stats.sm_kge(column=column)
  79
  80            self.misfit_stats.pearson(nodata=nodata, column=column)
  81
  82        if ret:
  83            return self.misfit.results
  84
  85    def fspatial_stats(self, ret=False):
  86        """
  87        Compute basics statistics over the time (mean, median, q20, q80, maximum, minimum)
  88         on the spatial discharges field.
  89        :param ret: return the result, defaults to False
  90        :type ret: TYPE, optional
  91        :return: object with attributes with different statistics.
  92        :rtype: class src.mystats.spatial_stats.results()
  93
  94        """
  95
  96        if not hasattr(self._parent_class.extra_smash_results, "q_domain"):
  97            raise ValueError(
  98                "No smash extra results  'q_domain' found. Run forward_run()"
  99                "with return_options={'q_domain': True}"
 100            )
 101
 102        self.spatial_stats.mean()
 103        self.spatial_stats.median()
 104        self.spatial_stats.q20()
 105        self.spatial_stats.q80()
 106        self.spatial_stats.maximum()
 107        self.spatial_stats.minimum()
 108        self.spatial_stats.var()
 109
 110        if ret:
 111            return self.spatial_stats.results
 112
 113    def foutlets_stats(self, ret=False):
 114        """
 115        Compute basices statistics over the time (mean, median, q20, q80, maximum, minimum)
 116         at every outlets.
 117        :param ret: return the result, defaults to False
 118        :type ret: TYPE, optional
 119        :return: object with attributes with different statistics.
 120        :rtype: class src.mystats.outlets_stats.results()
 121
 122        """
 123        if self._parent_class.smash is None:
 124            raise ValueError(
 125                "Attribut smash is None. Perhaps, you forget to buil and run the"
 126                "smash model..."
 127            )
 128
 129        self.outlets_stats.mean()
 130        self.outlets_stats.median()
 131        self.outlets_stats.q20()
 132        self.outlets_stats.q80()
 133        self.outlets_stats.maximum()
 134        self.outlets_stats.minimum()
 135        self.outlets_stats.var()
 136
 137        if ret:
 138            return self.outlets_stats.results
 139
 140    @tools.autocast_args
 141    def fmaxima_stats(
 142        self,
 143        t_axis: int = 2,
 144        nb_minimum_chunks: int = 4,
 145        chunk_size: int = 365,
 146        quantile_duration: list | tuple = [1, 2, 3, 4, 6, 12, 24, 48, 72],
 147        cumulated_maxima: bool = True,
 148    ):
 149        """
 150        Compute the maximum discharge values of a 3D array by chunk, corresponding to `chunk_size` (days), along axis `t_axis` (time) for each pixel of the array.(nbX, nbY).
 151
 152        Parameters
 153        ----------
 154
 155        t_axis : int
 156            The axis along with the maximum will be computed. This axis should correspond
 157            to the time.
 158        nb_minimum_chunks: int
 159            number of minimum chunks required to compute the maxima along
 160        the t_axis. Default is set to 4. If the number of chunks is lower, the function
 161        will return None.
 162        chunk_size: int
 163            Size of the chunks in days. Default is 365 days.
 164        quantile_duration: list | tuple
 165            The duration of every quantile (hours). The discharges will be resampled for every duration. Default is [1, 2, 3, 4, 6, 12, 24, 48, 72] hours.
 166        cumulated_maxima: bool
 167            For each call of the function, the maxima will accumulate in a matrix for
 168            each quantil duration. This provide a convient way to compute the quantile
 169            (fit gumbel/gev) after many successive simulations.
 170
 171
 172        Examples
 173        --------
 174        >>> import smashbox
 175        >>> import numpy as np
 176        >>>
 177        >>> graffas_prcp = np.zeros(shape=(142, 166, 1300))
 178        >>> rng = np.random.default_rng()
 179        >>> graffas_prcp = (
 180        >>>     graffas_prcp
 181        >>>     + rng.multinomial(20, [0.2, 0.8], size=graffas_prcp.shape)[:, :, :, 0]
 182        >>> )
 183        >>>
 184        >>> es=smashbox.SmashBox()
 185        >>> sb.newmodel("graffas_zone")
 186        >>> sb.graffas_zone.generate_mesh()
 187        >>> sb.graffas_zone.atmos_data_connector(input_prcp=graffas_prcp)
 188        >>> sb.graffas_zone.model()
 189        >>> sb.graffas_zone.forward_run(invert_states=True, return_options={"q_domain": True})
 190        >>>
 191        >>> sb.graffas_zone.mysmashmodel.mystats.stats_maxima(chunk_size=10)
 192        >>> sb.graffas_zone.mysmashmodel.mystats.stats_maxima(chunk_size=10)
 193        >>>
 194        >>> results = stats.fit_quantile(
 195        >>>    maxima=es.graffas_zone.mysmashmodel.mystats.spatial_quantile.spatial_cumulated_maxima[
 196        >>>        :, :, :, 0
 197        >>>    ],
 198        >>>    t_axis=2,
 199        >>>    return_periods=[2, 5, 10, 20, 50, 100],
 200        >>>    fit="gumbel",
 201        >>>    estimate_method="MLE",
 202        >>>    quantile_duration=1,
 203        >>>    ncpu=6,
 204        >>>)
 205
 206        """
 207        if pd.Timedelta(
 208            hours=max(quantile_duration),
 209        ) > pd.Timedelta(
 210            days=chunk_size,
 211        ):
 212            raise ValueError(
 213                f"The chunk_size {chunk_size} (days) must be"
 214                f" greater or equal than the quantile duration {max(quantile_duration)} (hours)"
 215            )
 216
 217        if not hasattr(self._parent_class.extra_smash_results, "q_domain"):
 218            raise ValueError(
 219                "No smash extra results  'q_domain' found. Run forward_run()"
 220                "with return_options={'q_domain': True}"
 221            )
 222
 223        all_maxima = None
 224        for id_dur, duration in enumerate(quantile_duration):
 225            array = stats.time_resample_array(
 226                array=self._parent_class.extra_smash_results.q_domain,
 227                quantile_duration=duration,
 228                model_time_step=self._parent_class.smash.setup.dt,
 229                quantile_chunk_size=chunk_size,
 230                t_axis=t_axis,
 231            )
 232
 233            maxima = stats.compute_maxima(
 234                array=array,
 235                t_axis=t_axis,
 236                nb_minimum_chunks=nb_minimum_chunks,
 237                chunk_size=chunk_size,
 238                quantile_duration=duration,
 239            )
 240
 241            if all_maxima is None:
 242                all_maxima = (
 243                    np.zeros(shape=(*maxima.shape, len(quantile_duration)))
 244                    + np.nan
 245                )
 246                all_maxima_outlets = (
 247                    np.zeros(
 248                        shape=(
 249                            len(self._parent_class.smash.mesh.code),
 250                            maxima.shape[t_axis],
 251                            len(quantile_duration),
 252                        )
 253                    )
 254                    + np.nan
 255                )
 256
 257            all_maxima[:, :, :, id_dur] = maxima
 258
 259            for i in range(len(self._parent_class.smash.mesh.code)):
 260                coords = self._parent_class.smash.mesh.gauge_pos[i]
 261                all_maxima_outlets[i, :, id_dur] = all_maxima[
 262                    coords[0], coords[1], :, id_dur
 263                ]
 264
 265        if cumulated_maxima:
 266            # if hasattr(self.quantile_stats, "spatial_cumulated_maxima"):
 267            #     self.quantile_stats.spatial_cumulated_maxima = np.concat(
 268            #         (self.quantile_stats.spatial_cumulated_maxima, all_maxima),
 269            #         axis=t_axis,
 270            #     )
 271            #     self.quantile_stats.spatial_cumulated_maxima_outlets = np.concat(
 272            #         (
 273            #             self.quantile_stats.spatial_cumulated_maxima_outlets,
 274            #             all_maxima_outlets,
 275            #         ),
 276            #         axis=1,
 277            #     )
 278            # else:
 279            #     setattr(self.quantile_stats, "spatial_cumulated_maxima", all_maxima)
 280            #     setattr(
 281            #         self.quantile_stats,
 282            #         "spatial_cumulated_maxima_outlets",
 283            #         all_maxima_outlets,
 284            #     )
 285
 286            if self.quantile_stats.spatial_cumulated_maxima is None:
 287                self.quantile_stats.spatial_cumulated_maxima = all_maxima
 288                self.quantile_stats.spatial_cumulated_maxima_outlets = (
 289                    all_maxima_outlets
 290                )
 291            else:
 292                self.quantile_stats.spatial_cumulated_maxima = np.concat(
 293                    (self.quantile_stats.spatial_cumulated_maxima, all_maxima),
 294                    axis=t_axis,
 295                )
 296                self.quantile_stats.spatial_cumulated_maxima_outlets = np.concat(
 297                    (
 298                        self.quantile_stats.spatial_cumulated_maxima_outlets,
 299                        all_maxima_outlets,
 300                    ),
 301                    axis=1,
 302                )
 303
 304        # setattr(self.quantile_stats, "spatial_maxima", all_maxima)
 305        # setattr(self.quantile_stats, "spatial_maxima_outlets", all_maxima_outlets)
 306        self.quantile_stats.spatial_maxima = all_maxima
 307        self.quantile_stats.spatial_maxima_outlets = all_maxima_outlets
 308
 309    @tools.autocast_args
 310    def fquantile_stats2(
 311        self,
 312        t_axis: int = 2,
 313        return_periods: list | tuple = [2, 5, 10, 20, 50, 100],
 314        fit: str = "gumbel",
 315        nb_minimum_chunks: int = 4,
 316        quantile_duration: list | tuple = [1, 2, 3, 4, 6, 12, 24, 48, 72],
 317        estimate_method: str = "MLE",
 318        chunk_size: int = 365,
 319        ncpu: int | None = None,
 320        # from_maxima: bool = False,
 321        compute_uncertainties: bool = False,
 322        bootstrap_sample: int = 100,
 323    ):
 324        """
 325        Compute the discharge quantile of an 3D array by chunk, corresponding to
 326         `chunk_size` (days), along axis `t_axis` (time) for each pixel of the array.
 327         (nbX, nbY), for each duration `quantile_duration` and for each return period
 328         `return_period`. This function perform the computation in parallel respect
 329         to the the list of quantile duration.
 330
 331        Parameters
 332        ----------
 333
 334        t_axis : int
 335            The axis along with the maximum will be computed. This axis should correspond
 336            to the time.
 337        return_periods: list | tuple
 338            The duration of every return period of unit `chunk_size`.
 339            Default is [2, 5, 10, 20, 50, 100] with a chunk_size=365 days.
 340        fit: str
 341            The extrem law to use to compute the quantile. Choice are
 342            'gumbel' | 'gev'. Default is 'gumbel'.
 343        nb_minimum_chunks: int
 344            number of minimum chunks required to compute the maxima along
 345        the t_axis. Default is set to 4. If the number of chunks is lower, the function
 346        will return None.
 347        estimate_method: str
 348            The method to use to fit rhe Gumbel or GEV law. Choice are `MLE`
 349            (Maximum Likelihood Estimate) or `MM` (Method of Moments). Default is `MLE`.
 350        chunk_size: int
 351            Size of the chunks in days. Default is 365 days.
 352        quantile_duration: list | tuple
 353            The duration of every quantile (hours). The discharges will be resampled for
 354            every duration. Default is [1, 2, 3, 4, 6, 12, 24, 48, 72] hours.
 355        ncpu: int
 356            Number of cpu to use to parrallelize the computation. Default is set to
 357            int(os.cpu_count() / 2).
 358        compute_uncertainties: bool
 359            Compute the uncertainties using the parametric bootstrap method
 360        bootstrap_sample: int
 361            Number of sample using by the bootstrap method, default is 100
 362
 363        Return:
 364        ------
 365        Results are stored in the class spatial_quantile wit different attributes:
 366            - spatial_quantile_matrix : matrix of the spatial quantile for each duration
 367            and each return period.
 368            - spatial_maxima_matrix : matrix of the spatial maxima for each duration and
 369            each return period.
 370            - spatial_cumulated_maxima_matrix : matrix of the spatial accumulated maxima
 371            for each duration and each return period.
 372            - Quantile_{`duration`}h : class of spatial_quantile_results() with attributes:
 373                - self.T : the return periods
 374                - self.Q_th : the quantile matrix for every return periods
 375                - self.T_emp : the empirical return period for each maximum
 376                - self.maxima : the matrix of the maxima
 377                - self.nb_chunks : nb of chunk, i.e data for each pixel
 378                - self.fit : fitting law
 379                - self.fit_shape : matrix of the shape coefficient
 380                - self.fit_scale : matrix of the scale coefficient
 381                - self.fit_loc : matrix of the localisation coefficient
 382                - self.duration : duration of the quantile
 383
 384        Examples
 385        --------
 386        >>> import smashbox
 387        >>> import numpy as np
 388        >>>
 389        >>> graffas_prcp = np.zeros(shape=(142, 166, 1300))
 390        >>> rng = np.random.default_rng()
 391        >>> graffas_prcp = (
 392        >>>     graffas_prcp
 393        >>>     + rng.multinomial(20, [0.2, 0.8], size=graffas_prcp.shape)[:, :, :, 0]
 394        >>> )
 395        >>>
 396        >>> es=smashbox.SmashBox()
 397        >>> sb.newmodel("graffas_zone")
 398        >>> sb.graffas_zone.generate_mesh()
 399        >>> sb.graffas_zone.atmos_data_connector(input_prcp=graffas_prcp)
 400        >>> sb.graffas_zone.model()
 401        >>> sb.graffas_zone.forward_run(invert_states=True, return_options={"q_domain": True})
 402        >>>
 403        >>> sb.graffas_zone.mysmashmodel.mystats.stats_quantile(chunk_size=10, ncpu=6)
 404
 405        """
 406
 407        if pd.Timedelta(
 408            hours=max(quantile_duration),
 409        ) > pd.Timedelta(
 410            days=chunk_size,
 411        ):
 412            raise ValueError(
 413                f"The chunk_size {chunk_size} (days) must be"
 414                f"greater or equal than the quantile duration {max(quantile_duration)} (hours)"
 415            )
 416
 417        if not hasattr(self._parent_class.extra_smash_results, "q_domain"):
 418            raise ValueError(
 419                "No smash extra results  'q_domain' found. Run forward_run()"
 420                "with return_options={'q_domain': True}"
 421            )
 422
 423        if ncpu is None:
 424            ncpu = int(os.cpu_count() / 2)
 425        else:
 426            ncpu = int(min(ncpu, os.cpu_count() - 1))
 427
 428        model_time_step = self._parent_class.smash.setup.dt
 429
 430        shape = list(self._parent_class.extra_smash_results.q_domain.shape)
 431        shape.insert(0, shape.pop(t_axis))
 432
 433        spatial_quantile_matrix = np.zeros(
 434            shape=(
 435                shape[1],
 436                shape[2],
 437                len(quantile_duration),
 438                len(return_periods),
 439            )
 440        )
 441
 442        spatial_quantile_matrix_outlets = np.zeros(
 443            shape=(
 444                len(self._parent_class.smash.mesh.code),
 445                len(quantile_duration),
 446                len(return_periods),
 447            )
 448        )
 449
 450        spatial_maxima = None
 451        spatial_maxima_outlets = None
 452
 453        q_domain = self._parent_class.extra_smash_results.q_domain
 454
 455        partial_sp_quantile = partial(
 456            stats.spatial_quantiles_unparallel,
 457            q_domain,
 458            t_axis,
 459            return_periods,
 460            fit,
 461            nb_minimum_chunks,
 462            model_time_step,
 463            estimate_method,
 464            chunk_size,
 465            compute_uncertainties,
 466            bootstrap_sample,
 467        )
 468
 469        if hasattr(self.quantile_stats, "spatial_cumulated_maxima"):
 470            imap_args = []
 471            for id_dur, duration in enumerate(quantile_duration):
 472                imap_args.append(
 473                    [
 474                        duration,
 475                        self.quantile_stats.spatial_cumulated_maxima[
 476                            :, :, :, id_dur
 477                        ],
 478                    ]
 479                )
 480        else:
 481            imap_args = [
 482                [duration] for id_dur, duration in enumerate(quantile_duration)
 483            ]
 484
 485        with multiprocessing.Pool(ncpu) as p, tqdm(
 486            total=len(quantile_duration)
 487        ) as pbar:
 488
 489            for res in p.starmap(
 490                partial_sp_quantile,
 491                imap_args,
 492                chunksize=1,
 493            ):
 494                pbar.update()
 495                pbar.refresh()
 496                pos = quantile_duration.index(res["duration"])
 497                spatial_quantile_matrix[:, :, pos, :] = res["Q_th"]
 498
 499                for i in range(len(self._parent_class.smash.mesh.code)):
 500                    coords = self._parent_class.smash.mesh.gauge_pos[i]
 501                    spatial_quantile_matrix_outlets[i, pos, :] = (
 502                        spatial_quantile_matrix[coords[0], coords[1], pos, :]
 503                    )
 504
 505                if spatial_maxima is None:
 506                    spatial_maxima = (
 507                        np.zeros(
 508                            shape=(
 509                                *res["maxima"].shape,
 510                                len(quantile_duration),
 511                            )
 512                        )
 513                        + np.nan
 514                    )
 515                    spatial_maxima_outlets = (
 516                        np.zeros(
 517                            shape=(
 518                                len(self._parent_class.smash.mesh.code),
 519                                spatial_quantile["maxima"].shape[t_axis],
 520                                len(quantile_duration),
 521                            )
 522                        )
 523                        + np.nan
 524                    )
 525
 526                spatial_maxima[:, :, :, pos] = res["maxima"]
 527
 528                for i in range(len(self._parent_class.smash.mesh.code)):
 529                    coords = self._parent_class.smash.mesh.gauge_pos[i]
 530                    spatial_maxima_outlets[i, :, pos] = spatial_maxima[
 531                        coords[0], coords[1], :, pos
 532                    ]
 533
 534                if not hasattr(
 535                    self.quantile_stats, f"Quantile_{res['duration']}h"
 536                ):
 537                    setattr(
 538                        self.quantile_stats,
 539                        f"Quantile_{res['duration']}h",
 540                        spatial_quantile_results(),
 541                    )
 542
 543                eval(
 544                    f"self.quantile_stats.Quantile_{res['duration']}h."
 545                    f"fill_attribute(res)"
 546                )
 547
 548                # setattr(self.quantile_stats, "spatial_quantile", spatial_quantile_matrix)
 549                self.quantile_stats.spatial_quantile = spatial_quantile_matrix
 550                self.quantile_stats.spatial_quantile_outlets = (
 551                    spatial_quantile_matrix_outlets
 552                )
 553
 554                # if not from_maxima:
 555                # setattr(self.quantile_stats, "spatial_maxima", spatial_maxima)
 556                self.quantile_stats.spatial_maxima = spatial_maxima
 557                self.quantile_stats.spatial_maxima_outlets = (
 558                    spatial_maxima_outlets
 559                )
 560
 561    @tools.autocast_args
 562    def fquantile_stats(
 563        self,
 564        t_axis: int = 2,
 565        return_periods: list | tuple = [2, 5, 10, 20, 50, 100],
 566        fit: str = "gumbel",
 567        nb_minimum_chunks: int = 4,
 568        quantile_duration: list | tuple = [1, 2, 3, 4, 6, 12, 24, 48, 72],
 569        estimate_method: str = "MLE",
 570        chunk_size: int = 365,
 571        ncpu: int | None = None,
 572        # from_maxima: bool = False,
 573        compute_uncertainties: bool = False,
 574        bootstrap_sample: int = 100,
 575    ):
 576        """
 577        Compute the discharge quantile of an 3D array by chunk, corresponding to
 578        `chunk_size` (days), along axis `t_axis` (time) for each pixel of the array
 579        (nbX, nbY), for each duration `quantile_duration` and for each return period
 580        `return_period`.
 581        This function fit the coefficient of the extrem law in parallel along the Y axis
 582        of the input maxima array (shape=(nbX,nbY,nbchunks)).
 583
 584        Parameters
 585        ----------
 586
 587        t_axis : int
 588            The axis along with the maximum will be computed. This axis should correspond
 589            to the time.
 590        return_periods: list | tuple
 591            The duration of every return period of unit `chunk_size`.
 592            Default is [2, 5, 10, 20, 50, 100] with a chunk_size=365 days.
 593        fit: str
 594            The extrem law to use to compute the quantile. Choice are
 595            'gumbel' | 'gev'. Default is 'gumbel'.
 596        nb_minimum_chunks: int
 597            number of minimum chunks required to compute the maxima along
 598        the t_axis. Default is set to 4. If the number of chunks is lower, the function
 599        will return None.
 600        estimate_method: str
 601            The method to use to fit rhe Gumbel or GEV law. Choice are `MLE`
 602            (Maximum Likelihood Estimate) or `MM` (Method of Moments). Default is `MLE`.
 603        chunk_size: int
 604            Size of the chunks in days. Default is 365 days.
 605        quantile_duration: list | tuple
 606            The duration of every quantile (hours). The discharges will be resampled for
 607            every duration. Default is [1, 2, 3, 4, 6, 12, 24, 48, 72] hours.
 608        ncpu: int
 609            Number of cpu to use to parrallelize the computation. Default is set to
 610            int(os.cpu_count() / 2).
 611        compute_uncertainties: bool
 612            Compute the uncertainties using the parametric bootstrap method
 613        bootstrap_sample: int
 614            Number of sample using by the bootstrap method, default is 100
 615
 616        Return:
 617        ------
 618        Results are stored in the class spatial_quantile wit different attributes:
 619            - spatial_quantile_matrix : matrix of the spatial quantile for each duration
 620            and each return period.
 621            - spatial_maxima_matrix : matrix of the spatial maxima for each duration and
 622            each return period.
 623            - spatial_cumulated_maxima_matrix : matrix of the spatial accumulated maxima
 624            for each duration and each return period.
 625            - Quantile_{`duration`}h : class of spatial_quantile_results() with attributes:
 626                - self.T : the return periods
 627                - self.Q_th : the quantile matrix for every return periods
 628                - self.T_emp : the empirical return period for each maximum
 629                - self.maxima : the matrix of the maxima
 630                - self.nb_chunks : nb of chunk, i.e data for each pixel
 631                - self.fit : fitting law
 632                - self.fit_shape : matrix of the shape coefficient
 633                - self.fit_scale : matrix of the scale coefficient
 634                - self.fit_loc : matrix of the localisation coefficient
 635                - self.duration : duration of the quantile
 636
 637        Examples
 638        --------
 639        >>> import smashbox
 640        >>> import numpy as np
 641        >>>
 642        >>> graffas_prcp = np.zeros(shape=(142, 166, 1300))
 643        >>> rng = np.random.default_rng()
 644        >>> graffas_prcp = (
 645        >>>     graffas_prcp
 646        >>>     + rng.multinomial(20, [0.2, 0.8], size=graffas_prcp.shape)[:, :, :, 0]
 647        >>> )
 648        >>>
 649        >>> es=smashbox.SmashBox()
 650        >>> sb.newmodel("graffas_zone")
 651        >>> sb.graffas_zone.generate_mesh()
 652        >>> sb.graffas_zone.atmos_data_connector(input_prcp=graffas_prcp)
 653        >>> sb.graffas_zone.model()
 654        >>> sb.graffas_zone.forward_run(invert_states=True, return_options={"q_domain": True})
 655        >>>
 656        >>> sb.graffas_zone.mysmashmodel.mystats.stats_quantile(chunk_size=10, ncpu=6)
 657
 658        """
 659
 660        if pd.Timedelta(
 661            hours=max(quantile_duration),
 662        ) > pd.Timedelta(
 663            days=chunk_size,
 664        ):
 665            raise ValueError(
 666                f"The chunk_size {chunk_size} (days) must be"
 667                f" greater or equal than the quantile duration {max(quantile_duration)} (hours)"
 668            )
 669
 670        if not hasattr(self._parent_class.extra_smash_results, "q_domain"):
 671            raise ValueError(
 672                "No smash extra results  'q_domain' found. Run forward_run()"
 673                "with return_options={'q_domain': True}"
 674            )
 675
 676        model_time_step = self._parent_class.smash.setup.dt
 677
 678        shape = list(self._parent_class.extra_smash_results.q_domain.shape)
 679        shape.insert(0, shape.pop(t_axis))
 680
 681        spatial_quantile_matrix = np.zeros(
 682            shape=(
 683                shape[1],
 684                shape[2],
 685                len(quantile_duration),
 686                len(return_periods),
 687            )
 688        )
 689        spatial_quantile_matrix_outlets = np.zeros(
 690            shape=(
 691                len(self._parent_class.smash.mesh.code),
 692                len(quantile_duration),
 693                len(return_periods),
 694            )
 695        )
 696
 697        spatial_maxima = None
 698        spatial_maxima_outlets = None
 699
 700        for id_dur, duration in tqdm(enumerate(quantile_duration)):
 701
 702            print(f"</> Computing spatial quantile for duration {duration}h")
 703
 704            # if hasattr(self.quantile_stats, "spatial_cumulated_maxima"):
 705            if self.quantile_stats.spatial_cumulated_maxima is not None:
 706                maxima = self.quantile_stats.spatial_cumulated_maxima[
 707                    :, :, :, id_dur
 708                ]
 709            else:
 710                maxima = None
 711
 712            spatial_quantile = stats.spatial_quantiles(
 713                array=self._parent_class.extra_smash_results.q_domain,
 714                t_axis=t_axis,
 715                return_periods=return_periods,
 716                fit=fit,
 717                nb_minimum_chunks=nb_minimum_chunks,
 718                model_time_step=model_time_step,
 719                quantile_duration=duration,
 720                estimate_method=estimate_method,
 721                chunk_size=chunk_size,
 722                ncpu=ncpu,
 723                compute_uncertainties=compute_uncertainties,
 724                bootstrap_sample=bootstrap_sample,
 725                maxima=maxima,
 726            )
 727            # else:
 728            #     spatial_quantile = stats.spatial_quantiles(
 729            #         array=self._parent_class.extra_smash_results.q_domain,
 730            #         t_axis=t_axis,
 731            #         return_periods=return_periods,
 732            #         fit=fit,
 733            #         nb_minimum_chunks=nb_minimum_chunks,
 734            #         model_time_step=model_time_step,
 735            #         quantile_duration=duration,
 736            #         estimate_method=estimate_method,
 737            #         chunk_size=chunk_size,
 738            #         ncpu=ncpu,
 739            #         compute_uncertainties=compute_uncertainties,
 740            #         bootstrap_sample=bootstrap_sample,
 741            #         maxima=None,
 742            #     )
 743
 744            print("</>")
 745
 746            pos = quantile_duration.index(spatial_quantile["duration"])
 747            spatial_quantile_matrix[:, :, pos, :] = spatial_quantile["Q_th"]
 748
 749            for i in range(len(self._parent_class.smash.mesh.code)):
 750                coords = self._parent_class.smash.mesh.gauge_pos[i]
 751                spatial_quantile_matrix_outlets[i, :, :] = (
 752                    spatial_quantile_matrix[coords[0], coords[1], :, :]
 753                )
 754
 755            if spatial_maxima is None:
 756                spatial_maxima = (
 757                    np.zeros(
 758                        shape=(
 759                            *spatial_quantile["maxima"].shape,
 760                            len(quantile_duration),
 761                        )
 762                    )
 763                    + np.nan
 764                )
 765                spatial_maxima_outlets = (
 766                    np.zeros(
 767                        shape=(
 768                            len(self._parent_class.smash.mesh.code),
 769                            spatial_quantile["maxima"].shape[t_axis],
 770                            len(quantile_duration),
 771                        )
 772                    )
 773                    + np.nan
 774                )
 775
 776            spatial_maxima[:, :, :, pos] = spatial_quantile["maxima"]
 777
 778            for i in range(len(self._parent_class.smash.mesh.code)):
 779                coords = self._parent_class.smash.mesh.gauge_pos[i]
 780                spatial_maxima_outlets[i, :, id_dur] = spatial_maxima[
 781                    coords[0], coords[1], :, id_dur
 782                ]
 783
 784            if not hasattr(self.quantile_stats, f"Quantile_{duration}h"):
 785                setattr(
 786                    self.quantile_stats,
 787                    f"Quantile_{duration}h",
 788                    spatial_quantile_results(),
 789                )
 790
 791            eval(
 792                f"self.quantile_stats.Quantile_{duration}h."
 793                f"fill_attribute(spatial_quantile)"
 794            )
 795
 796        # setattr(self.quantile_stats, "spatial_quantile", spatial_quantile_matrix)
 797        self.quantile_stats.spatial_quantile = spatial_quantile_matrix
 798        self.quantile_stats.spatial_quantile_outlets = (
 799            spatial_quantile_matrix_outlets
 800        )
 801
 802        # if not from_maxima:
 803        # setattr(self.quantile_stats, "spatial_maxima", spatial_maxima)
 804        self.quantile_stats.spatial_maxima = spatial_maxima
 805        self.quantile_stats.spatial_maxima_outlets = spatial_maxima_outlets
 806
 807
 808class misfit_results:
 809    """
 810    The class misfit_results stores the results of the misfits criterium
 811    """
 812
 813    def __init__(self):
 814        self.mse = None
 815        """MSE for each outlets of the Smash model.
 816         mse = (1.0 / nb_valid_data)* np.sum((obs - sim) ** 2.0)
 817        """
 818        self.rmse = None
 819        """RMSE for each outlets of the Smash model.
 820         rmse = np.sqrt((1.0 / nb_valid_data)* np.sum((obs - sim) ** 2.0))
 821        """
 822        self.nrmse = None
 823        """Normalized-RMSE for each outlets of the Smash model.
 824         nrmse = res_rmse / mean_obs
 825        """
 826        self.se = None
 827        """SE for each outlets of the Smash model. se =
 828            np.sum((obs - sim) ** 2.0)
 829        )"""
 830        self.mae = None
 831        """MAE for each outlets of the Smash model. 
 832            mae = np.sqrt(np.sum(abs(obs - sim))
 833        )"""
 834        self.mape = None
 835        """MAPE for each outlets of the Smash model. 
 836            mape = np.sqrt(
 837                np.sum(abs((obs - sim) / obs))
 838            )
 839        )"""
 840        self.lgrm = None
 841        """LGRM for each outlets of the Smash model. 
 842            lgrm = np.sum(
 843                obs * (np.log((obs / sim) ** 2.0)), axis=t_axis, where=mask_nodata
 844            )
 845        )"""
 846        self.nse = None
 847        """NSE for each outlets of the Smash model."""
 848        self.nnse = None
 849        """NNSE for each outlets of the Smash model. nnse = 1.0 / (2.0 - nse)"""
 850        self.kge = None
 851        """KGE for each outlets of the Smash model."""
 852        self.pearson = None
 853        """Pearson coefficient for each outlets of the Smash model."""
 854
 855
 856class spatial_stats_results:
 857    """
 858    The class spatial_stats_results stores the results of the spatial statistics
 859    on discharges at every pixel.
 860    """
 861
 862    def __init__(self):
 863        self.min = None
 864        """The minimum values for each pixel"""
 865        self.max = None
 866        """The maximum values for each pixel"""
 867        self.mean = None
 868        """The mean value for each pixel"""
 869        self.median = None
 870        """The median values for each pixel"""
 871        self.q20 = None
 872        """The percentile 20% for each pixel"""
 873        self.q80 = None
 874        """The percentile 80% for each pixel"""
 875        self.var = None
 876        """The variance for each pixel"""
 877
 878
 879class outlets_stats_results:
 880    """
 881    The class outlets_stats_results stores the results of the statistics
 882    on discharges at every outlets.
 883    """
 884
 885    def __init__(self):
 886        self.min = None
 887        """The minimum values for each outlets"""
 888        self.max = None
 889        """The maximum values for each outlets"""
 890        self.mean = None
 891        """The mean values for each outlets"""
 892        self.median = None
 893        """The median values for each outlets"""
 894        self.q20 = None
 895        """The percentile 20%  values for each outlets"""
 896        self.q80 = None
 897        """The percentile 80%  values for each outlets"""
 898        self.var = None
 899        """The variance for each pixel"""
 900
 901
 902class spatial_quantile_results:
 903    """
 904    The class spatial_quantile_results stores the results of the quantiles computed for
 905     different return period. Results include the quantile, the empirical return period,
 906     the maxima for each chunk of chunk_size, the fit parameters of the `fit` extrem law.
 907    """
 908
 909    def __init__(self):
 910        self.T = None
 911        """The returns periods"""
 912        self.Q_th = None
 913        """The theorical discharge quantiles for each return period"""
 914        self.T_emp = None
 915        """The empirical return period for each maxima"""
 916        self.maxima = None
 917        """The maxima for each chunk of size chunk_size (default is annual)"""
 918        self.nb_chunks = None
 919        """The number of chunk (default is number of years)"""
 920        self.fit = None
 921        """The extrem law used to fit the maxima and the empirical quantile"""
 922        self.fit_shape = None
 923        """The shape coefficient of the extrem law"""
 924        self.fit_scale = None
 925        """The scale coefficient of the scale law"""
 926        self.fit_loc = None
 927        """The fit coefficient of the extream law"""
 928        self.duration = None
 929        """The duration of the quantile (hours)"""
 930        self.chunk_size = None
 931        """The size of the chunk on which the maxima are computed (unit of the return
 932         period, default is 365 days (1 year))"""
 933        self.Umin = None
 934        """Uncertainties minimum values"""
 935        self.Umax = None
 936        """Uncertainties maximum values"""
 937
 938    def fill_attribute(self, stats_spatial_quantile: dict = None):
 939        """
 940        Fill the attribute of the class spatial_quantile_results.
 941
 942        :param stats_spatial_quantile: Dict of the spatial quantile results, defaults to None
 943        :type stats_spatial_quantile: dict, optional
 944
 945        """
 946
 947        self.T = stats_spatial_quantile["T"]
 948        self.Q_th = stats_spatial_quantile["Q_th"]
 949        self.T_emp = stats_spatial_quantile["T_emp"]
 950        self.maxima = stats_spatial_quantile["maxima"]
 951        self.nb_chunks = stats_spatial_quantile["nb_chunks"]
 952        self.fit = stats_spatial_quantile["fit"]
 953        self.fit_shape = stats_spatial_quantile["fit_shape"]
 954        self.fit_scale = stats_spatial_quantile["fit_scale"]
 955        self.fit_loc = stats_spatial_quantile["fit_loc"]
 956        self.duration = stats_spatial_quantile["duration"]
 957        self.chunk_size = stats_spatial_quantile["chunk_size"]
 958        if "Umin" in stats_spatial_quantile:
 959            self.Umin = stats_spatial_quantile["Umin"]
 960        if "Umax" in stats_spatial_quantile:
 961            self.Umax = stats_spatial_quantile["Umax"]
 962
 963
 964class spatial_quantile:
 965    """Parent class spatial quantile. Class to store results of the spatial quantile."""
 966
 967    def __init__(self):
 968        self.spatial_maxima = None
 969        """Spatial matrix of the maximum discharges computed on period long of `chunk_size`. If chunk_size is 365, the matrix store the maximum annual discharges. Shape=(nbx, nby, nb_chunk, duration)"""
 970        self.spatial_maxima_outlets = None
 971        """Outlets matrix of the maximum discharges computed on period long of `chunk_size`. If chunk_size is 365, the matrix store the maximum annual discharges.Shape=(nb_gauge, nb_chunk, duration)"""
 972        self.spatial_cumulated_maxima = None
 973        """Spatial matrix of the maximum discharges computed on period long of `chunk_size` and accumulated over several simulation. If chunk_size is 365, the matrix store the maximum annual discharges. These maximal values are stacked to the matrix through every simulations.Shape=(nbx, nby, nb_chunk*nb_simulations, duration)"""
 974        self.spatial_cumulated_maxima_outlets = None
 975        """Outlets matrix of the maximum discharges computed on period long of `chunk_size`. If chunk_size is 365, the matrix store the maximum annual discharges. These maximal values are stacked to the matrix through every simulations.Shape=(nb_gauge, nb_chunk*nb_simulations, duration)"""
 976        self.spatial_quantile = None
 977        """Spatial matrix of the discharges quantile computed from the maximum discharges.Shape=(nbx, nby, duration, return_period)"""
 978        self.spatial_quantile_outlets = None
 979        """Outlets matrix of the discharges quantile computed from the maximum discharges.Shape=(nb_gauge, duration, return_period)"""
 980        # pass
 981
 982
 983class spatial_stats:
 984    """Class for computing the spatial statistics on the discharges."""
 985
 986    def __init__(self, parent_class):
 987        self._parent_class = parent_class
 988        """The parent class in order to access to the results of the smash simulation"""
 989        self.results = spatial_stats_results()
 990        """The results of the spatial statistics"""
 991
 992    def mean(self):
 993        """Compute the spatial mean of the discharges."""
 994        self.results.mean = np.nanmean(
 995            self._parent_class._parent_class.extra_smash_results.q_domain,
 996            axis=2,
 997        )
 998
 999    def minimum(self):
1000        """Compute the spatial minimum of the discharges."""
1001        self.results.min = np.nanmin(
1002            self._parent_class._parent_class.extra_smash_results.q_domain,
1003            axis=2,
1004        )
1005
1006    def maximum(self):
1007        """Compute the spatial maximum of the discharges."""
1008        self.results.max = np.nanmax(
1009            self._parent_class._parent_class.extra_smash_results.q_domain,
1010            axis=2,
1011        )
1012
1013    def median(self):
1014        """Compute the spatial median of the discharges."""
1015        self.results.median = np.quantile(
1016            self._parent_class._parent_class.extra_smash_results.q_domain,
1017            0.5,
1018            axis=2,
1019        )
1020
1021    def q20(self):
1022        """Compute the spatial percentile 20% of the discharges."""
1023        self.results.q20 = np.quantile(
1024            self._parent_class._parent_class.extra_smash_results.q_domain,
1025            0.2,
1026            axis=2,
1027        )
1028
1029    def q80(self):
1030        """Compute the spatial percentile 80% of the discharges."""
1031        self.results.q80 = np.quantile(
1032            self._parent_class._parent_class.extra_smash_results.q_domain,
1033            0.8,
1034            axis=2,
1035        )
1036
1037    def var(self):
1038        """Compute the variance of the discharges for each pixel."""
1039        self.results.var = np.var(
1040            self._parent_class._parent_class.extra_smash_results.q_domain,
1041            axis=2,
1042        )
1043
1044
1045class outlets_stats:
1046    """Class for computing the statistics on the discharges at every outlets."""
1047
1048    def __init__(self, parent_class):
1049        self._parent_class = parent_class
1050        """The parent class in order to access to the results of the smash simulation"""
1051        self.results_sim = outlets_stats_results()
1052        """The results of the simulated dicharges statistics at every outlets"""
1053        self.results_obs = outlets_stats_results()
1054        """The results of the observed discharges statistics at every outlets"""
1055
1056    def mean(self):
1057        """Compute the mean of the discharges at every outlets."""
1058        qobs = self._parent_class._parent_class.smash.response_data.q
1059        qobs = np.where(qobs < 0, np.nan, qobs)
1060
1061        self.results_sim.mean = np.nanmean(
1062            self._parent_class._parent_class.smash.response.q, axis=1
1063        )
1064        self.results_obs.mean = np.nanmean(qobs, axis=1)
1065
1066    def minimum(self):
1067        """Compute the minimum of the discharges at every outlets."""
1068        qobs = self._parent_class._parent_class.smash.response_data.q
1069        qobs = np.where(qobs < 0, np.nan, qobs)
1070
1071        self.results_sim.min = np.nanmin(
1072            self._parent_class._parent_class.smash.response.q, axis=1
1073        )
1074        self.results_obs.min = np.nanmin(qobs, axis=1)
1075
1076    def maximum(self):
1077        """Compute the maximum of the discharges at every outlets."""
1078        qobs = self._parent_class._parent_class.smash.response_data.q
1079        qobs = np.where(qobs < 0, np.nan, qobs)
1080
1081        self.results_sim.max = np.nanmax(
1082            self._parent_class._parent_class.smash.response.q, axis=1
1083        )
1084        self.results_obs.max = np.nanmax(qobs, axis=1)
1085
1086    def median(self):
1087        """Compute the median of the discharges at every outlets."""
1088        qobs = self._parent_class._parent_class.smash.response_data.q
1089        qobs = np.where(qobs < 0, np.nan, qobs)
1090
1091        self.results_sim.median = np.nanquantile(
1092            self._parent_class._parent_class.smash.response.q,
1093            0.5,
1094            axis=1,
1095        )
1096        self.results_obs.median = np.nanquantile(qobs, 0.5, axis=1)
1097
1098    def q20(self):
1099        """Compute the percentile 20% of the discharges at every outlets."""
1100        qobs = self._parent_class._parent_class.smash.response_data.q
1101        qobs = np.where(qobs < 0, np.nan, qobs)
1102
1103        self.results_sim.q20 = np.nanquantile(
1104            self._parent_class._parent_class.smash.response.q,
1105            0.2,
1106            axis=1,
1107        )
1108        self.results_obs.q20 = np.nanquantile(qobs, 0.2, axis=1)
1109
1110    def q80(self):
1111        """Compute the percentile 80% of the discharges at every outlets."""
1112        qobs = self._parent_class._parent_class.smash.response_data.q
1113        qobs = np.where(qobs < 0, np.nan, qobs)
1114
1115        self.results_sim.q80 = np.nanquantile(
1116            self._parent_class._parent_class.smash.response.q,
1117            0.8,
1118            axis=1,
1119        )
1120        self.results_obs.q80 = np.nanquantile(qobs, 0.8, axis=1)
1121
1122    def var(self):
1123        """Compute the variance of the discharges at every outlets"""
1124        qobs = self._parent_class._parent_class.smash.response_data.q
1125        qobs = np.where(qobs < 0, np.nan, qobs)
1126
1127        self.results_sim.var = np.var(
1128            self._parent_class._parent_class.smash.response.q, axis=1
1129        )
1130        self.results_obs.var = np.var(qobs, axis=1)
1131
1132
1133class misfit_stats:
1134    """Class for computing the misfit criterium on the discharges at every outlets."""
1135
1136    def __init__(self, parent_class):
1137        self._parent_class = parent_class
1138        """The parent class in order to access to the results of the smash simulation"""
1139        self.results = misfit_results()
1140        """The results of the misfit criterium at every outlets"""
1141
1142    def _update_column(self, column, outlets_name):
1143        if len(outlets_name) > 0:
1144            column = tools.array_isin(
1145                self._parent_class._parent_class.smash.mesh.code,
1146                np.array(outlets_name),
1147            )
1148
1149        if len(column) == 0:
1150            column = list(
1151                range(
1152                    0,
1153                    self._parent_class._parent_class.smash.response.q.shape[0],
1154                )
1155            )
1156
1157        return column
1158
1159    def mse(self, nodata=-99.0, column: list = [], outlets_name: list = []):
1160        """
1161        Compute the mse between the oberved and simulated discharges.
1162        :param nodata: The no data value, defaults to -99.0
1163        :type nodata: float, optional
1164        :param column: the column nuber on which we want to compute the misfit,
1165         defaults to []
1166        :type column: list, optional
1167        :param outlets_name: The names of the outlets for which we want to compute
1168         the misfit, defaults to []
1169        :type outlets_name: list, optional
1170
1171        """
1172
1173        column = self._update_column(column, outlets_name)
1174
1175        self.results.mse = stats.mse(
1176            self._parent_class._parent_class.smash.response_data.q[column, :],
1177            self._parent_class._parent_class.smash.response.q[column, :],
1178            nodata=-99.0,
1179            t_axis=1,
1180        )
1181
1182    def sm_mse(self, column: list = [], outlets_name: list = []):
1183        """
1184        Compute the mse between the oberved and simulated discharges.
1185        :param nodata: The no data value, defaults to -99.0
1186        :type nodata: float, optional
1187        :param column: the column nuber on which we want to compute the misfit,
1188         defaults to []
1189        :type column: list, optional
1190        :param outlets_name: The names of the outlets for which we want to compute
1191         the misfit, defaults to []
1192        :type outlets_name: list, optional
1193
1194        """
1195
1196        column = self._update_column(column, outlets_name)
1197
1198        metric = np.zeros(shape=(len(column))) + np.nan
1199
1200        for i in range(len(column)):
1201            qobs = self._parent_class._parent_class.smash.response_data.q[i, :]
1202            qsim = self._parent_class._parent_class.smash.response.q[i, :]
1203
1204            metric[i] = smash_metrics.mse(
1205                qobs,
1206                qsim,
1207            )
1208
1209        self.results.mse = metric
1210
1211    def rmse(self, nodata=-99.0, column: list = [], outlets_name: list = []):
1212        """
1213        Compute the rmse between the oberved and simulated discharges.
1214        :param nodata: The no data value, defaults to -99.0
1215        :type nodata: float, optional
1216        :param column: the column nuber on which we want to compute the misfit,
1217         defaults to []
1218        :type column: list, optional
1219        :param outlets_name: The names of the outlets for which we want to compute
1220         the misfit, defaults to []
1221        :type outlets_name: list, optional
1222
1223        """
1224
1225        if len(outlets_name) > 0:
1226            column = tools.array_isin(
1227                self._parent_class._parent_class.smash.mesh.code,
1228                np.array(outlets_name),
1229            )
1230
1231        if len(column) == 0:
1232            column = list(
1233                range(
1234                    0,
1235                    self._parent_class._parent_class.smash.response.q.shape[0],
1236                )
1237            )
1238
1239        self.results.rmse = stats.rmse(
1240            self._parent_class._parent_class.smash.response_data.q[column, :],
1241            self._parent_class._parent_class.smash.response.q[column, :],
1242            nodata=-99.0,
1243            t_axis=1,
1244        )
1245
1246    def sm_rmse(self, column: list = [], outlets_name: list = []):
1247        """
1248        Compute the rmse between the oberved and simulated discharges.
1249        :param nodata: The no data value, defaults to -99.0
1250        :type nodata: float, optional
1251        :param column: the column nuber on which we want to compute the misfit,
1252         defaults to []
1253        :type column: list, optional
1254        :param outlets_name: The names of the outlets for which we want to compute
1255         the misfit, defaults to []
1256        :type outlets_name: list, optional
1257
1258        """
1259
1260        column = self._update_column(column, outlets_name)
1261
1262        metric = np.zeros(shape=(len(column))) + np.nan
1263
1264        for i in range(len(column)):
1265            qobs = self._parent_class._parent_class.smash.response_data.q[i, :]
1266            qsim = self._parent_class._parent_class.smash.response.q[i, :]
1267
1268            metric[i] = smash_metrics.rmse(
1269                qobs,
1270                qsim,
1271            )
1272
1273        self.results.rmse = metric
1274
1275    def nrmse(self, nodata=-99.0, column=[], outlets_name=[]):
1276        """
1277        Compute the nrmse between the oberved and simulated discharges.
1278        :param nodata: The no data value, defaults to -99.0
1279        :type nodata: float, optional
1280        :param column: the column nuber on which we want to compute the misfit,
1281         defaults to []
1282        :type column: list, optional
1283        :param outlets_name: The names of the outlets for which we want to compute
1284         the misfit, defaults to []
1285        :type outlets_name: list, optional
1286
1287        """
1288
1289        if len(outlets_name) > 0:
1290            column = tools.array_isin(
1291                self._parent_class._parent_class.smash.mesh.code,
1292                np.array(outlets_name),
1293            )
1294
1295        if len(column) == 0:
1296            column = list(
1297                range(
1298                    0,
1299                    self._parent_class._parent_class.smash.response.q.shape[0],
1300                )
1301            )
1302
1303        self.results.nrmse = stats.nrmse(
1304            self._parent_class._parent_class.smash.response_data.q[column, :],
1305            self._parent_class._parent_class.smash.response.q[column, :],
1306            nodata=-99.0,
1307            t_axis=1,
1308        )
1309
1310    def sm_nrmse(self, column: list = [], outlets_name: list = []):
1311        """
1312        Compute the nrmse between the oberved and simulated discharges.
1313        :param nodata: The no data value, defaults to -99.0
1314        :type nodata: float, optional
1315        :param column: the column nuber on which we want to compute the misfit,
1316         defaults to []
1317        :type column: list, optional
1318        :param outlets_name: The names of the outlets for which we want to compute
1319         the misfit, defaults to []
1320        :type outlets_name: list, optional
1321
1322        """
1323
1324        column = self._update_column(column, outlets_name)
1325
1326        metric = np.zeros(shape=(len(column))) + np.nan
1327
1328        for i in range(len(column)):
1329            qobs = self._parent_class._parent_class.smash.response_data.q[i, :]
1330            qsim = self._parent_class._parent_class.smash.response.q[i, :]
1331
1332            mean_qobs = np.mean(qobs)
1333            metric[i] = smash_metrics.rmse(qobs, qsim) / mean_qobs
1334
1335        self.results.nrmse = metric
1336
1337    def se(self, nodata=-99.0, column=[], outlets_name=[]):
1338        """
1339        Compute the se between the oberved and simulated discharges.
1340        :param nodata: The no data value, defaults to -99.0
1341        :type nodata: float, optional
1342        :param column: the column nuber on which we want to compute the misfit,
1343         defaults to []
1344        :type column: list, optional
1345        :param outlets_name: The names of the outlets for which we want to compute
1346         the misfit, defaults to []
1347        :type outlets_name: list, optional
1348
1349        """
1350        column = self._update_column(column, outlets_name)
1351
1352        self.results.se = stats.se(
1353            self._parent_class._parent_class.smash.response_data.q[column, :],
1354            self._parent_class._parent_class.smash.response.q[column, :],
1355            nodata=-99.0,
1356            t_axis=1,
1357        )
1358
1359    def sm_se(self, column: list = [], outlets_name: list = []):
1360        """
1361        Compute the se between the oberved and simulated discharges.
1362        :param nodata: The no data value, defaults to -99.0
1363        :type nodata: float, optional
1364        :param column: the column nuber on which we want to compute the misfit,
1365         defaults to []
1366        :type column: list, optional
1367        :param outlets_name: The names of the outlets for which we want to compute
1368         the misfit, defaults to []
1369        :type outlets_name: list, optional
1370
1371        """
1372        column = self._update_column(column, outlets_name)
1373        metric = np.zeros(shape=(len(column))) + np.nan
1374
1375        for i in range(len(column)):
1376            qobs = self._parent_class._parent_class.smash.response_data.q[i, :]
1377            qsim = self._parent_class._parent_class.smash.response.q[i, :]
1378
1379            if not np.all(qobs < 0):
1380                metric[i] = smash_metrics.se(
1381                    qobs,
1382                    qsim,
1383                )
1384
1385        self.results.se = metric
1386
1387    def mae(self, nodata=-99.0, column=[], outlets_name=[]):
1388        """
1389        Compute the mae between the oberved and simulated discharges.
1390        :param nodata: The no data value, defaults to -99.0
1391        :type nodata: float, optional
1392        :param column: the column nuber on which we want to compute the misfit,
1393         defaults to []
1394        :type column: list, optional
1395        :param outlets_name: The names of the outlets for which we want to compute
1396         the misfit, defaults to []
1397        :type outlets_name: list, optional
1398
1399        """
1400        column = self._update_column(column, outlets_name)
1401
1402        self.results.mae = stats.mae(
1403            self._parent_class._parent_class.smash.response_data.q[column, :],
1404            self._parent_class._parent_class.smash.response.q[column, :],
1405            nodata=-99.0,
1406            t_axis=1,
1407        )
1408
1409    def sm_mae(self, column=[], outlets_name=[]):
1410        """
1411        Compute the mae between the oberved and simulated discharges.
1412        :param nodata: The no data value, defaults to -99.0
1413        :type nodata: float, optional
1414        :param column: the column nuber on which we want to compute the misfit,
1415         defaults to []
1416        :type column: list, optional
1417        :param outlets_name: The names of the outlets for which we want to compute
1418         the misfit, defaults to []
1419        :type outlets_name: list, optional
1420
1421        """
1422        column = self._update_column(column, outlets_name)
1423        metric = np.zeros(shape=(len(column))) + np.nan
1424
1425        for i in range(len(column)):
1426            qobs = self._parent_class._parent_class.smash.response_data.q[i, :]
1427            qsim = self._parent_class._parent_class.smash.response.q[i, :]
1428
1429            metric[i] = smash_metrics.mae(
1430                qobs,
1431                qsim,
1432            )
1433
1434        self.results.mae = metric
1435
1436    def mape(self, nodata=-99.0, column=[], outlets_name=[]):
1437        """
1438        Compute the mape between the oberved and simulated discharges.
1439        :param nodata: The no data value, defaults to -99.0
1440        :type nodata: float, optional
1441        :param column: the column nuber on which we want to compute the misfit,
1442         defaults to []
1443        :type column: list, optional
1444        :param outlets_name: The names of the outlets for which we want to compute
1445         the misfit, defaults to []
1446        :type outlets_name: list, optional
1447
1448        """
1449        column = self._update_column(column, outlets_name)
1450
1451        self.results.mape = stats.mape(
1452            self._parent_class._parent_class.smash.response_data.q[column, :],
1453            self._parent_class._parent_class.smash.response.q[column, :],
1454            nodata=-99.0,
1455            t_axis=1,
1456        )
1457
1458    def sm_mape(self, column=[], outlets_name=[]):
1459        """
1460        Compute the mape between the oberved and simulated discharges.
1461        :param nodata: The no data value, defaults to -99.0
1462        :type nodata: float, optional
1463        :param column: the column nuber on which we want to compute the misfit,
1464         defaults to []
1465        :type column: list, optional
1466        :param outlets_name: The names of the outlets for which we want to compute
1467         the misfit, defaults to []
1468        :type outlets_name: list, optional
1469
1470        """
1471        column = self._update_column(column, outlets_name)
1472        metric = np.zeros(shape=(len(column))) + np.nan
1473
1474        for i in range(len(column)):
1475            qobs = self._parent_class._parent_class.smash.response_data.q[i, :]
1476            qsim = self._parent_class._parent_class.smash.response.q[i, :]
1477
1478            metric[i] = smash_metrics.mape(
1479                qobs,
1480                qsim,
1481            )
1482
1483        self.results.mape = metric
1484
1485    def lgrm(self, nodata=-99.0, column=[], outlets_name=[]):
1486        """
1487        Compute the lgrm between the oberved and simulated discharges.
1488        :param nodata: The no data value, defaults to -99.0
1489        :type nodata: float, optional
1490        :param column: the column nuber on which we want to compute the misfit,
1491         defaults to []
1492        :type column: list, optional
1493        :param outlets_name: The names of the outlets for which we want to compute
1494         the misfit, defaults to []
1495        :type outlets_name: list, optional
1496
1497        """
1498        column = self._update_column(column, outlets_name)
1499
1500        self.results.lgrm = stats.lgrm(
1501            self._parent_class._parent_class.smash.response_data.q[column, :],
1502            self._parent_class._parent_class.smash.response.q[column, :],
1503            nodata=-99.0,
1504            t_axis=1,
1505        )
1506
1507    def sm_lgrm(self, column=[], outlets_name=[]):
1508        """
1509        Compute the lgrm between the oberved and simulated discharges.
1510        :param nodata: The no data value, defaults to -99.0
1511        :type nodata: float, optional
1512        :param column: the column nuber on which we want to compute the misfit,
1513         defaults to []
1514        :type column: list, optional
1515        :param outlets_name: The names of the outlets for which we want to compute
1516         the misfit, defaults to []
1517        :type outlets_name: list, optional
1518
1519        """
1520        column = self._update_column(column, outlets_name)
1521        metric = np.zeros(shape=(len(column))) + np.nan
1522
1523        for i in range(len(column)):
1524            qobs = self._parent_class._parent_class.smash.response_data.q[i, :]
1525            qsim = self._parent_class._parent_class.smash.response.q[i, :]
1526
1527            if not np.all(qobs < 0):
1528                metric[i] = smash_metrics.lgrm(
1529                    qobs,
1530                    qsim,
1531                )
1532
1533        self.results.lgrm = metric
1534
1535    def nse(self, nodata=-99.0, column=[], outlets_name=[]):
1536        """
1537        Compute the nse between the oberved and simulated discharges.
1538        :param nodata: The no data value, defaults to -99.0
1539        :type nodata: float, optional
1540        :param column: the column nuber on which we want to compute the misfit,
1541         defaults to []
1542        :type column: list, optional
1543        :param outlets_name: The names of the outlets for which we want to compute
1544         the misfit, defaults to []
1545        :type outlets_name: list, optional
1546
1547        """
1548        column = self._update_column(column, outlets_name)
1549
1550        self.results.nse = stats.nse(
1551            self._parent_class._parent_class.smash.response_data.q[column, :],
1552            self._parent_class._parent_class.smash.response.q[column, :],
1553            nodata=-99.0,
1554            t_axis=1,
1555        )
1556
1557    def sm_nse(self, column=[], outlets_name=[]):
1558        """
1559        Compute the nse between the oberved and simulated discharges.
1560        :param nodata: The no data value, defaults to -99.0
1561        :type nodata: float, optional
1562        :param column: the column nuber on which we want to compute the misfit,
1563         defaults to []
1564        :type column: list, optional
1565        :param outlets_name: The names of the outlets for which we want to compute
1566         the misfit, defaults to []
1567        :type outlets_name: list, optional
1568
1569        """
1570        column = self._update_column(column, outlets_name)
1571        metric = np.zeros(shape=(len(column))) + np.nan
1572
1573        for i in range(len(column)):
1574            qobs = self._parent_class._parent_class.smash.response_data.q[i, :]
1575            qsim = self._parent_class._parent_class.smash.response.q[i, :]
1576
1577            metric[i] = smash_metrics.nse(
1578                qobs,
1579                qsim,
1580            )
1581
1582        self.results.nse = metric
1583
1584    def nnse(self, nodata=-99.0, column=[], outlets_name=[]):
1585        """
1586        Compute the nnse between the oberved and simulated discharges.
1587        :param nodata: The no data value, defaults to -99.0
1588        :type nodata: float, optional
1589        :param column: the column nuber on which we want to compute the misfit,
1590         defaults to []
1591        :type column: list, optional
1592        :param outlets_name: The names of the outlets for which we want to compute
1593         the misfit, defaults to []
1594        :type outlets_name: list, optional
1595
1596        """
1597        column = self._update_column(column, outlets_name)
1598
1599        self.results.nnse = stats.nnse(
1600            self._parent_class._parent_class.smash.response_data.q[column, :],
1601            self._parent_class._parent_class.smash.response.q[column, :],
1602            nodata=-99.0,
1603            t_axis=1,
1604        )
1605
1606    def sm_nnse(self, column=[], outlets_name=[]):
1607        """
1608        Compute the nnse between the oberved and simulated discharges.
1609        :param nodata: The no data value, defaults to -99.0
1610        :type nodata: float, optional
1611        :param column: the column nuber on which we want to compute the misfit,
1612         defaults to []
1613        :type column: list, optional
1614        :param outlets_name: The names of the outlets for which we want to compute
1615         the misfit, defaults to []
1616        :type outlets_name: list, optional
1617
1618        """
1619        column = self._update_column(column, outlets_name)
1620        metric = np.zeros(shape=(len(column))) + np.nan
1621
1622        for i in range(len(column)):
1623            qobs = self._parent_class._parent_class.smash.response_data.q[i, :]
1624            qsim = self._parent_class._parent_class.smash.response.q[i, :]
1625
1626            metric[i] = smash_metrics.nnse(
1627                qobs,
1628                qsim,
1629            )
1630
1631        self.results.nnse = metric
1632
1633    def kge(self, nodata=-99.0, column=[], outlets_name=[]):
1634        """
1635        Compute the kge between the oberved and simulated discharges.
1636        :param nodata: The no data value, defaults to -99.0
1637        :type nodata: float, optional
1638        :param column: the column nuber on which we want to compute the misfit,
1639         defaults to []
1640        :type column: list, optional
1641        :param outlets_name: The names of the outlets for which we want to compute
1642         the misfit, defaults to []
1643        :type outlets_name: list, optional
1644
1645        """
1646        column = self._update_column(column, outlets_name)
1647
1648        self.results.kge = stats.kge(
1649            self._parent_class._parent_class.smash.response_data.q[column, :],
1650            self._parent_class._parent_class.smash.response.q[column, :],
1651            nodata=-99.0,
1652            t_axis=1,
1653        )
1654
1655    def sm_kge(self, column=[], outlets_name=[]):
1656        """
1657        Compute the kge between the oberved and simulated discharges.
1658        :param nodata: The no data value, defaults to -99.0
1659        :type nodata: float, optional
1660        :param column: the column nuber on which we want to compute the misfit,
1661         defaults to []
1662        :type column: list, optional
1663        :param outlets_name: The names of the outlets for which we want to compute
1664         the misfit, defaults to []
1665        :type outlets_name: list, optional
1666
1667        """
1668        column = self._update_column(column, outlets_name)
1669        metric = np.zeros(shape=(len(column))) + np.nan
1670
1671        for i in range(len(column)):
1672            qobs = self._parent_class._parent_class.smash.response_data.q[i, :]
1673            qsim = self._parent_class._parent_class.smash.response.q[i, :]
1674
1675            metric[i] = smash_metrics.kge(
1676                qobs,
1677                qsim,
1678            )
1679
1680        self.results.kge = metric
1681
1682    def pearson(self, nodata=-99.0, column=[], outlets_name=[]):
1683        """
1684        Compute the pearson coefficient between the oberved and simulated discharges.
1685        :param nodata: The no data value, defaults to -99.0
1686        :type nodata: float, optional
1687        :param column: the column nuber on which we want to compute the misfit,
1688         defaults to []
1689        :type column: list, optional
1690        :param outlets_name: The names of the outlets for which we want to compute
1691         the misfit, defaults to []
1692        :type outlets_name: list, optional
1693
1694        """
1695        column = self._update_column(column, outlets_name)
1696
1697        self.results.pearson = stats.pearson(
1698            self._parent_class._parent_class.smash.response_data.q[column, :],
1699            self._parent_class._parent_class.smash.response.q[column, :],
1700            nodata=-99.0,
1701            t_axis=1,
1702        )