text:
function _initial_Σrp(mcmc::MCMC{S, M}) where {
                      S <: DiseaseStateSequence,
                      M <: ILM}
  rp = mcmc.risk_priors
  return Float64[i == j ? var(rp[i]) : 0.0 for i in 1:length(rp), j in 1:length(rp)]/(10 * length(rp))
end

function _initial_Σsm(mcmc::MCMC{S, M}) where {
                      S <: DiseaseStateSequence,
                      M <: PhyloILM}
  smp = mcmc.substitution_model_prior
  return Float64[i == j ? var(smp[i]) : 0.0 for i in 1:length(smp), j in 1:length(smp)]/(10 * length(smp))
end

function Base.convert(::Type{Array{Float64, 2}}, sm::Vector{NucleicAcidSubstitutionModel})
  prop_names = [propertynames(sm[1])...]
  return [getproperty(sm[i], θ) for i = 1:length(sm), θ in prop_names]
end

function iterate!(
  mc::MarkovChain{S, M},
  mcmc::MCMC{S, M},
  n::Int64,
  Σrp::Array{Float64, 2},
  σ::Float64;
  condition_on_network::Bool=false,
  event_batches::Int64=1,
  adapt_cov::Int64=100) where {
    S <: DiseaseStateSequence,
    M <: TNILM}
  if adapt_cov < 0
    @warn "`adapt_cov` argument indicates the increment in iterations in which the covariance matrix is updated and must be ≧ 0. Setting to 0 for non-adaptive Metropolis-Hastings."
    adapt_cov = 0
  end
  pmeter = Progress(n, "MCMC progress")
  use_adapted_Σrp = false
  local adapted_Σrp
  if (adapt_cov != 0 & mc.iterations > adapt_cov) && LinearAlgebra.isposdef(value(mc.Σrp))
    use_adapted_Σrp = true
    adapted_Σrp = OnlineStats.value(mc.Σrp) * 2.38^2 / length(mc.risk_parameters[1])
  end
  tr_cache = TransmissionRateCache(individuals(mcmc.population))
  for i = 1:n
    if use_adapted_Σrp
      update!(mc, mcmc, adapted_Σrp, σ, condition_on_network = condition_on_network, event_batches = event_batches, transmission_rate_cache=tr_cache)
    else
      update!(mc, mcmc, Σrp, σ, condition_on_network = condition_on_network, event_batches = event_batches, transmission_rate_cache=tr_cache)
    end
    next!(pmeter)
    @logmsg LogLevel(-5000) "MCMC progress" progress = i/n
    if adapt_cov > 0 && mod(i, adapt_cov) == 0
      OnlineStats.fit!(mc.Σrp, eachrow(convert(Array{Float64, 2}, mc.risk_parameters[end-adapt_cov:end])))
      if LinearAlgebra.isposdef(value(mc.Σrp))
        if !use_adapted_Σrp
          use_adapted_Σrp = true
          @debug "Now using adapted covariance matrix for adaptive MCMC on core $(Distributed.myid()) (updated every $adapt_cov iterations)"
        end
        adapted_Σrp = OnlineStats.value(mc.Σrp) * 2.38^2 / length(mc.risk_parameters[1])
        @debug "Covariance matrix updated for Adaptive Metropolis-Hastings MCMC" Covariance = adapted_Σrp
      end
    end
  end
  return mc
end

function iterate!(
  mc::MarkovChain{S, M},
  mcmc::MCMC{S, M},
  n::Int64,
  Σrp::Array{Float64, 2},
  Σsm::Array{Float64, 2},
  σ::Float64;
  event_batches::Int64=1,
  adapt_cov::Int64=100) where {
    S <: DiseaseStateSequence,
    M <: PhyloILM}
  if adapt_cov < 0
    @warn "`adapt_cov` argument indicates the increment in iterations in which the covariance matrix is updated and must be ≧ 0. Setting to 0 for non-adaptive Metropolis-Hastings."
    adapt_cov = 0
  end
  pmeter = Progress(n, "MCMC progress")
  use_adapted_Σrp = false
  use_adapted_Σsm = false
  local adapted_Σrp
  local adapted_Σsm

  if (adapt_cov != 0 & mc.iterations > adapt_cov)
    if LinearAlgebra.isposdef(value(mc.Σrp))
      use_adapted_Σrp = true
      adapted_Σrp = OnlineStats.value(mc.Σrp) * 2.38^2 / length(mc.risk_parameters[1])
    end
    if LinearAlgebra.isposdef(value(mc.Σsm))
      use_adapted_Σsm = true
      adapted_Σsm = OnlineStats.value(mc.Σsm) * 2.38^2 / length(propertynames(mcmc.substitution_model))
    end
  end
  tr_cache = TransmissionRateCache(individuals(mcmc.population))
  for i = 1:n
    if use_adapted_Σrp & use_adapted_Σsm
      update!(mc, mcmc, adapted_Σrp, adapted_Σsm, σ, event_batches = event_batches, transmission_rate_cache=tr_cache)
    elseif use_adapted_Σrp & !use_adapted_Σsm
      update!(mc, mcmc, adapted_Σrp, Σsm, σ, event_batches = event_batches, transmission_rate_cache=tr_cache)
    elseif !use_adapted_Σrp & use_adapted_Σsm
      update!(mc, mcmc, Σrp, adapted_Σsm, σ, event_batches = event_batches, transmission_rate_cache=tr_cache)
    elseif !use_adapted_Σrp & !use_adapted_Σsm
      update!(mc, mcmc, Σrp, Σsm, σ, event_batches = event_batches, transmission_rate_cache=tr_cache)
    end
    next!(pmeter)
    @logmsg LogLevel(-5000) "MCMC progress" progress = i/n
    if adapt_cov > 0 && mod(i, adapt_cov) == 0
      OnlineStats.fit!(mc.Σrp, eachrow(convert(Array{Float64, 2}, mc.risk_parameters[end-adapt_cov:end])))
      OnlineStats.fit!(mc.Σsm, eachrow(convert(Array{Float64, 2}, mc.substitution_model[end-adapt_cov:end])))
      if LinearAlgebra.isposdef(value(mc.Σrp))
        if !use_adapted_Σrp
          use_adapted_Σrp = true
          @debug "Now using adapted covariance matrix for adaptive MCMC for Risk Function Parameters" core = Distributed.myid() UpdateFreq = adapt_cov
        end
        adapted_Σrp = OnlineStats.value(mc.Σrp) * 2.38^2 / length(mc.risk_parameters[1])
        @debug "Covariance matrix updated for Adaptive Metropolis-Hastings MCMC for Risk Function Parameters" Σrp = adapted_Σrp
      end
      if LinearAlgebra.isposdef(value(mc.Σsm))
        if !use_adapted_Σsm
          use_adapted_Σsm = true
          @debug "Now using adapted covariance matrix for adaptive MCMC on core $(Distributed.myid()) (updated every $adapt_cov iterations)"
        end
        adapted_Σsm = OnlineStats.value(mc.Σsm) * 2.38^2 / length(propertynames(mcmc.substitution_model))
        @debug "Covariance matrix updated for Adaptive Metropolis-Hastings MCMC" Σsm = adapted_Σsm
      end
    end
  end
  return mc
end

function iterate!(mc::MarkovChain{S, M},
                  mcmc::MCMC{S, M},
                  n::Int64,
                  Σrp::Array{Float64, 2},
                  σ::Float64,
                  progress_channel::RemoteChannel;
                  condition_on_network::Bool=false,
                  event_batches::Int64=1,
                  adapt_cov::Int64=100) where {
                  S <: DiseaseStateSequence,
                  M <: TNILM}
  if adapt_cov < 0
    @warn "`adapt_cov` argument indicates the increment in iterations in which the covariance matrix is updated and must be ≧ 0. Setting to 0 for non-adaptive Metropolis-Hastings."
    adapt_cov = 0
  end
  use_adapted_Σrp = false
  local adapted_Σrp
  if (adapt_cov != 0 & mc.iterations > adapt_cov) && LinearAlgebra.isposdef(value(mc.Σrp))
    use_adapted_Σrp = true
    adapted_Σrp = OnlineStats.value(mc.Σrp) * 2.38^2 / length(mc.risk_parameters[1])
  end
  tr_cache = TransmissionRateCache(individuals(mcmc.population))
  for i = 1:n
    if use_adapted_Σrp
      update!(mc, mcmc, adapted_Σrp, σ, condition_on_network = condition_on_network, event_batches = event_batches, transmission_rate_cache=tr_cache)
    else
      update!(mc, mcmc, Σrp, σ, condition_on_network = condition_on_network, event_batches = event_batches, transmission_rate_cache=tr_cache)
    end
    put!(progress_channel, true)
    @logmsg LogLevel(-5000) "MCMC progress" progress = i/n
    if adapt_cov > 0 && mod(i, adapt_cov) == 0
      OnlineStats.fit!(mc.Σrp, eachrow(convert(Array{Float64, 2}, mc.risk_parameters[end-adapt_cov:end])))
      if LinearAlgebra.isposdef(value(mc.Σrp))
        if !use_adapted_Σrp
          use_adapted_Σrp = true
          @debug "Now using adapted covariance matrix for adaptive MCMC on core $(Distributed.myid()) (updated every $adapt_cov iterations)"
        end
        adapted_Σrp = OnlineStats.value(mc.Σrp) * 2.38^2 / length(mc.risk_parameters[1])
        @debug "Covariance matrix updated for Adaptive Metropolis-Hastings MCMC" Covariance = adapted_Σrp
      end
    end
  end
  return mc
end

function iterate!(
  mc::MarkovChain{S, M},
  mcmc::MCMC{S, M},
  n::Int64,
  Σrp::Array{Float64, 2},
  Σsm::Array{Float64, 2},
  σ::Float64,
  progress_channel::RemoteChannel;
  event_batches::Int64=1,
  adapt_cov::Int64=100) where {
    S <: DiseaseStateSequence,
    M <: PhyloILM}
  if adapt_cov < 0
    @warn "`adapt_cov` argument indicates the increment in iterations in which the covariance matrix is updated and must be ≧ 0. Setting to 0 for non-adaptive Metropolis-Hastings."
    adapt_cov = 0
  end
  use_adapted_Σrp = false
  use_adapted_Σsm = false
  local adapted_Σrp
  local adapted_Σsm

  if (adapt_cov != 0 & mc.iterations > adapt_cov)
    if LinearAlgebra.isposdef(value(mc.Σrp))
      use_adapted_Σrp = true
      adapted_Σrp = OnlineStats.value(mc.Σrp) * 2.38^2 / length(mc.risk_parameters[1])
    end
    if LinearAlgebra.isposdef(value(mc.Σsm))
      use_adapted_Σsm = true
      adapted_Σsm = OnlineStats.value(mc.Σsm) * 2.38^2 / length(propertynames(mcmc.substitution_model))
    end
  end
  tr_cache = TransmissionRateCache(individuals(mcmc.population))
  for i = 1:n
    if use_adapted_Σrp & use_adapted_Σsm
      update!(mc, mcmc, adapted_Σrp, adapted_Σsm, σ, event_batches = event_batches, transmission_rate_cache=tr_cache)
    elseif use_adapted_Σrp & !use_adapted_Σsm
      update!(mc, mcmc, adapted_Σrp, Σsm, σ, event_batches = event_batches, transmission_rate_cache=tr_cache)
    elseif !use_adapted_Σrp & use_adapted_Σsm
      update!(mc, mcmc, Σrp, adapted_Σsm, σ, event_batches = event_batches, transmission_rate_cache=tr_cache)
    elseif !use_adapted_Σrp & !use_adapted_Σsm
      update!(mc, mcmc, Σrp, Σsm, σ, event_batches = event_batches, transmission_rate_cache=tr_cache)
    end
    put!(progress_channel, true)
    @logmsg LogLevel(-5000) "MCMC progress" progress = i/n
    if adapt_cov > 0 && mod(i, adapt_cov) == 0
      OnlineStats.fit!(mc.Σrp, eachrow(convert(Array{Float64, 2}, mc.risk_parameters[end-adapt_cov:end])))
      OnlineStats.fit!(mc.Σsm, eachrow(convert(Array{Float64, 2}, mc.substitution_model[end-adapt_cov:end])))
      if LinearAlgebra.isposdef(value(mc.Σrp))
        if !use_adapted_Σrp
          use_adapted_Σrp = true
          @debug "Now using adapted covariance matrix for adaptive MCMC for Risk Function Parameters" core = Distributed.myid() UpdateFreq = adapt_cov
        end
        adapted_Σrp = OnlineStats.value(mc.Σrp) * 2.38^2 / length(mc.risk_parameters[1])
        @debug "Covariance matrix updated for Adaptive Metropolis-Hastings MCMC for Risk Function Parameters" Σrp = adapted_Σrp
      end
      if LinearAlgebra.isposdef(value(mc.Σsm))
        if !use_adapted_Σsm
          use_adapted_Σsm = true
          @debug "Now using adapted covariance matrix for adaptive MCMC on core $(Distributed.myid()) (updated every $adapt_cov iterations)"
        end
        adapted_Σsm = OnlineStats.value(mc.Σsm) * 2.38^2 / length(propertynames(mcmc.substitution_model))
        @debug "Covariance matrix updated for Adaptive Metropolis-Hastings MCMC" Σsm = adapted_Σsm
      end
    end
  end
  return mc
end

