text:
#= Custom calibration configuration file. =#

using Distributions
using StatsBase
using LinearAlgebra
using Random
using CalibrateEDMF
using CalibrateEDMF.ModelTypes
using CalibrateEDMF.DistributionUtils
using CalibrateEDMF.ReferenceModels
using CalibrateEDMF.ReferenceStats
using CalibrateEDMF.LESUtils
using CalibrateEDMF.TurbulenceConvectionUtils
using CalibrateEDMF.ModelTypes
using CalibrateEDMF.HelperFuncs
# Import EKP modules
using JLD2
using EnsembleKalmanProcesses
using EnsembleKalmanProcesses.ParameterDistributions

# Cases defined as structs for quick access to default configs
struct SCT2Train end
struct SCT2Val end

function get_config()
    config = Dict()
    # Flags for saving output data
    config["output"] = get_output_config()
    # Define regularization of inverse problem
    config["regularization"] = get_regularization_config()
    # Define reference used in the inverse problem
    config["reference"] = get_reference_config(SCT2Train())
    # Define reference used for validation
    config["validation"] = get_reference_config(SCT2Val())
    # Define the parameter priors
    config["prior"] = get_prior_config()
    # Define the kalman process
    config["process"] = get_process_config()
    # Define the SCM static configuration
    config["scm"] = get_scm_config()
    return config
end

function get_output_config()
    config = Dict()
    config["outdir_root"] = pwd()
    config["overwrite_scm_file"] = false # Flag for overwritting SCM input file
    return config
end

function get_regularization_config()
    config = Dict()
    # Regularization of observations: mean and covariance
    config["perform_PCA"] = true # Performs PCA on data
    config["variance_loss"] = 1.0e-2 # Variance truncation level in PCA
    config["normalize"] = true  # whether to normalize data by pooled variance
    config["tikhonov_mode"] = "relative" # Tikhonov regularization
    config["tikhonov_noise"] = 1.0e-6 # Tikhonov regularization
    config["dim_scaling"] = true # Dimensional scaling of the loss

    # Parameter regularization: L2 regularization with respect to prior mean.
    #  - Set to `nothing` to use prior covariance as regularizer,
    #  - Set to a float for isotropic parameter regularization.
    #  - Pass a dictionary of lists similar to config["prior_mean"] for
    #       anisotropic regularization. The dictionary must be complete.
    #       If you want to avoid regularizing a certain parameter, set the entry
    #       to [0].
    # To turn off regularization, set config["process"]["augmented"] to false.
    #
    #
    # Defaults set to batch_size/total_size to match total dataset uncertainty
    # in UKI. Feel free to set treat these as hyperparameters.
    config["l2_reg"] = Dict(
        # entrainment parameters
        "nn_ent_params" => repeat([0.0], 58),
        "turbulent_entrainment_factor" => [5.0 / 60.0],

        # diffusion parameters
        "tke_ed_coeff" => [5.0 / 60.0],
        "tke_diss_coeff" => [5.0 / 60.0],
        "static_stab_coeff" => [5.0 / 60.0],
        "tke_surf_scale" => [5.0 / 60.0],
        "Prandtl_number_0" => [5.0 / 60.0],

        # momentum exchange parameters
        "pressure_normalmode_adv_coeff" => [5.0 / 60.0],
        "pressure_normalmode_buoy_coeff1" => [5.0 / 60.0],
        "pressure_normalmode_drag_coeff" => [5.0 / 60.0],

        # surface
        "surface_area" => [5.0 / 60.0],
    )
    return config
end

function get_process_config()
    config = Dict()
    config["N_iter"] = 50
    config["N_ens"] = 50 # Must be 2p+1 when algorithm is "Unscented"
    config["algorithm"] = "Inversion" # "Sampler", "Unscented", "Inversion"
    config["noisy_obs"] = false # Choice of covariance in evaluation of y_{j+1} in EKI. True -> Γy, False -> 0
    # Artificial time stepper of the EKI.
    config["Δt"] = 60.0 / 5.0 # 1.0
    # Whether to augment the outputs with the parameters for regularization
    config["augmented"] = true
    config["failure_handler"] = "sample_succ_gauss" #"high_loss" #"sample_succ_gauss"
    return config
end

function get_reference_config(::SCT2Train)
    config = Dict()

    # AMIP data: July
    cfsite_numbers = (2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 17, 18, 19, 20, 21, 23)
    les_kwargs = (forcing_model = "HadGEM2-A", month = 7, experiment = "amip")
    ref_dirs = [get_cfsite_les_dir(cfsite_number; les_kwargs...) for cfsite_number in cfsite_numbers]
    suffixes = [get_gcm_les_uuid(cfsite_number; les_kwargs...) for cfsite_number in cfsite_numbers]
    # AMIP data: October
    cfsite_numbers = (2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 17, 18, 19, 20, 21, 23)
    les_kwargs = (forcing_model = "HadGEM2-A", month = 10, experiment = "amip")
    append!(ref_dirs, [get_cfsite_les_dir(cfsite_number; les_kwargs...) for cfsite_number in cfsite_numbers])
    append!(suffixes, [get_gcm_les_uuid(cfsite_number; les_kwargs...) for cfsite_number in cfsite_numbers])
    # AMIP data: January
    cfsite_numbers = (2, 3, 5, 7, 9, 11, 13, 21, 22, 23)
    les_kwargs = (forcing_model = "HadGEM2-A", month = 1, experiment = "amip")
    append!(ref_dirs, [get_cfsite_les_dir(cfsite_number; les_kwargs...) for cfsite_number in cfsite_numbers])
    append!(suffixes, [get_gcm_les_uuid(cfsite_number; les_kwargs...) for cfsite_number in cfsite_numbers])
    # AMIP data: April
    cfsite_numbers = (2, 3, 5, 7, 9, 11, 13, 19, 21, 23)
    les_kwargs = (forcing_model = "HadGEM2-A", month = 4, experiment = "amip")
    append!(ref_dirs, [get_cfsite_les_dir(cfsite_number; les_kwargs...) for cfsite_number in cfsite_numbers])
    append!(suffixes, [get_gcm_les_uuid(cfsite_number; les_kwargs...) for cfsite_number in cfsite_numbers])

    n_repeat = length(ref_dirs)
    config["case_name"] = repeat(["LES_driven_SCM"], n_repeat)
    # Flag to indicate whether reference data is from a perfect model (i.e. SCM instead of LES)
    config["y_reference_type"] = LES()
    config["Σ_reference_type"] = LES()
    config["y_names"] =
        repeat([["s_mean", "ql_mean", "qt_mean", "total_flux_qt", "total_flux_s", "u_mean", "v_mean"]], n_repeat)
    config["y_dir"] = ref_dirs
    config["scm_suffix"] = suffixes
    config["scm_parent_dir"] = repeat(["scm_init"], n_repeat)
    config["t_start"] = repeat([3.0 * 3600], n_repeat)
    config["t_end"] = repeat([6.0 * 3600], n_repeat)
    # Use full LES timeseries for covariance
    config["Σ_t_start"] = repeat([-5.75 * 24 * 3600], n_repeat)
    config["Σ_t_end"] = repeat([6.0 * 3600], n_repeat)
    config["batch_size"] = 5
    config["write_full_stats"] = false
    return config
end

function get_reference_config(::SCT2Val)
    config = Dict()

    # AMIP4K data: July
    cfsite_numbers = (3, 5, 8, 11, 14)
    les_kwargs = (forcing_model = "HadGEM2-A", month = 7, experiment = "amip4K")
    ref_dirs = [get_cfsite_les_dir(cfsite_number; les_kwargs...) for cfsite_number in cfsite_numbers]
    suffixes = [get_gcm_les_uuid(cfsite_number; les_kwargs...) for cfsite_number in cfsite_numbers]

    n_repeat = length(ref_dirs)
    config["case_name"] = repeat(["LES_driven_SCM"], n_repeat)
    # Flag to indicate whether reference data is from a perfect model (i.e. SCM instead of LES)
    config["y_reference_type"] = LES()
    config["Σ_reference_type"] = LES()
    config["y_names"] =
        repeat([["s_mean", "ql_mean", "qt_mean", "total_flux_qt", "total_flux_s", "u_mean", "v_mean"]], n_repeat)
    config["y_dir"] = ref_dirs
    config["scm_suffix"] = suffixes
    config["scm_parent_dir"] = repeat(["scm_init"], n_repeat)
    config["t_start"] = repeat([3.0 * 3600], n_repeat)
    config["t_end"] = repeat([6.0 * 3600], n_repeat)
    # Use full LES timeseries for covariance
    config["Σ_t_start"] = repeat([-5.75 * 24 * 3600], n_repeat)
    config["Σ_t_end"] = repeat([6.0 * 3600], n_repeat)
    config["write_full_stats"] = false
    return config
end

function get_prior_config()
    config = Dict()
    config["constraints"] = Dict(
        # entrainment parameters
        "nn_ent_params" => [repeat([no_constraint()], 58)...],
        "turbulent_entrainment_factor" => [bounded(0.0, 10.0)],

        # diffusion parameters
        "tke_ed_coeff" => [bounded(0.01, 1.0)],
        "tke_diss_coeff" => [bounded(0.01, 1.0)],
        "static_stab_coeff" => [bounded(0.01, 1.0)],
        "tke_surf_scale" => [bounded(1.0, 16.0)],
        "Prandtl_number_0" => [bounded(0.5, 1.5)],

        # momentum exchange parameters
        "pressure_normalmode_adv_coeff" => [bounded(0.0, 100.0)],
        "pressure_normalmode_buoy_coeff1" => [bounded(0.0, 10.0)],
        "pressure_normalmode_drag_coeff" => [bounded(0.0, 50.0)],

        # surface
        "surface_area" => [bounded(0.01, 0.5)],
    )

    # TC.jl prior mean
    config["prior_mean"] = Dict(
        # entrainment parameters
        "nn_ent_params" => 0.1 .* (rand(58) .- 0.5),
        "turbulent_entrainment_factor" => [0.075],

        # diffusion parameters
        "tke_ed_coeff" => [0.14],
        "tke_diss_coeff" => [0.22],
        "static_stab_coeff" => [0.4],
        "tke_surf_scale" => [3.75],
        "Prandtl_number_0" => [0.74],

        # momentum exchange parameters
        "pressure_normalmode_adv_coeff" => [0.001],
        "pressure_normalmode_buoy_coeff1" => [0.12],
        "pressure_normalmode_drag_coeff" => [10.0],

        # surface
        "surface_area" => [0.1],
    )

    config["unconstrained_σ"] = 1.0
    # Tight initial prior for Unscented
    # config["unconstrained_σ"] = 0.25
    return config
end

function get_scm_config()
    config = Dict()
    config["namelist_args"] = [
        ("time_stepping", "dt_min", 1.0),
        ("time_stepping", "dt_max", 2.0),
        ("stats_io", "frequency", 60.0),
        ("turbulence", "EDMF_PrognosticTKE", "entrainment", "NN"),
        ("turbulence", "EDMF_PrognosticTKE", "area_limiter_power", 0.0),
        ("turbulence", "EDMF_PrognosticTKE", "entr_dim_scale", "inv_z"),
    ]
    return config
end


text:
#-----------------------------------------------------------------------------# Part
struct Part{D, O<:OnlineStat} <: OnlineStat{TwoThings}
    domain::D
    stat::O
end
Part(stat::OnlineStat, domain) = Part(domain, stat)

value(o::Part) = (domain=o.domain, stat=o.stat)
Base.in(x, o::Part) = x ∈ o.domain
Base.isless(a::Part, b::Part) = isless(a.domain, b.domain)

nobs(o::Part) = nobs(o.stat)

function _fit!(o::Part, xy)
    first(xy) in o.domain || error("$(first(xy)) ∉ $(o.domain)")
    _fit!(o.stat, last(xy))
end

function _merge!(a::Part, b::Part)
    merge!(a.stat, b.stat)
    merge!(a.domain, b.domain)
end

Base.diff(a::Part, b::Part) = diff(a.domain, b.domain)

#-----------------------------------------------------------------------------#
#-----------------------------------------------------------------------------# Domains
#-----------------------------------------------------------------------------#
# Required methods:
#   - Base.merge!
#   - Base.in
#   - Base.isless
#   - Base.diff

