import os
import numpy as np
import pandas as pd

import matplotlib as mpl
import matplotlib.pyplot as plt

plt.ioff()
mpl.rcParams.update({'font.size': 18})
import statsmodels.api as sm
from sklearn.linear_model import LinearRegression
from statsmodels.sandbox.regression.predstd import wls_prediction_std
from scipy import odr as odr
from scipy import stats as scistats

import RTxploitation as rt
iopw = rt.auxdata.iopw().get_iopw

opj= os.path.join
idir = os.path.abspath('/DATA/OBS2CO/data/rogerio')
figdir = opj(idir,'fig')
file = opj(idir,'data/Simulation_Rrs_OSOAA_TSS_COD_aCDOM.xlsx')

data = pd.read_excel(file)
data.sort_values(by="TSS (mg L)", inplace=True)
data = data.set_index(['Station', 'Date'])

acdom = data.iloc[:, 17]
acdom_sd = data.iloc[:, 18]
#remove sparse data of CDOM
data = data.iloc[:,:-2]
params=['SPM','DOC']
param=params[0]
if param == 'DOC':
    data.dropna(inplace=True)
month = data.index.get_level_values(1).month
spm = data.iloc[:, 0]  # .values.reshape(-1, 1)
spm_sd = data.iloc[:, 1]
doc = data.iloc[:, 15]
doc_sd = data.iloc[:, 16]

def Rrs2rrs(Rrs):
        return Rrs / (0.52 + 1.7 * Rrs)

def inv_gordon88(rrs,g0=0.089,g1=0.125):
        deltas = g0**2 + 4 * g1 * rrs
        return (-g0 + np.sqrt(deltas)) / (2*g1)

def black_water_model_old(B, x , wl=550):
    aw,bbw = iopw(wl)
    rrs = Rrs2rrs(x)
    u = inv_gordon88(rrs)
    return (B[0])**2 * (u - np.abs(B[1])  )

def black_water_model(B, x, wl=550):

    rrs = Rrs2rrs(x)
    u = inv_gordon88(rrs)
    aw,bbw = iopw(wl)
    N = (u*(aw+bbw+B[0])-bbw)/ (B[1]- u *(B[1]+B[2]))

    return B[3] * N

def linear(B, x):
    '''Linear function y = m*x + b'''
    return B[0] + B[1] * x

def exponential(B, x):
    #return B[0] - (B[0]-B[1]) * np.exp(-B[2] * x)
    return B[0] * np.exp(B[1] * x) + B[2]

def hyperbolic(B,x,b=1):
    return B[1]/(1+B[0]*x)**(1/b)

def poly(B, x):
    return B[0] + B[1] * x + B[2] * x ** 2  # + B[3]*x**3

def confidence_interval(xn, model, res, nstd=1):
    '''

    :param xn: x-axis data for computation
    :param model: numerical model used for the fit
    :param res: output from scipy.odr.run
    :param nstd: number of sigma to compute confidence interval
    :return: data up and data down
    '''
    '''
    
    :param res: output from scipy.odr.run
    :param nstd: number of sigma to compute confidence interval
    :return: 
    '''

    popt_up = res.beta + nstd * res.sd_beta
    popt_dw = res.beta - nstd * res.sd_beta

    return model(popt_up, xn), model(popt_dw, xn)

if param == 'SPM':
    y=spm
    y_sd=spm_sd
    fit_eq='$g_0={:.1f};\ g_1={:.1f};\ C={:.1f}$'
    ylabel='SPM (mg/L)'
    # gordon values
    g0 = 0.089
    g1 = 0.125
    model,beta0=black_water_model_old,[150,10]
    #model,beta0=black_water_model,[1,1,1,10]
else:
    y=doc
    y_sd=doc_sd
    fit_eq='y = {:.1f}/(1+{:.1f}$x)$\n res_var = {:.1f}'
    ylabel='DOC (mg/L)'
    model,beta0=hyperbolic,[10,10]

c=doc
cmap = plt.cm.get_cmap("Spectral").reversed()
stats = pd.DataFrame(columns=['band','algo','b0','b1','b2','sig_b0','sig_b1','sig_b2'])


