text:
module Gen

using ProtoBuf.GoogleProtoBuf
using ProtoBuf.GoogleProtoBufCompiler

using ProtoBuf
using Compat

import ProtoBuf: meta, @logmsg, DEF_REQ, DEF_FNUM, DEF_VAL, DEF_PACK

export gen

const JTYPES              = [Float64, Float32, Int64, UInt64, Int32, UInt64,  UInt32,  Bool, AbstractString, Any, Any, Array{UInt8,1}, UInt32, Int32, Int32, Int64, Int32, Int64]
const JTYPE_DEFAULTS      = [0,       0,       0,     0,      0,     0,       0,       false, "",    nothing, nothing, UInt8[], 0,     0,     0,       0,       0,     0;]

isprimitive(fldtype) = (1 <= fldtype <= 8) || (13 <= fldtype <= 18)

# maps protofile name to package name
const _packages = Dict{AbstractString,AbstractString}()

# maps protofile name to array of imported protofiles
const protofile_imports = Dict()

const _keywords = [
    "if", "else", "elseif", "while", "for", "begin", "end", "quote",
    "try", "catch", "return", "local", "abstract", "function", "macro",
    "ccall", "finally", "typealias", "break", "continue", "type",
    "global", "module", "using", "import", "export", "const", "let",
    "bitstype", "do", "baremodule", "importall", "immutable",
    "Type", "Enum", "Any", "DataType", "Base"
]

_module_postfix = false
_map_as_array = false

type Scope
    name::AbstractString
    syms::Array{AbstractString,1}
    files::Array{AbstractString,1}
    is_module::Bool
    children::Array{Scope,1}
    parent::Scope

    Scope(name::AbstractString) = new(name, AbstractString[], AbstractString[], false, Scope[])
    function Scope(name::AbstractString, parent::Scope)
        s = new(name, AbstractString[], AbstractString[], false, Scope[], parent)
        push!(parent.children, s)
        s
    end
end

function fullname(s::Scope)
    if isempty(s.name) || isempty(s.parent.name)
        return s.name
    end
    string(fullname(s.parent), s.parent.is_module ? "." : "_", s.name)
end

const top_scope = Scope("")
top_scope.is_module = false

function get_module_scope(parent::Scope, newname::AbstractString)
    newname = _module_postfix ? newname * "_pb" : newname
    for s in parent.children
        if s.name == newname
            return s
        end
    end
    s = Scope(newname, parent)
    s.is_module = true
    s
end

function qualify(name::AbstractString, scope::Scope)
    if name in scope.syms
        return pfx(name, scope)
    elseif isdefined(scope, :parent)
        return qualify(name, scope.parent)
    else
        error("unresolved name $name at scope $(scope.name)")
    end
end

function readreq(srcio::IO)
    req = ProtoBuf.instantiate(CodeGeneratorRequest)
    readproto(srcio, req)
    req
end

pfx(name::AbstractString, scope::Scope) = isempty(scope.name) ? name : (fullname(scope) * (scope.is_module ? "." : "_") * name)
splitmodule(name::AbstractString) = rsplit(name, '.'; limit=2)
function findmodule(name::AbstractString)
    mlen = 0
    mpkg = ""
    for pkg in values(_packages)
        orig = _module_postfix ? join(map(m->m[1:end-3], split(pkg, ".")), ".") : pkg
        if (length(pkg) > mlen) && startswith(name, orig)
            mlen = length(orig)
            mpkg = pkg
        end
    end
    (mpkg, replace((0 == mlen) ? name : name[(mlen+2):end], '.', '_'))
end

type DeferredWrite
    iob::IOBuffer
    depends::Array{AbstractString,1}
end
const _deferred = Dict{AbstractString,DeferredWrite}()
# set of fully-qualified names we've resolved, to handle dependencies in other files
const _all_resolved = Set{AbstractString}()

function defer(name::AbstractString, iob::IOBuffer, depends::AbstractString)
    @logmsg("defer $name due to $depends")
    if isdeferred(name)
        depsnow = _deferred[name].depends
        !(depends in depsnow) && push!(depsnow, depends)
        return
    end
    _deferred[name] = DeferredWrite(iob, AbstractString[depends])
    nothing
end

isdeferred(name::AbstractString) = haskey(_deferred, name)
function isresolved(dtypename::AbstractString, referenced_name::AbstractString, full_referenced_name::AbstractString, exports::Array{AbstractString,1})
    (dtypename == referenced_name) && return true
    for jtype in JTYPES
        (referenced_name == string(jtype)) && return true
    end
    if '.' in referenced_name
        return !isdeferred(referenced_name)
    elseif referenced_name in exports
        return !(isdeferred(referenced_name) || isdeferred(full_referenced_name))
    end
    false
end

function resolve(iob::IOBuffer, name::AbstractString)
    fully_resolved = AbstractString[]
    for (typ,dw) in _deferred
        idx = findfirst(dw.depends, name)
        (idx == 0) && continue
        splice!(dw.depends, idx)
        isempty(dw.depends) && push!(fully_resolved, typ)
    end

    # write all fully resolved entities
    for typ in fully_resolved
        @logmsg("resolved $typ")
        print(iob, takebuf_string(_deferred[typ].iob))
        delete!(_deferred, typ)
        push!(_all_resolved, typ)
    end

    # mark them resolved as well
    for typ in fully_resolved
        resolve(iob, typ)
    end
end

chk_keyword(name) = (name in _keywords) ? string('_', name) : name

function generate(io::IO, errio::IO, enumtype::EnumDescriptorProto, scope::Scope, exports::Array{AbstractString,1})
    enumname = pfx(enumtype.name, scope)
    sm = splitmodule(enumname)
    (length(sm) > 1) && (enumname = sm[2])
    enumname = chk_keyword(enumname)
    push!(scope.syms, enumname)

    @logmsg("begin enum $(enumname)")
    println(io, "type __enum_$(enumname) <: ProtoEnum")
    values = Int32[]
    for value::EnumValueDescriptorProto in enumtype.value
        # If we find that the field name is a keyword prepend it with `_`
        fldname = chk_keyword(value.name)
        println(io, "    $(fldname)::Int32")
        push!(values, value.number)
    end
    println(io, "    __enum_$(enumname)() = new($(join(values,',')))")
    println(io, "end #type __enum_$(enumname)")
    println(io, "const $(enumname) = __enum_$(enumname)()")
    println(io, "")
    push!(exports, enumname)
    @logmsg("end enum $(enumname)")
    nothing
end

function short_type_name(full_type_name::AbstractString, depends::Array{AbstractString,1})
    comps = split(full_type_name, '.')
    type_name = pop!(comps)
    while !isempty(comps) && !(join(comps, '.') in depends)
        type_name = (pop!(comps) * "." * type_name)
    end
    @logmsg("check $full_type_name against $depends || found $type_name")
    return type_name
end

function generate(outio::IO, errio::IO, dtype::DescriptorProto, scope::Scope, syntax::AbstractString, exports::Vector{AbstractString}, depends::Vector{AbstractString}, mapentries::Dict)
    io = IOBuffer()
    full_dtypename = dtypename = pfx(dtype.name, scope)
    sm = splitmodule(dtypename)
    modul,dtypename = (length(sm) > 1) ? (sm[1],sm[2]) : ("",dtypename)
    dtypename = chk_keyword(dtypename)
    full_dtypename = (modul=="") ? dtypename : "$(modul).$(dtypename)"
    @logmsg("begin type $(full_dtypename)")

    scope = Scope(dtype.name, scope)

    # check oneof
    oneof_names = Compat.String[]
    if isfilled(dtype, :oneof_decl)
        for oneof_decl in dtype.oneof_decl
            if isfilled(oneof_decl, :name)
                push!(oneof_names,  "@compat(Symbol(\"$(oneof_decl.name)\"))")
            end
        end
    end

    # generate enums
    if isfilled(dtype, :enum_type)
        for enum_type in dtype.enum_type
            generate(io, errio, enum_type, scope, exports)
            (errio.size > 0) && return
        end
    end

    # generate nested types
    if isfilled(dtype, :nested_type)
        for nested_type::DescriptorProto in dtype.nested_type
            generate(io, errio, nested_type, scope, syntax, exports, depends, mapentries)
            (errio.size > 0) && return
        end
    end

    # generate this type
    println(io, "type $(dtypename)")
    reqflds = AbstractString[]
    packedflds = AbstractString[]
    fldnums = Int[]
    defvals = AbstractString[]
    wtypes = AbstractString[]
    oneofs = isempty(oneof_names) ? Int[] : zeros(Int, length(dtype.field))
    ismapentry = isfilled(dtype, :options) && isfilled(dtype.options, :map_entry) && dtype.options.map_entry
    map_key_type = ""
    map_val_type = ""

    if isfilled(dtype, :field)
        for fld_idx in 1:length(dtype.field)
            field = dtype.field[fld_idx]
            if isfilled(field, :oneof_index) && !isfilled_default(field, :oneof_index)
                oneof_idx = field.oneof_index + 1
                oneofs[fld_idx] = oneof_idx
            end

            # If we find that the field name is a keyword prepend it with _type
            fldname = chk_keyword(field.name)
            if field._type == FieldDescriptorProto_Type.TYPE_GROUP
                println(errio, "Groups are not supported")
                return
            end

            full_typ_name = ""
            if (field._type == FieldDescriptorProto_Type.TYPE_MESSAGE) || (field._type == FieldDescriptorProto_Type.TYPE_ENUM)
                typ_name = field.type_name
                if startswith(typ_name, '.')
                    (m,t) = findmodule(typ_name[2:end])
                    t = chk_keyword(t)
                    full_typ_name = m=="" ? t : "$(m).$(t)"
                    typ_name = (m == modul) ? t : full_typ_name
                else
                    full_typ_name = qualify(typ_name, scope)
                    typ_name = full_typ_name
                    m,t = splitmodule(typ_name)
                    t = chk_keyword(t)
                    (m == modul) && (typ_name = t)
                    full_typ_name = m=="" ? t : "$(m).$(t)"
                end
            elseif field._type == FieldDescriptorProto_Type.TYPE_SINT32
                push!(wtypes, ":$fldname => :sint32")
            elseif field._type == FieldDescriptorProto_Type.TYPE_SINT64
                push!(wtypes, ":$fldname => :sint64")
            elseif field._type == FieldDescriptorProto_Type.TYPE_FIXED32
                push!(wtypes, ":$fldname => :fixed32")
            elseif field._type == FieldDescriptorProto_Type.TYPE_SFIXED32
                push!(wtypes, ":$fldname => :sfixed32")
            elseif field._type == FieldDescriptorProto_Type.TYPE_FIXED64
                push!(wtypes, ":$fldname => :fixed64")
            elseif field._type == FieldDescriptorProto_Type.TYPE_SFIXED64
                push!(wtypes, ":$fldname => :sfixed64")
            end
            enum_typ_name = (field._type == FieldDescriptorProto_Type.TYPE_ENUM) ? typ_name : ""
            (field._type != FieldDescriptorProto_Type.TYPE_MESSAGE) && (typ_name = "$(JTYPES[field._type])")

            push!(fldnums, field.number)
            (FieldDescriptorProto_Label.LABEL_REQUIRED == field.label) && push!(reqflds, ":"*fldname)

            if isfilled(field, :default_value) && !isempty(field.default_value)
                if field._type == FieldDescriptorProto_Type.TYPE_STRING
                    push!(defvals, ":$fldname => \"$(escape_string(field.default_value))\"")
                elseif field._type == FieldDescriptorProto_Type.TYPE_MESSAGE
                    println(errio, "Default values for message types are not supported. Field: $(dtypename).$(fldname) has default value [$(field.default_value)]")
                    return
                elseif field._type == FieldDescriptorProto_Type.TYPE_BYTES
                    println(errio, "Default values for byte array types are not supported. Field: $(dtypename).$(fldname) has default value [$(field.default_value)]")
                    return
                else
                    defval = (field._type == FieldDescriptorProto_Type.TYPE_ENUM) ? "$(short_type_name(enum_typ_name, depends)).$(field.default_value)" : "$(field.default_value)"
                    push!(defvals, ":$fldname => $defval")
                end
            end

            packed = (isfilled(field, :options) && field.options.packed) ||
                     ((syntax == "proto3") && (FieldDescriptorProto_Label.LABEL_REPEATED == field.label) && isprimitive(field._type))
            packed && push!(packedflds, ":"*fldname)

            if !(isresolved(dtypename, typ_name, full_typ_name, exports) || full_typ_name in _all_resolved)
                defer(full_dtypename, io, full_typ_name)
            end

            typ_name = short_type_name(typ_name, depends)
            is_typ_mapentry = typ_name in keys(mapentries)
            if is_typ_mapentry && !_map_as_array
                k,v = mapentries[typ_name]
                typ_name = "Dict{$k,$v}"
            elseif FieldDescriptorProto_Label.LABEL_REPEATED == field.label
                typ_name = "Array{$typ_name,1}"
            end
            println(io, "    $(fldname)::$typ_name", is_typ_mapentry ? " # map entry" : "")

            if ismapentry
                if field.number == 1
                    map_key_type = typ_name
                elseif field.number == 2
                    map_val_type = typ_name
                end
            end
        end
    end
    println(io, "    $(dtypename)(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)")
    println(io, "end #type $(dtypename)", ismapentry ? " (mapentry)" : "")

    # generate the meta for this type if required
    _d_fldnums = [1:length(fldnums);]
    !isempty(reqflds) && println(io, "const __req_$(dtypename) = Symbol[$(join(reqflds, ','))]")
    !isempty(defvals) && println(io, "const __val_$(dtypename) = @compat Dict($(join(defvals, ", ")))")
    (fldnums != _d_fldnums) && println(io, "const __fnum_$(dtypename) = Int[$(join(fldnums, ','))]")
    !isempty(packedflds) && println(io, "const __pack_$(dtypename) = Symbol[$(join(packedflds, ','))]")
    !isempty(wtypes) && println(io, "const __wtype_$(dtypename) = @compat Dict($(join(wtypes, ", ")))")
    if !isempty(oneofs)
        println(io, "const __oneofs_$(dtypename) = Int[$(join(oneofs, ','))]")
        println(io, "const __oneof_names_$(dtypename) = [$(join(oneof_names, ','))]")
    end
    if !isempty(reqflds) || !isempty(defvals) || (fldnums != [1:length(fldnums);]) || !isempty(packedflds) || !isempty(wtypes) || !isempty(oneofs)
        @logmsg("generating meta for type $(dtypename)")
        print(io, "meta(t::Type{$dtypename}) = meta(t, ")
        print(io, isempty(reqflds) ? "ProtoBuf.DEF_REQ, " : "__req_$(dtypename), ")
        print(io, (fldnums == _d_fldnums) ? "ProtoBuf.DEF_FNUM, " : "__fnum_$(dtypename), ")
        print(io, isempty(defvals) ? "ProtoBuf.DEF_VAL, " : "__val_$(dtypename), ")
        print(io, "true, ")
        print(io, isempty(packedflds) ? "ProtoBuf.DEF_PACK, " : "__pack_$(dtypename), ")
        print(io, isempty(wtypes) ? "ProtoBuf.DEF_WTYPES, " : "__wtype_$(dtypename), ")
        print(io, isempty(oneofs) ? "ProtoBuf.DEF_ONEOFS, " : "__oneofs_$(dtypename), ")
        print(io, isempty(oneofs) ? "ProtoBuf.DEF_ONEOF_NAMES" : "__oneof_names_$(dtypename)")
        println(io, ")")
    end
    # generate hash, equality and isequal methods
    println(io, "hash(v::$(dtypename)) = ProtoBuf.protohash(v)")
    println(io, "isequal(v1::$(dtypename), v2::$(dtypename)) = ProtoBuf.protoisequal(v1, v2)")
    println(io, "==(v1::$(dtypename), v2::$(dtypename)) = ProtoBuf.protoeq(v1, v2)")

    println(io, "")
    push!(exports, dtypename)
    ismapentry && (mapentries[dtypename] = (map_key_type, map_val_type))

    if !isdeferred(full_dtypename)
        @logmsg("resolved $full_dtypename")
        print(outio, takebuf_string(io))
        resolve(outio, full_dtypename)
        push!(_all_resolved, full_dtypename)
    end

    @logmsg("end type $(full_dtypename)")
    nothing