function iterate!(
  mcmc::MCMC{S, M},
  n::Int64,
  Σrp::Array{Float64, 2},
  σ::Float64;
  condition_on_network::Bool=false,
  event_batches::Int64=1,
  adapt_cov::Int64=100) where {
    S <: DiseaseStateSequence,
    M <: TNILM}
  if length(mcmc.markov_chains) == 1
    iterate!(mcmc.markov_chains[1], mcmc, n, Σrp, σ,
             condition_on_network = condition_on_network,
             event_batches = event_batches,
             adapt_cov = adapt_cov)
    return mcmc
  else
    if adapt_cov < 0
      @warn "`adapt_cov` argument indicates the increment in iterations in which the covariance matrix is updated and must be ≧ 0. Setting to 0 for non-adaptive Metropolis-Hastings."
      adapt_cov = 0
    end
    pmeter = Progress(n*length(mcmc.markov_chains), "MCMC progress")
    pchannel = RemoteChannel(()->Channel{Bool}(10), 1)
    mc_Futures = Future[]
    @debug "Starting MCMC..."
    for mc in mcmc.markov_chains
      push!(mc_Futures, @spawn iterate!(mc, mcmc, n, Σrp, σ, pchannel, condition_on_network = condition_on_network, event_batches = event_batches, adapt_cov = adapt_cov))
    end
    @debug "MCMC in progress..."
    while !all(isready.(mc_Futures))
      take!(pchannel) && next!(pmeter)
    end
    close(pchannel)
    finish!(pmeter)
    mcmc.markov_chains = [fetch(i) for i in mc_Futures]
    @debug "MCMC complete..."
    return mcmc
  end
end

function iterate!(
  mcmc::MCMC{S, M},
  n::Int64,
  Σrp::Array{Float64, 2},
  Σsm::Array{Float64, 2},
  σ::Float64;
  condition_on_network::Bool=false,
  event_batches::Int64=1,
  adapt_cov::Int64=100) where {
    S <: DiseaseStateSequence,
    M <: PhyloILM}
  if length(mcmc.markov_chains) == 1
    iterate!(mcmc.markov_chains[1], mcmc, n, Σrp, Σsm, σ,
             event_batches = event_batches,
             adapt_cov = adapt_cov)
    return mcmc
  else
    if adapt_cov < 0
      @warn "`adapt_cov` argument indicates the increment in iterations in which the covariance matrix is updated and must be ≧ 0. Setting to 0 for non-adaptive Metropolis-Hastings."
      adapt_cov = 0
    end
    pmeter = Progress(n*length(mcmc.markov_chains), "MCMC progress")
    pchannel = RemoteChannel(()->Channel{Bool}(10), 1)
    mc_Futures = Future[]
    @debug "Starting MCMC..."
    for mc in mcmc.markov_chains
      push!(mc_Futures, @spawn iterate!(mc, mcmc, n, Σrp, Σsm, σ, pchannel, event_batches = event_batches, adapt_cov = adapt_cov))
    end
    @debug "MCMC in progress..."
    while !all(isready.(mc_Futures))
      take!(pchannel) && next!(pmeter)
    end
    close(pchannel)
    finish!(pmeter)
    mcmc.markov_chains = [fetch(i) for i in mc_Futures]
    @debug "MCMC complete..."
    return mcmc
  end
end

function iterate!(
  mc::MarkovChain{S, M},
  mcmc::MCMC{S, M},
  n::Int64,
  σ::Float64;
  condition_on_network::Bool=false,
  event_batches::Int64=1,
  adapt_cov::Int64=100) where {
    S <: DiseaseStateSequence,
    M <: TNILM}
  return iterate!(
    mc,
    mcmc,
    n,
    _initial_Σrp(mcmc),
    σ,
    condition_on_network = condition_on_network,
    event_batches = event_batches,
    adapt_cov = adapt_cov)
end

function iterate!(
  mc::MarkovChain{S, M},
  mcmc::MCMC{S, M},
  n::Int64,
  σ::Float64;
  event_batches::Int64=1,
  adapt_cov::Int64=100,
  DA_couple::Bool=false) where {
    S <: DiseaseStateSequence,
    M <: PhyloILM}
  return iterate!(
    mc,
    mcmc,
    n,
    _initial_Σrp(mcmc),
    _initial_Σsm(mcmc),
    σ,
    event_batches = event_batches,
    adapt_cov = adapt_cov,
    DA_couple = DA_couple)
end

function iterate!(
  mcmc::MCMC{S, M},
  n::Int64,
  σ::Float64;
  condition_on_network::Bool=false,
  event_batches::Int64=1,
  adapt_cov::Int64=100) where {
    S <: DiseaseStateSequence,
    M <: TNILM}
  return iterate!(
    mcmc,
    n,
    _initial_Σrp(mcmc),
    σ,
    condition_on_network = condition_on_network,
    event_batches = event_batches,
    adapt_cov = adapt_cov)
end

function iterate!(
  mcmc::MCMC{S, M},
  n::Int64,
  σ::Float64;
  event_batches::Int64=1,
  adapt_cov::Int64=100,
  DA_couple::Bool=false) where {
    S <: DiseaseStateSequence,
    M <: PhyloILM}
  return iterate!(
    mcmc,
    n,
    _initial_Σrp(mcmc),
    _initial_Σsm(mcmc),
    σ,
    event_batches = event_batches,
    adapt_cov = adapt_cov)
end


text:
function read_bgen(bgen_file::String)
    kwargs = Dict{Symbol, Any}(:sample_path => nothing, :idx_path => nothing)
    if bgen_file[end-3:end] == "bgen"
        base = bgen_file[1:end-4]

        samplefile = base * "sample"
        isfile(samplefile) ? kwargs[:sample_path] = samplefile : nothing

        bgifile = bgen_file * ".bgi"
        isfile(bgifile) ? kwargs[:idx_path] = bgifile : nothing
    end
    return Bgen(bgen_file; kwargs...)
end

function samples_genotype(probabilities, variant_genotypes, threshold=0.9)
    n = size(probabilities)[2]
    # The default value is missing
    t = Vector{Union{String, Missing}}(missing, n)
    for i in 1:n
        # If no allele has been annotated with sufficient confidence
        # the sample is declared as missing for this variant
        sample_gen_index = findfirst(x -> x >= threshold, probabilities[:, i])
        sample_gen_index isa Nothing || (t[i] = variant_genotypes[sample_gen_index])
    end
    return t
end

"""

A heterozygous genotype can be specified as (ALLELE₁, ALLELE₂) or (ALLELE₂, ALLELE₁).
Here we align this heterozygous specification on the query and default to
(ALLELE₁, ALLELE₂) provided in the BGEN file if nothing is specified in the query.
"""
function variant_genotypes(variant::Variant, queries)
    all₁, all₂ = alleles(variant)
    # Either (ALLELE₂, ALLELE₁) is provided in the query
    # and we return it as the heterozygous genotype.
    # Or the other specification will do in all other cases.
    queried_alleles = [q[Symbol(variant.rsid)] for (qn, q) in queries]
    if all₂*all₁ in collect(Iterators.flatten(queried_alleles))
        return [all₁*all₁, all₂*all₁, all₂*all₂]
    end
    return [all₁*all₁, all₁*all₂, all₂*all₂]
end

function UKBBGenotypes(queryfile, queries)
    config = TOML.parsefile(queryfile)
    snps = config["SNPS"]
    threshold = config["threshold"]
    # Let's load the variants by the files they are in
    bgen_groups = Dict()
    for (rsid, path) in snps
        haskey(bgen_groups, path) ? push!(bgen_groups[path], rsid) : bgen_groups[path] = [rsid]
    end

    genotypes = nothing
    for (path, rsids) in bgen_groups
        b = read_bgen(path)
        chr_genotypes = DataFrame(SAMPLE_ID=b.samples)

        # Iterate over variants in this chromosome
        for rsid in rsids
            v = variant_by_rsid(b, rsid)
            variant_gens = variant_genotypes(v, queries)
            probabilities = probabilities!(b, v)
            chr_genotypes[!, rsid] = samples_genotype(probabilities, variant_gens, threshold)
        end
        # I think concatenating should suffice but I still join as a safety
        genotypes isa Nothing ? genotypes = chr_genotypes :
            genotypes = innerjoin(genotypes, chr_genotypes, on=:SAMPLE_ID)
    end
    return genotypes
end

function preprocess(genotypes, confounders, phenotypes;
                    verbosity=1)

    # Make sure data SAMPLE_ID types coincide
    genotypes.SAMPLE_ID = string.(genotypes.SAMPLE_ID)
    confounders.SAMPLE_ID = string.(confounders.SAMPLE_ID)
    phenotypes.eid = string.(phenotypes.eid)

    # Join all elements together
    data = innerjoin(
            innerjoin(genotypes, confounders, on=:SAMPLE_ID),
            phenotypes,
            on = :SAMPLE_ID =>:eid
            )

    # Check no data has been lost at this stage
    any(nrows(array) != nrows(data) for array in (genotypes, confounders, phenotypes)) &&
        @warn """The number of matching samples is different in data sources: \n
            - Genotypes: $(nrows(genotypes)) \n
            - Confounders: $(nrows(confounders)) \n
            - Phenotypes: $(nrows(phenotypes)) \n
            - After Join: $(nrows(data)) \n
        """

    # Filter any line where an element is missing
    filtered_data = dropmissing(data)
    verbosity >= 1 && @info "Samples size after missing removed: $(nrows(filtered_data))"

    # Retrieve T and convert to categorical data
    T = filtered_data[!, filter(!=("SAMPLE_ID"), names(genotypes))]
    for name in names(T)
        T[!, name] = categorical(T[:, name])
    end

    # Retrieve W
    W = filtered_data[!, filter(!=("SAMPLE_ID"), names(confounders))]

    # Retrieve y which is assumed to be the first column after SAMPLE_ID
    # and convert to categorical array if needed
    y = filtered_data[!, filter(!=("eid"), names(phenotypes))][!, 1]

    # Only support binary and continuous traits for now
    is_binary(y) ? y = categorical(convert(Vector{Bool}, y)) : nothing

    return T, W, y
end

function PhenotypeTMLEEpistasis(tmles::Dict, T, W, y, queries; verbosity=1, phenotypename=nothing)
    # Build TMLE machine, can I update Y withoug altering the machine state
    # so that only the Q mach is fit again?
    tmle = is_binary(y) ? tmles["binary"] : tmles["continuous"]
    mach = machine(tmle, T, W, y)
    # Run TMLE
    fit!(mach; verbosity=verbosity-1)

    reports = briefreport(mach)
    results = DataFrame(
        PHENOTYPE=Symbol[],
        QUERYNAME=String[],
        QUERYSTRING=String[],
        ESTIMATE=Float64[],
        PVALUE=Float64[],
        LOWER_BOUND=Float64[],
        UPPER_BOUND=Float64[],
        STD_ERROR=Float64[],
        QSTACK_COEFS=String[],
        NROWS=Int[],
        TCOUNTS=String[]
        )

    tcounts = TreatmentCountsRepr(T)
    n = nrows(y)
    for (i, (queryname, query)) in enumerate(queries)
        querystring_ = querystring(query)

        queryreport = get_query_report(reports, i)
        pvalue = queryreport.pvalue
        lwb, upb = queryreport.confint
        estimate = queryreport.estimate
        stderror = queryreport.stderror
        qstack_coefs = repr_Qstack_coefs(mach)

        push!(results, (phenotypename, queryname, querystring_, estimate, pvalue, lwb, upb, stderror, qstack_coefs, n, tcounts))
    end
    return results
end

function PhenotypeCrossValidation(library::Dict, T, W, y, queries; verbosity=1, phenotypename=nothing)
    # Cross validate models for Q̅
    Qsettings = is_binary(y) ? library["Qcat"] : library["Qcont"]
    Qmeasure = is_binary(y) ? log_loss : l2
    hot_mach = machine(OneHotEncoder(drop_last=true), T)
    fit!(hot_mach, verbosity=verbosity)
    T_hot = transform(hot_mach, T)
    X = hcat(T_hot, W)
    Qresult_string = ""
    for (modelname, model) in Qsettings[2]
        mach = machine(model, X, y)
        result = evaluate!(mach, resampling=Qsettings[1], measure=Qmeasure, verbosity=verbosity)
        Qresult_string *= "$modelname: m=$(result.measurement[1]) std=$(std(result.per_fold[1])) | "
    end

    # Cross validate models for G
    Gsettings = library["G"]
    Gresult_string = ""
    t_target = encode(T)
    for (modelname, model) in Gsettings[2]
        mach = machine(model, W, t_target)
        result = evaluate!(mach, resampling=Gsettings[1], measure=log_loss, operation=predict, verbosity=verbosity)
        Gresult_string *= "$modelname: m=$(result.measurement[1]) std=$(std(result.per_fold[1])) | "
    end

    variants = join(keys(queries[1][2]), " && ")

    return DataFrame(
        PHENOTYPE=[phenotypename],
        VARIANTS=[variants],
        Q_RESULTSTRING=[Qresult_string[1:length(Qresult_string)-3]],
        G_RESULTSTRING=[Gresult_string[1:length(Gresult_string)-3]]
        )

end

