Phase2L1ParticleFlow/interface/conifer.h

0001 #ifndef L1TRIGGER_PHASE2L1PARTICLEFLOW_CONNIFER_H
0002 #define L1TRIGGER_PHASE2L1PARTICLEFLOW_CONNIFER_H
0003 #include "nlohmann/json.hpp"
0004 #include <fstream>
0005 #ifdef CMSSW_GIT_HASH
0006 #include "FWCore/Utilities/interface/Exception.h"
0007 #else
0008 #include <stdexcept>
0009 #endif
0010
0011 namespace conifer {
0012
0013   /* ---
0014 * Balanced tree reduce implementation.
0015 * Reduces an array of inputs to a single value using the template binary operator 'Op',
0016 * for example summing all elements with Op_add, or finding the maximum with Op_max
0017 * Use only when the input array is fully unrolled. Or, slice out a fully unrolled section
0018 * before applying and accumulate the result over the rolled dimension.
0019 * Required for emulation to guarantee equality of ordering.
0020 * --- */
0021   constexpr int floorlog2(int x) { return (x < 2) ? 0 : 1 + floorlog2(x / 2); }
0022
0023   template <int B>
0024   constexpr int pow(int x) {
0025     return x == 0 ? 1 : B * pow<B>(x - 1);
0026   }
0027
0028   constexpr int pow2(int x) { return pow<2>(x); }
0029
0030   template <class T, class Op>
0031   T reduce(std::vector<T> x, Op op) {
0032     int N = x.size();
0033     int leftN = pow2(floorlog2(N - 1)) > 0 ? pow2(floorlog2(N - 1)) : 0;
0034     //static constexpr int rightN = N - leftN > 0 ? N - leftN : 0;
0035     if (N == 1) {
0036       return x.at(0);
0037     } else if (N == 2) {
0038       return op(x.at(0), x.at(1));
0039     } else {
0040       std::vector<T> left(x.begin(), x.begin() + leftN);
0041       std::vector<T> right(x.begin() + leftN, x.end());
0042       return op(reduce<T, Op>(left, op), reduce<T, Op>(right, op));
0043     }
0044   }
0045
0046   template <class T>
0047   class OpAdd {
0048   public:
0049     T operator()(T a, T b) { return a + b; }
0050   };
0051
0052   template <class T, class U>
0053   class DecisionTree {
0054   private:
0055     std::vector<int> feature;
0056     std::vector<int> children_left;
0057     std::vector<int> children_right;
0058     std::vector<T> threshold_;
0059     std::vector<U> value_;
0060     std::vector<double> threshold;
0061     std::vector<double> value;
0062
0063   public:
0064     U decision_function(std::vector<T> x) const {
0065       /* Do the prediction */
0066       int i = 0;
0067       while (feature[i] != -2) {  // continue until reaching leaf
0068         bool comparison = x[feature[i]] <= threshold_[i];
0069         i = comparison ? children_left[i] : children_right[i];
0070       }
0071       return value_[i];
0072     }
0073
0074     void init_() {
0075       /* Since T, U types may not be readable from the JSON, read them to double and the cast them here */
0076       std::transform(
0077           threshold.begin(), threshold.end(), std::back_inserter(threshold_), [](double t) -> T { return (T)t; });
0078       std::transform(value.begin(), value.end(), std::back_inserter(value_), [](double v) -> U { return (U)v; });
0079     }
0080
0081     // Define how to read this class to/from JSON
0082     NLOHMANN_DEFINE_TYPE_INTRUSIVE(DecisionTree, feature, children_left, children_right, threshold, value);
0083
0084   };  // class DecisionTree
0085
0086   template <class T, class U, bool useAddTree = false>
0087   class BDT {
0088   private:
0089     unsigned int n_classes;
0090     unsigned int n_trees;
0091     unsigned int n_features;
0092     std::vector<double> init_predict;
0093     std::vector<U> init_predict_;
0094     // vector of decision trees: outer dimension tree, inner dimension class
0095     std::vector<std::vector<DecisionTree<T, U>>> trees;
0096     OpAdd<U> add;
0097
0098   public:
0099     // Define how to read this class to/from JSON
0100     NLOHMANN_DEFINE_TYPE_INTRUSIVE(BDT, n_classes, n_trees, n_features, init_predict, trees);
0101
0102     BDT(std::string filename) {
0103       /* Construct the BDT from conifer cpp backend JSON file */
0104       std::ifstream ifs(filename);
0105       nlohmann::json j = nlohmann::json::parse(ifs);
0106       from_json(j, *this);
0107       /* Do some transformation to initialise things into the proper emulation T, U types */
0108       if (n_classes == 2)
0109         n_classes = 1;
0110       std::transform(init_predict.begin(), init_predict.end(), std::back_inserter(init_predict_), [](double ip) -> U {
0111         return (U)ip;
0112       });
0113       for (unsigned int i = 0; i < n_trees; i++) {
0114         for (unsigned int j = 0; j < n_classes; j++) {
0115           trees.at(i).at(j).init_();
0116         }
0117       }
0118     }
0119
0120     std::vector<U> decision_function(std::vector<T> x) const {
0121       /* Do the prediction */
0122 #ifdef CMSSW_GIT_HASH
0123       if (x.size() != n_features) {
0124         throw cms::Exception("RuntimeError")
0125             << "Conifer : Size of feature vector mismatches expected n_features" << std::endl;
0126       }
0127 #else
0128       if (x.size() != n_features) {
0129         throw std::runtime_error("Conifer : Size of feature vector mismatches expected n_features");
0130       }
0131 #endif
0132       std::vector<U> values;
0133       std::vector<std::vector<U>> values_trees;
0134       values_trees.resize(n_classes);
0135       values.resize(n_classes, U(0));
0136       for (unsigned int i = 0; i < n_classes; i++) {
0137         std::transform(trees.begin(),
0138                        trees.end(),
0139                        std::back_inserter(values_trees.at(i)),
0140                        [&i, &x](std::vector<DecisionTree<T, U>> tree_v) { return tree_v.at(i).decision_function(x); });
0141         if (useAddTree) {
0142           values.at(i) = init_predict_.at(i);
0143           values.at(i) += reduce<U, OpAdd<U>>(values_trees.at(i), add);
0144         } else {
0145           values.at(i) = std::accumulate(values_trees.at(i).begin(), values_trees.at(i).end(), U(init_predict_.at(i)));
0146         }
0147       }
0148
0149       return values;
0150     }
0151
0152     std::vector<double> _decision_function_double(std::vector<double> x) const {
0153       /* Do the prediction with data in/out as double, cast to T, U before prediction */
0154       std::vector<T> xt;
0155       std::transform(x.begin(), x.end(), std::back_inserter(xt), [](double xi) -> T { return (T)xi; });
0156       std::vector<U> y = decision_function(xt);
0157       std::vector<double> yd;
0158       std::transform(y.begin(), y.end(), std::back_inserter(yd), [](U yi) -> double { return (double)yi; });
0159       return yd;
0160     }
0161
0162   };  // class BDT
0163
0164 }  // namespace conifer
0165
0166 #endif