BCIntegrate.h

#ifndef __BCINTEGRATE__H
#define __BCINTEGRATE__H

/*
 * Copyright (C) 2007-2018, the BAT core developer team
 * All rights reserved.
 *
 * For the licensing terms see doc/COPYING.
 * For documentation see http://mpp.mpg.de/bat
 */

// ---------------------------------------------------------

#include "BCEngineMCMC.h"

#include <TMinuitMinimizer.h>

#include <string>

// forward declarations
class BCIntegrate;

class TFile;
class TH1;
class TH2;
class TTree;

namespace BCCubaOptions
{

struct General {
    int ncomp, flags, nregions, neval, fail;
    double error, prob;

    General();
};

struct Vegas : public General {
    int nstart, nincrease, nbatch, gridno;

    Vegas();
};

struct Suave : public General {
    int nmin, nnew;
    double flatness;

    Suave();
};

struct Divonne : public General {
    int key1, key2, key3, maxpass;
    double border, maxchisq, mindeviation;

    Divonne();
};

struct Cuhre : public General {
    int key;

    Cuhre();
};
}

namespace BCMinimizer
{

class Adapter : public ROOT::Math::IMultiGenFunction
{
public:
    Adapter(BCIntegrate& m);
    virtual unsigned int NDim() const;
    virtual ROOT::Math::IMultiGenFunction* Clone() const;

    BCIntegrate* m;
    mutable std::vector<double> par;
private:
    virtual double DoEval (const double* x) const;
};

class Wrapper
{
public:
    Wrapper(BCIntegrate& m);
    void Init();
    void Init(const std::vector<double>& start, int printlevel);
    void Reset(BCIntegrate& m);

    TMinuitMinimizer min;
    Adapter adapt;
};
}

// ---------------------------------------------------------

class BCIntegrate : public BCEngineMCMC
{

public:

    enum BCOptimizationMethod {
        kOptEmpty,                                
        kOptSimAnn,                               
        kOptMetropolis,                           
        kOptMinuit,                               
        kOptDefault,                              
        NOptMethods                               
    };

    enum BCIntegrationMethod {
        kIntEmpty,                                
        kIntMonteCarlo,                           
        kIntCuba,                                 
        kIntGrid,                                 
        kIntLaplace,                              
        kIntDefault,                              
        NIntMethods                               
    };

    enum BCMarginalizationMethod {
        kMargEmpty,                               
        kMargMetropolis,                          
        kMargMonteCarlo,                          
        kMargGrid,                                
        kMargDefault,                             
        NMargMethods                              
    };

    enum BCSASchedule {
        kSACauchy,                                
        kSABoltzmann,                             
        kSACustom,                                
        NSAMethods                                
    };

    enum BCCubaMethod {
        kCubaVegas,                               
        kCubaSuave,                               
        kCubaDivonne,                             
        kCubaCuhre,                               
        kCubaDefault,                             
        NCubaMethods                              
    };

    typedef void (BCIntegrate::*tRandomizer)(std::vector<double>&) const;

    typedef double (BCIntegrate::*tEvaluator)(std::vector<double>&, const std::vector<double>&, bool&);

    typedef void (*tIntegralUpdater)(const std::vector<double>&, const int&, double&, double&);

    BCIntegrate(const std::string& name = "model");

    BCIntegrate(const std::string& filename, const std::string& name, bool loadObservables = true);

    BCIntegrate(const BCIntegrate& other);

    BCIntegrate& operator=(const BCIntegrate&);

    virtual ~BCIntegrate() {};

    double GetIntegral() const
    { return fIntegral; }

    BCIntegrate::BCOptimizationMethod GetOptimizationMethod() const
    { return fOptimizationMethodCurrent; }

    BCIntegrate::BCIntegrationMethod GetIntegrationMethod() const
    { return fIntegrationMethodCurrent; }

    BCIntegrate::BCMarginalizationMethod GetMarginalizationMethod() const
    { return fMarginalizationMethodCurrent; }

    BCIntegrate::BCSASchedule GetSASchedule() const
    { return fSASchedule; }

    void GetRandomVectorInParameterSpace(std::vector<double>& x) const;

    double GetRandomPoint(std::vector<double>& x);

    int GetNIterationsMin() const
    { return fNIterationsMin; }

    int GetNIterationsMax() const
    { return fNIterationsMax; }

    int GetNIterationsPrecisionCheck() const
    { return fNIterationsPrecisionCheck; }

    int GetNIterations() const
    { return fNIterations; }

    double GetRelativePrecision() const
    { return fRelativePrecision; }

    double GetAbsolutePrecision() const
    { return fAbsolutePrecision; }

    BCCubaMethod GetCubaIntegrationMethod() const
    { return fCubaIntegrationMethod; }

    const BCCubaOptions::Vegas& GetCubaVegasOptions() const
    { return fCubaVegasOptions; }

    const BCCubaOptions::Suave& GetCubaSuaveOptions() const
    { return fCubaSuaveOptions; }

    const BCCubaOptions::Divonne& GetCubaDivonneOptions() const
    { return fCubaDivonneOptions; }

