new tests + fixed docstrings

src/fitting.jl

@@ -19,14 +19,17 @@ function get_bkg_info(config)
 end
 
 
-function norm_linear(x::Float64,p::NamedTuple,b_name::Symbol,fit_range)
 """
+    norm_linear(x::Float64,p::NamedTuple,b_name::Symbol,fit_range)
+
 Normalised linear function defined by (1+slope*(x-center)/260)/norm.
+
 Parameters
 ----------
     - slope::Real, the slope of the background
     - x::Real,     the x value to evaluate at
 """
+function norm_linear(x::Float64,p::NamedTuple,b_name::Symbol,fit_range)
     range_l = [arr[1] for arr in fit_range]
     range_h = [arr[2] for arr in fit_range]
     center = range_l[1]
@@ -43,13 +46,16 @@ end
 
 
 
-function norm_uniform(x::Real,p::NamedTuple,b_name::Symbol,fit_range)
 """
+    norm_uniform(x::Real,p::NamedTuple,b_name::Symbol,fit_range)
+
 Normalised linear function defined by (1+slope*(x-center)/260)/norm.
+
 Parameters
 ----------
     - x::Real,     the x value to evaluate at
 """    
+function norm_uniform(x::Real,p::NamedTuple,b_name::Symbol,fit_range)
     range_l = [arr[1] for arr in fit_range]
     range_h = [arr[2] for arr in fit_range]
     center = range_l[1]
@@ -67,14 +73,19 @@ function exp_stable(x::Float64)
         return exp(x)
     end
 end
-function norm_exponential(x::Float64,p::NamedTuple,b_name::Symbol,fit_range)
+
+
 """
+    norm_exponential(x::Float64,p::NamedTuple,b_name::Symbol,fit_range)
+
 Normalised linear function defined by (1+slope*(x-center)/fit_range)/norm.
+
 Parameters
 ----------
     - slope::Real, the slope of the background
     - x::Real,     the x value to evaluate at
 """
+function norm_exponential(x::Float64,p::NamedTuple,b_name::Symbol,fit_range)
     range_l = [arr[1] for arr in fit_range]
     range_h = [arr[2] for arr in fit_range]
     center = range_l[1]
@@ -94,10 +105,13 @@ Parameters
 
 end
 
-function gaussian_plus_lowEtail(evt_energy::Float64,Qbb::Float64,bias::Float64,reso::Float64,part_k::NamedTuple)
+
 """
+    gaussian_plus_lowEtail(evt_energy::Float64,Qbb::Float64,bias::Float64,reso::Float64,part_k::NamedTuple)
+
 Signal model based on the peak shape used for the MJD analysis. The peak shape derives from considerations made in [S. I. Alvis et al., Phys. Rev. C 100, 025501 (2019)].
 """
+function gaussian_plus_lowEtail(evt_energy::Float64,Qbb::Float64,bias::Float64,reso::Float64,part_k::NamedTuple)
     γ = reso
     # following params are ALWAYS fixed
     f = part_k.frac
@@ -117,10 +131,12 @@ end
 ##############################################
 ##############################################
 ##############################################
-function get_stat_blocks(partitions,events::Array{Vector{Float64}},part_event_index,fit_ranges;config,bkg_only)
 """
+    get_stat_blocks(partitions,events::Array{Vector{Float64}},part_event_index,fit_ranges;config,bkg_only)
+
 Function to retrieve useful pieces (prior, likelihood, posterior), also in saving values
 """
+function get_stat_blocks(partitions,events::Array{Vector{Float64}},part_event_index,fit_ranges;config,bkg_only)
     settings=get_settings(config)
 
 
@@ -160,10 +176,12 @@ Function to retrieve useful pieces (prior, likelihood, posterior), also in savin
 end
 
 
-function run_fit_over_partitions(partitions,events::Array{Vector{Float64}},part_event_index::Vector{Int}, config,fit_ranges)
 """
+    run_fit_over_partitions(partitions,events::Array{Vector{Float64}},part_event_index::Vector{Int}, config,fit_ranges)
+
 Function to run the fit looping over partitions
 """
