#include <stdio.h> //printf
#include <stdlib.h>
#include <stdint.h>
#include <iostream> //cout
#include "opencv2/opencv.hpp" //cv::
#include "opencv2/gpu/gpu.hpp" //gpu::
#include "opencv2/gpu/gpumat.hpp" //gpu::GpuMat
#include <fstream> //streampos
#include <ctime> //time
#include <math.h> //exp,log,sqrt,pow
using namespace std;
using namespace cv;
using namespace gpu;

#define time double(std::clock())/CLOCKS_PER_SEC

double simpleKhatBeta(double N,double* distance2,double alpha2, double beta, double* nohhldij,double nohhldij_diagsum){
    double simplekhat;
    double sum = 0;
    for (int i = 0; i < 9250572; i++)
    {
        sum += (nohhldij[i] * alpha2)/(distance2[i] + alpha2);
    }
    simplekhat = N * beta/(sum + nohhldij_diagsum);
    return(simplekhat);
}

double calcLambdaSimpleKhat(double alpha2,double beta,double DS[21743733],double DKAOarray[977][61],double N,double* distance2,double* nohhldij,double nohhldij_diagsum){
    double S_LL = 0, KAO_LL = 0;
    double simplekhatbeta = simpleKhatBeta(N,distance2,alpha2,beta,nohhldij,nohhldij_diagsum);
    for (int i = 0; i < 21743733; i++)
    {
        S_LL += alpha2/(DS[i] + alpha2);
    }
    S_LL = simplekhatbeta * S_LL;
    double rowtotal[977] = {0};
    for (int i = 0; i < 977; i++)
    {
        for (int j = 0; j < 60; j++)
        {
            rowtotal[i] += (alpha2 / (DKAOarray[i][j] + alpha2));
        }
        KAO_LL += log(1 - exp(simplekhatbeta * (rowtotal[i] + DKAOarray[i][60]))); //take off NA contributions that are represented as zeroes within DKAOarray
    }
    //cout << "S_LL = " << S_LL << " KAO_LL = " << KAO_LL << "\n";
    return(S_LL+KAO_LL);
}

