ProfileFitter.h

Go to the documentation of this file.

//***************************************************************************************
//***************************************************************************************
 
#pragma once
 
#include "NAGASH/Tool.h"
#include "NAGASH/Global.h"
#include "NAGASH/ThreadPool.h"
 
#include "Minuit2/FCNBase.h"
#include "Minuit2/FCNGradientBase.h"
#include "Minuit2/FunctionMinimum.h"
#include "Minuit2/MnUserParameterState.h"
#include "Minuit2/MnPrint.h"
#include "Minuit2/MnMigrad.h"
#include "Minuit2/MnSimplex.h"
#include "Minuit2/MnHesse.h"
#include "Minuit2/MnMinos.h"
#include "Minuit2/MnScan.h"
#include "Minuit2/MnContours.h"
#include "Minuit2/MnPlot.h"
 
#include "Eigen/Dense"
 
namespace NAGASH
{
    class ProfileFitter : public Tool
    {
    public:
        enum class ParameterType
        {
            Nuisance,      
            Interest,      
            StatGamma,     
            Normalization, 
            Unknown
        };
        enum class FitMethod
        {
            ProfileLikelihood,                 
            AnalyticalChiSquare,               
            AnalyticalBoostedProfileLikelihood 
        };
 
        class Parameter
        {
        public:
            TString Name;
            double Value;
            double Error;
            double Min;
            double Max;
            double Init;
            ParameterType Type;
 
            double Nuisance_Gaus_Sigma = 1;
            double Nuisance_Gaus_Mean = 0;
 
            enum class ExtrapStrategy
            {
                FULLLINEAR,
                POLY6EXP
            };
 
            ExtrapStrategy Strategy = ExtrapStrategy::FULLLINEAR;
            bool isFitted = false;
            bool isFixed = false;
 
            double GetLikelihood();
            void SetStrategy(ExtrapStrategy str) { Strategy = str; }
            static double GetLikelihood_Thread(Parameter *tpar) { return tpar->GetLikelihood(); }
        };
 
        class ProfileCache
        {
        public:
            std::map<TString, double> Num_Map;
 
            double GetCache(const TString &name) { return Num_Map[name]; }
            void SaveCache(const TString &name, double val) { Num_Map.emplace(std::pair<TString, double>(name, val)); }
        };
 
        class Observable
        {
        public:
            TString Name;
            std::vector<double> Data;
            std::vector<double> Data_Stat;
            std::vector<double> Profile;
 
            class Sample
            {
            public:
                TString Name;
                std::vector<double> Nominal;
                std::vector<double> Nominal_Stat;
                std::vector<double> Profile;
 
                std::map<Parameter *, bool> NeedSymmetrization_Map;
                std::map<Parameter *, int> SmoothLevel_Map;
 
                std::map<Parameter *, std::vector<double>> Variation_Nominal_Map;
                std::map<Parameter *, std::vector<double>> Variation_Up_Map;
                std::map<Parameter *, std::vector<double>> Variation_Down_Map;
 
                std::map<Parameter *, std::vector<double>> Variation_Sigma_Map;
                std::map<Parameter *, std::vector<std::vector<double>>> Variation_Map;
 
                std::map<Parameter *, std::vector<std::shared_ptr<ProfileCache>>> Variation_Cache;
 
                std::vector<Parameter *> Gamma_Vec;
                Parameter *Norm_Par = nullptr;
 
                bool isnominal_booked = false;
 
                int Start_Index = -1;
                int End_Index = -1;
 
                int Rebin_Index = 1;
 
                int N_Variations = 0;
 
                void SortVariations();
                void GetProfile(std::map<TString, Parameter> &m_par);
                static NAGASH::StatusCode GetProfile_Thread(Sample *tsample, std::map<TString, Parameter> &m_par)
                {
                    tsample->GetProfile(m_par);
                    return NAGASH::StatusCode::SUCCESS;
                }
 
            private:
                double GetProfileBinVariation(int i, double aim_sigma, Parameter *par);
            };
 
            std::map<TString, Sample> Sample_Map;
 
            bool isdata_booked = false;
            bool isgamma_booked = false;
 
