text:
# Rendering SQL.

mutable struct RenderContext <: IO
    dialect::SQLDialect
    io::IOBuffer
    level::Int
    nested::Bool
    vars::Vector{Symbol}

    RenderContext(dialect) =
        new(dialect, IOBuffer(), 0, false, Symbol[])
end

Base.write(ctx::RenderContext, octet::UInt8) =
    write(ctx.io, octet)

Base.unsafe_write(ctx::RenderContext, input::Ptr{UInt8}, nbytes::UInt) =
    unsafe_write(ctx.io, input, nbytes)

function newline(ctx::RenderContext)
    print(ctx, "\n")
    for k = 1:ctx.level
        print(ctx, "  ")
    end
end

function render(c::AbstractSQLClause; dialect = :default)
    ctx = RenderContext(dialect)
    render(ctx, convert(SQLClause, c))
    sql = String(take!(ctx.io))
    SQLStatement(sql = sql, dialect = ctx.dialect, vars = ctx.vars)
end

function render(ctx, name::Symbol)
    s = string(name)
    lq, rq = ctx.dialect.identifier_quotes
    if rq in s
        s = replace(s, rq => rq * rq)
    end
    print(ctx, lq, s, rq)
end

function render(ctx, names::Vector{Symbol})
    first = true
    for name in names
        if !first
            print(ctx, ", ")
        else
            first = false
        end
        render(ctx, name)
    end
end

render(ctx, ::Missing) =
    print(ctx, "NULL")

render(ctx, val::Bool) =
    if ctx.dialect.has_boolean_literals
        print(ctx, val ? "TRUE" : "FALSE")
    else
        print(ctx, val ? "(1 = 1)" : "(1 = 0)")
    end

render(ctx, val::Number) =
    print(ctx, val)

function render(ctx, val::AbstractString)
    if '\'' in val
        val = replace(val, '\'' => "''")
    end
    print(ctx, '\'', val, '\'')
end

render(ctx, val::Dates.AbstractTime) =
    print(ctx, '\'', val, '\'')

function render(ctx, c::SQLClause)
    render(ctx, c[])
    nothing
end

function render(ctx, cs::AbstractVector{SQLClause}, sep = nothing)
    first = true
    for c in cs
        if !first
            if @dissect(c, KW())
                print(ctx, ' ')
            elseif sep === nothing
                print(ctx, ", ")
            else
                print(ctx, ' ', sep, ' ')
            end
        else
            first = false
        end
        render(ctx, c)
    end
end

function render_lines(ctx, cs::AbstractVector{SQLClause}, sep = nothing)
    !isempty(cs) || return
    if length(cs) == 1
        print(ctx, ' ')
        render(ctx, cs[1])
    else
        ctx.level += 1
        newline(ctx)
        first = true
        for c in cs
            if !first
                if sep === nothing
                    print(ctx, ',')
                    newline(ctx)
                else
                    print(ctx, ' ', sep)
                    newline(ctx)
                end
            else
                first = false
            end
            render(ctx, c)
        end
        ctx.level -= 1
    end
end

function render(ctx, c::AggregateClause)
    if c.filter !== nothing || c.over !== nothing
        print(ctx, '(')
    end
    print(ctx, c.name, '(')
    if c.distinct
        print(ctx, "DISTINCT ")
    end
    render(ctx, c.args)
    print(ctx, ')')
    if c.filter !== nothing
        print(ctx, " FILTER (WHERE ")
        render(ctx, c.filter)
        print(ctx, ')')
    end
    if c.over !== nothing
        print(ctx, " OVER (")
        render(ctx, c.over)
        print(ctx, ')')
    end
    if c.filter !== nothing || c.over !== nothing
        print(ctx, ')')
    end
end

function render(ctx, c::AsClause)
    over = c.over
    if @dissect(over, PARTITION())
        render(ctx, c.name)
        print(ctx, " AS (")
        render(ctx, over)
        print(ctx, ')')
    elseif over !== nothing
        render(ctx, over)
        print(ctx, " AS ")
        render(ctx, c.name)
    end
end

function render(ctx, c::CaseClause)
    print(ctx, "(CASE")
    nargs = length(c.args)
    for (i, arg) in enumerate(c.args)
        if isodd(i)
            print(ctx, i < nargs ? " WHEN " : " ELSE ")
        else
            print(ctx, " THEN ")
        end
        render(ctx, arg)
    end
    print(ctx, " END)")
end

function render(ctx, c::CTEClause)
    over = c.over
    columns = c.columns
    render(ctx, c.name)
    if columns !== nothing
        print(ctx, " (")
        render(ctx, columns)
        print(ctx, ')')
    end
    if over !== nothing
        print(ctx, " AS ")
        if c.materialized === true
            print(ctx, "MATERIALIZED ")
        elseif c.materialized === false
            print(ctx, "NOT MATERIALIZED ")
        end
        nested = ctx.nested
        ctx.nested = true
        render(ctx, over)
        ctx.nested = nested
    end
end

function render(ctx, c::FromClause)
    newline(ctx)
    print(ctx, "FROM")
    over = c.over
    if over !== nothing
        print(ctx, ' ')
        render(ctx, over)
    end
end

function render(ctx, c::FunctionClause)
    print(ctx, c.name, '(')
    render(ctx, c.args)
    print(ctx, ')')
end

function render(ctx, c::GroupClause)
    over = c.over
    if over !== nothing
        render(ctx, over)
    end
    !isempty(c.by) || return
    newline(ctx)
    print(ctx, "GROUP BY")
    render_lines(ctx, c.by)
end

function render(ctx, c::HavingClause)
    over = c.over
    if over !== nothing
        render(ctx, over)
    end
    newline(ctx)
    print(ctx, "HAVING")
    if @dissect(c.condition, OP(name = :AND, args = args)) && length(args) >= 2
        render_lines(ctx, args, "AND")
    else
        print(ctx, ' ')
        render(ctx, c.condition)
    end
end

function render(ctx, c::IdentifierClause)
    over = c.over
    if over !== nothing
        render(ctx, over)
        print(ctx, '.')
    end
    render(ctx, c.name)
end

function render(ctx, c::JoinClause)
    over = c.over
    if over !== nothing
        render(ctx, over)
    end
    newline(ctx)
    cross = !c.left && !c.right && @dissect(c.on, LIT(val = true))
    if cross
        print(ctx, "CROSS JOIN ")
    elseif c.left && c.right
        print(ctx, "FULL JOIN ")
    elseif c.left
        print(ctx, "LEFT JOIN ")
    elseif c.right
        print(ctx, "RIGHT JOIN ")
    else
        print(ctx, "JOIN ")
    end
    if c.lateral
        print(ctx, "LATERAL ")
    end
    render(ctx, c.joinee)
    if !cross
        print(ctx, " ON ")
        render(ctx, c.on)
    end
end

function render(ctx, c::KeywordClause)
    print(ctx, c.name)
    over = c.over
    if over !== nothing
        print(ctx, ' ')
        render(ctx, over)
    end
end

function render(ctx, c::LimitClause)
    over = c.over
    if over !== nothing
        render(ctx, over)
    end
    start = c.offset
    count = c.limit
    start !== nothing || count !== nothing || return
    if ctx.dialect.limit_style === LIMIT_STYLE.MYSQL
        newline(ctx)
        print(ctx, "LIMIT ")
        if start !== nothing
            print(ctx, start, ", ")
        end
        print(ctx, count !== nothing ? count : "18446744073709551615")
    elseif ctx.dialect.limit_style === LIMIT_STYLE.SQLITE
        newline(ctx)
        print(ctx, "LIMIT ", count !== nothing ? count : -1)
        if start !== nothing
            newline(ctx)
            print(ctx, "OFFSET ", start)
        end
    else
        if start !== nothing
            newline(ctx)
            print(ctx, "OFFSET ", start, start == 1 ? " ROW" : " ROWS")
        end
        if count !== nothing
            newline(ctx)
            print(ctx, start === nothing ? "FETCH FIRST " : "FETCH NEXT ",
                       count, count == 1 ? " ROW" : " ROWS",
                       c.with_ties ? " WITH TIES" : " ONLY")
        end
    end
end

render(ctx, c::LiteralClause) =
    render(ctx, c.val)

function render(ctx, c::OperatorClause)
    if isempty(c.args)
        print(ctx, c.name)
    elseif length(c.args) == 1
        print(ctx, '(', c.name, ' ')
        render(ctx, c.args[1])
        print(ctx, ')')
    else
        print(ctx, '(')
        render(ctx, c.args, c.name)
        print(ctx, ')')
    end
end

function render(ctx, c::OrderClause)
    over = c.over
    if over !== nothing
        render(ctx, over)
    end
    !isempty(c.by) || return
    newline(ctx)
    print(ctx, "ORDER BY")
    render_lines(ctx, c.by)
end

function render(ctx, m::FrameMode)
    if m == FRAME_MODE.RANGE
        print(ctx, "RANGE")
    elseif m == FRAME_MODE.ROWS
        print(ctx, "ROWS")
    elseif m == FRAME_MODE.GROUPS
        print(ctx, "GROUPS")
    else
        throw(DomainError(m))
    end
end

function render(ctx, e::FrameExclusion)
    if e == FRAME_EXCLUSION.NO_OTHERS
        print(ctx, "EXCLUDE NO OTHERS")
    elseif e == FRAME_EXCLUSION.CURRENT_ROW
        print(ctx, "EXCLUDE CURRENT ROW")
    elseif e == FRAME_EXCLUSION.GROUP
        print(ctx, "EXCLUDE GROUP")
    elseif e == FRAME_EXCLUSION.TIES
        print(ctx, "EXCLUDE TIES")
    else
        throw(DomainError(e))
    end
end

function render_frame_endpoint(ctx, val)
    if iszero(val)
        print(ctx, "CURRENT ROW")
    else
        pos = val > zero(val)
        inf = isinf(val)
        if pos && inf
            print(ctx, "UNBOUNDED FOLLOWING")
        elseif pos
            render(ctx, val)
            print(ctx, " FOLLOWING")
        elseif !inf
            render(ctx, - val)
            print(ctx, " PRECEDING")
        else
            print(ctx, "UNBOUNDED PRECEDING")
        end
    end
end

function render(ctx, f::PartitionFrame)
    render(ctx, f.mode)
    print(ctx, ' ')
    if f.finish === nothing
        render_frame_endpoint(ctx, something(f.start, -Inf))
    else
        print(ctx, "BETWEEN ")
        render_frame_endpoint(ctx, something(f.start, -Inf))
        print(ctx, " AND ")
        render_frame_endpoint(ctx, f.finish)
    end
    if f.exclude !== nothing
        print(ctx, ' ')
        render(ctx, f.exclude)
    end
end

function render(ctx, c::PartitionClause)
    need_space = false
    if c.over !== nothing
        render(ctx, c.over)
        need_space = true
    end
    if !isempty(c.by)
        if need_space
            print(ctx, ' ')
        end
        print(ctx, "PARTITION BY ")
        render(ctx, c.by)
        need_space = true
    end
    if !isempty(c.order_by)
        if need_space
            print(ctx, ' ')
        end
        print(ctx, "ORDER BY ")
        render(ctx, c.order_by)
        need_space = true
    end
    if c.frame !== nothing
        if need_space
            print(ctx, ' ')
        end
        render(ctx, c.frame)
        need_space = true
    end
end

function render(ctx, c::SelectClause)
    nested = ctx.nested
    if nested
        ctx.level += 1
        print(ctx, '(')
        newline(ctx)
    end
    ctx.nested = true
    print(ctx, "SELECT")
    top = c.top
    if top !== nothing
        print(ctx, " TOP ", top.limit)
        if top.with_ties
            print(ctx, " WITH TIES")
        end
    end
    if c.distinct
        print(ctx, " DISTINCT")
    end
    render_lines(ctx, c.args)
    over = c.over
    if over !== nothing
        render(ctx, over)
    end
    ctx.nested = nested
    if nested
        ctx.level -= 1
        newline(ctx)
        print(ctx, ')')
    end
end

function render(ctx, c::SortClause)
    over = c.over
    if over !== nothing
        render(ctx, over)
    end
    if c.value == VALUE_ORDER.ASC
        print(ctx, " ASC")
    elseif c.value == VALUE_ORDER.DESC
        print(ctx, " DESC")
    end
    if c.nulls == NULLS_ORDER.NULLS_FIRST
        print(ctx, " NULLS FIRST")
    elseif c.nulls == NULLS_ORDER.NULLS_LAST
        print(ctx, " NULLS LAST")
    end
end

function render(ctx, c::UnionClause)
    nested = ctx.nested
    if nested
        ctx.level += 1
        print(ctx, '(')
        newline(ctx)
    end
    ctx.nested = false
    over = c.over
    if over !== nothing
        render(ctx, over)
    end
    for arg in c.args
        newline(ctx)
        print(ctx, "UNION")
        if c.all
            print(ctx, " ALL")
        end
        newline(ctx)
        render(ctx, arg)
    end
    ctx.nested = nested
    if nested
        ctx.level -= 1
        newline(ctx)
        print(ctx, ')')
    end
end

function render(ctx, c::VariableClause)
    style = ctx.dialect.variable_style
    prefix = ctx.dialect.variable_prefix
    pos = nothing
    if style != VARIABLE_STYLE.POSITIONAL
        pos = findfirst(==(c.name), ctx.vars)
    end
    if pos === nothing
        push!(ctx.vars, c.name)
        pos = length(ctx.vars)
    end
    print(ctx, prefix)
    if style == VARIABLE_STYLE.NAMED
        print(ctx, c.name)
    elseif style == VARIABLE_STYLE.NUMBERED
        print(ctx, pos)
    end
end

function render(ctx, c::WhereClause)
    over = c.over
    if over !== nothing
        render(ctx, over)
    end
    newline(ctx)
    print(ctx, "WHERE")
    if @dissect(c.condition, OP(name = :AND, args = args)) && length(args) >= 2
        render_lines(ctx, args, "AND")
    else
        print(ctx, ' ')
        render(ctx, c.condition)
    end
end

function render(ctx, c::WindowClause)
    over = c.over
    if over !== nothing
        render(ctx, over)
    end
    !isempty(c.args) || return
    newline(ctx)
    print(ctx, "WINDOW")
    render_lines(ctx, c.args)
