text:
export AbstractSGD, SimpleSGD, AdaDelta, AdaGrad, AveragedSGD

abstract AbstractSGD
#should implement a update! method:
#update!(obj::AbstractSGD, weights::Vector{Float64}, gr::Vector{Float64})

#This function returns the vector of weights at which to evaluate the gradient
#It should be overridden by SGD implementations that need the gradient evaluated at
#values other than those stored in the Learner (AveragedSGD)
which_weights(obj::AbstractSGD, weights) = weights

type SimpleSGD <: AbstractSGD
    alpha1::Float64
    alpha2::Float64
    t::Int
    function SimpleSGD(alpha1::Float64, alpha2::Float64)
        alpha1 <= 0.0 && error("alpha1 should be positive")
        alpha2 < 0.0 && error("alpha2 should be non-negative")
        new(alpha1, alpha2, 0)
    end
end

function update!(obj::SimpleSGD, weights::Vector{Float64}, gr::Vector{Float64})
    obj.t += 1

    stepsize = - obj.alpha1 / (1.0 + obj.alpha1 * obj.alpha2 * obj.t)
    fma!(weights, gr, stepsize)
    weights
end

type AdaDelta <: AbstractSGD
    rho::Float64
    eps::Float64
    sqgr::Vector{Float64}
    squp::Vector{Float64}
    up::Vector{Float64}
    initialized::Bool
    function AdaDelta(rho::Float64, eps::Float64)
        (rho <= 0.0 || eps <= 0.0) && error("rho and epsilon should be positive")
        obj = new()
        obj.rho = rho
        obj.eps = eps
        obj.initialized = false
        obj
    end
end

Base.show(io::IO, obj::AdaDelta) = print(io, "AdaDelta(ρ=$(obj.rho), ε=$(obj.eps))")

function init!(obj::AdaDelta, weights)
    obj.initialized && error("already initialized")
    obj.sqgr = zeros(weights)
    obj.squp = zeros(weights)
    obj.up = similar(weights)
    obj.initialized = true
end

function update!(obj::AdaDelta, weights::Vector{Float64}, gr::Vector{Float64})
    obj.initialized || init!(obj, weights)
    @devec obj.sqgr[:] = obj.rho .* obj.sqgr + (1.0 - obj.rho) .* gr .* gr #line 4
    @devec obj.up[:] = - sqrt(obj.squp .+ obj.eps) ./ sqrt(obj.sqgr .+ obj.eps) .* gr #line 5
    @devec obj.squp[:] = obj.rho .* obj.squp + (1.0 - obj.rho) .* obj.up  .* obj.up #line 6
    add!(weights, obj.up) #line 7
end

type AdaGrad <: AbstractSGD
    eta::Float64
    sqgr::Vector{Float64}
    initialized::Bool
    function AdaGrad(eta::Float64)
        eta > 0.0 || error("eta should be positive")
        obj = new()
        obj.eta = eta
        obj.initialized = false
        obj
    end
end

Base.show(io::IO, obj::AdaGrad) = print(io, "AdaGrad(η=$(obj.eta))")

function init!(obj::AdaGrad, weights)
    obj.initialized && error("already initialized")
    obj.sqgr = fill!(similar(weights), 1.0e-8)
    obj.initialized = true
end

function update!(obj::AdaGrad, weights::Vector{Float64}, gr::Vector{Float64})
    obj.initialized || init!(obj, weights)
    @devec obj.sqgr[:] += gr .* gr
    @devec weights[:] -= obj.eta ./ sqrt(obj.sqgr) .* gr
    weights
end

type AveragedSGD <: AbstractSGD
    alpha1::Float64
    alpha2::Float64
    unaveraged_weights::Vector
    t0::Int
    t::Int
    initialized::Bool
    function AveragedSGD(alpha1::Float64, alpha2::Float64, t0::Int)
        alpha1 <= 0.0 && error("alpha1 should be positive")
        alpha2 < 0.0 && error("alpha2 should be non-negative")
        obj = new()
        obj.alpha1 = alpha1
        obj.alpha2 = alpha2
        obj.t0 = t0
        obj.t = 0
        obj.initialized = false
        obj
    end
end

Base.show(io::IO, obj::AveragedSGD) = print(io, "AveragedSGD(α1=$(obj.alpha1), α2=$(obj.alpha2), t0=$(obj.t0))")

function init!(obj::AveragedSGD, weights)
    obj.initialized && error("already initialized")
    obj.unaveraged_weights = zeros(weights)
    obj.initialized = true
end

which_weights(obj::AveragedSGD, weights) = obj.initialized? obj.unaveraged_weights: weights

function update!(obj::AveragedSGD, weights::Vector{Float64}, gr::Vector{Float64})
    obj.initialized || init!(obj, weights)
    obj.t += 1
    stepsize = - obj.alpha1 * (1.0 + obj.alpha1 * obj.alpha2 * obj.t)^-0.75
    fma!(obj.unaveraged_weights, gr, stepsize)
    mu = 1.0 / max(1.0, obj.t - obj.t0)
    @devec weights[:] += mu .* (obj.unaveraged_weights .- weights)
    weights
end


text:
export AbstractHook,
    ComposedHook,
    EmptyHook,
    StepsPerEpisode,
    RewardsPerEpisode,
    TotalRewardPerEpisode,
    TotalBatchRewardPerEpisode,
    BatchStepsPerEpisode,
    TimePerStep,
    DoEveryNEpisode,
    DoEveryNStep,
    UploadTrajectoryEveryNStep,
    MultiAgentHook

using UnicodePlots:lineplot

"""
A hook is called at different stage duiring a [`run`](@ref) to allow users to inject customized runtime logic.
By default, a `AbstractHook` will do nothing. One can override the behavior by implementing the following methods:

- `(hook::YourHook)(::PreActStage, agent, env, action)`, note that there's an extra argument of `action`.
- `(hook::YourHook)(::PostActStage, agent, env)`
- `(hook::YourHook)(::PreEpisodeStage, agent, env)`
- `(hook::YourHook)(::PostEpisodeStage, agent, env)`
"""
abstract type AbstractHook end

(hook::AbstractHook)(args...) = nothing

# https://github.com/JuliaLang/julia/issues/14919
# function (f::Function)(stage::T, args...;kw...) where T<: AbstractStage end

#####
# ComposedHook
#####

"""
    ComposedHook(hooks::AbstractHook...)

Compose different hooks into a single hook.
"""
struct ComposedHook{T<:Tuple} <: AbstractHook
    hooks::T
    ComposedHook(hooks...) = new{typeof(hooks)}(hooks)
end

function (hook::ComposedHook)(stage::AbstractStage, args...; kw...)
    for h in hook.hooks
        h(stage, args...; kw...)
    end
end

Base.getindex(hook::ComposedHook, inds...) = getindex(hook.hooks, inds...)

#####
# EmptyHook
#####

"""
Do nothing
"""
struct EmptyHook <: AbstractHook end

const EMPTY_HOOK = EmptyHook()

#####
# display
#####

Base.display(::AbstractStage, agent, env, args...; kwargs...) = display(env)

#####
# StepsPerEpisode
#####

"""
    StepsPerEpisode(; steps = Int[], count = 0)

Store steps of each episode in the field of `steps`.
"""
Base.@kwdef mutable struct StepsPerEpisode <: AbstractHook
    steps::Vector{Int} = Int[]
    count::Int = 0
end

Base.getindex(h::StepsPerEpisode) = h.steps

(hook::StepsPerEpisode)(::PostActStage, args...) = hook.count += 1

function (hook::StepsPerEpisode)(::Union{PostEpisodeStage,PostExperimentStage}, agent, env)
    push!(hook.steps, hook.count)
    hook.count = 0
end

#####
# RewardsPerEpisode
#####

"""
    RewardsPerEpisode(; rewards = Vector{Vector{Float64}}())

Store each reward of each step in every episode in the field of `rewards`.
"""
Base.@kwdef mutable struct RewardsPerEpisode <: AbstractHook
    rewards::Vector{Vector{Float64}} = Vector{Vector{Float64}}()
end

Base.getindex(h::RewardsPerEpisode) = h.rewards

function (hook::RewardsPerEpisode)(::PreEpisodeStage, agent, env)
    push!(hook.rewards, [])
end

function (hook::RewardsPerEpisode)(::PostActStage, agent, env)
    push!(hook.rewards[end], reward(env))
end

function (hook::RewardsPerEpisode)(::PostActStage, agent::NamedPolicy, env)
    push!(hook.rewards[end], reward(env, nameof(agent)))
end

#####
# TotalRewardPerEpisode
#####

"""
    TotalRewardPerEpisode(; rewards = Float64[], reward = 0.0, is_display_on_exit = true)

Store the total reward of each episode in the field of `rewards`. If
`is_display_on_exit` is set to `true`, a unicode plot will be shown at the [`PostExperimentStage`](@ref).
"""
Base.@kwdef mutable struct TotalRewardPerEpisode <: AbstractHook
    rewards::Vector{Float64} = Float64[]
    reward::Float64 = 0.0
    is_display_on_exit::Bool = true
end

Base.getindex(h::TotalRewardPerEpisode) = h.rewards

function (hook::TotalRewardPerEpisode)(::PostActStage, agent, env)
    hook.reward += reward(env)
end

function (hook::TotalRewardPerEpisode)(::PostActStage, agent::NamedPolicy, env)
    hook.reward += reward(env, nameof(agent))
end

function (hook::TotalRewardPerEpisode)(::PostEpisodeStage, agent, env)
    push!(hook.rewards, hook.reward)
    hook.reward = 0
end

function (hook::TotalRewardPerEpisode)(::PostExperimentStage, agent, env)
    println(lineplot(hook.rewards, title="Total reward per episode", xlabel="Episode", ylabel="Score"))
end

#####
# TotalBatchRewardPerEpisode
#####
struct TotalBatchRewardPerEpisode <: AbstractHook
    rewards::Vector{Vector{Float64}}
    reward::Vector{Float64}
end

Base.getindex(h::TotalBatchRewardPerEpisode) = h.rewards

"""
    TotalBatchRewardPerEpisode(batch_size::Int)

Similar to [`TotalRewardPerEpisode`](@ref), but is specific to environments
which return a `Vector` of rewards (a typical case with `MultiThreadEnv`).
"""
function TotalBatchRewardPerEpisode(batch_size::Int)
    TotalBatchRewardPerEpisode([Float64[] for _ in 1:batch_size], zeros(batch_size))
end

function (hook::TotalBatchRewardPerEpisode)(::PostActStage, agent, env)
    R = agent isa NamedPolicy ? reward(env, nameof(agent)) : reward(env)
    for (i, (t, r)) in enumerate(zip(is_terminated(env), R))
        hook.reward[i] += r
        if t
            push!(hook.rewards[i], hook.reward[i])
            hook.reward[i] = 0.0
        end
    end
end

#####
# BatchStepsPerEpisode
#####

struct BatchStepsPerEpisode <: AbstractHook
    steps::Vector{Vector{Int}}
    step::Vector{Int}
end

Base.getindex(h::BatchStepsPerEpisode) = h.steps

"""
    BatchStepsPerEpisode(batch_size::Int; tag = "TRAINING")

Similar to [`StepsPerEpisode`](@ref), but is specific to environments
which return a `Vector` of rewards (a typical case with `MultiThreadEnv`).
"""
function BatchStepsPerEpisode(batch_size::Int)
    BatchStepsPerEpisode([Int[] for _ in 1:batch_size], zeros(Int, batch_size))
end

function (hook::BatchStepsPerEpisode)(::PostActStage, agent, env)
    for (i, t) in enumerate(is_terminated(env))
        hook.step[i] += 1
        if t
            push!(hook.steps[i], hook.step[i])
            hook.step[i] = 0
        end
    end
end

#####
# TimePerStep
#####

"""
    TimePerStep(;max_steps=100)
    TimePerStep(times::CircularArrayBuffer{Float64}, t::UInt64)

Store time cost of the latest `max_steps` in the `times` field.
"""
mutable struct TimePerStep <: AbstractHook
    times::CircularArrayBuffer{Float64,1}
    t::UInt64
end

Base.getindex(h::TimePerStep) = h.times

TimePerStep(; max_steps = 100) =
    TimePerStep(CircularArrayBuffer{Float64}(max_steps), time_ns())

function (hook::TimePerStep)(::PostActStage, agent, env)
    push!(hook.times, (time_ns() - hook.t) / 1e9)
    hook.t = time_ns()
end

"""
    DoEveryNStep(f; n=1, t=0)

Execute `f(t, agent, env)` every `n` step.
`t` is a counter of steps.
"""
mutable struct DoEveryNStep{F} <: AbstractHook
    f::F
    n::Int
    t::Int
end

DoEveryNStep(f; n = 1, t = 0) = DoEveryNStep(f, n, t)

function (hook::DoEveryNStep)(::PostActStage, agent, env)
    hook.t += 1
    if hook.t % hook.n == 0
        hook.f(hook.t, agent, env)
    end
