text:
using Pkg
Pkg.activate("examples")
Pkg.instantiate()

# for det()
using LinearAlgebra

# for importing the data from CSV
using DataFrames
using CSV

# for GLOM
import GPLinearODEMaker; GLOM = GPLinearODEMaker

# for this module
using GLOM_RV_Example; GLOM_RV = GLOM_RV_Example

# for units in orbit fitting functions
using UnitfulAstro, Unitful

# for getting times from HDF5 file the SAFE stats were generated from
using HDF5

# for PCA and std()
using MultivariateStats
using Statistics

#####################################################################
# CTRL+F "CHANGE"TO FIND PLACES WHERE YOU SHOULD MAKE MODIFICATIONS #
#####################################################################

## Problem setup
length(ARGS) > 0 ? n_bs = parse(Int, ARGS[1]) : n_bs = 5
@assert 1 < n_bs < 6

# CHANGE: whether you want to look for a planet that you inject
inject_planet = false
# inject_ks = GLOM_RV.kep_signal(; K=1.0u"m/s", P=sqrt(2)*5u"d", M0=rand()*2*π, ω_or_k=rand()*2*π, e_or_h=0.1)
inject_ks = GLOM_RV.kep_signal(; K=1.0u"m/s", P=sqrt(2)*5u"d", M0=3.2992691080593275, ω_or_k=4.110936513051912, e_or_h=0.1)

# CHANGE: choose a kernel, I suggest 3 for Matern 5/2 or 4 for Quasi-periodic
# kernel
kernel_choice = 3
kernel_names = ["pp", "se", "m52", "qp", "m52_m52", "se_se"]
kernel_name = kernel_names[kernel_choice]
kernel_function, num_kernel_hyperparameters = GLOM.include_kernel(kernel_name)

# CHANGE: the stars rotation rate which is used as the first guess for some GLOM
# hyperparameters and starting point for priors
star_rot_rate = 25.05  # days

# importing SOAP SAFE data
data = CSV.read("timecorssim_SAFE_snr600.csv", DataFrame)[1:2:end, :]

# CHANGE: observation times go here
obs_xs = h5open("D:/Christian/Downloads/res-1000-1years_full_id1.h5")["phases"][1:2:366] * star_rot_rate
# taking out the mean observation times makes the optimization easier for models
# with periodic parameters
GLOM_RV.remove_mean!(obs_xs)

# CHANGE: rvs and their errors go here
obs_rvs = collect(data[!, "b1"])
if inject_planet
    obs_rvs[:] .+= ustrip.(inject_ks.(obs_xs.*u"d"))
end
obs_rvs_err = data[!, "b1"] ./ data[!, "t1"]

## PCA for indicators

find_data_col_ind(col_name) = findfirst(x -> x==col_name, names(data))
b_inds = find_data_col_ind("b0"):(find_data_col_ind("b0") + n_bs - 1)
t_inds = find_data_col_ind("t0") - find_data_col_ind("b0") .+ b_inds

Xte = Matrix(data[!, b_inds])'

# #PPCA version
# M = fit(PPCA, Xte; maxoutdim=2)
# # modifed from transform() at from https://github.com/JuliaStats/MultivariateStats.jl/blob/master/src/ppca.jl
# function transforming_mat(m::PPCA{T}) where {T<:Real}
#     W  = m.W
#     n = outdim(m)
#     M = W'W .+ m.σ² * Matrix{T}(I, n, n)
#     return inv(M)*m.W'
# end
# Z = transforming_mat(M)

# PCA version
M = fit(PCA, Xte; maxoutdim=2)
Z = M.proj'
println("Z: ", Z)