end

function render(ctx, c::WithClause)
    if !isempty(c.args)
        print(ctx, "WITH ")
        if c.recursive && ctx.dialect.has_with_recursive
            print(ctx, "RECURSIVE ")
        end
        first = true
        for arg in c.args
            if !first
                print(ctx, ",")
                newline(ctx)
            else
                first = false
            end
            render(ctx, arg)
        end
        newline(ctx)
    end
    if c.over !== nothing
        render(ctx, c.over)
    end
end


text:
using CUDAdrv

using Test

@testset "CUDAdrv" begin

include("util.jl")

include("base.jl")

if CUDAdrv.configured
    @test length(devices()) > 0
    if length(devices()) > 0
        @test CuCurrentContext() == nothing

        # pick most recent device (based on compute capability)
        global dev = nothing
        for newdev in devices()
            if dev == nothing || capability(newdev) > capability(dev)
                dev = newdev
            end
        end
        @info "Testing using device $(name(dev))"

        global ctx = CuContext(dev)
        @test CuCurrentContext() != nothing

        @testset "API wrappers" begin
            include("errors.jl")
            include("version.jl")
            include("devices.jl")
            include("context.jl")
            include("module.jl")
            include("memory.jl")
            include("stream.jl")
            include("execution.jl")
            include("events.jl")
            include("profile.jl")
        end

        include("gc.jl")

        include("examples.jl")
    end
else
    @warn "CUDAdrv.jl has not been configured; skipping most tests."
end

end

ext: jl


text:
function add_constrained_var(model, set::MOI.AbstractScalarSet)
    vi, con_idx = MOI.add_constrained_variable(model, set)
    return MOI.SingleVariable(vi)
end

function add_constrained_var(model, set::MOI.AbstractVectorSet)
    vis, con_idx = MOI.add_constrained_variables(model, set)
    return MOI.SingleVariable.(vis)
end

"""
    generate_mip(optimizer, m, nlp_sol, tabu_list, start_fpump)

Generate a mip using the linear constraints of the original model
TODO: This can include quadratic constraints when the mip_solver supports them
Minimize the distance to nlp_sol and avoid using solutions inside the tabu list
"""
function generate_mip(optimizer, m, nlp_sol, tabu_list, start_fpump)
    mip_optimizer = MOI.instantiate(m.mip_solver, with_bridge_type=Float64)
    index_map = MOI.copy_to(mip_optimizer, NoObjectiveFilter(LinearFilter(optimizer)); copy_names=false)
    x = [index_map[vi] for vi in MOI.get(optimizer, MOI.ListOfVariableIndices())]
    mx = MOI.SingleVariable.(x)
    for i in 1:m.num_var
        if m.var_type[i] == :Bin
            if m.l_var[i] > 0
                set_bounds(mip_optimizer, x[i], 1.0, 1.0)
            elseif m.u_var[i] < 1
                set_bounds(mip_optimizer, x[i], 0.0, 0.0)
            end
        end
    end

    mabsx = add_constrained_var(mip_optimizer, MOI.Nonnegatives(m.num_disc_var))
    for (mabsxi, vi) in zip(mabsx, m.disc2var_idx)
        MOI.add_constraint(
            mip_optimizer,
            MOIU.operate(vcat, Float64, mabsxi, mx[vi] - nlp_sol[vi]),
            MOI.NormOneCone(2)
        )
    end

    # How long is the tabu list
    num_sols = 0
    for i=1:tabu_list.length
        if !isnan(tabu_list.sols[i][1])
            num_sols += 1
        else
            break
        end
    end

    # If there solutions in the tabu list => avoid them
    if num_sols > 0
        z1 = [add_constrained_var(mip_optimizer, MOI.ZeroOne()) for _=1:m.num_disc_var, _=1:num_sols]
        z2 = [add_constrained_var(mip_optimizer, MOI.ZeroOne()) for _=1:m.num_disc_var, _=1:num_sols]
        v = tabu_list.sols
        for k=1:num_sols, j=1:m.num_disc_var
            i = m.disc2var_idx[j]
            lbi = m.l_var[i] > typemin(Int64) ? m.l_var[i] : typemin(Int64)
            ubi = m.u_var[i] < typemax(Int64) ? m.u_var[i] : typemax(Int64)
            MOI.add_constraint(mip_optimizer, 1.0z1[j,k] + 1.0z2[j,k], MOI.LessThan(1.0))
            MOI.add_constraint(mip_optimizer, MA.@rewrite((lbi - v[k][j]) * z1[j,k] + z2[j,k] - mx[i]), MOI.LessThan(-v[k][j]))
            MOI.add_constraint(mip_optimizer, MA.@rewrite((ubi - v[k][j]) * z2[j,k] - z1[j,k] - mx[i]), MOI.GreaterThan(-v[k][j]))
        end
        for k=1:num_sols
            MOI.add_constraint(mip_optimizer, sum(1.0z1) + sum(1.0z2), MOI.GreaterThan(1.0))
        end
    end

    MOI.set(mip_optimizer, MOI.ObjectiveSense(), MOI.MIN_SENSE)
    obj = sum(1.0mabsx)
    MOI.set(mip_optimizer, MOI.ObjectiveFunction{typeof(obj)}(), obj)

    # Break the mip solver if it takes too long or throw a warning when this option isn't available
    current_time = time()-start_fpump
    time_left = m.options.feasibility_pump_time_limit-current_time
    time_left < 0 && (time_left = 1.0)

    # set time limit if supported
    old_time_limit = set_time_limit!(mip_optimizer, time_left)

    MOI.optimize!(mip_optimizer)
    status = MOI.get(mip_optimizer, MOI.TerminationStatus())

    # reset time limit
    set_time_limit!(mip_optimizer, old_time_limit)

    obj_val = NaN
    values = fill(NaN,m.num_var)
    if state_is_optimal(status; allow_almost=m.options.allow_almost_solved)
        obj_val = MOI.get(mip_optimizer, MOI.ObjectiveValue())

        # round mip values
        values = MOI.get(mip_optimizer, MOI.VariablePrimal(), x)
        for i=1:m.num_disc_var
            vi = m.disc2var_idx[i]
            values[vi] = round(values[vi])
        end
    end

    return status, values, obj_val
end

"""
    generate_nlp(optimizer, m, mip_sol, start_fpump; random_start=false)

Generates the original nlp but changes the objective to minimize the distance to the mip solution
"""
function generate_nlp(optimizer, m, mip_sol, start_fpump; random_start=false)
    nlp_optimizer = MOI.instantiate(m.nl_solver, with_bridge_type=Float64)
    index_map = MOI.copy_to(nlp_optimizer, NoObjectiveFilter(IntegerRelaxation(optimizer)); copy_names=false)
    x = [index_map[vi] for vi in MOI.get(optimizer, MOI.ListOfVariableIndices())]
    nx = MOI.SingleVariable.(x)
    if random_start
        restart_values = generate_random_restart(m)
    else
        restart_values = mip_sol
    end
    MOI.set(nlp_optimizer, MOI.VariablePrimalStart(), x, restart_values)

    MOI.set(nlp_optimizer, MOI.ObjectiveSense(), MOI.MIN_SENSE)
    obj = sum((nx[m.disc2var_idx[i]] - mip_sol[m.disc2var_idx[i]])^2 for i=1:m.num_disc_var)
    MOI.set(nlp_optimizer, MOI.ObjectiveFunction{typeof(obj)}(), obj)

    current_time = time()-start_fpump
    time_left = m.options.feasibility_pump_time_limit-current_time
    time_left < 0 && (time_left = 1.0)

    # set time limit if supported
    old_time_limit = set_time_limit!(nlp_optimizer, time_left)

    MOI.optimize!(nlp_optimizer)
    status = MOI.get(nlp_optimizer, MOI.TerminationStatus())

    set_time_limit!(nlp_optimizer, old_time_limit)

    nlp_obj = NaN
    nlp_sol = fill(NaN, m.num_var)
    if state_is_optimal(status; allow_almost=m.options.allow_almost_solved)
        nlp_obj = MOI.get(nlp_optimizer, MOI.ObjectiveValue())
        nlp_sol = MOI.get(nlp_optimizer, MOI.VariablePrimal(), x)
    end

    return status, nlp_sol, nlp_obj
end

"""
    generate_real_nlp(optimizer, m, sol; random_start=false)

Generate the original nlp and get the objective for that
"""
function generate_real_nlp(optimizer, m, sol; random_start=false)
    if m.num_var == m.num_disc_var
        return MOI.OPTIMAL, sol, evaluate_objective(optimizer, m, sol)
    end

    nlp_optimizer = MOI.instantiate(m.nl_solver, with_bridge_type=Float64)
    vis = collect(MOI.get(optimizer, MOI.ListOfVariableIndices()))
    disc_vals = Dict(vis[i] => sol[i] for i in m.disc2var_idx)
    index_map = MOI.copy_to(nlp_optimizer, FixVariables(IntegerRelaxation(optimizer), disc_vals); copy_names=false)
    x = [index_map[vi] for vi in MOI.get(optimizer, MOI.ListOfVariableIndices())]
    rx = MOI.SingleVariable.(x)

    if random_start
        restart_values = generate_random_restart(m)
        for i=1:m.num_var
            if m.var_type[i] == :Cont
                MOI.set(nlp_optimizer, MOI.VariablePrimalStart(), x[i], restart_values[i])
            else
                # discrete are fixed anyway
                MOI.set(nlp_optimizer, MOI.VariablePrimalStart(), x[i], nothing)
            end
        end
    end

    MOI.optimize!(nlp_optimizer)
    status = MOI.get(nlp_optimizer, MOI.TerminationStatus())

    obj_val = NaN
    real_sol = fill(NaN,m.num_var)
    if state_is_optimal(status; allow_almost=m.options.allow_almost_solved)
        obj_val = MOI.get(nlp_optimizer, MOI.ObjectiveValue())
        real_sol = MOI.get(nlp_optimizer, MOI.VariablePrimal(), x)
    end

    return status, real_sol, obj_val
end

"""
    add!(t::TabuList, m, sol)

Add a solution to the tabu list (includes only the discrete variables)
"""
function add!(t::TabuList, m, sol)
    t.sols[t.pointer] = [sol[m.disc2var_idx[i]] for i=1:m.num_disc_var]
    t.pointer += 1
    if t.pointer > t.length
        t.pointer = 1
    end
end

function get_fp_table(mip_obj,nlp_obj,t, fields, field_chars, catol)
    ln = ""
    i = 1
    arr = []
    for f in fields
        val = ""
        if f == :MIPobj
            if isnan(mip_obj)
                val = "-"
            else
                digits = 4
                if mip_obj < 1e-4
                    digits = convert(Int,floor(log10(1/catol)))
                end
                val = string(round(mip_obj; digits=digits))
            end
        elseif f == :NLPobj
            if isnan(nlp_obj)
                val = "-"
            else
                digits = 4
                if nlp_obj < 1e-4
                    digits = convert(Int,floor(log10(1/catol)))
                end
                val = string(round(nlp_obj; digits=digits))
            end
        elseif f == :Time
            val = string(round(t; digits=1))
        end

        if length(val) > field_chars[i]
            # too long to display shouldn't happen normally but is better than error
            # if it happens
            val = "t.l."
        end

        padding = field_chars[i]-length(val)
        ln *= repeat(" ",trunc(Int, floor(padding/2)))
        ln *= val
        ln *= repeat(" ",trunc(Int, ceil(padding/2)))
        push!(arr,val)
        i += 1
    end
    return ln, arr
end

