26 const std::vector<double> &pdir,
double tlr,
unsigned int maxcalls)
const
34 unsigned int npar = par.size();
35 unsigned int nfcn = 0;
38 double mgr_tlr = 0.5 * tlr;
45 double tlf = tlr * up;
47 unsigned int maxitr = 15;
49 double aminsv =
fFval;
50 double aim = aminsv + up;
54 std::vector<double> alsb(3, 0.), flsb(3, 0.);
56 MnPrint print(
"MnFunctionCross");
58 print.
Debug([&](std::ostream &os) {
59 for (
unsigned int i = 0; i < par.size(); ++i)
60 os <<
"Parameter " << par[i] <<
" value " << pmid[i] <<
" dir " << pdir[i] <<
" function min = " << aminsv
61 <<
" contur value aim = (fmin + up) = " << aim;
67 for (
unsigned int i = 0; i < par.size(); i++) {
68 unsigned int kex = par[i];
70 double zmid = pmid[i];
71 double zdir = pdir[i];
81 aulim = std::min(aulim, (zlim - zmid) / zdir);
90 aulim = std::min(aulim, (zlim - zmid) / zdir);
95 print.
Debug(
"Largest allowed aulim", aulim);
98 if (limset && npar == 1) {
99 print.
Warn(
"Parameter is at limit", pmid[0],
"delta", pdir[0]);
103 if (aulim < aopt + tla)
108 print.
Info([&](std::ostream &os) {
109 os <<
"Run Migrad with fixed parameters:";
110 for (
unsigned i = 0; i < npar; ++i)
111 os <<
"\n Pos " << par[i] <<
": " <<
fState.
Name(par[i]) <<
" = " << pmid[i];
114 for (
unsigned int i = 0; i < npar; i++)
127 print.
Warn(
"New minimum found while scanning parameter", par.front(),
"new value =", min0.
Fval(),
128 "old value =",
fFval);
135 if (limset ==
true && min0.
Fval() < aim)
140 flsb[0] = min0.
Fval();
141 flsb[0] = std::max(flsb[0], aminsv + 0.1 * up);
142 aopt =
std::sqrt(up / (flsb[0] - aminsv)) - 1.;
156 print.
Debug(
"flsb[0]", flsb[0],
"aopt", aopt);
158 print.
Info([&](std::ostream &os) {
159 os <<
"Run Migrad again (2nd) with fixed parameters:";
160 for (
unsigned i = 0; i < npar; ++i)
161 os <<
"\n Pos " << par[i] <<
": " <<
fState.
Name(par[i]) <<
" = " << pmid[i] + (aopt)*pdir[i];
164 for (
unsigned int i = 0; i < npar; i++)
165 migrad.
SetValue(par[i], pmid[i] + (aopt)*pdir[i]);
178 if (limset ==
true && min1.
Fval() < aim)
183 flsb[1] = min1.
Fval();
184 double dfda = (flsb[1] - flsb[0]) / (alsb[1] - alsb[0]);
186 print.
Debug(
"aopt", aopt,
"min1Val", flsb[1],
"dfda", dfda);
191 print.
Debug(
"dfda < 0 - iterate from", ipt,
"to max of", maxitr);
194 unsigned int maxlk = maxitr - ipt;
195 for (
unsigned int it = 0; it < maxlk; it++) {
199 aopt = alsb[0] + 0.2 * (it + 1);
206 print.
Info([&](std::ostream &os) {
207 os <<
"Run Migrad again (iteration " << it <<
" ) :";
208 for (
unsigned i = 0; i < npar; ++i)
209 os <<
"\n parameter " << par[i] <<
" (" <<
fState.
Name(par[i]) <<
") fixed to "
210 << pmid[i] + (aopt)*pdir[i];
213 for (
unsigned int i = 0; i < npar; i++)
214 migrad.
SetValue(par[i], pmid[i] + (aopt)*pdir[i]);
216 min1 = migrad(maxcalls, mgr_tlr);
227 if (limset ==
true && min1.
Fval() < aim)
231 flsb[1] = min1.
Fval();
232 dfda = (flsb[1] - flsb[0]) / (alsb[1] - alsb[0]);
235 print.
Debug(
"aopt", aopt,
"min1Val", flsb[1],
"dfda", dfda);
248 aopt = alsb[1] + (aim - flsb[1]) / dfda;
250 print.
Debug(
"dfda > 0 : aopt", aopt);
257 if (adist < tla && fdist < tlf)
261 double bmin = std::min(alsb[0], alsb[1]) - 1.;
264 double bmax = std::max(alsb[0], alsb[1]) + 1.;
274 print.
Info([&](std::ostream &os) {
275 os <<
"Run Migrad again (3rd) with fixed parameters:";
276 for (
unsigned i = 0; i < npar; ++i)
277 os <<
"\n Pos " << par[i] <<
": " <<
fState.
Name(par[i]) <<
" = " << pmid[i] + (aopt)*pdir[i];
280 for (
unsigned int i = 0; i < npar; i++)
281 migrad.
SetValue(par[i], pmid[i] + (aopt)*pdir[i]);
294 if (limset ==
true && min2.
Fval() < aim)
299 flsb[2] = min2.
Fval();
303 double ecarmn =
std::fabs(flsb[2] - aim);
305 unsigned int ibest = 2;
306 unsigned int iworst = 0;
307 unsigned int noless = 0;
309 for (
unsigned int i = 0; i < 3; i++) {
311 if (ecart > ecarmx) {
315 if (ecart < ecarmn) {
323 print.
Debug(
"have three points : noless < aim; noless", noless,
"ibest", ibest,
"iworst", iworst);
328 if (noless == 1 || noless == 2)
331 if (noless == 0 && ibest != 2)
335 if (noless == 3 && ibest != 2) {
339 print.
Debug(
"All three points below - look again fir positive slope");
345 flsb[iworst] = flsb[2];
346 alsb[iworst] = alsb[2];
347 dfda = (flsb[1] - flsb[0]) / (alsb[1] - alsb[0]);
349 print.
Debug(
"New straight line using point 1-2; dfda", dfda);
363 print.
Debug(
"Parabola fit: iteration", ipt);
365 double coeff1 = parbol.
C();
366 double coeff2 = parbol.
B();
367 double coeff3 = parbol.
A();
368 double determ = coeff2 * coeff2 - 4. * coeff3 * (coeff1 - aim);
370 print.
Debug(
"Parabola fit: a =", coeff3,
"b =", coeff2,
"c =", coeff1,
"determ =", determ);
373 if (determ < prec.
Eps())
376 double x1 = (-coeff2 + rt) / (2. * coeff3);
377 double x2 = (-coeff2 - rt) / (2. * coeff3);
378 double s1 = coeff2 + 2. *
x1 * coeff3;
379 double s2 = coeff2 + 2. *
x2 * coeff3;
381 print.
Debug(
"Parabola fit: x1",
x1,
"x2",
x2,
"s1",
s1,
"s2", s2);
384 print.
Warn(
"Problem 1");
394 print.
Debug(
"Parabola fit: aopt", aopt,
"slope", slope);
401 print.
Debug(
"Delta(aopt)",
std::fabs(aopt - alsb[ibest]),
"tla", tla,
"Delta(F)",
std::fabs(flsb[ibest] - aim),
412 unsigned int ileft = 3;
413 unsigned int iright = 3;
414 unsigned int iout = 3;
418 for (
unsigned int i = 0; i < 3; i++) {
420 if (ecart < ecarmn) {
429 else if (flsb[i] > flsb[iright])
435 }
else if (ileft == 3)
437 else if (flsb[i] < flsb[ileft])
445 print.
Debug(
"ileft", ileft,
"iright", iright,
"iout", iout,
"ibest", ibest);
449 if (ecarmx > 10. *
std::fabs(flsb[iout] - aim))
450 aopt = 0.5 * (aopt + 0.5 * (alsb[iright] + alsb[ileft]));
453 double smalla = 0.1 * tla;
454 if (slope * smalla > tlf)
455 smalla = tlf / slope;
456 double aleft = alsb[ileft] + smalla;
457 double aright = alsb[iright] - smalla;
465 aopt = 0.5 * (aleft + aright);
475 print.
Info([&](std::ostream &os) {
476 os <<
"Run Migrad again at new point (#iter = " << ipt <<
" ):";
477 for (
unsigned i = 0; i < npar; ++i)
478 os <<
"\n\t - parameter " << i <<
" fixed to " << pmid[i] + (aopt)*pdir[i];
481 for (
unsigned int i = 0; i < npar; i++)
482 migrad.
SetValue(par[i], pmid[i] + (aopt)*pdir[i]);
484 min2 = migrad(maxcalls, mgr_tlr);
495 if (limset ==
true && min2.
Fval() < aim)
501 flsb[iout] = min2.
Fval();
503 }
while (ipt < maxitr);