z = MultivariateStats.transform(M, Xte)  # = Z * (Xte .- M.mean)
var_β = ((Matrix(data[!, b_inds]) ./ Matrix(data[!, t_inds]))').^2
σ_z = zeros(size(z))
for i in 1:size(z, 2)
    σ_z[:, i] = sqrt.(diag(Z * Diagonal(var_β[:, 1]) * Z'))
end

## adding indicators

# CHANGE: activity indicators and thier errors go here
# you can actually have as many as you want, but obviously it will take longer
# to fit
obs_indicator1 = z[1, :]
obs_indicator1_err = σ_z[1, :]
obs_indicator2 = z[2, :]
obs_indicator2_err = σ_z[2, :]

# removing means as the GP model assumes zero mean
GLOM_RV.remove_mean!(obs_rvs)
GLOM_RV.remove_mean!(obs_indicator1)
GLOM_RV.remove_mean!(obs_indicator2)

# CHANGE: change these lines if you add more than 2 indicators
# this takes the data and riffles it together so it takes the form
# [rv_1, ind1_1, ind2_1, rv_2, ind1_2, ind2_2, ...]
n_out = 3  # number of indicators + 1
obs_ys = collect(Iterators.flatten(zip(obs_rvs, obs_indicator1, obs_indicator2)))
obs_noise = collect(Iterators.flatten(zip(obs_rvs_err, obs_indicator1_err, obs_indicator2_err)))

# How many differention orders we want in the GLOM model
n_dif = 2 + 1

# CHANGE: consider changing a (the GLOM coefficients that are used, see
# commented lines below)
# If all a's active:
glo = GLOM.GLO(kernel_function, num_kernel_hyperparameters, n_dif, n_out, obs_xs, copy(obs_ys); noise=copy(obs_noise), a=[[1. 1 1];[1 1 1];[1 1 1]])
# glo = GLOM.GLO(kernel_function, num_kernel_hyperparameters, n_dif, n_out, obs_xs, copy(obs_ys); noise=copy(obs_noise), a=[[1. 1 0];[1 0 1];[1 0 1]])

# Makes the std of each output equal to 1, improves fitting stability
# the normalizations are stored in glo.normals
GLOM.normalize_GLO!(glo)

# CHANGE: Setting initial fit values
initial_total_hyperparameters = collect(Iterators.flatten(glo.a))
initial_hypers = [[star_rot_rate], [star_rot_rate], [star_rot_rate], [star_rot_rate, 2 * star_rot_rate, 1], [star_rot_rate, 2 * star_rot_rate, 1], [star_rot_rate, 2 * star_rot_rate, 1]]
append!(initial_total_hyperparameters, initial_hypers[kernel_choice])

# For RV, b2, and b3
# initial_total_hyperparameters[:] = [0.0023761507506676236, -0.3054296895488315, -0.11593928780248865, 0.8755168285072893, 11.8214828482746, -2.0796862114288186, 27.485528857812508, -30.289234793914247, -5.34066485785403, 16.17994314368127]
## Fitting GLOM Model

# CHANGE: Setting kernel hyperparameter priors and kick function
# kick functions help avoid saddle points
tighten_lengthscale_priors = 5
if kernel_name in ["pp", "se", "m52"]
    kernel_hyper_priors(hps::Vector{<:Real}, d::Integer) =
        GLOM_RV.kernel_hyper_priors_1λ(hps, d, star_rot_rate, star_rot_rate / 2 / tighten_lengthscale_priors)
    add_kick!(hps::Vector{<:Real}) = GLOM_RV.add_kick_1λ!(hps)
elseif kernel_name == "qp"
    kernel_hyper_priors(hps::Vector{<:Real}, d::Integer) =
        GLOM_RV.kernel_hyper_priors_qp(hps, d, [star_rot_rate, 2 * star_rot_rate, 1], [star_rot_rate / 2, star_rot_rate / 2, 0.4] ./ tighten_lengthscale_priors)
    add_kick!(hps::Vector{<:Real}) = GLOM_RV.add_kick_qp!(hps)
elseif kernel_name in ["se_se", "m52_m52"]
    kernel_hyper_priors(hps::Vector{<:Real}, d::Integer) =
        GLOM_RV.kernel_hyper_priors_2λ(hps, d, [star_rot_rate, 2 * star_rot_rate, 1], [star_rot_rate / 2, star_rot_rate / 2, 1] ./ tighten_lengthscale_priors)
    add_kick!(hps::Vector{<:Real}) = GLOM_RV.add_kick_2λ!(hps)
else
    # kernel_hyper_priors(hps::Vector{<:Real}, d::Integer) = custom function
end

fit1_total_hyperparameters, result = GLOM_RV.fit_GLOM(glo, initial_total_hyperparameters, kernel_hyper_priors, add_kick!)
# fit_GLOM returns a vector of num_kernel_hyperparameters gp hyperparameters
# followed by the GLOM coefficients and the Optim result object
workspace = GLOM.nlogL_matrix_workspace(glo, fit1_total_hyperparameters)

## Plotting initial results

plot_xs = collect(LinRange(obs_xs[1]-10, obs_xs[end]+10, 1000))
rv_std = std(obs_rvs)
noise_std = norm(holder[1:n_out:end])/sqrt(length(obs_rvs))

function posts(glo, hyperparameters)
    post, post_err, post_obs = GLOM_RV.GLOM_posteriors(glo, plot_xs, hyperparameters)
    GLOM_rvs_at_plot_xs, GLOM_ind1_at_plot_xs, GLOM_ind2_at_plot_xs = post
    GLOM_rvs_err_at_plot_xs, GLOM_ind1_err_at_plot_xs, GLOM_ind2_err_at_plot_xs = post_err
    GLOM_rvs_at_obs_xs, GLOM_ind1_at_obs_xs, GLOM_ind2_at_obs_xs = post_obs
    return GLOM_rvs_at_plot_xs, GLOM_ind1_at_plot_xs, GLOM_ind2_at_plot_xs,
        GLOM_rvs_err_at_plot_xs, GLOM_ind1_err_at_plot_xs, GLOM_ind2_err_at_plot_xs,
        GLOM_rvs_at_obs_xs, GLOM_ind1_at_obs_xs, GLOM_ind2_at_obs_xs
end
function std_print()
    println("\nstarting std:           ", rv_std)
    println("PPCA SAFE feat new std: ", after_model_std)
    println("std from noise ~        ", noise_std)
    println("std from activity: ", sqrt(rv_std^2 - noise_std^2), " -> ", sqrt(maximum([0, after_model_std^2 - noise_std^2])))
end
using Plots
function plots(name)
    plt = scatter(obs_xs, obs_rvs, yerror=obs_rvs_err; label="rv obs")
    plot!(plt, plot_xs, GLOM_rvs_at_plot_xs, ribbons=GLOM_rvs_err_at_plot_xs, fillalpha=0.3; label="model")
    savefig("examples/figs/" * name * "_0")
    plt = scatter(obs_xs, obs_indicator1, yerror=obs_indicator1_err; label="ind1 obs")
    plot!(plt, plot_xs, GLOM_ind1_at_plot_xs, ribbons=GLOM_ind1_err_at_plot_xs, fillalpha=0.3; label="model")
    savefig("examples/figs/" * name * "_1")
    plt = scatter(obs_xs, obs_indicator2, yerror=obs_indicator2_err; label="ind2 obs")
    plot!(plt, plot_xs, GLOM_ind2_at_plot_xs, ribbons=GLOM_ind2_err_at_plot_xs, fillalpha=0.3; label="model")
    savefig("examples/figs/" * name * "_2")
end

holder = copy(glo.noise)
glo.noise[1:n_out:end] .*= 100
GLOM_rvs_at_plot_xs, GLOM_ind1_at_plot_xs, GLOM_ind2_at_plot_xs,
    GLOM_rvs_err_at_plot_xs, GLOM_ind1_err_at_plot_xs, GLOM_ind2_err_at_plot_xs,
    GLOM_rvs_at_obs_xs, GLOM_ind1_at_obs_xs, GLOM_ind2_at_obs_xs =
    posts(glo, fit1_total_hyperparameters)
after_model_std = std(GLOM_rvs_at_obs_xs - obs_rvs)
std_print()
plots("comp_$(n_bs)_norv_snr600")
glo.noise[:] = holder

GLOM_rvs_at_plot_xs, GLOM_ind1_at_plot_xs, GLOM_ind2_at_plot_xs,
    GLOM_rvs_err_at_plot_xs, GLOM_ind1_err_at_plot_xs, GLOM_ind2_err_at_plot_xs,
    GLOM_rvs_at_obs_xs, GLOM_ind1_at_obs_xs, GLOM_ind2_at_obs_xs =
    posts(glo, fit1_total_hyperparameters)
after_model_std = std(GLOM_rvs_at_obs_xs - obs_rvs)
std_print()
plots("comp_$(n_bs)_norm_snr600")

## Coefficient exploration

#=
# Could use this to iterate through all of the possible combinations of GLOM
# coefficients
possible_as = Matrix[]
GLOM_RV.valid_as!(possible_as, zeros(n_out, n_dif))
append!(possible_as, [ones(n_out, n_dif)])
reverse!(possible_as)  # the expensive, filled out possibilities go first
nℓs = Float64[]
all_fit_total_hyperparameters = Vector[]
glos = GLOM.GLO[]
t0 = Libc.time()
for i in 1:length(possible_as)
    a = possible_as[i]
    glo = GLOM.GLO(kernel_function, num_kernel_hyperparameters, n_dif, n_out, obs_xs, copy(obs_ys); noise=copy(obs_noise), a=a)
    GLOM.normalize_GLO!(glo)
    initial_total_hyperparameters = collect(Iterators.flatten(glo.a))
    append!(initial_total_hyperparameters, initial_hypers[kernel_choice])
    fit_total_hyperparameters, result = GLOM_RV.fit_GLOM(glo, initial_total_hyperparameters, kernel_hyper_priors, add_kick!; print_stuff=false)
    append!(glos, [glo])
    append!(nℓs, result.minimum)
    append!(all_fit_total_hyperparameters, [fit_total_hyperparameters])
    println("\nDone with $(round(100 * i / length(possible_as); digits=2))% of a possibilities")
    t = Libc.time() - t0
    println("t: $(Int(round(t)))s t_left?: $(round(t * (length(possible_as) - i) / i / 60; digits=1)) mins")
end

best_fits = sortperm(nℓs)

plot_xs = collect(LinRange(obs_xs[1]-10, obs_xs[end]+10, 300))
post, post_err, post_obs = GLOM_RV.GLOM_posteriors(glos[best_fits[1]], plot_xs, all_fit_total_hyperparameters[best_fits[1]])
GLOM_rvs_at_plot_xs, GLOM_ind1_at_plot_xs, GLOM_ind2_at_plot_xs = post
GLOM_rvs_err_at_plot_xs, GLOM_ind1_err_at_plot_xs, GLOM_ind2_err_at_plot_xs = post_err
GLOM_rvs_at_obs_xs, GLOM_ind1_at_obs_xs, GLOM_ind2_at_obs_xs = post_obs

using Plots
plt = scatter(obs_xs, obs_rvs, yerror=obs_rvs_err)
plot!(plt, plot_xs, GLOM_rvs_at_plot_xs, ribbons=GLOM_rvs_err_at_plot_xs, fillalpha=0.3)

plt = scatter(obs_xs, obs_indicator1, yerror=obs_indicator1_err)
plot!(plt, plot_xs, GLOM_ind1_at_plot_xs, ribbons=GLOM_ind1_err_at_plot_xs, fillalpha=0.3)

plt = scatter(obs_xs, obs_indicator2, yerror=obs_indicator2_err)
plot!(plt, plot_xs, GLOM_ind2_at_plot_xs, ribbons=GLOM_ind2_err_at_plot_xs, fillalpha=0.3)
=#

## Finding a planet?

@assert inject_planet
nlogprior_hyperparameters(total_hyper::Vector, d::Int) = GLOM_RV.nlogprior_hyperparameters(kernel_hyper_priors, glo.n_kern_hyper, total_hyper, d)
glo_rv = GLO_RV(glo, 1u"d", glo.normals[1]u"m/s")

############
# Post fit #
############

println("starting hyperparameters")
println(initial_total_hyperparameters)
initial_nlogL = GLOM.nlogL_GLOM(glo, initial_total_hyperparameters)
initial_uE = -initial_nlogL - nlogprior_hyperparameters(initial_total_hyperparameters, 0)
println(initial_uE, "\n")

println("ending hyperparameters")
println(fit1_total_hyperparameters)
fit_nlogL1 = GLOM.nlogL_GLOM(glo, fit1_total_hyperparameters)
uE1 = -fit_nlogL1 - nlogprior_hyperparameters(fit1_total_hyperparameters, 0)
println(uE1, "\n")

#########################
# Keplerian periodogram #
#########################

# sample linearly in frequency space so that we get periods from the 1 / uneven Nyquist
freq_grid = GLOM_RV.autofrequency(glo_rv.time; samples_per_peak=11)
period_grid = 1 ./ reverse(freq_grid)
amount_of_periods = length(period_grid)

Σ_obs = GLOM.Σ_observations(glo, fit1_total_hyperparameters)

# making necessary variables local to all workers
fit1_total_hyperparameters_nz = GLOM.remove_zeros(fit1_total_hyperparameters)
nlogprior_kernel = nlogprior_hyperparameters(fit1_total_hyperparameters_nz, 0)

using Distributed

# CHANGE: can change full_fit to false to just do the epicyclic fit which is
# ~40x faster. Distributing the compute only gives me a factor of 2, so only
# try if the keplerian fitting takes a long time
use_distributed = true
full_fit = true

# concurrency is weird so you may have to run this twice
if use_distributed
    GLOM.auto_addprocs()
    @everywhere using Pkg; @everywhere Pkg.activate("examples"); @everywhere Pkg.instantiate()
    @everywhere using GLOM_RV_Example; @everywhere GLOM_RV = GLOM_RV_Example
    @everywhere import GPLinearODEMaker; @everywhere GLOM = GPLinearODEMaker
    @everywhere using UnitfulAstro, Unitful

    GLOM.sendto(workers(), kernel_name=kernel_name)
    @everywhere GLOM.include_kernel(kernel_name)

    GLOM.sendto(workers(), glo_rv=glo_rv, fit1_total_hyperparameters=fit1_total_hyperparameters, Σ_obs=Σ_obs, nlogprior_kernel=nlogprior_kernel, full_fit=full_fit)
end

@everywhere function fit_kep_hold_P(P::Unitful.Time; fast::Bool=false, kwargs...)
    #initialize with fast epicyclic fit
    ks = GLOM_RV.fit_kepler(glo_rv, Σ_obs, GLOM_RV.kep_signal_epicyclic(P=P))
    if !fast
        ks = GLOM_RV.fit_kepler(glo_rv, Σ_obs, GLOM_RV.kep_signal_wright(0u"m/s", P, ks.M0, minimum([ks.e, 0.3]), 0, 0u"m/s"); hold_P=true, avoid_saddle=false, print_stuff=false, kwargs...)
        return ks
    end
    if ks == nothing
        ks = GLOM_RV.fit_kepler(glo_rv, Σ_obs, GLOM_RV.kep_signal_wright(0u"m/s", P, 2 * π * rand(), 0.1, 0, 0u"m/s"); hold_P=true, avoid_saddle=false, print_stuff=false, kwargs...)
        return ks
    end
    return ks
end
@everywhere function kep_unnormalized_posterior_distributed(P::Unitful.Time; kwargs...)
    ks = fit_kep_hold_P(P; kwargs...)
    if ks == nothing
        return [-Inf, -Inf]
    else
        val = GLOM.nlogL_GLOM(
            glo_rv.GLO,
            fit1_total_hyperparameters;
            Σ_obs=Σ_obs,
            y_obs=GLOM_RV.remove_kepler(glo_rv, ks))
        return [-val, GLOM_RV.logprior_kepler(ks; use_hk=false) - nlogprior_kernel - val]
    end
end
@everywhere kep_unnormalized_posterior_distributed(P::Unitful.Time) = kep_unnormalized_posterior_distributed(P; fast=!full_fit)

likelihoods = zeros(amount_of_periods)
unnorm_posteriors = zeros(amount_of_periods)

@time if use_distributed
    # takes around 115s for 1/2 of SAFE data with 6 workers and 5000 periods
    holder = pmap(x->kep_unnormalized_posterior_distributed(x), period_grid, batch_size=Int(floor(amount_of_periods / (nworkers() + 1)) + 1))
    likelihoods[:] = [holder[i][1] for i in 1:length(holder)]
    unnorm_posteriors[:] = [holder[i][2] for i in 1:length(holder)]
else
    # takes around 15s for 101501 data and 1200 periods
    for i in 1:amount_of_periods
        likelihoods[i], unnorm_posteriors[i] = kep_unnormalized_posterior_distributed(period_grid[i])
    end
end

best_periods = period_grid[GLOM_RV.find_modes(unnorm_posteriors; amount=10)]
best_period = best_periods[1]

println("found period:    $(ustrip(best_period)) days")

# using Plots
# plot(ustrip.(period_grid), likelihoods; xaxis=:log, leg=false)
# plot(ustrip.(period_grid), unnorm_posteriors; xaxis=:log, leg=false)

####################################################################################################
# Refitting GP with full planet signal at found period subtracted (K,ω,γ-linear, P,M0,e-nonlinear) #
####################################################################################################
remainder(vec, x) = [i > 0 ? i % x : (i % x) + x for i in vec]

current_ks = fit_kep_hold_P(best_period; print_stuff=true)
println("before wright fit: ", GLOM_RV.kep_parms_str(current_ks))

#=
plot_kep_xs = collect(LinRange(0, ustrip(best_period), 1000))
# scatter(remainder(glo.x_obs, ustrip(best_period)), ustrip.(glo_rv.rv); yerror=ustrip.(glo_rv.rv_noise), label="data")
scatter(remainder(glo.x_obs, ustrip(best_period)), clean_rvs; yerror=GLOM_rvs_err_at_obs_xs, label="\"clean\" data")
plot!(plot_kep_xs, ustrip.(current_ks.(plot_kep_xs.*u"d")); label="kep")
=#
fit2_total_hyperparameters, current_ks = GLOM_RV.fit_GLOM_and_kep!(workspace,
    glo_rv, fit1_total_hyperparameters, kernel_hyper_priors,
    add_kick!, current_ks; avoid_saddle=false)

# # these should be near 0
# GLOM_RV.test_∇nlogL_kep(glo_rv, workspace.Σ_obs, current_ks; include_priors=true)
###########################################################################################
# Refitting GP with full planet signal at found period subtracted (K,P,M0,e,ω,γ-nonlinear)#
###########################################################################################

current_ks = GLOM_RV.kep_signal(current_ks)
println("\nbefore full fit: ", GLOM_RV.kep_parms_str(current_ks))
fit3_total_hyperparameters, current_ks = GLOM_RV.fit_GLOM_and_kep!(workspace,
    glo_rv, fit2_total_hyperparameters, kernel_hyper_priors,
    add_kick!, current_ks; avoid_saddle=true)

full_ks = GLOM_RV.kep_signal(current_ks)

# # these should be near 0
# GLOM_RV.test_∇nlogL_kep(glo_rv, workspace.Σ_obs, full_ks; include_priors=true)

###################
# Post planet fit #
###################

println("fit hyperparameters")
println(fit1_total_hyperparameters)
println(uE1, "\n")

println("kepler hyperparameters")
println(fit3_total_hyperparameters)
fit_nlogL2 = GLOM.nlogL_GLOM!(workspace, glo, fit3_total_hyperparameters; y_obs=GLOM_RV.remove_kepler(glo_rv, full_ks))
uE2 = -fit_nlogL2 - nlogprior_hyperparameters(fit3_total_hyperparameters, 0) + GLOM_RV.logprior_kepler(full_ks; use_hk=true)
println(uE2, "\n")

println("best fit keplerian")
println(GLOM_RV.kep_parms_str(full_ks))

##################################################################################
# refitting noise model to see if a better model was found during planet fitting #
##################################################################################

fit1_total_hyperparameters_temp, result = GLOM_RV.fit_GLOM(
    glo,
    fit3_total_hyperparameters,
    kernel_hyper_priors,
    add_kick!)

println(result)

println("first fit hyperparameters")
println(fit1_total_hyperparameters)
println(uE1, "\n")

println("fit after planet hyperparameters")
println(fit1_total_hyperparameters_temp)
fit_nlogL1_temp = GLOM.nlogL_GLOM!(workspace, glo, fit1_total_hyperparameters_temp)
uE1_temp = -fit_nlogL1_temp - nlogprior_hyperparameters(fit1_total_hyperparameters_temp, 0)
println(uE1_temp, "\n")

if uE1_temp > uE1
    println("new fit is better, switching hps")
    fit1_total_hyperparameters[:] = fit1_total_hyperparameters_temp
    fit_nlogL1 = fit_nlogL1_temp
    uE1 = uE1_temp
end

##########################
# Evidence approximation #
##########################

# no planet
H1 = (GLOM.∇∇nlogL_GLOM(glo, fit1_total_hyperparameters)
    + nlogprior_hyperparameters(GLOM.remove_zeros(fit1_total_hyperparameters), 2))
try
    global E1 = GLOM.log_laplace_approximation(H1, -uE1, 0)
catch err
    if isa(err, DomainError)
        println("Laplace approximation failed for initial GP fit")
        println("det(H1): $(det(H1)) (should've been positive)")
        global E1 = 0
    else
        rethrow()
    end
end

# planet
H2 = Matrix(GLOM_RV.∇∇nlogL_GLOM_and_planet!(workspace, glo_rv, fit3_total_hyperparameters, full_ks; include_kepler_priors=true))
n_hyper = length(GLOM.remove_zeros(fit3_total_hyperparameters))
H2[1:n_hyper, 1:n_hyper] += nlogprior_hyperparameters(GLOM.remove_zeros(fit3_total_hyperparameters), 2)
try
    global E2 = GLOM.log_laplace_approximation(Symmetric(H2), -uE2, 0)
catch err
    if isa(err, DomainError)
        println("Laplace approximation failed for planet fit")
        println("det(H2): $(det(H2)) (should've been positive)")
        global E2 = 0
    else
        rethrow()
    end
end

println("\nlog likelihood for GLOM model: " * string(-fit_nlogL1))
println("log likelihood for GLOM + planet model: " * string(-fit_nlogL2))

println("\nunnormalized posterior for GLOM model: " * string(uE1))
println("unnormalized posterior for GLOM + planet model: " * string(uE2))

println("\nevidence for GLOM model: " * string(E1))
println("evidence for GLOM + planet model: " * string(E2))

# # This should be pretty close to true
# GLOM_RV.test_∇∇nlogL_kep(glo_rv, workspace.Σ_obs, current_ks; include_priors=true)

## Plotting final results

post, post_err, post_obs = GLOM_RV.GLOM_posteriors(glo, plot_xs, fit3_total_hyperparameters; y_obs=GLOM_RV.remove_kepler(glo_rv, full_ks))
GLOM_rvs_at_plot_xs, GLOM_ind1_at_plot_xs, GLOM_ind2_at_plot_xs = post
GLOM_rvs_err_at_plot_xs, GLOM_ind1_err_at_plot_xs, GLOM_ind2_err_at_plot_xs = post_err
GLOM_rvs_at_obs_xs, GLOM_ind1_at_obs_xs, GLOM_ind2_at_obs_xs = post_obs

using Plots
plt = scatter(obs_xs, GLOM_RV.remove_kepler(glo_rv, full_ks)[1:3:end], yerror=obs_rvs_err)
plot!(plt, plot_xs, GLOM_rvs_at_plot_xs, ribbons=GLOM_rvs_err_at_plot_xs, fillalpha=0.3)

plt = scatter(obs_xs, obs_indicator1, yerror=obs_indicator1_err)
plot!(plt, plot_xs, GLOM_ind1_at_plot_xs, ribbons=GLOM_ind1_err_at_plot_xs, fillalpha=0.3)

plt = scatter(obs_xs, obs_indicator2, yerror=obs_indicator2_err)
plot!(plt, plot_xs, GLOM_ind2_at_plot_xs, ribbons=GLOM_ind2_err_at_plot_xs, fillalpha=0.3)

println(fit3_total_hyperparameters)


text:
#
# Copyright (c) 2021 Tobias Thummerer, Lars Mikelsons, Josef Kircher
# Licensed under the MIT license. See LICENSE file in the project root for details.
#

pathToFMU = joinpath(dirname(@__FILE__), "..", "model", ENV["EXPORTINGTOOL"], "SpringPendulum1D.fmu")

myFMU = fmiLoad(pathToFMU)

comp = fmiInstantiate!(myFMU; loggingOn=false)
@test comp != 0

# choose FMU or FMUComponent
fmuStruct = nothing
envFMUSTRUCT = ENV["FMUSTRUCT"]
if envFMUSTRUCT == "FMU"
    fmuStruct = myFMU
elseif envFMUSTRUCT == "FMUCOMPONENT"
    fmuStruct = comp
end
@assert fmuStruct != nothing "Unknwon fmuStruct, environment variable `FMUSTRUCT` = `$envFMUSTRUCT`"

t_start = 0.0
t_stop = 8.0

data = fmiSimulateCS(fmuStruct, t_start, t_stop; recordValues=["mass.s", "mass.v"])
@test length(data.dataPoints) == 100
@test length(data.dataPoints[1]) == 3

t = collect(d[1] for d in data.dataPoints)
s = collect(d[2] for d in data.dataPoints)
v = collect(d[3] for d in data.dataPoints)
@test t[1] == t_start
@test t[end] == t_stop

# reference values from Simulation in Dymola2020x (Dassl)
@test s[1] == 0.5
@test v[1] == 0.0
#@test abs(s[end] - 0.509219) < 0.01
#@test abs(v[end] - 0.314074) < 0.01

data = fmiSimulateCS(fmuStruct, t_start, t_stop)
@test data == nothing   # nothing was recorded, recordValues=[]

fmiUnload(myFMU)

ext: jl


text:
#==============================================================================#
# KinesisVideoArchivedMedia.jl
#
# This file is generated from:
# https://github.com/aws/aws-sdk-js/blob/master/apis/kinesis-video-archived-media-2017-09-30.normal.json
#==============================================================================#

__precompile__()

module KinesisVideoArchivedMedia

using AWSCore

"""
    using AWSSDK.KinesisVideoArchivedMedia.get_hlsstreaming_session_url
    get_hlsstreaming_session_url([::AWSConfig], arguments::Dict)
    get_hlsstreaming_session_url([::AWSConfig]; <keyword arguments>)

    using AWSCore.Services.kinesis_video_archived_media
    kinesis_video_archived_media([::AWSConfig], "POST", "/getHLSStreamingSessionURL", arguments::Dict)
    kinesis_video_archived_media([::AWSConfig], "POST", "/getHLSStreamingSessionURL", <keyword arguments>)

# GetHLSStreamingSessionURL Operation

Retrieves an HTTP Live Streaming (HLS) URL for the stream. You can then open the URL in a browser or media player to view the stream contents.

You must specify either the `StreamName` or the `StreamARN`.

An Amazon Kinesis video stream has the following requirements for providing data through HLS:

*   The media type must be `video/h264`.

*   Data retention must be greater than 0.

*   The fragments must contain codec private data in the AVC (Advanced Video Coding) for H.264 format ([MPEG-4 specification ISO/IEC 14496-15](https://www.iso.org/standard/55980.html)). For information about adapting stream data to a given format, see [NAL Adaptation Flags](http://docs.aws.amazon.com/kinesisvideostreams/latest/dg/latest/dg/producer-reference-nal.html).

Kinesis Video Streams HLS sessions contain fragments in the fragmented MPEG-4 form (also called fMP4 or CMAF), rather than the MPEG-2 form (also called TS chunks, which the HLS specification also supports). For more information about HLS fragment types, see the [HLS specification](https://tools.ietf.org/html/draft-pantos-http-live-streaming-23).

The following procedure shows how to use HLS with Kinesis Video Streams:

1.  Get an endpoint using [GetDataEndpoint](http://docs.aws.amazon.com/kinesisvideostreams/latest/dg/API_GetDataEndpoint.html), specifying `GET_HLS_STREAMING_SESSION_URL` for the `APIName` parameter.

2.  Retrieve the HLS URL using `GetHLSStreamingSessionURL`. Kinesis Video Streams creates an HLS streaming session to be used for accessing content in a stream using the HLS protocol. `GetHLSStreamingSessionURL` returns an authenticated URL (that includes an encrypted session token) for the session's HLS *master playlist* (the root resource needed for streaming with HLS).

    **Note**
    > Don't share or store this token where an unauthorized entity could access it. The token provides access to the content of the stream. Safeguard the token with the same measures that you would use with your AWS credentials.

    The media that is made available through the playlist consists only of the requested stream, time range, and format. No other media data (such as frames outside the requested window or alternate bit rates) is made available.

3.  Provide the URL (containing the encrypted session token) for the HLS master playlist to a media player that supports the HLS protocol. Kinesis Video Streams makes the HLS media playlist, initialization fragment, and media fragments available through the master playlist URL. The initialization fragment contains the codec private data for the stream, and other data needed to set up the video decoder and renderer. The media fragments contain H.264-encoded video frames and time stamps.

4.  The media player receives the authenticated URL and requests stream metadata and media data normally. When the media player requests data, it calls the following actions:

    *   **GetHLSMasterPlaylist:** Retrieves an HLS master playlist, which contains a URL for the `GetHLSMediaPlaylist` action, and additional metadata for the media player, including estimated bit rate and resolution.

    *   **GetHLSMediaPlaylist:** Retrieves an HLS media playlist, which contains a URL to access the MP4 initialization fragment with the `GetMP4InitFragment` action, and URLs to access the MP4 media fragments with the `GetMP4MediaFragment` actions. The HLS media playlist also contains metadata about the stream that the player needs to play it, such as whether the `PlaybackMode` is `LIVE` or `ON_DEMAND`. The HLS media playlist is typically static for sessions with a `PlaybackType` of `ON_DEMAND`. The HLS media playlist is continually updated with new fragments for sessions with a `PlaybackType` of `LIVE`.

    *   **GetMP4InitFragment:** Retrieves the MP4 initialization fragment. The media player typically loads the initialization fragment before loading any media fragments. This fragment contains the "`fytp`" and "`moov`" MP4 atoms, and the child atoms that are needed to initialize the media player decoder.

        The initialization fragment does not correspond to a fragment in a Kinesis video stream. It contains only the codec private data for the stream, which the media player needs to decode video frames.

    *   **GetMP4MediaFragment:** Retrieves MP4 media fragments. These fragments contain the "`moof`" and "`mdat`" MP4 atoms and their child atoms, containing the encoded fragment's video frames and their time stamps.

        **Note**
        > After the first media fragment is made available in a streaming session, any fragments that don't contain the same codec private data are excluded in the HLS media playlist. Therefore, the codec private data does not change between fragments in a session.

        Data retrieved with this action is billable. See [Pricing](aws.amazon.comkinesis/video-streams/pricing/) for details.

**Note**
> The following restrictions apply to HLS sessions:

*   A streaming session URL should not be shared between players. The service might throttle a session if multiple media players are sharing it. For connection limits, see [Kinesis Video Streams Limits](http://docs.aws.amazon.com/kinesisvideostreams/latest/dg/limits.html).

*   A Kinesis video stream can have a maximum of five active HLS streaming sessions. If a new session is created when the maximum number of sessions is already active, the oldest (earliest created) session is closed. The number of active `GetMedia` connections on a Kinesis video stream does not count against this limit, and the number of active HLS sessions does not count against the active `GetMedia` connection limit.

You can monitor the amount of data that the media player consumes by monitoring the `GetMP4MediaFragment.OutgoingBytes` Amazon CloudWatch metric. For information about using CloudWatch to monitor Kinesis Video Streams, see [Monitoring Kinesis Video Streams](http://docs.aws.amazon.com/kinesisvideostreams/latest/dg/monitoring.html). For pricing information, see [Amazon Kinesis Video Streams Pricing](https://aws.amazon.com/kinesis/video-streams/pricing/) and [AWS Pricing](https://aws.amazon.com/pricing/). Charges for both HLS sessions and outgoing AWS data apply.

For more information about HLS, see [HTTP Live Streaming](https://developer.apple.com/streaming/) on the [Apple Developer site](https://developer.apple.com).

# Arguments

## `StreamName = ::String`
The name of the stream for which to retrieve the HLS master playlist URL.

You must specify either the `StreamName` or the `StreamARN`.

## `StreamARN = ::String`
The Amazon Resource Name (ARN) of the stream for which to retrieve the HLS master playlist URL.

You must specify either the `StreamName` or the `StreamARN`.

## `PlaybackMode = "LIVE" or "ON_DEMAND"`
Whether to retrieve live or archived, on-demand data.

Features of the two types of session include the following:

*   **`LIVE`** : For sessions of this type, the HLS media playlist is continually updated with the latest fragments as they become available. We recommend that the media player retrieve a new playlist on a one-second interval. When this type of session is played in a media player, the user interface typically displays a "live" notification, with no scrubber control for choosing the position in the playback window to display.

    **Note**
    > In `LIVE` mode, the newest available fragments are included in an HLS media playlist, even if there is a gap between fragments (that is, if a fragment is missing). A gap like this might cause a media player to halt or cause a jump in playback. In this mode, fragments are not added to the HLS media playlist if they are older than the newest fragment in the playlist. If the missing fragment becomes available after a subsequent fragment is added to the playlist, the older fragment is not added, and the gap is not filled.

*   **`ON_DEMAND`** : For sessions of this type, the HLS media playlist contains all the fragments for the session, up to the number that is specified in `MaxMediaPlaylistFragmentResults`. The playlist must be retrieved only once for each session. When this type of session is played in a media player, the user interface typically displays a scrubber control for choosing the position in the playback window to display.

In both playback modes, if `FragmentSelectorType` is `PRODUCER_TIMESTAMP`, and if there are multiple fragments with the same start time stamp, the fragment that has the larger fragment number (that is, the newer fragment) is included in the HLS media playlist. The other fragments are not included. Fragments that have different time stamps but have overlapping durations are still included in the HLS media playlist. This can lead to unexpected behavior in the media player.

The default is `LIVE`.

## `HLSFragmentSelector = [ ... ]`
The time range of the requested fragment, and the source of the time stamps.

This parameter is required if `PlaybackMode` is `ON_DEMAND`. This parameter is optional if `PlaybackMode` is `LIVE`. If `PlaybackMode` is `LIVE`, the `FragmentSelectorType` can be set, but the `TimestampRange` should not be set. If `PlaybackMode` is `ON_DEMAND`, both `FragmentSelectorType` and `TimestampRange` must be set.
```
 HLSFragmentSelector = [
        "FragmentSelectorType" =>  "PRODUCER_TIMESTAMP" or "SERVER_TIMESTAMP",
        "TimestampRange" =>  [
            "StartTimestamp" =>  timestamp,
            "EndTimestamp" =>  timestamp
        ]
    ]
```

## `DiscontinuityMode = "ALWAYS" or "NEVER"`
Specifies when flags marking discontinuities between fragments will be added to the media playlists. The default is `ALWAYS` when [HLSFragmentSelector](@ref) is `SERVER_TIMESTAMP`, and `NEVER` when it is `PRODUCER_TIMESTAMP`.

Media players typically build a timeline of media content to play, based on the time stamps of each fragment. This means that if there is any overlap between fragments (as is typical if [HLSFragmentSelector](@ref) is `SERVER_TIMESTAMP`), the media player timeline has small gaps between fragments in some places, and overwrites frames in other places. When there are discontinuity flags between fragments, the media player is expected to reset the timeline, resulting in the fragment being played immediately after the previous fragment. We recommend that you always have discontinuity flags between fragments if the fragment time stamps are not accurate or if fragments might be missing. You should not place discontinuity flags between fragments for the player timeline to accurately map to the producer time stamps.

## `Expires = ::Int`
The time in seconds until the requested session expires. This value can be between 300 (5 minutes) and 43200 (12 hours).

When a session expires, no new calls to `GetHLSMasterPlaylist`, `GetHLSMediaPlaylist`, `GetMP4InitFragment`, or `GetMP4MediaFragment` can be made for that session.

The default is 300 (5 minutes).

## `MaxMediaPlaylistFragmentResults = ::Int`
The maximum number of fragments that are returned in the HLS media playlists.

When the `PlaybackMode` is `LIVE`, the most recent fragments are returned up to this value. When the `PlaybackMode` is `ON_DEMAND`, the oldest fragments are returned, up to this maximum number.

When there are a higher number of fragments available in a live HLS media playlist, video players often buffer content before starting playback. Increasing the buffer size increases the playback latency, but it decreases the likelihood that rebuffering will occur during playback. We recommend that a live HLS media playlist have a minimum of 3 fragments and a maximum of 10 fragments.

The default is 5 fragments if `PlaybackMode` is `LIVE`, and 1,000 if `PlaybackMode` is `ON_DEMAND`.

The maximum value of 1,000 fragments corresponds to more than 16 minutes of video on streams with 1-second fragments, and more than 2 1/2 hours of video on streams with 10-second fragments.

# Returns

`GetHLSStreamingSessionURLOutput`

# Exceptions

`ResourceNotFoundException`, `InvalidArgumentException`, `ClientLimitExceededException`, `NotAuthorizedException`, `UnsupportedStreamMediaTypeException`, `NoDataRetentionException`, `MissingCodecPrivateDataException` or `InvalidCodecPrivateDataException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/kinesis-video-archived-media-2017-09-30/GetHLSStreamingSessionURL)
"""
@inline get_hlsstreaming_session_url(aws::AWSConfig=default_aws_config(); args...) = get_hlsstreaming_session_url(aws, args)

@inline get_hlsstreaming_session_url(aws::AWSConfig, args) = AWSCore.Services.kinesis_video_archived_media(aws, "POST", "/getHLSStreamingSessionURL", args)

@inline get_hlsstreaming_session_url(args) = get_hlsstreaming_session_url(default_aws_config(), args)

"""
    using AWSSDK.KinesisVideoArchivedMedia.get_media_for_fragment_list
    get_media_for_fragment_list([::AWSConfig], arguments::Dict)
    get_media_for_fragment_list([::AWSConfig]; StreamName=, Fragments=)

    using AWSCore.Services.kinesis_video_archived_media
    kinesis_video_archived_media([::AWSConfig], "POST", "/getMediaForFragmentList", arguments::Dict)
    kinesis_video_archived_media([::AWSConfig], "POST", "/getMediaForFragmentList", StreamName=, Fragments=)

# GetMediaForFragmentList Operation

Gets media for a list of fragments (specified by fragment number) from the archived data in an Amazon Kinesis video stream.

The following limits apply when using the `GetMediaForFragmentList` API:

*   A client can call `GetMediaForFragmentList` up to five times per second per stream.

*   Kinesis Video Streams sends media data at a rate of up to 25 megabytes per second (or 200 megabits per second) during a `GetMediaForFragmentList` session.

# Arguments

## `StreamName = ::String` -- *Required*
The name of the stream from which to retrieve fragment media.

## `Fragments = [::String, ...]` -- *Required*
A list of the numbers of fragments for which to retrieve media. You retrieve these values with [ListFragments](@ref).

# Returns

`GetMediaForFragmentListOutput`

# Exceptions

`ResourceNotFoundException`, `InvalidArgumentException`, `ClientLimitExceededException` or `NotAuthorizedException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/kinesis-video-archived-media-2017-09-30/GetMediaForFragmentList)
"""
@inline get_media_for_fragment_list(aws::AWSConfig=default_aws_config(); args...) = get_media_for_fragment_list(aws, args)

@inline get_media_for_fragment_list(aws::AWSConfig, args) = AWSCore.Services.kinesis_video_archived_media(aws, "POST", "/getMediaForFragmentList", args)

@inline get_media_for_fragment_list(args) = get_media_for_fragment_list(default_aws_config(), args)

"""
    using AWSSDK.KinesisVideoArchivedMedia.list_fragments
    list_fragments([::AWSConfig], arguments::Dict)
    list_fragments([::AWSConfig]; StreamName=, <keyword arguments>)

    using AWSCore.Services.kinesis_video_archived_media
    kinesis_video_archived_media([::AWSConfig], "POST", "/listFragments", arguments::Dict)
    kinesis_video_archived_media([::AWSConfig], "POST", "/listFragments", StreamName=, <keyword arguments>)

# ListFragments Operation

Returns a list of [Fragment](@ref) objects from the specified stream and start location within the archived data.

# Arguments

## `StreamName = ::String` -- *Required*
The name of the stream from which to retrieve a fragment list.

## `MaxResults = ::Int`
The total number of fragments to return. If the total number of fragments available is more than the value specified in `max-results`, then a [ListFragmentsOutput\$NextToken](@ref) is provided in the output that you can use to resume pagination.

## `NextToken = ::String`
A token to specify where to start paginating. This is the [ListFragmentsOutput\$NextToken](@ref) from a previously truncated response.

## `FragmentSelector = [ ... ]`
Describes the time stamp range and time stamp origin for the range of fragments to return.
```
 FragmentSelector = [
        "FragmentSelectorType" => <required> "PRODUCER_TIMESTAMP" or "SERVER_TIMESTAMP",
        "TimestampRange" => <required> [
            "StartTimestamp" => <required> timestamp,
            "EndTimestamp" => <required> timestamp
        ]
    ]
```

# Returns

`ListFragmentsOutput`

# Exceptions

`ResourceNotFoundException`, `InvalidArgumentException`, `ClientLimitExceededException` or `NotAuthorizedException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/kinesis-video-archived-media-2017-09-30/ListFragments)
"""
@inline list_fragments(aws::AWSConfig=default_aws_config(); args...) = list_fragments(aws, args)

@inline list_fragments(aws::AWSConfig, args) = AWSCore.Services.kinesis_video_archived_media(aws, "POST", "/listFragments", args)

@inline list_fragments(args) = list_fragments(default_aws_config(), args)

end # module KinesisVideoArchivedMedia

#==============================================================================#
# End of file
#==============================================================================#


text:
# grep "#define" cpxconst.h | grep "CPX_PARAM_|CPXPARAM" | awk '{ printf("\"%s\" => %s,\n",$2,$3) }'
const paramName2Indx = Dict(
"CPX_AUTO_SOLN" => 0,
"CPX_BASIC_SOLN" => 1,
"CPX_NONBASIC_SOLN" => 2,
"CPX_PRIMAL_SOLN" => 3,
"CPX_PARAM_ALL_MIN" => 1000,
"CPX_PARAM_ALL_MAX" => 6000,
"CPXPARAM_H" => 1,
"CPXPARAM_Advance" => 1001,
"CPXPARAM_Barrier_Algorithm" => 3007,
"CPXPARAM_Barrier_ColNonzeros" => 3009,
"CPXPARAM_Barrier_ConvergeTol" => 3002,
"CPXPARAM_Barrier_Crossover" => 3018,
"CPXPARAM_Barrier_Display" => 3010,
"CPXPARAM_Barrier_Limits_Corrections" => 3013,
"CPXPARAM_Barrier_Limits_Growth" => 3003,
"CPXPARAM_Barrier_Limits_Iteration" => 3012,
"CPXPARAM_Barrier_Limits_ObjRange" => 3004,
"CPXPARAM_Barrier_Ordering" => 3014,
"CPXPARAM_Barrier_QCPConvergeTol" => 3020,
"CPXPARAM_Barrier_StartAlg" => 3017,
"CPXPARAM_Benders_Strategy" => 1501,
"CPXPARAM_Benders_Tolerances_feasibilitycut" => 1509,
"CPXPARAM_Benders_Tolerances_optimalitycut" => 1510,
"CPXPARAM_Benders_WorkerAlgorithm" => 1500,
"CPXPARAM_ClockType" => 1006,
"CPXPARAM_Conflict_Algorithm" => 1073,
"CPXPARAM_Conflict_Display" => 1074,
"CPXPARAM_CPUmask" => 1144,
"CPXPARAM_DetTimeLimit" => 1127,
"CPXPARAM_DistMIP_Rampup_DetTimeLimit" => 2164,
"CPXPARAM_DistMIP_Rampup_Duration" => 2163,
"CPXPARAM_DistMIP_Rampup_TimeLimit" => 2165,
"CPXPARAM_Emphasis_Memory" => 1082,
"CPXPARAM_Emphasis_MIP" => 2058,
"CPXPARAM_Emphasis_Numerical" => 1083,
"CPXPARAM_Feasopt_Mode" => 1084,
"CPXPARAM_Feasopt_Tolerance" => 2073,
"CPXPARAM_LPMethod" => 1062,
"CPXPARAM_MIP_Cuts_BQP" => 2195,
"CPXPARAM_MIP_Cuts_Cliques" => 2003,
"CPXPARAM_MIP_Cuts_Covers" => 2005,
"CPXPARAM_MIP_Cuts_Disjunctive" => 2053,
"CPXPARAM_MIP_Cuts_FlowCovers" => 2040,
"CPXPARAM_MIP_Cuts_Gomory" => 2049,
"CPXPARAM_MIP_Cuts_GUBCovers" => 2044,
"CPXPARAM_MIP_Cuts_Implied" => 2041,
"CPXPARAM_MIP_Cuts_LiftProj" => 2152,
"CPXPARAM_MIP_Cuts_LocalImplied" => 2181,
"CPXPARAM_MIP_Cuts_MCFCut" => 2134,
"CPXPARAM_MIP_Cuts_MIRCut" => 2052,
"CPXPARAM_MIP_Cuts_PathCut" => 2051,
"CPXPARAM_MIP_Cuts_RLT" => 2196,
"CPXPARAM_MIP_Cuts_ZeroHalfCut" => 2111,
"CPXPARAM_MIP_Display" => 2012,
"CPXPARAM_MIP_Interval" => 2013,
"CPXPARAM_MIP_Limits_AggForCut" => 2054,
"CPXPARAM_MIP_Limits_AuxRootThreads" => 2139,
"CPXPARAM_MIP_Limits_CutPasses" => 2056,
"CPXPARAM_MIP_Limits_CutsFactor" => 2033,
"CPXPARAM_MIP_Limits_EachCutLimit" => 2102,
"CPXPARAM_MIP_Limits_GomoryCand" => 2048,
"CPXPARAM_MIP_Limits_GomoryPass" => 2050,
"CPXPARAM_MIP_Limits_Nodes" => 2017,
"CPXPARAM_MIP_Limits_PolishTime" => 2066,
"CPXPARAM_MIP_Limits_Populate" => 2108,
"CPXPARAM_MIP_Limits_ProbeDetTime" => 2150,
"CPXPARAM_MIP_Limits_ProbeTime" => 2065,
"CPXPARAM_MIP_Limits_RepairTries" => 2067,
"CPXPARAM_MIP_Limits_Solutions" => 2015,
"CPXPARAM_MIP_Limits_StrongCand" => 2045,
"CPXPARAM_MIP_Limits_StrongIt" => 2046,
"CPXPARAM_MIP_Limits_TreeMemory" => 2027,
"CPXPARAM_MIP_OrderType" => 2032,
"CPXPARAM_MIP_PolishAfter_AbsMIPGap" => 2126,
"CPXPARAM_MIP_PolishAfter_DetTime" => 2151,
"CPXPARAM_MIP_PolishAfter_MIPGap" => 2127,
"CPXPARAM_MIP_PolishAfter_Nodes" => 2128,
"CPXPARAM_MIP_PolishAfter_Solutions" => 2129,
"CPXPARAM_MIP_PolishAfter_Time" => 2130,
"CPXPARAM_MIP_Pool_AbsGap" => 2106,
"CPXPARAM_MIP_Pool_Capacity" => 2103,
"CPXPARAM_MIP_Pool_Intensity" => 2107,
"CPXPARAM_MIP_Pool_RelGap" => 2105,
"CPXPARAM_MIP_Pool_Replace" => 2104,
"CPXPARAM_MIP_Strategy_Backtrack" => 2002,
"CPXPARAM_MIP_Strategy_BBInterval" => 2039,
"CPXPARAM_MIP_Strategy_Branch" => 2001,
"CPXPARAM_MIP_Strategy_CallbackReducedLP" => 2055,
"CPXPARAM_MIP_Strategy_Dive" => 2060,
"CPXPARAM_MIP_Strategy_File" => 2016,
"CPXPARAM_MIP_Strategy_FPHeur" => 2098,
"CPXPARAM_MIP_Strategy_HeuristicEffort" => 2120,
"CPXPARAM_MIP_Strategy_HeuristicFreq" => 2031,
"CPXPARAM_MIP_Strategy_KappaStats" => 2137,
"CPXPARAM_MIP_Strategy_LBHeur" => 2063,
"CPXPARAM_MIP_Strategy_MIQCPStrat" => 2110,
"CPXPARAM_MIP_Strategy_NodeSelect" => 2018,
"CPXPARAM_MIP_Strategy_Order" => 2020,
"CPXPARAM_MIP_Strategy_PresolveNode" => 2037,
"CPXPARAM_MIP_Strategy_Probe" => 2042,
"CPXPARAM_MIP_Strategy_RINSHeur" => 2061,
"CPXPARAM_MIP_Strategy_Search" => 2109,
"CPXPARAM_MIP_Strategy_StartAlgorithm" => 2025,
"CPXPARAM_MIP_Strategy_SubAlgorithm" => 2026,
"CPXPARAM_MIP_Strategy_VariableSelect" => 2028,
"CPXPARAM_MIP_SubMIP_StartAlg" => 2205,
"CPXPARAM_MIP_SubMIP_SubAlg" => 2206,
"CPXPARAM_MIP_SubMIP_NodeLimit" => 2212,
"CPXPARAM_MIP_SubMIP_Scale" => 2207,
"CPXPARAM_MIP_Tolerances_AbsMIPGap" => 2008,
"CPXPARAM_MIP_Tolerances_Linearization" => 2068,
"CPXPARAM_MIP_Tolerances_Integrality" => 2010,
"CPXPARAM_MIP_Tolerances_LowerCutoff" => 2006,
"CPXPARAM_MIP_Tolerances_MIPGap" => 2009,
"CPXPARAM_MIP_Tolerances_ObjDifference" => 2019,
"CPXPARAM_MIP_Tolerances_RelObjDifference" => 2022,
"CPXPARAM_MIP_Tolerances_UpperCutoff" => 2007,
"CPXPARAM_MultiObjective_Display" => 1600,
"CPXPARAM_Network_Display" => 5005,
"CPXPARAM_Network_Iterations" => 5001,
"CPXPARAM_Network_NetFind" => 1022,
"CPXPARAM_Network_Pricing" => 5004,
"CPXPARAM_Network_Tolerances_Feasibility" => 5003,
"CPXPARAM_Network_Tolerances_Optimality" => 5002,
"CPXPARAM_OptimalityTarget" => 1131,
"CPXPARAM_Output_CloneLog" => 1132,
"CPXPARAM_Output_IntSolFilePrefix" => 2143,
"CPXPARAM_Output_MPSLong" => 1081,
"CPXPARAM_Output_WriteLevel" => 1114,
"CPXPARAM_Parallel" => 1109,
"CPXPARAM_ParamDisplay" => 1163,
"CPXPARAM_Preprocessing_Aggregator" => 1003,
"CPXPARAM_Preprocessing_BoundStrength" => 2029,
"CPXPARAM_Preprocessing_CoeffReduce" => 2004,
"CPXPARAM_Preprocessing_Dependency" => 1008,
"CPXPARAM_Preprocessing_Dual" => 1044,
"CPXPARAM_Preprocessing_Fill" => 1002,
"CPXPARAM_Preprocessing_Folding" => 1164,
"CPXPARAM_Preprocessing_Linear" => 1058,
"CPXPARAM_Preprocessing_NumPass" => 1052,
"CPXPARAM_Preprocessing_Presolve" => 1030,
"CPXPARAM_Preprocessing_QCPDuals" => 4003,
"CPXPARAM_Preprocessing_QPMakePSD" => 4010,
"CPXPARAM_Preprocessing_QToLin" => 4012,
"CPXPARAM_Preprocessing_Reduce" => 1057,
"CPXPARAM_Preprocessing_Relax" => 2034,
"CPXPARAM_Preprocessing_RepeatPresolve" => 2064,
"CPXPARAM_Preprocessing_Symmetry" => 2059,
"CPXPARAM_QPMethod" => 1063,
"CPXPARAM_RandomSeed" => 1124,
"CPXPARAM_Read_APIEncoding" => 1130,
"CPXPARAM_Read_Constraints" => 1021,
"CPXPARAM_Read_DataCheck" => 1056,
"CPXPARAM_Read_FileEncoding" => 1129,
"CPXPARAM_Read_Nonzeros" => 1024,
"CPXPARAM_Read_QPNonzeros" => 4001,
"CPXPARAM_Read_Scale" => 1034,
"CPXPARAM_Read_Variables" => 1023,
"CPXPARAM_Read_WarningLimit" => 1157,
"CPXPARAM_Record" => 1162,
"CPXPARAM_ScreenOutput" => 1035,
"CPXPARAM_Sifting_Algorithm" => 1077,
"CPXPARAM_Sifting_Simplex" => 1158,
"CPXPARAM_Sifting_Display" => 1076,
"CPXPARAM_Sifting_Iterations" => 1078,
"CPXPARAM_Simplex_Crash" => 1007,
"CPXPARAM_Simplex_DGradient" => 1009,
"CPXPARAM_Simplex_Display" => 1019,
"CPXPARAM_Simplex_DynamicRows" => 1161,
"CPXPARAM_Simplex_Limits_Iterations" => 1020,
"CPXPARAM_Simplex_Limits_LowerObj" => 1025,
"CPXPARAM_Simplex_Limits_Perturbation" => 1028,
"CPXPARAM_Simplex_Limits_Singularity" => 1037,
"CPXPARAM_Simplex_Limits_UpperObj" => 1026,
"CPXPARAM_Simplex_Perturbation_Constant" => 1015,
"CPXPARAM_Simplex_Perturbation_Indicator" => 1027,
"CPXPARAM_Simplex_PGradient" => 1029,
"CPXPARAM_Simplex_Pricing" => 1010,
"CPXPARAM_Simplex_Refactor" => 1031,
"CPXPARAM_Simplex_Tolerances_Feasibility" => 1016,
"CPXPARAM_Simplex_Tolerances_Markowitz" => 1013,
"CPXPARAM_Simplex_Tolerances_Optimality" => 1014,
"CPXPARAM_SolutionType" => 1147,
"CPXPARAM_Threads" => 1067,
"CPXPARAM_TimeLimit" => 1039,
"CPXPARAM_Tune_DetTimeLimit" => 1139,
"CPXPARAM_Tune_Display" => 1113,
"CPXPARAM_Tune_Measure" => 1110,
"CPXPARAM_Tune_Repeat" => 1111,
"CPXPARAM_Tune_TimeLimit" => 1112,
"CPXPARAM_WorkDir" => 1064,
"CPXPARAM_WorkMem" => 1065,
"CPX_PARAM_ADVIND" => 1001,
"CPX_PARAM_AGGFILL" => 1002,
"CPX_PARAM_AGGIND" => 1003,
"CPX_PARAM_CLOCKTYPE" => 1006,
"CPX_PARAM_CRAIND" => 1007,
"CPX_PARAM_DEPIND" => 1008,
"CPX_PARAM_DPRIIND" => 1009,
"CPX_PARAM_PRICELIM" => 1010,
"CPX_PARAM_EPMRK" => 1013,
"CPX_PARAM_EPOPT" => 1014,
"CPX_PARAM_EPPER" => 1015,
"CPX_PARAM_EPRHS" => 1016,
"CPX_PARAM_SIMDISPLAY" => 1019,
"CPX_PARAM_ITLIM" => 1020,
"CPX_PARAM_ROWREADLIM" => 1021,
"CPX_PARAM_NETFIND" => 1022,
"CPX_PARAM_COLREADLIM" => 1023,
"CPX_PARAM_NZREADLIM" => 1024,
"CPX_PARAM_OBJLLIM" => 1025,
"CPX_PARAM_OBJULIM" => 1026,
"CPX_PARAM_PERIND" => 1027,
"CPX_PARAM_PERLIM" => 1028,
"CPX_PARAM_PPRIIND" => 1029,
"CPX_PARAM_PREIND" => 1030,
"CPX_PARAM_REINV" => 1031,
"CPX_PARAM_SCAIND" => 1034,
"CPX_PARAM_SCRIND" => 1035,
"CPX_PARAM_SINGLIM" => 1037,
"CPX_PARAM_TILIM" => 1039,
"CPX_PARAM_PREDUAL" => 1044,
"CPX_PARAM_PREPASS" => 1052,
"CPX_PARAM_DATACHECK" => 1056,
"CPX_PARAM_REDUCE" => 1057,
"CPX_PARAM_PRELINEAR" => 1058,
"CPX_PARAM_LPMETHOD" => 1062,
"CPX_PARAM_QPMETHOD" => 1063,
"CPX_PARAM_WORKDIR" => 1064,
"CPX_PARAM_WORKMEM" => 1065,
"CPX_PARAM_THREADS" => 1067,
"CPX_PARAM_CONFLICTALG" => 1073,
"CPX_PARAM_CONFLICTDISPLAY" => 1074,
"CPX_PARAM_SIFTDISPLAY" => 1076,
"CPX_PARAM_SIFTALG" => 1077,
"CPX_PARAM_SIFTITLIM" => 1078,
"CPX_PARAM_MPSLONGNUM" => 1081,
"CPX_PARAM_MEMORYEMPHASIS" => 1082,
"CPX_PARAM_NUMERICALEMPHASIS" => 1083,
"CPX_PARAM_FEASOPTMODE" => 1084,
"CPX_PARAM_PARALLELMODE" => 1109,
"CPX_PARAM_TUNINGMEASURE" => 1110,
"CPX_PARAM_TUNINGREPEAT" => 1111,
"CPX_PARAM_TUNINGTILIM" => 1112,
"CPX_PARAM_TUNINGDISPLAY" => 1113,
"CPX_PARAM_WRITELEVEL" => 1114,
"CPX_PARAM_RANDOMSEED" => 1124,
"CPX_PARAM_DETTILIM" => 1127,
"CPX_PARAM_FILEENCODING" => 1129,
"CPX_PARAM_APIENCODING" => 1130,
"CPX_PARAM_OPTIMALITYTARGET" => 1131,
"CPX_PARAM_CLONELOG" => 1132,
"CPX_PARAM_TUNINGDETTILIM" => 1139,
"CPX_PARAM_CPUMASK" => 1144,
"CPX_PARAM_SOLUTIONTYPE" => 1147,
"CPX_PARAM_WARNLIM" => 1157,
"CPX_PARAM_SIFTSIM" => 1158,
"CPX_PARAM_DYNAMICROWS" => 1161,
"CPX_PARAM_RECORD" => 1162,
"CPX_PARAM_PARAMDISPLAY" => 1163,
"CPX_PARAM_FOLDING" => 1164,
"CPX_PARAM_WORKERALG" => 1500,
"CPX_PARAM_BENDERSSTRATEGY" => 1501,
"CPX_PARAM_BENDERSFEASCUTTOL" => 1509,
"CPX_PARAM_BENDERSOPTCUTTOL" => 1510,
"CPX_PARAM_MULTIOBJDISPLAY" => 1600,
"CPX_PARAM_BRDIR" => 2001,
"CPX_PARAM_BTTOL" => 2002,
"CPX_PARAM_CLIQUES" => 2003,
"CPX_PARAM_COEREDIND" => 2004,
"CPX_PARAM_COVERS" => 2005,
"CPX_PARAM_CUTLO" => 2006,
"CPX_PARAM_CUTUP" => 2007,
"CPX_PARAM_EPAGAP" => 2008,
"CPX_PARAM_EPGAP" => 2009,
"CPX_PARAM_EPINT" => 2010,
"CPX_PARAM_MIPDISPLAY" => 2012,
"CPX_PARAM_MIPINTERVAL" => 2013,
"CPX_PARAM_INTSOLLIM" => 2015,
"CPX_PARAM_NODEFILEIND" => 2016,
"CPX_PARAM_NODELIM" => 2017,
"CPX_PARAM_NODESEL" => 2018,
"CPX_PARAM_OBJDIF" => 2019,
"CPX_PARAM_MIPORDIND" => 2020,
"CPX_PARAM_RELOBJDIF" => 2022,
"CPX_PARAM_STARTALG" => 2025,
"CPX_PARAM_SUBALG" => 2026,
"CPX_PARAM_TRELIM" => 2027,
"CPX_PARAM_VARSEL" => 2028,
"CPX_PARAM_BNDSTRENIND" => 2029,
"CPX_PARAM_HEURFREQ" => 2031,
"CPX_PARAM_MIPORDTYPE" => 2032,
"CPX_PARAM_CUTSFACTOR" => 2033,
"CPX_PARAM_RELAXPREIND" => 2034,
"CPX_PARAM_PRESLVND" => 2037,
"CPX_PARAM_BBINTERVAL" => 2039,
"CPX_PARAM_FLOWCOVERS" => 2040,
"CPX_PARAM_IMPLBD" => 2041,
"CPX_PARAM_PROBE" => 2042,
"CPX_PARAM_GUBCOVERS" => 2044,
"CPX_PARAM_STRONGCANDLIM" => 2045,
"CPX_PARAM_STRONGITLIM" => 2046,
"CPX_PARAM_FRACCAND" => 2048,
"CPX_PARAM_FRACCUTS" => 2049,
"CPX_PARAM_FRACPASS" => 2050,
"CPX_PARAM_FLOWPATHS" => 2051,
"CPX_PARAM_MIRCUTS" => 2052,
"CPX_PARAM_DISJCUTS" => 2053,
"CPX_PARAM_AGGCUTLIM" => 2054,
"CPX_PARAM_MIPCBREDLP" => 2055,
"CPX_PARAM_CUTPASS" => 2056,
"CPX_PARAM_MIPEMPHASIS" => 2058,
"CPX_PARAM_SYMMETRY" => 2059,
"CPX_PARAM_DIVETYPE" => 2060,
"CPX_PARAM_RINSHEUR" => 2061,
"CPX_PARAM_LBHEUR" => 2063,
"CPX_PARAM_REPEATPRESOLVE" => 2064,
"CPX_PARAM_PROBETIME" => 2065,
"CPX_PARAM_POLISHTIME" => 2066,
"CPX_PARAM_REPAIRTRIES" => 2067,
"CPX_PARAM_EPLIN" => 2068,
"CPX_PARAM_EPRELAX" => 2073,
"CPX_PARAM_FPHEUR" => 2098,
"CPX_PARAM_EACHCUTLIM" => 2102,
"CPX_PARAM_SOLNPOOLCAPACITY" => 2103,
"CPX_PARAM_SOLNPOOLREPLACE" => 2104,
"CPX_PARAM_SOLNPOOLGAP" => 2105,
"CPX_PARAM_SOLNPOOLAGAP" => 2106,
"CPX_PARAM_SOLNPOOLINTENSITY" => 2107,
"CPX_PARAM_POPULATELIM" => 2108,
"CPX_PARAM_MIPSEARCH" => 2109,
"CPX_PARAM_MIQCPSTRAT" => 2110,
"CPX_PARAM_ZEROHALFCUTS" => 2111,
"CPX_PARAM_POLISHAFTEREPAGAP" => 2126,
"CPX_PARAM_POLISHAFTEREPGAP" => 2127,
"CPX_PARAM_POLISHAFTERNODE" => 2128,
"CPX_PARAM_POLISHAFTERINTSOL" => 2129,
"CPX_PARAM_POLISHAFTERTIME" => 2130,
"CPX_PARAM_MCFCUTS" => 2134,
"CPX_PARAM_MIPKAPPASTATS" => 2137,
"CPX_PARAM_AUXROOTTHREADS" => 2139,
"CPX_PARAM_INTSOLFILEPREFIX" => 2143,
"CPX_PARAM_PROBEDETTIME" => 2150,
"CPX_PARAM_POLISHAFTERDETTIME" => 2151,
"CPX_PARAM_LANDPCUTS" => 2152,
"CPX_PARAM_RAMPUPDURATION" => 2163,
"CPX_PARAM_RAMPUPDETTILIM" => 2164,
"CPX_PARAM_RAMPUPTILIM" => 2165,
"CPX_PARAM_LOCALIMPLBD" => 2181,
"CPX_PARAM_BQPCUTS" => 2195,
"CPX_PARAM_RLTCUTS" => 2196,
"CPX_PARAM_SUBMIPSTARTALG" => 2205,
"CPX_PARAM_SUBMIPSUBALG" => 2206,
"CPX_PARAM_SUBMIPSCAIND" => 2207,
"CPX_PARAM_SUBMIPNODELIMIT" => 2212,
"CPX_PARAM_BAREPCOMP" => 3002,
"CPX_PARAM_BARGROWTH" => 3003,
"CPX_PARAM_BAROBJRNG" => 3004,
"CPX_PARAM_BARALG" => 3007,
"CPX_PARAM_BARCOLNZ" => 3009,
"CPX_PARAM_BARDISPLAY" => 3010,
"CPX_PARAM_BARITLIM" => 3012,
"CPX_PARAM_BARMAXCOR" => 3013,
"CPX_PARAM_BARORDER" => 3014,
"CPX_PARAM_BARSTARTALG" => 3017,
"CPX_PARAM_BARCROSSALG" => 3018,
"CPX_PARAM_BARQCPEPCOMP" => 3020,
"CPX_PARAM_QPNZREADLIM" => 4001,
"CPX_PARAM_CALCQCPDUALS" => 4003,
"CPX_PARAM_QPMAKEPSDIND" => 4010,
"CPX_PARAM_QTOLININD" => 4012,
"CPX_PARAM_NETITLIM" => 5001,
"CPX_PARAM_NETEPOPT" => 5002,
"CPX_PARAM_NETEPRHS" => 5003,
"CPX_PARAM_NETPPRIIND" => 5004,
"CPX_PARAM_NETDISPLAY" => 5005,
)


text:
@testset "minor allele dosage" begin
@testset "slow" begin
    path = BGEN.datadir("example.16bits.zstd.bgen")
    b = Bgen(path)
    for v in iterator(b)
        probs = probabilities!(b, v)
        dose = minor_allele_dosage!(b, v)

        # dosages for each allele
        a1 = probs[1, :] .* 2 .+ probs[2, :]
        a2 = probs[3, :] .* 2 .+ probs[2, :]

        dose_correct = sum(a1[.!isnan.(a1)]) < sum(a2[.!isnan.(a2)]) ? a1 : a2
        minor_allele_index = sum(a1[.!isnan.(a1)]) < sum(a2[.!isnan.(a2)]) ? 1 : 2
        @test v.genotypes.minor_idx == minor_allele_index
        @test all(isapprox.(dose, dose_correct; atol=2e-7, nans=true))
    end
end

@testset "fast" begin
    path = BGEN.datadir("example.8bits.bgen")
    b = Bgen(path)
    for v in iterator(b)
        probs = probabilities!(b, v)
        dose = minor_allele_dosage!(b, v)

        # dosages for each allele
        a1 = probs[1, :] .* 2 .+ probs[2, :]
        a2 = probs[3, :] .* 2 .+ probs[2, :]

        dose_correct = sum(a1[.!isnan.(a1)]) < sum(a2[.!isnan.(a2)]) ? a1 : a2
        minor_allele_index = sum(a1[.!isnan.(a1)]) < sum(a2[.!isnan.(a2)]) ? 1 : 2
        @test v.genotypes.minor_idx == minor_allele_index
        @test all(isapprox.(dose, dose_correct; atol=2e-7, nans=true))
    end
end

@testset "v11" begin
    path = BGEN.datadir("example.v11.bgen")
    b = Bgen(path)
    for v in iterator(b)
        probs = probabilities!(b, v)
        dose = minor_allele_dosage!(b, v)

        # dosages for each allele
        a1 = probs[1, :] .* 2 .+ probs[2, :]
        a2 = probs[3, :] .* 2 .+ probs[2, :]

        dose_correct = sum(a1[.!isnan.(a1)]) < sum(a2[.!isnan.(a2)]) ? a1 : a2
        minor_allele_index = sum(a1[.!isnan.(a1)]) < sum(a2[.!isnan.(a2)]) ? 1 : 2
        @test v.genotypes.minor_idx == minor_allele_index
        @test all(isapprox.(dose, dose_correct; atol=7e-5, nans=true))
    end
end

@testset "mean_impute" begin
    path = BGEN.datadir("example.8bits.bgen")
    b = Bgen(path)
    v = first(iterator(b; from_bgen_start=true))
    m = minor_allele_dosage!(b, v; mean_impute=true)
    @test isapprox(m[1], 0.3958112303037447)
end

@testset "haplotypes" begin
    path = BGEN.datadir("haplotypes.bgen")
    b = Bgen(path)
    for v in iterator(b)
        probs = probabilities!(b, v)
        dose = ref_allele_dosage!(b, v)
        dose_correct = probs[1, :] .+ probs[3, :]
        @test all(isapprox.(dose, dose_correct))
    end
end
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct LogAnalyticsOperationResult <: SwaggerModel
    properties::Any # spec type: Union{ Nothing, LogAnalyticsOutput } # spec name: properties

    function LogAnalyticsOperationResult(;properties=nothing)
        o = new()
        validate_property(LogAnalyticsOperationResult, Symbol("properties"), properties)
        setfield!(o, Symbol("properties"), properties)
        o
    end
end # type LogAnalyticsOperationResult

const _property_map_LogAnalyticsOperationResult = Dict{Symbol,Symbol}(Symbol("properties")=>Symbol("properties"))
const _property_types_LogAnalyticsOperationResult = Dict{Symbol,String}(Symbol("properties")=>"LogAnalyticsOutput")
Base.propertynames(::Type{ LogAnalyticsOperationResult }) = collect(keys(_property_map_LogAnalyticsOperationResult))
Swagger.property_type(::Type{ LogAnalyticsOperationResult }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_LogAnalyticsOperationResult[name]))}
Swagger.field_name(::Type{ LogAnalyticsOperationResult }, property_name::Symbol) =  _property_map_LogAnalyticsOperationResult[property_name]

function check_required(o::LogAnalyticsOperationResult)
    true
end

function validate_property(::Type{ LogAnalyticsOperationResult }, name::Symbol, val)
end


text:
using Test, Distributions, SpecialFunctions
using QuadGK

# `numeric_norm` is a helper function to compute numerically the squared L2
# norms of the distributions.  These methods aren't very robust because can't
# deal with divergent norms, or discrete distributions with infinite support.
numeric_norm(d::ContinuousUnivariateDistribution) =
    quadgk(x -> pdf(d, x) ^ 2, support(d).lb, support(d).ub)[1]

function numeric_norm(d::DiscreteUnivariateDistribution)
    # When the distribution has infinite support, sum up to an arbitrary large
    # value.
    upper = isfinite(maximum(d)) ? round(Int, maximum(d)) : 100
    return sum(pdf(d, k) ^ 2 for k in round(Int, minimum(d)):upper)
end

@testset "pdf L2 norm" begin
    # Test error on a non implemented norm.
    @test_throws MethodError pdfsquaredL2norm(Gumbel())

    @testset "Bernoulli" begin
        for d in (Bernoulli(0.5), Bernoulli(0), Bernoulli(0.25))
            @test pdfsquaredL2norm(d) ≈ numeric_norm(d)
        end
        # The norm is the same for complementary probabilities
        @test pdfsquaredL2norm(Bernoulli(0)) == pdfsquaredL2norm(Bernoulli(1))
        @test pdfsquaredL2norm(Bernoulli(0.25)) == pdfsquaredL2norm(Bernoulli(0.75))
    end

    @testset "Beta" begin
        for d in (Beta(1, 1), Beta(2, 2))
            @test pdfsquaredL2norm(d) ≈ numeric_norm(d)
        end
        @test pdfsquaredL2norm(Beta(0.25, 1)) ≈ Inf
        @test pdfsquaredL2norm(Beta(1, 0.25)) ≈ Inf
    end

    @testset "Categorical" begin
        for n in (1, 2, 5, 10)
            d = Categorical(collect(1 / n for _ in 1:n))
            @test pdfsquaredL2norm(d) ≈ numeric_norm(d)
        end
        for d in (Categorical([0.25, 0.75]), Categorical([1 / 6, 1 / 3, 1 / 2]))
            @test pdfsquaredL2norm(d) ≈ numeric_norm(d)
        end
    end

    @testset "Cauchy" begin
        for d in (Cauchy(0, 1), Cauchy(0, 2))
            @test pdfsquaredL2norm(d) ≈ numeric_norm(d)
        end
        # The norm doesn't depend on the mean
        @test pdfsquaredL2norm(Cauchy(100, 1)) == pdfsquaredL2norm(Cauchy(-100, 1)) == pdfsquaredL2norm(Cauchy(0, 1))
    end

    @testset "Chi" begin
        @test pdfsquaredL2norm(Chi(2)) ≈ numeric_norm(Chi(2))
        @test pdfsquaredL2norm(Chi(0.25)) ≈ Inf
    end

    @testset "Chisq" begin
        @test pdfsquaredL2norm(Chisq(2)) ≈ numeric_norm(Chisq(2))
        @test pdfsquaredL2norm(Chisq(1)) ≈ Inf
    end

    @testset "DiscreteUniform" begin
        for d in (DiscreteUniform(-1, 1), DiscreteUniform(1, 2))
            @test pdfsquaredL2norm(d) ≈ numeric_norm(d)
        end
    end

    @testset "Exponential" begin
        for d in (Exponential(1), Exponential(2))
            @test pdfsquaredL2norm(d) ≈ numeric_norm(d)
        end
    end

    @testset "Gamma" begin
        for d in (Gamma(1, 1), Gamma(1, 2), Gamma(2, 2), Gamma(1, 0.25))
            @test pdfsquaredL2norm(d) ≈ numeric_norm(d)
        end
        @test pdfsquaredL2norm(Gamma(0.5, 1)) ≈ Inf
    end

    @testset "Geometric" begin
        for d in (Geometric(0.20), Geometric(0.25), Geometric(0.50), Geometric(0.75), Geometric(0.80))
            @test pdfsquaredL2norm(d) ≈ numeric_norm(d)
        end
    end

    @testset "Logistic" begin
        for d in (Logistic(0, 1), Logistic(0, 2))
            @test pdfsquaredL2norm(d) ≈ numeric_norm(d)
        end
        # The norm doesn't depend on the mean
        @test pdfsquaredL2norm(Logistic(100, 1)) == pdfsquaredL2norm(Logistic(-100, 1)) == pdfsquaredL2norm(Logistic(0, 1))
    end

    @testset "Normal" begin
        for d in (Normal(0, 1), Normal(0, 2))
            @test pdfsquaredL2norm(d) ≈ numeric_norm(d)
        end
        @test pdfsquaredL2norm(Normal(1, 0)) ≈ Inf
        # The norm doesn't depend on the mean
        @test pdfsquaredL2norm(Normal(100, 1)) == pdfsquaredL2norm(Normal(-100, 1)) == pdfsquaredL2norm(Normal(0, 1))
    end

    @testset "Poisson" begin
        for d in (Poisson(0), Poisson(1), Poisson(pi))
            @test pdfsquaredL2norm(d) ≈ numeric_norm(d)
        end
    end

    @testset "Uniform" begin
        for d in (Uniform(-1, 1), Uniform(1, 2))
            @test pdfsquaredL2norm(d) ≈ numeric_norm(d)
        end
    end
end


text:
module TestLP

import VecMathOptInterface
using Test

const MOI = VecMathOptInterface
const MOIU = MOI.Utilities
const LP = MOI.FileFormats.LP
const LP_TEST_FILE = "test.lp"

function test_show()
    @test sprint(show, LP.Model()) == "A .LP-file model"
end

function test_comprehensive_write()
    model = LP.Model()
    MOIU.loadfromstring!(
        model,
        """
variables: a, x, y, z
minobjective: x
x >= -1.0
x <= 2.0
y == 3.0
z in Interval(4.0, 5.0)
c5: 1.1x + 0.0 <= 5.1
c6: 1.3x + -1.4 >= -0.1
c7: 1.5a + 1.6 == 0.2
c8: 1.7a + 1.8 in Interval(0.3, 0.4)
x in ZeroOne()
y in Integer()
""",
    )
    MOI.write_to_file(model, LP_TEST_FILE)
    @test read(LP_TEST_FILE, String) ==
          "minimize\n" *
          "obj: x\n" *
          "subject to\n" *
          "c5: 1.1 x <= 5.1\n" *
          "c6: -1.4 + 1.3 x >= -0.1\n" *
          "c7: 1.6 + 1.5 a = 0.2\n" *
          "c8: 0.3 <= 1.8 + 1.7 a <= 0.4\n" *
          "Bounds\n" *
          "x <= 2\n" *
          "x >= -1\n" *
          "y = 3\n" *
          "4 <= z <= 5\n" *
          "a free\n" *
          "General\n" *
          "y\n" *
          "Binary\n" *
          "x\n" *
          "End\n"

    @test !MOI.is_empty(model)
    MOI.empty!(model)
    @test MOI.is_empty(model)
end

function test_name_sanitization_start()
    for starting_letter in [".", "0", "E", "e"]
        model = LP.Model()
        MOI.Utilities.loadfromstring!(
            model,
            """
variables: x
minobjective: x
c1: 1.0 * x >= -1.0
""",
        )
        x = MOI.get(model, MOI.VariableIndex, "x")
        MOI.set(model, MOI.VariableName(), x, starting_letter * "x")
        c = MOI.get(model, MOI.ConstraintIndex, "c1")
        MOI.set(model, MOI.ConstraintName(), c, starting_letter * "c1")
        MOI.write_to_file(model, LP_TEST_FILE)
        @test read(LP_TEST_FILE, String) ==
              "minimize\n" *
              "obj: _$(starting_letter)x\n" *
              "subject to\n" *
              "_$(starting_letter)c1: 1 _$(starting_letter)x >= -1\n" *
              "Bounds\n" *
              "_$(starting_letter)x free\n" *
              "End\n"
    end
end

function test_name_sanitization_illegal()
    for illegal_letter in ["[", "]", "*", "^"]
        model = LP.Model()
        MOI.Utilities.loadfromstring!(
            model,
            """
variables: x
minobjective: x
c1: 1.0 * x >= -1.0
""",
        )
        x = MOI.get(model, MOI.VariableIndex, "x")
        MOI.set(model, MOI.VariableName(), x, "x$(illegal_letter)y")
        c = MOI.get(model, MOI.ConstraintIndex, "c1")
        MOI.set(model, MOI.ConstraintName(), c, "c$(illegal_letter)d")
        MOI.write_to_file(model, LP_TEST_FILE)
        @test read(LP_TEST_FILE, String) ==
              "minimize\n" *
              "obj: x_y\n" *
              "subject to\n" *
              "c_d: 1 x_y >= -1\n" *
              "Bounds\n" *
              "x_y free\n" *
              "End\n"
    end
end

function test_name_sanitization_duplicate()
    model = LP.Model()
    MOI.Utilities.loadfromstring!(
        model,
        """
variables: a, b, c, d
minobjective: a + b + c + d
c1: a + b + c + d >= -1.0
""",
    )
    variables = MOI.get(model, MOI.ListOfVariableIndices())
    MOI.set.(model, MOI.VariableName(), variables, ["a[", "a]", "a*", "a^"])
    MOI.write_to_file(model, LP_TEST_FILE)
    @test read(LP_TEST_FILE, String) ==
          "minimize\n" *
          "obj: 1 a_ + 1 a__1 + 1 a__2 + 1 a__3\n" *
          "subject to\n" *
          "c1: 1 a_ + 1 a__1 + 1 a__2 + 1 a__3 >= -1\n" *
          "Bounds\n" *
          "a_ free\n" *
          "a__1 free\n" *
          "a__2 free\n" *
          "a__3 free\n" *
          "End\n"
end

function test_name_sanitization_length()
    model = LP.Model(maximum_length = 3)
    x = MOI.add_variable(model)
    MOI.set(model, MOI.VariableName(), x, "abcdefg")
    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)
    MOI.set(model, MOI.ObjectiveFunction{MOI.VariableIndex}(), x)
    MOI.write_to_file(model, LP_TEST_FILE)
    @test read(LP_TEST_FILE, String) ==
          "minimize\n" *
          "obj: abc\n" *
          "subject to\n" *
          "Bounds\n" *
          "abc free\n" *
          "End\n"
