/***************************************************/
/* gccSIMPLEX.c */
/***************************************************/
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#define MPmax 5
#define NDmax 5000
/*gnuplot用構造体宣言*/
/*scom,sxjiku,syjiku :タイトル(凡例),x軸ラベル,ｙ軸ラベル*/
/*slogx,xmin,xmax,dx,sxfmt :x軸種,x軸最小値,x軸最大値,x軸増分,x軸書式*/
/*slogy,ymin,ymax,dy,syfmt :y軸種,y軸最小値,y軸最大値,y軸増分,y軸書式*/
struct GPdata{
    char scom[100];char sxjiku[100];char syjiku[100];
    char slogx[3];double xmin;double xmax;double dx;char sxfmt[20];
    char slogy[3];double ymin;double ymax;double dy;char syfmt[20];
};
/*関数プロトタイプ*/
int MPARA_SET();
double Vfunc(double x,double parameter[MPmax]);
void INIVAL(int mpara,double vertex[MPmax+1][MPmax],
    double stepsize[MPmax],double tolerance[MPmax]);
void main_part(char fnameR[50],char fnameW[50],char fnameP[50]);
void SIMPLEX(int mpara,double parameter[MPmax],
    int ndata,double datax[NDmax],double datay[NDmax],
    int *iteration,int *iflag);
void GPLOT(char fnameP[50],struct GPdata gp);
/*---------------------------------------------------------------------------*/

int MPARA_SET()
{
    /***********************************************/
    /*パラメータ数：mparaのセット*/
    /***********************************************/
    int mpara;

    /*mpara=5;*/
    mpara=4;

    return mpara;
}

/*---------------------------------------------------------------------------*/

double Vfunc(double x,double parameter[MPmax])
{
    /***********************************************/
    /*関数値Vfuncの定義・計算*/
    /*parameter[i]：回帰係数(i=0〜mpara-1)*/
    /***********************************************/
    double fx;

    /*パソコンによるデータ解析入門p236(mpara=5)*/
    /*
    double gamma2,dx;
    gamma2=parameter[1]*parameter[1];
    dx=x-parameter[0];
    fx=parameter[2]*gamma2/(dx*dx+gamma2)+parameter[3]+parameter[4]*x;
    */

    /*減衰振動(mpara=4)*/
    fx=parameter[0]*exp(-parameter[1]*x)*sin(parameter[2]*x+parameter[3]);
    
    return fx;
}
/*---------------------------------------------------------------------------*/

void INIVAL(int mpara,double vertex[MPmax+1][MPmax],
    double stepsize[MPmax],double tolerance[MPmax])
{
    /********************************************/
    /*解析ﾊﾟﾗﾒ-ﾀ初期値等入力(0〜mpara-1)        */
    /*vertex：初期値                            */
    /*stepsize：初期の検索ステップ              */
    /*tolerance：収束判定基準                   */
    /********************************************/
    int j;
    for(j=0;j<=mpara-1;j++){
        vertex[0][j]=1.0;
        stepsize[j]=0.1*vertex[0][j];
        tolerance[j]=0.000001;
    }
}
/*---------------------------------------------------------------------------*/

int main(int argc,char *argv[])
{
    char fnameR[50]; /*入力ファイル名*/
    char fnameW[50]; /*出力ファイル名*/
    char fnameP[50]; /*プロットファイル名*/

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

    main_part(fnameR,fnameW,fnameP);

    printf("リターンキーを押してください");
    getchar();
    return 0;
}
/*---------------------------------------------------------------------------*/

