text:
export run_until!, next_event!

function  gillespie_step(pop::PopState, par::PopRates)
    total = pop.size * par.λ # total rate
    Δt = randexp() / total
    ind = rand(1:pop.size) # choose an individual at random
    if  (2*rand() - 1)*par.λ + pop.fitness[ind] - pop.offset > 0
            process = ind
    else
            process = -ind
    end
    return (process, Δt)
    # sign of process tells us whether it's a birth or death
end

# Mutation

function flipat!(bits::BitArray, inds)
    for ii in inds
        bits[ii] = !bits[ii]
    end
end

"""
Mutates the genome at a rate given by μ
Assumes the action of a markov chain so that
mutate!(X, 2μ) = mutate!(mutate!(X, μ), μ).
mutate!(Inf) should randomize the genome
"""

function mutate!(genome::BitArray, ν::Real)
    N = length(genome)
    p = -expm1(-2ν)/2
    flips = randsubseq(1:N, p)
    flipat!(genome, flips)
end

"""
Recombines two genomes at a rate given by ν
Assumes the action of a markov chain so that
mutate!(X, 2ν) = mutate!(mutate!(X, ν), ν).
mutate!(Inf) should randomize the genome

To undertand the code, here we have randomly chosen cuts at [2,5].
From positions 3-5 the codes are swapped as in this diagram.
*--. ,---. ,----->
    /     /
*--' `---' `----->

"""

function recombine!(v1::T, v2::T, ρ::Real) where T <: AbstractArray
    N = length(v1)
    cut = sort(randsubseq(0:(N-1),ρ))
    if isodd(length(cut))
        push!(cut,N)
    end
    for ii in 1:2:length(cut)
         for jj in cut[ii]+1:cut[ii+1]
             (v2[jj], v1[jj]) = (v1[jj], v2[jj])
         end
    end
end

"""
Add an individual to the population.
Adjust the birth and death rates accordingly
"""
function add_individual!(pop::PopState, child::BitArray{1}, par::PopRates)
    # add to the fitness
    push!(pop.individuals, child)
    push!(pop.fitness, par.f(child))
    pop.offset += last(pop.fitness)/par.κ      # add to competition
    pop.size += 1
end

"""
Delete an individual to the population.
Adjust the birth and death rates accordingly
"""
function delete_individual!(pop::PopState, index::Int, par::PopRates)
    # add to the death rate
    pop.offset +=  -pop.fitness[index]/par.κ
    deleteat!(pop.fitness, index)
    deleteat!(pop.individuals, index)
    pop.size += -1
end

function birth_event!(pop::PopState, parent_index::Int, par::PopRates)
    child = copy(pop.individuals[parent_index])
    mutate!(child, par.ν)
    # add to the population
    add_individual!(pop, child, par)
end

death_event!(pop::PopState, index::Int, par::PopRates) = delete_individual!(pop, index, par)

"""
Choose an individual from the population (daddy) to recombine with the individual at
mommy_index.
"""

function recombination_event!(pop::PopState, mommy_index::Int, par::PopRates)
    daddy_index = rand(1:pop.size)
    # update their genomes
    @views recombine!(pop.individuals[mommy_index], pop.individuals[daddy_index], par.ρ)
    # the fitness differences
    Δmfit = par.f(pop.individuals[mommy_index]) - pop.fitness[mommy_index]
    Δdfit = par.f(pop.individuals[daddy_index]) - pop.fitness[daddy_index]
    pop.fitness[mommy_index] += Δmfit
    pop.fitness[daddy_index] += Δdfit
    # update offset
    pop.offset += (Δmfit+Δdfit)/par.κ
end

function next_event!(pop::PopState, par::PopRates)
    # step forward in time
    (index, Δt) = gillespie_step(pop, par)
    if index > 0
        birth_event!(pop, index, par)
        if rand(Bernoulli(par.χ))
            recombination_event!(pop, index, par)
        end
    else
        death_event!(pop, abs(index), par)
    end
    pop.time += Δt
end

"""
Measure diversity of neutral sites
    data is a matrix of counts ( 2 x n)
"""

function run_until!(pop::PopState, par::PopRates, t)
    renew_fitness(pop,par)
        # really only neccesary if the parameters have changed
        # slows things down if timesteps are recorded at order generation time
    while pop.time < t
        next_event!(pop, par)
    end
end

# function run_sim(pop::PopState, par::PopRates, timepoints;
#     stats = [freq, D, mean_fit, var_fit],
#     names = [:freq, :D, :fit, :fitvar])
#     rows = []
#     time = 0
#     for t in timepoints
#         while time < t
#             time += next_event!(pop, par)
#         end
#         push!(rows,[t;[s(pop,par) for s in stats]])
#     end
#     columns = Vector{Vector}(undef, 0)
#     for (i,r) in enumerate(first(rows))
#         column = Vector{typeof(r)}(undef, length(rows))
#         column .= getindex.(rows, i)
#         push!(columns, column)
#        end
#     df = DataFrame(columns, [:t;names], copycols=false)
# end

##

max_stars_repo_name: StatPhysBio/Tomoko

max_issues_repo_name: StatPhysBio/Tomoko

max_forks_repo_name: StatPhysBio/Tomoko


text:
module Dictionaries

using Random

using Indexing

using Base: @propagate_inbounds, Callable

export getindices, setindices!, mapview

export AbstractDictionary, AbstractIndices, IndexError, Indices, HashIndices, HashDictionary, Dictionary, MappedDictionary, DictionaryView, FilteredDictionary, FilteredIndices, BroadcastedDictionary

export issettable, isinsertable, set!, unset!, dictionary
export istokenizable, tokentype, tokens, tokenized, gettoken, gettokenvalue, istokenassigned, settokenvalue!, gettoken!, deletetoken!, sharetokens

export filterview # TODO move to SplitApplyCombine.jl (and re-order project dependencies?)

include("AbstractDictionary.jl")
include("AbstractIndices.jl")

include("tokens.jl")
include("indexing.jl")
include("insertion.jl")
include("filter.jl")
include("foreach.jl")
include("map.jl")
include("broadcast.jl")
include("find.jl")
include("show.jl")

include("PairDictionary.jl")
include("Indices.jl")
include("Dictionary.jl")
include("HashIndices.jl")
include("HashDictionary.jl")
include("MappedDictionary.jl")

end # module

# # TODO
#
# * Improved printing - don't calculate length (beyond some cutoff) if it is `SizeUnknown` and limit=true.
# * `hash` and `isless`
# * A manually-ordered dictionary would be quite useful, like [OrderedCollections.jl](https://github.com/JuliaCollections/OrderedCollections.jl).
# * A surface interface for updates like https://github.com/JuliaLang/julia/pull/31367
# * Soon we will have the concept of "ordered" indices/sets (sort-based dictionaries and
#   B-trees). We can probably formalize an interface around a trait here. Certain operations
#   like slicing out an interval or performing a sort-merge co-iteration for `merge` become
#   feasible.

ext: jl


text:
import Distributions

function effective_sample_size(log_normalized_weights::Vector{Float64})
    log_ess = -logsumexp(2. * log_normalized_weights)
    return exp(log_ess)
end

function normalize_weights(log_weights::Vector{Float64})
    log_total_weight = logsumexp(log_weights)
    log_normalized_weights = log_weights .- log_total_weight
    return (log_total_weight, log_normalized_weights)
end

#######################################
# building blocks for particle filter #
#######################################

mutable struct ParticleFilterState{U}
    traces::Vector{U}
    new_traces::Vector{U}
    log_weights::Vector{Float64}
    log_ml_est::Float64
    parents::Vector{Int}
end

"""
    traces = get_traces(state::ParticleFilterState)

Return the vector of traces in the current state, one for each particle.
"""
function get_traces(state::ParticleFilterState)
    state.traces

end

"""
    log_weights = get_log_weights(state::ParticleFilterState)

Return the vector of log weights for the current state, one for each particle.

The weights are not normalized, and are in log-space.
"""
function get_log_weights(state::ParticleFilterState)
    state.log_weights
end

"""
    estimate = log_ml_estimate(state::ParticleFilterState)

Return the particle filter's current estimate of the log marginal likelihood.
"""
function log_ml_estimate(state::ParticleFilterState)
    num_particles = length(state.traces)
    return state.log_ml_est + logsumexp(state.log_weights) - log(num_particles)
end

"""
    traces::Vector = sample_unweighted_traces(state::ParticleFilterState, num_samples::Int)

Sample a vector of `num_samples` traces from the weighted collection of traces in the given particle filter state.
"""
function sample_unweighted_traces(state::ParticleFilterState{U}, num_samples::Int) where {U}
    (_, log_normalized_weights) = normalize_weights(state.log_weights)
    weights = exp.(log_normalized_weights)
    traces = Vector{U}(undef, num_samples)
    for i=1:num_samples
        traces[i] = state.traces[categorical(weights)]
    end
    return traces
end

"""
    state = initialize_particle_filter(model::GenerativeFunction, model_args::Tuple,
        observations::ChoiceMap, proposal::GenerativeFunction, proposal_args::Tuple,
        num_particles::Int)

Initialize the state of a particle filter using a custom proposal for the initial latent state.
"""
function initialize_particle_filter(model::GenerativeFunction{T,U}, model_args::Tuple,
        observations::ChoiceMap, proposal::GenerativeFunction, proposal_args::Tuple,
        num_particles::Int) where {T,U}
    traces = Vector{Any}(undef, num_particles)
    log_weights = Vector{Float64}(undef, num_particles)
    for i=1:num_particles
        (prop_choices, prop_weight, _) = propose(proposal, proposal_args)
        (traces[i], model_weight) = generate(model, model_args, merge(observations, prop_choices))
        log_weights[i] = model_weight - prop_weight
    end
    ParticleFilterState{U}(traces, Vector{U}(undef, num_particles),
        log_weights, 0., collect(1:num_particles))
end

"""
    state = initialize_particle_filter(model::GenerativeFunction, model_args::Tuple,
        observations::ChoiceMap, num_particles::Int)

Initialize the state of a particle filter, using the default proposal for the initial latent state.
"""
function initialize_particle_filter(model::GenerativeFunction{T,U}, model_args::Tuple,
        observations::ChoiceMap, num_particles::Int) where {T,U}
    traces = Vector{Any}(undef, num_particles)
    log_weights = Vector{Float64}(undef, num_particles)
    for i=1:num_particles
        (traces[i], log_weights[i]) = generate(model, model_args, observations)
    end
    ParticleFilterState{U}(traces, Vector{U}(undef, num_particles),
        log_weights, 0., collect(1:num_particles))
end

