#ifndef GA_GPU_H
#define GA_GPU_H

#include <iostream>
#include <thread>
#include <complex>
#include <cv.h>
#include <stdio.h>
#include <stdlib.h>
#include <iostream>

#include "timer.h"

#include "basic_functions.h"
//LAPACKE support for Visual Studio

#ifndef LAPACK_COMPLEX_CUSTOM
#define LAPACK_COMPLEX_CUSTOM
#define lapack_complex_float std::complex<float>
#define lapack_complex_double std::complex<double>
#include "lapacke.h"
#endif


#define LAPACK_ROW_MAJOR               101
#define LAPACK_COL_MAJOR               102

//CUDA functions
void gpuIntialization(unsigned int** gpuP, size_t p, size_t f,													//variables required for the population allocation
	float** gpuCM, float* cpuCM, size_t nC, unsigned int ub,
	float** gpuM, float* cpuM, unsigned int** gpu_nPxInCls,
	float** gpuSb, float** gpuSw,
	float** gpuF, float* cpuF,
	unsigned int** gpuT, unsigned int* cpuT, size_t tP, unsigned int* cpu_nPxInCls);
void gpucomputeSbSw(unsigned int* gpuP, unsigned int* cpuP, size_t p, size_t f,
	float* gpuSb, float* cpuSb,
	float* gpuSw, float* cpuSw,
	float* gpuF, unsigned int* T, float* gpuM, float* gpuCM,
	size_t nC, size_t tP, cudaDeviceProp props, size_t gen, size_t gnrtn, size_t ub, unsigned int* gpu_nPxInCls, std::ofstream& profilefile);
void gpuDestroy(unsigned int* gpuP, float* gpuCM, float* gpuM, unsigned int* gpu_nPxInCls, float* gpuSb, float* gpuSw, float* gpuF, unsigned int* gpuT);

struct _fcomplex { float re, im; };
typedef struct _fcomplex fcomplex;

Timer timer;

class ga_gpu {

public:
	float* F;					//pointer to the raw data in host memory
	unsigned int* T;			//pointer to the class labels in host memory
	size_t gnrtn;				//total number of generations
	size_t p;					//population size
	size_t f;					// number of features to be selected

	unsigned int* P;			//pointer to population of current generation genotype matrix (p x f)
	float* S;					//pointer to score(fitness value) of each gnome from current population matric P
	unsigned int* i_guess;		//initial guess of features if mentioined in args add to initial population
	unsigned int ub;			//upper bound for gnome value (maximum feature index from raw feature matrix F)
	unsigned int lb;			//lower bound for gnome value (minimum feature index from raw feature matrix F = 0)
	float uniformRate;
	float mutationRate;
	size_t tournamentSize;  	//number of potential gnomes to select parent for crossover
	bool elitism;        	//if true then passes best gnome to next generation

							//declare gpu pointers
	float* gpuF;				//Feature matrix	
	unsigned int* gpuT;			//target responses of entire feature matrix	
	unsigned int* gpuP;		//population matrix
	unsigned int* gpu_nPxInCls;
	float* gpuCM;				//class mean of entire feature matrix
	float* gpuM;				//total mean of entire feature matrix
	float* gpuSb;				//between class scatter for all individuals of current population
	float* gpuSw;				//within class scatter for all individuals of current population

								//constructor
	ga_gpu() {}

	//==============================generate initial population

		void initialize_population(std::vector<unsigned int> i_guess, bool debug) {
			if (debug) {
				std::cout << std::endl;
				std::cout << "initial populatyion is: " << std::endl;
			}
			
			lb = 0;
			P = (unsigned int*)calloc(p * f, sizeof(unsigned int));					//allcate memory for genetic population(indices of features from F), p number of gnomes of size f
			S = (float*)calloc(p, sizeof(float));							//allcate memory for scores(fitness value) of each gnome from P 

			srand(1);
			//add intial guess to the population if specified by user as a output of other algorithm or by default just random guess
			std::memcpy(P, i_guess.data(), f * sizeof(unsigned int));

			//generate random initial population
			for (size_t i1 = 1; i1 < p; i1++) {
				for (size_t i2 = 0; i2 < f; i2++) {
					P[i1 * f + i2] = rand() % ub + lb;						//select element of gnome as random feature index within lower bound(0) and upper bound(B) 
					if (debug)	std::cout << P[i1 * f + i2] << "\t";
				}
				if (debug) std::cout << std::endl;
			}
		}