# plot Rrs = f(SPM)
wls = [560, 665, 704,740, 782, 864]
fig, axs = plt.subplots(nrows=2, ncols=3, figsize=(20, 12))
for i, ax in enumerate(axs.ravel()):
    wl=wls[i]

    Rrs = data.iloc[:, 2 * i + 2]
    print(Rrs.name,wl)
    Rrs_sd = data.iloc[:, 2 * i + 3]

    x,x_sd = Rrs, Rrs_sd

    ax.errorbar(x, y, xerr=x_sd, yerr=y_sd, fmt='o', color='grey', ecolor='black', alpha=0.8)
    # to put color in scatter dots
    #ax.scatter(x, y, c=c, s=55, cmap=cmap)
    ax.set_title(x.name)
    ax.set_xlabel(r'$R_{rs}\ (sr^{-1})$')
    ax.set_ylabel(ylabel)

    testdata = odr.RealData(x, y, sx=x_sd, sy=y_sd)
    _odr = odr.ODR(testdata, odr.Model(model, extra_args=[wl]), beta0=beta0) #

    # fit with ordinary least square (OLS, fit_type=2 )
    _odr.set_job(fit_type=2)
    res = _odr.run()

    res.pprint()
    xn = np.linspace(0, np.max(x) * 1.25, 500)
    yn = model(res.beta, xn)
    fit_up, fit_dw = confidence_interval(xn, model, res)

    #stats.loc[2*i]=np.concatenate([[Rrs.name, 'OLS'], res.beta, res.sd_beta])


    ax.plot(xn, yn, 'b--', label='OLS',#; '+fit_eq.format(*res.beta),
            linewidth=2)
    ax.fill_between(xn, fit_up, fit_dw, alpha=.25, facecolor="b")  # ,label="1-sigma interval")

    # fit with Orthogonal distance regression (ODR, fit_type = 0)
    _odr.set_job(fit_type=0)
    res = _odr.run()
    res.pprint()
    xn = np.linspace(0, np.max(x) * 1.25, 500)
    yn = model(res.beta, xn)
    fit_up, fit_dw = confidence_interval(xn, model, res)
    print()
    #stats.loc[2*i+1]=np.concatenate([[Rrs.name, 'ODR'], res.beta, res.sd_beta])

    ax.plot(xn, yn, 'r-', label='ODR',#; '+fit_eq.format(*res.beta),
            linewidth=2)
    ax.fill_between(xn, fit_up, fit_dw, alpha=.25, facecolor="r")  # ,label="1-sigma interval")
    ax.legend()
    ax.set_ylim([0,40])
plt.suptitle('Fit based on modified Gordon model')
stats.to_csv(opj(idir,'stats_SPM_Rrs.csv'), float_format='%.2f',index=False)
plt.tight_layout(rect=[0.0, 0.0, 0.99, 0.94])
fig.savefig(opj(figdir,'blackwater_'+param+'_vs_Rrs.png'), dpi=300)
plt.show()

plt.close()
stats= stats.set_index(['band','algo'])
stats = stats.astype('float')

#------------------------------------------------
# Error histograms
#------------------------------------------------

fig, axs = plt.subplots(nrows=2, ncols=3, figsize=(20, 12))
axs=axs.ravel()
i=0
bandwidth=0.4
for g,d in data.loc[:,['Rrs_sim(B3)','Rrs_sim(B4)','Rrs_sim(B5)','Rrs_sim(B6)','Rrs_sim(B7)','Rrs_sim(B8a)']].iteritems():
    print(g)
    ax=axs[i]
    for method in ('OLS','ODR'):
        c='red'
        if method == 'OLS':
            c='blue'
        stat=stats.loc[(g,method)]

        err = model([stat.b0, stat.b1,stat.b2], d) - y.values.ravel()

        kde = scistats.gaussian_kde(err,bw_method=bandwidth)# / x.std(ddof=1))
        xx = np.linspace(-12,12, 1000)

        ax.hist(err, bins=np.arange(-15,15,bandwidth), color=c,label=method, alpha=0.5, density=True,rwidth=0.8)
        ax.plot(xx, kde.evaluate(xx),color=c,label=method,lw=2)
    ax.axes.get_yaxis().set_visible(False)
    ax.set_title(g)
    ax.set_xlim(-12,12)
    ax.set_xlabel(r'Retrieved - Measured SSSC (mg/L)')
    ax.set_ylabel(r'N')
    ax.legend()
    print(y,)
    i+=1
plt.suptitle('Normalized histogram of errors')
plt.tight_layout(rect=[0.05, 0.05, 0.99, 0.95])
#fig.savefig(opj(figdir,'blackwater_histogram_SPM-err-retrieval.png'), dpi=200)
plt.show()