"""
    (log_incremental_weights,) = particle_filter_step!(
        state::ParticleFilterState, new_args::Tuple, argdiffs,
        observations::ChoiceMap, proposal::GenerativeFunction, proposal_args::Tuple)

Perform a particle filter update, where the model arguments are adjusted, new
observations are added, and some combination of a custom proposal and the
model's internal proposal is used for proposing new latent state.  That is, for
each particle,

* The proposal function `proposal` is evaluated with arguments `Tuple(t_old,
  proposal_args...)` (where `t_old` is the old model trace), and produces its
  own trace (call it `proposal_trace`); and
* The old model trace is replaced by a new model trace (call it `t_new`).

The choicemap of `t_new` satisfies the following conditions:

1. `get_choices(t_old)` is a subset of `get_choices(t_new)`;
2. `observations` is a subset of `get_choices(t_new)`;
3. `get_choices(proposal_trace)` is a subset of `get_choices(t_new)`.

Here, when we say one choicemap `a` is a "subset" of another choicemap `b`, we
mean that all keys that occur in `a` also occur in `b`, and the values at those
addresses are equal.

It is an error if no trace `t_new` satisfying the above conditions exists in
the support of the model (with the new arguments). If such a trace exists, then
the random choices not determined by the above requirements are sampled using
the internal proposal.
"""
function particle_filter_step!(state::ParticleFilterState{U}, new_args::Tuple, argdiffs::Tuple,
        observations::ChoiceMap, proposal::GenerativeFunction, proposal_args::Tuple) where {U}
    trace_translator = SimpleExtendingTraceTranslator(new_args, argdiffs, observations, proposal, proposal_args)
    num_particles = length(state.traces)
    log_incremental_weights = Vector{Float64}(undef, num_particles)
    for i=1:num_particles
        (state.new_traces[i], log_weight) = trace_translator(state.traces[i])
        log_incremental_weights[i] = log_weight
        state.log_weights[i] += log_weight
    end

    # swap references
    tmp = state.traces
    state.traces = state.new_traces
    state.new_traces = tmp

    return (log_incremental_weights,)
end

"""
    (log_incremental_weights,) = particle_filter_step!(
        state::ParticleFilterState, new_args::Tuple, argdiffs,
        observations::ChoiceMap)

Perform a particle filter update, where the model arguments are adjusted, new observations are added, and the default proposal is used for new latent state.
"""
function particle_filter_step!(state::ParticleFilterState{U}, new_args::Tuple, argdiffs::Tuple,
        observations::ChoiceMap) where {U}
    num_particles = length(state.traces)
    log_incremental_weights = Vector{Float64}(undef, num_particles)
    for i=1:num_particles
        (state.new_traces[i], increment, _, discard) = update(
            state.traces[i], new_args, argdiffs, observations)
        if !isempty(discard)
            error("Choices were updated or deleted inside particle filter step: $discard")
        end
        log_incremental_weights[i] = increment
        state.log_weights[i] += increment
    end

    # swap references
    tmp = state.traces
    state.traces = state.new_traces
    state.new_traces = tmp

    return (log_incremental_weights,)
end

"""
    did_resample::Bool = maybe_resample!(state::ParticleFilterState;
        ess_threshold::Float64=length(state.traces)/2, verbose=false)

Do a resampling step if the effective sample size is below the given threshold.
Return `true` if a resample thus occurred, `false` otherwise.
"""
function maybe_resample!(state::ParticleFilterState{U};
                        ess_threshold::Real=length(state.traces)/2, verbose=false) where {U}
    num_particles = length(state.traces)
    (log_total_weight, log_normalized_weights) = normalize_weights(state.log_weights)
    ess = effective_sample_size(log_normalized_weights)
    do_resample = ess < ess_threshold
    if verbose
        println("effective sample size: $ess, doing resample: $do_resample")
    end
    if do_resample
        weights = exp.(log_normalized_weights)
        Distributions.rand!(Distributions.Categorical(weights / sum(weights)), state.parents)
        state.log_ml_est += log_total_weight - log(num_particles)
        for i=1:num_particles
            state.new_traces[i] = state.traces[state.parents[i]]
            state.log_weights[i] = 0.
        end

        # swap references
        tmp = state.traces
        state.traces = state.new_traces
        state.new_traces = tmp
    end
    return do_resample
end

export initialize_particle_filter, particle_filter_step!, maybe_resample!
export get_traces, get_log_weights, log_ml_estimate, sample_unweighted_traces


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

mutable struct UpgradeOperationHistoricalStatusInfo <: SwaggerModel
    properties::Any # spec type: Union{ Nothing, UpgradeOperationHistoricalStatusInfoProperties } # spec name: properties
    type::Any # spec type: Union{ Nothing, String } # spec name: type
    location::Any # spec type: Union{ Nothing, String } # spec name: location

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

const _property_map_UpgradeOperationHistoricalStatusInfo = Dict{Symbol,Symbol}(Symbol("properties")=>Symbol("properties"), Symbol("type")=>Symbol("type"), Symbol("location")=>Symbol("location"))
const _property_types_UpgradeOperationHistoricalStatusInfo = Dict{Symbol,String}(Symbol("properties")=>"UpgradeOperationHistoricalStatusInfoProperties", Symbol("type")=>"String", Symbol("location")=>"String")
Base.propertynames(::Type{ UpgradeOperationHistoricalStatusInfo }) = collect(keys(_property_map_UpgradeOperationHistoricalStatusInfo))
Swagger.property_type(::Type{ UpgradeOperationHistoricalStatusInfo }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_UpgradeOperationHistoricalStatusInfo[name]))}
Swagger.field_name(::Type{ UpgradeOperationHistoricalStatusInfo }, property_name::Symbol) =  _property_map_UpgradeOperationHistoricalStatusInfo[property_name]

function check_required(o::UpgradeOperationHistoricalStatusInfo)
    true
end

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


text:
using BSON
using Dates
using Distributed
using FilePathsBase: SystemPath
using InteractiveUtils
using Memento
using Pkg
using Random
using Serialization
using Suppressor
using Test

using JLSO
using JLSO: JLSOFile, LOGGER, upgrade_jlso

# To test different types from common external packages
using DataFrames
using Distributions
using TimeZones

include("test_data_setup.jl")
@testset "JLSO" begin
    include("backwards_compat.jl")
    include("JLSOFile.jl")
    include("file_io.jl")
end


text:
@kwdef mutable struct DAddHandler <: AbstractHandler
    id = 0
end

@register(DAddHandler, +, [DBit, DBit] => [Number])

function (handler::DAddHandler)(netlist::Netlist,
                                inputs::Vector{Variable},
                                outputs::Vector{Variable})
    # check size
    if getsize(netlist, getname(inputs[1])) != (1, 1) && getsize(netlist, getname(inputs[1])) != (1, 1)
        error("Deterministic add does not support matrix hardware generation yet.")
    end

    outstring = """
        $stdcomment
        // BEGIN d_add$(handler.id)
        determ_add #(
                .BIT_WIDTH(BIT_WIDTH),
                .INT_WIDTH(INT_WIDTH)
            ) d_add$(handler.id) (
                .a($(inputs[1].name)),
                .b($(inputs[2].name)),
                .y($(outputs[1].name))
            );
        // END d_add$(handler.id)
        \n"""

    handler.id += 1

    return outstring
end

@kwdef mutable struct DAddFxpHandler <: AbstractHandler
    id = 0
end

@register(DAddFxpHandler, +, begin
    [DBit, Number] => [Number]
    [Number, DBit] => [Number]
end)

function (handler::DAddFxpHandler)(netlist::Netlist,
                                   inputs::Vector{Variable},
                                   outputs::Vector{Variable})
    # check size
    if getsize(netlist, getname(inputs[1])) != (1, 1) && getsize(netlist, getname(inputs[1])) != (1, 1)
        error("Deterministic add does not support matrix hardware generation yet.")
    end

    outstring = """
        $stdcomment
        // BEGIN d_add_fxp$(handler.id)
        determ_add_fxp #(
                .BIT_WIDTH(BIT_WIDTH),
                .INT_WIDTH(INT_WIDTH)
            ) d_add_fxp$(handler.id) (
                .a($(inputs[1].name)),
                .b($(inputs[2].name)),
                .y($(outputs[1].name))
            );
        // END d_add_fxp$(handler.id)
        \n"""

    handler.id += 1

    return outstring
end


text:
"""
    EmulatedABCRejection(
        reference_data,
        simulator_function,
        priors,
        threshold,
        n_particles,
        n_design_points;
        summary_statistic               = "keep_all",
        distance_function               = Distances.euclidean,
        batch_size                      = 10*n_particles,
        max_iter                        = 1000,
        emulator_training               = DefaultEmulatorTraining(),
        emulated_particle_selection     = MeanEmulatedParticleSelection(),
        write_progress                  = true,
        progress_every                  = 1000,
        )

Run emulation-based rejection ABC algorithm. Model simulation results are used to train the
emulator, which is then used to get the approximated distance for each particle. The rest
of the workflow is similar to [`SimulatedABCRejection`](@ref).

See [ABC Overview](@ref abc-overview) for more details.

# Mandatory arguments
- `reference_data::AbstractArray{Float,2}`: Observed data to which the simulated model output will be compared. Array dimensions sould match that of the simulator function result.
- `simulator_function::Function`: A function that takes a parameter vector as an argument and outputs model results.
- `priors::AbstractArray{ContinuousUnivariateDistribution,1}`: Continuous univariate distributions, from which candidate parameters will be sampled. Array size should match the number of parameters.
- `threshold::Float`: The ``\\varepsilon`` threshold to be used in ABC algorithm. Only those particles that produce emulated results that are within this threshold from the reference data are included into the posterior.
- `n_particles::Int`: The number of parameter vectors (particles) that will be included in the posterior.
- `n_design_points::Int`: The number of design particles that will be simulated to traing the emulator.

# Optional keyword arguments
- `summary_statistic::Union{String,AbstractArray{String,1},Function}`: Summary statistics that will be applied to the data before computing the distances. Defaults to `keep_all`. See [detailed documentation of summary statistics](@ref summary_stats).
- `distance_function::Function`: A function that will be used to compute the distance between the summary statistic of the simulated data and that of reference data. Defaults to `Distances.euclidean`.
- `batch_size::Int`: The number of particles that will be emulated on each iteration. Defaults to `1000 * n_particles`.
- `max_iter::Int`: The maximum number of emulations that will be run. Defaults to 1000.
- `emulator_training<:AbstractEmulatorTraining`: This determines how the emulator will be trained. See [`AbstractEmulatorTraining`](@ref) for more details.
- `emulated_particle_selection<:AbstractEmulatedParticleSelection`: This determines how the particles that will be added to the posterior are selected after each emulation run. See [`AbstractEmulatedParticleSelection`](@ref) for details. Defaults to [`MeanEmulatedParticleSelection`](@ref).
- `write_progress::Bool`: Whether algorithm progress should be printed on standard output. Defaults to `true`.

# Returns
An [`ABCRejectionOutput`](@ref) object.
"""
function EmulatedABCRejection(reference_data::AbstractArray{AF,2},
    simulator_function::Function,
    priors::AbstractArray{D,1},
    threshold::AF,
    n_particles::Int,
    n_design_points::Int;
    summary_statistic::Union{String,AbstractArray{String,1},Function}="keep_all",
    distance_function::Function=Distances.euclidean,
    emulation_type::AbstractEmulatorTraining = DefaultEmulatorTraining(),
    batch_size::Int=10*n_particles,
    max_iter::Int=1000,
    emulator_training::ET = DefaultEmulatorTraining(),
    emulated_particle_selection::EPS = MeanEmulatedParticleSelection(),
    kwargs...) where {
    AF<:AbstractFloat,
    D<:ContinuousUnivariateDistribution,
    EPS<:AbstractEmulatedParticleSelection,
    ET<:AbstractEmulatorTraining
    }

    summary_statistic = build_summary_statistic(summary_statistic)
    reference_summary_statistic = summary_statistic(reference_data)

    emulator_training_input = EmulatorTrainingInput(
        DistanceSimulationInput(reference_summary_statistic, simulator_function, summary_statistic, distance_function),
        n_design_points,
        emulator_training)

    input = EmulatedABCRejectionInput(length(priors), n_particles, threshold,
        priors, batch_size, max_iter, emulator_training_input, emulated_particle_selection)

    return ABCrejection(input; kwargs...)