"""
    fpump(optimizer, m)

Run the feasibility pump
"""
function fpump(optimizer, m)

    Random.seed!(JUNIPER_RNG, m.options.seed)

    if are_type_correct(m.relaxation_solution, m.var_type, m.disc2var_idx, m.options.atol)
        return Incumbent(m.relaxation_objval, m.relaxation_solution, only_almost_solved(m.status))
    end

    start_fpump = time()
    nlp_sol = m.relaxation_solution
    nlp_obj = 1 # should be not 0 for while
    c = 1
    tabu_list = TabuList()
    mip_sols = Dict{UInt64,Bool}()
    tabu_list.length = m.options.tabu_list_length
    tabu_list.pointer = 1
    tabu_list.sols = []
    for i=1:tabu_list.length
        push!(tabu_list.sols, NaN*ones(m.num_disc_var))
    end

    last_table_arr = []
    fields = []
    field_chars = []
    ps = m.options.log_levels
    # Print table init
    if check_print(ps,[:Table])
        fields, field_chars = [:MIPobj,:NLPobj,:Time], [20,20,5]
        print_table_header(fields,field_chars)
    end

    real_status = MOI.OPTIMIZE_NOT_CALLED
    fix = false
    nlp_status = :Error
    iscorrect = false
    tl = m.options.feasibility_pump_time_limit
    # the tolerance can be changed => current atol
    catol = m.options.atol
    atol_counter = 0
    while !are_type_correct(nlp_sol, m.var_type, m.disc2var_idx, catol) && time()-start_fpump < tl &&
        time()-m.start_time < m.options.time_limit

        # generate a mip or just round if no linear constraints
        if any(FS -> FS[1] != MOI.SingleVariable, MOI.get(LinearFilter(optimizer), MOI.ListOfConstraints()))
            mip_status, mip_sol, mip_obj = generate_mip(optimizer, m, nlp_sol, tabu_list, start_fpump)
        else
            # if no linear constraints just round the discrete variables
            mip_obj = NaN
            mip_sol = copy(nlp_sol)
            mip_status = MOI.OPTIMAL
            for vi=1:m.num_disc_var
                vidx = m.disc2var_idx[vi]
                mip_sol[vidx] = round(mip_sol[vidx])
            end
        end
        if mip_status != MOI.OPTIMAL
            @warn "MIP couldn't be solved to optimality. Terminated with status: "*string(mip_status)
            break
        end

        # If a cycle is detected which wasn't able to prevent by the tabu list (maybe too short)
        if haskey(mip_sols, hash(mip_sol))
            @warn "Cycle detected"
            break
        end
        add!(tabu_list, m, mip_sol)
        mip_sols[hash(mip_sol)] = true

        nlp_status, nlp_sol, nlp_obj = generate_nlp(optimizer, m, mip_sol, start_fpump)
        if !state_is_optimal(nlp_status; allow_almost=m.options.allow_almost_solved)
            cnlpinf = 0
            while cnlpinf < m.options.num_resolve_nlp_feasibility_pump && !state_is_optimal(nlp_status; allow_almost=m.options.allow_almost_solved) &&
                time()-start_fpump < tl && time()-m.start_time < m.options.time_limit
                nlp_status, nlp_sol, nlp_obj = generate_nlp(optimizer, m, mip_sol, start_fpump; random_start=true)
                cnlpinf += 1
            end
            if !state_is_optimal(nlp_status; allow_almost=m.options.allow_almost_solved)
                @warn "NLP couldn't be solved to optimality"
                if check_print(ps,[:Table])
                    print_fp_table(mip_obj, NaN, time()-start_fpump, fields, field_chars, catol)
                end
                break
            end
        end

        if check_print(ps,[:Table])
            print_fp_table(mip_obj, nlp_obj, time()-start_fpump, fields, field_chars, catol)
        end

        # if the current tolerance was nearly reached 5 times
        # => If reasonable should be an option
        if atol_counter >= m.options.feasibility_pump_tolerance_counter
            catol *= 10
            @warn "FPump tolerance changed to: ",catol
            atol_counter = 0
        end

        # if the difference is near 0 => try to improve the obj by using the original obj
        # set atol for type correct to a low value as it is checked with real_nlp anyway
        if are_type_correct(nlp_sol, m.var_type, m.disc2var_idx, catol*1000) || isapprox(nlp_obj, 0.0; atol=catol)
            real_status,real_sol, real_obj = generate_real_nlp(optimizer, m, mip_sol)
            cnlpinf = 0
            while cnlpinf < m.options.num_resolve_nlp_feasibility_pump && !state_is_optimal(real_status; allow_almost=m.options.allow_almost_solved) &&
                time()-start_fpump < tl && time()-m.start_time < m.options.time_limit
                real_status,real_sol, real_obj = generate_real_nlp(optimizer, m, mip_sol; random_start=true)
                cnlpinf += 1
            end
            if state_is_optimal(real_status) || (only_almost_solved(real_status) && m.options.allow_almost_solved_integral)
                if only_almost_solved(real_status)
                    @warn "Integral feasible point only almost solved. Disable with `allow_almost_solved_integral=false`"
                end
                nlp_obj = real_obj
                nlp_sol = real_sol
                iscorrect = true
                break
            elseif are_type_correct(nlp_sol, m.var_type, m.disc2var_idx, catol)
                real_status = MOI.LOCALLY_SOLVED
                nlp_obj = evaluate_objective(optimizer, m, nlp_sol)
                iscorrect = true
                @warn "Real objective wasn't solved to optimality"
                break
            end
        end
        if !isapprox(nlp_obj, 0.0; atol=catol) && isapprox(nlp_obj, 0.0; atol=10*catol)
            atol_counter += 1
        else
            atol_counter = 0
        end
        c += 1
    end

    if check_print(ps,[:Table])
        println()
    end

    if check_print(ps,[:Info])
        println("FP: ", time()-start_fpump, " s")
        println("FP: ", c == 1 ? "$c round" : "$c rounds")
    end
    m.fpump_info = Dict{Symbol,Float64}()
    m.fpump_info[:time] = time()-start_fpump
    m.fpump_info[:rounds] = c

    if iscorrect
        check_print(ps,[:Info]) && println("FP: Obj: ", nlp_obj)
        m.fpump_info[:obj] = nlp_obj
        m.fpump_info[:gap] = abs(m.relaxation_objval-nlp_obj)/abs(nlp_obj)
        return Incumbent(nlp_obj, nlp_sol, only_almost_solved(real_status))
    end

    m.fpump_info[:obj] = NaN
    m.fpump_info[:gap] = NaN
    check_print(ps,[:Info]) && println("FP: No integral solution found")

    return nothing
end


text:
#####
##### `AbstractMsgPackType`
#####

"""
    AbstractMsgPackType

An abstract type whose subtypes define a MessagePack <--> Julia type interface.

The subtypes of `AbstractMsgPackType` are:

- [`IntegerType`](@ref)
- [`NilType`](@ref)
- [`BooleanType`](@ref)
- [`FloatType`](@ref)
- [`StringType`](@ref)
- [`BinaryType`](@ref)
- [`ArrayType`](@ref)
- [`MapType`](@ref)
- [`ExtensionType`](@ref)
- [`AnyType`](@ref)
- [`ImmutableStructType`](@ref)
- [`MutableStructType`](@ref)
"""
abstract type AbstractMsgPackType end

"""
    IntegerType <: AbstractMsgPackType

A Julia type corresponding to the MessagePack Integer type.

If `msgpack_type(T)` is defined to return `IntegerType()`, then `T` must support:

- `to_msgpack(::IntegerType, ::T)::S`
- `from_msgpack(::Type{T}, ::S)::T`
- standard numeric comparators (`>`, `<`, `==`, etc.) against values of type `S`

where `S` may be one of the following types:

- `UInt8`
- `UInt16`
- `UInt32`
- `UInt64`
- `Int8`
- `Int16`
- `Int32`
- `Int64`
"""
struct IntegerType <: AbstractMsgPackType end

"""
    NilType <: AbstractMsgPackType

A Julia type corresponding to the MessagePack Nil type.

If `msgpack_type(T)` is defined to return `NilType()`, then `T` must support:

- `from_msgpack(::Type{T}, ::Nothing)::T`
"""
struct NilType <: AbstractMsgPackType end

"""
    BooleanType <: AbstractMsgPackType

A Julia type corresponding to the MessagePack Boolean type.

If `msgpack_type(T)` is defined to return `BooleanType()`, then `T` must support:

- `to_msgpack(::BooleanType, ::T)::Bool`
- `from_msgpack(::Type{T}, ::Bool)::T`
"""
struct BooleanType <: AbstractMsgPackType end

"""
    FloatType <: AbstractMsgPackType

A Julia type corresponding to the MessagePack Float type.

If `msgpack_type(T)` is defined to return `FloatType()`, then `T` must support:

- `to_msgpack(::FloatType, ::T)::S`
- `from_msgpack(::Type{T}, ::S)::T`
- standard numeric comparators (`>`, `<`, `==`, etc.) against values of type `S`

where `S` may be one of the following types:

- `Float32`
- `Float64`
"""
struct FloatType <: AbstractMsgPackType end

"""
    StringType <: AbstractMsgPackType

A Julia type corresponding to the MessagePack String type.

If `msgpack_type(T)` is defined to return `StringType()`, then `T` must support:

- `to_msgpack(::StringType, ::T)::String`
- `from_msgpack(::Type{T}, ::String)::T`
"""
struct StringType <: AbstractMsgPackType end

"""
    BinaryType <: AbstractMsgPackType

A Julia type corresponding to the MessagePack Binary type.

If `msgpack_type(T)` is defined to return `BinaryType()`, then `T` must support:

- `to_msgpack(::BinaryType, ::T)::Vector{UInt8}`
- `from_msgpack(::Type{T}, ::Vector{UInt8})::T`
"""
struct BinaryType <: AbstractMsgPackType end

"""
    ArrayType <: AbstractMsgPackType

A Julia type corresponding to the MessagePack Array type.

If `msgpack_type(T)` is defined to return `ArrayType()`, then `T` must support
the Julia `AbstractArray` interface, and/or must support:

- `to_msgpack(::ArrayType, ::T)::AbstractArray`
- `from_msgpack(::Type{T}, ::Vector)::T`
"""
struct ArrayType <: AbstractMsgPackType end

"""
    MapType <: AbstractMsgPackType

A Julia type corresponding to the MessagePack Map type.

If `msgpack_type(T)` is defined to return `MapType()`, then `T` must support
the Julia `AbstractDict` interface, and/or must support:

- `to_msgpack(::ArrayType, ::T)::AbstractDict`
- `from_msgpack(::Type{T}, ::Dict)::T`
"""
struct MapType <: AbstractMsgPackType end

"""
    ExtensionType <: AbstractMsgPackType

A Julia type corresponding to the MessagePack Extension type.

If `msgpack_type(T)` is defined to return `ExtensionType()`, then `T` must support:

- `to_msgpack(::ExtensionType, ::T)::Extension`
- `from_msgpack(::Type{T}, ::Extension)::T`

See also: [`Extension`](@ref)
"""
struct ExtensionType <: AbstractMsgPackType end

"""
    AnyType <: AbstractMsgPackType

The fallback return type of `msgpack_type(::Type)`, indicating that the given
Julia type does not have a known corresponding MessagePack type.
"""
struct AnyType <: AbstractMsgPackType end

"""
    ImmutableStructType <: AbstractMsgPackType

If `msgpack_type(T)` is defined to return `ImmutableStructType`, `T` will be
(de)serialized as a MessagePack Map type assuming certain constraints that
enable additional optimizations:

- `T` supports `fieldcount`, `fieldtype`, `fieldname`, `getfield`, and a constructor
that can be called as `T((getfield(x::T, i) for i in 1:fieldcount(T))...)` for any
`T` instance `x`.

- [`unpack`](@ref) will assume that incoming bytes to be deserialized to `T`
will always be formmatted as a MessagePack Map whose fields correspond exactly
to the fields of `T`. In other words, the `i`th key in the Map must correspond
to `fieldname(T, i)`, and the `i`th value must correspond to `getfield(::T, i)`.

This type is similar to [`MutableStructType`](@ref), but generally achieves
greater (de)serialization performance by imposing tighter constraints.
"""
struct ImmutableStructType <: AbstractMsgPackType end

"""
    MutableStructType <: AbstractMsgPackType

If `msgpack_type(T)` is defined to return `MutableStructType`, `T` will be
(de)serialized as a MessagePack Map type assuming certain constraints that
enable additional optimizations:

- `T` supports `fieldcount`, `fieldtype`, `fieldname`, `getfield`, `setfield!`,
and has the inner constructor `T() = new()`.

- [`unpack`](@ref) will assume that incoming bytes to be deserialized to `T`
will always be formmatted as a MessagePack Map whose fields are an unordered
subset of the fields of `T`. If a given field is not present in the MessagePack
Map, the corresponding field of the returned `T` instance will be left
uninitialized.

This type is similar to [`ImmutableStructType`](@ref), but imposes fewer
constraints at the cost of (de)serialization performance.
"""
struct MutableStructType <: AbstractMsgPackType end

#####
##### `msgpack_type`, `to_msgpack`, `from_msgpack` defaults
#####

"""
    msgpack_type(::Type{T}) where {T}

Return an instance of the [`AbstractMsgPackType`](@ref) subtype corresponding
to `T`'s intended MessagePack representation. For example:

```
msgpack_type(::Type{UUID}) = StringType()
```

If this method is overloaded such that `msgpack_type(T) === M()`, then `to_msgpack(::M, ::T)`
and `from_msgpack(::Type{T}, x)` should also be overloaded to handle conversion of `T`
instances to/from MsgPack2-compatible types.

By default, this method returns `AnyType()`. While this fallback method need not be overloaded to
support deserialization of `T` instances via `unpack`, `msgpack_type(T)` must be overloaded to
return a non-`AnyType` `AbstractMsgPackType` instance in order to support serializing `T` instances
via `pack`.

See also: [`from_msgpack`](@ref), [`to_msgpack`](@ref)
"""
msgpack_type(::Type) = AnyType()

"""
    to_msgpack(::M, value_to_serialize::T) where {M<:AbstractMsgPackType,T}

Return an "`M`-compatible" representation of `value_to_serialize` (for compatibility
definitions, see the docstrings for subtypes of [`AbstractMsgPackType`](@ref)).

By default, `to_msgpack` simply returns `value_to_serialize` directly.

The implementation of [`pack`](@ref) utilizes this function for every value
encountered during serialization, calling it in a manner similar to the
following psuedocode:

```
t = msgpack_type(T)
value_in_compatible_representation = to_msgpack(t, value_to_serialize::T)
_serialize_in_msgpack_format(t, value_in_compatible_representation)
```

For example, if `msgpack_type(UUID)` was defined to return `StringType()`, an
appropriate `to_msgpack` implementation might be:

```
to_msgpack(::StringType, uuid::UUID) = string(uuid)
```

See also: [`from_msgpack`](@ref), [`msgpack_type`](@ref), [`AbstractMsgPackType`](@ref)
"""
to_msgpack(::AbstractMsgPackType, x) = x

"""
    from_msgpack(::Type{T}, value_deserialized_by_msgpack) where {T}

Return the `value_deserialized_by_msgpack` converted to type `T`. By default,
this method simply calls `convert(T, value_deserialized_by_msgpack)`.

The implementation of [`unpack`](@ref) calls this function on every deserialized
value; in this case, `T` is generally derived from the type specified by
the caller of `unpack`.

For example, if `msgpack_type(UUID)` was defined to return `StringType()`, an
appropriate `from_msgpack` implementation might be:

```
from_msgpack(::Type{UUID}, uuid::AbstractString) = UUID(uuid)
```

See also: [`to_msgpack`](@ref), [`msgpack_type`](@ref), [`AbstractMsgPackType`](@ref)
"""
from_msgpack(T::Type, x) = convert(T, x)

# int-y things

msgpack_type(::Type{<:Integer}) = IntegerType()

# nil-y things

msgpack_type(::Type{Nothing}) = NilType()
msgpack_type(::Type{Missing}) = NilType()

from_msgpack(::Type{Missing}, ::Nothing) = missing

# bool-y things

msgpack_type(::Type{Bool}) = BooleanType()

# float-y things

msgpack_type(::Type{<:AbstractFloat}) = FloatType()

# string-y things

msgpack_type(::Type{<:AbstractString}) = StringType()
msgpack_type(::Type{Symbol}) = StringType()
msgpack_type(::Type{Char}) = StringType()

to_msgpack(::StringType, x::Char) = string(x)
to_msgpack(::StringType, x::Symbol) = _symbol_to_string(x)

from_msgpack(::Type{Char}, x::AbstractString) = first(x)
from_msgpack(::Type{Symbol}, x::AbstractString) = Symbol(x)