            int N_Samples = 0;
 
            int Start_Index = -1;
            int End_Index = -1;
 
            int Rebin_Index = 1;
 
            void SortVariations();
            double GetLikelihood();
            int GetNbins() { return Data.size(); }
            static double GetLikelihood_Thread(Observable *tobs)
            {
                for (int i = 0; i < tobs->Profile.size(); i++)
                {
                    double sum = 0;
                    for (auto &samp : tobs->Sample_Map)
                        sum = sum + samp.second.Profile[i];
                    tobs->Profile[i] = sum;
                }
                return tobs->GetLikelihood();
            }
            void GetProfile(std::map<TString, Parameter> &m_par);
 
        private:
            bool USE_POISSON = true;
        };
 
        class PLH_FCN : public ROOT::Minuit2::FCNBase
        {
        public:
            PLH_FCN(std::shared_ptr<ThreadPool> tpool,
                    std::shared_ptr<MSGTool> _msg_,
                    std::map<TString, Parameter> *par_map,
                    std::map<TString, Observable> *obs_map);
            // function value for negative log likelihood
            double operator()(const std::vector<double> &) const override;
            // error defination
            double Up() const { return 0.5; }
 
            std::shared_ptr<MSGTool> MSGUser() const { return msg; }
            void Reset()
            {
                NCall = 0;
                Max_Log_Likelihood = -999999;
            }
 
        private:
            std::map<TString, Parameter> *Par_Map;
            std::map<TString, Observable> *Obs_Map;
 
            std::shared_ptr<ThreadPool> pool;
            std::shared_ptr<MSGTool> msg;
 
            inline static std::recursive_mutex PrintMutex = std::recursive_mutex();
 
            mutable double Max_Log_Likelihood = -999999;
            mutable int NCall = 0;
        };
 
        ProfileFitter(std::shared_ptr<MSGTool> msg, int nthread = 1);
        void FixParameter(TString name, double val);
        void ReleaseParameter(TString name);
        void BookParameter(TString name, double init, double min, double max, ParameterType type);
        void BookParameterOfInterest(TString name, double init, double min, double max);
        void BookNuisanceParameter(TString name);
        void BookNuisanceParameter(TString name, TString gname);
        void BookNormalizationParameter(TString name, double init, double min, double max);
 
        void BookObservableVariation(TString name, TString samplename, const std::vector<double> &vec, TString parname, double sigma);
        void BookObservableVariationFromUnc(TString name, TString samplename, const std::vector<double> &vec, const std::vector<double> &vec_unc, TString parname);
        void BookObservableData(TString name, const std::vector<double> &vec, const std::vector<double> &vec_stat);
        void BookObservableNominal(TString name, TString samplename, const std::vector<double> &vec, const std::vector<double> &vec_stat);
 
        void BookObservableVariationNominal(TString name, TString samplename, const std::vector<double> &vec, TString parname)
        {
            BookObservableVariation(name, samplename, vec, parname, 0);
        }
        void BookObservableVariationUp(TString name, TString samplename, const std::vector<double> &vec, TString parname)
        {
            BookObservableVariation(name, samplename, vec, parname, 1);
        }
        void BookObservableVariationDown(TString name, TString samplename, const std::vector<double> &vec, TString parname)
        {
            BookObservableVariation(name, samplename, vec, parname, -1);
        }
 
        void SetGammaPrefix(TString samplename, TString prefix);
        void LinkNormalizationSample(TString name, TString samplename);
 
        void SetParameterStrategy(TString name, Parameter::ExtrapStrategy str) { GetParameter(name)->SetStrategy(str); }
 
        void SetSymmetrization(TString obsname, TString samplename, TString parname, bool need = true);
        void SymmetrizeVariation(std::vector<double> &up, std::vector<double> &down, std::vector<double> &norm);
 
        void SetSmooth(TString obsname, TString samplename, TString parname, int slevel);
        void SmoothVariation(std::vector<double> &var, std::vector<double> &norm, int level);
 
        void SetObservableRange(TString name, int start, int end);
        void RebinObservable(TString name, int index);
 