end

"""
    EmulatedABCSMC(
        reference_data,
        simulator_function,
        priors,
        threshold_schedule,
        n_particles,
        n_design_points;
        summary_statistic               = "keep_all",
        distance_function               = Distances.euclidean,
        batch_size                      = 10*n_particles,
        max_iter                        = 1000,
        emulator_training               = DefaultEmulatorTraining(),
        emulator_retraining             = NoopRetraining(),
        emulated_particle_selection     = MeanEmulatedParticleSelection(),
        write_progress                  = true,
        progress_every                  = 1000,
        )

Run emulation-based ABC-SMC algorithm. This is similar to [`EmulatedABCRejection`](@ref),
the main difference being that an array of thresholds is provided instead of a single threshold.
It is assumed that thresholds are sorted in decreasing order.

An emulation based ABC iteration is executed for each threshold. For the first threshold, the provided prior is used.
For each subsequent threshold, the posterior from the previous iteration is used as a prior.

See [ABC Overview](@ref abc-overview) for more details.

# Mandatory arguments
- `reference_data::AbstractArray{Float,2}`: Observed data to which the simulated model output will be compared. Array dimensions sould match that of the simulator function result.
- `simulator_function::Function`: A function that takes a parameter vector as an argument and outputs model results.
- `priors::AbstractArray{ContinuousUnivariateDistribution,1}`: Continuous univariate distributions, from which candidate parameters will be sampled during the first iteration. Array size should match the number of parameters.
- `threshold_schedule::AbstractArray{Float,1}`: The threshold schedule to be used in ABC algorithm. An ABC iteration is executed for each threshold. It is assumed that thresholds are sorted in decreasing order.
- `n_particles::Int`: The number of parameter vectors (particles) that will be included in the posterior.
- `n_design_points::Int`: The number of design particles that will be simulated to traing the emulator.

# Optional keyword arguments
- `summary_statistic::Union{String,AbstractArray{String,1},Function}`: Summary statistics that will be applied to the data before computing the distances. Defaults to `keep_all`. See [detailed documentation of summary statistics](@ref summary_stats).
- `distance_function::Function`: A function that will be used to compute the distance between the summary statistic of the simulated data and that of reference data. Defaults to `Distances.euclidean`.
- `batch_size::Int`: The number of particles that will be emulated on each iteration. Defaults to `1000 * n_particles`.
- `max_iter::Int`: The maximum number of emulations that will be run. Defaults to 1000.
- `emulator_training<:AbstractEmulatorTraining`: This determines how the emulator will be trained for each iteration. See [`AbstractEmulatorTraining`](@ref) for more details.
- `emulator_retraining<:AbstractEmulatorRetraining`: This is used to specify parameters of additional emulator retraining that can be done for each iteration. By default this retraining is switched off ([`NoopRetraining`](@ref)). See [`AbstractEmulatorRetraining`] for more details.
- `emulated_particle_selection<:AbstractEmulatedParticleSelection`: This determines how the particles that will be added to the posterior are selected after each emulation run. See [`AbstractEmulatedParticleSelection`](@ref) for details. Defaults to [`MeanEmulatedParticleSelection`](@ref).
- `write_progress::Bool`: Whether algorithm progress should be printed on standard output. Defaults to `true`.

# Returns
An [`ABCSMCOutput`](@ref) object.
"""
function EmulatedABCSMC(reference_data::AbstractArray{AF,2},
    simulator_function::Function,
    priors::AbstractArray{D,1},
    threshold_schedule::AbstractArray{AF,1},
    n_particles::Int,
    n_design_points::Int;
    summary_statistic::Union{String,AbstractArray{String,1},Function}="keep_all",
    emulator_training::ET = DefaultEmulatorTraining(),
    distance_metric::Function=Distances.euclidean,
    batch_size::Int=10*n_particles,
    max_iter::Int=20,
    emulator_retraining::ER = NoopRetraining(),
    emulated_particle_selection::EPS = MeanEmulatedParticleSelection(),
    kwargs...) where {
    AF<:AbstractFloat,
    ET<:AbstractEmulatorTraining,
    ER<:AbstractEmulatorRetraining,
    EPS<:AbstractEmulatedParticleSelection,
    D<:ContinuousUnivariateDistribution
    }

    summary_statistic = build_summary_statistic(summary_statistic)
    reference_summary_statistic = summary_statistic(reference_data)

    emulator_training_input = EmulatorTrainingInput(
        DistanceSimulationInput(reference_summary_statistic, simulator_function, summary_statistic, distance_metric),
        n_design_points,
        emulator_training)

    input = EmulatedABCSMCInput(length(priors), n_particles, threshold_schedule,
        priors, batch_size, max_iter, emulator_training_input, emulator_retraining, emulated_particle_selection)

    return ABCSMC(input; kwargs...)
end

"""
    EmulatedModelSelection

Perform model selection using emulation-based ABC.

# Arguments
- `n_design_points::Int64`: The number of parameter vectors used to train the Gaussian process emulator.
- `reference_data::AbstractArray{Float64,2}`: The observed data to which the simulated model output will be compared. Size: (n_model_trajectories, n_time_points)
- `n_particles::Int64`: The number of parameter vectors (particles) that will be included in the final posterior.
- `threshold_schedule::AbstractArray{Float64}`: A set of maximum distances from the summarised model output to summarised observed data for a parameter vector to be included in the posterior. Each distance will be used in a single run of the ABC-SMC algorithm.
- `parameter_priors::AbstractArray{AbstractArray{ContinuousUnivariateDistribution},1}`: Priors for the parameters of each model. The length of the outer array is the number of models.
- `summary_statistic::Union{String,AbstractArray{String,1},Function}`: Either: 1. A `String` or 1D Array of strings that Or 2. A function that outputs a 1D Array of Floats that summarises model output. Defaults to `keep_all`. See [detailed documentation of summary statistics](@ref summary_stats).
- `simulator_functions::AbstractArray{Function,1}`: An array of functions that take a parameter vector as an argument and outputs model results (one per model).
- 'model_prior::DiscreteUnivariateDistribution': The prior from which models are sampled. Default is a discrete, uniform distribution.
- `distance_function::Function`: Any function that computes the distance between 2 1D Arrays. Optional argument (default is to use the Euclidean distance).
- `max_iter::Int`: The maximum number of simulations that will be run. The default is 1000*`n_particles`. Each iteration samples a single model and performs ABC using a single particle.
- `max_batch_size::Int`: The maximum batch size for the emulator when making predictions.

# Returns
A [`ModelSelectionOutput`](@ref) object that contains which models are supported by the observed data.
"""
function EmulatedModelSelection(
    reference_data::AbstractArray{AF,2},
    simulator_functions::AbstractArray{Function,1},
    parameter_priors::AbstractArray{AD,1},
    threshold_schedule::AbstractArray{AF,1},
    n_particles::Int,
    n_design_points::Int,
    model_prior::DiscreteUnivariateDistribution=Distributions.DiscreteUniform(1,length(parameter_priors));
    summary_statistic::Union{String,AbstractArray{String,1},Function}="keep_all",
    distance_function::Function=Distances.euclidean,
    max_iter::Int=1000,
    max_batch_size::Int=1000,
    emulator_training::ET = DefaultEmulatorTraining(),
    emulator_retraining::ER = NoopRetraining(),
    emulated_particle_selection::EPS = MeanEmulatedParticleSelection(),
    ) where {
    D<:ContinuousUnivariateDistribution,
    AD<:AbstractArray{D,1},
    AF<:AbstractFloat,
    ET<:AbstractEmulatorTraining,
    ER<:AbstractEmulatorRetraining,
    EPS<:AbstractEmulatedParticleSelection,
    }

    summary_statistic = build_summary_statistic(summary_statistic)
    reference_summary_statistic = summary_statistic(reference_data)
    emulator_training_input = [EmulatorTrainingInput(
        DistanceSimulationInput(reference_summary_statistic, simulator_functions[m], summary_statistic, distance_function),
        n_design_points, emulator_training) for m in 1:length(parameter_priors)]

    input = EmulatedModelSelectionInput(length(parameter_priors),
        n_particles,
        threshold_schedule,
        model_prior,
        parameter_priors,
        emulator_training_input,
        emulator_retraining,
        emulated_particle_selection,
        max_batch_size,
        max_iter)

    return model_selection(input)
end


text:
module test_stdlib_interactiveutils

using InteractiveUtils, Test

@test @which(+) === Base

meth = @which(1+2)
@test meth.file === Symbol("int.jl")

end # module test_stdlib_interactiveutils


text:
using LibGit2

mutable struct Server
   storedir::String
   context::Pkg.Types.Context
   depot::Dict
end

struct PackageID
    name::String
    uuid::String
end

abstract type SymStore end

mutable struct ModuleStore <: SymStore
    name::String
    vals::Dict{String,Any}
    exported::Set{String}
    doc::String
    ver::String
    sha
end
ModuleStore(name::String) = ModuleStore(name, Dict{String,Any}(), Set{String}(), "", "", nothing)

struct MethodStore <: SymStore
    file::String
    line::Int
    args::Vector{Tuple{String,String}}
end

struct FunctionStore <: SymStore
    methods::Vector{MethodStore}
    doc::String
end

struct abstractStore <: SymStore
    params::Vector{String}
    doc::String
end

struct primitiveStore <: SymStore
    params::Vector{String}
    doc::String
end

struct structStore <: SymStore
    params::Vector{String}
    fields::Vector{String}
    ts::Vector{String}
    methods::Vector{MethodStore}
    doc::String
end

struct genericStore <: SymStore
    t::String
    params::Vector{String}
    doc::String
end

function _getdoc(x)
    # Packages can add methods to Docs.doc, and those can have a bug,
    # and we don't want that to kill the symbol server process
    try
        return string(Docs.doc(x))
    catch err
        @warn "Couldn't retrieve docs."
        return ""
    end
end

function read_methods(x)
    ms = methods(x)
    map(ms) do m
        path = isabspath(String(m.file)) ? String(m.file) : Base.find_source_file(String(m.file))
        if path === nothing
            path = ""
        end
        args = Base.arg_decl_parts(m)[2][2:end]
        if isdefined(ms.mt, :kwsorter)
            kws = Base.kwarg_decl(m, typeof(ms.mt.kwsorter))
            for kw in kws
                push!(args, (string(kw), ".KW"))
            end
        end
        for i = 1:length(args)
            if isempty(args[i][2])
                args[i] = (args[i][1], "Any")
            end
        end
        MethodStore(path,
                    m.line,
                    args)
    end
end

function collect_params(t, params = [])
    if t isa UnionAll
        push!(params, t.var)
        return collect_params(t.body, params)
    else
        return t, params
    end
end

function import_package_names(pkg::PackageID, depot, c, m = nothing)
    if pkg.uuid in keys(depot)
        return depot[pkg.uuid]
    else
        depot[pkg.uuid] = ModuleStore(pkg.name)
        depot[pkg.uuid].ver = pkg_ver(pkg, c)
    end
    path = pkg_path(pkg, c)
    if path isa String && isdir(path) && isgitrepo(path)
        depot[pkg.uuid].sha = getgithash(path)
    end
    if m isa Module
    elseif Symbol(pkg.name) in names(Main, all = true)
        m = getfield(Main, Symbol(pkg.name))
    else
        m = try
            Main.eval(:(import $(Symbol(pkg.name))))
            m = getfield(Main, Symbol(pkg.name))
        catch
            nothing
        end
    end
    if m isa Module
        get_module_names(m, pkg, depot, depot[pkg.uuid], c)
    end
    for dep in pkg_deps(pkg, c)
        depid = PackageID(first(dep), string(last(dep)))
        if !haskey(depot, depid.uuid)
            import_package_names(depid, depot, c)
        end
    end
    return depot[pkg.uuid]
