text:
function moveSide(r::Robot, side::HorizonSide, needMark=false)
    # Move untill you find pr
    steps = 0
    while isborder(r, side) != true
        moveStep(r, side, needMark)
        steps += 1
    end
    return steps
end

function followByCount(r::Robot, side::HorizonSide, count::Int)
    for i = 1:count
        moveStep(r, side, false)
    end
end

function moveStep(r::Robot, side::HorizonSide, needMark::Bool, check=false)
    # Take step in one direction
    move!(r, side)
    if needMark
        if check
            if ismarker(r)
                return 0
            end
        end
        putmarker!(r)
    end
end

function inverseSide(side::HorizonSide)
    # Return inverse side
    return HorizonSide(mod(Int(side)+2, 4))
end

function find_border!(r::Robot, direction_to_border::HorizonSide, direction_of_movement::HorizonSide)
    while isborder(r,direction_to_border)==false
        if isborder(r,direction_of_movement)==false
            move!(r,direction_of_movement)
        else
            move!(r,direction_to_border)
            direction_of_movement=inverse(direction_of_movement)
        end
    end
    #УТВ: непосредственно справа от Робота - внутренняя пергородка
end

function measureBoard(r::Robot)
    stepsToY = moveSide(r, Nord)
    stepsToX = moveSide(r, West)

    width = moveSide(r, Ost)
    height = moveSide(r, Sud)

    moveSide(r, Nord)
    moveSide(r, West)

    followByCount(r, Sud, stepsToY)
    followByCount(r, Ost, stepsToX)

    return (width, height, stepsToX, stepsToY)
end

max_stars_repo_name: TheDarkMoon/mirea-progs

max_issues_repo_name: TheDarkMoon/mirea-progs

max_forks_repo_name: TheDarkMoon/mirea-progs


text:
function hermite_polynomial_coefficient_matrix(left,right)
    [1.0   left   left^2   left^3
     0.0   1.0   2.0*left  3.0*(left^2)
     1.0   right  right^2  right^3
     0.0   1.0   2.0*right 3.0*right^2]
end

function hermite_polynomial_coefficients(left,right)
    M = hermite_polynomial_coefficient_matrix(left,right)
    coeffs = M\I
end

function hermite_polynomials(x,left,right)
    coeffs = hermite_polynomial_coefficients(left,right)
    order,numbasis = size(coeffs)
    return [polynomial_from_coefficients(x,coeffs[:,i]) for i = 1:numbasis]
end

ext: jl


text:
using CSV, DataFrames, Arrow, sem, ModelingToolkit,
    LinearAlgebra, SparseArrays, DataFrames, Optim, LineSearches,
    BenchmarkTools

function get_data_paths(config)
    data_paths = []
    for i = 1:nrow(config)
        row = config[i, :]
        file = string("nfact_", row.nfact_vec, "_nitem_", row.nitem_vec, ".arrow")
        push!(data_paths, file)
    end
    return data_paths
end

function read_files(dir, data_paths)
    data = Vector{Any}()
    for i in 1:length(data_paths)
        push!(data, DataFrames.DataFrame(Arrow.Table(dir*"/"*data_paths[i])))
    end
    return data
end

function gen_model(nfact, nitem, data, start_val)
    nfact = Int64(nfact)
    nitem = Int64(nitem)
    # observed
    #semobserved = SemObsCommon(data = Matrix{Float64}(data))

    semobserved = SemObsMissing(data)

    ## Model definition
    nobs = nfact*nitem
    nnod = nfact+nobs
    @ModelingToolkit.variables x[1:Int64(nobs + nobs)], m[1:nobs]

    #F
    Ind = collect(1:nobs)
    Jnd = collect(1:nobs)
    V = fill(1,nobs)
    F = sparse(Ind, Jnd, V, nobs, nnod)

    #S
    Ind = collect(1:nnod)
    Jnd = collect(1:nnod)
    V = [x[1:nobs]; fill(1.0, nfact)]
    S = sparse(Ind, Jnd, V, nnod, nnod)

    #A
    Ind = collect(1:nobs)
    Jnd = vcat([fill(nobs+i, nitem) for i in 1:nfact]...)
    V = x[(nobs+1):(2*nobs)]
    A = sparse(Ind, Jnd, V, nnod, nnod)

    M = [m..., fill(0.0, nfact)...]

    res_ind = (nobs+1):(2*nobs)
    load_ind = 1:nobs
    n_est = length(start_val)
    mean_ind = (n_est-nobs-nfact+1):(n_est-nfact)

    ind = [collect(res_ind)..., collect(load_ind)..., collect(mean_ind)...]

    start_val = start_val[ind]

    # imply
    imply = ImplySymbolic(A, S, F, [x..., m...], start_val; M = M)

    # loss
    #loss = Loss([SemML(semobserved, [0.0], similar(start_val))])
    loss = Loss([SemFIML(semobserved, imply, 0.0, 0.0)])

    #grad_ml = sem.∇SemML(A, S, F, x, start_val)

#=     diff_ana =
        SemAnalyticDiff(
            LBFGS(),
            Optim.Options(
                ;f_tol = 1e-10,
                x_tol = 1.5e-8),
                (grad_ml,))   =#

    grad_fiml = sem.∇SemFIML(semobserved, imply, A, S, F, [x..., m...], start_val; M = M)

#=     diff_ana =
        SemAnalyticDiff(
            LBFGS(
                alphaguess = LineSearches.InitialHagerZhang(),
                linesearch = LineSearches.HagerZhang()
            ),
            Optim.Options(
                ;f_tol = 1e-10,
                x_tol = 1.5e-8,
                show_trace=true),
            (grad_fiml,)) =#

    diff_ana =
            SemAnalyticDiff(
                :LD_LBFGS,
                nothing,
                (grad_fiml,))

#=     diff_fin = SemFiniteDiff(
        LBFGS(
            ;alphaguess = LineSearches.InitialHagerZhang(),
            linesearch = LineSearches.HagerZhang()
        ),
        Optim.Options(;
            f_tol = 1e-10,
            x_tol = 1.5e-8)) =#

    return (Sem(semobserved, imply, loss, diff_ana), ind)
end

function gen_models(config, data_vec, start_vec)
    models = []
    for i = 1:nrow(config)
        row = config[i, :]
        model = gen_model(row.nfact_vec, row.nitem_vec, Matrix(data_vec[i]),
        convert(Vector{Float64}, start_vec[i].est))
        push!(models, model)
    end
    return models
end

function benchmark_models(models)
    benchmarks = []
    for model in models
        bm = @benchmark sem_fit($model)
        push!(benchmarks, bm)
    end
    return benchmarks
end

function get_fits(models)
    fits = []
    for model in models
        fit = sem_fit(model)
        push!(fits, fit)
    end
    return fits
end

function check_solution(fits, par_vec, par_order; tol = 0.05)
    approx_equal = []
    for (fit, par, ind) in zip(fits, par_vec, par_order)
        is_equ =
            all(
                abs.(fit.minimizer .- par.est[ind]) .<
                tol*abs.(par.est[ind]))
        push!(approx_equal, is_equ)
    end
    return approx_equal
end


text:
@recipe(TeXImg, origin, tex) do scene
    merge(
        default_theme(scene),
        Attributes(
            color = AbstractPlotting.automatic,
            implant = true,
            dpi = 3000.0,
            align = (:left, :center)
        )
    )
end

function AbstractPlotting.convert_arguments(::Type{<: TeXImg}, bbox::T, x::AbstractString) where T <: Rect2D
    println("as")
    return (bbox, CachedTeX(implant_math(x)))
end

function AbstractPlotting.convert_arguments(::Type{<: TeXImg}, pos::T, x::AbstractString) where T <: Point2
    println("sa")
    doc = CachedTeX(implant_math(x))
    dims = doc.raw_dims
    return (Rect(pos..., dims.width, dims.height), CachedTeX(implant_math(x)))
end

function AbstractPlotting.convert_arguments(::Type{<: TeXImg}, bbox::RT, x::LaTeXString) where RT <: Rect2D
    if first(x) == "\$" && last(x) == "\$"
        return (bbox, CachedTeX(implant_math(x)))
    else
        return (bbox, CachedTeX(implant_text(x)))
    end
end

function AbstractPlotting.convert_arguments(::Type{<: TeXImg}, bbox::RT, doc::TeXDocument) where RT <: Rect2D
    println("doc")
    return (bbox, CachedTeX(doc))
end

function AbstractPlotting.convert_arguments(::Type{<: TeXImg}, pos::Point2, doc::CachedTeX)
    println("hi")
    dims = doc.raw_dims
    return (Rect(pos..., dims.width, dims.height), doc)
end

function AbstractPlotting.plot!(img::T) where T <: TeXImg

    bbox   = img[1]
    texdoc = img[2][]

    halign, valign = lift(first, img.align), lift(last, img.align)

    png = rsvg2img(texdoc.handle)

    xr, yr = Node{LinRange}(LinRange(0, 1, 10)), Node{LinRange}(LinRange(0, 1, 10))

    lift(bbox, halign, valign) do bbox, halign, valign
        x0, y0 = origin(bbox)
        w, h   = widths(bbox)

        if halign == :left
            x0 -= 0
        elseif halign == :center
            x0 -= w / 2
        elseif halign == :right
            x0 -= w
        end

        if valign == :top
            y0 -= h
        elseif valign == :center
            y0 -= h / 2
        elseif valign == :bottom
            y0 -= 0
        end

        xr[] = LinRange(x0, x0 + w, size(png, 1))
        yr[] = LinRange(y0, y0 + h, size(png, 2))
    end

    image!(img, xr, yr, png)
end

function get_ink_extents(surf::CairoSurface)
    dims = zeros(Float64, 4)

    ccall(
        (:cairo_recording_surface_ink_extents, CairoMakie.LIB_CAIRO),
        Cvoid,
        (Ptr{Cvoid}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}),
        surf.ptr, Ref(dims, 1), Ref(dims, 2), Ref(dims, 3), Ref(dims, 4)
    )

    return dims
end

function CairoMakie.draw_plot(scene::Scene, screen::CairoMakie.CairoScreen, img::T) where T <: MakieTeX.TeXImg

    bbox = AbstractPlotting.boundingbox(img)

    ctx = screen.context
    tex = img[2][]
    halign, valign = img.align[]

    handle = svg2rsvg(tex.svg)
    dims = tex.raw_dims

    surf, rctx = rsvg2recordsurf(handle)
    x0, y0, w, h = get_ink_extents(surf)

    pos = CairoMakie.project_position(scene, origin(img[1][]), img.model[])
    scale_factor = CairoMakie.project_scale(scene, widths(bbox), img.model[])

    pos = if halign == :left
        pos .+ (0, 0)
    elseif halign == :center
        pos .+ (scale_factor[1] / w / 2, 0)
    elseif halign == :right
        pos .+ (scale_factor[1] / w, 0)
    end

    pos = if valign == :top
        pos .- (0, scale_factor[2] / h)
    elseif valign == :center
        pos .- (0, scale_factor[2] / h / 2)
    elseif valign == :bottom
        pos .- (0, 0)
    end

    Cairo.save(ctx)
    Cairo.translate(
        ctx,
        pos[1],
        pos[2] - (h + y0) * scale_factor[2] / h
    )
    Cairo.scale(
        ctx,
        scale_factor[1] / w,
        scale_factor[2] / h
    )
    render_surface(ctx, surf)
    Cairo.restore(ctx)
end


text:
using BenchmarkTools
const SUITE = BenchmarkGroup()
include("bench_van_der_pol.jl")
include("bench_morris_lecar.jl")


text:
"""
    BenchmarkConfig

A `BenchmarkConfig` contains the configuration for the benchmarks to be executed
by [`benchmarkpkg`](@ref).

This includes the following:

* The commit of the package the benchmarks are run on.
* What julia command should be run, i.e. the path to the Julia executable and
  the command flags used (e.g. optimization level with `-O`).
* Custom environment variables (e.g. `JULIA_NUM_THREADS`).
"""
struct BenchmarkConfig
    id::Union{String,Nothing}
    juliacmd::Cmd
    env::Dict{String,Any}
end

