src/Utilities/AxisSpec.cpp

0001 // This file is part of the ACTS project.
0002 //
0003 // Copyright (C) 2016 CERN for the benefit of the ACTS project
0004 //
0005 // This Source Code Form is subject to the terms of the Mozilla Public
0006 // License, v. 2.0. If a copy of the MPL was not distributed with this
0007 // file, You can obtain one at https://mozilla.org/MPL/2.0/.
0008
0009 #include "Acts/Utilities/AxisSpec.hpp"
0010
0011 #include <algorithm>
0012 #include <ostream>
0013 #include <sstream>
0014 #include <stdexcept>
0015
0016 namespace Acts {
0017
0018 namespace {
0019
0020 /// Take the value from whichever side gives it, requiring both to agree if
0021 /// both do
0022 template <typename T>
0023 std::optional<T> mergeProperty(const std::optional<T>& specified,
0024                                const std::optional<T>& supplied,
0025                                const std::string& what) {
0026   if (specified.has_value() && supplied.has_value() &&
0027       *specified != *supplied) {
0028     throw std::invalid_argument("AxisSpec: the specified axis " + what +
0029                                 " does not match the one supplied");
0030   }
0031   return specified.has_value() ? specified : supplied;
0032 }
0033
0034 /// Unwrap a property that has to be known by now
0035 template <typename T>
0036 const T& requireProperty(const std::optional<T>& value,
0037                          const std::string& what) {
0038   if (!value.has_value()) {
0039     throw std::domain_error("AxisSpec: the axis " + what +
0040                             " is neither specified nor supplied");
0041   }
0042   return *value;
0043 }
0044
0045 /// Range of a spec, unset where it is left to the consumer
0046 std::optional<double> minOf(const AxisSpec::EquidistantParams& params) {
0047   return params.min;
0048 }
0049 std::optional<double> maxOf(const AxisSpec::EquidistantParams& params) {
0050   return params.max;
0051 }
0052 std::optional<double> minOf(const AxisSpec::VariableParams& params) {
0053   return params.edges.front();
0054 }
0055 std::optional<double> maxOf(const AxisSpec::VariableParams& params) {
0056   return params.edges.back();
0057 }
0058 std::optional<double> minOf(
0059     const AxisSpec::DeferredVariableParams& /*params*/) {
0060   return std::nullopt;
0061 }
0062 std::optional<double> maxOf(
0063     const AxisSpec::DeferredVariableParams& /*params*/) {
0064   return std::nullopt;
0065 }
0066
0067 void checkStrictlyIncreasing(const std::vector<double>& edges,
0068                              const std::string& context) {
0069   if (edges.size() < 2) {
0070     throw std::invalid_argument(context + ": at least two edges are required");
0071   }
0072   if (!std::ranges::is_sorted(edges, std::less_equal<double>())) {
0073     throw std::invalid_argument(context +
0074                                 ": edges must be strictly increasing");
0075   }
0076 }
0077
0078 }  // namespace
0079
0080 AxisSpec::AxisSpec(Variant variant,
0081                    std::optional<AxisBoundaryType> boundaryType,
0082                    std::optional<AxisDirection> direction)
0083     : m_variant(std::move(variant)),
0084       m_boundaryType(boundaryType),
0085       m_direction(direction) {}
0086
0087 AxisSpec AxisSpec::Equidistant(std::size_t nBins, std::optional<double> min,
0088                                std::optional<double> max,
0089                                std::optional<AxisBoundaryType> boundaryType,
0090                                std::optional<AxisDirection> direction) {
0091   if (min.has_value() && max.has_value() && *min >= *max) {
0092     throw std::invalid_argument("AxisSpec::Equidistant: min must be < max");
0093   }
0094   if (nBins == 0) {
0095     throw std::invalid_argument(
0096         "AxisSpec::Equidistant: at least one bin is required");
0097   }
0098   return AxisSpec(EquidistantParams{nBins, min, max}, boundaryType, direction);
0099 }
0100
0101 AxisSpec AxisSpec::DeferredEquidistant(std::size_t nBins,
0102                                        std::optional<AxisDirection> direction) {
0103   return Equidistant(nBins, std::nullopt, std::nullopt, std::nullopt,
0104                      direction);
0105 }
0106
0107 AxisSpec AxisSpec::Variable(std::vector<double> edges,
0108                             std::optional<AxisBoundaryType> boundaryType,
0109                             std::optional<AxisDirection> direction) {
0110   checkStrictlyIncreasing(edges, "AxisSpec::Variable");
0111   return AxisSpec(VariableParams{std::move(edges)}, boundaryType, direction);
0112 }
0113
0114 AxisSpec AxisSpec::DeferredVariable(
0115     std::vector<double> normalizedEdges,
0116     std::optional<AxisBoundaryType> boundaryType,
0117     std::optional<AxisDirection> direction) {
0118   checkStrictlyIncreasing(normalizedEdges, "AxisSpec::DeferredVariable");
0119   if (normalizedEdges.front() != 0. || normalizedEdges.back() != 1.) {
0120     throw std::invalid_argument(
0121         "AxisSpec::DeferredVariable: edges must be normalized to [0, 1]");
0122   }
0123   return AxisSpec(DeferredVariableParams{std::move(normalizedEdges)},
0124                   boundaryType, direction);
0125 }
0126
0127 AxisSpec AxisSpec::FromAxis(const IAxis& axis) {
0128   if (axis.getType() == AxisType::Equidistant) {
0129     return Equidistant(axis.getNBins(), axis.getMin(), axis.getMax(),
0130                        axis.getBoundaryType(), axis.getDirection());
0131   }
0132   return Variable(axis.getBinEdges(), axis.getBoundaryType(),
0133                   axis.getDirection());
0134 }
0135
0136 AxisSpec AxisSpec::withDirection(AxisDirection direction) const {
0137   return AxisSpec(m_variant, m_boundaryType, direction);
0138 }
0139
0140 AxisSpec AxisSpec::toDeferred() const {
0141   return std::visit(
0142       [this]<typename T>(const T& params) -> AxisSpec {
0143         if constexpr (std::is_same_v<T, EquidistantParams>) {
0144           return AxisSpec(EquidistantParams{params.nBins}, std::nullopt,
0145                           m_direction);
0146         } else if constexpr (std::is_same_v<T, VariableParams>) {
0147           std::vector<double> normalizedEdges = params.edges;
0148           double min = normalizedEdges.front();
0149           double max = normalizedEdges.back();
0150           for (double& edge : normalizedEdges) {
0151             edge = (edge - min) / (max - min);
0152           }
0153           // Force exact endpoints against floating point round-off
0154           normalizedEdges.front() = 0.;
0155           normalizedEdges.back() = 1.;
0156           return AxisSpec(DeferredVariableParams{std::move(normalizedEdges)},
0157                           std::nullopt, m_direction);
0158         } else {
0159           return AxisSpec(DeferredVariableParams{params.normalizedEdges},
0160                           std::nullopt, m_direction);
0161         }
0162       },
0163       m_variant);
0164 }
0165
0166 bool AxisSpec::isDeferred() const {
0167   if (!m_boundaryType.has_value()) {
0168     return true;
0169   }
0170   if (std::holds_alternative<DeferredVariableParams>(m_variant)) {
0171     return true;
0172   }
0173   if (const auto* eq = std::get_if<EquidistantParams>(&m_variant);
0174       eq != nullptr) {
0175     return !eq->min.has_value() || !eq->max.has_value();
0176   }
0177   return false;
0178 }
0179
0180 bool AxisSpec::isEquidistant() const {
0181   return std::holds_alternative<EquidistantParams>(m_variant);
0182 }
0183
0184 bool AxisSpec::isVariable() const {
0185   return !isEquidistant();
0186 }
0187
0188 bool AxisSpec::isDeferredVariable() const {
0189   return std::holds_alternative<DeferredVariableParams>(m_variant);
0190 }
0191
0192 std::optional<AxisDirection> AxisSpec::direction() const {
0193   return m_direction;
0194 }
0195
0196 std::optional<AxisBoundaryType> AxisSpec::boundaryType() const {
0197   return m_boundaryType;
0198 }
0199
0200 std::size_t AxisSpec::nBins() const {
0201   return std::visit(
0202       []<typename T>(const T& params) -> std::size_t {
0203         if constexpr (std::is_same_v<T, EquidistantParams>) {
0204           return params.nBins;
0205         } else if constexpr (std::is_same_v<T, VariableParams>) {
0206           return params.edges.size() - 1;
0207         } else {
0208           return params.normalizedEdges.size() - 1;
0209         }
0210       },
0211       m_variant);
0212 }
0213
0214 const AxisSpec::EquidistantParams& AxisSpec::asEquidistant() const {
0215   return std::get<EquidistantParams>(m_variant);
0216 }
0217
0218 const AxisSpec::VariableParams& AxisSpec::asVariable() const {
0219   return std::get<VariableParams>(m_variant);
0220 }
0221
0222 const AxisSpec::DeferredVariableParams& AxisSpec::asDeferredVariable() const {
0223   return std::get<DeferredVariableParams>(m_variant);
0224 }
0225
0226 std::unique_ptr<IAxis> AxisSpec::buildAxis(const Options& options) const {
0227   AxisBoundaryType boundaryType = requireProperty(
0228       mergeProperty(m_boundaryType, options.boundaryType, "boundary type"),
0229       "boundary type");
0230   std::optional<AxisDirection> direction =
0231       mergeProperty(m_direction, options.direction, "direction");
0232
0233   return std::visit(
0234       [&]<typename T>(const T& params) -> std::unique_ptr<IAxis> {
0235         // For absolute edges the range is implied by them, so merging only
0236         // validates a supplied one
0237         std::optional<double> min =
0238             mergeProperty(minOf(params), options.min, "minimum");
0239         std::optional<double> max =
0240             mergeProperty(maxOf(params), options.max, "maximum");
0241
0242         if constexpr (std::is_same_v<T, VariableParams>) {
0243           return IAxis::createVariable(boundaryType, params.edges, direction);
0244         } else {
0245           double minValue = requireProperty(min, "minimum");
0246           double maxValue = requireProperty(max, "maximum");
0247           if (minValue >= maxValue) {
0248             throw std::invalid_argument(
0249                 "AxisSpec::buildAxis: min must be < max");
0250           }
0251
0252           if constexpr (std::is_same_v<T, EquidistantParams>) {
0253             return IAxis::createEquidistant(boundaryType, minValue, maxValue,
0254                                             params.nBins, direction);
0255           } else {
0256             std::vector<double> edges = params.normalizedEdges;
0257             for (double& edge : edges) {
0258               edge = minValue + edge * (maxValue - minValue);
0259             }
0260             return IAxis::createVariable(boundaryType, edges, direction);
0261           }
0262         }
0263       },
0264       m_variant);
0265 }
0266
0267 std::string AxisSpec::toString() const {
0268   std::stringstream ss;
0269   ss << "AxisSpec: " << nBins() << " bins";
0270   if (m_direction.has_value()) {
0271     ss << " in " << axisDirectionName(*m_direction);
0272   }
0273   ss << (isEquidistant() ? ", equidistant" : ", variable");
0274   std::visit(
0275       [&ss](const auto& params) {
0276         std::optional<double> min = minOf(params);
0277         std::optional<double> max = maxOf(params);
0278         if (min.has_value() && max.has_value()) {
0279           ss << " within [" << *min << ", " << *max << "]";
0280         } else {
0281           ss << " within deferred range";
0282         }
0283       },
0284       m_variant);
0285   if (m_boundaryType.has_value()) {
0286     ss << ", " << *m_boundaryType;
0287   } else {
0288     ss << ", deferred boundary type";
0289   }
0290   return ss.str();
0291 }
0292
0293 }  // namespace Acts