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

#include "MersenneTwister.h"
#include "book.h"
using namespace std;


#define N 100 // 
#define P 10 // 
#define Z 2 //
#define TPB 100 // threads per block


// 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;

__device__ intMatrix devNC2;
__device__ intMatrix devJK;

// row-major storing
#define IDX2(i,j) (((i)*(P))+(j))




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




__global__ void GPU_ARD(Matrix X, Matrix ARD, intMatrix NC2, intMatrix JK){

  int bidx=blockIdx.x; // 0 .. 
  int bidy=blockIdx.y; // 0.. 2^P-1
  int tidx=threadIdx.x;

  int thread = bidx * TPB + tidx ;

  if (thread < (N*(N-1)/2)){
      // look up thread col in NC2 matrix
      int row1, row2;
      
      row1=NC2.elements[1*NC2.width+thread];
      row2=NC2.elements[2*NC2.width+thread];
      
      
      // look up bidy col in JK matrix
      double xx=0.0;
      int i;
      int j=0;
      for(i=0;i<P;i++){
	if(JK.elements[i*JK.width+bidy]>0){
	  j++;
	  xx += (X.elements[row1*X.width + i] - X.elements[row2*X.width+i])*(X.elements[row1*X.width + i] - X.elements[row2*X.width+i]);
	}
      }
      
      ARD.elements[bidy*ARD.width+thread] = sqrt(j/xx);
    }
}



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]=col;
	NC2.elements[1*NC2.width+col]=i;
	NC2.elements[2*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 aa,qt,rem;
  for(j=1;j<=JK.width;j++){
    for(i=0;i<P;i++) JK.elements[i*JK.width+(j-1)]=0;
    
    aa=j;
    i=0;
    while(aa>0){
      qt = aa/2; 
      rem= aa%2;
      if(rem==0){
	JK.elements[i*JK.width+(j-1)]=0;
      }
      else{
	JK.elements[i*JK.width+(j-1)]=1;
      }
      i++;
      aa=qt;
    }
  }
  
  
  
  
  double ard1=0.0;
  int roww1, roww2,k,jj;
  double xxx=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;
      roww1=NC2.elements[1*NC2.width+j];
      roww2=NC2.elements[2*NC2.width+j];
      
      
      //look up column i in JK matrix;
      xxx=0.0;
      jj=0;
      for(k=0;k<P;k++){
	if ( JK.elements[k*JK.width+i]>0 ){
	  jj++;
	  xxx += (X.elements[roww1*X.width+k]-X.elements[roww2*X.width+k])*(X.elements[roww1*X.width+k]-X.elements[roww2*X.width+k]);
	}
      }
      
      ard1 += sqrt(jj/xxx);
      
    }
  
  printf("Time to compute on CPU : %f seconds.\n ", ((double)clock() - startcpu)/CLOCKS_PER_SEC );
  
  printf("ARD distance is %f \n", ard1/TT);
  
  

  Matrix devX;
  devX.width=X.width;
  devX.height=X.height;
  HANDLE_ERROR(cudaMalloc((void **)&devX.elements, devX.width*devX.height*sizeof(double)));
  
  intMatrix devNC2;
  intMatrix devJK;

  devNC2.width=NC2.width;
  devNC2.height=NC2.height;
  HANDLE_ERROR(cudaMalloc((void **)&devNC2.elements, devNC2.width*devNC2.height*sizeof(int)));
  
  devJK.width=JK.width;
  devJK.height=JK.height;
  HANDLE_ERROR(cudaMalloc((void **)&devJK.elements, devJK.width*devJK.height*sizeof(int)));
  
  Matrix devARD;
  devARD.width=N*(N-1)/2;
  devARD.height=mypower(2,P)-1;
  HANDLE_ERROR(cudaMalloc((void **)&devARD.elements, devARD.width*devARD.height*sizeof(double)));
  

  Matrix ARD;
  ARD.width=N*(N-1)/2;
  ARD.height=mypower(2,P)-1;
  ARD.elements=(double*) malloc(ARD.width*ARD.height*sizeof(double));

  clock_t startgpu=clock();


  HANDLE_ERROR(cudaMemcpy(devX.elements, X.elements, X.width*X.height*sizeof(double), cudaMemcpyHostToDevice)) ;
  HANDLE_ERROR(cudaMemcpy(devNC2.elements, NC2.elements, NC2.width*NC2.height*sizeof(int), cudaMemcpyHostToDevice)) ;
  HANDLE_ERROR(cudaMemcpy(devJK.elements, JK.elements, JK.width*JK.height*sizeof(int), cudaMemcpyHostToDevice)) ;
 

  
  int xdim, ydim;
  xdim = (N*(N-1)/2)/TPB + 1;
  xdim=N*(N-1)/2;
  ydim=mypower(2,P)-1;

  dim3 grid(xdim,ydim);
  GPU_ARD<<<grid,TPB>>>(devX, devARD, devNC2, devJK);
  

  HANDLE_ERROR(cudaMemcpy(ARD.elements, devARD.elements, devARD.width*devARD.height*sizeof(double), cudaMemcpyDeviceToHost)) ;
  
  double ard=0.0;
  for(i=0;i<ARD.height;i++)
    for(j=0;j<ARD.width;j++)
      ard += ARD.elements[i*ARD.width+j];
  
  printf("Time to compute on GPU : %f seconds.\n ", ((double)clock() - startgpu)/CLOCKS_PER_SEC );
  
  printf("ARD distance is %f \n", ard/TT);
  

}