"""
    BenchmarkConfig(;id::Union{String, Nothing} = nothing,
                     juliacmd::Cmd = `joinpath(Sys.BINDIR, Base.julia_exename())`,
                     env::Dict{String, Any} = Dict{String, Any}())

Creates a `BenchmarkConfig` from the following keyword arguments:
* `id` - A git identifier like a commit, branch, tag, "HEAD", "HEAD~1" etc.
         If `id == nothing` then benchmark will be done on the current state
         of the repo (even if it is dirty).
* `juliacmd` - Used to exectue the benchmarks, defaults to the julia executable
               that the Pkgbenchmark-functions are called from. Can also include command flags.
* `env` - Contains custom environment variables that will be active when the
          benchmarks are run.

# Examples
```julia
julia> using Pkgbenchmark

julia> BenchmarkConfig(id = "performance_improvements",
                       juliacmd = `julia -O3`,
                       env = Dict("JULIA_NUM_THREADS" => 4))
BenchmarkConfig:
    id: performance_improvements
    juliacmd: `julia -O3`
    env: JULIA_NUM_THREADS => 4
```
"""
function BenchmarkConfig(;id::Union{String,Nothing} = nothing,

                 juliacmd::Cmd = `$(joinpath(Sys.BINDIR, Base.julia_exename()))`,
                 env::Dict = Dict{String,Any}())
    BenchmarkConfig(id, juliacmd, env)
end

BenchmarkConfig(cfg::BenchmarkConfig) = cfg
BenchmarkConfig(str::String) = BenchmarkConfig(id = str)
BenchmarkConfig(::Nothing) = BenchmarkConfig()

function BenchmarkConfig(d::Dict)
    BenchmarkConfig(
        d["id"],
        Cmd(d["juliacmd"]),
        d["env"]
    )
end

# Arr!...
function Base.Cmd(d::Dict)
    Cmd(
        Cmd(convert(Vector{String}, d["exec"])),
        d["ignorestatus"],
        d["flags"],
        d["env"],
        d["dir"],
    )
end

const _INDENT = "    "

function Base.show(io::IO, bcfg::BenchmarkConfig)
    println(io, "BenchmarkConfig:")
    println(io, _INDENT, "id: ", bcfg.id)
    println(io, _INDENT, "juliacmd: ", bcfg.juliacmd)
    print(io, _INDENT, "env: ")
    if !isempty(bcfg.env)
        first = true
        for (k, v) in bcfg.env
            if !first
                println(io)
                print(io, _INDENT, " "^strwidth("env: "))
            end
            first = false
            print(io, k, " => ", v)
        end
    end
end


text:
# TODO

struct PhilipsFile <: MRIFile
  filename::String
end

function MRIBase.RawAcquisitionData(f::PhilipsFile)
  error("Not yet implemented!")
end


text:
# This file is part of the ValidatedNumerics.jl package; MIT licensed

## Powers

# CRlibm does not contain a correctly-rounded ^ function for Float64
# Use the BigFloat version from MPFR instead, which is correctly-rounded:

# Write explicitly like this to avoid ambiguity warnings:
for T in (:Integer, :Rational, :Float64, :BigFloat, :Interval)

    @eval ^(a::Interval{Float64}, x::$T) = convert(Interval{Float64}, big53(a)^x)
end

# Integer power:

function ^(a::Interval{BigFloat}, n::Integer)
    isempty(a) && return a
    n == 0 && return one(a)
    n == 1 && return a
    n == 2 && return sqr(a)
    n < 0 && a == zero(a) && return emptyinterval(a)

    T = BigFloat

    if isodd(n) # odd power
        isentire(a) && return a
        if n > 0
            a.lo == zero(T) && return @round(T, zero(T), a.hi^n)
            a.hi == zero(T) && return @round(T, a.lo^n, zero(T))
            return @round(T, a.lo^n, a.hi^n)
        else
            if a.lo ≥ zero(T)
                a.lo == zero(T) && return @round(T, a.hi^n, convert(T, Inf))
                return @round(T, a.hi^n, a.lo^n)

            elseif a.hi ≤ zero(T)
                a.hi == zero(T) && return @round(T, convert(T, -Inf), a.lo^n)
                return @round(T, a.hi^n, a.lo^n)
            else
                return entireinterval(a)
            end
        end

    else # even power
        if n > 0
            if a.lo ≥ zero(T)
                return @round(T, a.lo^n, a.hi^n)
            elseif a.hi ≤ zero(T)
                return @round(T, a.hi^n, a.lo^n)
            else
                return @round(T, mig(a)^n, mag(a)^n)
            end

        else
            if a.lo ≥ zero(T)
                return @round(T, a.hi^n, a.lo^n)
            elseif a.hi ≤ zero(T)
                return @round(T, a.lo^n, a.hi^n)
            else
                return @round(T, mag(a)^n, mig(a)^n)
            end
        end
    end
end

function sqr{T<:Real}(a::Interval{T})
    isempty(a) && return a
    if a.lo ≥ zero(T)
        return @round(T, a.lo^2, a.hi^2)

    elseif a.hi ≤ zero(T)
        return @round(T, a.hi^2, a.lo^2)
    end

    return @round(T, mig(a)^2, mag(a)^2)
end

# Floating-point power of a BigFloat interval:
function ^(a::Interval{BigFloat}, x::AbstractFloat)
    T = BigFloat
    domain = Interval(zero(T), convert(T, Inf))

    if a == zero(a)
        a = a ∩ domain
        x > zero(x) && return zero(a)
        return emptyinterval(a)
    end

    isinteger(x) && return a^(round(Int, x))
    x == one(T)/2 && return sqrt(a)

    a = a ∩ domain
    (isempty(x) || isempty(a)) && return emptyinterval(a)

    xx = convert(Interval{BigFloat}, x)

    lo = @round(T, a.lo^xx.lo, a.lo^xx.lo)
    lo1 = @round(T, a.lo^xx.hi, a.lo^xx.hi)
    hi = @round(T, a.hi^xx.lo, a.hi^xx.lo)
    hi1 = @round(T, a.hi^xx.hi, a.hi^xx.hi)

    lo = hull(lo, lo1)
    hi = hull(hi, hi1)

    return hull(lo, hi)
end

function ^{T<:Integer,}(a::Interval{Rational{T}}, x::AbstractFloat)
    a = Interval(a.lo.num/a.lo.den, a.hi.num/a.hi.den)
    a = a^x
    convert(Interval{Rational{T}}, a)
end

# Rational power
function ^{S<:Integer}(a::Interval{BigFloat}, r::Rational{S})
    T = BigFloat
    domain = Interval(zero(a.lo), convert(T, Inf))

    if a == zero(a)
        a = a ∩ domain
        r > zero(r) && return zero(a)
        return emptyinterval(a)
    end

    isinteger(r) && return convert(Interval{T}, a^round(S,r))
    r == one(S)//2 && return sqrt(a)

    a = a ∩ domain
    (isempty(r) || isempty(a)) && return emptyinterval(a)

    y = convert(Interval{BigFloat}, r)

    a^y
end

# Interval power of an interval:
function ^(a::Interval{BigFloat}, x::Interval)
    T = BigFloat
    domain = Interval(zero(T), convert(T, Inf))

    a = a ∩ domain

    (isempty(x) || isempty(a)) && return emptyinterval(a)

    lo1 = a^x.lo
    lo2 = a^x.hi
    lo1 = hull(lo1, lo2)

    hi1 = a^x.lo
    hi2 = a^x.hi
    hi1 = hull(hi1, hi2)

    hull(lo1, hi1)
end

function sqrt{T}(a::Interval{T})
    domain = Interval(zero(T), convert(T, Inf))
    a = a ∩ domain

    isempty(a) && return a

    @round(T, sqrt(a.lo), sqrt(a.hi))  # `sqrt` is correctly-rounded
end

for f in (:exp, :expm1)
    @eval begin
        function ($f){T}(a::Interval{T})
            isempty(a) && return a
            Interval( ($f)(a.lo, RoundDown), ($f)(a.hi, RoundUp) )
        end
    end
end

for f in (:exp2, :exp10)

    @eval function ($f)(x::BigFloat, r::RoundingMode)  # add BigFloat functions with rounding:
            setrounding(BigFloat, r) do
                ($f)(x)
            end
        end

    @eval ($f)(a::Interval{Float64}) = float($f(big53(a)))  # no CRlibm version

    @eval function ($f)(a::Interval{BigFloat})
            isempty(a) && return a
            Interval( ($f)(a.lo, RoundDown), ($f)(a.hi, RoundUp) )
        end
end

for f in (:log, :log2, :log10, :log1p)

    @eval function ($f){T}(a::Interval{T})
            domain = Interval(zero(T), convert(T, Inf))
            a = a ∩ domain

            (isempty(a) || a.hi ≤ zero(T)) && return emptyinterval(a)

            Interval( ($f)(a.lo, RoundDown), ($f)(a.hi, RoundUp) )
        end
end


text:
function parse_commands(filename::AbstractString)
    @assert splitext(filename)[2] == ".json"
    parse_commands(JSON.parsefile(filename))
end

function parse_commands(d::Dict)
    output = ImagineSignal[]
    rig = d[METADATA_KEY]["rig"]
    for typ_key in (ANALOG_KEY, DIGITAL_KEY)
        for k in keys(d[typ_key])
            v = d[typ_key][k]
            seq = haskey(v, "sequence") ? convert(RLEVector, v["sequence"]) : convert(RLEVector, [])
            push!(output, _parse_command(rig, k, v["daq channel"], seq, d[COMPONENT_KEY], Int(d[METADATA_KEY]["samples per second"])*Unitful.s^-1))
        end
    end
    return output
end

function parse_command(filename::AbstractString, comname::AbstractString)
    @assert splitext(filename)[2] == ".json"
    parse_command(JSON.parsefile(filename; dicttype=Dict, use_mmap=true), comname)
end

function parse_command(d::Dict, comname::AbstractString)
    rig = d[METADATA_KEY]["rig"]
    for typ_key in (ANALOG_KEY, DIGITAL_KEY)
        ad = d[typ_key]
        if haskey(ad, comname)
            seq = haskey(ad[comname], "sequence") ? convert(RLEVector, ad[comname]["sequence"]) : convert(RLEVector, [])
            return _parse_command(rig, comname, ad[comname]["daq channel"], seq, d[COMPONENT_KEY], Int(d[METADATA_KEY]["samples per second"])*Unitful.s^-1)
        end
    end
    error("Command signal name not found")
end

#In the JSON arrays, waveforms and counts-of-waves are specified in alternating order: count,wave,count,wave...
function _parse_command(rig::AbstractString, chan_name::AbstractString, daq_chan_name::AbstractString, seqs_compressed::RLEVector, seqs_lookup::Dict, sample_rate::HasInverseTimeUnits)
    sampmapper = default_samplemapper(rig, daq_chan_name; sample_rate = sample_rate)
    rawtyp = rawtype(sampmapper)
    vectype = isoutput(daq_chan_name, rig) ? RLEVector{rawtyp} : AbstractVector{rawtyp}
    seqlist = Vector{vectype}(undef, 0)
    seqnames = String[]
    if !isoutput(daq_chan_name, rig) && !isempty(seqs_compressed)
        @warn "Found samples written to the .json file for input channel $(chan_name).  Loading anyway, but this probably indicates a problem"
    end
    for s in seqs_compressed
        for c = 1:s.n
            push!(seqlist, seqs_lookup[s.value])
            push!(seqnames, s.value)
        end
    end
    cumlen = zeros(Int, length(seqlist))
    calc_cumlength!(cumlen, seqlist)
    for s in unique(seqnames) #Convert from Any vectors to RLEVectors
        seqs_lookup[s] = convert(vectype, seqs_lookup[s])
    end
    return ImagineSignal{vectype}(chan_name, daq_chan_name, rig, seqlist, seqnames, seqs_lookup, cumlen, sampmapper)
end

"""
    ai = parse_ai(ainame; imaginename=splitext(ainame)[1]*".imagine")
    ai = parse_ai(ainame, header::Dict{String})

Parse the imagine header and specified `.ai` file to extract the signals used to
represent timing information in an Imagine recording.
`ai` is a vector of memory-mapped signals.
See `examples/reading_ai.jl` for a usage demonstration.
"""
function parse_ai(ai_name::AbstractString; imagine_header::AbstractString = splitext(ai_name)[1]*".imagine")
    if !isfile(ai_name) && isfile(ai_name*".ai")
        ai_name = ai_name*".ai"
    end
    if !isfile(imagine_header)
        error(".imagine header not found.  Please specify header file name by keyword argument.")
    end
    hdr = ImagineFormat.parse_header(imagine_header)
    parse_ai(ai_name, hdr)
end

function parse_ai(ai_name::AbstractString, hdr::Dict{String})
    chns = hdr["channel list"] #note: zero-based, need to add a constant for mapping these indices to DAQ channels
    chns = map(x-> "AI$(x)", chns)
    labs = split(hdr["label list"], "\$")
    rig = hdr["rig"]
    samp_rate = convert(Int, hdr["scan rate"]) * Unitful.s^-1
    if length(chns) != length(labs)
        error("Invalid .imagine header: The number of channels does not match the number of channel labels")
    end
    parse_ai(ai_name, chns, labs, rig, samp_rate)