#-----------------------------------------------------------------------------# Centroid
mutable struct Centroid{T}
    center::T
    function Centroid(x)
        x2 = x ./ one(eltype(x))
        new{typeof(x2)}(x2)
    end
end

Base.in(x::Number, c::Centroid{<:Number}) = true
Base.in(x::T, c::Centroid{T}) where {T} = true
Base.in(x, c::Centroid) = false
Base.isless(a::Centroid, b::Centroid) = isless(a.center, b.center)
Base.show(io::IO, c::Centroid) = print(io, "Centroid: $(c.center)")
Base.diff(a::Centroid, b::Centroid) = norm(a.center - b.center)

function Base.merge!(a::Part{Centroid{T}, O}, b::Part{Centroid{T}, O}) where {T <: Number, O}
    merge!(a.stat, b.stat)
    a.domain.center = smooth(a.domain.center, b.domain.center, nobs(b) / nobs(a))
    a
end
function Base.merge!(a::Part{Centroid{T}, O}, b::Part{Centroid{T}, O}) where {T, O}
    merge!(a.stat, b.stat)
    smooth!(a.domain.center, b.domain.center, nobs(b) / nobs(a))
    a
end

#-----------------------------------------------------------------------------# ClosedInterval
mutable struct ClosedInterval{T}
    first::T
    last::T
    ClosedInterval(a::T, b::T) where {T} = a ≤ b ? new{T}(a,b) : error("Arguments must be ordered: [$a, $b]")
end
Base.show(io::IO, b::ClosedInterval) = print(io, "ClosedInterval: [$(b.first), $(b.last)]")
Base.in(x, bucket::ClosedInterval) = bucket.first ≤ x ≤ bucket.last
Base.isless(a::ClosedInterval, b::ClosedInterval) = isless(a.first, b.first)
Base.diff(a::ClosedInterval, b::ClosedInterval) = a < b ? b.first - a.last : a.first - b.last
function Base.diff(a::ClosedInterval{T}, b::ClosedInterval{T}) where {T<:Dates.TimeType}
    a < b ? value(b.first) - value(a.last) : value(a.first) - value(b.last)
end

function Base.merge!(a::ClosedInterval, b::ClosedInterval)
    a.first = min(a.first, b.first)
    a.last = max(a.last, b.last)
    a
end

ext: jl


text:
function part1(dataset, wvar)
	display(cor(wvar););
	#= +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
	[1.3 Figure 2]
		In this part we are just interested in some graphics, then
		I will just consider a light version of each model
	+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ =#
	p = 2;
	h = 40;
	nc = nfields(dataset)
	tt = size(dataset[1])[1] - p
	U = fill(NaN,tt,nc);
	j = 1
	for y in dataset
		_a1,U[:,j], _a2= GShock.GSstimation(y, p, h, VarGS = 2, NF = false, nmodls = 1000);
		j+=1
	end
	umean = mean(U, dims=2)
	UB = maximum(U, dims=2) - umean
	LB = umean - minimum(U, dims=2)
	umean = umean .- mean(umean);
	gd1 = wvar[end-tt+1:end,1]
	gd1 = gd1 .- mean(gd1);
	pp = plot(1:tt, umean, ribbon = (LB, UB),c = :teal, fillalpha=0.1, w = 1.5,
	    	fg_legend = :transparent, label = "News-augmented comm. shocks",
	    	title = "Figure 2: News-augmented commodity shocks vs G20 growth" ,
	    	titlefontsize = 10,xticks = (1:12:tt, qlabel[end+1-tt:12:end]),
	    	legendfontsize=7, legend=:bottomleft);
	plot!(1:tt,gd1, c=:red, w = 1.25, alpha = 0.85, label= "G20 output");
	savefig(".//Figures//intro2.svg");

	#= +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
	[1.5 Correlations]
	+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ =#
	display(cor([umean gd1]));
	comprice = Array{Float64,2}(undef,tt+p,nc)
	[comprice[:,i] = dataset[i][:,2] for i in 1:10]
	corvec_a = vec(cor(comprice));
	corvec   = corvec_a[corvec_a .<1]
	display(sum(corvec .> 0.9))
	pp= histogram(corvec, fillalpha = 0.4, legend = :false,
		xlabel = "Correlation", lc = :black);
	savefig(".//Figures//histogram.svg");
end
#=	for i in 1:nc
	    start = 1 + nv*(i-1);
	    y = dataset[:, start:nv*i];
		model1, Uaux, _tup1 =
	    GShock.GSstimation(y, p, h, VarGS = 2, NF = false, nmodls = 1000);
	    tt = length(Uaux);
	    U[end+1-tt:end,i] = Uaux;
	end

	umean = [mean(U[i,:][.!isnan.(U[i,:])]) for i in 4:nq];
	umax  = [maximum(U[i,:][.!isnan.(U[i,:])]) for i in 4:nq];
	umin  = [minimum(U[i,:][.!isnan.(U[i,:])]) for i in 4:nq];
=#

max_stars_repo_name: ChrisVeJa/GlobalShock

max_issues_repo_name: ChrisVeJa/GlobalShock

max_forks_repo_name: ChrisVeJa/GlobalShock


text:
"""
`InitialGuessODE`: Initial guess for ordinary differential equations

### Fields

* `int`: interpolation structure
* `v`:   vector  field
* `Δt`:  time step
* `s`:   number of extrapolation stages (for initialisation)
"""
struct InitialGuessODE{TT, VT, IT <: Extrapolation}
    int::IT
    v::VT
    Δt::TT
    s::Int

    function InitialGuessODE(int::IT, v::VT, Δt::TT) where {TT <: Real, VT <: Function, IT <: Extrapolation}
        new{TT,VT,IT}(int, v, Δt, get_config(:ig_extrapolation_stages))
    end
end

function InitialGuessODE(interp::Type{<:Extrapolation}, v::Function, Δt)
    int = interp(zero(Δt), Δt)
    InitialGuessODE(int, v, Δt)
end

function InitialGuessODE(interp::Type{<:Extrapolation}, equation::Union{ODE,DAE,IODE,LODE}, Δt)
    InitialGuessODE(interp, _get_v̄(equation), Δt)
end

InitialGuess(interp, equation::Union{ODE,DAE}, Δt) = InitialGuessODE(interp, equation, Δt)
InitialGuess(equation::Union{ODE,DAE}, Δt) = InitialGuessODE(get_config(:ig_extrapolation), equation, Δt)

InitialGuess(interp, equation::ODE, Δt) = InitialGuessODE(interp, equation, Δt)
InitialGuess(interp, equation::DAE, Δt) = InitialGuessODE(interp, equation, Δt)

Base.:(==)(ig1::InitialGuessODE{TT1}, ig2::InitialGuessODE{TT2}) where {TT1,TT2}= (
                                TT1 == TT2
                             && ig1.int == ig2.int
                             && ig1.v   == ig2.v
                             && ig1.Δt  == ig2.Δt
                             && ig1.s   == ig2.s)

"Initialise initial guess, i.e., given t₀, t₁, q₁ compute q₀, v₀, v₁."
function initialize!(ig::InitialGuessODE{TT},
                t₁::TT,
                q₁::SolutionVector{DT},
                v₁::SolutionVector{DT},
                t₀::TT,
                q₀::SolutionVector{DT},
                v₀::SolutionVector{DT}) where {DT,TT}
    _midpoint_extrapolation_ode!(ig.v, t₁, t₀, q₁, q₀, ig.s)
    ig.v(t₀, q₀, v₀)
    ig.v(t₁, q₁, v₁)
end

"compute vector field of new solution"
function update_vector_fields!(ig::InitialGuessODE{TT}, t₁::TT,
                               q₁::SolutionVector{DT},
                               v₁::Vector{DT}) where {DT,TT}
    ig.v(t₁, q₁, v₁)
end

function GeometricBase.evaluate!(ig::InitialGuessODE{TT},
                q₀::SolutionVector{DT},
                v₀::SolutionVector{DT},
                q₁::SolutionVector{DT},
                v₁::SolutionVector{DT},
                guess::SolutionVector{DT},
                c::TT=one(TT)) where {DT,TT}

    if q₀ == q₁
        @warn "q₀ and q₁ in initial guess are identical! Setting q=q₁."
        guess .= q₁
    else
        evaluate!(ig.int, q₀, q₁, v₀, v₁, (one(TT)+c)*ig.Δt, guess)
    end
end

function GeometricBase.evaluate!(ig::InitialGuessODE{TT},
                q₀::SolutionVector{DT},
                v₀::SolutionVector{DT},
                q₁::SolutionVector{DT},
                v₁::SolutionVector{DT},
                guess_q::SolutionVector{DT},
                guess_v::SolutionVector{DT},
                c::TT=one(TT)) where {DT,TT}
                @assert length(guess_q) == length(guess_v)

    if q₀ == q₁
        @warn "q₀ and q₁ in initial guess are identical! Setting q=q₁ and v=0."
        guess_q .= q₁
        guess_v .= 0
    else
        evaluate!(ig.int, q₀, q₁, v₀, v₁, (one(TT)+c)*ig.Δt, guess_q, guess_v)
    end
end


text:
using CSV
using DataFrames
using Plots
using JLD2
using DiscretisedFluidQueues
import Distributions
include("/Users/anguslewis/Documents/SFFMProject/DFQExamples.jl/plot_utils.jl")
include("/Users/anguslewis/Documents/SFFMProject/DFQExamples.jl/error_metrics.jl")

@halfwidth_plot_defaults()
ic = "point_mass"
plotlyjs()
function p()
ks_data = CSV.read((@__DIR__)*"/../errors/ks_point_mass.csv",DataFrame)
ks_data.DG_lin_lim = CSV.read("/Users/anguslewis/Documents/SFFMProject/DFQExamples.jl/hitting_times_linearDGwLimiter/reflecting_model/data/"*ic*"/data2/errors/ks_error.csv",DataFrame).DG_limit

ks_data_lwr = CSV.read((@__DIR__)*"/../errors/ks_lwr_point_mass.csv",DataFrame)
ks_data_lwr.DG_lin_lim = CSV.read("/Users/anguslewis/Documents/SFFMProject/DFQExamples.jl/hitting_times_linearDGwLimiter/reflecting_model/data/"*ic*"/data2/errors/ks_error_lwr.csv",DataFrame).DG_limit

ks_data_upr = CSV.read((@__DIR__)*"/../errors/ks_upr_point_mass.csv",DataFrame)
ks_data_upr.DG_lin_lim = CSV.read("/Users/anguslewis/Documents/SFFMProject/DFQExamples.jl/hitting_times_linearDGwLimiter/reflecting_model/data/"*ic*"/data2/errors/ks_error_upr.csv",DataFrame).DG_limit

lwr = [log10(ks_data[i,j]) .- log10(ks_data_lwr[i,j]) for i in 1:size(ks_data,1), j in 1:size(ks_data,2)]
upr = [log10(ks_data_upr[i,j]) .- log10(ks_data[i,j]) for i in 1:size(ks_data,1), j in 1:size(ks_data,2)]
ks_data = log10.(ks_data)
plot()
linestyles_vec = [:solid,:dash,:dashdot, :dot,:dash]
markerstyles_vec = [:false,:false,:false,:false,:false]
_names = ["DG"; "DG lim"; "Unif"; "QBD-RAP"; "DG lin lim"]
colours = [1;3;2;4;5]
for (c,col) in enumerate(names(ks_data))
    os = collect(1:2:21)
    (c==5)&&(os.+=1)
    plot!(log10.(os),ks_data[:,col], ribbon=(lwr[:,c],upr[:,c]),
        label=_names[c],linestyle=linestyles_vec[c],
        marker=markerstyles_vec[c],
        linewidth=2,color=colours[c])
end
plot!(xlabel="Dimension")
plot!(ylabel="Error")
plot!(title="KS error")
# plot!(legend=(0.225,0.55))
error_ticks!(plot!())
@add_lines!(ks_data,("Unif","QBDRAP"),@__DIR__)

savefig((@__DIR__)*"/ks_error_formatted.pdf")
end
p()

function p()
@halfwidth_plot_defaults()

