text:
################################################################################
# Julia C wrapper generator using libclang from the LLVM project               #
################################################################################

module wrap_c
    version = v"0.0.0"

using Clang.cindex

export wrap_c_headers
export WrapContext

### Reserved Julia identifiers to prepend with "_"
reserved_words = ["type", "begin", "end", "global", "local"]
function name_safe(c::CLCursor)
    cur_name = name(c)
    return (cur_name in reserved_words) ? "_"*cur_name : cur_name
end

### InternalOptions
type InternalOptions
    wrap_structs::Bool
end
InternalOptions() = InternalOptions(true)

### WrapContext
# stores shared information about the wrapping session
type WrapContext
    index::cindex.CXIndex
    output_file::ASCIIString
    common_file::ASCIIString
    clang_includes::Array{ASCIIString,1}         # clang include paths
    clang_args::Array{ASCIIString,1}             # additional {"-Arg", "value"} pairs for clang
    header_wrapped::Function                     # called to determine header inclusion status
    header_library::Function                     # called to determine shared library for given header
    header_outfile::Function                     # called to determine output file group for given header
    cursor_wrapped::Function                     # called to determine cursor inclusion statusk
    common_stream
    cache_wrapped::Set{ASCIIString}
    output_streams::Dict{ASCIIString, IO}
    options::InternalOptions
    anon_count::Int
end
WrapContext(idx,outfile,cmnfile,clanginc,clangextra,hwrap,hlib,hout,cwrap) =
    WrapContext(idx,outfile,cmnfile,convert(Array{ASCIIString,1},clanginc),convert(Array{ASCIIString,1}, clangextra),hwrap,hlib,hout,cwrap,
                None,Set{ASCIIString}(), Dict{ASCIIString,IO}(), InternalOptions(),0)

### Convenience function to initialize wrapping context with defaults
function init(;
            index                           = None,
            output_file::ASCIIString        = "",
            common_file::ASCIIString        = "",
            clang_args::Array{ASCIIString,1}
                                            = ASCIIString[],
            clang_includes::Array{ASCIIString,1}
                                            = ASCIIString[],
            clang_diagnostics::Bool         = true,
            header_wrapped                  = (header, cursorname) -> true,
            header_library                  = None,
            header_outputfile               = None,
            cursor_wrapped                  = (cursorname, cursor) -> true)

    # Set up some optional args if they are not explicitly passed.

    (index == None)         && ( index = cindex.idx_create(0, (clang_diagnostics ? 1 : 0)) )

    if (output_file == "" && header_outputfile == None)
        header_outputfile = x->joinpath(strip(splitext(basename(x))[1]) * ".jl")
    end

    (common_file == "")    && ( common_file = output_file )

    if (header_library == None)
        header_library = x->strip(splitext(basename(x))[1])
    elseif isa(header_library, ASCIIString)
        header_library = x->header_library
    end
    if (header_outputfile == None)
        header_outputfile = x->output_file
    end

    # Instantiate and return the WrapContext
    global context = WrapContext(index, output_file, common_file, clang_includes, clang_args, header_wrapped, header_library,header_outputfile,cursor_wrapped)
    return context
end

###############################################################################
#
# Mapping from libclang types to Julia types
#
#   This is primarily used with CXTypeKind enums (aka: cindex.TypeKind)
#   but can also be used to map arbitrary type names to corresponding
#   Julia entities, for example "size_t" -> :Csize_t
#
##############################################################################
cl_to_jl = {
    cindex.VoidType         => Void,
    cindex.BoolType         => Bool,
    cindex.Char_U           => Uint8,
    cindex.UChar            => :Cuchar,
    cindex.Char16           => Uint16,
    cindex.Char32           => Uint32,
    cindex.UShort           => Uint16,
    cindex.UInt             => Uint32,
    cindex.ULong            => :Culong,
    cindex.ULongLong        => :Culonglong,
    cindex.Char_S           => Uint8,
    cindex.SChar            => Uint8,
    cindex.WChar            => Char,
    cindex.Short            => Int16,
    cindex.IntType          => :Cint,
    cindex.Long             => :Clong,
    cindex.LongLong         => :Clonglong,
    cindex.Float            => :Cfloat,
    cindex.Double           => :Cdouble,
    cindex.LongDouble       => Float64,
    cindex.Enum             => :Cint,
    cindex.NullPtr          => C_NULL,
    cindex.UInt128          => Uint128,
    cindex.FirstBuiltin     => Void,
    "size_t"                => :Csize_t,
    "ptrdiff_t"             => :Cptrdiff_t
    }

################################################################################
#
# libclang objects to Julia representation
#
# each repr_jl function takes one or more CLCursor or CLType objects,
# and returns the appropriate string representation.
#
################################################################################

function repr_jl(t::Union(cindex.Record, cindex.Typedef))
    tname = spelling(cindex.getTypeDeclaration(t))
    return get(cl_to_jl, tname, tname)
end

function repr_jl(t::TypeRef)
    reftype = cindex.getCursorReferenced(t)
    refdef = cindex.getCursorDefinition(reftype)
    if isa(refdef, cindex.InvalidFile) ||
       isa(refdef, cindex.FirstInvalid) ||
       isa(refdef, cindex.Invalid)
        return "Void"
    else
        return spelling(reftype)
    end
end

function repr_jl(ptr::cindex.Pointer)
    ptee = pointee_type(ptr)
    return string("Ptr{", ((r = repr_jl(ptee)) == "" ? "Void" : r), "}")
end

function repr_jl(parm::cindex.ParmDecl)
    return repr_jl(cu_type(parm))
end

function repr_jl(unxp::Unexposed)
    return spelling(cindex.getTypeDeclaration(unxp))
end

function repr_jl(arg::cindex.CLType)
    return string(cl_to_jl[typeof(arg)])
end
function repr_jl(t::ConstantArray)
    # For ConstantArray declarations, we make an immutable
    # array with that many members of the appropriate type.

    # Override pointer representation so pointee type is not written
    repr_short(t::CLType) = isa(t, Pointer) ? "Ptr" : repr_jl(t)

    # Grab the WrapContext in order to add bitstype
    global context::WrapContext
    buf = context.common_stream

    arrsize = cindex.getArraySize(t)
    eltype = cindex.getArrayElementType(t)
    typename = string("Array_", arrsize, "_", repr_short(eltype))
    if !(typename in context.cache_wrapped)
        println(buf, "immutable ", typename)
        for i = 1:arrsize
            println(buf, "    d", i, "::", repr_jl(eltype))
        end
        println(buf, "end")
    end
    push!(context.cache_wrapped, typename)
    return typename
end

function repr_jl(t::UnionType)
    tdecl = cindex.getTypeDeclaration(t)
    maxelem = largestfield(tdecl)
    return repr_jl(maxelem)
end

###############################################################################
# Get field decl sizes
#   - used for hacky union inclusion: we find the largest union field and
#     declare a block of bytes to match.
###############################################################################

typesize(cu::CLCursor)     = typesize(cu_type(cu))
typesize(t::ConstantArray) = cindex.getArraySize(t)

function largestfield(cu::UnionDecl)
    maxsize,maxelem = 0
    fields = children(cu)
    for i in 1:length(fields)
        maxelem = ( (maxsize > (typesize(fields[i]))) ? maxelem : i )
    end
    fields[maxelem]
end

function fieldsize(cu::FieldDecl)
    fieldsize(children(cu)[1])
end

################################################################################
# Handle declarations
################################################################################

function wrap(buf::IO, cursor::EnumDecl; usename="")
    if (usename == "" && (usename = name(cursor)) == "")
        usename = name_anon()
    end
    enumname = usename
    println(buf, "# begin enum $enumname")
    println(buf, "typealias $enumname ", repr_jl(cindex.getEnumDeclIntegerType(cursor)))
    for enumitem in children(cursor)
        cur_name = cindex.spelling(enumitem)
        if (length(cur_name) < 1) continue end

        println(buf, "const ", cur_name, " = ", value(enumitem))
    end
    println(buf, "# end enum $enumname")
end

function wrap (buf::IO, sd::StructDecl; usename = "")
    if (usename == "" && (usename = name(sd)) == "")
        warn("Skipping unnamed StructDecl")
        return
    end
    global context::WrapContext
    if (!context.options.wrap_structs)
        return
    end

    prebuf = IOBuffer()
    outbuf = IOBuffer()

    # Generate type declaration
    println(outbuf, "type $usename")
    ccl = children(sd)
    if (length(ccl) < 1)
        warn("Skipping empty struct: \"$usename\"")
        return
    end
    for cu in children(sd)
        cur_name = spelling(cu)
        if (isa(cu, StructDecl) || isa(cu, UnionDecl))
            continue
        elseif !(isa(cu, FieldDecl) || isa(cu, TypeRef))
            warn("Skipping struct: \"$usename\" due to unsupported field: $cur_name")
            return
        elseif (length(cur_name) < 1)
            error("Unnamed struct member in: $usename ... cursor: ", string(cu))
        end
        if ((cur_name in reserved_words)) cur_name = "_"*cur_name end

        println(outbuf, "    ", cur_name, "::", repr_jl(cu_type(cu)))
        #if ((hlp = help_type(ty,cu)) != None)
        #    (hlp[1] != "") && println(prebuf, hlp[1])
        #    println(outbuf, "    ", cur_name, "::", hlp[2])
        #else
        #    println(outbuf, "    ", cur_name, "::", repr_jl(ty))
        #end
    end
    println(outbuf, "end")

    # print the type helpers to real output
    print(buf, takebuf_string(prebuf))
    # print the type declaration to real output
    print(buf, takebuf_string(outbuf))
end

function wrap(buf::IO, funcdecl::FunctionDecl, libname::ASCIIString)
    function print_args(buf::IO, cursors, types)
        i = 1
        for (c,t) in zip(cursors,types)
            print(buf, name_safe(c), "::", t)
            (i < length(cursors)) && print(buf, ", ")
            i += 1
        end
    end

    cu_spelling = spelling(funcdecl)

    funcname = spelling(funcdecl)
    arg_types = cindex.function_args(funcdecl)
    args = [x for x in search(funcdecl, ParmDecl)]
    arg_list = tuple( [repr_jl(x) for x in arg_types]... )
    ret_type = repr_jl(return_type(funcdecl))

    print(buf, "function ")
    print(buf, spelling(funcdecl))
    print(buf, "(")
    print_args(buf, args, [repr_jl(x) for x in arg_types])
    println(buf, ")")
    print(buf, "  ")
    print(buf, "ccall( (:", funcname, ", ", libname, "), ")
    print(buf, rep_type(ret_type))
    print(buf, ", ")
    print(buf, rep_args(arg_list), ", ")
    for (i,arg) in enumerate(args)
        print(buf, name_safe(arg))
        (i < length(args)) && print(buf, ", ")
    end
    println(buf, ")")
    println(buf, "end")
end

function wrap(buf::IO, tref::TypeRef; usename="")
    println("Wrap: ", tref)
end

function wrap(buf::IO, tdecl::TypedefDecl; usename="")
    function wrap_td(out::IO, t::CLType)
        println(buf, "typealias ",    spelling(t), " ", repr_jl(t))
    end

    cursor_type = cindex.cu_type(tdecl)
    td_type = cindex.getTypedefDeclUnderlyingType(tdecl)

    if isa(td_type, Unexposed)
        tdunxp = children(tdecl)[1]
        if isa(tdunxp, TypeRef)
            td_type = tdunxp
        else
            wrap(buf, tdunxp; usename=name(tdecl))
            return # TODO.. ugly flow
        end
    elseif isa(td_type, FunctionProto)
        return string("# Skipping Typedef: FunctionProto", spelling(tdecl))
    end

    println(buf, "typealias ",    spelling(tdecl), " ", repr_jl(td_type) )
end

################################################################################
# Handler for macro definitions
################################################################################
#
# For handling of #define'd constants, allows basic expressions
# but bails out quickly.

function lex_exprn(tokens::TokenList, pos::Int)
    function trans(tok)
        ops = ["+" "-" "*" ">>" "<<" "/" "\\" "%" "|" "||" "^" "&" "&&"]
        if (isa(tok, cindex.Literal) ||
            (isa(tok,cindex.Identifier))) return 0
        elseif (isa(tok, cindex.Punctuation) && tok.text in ops) return 1
        else return -1
        end
    end

    # normalize literal with a size suffix
    function literally(tok)
        # note: put multi-character first, or it will break out too soon for those!
        literalsuffixes = ["UL" "Ul" "uL" "ul" "LU" "Lu" "lU" "lu" "U" "u" "L" "l"]
        txt = tok.text
        if isa(tok,cindex.Identifier) || isa(tok,cindex.Punctuation)
            # pass
        elseif isa(tok,cindex.Literal)
            txt = strip(tok.text)
            for sfx in literalsuffixes
                if endswith(txt, sfx)
                    txt = txt[1:end-length(sfx)]
                    break
                end
            end
        end
        return txt
    end

    # check whether identifiers and literals alternate
    # with punctuation
    exprn = ""
    prev = 1 >> trans(tokens[pos])
    for pos = pos:length(tokens)
        tok = tokens[pos]
        state = trans(tok)
        if ( state $ prev  == 1)
            prev = state
        else
            break
        end
        exprn = exprn * literally(tok)
    end
    return (exprn,pos)