end

function test_name_sanitization_too_long()
    model = LP.Model(maximum_length = 3)
    MOI.Utilities.loadfromstring!(
        model,
        """
variables: abcd, abce
minobjective: abcd + abce
c1: abcd + abce >= -1.0
""",
    )
    MOI.write_to_file(model, LP_TEST_FILE)
    @test read(LP_TEST_FILE, String) ==
          "minimize\n" *
          "obj: 1 abc + 1 abc_1\n" *
          "subject to\n" *
          "c1: 1 abc + 1 abc_1 >= -1\n" *
          "Bounds\n" *
          "abc free\n" *
          "abc_1 free\n" *
          "End\n"
end

function test_name_sanitization_other()
    model = LP.Model()
    MOIU.loadfromstring!(
        model,
        """
variables: x
maxobjective: 2.0 * x + -1.0
""",
    )
    MOI.write_to_file(model, LP_TEST_FILE)
    @test read(LP_TEST_FILE, String) ==
          "maximize\n" *
          "obj: -1 + 2 x\n" *
          "subject to\n" *
          "Bounds\n" *
          "x free\n" *
          "End\n"
end

function test_free_variables()
    model = LP.Model()
    MOIU.loadfromstring!(
        model,
        """
variables: x, y, z
maxobjective: x
x in ZeroOne()
y in Integer()
""",
    )
    MOI.write_to_file(model, LP_TEST_FILE)
    @test read(LP_TEST_FILE, String) ==
          "maximize\n" *
          "obj: x\n" *
          "subject to\n" *
          "Bounds\n" *
          "y free\n" *
          "z free\n" *
          "General\n" *
          "y\n" *
          "Binary\n" *
          "x\n" *
          "End\n"
