/***************************************************/
/* gccBTS.c */
/***************************************************/
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <time.h>
#define NDmax 1000
#define NED 10000
#define NYR 14
int COMP_IS(const void *a,const void *b);
int COMP_DS(const void *a,const void *b);
void IRND(int nd,int idata[NDmax]);
void main_part(int nmethod,char fnameR[50],char fnameW[50],char fnameH[50]);
void GUM_COEF(int nd,double datat[NDmax],double *a_gum,double *c_gum);
void GEV_COEF(int nd,double datat[NDmax],double *k_gev,double *a_gev,double *c_gev);
void SQR_COEF(int nd,double datat[NDmax],double *a_sqr,double *b_sqr);
void LN3_COEF(int nd,double datat[NDmax],double *a_ln3,double *my_ln3,double *sy_ln3);
void LP3_COEF(int nd,double datat[NDmax],double *a_lp3,double *b_lp3,double *c_lp3);
double XP_GUM(double p,double a_gum,double c_gum);
double XP_GEV(double p,double k_gev,double a_gev,double c_gev);
double XP_SQR(double p,double a_sqr,double b_sqr,double xmax);
double XP_LN3(double p,double a_ln3,double my_ln3,double sy_ln3);
double XP_LP3(double p,double a_lp3,double b_lp3,double c_lp3);
double XP_LP3_SP(int nd,double datat[NDmax],double p);
double FSQR(int nd,double datat[NDmax],double bb,double *a1,double *a2);
double GAMMAL(double x);
double TODA(double rx);
double GDISTPP(double eta,double p);
double LGAMMA1(double xx);
double LGAMMA2(double x);
double LINCGF(double y1,double z1);
/*---------------------------------------------------------------------------*/
int main(int argc,char *argv[])
{
    char fnameR[50]; /*Input file name*/
    char fnameW[50]; /*Output file name*/
    char fnameH[50]; /*Output file name for histogram*/
    int nmethod;

    nmethod=atoi(argv[1]);
    strcpy(fnameR,argv[2]);
    strcpy(fnameW,argv[3]);
    strcpy(fnameH,argv[4]);

    main_part(nmethod,fnameR,fnameW,fnameH);

    return 0;
}
/*---------------------------------------------------------------------------*/
void main_part(int nmethod,char fnameR[50],char fnameW[50],char fnameH[50])
{
    int i,j,iii,jjj,nd,kl,ku;
    int idata[NDmax+1];
    double datax[NDmax+1],datay[NDmax+1],edata[NED+1];
    double p,yy,eave;
    double a_gum,c_gum;
    double a_gev,c_gev,k_gev;
    double a_ln3,my_ln3,sy_ln3;
    double a_lp3,b_lp3,c_lp3;
    double a_sqr,b_sqr;
    int nnn=NYR;
    double Ty[]={0,2,5,10,20,30,50,100,200,300,500,1000,2000,3000,5000};
    double res_org[NYR+1],wk[NYR+1];
    char dat[500]="";
    char sv[50]="";
    char fnameT[]="_temp.csv";

    FILE *fin,*fout,*fhis;

    srand((unsigned int)time(NULL));
    for(i=0;i<=NDmax;i++){datax[i]=0.0;datay[i]=0.0;idata[i]=0;}
    for(i=0;i<=NED;i++){edata[i]=0.0;}
    /*============*/
    /* Data input */
    /*============*/
    fin=fopen(fnameR,"r");
    fgets(dat,sizeof dat,fin);
    i=0;
    while(fgets(dat,sizeof dat,fin)!=NULL){i=i+1;datax[i]=atof(dat);}
    fclose(fin);
    nd=i;
    for(i=0;i<=nd;i++){idata[i]=i;}
    /*========================*/
    /* Sorting in small order */
    /*========================*/
    datax[0]=-99.9;
    qsort(datax,nd+1,sizeof(datax[0]),COMP_DS);
    /*========================*/
    /* Calc. of Org.data      */
    /*========================*/
    switch(nmethod){
        case 1:GUM_COEF(nd,datax,&a_gum,&c_gum);         break;
        case 2:GEV_COEF(nd,datax,&k_gev,&a_gev,&c_gev);  break;
        case 3:SQR_COEF(nd,datax,&a_sqr,&b_sqr);         break;
        case 4:LN3_COEF(nd,datax,&a_ln3,&my_ln3,&sy_ln3);break;
        case 5:LP3_COEF(nd,datax,&a_lp3,&b_lp3,&c_lp3);  break;
        }
    for(iii=1;iii<=nnn;iii++){
        p=1.0-1.0/Ty[iii];
        switch(nmethod){
            case 1:yy=XP_GUM(p,a_gum,c_gum);          break;
            case 2:yy=XP_GEV(p,k_gev,a_gev,c_gev);    break;
            case 3:yy=XP_SQR(p,a_sqr,b_sqr,datax[nd]);break;
            case 4:yy=XP_LN3(p,a_ln3,my_ln3,sy_ln3);  break;
            case 5:
                if(b_lp3<1e4){
                    yy=XP_LP3(p,a_lp3,b_lp3,c_lp3);
                }else{
                    yy=XP_LP3_SP(nd,datax,p);
                }
                break;
        }
        res_org[iii]=yy;
    }
    /*========================*/
    /* Calc. of BTS           */
    /*========================*/
    fout=fopen(fnameT,"w");
    for(i=1;i<=NED;i++){
        IRND(nd,idata);
        for(j=1;j<=nd;j++){datay[j]=datax[idata[j]];}
        switch(nmethod){
            case 1:GUM_COEF(nd,datay,&a_gum,&c_gum);          break;
            case 2:GEV_COEF(nd,datay,&k_gev,&a_gev,&c_gev);   break;
            case 3:SQR_COEF(nd,datay,&a_sqr,&b_sqr);          break;
            case 4:LN3_COEF(nd,datay,&a_ln3,&my_ln3,&sy_ln3); break;
            case 5:LP3_COEF(nd,datay,&a_lp3,&b_lp3,&c_lp3);   break;
        }
        sprintf(sv,"%d",i);strcpy(dat,sv);
        for(iii=1;iii<=nnn;iii++){
            p=1.0-1.0/Ty[iii];
            switch(nmethod){
                case 1:yy=XP_GUM(p,a_gum,c_gum);           break;
                case 2:yy=XP_GEV(p,k_gev,a_gev,c_gev);     break;
                case 3:yy=XP_SQR(p,a_sqr,b_sqr,datax[nd]); break;
                case 4:yy=XP_LN3(p,a_ln3,my_ln3,sy_ln3);   break;
                case 5:
                    if(b_lp3<1e4){
                        yy=XP_LP3(p,a_lp3,b_lp3,c_lp3);
                    }else{
                        yy=XP_LP3_SP(nd,datay,p);
                    }
                    break;
            }
             sprintf(sv,"%.3f",yy);strcat(dat,",");strcat(dat,sv);
        }
        strcat(dat,"\n");fputs(dat,fout);
    }
    fclose(fout);

    fout=fopen(fnameW,"w");
    switch(nmethod){
        case 1:fprintf(fout,"Gumbel\n");break;
        case 2:fprintf(fout,"GEV\n");break;
        case 3:fprintf(fout,"SQRT-ET\n");break;
        case 4:fprintf(fout,"LN3\n");break;
        case 5:fprintf(fout,"LP3\n");break;
    }
    fprintf(fout,"year,prob,org,bts,lcl,ucl\n");
    for(iii=1;iii<=nnn;iii++){
        edata[0]=-99.9;
        fin=fopen(fnameT,"r");
        for(i=1;i<=NED;i++){
            fgets(dat,sizeof dat,fin);
            yy=atof(strtok(dat,","));
            for(jjj=1;jjj<=nnn;jjj++){
                wk[jjj]=atof(strtok(NULL,","));
            }
            edata[i]=wk[iii];
        }
        fclose(fin);
        qsort(edata,NED+1,sizeof(edata[0]),COMP_DS);
        eave=0.0;
        for(i=1;i<=NED;i++){eave=eave+edata[i];}
        eave=eave/(double)NED;
        kl=(int)(0.025*(double)NED);
        ku=NED-kl;
        printf("%4.0f %7.5f %4.1f %4.1f %4.1f %4.1f\n",Ty[iii],p,res_org[iii],eave,edata[kl],edata[ku]);
        fprintf(fout,"%.0f,%.5f,%.3f,%.3f,%.3f,%.3f\n",Ty[iii],p,res_org[iii],eave,edata[kl],edata[ku]);
        /*========================*/
        /* Output for histogram   */
        /*========================*/
        if(iii==nnn){
            fhis=fopen(fnameH,"w");
            switch(nmethod){
                case 1:fprintf(fhis,"Gumbel 5000year bootstrap\n");break;
                case 2:fprintf(fhis,"GEV 5000year bootstrap\n");break;
                case 3:fprintf(fhis,"SQRT-ET 5000year bootstrap\n");break;
                case 4:fprintf(fhis,"LN3 5000year bootstrap\n");break;
                case 5:fprintf(fhis,"LP3 5000year bootstrap\n");break;
            }
            for(i=1;i<=NED;i++){fprintf(fhis,"%g\n",edata[i]);}
            fclose(fhis);
        }
    }
    fclose(fout);
}
/*---------------------------------------------------------------------------*/
int COMP_IS(const void *a,const void *b)
{
    /*小さい順*/
    return *(int*)a-*(int*)b;
}
/*---------------------------------------------------------------------------*/
int COMP_DS(const void *a,const void *b)
{
    /*比較関数は戻り値が整数型なので*/
    /*整数値を返すようにする*/
    /*小さい順*/
    if(*(double*)a> *(double*)b)return  1;
    if(*(double*)a==*(double*)b)return  0;
    if(*(double*)a< *(double*)b)return -1;
}
/*---------------------------------------------------------------------------*/
void IRND(int nd,int idata[NDmax])
{
    int i;
    double r,a,b;

    idata[0]=0;
    a=1.0;
    b=(double)nd;
    for(i=1;i<=nd;i++){
        r=a+(b-a)*(double)rand()/(double)RAND_MAX;
        idata[i]=(int)(r+0.5);
    }
    qsort(idata,nd+1,sizeof(idata[0]),COMP_IS);
}
/*---------------------------------------------------------------------------*/
void GUM_COEF(int nd,double datat[NDmax],double *a_gum,double *c_gum)
{
    /* Gumbel distribution  */
    int i;
    double b0,b1,b2;
    double lam1,lam2,lam3;

    /* datat[] is sorted data */
    b0=0;b1=0;b2=0;
    for(i=1;i<=nd;i++){
        b0=b0+datat[i];
        b1=b1+(double)(i-1)*datat[i];
        b2=b2+(double)((i-1)*(i-2))*datat[i];
    }
    b0=b0/(double)nd; /*Average*/
    b1=b1/(double)(nd*(nd-1));
    b2=b2/(double)(nd*(nd-1)*(nd-2));
    lam1=b0;
    lam2=2.0*b1-b0;
    lam3=6.0*b2-6.0*b1+b0;

    /*Gumbel F(x)=exp{-exp[-(x-c)/a]} T=1/[1-F(x)] */
    *a_gum=lam2/log(2.0);
    *c_gum=lam1-0.5772*(*a_gum);
}
/*---------------------------------------------------------------------------*/
void GEV_COEF(int nd,double datat[NDmax],double *k_gev,double *a_gev,double *c_gev)
{
    /* Generalized Extreme Value distribution  */
    int i;
    double b0,b1,b2,d;
    double lam1,lam2,lam3;

    /* datat[] is sorted data */
    b0=0;b1=0;b2=0;
    for(i=1;i<=nd;i++){
        b0=b0+datat[i];
        b1=b1+(double)(i-1)*datat[i];
        b2=b2+(double)((i-1)*(i-2))*datat[i];
    }
    b0=b0/(double)nd; /*Average*/
    b1=b1/(double)(nd*(nd-1));
    b2=b2/(double)(nd*(nd-1)*(nd-2));
    lam1=b0;
    lam2=2.0*b1-b0;
    lam3=6.0*b2-6.0*b1+b0;

    /*GEV F(x)=exp{-[1-k*(x-c)/a]^(1/k)} T=1/[1-F(x)] */
    d=2.0*lam2/(lam3+3.0*lam2)-log(2.0)/log(3.0);
    *k_gev=7.8590*d+2.9554*d*d;
    *a_gev=(*k_gev)*lam2/(1.0-pow(2.0,-(*k_gev)))/GAMMAL(1.0+(*k_gev));
    *c_gev=lam1-((*a_gev)/(*k_gev))*(1.0-GAMMAL(1.0+(*k_gev)));
}
/*---------------------------------------------------------------------------*/
void SQR_COEF(int nd,double datat[NDmax],double *a_sqr,double *b_sqr)
{
    int i;
    double a1,a2,b1,b2,bb,f1,f2,ff;

    /* Condition for a1>0 */
    bb=0.0;
    for(i=1;i<=nd;i++){bb=bb+sqrt(datat[i]);}
    bb=(2.0*(double)nd/bb)*(2.0*(double)nd/bb)+0.001;

    /* Bisection method */
    b1=bb;
    b2=b1+0.5;
    bb=0.5*(b1+b2);
    f1=FSQR(nd,datat,b1,&a1,&a2);
    f2=FSQR(nd,datat,b2,&a1,&a2);
    ff=FSQR(nd,datat,bb,&a1,&a2);
    do{
        if(f1*ff<0.0)b2=bb;
        if(ff*f2<0.0)b1=bb;
        if(ff==0.0)break;
        if(0.0<f1*ff&&0.0<ff*f2){b1=b2;b2=b1+0.5;}
        bb=0.5*(b1+b2);
        f1=FSQR(nd,datat,b1,&a1,&a2);
        f2=FSQR(nd,datat,b2,&a1,&a2);
        ff=FSQR(nd,datat,bb,&a1,&a2);
    }while(0.001<fabs(a1-a2));

    *a_sqr=0.5*(a1+a2);
    *b_sqr=bb;
}
/*---------------------------------------------------------------------------*/
double FSQR(int nd,double datat[NDmax],double bb,double *a1,double *a2)
{
    int i;
    double bxr,bxe,bx1e,bxre;

    bxr=0.0;bxe=0.0;bx1e=0.0;bxre=0.0;
    for(i=1;i<=nd;i++){
        bxr=bxr+sqrt(bb*datat[i]);
        bxe=bxe+exp(-sqrt(bb*datat[i]));
        bx1e=bx1e+bb*datat[i]*exp(-sqrt(bb*datat[i]));
        bxre=bxre+sqrt(bb*datat[i])*exp(-sqrt(bb*datat[i]));
    }
    *a1=(bxr-2.0*(double)nd)/bx1e;
    *a2=(double)nd/(bxe+bxre);
    return (*a1)-(*a2);
}
/*---------------------------------------------------------------------------*/
void LN3_COEF(int nd,double datat[NDmax],double *a_ln3,double *my_ln3,double *sy_ln3)
{
    /* Log-normal distribution 3-parameter quantile-method */
    /* search of median xm*/
    int i;
    double xm,w;

    xm=0.5*(datat[(int)(nd/2)+1]+datat[nd-(int)(nd/2)]);
    *a_ln3=(datat[1]*datat[nd]-xm*xm)/(datat[1]+datat[nd]-2.0*xm);
    w=0.0;
    for(i=1;i<=nd;i++){w=w+log(datat[i]-(*a_ln3));}
    *my_ln3=w/(double)nd;
    w=0.0;
    for(i=1;i<=nd;i++){w=w+(log(datat[i]-(*a_ln3))-(*my_ln3))*(log(datat[i]-(*a_ln3))-(*my_ln3));}
    *sy_ln3=sqrt(w/(double)nd);
}
/*---------------------------------------------------------------------------*/
void LP3_COEF(int nd,double datat[NDmax],double *a_lp3,double *b_lp3,double *c_lp3)
{
    /* Log Pearson III distribution (Gamma distribution) */
    int i;
    double wmy,wsy,wry,ws,wc,w1,w2;
    
    wmy=0.0;
    for(i=1;i<=nd;i++){wmy=wmy+log(datat[i]);}
    wmy=wmy/(double)nd;
    ws=0.0;
    for(i=1;i<=nd;i++){ws=ws+(log(datat[i])-wmy)*(log(datat[i])-wmy);}
    ws=sqrt(ws/(double)nd);
    wsy=sqrt((double)nd/(double)(nd-1))*ws;
    wc=0.0;
    for(i=1;i<=nd;i++){wc=wc+pow((log(datat[i])-wmy)/ws,3.0);}
    wc=wc/(double)nd;
    w1=1.0+6.51/(double)nd+20.2/(double)nd/(double)nd;
    w2=1.48/(double)nd+6.77/(double)nd/(double)nd;
    wry=wc*(w1+w2*wc*wc);

    *b_lp3=4.0/(wry*wry);
    *a_lp3=wsy/sqrt(*b_lp3);
    if(wry<0.0)*a_lp3=-(*a_lp3);
    *c_lp3=wmy-(*a_lp3)*(*b_lp3);
}
/*---------------------------------------------------------------------------*/
double XP_GUM(double p,double a_gum,double c_gum)
{
    return c_gum-a_gum*log(-log(p));
}
/*---------------------------------------------------------------------------*/
double XP_GEV(double p,double k_gev,double a_gev,double c_gev)
{
    return c_gev+(a_gev/k_gev)*(1.0-pow(-log(p),k_gev));
}
/*---------------------------------------------------------------------------*/
double XP_SQR(double p,double a_sqr,double b_sqr,double xmax)
{
    double x1,x2,ff,df,tp;

    tp=sqrt(b_sqr*xmax); /* Initial value of tp */
        ff=log(1.0+tp)-tp-log(-1.0/a_sqr*log(p));
        df=1.0/(1.0+tp)-1.0;
        tp=tp-ff/df;
        x2=tp*tp/b_sqr; /* Iniial value of x2 */
    do{
        x1=x2;
        ff=log(1.0+tp)-tp-log(-1.0/a_sqr*log(p));
        df=1.0/(1.0+tp)-1.0;
        tp=tp-ff/df;
        x2=tp*tp/b_sqr;
    }while(0.001<fabs(x2-x1));

    return x2;
}
/*---------------------------------------------------------------------------*/
double XP_LN3(double p,double a_ln3,double my_ln3,double sy_ln3)
{
    return a_ln3+exp(my_ln3+sy_ln3*TODA(p));
}
/*-----------------------------------------------------------------------*/
double XP_LP3(double p,double a_lp3,double b_lp3,double c_lp3)
{
    double q;

    q=p;
    if(a_lp3<0.0)q=1.0-p;
    return exp(c_lp3+a_lp3*GDISTPP(b_lp3,q));
}
/*---------------------------------------------------------------------------*/
double XP_LP3_SP(int nd,double datat[NDmax],double p)
{
    /* Log Pearson III distribution (Gamma distribution) */
    int i;
    double wmy,wsy,wry,ws,wc,w1,w2;
    double q,kp,ry;

    wmy=0.0;
    for(i=1;i<=nd;i++){wmy=wmy+log(datat[i]);}
    wmy=wmy/(double)nd;
    ws=0.0;
    for(i=1;i<=nd;i++){ws=ws+(log(datat[i])-wmy)*(log(datat[i])-wmy);}
    ws=sqrt(ws/(double)nd);
    wsy=sqrt((double)nd/(double)(nd-1))*ws;
    wc=0.0;
    for(i=1;i<=nd;i++){wc=wc+pow((log(datat[i])-wmy)/ws,3.0);}
    wc=wc/(double)nd;
    w1=1.0+6.51/(double)nd+20.2/(double)nd/(double)nd;
    w2=1.48/(double)nd+6.77/(double)nd/(double)nd;
    wry=wc*(w1+w2*wc*wc);

    q=p;
    ry=wry;
    kp=2.0/ry*pow((1.0+ry*TODA(q)/6.0-ry*ry/36.0),3.0)-2.0/ry;

    return exp(wmy+wsy*kp);
}
/*---------------------------------------------------------------------------*/
double GAMMAL(double x)
{
    return exp(LGAMMA1(x));
}
/*---------------------------------------------------------------------------*/
double TODA(double rx)
{
    /* %-point of Non-exceeding probability */
    int i;
    double ry,ay,uay;
    double b[11],sum;

    if(rx==0.5){
        uay=0.0;
    }else{
        if(rx<0.5){
            ry=rx;
        }else{
            ry=1-rx;
        }
        ay=-log(4.0*ry*(1.0-ry));
        b[0] = 1.570796288;
        b[1] = 0.03706987906;
        b[2] =-0.0008364353589;
        b[3] =-0.0002250947176;
        b[4] = 0.000006841218299;
        b[5] = 0.000005824238515;
        b[6] =-0.00000104527497;
        b[7] = 0.00000008360937017;
        b[8] =-0.000000003231081277;
        b[9] = 0.00000000003657763036;
        b[10]= 0.0000000000006936233982;
        sum=b[0];
        for(i=1;i<=10;i++){sum=sum+b[i]*pow(ay,(double)i);}
        uay=sqrt(ay*sum);
        if(rx<0.5)uay=-uay;
    }
    return uay;
}
/*-----------------------------------------------------------------------*/
double GDISTPP(double eta,double p)
{
    /* eta=b_lp3,p=probability */
    int k;
    double etad,bd,et,g1,x,dk,xx,a,b,expa,q,qd,er;

    etad=eta;
    bd=p;
    et=log(etad);
    g1=log(bd)+LGAMMA1(etad);
    x=(et+g1)*1.0/etad;
    k=0;
    do{
        k=k+1;
        dk=(double)k;
        xx=x+log(dk);
    }while(LINCGF(xx,et)<g1);
    a=0.0;
    if(2<=k)a=x+log(dk-1.0);
    b=x+log(dk);
    expa=exp(a);
    if(k<=1)expa=0.0;
    k=0;
    do{
        k=k+1;
        q=0.5*(expa+exp(b));
        qd=log(q);
        er=LINCGF(qd,et)-g1;
        if(er<0.0){a=qd;expa=exp(a);}
        if(fabs(er)<1e-10)break;
        if(0.0<er)b=qd;
    }while(k<100000);

    return q;
}
/*-----------------------------------------------------------------------*/
double LGAMMA1(double xx)
{
    int i,m;
    double pi=M_PI;
    double x,gs,y,z,y1;

    x=xx;
    if(30.0<=x){
        gs=(x-0.5)*log(x)-x+0.5*log(2.0*pi)+1.0/(12.0*x)-1.0/(360.0*pow(x,3.0))+1.0/(1260.0*pow(x,5.0))-1.0/(1680.0*pow(x,7.0));
        if(10000.0<x)return gs;
        gs=gs+5.0/(5940.0*pow(x,9.0))-691.0/(360360.0*pow(x,11.0));
        return gs;
    }else{
        m=(int)x;
        y=x-(double)m;
        gs=LGAMMA2(y);
        if(x<=2.0){
            if(m<1)gs=gs-log(y);
            return gs;
        }else{
            z=0;
            for(i=1;i<=m-1;i++){y1=y+(double)i;z=z+log(y1);}
            gs=gs+z;
            return gs;
        }
    }
}
/*-----------------------------------------------------------------------*/
double LGAMMA2(double x)
{
    int i;
    double c[9],g;

    c[1]=-0.577191652;
    c[2]= 0.988205891;
    c[3]=-0.897056937;
    c[4]= 0.918206857;
    c[5]=-0.756704078;
    c[6]= 0.482199394;
    c[7]=-0.193527818;
    c[8]= 0.035868343;
    g=1.0;
    for(i=1;i<=8;i++){g=g+c[i]*pow(x,(double)i);}
    g=log(g);
    return g;
}
/*-----------------------------------------------------------------------*/
double LINCGF(double y1,double z1)
{
    int j;
    double z,y,f,d1,d2;

    z=exp(z1);
    y=exp(y1);
    f=1.0/z;
    d1=1.0/z;
    j=0;
    do{
        j=j+1;
        d2=d1*y/(z+(double)j);
        f=f+d2;
        if(y/(z+(double)j+1.0-y)*d2<f*1e-10&&y/(z+(double)j+1.0)<1.0)break;
        d1=d2;
    }while(j<100000);

    return z*y1-y+log(f);
}
/*-----------------------------------------------------------------------*/