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 a762c7328d61e..01f9687b60501 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, @@ -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); } @@ -1269,16 +1492,14 @@ 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)) { - duplicateModifierError( - channel, sample, "normsys", merged.name, - "multiplicative combination is only valid for log-space interpolation codes"); - } - merged.low *= modifier.low; - merged.high *= modifier.high; - }); + mergeDuplicateModifiers(channel, sample, sample.normsys, "normsys", [&](NormSys &merged, const NormSys &modifier) { + if (!normSysSupportsMultiplicativeMerge(merged.interpolation)) { + duplicateModifierError(channel, sample, "normsys", merged.name, + "multiplicative combination is only valid for log-space interpolation"); + } + merged.low *= modifier.low; + merged.high *= modifier.high; + }); } void mergeDuplicateHistoSys(const Channel &channel, Sample &sample) @@ -1286,10 +1507,11 @@ 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) { + if (merged.interpolation != additivePolynomialLinear) { duplicateModifierError( channel, sample, "histosys", merged.name, - "non-default interpolation cannot currently be represented by the HS3 exporter"); + "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) { @@ -1358,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 @@ -1368,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"]; @@ -1396,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(); @@ -1411,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 2f2a71b7835c5..78dafe999aacc 100644 --- a/roofit/hs3/test/testRooFitHS3.cxx +++ b/roofit/hs3/test/testRooFitHS3.cxx @@ -39,6 +39,7 @@ #include #include +#include #include #include @@ -212,6 +213,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) { @@ -230,6 +282,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). @@ -805,8 +873,7 @@ TEST(RooFitHS3, RooGenericPdf) TEST(RooFitHS3, GenericExpressionCleanup) { RooRealVar x{"x", "x", 0.5, -1.0, 1.0}; - RooFormulaVar formula{"formula", - "formula", + RooFormulaVar formula{"formula", "formula", "TMath::Floor(x) + TMath::Ceil(x) + TMath::Abs(x) + TMath::Tan(x) + " "TMath::ASin(x / 2.) + TMath::ACos(x / 2.) + TMath::ATan(x) + TMath::Pi() + TMath::E()", RooArgList{x}}; @@ -1329,6 +1396,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 @@ -1401,6 +1715,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": [ @@ -1439,6 +1755,7 @@ TEST(RooFitHS3, HistFactoryDuplicateModifiersAreCombined) "name": "poly", "parameter": "alpha_poly", "type": "normsys", + "interpolation": 5, "data": {"lo": 0.75, "hi": 1.25} }, { @@ -1470,9 +1787,8 @@ TEST(RooFitHS3, HistFactoryDuplicateModifiersAreCombined) ASSERT_TRUE(RooJSONFactoryWSTool{ws1}.importJSONfromString(jsonStr)); std::string exported; - const std::string warnings = captureMessages(RooFit::WARNING, RooFit::IO, [&] { - exported = RooJSONFactoryWSTool{ws1}.exportJSONtoString(); - }); + const std::string warnings = + captureMessages(RooFit::WARNING, RooFit::IO, [&] { exported = RooJSONFactoryWSTool{ws1}.exportJSONtoString(); }); EXPECT_NE(warnings.find("combined 2 duplicate modifiers named 'norm' of type 'normsys'"), std::string::npos) << warnings; @@ -1550,9 +1866,8 @@ TEST(RooFitHS3, HistFactoryDuplicateModifiersAreCombined) } // The operation must be logged as a warning, not as an error. - const std::string sentinelErrors = captureMessages(RooFit::ERROR, RooFit::IO, [] { - RooJSONFactoryWSTool::warning("duplicate-modifier warning-level sentinel"); - }); + const std::string sentinelErrors = captureMessages( + RooFit::ERROR, RooFit::IO, [] { RooJSONFactoryWSTool::warning("duplicate-modifier warning-level sentinel"); }); EXPECT_TRUE(sentinelErrors.empty()) << sentinelErrors; } @@ -1646,7 +1961,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"( { @@ -1692,7 +2007,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)