# array-y things

msgpack_type(::Type{<:AbstractArray}) = ArrayType()
msgpack_type(::Type{<:AbstractSet}) = ArrayType()
msgpack_type(::Type{<:Tuple}) = ArrayType()

from_msgpack(::Type{T}, x::Vector) where {T<:AbstractSet} = convert(T, Set(x))
from_msgpack(::Type{T}, x::Vector) where {T<:Tuple} = convert(T, (x...,))

# map-y things

msgpack_type(::Type{<:AbstractDict}) = MapType()
msgpack_type(::Type{<:NamedTuple}) = MapType()

@generated function to_msgpack(::MapType, x::NamedTuple)
    fields = Any[]
    for i in 1:fieldcount(x)
        push!(fields, Expr(:tuple, Expr(:quote, fieldname(x, i)), :(getfield(x, $i))))
    end
    return Expr(:tuple, fields...)
end

@generated function from_msgpack(::Type{T}, x::Dict) where {T<:NamedTuple}
    fields = Any[]
    for i in 1:fieldcount(T)
        name = fieldname(T, i)
        strname = String(name)
        push!(fields, :($name = x[$strname]))
    end
    return Expr(:tuple, fields...)
end

#####
##### `Extension`
#####

"""
    struct Extension
        type::Int8
        data::Vector{UInt8}
    end

A wrapper for bytes formatted to the MessagePack Extension type specification.

See also: [`extserialize`](@ref), [`extdeserialize`](@ref)
"""
struct Extension
    type::Int8
    data::Vector{UInt8}
end

Base.:(==)(a::Extension, b::Extension) = a.type == b.type && a.data == b.data

msgpack_type(::Type{Extension}) = ExtensionType()

@deprecate Ext(type, data) MsgPack2.Extension(type, data)

"""
    extserialize(type, x)

Return `Extension(type, data)` where `data` is the result of calling Julia's
built-in `serialize` function on `x`, facilitating direct serialization Julia
values to MessagePack format:

```
julia> struct Point{T}
           x::T
           y::T
       end

julia> val = [Point(rand(), rand()) for _ in 1:100];

julia> bytes = MsgPack2.pack(MsgPack2.extserialize(123, x));

julia> type, new_val = MsgPack2.extdeserialize(MsgPack2.unpack(bytes));

julia> type == 123
true

julia> new_val == val
true
```

Note that `extserialize`/[`extdeserialize`](@ref) are subject to the same caveat
as Julia's built-in `serialize`/`deserialize` functions: "In general, this
process will not work if the reading and writing are done by different versions
of Julia, or an instance of Julia with a different system image.".

See also: [`Extension`](@ref), [`extdeserialize`](@ref)
"""
function extserialize(type, x)
    io = IOBuffer()
    serialize(io, x)
    return Extension(type, take!(io))
end

"""
    extdeserialize(x::Extension)

Return `(x.type, value)` where `value` is the result of calling Julia's built-in
`deserialize` function on `x.data`.

This function is meant to be used in conjuction with [`extserialize`](@ref); see
that function's docstring for more details.

See also: [`Extension`](@ref), [`extserialize`](@ref)
"""
extdeserialize(x::Extension) = (x.type, deserialize(IOBuffer(x.data)))

#####
##### `Skip`
#####

struct Skip{T} end

msgpack_type(::Type{Skip{T}}) where {T} = msgpack_type(T)

from_msgpack(::Type{T}, x) where {T<:Skip} = T()

#####
##### `PointerString`
#####

struct PointerString
    ptr::Ptr{UInt8}
    len::UInt64
end

Base.length(s::PointerString) = s.len

Base.write(io::IO, s::PointerString) = Base.unsafe_write(io, s.ptr, s.len)

from_msgpack(::Type{Symbol}, x::PointerString) = ccall(:jl_symbol_n, Ref{Symbol}, (Ptr{UInt8}, Int), x.ptr, x.len)
from_msgpack(::Type{PointerString}, x::PointerString) = x
from_msgpack(S::Type, x::PointerString) = from_msgpack(S, unsafe_string(x.ptr, x.len))

_symbol_to_string(x::Symbol) = PointerString(Base.unsafe_convert(Ptr{UInt8}, x), sizeof(x))


text:
using MISAlgorithms
using Test

@testset "mis1" begin
    include("mis1.jl")
end

@testset "mis2" begin
    include("mis2.jl")
end


text:
using Oceananigans.Fields: FieldSlicer

#####
##### JLD2OutputWriter tests
#####

function jld2_field_output(model)

    model.clock.iteration = 0
    model.clock.time = 0.0

    set!(model, u = (x, y, z) -> rand(),
                v = (x, y, z) -> rand(),
                w = (x, y, z) -> rand(),
                T = 0,
                S = 0)

    simulation = Simulation(model, Δt=1.0, stop_iteration=1)

    simulation.output_writers[:velocities] = JLD2OutputWriter(model, model.velocities,
                                                                    schedule = TimeInterval(1),
                                                                        dir = ".",
                                                                     prefix = "test",
                                                                      force = true)

    u₀ = CUDA.@allowscalar data(model.velocities.u)[3, 3, 3]
    v₀ = CUDA.@allowscalar data(model.velocities.v)[3, 3, 3]
    w₀ = CUDA.@allowscalar data(model.velocities.w)[3, 3, 3]

    run!(simulation)

    file = jldopen("test.jld2")

    # Data is saved without halos by default
    u₁ = file["timeseries/u/0"][3, 3, 3]
    v₁ = file["timeseries/v/0"][3, 3, 3]
    w₁ = file["timeseries/w/0"][3, 3, 3]

    close(file)

    rm("test.jld2")

    FT = typeof(u₁)

    return FT(u₀) == u₁ && FT(v₀) == v₁ && FT(w₀) == w₁
end

function jld2_sliced_field_output(model, outputs=model.velocities)

    model.clock.iteration = 0
    model.clock.time = 0.0

    set!(model, u = (x, y, z) -> rand(),
                v = (x, y, z) -> rand(),
                w = (x, y, z) -> rand())

    simulation = Simulation(model, Δt=1.0, stop_iteration=1)

    simulation.output_writers[:velocities] =
        JLD2OutputWriter(model, outputs,
                             schedule = TimeInterval(1),
                         field_slicer = FieldSlicer(i=1:2, j=1:2:4, k=:),
                                  dir = ".",
                               prefix = "test",
                                force = true)

    run!(simulation)

    file = jldopen("test.jld2")

    u₁ = file["timeseries/u/0"]
    v₁ = file["timeseries/v/0"]
    w₁ = file["timeseries/w/0"]

    close(file)

    rm("test.jld2")

    return size(u₁) == (2, 2, 4) && size(v₁) == (2, 2, 4) && size(w₁) == (2, 2, 5)
end

function run_jld2_file_splitting_tests(arch)
    model = NonhydrostaticModel(architecture=arch, grid=RegularRectilinearGrid(size=(16, 16, 16), extent=(1, 1, 1)))
    simulation = Simulation(model, Δt=1, stop_iteration=10)

    function fake_bc_init(file, model)
        file["boundary_conditions/fake"] = π
    end

    ow = JLD2OutputWriter(model, (u=model.velocities.u,); dir=".", prefix="test", schedule=IterationInterval(1),
                          init=fake_bc_init, including=[:grid],
                          field_slicer=nothing, array_type=Array{Float64},
                          max_filesize=200KiB, force=true)

    push!(simulation.output_writers, ow)

    # 531 KiB of output will be written which should get split into 3 files.
    run!(simulation)

    # Test that files has been split according to size as expected.
    @test filesize("test_part1.jld2") > 200KiB
    @test filesize("test_part2.jld2") > 200KiB
    @test filesize("test_part3.jld2") < 200KiB

    for n in string.(1:3)
        filename = "test_part$n.jld2"
        jldopen(filename, "r") do file
            # Test to make sure all files contain structs from `including`.
            @test file["grid/Nx"] == 16

            # Test to make sure all files contain info from `init` function.
            @test file["boundary_conditions/fake"] == π
        end

        # Leave test directory clean.
        rm(filename)
    end
end

function jld2_time_averaging_of_horizontal_averages(model)

    model.clock.iteration = 0
    model.clock.time = 0.0

    set!(model, u = (x, y, z) -> 1,
                v = (x, y, z) -> 2,
                w = (x, y, z) -> 0,
                T = (x, y, z) -> 4)

    u, v, w = model.velocities
    T, S = model.tracers

    simulation = Simulation(model, Δt=1.0, stop_iteration=5)

    u, v, w = model.velocities
    T, S = model.tracers

    average_fluxes = (wu = AveragedField(w * u, dims=(1, 2)),
                      uv = AveragedField(u * v, dims=(1, 2)),
                      wT = AveragedField(w * T, dims=(1, 2)))

    simulation.output_writers[:fluxes] = JLD2OutputWriter(model, average_fluxes,
                                                          schedule = AveragedTimeInterval(4, window=2),
                                                               dir = ".",
                                                            prefix = "test",
                                                             force = true)

    run!(simulation)

    file = jldopen("test.jld2")

    # Data is saved without halos by default
    wu = file["timeseries/wu/4"][1, 1, 3]
    uv = file["timeseries/uv/4"][1, 1, 3]
    wT = file["timeseries/wT/4"][1, 1, 3]

    close(file)

    rm("test.jld2")

    FT = eltype(model.grid)

    return wu == zero(FT) && wT == zero(FT) && uv == FT(2)
end

for arch in archs
    # Some tests can reuse this same grid and model.
    topo =(Periodic, Periodic, Bounded)
    grid = RegularRectilinearGrid(topology=topo, size=(4, 4, 4), extent=(1, 1, 1))
    model = NonhydrostaticModel(architecture=arch, grid=grid)

    @testset "JLD2 output writer [$(typeof(arch))]" begin
        @info "  Testing JLD2 output writer [$(typeof(arch))]..."

        @test jld2_field_output(model)
        @test jld2_sliced_field_output(model)
        @test jld2_sliced_field_output(model, (u = model -> model.velocities.u,
                                               v = model -> model.velocities.v,
                                               w = model -> model.velocities.w))

        run_jld2_file_splitting_tests(arch)

        @test jld2_time_averaging_of_horizontal_averages(model)
    end
end


text:
### A Pluto.jl notebook ###
# v0.12.21

using Markdown
using InteractiveUtils

# ╔═╡ 09988650-82cc-11eb-33d6-f581c3a7f16d
using FileIO,Images,ImageCore,ImageFeatures,ImageView,ImageDraw,ColorVectorSpace,ImageMorphology,ImageBinarization

# ╔═╡ fed405b6-828d-11eb-2d86-977a78c747ee
md"""
# Correct Nonuniform Illumination and Analyze Foreground Objects
"""

# ╔═╡ c5e83452-8283-11eb-21e5-4d3d301bf70a
begin
img = load("assets/rice.png"); #second way
img_gray = Gray.(img)
end

# ╔═╡ b6c1ef02-8302-11eb-05df-3d28ed8ddf5d
opening(img_gray)

# ╔═╡ baaa035e-8302-11eb-1db8-1be228d4ddfa

# ╔═╡ ba8e5dd4-8302-11eb-073e-f38cb0d923d8

# ╔═╡ ba399810-8302-11eb-3971-23a86794a1f9

# ╔═╡ Cell order:
# ╠═fed405b6-828d-11eb-2d86-977a78c747ee
# ╠═09988650-82cc-11eb-33d6-f581c3a7f16d
# ╠═c5e83452-8283-11eb-21e5-4d3d301bf70a
# ╠═b6c1ef02-8302-11eb-05df-3d28ed8ddf5d
# ╠═baaa035e-8302-11eb-1db8-1be228d4ddfa
# ╠═ba8e5dd4-8302-11eb-073e-f38cb0d923d8
# ╠═ba399810-8302-11eb-3971-23a86794a1f9

max_stars_repo_name: ashwani-rathee/Demos

max_issues_repo_name: ashwani-rathee/Demos

max_forks_repo_name: ashwani-rathee/Demos


text:
"""
Create the base routes for this webservice.
The basic interface is ```/experiments```.
Issuing a POST request to ```/experiments```, will create a new experiment and report back the experiment id.

For accessing previously created experiments, append the experiment id to the URL.
I.e. to access experiment1, aim your requests at ```/experiments/experiment1```.
A simple GET request to this URL with give back a full JSON representation of this experiment (carefull, this can be large).

Access to the different parts of data in an experiment is given by again appending to the URL.
Currently, the different data types that can be accessed or worked with are:
 - genes
 - proteins
 - cell_lines
 - drugs
 - pathways
 - outcomes

POST request to the according URLs will create the data inside of an experiment.
I.e. to import gene representations into experiment1, issue a POST request to ```/experiments/experiment1/genes``` and provide the
JSON representation of a list of genes.

"""
function create_base_router()
    r = HTTP.Handlers.Router()
    HTTP.@register(r, "/", HTTP.Handlers.HandlerFunction(handle_base_request))
	# load sample data
	HTTP.@register(r, "/load_iorio_data", load_iorio_data)
	HTTP.@register(r, "/load_dream_challenge_data", load_dream_challenge_data)
    # experiment methods
    # create experiments
    HTTP.@register(r, "POST", "/experiments", handle_create_experiment_request)

    # query for whole experiment
    HTTP.@register(r, "GET", "/experiments/*", handle_get_experiment_request)

    # delete experiment
    HTTP.@register(r, "DELETE", "/experiments/*", handle_delete_experiment_request)

    # data methods
    # create collections, can be genes, proteins, cell_lines, drugs, pathways and outcomes
    HTTP.@register(r, "POST", "/experiments/*/data/*", handle_create_collection_request)

    # query for collections
    HTTP.@register(r, "GET", "/experiments/*/data/*", handle_get_collection_request)

    # query single objects
    HTTP.@register(r, "GET", "/experiments/*/data/*/*", handle_get_object_request)

    # delete collections
    HTTP.@register(r, "DELETE", "/experiments/*/data/*", handle_delete_collection_request)

    # delete single objects
    HTTP.@register(r, "DELETE", "/experiments/*/data/*/*", handle_delete_object_request)

    # computation methods
    HTTP.@register(r, "POST", "/experiments/*/train", train_request)
    HTTP.@register(r, "GET", "/experiments/*/progress", progress_request)
    HTTP.@register(r, "GET", "/experiments/*/results", results_request)
    HTTP.@register(r, "POST", "/experiments/*/predict", predict_request)
	HTTP.@register(r, "POST", "/debug", debug_request)
    r