L1_cell_probs_data = log10.(CSV.read((@__DIR__)*"/../errors/cell_probs_point_mass.csv",DataFrame))
L1_cell_probs_data.DG_lin_lim = log10.(CSV.read("/Users/anguslewis/Documents/SFFMProject/DFQExamples.jl/hitting_times_linearDGwLimiter/reflecting_model/data/"*ic*"/data2/errors/l1_cell_probs_error.csv",DataFrame).DG_limit)

L1_cell_probs_data_lwr = log10.(CSV.read((@__DIR__)*"/../errors/cell_probs_lwr_point_mass.csv",DataFrame))
L1_cell_probs_data_lwr.DG_lin_lim = log10.(CSV.read("/Users/anguslewis/Documents/SFFMProject/DFQExamples.jl/hitting_times_linearDGwLimiter/reflecting_model/data/"*ic*"/data2/errors/l1_cell_probs_error_lwr.csv",DataFrame).DG_limit)

L1_cell_probs_data_upr = log10.(CSV.read((@__DIR__)*"/../errors/cell_probs_upr_point_mass.csv",DataFrame))
L1_cell_probs_data_upr.DG_lin_lim = log10.(CSV.read("/Users/anguslewis/Documents/SFFMProject/DFQExamples.jl/hitting_times_linearDGwLimiter/reflecting_model/data/"*ic*"/data2/errors/l1_cell_probs_error_upr.csv",DataFrame).DG_limit)

lwr = [(L1_cell_probs_data[i,j]) .- (L1_cell_probs_data_lwr[i,j]) for i in 1:size(L1_cell_probs_data,1), j in 1:size(L1_cell_probs_data,2)]
upr = [(L1_cell_probs_data_upr[i,j]) .- (L1_cell_probs_data[i,j]) for i in 1:size(L1_cell_probs_data,1), j in 1:size(L1_cell_probs_data,2)]
plot()
linestyles_vec = [:solid,:dash,:dashdot,:dot,:dash]
markerstyles_vec = [:false,:false,:false,:false,:false]
_names = ["DG"; "DG lim"; "Unif"; "QBD-RAP"; "DG lin lim"]
colours = [1;3;2;4;5]
for (c,col) in enumerate(names(L1_cell_probs_data))
    os = collect(1:2:21)
    (c==5)&&(os.+=1)
    plot!(log10.(os),L1_cell_probs_data[:,col], ribbon=(lwr[:,c],upr[:,c]),
        label=_names[c],linestyle=linestyles_vec[c],
        marker=markerstyles_vec[c],
        linewidth=2,color=colours[c])
end
plot!(xlabel="Dimension")
plot!(ylabel="Error"); error_ticks!(plot!())
plot!(title="Cell-wise error")
# plot!(legend=(0.225,0.55))
plot!()
@add_lines!(L1_cell_probs_data,("Unif","QBDRAP"),@__DIR__)

savefig((@__DIR__)*"/L1_cell_probs_error_formatted.pdf")
end
p()

function p()
l1_cdf_data = log10.(CSV.read((@__DIR__)*"/../errors/l1_point_mass.csv",DataFrame))
l1_cdf_data.DG_lin_lim = log10.(CSV.read("/Users/anguslewis/Documents/SFFMProject/DFQExamples.jl/hitting_times_linearDGwLimiter/reflecting_model/data/"*ic*"/data2/errors/l1_error.csv",DataFrame).DG_limit)

l1_cdf_data_lwr = log10.(CSV.read((@__DIR__)*"/../errors/l1_lwr_point_mass.csv",DataFrame))
l1_cdf_data_lwr.DG_lin_lim = log10.(CSV.read("/Users/anguslewis/Documents/SFFMProject/DFQExamples.jl/hitting_times_linearDGwLimiter/reflecting_model/data/"*ic*"/data2/errors/l1_error_lwr.csv",DataFrame).DG_limit)

l1_cdf_data_upr = log10.(CSV.read((@__DIR__)*"/../errors/l1_upr_point_mass.csv",DataFrame))
l1_cdf_data_upr.DG_lin_lim = log10.(CSV.read("/Users/anguslewis/Documents/SFFMProject/DFQExamples.jl/hitting_times_linearDGwLimiter/reflecting_model/data/"*ic*"/data2/errors/l1_error_upr.csv",DataFrame).DG_limit)

lwr = [(l1_cdf_data[i,j]) .- (l1_cdf_data_lwr[i,j]) for i in 1:size(l1_cdf_data,1), j in 1:size(l1_cdf_data,2)]
upr = [(l1_cdf_data_upr[i,j]) .- (l1_cdf_data[i,j]) for i in 1:size(l1_cdf_data,1), j in 1:size(l1_cdf_data,2)]
@halfwidth_plot_defaults()
plot()
markerstyles_vec = [:false,:false,:false,:false,:false]
_names = ["DG"; "DG lim"; "Unif"; "QBD-RAP"; "DG lin lim"]
colours = [1;3;2;4;5]
linestyles_vec = [:solid,:dash,:dashdot,:dot,:dash]
for (c,col) in enumerate(names(l1_cdf_data))
    os = collect(1:2:21)
    (c==5)&&(os.+=1)
    plot!(log10.(os),l1_cdf_data[:,col],ribbon=(lwr[:,c],upr[:,c]),
        label=_names[c],linestyle=linestyles_vec[c],
        marker=markerstyles_vec[c],
        linewidth=2,color=colours[c])
end
plot!(xlabel="Dimension")
plot!(ylabel="Error"); error_ticks!(plot!())
plot!(title="L¹ error - CDF")
# plot!(legend=(0.225,0.55))
plot!()
@add_lines!(l1_cdf_data,("Unif","QBDRAP"),@__DIR__)

savefig((@__DIR__)*"/l1_cdf_error_formatted.pdf")
end
p()

function p()
    @fullwidth_plot_defaults

    dg_lim_lin = CSV.read("/Users/anguslewis/Documents/SFFMProject/DFQExamples.jl/hitting_times_linearDGwLimiter/reflecting_model/data/point_mass/data2/cdfs/order_21_model_limiter.csv", DataFrame)
    plot(dg_lim_lin.x,dg_lim_lin.phase_1,label="DG lim lin",color=4,linewidth=2)

    qbdrap = CSV.read("/Users/anguslewis/Documents/SFFMProject/DFQExamples.jl/hitting_times_model/reflecting_model/transient_distribution/point_mass/datat2/approx_cdf/order_21model_qbdrap4.csv",DataFrame)
    plot!(qbdrap.x,qbdrap.phase_1,label="QBD-RAP",color=5,linewidth=2)

    sims = CSV.read("/Users/anguslewis/Documents/SFFMProject/DFQExamples.jl/hitting_times_model/reflecting_model/transient_distribution/point_mass/datat2/sims_cdf_evaluated.csv",DataFrame)
    plot!(sims.x,sims.phase_1,label="Simulation",color=:grey,linewidth=2)

    plot!(title="CDF",xlabel="x",ylabel="Probability",xlims = (1.5,2.6),ylims = (0.4,0.65))

    savefig((@__DIR__)*"/order_21_CDFs_formatted.pdf")
end
p()


text:
Network = ResNet

netparams = ResNetHP(
  num_filters=128,
  num_blocks=10,
  conv_kernel_size=(3, 3),
  num_policy_head_filters=32,
  num_value_head_filters=32,
  batch_norm_momentum=0.1)

self_play = SelfPlayParams(
  num_games=4000,
  num_workers=128,
  use_gpu=true,
  reset_mcts_every=4,
  mcts=MctsParams(
    num_iters_per_turn=600,
    cpuct=3.0,
    prior_temperature=1.0,
    temperature=PLSchedule([0, 20], [1.0, 0.3]),
    dirichlet_noise_ϵ=0.2,
    dirichlet_noise_α=1.0))

arena = ArenaParams(
  num_games=128,
  num_workers=128,
  use_gpu=true,
  reset_mcts_every=2,
  flip_probability=0.5,
  update_threshold=0.05,
  mcts=MctsParams(
    self_play.mcts,
    temperature=ConstSchedule(0.2),
    dirichlet_noise_ϵ=0.05))

learning = LearningParams(
  use_gpu=true,
  use_position_averaging=true,
  samples_weighing_policy=LOG_WEIGHT,
  batch_size=512,
  loss_computation_batch_size=1024,
  optimiser=Adam(lr=5e-4),
  l2_regularization=1e-4,
  nonvalidity_penalty=1.,
  min_checkpoints_per_epoch=1,
  max_batches_per_checkpoint=4000,
  num_checkpoints=1)

params = Params(
  arena=arena,
  self_play=self_play,
  learning=learning,
  num_iters=50,
  ternary_rewards=true,
  use_symmetries=true,
  memory_analysis=MemAnalysisParams(
    num_game_stages=4),
  mem_buffer_size=PLSchedule(
  [      0,        40],
  [400_000, 2_000_000]))

mcts_baseline =
  Benchmark.MctsRollouts(
    MctsParams(
      arena.mcts,
      num_iters_per_turn=1000,
      cpuct=1.))

minmax_baseline = Benchmark.MinMaxTS(depth=5, amplify_rewards=true, τ=0.2)

players = [
  Benchmark.Full(arena.mcts),
  Benchmark.Full(arena.mcts),
  Benchmark.NetworkOnly(τ=0.5)]

baselines = [
  mcts_baseline,
  minmax_baseline,
  mcts_baseline]

make_duel(player, baseline) =
  Benchmark.Duel(
    player,
    baseline,
    num_games=128,
    num_workers=128,
    use_gpu=true,
    flip_probability=0.5,
    color_policy=CONTENDER_WHITE)

benchmark = [make_duel(p, b) for (p, b) in zip(players, baselines)]


text:
module QueueUnion

import GraphIdx.Tree: ChildrenIndex, reset!, dfs_walk_rev

"""
    Range(start, stop)

A range in a `Vector`.
In our algorithm is is used for the nodes of [`Queues`](@ref) of [`Event`](@ref)s.

!!! warning

    unlike in Python, `start` and `stop` are inclusive!
"""
struct Range
    start::Int
    stop::Int
end

Base.length(q::Range) = q.stop - q.start + 1

Base.show(io::IO, q::Range) = print(io, q.start, ":", q.stop)

Base.convert(::Type{Range}, r::UnitRange{Int}) =
    Range(r.start, r.stop)

function merge(elements::Vector{E}, parent::Range, child::Range)::Range where {E}
    if parent.start <= parent.stop
        gap = child.start - parent.stop - 1
        old_stop = parent.stop
        res = Range(parent.start, child.stop - gap)
        if gap > 0
            for i in old_stop+1:res.stop
                elements[i] = elements[i+gap]
            end
        end
        sort!(elements, Int(res.start), Int(res.stop),
              PartialQuickSort(Int(res.start):Int(res.stop)),
              Base.Sort.Forward)
        return res
    end
    return child
end

function merge(elements::Vector{E}, pq::Vector{Range}, p::I, i::I) where {E,I}
    if pq[p].start <= pq[p].stop
        gap = pq[i].start - pq[p].stop - 1
        old_stop = pq[p].stop
        res = Range(pq[p].start, pq[i].stop - gap)
        if gap > 0
            for i in old_stop+1:res.stop
                elements[i] = elements[i+gap]
            end
        end
        sort!(elements, Int(res.start), Int(res.stop),
              PartialQuickSort(Int(res.start):Int(res.stop)),
              Base.Sort.Forward)
        pq[p] = res
    else
        pq[p] = pq[i]
    end
end

"""
Set of double ended queues (DeQue) in a tree.
Initially all leaf nodes have its own queue.
Children's queues will be merged into the parent node.
"""
struct Queues{E}
    events::Vector{E}
    pq::Vector{Range}
end

Queues{E}(n::Integer) where {E} =
    Queues{E}(Vector{E}(undef, 2n),
              Vector{Range}(undef,  n))

Base.getindex(q::Queues, i) = q.pq[i]

function merge!(qs::Queues{E}, i::I, j::I) where {I,E}
    qs.pq[i] = merge(qs.events, qs.pq[i], qs.pq[j])
end