int main(){
    ifstream fileDS("DSvector21743733.bin", ios::in|ios::binary|ios::ate);
    
    // allocate enough space for DS due to its large array
    double* DS;
    long idx=21743733;
    DS = (double*)malloc((idx+1) * sizeof(DS[0]));
    if(DS == NULL)
        {
           printf("not enough memory");
           return 1;
        } 

    if (fileDS.is_open())
    {
        streampos size = fileDS.tellg();
        cout << "size=" << size << " bytes\n";
        char * memblock = new char [size];
        fileDS.seekg (0, ios::beg);
        fileDS.read (memblock, size);
        fileDS.close();
        cout << "DS is in memory \n";
        DS = (double*)memblock; //reinterpret as doubles
    }
    else 
    {
        cout << "Unable to open file";
    }
    printf("DS[0] = %f\n", DS[0]);

    ifstream filedistance2("distance2_9250572.bin", ios::in|ios::binary|ios::ate);
    double* distance2;
    if (filedistance2.is_open())
    {
        streampos size = filedistance2.tellg();
        cout << "size=" << size << " bytes\n";
        char * memblock = new char [size];
        filedistance2.seekg (0, ios::beg);
        filedistance2.read (memblock, size);
        filedistance2.close();
        cout << "distance2 is in memory \n";
        distance2 = (double*)memblock; //reinterpret as doubles
    }
    else 
    {
        cout << "Unable to open file";
    }
    printf("distance2[0] = %f\n", distance2[0]);

    ifstream filenohhldij("nohhldij_9250572.bin", ios::in|ios::binary|ios::ate);
    double* nohhldij;
    if (filenohhldij.is_open())
    {
        streampos size = filenohhldij.tellg();
        cout << "size=" << size << " bytes\n";
        char * memblock = new char [size];
        filenohhldij.seekg (0, ios::beg);
        filenohhldij.read (memblock, size);
        filenohhldij.close();
        cout << "nohhldij is in memory \n";
        nohhldij = (double*)memblock; //reinterpret as doubles
    }
    else 
    {
        cout << "Unable to open file";
    }
    printf("nohhldij[0] = %f\n", nohhldij[0]);

    ifstream fileDKAO("DKAOvector59597.bin", ios::in|ios::binary|ios::ate);
    double* DKAO;
    if (fileDKAO.is_open())
    {
        streampos size = fileDKAO.tellg();
        cout << "size=" << size << " bytes\n";
        char * memblock = new char [size];
        fileDKAO.seekg (0, ios::beg);
        fileDKAO.read (memblock, size);
        fileDKAO.close();
        cout << "DKAO is in memory \n";
        DKAO = (double*)memblock; //reinterpret as floats
    }
    else 
    {
        cout << "Unable to open file";
    }
    printf("DKAO[0] = %f\n", DKAO[0]);

        // read DKAO(a vector) into an array
    
    double DKAOarray[977][61];
    int vectorcount = 0;
    for (int i = 0; i < 977; i++)
    {
        for (int j = 0; j < 61; j++)
        {
            DKAOarray[i][j] = DKAO[vectorcount];
            vectorcount = vectorcount + 1;
        }
    }

    double DKAOarray60[977][60];
    double DKAOarray1[977][1];
    for (int i = 0; i < 977; i++)
    {
        for (int j = 0; j < 60; j++)
        {
            DKAOarray60[i][j] = DKAOarray[i][j];
        }
        for (int j = 60; j < 61; j++)
        {
            DKAOarray1[i][0] = DKAOarray[i][60];
        }
    }

    double sum = 0;
    for (int i = 0; i < 21743733; i++) sum += DS[i];
    printf("\nCPU SUMS\nDSsum = %f\n",sum);
    sum = 0;
    for (int i = 0; i < 9250572; i++) sum += distance2[i];
    printf("distance2sum = %f\n",sum);
    sum = 0;
    for (int i = 0; i < 9250572; i++) sum += nohhldij[i];
    printf("nohhldsum = %f\n",sum);
    sum = 0;
    for (int i = 0; i < 59597; i++) sum += DKAO[i];
    printf("DKAOsum = %f\n\n",sum);
    
    // init values
    double logL;
    double alpha = 1.2;
    double alpha2 = alpha*alpha;
    double beta = 0.6;
    beta = -1* beta; //taking negative beta to save computation for later
    double minussimplekhatbeta;

    //eval logL & print output
    double N = 25505;
    double nohhldij_diagsum = 126712;
    double cputime;
    double simplekhatcputime;
    double calclambdacputime;
    double sumcputime = 0;
    double sumsimplekhatcputime = 0;
    double sumcalclambdacputime = 0;
    double sum2cputime = 0;
    double sum2simplekhatcputime = 0;
    double sum2calclambdacputime = 0;
    for (int i = 0; i < 100; ++i)
    {
        cputime = time;
        simplekhatcputime = time;
        minussimplekhatbeta = simpleKhatBeta(N, distance2, alpha2, beta, nohhldij, nohhldij_diagsum);
        simplekhatcputime = time - simplekhatcputime;
        calclambdacputime = time;
        logL = calcLambdaSimpleKhat(alpha2, beta, DS, DKAOarray, N, distance2, nohhldij, nohhldij_diagsum);
        calclambdacputime = time - calclambdacputime;
        cputime = time - cputime;
        sumcputime += cputime;
        sumsimplekhatcputime += simplekhatcputime;
        sumcalclambdacputime += calclambdacputime;
        sum2cputime += cputime*cputime;
        sum2simplekhatcputime += simplekhatcputime*simplekhatcputime;
        sum2calclambdacputime += calclambdacputime*calclambdacputime;
    }
    cout << "sumcputime over 100 trials " << sumcputime/100 << "s with CIs (" << (sumcputime/100)-(1.96/10*(sqrt((sum2cputime/100)-(pow(sumcputime/100,2.0))))) << "," << (sumcputime/100)+(1.96/10*(sqrt((sum2cputime/100)-(pow(sumcputime/100,2.0))))) << ")" << "\n";
    cout << "sumsimplekhatcputime over 100 trials " << sumsimplekhatcputime/100 << "s with CIs (" << (sumsimplekhatcputime/100)-(1.96/10*(sqrt((sum2simplekhatcputime/100)-(pow(sumsimplekhatcputime/100,2.0))))) << "," << (sumsimplekhatcputime/100)+(1.96/10*(sqrt((sum2simplekhatcputime/100)-(pow(sumsimplekhatcputime/100,2.0))))) << ")" << "\n";
    cout << "sumcalclambdacputime over 100 trials " << sumcalclambdacputime/100 << "s with CIs (" << (sumcalclambdacputime/100)-(1.96/10*(sqrt((sum2calclambdacputime/100)-(pow(sumcalclambdacputime/100,2.0))))) << "," << (sumcalclambdacputime/100)+(1.96/10*(sqrt((sum2calclambdacputime/100)-(pow(sumcalclambdacputime/100,2.0))))) << ")" << "\n";

    int gpusupport = gpu::getCudaEnabledDeviceCount();
    printf("If OpenCV compiled with GPU support then next digit will be a one...%i\n",gpusupport);
    
    gpu::resetDevice(); // the next API call will re-initialise it
    gpu::DeviceInfo x;
    size_t freememory = x.freeMemory();
    cout << "free memory at initialisation = " << freememory/1000000 << "MB\n";

    //convert vector objs to matrices
    static double DS_matrix[4664][4664];
    long k = 0;
    for (long i = 0; i < 4664; i++){
        for (long j = 0; j < 4664; j++)
        {
            if (k < 21743733)
            {
                DS_matrix[i][j] = DS[k];
            }
            else{
                DS_matrix[i][j] = 0;
            }
            k = k + 1;
        }
    }
    //DS_matrix is 21,752,896 but I only fill it with 21743733 so 9163 are zero so need to take off -9163

    static double distance2_matrix[3042][3042];
    k = 0;
    for (int i = 0; i < 3042; i++){
        for (int j = 0; j < 3042; j++)
        {
            if (k < 9250572)
            {
                distance2_matrix[i][j] = distance2[k];
            }
            else{
                distance2_matrix[i][j] = 0;
            }
            k = k + 1;
        }   
    }

    static double nohhldij_matrix[3042][3042];
    k = 0;
    for (int i = 0; i < 3042; i++){
        for (int j = 0; j < 3042; j++)
        {
            if (k < 9250572)
            {
                nohhldij_matrix[i][j] = nohhldij[k];
            }
            else{
                nohhldij_matrix[i][j] = 0;
            }
            k = k + 1;
        }   
    }

    sum = 0;
    for (int i = 0; i < 3042; i++)
    {
        for (int j = 0; j < 3042; j++)
        {
            sum += nohhldij_matrix[i][j];
        }
    }
    printf("nohhldij_matrix sum = %f\n", sum);
    sum = 0;
    for (int i = 0; i < 9250572; i++)
    {
        sum += nohhldij[i];
    }
    printf("nohhldij  sum = %f\n", sum);

    double transfCPUGPUtime = 0;
    double gpucomputetime = 0;
    double transfGPUCPUtime = 0;
    double simplekhatgputime = 0;
    double calclambdasimplekhatgputime = 0;
    
    double sumtransfCPUGPUtime = 0;
    double sumgpucomputetime = 0;
    double sumtransfGPUCPUtime = 0;
    double sumsimplekhatgputime = 0;
    double sumcalclambdasimplekhatgputime = 0;
    double sum2transfCPUGPUtime = 0;
    double sum2gpucomputetime = 0;
    double sum2transfGPUCPUtime = 0;
    double sum2simplekhatgputime = 0;
    double sum2calclambdasimplekhatgputime = 0;

    for (int i = 0; i < 100; ++i) //time trial
    {
        try
        {
        // import data and init host(CPU) _h objs
        Mat alpha2_h(3042,3042,DataType<double>::type,alpha2);
        Mat alpha3_h(4664,4664,DataType<double>::type,alpha2);
        Mat alpha4_h(977,60,DataType<double>::type,alpha2);
        Mat beta_h(1,1,DataType<double>::type,beta);
        Mat DS_h(4664,4664,DataType<double>::type,DS_matrix);
        Mat DKAOarray_h(977,60,DataType<double>::type,DKAOarray60);
        Mat DKAOarraycorr_h(1,977,DataType<double>::type,DKAOarray1);
        Mat N_h(1,1,DataType<double>::type,N);
        Mat distance2_h(3042,3042,DataType<double>::type,distance2_matrix);
        Mat nohhldij_h(3042,3042,DataType<double>::type,nohhldij_matrix);
        Mat nohhldij_diagsum_h(1,1,DataType<double>::type,nohhldij_diagsum);
        Mat result_h(1,1,DataType<double>::type);
        Mat simplekhat_h(1,1,DataType<double>::type);
        Mat DSsum_h(1,1,DataType<double>::type);
        Mat distance2sum_h(1,1,DataType<double>::type);
        Mat nohhldsum_h(1,1,DataType<double>::type);
        Mat DKAOsum_h(1,1,DataType<double>::type);
        Mat S_LL_h(1,1,DataType<double>::type);
        Mat KAO_LL_h(1,1,DataType<double>::type);
        Mat rowtotal_h(977,1,DataType<double>::type);
        
        // cout << "\nGPU\nalpha2 " << alpha2 << "\n";
        // cout << "beta_h " << beta_h << "\n";
        // cout << "N_h " << N_h << "\n";
        // cout << "nohhldij_diagsum " << nohhldij_diagsum_h << "\n";

        gpu::GpuMat beta_g(1,1,DataType<double>::type);
        gpu::GpuMat simplekhat_g(1,1,DataType<double>::type);
        gpu::GpuMat N_g(1,1,DataType<double>::type);
        gpu::GpuMat nohhldij_diagsum_g(1,1,DataType<double>::type);
        gpu::GpuMat result_g(3042,3042,DataType<double>::type);
        gpu::GpuMat alpha2_g(3042,3042,DataType<double>::type);
        gpu::GpuMat nohhldij_g(3042,3042,DataType<double>::type);
        gpu::GpuMat distance2_g(3042,3042,DataType<double>::type);
        gpu::GpuMat alpha3_g(4664,4664,DataType<double>::type);
        gpu::GpuMat DS_g(4664,4664,DataType<double>::type);
        gpu::GpuMat result2_g(4664,4664,DataType<double>::type);
        gpu::GpuMat alpha4_g(977,60,DataType<double>::type);
        gpu::GpuMat DKAOarray_g(977,60,DataType<double>::type);
        gpu::GpuMat DKAOarraycorr_g(1,977,DataType<double>::type);
        gpu::GpuMat result3_g(977,60,DataType<double>::type);
        gpu::GpuMat DSsum_g(1,1,DataType<double>::type);
        gpu::GpuMat distance2sum_g(1,1,DataType<double>::type);
        gpu::GpuMat nohhldijsum_g(1,1,DataType<double>::type);
        gpu::GpuMat DKAOsum_g(1,1,DataType<double>::type);

        // //copy objs host->GPU
        transfCPUGPUtime = time;
        beta_g.upload(beta_h);
        N_g.upload(N_h);
        nohhldij_diagsum_g.upload(nohhldij_diagsum_h);
        alpha2_g.upload(alpha2_h);
        nohhldij_g.upload(nohhldij_h);
        distance2_g.upload(distance2_h);
        alpha3_g.upload(alpha3_h);
        alpha4_g.upload(alpha4_h);
        DS_g.upload(DS_h);
        DKAOarray_g.upload(DKAOarray_h);
        DKAOarraycorr_g.upload(DKAOarraycorr_h);
        result3_g.upload(alpha4_h); //sic
        transfCPUGPUtime = time - transfCPUGPUtime;
        //cout << "transfer time host -> GPU = " << transfCPUGPUtime << "s \n";

        //GPU sum checks
        DSsum_g = gpu::sum(DS_g);
        distance2sum_g = gpu::sum(distance2_g);
        nohhldijsum_g = gpu::sum(nohhldij_g);
        DKAOsum_g = gpu::sum(DKAOarray_g);
        DSsum_g.download(DSsum_h);
        distance2sum_g.download(distance2sum_h);
        nohhldijsum_g.download(nohhldsum_h);
        DKAOsum_g.download(DKAOsum_h);
        // cout << "DSsum_h " << DSsum_h << "\n";
        // cout << "distance2sum_h " << distance2sum_h << "\n";
        // cout << "nohhldsum_h " << nohhldsum_h << "\n";
        // cout << "DKAOsum_h " << DKAOsum_h << "\n";

        //simpleKhatBeta(){}
        gpucomputetime = time;
        simplekhatgputime = time;
        gpu::add(distance2_g,alpha2_g,result_g);
        gpu::divide(nohhldij_g,result_g,result_g); 
        gpu::multiply(result_g,alpha2_g,result_g);
        simplekhat_g = gpu::sum(result_g);
        gpu::add(simplekhat_g,nohhldij_diagsum_g,simplekhat_g);
        gpu::divide(N_g,simplekhat_g,simplekhat_g);
        gpu::multiply(simplekhat_g,beta_g,simplekhat_g);
        double simplekhat_g_scalar = gpu::sum(simplekhat_g)[0]; //yuck but it works to output a Scalar object. unclear if cv::Scalar objs can be created on the GPU or if this involves GPU<>host communication
        simplekhatgputime = time - simplekhatgputime;
        //cout << "simplekhatBeta computation time on GPU = " << simplekhatgputime << "s \n";
        simplekhat_g.download(simplekhat_h);

        //calcLambdaSimpleKhat(){}
        calclambdasimplekhatgputime = time;
        gpu::add(DS_g,alpha3_g,result2_g);
        gpu::divide(alpha3_g,result2_g,result2_g);   
        double result2_g_sum = gpu::sum(result2_g)[0];
        gpu::GpuMat S_LL_g(1,1,DataType<double>::type,result2_g_sum);
        gpu::subtract(S_LL_g,9163,S_LL_g); //correction to DS for padded zeroes
        gpu::multiply(S_LL_g,simplekhat_g,S_LL_g); 
        S_LL_g.download(S_LL_h);

        gpu::GpuMat rowtotal_g(977,1,DataType<double>::type,0.0);
        gpu::GpuMat one_g(1,977,DataType<double>::type,1.0);
        gpu::add(DKAOarray_g,result3_g,result3_g);
        gpu::divide(alpha4_g,result3_g,result3_g);     
        gpu::reduce(result3_g,rowtotal_g,1,CV_REDUCE_SUM,-1);
        gpu::add(rowtotal_g,DKAOarraycorr_g,rowtotal_g);
        gpu::multiply(rowtotal_g,simplekhat_g_scalar,rowtotal_g);
        gpu::exp(rowtotal_g,rowtotal_g);
        gpu::subtract(one_g,rowtotal_g,rowtotal_g);
        gpu::log(rowtotal_g,rowtotal_g);
        double KAO_LL_g_scalar = gpu::sum(rowtotal_g)[0]; //yuck
        gpu::GpuMat KAO_LL_g(1,1,DataType<double>::type,KAO_LL_g_scalar);
        KAO_LL_g.download(KAO_LL_h);
        gpu::add(S_LL_g,KAO_LL_g,S_LL_g);
        calclambdasimplekhatgputime = time - calclambdasimplekhatgputime;
        gpucomputetime = time - gpucomputetime;
        // cout << "calcLambdaSimpleKhat computation time on GPU = " << calclambdasimplekhatgputime << "s \n";
        // cout << "total gpucomputetime on GPU = " << gpucomputetime << "s \n";

        transfGPUCPUtime = time;
        S_LL_g.download(result_h);
        transfGPUCPUtime = time - transfGPUCPUtime;
        // cout << "transfer time GPU -> host = " << transfGPUCPUtime << "s \n";
        // cout << "simplekhat_h(from GPU) = " << simplekhat_h << "\n";
        // cout << "S_LL_h(from GPU) = " << S_LL_h << "\n";
        // cout << "KAO_LL_h(from GPU) = " << KAO_LL_h << "\n";
        // cout << "result_h(from GPU) = " << result_h << "\n";

        freememory = x.freeMemory();
        //cout << "free memory afterwards = " << freememory/1000000 << "MB\n";

        // // float data[3][1] = {{1.3},{1.3},{1.3}};
        // // float data2[3][3] = {{1,2,3},{4,5,6},{7,8,9}};
        // // Mat vector(3,1,DataType<float>::type,data);
        // // Mat array2(3,3,DataType<float>::type,data2);
        
        // //init objs on host
        // // Mat vector_host;
        // // Mat vector_add;
        // // Mat vectorscal_mul;
        // // Mat vector_exp;
        // // Mat matrix_host;
        // // Mat vector_subtract;
        // // Mat vector_log;
        // // Mat vector_sum;
        // // Mat matrix_rowsumreduce;
        // // Mat matrix_add;
        // // Mat matrixscal_mul;
        // // Mat matrix_divide;
        // // Mat matrix_sum;
        // // Mat scalar_prod;

        // // //print input array
        // // cout << "vector(init on host) = " << endl << " " << vector << endl;
        // // cout << "array2(init on host) = " << endl << " " << array2 << endl;
        
        // //init objs on GPU
        // // gpu::GpuMat vector_gpu(3,1,DataType<float>::type);
        // // gpu::GpuMat vectoradd_gpu(3,1,DataType<float>::type);
        // // gpu::GpuMat vectorscalmul_gpu(3,1,DataType<float>::type);
        // // gpu::GpuMat vectorexp_gpu(3,1,DataType<float>::type);
        // // gpu::GpuMat vectorsubtr_gpu(3,1,DataType<float>::type);
        // // gpu::GpuMat vectorlog_gpu(3,1,DataType<float>::type);
        // // gpu::GpuMat vectorsum_gpu(1,1,DataType<float>::type); //note even scalars are 'matrix' types as .download() only works on GpuMat types
        // // gpu::GpuMat matrix_gpu(3,3,DataType<float>::type); 
        // // gpu::GpuMat matrixrowsumreduce_gpu(3,1,DataType<float>::type); //note this is a vector to receive the result
        // // gpu::GpuMat matrixadd_gpu(3,3,DataType<float>::type);
        // // gpu::GpuMat matrixscalmul_gpu(3,1,DataType<float>::type);
        // // gpu::GpuMat matrixdiv_gpu(3,3,DataType<float>::type);
        // // gpu::GpuMat matrixsum_gpu(1,1,DataType<float>::type);
        // // gpu::GpuMat scalarprod_gpu(1,1,DataType<float>::type);
        // // scalarprod_gpu = 3.14;

        // //copy objs host->GPU
        // // a = time;
        // // vector_gpu.upload(vector);
        // // matrix_gpu.upload(array2);
        // // cout << "transfer time host -> GPU = " << time - a << "s \n";
        
        // // //do GPU operations
        // // a = time;
        // // gpu::add(vector_gpu,vector_gpu,vectoradd_gpu);
        // // gpu::multiply(vector_gpu,3.01,vectorscalmul_gpu);
        // // gpu::exp(vector_gpu,vectorexp_gpu);
        // // gpu::subtract(vector_gpu,0.09,vectorsubtr_gpu);
        // // gpu::log(vector_gpu,vectorlog_gpu);
        // // vectorsum_gpu = gpu::sum(vector_gpu);
        // // gpu::reduce(matrix_gpu,matrixrowsumreduce_gpu,1,CV_REDUCE_SUM,-1);
        // // gpu::add(matrix_gpu,matrix_gpu,matrixadd_gpu);
        // // gpu::multiply(matrix_gpu,2.4,matrixscalmul_gpu);
        // // gpu::divide(matrix_gpu,matrix_gpu,matrixdiv_gpu);
        // // matrixsum_gpu = gpu::sum(matrix_gpu);
        // // gpu::multiply(scalarprod_gpu,scalarprod_gpu,scalarprod_gpu); //note overwrites scalarprod_gpu
        // // cout << "computation time on GPU = " << time - a << "\n";

        // // //copy objs GPU->host
        // // a = time;
        // // vector_gpu.download(vector_host);
        // // vectoradd_gpu.download(vector_add);
        // // vectorscalmul_gpu.download(vectorscal_mul);
        // // vectorexp_gpu.download(vector_exp);
        // // vectorsubtr_gpu.download(vector_subtract);
        // // vectorlog_gpu.download(vector_log);
        // // vectorsum_gpu.download(vector_sum);
        // // matrixrowsumreduce_gpu.download(matrix_rowsumreduce);
        // // matrixadd_gpu.download(matrix_add);
        // // matrixscalmul_gpu.download(matrixscal_mul);
        // // matrixdiv_gpu.download(matrix_divide);
        // // matrixsum_gpu.download(matrix_sum);
        // // scalarprod_gpu.download(scalar_prod);
        // // cout << "transfer time GPU -> host = " << time - a << "\n";

        // // //print output array
        // // cout << "vector_host(unaltered (vector) back from GPU) = " << endl << " " << vector_host << endl;
        // // cout << "vector_add(result of GPU addition) = " << endl << " " << vector_add << endl;
        // // cout << "vectorscal_mul(result of GPU V*3.01[scalar]) = " << endl << " " << vectorscal_mul << endl;
        // // cout << "vector_exp(result of exp(vector)) = " << endl << " " << vector_exp << endl;
        // // cout << "vector_subtr(result of (vector-0.09)) = " << endl << " " << vector_subtract << endl;
        // // cout << "vector_log(result of natural log(vector)) = " << endl << " " << vector_log << endl;
        // // cout << "vector_sum(result of sum(vector)) = " << endl << " " << vector_sum << endl;
        // // cout << "matrix_rowsumreduce(result of rowsum(array2)) = " << endl << " " << matrix_rowsumreduce << endl;
        // // cout << "matrix_add(result of array2+array2) = " << endl << " " << matrix_add << endl;
        // // cout << "matrixscal_mul(result of array2*2.4) = " << endl << " " << matrixscal_mul << endl;
        // // cout << "matrix_divide(result of array2/array2) = " << endl << " " << matrix_divide << endl;
        // // cout << "matrix_sum(result of sum(array2)) = " << endl << " " << matrix_sum << endl;
        // // cout << "scalar_prod(result of 3.14*3.14) = " << endl << " " << scalar_prod << endl;
        sumtransfCPUGPUtime += transfCPUGPUtime;
        sumgpucomputetime += gpucomputetime;
        sumtransfGPUCPUtime += transfGPUCPUtime;
        sumsimplekhatgputime += simplekhatgputime;
        sumcalclambdasimplekhatgputime += calclambdasimplekhatgputime;
        sum2transfCPUGPUtime += transfCPUGPUtime*transfCPUGPUtime;
        sum2gpucomputetime += gpucomputetime*gpucomputetime;
        sum2transfGPUCPUtime += transfGPUCPUtime*transfGPUCPUtime;
        sum2simplekhatgputime += simplekhatgputime*simplekhatgputime;
        sum2calclambdasimplekhatgputime += calclambdasimplekhatgputime*calclambdasimplekhatgputime;
        }

        catch(const cv::Exception& ex)
        {
            cout << "Error: " << ex.what() << "Reason: " << ex.msg << endl;
        }
    }
    cout << "sumtransfCPUGPUtime over 100 trials " << sumtransfCPUGPUtime/100 << "s with CIs (" << (sumtransfCPUGPUtime/100)-(1.96/10*(sqrt((sum2transfCPUGPUtime/100)-(pow(sumtransfCPUGPUtime/100,2.0))))) << "," << (sumtransfCPUGPUtime/100)+(1.96/10*(sqrt((sum2transfCPUGPUtime/100)-(pow(sumtransfCPUGPUtime/100,2.0))))) << ")" << "\n";
    cout << "sumgpucomputetime over 100 trials " << sumgpucomputetime/100 << "s with CIs (" << (sumgpucomputetime/100)-(1.96/10*(sqrt((sum2gpucomputetime/100)-(pow(sumgpucomputetime/100,2.0))))) << "," << (sumgpucomputetime/100)+(1.96/10*(sqrt((sum2gpucomputetime/100)-(pow(sumgpucomputetime/100,2.0))))) << ")" << "\n";
    cout << "sumtransfGPUCPUtime over 100 trials " << sumtransfGPUCPUtime/100 << "s with CIs (" << (sumtransfGPUCPUtime/100)-(1.96/10*(sqrt((sum2transfGPUCPUtime/100)-(pow(sumtransfGPUCPUtime/100,2.0))))) << "," << (sumtransfGPUCPUtime/100)+(1.96/10*(sqrt((sum2transfGPUCPUtime/100)-(pow(sumtransfGPUCPUtime/100,2.0))))) << ")" << "\n";
    cout << "sumsimplekhatgputime over 100 trials " << sumsimplekhatgputime/100 << "s with CIs (" << (sumsimplekhatgputime/100)-(1.96/10*(sqrt((sum2simplekhatgputime/100)-(pow(sumsimplekhatgputime/100,2.0))))) << "," << (sumsimplekhatgputime/100)+(1.96/10*(sqrt((sum2simplekhatgputime/100)-(pow(sumsimplekhatgputime/100,2.0))))) << ")" << "\n";
    cout << "sumcalclambdasimplekhatgputime over 100 trials " << sumcalclambdasimplekhatgputime/100 << "s with CIs (" << (sumcalclambdasimplekhatgputime/100)-(1.96/10*(sqrt((sum2calclambdasimplekhatgputime/100)-(pow(sumcalclambdasimplekhatgputime/100,2.0))))) << "," << (sumcalclambdasimplekhatgputime/100)+(1.96/10*(sqrt((sum2calclambdasimplekhatgputime/100)-(pow(sumcalclambdasimplekhatgputime/100,2.0))))) << ")" << "\n";

    free(DS);
    free(distance2);
    free(nohhldij);
    free(DKAO);
    return 0;
}