void main_part(char fnameR[50],char fnameW[50],char fnameP[50])
{
    int i;
    int ndata,mpara;
    int iflag,iteration;
    double parameter[MPmax];
    double datax[NDmax],datay[NDmax];
    double x,y;
    char dat[256]="";
    char sv[10]="";
    struct GPdata gp={
        "","","",
        "",0.0,0.0,0.0,"",
        "",0.0,0.0,0.0,""
    };
    FILE *fin,*fout;

    /**********************/
    /*パラメータ数定義*/
    /**********************/
    mpara=MPARA_SET();

    /*データ入力*/
    fin=fopen(fnameR,"r");
    fgets(dat,sizeof dat,fin);
        strcpy(sv,strtok(dat,","));
        strncat(gp.scom,sv,strlen(sv)-1);
    fgets(dat,sizeof dat,fin);
        strcpy(sv,strtok(dat,","));
        strncat(gp.sxjiku,sv,strlen(sv));
        strcpy(sv,strtok(NULL,","));
        strncat(gp.syjiku,sv,strlen(sv)-1);
    fgets(dat,sizeof dat,fin);
        strcpy(gp.slogx,strtok(dat,","));
        gp.xmin=atof(strtok(NULL,","));
        gp.xmax=atof(strtok(NULL,","));
        if(strncmp(gp.slogx,"N",1)==0){
            gp.dx=atof(strtok(NULL,","));
            strcpy(sv,strtok(NULL,","));
            strncat(gp.sxfmt,sv,strlen(sv)-1);
        }
    fgets(dat,sizeof dat,fin);
        strcpy(gp.slogy,strtok(dat,","));
        gp.ymin=atof(strtok(NULL,","));
        gp.ymax=atof(strtok(NULL,","));
        if(strncmp(gp.slogy,"N",1)==0){
            gp.dy=atof(strtok(NULL,","));
            strcpy(sv,strtok(NULL,","));
            strncat(gp.syfmt,sv,strlen(sv)-1);
        }
    fgets(dat,sizeof dat,fin);
        ndata=atof(strtok(dat,","))-1;
        printf("ndata=%d\n",ndata+1);
    for(i=0;i<=ndata;i++){
        fgets(dat,sizeof dat,fin);
            datax[i]=atof(strtok(dat,","));
            datay[i]=atof(strtok(NULL,","));
            printf("%d  %g  %g\n",i,datax[i],datay[i]);
    }
    fclose(fin);

    /*回帰係数計算*/
    SIMPLEX(mpara,parameter,ndata,datax,datay,&iteration,&iflag);

    /*結果出力*/
    fout=fopen(fnameW,"w");
        strcpy(dat,"iflag,");sprintf(sv,"%d",iflag);strcat(dat,sv);
        strcat(dat,"\n");
        fputs(dat,fout);
        strcpy(dat,"iteration,");sprintf(sv,"%d",iteration);strcat(dat,sv);
        strcat(dat,"\n");
        fputs(dat,fout);
        for(i=0;i<=mpara-1;i++){
            strcpy(dat,"parameter");
            sprintf(sv,"%d",i);strcat(dat,"_");strcat(dat,sv);
            sprintf(sv,"%g",parameter[i]);strcat(dat,",");strcat(dat,sv);
            strcat(dat,"\n");
            fputs(dat,fout);
        }
        strcpy(dat,"datax,datay,y_est\n");
        fputs(dat,fout);
        for(i=0;i<=ndata;i++){
            sprintf(sv,"%g",datax[i]);strcpy(dat,sv);
            sprintf(sv,"%g",datay[i]);strcat(dat,",");strcat(dat,sv);
            sprintf(sv,"%g",Vfunc(datax[i],parameter));strcat(dat,",");strcat(dat,sv);
            strcat(dat,"\n");
            fputs(dat,fout);
        }
    fclose(fout);

    /*入力データプロット用一時ファイル書き出し*/
    fout=fopen("temp0.prn","w");
    for(i=0;i<=ndata;i++){
        fprintf(fout,"%f %f\n",datax[i],datay[i]);
    }
    fclose(fout);

    /*回帰曲線プロット用一時ファイル書き出し*/
    fout=fopen("temp1.prn","w");
    for(i=0;i<=NDmax-1;i++){
        x=gp.xmin+(gp.xmax-gp.xmin)/(double)(NDmax-1)*(double)i;
        y=Vfunc(x,parameter);
        fprintf(fout,"%f %f\n",x,y);
    }
    fclose(fout);

    /*gnuplotによる作図*/
    GPLOT(fnameP,gp);

}
/*---------------------------------------------------------------------------*/