		//===================generation of new population==========================================

		size_t evolvePopulation(unsigned int* newPop, float* M, bool debug) {

			//gget index of best gnome in the current population
			size_t bestG_Indx = gIdxbestGnome();
			//-------------(reproduction)-------
			if (elitism) {
				saveGnomeIdx(0, bestG_Indx, newPop);							//keep best gnome from previous generation to new generation
			}
			// ------------Crossover population---------------
			int elitismOffset;
			if (elitism) {
				elitismOffset = 1;
			}
			else {
				elitismOffset = 0;
			}

			//Do crossover for rest of population size
			for (int i = elitismOffset; i <p; i++) {
				//	std::cout<<"crossover of gnome "<<i<<std::endl;
				std::vector<unsigned int>gnome1;
				gnome1.reserve(f);
				gnome1 = tournamentSelection(5);								//select first parent for crossover from tournament selection of 5 gnomes
																				//	displaygnome(gnome1);
				std::vector<unsigned int>gnome2;
				gnome2.reserve(f);
				gnome2 = tournamentSelection(5);								//select first parent for crossover from tournament selection of 5 gnomes
																				//	displaygnome(gnome2);
				std::vector<unsigned int>gnome;
				gnome.reserve(f);
				gnome = crossover(gnome1, gnome2, M);								//Do crossover of above parent gnomes to produce new gnome
																					//	displaygnome(gnome);
				saveGnome(i, gnome, newPop);									//save crosseover result to new population
			}

			//--------------Mutate population------------
			// introduce some mutation in new population
			for (int i = elitismOffset; i <p; i++) {
				//std::cout<<"mutation of gnome"<<std::endl;
				std::vector<unsigned int>gnome;
				gnome.reserve(f);

				for (size_t n = 0; n < f; n++)
					gnome.push_back(newPop[i*f + n]);
				//std::cout<<"\n starting address "<<(&newPop[0] + i*f)<<"\t end address is "<<(&newPop[0] + i*f + f-1) <<std::endl;
				//std::copy((&newPop[0] + i*f), (&newPop[0] + i*f +f-1), gnome.begin());
				//	displaygnome(gnome);
				mutate(gnome);
				//	displaygnome(gnome);
				saveGnome(i, gnome, newPop);								//save new gnome to new population at position i
			}
			return bestG_Indx;
		}

		//============================== functions for population evolution ===========================================================================
		std::vector<unsigned int> tournamentSelection(size_t tSize) {
			// Create a tournament population
			unsigned int* tournamentP = (unsigned int*)malloc(tSize * f * sizeof(unsigned int));
			std::vector<float>tournamentS;

			// For each place in the tournament get a random individual
			for (size_t i = 0; i < tSize; i++) {
				size_t rndmIdx = rand() % p + lb;
				tournamentS.push_back(S[rndmIdx]);
				//for (size_t n = 0; n <f; n++) 
				//tournamentP[i * f + n] = (getGnome(rndmIdx)).at(n);
				std::vector<unsigned int> temp_g(getGnome(rndmIdx));
				std::copy(temp_g.begin(), temp_g.end(), tournamentP + i*f);
			}
			// Get the fittest
			std::vector<unsigned int>fittestgnome;
			fittestgnome.reserve(f);

			//select index of best gnome from fitness score
			size_t bestSIdx = 0;
			for (size_t i = 0; i < tSize; i++) {
				if (tournamentS[i] < tournamentS[bestSIdx])
					bestSIdx = i;    //float check : it was like this b(&idx[i], &idx[j]) but gave me error
			}

			for (size_t n = 0; n < f; n++)
				fittestgnome.push_back(tournamentP[bestSIdx * f + n]);
			return fittestgnome;
		} //end of tournament selection 


