/***************************************************/
/* gccWEI.c */
/***************************************************/
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#define NDmax 5000
void main_part(char fnameR[50],char fnameW[50]);
double COC(int nd,double datax[NDmax],double datay[NDmax]);
void WLM_COEF(int nd,double datax[NDmax],double *k_wlm,double *a_wlm,double *c_wlm);
void WML_COEF(int nd,double datax[NDmax],double datay[NDmax],double *k_t,double *a_t,double *c_t,double *k_wml,double *a_wml,double *c_wml);
double XP_wei(double p,double k_wei,double a_wei,double c_wei);
double GAMMAF(double x);
void RSM(int nd,double vecx[NDmax],double vecy[NDmax],double *aa,double *bb,double *rr);
/*---------------------------------------------------------------------------*/

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],datay[NDmax];
    double alpha;
    double p,yy1,yy2,yy3,yy4,yy5,work;
    double rr1,rr2,rr3,rr4,rr5;
    double a_wlm,c_wlm,k_wlm;
    double a_wrs,c_wrs,k_wrs;
    double a_wml,c_wml,k_wml;
    char dat[256]="";
    char fnameT[50]="_temp.txt";

    FILE *fin,*fout;

    /*============*/
    /* 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 (Mean Rank)         */
    /*   (Median Rank : alpha=0.3 )             */
    /*    Cunnane : alpha=0.4                   */
    /*    Hazen   : alpha=0.5 (Symmetrical CDF) */
    /*------------------------------------------*/
    alpha=0.4;
    for(i=1;i<=nd;i++){datay[i]=(i-alpha)/(nd+1-2.0*alpha);}
    /*=============================*/
    /* Calculation of coefficients */
    /*=============================*/
    WLM_COEF(nd,datax,&k_wlm,&a_wlm,&c_wlm);
    WML_COEF(nd,datax,datay,&k_wrs,&a_wrs,&c_wrs,&k_wml,&a_wml,&c_wml);
    /*===============================*/
    /* Calculation of quantile value */
    /*===============================*/
    fout=fopen(fnameT,"w");
    for(i=1;i<=nd;i++){
        p=datay[i];
        yy1=XP_wei(p,k_wlm,a_wlm,c_wlm);
        yy2=XP_wei(p,k_wrs,a_wrs,c_wrs);
        yy3=XP_wei(p,k_wml,a_wml,c_wml);
        fprintf(fout,"%.4f,%.1f,%.1f,%.1f,%.1f\n",p,datax[i],yy1,yy2,yy3);
    }
    fclose(fout);
    /*============================*/
    /* Calculation of correlation */
    /*============================*/
    for(nnn=1;nnn<=3;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,","));
            switch(nnn){
                case 1:datay[i]=yy1;break;
                case 2:datay[i]=yy2;break;
                case 3:datay[i]=yy3;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;
        }
    }
    fout=fopen(fnameW,"w");
    fprintf(fout,"Input file : %s\n",fnameR);
    fprintf(fout,"LMM(k,a,c)=(%.4f,%.4f %.4f)\n",k_wlm,a_wlm,c_wlm);
    fprintf(fout,"RSM(k,a,c)=(%.4f,%.4f %.4f)\n",k_wrs,a_wrs,c_wrs);
    fprintf(fout,"MLM(k,a,c)=(%.4f,%.4f %.4f)\n",k_wml,a_wml,c_wml);
    fprintf(fout,"rr,%d,%.3f,%.3f,%.3f\n",nd,rr1,rr2,rr3);
        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 WLM_COEF(int nd,double datax[NDmax],double *k_wlm,double *a_wlm,double *c_wlm)
{
    /* Generalized Extreme Value distribution  */
    int i;
    double b0,b1,b2,t;
    double lam1,lam2,lam3;

    b0=0;b1=0;b2=0;
    for(i=1;i<=nd;i++){
        b0=b0+datax[i];
        b1=b1+(double)(i-1)*datax[i];
        b2=b2+(double)((i-1)*(i-2))*datax[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;

    /*Weibull F(x)=1-exp{-[(x-c)/a]^k} T=1/[1-F(x)] */
    t=lam3/lam2;
    *k_wlm=285.3*pow(t,6.0)-658.6*pow(t,5.0)+622.8*pow(t,4.0)-317.2*pow(t,3.0)+98.52*pow(t,2.0)-21.256*t+3.5160;
    *a_wlm=lam2/(1.0-pow(2.0,-1.0/ *k_wlm))/GAMMAF(1.0+1.0/ *k_wlm);
    *c_wlm=lam1- *a_wlm*GAMMAF(1.0+1.0/ *k_wlm);
}
/*---------------------------------------------------------------------------*/
double XP_wei(double p,double k_wei,double a_wei,double c_wei)
{
    return c_wei+a_wei*pow(-log(1.0-p),1.0/k_wei);
}
/*---------------------------------------------------------------------------*/
double GAMMAF(double x)
{
    double a,b,p,q;

    b=x;
    a=1.0/(x*(x+1.0));
    b=b-1.0;
    if(0.0<=b){
        do{
            a=a*(b+2.0);
            b=b-1.0;
        }while(0.0<=b);
    }
    b=x-(int)x;
    p=(((-45.30104765151624*b-261.0284950616458)*b-1817.799191156434)*b-4954.768972677267)*b-13572.65011103686;
    q=((((((1.0*b-16.16160629773462)*b+93.7769173224393)*b-124.7947256227383)*b-931.1785045593485)*b+3440.699241346454)*b+783.5348800559552)*b-13572.65011103691;
    return a*p/q;
}
/*---------------------------------------------------------------------------*/
void WML_COEF(int nd,double datax[NDmax],double datay[NDmax],double *k_t,double *a_t,double *c_t,double *k_wml,double *a_wml,double *c_wml)
{
    int i;
    double c,vecx[NDmax],vecy[NDmax];
    double aa,bb,rr1,rr2;
    double kk,ff,df,tt,t0,t1,t2,t3;

    /* Culculation of coefficient (c) */
    c=datax[1];rr1=0.0;
    do{
        c=c-0.001;
        for(i=1;i<=nd;i++){
            vecx[i]=log(datax[i]-c);
            vecy[i]=log(-log(1.0-datay[i]));
        }
        RSM(nd,vecx,vecy,&aa,&bb,&rr2);
        if(fabs(rr2)<fabs(rr1))break;
        rr1=rr2;
    }while(0.0<=c);
    *k_t=aa;
    *a_t=exp(-bb/aa);
    *c_t=c;

    /* Newton-Raphson method */
    kk=*k_t;
    do{
        t0=0.0;t1=0.0;t2=0.0;t3=0.0;
        for(i=1;i<=nd;i++){
            tt=datax[i]- *c_t;
            t0=t0+log(tt);
            t1=t1+pow(tt,kk);
            t2=t2+log(tt)*pow(tt,kk);
            t3=t3+log(tt)*log(tt)*pow(tt,kk);
        }
        ff=1.0/kk+t0/(double)nd-t2/t1;
        df=-1.0/kk/kk-(t3*t1-t2*t2)/t1/t1;
        kk=kk-ff/df;
    }while(1e-5<fabs(ff/df));

    *k_wml=kk;
    *a_wml=pow(t1/(double)nd,1.0/kk);
    *c_wml=*c_t;
}
/*---------------------------------------------------------------------------*/
void RSM(int nd,double vecx[NDmax],double vecy[NDmax],double *aa,double *bb,double *rr)
{
    /*--------------------------*/
    /* Regression (y=aa*x+bb)   */
    /*--------------------------*/
    int i;
    double x1,y1,x2,y2,xy,xm,ym,c1,c2,c3;
    
    x1=0.0;y1=0.0;x2=0.0;xy=0.0;
    for(i=1;i<=nd;i++){
        x1=x1+vecx[i];
        y1=y1+vecy[i];
        x2=x2+vecx[i] * vecx[i];
        xy=xy+vecx[i] * vecy[i];
    }
    xm=x1/(double)nd;ym=y1/(double)nd;
    *aa=((double)nd*xy-x1*y1)/((double)nd*x2-x1*x1);
    *bb=(x2*y1-x1*xy)/((double)nd*x2-x1*x1);
    c1=0.0;c2=0.0;c3=0.0;
    for(i=1;i<=nd;i++){
        c1=c1+(vecx[i]-xm)*(vecy[i]-ym);
        c2=c2+(vecx[i]-xm)*(vecx[i]-xm);
        c3=c3+(vecy[i]-ym)*(vecy[i]-ym);
    }
    *rr=c1/sqrt(c2*c3);
}
/*---------------------------------------------------------------------------*/