Heavy ion Analysis Libriares
Loading...
Searching...
No Matches
Minimizer.cxx
1/*
2 * HalMinimizer.cxx
3 *
4 * Created on: 7 lis 2019
5 * Author: Daniel Wielanek
6 * E-mail: daniel.wielanek@gmail.com
7 * Warsaw University of Technology, Faculty of Physics
8 */
9
10#include "Minimizer.h"
11
12#include "Array.h"
13#include "Cout.h"
14#include "Std.h"
15#include <Math/IFunctionfwd.h>
16#include <RtypesCore.h>
17#include <TMath.h>
18#include <TMatrixD.h>
19#include <TString.h>
20#include <fstream>
21#include <string>
22
23namespace Hal {
24 Minimizer* Minimizer::fgInstance = NULL;
25 Minimizer::Minimizer() :
26 fFunc(nullptr),
27 fMinimizeType(kScan),
28 fNCalls(0),
29 fDiscreteFit(kFALSE),
30 fTrace(kFALSE),
31 fGlobMin(1.),
32 fNDF(0),
33 fQuantumFits(nullptr),
34 fSmoothFits(nullptr),
35 fErrors(nullptr),
36 fTempParams(nullptr) {}
37
38 // set param methods
39
40 bool Minimizer::SetFixedVariable(unsigned int ivar, const std::string& name, double val) {
41 if (fParameters.size() <= ivar) { fParameters.resize(ivar + 1); }
42 fParameters[ivar].SetParName(name);
43 fParameters[ivar].SetRange(val, val);
44 fParameters[ivar].SetIsFixed(kTRUE);
45 return true;
46 }
47
48 bool Minimizer::SetVariableLimits(unsigned int ivar, double lower, double upper) {
49 if (lower > upper) {
50 double temp = upper;
51 upper = lower;
52 lower = temp;
53 }
54 if (fParameters.size() <= ivar) { fParameters.resize(ivar + 1); }
55 fParameters[ivar].SetRange(lower, upper);
56 fParameters[ivar].SetIsFixed(kFALSE);
57 return true;
58 }
59
60 bool
61 Minimizer::SetLimitedVariable(unsigned int ivar, const std::string& name, double val, double step, double lower, double upper) {
62 if (fParameters.size() <= ivar) { fParameters.resize(ivar + 1); }
63 fParameters[ivar].SetParName(name);
64 fParameters[ivar].SetRange(lower, upper);
65 fParameters[ivar].SetStartVal(val);
66 fParameters[ivar].SetMapRange(lower, upper, TMath::Nint((upper - lower) / step + 1));
67 fParameters[ivar].SetIsFixed(kFALSE);
68 return false;
69 }
70
71 // recalc methods
72
73 void Minimizer::InitFit() {
74 fNonConstMap.clear();
75 fAllowedValues.clear();
76 fStateVector.resize(GetNParams());
77 fStateVectorMin.resize(GetNParams());
78 for (auto& x : fStateVector)
79 x = 0;
80 for (auto& x : fStateVectorMin)
81 x = 0;
82 fNCalls = 0;
83 fFreePars = 0;
84 for (unsigned int i = 0; i < fParameters.size(); i++) {
85 if (!fParameters[i].IsFixed()) {
86 fNonConstMap.push_back(i);
87 } else {
88 fFreePars++;
89 }
90 }
91 auto lambda = [](Int_t a, Int_t b) -> bool { return a < b; };
92 std::sort(fNonConstMap.begin(), fNonConstMap.end(), lambda);
93 Cout::PrintInfo("BEFORE INIT ", EInfo::kLowWarning);
94 std::cout << Cout::GetColor(kBlue);
95 Cout::Database({"ParName", "MinMap", "MaxMap", "Points", "Step", "Min", "Max"});
96 std::cout << Cout::GetDisableColor();
97 for (unsigned int i = 0; i < fParameters.size(); i++) {
98 if (fParameters[i].IsFixed()) { std::cout << Cout::GetColor(kOrange); }
99 Cout::Database({fParameters[i].GetParName().Data(),
100 Form("%4.4f", fParameters[i].GetMapMin()),
101 Form("%4.4f", fParameters[i].GetMapMax()),
102 Form("%d", fParameters[i].GetNPoints()),
103 Form("%4.4f", fParameters[i].GetDParam()),
104 Form("%4.4f", fParameters[i].GetMin()),
105 Form("%4.4f", fParameters[i].GetMax())});
106 if (fParameters[i].IsFixed()) std::cout << Cout::GetDisableColor();
107 }
108
109 // round minimum
110 for (unsigned int i = 0; i < fParameters.size(); i++) {
111 fParameters[i].SetIsDiscrete(kTRUE); // all parameters must be discrete
112 fParameters[i].Init();
113 // fParameters[i].Print();
114 }
115 Cout::PrintInfo("AFTER INIT ", EInfo::kLowWarning);
116 std::cout << Cout::GetColor(kBlue);
117 Cout::Database({"ParName", "MinMap", "MaxMap", "Points", "Step", "Min", "Max"});
118 std::cout << Cout::GetDisableColor();
119 for (unsigned int i = 0; i < fParameters.size(); i++) {
120 if (fParameters[i].IsFixed()) { std::cout << Cout::GetColor(kOrange); }
121 Cout::Database({fParameters[i].GetParName().Data(),
122 Form("%4.4f", fParameters[i].GetMapMin()),
123 Form("%4.4f", fParameters[i].GetMapMax()),
124 Form("%d", fParameters[i].GetNPoints()),
125 Form("%4.4f", fParameters[i].GetDParam()),
126 Form("%4.4f", fParameters[i].GetMin()),
127 Form("%4.4f", fParameters[i].GetMax())});
128 if (fParameters[i].IsFixed()) std::cout << Cout::GetDisableColor();
129 }
130 if (fQuantumFits) {
131 delete[] fQuantumFits;
132 delete[] fSmoothFits;
133 delete[] fErrors;
134 delete[] fTempParams;
135 }
136 fQuantumFits = new double[GetNParams()];
137 fSmoothFits = new double[GetNParams()];
138 fErrors = new double[GetNParams()];
139 fParamsMin.resize(GetNParams());
140 fTempParams = new double[GetNParams()];
141 for (int i = 0; i < GetNParams(); i++) {
142 if (fParameters[i].IsFixed()) { fTempParams[i] = fParameters[i].GetStartVal(); }
143 }
144 fGlobMin = DBL_MAX;
145 for (int i = 0; i < GetNParams(); i++) {
146 fAllowedValues.push_back(fParameters[i].GetValuesArray());
147 }
148 }
149
150 void Minimizer::Reset() { fParameters.clear(); }
151
152
153 // minimize methods
154
155
156 bool Minimizer::Minimize() {
157 if (GetNParams() == 0) return false;
158 InitFit();
159 switch (fMinimizeType) {
160 case eMinimizeType::kAnt: {
161 /* LoopOverParameter(0);
162 for (unsigned int i = 0; i < GetNParams(); i++) {
163 fTempParams[i] = fParamsMin[i];
164 fSmoothFits[i] = fParamsMin[i];
165 }
166 std::cout << "*************************" << std::endl;
167 EstimateErrors();
168 FinishFit();*/
169 std::cout << "Ant Minimiser" << std::endl;
170 MinimizeAnt();
171 } break;
172 case eMinimizeType::kScan: {
173 // MinimizeAnt();
174 std::cout << "Scan minimizer" << std::endl;
175 MinimizeScan();
176 } break;
177 }
178 FinishFit();
179 return true;
180 }
181
182 const double* Minimizer::X() const {
183 if (!fDiscreteFit) return fSmoothFits;
184 return fQuantumFits;
185 }
186
187 void Minimizer::SetFunction(const ROOT::Math::IMultiGenFunction& func) {
188 fFunc = const_cast<ROOT::Math::IMultiGenFunction*>(&func);
189 }
190
191 void Minimizer::SetParamConf(const MinimizerStepConf& conf, Bool_t overwrite) { conf.SetParameters(fParameters, overwrite); }
192
193 Minimizer::~Minimizer() {
194 if (fQuantumFits) {
195 delete[] fQuantumFits;
196 delete[] fSmoothFits;
197 delete[] fErrors;
198 delete[] fTempParams;
199 }
200 }
201
202 Bool_t Minimizer::LoopOverParameter(Int_t param) {
204 if ((UInt_t) param == fNonConstMap.size()) {
205 Double_t chi = GetChi2();
206 if (chi < fGlobMin) {
207 fGlobMin = chi;
208 fStateVectorMin = fStateVector;
209 }
210 return kTRUE;
211 }
213 Int_t pairID = fNonConstMap[param];
214 for (unsigned int i = 0; i < fAllowedValues[pairID].size(); i++) {
215 fStateVector[pairID] = i;
216 Bool_t val = LoopOverParameter(param + 1);
217 if (val == kFALSE) return kFALSE;
218 }
219 return kTRUE;
220 }
221
222 void Minimizer::EstimateErrors() {
223 for (int i = 0; i < GetNParams(); i++) {
224 if (IsFixed(i)) {
225 fQuantumFits[i] = fParamsMin[i];
226 fSmoothFits[i] = fParamsMin[i];
227 fErrors[i] = 0;
228 } else {
229 Double_t par, err, parQ;
230 EstimateError(i, par, parQ, err);
231
232 if (fParameters[i].IsDiscrete()) {
233 fQuantumFits[i] = fParamsMin[i];
234 fSmoothFits[i] = fParamsMin[i];
235 } else {
236 fQuantumFits[i] = parQ;
237 fSmoothFits[i] = par;
238 }
239
240 fErrors[i] = err;
241 }
242 }
243 }
244
245 void Minimizer::EstimateError(Int_t par, Double_t& min, Double_t& quantumMin, Double_t& error) {
246
247 Double_t x2 = fParamsMin[par];
248 Double_t step = fParameters[par].GetDParam();
249 Double_t x1 = x2 - step;
250 Double_t x3 = x2 + step;
251 Double_t y2 = fGlobMin / fNDF;
252 Double_t y3, y1 = 0;
253 if (fParameters[par].IsFixed()) {
254 min = x2;
255 quantumMin = x2;
256 error = 0;
257 return;
258 }
259
260 for (int i = 0; i < GetNParams(); i++)
261 fTempParams[i] = fParamsMin[i];
262 if (x1 >= fParameters[par].GetMapMin()) {
263 fTempParams[par] = x1;
264 y1 = (*fFunc)(fTempParams) / fNDF;
265 }
266 if (x3 > fParameters[par].GetMapMax()) {
267 y3 = y1;
268 } else {
269 fTempParams[par] = x3;
270 y3 = (*fFunc)(fTempParams) / fNDF;
271 }
272 if (x1 < fParameters[par].GetMapMin()) { y1 = y3; }
273 fTempParams[par] = fParamsMin[par];
274 Double_t a, b, c;
275 Hal::Std::FitParabola(x1, x2, x3, y1, y2, y3, a, b, c);
276 Double_t X1, X2;
277 Hal::Std::SolveParabola(a, b, c - 0.9 * fGlobMin / fNDF, X1, X2);
278 Double_t min_Fit = -b / (2.0 * a);
279
280
281 // std::cout << " pAR " << x1 << " " << x2 << " " << x3 << " " << y1 << " " << y2 << " " << y3 << " " << a << std::endl;
282
283 quantumMin = fParamsMin[par];
284 min = min_Fit;
285 error = TMath::Max(TMath::Abs(min_Fit - X1), TMath::Abs(min_Fit - X2));
286 // hessian-like method
287 error = TMath::Sqrt(1.0 / a) * y2;
288 return;
289 /*
290 // SolveEquation(x1, x2, x3, y1, y2, y3, par);
291 Double_t C = fDeriveratives[par][0];
292 Double_t B = fDeriveratives[par][1];
293 Double_t A = fDeriveratives[par][2];
294 Double_t min_x = -B / (2.0 * A);
295 Double_t min_chi = min_x * B + min_x * min_x * A + C;
296 min = min_x;
297 C = C - min_chi * 1.1;
298 // ten minimizer daje bledy o polowe mniejsze dlaczego?
299
300 * przykladowo w pewnej analizie robilem skok o 1, mama chi2 znajduje minimum
301 * na 4 ale jak to sie fituje parabola (jak robi to minimizer) to minimum jest
302 * na 4.15 w efekcie najmniejsza wartoscia z paraboli jest 4.15 z dwoma
303 * miejscami zerowymi np. 4 i 4.3 chi2 liczy odleglosc do min. w mapie i daje
304 * 4 ten kod liczy odleglosc do min. w ficie parabolicznym i daje 4.15 w
305 * efekcie rozrzut w przypadku 1 to 0.3 a w drugim to 0.15
306
307 Double_t delta = TMath::Sqrt(B * B - 4. * A * C);
308 Double_t xs = (-B + delta) / (2. * A);
309 error = TMath::Abs(xs - min_x);*/
310 }
311
312 Bool_t Minimizer::IsFixed(Int_t i) const { return fParameters[i].IsFixed(); }
313
314 Minimizer* Minimizer::Instance() {
315 if (fgInstance == NULL) fgInstance = new Minimizer();
316 // fgInstance->Reset();
317 return fgInstance;
318 }
319
320 void Minimizer::FinishFit() {
321 for (int i = 0; i < GetNParams(); i++) {
322 fParamsMin[i] = fAllowedValues[i][fStateVectorMin[i]];
323 }
324 EstimateErrors();
325 delete[] fTempParams;
326 fTempParams = nullptr;
327 }
328
329 void Minimizer::SetMinimizerType(TString opt) {
330 opt.ToLower();
331 if (opt == "scan") {
332 fMinimizeType = kScan;
333 } else if (opt == "ant") {
334 fMinimizeType = kAnt;
335 }
336 }
337
338 void Minimizer::MinimizeAnt() {
339 std::vector<std::vector<int>> directionIndex;
340 for (int i = 0; i < GetNFree(); i++) {
341 std::vector<int> poz;
342 for (int j = 0; j < GetNFree(); j++) {
343 poz.push_back(0);
344 }
345 poz[i] = 1;
346 directionIndex.push_back(poz);
347 }
348 for (int i = 0; i < GetNFree(); i++) {
349 std::vector<int> poz;
350 for (int j = 0; j < GetNFree(); j++) {
351 poz.push_back(0);
352 }
353 poz[i] = -1;
354 directionIndex.push_back(poz);
355 }
356 int direction = 0;
357 int freePars = GetNFree();
358 std::vector<int> tmpVec(freePars, 0);
359 Bool_t cont = kTRUE;
360
361 auto applyVector = [&]() {
362 for (int p = 0; p < freePars; p++) {
363 int constToPar = fNonConstMap[p];
364 tmpVec[constToPar] += directionIndex[direction][p];
365 if (tmpVec[constToPar] < 0 || tmpVec[constToPar] > fParameters[constToPar].GetNPoints()) { return kFALSE; }
366 }
367 /*
368 std::cout << "VEC ";
369 for (auto i : directionIndex[direction]) {
370 std::cout << i << " ";
371 }
372 std::cout << std::endl;
373 std::cout << "STATE ";
374 for (auto i : fStateVector) {
375 std::cout << i << " ";
376 }
377 std::cout << std::endl;*/
378 for (int p = 0; p < freePars; p++) {
379 fStateVector[fNonConstMap[p]] = tmpVec[p];
380 }
381
382 return kTRUE;
383 };
384 double old_chi = 1E+9;
385 int bad = 0;
386 while (cont) {
387 cont = applyVector();
388 if (cont) {
389 double chi = GetChi2();
390 if (chi > old_chi) { // larger chi, change direction
391 direction = (direction + 1) % freePars;
392 bad++;
393 }
394 if (chi < fGlobMin) { // smaller chi
395 fGlobMin = chi;
396 fStateVectorMin = fStateVector;
397 bad = 0;
398 old_chi = chi;
399 }
400 } else { // out of bounds, whatever increase bad attemps counter
401 direction = (direction + 1) % freePars;
402 bad++;
403 }
404 fNCalls++;
405 if (bad == freePars) { // to many bad attempts
406 // std::cout << "Cannot find better minimum" << std::endl;
407 break;
408 }
409 if (fNCalls > 1E+6) break; // to many call
410 }
411 }
412
413 void Minimizer::MinimizeScan() {
414 // set start parameters
415 for (auto& x : fStateVector) {
416 x = 0;
417 }
418 SetTempParams();
419 for (unsigned int i = 0; i < fParameters.size(); i++) {
420 fTempParams[i] = fParameters[i].GetMin();
421 }
422 LoopOverParameter(0);
423 for (int i = 0; i < GetNParams(); i++) {
424 fTempParams[i] = fParamsMin[i];
425 fSmoothFits[i] = fParamsMin[i];
426 fQuantumFits[i] = fParamsMin[i];
427 }
428 }
429
430 void Minimizer::MinimizeNelderMead() {
431 InitFit();
432 // meat
433 Hal::Array_2<Double_t> vertices;
434 vertices.MakeBigger(fFreePars + 1, fFreePars);
435 Hal::Array_1<Double_t> chiMap;
436 chiMap.MakeBigger(fFreePars + 1);
437 Hal::Array_1<Double_t> centroid, futurePar, worstPar;
438 centroid.MakeBigger(fFreePars);
439 futurePar.MakeBigger(fFreePars);
440 worstPar.MakeBigger(fFreePars);
441 for (int i = 0; i < fFreePars; i++) {
442 for (int j = 0; j < fFreePars; j++) {
443 vertices[i][j] = fParameters[fNonConstMap[i]].GetMin();
444 }
445 }
446 for (int j = 0; j < fFreePars; j++) {
447 vertices[fFreePars][j] = fParameters[fNonConstMap[0]].GetMax();
448 }
449
450 auto copyPars = [&](const Hal::Array_1<Double_t>& /*val*/) {
451 for (int i = 0; i < fFreePars; i++) {
452 // fTempParams[fNonConstMap[i]] = fParameters[i].FindClosest(val[i]);
453 }
454 };
455
456 // calcualte first vertices
457 for (int i = 0; i < fFreePars + 1; i++) {
458 for (int j = 0; j < fFreePars; j++) {
459 copyPars(vertices[i]);
460 }
461 chiMap[i] = (*fFunc)(fTempParams);
462 }
463
464 Bool_t found = kTRUE;
465 while (found) {
466 int largest_chi, lowest_chi;
467 double max_chi = chiMap.FindMax(largest_chi);
468 double min_chi = chiMap.FindMin(lowest_chi);
469 // find worse vertice
470 double almost_max_chi = 0;
471 for (int i = 0; i < fFreePars + 1; i++) {
472 if (chiMap[i] > almost_max_chi && chiMap[i] != max_chi) { almost_max_chi = chiMap[i]; }
473 }
474 worstPar = vertices[largest_chi];
475
476 // double calc centroid
477
478 for (int i = 0; i < fFreePars; i++) {
479 centroid[i] = 0;
480 for (int j = 0; j < fFreePars + 1; j++) {
481 if (j == largest_chi) continue;
482 centroid += vertices[i];
483 }
484 centroid = centroid / double(fFreePars);
485 }
486 // do reflection
487 futurePar = centroid + (centroid - worstPar);
488 copyPars(futurePar);
489 double chi_r = (*fFunc)(fTempParams);
490
491 if (chi_r < almost_max_chi && chi_r > min_chi) {
492 vertices[largest_chi] = futurePar;
493 continue;
494 }
495 if (chi_r <= min_chi) { // expand
496 }
497
498 // double second_worst_chi = 0; //ntot used warning
499 for (int i = 0; i < fFreePars; i++) {
500 centroid[i] = 0;
501 for (int j = 0; j < fFreePars + 1; j++) {
502 if (j == largest_chi) continue;
503 centroid += vertices[i];
504 centroid += vertices[i];
505 }
506 centroid = centroid / double(fFreePars);
507 }
508 }
509
510 futurePar = centroid + (centroid - worstPar);
511 copyPars(futurePar);
512 // double chi_r = (*fFunc)(fTempParams); not used warning
513 // TO DO complete this https://codesachin.wordpress.com/2016/01/16/nelder-mead-optimization/
514
515 EstimateErrors();
516 FinishFit();
517 }
518
519 void Minimizer::ChangeStateVector(std::vector<Int_t>& vec) {
520 std::vector<Int_t> temp = vec;
521 for (unsigned int i = 0; i < vec.size(); i++) {
522 vec[i + 1] = temp[i];
523 }
524 vec[0] = temp[temp.size() - 1];
525 }
526
527 bool Minimizer::SetParLimits(unsigned int ivar, double min, double max) {
528 if (min > max) {
529 double temp = max;
530 max = min;
531 min = temp;
532 }
533 fParameters[ivar].SetMin(min);
534 fParameters[ivar].SetMax(max);
535 return true;
536 }
537
538 void Minimizer::SetTempParams() {
539 for (int i = 0; i < GetNParams(); i++) {
540 fTempParams[i] = fAllowedValues[i][fStateVector[i]];
541 }
542 }
543
544 Double_t Minimizer::GetChi2() {
545 SetTempParams();
546 return (*fFunc)(fTempParams);
547 }
548
549} // namespace Hal
Hal::Array_1
Definition Array.h:19
Hal::Array_1::MakeBigger
void MakeBigger(Int_t new_dim)
Definition Array.cxx:5
Hal::Array_2
Definition Array.h:127
Hal::Array_2::MakeBigger
void MakeBigger(Int_t sizeA, Int_t sizeB)
Definition Array.cxx:197