end

function parse_ai(ai_name, chns, labels, rig, sample_rate::HasInverseTimeUnits)
    nchannels = length(chns)
    tmp = rigtemplate(rig; sample_rate = sample_rate)
    incoms = getanalog(getinputs(tmp))
    aitype = rawtype(incoms[1])
    nbytes = filesize(ai_name)
    nsamples = nbytes/nchannels/sizeof(aitype)
    if !isapprox(nsamples, round(Int, nsamples))
        error("The .ai file doesn't have the correct number of samples.  The recording may have been interrupted.")
    else
        nsamples = round(Int, nsamples)
    end
    f = open(ai_name, "r")
    A = Mmap.mmap(f, Matrix{aitype}, (nchannels,nsamples))

    parse_ai(A, chns, labels, rig, sample_rate)
end

function parse_ai(ai_name, chns, rig, sample_rate::HasInverseTimeUnits)
    labels = map(x->DEFAULT_DAQCHANS_TO_NAMES[rig][x], chns)
    parse_ai(ai_name, chns, labels, rig, sample_rate)
end

function parse_ai(A::Matrix{T}, chns, labels, rig, sample_rate::HasInverseTimeUnits) where T
    nchannels = size(A,1)
    tmp = rigtemplate(rig; sample_rate = sample_rate)
    incoms = getanalog(getinputs(tmp))
    output = ImagineSignal[]
    for i = 1:length(chns)
        comi = finddaqchan(incoms, chns[i])
        if comi == 0
            @warn "DAQ channel $(chns[i]) is not an analog input channel for this rig.  Attempting to load it anyway."
        end
        com = incoms[comi]
        samps = view(A, i, :)
        #sampsarr = Array{AbstractVector}(undef, 0)
        sampsarr = Array{typeof(samps)}(undef, 0)
        push!(sampsarr, samps)
        lookup_nm = string(hash(chns[i]))
        #mon = ImagineSignal{AbstractVector}(labels[i], chns[i], rig, sampsarr, [lookup_nm;], Dict{String,Any}(lookup_nm=>samps), [length(samps);], mapper(com))
        mon = ImagineSignal(labels[i], chns[i], rig, sampsarr, [lookup_nm;], Dict{String,Any}(lookup_nm=>samps), [length(samps);], mapper(com))
        push!(output, mon)
    end
    return output
end

function parse_ai(A::Matrix, chns, rig, sample_rate::HasInverseTimeUnits)
    labels = map(x->DEFAULT_DAQCHANS_TO_NAMES[rig][x], chns)
    parse_ai(A, chns, labels, rig, sample_rate)
end

"""
    di = parse_di(diname; imaginename=splitext(diname)[1]*".imagine", kwargs...)

Parse the imagine header and specified `.di` file to extract the signals used to
represent timing information in an Imagine recording.
`di` is a vector of memory-mapped signals.
"""
function parse_di(di_name::AbstractString; imagine_header = splitext(di_name)[1]*".imagine", kwargs...)
    if !isfile(di_name) && isfile(di_name*".di")
        di_name = di_name*".di"
    end
    if !isfile(imagine_header)
        error(".imagine header not found.  Please specify header file name by keyword argument.")
    end
    hdr = ImagineFormat.parse_header(imagine_header)
    return parse_di(di_name, hdr; kwargs...)
end

"""
    di = parse_di(diname, header::Dict{String}; load_unused=false)

Parse the specified `.di` file to extract the signals used to
represent timing information in an Imagine recording.
`di` is a vector of memory-mapped signals. `header` is an Imagine header.
"""
function parse_di(di_name::AbstractString, hdr::Dict{String,Any}; load_unused = false)
    chns = hdr["di channel list"] #note: zero-based, need to add a constant for mapping these indices to DAQ channels
    chns = map(x->"P0.$(x)", chns)
    labs = split(hdr["di label list"], "\$")
    rig = hdr["rig"]
    samp_rate = convert(Int, hdr["di scan rate"]) * Unitful.s^-1
    nchannels = length(chns)
    if nchannels != length(labs)
        error("Invalid .imagine header: The number of channels does not match the number of channel labels")
    end
    for i = 1:length(chns)
        biti = findfirst(x->x==parse(Int, split(chns[i], ".")[2]), hdr["di channel list"])
        if biti != findfirst(x->x==chns[i], DI_CHANS[rig])
            @warn "DI channel list entry #$(biti) found in the .imagine header does not match the expected entry for this rig. Attempting to load anyway, but please report this issue"
        end
    end
    parse_di(di_name, chns, labs, rig, samp_rate::HasInverseTimeUnits; load_unused = load_unused)
end

#Does not require Imagine file, but allows a loading a subset of the 8 channels based on names
function parse_di(di_name, chns, labels, rig, sample_rate::HasInverseTimeUnits; load_unused = false)
    nchannels = length(labels)
    di_sigs = parse_di(di_name, rig, sample_rate)
    di_sigs_used = ImagineSignal[]
    #assign labels to di_sigs
    for i = 1:length(chns)
        sig = getdaqchan(di_sigs, chns[i])
        if isfree(sig)
            rename!(sig, String(labels[i]))
        elseif String(labels[i]) != name(sig) && labels[i] != "unused"
            @warn "The name $(String(labels[i])) was provided for a di channel that cannot be renamed.  Defaulting to its allowed name"
        end
        if labels[i] != "unused"
            push!(di_sigs_used, sig)
        end
    end
    if load_unused
        return di_sigs
    else
        return di_sigs_used
    end
end

#Does not require imagine file, loads all 8 channels
function parse_di(di_name, rig, sample_rate::HasInverseTimeUnits)
    insigs = getdigital(getinputs(rigtemplate(rig; sample_rate = sample_rate)))
    nsamples = filesize(di_name) #each sample is one byte
    f = open(di_name, "r")
    rd_nsamples = floor(Int, nsamples/64)*64 #to truncate excessive samples that doesn't fit into page size chunks
    A = Mmap.mmap(f, BitArray, (8,rd_nsamples)) #we need all 8 bits even if the user didn't ask for 8 channels of DI
    output = ImagineSignal[]
    close(f)
    for i = 1:length(insigs)
        sig = insigs[i]
        daq_chan_str = daq_channel(sig)
        biti = findfirst(x->x==daq_chan_str, DI_CHANS[rig])
        samps = mappedarray(UInt8, view(A, biti, :)) #feels a bit dishonest to call UInt8 the raw type, but there is no Bit type...
        sampsarr = Array{AbstractVector}(undef, 0)
        push!(sampsarr, samps)
        lookup_nm = string(hash(daq_chan_str))
        mpr = SampleMapper(UInt8(false), UInt8(true), 0.0*Unitful.V, 3.3*Unitful.V, false, true, sample_rate) #unlike for outputs, the raw type is Bool
        mon = ImagineSignal{AbstractVector}(name(sig), daq_chan_str, rig, sampsarr, [lookup_nm;], Dict{String,Any}(lookup_nm=>samps), [length(samps);], mpr)
        push!(output, mon)
    end
    return output
end

function append_or_replace!(coms::Vector{ImagineSignal}, newcoms)
    for c in newcoms
        i = findname(coms, name(c))
        if i != 0
            coms[i] = c
        else
            append!(coms, c)
        end
    end
    return coms
end

function load_signals(any_name::AbstractString)
    if !isfile(any_name)
        error("File not found")
    end
    extn = splitext(any_name)[2]
    basenm = splitext(any_name)[1]
    if extn != ".json" && !isfile(basenm * ".json") && !isfile(basenm * ".imagine")
        error("A matching .json or .imagine header was not found in the supplied directory, so the experiment cannot be loaded")
    end

    coms = ImagineSignal[]
    comnm = basenm * ".json"
    has_comfile = false
    if isfile(comnm)
        has_comfile = true
        append!(coms, parse_commands(comnm))
    else
        @warn "A matching .json command file was not found in the supplied directory.  Attempting to load .ai and .di file anyway"
    end
    ainm = basenm * ".ai"
    dinm = basenm * ".di"
    if isempty(coms) && !isfile(ainm) && !isfile(dinm)
        error("Cannot load anything: no .json file was found, and no .ai or .di file was found to match the .imagine file")
    end

    if isfile(ainm)
        ai_sigs = parse_ai(ainm)
        if has_comfile
            for s in ai_sigs
                if !is_similar(s, getname(coms, name(s)))
                    error("AI file signal $(name(s)) does not match the entry found in the .json command file")
                end
            end
        end
        append_or_replace!(coms, ai_sigs)
    else
        @warn "Analog input file (with .ai extension) was not found in the supplied directory"
    end

    if isfile(dinm)
        di_sigs = parse_di(dinm)
        if has_comfile
            for s in di_sigs
                if !is_similar(s, getname(coms, name(s)))
                    error("DI file signal $(name(s)) does not match the entry found in the .json command file")
                end
            end
        end
        append_or_replace!(coms, di_sigs)
    else
        @warn "Digital input file (with .di extension) was not found in the supplied directory"
    end
    return coms
end

max_forks_repo_name: HolyLab/ImagineInterface


text:
# Parsing expressions
Meta.parse("x + y")
#=
julia> Meta.parse("x + y")
:(x + y)
=#

# Expr type
Meta.parse("x + y") |> typeof
#=
julia> Meta.parse("x + y") |> typeof
Expr
=#

# Expr object dump
Meta.parse("x + y") |> dump
#=
julia> Meta.parse("x + y") |> dump
Expr
  head: Symbol call
  args: Array{Any}((3,))
    1: Symbol +
    2: Symbol x
    3: Symbol y
=#

# -----------------------------------------------------------------------------

# Single variable
Meta.parse("x") |> dump
#=
julia> Meta.parse("x") |> dump
Symbol x
=#

# Function calls with keyword arguments
Meta.parse("""open("/tmp/test.txt", read = true, write = true)""") |> dump
#=
julia> Meta.parse("""open("/tmp/test.txt", read = true, write = true)""") |> dump
Expr
  head: Symbol call
  args: Array{Any}((4,))
    1: Symbol open
    2: String "/tmp/test.txt"
    3: Expr
      head: Symbol kw
      args: Array{Any}((2,))
        1: Symbol read
        2: Bool true
    4: Expr
      head: Symbol kw
      args: Array{Any}((2,))
        1: Symbol write
        2: Bool true
=#

# Nested functions
Meta.parse("cos(sin(x+1))") |> dump
#=
julia> Meta.parse("cos(sin(x+1))") |> dump
Expr
  head: Symbol call
  args: Array{Any}((2,))
    1: Symbol cos
    2: Expr
      head: Symbol call
      args: Array{Any}((2,))
        1: Symbol sin
        2: Expr
          head: Symbol call
          args: Array{Any}((3,))
            1: Symbol +
            2: Symbol x
            3: Int64 1
=#

# -----------------------------------------------------------------------------

# constructing Expr manually
Expr(:call, :+, :x, :y)
#=
julia> Expr(:call, :+, :x, :y)
:(x + y)
=#

# nested one
Expr(:call, :sin, Expr(:call, :+, :x, :y))
#=
julia> Expr(:call, :sin, Expr(:call, :+, :x, :y))
:(sin(x + y))
=#

# using :() syntax
ex = :(x + y)
dump(ex)
#=
julia> ex = :(x + y)
:(x + y)

julia> dump(ex)
Expr
  head: Symbol call
  args: Array{Any}((3,))
    1: Symbol +
    2: Symbol x
    3: Symbol y
=#

# quoted block
:( begin
        x = 1
        y = 2
   end )
#=
julia> :(begin
       x = 1
       y = 2
       end)
quote
    #= REPL[24]:2 =#
    x = 1
    #= REPL[24]:3 =#
    y = 2
end
=#

# using quote block directly
quote
    x = 1
    y = 2
end
#=
julia> quote
           x = 1
           y = 2
       end
quote
    #= REPL[25]:2 =#
    x = 1
    #= REPL[25]:3 =#
    y = 2
end
=#

# -----------------------------------------------------------------------------
# more complex expression

# Variable Assignment
:(x = 1 + 1) |> dump
#=
julia> :(x = 1 + 1) |> dump
Expr
  head: Symbol =
  args: Array{Any}((2,))
    1: Symbol x
    2: Expr
      head: Symbol call
      args: Array{Any}((3,))
        1: Symbol +
        2: Int64 1
        3: Int64 1
=#

# Code Block
:(begin
        println("hello")
        println("world")
    end) |> dump
#=
julia> :(begin
           println("hello")
           println("world")
       end) |> dump