end

function test_quadratic_objective()
    model = LP.Model()
    @test_throws(
        MOI.UnsupportedAttribute,
        MOIU.loadfromstring!(
            model,
            """
variables: x
minobjective: 1.0*x*x
""",
        )
    )
end

function test_read()
    model = LP.Model()
    exception = ErrorException("read! is not implemented for LP files.")
    @test_throws exception MOI.read_from_file(model, LP_TEST_FILE)
end

function runtests()
    for name in names(@__MODULE__, all = true)
        if startswith("$(name)", "test_")
            @testset "$name" begin
                getfield(@__MODULE__, name)()
            end
        end
    end
    sleep(1.0)  # Allow time for unlink to happen.
    rm(LP_TEST_FILE, force = true)
    return
end

end

TestLP.runtests()


text:
struct OBDD
  ordering::Ordering
  root::BDDNode

  function OBDD(ordering::Array{String,1},bdd::BDDNode)
    try
      O=ListOrdering(ordering)
      return OBDD(O,bdd)
    catch e
      throw(typeof(e)(e.msg))
    end
  end

  function OBDD(O::Ordering,bdd::BDDNode)
    global dynamicOBDDordering

    if !respect(bdd,O)
      if dynamicOBDDordering
        return OBDD(O,string(bdd))
      else
        throw(ArgumentError("$(bdd) does not respect $(O)"))
      end
    end

    try
      return new(O,bdd)
    catch e
      throw(typeof(e)(e.msg))
    end
  end

  function OBDD(O::Ordering,value::BinBoolType)
    try
      return new(O,BDD(value))
    catch e
      throw(typeof(e)(e.msg))
    end
  end
