interface/bdt/LossFunctions.h

0001 // LossFunctions.h
0002 // Here we define the different loss functions that can be used
0003 // with the BDT system.
0004
0005 #ifndef L1Trigger_L1TMuonEndCap_emtf_LossFunctions
0006 #define L1Trigger_L1TMuonEndCap_emtf_LossFunctions
0007
0008 #include "Event.h"
0009 #include <string>
0010 #include <algorithm>
0011 #include <cmath>
0012
0013 // ========================================================
0014 // ================ Define the Interface ==================
0015 //=========================================================
0016
0017 namespace emtf {
0018
0019   // Define the Interface
0020   class LossFunction {
0021   public:
0022     // The gradient of the loss function.
0023     // Each tree is a step in the direction of the gradient
0024     // towards the minimum of the Loss Function.
0025     virtual double target(Event* e) = 0;
0026
0027     // The fit should minimize the loss function in each
0028     // terminal node at each iteration.
0029     virtual double fit(std::vector<Event*>& v) = 0;
0030     virtual std::string name() = 0;
0031     virtual int id() = 0;
0032     virtual ~LossFunction() = default;
0033   };
0034
0035   // ========================================================
0036   // ================ Least Squares =========================
0037   // ========================================================
0038
0039   class LeastSquares : public LossFunction {
0040   public:
0041     LeastSquares() {}
0042     ~LeastSquares() override {}
0043
0044     double target(Event* e) override {
0045       // Each tree fits the residuals when using LeastSquares.
0046       return e->trueValue - e->predictedValue;
0047     }
0048
0049     double fit(std::vector<Event*>& v) override {
0050       // The average of the residuals minmizes the Loss Function for LS.
0051
0052       double SUM = 0;
0053       for (unsigned int i = 0; i < v.size(); i++) {
0054         Event* e = v[i];
0055         SUM += e->trueValue - e->predictedValue;
0056       }
0057
0058       return SUM / v.size();
0059     }
0060     std::string name() override { return "Least_Squares"; }
0061     int id() override { return 1; }
0062   };
0063
0064   // ========================================================
0065   // ============== Absolute Deviation    ===================
0066   // ========================================================
0067
0068   class AbsoluteDeviation : public LossFunction {
0069   public:
0070     AbsoluteDeviation() {}
0071     ~AbsoluteDeviation() override {}
0072
0073     double target(Event* e) override {
0074       // The gradient.
0075       if ((e->trueValue - e->predictedValue) >= 0)
0076         return 1;
0077       else
0078         return -1;
0079     }
0080
0081     double fit(std::vector<Event*>& v) override {
0082       // The median of the residuals minimizes absolute deviation.
0083       if (v.empty())
0084         return 0;
0085       std::vector<double> residuals(v.size());
0086
0087       // Load the residuals into a vector.
0088       for (unsigned int i = 0; i < v.size(); i++) {
0089         Event* e = v[i];
0090         residuals[i] = (e->trueValue - e->predictedValue);
0091       }
0092
0093       // Get the median and return it.
0094       int median_loc = (residuals.size() - 1) / 2;
0095
0096       // Odd.
0097       if (residuals.size() % 2 != 0) {
0098         std::nth_element(residuals.begin(), residuals.begin() + median_loc, residuals.end());
0099         return residuals[median_loc];
0100       }
0101
0102       // Even.
0103       else {
0104         std::nth_element(residuals.begin(), residuals.begin() + median_loc, residuals.end());
0105         double low = residuals[median_loc];
0106         std::nth_element(residuals.begin() + median_loc + 1, residuals.begin() + median_loc + 1, residuals.end());
0107         double high = residuals[median_loc + 1];
0108         return (high + low) / 2;
0109       }
0110     }
0111     std::string name() override { return "Absolute_Deviation"; }
0112     int id() override { return 2; }
0113   };
0114
0115   // ========================================================
0116   // ============== Huber    ================================
0117   // ========================================================
0118
0119   class Huber : public LossFunction {
0120   public:
0121     Huber() {}
0122     ~Huber() override {}
0123
0124     double quantile;
0125     double residual_median;
0126
0127     double target(Event* e) override {
0128       // The gradient of the loss function.
0129
0130       if (std::abs(e->trueValue - e->predictedValue) <= quantile)
0131         return (e->trueValue - e->predictedValue);
0132       else
0133         return quantile * (((e->trueValue - e->predictedValue) > 0) ? 1.0 : -1.0);
0134     }
0135
0136     double fit(std::vector<Event*>& v) override {
0137       // The constant fit that minimizes Huber in a region.
0138
0139       quantile = calculateQuantile(v, 0.7);
0140       residual_median = calculateQuantile(v, 0.5);
0141
0142       double x = 0;
0143       for (unsigned int i = 0; i < v.size(); i++) {
0144         Event* e = v[i];
0145         double residual = e->trueValue - e->predictedValue;
0146         double diff = residual - residual_median;
0147         x += ((diff > 0) ? 1.0 : -1.0) * std::min(quantile, std::abs(diff));
0148       }
0149
0150       return (residual_median + x / v.size());
0151     }
0152
0153     std::string name() override { return "Huber"; }
0154     int id() override { return 3; }
0155
0156     double calculateQuantile(std::vector<Event*>& v, double whichQuantile) {
0157       // Container for the residuals.
0158       std::vector<double> residuals(v.size());
0159
0160       // Load the residuals into a vector.
0161       for (unsigned int i = 0; i < v.size(); i++) {
0162         Event* e = v[i];
0163         residuals[i] = std::abs(e->trueValue - e->predictedValue);
0164       }
0165
0166       std::sort(residuals.begin(), residuals.end());
0167       unsigned int quantile_location = whichQuantile * (residuals.size() - 1);
0168       return residuals[quantile_location];
0169     }
0170   };
0171
0172   // ========================================================
0173   // ============== Percent Error ===========================
0174   // ========================================================
0175
0176   class PercentErrorSquared : public LossFunction {
0177   public:
0178     PercentErrorSquared() {}
0179     ~PercentErrorSquared() override {}
0180
0181     double target(Event* e) override {
0182       // The gradient of the squared percent error.
0183       return (e->trueValue - e->predictedValue) / (e->trueValue * e->trueValue);
0184     }
0185
0186     double fit(std::vector<Event*>& v) override {
0187       // The average of the weighted residuals minimizes the squared percent error.
0188       // Weight(i) = 1/true(i)^2.
0189
0190       double SUMtop = 0;
0191       double SUMbottom = 0;
0192
0193       for (unsigned int i = 0; i < v.size(); i++) {
0194         Event* e = v[i];
0195         SUMtop += (e->trueValue - e->predictedValue) / (e->trueValue * e->trueValue);
0196         SUMbottom += 1 / (e->trueValue * e->trueValue);
0197       }
0198
0199       return SUMtop / SUMbottom;
0200     }
0201     std::string name() override { return "Percent_Error"; }
0202     int id() override { return 4; }
0203   };
0204
0205 }  // namespace emtf
0206
0207 #endif