"""
Module containing classical generalistic models

Gumbel:
    To be used applying the Block Maxima approach
    
Generalised extreme value distribution (GEV):
    To be used applying the Block Maxima approach
    
Generalised Pareto Distribution (GPD):
    To be used applying the Peak-Over-Threshold approach
    TODO
"""

from collections import OrderedDict

from scipy import stats as _st
from scipy import optimize as _op
from lmoments3 import distr as _lmdistr
import numpy as _np
import matplotlib.pyplot as _plt
import numdifftools as _ndt

from ..utils import bootstrap_ci as _bsci
from ..utils import gev_momfit as _gev_momfit
from ..utils import gum_momfit as _gum_momfit

class _Base:
    
    def __init__(self, data, fit_method = 'mle', 
                       ci = 0, ci_method = None,
                       return_periods = None, frec = 1):        
        # Data to be used for the fit
        self.data = data
        
        # Fit method to be used
        if fit_method in ['mle', 'mom', 'lmoments']:
            self.fit_method = fit_method
        else:
            raise ValueError(
                ("fit methods accepted are:\n"
                 "    mle (Maximum Likelihood Estimation)\n"
                 "    lmoments\n"
                 "    mom (method of moments)\n")
            )
        
        # Calculate shape, location, scale and a frozen distribution
        # with the calculated estimators (shape, location, scale)
        self._fit()
        
        # Check for calculations of return periods and return values.
        self.frec = frec
        if return_periods:
            self.return_periods = _np.array(return_periods)
            self.return_values = self.distr.isf(self.frec / 
                                                self.return_periods)
        else:
            self.return_periods = _np.array([])
            self.return_values = _np.array([])
        
        # Check for the estimation of confidence intervals
        if ci  == 0 or 0 < ci < 1:
            self.ci = ci
        else:
            raise ValueError("ci should be a value in the interval 0 < ci < 1")
        if self.ci:
            if (ci_method and
                fit_method == 'mle' and 
                ci_method in ['delta', 'bootstrap']):
                self.ci_method = ci_method
                self._ci()
            elif (ci_method and
                fit_method == 'lmoments' and 
                ci_method in ['bootstrap']):
                self.ci_method = ci_method
                self._ci()
            elif (ci_method and
                fit_method == 'mom' and 
                ci_method in ['bootstrap']):
                self.ci_method = ci_method
                self._ci()
            else:
                raise ValueError(
                ("You should provide a valid value for the confidence\n"
                 "interval calculation, 'ci_method'\n"))
         
    
    def _fit(self):
        # This is a base class and shouldn't be used as it. 
        # This method should be implemented in the subclass.
        raise NotImplementedError("Subclasses should implement this!")
    
    def _ci(self):
        # This is a base class and shouldn't be used as it. 
        # This method should be implemented in the subclass.
        raise NotImplementedError("Subclasses should implement this!")
    
    def pdf(self, quantiles):
        # A shortcut to the frozen distribution pdf as provided by scipy.
        """
        Probability density function at x of the given frozen RV.
        
        **Parameters**
        
        x : array_like
            quantiles
            
        **Returns**
        
        pdf : ndarray
            Probability density function evaluated at x
        """
        
        return self.distr.pdf(quantiles)
    
    def cdf(self, quantiles):
        # A shortcut to the frozen distribution cdf as provided by scipy.
        """
        Cumulative distribution function of the given frozen RV.
        
        **Parameters**
        
        x : array_like
            quantiles

        **Returns**
        
        cdf : ndarray
            Cumulative distribution function evaluated at `x`
        """
        
        return self.distr.cdf(quantiles)
        
    def ppf(self, q):
        # A shortcut to the frozen distribution ppf as provided by scipy.
        """
        Percent point function (inverse of cdf) at q of the given frozen RV.

        **Parameters**
        
        q : array_like
            lower tail probability
        
        **Returns**
        
        x : array_like
            quantile corresponding to the lower tail probability q.
        """
        
        return self.distr.ppf(q)
    
    def stats(self, moments):
        # A shortcut to the frozen distribution stats as provided by scipy.
        """
        Some statistics of the given RV.

        **Parameters**
        
        moments : str, optional
            composed of letters ['mvsk'] defining which moments to compute:
            'm' = mean,
            'v' = variance,
            's' = (Fisher's) skew,
            'k' = (Fisher's) kurtosis.
            (default='mv')

        **Returns**
        
        stats : sequence
            of requested moments.
        """
        
        return self.distr.stats(moments)    
    
    def _plot(self, ax, title, xlabel, ylabel):
        # helper function for:
        #     self.plot_density()
        #     self.plot_pp()
        #     self.plot_qq()
        #     self.plot_return_values()
        #     self.plot_summary()
        ax.set_facecolor((0.95, 0.95, 0.95))
        _plt.setp(ax.lines, linewidth = 2, color = 'magenta')
        ax.set_title(title)
        ax.set_xlabel(xlabel)
        ax.set_ylabel(ylabel)
        ax.grid(True)
        return ax
        
    def plot_density(self):
        """
        Histogram of the empirical pdf data and the pdf plot of the 
        fitted distribution.
        All parameters are predefined from the frozen fitted model and empirical
        data available.

        **Returns**
        
        Density plot.
        """
        
        fig, ax = _plt.subplots(figsize=(8, 6))
        
        # data
        x = _np.linspace(self.distr.ppf(0.001), self.distr.ppf(0.999), 100)
        
        # plot
        ax.plot(x, self.distr.pdf(x), label = 'Fitted', color = 'k')
        ax.hist(self.data, density = True, 
                color = 'yellow', alpha = 0.75, label = "Empirical")
        ax = self._plot(ax, 'Density Plot', 'x', 'f(x)')
        ax.legend(loc='best', frameon=False)     
    
    def plot_pp(self):
        """
        PP (probability) plot between empirical and fitted data.
        All parameters are predefined from the frozen fitted model and empirical
        data available.
        
        **Returns**
        
        PP plot. 
        """
        
        fig, ax = _plt.subplots(figsize=(8, 6))
        
        # data
        data = _np.sort(self.data)
        N = len(data)
        y = _np.arange(1, N + 1) / (N + 1)
        x = self.distr.cdf(data)
        
        # plot
        ax.scatter(x, y, color = 'darkcyan')
        ax.plot([0, 1], [0, 1])
        ax.set_xlim(0, 1)
        ax.set_ylim(0, 1)
        ax = self._plot(ax, 'P-P Plot', 'model', 'empirical')
        
    def plot_qq(self):
        """
        QQ (Quantile-Quantile) plot between empirical and fitted data.
        All parameters are predefined from the frozen fitted model and empirical
        data available.
        
        **Returns**
        
        QQ plot. 
        """
        
        fig, ax = _plt.subplots(figsize=(8, 6))
        
        # data
        y = _np.sort(self.data)
        N = len(y)
        x = _np.arange(1, N + 1) / (N + 1)
        x = self.distr.ppf(x)
        
        # plot
        ax = self._plot(ax, 'Q-Q Plot', 'model', 'empirical')
        ax.scatter(x, y, color = 'forestgreen')
        low_lim = _np.min([x, y]) * 0.95
        high_lim = _np.max([x, y]) * 1.05
        ax.plot([low_lim, high_lim], [low_lim, high_lim], c='k')
        ax.set_xlim(low_lim, high_lim)
        ax.set_ylim(low_lim, high_lim)
        
    def plot_return_values(self):
        """
        Return values and return periods of data. If confidence interval 
        information has been provided it will show the confidence interval 
        values.
        
        **Returns**
        
        Return values and return periods plot. 
        """
        
        fig, ax = _plt.subplots(figsize=(8, 6))
        
        # data
        T = _np.arange(0.1, 500.1, 0.1)
        sT = self.distr.isf(self.frec * 1./T)
        N = _np.r_[1:len(self.data)+1] * self.frec
        Nmax = max(N)
        
        # plot
        ax = self._plot(ax, 'Return Level Plot', 'Return period', 'Return level')
        ax.semilogx(T, sT)
        ax.scatter(self.frec * Nmax/N, sorted(self.data)[::-1], color = 'orangered')
        
        # plot confidence intervals if available
        if self.ci:
            #y1 = sT - st.norm.ppf(1 - self.ci / 2) * np.sqrt(self._ci_se)
            #y2 = sT + st.norm.ppf(1 - self.ci / 2) * np.sqrt(self._ci_se)
            ax.semilogx(T, self._ci_Td, '--')
            ax.semilogx(T, self._ci_Tu, '--')
            ax.fill_between(T, self._ci_Td, self._ci_Tu, color = '0.75', alpha = 0.5)
        
    def plot_summary(self):
        """
        Summary plot including PP plot, QQ plot, empirical and fitted pdf and
        return values and periods.
        
        **Returns**
        
        4-panel plot including PP, QQ, pdf and return level plots
        """
        
        fig, ((ax3, ax2), (ax4, ax1)) = _plt.subplots(2, 2, figsize=(8, 6))
        
        # PDF plot
        x = _np.linspace(self.distr.ppf(0.001), 
                        self.distr.ppf(0.999), 
                        100)
        ax1.plot(x, self.distr.pdf(x), label = 'Fitted')
        ax1.hist(self.data, density = True, 
                color = 'yellow', alpha = 0.75, label = "Empirical")
        ax1 = self._plot(ax1, 'Density Plot', 'x', 'f(x)')
        ax1.legend(loc='best', frameon=False)
        
        # QQ plot
        data = _np.sort(self.data)
        N = len(data)
        y = _np.arange(1, N + 1) / (N + 1)
        x = self.distr.cdf(data)
        ax2.plot([0, 1], [0, 1])
        ax2.set_xlim(0, 1)
        ax2.set_ylim(0, 1)
        ax2 = self._plot(ax2, 'P-P Plot', 'model', 'empirical')
        ax2.scatter(x, y, color = 'darkcyan')
        
        # PP Plot
        y = _np.sort(self.data)
        N = len(y)
        x = _np.arange(1, N + 1) / (N + 1)
        x = self.distr.ppf(x)
        ax3.scatter(x, y, color = 'forestgreen')
        low_lim = _np.min([x, y]) * 0.95
        high_lim = _np.max([x, y]) * 1.05
        ax3.plot([low_lim, high_lim], [low_lim, high_lim])
        ax3.set_xlim(low_lim, high_lim)
        ax3.set_ylim(low_lim, high_lim)
        ax3 = self._plot(ax3, 'Q-Q Plot', 'model', 'empirical')
        
        # Return levels plot
        T = _np.arange(0.1, 500.1, 0.1)
        sT = self.distr.isf(self.frec/T)
        N = _np.r_[1:len(self.data)+1] * self.frec
        Nmax=max(N)
        ax4 = self._plot(ax4, 'Return Level Plot', 
                              'Return period', 
                              'Return level')
        ax4.semilogx(T, sT, 'k')
        ax4.scatter(self.frec * Nmax/N, sorted(self.data)[::-1], color = 'orangered')
        if self.ci:
            #y1 = sT - st.norm.ppf(1 - self.ci / 2) * np.sqrt(self._ci_se)
            #y2 = sT + st.norm.ppf(1 - self.ci / 2) * np.sqrt(self._ci_se)
            ax4.semilogx(T, self._ci_Td, '--')
            ax4.semilogx(T, self._ci_Tu, '--')
            ax4.fill_between(T, self._ci_Td, self._ci_Tu, color = '0.75', alpha = 0.5)
        
        # I love matplotlib for stuff like this, thanks, guys!!!
        _plt.tight_layout()
        return fig, ax1, ax2, ax3, ax4
        
        