end


text:
# This code is part of QuantumCircuits.
#
# (C) Copyright Rafał Pracht 2022.
#
# This code is licensed under the Apache License, Version 2.0. You may
# obtain a copy of this license in the LICENSE.txt file in the root directory
# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
#
# Any modifications or derivative works of this code must retain this
# copyright notice, and modified files need to carry a notice indicating
# that they have been altered from the originals.

using Test

using QuantumCircuits
using QuantumCircuits.QML

using QuantumCircuits.QCircuits.Math
using QuantumCircuits.QCircuits.Circuit: toQiskit
using QuantumCircuits.QCircuits.Gates: Rx, Ry, Rz, X, Y, Z, S, Sd, T, Td, H, CX, U3

using Zygote

function test_derivate(qc, explen)
    pa = getRandParameters(qc)
    @test length(pa) == explen
    loss(params) = real(matrix_norm(tomatrix(qc, params)))
    loss(pa)
    loss'(pa)
end

################################################################################
#  1 qubit, rotation gate                                                      #
################################################################################
gates = [Rx, Ry, Rz]
θs = [0, π/4, π/2, 3π/4, π, 5π/4, 3π/2, 7π/4, 2π]

for g in gates
    for θ in θs
        circ = QCircuit(1)
        add!(circ, g, 0, ParameterT(θ))
        err = unitary_error(tomatrix(toQiskit(circ)), tomatrix(circ))
        if err >= 1e-8
            println(circ)
        end
        @test err < 1e-8
        test_derivate(circ, 1)
    end
end

################################################################################
#  1 qubit, 1 gate unitary                                                     #
################################################################################
gates = [X, Y, Z, S, Sd, T, Td, H]

for g in gates
    circ = QCircuit(1)
    add!(circ, g, 0)
    @test unitary_error(tomatrix(toQiskit(circ)), tomatrix(circ)) < 1e-8
    test_derivate(circ, 0)
end

################################################################################
#  1 qubit, 2/3 gate unitary                                                   #
################################################################################
for g in gates
    for g2 in gates
        circ = QCircuit(1)
        add!(circ, g, 0)
        add!(circ, g2, 0)
        @test unitary_error(tomatrix(toQiskit(circ)), tomatrix(circ)) < 1e-8
        test_derivate(circ, 0)
    end
end
for g in gates
    for g2 in gates
        for g3 in gates
            circ = QCircuit(1)
            add!(circ, g, 0)
            add!(circ, g2, 0)
            add!(circ, g3, 0)
            @test unitary_error(tomatrix(toQiskit(circ)), tomatrix(circ)) < 1e-8
            test_derivate(circ, 0)
        end
    end
end

qc = QCircuit(1)
qc.u3(0, ParameterT(0), ParameterT(1.0), ParameterT(2.0))
qc.x(0)
qc.u3(0, ParameterT(2.0), ParameterT(3.0), ParameterT(4.0))
qiskitmatrix = tomatrix(toQiskit(qc))
@test unitary_error(qiskitmatrix, tomatrix(qc)) < 1e-8
params = [4.0, 3.0, 2.0, 1.0, 0.0, -1.0]
gcparammatrix = tomatrix(qc, params)
@test unitary_error(qiskitmatrix, gcparammatrix) > 1.0
setparameters!(qc, params)
qiskitmatrix = tomatrix(toQiskit(qc))
@test unitary_error(qiskitmatrix, gcparammatrix) < 1e-8

################################################################################
#  2 qubit, 2 gate unitary                                                     #
################################################################################
for g in gates
    for g2 in gates
        circ = QCircuit(2)
        add!(circ, g, 0)
        add!(circ, g2, 1)
        @test unitary_error(tomatrix(toQiskit(circ)), tomatrix(circ)) < 1e-8
        test_derivate(circ, 0)
    end
end

for g in gates
    for g2 in gates
        for g3 in gates
            circ = QCircuit(2)
            add!(circ, g, 0)
            add!(circ, g2, 1)
            add!(circ, g3, 1)
            @test unitary_error(tomatrix(toQiskit(circ)), tomatrix(circ)) < 1e-8
            test_derivate(circ, 0)
        end
    end
end

for g in gates
    for g2 in gates
        for g3 in gates
            circ = QCircuit(3)
            add!(circ, g, 0)
            add!(circ, g2, 1)
            add!(circ, g3, 2)
            @test unitary_error(tomatrix(toQiskit(circ)), tomatrix(circ)) < 1e-8
            test_derivate(circ, 0)
        end
    end
end

################################################################################
#  CX                                                                          #
################################################################################
function test_unitary(qubits, gates, expected_param=0)
    circ = QCircuit(qubits)
    for (g, args) in gates
        add!(circ, g, args...)
    end
    @test unitary_error(tomatrix(toQiskit(circ)), tomatrix(circ)) < 1e-8
    test_derivate(circ, expected_param)
end

test_unitary(2, [(CX, [0, 1])])
test_unitary(2, [(CX, [1, 0])])
test_unitary(2, [(H, [0]), (H, [1]), (CX, [0, 1]), (X, [0]), (Z, [1])])
test_unitary(2, [(H, [0]), (H, [1]), (CX, [1, 0]), (X, [0]), (Z, [1])])
test_unitary(6, [(H, [0]), (H, [1]), (Z, [2]), (CX, [2, 3]), (S, [0]), (Y, [1]),
                 (Td, [5])])
test_unitary(6, [(H, [0]), (H, [1]), (CX, [3, 2]), (X, [0]), (Sd, [1]), (H, [4])])

test_unitary(3, [(CX, [1, 2])])

test_unitary(4, [(CX, [0, 3])])
test_unitary(4, [(CX, [0, 3]), (H, [0]), (Y, [1])])
test_unitary(4, [(CX, [3, 0])])
test_unitary(5, [(CX, [1, 3])])
test_unitary(5, [(CX, [3, 1])])

test_unitary(6, [(H, [0]), (H, [1]), (CX, [2, 3]), (S, [0]), (Y, [1]),
                 (Td, [5]), (Z, [2]), (CX, [0, 3]), (H, [1]),
                 (Z, [2]), (S, [5]), (CX, [0, 3]), (Td, [3]), (H, [4]),
                 (Y, [1]), (Y, [2]), (CX, [1, 3]), (H, [5]), (H, [3]),
                 (CX, [3, 1]), (Z, [2])])

################################################################################
#  U3                                                                          #
################################################################################
test_unitary(1, [(U3, (0, ParameterT(π/4), ParameterT(π/2), ParameterT(π)))], 3)
test_unitary(1, [(U3, (0, ParameterT(-π/4), ParameterT(-π/2), ParameterT(-π)))], 3)
test_unitary(1, [(U3, (0, ParameterT(π/4), ParameterT(π/2), ParameterT(π))), (U3, (0, ParameterT(-π/4), ParameterT(-π/2), ParameterT(-π)))], 6)
test_unitary(2, [(U3, (0, ParameterT(π/4), ParameterT(π/2), ParameterT(π))), (CX, [1, 0]), (U3, (1, ParameterT(-π/4), ParameterT(-π/2), ParameterT(-π)))], 6)


text:
using NNlib
using Base.Test

@testset "NNlib" begin

include("activation.jl")
include("conv.jl")

xs = rand(5)
@test softmax(xs) ≈ exp.(xs) ./ sum(exp.(xs))
@test logsoftmax(xs) ≈ log.(softmax(xs))
@test logsigmoid.(xs) ≈ log.(sigmoid.(xs))

xs = rand(5,10)
@test softmax(xs) ≈ exp.(xs) ./ sum(exp.(xs),1)
@test logsoftmax(xs) ≈ log.(softmax(xs))
@test logsigmoid.(xs) ≈ log.(sigmoid.(xs))

for T in [:Float32, :Float64]
  @eval @test logsigmoid.($T[-100_000, 100_000.]) ≈ $T[-100_000, 0.]
end

## compare the outputs with the PyTorch nn.LogSoftmax returns
xs = Float32[1, 2, 3000.]
@test logsoftmax(xs) ≈ [-2999, -2998, 0]

xs = Float32[1 2 3; 1000 2000 3000]
@test logsoftmax(xs) ≈ [-999 -1998 -2997; 0 0 0.]
end


text:
# Model
include_model("quadruped")

include(joinpath(pwd(), "models/visualize.jl"))
vis = Visualizer()
render(vis)

# Horizon
T = 5

# Time step
h = 0.1
tf = h * (T - 1)

function initial_configuration(model::Quadruped, θ1, θ2, θ3)
    q1 = zeros(model.nq)
    q1[3] = pi / 2.0
    q1[4] = -θ1
    q1[5] = θ2

    q1[8] = -θ1
    q1[9] = θ2

    q1[2] = model.l_thigh1 * cos(q1[4]) + model.l_calf1 * cos(q1[5])

    q1[10] = -θ3
    q1[11] = acos((q1[2] - model.l_thigh2 * cos(q1[10])) / model.l_calf2)

    q1[6] = -θ3
    q1[7] = acos((q1[2] - model.l_thigh2 * cos(q1[6])) / model.l_calf2)

    return q1
end

θ1 = pi / 4.0
θ2 = pi / 4.0
θ3 = pi / 3.0

q1 = initial_configuration(model, θ1, θ2, θ3)
visualize!(vis, model, [q1])

# control
# u1 = initial_torque(model, q1, h)[model.idx_u]
# ul <= u <= uu
_uu = Inf * ones(model.m)
_uu[model.idx_u] .= 0.0
_ul = zeros(model.m)
_ul[model.idx_u] .= -0.0
ul, uu = control_bounds(model, T, _ul, _uu)

xl, xu = state_bounds(model, T,
    x1 = [q1; q1])

# Objective
include_objective(["velocity", "nonlinear_stage"])
q_ref = linear_interpolation(q1, q1, T)
render(vis)
visualize!(vis, model, q_ref)
x0 = configuration_to_state(q_ref)

# penalty on slack variable
obj_penalty = PenaltyObjective(1.0e4, model.m)

# quadratic tracking objective
# Σ (x - xref)' Q (x - x_ref) + (u - u_ref)' R (u - u_ref)
obj_control = quadratic_time_tracking_objective(
    [1.0 * Diagonal(1.0e-5 * ones(model.n)) for t = 1:T],
    [1.0 * Diagonal([1.0e-5 * ones(model.nu)..., zeros(model.m - model.nu)...]) for t = 1:T-1],
    [[q1; q1] for t = 1:T],
    [zeros(model.m) for t = 1:T],
    1.0)

obj = MultiObjective([obj_penalty, obj_control])
# Constraints
include_constraints(["stage", "contact", "free_time", "loop"])

con_contact = contact_constraints(model, T)
con = multiple_constraints([con_contact])#,
    # con_free_time])#, con_loop,
	# con_pinned1, con_pinned2])

# cc = zeros(4)
# xx = rand(model.n)
# uu = rand(model.m)

# pinned2!(cc, xx, uu)

# Problem
prob = trajectory_optimization_problem(model,
               obj,
               T,
               xl = xl,
               xu = xu,
               ul = ul,
               uu = uu,
			   h = h,
               con = con)

# trajectory initialization
u0 = [1.0e-2 * rand(model.m) for t = 1:T-1] # random controls

# Pack trajectories into vector
z0 = pack(x0, u0, prob)
# z0 .+= 0.001 * randn(prob.num_var)

# Solve
include_snopt()
# @load joinpath(@__DIR__, "quadruped_gait.jld2") z̄ q̄ ū τ̄ λ̄ b̄ h̄
@time z̄, info = solve(prob, copy(z0),
    nlp = :SNOPT7,
    tol = 1.0e-3, c_tol = 1.0e-3,
	max_iter = 5000,
    time_limit = 60 * 2, mapl = 5)

check_slack(z̄, prob)
x̄, ū = unpack(z̄, prob)
q̄ = state_to_configuration(x̄)
τ̄ = [u[model.idx_u] for u in ū]
λ̄ = [u[model.idx_λ] for u in ū]
b̄ = [u[model.idx_b] for u in ū]

[norm(fd(model, x̄[t+1], x̄[t], ū[t], zeros(model.d), h, t)) for t = 1:T-1]
# Visualize
vis = Visualizer()
render(vis)
# open(vis)
visualize!(vis, model, state_to_configuration(x̄), Δt = h)


text:
module LyceumDevTools

using Base: UUID
using Base64

using GitHub: GitHub, GitHubAPI, Repo, AnonymousAuth, Authorization
using GitHub: authenticate_headers!, authenticate, github2json, api_uri

using HTTP: HTTP
using JuliaFormatter: format_file, format
using LibGit2: LibGit2

using Pkg
using Pkg.Types: projectfile_path, read_project, write_project
using Pkg.Types: manifestfile_path, read_manifest, write_manifest

using PkgTemplates: Template, generate
using Registrator: Registrator
using RegistryTools: register

include("configs.jl")

export treehash, repository_dispatch
include("git.jl")

export ly_format_file, ly_format
export with_tempdir, genkeys, with_sandbox_env, cpinto, flattendir
include("misc.jl")

export ly_generate, ly_register, incversion!
include("packaging.jl")

include("Compat/Compat.jl")
using .Compat

include("TestUtil/TestUtil.jl")
using .TestUtil

end # module


text:
"""
$SIGNATURES

Find footnotes refs and defs and eliminate the ones that don't verify the appropriate regex.
For a footnote ref: `\\[\\^[a-zA-Z0-0]+\\]` and `\\[\\^[a-zA-Z0-0]+\\]:` for the def.
"""
function validate_footnotes!(tokens::Vector{Token})
    fn_refs = Vector{Token}()
    rm      = Int[]
    for (i, τ) in enumerate(tokens)
        if τ.name == :FOOTNOTE_REF
            # footnote ref [^1]:
            m = match(r"^\[\^[a-zA-Z0-9]+\](:)?$", τ.ss)
            if m !== nothing
                if m.captures[1] !== nothing
                    # it's a def
                    tokens[i] = Token(:FOOTNOTE_DEF, τ.ss)
                else
                    # it's a ref, take and delete
                    push!(fn_refs, τ)
                    push!(rm, i)
                end
            else
                # delete
                push!(rm, i)
            end
        end
    end
    deleteat!(tokens, rm)
    return fn_refs