#------------------------------------------------
# compare with uncertainty
#------------------------------------------------
def sig_param(Rrs,sig_a,sig_b,cov_ab):
    return ((Rrs*sig_a)**2 + sig_b**2 + 2*Rrs*cov_ab)**0.5

fig, axs = plt.subplots(nrows=2, ncols=3, figsize=(20, 12))
axs=axs.ravel()
i=0

for g,d in data.loc[:,['Rrs_sim(B3)','Rrs_sim(B4)','Rrs_sim(B5)','Rrs_sim(B6)','Rrs_sim(B7)','Rrs_sim(B8a)']].iteritems():
    print(g)
    ax=axs[i]
    for method in ('OLS','ODR'):
        c='red'
        if method == 'OLS':
            c='blue'
        stat=stats.loc[(g,method)]

        est = model([stat.b0, stat.b1,stat.b2], d)
        slope, intercept, r_value, p_value, std_err = scistats.linregress(y, est)
        ax.plot([-10, 100], intercept + slope * np.array([-10, 100]), color=c, lw=1.6)
        #sig_est = sig_param(d,stat.sig_a,stat.sig_b,stat.cov_ab)
        #yerr = sig_est,
        ax.errorbar(y, est, xerr=y_sd,  fmt='o', color=c, ecolor=c, label=method+'; y={:.1f}x+{:.1f}; $R^2$={:.3f}'.format(slope, intercept,r_value**2),alpha=0.8)


    ax.set_title(g)

    ax.set_xlabel(r'Measured SSSC (mg/L)')
    ax.set_ylabel(r'Retrieved SSSC (mg/L)')
    ax.plot([-10, 100], [-10, 100], '--', color="grey", lw=1.6)
    ax.set_xlim(0,30)
    ax.set_ylim(0,30)

    ax.legend()
    print(y,)
    i+=1
plt.suptitle('Comparison retrievals')
plt.tight_layout(rect=[0.05, 0.05, 0.99, 0.95])
fig.savefig(opj(figdir,'blackwater_compar_SPM-retrieval.png'), dpi=200)
plt.show()


fa841e375ea6a23c1c6b8fbc54456e78d5a0147f
<mask token>


def inv_gordon88(rrs, g0=0.089, g1=0.125):
    deltas = g0 ** 2 + 4 * g1 * rrs
    return (-g0 + np.sqrt(deltas)) / (2 * g1)


def black_water_model_old(B, x, wl=550):
    aw, bbw = iopw(wl)
    rrs = Rrs2rrs(x)
    u = inv_gordon88(rrs)
    return B[0] ** 2 * (u - np.abs(B[1]))


def black_water_model(B, x, wl=550):
    rrs = Rrs2rrs(x)
    u = inv_gordon88(rrs)
    aw, bbw = iopw(wl)
    N = (u * (aw + bbw + B[0]) - bbw) / (B[1] - u * (B[1] + B[2]))
    return B[3] * N


<mask token>


def exponential(B, x):
    return B[0] * np.exp(B[1] * x) + B[2]


def hyperbolic(B, x, b=1):
    return B[1] / (1 + B[0] * x) ** (1 / b)


def poly(B, x):
    return B[0] + B[1] * x + B[2] * x ** 2


<mask token>


def sig_param(Rrs, sig_a, sig_b, cov_ab):
    return ((Rrs * sig_a) ** 2 + sig_b ** 2 + 2 * Rrs * cov_ab) ** 0.5


<mask token>

<mask token>


def Rrs2rrs(Rrs):
    return Rrs / (0.52 + 1.7 * Rrs)


def inv_gordon88(rrs, g0=0.089, g1=0.125):
    deltas = g0 ** 2 + 4 * g1 * rrs
    return (-g0 + np.sqrt(deltas)) / (2 * g1)


def black_water_model_old(B, x, wl=550):
    aw, bbw = iopw(wl)
    rrs = Rrs2rrs(x)
    u = inv_gordon88(rrs)
    return B[0] ** 2 * (u - np.abs(B[1]))


def black_water_model(B, x, wl=550):
    rrs = Rrs2rrs(x)
    u = inv_gordon88(rrs)
    aw, bbw = iopw(wl)
    N = (u * (aw + bbw + B[0]) - bbw) / (B[1] - u * (B[1] + B[2]))
    return B[3] * N


<mask token>