function UKBBVariantRun(parsed_args; run_fn=PhenotypeTMLEEpistasis)
    v = parsed_args["verbosity"]

    # Parse queries
    queries = parse_queries(parsed_args["queries"])

    # Build estimators
    tmle_config = TOML.parsefile(parsed_args["estimator"])
    estimators = estimators_from_toml(tmle_config, queries, run_fn)

    v >= 1 && @info "Loading Genotypes and Confounders."
    # Build Genotypes
    genotypes = UKBBGenotypes(parsed_args["queries"], queries)

    # Read Confounders
    confounders = CSV.File(parsed_args["confounders"]) |> DataFrame

    # Generate phenotypes range
    all_phenotype_names = phenotypesnames(parsed_args["phenotypes"])
    done_phenotypes = init_or_retrieve_results(parsed_args["output"], run_fn)
    phenotypes_range = phenotypes_list(
        parsed_args["phenotypes-list"],
        done_phenotypes,
        all_phenotype_names
        )

    for phenotypename in phenotypes_range
        v >= 1 && @info "Running procedure with phenotype: $phenotypename."
        # Read Target
        phenotype = CSV.File(parsed_args["phenotypes"], select=[:eid, phenotypename]) |> DataFrame
        # Preprocess all data together
        T, W, y = preprocess(genotypes,
                            confounders,
                            phenotype;
                            verbosity=v)
        # Run the estimation procedure
        results = run_fn(estimators, T, W, y, queries; verbosity=v, phenotypename=phenotypename)
        # Update the results
        CSV.write(parsed_args["output"], results, append=true)

    end

    v >= 1 && @info "Done."
end

ext: jl


text:
@testset "maxintfloat $T" for T in (Double32, Double64)
    @test isinteger(maxintfloat(T))

    # Previous integer is representable, next integer is not
    @test maxintfloat(T) == (maxintfloat(T) - one(T)) + one(T)
    @test maxintfloat(T) != (maxintfloat(T) + one(T)) - one(T)
end

# Double16 is a special case
@test isinteger(maxintfloat(Double16))
@test maxintfloat(Double16) == floatmax(Double16)

@testset "Inf and NaN generation $T" for T in (Double16, Double32, Double64)
    @test T(Inf) == inf(T)
    @test T(Inf) == posinf(T)
    @test T(-Inf) == neginf(T)
    @test T(NaN) == nan(T)

    @test T(Inf32) == inf(T)
    @test T(Inf32) == posinf(T)
    @test T(-Inf32) == neginf(T)
    @test T(NaN32) == nan(T)

    @test T(Inf16) == inf(T)
    @test T(Inf16) == posinf(T)
    @test T(-Inf16) == neginf(T)
    @test T(NaN16) == nan(T)
end


text:
"""
	Box3d(p1::Point3d,p2::Point3d)

Mutable `struct`, with two fields `p1::Point3d` and `p2::Point3d`.
Two constructors `Box3d()` and `Box3d(p1::Point3d,p2::Point3d)`
# Example
```

```
"""
mutable struct Box3d
	p1::Point3d
	p2::Point3d

	# constyructor
	function Box3d()
		new(Point3d(0,0,0),Point3d(0,0,0))
	end

	# constructor
	function Box3d(p1::Point3d,p2::Point3d)
		new(p1,p2)
	end

end

Base.:(==)(a::Box3d, b::Box3d) = a.p1 == b.p1 && a.p2 == b.p2

"""
	invalidBox()

Returns an object of type `ViewerGL.Box3d`.
# Example
```
julia> invalidBox().p1
3-element StaticArrays.MArray{Tuple{3},Float64,1,3}:
 Inf
 Inf
 Inf

 julia> invalidBox().p2
 3-element StaticArrays.MArray{Tuple{3},Float64,1,3}:
  -Inf
  -Inf
  -Inf
```
"""
function invalidBox()
	m,M=typemin(Float64),typemax(Float64)
	return Box3d(Point3d(M,M,M),Point3d(m,m,m))
end

"""
	addPoint(box::GL.Box3d,point::GL.Point3d)

Update a mutable `Box3d`.
# Example
```
julia> box = GL.addPoint(GL.invalidBox(), GL.Point3d(1.0,2.0,0.0))
ViewerGL.Box3d([1.0, 2.0, 0.0], [1.0, 2.0, 0.0])
```
"""
function addPoint(box::Box3d,point::Point3d)
	for i in 1:3
		box.p1[i]=min(box.p1[i],point[i])
		box.p2[i]=max(box.p2[i],point[i])
	end
	return box
end

"""
	getPoints(box::Box3d)

Return the 8 vertices of an object of `Box3d` type.
# Example
```
julia> box = GL.getPoints(GL.Box3d(GL.Point3d(0,0,0),GL.Point3d(1,1,1)))
8-element Array{StaticArrays.MArray{Tuple{3},Float64,1,3},1}:
 [0.0, 0.0, 0.0]
 [1.0, 0.0, 0.0]
 [1.0, 1.0, 0.0]
 [0.0, 1.0, 0.0]
 [0.0, 0.0, 1.0]
 [1.0, 0.0, 1.0]
 [1.0, 1.0, 1.0]
 [0.0, 1.0, 1.0]
```
"""
function getPoints(box::Box3d)
	return [
		Point3d(box.p1[1],box.p1[2],box.p1[3]),
		Point3d(box.p2[1],box.p1[2],box.p1[3]),
		Point3d(box.p2[1],box.p2[2],box.p1[3]),
		Point3d(box.p1[1],box.p2[2],box.p1[3]),
		Point3d(box.p1[1],box.p1[2],box.p2[3]),
		Point3d(box.p2[1],box.p1[2],box.p2[3]),
		Point3d(box.p2[1],box.p2[2],box.p2[3]),
		Point3d(box.p1[1],box.p2[2],box.p2[3])
	]
end


text:
#Pkg.add("Ipopt");Pkg.build("Ipopt")
Pkg.add("PyPlot");Pkg.build("PyPlot")
#Pkg.add("PGFPlots");Pkg.build("PGFPlots")
#Pkg.add("GR");Pkg.build("GR")

using Documenter,NLOptControl
makedocs(modules=[NLOptControl],
        doctest=false, clean=true,
        format =:html,
        authors="Huckleberry Febbo",
        assets = ["assets/style.css"],
        sitename="NLOptControl.jl",
        pages = Any[
        "Home" => "index.md",
        "Tutorials"=>Any[
              "tutorials/Brachistochrone/main.md",
              "tutorials/MoonLander/main.md",
              "tutorials/BrysonDenham/main.md",
              "tutorials/Beam/main.md",
              "tutorials/HyperSensitive/main.md",
              "tutorials/RobotArm/main.md",
              "tutorials/Rocket/main.md",
              "tutorials/Unicycle/main.md",
               ],
       "MPC"=>Any[
           "mpc/index.md"
           ]
               ]
               )

deploydocs(
    deps=nothing,
    repo="github.com/JuliaMPC/NLOptControl.jl.git",
    target="build",
    osname="linux",
    julia="1.0",
    make=nothing)

#=

     "Exported Functions"=>Any[
     "functions/NLOptControl.md",
     "functions/PrettyPlots.md"
     ]
     ]=#


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

function rectangle_path(lo::Vector{<:Real}, hi::Vector{<:Real})
    [
        lo[1] lo[2];
        hi[1] lo[2];
        hi[1] hi[2];
        lo[1] hi[2];
        lo[1] lo[2];
    ]
end

function err_ellipsis_path(μ::Vector{<:Real}, Σ::Matrix{<:Real}, confidence::Real = 0.68, npts = 256)
    σ_sqr, A = eigen(Hermitian(Σ))
    σ = sqrt.(σ_sqr)
    ϕ = range(0, stop = 2π, length = 100)
    σ_scaled = σ .* sqrt(invlogcdf(Chisq(2), log(confidence)))
    xy = hcat(σ_scaled[1] * cos.(ϕ), σ_scaled[2] * sin.(ϕ)) * [A[1,1] A[1,2]; A[2,1] A[2,2]]
    xy .+ μ'
end

@recipe function f(stats::MCMCBasicStats, parsel::NTuple{2,Integer})
    pi_x, pi_y = parsel

    Σ_all = stats.param_stats.cov
    Σ = [Σ_all[pi_x, pi_x] Σ_all[pi_x, pi_y]; Σ_all[pi_y, pi_x] Σ_all[pi_y, pi_y]]

    μ = stats.param_stats.mean[[pi_x, pi_y]]
    mode_xy = stats.mode[[pi_x, pi_y]]

    conf = standard_confidence_vals

    color --> :darkviolet
    alpha --> 0.75
    linewidth --> 2

    for i in eachindex(conf)
        xy = err_ellipsis_path(μ, Σ, conf[i])
        @series begin
            seriestype := :path
            label --> "$(100 *conf[i])%"
            (xy[:,1], xy[:,2])
        end
    end

    markercolor --> :white
    markeralpha --> 0
    markersize --> 7
    markerstrokecolor --> :black
    markerstrokealpha --> 1
    markerstrokewidth --> 2

    # @series begin
    #     seriestype := :scatter
    #     label := "mean"
    #     markershape := :circle
    #     ([μ[1]], [μ[2]])
    # end

    # @series begin
    #     seriestype := :scatter
    #     label := "mode"
    #     markershape := :rect
    #     ([mode_xy[1]], [mode_xy[2]])
    # end

    vlo = stats.param_stats.minimum[[pi_x, pi_y]]
    vhi = stats.param_stats.maximum[[pi_x, pi_y]]
    rect_xy = rectangle_path(vlo, vhi)

    @series begin
        seriestype := :path
        label --> "bbox"
        (rect_xy[:,1], rect_xy[:,2])
    end

    nothing
end


text:
import ..Abstract.AbstractMolecularContainer

type Molecule <: AbstractMolecularContainer
    atoms::Array{Atom,1}
end

include("constructor.jl")


text:
operators = [:-
            :+
            :*
            :/
            :^
            :exp
            :log
            :Expect
            :max
            :min
            :∫
            :≤]

function subs!(x::Expr,list::Dict)
    for i = 1:length(x.args)
        if in(x.args[i],keys(list))
            x.args[i] = list[x.args[i]]
        elseif isa(x.args[i],Expr)
            subs!(x.args[i],list)
        end
    end
end

function subs(x1::Expr,list::Dict)
    x = deepcopy(x1)
    subs!(x,list)
    return x
end

function subs!(x::Expr,s::Pair)
    for i = 1:length(x.args)
        if x.args[i]==s.first
            x.args[i] = s.second
        elseif isa(x.args[i],Expr)
            subs!(x.args[i],s)
        end
    end
end

function subs(x1::Expr,s::Pair)
    x = deepcopy(x1)
    subs!(x,s)
    return x
end

function addindex!(x,ignore=operators)
    if typeof(x) == Expr
        if x.head == :ref
            return x
        else
            for i = 1:length(x.args)
                x.args[i]=addindex!(x.args[i],ignore)
            end
        end
    elseif typeof(x) == Symbol
        if !in(x,ignore)
            x= :($x[0])
            return x
        end
    end
    return x
end

addindex(x,ignore=operators) = addindex!(deepcopy(x),ignore)

function removeindex!(x)
    if isa(x,Expr)
        if x.head == :ref
            return x.args[1]
        else
            for i = 1:length(x.args)
                x.args[i]=removeindex!(x.args[i])
            end
        end
    end
    return x
end

removeindex(x,ignore=operators) = removeindex!(deepcopy(x),ignore)

function tchange!(x::Expr,t::Int,ignore=operators)
    if x.head == :ref
        if (x.args[2] == -1 && t==1) || (x.args[2] == 1 && t==-1)
            x = :($(x.args[1])[0])
        elseif x.args[2]==0
            x.args[2] = t
        end
    else
        for i = 1:length(x.args)
            if isa(x.args[i],Expr)
                x.args[i] = tchange!(x.args[i],t,ignore)
            end
        end
    end
    x
end

tchange(x::Expr,t::Int,ignore=operators) = tchange!(deepcopy(x),t,ignore)

function getv(x,list=Expr[],ignore=operators)
    if isa(x,Expr)
        if x.head==:ref
            push!(list,x)
        else
            for i = 1:length(x.args)
                list = getv(x.args[i],list,ignore)
            end
        end
    end
    return unique(list)
end


text:
export iauCr
"""
Copy an r-matrix.

This function is part of the International Astronomical Union's
SOFA (Standards Of Fundamental Astronomy) software collection.

Status:  vector/matrix support function.

Given:
   r        double[3][3]    r-matrix to be copied

Returned:
   c        double[3][3]    copy

Called:
   iauCp        copy p-vector

This revision:  2016 May 19

SOFA release 2018-01-30

Copyright (C) 2018 IAU SOFA Board.  See notes at end.
"""
# void iauCr(double r[3][3], double c[3][3])

