LossFunctions.h

Go to the documentation of this file.

// LossFunctions.h
// Here we define the different loss functions that can be used
// with the BDT system. 

#ifndef L1Trigger_L1TMuonEndCap_L1TLossFunctions
#define L1Trigger_L1TMuonEndCap_L1TLossFunctions

#include "L1Trigger/L1TMuonEndCap/interface/Event.h"
#include <string>
#include <algorithm>

// ========================================================
// ================ Define the Interface ==================
//=========================================================

// Define the Interface
class L1TLossFunction
{
 public:
  
  // The gradient of the loss function.
  // Each tree is a step in the direction of the gradient
  // towards the minimum of the Loss Function.
  virtual Double_t target(emtf::Event* e) = 0;
  
  // The fit should minimize the loss function in each
  // terminal node at each iteration.
  virtual Double_t fit(std::vector<emtf::Event*>& v) = 0;
  virtual std::string name() = 0;
  virtual int id() = 0;
};

// ========================================================
// ================ Least Squares =========================
// ========================================================

class LeastSquares : public L1TLossFunction
{
 public:
  LeastSquares(){}
  ~LeastSquares(){}
  
  Double_t target(emtf::Event* e)
  {
    // Each tree fits the residuals when using LeastSquares.
    return e->trueValue - e->predictedValue;
  }
  
  Double_t fit(std::vector<emtf::Event*>& v)
  {
    // The average of the residuals minmizes the Loss Function for LS.
    
    Double_t SUM = 0;
    for(unsigned int i=0; i<v.size(); i++)
      {
        emtf::Event* e = v[i];
    SUM += e->trueValue - e->predictedValue;
      }
    
    return SUM/v.size();
  }
  std::string name() { return "Least_Squares"; }
  int id(){ return 1; }
  
};

// ========================================================
// ============== Absolute Deviation    ===================
// ========================================================

class AbsoluteDeviation : public L1TLossFunction
{
 public:
  AbsoluteDeviation(){}
  ~AbsoluteDeviation(){}
  
  Double_t target(emtf::Event* e)
  {
    // The gradient.
    if ((e->trueValue - e->predictedValue) >= 0)
      return 1;
    else
      return -1;
  }
  
  Double_t fit(std::vector<emtf::Event*>& v)
  {
    // The median of the residuals minimizes absolute deviation.
    if(v.size()==0) return 0;
    std::vector<Double_t> residuals(v.size());
    
    // Load the residuals into a vector. 
    for(unsigned int i=0; i<v.size(); i++)
      {
        emtf::Event* e = v[i];
    residuals[i] = (e->trueValue - e->predictedValue);
      }
    
    // Get the median and return it.
    int median_loc = (residuals.size()-1)/2;
    
    // Odd.
    if(residuals.size()%2 != 0)
      {
    std::nth_element(residuals.begin(), residuals.begin()+median_loc, residuals.end());
    return residuals[median_loc];
      }
    
    // Even.
    else
      {
    std::nth_element(residuals.begin(), residuals.begin()+median_loc, residuals.end());
    Double_t low = residuals[median_loc];
    std::nth_element(residuals.begin()+median_loc+1, residuals.begin()+median_loc+1, residuals.end());
    Double_t high = residuals[median_loc+1];
    return (high + low)/2;
      }
  }
  std::string name() { return "Absolute_Deviation"; }
  int id(){ return 2; }
};

// ========================================================
// ============== Huber    ================================
// ========================================================

class Huber : public L1TLossFunction
{
 public:
  Huber(){}
  ~Huber(){}
  
  double quantile;
  double residual_median;
  
  Double_t target(emtf::Event* e)
  {
    // The gradient of the loss function.
    
    if (TMath::Abs(e->trueValue - e->predictedValue) <= quantile)
      return (e->trueValue - e->predictedValue);
    else
      return quantile*(((e->trueValue - e->predictedValue) > 0)?1.0:-1.0);
  }
  
  Double_t fit(std::vector<emtf::Event*>& v)
  {
    // The constant fit that minimizes Huber in a region.
    
    quantile = calculateQuantile(v, 0.7);
    residual_median = calculateQuantile(v, 0.5); 
    
    double x = 0;
    for(unsigned int i=0; i<v.size(); i++)
      {
        emtf::Event* e = v[i];
    double residual = e->trueValue - e->predictedValue;
    double diff = residual - residual_median; 
    x += ((diff > 0)?1.0:-1.0)*std::min(quantile, TMath::Abs(diff));
      }
    
    return (residual_median + x/v.size());
    
  }
  
  std::string name() { return "Huber"; }
  int id(){ return 3; }
  
  double calculateQuantile(std::vector<emtf::Event*>& v, double whichQuantile)
  {
    // Container for the residuals.
    std::vector<Double_t> residuals(v.size());
    
    // Load the residuals into a vector. 
    for(unsigned int i=0; i<v.size(); i++)
      {
        emtf::Event* e = v[i];
    residuals[i] = TMath::Abs(e->trueValue - e->predictedValue);
      }
    
    std::sort(residuals.begin(), residuals.end());             
    unsigned int quantile_location = whichQuantile*(residuals.size()-1);
    return residuals[quantile_location];
  }        
};

// ========================================================
// ============== Percent Error ===========================
// ========================================================

class PercentErrorSquared : public L1TLossFunction
{
 public:
  PercentErrorSquared(){}
  ~PercentErrorSquared(){}
  
  Double_t target(emtf::Event* e)
  {   
    // The gradient of the squared percent error.
    return (e->trueValue - e->predictedValue)/(e->trueValue * e->trueValue);
  }   
  
  Double_t fit(std::vector<emtf::Event*>& v)
  {   
    // The average of the weighted residuals minimizes the squared percent error.
    // Weight(i) = 1/true(i)^2. 
    
    Double_t SUMtop = 0;
    Double_t SUMbottom = 0;
    
    for(unsigned int i=0; i<v.size(); i++)
      {   
        emtf::Event* e = v[i];
    SUMtop += (e->trueValue - e->predictedValue)/(e->trueValue*e->trueValue); 
    SUMbottom += 1/(e->trueValue*e->trueValue);
      }   
    
    return SUMtop/SUMbottom;
  }   
  std::string name() { return "Percent_Error"; }
  int id(){ return 4; }
};

#endif