end

function protofile_name_to_module_name(n::AbstractString)
    name = splitext(basename(n))[1]
    name = replace(name, '.', '_')
    name = string(name, "_pb")
    return name
end

function has_gen_services(opt::FileOptions)
    isfilled(opt, :cc_generic_services) && opt.cc_generic_services && return true
    isfilled(opt, :py_generic_services) && opt.py_generic_services && return true
    isfilled(opt, :java_generic_services) && opt.java_generic_services && return true
    return false
end

function generate(io::IO, errio::IO, stype::ServiceDescriptorProto, svcidx::Int, exports::Array{AbstractString,1})
    nmethods = isfilled(stype, :method) ? length(stype.method) : 0

    # generate method and service descriptors
    println(io, "# service methods for $(stype.name)")
    println(io, "const _$(stype.name)_methods = MethodDescriptor[")
    for idx in 1:nmethods
        method = stype.method[idx]
        sm = splitmodule(method.input_type)
        _modul,in_typ_name = (length(sm) > 1) ? (sm[1],sm[2]) : ("",method.input_type)
        sm = splitmodule(method.output_type)
        _modul,out_typ_name = (length(sm) > 1) ? (sm[1],sm[2]) : ("",method.output_type)
        elem_sep = (idx < nmethods) ? "," : ""
        println(io, "        MethodDescriptor(\"$(method.name)\", $(idx), $(in_typ_name), $(out_typ_name))$(elem_sep)")
    end
    println(io, "    ] # const _$(stype.name)_methods")
    println(io, "const _$(stype.name)_desc = ServiceDescriptor(\"$(stype.name)\", $(svcidx), _$(stype.name)_methods)")
    println(io, "")

    # generate service
    println(io, "$(stype.name)(impl::Module) = ProtoService(_$(stype.name)_desc, impl)")
    println(io, "")
    push!(exports, stype.name)

    # generate stubs
    stub = "$(stype.name)Stub"
    println(io, "type $(stub) <: AbstractProtoServiceStub{false}")
    println(io, "    impl::ProtoServiceStub")
    println(io, "    $(stub)(channel::ProtoRpcChannel) = new(ProtoServiceStub(_$(stype.name)_desc, channel))")
    println(io, "end # type $(stub)")
    println(io, "")
    push!(exports, stub)

    nbstub = "$(stype.name)BlockingStub"
    println(io, "type $(nbstub) <: AbstractProtoServiceStub{true}")
    println(io, "    impl::ProtoServiceBlockingStub")
    println(io, "    $(nbstub)(channel::ProtoRpcChannel) = new(ProtoServiceBlockingStub(_$(stype.name)_desc, channel))")
    println(io, "end # type $(nbstub)")
    println(io, "")
    push!(exports, nbstub)

    for idx in 1:nmethods
        method = stype.method[idx]
        sm = splitmodule(method.input_type)
        _modul,in_typ_name = (length(sm) > 1) ? (sm[1],sm[2]) : ("",method.input_type)
        println(io, "$(method.name)(stub::$(stub), controller::ProtoRpcController, inp::$(in_typ_name), done::Function) = call_method(stub.impl, _$(stype.name)_methods[$(idx)], controller, inp, done)")
        println(io, "$(method.name)(stub::$(nbstub), controller::ProtoRpcController, inp::$(in_typ_name)) = call_method(stub.impl, _$(stype.name)_methods[$(idx)], controller, inp)")
        println(io, "")
        push!(exports, method.name)
    end
    nothing
end

function generate(io::IO, errio::IO, protofile::FileDescriptorProto)
    @logmsg("generate begin for $(protofile.name), package $(protofile.package)")

    svcs = isfilled(protofile, :options) ? has_gen_services(protofile.options) : false
    @logmsg("generate services: $svcs")

    scope = top_scope
    if !isempty(protofile.package)
        for pkgname in split(protofile.package, '.')
            scope = get_module_scope(scope, pkgname)
        end
    end

    push!(scope.files, protofile.name)
    @logmsg("generated scope for $(protofile.name), package $(protofile.package)")

    # generate module begin
    if !isempty(scope.name)
        scope.is_module = true
        _packages[protofile.name] = fullname(scope)
    end

    # generate syntax version
    syntax = (isfilled(protofile, :syntax) && !isempty(protofile.syntax)) ? protofile.syntax : "proto2"
    println(io, "# syntax: $(syntax)")

    depends = AbstractString[]
    @logmsg("generating imports")
    # generate imports
    println(io, "using Compat")
    println(io, "using ProtoBuf")
    println(io, "import ProtoBuf.meta")
    println(io, "import Base: hash, isequal, ==")
    if isfilled(protofile, :dependency)
        protofile_imports[protofile.name] = protofile.dependency
        using_pkgs = Set{AbstractString}()
        for dependency in protofile.dependency
            if haskey(_packages, dependency)
                push!(using_pkgs, _packages[dependency])
                push!(depends, _packages[dependency])
            end
        end

        fullscopename = scope.is_module ? fullname(scope) : ""
        parentscope = (isdefined(scope, :parent) && scope.parent.is_module) ? fullname(scope.parent) : ""
        for dependency in using_pkgs
            (fullscopename == dependency) && continue
            !isempty(parentscope) && startswith(dependency, parentscope) && (dependency = ".$(dependency[length(parentscope)+1:end])")
            println(io, "using $dependency")
        end
    end
    println(io, "")

    exports = AbstractString[]
    mapentries = Dict{AbstractString, Tuple{AbstractString,AbstractString}}()

    # generate top level enums
    @logmsg("generating enums")
    if isfilled(protofile, :enum_type)
        for enum_type in protofile.enum_type
            generate(io, errio, enum_type, scope, exports)
            (errio.size > 0) && return
        end
    end

    # generate message types
    @logmsg("generating types")
    if isfilled(protofile, :message_type)
        for message_type in protofile.message_type
            generate(io, errio, message_type, scope, syntax, exports, depends, mapentries)
            (errio.size > 0) && return
        end
    end

    # generate service stubs
    @logmsg("generating services")
    if svcs && isfilled(protofile, :service)
        nservices = length(protofile.service)
        for idx in 1:nservices
            service = protofile.service[idx]
            generate(io, errio, service, idx, exports)
            (errio.size > 0) && return
        end
    end

    # check if everything is resolved
    if !isempty(_deferred)
        println(errio, "All types could not be fully resolved:")
        for (k,v) in _deferred
            println(errio, "$k depends on $(join(v.depends, ','))")
        end
        return
    end

    # generate exports
    !isempty(exports) && println(io, "export ", join(exports, ", "))

    # mention mapentries
    !isempty(mapentries) && println(io, "# mapentries: ", join(mapentries, ", "))

    @logmsg("generate end for $(protofile.name)")
    nothing
end

function append_response(resp::CodeGeneratorResponse, protofile::FileDescriptorProto, io::IOBuffer)
    outdir = dirname(protofile.name)
    filename = protofile_name_to_module_name(protofile.name)
    filename = string(filename, ".jl")

    append_response(resp, filename, io)
end

function append_response(resp::CodeGeneratorResponse, filename::AbstractString, io::IOBuffer)
    jfile = ProtoBuf.instantiate(CodeGeneratorResponse_File)

    jfile.name = filename
    jfile.content = takebuf_string(io)

    !isdefined(resp, :file) && (resp.file = CodeGeneratorResponse_File[])
    push!(resp.file, jfile)
    resp
end

function err_response(errio::IOBuffer)
    resp = ProtoBuf.instantiate(CodeGeneratorResponse)
    resp.error = takebuf_string(errio)
    resp
end

scope_has_file(s::Scope, f::AbstractString) = (f in s.files || any(x->scope_has_file(x,f), s.children))

function print_package(io::IO, s::Scope, indent="")
    s.is_module || return
    println(io, "$(indent)module $(s.name)")
    nested = indent*"  "
    children = Set(s.children)
    for f in s.files
        if haskey(protofile_imports,f)
            deps = protofile_imports[f]
            # if the file we're about to `include` depends on any of our child modules,
            # generate those first.
            for d in deps
                for c in filter(c->scope_has_file(c,d), children)
                    delete!(children, c)
                    print_package(io, c, nested)
                end
            end
        end

        fname = string(protofile_name_to_module_name(f), ".jl")
        println(io, "$(nested)include(\"$fname\")")
    end
    for c in s.children
        # check if already included
        (c in children) || continue
        print_package(io, c, nested)
    end
    println(io, "$(indent)end")
end

function codegen(srcio::IO)
    errio = IOBuffer()
    resp = ProtoBuf.instantiate(CodeGeneratorResponse)
    @logmsg("generate begin")
    while !eof(srcio)
        req = readreq(srcio)

        if !isfilled(req, :file_to_generate)
            @logmsg("no files to generate!!")
            continue
        end

        @logmsg("generate request for $(length(req.file_to_generate)) proto files")
        @logmsg("$(req.file_to_generate)")

        #isfilled(req, :parameter) && @logmsg("parameter $(req.parameter)")

        for protofile in req.proto_file
            io = IOBuffer()
            generate(io, errio, protofile)
            (errio.size > 0) && return err_response(errio)
            append_response(resp, protofile, io)
        end
    end
    if !isempty(top_scope.children)
        for pkg in top_scope.children
            if pkg.is_module
                io = IOBuffer()
                print_package(io, pkg)
                append_response(resp, "$(pkg.name).jl", io)
            end
        end
    end
    @logmsg("generate end")
    resp
end

##
# the main read - write method
function gen()
    try
        global _module_postfix = in("--module-postfix-enabled", ARGS)
        global _map_as_array = in("--map-as-array", ARGS)
        writeproto(STDOUT, codegen(STDIN))
    catch ex
        println(STDERR, "Exception while generating Julia code")
        rethrow()
    end
end

end # module Gen


text:
#
# Date created: 2022-03-08
# Author: aradclif
#
#
############################################################################################
# Cleaned version of experimental.jl
using Static
using LoopVectorization

function loopgen(N::Int) #A::AbstractArray{T, N} where {T, N}
    if N == 1
        loops = Expr(:for)
        ex = Expr(:(=), Symbol(:i_, 1), Expr(:call, :axes, :A, 1))
        push!(loops.args, ex)
    else
        loops = Expr(:for)
        block = Expr(:block)
        for d = N:-1:1
            ex = Expr(:(=), Symbol(:i_, d), Expr(:call, :axes, :A, d))
            push!(block.args, ex)
        end
        push!(loops.args, block)
    end
    return loops
end
loopgen(A::AbstractArray{T, N}) where {T, N} = loopgen(N)

function reducebody(f, N, D)
    body = Expr(:block)
    a = Expr(:ref, :A, ntuple(d -> Symbol(:i_, d), N)...)
    b = Expr(:ref, :B, ntuple(d -> D[d] == 1 ? 1 : Symbol(:i_, d), N)...)
    e = Expr(:(=), b, Expr(:call, Symbol(f), b, a))
    push!(body.args, e)
    body