function iauCr(r::AbstractMatrix{<:Real})

   # Allocate return value
   c = zeros(Float64, 3, 3)

   ccall((:iauCr, libsofa_c), Cvoid,
         (Ptr{Cdouble}, Ptr{Cdouble}),
         convert(Array{Float64, 2}, r'), c)

   return SMatrix{3,3}(c')
end


text:
using LinearAlgebra

@testset "generators" begin
    @testset "erdos_renyi" begin
        @test_throws DomainError erdos_renyi_undirected(10,11.)
        @test_throws DomainError erdos_renyi_directed(10,11.)

        n = 100
        avgdegs = range(1., stop=2*log(n), length=100)
        compsizes = map( (dbar) ->
                maximum(scomponents(erdos_renyi_undirected(n,dbar)).sizes),
            avgdegs )

        @test is_undirected(erdos_renyi_undirected(10,2.))

        erdos_renyi_undirected(0, 0.)
        erdos_renyi_undirected(0, 0)
        erdos_renyi_undirected(5, 0)
        erdos_renyi_undirected(5, 0.)
        erdos_renyi_undirected(5, 1.)
        erdos_renyi_undirected(5, 1)
        @test all(diag(sparse_transpose(erdos_renyi_undirected(10, 0.5))) .== 0.)
        erdos_renyi_directed(0, 0)
        erdos_renyi_directed(0, 0)
        erdos_renyi_directed(5, 0.)
        erdos_renyi_directed(5, 1.)
        @test all(diag(sparse_transpose(erdos_renyi_directed(10, 0.5))) .== 0.)
    end

    @testset "chung_lu" begin
        @test_throws ArgumentError chung_lu_undirected([1,2],2)
        @test_throws ArgumentError chung_lu_undirected([1,3])
        @test size(chung_lu_undirected([1,2]),1) == 2
    end

    @testset "havel_hakimi" begin
        @test is_graphical_sequence([1,1]) == true
        @test is_graphical_sequence([0,1,1]) == true
        @test is_graphical_sequence([2,2,2]) == true
        @test is_graphical_sequence([1,1,2]) == true
        @test is_graphical_sequence([1,2,2]) == false
        @test is_graphical_sequence([1,1,3]) == false
        @test is_graphical_sequence([0]) == true
        @test is_graphical_sequence(Int[]) == true
        @test is_graphical_sequence([1,1,1,1]) == true

        @test Matrix(sparse_transpose(havel_hakimi_graph(Int[]))) == zeros(Bool,0,0)
        @test Matrix(sparse_transpose(havel_hakimi_graph([0]))) == zeros(Bool,1,1)
        @test Matrix(sparse_transpose(havel_hakimi_graph([1,1]))) == [0 1; 1 0]
        @test Matrix(sparse_transpose(havel_hakimi_graph([2,2,2]))) == [0 1 1; 1 0 1; 1 1 0]
        six_star = [0 1 1 1 1 1; 1 0 0 0 0 0; 1 0 0 0 0 0; 1 0 0 0 0 0; 1 0 0 0 0 0; 1 0 0 0 0 0]
        @test Matrix(sparse_transpose(havel_hakimi_graph([5,1,1,1,1,1]))) == six_star

        @test_throws ArgumentError is_graphical_sequence([-1])
        @test_throws ArgumentError havel_hakimi_graph([1,2,2])
    end

    @testset "pa_graph" begin
        @test typeof(pa_graph(10,5,5)) == MatrixNetwork{Bool}
        @test typeof(pa_edges!(2,1,[(1,1)])) == Vector{Tuple{Int,Int}}
        @test_throws ArgumentError pa_edges!(5,2,Vector{Tuple{Int,Int}}())
        @test is_empty(pa_graph(0,0,0)) == true
        @test all(diag(sparse_transpose(pa_graph(10, 2, 3))) .== 0.)
        @test is_undirected(pa_graph(10, 12, 3))
        @test_throws ArgumentError pa_graph(-1,10,3)
        @test_throws ArgumentError pa_graph(5,10,-3)
        @test maximum(map(first, pa_edges!(5,1,[(1,2),(2,1)],2))) == 7
        @test size(pa_graph(10,5,5),1) == 10
        @test size(pa_graph(10,-5,5),1) == 10
        @test nnz(sparse_transpose(pa_graph(10,-5,5))) == 20
    end

    @testset "gpa_graph" begin
        @test typeof(gpa_graph(10,0.5,0.3,2)) == MatrixNetwork{Bool}
        @test typeof(gpa_graph(10,0.5,0.3,2,Val{true})) == MatrixNetwork{Bool}
        @test typeof(gpa_edges!(4,.75,.25,[(1,2)],1)) == Vector{Tuple{Int,Int}}
        @test_throws ArgumentError gpa_edges!(4,.75,.25,[(1,1)],1)
        @test typeof(gpa_edges!(4,.75,.25,[(1,1)],1,Val{true})) == Vector{Tuple{Int,Int}}
        @test is_empty(gpa_graph(0,0.0,0.0,0)) == true
        @test all(diag(sparse_transpose(gpa_graph(10, .3,.4, 3))) .== 0.)
        @test is_undirected(gpa_graph(10, .2, .8,4))
        @test_throws ArgumentError gpa_graph(-1,.3,.4,5)
        @test_throws ArgumentError gpa_graph(8,.8,.4,2)
        @test_throws ArgumentError gpa_graph(10,.3,.4,-2)
        @test size(gpa_graph(12,.5,.4,2),1) == 12
        @test size(gpa_graph(21,.1,.1,3),1) == 21
    end

    @testset "roach_graph" begin
        @test_throws ArgumentError roach_graph(-1)
        @test is_empty(roach_graph(0))
        @test is_connected(roach_graph(5))
        @test is_connected(roach_graph(5, Val{true})[1])
        @test maximum(bfs(roach_graph(5, Val{false}),20)[1]) == 11
        @test maximum(bfs(roach_graph(10),40)[1]) == 21
        @test (bfs(roach_graph(10),40)[1])[1] == 20

        G,xy = roach_graph(5, Val{true})
        filt = vec(all(xy .>= 0,dims = 2))
        Asub = sparse_transpose(G)[filt,filt]
        # Asub should be an antennae...
        @test bfs(Asub,1)[1][5] == 4
        filt = vec(all(xy .<= 0,dims = 2))
        Asub = sparse_transpose(G)[filt,filt]
        # Asub should be an line...
        @test bfs(Asub,1)[1][5] == 4
    end

    @testset "lollipop_graph" begin
        @test_throws ArgumentError lollipop_graph(-1)
        @test_throws ArgumentError lollipop_graph(-1,-1)
        @test_throws ArgumentError lollipop_graph(5,-1,Val{true})

        @test is_connected(lollipop_graph(5))
        @test is_connected(lollipop_graph(5, Val{true})[1])
        @test is_connected(lollipop_graph(5, 10))
        @test is_connected(lollipop_graph(10, 5))

        @test maximum(bfs(lollipop_graph(5,10),1)[1]) == 6
        @test maximum(bfs(lollipop_graph(10,5),1)[1]) == 11

        G,xy = lollipop_graph(6, 5, Val{true})
        filt = vec(all(xy .<= 0, dims = 2))
        Asub = sparse_transpose(G)[filt,filt]
        # Asub should be an tail (or a line graph
        @test bfs(Asub,1)[1][6] == 5
        Asub = sparse_transpose(G)[.~filt,.~filt]
        @test nnz(Asub) == 5*4

    end
end


text:
module Solutions

    using HDF5
    using OffsetArrays
    using SharedArrays

    using GeometricBase
    using GeometricBase.Config
    using GeometricBase.Utils
    using GeometricEquations

    using GeometricBase: fromarray

    export DEFAULT_NSAVE, DEFAULT_NWRITE

    const DEFAULT_NSAVE = 1
    const DEFAULT_NWRITE = 0

    export DataSeries, TimeSeries, Solution, DeterministicSolution
    export get_solution, get_solution!, set_solution!

    include("solutions/solution.jl")

    export AtomicSolution,
           AtomicSolutionODE, AtomicSolutionPODE,
           AtomicSolutionDAE, AtomicSolutionPDAE
    export update!, cut_periodic_solution!

    include("solutions/atomic_solution.jl")
    include("solutions/atomic_solution_ode.jl")
    include("solutions/atomic_solution_pode.jl")
    include("solutions/atomic_solution_dae.jl")
    include("solutions/atomic_solution_pdae.jl")

    export SolutionODE, SSolutionODE, PSolutionODE, SolutionPODE, SSolutionPODE, PSolutionPODE
    export SolutionDAE, SSolutionDAE, PSolutionDAE, SolutionPDAE, SSolutionPDAE, PSolutionPDAE
    export get_initial_conditions, get_initial_conditions!, set_initial_conditions!
    export create_hdf5, create_hdf5!, write_to_hdf5, hdf5

    include("solutions/solution_ode.jl")
    include("solutions/solution_pode.jl")
    include("solutions/solution_dae.jl")
    include("solutions/solution_pdae.jl")

    include("solutions/solution_hdf5.jl")

    include("solutions/atomic_solution_constructors.jl")

end


text:
using DelimitedFiles
using Random
using Serialization

export Basket, convertSparseCsvToBaskets, convertBasketsToSparseCsv

struct Basket
  basketItems::Vector{Int}
  heldOutItem::Int
  numItemsInCatalog::Int
end

# Converts data from a file in sparse CSV basket format to the data format
# (sparse vector representation) expected by the DPP learning implementation
function convertSparseCsvToBaskets(csvBasketDataFileName,
  trainingBasketsDictFileName, testBasketsDictFileName, trainingSetSizePercent)
  Random.seed!(1234)

  csvBasketData = readdlm(csvBasketDataFileName, ',')
  numCsvRecords = size(csvBasketData, 1)
  numCsvColumns = size(csvBasketData, 2)
  numItems = 0
  basketItemsDict = Dict{Int, Vector{Int}}()

  # Populate basketItemsDict with baskets from csvBasketData
  basketCounter = 1
  for i = 1:numCsvRecords
    j = 1
    basket = Int[]

    for j = 1:numCsvColumns
      if csvBasketData[i, j] != ""
        itemId = csvBasketData[i, j]
        push!(basket, itemId)

        # Update the number of items in the catalog
        if itemId > numItems
          numItems = itemId
        end
      else
        # We've reached the end of the item entries for this basket
        break
      end
    end

    basketItemsDict[basketCounter] = basket

    basketCounter += 1
  end

  println("num items in catalog: $numItems")

  # Generate histogram of number of baskets for each basket length (item count)
  numBasketsPerBasketLengthHistDict = Dict{Int, Int}()
  for basketId in keys(basketItemsDict)
    itemsInBasket = basketItemsDict[basketId]
    itemCountForBasket = length(itemsInBasket)

    if haskey(numBasketsPerBasketLengthHistDict, itemCountForBasket)
      # There is already at least one basket with this item count
      numBasketsWithItemCount = numBasketsPerBasketLengthHistDict[itemCountForBasket]
      numBasketsWithItemCount += 1
      numBasketsPerBasketLengthHistDict[itemCountForBasket] = numBasketsWithItemCount
    else
      numBasketsPerBasketLengthHistDict[itemCountForBasket] = 1
    end
  end

  for itemCountForBasket in sort(collect(keys(numBasketsPerBasketLengthHistDict)))
    println("Item count $itemCountForBasket: $(numBasketsPerBasketLengthHistDict[itemCountForBasket])")
  end

  # Remove baskets of size 1 from basketItemsDict, since we are not interested in them for a DPP
  for basketId in keys(basketItemsDict)
    itemsInBasket = basketItemsDict[basketId]
    if length(itemsInBasket) == 1
      delete!(basketItemsDict, basketId)
    end
  end

  # Generate training and test sets
  trainingBasketsDict = Dict{Int, Basket}()
  testBasketsDict = Dict{Int, Basket}()
  numBaskets = length(collect(keys(basketItemsDict)))
  trainingSetNumBaskets = convert(Int, round(numBaskets * trainingSetSizePercent));
  testSetNumBaskets = numBaskets - trainingSetNumBaskets
  shuffledBasketItemsDictKeys = shuffle(collect(keys(basketItemsDict)))
  trainingSetBaskets = shuffledBasketItemsDictKeys[1:trainingSetNumBaskets]
  testSetBaskets = shuffledBasketItemsDictKeys[(trainingSetNumBaskets + 1):numBaskets]
  for basketId in trainingSetBaskets
    basketItems = shuffle(basketItemsDict[basketId])
    heldOutItem = pop!(basketItems)
    basket = Basket(basketItems, heldOutItem, numItems)
    trainingBasketsDict[basketId] = basket
  end
  for basketId in testSetBaskets
    basketItems = shuffle(basketItemsDict[basketId])
    heldOutItem = pop!(basketItems)
    basket = Basket(basketItems, heldOutItem, numItems)
    testBasketsDict[basketId] = basket
  end

  println("Num baskets in training set: $(length(collect(keys(trainingBasketsDict))))")
  println("Num baskets in test set: $(length(collect(keys(testBasketsDict))))")

  open(f -> serialize(f, trainingBasketsDict), trainingBasketsDictFileName, "w");
  println("Saved trainingBasketsDict")

  open(f -> serialize(f, testBasketsDict), testBasketsDictFileName, "w");
  println("Saved testBasketsDict")
end

# Converts data from a file in sparse vector (DPP implementation) format to
# sparse CSV format
function convertBasketsToSparseCsv(basketsDictFileName, basketsDictObjectName,
                                   csvBasketDataFileName)
  # Load basket data
  basketsDict = load(basketsDictFileName, basketsDictObjectName)
  println("Loaded $basketsDictFileName")

  # Write basket data to csvBasketDataFileName
  csvBasketDataFile = open(csvBasketDataFileName, "w")
  for basketId in collect(keys(basketsDict))
    basket = basketsDict[basketId]
    basketArray = basket.basketItems
    push!(basketArray, basket.heldOutItem)

    println("Contents for basket $basketId: $basketArray")
    write(csvBasketDataFile, join(sort(basketArray), " "), "\n")
  end
  close(csvBasketDataFile)
end


text:
# AUTO GENERATED FILE - DO NOT EDIT

export 'feffery'_antdcountdown

"""
    'feffery'_antdcountdown(;kwargs...)

An AntdCountdown component.

Keyword arguments:
- `id` (String; optional)
- `className` (String; optional)
- `format` (String; optional)
- `loading_state` (optional): . loading_state has the following type: lists containing elements 'is_loading', 'prop_name', 'component_name'.
Those elements have the following types:
  - `is_loading` (Bool; optional): Determines if the component is loading or not
  - `prop_name` (String; optional): Holds which property is loading
  - `component_name` (String; optional): Holds the name of the component that is loading
- `prefix` (optional): . prefix has the following type: lists containing elements 'mode', 'content'.
Those elements have the following types:
  - `mode` (a value equal to: 'text', 'icon'; optional)
  - `content` (String; optional)
- `style` (Dict; optional)
- `suffix` (optional): . suffix has the following type: lists containing elements 'mode', 'content'.
Those elements have the following types:
  - `mode` (a value equal to: 'text', 'icon'; optional)
  - `content` (String; optional)
- `title` (String; optional)
- `value` (String; optional)
- `valueFormat` (String; optional)
- `valueStyle` (Dict; optional)
"""
function 'feffery'_antdcountdown(; kwargs...)
        available_props = Symbol[:id, :className, :format, :loading_state, :prefix, :style, :suffix, :title, :value, :valueFormat, :valueStyle]
        wild_props = Symbol[]
        return Component("'feffery'_antdcountdown", "AntdCountdown", "feffery_antd_components", available_props, wild_props; kwargs...)
end

max_stars_repo_name: ReilYoung/feffery-antd-components

max_issues_repo_name: ReilYoung/feffery-antd-components

max_forks_repo_name: ReilYoung/feffery-antd-components


text:
# Nest
#
# If the line exceeds the print width it will be nested.
#
# This is done by replacing `PLACEHOLDER` nodes with `NEWLINE`
# nodes and updating the PTree's indent.
#
# `extra_width` provides additional width to consider from
# the top-level node.
#
# Example:
#
#     arg1 op arg2
#
# the length of " op" will be considered when nesting arg1

function skip_indent(x::PTree)
    if x.typ === CSTParser.LITERAL && x.val == ""
        return true
    elseif x.typ === NEWLINE || x.typ === NOTCODE
        return true
    end
    false
end

function walk(f, x::PTree, s::State)
    f(x, s)
    is_leaf(x) && return
    for (i, n) in enumerate(x.nodes)
        if n.typ === NEWLINE && i < length(x.nodes)
            if is_closer(x.nodes[i+1])
                s.line_offset = x.nodes[i+1].indent
            elseif !skip_indent(x.nodes[i+1])
                s.line_offset = x.indent
            end
        else
            walk(f, n, s)
        end
    end
end

function reset_line_offset!(x::PTree, s::State)
    !is_leaf(x) && return
    s.line_offset += length(x)
end

function add_indent!(x::PTree, s::State, indent)
    indent == 0 && return
    lo = s.line_offset
    f = (x::PTree, s::State) -> x.indent += indent
    walk(f, x, s)
    s.line_offset = lo
end

function nest!(nodes::Vector{PTree}, s::State, indent; extra_width = 0)
    for (i, n) in enumerate(nodes)
        if n.typ === NEWLINE && nodes[i+1].typ === CSTParser.Block
            s.line_offset = nodes[i+1].indent
        elseif n.typ === NOTCODE && nodes[i+1].typ === CSTParser.Block
            s.line_offset = nodes[i+1].indent
        elseif n.typ === NEWLINE
            s.line_offset = indent
        else
            nest!(n, s, extra_width = extra_width)
        end
    end
end

function nest!(x::PTree, s::State; extra_width = 0)
    if is_leaf(x)
        s.line_offset += length(x)
        return
    end

    if x.typ === CSTParser.Import
        n_import!(x, s, extra_width = extra_width)
    elseif x.typ === CSTParser.Export
        n_import!(x, s, extra_width = extra_width)
    elseif x.typ === CSTParser.Using
        n_import!(x, s, extra_width = extra_width)
    elseif x.typ === CSTParser.ImportAll
        n_import!(x, s, extra_width = extra_width)
    elseif x.typ === CSTParser.WhereOpCall
        n_wherecall!(x, s, extra_width = extra_width)
    elseif x.typ === CSTParser.ConditionalOpCall
        n_condcall!(x, s, extra_width = extra_width)
    elseif x.typ === CSTParser.BinaryOpCall
        n_binarycall!(x, s, extra_width = extra_width)
    elseif x.typ === CSTParser.Curly
        n_call!(x, s, extra_width = extra_width)
    elseif x.typ === CSTParser.Call
        n_call!(x, s, extra_width = extra_width)
    elseif x.typ === CSTParser.MacroCall
        n_call!(x, s, extra_width = extra_width)
    elseif x.typ === CSTParser.Ref
        n_call!(x, s, extra_width = extra_width)
    elseif x.typ === CSTParser.ChainOpCall
        n_chainopcall!(x, s, extra_width = extra_width)
    elseif x.typ === CSTParser.Comparison
        n_comparison!(x, s, extra_width = extra_width)
    elseif x.typ === CSTParser.TupleH
        n_tuple!(x, s, extra_width = extra_width)
    elseif x.typ === CSTParser.Vect
        n_tuple!(x, s, extra_width = extra_width)
    elseif x.typ === CSTParser.Vcat
        n_tuple!(x, s, extra_width = extra_width)
    elseif x.typ === CSTParser.Block
        n_block!(x, s, extra_width = extra_width)
    elseif x.typ === CSTParser.For
        n_for!(x, s, extra_width = extra_width)
    elseif x.typ === CSTParser.Let
        n_for!(x, s, extra_width = extra_width)
    elseif x.typ === CSTParser.TypedVcat
        n_call!(x, s, extra_width = extra_width)
    elseif x.typ === CSTParser.Parameters
        n_tuple!(x, s, extra_width = extra_width)
    elseif x.typ === CSTParser.Braces
        n_tuple!(x, s, extra_width = extra_width)
    elseif x.typ === CSTParser.InvisBrackets
        n_tuple!(x, s, extra_width = extra_width)
    elseif x.typ === CSTParser.UnaryOpCall && x.nodes[2].typ === CSTParser.OPERATOR
        nest!(x.nodes[1], s, extra_width = extra_width + length(x.nodes[2]))
        nest!(x.nodes[2], s)
    elseif x.typ === CSTParser.Do
        n_do!(x, s, extra_width = extra_width)
    elseif x.typ === CSTParser.Generator
        n_gen!(x, s, extra_width = extra_width)
    elseif x.typ === CSTParser.Filter
        n_gen!(x, s, extra_width = extra_width)
    elseif x.typ === CSTParser.Comprehension
        n_tuple!(x, s, extra_width = extra_width)
    else
        nest!(x.nodes, s, x.indent, extra_width = extra_width)
    end
end

function n_do!(x, s; extra_width = 0)
    ew = sum(length.(x.nodes[2:3]))
    # make sure there are nodes after "do"
    if x.nodes[4].typ === WHITESPACE
        ew += length(x.nodes[4])
        ew += length(x.nodes[5])
    end
    nest!(x.nodes[1], s, extra_width = extra_width + ew)
    nest!(x.nodes[2:end], s, x.indent, extra_width = extra_width)
end

# Import,Using,Export,ImportAll
function n_import!(x, s; extra_width = 0)
    line_margin = s.line_offset + length(x) + extra_width
    idx = findfirst(n -> n.typ === PLACEHOLDER, x.nodes)
    if idx !== nothing && (line_margin > s.margin || x.force_nest)
        sidx = findfirst(is_colon, x.nodes)
        if sidx === nothing
            sidx = 2
        else
            # for the WHITESPACE node after the colon
            sidx += 1
        end
        x.indent = s.line_offset + sum(length.(x.nodes[1:sidx]))
        s.line_offset = x.indent

        for (i, n) in enumerate(x.nodes)
            if n.typ === NEWLINE
                s.line_offset = x.indent
            elseif n.typ === PLACEHOLDER
                x.nodes[i] = Newline()
                s.line_offset = x.indent
            else
                nest!(n, s)
            end
        end
    else
        nest!(x.nodes, s, x.indent, extra_width = extra_width)
    end
end

function n_tuple!(x, s; extra_width = 0)
    line_margin = s.line_offset + length(x) + extra_width
    idx = findlast(n -> n.typ === PLACEHOLDER, x.nodes)
    # @info "ENTERING" idx x.typ s.line_offset length(x) extra_width
    opener = length(x.nodes) > 0 ? is_opener(x.nodes[1]) : false
    if idx !== nothing && (line_margin > s.margin || x.force_nest)
        if opener
            x.nodes[end].indent = x.indent
        end
        line_offset = s.line_offset

        x.indent += s.indent_size
        if x.indent - s.line_offset > 1
            x.indent = s.line_offset
            opener && (x.indent += 1)
        end

        # @debug "DURING" x.indent s.line_offset x.typ
        for (i, n) in enumerate(x.nodes)
            if n.typ === NEWLINE
                s.line_offset = x.indent
            elseif n.typ === PLACEHOLDER
                x.nodes[i] = Newline()
                s.line_offset = x.indent
            elseif n.typ === TRAILINGCOMMA
                x.nodes[i].val = ","
                x.nodes[i].len = 1
                nest!(x.nodes[i], s)
            elseif n.typ === TRAILINGSEMICOLON
                x.nodes[i].val = ""
                x.nodes[i].len = 0
                nest!(x.nodes[i], s)
            elseif opener && (i == 1 || i == length(x.nodes))
                nest!(n, s)
            else
                diff = x.indent - x.nodes[i].indent
                add_indent!(n, s, diff)
                nest!(n, s, extra_width = 1)
            end
        end

        if opener
            s.line_offset = x.nodes[end].indent + 1
        end
        # @info "EXITING" x.typ s.line_offset x.indent x.nodes[end].indent
    else
        opener && (extra_width += 1)
        nest!(x.nodes, s, x.indent, extra_width = extra_width)
    end
end

function n_for!(x, s; extra_width = 0)
    block_idx = findfirst(n -> !is_leaf(n), x.nodes)
    if block_idx === nothing
        nest!(x.nodes, s, x.indent, extra_width = extra_width)
        return
    end
    ph_idx = findfirst(n -> n.typ === PLACEHOLDER, x.nodes[block_idx].nodes)
    nest!(x.nodes, s, x.indent, extra_width = extra_width)

    # return if the argument block was nested
    ph_idx !== nothing && x.nodes[3].nodes[ph_idx].typ === NEWLINE && return

    idx = 5
    n = x.nodes[idx]
    if n.typ === NOTCODE && n.startline == n.endline
        res = get(s.doc.comments, n.startline, (0, ""))
        res == (0, "") && deleteat!(x.nodes, idx - 1)
    end
end

function n_call!(x, s; extra_width = 0)
    line_margin = s.line_offset + length(x) + extra_width
    idx = findlast(n -> n.typ === PLACEHOLDER, x.nodes)
    # @info "ENTERING" x.typ s.line_offset length(x) extra_width s.margin x.nodes[1].val
    if idx !== nothing && (line_margin > s.margin || x.force_nest)
        x.nodes[end].indent = x.indent
        line_offset = s.line_offset

        caller_len = length(x.nodes[1])

        x.indent += s.indent_size
        # @debug "ENTERING" x.indent s.line_offset x.typ
        if x.indent - s.line_offset > caller_len + 1
            x.indent = s.line_offset + caller_len + 1
            x.nodes[end].indent = s.line_offset
        end

        # @debug "DURING" x.indent s.line_offset x.typ
        for (i, n) in enumerate(x.nodes)
            if n.typ === NEWLINE
                s.line_offset = x.indent
            elseif n.typ === PLACEHOLDER
                x.nodes[i] = Newline()
                s.line_offset = x.indent
            elseif n.typ === TRAILINGCOMMA
                x.nodes[i].val = ","
                x.nodes[i].len = 1
                nest!(x.nodes[i], s)
            elseif n.typ === TRAILINGSEMICOLON
                x.nodes[i].val = ""
                x.nodes[i].len = 0
                nest!(x.nodes[i], s)
            elseif i == 1 || i == length(x.nodes)
                n.typ === CSTParser.Parameters && (n.force_nest = true)
                nest!(n, s, extra_width = 1)
            else
                diff = x.indent - x.nodes[i].indent
                add_indent!(n, s, diff)
                n.typ === CSTParser.Parameters && (n.force_nest = true)
                nest!(n, s, extra_width = 1)
            end
        end

        s.line_offset = x.nodes[end].indent + 1
    else
        is_closer(x.nodes[end]) && (extra_width += 1)
        nest!(x.nodes, s, x.indent, extra_width = extra_width)
    end
end

# "A where B"
function n_wherecall!(x, s; extra_width = 0)
    line_margin = s.line_offset + length(x) + extra_width
    # @info "" s.line_offset x.typ line_margin extra_width length(x) x.force_nest line_margin > s.margin
    if line_margin > s.margin || x.force_nest
        line_offset = s.line_offset
        # after "A where "
        idx = findfirst(n -> n.typ === PLACEHOLDER, x.nodes)
        Blen = sum(length.(x.nodes[idx+1:end]))

        # A, +7 is the length of " where "
        ew = Blen + 7 + extra_width
        nest!(x.nodes[1], s, extra_width = ew)

        for (i, n) in enumerate(x.nodes[2:idx-1])
            nest!(n, s)
        end

        # "B"

        has_braces = is_closer(x.nodes[end])
        if has_braces
            x.nodes[end].indent = x.indent
        end

        over = (s.line_offset + Blen + extra_width > s.margin) || x.force_nest
        x.indent += s.indent_size

        last_typ::Union{CSTParser.Head,Nothing} = nothing
        for (i, n) in enumerate(x.nodes[idx+1:end])
            if n.typ === NEWLINE
                s.line_offset = x.indent
            elseif is_opener(n)
                if x.indent - s.line_offset > 1
                    x.indent = s.line_offset + 1
                    x.nodes[end].indent = s.line_offset
                end
                nest!(n, s)
            elseif n.typ === PLACEHOLDER && over
                x.nodes[i+idx] = Newline()
                s.line_offset = x.indent
            elseif n.typ === TRAILINGCOMMA && over
                x.nodes[i+idx].val = ","
                x.nodes[i+idx].len = 1
                nest!(x.nodes[i+idx], s)
            elseif has_braces
                nest!(n, s, extra_width = 1 + extra_width)
            else
                nest!(n, s, extra_width = extra_width)
            end
            (!is_leaf(n) || n.typ === CSTParser.IDENTIFIER) && (last_typ = n.typ)
        end

        s.line_offset = line_offset
        walk(reset_line_offset!, x, s)
    else
        nest!(x.nodes, s, x.indent, extra_width = extra_width)
    end
end

function n_condcall!(x, s; extra_width = 0)
    line_margin = s.line_offset + length(x) + extra_width
    # @info "ENTERING" x.typ s.line_offset line_margin x.force_nest
    if line_margin > s.margin || x.force_nest
        line_offset = s.line_offset
        x.indent = s.line_offset
        phs = reverse(findall(n -> n.typ === PLACEHOLDER, x.nodes))
        for (i, idx) in enumerate(phs)
            if i == 1
                x.nodes[idx] = Newline()
            else
                nidx = phs[i-1]
                l1 = sum(length.(x.nodes[1:idx-1]))
                l2 = sum(length.(x.nodes[idx:nidx-1]))
                width = line_offset + l1 + l2
                if x.force_nest || width > s.margin
                    x.nodes[idx] = Newline()
                end
            end
        end

        for (i, n) in enumerate(x.nodes)
            if i == length(x.nodes)
                nest!(n, s, extra_width = extra_width)
            elseif x.nodes[i+1].typ === WHITESPACE
                ew = length(x.nodes[i+1]) + length(x.nodes[i+2])
                nest!(n, s, extra_width = ew)
            elseif n.typ === NEWLINE
                s.line_offset = x.indent
            else
                nest!(n, s)
            end
        end

        s.line_offset = line_offset
        for (i, n) in enumerate(x.nodes)
            if n.typ === NEWLINE && !is_comment(x.nodes[i+1]) && !is_comment(x.nodes[i-1])
                # +1 for newline to whitespace conversion
                width = s.line_offset + 1
                if i == length(x.nodes) - 1
                    width += sum(length.(x.nodes[i+1:end])) + extra_width
                else
                    width += sum(length.(x.nodes[i+1:i+3]))
                end
                # @debug "" s.line_offset l  s.margin
                if width <= s.margin
                    x.nodes[i] = Whitespace(1)
                else
                    s.line_offset = x.indent
                end
            end
            s.line_offset += length(x.nodes[i])
        end

        s.line_offset = line_offset
        walk(reset_line_offset!, x, s)
    else
        nest!(x.nodes, s, x.indent, extra_width = extra_width)
    end
end
n_comparison!(x, s; extra_width = 0) = n_block!(x, s, extra_width = extra_width)
n_chainopcall!(x, s; extra_width = 0) = n_block!(x, s, extra_width = extra_width)

function dedent!(x::PTree, s::State, pindent::Int, line_margin::Int)
    is_leaf(x) && return
    x.typ === CSTParser.ConditionalOpCall && return
    x.typ === CSTParser.Comparison && return
    x.typ === CSTParser.ChainOpCall && return
    x.typ === CSTParser.BinaryOpCall && return

    if is_iterable(x)
        indent = x.indent
        args = get_args(x)
        while length(args) == 1 && is_iterable(args[1])
            indent = args[1].indent
            args = get_args(args[1])
        end
        diff = min(line_margin, pindent + s.indent_size) - indent
        # @info "" pindent diff x.indent line_margin indent
        # @info "" pindent diff x.indent x.typ
        add_indent!(x, s, diff)
        x.nodes[end].indent = pindent
    else
        add_indent!(x, s, -s.indent_size)
    end
end

# arg1 op arg2
#
# nest in order of
#
# arg1 op
# arg2
function n_binarycall!(x, s; extra_width = 0)
    # If there's no placeholder the binary call is not nestable
    idxs = findall(n -> n.typ === PLACEHOLDER, x.nodes)
    line_margin = s.line_offset + length(x) + extra_width
    # @info "ENTERING" x.typ extra_width s.line_offset length(x) idxs x.ref[][2]
    if length(idxs) == 2 && (line_margin > s.margin || x.force_nest)
        line_offset = s.line_offset
        i1 = idxs[1]
        i2 = idxs[2]
        x.nodes[i1] = Newline()
        cst = x.ref[]

        has_eq = CSTParser.defines_function(cst) || nest_assignment(cst)
        if has_eq
            s.line_offset = x.indent + s.indent_size
            x.nodes[i2] = Whitespace(s.indent_size)
            add_indent!(x.nodes[end], s, s.indent_size)
        else
            x.indent = s.line_offset
        end

        # arg2
        nest!(x.nodes[end], s, extra_width = extra_width)

        # arg1 op
        s.line_offset = line_offset
        # extra_width for " op"
        ew = length(x.nodes[2]) + length(x.nodes[3])
        # @debug "DURING" ew s.line_offset x.typ
        nest!(x.nodes[1], s, extra_width = ew)
        for n in x.nodes[2:i1]
            nest!(n, s)
        end

        # Undo nest if possible
        if !x.force_nest
            arg2 = x.nodes[end]
            arg2.typ === CSTParser.Block && (arg2 = arg2.nodes[1])
            cst = arg2.ref[]

            if (
                arg2.typ === CSTParser.BinaryOpCall && (
                    !(is_lazy_op(cst) && !has_eq) && cst[2].kind !== Tokens.IN
                )
            ) || arg2.typ === CSTParser.UnaryOpCall || is_block(cst)
                line_margin = s.line_offset + 1 + length(x.nodes[end])
            else
                rw, _ = length_to(x, [NEWLINE], start = i2 + 1)
                line_margin = s.line_offset + 1 + rw
            end
            # @info "" can_unnest arg2.typ has_eq

            if line_margin + extra_width <= s.margin
                x.nodes[i1] = Whitespace(1)
                if has_eq
                    x.nodes[i2] = Placeholder(0)
                    # x.indent -= s.indent_size
                    # recursive dedent
                    if arg2.typ === CSTParser.BinaryOpCall
                        dedent!(arg2.nodes[1], s, x.indent, line_margin)
                        dedent!(arg2.nodes[end], s, x.indent, line_margin)
                    else
                        dedent!(arg2, s, x.indent, line_margin)
                    end
                end
            end
        end

        # @info "before reset" line_offset s.line_offset

        s.line_offset = line_offset
        walk(reset_line_offset!, x, s)

        # @info "after reset" line_offset s.line_offset
    else
        # Handles the case of a function def defined
        # as "foo(a)::R where {A,B} = body".
        #
        # In this case instead of it being parsed as:
        #
        # CSTParser.BinaryOpCall
        #  - CSTParser.WhereOpCall
        #  - OP
        #  - arg2
        #
        # It's parsed as:
        #
        # CSTParser.BinaryOpCall
        #  - CSTParser.BinaryOpCall
        #   - arg1
        #   - OP
        #   - CSTParser.WhereOpCall
        #    - R
        #    - ...
        #  - OP
        #  - arg2
        #
        # The result being extra width is trickier to deal with.
        idx = findfirst(n -> n.typ === CSTParser.WhereOpCall, x.nodes)
        return_width = 0
        if idx !== nothing && idx > 1
            return_width = length(x.nodes[idx].nodes[1]) + length(x.nodes[2])
        elseif idx === nothing
            return_width, _ = length_to(x, [NEWLINE], start = 2)
        end

        # @info "" return_width
        # @info "" s.line_offset return_width length(x.nodes[1])

        for (i, n) in enumerate(x.nodes)
            if n.typ === NEWLINE
                s.line_offset = x.indent
            elseif i == 1
                nest!(n, s, extra_width = return_width)
            elseif i == idx
                nest!(n, s, extra_width = extra_width)
            else
                nest!(n, s, extra_width = extra_width)
            end
        end
    end
end

function n_block!(x, s; extra_width = 0, custom_indent = 0)
    line_margin = s.line_offset + length(x) + extra_width
    idx = findfirst(n -> n.typ === PLACEHOLDER, x.nodes)
    # @info "ENTERING" idx x.typ s.line_offset length(x) extra_width
    if idx !== nothing && (line_margin > s.margin || x.force_nest)
        line_offset = s.line_offset

        x.indent = custom_indent > 0 ? custom_indent : s.line_offset

        # @info "DURING" x.indent s.line_offset x.typ
        for (i, n) in enumerate(x.nodes)
            if n.typ === NEWLINE
                s.line_offset = x.indent
            elseif n.typ === PLACEHOLDER
                x.nodes[i] = Newline()
                s.line_offset = x.indent
            elseif n.typ === TRAILINGCOMMA
                x.nodes[i].val = ","
                x.nodes[i].len = 1
                nest!(x.nodes[i], s)
            elseif n.typ === TRAILINGSEMICOLON
                x.nodes[i].val = ""
                x.nodes[i].len = 0
                nest!(x.nodes[i], s)
            else
                diff = x.indent - x.nodes[i].indent
                add_indent!(n, s, diff)
                nest!(n, s, extra_width = 1)
            end
        end

        # @info "EXITING" x.typ s.line_offset x.indent x.nodes[end].indent
    else
        nest!(x.nodes, s, x.indent, extra_width = extra_width)
    end
end

n_gen!(x, s; extra_width = 0) =
    n_block!(x, s, extra_width = extra_width, custom_indent = x.indent)


text:
function K_MITC4(elm::Quad)::Matrix{Float64}
    E₀,ν₀=elm.material.E,elm.material.ν
    center=elm.center
    t=elm.t
    T=elm.T[1:3,1:3]
    x₁,y₁,z₁=T*(elm.node1.loc-center)
    x₂,y₂,z₂=T*(elm.node2.loc-center)
    x₃,y₃,z₃=T*(elm.node3.loc-center)
    x₄,y₄,z₄=T*(elm.node4.loc-center)
    K=Matrix{Float64}(undef,12,12)
    J=Matrix{Float64}(undef,2,2)
    return K
end


text:
using LAP_julia

viewimg(img) = colorview(Gray, img)

img_size = (250,255)

view(img)
view(imgw)

timer = TimerOutput(string(args_dict["reg_alg"]))
method_kwargs = Dict(:display => true, :timer => timer)
img, imgw, flow = gen_init(:spaghetti, img_args=[img_size], flow_args=[25])
flow_est, source_reg, timer, time_in_secs = test_registration_alg(sparse_pflap,
                                                                  img,
                                                                  imgw,
                                                                  flow,
                                                                  timer=method_kwargs[:timer],
                                                                  method_kwargs=method_kwargs)

showflow(flow_est)
showflow(flow.*(-1))

count(isnan.([real.(flow_est) imag.(flow_est)]))

LAP_julia.vec_len.(flow_est)
view(source_reg)
view(img)
view(imgw)

##
base_path_bikes = "/Users/MrTrololord/Downloads/bikes"

img = load(joinpath(base_path_bikes, "img1.png"))
for k in 2:2
    imgw = load(joinpath(base_path_bikes, "img$(k).png"))
    global imgw
end

imgw

img = Float64.(Gray.(img))
imgw = Float64.(Gray.(imgw))

timer = TimerOutput("sparse_pflap");
method_kwargs = Dict(:display => false, :timer => timer)
flow_est, source_reg, timer, time_in_secs = time_reg_alg(pflap, img, imgw, timer=method_kwargs[:timer], method_kwargs=method_kwargs)

view(imgw)
view(source_reg)
view(img)

showflow(flow_est)

##

base_path_zhang = "/Users/MrTrololord/Downloads/dataset/"
fnames = ["yard", "PGH", "testchart"]
imgs = Array{Image}(undef, size(fnames))
imgws = Array{Image}(undef, size(fnames))

for (k, name) in zip(1:length(fnames), fnames)
    imgs[k] = load(joinpath(base_path_zhang, "$(name)_t.jpg"))
    imgws[k] = load(joinpath(base_path_zhang, "$(name)_s.jpg"))
    global img, imgw
end

img = imgs[3]
imgw = imgws[3]

img = Float64.(Gray.(img))
imgw = Float64.(Gray.(imgw))

##
timer = TimerOutput("pflap");
method_kwargs = Dict(:display => false, :timer => timer)
flow_est, source_reg, timer, time_in_secs = time_reg_alg(pflap, img, imgw, timer=method_kwargs[:timer], method_kwargs=method_kwargs)

source_reg_zeros = warp_img(imgw, -real.(flow_est), -imag.(flow_est), border_strat=:zeros)
view(source_reg_zeros)

view(imgw)
view(source_reg)
viewimg(img)

showflow(imresize(flow_est, (200,200)))

##
# @info "size from $(size)"
(img, ratio), (imgw, ratio) = resize_to_diag_size(img, 400), resize_to_diag_size(imgw, 400)

img, imgw = Gray.(img), Gray.(imgw)

timer = TimerOutput("sparse_pflap");
method_kwargs = Dict(:display => false, :timer => timer)
flow_est, source_reg, timer, time_in_secs = time_reg_alg(sparse_pflap_psnr, imgw, img, timer=method_kwargs[:timer], method_kwargs=method_kwargs)

source_reg_zeros = warp_img(img, -real.(flow_est), -imag.(flow_est), border_strat=:zeros)
viewimg(source_reg_zeros)
viewimg(imgw)

flow_est

##
# load float from text

function parse_numbers(s)
    pieces = split(s, ' ', keepempty=false)
    map(pieces) do piece
        parse(Float64, piece)
    end
end

function load_H(path)
    lines = readlines(path)
    H = Array{Float64}(undef, (3,3))
    for (k, line) in zip(1:3, lines)
        H[k,:] = parse_numbers(line)
    end
    return H
end

function make_flow_from_H(H, size)
    warp_perspective(y, x) = (H[1,1]*x + H[1,2]*y + H[1,3])/(H[3,1]*x + H[3,2]*y + H[3,3]),
                             (H[2,1]*x + H[2,2]*y + H[2,3])/(H[3,1]*x + H[3,2]*y + H[3,3])
    flow = Array{Complex{Float64},2}(undef, size)
    X = ones(size[1]) * collect(range(0,1,length=size[2]))'
    Y = collect(range(0,1,length=size[1])) * ones(size[2])'
    for ind in CartesianIndices(flow)
        calc = [Y[ind], X[ind]] .- warp_perspective(Y[ind], X[ind])
        flow[ind] = calc[1] + im*calc[2]
    end
    return flow
end

H = load_H(joinpath(base_path_bikes, "H1to2p"))
out = make_flow_from_H(H, size(img))

showflow(out)

img
imgw

asdf = warp_img(img, -imag.(out), -real.(out), border_strat=:zeros)

Gray.(img) .- Gray.(asdf)

##
using CoordinateTransformations, StaticArrays, ImageTransformations, LinearAlgebra;
struct Homography{T} <: AbstractAffineMap
    m::SMatrix{3, 3, T, 9}
end

# having T in the RHS forces calling the constructor again
Homography(m::AbstractMatrix{}) = Homography(m);
function (trans::Homography{M})(x::SVector{3}) where M <: Real
    out = trans.m * x;
    out = out / out[end];
    SVector{2}(out[1:2])
end

function (trans::Homography{M})(x::SVector{2}) where M <: Real
    trans(SVector{3}([x[1], x[2], 1.0]))
end

function (trans::Homography{M})(x::CartesianIndex{2}) where M <: Real
    trans(SVector{3}([collect(x.I); 1]))
end

function (trans::Homography{M})(x::Tuple{Int, Int}) where M <: Real
    trans(CartesianIndex{2}(x))
end

function (trans::Homography{M})(x::Array{CartesianIndex{2}, 1}) where M <: Real
    CartesianIndex{2}.([tuple(y...) for y in trunc.(Int, collect.(trans.(x)))])
end

function Base.inv(trans::Homography)
    i = inv(trans.m);
    Homography(i ./ i[end])
end

Homo = Homography{Float64}(H)

img1_warp = ImageTransformations.warp(imgw, Homo)

img

##

## posun homografii

base_path_bikes = "/Users/MrTrololord/Downloads/bikes"

img = load(joinpath(base_path_bikes, "img1.png"))
k = 2
imgw = load(joinpath(base_path_bikes, "img$(k).png"))

function parse_numbers(s)
    pieces = split(s, ' ', keepempty=false)
    map(pieces) do piece
        parse(Float64, piece)
    end
end

function load_H(path)
    lines = readlines(path)
    H = Array{Float64}(undef, (3,3))
    for (k, line) in zip(1:3, lines)
        H[k,:] = parse_numbers(line)
    end
    return H
end

function make_flow_from_H(H, size)
    warp_perspective(y, x) = (H[1,1]*x + H[1,2]*y + H[1,3])/(H[3,1]*x + H[3,2]*y + H[3,3]),
                             (H[2,1]*x + H[2,2]*y + H[2,3])/(H[3,1]*x + H[3,2]*y + H[3,3])
    flow = Array{Complex{Float64},2}(undef, size)
    X = ones(size[1]) * collect(range(0,1,length=size[2]))'
    Y = collect(range(0,1,length=size[1])) * ones(size[2])'
    for ind in CartesianIndices(flow)
        calc = [Y[ind], X[ind]] .- warp_perspective(Y[ind], X[ind])
        flow[ind] = calc[1] + im*calc[2]
    end
    return flow
end

H = load_H(joinpath(base_path_bikes, "H1to2p"))
out = make_flow_from_H(H, size(img))

warped_img = warp_img(img, -imag.(out), -real.(out), border_strat=:zeros)
imgw

## posun mym alg

viewimg(img) = colorview(Gray, img)

img, imgw = Gray.(img), Gray.(imgw)

timer = TimerOutput("sparse_pflap");
method_kwargs = Dict(:display => false, :timer => timer)
flow_est, source_reg, timer, time_in_secs = time_reg_alg(sparse_pflap_psnr, imgw, img, timer=method_kwargs[:timer], method_kwargs=method_kwargs)

source_reg_zeros = warp_img(img, -real.(flow_est), -imag.(flow_est), border_strat=:zeros)
viewimg(source_reg_zeros)
viewimg(imgw)

# posun mym alg je ma lepsi vysledky?? nesmysl -> neco spatne


text:
using Turing, ContinuousBenchmarks

include(joinpath(ContinuousBenchmarks.STAN_DATA_DIR, "negative_binomial.data.jl"))
include(joinpath(ContinuousBenchmarks.STAN_MODELS_DIR, "negative_binomial.model.jl"))

# Produce an error.
sample(negbinmodel(negbindata), HMC(1000, 0.02, 1));


text:
#  _____ _____ ____ _____ ____        ____        ____  ___   ___  ____
# |_   _| ____/ ___|_   _/ ___|      |  _ \      / ___|/ _ \ / _ \|  _ \
#   | | |  _| \___ \ | | \___ \ _____| |_) |____| |  _| | | | | | | | | |
#   | | | |___ ___) || |  ___) |_____|  _ <_____| |_| | |_| | |_| | |_| |
#   |_| |_____|____/ |_| |____/      |_| \_\     \____|\___/ \___/|____/

# NOTE TO READERS:
# This is a test case, so we'll be exercising advanced features
# provided by DynamoDB. If you're looking for simple examples to start
# from, please check out example.jl in the doc/ folder of this repo.

# see runtests.jl for Foo's definition

foo_basic = dynamo_table(Foo, "foo_basic", :a)
foo_basic_lsi = dynamo_local_index(foo_basic, "foo_basic_indexed_on_b", :b)
foo_basic_gsi = dynamo_global_index(foo_basic, "foo_basic_global_index_on_b_a", :b, :a)
foo_basic_gsi_no_idx = dynamo_global_index(foo_basic, "foo_basic_global_index_on_b", :b)

foo_range = dynamo_table(Foo, "foo_range", :a, :b)
foo_range_gsi = dynamo_global_index(foo_range, "foo_range_global_index_on_b_a", :b, :a)

# keydict in all its glory
@test_throws ErrorException DynamoDB.keydict(foo_basic, 1, 2)
@test DynamoDB.keydict(foo_basic, 1) == Dict("a" => Dict("N" => "1"))

@test DynamoDB.keydict(foo_basic_lsi, 1, 2) == Dict("a" => Dict("N" => "1"), "b" => Dict("N" => "2"))
@test_throws ErrorException DynamoDB.keydict(foo_basic_lsi, 1)

@test DynamoDB.keydict(foo_basic_gsi, 1, 2) == Dict("b" => Dict("N" => "1"), "a" => Dict("N" => "2"))
@test DynamoDB.keydict(foo_basic_gsi, 1) == Dict("b" => Dict("N" => "1"), "a" => Dict("NULL" => true))

@test_throws ErrorException DynamoDB.keydict(foo_basic_gsi_no_idx, 1, 2)
@test DynamoDB.keydict(foo_basic_gsi, 1) == Dict("b" => Dict("N" => "1"), "a" => Dict("NULL" => true))

@test DynamoDB.keydict(foo_range, 1, 2) == Dict("a" => Dict("N" => "1"), "b" => Dict("N" => "2"))
@test DynamoDB.keydict(foo_range, 1) == Dict("a" => Dict("N" => "1"), "b" => Dict("NULL" => true))

@test DynamoDB.keydict(foo_range_gsi, 1, 2) == Dict("b" => Dict("N" => "1"), "a" => Dict("N" => "2"))
@test DynamoDB.keydict(foo_range_gsi, 1) == Dict("b" => Dict("N" => "1"), "a" => Dict("NULL" => true))

@test_throws ErrorException DynamoDB.check_status(404, "test")
@test DynamoDB.check_status(200, "test") == nothing


text:
# ---
# title: 97. Interleaving String
# id: problem97
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: String, Dynamic Programming
# link: <https://leetcode.com/problems/interleaving-string/description/>
# hidden: true
# ---
#
# Given strings `s1`, `s2`, and `s3`, find whether `s3` is formed by an
# **interleaving** of `s1` and `s2`.
#
# An **interleaving** of two strings `s` and `t` is a configuration where they
# are divided into **non-empty** substrings such that:
#
#   * `s = s1 + s2 + ... + sn`
#   * `t = t1 + t2 + ... + tm`
#   * `|n - m| <= 1`
#   * The **interleaving** is `s1 + t1 + s2 + t2 + s3 + t3 + ...` or `t1 + s1 + t2 + s2 + t3 + s3 + ...`
#
# **Note:** `a + b` is the concatenation of strings `a` and `b`.
#
#
#
# **Example 1:**
#
# ![](https://assets.leetcode.com/uploads/2020/09/02/interleave.jpg)
#
#
#
#     Input: s1 = "aabcc", s2 = "dbbca", s3 = "aadbbcbcac"
#     Output: true
#
#
# **Example 2:**
#
#
#
#     Input: s1 = "aabcc", s2 = "dbbca", s3 = "aadbbbaccc"
#     Output: false
#
#
# **Example 3:**
#
#
#
#     Input: s1 = "", s2 = "", s3 = ""
#     Output: true
#
#
#
#
# **Constraints:**
#
#   * `0 <= s1.length, s2.length <= 100`
#   * `0 <= s3.length <= 200`
#   * `s1`, `s2`, and `s3` consist of lower-case English letters.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
export joLagrangeInterp1D

# 1D Cubic Lagrange Interpolation
#
# Parameters:
#   xin  - input grid
#   xout - output grid
#   T    - vector data type
#
# Note: the interval [minimum(xout),maximum(xout)] must be contained
# in the interval [minimum(xin),maximum(xin)]

function joLagrangeInterp1D(xin,xout,T)
    ik = find((xout .>= xin[1]) .& (xout .<= xin[end]))
    nout = length(xout)
    nin = length(xin)
    nk = length(ik)
    nk > 0 || error("The set generated by xout must be contained within the interval [xin[1],xin[end]]")
    I = zeros(Int64,4*nk)
    J = zeros(Int64,4*nk)
    S = zeros(T,4*nk)
    a,b,c,d = 1,2,3,4
    l=1
    for i = 1:nk
        k = ik[i];
        if xout[k]<xin[b]
            while (xout[k]<xin[b])&&(b-1>1)
                b=b-1;
            end
            a=b-1;c=b+1;d=c+1;
        elseif xout[k]>xin[c]
            while (xout[k]>xin[c])&&(c+1<nin)
                c=c+1;
            end
            a=c-2;b=c-1;d=c+1;
        end
        I[l:l+3] = k;
        J[l]   = a;
        S[l]   = ((xout[k]-xin[b])*(xout[k]-xin[c])*(xout[k]-xin[d]))/((xin[a]-xin[b])*(xin[a]-xin[c])*(xin[a]-xin[d]));

        J[l+1] = b;
        S[l+1] = ((xout[k]-xin[a])*(xout[k]-xin[c])*(xout[k]-xin[d]))/((xin[b]-xin[a])*(xin[b]-xin[c])*(xin[b]-xin[d]));

        J[l+2] = c;
        S[l+2] = ((xout[k]-xin[b])*(xout[k]-xin[a])*(xout[k]-xin[d]))/((xin[c]-xin[b])*(xin[c]-xin[a])*(xin[c]-xin[d]));

        J[l+3] = d;
        S[l+3] = ((xout[k]-xin[b])*(xout[k]-xin[c])*(xout[k]-xin[a]))/((xin[d]-xin[b])*(xin[d]-xin[c])*(xin[d]-xin[a]));

        l += 4;
    end

    return joMatrix(sparse(I[1:l-1],J[1:l-1],S[1:l-1],nout,nin),DDT=T,name="joLagrangeInterp1D")
end


text:
# initialize world
function initialize_world(dim, pop_density)
    world = rand(dim, dim)
    # distribute agents probabilistically
    for n in 1:dim, m in 1:dim
        if world[n, m] < pop_density
            if rand() < 0.5
                world[n, m] = 1.0
            else
                world[n, m] = 2.0
            end
        else
            world[n, m] = 0.0
        end
    end
    return world
end

# is agent happy?
function agent_happy(world, i, j, threshold)
    similar = 0
    different = 0
    positions = [
        (i-1, j-1), (i-1, j  ), (i-1, j+1),
        (i,   j-1),             (i,   j+1),
        (i+1, j-1), (i+1, j  ), (i+1, j+1)
    ]
    for pos in positions
        # case 1: row out of bounds (0)
        if pos[1] == 0
            # case 1a: column out of bounds (0)
            if pos[2] == 0
                world[size(world)[1], size(world)[1]] == world[i, j] ? similar += 1 : different += 1
            # case 1b: column out of bounds (dim)
            elseif pos[2] == (size(world)[1] + 1)
                world[size(world)[1], 1] == world[i, j] ? similar += 1 : different += 1
            # case 1c: column not out of bounds
            else
                world[size(world)[1], pos[2]] == world[i, j] ? similar += 1 : different += 1
            end
        # case 2: row out of bounds (dim)
        elseif pos[1] == (size(world)[1] + 1)
            # case 2a: column out of bounds (0)
            if pos[2] == 0
                world[1, size(world)[1]] == world[i, j] ? similar += 1 : different += 1
            # case 2b: column oht of bounds (dim)
            elseif pos[2] == (size(world)[1] + 1)
                world[1, 1] == world[i, j] ? similar += 1 : different += 1
            # case 2c: column not out of bounds
            else
                world[1, pos[2]] == world[i, j] ? similar += 1 : different += 1
            end
        # case 3: row not out of bounds
        else
            # case 3a: column out of bounds (0)
            if pos[2] == 0
                world[pos[1], size(world)[1]] == world[i, j] ? similar += 1 : different += 1
            # case 3b: column out of bounds (dim)
            elseif pos[2] == (size(world)[1] + 1)
                world[pos[1], 1] == world[i, j] ? similar += 1 : different += 1
            # case 3c: column not out of bounds
            else
                world[pos[1], pos[2]] == world[i, j] ? similar += 1 : different += 1
            end
        end
    end
    # return true if relative count of similar agents in neighborhood is > threshold, else false
    return similar/different >= threshold ? true : false
end

# find a new place to live
function relocate(world, i, j)
    new_row = convert(Int64, ceil(rand() * size(world)[1]))
    new_col = convert(Int64, ceil(rand() * size(world)[1]))
    if world[new_row, new_col] == 0.0
        world[new_row, new_col] = world[i, j]
        world[i, j] = 0.0
    else
        relocate(world, i, j)
    end
    return world
end

# simulation step
function tick(world, threshold)
    for n in 1:size(world)[1], m in 1:size(world)[1]
        if !agent_happy(world, n, m, threshold)
            relocate(world, n, m)
        end
    end
    return world
end

# run the simulation for max_ticks steps at given threshold
function simulate(world, max_ticks, threshold)
    for t in 1:max_ticks
        world = tick(world, threshold)
    end
    return world
end

max_stars_repo_name: JohannesNakayama/schelling-model

max_issues_repo_name: JohannesNakayama/schelling-model

max_forks_repo_name: JohannesNakayama/schelling-model


text:
module Web

export Backend

include("backend/backend.jl")
using Web.Backend

end # module

max_stars_repo_name: Gollor/JuliaAutomation

max_issues_repo_name: Gollor/JuliaAutomation

max_forks_repo_name: Gollor/JuliaAutomation


text:
# Some utilities to assist in defining BNF grammars.

which_BNF_grammar = nothing

deferred_bnf_strs = []

function do_bnf_str(str::String, grammar_name::Symbol)
    if !haskey(AllGrammars, grammar_name)
        BNFGrammar(grammar_name)
    end
    g::BNFGrammar = AllGrammars[Symbol(grammar_name)]
    bnf = AllGrammars[which_BNF_grammar]
    recognize(bnf["<syntax>"], str)
end

"""
Parse `str` as BNF and add those productions to the grammar named `grammar_name`.
"""
macro bnf_str(str, grammar_name)
    # HOW TO CAPTURE SOURCE LOCATION? SEE cl_str.
    grammar_name = Symbol(grammar_name)
    if which_BNF_grammar === nothing
        push!(deferred_bnf_strs, (str, grammar_name))
        return str
    end
    do_bnf_str(str, grammar_name)
end


text:
using JWebImageDemo
using Test

import HTTP

const ON_HEROKU = "true"
const PORT = rand(8090:9999)
const HOST = "127.0.0.1"
const VERBOSE = 2

const TEST_IMAGE = joinpath(@__DIR__, "resources", "green_python.jpg")
const NUM = 5

const TEST_REQUESTS = Dict(
    "f_image" => HTTP.Multipart("test.jpg", open(TEST_IMAGE)),
    "num" => HTTP.Multipart(nothing, IOBuffer("5"), "plain/text")
)

@info "Running web service"
service = @async JWebImageDemo.AppServer.start_server(;host = string(HOST), port = PORT)

sleep(10)
@testset "web service" begin
   for i = 1:10
        try
            r = HTTP.get("http://$(HOST):$(PORT)/"; verbose = VERBOSE)
            # dump(r)
            @test r.status == 200
            @test contains(String(r.body), "form")

            r = HTTP.post(
                "http://$(HOST):$(PORT)/api/process_image/",
                HTTP.Form(TEST_REQUESTS);
                verbose = VERBOSE,
            )
            # dump(r)
            @info "Status: $(r.status)"
            @test r.status == 200
            @test any(t -> contains(t.first, "Content-Type") &&
                           contains(t.second, "image/png"),
                      r.headers)

           break
        catch
        end
        sleep(30)
   end
end

ex = InterruptException()
@async Base.throwto(service, ex)
sleep(5)


text:
"""
`Basics.provide()`  ... provides various data and transformations.

  + `("subshell states: antisymmetric, seniority", sh::Subshell, occ::Int64)`
    ... to provide all antisymmetric SubshellStates within the relativistic seniority scheme for the given
        shell and occupation; an Array{SubshellStateR,1} is returned.
"""
function Basics.provide(sa::String, sh::Subshell, occ::Int64)

    !(sa == "subshell states: antisymmetric, seniority")   &&   error("Unsupported keystring = $sa")

    # Use the j-value = wa[4] of the subshell to select the allowed SubshellStateR's
    wb = Basics.SubshellStateR[]

    if       Basics.subshell_2j(sh) == 1
        if       occ == 0  || occ == 2      push!( wb, Basics.SubshellStateR( sh, occ, 0, 0) );       return(wb)
        elseif   occ == 1                   push!( wb, Basics.SubshellStateR( sh, occ, 1, 1) );       return(wb)
        end
    elseif   Basics.subshell_2j(sh) == 3
        if       occ == 0  || occ == 4      push!( wb, Basics.SubshellStateR( sh, occ, 0, 0) );       return(wb)
        elseif   occ == 1  || occ == 3      push!( wb, Basics.SubshellStateR( sh, occ, 1, 3) );       return(wb)
        elseif   occ == 2                   push!( wb, Basics.SubshellStateR( sh, occ, 0, 0) )
                                            push!( wb, Basics.SubshellStateR( sh, occ, 2, 4) );       return(wb)
        end
    elseif   Basics.subshell_2j(sh) == 5
        if       occ == 0  || occ == 6      push!( wb, Basics.SubshellStateR( sh, occ, 0, 0) );       return(wb)
        elseif   occ == 1  || occ == 5      push!( wb, Basics.SubshellStateR( sh, occ, 1, 5) );       return(wb)
        elseif   occ == 2  || occ == 4      push!( wb, Basics.SubshellStateR( sh, occ, 0, 0) )
                                            push!( wb, Basics.SubshellStateR( sh, occ, 2, 4) )
                                            push!( wb, Basics.SubshellStateR( sh, occ, 2, 8) );       return(wb)
        elseif   occ == 3                   push!( wb, Basics.SubshellStateR( sh, occ, 1, 5) )
                                            push!( wb, Basics.SubshellStateR( sh, occ, 3, 3) )
                                            push!( wb, Basics.SubshellStateR( sh, occ, 3, 9) );       return(wb)
        end
    elseif   Basics.subshell_2j(sh) == 7
        if       occ == 0  || occ == 8      push!( wb, Basics.SubshellStateR( sh, occ, 0, 0) );       return(wb)
        elseif   occ == 1  || occ == 7      push!( wb, Basics.SubshellStateR( sh, occ, 1, 7) );       return(wb)
        elseif   occ == 2  || occ == 6      push!( wb, Basics.SubshellStateR( sh, occ, 0, 0) )
                                            push!( wb, Basics.SubshellStateR( sh, occ, 2, 4) )
                                            push!( wb, Basics.SubshellStateR( sh, occ, 2, 8) )
                                            push!( wb, Basics.SubshellStateR( sh, occ, 2,12) );       return(wb)
        elseif   occ == 3  || occ == 5      push!( wb, Basics.SubshellStateR( sh, occ, 1, 7) )
                                            push!( wb, Basics.SubshellStateR( sh, occ, 3, 3) )
                                            push!( wb, Basics.SubshellStateR( sh, occ, 3, 5) )
                                            push!( wb, Basics.SubshellStateR( sh, occ, 3, 9) )
                                            push!( wb, Basics.SubshellStateR( sh, occ, 3,11) )
                                            push!( wb, Basics.SubshellStateR( sh, occ, 3,15) );       return(wb)
        elseif   occ == 4                   push!( wb, Basics.SubshellStateR( sh, occ, 0, 0) )
                                            push!( wb, Basics.SubshellStateR( sh, occ, 2, 4) )
                                            push!( wb, Basics.SubshellStateR( sh, occ, 2, 8) )
                                            push!( wb, Basics.SubshellStateR( sh, occ, 2,12) )
                                            push!( wb, Basics.SubshellStateR( sh, occ, 4, 4) )
                                            push!( wb, Basics.SubshellStateR( sh, occ, 4, 8) )
                                            push!( wb, Basics.SubshellStateR( sh, occ, 4,10) )
                                            push!( wb, Basics.SubshellStateR( sh, occ, 4,16) );       return(wb)
        end
    elseif   Basics.subshell_2j(sh) in [ 9, 11, 13, 15]
        if       occ == 0                   push!( wb, Basics.SubshellStateR( sh, occ, 0, 0) );                 return(wb)
        elseif   occ == 1                   push!( wb, Basics.SubshellStateR( sh, occ, 1, subshell_2j(sh)) );   return(wb)
        else     error("stop b")
        end

    # For an extension to other subshell-occupations, see the Racah-manual; Table ...
    else  error("stop a")
    end
end

"""
  + `("binding energy", Z::Float64, sh::Subshell)`
    ... to provide the binding energy of a subshell electron, taken from a semi-empirical tabulations by Williams et al.,
        https://userweb.jlab.org/~gwyn/ebindene.html. A energy::Float64 in  Hartree is returned.
"""
function Basics.provide(sa::String, Z::Int64, sh::Subshell)
    if     sa == "binding energy"
        wa = JAC.store_Williams2000(Z)
        if      sh == Subshell("1s_1/2")    wb = wa[1]
        elseif  sh == Subshell("2s_1/2")    wb = wa[2]
        elseif  sh == Subshell("2p_1/2")    wb = wa[3]
        elseif  sh == Subshell("2p_3/2")    wb = wa[4]
        elseif  sh == Subshell("3s_1/2")    wb = wa[5]
        elseif  sh == Subshell("3p_1/2")    wb = wa[6]
        elseif  sh == Subshell("3p_3/2")    wb = wa[7]
        elseif  sh == Subshell("3d_3/2")    wb = wa[8]
        elseif  sh == Subshell("3d_5/2")    wb = wa[9]
        elseif  sh == Subshell("4s_1/2")    wb = wa[10]
        elseif  sh == Subshell("4p_1/2")    wb = wa[11]
        elseif  sh == Subshell("4p_3/2")    wb = wa[12]
        else    error("No binding energy available for Z = $Z and subshell $sh ")
        end
    else   error("Unsupported keystring")
    end
    #
    if     wb == -1.   error("No binding energy available for Z = $Z and subshell $sh ")
    else   wb = Defaults.convertUnits("energy: from eV to atomic", wb)
    end
    #
    return( wb )
end

"""
  + `("binding energy", Z::Float64, conf::Configuration)`
    ... to provide an approximate binding energy of a given electron configuration. This estimate adds the binding
        energies of all subshell, taken frogm a semi-empirical tabulations by Williams et al.,
        https://userweb.jlab.org/~gwyn/ebindene.html. No relaxation effects are included if several hole states
        occur with regard to the neutral atom. An energy::Float64 in  Hartree is returned.
"""
function Basics.provide(sa::String, Z::Int64, conf::Configuration)
    if     sa == "binding energy"
        wa = JAC.PeriodicTable.bindingEnergies_Williams2000(Z);    wb = 0.
        for (sh,v) in  conf.shells
            if      sh == Shell("1s")    wb = wb + v * wa[1];    if wa[1]  == -1.   error("stop aa")   end
            elseif  sh == Shell("2s")    wb = wb + v * wa[2];    if wa[2]  == -1.   error("stop ab")   end
            elseif  sh == Shell("2p")    wb = wb + v * wa[3];    if wa[3]  == -1.   error("stop ac")   end
            elseif  sh == Shell("3s")    wb = wb + v * wa[5];    if wa[5]  == -1.   error("stop ae")   end
            elseif  sh == Shell("3p")    wb = wb + v * wa[6];    if wa[6]  == -1.   error("stop af")   end
            elseif  sh == Shell("3d")    wb = wb + v * wa[8];    if wa[8]  == -1.   error("stop ah")   end
            elseif  sh == Shell("4s")    wb = wb + v * wa[10];   if wa[10] == -1.   error("stop aj")   end
            elseif  sh == Shell("4p")    wb = wb + v * wa[11];   if wa[11] == -1.   error("stop ak")   end
            else    error("No binding energy available for Z = $Z and subshell $sh ")
            end
        end
    else   error("Unsupported keystring")
    end
    #
    wb = Defaults.convertUnits("energy: from eV to atomic", wb)
    #
    return( wb )
end