end

"""
$SIGNATURES

Verify that a given string corresponds to a well formed html entity.
"""
function validate_html_entity(ss::AS)
    match(r"&(?:[a-z0-9]+|#[0-9]{1,6}|#x[0-9a-f]{1,6});", ss) !== nothing
end

"""
$(SIGNATURES)

Given a candidate header block, check that the opening `#` is at the start of a line, otherwise
ignore the block.
"""
function validate_header_block(β::OCBlock)::Bool
    # skip non-header blocks
    β.name ∈ MD_HEADER || return true
    # if it's a header block, have a look at the opening token
    τ = otok(β)
    # check if it overlaps with the first character
    from(τ) == 1 && return true
    # otherwise check if the previous character is a linereturn
    s = str(β.ss) # does not allocate
    prevc = s[prevind(str(β.ss), from(τ))]
    prevc == '\n' && return true
    return false
end

"""
$(SIGNATURES)

Keep track of link defs.
"""
function validate_and_store_link_defs!(blocks::Vector{OCBlock})::Nothing
    isempty(blocks) && return
    rm = Int[]
    parent = str(blocks[1])
    for (i, β) in enumerate(blocks)
        if β.name == :LINK_DEF
            # incremental backward look until we find a `[` or a `\n` if `\n` first, discard
            ini  = prevind(parent, from(β))
            k    = ini
            char = '\n'
            while k ≥ 1
                char = parent[k]
                char ∈ ('[','\n') && break
                k = prevind(parent, k)
            end
            if char == '['
                # redefine the full block
                ftk = Token(:FOO,subs(""))
                # we have a [id]: lk add it to PAGE_LINK_DEFS
                id = subs(parent, nextind(parent, k), ini) |> htmlesc
                lk = β |> content |> strip |> string
                PAGE_LINK_DEFS[id] = lk
                # replace the block by a full one so that it can be fully
                # discarded in the process of md blocks
                blocks[i] = OCBlock(:LINK_DEF, ftk=>ftk, subs(parent, k, to(β)))
            else
                # discard
                push!(rm, i)
            end
        end
    end
    deleteat!(blocks, rm)
    return nothing
end


text:
"""
    UncertainIndexDataset

Generic dataset containing uncertain indices.

## Fields
- **`indices::AbstractVector{AbstractUncertainValue}`**: The uncertain values.
"""
struct UncertainIndexDataset <: AbstractUncertainIndexDataset
    indices::AbstractVector{<:AbstractUncertainValue}
end

"""
    ConstrainedUncertainIndexDataset

Constrained dataset containing uncertain indices.

## Fields
- **`indices::AbstractVector{AbstractUncertainValue}`**: The uncertain indices.
"""
struct ConstrainedUncertainIndexDataset <: AbstractUncertainIndexDataset
    indices::AbstractVector{<:AbstractUncertainValue}
end

function UncertainIndexDataset(x::AbstractArray{T, 1}) where T
    UncertainIndexDataset(CertainValue.(x))
end

export
UncertainIndexDataset,
ConstrainedUncertainIndexDataset,
distributions


text:
using JuliaGrid
using HDF5
using SparseArrays
using Test

include("flow/flow.jl")
include("generator/generator.jl")
include("estimation/estimation.jl")

max_stars_repo_name: mcosovic/JuliaGrid

max_issues_repo_name: mcosovic/JuliaGrid

max_forks_repo_name: mcosovic/JuliaGrid


text:
using Rubi
using Test

@testset "Rubi.jl" begin
    # Write your tests here.
end


text:
export @default_compute_rf_inner_solver
# Try to automatically determine an appropriate scalar compute_rf solver.
macro default_compute_rf_inner_solver(nep)
    return (nep isa PEP) ? PolyeigInnerSolver() : ScalarNewtonInnerSolver();
end

"""
    struct ScalarNewtonInnerSolver <: InnerSolver
    function ScalarNewtonInnerSolver(;tol=eps()*100,maxit=80,bad_solution_allowed=true)

This type can be used with  `compute_rf`. It solves the scalar nonlinear problem
using a Rayleig functional. In contrast to the `NewtonInnerSolver`, this
function does not require an SPMF (or precomputations associated with SPMFs).
"""
struct ScalarNewtonInnerSolver <: InnerSolver
    tol::Float64
    maxit::Int
    bad_solution_allowed::Bool
    function ScalarNewtonInnerSolver(;tol=eps()*100,maxit=80,bad_solution_allowed=true)
        return new(tol,maxit,bad_solution_allowed);
    end
end

function compute_rf(T0::Type{T}, nep::NEP, x, inner_solver::ScalarNewtonInnerSolver;
                        y=x, target=zero(T), λ0=target,
                        kwargs...) where T
    λ_iter = T(λ0);
    Δλ = T(Inf)
    count = 0
    while (abs(Δλ)>inner_solver.tol) & (count<inner_solver.maxit)
        count = count+1
        # compute function value and derivative
        z1 = compute_Mlincomb(nep, λ_iter, reshape(x,size(nep,1),1))
        z2 = compute_Mlincomb(nep, λ_iter, reshape(x,size(nep,1),1),[T(1)],1)

        Δλ = -dot(y,z1)/dot(y,z2);
        λ_iter += Δλ
    end

    if ((count==inner_solver.maxit) && (! (inner_solver.bad_solution_allowed )))
        throw(NoConvergenceException());
    end

    # Return type is a vector of correct type
    λ_star::Array{T,1} = Array{T,1}(undef, 1)
    if (T <: Real) && (typeof(λ_iter) != T) && (imag(λ_iter)/real(λ_iter) < inner_solver.tol)
        # Looking for a real quantity (AND) iterate is not real (AND) complex part is negligible
        λ_star[1] = real(λ_iter) # Truncate to real
    else
        λ_star[1] = λ_iter
    end
    return λ_star
end

"""
    compute_rf(eltype::Type,nep::NEP,x,inner_solver::InnerSolver;
               y=x, target=zero(T), λ0=target,TOL=eps(real(T))*1e3,max_iter=10)

Computes the Rayleigh functional of the `nep`, i.e., computes a vector
``Λ`` of values ``λ``
such that ``y^TM(λ)x=0``, using the procedure specified in `inner_solver`.
The default behaviour consists of a scalar valued
Newton-iteration, and the returned vector has only one element.

The given eltype<:Number is the type of the returned vector.

# Example

This uses `iar` to solve the (scalar) nonlinear problem.

```julia-repl
julia> nep=nep_gallery("dep0");
julia> x=ones(size(nep,1));
julia> s=compute_rf(ComplexF64,nep,x,IARInnerSolver())[1] # Take just first element
-1.186623627962043 - 1.5085094961223182im
julia> x'*compute_Mlincomb(nep,s,x)
-8.881784197001252e-16 + 1.0880185641326534e-14im
```
"""
    function compute_rf(T0::Type{T}, nep::NEP, x, inner_solver::InnerSolver;
                        y=x, target=zero(T), λ0=target,
                        TOL=eps(real(T))*1e3, max_iter=30,kwargs...) where T

        pnep=create_proj_NEP(nep);
        n=size(nep,1);
        set_projectmatrices!(pnep,reshape(y,n,1),reshape(x,n,1));

        local λv
        try
            (λv,xv)=inner_solve(inner_solver,T0,pnep,σ=target,λv=[λ0]);
        catch e
            # Even return the approximation upon failure
            if (e isa NoConvergenceException)
                λv=e.λ
            else
                rethrow(e)
            end
        end
        if (λv isa Vector)
            return λv
        else
            return [λv];
        end
    end

function compute_rf(::Type{T},nep::NEP,x, inner_solver::PolyeigInnerSolver; y=x, target=zero(T), λ0=target,
                    TOL=eps(real(T))*1e3,max_iter=10) where T

    @show λ0
    a=zeros(T,size(nep.A,1))
    for i=1:size(nep.A,1)
        a[i]=dot(y,nep.A[i]*x);
    end
    m=size(nep.A,1);
    # Build a companion matrix (maybe bigfloat)
    A=zeros(T,m-1,m-1);
    for k=1:m-1
        if (k<m-1)
            A[k,k+1]=1;
        end
        A[end,k]=-a[k]/a[end]
    end
    pp=eigvals(A);
    if (T <: Real) # If we specify real arithmetic, return real eigvals
        pp=real(pp);
    end

    II=sortperm(abs.(pp .- target))
    return pp[II];
end

# Overload a version of eigvals(::Matrix{BigFloat}) for compute_rf
import LinearAlgebra.eigvals
function eigvals(A::Union{Matrix{BigFloat},Matrix{Complex{BigFloat}}})
    T=eltype(A);
    Tfloat= (T <: Complex) ? (ComplexF64) : (Float64)
    Afloat=Matrix{Tfloat}(A);
    λv,X=eigen(Afloat);
    local λv_bigfloat=Vector{Complex{BigFloat}}(λv);
    # Some Rayleigh quotient iteration steps in bigfloat to improve accuracy:
    for j=1:size(λv,1)
        local s=λv_bigfloat[j];
        z=Vector{Complex{BigFloat}}(X[:,j])
        for f=1:3
            z=(A-s*I)\z; normalize!(z);
            s=z'*A*z;
        end
        λv_bigfloat[j]=s
    end
    if (maximum(abs.((λv_bigfloat-λv)./λv_bigfloat)) > eps()*100)
        @warn("Correction iteration of eigvals for bigfloat, made the eigenvalues move considerably")
    end

    return λv_bigfloat
end


text:
export Clock, TICK, SLEEP_NSEC, SLEEP_SEC
export tick, job!, oneshot!

"""
```julia
mutable struct Clock
  start_event::Base.Event
  is_stopped::Bool
  frequency::AbstractFloat
  Clock() = new(Base.Event(), true, 0.01667)
end
```
# Fields
- `start_event`::Base.Event Whether awake! has been called
- `is_stopped`::Bool Whether shutdown! has been called
- `frequency`::AbstractFloat Tick frequency
"""
mutable struct Clock
  start_event::Base.Event
  is_stopped::Bool
  frequency::AbstractFloat
  Clock() = new(Base.Event(), true, 0.01667)
end

"""
```julia
struct TICK
  Δ::AbstractFloat # seconds
end
```
"""
struct TICK
  Δ::AbstractFloat # seconds
end

"""
```julia
struct SLEEP_SEC
  Δ::AbstractFloat
end
```
"""
struct SLEEP_SEC
  Δ::AbstractFloat
end

"""
```julia
struct SLEEP_NSEC
  Δ::UInt
end
```
"""
struct SLEEP_NSEC
  Δ::UInt
end

"""
```julia
function Base.sleep(s::SLEEP_SEC) ::Float64
  t1 = time()
  while true
    if time() - t1 >= s.Δ break end
    yield()
  end
  return time() - t1
end
```

`sleep(::SLEEP_SEC) -> Float64`

Sleep the specified number of seconds.
"""
function Base.sleep(s::SLEEP_SEC) ::Float64
  t1 = time()
  while true
    if time() - t1 >= s.Δ break end
    yield()
  end
  return time() - t1
end

"""
```julia
function Base.sleep(s::SLEEP_NSEC)
  t1 = time_ns()
  while true
    if time_ns() - t1 >= s.Δ break end
    yield()
  end
  return time_ns() - t1
end
```

`sleep(::SLEEP_NSEC) -> UInt64`

Sleep the specified number of nanoseconds
"""
function Base.sleep(s::SLEEP_NSEC) ::UInt64
  t1 = time_ns()
  while true
    if time_ns() - t1 >= s.Δ break end
    yield()
  end
  return time_ns() - t1
end

"""
```julia
function tick(Δ)
  δ = sleep(SLEEP_NSEC(Δ * 1e9))
  sendMessage(TICK(δ / 1e9))
end
```

`tick(::Any) -> Nothing`

Sleep Δ seconds, the raise a TICK event with the actual amount of time slept.
Called in a background task in an infinite loop. Primary purpose is to trigger subsystem TICK handlers.
"""
function tick(Δ) ::Nothing
  δ = sleep(SLEEP_NSEC(Δ * 1e9))
  sendMessage(TICK(δ / 1e9))
  Log.@debug "tick"

  return nothing
end

"""
```julia
function job!(c::Clock, f, arg=1)
  function job()
    Base.wait(c.start_event)
    while !c.is_stopped
      f(arg)
    end
  end
  schedule(Task(job))
end
```

`job!(::Clock, f::Any, arg::Any) -> Task`

Schedule a background task to be run and synchronized with clock state.
Used internally to run dispatchMessage and tick.
"""
function job!(c::Clock, f, arg=1) ::Task
  function job()
    Base.wait(c.start_event)
    while !c.is_stopped
      f(arg)
    end
  end
  schedule(Task(job))
end

"""
```julia
function oneshot!(c::Clock, f, arg=1)
  function oneshot()
    Base.wait(c.start_event)
    f(arg)
  end
  schedule(Task(oneshot))
end
```

`oneshot!(::Clock, f::Any, arg::Any) -> Task`

Schedule a background task to be run once, synchronized with Clock state.
"""
function oneshot!(c::Clock, f, arg=1) ::Task
  function oneshot()
    Base.wait(c.start_event)
    f(arg)
  end
  schedule(Task(oneshot))
end

"""
```julia
function awake!(c::Clock)
  job!(c, tick, c.frequency)

  c.is_stopped = false

  Base.notify(c.start_event)
end
```

`awake!(::Clock) -> Nothing`

Starts the tick job and signals all waiting Tasks that the clock has started.
"""
function awake!(c::Clock) ::Nothing
  Log.@debug "Clock awake!"

  job!(c, tick, c.frequency)
  c.is_stopped = false
  Base.notify(c.start_event)
end