		std::vector<unsigned int>  crossover(std::vector<unsigned int>  gnome1, std::vector<unsigned int>  gnome2, float* M) {
			std::vector<unsigned int> gnome;
			for (size_t i = 0; i < f; i++) {
				// Crossover
				float r = static_cast <float> (rand()) / static_cast <float> (RAND_MAX);
				if (r <= uniformRate) {
					gnome.push_back(gnome1.at(i));
				}
				else {
					gnome.push_back(gnome2.at(i));
				}
			}

			//check new gnome for all zero bands and duplicated values
			std::vector<unsigned int> gnomeunique;
			int flag = 0;
			std::sort(gnome.begin(), gnome.end()); // 1 1 2 2 3 3 3 4 4 5 5 6 7 
			std::unique_copy(gnome.begin(), gnome.end(), std::back_inserter(gnomeunique));
			/*	if(gnomeunique.size()< gnome.size()){
			flag = 1;
			std::cout<<"gnome:["<<g<<"] "<<"\t duplications are "<< (gnome.size() - gnomeunique.size())<<std::endl;
			}*/
			unsigned int featureband, featureband1, featureband2;
			if (gnomeunique.size() < f) {
				for (size_t k = gnomeunique.size(); k < f; k++) {
					featureband = rand() % ub + lb;
					for (size_t i = 0; i < f; i++) {
						featureband1 = gnome1.at(i);
						featureband2 = gnome2.at(i);
						for (size_t j = 0; j < gnomeunique.size(); j++) {
							if (gnomeunique.at(j) != featureband1) {
								featureband = featureband1;
							}
							else if (gnomeunique.at(j) != featureband2) {
								featureband = featureband2;
							}
							else if (gnomeunique.at(j) == featureband) {
								featureband = rand() % ub + lb;
								while (M[featureband] == 0) {
									featureband = rand() % ub + lb;
								}
							}
						}
					}
					gnomeunique.push_back(featureband);
				}
			}
			//if(flag ==1){
			//	std::cout<<"\n original gnome "<<g<<" are "<<std::endl;
			//	for(int k = 0; k < gnome.size(); k++)
			//		std::cout<<gnome[k]<<"\t";
			//	std::cout<<"\n unique results in cpp for gnome "<<g<<" are "<<std::endl;
			//	for(int k = 0; k < gnomeunique.size(); k++)
			//		std::cout<<gnomeunique[k]<<"\t";
			//}

			return gnomeunique;
		}

		void mutate(std::vector<unsigned int> gnome) {
			for (size_t i = 0; i < f; i++) {
				float LO = (float)0.01;
				float HI = 1;
				float r3 = LO + static_cast <float> (rand()) / (static_cast <float> (RAND_MAX / (HI - LO)));
				//if random value is less than mutationRate then mutate this gnome
				if (r3 <= mutationRate) {
					gnome.at(i) = (rand() % ub + lb);
					gnome.push_back(rand() % ub + lb);
				}
			}
		}

		///returns gnome of given index
		std::vector<unsigned int> getGnome(size_t idx) {
			std::vector<unsigned int> gnome;
			gnome.reserve(f);
			//pulling gnome idx from population P
			for (size_t n = 0; n < f; n++)
				gnome.push_back(P[idx * f + n]);
			//memcpy(&gnome[0], P+idx*f, f*sizeof(size_t));
			return gnome;
		}

		//save gnome of index gIdx from previous population at position i in the new population
		void saveGnomeIdx(size_t i, size_t gIdx, unsigned int* newPop) {
			for (size_t n = 0; n < f; n++)
				newPop[i * f + n] = P[gIdx * f + n];
		}

		void saveGnome(size_t idx, std::vector<unsigned int>gnome, unsigned int* newPop) {
			std::copy(gnome.begin(), gnome.end(), newPop + idx*f);
		}

		size_t gIdxbestGnome() {
			//std::cout<<"best gnome indes is: "<<sortSIndx()[0];
			return sortSIndx()[0];
		}

		void displaygnome(std::vector<unsigned int> gnome) {
			std::cout << "\t gnome: ";
			for (int i = 0; i<gnome.size(); ++i)
				std::cout << gnome[i] << ' ';
			std::cout << std::endl;
		}

		//---------------------post processing of score-------------------------------------
		void Snorm() { //normalize gnome scores 
			double s;
			for (size_t i = 0; i < p; i++) {
				s += S[i];			//sum of all gnome score in population
			}
			//std::cout<<"mean Score is: "<<(double) s/p;
			for (size_t i = 0; i <p; i++)
				S[i] = S[i] / s;
		}