function reset!(
    pq::Vector{Range},
    childs::ChildrenIndex,
    stack::Vector{I} = I[],
    proc_order::Vector{I} = I[],
) where {I}
    @assert length(pq) == length(childs)
    resize!(proc_order, length(childs))
    empty!(proc_order)
    @assert isempty(proc_order)
    empty!(stack)

    local t::Ref{Int} = Ref{Int}(0)
    dfs_walk_rev(childs, stack) do v::Int
        t[] += 1
        if v >= 0
            push!(proc_order, v)
            pq[v] = Range(t[], t[]-1)
        end
    end

    return proc_order
end

end

max_stars_repo_name: EQt/treelas

max_issues_repo_name: EQt/treelas

max_forks_repo_name: EQt/treelas


text:
using FactCheck
using LLVM.Ast

using DataStructures

facts("Data Layout") do
    context("litle-endian") do
        str = "e"
        @fact string(Ast.DataLayout(endianness=Ast.LittleEndian())) => str
        @fact string(LLVM.parse_datalayout(str)) => str
    end

    context("big-endian") do
        str = "E"
        @fact string(Ast.DataLayout(endianness=Ast.BigEndian())) => str
        @fact string(LLVM.parse_datalayout(str)) => str
    end

    context("native") do
        str = "n8:32"
        @fact string(Ast.DataLayout(nativesizes=Set(8,32))) => str
        @fact string(LLVM.parse_datalayout(str)) => str
    end

    context("pointerlayouts") do
        str = "p:8:64"
        as = Ast.AddrSpace(0)
        ai = Ast.AlignmentInfo(64, nothing)
        @fact string(Ast.DataLayout(pointerlayouts=[as=>(8, ai)])) => str
        @fact string(LLVM.parse_datalayout(str)) => str
    end

    context("pointerlayouts") do
        str = "p1:8:32:64"
        as = Ast.AddrSpace(1)
        ai = Ast.AlignmentInfo(32, 64)
        @fact string(Ast.DataLayout(pointerlayouts=[as=>(8, ai)])) => str
        @fact string(LLVM.parse_datalayout(str)) => str
    end

    context("big") do
        str = "e-m:e-S128-p:64:64:64-i1:8:8-i8:8:8-i16:16:16-i32:32:32-i64:64:64-v64:64:64-v128:128:128-f32:32:32-f64:64:64-f80:128:128-a0:0:64-n8:16:32:64"
        dl = Ast.DataLayout(
                endianness=Ast.LittleEndian(),
                mangling=Ast.ELFMangling(),
                stackalignment=128,
                pointerlayouts=[Ast.AddrSpace(0) => (64, Ast.AlignmentInfo(64, 64))],
                typelayouts=OrderedDict(
                             [((Ast.IntegerAlign(),  1), Ast.AlignmentInfo(8, 8)),
                              ((Ast.IntegerAlign(),  8), Ast.AlignmentInfo(8, 8)),
                              ((Ast.IntegerAlign(), 16), Ast.AlignmentInfo(16, 16)),
                              ((Ast.IntegerAlign(), 32), Ast.AlignmentInfo(32, 32)),
                              ((Ast.IntegerAlign(), 64), Ast.AlignmentInfo(64, 64)),
                              ((Ast.VectorAlign(),  64), Ast.AlignmentInfo(64, 64)),
                              ((Ast.VectorAlign(), 128), Ast.AlignmentInfo(128, 128)),
                              ((Ast.FloatAlign(),   32), Ast.AlignmentInfo(32, 32)),
                              ((Ast.FloatAlign(),   64), Ast.AlignmentInfo(64, 64)),
                              ((Ast.FloatAlign(),   80), Ast.AlignmentInfo(128, 128)),
                              ((Ast.AggregateAlign(),0), Ast.AlignmentInfo(0, 64))]),
                nativesizes=Set(8,16,32,64))
        @fact string(dl) => str
        @fact string(LLVM.parse_datalayout(str)) => str
    end
end


text:
using BenchmarkTools
const SUITE = BenchmarkGroup()

using PNML

bench = SUITE["label"] = BenchmarkGroup()

# FROM CATLAB
function benchmark_pullback(suite::BenchmarkGroup, name, arg)
  for alg in (NestedLoopJoin(), SortMergeJoin(), HashJoin())
    suite["$name:$(nameof(typeof(alg)))"] =
      @benchmarkable limit($arg, alg=$alg)
  end
end

sizes = (100, 150)
f, g = (FinFunction(ones(Int, size), 1) for size in sizes)
benchmark_pullback(bench, "pullback_terminal", Cospan(f, g))


text:
using PyAMG
using Base.Test

L1d(N) = spdiagm((-ones(N-1), 2*ones(N), -ones(N-1)), (-1,0,1), N, N) * N^2
function L2d(N)
  L1 = L1d(N)
  A = kron(speye(N), L1) + kron(L1, speye(N))
  b = ones(size(A,1))
  return A, b
end
function L3d(N)
  L1 = L1d(N)
  I = speye(N)
  A = kron(L1, I, I) + kron(I, L1, I) + kron(I, I, L1)
  b = ones(size(A,1))
  return A, b
end

println("=================================================")
println("Test 1: Julia `\\` vs PyAMG Blackbox `solve` (2D) Laplacian")
A, b = L2d(100)
print(" \\ : ")
@time x1 = A \ b
print("PyAMG-Warmup: ")
@time x2 = PyAMG.solve(A, b; tol=1e-12, verb=false);
print("PyAMG: ")
@time x3 = PyAMG.solve(A, b; tol=1e-12, verb=false)
println("|x_\\ - x_amg|_∞ = ", norm(x3-x1, Inf))
println("|A x_amg - b|_∞ = ", norm(A*x2 - b, Inf))
@assert norm(x3-x1, Inf) < 1e-10

println("=================================================")
println("Test 2: Julia `\\` vs PyAMG Blackbox `solve` (3D)")
A, b = L3d(30)
print(" \\ : ")
@time x1 = A \ b
print("PyAMG-Warmup: ")
@time x2 = PyAMG.solve(A, b; tol=1e-12, verb=false)
print("PyAMG: ")
@time x3 = PyAMG.solve(A, b; tol=1e-12, verb=false)
println("|x_\\ - x_amg|_∞ = ", norm(x3-x1, Inf))
println("|A x_amg - b|_∞ = ", norm(A*x2 - b, Inf))
@assert norm(x3-x1, Inf) < 1e-9

println("=================================================")
println("Test 3: RugeStubenSolver (3D)")
println("        50 x 50 x 50 grid = 125k dofs ")
A, b = L3d(50)
print("Create solver: ")
@time amg = RugeStubenSolver(A)
print("First solve: ")
@time x1 = solve(amg, b)
print("Second solve: (tol 1e-9)")
@time x2 = solve(amg, b; tol=1e-9)
print("Third solve: (tol 1e-6)")
@time x3 = solve(amg, b; tol=1e-6)
println("|A x_amg - b|_∞ = ", norm(A * x3 - b, Inf))
@assert norm(A * x2 - b, Inf) < 1e-8

println("=================================================")
println("Test 4: \\, ldiv, and kwargs")
A, b = L2d(50)
amg = RugeStubenSolver(A)
x_solve = solve(amg, b, tol=1e-6, cycle="V", accel="cg")
set_kwargs!(amg, tol=1e-6, cycle="V", accel="cg")
x_bs = amg \ b
x_ldiv = A_ldiv_B!(zeros(x_solve), amg, b)
@assert x_solve == x_bs
@assert x_solve == x_ldiv

println("=================================================")
println("Test 5: AMG as a preconditioner")
println("        100 x 100 Dirichlet problem, TOL = 1e-4")
println("        PyAMG vs CG vs PCG  (using IterativeSolvers)")
A, b = L2d(100)
amg = RugeStubenSolver(A)
M = aspreconditioner(amg)
using IterativeSolvers
TOL = 1e-4
println("Plain CG:")
@time x_cg, _ = IterativeSolvers.cg(A, b, 1; tol=TOL)
@time x_cg, _ = IterativeSolvers.cg(A, b, 1; tol=TOL)
println("PyAMG-preconditionerd CG:  (see `aspreconditioner`)")
@time x_pcg, _ = IterativeSolvers.cg(A, b, M; tol=TOL)
@time x_pcg, _ = IterativeSolvers.cg(A, b, M; tol=TOL)
println("PyAMG solver")
@time x_pyamg = PyAMG.solve(amg, b; tol=TOL*1e-2, accel="cg")
@time x_pyamg = PyAMG.solve(amg, b; tol=TOL*1e-2, accel="cg")
x = A \ b
println("|x_cg-x| = ", norm(x_cg - x), " \n",
         "|x_pcg-x| = ", norm(x_pcg - x), "\n",
         "|x_pyamg-x| = ", norm(x_pyamg - x) )
@assert norm(x_cg - x) < 1e-5
@assert norm(x_pcg - x) < 1e-5
@assert norm(x_pyamg - x) < 1e-5


text:
using DataFrames
using ProgressMeter

abstract type SimulationResult end

function _logmeanexp(x::AbstractArray)
    x_max = maximum(x)
    log(mean(xi -> exp(xi - x_max), x)) + x_max
end

"""
    logmeanexp(x[; dims=nothing])

Compute log(mean(exp(x))) in a numerically stable way.
"""
logmeanexp(x::AbstractArray; dims=nothing) = if dims === nothing _logmeanexp(x) else mapslices(_logmeanexp, x, dims=dims) end

include("ThermodynamicIntegration.jl")
include("AIS.jl")
include("DirectMC.jl")
include("SMC.jl")

"""
    mutual_information(system, algorithm; num_samples=1, progress=true)

Perform a simulation to compute the mutual information between input
and output trajectories of `system`.

# Arguments

The required marginalization integrals to obtain
the marginal probability ``\\mathcal{P}[\\bm{x}]`` are performed using the
specified `algorithm`.

Overall, `num_samples` Monte Carlo samples are performed. For each
individual sample, one or mupltiple marginalization operations need to be performed.

If `progress == true`, a progress bar will be shown during the computation.

# Returns

Returns a `DataFrame` containing the results of the simulation. This resulting
`DataFrame` has 3 columns. Assuming, the returned value has been named `result`
the columns can be accessed by:

- `result.MutualInformation`: A vector of vectors that contains the results of the simulation. Each element of the outer vector is the result of a single Monte Carlo sample. Each element is a vector containing the trajectory mutual information estimates for each time specified in `system.dtimes`.
- `result.TimeMarginal`: A vector containing, for each sample, the CPU time in seconds used for the computation of the marginal entropy.
- `result.TimeConditional`: A vector containing, for each sample, the CPU time in seconds used for the computation of the conditional entropy.
"""
function mutual_information(system, algorithm; num_samples::Integer=1, progress=true, compile_args = (;))
    # initialize the ensembles
    compiled_system = compile(system; compile_args...)

    # this is the outer Direct Monte-Carlo loop
    # result = Base.invokelatest(_mi_inner, compiled_system, algorithm, num_samples)
    result = _mi_inner(compiled_system, algorithm, num_samples, progress)

    result
end

function _mi_inner(compiled_system, algorithm, num_samples, show_progress)
    stats = DataFrame(
        TimeConditional=zeros(Float64, num_samples),
        TimeMarginal=zeros(Float64, num_samples),
    )

    p = Progress(num_samples; showspeed=true, enabled=show_progress)
    result = progress_map(1:num_samples, progress=p) do i
        # draw an independent sample
        sample = generate_configuration(compiled_system.system)

        # compute P(x|s) and P(x)
        cond_result = @timed conditional_density(compiled_system, algorithm, sample)
        marg_result = @timed marginal_density(compiled_system, algorithm, sample)

        # record the simulation time
        stats.TimeConditional[i] = cond_result.time
        stats.TimeMarginal[i] = marg_result.time

        # summary of the configuration
        dtimes = vec(compiled_system.system.dtimes)
        summary = vcat(sample.s_traj(dtimes), sample.x_traj(dtimes))

        # compute a sample for the mutual information using
        # ln [P(x,s)/(P(x)P(s))] = ln [P(x|s)/P(x)] = ln P(x|s) - ln P(x)
        cond_result.value - marg_result.value, summary
    end

    stats[!, :MutualInformation] = getindex.(result, 1)
    stats[!, :Trajectory] = getindex.(result, 2)

    stats
end


text:
using SemiDiscretizationMethod
using Test