end

"""
    DoEveryNEpisode(f; n=1, t=0)

Execute `f(agent, env)` every `n` episode.
`t` is a counter of steps.
"""
Base.@kwdef mutable struct DoEveryNEpisode{S<:Union{PreEpisodeStage,PostEpisodeStage},F} <:
                           AbstractHook
    f::F
    n::Int = 1
    t::Int = 0
end

DoEveryNEpisode(f::F, n = 1, t = 0; stage::S = POST_EPISODE_STAGE) where {S,F} =
    DoEveryNEpisode{S,F}(f, n, t)

function (hook::DoEveryNEpisode{S})(::S, agent, env) where {S}
    hook.t += 1
    if hook.t % hook.n == 0
        hook.f(hook.t, agent, env)
    end
end

"""
    UploadTrajectoryEveryNStep(;mailbox, n, sealer=deepcopy)
"""
Base.@kwdef mutable struct UploadTrajectoryEveryNStep{M,S} <: AbstractHook
    mailbox::M
    n::Int
    t::Int = -1
    sealer::S = deepcopy
end

function (hook::UploadTrajectoryEveryNStep)(::PostActStage, agent::Agent, env)
    hook.t += 1
    if hook.t > 0 && hook.t % hook.n == 0
        put!(hook.mailbox, hook.sealer(agent.trajectory))
    end
end

"""
    MultiAgentHook(player=>hook...)
"""
struct MultiAgentHook <: AbstractHook
    hooks::Dict{Any,Any}
end

MultiAgentHook(player_hook_pair::Pair...) = MultiAgentHook(Dict(player_hook_pair...))

Base.getindex(h::MultiAgentHook, p) = getindex(h.hooks, p)

function (hook::MultiAgentHook)(
    s::AbstractStage,
    m::MultiAgentManager,
    env::AbstractEnv,
    args...,
)
    for (p, h) in zip(values(m.agents), values(hook.hooks))
        h(s, p, env, args...)
    end
end

ext: jl


text:
getalg(alg0::DEAlgorithm) = alg0
getalg(alg0::MethodOfSteps) = alg0.alg

"""
    has_constant_lags(integrator::DDEIntegrator)

Return if the DDE problem of the `integrator` contains constant delays.
"""
has_constant_lags(integrator::SDDEIntegrator) = has_constant_lags(integrator.sol.prob)

"""
    has_dependent_lags(integrator::DDEIntegrator)

Return if the DDE problem of the `integrator` contains dependent delays.
"""
has_dependent_lags(integrator::SDDEIntegrator) = has_dependent_lags(integrator.sol.prob)

"""
    has_constant_lags(prob::DDEProblem)

Return if the DDE problem `prob` contains constant delays.
"""
has_constant_lags(prob::SDDEProblem) =
  prob.constant_lags !== nothing && !isempty(prob.constant_lags)

"""
    has_dependent_lags(prob::DDEProblem)

Return if the DDE problem `prob` contains dependent delays.
"""
has_dependent_lags(prob::SDDEProblem) =
  prob.dependent_lags !== nothing && !isempty(prob.dependent_lags)

function u_uprev(u0; alias_u0 = false)
  if  typeof(u0) <: Tuple
    u = ArrayPartition(prob.u0,Val{true})
  else
    if alias_u0
      u = u0
    else
      u = recursivecopy(u0)
    end
  end
  uprev = recursivecopy(u)
  u, uprev
end

"""
    solution_arrays(u, tspan, rate_prototype; kwargs...)

Return arrays of saved time points, states, and rates, initialized with the solution at the
first time point if `save_start = true` (the default).
"""
function solution_arrays(u::uType, tspan, rate_prototype;
  timeseries_init = typeof(u)[],
  ts_init = eltype(tspan)[],
  save_idxs = nothing,
  save_start = true) where {uType}
# determine types of time and state
  # uType = typeof(u)
  tType = eltype(tspan)

  # initialize vector of saved time points
  ts = convert(Vector{tType}, ts_init)
  alg_choice = Int[]

  # initialize vector of saved states
  if save_idxs === nothing
    timeseries = convert(Vector{uType}, timeseries_init)
  else
    u_initial = u[save_idxs]
    timeseries = convert(Vector{typeof(u_initial)}, timeseries_init)
  end

  # save solution at initial time point
  if save_start
    saveiter=1
    copyat_or_push!(ts, 1, first(tspan))
    if save_idxs === nothing
      copyat_or_push!(timeseries, 1, u)
    else
      copyat_or_push!(timeseries, 1, u_initial, Val{false})
    end
  else
    saveiter=0;
  end

  ts, timeseries, saveiter
end
"""
    sizehint!(sol::DESolution, n)

Suggest that solution `sol` reserves capacity for at least `n` elements.
"""
function Base.sizehint!(sol::RODESolution, n)
  sizehint!(sol.u, n)
  sizehint!(sol.t, n)

  nothing
end

"""
    sizehint!(sol::DESolution, alg, tspan, tstops, saveat; kwargs...)

Suggest that solution `sol` reserves capacity for a number of elements that
depends on the parameter settings of the numerical solver.
"""
function Base.sizehint!(sol::RODESolution, alg, tspan, tstops, saveat;
                        save_everystep = isempty(saveat),
                        adaptive = StochasticDiffEq.isadaptive(getalg(alg)),
                        internalnorm = DiffEqBase.ODE_DEFAULT_NORM,
                        dt = zero(eltype(tspan)))
  # obtain integration time
  t0 = first(tspan)
  integrationtime = last(tspan) - t0

  if !adaptive && save_everystep && tspan[2]-tspan[1] != Inf
    # determine number of steps if known a priori
    if iszero(dt)
      steps = length(tstops)
    else
      steps = ceil(Int, internalnorm(integrationtime / dt, t0))
    end

    sizehint!(sol,steps+1)
  elseif save_everystep
    sizehint!(sol,50)
  elseif !isempty(saveat)
    sizehint!(sol,length(saveat)+1)
  else
    sizehint!(sol,2)
  end

  nothing
end

function build_history_function(prob, alg, reltol, rate_prototype, noise_rate_prototype, W, _seed, dense;
  dt = zero(eltype(prob.tspan)),
  adaptive = StochasticDiffEq.isadaptive(getalg(alg)),
  calck = false,
  internalnorm = DiffEqBase.ODE_DEFAULT_NORM)
  @unpack f, g, h, u0, tspan, p = prob

  t0 = first(tspan)
  tType = eltype(tspan)
  tTypeNoUnits = typeof(one(tType))
  tdir = sign(last(tspan) - t0)

  uEltypeNoUnits = recursive_unitless_eltype(u0)
  uBottomEltypeNoUnits = recursive_unitless_bottom_eltype(u0)

  # bootstrap an SDE integrator
  # - whose solution captures the dense history of the simulation
  # - that is used for extrapolation of the history for time points past the
  #   already fixed history
  # - that is used for interpolation of the history for time points in the
  #   current integration step (so the interpolation is fixed while updating the stages)
  # we wrap the user-provided history function such that function calls during the setup
  # of the integrator do not fail
  # sde_f = SDEFunctionWrapper(f, prob.h)
  # sde_g = SDEDiffusionTermWrapper(g,prob.h)

  sde_f,sde_g = wrap_functions_and_history(f, g, h)
  sde_prob = SDEProblem{isinplace(prob)}(sde_f, sde_g, u0, tspan, p)

  # get states of ODE integrator (do not alias uprev)
  sde_u, sde_uprev = u_uprev(u0; alias_u0 = false)

  # # initialize output arrays
  sde_ts, sde_timeseries, sde_saveiter = solution_arrays(sde_u, tspan, rate_prototype,
                                      save_idxs = nothing,
                                      save_start = true)

  # # obtain cache (we alias uprev2 and uprev)
  sde_cache = StochasticDiffEq.alg_cache(getalg(alg),prob, sde_u, W.dW, W.dZ, p, rate_prototype,noise_rate_prototype, uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits, sde_uprev, f, t0, dt, Val{isinplace(prob)})

  # build dense interpolation of history
  id = StochasticDiffEq.LinearInterpolationData(sde_timeseries,sde_ts)
  if typeof(getalg(alg)) <: StochasticDiffEq.StochasticDiffEqCompositeAlgorithm
    alg_choice = Int[]
    sde_sol = DiffEqBase.build_solution(prob,alg,sde_ts,sde_timeseries,W=W,
                                    destats = DiffEqBase.DEStats(0),
                                    calculate_error = false, alg_choice=alg_choice,
                                    interp = id, dense = dense, seed = _seed)
  else
    sde_sol = DiffEqBase.build_solution(prob,alg,sde_ts,sde_timeseries,W=W,
                                    destats = DiffEqBase.DEStats(0),
                                    calculate_error = false,
                                    interp = id, dense = dense, seed = _seed)
  end

  # # reserve capacity
  sizehint!(sde_sol, getalg(alg), tspan, (), ();
    save_everystep = true, adaptive = adaptive, internalnorm = internalnorm, dt = tType(dt))

  # # create simple integrator
  sde_integrator = HistorySDEIntegrator{typeof(getalg(alg)),isinplace(prob),typeof(prob.u0),
            tType, typeof(sde_sol),typeof(sde_cache)}(
              sde_sol, sde_u, t0, zero(tType), sde_uprev,
              t0, getalg(alg), zero(tType), tdir, 1, sde_cache)

  # # combine the user-provided history function and the ODE integrator with dense solution
  # # to a joint dense history of the DDE
  # # we use this history information to create a problem function of the DDE with all
  # # available history information that is of the form f(du,u,p,t) or f(u,p,t) such that
  # # SDE algorithms can be applied
  HistoryFunction(prob.h, sde_integrator)
end

"""
    initialize_solution!(integrator::SDDEIntegrator)

Initialize the solution of an integrator by adjusting the cache for composite algorithms.
"""
function initialize_solution!(integrator::SDDEIntegrator)
  if iscomposite(getalg(integrator.alg))
    copyat_or_push!(integrator.integrator.sol.alg_choice, 1, integrator.cache.current)
    if integrator.opts.save_start
      copyat_or_push!(integrator.sol.alg_choice, 1, integrator.cache.current)
    end
  end

  nothing
end

DiffEqBase.nlsolve_f(integrator::SDDEIntegrator) =
  DiffEqBase.nlsolve_f(integrator.f, unwrap_alg(integrator, true))

function unwrap_alg(integrator::SDDEIntegrator, is_stiff)
  alg = integrator.alg
  iscomp = typeof(alg) <: StochasticDiffEq.StochasticCompositeAlgorithm
  if !iscomp
    return alg
  elseif typeof(alg.choice_function) <: DiffEqBase.AutoSwitch
    num = is_stiff ? 2 : 1
    return alg.algs[num]
  else
    return alg.algs[integrator.cache.current]
  end
end

function DiffEqBase.unwrap_cache(integrator::SDDEIntegrator, is_stiff)
  alg   = getalg(integrator.alg)
  cache = integrator.cache
  iscomp = alg isa StochasticDiffEq.StochasticCompositeAlgorithm
  if !iscomp
    return cache
  elseif alg.choice_function isa DiffEqBase.AutoSwitch
    num = is_stiff ? 2 : 1
    return cache.caches[num]
  else
    return cache.caches[integrator.cache.current]
  end
end


text:
using BinaryProvider
using CMakeWrapper: cmake_executable
using Libdl: dlext

function compile(libname, tarball_url, hash; prefix=BinaryProvider.global_prefix, verbose=false)
    # download to tarball_path
    tarball_path = joinpath(prefix, "downloads", "libmakef.tar.gz")
    download_verify(tarball_url, hash, tarball_path; force=true, verbose=verbose)

    # unpack into source_path
    tarball_dir = joinpath(prefix, "downloads", dirname(first(list_tarball_files(tarball_path))))
    source_path = joinpath(prefix, "downloads", "src")

    verbose && @info("Unpacking $tarball_path into $source_path")
    rm(tarball_dir, force=true, recursive=true)
    rm(source_path, force=true, recursive=true)
    unpack(tarball_path, dirname(tarball_dir); verbose=verbose)
    mv(tarball_dir, source_path)

    build_dir = joinpath(source_path, "build")
    src_dir   = joinpath(build_dir, "lib")

    mkdir(build_dir)
    verbose && @info("Compiling in $build_dir...")
    cd(build_dir) do
        run(`$cmake_executable ..`)
        run(`$cmake_executable --build ./lib`)
        mkpath(libdir(prefix))

        src_file = joinpath(src_dir, libname*"."*dlext)
        dst_file = joinpath(libdir(prefix), libname*"."*dlext)
        cp(src_file, dst_file, force=true, follow_symlinks=true)
    end
end


