StatisticManager.h

Go to the documentation of this file.

#ifndef STATISTIC_MANAGER_H
#define STATISTIC_MANAGER_H
 
#include <vector>
#include <iomanip>
#include <iostream>
#include <filesystem>
#include <optional>
#include <limits>
#include <numeric>
#include <set>
#include <random>
#include <map>
 
#include "Include.h"
#include "IModel.h"
#include "IStatCorrelationProxy.h"
#include "IStatParameterProxy.h"
#include "IStatSourcesProxy.h"
#include "IStatDependencyPruner.h"
#include "CovarianceTransformer.h"
#include "JointDistribution.h"
#include "RvgNuisanceSampler.h"
#include "MCEngine.h"
#include "MarginalFactory.h"
#include "CopulaFactory.h"
#include "Fit.h"
#include "ChiSquaredLikelihood.h"
#include "MarginalConfigFactory.h"
#include "INuisanceReader.h"
#include "NuisanceSpec.h"
#include "IStatParamOptimizerProxy.h"
#include "ChiSquaredLikelihood.h"
 
enum class StatisticLikelihoodMode {
    PROFILED_NUISANCE,  
    CHI2_MC_COVARIANCE  
};
 
struct AdvancedStatisticConfig {
    std::map<ParamId, MarginalType> override_nuisance_marginals {};      
    std::map<ExperimentObs, MarginalType> override_exp_data_marginals {};
    CopulaType nuisance_copula_type = CopulaType::GAUSSIAN;              
    CopulaType exp_data_copula_type = CopulaType::GAUSSIAN;              
 
    std::size_t MLE_max_iter = 500;                                      
    double MLE_tol = 1e-8;                                               
    unsigned MLE_strategy = 2;                                           
    bool MLE_run_hesse = true;                                           
    bool MLE_request_minos = false;                                      
    bool MLE_verbose = false;                                            
 
    double nuisance_relevance_cutoff = 1e-8;                             
    bool nuisance_sensitivity_pruning = true;                            
    double nuisance_sensitivity_probe_sigmas = 1.0;                      
    double nuisance_sensitivity_rel_cutoff = 1e-6;                       
    double nuisance_sensitivity_abs_cutoff = 1e-12;                      
    double nuisance_sensitivity_scale_floor = 1e-3;                      
    int nuisance_sensitivity_contexts = 2;                               
    double nuisance_sensitivity_context_sigma = 0.35;                    
    unsigned nuisance_sensitivity_seed = 12345;                          
    bool nuisance_sensitivity_keep_on_failure = true;                    
 
    bool fit_parameter_sensitivity_check = true;                         
    double fit_parameter_sensitivity_probe_fraction = 0.05;              
    double fit_parameter_sensitivity_rel_cutoff = 1e-10;                 
    double fit_parameter_sensitivity_abs_cutoff = 1e-12;                 
    bool fit_parameter_sensitivity_keep_on_failure = true;               
 
    bool MLE_trace_first_evals = false;                                  
    std::size_t MLE_trace_max_evals = 25;                                
    bool MLE_allow_profile_hessian_fallback = true;                      
    double MLE_profile_hessian_step_scale = 1.0;                         
    double MLE_profile_hessian_eig_floor_rel = 1e-8;                     
 
    StatisticLikelihoodMode likelihood_mode = StatisticLikelihoodMode::CHI2_MC_COVARIANCE; 
    double chi2_covariance_ridge_rel = 1e-8;                             
    double chi2_covariance_ridge_abs = 1e-12;                            
 
    bool MC_force_decay_threads_to_one = true;                           
    std::size_t MC_forced_decay_threads = 1;                             
};
 
struct StatisticConfig {
    std::size_t MC_draws = 100;                 
    std::size_t MC_threads = 1;                 
    unsigned int MC_seed = 123456u;             
    double skew_abs_threshold = 0.2;            
 