Expr
  head: Symbol block
  args: Array{Any}((4,))
    1: LineNumberNode
      line: Int64 2
      file: Symbol REPL[10]
    2: Expr
      head: Symbol call
      args: Array{Any}((2,))
        1: Symbol println
        2: String "hello"
    3: LineNumberNode
      line: Int64 3
      file: Symbol REPL[10]
    4: Expr
      head: Symbol call
      args: Array{Any}((2,))
        1: Symbol println
        2: String "world"
=#

# Conditional
:(if 2 > 1
        "good"
    else
        "bad"
    end) |> dump
#=
julia> :(if 2 > 1
           "good"
       else
           "bad"
       end) |> dump
Expr
  head: Symbol if
  args: Array{Any}((3,))
    1: Expr
      head: Symbol call
      args: Array{Any}((3,))
        1: Symbol >
        2: Int64 2
        3: Int64 1
    2: Expr
      head: Symbol block
      args: Array{Any}((2,))
        1: LineNumberNode
          line: Int64 2
          file: Symbol REPL[12]
        2: String "good"
    3: Expr
      head: Symbol block
      args: Array{Any}((2,))
        1: LineNumberNode
          line: Int64 4
          file: Symbol REPL[12]
        2: String "bad"
=#

# Loop
:(for i in 1:5
        println("hello world")
    end) |> dump
#=
julia> :(for i in 1:5
           println("hello world")
       end) |> dump
Expr
  head: Symbol for
  args: Array{Any}((2,))
    1: Expr
      head: Symbol =
      args: Array{Any}((2,))
        1: Symbol i
        2: Expr
          head: Symbol call
          args: Array{Any}((3,))
            1: Symbol :
            2: Int64 1
            3: Int64 5
    2: Expr
      head: Symbol block
      args: Array{Any}((2,))
        1: LineNumberNode
          line: Int64 2
          file: Symbol REPL[13]
        2: Expr
          head: Symbol call
          args: Array{Any}((2,))
            1: Symbol println
            2: String "hello world"
=#

# Function Definition
:(function foo(x; y = 1)
      return x + y
  end) |> dump
#=
julia> :(function foo(x; y = 1)
           return x + y
       end) |> dump
Expr
  head: Symbol function
  args: Array{Any}((2,))
    1: Expr
      head: Symbol call
      args: Array{Any}((3,))
        1: Symbol foo
        2: Expr
          head: Symbol parameters
          args: Array{Any}((1,))
            1: Expr
              head: Symbol kw
              args: Array{Any}((2,))
                1: Symbol y
                2: Int64 1
        3: Symbol x
    2: Expr
      head: Symbol block
      args: Array{Any}((2,))
        1: LineNumberNode
          line: Int64 2
          file: Symbol REPL[23]
        2: Expr
          head: Symbol return
          args: Array{Any}((1,))
            1: Expr
              head: Symbol call
              args: Array{Any}((3,))
                1: Symbol +
                2: Symbol x
                3: Symbol y
=#

# -----------------------------------------------------------------------------

# eval expressions
eval(:(x = 1))
#=
julia> eval(:(x = 1))
1

julia> x
1
=#

# eval from within a function
function foo()
    eval(:(y = 1))
end
#=
julia> function foo()
           eval(:(y = 1))
       end

julia> foo()
1

julia> y
1
=#

# interpolation
x = 2
:(sqrt($x))
#=
julia> x = 2
2

julia> :(sqrt($x))
:(sqrt(2))
=#

# splatting interpolation
v = [1, 2, 3]
quote
    max($(v...))
end
#=
julia> v = [1, 2, 3]
3-element Array{Int64,1}:
 1
 2
 3

julia> quote
           max($(v...))
       end
quote
    #= REPL[47]:2 =#
    max(1, 2, 3)
end
=#

# incorrect splatting usage
quote
    max($v...)
end
#=
julia> quote
           max($v...)
       end
quote
    #= REPL[46]:2 =#
    max([1, 2, 3]...)
end
=#

# -----------------------------------------------------------------------------

# what is QuoteNode
:( x = :hello ) |> dump
#=
julia> :( x = :hello ) |> dump
Expr
  head: Symbol =
  args: Array{Any}((2,))
    1: Symbol x
    2: QuoteNode
      value: Symbol hello
=#

# variable assignment
:( x = hello ) |> dump
#=
julia> :( x = hello ) |> dump
Expr
  head: Symbol =
  args: Array{Any}((2,))
    1: Symbol x
    2: Symbol hello
=#

# variable assignment (using interpolation)
sym = :hello
:( x = $sym)
#=
julia> sym = :hello
:hello

julia> :( x = $sym)
:(x = hello)
=#

# using QuoteNode
sym = QuoteNode(:hello)
:( x = $sym)
#=
julia> sym = QuoteNode(:hello)
:(:hello)

julia> :( x = $sym)
:(x = :hello)
=#

# -----------------------------------------------------------------------------

# interpolation into nested expression
:( :( x = 1 ) ) |> dump
#=
julia> :( :( x = 1 ) ) |> dump
Expr
  head: Symbol quote
  args: Array{Any}((1,))
    1: Expr
      head: Symbol =
      args: Array{Any}((2,))
        1: Symbol x
        2: Int64 1

julia> :( :( x = $($v) ) )
:($(Expr(:quote, :(x = $(Expr(:$, 2))))))
=#

:( :( x = $v ) ) |> dump
#=
julia> :( :( x = $v ) ) |> dump
Expr
  head: Symbol quote
  args: Array{Any}((1,))
    1: Expr
      head: Symbol =
      args: Array{Any}((2,))
        1: Symbol x
        2: Expr
          head: Symbol $
          args: Array{Any}((1,))
            1: Symbol v
=#


text:
type CallableFixedRateBond{DC <: DayCount, P <: PricingEngine, P2 <: PricingEngine} <: AbstractCallableBond
  lazyMixin::LazyMixin
  bondMixin::BondMixin
  faceAmount::Float64
  schedule::Schedule
  cashflows::FixedRateLeg
  dc::DC
  redemption::Float64
  startDate::Date
  pricingEngine::P
  settlementValue::Float64
  putCallSchedule::CallabilitySchedule
  blackEngine::P2
  blackVolQuote::Quote
end

function CallableFixedRateBond{DC <: DayCount, P <: PricingEngine}(settlementDays::Int, faceAmount::Float64, schedule::Schedule, coupons::Union{Vector{Float64}, Float64},
                              accrualDayCounter::DC, paymentConvention::BusinessDayConvention, redemption::Float64, issueDate::Date,
                              putCallSchedule::CallabilitySchedule, pe::P)
  maturityDate = schedule.dates[end]

  # build bond
  isZeroCouponBond = length(coupons) == 1 && is_close(coupons[1], 0.0)

  if ~isZeroCouponBond
    coups = FixedRateLeg(schedule, faceAmount, coupons, schedule.cal, paymentConvention, accrualDayCounter)
  else
    # build redemption CashFlow
    redemptionDate = adjust(schedule.cal, paymentConvention, maturityDate)
    coups = ZeroCouponLeg(SimpleCashFlow(redemption, redemptionDate))
  end

  blackVolQuote = Quote(0.0)
  blackEngine = BlackCallableFixedRateBondEngine(blackVolQuote)

  return CallableFixedRateBond{DC, P, BlackCallableFixedRateBondEngine}(LazyMixin(), BondMixin(settlementDays, issueDate, maturityDate), faceAmount, schedule, coups, accrualDayCounter, redemption,
                              schedule.dates[1], pe, 0.0, putCallSchedule, blackEngine, blackVolQuote)
end

function accrued(bond::CallableFixedRateBond, d::Date)
  if d == Date()
    d = settlement_date(bond)
  end

  for i in eachindex(bond.cashflows.coupons)
    # the first coupon paying after d is the one we are after
    if ~has_occurred(bond.cashflows[i], d, false)
      if isa(bond.cashflows[i], Coupon)
        return accrued_amount(bond.cashflows[i], d) / notional(bond, d) * 100.0
      else
        return 0.0
      end
    end
  end

  return 0.0
end

type CallableBondArgs{F <: Frequency}
  couponDates::Vector{Date}
  couponAmounts::Vector{Float64}
  redemption::Float64
  redemptionDate::Date
  freq::F
  putCallSchedule::CallabilitySchedule
  callabilityPrices::Vector{Float64}
  callabilityDates::Vector{Date}
end

function CallableBondArgs(bond::CallableFixedRateBond)
  redemption = get_redemption(bond)
  settlement = settlement_date(bond)
  coupDates = Vector{Date}()
  coupAmts = Vector{Float64}()
  for i = 1:length(bond.cashflows.coupons)-1
    if ~has_occurred(bond.cashflows[i], settlement, false)
      push!(coupDates, date(bond.cashflows[i]))
      push!(coupAmts, amount(bond.cashflows[i]))
    end
  end

  callabilityPrices = Vector{Float64}()
  callabilityDates = Vector{Date}()
  for i in eachindex(bond.putCallSchedule)
    if ~has_occurred(bond.putCallSchedule[i], settlement, false)
      push!(callabilityDates, bond.putCallSchedule[i].date)
      push!(callabilityPrices, bond.putCallSchedule[i].price.amount)

      # calling accrued() forces accrued interest to be zero if future option date
      # is also coupon date, so that dirty price = clean price.  Used here because
      # callability is always applied before coupon in tree engine
      if isa(bond.putCallSchedule[i].price.callType, CleanCall)
        callabilityPrices[end] += accrued(bond, bond.putCallSchedule[i].date)
      end
    end
  end

  return CallableBondArgs(coupDates, coupAmts, amount(redemption), date(redemption), get_frequency(bond), bond.putCallSchedule,
                          callabilityPrices, callabilityDates)
end


text:
using DataStructures
using Base.Test

# construction

@test is(typeof(OrderedSet()), OrderedSet{Any})
@test is(typeof(OrderedSet([1,2,3])), OrderedSet{Int})
@test is(typeof(OrderedSet{Int}([3])), OrderedSet{Int})
data_in = (1, "banana", ())
s = OrderedSet(data_in)
data_out = collect(s)
@test is(typeof(data_out), Array{Any,1})
@test is(tuple(data_out...), data_in)
@test is(tuple(data_in...), tuple(s...))
@test length(data_out) == length(data_in)

# hash
s1 = OrderedSet{ASCIIString}(["bar", "foo"])
s2 = OrderedSet{ASCIIString}(["foo", "bar"])
s3 = OrderedSet{ASCIIString}(["baz"])
@test hash(s1) != hash(s2)
@test hash(s1) != hash(s3)

# isequal
@test  isequal(OrderedSet(), OrderedSet())
@test !isequal(OrderedSet(), OrderedSet([1]))
@test  isequal(OrderedSet{Any}(Any[1,2]), OrderedSet{Int}([1,2]))
@test !isequal(OrderedSet{Any}(Any[1,2]), OrderedSet{Int}([1,2,3]))

@test  isequal(OrderedSet{Int}(), OrderedSet{AbstractString}())
@test !isequal(OrderedSet{Int}(), OrderedSet{AbstractString}([""]))
@test !isequal(OrderedSet{AbstractString}(), OrderedSet{Int}([0]))
@test !isequal(OrderedSet{Int}([1]), OrderedSet{AbstractString}())
@test  isequal(OrderedSet{Any}([1,2,3]), OrderedSet{Int}([1,2,3]))
@test  isequal(OrderedSet{Int}([1,2,3]), OrderedSet{Any}([1,2,3]))
@test !isequal(OrderedSet{Any}([1,2,3]), OrderedSet{Int}([1,2,3,4]))
@test !isequal(OrderedSet{Int}([1,2,3]), OrderedSet{Any}([1,2,3,4]))
@test !isequal(OrderedSet{Any}([1,2,3,4]), OrderedSet{Int}([1,2,3]))
@test !isequal(OrderedSet{Int}([1,2,3,4]), OrderedSet{Any}([1,2,3]))

# eltype, similar
s1 = similar(OrderedSet([1,"hello"]))
@test isequal(s1, OrderedSet())
@test is(eltype(s1), Any)
s2 = similar(OrderedSet{Float32}([2.0f0,3.0f0,4.0f0]))
@test isequal(s2, OrderedSet())
@test is(eltype(s2), Float32)

# show
@test endswith(sprint(show, OrderedSet()), "OrderedSet{Any}()")
@test endswith(sprint(show, OrderedSet(['a'])), "OrderedSet{Char}(['a'])")

s = OrderedSet(); push!(s,1); push!(s,2); push!(s,3)
@test !isempty(s)
@test in(1,s)
@test in(2,s)
@test length(s) == 3
push!(s,1); push!(s,2); push!(s,3)
@test length(s) == 3
@test pop!(s,1) == 1
@test !in(1,s)
@test in(2,s)
@test length(s) == 2
@test_throws KeyError pop!(s,1)
@test pop!(s,1,:foo) == :foo
@test length(delete!(s,2)) == 1
@test !in(1,s)
@test !in(2,s)
@test pop!(s) == 3
@test length(s) == 0
@test isempty(s)