end

function wrap(strm::IO, md::cindex.MacroDefinition)
    tokens = tokenize(md)
    # Skip any empty definitions
    if(tokens.size < 2) return end
    if(beginswith(name(md), "_")) return end

    pos = 1; exprn = ""
    if(tokens[2].text == "(")
        exprn,pos = lex_exprn(tokens, 3)
        if (pos != endof(tokens) || tokens[pos].text != ")")
            print(strm, "# Skipping MacroDefinition: ", join([c.text for c in tokens]), "\n")
            return
        end
        exprn = "(" * exprn * ")"
    else
        (exprn,pos) = lex_exprn(tokens, 2)
    end
    print(strm, "const " * string(tokens[1].text) * " = " * exprn * "\n")
end

function wrap(buf::IO, cursor::TypeRef)
    println("Printing typeref: ", cursor)
    print(buf, name(cursor))
end

################################################################################
# Wrapping driver
################################################################################

function wrap_header(wc::WrapContext, topcu::CLCursor, top_hdr, ostrm::IO)
    println("WRAPPING HEADER: $top_hdr")

    topcl = children(topcu)

    # Loop over all of the child cursors and wrap them, if appropriate.
    for i=1:topcl.size
        cursor = topcl[i]
        cursor_hdr = cu_file(cursor)
        cursor_name = name(cursor)

        # Heuristic to decide what should be wrapped:
        #    1. always wrap things in the current top header (ie not includes)
        #    2. everything else is from includes, wrap if:
        #         - the client wants it wc.header_wrapped == True
        #         - the item has not already been wrapped (ie not in wc.cache_wrapped)
        if (cursor_hdr == top_hdr)
            # pass
        elseif (!wc.header_wrapped(top_hdr, cu_file(cursor)) ||
                        (cursor_name in wc.cache_wrapped) )
            continue
        elseif (!wc.cursor_wrapped(name(cursor), cursor))
            continue
        elseif (beginswith(cursor_name, "__"))
            # skip compiler definitions
            continue
        end

        if (isa(cursor, FunctionDecl))
            wrap(ostrm, cursor, wc.header_library(cu_file(cursor)))
        elseif (isa(cursor, EnumDecl))
            wrap(wc.common_stream, cursor)
        elseif (isa(cursor, TypedefDecl))
            wrap(wc.common_stream,cursor)
        elseif (isa(cursor, MacroDefinition))
            wrap(wc.common_stream,cursor)
        elseif (wc.options.wrap_structs && isa(cursor, StructDecl))
            wrap(wc.common_stream,cursor)
        else
            continue
        end
    end
    cindex.cl_dispose(topcl)
end

function header_output_stream(wc::WrapContext, hfile)
    jloutfile = wc.header_outfile(hfile)
    if (_x = get(wc.output_streams, jloutfile, None)) != None
        return _x
    else
        strm = IOStream("")
        try strm = open(jloutfile, "w")
        catch error("Unable to create output: $jloutfile for header: $hfile") end
        wc.output_streams[jloutfile] = strm
    end
    return strm
end

################################################################################
# Wrapping driver
################################################################################

function wrap_c_headers(wc::WrapContext, headers)

    println(wc.clang_includes)

    # Check headers!
    for h in headers
        if !isfile(h)
            error("Header file: ", h, " cannot be found")
        end
    end
    for d in wc.clang_includes
        if !isdir(d)
            error("Include file: ", d, " cannot be found")
        end
    end

    # Output stream for common items: typedefs, enums, etc.
    wc.common_stream = IOBuffer()

    # Generate the wrappings
    try
        for hfile in headers
            ostrm = header_output_stream(wc, hfile)
            println(ostrm, "# Julia wrapper for header: $hfile")
            println(ostrm, "# Automatically generated using Clang.jl wrap_c, version $version\n")

            topcu = cindex.parse_header(hfile;
                                 index = wc.index,
                                 args  = wc.clang_args,
                                 includes = wc.clang_includes,
                                 flags = TranslationUnit_Flags.DetailedPreprocessingRecord |
                                              TranslationUnit_Flags.SkipFunctionBodies)
            wrap_header(wc, topcu, hfile, ostrm)

            println(ostrm)
        end
    finally
        [close(os) for os in values(wc.output_streams)]
    end

    # Sort the common includes so that things aren't used out-of-order
    incl_lines = sort_common_includes(wc.common_stream)
    open(wc.common_file, "w") do strm
        [print(strm, l) for l in incl_lines]
    end

    close(wc.common_stream)
end

###############################################################################
# Utilities
###############################################################################

function sort_common_includes(strm::IOBuffer)
    # TODO: too many temporaries
    seek(strm,0)
    col1 = Dict{ASCIIString,Int}()
    col2 = Dict{ASCIIString,Int}()
    tmp = Dict{Int,ASCIIString}()
    pos = Int[]
    fnl = ASCIIString[]

    for (i,ln) in enumerate(readlines(strm))
        tmp[i] = ln
        if (m = match(r"@ctypedef (\w+) (?:Ptr{)?(\w+)(?:}?)", ln)) != nothing
            col1[m.captures[1]] = i
            col2[m.captures[2]] = i
        end
    end
    for s in sort(collect(keys(col2)))
        if( (m = get(col1, s, None))!=None)
            push!(fnl, tmp[m])
            push!(pos, m)
        end
    end

    kj = setdiff(1:length(keys(tmp)), pos)
    vcat(fnl,[tmp[i] for i in kj])
end

# eliminate symbol from the final type representation
function rep_type(t)
    replace(string(t), ":", "")
end

# Tuple representation without quotes
function rep_args(v)
    o = IOBuffer()
    print(o, "(")
    s = first = start(v)
    r = (!done(v,s))
    while r
        x,s = next(v,s)
        print(o, x)
        done(v,s) && break
        print(o, ", ")
    end
    r && (s == next(v,first)[2]) && print(o, ",")
    print(o, ")")
    seek(o, 0)
    readall(o)
end

function name_anon()
    global context::WrapContext
    "ANONYMOUS_"*string(context.anon_count += 1)
end

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

end # module wrap_c


text:
## PARTE 1
cd(@__DIR__())
using IterTools: groupby
keys = Set(("byr", "iyr", "eyr", "hgt", "hcl", "ecl", "pid"))
colors = Set(("amb", "blu", "brn", "gry", "grn", "hzl", "oth"))
rules = Dict(
    "byr" => x -> begin
        m = match(r"^([0-9][0-9][0-9][0-9])$", x)
        isnothing(m) && return false
        1920 <= parse(Int, m[1]) <= 2002
    end,
    "iyr" => x -> begin
        m = match(r"^([0-9][0-9][0-9][0-9])$", x)
        isnothing(m) && return false
        2010 <= parse(Int, m[1]) <= 2020
    end,
    "eyr" => x -> begin
        m = match(r"^([0-9][0-9][0-9][0-9])$", x)
        isnothing(m) && return false
        2020 <= parse(Int, m[1]) <= 2030
    end,
    "hgt" => x -> begin
        m = match(r"^([1-9][0-9]+)(cm|in)$", x)
        if isnothing(m)
            return false
        elseif m[2] == "cm"
            return 150 <= parse(Int, m[1]) <= 193
        elseif m[2] == "in"
            return 59 <= parse(Int, m[1]) <= 76
        end
    end,
    "hcl" => contains(r"^#[0-9a-f][0-9a-f][0-9a-f][0-9a-f][0-9a-f][0-9a-f]$"),
    "ecl" => x -> x ∈ colors,
    "pid" => contains(r"^[0-9][0-9][0-9][0-9][0-9][0-9][0-9][0-9][0-9]$"),
    "cid" => _ -> true,
)

##
sum(map(x->map(a->a[1:3],foldl(∪, split.(x))),Iterators.filter(x->x[1]!="", groupby(!isempty, eachline("input.txt"))))) do p
    keys ⊆ p
end

## PARTE 2
_validate_fields(p) = all(pairs(p)) do (k, f)
    rules[k](f)
end

##
res = filter(map(x->map(a->split(a,":"),foldl(∪, split.(x)))|>Dict,Iterators.filter(x->x[1]!="", groupby(!isempty, eachline("input.txt"))))) do p
    keys ⊆ Base.keys(p) && _validate_fields(p)
end

ext: jl


text:
function debugger(arena::Honeycomb,gstack::Array{UInt64,1},list::Array{Pointer,1},deb::Debugstate,ticks::Time)
    if deb.d==1
        println("\nticks:",ticks.t)
        @inbounds for ind=1:length(list)
            #println(ind,"  instr: ",arena.a[list[ind].loc[1],list[ind].loc[2]],"  dir: ",showdir(list[ind].dir),"  loc: ",list[ind].loc,"  ",transpose(list[ind].lstack),"•")
            println(ind,"  instr: ",arena.a[list[ind].loc[1],list[ind].loc[2]],"  dir: ",showdir(list[ind].dir),"  loc: ",list[ind].loc,"  ",transpose(convert(Array{BigInt},list[ind].lstack)),"•")
        end
        print("\n[")
        @inbounds for n=1:length(gstack)
            print(Int128(gstack[n]))
            n<length(gstack) ? print(" ") : nothing
        end
        println("]•")
        print("\n")
        println("———————")
    elseif deb.d >= 2
        println("\nticks: ",ticks.t)
        gdebug(arena,gstack,list,deb);print("\n")
        println("———————")
    end
end

function gdebug(arena::Honeycomb,gstack::Array{UInt64,1},list::Array{Pointer,1},deb::Debugstate)
    rows=maximum(length(arena.a[:,1]))
    cols=maximum(length(arena.a[1,:]))
    comb=deepcopy(arena.a)
    @inbounds for ind=1:length(list)
        comb[list[ind].loc[1],list[ind].loc[2]]=Char(ind+944)#translates IP numbers to greek letters
    end

    @inbounds for ind=1:length(list)
        #println(Char(ind+944),"  ",ind,"  instr: ",arena.a[list[ind].loc[1],list[ind].loc[2]]," wait: ",showwait(list[ind].wait)," dir: ",showdir(list[ind].dir),"  loc: ",list[ind].loc,"  ",transpose(list[ind].lstack),"•")
        println(Char(ind+944),"  ",ind,"  instr: ",arena.a[list[ind].loc[1],list[ind].loc[2]]," wait: ",showwait(list[ind].wait)," dir: ",showdir(list[ind].dir),"  loc: ",list[ind].loc,"  ",transpose(convert(Array{BigInt},list[ind].lstack)),"•")
    end

    println("\ngstack: ",transpose(convert(Array{BigInt},gstack)),"•\n\n")

    @inbounds for r=1:rows
        @inbounds for c=1:cols
            if Int(comb[r,c])>944 && Int(comb[r,c])<976 #greek letters unicode range
                if deb.d==3
                    c<cols ? print_with_color(:red,"$(comb[r,c])") : print_with_color(:red,"$(comb[r,c])\n")
                elseif deb.d==2
                    c<cols ? print(comb[r,c]) : println(comb[r,c])
                end
            else c<cols ? print(comb[r,c]) : println(comb[r,c])
            end
        end
    end
end

function showwait(wait::Int)
    wait==0 ? (return "∅") :
    wait==1 ? (return "∫") :
    wait==2 ? (return "∬") : nothing
end

function showdir(dir::Int)
    dir==0 ? (return "→") :
    dir==1 ? (return "↗") :
    dir==2 ? (return "↖") :
    dir==3 ? (return "←") :
    dir==4 ? (return "↙") :
    dir==5 ? (return "↘") : nothing
end


text:
module Coord

using Gadfly
using Compat
using Compose
using DataArrays
import Gadfly.Maybe
import Iterators: cycle

export cartesian

# Cartesian coordinates with position given by the x and y (and similar)
# aesthetics.
#
# Args:
#   xvars: Aesthetics to consider when choosing x bounds.
#   yvars: Aesthetics to consider when choosing y bounds.
#   xmin, xmax, ymin, ymax:
#     Force a particular x or y bound, rather than trying to choose it
#     from the data.
#
#
immutable Cartesian <: Gadfly.CoordinateElement
    xvars::Vector{Symbol}
    yvars::Vector{Symbol}
    xmin
    xmax
    ymin
    ymax
    xflip::Bool
    yflip::Bool
    fixed::Bool
    aspect_ratio::@compat(Union{(@compat Void), Float64})
    raster::Bool

    function Cartesian(
            xvars=[:x, :xmin, :xmax, :xintercept],
            yvars=[:y, :ymin, :ymax, :yintercept, :middle,
                   :lower_hinge, :upper_hinge, :lower_fence, :upper_fence, :outliers];
            xflip::Bool=false, yflip::Bool=false,
            xmin=nothing, xmax=nothing, ymin=nothing, ymax=nothing,
            fixed=false, aspect_ratio=nothing, raster=false)
        if isa(aspect_ratio, Real)
            aspect_ratio = convert(Float64, aspect_ratio)
        end
        new(xvars, yvars, xmin, xmax, ymin, ymax, xflip, yflip, fixed,
            aspect_ratio, raster)
    end
