/*
 *  Accuracy assessment for Error Matrices at most 20 x 20
 */


#include <stdio.h>
#include <string.h>
#include <math.h>
#define M 20  /* Maximum size of Error matrices */
#define S 400 /* = M^2                          */

     FILE     *fopen(),*data_r1,*data_r2,*data_w;
     char     *f_name[4],tmp1[10],tmp2[10],tmp3[10],cc[3];
              /* Variables used for files */

     int      c,num,i,j,r,r2,X[M][M],tmp[S];
     long int xc[M],xr[M],N,ni,sum_a;

     double   G1[M][M],G2[M][M],
              Pai[M],CX,AX,KapX,CBX[M],sum_c,K_hat,
              J_u,J_p,Ja_u,Ja_p,R_u,R_p,Tau_u[3],Tau_p[3],p[M],
              Var_u,Var_p,
              buf1[4],buf2[4],buf3[4];
     /* Variables used for Accuracy Coefficients */

main()
{
     void     Input_X_G(),
              See_matrices(),
              Calculate(),
              See_Result(),
              Output2file(),
              See_Result1();

     Input_X_G();
     /* Reading of elements of Error Matrix */
     if(r2==0) return;

     printf("\n\n");
     printf(" Do you want see to the input matrices? (y/n) : ");
     gets(cc);
     if(cc[0] == 'y') See_matrices();
     /* Confirmation of Input Matrix */

     printf("\n\n");

     Calculate();
     /* Evaluation of Coefficients */


     printf(" Do you want to see the result of calculation? (y/n) : ");
     gets(cc);
     printf("\n\n");
     if(cc[0] == 'y') See_Result();
     /* Confirmation of Evaluation Results on Screen */

     printf(" Do you want to output the result to file?  (y/n) : ");
     gets(cc);
     printf("\n\n");
     if(cc[0] == 'y'){
          Output2file();
          printf("  Save data -----> %s \n\n",f_name[2]);
      }
     /* Output Evaluation Results to File */
 }



/*-----------------Begin of Func. Input_X_G()--------------------*/
void Input_X_G()
{

     /* Input File Name ( X and Gamma ) */
     printf("\n Length of filename should be TWO");
     printf("\n Input filename of error matrix, (e.g., X1) : ");
     gets(tmp1);
     for(i=0;i<20;i++) tmp3[i]=tmp1[i];
     f_name[0]=tmp1;

     printf("\n Enter the file name of matrix,  (e.g., Ga) : ");
     gets(tmp2);
     f_name[1]=tmp2;

     /* Reading of Error Matrix X */
     data_r1=fopen(f_name[0],"r");
     num=0;
     while(feof(data_r1)==0){
          while ( ( c = fgetc(data_r1) ) == ' ' || c == '\t' );
               ungetc(c,data_r1);
          fscanf(data_r1,"%d",&tmp[num]);
          num++;
     /* Each element is stored in tmp */
      }
     fclose(data_r1);

     r=sqrt((float)(num-1));
     /* Find Category Number */

     num=0;
     for(i=0;i<r;i++)
     for(j=0;j<r;j++){
          X[i][j]=tmp[num];
          num++;
     /* Vector tmp is rearranged into squared matrix X */
      }

     /* Input of Loss matrix Gamma */
     data_r2=fopen(f_name[1],"r");
     num=0;
     while(feof(data_r2)==0){
          while ( ( c = fgetc(data_r2) ) == ' ' || c == '\t' );
               ungetc(c,data_r2);
          fscanf(data_r2,"%d",&tmp[num]);
          num++;
     /* Each element is stored in tmp */
      }
     fclose(data_r2);

     r2=sqrt((float)(num-1));
    /* Find Category Number */

     if(r != r2){
          printf("\nWarning : Input Matrix is not squared (%s & %s)\n\n"
                  ,tmp3,tmp2);
          r2=0;
          return;
      }
     /* Warning for X and Gamma are not squared matrices */

     num=0;
     for(i=0;i<r;i++)
     for(j=0;j<r;j++){
          G1[i][j]=tmp[num];
          num++;
      }
     /* Vector tmp is rearranged into squared matrix X */
     fclose(data_r2);

     for(i=0;i<r;i++)
     for(j=0;j<r;j++){
          if(i==j) G2[i][j]=0.0;
          else G2[i][j]=1.0;
      }
     /* G2 is a rxr matrix whose diagonal elements are 0 and others are 1 */
     return;

 }
