Commit 8967610271d60df76928a31790ca0bc61cb29a94
1 parent
34f743be
re-added main.cu
Showing
1 changed file
with
513 additions
and
0 deletions
Show diff stats
| 1 | +#include <iostream> | ||
| 2 | + | ||
| 3 | +//stim libraries | ||
| 4 | +#include <stim/envi/envi.h> | ||
| 5 | +#include <stim/image/image.h> | ||
| 6 | +#include <stim/ui/progressbar.h> | ||
| 7 | +#include <stim/parser/filename.h> | ||
| 8 | +#include <stim/parser/table.h> | ||
| 9 | +#include <stim/parser/arguments.h> | ||
| 10 | +//input arguments | ||
| 11 | +stim::arglist args; | ||
| 12 | +#include <fstream> | ||
| 13 | +#include <thread> | ||
| 14 | +#include <random> | ||
| 15 | +#include <vector> | ||
| 16 | +#include <math.h> | ||
| 17 | +#include <limits> | ||
| 18 | + | ||
| 19 | +#define NOMINMAX | ||
| 20 | + | ||
| 21 | + | ||
| 22 | + | ||
| 23 | +//GA | ||
| 24 | +#include "ga_gpu.h" | ||
| 25 | +#include "enviload.h" | ||
| 26 | + | ||
| 27 | + | ||
| 28 | +//envi input file and associated parameters | ||
| 29 | +stim::envi E; //ENVI binary file object | ||
| 30 | +unsigned int B; //shortcuts storing the spatial and spectral size of the ENVI image | ||
| 31 | +//mask and class information used for training | ||
| 32 | +//std::vector< stim::image<unsigned char> > C; //2D array used to access each mask C[m][p], where m = mask# and p = pixel# | ||
| 33 | +std::vector<unsigned int> nP; //array holds the number of pixels in each mask: nP[m] is the number of pixels in mask m | ||
| 34 | +size_t nC = 0; //number of classes | ||
| 35 | +size_t tP = 0; //total number of pixels in all masks: tP = nP[0] + nP[1] + ... + nP[nC] | ||
| 36 | +float* fea; | ||
| 37 | + | ||
| 38 | +//ga_gpu class object | ||
| 39 | +ga_gpu ga; | ||
| 40 | +bool debug; | ||
| 41 | +bool binaryClass; | ||
| 42 | +int binClassOne; | ||
| 43 | + | ||
| 44 | +//creating struct to pass to thread functions as it limits number of arguments to 3 | ||
| 45 | +typedef struct { | ||
| 46 | + float* S; | ||
| 47 | + float* Sb; | ||
| 48 | + float* Sw; | ||
| 49 | + float* lda; | ||
| 50 | +}gnome; | ||
| 51 | +gnome gnom; | ||
| 52 | + | ||
| 53 | + | ||
| 54 | +void gpuComputeEignS( size_t g, size_t fea){ | ||
| 55 | + //eigen value computation will return r = (nC-1) eigen vectors so new projected data will have dimension of r rather than f | ||
| 56 | + // std::thread::id this_id = std::this_thread::get_id(); | ||
| 57 | + // std::cout<<"thread id is "<< this_id<<std::endl; | ||
| 58 | + size_t f = fea; | ||
| 59 | + //std::thread::id g = std::this_thread::get_id(); | ||
| 60 | + float* LeftEigVectors_a = (float*) malloc(f * f * sizeof(float)); | ||
| 61 | + float* gSw_a = (float*) malloc(f * f * sizeof(float)); //copy of between class scatter | ||
| 62 | + std::memcpy(gSw_a, &gnom.Sw[g * f * f], f * f *sizeof(float)); | ||
| 63 | + if(debug){ | ||
| 64 | + std::cout<<"From Eigen function: Sb and Sw "<<std::endl; | ||
| 65 | + displayS(gSw_a, f); //display Sb | ||
| 66 | + displayS(&gnom.Sb[g * f * f], f); //display Sw | ||
| 67 | + std::cout<<std::endl; | ||
| 68 | + } | ||
| 69 | + | ||
| 70 | + std::vector<unsigned int> features = ga.getGnome(g); | ||
| 71 | + std::vector<unsigned int> featuresunique; | ||
| 72 | + int flag = 0; | ||
| 73 | + std::sort(features.begin(), features.end()); // 1 1 2 2 3 3 3 4 4 5 5 6 7 | ||
| 74 | + std::unique_copy(features.begin(), features.end(), std::back_inserter(featuresunique)); | ||
| 75 | + if(featuresunique.size()< features.size()){ | ||
| 76 | + f = featuresunique.size(); | ||
| 77 | + } | ||
| 78 | + | ||
| 79 | + size_t r = nC-1; //LDA projected dimension (limited to number of classes - 1 by rank) | ||
| 80 | + if(r > f){ | ||
| 81 | + r = f; | ||
| 82 | + } | ||
| 83 | + | ||
| 84 | + int info; | ||
| 85 | + float* EigenvaluesI_a = (float*)malloc(f * sizeof(float)); | ||
| 86 | + float* Eigenvalues_a = (float*)malloc(f * sizeof(float)); | ||
| 87 | + int *IPIV = (int*) malloc(sizeof(int) * f); | ||
| 88 | + //computing inverse of matrix Sw | ||
| 89 | + memset(IPIV, 0, f * sizeof(int)); | ||
| 90 | + LAPACKE_sgetrf(LAPACK_COL_MAJOR, (int)f, (int)f, gSw_a, (int)f, IPIV); | ||
| 91 | + // DGETRI computes the inverse of a matrix using the LU factorization computed by DGETRF. | ||
| 92 | + LAPACKE_sgetri(LAPACK_COL_MAJOR, (int)f, gSw_a, (int)f, IPIV); | ||
| 93 | + | ||
| 94 | + float* gSbSw_a = (float*)calloc(f * f, sizeof(float)); | ||
| 95 | + //mtxMul(gSbSw_a, gSw_a, &gnom.Sb[g * f * f * sizeof(float)], f, f, f,f); | ||
| 96 | + mtxMul(gSbSw_a, gSw_a, &gnom.Sb[g * f * f], f, f, f,f); | ||
| 97 | + if(debug){ | ||
| 98 | + std::cout<<"From Eigen function: inverse of sw and ratio of sb and sw (Sb/Sw)"; | ||
| 99 | + displayS(gSw_a, f); //display inverse of Sw (1/Sw) | ||
| 100 | + displayS(gSbSw_a, f); //display ratio of Sb and Sw (Sb/Sw) | ||
| 101 | + } | ||
| 102 | + | ||
| 103 | + //compute left eigenvectors for current gnome from ratio of between class scatter and within class scatter: Sb/Sw | ||
| 104 | + info = LAPACKE_sgeev(LAPACK_COL_MAJOR, 'V', 'N', (int)f, gSbSw_a, (int)f, Eigenvalues_a, EigenvaluesI_a, LeftEigVectors_a, (int)f, 0, (int)f); | ||
| 105 | + //sort eignevalue indices in descending order | ||
| 106 | + size_t* sortedindx = sortIndx(Eigenvalues_a, f); | ||
| 107 | + //displayS(LeftEigVectors_a, f); //display Eignevectors (Note these are -1 * matlab eigenvectors does not change fitness score results but keep in mind while projecting data on it) | ||
| 108 | + //sorting left eigenvectors (building forward transformation matrix As) | ||
| 109 | + for (size_t rowE = 0; rowE < r; rowE++){ | ||
| 110 | + for (size_t colE = 0; colE < f; colE++){ | ||
| 111 | + size_t ind1 = g * r * f + rowE * f + colE; | ||
| 112 | + //size_t ind1 = rowE * f + colE; | ||
| 113 | + size_t ind2 = sortedindx[rowE] * f + colE; //eigenvector as row vector | ||
| 114 | + gnom.lda[ind1] = LeftEigVectors_a[ind2]; | ||
| 115 | + } | ||
| 116 | + } | ||
| 117 | + | ||
| 118 | + if(debug){ | ||
| 119 | + std::cout<<"Eigenvalues are"<<std::endl; | ||
| 120 | + for(size_t n = 0 ; n < f; n ++){ | ||
| 121 | + std::cout << Eigenvalues_a[n] << ", " ; | ||
| 122 | + } | ||
| 123 | + std::cout<< std::endl; | ||
| 124 | + std::cout<<"From Eigen function: Eignevector"<<std::endl; | ||
| 125 | + | ||
| 126 | + std::cout<<"LDA basis is "<<std::endl; | ||
| 127 | + std::cout << "r is " << r << std::endl; | ||
| 128 | + for(size_t l = 0 ; l < r; l++){ | ||
| 129 | + for(size_t n = 0 ; n < f; n ++){ | ||
| 130 | + std::cout << gnom.lda[g * l * f + l * f + n] << ", " ; | ||
| 131 | + } | ||
| 132 | + std::cout<<std::endl; | ||
| 133 | + } | ||
| 134 | + | ||
| 135 | + } | ||
| 136 | + //Extract only r eigne vectors as a LDA projection basis | ||
| 137 | + float* tempgSb = (float*)calloc(r * f, sizeof(float)); | ||
| 138 | + //mtxMul(tempgSb, &gnom.lda[g * r * f * sizeof(float)], &gnom.Sb[g * f * f * sizeof(float)], r, f, f,f); | ||
| 139 | + //mtxMul(tempgSb, &lda[g * r * f ], gSb, r, f, f,f); | ||
| 140 | + mtxMul(tempgSb, &gnom.lda[g * r * f], &gnom.Sb[g * f * f], r, f, f,f); | ||
| 141 | + float* nSb = (float*)calloc(r * r, sizeof(float)); | ||
| 142 | + mtxMultranspose(nSb, tempgSb, &gnom.lda[g * r * f], r, f, r, f); | ||
| 143 | + | ||
| 144 | + float* tempgSw = (float*)calloc(r * f, sizeof(float)); | ||
| 145 | + //mtxMul(tempgSw, &gnom.lda[g * r * f * sizeof(float)], &gnom.Sw[g * f * f * sizeof(float)], r, f, f,f); | ||
| 146 | + mtxMul(tempgSw, &gnom.lda[g * r * f], &gnom.Sw[g * f * f], r, f, f,f); | ||
| 147 | + float* nSw = (float*)calloc(r * r, sizeof(float)); | ||
| 148 | + mtxMultranspose(nSw, tempgSw, &gnom.lda[g * r * f], r, f, r, f); | ||
| 149 | + if(debug){ | ||
| 150 | + std::cout<<"From Eigen function: projected Sb sn Sw"<<std::endl; | ||
| 151 | + displayS(nSb, r); //display Sb | ||
| 152 | + displayS(nSw, r); //display Sw | ||
| 153 | + std::cout<<std::endl; | ||
| 154 | + } | ||
| 155 | + | ||
| 156 | + ///DETERMINANT CURRENTLY REQUIRES OPENCV | ||
| 157 | + std::cout<<"ERROR: This code requires a fix to remove an OpenCV dependence."<<std::endl; | ||
| 158 | + mtxOutputFile("newSw.csv", nSw, r, r); | ||
| 159 | + exit(1); | ||
| 160 | + //FIX BY REPLACING THE FOLLOWING THREE LINES OF CODE USING LAPACK | ||
| 161 | + //cv::Mat newSw = cv::Mat((int)r, (int)r, CV_32FC1, nSw); //within scatter of gnome g in the population | ||
| 162 | + //cv::Mat newSb = cv::Mat((int)r, (int)r, CV_32FC1, nSb); //within scatter of gnome g in the population | ||
| 163 | + | ||
| 164 | + //fisher's ratio from ratio of projected sb and sw | ||
| 165 | + float fisherRatio = 0;// = cv::determinant(newSb) /cv::determinant(newSw); | ||
| 166 | + gnom.S[g] = 1/fisherRatio; | ||
| 167 | + if (debug) { | ||
| 168 | + std::cout<<"Score["<<g<<"]: "<< gnom.S[g]<<std::endl; | ||
| 169 | + | ||
| 170 | + std::cout << "best gnoem is " << std::endl; | ||
| 171 | + for (size_t i = 0; i < f; i++) | ||
| 172 | + std::cout << ga.P[ga.f * g + i] << ", "; | ||
| 173 | + std::cout << std::endl; | ||
| 174 | + } | ||
| 175 | + | ||
| 176 | + | ||
| 177 | + if(IPIV!= NULL) std::free(IPIV); | ||
| 178 | + if(gSw_a!= NULL) std::free(gSw_a); | ||
| 179 | + if(gSbSw_a!= NULL) std::free(gSbSw_a); | ||
| 180 | + if(Eigenvalues_a!= NULL) std::free(Eigenvalues_a); | ||
| 181 | + if(EigenvaluesI_a!= NULL) std::free(EigenvaluesI_a); | ||
| 182 | + if(tempgSb!= NULL) std::free(tempgSb); | ||
| 183 | + if(tempgSw!= NULL) std::free(tempgSw); | ||
| 184 | + | ||
| 185 | +} | ||
| 186 | + | ||
| 187 | + | ||
| 188 | +void fitnessFunction( float* sb, float* sw, float* lda, float* M, float* cM, size_t f, cudaDeviceProp props, size_t gen, std::ofstream& profilefile){ | ||
| 189 | + | ||
| 190 | + size_t tP = 0; //total number of pixels | ||
| 191 | + std::for_each(nP.begin(), nP.end(), [&] (size_t n) { | ||
| 192 | + tP += n; | ||
| 193 | + }); | ||
| 194 | + size_t nC = nP.size(); //total number of classes | ||
| 195 | + | ||
| 196 | + //--------------Compute between class scatter | ||
| 197 | + | ||
| 198 | + ga.gpu_computeSbSw(sb, sw, nC, tP, props, gen, debug, profilefile); | ||
| 199 | + | ||
| 200 | + if(debug){ | ||
| 201 | + std::cout<<"From fitness function: gpu results of Sb sn Sw"<<std::endl; | ||
| 202 | + displayS(sb, ga.f); //display Sb | ||
| 203 | + displayS(sw, ga.f); //display Sw | ||
| 204 | + std::cout<<std::endl; | ||
| 205 | + } | ||
| 206 | + | ||
| 207 | + // ----------------------- Linear discriminant Analysis -------------------------------------- | ||
| 208 | + gnom.S = ga.S; | ||
| 209 | + gnom.Sw = sw; | ||
| 210 | + gnom.Sb = sb; | ||
| 211 | + gnom.lda = lda; | ||
| 212 | + | ||
| 213 | + //calling function without using threads | ||
| 214 | + for (size_t i = 0; i<ga.p; i++){ | ||
| 215 | + //calling function for eigencomputation | ||
| 216 | + gpuComputeEignS(i, f); | ||
| 217 | + } | ||
| 218 | + | ||
| 219 | + const auto elapsed1 = timer.time_elapsed(); | ||
| 220 | + if(gen > ga.gnrtn - 2){ | ||
| 221 | + std::cout << "gpu_eigen time "<<std::chrono::duration_cast<std::chrono::microseconds>(elapsed1).count() << "us" << std::endl; | ||
| 222 | + profilefile<< "gpu_eigen time "<<std::chrono::duration_cast<std::chrono::microseconds>(elapsed1).count() << "us" << std::endl; | ||
| 223 | + } | ||
| 224 | + | ||
| 225 | +}//end of fitness function | ||
| 226 | + | ||
| 227 | +void binaryclassifier(int classnum){ | ||
| 228 | + unsigned int* target = (unsigned int*) calloc(tP, sizeof(unsigned int)); | ||
| 229 | + memcpy(target, ga.T, tP * sizeof(unsigned int)); | ||
| 230 | + for(int i = 0 ; i < tP; i++){ | ||
| 231 | + if(target[i]==classnum){ | ||
| 232 | + ga.T[i] = 1; | ||
| 233 | + | ||
| 234 | + }else | ||
| 235 | + ga.T[i] = 0; | ||
| 236 | + } | ||
| 237 | +} | ||
| 238 | + | ||
| 239 | + | ||
| 240 | + | ||
| 241 | +void advertisement() { | ||
| 242 | + std::cout << std::endl; | ||
| 243 | + std::cout << "=========================================================================" << std::endl; | ||
| 244 | + std::cout << "Thank you for using the GA-GPU features selection for spectroscopic image!" << std::endl; | ||
| 245 | + std::cout << "=========================================================================" << std::endl << std::endl; | ||
| 246 | +} | ||
| 247 | + | ||
| 248 | +int main(int argc, char* argv[]){ | ||
| 249 | + | ||
| 250 | +//Add the argument options and set some of the default parameters | ||
| 251 | + args.add("help", "print this help"); | ||
| 252 | + args.section("Genetic Algorithm"); | ||
| 253 | + args.add("features", "select features selection algorithm parameters","10", "number of features to be selected"); | ||
| 254 | + args.add("classes", "image masks used to specify classes", "", "class1.bmp class2.bmp class3.bmp"); | ||
| 255 | + args.add("population", "total number of feature subsets in puplation matrix", "1000"); | ||
| 256 | + args.add("generations", "number of generationsr", "50"); | ||
| 257 | + args.add("initial_guess", "initial guess of featues", ""); | ||
| 258 | + args.add("debug", "display intermediate data for debugging"); | ||
| 259 | + args.add("binary", "Calculate features based on class1 vs. all other classes", ""); | ||
| 260 | + args.add("trim", "this gives wavenumber to use in trim option of siproc which trims all bands from envi file except gagpu selected bands"); | ||
| 261 | + | ||
| 262 | + args.parse(argc,argv); //parse the command line arguments | ||
| 263 | + | ||
| 264 | +//Print the help text if set | ||
| 265 | + if(args["help"].is_set()){ //display the help text if requested | ||
| 266 | + advertisement(); | ||
| 267 | + std::cout<<std::endl<<"usage: ga-gpu input_ENVI output.txt --classes class1.bmp class2.bmp ... --option [A B C ...]"<<std::endl; | ||
| 268 | + std::cout<<std::endl<<std::endl; | ||
| 269 | + std::cout<<args.str()<<std::endl; | ||
| 270 | + exit(1); | ||
| 271 | + } | ||
| 272 | + if (args.nargs() < 2) { //if the user doesn't provide input and output files | ||
| 273 | + std::cout << "ERROR: GA-GPU requires an input (ENVI) file and an output (features, text) file." << std::endl; | ||
| 274 | + return 1; | ||
| 275 | + } | ||
| 276 | + if (args["classes"].nargs() < 2) { //if the user doesn't specify at least two class images | ||
| 277 | + std::cout << "ERROR: GA-GPU requires at least two class images to be specified using the --classes option" << std::endl; | ||
| 278 | + return 1; | ||
| 279 | + } | ||
| 280 | + | ||
| 281 | + std::string outfile = args.arg(1); //outfile is text file where bnad index, LDA-basis, wavelength and if --trim option is set then trim wavelengths are set respectively | ||
| 282 | + std::string profile_file = "profile_" + outfile ; | ||
| 283 | + std::ofstream profilefile(profile_file.c_str(), std::ios::out); //open outfstream for outfile | ||
| 284 | + | ||
| 285 | + time_t t_start = time(NULL); //start a timer for file reading | ||
| 286 | + E.open(args.arg(0), std::string(args.arg(0)) + ".hdr"); //open header file | ||
| 287 | + size_t X = E.header.samples; //total number of pixels in X dimension | ||
| 288 | + size_t Y = E.header.lines; //total number of pixels in Y dimension | ||
| 289 | + B = (unsigned int)E.header.bands; //total number of bands (features) | ||
| 290 | + std::vector<double> wavelengths = E.header.wavelength; //wavelengths of each band | ||
| 291 | + | ||
| 292 | + if(E.header.interleave != stim::envi_header::BIP){ //this code can only load bip files and hence check that in header file | ||
| 293 | + std::cout<<"this code works for only bip files. please convert file to bip file"<<std::endl; | ||
| 294 | + exit(1); //if file is not bip file exit code execution | ||
| 295 | + } | ||
| 296 | + | ||
| 297 | +///--------------------------Load features--------------------------------------------- | ||
| 298 | + nP = ga_load_class_images(argc, args, &nC, &tP); //load supervised class images | ||
| 299 | + ga.F = load_features( nC, tP, B, E, nP); //generate the feature matrix | ||
| 300 | + ga.T = ga_load_responses(tP, nC, nP); //load the responses for RF training | ||
| 301 | + E.close(); //close the hyperspectral file | ||
| 302 | + time_t t_end = time(NULL); | ||
| 303 | + std::cout<<"Total time: "<<t_end - t_start<<" s"<<std::endl; | ||
| 304 | + | ||
| 305 | +///--------------------------Genetic algorith configurations with defult paramets and from argument values--------------------- | ||
| 306 | + ga.f = args["features"].as_int(0); //number of features to be selected by user default value is 10 | ||
| 307 | + ga.p = args["population"].as_int(0); //population size to be selected by user default value is 1000 | ||
| 308 | + ga.gnrtn = args["generations"].as_int(0); //number of generations to be selected by user default value is 50 | ||
| 309 | + if(args["binary"]) { //set this option when features are to be selected as binary clas features (class vs stroma) | ||
| 310 | + binClassOne = args["binary"].as_int(0); //sel class number here, if 2 then features are selected for (class-2 vs stroma) | ||
| 311 | + //feture selection for class selected by user with user arguments (make it binary class data by making chosen class label as 1 and al other class labels 0 from multiclass data ) | ||
| 312 | + //to select feature for all classes in joint class data using binary class system need to write a script with loop covering all classes | ||
| 313 | + binaryclassifier(binClassOne); | ||
| 314 | + } ///not fully implemented yet | ||
| 315 | + | ||
| 316 | + ga.ub = B; //upper bound is number of bands (i.e. size of z dimension) Note: for this particular application and way code is written lower bound is 0 and upper bound is size of z dimension | ||
| 317 | + ga.uniformRate = 0.5; //uniform rate is used in crossover | ||
| 318 | + ga.mutationRate = 0.5f; //in percentage for mutation operation on gnome | ||
| 319 | + ga.tournamentSize = 5; //for crossover best parents are selected from tournament of gnomes | ||
| 320 | + ga.elitism = true; // if it is true then best gnome of current generation is passed to next generation | ||
| 321 | + //initial guess of population | ||
| 322 | + ga.i_guess = (unsigned int*) calloc(ga.f, sizeof(unsigned int)); | ||
| 323 | + debug = args["debug"]; | ||
| 324 | + | ||
| 325 | +//==================Generate intial population ================================================= | ||
| 326 | + std::vector<unsigned int> i_guess(ga.f); | ||
| 327 | + for (size_t b = 0; b < ga.f; b++) //generate default initial guess | ||
| 328 | + i_guess[b] = rand() % B + 0; | ||
| 329 | + | ||
| 330 | + if (args["initial_guess"].is_set()) {//if the user specifies the --initialguess option & provides feature indices as initial guess | ||
| 331 | + size_t nf = args["initial_guess"].nargs(); //get the number of arguments after initial_guess | ||
| 332 | + if (nf == 1 || nf == ga.f) { //check if file with initial guessed indices is given or direct indices are given as argument | ||
| 333 | + if (nf == 1) { //if initial guessed feature indices are given in file | ||
| 334 | + std::ifstream in; //open the file containing the baseline points | ||
| 335 | + in.open(args["initial_guess"].as_string().c_str()); | ||
| 336 | + if (in.is_open()){ //if file is present and can be opened then read it | ||
| 337 | + unsigned int b_ind; | ||
| 338 | + while (in >> b_ind) //get each band index and push it into the vector | ||
| 339 | + i_guess.push_back(b_ind); | ||
| 340 | + } | ||
| 341 | + else | ||
| 342 | + std::cout << "cannot open file of initial_guess indices" << std::endl; | ||
| 343 | + } | ||
| 344 | + else if (nf == ga.f) { //if direct indices are given as argument | ||
| 345 | + for (size_t b = 0; b < nf; b++) //for each band given by the user | ||
| 346 | + i_guess[b] = args["initial_guess"].as_int(b); //store that band in the i_guess array | ||
| 347 | + } | ||
| 348 | + } | ||
| 349 | + } | ||
| 350 | + | ||
| 351 | + ga.initialize_population(i_guess, debug); //initialize first population set for first generation, user can pass manually preferred features from command line | ||
| 352 | + //display_gnome(0); | ||
| 353 | + | ||
| 354 | +//------------------Calculate class means and total mean of features---------------------------- | ||
| 355 | + float* M = (float*)calloc( B , sizeof(float)); //total mean of entire feature martix for all features (bands B) | ||
| 356 | + ga.ttlMean(M, tP, B); //calculate total mean, ga.F is entire feature matrix, M is mean for all bands B(features) | ||
| 357 | + if(debug) ga.dispalymean(M); //if option --debug is used display all bands mean | ||
| 358 | + | ||
| 359 | + //display band index of bands with mean zero, this indicates that band has all zero values | ||
| 360 | + std::cout<<"Display features indices with zero mean "<<std::endl; | ||
| 361 | + for(unsigned int i = 0; i < B; i++){ | ||
| 362 | + if(M[i]== 0) | ||
| 363 | + std::cout<<"\t"<<i; | ||
| 364 | + } | ||
| 365 | + std::cout<<std::endl; | ||
| 366 | +// std::cout << "pixel target is " << ga.T[0] << " " << ga.T[1] << " " << ga.T[tP - 2] << " " << ga.T[tP - 1]<<std::endl; | ||
| 367 | + float* cM = (float*)calloc(nC * B , sizeof(float)); //cM is nC X B matrix with each row as mean of all samples in one class for all features (bands B) | ||
| 368 | + ga.classMean(cM, tP, nC, B, nP); //calculate class mean, ga.F is entire feature matrix, M is mean for all bands B(features) | ||
| 369 | + if(debug) ga.dispalyClassmean(cM, nC); | ||
| 370 | + | ||
| 371 | +//------------------------------------GPU init---------------------------------------------------- | ||
| 372 | + //checking for cuda device | ||
| 373 | + int count; | ||
| 374 | + HANDLE_ERROR(cudaGetDeviceCount(&count)); | ||
| 375 | + if(count < 1){ | ||
| 376 | + std::cout<<"no cuda device is available"<<std::endl; | ||
| 377 | + return 1; | ||
| 378 | + } | ||
| 379 | + cudaDeviceProp props; | ||
| 380 | + HANDLE_ERROR(cudaGetDeviceProperties(&props, 0)); | ||
| 381 | + ga.gpuInitializationfrommain(M, cM, nP, tP, nC); | ||
| 382 | + | ||
| 383 | + | ||
| 384 | +//============================= GA evolution by generations ==================================================== | ||
| 385 | + std::vector<unsigned int> bestgnome; //holds best gnome after each generation evaluation | ||
| 386 | + size_t bestG_Indx; //This gives index of best gnome in the current population to get best gnome and its fitness value | ||
| 387 | + unsigned int* newPop = (unsigned int*) calloc(ga.p * ga.f, sizeof(unsigned int)); //temprory storage of new population | ||
| 388 | + double* best_S = (double*) calloc (ga.gnrtn, sizeof(double)); //stores fitness value of best gnome at each iteration | ||
| 389 | + float* lda = (float*) calloc (ga.p * (nC-1) * ga.f, sizeof(float)); //stores LDA basis for each gnome so that we can have best gnome's LDA basis | ||
| 390 | + float* sb = (float*) calloc( ga.p * ga.f * ga.f , sizeof(float)) ; //3d matrix for between class scatter (each 2d matrix between class scatter for one gnome) | ||
| 391 | + float* sw = (float*) calloc( ga.p * ga.f * ga.f , sizeof(float)) ; //3d matrix for within class scatter (each 2d matrix within class scatter for one gnome) | ||
| 392 | + ga.zerobandcheck(M, true); //checking bands with all zeros and duplicated bands in a gnome replacing them with other bands avoiding duplication and zero mean | ||
| 393 | + ga.zerobandcheck(M, true); //Repeating zeroband cheack as some of these bands are not replaced in previous run and gave random results | ||
| 394 | + time_t gpu_timestart = time(NULL); //start a timer for total evoluation | ||
| 395 | + | ||
| 396 | + for (size_t gen = 0; gen < ga.gnrtn; gen++){ //for each generation find fitness value of all gnomes in population matrix and generate population for next generation | ||
| 397 | + //std::cout<<"Generation: "<<gen<<std::endl; | ||
| 398 | + fitnessFunction(sb, sw, lda, M , cM, ga.f, props, gen, profilefile); //Evaluate phe(feature matrix for current population) for fitness of all gnomes in current population | ||
| 399 | + timer.start(); //start timer for new population generation | ||
| 400 | + bestG_Indx = ga.evolvePopulation(newPop, M, debug); //evolve population to generate new generation population | ||
| 401 | + const auto pop_generation = timer.time_elapsed(); // end timer for new population generation | ||
| 402 | + if(gen >ga.gnrtn -2){ | ||
| 403 | + std::cout << "population evolution time "<<std::chrono::duration_cast<std::chrono::microseconds>(pop_generation).count() << "us" << std::endl; | ||
| 404 | + profilefile<<"population evolution time "<<std::chrono::duration_cast<std::chrono::microseconds>(pop_generation).count() << "us" <<std::endl; | ||
| 405 | + } | ||
| 406 | + | ||
| 407 | + best_S[gen] = ga.S[bestG_Indx]; //score of best gnome in current generation | ||
| 408 | + bestgnome = ga.getGnome(bestG_Indx); //Best gnome of current populaation | ||
| 409 | + ga.generateNewP(newPop); //replace current population with new populaiton in the ga classs object | ||
| 410 | + ga.zerobandcheck(M, false); //checking bands with all zeros and duplicated bands in a gnome replacing them with other bands avoiding duplication and zero mean | ||
| 411 | + ga.zerobandcheck(M, false); //Repeating zeroband cheack as some of these bands are not replaced in previous run and gave random results | ||
| 412 | + }//end generation | ||
| 413 | + | ||
| 414 | + time_t gpu_timeend = time(NULL); //end a timer for total evoluation | ||
| 415 | + std::cout<<"Total gpu time: "<<gpu_timeend - gpu_timestart<<" s"<<std::endl; | ||
| 416 | + profilefile<<"Total gpu time: "<<gpu_timeend - gpu_timestart<<" s"<<std::endl; | ||
| 417 | + | ||
| 418 | +//================================ Results of GA =============================================================== | ||
| 419 | + std::cout<<"best gnome's fitness value is "<<best_S[ga.gnrtn-1]<<std::endl; | ||
| 420 | + std::cout<<"best gnome is: "; | ||
| 421 | + for(size_t i = 0; i < ga.f; i++){ | ||
| 422 | + std::cout<<" "<<(bestgnome.at(i)); | ||
| 423 | + } | ||
| 424 | + std::cout<<std::endl; | ||
| 425 | + | ||
| 426 | + //create a text file to store the LDA stats (features subset and LDA-basis) | ||
| 427 | + ////format of CSV file is: 1st row - band index, 2nd LDA basis depending on number of classes, 3rd - wavenumber corresponding to band index and it --trim is selected then trim wavnumbersare also given | ||
| 428 | + std::ofstream csv(outfile.c_str(), std::ios::out); //open outfstream for outfile | ||
| 429 | + size_t ldaindx = bestG_Indx * (nC-1) * ga.f ; //Compute LDA basis index of best gnome | ||
| 430 | + | ||
| 431 | + | ||
| 432 | + //fitness values of best gnome is | ||
| 433 | + csv<<"best gnome's fitness value is "<<best_S[ga.gnrtn-1]<<std::endl; //output fitness value of best gnome in last generation | ||
| 434 | + //output gnome i.e. band index of selected featurs | ||
| 435 | + csv<<(bestgnome.at(0)); //output feature subset | ||
| 436 | + for(size_t i = 1; i < ga.f; i++) | ||
| 437 | + csv<<","<<(bestgnome.at(i)); | ||
| 438 | + csv<<std::endl; | ||
| 439 | + | ||
| 440 | + //output LDA basis of size r X f, r is nC - 1 as LDA projection is rank limited by number of classes - 1 | ||
| 441 | + for (size_t i = 0; i < nC-1; i++){ | ||
| 442 | + csv<<lda[ldaindx + i * ga.f ]; | ||
| 443 | + for (size_t j = 1; j < ga.f; j++){ | ||
| 444 | + csv<<","<<lda[ldaindx + i * ga.f +j]; | ||
| 445 | + } | ||
| 446 | + csv << std::endl; | ||
| 447 | + } | ||
| 448 | + //output actual wavelenths corresponding to those band indices | ||
| 449 | + csv << (wavelengths[bestgnome.at(0)]); | ||
| 450 | + for (size_t i = 1; i < ga.f; i++) | ||
| 451 | + csv << "," << (wavelengths[bestgnome.at(i)]); | ||
| 452 | + csv << std::endl; | ||
| 453 | + | ||
| 454 | + | ||
| 455 | + if (args["trim"].is_set()) { | ||
| 456 | + csv << "trim info" << std::endl; | ||
| 457 | + std::sort(bestgnome.begin(), bestgnome.end()); //sort features index in best gnome | ||
| 458 | + | ||
| 459 | + std::vector<unsigned int> trimindex(ga.f); //create a vector to store temprory trim index bounds | ||
| 460 | + std::vector<unsigned int> finaltrim_ind; //create a vector to store final trim index bounds | ||
| 461 | + std::vector<unsigned int> trim_wv; //create a vector to store final trim wavelength bounds | ||
| 462 | + | ||
| 463 | + //trim index | ||
| 464 | + trimindex.push_back(1); //1st trimming band index is 1 | ||
| 465 | + for (size_t i = 0; i < ga.f; i++) { // for each feature find its bound indexes | ||
| 466 | + trimindex[i * 2] = bestgnome.at(i) - 1; | ||
| 467 | + trimindex[i * 2 + 1] = bestgnome.at(i) + 1; | ||
| 468 | + } | ||
| 469 | + trimindex.push_back(B); //last bound index is B | ||
| 470 | + | ||
| 471 | + //organize trim index | ||
| 472 | + int k = 0; | ||
| 473 | + for (size_t i = 0; i < ga.f + 1; i++) { // find valid pair of trim indices bound excluding adjacent trim indices | ||
| 474 | + if (trimindex[2 * i] < trimindex[2 * i + 1]) { | ||
| 475 | + finaltrim_ind.push_back(trimindex[2 * i]); //this is left bound | ||
| 476 | + finaltrim_ind.push_back(trimindex[2 * i + 1]); | ||
| 477 | + k = k + 2; | ||
| 478 | + } | ||
| 479 | + } | ||
| 480 | + //add duplicated trim indices as single index to final trim index | ||
| 481 | + for (size_t i = 0; i < ga.f + 1; i++) { //check each pair of trim indices for duplications | ||
| 482 | + if (trimindex[2 * i] == trimindex[2 * i + 1]) { // (duplication caused due to adjacent features) | ||
| 483 | + finaltrim_ind.push_back(trimindex[2 * i]); // remove duplicated trim indices replace by one | ||
| 484 | + k = k + 1; | ||
| 485 | + } | ||
| 486 | + } | ||
| 487 | + | ||
| 488 | + | ||
| 489 | + ////output actual wavelenths corresponding to those trim indices | ||
| 490 | + ////these wavenumber are grouped in pairs, check each pair, if duplicated numbers are there in pair delete one and keep other as band to trim, if 2nd wavnumber is smaller than 1st in a pair ignore that pair | ||
| 491 | + ////e.g [1, 228, 230, 230, 232, 350,352, 351, 353, 1200] pairas [1(start)-228,230-230, 232-350, 352-351, 353-1200(end)], trimming wavenumbers are [1-228, 230, 233-350, 353-1200] | ||
| 492 | + csv << (wavelengths[finaltrim_ind.at(0)]); | ||
| 493 | + for (size_t i = 1; i < ga.f * + 2 ; i++) | ||
| 494 | + csv << "," << (wavelengths[finaltrim_ind.at(i)]); | ||
| 495 | + csv << std::endl; | ||
| 496 | + } //end trim option | ||
| 497 | + | ||
| 498 | + | ||
| 499 | + //free gpu pointers | ||
| 500 | + ga.gpu_Destroy(); | ||
| 501 | + | ||
| 502 | + //free pointers | ||
| 503 | + std::free(sb); | ||
| 504 | + std::free(sw); | ||
| 505 | + std::free(M); | ||
| 506 | + std::free(cM); | ||
| 507 | + std::free(best_S); | ||
| 508 | + std::free(lda); | ||
| 509 | + std::free(newPop); | ||
| 510 | + | ||
| 511 | +}//end main | ||
| 512 | + | ||
| 513 | + |