    const BCCubaOptions::Cuhre& GetCubaCuhreOptions() const
    { return fCubaCuhreOptions; }

    TH1* GetSlice(std::vector<unsigned> indices, unsigned& nIterations, double& log_max_val, const std::vector<double> parameters = std::vector<double>(0), int nbins = 0, bool normalize = true);

    TH1* GetSlice(const std::string& name, unsigned& nIterations, double& log_max_val, const std::vector<double> parameters = std::vector<double>(0), int nbins = 0, bool normalize = true)
    { return GetSlice(fParameters.Index(name), nIterations, log_max_val, parameters, nbins, normalize); }

    TH1* GetSlice(unsigned index, unsigned& nIterations, double& log_max_val, const std::vector<double> parameters = std::vector<double>(0), int nbins = 0, bool normalize = true)
    { return GetSlice(std::vector<unsigned>(1, index), nIterations, log_max_val, parameters, nbins, normalize); }

    TH2* GetSlice(const std::string& name1, const std::string& name2, unsigned& nIterations, double& log_max_val, const std::vector<double> parameters = std::vector<double>(0), int nbins = 0, bool normalize = true)
    { return GetSlice(fParameters.Index(name1), fParameters.Index(name2), nIterations, log_max_val, parameters, nbins, normalize); }

    TH2* GetSlice(unsigned index1, unsigned index2, unsigned& nIterations, double& log_max_val, const std::vector<double> parameters = std::vector<double>(0), int nbins = 0, bool normalize = true);

    double GetError() const
    { return fError; }

    TMinuitMinimizer& GetMinuit()
    { return fMinimizer.min; }

    double GetSAT0() const
    { return fSAT0; }

    double GetSATmin() const
    { return fSATmin; }

    virtual const std::vector<double>& GetBestFitParameters() const;

    const std::vector<double>& GetBestFitParameterErrors() const;

    double GetLogMaximum() const
    { return fLogMaximum; };

    void SetFlagIgnorePrevOptimization(bool flag)
    { fFlagIgnorePrevOptimization = flag; }

    void SetOptimizationMethod(BCIntegrate::BCOptimizationMethod method)
    { fOptimizationMethodCurrent = method; }

    void SetIntegrationMethod(BCIntegrate::BCIntegrationMethod method);

    void SetMarginalizationMethod(BCIntegrate::BCMarginalizationMethod method)
    { fMarginalizationMethodCurrent = method; }

    void SetSASchedule(BCIntegrate::BCSASchedule schedule)
    { fSASchedule = schedule; }

    void SetNIterationsMin(int niterations)
    { fNIterationsMin = niterations; }

    void SetNIterationsMax(int niterations)
    { fNIterationsMax = niterations; }

    void SetNIterationsPrecisionCheck(int niterations)
    { fNIterationsPrecisionCheck = niterations; }

    void SetRelativePrecision(double relprecision)
    { fRelativePrecision = relprecision; }

    void SetAbsolutePrecision(double absprecision)
    { fAbsolutePrecision = absprecision; }

    void SetCubaIntegrationMethod(BCCubaMethod type);

    void SetCubaOptions(const BCCubaOptions::Vegas& options)
    { fCubaVegasOptions = options; }

    void SetCubaOptions(const BCCubaOptions::Suave& options)
    { fCubaSuaveOptions = options; }

    void SetCubaOptions(const BCCubaOptions::Divonne& options)
    { fCubaDivonneOptions = options; }

    void SetCubaOptions(const BCCubaOptions::Cuhre& options)
    { fCubaCuhreOptions = options; }

    void SetSAT0(double T0)
    { fSAT0 = T0; }

    void SetSATmin(double Tmin)
    { fSATmin = Tmin; }
    virtual double Eval(const std::vector<double>& x) = 0;

    virtual double LogEval(const std::vector<double>& x)
    { return log(Eval(x)); }

    double Normalize()
    { return Integrate(); };

    double Integrate(BCIntegrationMethod intmethod);

    double Integrate();

    double Integrate(BCIntegrationMethod type, tRandomizer randomizer, tEvaluator evaluator, tIntegralUpdater updater, std::vector<double>& sums);

    double EvaluatorMC(std::vector<double>& sums, const std::vector<double>& point, bool& accepted);

    static void IntegralUpdaterMC(const std::vector<double>& sums, const int& nIterations, double& integral, double& absprecision);

    int MarginalizeAll();

    int MarginalizeAll(BCMarginalizationMethod margmethod);

    virtual void MarginalizePreprocess()
    {};

    virtual void MarginalizePostprocess()
    {};

    std::vector<double> FindMode(std::vector<double> start = std::vector<double>());

    std::vector<double> FindMode(BCIntegrate::BCOptimizationMethod optmethod, std::vector<double> start = std::vector<double>());

    double SATemperature(double t) const;

    double SATemperatureBoltzmann(double t) const;

    double SATemperatureCauchy(double t) const;

    virtual double SATemperatureCustom(double t) const;

    std::vector<double> GetProposalPointSA(const std::vector<double>& x, int t) const;