    bool print_mc_progress = true;              
    bool print_chi2_pipeline_progress = false;  
    bool print_mc_config = false;               
    bool print_fit_summary = false;             
    bool print_scan_summary = false;            
    bool print_cache_summary = false;           
    bool print_debug = false;                   
 
    bool write_mc_samples_csv = false;          
    std::string mc_samples_csv_path = "obs_samples.csv"; 
    std::size_t mc_progress_probe_draws = 5;    
    std::size_t mc_progress_update_every = 1;   
 
    std::shared_ptr<StatisticProgressMonitor> progress_monitor {}; 
 
    std::map<ParamId, std::pair<double, double>> fit_parameter_bounds {}; 
    std::map<ParamId, double> fit_parameter_offsets {}; 
 
    AdvancedStatisticConfig advanced {};        
};
 
struct LikelihoodScanPoint {
    double x = 0.0;             
    double y = 0.0;             
    double nll = 0.0;           
    double delta_nll = 0.0;     
};
 
struct LikelihoodScanGrid {
    ParamId x_param;                            
    ParamId y_param;                            
    double x_center = 0.0;                      
    double y_center = 0.0;                      
    std::size_t nx = 0;                         
    std::size_t ny = 0;                         
    std::vector<LikelihoodScanPoint> points;    
};
 
struct FitResultWithMaps {
    bool fit_ok {false};                                             
    std::map<ParamId, double> p_hat;                                 
    std::map<ParamId, double> eta_hat;                               
    std::map<ParamId, double> p_hat_std;                             
    std::map<ParamId, std::map<ParamId, double>> p_correlations;     
    double ell_hat {0.0};                                            
};
 
struct StatCache {
    std::map<ExperimentObs, double> exp_obs;                         
    std::map<ParamId, double> eta_specs_real;                        
    std::map<ParamId, double> p_specs;                               
    std::map<ParamId, std::map<ParamId, double>> SigmaEta;           
    std::map<ExperimentObs, std::map<ExperimentObs, double>> SigmaObs;
 
    FitResultWithMaps mle_result;                                    
};
 
class StatisticManager {
public:
    StatisticManager(StatisticConfig config,
                     std::shared_ptr<IModel> obs_int,
                     std::shared_ptr<IStatCorrelationProxy> pscp,
                     std::shared_ptr<IStatParameterProxy> pspp,
                     std::shared_ptr<IStatSourcesProxy> sp,
                     std::shared_ptr<IStatDependencyPruner> dp,
                     std::shared_ptr<INuisanceReader> nuisance_reader,
                     std::shared_ptr<IStatParamOptimizerProxy> spop);
    
    std::vector<std::unique_ptr<IMarginalDistribution>> build_nuisance_marginal_distributions();
 
    std::unique_ptr<JointDistribution> build_nuisance_distribution();
 
    std::unique_ptr<JointDistribution> build_exp_data_distribution();
 
    void reload_nuisance_specs();
 
    void set_nuisance_user_file(const fs::path& user_yaml_path);
 
    void clear_nuisance_user_file();
 
    const NuisanceRegistry& default_nuisance_specs() const { return default_nuisance_specs_; }
 
    const NuisanceRegistry& user_nuisance_specs() const { return user_nuisance_specs_; }
 
    const NuisanceRegistry& merged_nuisance_specs() const { return merged_nuisance_specs_; }
 
    std::map<BinnedObservableId, GaussianSummary> compute_uncertainties();
 
    MCResult compute_uncertainties_and_sampling();
    
    FitResultWithMaps compute_MLE(const std::vector<ParamId>& p_specs);
 
    Contour confidence_contour(ParamId p1, ParamId p2, double z, std::array<double, 4> bounds, ContourOptions options);
 
    void prepare_likelihood_for_scan(const std::vector<ParamId>& p_specs);
 
    void set_manual_scan_point(const std::map<ParamId, double>& p_hat,
                            const std::map<ParamId, double>& eta_hat);
 