		size_t* sortSIndx() {     //sort gnome index according to gnome scores
			//sort indices of score in ascending order (fitness value)
			size_t *idx = (size_t*)malloc(p * sizeof(size_t));	//array to hold sorted gnome index
			for (size_t i = 0; i < p; i++) {					//initialize index array from 1 to p(population size) in an ascending order
				idx[i] = i;
			}

			for (size_t i = 0; i<p; i++) {				//sort gnome indices according to score values using bubble sort 
				for (size_t j = i + 1; j<p; j++) {
					if (S[idx[i]] > S[idx[j]]) {
						std::swap(idx[i], idx[j]);		//float check : it was like this b(&idx[i], &idx[j]) but gave me error
					}
				}
			}

			//display best gnome
            //std::cout << "best fitness value: " << S[idx[0]] << std::endl;
			/*if (S[idx[0]] < 0) {				
				std::cout << "best gnome is " << std::endl;
				for (size_t i = 0; i < f; i++)
					std::cout << P[f * idx[0] + i] << ", ";
				std::cout << std::endl;
			}*/

			return idx; //use as sortSIdx in selection
		}


		//size_t* sortIndx(float* input, size_t size) {
		//	//sort indices of score in ascending order (fitness value)
		//	size_t *idx;
		//	idx = (size_t*)malloc(size * sizeof(size_t));
		//	for (size_t i = 0; i < size; i++)
		//		idx[i] = i;

		//	for (size_t i = 0; i<size; i++) {
		//		for (size_t j = i + 1; j<size; j++) {
		//			if (input[idx[i]] < input[idx[j]]) {
		//				std::swap(idx[i], idx[j]);				   //float check : it was like this b(&idx[i], &idx[j]) but gave me error
		//			}
		//		}
		//	}
		//	return idx; //use as sortSIdx in selection

		//}

		void generateNewP(unsigned int* newPop) {
			//std::memcpy(P, 0 , p * f *sizeof(unsigned int));			   //copy sb of gnome 'g' into bufferarray tempg_s	
			std::memcpy(P, newPop, p * f * sizeof(unsigned int));		   //copy sb of gnome 'g' into bufferarray tempg_s	
		}

		//============================== functions for fitness function  ===========================================================================
		//compute total mean M (1 X B) of all features (tP X B)
		void ttlMean(float* M, size_t tP, size_t B) {
			//std::cout<<"total number of pixels are "<<tP<<std::endl;
			for (int k = 0; k < tP; k++) {								  //total number of pixel in feature matrix
				for (size_t n = 0; n < B; n++) {							  // index of feature in ith gnome
					M[n] += F[k * B + n];
				}
			}
			for (size_t n = 0; n < B; n++)								 //take an avarage of above summation 
				M[n] = M[n] / (float)tP;
		}

		void dispalymean(float* M) {				//display mean
			std::cout << std::endl;
			std::cout << "Total mean of gnome 1 features are is " << std::endl;

			for (size_t i = 0; i < 1; i++) {
				for (size_t j = 0; j < f; j++) {
					size_t index = P[i*f + j];
					std::cout << "feature index " << index << "\t total mean" << M[index] << std::endl;
				}
			}
			std::cout << std::endl;
		}

		//Compute class means cM (p x nC x f) of all gnome features phe(tP x f)
		void classMean(float* cM, size_t tP, size_t nC, size_t B, std::vector<unsigned int> nPxInCls) {
			for (size_t c = 0; c < nC; c++) {									//index of class feature matrix responses
				float* tempcM = (float*)calloc(B, sizeof(float));			//tempcM holds classmean vector for current gnome 'i', class 'c'
				for (size_t k = 0; k < tP; k++) {							//total number of pixel in feature matrix
					if (T[k] == c + 1) {									//class numbers start from 1 not 0
						for (size_t n = 0; n < B; n++) {					//total number of features in a gnome
							tempcM[n] += F[k * B + n];					//add phe value for feature n of class 'c' in  ith gnome
						}
					}
				}
				for (size_t n = 0; n < B; n++)
					cM[c * B + n] = tempcM[n] / (float)nPxInCls[c];			//divide by number of pixels from class 'c'

			}

		}

