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