function createHayesProblem(a,b)
    AMx =  ProportionalMX(a*ones(1,1));
    τ1=1.
    BMx1 = DelayMX(τ1,b*ones(1,1));
    cVec = Additive(ones(1))
    LDDEProblem(AMx,[BMx1],cVec)
end

function createMathieuProblem(δ,ε,b0,a1;T=2π)
    AMx =  ProportionalMX(t->@SMatrix [0. 1.; -δ-ε*cos(2π/T*t) -a1]);
    τ1=2π # if function is needed, the use τ1 = t->foo(t)
    BMx1 = DelayMX(τ1,t->@SMatrix [0. 0.; b0 0.]);
    cVec = Additive(t->@SVector [0.,sin(4π/T*t)])
    LDDEProblem(AMx,[BMx1],cVec)
end

function tests()
    @testset "Testing the SemiDiscretizationMethod package with the examples" begin
		#Hayes
        @test begin
            hayes_lddep=createHayesProblem(-1.,-1.); # LDDE problem for Hayes equation
			method=SemiDiscretization(1,0.1) # 3rd order semi discretization with Δt=0.1
			τmax=1. # the largest τ of the system
			mapping=DiscreteMapping(hayes_lddep,method,τmax,n_steps=1,calculate_additive=true); #The discrete mapping of the system
			spectralRadiusOfMapping(mapping)
			fixPointOfMapping(mapping)
            true
        end
		#Delay Mathieu
        @test begin
            τmax=2π # the largest τ of the system
			T=2π #Principle period of the system (sin(t)=sin(t+P))
			mathieu_lddep=createMathieuProblem(3.,2.,-0.15,0.1,T=T); # LDDE problem for Hayes equation
			method=SemiDiscretization(1,0.01) # 3rd order semi discretization with Δt=0.1
			# if P = τmax, then n_steps is automatically calculated
			mapping=DiscreteMapping(mathieu_lddep,method,τmax,
		    n_steps=Int((T+100eps(T))÷method.Δt),calculate_additive=true); #The discrete mapping of the system

			spectralRadiusOfMapping(mapping); # spectral radius ρ of the mapping matrix (ρ>1 unstable, ρ<1 stable)
			# spectralRadiusOfMapping(mapping) = 0.5131596340374617
			fp=fixPointOfMapping(mapping); # stationary solution of the hayes equation (equilibrium position)
			true
        end
    end
end
tests()


text:
abstract type AbstractDomain{N} <: AbstractVector{Union{Float64, Vector{Float64}}} end

struct Domain{N} <: AbstractDomain{N}
    coords::AbstractMatrix{Float64}

    Domain(coords::Matrix{T}) where {T <: Real} =
        new{size(coords)[2]}(coords)

    function Domain(coords::Vararg{AbstractVector{T},N}) where {T <: Real, N}
        (N > 1)  &&  (@assert all(length(coords[1]) .== [length.(coords)...]))
        mat = deepcopy(coords[1])
        if N > 1
            for i in 2:N
                mat = hcat(mat, coords[i])
            end
        else
            mat = reshape(mat, length(mat), 1)
        end
        return new{N}(mat)
    end

    function Domain(length::Integer)
        @assert length >= 1
        return new{1}(reshape(collect(1.:length), length, 1))
    end
end

ndims(d::Domain{N}) where N = N
length(d::Domain) = size(d.coords)[1]
size(d::Domain{N}) where N = (length(d),)

getindex(d::Domain{1}, ::Colon) = d.coords[:, 1]
getindex(d::Domain{1}, index::Integer) = d.coords[index, 1]
getindex(d::Domain, dim::Integer) = d.coords[:, dim]

function iterate(d::Domain{1}, ii=1)
    (ii > size(d.coords)[1])  &&  (return nothing)
    return (d.coords[ii, 1], ii+1)
end
function iterate(d::Domain, ii=1)
    (ii > size(d.coords)[1])  &&  (return nothing)
    return (d.coords[ii, :], ii+1)
end

coords(d::Domain{N}) where N = [d.coords[:, i] for i in 1:N]

struct CartesianDomain{N} <: AbstractDomain{N}
    axis::NTuple{N, Vector{Float64}}
    roi::Vector{Int}
    ldomain::Domain{N}

    function CartesianDomain(axis::Vararg{AbstractVector{T},N}; roi=nothing) where {T <: Real, N}
        @assert N >= 2 "A cartesian domain requires at least 2 dimensions"
        len = prod(length.(axis))
        ss = tuple(length.(axis)...)
        isnothing(roi)  &&  (roi = collect(1:len))

        # Prepare corresponding linear domain
        mat = Matrix{Float64}(undef, length(roi), N)
        ci = Tuple.(CartesianIndices(ss))[roi]
        for i = 1:N
            mat[:, i] .= axis[i][getindex.(ci, i)]
        end
        ldomain = Domain(mat)
        return new{N}(deepcopy(axis), roi, ldomain)
    end

    function CartesianDomain(lengths::Vararg{T,N}; kw...) where {T <: Integer, N}
        @assert all(lengths .>= 1)
        axis = [collect(1.:lengths[i]) for i in 1:N]
        return CartesianDomain(axis...; kw...)
    end
end

# Forward methods to ldomain field
ndims(d::CartesianDomain{N}) where N = N
length(d::CartesianDomain) = length(d.ldomain)
size(d::CartesianDomain) = size(d.ldomain)
getindex(d::CartesianDomain, ii) = getindex(d.ldomain, ii)
iterate(d::CartesianDomain, args...) = iterate(d.ldomain, args...)
coords(d::CartesianDomain) = coords(d.ldomain)

# Cartesian-only methods
orig_size(d::CartesianDomain) = tuple(length.(d.axis)...)
axis(d::CartesianDomain, dim) = d.axis[dim]
roi(d::CartesianDomain) = d.roi

# ====================================================================
# Measures and Counts types
#
abstract type AbstractData{T,N} end

struct Measures{N} <: AbstractData{Float64,N}
    domain::AbstractDomain{N}
    val::Array{Float64,N}
    unc::Array{Float64,N}

    function Measures(domain::Domain{N}, val::AbstractArray{T, N}, unc::AbstractArray{T, N}) where {T <: Real, N}
        @assert length(domain) == length(val) "Domain and dataset have incompatible lengths"
        @assert length(unc) == length(val) "Dataset values and uncertainties have incompatible lengths"
        return new{N}(deepcopy(domain), deepcopy(val), deepcopy(unc))
    end

    function Measures(domain::CartesianDomain{N}, val::AbstractArray{T, N}, unc::AbstractArray{T, N}) where {T <: Real, N}
        @assert prod(orig_size(domain)) == length(val) "Domain and dataset have incompatible lengths"
        return new{N}(deepcopy(domain), deepcopy(val), deepcopy(unc))
    end
end

Measures(domain::AbstractDomain{N}, val::AbstractArray{T, N}, unc::T) where {T <: Real, N} =
    Measures(domain, val, fill(unc, size(val)))

Measures(domain::AbstractDomain{N}, val::AbstractArray{T, N}) where {T <: Real, N} =
    Measures(domain, val, one(Float64))

struct Counts{N} <: AbstractData{Int,N}
    domain::AbstractDomain{N}
    val::Array{Int,N}

    function Counts(domain::Domain{N}, val::AbstractArray{T, N}) where {T <: Integer, N}
        @assert length(domain) == length(val) "Domain and dataset have incompatible lengths"
        new{N}(deepcopy(domain), deepcopy(val))
    end

    function Counts(domain::CartesianDomain{N}, val::AbstractArray{T, N}) where {T <: Real, N}
        @assert prod(orig_size(domain)) == length(val) "Domain and dataset have incompatible lengths"
        return new{N}(deepcopy(domain), deepcopy(val))
    end
end

ndims(d::Measures{N}) where N = N
size(d::Measures) = size(d.val)
length(d::Measures) = length(d.val)
uncerts(d::Measures) = d.uncerts

ndims(d::Counts{N}) where N = N
size(d::Counts) = size(d.val)
iterate(d::Counts, args...) = iterate(d.val, args...)

# ====================================================================
# Methods to "flatten" a multidimensional object into a 1D one
#
function flatten(data::Measures{N}) where N
    (N == 1)  &&  return data
    if isa(data.domain, CartesianDomain)
        return Measures(data.val[roi(data.domain)], data.unc[roi(data.domain)])
    end
    return Measures(data.val[:], data.unc[:])
end

function flatten(data::Counts{N}) where N
    (N == 1)  &&  return data
    if isa(data.domain, CartesianDomain)
        return Counts(data.val[roi(data.domain)])
    end
    return Counts(data.val[:])
end


text:
## header

mutable struct pthread_attr_t
    _pthread_attr_t::NTuple{56, UInt8}
end

const pthread_attr_t = pthread_attr_t

mutable struct pthread_mutex_t
    _pthread_mutex_t::NTuple{40, UInt8}
end

mutable struct pthread_mutexattr_t
    _pthread_mutexattr_t::NTuple{4, UInt8}
end

mutable struct pthread_cond_t
    _pthread_cond_t::NTuple{48, UInt8}
end

mutable struct pthread_condattr_t
    _pthread_condattr_t::NTuple{4, UInt8}
end

const pthread_key_t = UInt32
const pthread_once_t = Cint

mutable struct pthread_rwlock_t
    _pthread_rwlock_t::NTuple{56, UInt8}
end

mutable struct pthread_rwlockattr_t
    _pthread_rwlockattr_t::NTuple{8, UInt8}
end

const pthread_spinlock_t = Cint

mutable struct pthread_barrier_t
    _pthread_barrier_t::NTuple{32, UInt8}
end

mutable struct pthread_barrierattr_t
    _pthread_barrierattr_t::NTuple{4, UInt8}
end

mutable struct random_data
    fptr::Ptr{Int32}
    rptr::Ptr{Int32}
    state::Ptr{Int32}
    rand_type::Cint
    rand_deg::Cint
    rand_sep::Cint
    end_ptr::Ptr{Int32}
end

mutable struct drand48_data
    __x::NTuple{3, UInt16}
    __old_x::NTuple{3, UInt16}
    __c::UInt16
    __init::UInt16
    __a::Culonglong
end

mutable struct matd_t
    nrows::UInt32
    ncols::UInt32
    data::Ptr{Cdouble}
end

mutable struct matd_svd_t
    U::Ptr{matd_t}
    S::Ptr{matd_t}
    V::Ptr{matd_t}
end

mutable struct matd_plu_t
    singular::Cint
    piv::Ptr{UInt32}
    pivsign::Cint
    lu::Ptr{matd_t}
end

mutable struct matd_chol_t
    is_spd::Cint
    u::Ptr{matd_t}
end

const uint8_t = Cuchar
const uint16_t = UInt16
const uint32_t = UInt32
const uint64_t = Culong
const int_least8_t = UInt8
const int_least16_t = Int16
const int_least32_t = Cint
const int_least64_t = Clong
const uint_least8_t = Cuchar
const uint_least16_t = UInt16
const uint_least32_t = UInt32
const uint_least64_t = Culong
const int_fast8_t = UInt8
const int_fast16_t = Clong
const int_fast32_t = Clong
const int_fast64_t = Clong
const uint_fast8_t = Cuchar
const uint_fast16_t = Culong
const uint_fast32_t = Culong
const uint_fast64_t = Culong
const intptr_t = Clong
const uintptr_t = Culong
const intmax_t = Clong
const uintmax_t = Culong

mutable struct image_u8
    width::Int32
    height::Int32
    stride::Int32
    buf::Ptr{UInt8}
end

const image_u8_t = image_u8

mutable struct image_u8x3
    width::Int32
    height::Int32
    stride::Int32
    buf::Ptr{UInt8}
end

const image_u8x3_t = image_u8x3

mutable struct image_u8x4
    width::Int32
    height::Int32
    stride::Int32
    buf::Ptr{UInt8}
end

const image_u8x4_t = image_u8x4

mutable struct image_f32
    width::Int32
    height::Int32
    stride::Int32
    buf::Ptr{Cfloat}
end

const image_f32_t = image_f32

mutable struct image_u32
    width::Int32
    height::Int32
    stride::Int32
    buf::Ptr{UInt32}
end

const image_u32_t = image_u32

mutable struct image_u8_lut
    scale::Cfloat
    nvalues::Cint
    values::Ptr{UInt8}
end

const image_u8_lut_t = image_u8_lut