        Parameter *GetParameter(TString name);
        TString GetParameterName(int index);
        int GetNParameters(ParameterType type = ParameterType::Unknown);
        void SetParameterGroup(TString gname, TString parname);
        // Call This Function to Fit
        void Fit(FitMethod method);
 
        std::map<TString, double> EstimateUnc(const std::vector<TString> &par_vec, FitMethod method, int toy_num);
        std::map<TString, double> EstimateUnc(TString g_name, FitMethod method, int toy_num);
        // std::map<TString, double> EstimateUncParallel(TString g_name, FitMethod method, int toy_num);
 
        double GetChi2(const std::map<TString, double> &par_val);
        double GetChi2();
 
    private:
        std::map<TString, Parameter> Par_Map;
        std::map<TString, Observable> Obs_Map;
        std::map<TString, std::vector<TString>> Gamma_Map;
        std::map<TString, std::vector<TString>> Norm_Map;
        std::map<TString, std::vector<TString>> ParGroup_Map;
 
        std::shared_ptr<ThreadPool> pool = nullptr;
        int NTHREAD = 1;
 
        double MAX_NUISANCE_SIGMA = 3.0;
 
        Parameter::ExtrapStrategy POI_DEFAULT_STRATEGY = Parameter::ExtrapStrategy::FULLLINEAR;
        Parameter::ExtrapStrategy NUISANCE_DEFAULT_STRATEGY = Parameter::ExtrapStrategy::POLY6EXP;
 
        Observable *InitializeNewObservable(TString name);
        void PrintInfo();
        bool Check();
        void Prepare(FitMethod method);
 
        void FitPLH(); // Profile Likelihood (PLH)
        void FitACS(); // Analytical Chi-square (ACS)
 
        // Eigen Elements for ACS
        Eigen::MatrixXd Gamma_Matrix;
        Eigen::MatrixXd V_Matrix;
        Eigen::MatrixXd Q_Matrix;
        Eigen::MatrixXd S_Matrix;
        Eigen::MatrixXd C_Matrix;
        Eigen::MatrixXd H_Matrix;
        Eigen::MatrixXd Y0_Matrix;
        Eigen::MatrixXd Lambda_Matrix;
        Eigen::MatrixXd Rho_Matrix;
        Eigen::MatrixXd Epsilon_Matrix;
 
        Eigen::MatrixXd POI_Shift;
        Eigen::MatrixXd Nuisance_Shift;
 
        Eigen::MatrixXd POI_Cov_Matrix;
        Eigen::MatrixXd Nuisance_Cov_Matrix;
 
        void PrintResult();
    };
 
    inline void ProfileFitter::FixParameter(TString name, double val)
    {
        Parameter *tpar = GetParameter(name);
        tpar->Value = val;
        tpar->isFixed = true;
    }
 
    inline void ProfileFitter::ReleaseParameter(TString name)
    {
        Parameter *tpar = GetParameter(name);
        tpar->isFixed = false;
    }
 
    inline void ProfileFitter::BookNormalizationParameter(TString name, double init, double min, double max)
    {
        BookParameter(name, init, min, max, ParameterType::Normalization);
    }
 
    inline void ProfileFitter::BookParameterOfInterest(TString name, double init, double min, double max)
    {
        BookParameter(name, init, min, max, ParameterType::Interest);
    }
 
    inline void ProfileFitter::BookNuisanceParameter(TString name)
    {
        BookParameter(name, 0, -MAX_NUISANCE_SIGMA, MAX_NUISANCE_SIGMA, ParameterType::Nuisance);
    }
    inline void ProfileFitter::BookNuisanceParameter(TString name, TString gname)
    {
        BookNuisanceParameter(name);
        SetParameterGroup(gname, name);
    }
 
    inline ProfileFitter::PLH_FCN::PLH_FCN(std::shared_ptr<ThreadPool> tpool,
                                           std::shared_ptr<MSGTool> _msg_,
                                           std::map<TString, Parameter> *par_map,
                                           std::map<TString, Observable> *obs_map)
    {
        pool = tpool;
        msg = _msg_;
        Par_Map = par_map;
        Obs_Map = obs_map;
    }
}