19 #include "analysis/OrcaKinFit/NewtonFitterGSL.h"
25 #include "analysis/OrcaKinFit/BaseFitObject.h"
26 #include "analysis/OrcaKinFit/BaseHardConstraint.h"
27 #include "analysis/OrcaKinFit/BaseSoftConstraint.h"
28 #include "analysis/OrcaKinFit/BaseTracer.h"
29 #include <framework/logging/Logger.h>
31 #include <gsl/gsl_block.h>
32 #include <gsl/gsl_vector.h>
33 #include <gsl/gsl_matrix.h>
34 #include <gsl/gsl_permutation.h>
35 #include <gsl/gsl_linalg.h>
36 #include <gsl/gsl_blas.h>
37 #include <gsl/gsl_cdf.h>
41 static int nitdebug = 100;
42 static int nitcalc = 0;
43 static int nitsvd = 0;
50 namespace OrcaKinFit {
54 : npar(0), ncon(0), nsoft(0), nunm(0), ierr(0), nit(0), fitprob(0), chi2(0),
56 x(nullptr), xold(nullptr), xbest(nullptr), dx(nullptr), dxscal(nullptr), grad(nullptr), y(nullptr), yscal(nullptr),
57 perr(nullptr), v1(nullptr), v2(nullptr), Meval(nullptr),
58 M(nullptr), Mscal(nullptr), M1(nullptr), M2(nullptr), M3(nullptr), M4(nullptr), M5(nullptr), Mevec(nullptr),
59 CC(nullptr), CC1(nullptr), CCinv(nullptr), permM(nullptr), ws(nullptr),
60 wsdim(0), chi2best(0),
63 fvalbest(0), scale(0), scalebest(0), stepsize(0), stepbest(0),
73 if (x) gsl_vector_free(x);
74 if (xold) gsl_vector_free(xold);
75 if (xbest) gsl_vector_free(xbest);
76 if (dx) gsl_vector_free(dx);
77 if (dxscal) gsl_vector_free(dxscal);
78 if (grad) gsl_vector_free(grad);
79 if (y) gsl_vector_free(y);
80 if (yscal) gsl_vector_free(yscal);
81 if (perr) gsl_vector_free(perr);
82 if (v1) gsl_vector_free(v1);
83 if (v2) gsl_vector_free(v2);
84 if (Meval) gsl_vector_free(Meval);
85 if (M) gsl_matrix_free(M);
86 if (Mscal) gsl_matrix_free(Mscal);
87 if (M1) gsl_matrix_free(M1);
88 if (M2) gsl_matrix_free(M2);
89 if (M3) gsl_matrix_free(M3);
90 if (M4) gsl_matrix_free(M4);
91 if (M5) gsl_matrix_free(M5);
92 if (Mevec) gsl_matrix_free(Mevec);
93 if (CC) gsl_matrix_free(CC);
94 if (CC1) gsl_matrix_free(CC1);
95 if (CCinv) gsl_matrix_free(CCinv);
96 if (permM) gsl_permutation_free(permM);
97 if (ws) gsl_eigen_symmv_free(ws);
123 ws =
nullptr; wsdim = 0;
135 assert(x && (
unsigned int)x->size == idim);
136 assert(dx && (
unsigned int)dx->size == idim);
137 assert(y && (
unsigned int)y->size == idim);
138 assert(perr && (
unsigned int)perr->size == idim);
139 assert(v1 && (
unsigned int)v1->size == idim);
140 assert(v2 && (
unsigned int)v2->size == idim);
141 assert(Meval && (
unsigned int)Meval->size == idim);
142 assert(M && (
unsigned int)M->size1 == idim);
143 assert(M1 && (
unsigned int)M1->size1 == idim);
144 assert(Mevec && (
unsigned int)Mevec->size1 == idim);
145 assert(permM && (
unsigned int)permM->size == idim);
147 gsl_vector_set_zero(x);
148 gsl_vector_set_zero(y);
149 gsl_vector_set_all(perr, 1);
157 gsl_vector_memcpy(x, xold);
166 bool converged =
false;
172 B2INFO(
"Fit objects:\n");
173 for (
auto fo : fitobjects) {
175 B2INFO(fo->getName() <<
": " << *fo <<
", chi2=" << fo->getChi2());
176 for (
int ilocal = 0; ilocal < fo->getNPar(); ++ilocal) {
177 if (!fo->isParamFixed(ilocal)) {
178 int iglobal = fo->getGlobalParNum(ilocal);
179 B2INFO(
" par " << fo->getParamName(ilocal) <<
": local: " << ilocal <<
": global: " << iglobal
180 <<
" value=" << fo->getParam(ilocal) <<
" +- " << fo->getError(ilocal));
181 if (fo->isParamMeasured(ilocal))
182 B2INFO(
" measured: " << fo->getMParam(ilocal));
184 B2INFO(
" par " << fo->getParamName(ilocal) <<
": local: " << ilocal <<
" -- fixed -- "
185 <<
" value=" << fo->getParam(ilocal) <<
" +- " << fo->getError(ilocal));
190 B2INFO(
"constraints:\n");
191 for (
auto i = constraints.begin(); i != constraints.end(); ++i) {
194 B2INFO(i - constraints.begin() <<
" " << c->getName() <<
": " << c->getValue() <<
"+-" << c->getError());
195 int kglobal = c->getGlobalNum();
196 B2INFO(
" global number: " << kglobal);
198 B2INFO(
"soft constraints:\n");
199 for (
auto i = softconstraints.begin(); i != softconstraints.end(); ++i) {
202 B2INFO(i - softconstraints.begin() <<
" " << c->getName() <<
": " << c->getValue() <<
"+-" << c->getError()
203 <<
", chi2: " << c->getChi2());
211 if (
nit == 0 ||
nit < nitdebug) {
212 B2DEBUG(11,
"===================\nStarting iteration " <<
nit);
213 std::string printout =
"\n Fit objects:\n";
214 for (
auto fo : fitobjects) {
216 printout += fo->getName();
217 printout +=
", chi2=" ;
218 printout += fo->getChi2();
221 printout +=
"constraints:\n";
222 for (
auto c : constraints) {
224 printout += c->getName() ;
226 printout += c->getValue();
228 printout += c->getError() ;
231 printout +=
"soft constraints:\n";
232 for (
auto c : softconstraints) {
234 printout += c->getName() ;
236 printout += c->getValue() ;
238 printout += c->getError();
239 printout +=
", chi2: ";
240 printout += c->getChi2() ;
243 B2DEBUG(12, printout);
258 if (
nit == 0 ||
nit < nitdebug) {
259 B2DEBUG(13,
"After setting up equations: \n");
261 debug_print(Mscal,
"Mscal");
263 debug_print(yscal,
"yscal");
264 debug_print(perr,
"perr");
266 debug_print(xold,
"xold");
270 fvals[0] = 0.5 * pow(gsl_blas_dnrm2(yscal), 2);
278 B2INFO(
"NewtonFitterGSL::fit: calcDx: ifail=" << ifail);
286 debug_print(xbest,
"new parameters");
290 B2DEBUG(13,
"New fval: " << 0.5 * pow(gsl_blas_dnrm2(yscal), 2));
292 B2DEBUG(13,
"chi2: " << chi2old <<
" -> " << chi2new);
294 if (
nit == 0 ||
nit < nitdebug) {
295 B2DEBUG(13,
"After solving equations: \n");
296 debug_print(xbest,
"xbest");
302 if (
nit < nitdebug) {
303 B2DEBUG(11,
"old chi2: " << chi2old <<
", new chi2: " << chi2new <<
", diff=" << chi2old - chi2new);
308 converged = (abs(chi2new - chi2old) < 0.001 && fvalbest < 1E-3 &&
309 (fvalbest < 1E-6 || abs(fvals[0] - fvalbest) < 0.2 * fvalbest));
322 B2DEBUG(10,
"abs (chi2new - chi2old)=" << abs(chi2new - chi2old) <<
"\n"
323 <<
"fvalbest=" << fvalbest <<
"\n"
324 <<
"abs(fvals[0]-fvalbest)=" << abs(fvals[0] - fvalbest) <<
"\n");
328 if (tracer) tracer->
step(*
this);
331 }
while (!(converged ||
ierr));
342 for (
auto& fitobject : fitobjects) {
343 for (
int ilocal = 0; ilocal < fitobject->getNPar(); ++ilocal) {
344 int iglobal = fitobject->getGlobalParNum(ilocal);
345 for (
int jlocal = ilocal; jlocal < fitobject->getNPar(); ++jlocal) {
346 int jglobal = fitobject->getGlobalParNum(jlocal);
347 if (iglobal >= 0 && jglobal >= 0)
348 fitobject->setCov(ilocal, jlocal, gsl_matrix_get(CCinv, iglobal, jglobal));
355 std::string printout =
"\n ========= END =========\n";
356 printout +=
"Fit objects:\n";
357 for (
auto fo : fitobjects) {
359 printout += fo->getName() ;
361 printout +=
", chi2=" ;
362 printout += fo->getChi2() ;
365 printout +=
"constraints:\n";
366 for (
auto c : constraints) {
368 printout += c->getName();
370 printout += c->getValue();
372 printout += c->getError();
375 printout +=
"=============================================\n";
376 B2DEBUG(11, printout);
382 if (tracer) tracer->
finish(*
this);
385 B2DEBUG(10,
"NewtonFitterGSL::fit: converged=" << converged
386 <<
", nit=" <<
nit <<
", fitprob=" <<
fitprob);
402 for (
auto& fitobject : fitobjects) {
403 for (
int ilocal = 0; ilocal < fitobject->getNPar(); ++ilocal) {
404 if (!fitobject->isParamFixed(ilocal)) {
405 B2DEBUG(13,
"NewtonFitterGSL::initialize: parameter " << ilocal
406 <<
" of fitobject " << fitobject->getName()
407 <<
" gets global number " <<
npar);
408 fitobject->setGlobalParNum(ilocal,
npar);
410 if (!fitobject->isParamMeasured(ilocal)) ++
nunm;
416 ncon = constraints.size();
418 for (
int icon = 0; icon <
ncon; ++icon) {
421 B2DEBUG(13,
"NewtonFitterGSL::initialize: constraint " << c->getName()
422 <<
" gets global number " <<
npar + icon);
423 c->setGlobalNum(
npar + icon);
424 B2DEBUG(14,
"Constraint " << icon <<
" -> global " << c->getGlobalNum());
427 nsoft = softconstraints.size();
430 B2ERROR(
"NewtonFitterGSL::initialize: nunm=" <<
nunm <<
" > ncon+nsoft="
436 ini_gsl_vector(x, idim);
437 ini_gsl_vector(xold, idim);
438 ini_gsl_vector(xbest, idim);
439 ini_gsl_vector(dx, idim);
440 ini_gsl_vector(dxscal, idim);
441 ini_gsl_vector(grad, idim);
442 ini_gsl_vector(y, idim);
443 ini_gsl_vector(yscal, idim);
444 ini_gsl_vector(perr, idim);
445 ini_gsl_vector(v1, idim);
446 ini_gsl_vector(v2, idim);
447 ini_gsl_vector(Meval, idim);
449 ini_gsl_matrix(M, idim, idim);
450 ini_gsl_matrix(Mscal, idim, idim);
451 ini_gsl_matrix(M1, idim, idim);
452 ini_gsl_matrix(M2, idim, idim);
453 ini_gsl_matrix(M3, idim, idim);
454 ini_gsl_matrix(M4, idim, idim);
455 ini_gsl_matrix(M5, idim, idim);
456 ini_gsl_matrix(Mevec, idim, idim);
457 ini_gsl_matrix(CC, idim, idim);
458 ini_gsl_matrix(CC1, idim, idim);
459 ini_gsl_matrix(CCinv, idim, idim);
461 ini_gsl_permutation(permM, idim);
463 if (ws && wsdim != idim) {
464 gsl_eigen_symmv_free(ws);
467 if (ws ==
nullptr) ws = gsl_eigen_symmv_alloc(idim);
477 for (
auto fo : fitobjects) {
479 chi2 += fo->getChi2();
481 for (
auto bsc : softconstraints) {
483 chi2 += bsc->getChi2();
490 for (
int i = 0; i < idime; ++i) {
491 B2INFO(i <<
" [ " << Ma[idime * i + 0]);
492 for (
int j = 1; j < idime; ++j) B2INFO(
", " << Ma[idime * i + j]);
493 B2INFO(
"] [" << yo[i] <<
"]\n");
505 B2DEBUG(11,
"entering calcDx");
510 gsl_matrix_memcpy(M1, Mscal);
516 int result = gsl_linalg_LU_decomp(M1, permM, &signum);
517 B2DEBUG(11,
"calcDx: gsl_linalg_LU_decomp result=" << result);
519 ifail = gsl_linalg_LU_solve(M1, permM, yscal, dxscal);
520 B2DEBUG(11,
"calcDx: gsl_linalg_LU_solve result=" << ifail);
523 B2ERROR(
"NewtonFitter::calcDx: ifail from gsl_linalg_LU_solve=" << ifail);
527 stepsize = std::abs(gsl_vector_get(dxscal, gsl_blas_idamax(dxscal)));
530 gsl_vector_memcpy(dx, dxscal);
531 gsl_vector_mul(dx, perr);
535 gsl_vector_memcpy(xbest, xold);
540 if (scalebest < 0.01) {
541 B2DEBUG(11,
"NewtonFitter::calcDx: reverting to calcDxSVD\n");
555 for (
unsigned int i = 0; i < idim; ++i) assert(gsl_vector_get(perr, i) > 0);
559 gsl_matrix_memcpy(M1, Mscal);
560 B2DEBUG(13,
"NewtonFitterGSL::calcDxSVD: Calling gsl_eigen_symmv");
561 ierrc = gsl_eigen_symmv(M1, Meval, Mevec, ws);
562 B2DEBUG(13,
"NewtonFitterGSL::calcDxSVD: result of gsl_eigen_symmv: ");
564 B2ERROR(
"NewtonFitter::calcDxSVD: ierr=" << ierrc <<
"from gsl_eigen_symmv!\n");
567 ierrc = gsl_eigen_symmv_sort(Meval, Mevec, GSL_EIGEN_SORT_ABS_DESC);
569 B2ERROR(
"NewtonFitter::calcDxSVD: ierr=" << ierrc <<
"from gsl_eigen_symmv_sort!\n");
592 unsigned int ndim = 0;
593 while (ndim < idim && gsl_vector_get(Meval, ndim) != 0) ++ndim;
596 B2INFO(
"calcDxSVD: idim = " << idim <<
" > ndim = " << ndim);
601 gsl_blas_dgemv(CblasTrans, 1, Mevec, yscal, 0, v2);
604 for (
unsigned int i = 0; i < idim; ++i) {
605 if (
double e = gsl_vector_get(Meval, i)) gsl_vector_set(v2, i, gsl_vector_get(v2, i) / e);
606 else gsl_vector_set(v2, i, 0);
612 gsl_vector_view v2part = gsl_vector_subvector(v2, 0, ndim);
613 gsl_matrix_view Mevecpart = gsl_matrix_submatrix(Mevec, 0, 0, idim, ndim);
616 gsl_blas_dgemv(CblasNoTrans, 1, &Mevecpart.matrix, &v2part.vector, 0, dxscal);
622 stepsize = std::abs(gsl_vector_get(dxscal, gsl_blas_idamax(dxscal)));
625 gsl_vector_memcpy(dx, dxscal);
626 gsl_vector_mul(dx, perr);
629 B2INFO(
"calcDxSVD: Optimizing scale for ndim=" << ndim);
630 debug_print(dxscal,
"dxscal");
637 if (scalebest < 0.01 || ndim < idim - 1) {
638 B2DEBUG(11,
"ndim=" << ndim <<
", scalebest=" << scalebest);
640 }
while (ndim > 0 && scalebest < 0.01);
646 void NewtonFitterGSL::ini_gsl_permutation(gsl_permutation*& p,
unsigned int size)
649 if (p->size != size) {
650 gsl_permutation_free(p);
651 if (size > 0) p = gsl_permutation_alloc(size);
654 }
else if (size > 0) p = gsl_permutation_alloc(size);
657 void NewtonFitterGSL::ini_gsl_vector(gsl_vector*& v,
unsigned int size)
661 if (v->size != size) {
663 if (size > 0) v = gsl_vector_alloc(size);
666 }
else if (size > 0) v = gsl_vector_alloc(size);
669 void NewtonFitterGSL::ini_gsl_matrix(gsl_matrix*& m,
unsigned int size1,
unsigned int size2)
672 if (m->size1 != size1 || m->size2 != size2) {
674 if (size1 * size2 > 0) m = gsl_matrix_alloc(size1, size2);
677 }
else if (size1 * size2 > 0) m = gsl_matrix_alloc(size1, size2);
680 void NewtonFitterGSL::debug_print(gsl_matrix* m,
const char* name)
682 for (
unsigned int i = 0; i < m->size1; ++i)
683 for (
unsigned int j = 0; j < m->size2; ++j)
684 if (gsl_matrix_get(m, i, j) != 0)
685 B2INFO(name <<
"[" << i <<
"][" << j <<
"]=" << gsl_matrix_get(m, i, j));
688 void NewtonFitterGSL::debug_print(gsl_vector* v,
const char* name)
690 for (
unsigned int i = 0; i < v->size; ++i)
691 if (gsl_vector_get(v, i) != 0)
692 B2INFO(name <<
"[" << i <<
"]=" << gsl_vector_get(v, i));
700 bool NewtonFitterGSL::updateParams(gsl_vector* xnew)
703 assert(xnew->size == idim);
704 bool significant =
false;
705 for (
auto i = fitobjects.begin(); i != fitobjects.end(); ++i) {
706 BaseFitObject* fo = *i;
708 bool s = fo->updateParams(xnew->block->data, xnew->size);
710 if (
nit < nitdebug && s) {
711 B2DEBUG(11,
"Significant update for FO " << i - fitobjects.begin() <<
" ("
712 << fo->getName() <<
")\n");
719 void NewtonFitterGSL::fillxold()
722 assert(xold->size == idim);
723 for (
auto fo : fitobjects) {
725 for (
int ilocal = 0; ilocal < fo->getNPar(); ++ilocal) {
726 if (!fo->isParamFixed(ilocal)) {
727 int iglobal = fo->getGlobalParNum(ilocal);
728 assert(iglobal >= 0 && iglobal <
npar);
729 gsl_vector_set(xold, iglobal, fo->getParam(ilocal));
733 for (
auto c : constraints) {
735 int iglobal = c->getGlobalNum();
736 assert(iglobal >= 0 && iglobal < (
int)idim);
737 gsl_vector_set(xold, iglobal, gsl_vector_get(x, iglobal));
741 void NewtonFitterGSL::fillperr()
744 assert(perr->size == idim);
745 gsl_vector_set_all(perr, 1);
746 for (
auto fo : fitobjects) {
748 for (
int ilocal = 0; ilocal < fo->getNPar(); ++ilocal) {
749 if (!fo->isParamFixed(ilocal)) {
750 int iglobal = fo->getGlobalParNum(ilocal);
751 assert(iglobal >= 0 && iglobal <
npar);
752 double e = std::abs(fo->getError(ilocal));
753 gsl_vector_set(perr, iglobal, e ? e : 1);
757 for (
auto c : constraints) {
759 int iglobal = c->getGlobalNum();
760 assert(iglobal >= 0 && iglobal < (
int)idim);
761 double e = c->getError();
762 gsl_vector_set(perr, iglobal, e ? 1 / e : 1);
766 int NewtonFitterGSL::calcM(
bool errorpropagation)
769 assert(M->size1 == idim && M->size2 == idim);
771 gsl_matrix_set_zero(M);
774 for (
auto fo : fitobjects) {
776 fo->addToGlobalChi2DerMatrix(M->block->data, M->tda);
779 B2INFO(
"After adding covariances from fit objects:\n");
786 for (
auto c : constraints) {
788 int kglobal = c->getGlobalNum();
789 assert(kglobal >= 0 && kglobal < (
int)idim);
790 c->add1stDerivativesToMatrix(M->block->data, M->tda);
792 B2INFO(
"After adding first derivatives of constraint " << c->getName());
794 B2INFO(
"errorpropagation = " << errorpropagation);
798 if (!errorpropagation) c->add2ndDerivativesToMatrix(M->block->data, M->tda, gsl_vector_get(x, kglobal));
801 B2INFO(
"After adding derivatives of constraints::\n");
803 B2INFO(
"===========================================::\n");
809 for (
auto bsc : softconstraints) {
811 bsc->add2ndDerivativesToMatrix(M->block->data, M->tda);
814 B2INFO(
"After adding soft constraints::\n");
816 B2INFO(
"===========================================::\n");
820 for (
unsigned int i = 0; i < idim; ++i)
821 for (
unsigned int j = 0; j < idim; ++j)
822 gsl_matrix_set(Mscal, i, j,
823 gsl_vector_get(perr, i)*gsl_vector_get(perr, j)*gsl_matrix_get(M, i, j));
829 int NewtonFitterGSL::calcy()
832 assert(y->size == idim);
833 gsl_vector_set_zero(y);
835 for (
auto fo : fitobjects) {
837 fo->addToGlobalChi2DerVector(y->block->data, y->size);
842 for (
auto c : constraints) {
844 int kglobal = c->getGlobalNum();
845 assert(kglobal >= 0 && kglobal < (
int)idim);
846 c->addToGlobalChi2DerVector(y->block->data, y->size, gsl_vector_get(x, kglobal));
847 gsl_vector_set(y, kglobal, c->getValue());
852 for (
auto bsc : softconstraints) {
854 bsc->addToGlobalChi2DerVector(y->block->data, y->size);
856 gsl_vector_memcpy(yscal, y);
857 gsl_vector_mul(yscal, perr);
861 int NewtonFitterGSL::optimizeScale()
866 B2INFO(
"NewtonFitterGSL::optimizeScale");
867 debug_print(xold,
"xold");
868 debug_print(yscal,
"yscal");
869 debug_print(dx,
"dx");
870 debug_print(dxscal,
"dxscal");
873 fvals[0] = 0.5 * pow(gsl_blas_dnrm2(yscal), 2);
874 B2DEBUG(11,
"NewtonFitterGSL::optimizeScale: fvals[0] = " << fvals[0]);
881 gsl_blas_dgemv(CblasNoTrans, 1, Mscal, yscal, 0, grad);
883 debug_print(grad,
"grad");
891 static const double ALF = 1E-4;
893 stepsize = std::abs(gsl_vector_get(dxscal, gsl_blas_idamax(dxscal)));
894 static const double maxstepsize = 5;
895 double scalefactor = maxstepsize / stepsize;
896 if (stepsize > maxstepsize) {
897 gsl_vector_scale(dxscal, maxstepsize / stepsize);
898 B2DEBUG(12,
"NewtonFitterGSL::optimizeScale: Rescaling dxscal by factor " << scalefactor);
899 stepsize = std::abs(gsl_vector_get(dxscal, gsl_blas_idamax(dxscal)));
901 debug_print(dxscal,
"dxscal");
906 gsl_blas_ddot(dxscal, grad, &slope);
908 B2DEBUG(12,
"NewtonFitterGSL::optimizeScale: slope=" << slope
909 <<
", 2*fvals[0]*factor=" << 2 * fvals[0]*scalefactor);
912 double tmpscale, scaleold = 1;
915 gsl_vector_memcpy(x, xold);
918 debug_print(x,
"x(1)");
920 gsl_blas_daxpy(-scale, dx, x);
922 debug_print(x,
"x(2)");
928 debug_print(x,
"x(3)");
929 debug_print(yscal,
"yscal");
932 scalevals[nite] = scale;
933 fvals[nite] = 0.5 * pow(gsl_blas_dnrm2(yscal), 2);
941 if ((fvals[nite] < fvalbest)) {
942 B2DEBUG(13,
"new best value: "
943 <<
" scale " << scalevals[nite] <<
" -> |y|^2 = " << fvals[nite]
944 <<
", chi2=" << chi2new <<
", old best chi2: " << chi2best);
945 gsl_vector_memcpy(xbest, x);
947 fvalbest = fvals[nite];
949 stepbest = scale * stepsize;
952 if (fvals[nite] < fvals[0] + ALF * scale * slope)
break;
954 tmpscale = -slope / (2 * (fvals[nite] - fvals[0] - slope));
955 B2DEBUG(13,
"quadratic estimate for best scale: " << tmpscale);
957 double rhs1 = fvals[nite] - fvals[0] - scale * slope;
958 double rhs2 = fvals[nite - 1] - fvals[0] - scaleold * slope;
959 double a = (rhs1 / (scale * scale) - rhs2 / (scaleold * scaleold)) / (scale - scaleold);
960 double b = (-scaleold * rhs1 / (scale * scale) + scale * rhs2 / (scaleold * scaleold)) / (scale - scaleold);
961 if (a == 0) tmpscale = -slope / (2 * b);
963 double disc = b * b - 3 * a * slope;
964 if (disc < 0) tmpscale = 0.5 * scale;
965 else if (b <= 0) tmpscale = (-b + sqrt(disc)) / (3 * a);
966 else tmpscale = -slope / (b + sqrt(disc));
968 B2DEBUG(13,
"cubic estimate for best scale: " << tmpscale);
969 if (tmpscale > 0.5 * scale) tmpscale = 0.5 * scale;
972 scale = (tmpscale < 0.1 * scale) ? 0.1 * scale : tmpscale;
973 B2DEBUG(11,
"New scale: " << scale);
975 }
while (nite < NITMAX && scale > 0.0001);
978 for (
int it = 0; it <= nite; ++it) {
979 B2INFO(
" scale " << scalevals[it] <<
" -> |y|^2 = " << fvals[it]
980 <<
" should be " << fvals[0] + ALF * scale * slope);
983 return chi2best < chi2old;
986 int NewtonFitterGSL::invertM()
988 gsl_matrix_memcpy(M1, M);
994 int result = gsl_linalg_LU_decomp(M1, permM, &signum);
995 B2DEBUG(11,
"invertM: gsl_linalg_LU_decomp result=" << result);
997 ifail = gsl_linalg_LU_invert(M1, permM, M);
998 B2DEBUG(11,
"invertM: gsl_linalg_LU_solve result=" << ifail);
1001 B2ERROR(
"NewtonFitter::invert: ifail from gsl_linalg_LU_invert=" << ifail);
1012 void NewtonFitterGSL::calcCovMatrix()
1023 gsl_matrix_set_zero(M1);
1024 gsl_matrix_set_zero(M2);
1026 for (
auto fo : fitobjects) {
1028 fo->addToGlobalChi2DerMatrix(M1->block->data, M1->tda);
1029 fo->addToGlobCov(M2->block->data, M2->tda);
1033 gsl_matrix_scale(M1, -1);
1038 gsl_matrix_view dydeta = gsl_matrix_submatrix(M1, 0, 0, idim,
npar);
1039 gsl_matrix_view Cov_eta = gsl_matrix_submatrix(M2, 0, 0,
npar,
npar);
1042 B2INFO(
"NewtonFitterGSL::calcCovMatrix\n");
1043 debug_print(&dydeta.matrix,
"dydeta");
1044 debug_print(&Cov_eta.matrix,
"Cov_eta"); \
1053 debug_print(M,
"M");
1060 int result = gsl_linalg_LU_decomp(M, permM, &signum);
1063 B2INFO(
"invertM: gsl_linalg_LU_decomp result=" << result);
1064 debug_print(M,
"M_LU");
1068 int ifail = gsl_linalg_LU_invert(M, permM, M3);
1071 B2INFO(
"invertM: gsl_linalg_LU_invert ifail=" << ifail);
1072 debug_print(M3,
"Minv");
1076 gsl_matrix_set_zero(M4);
1077 gsl_matrix_view dadeta = gsl_matrix_submatrix(M4, 0, 0, idim,
npar);
1080 debug_print(&dadeta.matrix,
"dadeta");
1084 gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1, M3, &dydeta.matrix, 0, &dadeta.matrix);
1090 gsl_matrix_view M3part = gsl_matrix_submatrix(M3, 0, 0,
npar, idim);
1091 gsl_blas_dgemm(CblasNoTrans, CblasTrans, 1, &Cov_eta.matrix, &dadeta.matrix, 0, &M3part.matrix);
1093 gsl_matrix_set_zero(M5);
1094 gsl_matrix_view Cov_a = gsl_matrix_submatrix(M5, 0, 0,
npar,
npar);
1095 gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1, &dadeta.matrix, &M3part.matrix, 0, M5);
1096 gsl_matrix_memcpy(CCinv, M5);
1099 debug_print(&Cov_a.matrix,
"Cov_a");
1100 debug_print(CCinv,
"full Cov from err prop");
1101 debug_print(M1,
"uncorr Cov from err prop");
1111 for (
int i = 0; i <
covDim; ++i) {
1112 for (
int j = 0; j <
covDim; ++j) {
1113 cov[i *
covDim + j] = gsl_matrix_get(&Cov_a.matrix, i, j);
1120 double NewtonFitterGSL::meritFunction(
double mu)
1126 for (
auto c : constraints) {
1128 int kglobal = c->getGlobalNum();
1129 assert(kglobal >= 0 && kglobal < (
int)idim);
1130 result += mu * std::fabs(c->getValue());
1138 double NewtonFitterGSL::meritFunctionDeriv()