/***************************************************/
/* gccSTA.c */
/***************************************************/
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#define NDmax 1000
void main_part(char fnameR[50],char fnameW[50]);
double COC(int nd,double datat[NDmax],double datay[NDmax]);
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 fnameP[50]; /*Plotting file name*/

    strcpy(fnameR,argv[1]);
    strcpy(fnameW,argv[2]);
/*    strcpy(fnameP,argv[3]);*/

    main_part(fnameR,fnameW);

    return 0;
}
/*---------------------------------------------------------------------------*/

void main_part(char fnameR[50],char fnameW[50])
{
    int i,j,nd,nnn;
    double datax[NDmax+1],datay[NDmax+1];
    double alpha;
    double p,yy1,yy2,yy3,yy4,yy5,work;
    double rr1,rr2,rr3,rr4,rr5;
    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;
    char dat[256]="";
    char fnameT[50]="_temp.txt";

    FILE *fin,*fout;

    for(i=0;i<=NDmax;i++){datax[i]=0.0;datay[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;
    /*========================*/
    /* Sorting in small order */
    /*========================*/
    for(i=1;i<=nd;i++){
        for(j=i+1;j<=nd;j++){
            if(datax[i]>datax[j]){
                work=datax[i];datax[i]=datax[j];datax[j]=work;
            }
        }
    }
    /*----------------------------------*/
    /* Non-exceeding probability        */
    /* Plotting and positioned formula  */
    /*    Weibull : alpha=0             */
    /*    Cunnane : alpha=0.4           */
    /*    Hazen   : alpha=0.5           */
    /*----------------------------------*/
    alpha=0.4;
    for(i=1;i<=nd;i++){datay[i]=(i-alpha)/(nd+1-2.0*alpha);}
    /*=============================*/
    /* Calculation of coefficients */
    /*=============================*/
    GUM_COEF(nd,datax,&a_gum,&c_gum);
    GEV_COEF(nd,datax,&k_gev,&a_gev,&c_gev);
    SQR_COEF(nd,datax,&a_sqr,&b_sqr);
    LN3_COEF(nd,datax,&a_ln3,&my_ln3,&sy_ln3);
    LP3_COEF(nd,datax,&a_lp3,&b_lp3,&c_lp3);
    /*===============================*/
    /* Calculation of quantile value */
    /*===============================*/
    fout=fopen(fnameT,"w");
    for(i=1;i<=nd;i++){
        p=datay[i];
        yy1=XP_GUM(p,a_gum,c_gum);
        yy2=XP_GEV(p,k_gev,a_gev,c_gev);
        yy3=XP_SQR(p,a_sqr,b_sqr,datax[nd]);
        yy4=XP_LN3(p,a_ln3,my_ln3,sy_ln3);
        if(b_lp3<1e4){
            yy5=XP_LP3(p,a_lp3,b_lp3,c_lp3);
        }else{
            yy5=XP_LP3_SP(nd,datax,p);
        }
        fprintf(fout,"%.4f,%.1f,%.1f,%.1f,%.1f,%.1f,%.1f\n",p,datax[i],yy1,yy2,yy3,yy4,yy5);
        printf("%.4f,%.1f,%.1f,%.1f,%.1f,%.1f,%.1f\n",p,datax[i],yy1,yy2,yy3,yy4,yy5);
    }
    fclose(fout);
    /*============================*/
    /* Calculation of correlation */
    /*============================*/
    for(nnn=1;nnn<=5;nnn++){
        fin=fopen(fnameT,"r");
        i=0;
        while(fgets(dat,sizeof dat,fin)!=NULL){
            i=i+1;
            work=atof(strtok(dat,","));
            datax[i]=atof(strtok(NULL,","));
            yy1=atof(strtok(NULL,","));
            yy2=atof(strtok(NULL,","));
            yy3=atof(strtok(NULL,","));
            yy4=atof(strtok(NULL,","));
            yy5=atof(strtok(NULL,","));
            switch(nnn){
                case 1:datay[i]=yy1;break;
                case 2:datay[i]=yy2;break;
                case 3:datay[i]=yy3;break;
                case 4:datay[i]=yy4;break;
                case 5:datay[i]=yy5;break;
            }
        }
        fclose(fin);
        nd=i;
        switch(nnn){
            case 1:rr1=COC(nd,datax,datay);break;
            case 2:rr2=COC(nd,datax,datay);break;
            case 3:rr3=COC(nd,datax,datay);break;
            case 4:rr4=COC(nd,datax,datay);break;
            case 5:rr5=COC(nd,datax,datay);break;
        }
    }
    fout=fopen(fnameW,"w");
    fprintf(fout,"Input file : %s\n",fnameR);
    fprintf(fout,"Gumbel(a,c)=(%.4f,%.4f)\n",a_gum,c_gum);
    fprintf(fout,"GEV(k,a,c)=(%.4f,%.4f %.4f)\n",k_gev,a_gev,c_gev);
    fprintf(fout,"SQRT-ET(a,b)=(%.4f,%.4f)\n",a_sqr,b_sqr);
    fprintf(fout,"LN3(a,mu,sig)=(%.4f,%.4f,%.4f)\n",a_ln3,my_ln3,sy_ln3);
    fprintf(fout,"LP3(a,b,c)=(%.4f,%.4f,%.4f)\n",a_lp3,b_lp3,c_lp3);
    fprintf(fout,"Prob,Obs,Gumbel,GEV,SQRT-ET,LN3,LP3\n");
    fprintf(fout,"rr,%d,%.3f,%.3f,%.3f,%.3f,%.3f\n",nd,rr1,rr2,rr3,rr4,rr5);
        fin=fopen(fnameT,"r");
        while(fgets(dat,sizeof dat,fin)!=NULL){fprintf(fout,"%s",dat);}
        fclose(fin);
    fclose(fout);

}
/*---------------------------------------------------------------------------*/
double COC(int nd,double datax[NDmax],double datay[NDmax])
{
    /*Coefficient of correlation*/
    int i;
    double x1,y1,x2,y2,xy,xm,ym,sx,sy,cv,rr;

    x1=0.0;y1=0.0;x2=0.0;y2=0.0;xy=0.0;
    for(i=1;i<=nd;i++){
        x1=x1+datax[i];
        y1=y1+datay[i];
        x2=x2+datax[i]*datax[i];
        y2=y2+datay[i]*datay[i];
        xy=xy+datax[i]*datay[i];
    }
    xm=x1/(double)nd;
    ym=y1/(double)nd;
    sx=sqrt(x2/(double)nd-xm*xm);
    sy=sqrt(y2/(double)nd-ym*ym);
    cv=xy/(double)nd-xm*ym;

    return cv/sx/sy;
}
/*---------------------------------------------------------------------------*/
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);
}
/*-----------------------------------------------------------------------*/