text:
function test_combination_combinations()
    @test collect(combinations(Int[], 0, Val{:unique})) == Array{Int,1}[[]]
    @test length(combinations(Int[], 0, Val{:unique})) == 1
    @test eltype(combinations(Int[], 0, Val{:unique})) == Array{Int,1}

    @test collect(combinations(Int[], 1, Val{:unique})) == Array{Int,1}[]
    @test length(combinations(Int[], 1, Val{:unique})) == 0
    @test eltype(combinations(Int[], 1, Val{:unique})) == Array{Int,1}

    @test collect(combinations(["a", "b", "c", "d"], 2, Val{:unique})) == Array[
        ["a", "b"],
        ["a", "c"],
        ["a", "d"],
        ["b", "c"],
        ["b", "d"],
        ["c", "d"],
    ]
    @test length(combinations(["a", "b", "c", "d"], 2, Val{:unique})) == 6
    @test eltype(combinations(["a", "b", "c", "d"], 2, Val{:unique})) == Array{String,1}

    @test collect(combinations(["a", "b"], 3, Val{:unique})) == Array{String,1}[]
    @test length(combinations(["a", "b"], 3, Val{:unique})) == 0
end

function test_combination_combinations_with_repetition()
    @test collect(combinations(Int[], 0, Val{:repeated})) == Array{Int,1}[[]]
    @test length(combinations(Int[], 0, Val{:repeated})) == 1
    @test eltype(combinations(Int[], 0, Val{:repeated})) == Array{Int,1}

    @test collect(combinations(Int[], 1, Val{:repeated})) == Array{Int,1}[]
    @test length(combinations(Int[], 1, Val{:repeated})) == 0
    @test eltype(combinations(Int[], 1, Val{:repeated})) == Array{Int,1}

    @test collect(combinations(["a", "b", "c", "d"], 2, Val{:repeated})) == Array[
        ["a", "a"],
        ["a", "b"],
        ["a", "c"],
        ["a", "d"],
        ["b", "b"],
        ["b", "c"],
        ["b", "d"],
        ["c", "c"],
        ["c", "d"],
        ["d", "d"],
    ]
    @test length(combinations(["a", "b", "c", "d"], 2, Val{:repeated})) == 10
    @test eltype(combinations(["a", "b", "c", "d"], 2, Val{:repeated})) == Array{String,1}

    @test collect(combinations(["a", "b"], 3, Val{:repeated})) == Array[
        ["a", "a", "a"],
        ["a", "a", "b"],
        ["a", "b", "b"],
        ["b", "b", "b"],
    ]
    @test length(combinations(["a", "b"], 3, Val{:repeated})) == 4
end

function test_combination_combinations_unknown_mode()
    @test_throws MethodError combinations(["a", "b", "c", "d"], 2, Val{:unknown})
end

function test_combination_all()
    print(rpad("Math.Combinatorics.Combination...", 50, ' '))

    test_combination_combinations()
    test_combination_combinations_with_repetition()
    test_combination_combinations_unknown_mode()

    println("PASS")
end


text:
using Parameters
import Roots
import Calculus

export ZeroTemperature,
      ZeroTemperatureGrand,
      MerminFreeEnergy,
      GrandPotential,
      FermiDiracSmearing,
      ZeroT,
      ZeroTGrand

# a smearing function can be
#  [1] 0T or finite T
#  [2] fixed eF or variable eF
# TODO: ZeroTemperature models, MerminFreeEnergy

# ChemicalPotentials define:
# * occupancy   n(ϵ)
# * energy   f(ϵ)
# * grad     f'(ϵ)
#     NB: grad is defined via ForwardDiff, but for an efficient implemention
#         this can of course be overloaded

# AD derivative of potential and occupancy
grad(f::FiniteTPotential, epsn::Real) =
      ForwardDiff.derivative(s -> energy(f, s), epsn)
occupancy_d(f::FiniteTPotential, epsn::Real, args...) =
      ForwardDiff.derivative(s -> occupancy(f, s, args...), epsn)

# vectorized versions
energy(f::ChemicalPotential, epsn::AbstractVector, args...) =
      [energy(f, es, args...) for es in epsn]
occupancy(f::ChemicalPotential, epsn::AbstractVector, args...) =
      [occupancy(f, es, args...) for es in epsn]
grad(f::ChemicalPotential, epsn::AbstractVector, args...) =
      [grad(f, s, args...) for s in epsn]
occupancy_d(f::FiniteTPotential, epsn::AbstractVector, args...) =
      [occupancy_d(f, s, args...) for s in epsn]

# default implementations for extracting beta and eF
beta(f::FiniteTPotential) = f.beta
get_eF(f::FiniteTPotential) = f.eF

function set_eF!(fd::ChemicalPotential, eF)
   fd.eF = eF
end

function set_beta!(fd::FiniteTPotential, beta)
   fd.beta = beta
end

# total occupancy calculation

occupancy(at::AbstractAtoms, f::ChemicalPotential, tbm, μ = get_eF(f)) =
   sum( w * occupancy(f, epsn, μ) for (w, _k, epsn, _ψ) in BZiter(tbm, at) )

occupancy_d(at::AbstractAtoms, f::ChemicalPotential, tbm, μ = get_eF(f)) =
   sum( w * occupancy_d(f, epsn, μ) for (w, _k, epsn, _ψ) in BZiter(tbm, at) )

# ================= Fermi-Dirac Distribution  ========================

# symbolic differentiation of the Fermi-dirac distribution
_fd0_ = :(1.0 / ( 1.0 + exp((epsn - eF) * beta) ))
_fd1_ = Calculus.differentiate(_fd0_, :epsn)
_fd2_ = Calculus.differentiate(_fd1_, :epsn)
_fd3_ = Calculus.differentiate(_fd2_, :epsn)

eval( :( fermidirac(epsn::Number, eF, beta) = $_fd0_ ) )
eval( :( fermidirac1(epsn, eF, beta) = $_fd1_ ) )
eval( :( fermidirac2(epsn, eF, beta) = $_fd2_ ) )
eval( :( fermidirac3(epsn, eF, beta) = $_fd3_ ) )

fermidirac(epsn::AbstractVector, eF, beta) = [fermidirac(e, eF, beta) for e in epsn]

# ================= Grand Potential ============================

mutable struct GrandPotential <: FiniteTPotential
   beta::Float64
   eF::Float64
end

@pot GrandPotential

"""
Finite temperature grand-potential, e ↦ 2/β log(1 - f(e)),
where f is the Fermi-dirac function.
"""
GrandPotential

GrandPotential(; β=30.0, beta=β, eF = 0.0) = GrandPotential(beta, eF)

occupancy(f::GrandPotential, epsn::Number, eF) = fermidirac(epsn, eF, f.beta)
occupancy(f::GrandPotential, epsn::Number) = occupancy(f, epsn, f.eF)

# _logfd(z) = real(z) < 0 ? z - log(1+exp(z)) : log( exp(z) / (1+exp(z)) )
_logfd(z) = real(z) < 0 ? z - log1p(exp(z)) : -log1p(exp(-z));
# _logfd(z) = log( exp(z) / (1+exp(z)) )
# _logfd(z) = z - log(1+exp(z))
grand(epsn, eF, beta) = 2.0/beta * _logfd(beta*(epsn-eF))

energy(f::GrandPotential, epsn::Number) = grand(epsn, f.eF, f.beta)

update!(at::AbstractAtoms, f::GrandPotential, tbm::TBModel) = nothing

set_Nel!(f::GrandPotential, tbm, at, Nel) = set_eF!(f, eF_solver(at, f, tbm, Nel))

# get_Ne(f::GrandPotential) = occupancy(at::AbstractAtoms, f::ChemicalPotential, tbm, μ = get_eF(f))

fixed_eF(::GrandPotential) = true

# ================= The Old Smearing Function ============================
# should still implement it, then deprecate and remove once
# we have fixed the Atoms.jl - SKTB.jl equivalence

mutable struct FermiDiracSmearing  <: FiniteTPotential
    beta::Float64
    eF::Float64
    Ne::Float64
    fixed_eF::Bool
end

@pot FermiDiracSmearing

get_Ne(f::FermiDiracSmearing) = f.Ne
fixed_eF(f::FermiDiracSmearing) = f.fixed_eF

"""`FermiDiracSmearing`:

f(e) = ( 1 + exp( beta (e - eF) ) )^{-1}

Constructor: `FermiDiracSmearing(beta; eF=0.0, Ne = 0.0, fixed_eF = true)`
"""
FermiDiracSmearing

FermiDiracSmearing(beta; eF=0.0, Ne = 0.0, fixed_eF = true) =
      FermiDiracSmearing(beta, eF, Ne, fixed_eF)

occupancy(fd::FermiDiracSmearing, epsn::Number) = fermidirac(epsn, fd.eF, fd.beta)
occupancy(fd::FermiDiracSmearing, epsn::Number, eF) = fermidirac(epsn, eF, fd.beta)
energy(fd::FermiDiracSmearing, epsn::Number) = fermidirac(epsn, fd.eF, fd.beta) * epsn
# avoid NaN
function grad(f::FermiDiracSmearing, epsn::Real)
   exp_num = (epsn - f.eF) * f.beta
   exp_item = exp(exp_num)
   if exp_num > 0
      exp_item_inv = exp(-exp_num)
      return (exp_item_inv + 1 - epsn * f.beta) / (exp_item + 2 + exp_item_inv)
   else
      return (1 + (1 - epsn * f.beta) * exp_item) / (exp_item ^ 2 + 2 * exp_item + 1)
   end
end

function update!(at::AbstractAtoms, f::FermiDiracSmearing, tbm::TBModel)
   if !f.fixed_eF
      f.eF = eF_solver(at, f, tbm)
   end
   return nothing
end

"""
`eF_solver(at, f, tbm, Ne = get_Ne(f))`

* `at`: `AbstractAtoms`
* `f`: potential (e.g. `FermiDiracSmearing`, `GrandPotential`, etc)
* `tbm`: `SKTBModel`
* `Ne` : number of electrons in the system

Given an electron number `Ne` compute the fermi-level (chemical potential)
such that `∑_s f_s = Ne` (where `f_s` is the occupation number)
"""
function eF_solver(at, f, tbm, Ne = get_Ne(f))
   update!(at, tbm)
   # guess-timate an initial eF (this assumes Fermi-Dirac like behaviour)
   nf = ceil(Int, Ne)
   μ = 0.0
   for (w, k) in tbm.bzquad
      epsn_k = get_k_array(at, :epsn, k)
      μ += w * (epsn_k[nf-1] + epsn_k[nf]) / 2
   end
   # call the Roots package
   μ = Roots.fzero( _μ -> Ne - occupancy(at, f, tbm, _μ), μ )
   @assert abs(Ne - occupancy(at, f, tbm, μ)) < 1e-6

   return μ
end

"""
`set_δNel!(tbm::TBModel, at::AbstractAtoms, δNel = 0.0)`

Set the particle number (number of electrons) to `(ndofs(at) + δNel) / 2`.
"""
set_δNel!(tbm::TBModel, at::AbstractAtoms, δNel = 0.0) =
   set_Nel!(tbm, at, (ndofs(tbm.H, at)+δNel)/2)

set_Nel!(tbm::TBModel, at::AbstractAtoms, Nel) =
   set_Nel!(tbm.potential, tbm, at, Nel)

function set_Nel!(f::FermiDiracSmearing, tbm, at, Nel)
   f.Ne = Nel
   f.eF = eF_solver(at, f, tbm, Nel)
   return nothing
end

# ================= Zero-Temperature Models  ============================

function fermilevel(tbm::TBModel, at::AbstractAtoms, Nel)
   e = spectrum(tbm, at) |> sort
   # Nel = length(e) ÷ 2
   nf = floor(Int, Nel * length(tbm.bzquad) / 2)
   return 0.5 * (e[nf] + e[nf + 1])
end

abstract type ZeroTPotential <: ChemicalPotential end

"""
`ZeroT`: 0T canonical model (Nel fixed)
"""
@with_kw mutable struct ZeroT <: ZeroTPotential
   Nel::Float64 = 0.0
   eF::Float64 = 0.0
end

"""
`ZeroTGrand`: 0T Grand-canonical model (eF fixed)
"""
@with_kw mutable struct ZeroTGrand <: ZeroTPotential
   Nel::Float64 = 0.0
   eF::Float64 = 0.0
end

beta(f::ZeroTPotential) = 0.0

get_eF(f::ZeroTPotential) = f.eF

get_Nel(f::ZeroTPotential) = f.Nel

occupancy(f::ZeroTPotential, epsn::Number) = occupancy(f, epsn, f.eF)
occupancy(f::ZeroTPotential, epsn::Number, μ) = epsn < μ ? 1.0 : 0.0

energy(f::ZeroTPotential, epsn::Number) = occupancy(f, epsn) * epsn

grad(f::ZeroTPotential, epsn::Number) = occupancy(f, epsn)

function set_Nel!(f::ZeroTPotential, tbm, at, Nel)
   f.Nel = Nel
   f.eF = fermilevel(tbm, at, Nel)
   return nothing
end