		//display class mean
		void dispalyClassmean(float* cM, size_t nC) {
			std::cout << std::endl;
			std::cout << "class mean of gnome 1 with total classes " << nC << " is :" << std::endl;
			for (size_t i = 0; i < 1; i++) {
				for (size_t c = 0; c < nC; c++) {
					for (size_t j = 0; j < f; j++) {
						size_t index = P[i*f + j];

						std::cout << "class index: " << c << "\t feature index " << index << "\t  class mean " << cM[c * ub + index] << std::endl;
					}
				}
			}
			std::cout << std::endl;
		}

		//-----------------------------------------between and within class Scattering computation---------------------------------------------------------------
		//computation on CPU
		void cpu_computeSbSw(float* sb, float* sw, float* M, float* cM, size_t nC, size_t tP, std::vector<unsigned int> nPxInCls) {
			timer.start();
			computeSb(sb, M, cM, nC, nPxInCls);					//compute between class scatter on CPU
			const auto elapsed = timer.time_elapsed();
			std::cout << "Sb CPU time " << std::chrono::duration_cast<std::chrono::microseconds>(elapsed).count() << "us" << std::endl;

			timer.start();
			computeSw(sw, cM, nC, tP);							//compute within class scatter on CPU
			const auto elapsed1 = timer.time_elapsed();
			std::cout << "Sw CPU time " << std::chrono::duration_cast<std::chrono::microseconds>(elapsed1).count() << "us" << std::endl;
		}

		//display between class scatter
		void displaySb(float* sb) {
			std::cout << "between scatter is  " << std::endl;
			for (size_t g = 0; g<1; g++) {
				std::cout << std::endl;
				for (size_t j = 0; j < f; j++) {								//total number of features in a gnome
					for (size_t k = 0; k < f; k++) {							//total number of features in a gnome	
						std::cout << sb[g * f * f + j * f + k] << "  ";
					}
					std::cout << std::endl;
				}
			}
			std::cout << std::endl;
		}

		//Compute between class scatter sb (p x f x f) of all gnome features phe(tP x f)
		void computeSb(float* sb, float* M, float* cM, size_t nC, std::vector<unsigned int> nPxInCls) {
			float tempsbval;
			size_t n1;
			size_t n2;
			size_t classIndx;										//class index in class mean matrix
			/*std::cout <<"population of computation of cpusb "<< std::endl;
			for (size_t i2 = 0; i2 < f; i2++) {
				std::cout << P[i2] << "\t";
			}*/

			for (size_t gnomeIndx = 0; gnomeIndx < p; gnomeIndx++) {
				for (size_t c = 0; c < nC; c++) {
					for (size_t i = 0; i < f; i++) {
						for (size_t j = 0; j < f; j++) {
							tempsbval = 0;
							classIndx = c * ub;
							n1 = P[gnomeIndx * f + i];						//actual feature index in original feature matrix 
							n2 = P[gnomeIndx * f + j];						//actual feature index in original feature matrix 
						//	std::cout << "i: " << i << " j: " <<j<< " n1: " << n1 << " n2:" << n2 << std::endl;
							tempsbval = ((cM[classIndx + n1] - M[n1]) *(cM[classIndx + n2] - M[n2]));							
							sb[gnomeIndx * f * f + i * f + j] += tempsbval * (float)nPxInCls[c];  // compute tempsb[j][k] element of class 'c' of gnome 'i'
						}
					}
				}
			}

		}

