53 template<
class GradFunc = IGradModelFunction>
99 deriv = 0.5 * ( f2 -
f1 )/
h;
102 deriv = (
f1 - f0 )/
h;
127 for (
unsigned int k = 0; k <
fN; ++k) {
133 void Gradient(
const double *
x,
const double *
p,
double f0,
double *
g) {
136 for (
unsigned int k = 0; k <
fN; ++k) {
146 g[k] = 0.5 * ( f2 -
f1 )/
h;
149 g[k] = (
f1 - f0 )/
h;
162 mutable std::vector<double>
fVec;
169 if (rval > - std::numeric_limits<double>::max() && rval < std::numeric_limits<double>::max() )
173 return -std::numeric_limits<double>::max();
176 return + std::numeric_limits<double>::max();
183 if (rval > - std::numeric_limits<double>::max() && rval < std::numeric_limits<double>::max() )
187 rval = -std::numeric_limits<double>::max();
192 rval = + std::numeric_limits<double>::max();
201 template <
class GFunc>
203 const double *
x1,
const double *
x2,
const double *
p,
double *
g) {
209 unsigned int npar = gfunc.NPar();
210 for (
unsigned int k = 0; k < npar; ++k ) {
212 g[k] = igDerEval(
x1,
x2 );
234 unsigned int n =
data.Size();
246 bool useBinIntegral = fitOpt.
fIntegral &&
data.HasBinEdges();
249 bool isWeighted =
data.IsWeighted();
252 std::cout <<
"\n\nFit data size = " <<
n << std::endl;
253 std::cout <<
"evaluate chi2 using function " << &func <<
" " <<
p << std::endl;
254 std::cout <<
"use empty bins " << fitOpt.
fUseEmpty << std::endl;
255 std::cout <<
"use integral " << fitOpt.
fIntegral << std::endl;
256 std::cout <<
"use all error=1 " << fitOpt.
fErrors1 << std::endl;
257 if (isWeighted) std::cout <<
"Weighted data set - sumw = " <<
data.SumOfContent() <<
" sumw2 = " <<
data.SumOfError2() << std::endl;
270 double maxResValue = std::numeric_limits<double>::max() /
n;
271 double wrefVolume = 1.0;
278 auto mapFunction = [&](
const unsigned i){
283 const auto x1 =
data.GetCoordComponent(i, 0);
284 const auto y =
data.Value(i);
285 auto invError =
data.InvError(i);
289 const double *
x =
nullptr;
290 std::vector<double> xc;
291 double binVolume = 1.0;
293 unsigned int ndim =
data.NDim();
294 xc.resize(
data.NDim());
295 for (
unsigned int j = 0; j < ndim; ++j) {
296 double xx = *
data.GetCoordComponent(i, j);
297 double x2 =
data.GetBinUpEdgeComponent(i, j);
299 xc[j] = 0.5*(
x2 + xx);
303 binVolume *= wrefVolume;
304 }
else if(
data.NDim() > 1) {
306 xc.resize(
data.NDim());
308 for (
unsigned int j = 1; j <
data.NDim(); ++j)
309 xc[j] = *
data.GetCoordComponent(i, j);
316 if (!useBinIntegral) {
320 fval = func (
x,
p );
326 std::vector<double>
x2(
data.NDim());
327 data.GetBinUpEdgeCoordinates(i,
x2.data());
328 fval = igEval(
x,
x2.data());
331 if (useBinVolume) fval *= binVolume;
335 double invWeight = 1.0;
342 invWeight =
data.SumOfContent()/
data.SumOfError2();
348 double invError2 = (fval > 0) ? invWeight / fval : 0.0;
355 std::cout <<
x[0] <<
" " <<
y <<
" " << 1./invError <<
" params : ";
356 for (
unsigned int ipar = 0; ipar < func.
NPar(); ++ipar)
357 std::cout <<
p[ipar] <<
"\t";
358 std::cout <<
"\tfval = " << fval <<
" bin volume " << binVolume <<
" ref " << wrefVolume << std::endl;
364 double tmp = (
y -fval )* invError;
365 double resval = tmp * tmp;
369 if ( resval < maxResValue )
380 auto redFunction = [](
const std::vector<double> & objs){
381 return std::accumulate(objs.begin(), objs.end(),
double{});
388 Warning(
"FitUtil::EvaluateChi2",
"Multithread execution policy requires IMT, which is disabled. Changing "
389 "to ROOT::EExecutionPolicy::kSequential.");
396 for (
unsigned int i=0; i<
n; ++i) {
397 res += mapFunction(i);
409 Error(
"FitUtil::EvaluateChi2",
"Execution policy unknown. Avalaible choices:\n ROOT::EExecutionPolicy::kSequential (default)\n ROOT::EExecutionPolicy::kMultiThread (requires IMT)\n");
428 unsigned int n =
data.Size();
431 std::cout <<
"\n\nFit data size = " <<
n << std::endl;
432 std::cout <<
"evaluate effective chi2 using function " << &func <<
" " <<
p << std::endl;
435 assert(
data.HaveCoordErrors() ||
data.HaveAsymErrors());
443 unsigned int ndim = func.
NDim();
448 double maxResValue = std::numeric_limits<double>::max() /
n;
452 for (
unsigned int i = 0; i <
n; ++ i) {
456 const double *
x =
data.GetPoint(i,
y);
458 double fval = func(
x,
p );
460 double delta_y_func =
y - fval;
464 const double *
ex = 0;
465 if (!
data.HaveAsymErrors() )
468 double eylow, eyhigh = 0;
469 ex =
data.GetPointError(i, eylow, eyhigh);
470 if ( delta_y_func < 0)
478 while ( j < ndim &&
ex[j] == 0.) { j++; }
485 double kPrecision = 1.E-8;
486 for (
unsigned int icoord = 0; icoord < ndim; ++icoord) {
488 if (
ex[icoord] > 0) {
492 double x0=
x[icoord];
493 double h = std::max( kEps*
std::abs(
ex[icoord]), 8.0*kPrecision*(
std::abs(x0) + kPrecision) );
495 double edx =
ex[icoord] * deriv;
498 std::cout <<
"error for coord " << icoord <<
" = " <<
ex[icoord] <<
" deriv " << deriv << std::endl;
503 double w2 = (e2 > 0) ? 1.0/e2 : 0;
504 double resval = w2 * (
y - fval ) * (
y - fval);
507 std::cout <<
x[0] <<
" " <<
y <<
" ex " <<
ex[0] <<
" ey " <<
ey <<
" params : ";
508 for (
unsigned int ipar = 0; ipar < func.
NPar(); ++ipar)
509 std::cout <<
p[ipar] <<
"\t";
510 std::cout <<
"\tfval = " << fval <<
"\tresval = " << resval << std::endl;
515 if ( resval < maxResValue )
528 std::cout <<
"chi2 = " << chi2 <<
" n = " << nPoints << std::endl;
545 MATH_ERROR_MSG(
"FitUtil::EvaluateChi2Residual",
"Error on the coordinates are not used in calculating Chi2 residual");
552 double y, invError = 0;
555 bool useBinIntegral = fitOpt.
fIntegral &&
data.HasBinEdges();
559 const double *
x1 =
data.GetPoint(i,
y, invError);
563 unsigned int ndim =
data.NDim();
564 double binVolume = 1.0;
565 const double *
x2 = 0;
566 if (useBinVolume || useBinIntegral)
x2 =
data.BinUpEdge(i);
571 xc =
new double[ndim];
572 for (
unsigned int j = 0; j < ndim; ++j) {
574 xc[j] = 0.5*(
x2[j]+
x1[j]);
577 binVolume /=
data.RefVolume();
580 const double *
x = (useBinVolume) ? xc :
x1;
582 if (!useBinIntegral) {
583 fval = func (
x,
p );
588 fval = igEval(
x1,
x2) ;
591 if (useBinVolume) fval *= binVolume;
602 double invError2 = (fval > 0) ? 1.0 / fval : 0.0;
607 double resval = (
y -fval )* invError;
615 unsigned int npar = func.
NPar();
621 if (!
h ) useFullHessian =
false;
622 if (useFullHessian && (!gfunc || useBinIntegral || (gfunc && !gfunc->
HasParameterHessian())))
623 return std::numeric_limits<double>::quiet_NaN();
627 if (!useBinIntegral ) {
629 if (useFullHessian) {
640 if (!useBinIntegral ) {
641 gc.ParameterGradient(
x,
p, fval,
g);
648 for (
unsigned int k = 0; k < npar; ++k) {
650 if (useBinVolume)
g[k] *= binVolume;
652 for (
unsigned int l = 0;
l <= k;
l++) {
653 unsigned int idx =
l + k * (k + 1) / 2;
654 if (useFullHessian) {
655 h[idx] *= 2.* resval * (-invError);
656 if (useBinVolume)
h[idx] *= binVolume;
662 h[idx] += 2. *
g[k]*
g[
l];
668 if (useBinVolume)
delete [] xc;
683 if (
data.HaveCoordErrors()) {
685 "Error on the coordinates are not used in calculating Chi2 gradient");
690 assert(fg !=
nullptr);
695 std::cout <<
"\n\nFit data size = " <<
n << std::endl;
696 std::cout <<
"evaluate chi2 using function gradient " << &func <<
" " <<
p << std::endl;
699 const DataOptions &fitOpt =
data.Opt();
700 bool useBinIntegral = fitOpt.fIntegral &&
data.HasBinEdges();
701 bool useBinVolume = (fitOpt.fBinVolume &&
data.HasBinEdges());
703 double wrefVolume = 1.0;
705 if (fitOpt.fNormBinVolume) wrefVolume /=
data.RefVolume();
713 IntegralEvaluator<> igEval(func,
p, useBinIntegral,igType);
715 unsigned int npar = func.NPar();
716 unsigned initialNPoints =
data.Size();
718 std::vector<bool> isPointRejected(initialNPoints);
720 auto mapFunction = [&](
const unsigned int i) {
722 std::vector<double> gradFunc(npar);
723 std::vector<double> pointContribution(npar);
725 const auto x1 =
data.GetCoordComponent(i, 0);
726 const auto y =
data.Value(i);
727 auto invError =
data.Error(i);
729 invError = (invError != 0.0) ? 1.0 / invError : 1;
733 const double *
x =
nullptr;
734 std::vector<double> xc;
736 unsigned int ndim =
data.NDim();
737 double binVolume = 1;
740 for (
unsigned int j = 0; j < ndim; ++j) {
741 double x1_j = *
data.GetCoordComponent(i, j);
742 double x2_j =
data.GetBinUpEdgeComponent(i, j);
744 xc[j] = 0.5 * (x2_j + x1_j);
750 binVolume *= wrefVolume;
751 }
else if (ndim > 1) {
754 for (
unsigned int j = 1; j < ndim; ++j)
755 xc[j] = *
data.GetCoordComponent(i, j);
761 if (!useBinIntegral) {
763 func.ParameterGradient(
x,
p, &gradFunc[0]);
765 std::vector<double>
x2(
data.NDim());
766 data.GetBinUpEdgeCoordinates(i,
x2.data());
769 fval = igEval(
x,
x2.data());
776 std::cout <<
x[0] <<
" " <<
y <<
" " << 1. / invError <<
" params : ";
777 for (
unsigned int ipar = 0; ipar < npar; ++ipar)
778 std::cout <<
p[ipar] <<
"\t";
779 std::cout <<
"\tfval = " << fval << std::endl;
782 isPointRejected[i] =
true;
784 return pointContribution;
788 unsigned int ipar = 0;
789 for (; ipar < npar; ++ipar) {
793 gradFunc[ipar] *= binVolume;
797 double dfval = gradFunc[ipar];
803 pointContribution[ipar] = -2.0 * (
y - fval) * invError * invError * gradFunc[ipar];
808 isPointRejected[i] =
true;
811 return pointContribution;
815 auto redFunction = [&](
const std::vector<std::vector<double>> &pointContributions) {
816 std::vector<double>
result(npar);
818 for (
auto const &pointContribution : pointContributions) {
819 for (
unsigned int parameterIndex = 0; parameterIndex < npar; parameterIndex++)
820 result[parameterIndex] += pointContribution[parameterIndex];
826 std::vector<double>
g(npar);
831 Warning(
"FitUtil::EvaluateChi2Gradient",
"Multithread execution policy requires IMT, which is disabled. Changing "
832 "to ROOT::EExecutionPolicy::kSequential.");
838 std::vector<std::vector<double>> allGradients(initialNPoints);
839 for (
unsigned int i = 0; i < initialNPoints; ++i) {
840 allGradients[i] = mapFunction(i);
842 g = redFunction(allGradients);
856 Error(
"FitUtil::EvaluateChi2Gradient",
857 "Execution policy unknown. Available choices:\n 0: Serial (default)\n 1: MultiThread (requires IMT)\n");
866 nPoints = initialNPoints;
868 if (std::any_of(isPointRejected.begin(), isPointRejected.end(), [](
bool point) { return point; })) {
869 unsigned nRejected = std::accumulate(isPointRejected.begin(), isPointRejected.end(), 0);
870 assert(nRejected <= initialNPoints);
871 nPoints = initialNPoints - nRejected;
875 "Error - too many points rejected for overflow in gradient calculation");
879 std::copy(
g.begin(),
g.end(), grad);
898 const double *
x =
data.Coords(i);
899 double fval = func (
x,
p );
902 if (
g == 0)
return logPdf;
916 gc.ParameterGradient(
x,
p, fval,
g );
919 for (
unsigned int ipar = 0; ipar < func.
NPar(); ++ipar) {
924 std::cout <<
x[i] <<
"\t";
925 std::cout <<
"\tpar = [ " << func.
NPar() <<
" ] = ";
926 for (
unsigned int ipar = 0; ipar < func.
NPar(); ++ipar)
927 std::cout <<
p[ipar] <<
"\t";
928 std::cout <<
"\tfval = " << fval;
929 std::cout <<
"\tgrad = [ ";
930 for (
unsigned int ipar = 0; ipar < func.
NPar(); ++ipar)
931 std::cout <<
g[ipar] <<
"\t";
932 std::cout <<
" ] " << std::endl;
940 int iWeight,
bool extended,
unsigned int &nPoints,
945 unsigned int n =
data.Size();
949 bool normalizeFunc =
false;
953 (
const_cast<IModelFunctionTempl<double> &
>(func)).SetParameters(
p);
956 nPoints =
data.Size();
961 if (!normalizeFunc) {
962 if (
data.NDim() == 1) {
963 const double *
x =
data.GetCoordComponent(0,0);
967 std::vector<double>
x(
data.NDim());
968 for (
unsigned int j = 0; j <
data.NDim(); ++j)
969 x[j] = *
data.GetCoordComponent(0, j);
978 std::vector<double>
xmin(
data.NDim());
979 std::vector<double>
xmax(
data.NDim());
980 IntegralEvaluator<> igEval(func,
p,
true);
982 if (
data.Range().Size() > 0) {
984 for (
unsigned int ir = 0; ir <
data.Range().
Size(); ++ir) {
986 norm += igEval.Integral(
xmin.data(),
xmax.data());
992 if (func(
xmin.data(),
p) != 0 || func(
xmax.data(),
p) != 0) {
994 "A range has not been set and the function is not zero at +/- inf");
997 norm = igEval.Integral(&
xmin[0], &
xmax[0]);
1003 auto mapFunction = [&](
const unsigned i) {
1008 if (
data.NDim() > 1) {
1009 std::vector<double>
x(
data.NDim());
1010 for (
unsigned int j = 0; j <
data.NDim(); ++j)
1011 x[j] = *
data.GetCoordComponent(i, j);
1012#ifdef USE_PARAMCACHE
1013 fval = func(
x.data());
1015 fval = func(
x.data(),
p);
1020 const auto x =
data.GetCoordComponent(i, 0);
1021#ifdef USE_PARAMCACHE
1029 fval = fval * (1 / norm);
1034 double weight =
data.Weight(i);
1041 W2 = weight * weight;
1045 return LikelihoodAux<double>(logval, W, W2);
1056 auto redFunction = [](
const std::vector<LikelihoodAux<double>> & objs){
1057 auto l0 = LikelihoodAux<double>(0.0,0.0,0.0);
1058 for (
auto &
l : objs ) {
1068 Warning(
"FitUtil::EvaluateLogL",
"Multithread execution policy requires IMT, which is disabled. Changing "
1069 "to ROOT::EExecutionPolicy::kSequential.");
1078 for (
unsigned int i=0; i<
n; ++i) {
1079 auto resArray = mapFunction(i);
1080 logl+=resArray.logvalue;
1081 sumW+=resArray.weight;
1082 sumW2+=resArray.weight2;
1089 logl=resArray.logvalue;
1090 sumW=resArray.weight;
1091 sumW2=resArray.weight2;
1097 Error(
"FitUtil::EvaluateLogL",
"Execution policy unknown. Avalaible choices:\n ROOT::EExecutionPolicy::kSequential (default)\n ROOT::EExecutionPolicy::kMultiThread (requires IMT)\n");
1102 double extendedTerm = 0;
1106 if (!normalizeFunc) {
1107 IntegralEvaluator<> igEval( func,
p,
true);
1108 std::vector<double>
xmin(
data.NDim());
1109 std::vector<double>
xmax(
data.NDim());
1112 if (
data.Range().Size() > 0 ) {
1114 for (
unsigned int ir = 0; ir <
data.Range().
Size(); ++ir) {
1116 nuTot += igEval.Integral(
xmin.data(),
xmax.data());
1122 if (func(
xmin.data(),
p) != 0 || func(
xmax.data(),
p) != 0) {
1123 MATH_ERROR_MSG(
"FitUtil::EvaluateLogLikelihood",
"A range has not been set and the function is not zero at +/- inf");
1126 nuTot = igEval.Integral(&
xmin[0],&
xmax[0]);
1132 extendedTerm = - nuTot;
1136 extendedTerm = - (sumW2 / sumW) * nuTot;
1145 logl += extendedTerm;
1150 std::cout <<
"Evaluated log L for parameters (";
1151 for (
unsigned int ip = 0; ip < func.NPar(); ++ip)
1152 std::cout <<
" " <<
p[ip];
1153 std::cout <<
") fval = " << -logl << std::endl;
1165 assert(fg !=
nullptr);
1169 unsigned int npar = func.
NPar();
1170 unsigned initialNPoints =
data.Size();
1175 std::cout <<
"\n===> Evaluate Gradient for parameters ";
1176 for (
unsigned int ip = 0; ip < npar; ++ip)
1177 std::cout <<
" " <<
p[ip];
1181 const double kdmax1 =
std::sqrt(std::numeric_limits<double>::max());
1182 const double kdmax2 = std::numeric_limits<double>::max() / (4 * initialNPoints);
1184 auto mapFunction = [&](
const unsigned int i) {
1185 std::vector<double> gradFunc(npar);
1186 std::vector<double> pointContribution(npar);
1189 const double *
x =
nullptr;
1190 std::vector<double> xc;
1191 if (
data.NDim() > 1) {
1192 xc.resize(
data.NDim() );
1193 for (
unsigned int j = 0; j <
data.NDim(); ++j)
1194 xc[j] = *
data.GetCoordComponent(i, j);
1197 x =
data.GetCoordComponent(i, 0);
1200 double fval = func(
x,
p);
1201 func.ParameterGradient(
x,
p, &gradFunc[0]);
1206 if (i < 5 || (i >
data.Size()-5) ) {
1207 if (
data.NDim() > 1) std::cout << i <<
" x " <<
x[0] <<
" y " <<
x[1] <<
" func " << fval
1208 <<
" gradient " << gradFunc[0] <<
" " << gradFunc[1] <<
" " << gradFunc[3] << std::endl;
1209 else std::cout << i <<
" x " <<
x[0] <<
" gradient " << gradFunc[0] <<
" " << gradFunc[1] <<
" " << gradFunc[3] << std::endl;
1214 for (
unsigned int kpar = 0; kpar < npar; ++kpar) {
1216 pointContribution[kpar] = -1. / fval * gradFunc[kpar];
1217 else if (gradFunc[kpar] != 0) {
1218 double gg = kdmax1 * gradFunc[kpar];
1220 gg = std::min(gg, kdmax2);
1222 gg = std::max(gg, -kdmax2);
1223 pointContribution[kpar] = -gg;
1228 return pointContribution;
1232 auto redFunction = [&](
const std::vector<std::vector<double>> &pointContributions) {
1233 std::vector<double>
result(npar);
1235 for (
auto const &pointContribution : pointContributions) {
1236 for (
unsigned int parameterIndex = 0; parameterIndex < npar; parameterIndex++)
1237 result[parameterIndex] += pointContribution[parameterIndex];
1243 std::vector<double>
g(npar);
1248 Warning(
"FitUtil::EvaluateLogLGradient",
"Multithread execution policy requires IMT, which is disabled. Changing "
1249 "to ROOT::EExecutionPolicy::kSequential.");
1255 std::vector<std::vector<double>> allGradients(initialNPoints);
1256 for (
unsigned int i = 0; i < initialNPoints; ++i) {
1257 allGradients[i] = mapFunction(i);
1259 g = redFunction(allGradients);
1269 Error(
"FitUtil::EvaluateLogLGradient",
"Execution policy unknown. Avalaible choices:\n "
1270 "ROOT::EExecutionPolicy::kSequential (default)\n "
1271 "ROOT::EExecutionPolicy::kMultiThread (requires IMT)\n");
1280 std::copy(
g.begin(),
g.end(), grad);
1281 nPoints =
data.Size();
1284 std::cout <<
"FitUtil.cxx : Final gradient ";
1285 for (
unsigned int param = 0; param < npar; param++) {
1286 std::cout <<
" " << grad[param];
1299 const double *
x1 =
data.GetPoint(i,
y);
1302 bool useBinIntegral = fitOpt.
fIntegral &&
data.HasBinEdges();
1307 const double *
x2 = 0;
1309 double binVolume = 1;
1310 std::vector<double> xc;
1312 unsigned int ndim =
data.NDim();
1314 for (
unsigned int j = 0; j < ndim; ++j) {
1315 double x2j =
data.GetBinUpEdgeComponent(i, j);
1317 xc[j] = 0.5*(x2j+
x1[j]);
1320 binVolume /=
data.RefVolume();
1323 const double *
x = (useBinVolume) ? &xc.front() :
x1;
1325 if (!useBinIntegral ) {
1326 fval = func (
x,
p );
1330 std::vector<double> vx2(
data.NDim());
1331 data.GetBinUpEdgeCoordinates(i, vx2.data());
1332 fval = igEval(
x1, vx2.data() ) ;
1334 if (useBinVolume) fval *= binVolume;
1337 fval = std::max(fval, 0.0);
1338 double nlogPdf = fval;
1344 if (
g ==
nullptr)
return nlogPdf;
1346 unsigned int npar = func.
NPar();
1351 if (useFullHessian && (!gfunc || useBinIntegral || (gfunc && !gfunc->
HasParameterHessian())))
1352 return std::numeric_limits<double>::quiet_NaN();
1357 if (!useBinIntegral ) {
1359 if (useFullHessian &&
h) {
1361 return std::numeric_limits<double>::quiet_NaN();
1363 if (!goodHessFunc) {
1364 return std::numeric_limits<double>::quiet_NaN();
1376 if (!useBinIntegral )
1377 gc.ParameterGradient(
x,
p, fval,
g);
1384 double coeffGrad = (fval > 0) ? (1. -
y/fval) : ( (
y > 0) ?
std::sqrt( std::numeric_limits<double>::max() ) : 1. );
1385 double coeffHess = (fval > 0) ?
y/(fval*fval) : ( (
y > 0) ?
std::sqrt( std::numeric_limits<double>::max() ) : 0. );
1387 coeffGrad *= binVolume;
1388 coeffHess *= binVolume*binVolume;
1390 for (
unsigned int k = 0; k < npar; ++k) {
1393 for (
unsigned int l = k;
l < npar; ++
l) {
1394 unsigned int idx = k +
l * (
l + 1) / 2;
1395 if (useFullHessian) {
1396 h[idx] *= coeffGrad;
1402 h[idx] += coeffHess *
g[k]*
g[
l];
1413 std::cout <<
"x = " <<
x[0] <<
" y " <<
y <<
" fval " << fval <<
" logPdf = " << nlogPdf <<
" gradient : ";
1414 for (
unsigned int ipar = 0; ipar < npar; ++ipar)
1415 std::cout <<
g[ipar] <<
"\t";
1417 std::cout <<
"\thessian : ";
1418 for (
unsigned int ipar = 0; ipar < npar; ++ipar) {
1420 for (
unsigned int jpar = 0; jpar <= ipar; ++jpar) {
1421 std::cout <<
h[ipar + jpar * (jpar + 1) / 2] <<
"\t";
1426 std::cout << std::endl;
1454 unsigned int n =
data.Size();
1456#ifdef USE_PARAMCACHE
1460 nPoints =
data.Size();
1465 bool useBinIntegral = fitOpt.
fIntegral &&
data.HasBinEdges();
1467 bool useW2 = (iWeight == 2);
1470 double wrefVolume = 1.0;
1477 std::cout <<
"Evaluate PoissonLogL for params = [ ";
1478 for (
unsigned int j = 0; j < func.NPar(); ++j) std::cout <<
p[j] <<
" , ";
1479 std::cout <<
"] - data size = " <<
n <<
" useBinIntegral " << useBinIntegral <<
" useBinVolume "
1480 << useBinVolume <<
" useW2 " << useW2 <<
" wrefVolume = " << wrefVolume << std::endl;
1489#ifdef USE_PARAMCACHE
1490 IntegralEvaluator<> igEval(func, 0, useBinIntegral, igType);
1492 IntegralEvaluator<> igEval(func,
p, useBinIntegral, igType);
1495 auto mapFunction = [&](
const unsigned i) {
1496 auto x1 =
data.GetCoordComponent(i, 0);
1497 auto y = *
data.ValuePtr(i);
1499 const double *
x =
nullptr;
1500 std::vector<double> xc;
1502 double binVolume = 1.0;
1505 unsigned int ndim =
data.NDim();
1506 xc.resize(
data.NDim());
1507 for (
unsigned int j = 0; j < ndim; ++j) {
1508 double xx = *
data.GetCoordComponent(i, j);
1509 double x2 =
data.GetBinUpEdgeComponent(i, j);
1511 xc[j] = 0.5 * (
x2 + xx);
1515 binVolume *= wrefVolume;
1516 }
else if (
data.NDim() > 1) {
1517 xc.resize(
data.NDim());
1519 for (
unsigned int j = 1; j <
data.NDim(); ++j) {
1520 xc[j] = *
data.GetCoordComponent(i, j);
1527 if (!useBinIntegral) {
1528#ifdef USE_PARAMCACHE
1536 std::vector<double>
x2(
data.NDim());
1537 data.GetBinUpEdgeCoordinates(i,
x2.data());
1538 fval = igEval(
x,
x2.data());
1540 if (useBinVolume) fval *= binVolume;
1546 if (i % NSAMPLE == 0) {
1547 std::cout <<
"evt " << i <<
" x = [ ";
1548 for (
unsigned int j = 0; j < func.NDim(); ++j) std::cout <<
x[j] <<
" , ";
1551 std::cout <<
"x2 = [ ";
1552 for (
unsigned int j = 0; j < func.NDim(); ++j) std::cout <<
data.GetBinUpEdgeComponent(i, j) <<
" , ";
1555 std::cout <<
" y = " <<
y <<
" fval = " << fval << std::endl;
1562 fval = std::max(fval, 0.0);
1564 double nloglike = 0;
1573 double weight = 1.0;
1575 double error =
data.Error(i);
1576 weight = (error * error) /
y;
1581 weight =
data.SumOfError2()/
data.SumOfContent();
1584 nloglike += weight * ( fval -
y);
1592 if (extended) nloglike = fval -
y;
1601 std::cout <<
" nll = " << nloglike << std::endl;
1608 auto redFunction = [](
const std::vector<double> &objs) {
1609 return std::accumulate(objs.begin(), objs.end(),
double{});
1616 Warning(
"FitUtil::EvaluatePoissonLogL",
"Multithread execution policy requires IMT, which is disabled. Changing "
1617 "to ROOT::EExecutionPolicy::kSequential.");
1624 for (
unsigned int i = 0; i <
n; ++i) {
1625 res += mapFunction(i);
1637 Error(
"FitUtil::EvaluatePoissonLogL",
1638 "Execution policy unknown. Avalaible choices:\n ROOT::EExecutionPolicy::kSequential (default)\n ROOT::EExecutionPolicy::kMultiThread (requires IMT)\n");
1642 std::cout <<
"Loglikelihood = " << res << std::endl;
1654 assert(fg !=
nullptr);
1658#ifdef USE_PARAMCACHE
1662 const DataOptions &fitOpt =
data.Opt();
1663 bool useBinIntegral = fitOpt.fIntegral &&
data.HasBinEdges();
1664 bool useBinVolume = (fitOpt.fBinVolume &&
data.HasBinEdges());
1666 double wrefVolume = 1.0;
1667 if (useBinVolume && fitOpt.fNormBinVolume)
1668 wrefVolume /=
data.RefVolume();
1676 IntegralEvaluator<> igEval(func,
p, useBinIntegral, igType);
1678 unsigned int npar = func.NPar();
1679 unsigned initialNPoints =
data.Size();
1681 auto mapFunction = [&](
const unsigned int i) {
1683 std::vector<double> gradFunc(npar);
1684 std::vector<double> pointContribution(npar);
1686 const auto x1 =
data.GetCoordComponent(i, 0);
1687 const auto y =
data.Value(i);
1688 auto invError =
data.Error(i);
1690 invError = (invError != 0.0) ? 1.0 / invError : 1;
1694 const double *
x =
nullptr;
1695 std::vector<double> xc;
1697 unsigned ndim =
data.NDim();
1698 double binVolume = 1.0;
1703 for (
unsigned int j = 0; j < ndim; ++j) {
1704 double x1_j = *
data.GetCoordComponent(i, j);
1705 double x2_j =
data.GetBinUpEdgeComponent(i, j);
1706 binVolume *=
std::abs(x2_j - x1_j);
1707 xc[j] = 0.5 * (x2_j + x1_j);
1713 binVolume *= wrefVolume;
1714 }
else if (ndim > 1) {
1717 for (
unsigned int j = 1; j < ndim; ++j)
1718 xc[j] = *
data.GetCoordComponent(i, j);
1724 if (!useBinIntegral) {
1726 func.ParameterGradient(
x,
p, &gradFunc[0]);
1730 std::vector<double>
x2(
data.NDim());
1731 data.GetBinUpEdgeCoordinates(i,
x2.data());
1732 fval = igEval(
x,
x2.data());
1741 if (i < 5 || (i >
data.Size()-5) ) {
1742 if (
data.NDim() > 1) std::cout << i <<
" x " <<
x[0] <<
" y " <<
x[1] <<
" func " << fval
1743 <<
" gradient " << gradFunc[0] <<
" " << gradFunc[1] <<
" " << gradFunc[3] << std::endl;
1744 else std::cout << i <<
" x " <<
x[0] <<
" gradient " << gradFunc[0] <<
" " << gradFunc[1] <<
" " << gradFunc[3] << std::endl;
1750 for (
unsigned int ipar = 0; ipar < npar; ++ipar) {
1754 gradFunc[ipar] *= binVolume;
1758 pointContribution[ipar] = gradFunc[ipar] * (1. -
y / fval);
1759 else if (gradFunc[ipar] != 0) {
1760 const double kdmax1 =
std::sqrt(std::numeric_limits<double>::max());
1761 const double kdmax2 = std::numeric_limits<double>::max() / (4 * initialNPoints);
1762 double gg = kdmax1 * gradFunc[ipar];
1764 gg = std::min(gg, kdmax2);
1766 gg = std::max(gg, -kdmax2);
1767 pointContribution[ipar] = -gg;
1772 return pointContribution;
1776 auto redFunction = [&](
const std::vector<std::vector<double>> &pointContributions) {
1777 std::vector<double>
result(npar);
1779 for (
auto const &pointContribution : pointContributions) {
1780 for (
unsigned int parameterIndex = 0; parameterIndex < npar; parameterIndex++)
1781 result[parameterIndex] += pointContribution[parameterIndex];
1787 std::vector<double>
g(npar);
1792 Warning(
"FitUtil::EvaluatePoissonLogLGradient",
1793 "Multithread execution policy requires IMT, which is disabled. Changing "
1794 "to ROOT::EExecutionPolicy::kSequential.");
1800 std::vector<std::vector<double>> allGradients(initialNPoints);
1801 for (
unsigned int i = 0; i < initialNPoints; ++i) {
1802 allGradients[i] = mapFunction(i);
1804 g = redFunction(allGradients);
1819 Error(
"FitUtil::EvaluatePoissonLogLGradient",
1820 "Execution policy unknown. Avalaible choices:\n 0: Serial (default)\n 1: MultiThread (requires IMT)\n");
1829 std::copy(
g.begin(),
g.end(), grad);
1832 std::cout <<
"***** Final gradient : ";
1833 for (
unsigned int ii = 0; ii< npar; ++ii) std::cout << grad[ii] <<
" ";
1842 if (nEvents/ncpu < 1000)
return ncpu;
1843 return nEvents/1000;
#define MATH_ERROR_MSG(loc, str)
void Error(const char *location, const char *msgfmt,...)
Use this function in case an error occurred.
void Warning(const char *location, const char *msgfmt,...)
Use this function in warning situations.
winID h TVirtualViewer3D TVirtualGLPainter p
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void data
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t Float_t Float_t Int_t Int_t UInt_t UInt_t Rectangle_t result
Option_t Option_t TPoint TPoint const char x2
Option_t Option_t TPoint TPoint const char x1
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t Float_t g
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void gc
R__EXTERN TVirtualMutex * gROOTMutex
#define R__LOCKGUARD(mutex)
Class describing the binned data sets : vectors of x coordinates, y values and optionally error on y ...
SimpleGradientCalculator(int gdim, const IModelFunction &func, double eps=2.E-8, int istrat=1)
const IModelFunction & fFunc
std::vector< double > fVec
double ParameterDerivative(const double *x, const double *p, int ipar) const
void Gradient(const double *x, const double *p, double f0, double *g)
unsigned int NDim() const
void ParameterGradient(const double *x, const double *p, double f0, double *g)
unsigned int NPar() const
double DoParameterDerivative(const double *x, const double *p, double f0, int k) const
Class describing the un-binned data sets (just x coordinates values) of any dimensions.
virtual unsigned int NDim() const =0
Retrieve the dimension of the function.
virtual void SetParameters(const double *p)=0
Set the parameter values.
virtual unsigned int NPar() const =0
Return the number of Parameters.
Interface (abstract class) for parametric gradient multi-dimensional functions providing in addition ...
virtual void ParameterGradient(const T *x, const double *p, T *grad) const
Evaluate the all the derivatives (gradient vector) of the function with respect to the parameters at ...
virtual bool HasParameterHessian() const
virtual bool ParameterHessian(const T *, const double *, T *) const
Evaluate the all the Hessian (second derivatives matrix) of the function with respect to the paramete...
OneDimMultiFunctionAdapter class to wrap a multidimensional function in one dimensional one.
void SetCoord(int icoord)
User class for calculating the derivatives of a function.
double Derivative1(double x)
Returns the first derivative of the function at point x, computed by Richardson's extrapolation metho...
A pseudo container class which is a generator of indices.
This class provides a simple interface to execute the same task multiple times in parallel threads,...
auto MapReduce(F func, unsigned nTimes, R redfunc) -> InvokeResult_t< F >
Execute a function nTimes in parallel (Map) and accumulate the results into a single value (Reduce).
Type
enumeration specifying the integration types.
@ kGAUSS
simple Gauss integration method with fixed rule
@ kDEFAULT
default type specified in the static options
RVec< PromoteType< T > > abs(const RVec< T > &v)
TFitResultPtr Fit(FitObject *h1, TF1 *f1, Foption_t &option, const ROOT::Math::MinimizerOptions &moption, const char *goption, ROOT::Fit::DataRange &range)
ROOT::Math::IParamMultiGradFunction IGradModelFunction
ROOT::Math::IParamMultiFunction IModelFunction
void CalculateGradientIntegral(const GFunc &gfunc, const double *x1, const double *x2, const double *p, double *g)
double EvaluatePoissonBinPdf(const IModelFunction &func, const BinData &data, const double *x, unsigned int ipoint, double *g=nullptr, double *h=nullptr, bool hasGrad=false, bool fullHessian=false)
evaluate the pdf contribution to the Poisson LogL given a model function and the BinPoint data.
double CorrectValue(double rval)
void EvaluatePoissonLogLGradient(const IModelFunction &func, const BinData &data, const double *p, double *grad, unsigned int &nPoints, ::ROOT::EExecutionPolicy executionPolicy=::ROOT::EExecutionPolicy::kSequential, unsigned nChunks=0)
evaluate the Poisson LogL given a model function and the data at the point p.
double EvaluateChi2Residual(const IModelFunction &func, const BinData &data, const double *p, unsigned int ipoint, double *g=nullptr, double *h=nullptr, bool hasGrad=false, bool fullHessian=false)
evaluate the residual contribution to the Chi2 given a model function and the BinPoint data and if th...
double EvaluatePoissonLogL(const IModelFunction &func, const BinData &data, const double *p, int iWeight, bool extended, unsigned int &nPoints, ::ROOT::EExecutionPolicy executionPolicy, unsigned nChunks=0)
evaluate the Poisson LogL given a model function and the data at the point p.
double EvaluateChi2Effective(const IModelFunction &func, const BinData &data, const double *x, unsigned int &nPoints)
evaluate the effective Chi2 given a model function and the data at the point x.
void EvaluateLogLGradient(const IModelFunction &func, const UnBinData &data, const double *p, double *grad, unsigned int &nPoints, ::ROOT::EExecutionPolicy executionPolicy=::ROOT::EExecutionPolicy::kSequential, unsigned nChunks=0)
evaluate the LogL gradient given a model function and the data at the point p.
double EvaluatePdf(const IModelFunction &func, const UnBinData &data, const double *p, unsigned int ipoint, double *g=nullptr, double *h=nullptr, bool hasGrad=false, bool fullHessian=false)
evaluate the pdf contribution to the LogL given a model function and the BinPoint data.
bool CheckInfNaNValue(double &rval)
unsigned setAutomaticChunking(unsigned nEvents)
double EvaluateLogL(const IModelFunction &func, const UnBinData &data, const double *p, int iWeight, bool extended, unsigned int &nPoints, ::ROOT::EExecutionPolicy executionPolicy, unsigned nChunks=0)
evaluate the LogL given a model function and the data at the point x.
double EvaluateChi2(const IModelFunction &func, const BinData &data, const double *p, unsigned int &nPoints, ::ROOT::EExecutionPolicy executionPolicy, unsigned nChunks=0)
Chi2 Functions.
void EvaluateChi2Gradient(const IModelFunction &func, const BinData &data, const double *p, double *grad, unsigned int &nPoints, ::ROOT::EExecutionPolicy executionPolicy=::ROOT::EExecutionPolicy::kSequential, unsigned nChunks=0)
evaluate the Chi2 gradient given a model function and the data at the point p.
void(off) SmallVectorTemplateBase< T
T EvalLog(T x)
safe evaluation of log(x) with a protections against negative or zero argument to the log smooth line...
VecExpr< UnaryOp< Sqrt< T >, VecExpr< A, T, D >, T >, T, D > sqrt(const VecExpr< A, T, D > &rhs)
This file contains a specialised ROOT message handler to test for diagnostic in unit tests.
UInt_t GetThreadPoolSize()
Returns the size of ROOT's thread pool.
DataOptions : simple structure holding the options on how the data are filled.
bool fErrors1
use all errors equal to 1, i.e. fit without errors (default is false)
bool fNormBinVolume
normalize data by a normalized the bin volume (bin volume divided by a reference value)
bool fUseEmpty
use empty bins (default is false) with a fixed error of 1
bool fIntegral
use integral of bin content instead of bin center (default is false)
bool fExpErrors
use expected errors from the function and not from the data
bool fBinVolume
normalize data by the bin volume (it is used in the Poisson likelihood fits)
double operator()(const double *x, const double *p) const
void SetDerivComponent(unsigned int ipar)
unsigned int NDim() const
ParamDerivFunc(const GradFunc &f)