# copy
data_in = (1,2,9,8,4)
s = OrderedSet(data_in)
c = copy(s)
@test isequal(s,c)
v = pop!(s)
@test !in(v,s)
@test  in(v,c)
push!(s,100)
push!(c,200)
@test !in(100,c)
@test !in(200,s)

# sizehint, empty
s = OrderedSet([1])
@test isequal(sizehint(s, 10), OrderedSet([1]))
@test isequal(empty!(s), OrderedSet())
# TODO: rehash

# start, done, next
for data_in in ((7,8,4,5),
                ("hello", 23, 2.7, (), [], (1,8)))
    s = OrderedSet(data_in)

    s_new = OrderedSet()
    for el in s
        push!(s_new, el)
    end
    @test isequal(s, s_new)

    t = tuple(s...)

    @test is(t, data_in)
    @test length(t) == length(s)
    for (e,f) in zip(t,s)
        @test is(e,f)
    end
end

# union
@test isequal(union(OrderedSet([1])),OrderedSet([1]))
s = ∪(OrderedSet([1,2]), OrderedSet([3,4]))
@test isequal(s, OrderedSet([1,2,3,4]))
s = union(OrderedSet([5,6,7,8]), OrderedSet([7,8,9]))
@test isequal(s, OrderedSet([5,6,7,8,9]))
s = OrderedSet([1,3,5,7])
union!(s,(2,3,4,5))
# TODO: order is not the same, so isequal should return false...
@test isequal(s,OrderedSet([1,2,3,4,5,7]))

# intersect
@test isequal(intersect(OrderedSet([1])),OrderedSet([1]))
s = ∩(OrderedSet([1,2]), OrderedSet([3,4]))
@test isequal(s, OrderedSet())
s = intersect(OrderedSet([5,6,7,8]), OrderedSet([7,8,9]))
@test isequal(s, OrderedSet([7,8]))
@test isequal(intersect(OrderedSet([2,3,1]), OrderedSet([4,2,3]), OrderedSet([5,4,3,2])), OrderedSet([2,3]))

# setdiff
@test isequal(setdiff(OrderedSet([1,2,3]), OrderedSet()),        OrderedSet([1,2,3]))
@test isequal(setdiff(OrderedSet([1,2,3]), OrderedSet([1])),     OrderedSet([2,3]))
@test isequal(setdiff(OrderedSet([1,2,3]), OrderedSet([1,2])),   OrderedSet([3]))
@test isequal(setdiff(OrderedSet([1,2,3]), OrderedSet([1,2,3])), OrderedSet())
@test isequal(setdiff(OrderedSet([1,2,3]), OrderedSet([4])),     OrderedSet([1,2,3]))
@test isequal(setdiff(OrderedSet([1,2,3]), OrderedSet([4,1])),   OrderedSet([2,3]))
s = OrderedSet([1,3,5,7])
setdiff!(s,(3,5))
@test isequal(s,OrderedSet([1,7]))
s = OrderedSet([1,2,3,4])
setdiff!(s, OrderedSet([2,4,5,6]))
@test isequal(s,OrderedSet([1,3]))

# ordering
@test OrderedSet() < OrderedSet([1])
@test OrderedSet([1]) < OrderedSet([1,2])
@test !(OrderedSet([3]) < OrderedSet([1,2]))
@test !(OrderedSet([3]) > OrderedSet([1,2]))
@test OrderedSet([1,2,3]) > OrderedSet([1,2])
@test !(OrderedSet([3]) <= OrderedSet([1,2]))
@test !(OrderedSet([3]) >= OrderedSet([1,2]))
@test OrderedSet([1]) <= OrderedSet([1,2])
@test OrderedSet([1,2]) <= OrderedSet([1,2])
@test OrderedSet([1,2]) >= OrderedSet([1,2])
@test OrderedSet([1,2,3]) >= OrderedSet([1,2])
@test !(OrderedSet([1,2,3]) >= OrderedSet([1,2,4]))
@test !(OrderedSet([1,2,3]) <= OrderedSet([1,2,4]))

# issubset, symdiff
for (l,r) in ((OrderedSet([1,2]),     OrderedSet([3,4])),
              (OrderedSet([5,6,7,8]), OrderedSet([7,8,9])),
              (OrderedSet([1,2]),     OrderedSet([3,4])),
              (OrderedSet([5,6,7,8]), OrderedSet([7,8,9])),
              (OrderedSet([1,2,3]),   OrderedSet()),
              (OrderedSet([1,2,3]),   OrderedSet([1])),
              (OrderedSet([1,2,3]),   OrderedSet([1,2])),
              (OrderedSet([1,2,3]),   OrderedSet([1,2,3])),
              (OrderedSet([1,2,3]),   OrderedSet([4])),
              (OrderedSet([1,2,3]),   OrderedSet([4,1])))
    @test issubset(intersect(l,r), l)
    @test issubset(intersect(l,r), r)
    @test issubset(l, union(l,r))
    @test issubset(r, union(l,r))
    @test isequal(union(intersect(l,r),symdiff(l,r)), union(l,r))
end
@test ⊆(OrderedSet([1]), OrderedSet([1,2]))

## TODO: not implemented for OrderedSets
#@test ⊊(OrderedSet([1]), OrderedSet([1,2]))
#@test !⊊(OrderedSet([1]), OrderedSet([1]))
#@test ⊈(OrderedSet([1]), OrderedSet([2]))

# TODO: returns false!
#       == is not properly defined for OrderedSets
#@test symdiff(OrderedSet([1,2,3,4]), OrderedSet([2,4,5,6])) == OrderedSet([1,3,5,6])
@test isequal(symdiff(OrderedSet([1,2,3,4]), OrderedSet([2,4,5,6])), OrderedSet([1,3,5,6]))

# filter
s = OrderedSet([1,2,3,4])
@test isequal(filter(isodd,s), OrderedSet([1,3]))
filter!(isodd, s)
@test isequal(s, OrderedSet([1,3]))

# first
# TODO: throws BoundsError in v0.3, ArgumentError in v0.4
#@test_throws ArgumentError first(OrderedSet())
@test first(OrderedSet([2])) == 2

#####################

# empty set
d = OrderedSet{Char}()
@test length(d) == 0
@test isempty(d)
@test !('c' in d)
push!(d, 'c')
@test !isempty(d)
empty!(d)
@test isempty(d)

# access, modification

for c in 'a':'z'
    push!(d, c)
end

for c in 'a':'z'
    @test c in d
end

@test collect(d) == collect('a':'z')


text:
module WellTesting

using RecipesBase
using Plots
import LinearSolvers.Reg
using Statistics

export CentralDiff

include("CentralDiff.jl")
include("BuildUpRecipes.jl")
include("DrawDownRecipes.jl")

end # module


text:
facts("Containers") do
  context("incorrect type creation") do
    @fact_throws Exception Container("this"; volume="/a")
  end

  context("type creation") do
    container = Container("this"; volume=("/this", "/a"), volume=("/that", "/b"))
    context("standard") do
      @fact container.image.name => "this"
      @fact length(container.volumes) => 2
      @fact container.volumes => x -> haskey(x, "/a")
      @fact container.volumes => x -> haskey(x, "/b")
      @fact container.volumes["/a"] => "/this"
      @fact container.volumes["/b"] => "/that"

      for name in names(RudeOil.Container)
        if name != :image && name != :volumes
          @fact getfield(container, name) => RudeOil.DEFAULT_CONTAINER[name]
        end
      end

      @fact Container().image.name => RudeOil.DEFAULT_CONTAINER[:image].name
    end

    context("copy creation") do
      copy = Container(container; volume=("/t", "/a"), interactive=true)

      @fact copy.image.name => "this"
      @fact copy.interactive => true
      @fact length(copy.volumes) => 2
      @fact copy.volumes => x -> haskey(x, "/a")
      @fact copy.volumes => x -> haskey(x, "/b")
      @fact copy.volumes["/a"] => "/t"
      @fact copy.volumes["/b"] => "/that"

      for name in names(RudeOil.Container)
        if name ∉ (:image, :volumes, :interactive)
          @fact getfield(copy, name) => RudeOil.DEFAULT_CONTAINER[name]
        end
      end
    end
  end

  context("Run container") do
    vm = activate(machine)

    context("Hello World") do
      result = vm |> Container() |> `echo "Hello World"` |> readchomp
      @fact result => "Hello World"
      result = vm |> `echo "Hello World"` |> readchomp
      @fact result => "Hello World"
    end

    context("With Image") do
      python = image("python", packages=["python"])
      result = vm |> python |> `python -c "print('Hello World')"` |> readchomp
      @fact result => "Hello World"
    end

    context("With volumes") do
      volume = image("volume", volume="/a")
      dirname = RudeOil.Util.heretempfile()
      try
        mkdir(dirname)
        vm |> volume |> Container(volume=(dirname, "/a"), workdir="/a") |> `touch hello` |> run
        @fact joinpath(dirname, "hello") => isfile
      finally
        rm(joinpath(dirname, "hello"))
        rm(dirname)
      end
    end

    context("With array of commands") do
      volume = image("volume", volume="/a")
      dirname = RudeOil.Util.heretempfile()
      try
        mkdir(dirname)
        vm |> volume |> Container(volume=(dirname, "/a"), workdir="/a") |> [
          `touch hello`,
          `mkdir small`,
          setenv(`touch world`, dir="small"),
        ] |> run
        @fact joinpath(dirname, "hello") => isfile
        @fact joinpath(dirname, "small", "world") => isfile
      finally
        rm(dirname, recursive=true)
      end
    end
  end
end


text:
"""
    function make_mc_time_converter(event_info::Union{MCEventInfo,DAQEventInfo})

Returns a function which converts MC time to JTE time.
"""
function make_mc_time_converter(event_info::Union{MCEventInfo,DAQEventInfo})
    function time_converter(time)
        return time - (event_info.timestamp * 1e9 + event_info.nanoseconds) + event_info.mc_time
    end
    return time_converter
end

"""
    function cherenkov_origin(pos, t::Track)

Calculate the origin of the Cherenkov photon on a track.
"""
function cherenkov_origin(pos, t::Track)
    θ = acos(1/N_SEAWATER)
    P = project(pos, t)
    dir = -normalize(t.dir)
    track_distance = pld3(pos, t.pos, t.dir)
    distance = track_distance / tan(θ)
    Position(P + distance*dir)
end


text:
module ModalC0RefFEsTests

using Test
using Gridap
using Gridap.ReferenceFEs
using Gridap.FESpaces
using Gridap.Fields
using Gridap.CellData

# using BenchmarkTools

order = 1
p = QUAD
T = VectorValue{2,Float64}

m = ReferenceFE(p,modalC0,T,order)
l = ReferenceFE(p,lagrangian,T,order)

test_reference_fe(m)

@test num_dofs(m) == num_dofs(l)
@test Conformity(m) === Conformity(l)
@test get_face_own_dofs(m,Conformity(m)) == get_face_own_dofs(l,Conformity(l))
@test get_face_own_dofs_permutations(m,Conformity(m)) == get_face_own_dofs_permutations(l,Conformity(l))

# domain = (0,1,0,1)
# partition = (2,2)
# model = CartesianDiscreteModel(domain,partition)

function _test_function_interpolation(order,u,V)
  test_single_field_fe_space(V)
  uh = interpolate(u,V)
  Ω = Triangulation(model)
  degree = 2*order
  dΩ = Measure(Ω,degree)
  l2(u) = sqrt(sum( ∫( u⊙u )*dΩ ))
  e = u - uh
  el2 = l2(e)
  @test el2 < 1.0e-9
  # writevtk(Ω,"results",nsubcells=20,cellfields=["uh"=>uh])
  # free_vals = zeros(num_free_dofs(V)); free_vals[7] = 1
  # uh = FEFunction(V,free_vals)
  # writevtk(Ω,"shape_7",nsubcells=20,cellfields=["s7"=>uh])
  # free_vals = zeros(num_free_dofs(V)); free_vals[9] = 1
  # uh = FEFunction(V,free_vals)
  # writevtk(Ω,"shape_9",nsubcells=20,cellfields=["s9"=>uh])
  # free_vals = zeros(num_free_dofs(V)); free_vals[11] = 1
  # uh = FEFunction(V,free_vals)
  # writevtk(Ω,"shape_11",nsubcells=20,cellfields=["s11"=>uh])
  # free_vals = zeros(num_free_dofs(V)); free_vals[13] = 1
  # uh = FEFunction(V,free_vals)
  # writevtk(Ω,"shape_13",nsubcells=20,cellfields=["s13"=>uh])
