Skip to content

Commit 111257e

Browse files
guitargeekguitargeek
committed
[math] Support multithreaded fits of TGraphErrors with x errors
Make FitUtil::EvaluateChi2Effective execution-policy aware via TThreadExecutor (same pattern as EvaluateChi2) and thread the policy through Chi2FCN. Point coordinates/errors are read via the thread-safe component accessors, since GetPoint/GetPointError use shared temp buffers. Also honour the SERIAL/MULTITHREAD options for graph fits. Adds exec-policy tests to testFitter. Closes #10021. 🤖 Done with the help of AI.
1 parent 5bcfc92 commit 111257e

7 files changed

Lines changed: 171 additions & 40 deletions

File tree

‎hist/hist/src/HFitImpl.cxx‎

Lines changed: 12 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -695,6 +695,18 @@ void ROOT::Fit::FitOptionsMake(EFitObjectType type, const char *option, Foption_
695695
TString opt = option;
696696
opt.ToUpper();
697697

698+
// Parse the execution policy options first and strip them from the option
699+
// string, so that the remaining letters (e.g. the "T", "E", "R" in
700+
// "MULTITHREAD") are not mistaken for single-letter options below.
701+
if (opt.Contains("SERIAL")) {
702+
fitOption.ExecPolicy = ROOT::EExecutionPolicy::kSequential;
703+
opt.ReplaceAll("SERIAL", "");
704+
}
705+
if (opt.Contains("MULTITHREAD")) {
706+
fitOption.ExecPolicy = ROOT::EExecutionPolicy::kMultiThread;
707+
opt.ReplaceAll("MULTITHREAD", "");
708+
}
709+
698710
// parse firt the specific options
699711
if (type == EFitObjectType::kHistogram) {
700712

@@ -715,16 +727,6 @@ void ROOT::Fit::FitOptionsMake(EFitObjectType type, const char *option, Foption_
715727
// opt.ReplaceAll("MULTIPROC","");
716728
// }
717729

718-
if (opt.Contains("SERIAL")) {
719-
fitOption.ExecPolicy = ROOT::EExecutionPolicy::kSequential;
720-
opt.ReplaceAll("SERIAL","");
721-
}
722-
723-
if (opt.Contains("MULTITHREAD")) {
724-
fitOption.ExecPolicy = ROOT::EExecutionPolicy::kMultiThread;
725-
opt.ReplaceAll("MULTITHREAD","");
726-
}
727-
728730
if (opt.Contains("I")) fitOption.Integral= 1; // integral of function in the bin (no sense for graph)
729731
if (opt.Contains("W")) fitOption.W1 = 1; // all non-empty bins or points have weight =1 (for chi2 fit)
730732
if (opt.Contains("WW")) fitOption.W1 = 2; //all bins have weight=1, even empty bins

‎hist/hist/src/TGraph.cxx‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1181,6 +1181,8 @@ TObject *TGraph::FindObject(const TObject *obj) const
11811181
/// "G" | Uses the gradient implemented in `TF1::GradientPar` for the minimization. This allows to use Automatic Differentiation when it is supported by the provided TF1 function.
11821182
/// "EX0" | When fitting a TGraphErrors or TGraphAsymErrors do not consider errors in the X coordinates
11831183
/// "ROB" | In case of linear fitting, compute the LTS regression coefficients (robust (resistant) regression), using the default fraction of good points "ROB=0.x" - compute the LTS regression coefficients, using 0.x as a fraction of good points
1184+
/// "SERIAL" | Runs in serial mode. By default, if ROOT is built with MT support and MT is enabled, the fit is performed in multi-thread.
1185+
/// "MULTITHREAD" | Forces usage of multi-thread execution whenever possible.
11841186
///
11851187
///
11861188
/// This function is used for fitting also the derived TGraph classes such as TGraphErrors or TGraphAsymmErrors.

‎math/mathcore/inc/Fit/Chi2FCN.h‎

Lines changed: 6 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -141,12 +141,15 @@ class Chi2FCN : public BasicFCN<DerivFunType, ModelFunType, BinData> {
141141
/**
142142
Evaluation of the function (required by interface)
143143
*/
144-
double DoEval (const double * x) const override {
144+
double DoEval(const double *x) const override
145+
{
145146
this->UpdateNCalls();
146147
if (BaseFCN::Data().HaveCoordErrors() || BaseFCN::Data().HaveAsymErrors())
147-
return FitUtil::Evaluate<T>::EvalChi2Effective(BaseFCN::ModelFunction(), BaseFCN::Data(), x, fNEffPoints);
148+
return FitUtil::Evaluate<T>::EvalChi2Effective(BaseFCN::ModelFunction(), BaseFCN::Data(), x, fNEffPoints,
149+
fExecutionPolicy);
148150
else
149-
return FitUtil::Evaluate<T>::EvalChi2(BaseFCN::ModelFunction(), BaseFCN::Data(), x, fNEffPoints, fExecutionPolicy);
151+
return FitUtil::Evaluate<T>::EvalChi2(BaseFCN::ModelFunction(), BaseFCN::Data(), x, fNEffPoints,
152+
fExecutionPolicy);
150153
}
151154

152155
// for derivatives

‎math/mathcore/inc/Fit/FitUtil.h‎

Lines changed: 8 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -263,7 +263,9 @@ namespace FitUtil {
263263
The effective chi2 uses the errors on the coordinates : W = 1/(sigma_y**2 + ( sigma_x_i * df/dx_i )**2 )
264264
return also nPoints as the effective number of used points in the Chi2 evaluation
265265
*/
266-
double EvaluateChi2Effective(const IModelFunction &func, const BinData &data, const double *x, unsigned int &nPoints);
266+
double EvaluateChi2Effective(const IModelFunction &func, const BinData &data, const double *x, unsigned int &nPoints,
267+
::ROOT::EExecutionPolicy executionPolicy = ::ROOT::EExecutionPolicy::kSequential,
268+
unsigned nChunks = 0);
267269

268270
/**
269271
evaluate the Chi2 gradient given a model function and the data at the point p.
@@ -363,7 +365,8 @@ namespace FitUtil {
363365
unsigned nChunks = 0);
364366

365367
static double
366-
EvalChi2Effective(const IModelFunctionTempl<Double_v> &, const BinData &, const double *, unsigned int &)
368+
EvalChi2Effective(const IModelFunctionTempl<Double_v> &, const BinData &, const double *, unsigned int &,
369+
::ROOT::EExecutionPolicy = ::ROOT::EExecutionPolicy::kSequential, unsigned = 0)
367370
{
368371
Error("FitUtil::Evaluate<T>::EvalChi2Effective", "The vectorized evaluation of the Chi2 with coordinate errors is still not supported");
369372
return -1.;
@@ -441,9 +444,10 @@ namespace FitUtil {
441444
return FitUtil::EvaluatePoissonLogL(func, data, p, iWeight, extended, nPoints, executionPolicy, nChunks);
442445
}
443446

444-
static double EvalChi2Effective(const IModelFunctionTempl<double> &func, const BinData & data, const double * p, unsigned int &nPoints)
447+
static double EvalChi2Effective(const IModelFunctionTempl<double> &func, const BinData & data, const double * p, unsigned int &nPoints,
448+
::ROOT::EExecutionPolicy executionPolicy = ::ROOT::EExecutionPolicy::kSequential, unsigned nChunks = 0)
445449
{
446-
return FitUtil::EvaluateChi2Effective(func, data, p, nPoints);
450+
return FitUtil::EvaluateChi2Effective(func, data, p, nPoints, executionPolicy, nChunks);
447451
}
448452
static void EvalChi2Gradient(const IModelFunctionTempl<double> &func, const BinData &data, const double *p,
449453
double *g, unsigned int &nPoints,

‎math/mathcore/src/FitUtil.cxx‎

Lines changed: 52 additions & 22 deletions
Original file line numberDiff line numberDiff line change
@@ -425,7 +425,8 @@ double FitUtil::EvaluateChi2(const IModelFunction &func, const BinData &data, co
425425

426426
//___________________________________________________________________________________________________________________________
427427

428-
double FitUtil::EvaluateChi2Effective(const IModelFunction & func, const BinData & data, const double * p, unsigned int & nPoints) {
428+
double FitUtil::EvaluateChi2Effective(const IModelFunction & func, const BinData & data, const double * p, unsigned int & nPoints,
429+
::ROOT::EExecutionPolicy executionPolicy, unsigned nChunks) {
429430
// evaluate the chi2 given a function reference , the data and returns the value and also in nPoints
430431
// the actual number of used points
431432
// method using the error in the coordinates
@@ -440,39 +441,42 @@ double FitUtil::EvaluateChi2Effective(const IModelFunction & func, const BinData
440441

441442
assert(data.HaveCoordErrors() || data.HaveAsymErrors());
442443

443-
double chi2 = 0;
444-
//int nRejected = 0;
445-
446-
447444
//func.SetParameters(p);
448445

449446
unsigned int ndim = func.NDim();
450447

451-
// use Richardson derivator
452-
ROOT::Math::RichardsonDerivator derivator;
453-
454448
double maxResValue = std::numeric_limits<double>::max() /n;
455449

450+
auto mapFunction = [&](const unsigned i) {
456451

452+
// use a Richardson derivator local to the iteration since it is not thread safe
453+
ROOT::Math::RichardsonDerivator derivator;
457454

458-
for (unsigned int i = 0; i < n; ++ i) {
459-
455+
// copy the coordinates and the coordinate errors of the point into thread-local
456+
// buffers: the BinData::GetPoint and GetPointError accessors return pointers into
457+
// shared temporary storage and are not thread safe (in contrast to the per-component
458+
// accessors used here)
459+
// TODO: add threadsafe getters to FitData and BinData!
460+
std::vector<double> x(ndim);
461+
std::vector<double> ex(ndim);
462+
for (unsigned int icoord = 0; icoord < ndim; ++icoord) {
463+
x[icoord] = *data.GetCoordComponent(i, icoord);
464+
ex[icoord] = data.GetCoordErrorComponent(i, icoord);
465+
}
460466

461-
double y = 0;
462-
const double * x = data.GetPoint(i,y);
467+
double y = data.Value(i);
463468

464-
double fval = func( x, p );
469+
double fval = func( x.data(), p );
465470

466471
double delta_y_func = y - fval;
467472

468473

469474
double ey = 0;
470-
const double * ex = nullptr;
471475
if (!data.HaveAsymErrors() )
472-
ex = data.GetPointError(i, ey);
476+
ey = data.Error(i);
473477
else {
474-
double eylow, eyhigh = 0;
475-
ex = data.GetPointError(i, eylow, eyhigh);
478+
double eylow = 0, eyhigh = 0;
479+
data.GetAsymError(i, eylow, eyhigh);
476480
if ( delta_y_func < 0)
477481
ey = eyhigh; // function is higher than points
478482
else
@@ -485,7 +489,7 @@ double FitUtil::EvaluateChi2Effective(const IModelFunction & func, const BinData
485489
// if j is less ndim some elements are not zero
486490
if (j < ndim) {
487491
// need an adapter from a multi-dim function to a one-dimensional
488-
ROOT::Math::OneDimMultiFunctionAdapter<const IModelFunction &> f1D(func,x,0,p);
492+
ROOT::Math::OneDimMultiFunctionAdapter<const IModelFunction &> f1D(func,x.data(),0,p);
489493
// select optimal step size (use 10--2 by default as was done in TF1:
490494
double kEps = 0.01;
491495
double kPrecision = 1.E-8;
@@ -518,12 +522,38 @@ double FitUtil::EvaluateChi2Effective(const IModelFunction & func, const BinData
518522

519523
// avoid (infinity and nan ) in the chi2 sum
520524
// eventually add possibility of excluding some points (like singularity)
521-
if ( resval < maxResValue )
522-
chi2 += resval;
523-
else
524-
chi2 += maxResValue;
525+
return ( resval < maxResValue ) ? resval : maxResValue;
525526
//nRejected++;
527+
};
528+
529+
#ifdef R__USE_IMT
530+
auto redFunction = [](const std::vector<double> & objs) {
531+
return std::accumulate(objs.begin(), objs.end(), double{});
532+
};
533+
#else
534+
(void)nChunks;
526535

536+
// If IMT is disabled, force the execution policy to the serial case
537+
if (executionPolicy == ROOT::EExecutionPolicy::kMultiThread) {
538+
Warning("FitUtil::EvaluateChi2Effective", "Multithread execution policy requires IMT, which is disabled. Changing "
539+
"to ROOT::EExecutionPolicy::kSequential.");
540+
executionPolicy = ROOT::EExecutionPolicy::kSequential;
541+
}
542+
#endif
543+
544+
double chi2 = 0;
545+
if (executionPolicy == ROOT::EExecutionPolicy::kSequential) {
546+
for (unsigned int i = 0; i < n; ++i) {
547+
chi2 += mapFunction(i);
548+
}
549+
#ifdef R__USE_IMT
550+
} else if (executionPolicy == ROOT::EExecutionPolicy::kMultiThread) {
551+
ROOT::TThreadExecutor pool;
552+
auto chunks = nChunks != 0 ? nChunks : setAutomaticChunking(data.Size());
553+
chi2 = pool.MapReduce(mapFunction, ROOT::TSeq<unsigned>(0, n), redFunction, chunks);
554+
#endif
555+
} else {
556+
Error("FitUtil::EvaluateChi2Effective", "Execution policy unknown. Available choices:\n ROOT::EExecutionPolicy::kSequential (default)\n ROOT::EExecutionPolicy::kMultiThread (requires IMT)\n");
527557
}
528558

529559
// reset the number of fitting data points

‎math/mathcore/test/CMakeLists.txt‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -77,7 +77,7 @@ ROOT_ADD_GTEST(MulmodUnitOpt mulmod_opt.cxx)
7777
ROOT_ADD_GTEST(RanluxLCGUnit ranlux_lcg.cxx)
7878
ROOT_ADD_GTEST(RanluxppEngineTests RanluxppEngine.cxx LIBRARIES Core MathCore)
7979
ROOT_ADD_GTEST(testDelaunay2D testDelaunay2D.cxx LIBRARIES Core MathCore)
80-
ROOT_ADD_GTEST(testFitter testFitter.cxx LIBRARIES Core MathCore)
80+
ROOT_ADD_GTEST(testFitter testFitter.cxx LIBRARIES Core MathCore Hist)
8181
ROOT_ADD_GTEST(testKNNDensity testKNNDensity.cxx LIBRARIES Core MathCore Hist)
8282
ROOT_ADD_GTEST(testKahan testKahan.cxx LIBRARIES Core MathCore)
8383
ROOT_ADD_GTEST(testLFSR testLFSR.cxx LIBRARIES Core MathCore)

‎math/mathcore/test/testFitter.cxx‎

Lines changed: 90 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -2,9 +2,18 @@
22

33
#include <Fit/Fitter.h>
44
#include <Math/Functor.h>
5+
#include <Math/MinimizerOptions.h>
6+
7+
#include <TGraphErrors.h>
8+
#include <TGraphAsymmErrors.h>
9+
#include <TF1.h>
10+
#include <TFitResult.h>
511

612
#include <gtest/gtest.h>
713

14+
#include <cmath>
15+
#include <vector>
16+
817
namespace {
918

1019
double gauss(double *x, double *p)
@@ -16,6 +25,66 @@ double gauss(double *x, double *p)
1625
return A * exp(-(xx - mu) * (xx - mu) / (2 * sigma * sigma));
1726
}
1827

28+
#ifdef R__USE_IMT
29+
30+
// Build a deterministic dataset so the exec-policy tests are reproducible.
31+
TGraphErrors makeGraphErrors(int n)
32+
{
33+
std::vector<double> x(n), y(n), ex(n), ey(n);
34+
for (int i = 0; i < n; ++i) {
35+
double xi = 0.01 * i;
36+
x[i] = xi;
37+
y[i] = 3.0 + 2.0 * xi + 0.1 * std::sin(xi);
38+
ex[i] = 0.05;
39+
ey[i] = 0.2;
40+
}
41+
return TGraphErrors(n, x.data(), y.data(), ex.data(), ey.data());
42+
}
43+
44+
TGraphAsymmErrors makeGraphAsymmErrors(int n)
45+
{
46+
std::vector<double> x(n), y(n), exl(n), exh(n), eyl(n), eyh(n);
47+
for (int i = 0; i < n; ++i) {
48+
double xi = 0.01 * i;
49+
x[i] = xi;
50+
y[i] = 3.0 + 2.0 * xi + 0.1 * std::sin(xi);
51+
exl[i] = 0.04;
52+
exh[i] = 0.06;
53+
eyl[i] = 0.15;
54+
eyh[i] = 0.25;
55+
}
56+
return TGraphAsymmErrors(n, x.data(), y.data(), exl.data(), exh.data(), eyl.data(), eyh.data());
57+
}
58+
59+
// Fit a graph with errors on the x coordinate (which uses the "effective chi2",
60+
// FitUtil::EvaluateChi2Effective) both sequentially and multi-threaded, and check
61+
// the two execution policies converge to the same minimum.
62+
template <class Graph>
63+
void checkEffectiveChi2ExecPolicy(Graph &g)
64+
{
65+
ROOT::Math::MinimizerOptions::SetDefaultMinimizer("Minuit2");
66+
67+
TF1 fSeq("fSeq", "[0]+[1]*x", 0, 20);
68+
fSeq.SetParameters(1, 1);
69+
auto rSeq = g.Fit(&fSeq, "S Q N SERIAL");
70+
ASSERT_TRUE(rSeq.Get() != nullptr);
71+
ASSERT_TRUE(rSeq->IsValid());
72+
73+
TF1 fMT("fMT", "[0]+[1]*x", 0, 20);
74+
fMT.SetParameters(1, 1);
75+
auto rMT = g.Fit(&fMT, "S Q N MULTITHREAD");
76+
ASSERT_TRUE(rMT.Get() != nullptr);
77+
ASSERT_TRUE(rMT->IsValid());
78+
79+
// The tiny residual difference is only due to the order of the floating point
80+
// sum in the parallel reduction.
81+
EXPECT_NEAR(rSeq->Parameter(0), rMT->Parameter(0), 1e-6);
82+
EXPECT_NEAR(rSeq->Parameter(1), rMT->Parameter(1), 1e-6);
83+
EXPECT_NEAR(rSeq->MinFcnValue(), rMT->MinFcnValue(), 1e-4 * std::abs(rSeq->MinFcnValue()));
84+
}
85+
86+
#endif // R__USE_IMT
87+
1988
} // namespace
2089

2190
/// Check if we can fix and release parameters, and that this will be correctly
@@ -75,3 +144,24 @@ TEST(Fitter, FitAndReleaseParams)
75144
}
76145
}
77146
}
147+
148+
#ifdef R__USE_IMT
149+
150+
/// Fitting a TGraphErrors with errors on the x coordinate (effective chi2) must
151+
/// give the same result sequentially and multi-threaded.
152+
/// Covers https://github.com/root-project/root/issues/10021
153+
/// See also https://root-forum.cern.ch/t/multithreading-on-minuit/46225
154+
TEST(Fitter, EffectiveChi2ExecPolicy)
155+
{
156+
auto g = makeGraphErrors(2000);
157+
checkEffectiveChi2ExecPolicy(g);
158+
}
159+
160+
/// Same as above but with asymmetric errors (also goes through the effective chi2 path).
161+
TEST(Fitter, EffectiveChi2ExecPolicyAsymm)
162+
{
163+
auto g = makeGraphAsymmErrors(1500);
164+
checkEffectiveChi2ExecPolicy(g);
165+
}
166+
167+
#endif

0 commit comments

Comments
 (0)