def exponential(B, x):
    return B[0] * np.exp(B[1] * x) + B[2]


def hyperbolic(B, x, b=1):
    return B[1] / (1 + B[0] * x) ** (1 / b)


def poly(B, x):
    return B[0] + B[1] * x + B[2] * x ** 2


def confidence_interval(xn, model, res, nstd=1):
    """

    :param xn: x-axis data for computation
    :param model: numerical model used for the fit
    :param res: output from scipy.odr.run
    :param nstd: number of sigma to compute confidence interval
    :return: data up and data down
    """
    """
    
    :param res: output from scipy.odr.run
    :param nstd: number of sigma to compute confidence interval
    :return: 
    """
    popt_up = res.beta + nstd * res.sd_beta
    popt_dw = res.beta - nstd * res.sd_beta
    return model(popt_up, xn), model(popt_dw, xn)


<mask token>


def sig_param(Rrs, sig_a, sig_b, cov_ab):
    return ((Rrs * sig_a) ** 2 + sig_b ** 2 + 2 * Rrs * cov_ab) ** 0.5


<mask token>

<mask token>


def Rrs2rrs(Rrs):
    return Rrs / (0.52 + 1.7 * Rrs)


def inv_gordon88(rrs, g0=0.089, g1=0.125):
    deltas = g0 ** 2 + 4 * g1 * rrs
    return (-g0 + np.sqrt(deltas)) / (2 * g1)


def black_water_model_old(B, x, wl=550):
    aw, bbw = iopw(wl)
    rrs = Rrs2rrs(x)
    u = inv_gordon88(rrs)
    return B[0] ** 2 * (u - np.abs(B[1]))


def black_water_model(B, x, wl=550):
    rrs = Rrs2rrs(x)
    u = inv_gordon88(rrs)
    aw, bbw = iopw(wl)
    N = (u * (aw + bbw + B[0]) - bbw) / (B[1] - u * (B[1] + B[2]))
    return B[3] * N


def linear(B, x):
    """Linear function y = m*x + b"""
    return B[0] + B[1] * x


def exponential(B, x):
    return B[0] * np.exp(B[1] * x) + B[2]


def hyperbolic(B, x, b=1):
    return B[1] / (1 + B[0] * x) ** (1 / b)


def poly(B, x):
    return B[0] + B[1] * x + B[2] * x ** 2


def confidence_interval(xn, model, res, nstd=1):
    """

    :param xn: x-axis data for computation
    :param model: numerical model used for the fit
    :param res: output from scipy.odr.run
    :param nstd: number of sigma to compute confidence interval
    :return: data up and data down
    """
    """
    
    :param res: output from scipy.odr.run
    :param nstd: number of sigma to compute confidence interval
    :return: 
    """
    popt_up = res.beta + nstd * res.sd_beta
    popt_dw = res.beta - nstd * res.sd_beta
    return model(popt_up, xn), model(popt_dw, xn)


<mask token>


def sig_param(Rrs, sig_a, sig_b, cov_ab):
    return ((Rrs * sig_a) ** 2 + sig_b ** 2 + 2 * Rrs * cov_ab) ** 0.5


<mask token>

<mask token>
plt.ioff()
mpl.rcParams.update({'font.size': 18})
<mask token>
data.sort_values(by='TSS (mg L)', inplace=True)
<mask token>
if param == 'DOC':
    data.dropna(inplace=True)
<mask token>


def Rrs2rrs(Rrs):
    return Rrs / (0.52 + 1.7 * Rrs)


def inv_gordon88(rrs, g0=0.089, g1=0.125):
    deltas = g0 ** 2 + 4 * g1 * rrs
    return (-g0 + np.sqrt(deltas)) / (2 * g1)


def black_water_model_old(B, x, wl=550):
    aw, bbw = iopw(wl)
    rrs = Rrs2rrs(x)
    u = inv_gordon88(rrs)
    return B[0] ** 2 * (u - np.abs(B[1]))


def black_water_model(B, x, wl=550):
    rrs = Rrs2rrs(x)
    u = inv_gordon88(rrs)
    aw, bbw = iopw(wl)
    N = (u * (aw + bbw + B[0]) - bbw) / (B[1] - u * (B[1] + B[2]))
    return B[3] * N


def linear(B, x):
    """Linear function y = m*x + b"""
    return B[0] + B[1] * x


def exponential(B, x):
    return B[0] * np.exp(B[1] * x) + B[2]