end

function test_function_interpolation(::Type{T},order,C,u) where T
  reffe = ReferenceFE(modalC0,T,order)
  V = FESpace(model,reffe,conformity=C)
  _test_function_interpolation(order,u,V)
end

function test_function_interpolation(::Type{T},order,C,u,bboxes,space) where T
  reffe = ReferenceFE(modalC0,T,order,bboxes,space=space)
  V = FESpace(model,reffe,conformity=C)
  _test_function_interpolation(order,u,V)
end

domain = (0,1)
partition = (1,)
model = CartesianDiscreteModel(domain,partition)
trian = Triangulation(model)

T = Float64; order = 3; C = :H1; u(x) = x[1]^3
bboxes = [Point(0.0),Point(1.0),Point(0.0),Point(1.0),Point(-3.0),Point(1.0)]
bboxes = CellPoint(fill(bboxes,1),trian,PhysicalDomain())
test_function_interpolation(T,order,C,u,bboxes.cell_ref_point,:Q)

domain = (0,1)
partition = (4,)
model = CartesianDiscreteModel(domain,partition)

T = Float64; order = 3; C = :H1; u(x) = x[1]^3
bboxes = [ [Point(0.0),Point(1.0),Point(0.0),Point(1.0),Point(-1.0),Point(3.0)],
           [Point(0.0),Point(1.0),Point(0.0),Point(1.0),Point(-1.5),Point(3.2)],
           [Point(0.0),Point(1.0),Point(0.0),Point(1.0),Point(-0.2),Point(1.0)],
           [Point(0.0),Point(1.0),Point(0.0),Point(1.0),Point(-1.2),Point(2.0)] ]
test_function_interpolation(T,order,C,u,bboxes,:Q)

domain = (0,4)
partition = (4,)
model = CartesianDiscreteModel(domain,partition)

T = Float64; order = 3; C = :H1; u(x) = x[1]^3
bboxes = [ [Point(0.0),Point(1.0),Point(0.0),Point(1.0),Point(-1.0),Point(3.0)],
           [Point(0.0),Point(1.0),Point(0.0),Point(1.0),Point(-1.5),Point(3.2)],
           [Point(0.0),Point(1.0),Point(0.0),Point(1.0),Point(-0.2),Point(1.0)],
           [Point(0.0),Point(1.0),Point(0.0),Point(1.0),Point(-1.2),Point(2.0)] ]
test_function_interpolation(T,order,C,u,bboxes,:Q)

domain = (0,1,0,1)
partition = (2,2,)
model = CartesianDiscreteModel(domain,partition)
trian = Triangulation(model)

T = Float64; order = 3; C = :H1; u(x) = (x[1]+x[2])^3
bboxes = reshape( [Point(0.0,0.0),Point(1.0,1.0),Point(0.0,0.0),
                   Point(1.0,1.0),Point(0.0,0.0),Point(1.0,1.0),
                   Point(0.0,0.0),Point(1.0,1.0),Point(-0.5,2.5),
                   Point(2.5,2.5),Point(-0.5,2.5),Point(2.5,2.5),
                   Point(-1.0,-1.0),Point(-1.0,1.25),Point(-1.0,-1.0),
                   Point(-1.0,1.25),Point(0.0,0.0),Point(1.0,1.0)], 18 )
bboxes = CellPoint(fill(bboxes,4),trian,PhysicalDomain())
test_function_interpolation(T,order,C,u,bboxes.cell_ref_point,:Q)

T = Float64; order = 5; C = :H1; u(x) = (x[1]+x[2])^5
bboxes = reshape( [Point(0.0,0.0),Point(1.0,1.0),Point(0.0,0.0),
                   Point(1.0,1.0),Point(0.0,0.0),Point(1.0,1.0),
                   Point(0.0,0.0),Point(1.0,1.0),Point(-0.5,2.5),
                   Point(2.5,2.5),Point(-0.5,2.5),Point(2.5,2.5),
                   Point(-1.0,-1.0),Point(-1.0,1.25),Point(-1.0,-1.0),
                   Point(-1.0,1.25),Point(0.0,0.0),Point(1.0,1.0)], 18 )
bboxes = CellPoint(fill(bboxes,4),trian,PhysicalDomain())
test_function_interpolation(T,order,C,u,bboxes.cell_ref_point,:S)

order = 1; T = Float64; C = :L2; u(x) = x[1]+x[2]
test_function_interpolation(T,order,C,u)

order = 1; T = VectorValue{2,Float64}; C = :H1
u(x) = VectorValue(x[1]+x[2],x[2])
test_function_interpolation(T,order,C,u)

domain = (0,1,0,1,0,1)
partition = (2,2,2)
model = CartesianDiscreteModel(domain,partition)

order = 1; T = Float64; C = :H1; u(x) = x[1]+x[2]+x[3]
test_function_interpolation(T,order,C,u)

# Inspect operator matrix to check if L2-scalar product of
# gradients of bubble functions satisfy Kronecker's delta
# domain = (0,1)
# partition = (1)
# model = CartesianDiscreteModel(domain,partition)
# order = 6; T = Float64; C = :H1;
# reffe = ReferenceFE(modalC0,T,order)
# V = FESpace(model,reffe,conformity=C)
# Ω = Triangulation(model)
# degree = 2*order
# dΩ = LebesgueMeasure(Ω,degree)
# a(u,v) = ∫( ∇(v)⊙∇(u) )*dΩ
# b(v) = 0.0
# op = AffineFEOperator(bboxes,V,V)

end # module


text:
using OrderedCollections, Test

@testset "OrderedDict" begin

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

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

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

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

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

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

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

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

        @test od60[14] == 15
    end

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

    # OrderedDict

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

    @testset "Issue #30" begin
        d = OrderedDict(:a=>1, :b=>2)
        d1 = OrderedDict(k=>v for (k,v) in d)
        @test keytype(d1) == keytype(d)
        @test valtype(d1) == valtype(d)
    end

    @testset "Sorting" begin
        d = Dict(i=>Char(123-i) for i = 1:10)
        @test collect(keys(d)) != 1:10
        sd = sort!(OrderedDict(d))
        @test collect(keys(sd)) == 1:10
        @test collect(values(sd)) == collect('z':-1:'q')
        @test sort(sd) == sd
        sdv = sort!(OrderedDict(d); byvalue=true)
        @test collect(keys(sdv)) == 10:-1:1
        @test collect(values(sdv)) == collect('q':'z')
    end

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

    @testset "map!(f, values(OrderedDict))" begin
        testdict = OrderedDict(:a=>1, :b=>2)
        map!(v->v-1, values(testdict))
        @test testdict[:a] == 0
        @test testdict[:b] == 1
    end

    @testset "Issue #47" begin
        @test eltype(OrderedDict(String => :string, SubString => :substring)) == Pair{Type,Symbol}
        @test eltype(OrderedDict(:string => String, :substring => SubString)) == Pair{Symbol,Type}
        @test eltype(OrderedDict(String => String, SubString => SubString)) == Pair{Type,Type}

        @test eltype(OrderedDict(tuple(String => :string, SubString => :substring))) == Pair{Type,Symbol}
        @test eltype(OrderedDict(tuple(:string => String, :substring => SubString))) == Pair{Symbol,Type}
        @test eltype(OrderedDict(tuple(String => String, SubString => SubString))) == Pair{Type,Type}
    end

    @testset "Issue #71" begin
        od = OrderedDict(Dict(i=>0 for i=1:158))
        sort!(od)
        @test od[158] == 0
    end

    @testset "Issue #71b" begin
        # This is actually a simplified version of #60, which was triggered while fixing #71
        # It doesn't actually fail on previous versions of OrderedCollections
        od = OrderedDict{Int,Int}(13=>13)
        delete!( od, 13 )
        od[14]=14
        @test od[14] == 14
    end

end # @testset OrderedDict


text:
using Comonicon
using GraphUtilities, GenericTensorNetworks, CUDA, Graphs
using DelimitedFiles

########################### GENERATE ###################################
function mkdir_and_dumpcode(datafolder, config; ntrials, niters, nslices, sc_weight, rw_weight, sc_target, prefix)
    folder = foldername(datafolder, config; create=true, prefix)
    instance = problem_instance(config;
        optimizer=TreeSA(; sc_target, sc_weight, ntrials, niters, rw_weight, nslices, βs=0.01:0.05:30),
        simplifier=MergeGreedy()
        )
    println("graph size = $(nv(instance.graph))")
    println("time, space, RW complexitys are $(timespacereadwrite_complexity(instance))")
    GraphUtilities.save_code(folder, instance)
end

function get_graph(graphname, size; degree, filling, seed)
    if graphname=="diag"
        return DiagGraphConfig(; filling, n=size, m=size, seed=seed)
    elseif graphname == "regular"
        return RegularGraphConfig(; degree, size=size, seed=seed)
    elseif graphname == "square"
        return SquareGraphConfig(; filling, m=size, n=size, seed=seed)
    elseif graphname == "misproject"
        return MISProjectGraphConfig(; n=size, index=seed)
    elseif graphname == "mapped-regular"
        return MappedRegularGraphConfig(; degree, size, seed)
    else
        return SmallGraphConfig(graphname)
    end
end

@cast function generate(graphname::String, problem::String, size::Int;
    sizestop::Int=size, sizestep::Int=1,
    degree::Int=3, seed::Int=1, seedstop::Int=seed,
    filling::Float64=0.8,
    weights::Union{Nothing,Vector{Float64}}=nothing,
    datafolder::String="data", prefix::String="",
    ntrials::Int=1, niters::Int=10, nslices::Int=0, sc_target::Int=25,
    sc_weight::Float64=1.0, rw_weight::Float64=1.0,  # config TreeSA
    )
    for sz in size:sizestep:sizestop, sd in seed:seedstop
        println("seed = $sd, size = $(sz)")
        graph = get_graph(graphname, sz; degree, filling, seed=sd)
        config = GraphProblemConfig(; problem, graph, weights, openvertices=Int[])
        mkdir_and_dumpcode(datafolder, config; ntrials, niters, nslices, sc_target, sc_weight, rw_weight, prefix)
    end
end

########################### COMPUTE ###################################
function load_and_compute(datafolder, config, property; prefix, cudadevice)
    folder = foldername(datafolder, config; create=false, prefix)
    instance = GraphUtilities.load_code(config, folder)
    println("time, space, RW complexitys are $(timespacereadwrite_complexity(instance))")
    if cudadevice >=0
        CUDA.device!(cudadevice)
    end
    res = solve(instance, property; usecuda=cudadevice>=0)[]
    save_property(folder, property, res)
end

@cast function compute(graphname::String, problem::String, size::Int, property::String;
    sizestop::Int=size, sizestep::Int=1,
    degree::Int=3, seed::Int=1, seedstop::Int=seed,
    filling::Float64=0.8,
    weights::Union{Nothing,Vector{Float64}}=nothing,
    datafolder::String="data", prefix::String="",
    cudadevice::Int=-1)
    for sz in size:sizestep:sizestop, sd in seed:seedstop
        println("seed = $sd, size = $(sz)")
        graph = get_graph(graphname, sz; degree, filling, seed=sd)
        config = GraphProblemConfig(; problem, graph, weights, openvertices=Int[])
        load_and_compute(datafolder, config, GraphUtilities.parseproperty(property); prefix, cudadevice)
    end
end

########################### OGP ###################################

@cast function ogp(graphname::String, size::Int;
    sizestop::Int=size, sizestep::Int=1,
    filling::Float64=0.8,
    degree::Int=3, seed::Int=1, seedstop::Int=seed,
    alpha::Float64=0.1,
    datafolder::String="data",
    prefix::String="", cleanup::Bool=false, overwrite::Bool=false)
    for sz in size:sizestep:sizestop, sd in seed:seedstop
        println("seed = $sd, size = $(sz)")
        graph = get_graph(graphname, sz; degree, filling, seed=sd)
        config = GraphProblemConfig(; problem="IndependentSet", graph, weights=nothing, openvertices=Int[])
        maxn = load_property(foldername(datafolder, config; create=false, prefix), SizeMax()).n
        K = ceil(Int, maxn * alpha)

        prop = ConfigsMax(K; tree_storage=true, bounded=true)
        folder = foldername(datafolder, config; create=false, prefix)
        fd = joinpath(folder, "$(GraphUtilities.unique_string(prop)).dat")
        if cleanup
            if ispath(fd)
                println("Removing MISs of size `$(maxn-K+1):$(maxn)`: $(fd)")
                rm(fd, recursive=true)
            else
                println("Did not find files for MISs of size `$(maxn-K+1):$(maxn)`: $(fd)")
            end
        else
            if !overwrite && ispath(fd)
                println("Find existing data for MISs of size `$(maxn-K+1):$(maxn): $(fd)`, pass")
            else
                println("Computing MISs of size `$(maxn-K+1):$(maxn)`")
                load_and_compute(datafolder, config, prop; prefix, cudadevice=-1)
            end
        end
    end