+function run_fit_over_partitions(partitions,events::Array{Vector{Float64}},part_event_index::Vector{Int}, config,fit_ranges)
     bkg_only = config["bkg_only"]
     prior,likelihood,posterior,par_names,nuisance_info = get_stat_blocks(partitions,events,part_event_index,fit_ranges,config=config,bkg_only=bkg_only)
 
@@ -186,11 +204,6 @@ end
 
 
 
-
-
-##############################################
-##############################################
-##############################################
 function get_qbb_posterior(fit_function,samples)
     qbb=[]
 
@@ -204,9 +217,6 @@ function get_qbb_posterior(fit_function,samples)
     return qbb
 end
 
-##############################################
-##############################################
-##############################################
 function get_par_posterior(samples,par;idx)
 
     pars=[]
@@ -227,9 +237,6 @@ function get_par_posterior(samples,par;idx)
 end
 
 
-##############################################
-##############################################
-##############################################
 function get_n_posterior(samples)
     ns=[]
 
@@ -244,8 +251,6 @@ function get_n_posterior(samples)
 end
 
 
-
-
 struct LogFlat <: ContinuousUnivariateDistribution
         a::Float64
         b::Float64

src/likelihood.jl

@@ -36,10 +36,13 @@ function get_signal_pdf(evt_energy::Float64,Qbb::Float64,part_k::NamedTuple)
 
 end
 
+
+""" 
+    get_energy_scale_pars(part_k::NamedTuple,p::NamedTuple,settings::Dict,idx_part_with_events)
+
+Get the resolution and bias
+"""
 function get_energy_scale_pars(part_k::NamedTuple,p::NamedTuple,settings::Dict,idx_part_with_events)
-    """ 
-    Get the resolution and bias
-    """
     if (settings[:energy_scale_fixed]==true || idx_part_with_events==0)
         reso = part_k.width
         bias = part_k.bias
@@ -60,10 +63,13 @@ return reso*1.0,bias*1.0
 
 end
 
+
+"""
+    get_mu_s_b(p::NamedTuple,part_k::NamedTuple,idx_part_with_events::Int,settings::Dict,fit_range)
+
+Get the expected number of signal and background counts in a partition
+"""
 function get_mu_s_b(p::NamedTuple,part_k::NamedTuple,idx_part_with_events::Int,settings::Dict,fit_range)
-    """
-    Get the expected number of signal and background counts in a partition
-    """
     N_A = constants.N_A
     m_76 = constants.m_76
     sig_units = constants.sig_units