end

function SAT_assignments(T::OBDD)::Array{Dict{String, BinBoolType}}
  return SAT_assignments(T.root)
end

function string(A::OBDD)
  if typeof(A.ordering) == ListOrdering
    ordstr=string(A.ordering)
  else
    ordstr="(?)"
  end
  return string(ordstr,"->",string(A.root))
end

print(io::IO, A::OBDD) = print(io, string(A))
println(io::IO, A::OBDD) = println(io, string(A))

show(io::IO, A::OBDD) = show(io, string(A))

function variables(A::OBDD)
  return variables(A.root)
end

function ~(A::OBDD)
  return OBDD(A.ordering,~A.root)
end

function mergeordering(A::OBDD,B::OBDD)
  if definedfor(A.ordering,variablesin(B.ordering))
    B=changeordering(B,A.ordering)
    return (A,B)
  end

  if definedfor(B.ordering,variablesin(A.ordering))
    A=changeordering(A,B.ordering)
    return (A,B)
  end

  new_ordering=merge(A.ordering,B.ordering)

  A=changeordering(A,new_ordering)
  B=changeordering(B,new_ordering)

  return (A,B)
end

function applyoperator(operator::Function,A::OBDD,B::OBDD)
  global dynamicOBDDordering

  if A.ordering!=B.ordering
    if dynamicOBDDordering
      (A,B)=mergeordering(A,B)
    else
      throw(ArgumentError("$(A) and $(B) do not share the same ordering"))
    end
  end

  return OBDD(A.ordering,applyoperator(operator,A.root,B.root,A.ordering))