/*----------------End of Func. Input_X_G()-------------------*/


/*------------Begin of Func. See_matrices()------------------*/
void See_matrices()
{
     printf("\n %s = \n",f_name[0]);
     for(i=0;i<r;i++){
          for(j=0;j<r;j++)
          printf("%7d",X[i][j]);
          printf("\n");
      }
     /* Dipplay X on Screen */
     printf("\n\n");
     printf("  Press <return> Key ");
     gets(cc);
     printf("\n\n");

     printf(" %s = \n",f_name[1]);
     for(i=0;i<r;i++){
          for(j=0;j<r;j++)
          printf("%7.2lf",G1[i][j]);
          printf("\n");
      }
     /* Dipplay Gamma on Screen */
     printf("\n");
     return;
 }
/*-------------End of Func. See_matrices()---------------*/


/*---------------Begin of Func. Calculate()-----------------*/
void Calculate()
{
     for(i=0;i<r;i++){
          xc[i]=0;xr[i]=0;
      }
     N=0;
     sum_c=0;
     sum_a=0;

     for(i=0;i<r;i++)
     for(j=0;j<r;j++){
          N=N+X[i][j];                  /* Sample size   */
          xc[i]=xc[i]+X[i][j];          /* sum of column */
          xr[j]=xr[j]+X[i][j];          /* sum of row    */
          if(i==j) sum_a=sum_a+X[i][j]; /* sum of diagonal elements */
      }

     for(i=0;i<r;i++){
          CBX[i]=X[i][i]/(float)xc[i]; /* Class-averaged accuracy */
          sum_c=sum_c+CBX[i];          
                                       
      }

     for(i=0;i<r;i++) Pai[i]=xc[i]/(float)N;   
     /* Proportional] Prior Probability */

     AX=sum_a/(float)N;    /* Overall accuracy A(X) */
     CX=sum_c/(float)r;    /* Class-averaged accuracy C(X) */

     buf1[0]=0;buf2[0]=0;
     for(i=0;i<r;i++){
          buf1[0] += X[i][i]-(xc[i]*xr[i])/(float)N;
          buf2[0] += (xc[i]*xr[i])/(float)N;
          KapX=buf1[0]/((float)N-buf2[0]);     /* Kappa(X) */
      }

     buf1[0]=0;
     for(i=0;i<r;i++) buf1[0] += xc[i]*xr[i];
     K_hat=(sum_a-buf1[0])/(float)(N*N-buf1[0]);/* K_hat */

     for(i=0;i<r;i++) p[i] = (2*X[i][i]+1)/(float)(2*xc[i]+1);
                                                /* p_ i|i */

     /* --------- Uniform Prior Probability --------- */

     buf1[0]=1.0;
     for(i=0;i<r;i++) buf1[0] *= (2*X[i][i]+1)/(float)(2*xc[i]+1);
     J_u = pow(buf1[0],1/(float)r);
     /* J_uni(X) */

     buf1[0]=1.0;
     for(i=0;i<r;i++) buf1[0] *= X[i][i]/(float)xc[i];
     Ja_u = pow(buf1[0],1/(float)r);
     /* J*_uni(X) */

     R_u=0;
     for(i=0;i<r;i++){
          buf1[0]=(1/(float)r)/(float)xc[i];
          buf1[1]=0;
          for(j=0;j<r;j++) buf1[1] += (G1[i][j] * X[i][j]);
          R_u += (buf1[0] * buf1[1]);
      }
     /* R(X|G) with uniform priors */

     Tau_u[1]=0;
     Tau_u[2]=0;
     for(i=0;i<r;i++){
          buf1[0] = (1/(float)r)/(float)(xc[i]*xc[i]);
          buf1[1]=0;buf1[2]=0;
          buf2[1]=0;buf2[2]=0;
          for(j=0;j<r;j++){
               buf1[1] += pow(G1[i][j],2.0)*X[i][j];
               buf1[2] += G1[i][j]*X[i][j];
               buf2[1] += G2[i][j]*G2[i][j]*X[i][j];
               buf2[2] += G2[i][j]*X[i][j];
	   }
          Tau_u[1] += buf1[0]*(buf1[1]-(1/(float)xc[i])*pow(buf1[2],2.0));
          Tau_u[2] += buf1[0]*(buf2[1]-(1/(float)xc[i])*pow(buf2[2],2.0));
      }
     /* tau^2 with uniform priors */

     Var_u=0;
     for(i=0;i<r;i++)
     Var_u += (1/(float)(r*r))*((1.0-p[i])/(float)(p[i]*xc[i]));
     /* Variance sigma^2 with uniform priors */


     /* --------- Proportional Prior Probabilities --------- */

     J_p=1.0;
     for(i=0;i<r;i++)
     J_p *= pow((2*X[i][i]+1)/(float)(2*xc[i]+1),xc[i]/(float)N);
     /* J_pro(X) */

     Ja_p=1.0;
     for(i=0;i<r;i++)
     Ja_p *= pow(X[i][i]/(float)xc[i],xc[i]/(float)N);
     /* J*_pro(X) */

     R_p=0;
     for(i=0;i<r;i++){
          buf1[0] = Pai[i]/(float)xc[i];
          buf1[1]=0;
          for(j=0;j<r;j++){
               buf1[1] += (G1[i][j]* X[i][j]);
	   }
          R_p += (buf1[0] * buf1[1]);
      }
     /* R(X|G) with proportinal prios */

     Tau_p[1]=0;
      Tau_p[2]=0;
     for(i=0;i<r;i++){
          buf1[0] = Pai[i]/(float)(xc[i]*xc[i]);
          buf1[1]=0;buf1[2]=0;
          buf2[1]=0;buf2[2]=0;
          for(j=0;j<r;j++){
               buf1[1] += pow(G1[i][j],2.0)*X[i][j];
               buf1[2] += G1[i][j]*X[i][j];
               buf2[1] += pow(G2[i][j],2.0)*X[i][j];
               buf2[2] += G2[i][j]*X[i][j];
	   }
          Tau_p[1] += buf1[0]*(buf1[1]-(1/(float)xc[i])*pow(buf1[2],2.0));
          Tau_p[2] += buf1[0]*(buf2[1]-(1/(float)xc[i])*pow(buf2[2],2.0));
      }
     /* tau^2 with proportinal prios */

     Var_p=0;
     for(i=0;i<r;i++)
     Var_p += pow(Pai[i],2.0)*((1.0-p[i])/(float)(p[i]*xc[i]));
     /* Variance sigma^2 with proportinal prios */
     return;
 }