@@ -101,10 +107,12 @@ function get_mu_s_b(p::NamedTuple,part_k::NamedTuple,idx_part_with_events::Int,s
 end
 
 
+"""
+    build_likelihood_zero_obs_evts(part_k::NamedTuple, p::NamedTuple,settings::Dict,fit_range)
+
+Function to calculate the partial likelihood for a partition with 0 events
+"""
 function build_likelihood_zero_obs_evts(part_k::NamedTuple, p::NamedTuple,settings::Dict,fit_range)
-    """
-    Function to calculate the partial likelihood for a partition with 0 events
-    """
 
     ll_value = 0
     model_s_k,model_b_k = get_mu_s_b(p,part_k,0,settings,fit_range)
@@ -116,12 +124,14 @@ function build_likelihood_zero_obs_evts(part_k::NamedTuple, p::NamedTuple,settin
 end
 
 
-function build_likelihood_per_partition(idx_k::Int, idx_part_with_events::Int,part_k::NamedTuple, events_k::Vector{Union{Float64}}, 
-    p::NamedTuple,settings::Dict,bkg_shape::Symbol,fit_range)
 """
+    build_likelihood_per_partition(idx_k::Int, idx_part_with_events::Int,part_k::NamedTuple, events_k::Vector{Union{Float64}},p::NamedTuple,settings::Dict,bkg_shape::Symbol,fit_range)
+
 Function which computes the partial likelihood for a single data partiton
 free parameters: signal (S), background (B), energy bias (biask) and resolution per partition (resk)
 """
+function build_likelihood_per_partition(idx_k::Int, idx_part_with_events::Int,part_k::NamedTuple, events_k::Vector{Union{Float64}}, 
+    p::NamedTuple,settings::Dict,bkg_shape::Symbol,fit_range)
     Qbb = constants.Qbb
 
     ll_value = 0
@@ -161,19 +171,23 @@ end
 
 
 
-function build_likelihood_looping_partitions(partitions::TypedTables.Table, events::Array{Vector{Float64}},
-    part_event_index::Vector{Int},settings::Dict,sqrt_prior::Bool,s_max::Union{Float64,Nothing},fit_ranges;bkg_shape::Symbol=:uniform)
 """
+    build_likelihood_looping_partitions(partitions::TypedTables.Table,events::Array{Vector{Float64}},part_event_index::Vector{Int},settings::Dict,sqrt_prior::Bool,s_max::Union{Float64,Nothing},fit_ranges;bkg_shape::Symbol=:uniform)
+
 Function which creates the likelihood function for the fit (looping over partitions)
+
 Parameters:
 -----------
     -partitions: Table - partitions input file
     -events: Array      - list of events in each partitions (with their energy)
     -nuis_prior:bool     - true if we want to include priors for nuisance parameters (bias, res, eff)
+
 Returns:
 --------
     DensityInterface.logfuncdensity - the likelihood function
 """
+function build_likelihood_looping_partitions(partitions::TypedTables.Table, events::Array{Vector{Float64}},
+    part_event_index::Vector{Int},settings::Dict,sqrt_prior::Bool,s_max::Union{Float64,Nothing},fit_ranges;bkg_shape::Symbol=:uniform)
     @debug part_event_index
     return DensityInterface.logfuncdensity( function(p::NamedTuple)
             total_ll = 0.0
@@ -202,17 +216,19 @@ end
 
 
 
-function generate_data(samples::BAT.DensitySampleVector,partitions::TypedTables.Table,part_event_index::Vector{Int},settings::Dict,fit_ranges;
-    best_fit::Bool=false,seed=nothing,bkg_only=false)
 """
+    generate_data(samples::BAT.DensitySampleVector,partitions::TypedTables.Table,part_event_index::Vector{Int},settings::Dict,fit_ranges;best_fit::Bool=false,seed=nothing,bkg_only=false)
+
 Generates data from a posterior distribution.
 This is based on the posterior predictive distributions. 
 Given a model with some parameters `theta_i`, the posterior predictive distribution,
 or the distribution of data generated according to the posterior distribution of theta
 and the likelihood is:
+
 ```math
 p(y|D) =int p(y|theta)p(theta|D)dtheta
 ```
+
 Or in terms of sampling we first draw samples of `theta` from the posterior and then generate,
 datasets based on the likelihood.
 We also give the options to fix the posterior distribution to the best fit,
@@ -223,15 +239,20 @@ Parameters
     - samples::DensitySamplesVector the samples of a past fit or a NamedTuple of best fit
     - partitions::Table of the partition info
     - part_event_index: index for the parameters for partitions with events
+
 Keyword arguments
 -----------------
     - best_fit::Bool where to fix the paramaters to the best fit
     - nuis_prior::Bool whether only statistical parameters were included in the posterior
     - bkg_only::Bool where the fit was without signal,
     - seed::Int random seed
+
 Returns
+-------
     OrderedDict of the data
 """
+function generate_data(samples::BAT.DensitySampleVector,partitions::TypedTables.Table,part_event_index::Vector{Int},settings::Dict,fit_ranges;
+    best_fit::Bool=false,seed=nothing,bkg_only=false)
     Qbb = constants.Qbb
 
     # seed the seed
@@ -294,12 +315,16 @@ end
 
 
 
-function get_signal_bkg_priors(config)
 """
+    get_signal_bkg_priors(config)
+
 Defines specific priors for signal and background contributions
+
 Parameters
+----------
     - config: the Dict of the fit config
 """
+function get_signal_bkg_priors(config)
 
     uppS = config["signal"]["upper_bound"]
     uppB = config["bkg"]["upper_bound"]
@@ -321,15 +346,19 @@ Parameters
     return distrS, distrB
 end
 
-function build_prior(partitions,part_event_index,config,settings::Dict;hierachical=false,hierachical_mode=nothing,hierachical_range=nothing,bkg_shape=:uniform,shape_pars=nothing)
+
 """
+    build_prior(partitions,part_event_index,config,settings::Dict;hierachical=false,hierachical_mode=nothing,hierachical_range=nothing,bkg_shape=:uniform,shape_pars=nothing)
+
 Builds the priors for use in the fit
-----------
+
 Parameters
+----------
     - partitions:Table of the partition info
     - config: the Dict of the fit config
     - nuis_prior; true if we want to include priors for nuisance parameters (bias, res, eff)
 """
+function build_prior(partitions,part_event_index,config,settings::Dict;hierachical=false,hierachical_mode=nothing,hierachical_range=nothing,bkg_shape=:uniform,shape_pars=nothing)
 
     # bkg indexs
     list_names = partitions.bkg_name
@@ -531,6 +560,7 @@ Parameters
     end
 end
 
+
 function print_names(priors,pretty_names)
 
     for (key, value) in pairs(priors)

src/plotting.jl

@@ -26,9 +26,6 @@ color_schemes=Dict(:blue =>[tol_colors[1],tol_colors[3],tol_colors[2]],
                    :muted=>[:olivedrab,:goldenrod,:indianred1]
 )
 
-##############################################
-##############################################
-##############################################
 function fit_model(idx_part_with_events::Int,part_k::NamedTuple, p::NamedTuple,settings::Dict,bkg_shape::Symbol,fit_range,x)
     Qbb = constants.Qbb
     N_A = constants.N_A
@@ -72,13 +69,12 @@ function fit_model(idx_part_with_events::Int,part_k::NamedTuple, p::NamedTuple,s
 end
 
 
-##############################################
-##############################################
-##############################################
-function plot_data(hist::Histogram,name,partitions,part_event_index,pars,samples,posterior,plotflag,settings::Dict,bkg_shape::Symbol,fit_ranges)
 """
+    plot_data(hist::Histogram,name,partitions,part_event_index,pars,samples,posterior,plotflag,settings::Dict,bkg_shape::Symbol,fit_ranges)
+
 Function to plot events in the Qbb analysis window and BAT fit results
 """
+function plot_data(hist::Histogram,name,partitions,part_event_index,pars,samples,posterior,plotflag,settings::Dict,bkg_shape::Symbol,fit_ranges)
 
     counts=sum(hist.weights)
     p = plot() 
@@ -145,9 +141,6 @@ end
 
 
 
-##############################################
-##############################################
-##############################################
 function plot_fit_and_data(partitions, events, part_event_index, samples, posterior, pars, output, config, fit_ranges; toy_idx=nothing)
 
     plotflag=config["plot"]
@@ -177,14 +170,13 @@ function plot_fit_and_data(partitions, events, part_event_index, samples, poster
 
 end
 
-##############################################
-##############################################
-##############################################
 
-function plot_correlation_matrix(samples,output;par_names=nothing,toy_idx=nothing)
 """
+    plot_correlation_matrix(samples,output;par_names=nothing,toy_idx=nothing)
+
 Plots the correlation matrixs
 """
+function plot_correlation_matrix(samples,output;par_names=nothing,toy_idx=nothing)
     unshaped_samples, f_flatten = bat_transform(Vector, samples)
     covariance_matrix = cov(unshaped_samples)
     var = std(unshaped_samples)
@@ -233,11 +225,13 @@ end
 ##############################################
 ##############################################
 
-function plot_marginal_distr(partitions,samples,pars,output;sqrt_prior=false,
-    priors=nothing,par_names=nothing,plot_config=nothing,s_max=nothing,hier=false,toy_idx=nothing)    
 """
+    plot_marginal_distr(partitions,samples,pars,output;sqrt_prior=false,priors=nothing,par_names=nothing,plot_config=nothing,s_max=nothing,hier=false,toy_idx=nothing) 
+
 Function to plot 1D and 2D marginalized distributions (and priors)
 """
+function plot_marginal_distr(partitions,samples,pars,output;sqrt_prior=false,
+    priors=nothing,par_names=nothing,plot_config=nothing,s_max=nothing,hier=false,toy_idx=nothing)    
 
     name = split(output, "output/")[end]
     first_sample = samples.v[1]