end
function mapreducebody(f, op, N, D)
    body = Expr(:block)
    a = Expr(:ref, :A, ntuple(d -> Symbol(:i_, d), N)...)
    b = Expr(:ref, :B, ntuple(d -> D[d] == 1 ? 1 : Symbol(:i_, d), N)...)
    e = Expr(:(=), b, Expr(:call, Symbol(op), b, Expr(:call, Symbol(f), a)))
    push!(body.args, e)
    body
end

#### quote for use in @generated
function reduce_quote(F, N::Int, D)
    ls = loopgen(N)
    body = Expr(:block)
    params = D.parameters
    f = F.instance
    a = Expr(:ref, :A, ntuple(d -> Symbol(:i_, d), N)...)
    b = Expr(:ref, :B, ntuple(d -> params[d] == Static.One ? 1 : Symbol(:i_, d), N)...)
    e = Expr(:(=), b, Expr(:call, Symbol(f), b, a))
    push!(body.args, e)
    push!(ls.args, body)
    return quote
        @turbo $ls
        return B
    end
end

function lvreduce(f, A::AbstractArray{T, N}, dims::NTuple{M, Int}) where {T, N, M}
    if ntuple(identity, Val(N)) ⊆ dims
        B = hvncat(ntuple(i -> 1, Val(N)), true, lvreduce1(f, A))
    else
        Dᴬ = size(A)
        Dᴮ = ntuple(d -> d ∈ dims ? 1 : Dᴬ[d], N)
        B = zeros(Base.promote_op(f, T), Dᴮ)
        # Dᴮ′ = ntuple(d -> StaticInt(Dᴮ[d]), N)
        Dᴮ′ = ntuple(d -> d ∈ dims ? StaticInt(1) : Dᴮ[d], N)
        _lvreduce!(f, B, A, Dᴮ′)
    end
    return B
end

@generated function _lvreduce!(f::F, B::AbstractArray{Tₒ, N}, A::AbstractArray{T, N}, dims::D) where {F, Tₒ, T, N, D}
    reduce_quote(F, N, D)
end

################
# Handle scalar dims by wrapping in Tuple
lvreduce(f, A::AbstractArray{T, N}, dims::Int) where {T, N} = lvreduce(f, A, (dims,))
# Convenience dispatches to match JuliaBase
lvreduce(f, A::AbstractArray{T, N}; dims=:) where {T, N} = lvreduce(f, A, dims)
lvreduce(f, A::AbstractArray{T, N}) where {T, N} = lvreduce1(f, A)
lvreduce(f, A::AbstractArray{T, N}, ::Colon) where {T, N} = lvreduce1(f, A)

# When dimensions unspecified, treat as vector. Aside from special cases,
# one must usually provide an initial value in order for the reduce to be defined.
# In lieu of a proper approach, zero seems to be a decent guess.
@generated function lvreduce1(f::F, A::AbstractArray{T, N}) where {F, T, N}
    f = F.instance
    quote
        s = zero($T)
        @turbo for i ∈ eachindex(A)
            s = $f(s, A[i])
        end
        s
    end
end
for (op, init) ∈ zip((:+, :-, :*, :max, :min), (:zero, :zero, :one, :typemin, :typemax))
    @eval function lvreduce(::typeof($op), A::AbstractArray{T, N}, dims::NTuple{M, Int}) where {T, N, M}
        if ntuple(identity, Val(N)) ⊆ dims
            B = hvncat(ntuple(i -> 1, Val(N)), true, lvreduce1($op, A))
        else
            Dᴬ = size(A)
            Dᴮ = ntuple(d -> d ∈ dims ? 1 : Dᴬ[d], N)
            B = fill($init(T), Dᴮ)
            # Dᴮ′ = ntuple(d -> StaticInt(Dᴮ[d]), N)
            Dᴮ′ = ntuple(d -> d ∈ dims ? StaticInt(1) : Dᴮ[d], N)
            _lvreduce!($op, B, A, Dᴮ′)
        end
        return B
    end
    @eval function lvreduce1(::typeof($op), A::AbstractArray{T, N}) where {T, N}
        s = $init(T)
        @turbo for i ∈ eachindex(A)
            s = $op(s, A[i])
        end
        s
    end
end

# Convenience definitions
lvsum(A::AbstractArray{T, N}; dims=:) where {T, N} = lvreduce(+, A, dims=dims)
lvprod(A::AbstractArray{T, N}; dims=:) where {T, N} = lvreduce(*, A, dims=dims)
lvmaximum(A::AbstractArray{T, N}; dims=:) where {T, N} = lvreduce(max, A, dims=dims)
lvminimum(A::AbstractArray{T, N}; dims=:) where {T, N} = lvreduce(min, A, dims=dims)
lvextrema(A::AbstractArray{T, N}; dims=:) where {T, N} =
    collect(zip(lvminimum(A, dims=dims), lvmaximum(A, dims=dims)))
################

function mapreduce_quote(F, OP, N::Int, D)
    ls = loopgen(N)
    body = Expr(:block)
    params = D.parameters
    f = F.instance
    op = OP.instance
    a = Expr(:ref, :A, ntuple(d -> Symbol(:i_, d), N)...)
    b = Expr(:ref, :B, ntuple(d -> params[d] == Static.One ? 1 : Symbol(:i_, d), N)...)
    e = Expr(:(=), b, Expr(:call, Symbol(op), b, Expr(:call, Symbol(f), a)))
    push!(body.args, e)
    push!(ls.args, body)
    return quote
        @turbo $ls
        return B
    end
end

function lvmapreduce(f, op, A::AbstractArray{T, N}, dims::NTuple{M, Int}) where {T, N, M}
    if ntuple(identity, Val(N)) ⊆ dims
        B = hvncat(ntuple(i -> 1, Val(N)), true, lvmapreduce1(f, op, A))
    else
        Dᴬ = size(A)
        Dᴮ = ntuple(d -> d ∈ dims ? 1 : Dᴬ[d], N)
        B = zeros(Base.promote_op(f, T), Dᴮ)
        # Dᴮ′ = ntuple(d -> StaticInt(Dᴮ[d]), N)
        Dᴮ′ = ntuple(d -> d ∈ dims ? StaticInt(1) : Dᴮ[d], N)
        _lvmapreduce!(f, op, B, A, Dᴮ′)
    end
    return B
end

@generated function _lvmapreduce!(f::F, op::OP, B::AbstractArray{Tₒ, N}, A::AbstractArray{T, N}, dims::D) where {F, OP, Tₒ, T, N, D}
    mapreduce_quote(F, OP, N, D)
end

################
# Handle scalar dims by wrapping in Tuple
lvmapreduce(f, op, A::AbstractArray{T, N}, dims::Int) where {T, N} = lvmapreduce(f, op, A, (dims,))
# Convenience dispatches to match JuliaBase
lvmapreduce(f, op, A::AbstractArray{T, N}; dims=:) where {T, N} = lvmapreduce(f, op, A, dims)
lvmapreduce(f, op, A::AbstractArray{T, N}) where {T, N} = lvmapreduce1(f, op, A)
lvmapreduce(f, op, A::AbstractArray{T, N}, ::Colon) where {T, N} = lvmapreduce1(f, op, A)

# When dimensions unspecified, treat as vector. Aside from special cases,
# one must usually provide an initial value in order for the reduce to be defined.
# In lieu of a proper approach, zero seems to be a decent guess.
@generated function lvmapreduce1(f::F, op::OP, A::AbstractArray{T, N}) where {F, OP, T, N}
    op = OP.instance
    f = F.instance
    quote
        s = zero(Base.promote_op($f, $T))
        @turbo for i ∈ eachindex(A)
            s = $op(s, $f(A[i]))
        end
        s
    end
end
for (op, init) ∈ zip((:+, :-, :*, :max, :min), (:zero, :zero, :one, :typemin, :typemax))
    @eval function lvmapreduce(f, ::typeof($op), A::AbstractArray{T, N}, dims::NTuple{M, Int}) where {T, N, M}
        if ntuple(identity, Val(N)) ⊆ dims
            B = hvncat(ntuple(i -> 1, Val(N)), true, lvmapreduce1(f, $op, A))
        else
            Dᴬ = size(A)
            Dᴮ = ntuple(d -> d ∈ dims ? 1 : Dᴬ[d], N)
            B = fill($init(Base.promote_op(f, T)), Dᴮ)
            # Dᴮ′ = ntuple(d -> StaticInt(Dᴮ[d]), N)
            Dᴮ′ = ntuple(d -> d ∈ dims ? StaticInt(1) : Dᴮ[d], N)
            _lvmapreduce!(f, $op, B, A, Dᴮ′)
        end
        return B
    end
    @eval function lvmapreduce1(f, ::typeof($op), A::AbstractArray{T, N}) where {T, N}
        s = $init(Base.promote_op(f, T))
        @turbo for i ∈ eachindex(A)
            s = $op(s, f(A[i]))
        end
        s
    end
end

# Convenience definitions
lvsum(f, A::AbstractArray{T, N}; dims=:) where {T, N} = lvmapreduce(f, +, A, dims=dims)
lvprod(f, A::AbstractArray{T, N}; dims=:) where {T, N} = lvmapreduce(f, *, A, dims=dims)
lvmaximum(f, A::AbstractArray{T, N}; dims=:) where {T, N} = lvmapreduce(f, max, A, dims=dims)
lvminimum(f, A::AbstractArray{T, N}; dims=:) where {T, N} = lvmapreduce(f, min, A, dims=dims)
lvextrema(f, A::AbstractArray{T, N}; dims=:) where {T, N} =
    collect(zip(lvminimum(f, A, dims=dims), lvmaximum(f, A, dims=dims)))
############################################################################################
#### threaded versions
function treduce_quote(F, N::Int, D)
    ls = loopgen(N)
    body = Expr(:block)
    params = D.parameters
    f = F.instance
    a = Expr(:ref, :A, ntuple(d -> Symbol(:i_, d), N)...)
    b = Expr(:ref, :B, ntuple(d -> params[d] == Static.One ? 1 : Symbol(:i_, d), N)...)
    e = Expr(:(=), b, Expr(:call, Symbol(f), b, a))
    push!(body.args, e)
    push!(ls.args, body)
    return quote
        @tturbo $ls
        return B
    end
end
function lvtreduce(f, A::AbstractArray{T, N}, dims::NTuple{M, Int}) where {T, N, M}
    if ntuple(identity, Val(N)) ⊆ dims
        B = hvncat(ntuple(i -> 1, Val(N)), true, lvtreduce1(f, A))
    else
        Dᴬ = size(A)
        Dᴮ = ntuple(d -> d ∈ dims ? 1 : Dᴬ[d], N)
        B = zeros(Base.promote_op(f, T), Dᴮ)
        # Dᴮ′ = ntuple(d -> StaticInt(Dᴮ[d]), N)
        Dᴮ′ = ntuple(d -> d ∈ dims ? StaticInt(1) : Dᴮ[d], N)
        _lvtreduce!(f, B, A, Dᴮ′)
    end
    return B
end
@generated function _lvtreduce!(f::F, B::AbstractArray{T, N}, A::AbstractArray{T, N}, dims::D) where {F, T, N, D}
    treduce_quote(F, N, D)
end

################
# Handle scalar dims by wrapping in Tuple
lvtreduce(f, A::AbstractArray{T, N}, dims::Int) where {T, N} = lvtreduce(f, A, (dims,))
# Convenience dispatches to match JuliaBase
lvtreduce(f, A::AbstractArray{T, N}; dims=:) where {T, N} = lvtreduce(f, A, dims)
lvtreduce(f, A::AbstractArray{T, N}) where {T, N} = lvtreduce1(f, A)
lvtreduce(f, A::AbstractArray{T, N}, ::Colon) where {T, N} = lvtreduce1(f, A)

# When dimensions unspecified, treat as vector. Aside from special cases,
# one must usually provide an initial value in order for the reduce to be defined.
# In lieu of a proper approach, zero seems to be a decent guess.
@generated function lvtreduce1(f::F, A::AbstractArray{T, N}) where {F, T, N}
    f = F.instance
    quote
        s = zero($T)
        @tturbo for i ∈ eachindex(A)
            s = $f(s, A[i])
        end
        s
    end
end
for (op, init) ∈ zip((:+, :-, :*, :max, :min), (:zero, :zero, :one, :typemin, :typemax))
    @eval function lvtreduce(::typeof($op), A::AbstractArray{T, N}, dims::NTuple{M, Int}) where {T, N, M}
        if ntuple(identity, Val(N)) ⊆ dims
            B = hvncat(ntuple(i -> 1, Val(N)), true, lvtreduce1($op, A))
        else
            Dᴬ = size(A)
            Dᴮ = ntuple(d -> d ∈ dims ? 1 : Dᴬ[d], N)
            B = fill($init(T), Dᴮ)
            # Dᴮ′ = ntuple(d -> StaticInt(Dᴮ[d]), N)
            Dᴮ′ = ntuple(d -> d ∈ dims ? StaticInt(1) : Dᴮ[d], N)
            _lvtreduce!($op, B, A, Dᴮ′)
        end
        return B
    end
    @eval function lvtreduce1(::typeof($op), A::AbstractArray{T, N}) where {T, N}
        s = $init(T)
        @tturbo for i ∈ eachindex(A)
            s = $op(s, A[i])
        end
        s
    end
end

