#include<stdio.h>
#include<stdlib.h>
#include<math.h>

/*
이 프로그램은 전산통계 과목에 사용된 프로그램으로
주어진 행렬의 HalfMatrix를 구 하고
HalfMatrix의 역행렬인 B행렬을 구하는 프로그램입니다.

목적 : 전치행렬 및 행렬의 곱을 구하는 프로그램.

*/

int indx(int i, int j);

void main() {
 int
  n, p, //det, /* for size of matrix */
  ij, i, j, k, fp,/* index */
  df, dfn, dfp;
 double
  *x, *z, *h, //*b, /* for matrix */
  sum,  //각각의 합과 역행렬의 determinant
  SSR=0, SSE=0, SST=0,
  MSE=0, MSR=0, RMSE=0,
  RR=0, F=0, DPM=0;
 printf("자료인 X행렬을 통해 Z행렬을 찾아내고\\nZ행렬의 Half Maxrix를 이용해 ANOVA Table을 구하는 프로그램입니다.\\n\\n"); 
 printf("먼저, X행렬을 출력하겠습니다.\\n\\n  =X Matrix=\\n");

 FILE *fp;
 /* 입력 파일 정하기 */
 //현재 소스 폴더에 있는 파일을 읽을때 그냥 " "안에만 넣으면 읽히는것으로 알고 있습니다.

 if((fp=fopen("C:\\\\reg.txt","r"))==NULL){
  printf("Error Opening file \\n");
  exit(1);
 }

 /* 행의 개수와 열의 개수 읽기 */
 fscanf(fp,"%d %d",&n,&p);

 /* 메모리 할당 */
 x=(double *) malloc(n*(p+1)*sizeof(double));

 /* 자료 읽기 */
 for (i=0; i<n; i++) {    
  for (j=0; j<(p+1); j++) {//X행렬의 가장 앞에 1이 오므로 열이 하나 늘어나게 되어 p+1까지 계산합니다.
   if(j==0) x[i*(p+1)+j]=1; //X행렬의 가장 앞은 1이 저장되도록 했습니다.
   else{
   fscanf(fp,"%lf",&x[i*(p+1)+j]);
  }}
 }

 //Matrix X의 출력
 for (i=0; i<n; i++) {
  printf("│");
  for (j=0; j<(p+1); j++) {
   printf("%7.2f ",x[(p+1)*i+j]);
  }
  printf("│\\n");
 }
 printf("\\n다음, X\'X행렬인 Z를 출력하겠습니다.\\n\\n  =Z Matrix=\\n");
 //X\'X인 Z행렬의 저장 및 출력.

  z=(double *) malloc((p+1)*(p+2)/2*sizeof(double)); //행렬의 크기가 늘어나므로 동적할당에도 용량을 늘렸습니다.
 for(i=0; i<p+1; i++){
  for(j=0; j<=i;j++){
   ij=indx(i,j);//상삼각행렬이 되도록 ij값을 indx를 이용하여 적절하게 바꿔줍니다.
   z[ij]=0; //행렬의 초기화
   for (k=0; k<n; k++) z[ij] += x[k*(p+1)+i]*x[k*(p+1)+j]; //X\'X 의 계산
  }
 }
 //z출력
 for (i=0; i<p+1; i++) {
  printf("│");
  for (j=0; j<p+1; j++) {
   printf("%9.2f ",z[indx(i,j)]);
  }
  printf("│");
  printf("\\n");
 }

        
 //Half Metrix//
 
 printf("\\n여기서 Z의 Half Matrix인 H행렬을 구해 보겠습니다.\\n");
  printf("  =H Matrix=\\n");
  h=(double *) malloc((p+1)*(p+2)/2*sizeof(double));
 for (i=0; i<n; i++){
  for(j=0; j<p+1; j++){
   if (i==j){  //i와 j가 같을 경우입니다.
     sum=0;  
     //sigma를 계산할때 사용하기 위한 변수입니다. 실행될때마다 초기화합니다.
     for(k=0; k<i; k++){ // sigma 계산
       sum += h[(indx(k,i))]*h[(indx(k,i))];  
       //uppertriangle 이므로 indx를 이용하여 불러옵니다
     }     
    h[(indx(i,j))]=sqrt(z[indx(i,j)]-sum);
    //uppertriangle이므로 h의 저장을 이렇게 합니다.
    //a메트릭스는 triangle이 아니므로 일반적인 방법으로 불러옵니다.  
   }
   else if(i<j) {  //j가 i보다 큰 경우입니다.
     sum=0;
     //sigma를 계산할때 사용하기 위한 변수입니다. 실행될때마다 초기화합니다.
     for(k=0; k<i; k++){      
       sum += h[(indx(k,i))]*h[(indx(k,j))];
     }
     //h(ii)값이 0에 가까우면 오류 메시지 출력
     if(h[(indx(i,i))] > -0.00001 && h[(indx(i,i))] < 0.00001){
      printf("\\nh(ii)값이 너무 작아 진행할 수 없습니다.\\n프로그램을 종료합니다.\\n\\n");
      exit(1);
     }
     else h[(indx(i,j))]=(z[indx(i,j)]-sum)/h[(indx(i,i))];
   }
   //이 부분에 i가 j보다 큰 경우를 적어주어야 하는데 h메트릭스의 저장공간 절약을 위해 0을 저장하지 않으므로
   //이 경우에는 그냥 건너 뛰도록 아무문장도 넣지 않았습니다. 즉 j>=i일때만 h에 저장하도록 했습니다.
  }
 }
 //메트릭스의 출력.
 for (i=0; i<p+1; i++) {
  printf("│");
  for (j=0; j<p+1; j++) {
    if(j>=i) printf("%9.4f ",h[(indx(i,j))]);//j가i보다 같거나 더 클 경우만 출력되도록 했습니다.
               //이렇게 하지 않으면 symmatric한 모습으로 출력되겠지요.
    else printf("%9.4f ",0);     //i가 더 클경우에는 0이 출력되도록 했습니다.
  }
  printf("│\\n");
 } 
/*  B메트릭스를 구하는 소스는 숨겨놓았습니다.

  printf("\\n그럼 H행렬의 역행렬인 B행렬을 구해 보겠습니다.\\n\\n");
  printf("  =B Matrix=\\n");
//B Matrix
  b=(double *) malloc((p+1)*(p+2)/2*sizeof(double));
 det=1;
 for(i=0;i<p;i++){
  det = h[indx(i,i)]*det;
 }
 if(det > -0.00001 && det < 0.00001){
  printf("\\n행렬식의 크기가 너무작아 진행할 수 없습니다.\\n프로그램을 종료합니다.\\n\\n");
  exit(1);
 }//행렬식의 크기가 0에 가까우면 실행되지 않도록 했습니다.
 for(i=0; i<p+1; i++){
  for(j=0; j<p+1; j++){
   if(i==j){ //i와 j가 같은 경우입니다.
    b[(indx(i,j))] = 1/h[(indx(i,j))];
    //b와 h행렬이 모두 Uppertriangle이므로 indx문을 사용해서 입력
   }
   else if(i<j){  //j가 i보다 큰 경우입니다.
    sum=0;
     for(k=i; k<j; k++){  //sigma의 계산
      sum += b[(indx(i,k))]*h[(indx(k,j))];
     }
    if(h[(indx(j,j))] > -0.00001 && h[(indx(j,j))] < 0.00001){
     printf("\\nh(jj)값이 너무 작아 진행할 수 없습니다.\\n프로그램을 종료합니다.\\n\\n");
     exit(1);
    }//h(jj)가 0의 값에 가까우면 오류가 나도록 설정하였습니다.
    else b[(indx(i,j))] = -sum/h[(indx(j,j))];
    //b행렬 구하는 식
   }
   //이곳에 i가j보다 큰 경우를 적어주어야 하는데 b메트릭스의 저장공간의 절약을 위해 0을 저장하지 않으므로
   //i가 j보다 큰 경우엔 그냥 건너뛰도록 아무 문장도 넣지 않았습니다.
  }
 }

 for (i=0; i<p+1; i++) {
  printf("│");
  for (j=0; j<p+1; j++) {
    if(j>=i) printf("%9.4f ",b[(indx(i,j))]);
    //j가i보다 같거나 더 클 경우만 출력되도록 했습니다
    else printf("%9.4f ",0);
  }
  printf("│\\n");
 }
*/

 for(i=1 ; i<p ; i++){ //SSR입니다. 제일 위의 행과 아래 행을 제외하고 제곱하여 구합니다.
  SSR += h[indx(i,p)]*h[indx(i,p)];
 }
 for(k=0 ; k<n ; k++){ //Dependent Mean 을 구하기 위해서 사용한 반복문입니다.
  DPM += x[k*(p+1)+p];
 }
 dfp = p-1;
 dfn = n-1;
 df = dfn-dfp; // 각각의 자유도들 입니다.
 SSE = h[indx(p,p)]*h[indx(p,p)]; //SSE는 h행렬의 마지막줄 마지막열의 곱으로 나타납니다.
 SST = SSR+SSE; // SST는 SSR+SSE로 나타 낼 수 있습니다.
 RR = SSR/SST; //R-Square입니다.
 MSE = SSE/df; //MSE는 SSE를 오차의 자유도로 나눈값입니다.
 MSR = SSR/dfp; //MSR은 SSR을 변수의 자유도로 나눈 값입니다.
 F =  MSR/MSE; //F값은 MSR/MSE로 나타냅니다.
 RMSE = sqrt(MSE); //Root MSE입니다. MSE값에 제곱근을 한 값입니다.
 DPM = DPM/n; //위의 DPM값을 n으로 나누어서 구한 y값의 평균입니다.
 printf("\\nANOVA TABLE을 출력하겠습니다.");
 
 //아노바 테이블의 출력입니다.
 
 printf("\\n\\n              Analysis of Variance(ANOVA TABLE)\\n\\n");
 printf("┌────────┬─┬───────┬──────┬────┐\\n");
 printf("│Source          │DF│Sum of Squares│ Mean Square│ F-Value│\\n");
 printf("├────────┼─┼───────┼──────┼────┤\\n");
 printf("│Model           │%2d│%14.0f│%12.5f│%8.2f│\\n",dfp,SSR,MSR,F);
 printf("│Error           │%2d│%14.5f│%12.5f│        │\\n",df,SSE,MSE);
 printf("│Corrected Total │%2d│%14.0f│            │        │\\n",dfn,SST);
 printf("├────────┼─┴───┬───┴──────┴────┘\\n");
 printf("│Root MSE        │%10.5f│\\n│R-Square        │%10.4f│\\n",RMSE,RR);
 printf("│Dependent Mean  │%10.5f│\\n",DPM);
 printf("└────────┴─────┘\\n");
 printf("\\n이상으로 프로그램을 종료합니다.\\n");

}

int indx(int i, int j) //upper triangle로 저장되도록 약간 변경 하였습니다.
{
 if (j>=i) return (j*(j+1)/2+i);
 else return (i*(i+1)/2+j);
}