From 297baa63644ca043f5f43b91017ca05d09e1e70c Mon Sep 17 00:00:00 2001 From: Phmonski Date: Fri, 17 Jul 2026 14:32:06 +0200 Subject: [PATCH 1/6] [RF][HS3] Add interpolation export for HistFactory modifiers --- .../HistFactory/PiecewiseInterpolation.h | 3 + .../src/PiecewiseInterpolation.cxx | 27 +- .../test/testPiecewiseInterpolation.cxx | 19 + roofit/hs3/src/JSONFactories_HistFactory.cxx | 357 +++++++++++++++--- roofit/hs3/test/testRooFitHS3.cxx | 322 +++++++++++++++- 5 files changed, 676 insertions(+), 52 deletions(-) diff --git a/roofit/histfactory/inc/RooStats/HistFactory/PiecewiseInterpolation.h b/roofit/histfactory/inc/RooStats/HistFactory/PiecewiseInterpolation.h index 4d6fc9e66d67d..03255d6257d07 100644 --- a/roofit/histfactory/inc/RooStats/HistFactory/PiecewiseInterpolation.h +++ b/roofit/histfactory/inc/RooStats/HistFactory/PiecewiseInterpolation.h @@ -33,6 +33,9 @@ class PiecewiseInterpolation : public RooAbsReal { PiecewiseInterpolation() ; PiecewiseInterpolation(const char *name, const char *title, const RooAbsReal &nominal, const RooArgList &lowSet, const RooArgList &highSet, const RooArgList ¶mSet); + PiecewiseInterpolation(const char *name, const char *title, const RooAbsReal &nominal, const RooArgList &lowSet, + const RooArgList &highSet, const RooArgList ¶mSet, + const std::vector &interpolationCodes); ~PiecewiseInterpolation() override ; PiecewiseInterpolation(const PiecewiseInterpolation& other, const char *name = nullptr); diff --git a/roofit/histfactory/src/PiecewiseInterpolation.cxx b/roofit/histfactory/src/PiecewiseInterpolation.cxx index e52fafb461804..fdf117d91f3f6 100644 --- a/roofit/histfactory/src/PiecewiseInterpolation.cxx +++ b/roofit/histfactory/src/PiecewiseInterpolation.cxx @@ -131,7 +131,32 @@ PiecewiseInterpolation::PiecewiseInterpolation(const char *name, const char *tit TRACE_CREATE; } - +//////////////////////////////////////////////////////////////////////////////// +/// Construct a new interpolation and set the interpolation code for each +/// parameter by position. +/// \param name Name of the object. +/// \param title Title (for e.g. plotting). +/// \param nominal Nominal value of the function. +/// \param lowSet Set of down variations. +/// \param highSet Set of up variations. +/// \param paramSet Parameters that control the interpolation. +/// \param interpolationCodes Interpolation code for each parameter. +PiecewiseInterpolation::PiecewiseInterpolation(const char *name, const char *title, const RooAbsReal &nominal, + const RooArgList &lowSet, const RooArgList &highSet, + const RooArgList ¶mSet, const std::vector &interpolationCodes) + : PiecewiseInterpolation(name, title, nominal, lowSet, highSet, paramSet) +{ + if (interpolationCodes.size() != _paramSet.size()) { + coutE(InputArguments) << "PiecewiseInterpolation::ctor(" << GetName() + << ") ERROR: interpolation code vector should have the same length as the parameter list" + << std::endl; + RooErrorHandler::softAbort(); + return; + } + for (std::size_t i = 0; i < interpolationCodes.size(); ++i) { + setInterpCodeForParam(i, interpolationCodes[i]); + } +} //////////////////////////////////////////////////////////////////////////////// /// Copy constructor diff --git a/roofit/histfactory/test/testPiecewiseInterpolation.cxx b/roofit/histfactory/test/testPiecewiseInterpolation.cxx index a51bad25fc645..fad4847965ea0 100644 --- a/roofit/histfactory/test/testPiecewiseInterpolation.cxx +++ b/roofit/histfactory/test/testPiecewiseInterpolation.cxx @@ -101,3 +101,22 @@ TEST(PiecewiseInterpolation, AdditiveOrMultiplicative) } } } + +TEST(PiecewiseInterpolation, ConstructWithPerParameterCodes) +{ + RooRealVar nominal{"nominal", "nominal", 1.0}; + RooRealVar parameter{"parameter", "parameter", 0.0, -2.0, 2.0}; + RooRealVar low1{"low1", "low1", 0.8}; + RooRealVar low2{"low2", "low2", 0.7}; + RooRealVar high1{"high1", "high1", 1.2}; + RooRealVar high2{"high2", "high2", 1.4}; + + // Repeating a parameter is supported by PiecewiseInterpolation, so codes + // need to be assigned by position rather than by looking up the parameter. + RooArgList parameters{parameter, parameter}; + RooArgList lows{low1, low2}; + RooArgList highs{high1, high2}; + PiecewiseInterpolation interpolation{"interpolation", "", nominal, lows, highs, parameters, {0, 2}}; + + EXPECT_EQ(interpolation.interpolationCodes(), (std::vector{0, 2})); +} diff --git a/roofit/hs3/src/JSONFactories_HistFactory.cxx b/roofit/hs3/src/JSONFactories_HistFactory.cxx index cc0d670ee99a8..4630f4e6b6989 100644 --- a/roofit/hs3/src/JSONFactories_HistFactory.cxx +++ b/roofit/hs3/src/JSONFactories_HistFactory.cxx @@ -33,7 +33,12 @@ #include #include +#include +#include +#include +#include #include +#include #include "static_execute.h" #include "JSONIOUtils.h" @@ -63,6 +68,209 @@ namespace Literals { constexpr auto staterror = "staterror"; } +struct Interpolation { + std::string type; + std::string in; + std::optional out; + + bool operator==(const Interpolation &other) const + { + return std::tie(type, in, out) == std::tie(other.type, other.in, other.out); + } + + bool operator!=(const Interpolation &other) const { return !(*this == other); } + + bool operator<(const Interpolation &other) const + { + return std::tie(type, in, out) < std::tie(other.type, other.in, other.out); + } +}; + +const Interpolation additivePiecewiseLinear{"add", "poly1", std::nullopt}; +const Interpolation multiplicativePiecewiseExponential{"mult", "exp", std::nullopt}; +const Interpolation additiveQuadraticLinear{"add", "poly2", "poly1"}; +const Interpolation additivePolynomialLinear{"add", "poly6", "poly1"}; +const Interpolation multiplicativePolynomialExponential{"mult", "poly6", "exp"}; +const Interpolation multiplicativePolynomialLinear{"mult", "poly6", "poly1"}; + +std::string interpolationString(const Interpolation &interpolation) +{ + std::stringstream ss; + ss << R"({"type":")" << interpolation.type << R"(","in":")" << interpolation.in << R"(","out":)"; + if (interpolation.out) { + ss << '"' << *interpolation.out << '"'; + } else { + ss << "null"; + } + ss << '}'; + return ss.str(); +} + +bool isInterpolationFunction(std::string_view function) +{ + return function == "poly1" || function == "poly2" || function == "poly6" || function == "exp"; +} + +Interpolation readInterpolation(const JSONNode &node, const std::string &context) +{ + if (!node.is_map()) { + RooJSONFactoryWSTool::error(context + " must be a struct with components 'type', 'in', and 'out'"); + } + for (const char *component : {"type", "in", "out"}) { + if (!node.has_child(component)) { + RooJSONFactoryWSTool::error(context + " does not define the required '" + component + "' component"); + } + } + + const auto &typeNode = node["type"]; + const auto &inNode = node["in"]; + const auto &outNode = node["out"]; + if (typeNode.is_container() || typeNode.is_null() || inNode.is_container() || inNode.is_null()) { + RooJSONFactoryWSTool::error(context + " components 'type' and 'in' must be strings"); + } + + Interpolation interpolation{typeNode.val(), inNode.val(), std::nullopt}; + if (interpolation.type != "add" && interpolation.type != "mult") { + RooJSONFactoryWSTool::error(context + " has unknown interpolation type '" + interpolation.type + "'"); + } + if (!isInterpolationFunction(interpolation.in)) { + RooJSONFactoryWSTool::error(context + " has unknown interpolation function '" + interpolation.in + "'"); + } + + if (!outNode.is_null()) { + if (outNode.is_container()) { + RooJSONFactoryWSTool::error(context + " component 'out' must be a string or null"); + } + interpolation.out = outNode.val(); + if (!isInterpolationFunction(*interpolation.out)) { + RooJSONFactoryWSTool::error(context + " has unknown extrapolation function '" + *interpolation.out + "'"); + } + } + + return interpolation; +} + +void writeInterpolation(JSONNode &node, const Interpolation &interpolation) +{ + node.set_map(); + node["type"] << interpolation.type; + node["in"] << interpolation.in; + if (interpolation.out) { + node["out"] << *interpolation.out; + } else { + node["out"].set_null(); + } +} + +int readLegacyInterpolationCode(const JSONNode &node, const std::string &context) +{ + if (node.is_container() || node.is_null() || !node.has_val()) { + RooJSONFactoryWSTool::error(context + " must be a structured interpolation or a legacy integer code"); + } + + const std::string value = node.val(); + int code = 0; + const auto result = std::from_chars(value.data(), value.data() + value.size(), code); + if (result.ec != std::errc{} || result.ptr != value.data() + value.size()) { + RooJSONFactoryWSTool::error(context + " has invalid legacy interpolation code '" + value + "'"); + } + return code; +} + +Interpolation interpolationFromPiecewiseCode(int code, const std::string &context) +{ + switch (code) { + case 0: return additivePiecewiseLinear; + case 1: return multiplicativePiecewiseExponential; + case 2: + case 3: return additiveQuadraticLinear; + case 4: return additivePolynomialLinear; + case 5: return multiplicativePolynomialExponential; + case 6: return multiplicativePolynomialLinear; + default: + RooJSONFactoryWSTool::error(context + " has unsupported PiecewiseInterpolation code " + std::to_string(code)); + } +} + +Interpolation interpolationFromFlexibleCode(int code, const std::string &context) +{ + switch (code) { + case 0: return additivePiecewiseLinear; + case 1: return multiplicativePiecewiseExponential; + case 2: + case 3: return additiveQuadraticLinear; + case 4: + case 5: return multiplicativePolynomialExponential; + default: RooJSONFactoryWSTool::error(context + " has unsupported FlexibleInterpVar code " + std::to_string(code)); + } +} + +int piecewiseCodeFromInterpolation(const Interpolation &interpolation, const std::string &context) +{ + if (interpolation == additivePiecewiseLinear) + return 0; + if (interpolation == multiplicativePiecewiseExponential) + return 1; + if (interpolation == additiveQuadraticLinear) + return 2; + if (interpolation == additivePolynomialLinear) + return 4; + if (interpolation == multiplicativePolynomialExponential) + return 5; + if (interpolation == multiplicativePolynomialLinear) + return 6; + RooJSONFactoryWSTool::error(context + " " + interpolationString(interpolation) + + " cannot be represented by PiecewiseInterpolation"); +} + +int flexibleCodeFromInterpolation(const Interpolation &interpolation, const std::string &context) +{ + if (interpolation == additivePiecewiseLinear) + return 0; + if (interpolation == multiplicativePiecewiseExponential) + return 1; + if (interpolation == additiveQuadraticLinear) + return 2; + if (interpolation == multiplicativePolynomialExponential) + return 4; + RooJSONFactoryWSTool::error(context + " " + interpolationString(interpolation) + + " cannot be represented by FlexibleInterpVar"); +} + +enum class InterpolationClass { + Piecewise, + Flexible +}; + +int interpolationCode(const JSONNode &modifier, const std::optional &defaultInterpolation, + InterpolationClass interpolationClass, const std::string &context) +{ + const auto toCode = [&](const Interpolation &interpolation) { + return interpolationClass == InterpolationClass::Piecewise + ? piecewiseCodeFromInterpolation(interpolation, context) + : flexibleCodeFromInterpolation(interpolation, context); + }; + + if (const auto *interpolationNode = modifier.find("interpolation")) { + if (interpolationNode->is_map()) { + return toCode(readInterpolation(*interpolationNode, context)); + } + const int legacyCode = readLegacyInterpolationCode(*interpolationNode, context); + const Interpolation interpolation = interpolationClass == InterpolationClass::Piecewise + ? interpolationFromPiecewiseCode(legacyCode, context) + : interpolationFromFlexibleCode(legacyCode, context); + return toCode(interpolation); + } + if (defaultInterpolation) { + return toCode(*defaultInterpolation); + } + + // Before structured interpolation was introduced, both modifier classes + // used the integer code 4 as their implicit default. The meaning of code 4 + // is class-dependent. + return 4; +} + void erasePrefix(std::string &str, std::string_view prefix) { if (startsWith(str, prefix)) { @@ -347,7 +555,7 @@ double constraintRelError(RooAbsPdf const &constraint, RooAbsArg const &gamma) bool importHistSample(RooJSONFactoryWSTool &tool, RooDataHist &dh, RooArgSet const &varlist, RooAbsArg const *mcStatObject, const std::string &fprefix, const JSONNode &p, - RooArgSet &constraints) + const std::optional &defaultInterpolation, RooArgSet &constraints) { RooWorkspace &ws = *tool.workspace(); @@ -382,6 +590,7 @@ bool importHistSample(RooJSONFactoryWSTool &tool, RooDataHist &dh, RooArgSet con RooArgList histNps; RooArgList histoLo; RooArgList histoHi; + std::vector histoInterp; int idx = 0; for (const auto &mod : p["modifiers"].children()) { @@ -408,17 +617,16 @@ bool importHistSample(RooJSONFactoryWSTool &tool, RooDataHist &dh, RooArgSet con auto &par = getOrCreate(ws, parname, 0., -5, 5); overall_nps.add(par); auto &data = mod["data"]; - int interp = 4; - if (mod.has_child("interpolation")) { - interp = mod["interpolation"].val_int(); - } + const std::string context = "interpolation for normsys modifier '" + sysname + "' in sample '" + + sampleName + "' of channel '" + channelName + "'"; + const int interp = interpolationCode(mod, defaultInterpolation, InterpolationClass::Flexible, context); double low = data["lo"].val_double(); double high = data["hi"].val_double(); // the below contains a a hack to cut off variations that go below 0 - // this is needed because with interpolation code 4, which is the default, interpolation is done in - // log-space. hence, values <= 0 result in NaN which propagate throughout the model and cause evaluations to - // fail if you know a nicer way to solve this, please go ahead and fix the lines below + // This is needed because FlexibleInterpVar code 4 interpolates in log-space. Hence, values <= 0 result in + // NaN, which propagates throughout the model and causes evaluations to fail. If you know a nicer way to + // solve this, please go ahead and fix the lines below. if (interp == 4 && low <= 0) low = std::numeric_limits::epsilon(); if (interp == 4 && high <= 0) @@ -442,6 +650,9 @@ bool importHistSample(RooJSONFactoryWSTool &tool, RooDataHist &dh, RooArgSet con histoHi.add(tool.wsEmplace( sysname + "High_" + prefixedName, varlist, RooJSONFactoryWSTool::readBinnedData(data["hi"], sysname + "High_" + prefixedName, varlist))); + const std::string context = "interpolation for histosys modifier '" + sysname + "' in sample '" + + sampleName + "' of channel '" + channelName + "'"; + histoInterp.push_back(interpolationCode(mod, defaultInterpolation, InterpolationClass::Piecewise, context)); constraints.add(getOrCreateConstraint(tool, mod, par, sampleName)); } else if (modtype == "shapesys" || modtype == "shapefactor") { std::string funcName = channelName + "_" + sysname + "_ShapeSys"; @@ -527,9 +738,9 @@ bool importHistSample(RooJSONFactoryWSTool &tool, RooDataHist &dh, RooArgSet con normElems.add(v); } if (!histNps.empty()) { - auto &v = tool.wsEmplace("histoSys_" + prefixedName, hf, histoLo, histoHi, histNps); + auto &v = tool.wsEmplace("histoSys_" + prefixedName, hf, histoLo, histoHi, histNps, + histoInterp); v.setPositiveDefinite(); - v.setAllInterpCodes(4); // default interpCode for HistFactory shapeElems.add(v); } else { shapeElems.add(hf); @@ -556,6 +767,11 @@ class HistFactoryImporter : public RooFit::JSONIO::Importer { } double statErrThresh = 0; std::string statErrType = "Poisson"; + std::optional defaultInterpolation; + if (p.has_child("default_interpolation")) { + defaultInterpolation = + readInterpolation(p["default_interpolation"], "default_interpolation of channel '" + name + "'"); + } if (p.has_child(::Literals::staterror)) { auto &staterr = p[::Literals::staterror]; if (staterr.has_child("relThreshold")) @@ -622,7 +838,8 @@ class HistFactoryImporter : public RooFit::JSONIO::Importer { RooArgList funcs; RooArgList coefs; for (const auto &comp : p["samples"].children()) { - importHistSample(*tool, *data[idx], observables, mcStatObject, fprefix, comp, constraints); + importHistSample(*tool, *data[idx], observables, mcStatObject, fprefix, comp, defaultInterpolation, + constraints); ++idx; std::string const &compName = RooJSONFactoryWSTool::name(comp); @@ -769,11 +986,11 @@ struct NormSys { RooAbsReal const *param = nullptr; double low = 1.; double high = 1.; - int interpolationCode = 4; + Interpolation interpolation = multiplicativePolynomialExponential; RooAbsPdf const *constraint = nullptr; NormSys() {}; - NormSys(const std::string &n, RooAbsReal *const p, double h, double l, int i, const RooAbsPdf *c) - : name(n), param(p), low(l), high(h), interpolationCode(i), constraint(c) + NormSys(const std::string &n, RooAbsReal *const p, double h, double l, Interpolation i, const RooAbsPdf *c) + : name(n), param(p), low(l), high(h), interpolation(std::move(i)), constraint(c) { } }; @@ -783,12 +1000,11 @@ struct HistoSys { RooAbsReal const *param = nullptr; std::vector low; std::vector high; - // Used to validate duplicate RooFit modifiers. This is intentionally not serialized until HS3 defines the - // structured histosys interpolation representation. - int interpolationCode = 4; + Interpolation interpolation = additivePolynomialLinear; RooAbsPdf const *constraint = nullptr; - HistoSys(const std::string &n, RooAbsReal *const p, RooHistFunc *l, RooHistFunc *h, int i, const RooAbsPdf *c) - : name(n), param(p), interpolationCode(i), constraint(c) + HistoSys(const std::string &n, RooAbsReal *const p, RooHistFunc *l, RooHistFunc *h, Interpolation i, + const RooAbsPdf *c) + : name(n), param(p), interpolation(std::move(i)), constraint(c) { low.assign(l->dataHist().weightArray(), l->dataHist().weightArray() + l->dataHist().numEntries()); high.assign(h->dataHist().weightArray(), h->dataHist().weightArray() + h->dataHist().numEntries()); @@ -913,8 +1129,8 @@ NormSys parseOverallModifierFormula(const std::string &s, RooFormulaVar *formula sys.low = -sign * p2->getVal(); } - // interpolation code 1 means linear, which is what we have here - sys.interpolationCode = 1; + // Preserve the legacy export behaviour for recognized explicit formulae. + sys.interpolation = multiplicativePiecewiseExponential; erasePrefix(sys.name, "alpha_"); } @@ -1045,8 +1261,11 @@ Channel readChannel(RooJSONFactoryWSTool *tool, const std::string &pdfname, cons if (!constraint && !var->isConstant()) { RooJSONFactoryWSTool::error("cannot find constraint for " + std::string(var->GetName())); } else { + const std::string context = "normsys modifier '" + sysname + "' in sample '" + sample.name + + "' of channel '" + channel.name + "'"; sample.normsys.emplace_back(sysname, var, fip->high()[i], fip->low()[i], - fip->interpolationCodes()[i], constraint); + interpolationFromFlexibleCode(fip->interpolationCodes()[i], context), + constraint); } } } else if (!pip && (pip = dynamic_cast(e))) { @@ -1124,7 +1343,11 @@ Channel readChannel(RooJSONFactoryWSTool *tool, const std::string &pdfname, cons if (!constraint && !var->isConstant()) { RooJSONFactoryWSTool::error("cannot find constraint for " + std::string(var->GetName())); } else { - sample.histosys.emplace_back(sysname, var, lo, hi, pip->interpolationCodes()[i], constraint); + const std::string context = "histosys modifier '" + sysname + "' in sample '" + sample.name + + "' of channel '" + channel.name + "'"; + sample.histosys.emplace_back(sysname, var, lo, hi, + interpolationFromPiecewiseCode(pip->interpolationCodes()[i], context), + constraint); } } } @@ -1221,13 +1444,13 @@ void warnDuplicateModifiersCombined(const Channel &channel, const Sample &sample // the piecewise-exponential code 1 (exact everywhere) and for the default code 4 (exact at the +-1 sigma anchors and in // the exponential extrapolation region). The linear-space codes (e.g. 0 and 2) would turn the product into a shape that // cannot be represented by a single normsys, so those must not be merged. -bool normSysSupportsMultiplicativeMerge(int interpolationCode) +bool normSysSupportsMultiplicativeMerge(const Interpolation &interpolation) { - return interpolationCode == 1 || interpolationCode == 4; + return interpolation == multiplicativePiecewiseExponential || interpolation == multiplicativePolynomialExponential; } // Combines runs of adjacent modifiers that share the same name (the container is sorted by name beforehand) into a -// single modifier. The shared metadata (constraint, parameter and interpolation code) must be identical across the +// single modifier. The shared metadata (constraint, parameter and interpolation behaviour) must be identical across the // duplicates; the type-specific `combine` callable performs the actual merge and any additional validation. template void mergeDuplicateModifiers(const Channel &channel, const Sample &sample, Modifiers &modifiers, std::string_view type, @@ -1251,8 +1474,8 @@ void mergeDuplicateModifiers(const Channel &channel, const Sample &sample, Modif if (!hasSameMetadata(merged.param, modifier.param)) { duplicateModifierError(channel, sample, type, merged.name, "parameter metadata differs"); } - if (merged.interpolationCode != modifier.interpolationCode) { - duplicateModifierError(channel, sample, type, merged.name, "interpolation codes differ"); + if (merged.interpolation != modifier.interpolation) { + duplicateModifierError(channel, sample, type, merged.name, "interpolation behaviours differ"); } combine(merged, modifier); } @@ -1270,9 +1493,9 @@ void mergeDuplicateModifiers(const Channel &channel, const Sample &sample, Modif void mergeDuplicateNormSys(const Channel &channel, Sample &sample) { mergeDuplicateModifiers(channel, sample, sample.normsys, "normsys", [&](NormSys &merged, const NormSys &modifier) { - if (!normSysSupportsMultiplicativeMerge(merged.interpolationCode)) { + if (!normSysSupportsMultiplicativeMerge(merged.interpolation)) { duplicateModifierError(channel, sample, "normsys", merged.name, - "multiplicative combination is only valid for log-space interpolation codes"); + "multiplicative combination is only valid for log-space interpolation"); } merged.low *= modifier.low; merged.high *= modifier.high; @@ -1282,21 +1505,24 @@ void mergeDuplicateNormSys(const Channel &channel, Sample &sample) void mergeDuplicateHistoSys(const Channel &channel, Sample &sample) { const std::size_t nBins = sample.hist.size(); - mergeDuplicateModifiers( - channel, sample, sample.histosys, "histosys", [&](HistoSys &merged, const HistoSys &modifier) { - if (merged.interpolationCode != 4) { - duplicateModifierError(channel, sample, "histosys", merged.name, - "non-default interpolation cannot currently be represented by the HS3 exporter"); - } - if (merged.low.size() != nBins || merged.high.size() != nBins || modifier.low.size() != nBins || - modifier.high.size() != nBins) { - duplicateModifierError(channel, sample, "histosys", merged.name, "histogram binning differs"); - } - for (std::size_t bin = 0; bin < nBins; ++bin) { - merged.low[bin] += modifier.low[bin] - sample.hist[bin]; - merged.high[bin] += modifier.high[bin] - sample.hist[bin]; - } - }); + mergeDuplicateModifiers(channel, sample, sample.histosys, "histosys", + [&](HistoSys &merged, const HistoSys &modifier) { + if (merged.interpolation != additivePolynomialLinear) { + duplicateModifierError( + channel, sample, "histosys", merged.name, + "this interpolation cannot currently be combined for duplicate histosys " + "modifiers"); + } + if (merged.low.size() != nBins || merged.high.size() != nBins || + modifier.low.size() != nBins || modifier.high.size() != nBins) { + duplicateModifierError(channel, sample, "histosys", merged.name, + "histogram binning differs"); + } + for (std::size_t bin = 0; bin < nBins; ++bin) { + merged.low[bin] += modifier.low[bin] - sample.hist[bin]; + merged.high[bin] += modifier.high[bin] - sample.hist[bin]; + } + }); } void ensureUniqueModifiers(const Channel &channel, const Sample &sample) @@ -1354,6 +1580,32 @@ void configureStatError(Channel &channel) } } +std::optional defaultInterpolation(const Channel &channel) +{ + std::map counts; + for (const auto &sample : channel.samples) { + for (const auto &modifier : sample.normsys) { + ++counts[modifier.interpolation]; + } + for (const auto &modifier : sample.histosys) { + ++counts[modifier.interpolation]; + } + } + if (counts.empty()) { + return std::nullopt; + } + + auto best = counts.begin(); + for (auto current = std::next(counts.begin()); current != counts.end(); ++current) { + if (current->second > best->second || + (current->second == best->second && current->first == multiplicativePolynomialExponential && + best->first != multiplicativePolynomialExponential)) { + best = current; + } + } + return best->first; +} + bool exportChannel(RooJSONFactoryWSTool *tool, const Channel &channel, JSONNode &elem) { // Write the constraint reference for any modifier that supports an @@ -1364,11 +1616,15 @@ bool exportChannel(RooJSONFactoryWSTool *tool, const Channel &channel, JSONNode } }; + elem["type"] << "histfactory_dist"; + const auto channelDefaultInterpolation = defaultInterpolation(channel); + if (channelDefaultInterpolation) { + writeInterpolation(elem["default_interpolation"], *channelDefaultInterpolation); + } + bool observablesWritten = false; for (const auto &sample : channel.samples) { - elem["type"] << "histfactory_dist"; - auto &s = RooJSONFactoryWSTool::appendNamedChild(elem["samples"], sample.name); auto &modifiers = s["modifiers"]; @@ -1392,8 +1648,8 @@ bool exportChannel(RooJSONFactoryWSTool *tool, const Channel &channel, JSONNode mod["name"] << sys.name; mod["type"] << "normsys"; mod["parameter"] << sys.param->GetName(); - if (sys.interpolationCode != 4) { - mod["interpolation"] << sys.interpolationCode; + if (!channelDefaultInterpolation || sys.interpolation != *channelDefaultInterpolation) { + writeInterpolation(mod["interpolation"], sys.interpolation); } writeConstraint(mod, sys); auto &data = mod["data"].set_map(); @@ -1407,6 +1663,9 @@ bool exportChannel(RooJSONFactoryWSTool *tool, const Channel &channel, JSONNode mod["name"] << sys.name; mod["type"] << "histosys"; mod["parameter"] << sys.param->GetName(); + if (!channelDefaultInterpolation || sys.interpolation != *channelDefaultInterpolation) { + writeInterpolation(mod["interpolation"], sys.interpolation); + } writeConstraint(mod, sys); auto &data = mod["data"].set_map(); if (channel.nBins != sys.low.size() || channel.nBins != sys.high.size()) { diff --git a/roofit/hs3/test/testRooFitHS3.cxx b/roofit/hs3/test/testRooFitHS3.cxx index 36ad1233fa4c3..003b9a0105875 100644 --- a/roofit/hs3/test/testRooFitHS3.cxx +++ b/roofit/hs3/test/testRooFitHS3.cxx @@ -41,6 +41,7 @@ #include #include +#include #include #include #include @@ -233,6 +234,57 @@ std::size_t countOccurrences(std::string_view haystack, std::string_view needle) return result; } +std::string makeHistFactoryJSON(std::string const &modifiers, std::string const &defaultInterpolation = "") +{ + const std::string defaultNode = + defaultInterpolation.empty() ? "" : "\"default_interpolation\":" + defaultInterpolation + ","; + return R"({ + "metadata": {"hs3_version": "0.2"}, + "distributions": [ + { + "name": "model_channel0", + "type": "histfactory_dist", + "axes": [{"name": "x", "min": 0.0, "max": 2.0, "nbins": 2}],)" + + defaultNode + R"( + "samples": [ + { + "name": "sig", + "data": {"contents": [10.0, 20.0]}, + "modifiers": [)" + + modifiers + R"(] + } + ] + } + ] + })"; +} + +void expectInterpolation(const RooFit::Detail::JSONNode &node, std::string_view type, std::string_view in, + std::optional out) +{ + ASSERT_TRUE(node.is_map()); + ASSERT_TRUE(node.has_child("type")); + ASSERT_TRUE(node.has_child("in")); + ASSERT_TRUE(node.has_child("out")); + EXPECT_EQ(node["type"].val(), type); + EXPECT_EQ(node["in"].val(), in); + if (out) { + EXPECT_FALSE(node["out"].is_null()); + EXPECT_EQ(node["out"].val(), *out); + } else { + EXPECT_TRUE(node["out"].is_null()); + } +} + +const RooFit::Detail::JSONNode *findModifier(const RooFit::Detail::JSONNode &channel, std::string const &name) +{ + const auto *sample = RooJSONFactoryWSTool::findNamedChild(channel["samples"], "sig"); + if (!sample) { + return nullptr; + } + return RooJSONFactoryWSTool::findNamedChild((*sample)["modifiers"], name); +} + // Asserts that exporting `ws` to HS3 throws and logs an error message containing `expectedReason`. void expectExportThrowsWithError(RooWorkspace &ws, std::string const &expectedReason) { @@ -251,6 +303,22 @@ void expectExportThrowsWithError(RooWorkspace &ws, std::string const &expectedRe EXPECT_NE(errors.find(expectedReason), std::string::npos) << errors; } +void expectImportThrowsWithError(std::string const &json, std::string const &expectedReason) +{ + RooWorkspace ws; + bool threw = false; + const std::string errors = captureMessages(RooFit::ERROR, RooFit::IO, [&] { + try { + RooJSONFactoryWSTool{ws}.importJSONfromString(json); + } catch (const std::runtime_error &) { + threw = true; + } + }); + + EXPECT_TRUE(threw); + EXPECT_NE(errors.find(expectedReason), std::string::npos) << errors; +} + // Builds a single-bin-observable RooDataHist filled with the given bin contents. Returns an owning pointer so that it // can be handed to the RooHistFunc constructor that takes ownership (the const-reference overload keeps only an unowned // pointer, which would dangle for a temporary). @@ -1615,6 +1683,253 @@ TEST(RooFitHS3, UnbinnedDatasetAxisRange) EXPECT_EQ(axesNode.find("\"value\""), std::string::npos) << axesNode; } +TEST(RooFitHS3, HistFactoryInterpolationCodeMapping) +{ + struct TestCase { + const char *modifierType; + int inputCode; + const char *type; + const char *in; + const char *out; + int canonicalCode; + }; + const TestCase testCases[]{ + {"normsys", 0, "add", "poly1", nullptr, 0}, {"normsys", 1, "mult", "exp", nullptr, 1}, + {"normsys", 2, "add", "poly2", "poly1", 2}, {"normsys", 3, "add", "poly2", "poly1", 2}, + {"normsys", 4, "mult", "poly6", "exp", 4}, {"normsys", 5, "mult", "poly6", "exp", 4}, + {"histosys", 0, "add", "poly1", nullptr, 0}, {"histosys", 1, "mult", "exp", nullptr, 1}, + {"histosys", 2, "add", "poly2", "poly1", 2}, {"histosys", 3, "add", "poly2", "poly1", 2}, + {"histosys", 4, "add", "poly6", "poly1", 4}, {"histosys", 5, "mult", "poly6", "exp", 5}, + {"histosys", 6, "mult", "poly6", "poly1", 6}, + }; + + for (const auto &testCase : testCases) { + SCOPED_TRACE(std::string(testCase.modifierType) + " code " + std::to_string(testCase.inputCode)); + const bool isNormSys = std::string_view{testCase.modifierType} == "normsys"; + const std::string data = + isNormSys ? R"({"lo":0.8,"hi":1.2})" : R"({"lo":{"contents":[8.0,18.0]},"hi":{"contents":[12.0,22.0]}})"; + const std::string modifier = R"({"name":"modifier","parameter":"alpha_modifier","type":")" + + std::string(testCase.modifierType) + R"(","interpolation":)" + + std::to_string(testCase.inputCode) + R"(,"data":)" + data + "}"; + + RooWorkspace ws; + ASSERT_TRUE(RooJSONFactoryWSTool{ws}.importJSONfromString(makeHistFactoryJSON(modifier))); + const std::string exported = RooJSONFactoryWSTool{ws}.exportJSONtoString(); + auto tree = RooFit::Detail::JSONTree::create(exported); + const auto *channel = RooJSONFactoryWSTool::findNamedChild(tree->rootnode()["distributions"], "model_channel0"); + ASSERT_NE(channel, nullptr); + ASSERT_TRUE(channel->has_child("default_interpolation")); + expectInterpolation((*channel)["default_interpolation"], testCase.type, testCase.in, + testCase.out ? std::optional{testCase.out} : std::nullopt); + + const auto *exportedModifier = findModifier(*channel, "modifier"); + ASSERT_NE(exportedModifier, nullptr); + EXPECT_FALSE(exportedModifier->has_child("interpolation")); + + RooWorkspace roundTripped; + ASSERT_TRUE(RooJSONFactoryWSTool{roundTripped}.importJSONfromString(exported)); + if (isNormSys) { + auto *interpolation = + dynamic_cast(roundTripped.function("sig_channel0_epsilon")); + ASSERT_NE(interpolation, nullptr); + ASSERT_EQ(interpolation->interpolationCodes().size(), 1u); + EXPECT_EQ(interpolation->interpolationCodes()[0], testCase.canonicalCode); + } else { + auto *interpolation = + dynamic_cast(roundTripped.function("histoSys_model_channel0_sig")); + ASSERT_NE(interpolation, nullptr); + ASSERT_EQ(interpolation->interpolationCodes().size(), 1u); + EXPECT_EQ(interpolation->interpolationCodes()[0], testCase.canonicalCode); + } + } +} + +TEST(RooFitHS3, HistFactoryInterpolationDefaultsAndOverrides) +{ + const std::string modifiers = R"( + { + "name":"norm_default_1", "parameter":"alpha_norm_default_1", "type":"normsys", + "interpolation":4, "data":{"lo":0.8,"hi":1.2} + }, + { + "name":"norm_default_2", "parameter":"alpha_norm_default_2", "type":"normsys", + "interpolation":5, "data":{"lo":0.9,"hi":1.1} + }, + { + "name":"norm_linear", "parameter":"alpha_norm_linear", "type":"normsys", + "interpolation":0, "data":{"lo":0.85,"hi":1.15} + }, + { + "name":"shape_default", "parameter":"alpha_shape_default", "type":"histosys", + "interpolation":4, + "data":{"lo":{"contents":[8.0,18.0]},"hi":{"contents":[12.0,22.0]}} + }, + { + "name":"shape_quadratic", "parameter":"alpha_shape_quadratic", "type":"histosys", + "interpolation":2, + "data":{"lo":{"contents":[9.0,17.0]},"hi":{"contents":[11.0,23.0]}} + })"; + + RooWorkspace ws; + ASSERT_TRUE(RooJSONFactoryWSTool{ws}.importJSONfromString(makeHistFactoryJSON(modifiers))); + const std::string exported = RooJSONFactoryWSTool{ws}.exportJSONtoString(); + auto tree = RooFit::Detail::JSONTree::create(exported); + const auto *channel = RooJSONFactoryWSTool::findNamedChild(tree->rootnode()["distributions"], "model_channel0"); + ASSERT_NE(channel, nullptr); + + expectInterpolation((*channel)["default_interpolation"], "mult", "poly6", "exp"); + for (const char *name : {"norm_default_1", "norm_default_2"}) { + const auto *modifier = findModifier(*channel, name); + ASSERT_NE(modifier, nullptr); + EXPECT_FALSE(modifier->has_child("interpolation")); + } + + const auto *linear = findModifier(*channel, "norm_linear"); + ASSERT_NE(linear, nullptr); + ASSERT_TRUE(linear->has_child("interpolation")); + expectInterpolation((*linear)["interpolation"], "add", "poly1", std::nullopt); + + const auto *shapeDefault = findModifier(*channel, "shape_default"); + ASSERT_NE(shapeDefault, nullptr); + ASSERT_TRUE(shapeDefault->has_child("interpolation")); + expectInterpolation((*shapeDefault)["interpolation"], "add", "poly6", "poly1"); + + const auto *shapeQuadratic = findModifier(*channel, "shape_quadratic"); + ASSERT_NE(shapeQuadratic, nullptr); + ASSERT_TRUE(shapeQuadratic->has_child("interpolation")); + expectInterpolation((*shapeQuadratic)["interpolation"], "add", "poly2", "poly1"); + + RooWorkspace roundTripped; + ASSERT_TRUE(RooJSONFactoryWSTool{roundTripped}.importJSONfromString(exported)); + auto *shapeInterpolation = + dynamic_cast(roundTripped.function("histoSys_model_channel0_sig")); + ASSERT_NE(shapeInterpolation, nullptr); + ASSERT_EQ(shapeInterpolation->paramList().size(), 2u); + for (std::size_t i = 0; i < shapeInterpolation->paramList().size(); ++i) { + const std::string parameter = shapeInterpolation->paramList().at(i)->GetName(); + if (parameter == "alpha_shape_default") { + EXPECT_EQ(shapeInterpolation->interpolationCodes()[i], 4); + } else if (parameter == "alpha_shape_quadratic") { + EXPECT_EQ(shapeInterpolation->interpolationCodes()[i], 2); + } else { + FAIL() << "Unexpected histosys parameter " << parameter; + } + } + + const std::string implicitModifiers = R"( + {"name":"norm","parameter":"alpha_norm","type":"normsys","data":{"lo":0.8,"hi":1.2}}, + { + "name":"shape","parameter":"alpha_shape","type":"histosys", + "data":{"lo":{"contents":[8.0,18.0]},"hi":{"contents":[12.0,22.0]}} + })"; + RooWorkspace implicitWorkspace; + ASSERT_TRUE(RooJSONFactoryWSTool{implicitWorkspace}.importJSONfromString(makeHistFactoryJSON(implicitModifiers))); + const std::string implicitExported = RooJSONFactoryWSTool{implicitWorkspace}.exportJSONtoString(); + auto implicitTree = RooFit::Detail::JSONTree::create(implicitExported); + const auto *implicitChannel = + RooJSONFactoryWSTool::findNamedChild(implicitTree->rootnode()["distributions"], "model_channel0"); + ASSERT_NE(implicitChannel, nullptr); + expectInterpolation((*implicitChannel)["default_interpolation"], "mult", "poly6", "exp"); + const auto *implicitNorm = findModifier(*implicitChannel, "norm"); + const auto *implicitShape = findModifier(*implicitChannel, "shape"); + ASSERT_NE(implicitNorm, nullptr); + ASSERT_NE(implicitShape, nullptr); + EXPECT_FALSE(implicitNorm->has_child("interpolation")); + expectInterpolation((*implicitShape)["interpolation"], "add", "poly6", "poly1"); + + const std::string stableTieModifiers = R"( + {"name":"norm","parameter":"alpha_norm","type":"normsys","interpolation":0,"data":{"lo":0.8,"hi":1.2}}, + { + "name":"shape","parameter":"alpha_shape","type":"histosys","interpolation":2, + "data":{"lo":{"contents":[8.0,18.0]},"hi":{"contents":[12.0,22.0]}} + })"; + RooWorkspace stableTieWorkspace; + ASSERT_TRUE(RooJSONFactoryWSTool{stableTieWorkspace}.importJSONfromString(makeHistFactoryJSON(stableTieModifiers))); + auto stableTieTree = RooFit::Detail::JSONTree::create(RooJSONFactoryWSTool{stableTieWorkspace}.exportJSONtoString()); + const auto *stableTieChannel = + RooJSONFactoryWSTool::findNamedChild(stableTieTree->rootnode()["distributions"], "model_channel0"); + ASSERT_NE(stableTieChannel, nullptr); + expectInterpolation((*stableTieChannel)["default_interpolation"], "add", "poly1", std::nullopt); + + RooWorkspace noInterpolationWorkspace; + ASSERT_TRUE(RooJSONFactoryWSTool{noInterpolationWorkspace}.importJSONfromString( + makeHistFactoryJSON(R"({"name":"mu","parameter":"mu","type":"normfactor"})"))); + const std::string noInterpolationExported = RooJSONFactoryWSTool{noInterpolationWorkspace}.exportJSONtoString(); + auto noInterpolationTree = RooFit::Detail::JSONTree::create(noInterpolationExported); + const auto *noInterpolationChannel = + RooJSONFactoryWSTool::findNamedChild(noInterpolationTree->rootnode()["distributions"], "model_channel0"); + ASSERT_NE(noInterpolationChannel, nullptr); + EXPECT_FALSE(noInterpolationChannel->has_child("default_interpolation")); +} + +TEST(RooFitHS3, HistFactoryStructuredInterpolationRoundTrip) +{ + const std::string modifiers = R"( + { + "name":"norm", "parameter":"alpha_norm", "type":"normsys", + "data":{"lo":0.8,"hi":1.3} + }, + { + "name":"shape", "parameter":"alpha_shape", "type":"histosys", + "interpolation":{"type":"add","in":"poly2","out":"poly1"}, + "data":{"lo":{"contents":[8.0,17.0]},"hi":{"contents":[12.0,23.0]}} + })"; + const std::string defaultInterpolation = R"({"type":"mult","in":"poly6","out":"exp"})"; + + RooWorkspace before; + ASSERT_TRUE(RooJSONFactoryWSTool{before}.importJSONfromString(makeHistFactoryJSON(modifiers, defaultInterpolation))); + auto *normBefore = dynamic_cast(before.function("sig_channel0_epsilon")); + auto *shapeBefore = dynamic_cast(before.function("histoSys_model_channel0_sig")); + ASSERT_NE(normBefore, nullptr); + ASSERT_NE(shapeBefore, nullptr); + EXPECT_EQ(normBefore->interpolationCodes()[0], 4); + EXPECT_EQ(shapeBefore->interpolationCodes()[0], 2); + + const std::string exported = RooJSONFactoryWSTool{before}.exportJSONtoString(); + RooWorkspace after; + ASSERT_TRUE(RooJSONFactoryWSTool{after}.importJSONfromString(exported)); + + auto prediction = [](RooWorkspace &ws, int bin, double norm, double shape) { + ws.var("x")->setBin(bin); + ws.var("alpha_norm")->setVal(norm); + ws.var("alpha_shape")->setVal(shape); + RooArgSet observables{*ws.var("x")}; + return ws.function("model_channel0_sig_shapes")->getVal(observables) * + ws.function("model_channel0_sig_scaleFactors")->getVal(); + }; + + for (double theta : {-2.0, -1.0, -0.5, 0.0, 0.5, 1.0, 2.0}) { + for (int bin = 0; bin < 2; ++bin) { + const double expected = prediction(before, bin, theta, theta); + EXPECT_NEAR(prediction(after, bin, theta, theta), expected, std::abs(expected) * 1e-13) + << "theta=" << theta << ", bin=" << bin; + } + } +} + +TEST(RooFitHS3, HistFactoryStructuredInterpolationValidation) +{ + const std::string normsys = + R"({"name":"norm","parameter":"alpha_norm","type":"normsys","data":{"lo":0.8,"hi":1.2}})"; + + expectImportThrowsWithError(makeHistFactoryJSON(normsys, R"({"type":"mult","in":"poly6"})"), + "does not define the required 'out' component"); + expectImportThrowsWithError(makeHistFactoryJSON(normsys, R"({"type":"mult","in":"spline","out":null})"), + "unknown interpolation function 'spline'"); + expectImportThrowsWithError( + makeHistFactoryJSON( + R"({"name":"norm","parameter":"alpha_norm","type":"normsys","interpolation":{"type":"add","in":"poly6","out":"poly1"},"data":{"lo":0.8,"hi":1.2}})"), + "cannot be represented by FlexibleInterpVar"); + expectImportThrowsWithError( + makeHistFactoryJSON( + R"({"name":"shape","parameter":"alpha_shape","type":"histosys","interpolation":{"type":"add","in":"exp","out":null},"data":{"lo":{"contents":[8.0,18.0]},"hi":{"contents":[12.0,22.0]}}})"), + "cannot be represented by PiecewiseInterpolation"); + expectImportThrowsWithError( + makeHistFactoryJSON( + R"({"name":"norm","parameter":"alpha_norm","type":"normsys","interpolation":"four","data":{"lo":0.8,"hi":1.2}})"), + "invalid legacy interpolation code 'four'"); +} + // HistFactory channels with samples that have a zero-yield bin together with a // staterror modifier used to produce NaN gamma errors because the relative // error is computed as sqrt(sumW2)/sumW. Importing such a channel should now @@ -1687,6 +2002,8 @@ TEST(RooFitHS3, HistFactoryZeroYieldBin) TEST(RooFitHS3, HistFactoryDuplicateModifiersAreCombined) { + // The two "poly" modifiers deliberately use implicit code 4 and explicit + // code 5, which are semantic aliases in FlexibleInterpVar. const std::string jsonStr = R"({ "metadata": {"hs3_version": "0.1.90"}, "distributions": [ @@ -1725,6 +2042,7 @@ TEST(RooFitHS3, HistFactoryDuplicateModifiersAreCombined) "name": "poly", "parameter": "alpha_poly", "type": "normsys", + "interpolation": 5, "data": {"lo": 0.75, "hi": 1.25} }, { @@ -1930,7 +2248,7 @@ INSTANTIATE_TEST_SUITE_P( "name": "input_2", "parameter": "alpha_dup", "type": "normsys", "interpolation": 4, "data": {"lo": 0.9, "hi": 1.1} })", - "", "interpolation codes differ"}, + "", "interpolation behaviours differ"}, IncompatibleModifierCase{"Constraint", R"( { @@ -1976,7 +2294,7 @@ TEST(RooFitHS3, HistFactoryDuplicateHistoSysWithNonDefaultInterpolationFails) RooWorkspace ws{"ws_incompatible_histosys_interpolation"}; importDuplicateHistoSysModel(ws, /*interpCode=*/2, /*low2=*/{9.0, 16.0}, /*high2=*/{11.0, 24.0}); - expectExportThrowsWithError(ws, "non-default interpolation cannot currently be represented by the HS3 exporter"); + expectExportThrowsWithError(ws, "this interpolation cannot currently be combined for duplicate histosys modifiers"); } TEST(RooFitHS3, HistFactoryConstraintKeyMigration) From 74278ecd7fb808d759bbe61df791ca5fa4a623da Mon Sep 17 00:00:00 2001 From: Phmonski Date: Mon, 27 Jul 2026 11:00:48 +0200 Subject: [PATCH 2/6] [RF][HS3] Add interp. IO for FlexibleInterpVar + PiecewiseInterpolation --- roofit/hs3/src/JSONFactories_HistFactory.cxx | 159 ++++++++++++--- roofit/hs3/test/testRooFitHS3.cxx | 200 +++++++++++++++++++ 2 files changed, 331 insertions(+), 28 deletions(-) diff --git a/roofit/hs3/src/JSONFactories_HistFactory.cxx b/roofit/hs3/src/JSONFactories_HistFactory.cxx index 4630f4e6b6989..f7dd967c55c8c 100644 --- a/roofit/hs3/src/JSONFactories_HistFactory.cxx +++ b/roofit/hs3/src/JSONFactories_HistFactory.cxx @@ -242,13 +242,111 @@ enum class InterpolationClass { Flexible }; +Interpolation interpolationFromCode(int code, InterpolationClass interpolationClass, const std::string &context) +{ + return interpolationClass == InterpolationClass::Piecewise ? interpolationFromPiecewiseCode(code, context) + : interpolationFromFlexibleCode(code, context); +} + +int codeFromInterpolation(const Interpolation &interpolation, InterpolationClass interpolationClass, + const std::string &context) +{ + return interpolationClass == InterpolationClass::Piecewise ? piecewiseCodeFromInterpolation(interpolation, context) + : flexibleCodeFromInterpolation(interpolation, context); +} + +void writeInterpolations(JSONNode &node, const std::vector &codes, InterpolationClass interpolationClass, + const std::string &context) +{ + auto &interpolations = node.set_seq(); + if (codes.empty()) { + return; + } + + std::vector descriptors; + descriptors.reserve(codes.size()); + for (std::size_t i = 0; i < codes.size(); ++i) { + descriptors.push_back( + interpolationFromCode(codes[i], interpolationClass, context + " at parameter index " + std::to_string(i))); + } + + bool allEqual = true; + for (std::size_t i = 1; i < descriptors.size(); ++i) { + if (descriptors[i] != descriptors.front()) { + allEqual = false; + break; + } + } + + const std::size_t outputSize = allEqual ? 1 : descriptors.size(); + for (std::size_t i = 0; i < outputSize; ++i) { + writeInterpolation(interpolations.append_child(), descriptors[i]); + } +} + +std::vector readInterpolations(const JSONNode &object, std::size_t nParameters, + InterpolationClass interpolationClass, const std::string &context) +{ + if (const auto *interpolations = object.find("interpolations")) { + if (!interpolations->is_seq()) { + RooJSONFactoryWSTool::error(context + " component 'interpolations' must be an array"); + } + + const std::size_t size = interpolations->num_children(); + const bool validSize = nParameters == 0 ? size == 0 : size == 1 || size == nParameters; + if (!validSize) { + RooJSONFactoryWSTool::error(context + + " component 'interpolations' must contain either one descriptor or one " + "descriptor per parameter (got " + + std::to_string(size) + " for " + std::to_string(nParameters) + " parameters)"); + } + + std::vector codes; + codes.reserve(size); + std::size_t i = 0; + for (const auto &node : interpolations->children()) { + const std::string entryContext = context + " component 'interpolations' at index " + std::to_string(i); + codes.push_back( + codeFromInterpolation(readInterpolation(node, entryContext), interpolationClass, entryContext)); + ++i; + } + if (size == 1) { + codes.resize(nParameters, codes.front()); + } + return codes; + } + + std::vector codes(nParameters, 0); + if (const auto *legacyCodes = object.find("interpolationCodes")) { + if (!legacyCodes->is_seq()) { + RooJSONFactoryWSTool::error(context + " legacy component 'interpolationCodes' must be an array"); + } + if (legacyCodes->num_children() != nParameters) { + RooJSONFactoryWSTool::error(context + + " legacy component 'interpolationCodes' must contain one code per " + "parameter (got " + + std::to_string(legacyCodes->num_children()) + " for " + + std::to_string(nParameters) + " parameters)"); + } + + std::size_t i = 0; + for (const auto &node : legacyCodes->children()) { + const std::string entryContext = + context + " legacy component 'interpolationCodes' at index " + std::to_string(i); + const Interpolation interpolation = + interpolationFromCode(readLegacyInterpolationCode(node, entryContext), interpolationClass, entryContext); + codes[i] = codeFromInterpolation(interpolation, interpolationClass, entryContext); + ++i; + } + } + return codes; +} + int interpolationCode(const JSONNode &modifier, const std::optional &defaultInterpolation, InterpolationClass interpolationClass, const std::string &context) { const auto toCode = [&](const Interpolation &interpolation) { - return interpolationClass == InterpolationClass::Piecewise - ? piecewiseCodeFromInterpolation(interpolation, context) - : flexibleCodeFromInterpolation(interpolation, context); + return codeFromInterpolation(interpolation, interpolationClass, context); }; if (const auto *interpolationNode = modifier.find("interpolation")) { @@ -256,9 +354,7 @@ int interpolationCode(const JSONNode &modifier, const std::optional(func); + const std::size_t nParameters = fip->variables().size(); + if (fip->low().size() != nParameters || fip->high().size() != nParameters || + fip->interpolationCodes().size() != nParameters) { + RooJSONFactoryWSTool::error("FlexibleInterpVar '" + std::string{fip->GetName()} + + "' has non-matching parameter, variation, and interpolation lengths"); + } elem["type"] << key(); - elem["interpolationCodes"].fill_seq(fip->interpolationCodes()); + writeInterpolations(elem["interpolations"], fip->interpolationCodes(), InterpolationClass::Flexible, + "FlexibleInterpVar '" + std::string{fip->GetName()} + "'"); RooJSONFactoryWSTool::fillSeq(elem["vars"], fip->variables()); elem["nom"] << fip->nominal(); elem["high"].fill_seq(fip->high(), fip->variables().size()); @@ -890,8 +993,15 @@ class PiecewiseInterpolationStreamer : public RooFit::JSONIO::Exporter { bool exportObject(RooJSONFactoryWSTool *, const RooAbsArg *func, JSONNode &elem) const override { const PiecewiseInterpolation *pip = static_cast(func); + const std::size_t nParameters = pip->paramList().size(); + if (pip->lowList().size() != nParameters || pip->highList().size() != nParameters || + pip->interpolationCodes().size() != nParameters) { + RooJSONFactoryWSTool::error("PiecewiseInterpolation '" + std::string{pip->GetName()} + + "' has non-matching parameter, variation, and interpolation lengths"); + } elem["type"] << key(); - elem["interpolationCodes"].fill_seq(pip->interpolationCodes()); + writeInterpolations(elem["interpolations"], pip->interpolationCodes(), InterpolationClass::Piecewise, + "PiecewiseInterpolation '" + std::string{pip->GetName()} + "'"); elem["positiveDefinite"] << pip->positiveDefinite(); RooJSONFactoryWSTool::fillSeq(elem["vars"], pip->paramList()); elem["nom"] << pip->nominalHist()->GetName(); @@ -908,21 +1018,20 @@ class PiecewiseInterpolationFactory : public RooFit::JSONIO::Importer { std::string name(RooJSONFactoryWSTool::name(p)); RooArgList vars{tool->requestArgList(p, "vars")}; + RooArgList low{tool->requestArgList(p, "low")}; + RooArgList high{tool->requestArgList(p, "high")}; + if (vars.size() != low.size() || vars.size() != high.size()) { + RooJSONFactoryWSTool::error("PiecewiseInterpolation '" + name + + "' has non-matching lengths of 'vars', 'high' and 'low'"); + } + const std::vector codes = + readInterpolations(p, vars.size(), InterpolationClass::Piecewise, "PiecewiseInterpolation '" + name + "'"); - auto &pip = tool->wsEmplace(name, *tool->requestArg(p, "nom"), - tool->requestArgList(p, "low"), - tool->requestArgList(p, "high"), vars); + auto &pip = + tool->wsEmplace(name, *tool->requestArg(p, "nom"), low, high, vars, codes); pip.setPositiveDefinite(p["positiveDefinite"].val_bool()); - if (p.has_child("interpolationCodes")) { - std::size_t i = 0; - for (auto const &node : p["interpolationCodes"].children()) { - pip.setInterpCode(*static_cast(vars.at(i)), node.val_int(), true); - ++i; - } - } - return true; } }; @@ -956,16 +1065,10 @@ class FlexibleInterpVarFactory : public RooFit::JSONIO::Importer { RooJSONFactoryWSTool::error("FlexibleInterpVar '" + name + "' has non-matching lengths of 'vars', 'high' and 'low'!"); } + const std::vector codes = + readInterpolations(p, vars.size(), InterpolationClass::Flexible, "FlexibleInterpVar '" + name + "'"); - auto &fip = tool->wsEmplace(name, vars, nom, low, high); - - if (p.has_child("interpolationCodes")) { - size_t i = 0; - for (auto const &node : p["interpolationCodes"].children()) { - fip.setInterpCode(*static_cast(vars.at(i)), node.val_int()); - ++i; - } - } + tool->wsEmplace(name, vars, nom, low, high, codes); return true; } diff --git a/roofit/hs3/test/testRooFitHS3.cxx b/roofit/hs3/test/testRooFitHS3.cxx index 003b9a0105875..c1ae9e2906128 100644 --- a/roofit/hs3/test/testRooFitHS3.cxx +++ b/roofit/hs3/test/testRooFitHS3.cxx @@ -259,6 +259,29 @@ std::string makeHistFactoryJSON(std::string const &modifiers, std::string const })"; } +std::string makeStandaloneInterpolationJSON(std::string const &functions) +{ + return R"({ + "metadata": {"hs3_version": "0.2"}, + "parameter_points": [ + { + "name": "default_values", + "parameters": [ + {"name": "alpha", "value": 0.25}, + {"name": "beta", "value": -0.4}, + {"name": "nominal", "value": 10.0}, + {"name": "low1", "value": 8.0}, + {"name": "low2", "value": 7.0}, + {"name": "high1", "value": 12.0}, + {"name": "high2", "value": 14.0} + ] + } + ], + "functions": [)" + + functions + R"(] + })"; +} + void expectInterpolation(const RooFit::Detail::JSONNode &node, std::string_view type, std::string_view in, std::optional out) { @@ -1683,6 +1706,183 @@ TEST(RooFitHS3, UnbinnedDatasetAxisRange) EXPECT_EQ(axesNode.find("\"value\""), std::string::npos) << axesNode; } +TEST(RooFitHS3, StandaloneFlexibleInterpVarStructuredInterpolations) +{ + RooRealVar alpha{"alpha", "alpha", 0.25, -2.0, 2.0}; + RooRealVar beta{"beta", "beta", -0.4, -2.0, 2.0}; + RooArgList parameters{alpha, beta}; + RooStats::HistFactory::FlexibleInterpVar shared{"fiv_shared", "", parameters, 1.0, {0.8, 0.7}, {1.2, 1.4}, {4, 5}}; + RooStats::HistFactory::FlexibleInterpVar mixed{"fiv_mixed", "", parameters, 1.0, {0.8, 0.7}, {1.2, 1.4}, {0, 1}}; + + RooWorkspace source; + source.import(shared, RooFit::Silence()); + source.import(mixed, RooFit::RecycleConflictNodes(), RooFit::Silence()); + const std::string json = RooJSONFactoryWSTool{source}.exportJSONtoString(); + auto tree = RooFit::Detail::JSONTree::create(json); + + const auto *sharedNode = RooJSONFactoryWSTool::findNamedChild(tree->rootnode()["functions"], "fiv_shared"); + ASSERT_NE(sharedNode, nullptr); + EXPECT_FALSE(sharedNode->has_child("interpolationCodes")); + ASSERT_TRUE(sharedNode->has_child("interpolations")); + ASSERT_TRUE((*sharedNode)["interpolations"].is_seq()); + ASSERT_EQ((*sharedNode)["interpolations"].num_children(), 1u); + expectInterpolation((*sharedNode)["interpolations"].child(0), "mult", "poly6", "exp"); + + const auto *mixedNode = RooJSONFactoryWSTool::findNamedChild(tree->rootnode()["functions"], "fiv_mixed"); + ASSERT_NE(mixedNode, nullptr); + EXPECT_FALSE(mixedNode->has_child("interpolationCodes")); + ASSERT_TRUE(mixedNode->has_child("interpolations")); + ASSERT_EQ((*mixedNode)["interpolations"].num_children(), 2u); + expectInterpolation((*mixedNode)["interpolations"].child(0), "add", "poly1", std::nullopt); + expectInterpolation((*mixedNode)["interpolations"].child(1), "mult", "exp", std::nullopt); + + RooWorkspace imported; + ASSERT_TRUE(RooJSONFactoryWSTool{imported}.importJSONfromString(json)); + auto *importedShared = dynamic_cast(imported.function("fiv_shared")); + auto *importedMixed = dynamic_cast(imported.function("fiv_mixed")); + ASSERT_NE(importedShared, nullptr); + ASSERT_NE(importedMixed, nullptr); + EXPECT_EQ(importedShared->interpolationCodes(), (std::vector{4, 4})); + EXPECT_EQ(importedMixed->interpolationCodes(), (std::vector{0, 1})); + EXPECT_DOUBLE_EQ(importedShared->getVal(), source.function("fiv_shared")->getVal()); + EXPECT_DOUBLE_EQ(importedMixed->getVal(), source.function("fiv_mixed")->getVal()); +} + +TEST(RooFitHS3, StandalonePiecewiseInterpolationStructuredInterpolations) +{ + RooRealVar alpha{"alpha", "alpha", 0.25, -2.0, 2.0}; + RooRealVar beta{"beta", "beta", -0.4, -2.0, 2.0}; + RooRealVar nominal{"nominal", "nominal", 10.0}; + RooRealVar low1{"low1", "low1", 8.0}; + RooRealVar low2{"low2", "low2", 7.0}; + RooRealVar high1{"high1", "high1", 12.0}; + RooRealVar high2{"high2", "high2", 14.0}; + RooArgList lows{low1, low2}; + RooArgList highs{high1, high2}; + + PiecewiseInterpolation shared{"pip_shared", "", nominal, lows, highs, RooArgList{alpha, beta}, {4, 4}}; + PiecewiseInterpolation mixedRepeated{"pip_mixed_repeated", "", nominal, lows, highs, + RooArgList{alpha, alpha}, {0, 2}}; + mixedRepeated.setPositiveDefinite(); + + RooWorkspace source; + source.import(shared, RooFit::Silence()); + source.import(mixedRepeated, RooFit::RecycleConflictNodes(), RooFit::Silence()); + const std::string json = RooJSONFactoryWSTool{source}.exportJSONtoString(); + auto tree = RooFit::Detail::JSONTree::create(json); + + const auto *sharedNode = RooJSONFactoryWSTool::findNamedChild(tree->rootnode()["functions"], "pip_shared"); + ASSERT_NE(sharedNode, nullptr); + EXPECT_FALSE(sharedNode->has_child("interpolationCodes")); + ASSERT_TRUE(sharedNode->has_child("interpolations")); + ASSERT_EQ((*sharedNode)["interpolations"].num_children(), 1u); + expectInterpolation((*sharedNode)["interpolations"].child(0), "add", "poly6", "poly1"); + + const auto *mixedNode = RooJSONFactoryWSTool::findNamedChild(tree->rootnode()["functions"], "pip_mixed_repeated"); + ASSERT_NE(mixedNode, nullptr); + EXPECT_FALSE(mixedNode->has_child("interpolationCodes")); + ASSERT_TRUE(mixedNode->has_child("interpolations")); + ASSERT_EQ((*mixedNode)["interpolations"].num_children(), 2u); + expectInterpolation((*mixedNode)["interpolations"].child(0), "add", "poly1", std::nullopt); + expectInterpolation((*mixedNode)["interpolations"].child(1), "add", "poly2", "poly1"); + + RooWorkspace imported; + ASSERT_TRUE(RooJSONFactoryWSTool{imported}.importJSONfromString(json)); + auto *importedShared = dynamic_cast(imported.function("pip_shared")); + auto *importedMixed = dynamic_cast(imported.function("pip_mixed_repeated")); + ASSERT_NE(importedShared, nullptr); + ASSERT_NE(importedMixed, nullptr); + EXPECT_EQ(importedShared->interpolationCodes(), (std::vector{4, 4})); + EXPECT_EQ(importedMixed->interpolationCodes(), (std::vector{0, 2})); + ASSERT_EQ(importedMixed->paramList().size(), 2u); + EXPECT_STREQ(importedMixed->paramList().at(0)->GetName(), "alpha"); + EXPECT_STREQ(importedMixed->paramList().at(1)->GetName(), "alpha"); + EXPECT_TRUE(importedMixed->positiveDefinite()); + EXPECT_DOUBLE_EQ(importedShared->getVal(), source.function("pip_shared")->getVal()); + EXPECT_DOUBLE_EQ(importedMixed->getVal(), source.function("pip_mixed_repeated")->getVal()); +} + +TEST(RooFitHS3, StandaloneInterpolationLegacyAndDefaults) +{ + const std::string functions = R"( + { + "name":"legacy_fiv", "type":"interpolation0d", "vars":["alpha","beta"], "nom":1.0, + "low":[0.8,0.7], "high":[1.2,1.4], "interpolationCodes":[5,3] + }, + { + "name":"default_fiv", "type":"interpolation0d", "vars":["alpha","beta"], "nom":1.0, + "low":[0.8,0.7], "high":[1.2,1.4] + }, + { + "name":"precedence_fiv", "type":"interpolation0d", "vars":["alpha","beta"], "nom":1.0, + "low":[0.8,0.7], "high":[1.2,1.4], + "interpolations":[{"type":"add","in":"poly1","out":null}], "interpolationCodes":[4,4] + }, + { + "name":"legacy_pip", "type":"interpolation", "vars":["alpha","alpha"], "nom":"nominal", + "low":["low1","low2"], "high":["high1","high2"], "positiveDefinite":true, + "interpolationCodes":[0,2] + }, + { + "name":"default_pip", "type":"interpolation", "vars":["alpha","beta"], "nom":"nominal", + "low":["low1","low2"], "high":["high1","high2"], "positiveDefinite":false + })"; + + ScopedNoDomainConstVarImportFlag flagGuard{false}; + RooWorkspace imported; + ASSERT_TRUE(RooJSONFactoryWSTool{imported}.importJSONfromString(makeStandaloneInterpolationJSON(functions))); + + auto *legacyFiv = dynamic_cast(imported.function("legacy_fiv")); + auto *defaultFiv = dynamic_cast(imported.function("default_fiv")); + auto *precedenceFiv = dynamic_cast(imported.function("precedence_fiv")); + auto *legacyPip = dynamic_cast(imported.function("legacy_pip")); + auto *defaultPip = dynamic_cast(imported.function("default_pip")); + ASSERT_NE(legacyFiv, nullptr); + ASSERT_NE(defaultFiv, nullptr); + ASSERT_NE(precedenceFiv, nullptr); + ASSERT_NE(legacyPip, nullptr); + ASSERT_NE(defaultPip, nullptr); + EXPECT_EQ(legacyFiv->interpolationCodes(), (std::vector{4, 2})); + EXPECT_EQ(defaultFiv->interpolationCodes(), (std::vector{0, 0})); + EXPECT_EQ(precedenceFiv->interpolationCodes(), (std::vector{0, 0})); + EXPECT_EQ(legacyPip->interpolationCodes(), (std::vector{0, 2})); + EXPECT_EQ(defaultPip->interpolationCodes(), (std::vector{0, 0})); + EXPECT_TRUE(legacyPip->positiveDefinite()); + + const std::string exported = RooJSONFactoryWSTool{imported}.exportJSONtoString(); + EXPECT_EQ(exported.find("interpolationCodes"), std::string::npos) << exported; + auto tree = RooFit::Detail::JSONTree::create(exported); + const auto *legacyFivNode = RooJSONFactoryWSTool::findNamedChild(tree->rootnode()["functions"], "legacy_fiv"); + const auto *legacyPipNode = RooJSONFactoryWSTool::findNamedChild(tree->rootnode()["functions"], "legacy_pip"); + ASSERT_NE(legacyFivNode, nullptr); + ASSERT_NE(legacyPipNode, nullptr); + ASSERT_EQ((*legacyFivNode)["interpolations"].num_children(), 2u); + ASSERT_EQ((*legacyPipNode)["interpolations"].num_children(), 2u); +} + +TEST(RooFitHS3, StandaloneInterpolationRejectsMalformedLengths) +{ + ScopedNoDomainConstVarImportFlag flagGuard{false}; + const std::string descriptor = R"({"type":"add","in":"poly1","out":null})"; + + expectImportThrowsWithError( + makeStandaloneInterpolationJSON(R"({"name":"bad","type":"interpolation0d","vars":["alpha","beta"],"nom":1.0,)" + R"("low":[0.8,0.7],"high":[1.2,1.4],"interpolations":[]})"), + "component 'interpolations' must contain either one descriptor or one descriptor per parameter"); + + expectImportThrowsWithError( + makeStandaloneInterpolationJSON( + R"({"name":"bad","type":"interpolation","vars":["alpha","beta"],"nom":"nominal",)" + R"("low":["low1","low2"],"high":["high1","high2"],"positiveDefinite":false,"interpolations":[)" + + descriptor + "," + descriptor + "," + descriptor + "]}"), + "component 'interpolations' must contain either one descriptor or one descriptor per parameter"); + + expectImportThrowsWithError( + makeStandaloneInterpolationJSON(R"({"name":"bad","type":"interpolation0d","vars":["alpha","beta"],"nom":1.0,)" + R"("low":[0.8,0.7],"high":[1.2,1.4],"interpolationCodes":[4]})"), + "legacy component 'interpolationCodes' must contain one code per parameter"); +} + TEST(RooFitHS3, HistFactoryInterpolationCodeMapping) { struct TestCase { From 54d838539b3388418dde22c82b40b9844a7ea16f Mon Sep 17 00:00:00 2001 From: Phmonski Date: Mon, 27 Jul 2026 11:03:59 +0200 Subject: [PATCH 3/6] [RF][HS3] Apply changes suggested by clang-format --- roofit/hs3/src/JSONFactories_HistFactory.cxx | 34 +++++++++----------- 1 file changed, 16 insertions(+), 18 deletions(-) diff --git a/roofit/hs3/src/JSONFactories_HistFactory.cxx b/roofit/hs3/src/JSONFactories_HistFactory.cxx index f7dd967c55c8c..c9af587f8074a 100644 --- a/roofit/hs3/src/JSONFactories_HistFactory.cxx +++ b/roofit/hs3/src/JSONFactories_HistFactory.cxx @@ -1608,24 +1608,22 @@ void mergeDuplicateNormSys(const Channel &channel, Sample &sample) void mergeDuplicateHistoSys(const Channel &channel, Sample &sample) { const std::size_t nBins = sample.hist.size(); - mergeDuplicateModifiers(channel, sample, sample.histosys, "histosys", - [&](HistoSys &merged, const HistoSys &modifier) { - if (merged.interpolation != additivePolynomialLinear) { - duplicateModifierError( - channel, sample, "histosys", merged.name, - "this interpolation cannot currently be combined for duplicate histosys " - "modifiers"); - } - if (merged.low.size() != nBins || merged.high.size() != nBins || - modifier.low.size() != nBins || modifier.high.size() != nBins) { - duplicateModifierError(channel, sample, "histosys", merged.name, - "histogram binning differs"); - } - for (std::size_t bin = 0; bin < nBins; ++bin) { - merged.low[bin] += modifier.low[bin] - sample.hist[bin]; - merged.high[bin] += modifier.high[bin] - sample.hist[bin]; - } - }); + mergeDuplicateModifiers( + channel, sample, sample.histosys, "histosys", [&](HistoSys &merged, const HistoSys &modifier) { + if (merged.interpolation != additivePolynomialLinear) { + duplicateModifierError(channel, sample, "histosys", merged.name, + "this interpolation cannot currently be combined for duplicate histosys " + "modifiers"); + } + if (merged.low.size() != nBins || merged.high.size() != nBins || modifier.low.size() != nBins || + modifier.high.size() != nBins) { + duplicateModifierError(channel, sample, "histosys", merged.name, "histogram binning differs"); + } + for (std::size_t bin = 0; bin < nBins; ++bin) { + merged.low[bin] += modifier.low[bin] - sample.hist[bin]; + merged.high[bin] += modifier.high[bin] - sample.hist[bin]; + } + }); } void ensureUniqueModifiers(const Channel &channel, const Sample &sample) From 4fd7ed88ca68c5ade21dbc6af180afa605ae7fac Mon Sep 17 00:00:00 2001 From: Jonas Rembser Date: Wed, 29 Jul 2026 12:40:28 +0000 Subject: [PATCH 4/6] [RF][HS3] Consolidate HistFactory interpolation code mapping into a table The mapping between the structured HS3 interpolation descriptors and the RooFit integer codes was spread across four functions (two per direction, one pair per interpolation class) plus a pair of dispatchers. They all encoded the same table and had to be kept mutually consistent by hand. Replace them with a single static table that mirrors the documented mapping, plus two generic lookups. The class-specific quirks (legacy code 3 aliasing code 2, and FlexibleInterpVar treating code 4 and 5 identically) are now captured by two explicit normalization lines rather than being scattered across switch cases. No behavioural change; all HS3 interpolation tests still pass. --- roofit/hs3/src/JSONFactories_HistFactory.cxx | 130 +++++++++---------- 1 file changed, 63 insertions(+), 67 deletions(-) diff --git a/roofit/hs3/src/JSONFactories_HistFactory.cxx b/roofit/hs3/src/JSONFactories_HistFactory.cxx index c9af587f8074a..1627723fee8f1 100644 --- a/roofit/hs3/src/JSONFactories_HistFactory.cxx +++ b/roofit/hs3/src/JSONFactories_HistFactory.cxx @@ -177,82 +177,76 @@ int readLegacyInterpolationCode(const JSONNode &node, const std::string &context return code; } -Interpolation interpolationFromPiecewiseCode(int code, const std::string &context) -{ - switch (code) { - case 0: return additivePiecewiseLinear; - case 1: return multiplicativePiecewiseExponential; - case 2: - case 3: return additiveQuadraticLinear; - case 4: return additivePolynomialLinear; - case 5: return multiplicativePolynomialExponential; - case 6: return multiplicativePolynomialLinear; - default: - RooJSONFactoryWSTool::error(context + " has unsupported PiecewiseInterpolation code " + std::to_string(code)); - } -} +enum class InterpolationClass { + Piecewise, + Flexible +}; -Interpolation interpolationFromFlexibleCode(int code, const std::string &context) -{ - switch (code) { - case 0: return additivePiecewiseLinear; - case 1: return multiplicativePiecewiseExponential; - case 2: - case 3: return additiveQuadraticLinear; - case 4: - case 5: return multiplicativePolynomialExponential; - default: RooJSONFactoryWSTool::error(context + " has unsupported FlexibleInterpVar code " + std::to_string(code)); - } -} +// Single source of truth for the mapping between the structured HS3 interpolation +// descriptors and the RooFit integer codes of PiecewiseInterpolation and +// FlexibleInterpVar. A code of `kUnrepresentable` means the descriptor cannot be +// expressed by that class: FlexibleInterpVar internally remaps code 4 to code 5, +// so it has no additive-linear or multiplicative-linear poly6 variant. +struct InterpolationCodes { + const Interpolation &descriptor; + int piecewise; + int flexible; +}; + +constexpr int kUnrepresentable = -1; -int piecewiseCodeFromInterpolation(const Interpolation &interpolation, const std::string &context) +// clang-format off +const std::vector interpolationTable{ + {additivePiecewiseLinear, 0, 0}, + {multiplicativePiecewiseExponential, 1, 1}, + {additiveQuadraticLinear, 2, 2}, + {additivePolynomialLinear, 4, kUnrepresentable}, + {multiplicativePolynomialExponential, 5, 4}, + {multiplicativePolynomialLinear, 6, kUnrepresentable}, +}; +// clang-format on + +int codeForClass(const InterpolationCodes &row, InterpolationClass interpolationClass) { - if (interpolation == additivePiecewiseLinear) - return 0; - if (interpolation == multiplicativePiecewiseExponential) - return 1; - if (interpolation == additiveQuadraticLinear) - return 2; - if (interpolation == additivePolynomialLinear) - return 4; - if (interpolation == multiplicativePolynomialExponential) - return 5; - if (interpolation == multiplicativePolynomialLinear) - return 6; - RooJSONFactoryWSTool::error(context + " " + interpolationString(interpolation) + - " cannot be represented by PiecewiseInterpolation"); + return interpolationClass == InterpolationClass::Piecewise ? row.piecewise : row.flexible; } -int flexibleCodeFromInterpolation(const Interpolation &interpolation, const std::string &context) +const char *interpolationClassName(InterpolationClass interpolationClass) { - if (interpolation == additivePiecewiseLinear) - return 0; - if (interpolation == multiplicativePiecewiseExponential) - return 1; - if (interpolation == additiveQuadraticLinear) - return 2; - if (interpolation == multiplicativePolynomialExponential) - return 4; - RooJSONFactoryWSTool::error(context + " " + interpolationString(interpolation) + - " cannot be represented by FlexibleInterpVar"); + return interpolationClass == InterpolationClass::Piecewise ? "PiecewiseInterpolation" : "FlexibleInterpVar"; } -enum class InterpolationClass { - Piecewise, - Flexible -}; - Interpolation interpolationFromCode(int code, InterpolationClass interpolationClass, const std::string &context) { - return interpolationClass == InterpolationClass::Piecewise ? interpolationFromPiecewiseCode(code, context) - : interpolationFromFlexibleCode(code, context); + // Code 3 was historically an unimplemented alias of code 2 for both classes, + // and FlexibleInterpVar treats code 5 identically to its canonical code 4. + if (code == 3) { + code = 2; + } + if (interpolationClass == InterpolationClass::Flexible && code == 5) { + code = 4; + } + for (const auto &row : interpolationTable) { + const int rowCode = codeForClass(row, interpolationClass); + if (rowCode != kUnrepresentable && rowCode == code) { + return row.descriptor; + } + } + RooJSONFactoryWSTool::error(context + " has unsupported " + interpolationClassName(interpolationClass) + " code " + + std::to_string(code)); } int codeFromInterpolation(const Interpolation &interpolation, InterpolationClass interpolationClass, const std::string &context) { - return interpolationClass == InterpolationClass::Piecewise ? piecewiseCodeFromInterpolation(interpolation, context) - : flexibleCodeFromInterpolation(interpolation, context); + for (const auto &row : interpolationTable) { + const int rowCode = codeForClass(row, interpolationClass); + if (rowCode != kUnrepresentable && row.descriptor == interpolation) { + return rowCode; + } + } + RooJSONFactoryWSTool::error(context + " " + interpolationString(interpolation) + " cannot be represented by " + + interpolationClassName(interpolationClass)); } void writeInterpolations(JSONNode &node, const std::vector &codes, InterpolationClass interpolationClass, @@ -1366,9 +1360,10 @@ Channel readChannel(RooJSONFactoryWSTool *tool, const std::string &pdfname, cons } else { const std::string context = "normsys modifier '" + sysname + "' in sample '" + sample.name + "' of channel '" + channel.name + "'"; - sample.normsys.emplace_back(sysname, var, fip->high()[i], fip->low()[i], - interpolationFromFlexibleCode(fip->interpolationCodes()[i], context), - constraint); + sample.normsys.emplace_back( + sysname, var, fip->high()[i], fip->low()[i], + interpolationFromCode(fip->interpolationCodes()[i], InterpolationClass::Flexible, context), + constraint); } } } else if (!pip && (pip = dynamic_cast(e))) { @@ -1448,9 +1443,10 @@ Channel readChannel(RooJSONFactoryWSTool *tool, const std::string &pdfname, cons } else { const std::string context = "histosys modifier '" + sysname + "' in sample '" + sample.name + "' of channel '" + channel.name + "'"; - sample.histosys.emplace_back(sysname, var, lo, hi, - interpolationFromPiecewiseCode(pip->interpolationCodes()[i], context), - constraint); + sample.histosys.emplace_back( + sysname, var, lo, hi, + interpolationFromCode(pip->interpolationCodes()[i], InterpolationClass::Piecewise, context), + constraint); } } } From a12f24234268f2e2e79e1b290291f78d2f1fdc57 Mon Sep 17 00:00:00 2001 From: Jonas Rembser Date: Wed, 29 Jul 2026 12:44:46 +0000 Subject: [PATCH 5/6] [RF][HF] Remove dead commented-out code from PiecewiseInterpolation The analytical-integral routines carried several large commented-out blocks (old dirty-inhibit experiments, an abandoned numeric integral attempt, and two obsolete per-code interpolation loops), and the file ended with a fully commented-out printMetaArgs implementation. None of it is compiled; remove it. The git history preserves the old code for anyone who wants to revisit those ideas. No behavioural change. --- .../src/PiecewiseInterpolation.cxx | 191 ------------------ 1 file changed, 191 deletions(-) diff --git a/roofit/histfactory/src/PiecewiseInterpolation.cxx b/roofit/histfactory/src/PiecewiseInterpolation.cxx index fdf117d91f3f6..b78020fa8cd4a 100644 --- a/roofit/histfactory/src/PiecewiseInterpolation.cxx +++ b/roofit/histfactory/src/PiecewiseInterpolation.cxx @@ -283,18 +283,6 @@ bool PiecewiseInterpolation::setBinIntegrator(RooArgSet& allVars) Int_t PiecewiseInterpolation::getAnalyticalIntegralWN(RooArgSet& allVars, RooArgSet& analVars, const RooArgSet* normSet, const char* /*rangeName*/) const { - /* - std::cout << "---------------------------\nin PiecewiseInterpolation get analytic integral " << std::endl; - std::cout << "all vars = "<< std::endl; - allVars.Print("v"); - std::cout << "anal vars = "<< std::endl; - analVars.Print("v"); - std::cout << "normset vars = "<< std::endl; - if(normSet2) - normSet2->Print("v"); - */ - - // Handle trivial no-integration scenario if (allVars.empty()) return 0 ; if (_forceNumInt) return 0 ; @@ -356,75 +344,6 @@ Int_t PiecewiseInterpolation::getAnalyticalIntegralWN(RooArgSet& allVars, RooArg double PiecewiseInterpolation::analyticalIntegralWN(Int_t code, const RooArgSet* /*normSet2*/,const char* /*rangeName*/) const { - /* - std::cout <<"Enter analytic Integral"<< std::endl; - printDirty(true); - // _nominal.arg().setDirtyInhibit(true) ; - _nominal.arg().setShapeDirty() ; - RooAbsReal* temp ; - RooFIter lowIter(_lowSet.fwdIterator()) ; - while((temp=(RooAbsReal*)lowIter.next())) { - // temp->setDirtyInhibit(true) ; - temp->setShapeDirty() ; - } - RooFIter highIter(_highSet.fwdIterator()) ; - while((temp=(RooAbsReal*)highIter.next())) { - // temp->setDirtyInhibit(true) ; - temp->setShapeDirty() ; - } - */ - - /* - RooAbsArg::setDirtyInhibit(true); - printDirty(true); - std::cout <<"done setting dirty inhibit = true"<< std::endl; - - // old integral, only works for linear and not positive definite - CacheElem* cache = (CacheElem*) _normIntMgr.getObjByIndex(code-1) ; - - - std::unique_ptr vars2( getParameters(RooArgSet()) ); - std::unique_ptr iset( _normIntMgr.nameSet2ByIndex(code-1)->select(*vars2) ); - std::cout <<"iset = "<< std::endl; - iset->Print("v"); - - double sum = 0; - RooArgSet* vars = getVariables(); - vars->remove(_paramSet); - _paramSet.Print("v"); - vars->Print("v"); - if(vars->size()==1){ - RooRealVar* obs = (RooRealVar*) vars->first(); - for(int i=0; inumBins(); ++i){ - obs->setVal( obs->getMin() + (.5+i)*(obs->getMax()-obs->getMin())/obs->numBins()); - sum+=evaluate()*(obs->getMax()-obs->getMin())/obs->numBins(); - std::cout << "obs = " << obs->getVal() << " sum = " << sum << std::endl; - } - } else{ - std::cout <<"only know how to deal with 1 observable right now"<< std::endl; - } - */ - - /* - _nominal.arg().setDirtyInhibit(false) ; - RooFIter lowIter2(_lowSet.fwdIterator()) ; - while((temp=(RooAbsReal*)lowIter2.next())) { - temp->setDirtyInhibit(false) ; - } - RooFIter highIter2(_highSet.fwdIterator()) ; - while((temp=(RooAbsReal*)highIter2.next())) { - temp->setDirtyInhibit(false) ; - } - */ - - /* - RooAbsArg::setDirtyInhibit(false); - printDirty(true); - std::cout <<"done"<< std::endl; - std::cout << "sum = " <(_normIntMgr.getObjByIndex(code-1)) ; if( cache==nullptr ) { @@ -464,70 +383,6 @@ double PiecewiseInterpolation::analyticalIntegralWN(Int_t code, const RooArgSet* ++i; } - /* // MB : old bit of interpolation code - while( (param=(RooAbsReal*)_paramIter->Next()) ) { - low = (RooAbsReal*)lowIntIter->Next() ; - high = (RooAbsReal*)highIntIter->Next() ; - - if(param->getVal()>0) { - value += param->getVal()*(high->getVal() - nominal ); - } else { - value += param->getVal()*(nominal - low->getVal()); - } - ++i; - } - */ - - /* KC: the code below is wrong. Can't pull out a constant change to a non-linear shape deformation. - while( (param=(RooAbsReal*)paramIter.next()) ) { - low = (RooAbsReal*)lowIntIter.next() ; - high = (RooAbsReal*)highIntIter.next() ; - - if(_interpCode.empty() || _interpCode.at(i)==0){ - // piece-wise linear - if(param->getVal()>0) - value += param->getVal()*(high->getVal() - nominal ); - else - value += param->getVal()*(nominal - low->getVal()); - } else if(_interpCode.at(i)==1){ - // piece-wise log - if(param->getVal()>=0) - value *= pow(high->getVal()/nominal, +param->getVal()); - else - value *= pow(low->getVal()/nominal, -param->getVal()); - } else if(_interpCode.at(i)==2){ - // parabolic with linear - double a = 0.5*(high->getVal()+low->getVal())-nominal; - double b = 0.5*(high->getVal()-low->getVal()); - double c = 0; - if(param->getVal()>1 ){ - value += (2*a+b)*(param->getVal()-1)+high->getVal()-nominal; - } else if(param->getVal()<-1 ) { - value += -1*(2*a-b)*(param->getVal()+1)+low->getVal()-nominal; - } else { - value += a*pow(param->getVal(),2) + b*param->getVal()+c; - } - } else if(_interpCode.at(i)==3){ - //parabolic version of log-normal - double a = 0.5*(high->getVal()+low->getVal())-nominal; - double b = 0.5*(high->getVal()-low->getVal()); - double c = 0; - if(param->getVal()>1 ){ - value += (2*a+b)*(param->getVal()-1)+high->getVal()-nominal; - } else if(param->getVal()<-1 ) { - value += -1*(2*a-b)*(param->getVal()+1)+low->getVal()-nominal; - } else { - value += a*pow(param->getVal(),2) + b*param->getVal()+c; - } - - } else { - coutE(InputArguments) << "PiecewiseInterpolation::analyticalIntegralWN ERROR: " << param->GetName() - << " with unknown interpolation code" << std::endl ; - } - ++i; - } - */ - // std::cout << "value = " << value << std::endl; return value; } @@ -619,49 +474,3 @@ void PiecewiseInterpolation::Streamer(TBuffer &R__b) R__b.WriteClassBuffer(PiecewiseInterpolation::Class(),this); } } - - -/* -//////////////////////////////////////////////////////////////////////////////// -/// Customized printing of arguments of a PiecewiseInterpolation to more intuitively reflect the contents of the -/// product operator construction - -void PiecewiseInterpolation::printMetaArgs(ostream& os) const -{ - _lowIter->Reset() ; - if (_highIter) { - _highIter->Reset() ; - } - - bool first(true) ; - - RooAbsArg* arg1, *arg2 ; - if (_highSet.size()!=0) { - - while((arg1=(RooAbsArg*)_lowIter->Next())) { - if (!first) { - os << " + " ; - } else { - first = false ; - } - arg2=(RooAbsArg*)_highIter->Next() ; - os << arg1->GetName() << " * " << arg2->GetName() ; - } - - } else { - - while((arg1=(RooAbsArg*)_lowIter->Next())) { - if (!first) { - os << " + " ; - } else { - first = false ; - } - os << arg1->GetName() ; - } - - } - - os << " " ; -} - -*/ From 7d3dd58f5db764daa322991291e96d2d3189dd4f Mon Sep 17 00:00:00 2001 From: Jonas Rembser Date: Wed, 29 Jul 2026 12:51:20 +0000 Subject: [PATCH 6/6] [RF][HF] Share interpolation-code validation between the two classes FlexibleInterpVar and PiecewiseInterpolation had an almost identical setInterpCodeForParam(): the same range check, the same code-3 -> code-2 aliasing (with the same explanatory comment), and the same error messages, differing only in the highest supported code and the class name in the log output. Keeping the code-3 rule correct meant editing it in two places in lockstep. Extract the shared logic into a small internal helper, Detail::setInterpolationCode(), used by both classes. Each override now just forwards its parameter list and maximum code and flags itself value-dirty when a code was stored. No behavioural change; the error messages are preserved verbatim. --- roofit/histfactory/src/FlexibleInterpVar.cxx | 21 ++------ .../src/HistFactoryInterpolationCodeUtils.h | 53 +++++++++++++++++++ .../src/PiecewiseInterpolation.cxx | 21 ++------ 3 files changed, 63 insertions(+), 32 deletions(-) create mode 100644 roofit/histfactory/src/HistFactoryInterpolationCodeUtils.h diff --git a/roofit/histfactory/src/FlexibleInterpVar.cxx b/roofit/histfactory/src/FlexibleInterpVar.cxx index 254f517909fe5..9a8bbf01e6192 100644 --- a/roofit/histfactory/src/FlexibleInterpVar.cxx +++ b/roofit/histfactory/src/FlexibleInterpVar.cxx @@ -21,6 +21,8 @@ #include #include +#include "HistFactoryInterpolationCodeUtils.h" + #include #include @@ -133,23 +135,10 @@ void FlexibleInterpVar::setAllInterpCodes(int code) void FlexibleInterpVar::setInterpCodeForParam(int iParam, int code) { - RooAbsArg const ¶m = _paramList[iParam]; - if (code < 0 || code > 5) { - coutE(InputArguments) << "FlexibleInterpVar::setInterpCode ERROR: " << param.GetName() - << " with unknown interpolation code " << code << ", keeping current code " - << _interpCode[iParam] << std::endl; - return; - } - if (code == 3) { - // In the past, code 3 was equivalent to code 2, which confused users. - // Now, we just say that code 3 doesn't exist and default to code 2 in - // that case for backwards compatible behavior. - coutE(InputArguments) << "FlexibleInterpVar::setInterpCode ERROR: " << param.GetName() - << " with unknown interpolation code " << code << ", defaulting to code 2" << std::endl; - code = 2; + if (Detail::setInterpolationCode(*this, "FlexibleInterpVar", _paramList[iParam], _interpCode, iParam, code, + /*maxCode=*/5)) { + setValueDirty(); } - _interpCode.at(iParam) = code; - setValueDirty(); } //////////////////////////////////////////////////////////////////////////////// diff --git a/roofit/histfactory/src/HistFactoryInterpolationCodeUtils.h b/roofit/histfactory/src/HistFactoryInterpolationCodeUtils.h new file mode 100644 index 0000000000000..a0eb2461f1c11 --- /dev/null +++ b/roofit/histfactory/src/HistFactoryInterpolationCodeUtils.h @@ -0,0 +1,53 @@ +/// \cond ROOFIT_INTERNAL + +#ifndef ROOSTATS_HISTFACTORY_HISTFACTORYINTERPOLATIONCODEUTILS_H +#define ROOSTATS_HISTFACTORY_HISTFACTORYINTERPOLATIONCODEUTILS_H + +#include +#include + +#include +#include + +namespace RooStats { +namespace HistFactory { +namespace Detail { + +/// Validate an interpolation code and store it in `codes[iParam]`. +/// +/// Shared implementation for FlexibleInterpVar and PiecewiseInterpolation, which +/// only differ in the highest supported code. Codes outside [0, maxCode] are +/// rejected and the current code is kept. Code 3 is mapped to code 2 to preserve +/// the historical behaviour where the two were equivalent. `self` and `param` are +/// only used for the error messages. +/// +/// \return true if a code was stored, so the caller can flag itself value-dirty. +inline bool setInterpolationCode(RooAbsReal const &self, const char *className, RooAbsArg const ¶m, + std::vector &codes, std::size_t iParam, int code, int maxCode) +{ + if (code < 0 || code > maxCode) { + oocoutE(&self, InputArguments) << className << "::setInterpCode ERROR: " << param.GetName() + << " with unknown interpolation code " << code << ", keeping current code " + << codes[iParam] << std::endl; + return false; + } + if (code == 3) { + // In the past, code 3 was equivalent to code 2, which confused users. + // Now, we just say that code 3 doesn't exist and default to code 2 in + // that case for backwards compatible behavior. + oocoutE(&self, InputArguments) << className << "::setInterpCode ERROR: " << param.GetName() + << " with unknown interpolation code " << code << ", defaulting to code 2" + << std::endl; + code = 2; + } + codes.at(iParam) = code; + return true; +} + +} // namespace Detail +} // namespace HistFactory +} // namespace RooStats + +#endif + +/// \endcond diff --git a/roofit/histfactory/src/PiecewiseInterpolation.cxx b/roofit/histfactory/src/PiecewiseInterpolation.cxx index b78020fa8cd4a..9e4ae37a18e06 100644 --- a/roofit/histfactory/src/PiecewiseInterpolation.cxx +++ b/roofit/histfactory/src/PiecewiseInterpolation.cxx @@ -40,6 +40,8 @@ #include +#include "HistFactoryInterpolationCodeUtils.h" + #include "Riostream.h" #include "TBuffer.h" @@ -407,23 +409,10 @@ void PiecewiseInterpolation::setAllInterpCodes(int code) void PiecewiseInterpolation::setInterpCodeForParam(int iParam, int code) { - RooAbsArg const ¶m = _paramSet[iParam]; - if (code < 0 || code > 6) { - coutE(InputArguments) << "PiecewiseInterpolation::setInterpCode ERROR: " << param.GetName() - << " with unknown interpolation code " << code << ", keeping current code " - << _interpCode[iParam] << std::endl; - return; - } - if (code == 3) { - // In the past, code 3 was equivalent to code 2, which confused users. - // Now, we just say that code 3 doesn't exist and default to code 2 in - // that case for backwards compatible behavior. - coutE(InputArguments) << "PiecewiseInterpolation::setInterpCode ERROR: " << param.GetName() - << " with unknown interpolation code " << code << ", defaulting to code 2" << std::endl; - code = 2; + if (RooStats::HistFactory::Detail::setInterpolationCode(*this, "PiecewiseInterpolation", _paramSet[iParam], + _interpCode, iParam, code, /*maxCode=*/6)) { + setValueDirty(); } - _interpCode.at(iParam) = code; - setValueDirty(); } ////////////////////////////////////////////////////////////////////////////////