# Convenience definitions
lvtsum(A::AbstractArray{T, N}; dims=:) where {T, N} = lvtreduce(+, A, dims=dims)
lvtprod(A::AbstractArray{T, N}; dims=:) where {T, N} = lvtreduce(*, A, dims=dims)
lvtmaximum(A::AbstractArray{T, N}; dims=:) where {T, N} = lvtreduce(max, A, dims=dims)
lvtminimum(A::AbstractArray{T, N}; dims=:) where {T, N} = lvtreduce(min, A, dims=dims)
lvtextrema(A::AbstractArray{T, N}; dims=:) where {T, N} =
    collect(zip(lvtminimum(A, dims=dims), lvtmaximum(A, dims=dims)))
################

function tmapreduce_quote(F, OP, N::Int, D)
    ls = loopgen(N)
    body = Expr(:block)
    params = D.parameters
    f = F.instance
    op = OP.instance
    a = Expr(:ref, :A, ntuple(d -> Symbol(:i_, d), N)...)
    b = Expr(:ref, :B, ntuple(d -> params[d] == Static.One ? 1 : Symbol(:i_, d), N)...)
    e = Expr(:(=), b, Expr(:call, Symbol(op), b, Expr(:call, Symbol(f), a)))
    push!(body.args, e)
    push!(ls.args, body)
    return quote
        @tturbo $ls
        return B
    end
end
function lvtmapreduce(f, op, A::AbstractArray{T, N}, dims::NTuple{M, Int}) where {T, N, M}
    if ntuple(identity, Val(N)) ⊆ dims
        B = hvncat(ntuple(i -> 1, Val(N)), true, lvtmapreduce1(f, op, A))
    else
        Dᴬ = size(A)
        Dᴮ = ntuple(d -> d ∈ dims ? 1 : Dᴬ[d], N)
        B = zeros(Base.promote_op(f, T), Dᴮ)
        # Dᴮ′ = ntuple(d -> StaticInt(Dᴮ[d]), N)
        Dᴮ′ = ntuple(d -> d ∈ dims ? StaticInt(1) : Dᴮ[d], N)
        _lvtmapreduce!(f, op, B, A, Dᴮ′)
    end
    return B
end

@generated function _lvtmapreduce!(f::F, op::OP, B::AbstractArray{T, N}, A::AbstractArray{T, N}, dims::D) where {F, OP, T, N, D}
    tmapreduce_quote(F, OP, N, D)
end

################
# Handle scalar dims by wrapping in Tuple
lvtmapreduce(f, op, A::AbstractArray{T, N}, dims::Int) where {T, N} = lvtmapreduce(f, op, A, (dims,))
# Convenience dispatches to match JuliaBase
lvtmapreduce(f, op, A::AbstractArray{T, N}; dims=:) where {T, N} = lvtmapreduce(f, op, A, dims)
lvtmapreduce(f, op, A::AbstractArray{T, N}) where {T, N} = lvtmapreduce1(f, op, A)
lvtmapreduce(f, op, A::AbstractArray{T, N}, ::Colon) where {T, N} = lvtmapreduce1(f, op, A)

# When dimensions unspecified, treat as vector. Aside from special cases,
# one must usually provide an initial value in order for the reduce to be defined.
# In lieu of a proper approach, zero seems to be a decent guess.
@generated function lvtmapreduce1(f::F, op::OP, A::AbstractArray{T, N}) where {F, OP, T, N}
    f = F.instance
    op = OP.instance
    quote
        s = zero(Base.promote_op($f, $T))
        @tturbo for i ∈ eachindex(A)
            s = $op(s, $f(A[i]))
        end
        s
    end
end
for (op, init) ∈ zip((:+, :-, :*, :max, :min), (:zero, :zero, :one, :typemin, :typemax))
    @eval function lvtmapreduce(f, ::typeof($op), A::AbstractArray{T, N}, dims::NTuple{M, Int}) where {T, N, M}
        if ntuple(identity, Val(N)) ⊆ dims
            B = hvncat(ntuple(i -> 1, Val(N)), true, lvtmapreduce1(f, $op, A))
        else
            Dᴬ = size(A)
            Dᴮ = ntuple(d -> d ∈ dims ? 1 : Dᴬ[d], N)
            B = fill($init(Base.promote_op(f, T)), Dᴮ)
            # Dᴮ′ = ntuple(d -> StaticInt(Dᴮ[d]), N)
            Dᴮ′ = ntuple(d -> d ∈ dims ? StaticInt(1) : Dᴮ[d], N)
            _lvtmapreduce!(f, $op, B, A, Dᴮ′)
        end
        return B
    end
    @eval function lvtmapreduce1(f, ::typeof($op), A::AbstractArray{T, N}) where {T, N}
        s = $init(Base.promote_op(f, T))
        @tturbo for i ∈ eachindex(A)
            s = $op(s, f(A[i]))
        end
        s
    end
end

# Convenience definitions
lvtsum(f, A::AbstractArray{T, N}; dims=:) where {T, N} = lvtmapreduce(f, +, A, dims=dims)
lvtprod(f, A::AbstractArray{T, N}; dims=:) where {T, N} = lvtmapreduce(f, *, A, dims=dims)
lvtmaximum(f, A::AbstractArray{T, N}; dims=:) where {T, N} = lvtmapreduce(f, max, A, dims=dims)
lvtminimum(f, A::AbstractArray{T, N}; dims=:) where {T, N} = lvtmapreduce(f, min, A, dims=dims)
lvtextrema(f, A::AbstractArray{T, N}; dims=:) where {T, N} =
    collect(zip(lvtminimum(f, A, dims=dims), lvtmaximum(f, A, dims=dims)))
################

ext: jl


text:
using PastaQ
using PastaQ.ITensors
using Test
using LinearAlgebra
using Random
using Observers
using JLD2
using HDF5
using Optimisers

@testset "write and read samples" begin
  N = 3
  depth = 4
  nshots = 100
  circuit = randomcircuit(N; depth = depth)
  path = "test_data_writesamples.h5"

  X = runcircuit(circuit)
  data = getsamples(X, nshots)
  writesamples(data, path)
  datatest = readsamples(path)
  @test datatest ≈ data
  writesamples(data, X, path)
  datatest, Xtest = readsamples(path)
  @test X ≈ Xtest
  @test datatest ≈ data

  X = runcircuit(circuit; noise = ("DEP",(p=0.01,)))
  data = getsamples(X, nshots)
  writesamples(data, path)
  datatest = readsamples(path)
  @test datatest ≈ data
  writesamples(data, X, path)
  datatest, Xtest = readsamples(path)
  @test X ≈ Xtest
  @test datatest ≈ data

  X = randomstate(N; ξ = 2, χ = 3)
  writesamples(data, X, path)
  datatest, Xtest = readsamples(path)
  @test X.X ≈ Xtest.X
  @test datatest ≈ data

  X = runcircuit(circuit)
  bases = randombases(N, 10)
  data = getsamples(X, bases, nshots)
  writesamples(data, path)
  datatest = readsamples(path)
  @test all(data .== datatest)
  writesamples(data, X, path)
  datatest, Xtest = readsamples(path)
  @test X ≈ Xtest
  @test all(data .== datatest)

  X = runcircuit(circuit; noise = ("DEP",(p=0.01,)))
  bases = randombases(N, 10)
  data = getsamples(X, bases, nshots)
  writesamples(data, path)
  datatest = readsamples(path)
  @test all(data .== datatest)
  writesamples(data, X, path)
  datatest, Xtest = readsamples(path)
  @test X ≈ Xtest
  @test all(data .== datatest)

  X = randomstate(N; ξ = 2, χ = 3)
  writesamples(data, X, path)
  datatest, Xtest = readsamples(path)
  @test X.X ≈ Xtest.X
  @test all(data .== datatest)

  X = runcircuit(circuit; process = true)
  preps = randompreparations(N, 8)
  bases = randombases(N, 10)
  data = getsamples(X, preps, bases, nshots)
  writesamples(data, path)
  datatest = readsamples(path)
  @test all(data .== datatest)
  writesamples(data, X, path)
  datatest, Xtest = readsamples(path)
  @test X ≈ Xtest
  @test all(data .== datatest)

  X = runcircuit(circuit; noise = ("DEP",(p=0.01,)), process = true)
  preps = randompreparations(N, 8)
  bases = randombases(N, 10)
  data = getsamples(X, preps, bases, nshots)
  writesamples(data, path)
  datatest = readsamples(path)
  @test all(data .== datatest)
  writesamples(data, X, path)
  datatest, Xtest = readsamples(path)
  @test all(data .== datatest)

  X = randomstate(N; ξ = 2, χ = 3)
  writesamples(data, X, path)
  datatest, Xtest = readsamples(path)
  @test X.X ≈ Xtest.X
  @test all(data .== datatest)
end

@testset "circuit observer: MPS" begin
  N = 8
  depth = 9
  R = 5
  Random.seed!(1234)
  circuit = randomcircuit(N; depth = depth)
  layer = Tuple[]
  push!(circuit, [("CX",(1,N)),])

  sites = siteinds("Qubit", N)

  outputpath = "simulation"
  ϕ = randomstate(sites; χ = 10, normalize=true)
  ψ = runcircuit(sites, circuit)
  Ftest = fidelity(ψ,ϕ)
  f(ψ::MPS; kwargs...) = fidelity(ψ, ϕ)#; kwargs...) = fidelity(ψ, ϕ)
  obs = Observer(["f" => f])
  ψ = runcircuit(sites, circuit; (observer!)=obs, move_sites_back_before_measurements=true,
                 outputpath = outputpath, savestate = true, outputlevel = 0)
  @test Ftest ≈ results(obs, "f")[end]
  @test length(results(obs, "f")) == depth+1

  obs2 = load("simulation_observer.jld2")
  for (k,v) in obs2
    @test last(v) ≈ results(obs, k)
  end

  fin = h5open("simulation_state.h5","r")
  M = read(fin, "state", MPS)
  close(fin)
  @test M ≈ ψ
end

@testset "circuit observer: MPO" begin
  N = 4
  depth = 5
  #circuit = Vector{Vector{<:Any}}(undef, depth)
  sites = siteinds("Qubit", N)
  circuit = randomcircuit(N; depth = depth)

  @disable_warn_order begin
    L = randomstate(sites; χ = 10, ξ = 3, normalize=true)
    ϱ = MPO(L)
    ρ = runcircuit(sites, circuit; noise = ("DEP",(p=0.001,)))
    Ftest = fidelity(ϱ,ρ)
    g(ρ::MPO; kwargs...) = fidelity(ρ, ϱ)#; kwargs...) = fidelity(ψ, ϕ)
    obs = Observer(["g" => g])
    outputpath = "simulation"
    ρ = runcircuit(sites, circuit;
                   (observer!)=obs,
                   move_sites_back_before_measurements=true,
                   outputlevel = 0,
                   noise = ("DEP",(p=0.001,)),
                   outputpath = outputpath, savestate = true)
 end
  @test Ftest ≈ results(obs, "g")[end]
  @test length(results(obs, "g")) == depth
  obs2 = load("simulation_observer.jld2")
  for (k,v) in obs2
    @test last(v) ≈ results(obs, k)
  end

  fin = h5open("simulation_state.h5","r")
  M = read(fin, "state", MPO)
  close(fin)
  @test M ≈ ρ
end

@testset "state tomography observer output: MPS" begin
  Random.seed!(1234)
  N = 4
  depth = 4
  nshots = 100
  circuit = randomcircuit(N; depth = depth)
  Ψ = runcircuit(circuit)
  bases = randombases(N,2)
  data = getsamples(Ψ, bases, nshots)
  test_data = copy(data[1:10, :])

  N = length(Ψ)     # Number of productstate
  χ = maxlinkdim(Ψ) # Bond dimension of variational MPS
  ψ0 = randomstate(Ψ; χ=χ, σ=0.1)
  opt = Optimisers.Descent(0.01)

  F(ψ::MPS; kwargs...) = fidelity(ψ, Ψ)
  obs = Observer(["F" => F])
  epochs = 18

  batchsize = 10
  observe_step = 3

  outputpath = "simulation"
  ψ = tomography(
    data,
    ψ0;
    test_data=test_data,
    batchsize=10,
    epochs=epochs,
    (observer!)=obs,
    observe_step = observe_step,
    outputpath = outputpath,
    savestate = true,
    outputlevel = 0
   )
  @test length(results(obs, "F")) == epochs ÷ observe_step

  obs2 = load("simulation_observer.jld2")
  for (k,v) in obs2
    @test last(v) ≈ results(obs, k)
  end

  fin = h5open("simulation_state.h5","r")
  M = read(fin, "state", MPS)
  close(fin)
  @test M ≈ ψ
end

@testset "state tomography observer output: LPDO" begin
  Random.seed!(1234)
  N = 4
  depth = 4
  nshots = 100
  circuit = randomcircuit(N; depth = depth)
  ϱ = runcircuit(circuit; noise = ("DEP",(p=0.01,)))
  bases = randombases(N,2)
  data = getsamples(ϱ, bases, nshots)
  test_data = copy(data[1:10, :])

  N = length(ϱ)     # Number of productstate
  χ = maxlinkdim(ϱ) # Bond dimension of variational MPS
  ρ = randomstate(ϱ; χ=χ÷2, ξ=2)
  opt = Optimisers.Descent(0.01)

  F(ρ::LPDO; kwargs...) = fidelity(ρ, ϱ)
  obs = Observer(["F" => F])
  epochs = 9

  batchsize = 10
  observe_step = 3

  outputpath = "simulation"
  ρ = tomography(
    data,
    ρ;
    test_data=test_data,
    batchsize=10,
    epochs=epochs,
    (observer!)=obs,
    observe_step = observe_step,
    outputpath = outputpath,
    savestate = true,
    outputlevel = 0
  )
  @test length(results(obs, "F")) == epochs ÷ observe_step

  obs2 = load("simulation_observer.jld2")
  for (k,v) in obs2
    @test last(v) ≈ results(obs, k)
  end

  fin = h5open("simulation_state.h5","r")
  M = read(fin, "state", LPDO{MPO})
  close(fin)
  @test M.X ≈ ρ.X