mutable struct zarray
    el_sz::Csize_t
    size::Cint
    alloc::Cint
    data::Ptr{UInt8}
    zarray() = new()
end

const zarray_t = zarray

mutable struct workerpool
end

const workerpool_t = Nothing

const va_list = Cint

mutable struct timezone
    tz_minuteswest::Cint
    tz_dsttime::Cint
end

mutable struct _pthread_cleanup_buffer
    __routine::Ptr{Nothing}
    __arg::Ptr{Nothing}
    __canceltype::Cint
    __prev::Ptr{_pthread_cleanup_buffer}
end

mutable struct quad
    p::NTuple{4, NTuple{2, Cfloat}}
    reversed_border::Bool
    H::Ptr{matd_t}
    Hinv::Ptr{matd_t}
end

mutable struct timeprofile
    utime::Int64
    stamps::Ptr{zarray_t}
end

const timeprofile_t = timeprofile

#TODO: untested
mutable struct apriltag_family
    ncodes::UInt32
    codes::Ptr{UInt64}
    width_at_border::Cint
    total_width::Cint
    reversed_border::Bool
    nbits::UInt32
    bit_x::Ptr{UInt32}
    bit_y::Ptr{UInt32}
    h::UInt32
    name::Cstring
    impl::Ptr{Nothing}
end

const apriltag_family_t = apriltag_family

mutable struct apriltag_quad_thresh_params
    min_cluster_pixels::Cint
    max_nmaxima::Cint
    cos_critical_rad::Cfloat
    max_line_fit_mse::Cfloat
    min_white_black_diff::Cint
    deglitch::Cint
end

mutable struct apriltag_detector
# #############################################################################
# User-configurable parameters.
    # How many threads should be used?
    nthreads::Cint
    # // detection of quads can be done on a lower-resolution image,
    # // improving speed at a cost of pose accuracy and a slight
    # // decrease in detection rate. Decoding the binary payload is
    # // still done at full resolution. .
    quad_decimate::Cfloat
    #  What Gaussian blur should be applied to the segmented image
    #  (used for quad detection?)  Parameter is the standard deviation
    #  in pixels.  Very noisy images benefit from non-zero values
    #  (e.g. 0.8).
    quad_sigma::Cfloat
    #  When non-zero, the edges of the each quad are adjusted to "snap
    #  to" strong gradients nearby. This is useful when decimation is
    #  employed, as it can increase the quality of the initial quad
    #  estimate substantially. Generally recommended to be on (1).
    #  Very computationally inexpensive. Option is ignored if quad_decimate = 1.
    refine_edges::Cint

    # How much sharpening should be done to decoded images? This
    # can help decode small tags but may or may not help in odd
    # lighting conditions or low light conditions.
    #
    # The default value is 0.25.
    decode_sharpening::Cdouble

    #  When non-zero, write a variety of debugging images to the
    #  current working directory at various stages through the
    #  detection process. (Somewhat slow).
    # int debug;
    debug::Cint

    # qtp::apriltag_quad_thresh_params
        min_cluster_pixels::Cint
        max_nmaxima::Cint
        critical_rad::Cfloat
        max_line_fit_mse::Cfloat
        min_white_black_diff::Cint
        deglitch::Cint

    tp::Ptr{timeprofile_t}
    nedges::UInt32
    nsegments::UInt32
    nquads::UInt32
    tag_families::Ptr{zarray_t}
    wp::Ptr{workerpool_t}

    # mutex::pthread_mutex_t
        mutex::NTuple{40, UInt8}

    apriltag_detector() = new()
end

const apriltag_detector_t = apriltag_detector

mutable struct apriltag_detection
    family::Ptr{apriltag_family_t}
    id::Cint
    hamming::Cint
    # goodness::Cfloat
    decision_margin::Cfloat
    H::Ptr{matd_t}
    c::NTuple{2, Cdouble}
    p::NTuple{4, NTuple{2, Cdouble}}
    apriltag_detection() = new()
end

const apriltag_detection_t = apriltag_detection

# Tag families constructors and destructors
"""
	tag36h11_create()
Create a AprilTag family object for tag36h11 with all fields set to default value.
"""
function tag36h11_create()
    ccall((:tag36h11_create, libapriltag), Ptr{apriltag_family_t}, ())
end

"""
	tag36h11_destroy(tf)
Destroy the AprilTag family object.
"""
function tag36h11_destroy(tf)
    ccall((:tag36h11_destroy, libapriltag), Nothing, (Ptr{apriltag_family_t},), tf)
end

#NOTE Not in apriltag 3
# function tag36h10_create()
#     ccall((:tag36h10_create, libapriltag), Ptr{apriltag_family_t}, ())
# end
#
# function tag36h10_destroy(tf)
#     ccall((:tag36h10_destroy, libapriltag), Nothing, (Ptr{apriltag_family_t},), tf)
# end

function tag25h9_create()
    ccall((:tag25h9_create, libapriltag), Ptr{apriltag_family_t}, ())
end

function tag25h9_destroy(tf)
    ccall((:tag25h9_destroy, libapriltag), Nothing, (Ptr{apriltag_family_t},), tf)
end

#NOTE Not in apriltag 3
# function tag25h7_create()
#     ccall((:tag25h7_create, libapriltag), Ptr{apriltag_family_t}, ())
# end

# function tag25h7_destroy(tf)
#     ccall((:tag25h7_destroy, libapriltag), Nothing, (Ptr{apriltag_family_t},), tf)
# end

function tag16h5_create()
    ccall((:tag16h5_create, libapriltag), Ptr{apriltag_family_t}, ())
end

function tag16h5_destroy(tf)
    ccall((:tag16h5_destroy, libapriltag), Nothing, (Ptr{apriltag_family_t},), tf)
end

"""
	apriltag_detector_create()
Create a AprilTag Detector object with all fields set to default value.
"""
function apriltag_detector_create()
    ccall((:apriltag_detector_create, libapriltag), Ptr{apriltag_detector_t}, ())
end

function apriltag_detector_add_family_bits(td, fam, bits_corrected::Cint)
    ccall((:apriltag_detector_add_family_bits, libapriltag), Nothing, (Ptr{apriltag_detector_t}, Ptr{apriltag_family_t}, Cint), td, fam, bits_corrected)
end

"""
	apriltag_detector_add_family(tag_detector, tag_family)
Add a tag family to an AprilTag Detector object.
The caller still "owns" the family and a single instance should only be provided to one apriltag detector instance.
"""
function apriltag_detector_add_family(td, fam)
    # ccall((:apriltag_detector_add_family, libapriltag), Nothing, (Ptr{apriltag_detector_t}, Ptr{apriltag_family_t}), td, fam)
    apriltag_detector_add_family_bits(td, fam, Int32(2))
end

function apriltag_detector_remove_family(td, fam)
    ccall((:apriltag_detector_remove_family, libapriltag), Nothing, (Ptr{apriltag_detector_t}, Ptr{apriltag_family_t}), td, fam)
end

function apriltag_detector_clear_families(td)
    ccall((:apriltag_detector_clear_families, libapriltag), Nothing, (Ptr{apriltag_detector_t},), td)
end

function apriltag_detector_destroy(td)
    ccall((:apriltag_detector_destroy, libapriltag), Nothing, (Ptr{apriltag_detector_t},), td)
end

"""
	apriltag_detector_detect(tag_detector, image)
Detect tags from an image and return an array of apriltag_detection_t*.
You can use apriltag_detections_destroy to free the array and the detections it contains, or call
detection_destroy and zarray_destroy yourself.
"""
function apriltag_detector_detect(td, im_orig)
    ccall((:apriltag_detector_detect, libapriltag), Ptr{zarray_t}, (Ptr{apriltag_detector_t}, Ptr{image_u8_t}), td, Ref(im_orig))
end

"""
	threadcall_apriltag_detector_detect(tag_detector, image)
*Experimental* call apriltag_detector_detect in a seperate thread using the experimantal `@threadcall`
Detect tags from an image and return an array of apriltag_detection_t*.
You can use apriltag_detections_destroy to free the array and the detections it contains, or call
detection_destroy and zarray_destroy yourself.
"""
function threadcall_apriltag_detector_detect(td, im_orig)
    @threadcall((:apriltag_detector_detect, libapriltag), Ptr{zarray_t}, (Ptr{apriltag_detector_t}, Ptr{image_u8_t}), td, Ref(im_orig))
end

function apriltag_detection_destroy(det)
    ccall((:apriltag_detection_destroy, libapriltag), Nothing, (Ptr{apriltag_detection_t},), det)
end

function apriltag_detections_destroy(detections)
    ccall((:apriltag_detections_destroy, libapriltag), Nothing, (Ptr{zarray_t},), detections)
end

function apriltag_to_image(fam, idx::Cint)
    ccall((:apriltag_to_image, libapriltag), Ptr{image_u8_t}, (Ptr{apriltag_family_t}, Cint), fam, idx)
end

#common

# matd_t *homography_to_pose(const matd_t *H, double fx, double fy, double cx, double cy)
function homography_to_pose(H, fx, fy, cx, cy)
    ccall((:homography_to_pose, libapriltag), Ptr{matd_t}, (Ptr{matd_t}, Cdouble, Cdouble, Cdouble, Cdouble), H, fx, fy, cx, cy)
end

"""
    $TYPEDEF

new pose estimation from apriltag 3

Notes:
- This low-level `ccall` wrapped C-library type uses the convention (i.e. the camera-frame):
  - `fx == f_width`,
  - `cy == c_height`, and
  - C-library camara coordinate system: camera looking along positive Z axis with `x` to the right and `y` down.
    - C-library internally follows: https://docs.opencv.org/3.4/d9/d0c/group__calib3d.html
- Images.jl uses `::Array` in Julia as column-major (i.e. vertical major) convention, that is `size(img) == (480, 640)`
  - Axes start top left-corner of the image plane (i.e. the image-frame):
  - `width` is from left to right,
  - `height` is from top downward.
"""
mutable struct apriltag_detection_info_t
    det::Ptr{apriltag_detection_t}
    tagsize::Cdouble
    fx::Cdouble
    fy::Cdouble
    cx::Cdouble
    cy::Cdouble
end

mutable struct apriltag_pose_t
    R::Ptr{matd_t}
    t::Ptr{matd_t}
    apriltag_pose_t() = new()
end