end

const cartesian = Cartesian

# Return the first concrete aesthetic value in one of the given aesthetics
#
# Args:
#   aess: An array of Aesthetics to search through.
#   vars: Aesthetic variables to search though.
#
# Returns:
#   A concrete value if one is found, otherwise nothing.
#
function first_concrete_aesthetic_value(aess::Vector{Gadfly.Aesthetics},
                                        vars::Vector{Symbol})
    T = aesthetics_type(aess, vars)
    for var in vars
        for aes in aess
            vals = getfield(aes, var)
            if vals === nothing
                continue
            end

            if !isa(vals, AbstractArray)
                vals = [vals]
            end

            if var == :outliers
                for outlier_vals in aes.outliers
                    for val in outlier_vals
                        if Gadfly.isconcrete(val)
                            return convert(T, val)
                        end
                    end
                end
                continue
            end

            for val in vals
                if Gadfly.isconcrete(val)
                    return convert(T, val)
                end
            end
        end
    end

    return nothing
end

# Find a common type among a group of aesthetics.
#
# Args:
#   aess: An array of Aesthetics to search through.
#   vars: Aesthetic variables to search though.
#
# Returns:
#   A common type.
function aesthetics_type(aess::Vector{Gadfly.Aesthetics},
                              vars::Vector{Symbol})
    T = @compat(Union{})
    for var in vars
        for aes in aess
            vals = getfield(aes, var)
            if vals === nothing
                continue
            end

            if var == :outliers
                if !isempty(vals)
                    T = promote_type(T, eltype(first(vals)))
                end
            else
                T = promote_type(T, eltype(vals))
            end
        end
    end

    return T
end

# Produce a context with suitable cartesian coordinates.
#
# Args:
#   coord: cartesian coordinate instance.
#
# Returns:
#   A compose Canvas.
#
function apply_coordinate(coord::Cartesian, aess::Vector{Gadfly.Aesthetics},
                          scales::Dict{Symbol, Gadfly.ScaleElement})
    pad_categorical_x = Nullable{Bool}()
    pad_categorical_y = Nullable{Bool}()
    for aes in aess
        if !isnull(aes.pad_categorical_x)
            if isnull(pad_categorical_x)
                pad_categorical_x = aes.pad_categorical_x
            else
                pad_categorical_x = Nullable(get(pad_categorical_x) || get(aes.pad_categorical_x))
            end
        end
        if !isnull(aes.pad_categorical_y)
            if isnull(pad_categorical_y)
                pad_categorical_y = aes.pad_categorical_y
            else
                pad_categorical_y = Nullable(get(pad_categorical_y) || get(aes.pad_categorical_y))
            end
        end
    end

    xmin = xmax = first_concrete_aesthetic_value(aess, coord.xvars)

    if xmin != nothing
        for var in coord.xvars
            for aes in aess
                vals = getfield(aes, var)
                if vals === nothing
                    continue
                end

                if !isa(vals, AbstractArray)
                    vals = [vals]
                end

                xmin, xmax = Gadfly.concrete_minmax(vals, xmin, xmax)
            end
        end
    end

    ymin = ymax = first_concrete_aesthetic_value(aess, coord.yvars)
    if ymin != nothing
        for var in coord.yvars
            for aes in aess
                vals = getfield(aes, var)
                if vals === nothing
                    continue
                end

                # Outliers is an odd aesthetic that needs special treatment.
                if var == :outliers
                    for outlier_vals in aes.outliers
                        ymin, ymax = Gadfly.concrete_minmax(outlier_vals, ymin, ymax)
                    end

                    continue
                end

                if !isa(vals, AbstractArray)
                    vals = [vals]
                end

                ymin, ymax = Gadfly.concrete_minmax(vals, ymin, ymax)
            end
        end
    end

    xviewmin = xviewmax = yviewmin = yviewmax = nothing

    # viewmin/max that is set explicitly should override min/max
    for aes in aess
        if aes.xviewmin != nothing
            xviewmin = xviewmin === nothing ? aes.xviewmin : min(xviewmin, aes.xviewmin)
        end

        if aes.xviewmax != nothing
            xviewmax = xviewmax === nothing ? aes.xviewmax : max(xviewmax, aes.xviewmax)
        end

        if aes.yviewmin != nothing
            yviewmin = yviewmin === nothing ? aes.yviewmin : min(yviewmin, aes.yviewmin)
        end

        if aes.yviewmax != nothing
            yviewmax = yviewmax === nothing ? aes.yviewmax : max(yviewmax, aes.yviewmax)
        end
    end

    xmax = xviewmax === nothing ? xmax : max(xmax, xviewmax)
    xmin = xviewmin === nothing ? xmin : min(xmin, xviewmin)
    ymax = yviewmax === nothing ? ymax : max(ymax, yviewmax)
    ymin = yviewmin === nothing ? ymin : min(ymin, yviewmin)

    # Hard limits set in Coord should override everything else
    xmin = coord.xmin === nothing ? xmin : coord.xmin
    xmax = coord.xmax === nothing ? xmax : coord.xmax
    ymin = coord.ymin === nothing ? ymin : coord.ymin
    ymax = coord.ymax === nothing ? ymax : coord.ymax

    if xmin === nothing || !isfinite(xmin)
        xmin = 0.0
        xmax = 1.0
    end

    if ymin === nothing || !isfinite(ymin)
        ymin = 0.0
        ymax = 1.0
    end

    xpadding = Scale.iscategorical(scales, :x) ? 0mm : 2mm
    ypadding = Scale.iscategorical(scales, :y) ? 0mm : 2mm

    if Scale.iscategorical(scales, :x) && (isnull(pad_categorical_x) || get(pad_categorical_x))
        xmin -= 0.5
        xmax += 0.5
    end

    if Scale.iscategorical(scales, :y) && (isnull(pad_categorical_y) || get(pad_categorical_y))
        ymin -= 0.5
        ymax += 0.5
    end

    width  = xmax - xmin
    height = ymax - ymin

    compose!(
        context(units=UnitBox(
            coord.xflip ? xmax : xmin,
            coord.yflip ? ymin : ymax,
            coord.xflip ? -width : width,
            coord.yflip ? height : -height,
            leftpad=xpadding,
            rightpad=xpadding,
            toppad=ypadding,
            bottompad=ypadding),
            raster=coord.raster),
        svgclass("plotpanel"))
end

immutable SubplotGrid <: Gadfly.CoordinateElement
end

function apply_coordinate(coord::SubplotGrid, aess::Vector{Gadfly.Aesthetics},
                          scales::Dict{Symbol, Gadfly.ScaleElement})
    compose!(context(), svgclass("plotpanel"))
end

const subplot_grid = SubplotGrid

end # module Coord


text:
using LinearAlgebra
using LightGraphs
using MetaGraphs
using NPZ
using Plots
using SpinGlassPEPS
using Logging
using ArgParse
using CSV
using Test

disable_logging(LogLevel(0))
δH = 0.9
# this is axiliary function for npz write

function vecvec2matrix(v::Vector{Vector{Int}})
    M = v[1]
    for i in 2:length(v)
        M = hcat(M, v[i])
    end
    M
end

s = ArgParseSettings("description")
  @add_arg_table! s begin
    "--file", "-i"
    arg_type = String
    help = "the file name"
    "--beta", "-b"
    default = 2.
    arg_type = Float64
    help = "beta value"
    "--chi", "-c"
    default = typemax(Int)
    arg_type = Int
    help = "cutoff size"
    "--n_sols", "-n"
    default = 10
    arg_type = Int
    help = "maximal number of solutions"
    "--node_size1", "-r"
    default = 1
    arg_type = Int
    help = "chimera node size in rows"
    "--node_size2", "-o"
    default = 1
    arg_type = Int
    help = "chimera node size in columns"
    "--spectrum_cutoff", "-u"
    default = 256
    arg_type = Int
    help = "maximal size of the lower spectrum"
    "--lower_cutoff", "-v"
    default = 1
    arg_type = Int
    help = "minimal size of the lower spectrum"
  end

fi = parse_args(s)["file"]
file = split(fi, "/")[end]
folder = fi[1:end-length(file)]
println(file)
println(folder)
ig = ising_graph(fi)

β = parse_args(s)["beta"]
χ = parse_args(s)["chi"]

spectrum_cutoff = parse_args(s)["spectrum_cutoff"]
lower_cutoff = parse_args(s)["lower_cutoff"]

n_sols = parse_args(s)["n_sols"]
node_size = (parse_args(s)["node_size1"], parse_args(s)["node_size2"])
println(node_size)

s1 = isqrt(div( nv(ig), 8))
n = div(s1, node_size[1])
m = div(s1, node_size[2])

fil = folder*"groundstates_otn2d.txt"
data = split.(readlines(open(fil)))

i = findall(x->x[1]==file, data)[1]
ground_ref = [parse(Int, el) for el in data[i][4:end]]

f(i) = (i == 0)  ? -1 : 1
ground_spins = map(f, ground_ref)
energy_ref = energy(ground_spins, ig)

ses = spectrum_cutoff:-1:lower_cutoff
step = 10
n_s = n_sols:-step:1

delta_e = ones(length(ses), length(n_s))
cut = ones(Int, length(ses), length(n_s))
rule = super_square_lattice((m, node_size[1], n, node_size[2], 8))

function proceed()
  j = 1
  for n_sol in n_s
    i = 1
    for sc in ses
      println(sc)
      D = Dict{Int, Int}()
      for v in vertices(ig)
        push!(D, (v => sc))
      end
      fg = factor_graph(
            ig,
            D,
            cluster_assignment_rule=rule,
            spectrum=brute_force,
        )

      println("..............")
      for v in vertices(fg)
        println(v, " ", length(props(fg, v)[:spectrum].energies))
      end
      println("..............")

      control_params = Dict(
           "bond_dim" => χ,
           "var_tol" => 1E-12,
           "sweeps" => 4.
       )

      peps = PEPSNetwork(m, n, fg, β, :NW, control_params)

      @time sols = low_energy_spectrum(peps, n_sol)

      en = minimum(sols.energies)

      cut[i,j] = sc
      delta_e[i,j] = (en-energy_ref)/abs(energy_ref)
      i = i+1

      println("spectrum cutoff = ", sc)
      println("no sols = ", n_sol)
      println("percentage E = ", (en-energy_ref)/abs(energy_ref))
    end
    j = j+1
  end

  d = Dict("percentage_delta_e" => delta_e, "spectrum_size" => cut, "chi" => χ, "beta" => β, "chimera_size" => nv(ig), "no_solutions" => n_s)
  npzwrite(folder*"low_energy_sols/"*file[1:3]*"approx$(nv(ig)).npz", d)

  p = plot(cut[:,1], delta_e[:,1], title = "low spectr. approx., beta = $β, chi = $χ, chimera_s. = $(nv(ig))", label = "n_sols = $(n_s[1])", lw = 1.5, left_margin = 10Plots.mm, bottom_margin = 10Plots.mm)
  xlabel!("size of the low energy spectrum")
  ylabel!("percentage delta E")
  for j in 2:size(cut,2)
    plot!(p, cut[:,j], delta_e[:,j], label = "n_sols = $(n_s[j])", lw = 1.5)
  end
  savefig(folder*"low_energy_sols/"*file[1:3]*"approx$(nv(ig)).pdf")
end

proceed()


text:
using GridapGPU
using Test

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


text:
# CUDA module management

immutable CuModule
    handle::Ptr{Void}

    function CuModule(filename::ASCIIString)
        a = Array(Ptr{Void}, 1)
        @cucall(cuModuleLoad, (Ptr{Ptr{Void}}, Ptr{Cchar}), a, filename)
        new(a[1])
    end
end

function unload(md::CuModule)
    @cucall(cuModuleUnload, (Ptr{Void},), md.handle)
end

immutable CuFunction
    handle::Ptr{Void}

    function CuFunction(md::CuModule, name::ASCIIString)
        a = Array(Ptr{Void}, 1)
        @cucall(cuModuleGetFunction, (Ptr{Ptr{Void}}, Ptr{Void}, Ptr{Cchar}),
            a, md.handle, name)
        new(a[1])
    end
end


text:
using GenomicFeatures
using Test
using Distributions
import Random
import BGZFStreams
import YAML
import ColorTypes: RGB
import BioSequences: @dna_str, FASTA
import FixedPointNumbers: N0f8
import BioCore.Exceptions: MissingFieldException

# Test that an array of intervals is well ordered
function is_all_ordered(intervals::Vector{I}) where I <: Interval
    for i = 2:length(intervals)
        if !GenomicFeatures.isordered(intervals[i-1], intervals[i])
            return false
        end
    end
    return true
end

