import numpy as np
from math import *
import sys
from scipy.stats import gamma
from scipy.special import gamma as Gamma
import matplotlib.pyplot as plt

def MSR(n,z):
    mz=np.sum(z)/float(n)
    sz=np.sqrt(np.sum((z-mz)**2)/float(n-1))
    _s=np.sqrt(np.sum((z-mz)**2)/float(n))
    _c=np.sum(((z-mz)/_s)**3)/float(n)
    rz=sqrt(float(n)*float(n-1))/float(n-2)*_c
    return mz,sz,rz

def X_LP3(n,x,alpha,p):
    y=np.log(np.array(x))
    ym,ys,yr=MSR(n,y)
    bb=4.0/yr**2
    aa=ys/sqrt(bb)
#    if yr<0.0: aa=-ys/sqrt(bb)
    cc=ym-aa*bb
    xp=exp(cc+aa*gamma.ppf(p,bb))
    return xp

def DRAWHIST(data1,data2,data3,nn,str_t):
    xmin=0; xmax=1000
    ymin=0; ymax=int(nn/5)
    fig=plt.figure()
    ax1=fig.add_subplot(2,1,1)
    ax1.set_xlabel('Maximum daily rainfall (mm/day)', size=14)
    ax1.set_ylabel('Frequency', size=14)
    ax1.set_xlim(xmin,xmax)
    ax1.set_ylim(ymin,ymax)
    ax1.grid(True)
    ndx=int((xmax-xmin)/20)
    ax1.hist(data1,bins=ndx,range=(xmin,xmax),alpha=0.2, histtype='stepfilled', color='red',label='100years')
    ax1.hist(data2,bins=ndx,range=(xmin,xmax),alpha=0.2, histtype='stepfilled', color='yellow',label='1000years')
    ax1.hist(data3,bins=ndx,range=(xmin,xmax),alpha=0.2, histtype='stepfilled', color='blue',label='10000years')
    ax1.legend(fontsize=10,title='LP3',loc='upper right')
    ax2=fig.add_subplot(2,1,2)
    ax2.set_xlabel('Maximum daily rainfall (mm/day)', size=14)
    ax2.set_ylabel('Cumulative probability', size=14)
    ax2.set_ylim(0.0,1.0)
    ax2.grid(True)
    ax2.hist(data1, bins = nn, range=(xmin,xmax), normed=True, cumulative=True, histtype='step',color='red',label='100years')
    ax2.hist(data2, bins = nn, range=(xmin,xmax), normed=True, cumulative=True, histtype='step',color='lime',label='1000years')
    ax2.hist(data3, bins = nn, range=(xmin,xmax), normed=True, cumulative=True, histtype='step',color='blue',label='10000years')
    ax2.axvline(x=np.median(data1),linestyle='--',color='red')
    ax2.axvline(x=np.median(data2),linestyle='--',color='lime')
    ax2.axvline(x=np.median(data3),linestyle='--',color='blue')
    txt1='{0:5.0f}'.format(np.median(data1))
    txt2='{0:5.0f}'.format(np.median(data2))
    txt3='{0:5.0f}'.format(np.median(data3))
    ax2.text(np.median(data1),0.95,txt1,fontname='monospace',fontsize=12,color='red', ha='right',va='top',rotation=90)
    ax2.text(np.median(data2),0.95,txt2,fontname='monospace',fontsize=12,color='lime',ha='right',va='top',rotation=90)
    ax2.text(np.median(data3),0.95,txt3,fontname='monospace',fontsize=12,color='blue',ha='right',va='top',rotation=90)
    ax2.legend(fontsize=10,title='LP3',loc='lower right',handleheight=0.1)
    plt.tight_layout()
    plt.savefig('fig_bs_LP3_M.png',dpi=200)
    plt.show()
    plt.close()


# Main routine
alpha=0.4
nn=10000
param=sys.argv
fnameR=param[1]

# Input data
x_org=[]
fin=open(fnameR,'r')
text=fin.readline()
while 1:
    text=fin.readline()
    if not text: break
    text=text.strip()
    text=text.split()
    x_org=x_org+[float(text[0])]
fin.close()
n=len(x_org)

str_t='LP3'
years=[100,1000,10000]
for k in range(0,3):
    p=1.0-1.0/years[k]
    xp=[]
    for ii in range(0,nn):
        num=np.random.randint(0,n,size=n)
        bs=[]
        for i in range(0,n):
            bs=bs+[x_org[num[i]]]
        bs.sort()
        xp=xp+[X_LP3(n,bs,alpha,p)]
        del bs
    if k==0: data1=xp
    if k==1: data2=xp
    if k==2: data3=xp
    del xp
DRAWHIST(data1,data2,data3,nn,str_t)