end

function applyoperator(operator::Function,a::BinBoolType,B::OBDD)
  return applyoperator(operator,BDD(a),B)
end

function applyoperator(operator::Function,A::OBDD,b::BinBoolType)
  return applyoperator(operator,A,BDD(b))
end

function applyoperator(operator::Function,A::OBDD,B::BDDNode)
  return applyoperator(operator,A,OBDD(A.ordering,B))
end

function applyoperator(operator::Function,A::BDDNode,B::OBDD)
  return applyoperator(operator,OBDD(B.ordering,A),B)
end

function (&)(A::OBDD,B::OBDD)
  return applyoperator(((v1,v2)->v1&v2),A,B)
end

function (&)(a::BinBoolType,B::OBDD)
  return BDD(a)&B
end

function (&)(A::OBDD,b::BinBoolType)
  return A&BDD(b)
end

function (&)(A::OBDD,B::BDDNode)
  return A&OBDD(A.ordering,B)
end

function (&)(A::BDDNode,B::OBDD)
  return OBDD(B.ordering,A)&B
end

function (|)(A::OBDD,B::OBDD)
  return applyoperator(((v1,v2)->v1|v2),A,B)
end

function (|)(a::BinBoolType,B::OBDD)
  return BDD(a)|B
end

function (|)(A::OBDD,b::BinBoolType)
  return A|BDD(b)