# Generate an array of n random Interval{Int} object. With sequence names
# samples from seqnames, and intervals drawn to lie in [1, maxpos].
function random_intervals(seqnames, maxpos::Int, n::Int)
    seq_dist = Categorical(length(seqnames))
    strand_dist = Categorical(2)
    length_dist = Normal(1000, 1000)
    intervals = Vector{Interval{Int}}(undef, n)
    for i in 1:n
        intlen = maxpos
        while intlen >= maxpos || intlen <= 0
            intlen = ceil(Int, rand(length_dist))
        end
        first = rand(1:maxpos-intlen)
        last = first + intlen - 1
        strand = rand(strand_dist) == 1 ? STRAND_POS : STRAND_NEG
        intervals[i] = Interval{Int}(seqnames[rand(seq_dist)],
                                     first, last, strand, i)
    end
    return intervals
end

# A simple interval intersection implementation to test against.
function simple_intersection(intervals_a, intervals_b; filter=(a,b)->true)
    sort!(intervals_a)
    sort!(intervals_b)
    intersections = Any[]
    i = 1
    j = 1
    while i <= length(intervals_a) && j <= length(intervals_b)
        ai = intervals_a[i]
        bj = intervals_b[j]
        if isless(ai.seqname, bj.seqname) ||
           (ai.seqname == bj.seqname && ai.last < bj.first)
            i += 1
        elseif isless(bj.seqname, ai.seqname) ||
               (ai.seqname == bj.seqname && bj.last < ai.first)
            j += 1
        else
            k = j
            while k <= length(intervals_b) && intervals_b[k].first <= ai.last
                if isoverlapping(ai, intervals_b[k]) && filter(ai, intervals_b[k])
                    push!(intersections, (ai, intervals_b[k]))
                end
                k += 1
            end
            i += 1
        end
    end
    return intersections
end

function simple_coverage(intervals)
    seqlens = Dict{String, Int}()
    for interval in intervals
        if get(seqlens, interval.seqname, -1) < interval.last
            seqlens[interval.seqname] = interval.last
        end
    end

    covarrays = Dict{String, Vector{Int}}()
    for (seqname, seqlen) in seqlens
        covarrays[seqname] = zeros(Int, seqlen)
    end

    for interval in intervals
        arr = covarrays[interval.seqname]
        for i in interval.first:interval.last
            arr[i] += 1
        end
    end

    covintervals = Interval{UInt32}[]
    for (seqname, arr) in covarrays
        i = j = 1
        while i <= length(arr)
            if arr[i] > 0
                j = i + 1
                while j <= length(arr) && arr[j] == arr[i]
                    j += 1
                end
                push!(covintervals,
                      Interval{UInt32}(seqname, i, j - 1, STRAND_BOTH, arr[i]))
                i = j
            else
                i += 1
            end
        end
    end

    return covintervals
end

function get_bio_fmt_specimens(commit="222f58c8ef3e3480f26515d99d3784b8cfcca046")
    path = joinpath(dirname(@__FILE__), "BioFmtSpecimens")
    if !isdir(path)
        run(`git clone https://github.com/BioJulia/BioFmtSpecimens.git $(path)`)
    end
    cd(path) do
        #run(`git checkout $(commit)`)
    end
    return path
end

@testset "Strand" begin
    @testset "Constructor" begin
        @test Strand('?') === STRAND_NA
        @test Strand('+') === STRAND_POS
        @test Strand('-') === STRAND_NEG
        @test Strand('.') === STRAND_BOTH
        @test_throws Exception Strand('x')
    end

    @testset "Conversion" begin
        @test convert(Strand, '?') === STRAND_NA
        @test convert(Strand, '+') === STRAND_POS
        @test convert(Strand, '-') === STRAND_NEG
        @test convert(Strand, '.') === STRAND_BOTH

        @test convert(Char, STRAND_NA) === '?'
        @test convert(Char, STRAND_POS) === '+'
        @test convert(Char, STRAND_NEG) === '-'
        @test convert(Char, STRAND_BOTH) === '.'
    end

    @testset "Order" begin
        @test STRAND_NA < STRAND_POS < STRAND_NEG < STRAND_BOTH
    end

    @testset "IO" begin
        # show
        buf = IOBuffer()
        for s in [STRAND_NA, STRAND_POS, STRAND_NEG, STRAND_BOTH]
            show(buf, s); print(buf, " ")
        end
        @test String(take!(buf)) == "STRAND_NA STRAND_POS STRAND_NEG STRAND_BOTH "

        # print
        buf = IOBuffer()
        for s in [STRAND_NA, STRAND_POS, STRAND_NEG, STRAND_BOTH]
            print(buf, s)
        end
        @test String(take!(buf)) == "?+-."

        # read and write
        buf = IOBuffer()
        @test write(buf, STRAND_NA)   == 1
        @test write(buf, STRAND_POS)  == 1
        @test write(buf, STRAND_NEG)  == 1
        @test write(buf, STRAND_BOTH) == 1
        seekstart(buf)
        @test read(buf, Strand) === STRAND_NA
        @test read(buf, Strand) === STRAND_POS
        @test read(buf, Strand) === STRAND_NEG
        @test read(buf, Strand) === STRAND_BOTH
        @test eof(buf)
    end
end

@testset "Interval" begin
    @testset "Constructor" begin
        i = Interval("chr1", 10, 20)
        @test seqname(i) == "chr1"
        @test leftposition(i) == 10
        @test rightposition(i) == 20
        @test strand(i) == STRAND_BOTH
        @test i == Interval("chr1", 10:20)

        i1 = Interval("chr1", 10, 20, '+')
        i2 = Interval("chr1", 10, 20, STRAND_POS)
        @test i1 == i2
        @test i1 == Interval("chr1", 10:20, '+')

        i1 = Interval("chr2", 5692667, 5701385, '+',        "SOX11")
        i2 = Interval("chr2", 5692667, 5701385, STRAND_POS, "SOX11")
        @test i1 == i2
        @test i1 == Interval("chr2", 5692667:5701385, '+', "SOX11")
    end
end

@testset "IntervalCollection" begin
    @testset "Insertion/Iteration" begin
        n = 100000
        intervals = random_intervals(["one", "two", "three"], 1000000, n)
        ic = IntervalCollection{Int}()

        @test isempty(ic)
        @test collect(Interval{Int}, ic) == Interval{Int}[]

        for interval in intervals
            push!(ic, interval)
        end
        @test is_all_ordered(collect(Interval{Int}, ic))
    end

    @testset "Intersection" begin
        n = 1000
        Random.seed!(1234)
        intervals_a = random_intervals(["one", "two", "three"], 1000000, n)
        intervals_b = random_intervals(["one", "three", "four"], 1000000, n)

        # empty versus empty
        ic_a = IntervalCollection{Int}()
        ic_b = IntervalCollection{Int}()
        @test collect(eachoverlap(ic_a, ic_b)) == Any[]

        # empty versus non-empty
        for interval in intervals_a
            push!(ic_a, interval)
        end

        @test collect(eachoverlap(ic_a, ic_b)) == Any[]
        @test collect(eachoverlap(ic_b, ic_a)) == Any[]

        # non-empty versus non-empty
        for interval in intervals_b
            push!(ic_b, interval)
        end

        @test sort(collect(eachoverlap(ic_a, ic_b))) ==
              sort(simple_intersection(intervals_a, intervals_b))
        @test sort(collect(eachoverlap(ic_a, ic_b, filter=(a,b) -> isodd(first(a))))) ==
              sort(simple_intersection(intervals_a, intervals_b, filter=(a,b) -> isodd(first(a))))
    end

    @testset "Show" begin
        ic = IntervalCollection{Int}()
        show(devnull, ic)

        push!(ic, Interval{Int}("one", 1, 1000, STRAND_POS, 0))
        show(devnull, ic)

        intervals = random_intervals(["one", "two", "three"], 1000000, 100)
        for interval in intervals
            push!(ic, interval)
        end
        show(devnull, ic)

        show(devnull, STRAND_NA)
        show(devnull, STRAND_POS)
        show(devnull, STRAND_NEG)
        show(devnull, STRAND_BOTH)
    end
end

@testset "IntervalStream" begin
    @testset "Intersection" begin
        n = 1000
        Random.seed!(1234)
        intervals_a = random_intervals(["one", "two", "three"], 1000000, n)
        intervals_b = random_intervals(["one", "three", "four"], 1000000, n)

        ic_a = IntervalCollection{Int}()
        ic_b = IntervalCollection{Int}()

        for interval in intervals_a
            push!(ic_a, interval)
        end

        for interval in intervals_b
            push!(ic_b, interval)
        end

        # non-empty versus non-empty, stream intersection
        it = eachoverlap(ic_a, ic_b, isless)
        @test sort(collect(it)) == sort(simple_intersection(intervals_a, intervals_b))

        it = eachoverlap(ic_a, ic_b, isless, filter=(a,b) -> isodd(first(a)))
        @test sort(collect(it)) ==
              sort(simple_intersection(intervals_a, intervals_b, filter=(a,b) -> isodd(first(a))))

        it = eachoverlap(
            [Interval("a", 1, 100, STRAND_POS, nothing), Interval("c", 1, 100, STRAND_POS, nothing)],
            [Interval("a", 1, 100, STRAND_POS, nothing), Interval("b", 1, 100, STRAND_POS, nothing)],
            isless)
        @test length(collect(it)) == 1

        it = eachoverlap(
            [Interval("c", 1, 100, STRAND_POS, nothing), Interval("d", 1, 100, STRAND_POS, nothing)],
            [Interval("b", 1, 100, STRAND_POS, nothing), Interval("d", 1, 100, STRAND_POS, nothing)],
            isless)
        @test length(collect(it)) == 1

        # unsorted streams are not allowed
        @test_throws Exception begin
            it = eachoverlap(
                [Interval("b", 1, 1000, STRAND_POS, nothing),
                 Interval("a", 1, 1000, STRAND_POS, nothing)],
                [Interval("a", 1, 1000, STRAND_POS, nothing),
                 Interval("b", 1, 1000, STRAND_POS, nothing)], isless)
            collect(it)
        end

        @test_throws Exception begin
            it = eachoverlap(
                [Interval("a", 1, 1000, STRAND_POS, nothing),
                 Interval("a", 500, 1000, STRAND_POS, nothing),
                 Interval("a", 400, 2000, STRAND_POS, nothing)],
                [Interval("a", 1, 1000, STRAND_POS, nothing),
                 Interval("b", 1, 1000, STRAND_POS, nothing)], isless)
            collect(it)
        end
    end

    @testset "IntervalStream Intersection" begin
        n = 1000
        Random.seed!(1234)
        intervals_a = random_intervals(["one", "two", "three"], 1000000, n)
        intervals_b = random_intervals(["one", "two", "three"], 1000000, n)

        ic_a = IntervalCollection{Int}()
        ic_b = IntervalCollection{Int}()

        for interval in intervals_a
            push!(ic_a, interval)
        end

        for interval in intervals_b
            push!(ic_b, interval)
        end

        @test sort(collect(eachoverlap(ic_a, ic_b, isless))) ==
              sort(simple_intersection(intervals_a, intervals_b))
        @test sort(collect(eachoverlap(ic_a, ic_b, isless, filter=(a,b) -> isodd(first(a))))) ==
              sort(simple_intersection(intervals_a, intervals_b, filter=(a,b) -> isodd(first(a))))
    end

    @testset "IntervalStream Coverage" begin
        n = 10000
        Random.seed!(1234)
        intervals = random_intervals(["one", "two", "three"], 1000000, n)

        ic = IntervalCollection{Int}()
        for interval in intervals
            push!(ic, interval)
        end

        @test sort(simple_coverage(intervals)) == sort(collect(coverage(ic)))
    end

    @testset "eachoverlap" begin
        i = Interval("chr1", 1, 10)
        intervals_a = typeof(i)[]
        @test length(collect(eachoverlap(intervals_a, intervals_a))) == 0

        intervals_a = [Interval("chr1", 1, 10)]
        intervals_b = eltype(intervals_a)[]
        @test length(collect(eachoverlap(intervals_a, intervals_b))) == 0
        @test length(collect(eachoverlap(intervals_b, intervals_a))) == 0

        intervals_a = [Interval("chr1", 1, 10)]
        @test length(collect(eachoverlap(intervals_a, intervals_a))) == 1

        intervals_a = [Interval("chr1", 1, 10)]
        intervals_b = [Interval("chr2", 1, 10)]
        @test length(collect(eachoverlap(intervals_a, intervals_b))) == 0
        @test length(collect(eachoverlap(intervals_b, intervals_a))) == 0

        intervals_a = [Interval("chr1", 1, 10)]
        intervals_b = [Interval("chr1", 11, 15)]
        @test length(collect(eachoverlap(intervals_a, intervals_b))) == 0
        @test length(collect(eachoverlap(intervals_b, intervals_a))) == 0

        intervals_a = [Interval("chr1", 11, 15)]
        intervals_b = [Interval("chr1", 1, 10), Interval("chr1", 12, 13)]
        @test length(collect(eachoverlap(intervals_a, intervals_b))) == 1
        @test length(collect(eachoverlap(intervals_b, intervals_a))) == 1

        intervals_a = [Interval("chr1", 1, 2), Interval("chr1", 2, 5), Interval("chr2", 1, 10)]
        intervals_b = [Interval("chr1", 1, 2), Interval("chr1", 2, 3), Interval("chr2", 1, 2)]
        @test length(collect(eachoverlap(intervals_a, intervals_b))) == 5
        @test length(collect(eachoverlap(intervals_b, intervals_a))) == 5

        intervals_a = [Interval("chr1", 1, 2), Interval("chr1", 3, 5), Interval("chr2", 1, 10)]
        @test length(collect(eachoverlap(intervals_a, intervals_a))) == 3

        # compare generic and specific eachoverlap methods
        intervals_a = [Interval("chr1", 1, 2), Interval("chr1", 1, 3), Interval("chr1", 5, 9),
                       Interval("chr2", 1, 5), Interval("chr2", 6, 6), Interval("chr2", 6, 8)]
        intervals_b = intervals_a
        ic_a = IntervalCollection(intervals_a)
        ic_b = IntervalCollection(intervals_b)
        iter1 = eachoverlap(intervals_a, intervals_b)
        iter2 = eachoverlap(intervals_a, ic_b)
        iter3 = eachoverlap(ic_a, intervals_b)
        iter4 = eachoverlap(ic_a, ic_b)
        @test collect(iter1) == collect(iter2) == collect(iter3) == collect(iter4)
    end
