initial commit

Author: Cornelius-G

README.md

@@ -0,0 +1,44 @@
+# CombinationTool
+---
+## Idea:
+Tool to constrain the parameters of physics models using Bayesian inference by combining measurements of (different) observables. Especially suitable for EFT models. Bayesian inference is performed using the <a target="_blank" href="https://github.com/bat/BAT.jl">*Bayesian Analysis Toolkit - BAT.jl*</a>.
+
+The likelihood used in the CombinationTool is based on the
+<a target="_blank" href="https://link.springer.com/article/10.1140/epjc/s10052-016-4280-9">EFT*fitter* </a>
+
+*  **Assumption**: Measurements of physical quantities are approximately gaussian. This allows to combine the measurements using the following likelihood:
+
+<img src="http://latex2png.com/pngs/9fcbfbebe11bc63f39c64ac44e6bf790.png" width="1000" height="55" />
+
+* **Likelihood**:
+
+   <img src="http://latex2png.com/output//latex_9bcf1dbba64bc349f0c0d81e817d269c.png" width="600" height="65" />
+
+with the **covariance matrix**:
+
+  <img src="http://latex2png.com/pngs/5ceee38a8fa4eefd3b91ec36ce075746.png" width="550" height="65" />
+
+
+
+## Input needed:
+* **Model:** Predicting physical observables as a function of the model parameters
+* **Measurements:** Measured values of the observables, including (various types of) uncertainties
+* **Correlations:** Correlation matrices for all types of uncertainties
+
+
+## Current structure of the Julia package *EFTfitter_julia.jl*
+* /
+  * *runCombinationTool.jl:* Main routine for performing a run of the CombinationTool
+  * *inputs.jl:* Input file. Providing observables, measurements, uncertainties and correlations.
+
+* /src
+  * *CombinationTool.jl:* module definition
+
+  * *datahandling.jl:* type definitions for observables, measurements, uncertainties and correlations. Functions to add these objects to arrays.
+
+  * *combination.jl:* Calculation of the covariance matrix and the combination likelihood
+
+  * *callingBAT.jl:* passing the combination likelihood to BAT.jl
+
+* /test:
+   * Unit tests for parts of the code

examples/inputs.jl

@@ -0,0 +1,30 @@
+#============= Observables =============================================#
+function obs2(params)
+    2*params.p1.^2 + 4*params.p2
+end
+
+observables = [
+    Observable("Obs1", params -> params.p1.^4-params.p2.^2),
+    Observable("Obs2", obs2)
+]
+
+
+#============= Measurements ============================================#
+measurements = [
+    Measurement("Meas1","Obs1", 42.0, Uncertainties("stat"=>1.1,
+                                                    "syst"=>1.2)),
+
+    Measurement("Meas2","Obs2", 22.0, Uncertainties("stat"=>2.1,
+                                                    "syst"=>2.2))
+]
+
+
+#============= Correlations ============================================#
+correlations = [
+    Correlation("stat", [1.0 0.0
+                         0.0 1.0], active=false),
+
+    Correlation("syst", [1.0 0.5
+                         0.5 1.0])
+]
+#=======================================================================#

examples/runCombination.jl

@@ -0,0 +1,29 @@
+#using CombinationTool
+using BAT
+using LinearAlgebra
+using ValueShapes
+using IntervalSets
+using Plots
+
+push!(LOAD_PATH, "1_BAT2.0")
+using CombinationTool
+
+include("inputs.jl")
+
+m = createmodel(observables, measurements, correlations)
+#CombinationTool.printmodel(m)
+
+likelihood = CombinationDensity(m)
+
+prior = NamedTupleDist(
+    p1 = [-10.0..10.0],
+    p2 = [-20.0..20.0]
+)
+
+posterior = PosteriorDensity(likelihood, prior)
+
+nchains = 4
+nsamples = 10^4
+algorithm = MetropolisHastings()
+
+samples, stats = bat_sample(posterior, (nsamples, nchains), algorithm)

src/CombinationTool.jl

@@ -0,0 +1,14 @@
+__precompile__(true)
+
+module CombinationTool
+
+using BAT, StatsBase, IntervalSets
+using LinearAlgebra
+
+
+include("datahandling.jl")
+include("combination.jl")
+include("callBAT.jl")
+include("utils.jl")
+
+end # module

src/callBAT.jl

@@ -0,0 +1,29 @@
+export CombinationDensity
+
+
+struct CombinationDensity <: AbstractDensity
+    model::Model
+    invcov::Matrix{Real}
+    sqrtdetcov::Real
+end
+
+function CombinationDensity(m::Model)
+    covariance = calculate_covariancematrix(m.uncertainties, m.correlationmatrices)
+    sqrtdetcov = sqrt(abs(det(covariance)))
+
+    CombinationDensity(m, inv(covariance), sqrtdetcov)
+end
+
+
+function BAT.density_logval(
+    dens::CombinationDensity,
+    params
+)
+
+    combinemeasurements(
+        dens.model,
+        dens.invcov,
+        dens.sqrtdetcov,
+        params
+    ) 
+end 

src/combination.jl

@@ -0,0 +1,50 @@
+function calculate_covariancematrix(
+    uncertainties::AbstractArray{Real, 2},
+    correlationmatrices::Vector{Array{Real,2}}
+)
+
+    nmeas= size(uncertainties, 1)
+    nunc = size(uncertainties, 2)
+    covariance = zeros(nmeas, nmeas)
+
+    for i in 1:nunc
+        for j in 1:nmeas
+            for k in j:nmeas
+                if(j == k)
+                    covariance[j, k] += uncertainties[j, i]^2
+                else
+                    covariance[j, k] += correlationmatrices[i][j, k] * uncertainties[j, i] * uncertainties[k, i]
+                    covariance[k, j]  = covariance[j, k]
+                end
+            end
+        end
+    end
+
+    return covariance
+end
+
+
+
+function combinemeasurements(
+    m::Model,
+    invcov::Matrix{<:Real},
+    sqrtdetcov::Real,
+    parameters
+)
+
+    result::Float64=0.0
+
+    nmeas=length(m.measurement_values)
+
+    for i in 1:nmeas
+        r1 = m.measurement_values[i] - m.observable_functions[m.measurement_observables[i]].f.obj.x(parameters)[1]
+        for j in (i+1):nmeas
+            r2 = m.measurement_values[j] - m.observable_functions[m.measurement_observables[j]].f.obj.x(parameters)[1]
+            result += r1 * invcov[i,j] * r2
+        end
+        result += 0.5 * r1 * invcov[i,i] * r1
+    end
+    
+    final_result::Float64 = -result - log((2*π)^(0.5*nmeas) * sqrtdetcov) 
+    return  final_result
+end