"""
```julia
function shutdown!(c::Clock)
  c.is_stopped = true

  c.start_event = Base.Event()
end
```

`shutdown!(::Clock) -> Nothing`

Signal Tasks that the Clock has stopped.
"""
function shutdown!(c::Clock) ::Nothing
  Log.@debug "Clock shutdown!"

  c.is_stopped = true
  c.start_event = Base.Event() # old one remains signaled no matter what, replace
  return nothing
end


text:
"""
    plot_cyclotron(p::MagneticParticle; use_cell=true, kwargs...)
Plot the circle traced by the free particle motion. Optionally use `p.current_cell` for
the particle's position. The user provided `kwargs...`
are passed onto `PyPlot.plt.Circle`.
"""
function plot_cyclotron(p::MagneticParticle; use_cell=true, kwargs...)
  ω = p.omega
  pos = use_cell ? p.pos + p.current_cell : p.pos
  c = pos - (1/ω)*[p.vel[2], -p.vel[1]]
  r = abs(1/ω)
  circle1 = PyPlot.plt.Circle(c, r,
  edgecolor = (0.0,0.0,0.8, 0.5), linewidth = 2.0, facecolor = (0., 0.0, 0.8, 0.05),
  kwargs...)
  PyPlot.gca().add_artist(circle1)
end

function plot_particle(x,y,vx,vy; ax = PyPlot.gca(), c = "k", kwargs...)
    # Plot position:
    s1 = ax.scatter(x, y; color=c, s = 20.0, kwargs..., zorder = 99)
    # Plot velocity:
    q1 = ax.quiver(x, y, 0.04vx, 0.04vy; angles = "xy", scale = 1, width = 0.004, color=c, zorder = 99)

    return s1, q1
end

"""
    plot(p::AbstractParticle [, cyclotron=false]; use_cell=true, ax=gca(), kwargs...)
Plot given particle on `ax`. Optionally use `p.current_cell` for
the particle's position. Given `kwargs...` are passed onto `PyPlot.scatter`.

The particle is represented as a small ball (`PyPlot.scatter`) and a small arrow
(`PyPlot.quiver`).
All `kwargs...` are given to `scatter` but if a keyword argument `color` is given,
it is also passed to `quiver`.

Optionally you can plot the cyclotron traced by a `MagneticParticle` by giving
`true` as second argument.
"""
function plot(p::AbstractParticle{T}, cycl::Bool = false; use_cell=true, kwargs...) where {T}
    if typeof(p) <: MagneticParticle && cycl
        plot_cyclotron(p; use_cell = use_cell)
    end
    pos = use_cell ? p.pos + p.current_cell : p.pos
    kwargs = Dict(kwargs)
    # Set same color for arrow and point:
    if haskey(kwargs, :color)
        c = kwargs[:color]
    else
        c = (0,0,0)
    end

    return plot_particle(pos..., p.vel...; c = c, kwargs...)
end


text:
function process_reload_hash(request::HTTP.Request, state::HandlerState)
    reload_tuple = (
        reloadHash = state.reload.hash,
        hard = state.reload.hard,
        packages = keys(state.cache.resources.files),
        files = state.reload.changed_assets
    )
    state.reload.hard = false
    state.reload.changed_assets = []
    return HTTP.Response(200, ["Content-Type" => "application/json"], body = JSON3.write(reload_tuple))

end


text:
# This file was generated, do not modify it. # hide
print(
    "Optimal hyper-parameters: \n",
    "  max_depth: ", best_booster.max_depth, "\n",
    "  η:         ", best_booster.η
)

max_stars_repo_name: alan-turing-institute/MLJTutorials

max_issues_repo_name: alan-turing-institute/MLJTutorials

max_forks_repo_name: alan-turing-institute/MLJTutorials


text:
function menfromboys(arr)
    menfirst = (a,b) -> ( iseven(a) && !iseven(b) ) ||
                        ( iseven(a) == iseven(b) == (a < b) )
    result = sort(arr,lt=menfirst)
    unique!(result)
end

max_stars_repo_name: UlrichBerntien/Codewars-Katas

max_issues_repo_name: UlrichBerntien/Codewars-Katas

max_forks_repo_name: UlrichBerntien/Codewars-Katas


text:
#
# Factors Indexing Tests
#

@testset "Factors Index" begin

@testset "basics" begin
    x = Dimension(:X, [2, 5, 7])
    y = Dimension(:Y, [2, 3])
    z = Dimension(:Z, ['a', 'c'])

    ft = Factor([x, y, z], [1, 2, 3, 2, 3, 4, 4, 6, 7, 8, 10, 16])
    df = DataFrame(ft)

    a = Assignment(:Y=> [3, 2], :K => 16, :Z => 'a')
    @test ft[a].potential == [2 1; 3 2; 4 3]

    ft[Assignment(:X=>2, :Y=> 2, :Z=>'c')] = 1600.0
    @test ft.potential[1, 1, 2] == 1600.0
    @test DataFrame(ft)[sub2ind(ft, 1, 1, 2), :potential] == 1600.0
end

@testset "basics 2" begin
    ft = Factor([1, 2, 3, 2, 3, 4, 4, 6, 7, 8, 10, 16], :X=>3, :Y =>2, :Z=>2)
    df = DataFrame(ft)

    @test_throws MethodError ft[:Y => "waldo"]
    @test_throws ArgumentError ft[:Y => 16]

    a = Assignment(:Y=> 2, :K => 16, :Z => 1)
    @test ft[a].potential == [2, 3, 4]

    ft[:X => 2, :Y => 1, :Z => 2] = 1600.0
    @test ft.potential[2, 1, 2] == 1600.0
    @test DataFrame(ft)[sub2ind(ft, 2, 1, 2), :potential] == 1600.0

    ft[:X => 1, :Y => 2] = 2016
    @test ft.potential[1, 2, :] == [2016, 2016]
end

@testset "dims updating" begin
    ft = Factor([1 4; 2 5; 3 6], :X => 2:4, :Y=>['a', 'b'])
    ft2 = ft[:X=>[3, 2], :Y=>'a']
    @test names(ft2) == [:X]
    @test first(scope(ft2)) == Dimension(:X, [3, 2])
end

@testset "bits" begin
    ft = Factor([1 4; 2 5; 3 6], :X => 2:4, :Y=>['a', 'b'])
    @test_throws BoundsError ft[:X=>[true, false]]

    ft2 = ft[:X=>[true, false, true]]
    @test names(ft2) == [:X, :Y]
    @test getdim(ft2, :X) == Dimension(:X, [2, 4])
    @test ft2.potential == [1 4; 3 6]
end

end


text:
using ExampleServer
using Test

@testset "ExampleServer.jl" begin
    # Write your tests here.
end


text:
"""
    BoseFS2C{NA,NB,M,AA,AB} <: AbstractFockAddress
    BoseFS2C(onr_a, onr_b)

Address type that constructed with two [`BoseFS{N,M,S}`](@ref). It represents a
Fock state with two components, e.g. two different species of bosons with particle
number `NA` from species S and particle number `NB` from species B. The number of
modes `M` is expected to be the same for both components.
"""
struct BoseFS2C{NA,NB,M,SA,SB,N} <: AbstractFockAddress{N,M}
    bsa::BoseFS{NA,M,SA}
    bsb::BoseFS{NB,M,SB}

end

function BoseFS2C(bsa::BoseFS{NA,M,SA}, bsb::BoseFS{NB,M,SB}) where {NA,NB,M,SA,SB}
    N = NA + NB
    return BoseFS2C{NA,NB,M,SA,SB,N}(bsa, bsb)
end
BoseFS2C(onr_a::Tuple, onr_b::Tuple) = BoseFS2C(BoseFS(onr_a),BoseFS(onr_b))

function Base.show(io::IO, b::BoseFS2C)
    print(io, "BoseFS2C(")
    Base.show(io,b.bsa)
    print(io, ",")
    Base.show(io,b.bsb)
    print(io, ")")
end

num_components(::Type{<:BoseFS2C}) = 2
Base.isless(a::T, b::T) where {T<:BoseFS2C} = isless((a.bsa, a.bsb), (b.bsa, b.bsb))
onr(b2::BoseFS2C) = (onr(b2.bsa), onr(b2.bsb))

function near_uniform(::Type{<:BoseFS2C{NA,NB,M}}) where {NA,NB,M}
    return BoseFS2C(near_uniform(BoseFS{NA,M}), near_uniform(BoseFS{NB,M}))
end

"""
    CompositeFS(addresses::SingleComponentFockAddress...) <: AbstractFockAddress

Used to encode addresses for multi-component models. All component addresses
are expected have the same number of modes.

See also: [`BoseFS`](@ref), [`FermiFS`](@ref), [`SingleComponentFockAddress`](@ref),
[`num_modes`](@ref), [`FermiFS2C`](@ref).
"""
struct CompositeFS{C,N,M,T} <: AbstractFockAddress{N,M}
    components::T
    # C: components, N: total particles, M: modes in each component, T: tuple type with constituent address types
    function CompositeFS{C,N,M,T}(adds::T) where {C,N,M,T}
        return new{C,N,M,T}(adds)
    end
end

# Slow constructor - not to be used internallly
function CompositeFS(adds::Vararg{AbstractFockAddress}) # CompositeFS(adds)
    N = sum(num_particles, adds)
    M1, M2 = extrema(num_modes, adds)
    if M1 ≠ M2
        error("all addresses must have the same number of modes")
    end
    return CompositeFS{length(adds),N,M1,typeof(adds)}(adds)
end

num_components(::CompositeFS{C}) where {C} = C
Base.hash(c::CompositeFS, u::UInt) = hash(c.components, u)

function Base.show(io::IO, c::CompositeFS{C}) where {C}
    println(io, "CompositeFS(")
    for add in c.components
        println(io, "  ", add, ",")
    end
    print(io, ")")
end

"""
    update_component(c::CompositeFS, new, ::Val{i})

Replace the `i`-th component in `c` with `new`. Used for updating a single component in the
address.
"""
function update_component(c::CompositeFS, new, ::Val{I}) where {I}
    return typeof(c)(_update_component(c.components, new, Val(I)))
end

@inline _update_component((a, as...), new, ::Val{1}) = (new, as...)
@inline function _update_component((a, as...), new, ::Val{I}) where {I}
    return (a, _update_component(as, new, Val(I - 1))...)
end

Base.isless(a::T, b::T) where {T<:CompositeFS} = isless(a.components, b.components)

function onr(a::CompositeFS)
    map(onr, a.components)
end

# Convenience
"""
    FermiFS2C <: AbstractFockAddress
    FermiFS2C(onr_a, onr_b)

Fock state address with two fermionic (spin) components. Alias for [`CompositeFS`](@ref)
with two [`FermiFS`](@ref) components.
Construct by specifying either two compatible [`FermiFS`](@ref)s or [`onr`](@ref)s.
```jldoctest
julia> f2c = FermiFS2C((1,0,0),(0,1,1))
CompositeFS(
  FermiFS{1,3}((1, 0, 0)),
  FermiFS{2,3}((0, 1, 1)),
)

julia> f2c isa FermiFS2C
true
```
"""
const FermiFS2C{N1,N2,M,N,F1,F2} =
    CompositeFS{2,N,M,Tuple{F1,F2}} where {F1<:FermiFS{N1,M},F2<:FermiFS{N2,M}}
FermiFS2C(f1::FermiFS{<:Any,M}, f2::FermiFS{<:Any,M}) where {M} = CompositeFS(f1, f2)
FermiFS2C(onr_a, onr_b) = FermiFS2C(FermiFS(onr_a), FermiFS(onr_b))


text:
### A Pluto.jl notebook ###
# v0.16.4

using Markdown
using InteractiveUtils

# ╔═╡ d80342ae-522b-11eb-19c0-cf931fb4ae81
begin
	import Pkg
	Pkg.activate(Base.current_project())
	using Statistics
end

# ╔═╡ 9c2ef32a-522d-11eb-29e6-29da29e4dec7
using Plots

# ╔═╡ 96096b41-fc25-43f9-ab75-50c5196e9604
md"# Chapter 8.5 Expected vs. Sample Updates"

# ╔═╡ 88b47ee6-522d-11eb-1a98-ed4dd43dfd11
begin
	Base.@kwdef mutable struct SampleAvg
		t::Int = 0
		avg::Float64 = 0.0
	end
	function (s::SampleAvg)(x)
		s.t += 1
		s.avg += (x - s.avg) / s.t
		s.avg
	end
end

# ╔═╡ 8fde026e-522d-11eb-2aa1-f3e28c238ea6
function run_once(b)
    rms = Float64[]
    distribution = randn(b)
    expectation = mean(distribution)
    sample_avg = SampleAvg()

    for i in 1:2*b
        avg = sample_avg(distribution[rand(1:b)])
        push!(rms, abs(avg - expectation))
    end
    rms
end

# ╔═╡ 39865f6a-8500-4277-884f-81b7de151fc6
md"## Figure 8.7"

# ╔═╡ 969b48aa-522d-11eb-2039-61916ddf87d4
begin
	n_runs = 1000
	p = plot(legend=:topright)

	for b in [2, 10, 100, 1000]
		rms = mean(run_once(b) for _ in 1:n_runs)
		xs = (1:2*b) ./ b
		plot!(p, xs, rms, label="b=$b", ylabel="RMS error in value estimate")
	end

	p
end

# ╔═╡ Cell order:
# ╟─96096b41-fc25-43f9-ab75-50c5196e9604
# ╠═d80342ae-522b-11eb-19c0-cf931fb4ae81
# ╠═88b47ee6-522d-11eb-1a98-ed4dd43dfd11
# ╠═8fde026e-522d-11eb-2aa1-f3e28c238ea6
# ╟─39865f6a-8500-4277-884f-81b7de151fc6
# ╠═9c2ef32a-522d-11eb-29e6-29da29e4dec7
# ╠═969b48aa-522d-11eb-2039-61916ddf87d4


text:
export selu