function update!(at::AbstractAtoms, f::ZeroT, tbm::TBModel)
   f.eF = fermilevel(tbm, at, f.Nel)
   return nothing
end

function update!(at::AbstractAtoms, f::ZeroTGrand, tbm::TBModel)
   f.Nel = occupancy(at, f, tbm)
   return nothing
end

# # ================= Canonical Ensemble (Mermin) ============================
#
# @pot mutable struct MerminFreeEnergy <: ChemicalPotential
#    Nel::Float64
#    beta::Float64
#    eF::Float64
# end
#
# """
# Mermin Free energy is given by
# `e ↦  2 e f(e - μ) + 2 β⁻¹ entropy(f(e-μ))`
#
# where μ is determined from  `2 ∑_s f(e_s - μ) = Nel.`
# """
# MerminFreeEnergy
#
# _en0_ = :( f * log(f) - (1-f) * log(1-f) )
# _en1_ = Calculus.differentiate(_en0_, :f)
# _en2_ = Calculus.differentiate(_en1_, :f)
# _en3_ = Calculus.differentiate(_en2_, :f)
#
# eval( :( entropy(f) = $_en0_ ) )
# eval( :( entropy1(f) = $_en1_ ) )
# eval( :( entropy2(f) = $_en2_ ) )
# eval( :( entropy3(f) = $_en3_ ) )
#
#
# function update!(at::AbstractAtoms, f::MerminFreeEnergy)
#    error("`update!` for `MerminFreeEnergy` has not been implemented yet.")
#    # need to solve the nonlinear system that ensures
#    # ∑_k ∑_s f(ϵ_sk) = Nel/2    (or whatever)
# end


text:
const PlusMinus = "±"
const Exactly   = "⌁"
const Inexactly = "~"
const OpenedOpened = "⍿"
const ClosedOpened = "⫯"
const OpenedClosed = "⫰"
const ClosedClosed = "∣"

const PosInfinity = "∞"
const NegInfinity = "-∞"
const Empty = "∅"
const Full  = "▣"

#                        ClCl          ClOp            OpCl      OpOp
const Postfixes = ( ClosedClosed, ClosedOpened, OpenedClosed, OpenedOpened );
@inline postfixes(loIsOpen::Bool, hiIsOpen::Bool) =
      Postfixes[ one(Int8)+(reinterpret(Int8,loIsOpen)<<1)+hiIsOpen ]

# select characters corresponding to boundries as specified (lo,hi)
#   true for Open, false for Closed
#                     ClCl      ClOp        OpCl      OpOp
# const Delimiters = [ ("▪","▪"), ("▪","▫"), ("▫","▪"), ("▫","▫") ];
const Delimiters = [ ( "⟨" , "⟩" ) , ( "⟨" , "⟫" ) ,  ( "⟪" , "⟩" ), ( "⟪" , "⟫" )  ];
@inline delimiters(loIsOpen::Bool, hiIsOpen::Bool) =
      Delimiters[ one(Int8)+(reinterpret(Int8,loIsOpen)<<1)+hiIsOpen ]

function show{S<:Sculpt, Q<:Qualia, C<:Clay}(io::IO, x::Flex{S,Q,C})
    tiesym = (Q==INEXACT) ? Inexactly : Exactly
    loDelim,hiDelim = delimiters(boundries(x)...)
    s = (x.lo != x.hi) ? string(loDelim, x.lo, tiesym, x.hi, hiDelim) : string(loDelim, x.lo, tiesym, hiDelim)
    print(io, s)
end

function showcompact{S<:Sculpt, Q<:Qualia, C<:Clay}(io::IO, x::Flex{S,Q,C})
    tiesym = (Q==INEXACT) ? Inexactly : Exactly
    loDelim,hiDelim = delimiters(boundries(x)...)
    lo = @sprintf("%7.5g", x.lo); lo = strip(lo)
    hi = @sprintf("%7.5g", x.hi); hi = strip(hi)
    s = (x.lo != x.hi) ? string(loDelim, lo, tiesym, hi, hiDelim) : string(lo, tiesym)
    print(io, s)
end

function showmidrad{S<:Sculpt, Q<:Qualia, C<:Clay}(io::IO, x::Flex{S,Q,C})
    tiesym = (Q==INEXACT) ? Inexactly : Exactly
    md = @sprintf("%7.5g", mid(x)); md = strip(md)
    rd = @sprintf("%7.5g", rad(x)); rd = strip(rd)
    s = string(md, PlusMinus, rd)
    print(io, s)
end

showmidrad{S<:Sculpt, Q<:Qualia, C<:Clay}(x::Flex{S,Q,C}) = showmidrad(STDOUT,x)

function showmidpointradius{S<:Sculpt, Q<:Qualia, C<:Clay}(io::IO, x::Flex{S,Q,C})
    tiesym = (Q==INEXACT) ? Inexactly : Exactly
    md = mid(x)
    rd = rad(x)
    s  = string(md, PlusMinus, rd)
    print(io, s)
end

showmidpointradius{S<:Sculpt, Q<:Qualia, C<:Clay}(x::Flex{S,Q,C}) = showmidpointradius(STDOUT,x)


text:
@testset "OrderedRobinDict" begin

    @testset "Constructors" begin
        @test isa(@inferred(OrderedRobinDict()), OrderedRobinDict{Any,Any})
        @test isa(@inferred(OrderedRobinDict([(1,2.0)])), OrderedRobinDict{Int,Float64})
        @test isa(@inferred(OrderedRobinDict([("a",1),("b",2)])), OrderedRobinDict{String,Int})
        @test isa(@inferred(OrderedRobinDict(Pair(1, 1.0))), OrderedRobinDict{Int,Float64})
        @test isa(@inferred(OrderedRobinDict(Pair(1, 1.0), Pair(2, 2.0))), OrderedRobinDict{Int,Float64})
        @test isa(@inferred(OrderedRobinDict{Int,Float64}(Pair(1, 1), Pair(2, 2))), OrderedRobinDict{Int,Float64})
        @test isa(@inferred(OrderedRobinDict(Pair(1, 1.0), Pair(2, 2.0), Pair(3, 3.0))), OrderedRobinDict{Int,Float64})
        @test OrderedRobinDict(()) == OrderedRobinDict{Any,Any}()
        @test isa(@inferred(OrderedRobinDict([Pair(1, 1.0), Pair(2, 2.0)])), OrderedRobinDict{Int,Float64})
        @test_throws ArgumentError OrderedRobinDict([1,2,3,4])
        iter = Iterators.filter(x->x.first>1, [Pair(1, 1.0), Pair(2, 2.0), Pair(3, 3.0)])
        @test @inferred(OrderedRobinDict(iter)) == OrderedRobinDict{Int,Float64}(2=>2.0, 3=>3.0)
        iter = Iterators.drop(1:10, 1)
        @test_throws ArgumentError OrderedRobinDict(iter)
    end

    @testset "empty dictionary" begin
        d = OrderedRobinDict{Char, Int}()
        @test length(d) == 0
        @test isempty(d)
        @test_throws KeyError d['c'] == 1
        d['c'] = 1
        @test !isempty(d)
        @test_throws KeyError d[0.01]
        @test isempty(empty(d))
        empty!(d)
        @test isempty(d)

        # access, modification
        for c in 'a':'z'
            d[c] = c - 'a' + 1
        end

        @test (d['a'] += 1) == 2
        @test 'a' in keys(d)
        @test haskey(d, 'a')
        @test get(d, 'B', 0) == 0
        @test getkey(d, 'b', nothing) == 'b'
        @test getkey(d, 'B', nothing) == nothing
        @test !('B' in keys(d))
        @test !haskey(d, 'B')
        @test pop!(d, 'a') == 2

        @test collect(keys(d)) == collect('b':'z')
        @test collect(values(d)) == collect(2:26)
        @test collect(d) == [Pair(a,i) for (a,i) in zip('b':'z', 2:26)]
    end

    @testset "convert" begin
        d = OrderedRobinDict{Int,Float32}(i=>Float32(i) for i = 1:10)
        @test convert(OrderedRobinDict{Int,Float32}, d) === d
        dc = convert(OrderedRobinDict{Int,Float64}, d)
        @test dc !== d
        @test keytype(dc) == Int
        @test valtype(dc) == Float64
        @test keys(dc) == keys(d)
        @test collect(values(dc)) == collect(values(d))
    end

    @testset "Issue #60" begin
        od60 = OrderedRobinDict{Int,Int}()
        od60[1] = 2

        ranges = [2:5, 6:9, 10:13]
        for range in ranges
            for i = range
                od60[i] = i+1
            end
            for i = range
                delete!( od60, i )
            end
        end
        od60[14]=15

        @test od60[14] == 15
    end

    # #############################
    # Copied and modified from Base/test/dict.jl

    # OrderedRobinDict

    @testset "OrderedRobinDict{Int,Int}" begin
        h = OrderedRobinDict{Int,Int}()
        for i=1:10000
            h[i] = i+1
        end

        @test collect(h) == [Pair(x,y) for (x,y) in zip(1:10000, 2:10001)]

        for i=1:2:10000
            delete!(h, i)
        end
        for i=1:2:10000
            h[i] = i+1
        end

        for i=1:10000
            @test h[i]==i+1
        end

        for i=1:10000
            delete!(h, i)
        end
        @test isempty(h)

        h[77] = 100
        @test h[77]==100

        for i=1:10000
            h[i] = i+1
        end

        for i=1:2:10000
            delete!(h, i)
        end

        for i=10001:20000
            h[i] = i+1
        end

        for i=2:2:10000
            @test h[i]==i+1
        end

        for i=10000:20000
            @test h[i]==i+1
        end
    end

    @testset "OrderedRobinDict{Any,Any}" begin
        h = OrderedRobinDict{Any,Any}([("a", 3)])
        @test h["a"] == 3
        h["a","b"] = 4
        @test h["a","b"] == h[("a","b")] == 4
        h["a","b","c"] = 4
        @test h["a","b","c"] == h[("a","b","c")] == 4
    end

    @testset "KeyError" begin
        z = OrderedRobinDict()
        get_KeyError = false
        try
            z["a"]
        catch _e123_
            get_KeyError = isa(_e123_, KeyError)
        end
        @test get_KeyError
    end

    @testset "filter" begin
        _d = OrderedRobinDict([("a", 0)])
        v = [k for k in filter(x->length(x)==1, collect(keys(_d)))]
        @test isa(v, Vector{String})
    end

    @testset "from tuple/vector/pairs/tuple of pair 1" begin
        d = OrderedRobinDict(((1, 2), (3, 4)))
        d2 = OrderedRobinDict([(1, 2), (3, 4)])
        d3 = OrderedRobinDict(1 => 2, 3 => 4)
        d4 = OrderedRobinDict((1 => 2, 3 => 4))

        @test d[1] === 2
        @test d[3] === 4

        @test d == d2 == d3 == d4
        @test isa(d, OrderedRobinDict{Int,Int})
        @test isa(d2, OrderedRobinDict{Int,Int})
        @test isa(d3, OrderedRobinDict{Int,Int})
        @test isa(d4, OrderedRobinDict{Int,Int})

        h = OrderedRobinDict{Int, Char}(1=>'a')
        @test h[1] == 'a'
    end

    @testset "from tuple/vector/pairs/tuple of pair 2" begin
        d = OrderedRobinDict(((1, 2), (3, "b")))
        d2 = OrderedRobinDict([(1, 2), (3, "b")])
        d3 = OrderedRobinDict(1 => 2, 3 => "b")
        d4 = OrderedRobinDict((1 => 2, 3 => "b"))

        @test d2[1] === 2
        @test d2[3] == "b"

        ## TODO: tuple of tuples doesn't work for mixed tuple types
        # @test d == d2 == d3 == d4
        # @test isa(d, OrderedRobinDict{Int,Any})
        @test d2 == d3 == d4
        @test isa(d2, OrderedRobinDict{Int,Any})
        @test isa(d3, OrderedRobinDict{Int,Any})
        @test isa(d4, OrderedRobinDict{Int,Any})
    end

    @testset "from tuple/vector/pairs/tuple of pair 3" begin
        d = OrderedRobinDict(((1, 2), ("a", 4)))
        d2 = OrderedRobinDict([(1, 2), ("a", 4)])
        d3 = OrderedRobinDict(1 => 2, "a" => 4)
        d4 = OrderedRobinDict((1 => 2, "a" => 4))

        @test d2[1] === 2
        @test d2["a"] === 4

        ## TODO: tuple of tuples doesn't work for mixed tuple types
        # @test d == d2 == d3 == d4
        @test d2 == d3 == d4
        # @test isa(d, OrderedRobinDict{Any,Int})
        @test isa(d2, OrderedRobinDict{Any,Int})
        @test isa(d3, OrderedRobinDict{Any,Int})
        @test isa(d4, OrderedRobinDict{Any,Int})
    end

    @testset "from tuple/vector/pairs/tuple of pair 4" begin
        d = OrderedRobinDict(((1, 2), ("a", "b")))
        d2 = OrderedRobinDict([(1, 2), ("a", "b")])
        d3 = OrderedRobinDict(1 => 2, "a" => "b")
        d4 = OrderedRobinDict((1 => 2, "a" => "b"))

        @test d[1] === 2
        @test d["a"] == "b"

        @test d == d2 == d3 == d4
        @test isa(d, OrderedRobinDict{Any,Any})
        @test isa(d2, OrderedRobinDict{Any,Any})
        @test isa(d3, OrderedRobinDict{Any,Any})
        @test isa(d4, OrderedRobinDict{Any,Any})
    end

    @testset "first" begin
        @test_throws ArgumentError first(OrderedRobinDict())
        @test first(OrderedRobinDict([(:f, 2)])) == Pair(:f,2)
    end

    @testset "Issue #1821" begin
        d = OrderedRobinDict{String, Vector{Int}}()
        d["a"] = [1, 2]
        @test_throws MethodError d["b"] = 1
        @test isa(repr(d), AbstractString)  # check that printable without error
    end

    @testset "Issue #2344" begin
        bestkey(d, key) = key
        bestkey(d::AbstractDict{K,V}, key) where {K<:AbstractString,V} = string(key)
        bar(x) = bestkey(x, :y)
        @test bar(OrderedRobinDict([(:x, [1,2,5])])) == :y
        @test bar(OrderedRobinDict([("x", [1,2,5])])) == "y"
    end

    @testset "isequal" begin
        @test  isequal(OrderedRobinDict(), OrderedRobinDict())
        @test  isequal(OrderedRobinDict([(1, 1)]), OrderedRobinDict([(1, 1)]))
        @test !isequal(OrderedRobinDict([(1, 1)]), OrderedRobinDict())
        @test !isequal(OrderedRobinDict([(1, 1)]), OrderedRobinDict([(1, 2)]))
        @test !isequal(OrderedRobinDict([(1, 1)]), OrderedRobinDict([(2, 1)]))

        @test isequal(OrderedRobinDict(), sizehint!(OrderedRobinDict(),96))

        # Here is what currently happens when dictionaries of different types
        # are compared. This is not necessarily desirable. These tests are
        # descriptive rather than proscriptive.
        @test !isequal(OrderedRobinDict([(1, 2)]), OrderedRobinDict([("dog", "bone")]))
        @test isequal(OrderedRobinDict{Int,Int}(), OrderedRobinDict{AbstractString,AbstractString}())
    end

    @testset "data_in" begin
        # Generate some data to populate dicts to be compared
        data_in = [ (rand(1:1000), randstring(2)) for _ in 1:1001 ]

        # Populate the first dict
        d1 = OrderedRobinDict{Int, String}()
        for (k,v) in data_in
            d1[k] = v
        end
        data_in = collect(d1)
        # shuffle the data
        for i in 1:length(data_in)
            j = rand(1:length(data_in))
            data_in[i], data_in[j] = data_in[j], data_in[i]
        end
        # Inserting data in different (shuffled) order should result in
        # equivalent dict.
        d2 = OrderedRobinDict{Int, AbstractString}()
        for (k,v) in data_in
            d2[k] = v
        end

        @test  isequal(d1, d2)
        d3 = copy(d2)
        d4 = copy(d2)
        # Removing an item gives different dict
        delete!(d1, data_in[rand(1:length(data_in))][1])
        @test !isequal(d1, d2)
        # Changing a value gives different dict
        d3[data_in[rand(1:length(data_in))][1]] = randstring(3)
        !isequal(d1, d3)
        # Adding a pair gives different dict
        d4[1001] = randstring(3)
        @test !isequal(d1, d4)
    end

    @testset "get!" begin
        # get! (get with default values assigned to the given location)
        f(x) = x^2
        d = OrderedRobinDict(8 => 19)

        @test get!(d, 8, 5) == 19
        @test get!(d, 19, 2) == 2

        @test get!(d, 42) do  # d is updated with f(2)
            f(2)
        end == 4

        @test get!(d, 42) do  # d is not updated
            f(200)
        end == 4

        @test get(d, 13) do   # d is not updated
            f(4)
        end == 16

        @test d == OrderedRobinDict(8=>19, 19=>2, 42=>4)
    end

    @testset "Issue #5886" begin
        d5886 = OrderedRobinDict()
        for k5886 in 1:11
            d5886[k5886] = 1
        end
        for k5886 in keys(d5886)
            # undefined ref if not fixed
            d5886[k5886] += 1
        end
    end

    @testset "Issue #216" begin
        @test !isordered(Dict{Int, String})
        @test isordered(OrderedRobinDict{Int, String})
    end

    @testset "Test merging" begin
        a = OrderedRobinDict("foo"  => 0.0, "bar" => 42.0)
        b = OrderedRobinDict("フー" => 17, "バー" => 4711)
        @test isa(merge(a, b), OrderedRobinDict{String,Float64})
    end

    @testset "Issue #9295" begin
        d = OrderedRobinDict()
        @test push!(d, 'a'=> 1) === d
        @test d['a'] == 1
        @test push!(d, 'b' => 2, 'c' => 3) === d
        @test d['b'] == 2
        @test d['c'] == 3
        @test push!(d, 'd' => 4, 'e' => 5, 'f' => 6) === d
        @test d['d'] == 4
        @test d['e'] == 5
        @test d['f'] == 6
        @test length(d) == 6
    end

    @testset "Serialization" begin
        s = IOBuffer()
        od = OrderedRobinDict{Char,Int64}()
        for c in 'a':'e'
            od[c] = c-'a'+1
        end
        serialize(s, od)
        seek(s, 0)
        dd = deserialize(s)
        @test isa(dd, OrderedRobinDict{Char,Int64})
        @test dd == od
        close(s)
    end

    @testset "Issue #148" begin
        d148 = OrderedRobinDict(
                    :gps => [],
                    :direction => 1:8,
                    :weather => 1:10
            )

        d148_2 = OrderedRobinDict(
            :time => 1:10,
            :features => OrderedRobinDict(
                :gps => 1:5,
                :direction => 1:8,
                :weather => 1:10
            )
        )
    end

    @testset "Issue #400" begin
        @test filter(p->first(p) > 1, OrderedRobinDict(1=>2, 3=>4)) isa OrderedRobinDict
    end

    @testset "Test that OrderedRobinDict merge with combiner returns type OrderedRobinDict" begin
        @test merge(+, OrderedRobinDict(:a=>1, :b=>2), OrderedRobinDict(:b=>7, :c=>4)) == OrderedRobinDict(:a=>1, :b=>9, :c=>4)
        @test merge(+, OrderedRobinDict(:a=>1, :b=>2), Dict(:b=>7, :c=>4)) isa OrderedRobinDict
    end

    @testset "pop!" begin
        h = OrderedRobinDict{Int, Int}()
        for i in 1:10
            h[i] = i+1
        end
        for i in 10:-1:1
            @test pop!(h) == (i => i+1)
        end

        d = OrderedRobinDict(1=>'a', 2=>'b', 3=>'c')
        @test pop!(d, 4, 'e') == 'e'
    end

    @testset "filter!" begin
        h = OrderedRobinDict{Int, Int}()
        for i in 1:100
            h[i] = i
        end
        filter!(x->isodd(x.first), h)
        for i in 1:2:100
            @test h[i] == i
        end
    end