end

function (|)(A::OBDD,B::BDDNode)
  return A|OBDD(A.ordering,B)
end

function (|)(A::BDDNode,B::OBDD)
  return OBDD(B.ordering,A)|B
end

function Base.xor(A::OBDD,B::OBDD)
  return applyoperator(((v1,v2)->v1 ⊻ v2),A,B)
end

function Base.xor(a::BinBoolType,B::OBDD)
  return BDD(a) ⊻ B
end

function Base.xor(A::OBDD,b::BinBoolType)
  return A ⊻ BDD(b)
end

function Base.xor(A::OBDD,B::BDDNode)
  return A ⊻ OBDD(A.ordering,B)
end

function Base.xor(A::BDDNode,B::OBDD)
  return OBDD(B.ordering,A) ⊻ B
end

function ==(A::OBDD,B::OBDD)
  global dynamicOBDDordering

  if A.ordering!=B.ordering
    if dynamicOBDDordering
      (A,B)=mergeordering(A,B)
    else
      throw(ArgumentError("$(A) and $(B) do not share the same ordering"))
    end
  end

  return A.root == B.root
end

function ==(A::OBDD,b::BDDNode)
  return A.root == b
end

function ==(a::BDDNode,B::OBDD)
  return a == B.root
end

function ==(A::OBDD,b::BinBoolType)
  return A.root == BDD(b)
end

function ==(a::BinBoolType,B::OBDD)
  return B.root == BDD(a)
end

function ==(A::OBDD,b::String)
  return A == OBDD(A.ordering,b)
end

function ==(a::String,B::OBDD)
  return OBDD(B.ordering,a) == B
end

function restrict(A::OBDD,var::String,value)
  root=restrict(A.root,var,value)
  return OBDD(A.ordering,root)
end

function restrict(A::OBDD, assignment::Dict{String, BinBoolType})
  root=restrict(A.root,assignment)
  return OBDD(A.ordering,root)
end


text:
# ===========
println("Test de l'analyse des logging... ")

include("../src/utils/log.jl")

args = Args.parse_commandline()
Args.set(:level, 2)

# EXEMPLE DE CAPTURE ET RESTAURATION DU FLUS DE SORTIE STANDARD
#
# original_stdout = stdout
# (read_pipe, write_pipe) = redirect_stdout();
#
# redirect_stdout(original_stdout);
# close(write_pipe)

##########################
# CAPTURE DE STDOUT
original_stdout = stdout
(read_pipe, write_pipe) = redirect_stdout();

# --------------------------------------------
lg2("A") # Affiche A sans saut de ligne
lg2("B") # Affiche B sans saut de ligne
ln2("")  # Affiche un saut de ligne

# --------------------------------------------
ln2("C") # Affiche C avec saut de ligne

# --------------------------------------------
# Affiche des valeura de type varié sans saut de ligne malgré le ln3()
# car suffix est positionné pour annuler le "\n" par défaut.
n=5
f=3.1
r=22//7
ln2("n=", n, " f=", f, " r=", r,
    prefix="UN_PREFIX==",
    suffix="==UN_SUFFIX"
   )

# --------------------------------------------
lg2("-1-") # On ajoute à la fin de la ligne précédente
ln2() # n'affiche rien, même pas un saut de ligne
lg2("-2-") # On ajoute à la fin de la ligne précédente
ln2("") # affiche un saut de ligne

# --------------------------------------------
ln2("FIN") # affiche "FIN" puis un saut de ligne

##########################
# RESTAURATION DE STDOUT
redirect_stdout(original_stdout)
close(write_pipe)

# On exploite l'ecriture capturée présendant pour vérification
@testset "Test des logging via capture de stdout par redirect_stdout" begin

    # --------------------------------------------
    # txt = read(read_pipe, String)
    line =  readline(read_pipe, keep=true)
    # @show line
    @test line == "AB\n"

    # --------------------------------------------
    line =  readline(read_pipe, keep=true)
    # @show line
    @test line == "C\n"

    # --------------------------------------------
    line =  readline(read_pipe, keep=true)
    # @show line
    # @test line == "UN_PREFIX==n=5 f=3.1 r=22//7==UN_SUFFIX" # NON car suffix écrasé !
    @test line == "UN_PREFIX==n=5 f=3.1 r=22//7\n"

    # --------------------------------------------
    #
    line =  readline(read_pipe, keep=true)
    # @show line
    @test line == "-1--2-\n"

    # --------------------------------------------
    line =  readline(read_pipe, keep=true)
    # @show line
    @test line == "FIN\n"

end
# Suppress unnecessary output when include this file.
return nothing


text:
# Manipulations for creating the QR code matrix

"""
Finder patterns in the corners of QR codes.
"""
const finderpattern =
    BitArray{2}([1 1 1 1 1 1 1;
                 1 0 0 0 0 0 1;
                 1 0 1 1 1 0 1;
                 1 0 1 1 1 0 1;
                 1 0 1 1 1 0 1;
                 1 0 0 0 0 0 1;
                 1 1 1 1 1 1 1])

"""
Alignment pattern in the center of large QR codes.
"""
const alignmentpattern =
    BitArray{2}([1 1 1 1 1;
                 1 0 0 0 1;
                 1 0 1 0 1;
                 1 0 0 0 1;
                 1 1 1 1 1])

"""
    emptymatrix(version::Int64)

Return a matrix for the QR code filled with data-independent elements. `nothing`
matrix elements act as a placeholder for the data.
"""
function emptymatrix(version::Int64)::Array{Union{Bool,Nothing},2}
    n = (version - 1) * 4 + 21
    # nothing is used as a placeholder for the data
    matrix = Array{Union{Bool,Nothing},2}(nothing, (n, n))

    # Finder Patterns
    matrix[1:7, 1:7]     = finderpattern
    matrix[1:7, n - 6:n] = finderpattern
    matrix[n - 6:n, 1:7] = finderpattern

    # Separators around the finder patterns
    matrix[8, 1:8]     .= 0
    matrix[8, n - 7:n] .= 0
    matrix[n - 7, 1:8] .= 0
    matrix[1:8, 8]     .= 0
    matrix[n - 7:n, 8] .= 0
    matrix[1:8, n - 7] .= 0

    # Alignment Location Patterns
    loc = alignmentlocation[version]
    for i in loc, j in loc
        if !( (i - 1 < 8 && j - 1 < 8)   ||
              (i - 1 < 8 && j + 3 > n - 7) ||
              (i + 3 > n - 7 && j - 1 < 8)    )
            matrix[i - 1:i + 3, j - 1:j + 3] = alignmentpattern
        end
    end

    # Timing Patterns
    timing = repeat(BitArray([0, 1]), (n - 15) ÷ 2)
    matrix[7, 8:n - 8] = timing
    matrix[8:n - 8, 7] = timing

    # Format information area + dark module
    matrix[9, 1:9]     .= 1
    matrix[9, n - 7:n] .= 1
    matrix[1:9, 9]     .= 1
    matrix[n - 7:n, 9] .= 1

    # Version information areas for version 7+
    if version >= 7
        matrix[1:6, n - 10:n - 8] .= 1
        matrix[n - 10:n - 8, 1:6] .= 1
    end

    return matrix
end

"""
    placedata!( matrix::Array{Union{Bool,Nothing},2}, data::BitArray{1})

Fill the matrix with the data.
"""
function placedata!( matrix::Array{Union{Bool,Nothing},2}
                   , data::BitArray{1}
                   )::BitArray{2}
    n = size(matrix, 1)
    col, row = n, n + 1
    while col > 0
        # Skip the column with the timing pattern
        if col == 7
            col -= 1
            continue
        end

        # path goes up and down
        if row > n
            row, δrow = n, -1
        else
            row, δrow = 1, 1
        end

        # place data if matrix element is nothing
        for _ in 1:n
            if isnothing(matrix[row, col])
                matrix[row, col] = popfirst!(data)
            end
            if isnothing(matrix[row, col - 1])
                matrix[row, col - 1] = popfirst!(data)
            end
            row += δrow
        end
        # go left
        col -= 2
    end

    return BitArray{2}(matrix)
end

"""
    makemasks(matrix::Array{Union{Bool,Nothing},2})

Create 8 bitmasks for a given matrix.
"""
function makemasks(matrix::Array{Union{Bool,Nothing},2})::Array{BitArray{2},1}
    n = size(matrix, 1)
    masks = [falses(size(matrix)) for _ in 1:8]

    # Weird indexing due to 0-based indexing in documentation
    for row in 0:n-1, col in 0:n-1
        if !isnothing(matrix[row+1, col+1])
            continue
        end
        if (row + col) % 2 == 0
            masks[1][row+1, col+1] = true
        end
        if row % 2 == 0
            masks[2][row+1, col+1] = true
        end
        if col % 3 == 0
            masks[3][row+1, col+1] = true
        end
        if (row + col) % 3 == 0
            masks[4][row+1, col+1] = true
        end
        if ((row ÷ 2) + (col ÷ 3)) % 2 == 0
            masks[5][row+1, col+1] = true
        end
        if ((row * col) % 2) + ((row * col) % 3) == 0
            masks[6][row+1, col+1] = true
        end
        if (((row * col) % 2) + ((row * col) % 3)) % 2 == 0
            masks[7][row+1, col+1] = true
        end
        if (((row + col) % 2) + ((row * col) % 3)) % 2 == 0
            masks[8][row+1, col+1] = true
        end
    end

    return masks
end

"""
    penalty(matrix::BitArray{2})

Calculate a penalty score in order to pick the best mask.
"""
function penalty(matrix::BitArray{2})::Int64
    rows, cols = size(matrix)

    # Condition 1: 5+ in a row of the same color
    score(c) = foldl(check, c, init = (0, !c[1], 0))[1]
    function check((tot, pb, cnt), b)
        if pb != b
            return (tot, b, 1)
        elseif cnt < 4
            return (tot, b, cnt + 1)
        elseif cnt == 4
            return (tot + 3, b, cnt + 1)
        else
            return (tot + 1, b, cnt)
        end
    end
    if VERSION < v"1.1"
        p1 = sum(score(matrix[i, :]) for i in axes(matrix, 1))
           + sum(score(matrix[:, i]) for i in axes(matrix, 2))
    else
        p1 = sum(map(score, eachrow(matrix)))
           + sum(map(score, eachcol(matrix)))
    end

    # Condition 2: number of 2x2 blocks of the same color
    p2 = 0
    for row in 1:rows - 1, col in 1:cols - 1
        block = matrix[row:row + 1, col:col + 1]
        if all(block) || all(.! block)
            p2 += 2
        end
    end

    # Condition 3: specific patterns in rows or columns
    p3 = 0
    patt1 = BitArray([1, 0, 1, 1, 1, 0, 1, 0, 0, 0 ,0])
    patt2 = BitArray([0, 0, 0 ,0, 1, 0, 1, 1, 1, 0, 1])
    for row in 1:rows, col in 1:cols - 10
        hline = matrix[row, col:col + 10]
        if all(.! xor.(hline, patt1)) || all(.! xor.(hline, patt2))
            p3 += 40
        end

        vline = matrix[col:col + 10, row]
        if all(.! xor.(vline, patt1)) || all(.! xor.(vline, patt2))
            p3 += 40
        end
    end

    # Condition 4: percentage of black and white
    t = (sum(matrix) / length(matrix) * 100) ÷ 5
    p4 = 10 * min(abs(t - 10), abs(t - 9))

    return p1 + p2 + p3 + p4
end