end

@testset "Indexes" begin
    # TODO
    @test GenomicFeatures.Indexes.Tabix === GenomicFeatures.Indexes.Tabix
end

@testset "BED" begin
    @testset "Record" begin
        record = BED.Record("chr1\t17368\t17436")
        @test BED.chrom(record) == "chr1"
        @test BED.chromstart(record) === 17369
        @test BED.chromend(record) === 17436
        @test !BED.hasname(record)
        @test !BED.hasscore(record)
        @test !BED.hasstrand(record)
        @test !BED.hasthickstart(record)
        @test !BED.hasthickend(record)
        @test !BED.hasitemrgb(record)
        @test !BED.hasblockcount(record)
        @test !BED.hasblocksizes(record)
        @test !BED.hasblockstarts(record)
        @test startswith(repr(record), "GenomicFeatures.BED.Record:\n")
        @test string(record) == "chr1\t17368\t17436"

        record = BED.Record("chrXIII\t854794\t855293\tYMR292W\t0\t+\t854794\t855293\t0\t2\t22,395,\t0,104,")
        @test BED.chrom(record) == "chrXIII"
        @test BED.chromstart(record) === 854795
        @test BED.chromend(record) === 855293
        @test BED.name(record) == "YMR292W"
        @test BED.score(record) === 0
        @test BED.strand(record) === STRAND_POS
        @test BED.thickstart(record) === 854795
        @test BED.thickend(record) === 855293
        @test BED.itemrgb(record) === RGB(0, 0, 0)
        @test BED.blockcount(record) === 2
        @test BED.blocksizes(record) == [22, 395]
        @test BED.blockstarts(record) == [1, 105]
        @test startswith(repr(record), "GenomicFeatures.BED.Record:\n")
        @test string(record) == "chrXIII\t854794\t855293\tYMR292W\t0\t+\t854794\t855293\t0\t2\t22,395,\t0,104,"

        record = BED.Record("chrX\t151080532\t151081699\tCHOCOLATE1\t0\t-\t151080532\t151081699\t255,127,36")
        @test BED.chrom(record) == "chrX"
        @test BED.chromstart(record) === 151080533
        @test BED.chromend(record) === 151081699
        @test BED.name(record) == "CHOCOLATE1"
        @test BED.score(record) === 0
        @test BED.strand(record) === STRAND_NEG
        @test BED.thickstart(record) === 151080533
        @test BED.thickend(record) === 151081699
        @test BED.itemrgb(record) === RGB(map(x -> reinterpret(N0f8, UInt8(x)), (255, 127, 36))...)
        @test !BED.hasblockcount(record)
        @test !BED.hasblocksizes(record)
        @test !BED.hasblockstarts(record)
        @test startswith(repr(record), "GenomicFeatures.BED.Record:\n")
        @test string(record) == "chrX\t151080532\t151081699\tCHOCOLATE1\t0\t-\t151080532\t151081699\t255,127,36"
    end

    function check_bed_parse(filename)
        # Reading from a stream
        for interval in BED.Reader(open(filename))
        end

        # Reading from a regular file
        for interval in open(BED.Reader, filename)
        end

        # in-place parsing
        stream = open(BED.Reader, filename)
        entry = eltype(stream)()
        while !eof(stream)
            read!(stream, entry)
        end
        close(stream)

        # Check round trip
        output = IOBuffer()
        writer = BED.Writer(output)
        expected_entries = BED.Record[]
        for interval in open(BED.Reader, filename)
            write(writer, interval)
            push!(expected_entries, interval)
        end
        flush(writer)

        seekstart(output)
        read_entries = BED.Record[]
        for interval in BED.Reader(output)
            push!(read_entries, interval)
        end

        return expected_entries == read_entries
    end

    path = joinpath(get_bio_fmt_specimens(), "BED")
    for specimen in YAML.load_file(joinpath(path, "index.yml"))
        if "gzip" ∈ split(get(specimen, "tags", ""))
            # skip compressed files
            continue
        end
        valid = get(specimen, "valid", true)
        if valid
            @test check_bed_parse(joinpath(path, specimen["filename"]))
        else
            @test_throws Exception check_bed_parse(joinpath(path, specimen["filename"]))
        end
    end

    # Testing strategy: there are two entirely separate intersection
    # algorithms for IntervalCollection and IntervalStream. Here we test
    # them both by checking that they agree by generating and intersecting
    # random BED files.

    function check_intersection(filename_a, filename_b)
        ic_a = IntervalCollection{BED.Record}()
        open(BED.Reader, filename_a) do reader
            for record in reader
                push!(ic_a, Interval(record))
            end
        end

        ic_b = IntervalCollection{BED.Record}()
        open(BED.Reader, filename_b) do reader
            for record in reader
                push!(ic_b, Interval(record))
            end
        end

        # This is refactored out to close streams
        fa = open(BED.Reader, filename_a)
        fb = open(BED.Reader, filename_b)
        xs = sort(collect(eachoverlap(fa, fb)))
        close(fa)
        close(fb)

        ys = sort(collect(eachoverlap(ic_a, ic_b)))

        return xs == ys
    end

    function write_intervals(filename, intervals)
        open(filename, "w") do out
            for interval in sort(intervals)
                println(out, interval.seqname, "\t", interval.first - 1,
                        "\t", interval.last, "\t", interval.metadata, "\t",
                        1000, "\t", interval.strand)
            end
        end

    end

    n = 10000
    Random.seed!(1234)
    intervals_a = random_intervals(["one", "two", "three", "four", "five"],
                                   1000000, n)
    intervals_b = random_intervals(["one", "two", "three", "four", "five"],
                                   1000000, n)

    mktempdir() do dir
        filename_a = joinpath(dir, "test_a.bed")
        filename_b = joinpath(dir, "test_b.bed")
        write_intervals(filename_a, intervals_a)
        write_intervals(filename_b, intervals_b)
        @test check_intersection(filename_a, filename_b)
    end

    @testset "eachoverlap" begin
        path = joinpath(get_bio_fmt_specimens(), "BED", "ws245Genes.WBGene.bed.bgz")
        stream = BGZFStreams.BGZFStream(path)
        reader = BED.Reader(stream, index=string(path, ".tbi"))
        for (interval, n_records) in [
                (Interval("chrII", 500_000:10_000_000), 38),
                (Interval("chrII", 1:70_000), 0),
                (Interval("chrIV", 500_000:10_000_000), 43),
                (Interval("chrI",  1_000_000:12_000_000), 27),]
            n = 0
            for record in eachoverlap(reader, interval)
                n += 1
            end
            @test n == n_records
        end
        @test isa(BED.Reader(path), BED.Reader)
    end
end

@testset "GFF3" begin
    record = GFF3.Record()
    @test !isfilled(record)
    @test repr(record) == "GenomicFeatures.GFF3.Record: <not filled>"

    record = GFF3.Record("CCDS1.1\tCCDS\tgene\t801943\t802434\t.\t-\t.\tNAME=LINC00115")
    @test isfilled(record)
    @test GFF3.isfeature(record)
    @test hasseqname(record)
    @test GFF3.hasseqid(record)
    @test seqname(record) == GFF3.seqid(record) == "CCDS1.1"
    @test GFF3.hassource(record)
    @test GFF3.source(record) == "CCDS"
    @test GFF3.hasfeaturetype(record)
    @test GFF3.featuretype(record) == "gene"
    @test GFF3.hasseqstart(record) === hasleftposition(record) === true
    @test GFF3.seqstart(record) === leftposition(record) === 801943
    @test GFF3.hasseqend(record) === hasrightposition(record) === true
    @test GFF3.seqend(record) === rightposition(record) === 802434
    @test !GFF3.hasscore(record)
    @test_throws MissingFieldException GFF3.score(record)
    @test GFF3.hasstrand(record)
    @test strand(record) === GFF3.strand(record) === STRAND_NEG
    @test !GFF3.hasphase(record)
    @test_throws MissingFieldException GFF3.phase(record)
    @test GFF3.attributes(record) == ["NAME" => ["LINC00115"]]
    @test GFF3.attributes(record, "NAME") == ["LINC00115"]
    @test GFF3.content(record) == "CCDS1.1\tCCDS\tgene\t801943\t802434\t.\t-\t.\tNAME=LINC00115"
    @test startswith(repr(record), "GenomicFeatures.GFF3.Record:\n")
    @test string(record) == "CCDS1.1\tCCDS\tgene\t801943\t802434\t.\t-\t.\tNAME=LINC00115"

    record = GFF3.Record("##gff-version 3")
    @test isfilled(record)
    @test GFF3.isdirective(record)
    @test GFF3.content(record) == "gff-version 3"
    @test convert(String, record) == "##gff-version 3"

    record = GFF3.Record("#comment")
    @test isfilled(record)
    @test GFF3.iscomment(record)
    @test GFF3.content(record) == "comment"
    @test convert(String, record) == "#comment"

    function check_gff3_parse(filename)
        # Reading from a stream
        num_intervals = 0
        for interval in GFF3.Reader(open(filename))
            num_intervals += 1
        end

        # Reading from a regular file
        for interval in open(GFF3.Reader, filename)
        end

        collection = IntervalCollection(open(GFF3.Reader, filename))
        @test length(collection) == num_intervals

        # in-place parsing
        stream = open(GFF3.Reader, filename)
        entry = eltype(stream)()
        while !eof(stream)
            try
                read!(stream, entry)
            catch ex
                if isa(ex, EOFError)
                    break
                end
            end
        end
        close(stream)

        # copy
        records = GFF3.Record[]
        reader = open(GFF3.Reader, filename)
        output = IOBuffer()
        writer = GFF3.Writer(output)
        for record in reader
            write(writer, record)
            push!(records, record)
        end
        close(reader)
        flush(writer)

        records2 = GFF3.Record[]
        for record in GFF3.Reader(IOBuffer(take!(output)))
            push!(records2, record)
        end
        return records == records2
    end

    path = joinpath(get_bio_fmt_specimens(), "GFF3")
    for specimen in YAML.load_file(joinpath(path, "index.yml"))
        if "gzip" ∈ split(get(specimen, "tags", ""))
            # skip compressed files
            continue
        end
        valid = get(specimen, "valid", true)
        if valid
            @test check_gff3_parse(joinpath(path, specimen["filename"]))
        else
            @test_throws Exception check_gff3_parse(joinpath(path, specimen["filename"]))
        end
    end

    # no fasta
    test_input = """
1	havana	exon	870086	870201	.	-	.	Parent=transcript:ENST00000432963;Name=ENSE00001791782;constitutive=0;ensembl_end_phase=-1;ensembl_phase=-1;exon_id=ENSE00001791782;rank=1;version=1
1	havana	lincRNA	868403	876802	.	-	.	ID=transcript:ENST00000427857;Parent=gene:ENSG00000230368;Name=FAM41C-002;biotype=lincRNA;havana_transcript=OTTHUMT00000007022;havana_version=1;transcript_id=ENST00000427857;transcript_support_level=3;version=1
"""
    stream = GFF3.Reader(IOBuffer(test_input))
    collect(stream)
    @test !GFF3.hasfasta(stream)
    @test_throws Exception GFF3.getfasta(stream)

    # implicit fasta
    test_input2 = string(test_input, """
>seq1
ACGTACGT
>seq2
TGCATGCA
""")
    stream = GFF3.Reader(IOBuffer(test_input2))
    collect(stream)
    @test GFF3.hasfasta(stream)
    @test collect(GFF3.getfasta(stream)) ==
        [FASTA.Record("seq1", dna"ACGTACGT")
         FASTA.Record("seq2", dna"TGCATGCA")]

    # explicit fasta
    test_input3 = string(test_input, """
##FASTA
>seq1
ACGTACGT
>seq2
TGCATGCA
""")
    stream = GFF3.Reader(IOBuffer(test_input3))
    collect(stream)
    @test GFF3.hasfasta(stream)
    @test collect(GFF3.getfasta(stream)) ==
        [FASTA.Record("seq1", dna"ACGTACGT")
         FASTA.Record("seq2", dna"TGCATGCA")]

    test_input4 = """
##directive1
#comment1
##directive2
1	havana	exon	869528	869575	.	-	.	Parent=transcript:ENST00000432963;Name=ENSE00001605362;constitutive=0;ensembl_end_phase=-1;ensembl_phase=-1;exon_id=ENSE00001605362;rank=2;version=1
1	havana	exon	870086	870201	.	-	.	Parent=transcript:ENST00000432963;Name=ENSE00001791782;constitutive=0;ensembl_end_phase=-1;ensembl_phase=-1;exon_id=ENSE00001791782;rank=1;version=1
##directive3
#comment2
##directive4
1	havana	lincRNA	868403	876802	.	-	.	ID=transcript:ENST00000427857;Parent=gene:ENSG00000230368;Name=FAM41C-002;biotype=lincRNA;havana_transcript=OTTHUMT00000007022;havana_version=1;transcript_id=ENST00000427857;transcript_support_level=3;version=1
##directive5
#comment3
##directive6
"""
    stream = GFF3.Reader(IOBuffer(test_input4), save_directives=true)
    read(stream)
    @test GFF3.directives(stream) == ["directive1", "directive2"]
    read(stream)
    @test isempty(GFF3.directives(stream))
    read(stream)
    @test GFF3.directives(stream) == ["directive3", "directive4"]
    @test_throws EOFError read(stream)
    @test eof(stream)
    close(stream)
    @test GFF3.directives(stream) == ["directive5", "directive6"]

    test_input5 = """
    ##directive1
    feature1\t.\t.\t.\t.\t.\t.\t.\t
    #comment1
    feature2\t.\t.\t.\t.\t.\t.\t.\t
    ##directive2
    feature3\t.\t.\t.\t.\t.\t.\t.\t
    """
    @test [r.kind for r in GFF3.Reader(IOBuffer(test_input5))] == [:feature, :feature, :feature]
    @test [r.kind for r in GFF3.Reader(IOBuffer(test_input5), skip_directives=false)] == [:directive, :feature, :feature, :directive, :feature]
    @test [r.kind for r in GFF3.Reader(IOBuffer(test_input5), skip_directives=false, skip_comments=false)] == [:directive, :feature, :comment, :feature, :directive, :feature]

    @testset "eachoverlap" begin
        path = joinpath(get_bio_fmt_specimens(), "GFF3", "TAIR10.part.gff.bgz")
        stream = BGZFStreams.BGZFStream(path)
        reader = GFF3.Reader(stream, index=string(path, ".tbi"))
        for (interval, n_records) in [
                (Interval("ChrC", 1:70_000), 68),
                (Interval("Chr4", 1:10_000),  2),
                (Interval("ChrM", 1:30_000), 10),]
            n = 0
            for record in eachoverlap(reader, interval)
                n += 1
            end
            @test n == n_records
        end
        @test isa(GFF3.Reader(path), GFF3.Reader)
    end