/*---------------End of Func. Calculate()-------------------*/



/*--------------Begin of Func. See_Result()-----------------*/
void See_Result()
{
     printf("   #========================================#\n");
     printf("   #                %s , %s                 #\n",tmp3,tmp2);
     printf("   #========================================#\n");
     printf("\n");
     printf("       C(X)=\t\t\t%9.4lf\n",CX);
     printf("       A(X)=\t\t\t%9.4lf\n",AX);
     printf("       Kappa(X)=\t\t%9.4lf\n",KapX);
     printf("       K_hat=\t\t\t%9.4lf\n",K_hat);
     printf("\n");
     printf("   #----------------------------------------#\n");
     printf("   #           proportional prob.           #\n");
     printf("   #----------------------------------------#\n");
     printf("       Jpro=\t\t\t%9.4lf\n",J_p);
     printf("       Jpro*=\t\t\t%9.4lf\n",Ja_p);
     printf("       R(%s|%s)=\t\t%9.4lf\n",tmp3,tmp2,R_p);
     printf("       tau^2(%s|%s)=\t\t%9.4lf\n",tmp3,tmp2,Tau_p[1]);
     printf("       Var(%s)=\t\t\t%9.4lf\n",tmp3,Var_p);
     printf("\n");
     printf("   #----------------------------------------#\n");
     printf("   #             uniform prob.              #\n");
     printf("   #----------------------------------------#\n");
     printf("       Juni=\t\t\t%9.4lf\n",J_u);
     printf("       Juni*=\t\t\t%9.4lf\n",Ja_u);
     printf("       R(%s|%s)=\t\t%9.4lf\n",tmp3,tmp2,R_u);
     printf("       tau^2(%s|%s)=\t\t%9.4lf\n",tmp3,tmp2,Tau_u[1]);
     printf("       Var(%s)=\t\t\t%9.4lf\n",tmp3,Var_u);
     printf("   #----------------------------------------#\n");
     printf("\n\n");

     return;
 }