end

@testset "process tomography observer output: MPO" begin
  Random.seed!(1234)
  N = 3
  depth = 4
  nshots = 100
  circuit = randomcircuit(N; depth =  depth)
  V = runcircuit(circuit; process = true)
  preps = randompreparations(N,2)
  bases = randombases(N,2)
  data = getsamples(V, preps, bases, nshots)
  test_data = copy(data[1:10, :])
  N = length(V)     # Number of productstate
  χ = maxlinkdim(V) # Bond dimension of variational MPS
  U0 = randomprocess(V; χ=χ, σ=0.1)
  opt = Optimisers.Descent(0.01)

  F(U::MPO; kwargs...) = fidelity(U, V; process=true)
  obs = Observer(["F" => F])
  epochs = 9

  batchsize = 10
  observe_step = 3

  outputpath = "simulation"
  U = tomography(
    data,
    U0;
    test_data=test_data,
    batchsize=10,
    epochs=epochs,
    (observer!)=obs,
    observe_step = observe_step,
    outputpath = outputpath,
    savestate = true,
    outputlevel = 0
   )

  @test length(results(obs, "F")) == epochs ÷ observe_step

  obs2 = load("simulation_observer.jld2")
  for (k,v) in obs2
    @test last(v) ≈ results(obs, k)
  end
  fin = h5open("simulation_state.h5","r")
  M = read(fin, "state", MPO)
  close(fin)
  @test M ≈ U
end

@testset "state tomography observer output: LPDO" begin
  Random.seed!(1234)
  N = 2
  depth = 4
  nshots = 100
  circuit = randomcircuit(N; depth = depth)
  Φ = runcircuit(circuit; process = true, noise = ("DEP",(p=0.01,)))
  preps = randompreparations(N,2)
  bases = randombases(N,2)
  data = getsamples(Φ, preps, bases, nshots)
  test_data = copy(data[1:10, :])

  N = length(Φ)     # Number of productstate
  χ = maxlinkdim(Φ) # Bond dimension of variational MPS
  Λ = randomprocess(Φ; χ=χ÷2, ξ=2)
  opt = Optimisers.Descent(0.01)

  F(Λ::LPDO; kwargs...) = fidelity(Λ, Φ)
  obs = Observer(["F" => F])
  epochs = 9

  batchsize = 10
  observe_step = 3

  outputpath = "simulation"
  Λ = tomography(
    data,
    Λ;
    test_data=test_data,
    batchsize=10,
    epochs=epochs,
    (observer!)=obs,
    observe_step = observe_step,
    outputpath = outputpath,
    savestate = true,
    outputlevel = 0
  )
  @test length(results(obs, "F")) == epochs ÷ observe_step

  obs2 = load("simulation_observer.jld2")
  for (k,v) in obs2
    @test last(v) ≈ results(obs, k)
  end

  fin = h5open("simulation_state.h5","r")
  M = read(fin, "state", LPDO{MPO})
  close(fin)
  @test M.X ≈ Λ.X
end


text:
# Run thermal storage sensitivity cases
# Note: the penalty for unmet loads should be L1 and not L2 for this file to produce expected results.
# An L1 norm will cause the algorithm to detect unfeasibility faster (because the unmet loads will all be on the same hour).

push!(LOAD_PATH, joinpath(ENV["REPO_PATH"], "Sesi.jl/src"))

using CSV
using JuMP
using IterTools
using JSON
using DataFrames

import Sesi # import vs using so I can reload this from the REPL

reload("Sesi")

include("jl_utils.jl")

save = true
verb = true

if verb
    println(repeat("#",80))
    println("Starting storage sensitivity")
    println(repeat("#",80))
end
folderIn = joinpath(ENV["REPO_PATH"], "data/jl_in/")
folderOut = joinpath(ENV["REPO_PATH"], "data/jl_out/")

results = Dict()

fileNm = folderIn * "jlInput_x0.csv"

# Increase resolution as needed here
hot_lst=linspace(0., 1100., 30)
cold_lst=linspace(0., 120.e3, 30)

results["hot_lst"] = hot_lst
results["cold_lst"] = cold_lst
results["cost"] = zeros(length(hot_lst), length(cold_lst))
results["feas"] = zeros(length(hot_lst), length(cold_lst))
results["peak"] = zeros(length(hot_lst), length(cold_lst))
results["demandCharge"] = zeros(length(hot_lst), length(cold_lst))
results["agg_dollar"] = zeros(length(hot_lst), length(cold_lst))
results["chillers"] = zeros(length(hot_lst), length(cold_lst))
results["umHot"] = zeros(length(hot_lst), length(cold_lst))
results["umCold"] = zeros(length(hot_lst), length(cold_lst))
results["bill"] = Dict()

for (ihot, hot) in enumerate(hot_lst)
for (icold, cold) in enumerate(cold_lst)
    nChiller = 4.
    umCost = 1e7
    while umCost > 1e-2
        if verb
            @printf("hot: %.1f mmbtu - cold: %.1f tons - chillers %.0f\n", hot,
                cold, nChiller)
        end
        sesi, p, f = build_model(fileNm; verb=false, startDate=Date(2016,1,1),
            endDate=Date(2016,12,31), hot_storage=hot, cold_storage=cold,
            nChiller=nChiller, umLoadsL2=false)
        status = solve(sesi)

        if status == :Optimal
            bill = Sesi.campus_bill(sesi, f, p, verb=0)
            results["feas"][ihot,icold] = 1.
            results["peak"][ihot,icold] = maximum(getvalue(sesi[:peak]))
            results["demandCharge"][ihot,icold] = bill["demandCharge"]
            results["agg_dollar"][ihot,icold] = bill["agg_dollar"]
            results["bill"][(ihot,icold)] = bill
            results["chillers"][ihot,icold] = nChiller-4.
            if nChiller == 4.
                results["umHot"][ihot,icold] = sum(getvalue(sesi[:umHot]))
                results["umCold"][ihot,icold] = sum(getvalue(sesi[:umCold]))
            end
            umCost = bill["umLoads"]
            @printf("unmet Load Cost: %.6f\n", umCost)
        else
            warn("Unknown status!")
            println(hot)
            println(cold)
        end

        nChiller += 1.
    end
end
end

if save
    json_data = json(results);
    open(folderOut * "storage_sens_json_data", "w") do fw
        write(fw, json_data)
    end
end

max_stars_repo_name: jdechalendar/sesi

max_issues_repo_name: jdechalendar/sesi

max_forks_repo_name: jdechalendar/sesi


text:
module ImageBinarization

using Base.Iterators: repeated
using LinearAlgebra
using Polynomials
using Statistics
using Reexport

using ImageContrastAdjustment
@reexport using HistogramThresholding

using ImageCore
using ImageCore.MappedArrays
using ImageCore: GenericGrayImage
using ColorVectorSpace

# TODO: port BinarizationAPI to ImagesAPI
include("BinarizationAPI/BinarizationAPI.jl")
import .BinarizationAPI: AbstractImageBinarizationAlgorithm,
                         binarize, binarize!
import HistogramThresholding: ThresholdAlgorithm

include("integral_image.jl")
include("util.jl")
include("compat.jl")

# Concrete binarization algorithms

# Balanced, Entropy, Intermodes, Minimum, MinimumError, Moments, Otsu, UnimodalRosin, Yen
include("algorithms/single_histogram_threshold.jl")
include("algorithms/adaptive_threshold.jl") # AdaptiveThreshold
include("algorithms/niblack.jl") # Niblack
include("algorithms/polysegment.jl") # Polysegment
include("algorithms/sauvola.jl") # Sauvola

include("deprecations.jl")

export
    # generic API
    binarize, binarize!,

    # Algorithms
    AdaptiveThreshold, recommend_size,
    Niblack,
    Polysegment,
    Sauvola,

    # also reexport algorithms in HistogramThresholding
    SingleHistogramThreshold

end # module ImageBinarization


text:
using BioSequences
using DataStructures: DefaultDict, OrderedDict
using TextWrap: wrap
import JLD: load, save
import Blosc
import JSON
import GZip
### Main calling function:
function run_calling_pipeline(args)
  # get samples/fastq files from command line parameters:
  sample_files = build_sample_files(args["sample_input_file"], args["1"], args["2"])

  # check if DB and fast(Q/A) files exist:
  check_files([args["db"];[fastq for fq_files in values(sample_files) for fastq in fq_files]])
  # array for saving each sample result:
  sample_results = []
  @info("Opening kmer database ... ")
  kmer_db, loci, loci2alleles, db_coverage, k, profile, build_args = open_db(args["db"])

  # prepend base path of the kmer database on to the relative paths for allele fasta files stored in the db
  for (index, file) in enumerate(build_args["fasta_files"])
    build_args["fasta_files"][index] = joinpath(dirname(args["db"]), file)
  end
  # check if scheme fasta files exist
    check_files(build_args["fasta_files"])
  # process each sample:
  for (sample, fq_files) in sample_files
    @info("Sample: $sample. Opening fastq file(s) and counting kmers ... ")
    kmer_count = count_kmers(DNAKmer{k}, fq_files, args["fasta"])
    @info("Voting for alleles ... ")
    votes, loci_votes = count_votes(kmer_count, kmer_db, loci2alleles, db_coverage)
    @info("Calling alleles and novel alleles ...")
    # If fasta are given as input, set kt as 1:
    kt = args["fasta"] ? 1 : args["kt"]
    allele_calls, voting_result = call_alleles(DNAKmer{k}, kmer_count, votes, loci_votes, loci, loci2alleles, build_args["fasta_files"], kt, args["mutation_threshold"], args["output_votes"])
    push!(sample_results, (sample, allele_calls, voting_result))
  end
  @info("Writing output ...")
  write_calls(sample_results, loci, loci2alleles, args["o"], profile, args["output_special"], args["output_votes"])

  @info("Done.")
end

## Helper structs:
struct AlleleCoverage
  allele::Int16
  votes::Int64
  depth::Int16
  covered_kmers::Int16
  uncovered_kmers::Int16
  gap_list::Vector{Tuple{Int16,Int16}}
end

mutable struct NovelAllele
  template_allele::Int
  sequence::String
  n_mutations::Int8
  mutations_list::Array{Any}
  depth::Int16
  uncorrected_gaps::Array{Tuple{Int32,Int32}}
end

function EmptyNovelAllele()
  return NovelAllele(-1,"",0,[],0,[])
end

struct AlleleCall
  allele::String
  flag::String
  coverage::Float16
  depth::Int
  report_txt::String
  novel_allele::NovelAllele
  alleles_to_check::Vector{Tuple{String,String, String}} # tuples of allele label, sequence, report.
end

function open_db(filename)
  # profile:
  profile = nothing
  if isfile("$filename.profile")
    types = Array{String}[]
    open("$filename.profile") do f
      header = split(readline(f),"\t")
      for l in eachline(f)
        values = split(strip(l),"\t")
        push!(types, values)
      end
      profile = Dict("header"=>header, "types"=>types)
    end
  end

  # Compressed database, open and decompress/decode in memory:
  d = load("$filename")
  @info("Finished the JLD load, building alleles list...")
  # alleles_list might be split into smaller parts:
  alleles_list = []
  if haskey(d, "alleles_list")
    alleles_list = Blosc.decompress(Int16, d["alleles_list"])
  else
    idx = 0
    while haskey(d, "alleles_list_$idx")
      append!(alleles_list,Blosc.decompress(Int16, d["alleles_list_$idx"]))
      idx += 1
    end
  end
  @info("Decompressing weight list...")
  build_args = JSON.parse(d["args"])
  k = build_args["k"]
  loci = d["loci"]
  loci_list = Blosc.decompress(Int32, d["loci_list"])
  # weight_list = Blosc.decompress(Int8, d["weight_list"])
  weight_list = Blosc.decompress(Int16, d["weight_list"])
  allele_ids_per_locus = Blosc.decompress(Int, d["allele_ids_per_locus"])
  kmer_str = d["kmer_list"]
  # build a dict to transform allele idx (1,2,...) to original allele ids:
  # TODO: loci2alleles can be a vector, since indexes are 1..length(loci)
  loci2alleles = Dict{Int16, Vector{Int16}}(idx => Int16[] for (idx,locus) in enumerate(loci))
  allele_ids_idx = 1
  locus_idx = 1
  while allele_ids_idx < length(allele_ids_per_locus)
    n_alleles = allele_ids_per_locus[allele_ids_idx]
    append!(loci2alleles[locus_idx], allele_ids_per_locus[allele_ids_idx+1:allele_ids_idx + n_alleles])
    allele_ids_idx += 1 + n_alleles
    locus_idx += 1
  end
  # build the kmer db in the usual format:
  @info("Building kmer index ...")
  kmer_classification = Dict{DNAKmer{k}, Vector{Tuple{Int16, Int16, Vector{Int16}}}}()
  # tuple is locus idx, weight, and list of alleles;
  loci_list_idx = 1
  allele_list_idx = 1
  kmer_idx = 1
  weight_idx = 1
  # for each loci, the number of kmers and the locus idx;
  while loci_list_idx < length(loci_list)
    n_kmers = loci_list[loci_list_idx]
    locus_idx = loci_list[loci_list_idx+1]
    loci_list_idx += 2
    for i in 1:n_kmers
      # get current kmer
      kmer = DNAKmer{k}(kmer_str[kmer_idx:kmer_idx+k-1])
      kmer_idx += k
      # get list of alleles for this kmer:
      n_alleles = alleles_list[allele_list_idx]
      current_allele_list = alleles_list[allele_list_idx+1:allele_list_idx+n_alleles]
      allele_list_idx += n_alleles + 1
      weight = weight_list[weight_idx]
      weight_idx += 1
      # save in db
      if !haskey(kmer_classification, kmer)
        kmer_classification[kmer] = Tuple{Int16, Int16, Vector{Int16}}[]
      end
      push!(kmer_classification[kmer], (locus_idx, weight, current_allele_list))
    end
  end
  # db_coverage: [i][j] -> number of kmers that cover allele j in locus i
  db_coverage = [Int16[] for (idx,locus) in enumerate(loci)]
  current_idx = 1
  if haskey(d, "allele_coverages") # might not exist for older dbs;
    allele_coverages_list = Blosc.decompress(Int16,d["allele_coverages"])
    for (idx,locus) in enumerate(loci)
      n_alleles = length(loci2alleles[idx])
      db_coverage[idx] = allele_coverages_list[current_idx:current_idx+n_alleles-1]
      current_idx += n_alleles
    end
  end
  return kmer_classification, loci, loci2alleles, db_coverage, k, profile, build_args