"""
    addformat(matrix::BitArray{2}, mask::Int64, version::Int64, eclevel::ErrCorrLevel)

Add information about the `version` and mask number in `matrix`.
"""
function addformat( matrix::BitArray{2}
                  , mask::Int64
                  , version::Int64
                  , eclevel::ErrCorrLevel)::BitArray{2}

    format = formatinfo[(eclevel, mask)]
    n = size(matrix, 1)

    # Info around the top left finder pattern
    matrix[9, 1:6]    = format[1:6]
    matrix[9, 8:9]    = format[7:8]
    matrix[8, 9]      = format[9]
    matrix[6:-1:1, 9] = format[10:15]

    # Info around the bottom and right finder pattern
    matrix[n:-1:n - 6, 9] = format[1:7]
    matrix[9, n - 7:n]    = format[8:15]

    # Extra version information for version 7+
    if version >= 7
        vinfo = reshape(versioninfo[version], (3, 6))

        matrix[n - 10 : n - 8, 1:6] = vinfo
        matrix[1:6, n - 10 : n - 8] = transpose(vinfo)
    end

    return matrix
end


text:
@testset "Context" begin
    let space = StateSpace([1,2], 2, [0.1, 1.0], 2)
        ctx = Context(0.5, space, [0,1])
        @test ctx.β == 0.5
        @test ctx.system_space === space
        @test ctx.system_state == [0,1]
        @test ctx.bath_space == space
        let a = e^-2.55, b = e^-1.1
            @test ctx.bath_state ≈ [2a/(2a + b), b/(2a + b)]
        end
    end
end

let rng = MersenneTwister(2018)
    @testset "Timestep" begin
        let ctx = Context(0.5, StateSpace([1,2], 4, [0.1, 1.0], 4), [1, 0, 0])
            re = relative_entropy(ctx.system_state, ctx.bath_state)
            for _ in 1:10000
                timestep(rng, ctx)
                r = relative_entropy(ctx.system_state, ctx.bath_state)
                @test isapprox(r, re; atol=1e-9) || r < re
                re = r
            end
            @test ctx.system_state ≈ ctx.bath_state
        end
        let ctx = Context(0.5, StateSpace(1:4, 4, [0.1, 1.0], 4), [1,0,0,0,0])
            re = relative_entropy(ctx.system_state, ctx.bath_state)
            for _ in 1:1000000
                timestep(rng, ctx)
                r = relative_entropy(ctx.system_state, ctx.bath_state)
                @test isapprox(r, re; atol=1e-9) || r < re
                re = r
            end
            @test ctx.system_state ≈ ctx.bath_state
        end
    end
end


text:
module DiscontinuousGherkins

export DiscontinuousGalerkin1DGrid

include("jacobi.jl")
include("vandermonde.jl")
include("lift.jl")
include("grids.jl")

end # module


text:
#Radiative decay pathways for 87Rb
#23/07/2017

function Rb87DecayPaths(ll1,jj1)

    D1trans = 1.56051536 #D1 line transition energy in eV
    D2trans = 1.58999085 #D2 line transition energy in eV
    n4l2j5_2en = 2.40116159
    n4l2j3_2en = 2.40121692
    n6l0j1_2en = 2.49759249
    n6l1j1_2en = 2.94203794
    n6l1j3_2en = 2.95165365
    n5l0j1_2en = 0
    loadparam = [5 1 0.5; 5 1 1.5; 6 1 1.5; 6 1 0.5; 5 0 0.5; 4 2 1.5; 6 0 0.5; 4 2 2.5] #Values to load based on allowed decay pathways
    stateenvec = [D1trans, D2trans, n6l1j3_2en, n6l1j1_2en, n5l0j1_2en, n4l2j3_2en, n6l0j1_2en, n4l2j5_2en] #Energy levels
    Rdecayratevec = zeros(8) #Create a vector to store decay rates in

    if ll1==0||(ll1==2&&jj1==1.5)
        qvec1 = 1
        qvec4 = 4
    elseif ll1==1&&jj1==1.5
        qvec1 = 5
        qvec4 = 8
    elseif ll1==1&&jj1==0.5
        qvec1 = 5
        qvec4 = 7
    elseif ll1==2&&jj1==2.5
        qvec1 = 2
        qvec4 = 3
    end

    return loadparam, stateenvec, Rdecayratevec, qvec1, qvec4
end

max_stars_repo_name: CSChisholm/Rydberg

max_issues_repo_name: CSChisholm/Rydberg

max_forks_repo_name: CSChisholm/Rydberg


text:
include("consts.jl")
include("rules.jl")
include("integral.jl")


text:
using Test
using ColorTypes, FileIO, ImageCore, ImageMagick, Dates, FixedPointNumbers, Statistics

import VideoIO

testdir = dirname(@__FILE__)
videodir = joinpath(testdir, "..", "videos")

VideoIO.TestVideos.available()
VideoIO.TestVideos.download_all()

swapext(f, new_ext) = "$(splitext(f)[1])$new_ext"

@noinline function isblank(img)
    all(c->green(c) == 0, img) || all(c->blue(c) == 0, img) || all(c->red(c) == 0, img) || maximum(rawview(channelview(img))) < 0xcf
end

# tesing that the executables run by testing the standard first line output
@testset "Executable functionality (ffmpeg & ffprobe)" begin
    withenv("PATH" => VideoIO.libpath, "LD_LIBRARY_PATH" => VideoIO.libpath, "DYLD_LIBRARY_PATH" => VideoIO.libpath) do
        out = VideoIO.collectexecoutput(`$(VideoIO.ffmpeg)`)
        @test occursin("ffmpeg version ",out[1])
    end
    withenv("PATH" => VideoIO.libpath, "LD_LIBRARY_PATH" => VideoIO.libpath, "DYLD_LIBRARY_PATH" => VideoIO.libpath) do
        out = VideoIO.collectexecoutput(`$(VideoIO.ffprobe)`)
        @test occursin("ffprobe version ",out[1])
    end
end

@testset "File reading" begin

    for name in VideoIO.TestVideos.names()
        Sys.isapple() && startswith(name, "crescent") && continue
        @testset "Reading $name" begin
            first_frame_file = joinpath(testdir, swapext(name, ".png"))
            fiftieth_frame_file = joinpath(testdir, swapext(name, "")*"50.png")
            first_frame = load(first_frame_file) # comment line when creating png files

            f = VideoIO.testvideo(name)
            v = VideoIO.openvideo(f)

            if size(first_frame, 1) > v.height
                first_frame = first_frame[1+size(first_frame,1)-v.height:end,:]
            end

            img = read(v)

            # Find the first non-trivial image
            while isblank(img)
                read!(v, img)
            end

            #save(first_frame_file,img)        # uncomment line when creating png files)

            @test img == first_frame               # comment line when creating png files
            for i in 1:50
                read!(v,img)
            end
            fiftieth_frame = img
            timebase = v.avin.video_info[1].stream.time_base
            tstamp = v.aVideoFrame[1].pkt_dts
            video_tstamp = v.avin.video_info[1].stream.first_dts
            fiftytime = (tstamp-video_tstamp)/(convert(Float64,timebase.den)/convert(Float64,timebase.num))

            while !eof(v)
                read!(v, img)
            end

            seek(v,float(fiftytime))
            read!(v,img)

            @test img == fiftieth_frame

            # read first frames again, and compare
            seekstart(v)

            read!(v, img)

            while isblank(img)
                read!(v, img)
            end

            @test img == first_frame

            close(v)
        end
    end
end

@testset "IO reading" begin
    for name in VideoIO.TestVideos.names()
        Sys.isapple() && startswith(name, "crescent") && continue
        # TODO: fix me?
        (startswith(name, "ladybird") || startswith(name, "NPS")) && continue
        @testset "Testing $name" begin
            first_frame_file = joinpath(testdir, swapext(name, ".png"))
            first_frame = load(first_frame_file) # comment line when creating png files

            filename = joinpath(videodir, name)
            v = VideoIO.openvideo(open(filename))

            if size(first_frame, 1) > v.height
                first_frame = first_frame[1+size(first_frame,1)-v.height:end,:]
            end

            img = read(v)

            # Find the first non-trivial image
            while isblank(img)
                read!(v, img)
            end

            #save(first_frame_file,img)        # uncomment line when creating png files

            @test img == first_frame               # comment line when creating png files

            while !eof(v)
                read!(v, img)
            end
        end
    end

    VideoIO.testvideo("ladybird") # coverage testing
    @test_throws ErrorException VideoIO.testvideo("rickroll")
    @test_throws ErrorException VideoIO.testvideo("")
end

@testset "Reading video duration and date/datetime" begin
    # tesing the duration and date & time functions:
    file = joinpath(videodir, "annie_oakley.ogg")
    @test VideoIO.get_duration(file) == Dates.Microsecond(24224200)
    @test VideoIO.get_start_time(file) == DateTime(1970, 1, 1)
    @test VideoIO.get_time_duration(file) == (DateTime(1970, 1, 1), Dates.Microsecond(24224200))
end

@testset "Video encoding (read, encode, read, compare)" begin
    file = joinpath(videodir, "annie_oakley.ogg")
    f = VideoIO.openvideo(file)
    imgstack_rgb = []
    imgstack_gray = []
    while !eof(f)
        img = collect(read(f))
        img_gray = convert(Array{Gray{N0f8}},img)
        push!(imgstack_rgb,img)
        push!(imgstack_gray,img_gray)
    end
    @testset "Lossless Grayscale encoding" begin
        file_lossless_gray_copy = joinpath(videodir, "annie_oakley_lossless_gray.mp4")
        prop = [:color_range=>2, :priv_data => ("crf"=>"0","preset"=>"medium")]
        codec_name="libx264"
        VideoIO.encodevideo(file_lossless_gray_copy,imgstack_gray,codec_name=codec_name,AVCodecContextProperties=prop, silent=true)

        fcopy = VideoIO.openvideo(file_lossless_gray_copy,target_format=VideoIO.AV_PIX_FMT_GRAY8)
        imgstack_gray_copy = []
        while !eof(fcopy)
            push!(imgstack_gray_copy,collect(read(fcopy)))
        end
        close(f)
        @test eltype(imgstack_gray) == eltype(imgstack_gray_copy)
        @test length(imgstack_gray) == length(imgstack_gray_copy)
        @test size(imgstack_gray[1]) == size(imgstack_gray_copy[1])
        @test !any(.!(imgstack_gray .== imgstack_gray_copy))
    end

    @testset "Lossless RGB encoding" begin
        file_lossless_rgb_copy = joinpath(videodir, "annie_oakley_lossless_rgb.mp4")
        prop = [:priv_data => ("crf"=>"0","preset"=>"medium")]
        codec_name="libx264rgb"
        VideoIO.encodevideo(file_lossless_rgb_copy,imgstack_rgb,codec_name=codec_name,AVCodecContextProperties=prop, silent=true)

        fcopy = VideoIO.openvideo(file_lossless_rgb_copy)
        imgstack_rgb_copy = []
        while !eof(fcopy)
            img = collect(read(fcopy))
            push!(imgstack_rgb_copy,img)
        end
        close(f)
        @test eltype(imgstack_rgb) == eltype(imgstack_rgb_copy)
        @test length(imgstack_rgb) == length(imgstack_rgb_copy)
        @test size(imgstack_rgb[1]) == size(imgstack_rgb_copy[1])
        @test !any(.!(imgstack_rgb .== imgstack_rgb_copy))
    end
end

#VideoIO.TestVideos.remove_all()


text:
# AST
"""
    get_identifier_table(x::ASTContext)
Return the [`IdentifierTable`](@ref) in the [`ASTContext`](@ref).
"""
get_identifier_table(x::ASTContext) = getIdents(x)

get_name(x::ASTContext, s::String) = get(getIdents(x), s)

function get_builtin_type(ctx::ASTContext, ::Type{T}) where {T<:AbstractBuiltinType}
    @check_ptrs ctx
    return T(ctx)
end

get_decl_type(x::ASTContext, decl) = getTypeDeclType(x, decl)
get_decl_type(x::ASTContext, decl, prev) = getTypeDeclType(x, decl, prev)

get_pointer_type(x::ASTContext, ty::AbstractQualType) = getPointerType(x, ty)
get_lvalue_reference_type(x::ASTContext, ty::AbstractQualType)= getLValueReferenceType(x, ty)

# Decl
get_ast_context(x::AbstractDecl) = getASTContext(x)
get_ast_context(x::DeclContext) = getParentASTContext(x)

get_name(x::AbstractNamedDecl) = getName(x)
get_name(x::DeclarationName) = getAsString(x)

get_begin_loc(x::AbstractDecl) = getBeginLoc(x)
get_end_loc(x::AbstractDecl) = getEndLoc(x)
get_loc(x::AbstractDecl) = getLocation(x)

dump(x::AbstractDecl) = dumpColor(x)

Base.isempty(x::DeclarationName) = isEmpty(x)
Base.string(x::DeclarationName) = get_as_string(x)

dump(x::CXXScopeSpec) = dump(getScopeRep(x))


text:
module Util
    export find_extremes, random_cities

    """
        find_extremes(array)
    Locate minimum and maximum in array.
    #Return value
    ((minimum, idx_of_min), (maximum, idx_of_max))
    """
    function find_extremes(arr::AbstractArray)
        min = max = arr[1]
        min_i = max_i = 1
        for i in 2:length(arr)
            if arr[i] > max
                max, max_i = arr[i], i
            elseif arr[i] < min
                min, min_i = arr[i], i
            end
        end
        return ((min, min_i), (max, max_i))
    end

end


text:
module BenchDescent

using NearestNeighborDescent
using NearestNeighborDescent: local_join!, get_neighbors!
using LightGraphs
using Distances
using BenchmarkTools
using Random

function two_local_join!(graph, data, metric, sample_rate)
    local_join!(graph, data, metric; sample_rate=sample_rate)
    local_join!(graph, data, metric; sample_rate=sample_rate)
end

suite = BenchmarkGroup()

knn_graph_types = [HeapKNNGraph, LockHeapKNNGraph]
num_points = [500]
dims = [50]
num_neighbs = [15]
metric = Euclidean

for graph in knn_graph_types
    for np in num_points
        for d in dims
            Random.seed!(0)
            data = [2 .* rand(d) .- 1 for _ in 1:np] # uniform dist in [-1, 1]
            for k in num_neighbs
                if k ≥ np
                    continue
                end
                for rate in [1, .25]
                    prepath = ["N="*string(np), "D="*string(d), "K="*string(k), string(graph)]
                    postpath = [string(metric), string(rate)]
                    g = graph(data, k, metric())

                    suite[vcat(prepath, ["get_neighbors!"], postpath)] =
                        @benchmarkable get_neighbors!(g_, $rate) setup=(g_ = deepcopy($g))
                    g2 = graph(data, k, metric())
                    suite[vcat(prepath, ["local_join! x 2"], postpath)] =
                        @benchmarkable two_local_join!(g_, $(data), $(metric()), $(rate)) setup=(g_ = deepcopy($g2))
                end
            end
        end
    end
end

end # module

BenchDescent.suite