SedSat3 1.1.6
Sediment Source Apportionment Tool - Advanced statistical methods for environmental pollution research
Loading...
Searching...
No Matches
GA.hpp
Go to the documentation of this file.
1
2
3// GA.cpp: implementation of the CGA class.
5#include "GA.h"
6#include <stdlib.h>
7#ifndef mac_version
8#include <omp.h>
9#endif
10
11#include "parameter.h"
12
13#ifdef Q_version
14 #include "ProgressWindow.h"
15#endif
16
17template<class T>
18CGA<T>::CGA()
19{
20 GA_params.maxpop = 100;
21 Ind.resize(GA_params.maxpop);
22 Ind_old.resize(GA_params.maxpop);
23 fitdist = GADistribution(GA_params.maxpop);
24 GA_params.N = 1; //MM N=1 by default
25 GA_params.pcross = 1;
26 GA_params.cross_over_type = 1;
27 MaxFitness = 0;
28 numberOfThreads = 16;
29}
30
31template<class T>
32CGA<T>::CGA(string filename, const T &model)
33{
34 Model = model;
35 ifstream file(filename);
36 GA_params.nParam = 0;
37 GA_params.pcross = 1;
38 GA_params.N = 1;
39 GA_params.RCGA = false;
40 loged.clear();
41 numberOfThreads = 20;
42 filenames.pathname = Model->OutputPath();
43 vector<string> s;
44 while (file.eof() == false)
45 {
46 s = aquiutils::getline(file);
47 if (s.size()>0)
48 { if (s[0] == "maxpop") GA_params.maxpop = atoi(s[1].c_str());
49 if (s[0] == "ngen") GA_params.nGen = atoi(s[1].c_str());
50 if (s[0] == "pcross") GA_params.pcross = atof(s[1].c_str());
51 if (s[0] == "pmute") GA_params.pmute = atof(s[1].c_str());
52 if (s[0] == "shakescale") GA_params.shakescale = atof(s[1].c_str());
53 if (s[0] == "shakescalered") GA_params.shakescalered = atof(s[1].c_str());
54 if (s[0] == "outputfile") filenames.outputfilename = s[1];
55 if (s[0] == "getfromfilename") filenames.getfromfilename = s[1].c_str();
56 if (s[0] == "initial_population") filenames.initialpopfilemame = s[1];
57 if (s[0] == "numthreads") numberOfThreads = atoi(s[1].c_str());
58 }
59 }
60
61 file.close();
62 for (int i=0; i<Model->Parameters().size(); i++)
63 {
64 GA_params.nParam++;
65 params.push_back(i);
66 if (Model->Parameters()[i]->GetPriorDistribution() == "lognormal")
67 { minval.push_back(log10(Model->Parameters()[i]->GetRange(_range::low)));
68 maxval.push_back(log10(Model->Parameters()[i]->GetRange(_range::high)));
69
70 }
71 else
72 {
73 minval.push_back(Model->Parameters()[i]->GetRange(_range::low));
74 maxval.push_back(Model->Parameters()[i]->GetRange(_range::high));
75 }
76 apply_to_all.push_back(false);
77 if (Model->Parameters()[i]->GetPriorDistribution() == "lognormal")
78 loged.push_back(1);
79 else
80 loged.push_back(0);
81
82 paramname.push_back(Model->Parameters().getKeyAtIndex(i));
83
84 }
85
86
87 Ind.resize(GA_params.maxpop);
88 Ind_old.resize(GA_params.maxpop);
89
90 fitdist = GADistribution(GA_params.maxpop);
91 GA_params.cross_over_type = 1;
92
93 for (int i=0; i<GA_params.maxpop; i++)
94 {
95 Ind[i] = CIndividual(GA_params.nParam);
96 Ind[i].fit_measures.resize(model->ObservationsCount()*3);
97 Ind_old[i] = CIndividual(GA_params.nParam);
98 Ind_old[i].fit_measures.resize(model->ObservationsCount()*3);
99 }
100
101 for (int i = 0; i<GA_params.nParam; i++)
102 Setminmax(i, minval[i], maxval[i],4);
103
104 MaxFitness = 0;
105}
106
107template<class T>
108CGA<T>::CGA(T *model)
109{
110 Model = model;
111 GA_params.nParam = 0;
112 GA_params.pcross = 1;
113 GA_params.N = 1;
114 GA_params.RCGA = false;
115 numberOfThreads = 20;
116 loged.clear();
117 filenames.pathname = Model->GetOutputPath();
118 GA_params.maxpop = max(1,GA_params.maxpop);
119 for (unsigned int i=0; i<Model->Parameters().size(); i++)
120 {
121 GA_params.nParam++;
122 params.push_back(i);
123 if (Model->parameter(i)->GetPriorDistribution().distribution == distribution_type::lognormal)
124 { minval.push_back(log10(Model->parameter(i)->GetVal("low")));
125 maxval.push_back(log10(Model->parameter(i)->GetVal("high")));
126
127 }
128 else
129 {
130 minval.push_back(Model->parameter(i)->GetVal("low"));
131 maxval.push_back(Model->parameter(i)->GetVal("high"));
132 }
133 apply_to_all.push_back(false);
134 if (Model->parameter(i)->GetPriorDistribution().distribution == distribution_type::lognormal)
135 loged.push_back(1);
136 else
137 loged.push_back(0);
138
139 paramname.push_back(Model->GetParameterName(i));
140
141 }
142
143
144 Ind.resize(GA_params.maxpop);
145 Ind_old.resize(GA_params.maxpop);
146
147 fitdist = GADistribution(GA_params.maxpop);
148 GA_params.cross_over_type = 1;
149
150 for (int i=0; i<GA_params.maxpop; i++)
151 {
152 Ind[i] = CIndividual(GA_params.nParam);
153 Ind[i].fit_measures.resize(model->ObservationsCount());
154 Ind_old[i] = CIndividual(GA_params.nParam);
155 Ind_old[i].fit_measures.resize(model->ObservationsCount());
156 }
157
158 for (int j = 0; j < GA_params.nParam; j++)
159 Setminmax(j, minval[j], maxval[j], 4);
160
161 MaxFitness = 0;
162}
163
164template<class T>
165void CGA<T>::InitiatePopulation()
166{
167 GA_params.nParam = Model->Parameters().size();
168 loged.clear();
169 for (unsigned int i=0; i<Model->Parameters().size(); i++)
170 {
171 params.push_back(i);
172 if (Model->parameter(i)->GetPriorDistribution().distribution == distribution_type::lognormal)
173 { minval.push_back(log10(Model->parameter(i)->GetVal("low")));
174 maxval.push_back(log10(Model->parameter(i)->GetVal("high")));
175
176 }
177 else
178 {
179 minval.push_back(Model->parameter(i)->GetVal("low"));
180 maxval.push_back(Model->parameter(i)->GetVal("high"));
181 }
182 apply_to_all.push_back(false);
183 if (Model->parameter(i)->GetPriorDistribution().distribution == distribution_type::lognormal)
184 loged.push_back(1);
185 else
186 loged.push_back(0);
187
188 paramname.push_back(Model->GetParameterName(i));
189
190 }
191
192
193 Ind.resize(GA_params.maxpop);
194 Ind_old.resize(GA_params.maxpop);
195
196 fitdist = GADistribution(GA_params.maxpop);
197 GA_params.cross_over_type = 1;
198
199 for (int i=0; i<GA_params.maxpop; i++)
200 {
201 Ind[i] = CIndividual(GA_params.nParam);
202 Ind[i].fit_measures.resize(Model->ObservationsCount());
203 Ind_old[i] = CIndividual(GA_params.nParam);
204 Ind_old[i].fit_measures.resize(Model->ObservationsCount());
205 }
206
207 for (int j = 0; j < GA_params.nParam; j++)
208 Setminmax(j, minval[j], maxval[j], 4);
209
210 MaxFitness = 0;
211}
212
213template<class T>
214void CGA<T>::setnparams(int n_params)
215{
216 Ind.resize(GA_params.maxpop);
217 Ind_old.resize(GA_params.maxpop);
218 for (int i=0; i<GA_params.maxpop; i++)
219 {
220 Ind[i] = CIndividual(n_params);
221 Ind_old[i] = CIndividual(n_params);
222
223 }
224
225
226}
227
228template<class T>
229void CGA<T>::setnumpop(int n)
230{
231 GA_params.maxpop = n;
232 CIndividual TempInd = Ind[0];
233
234 int nParam = Ind[0].nParams;
235 Ind.resize(GA_params.maxpop);
236 Ind_old.resize(GA_params.maxpop);
237 for (int i=0; i<n; i++)
238 {
239 Ind[i] = CIndividual(GA_params.nParam);
240 Ind_old[i] = CIndividual(GA_params.nParam);
241 for (int j = 0; j<nParam; j++)
242 {
243 Ind[i].minrange[j] = TempInd.minrange[j];
244 Ind[i].maxrange[j] = TempInd.maxrange[j];
245 Ind[i].precision[j] = TempInd.precision[j];
246 Ind_old[i].minrange[j] = TempInd.minrange[j];
247 Ind_old[i].maxrange[j] = TempInd.maxrange[j];
248 Ind_old[i].precision[j] = TempInd.precision[j];
249 }
250
251 }
252 fitdist = GADistribution(GA_params.maxpop);
253}
254
255template<class T>
256CGA<T>::CGA(const CGA<T> &C)
257{
258 GA_params.maxpop = C.GA_params.maxpop;
259 Ind.resize(GA_params.maxpop);
260 Ind_old.resize(GA_params.maxpop);
261 Ind = C.Ind;
262 GA_params = C.GA_params;
263 filenames = C.filenames;
264 params = C.params;
265 loged = C.loged;
266 fitdist = C.fitdist;
267 MaxFitness = C.MaxFitness;
268 paramname = C.paramname;
269
270}
271
272template<class T>
273CGA<T> CGA<T>::operator=(CGA<T> &C)
274{
275 GA_params = C.GA_params;
276 filenames = C.filenames;
277 Ind = C.Ind;
278 Ind_old = C.Ind_old;
279 params = C.params;
280 loged = C.loged;
281 fitdist = C.fitdist;
282 MaxFitness = C.MaxFitness;
283 paramname = C.paramname;
284
285 return *this;
286
287}
288
289template<class T>
290CGA<T>::~CGA()
291{
292
293}
294
295template<class T>
296void CGA<T>::initialize()
297{
298 for (int i=0; i<GA_params.maxpop; i++)
299 {
300 Ind[i].initialize();
301 }
302
303 if (filenames.initialpopfilemame!="")
304 {
305 getinifromoutput(filenames.pathname+filenames.initialpopfilemame);
306 for (int i=0; i<initial_pop.size(); i++)
307 for (int j=0; j<max(int(initial_pop[i].size()),GA_params.nParam); j++)
308 if (loged[j]==1)
309 Ind[i].x[j] = log10(initial_pop[i][j]);
310 else
311 Ind[i].x[j] = initial_pop[i][j];
312 }
313}
314
315template<class T>
316void CGA<T>::Setminmax(int a, double minrange, double maxrange, int prec)
317{
318 for (int i=0; i<GA_params.maxpop; i++)
319 {
320 Ind[i].maxrange[a] = maxrange;
321 Ind[i].minrange[a] = minrange;
322 Ind[i].precision[a] = prec;
323 }
324
325}
326
327template<class T>
328void CGA<T>::assignfitnesses()
329{
330 sumfitness = 0;
331
332 vector<vector<double>> inp;
333
334 inp.resize(GA_params.maxpop);
335
336
337 for (int k=0; k<GA_params.maxpop; k++)
338 inp[k].resize(GA_params.nParam);
339
340 vector<double> time_(GA_params.maxpop);
341 vector<int> epochs(GA_params.maxpop);
342 clock_t t0,t1;
343 //Models.clear();
344
345 for (int k = 0; k < GA_params.maxpop; k++)
346 {
347 for (int i = 0; i < GA_params.nParam; i++)
348 {
349 if (loged[i] != 1)
350 {
351 inp[k][i] = Ind[k].x[i]; //Ind
352 }
353 else
354 {
355 inp[k][i] = pow(10, Ind[k].x[i]);
356 }
357 }
358
359 Ind[k].actual_fitness = 0;
360
361 //Models[k] = *Model;
362 //Models[k].SetSilent(true);
363 //Models[k].SetRecordResults(false);
364 //Models[k].SetNumThreads(1);
365 for (int i = 0; i < GA_params.nParam; i++)
366 Models[k].SetParameterValue(i, inp[k][i]);
367 //Models[k].ApplyParameters();
368
369 }
370
371
372#ifndef NO_OPENMP
373 omp_set_num_threads(numberOfThreads);
374#endif
375int counter=0;
376#pragma omp parallel for //private(ts,l)
377 for (int k=0; k<GA_params.maxpop; k++)
378 {
379
380 if (GA_params.Steepest_Descent && k<min(GA_params.maxpop/10,1))
381 {
382 if (k==0)
383 qDebug()<<"Prior Likelihood: " <<Models[k].GetObjectiveFunctionValue();
384 CVector updated_params;
385 for (int j=0; j<5; j++)
386 updated_params = Models[k].GradientUpdate();
387 for (int i = 0; i < GA_params.nParam; i++)
388 {
389 if (loged[i] != 1)
390 {
391 inp[k][i] = updated_params[i];
392 Ind[k].x[i] = updated_params[i];
393
394 }
395 else
396 {
397 inp[k][i] = updated_params[i];
398 Ind[k].x[i] = log10(updated_params[i]);
399 }
400 }
401
402 if (k==0)
403 qDebug()<<"Posterior Likelihood: " <<Models[k].GetObjectiveFunctionValue();
404
405 }
406
407 FILE *FileOut;
408#pragma omp critical
409 {
410 FileOut = fopen((filenames.pathname+"detail_GA.txt").c_str(),"a");
411 fprintf(FileOut, "%i, ", k);
412 for (int l=0; l<Ind[0].nParams; l++)
413 if (loged[l]==1)
414 fprintf(FileOut, "%le, ", pow(10,Ind[k].x[l]));
415 else
416 fprintf(FileOut, "%le, ", Ind[k].x[l]);
417
418 //fprintf(FileOut, "%le, %le, %i, %e, %i, %i", Ind[k].actual_fitness, Ind[k].fitness, Ind[k].rank, time_[k], threads_num[k],num_threads[k]);
419 //fprintf(FileOut, "%le, %le, %i, %e", Ind[k].actual_fitness, Ind[k].fitness, Ind[k].rank, time_[k]);
420 fprintf(FileOut, "\n");
421 fclose(FileOut);
422 }
423 time_t t0 = time(nullptr);
424
425
426
427 Ind[k].actual_fitness = Models[k].GetObjectiveFunctionValue();
428 //for (unsigned int i=0; i<Models[k].fit_measures.size(); i++)
429 // Ind[k].fit_measures[i] = Models[k].fit_measures[i];
430
431// epochs[k] += Models[k].EpochCount();
432 time_[k] = time(nullptr)-t0;
433 counter++;
434#pragma omp critical
435 {
436#ifdef Q_GUI_SUPPORT
437 if (rtw != nullptr)
438 {
439#ifndef NO_OPENMP
440 if (omp_get_thread_num() == 0)
441#endif
442 {
443// rtw->SetProgress2(double(counter + 1) / GA_params.maxpop);
444// QCoreApplication::processEvents();
445 }
446 }
447#endif
448
449 {
450// FileOut = fopen((filenames.pathname+"detail_GA.txt").c_str(),"a");
451// fprintf(FileOut, "%i, fitness=%e, time=%e, internal_time=%e, failed=%i\n", k, Ind[k].actual_fitness, time_[k], double(Models[k].GetSimulationDuration()), Models[k].GetSolutionFailed());
452// fclose(FileOut);
453 }
454
455 }
456 }
457 Model_out = Models[maxfitness()];
458#ifdef Q_version
459// if (rtw != nullptr)
460// {
461// rtw->SetProgress(1);
462// QCoreApplication::processEvents();
463// }
464#endif
465 inp.clear();
466 assignfitness_rank(GA_params.N);
467
468}
469
470template<class T>
471void CGA<T>::crossover()
472{
473
474 Ind_old = Ind;
475 int a = maxfitness();
476 Ind[0] = Ind_old[a];
477 Ind[1] = Ind_old[a];
478 for (int i=2; i<GA_params.maxpop; i+=2)
479 {
481 int j1 = fitdist.GetRand();
482 int j2 = fitdist.GetRand();
483 double x = fitdist.GetRndUniF(0,1);
484 if (x<GA_params.pcross)
485 if (GA_params.cross_over_type == 1)
486 cross(Ind_old[j1], Ind_old[j2], Ind[i], Ind[min(i+1,GA_params.maxpop-1)]); //1 Breaking point
487 else
488 cross2p(Ind_old[j1], Ind_old[j2], Ind[i], Ind[min(i + 1, GA_params.maxpop - 1)]); //2 Breaking point
489 else
490 {
491 Ind[i] = Ind_old[j1];
492 Ind[i+1] = Ind_old[j2];
493 }
494
495 }
496
497}
498
499template<class T>
500void CGA<T>::crossoverRC()
501{
502
503 for (int i=0; i<GA_params.maxpop; i++)
504 Ind_old[i] = Ind[i];
505 int a = maxfitness();
506 Ind[0] = Ind_old[a];
507 Ind[1] = Ind_old[a];
508 for (int i=2; i<GA_params.maxpop; i+=2)
509 {
510 int j1 = fitdist.GetRand();
511 int j2 = fitdist.GetRand();
512 double x = GetRndUnif(0,1);
513 if (x<GA_params.pcross)
514 cross_RC_L(Ind_old[j1], Ind_old[j2], Ind[i], Ind[i+1]);
515 else
516 {
517 Ind[i] = Ind_old[j1];
518 Ind[i+1] = Ind_old[j2];
519 }
520 }
521}
522
523template<class T>
524bool CGA<T>::SetProperty(const string &varname, const string &value)
525{
526 if (aquiutils::tolower(varname) == "maxpop" || varname == "Population") {GA_params.maxpop = aquiutils::atoi(value); setnumpop(GA_params.maxpop); return true;}
527 if (aquiutils::tolower(varname) == "ngen" || varname == "Number of Generations") {GA_params.nGen = aquiutils::atoi(value); return true;}
528 if (aquiutils::tolower(varname) == "pcross" || varname == "Cross-over probability") {GA_params.pcross = aquiutils::atof(value); return true;}
529 if (aquiutils::tolower(varname) == "pmute" || varname == "Mutation probability") {GA_params.pmute = aquiutils::atof(value); return true;}
530 if (aquiutils::tolower(varname) == "shakescale" || varname == "Shake coefficient") {GA_params.shakescale = aquiutils::atof(value); return true;}
531 if (aquiutils::tolower(varname) == "shakescalered" || varname == "Shake coefficient reduction factor") {GA_params.shakescalered = aquiutils::atof(value); return true;}
532 if (aquiutils::tolower(varname) == "outputfile" || varname == "GA output file") {filenames.outputfilename = value; return true;}
533 if (aquiutils::tolower(varname) == "getfromfilename") {filenames.getfromfilename = value.c_str(); return true;}
534 if (aquiutils::tolower(varname) == "initial_population") {filenames.initialpopfilemame = value; return true;}
535 if (aquiutils::tolower(varname) == "numthreads" || varname == "Number of threads to be used") {numberOfThreads = aquiutils::atoi(value.c_str()); return true;}
536 if (aquiutils::tolower(varname) == "steepest descent") {
537
538 if (aquiutils::tolower(value) == "true")
539 GA_params.Steepest_Descent = true;
540 else
541 GA_params.Steepest_Descent = false;
542 return true;
543 }
544 last_error = "Property '" + varname + "' was not found!";
545 return false;
546}
547
548template<class T>
549bool CGA<T>::SetProperties(const map<string,string> &arguments)
550{
551 for (map<string,string>::const_iterator it=arguments.begin(); it!=arguments.end(); it++)
552 SetProperty(it->first, it->second);
553 return true;
554}
555
556template<class T>
557double CGA<T>::avgfitness()
558{
559 double sum=0;
560 for (int i=0; i<GA_params.maxpop; i++)
561 sum += Ind[i].fitness;
562 return sum/GA_params.maxpop;
563}
564
565template<class T>
566void CGA<T>::write_to_detailed_GA(string s)
567{
568 FILE *FileOut;
569
570 FileOut = fopen((filenames.pathname + "detail_GA.txt").c_str(), "a");
571 fprintf(FileOut, "%s\n", s.c_str());
572 fclose(FileOut);
573
574}
575
576template<class T>
577int CGA<T>::optimize()
578{
579 #ifdef Q_version
580 QCoreApplication::processEvents();
581 #endif // Q_version
582 string RunFileName = filenames.pathname + filenames.outputfilename;
583
584 FILE *FileOut;
585 FILE *FileOut1;
586
587 FileOut = fopen(RunFileName.c_str(),"w");
588 fclose(FileOut);
589 FileOut1 = fopen((filenames.pathname + "detail_GA.txt").c_str(), "w");
590 fclose(FileOut1);
591
592 double shakescaleini = GA_params.shakescale;
593
594 vector<double> X(Ind[0].nParams);
595
596 initialize();
597 double ininumenhancements = GA_params.numenhancements;
598 GA_params.numenhancements = 0;
599
600 CMatrix Fitness(GA_params.nGen, 3);
601
602 Models.clear();
603 Models.resize(GA_params.maxpop);
604 for (int k=0; k<GA_params.maxpop; k++)
605 Models[k] = *Model;
606 for (current_generation=0; current_generation<GA_params.nGen; current_generation++)
607 {
608
609 write_to_detailed_GA("Assigning fitnesses ...");
610
611 assignfitnesses();
612
613 write_to_detailed_GA("Assigning fitnesses done!");
614 FileOut = fopen(RunFileName.c_str(),"a");
615 printf("Generation: %i\n", current_generation);
616 fprintf(FileOut, "Generation: %i\n", current_generation);
617 fprintf(FileOut, "ID, ");
618 for (int k=0; k<Ind[0].nParams; k++)
619 fprintf(FileOut, "%s, ", paramname[k].c_str());
620
621 //fprintf(FileOut, "%s, %s, %s, ", "likelihood", "Fitness", "Rank");
622 for (unsigned int i=0; i<Model->ObservationsCount();i++)
623 {
624 fprintf(FileOut, "%s, %s, %s", (Model->observation(i)->GetName()+"_MSE").c_str(), (Model->observation(i)->GetName()+"_R2").c_str(), (Model->observation(i)->GetName()+"_NSE").c_str());
625 }
626 fprintf(FileOut, "\n");
627 write_to_detailed_GA("Generation: " + aquiutils::numbertostring(current_generation));
628 for (int j1=0; j1<GA_params.maxpop; j1++)
629 {
630
631 fprintf(FileOut, "%i, ", j1);
632
633 for (int k=0; k<Ind[0].nParams; k++)
634 if (loged[k] == 1)
635 fprintf(FileOut, "%le, ", pow(10, Ind[j1].x[k]));
636 else
637 fprintf(FileOut, "%le, ", Ind[j1].x[k]);
638
639 fprintf(FileOut, "%le, %le, %i, ", Ind[j1].actual_fitness, Ind[j1].fitness, Ind[j1].rank);
640 for (unsigned int i=0; i<Model->ObservationsCount();i++)
641 {
642 fprintf(FileOut, ",%le, %le, %le", Ind[j1].fit_measures[i]);
643 }
644 fprintf(FileOut, "\n");
645 }
646 fclose(FileOut);
647
648 int j = maxfitness();
649
650 Fitness[current_generation][0] = Ind[j].actual_fitness;
651
652#ifdef Q_GUI_SUPPORT
653 if (rtw)
654 { if (current_generation==0)
655 {
656 rtw->SetYRange(0,Ind[j].actual_fitness*1.1);
657 rtw->SetXRange(0,GA_params.nGen);
658 }
659 rtw->SetProgress(double(current_generation)/double(GA_params.nGen));
660 rtw->AppendPoint(current_generation+1,Ind[j].actual_fitness);
661 //rtw->Replot();
662 QCoreApplication::processEvents();
663 }
664#endif
665 if (current_generation>10)
666 {
667 if ((Fitness[current_generation][0] == Fitness[current_generation - 3][0]) && GA_params.shakescale>pow(10.0, -Ind[0].precision[0]))
668 GA_params.shakescale *= GA_params.shakescalered;
669
670
671 if ((Fitness[current_generation][0]>Fitness[current_generation - 1][0]) && (GA_params.shakescale<shakescaleini))
672 GA_params.shakescale /= GA_params.shakescalered;
673 GA_params.numenhancements = 0;
674 }
675
676 if (current_generation>50)
677 {
678 if (Fitness[current_generation][0] == Fitness[current_generation - 20][0])
679 {
680 GA_params.numenhancements *= 1.05;
681 if (GA_params.numenhancements == 0) GA_params.numenhancements = ininumenhancements;
682 }
683
684 if (Fitness[current_generation][0] == Fitness[current_generation - 50][0])
685 GA_params.numenhancements = ininumenhancements * 10;
686 }
687
688 Fitness[current_generation][1] = GA_params.shakescale;
689 Fitness[current_generation][2] = GA_params.pmute;
690
691 if (current_generation>20)
692 {
693 if (GA_params.shakescale == Fitness[current_generation - 20][1])
694 GA_params.shakescale = shakescaleini;
695 }
696
697
698 j = maxfitness();
699 MaxFitness = Ind[j].actual_fitness;
700
701 Fitness[current_generation][0] = Ind[j].actual_fitness;
702
703
704 fillfitdist();
705
706 write_to_detailed_GA("Cross-over ...");
707
708 if (GA_params.RCGA == true)
709 crossoverRC();
710 else
711 crossover();
712
713 write_to_detailed_GA("Cross-over done! ");
714
715 write_to_detailed_GA("Mutation ...");
716
717 mutate(GA_params.pmute);
718 write_to_detailed_GA("Mutation done!");
719 write_to_detailed_GA("Shake...!");
720 shake();
721 write_to_detailed_GA("Shake done!");
722
723
724 }
725
726 assignfitnesses();
727 FileOut = fopen(RunFileName.c_str(), "a");
728 fprintf(FileOut, "Final Enhancements\n");
729
730 int j = maxfitness();
731
732 MaxFitness = Ind[j].actual_fitness;
733 final_params.resize(GA_params.nParam);
734
735
736 for (int k = 0; k<Ind[0].nParams; k++)
737 {
738 if (loged[k] == 1) final_params[k] = pow(10, Ind[j].x[k]); else final_params[k] = Ind[j].x[k];
739 fprintf(FileOut, "%s, ", paramname[k].c_str());
740 fprintf(FileOut, "%le, ", final_params[k]);
741 fprintf(FileOut, "%le, %le\n", Ind[j].actual_fitness, Ind[j].fitness);
742 }
743 for (unsigned int i=0; i<Model->ObservationsCount();i++)
744 {
745 fprintf(FileOut, ",%le, %le, %le\n", Ind[j].fit_measures[i]);
746 }
747 fclose(FileOut);
748
749 assignfitnesses(final_params);
750
751#ifdef Q_GUI_SUPPORT
752 if (rtw)
753 {
754 rtw->SetProgress(1.0);
755 QCoreApplication::processEvents();
756 }
757#endif
758
759 Models.clear();
760
761 return maxfitness();
762}
763
764template<class T>
765double CGA<T>::assignfitnesses(vector<double> inp)
766{
767
768 double likelihood = 0;
769 T Model1;
770 Model1 = *Model;
771
772 for (int i = 0; i < GA_params.nParam; i++)
773 Model1.SetParameterValue(i, inp[i]);
774
775 //Model1.ApplyParameters();
776 //Model1.Solve();
777 likelihood -= Model1.GetObjectiveFunctionValue();
778
779 Model_out = Model1;
780 //Model_out.TransferResultsFrom(&Model1);
781 return likelihood;
782
783}
784
785/*
786template<class T>
787vector<T>& CGA<T>::assignfitnesses_p(vector<double> inp) //Generates an instance of the model with the provided parameters
788{
789 double likelihood = 0;
790 vector<T> Models(1);
791 for (int ts = 0; ts<1; ts++)
792 {
793 Models[ts] = Model;
794
795 int l = 0;
796 for (int i = 0; i<GA_params.nParam; i++)
797 Models[ts].SetParam(params[i], inp[getparamno(i, ts)]);
798 Models[ts].FinalizeSetParams();
799 likelihood -= Models[ts].EvaluateObjectiveFunction();
800 }
801 return Models;
802}
803
804template<class T>
805int CGA<T>::get_act_paramno(int i)
806{
807 int l = -1;
808 for (int j = 0; j<GA_params.nParam; j++)
809 {
810 if (apply_to_all[j]) l++; else l += 1;
811 if (l >= i)
812 {
813 if (apply_to_all[j]) l -= 1; else l--;
814 return j;
815 }
816 }
817}
818*/
819template<class T>
820double CGA<T>::getfromoutput(string filename)
821{
822 ifstream file(filename);
823 vector<string> s;
824 final_params.resize(GA_params.nParam);
825 while (file.eof() == false)
826 {
827 s = aquiutils::getline(file);
828 if (s.size()>0)
829 {
830 if (s[0] == "Final Enhancements")
831 for (int i = 0; i<GA_params.nParam; i++)
832 {
833 s = aquiutils::getline(file);
834 if (s.size() == 0)
835 write_to_detailed_GA("The number of parameters in GA output file does not match the number of unknown parameters");
836 else
837 final_params[i] = atof(s[1].c_str());
838 }
839 }
840 }
841 double ret = assignfitnesses(final_params);
842 return ret;
843}
844
845
846template<class T>
847int CGA<T>::getparamno(int i, int ts)
848{
849 int l = 0;
850 for (int j = 0; j<i; j++)
851 if (apply_to_all[j]) l++; else l += 1;
852
853 if (apply_to_all[i])
854 return l;
855 else
856 return l + ts;
857
858}
859template<class T>
860int CGA<T>::get_time_series(int i)
861{
862 int l = 0;
863 for (int j = 0; j<GA_params.nParam; j++)
864 {
865 if (apply_to_all[j]) l += 1; else l += 1;
866 if (l >= i)
867 {
868 if (apply_to_all[j]) l -= 1; else l--;
869 return i - l;
870 }
871 }
872}
873
874
875template<class T>
876void CGA<T>::shake()
877{
878 for (int i=1; i<GA_params.maxpop; i++)
879 Ind[i].shake(GA_params.shakescale);
880
881}
882
883template<class T>
884void CGA<T>::mutate(double mu)
885{
886 for (int i=2; i<GA_params.maxpop; i++)
887 Ind[i].mutate(mu);
888
889}
890
891template<class T>
892int CGA<T>::maxfitness()
893{
894 double max_fitness = 1E+308 ;
895 int i_max = 0;
896 for (int i=0; i<GA_params.maxpop; i++)
897 if (max_fitness>Ind[i].actual_fitness)
898 {
899 max_fitness = Ind[i].actual_fitness;
900 i_max = i;
901 }
902 return i_max;
903
904}
905
906template<class T>
907double CGA<T>::variancefitness()
908{
909 double sum = 0;
910 double a = avgfitness();
911 for (int i=0; i<GA_params.maxpop; i++)
912 sum += (a - Ind[i].fitness)*(a - Ind[i].fitness);
913 return sum;
914
915}
916
917template<class T>
918double CGA<T>::stdfitness()
919{
920 double sum = 0;
921 double a = avg_inv_actual_fitness();
922 for (int i=0; i<GA_params.maxpop; i++)
923 sum += (a - 1/Ind[i].actual_fitness)*(a - 1/Ind[i].actual_fitness);
924 return sqrt(sum)/GA_params.maxpop/a;
925
926}
927
928template<class T>
929double CGA<T>::avg_actual_fitness()
930{
931 double sum=0;
932 for (int i=0; i<GA_params.maxpop; i++)
933 sum += Ind[i].actual_fitness;
934 return sum/GA_params.maxpop;
935
936}
937
938template<class T>
939double CGA<T>::avg_inv_actual_fitness()
940{
941 double sum=0;
942 for (int i=0; i<GA_params.maxpop; i++)
943 sum += 1/Ind[i].actual_fitness;
944 return sum/GA_params.maxpop;
945
946}
947
948
949template<class T>
950void CGA<T>::assignrank()
951{
952 for (int i=0; i<GA_params.maxpop; i++)
953 {
954 int r = 1;
955 for (int j=0; j<GA_params.maxpop; j++)
956 {
957 if (Ind[i].actual_fitness > Ind[j].actual_fitness) r++;
958 }
959 Ind[i].rank = r;
960 }
961
962}
963
964template<class T>
965void CGA<T>::assignfitness_rank(double N)
966{
967 assignrank();
968 for (int i=0; i<GA_params.maxpop; i++)
969 {
970 Ind[i].fitness = pow(1.0/static_cast<double>(Ind[i].rank),GA_params.N);
971 }
972}
973
974
975template<class T>
976void CGA<T>::fillfitdist()
977{
978 double sum=0;
979 for (int i=0; i<GA_params.maxpop; i++)
980 {
981 sum+=Ind[i].fitness;
982 }
983
984 fitdist.s[0] = 0;
985 fitdist.e[0] = Ind[0].fitness/sum;
986 for (int i=1; i<GA_params.maxpop-1; i++)
987 {
988 fitdist.e[i] = fitdist.e[i-1] + Ind[i].fitness/sum;
989 fitdist.s[i] = fitdist.e[i-1];
990 }
991 fitdist.s[GA_params.maxpop-1] = fitdist.e[GA_params.maxpop-2];
992 fitdist.e[GA_params.maxpop-1] = 1;
993
994}
995
996template<class T>
997double CGA<T>::evaluateforward()
998{
999 vector<double> v(1);
1000 v[0] = 1;
1001 CVector out(1);
1002 int x_nParam = GA_params.nParam;
1003 vector<int> x_params = params;
1004 GA_params.nParam = 1;
1005 params.resize(1);
1006 params[0] = 100;
1007 out[0] = assignfitnesses(v);
1008
1009 out.writetofile(filenames.pathname + "likelihood.txt");
1010 params = x_params;
1011 GA_params.nParam = x_nParam;
1012 return out[0];
1013}
1014
1015template<class T>
1016void CGA<T>::getinifromoutput(string filename)
1017{
1018 ifstream file(filename);
1019 vector<string> s;
1020 initial_pop.resize(1);
1021 initial_pop[0].resize(GA_params.nParam);
1022 while (file.eof() == false)
1023 {
1024 s = aquiutils::getline(file);
1025 if (s.size()>0)
1026 { if (s[0] == "Final Enhancements")
1027 for (int i=0; i<GA_params.nParam; i++)
1028 {
1029 s = aquiutils::getline(file);
1030 if (loged[i]==1)
1031 initial_pop[0][i] = atof(s[1].c_str());
1032 else
1033 initial_pop[0][i] = atof(s[1].c_str());
1034
1035 }
1036 }
1037 }
1038
1039}
1040
1041template<class T>
1042void CGA<T>::getinitialpop(string filename)
1043{
1044 ifstream file(filename);
1045 vector<string> s;
1046
1047 while (file.eof() == false)
1048 {
1049 s = aquiutils::getline(file);
1050
1051 if (s.size()>0)
1052 {
1053 vector<double> x;
1054 for (int j=0; j<s.size(); j++)
1055 initial_pop.push_back(aquiutils::ATOF(s));
1056
1057 }
1058 }
1059 file.close();
1060}
1061
1062
cross2p
void cross2p(CBinary &B1, CBinary &B2, int p1, int p2)
Definition Binary.cpp:190
cross
void cross(CBinary &B1, CBinary &B2, int p)
Definition Binary.cpp:142
cross_RC_L
void cross_RC_L(const CIndividual &I1, const CIndividual &I2, CIndividual &IR1, CIndividual &IR2)
Definition Individual.cpp:188
GetRndUnif
double GetRndUnif(double xmin, double xmax)
Definition Individual.cpp:126
CGA
Genetic Algorithm optimizer for global parameter estimation.
Definition GA.h:350
CGA::MaxFitness
double MaxFitness
Maximum fitness value in current population.
Definition GA.h:366
CGA::fillfitdist
void fillfitdist()
Fill fitness distribution for roulette wheel selection.
Definition GA.cpp:754
CGA::SetProperty
bool SetProperty(const std::string &varname, const std::string &value)
Set GA property from string key-value pair.
Definition GA.hpp:524
CGA::params
std::vector< int > params
Parameter indices being optimized.
Definition GA.h:415
CGA::maxfitness
int maxfitness()
Find index of individual with maximum fitness.
Definition GA.cpp:670
CGA::assignrank
void assignrank()
Assign ranks to individuals based on fitness.
Definition GA.cpp:728
CGA::get_time_series
int get_time_series(int i)
Get time series index for parameter.
Definition GA.cpp:638
CGA::avg_actual_fitness
double avg_actual_fitness()
Calculate average of actual (untransformed) fitness values.
Definition GA.cpp:707
CGA::SetProperties
bool SetProperties(const std::map< std::string, std::string > &arguments)
Set multiple GA properties from map.
Definition GA.hpp:549
CGA::GA_params
GA_Tweaking_parameters GA_params
GA algorithm configuration parameters.
Definition GA.h:374
CGA::shake
void shake()
Apply shake/perturbation to population.
Definition GA.cpp:654
CGA::~CGA
virtual ~CGA()
Destructor.
Definition GA.cpp:203
CGA::Ind
std::vector< CIndividual > Ind
Current population of individuals.
Definition GA.h:475
CGA::fitdist
GADistribution fitdist
Fitness distribution for parent selection.
Definition GA.h:986
CGA::operator=
CGA operator=(CGA &C)
Assignment operator.
Definition GA.cpp:186
CGA::setnumpop
void setnumpop(int n)
Set population size.
Definition GA.cpp:142
CGA::evaluateforward
double evaluateforward()
Evaluate model forward (compute predictions)
Definition GA.cpp:775
CGA::paramname
std::vector< std::string > paramname
Names of parameters being optimized.
Definition GA.h:491
CGA::Setminmax
void Setminmax(int a, double minrange, double maxrange, int prec)
Set parameter bounds and precision.
Definition GA.cpp:229
CGA::avg_inv_actual_fitness
double avg_inv_actual_fitness()
Calculate average of inverse actual fitness.
Definition GA.cpp:717
CGA::crossoverRC
void crossoverRC()
Real-coded crossover operation.
Definition GA.cpp:351
CGA::CGA
CGA()
Default constructor.
Definition GA.cpp:12
CGA::assignfitnesses
void assignfitnesses()
Evaluate fitness for all individuals in population.
Definition GA.cpp:241
CGA::setnparams
void setnparams(int n)
Set number of parameters.
Definition GA.cpp:129
CGA::optimize
int optimize()
Main optimization loop.
Definition GA.cpp:394
CGA::getinifromoutput
void getinifromoutput(std::string filename)
Initialize population from previous output file.
Definition GA.cpp:794
CGA::stdfitness
double stdfitness()
Calculate standard deviation of population fitness.
Definition GA.cpp:696
CGA::getparamno
int getparamno(int i, int ts)
Get parameter index for variable and time series.
Definition GA.cpp:625
CGA::mutate
void mutate(double mu)
Apply mutation to population.
Definition GA.cpp:662
CGA::Ind_old
std::vector< CIndividual > Ind_old
Previous generation's population.
Definition GA.h:483
CGA::initialize
void initialize()
Initialize GA structures and allocate memory.
Definition GA.cpp:209
CGA::crossover
void crossover()
Perform crossover operation on population.
Definition GA.cpp:322
CGA::variancefitness
double variancefitness()
Calculate variance of population fitness.
Definition GA.cpp:685
CGA::loged
std::vector< int > loged
Flags indicating if parameters are log-transformed.
Definition GA.h:423
CGA::getfromoutput
double getfromoutput(std::string filename)
Read best solution from previous output file.
Definition GA.cpp:598
CGA::write_to_detailed_GA
void write_to_detailed_GA(std::string s)
Write detailed GA information to file.
Definition GA.cpp:383
CGA::InitiatePopulation
void InitiatePopulation()
Create initial population.
Definition GA.hpp:165
CGA::assignfitness_rank
void assignfitness_rank(double N)
Assign fitness using rank-based scheme.
Definition GA.cpp:743
CGA::getinitialpop
void getinitialpop(std::string filename)
Read initial population from file.
Definition GA.cpp:820
CGA::avgfitness
double avgfitness()
Calculate average fitness of population.
Definition GA.cpp:374
CGA::filenames
_filenames filenames
File paths for GA input/output.
Definition GA.h:382
CIndividual::precision
std::vector< int > precision
Definition Individual.h:22
CIndividual::minrange
std::vector< double > minrange
Definition Individual.h:23
CIndividual::maxrange
std::vector< double > maxrange
Definition Individual.h:23
distribution_type::lognormal
@ lognormal
Lognormal distribution: ln(x) ~ N(μ, σ²), for strictly positive variables.
_range::low
@ low
Lower bound of the parameter range.
_range::high
@ high
Upper bound of the parameter range.
GA_Tweaking_parameters::maxpop
int maxpop
Maximum population size (number of individuals)
Definition GA.h:45