end

function _find_profile(alleles, profile)
  if profile == nothing
    return 0,0
  end
  l = length(alleles)
    alleles_str = map(string,alleles)
  for genotype in profile["types"]
    # first element is the type id, all the rest is the actual genotype:
    if alleles_str == genotype[2:1+l]
      # assuming Clonal Complex after genotype, when present.
      return genotype[1], length(genotype) >= l + 2 ? genotype[l+2] : ""
    end
  end
  return 0,0
end

function read_alleles(fastafile, ids)
  alleles = Dict{Int, String}()
  idx = 1
  fh = open(fastafile, "r")
  reader = FASTA.Reader(fh)
  record = FASTA.Record()
  while !eof(reader)
    read!(reader, record)
    if in(idx, ids)
      alleles[idx] = FASTA.sequence(record)
    end
    idx += 1
  end
  close(fh)
  return alleles
end

function find_allele(fastafile, wanted_seq)
  idx = 1
  reader = FASTA.Reader(open(fastafile, "r"))
  record = FASTA.Record()
  while !eof(reader)
    read!(reader, record)
    if FASTA.sequence(fastafile) == wanted_seq
      return idx
    end
    idx += 1
  end
  return nothing
end

function LocusNotPresentCall()
  return AlleleCall("0", "", 0.0, 0, "Not present, no kmers found.", EmptyNovelAllele(), [])
end

function LocusPresentCall(call, depth, txt)
  return AlleleCall(call, "", 1.0, depth, txt, EmptyNovelAllele(), [])
end

function MultiplePresentCall(call, depth, txt, alleles_to_check)
  return AlleleCall(call, "+", 1.0, depth, txt, EmptyNovelAllele(), alleles_to_check)
end

function UncoveredAlleleCall(report_txt, alleles_to_check)
  # return AlleleCall("0", "?", 0.0, 0, report_txt, nothing, alleles_to_check)
  return AlleleCall("0", "?", 0.0, 0, report_txt, EmptyNovelAllele(), alleles_to_check)
end

function NovelAlleleCall(report_txt, novel_allele, alleles_to_check)
  return AlleleCall("N", "", 1.0, novel_allele.depth, report_txt, novel_allele, alleles_to_check)
end

function PartiallyCoveredAllele(call, coverage, depth, report_txt, alleles_to_check)
  return AlleleCall(call, "-", coverage, depth, report_txt, EmptyNovelAllele(), alleles_to_check)
end

N_SEARCH = 40 # number of alleles to search, in order of voting;
function call_alleles(::Type{DNAKmer{k}}, kmer_count, votes, loci_votes, loci, loci2alleles, fasta_files, kmer_thr, max_mutations, output_votes) where {k}

  # Function to call for each allele
  function call_allele(votes_per_allele, locus_votes, locus, locus2alleles, fasta_file)
    # given a subset of the voted alleles, find the allele call:
    function call_from_selected_allele_votes(selected_allele_votes)
      allele_set = Set([al for (al, votes) in selected_allele_votes])
      allele_seqs = read_alleles(fasta_file, allele_set)
      allele_coverage = [AlleleCoverage(al, votes, sequence_coverage(DNAKmer{k}, allele_seqs[al], kmer_count, kmer_thr)...) for (al, votes) in selected_allele_votes]
      # filter to find fully covered alleles:
      covered = [x for x in allele_coverage if (x.depth >= kmer_thr)] # if I remove alleles with negative votes, I might the reconstruct the same allele on the novel rebuild; better to flag output
      if length(covered) == 1
        # choose the only fully covered allele
        allele_label = locus2alleles[covered[1].allele]
        allele_votes  = covered[1].votes
        warning = allele_votes < 0 ? " Call with negative votes: $allele_votes." : ""
        return LocusPresentCall("$allele_label", covered[1].depth,  "Called allele $allele_label.$warning")
      end
      if length(covered) > 1
        # choose the most voted, but indicate that there are more fully covered that might be possible:
        allele_label = locus2alleles[covered[1].allele]
        multiple_alleles = join([locus2alleles[x.allele] for x in covered],", ")
        multiple_votes = join([x.votes for x in covered],", ")
        multiple_depth = join([x.depth for x in covered],", ")
        report_txt = "Multiple possible alleles:$multiple_alleles with depth $multiple_depth and votes $multiple_votes. Most voted ($allele_label) is chosen on call file."
        alleles_to_check = [("$(locus2alleles[x.allele])", allele_seqs[x.allele], "Multiple alleles present for this locus. Depth:$(x.depth) Votes:$(x.votes)") for x in covered]
        return MultiplePresentCall("$allele_label", covered[1].depth, report_txt, alleles_to_check)
      end
      # got here, no allele is covered; if the uncovered part is small enough, try novel allele
      # Try to find novel; sort by # of uncovered kmers, smallest to largest
      sorted_allele_coverage = sort(allele_coverage, by=x->x.uncovered_kmers)
      uncovered_kmers = sorted_allele_coverage[1].uncovered_kmers
      covered_kmers = sorted_allele_coverage[1].covered_kmers
      coverage = covered_kmers / (covered_kmers + uncovered_kmers)
      # each mutation gives around k uncovered kmers, so if smallest_uncovered > k * max_mutations, then declare missing;
      if uncovered_kmers > k * max_mutations # consider not present; TODO: better criteria?
        allele = sorted_allele_coverage[1].allele
        allele_label = "$(locus2alleles[allele])"
        depth = sorted_allele_coverage[1].depth
        report_txt = "Not present; allele $allele_label is the best covered but below threshold with $uncovered_kmers/$(covered_kmers+uncovered_kmers) missing kmers."
        return UncoveredAlleleCall(report_txt, [(allele_label, allele_seqs[allele], "Best covered allele, but declared not present, with $uncovered_kmers uncovered kmers and coverage $coverage.")])
      end

      # If gap is small, try to find a novel allele using an existing allele as template;
      # and filling the 'gaps' of coverage searching for mutations.
      # sort candidates by smaller # of gaps, then larger vote:
      sorted_by_gap_allele_coverage = sort(allele_coverage, by=x->(length(x.gap_list),-x.votes))
      best = sorted_by_gap_allele_coverage[1]
      # get all templates with same # of gaps:
      allele_templates = [x for x in sorted_by_gap_allele_coverage if length(x.gap_list) == length(best.gap_list)]
      candidate_list = [] # list with all candidate novel;
      for al_cov in allele_templates # al_cov is a struct AlleleCoverage
        template_seq = allele_seqs[al_cov.allele]
        novel_allele = correct_template(DNAKmer{k}, allele_seqs[al_cov.allele], al_cov.gap_list, kmer_count, kmer_thr, max_mutations)
        novel_allele.template_allele = locus2alleles[al_cov.allele]
        push!(candidate_list, (novel_allele, al_cov))
      end
      # sort by no gaps, then minimum mutations; then more votes
      sorted_candidates = sort(candidate_list,by=x->(length(x[1].uncorrected_gaps),x[1].n_mutations,-x[2].votes))
      # get best
      novel_allele, al_cov = sorted_candidates[1]
      template_allele_label = "$(novel_allele.template_allele)"

      if length(novel_allele.uncorrected_gaps) == 0
        # report:
        mutations_txt = novel_allele.n_mutations > 1 ? "$(novel_allele.n_mutations) mutations" : "$(novel_allele.n_mutations) mutation"
        mutation_desc = join([describe_mutation(ev) for ev in novel_allele.mutations_list], ", ")
        # save, closest and novel:
        alleles_to_check =[(template_allele_label, novel_allele.sequence, "Template for novel allele."),
          ("Novel", novel_allele.sequence, "$mutations_txt from allele $template_allele_label: $mutation_desc")]
        report_txt = "Novel, $mutations_txt from allele $template_allele_label: $mutation_desc"
        return NovelAlleleCall(report_txt, novel_allele, alleles_to_check)

      else # TODO: if all gaps have not been covered, just call directly uncovered? (no novel at all)
        # did not a fully corrected novel; possibly an uncovered allele, maybe partially corrected; output allele/N
        report_txt = "Partially covered allele or novel allele; Best allele $template_allele_label has $uncovered_kmers/$(covered_kmers+uncovered_kmers) missing kmers, and no novel was found. Gaps on positions: $(join(al_cov.gap_list, ','))"
        alleles_to_check =[("$(novel_allele.template_allele)", novel_allele.sequence, "Best allele $template_allele_label, but partially covered.")]
        return PartiallyCoveredAllele("$template_allele_label", round(coverage,digits=4), 0, report_txt, alleles_to_check) # if it's not 100% covered, then depth=0
      end
    end

    # if locus_votes == 0   # 1st - if no locus votes, no presence:
    #   return LocusNotPresentCall()
    # end
    # Votes: Dict{Int16,Dict{Int16,Int64}}; locus_idx -> { allele_idx -> votes};

    sorted_voted_alleles = sort(collect(votes_per_allele), by=x->-x[2]) # sort by max votes for this locus;
    tie_idx = 1
    while tie_idx < length(sorted_voted_alleles) && sorted_voted_alleles[tie_idx+1][2] == sorted_voted_alleles[1][2]
      tie_idx += 1
    end
    # tied_best = [x for x in sorted_voted_alleles if x[2] == sorted_voted_alleles[1][2]]
    tie_idx +=10 # add a bit more
    tied_best = sorted_voted_alleles[1:min(tie_idx,end)]
    allele_call = call_from_selected_allele_votes(tied_best)

    # if it is novel; check if it already exists in DB:
    # if allele_call.allele == "N"
    #   allele_idx = find_allele(fastafile, wanted_seq)
    #   if exists != nothing

    # If I did not get a full covered call, or novel, try looking at ALL alleles
    if (allele_call.depth == 0 || allele_call.allele == "N") && (tie_idx <= length(sorted_voted_alleles))
      second_call = call_from_selected_allele_votes(sorted_voted_alleles[tie_idx:end])
      if second_call.depth == 0
        return allele_call
      else
        return second_call
      end
    else
      return allele_call
    end
  end
  # Function to call all alleles by vote only:
  function call_by_vote(loci, votes, loci_votes, loci2alleles)
    ties = DefaultDict{String,Vector{Int16}}(() -> Int16[])
    vote_log = []
    best_voted_alleles = []
    for (idx,locus) in enumerate(loci)
      # if locus has no votes (no kmer present), output a zero:
      if loci_votes[idx] == 0
        push!(best_voted_alleles, "0")
        push!(vote_log, ("0","No kmers matches on this locus."))
        continue
      end
      # if it has votes, sort and pick the best;
      sorted_vote = sort(collect(votes[idx]), by=x->-x[2])
      push!(best_voted_alleles, loci2alleles[idx][sorted_vote[1][1]])
      votes_txt = join(["$(loci2alleles[idx][al])($votes)" for (al,votes) in sorted_vote[1:min(N_SEARCH,end)]],",")
      push!(vote_log, (loci_votes[idx], votes_txt))
      # find if there was a tie:
      if length(sorted_vote) > 1 && sorted_vote[1][2] == sorted_vote[2][2] # there is a tie;
        for i in 1:length(sorted_vote)
          if sorted_vote[i][2] < sorted_vote[1][2]
            break
          end
          push!(ties[locus], sorted_vote[i][1])
        end
      end
    end
    return best_voted_alleles, vote_log, ties
  end

  # get all calls:
  allele_calls = [call_allele(votes[idx], loci_votes[idx], loci[idx], loci2alleles[idx], fasta_files[idx]) for (idx, locus) in enumerate(loci)]

  # Also do voting output?
  if output_votes
    voting_result = call_by_vote(loci, votes, loci_votes, loci2alleles)
  else
    voting_result = nothing
  end

  return allele_calls, voting_result
end