end

function get_module_names(m::Module, pkg::PackageID, depot, out::ModuleStore, c::Pkg.Types.Context)
    out.doc = string(Docs.doc(m))
    out.exported = Set{String}(string.(names(m)))
    allnames = names(m, all = true, imported = true)
    for n in allnames
        !isdefined(m, n) && continue
        startswith(string(n), "#") && continue
        if Base.isdeprecated(m, n)
        else
            x = getfield(m, n)
            t, p = collect_params(x)
            if x isa Function
                out.vals[String(n)] = FunctionStore(read_methods(x), _getdoc(x))
            elseif t isa DataType
                if t.abstract
                    out.vals[String(n)] = abstractStore(string.(p), _getdoc(x))
                elseif t.isbitstype
                        out.vals[String(n)] = primitiveStore(string.(p), _getdoc(x))
                elseif !(isempty(t.types) || Base.isvatuple(t)) || t.mutable
                        out.vals[String(n)] = structStore(string.(p),
                                                     collect(string.(fieldnames(t))),
                                                     string.(collect(t.types)),
                                                     read_methods(x),
                                                     _getdoc(x))
                else
                    out.vals[String(n)] = genericStore("DataType", string.(p), _getdoc(x))
                end
            elseif x isa Module && x != m # include reference to current module
                if parentmodule(x) == m # load non-imported submodules
                    out.vals[String(n)] = ModuleStore(String(n))
                    get_module_names(x, pkg, depot, out.vals[String(n)], c)
                end
            else
                out.vals[String(n)] = genericStore(string(typeof(x)), [], _getdoc(x))
            end
        end
    end
    for dep in pkg_deps(pkg, c)
        depid = PackageID(first(dep), string(last(dep)))
        dep_module = can_access(m, Symbol(depid.name))
        if dep_module isa Module
            out.vals[depid.name] = depid.name
            if !haskey(depot, depid.uuid)
                import_package_names(depid, depot, c, dep_module)
            end
        end
    end
end

function load_core()
    c = Pkg.Types.Context()
    depot = Dict{String,Any}()
    for m in (Base,Core)
        depot[string(m)] = ModuleStore(string(m))
        get_module_names(m, PackageID(string(m), ""), depot, depot[string(m)], c)
    end

    # Add special cases
    push!(depot["Base"].exported, "include")
    depot["Base"].vals["@."] = depot["Base"].vals["@__dot__"]
    push!(depot["Base"].exported, "@.")
    delete!(depot["Core"].exported, "Main")
    # Add built-ins
    add_builtin_methods(depot["Core"])
    return depot
end

function add_builtin_methods(core_dep::ModuleStore)
    builtins = (split("=== typeof sizeof <: isa typeassert throw tuple getfield setfield! fieldtype nfields isdefined arrayref arrayset arraysize applicable invoke apply_type _apply _expr svec"))
    for f in builtins
        if haskey(core_dep.vals, f)
            push!(core_dep.vals[f].methods, MethodStore("built-in", 0, [("args...", "Any")]))
        else
            core_dep.vals[f] = FunctionStore(MethodStore[MethodStore("built-in", 0, [("args...", "Any")])], _getdoc(getfield(Core, Symbol(f))))
        end
    end
end

function save_store_to_disc(store, file)
    io = open(file, "w")
    serialize(io, store)
    close(io)
end

function pkg_deps(pkg::PackageID, c::Pkg.Types.Context)
    if VERSION < v"1.1.0-DEV.857"
        if haskey(c.env.manifest, pkg.name)
            return get(c.env.manifest[pkg.name][1], "deps", Dict{Any,Any}())
        else
            return Dict{Any,Any}()
        end
    else
        if !isempty(pkg.uuid) && haskey(c.env.manifest, Base.UUID(pkg.uuid))
            return c.env.manifest[Base.UUID(pkg.uuid)].deps
        else
            return Dict{String,Base.UUID}()
        end
    end
end

function context_deps(c::Pkg.Types.Context)
    if VERSION < v"1.1.0-DEV.857"
        c.env.project["deps"]
    else
        c.env.project.deps
    end
end

function pkg_ver(pkg::PackageID, c::Pkg.Types.Context)
    if VERSION < v"1.1.0-DEV.857"
        if haskey(c.env.manifest, pkg.name)
            return get(c.env.manifest[pkg.name][1], "version", "")
        else
            return ""
        end
    else
        if !isempty(pkg.uuid) && haskey(c.env.manifest, Base.UUID(pkg.uuid))
            return get(c.env.manifest[Base.UUID(pkg.uuid)].other, "version", "")
        else
            return ""
        end
    end
end

function pkg_path(pkg::PackageID, c::Pkg.Types.Context)
    if VERSION < v"1.1.0-DEV.857"
        if haskey(c.env.manifest, pkg.name)
            return get(c.env.manifest[pkg.name][1], "path", "")
        else
            return ""
        end
    else
        if !isempty(pkg.uuid) && haskey(c.env.manifest, Base.UUID(pkg.uuid))
            return get(c.env.manifest[Base.UUID(pkg.uuid)].other, "path", "")
        else
            return ""
        end
    end
end

function get_manifest(c::Pkg.Types.Context)
    out = PackageID[]
    if VERSION < v"1.1.0-DEV.857"
        for pkg in c.env.manifest
            push!(out, PackageID(pkg[1], get(pkg[2][1], "uuid", "")))
        end
    else
        for pkg in c.env.manifest
            push!(out, PackageID(pkg[2].name, string(pkg[1])))
        end
    end
    return out
end

function can_access(m::Module, s::Symbol)
    try
        return Base.eval(m, :($m.$s))
    catch
        return nothing
    end

end

function find_parent(c, uuid::String, out = Set{PackageID}())
    uuid = typeof(c.env.manifest).parameters[1](uuid)
    for pkg in c.env.manifest
        pkgname = VERSION < v"1.1.0-DEV.857" ? pkg[1] : pkg[2].name
        pkguuid = VERSION < v"1.1.0-DEV.857" ? pkg[2][1]["uuid"] : string(pkg[1])
        if uuid in values(pkg_deps(PackageID(pkgname, pkguuid), c))
            if pkgname in keys(context_deps(c))
                push!(out, PackageID(pkgname, pkguuid))
            else
                find_parent(c, pkguuid, out)
            end
        end
    end
    return out
end

function getgithash(path::String)
    repo = LibGit2.GitRepo(path)
    LibGit2.GitHash(LibGit2.GitObject(repo, LibGit2.name(LibGit2.head(repo))))
end

function isgitrepo(path::String)
    try
        LibGit2.GitRepo(path)
        return true
    catch err
        return false
    end
end


text:
module TestModelAPI

using Test
using MLJBase
import MLJModelInterface
using ..Models
using Distributions
using StableRNGs

rng = StableRNG(661)

@testset "predict_*" begin
    X = rand(rng, 5)
    yfinite = categorical(collect("abaaa"))
    ycont = float.(1:5)

    clf = ConstantClassifier()
    fitresult, _, _ = MLJBase.fit(clf, 1, X, yfinite)
    @test predict_mode(clf, fitresult, X)[1] == 'a'
    @test_throws ArgumentError predict_mean(clf, fitresult, X)
    @test_throws ArgumentError predict_median(clf, fitresult, X)

    rgs = ConstantRegressor()
    fitresult, _, _ = MLJBase.fit(rgs, 1, X, ycont)
    @test predict_mean(rgs, fitresult, X)[1] == 3
    @test predict_median(rgs, fitresult, X)[1] == 3
    @test_throws ArgumentError predict_mode(rgs, fitresult, X)
end

mutable struct UnivariateFiniteFitter <: MLJModelInterface.Probabilistic
    alpha::Float64
end
UnivariateFiniteFitter(;alpha=1.0) = UnivariateFiniteFitter(alpha)

@testset "models that fit a distribution" begin
    function MLJModelInterface.fit(model::UnivariateFiniteFitter,
                               verbosity, X, y)

        α = model.alpha
        N = length(y)
        _classes = classes(y)
        d = length(_classes)

        frequency_given_class = Distributions.countmap(y)
        prob_given_class =
            Dict(c => (frequency_given_class[c] + α)/(N + α*d) for c in _classes)

        fitresult = MLJBase.UnivariateFinite(prob_given_class)

        report = (params=Distributions.params(fitresult),)
        cache = nothing

        verbosity > 0 && @info "Fitted a $fitresult"

        return fitresult, cache, report
    end

    MLJModelInterface.predict(model::UnivariateFiniteFitter,
                              fitresult,
                              X) = fitresult

    MLJModelInterface.input_scitype(::Type{<:UnivariateFiniteFitter}) =
        Nothing
    MLJModelInterface.target_scitype(::Type{<:UnivariateFiniteFitter}) =
        AbstractVector{<:Finite}

    y = coerce(collect("aabbccaa"), Multiclass)
    X = nothing
    model = UnivariateFiniteFitter(alpha=0)
    mach = machine(model, X, y)
    fit!(mach, verbosity=0)

    ytest = y[1:3]
    yhat = predict(mach, nothing) # single UnivariateFinite distribution

    @test cross_entropy(fill(yhat, 3), ytest) ≈
        [-log(1/2), -log(1/2), -log(1/4)]

end

end
true


text:
function write_xdmf_2d(filename::String,
                       HRPM::HierarchicalRectangularlyPartitionedMesh)
    @info "Writing $filename"
    # Check valid filename
    if !occursin(".xdmf", filename)
        @error "Invalid filename. '.xdmf' does not occur in $filename"
    end

    # If there are materials, map all material names to an integer
    material_map = make_material_name_to_id_map(HRPM)

    # h5 filename
    h5_filename = replace(filename, ("xdmf" => "h5"))
    h5_file = h5open(h5_filename, "w")
    # XML
    xdoc = XMLDocument()
    # Xdmf
    xroot = create_root(xdoc, "Xdmf")
    set_attribute(xroot, "Version", "3.0")
    # Domain
    xdomain = new_child(xroot, "Domain")
    # Material names
    if material_map != Dict{String, Int64}()
        xmaterials = new_child(xdomain, "Information")
        set_attribute(xmaterials, "Name", "MaterialNames")
        add_text(xmaterials, join(keys(material_map), " "))
    end
    add_HRPM_xdmf!(xdomain, h5_filename, h5_file, HRPM, material_map)
    save_file(xdoc, filename)
    close(h5_file)
end