		//Compute within class scatter sw (p x f x f) of all gnome features phe(tP x f)
		void computeSw(float* sw, float* cM, size_t nC, size_t tP) {
			float tempswval;
			size_t n1;
			size_t n2;
			size_t cMclass;										//class index in class mean matrix
			size_t Pg;
			size_t swg;
			size_t pheg;
			for (size_t gnomeIndx = 0; gnomeIndx < p; gnomeIndx++) {
				Pg = gnomeIndx * f;
				swg = gnomeIndx * f * f;
				pheg = gnomeIndx * tP * f;;
				for (size_t c = 0; c < nC; c++) {
					cMclass = c * ub;

					for (size_t k = 0; k < tP; k++) {
						if (T[k] == (c + 1)) {
							for (size_t i = 0; i < f; i++) {
								for (size_t j = 0; j < f; j++) {
									n1 = P[Pg + i];						//actual feature index in original feature matrix 
									n2 = P[Pg + j];						//actual feature index in original feature matrix 

									tempswval = 0;
									tempswval = ((F[k * ub + n1] - cM[cMclass + n1]) * (F[k * ub + n2] - cM[cMclass + n2]));
									//tempswval = ((phe[gnomeIndx * tP * f + k * f + i] - cM[c * ub + P[gnomeIndx * f + i]]) * (phe[gnomeIndx * tP *f + k * f + j] - cM[c * ub + P[gnomeIndx * f + j]]));	
									sw[gnomeIndx * f * f + i * f + j] += tempswval;
								}
							}
						}
					}
				}

			}
		}
		//checking bands with all zeros and replacing duplicated bands in gnome but this function is only for initial population
		//void zerobandcheck(float* M, bool initial) {
		//	for (size_t g = 0; g < p; g++) {										// for each gnome			
		//		for (size_t i = 0; i < f; i++) {									//check each band (feature) index in that gnome
		//			while (M[P[g * f + i]] == 0) {									//if mean of band is zero then replace band index in population
		//				P[g * f + i] = rand() % ub + lb;
		//			}
		//		}
		//		//checking for duplicats in a gnome
		//		std::vector<unsigned int> gnome = getGnome(g);
		//		std::vector<unsigned int> gnomeunique;
		//		int flag = 0;							//flag will be set if gnome has duplicated band (feature) index
		//		std::sort(gnome.begin(), gnome.end());	// 1 1 2 2 3 3 3 4 4 5 5 6 7 
		//		std::unique_copy(gnome.begin(), gnome.end(), std::back_inserter(gnomeunique));		//keep only unique copies of indices and remove duplicate copies
		//		if (gnomeunique.size()< gnome.size()) {
		//			flag = 1;							//set flag for those if there are duplicated indices
		//												//std::cout<<"gnome:["<<g<<"] "<<"\t duplications are "<< (gnome.size() - gnomeunique.size())<<std::endl;
		//		}

		//		//adding extra random feature indices to unique copy of gnome to achive gnome size  = f
		//		if (gnomeunique.size() < f) {
		//			for (size_t k = gnomeunique.size(); k < f; k++) {
		//				unsigned int rnumber = rand() % ub + lb;
		//				//check if this randomaly generated number is already present in that gnome or not
		//				for (size_t j = 0; j < gnomeunique.size(); j++) {
		//					if (gnomeunique.at(j) == rnumber) {				//if new index is duplicated copy of any of previous gnome element replace it with another random number
		//						rnumber = rand() % ub + lb;
		//						j = 0;										//set j = 0 to start checking of duplication of feature index from the first element of gnome
		//					}
		//				}
		//				gnomeunique.push_back(rnumber);						//add feature index to gnomeunique 			
		//			}
		//		}
		//		std::copy(gnomeunique.begin(), gnomeunique.end(), P + g * f);
		//	}
		//}