function write_calls(sample_results, loci, loci2alleles, filename, profile, output_special_cases, output_votes)
# function write_calls(loci2alleles, allele_calls, loci, voting_result, sample, filename, profile, output_special_cases)
  # write the main call:
  open(filename, "w") do f
    header = join(vcat(["Sample"], loci, ["ST", "clonal_complex"]), "\t")
    write(f,  "$header\n")
    for (sample, allele_calls, voting_result) in sample_results
      allele_label_calls = ["$(call.allele)$(call.flag)" for call in allele_calls]
      st, clonal_complex = _find_profile(allele_label_calls, profile)
      calls = join(vcat([sample], allele_label_calls, ["$st", clonal_complex]), "\t")
      write(f, "$calls\n")
    end
  end
  # write the call report:
  open("$filename.coverage.txt", "w") do f
    write(f, "Sample\tLocus\tCoverage\tMinKmerDepth\tCall\n")
    for (sample, allele_calls, voting_result) in sample_results
      for (locus, call) in zip(loci, allele_calls)
        write(f, "$sample\t$locus\t$(call.coverage)\t$(call.depth)\t$(call.report_txt)\n")
      end
    end
  end
  # write the alleles with novel, missing, or multiple calls:
  alleles_to_check = [(locus, allele_to_check...) for (sample, allele_calls, voting_result) in sample_results for (locus, call) in zip(loci, allele_calls) for allele_to_check in call.alleles_to_check if length(call.alleles_to_check) > 0]
  if output_special_cases && length(alleles_to_check) > 0
    open("$filename.special_cases.fa", "w") do f
      # TODO: deal with repeated alleles;
      for (locus, al, seq, desc) in alleles_to_check
        write(f,">$(locus)_$al $desc\n$seq\n")
      end
    end
  end
  # write novel alleles
  # seen_novel = DefaultDict{String, Set{String}}(() -> Set{String}()) # locus to set of seen novel sequences
  times_seen = DefaultDict{String, Int}(0) # times that a novel DNA sequence was seen; (DNA -> count)
  novel_id = Dict{String, Int}() # novel id of a DNA seq (DNA -> id)
  current_novel_id = DefaultDict{String, Int}(0) # current novel id for a given locus. (locus -> current_id)
  novel_to_fasta = []
  open("$filename.novel.txt", "w") do text
    write(text, "Sample\tLocus\tNovel_id\tMinKmerDepth\tNmut\tDesc\n")
    for (sample, allele_calls, voting_result) in sample_results
      # println(novel_to_fasta)
      # println("------------")
      novel_alleles = [(locus, call.novel_allele) for (locus, call) in zip(loci, allele_calls) if call.novel_allele.template_allele != -1]
      # print("novel_allele: ")
      # println(novel_alleles)
      for (locus, novel_allele) in novel_alleles
        mutation_desc = join([describe_mutation(ev) for ev in novel_allele.mutations_list], ", ")
        # check if this novel was seen before:
        if (times_seen[novel_allele.sequence] == 0)
          # println("new novel!")
          current_novel_id[locus] += 1
          # println(locus)
          # println(current_novel_id[locus])
          push!(novel_to_fasta, (locus, novel_allele.sequence))
          novel_id[novel_allele.sequence] = current_novel_id[locus]
          write(text, "$sample\t$locus\tN$(novel_id[novel_allele.sequence])\t$(novel_allele.depth)\t$(novel_allele.n_mutations)\tFrom allele $(novel_allele.template_allele), $mutation_desc.\n")
        end
        times_seen[novel_allele.sequence] += 1
        # println(">>> Times we seen this allele")
        # println(times_seen[novel_allele.sequence])
        # println("--------")
        # println(novel_allele.sequence)
        # println(novel_id[novel_allele.sequence])

      end
    end
  end
  # now write the fasta:
  # println("<><><><><><><><><><><><><><><><><><>")
  # println(novel_to_fasta)
  open("$filename.novel.fa", "w") do f
    for (locus, seq) in novel_to_fasta
      desc = "Seen in $(times_seen[seq]) sample(s)"
      write(f, ">$(locus)_N$(novel_id[seq]) $desc.\n$(wrap(seq, width=120))\n")
    end
  end

  # TODO: with multiple samples, we have to check how many
  # novel_alleles = [(locus, call.novel_allele) for (locus, call) in zip(loci, allele_calls) if call.novel_allele.template_allele != -1]
  # if length(novel_alleles) > 0
  #   open("$filename.novel.fa", "w") do fasta
  #       for (locus, novel_allele) in novel_alleles
  #         write(fasta, ">$locus\n$(novel_allele.sequence)\n")
  #         mutation_desc = join([describe_mutation(ev) for ev in novel_allele.mutations_list], ", ")
  #         write(text, "$locus\t$(novel_allele.depth)\t$(novel_allele.n_mutations)\tFrom allele $(novel_allele.template_allele), $mutation_desc.\n")
  #       end
  #     end
  #   end
  # end
  # write also the votes if we got them:
  if output_votes
    f_vote_call = open("$filename.byvote", "w")
    f_details = open("$filename.votes.txt", "w")
    f_ties = open("$filename.ties.txt", "w")
    # headers:
    write(f_vote_call,  "$(join(vcat(["Sample"], loci, ["ST", "clonal_complex"]), "\t"))\n")
    write(f_details, "Sample\tLocus\tTotal locus votes\tAllele(relative votes),...\n")
    write(f_ties, "Sample\tLocus\tTied Alleles\n")
    # loop per sample:
    for (sample, allele_calls, voting_result) in sample_results
      # write the voting call:
      best_voted_alleles, vote_log, ties = voting_result
      st, clonal_complex = _find_profile(best_voted_alleles, profile)
      calls = join(vcat([sample], best_voted_alleles, ["$st", clonal_complex]), "\t")
      write(f_vote_call, "$calls\n")
      # write the detailed votes:
      for (locus, data) in zip(loci, vote_log)
        write(f_details, "$sample\t$locus\t$(join(data,'\t'))\n")
      end
      # write ties:
      if length(ties) > 0
        for (locus, tied_alleles) in sort(collect(ties), by=x->x[1])
          ties_txt = join(["$t" for t in sort(tied_alleles)],", ")
          write(f_ties, "$sample\t$locus\t$ties_txt\n")
        end
      end
    end
    close(f_vote_call)
    close(f_details)
    close(f_ties)
  end
end

function count_kmers_in_fasta_file(::Type{DNAKmer{k}}, fasta_file, kmer_count) where {k}
  reader = FASTA.Reader(open(fasta_file, "r"))
	record = FASTA.Record()
  while !eof(reader)
      read!(reader, record)
      for (pos, kmer) in each(DNAKmer{k}, FASTA.sequence(record))
        kmer_count[canonical(kmer)] += 1
      end
    end
end

function count_kmers_in_fastq_file(::Type{DNAKmer{k}}, fastq_file, kmer_count) where {k}
  reader = FASTQ.Reader(endswith(fastq_file, ".gz") ? GZip.open(fastq_file, "r") : open(fastq_file, "r"))
  record = FASTQ.Record()
  while !eof(reader)
      read!(reader, record)
      for (pos, kmer) in each(DNAKmer{k}, FASTQ.sequence(record))
        kmer_count[canonical(kmer)] += 1
      end
    end
end

function count_kmers(::Type{DNAKmer{k}}, files, is_fasta=false) where {k}
  # Count kmers:
  kmer_count = DefaultDict{DNAKmer{k},Int}(0)
  for f in files
    # check type:
    if is_fasta
      count_kmers_in_fasta_file(DNAKmer{k}, f, kmer_count)
    else
      count_kmers_in_fastq_file(DNAKmer{k}, f, kmer_count)
    end
  end
  return kmer_count
end

function count_votes(kmer_count, kmer_db, loci2alleles, db_coverage)
# function count_kmers_and_vote{k}(::Type{DNAKmer{k}}, files, kmer_db, loci2alleles)
  # Now count votes:
  # 0 votes for all alleles everyone at the start:
  # TODO: also list, instead of DICT?
  votes = Dict(locus_idx => Dict{Int16, Int}(i => 0 for i in 1:length(alleles)) for (locus_idx,alleles) in loci2alleles)
  # votes per locus, to decide presence/absence:
  loci_votes = DefaultDict{Int16, Int}(0)
  # for (locus,) in loci2alleles
  #     loci_votes[locus] = 0
  # end
  for (kmer, count) in kmer_count
    if haskey(kmer_db, kmer)
      for (locus, weight, alleles) in kmer_db[kmer]
        v = weight * count
        loci_votes[locus] += abs(v)
        for allele in alleles
          votes[locus][allele] += v
        end
      end
    end
  end
  # adjust for allele coverage:
  # for (locus, allele_votes) in votes
  #   for al in keys(allele_votes)
  #     allele_votes[al] = Int(round(allele_votes[al]/db_coverage[locus][al]))
  #   end
  # end
  return votes, loci_votes
end

# Calling:
function sequence_coverage(::Type{DNAKmer{k}}, sequence, kmer_count, kmer_thr=6) where {k}
  # TODO: not using the gap list anymore, just for debug. Eventually delete!
  # find the minimum coverage depth of all kmers of the given variant, and
  # how many kmers are covered/uncovered
  gap_list = Tuple{Int,Int}[]
  smallest_coverage::Int16 = 32766
  covered = 0
  uncovered = 0

  in_gap = false
  gap_start = 0

  for (pos, kmer) in each(DNAKmer{k}, DNASequence(sequence))
    can_kmer = canonical(kmer)
    cov = get(kmer_count, can_kmer, 0)
    if cov >= kmer_thr
      covered += 1
      # gap:
      if in_gap
        push!(gap_list, (gap_start, pos-1))
        in_gap = false
      end
    else
      uncovered += 1
      # gap
      if !in_gap
        gap_start = pos
        in_gap = true
      end
    end
    if cov < smallest_coverage
      smallest_coverage = cov
    end
  end
  if in_gap
    push!(gap_list, (gap_start, length(sequence)-k+1))
  end
  # merge too close gaps: if by chance we get a kmer match in the middle of a gap, this will break a gap in two, even
  # with just one mutation, giving errors. TODO: How to prevent that?
  # simple idea here, merge gaps (s1,e1) and (s2,e2) if s1 + k >= s2;
  merged_gap_list = []
  i = 1
  while i <= length(gap_list)
    # if too close to the next, merge:
    if i<length(gap_list) && gap_list[i][1] + k >= gap_list[i+1][1]
      push!(merged_gap_list, (gap_list[i][1], gap_list[i+1][2]))
      i += 2
    else
      push!(merged_gap_list, gap_list[i])
      i += 1
    end
  end
  return smallest_coverage, covered, uncovered, merged_gap_list
end

function cover_sequence_gap(::Type{DNAKmer{k}}, sequence, kmer_count, kmer_thr=8, max_mutations=4) where {k}
  gap_cover_list = []
  found = Set{String}()
  function scan_variant(n_mut, sequence, events, start_pos=1)
    if n_mut > max_mutations
      return
    end
    last_kmer_present = false
    for (pos, kmer) in each(DNAKmer{k}, DNASequence(sequence))
      if pos < start_pos
        continue
      end
      can_kmer = canonical(kmer)
      if get(kmer_count, can_kmer, 0) >= kmer_thr
        if pos > start_pos +1 && (!last_kmer_present)
          # found a kmer and the last was not there; check if I can find kmer variations for the previous position
          km = "$kmer"
          for base in ["A","C","T","G"]
            if get(kmer_count, canonical(DNAKmer{k}("$base$(km[1:end-1])")), 0) >= kmer_thr
              # found a kmer variant;
              # Test one more kmer before adding to the pool:
              # substitution at position pos-1
              substitution_seq = sequence[1:pos-2] * base * sequence[pos:end]
              push!(candidate_list, (n_mut+1, substitution_seq, [events; [("S",pos-1,"$(sequence[pos-1])->$base")]], start_pos))
              # insertion:
              insertion_seq = sequence[1:pos-1] * base * sequence[pos:end]
              push!(candidate_list, (n_mut+1, insertion_seq, [events; [("I",pos-1,base)]], start_pos))
              # deletion:
              for i in 0:2
                if pos-2-i < 1
                  break
                end
                if sequence[pos-2-i:pos-2-i] == base
                  deletion_seq = sequence[1:pos-2-i] * sequence[pos:end]
                  push!(candidate_list, (n_mut+1+i, deletion_seq, [events; [("D",pos-1,"$(i+1)")]], start_pos))
                  break
                end
              end
              # this variant has mismatches, so just return;
              return
            end
          end
        end
        last_kmer_present = true
      else
        if last_kmer_present
          # last_kmer_present was present, but this not; check kmer variations for nest position
          km = "$kmer"
          for base in ["A","C","T","G"]
            if get(kmer_count, canonical(DNAKmer{k}("$(km[1:end-1])$base")), 0) >= kmer_thr
              # println("Found a $base at pos $(pos+k-1)")
              # substitution at position pos+k-1;
              substitution_seq = sequence[1:pos+k-2] * base * sequence[pos+k:end]
              push!(candidate_list, (n_mut+1, substitution_seq, [events; [("S",pos+k-1,"$(sequence[pos+k-1])->$base")]], pos-1))
              # insertion:
              insertion_seq = sequence[1:pos+k-2] * base * sequence[pos+k-1:end]
              push!(candidate_list, (n_mut+1, insertion_seq, [events; [("I",pos+k-1,base)]], pos-1))
              # deletion:
              for i in 0:2
                # println("trying $(pos+k+i), found $(sequence[pos+k+i]) (want $base) $(sequence[pos+k+i:pos+k+i] == base).")
                if pos+k+i > length(sequence)
                  break
                end
                if sequence[pos+k+i:pos+k+i] == base
                  # println("Found a possible deletion, pos $(pos+k-1) to $(pos+k+i-1)")
                  deletion_seq = sequence[1:pos+k-2] * sequence[pos+k+i:end]
                  push!(candidate_list, (n_mut+1+i, deletion_seq, [events; [("D",pos+k-1,"$(i+1)")]], pos-1))
                  break
                end
              end
              # the original sequence has mismatches, so just return;
              return
            end
          end
        end
        last_kmer_present = false
      end
    end
    # passed without suggesting mutations, test this variant:
    seq_cov = sequence_coverage(DNAKmer{k}, sequence, kmer_count)[1]
    if seq_cov[1] >= kmer_thr
      if !in(sequence, found)
        push!(gap_cover_list, (n_mut, sequence, events, seq_cov))
        push!(found, sequence)
        max_mutations = n_mut
      end
    end
  end
  candidate_list = [(0, sequence, [], 0)]
  idx = 1
  while idx <= length(candidate_list)
    n_mut, candidate_sequence, mutations, start_pos = candidate_list[idx]
    scan_variant(n_mut, candidate_sequence, mutations, start_pos)
    idx += 1
  end

  if length(gap_cover_list) > 0
    # sort and return the one with minimum mutations: TODO: tie?
    return sort(gap_cover_list, by=x->x[1])[1]
  else
    return nothing
  end