function write_xdmf_2d(filename::String, mesh::UnstructuredMesh_2D)
    @info "Writing $filename"
    # Check valid filename
    if !occursin(".xdmf", filename)
        @error "Invalid filename. '.xdmf' does not occur in $filename"
    end

    # If there are materials, map all material names to an integer
    material_map = Dict{String, Int64}()
    if mesh.face_sets != Dict{String, Set{UInt32}}()
        material_map = make_material_name_to_id_map(mesh)
    end

    # h5 filename
    h5_filename = replace(filename, ("xdmf" => "h5"))
    h5_file = h5open(h5_filename, "w")
    # XML
    xdoc = XMLDocument()
    # Xdmf
    xroot = create_root(xdoc, "Xdmf")
    set_attribute(xroot, "Version", "3.0")
    # Domain
    xdomain = new_child(xroot, "Domain")
    # Material names
    if material_map != Dict{String, Int64}()
        xmaterials = new_child(xdomain, "Information")
        set_attribute(xmaterials, "Name", "MaterialNames")
        add_text(xmaterials, join(keys(material_map), " "))
    end

    # Add uniform grid
    add_uniform_grid_xdmf!(xdomain, h5_filename, h5_file, mesh, material_map)

    save_file(xdoc, filename)
    close(h5_file)
end

# Helper functions for write_xdmf
# -------------------------------------------------------------------------------------------------
function add_uniform_grid_xdmf!(xml::XMLElement,
                                h5_filename::String,
                                h5_mesh::Union{HDF5.Group, HDF5.File},
                                mesh::UnstructuredMesh_2D,
                                material_map::Dict{String, Int64})
    @debug "Adding uniform grid for $(mesh.name)"
    # Grid
    xgrid = new_child(xml, "Grid")
    set_attribute(xgrid, "Name", mesh.name)
    set_attribute(xgrid, "GridType", "Uniform")

    # h5
    h5_group = create_group(h5_mesh, mesh.name)

    # Geometry
    write_xdmf_geometry!(xgrid, h5_filename, h5_group, mesh)

    # Topology
    write_xdmf_topology!(xgrid, h5_filename, h5_group, mesh)

    # Non-material face sets
    if mesh.face_sets != Dict{String, Set{UInt64}}()
        write_xdmf_face_sets!(xgrid, h5_filename, h5_group, mesh)
    end

    # Materials
    if material_map != Dict{String, Int64}()
        write_xdmf_materials!(xgrid, h5_filename, h5_group, mesh, material_map)
    end
    return nothing
end

function write_xdmf_geometry!(xml::XMLElement,
                              h5_filename::String,
                              h5_mesh::HDF5.Group,
                              mesh::UnstructuredMesh_2D)
    @debug "Writing XDMF geometry"
    # Geometry
    xgeom = new_child(xml, "Geometry")
    set_attribute(xgeom, "GeometryType", "XYZ")
    # DataItem
    xdataitem = new_child(xgeom, "DataItem")
    set_attribute(xdataitem, "DataType", "Float")
    npoints = length(mesh.points)
    set_attribute(xdataitem, "Dimensions", "$npoints 3")
    set_attribute(xdataitem, "Format", "HDF")
    float_precision = 8
    set_attribute(xdataitem, "Precision", "$float_precision")
    h5_text_item = split(string(h5_mesh))[2]
    add_text(xdataitem, string(h5_filename, ":", h5_text_item, "/points"))
    # Convert the points into an array
    point_array = zeros(3, npoints)
    for i = 1:npoints
        point_array[1, i] = mesh.points[i][1]
        point_array[2, i] = mesh.points[i][2]
    end
    # Write the h5
    h5_mesh["points"] = point_array
    return nothing
end

function write_xdmf_topology!(xml::XMLElement,
                              h5_filename::String,
                              h5_mesh::HDF5.Group,
                              mesh::UnstructuredMesh_2D)
    @debug "Writing XDMF topology"
    # Topology
    xtopo = new_child(xml, "Topology")
    set_attribute(xtopo, "TopologyType", "Mixed")
    nelements = length(mesh.faces)
    set_attribute(xtopo, "NumberOfElements", "$nelements")
    # DataItem
    xdataitem = new_child(xtopo, "DataItem")
    set_attribute(xdataitem, "DataType", "Int")
    topo_length = mapreduce(x->length(x), +, mesh.faces) + nelements
    topo_array = Vector{UInt64}(undef, topo_length)
    convert_xdmf_faces_to_array!(topo_array, mesh.faces)
    ndimensions = length(topo_array)
    set_attribute(xdataitem, "Dimensions", "$ndimensions")
    set_attribute(xdataitem, "Format", "HDF")
    set_attribute(xdataitem, "Precision", "8")
    h5_text_item = split(string(h5_mesh))[2]
    add_text(xdataitem, string(h5_filename, ":", h5_text_item, "/cells"))
    # Write the h5
    h5_mesh["cells"] = topo_array
    return nothing
end

function convert_xdmf_faces_to_array!(topo_array::Vector{UInt64},
                                      faces::Vector{<:SArray{S, UInt32, 1, L} where {S<:Tuple, L}})
#    length_to_xdmf_type = Dict{UInt32, UInt32}(
#        # triangle
#        3  => 4,
#        # triangle6
#        6 => 36,
#        # quadrilateral
#        4 => 5,
#        # quad8
#        8 => 37
#       )
    topo_ctr = 1
    for face in faces
        # convert face to vector for mutability
        # adjust 1-based to 0-based indexing
        face_xdmf = collect(face) .- 1
        # add xdmf type
        face_length = length(face)
        if face_length === 3
            pushfirst!(face_xdmf, 4)
        elseif face_length === 4
            pushfirst!(face_xdmf, 5)
        elseif face_length === 6
            pushfirst!(face_xdmf, 36)
        elseif face_length === 8
            pushfirst!(face_xdmf, 37)
        else
            @error "Unknown face type"
        end
        topo_array[topo_ctr:topo_ctr + face_length] = face_xdmf
        topo_ctr += face_length + 1
    end
    return nothing
end

function make_material_name_to_id_map(mesh::UnstructuredMesh_2D)
    material_map = Dict{String, Int64}()
    mat_names = String[]
    max_length = 0
    for set_name in keys(mesh.face_sets)
        if occursin("MATERIAL", uppercase(set_name))
            material_map[set_name] = 0
            push!(mat_names, set_name)
            if length(set_name) > max_length
                max_length = length(set_name)
            end
        end
    end
    # Adjust material IDs so that IDs are assigned alphabetically
    sort!(mat_names)
    for (i, set_name) in enumerate(mat_names)
        material_map[set_name] = i - 1
    end
    if max_length < 13
        max_length = 13
    end
    println(string(rpad("Material Name", max_length, ' '), " : XDMF Material ID"))
    println(rpad("=", max_length + 19, '='))
    for set_name in mat_names
        if occursin("MATERIAL", uppercase(set_name))
            id = material_map[set_name]
            println(string(rpad(set_name, max_length, ' '), " : $id"))
        end
    end
    return material_map
end

function make_material_name_to_id_map(HRPM::HierarchicalRectangularlyPartitionedMesh)
    mesh_children = [HRPM]
    next_mesh_children = HierarchicalRectangularlyPartitionedMesh[]
    leaves_reached = false
    while !leaves_reached
        for child_mesh in mesh_children
            if length(child_mesh.children) > 0
                for child_ref in child_mesh.children
                    push!(next_mesh_children, child_ref[])
                end
            else
                leaves_reached = true
            end
        end
        if !leaves_reached
            mesh_children = next_mesh_children
            next_mesh_children = HierarchicalRectangularlyPartitionedMesh[]
        end
    end
    material_map = Dict{String, Int64}()
    mat_names = String[]
    max_length = 0
    for leaf_mesh in mesh_children
        UM = leaf_mesh.mesh[]
        for set_name in keys(UM.face_sets)
            if occursin("MATERIAL", uppercase(set_name))
                if set_name ∉  keys(material_map)
                    material_map[set_name] = 0
                    push!(mat_names, set_name)
                    if length(set_name) > max_length
                        max_length = length(set_name)
                    end
                end
            end
        end
    end
    # Adjust material IDs so that IDs are assigned alphabetically
    sort!(mat_names)
    for (i, set_name) in enumerate(mat_names)
        material_map[set_name] = i - 1
    end
    if max_length < 13
        max_length = 13
    end
    println(string(rpad("Material Name", max_length, ' '), " : XDMF Material ID"))
    println(rpad("=", max_length + 19, '='))
    for set_name in mat_names
        if occursin("MATERIAL", uppercase(set_name))
            id = material_map[set_name]
            println(string(rpad(set_name, max_length, ' '), " : $id"))
        end
    end
    return material_map
end

function write_xdmf_materials!(xml::XMLElement,
                               h5_filename::String,
                               h5_mesh::HDF5.Group,
                               mesh::UnstructuredMesh_2D,
                               material_map::Dict{String, Int64})
    @debug "Writing XDMF materials"
    # MaterialID
    xmaterial = new_child(xml, "Attribute")
    set_attribute(xmaterial, "Center", "Cell")
    set_attribute(xmaterial, "Name", "MaterialID")
    nelements = length(mesh.faces)
    mat_ID_array = fill(-1, nelements)
    for material_name in keys(material_map)
        if material_name ∈  keys(mesh.face_sets)
            material_ID = material_map[material_name]
            for cell in mesh.face_sets[material_name]
                if mat_ID_array[cell] === -1
                    mat_ID_array[cell] = material_ID
                else
                    @error "Mesh cell $cell has multiple materials assigned to it."
                end
            end
        end
    end
    if any(x->x === -1, mat_ID_array)
        @error "Some mesh cells do not have a material."
    end
    # DataItem
    xdataitem = new_child(xmaterial, "DataItem")
    set_attribute(xdataitem, "DataType", "Int")
    set_attribute(xdataitem, "Dimensions", "$nelements")
    set_attribute(xdataitem, "Format", "HDF")
    set_attribute(xdataitem, "Precision", "8")
    h5_text_item = split(string(h5_mesh))[2]
    add_text(xdataitem, string(h5_filename, ":", h5_text_item, "/material_id"))
    # Write the h5
    h5_mesh["material_id"] = mat_ID_array
    return nothing
end

function write_xdmf_face_sets!(xml::XMLElement,
                               h5_filename::String,
                               h5_mesh::HDF5.Group,
                               mesh::UnstructuredMesh_2D)
    @debug "Writing XDMF face_sets"
    for set_name in keys(mesh.face_sets)
        if occursin("MATERIAL", uppercase(set_name))
            continue
        end
        # Set
        xset = new_child(xml, "Set")
        set_attribute(xset, "Name", set_name)
        set_attribute(xset, "SetType", "Cell")
        # DataItem
        xdataitem = new_child(xset, "DataItem")
        set_attribute(xdataitem, "DataType", "Int")
        nelements = length(mesh.face_sets[set_name])
        set_attribute(xdataitem, "Dimensions", "$nelements")
        set_attribute(xdataitem, "Format", "HDF")
        set_attribute(xdataitem, "Precision", "8")
        h5_text_item = split(string(h5_mesh))[2]
        add_text(xdataitem, string(h5_filename, ":", h5_text_item, "/", set_name))
        # Write the h5
        ID_array = collect(mesh.face_sets[set_name]) .- 1
        h5_mesh[set_name] = ID_array
    end
    return nothing
end

function add_HRPM_xdmf!(xml::XMLElement,
                        h5_filename::String,
                        h5_mesh::Union{HDF5.Group, HDF5.File},
                        HRPM::HierarchicalRectangularlyPartitionedMesh,
                        material_map::Dict{String, Int64})
    @debug "Adding HRPM for $(HRPM.name)"
    if length(HRPM.children) > 0
        # Grid
        xgrid = new_child(xml, "Grid")
        set_attribute(xgrid, "Name", HRPM.name)
        set_attribute(xgrid, "GridType", "Tree")
        # h5_group
        h5_group = create_group(h5_mesh, HRPM.name)
        for child in HRPM.children
            add_HRPM_xdmf!(xgrid, h5_filename, h5_group, child[], material_map)
        end
    else
        add_uniform_grid_xdmf!(xml, h5_filename, h5_mesh, HRPM.mesh[], material_map)
    end