#special matd_t to help with construction and usage
mutable struct Matd4x4
    nrows::UInt32
    ncols::UInt32
    data::NTuple{16,Cdouble}
    Matd4x4() = new(4,4,NTuple{16,Cdouble}(Matrix{Float64}(undef,4,4)))
    Matd4x4(M::Matrix{Float64}) = new(4,4,NTuple{16,Cdouble}(M'))#TODO test if order is correct
end

mutable struct Matd3x3
    nrows::UInt32
    ncols::UInt32
    data::NTuple{9,Cdouble}
    Matd3x3() = new(3,3,NTuple{9,Cdouble}(Matrix{Float64}(undef,3,3)))
    Matd3x3(M::Matrix{Float64}) = new(3,3,NTuple{9,Cdouble}(M'))#TODO test if order is correct
end

mutable struct Matd3x1
    nrows::UInt32
    ncols::UInt32
    data::NTuple{3,Cdouble}
    Matd3x1() = new(3,1,NTuple{3,Cdouble}(Matrix{Float64}(undef,3,1)))
    Matd3x1(V::Vector{Float64}) = new(3,1,NTuple{3,Cdouble}(V))
end

# void estimate_pose_for_tag_homography(apriltag_detection_info_t* info, apriltag_pose_t* pose);
function estimate_pose_for_tag_homography(info, pose)
    ccall((:estimate_pose_for_tag_homography, libapriltag), Nothing, (Ptr{apriltag_detection_info_t}, Ptr{apriltag_pose_t}), Ref(info), Ref(pose))
end

# double orthogonal_iteration(matd_t** v, matd_t** p, matd_t** t, matd_t** R, int n_points, int n_steps)
function orthogonal_iteration(v::NTuple{4,Matd3x1}, p::NTuple{4,Matd3x1}, t::NTuple{1,Matd3x1}, R::NTuple{1,Matd3x3}, n_points::Int, n_steps::Int)
# function orthogonal_iteration(v, p, t, R, n_points::Int, n_steps::Int)

    vp = Base.unsafe_convert(Ptr{Nothing}, Ref(v))
    pp = Base.unsafe_convert(Ptr{Nothing}, Ref(p))
    tp = Base.unsafe_convert(Ptr{Nothing}, Ref(t))
    Rp = Base.unsafe_convert(Ptr{Nothing}, Ref(R))

    ccall((:orthogonal_iteration, libapriltag), Cdouble,
            (Ptr{Ptr{matd_t}}, Ptr{Ptr{matd_t}}, Ptr{Ptr{matd_t}}, Ptr{Ptr{matd_t}}, Cint, Cint),
            vp, pp, tp, Rp, n_points, n_steps)
end

# matd_t* fix_pose_ambiguities(matd_t** v, matd_t** p, matd_t* t, matd_t* R, int n_points) {
function fix_pose_ambiguities(v::NTuple{4,Matd3x1}, p::NTuple{4,Matd3x1}, t::Matd3x1, R::Matd3x3, n_points::Int)

    vp = Base.unsafe_convert(Ptr{Nothing}, Ref(v))
    pp = Base.unsafe_convert(Ptr{Nothing}, Ref(p))
    tp = Base.unsafe_convert(Ptr{Nothing}, Ref(t))
    Rp = Base.unsafe_convert(Ptr{Nothing}, Ref(R))

    ccall((:fix_pose_ambiguities, libapriltag), Ptr{Matd3x3},
            (Ptr{Ptr{matd_t}}, Ptr{Ptr{matd_t}}, Ptr{matd_t}, Ptr{matd_t}, Cint),
            vp, pp, tp, Rp, n_points)
end


text:
# Transcode
# =========

"""
    transcode(::Type{C}, data::Vector{UInt8})::Vector{UInt8} where C<:Codec

Transcode `data` by applying a codec `C()`.

Note that this method does allocation and deallocation of `C()` in every call,
which is handy but less efficient when transcoding a number of objects.
`transcode(codec, data)` is a recommended method in terms of performance.

Examples
--------

```julia
julia> using CodecZlib

julia> data = b"abracadabra";

julia> compressed = transcode(ZlibCompressor, data);

julia> decompressed = transcode(ZlibDecompressor, compressed);

julia> String(decompressed)
"abracadabra"

```
"""
function Base.transcode(::Type{C}, data::ByteData) where C<:Codec
    codec = C()
    initialize(codec)
    try
        return transcode(codec, data)
    finally
        finalize(codec)
    end
end

"""
    transcode(codec::Codec, data::Vector{UInt8})::Vector{UInt8}

Transcode `data` by applying `codec`.

Note that this method does not initialize or finalize `codec`. This is
efficient when you transcode a number of pieces of data, but you need to call
[`TranscodingStreams.initialize`](@ref) and
[`TranscodingStreams.finalize`](@ref) explicitly.

Examples
--------

```julia
julia> using CodecZlib

julia> data = b"abracadabra";

julia> codec = ZlibCompressor();

julia> TranscodingStreams.initialize(codec)

julia> compressed = transcode(codec, data);

julia> TranscodingStreams.finalize(codec)

julia> codec = ZlibDecompressor();

julia> TranscodingStreams.initialize(codec)

julia> decompressed = transcode(codec, compressed);

julia> TranscodingStreams.finalize(codec)

julia> String(decompressed)
"abracadabra"

```
"""
function Base.transcode(codec::Codec, data::ByteData)
    # Add `minoutsize` because `transcode` will be called at least two times.
    buffer2 = Buffer(
        expectedsize(codec, Memory(data)) + minoutsize(codec, Memory(C_NULL, 0)))
    mark!(buffer2)
    stream = TranscodingStream(codec, devnull, State(Buffer(data), buffer2); initialized=true)
    write(stream, TOKEN_END)
    return takemarked!(buffer2)
end


text:
##########################################################################
# Copyright 2017 Samuel Ridler.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
##########################################################################

using Random

@testset "grid" begin
	# create map
	mp = Map(0, 1, 0, 1, 10, 10) # Map(xMin, xMax, yMin, yMax, xScale, yScale)

	# create nodes
	numNodes = 100
	offRoadAccessProb = 0.8
	rng = Random.MersenneTwister(0)
	nodes = [Node() for i = 1:numNodes]
	for i = 1:numNodes
		nodes[i].index = i
		nodes[i].location = randLocation(mp; rng = rng)
		nodes[i].offRoadAccess = rand() < offRoadAccessProb
	end
	numOffRoadAccessNodes = sum(node -> node.offRoadAccess, nodes)
	@assert(numOffRoadAccessNodes < numNodes) # otherwise not all cases will be tested

	# create grids, add nodes to grid
	nx = ny = 10
	grid1 = Grid(mp, nx, ny) # grid with only nodes that have offRoadAccess == true
	JEMSS.gridPlaceNodes!(mp, grid1, nodes)
	grid2 = Grid(mp, nx, ny) # grid with all nodes
	JEMSS.gridPlaceNodes!(mp, grid2, nodes; offRoadAccessRequired = false)

	# test number of nodes added to each grid
	@test sum(gr -> length(gr.nodeIndices), grid1.rects) == numOffRoadAccessNodes
	@test sum(gr -> length(gr.nodeIndices), grid2.rects) == numNodes

	# each node should be closest to itself (for grid2, not grid1)
	@test all(node -> findNearestNode(mp, grid2, nodes, node.location)[1] == node.index, nodes)

	# check that grid rectangles have some different numbers of nodes, so that different cases will be tested
	@assert(all(grid -> any(rect -> length(rect.nodeIndices) == 0, grid.rects), (grid1, grid2)))
	@assert(all(grid -> any(rect -> length(rect.nodeIndices) == 1, grid.rects), (grid1, grid2)))
	@assert(all(grid -> any(rect -> length(rect.nodeIndices) >= 2, grid.rects), (grid1, grid2)))

	# compare findNearestNode function with brute-force method (no grid)
	grid1pass = grid2pass = true
	for x = range(0; length = 10*nx+1, stop = 1), y = range(0; length = 10*ny+1, stop = 1)
		location = Location(x,y)
		grid1pass &= findNearestNode(mp, grid1, nodes, location) == findNearestNode(mp, nodes, location; offRoadAccessRequired = true)
		grid2pass &= findNearestNode(mp, grid2, nodes, location) == findNearestNode(mp, nodes, location; offRoadAccessRequired = false)
	end
	@test grid1pass
	@test grid2pass

	# plot of nearest node from sample locations
	# this is not a test, but I wanted visual confirmation that the findNearestNode function was working
	false && begin
		# sample locations
		xs = range(0; length = 100*nx+1, stop = 1)
		ys = range(0; length = 100*ny+1, stop = 1)
		z = zeros(Int, length(ys), length(xs));
		for (ix, x) in enumerate(xs), (iy, y) in enumerate(ys)
			location = Location(x,y)
			z[iy,ix] = findNearestNode(mp, grid2, nodes, location)[1]
		end

		using Plots
		Plots.plotly()
		plt = heatmap(xs, ys, z, aspect_ratio=1);
		x = [node.location.x for node in nodes]
		y = [node.location.y for node in nodes]
		Plots.scatter!(plt,x,y)
	end
end


text:
"""
    KmeansAlg(m::Int, metric::SemiMetric=SqEuclidean(); nMarkov = 10, tol = 1e-3)

k-Means [1] initialization on the data `X` taking `m` inducing points.
The seeding is computed via [2], `nMarkov` gives the number of MCMC steps for the seeding.

## Arguments
- `k::Int` : Number of inducing points
- `metric::SemiMetric` : Metric used to compute the distances for the k-means algorithm

## Keyword Arguments
- `nMarkov::Int` : Number of random steps for the seeding
- `tol::Real` : Tolerance for the `kmeans` algorithm

[1] Arthur, D. & Vassilvitskii, S. k-means++: The advantages of careful seeding. in Proceedings of the eighteenth annual ACM-SIAM symposium on Discrete algorithms 1027–1035 (Society for Industrial and Applied Mathematics, 2007).
[2] Bachem, O., Lucic, M., Hassani, S. H. & Krause, A. Fast and Provably Good Seedings for k-Means. Advances in Neural Information Processing Systems 29 55--63 (2016) doi:10.1109/tmtt.2005.863818.
"""
struct KmeansAlg{Tm,T} <: OffIPSA
    m::Int
    metric::Tm
    nMarkov::Int
    tol::T
    function KmeansAlg(
        m::Int, metric::SemiMetric=SqEuclidean(); nMarkov::Int=10, tol::T=1e-3
    ) where {T<:Real}
        m > 0 || throw(ArgumentError("The number of inducing points m should be positive"))
        tol > 0 || throw(ArgumentError("The tolerance tol should be positive"))
        return new{typeof(metric),T}(m, metric, nMarkov, tol)
    end
end

function inducingpoints(
    rng::AbstractRNG, alg::KmeansAlg, X::AbstractVector; weights=nothing, kwargs...
)
    if alg.m >= length(X)
        return edge_case(alg.m, length(X), X)
    end
    return kmeans_ip(
        rng, X, alg.m, alg.metric; nMarkov=alg.nMarkov, weights=weights, tol=alg.tol
    )
end

function Base.show(io::IO, alg::KmeansAlg)
    return print(io, "k-Means Selection of Inducing Points (m : $(alg.m))")
end

#Return K inducing points from X, m being the number of Markov iterations for the seeding
function kmeans_ip(
    rng::AbstractRNG,
    X::AbstractVector,
    nC::Int,
    metric::SemiMetric;
    nMarkov::Int=10,
    weights=nothing,
    tol=1e-3,
)
    C = kmeans_seeding(rng, X, nC, metric, nMarkov)
    C = reduce(hcat, C)
    kmeans!(reduce(hcat, X), C; weights=weights, tol=tol, distance=metric)
    return ColVecs(C)
end

"""Fast and efficient seeding for KMeans based on [`Fast and Provably Good Seeding for k-Means](https://las.inf.ethz.ch/files/bachem16fast.pdf)"""
function kmeans_seeding(
    rng::AbstractRNG,
    X::AbstractVector, # Input data
    nC::Integer, # Number of centroids
    metric::SemiMetric, # Metric used
    nMarkov::Integer, # Number of Markov iterations
)
    nSamples = length(X) # Number of input samples
    # Preprocessing,
    init = rand(rng, 1:nSamples) # Sample first random center
    C = collect.(X[init:init]) # Initialize the collection of centroids
    q = vec(pairwise(metric, X, C)) # Create the pairwise values between the data and the first centroid
    q = Weights(q / sum(q) .+ 1.0 / (2 * nSamples), 1) # Create weights to work with
    for i in 2:nC
        x = X[sample(rng, q)] # weighted sampling,
        mindist = mindistance(metric, x, C) # Find the closest centroid to the random sample x
        for j in 2:nMarkov # Iterate over nMarkov iterations
            y = X[sample(rng, q)] # Draw a new sample
            dist = mindistance(metric, y, C) # Find the closest centroid to the random sample y
            if (dist / mindist) > rand(rng) # Acceptance step to see if y is better than x
                x = y
                mindist = dist
            end
        end
        push!(C, x)
    end
    return C
end


text:
module IntervalUnionArithmetic

using Reexport

@reexport using IntervalArithmetic

import Base: getindex, ∪, intersect, \, in, isequal, isempty
import Base: +, -, *, /, min, max, ^, log, <, >, exp, sin, cos, tan, sqrt, ==
import IntervalArithmetic: union, ∪, hull, bisect, intersect, ⊆, ⊂, setdiff

export
    IntervalUnion, intervalUnion, remove_empties, condense, closeGaps!,
    iscondensed, condense_weak, iscondensed_weak, env, hull, complement, left, right, bisect,
    intersect, setdiff, \, ⊆, ⊂, isequal, in, ==, isempty

global MAXINTS = [100]

include("interval_unions.jl")
include("arithmetic.jl")
include("utilities.jl")
include("set_operations.jl")

end # module


text:
"
Author: Daniel O'Connor
Email: dt.oconnor@outlook.com

File used to profile and benchmark the algorithms implmented.
"

include("samplers/tools.jl")
include("samplers/SA.jl")
include("samplers/SVMC.jl")
include("samplers/SVMC_TF.jl")
include("samplers/pimc.jl")
using Random
using ..SA
using ..SVMC
using ..SVMCTF
using ..PIMC
using ..SamplingTools
using LinearAlgebra
using Profile
using BenchmarkTools
using NPZ
using Interpolations

function main()
    # test 0 =  SA, test 1 = SVMC, test 2 = SVMCTF, test 3 = PIMC
    test = 3

    # simple PFC problem
    d = 0.1
    h = Dict([(1, -1.0), (2, 1.0 - d), (3, 1.0-d), (4, -1.0)])
    J = Dict([((1, 2), -1.0), ((2, 3), -1.0), ((3, 4), -1.0)])

    num_samples = 1000
    num_sweeps = 10000
    num_slices = 20
    nspins = length(h)
    temp = 0.01226 * 1.380649e-23 / (1e9 * 6.62607015e-34)

    a_sch= 1.0 .- collect(LinRange(0, 1, num_sweeps))
    b_sch= collect(LinRange(0, 1, num_sweeps))

    # schedule_data = npzread("data\\DWAVE_schedule.npz")
    # a_fn = LinearInterpolation(schedule_data[:, 1], 0.5.*schedule_data[:, 2])
    # b_fn = LinearInterpolation(schedule_data[:, 1], 0.5.*schedule_data[:, 3])
    # a_sch = a_fn.(collect(LinRange(0, 1, num_sweeps)))
    # b_sch = b_fn.(collect(LinRange(0, 1, num_sweeps)))

    a_sch_dict = Dict([(1, deepcopy(a_sch)), (2, 0.0 .* deepcopy(a_sch)),
                       (3, deepcopy(a_sch)), (4, deepcopy(a_sch)),])
    b_sch_dict = Dict([(1, deepcopy(b_sch)), (2, deepcopy(b_sch)),
                       (3, deepcopy(b_sch)), (4, deepcopy(b_sch)),])

    # initialize the problem parameter
    JM, M, degs = SamplingTools.hj_info(h, J)
    nbs_idxs, nbs_vals = SamplingTools.generate_nbr(nspins, M, maximum(degs) + 1)

    hot_t, cold_t = SamplingTools.default_temp_range(JM)
    temp_schedule = collect(LinRange(hot_t, cold_t, num_sweeps))

    # initialize state
    initial_state = Array{Float64, 1}([1.0, -1.0, -1.0, -1.0])
    if test == 0
        # spin_vectors = [deepcopy(initial_state) for i =1:num_samples]
        initial_state = rand([-1.0, 1.0], nspins, num_samples)
        spin_vectors = [initial_state[:, i] for i =1:num_samples]
        _fn = SA.anneal!
        args_in = [nbs_idxs, nbs_vals, temp_schedule]
    elseif test == 1 || test == 2
        # spin_vectors = [acos.(deepcopy(initial_state)) for i =1:num_reads]
        spin_vectors = fill(0.5 * pi, nspins, num_samples)
        spin_vectors = [spin_vectors[:, i] for i =1:num_samples]
        if test == 1
            _fn = SVMC.anneal!
        else
            _fn = SVMCTF.anneal!
        end
        args_in = [nbs_idxs, nbs_vals, a_sch, b_sch, temp]
    elseif  test == 3
        # spin_array = repeat(initial_state', num_slices)'
        # spin_vectors = [copy(spin_array) for i =1:num_samples]
        initial_state = rand([-1.0, 1.0], nspins, num_slices, num_samples)
        spin_vectors = [initial_state[:, :, i] for i =1:num_samples]
        _fn = PIMC.anneal!
        args_in = [nbs_idxs, nbs_vals, a_sch, b_sch, temp]
    end

    function run()
        @inbounds for i = 1:num_samples
            _fn(spin_vectors[i], args_in...)
        end
    end

    run()

    if test == 0
        print(discrete_aggregation(spin_vectors, num_samples, collect(UpperTriangular(JM))))
    elseif test == 1 || test == 2
        print(rotor_aggregation(spin_vectors, num_samples, collect(UpperTriangular(JM))))
    elseif  test == 3
        print(pimc_aggregation(spin_vectors, num_samples, collect(UpperTriangular(JM))))
    end

    p = @profiler run()
    Profile.print()

    return  @benchmark $run()
end

max_stars_repo_name: dtoconnor/Annealing-Samplers

max_issues_repo_name: dtoconnor/Annealing-Samplers

max_forks_repo_name: dtoconnor/Annealing-Samplers


text:
struct PCSAFTParam <: EoSParam
    Mw::SingleParam{Float64}
    segment::SingleParam{Float64}
    sigma::PairParam{Float64}
    epsilon::PairParam{Float64}
    epsilon_assoc::AssocParam{Float64}
    bondvol::AssocParam{Float64}
end

abstract type PCSAFTModel <: SAFTModel end
@newmodel PCSAFT PCSAFTModel PCSAFTParam

export PCSAFT
function PCSAFT(components::Array{String,1}; idealmodel=BasicIdeal, userlocations::Array{String,1}=String[], verbose=false)
    params = getparams(components, ["SAFT/PCSAFT","properties/molarmass.csv"]; userlocations=userlocations, verbose=verbose)
    segment = params["m"]
    k = params["k"]
    Mw = params["Mw"]
    params["sigma"].values .*= 1E-10
    sigma = sigma_LorentzBerthelot(params["sigma"])
    epsilon = epsilon_LorentzBerthelot(params["epsilon"], k)
    epsilon_assoc = params["epsilon_assoc"]
    bondvol = params["bondvol"]
    sites = SiteParam(Dict("e" => params["n_e"], "H" => params["n_H"]))

    packagedparams = PCSAFTParam(Mw,segment, sigma, epsilon, epsilon_assoc, bondvol)
    references = ["10.1021/ie0003887", "10.1021/ie010954d"]

    model = PCSAFT(packagedparams, sites, idealmodel; references=references, verbose=verbose)
    return model
end

function a_res(model::PCSAFTModel, V, T, z)
    return @f(a_hc) + @f(a_disp) + @f(a_assoc)
end

function a_hc(model::PCSAFTModel, V, T, z)
    x = z/∑(z)
    m = model.params.segment.values
    m̄ = ∑(x .* m)
    return m̄*@f(a_hs) - ∑(x[i]*(m[i]-1)*log(@f(g_hs,i,i)) for i ∈ @comps)
end

function a_disp(model::PCSAFTModel, V, T, z)
    x = z/∑(z)
    m = model.params.segment.values
    m̄ = ∑(x .* m)
    return -2*π*N_A*∑(z)/V*@f(I,1)*@f(m2ϵσ3, 1) - π*m̄*N_A*∑(z)/V*@f(C1)*@f(I,2)*@f(m2ϵσ3,2)
end

function d(model::PCSAFTModel, V, T, z, i)
    ϵii = model.params.epsilon.diagvalues[i]
    σii = model.params.sigma.diagvalues[i]
    return σii * (1 - 0.12exp(-3ϵii/T))
end

function ζ(model::PCSAFTModel, V, T, z, n)
    ∑z = ∑(z)
    x = z * (one(∑z)/∑z)
    m = model.params.segment.values
    res = N_A*∑z*π/6/V * ∑((x[i]*m[i]*@f(d,i)^n for i ∈ @comps))
end

function g_hs(model::PCSAFTModel, V, T, z, i, j)
    di = @f(d,i)
    dj = @f(d,j)
    ζ2 = @f(ζ,2)
    ζ3 = @f(ζ,3)
    return 1/(1-ζ3) + di*dj/(di+dj)*3ζ2/(1-ζ3)^2 + (di*dj/(di+dj))^2*2ζ2^2/(1-ζ3)^3
end

function a_hs(model::PCSAFTModel, V, T, z)
    ζ0 = @f(ζ,0)
    ζ1 = @f(ζ,1)
    ζ2 = @f(ζ,2)
    ζ3 = @f(ζ,3)
    return 1/ζ0 * (3ζ1*ζ2/(1-ζ3) + ζ2^3/(ζ3*(1-ζ3)^2) + (ζ2^3/ζ3^2-ζ0)*log(1-ζ3))
end

function C1(model::PCSAFTModel, V, T, z)
    x = z/∑(z)
    m = model.params.segment.values
    m̄ = ∑(x .* m)
    η = @f(ζ,3)
    return (1 + m̄*(8η-2η^2)/(1-η)^4 + (1-m̄)*(20η-27η^2+12η^3-2η^4)/((1-η)*(2-η))^2)^-1
end

function m2ϵσ3(model::PCSAFTModel, V, T, z, n = 1)
    x = z/∑(z)
    m = model.params.segment.values
    σ = model.params.sigma.values
    ϵ = model.params.epsilon.values
    return ∑(x[i]*x[j]*m[i]*m[j] * (ϵ[i,j]*(1)/T)^n * σ[i,j]^3 for i ∈ @comps, j ∈ @comps)
end

function I(model::PCSAFTModel, V, T, z, n)
    x = z/∑(z)
    m = model.params.segment.values
    m̄ = ∑(x .* m)
    σ = model.params.sigma.values
    η = @f(ζ,3)
    if n == 1
        corr = PCSAFTconsts.corr1
    elseif n == 2
        corr = PCSAFTconsts.corr2
    end
    return ∑((corr[i+1,1] + (m̄-1)/m̄*corr[i+1,2] + (m̄-1)/m̄*(m̄-2)/m̄*corr[i+1,3]) * η^i for i = 0:6)
end

function a_assoc(model::PCSAFTModel, V, T, z)
    x = z/∑(z)
    X_ = @f(X)
    n = model.allcomponentnsites
    return ∑(x[i]*∑(n[i][a] * (log(X_[i][a]) - X_[i][a]/2 + 0.5) for a ∈ @sites(i)) for i ∈ @comps)
end

function X(model::PCSAFTModel, V, T, z)
    _1 = one(V+T+first(z))
    Σz = ∑(z)
    x = z/ Σz
    ρ = N_A* Σz/V
    itermax = 100
    dampingfactor = 0.5
    error = 1.
    tol = model.absolutetolerance
    iter = 1
    X_ = [[_1 for a ∈ @sites(i)] for i ∈ @comps]
    X_old = deepcopy(X_)
    while error > tol
        iter > itermax && error("X has failed to converge after $itermax iterations")
        for i ∈ @comps, a ∈ @sites(i)
            rhs = 1/(1+∑(ρ*x[j]*∑(X_old[j][b]*@f(Δ,i,j,a,b) for b ∈ @sites(j)) for j ∈ @comps))
            X_[i][a] = (1-dampingfactor)*X_old[i][a] + dampingfactor*rhs
        end
        error = sqrt(∑(∑((X_[i][a] - X_old[i][a])^2 for a ∈ @sites(i)) for i ∈ @comps))
        for i = 1:length(X_)
            X_old[i] .= X_[i]
        end
        iter += 1
    end
    return X_
end

function Δ(model::PCSAFTModel, V, T, z, i, j, a, b)
    ϵ_assoc = model.params.epsilon_assoc.values
    κ = model.params.bondvol.values
    σ = model.params.sigma.values
    gij = @f(g_hs,i,j)
    return gij*σ[i,j]^3*(exp(ϵ_assoc[i,j][a,b]/T)-1)*κ[i,j][a,b]
end

const PCSAFTconsts = (
    corr1 =
    [0.9105631445 -0.3084016918 -0.0906148351;
    0.6361281449 0.1860531159 0.4527842806;
    2.6861347891 -2.5030047259 0.5962700728;
    -26.547362491 21.419793629 -1.7241829131;
    97.759208784 -65.255885330 -4.1302112531;
    -159.59154087 83.318680481 13.776631870;
    91.297774084 -33.746922930 -8.6728470368],

    corr2 =
    [0.7240946941 -0.5755498075 0.0976883116;
    2.2382791861 0.6995095521 -0.2557574982;
    -4.0025849485 3.8925673390 -9.1558561530;
    -21.003576815 -17.215471648 20.642075974;
    26.855641363 192.67226447 -38.804430052;
    206.55133841 -161.82646165 93.626774077;
    -355.60235612 -165.20769346 -29.666905585]
)


text:
getTemplate(x,y,z) =  string(x) * "時の" * string(y) * "は" * string(z)

println(getTemplate(12,"気温",22.4))