end # @testset OrderedRobinDict

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
include("transforms.jl")
include("kernel.jl")
include("model.jl")
include("prediction.jl")
include("statistics.jl")
using Plots
using Random
using Printf

Random.seed!(1234);
s=[-10 1 2 10 12 24 30 56 80 85 90 100]
s0 = [5]
K1 = K(s, s, 12, rbf)
z = sample(K1) #zero mean data with fixed covariance structure
mean = (minimum(z) < 0 ? abs(minimum(z)) + 12 : 12.1);
z = z .+  mean; # sample positive data with fixed mean
z = map(x -> invBoxCox(x, 2), z)
println(z)
X = ones(length(s), 1) #assume constant mean b/c of lack of covariate information
X0 = ones(length(s0), 1)
pλ = x -> 1 #let lambda be uniformly distributed (value of constant is immaterial, because of renormalization)
pθ = x -> 1 #let theta be uniformly distributed

if false #DEPRECATED
    pdf = z0 -> prediction(X, X0, s, s0, boxCox, boxCoxPrime, pθ, pλ, z, [z0], [10 14], [1 3])
    cdf = z0 ->  int1D(pdf, 0, z0)
    un = cdf(100)
    @printf("Unnormalized Integral of Bayesian Predictive Density Function: %f\n", un)
    cdfn = x -> cdf(x)/un #normalized cdf
    pdfn = x -> pdf(x)/un #normalized pdf
    @time begin
    med = bisection(x -> cdfn(x)-0.5, 1e-3, 12, 1e-3, 120)
    end
    @printf("Median of CDF: %f \n", med)
    x = range(.01, stop=20, length= 120)
    y= []
    for i = 1:120
        append!(y, prediction(X, X0, s, s0, boxCox, boxCoxPrime, pθ, pλ, z, [x[i]], [12 14], [1 3]))
    end
    display(plot(x, y, seriestype = :line))
    display(plot!([med], seriestype = :vline))

    @time begin
    h = confidence(pdfn, med)
    end
    @printf("interval width: %f \n: ", h)
    display(plot!([med-h, med+h], seriestype = :vline))
end
"""
Structure for storing statistics
"""
struct sts
    median::Float64
    intwidth::Float64
end
"""
Use single point deletion to validate model
"""
function cross_validate(X, X0, s, g, gprime, pθ, pλ, z, rangeθ, rangeλ)
    ss = Array{sts}(undef, length(z))
    num_pts = 400 #number of points used in plots
    endpt = 20 #right end-point of plot for the pdf
    x = collect(range(.001, stop = endpt, length = num_pts))
    X = repeat(x, outer = [1, length(z)]) #preallocate space for plotting data
    Y = Array{Float64}(undef, num_pts, length(z))
    for i=1:length(z)
        println(i)
        ind = [collect(1:i-1); collect(i+1:length(z))]
        @time begin
        pdf = z0 -> prediction(X[ind], [X[i]], s[ind], [s[i]], g, gprime, pθ, pλ, z[ind], [z0], rangeθ, rangeλ)#leave out one index
        cdf = z0 ->  int1D(pdf, 0, z0); un = cdf(50) #compute unnormalized cdf
        cdfn = x -> cdf(x)/un; pdfn = x -> pdf(x)/un  #compute normalized pdf and cdf
        end
        @time begin
        med = bisection(x -> cdfn(x)-0.5, 1e-3, 50, 1e-3, 15) #compute median
        end
        @time begin
        intwidth = confidence(pdfn, med) #compute width of credible interval corresponding to 95% confidence level
        end
        st = sts(med, intwidth); ss[i]=st
        @time begin
        for j = 1:num_pts
            Y[j, i] = pdfn(x[j])
        end
        end
    end
    return ss, X, Y
end