end

@cast function hamming(graphname::String, size::Int;
    sizestop::Int=size, sizestep::Int=1,
    filling::Float64=0.8,
    degree::Int=3, seed::Int=1, seedstop::Int=seed,
    datafolder::String="data", overwrite::Bool=false,
    alpha::Float64=0.1, alpha0::Float64=alpha,
    prefix::String="", nsample::Int = 10000)
    for sz in size:sizestep:sizestop, sd in seed:seedstop
        println("seed = $sd, size = $(sz)")
        graph = get_graph(graphname, sz; degree, filling, seed=sd)
        config = GraphProblemConfig(; problem="IndependentSet", graph, weights=nothing, openvertices=Int[])
        folder = foldername(datafolder, config; create=false, prefix)
        maxn = load_property(folder, SizeMax()).n
        K = ceil(Int, maxn * alpha0)
        Kreal = ceil(Int, maxn * alpha)
        fd = joinpath(folder, "hamming-K$(Kreal)-n$(nsample).dat")
        if !overwrite && ispath(fd)
            println("Hamming distance file exists for MISs of size `$(maxn-Kreal+1):$(maxn)`: $fd")
        end
        println("Hamming for MISs of size `$(maxn-Kreal+1):$(maxn)`")
        configs = load_property(folder, ConfigsMax(K; tree_storage=true, bounded=true))
        if K == 1
            tree = configs.c
        else
            tree = sum(configs.coeffs[end-Kreal+1:end])
        end
        samples = generate_samples(tree, 2*nsample);
        hd = hamming_distribution(samples[1:nsample], samples[nsample+1:2*nsample])
        writedlm(fd, hd)
    end
end

@cast function hamming1(graphname::String, size::Int, K::Int;
    sizestop::Int=size, sizestep::Int=1,
    filling::Float64=0.8,
    degree::Int=3, seed::Int=1, seedstop::Int=seed,
    datafolder::String="data", overwrite::Bool=false,
    prefix::String="", nsample::Int = 10000)
    for sz in size:sizestep:sizestop, sd in seed:seedstop
        println("seed = $sd, size = $(sz)")
        graph = get_graph(graphname, sz; degree, filling, seed=sd)
        config = GraphProblemConfig(; problem="IndependentSet", graph, weights=nothing, openvertices=Int[])
        folder = foldername(datafolder, config; create=false, prefix)
        fd = joinpath(folder, "hamming-K$(K)-n$(nsample).dat")
        if !overwrite && ispath(fd)
            println("Hamming distance file exists for MISs of `K = $K`: $fd")
        end
        configs = load_property(folder, ConfigsMax(K; tree_storage=true, bounded=true))
        println("Hamming for MISs of size `$(configs.maxorder-K+1)`")
        tree = configs.coeffs[1]
        samples = generate_samples(tree, 2*nsample);
        hd = hamming_distribution(samples[1:nsample], samples[nsample+1:2*nsample])
        writedlm(fd, hd)
    end
end

@main


text:
using StructDualDynProg, CutPruners
export appendtoSDDPLattice!, updatemaxncuts!

function MathProgBase.linprog(c::DualEntropy, h::EntropyCone, cut::DualEntropy)
    cuthrep = HalfSpace(cut.h, 1)
    MathProgBase.linprog(c.h, intersect(h.poly, cuthrep))
end

function getNLDS(c::DualEntropy, W, h, T, linset, solver, newcut::Symbol, pruningalgo::AbstractCutPruningAlgo)
    K = [(:NonNeg, collect(setdiff(BitSet(1:size(W, 1)), linset))), (:Zero, collect(linset))]
    C = [(:NonNeg, collect(1:size(W, 2)))]
    StructDualDynProg.StructProg.NLDS{Float64}(W, h, T, K, C, c.h, solver, pruningalgo, newcut)
end

function extractNLDS(c, h::EntropyConeLift, id, idp, solver, newcut, pruningalgo::AbstractCutPruningAlgo)
    hr = MixedMatHRep(hrep(h.poly))
    idx  = rangefor(h, id)
    idxp = rangefor(h, idp)
    W = hr.A[:,idx]
    T = hr.A[:,idxp]
    getNLDS(c, W, hr.b, T, hr.linset, solver, newcut, pruningalgo)
end

function next_perm(arr)
    # Find non-increasing suffix
    i = length(arr)
    while i > 1 && arr[i - 1] > arr[i]
        i -= 1
    end
    if i <= 1
        false
    else
        # Find successor to pivot
        j = length(arr)
        while arr[j] < arr[i - 1]
            j -= 1
        end
        arr[i - 1], arr[j] = arr[j], arr[i - 1]

        # Reverse suffix
        arr[i:end] = arr[length(arr):-1:i]
        true
    end
end

function addchildren!(sp::StructDualDynProg.StructProg.StochasticProgram{S}, node::Int, n::Int, old::Bool, oldnodes, newnodes, solver, max_n::Int, newcut::Symbol, pruningalgo::Vector) where S
    proba = 0.0 # It will be set at once the number of children is known (i.e. at the end of the function)
    transitions = []
    N = Int(ntodim(n))
    function addchild(J::EntropyIndex, K::EntropyIndex, adh::Symbol,
                      T = SparseMatrixCSC(1.0LinearAlgebra.I, N, N))
        if (n,J,K,adh) in keys(oldnodes)
            if !old
                child = oldnodes[(n,J,K,adh)]
                StructDualDynProg.SOI.add_scenario_transition!(sp, node, child, proba, T)
            end
        else
            child = create_node(sp, oldnodes, newnodes, n, J, K, adh, node, solver, max_n, newcut, pruningalgo)
            StructDualDynProg.SOI.add_scenario_transition!(sp, node, child, proba, T)
        end
    end

    if true
        #childT = Vector{AbstractMatrix{S}}()
        dn = max_n - n
        N = Int(ntodim(n))
        for i in 1:n-1
            for j in i+1:min(n, i+dn)
                I = set(1:i)
                J = set(1:j)
                push!(transitions, addchild(J,I,:Self))
                σ = collect(1:n)
                done = [(I,J)]
                while next_perm(σ)
                    Iperm = mymap(σ, I, n)
                    Jperm = mymap(σ, J, n)
                    if !((Iperm, Jperm) in done)
                        σinv = collect(1:n)[σ]
                        σi = map(I->mymap(σinv, I, n), indset(n))
                        push!(done, (Iperm, Jperm))
                        # Why transpose ?
                        P = sparse(Vector{Int}(σi), 1:N, 1, N, N)'
                        push!(transitions, addchild(J,I,:Self,P))
                    end
                end
            end
        end
    else
        #childT = nothing
        for adh in [:Self, :Inner]
            for J in indset(n)
                for K in indset(n)
                    if J == K ||
                        (adh == :Self && !(K ⊆ J)) ||
                        #(adh == :Self && !(fullset(n) ⊆ J)) || # Don't do that, Z-Y is not found with max_n = 5 and no inner otherwise
                        (adh == :Self && (fullset(n) ⊆ J)) || # Don't do that, Z-Y is not found with max_n = 5 and no inner otherwise
                        (adh == :Inner && K ⊆ J) ||
                        (adh == :Inner && J ⊆ K) ||
                        (adh == :Inner && (K ∩ J) == 0) ||
                        (adh == :Inner && !(fullset(n) ⊆ (K ∪ J))) ||
                        (adh == :Self && n <= 3) ||
                        (adh == :Inner && n <= 4) ||
                        adh == :Inner
                        continue
                    end
                    if nadh(n, J, K, adh) <= max_n
                        push!(transitions, addchild(J,K,adh))
                    end
                end
            end
        end
    end
    # no optimality cut needed since only the root has an objective
    # so the probably only influence sampling hence equal probability is a fair choice
    for tr in transitions
        StructDualDynProg.SOI.set!(sp, StructDualDynProg.SOI.Probability(), tr, 1 / length(transitions))
    end
end

# np: n for parent
function create_node(sp::StructDualDynProg.StructProg.StochasticProgram{S}, oldnodes, newnodes, np, Jp, Kp, adhp, parent, solver, max_n, newcut, pruningalgo::Vector) where S
    @assert !((np,Jp,Kp,adhp) in keys(oldnodes))
    if !((np,Jp,Kp,adhp) in keys(newnodes))
        n = nadh(np, Jp, Kp, adhp)
        #h = polymatroidcone(np)
        # h is for the parent.
        # We do not need the constraint of the parent so we want the full polyhedron with no constraint
        Np = Int(ntodim(np))
        fullpoly = polyhedron(hrep(HalfSpace{S, SparseVector{S, Int}}[]; d = Np))
        h = EntropyCone{S}(np, fullpoly)
        lift = adhesivelift(h, Jp, Kp, adhp)
        c = constdualentropy(n, 0)
        nlds = extractNLDS(c, lift, 2, 1, solver, newcut, pruningalgo[n])
        newnodedata = StructDualDynProg.StructProg.NodeData(nlds, Np)
        newnode = StructDualDynProg.SOI.add_scenario_node!(sp, newnodedata)
        # Only the root node has a non-zero objective so no need for optimality cuts
        StructDualDynProg.SOI.set!(sp, StructDualDynProg.StructProg.CutGenerator(),
                                   newnode, StructDualDynProg.StructProg.NoOptimalityCutGenerator())
        newnodes[(np,Jp,Kp,adhp)] = newnode
        addchildren!(sp, newnode, n, false, oldnodes, newnodes, solver, max_n, newcut, pruningalgo)
    end
    newnodes[(np,Jp,Kp,adhp)]
end

function fillroot!(sp::StructDualDynProg.StructProg.StochasticProgram, c::DualEntropy, H::EntropyCone, cut::DualEntropy, newnodes, solver, max_n::Integer, newcut::Symbol, pruningalgo::Vector)
    h = MixedMatHRep(hrep(H.poly ∩ HalfSpace(cut.h, 1)))
    @assert h.A isa AbstractSparseMatrix
    T = spzeros(Float64, size(h.A, 1), 0)
    hb = sparsevec(h.b) # FIXME it was done before but not sure it is useful as it the the rhs
    nlds = getNLDS(c, h.A, hb, T, h.linset, solver, newcut, AvgCutPruningAlgo(-1))
    rootdata = StructDualDynProg.StructProg.NodeData(nlds, 0)
    root = StructDualDynProg.SOI.add_scenario_node!(sp, rootdata)
    StructDualDynProg.SOI.set!(sp, StructDualDynProg.StructProg.CutGenerator(),
                               root, StructDualDynProg.StructProg.NoOptimalityCutGenerator())
    newnodes[(H.n,emptyset(),emptyset(),:NoAdh)] = root
    oldnodes = Dict{Tuple{Int,EntropyIndex,EntropyIndex,Symbol},Int}()
    addchildren!(sp, root, H.n, false, oldnodes, newnodes, solver, max_n, newcut, pruningalgo)
end

function StructDualDynProg.SOI.stochasticprogram(num_stages::Int, c::DualEntropy, h::EntropyCone, solver, max_n, cut::DualEntropy, newcut::Symbol, pruningalgo::Vector)
    # allnodes[n][J][K]: if K ⊆ J, it is self-adhesivity, otherwise it is inner-adhesivity
    allnodes = Dict{Tuple{Int,EntropyIndex,EntropyIndex,Symbol},Int}()
    sp = StructDualDynProg.StructProg.StochasticProgram{Float64}(num_stages)
    fillroot!(sp, c, h, cut, allnodes, solver, max_n, newcut, pruningalgo)
    sp, allnodes
end
function Base.append!(sp::StructDualDynProg.StructProg.StochasticProgram, oldnodes, solver, max_n, newcut, pruningalgo::Vector)
    newnodes = Dict{Tuple{Int,EntropyIndex,EntropyIndex,Symbol},Int}()
    for ((n,J,K,adh), node) in oldnodes
        addchildren!(sp, node, nadh(n,J,K,adh), true, oldnodes, newnodes, solver, max_n, newcut, pruningalgo)
    end
    merge!(oldnodes, newnodes)
end
function updatemaxncuts!(sp::StructDualDynProg.StructProg.StochasticProgram, allnodes, maxncuts::Vector{Int})
    for ((n,J,K,adh), node) in allnodes
        StructDualDynProg.updatemaxncuts!(nodedata(sp, node).nlds, maxncuts[nadh(n,J,K,adh)])
    end