/*--------------End of Func. See_Result()-----------------*/


/*--------------Begin of Func. Output2file()-----------------*/
void Output2file()
{
     f_name[3]=".acc";
     f_name[2]=strcat(f_name[0],f_name[1]);
     f_name[2]=strcat(f_name[2],f_name[3]);
     data_w=fopen(f_name[2],"w+");

     fprintf(data_w,"   #========================================#\n");
     fprintf(data_w,"   #                %s , %s                 #\n",
                    tmp3,tmp2);
     fprintf(data_w,"   #========================================#\n");
     fprintf(data_w,"\n");
     fprintf(data_w,"       C(X)=\t\t%18.14lf\n",CX);
     fprintf(data_w,"       A(X)=\t\t%18.14lf\n",AX);
     fprintf(data_w,"       Kappa(X)=\t% 18.14lf\n",KapX);
     fprintf(data_w,"       K_hat=\t\t%18.14lf\n",K_hat);
     fprintf(data_w,"\n");
     fprintf(data_w,"   #----------------------------------------#\n");
     fprintf(data_w,"   #           proportional prob.           #\n");
     fprintf(data_w,"   #----------------------------------------#\n");
     fprintf(data_w,"       Jpro=\t\t%18.14lf\n",J_p);
     fprintf(data_w,"       Jpro*=\t\t%18.14lf\n",Ja_p);
     fprintf(data_w,"       R(%s|%s)=\t%18.14lf\n",tmp3,tmp2,R_p);
     fprintf(data_w,"       tau^2(%s|%s)=\t%18.14lf\n",tmp3,tmp2,Tau_p[1]);
     fprintf(data_w,"       tau^2(%s|G_sp)=\t%18.14lf\n",tmp3,Tau_p[2]);

     fprintf(data_w,"       Var(%s)=\t\t%18.14lf\n",tmp3,Var_p);
     fprintf(data_w,"\n");
     fprintf(data_w,"   #----------------------------------------#\n");
     fprintf(data_w,"   #             uniform prob.              #\n");
     fprintf(data_w,"   #----------------------------------------#\n");
     fprintf(data_w,"       Juni=\t\t%18.14lf\n",J_u);
     fprintf(data_w,"       Juni*=\t\t%18.14lf\n",Ja_u);
     fprintf(data_w,"       R(%s|%s)=\t%18.14lf\n",tmp3,tmp2,R_u);
     fprintf(data_w,"       tau^2(%s|%s)=\t%18.14lf\n",tmp3,tmp2,Tau_u[1]);

     fprintf(data_w,"       tau^2(%s|G_sp)=\t%18.14lf\n",tmp3,Tau_u[2]);
     fprintf(data_w,"       Var(%s)=\t\t%18.14lf\n",tmp3,Var_u);
     fprintf(data_w,"   #----------------------------------------#\n");

     fclose(data_w);

     return;
 }
/*--------------End of Func. Output2file()-----------------*/