end
# struct CorrectedTemplate
#
# end
function correct_template(::Type{DNAKmer{k}}, template_seq, gap_list, kmer_count, kmer_thr, max_mutations) where {k}
  function find_next_gap(sequence, skip=1)
    in_gap = false
    gap_start = 0
    for (pos, kmer) in each(DNAKmer{k}, DNASequence(sequence))
      if pos < skip
        continue
      end
      can_kmer = canonical(kmer)
      cov = get(kmer_count, can_kmer, 0)
      if cov >= kmer_thr
        if in_gap
          return (gap_start, pos-1)
        end
      else
        if !in_gap
          gap_start = pos
          in_gap = true
        end
      end
    end # loop
    if in_gap
      return (gap_start, length(sequence)-k+1)
    end
    return nothing
  end # end auxiliar function find_next_gap

  # start correction:
  current_skip = 1
  uncorrected_gaps = []
  corrected_seq = template_seq
  total_mut = 0
  max_depth = 0
  mutations_list = []
  while current_skip < length(corrected_seq)
    gap = find_next_gap(corrected_seq, current_skip)
    if gap == nothing
      break
    end
    st_pos, end_pos = gap
    adj_start = max(st_pos-1, 1)
    adj_end = min(end_pos+k, length(corrected_seq))
    gap_seq = corrected_seq[adj_start:adj_end]
    gap_cover = cover_sequence_gap(DNAKmer{k}, gap_seq, kmer_count, kmer_thr, max_mutations)
    if gap_cover == nothing
      push!(uncorrected_gaps, (st_pos, end_pos))
      current_skip = end_pos + 1
    else
      n_mut, gap_cover_seq, mut_list, depth = gap_cover
      total_mut += n_mut
      if depth > max_depth
        max_depth = depth
      end
      for mut in mut_list
        # TODO: correct idx
        push!(mutations_list, (mut[1], mut[2] + adj_start - 1, mut[3]))
      end
      corrected_seq = corrected_seq[1:adj_start-1] * gap_cover_seq * corrected_seq[adj_end+1:end]
      current_skip = adj_start + length(gap_cover_seq) - k
    end
  end
  # 1st field, template_allele, is filled later:
  return NovelAllele(-1, corrected_seq, total_mut, mutations_list, max_depth, uncorrected_gaps)
end

function describe_mutation(mut)
  mut, pos, desc = mut
  if mut == "D"
    return "Del of len $desc at pos $pos"
  elseif mut == "I"
    return "Ins of base $desc at pos $pos"
  elseif mut == "S"
    return "Subst $desc at pos $pos"
  else
    return ""
  end
end

### input file helper calling_functions

function build_sample_files(input_file, forward_files, reverse_files)
  # First, sanity checks
  # at least one input option:
  if all([x == nothing for x in [input_file, forward_files, reverse_files]])
    exit_error("Please give at least one input file using the parameters -m, --1 and/or --2.")
  end
  #
  if length(reverse_files)>0 && length(forward_files) != length(reverse_files)
    exit_error("Forward and reverse input file does not match, got $(length(forward_files)) forward and $(length(reverse_files)) reverse.")
  end

  # Create an ordered dict to get all samples. Combinations of -m and -1,-2 are allowed.
  sample_files = OrderedDict{String, Vector{String}}()
  # Check
  if input_file != nothing
      check_files([input_file])
      open(input_file) do f
        for ln in eachline(f)
          sample, fastq = split(strip(ln))
          if !haskey(sample_files, sample)
            sample_files[sample] = String[]
          end
          push!(sample_files[sample], fastq)
        end
      end
  end
  if length(reverse_files) == 0 # single files only
    for fw_file in forward_files
      sample_files[remove_fastq_ext(fw_file)] = [fw_file]
    end
  end
  # build sample files for fw and rev:
  for (fw, rev) in zip(forward_files, reverse_files)
    sample = lcp([fw,rev]) # largest common prefix, removing special symbols at the end, will be the name of the sample.
    if sample == "" # if empty, fastq names dont match, probably wrong pairing.
      exit_error("No match between forward and reverse files $fw and $bw, please check the input options --1 and --2.")
    end
    sample_files[sample] = [fw,rev]
  end
  return sample_files
end

# longest common prefix
function lcp(str::Vector{String})
  r = IOBuffer()
  i = 1
  while all(i <= length(s) for s in str) && all(s == str[1][i] for s in getindex.(str, i))
    print(r, str[1][i])
    i += 1
  end
  return strip(basename(String(take!(r))),['_','.','-'])
end

# # Remove fastq
function remove_fastq_ext(str)
  return replace(str,  r"\.(fastq|fq)(\.gz)?" => "")
end

max_stars_repo_name: liberuser/MentaLiST

max_issues_repo_name: liberuser/MentaLiST

max_forks_repo_name: liberuser/MentaLiST


text:
#using Pkg
#Pkg.activate(ENV["HOME"] * "/.julia/environments/" * "covid19")
if Sys.iswindows() == false
    covid19dir = ENV["HOME"] * "/GitHub/covid19"
else
    covid19dir = ENV["HOMEDRIVE"] * ENV["HOMEPATH"] * "\\GitHub\\covid19"
end

if (covid19dir in LOAD_PATH) == false
    push!(LOAD_PATH, pwd());
end

#Pkg.update()
#ENV["PYTHON"] = ENV["CONDA_PYTHON_EXE"]
#]add CSV HTTP GitHub DataFrames PyCall PyPlot Plots Dates Statistics TimeSeries MLJ XGBoost Measurements LsqFit NCDatasets Glob JSON3 Tables Pandas FileIO CSVFiles Missings JLD2 OnlineStats CUDAnative CuArrays
#Pkg.build("PyCall")

using Dates, Statistics, Shapefile, DataFrames, Glob

## Load CENSUS TIGER/Line Shapefile Data
function load_county_data(countyshppath = "./tl_2019_us_county.shp")
    #countyshpdownloadurl = "https://liviandonglai.box.com/s/pv1zdejyy0ntsjikeiljqlfmoxiz8gba"
    #countyshpboxlocation = "/Users/gong/Box/Data/CENSUS/tl_2019_us_county/tl_2019_us_county.shp"
    return countytable = Shapefile.Table(countyshppath);
end

## Load COVID-19 Data
#import COVID19

(covid19g, covid19us, covid19nyt_state, covid19nyt_county) = COVID19.load_data("./");

country = [covid19us[i].country for i in 1:length(covid19us)];
fips = [covid19us[i].fips for i in 1:length(covid19us)];
county = [covid19us[i].county for i in 1:length(covid19us)];
state = [covid19us[i].province_state for i in 1:length(covid19us)];
clat = [covid19us[i].lat for i in 1:length(covid19us)];
clon = [covid19us[i].lon for i in 1:length(covid19us)];
cpop = [covid19us[i].population for i in 1:length(covid19us)];

cconfirmed = [covid19us[i].confirmed[end] for i in 1:length(covid19us)];
dcconfirmed = [mean(covid19us[i].confirmed[end-7:end] - covid19us[i].confirmed[end-8:end-1]) for i in 1:length(covid19us)];
dcconfirmedinst = [mean(covid19us[i].confirmed[end] - covid19us[i].confirmed[end-1]) for i in 1:length(covid19us)];
cdeath = [covid19us[i].death[end] for i in 1:length(covid19us)];
dcdeath = [mean(covid19us[i].death[end-7:end] - covid19us[i].death[end-8:end-1]) for i in 1:length(covid19us)];
dcdeathinst = [mean(covid19us[i].death[end] - covid19us[i].death[end-1]) for i in 1:length(covid19us)];

ustate = unique(state);
mfips = deepcopy(fips);
strfips = string.(fips);
terrind = findall(0 .<= mfips .< 100);
mfips[terrind] = mfips[terrind] .* 1000;
gind = findall(mfips .>= 0);
mfips[gind] = mfips[gind] .+ 100000;
strfips = [string(mfips[i])[2:end] for i in 1:length(mfips)];
statefips = [strfips[i][1:2] for i in 1:length(strfips)];

countytable = load_county_data();
ctcensusfips = countytable.STATEFP .* countytable.COUNTYFP;
ctaland = countytable.ALAND;
ctawater = countytable.AWATER;
ctname = countytable.NAME;

countyarea = Array{Float64}(undef,(length(strfips)));
countyarea .= NaN;
for i = 1:length(strfips)
    if parse(Int,strfips[i][1:2]) <= 59
        local fipsind = findall(strfips[i] .== ctcensusfips);
        if isempty(fipsind) == false
            countyarea[i] = ctaland[fipsind[1]];
        end
    end
end

## Calculating Cases per population/area
#strfips[1:end-110] .* state[1:end-110]
#sort(censusfips)[1:end-90]


text:
using GLFW
using ModernGL
using CSyntax
using CSyntax.CStatic

using CImGui
using CImGui.LibCImGui
using CImGui.GLFWBackend
using CImGui.OpenGLBackend

using ImPlot
import CImGui.LibCImGui: ImGuiCond_Always

@static if Sys.isapple()
    # OpenGL 3.2 + GLSL 150
    tmp_glsl_version = 150
    GLFW.WindowHint(GLFW.CONTEXT_VERSION_MAJOR, 3)
    GLFW.WindowHint(GLFW.CONTEXT_VERSION_MINOR, 2)
    GLFW.WindowHint(GLFW.OPENGL_PROFILE, GLFW.OPENGL_CORE_PROFILE) # 3.2+ only
    GLFW.WindowHint(GLFW.OPENGL_FORWARD_COMPAT, GL_TRUE) # required on Mac
else
    # OpenGL 3.0 + GLSL 130
    tmp_glsl_version = 130
    GLFW.WindowHint(GLFW.CONTEXT_VERSION_MAJOR, 3)
    GLFW.WindowHint(GLFW.CONTEXT_VERSION_MINOR, 0)
    # GLFW.WindowHint(GLFW.OPENGL_PROFILE, GLFW.OPENGL_CORE_PROFILE) # 3.2+ only
    # GLFW.WindowHint(GLFW.OPENGL_FORWARD_COMPAT, GL_TRUE) # 3.0+ only
    GLFW.WindowHint(GLFW.FOCUSED,GL_TRUE)
end
const glsl_version=tmp_glsl_version

# setup GLFW error callback
error_callback(err::GLFW.GLFWError) = @error "GLFW ERROR: code $(err.code) msg: $(err.description)"
GLFW.SetErrorCallback(error_callback)

# create window
window = GLFW.CreateWindow(1100, 720, "PlayGround")

@assert window != C_NULL
GLFW.MakeContextCurrent(window)
GLFW.SwapInterval(1)  # enable vsync

# setup Dear ImGui context
ctx = CImGui.CreateContext()
ctxp = ImPlot.CreateContext()

# setup Dear ImGui style
CImGui.StyleColorsClassic()

# setup Platform/Renderer bindings
ImGui_ImplGlfw_InitForOpenGL(window, true)
ImGui_ImplOpenGL3_Init(glsl_version)

try
    clear_color = Cfloat[0.45, 0.55, 0.60, 1.00]

    xmin=0.0
    xmax=1000.0

    while  !GLFW.WindowShouldClose(window)
        GLFW.PollEvents()

        # start the Dear ImGui frame
        ImGui_ImplOpenGL3_NewFrame()
        ImGui_ImplGlfw_NewFrame()
        CImGui.NewFrame()
        CImGui.SetNextWindowSize(CImGui.ImVec2(350.0,200.0),CImGui.ImGuiCond_Once);

        @c CImGui.Begin("Plot Window")
        y = rand(1000)
        ImPlot.SetNextPlotLimits(xmin,xmax,0.0,1.0, ImGuiCond_Once)
        if (ImPlot.BeginPlot("Foo", "x1", "y1", CImGui.ImVec2(-1,300)))
            ImPlot.PlotLine(y)
            ImPlot.EndPlot()
        end
        CImGui.End()

        # rendering
        CImGui.Render()
        GLFW.MakeContextCurrent(window)
        display_w, display_h = GLFW.GetFramebufferSize(window)
        glViewport(0, 0, display_w, display_h)
        glClearColor(clear_color...)
        glClear(GL_COLOR_BUFFER_BIT)
        ImGui_ImplOpenGL3_RenderDrawData(CImGui.GetDrawData())

        GLFW.MakeContextCurrent(window)
        GLFW.SwapBuffers(window)
    end

catch e
    @error "Error in renderloop!" exception=e
    Base.show_backtrace(stderr, catch_backtrace())
finally
    ImGui_ImplOpenGL3_Shutdown()
    ImGui_ImplGlfw_Shutdown()
    CImGui.DestroyContext(ctx)
    ImPlot.DestroyContext(ctxp)
    GLFW.DestroyWindow(window)
end