end


text:
println()
i = Float64(1.0)

while float(i*(float(1/i))) == 1
    global i = nextfloat(i)
end
println(i)
# 1.000000057228997

function nierowne(x = Float64)
    while float(x * float(1/x)) == 1
        x = nextfloat(x)
    end
    return x
end

println(nierowne(1.0))

max_stars_repo_path: II rok/I semestr/ANM/Cwiczenia 01/zad8.jl

max_stars_repo_name: antonitomaszewski/Studia

max_issues_repo_name: antonitomaszewski/Studia

max_forks_repo_name: antonitomaszewski/Studia


text:
include("../utils.jl")

@testset "SphereSymmetricMatrices" begin
    M = SphereSymmetricMatrices(3)
    M_complex = SphereSymmetricMatrices(3, ℂ)
    A = [1 2 3; 2 4 -5; 3 -5 6] / norm([1 2 3; 2 4 -5; 3 -5 6])
    B = [1 2; 2 4] / norm([1 2; 2 4])                                     #wrong dimensions
    C = [-3 -im 4; im 5 im; 4 -im 0] / norm([-3 -im 4; im 5 im; 4 -im 0]) #complex
    D = [1 2 3; 4 5 6; 7 8 9] / norm([1 2 3; 4 5 6; 7 8 9])               #not symmetric
    E = [1 2 3; 2 4 -5; 3 -5 6]                                           #not of unit norm
    J = [1.0 0.0 0.0; 0.0 0.0 0.0; 0.0 0.0 0.0]
    K = [0.0 1.0 0.0; 0.0 0.0 0.0; 0.0 0.0 0.0]
    @testset "Real Sphere Symmetric Matrices Basics" begin
        @test repr(M) == "SphereSymmetricMatrices(3, ℝ)"
        @test representation_size(M) == (3, 3)
        @test base_manifold(M) === M
        @test typeof(get_embedding(M)) === ArraySphere{Tuple{3,3},ℝ}
        @test check_point(M, A) === nothing
        @test_throws DomainError is_point(M, B, true)
        @test_throws DomainError is_point(M, C, true)
        @test_throws DomainError is_point(M, D, true)
        @test_throws DomainError is_point(M, E, true)
        @test check_vector(M, A, zeros(3, 3)) === nothing
        @test_throws DomainError is_vector(M, A, B, true)
        @test_throws DomainError is_vector(M, A, C, true)
        @test_throws DomainError is_vector(M, A, D, true)
        @test_throws DomainError is_vector(M, D, A, true)
        @test_throws DomainError is_vector(M, A, E, true)
        @test_throws DomainError is_vector(M, J, K, true)
        @test manifold_dimension(M) == 5
        A2 = similar(A)
        @test A == project!(M, A2, A)
        @test A == project(M, A)
        embed!(M, A2, A)
        A3 = embed(M, A)
        @test A2 == A
        @test A3 == A
        F = [2 4 -1; 4 9 7; -1 7 5] / norm([2 4 -1; 4 9 7; -1 7 5])
        G = [-10 9 -8; 9 7 6; -8 6 5] / norm([-10 9 -8; 9 7 6; -8 6 5])
        test_manifold(
            M,
            [A, F, G],
            test_injectivity_radius=false,
            test_forward_diff=false,
            test_reverse_diff=false,
            test_vector_spaces=true,
            test_project_tangent=true,
            test_musical_isomorphisms=true,
            test_default_vector_transport=true,
            is_tangent_atol_multiplier=2,
            test_inplace=true,
        )
    end
    @testset "Complex Sphere Symmetric Matrices Basics" begin
        @test repr(M_complex) == "SphereSymmetricMatrices(3, ℂ)"
        @test manifold_dimension(M_complex) == 8
        H =
            [7.0 -1.5im 0.0; 1.5im 3.0 -1.0im; 0.0 1.0im 4.5] /
            norm([7.0 -1.5im 0.0; 1.5im 3.0 -1.0im; 0.0 1.0im 4.5])
        I =
            [2.0 4.0 5.0im; 4.0 -1.0 -9.0; -5.0im -9.0 2.5] /
            norm([2.0 4.0 5.0im; 4.0 -1.0 -9.0; -5.0im -9.0 2.5])
        test_manifold(
            M_complex,
            [C, H, I],
            test_injectivity_radius=false,
            test_forward_diff=false,
            test_reverse_diff=false,
            test_vector_spaces=true,
            test_project_tangent=true,
            test_musical_isomorphisms=true,
            test_default_vector_transport=true,
            is_tangent_atol_multiplier=2,
            is_point_atol_multiplier=2,
            projection_atol_multiplier=2,
            exp_log_atol_multiplier=2,
            test_inplace=true,
        )
    end
end


text:
using HDF5, JLD, Test, ProgressMeter, Pkg

include("pkgversions.jl");

pkgsVersion()

arr1 = rand(1000,1000)
arr2 = rand(10,25)

h5write("testfile.h5", "arr1", arr1)
# @write "testfile.h5" arr1
# h5write("testfile.h5", "arr1", arr1, "arr2", arr2) for several datasets at the same time
# the first option is more specific (you can tell the name of the dataset to be saved), so I liked more.

arr_load = h5read("testfile.h5", "arr1")

farr = h5open("testfile.h5", "r")

dump(farr)

names(farr)

readmmap(farr["arr1"])

close(farr)

@time arr_load = h5read("testfile.h5", "arr1");

@time data = h5read("testfile.h5", "arr1", (2:15, 2:15))

h5write("testfile.h5", "arr2", arr2)

farr2 = h5open("testfile.h5", "r")

names(farr2)

read(farr2["arr2"])

close(farr2)

function szprint(szfile)
    if szfile <= 10^6
        println("File size = $(szfile/10^3) KB")
    elseif 10^6 < szfile <= 10^9
        println("File size = $(szfile/10^6) MB")
    elseif 10^9 < szfile <= 10^12
        println("File size = $(szfile/10^9) GB")
    elseif 10^12 < szfile <= 10^15
        println("File size = $(szfile/10^12) TB")
    else
        println("File size is really big. Petabytes or more!")
    end
end

szprint(700001000)

f = h5open("incrementalFilev1.h5", "w")
dset_A = d_create(f,"arr1",datatype(Float64), dataspace(10^6,1))
for i in 1:10^4:10^6
    dset_A[i:i + 10^4 - 1,1] = rand(10^4)
    szfile = filesize("incrementalFilev1.h5")
    szprint(szfile)
    sleep(.1)
    IJulia.clear_output(true)
end
flush(f)
close(f)

function szprintp(szfile)
    if szfile <= 10^6
        return (szfile/10^3, "KB")
    elseif 10^6 < szfile <= 10^9
        return (szfile/10^6, "MB")
    elseif 10^9 < szfile <= 10^12
        return (szfile/10^9, "GB")
    elseif 10^12 < szfile <= 10^15
        return(szfile/10^12, "TB")
    else
        return ("=> Petabytes")
    end
end

szprintp(1000)

h5open("incrementalFilev2.h5", "w") do file
    dset = d_create(file, "arr1k", datatype(Float64), dataspace(10^6,1))
    p = Progress(100, dt=0.1, desc="Saving...", color=:blue)
    for iter in 1:10^4:10^6
        dset[iter:iter + 10^4 - 1,1] = rand(10^4)
        szfile = filesize("incrementalFilev2.h5")
        ProgressMeter.next!(p; showvalues = [(:Filesize, szprintp(szfile))])
        sleep(0.1)
        IJulia.clear_output(true)
    end
    end;

f1 = h5open("incrementalFilev1.h5","r")

names(f1)

A_mmaped = readmmap(f1["arr1"])

#rand(100,100)

fm = h5open("MatrixFile.h5", "w")
dset_fm = d_create(fm,"arraym",datatype(Float64), dataspace(100,100,100)) #, "chunk", (100,100,1))
p = Progress(100, dt=0.1, desc="Saving...", color=:blue)
for i in 1:100
    dset_fm[:, :, i] = rand(100,100)
    szfilem = filesize("MatrixFile.h5")
    ProgressMeter.next!(p; showvalues = [(:File_Size, szprintp(szfilem))])
    sleep(0.1)
    IJulia.clear_output(true)
end
flush(fm)
close(fm)

fmatrix = h5open("MatrixFile.h5", "r")

names(fmatrix)

close(fmatrix)

h5read("MatrixFile.h5", "arraym", (:, :, 1))

save("testfile.jld", "arr1", arr1)
#@save "testfile.jld" arr1
# again, I prefer the first option, just for control.

openfile = jldopen("testfile.jld","r", mmaparrays=true)
#load("testfile.jld", "arr1")
#@load "testfile.jld" arr1
# eval(:arr1)
# these last options are also valid, but for big files the first option is better.

@time read(openfile["arr1"])[1:15, 2:5]

@load "testfile.jld"

names(openfile)

close(openfile)

max_stars_repo_name: lazarusA/helpfulJulia

max_issues_repo_name: lazarusA/helpfulJulia

max_forks_repo_name: lazarusA/helpfulJulia


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.11.2
#   kernelspec:
#     display_name: Julia 1.8.0-DEV
#     language: julia
#     name: julia-1.8
# ---

# %%
using Distributions, Plots
pyplot(fmt=:svg)
x ⪅ y = x < y || x ≈ y

chisq(a, b, c, d) = (a + b + c + d)*(a*d - b*c)^2/((a + b)*(c + d)*(a + c)*(b + d))
pval_chisq(a, b, c, d) = ccdf(Chisq(1), chisq(a, b, c, d))

function pval_fisher(a, b, c, d)
    hg = Hypergeometric(a + b, c + d, a + c)
    p0 = pdf(hg, a)
    sum(pdf(hg, x) for x in support(hg) if pdf(hg, x) ⪅ p0)
end

pdfbin(r, s, p, q, (x, y)) = pdf(Binomial(r, p), x) * pdf(Binomial(s, q), y)
supportbin(r, s, p, q) = Iterators.product(support(Binomial(r, p)), support(Binomial(s, q)))
prob_alphaerror(pval, α, r, s, p, q) = sum(pdfbin(r, s, p, q, (a, c))
    for (a, c) in supportbin(r, s, p, q) if pval(a, r-a, c, s-c) < α; init=0.0)

# %%
α = 0.05
p = 0:0.01:1
r, s = 25, 20
pae_chisq = prob_alphaerror.(pval_chisq, α, r, s, p, p)
pae_fisher = prob_alphaerror.(pval_fisher, α, r, s, p, p)

plot(; title="Model: Binomial($r, p)×Binomial($s, p)", titlefontsize=12)
plot!(p, pae_chisq; label="chisq")
plot!(p, pae_fisher; label="fisher")
hline!([α]; color=:black, label="α", alpha=0.5)
plot!(; ylim=(-0.05α, 1.4α), ylabel="probability of alpha errors")
plot!(; xtick=0:0.1:1, xlabel="p")

# %%
a, b, c, d = 1, 5, 24, 15
r, s = a+c, b+d
p = (a+b)*(a+c)/(a+b+c+d)^2
pp = round(p; digits=3)

α = 0:0.001:1
pae_chisq = prob_alphaerror.(pval_chisq, α, r, s, p, p)
pae_fisher = prob_alphaerror.(pval_fisher, α, r, s, p, p)

P = plot(; title="Model: Binomial($r, $pp)×Binomial($s, $pp)", titlefontsize=10)
plot!(α, pae_chisq; label="chisq")
plot!(α, pae_fisher; label="fisher")
plot!([0, 1], [0, 1]; color=:black, label="", alpha=0.5, legend=:topleft)
plot!(; ytick=0:0.1:1, ylabel="probability of alpta errors")
plot!(; xtick=0:0.1:1, xlabel="α")

α = 0:0.0001:0.1
pae_chisq = prob_alphaerror.(pval_chisq, α, r, s, p, p)
pae_fisher = prob_alphaerror.(pval_fisher, α, r, s, p, p)

Q = plot(; title="Model: Binomial($r, $pp)×Binomial($s, $pp)", titlefontsize=10)
plot!(α, pae_chisq; label="chisq")
plot!(α, pae_fisher; label="fisher")
plot!([0, 0.1], [0, 0.1]; color=:black, label="", alpha=0.5, legend=:topleft)
plot!(; ytick=0:0.01:0.1, ylabel="probability of alpha errors")
plot!(; xtick=0:0.01:0.1, xlabel="α")

plot(P, Q; size=(800, 420))

# %%

max_stars_repo_path: 0013/probability of alpha errors of 2x2 chisq test.jl

max_stars_repo_name: genkuroki/public

max_issues_repo_name: genkuroki/public

max_forks_repo_name: genkuroki/public