if true
    svals, X, Y = cross_validate(X, X0, s, boxCox, boxCoxPrime, pθ, pλ, z, [11.5 12.5], [1.95 2.05])
    ss = Array{Float64}(undef, 3, length(s))
    for j = 1:length(s)
        ss[1, j] = svals[j].median-svals[j].intwidth
        ss[2, j] = svals[j].median
        ss[3, j] = svals[j].median+svals[j].intwidth
    end
    display(plot(X, Y, layout=length(s), legend = nothing, ylims = (0, maximum(Y)+.03)))
    display(plot!(ss, seriestype = :vline, layout = length(s)))
    display(plot!(z', seriestype = :vline, layout=length(s)))
    savefig("validate18.png")
end

max_stars_repo_name: dbindel/btg

max_issues_repo_name: dbindel/btg

max_forks_repo_name: dbindel/btg


text:
module JFET

using FastGaussQuadrature
using MultiPoly

export ReferenceLine, ReferenceTriangle, ReferenceSquare,
       ReferenceTetrahedron, ReferenceCube, LagrangeFE

export dim, volume, make_quadrule, integrate, points,
        weights, vertices, edges, faces, facet, tabulate

#include("polybasis.jl")
include("refelements.jl")
include("finite_element.jl")
include("quadrature.jl")

end # module

max_stars_repo_name: JuliaTagBot/JFET

max_issues_repo_name: JuliaTagBot/JFET

max_forks_repo_name: JuliaTagBot/JFET


text:
for f in ["point", "bounds", "transform", "distance"]
    include("$f.jl")
end


text:
import Base: @deprecate, depwarn

abstract type MemoryLayout end

struct ColumnMajor <: MemoryLayout
    function ColumnMajor()
        depwarn("Use `Val(:col)` instead of `ColumnMajor()`", :ColumnMajor)
        new()
    end
end

struct RowMajor <: MemoryLayout
    function RowMajor()
        depwarn("Use `Val(:row)` instead of `RowMajor()`", :RowMajor)
        new()
    end
end

@deprecate RowMajor Val{:row}
@deprecate ColumnMajor Val{:col}

layout_map(orient) = typeof(orient) <: RowMajor ? Val(:row) : Val(:col)

col_warn = "Use `Val(:col)` instead of `ColumnMajor()`"
row_warn = "Use `Val(:row)` instead of `RowMajor()`"

function kernelmatrix!(
        σ::MemoryLayout,
        P::Matrix,
        κ::Kernel,
        X::AbstractMatrix,
        symmetrize::Bool = true
    )
    orientation = layout_map(σ)
    kernelmatrix!(orientation, P, κ, X, symmetrize)
end

function kernelmatrix!(
        σ::MemoryLayout,
        P::Matrix,
        κ::Kernel,
        X::AbstractMatrix,
        Y::AbstractMatrix
    )
    orientation = layout_map(σ)
    kernelmatrix!(orientation, P, κ, X, Y)
end

function kernelmatrix(
        σ::MemoryLayout,
        κ::Kernel,
        X::AbstractMatrix,
        symmetrize::Bool = true
    )
    orientation = layout_map(σ)
    kernelmatrix(orientation, κ, X, symmetrize)
end

function kernelmatrix(
        σ::MemoryLayout,
        κ::Kernel,
        X::AbstractMatrix,
        Y::AbstractMatrix
    )
    orientation = layout_map(σ)
    kernelmatrix(orientation, κ, X, Y)
end


text:
using MedicalTransforms
using Documenter

makedocs(;
    modules=[MedicalTransforms],
    authors="Dale <djblack@uci.edu> and contributors",
    repo="https://github.com/Dale-Black/MedicalTransforms.jl/blob/{commit}{path}#L{line}",
    sitename="MedicalTransforms.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://Dale-Black.github.io/MedicalTransforms.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/Dale-Black/MedicalTransforms.jl",
)


text:
using Delegate
using Base.Test

import Base: (==), (<), (<=), isequal, isless, abs, iseven, isodd, (-), (+), (*), hypot;

import Delegate: @delegate_1field1var, @delegate_1field2vars, @delegate_2fields1var,
                 @traject_1field1var,  @traject_1field2vars,  @traject_2fields1var;

type AnInt
   val::Int
end;

@delegate_1field1var(  AnInt, val, [ iseven, isodd ] );

@delegate_1field2vars( AnInt, val, [ (==), (<), (<=), isequal, isless ] );

@traject_1field1var(   AnInt, val, [ (-), abs ] );

@traject_1field2vars(  AnInt, val, [ (+), (-), (*) ] );

a = AnInt(3);
b = AnInt(7);

@test iseven(a) == false;
@test isodd(b)  == true;

@test -a != a;
@test -(-a) == a;
@test abs(-a) == a;

@test (a <  b) == true;
@test (b <= a) == false;

@test (a + b) == AnInt(10);
@test (a - b) == AnInt(-4);
@test (a * b) == AnInt(21);

typealias SysFloat Union{Float64,Float32};

function renormalize{T<:SysFloat}(a::T, b::T)
    hi = a + b
    t = hi - a
    lo = (a - (hi - t)) + (b - t)
    hi, lo
end;

type HiLo{T<:SysFloat}
   hi::T
   lo::T
end;

function (==){T<:SysFloat}(a::HiLo{T}, b::HiLo{T})
    (a.hi == b.hi) & (a.lo == b.lo)
end;

@traject_2fields1var( HiLo, hi, lo, [ renormalize, ] );

myHiLo = renormalize( HiLo(12.555555555, 8000.333333333) );
validHiLo = HiLo(8012.888888888, 4.440892098500626e-14);

@test myHiLo == validHiLo;


text:
"""
MasterOptions: A struct containing parameters for a MasterProblem

### Keyword Arguments
* `nlegs`: 1 (direct-route model) or 2 (single-transfer model)
* `nfreqs`: number of frequency levels to choose from
* `freqwts`: gamma in the paper; relative direct ridership coefficients
* `xfrwts`: lambda in the paper; relative transferring ridership coefficients
* `costwts`: rho in the paper; relative cost coefficients
* `constrainedg`: whether to constrain edges to be used at most once
* `ninitial`: the first `ninitial` lines in the `linelist` were used to start the master problem
* `angleparam`: old parameter; keep to -1.0
* `distparam`: Gamma in the paper; higher values allow longer commutes
"""
@with_kw mutable struct MasterOptions
    nlegs::Int                  = 1
    nfreqs::Int                 = 1
    freqwts::Vector{Float64}    = [1.0]
    xfrwts::Vector{Float64}     = [1.0]
    costwts::Vector{Float64}    = [1.0]
    constrainedg::Bool          = false
    ninitial::Int               = 0
    angleparam::Float64         = -1.0 # unused by default
    distparam::Float64          = 1.5
    modeltype::Symbol           = :lp
end

"""
MasterProblem

* `np`: a TransitNetwork object that contains the network information
* `linelist`: a list of lines to choose to operate
* `commutelines`: has length `options.nlegs`.
                  The first entry maps a commute (u, v) to the lines that contain both stops
                  The second entry is unused.
* `costs`: the cost for each line
* `model`: a JuMP model
* `budget`: the budget for the budget constraint, which is fixed by `NDM.optimize()`
* `x`: x[l, f] indicates whether a line l is operated at frequency level f
* `θ`: θ[u, v] indicates the fraction of commuters for commute (u, v) can be serviced
* `solver`: GurobiSolver by default
* `options`: struct containing options for the optimization
"""
mutable struct MasterProblem
    # Network information
    np::TransitNetwork
    linelist::Vector{Vector{Int}}
    commutelines::Vector{Dict{Tuple{Int,Int},Any}}
    costs::Vector{Float64}

    # Optimization information
    model::JuMP.Model
    budget::JuMP.Variable
    x::Array{JuMP.Variable,2}
    θ::JuMP.JuMPArray
    solver

    # Parameters
    options::MasterOptions
end

"""
Returns a MasterProblem object for the given TransitNetwork, initialized on `linelist`
with default parameters in `options`.
To start with an empty set of lines, use linelist = Vector{Vector{Int}}()
"""
function MasterProblem(
        np::TransitNetwork;
        linelist::Vector{Vector{Int}} = uniquelines(np.lines),
        options = MasterOptions(),
        solver = Gurobi.GurobiSolver(OutputFlag = 0)
    )
    @assert in(options.modeltype, [:lp, :ip])
    @assert in(np.gridtype, [:latlong, :euclidean])
    @assert abs(options.angleparam) <= 1.0
    @assert options.distparam >= 1.0
    @assert length(options.freqwts) == length(options.xfrwts) == length(options.costwts) == options.nfreqs

    commutelines = allcommutelines(np, options.nlegs, linelist)
    costs = Vector{Float64}([linecost(np, line) for line in linelist])
    rmp, budget, x, θ = mastermodel(
        np, linelist, commutelines, costs, solver, options
    )

    MasterProblem(
        np, linelist, commutelines, costs,
        rmp, budget, x, θ, solver, options
    )
end

"Builds base master problem model"
function mastermodel(
        np::TransitNetwork,
        linelist::Vector{Vector{Int}},
        commutelines::Vector{Dict{Tuple{Int,Int},Any}},
        costs::Vector{Float64},
        solver,
        options::MasterOptions
    )
    nlines = length(linelist)
    nfreqs = options.nfreqs
    freqwts = options.freqwts
    xfrwts = options.xfrwts
    costwts = options.costwts

    rmp = JuMP.Model(solver=solver)

    if options.modeltype == :lp
        JuMP.@variable(rmp, x[l=1:nlines,f=1:nfreqs] >= 0)
    elseif options.modeltype == :ip
        JuMP.@variable(rmp, x[l=1:nlines,f=1:nfreqs], Bin)
    end
    if options.nlegs == 2
        JuMP.@variable(rmp, twoline[
            (u,v)=commutes(np),
            w=np.xfrstns[(u,v)]
        ])
    end
    JuMP.@variable(rmp, budget)
    JuMP.@variable(rmp, θ[(u,v)=commutes(np)])
    # maximize ridership
    JuMP.@objective(rmp,
        Max,
        sum(demand(np,(u,v))*θ[(u,v)] for (u,v) in commutes(np))
    )
    # choice constraints
    JuMP.@constraint(rmp,
        choseub[(u,v)=commutes(np)],
        θ[(u,v)] <= 1
    )
    if options.nlegs == 1
        JuMP.@constraint(rmp,
            choseline[(u,v)=commutes(np)],
            θ[(u,v)] <= sum(freqwts[f]*sum(x[l,f] for l in commutelines[1][u,v])
                            for f in 1:nfreqs))
    else
        JuMP.@constraint(rmp,
            choseline[(u,v)=commutes(np)],
            θ[(u,v)] <=
                sum(freqwts[f]*sum(x[l,f] for l in commutelines[1][u,v]) for f in 1:nfreqs) +
                sum(twoline[(u,v),w] for w in np.xfrstns[(u,v)]))
        JuMP.@constraint(rmp,
            freq2a[(u,v)=commutes(np), w=np.xfrstns[(u,v)]],
            twoline[(u,v),w] <=
                sum(xfrwts[f]*sum(x[l,f] for l in setdiff(commutelines[1][u,w], commutelines[1][u,v]))
                    for f in 1:nfreqs))
        JuMP.@constraint(rmp,
            freq2b[(u,v)=commutes(np), w=np.xfrstns[(u,v)]],
            twoline[(u,v),w] <=
                sum(xfrwts[f]*sum(x[l,f] for l in setdiff(commutelines[1][w,v], commutelines[1][u,v]))
                    for f in 1:nfreqs))
    end

    if options.constrainedg
        edgelines = Dict{Tuple{Int,Int},Vector{Int}}()
        for l in 1:nlines, k in 2:length(linelist[l])
            u = linelist[l][k-1]
            v = linelist[l][k]
            edg = (min(u,v), max(u,v))
            if haskey(edgelines, edg)
                push!(edgelines[edg], l)
            else
                edgelines[edg] = [l]
            end
        end
        edges = keys(edgelines)
        coeffs = Dict{Tuple{Int,Int},Vector{Float64}}()
        for edg in edges
            if length(edgelines[edg]) == 1
                delete!(edgelines, edg)
            else
                coeffs[edg] = ones(length(edgelines[edg]))
                # allow for overlap between lines in original key network
                denom = length(intersect(1:options.ninitial, edgelines[edg]))
                if denom > 1
                    for l in 1:length(edgelines[edg])
                        if edgelines[edg][l] <= options.ninitial
                            coeffs[edg][l] = 1/denom
                        end
                    end
                end
            end
        end
        JuMP.@constraint(rmp,
            ccon[edg in keys(edgelines)],
            sum(sum(coeffs[edg][l]*x[edgelines[edg][l],f]
                    for l in 1:length(edgelines[edg]))
                for f in 1:nfreqs) <= 1)
    end

    # choose one frequency
    if nfreqs > 1
        JuMP.@constraint(rmp,
            [l=1:nlines],
            sum(x[l,f] for f in 1:nfreqs) <= 1)
    end

    # budget constraint
    JuMP.@constraint(rmp, bcon,
        sum(costwts[f]*costs[l]*x[l,f] for l in 1:nlines for f in 1:nfreqs) <= budget)

    rmp, budget, x, θ
end

"""
Adds a new column to the master problem.

It updates network information and creates a new model.
"""
function addcolumn!(
        rmp::MasterProblem,
        line::Vector{Int}
    )
    # recompute line information
    addline!(
        rmp.np, rmp.linelist, rmp.commutelines, line
    )
    push!(rmp.costs, linecost(rmp.np, line))
    initialbudget = JuMP.getvalue(rmp.budget)

    rmp.model, rmp.budget, rmp.x, rmp.θ =
        mastermodel(
            rmp.np, rmp.linelist, rmp.commutelines, rmp.costs, rmp.solver,
            rmp.options
        )
    JuMP.fix(rmp.budget, initialbudget)
end

"""
Returns a vector (with length `nlegs<=2`) of Dict{Tuple{Int,Int},Vector}, where

    commutelines[1][u,v] is a vector of the i lines that connect u to v, and
    commutelines[2][u,v] is empty (unused code)

"""
function allcommutelines(
        np::TransitNetwork,
        nlegs::Int,
        linelist::Vector{Vector{Int}}
    )
    @assert nlegs <= 2

    commutelines = [Dict{Tuple{Int,Int},Any}() for i in 1:nlegs]
    # 1. Initialize each (u,v) entry as an empty vector.
    for u=1:np.nstations, v=setdiff(1:np.nstations, u)
        commutelines[1][u,v] = Int[]
    end
    # 2. Populate the (u,v) entries in commutelines.
    #
    # We begin with an empty list, and iteratively add each line using the
    # addline!() method. At the end, we should recover the original linelist.
    linelistcopy = Vector{Int}[]
    for line in linelist
        linelistcopy = addline!(
            np, linelistcopy, commutelines, line
        )
    end
    @assert linelistcopy == linelist

    return commutelines
end

"""
Modify `oldlines` and `commutelines` in-place, adding information about `line`.
"""
function addline!(
        np::TransitNetwork,
        oldlines::Vector{Vector{Int}},
        commutelines::Vector{Dict{Tuple{Int,Int},Any}},
        line::Vector{Int}
    )
    l1 = length(oldlines) + 1 # We introduce a new `line` with index `l1`.
    # Single-leg commutes
    for u in line, v in line
        u != v && push!(commutelines[1][u,v], l1)
    end
    # Two-leg commutes are now handed directly in the master problem, so no need to update.
    # We update oldlines at the end to avoid conflicts with the earlier updates.
    push!(oldlines, line)
end

"""
Solves `mp::MasterProblem` with the given `budget`.

### Returns
The optimal solution `x` as a vector of floats.
"""
function optimize(mp::MasterProblem, budget::Real)
    JuMP.fix(mp.budget, budget)
    JuMP.solve(mp.model)
    JuMP.getvalue(mp.x)
end


text:
#Function for evaluating radial matrix elements
#23/07/2017

function radiel(atom,nn1,ll1,jj1,nn2,ll2,jj2)

  #Compute radial wavefunctions
    radial1, rscale1 = numerovfunc(atom,nn1,ll1,jj1)

    radial2, rscale2 = numerovfunc(atom,nn2,ll2,jj2)

  #Some preparation
  if nn1>=nn2
    Yk1 = radial1
    Yk2 = radial2
    rscale = rscale1
  else
    Yk1 = radial2
    Yk2 = radial1
    rscale = rscale2
  end

  #Resize smaller solution vector by attaching zeros to the end such that the two solution vectors are the same length
  if length(Yk1)≠length(Yk2)
    Yk2conc = append!(Yk2,zeros(length(Yk1)-length(Yk2)))
  else
    Yk2conc = Yk2
  end

  #Solve the matrix elements using method adapted from Zimmerman et al. (1979)
  deltar = rscale[2:end]-rscale[1:(end-1)]
  numervec = Yk1[2:end].*Yk2conc[2:end].*(rscale[2:end].^2).*deltar
  return matrixelement = sum(numervec)
end

max_stars_repo_name: CSChisholm/Rydberg

max_issues_repo_name: CSChisholm/Rydberg

max_forks_repo_name: CSChisholm/Rydberg


text:
####################################################################
function blrmodel(x, flux, σ2, idx; maxiterml = 10, maxitervi = 10, maxrandom=10, seed = 1, show_trace = true, S=S, Stest=0)
####################################################################

    M = numcentres = 30

    JITTER = 1e-8

    #--------------------------------------------------------
    # Auxiliary
    #--------------------------------------------------------

    rg = MersenneTwister(seed)

    Σdiag = Diagonal(σ2)

    numF = length(unique(idx))

    @assert(numF == length(galaxyvector(3000.0)))

    numparam = 2*numF + 2

    @assert(all(unique(idx) .== collect(1:numF)))

    @show N = length(x)

    @assert(N == length(flux) == length(idx) == length(σ2))

    @printf("Running blrmodel with %d and %d filters\n", N, numF)
    @printf("\tNumber of free parameters is %d\n", numparam)
    @printf("\tVI has S=%d samples and runs for %d iterations\n", S, maxitervi)

    #--------------------------------------------------------
    function unpack(param)
    #--------------------------------------------------------

        @assert(length(param) == numparam)

        local a       = param[1:numF]
        local b       = param[numF+1:2*numF]
        local logrinv = param[2*numF+1]
        local logα    = param[2*numF+2]

        a, b, logrinv, logα

    end

    #--------------------------------------------------------
    function objective(param)
    #--------------------------------------------------------

        local a, b, logrinv, logα = unpack(param)

        -1.0 * logl(a = a, logrinv = logrinv, b = b, logα = logα)

    end

    #--------------------------------------------------------
    function blockvector!(bvec, idx, b)
    #--------------------------------------------------------

        for i in 1:N
            bvec[i] = b[idx[i]]
        end

    end

    #--------------------------------------------------------
    function diagonalmatrix!(A, idx, a)
    #--------------------------------------------------------

        for i in 1:N
            A.diag[i] = a[idx[i]]
        end

    end

    #--------------------------------------------------------
    function logl(; a = a, b = b, logrinv = logrinv, logα = logα)
    #--------------------------------------------------------

        local Φ = PHI(x, centresongrid(x; numcentres=numcentres), exp(logrinv))

        local A = Diagonal(zeros(eltype(a), N))

        local bvec = zeros(eltype(a),N)

        blockvector!(bvec, idx, b)

        diagonalmatrix!(A, idx, a)

        local C = Σdiag + (A*Φ*Φ'*A) / exp(logα)

        logpdf(MvNormal(bvec, 0.5*(C+C')), flux)

    end

    #--------------------------------------------------------
    function loglΦ(; a = a, b = b, logα = logα, Φ = Φ)
    #--------------------------------------------------------

        local A = Diagonal(zeros(eltype(a), N))

        local bvec = zeros(eltype(a),N)

        blockvector!(bvec, idx, b)

        diagonalmatrix!(A, idx, a)

        local U = A*Φ / sqrt(exp(logα))

        local logd = logdet(Σdiag) + logdet(I + U'*(Σdiag\U))

        local Σinv = Σdiag\I

        local Cinv = Σinv - (Σinv*U)*((I + U'*Σinv*U)\(U'*Σinv))

         -0.5*logd -0.5*N*log(2π) -0.5*(bvec - flux)'*(Cinv*(bvec - flux))

    end

    #--------------------------------------------------------
    # Initialise GP hyperparameters with random search
    # and few iterations
    #--------------------------------------------------------

    initsolutions       = [randn(rg, numparam) for i=1:maxrandom]

    initialfitness      = @showprogress "Initial random search " map(objective, initsolutions)

    bestinitialsolution = initsolutions[argmin(initialfitness)]

    #--------------------------------------------------------
    # Call optimiser
    #--------------------------------------------------------

    opt    = Optim.Options(iterations = maxiterml, show_trace = show_trace, show_every=100)

    result = Optim.optimize(objective, bestinitialsolution, NelderMead(), opt)

    #--------------------------------------------------------
    # Retrieve optimised parameters
    #--------------------------------------------------------

    a, b, logrinv, logα = unpack(result.minimizer)

    @show exp(logα) 1/exp(logrinv)

    #--------------------------------------------------------
    # Call VI but fix length scale parameter, logα and basis Φ
    #--------------------------------------------------------

    Φ = PHI(x, centresongrid(x; numcentres=numcentres), exp(logrinv))

    A = Diagonal(a[idx])

    bvec = b[idx]

    # @show  logl(; a = a, b = b, logrinv = logrinv, logα = logα)
    # @show loglΦ(; a = a, b = b, logα = logα, Φ = Φ)
    # return 0,0

    function newobjective(θ)
        @assert(length(θ) == numF*2)
        local a = @view θ[1:numF]
        local b = @view θ[numF+1:2*numF]
        loglΦ(; a = a, b = b, logα = logα, Φ = Φ)
    end

    posterior, logevidence = VI(newobjective, [a;b], S=S, Stest=Stest, show_every=10, iterations=maxitervi, optimiser = :lbfgs)

    return posterior

    #--------------------------------------------------------
    # Reconstruct curves on test points
    #--------------------------------------------------------

    xtest  = collect(LinRange(minimum(x), maximum(x), 200))

    Φ      = PHI(x, centresongrid(x; numcentres=numcentres), exp(logrinv))

    Φxtest = PHI(xtest, centresongrid(x; numcentres=numcentres), exp(logrinv))

    Sninv  = exp(logα)*I + Φ'*A*inv(Σdiag)*A*Φ     # 3.81

    Mn     = Sninv \ (Φ'*A*inv(Σdiag)*(flux .- bvec)) # 3.84

    # predicted mean

    μ = vec(Φxtest*Mn)

    figure(101)

    cla()

    clr = ["b","g","r","m","k","c"]

    for (i,l) in enumerate(unique(idx))

        plot(xtest, a[i]*μ .+ b[i], clr[i]*"-")

        aux = findall(idx .== l)

        plot(x[aux], flux[aux], clr[i]*".")

    end

    return mu, Sigma

end


text:
@testset "Scenario Tree Errors" begin
    st = PH.ScenarioTree()
    @test_throws(ErrorException, PH.add_leaf(st, PH.root(st), -1.0))
end

@testset "Scenario Tree Functions" begin

    # First stage -- Creates root/first stage node automatically
    st = PH.ScenarioTree()

    # Second stage
    n1 = PH.add_node(st, PH.root(st))
    n2 = PH.add_node(st, PH.root(st))

    # Third stage
    sc0 = PH.add_leaf(st, n1, 0.25*0.75)
    sc1 = PH.add_leaf(st, n1, 0.75*0.75)
    sc2 = PH.add_leaf(st, n2, 0.25)

    @test length(n1.children) == 2
    (n3, n4) = n1.children
    if n3.id > n4.id
        temp = n4
        n4 = n3
        n3 = temp
    end
    @test length(n2.children) == 1
    n5 = first(n2.children) # only element in the set

    @test PH.last_stage(st) == PH.StageID(3)
    @test PH.is_leaf(st.root) == false
    @test PH.is_leaf(n1) == false
    @test PH.is_leaf(n2) == true # only having 1 child scenario means this is actually a leaf
    @test PH.is_leaf(n3) == true
    @test PH.is_leaf(n4) == true
    @test PH.is_leaf(n5) == true

    @test st.prob_map[sc0] == 0.25*0.75
    @test st.prob_map[sc1] == 0.75*0.75
    @test st.prob_map[sc2] == 0.25

    nodes = [st.root, n1, n2, n3, n4, n5]
    for n in nodes
        @test st.tree_map[n.id] == n
    end

    @test PH.scenarios(st) == Set([sc0, sc1, sc2])
    @test PH.scenario_bundle(st.root) == Set([sc0, sc1, sc2])
    @test PH.scenario_bundle(n1) == Set([sc0, sc1])
    @test PH.scenario_bundle(n2) == Set([sc2])
    @test PH.scenario_bundle(n3) == Set([sc0])
    @test PH.scenario_bundle(n4) == Set([sc1])
    @test PH.scenario_bundle(n5) == Set([sc2])

end

@testset "Scenario Tree Utilities" begin
    N = 10
    p = rand(N)
    p /= sum(p)
    st = PH.two_stage_tree(p)

    @test PH.last_stage(st) == PH.stid(2)
    @test length(st.tree_map) == N + 1
    for (scen, prob) in pairs(st.prob_map)
        @test prob == p[PH.value(scen) + 1]
    end

    st = PH.two_stage_tree(N)
    @test PH.last_stage(st) == PH.stid(2)
    @test length(st.tree_map) == N + 1
    for (scen, prob) in pairs(st.prob_map)
        @test prob == 1.0 / 10.0
    end
end


text:
module TestTransformer

using Test, MLJModels
using Tables, CategoricalArrays, Random
using ScientificTypes
using StatsBase

import MLJBase

#### FEATURE SELECTOR ####

@testset "Feat Selector" begin
    N = 100
    X = (Zn   = rand(N),
         Crim = rand(N),
         x3   = categorical(rand("YN", N)),
         x4   = categorical(rand("YN", N)))

    namesX   = Tables.schema(X).names |> collect
    selector = FeatureSelector()
    f,       = MLJBase.fit(selector, 1, X)

    @test f == namesX

    Xt = MLJBase.transform(selector, f, MLJBase.selectrows(X, 1:2))

    @test Set(Tables.schema(Xt).names) == Set(namesX)
    @test length(Xt.Zn) == 2

    selector = FeatureSelector(features=[:Zn, :Crim])
    f,       = MLJBase.fit(selector, 1, X)

    @test MLJBase.transform(selector, f, MLJBase.selectrows(X, 1:2)) ==
            MLJBase.select(X, 1:2, [:Zn, :Crim])

    infos = MLJBase.info_dict(selector)
    @test infos[:input_scitype]  == MLJBase.Table(MLJBase.Scientific)
    @test infos[:output_scitype] == MLJBase.Table(MLJBase.Scientific)
end

#  To be added with FeatureSelectorRule X = (n1=["a", "b", "a"], n2=["g", "g", "g"], n3=[7, 8, 9],
#               n4 =UInt8[3,5,10],  o1=[4.5, 3.6, 4.0], )
# MLJBase.schema(X)
# Xc = coerce(X,  :n1=>Multiclass, :n2=>Multiclass)

# t = Discretizer(features=[:o1, :n3, :n2, :n1])
# @test Xt.features == [:o1, :n3, :n2, :n1]
# @test Xt.is_ordinal == [true, false, false, false]
# @test Xt.A == [512 1 1 1; 1 2 1 2; 256 3 1 1]

#### UNIVARIATE DISCRETIZATION ####

@testset "U-Discr" begin
    v = randn(10000)
    t = UnivariateDiscretizer(n_classes=100);
    result, = MLJBase.fit(t, 1, v)
    w = MLJBase.transform(t, result, v)
    bad_values = filter(v - MLJBase.inverse_transform(t, result, w)) do x
        abs(x) > 0.05
    end
    @test length(bad_values)/length(v) < 0.06

    # scalars:
    @test MLJBase.transform(t, result, v[42]) == w[42]
    r =  MLJBase.inverse_transform(t, result, w)[43]
    @test MLJBase.inverse_transform(t, result, w[43]) ≈ r

    # test of permitted abuses of argument:
    @test MLJBase.inverse_transform(t, result, get(w[43])) ≈ r
    @test MLJBase.inverse_transform(t, result, map(get, w)) ≈
        MLJBase.inverse_transform(t, result, w)

    # all transformed vectors should have an identical pool (determined in
    # call to fit):
    v2 = v[1:3]
    w2 = MLJBase.transform(t, result, v2)
    @test levels(w2) == levels(w)

    #### UNIVARIATE STANDARDIZER ####

    stand = UnivariateStandardizer()
    f,    = MLJBase.fit(stand, 1, [0, 2, 4])

    @test round.(Int, MLJBase.transform(stand, f, [0,4,8])) == [-1.0,1.0,3.0]
    @test round.(Int, MLJBase.inverse_transform(stand, f, [-1, 1, 3])) == [0, 4, 8]

    infos = MLJBase.info_dict(stand)
end

#### STANDARDIZER ####

@testset "Standardizer" begin
    N = 5
    rand_char = rand("abcefgh", N)
    while length(unique(rand_char)) < 2
        rand_char = rand("abcefgh", N)
    end
    X = (OverallQual  = rand(UInt8, N),
         GrLivArea    = rand(N),
         Neighborhood = categorical(rand_char, ordered=true),
         x1stFlrSF    = sample(1:10, N, replace=false),
         TotalBsmtSF  = rand(N))

    # introduce a field of type `Char`:
    x1 = categorical(map(Char, (X.OverallQual |> collect)))

    X = (x1=x1, x2=X[2], x3=X[3], x4=X[4], x5=X[5])

    stand = Standardizer()
    f,    = MLJBase.fit(stand, 1, X)
    Xnew  = MLJBase.transform(stand, f, X)

    @test Xnew[1] == X[1]
    @test MLJBase.std(Xnew[2]) ≈ 1.0
    @test Xnew[3] == X[3]
    @test Xnew[4] == X[4]
    @test MLJBase.std(Xnew[5]) ≈ 1.0

    stand.features = [:x1, :x5]
    f,   = MLJBase.fit(stand, 1, X)
    Xnew = MLJBase.transform(stand, f, X)

    @test issubset(Set(keys(f)), Set(Tables.schema(X).names[[5,]]))

    Xt = MLJBase.transform(stand, f, X)

    @test Xnew[1] == X[1]
    @test Xnew[2] == X[2]
    @test Xnew[3] == X[3]
    @test Xnew[4] == X[4]
    @test MLJBase.std(Xnew[5]) ≈ 1.0

    # test on ignoring a feature, even if it's listed in the `features`
    stand.ignore = true
    f,   = MLJBase.fit(stand, 1, X)
    Xnew = MLJBase.transform(stand, f, X)

    @test issubset(Set(keys(f)), Set(Tables.schema(X).names[[2,]]))

    Xt = MLJBase.transform(stand, f, X)

    @test Xnew[1] == X[1]
    @test MLJBase.std(Xnew[2]) ≈ 1.0
    @test Xnew[3] == X[3]
    @test Xnew[4] == X[4]
    @test Xnew[5] == X[5]

    stand = Standardizer(features=[:x1, :mickey_mouse])
    @test_logs(
        (:warn, r"Some specified"),
        (:warn, r"No features left"),
        MLJBase.fit(stand, 1, X)
    )

    stand.ignore = true
    @test_logs (:warn, r"Some specified") MLJBase.fit(stand, 1, X)

    @test_throws ArgumentError Standardizer(ignore=true)

    stand = Standardizer(features=[:x3, :x4], count=true, ordered_factor=true)
    f,   = MLJBase.fit(stand, 1, X)
    Xnew = MLJBase.transform(stand, f, X)
    @test issubset(Set(keys(f)), Set(Tables.schema(X).names[3:4,]))

    Xt = MLJBase.transform(stand, f, X)

    @test Xnew[1] == X[1]
    @test Xnew[2] == X[2]
    @test elscitype(X[3]) <: OrderedFactor
    @test elscitype(Xnew[3]) <: Continuous
    @test MLJBase.std(Xnew[3]) ≈ 1.0
    @test elscitype(X[4]) == Count
    @test elscitype(Xnew[4]) <: Continuous
    @test MLJBase.std(Xnew[4]) ≈ 1.0
    @test Xnew[5] == X[5]

    stand = Standardizer(features= x-> x == (:x2))
    f,    = MLJBase.fit(stand, 1, X)
    Xnew  = MLJBase.transform(stand, f, X)

    @test Xnew[1] == X[1]
    @test MLJBase.std(Xnew[2]) ≈ 1.0
    @test Xnew[3] == X[3]
    @test Xnew[4] == X[4]
    @test Xnew[5] == X[5]

    infos = MLJBase.info_dict(stand)

    @test infos[:name] == "Standardizer"
    @test infos[:input_scitype] == MLJBase.Table(MLJBase.Scientific)
    @test infos[:output_scitype] == MLJBase.Table(MLJBase.Scientific)
end

#### UNIVARIATE BOX COX TRANSFORMER ####

@testset "U-boxcox" begin
    # create skewed non-negative vector with a zero value:
    Random.seed!(1551)
    v = abs.(randn(1000))
    v = v .- minimum(v)

    t  = UnivariateBoxCoxTransformer(shift=true)
    f, = MLJBase.fit(t, 2, v)

    e = v - MLJBase.inverse_transform(t, f, MLJBase.transform(t, f, v))
    @test sum(abs, e) <= 5000*eps()

    infos = MLJBase.info_dict(t)

    @test infos[:name] == "UnivariateBoxCoxTransformer"
    @test infos[:input_scitype] == AbstractVector{MLJBase.Continuous}
    @test infos[:output_scitype] == AbstractVector{MLJBase.Continuous}
end

#### ONE HOT ENCODER ####

@testset "One-Hot" begin
    X = (name   = categorical(["Ben", "John", "Mary", "John"], ordered=true),
         height = [1.85, 1.67, 1.5, 1.67],
         favourite_number = categorical([7, 5, 10, 5]),
         age    = [23, 23, 14, 23])

    t  = OneHotEncoder()
    f, _, report = @test_logs((:info, r"Spawning 3"),
                    (:info, r"Spawning 3"), MLJBase.fit(t, 1, X))

    Xt = MLJBase.transform(t, f, X)

    @test Xt.name__John == float.([false, true, false, true])
    @test Xt.height == X.height
    @test Xt.favourite_number__10 == float.([false, false, true, false])
    @test Xt.age == X.age
    @test MLJBase.schema(Xt).names == (:name__Ben, :name__John, :name__Mary,
                               :height, :favourite_number__5,
                               :favourite_number__7, :favourite_number__10,
                               :age)

    @test report.new_features == collect(MLJBase.schema(Xt).names)

    # test that *entire* pool of categoricals is used in fit, including
    # unseen levels:
    f, = @test_logs((:info, r"Spawning 3"), (:info, r"Spawning 3"),
                          MLJBase.fit(t, 1, MLJBase.selectrows(X,1:2)))
    Xtsmall = MLJBase.transform(t, f, X)
    @test Xt == Xtsmall

    # test that transform can be applied to subset of the data:
    @test MLJBase.transform(t, f, MLJBase.selectcols(X, [:name, :age])) ==
        MLJBase.selectcols(MLJBase.transform(t, f, X),
                           [:name__Ben, :name__John, :name__Mary, :age])

    # test ignore
    t = OneHotEncoder(features=[:name,], ignore=true)
    f, = MLJBase.fit(t, 0, X)
    Xt = MLJBase.transform(t, f, X)

    # test exclusion of ordered factors:
    t  = OneHotEncoder(ordered_factor=false)
    f, = MLJBase.fit(t, 0, X)
    Xt = MLJBase.transform(t, f, X)
    @test keys(Xt) == (:name, :height, :favourite_number__5,
                       :favourite_number__7, :favourite_number__10, :age)

    @test :name in Tables.schema(Xt).names
    @test :favourite_number__5 in Tables.schema(Xt).names
    @test MLJBase.schema(Xt).scitypes == (OrderedFactor{3}, Continuous,
                                          Continuous, Continuous,
                                          Continuous, Count)

    # test that one may not add new columns:
    X = (name = categorical(["Ben", "John", "Mary", "John"], ordered=true),
         height     = [1.85, 1.67, 1.5, 1.67],
         favourite_number = categorical([7, 5, 10, 5]),
         age        = [23, 23, 14, 23],
         gender     = categorical(['M', 'M', 'F', 'M']))
    @test_throws Exception MLJBase.transform(t, f, X)

    infos = MLJBase.info_dict(t)

    @test infos[:name] == "OneHotEncoder"
    @test infos[:input_scitype] == MLJBase.Table(MLJBase.Scientific)
    @test infos[:output_scitype] == MLJBase.Table(MLJBase.Scientific)
end

#### FILL IMPUTER ####

@testset "Imputer" begin
    X = (
        x = [missing,ones(10)...],
        y = [missing,ones(10)...],
        z = [missing,ones(10)...]
        )

    imp = FillImputer()
    f,  = MLJBase.fit(imp, 1, X)
    Xt  = MLJBase.transform(imp, f, X)
    @test all(.!ismissing.(Xt.x))
    @test Xt.x isa Vector{Float64} # no missing
    @test all(Xt.x .== 1.0)

    imp = FillImputer(features=[:x,:y])
    f,  = MLJBase.fit(imp, 1, X)
    @test_throws ErrorException MLJBase.transform(imp, f, X) # in X there's :z which we haven't trained on

    X = (x = categorical([missing, missing, missing, missing, "Old", "Young", "Middle", "Young",
                          "Old", "Young", "Middle", "Young"]), )

    mode_ = MLJBase.mode(["Old", "Young", "Middle", "Young", "Old", "Young", "Middle", "Young"])

    imp = FillImputer()
    f,  = MLJBase.fit(imp, 1, X)
    Xt  = MLJBase.transform(imp, f, X)
    @test all(.!ismissing.(Xt.x))
    @test all(Xt.x[ismissing.(X.x)] .== mode_)

    X  = (x = [missing, missing, 1, 1, 1, 1, 1, 5], )
    f, = MLJBase.fit(imp, 1, X)
    Xt = MLJBase.transform(imp, f, X)
    @test Xt.x == [1, 1, 1, 1, 1, 1, 1, 5]

    X = (x = categorical([missing, missing, missing, missing, "Old", "Young", "Middle", "Young",
                          "Old", "Young", "Middle", "Young"]),
         y = [missing, ones(11)...],
         z = [missing, missing, 1,1,1,1,1,5,1,1,1,1],
         a = rand("abc", 12))

    f, = MLJBase.fit(imp, 1, X)
    Xt = MLJBase.transform(imp, f, X)

    @test all(.!ismissing.(Xt.x))
    @test all(.!ismissing.(Xt.y))
    @test all(.!ismissing.(Xt.z))
    @test all(.!ismissing.(Xt.a))

    @test Xt.x[1] == mode_
    @test Xt.y[1] == 1
    @test Xt.z[1] == 1
end

#### CONTINUOUS ENCODER ####

@testset "Continuous encoder" begin

    X = (name  = categorical(["Ben", "John", "Mary", "John"], ordered=true),
         height = [1.85, 1.67, 1.5, 1.67],
         rubbish = ["a", "b", "c", "a"],
         favourite_number = categorical([7, 5, 10, 5]),
         age    = [23, 23, 14, 23])

    t  = ContinuousEncoder()
    f, _, _ = @test_logs((:info, r"Some.*dropped\:.*\:rubbish\]"),
                              MLJBase.fit(t, 1, X))

    Xt = MLJBase.transform(t, f, X)
    @test scitype(Xt) <: MLJBase.Table(MLJBase.Continuous)
    s = MLJBase.schema(Xt)
    @test s.names == (:name, :height, :favourite_number__5,
                      :favourite_number__7, :favourite_number__10, :age)

    t  = ContinuousEncoder(drop_last=true, one_hot_ordered_factors=true)
    f, _, r = MLJBase.fit(t, 0, X)
    Xt = MLJBase.transform(t, f, X)
    @test scitype(Xt) <: MLJBase.Table(MLJBase.Continuous)
    s = MLJBase.schema(Xt)
    @test s.names == (:name__Ben, :name__John, :height, :favourite_number__5,
                      :favourite_number__7, :age)

end

end
true