void SIMPLEX(int mpara,double parameter[MPmax],
    int ndata,double datax[NDmax],double datay[NDmax],
    int *iteration,int *iflag)
{
    double vertex[MPmax+1][MPmax];
    double criterion[MPmax+1];
    double stepsize[MPmax];
    double tolerance[MPmax];
    double s,v1,v2,cr,v,vmax,vmin,vsum,p,x,y;
    int i,j,k,kworst,kbest,ks;

    /********************************************/
    /*解析ﾊﾟﾗﾒ-ﾀ初期値等入力*/
    /********************************************/
    INIVAL(mpara,vertex,stepsize,tolerance);

    *iflag=0;
    ks=0;
    s=sqrt(0.5);
    v1=s*(sqrt((double)(mpara+1))-1)/(double)mpara;
    for (j=0;j<=mpara-1; j++){
        v2=vertex[0][j]+v1*stepsize[j];
        for(k=1;k<=mpara;k++){
            vertex[k][j]=v2;
        }
        vertex[j+1][j]=v2+s*stepsize[j];
    }
    for(k=0;k<=mpara;k++){
        for (j=0;j<=mpara-1;j++){
            parameter[j]=vertex[k][j];
        }
        cr=0; /*Vcrit-Start*/
        for(i=0;i<=ndata;i++){
            x=datax[i];
            y=Vfunc(x,parameter);
            cr=cr+(y-datay[i])*(y-datay[i]);
        }     /*Vcrit-End*/
        criterion[k]=cr;
    }
    *iteration=0;
    while(*iteration<=100000){
        *iteration=*iteration+1;
        kworst=0;
        kbest=0;
        for(k=1;k<=mpara;k++){
            if(criterion[k]>criterion[kworst])kworst=k;
            if(criterion[k]<criterion[kbest])kbest=k;
        }
        *iflag=0;
        for(j=0;j<=mpara-1;j++){
            vmax=vertex[0][j];
            vmin=vmax;
            for(k=1;k<=mpara;k++){
                v=vertex[k][j];
                if(v>=vmax)vmax=v;
                if(v<vmin)vmin=v;
            }
            if(vmax-vmin>tolerance[j])*iflag=1;
        }
        if(*iflag==0)break;

        for(j=0;j<=mpara-1;j++){
            vsum=0;
            for(k=0;k<=mpara;k++){
                if(k!=kworst)vsum=vsum+vertex[k][j];
            }
            parameter[j]=2.0*vsum/(double)mpara-vertex[kworst][j];
        }
        cr=0; /*Vcrit-Start*/
        for(i=0;i<=ndata;i++){
            x=datax[i];
            y=Vfunc(x,parameter);
            cr=cr+(y-datay[i])*(y-datay[i]);
        }     /*Vcrit-End*/
        if(cr<criterion[kbest])ks=0;
        if(criterion[kbest]<=cr&&cr<=criterion[kworst])ks=1;
        if(criterion[kworst]<cr)ks=2;
        switch(ks){
            case 0:
                for(j=0;j<=mpara-1;j++){
                    p=1.5*parameter[j]-0.5*vertex[kworst][j];
                    vertex[kworst][j]=parameter[j];
                    parameter[j]=p;
                }
                criterion[kworst]=cr;
                cr = 0; /*Vcrit-Start*/
                for(i=0;i<=ndata;i++){
                    x=datax[i];
                    y=Vfunc(x,parameter);
                    cr=cr+(y-datay[i])*(y-datay[i]);
                }     /*Vcrit-End*/
                if(cr<criterion[kworst]){
                    for(j=0;j<=mpara-1;j++){  /*REPLACE-Start*/
                        vertex[kworst][j] = parameter[j];
                    }
                    criterion[kworst]=cr; /*REPLACE-End*/
                }
                break;
            case 1:
                for(j=0;j<=mpara-1;j++){  /*REPLACE-Start*/
                    vertex[kworst][j]=parameter[j];
                }
                criterion[kworst]=cr; /*REPLACE-End*/
                break;
            case 2:
                for(j=0;j<=mpara-1;j++){
                    parameter[j]=0.75*vertex[kworst][j]+0.25*parameter[j];
                }
                cr=0; /*Vcrit-Start*/
                for(i=0;i<=ndata;i++){
                    x=datax[i];
                    y=Vfunc(x,parameter);
                    cr=cr+(y-datay[i])*(y-datay[i]);
                }     /*Vcrit-End*/
                if(cr<criterion[kworst]){
                    for(j=0;j<=mpara-1;j++){  /*REPLACE-Start*/
                        vertex[kworst][j]=parameter[j];
                    }
                    criterion[kworst]=cr; /*REPLACE-End*/
                }else{
                    for(k=0;k<=mpara;k++){
                        if(k!=kbest){
                            for(j=0;j<=mpara-1;j++){
                                parameter[j]=0.5*(vertex[k][j]+vertex[kbest][j]);
                                vertex[k][j]=parameter[j];
                            }
                            cr=0; /*Vcrit-Start*/
                            for(i=0;i<=ndata;i++){
                                x=datax[i];
                                y=Vfunc(x,parameter);
                                cr=cr+(y-datay[i])*(y-datay[i]);
                            }     /*Vcrit-End*/
                            criterion[k]=cr;
                        }
                    }
                }
                break;
        } /*switch*/
    } /*while*/
}
/*---------------------------------------------------------------------------*/

