//#include "gpu.hpp"
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
#include <CL/sycl.hpp>
#include <CL/sycl/INTEL/fpga_extensions.hpp>

typedef long long   int64;
typedef unsigned long long uint64;
#define FMT_INT64   "%lld"
#define FMT_UINT64   "%llu"
#define FMT_HEX64   "%llx"

struct genotype {
	uint64 *genocase;
	uint64 *genoctrl;
};

struct MarginalDistr {
	int MarginalDistrSNP[3];
	int MarginalDistrSNP_Y[3][2];
};
using namespace std;
using namespace cl::sycl;

double Abs(double a)
{
	return((a < 0) ? -a : a);
}
int main() {
	//void calulateJointResult(int n, int p, int ncase, int nctrl, int **GenoDistr) {

	queue Q{ cl::sycl::gpu_selector{} };
	//queue Q{ INTEL::fpga_emulator_selector{} };
	std::cout << "Running on GPU device: " << Q.get_device().get_info<info::device::name>() << "\n";
	time_t st, ed;
	time(&st);

	int  nlongintcase, nlongintctrl;
	int flag, i, ii, j, k, j1, j2;
	int c, tmp;
	int n = 2000, p = 1000, ncase = 987, nctrl = 1013;
	double thresholdRecord = 30;

	double *InteractionMeasureG = malloc_host<double>(p*(p - 1) / 2, Q);
	int InteractionCount = 0;

	FILE *fpRJiont = fopen("JiontB.txt", "r");
	int **GenoDistr_D = malloc_host<int*>(p*(p - 1) / 2, Q);    //连列表 3*3*2
	for (int t = 0; t < p*(p - 1) / 2; t++)GenoDistr_D[t] = malloc_host<int>(2 * 9, Q);    //连列表 3*3*2
	for (int m = 0; m < p * (p - 1) / 2; m++) {
		for (int t = 0; t < 18; t++)
			fscanf(fpRJiont, "%d", &GenoDistr_D[m][t]);
	}
	fclose(fpRJiont);
	printf(" load GenoDistr is over\n");

	Q.submit([&](auto &h) {
		h.parallel_for(range{ size_t(p - 1), size_t(p) }, [=](id<2> idx) {
			int j1 = idx[0]; int j2 = idx[1];
			if (j2 > j1) {
				int sum, i, j;
				double ptmp1, ptmp2;
				if (j1 == 0)sum = 0;
				else if (j1 % 2 == 0)sum = ((j1 - 1) / 2)*(j1 + 1) + 1;
				else sum = (j1 - 1) / 2 * j1;
				int loc = j1 * (p - 1) - sum + j2 - 1 - j1;
				double Pab[9], Pbc[6], Pca[6];   //conditional probability P(a|b), P(b|c), P(c|a)
				//计算pMarginalDistr
				int MarginalDistrJ1[3], MarginalDistrJ2[3];
				int MarginalDistrJ1_Y[3][2], MarginalDistrJ2_Y[3][2];
				for (i = 0; i < 3; i++) {
					for (int j = 0; j < 2; j++) {
						MarginalDistrJ1_Y[i][j] = GenoDistr_D[loc][j * 9 + i * 3] + GenoDistr_D[loc][j * 9 + i * 3 + 1] + GenoDistr_D[loc][j * 9 + i * 3 + 2];
						MarginalDistrJ2_Y[i][j] = GenoDistr_D[loc][j * 9 + i] + GenoDistr_D[loc][j * 9 + 1 * 3 + i] + GenoDistr_D[loc][j * 9 + 2 * 3 + i];
						Pbc[3 * j + i] = (double)MarginalDistrJ2_Y[i][j];
					}
					MarginalDistrJ1[i] = MarginalDistrJ1_Y[i][0] + MarginalDistrJ1_Y[i][1];
					MarginalDistrJ2[i] = MarginalDistrJ2_Y[i][0] + MarginalDistrJ2_Y[i][1];

					j = 0;
					Pca[2 * i + j] = (double)MarginalDistrJ1_Y[i][j] / MarginalDistrJ1[i];
					Pbc[3 * j + i] = Pbc[3 * j + i] / ncase;
					Pab[3 * i + j] = (double)(GenoDistr_D[loc][3 * j + i] + GenoDistr_D[loc][3 * j + i + 9]) / MarginalDistrJ2[i];

					j = 1;
					Pca[2 * i + j] = (double)MarginalDistrJ1_Y[i][j] / MarginalDistrJ1[i];
					Pbc[3 * j + i] = Pbc[3 * j + i] / nctrl;
					Pab[3 * i + j] = (double)(GenoDistr_D[loc][3 * j + i] + GenoDistr_D[loc][3 * j + i + 9]) / MarginalDistrJ2[i];

					j = 2;
					Pab[3 * i + j] = (double)(GenoDistr_D[loc][3 * j + i] + GenoDistr_D[loc][3 * j + i + 9]) / MarginalDistrJ2[i];
				}
				double tao = 0.0;
				double InteractionMeasure = 0.0;
				//double ptmp1, ptmp2;
				for (i = 0; i < 3; i++)// index for A
				{
					for (j = 0; j < 3; j++) //index for B
					{
						for (int k = 0; k < 2; k++)	//index for C
						{
							ptmp1 = (double)GenoDistr_D[loc][k * 9 + i * 3 + j] / n;
							if (ptmp1 > 0)
							{
								InteractionMeasure += ptmp1 * log(ptmp1);
							}

							ptmp2 = Pab[3 * j + i] * Pbc[3 * k + j] * Pca[2 * i + k];
							if (ptmp2 > 0)
							{
								InteractionMeasure += -ptmp1 * log(ptmp2);
								tao += ptmp2;
							}
						}
					}
				}

				InteractionMeasure = (InteractionMeasure + log(tao))*n * 2;
				InteractionMeasureG[loc] = InteractionMeasure;
			}//if(j2>j1)
	//}
		});
	});
	InteractionCount = 0;
	for (int j1 = 0; j1 < p*(p - 1) / 2; j1++) {
		if (InteractionMeasureG[j1] > thresholdRecord) {
			InteractionCount++;
		}
	}
	printf("the InteractionCount = %d\n", InteractionCount);


	time(&ed);
	printf("B2 cputime for run: %d seconds\n", (int)ed - st);
	free(InteractionMeasureG, Q);
	for (int t = 0; t < p*(p - 1) / 2; t++)free(GenoDistr_D[t], Q);    //连列表 3*3*2
	free(GenoDistr_D, Q);

	printf("RUN on GPU is over\n");
	return 0;
}