end

max_stars_repo_name: KyleVaughn/MOCNeutronTransport

max_issues_repo_name: KyleVaughn/MOCNeutronTransport

max_forks_repo_name: KyleVaughn/MOCNeutronTransport


text:
"""
    summarize(doc [, ns])

Summarizes the document and returns `ns` number of sentences.
By default `ns` is set to the value 5.

# Example
```julia-repl
julia> s = StringDocument("Assume this Short Document as an example. Assume this as an example summarizer. This has too foo sentences.")

julia> summarize(s, ns=2)
2-element Array{SubString{String},1}:
 "Assume this Short Document as an example."
 "This has too foo sentences."
```
"""
function summarize(d::AbstractDocument; ns=5)
    sentences = sentence_tokenize(language(d), text(d))
    num_sentences = length(sentences)
    s = StringDocument.(sentences)
    c = Corpus(s)
    prepare!(c, strip_case | strip_stopwords | stem_words )
    update_lexicon!(c)
    t = tf_idf(dtm(c))
    T = t * t'
    p = pagerank(T)
    return sentences[sort(sortperm(vec(p), rev=true)[1:min(ns, num_sentences)])]
end

function pagerank( A; Niter=20, damping=.15)
         Nmax = size(A, 1)
         r = rand(1,Nmax);              # Generate a random starting rank.
         r = r ./ norm(r,1);            # Normalize
         a = (1-damping) ./ Nmax;       # Create damping vector

         for i=1:Niter
             s = r * A
             rmul!(s, damping)
             r = s .+ (a * sum(r, dims=2));   # Compute PageRank.
         end

         r = r./norm(r,1);

         return r
end


text:
using LibPQ, Tables, DelimitedFiles
conn = LibPQ.Connection("dbname=postgres host=127.0.0.1 user=postgres password=test port=5432")
execute(conn, "create table if not exists hoge (id int)")
for i in 1:10
    execute(conn, "insert into hoge values ($i)")
end
result = execute(conn, "SELECT * FROM hoge")
writedlm("test.csv", result)


text:
#= JuTeX types =#

#= Basic types =#
const Boolean = Bool

"size_t"
const Pointer = UInt
const UChar = UInt8

#= [#18] Types in the outer block
## Sections
    18, 25, 38, 101, 109,
    113, 150, 212, 269, 300,
    [548], [594], 920, 925.

## All types
    #18: ASCIICode::UChar(8)
    #25: EightBits::UChar(8)
    #38: PackedASCIICode::UChar(8)
    #101: Scaled::Integer(64)
    #101: NonNegativeInteger::UInteger(64)
    #101: SmallNumber::UChar(8)
    #109: GlueRatio::Real(64)
    #113: QuarterWord(8=>16)
    #113: HalfWord(16=>32)
    #113: TwoHalves(64)
    #113: FourQuarters(64)
    #113: MemoryWord(64)
    #150: GlueOrd(8)
    #212: ListStateRecord(366)
    #269: GroupCode(8)
    #300: InStateRecord(1344)
    #920: TriePointer::Pointer(64)
    #925: HyphPointer::Pointer(64)
=#

@group "S25x37_P13x18" begin
const EightBits = missing
end # S25x37_P13x18

const PackedASCIICode = missing

""" [p38#101]: `Scaled::Int(64)` = [−2^{64}, 2^{64})
"""
const Scaled = Int64

const NonNegativeInteger = missing
""" [p38#101]: `SmallNumber::UChar(8)` = [0, 255]
only need: (6bit)[0, 63] .
"""
const SmallNumber = UChar
const GlueRatio = missing

# S110x114_P42x43
@group "S110x114_P42x43" begin
const QuarterWord = missing
const HalfWord = missing

struct TwoHalves end
struct FourQuarters end
struct MemoryWord end
end # S110x114_P42x43

struct GlueOrd end
struct GroupCode end
struct TriePointer end
struct HyphPointer end


text:
# Utilities

##### Buffer

struct XBuffer
    ptr::Xptr

    function XBuffer(bytes::Integer)
        p = ccall((:xmlBufferCreateSize,libxml2), Xptr, (Csize_t,), bytes)
        p != C_NULL || error("Failed to create buffer of $bytes bytes.")
        new(p)
    end
end

free(buf::XBuffer) = ccall((:xmlBufferFree,libxml2), Cvoid, (Xptr,), buf.ptr)

Base.length(buf::XBuffer) = int(ccall((:xmlBufferLength,libxml2), Cint, (Xptr,), buf.ptr))

content(buf::XBuffer) = unsafe_string(ccall((:xmlBufferContent,libxml2), Xstr, (Xptr,), buf.ptr))


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/MortarContact2D.jl/blob/master/LICENSE

mutable struct Contact2D <: BoundaryProblem
    master_elements :: Vector{Element}
    rotate_normals :: Bool
    dual_basis :: Bool
    max_distance :: Float64
    iteration :: Int
    contact_state_in_first_iteration :: Symbol
    store_fields :: Vector{String}
end

function Contact2D()
    return Contact2D([], false, true, Inf, 0, :AUTO, [])
end

function FEMBase.add_elements!(::Problem{Contact2D}, ::Any)
    error("use `add_slave_elements!` and `add_master_elements!` to add ",
          "elements to the Contact2D problem.")
end

function FEMBase.get_elements(problem::Problem{Contact2D})
    return [problem.elements; problem.properties.master_elements]
end

function FEMBase.add_slave_elements!(problem::Problem{Contact2D}, elements)
    for element in elements
        push!(problem.elements, element)
    end
end

function FEMBase.add_master_elements!(problem::Problem{Contact2D}, elements)
    for element in elements
        push!(problem.properties.master_elements, element)
    end
end

function FEMBase.get_slave_elements(problem::Problem{Contact2D})
    return problem.elements
end

function FEMBase.get_master_elements(problem::Problem{Contact2D})
    return problem.properties.master_elements
end

function get_slave_dofs(problem::Problem{Contact2D})
    dofs = Int64[]
    for element in get_slave_elements(problem)
        append!(dofs, get_gdofs(problem, element))
    end
    return sort(unique(dofs))
end

function get_master_dofs(problem::Problem{Contact2D})
    dofs = Int64[]
    for element in get_master_elements(problem)
        append!(dofs, get_gdofs(problem, element))
    end
    return sort(unique(dofs))
end

function create_rotation_matrix(element::Element{Seg2}, time::Float64)
    n = element("normal", time)
    R = [0.0 -1.0; 1.0 0.0]
    t1 = R'*n[1]
    t2 = R'*n[2]
    Q1 = [n[1] t1]
    Q2 = [n[2] t2]
    Z = zeros(2, 2)
    Q = [Q1 Z; Z Q2]
    return Q
end

function create_contact_segmentation(problem::Problem{Contact2D}, slave_element::Element{Seg2}, master_elements::Vector, time::Float64; deformed=false)

    result = []

    max_distance = problem.properties.max_distance

    x1 = slave_element("geometry", time)

    if deformed && haskey(slave_element, "displacement")
        x1 = map(+, x1, slave_element("displacement", time))
    end

    slave_midpoint = 1/2*(x1[1] + x1[2])
    slave_length = norm(x1[2]-x1[1])

    for master_element in master_elements

        x2 = master_element("geometry", time)

        if deformed && haskey(master_element, "displacement")
            x2 = map(+, x2, master_element("displacement", time))
        end

        master_midpoint = 1/2*(x2[1] + x2[2])
        master_length = norm(x2[2]-x2[1])
        # charasteristic length
        dist = norm(slave_midpoint - master_midpoint)
        cl = dist/max(slave_length, master_length)
        if cl > max_distance
            continue
        end

        # 3.1 calculate segmentation
        x21, x22 = x2
        xi1a = project_from_master_to_slave(slave_element, x21, time)
        xi1b = project_from_master_to_slave(slave_element, x22, time)
        xi1 = clamp.([xi1a; xi1b], -1.0, 1.0)
        l = 1/2*abs(xi1[2]-xi1[1])
        if isapprox(l, 0.0)
            continue # no contribution in this master element
        end
        push!(result, (master_element, xi1, l))
    end
    return result
end