def hyperbolic(B, x, b=1):
    return B[1] / (1 + B[0] * x) ** (1 / b)


def poly(B, x):
    return B[0] + B[1] * x + B[2] * x ** 2


def confidence_interval(xn, model, res, nstd=1):
    """

    :param xn: x-axis data for computation
    :param model: numerical model used for the fit
    :param res: output from scipy.odr.run
    :param nstd: number of sigma to compute confidence interval
    :return: data up and data down
    """
    """
    
    :param res: output from scipy.odr.run
    :param nstd: number of sigma to compute confidence interval
    :return: 
    """
    popt_up = res.beta + nstd * res.sd_beta
    popt_dw = res.beta - nstd * res.sd_beta
    return model(popt_up, xn), model(popt_dw, xn)


if param == 'SPM':
    y = spm
    y_sd = spm_sd
    fit_eq = '$g_0={:.1f};\\ g_1={:.1f};\\ C={:.1f}$'
    ylabel = 'SPM (mg/L)'
    g0 = 0.089
    g1 = 0.125
    model, beta0 = black_water_model_old, [150, 10]
else:
    y = doc
    y_sd = doc_sd
    fit_eq = 'y = {:.1f}/(1+{:.1f}$x)$\n res_var = {:.1f}'
    ylabel = 'DOC (mg/L)'
    model, beta0 = hyperbolic, [10, 10]
<mask token>
for i, ax in enumerate(axs.ravel()):
    wl = wls[i]
    Rrs = data.iloc[:, 2 * i + 2]
    print(Rrs.name, wl)
    Rrs_sd = data.iloc[:, 2 * i + 3]
    x, x_sd = Rrs, Rrs_sd
    ax.errorbar(x, y, xerr=x_sd, yerr=y_sd, fmt='o', color='grey', ecolor=
        'black', alpha=0.8)
    ax.set_title(x.name)
    ax.set_xlabel('$R_{rs}\\ (sr^{-1})$')
    ax.set_ylabel(ylabel)
    testdata = odr.RealData(x, y, sx=x_sd, sy=y_sd)
    _odr = odr.ODR(testdata, odr.Model(model, extra_args=[wl]), beta0=beta0)
    _odr.set_job(fit_type=2)
    res = _odr.run()
    res.pprint()
    xn = np.linspace(0, np.max(x) * 1.25, 500)
    yn = model(res.beta, xn)
    fit_up, fit_dw = confidence_interval(xn, model, res)
    ax.plot(xn, yn, 'b--', label='OLS', linewidth=2)
    ax.fill_between(xn, fit_up, fit_dw, alpha=0.25, facecolor='b')
    _odr.set_job(fit_type=0)
    res = _odr.run()
    res.pprint()
    xn = np.linspace(0, np.max(x) * 1.25, 500)
    yn = model(res.beta, xn)
    fit_up, fit_dw = confidence_interval(xn, model, res)
    print()
    ax.plot(xn, yn, 'r-', label='ODR', linewidth=2)
    ax.fill_between(xn, fit_up, fit_dw, alpha=0.25, facecolor='r')
    ax.legend()
    ax.set_ylim([0, 40])
plt.suptitle('Fit based on modified Gordon model')
stats.to_csv(opj(idir, 'stats_SPM_Rrs.csv'), float_format='%.2f', index=False)
plt.tight_layout(rect=[0.0, 0.0, 0.99, 0.94])
fig.savefig(opj(figdir, 'blackwater_' + param + '_vs_Rrs.png'), dpi=300)
plt.show()
plt.close()
<mask token>
for g, d in data.loc[:, ['Rrs_sim(B3)', 'Rrs_sim(B4)', 'Rrs_sim(B5)',
    'Rrs_sim(B6)', 'Rrs_sim(B7)', 'Rrs_sim(B8a)']].iteritems():
    print(g)
    ax = axs[i]
    for method in ('OLS', 'ODR'):
        c = 'red'
        if method == 'OLS':
            c = 'blue'
        stat = stats.loc[g, method]
        err = model([stat.b0, stat.b1, stat.b2], d) - y.values.ravel()
        kde = scistats.gaussian_kde(err, bw_method=bandwidth)
        xx = np.linspace(-12, 12, 1000)
        ax.hist(err, bins=np.arange(-15, 15, bandwidth), color=c, label=
            method, alpha=0.5, density=True, rwidth=0.8)
        ax.plot(xx, kde.evaluate(xx), color=c, label=method, lw=2)
    ax.axes.get_yaxis().set_visible(False)
    ax.set_title(g)
    ax.set_xlim(-12, 12)
    ax.set_xlabel('Retrieved - Measured SSSC (mg/L)')
    ax.set_ylabel('N')
    ax.legend()
    print(y)
    i += 1