end

@testset "BigWig" begin
    @testset "empty" begin
        buffer = IOBuffer()
        data = buffer.data
        writer = BigWig.Writer(buffer, [("chr1", 1000)])
        close(writer)
        reader = BigWig.Reader(IOBuffer(data))
        @test length(collect(reader)) == 0
    end

    @testset "small" begin
        buffer = IOBuffer()
        data = buffer.data
        writer = BigWig.Writer(buffer, [("chr1", 1000)])
        write(writer, ("chr1", 50, 100, 3.14))
        close(writer)
        reader = BigWig.Reader(IOBuffer(data))
        records = collect(reader)
        @test length(records) == 1
        @test BigWig.haschrom(records[1]) === hasseqname(records[1]) === true
        @test BigWig.chrom(records[1]) == seqname(records[1]) == "chr1"
        @test BigWig.haschromstart(records[1]) === hasleftposition(records[1]) === true
        @test BigWig.chromstart(records[1]) === leftposition(records[1]) === 50
        @test BigWig.haschromend(records[1]) === hasrightposition(records[1]) === true
        @test BigWig.chromend(records[1]) === rightposition(records[1]) === 100
        @test BigWig.hasvalue(records[1])
        @test BigWig.value(records[1]) === 3.14f0
        @test startswith(repr(records[1]), "GenomicFeatures.BigWig.Record:\n")
        interval = convert(Interval, records[1])
        @test seqname(interval) == "chr1"
        @test leftposition(interval) === 50
        @test rightposition(interval) === 100
        @test metadata(interval) === 3.14f0
        @test all(isnan.(BigWig.values(reader, "chr1", 1:49)))
        @test BigWig.values(reader, "chr1", 50:51) == [3.14f0, 3.14f0]
        @test BigWig.values(reader, "chr1", 99:100) == [3.14f0, 3.14f0]
        @test all(isnan.(BigWig.values(reader, "chr1", 101:200)))
        @test BigWig.values(reader, Interval("chr1", 55, 56)) == [3.14f0, 3.14f0]

        # bedgraph (default)
        buffer = IOBuffer()
        data = buffer.data
        writer = BigWig.Writer(buffer, [("chr1", 1000)]; datatype=:bedgraph)
        write(writer, ("chr1",  1, 10, 0.0))
        write(writer, ("chr1", 15, 15, 1.0))
        write(writer, ("chr1", 90, 99, 2.0))
        close(writer)
        reader = BigWig.Reader(IOBuffer(data))
        records = collect(reader)
        @test length(records) == 3
        @test BigWig.chrom.(records) == ["chr1", "chr1", "chr1"]
        @test BigWig.chromstart.(records) == [1, 15, 90]
        @test BigWig.chromend.(records) == [10, 15, 99]
        @test BigWig.value.(records) == [0.0, 1.0, 2.0]

        # varstep
        buffer = IOBuffer()
        data = buffer.data
        writer = BigWig.Writer(buffer, [("chr1", 1000)]; datatype=:varstep)
        write(writer, ("chr1",  1, 10, 0.0))
        write(writer, ("chr1", 15, 24, 1.0))
        write(writer, ("chr1", 90, 99, 2.0))
        close(writer)
        reader = BigWig.Reader(IOBuffer(data))
        records = collect(reader)
        @test length(records) == 3
        @test BigWig.chrom.(records) == ["chr1", "chr1", "chr1"]
        @test BigWig.chromstart.(records) == [1, 15, 90]
        @test BigWig.chromend.(records) == [10, 24, 99]
        @test BigWig.value.(records) == [0.0, 1.0, 2.0]

        # fixedstep
        buffer = IOBuffer()
        data = buffer.data
        writer = BigWig.Writer(buffer, [("chr1", 1000)]; datatype=:fixedstep)
        write(writer, ("chr1",  1,  5, 0.0))
        write(writer, ("chr1", 11, 15, 1.0))
        write(writer, ("chr1", 21, 25, 2.0))
        close(writer)
        reader = BigWig.Reader(IOBuffer(data))
        records = collect(reader)
        @test length(records) == 3
        @test BigWig.chrom.(records) == ["chr1", "chr1", "chr1"]
        @test BigWig.chromstart.(records) == [1, 11, 21]
        @test BigWig.chromend.(records) == [5, 15, 25]
        @test BigWig.value.(records) == [0.0, 1.0, 2.0]
    end

    @testset "large" begin
        buffer = IOBuffer()
        data = buffer.data
        binsize = 32
        writer = BigWig.Writer(buffer, [("chr1", 100_000), ("chr2", 90_000)], binsize=binsize)
        for i in 1:10_000
            write(writer, ("chr1", (i-1)*10+1, i*10, log(i)))
        end
        n = 0
        p = 1
        while p ≤ 90_000
            sz = min(rand(1:100), 90_000 - p)
            write(writer, ("chr2", p, p+sz, log(p)))
            n += 1
            p += sz + 1
        end
        close(writer)
        reader = BigWig.Reader(IOBuffer(data))
        records = collect(reader)
        @test length(records) == 10_000 + n
        records = collect(eachoverlap(reader, Interval("chr1", 50_001, 50_165)))
        @test length(records) == 17
        @testset for bin in [1, 5, 10, 51, 300]
            for scale in 1:2
                binsize_scaled = binsize * BigWig.ZOOM_SCALE_FACTOR^(scale-1)
                chromstart = (bin - 1) * binsize_scaled + 1
                chromend = bin * binsize_scaled
                @test BigWig.coverage(reader, "chr1", chromstart, chromend; usezoom=false) == BigWig.coverage(reader, "chr1", chromstart, chromend; usezoom=true)
                @test BigWig.mean(reader, "chr1", chromstart, chromend; usezoom=false)     ≈  BigWig.mean(reader, "chr1", chromstart, chromend; usezoom=true)
                @test BigWig.minimum(reader, "chr1", chromstart, chromend; usezoom=false)  ≈  BigWig.minimum(reader, "chr1", chromstart, chromend; usezoom=true)
                @test BigWig.maximum(reader, "chr1", chromstart, chromend; usezoom=false)  ≈  BigWig.maximum(reader, "chr1", chromstart, chromend; usezoom=true)
                # TODO: use more stable algorithm?
                #@test_approx_eq BigWig.std(reader, "chr1", chromstart, chromend; usezoom=false)     BigWig.std(reader, "chr1", chromstart, chromend; usezoom=true)
            end
        end
    end

    @testset "round trip" begin
        function test_round_trip(filepath)
            reader = open(BigWig.Reader, filepath)
            buffer = IOBuffer()
            data = buffer.data
            writer = BigWig.Writer(buffer, BigWig.chromlist(reader))
            original = []
            for record in reader
                t = (BigWig.chrom(record), BigWig.chromstart(record), BigWig.chromend(record), BigWig.value(record))
                write(writer, t)
                push!(original, t)
            end
            close(writer)
            close(reader)

            reader = BigWig.Reader(IOBuffer(data))
            copy = []
            for record in reader
                t = (BigWig.chrom(record), BigWig.chromstart(record), BigWig.chromend(record), BigWig.value(record))
                push!(copy, t)
            end
            close(reader)

            @test original == copy
        end

        dir = joinpath(get_bio_fmt_specimens(), "BBI")
        for specimen in YAML.load_file(joinpath(dir, "index.yml"))
            valid = get(specimen, "valid", true)
            bigwig = "bigwig" ∈ split(specimen["tags"])
            if valid && bigwig
                test_round_trip(joinpath(dir, specimen["filename"]))
            end
        end
    end
end