void GPLOT(char fnameP[50],struct GPdata gp)
{
    FILE *fp;

    /****************************/
    /*gnuplot制御用ファイル作成*/
    /****************************/
    fp=fopen("gpl_temp.txt","w");
        /*出力ファイル定義*/
        fprintf(fp,"set terminal png\n");
        fprintf(fp,"set output \"%s\"\n",fnameP);
        /*軸ラベル定義*/
        fprintf(fp,"set xlabel \"%s\"\n",gp.sxjiku);
        fprintf(fp,"set ylabel \"%s\"\n",gp.syjiku);
        /*ｘ軸定義*/
        if(strncmp(gp.slogx,"L",1)==0){
            fprintf(fp,"set logscale x\n");
            fprintf(fp,"set mxtics 10\n");
        }
        if(strncmp(gp.slogx,"N",1)==0){
            fprintf(fp,"set xtics %f\n",gp.dx);
            fprintf(fp,"set format x %s\n",gp.sxfmt);
        }
        /*ｙ軸定義*/
        if(strncmp(gp.slogy,"L",1)==0){
            fprintf(fp,"set logscale y\n");
            fprintf(fp,"set mytics 10\n");
        }
        if(strncmp(gp.slogy,"N",1)==0){
            fprintf(fp,"set ytics %f\n",gp.dy);
            fprintf(fp,"set format y %s\n",gp.syfmt);
        }
        /*軸範囲定義*/
        fprintf(fp,"set xrange [%f:%f]\n",gp.xmin,gp.xmax);
        fprintf(fp,"set yrange [%f:%f]\n",gp.ymin,gp.ymax);
        /*軸グリッド描画設定*/
        fprintf(fp,"set grid xtics ytics mxtics mytics\n");
        /*プロット実行（プロットは各種定義が終わって最後に実行）*/
        fprintf(fp,"plot \\\n");
        fprintf(fp,"\"temp0.prn\" with points pointtype 1 title \"Input data\", \\\n");
        fprintf(fp,"\"temp1.prn\" with lines linetype 3 title \"Estimated\"\n");
    fclose(fp);

    /****************************/
    /*gnuplot制御*/
    /****************************/
    fp=popen("gnuplot.exe","w");
        fprintf(fp,"load \"gpl_temp.txt\"\n");
        fflush(fp);
    pclose(fp);
}
/*---------------------------------------------------------------------------*/