/*===================================================
This code evaluates the AVERAGE RECIPROCAL DISTANCE on CPU
Author: Radu Herbei 
=====================================================*/

#include <iostream>
#include <iomanip>
#include <fstream>
#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <math.h>

#include "MersenneTwister.h"
using namespace std;


#define N 100 // 
#define P 10 // 


// Matrices are stored in row-major order
typedef struct {
    int width;
    int height;
    double* elements;
} Matrix;


// Matrices are stored in row-major order
typedef struct {
    int width;
    int height;
    int* elements;
} intMatrix;



int mypower(int a, int b)
{
     int c=a;
     for (int n=b; n>1; n--) c*=a;
     return c;
}





int main(void){


  MTRand mt;
  Matrix X; // this is the full design matrix
  int i,j,col;
 
  X.width=P;
  X.height=N;
  X.elements = (double*) malloc(X.width*X.height*sizeof(double));

  double TT;
  TT=(N*(N-1)/2)*(mypower(2,P)-1);

  FILE *file;
  float value;
  file = fopen("design100.txt", "r");
  for(i=0;i<N*P;i++){
    fscanf(file, "%f", &value);
    X.elements[i]=value/99;
  }
  fclose(file);




  clock_t startcpu=clock();


  //create the N choose 2 matrix
  intMatrix NC2;
  NC2.width=N*(N-1)/2;
  NC2.height=3;
  NC2.elements=(int *) malloc(NC2.height*NC2.width*sizeof(int));

  col=0;
  for(i=0;i<N;i++)
    for(j=0;j<N;j++){
      if(i>j){
	NC2.elements[0*NC2.width+col]=i;
	NC2.elements[1*NC2.width+col]=j;
	col++;
      }
    }
	

  intMatrix JK;
  JK.height=P;
  JK.width=mypower(2,P)-1;
  JK.elements=(int*)malloc(JK.height*JK.width*sizeof(int));
  // fill JK
  int temp,qt,rem;
  for(j=1;j<=JK.width;j++){
    for(i=0;i<P;i++) JK.elements[i*JK.width+(j-1)]=0;
    
    temp=j;
    i=0;
    while(temp>0){
      qt = temp/2; 
      rem= temp%2;
      if(rem==0){
	JK.elements[i*JK.width+(j-1)]=0;
      }
      else{
	JK.elements[i*JK.width+(j-1)]=1;
      }
      i++;
      temp=qt;
    }
  }




  double ard=0.0;
  int row1, row2,k,jj;
  double sum=0.0;


  for(i=0;i<(mypower(2,P)-1);i++)
    for(j=0;j<(N*(N-1)/2);j++){

      // look up column j in NC2 matrix;
      row1=NC2.elements[0*NC2.width+j];
      row2=NC2.elements[1*NC2.width+j];


      //look up column i in JK matrix;
      sum=0.0;
      jj=0;
      for(k=0;k<P;k++){
	if ( JK.elements[k*JK.width+i]>0 ){
	  jj++;
	  sum += (X.elements[row1*X.width+k]-X.elements[row2*X.width+k])*(X.elements[row1*X.width+k]-X.elements[row2*X.width+k]);
	}
      }

      ard += sqrt(jj/sum);
	
    }
      
  printf("ARD distance is %f \n", ard/TT);
  printf("Time to compute on CPU : %f seconds.\n ", ((double)clock() - startcpu)/CLOCKS_PER_SEC );
  
  free(X.elements);
  free(NC2.elements);
  free(JK.elements);



}