@testset "BigBed" begin
    @testset "empty" begin
        buffer = IOBuffer()
        data = buffer.data
        writer = BigBed.Writer(buffer, [("chr1", 1000)])
        close(writer)
        reader = BigBed.Reader(IOBuffer(data))
        @test length(collect(reader)) == 0
    end

    @testset "small" begin
        buffer = IOBuffer()
        data = buffer.data
        writer = BigBed.Writer(buffer, [("chr1", 1000)])
        write(writer, ("chr1", 50, 100, "name1"))
        close(writer)
        reader = BigBed.Reader(IOBuffer(data))
        records = collect(reader)
        @test length(records) == 1
        @test BigBed.chrom(records[1]) == "chr1"
        @test BigBed.chromstart(records[1]) === 50
        @test BigBed.chromend(records[1]) === 100
        @test BigBed.name(records[1]) == "name1"
        @test !BigBed.hasscore(records[1])
        @test BigBed.optionals(records[1]) == ["name1"]

        buffer = IOBuffer()
        data = buffer.data
        writer = BigBed.Writer(buffer, [("chr1", 1000)])
        write(writer, ("chr1", 1, 100, "some name", 100, '+', 10, 90, RGB(0.5, 0.1, 0.2), 2, [4, 10], [10, 20]))
        close(writer)
        reader = BigBed.Reader(IOBuffer(data))
        records = collect(reader)
        @test length(records) == 1
        @test BigBed.haschrom(records[1]) === hasseqname(records[1]) === true
        @test BigBed.chrom(records[1]) == seqname(records[1]) == "chr1"
        @test BigBed.haschromstart(records[1]) === hasleftposition(records[1]) === true
        @test BigBed.chromstart(records[1]) === leftposition(records[1]) === 1
        @test BigBed.haschromend(records[1]) === hasrightposition(records[1]) === true
        @test BigBed.chromend(records[1]) === rightposition(records[1]) === 100
        @test BigBed.hasname(records[1])
        @test BigBed.name(records[1]) == "some name"
        @test BigBed.hasscore(records[1])
        @test BigBed.score(records[1]) === 100
        @test BigBed.hasstrand(records[1])
        @test BigBed.strand(records[1]) === STRAND_POS
        @test BigBed.hasthickstart(records[1])
        @test BigBed.thickstart(records[1]) === 10
        @test BigBed.hasthickend(records[1])
        @test BigBed.thickend(records[1]) === 90
        @test BigBed.hasitemrgb(records[1])
        @test BigBed.itemrgb(records[1]) === convert(RGB{N0f8}, RGB(0.5, 0.1, 0.2))
        @test BigBed.hasblockcount(records[1])
        @test BigBed.blockcount(records[1]) === 2
        @test BigBed.hasblocksizes(records[1])
        @test BigBed.blocksizes(records[1]) == [4, 10]
        @test BigBed.hasblockstarts(records[1])
        @test BigBed.blockstarts(records[1]) == [10, 20]
        @test BigBed.optionals(records[1]) == ["some name", "100", "+", "9", "90", "128,26,51", "2", "4,10,", "9,19,"]
    end

    @testset "large" begin
        buffer = IOBuffer()
        data = buffer.data
        binsize = 32
        writer = BigBed.Writer(buffer, [("chr1", 100_000), ("chr2", 90_000)], binsize=binsize)
        for i in 1:10_000
            write(writer, ("chr1", (i-1)*10+1, i*10, string("name", i)))
        end
        n = 0
        p = 1
        while p ≤ 90_000
            sz = min(rand(1:100), 90_000 - p)
            write(writer, ("chr2", p, p+sz, string("name", n + 1)))
            n += 1
            p += sz + 1
        end
        close(writer)
        reader = BigBed.Reader(IOBuffer(data))
        records = collect(reader)
        @test length(records) == 10_000 + n
        records = collect(eachoverlap(reader, Interval("chr1", 50_001, 50_165)))
        @test length(records) == 17
    end

    @testset "round trip" begin
        function test_round_trip(filepath)
            reader = open(BigBed.Reader, filepath)
            buffer = IOBuffer()
            data = buffer.data
            writer = BigBed.Writer(buffer, BigBed.chromlist(reader))
            original = []
            for record in reader
                t = (BigBed.chrom(record), BigBed.chromstart(record), BigBed.chromend(record), BigBed.optionals(record)...)
                write(writer, t)
                push!(original, t)
            end
            close(writer)
            close(reader)

            reader = BigBed.Reader(IOBuffer(data))
            copy = []
            for record in reader
                t = (BigBed.chrom(record), BigBed.chromstart(record), BigBed.chromend(record), BigBed.optionals(record)...)
                push!(copy, t)
            end
            close(reader)

            @test original == copy
        end

        dir = joinpath(get_bio_fmt_specimens(), "BBI")
        for specimen in YAML.load_file(joinpath(dir, "index.yml"))
            valid = get(specimen, "valid", true)
            bigbed = "bigbed" ∈ split(specimen["tags"])
            if valid && bigbed
                test_round_trip(joinpath(dir, specimen["filename"]))
            end
        end
    end

    @testset "overlap" begin
        chromlen = 1_000_000
        Random.seed!(1234)
        chroms = ["one", "two", "three", "four", "five"]
        intervals = IntervalCollection(
             [Interval(i.seqname, i.first, i.last)
             for i in random_intervals(chroms, chromlen, 10_000)], true)

        buffer = IOBuffer()
        data = buffer.data
        writer = BigBed.Writer(buffer, [(chrom, chromlen) for chrom in chroms])
        for i in intervals
            write(writer, i)
        end
        close(writer)

        reader = BigBed.Reader(IOBuffer(data))
        queries = random_intervals(chroms, chromlen, 1000)
        triplet(x::Interval) = String(x.seqname), x.first, x.last
        triplet(x::BigBed.Record) = BigBed.chrom(x), BigBed.chromstart(x), BigBed.chromend(x)
        @test all(triplet.(collect(eachoverlap(intervals, q))) == triplet.(collect(eachoverlap(reader, q))) for q in queries)
        close(reader)
    end
end


text:
module KdbMode
    import Base: LineEdit, REPL
    import Q._k: k, kj, kp
    import Q: K_new, chkparens
    using Q

    const PROMPT = "q)"
    const PROMPT_COLOR = :light_blue
    const MODE_NAME = :q
    const MODE_KEY = '\\'

    function __init__()
        isdefined(Base, :active_repl) && install_kdb_mode(Base.active_repl)
    end

    function install_kdb_mode(repl)
        handle = Q.GOT_Q ? 0 : Q.open_default_kdb_handle(repl.t.err_stream)
        if handle >= 0
            install_kdb_mode(repl, handle)
            info(repl.t.err_stream, "Press \u005c\u20e3 for q) prompt.")
            true
        else
            false
        end
    end

    function on_enter_do(prompt_state)
        !isempty(prompt_state) && begin
            chkparens(LineEdit.buffer(prompt_state).data) == 0
        end
    end

    function on_done_do_function(repl, handle)
        function(cmd)
            out = repl.t.out_stream
            x = K[]
            n = count(c == '\n' for c in cmd)
            size = max.([10, 0], Int[displaysize(out)...] - [n, 0])
            try
                x = K(k(handle, "{.Q.S[x-3 0;y;]value z}",
                        K_new(size), kj(0), kp("get$(repr(cmd))")))
            catch error
                if error isa KdbException
                    print_with_color(:red, out, "'", error.s)
                    println(out)
                else
                    println(out, error)
                end
            end
            for line in x
                println(out, line)
            end
        end
    end

    function create_kdb_panel(repl, handle)
        on_done_do = on_done_do_function(repl, handle)
        main_mode = repl.interface.modes[1]
        # Setup q) panel
        panel = LineEdit.Prompt(PROMPT;
            prompt_prefix=Base.text_colors[PROMPT_COLOR],
            prompt_suffix=Base.text_colors[:white],
            on_enter=on_enter_do)

        panel.on_done = REPL.respond(on_done_do, repl, panel)

        push!(repl.interface.modes, panel)

        hp = main_mode.hist
        hp.mode_mapping[MODE_NAME] = panel
        panel.hist = hp

        search_prompt, skeymap = LineEdit.setup_search_keymap(hp)
        mk = REPL.mode_keymap(main_mode)

        b = Dict{Any,Any}[skeymap, mk,
                LineEdit.history_keymap,
                LineEdit.default_keymap,
                LineEdit.escape_defaults]
        panel.keymap_dict = LineEdit.keymap(b)

        panel
    end

    function install_kdb_mode(repl, handle)
        main_mode = repl.interface.modes[1]
        panel = create_kdb_panel(repl, handle)

        # Install this mode into the main mode
        q_keymap = Dict{Any,Any}(
            MODE_KEY => function (s, args...)
                if isempty(s) || position(LineEdit.buffer(s)) == 0
                    buf = copy(LineEdit.buffer(s))
                    LineEdit.transition(s, panel) do
                        LineEdit.state(s, panel).input_buffer = buf
                    end
                else
                    LineEdit.edit_insert(s, MODE_KEY)
                end
            end
        )
        LineEdit.keymap_merge!(main_mode.keymap_dict, q_keymap);
        nothing
    end
end


text:
mutable struct Tape
  tape::Array{Int, 1}
  pos::Int

  Tape() = new(zeros(Int, 1), 1)
end

inc(this::Tape) = (this.tape[this.pos] += 1)
dec(this::Tape) = (this.tape[this.pos] -= 1)
get(this::Tape) = this.tape[this.pos]

function advance(this::Tape)
  this.pos += 1
  if this.pos > length(this.tape)
    push!(this.tape, 0)
  end
end

function devance(this::Tape)
  if this.pos > 1
    this.pos -= 1
  end
end

validbfsymbol(x) = in(x, ['>', '<', '+', '-', '.', ',', '[', ']'])

struct Program
  code::String
  bracket_map::Dict{Int, Int}

  function Program(text)
    code = filter(validbfsymbol, text)
    bracket_map = Dict{Int, Int}()
    stack = Int[]
    pc = 1
    for ch in code
      if ch == '['
        push!(stack, pc)
      elseif ch == ']' && length(stack) > 0
        right = pop!(stack)
        bracket_map[pc] = right
        bracket_map[right] = pc
      end
      pc += 1
    end
    return new(code, bracket_map)
  end
end

function run(this::Program)
  code = this.code
  bracket_map = this.bracket_map
  tape = Tape()
  pc = 1
  while pc <= length(code)
    ch = code[pc]
    if ch == '+'
      inc(tape)
    elseif ch == '-'
      dec(tape)
    elseif ch == '>'
      advance(tape)
    elseif ch == '<'
      devance(tape)
    elseif ch == '['
      if get(tape) == 0
        pc = bracket_map[pc]
      end
    elseif ch == ']'
      if get(tape) != 0
        pc = bracket_map[pc]
      end
    elseif ch == '.'
      print(Char(get(tape)))
    end
    pc += 1
  end
end

function main(text)
  run(Program(text))
end

println("JIT warming up")
main(">++[<+++++++++++++>-]<[[>+>+<<-]>[<+>-]++++++++[>++++++++<-]>[-]<<>++++++++++[>++++++++++[>++++++++++[>++++++++++[>++++++++++[>++++++++++[>++++++++++[-]<-]<-]<-]<-]<-]<-]<-]++++++++++")

println("bench")
text = open(ARGS[1]) do file
  read(file, String)
end
x = @timed main(text)
println("Elapsed: $(x[2]), Allocated: $(x[3]), GC Time: $(x[4])")

max_stars_repo_name: wookay/benchmarks

max_issues_repo_name: wookay/benchmarks

max_forks_repo_name: wookay/benchmarks


text:
import Observables
import AbstractPlotting

q1 = Biquaternion(Quaternion(rand(4)), ℝ³(rand(3)))
scene = AbstractPlotting.Scene()
radius = rand()
segments = rand(5:10)
color = AbstractPlotting.RGBAf0(rand(4)...)
transparency = false
sphere = Sphere(q1,
                scene,
                radius = radius,
                segments = segments,
                color = color,
                transparency = transparency)
q2 = q1 * Biquaternion(Quaternion(rand(4)), ℝ³(rand(3)))
value1 = getsurface(sphere.observable, segments)
update(sphere, q2)
value2 = getsurface(sphere.observable, segments)

@test isapprox(sphere.q, q2)
@test isapprox(value1, value2) == false


text:
module DeZero

export Variable, backward, cleargrad, add, square

include("step16.jl")

end # module


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .jl
#       format_name: light
#       format_version: '1.3'
#       jupytext_version: 0.8.6
#   kernelspec:
#     display_name: Julia 1.0.3
#     language: julia
#     name: julia-1.0
# ---

using Random
using DifferentialEquations

# BM model for beta(t): d beta(t) = beta(t) d W(t)
# Function to compute rates of change
spsir_bm = function(du, u, p, t)
    #Susceptible Individuals
    S = u[1]
    #Infected Individuals
    I = u[2]
    #Recovered Individuals
    R = u[3]

    β = max(0., u[4])

    N = S + I + R
    γ = p[1]

    du[1] = -β * S * I / N
    du[2] = β * S * I / N - γ * I
    du[3] = γ * I
    du[4] = 0.
end

# Function to add noise
sigma_spsir = function(du, u, p, t )
    σ = p[2]

    du[1] = 0.
    du[2] = 0.
    du[3] = 0.
    du[4] = σ
end

# +
# BM for logbeta(t), with drift
spsir_logbm_drift = function(du, u, p, t)
    S = u[1]
    I = u[2]
    R = u[3]
    β = exp( u[4] )

    N = S + I + R
    γ = p[1]
    α = p[3]

    du[1] = -β * S * I / N
    du[2] = β * S * I / N - γ * I
    du[3] = γ * I
    du[4] = -α * I
end

# +
# set random seed
Random.seed!( 1111 )

## Simulation of BM model
# starting conditions
u0 = [50.;1.0;0.0;2.0]
tspan = (0.0,10.0)
# parameters gamma sigma
p = [1.; 1.]
spsir_bm_prob = SDEProblem(spsir_bm, sigma_spsir, u0, tspan, p)
spsir_bm_sol = solve(spsir_bm_prob)

# +
## Simulation of log-BM with drift model

# starting condtions
u0 = [50.;1.0;0.0;log(3.) ]
tspan = (0.0,10.0)
# parameters gamma sigma alpha
p = [1.; 1.; 0.1]
spsir_logbm_drift_prob = SDEProblem(spsir_logbm_drift, sigma_spsir, u0, tspan, p)
spsir_logbm_drift_sol = solve(spsir_logbm_drift_prob)

# -

## Plotting for BM model
using Plots
# Plot evolution of number infected
plot( spsir_bm_sol , vars = 2 )

# Plot evolution of transmission rate
plot( spsir_bm_sol , vars = 4 )

# +
## Plotting for BM with drift model
# Plot evolution of number infected
plot( spsir_logbm_drift_sol , vars = 2 )

# -

# Plot evolution of transmission rate
plot( spsir_logbm_drift_sol , vars = 4 )

max_stars_repo_name: epirecipes/epiRecipes

max_issues_repo_name: epirecipes/epiRecipes

max_forks_repo_name: epirecipes/epiRecipes


text:
import Sockets.connect, Sockets.TCPSocket, Base.StatusActive, Base.StatusOpen, Base.StatusPaused

abstract type RedisConnectionBase end
abstract type SubscribableConnection<:RedisConnectionBase end