plt.suptitle('Normalized histogram of errors')
plt.tight_layout(rect=[0.05, 0.05, 0.99, 0.95])
plt.show()


def sig_param(Rrs, sig_a, sig_b, cov_ab):
    return ((Rrs * sig_a) ** 2 + sig_b ** 2 + 2 * Rrs * cov_ab) ** 0.5


<mask token>
for g, d in data.loc[:, ['Rrs_sim(B3)', 'Rrs_sim(B4)', 'Rrs_sim(B5)',
    'Rrs_sim(B6)', 'Rrs_sim(B7)', 'Rrs_sim(B8a)']].iteritems():
    print(g)
    ax = axs[i]
    for method in ('OLS', 'ODR'):
        c = 'red'
        if method == 'OLS':
            c = 'blue'
        stat = stats.loc[g, method]
        est = model([stat.b0, stat.b1, stat.b2], d)
        slope, intercept, r_value, p_value, std_err = scistats.linregress(y,
            est)
        ax.plot([-10, 100], intercept + slope * np.array([-10, 100]), color
            =c, lw=1.6)
        ax.errorbar(y, est, xerr=y_sd, fmt='o', color=c, ecolor=c, label=
            method + '; y={:.1f}x+{:.1f}; $R^2$={:.3f}'.format(slope,
            intercept, r_value ** 2), alpha=0.8)
    ax.set_title(g)
    ax.set_xlabel('Measured SSSC (mg/L)')
    ax.set_ylabel('Retrieved SSSC (mg/L)')
    ax.plot([-10, 100], [-10, 100], '--', color='grey', lw=1.6)
    ax.set_xlim(0, 30)
    ax.set_ylim(0, 30)
    ax.legend()
    print(y)
    i += 1
plt.suptitle('Comparison retrievals')
plt.tight_layout(rect=[0.05, 0.05, 0.99, 0.95])
fig.savefig(opj(figdir, 'blackwater_compar_SPM-retrieval.png'), dpi=200)
plt.show()

import os
import numpy as np
import pandas as pd

import matplotlib as mpl
import matplotlib.pyplot as plt

plt.ioff()
mpl.rcParams.update({'font.size': 18})
import statsmodels.api as sm
from sklearn.linear_model import LinearRegression
from statsmodels.sandbox.regression.predstd import wls_prediction_std
from scipy import odr as odr
from scipy import stats as scistats

import RTxploitation as rt
iopw = rt.auxdata.iopw().get_iopw

opj= os.path.join
idir = os.path.abspath('/DATA/OBS2CO/data/rogerio')
figdir = opj(idir,'fig')
file = opj(idir,'data/Simulation_Rrs_OSOAA_TSS_COD_aCDOM.xlsx')

data = pd.read_excel(file)
data.sort_values(by="TSS (mg L)", inplace=True)
data = data.set_index(['Station', 'Date'])

acdom = data.iloc[:, 17]
acdom_sd = data.iloc[:, 18]
#remove sparse data of CDOM
data = data.iloc[:,:-2]
params=['SPM','DOC']
param=params[0]
if param == 'DOC':
    data.dropna(inplace=True)
month = data.index.get_level_values(1).month
spm = data.iloc[:, 0]  # .values.reshape(-1, 1)
spm_sd = data.iloc[:, 1]
doc = data.iloc[:, 15]
doc_sd = data.iloc[:, 16]

def Rrs2rrs(Rrs):
        return Rrs / (0.52 + 1.7 * Rrs)

def inv_gordon88(rrs,g0=0.089,g1=0.125):
        deltas = g0**2 + 4 * g1 * rrs
        return (-g0 + np.sqrt(deltas)) / (2*g1)

def black_water_model_old(B, x , wl=550):
    aw,bbw = iopw(wl)
    rrs = Rrs2rrs(x)
    u = inv_gordon88(rrs)
    return (B[0])**2 * (u - np.abs(B[1])  )

def black_water_model(B, x, wl=550):

    rrs = Rrs2rrs(x)
    u = inv_gordon88(rrs)
    aw,bbw = iopw(wl)
    N = (u*(aw+bbw+B[0])-bbw)/ (B[1]- u *(B[1]+B[2]))

    return B[3] * N

