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

/*
float Globals::calcQE(float *v1, float *v2){
  float res = 0;
  for (int i=0;i<vectorlength;i++) {
    res += (v1[i]-v2[i]) * (v1[i]-v2[i]);
  }
  return res;
}
*/


int main(int argc, char* argv[]) {

  FILE *f,*f2;
  long nv, vd;
  float *meanVec, *tempVec;
  char dummy[1024];
  long vecs,dims,i,j;
  float elem;
  int norm = 0;
  float absVector = 0;
  float mqe=0, MQE0=0;

  if (argc<3) {
    printf("Usage: calc_mqe0 inputvectorfile outputfile [-n]\n");
    exit(1);
  }

  if ((argc==4) && (strcmp(argv[3],"-n")==0)) {
    norm = 1;
    printf("norm vectors on\n");
  }

  if ((f = fopen(argv[1],"r")) == NULL) {
    printf("cannot open input vector file %s\n",argv[1]);
    exit(1);
  }

  /* TYPE */
  fscanf(f,"%s",&dummy);
  fscanf(f,"%s",&dummy);
  /* XDIM */
  fscanf(f,"%s",&dummy);
  if (strcmp(dummy,"$XDIM")!=0) {
    fclose(f);
    printf("cannot read $XDIM\n");
    exit(1);
  }
  fscanf(f,"%d",&nv);
  /* YDIM */
  fscanf(f,"%s",&dummy);
  fscanf(f,"%s",&dummy);
  /* VEC_DIM */
  fscanf(f,"%s",&dummy);
  if (strcmp(dummy,"$VEC_DIM")!=0) {
    fclose(f);
    printf("cannot read $VEC_DIM\n");
    exit(1);
  }
  fscanf(f,"%d",&vd);

  meanVec = (float *) calloc(vd,sizeof(float));
  for (i=0;i<vd;i++) {
    meanVec[i] = 0;
  }
  tempVec = (float *) calloc(vd,sizeof(float));
  for (i=0;i<vd;i++) {
    tempVec[i] = 0;
  }

  if (norm == 0) { // no normalization
    for (vecs=0;vecs<nv;vecs++) {
      for (dims=0;dims<vd;dims++) {
	fscanf(f,"%f",&elem);
	meanVec[dims] += elem;
      }
      fscanf(f,"%s",&dummy);
      printf("added %s\n",dummy);
    }
    
    for (i=0;i<vd;i++) {
      meanVec[i] = meanVec[i] / vd;
      //    printf("%f ",meanVec[i]);
    }
  } else { // norm vectors
    for (vecs=0;vecs<nv;vecs++) {
      for (dims=0;dims<vd;dims++) {
	fscanf(f,"%f",&elem);
	tempVec[dims] = elem;
      }
      fscanf(f,"%s",&dummy);
      printf("added %s\n",dummy);
      // *******************************
      absVector = 0;
      for (j=0;j<vd;j++) {
	absVector += tempVec[j] * tempVec[j];
      }
      absVector = sqrt(absVector);

      if (absVector>0) {
	for (j=0;j<vd;j++) {
	  tempVec[j] /= absVector;
	  meanVec[j] += tempVec[j];
	}
      } else {
	for (j=0;j<vd;j++) {
	  meanVec[j] += tempVec[j];
	}
      }
      // *******************************
    }
    for (j=0;j<vd;j++) {
      meanVec[j] = meanVec[j] / vd;
    }
  }

  fclose(f);

  // ************* MQE0 **************
  if ((f = fopen(argv[1],"r")) == NULL) {
    printf("cannot open input vector file %s\n",argv[1]);
    exit(1);
  }

  for (j=0;j<8;j++) {
    fscanf(f,"%s",&dummy);
  }
  
  if (norm == 0) { // no normalization
    MQE0 = 0;
    for (vecs=0;vecs<nv;vecs++) {
      mqe = 0;
      for (dims=0;dims<vd;dims++) {
	fscanf(f,"%f",&elem);
	tempVec[dims] = elem;
      }
      fscanf(f,"%s",&dummy);
      
      for (i=0;i<vd;i++) {
	mqe += (tempVec[i]-meanVec[i]) * (tempVec[i]-meanVec[i]);
      }
      MQE0 += mqe;
      printf("calculated %s\n",dummy);
    }
    //    MQE0 = MQE0 / (float) nv;
    
  } else { // norm vectors
    MQE0 = 0;
    for (vecs=0;vecs<nv;vecs++) {
      for (dims=0;dims<vd;dims++) {
	fscanf(f,"%f",&elem);
	tempVec[dims] = elem;
      }
      fscanf(f,"%s",&dummy);
      //      printf("added %s\n",dummy);
      // *******************************
      absVector = 0;
      for (j=0;j<vd;j++) {
	absVector += tempVec[j] * tempVec[j];
      }
      absVector = sqrt(absVector);
      
      if (absVector>0) {
	for (j=0;j<vd;j++) {
	  tempVec[j] /= absVector;
	}
      }

      mqe=0;
      for (i=0;i<vd;i++) {
	mqe += (tempVec[i]-meanVec[i]) * (tempVec[i]-meanVec[i]);
      }
      MQE0 += mqe;
      printf("calculated %s\n",dummy);
    
    }
    //    MQE0 = MQE0 / (float) nv;

    // *******************************

  }
  

  fclose(f);

  // ************* write file *************
  if ((f2 = fopen(argv[2],"w")) == NULL) {
    printf("cannot open ouput file %s\n",argv[2]);
    exit(1);
  }

  fprintf(f2,"$MEAN_VEC ");
  for (i=0;i<vd-1;i++) {
    fprintf(f2,"%f ",meanVec[i]);
  }
  fprintf(f2,"%f\n",meanVec[vd-1]);

  fprintf(f2,"$MQE0 %f\n",MQE0);

  fclose(f2);

  return 0;
}