class GEV(_Base):
    """
    Class to fit data to a Generalised extreme value (GEV) distribution.
    
    **Parameters**
        
    data : array_like
        1D array_like with the extreme values to be considered
    fit_method : str
        String indicating the method used to fit the distribution.
        Availalable values are 'mle' (default value), 'mom' and 'lmoments'.
    ci : float (optional)
        Float indicating the value to be used for the calculation of the 
        confidence interval. The returned values are (ci/2, 1-ci/2) 
        percentile confidence intervals. E.g., a value of 0.05 will 
        return confidence intervals at 0.025 and 0.975 percentiles.
    ci_method : str (optional)
        String indicating the method to be used to calculate the 
        confidence intervals. If ``ci`` is not supplied this parameter will 
        be ignored. Possible values depend of the fit method chosen. If 
        the fit method is 'mle' possible values for ci_method are 
        'delta' and 'bootstrap', if the fit method is 'mom' or 
        'lmoments' possible value for ci_method is 'bootstrap'.
            'delta' is for delta method.
            'bootstrap' is for parametric bootstrap.
    return_period : array_like (optional)
        1D array_like of values for the *return period*. Values indicate
        **years**. 
    frec : int or float
        Value indicating the frecuency of events per year. If frec is 
        not provided the data will be treated as yearly data (1 value per 
        year).

    **Attributes and Methods**
    
    params : OrderedDict
        Ordered dictionary with the values of the *shape*, *location* and
        *scale* parameters of the distribution.
    c : flt
        Float value for the *shape* parameter of the distribution.
    loc : flt
        Float value for the *location* parameter of the distribution.
    scale : flt
        Float value for the *scale* parameter of the distribution.
    distr : object
        Frozen RV object with the same methods of a continuous scipy
        distribution but holding the given *shape*, *location*, and *scale* 
        fixed. See http://docs.scipy.org/doc/scipy/reference/stats.html
        for more info.
    data : array_like
        Input data used for the fit
    fit_method : str
        String indicating the method used to fit the distribution,
        values can be 'mle', 'mom' or 'lmoments'.
    """
    
    def _fit(self):
        # Fit can be made using Maximum Likelihood Estimation (mle) or using
        # l-moments. 
        # L-moments is fast and accurate most of the time for the GEV
        # distribution. 
        
        # MLE FIT
        # In the case of the mle estimation, sometimes we get unstable values
        # if we don't provide an initial guess of the parameters. Loc and scale 
        # are more or less stable but shape can be quite unstable depending the 
        # input data. This is why we are using lmoments to obtain start values 
        # for the mle optimization. For mle we are using fmin_bfgs as it is 
        # faster than others and with the first guess provide accurate results.
        if self.fit_method == 'mle':
            # Initial guess to make the fit of GEV more stable
            # To do the initial guess we are using lmoments...
            _params0 = _lmdistr.gev.lmom_fit(self.data)
            # The mle fit will start with the initial estimators obtained
            # with lmoments above
            _params = _st.genextreme.fit(self.data, _params0['c'],
                                        loc = _params0['loc'],
                                        scale = _params0['scale'],
                                        optimizer = _op.fmin_bfgs)
            self.params = OrderedDict()
            # For the shape parameter the value provided by scipy 
            # is defined as negative as that obtained from other
            # packages in R, some textbooks, wikipedia,... ¿?
            self.params["shape"]    = _params[0]
            self.params["location"] = _params[1]
            self.params["scale"]    = _params[2]
        
        # L-MOMENTS FIT
        if self.fit_method == 'lmoments':
            _params = _lmdistr.gev.lmom_fit(self.data)
            self.params = OrderedDict()
            # For the shape parameter the value provided by lmoments3 
            # is defined as negative as that obtained from other
            # packages in R, some textbooks, wikipedia,... ¿?
            self.params["shape"]    = _params['c']
            self.params["location"] = _params['loc']
            self.params["scale"]    = _params['scale']
        
        # METHOD OF MOMENTS FIT
        if self.fit_method == 'mom':
            _params = _gev_momfit(self.data)
            self.params = OrderedDict()
            self.params["shape"]    = _params[0]
            self.params["location"] = _params[1]
            self.params["scale"]    = _params[2]
        
        # Estimators and a frozen distribution for the estimators
        self.c     = self.params['shape']      # shape
        self.loc   = self.params['location']   # location
        self.scale = self.params['scale']      # scale
        self.distr = _st.genextreme(self.c,     # frozen distribution
                                   loc = self.loc, 
                                   scale = self.scale)
                                   
    def _nnlf(self, theta):
        # This is used to calculate the variance-covariance matrix using the
        # Hessian from numdifftools
        # see self._ci_delta() method below
        
        x = self.data
        
        # Here we provide code for the GEV distribution and for the special
        # case when shape parameter is 0 (Gumbel distribution).
        if len(theta) == 3:
            c = theta[0]
            loc = theta[1]
            scale = theta[2]
        if len(theta) == 2:
            c = 0
            loc = theta[0]
            scale = theta[1]
        if c != 0:
            expr = 1. + c * ((x - loc) / scale)
            return (len(x) * _np.log(scale) + 
                   (1. + 1. / c) * _np.sum(_np.log(expr)) +
                   _np.sum(expr ** (-1. / c)))
        else:
            expr = (x - loc) / scale
            return (len(x) * _np.log(scale) + 
                   _np.sum(expr) +
                   _np.sum(_np.exp(-expr)))
     
    def _ci_delta(self):
        # Calculate the variance-covariance matrix using the
        # hessian from numdifftools
        # This is used to obtain confidence intervals for the estimators and
        # the return values for several return values.
        #
        # More info about the delta method can be found on:
        #     - Coles, Stuart: "An Introduction to Statistical Modeling of  
        #     Extreme Values", Springer (2001)
        #     - https://en.wikipedia.org/wiki/Delta_method
        
        # data
        c  = -self.c    # We negate the shape to avoid inconsistency problems!?
        loc = self.loc
        scale = self.scale
        hess = _ndt.Hessian(self._nnlf)
        T = _np.arange(0.1, 500.1, 0.1)
        sT = -_np.log(1.-self.frec/T)
        sT2 = self.distr.isf(self.frec/T)
        
        # VarCovar matrix and confidence values for estimators and return values
        # Confidence interval for return values (up values and down values)
        ci_Tu = _np.zeros(sT.shape)
        ci_Td = _np.zeros(sT.shape)
        if c:         # If c then we are calculating GEV confidence intervals
            varcovar = _np.linalg.inv(hess([c, loc, scale]))
            self.params_ci = OrderedDict()
            se = _np.sqrt(_np.diag(varcovar))
            self._se = se
            self.params_ci['shape']    = (self.c - _st.norm.ppf(1 - self.ci / 2) * se[0],
                                          self.c + _st.norm.ppf(1 - self.ci / 2) * se[0])
            self.params_ci['location'] = (self.loc - _st.norm.ppf(1 - self.ci / 2) * se[1],
                                          self.loc + _st.norm.ppf(1 - self.ci / 2) * se[1])
            self.params_ci['scale']    = (self.scale - _st.norm.ppf(1 - self.ci / 2) * se[2],
                                          self.scale + _st.norm.ppf(1 - self.ci / 2) * se[2])
            for i, val in enumerate(sT2):
                gradZ = [scale * (c**-2) * (1 - sT[i] ** (-c)) - scale * (c**-1) * (sT[i]**-c) * _np.log(sT[i]),
                         1, 
                         -(1 - sT[i] ** (-c)) / c]
                se = _np.dot(_np.dot(gradZ, varcovar), _np.array(gradZ).T)
                ci_Tu[i] = val + _st.norm.ppf(1 - self.ci / 2) * _np.sqrt(se)
                ci_Td[i] = val - _st.norm.ppf(1 - self.ci / 2) * _np.sqrt(se)
        else:         # else then we are calculating Gumbel confidence intervals
            varcovar = _np.linalg.inv(hess([loc, scale]))
            self.params_ci = OrderedDict()
            se = _np.sqrt(_np.diag(varcovar))
            self._se = se
            self.params_ci['shape']    = (0, 0)
            self.params_ci['location'] = (self.loc - _st.norm.ppf(1 - self.ci / 2) * se[0],
                                          self.loc + _st.norm.ppf(1 - self.ci / 2) * se[0])
            self.params_ci['scale']    = (self.scale - _st.norm.ppf(1 - self.ci / 2) * se[1],
                                          self.scale + _st.norm.ppf(1 - self.ci / 2) * se[1])
            for i, val in enumerate(sT2):
                gradZ = [1, -_np.log(sT[i])]
                se = _np.dot(_np.dot(gradZ, varcovar), _np.array(gradZ).T)
                ci_Tu[i] = val + _st.norm.ppf(1 - self.ci / 2) * _np.sqrt(se)
                ci_Td[i] = val - _st.norm.ppf(1 - self.ci / 2) * _np.sqrt(se)
        self._ci_Tu = ci_Tu
        self._ci_Td = ci_Td
    
    def _ci_bootstrap(self):
        # Calculate confidence intervals using parametric bootstrap and the
        # percentil interval method
        # This is used to obtain confidence intervals for the estimators and
        # the return values for several return values.
        # all the code in skextremes.utils.bootstrap_ci has been adapted and
        # simplified from that on https://github.com/cgevans/scikits-bootstrap.
        #
        # More info about bootstrapping can be found on:
        #     - https://github.com/cgevans/scikits-bootstrap
        #     - Efron: "An Introduction to the Bootstrap", Chapman & Hall (1993)
        #     - https://en.wikipedia.org/wiki/Bootstrapping_%28statistics%29
        
        # parametric bootstrap for return levels and parameters   
        
        # The function to bootstrap     
        def func(data):
            sample = _st.genextreme.rvs(self.c, 
                                       loc = self.loc, 
                                       scale = self.scale, 
                                       size = len(self.data))
            c, loc, scale = _st.genextreme.fit(sample, self.c, 
                                              loc = self.loc, 
                                              scale = self.scale,
                                              optimizer = _op.fmin_bfgs)
            T = _np.arange(0.1, 500.1, 0.1)
            sT = _st.genextreme.isf(self.frec/T, c, loc = loc, scale = scale)
            res = [c, loc, scale]
            res.extend(sT.tolist())
            return tuple(res)
        
        # the calculations itself
        out = _bsci(self.data, statfunction = func, n_samples = 500)
        self._ci_Td = out[0, 3:]
        self._ci_Tu = out[1, 3:]
        self.params_ci = OrderedDict()
        self.params_ci['shape']    = (out[0,0], out[1,0])
        self.params_ci['location'] = (out[0,1], out[1,1])
        self.params_ci['scale']    = (out[0,2], out[1,3])
    
    def _ci(self):
        # Method called internally to calculate confidence intervals if 
        # required. To see more info about available methods see comments on
        # self._ci_delta and self._ci_bootstrap methods.
        
        if self.ci_method == "delta":
            self._ci_delta()
        if self.ci_method == "bootstrap":
            self._ci_bootstrap()

