#                                               -*- Python -*-
#
#  @file  polynomial_qsa.py
#  @brief Validation of sensitivity evaluation for quantiles with the Polynomial model
#
#  Copyright (C) 2017 Airbus-IMACS
#
#  Written by Sofiane Haddad, haddad@imacs.polytechnique.fr
#             Nabil Rachdi, nabil.rachdi@airbus.com
#
#  This program is free software; you can redistribute it and/or
#  it under the terms of the GNU Lesser General Public License as published by
#  the Free Software Foundation, either version 3 of the License, or
#  (at your option) any later version.
#
#  This program is distributed in the hope that it will be useful
#  but WITHOUT ANY WARRANTY; without even the implied warranty of
#  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#  GNU Lesser General Public License for more details.
#
#  You should have received a copy of the GNU Lesser General Public
#  along with this library.  If not, see <http://www.gnu.org/licenses/>.
#

"""
    polynomial_qsa
    ==============
    Part of CHORUS project

"""
import openturns as ot
from math import sin, pi
import numpy as np
import matplotlib
matplotlib.style.use('ggplot')
matplotlib.rcParams["xtick.labelsize"] = "small"
import matplotlib.pylab as plt
from pygosa import plot_quantiles_sensitivities
from pygosa.tests import sensitivity_quantiles_evaluation
from pygosa.tests import polynomial_uc



# Size
size = 3000
qLevels = [1.0/size, 0.01, 0.05, 0.10, 0.15, 0.20, 0.25, 0.50, 0.75, 0.80, 0.85, 0.90, 0.95, 0.99, 1 - 1.0/size]
# UC
distribution_polynomial = polynomial_uc.define_polynomial_distribution()

s1 = []
s2 = []
Alpha = range(1,10)
for alpha in Alpha:
    polynomial_model = polynomial_uc.define_polynomial_function(alpha)

    #  Compute sensitivities using quantile contrast
    print("alpha=", alpha)
    ot.RandomGenerator.SetSeed(0)
    quantile_factors, qs =  sensitivity_quantiles_evaluation(qLevels, dist=distribution_polynomial, model=polynomial_model, size=size)
    # Add quantile factors to S1/S2
    s1.append(quantile_factors[:,0])
    s2.append(quantile_factors[:,1])

s1 = np.array(s1).T
s2 = np.array(s2).T

fig = plt.figure()
plt.plot(Alpha, s1[8, :])
plt.plot(Alpha, s2[8, :])
plt.plot(Alpha, s1[1, :])
plt.plot(Alpha, s2[1, :])
plt.savefig("sensitivity_f_alpha.png")
plt.close()
"""
# QAD ~ Quantile Absolute Deviation might writes k * sigma
A = np.sqrt( quantile_factors[:,0] ) + np.sqrt( quantile_factors[:,1] )
print(np.sqrt( quantile_factors[:,0] ) / A)
print(np.sqrt( quantile_factors[:,1] ) / A)

from pygosa import QuantileSensitivities
v = np.array( [np.mean( QuantileSensitivities._contrast_quantile(np.array(qs.mcsp.y_), qs.mcsp.y_.computeQuantile(q)[0], q) ) for q in qLevels] )
print(v)

from statsmodels.regression.quantile_regression import QuantReg
mod = QuantReg(y.ravel(),x)
mod = mod.fit(q=0.5)
print(mod.summary())
"""