doc"""
    selu(x::Var)

Scaled Exponential Linear Unit.

```math
f(x) = \lambda
\begin{cases}
x & x > 0 \\
\alpha e^{x}-\alpha & x\leq0
\end{cases}
```
where ``\lambda=1.0507`` and ``\alpha=1.6733``.

# References
Klambauer et al., ["Self-Normalizing Neural Networks"](https://arxiv.org/abs/1706.02515), NIPS 2017.
"""
selu(x::Var) = Var(selu.(x.data), (selu,x))
selu(x::T) where T = x > 0 ? T(1.0507)*x : T(1.0507)*T(1.6733)*(exp(x)-1)
selu(x::Node) = Node(selu, x)

function addgrad!(y::Var, ::typeof(selu), x::Var)
    isvoid(x.grad) || ∇selu!(y.data, y.grad, x.data, x.grad)
end

function ∇selu!(y::Array{T}, gy::Array{T}, x::Array{T}, gx::Array{T}) where T
    alpha = T(1.6733)
    lambda = T(1.0507)
    @inbounds for i = 1:length(x)
        gx[i] +=  x[i] > 0 ? gy[i]*lambda : gy[i]*(y[i]+lambda*alpha)
    end
end


text:
function convert_units(v::Vector{Float64})
    Float64[convert_units(v[i], MAP_IND2VAR_L[i]) for i=1:endof(v)]
end

function unconvert_units(v::Vector{Float64})
    Float64[unconvert_units(v[i], MAP_IND2VAR_L[i]) for i=1:endof(v)]
end

function convert_units(x::Float64, var::Symbol)
  if var == :v_d0 || var == :v_d1
    return x * 1.68780
  elseif var == :h_d0 || var == :h_d1
    return x / 60
  else
    return x
  end
end

function unconvert_units(x::Float64,var::Symbol)
  if var == :v_d0 || var == :v_d1
    return x / 1.68780
  elseif var == :h_d0 || var == :h_d1
    return x * 60
  else
    return x
  end
end

function val2ind(boundariesi, ri, value)
  if !isempty(boundariesi)
    index = findfirst(x -> (x > value), boundariesi) - 1

    if index == -1
      index = ri
    end
  else
    index = value
  end
  return index
end

function discretize(p::CorrAEMParameters, v::Vector{Float64})
  Int64[val2ind(p.boundaries[MAP_L2G[i]], p.r_transition[MAP_L2G[i]], val) for (i, val) in enumerate(v)]
end

function dediscretize(dval::Int64, boundaries::Vector{Float64}, zero_bin::Int64)
  val_min = boundaries[dval]
  val_max = boundaries[dval + 1]

  if dval == zero_bin
    val = 0.0
    prob = 1.0
  elseif val_max == val_min
    val = val_min
    prob = 1.0
  else
    val = randfloat(val_min, val_max)
    prob = 1.0 / (val_max - val_min)
    #this is a density so it won't be normalized to [0,1]
  end

  return (val, prob)
end

function select_random_cumweights(cweights::Vector{Float64})
  r = cweights[end] * rand()

  return findfirst(x -> (x >= r), cweights)
end

function alt_bounds(L::Int64)
    if L == 1
        return (1000.0, 3000.0)
    elseif L == 2
        return (3000.0, 10000.0)
    elseif L == 3
        return (10000.0, 18000.0)
    elseif L == 4
        return (18000.0, 29000.0)
    elseif L == 5
        return (29000.0, 40000.0)
    else
        error("invalid L")
    end
end

#mods x to the range [-b, b]
function to_plusminus_b(x::AbstractFloat, b::AbstractFloat)
  z = mod(x, 2 * b)

  return (z > b) ? (z - 2 * b) : z
end

to_plusminus_180(x::AbstractFloat) = to_plusminus_b(x, 180.0)
randfloat(xmin::Float64, xmax::Float64) = xmin + rand() * (xmax-xmin)


text:
module FileSystemCache

using ...SimpleTasks.Types

import SimpleTasks.Services.Bucket
import SimpleTasks.Services.Cache

export FileSystemCacheService

const FOLDER_SEPARATOR = "/"

type FileSystemCacheService <: CacheService
    base_directory::AbstractString

    FileSystemCacheService(base_directory::AbstractString) =
        isempty(strip(base_directory)) ?  throw(ArgumentError("Base directory can
            not be empty")) : new(base_directory)
end

const MEGABYTE = 1024 * 1024
const WRITE_BYTES = 50 * MEGABYTE

function sanitize(key::AbstractString)
    return replace(key, r"\.\.", "")
end

function to_filename(cache::FileSystemCacheService, key::AbstractString)
    return "$(cache.base_directory)$FOLDER_SEPARATOR$(sanitize(key))"
end

function create_path(full_path_file_name::AbstractString)
    path_end = rsearch(full_path_file_name, "/").start
    if path_end > 0
        mkpath(full_path_file_name[1:path_end])
    end
end

function Cache.haskey(cache::FileSystemCacheService, key::AbstractString)
    filename = to_filename(cache, key)
    return isfile(filename) && isreadable(filename)
end

function Cache.put!(cache::FileSystemCacheService, key::AbstractString,
        value_io::IO)
    filename = to_filename(cache, key)
    create_path(filename)
    filestream = open(filename, "w")
    try
        if !iswritable(filestream)
            error("Unable to write to $filename")
        end

        if typeof(value_io) <: IOBuffer && position(value_io) == value_io.size
            println("wARNING: trying to read from an IOBuffer with current " *
                "position at the end of the buffer")
        end
        while !eof(value_io)
            write(filestream, readbytes(value_io, WRITE_BYTES))
        end
        close(filestream)
    catch e
        # if there was an error in reading/writing the stream, we should
        # clear the file in case of transmission error
        truncate(filestream, 0)
        println(STDERR, "Clearing FSCache for $key in $filename due to " *
            "stream error")
        showerror(STDERR, e, catch_backtrace(); backtrace = true)
        # error is thrown later on because the file is empty
    end

    # if the outcome of the put resulted in an empty file, then remove it
    # from the cache
    close(filestream)
    if Cache.haskey(cache, key) && stat(filename).size == 0
        rm(filename)
        error("Tried to cache an empty file $key into $filename")
    end
end

function Cache.get(cache::FileSystemCacheService, key::AbstractString)
    if Cache.haskey(cache, key)
        return open(to_filename(cache, key), "r")
    else
        throw(KeyError(key))
    end
end

function Cache.delete!(cache::FileSystemCacheService, key::AbstractString)
    if Cache.haskey(cache, key)
        rm(to_filename(cache, key))
    end
    return cache
end

function Cache.clear!(cache::FileSystemCacheService)
    contents = readdir(cache.base_directory)
    for content in contents
        rm(to_filename(cache, content); recursive=true)
    end
end

end # module FileSystemCache

max_forks_repo_name: seung-lab/SimpleTasks


text:
#finitesemi.jl
#Algorithm 5.6 of Random Eigenvalues by Alan Edelman

#Experiment:    Eigenvalues of GUE matrices
#Plot:          Histogram of eigenvalues
#Theory:        Semicircle and finite semicircle

## Parameters
n = 3           # size of matrix
s = 10000       # number of samples
d = 0.1         # bin size

## Experiment
function finitesemi_experiment(n,s,d)
    e = Float64[]
    for i = 1:s
        a = randn(n,n)+im*randn(n,n)
        a = (a+a')/(2*sqrt(4n))
        v = eigvals(a)
        append!(e, v)
    end
    return hist(e, -1.5:d:1.5)
end
grid, count = finitesemi_experiment(n,s,d)

## Theory
x = -1:0.01:1
y = sqrt(1-x.^2)

## Plot
using Winston
p = FramedPlot()
h = Histogram(count*pi/(2*d*n*s), step(grid))
h.x0 = first(grid)
add(p, h)
add(p, Curve(x, y, "color", "red", "linewidth", 2))
include("levels.jl")
add(p, Curve(levels(n)..., "color", "blue", "linewidth", 2))
if isinteractive()
    Winston.display(p)
else
    file(p, "spherecomponent.png")
end


text: my_f(x,y) = x+y

max_stars_repo_name: zhaotianjing/MyExample

max_issues_repo_name: zhaotianjing/MyExample

max_forks_repo_name: zhaotianjing/MyExample


text:
function _template(; n, N)
    dirs = Vector{SingleEntryVector{Float64}}()
    d = sum(n^i for i in 1:N)

    for i in 1:n
        x = SingleEntryVector(i, d, 1.0)
        push!(dirs, x)
    end
    for i in 1:n
        x = SingleEntryVector(i, d, -1.0)
        push!(dirs, x)
    end
    return CustomDirections(dirs)
end

# TODO check if it can be removed
function _project(sol, vars)
    πsol_1n = Flowpipe([ReachSet(set(project(R, vars)), tspan(R)) for R in sol])
end

# TODO refactor to MathematicalSystems
import MathematicalSystems: statedim

struct CanonicalQuadraticForm{T, MT<:AbstractMatrix{T}} <: AbstractContinuousSystem
    F1::MT
    F2::MT
end
statedim(f::CanonicalQuadraticForm) = size(f.F1, 1)
canonical_form(s::CanonicalQuadraticForm) = s.F1, s.F2
export CanonicalQuadraticForm
#-

# TODO generalize to AbstractContinuousSystem using vector_field
function canonical_form(s::BlackBoxContinuousSystem)
    @requires Symbolics
    f = s.f
    # differentiate
end

function post(alg::CARLIN, ivp::IVP{<:AbstractContinuousSystem}, tspan; Δt0=interval(0), kwargs...)
    @unpack N, compress, δ, bloat, resets = alg

    # for now we assume there are no resets
    if haskey(kwargs, :NSTEPS)
        NSTEPS = kwargs[:NSTEPS]
        T = NSTEPS * δ
    else
        # get time horizon from the time span imposing that it is of the form (0, T)
        T = _get_T(tspan, check_zero=true, check_positive=true)
        NSTEPS = ceil(Int, T / δ)
    end

    # extract initial states
    X0 = initial_state(ivp)
    # compute canonical quadratic form from the given system
    F1, F2 = canonical_form(system(ivp))
    n = statedim(ivp)
    dirs = _template(n=n, N=N)
    alg = LGG09(δ=δ, template=dirs, approx_model=Forward())

    return reach_CARLIN_alg(X0, F1, F2; alg, resets, N, T, Δt0, bloat, compress)
end


text:
using AbstractAlgebra

@testset "AbstractAlgebra io consistency" begin
    R, x = PolynomialRing(GF(2^31 - 1), "x")
    @test Groebner.groebner([x^2 - 4, x + 2]) == [x + 2]

    R, x = PolynomialRing(QQ, "x")
    @test Groebner.groebner([x^2 - 4, x + 2]) == [x + 2]

    R, x = PolynomialRing(ZZ, "x")
    @test Groebner.groebner([x^2 - 4, x + 2]) == [x + 2]
end


text:
meta(T) = (usetyped ? IRTools.typed_meta : IRTools.meta)(T)

function code_ir(f, T)
  m = meta(Tuple{Typeof(f),T.parameters...})
  return IRCode(m)
end

function code_irm(ex)
  isexpr(ex, :call) || error("@code_ir f(args...)")
  f, args = ex.args[1], ex.args[2:end]
  :(code_ir($(esc(f)), typesof($(esc.(args)...))))
end

macro code_ir(ex)
  code_irm(ex)
end

function argnames!(meta, names...)
  meta.code.slotnames = [names...]
  meta.code.slotflags = [0x00 for name in names]
end

function spliceargs!(meta, ir::IRCode, args...)
  for i = 1:length(ir.stmts)
    ir[SSAValue(i)] = argmap(x -> Argument(x.n+length(args)), ir[SSAValue(i)])
  end
  for (name, T) in reverse(args)
    pushfirst!(ir.argtypes, T)
    pushfirst!(meta.code.slotnames, name)
  end
  return ir
end

# Behave as if the function signature is f(args...)
function varargs!(meta, ir::IRCode, n = 1)
  isva = meta.method.isva
  Ts = widenconst.(ir.argtypes[n+1:end])
  argtypes = !isva ?
    Any[ir.argtypes[1:n]..., Tuple{Ts...}] :
    Any[ir.argtypes[1:n]..., Tuple{Ts[1:end-1]...,Ts[end].parameters...}]
  empty!(ir.argtypes); append!(ir.argtypes, argtypes)
  ir = IncrementalCompact(ir)
  map = Dict{Argument,Any}()
  for i = 1:(length(Ts)-isva)
    map[Argument(i+n)] = insert_node_here!(ir, xcall(Base, :getfield, Argument(n+1), i), Ts[i], Int32(0))
  end
  if isva
    i = length(Ts)
    xs, T = Argument(n+1), argtypes[end]
    for _ = 1:i-1
      T = Tuple{T.parameters[2:end]...}
      xs = insert_node_here!(ir, xcall(Base, :tail, xs), T, Int32(0))
    end
    map[Argument(i+n)] = xs
  end
  for (i, x) in ir
    ir[i] = argmap(a -> get(map, a, a), x)
  end
  return finish_dc(ir)
end

function pis!(ir::IRCode)
  for i = 1:length(ir.stmts)
    ex = ir.stmts[i]
    ex isa PiNode || continue
    ir.stmts[i] = xcall(Core, :typeassert, ex.val, ex.typ)
  end
  return ir
end

function slots!(ir::IRCode)
  n = 0
  for b = 1:length(ir.cfg.blocks)
    i = first(ir.cfg.blocks[b].stmts)
    while (phi = ir[SSAValue(i)]) isa PhiNode
      slot = IRTools.Slot(Symbol(:phi, n+=1))
      ir[SSAValue(i)] = slot
      for (pred, val) in zip(phi.edges, phi.values)
        insert_blockend!(ir, pred, Any, :($slot = $val))
      end
      i += 1
    end
  end
  return compact!(ir)
end

@generated function roundtrip(f, args...)
  m = meta(Tuple{f,args...})
  ir = IRCode(m)
  ir = varargs!(m, ir)
  argnames!(m, :f, :args)
  ir = spliceargs!(m, ir, (Symbol("#self#"), typeof(roundtrip)))
  ir = slots!(pis!(ir))
  return IRTools.update!(m, ir)
end

function inlineable!(ir)
  insert_node!(ir, 1, Any, Expr(:meta, :inline))
  compact!(ir)
end

function log!(ir, msg)
  insert_node!(ir, 1, Any, xcall(Core, :println, msg))
  compact!(ir)
end