class Gumbel(GEV):
    __doc__ = GEV.__doc__.replace("Generalised extreme value (GEV) distribution.", 
                                  ("Gumbel distribution. Note that this is a "
                                   "special case of the ``GEV`` class where "
                                   "the 'shape' is fixed to 0."))
    
    def _fit(self):
        
        if self.fit_method == 'mle':
            _params = _st.gumbel_r.fit(self.data)
            self.params = OrderedDict()
            self.params["shape"]    = 0
            self.params["location"] = _params[0]
            self.params["scale"]    = _params[1]
        
        if self.fit_method == 'lmoments':
            _params = _lmdistr.gum.lmom_fit(self.data)
            self.params = OrderedDict()
            self.params["shape"]    = 0
            self.params["location"] = _params['loc']
            self.params["scale"]    = _params['scale']
        
        # METHOD OF MOMENTS FIT
        if self.fit_method == 'mom':
            _params = _gum_momfit(self.data)
            self.params = OrderedDict()
            self.params["shape"]    = _params[0]
            self.params["location"] = _params[1]
            self.params["scale"]    = _params[2]
        
        self.c     = self.params['shape']
        self.loc   = self.params['location']
        self.scale = self.params['scale']
        self.distr = _st.gumbel_r(loc = self.loc, 
                                   scale = self.scale)

class GPD(_Base):
    pass