def linear(B, x):
    '''Linear function y = m*x + b'''
    return B[0] + B[1] * x

def exponential(B, x):
    #return B[0] - (B[0]-B[1]) * np.exp(-B[2] * x)
    return B[0] * np.exp(B[1] * x) + B[2]

def hyperbolic(B,x,b=1):
    return B[1]/(1+B[0]*x)**(1/b)

def poly(B, x):
    return B[0] + B[1] * x + B[2] * x ** 2  # + B[3]*x**3

def confidence_interval(xn, model, res, nstd=1):
    '''

    :param xn: x-axis data for computation
    :param model: numerical model used for the fit
    :param res: output from scipy.odr.run
    :param nstd: number of sigma to compute confidence interval
    :return: data up and data down
    '''
    '''
    
    :param res: output from scipy.odr.run
    :param nstd: number of sigma to compute confidence interval
    :return: 
    '''

    popt_up = res.beta + nstd * res.sd_beta
    popt_dw = res.beta - nstd * res.sd_beta

    return model(popt_up, xn), model(popt_dw, xn)

if param == 'SPM':
    y=spm
    y_sd=spm_sd
    fit_eq='$g_0={:.1f};\ g_1={:.1f};\ C={:.1f}$'
    ylabel='SPM (mg/L)'
    # gordon values
    g0 = 0.089
    g1 = 0.125
    model,beta0=black_water_model_old,[150,10]
    #model,beta0=black_water_model,[1,1,1,10]
else:
    y=doc
    y_sd=doc_sd
    fit_eq='y = {:.1f}/(1+{:.1f}$x)$\n res_var = {:.1f}'
    ylabel='DOC (mg/L)'
    model,beta0=hyperbolic,[10,10]

c=doc
cmap = plt.cm.get_cmap("Spectral").reversed()
stats = pd.DataFrame(columns=['band','algo','b0','b1','b2','sig_b0','sig_b1','sig_b2'])


# plot Rrs = f(SPM)
wls = [560, 665, 704,740, 782, 864]
fig, axs = plt.subplots(nrows=2, ncols=3, figsize=(20, 12))
for i, ax in enumerate(axs.ravel()):
    wl=wls[i]

    Rrs = data.iloc[:, 2 * i + 2]
    print(Rrs.name,wl)
    Rrs_sd = data.iloc[:, 2 * i + 3]

    x,x_sd = Rrs, Rrs_sd

    ax.errorbar(x, y, xerr=x_sd, yerr=y_sd, fmt='o', color='grey', ecolor='black', alpha=0.8)
    # to put color in scatter dots
    #ax.scatter(x, y, c=c, s=55, cmap=cmap)
    ax.set_title(x.name)
    ax.set_xlabel(r'$R_{rs}\ (sr^{-1})$')
    ax.set_ylabel(ylabel)

    testdata = odr.RealData(x, y, sx=x_sd, sy=y_sd)
    _odr = odr.ODR(testdata, odr.Model(model, extra_args=[wl]), beta0=beta0) #

    # fit with ordinary least square (OLS, fit_type=2 )
    _odr.set_job(fit_type=2)
    res = _odr.run()

    res.pprint()
    xn = np.linspace(0, np.max(x) * 1.25, 500)
    yn = model(res.beta, xn)
    fit_up, fit_dw = confidence_interval(xn, model, res)

    #stats.loc[2*i]=np.concatenate([[Rrs.name, 'OLS'], res.beta, res.sd_beta])


    ax.plot(xn, yn, 'b--', label='OLS',#; '+fit_eq.format(*res.beta),
            linewidth=2)
    ax.fill_between(xn, fit_up, fit_dw, alpha=.25, facecolor="b")  # ,label="1-sigma interval")

    # fit with Orthogonal distance regression (ODR, fit_type = 0)
    _odr.set_job(fit_type=0)
    res = _odr.run()
    res.pprint()
    xn = np.linspace(0, np.max(x) * 1.25, 500)
    yn = model(res.beta, xn)
    fit_up, fit_dw = confidence_interval(xn, model, res)
    print()
    #stats.loc[2*i+1]=np.concatenate([[Rrs.name, 'ODR'], res.beta, res.sd_beta])

    ax.plot(xn, yn, 'r-', label='ODR',#; '+fit_eq.format(*res.beta),
            linewidth=2)
    ax.fill_between(xn, fit_up, fit_dw, alpha=.25, facecolor="r")  # ,label="1-sigma interval")
    ax.legend()
    ax.set_ylim([0,40])