function FEMBase.assemble_elements!(problem::Problem{Contact2D}, assembly::Assembly,
                                    elements::Vector{Element{Seg2}}, time::Float64)

    props = problem.properties
    props.iteration += 1

    slave_elements = get_slave_elements(problem)

    field_dim = get_unknown_field_dimension(problem)
    field_name = get_parent_field_name(problem)

    # 1. calculate nodal normals for slave element nodes j ∈ S
    normals = calculate_normals([slave_elements...], time; rotate_normals=props.rotate_normals)
    tangents = Dict(j => [t2, -t1] for (j, (t1,t2)) in normals)
    update!(slave_elements, "normal", time => normals)
    update!(slave_elements, "tangent", time => tangents)

    Rn = 0.0

    # 2. loop all slave elements
    for slave_element in slave_elements

        nsl = length(slave_element)
        X1 = slave_element("geometry", time)
        x1 = slave_element("geometry", time)
        if haskey(slave_element, "displacement")
            u1 = slave_element("displacement", time)
            x1 = map(+, x1, u1)
        end
        la1 = slave_element("lambda", time)
        n1 = slave_element("normal", time)
        t1 = slave_element("tangent", time)

        contact_area = 0.0
        contact_error = 0.0
        Q2 = create_rotation_matrix(slave_element, time)

        master_elements = get_master_elements(problem)
        segmentation = create_contact_segmentation(problem, slave_element, master_elements, time)
        if length(segmentation) == 0 # no overlapping in master and slave surfaces with this slave element
            continue
        end

        Ae = Matrix(I, nsl, nsl)
        if props.dual_basis
            De = zeros(nsl, nsl)
            Me = zeros(nsl, nsl)
            for (master_element, xi1, l) in segmentation
                for ip in get_integration_points(slave_element, 3)
                    detJ = slave_element(ip, time, Val{:detJ})
                    w = ip.weight*detJ*l
                    xi = ip.coords[1]
                    xi_s = dot([1/2*(1-xi); 1/2*(1+xi)], xi1)
                    N1 = vec(get_basis(slave_element, xi_s, time))
                    De += w*Matrix(Diagonal(N1))
                    Me += w*N1*N1'
                end
                Ae = De*inv(Me)
            end
        end

        # loop all segments
        for (master_element, xi1, l) in segmentation

            nm = length(master_element)
            X2 = master_element("geometry", time)
            x2 = master_element("geometry", time)
            if haskey(master_element, "displacement")
                u2 = master_element("displacement", time)
                x2 = map(+, x2, u2)
            end

            # 3.3. loop integration points of one integration segment and calculate
            # local mortar matrices
            De = zeros(nsl, nsl)
            Me = zeros(nsl, nsl)
            Ne = zeros(nsl, 2*nsl)
            Te = zeros(nsl, 2*nsl)
            He = zeros(nsl, 2*nsl)
            ce = zeros(nsl)
            ge = zeros(nsl)
            for ip in get_integration_points(slave_element, 3)
                detJ = slave_element(ip, time, Val{:detJ})
                w = ip.weight*detJ*l
                xi = ip.coords[1]
                xi_s = dot([1/2*(1-xi); 1/2*(1+xi)], xi1)
                N1 = vec(get_basis(slave_element, xi_s, time))
                Phi = Ae*N1

                # project gauss point from slave element to master element in direction n_s
                X_s = interpolate(N1, X1) # coordinate in gauss point
                n_s = interpolate(N1, n1) # normal direction in gauss point
                t_s = interpolate(N1, t1) # tangent condition in gauss point
                n_s /= norm(n_s)
                t_s /= norm(t_s)
                xi_m = project_from_slave_to_master(master_element, X_s, n_s, time)
                N2 = vec(get_basis(master_element, xi_m, time))
                X_m = interpolate(N2, X2)

                u_s = interpolate(N1, u1)
                u_m = interpolate(N2, u2)
                x_s = map(+, X_s, u_s)
                x_m = map(+, X_m, u_m)
                la_s = interpolate(Phi, la1)

                # contact virtual work
                De += w*Phi*N1'
                Me += w*Phi*N2'

                # contact constraints
                Ne += w*reshape(kron(N1, n_s, Phi), 2, 4)
                Te += w*reshape(kron(N2, n_s, Phi), 2, 4)
                He += w*reshape(kron(N1, t_s, Phi), 2, 4)
                ge += w*Phi*dot(n_s, x_m-x_s)
                ce += w*N1*dot(n_s, la_s)
                Rn += w*dot(n_s, la_s)

                contact_area += w
                contact_error += 1/2*w*dot(n_s, x_s-x_m)^2
            end

            sdofs = get_gdofs(problem, slave_element)
            mdofs = get_gdofs(problem, master_element)

            # add contribution to contact virtual work
            for i=1:field_dim
                lsdofs = sdofs[i:field_dim:end]
                lmdofs = mdofs[i:field_dim:end]
                add!(problem.assembly.C1, lsdofs, lsdofs, De)
                add!(problem.assembly.C1, lsdofs, lmdofs, -Me)
            end

            # add contribution to contact constraints
            add!(problem.assembly.C2, sdofs[1:field_dim:end], sdofs, Ne)
            add!(problem.assembly.C2, sdofs[1:field_dim:end], mdofs, -Te)
            add!(problem.assembly.D, sdofs[2:field_dim:end], sdofs, He)
            add!(problem.assembly.g, sdofs[1:field_dim:end], ge)
            add!(problem.assembly.c, sdofs[1:field_dim:end], ce)

        end # master elements done

        if "contact area" in props.store_fields
            update!(slave_element, "contact area", time => contact_area)
        end

        if "contact error" in props.store_fields
            update!(slave_element, "contact error", time => contact_error)
        end

    end # slave elements done, contact virtual work ready

    S = sort(collect(keys(normals))) # slave element nodes
    weighted_gap = Dict{Int64, Vector{Float64}}()
    contact_pressure = Dict{Int64, Vector{Float64}}()
    complementarity_condition = Dict{Int64, Vector{Float64}}()
    is_active = Dict{Int64, Int}()
    is_inactive = Dict{Int64, Int}()
    is_slip = Dict{Int64, Int}()
    is_stick = Dict{Int64, Int}()

    la = problem.assembly.la
    # FIXME: for matrix operations, we need to know the dimensions of the
    # final matrices
    ndofs = 2*length(S)
    ndofs = max(ndofs, length(la))
    ndofs = max(ndofs, size(problem.assembly.K, 2))
    ndofs = max(ndofs, size(problem.assembly.C1, 2))
    ndofs = max(ndofs, size(problem.assembly.C2, 2))
    ndofs = max(ndofs, size(problem.assembly.D, 2))
    ndofs = max(ndofs, size(problem.assembly.g, 2))
    ndofs = max(ndofs, size(problem.assembly.c, 2))

    C1 = sparse(problem.assembly.C1, ndofs, ndofs)
    C2 = sparse(problem.assembly.C2, ndofs, ndofs)
    D = sparse(problem.assembly.D, ndofs, ndofs)
    g = Vector(sparse(problem.assembly.g, ndofs, 1)[:])
    c = Vector(sparse(problem.assembly.c, ndofs, 1)[:])

    for j in S
        dofs = [2*(j-1)+1, 2*(j-1)+2]
        weighted_gap[j] = g[dofs]
    end

    state = problem.properties.contact_state_in_first_iteration
    if problem.properties.iteration == 1
        @info("First contact iteration, initial contact state = $state")

        if state == :AUTO
            avg_gap = mean([weighted_gap[j][1] for j in S])
            std_gap = std([weighted_gap[j][1] for j in S])
            if (avg_gap < 1.0e-12) && (std_gap < 1.0e-12)
                state = :ACTIVE
            else
                state = :UNKNOWN
            end
            @info("Average weighted gap = $avg_gap, std gap = $std_gap, automatically determined contact state = $state")
        end

    end

    # active / inactive node detection
    for j in S
        dofs = [2*(j-1)+1, 2*(j-1)+2]
        weighted_gap[j] = g[dofs]

        if length(la) != 0
            p = dot(normals[j], la[dofs])
            t = dot(tangents[j], la[dofs])
            contact_pressure[j] = [p, t]
        else
            contact_pressure[j] = [0.0, 0.0]
        end

        complementarity_condition[j] = contact_pressure[j] - weighted_gap[j]
        if complementarity_condition[j][1] < 0
            is_inactive[j] = 1
            is_active[j] = 0
            is_slip[j] = 0
            is_stick[j] = 0
        else
            is_inactive[j] = 0
            is_active[j] = 1
            is_slip[j] = 1
            is_stick[j] = 0
        end
    end

    if (problem.properties.iteration == 1) && (state == :ACTIVE)
        for j in S
            is_inactive[j] = 0
            is_active[j] = 1
            is_slip[j] = 1
            is_stick[j] = 0
        end
    end

    if (problem.properties.iteration == 1) && (state == :INACTIVE)
        for j in S
            is_inactive[j] = 1
            is_active[j] = 0
            is_slip[j] = 0
            is_stick[j] = 0
        end
    end

    if "weighted gap" in props.store_fields
        update!(slave_elements, "weighted gap", time => weighted_gap)
    end
    if "contact pressure" in props.store_fields
        update!(slave_elements, "contact pressure", time => contact_pressure)
    end
    if "complementarity condition" in props.store_fields
        update!(slave_elements, "complementarity condition", time => complementarity_condition)
    end
    if "active nodes" in props.store_fields
        update!(slave_elements, "active nodes", time => is_active)
    end
    if "inactive nodes" in props.store_fields
        update!(slave_elements, "inactive nodes", time => is_inactive)
    end
    if "stick nodes" in props.store_fields
        update!(slave_elements, "stick nodes", time => is_stick)
    end
    if "slip nodes" in props.store_fields
        update!(slave_elements, "slip nodes", time => is_slip)
    end

    @info("# | A | I | St | Sl | gap | pres | comp")
    for j in S
        str1 = "$j | $(is_active[j]) | $(is_inactive[j]) |  $(is_stick[j]) |  $(is_slip[j]) | "
        str2 = "$(round(weighted_gap[j][1]; digits=3)) | $(round(contact_pressure[j][1]; digits=3)) | $(round(complementarity_condition[j][1]; digits=3))"
        @info(str1 * str2)
    end

    # solve variational inequality

    # remove inactive nodes from assembly
    for j in S
        dofs = [2*(j-1)+1, 2*(j-1)+2]
        n, t = dofs
        if is_inactive[j] == 1
            @debug("$j is inactive, removing dofs $dofs")
            C1[dofs,:] .= 0.0
            C2[dofs,:] .= 0.0
            D[dofs,:] .= 0.0
            g[dofs,:] .= 0.0
        elseif (is_active[j] == 1) && (is_slip[j] == 1)
            @debug("$j is in active/slip, removing tangential constraint $t")
            C2[t,:] .= 0.0
            g[t] = 0.0
            D[t, dofs] .= tangents[j]
        end
    end

    problem.assembly.C1 = C1
    problem.assembly.C2 = C2
    problem.assembly.D = D
    problem.assembly.g = sparsevec(g)

end


text:
module AutoLandmarking

include("data_augmentation.jl")
include("data_preparation.jl")
include("handling_and_utility.jl")
include("models.jl")
include("outliers.jl")
include("custom_data_loading.jl")
# Write your package code here.
export
    # data loading
    load_imgs,
    read_landmarks,

    #data augmentation
    flip_volume_side,
    flip_volume_front,
    mirror_vol,
    flip_2D,
    flip_3D,
    jitter_2D,
    jitter_3D,
    rotate_images,
    rotate_volumes,

    #utility and data preparation
    swap_xy,
    empty_3rd_dimension,
    resize_2d_images,
    train_test_split_3d,
    regular_train_test_split_3d,
    regular_train_test_split_2d,
    landmark_to_surface,
    image_gradients,
    align_principal,
    resize_relevant,
    translate_lms_back,
    choose_dims,
    depth_map_all_sides,

    #evaluation and training
    cost_whole_data_2D,
    cost_whole_data_3D,
    run_model,
    avg_accuracy_per_point,
    predict_single,
    predict_set,

    #confidence testing and exporting
    save_vols_to_folder,
    response_distribution,
    procrustes_distance_list,
    writedlm

end

max_stars_repo_name: imi-diener/AutoLM

max_issues_repo_name: imi-diener/AutoLM

max_forks_repo_name: imi-diener/AutoLM


text:
using Parameters, Combinatorics, Distributions

"""
The purpose of this module is to define a flexible stochastic blockmodel for hypergraphs. At this stage, all we are aiming to do is *sample* from the model give a partition, a group intensity function, and a degree vector.
"""

Parameters.@with_kw mutable struct hypergraph
    """
    A very simple hypergraph composite type, designed to hold a node list N, an edge list E, a degree sequence D,
    """

    N::Vector{Int64}
    E::Dict{Int64, Dict}
    D::Array{Int64, 1} = Array{Int64, 1}()

end

function sampleEdge(S::Vector{Int64}, Z::Vector{Int64}, ϑ::Vector{Float64}, Ω::IntensityFunction; α::Any)
    """
    Sample a Poisson number of edges on a sequence of nodes S according to the hypergraph SBM law.
    # Arguments
    S::Vector{Int64},  an array of node indices
    Z::Vector{Int64}, an integer array of length n. Each integer is a group label.
    Ω::IntensityFunction, an IntensityFunction governing interaction rates between nodes of different groups.
    ϑ::Vector{Float64},  a nonnegative float array of length n
    # Returns
    x::Int64, an integer number of edge counts sampled according to a Poisson random variable with specified formula.
    """
    z = Z[S]
    θ = ϑ[S]

    # combinatorial factor associated with repeated indices. Equivalent to sampling a separate Poisson for each of the c distinct permutations of the node labels S
    c = counting_coefficient(S)
    X = Distributions.Poisson(prod(θ)*Ω.ω(Ω.P(z), α)*c)
    return(rand(X))
end

