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