    LikelihoodScanGrid scan_likelihood_around_current_point(
        ParamId p1,
        ParamId p2,
        double x_half_width,
        double y_half_width,
        std::size_t nx,
        std::size_t ny
    ) const;
 
    void save_likelihood_scan_csv(const std::string& path,
                                const LikelihoodScanGrid& grid) const;
 
    void update_cache(const std::vector<ParamId>& p_specs = std::vector<ParamId>());
 
    void select_experiment(const std::string& experiment);
 
    void select_experiments(const std::set<std::string>& experiments);
 
    void select_experiments(const std::vector<std::string>& experiments);
 
    void select_experiments_all();
 
    bool has_experiment_selection() const noexcept;
 
    std::set<std::string> selected_experiments() const;
 
    void select_experiment_observables(const std::set<ExperimentObs>& observables);
 
    void select_experiment_observables(const std::vector<ExperimentObs>& observables);
 
    void select_experiment_observables_all();
 
    bool has_experiment_observable_selection() const noexcept;
 
    std::set<ExperimentObs> selected_experiment_observables() const;
 
    std::map<ParamId, double> get_all_obss_deps();
 
    std::map<ParamId, double> get_p_specs(const std::vector<ParamId>& p_specs);
 
    std::map<ParamId, std::map<ParamId, double>> get_all_correlations();
 
    std::map<ExperimentObs, std::map<ExperimentObs, double>> get_all_obs_correlations();
 
    std::map<ExperimentObs, double> get_obs_exp();
 
    void print_cache();
 
private:
    void rebuild_merged_nuisance_specs();
 
    void invalidate_fit_state();
 
    const NuisanceSpec* find_nuisance_spec(const ParamId& pid) const;
 
    MarginalType resolve_nuisance_marginal_type(const ParamId& pid) const;
 
    fit_app::ParameterDefinition make_nuisance_parameter_definition(const ParamId& pid,
                                                                    double value,
                                                                    double sigma_hint) const;
 
    MarginalConfig make_nuisance_marginal_config(const ParamId& pid,
                                                 MarginalType mt) const;
 
    void validate_fit_parameter_sensitivity();
 
    bool accepts_experiment_observable(const ExperimentObs& exp_obs) const;
 
    std::shared_ptr<IModel> obs_int;                       
    std::shared_ptr<IStatCorrelationProxy> pscp;           
    std::shared_ptr<IStatParameterProxy> pspp;             
    MarginalConfigFactory marginal_config_factory_;        
    std::shared_ptr<IStatSourcesProxy> sp;                 
    std::shared_ptr<IStatDependencyPruner> dp;             
    std::shared_ptr<INuisanceReader> nuisance_reader_;     
    std::shared_ptr<IStatParamOptimizerProxy> spop;        
 
    StatisticConfig config;                                
    StatCache cache;                                       
 
    NuisanceRegistry default_nuisance_specs_;              
    NuisanceRegistry user_nuisance_specs_;                 
    NuisanceRegistry merged_nuisance_specs_;               
    std::optional<fs::path> current_user_nuisance_file_;   
 
    std::shared_ptr<LikelihoodContext> last_ctx_;          
    std::shared_ptr<BaseLikelihood> last_like_;            
    std::shared_ptr<MLFitter> last_fitter_;                
    FitResult last_fit_raw_;                               
 
    std::vector<ParamId> last_fit_param_ids_;              
    std::vector<ParamId> last_nuisance_ids_;               
    std::map<ParamId, std::size_t> last_fit_param_index_;  
 
    std::vector<double> last_scan_p_;                      
    std::vector<double> last_scan_eta_;                    
    bool has_manual_scan_point_ = false;                   
 
    std::vector<ParamId> last_detached_fit_params_;        
    std::vector<std::pair<ParameterType, std::string>> last_detached_fit_blocks_; 
 
    std::optional<std::set<std::string>> selected_experiments_; 
    std::optional<std::set<ExperimentObs>> selected_experiment_observables_; 
};
 
#endif