struct RedisConnection <: SubscribableConnection
    host::AbstractString
    port::Integer
    password::AbstractString
    db::Integer
    socket::TCPSocket
end

struct SentinelConnection <: SubscribableConnection
    host::AbstractString
    port::Integer
    password::AbstractString
    db::Integer
    socket::TCPSocket
end

struct TransactionConnection <: RedisConnectionBase
    host::AbstractString
    port::Integer
    password::AbstractString
    db::Integer
    socket::TCPSocket
end

mutable struct PipelineConnection <: RedisConnectionBase
    host::AbstractString
    port::Integer
    password::AbstractString
    db::Integer
    socket::TCPSocket
    num_commands::Integer
end

struct SubscriptionConnection <: RedisConnectionBase
    host::AbstractString
    port::Integer
    password::AbstractString
    db::Integer
    callbacks::Dict{AbstractString, Function}
    pcallbacks::Dict{AbstractString, Function}
    socket::TCPSocket
end

function RedisConnection(; host="127.0.0.1", port=6379, password="", db=0)
    try
        socket = connect(host, port)
        connection = RedisConnection(host, port, password, db, socket)
        on_connect(connection)
    catch
        throw(ConnectionException("Failed to connect to Redis server"))
    end
end

function SentinelConnection(; host="127.0.0.1", port=26379, password="", db=0)
    try
        socket = connect(host, port)
        sentinel_connection = SentinelConnection(host, port, password, db, socket)
        on_connect(sentinel_connection)
    catch
        throw(ConnectionException("Failed to connect to Redis sentinel"))
    end
end

function TransactionConnection(parent::RedisConnection)
    try
        socket = connect(parent.host, parent.port)
        transaction_connection = TransactionConnection(parent.host,
            parent.port, parent.password, parent.db, socket)
        on_connect(transaction_connection)
    catch
        throw(ConnectionException("Failed to create transaction"))
    end
end

function PipelineConnection(parent::RedisConnection)
    try
        socket = connect(parent.host, parent.port)
        pipeline_connection = PipelineConnection(parent.host,
            parent.port, parent.password, parent.db, socket, 0)
        on_connect(pipeline_connection)
    catch
        throw(ConnectionException("Failed to create pipeline"))
    end
end

function SubscriptionConnection(parent::SubscribableConnection)
    try
        socket = connect(parent.host, parent.port)
        subscription_connection = SubscriptionConnection(parent.host,
            parent.port, parent.password, parent.db, Dict{AbstractString, Function}(),
            Dict{AbstractString, Function}(), socket)
        on_connect(subscription_connection)
    catch
        throw(ConnectionException("Failed to create subscription"))
    end
end

function on_connect(conn::RedisConnectionBase)
    conn.password != "" && auth(conn, conn.password)
    conn.db != 0        && select(conn, conn.db)
    conn
end

function disconnect(conn::RedisConnectionBase)
    close(conn.socket)
end

function is_connected(conn::RedisConnectionBase)
    conn.socket.status == StatusActive || conn.socket.status == StatusOpen || conn.socket.status == StatusPaused
end

function send_command(conn::RedisConnectionBase, command::AbstractString)
    write(conn.socket, command)
end


text:
"""
	mutable struct DictEntry <: JSONEntry
		childs::Dict{String, Any}
		updated::Int
	end

keeps statistics about an object in json
- `childs` maintains key-value statistics of childrens. All values should be JSONEntries
- `updated` counts how many times the struct was updated.
"""

mutable struct DictEntry <: JSONEntry
	childs::Dict{Symbol, Any}
	updated::Int
end

DictEntry() = DictEntry(Dict{Symbol,Any}(),0)
Base.getindex(s::DictEntry, k::Symbol) = s.childs[k]
Base.setindex!(s::DictEntry, i, k::Symbol) = s.childs[k] = i
Base.get(s::Dict{Symbol, <:Any}, key::String, default) = get(s, Symbol(key), default)
Base.keys(e::DictEntry) = keys(e.childs)
Base.isempty(e::DictEntry) = isempty(e.childs)

function update!(s::DictEntry, d::Dict; path = "")
	s.updated += 1
	for (k,v) in d
		kc = Symbol(k)
		isnothing(v) && continue
		# isempty(v) && continue
		if haskey(s.childs, kc)
			s.childs[kc] = safe_update!(s.childs[kc], v, path="$path[$k]")
		else
			o = newentry!(v)
			if !isnothing(o)
				s.childs[kc] = o
			end
		end
	end
	return true
end

"""
	suggestextractor(e::DictEntry, settings = NamedTuple())

create convertor of json to tree-structure of `DataNode`

- `e` top-level of json hierarchy, typically returned by invoking schema
- `settings` can be any container supporting `get` function
- `settings.mincountkey` contains minimum repetition of the key to be included into the extractor
  (if missing it is equal to zero)
- `settings.key_as_field` of the number of keys exceeds this value, it is assumed that keys contains a value,
  which means that they will be treated as strings.
- `settings.scalar_extractors` contains rules for determining which extractor to use for leaves.
  Default value is return value of `default_scalar_extractor()`, it's array of pairs where first element is predicate
  and if it matches, second element, function which maps schema to specific extractor, is called.
"""
function suggestextractor(e::DictEntry, settings = NamedTuple(); path = "", child_less_than_parent = false)
	length(e.childs) >= get(settings, :key_as_field, 500) && return key_as_field(e, settings;
		path = path, child_less_than_parent = child_less_than_parent)

	for k in filter(k->!isnothing(e.childs[k]) && isempty(e.childs[k]), keys(e.childs))
		@warn "$(path): key $k contains empty array, skipping"
	end
	ks = filter(k->!isempty(e.childs[k]), keys(e.childs))
	mincount = get(settings, :mincountkey, 0)
	ks = filter(k -> updated(e.childs[k]) > mincount, ks)
	isempty(ks) && return nothing
	c = [(k,suggestextractor(e.childs[k], settings,
			path = path*"[:$(k)]",
			child_less_than_parent = child_less_than_parent || e.updated > e.childs[k].updated)
		) for k in ks]
	c = filter(s -> s[2] != nothing, c)
	isempty(c) && return nothing
	ExtractDict(Dict(c))
end

function key_as_field(e::DictEntry, settings; path = "", child_less_than_parent = false)
	@info "$(path) seems to store values in keys, therefore node is treated as bag with keys as extra values."
	child_schema = reduce(merge, collect(values(e.childs)), init = nothing)
	key_schema = Entry(String(first(keys(e))))
	for k in keys(e)
		update!(key_schema, k, path=path)
	end
	ExtractKeyAsField(ExtractString(3, 256, 2053, true), suggestextractor(child_schema, settings, path = path*"[:childs]"))
end

"""
Dispatch of Base.merge on JsonGrinder.JSONEntry structures.
Allows to merge multiple schemas to single one.

	merge(es::Entry...)
	merge(es::DictEntry...)
	merge(es::ArrayEntry...)
	merge(es::MultiEntry...)
	merge(es::JsonGrinder.JSONEntry...)

it can be used to distribute calculation of schema across multiple workers to merge their partial results into bigger one.

# Example

If we want to calculate schema from e.g. array of jsons in a distributed manner, if we have `jsons` array and ,
we can do it using
```julia
using ThreadsX
ThreadsX.mapreduce(schema, merge, Iterators.partition(jsons, length(jsons) ÷ Threads.nthreads()))
```
or
```julia
using ThreadTools
merge(tmap(schema, Threads.nthreads(), Iterators.partition(jsons, length(jsons) ÷ Threads.nthreads()))
```
or, if you like to split it into multiple jobs and having them processed by multiple threads, it can look like
```julia
using ThreadTools
merge(tmap(schema, Threads.nthreads(), Iterators.partition(jsons, 1_000))
```
where we split array to smaller array of size 1k and let all available threads create partial schemas.

If your data is too large to fit into ram, following approach works well also with filenames and similar other ways to process large data.
"""
function merge(es::DictEntry...)
	updates_merged = sum(updated.(es))
	childs_merged = merge(merge, childs.(es)...)
	DictEntry(childs_merged, updates_merged)
end

childs(s::T) where {T<:DictEntry} = s.childs
Base.hash(e::DictEntry, h::UInt) = hash((e.childs, e.updated), h)
Base.:(==)(e1::DictEntry, e2::DictEntry) = e1.updated === e2.updated && e1.childs == e2.childs

sample_synthetic(e::DictEntry) = Dict(k => sample_synthetic(v) for (k, v) in e.childs)


text:
using WindMixing: plot_loss

FILE_NAME = "NDE_3sim_diurnal_18simBFGST0.8nograd_divide1f5_gradient_smallNN_leakyrelu_rate_2e-4_T0.8"
DATA_PATH = joinpath(pwd(), "extracted_training_output", "$(FILE_NAME)_extracted.jld2")
OUTPUT_PATH = "C:\\Users\\xinle\\Documents\\OceanParameterizations.jl"

file = jldopen(DATA_PATH, "r")
losses = (
    u = file["losses/u"],
    v = file["losses/v"],
    T = file["losses/T"],
    ∂u∂z = file["losses/∂u∂z"],
    ∂v∂z = file["losses/∂v∂z"],
    ∂T∂z = file["losses/∂T∂z"],
)

train_files = file["training_info/train_files"]

diurnal = occursin("diurnal", train_files[1])

train_parameters = file["training_info/parameters"]
uw_NN = file["neural_network/uw"]
vw_NN = file["neural_network/vw"]
wT_NN = file["neural_network/wT"]

loss_scalings = (u=1f0, v=1f0, T=1f0, ∂u∂z=1f0, ∂v∂z=1f0, ∂T∂z=1f0)

close(file)

plot_loss(losses, joinpath(OUTPUT_PATH, "$(FILE_NAME)_loss.png"))


text:
abstract type RDataException <: Exception end

struct UnsupportedROBJ <: RDataException
    sxtype::SXType
    msg::String

    UnsupportedROBJ(sxtype::SXType, msg::AbstractString) = new(sxtype, msg)
    UnsupportedROBJ(sxtype::SXType) =
        UnsupportedROBJ(sxtype, "Unsupported R object (sxtype=$(sxtype))")
end

Base.showerror(io::IO, e::UnsupportedROBJ) = print(io, e.msg)

struct CodecMissingError <: RDataException
    formatName::Symbol
end

function Base.showerror(io::IO, e::CodecMissingError)
    codecpkg = "Codec" * titlecase(String(e.formatName))
    print(io, string(
        "$(typeof(e)): $codecpkg package is required to read ",
        "$(e.formatName)-compressed RData files. Run ",
        "Pkg.add(\"$codecpkg\") to install it. Then in your code call ",
        "\"using $codecpkg\" before \"using RData\"."
    ))
end


text:
import COBRA
import Utils # Pkg.clone("https://github.com/cossio/Utils.jl")

export readCobraModel, prepareXmlModel

"""
Read model in MATLAB format.
"""
function readCobraModel(path::String, matrix::String = "S", model::String = "model")
    COBRA.loadModel(path, matrix, model)
end

"""
Read model in MATLAB format.
Returned stoichiometric matrix is sparse.
"""
function readCobraModelSparse(path::String, matrix::String = "S", model::String = "model")
    model = readCobraModel(path, matrix, model)
    COBRA.LPproblem(sparse(model.S), model.b, model.c, model.lb, model.ub,
                           model.osense, model.csense, model.rxns, model.mets)
end

"""
Converts SBML (.xml) model to MATLAB (.mat) format.
"""
function prepareXmlModel(xmlpath::String)
    Utils.gunzip(xmlpath)
    matpath = Utils.removeExtension(xmlpath) * ".mat"
    @assert isfile(xmlpath) || isfile(matpath)
    if !isfile(matpath)
        pkgdir = Pkg.Dir.path() * "/MaxEntChemostat2018/"
        run(`python $(pkgdir * "/scripts/save_as_mat.py") $xmlpath`)
    end
    return matpath
end

function loadEColiTestModel()
    path = Pkg.dir() * "/COBRAUtils/metabolic_networks/ecoli_core_model.mat"
    model = COBRAUtils.readCobraModel(path)
    model.osense = 1
    return model
end


text:
#=--- TerminalGraphics / integrations / Luxor.jl --------------------------=#

#=-------------------------------------------------------------------------=#
#
#   Integrations for 'Luxor.jl'
#
#   When Luxor is loaded, automatically allow showing a drawing on the
#   console by drawing the underlying file.
#
#=-------------------------------------------------------------------------=#

import Base: display
display(d::TerminalGraphics.SixelDisplay, drw::Luxor.Drawing) = begin
    fn = drw.filename
    isfile(fn) &&
        error("The Luxor drawing holds a filename '$fn' that doesn't exist.")
    drw.surfacetype == "png" ||
        error("Currently, only 'png' Luxor Drawings can be displayed. Sorry.")
    # draw the 'png'
    title = string(typeof(drw))
    size = "size ($(drw.width), $(drw.height))"
    println(d.io, "Drawing $title of $size, at $fn")
    TerminalGraphics.draw(d.io, fn)
end

#=--- end of file ---------------------------------------------------------=#