		//checking bands with all zeros and replacing duplicated bands in gnome
		void zerobandcheck(float* M, bool initialPop) {
			size_t startgnome;
			if (initialPop) {
				startgnome = 0;    //for initial population check all gnomes 
			}
			else {
				startgnome = 1;    //for next generations start gnome check after elite children offset 
			}
			for (size_t g = startgnome; g < p; g++) {							// for each gnome except 		
				
				for (size_t i = 0; i < f; i++) {						//check each band (feature) index in that gnome
					while (M[P[g * f + i]] == 0) {						//if mean of band is zero then replace band index in population
						P[g * f + i] = rand() % ub + lb;
					}
				}
				//checking for duplicats in a gnome
				std::vector<unsigned int> gnome = getGnome(g);			//get current gnome g from population matrix P
				std::vector<unsigned int> gnomeunique;					//array to store only unique band indicies in a genome
				int flag = 0;											//flag will be set if gnome has duplicated band (feature) index
				std::sort(gnome.begin(), gnome.end());					//sort current gnome 
				std::unique_copy(gnome.begin(), gnome.end(), std::back_inserter(gnomeunique));		//remove duplicat copies of band indices and keep only unique in a gnome
				if (gnomeunique.size()< gnome.size()) {
					flag = 1;							//set flag for those if there are duplicated indices
					//std::cout<<"gnome:["<<g<<"] "<<"\t duplications are "<< (gnome.size() - gnomeunique.size())<<std::endl;
				}

				//adding extra random feature indices to unique copy of gnome to achive gnome size  = f
				if (gnomeunique.size() < f) {
					for (size_t k = gnomeunique.size(); k < f; k++) {
						unsigned int rnumber = rand() % ub + lb;
						//check if this randomaly generated number is already present in that gnome or not
						for (size_t j = 0; j < gnomeunique.size(); j++) {
							if (gnomeunique.at(j) == rnumber) {				//if new index is duplicated copy of any of previous gnome element replace it with another random number
								rnumber = rand() % ub + lb;					//generate random number between upper bound and lower bound (ub. lb)
								j = 0;										//set j = 0 to start checking of duplication of feature index from the first element of gnome
							}
						}
						gnomeunique.push_back(rnumber);						//add feature index to gnomeunique 			
					}
				}

				//diplay loop only if gnome has duplicated indices
				//if(flag ==1){
				//	std::cout<<"\n original gnome "<<g<<" are "<<std::endl;
				//	for(int k = 0; k < gnome.size(); k++)
				//		std::cout<<gnome[k]<<"\t";
				//	std::cout<<"\n unique results in cpp for gnome "<<g<<" are "<<std::endl;
				//	for(int k = 0; k < gnomeunique.size(); k++)
				//		std::cout<<gnomeunique[k]<<"\t";
				//}
				std::copy(gnomeunique.begin(), gnomeunique.end(), P + g * f);  //copy new gnome without any duplicate band index at current gnome location
			}
		}

		

		//gpu calling functions
		//gpu initialization (allocating space for all array on GPU)
		void gpuInitializationfrommain(float* cpuM, float* cpuCM, std::vector<unsigned int>cpu_nPxInCls, size_t tP, size_t nC) {
			//	call gpuInitialization(......) with all of the necessary parameters
			gpuIntialization(&gpuP, p, f, &gpuCM, cpuCM, nC, ub, &gpuM, cpuM, &gpu_nPxInCls, &gpuSb, &gpuSw, &gpuF, F, &gpuT, T, tP, &cpu_nPxInCls[0]);
			
		}

		//Computation of between class scatter and within class scatter in GPU
		void gpu_computeSbSw(float* cpuSb, float* cpuSw, size_t nC, size_t tP, cudaDeviceProp props, size_t gen, bool debug, std::ofstream& profilefile) {
			//calling function for SW and Sb computation and passing necessary arrays for computation
     // std::cout<<"gpu function calling"<<std::endl;
			gpucomputeSbSw(gpuP, P, p, f, gpuSb, cpuSb, gpuSw, cpuSw, gpuF, gpuT, gpuM, gpuCM, nC, tP, props, gen, gnrtn, ub, gpu_nPxInCls, profilefile);

			//display computed Sb and Sw if debug is set 
			if (debug) {
				std::cout << "From GA-GPU class: gpu results of Sb sn Sw" << std::endl;
				displayS(cpuSb, f);									//display Sb				
				displayS(cpuSw, f);									//display Sw
				std::cout << std::endl;
			}
		}

		//call function to free gpu pointers
		//free all gpu pointers
		void gpu_Destroy() {
			gpuDestroy(gpuP, gpuCM, gpuM, gpu_nPxInCls, gpuSb, gpuSw, gpuF, gpuT);
		}

		//Write a destructor here
		~ga_gpu() {

			if (F != NULL) 		std::free(F);		//not sure about this as it is only for 2nd constructor 
			if (T != NULL)		std::free(T);		//same as above
			if (P != NULL)		std::free(P);		//not sure about this as it is only for 2nd constructor 
			if (S != NULL)		std::free(S);		//same as above
														//if(i_guess!=NULL) 	std::free(i_guess);     //same as above
														//HANDLE_ERROR(cudaDeviceReset());

		}
	};

#endif