function sampleEdges(Z::Vector{Int64}, ϑ::Vector{Float64}, Ω::IntensityFunction; α::Vector{Float64}, kmax::Int64=3, kmin::Int64=2)
    """
    Sample edges of specified sizes according to the degree-corrected hypergraph configuration model (DCHSBM). Runs sampleEdge() for each possible distinct sequence of nodes of size between kmin and kmax, inclusive.
    # Arguments
    Z::Vector{Int64}, a vector of group labels
    ϑ::Vector{Float64}, a vector of degree parameters
    Ω::IntensityFunction, an intensity function governing the rate of interaction between nodes of various groups.
    α::Vector{Float64}, a parameter passed to Ω
    kmax::Int64, the size of the largest hyperedge to form
    kmin::Int64, the size of the smallest hyperedge to form
    # Returns
    E::Dict{Int64,Dict{Vector{Int64}, Int64}}, a dictionary keyed by edge sizes. E[k] is itself a dictionary whose keys are ordered sets of nodes and whose values are the number of edges on those node sets.
    N::Vector{Int64}, the list of nodes.
    """
    E = Dict{Int64, Dict}() # initialize empty edge list
    n = length(Z)           # number of nodes inferred from group labels

    for k in kmin:kmax
        T = Combinatorics.with_replacement_combinations(1:n, k) # all distinct combos of k nodes
        Ek = Dict{Vector{Int64}, Int64}()         # list of edges of size k
        for S in T
            X = sampleEdge(S, Z, ϑ, Ω; α=α)
            if X > 0                              # only store combos with at least one edge
                Ek[S] = X
            end
        end
        E[k] = Ek
    end
    N = 1:n
    return(E, N)
end

function sampleEdges(Z::Dict, ϑ::Dict, Ω::IntensityFunction; α::Any, kmax::Int64=3, kmin::Int64=2)
    """
    An alternative method when Z and ϑ are specified as dicts keyed by node rather than as arrays.
    """
    Z = [Z[i] for i in 1:length(Z)]
    ϑ = [ϑ[i] for i in 1:length(ϑ)]
    sampleEdges(Z, ϑ, Ω; α=α, kmax=kmax, kmin=kmin)
end

function computeDegrees(E::Dict{Int64, Dict}, N::Vector{Int64})
    """
    Compute the degree sequence of an edge list.
    """

    d = zeros(length(N))

    for k in keys(E)
        Ek = E[k]
        for e in keys(Ek)
            for i in e
                d[i] += 1
            end
        end
    end
    return(d)
end

function computeDegrees(H::hypergraph)
    return computeDegrees(H.E, H.N)
end

function computeDegrees!(H::hypergraph)
    """
    Compute the degree sequence of a hypergraph and store it as a field of the hypergraph.
    """
    H.D = computeDegrees(H)
end

function sampleSBM(args...;kwargs...)
    """
    Sample a hypergraph with specified parameters and return it with its degree sequence pre-computed. This is the primary user-facing function for sampling tasks.
    # Arguments
    args and kwargs: See documentation for sampleEdges()
    # Returns
    H::hypergraph, a hypergraph sampled according to the DCHSBM with specified args and kwargs.
    """
    E, N = sampleEdges(args...;kwargs...)
    H = hypergraph(E = E, N = N)

    computeDegrees!(H)
    return(H)
end

function countEdges(H::hypergraph)
    """
    count the number of edges in H
    """
    sum([length(H.E[k]) for k in keys(H.E)])
end

Base.copy(H::hypergraph) = hypergraph(H.N, H.E, H.D)


text:
import TexTables:   get_level, generate_schema, check_table_type,
                    default_sep, schema_lookup

@testset "Extracting Index Levels" begin
    # Standard Composite Tables
    Random.seed!(1234)
    x  = randn(10)
    y  = [Symbol(:key, i) for i=1:10]
    t1 = TableCol("test", y, x)
    t2 = TableCol("test2", y[2:9], x[2:9])
    t3 = TableCol("test3", y, x, randn(10) .|> abs .|> sqrt)
    t4 = TableCol("test" , Dict("Fixed Effects"=>"Yes"))
    t5 = TableCol("test2", Dict("Fixed Effects"=>"No"))
    t6 = TableCol("test3", Dict("Fixed Effects"=>"Yes"))
    c1 = append_table(t1, t4)
    c2 = append_table(t2, t5)
    c3 = append_table(t3, t6)

    # Check that the index levels are what I think they should be:
    @test map(x->get_level(x, 1), c1.row_index) == begin
        Tuple{Int64,Symbol}[(1, Symbol("")), (1, Symbol("")),
                            (1, Symbol("")), (1, Symbol("")),
                            (1, Symbol("")), (1, Symbol("")),
                            (1, Symbol("")), (1, Symbol("")),
                            (1, Symbol("")), (1, Symbol("")),
                            (2, Symbol(""))]
    end

    @test map(x->get_level(x, 2), c1.row_index) == begin
        Tuple{Int64,Symbol}[(1, :key1), (2, :key2), (3, :key3),
                            (4, :key4), (5, :key5), (6, :key6),
                            (7, :key7), (8, :key8), (9, :key9),
                            (10, :key10),
                            (1, Symbol("Fixed Effects"))]
    end

    @test map(x->get_level(x, 1), c1.col_index) == [tuple(1, :test)]
    @test map(x->get_level(x, 1), c2.col_index) == [tuple(1, :test2)]
    @test_throws BoundsError map(x->get_level(x, 2), c1.col_index)
end

@testset "Index Schemas" begin
    # Standard Composite Tables
    Random.seed!(1234)
    x  = randn(10)
    y  = [Symbol(:key, i) for i=1:10]
    t1 = TableCol("test", y, x)
    t2 = TableCol("test2", y[2:9], x[2:9])
    t3 = TableCol("test3", y, x, randn(10) .|> abs .|> sqrt)
    t4 = TableCol("test" , Dict("Fixed Effects"=>"Yes"))
    t5 = TableCol("test2", Dict("Fixed Effects"=>"No"))
    t6 = TableCol("test3", Dict("Fixed Effects"=>"Yes"))
    c1 = append_table(t1, t4)
    c2 = append_table(t2, t5)
    c3 = append_table(t3, t6)

    @test generate_schema(c1.row_index, 1) == Any[(1, Symbol(""))=>10,
                                                  (2, Symbol(""))=>1]
    @test generate_schema(c1.row_index, 2) == begin
        Any[(1, :key1)=>1, (2, :key2)=>1, (3, :key3)=>1, (4, :key4)=>1,
            (5, :key5)=>1, (6, :key6)=>1, (7, :key7)=>1, (8, :key8)=>1,
            (9, :key9)=>1, (10, :key10)=>1, (1, Symbol("Fixed Effects"))=>1]
    end

    @test generate_schema(c3.row_index,1)==generate_schema(c1.row_index,1)
    @test generate_schema(c3.row_index,2)==generate_schema(c1.row_index,2)
    @test_throws BoundsError generate_schema(c3.row_index, 3)

    s = generate_schema(c2.row_index, 1)
    @test schema_lookup(s, 1) == 8
    @test schema_lookup(s, 8) == 8
    @test schema_lookup(s, 9) == 1

    s = generate_schema(c2.row_index, 2)
    @test schema_lookup(s, 1) == 1
    @test schema_lookup(s, 8) == 1
    @test schema_lookup(s, 9) == 1

end

@testset "Argument Checking" begin
    @test check_table_type(:ascii) == nothing
    @test check_table_type(:latex) == nothing
    @test_throws ArgumentError check_table_type(:html)
    @test_throws ArgumentError check_table_type("ascii")
    @test_throws ArgumentError check_table_type("latex")
    @test_throws ArgumentError check_table_type(5)
    @test_throws ArgumentError check_table_type(5.0)
end

@testset "table_type Defaults" begin
    @test default_sep(:ascii) == "|"
    @test default_sep(:latex) == "&"
    @test_throws ArgumentError default_sep("ascii")
end

# c1_a = TableCol("Column 1",
#                 OrderedDict("Row 1"=>1,
#                          "Row 2"=>2.3,
#                          "Row 3"=>(2.3, .3),
#                          "Row 4"=>(832.1, 20.0)))
# c1_b = TableCol("Column 1",
#                 OrderedDict("Stat 1" => 232))
#
# c1   = append_table(c1_a, c1_b)


text:
"""
    DAMMplot(Tₛ::Array{Float64, 1}, θ::Array{Float64, 1}, Rₛ::Array{Float64, 1}, r::Int64)

Plot scatter of data and fitted DAMM surface

# Example
```julia-repl
julia> df = DAMMfdata(100)
julia> r = 50
julia> fig = DAMMplot(df.Tₛ, df.θ, df.Rₛ, r)
```
"""
function DAMMplot(Tₛ::Array{Float64, 1}, θ::Array{Float64, 1}, Rₛ::Array{Float64, 1}, r::Int64)
  fontsize_theme = Theme(fontsize = 20, font = "JuliaMono")
  set_theme!(fontsize_theme)
  fig = Figure()
  ax3D = Axis3(fig[1, 1])
  out = DAMMmat(Tₛ, θ, Rₛ, r)
  ax3D.xlabel = to_latex("T_{soil} (°C)");
  ax3D.ylabel = to_latex("\\theta (m^3 m^{-3})");
  ax3D.zlabel = to_latex("R_{soil} (\\mumol m^{-2} s^{-1})");
  data3D = Vec3f.(Tₛ, θ, Rₛ)
  p3D = scatter!(ax3D, data3D, markersize = 2500, strokewidth = 3,
	color = Rₛ, colormap = Reverse(:Spectral))
  s3D = surface!(ax3D, out.x, out.y, out.DAMM_Matrix, colormap = Reverse(:Spectral),
	transparency = true, alpha = 0.1, shading = false)
  w3D = wireframe!(ax3D, out.x, out.y, out.DAMM_Matrix, overdraw = true,
	transparency = true, color = (:black, 0.1));
  Colorbar(fig[1, 2], colormap = Reverse(:Spectral),
	   label = to_latex("R_{soil} (\\mumol m^{-2} s^{-1})"));
  fig
  return fig
end


text:
#Copyright (c) 2015 Rex Computing and Isaac Yonemoto

#see LICENSE.txt

#this work was supported in part by DARPA Contract D15PC00135

#log-table.jl

#log table describing logarithmic fractions for all 64-bit fraction bits...  starting
#with 1.1b, 1.01b, 1.001b, 1.0001b etc.

const lt64 = [
0x5269e12f346e2c00,
0x2b803473f7ad1000,
0x1663f6fac9131600,
0x0b5d69bac77ec380,
0x05b9e5a170b48a80,
0x02dfca16dde10a20,
0x01709c46d7aac770,
0x00b87c1ff853ab28,
0x005c4994dd0fd150,
0x002e27ac5ef2af86,
0x0017148ec2a1bfc9,
0x000b8a7588fd29b1,
0x0005c5464ec5f4d7,
0x0002e2a60a005c95,
0x00017153bda8f822,
0x0000b8aa0cfedcb1,
0x00005c55120a0c45,
0x00002e2a8be7ae56,
0x0000171546ac814f,
0x00000b8aa3846b33,
0x000005c551cdc03d,
0x000002e2a8e9c2c7,
0x0000017154759a0d,
0x000000b8aa3afb31,
0x0000005c551d8923,
0x0000002e2a8ec774,
0x0000001715476472,
0x0000000b8aa3b267,
0x00000005c551d93f,
]

max_stars_repo_name: REX-Computing/unumjl

max_issues_repo_name: REX-Computing/unumjl

max_forks_repo_name: REX-Computing/unumjl


text:
# Uses SSDP to locate the IP for the OpenBCI WiFi board if it is
# functioning on the local network.

using SSDPClient

const ARDUINOMATCH = "SERVER: Arduino.+(http://[^/]+)/description"
println("The OpenBCI board responds at URI: ", ssdpquery(ARDUINOMATCH)[1])