plt.suptitle('Fit based on modified Gordon model')
stats.to_csv(opj(idir,'stats_SPM_Rrs.csv'), float_format='%.2f',index=False)
plt.tight_layout(rect=[0.0, 0.0, 0.99, 0.94])
fig.savefig(opj(figdir,'blackwater_'+param+'_vs_Rrs.png'), dpi=300)
plt.show()

plt.close()
stats= stats.set_index(['band','algo'])
stats = stats.astype('float')

#------------------------------------------------
# Error histograms
#------------------------------------------------

fig, axs = plt.subplots(nrows=2, ncols=3, figsize=(20, 12))
axs=axs.ravel()
i=0
bandwidth=0.4
for g,d in data.loc[:,['Rrs_sim(B3)','Rrs_sim(B4)','Rrs_sim(B5)','Rrs_sim(B6)','Rrs_sim(B7)','Rrs_sim(B8a)']].iteritems():
    print(g)
    ax=axs[i]
    for method in ('OLS','ODR'):
        c='red'
        if method == 'OLS':
            c='blue'
        stat=stats.loc[(g,method)]

        err = model([stat.b0, stat.b1,stat.b2], d) - y.values.ravel()

        kde = scistats.gaussian_kde(err,bw_method=bandwidth)# / x.std(ddof=1))
        xx = np.linspace(-12,12, 1000)

        ax.hist(err, bins=np.arange(-15,15,bandwidth), color=c,label=method, alpha=0.5, density=True,rwidth=0.8)
        ax.plot(xx, kde.evaluate(xx),color=c,label=method,lw=2)
    ax.axes.get_yaxis().set_visible(False)
    ax.set_title(g)
    ax.set_xlim(-12,12)
    ax.set_xlabel(r'Retrieved - Measured SSSC (mg/L)')
    ax.set_ylabel(r'N')
    ax.legend()
    print(y,)
    i+=1
plt.suptitle('Normalized histogram of errors')
plt.tight_layout(rect=[0.05, 0.05, 0.99, 0.95])
#fig.savefig(opj(figdir,'blackwater_histogram_SPM-err-retrieval.png'), dpi=200)
plt.show()

#------------------------------------------------
# compare with uncertainty
#------------------------------------------------
def sig_param(Rrs,sig_a,sig_b,cov_ab):
    return ((Rrs*sig_a)**2 + sig_b**2 + 2*Rrs*cov_ab)**0.5

fig, axs = plt.subplots(nrows=2, ncols=3, figsize=(20, 12))
axs=axs.ravel()
i=0

for g,d in data.loc[:,['Rrs_sim(B3)','Rrs_sim(B4)','Rrs_sim(B5)','Rrs_sim(B6)','Rrs_sim(B7)','Rrs_sim(B8a)']].iteritems():
    print(g)
    ax=axs[i]
    for method in ('OLS','ODR'):
        c='red'
        if method == 'OLS':
            c='blue'
        stat=stats.loc[(g,method)]

        est = model([stat.b0, stat.b1,stat.b2], d)
        slope, intercept, r_value, p_value, std_err = scistats.linregress(y, est)
        ax.plot([-10, 100], intercept + slope * np.array([-10, 100]), color=c, lw=1.6)
        #sig_est = sig_param(d,stat.sig_a,stat.sig_b,stat.cov_ab)
        #yerr = sig_est,
        ax.errorbar(y, est, xerr=y_sd,  fmt='o', color=c, ecolor=c, label=method+'; y={:.1f}x+{:.1f}; $R^2$={:.3f}'.format(slope, intercept,r_value**2),alpha=0.8)


    ax.set_title(g)

    ax.set_xlabel(r'Measured SSSC (mg/L)')
    ax.set_ylabel(r'Retrieved SSSC (mg/L)')
    ax.plot([-10, 100], [-10, 100], '--', color="grey", lw=1.6)
    ax.set_xlim(0,30)
    ax.set_ylim(0,30)

    ax.legend()
    print(y,)
    i+=1
plt.suptitle('Comparison retrievals')
plt.tight_layout(rect=[0.05, 0.05, 0.99, 0.95])
fig.savefig(opj(figdir,'blackwater_compar_SPM-retrieval.png'), dpi=200)
plt.show()