    std::vector<double> GetProposalPointSABoltzmann(const std::vector<double>& x, int t) const;

    std::vector<double> GetProposalPointSACauchy(const std::vector<double>& x, int t) const;

    virtual std::vector<double> GetProposalPointSACustom(const std::vector<double>& x, int t) const;

    std::vector<double> SAHelperGetRandomPointOnHypersphere() const;

    double SAHelperGetRadialCauchy() const;

    double SAHelperSinusToNIntegral(int dim, double theta) const;

    virtual void ResetResults();

    std::string DumpIntegrationMethod(BCIntegrationMethod type) const;

    std::string DumpCurrentIntegrationMethod() const
    { return DumpIntegrationMethod(fIntegrationMethodCurrent); }

    std::string DumpUsedIntegrationMethod() const
    { return DumpIntegrationMethod(fIntegrationMethodUsed); }

    std::string DumpMarginalizationMethod(BCMarginalizationMethod type) const;

    std::string DumpCurrentMarginalizationMethod() const
    { return DumpMarginalizationMethod(fMarginalizationMethodCurrent); }

    std::string DumpUsedMarginalizationMethod() const
    { return DumpMarginalizationMethod(fMarginalizationMethodUsed); }

    std::string DumpOptimizationMethod(BCOptimizationMethod type) const;

    std::string DumpCurrentOptimizationMethod() const
    { return DumpOptimizationMethod(fOptimizationMethodCurrent); }

    std::string DumpUsedOptimizationMethod() const
    { return DumpOptimizationMethod(fOptimizationMethodUsed); }

    std::string DumpCubaIntegrationMethod(BCCubaMethod type) const;

    std::string DumpCubaIntegrationMethod() const
    { return DumpCubaIntegrationMethod(fCubaIntegrationMethod); }

    void SetBestFitParameters(const std::vector<double>& x);

    void SetBestFitParameters(const std::vector<double>& x, const double& new_value, double& old_value);

    bool CheckMarginalizationAvailability(BCMarginalizationMethod type);

    bool CheckMarginalizationIndices(TH1* hist, const std::vector<unsigned>& index);

    double IntegrateLaplace();

protected:
    unsigned IntegrationOutputFrequency() const;

    virtual void PrintBestFitSummary() const;

    virtual std::string GetBestFitSummary(unsigned i) const;

    virtual void PrintMarginalizationSummary() const;

    bool fFlagIgnorePrevOptimization;

    double fSAT0;

    double fSATmin;

    int fSANIterations;

    double fSATemperature;

    double fSALogProb;

    std::vector<double> fSAx;

    void LogOutputAtStartOfIntegration(BCIntegrationMethod type, BCCubaMethod cubatype);

    void LogOutputAtIntegrationStatusUpdate(BCIntegrationMethod type, double integral, double absprecision, int nIterations);

    void LogOutputAtEndOfIntegration(double integral, double absprecision, double relprecision, int nIterations);

    bool fFlagMarginalized;

private:

    BCMinimizer::Wrapper fMinimizer;

    std::vector<double> FindModeMinuit(std::vector<double>& mode, std::vector<double>& errors, std::vector<double> start = std::vector<double>(0), int printlevel = -1);

    std::vector<double> FindModeMCMC(std::vector<double>& mode, std::vector<double>& errors);

    std::vector<double> FindModeSA(std::vector<double>& mode, std::vector<double>& errors, std::vector<double> start = std::vector<double>(0));

    double IntegrateCuba()
    { return IntegrateCuba(fCubaIntegrationMethod); }

    double IntegrateCuba(BCCubaMethod cubatype);

    static int CubaIntegrand(const int* ndim, const double xx[], const int* ncomp, double ff[], void* userdata);

    double IntegrateSlice();

    BCIntegrate::BCOptimizationMethod fOptimizationMethodCurrent;

    BCIntegrate::BCOptimizationMethod fOptimizationMethodUsed;

    BCIntegrate::BCIntegrationMethod fIntegrationMethodCurrent;

    BCIntegrate::BCIntegrationMethod fIntegrationMethodUsed;

    BCIntegrate::BCMarginalizationMethod fMarginalizationMethodCurrent;

    BCIntegrate::BCMarginalizationMethod fMarginalizationMethodUsed;

    BCIntegrate::BCSASchedule fSASchedule;

    unsigned fNIterationsMin;

    unsigned fNIterationsMax;

    unsigned fNIterationsPrecisionCheck;

    int fNIterations;

    std::vector<double> fBestFitParameters;

    std::vector<double> fBestFitParameterErrors;

    double fLogMaximum;

    double fIntegral;

    double fRelativePrecision;

    double fAbsolutePrecision;

    double fError;

    BCCubaMethod fCubaIntegrationMethod;

    BCCubaOptions::Vegas fCubaVegasOptions;         
    BCCubaOptions::Suave fCubaSuaveOptions;         
    BCCubaOptions::Divonne fCubaDivonneOptions;     
    BCCubaOptions::Cuhre fCubaCuhreOptions;         
};

#endif