#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#define Nmax 100
#define Mmax 20
void main_part(char fnameR[],char fnameW[]);
void JACOBI(int n,double a[Mmax][Mmax],double ev[Mmax],double vec[Mmax][Mmax]);
/*-----------------------------------------------------------------------*/

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

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

    main_part(fnameR,fnameW);

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

void main_part(char fnameR[],char fnameW[])
{
    /* 主成分分析（principal component analysis）*/
    FILE *fin,*fout;
    char dat[256];
    char sv[20];
    double xd[Nmax+1][Mmax]; /* データ入力用配列 */
    double a[Mmax][Mmax];       /* 分散共分散行列格納配列 */
    double xm[Mmax];          /* 平均値格納配列 */
    double sd[Mmax];          /* 標準偏差格納配列 */
    double ev[Mmax];          /* 固有値 */
    double vec[Mmax][Mmax];     /* 固有ベクトル */
    double zd[Nmax+1][Mmax]; /* 主成分得点（標準化変数×固有ベクトル） */
    
    int i,j,k,nn;
    char scom[256]; /*コメント */
    int mm;         /*変数の数−１ */
    int ndata;      /*データ組数−１ */
    double tr;      /*固有値総和（分散共分散行列対角項総和） */

    /* データ読み込み */
    fin=fopen(fnameR,"r");
    fgets(dat,sizeof dat,fin);strcpy(scom,strtok(dat,","));
    fgets(dat,sizeof dat,fin);mm=atoi(strtok(dat,","));ndata=atoi(strtok(NULL,","))-1;
    i=0;
    while(fgets(dat,sizeof dat,fin)!=NULL){
        xd[i][0]=atof(strtok(dat,","));
        for(j=1;j<=mm-1;j++){
            xd[i][j]=atof(strtok(NULL,","));
        }
        i=i+1;
    }
    fclose(fin);
    ndata=i-1;

    /* 平均値の計算 */
    for(i=0;i<=mm-1;i++){
        xm[i]=0.0;
            for(k=0;k<=ndata;k++){
                xm[i]=xm[i]+xd[k][i];
            }
            xm[i]=xm[i]/(double)(ndata+1);
    }
    /* 標準偏差の計算 */
    for(i=0;i<=mm-1;i++){
        sd[i]=0.0;
        for(k=0;k<=ndata;k++){
            sd[i]=sd[i]+(xd[k][i]-xm[i])*(xd[k][i]-xm[i]);
        }
        sd[i]=sqrt(sd[i]/(double)(ndata));
    }
    /* 変数の標準化 */
    for(i=0;i<=mm-1;i++){
        for(k=0;k<=ndata;k++){
            xd[k][i]=(xd[k][i]-xm[i])/sd[i];
        }
    }
    /* 分散共分散行列(相関行列)の計算 */
    tr=0.0;
    for(i=0;i<=mm-1;i++){
        for(j=0;j<=mm-1;j++){
            a[i][j]=0.0;
            for(k=0;k<=ndata;k++){
                a[i][j]=a[i][j]+xd[k][i]*xd[k][j];
            }
            a[i][j]=a[i][j]/(double)ndata;
        }
        tr=tr+a[i][i]; /* 対角項の総和（=固有値の総和） */
    }
    /* 固有値・固有ベクトルの計算 */
    JACOBI(mm,a,ev,vec);
    /* 主成分得点の計算 */
    for(k=0;k<=ndata;k++){
        for(i=0;i<=mm-1;i++){
            zd[k][i]=0.0;
            for(j=0;j<=mm-1;j++){
                zd[k][i]=zd[k][i]+xd[k][j]*vec[j][i];
            }
        }
    }

    /* 画面書き出し */
    nn=4;
    if(nn<mm)nn=mm;
    printf("%s",scom);
    strcpy(dat,"No.");
    for(i=0;i<=nn-1;i++){
        sprintf(sv," %d",i+1);strcat(dat,",");strcat(dat,sv);
    }
    printf("%s\n",dat);
    strcpy(dat,"固有値");
    for(i=0;i<=nn-1;i++){
        sprintf(sv," %e",ev[i]);strcat(dat,",");strcat(dat,sv);
    }
    printf("%s\n",dat);
    strcpy(dat,"寄与率");
    for(i=0;i<=nn-1;i++){
        sprintf(sv," %e",ev[i]/tr);strcat(dat,",");strcat(dat,sv);
    }
    printf("%s\n",dat);
    for(j=0;j<=mm-1;j++){
        sprintf(sv,"[%d]",j+1);
        strcpy(dat,"固有ベクトル");strcat(dat,sv);
        for(i=0;i<=nn-1;i++){
            sprintf(sv," %e",vec[j][i]);strcat(dat,",");strcat(dat,sv);
        }
        printf("%s\n",dat);
    }

    /*結果書き出し*/
    fout=fopen(fnameW,"w");
    fprintf(fout,"%s",scom);
    strcpy(dat,"No.");
    for(i=0;i<=mm-1;i++){
        sprintf(sv," %d",i+1);strcat(dat,",");strcat(dat,sv);
    }
    fprintf(fout,"%s\n",dat);
    strcpy(dat,"固有値");
    for(i=0;i<=mm-1;i++){
        sprintf(sv," %e",ev[i]);strcat(dat,",");strcat(dat,sv);
    }
    fprintf(fout,"%s\n",dat);
    strcpy(dat,"寄与率");
    for(i=0;i<=mm-1;i++){
        sprintf(sv," %e",ev[i]/tr);strcat(dat,",");strcat(dat,sv);
    }
    fprintf(fout,"%s\n",dat);
    for(j=0;j<=mm-1;j++){
        sprintf(sv,"[%d]",j+1);
        strcpy(dat,"固有ベクトル");
        strcat(dat,sv);
        for(i=0;i<=mm-1;i++){
            sprintf(sv," %e",vec[j][i]);strcat(dat,",");strcat(dat,sv);
        }
        fprintf(fout,"%s\n",dat);
    }
    for(i=0;i<=ndata;i++){
        sprintf(sv,"[%d]",i+1);
        strcpy(dat,"data");
        strcat(dat,sv);
        for(j=0;j<=mm-1;j++){
            sprintf(sv," %e",zd[i][j]);strcat(dat,",");strcat(dat,sv);
        }
        fprintf(fout,"%s\n",dat);
    }
    fclose(fout);

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

void JACOBI(int n,double a[Mmax][Mmax],double ev[Mmax],double vec[Mmax][Mmax])
{
    /* Jacobi法による固有値・固有ベクトル計算（奥村先生）*/
    /* 固有値ev[k],固有ベクトルvec[i][k] */

    int i,j,k,im=0;
    double s,c,aoff,eps,t,aa=0.0,aj,ak,vj,vk,vmax,bb;

    for(i=0;i<=n-1;i++){
        for(j=0;j<=n-1;j++){
            vec[i][j]=0.0;
        }
        vec[i][i]=1.0;
    }
    aoff=0.0;
    for(j=0;j<=n-2;j++){
        for(k=j+1;k<=n-1;k++){
            aoff=aoff+a[j][k]*a[j][k];
        }
    }
    eps=aoff*2.0;
    for(j=0;j<=n-1;j++){
        eps=eps+a[j][j]*a[j][j];
    }
    eps=eps*0.00000000005;

    while(aoff>eps){
        for(j=0;j<=n-2;j++){
            for(k=j+1;k<=n-1;k++){
                if(a[j][k]==0.0){
                    t=0.0;
                }else{
                    aa=(a[k][k]-a[j][j])/(2.0*a[j][k]);
                }
                if(aa>=0){
                    t=1.0/(aa+sqrt(1.0+aa*aa));
                }else{
                    t=1.0/(aa-sqrt(1.0+aa*aa));
                }
                c=1.0/sqrt(1.0+t*t);
                s=t*c;
                a[j][j]=a[j][j]-t*a[j][k];
                a[k][k]=a[k][k]+t*a[j][k];
                a[j][k]=0.0;
                for(i=0;i<=j-1;i++){
                    aj=a[i][j];
                    ak=a[i][k];
                    a[i][j]=aj*c-ak*s;
                    a[i][k]=aj*s+ak*c;
                }
                for(i=j+1;i<=k-1;i++){
                    aj=a[j][i];
                    ak=a[i][k];
                    a[j][i]=aj*c-ak*s;
                    a[i][k]=aj*s+ak*c;
                }
                for(i=k+1;i<=n-1;i++){
                    aj=a[j][i];
                    ak=a[k][i];
                    a[j][i]=aj*c-ak*s;
                    a[k][i]=aj*s+ak*c;
                }
                for(i=0;i<=n-1;i++){
                    vj=vec[i][j];
                    vk=vec[i][k];
                    vec[i][j]=vj*c-vk*s;
                    vec[i][k]=vj*s+vk*c;
                }
            }
        }
        aoff=0.0;
        for(j=0;j<=n-2;j++){
            for(k=j+1;k<=n-1;k++){
                aoff=aoff+a[j][k]*a[j][k];
            }
        }
    }
            
    for(k=0;k<=n-1;k++){
        ev[k]=a[k][k];
    }
    for(k=0;k<=n-2;k++){
        vmax=-1.0E+30;
        for(i=k;i<=n-1;i++){
            if(ev[i]>=vmax){
                im=i;
                vmax=ev[i];
            }
        }
        aa=ev[k];ev[k]=ev[im];ev[im]=aa;
        for (i=0;i<=n-1;i++){
            bb=vec[i][k];
            vec[i][k]=vec[i][im];
            vec[i][im]=bb;
        }
    }

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