text:
# adapted from https://github.com/JuliaLang/julia/blob/cced577b79316ea38482d84b69e7be9666d14041/stdlib/REPL/src/REPLCompletions.jl
# by aviatesk, 2020/02/28

module FuzzyCompletions

export completions, shell_completions, bslash_completions, completion_text, score

using Base.Meta
using Base: propertynames, something
using REPL

# subtracting a small term proportional to levenshtein distance allows for case-sensitive matching
# without affecting other behaviours
const DISCOUNT_COEF_LEVENSTEIN = 1e-4

fuzzyscore(needle::String, haystack::String) =
    REPL.fuzzyscore(needle, haystack) - DISCOUNT_COEF_LEVENSTEIN * REPL.levenshtein(needle, haystack)
fuzzyscore(needle, haystack) = fuzzyscore(string(needle), string(haystack))

abstract type Completion end

struct KeywordCompletion <: Completion
    keyword::String
    score::Float64
end
KeywordCompletion(keyword::String) = KeywordCompletion(keyword, Inf)
KeywordCompletion(keyword::String, needle) = KeywordCompletion(keyword, fuzzyscore(needle, keyword))

struct PathCompletion <: Completion
    path::String
    score::Float64
end
PathCompletion(path::String) = PathCompletion(path, Inf)
PathCompletion(path::String, needle) = PathCompletion(path, fuzzyscore(needle, path))

struct ModuleCompletion <: Completion
    parent::Module
    mod::String
    score::Float64
end
ModuleCompletion(parent::Module, mod::String) = ModuleCompletion(parent, mod, Inf)
ModuleCompletion(parent::Module, mod::String, needle) = ModuleCompletion(parent, mod, fuzzyscore(needle, mod))

struct PackageCompletion <: Completion
    package::String
    score::Float64
end
PackageCompletion(package::String) = PackageCompletion(package, Inf)
PackageCompletion(package::String, needle) = PackageCompletion(package, fuzzyscore(needle, package))

struct PropertyCompletion <: Completion
    value
    property::Symbol
    score::Float64
end
PropertyCompletion(value, property::Symbol) = PropertyCompletion(value, property, Inf)
PropertyCompletion(value, property::Symbol, needle) = PropertyCompletion(value, property, fuzzyscore(needle, property))

struct FieldCompletion <: Completion
    typ::DataType
    field::Symbol
    score::Float64
end
FieldCompletion(typ::DataType, field::Symbol) = FieldCompletion(typ, field, Inf)
FieldCompletion(typ::DataType, field::Symbol, needle) = FieldCompletion(typ, field, fuzzyscore(needle, field))

# NOTE: no fuzzyness is needed for this
struct MethodCompletion <: Completion
    func
    input_types::Type
    method::Method
    orig_method::Union{Nothing,Method}
end

struct BslashCompletion <: Completion
    bslash::String
    score::Float64
end
BslashCompletion(bslash::String) = BslashCompletion(bslash, Inf)
BslashCompletion(bslash::String, needle) = BslashCompletion(bslash, fuzzyscore(needle, bslash))

struct ShellCompletion <: Completion
    text::String
end

struct DictCompletion <: Completion
    dict::AbstractDict
    key::String
    score::Float64
end
DictCompletion(dict::AbstractDict, key::String) = DictCompletion(dict, key, Inf)
DictCompletion(dict::AbstractDict, key::String, needle) = DictCompletion(dict, key, fuzzyscore(needle, key))

# interface definition
function Base.getproperty(c::Completion, name::Symbol)
    if name === :keyword
        return getfield(c, :keyword)::String
    elseif name === :path
        return getfield(c, :path)::String
    elseif name === :parent
        return getfield(c, :parent)::Module
    elseif name === :mod
        return getfield(c, :mod)::String
    elseif name === :package
        return getfield(c, :package)::String
    elseif name === :property
        return getfield(c, :property)::Symbol
    elseif name === :field
        return getfield(c, :field)::Symbol
    elseif name === :method
        return getfield(c, :method)::Method
    elseif name === :bslash
        return getfield(c, :bslash)::String
    elseif name === :text
        return getfield(c, :text)::String
    elseif name === :key
        return getfield(c, :key)::String
    end
    return getfield(c, name)
end

_completion_text(c::KeywordCompletion) = c.keyword
_completion_text(c::PathCompletion) = c.path
_completion_text(c::ModuleCompletion) = c.mod
_completion_text(c::PackageCompletion) = c.package
_completion_text(c::PropertyCompletion) = string(c.property)
_completion_text(c::FieldCompletion) = string(c.field)
_completion_text(c::MethodCompletion) = sprint(io -> show(io, isnothing(c.orig_method) ? c.method : c.orig_method::Method))
_completion_text(c::BslashCompletion) = c.bslash
_completion_text(c::ShellCompletion) = c.text
_completion_text(c::DictCompletion) = c.key

completion_text(c) = _completion_text(c)::String

const Completions = Tuple{Vector{Completion}, UnitRange{Int}, Bool}

score(c::Completion) = c.score
score(c::MethodCompletion) = 0.0
score(c::ShellCompletion) = 0.0

import REPL.REPLCompletions:
    appendmacro!

function filtered_mod_names(ffunc::Function, mod::Module, name::AbstractString, all::Bool = false, imported::Bool = false)
    ssyms = names(mod, all = all, imported = imported)
    filter!(ffunc, ssyms)
    syms = String[string(s) for s in ssyms]
    macros = filter(x -> startswith(x, "@" * name), syms)
    appendmacro!(syms, macros, "_str", "\"")
    appendmacro!(syms, macros, "_cmd", "`")
    filter!(sym -> '#' ∉ sym, syms)
    return [ModuleCompletion(mod, sym, name) for sym in syms]
end

# REPL Symbol Completions
function complete_symbol(sym, ffunc, context_module=Main)::Vector{Completion}
    mod = context_module
    name = sym

    lookup_module = true
    t = Union{}
    val = nothing
    if something(findlast(in(non_identifier_chars), sym), 0) < something(findlast(isequal('.'), sym), 0)
        # Find module
        lookup_name, name = rsplit(sym, ".", limit=2)

        ex = Meta.parse(lookup_name, raise=false, depwarn=false)

        b, found = get_value(ex, context_module)
        if found
            val = b
            if isa(b, Module)
                mod = b
                lookup_module = true
            elseif Base.isstructtype(typeof(b))
                lookup_module = false
                t = typeof(b)
            end
        else # If the value is not found using get_value, the expression contain an advanced expression
            lookup_module = false
            t, found = get_type(ex, context_module)
        end
        found || return Completion[]
        # Ensure REPLCompletion do not crash when asked to complete a tuple, #15329
        !lookup_module && t <: Tuple && return Completion[]
    end

    suggestions = Completion[]
    if lookup_module
        # We will exclude the results that the user does not want, as well
        # as excluding Main.Main.Main, etc., because that's most likely not what
        # the user wants
        p = s->(!Base.isdeprecated(mod, s) && s != nameof(mod) && ffunc(mod, s)::Bool)
        # Looking for a binding in a module
        if mod == context_module
            # Also look in modules we got through `using`
            mods = ccall(:jl_module_usings, Any, (Any,), context_module)::Vector
            for m in mods
                append!(suggestions, filtered_mod_names(p, m::Module, name))
            end
            append!(suggestions, filtered_mod_names(p, mod, name, true, true))
        else
            append!(suggestions, filtered_mod_names(p, mod, name, true, false))
        end
    elseif val !== nothing # looking for a property of an instance
        for property in propertynames(val, false)
            s = string(property)
            push!(suggestions, PropertyCompletion(val, property, name))
        end
    else
        # Looking for a member of a type
        if t isa DataType && t != Any
            # Check for cases like Type{typeof(+)}
            if t isa DataType && t.name === Base._TYPE_NAME
                t = typeof(t.parameters[1])
            end
            # Only look for fields if this is a concrete type
            if isconcretetype(t)
                fields = fieldnames(t)
                for field in fields
                    s = string(field)
                    push!(suggestions, FieldCompletion(t, field, name))
                end
            end
        end
    end
    suggestions
end

import REPL.REPLCompletions:
    sorted_keywords

# NOTE:
# I would like to be a bit strict on `KeywordCompletion`s:
# they are so common that it would look verbose if they appear in every completion.
# I want to restict their fuzzyness only after a strict match on the first character.
function complete_keyword(s::Union{String,SubString{String}})
    c = first(s, 1)
    filtered_keywords = filter(k -> startswith(k, c), sorted_keywords)
    Completion[KeywordCompletion(kw, s) for kw in filtered_keywords]
end

function complete_path(path::AbstractString, pos; use_envpath=false, shell_escape=false)
    # https://github.com/JunoLab/FuzzyCompletions.jl/issues/7
    if endswith(path, ' ')
        # in order to keep type stability, we just return empty completion for this case
        return Completion[], 0:-1, false
        # return REPL.REPLCompletions.complete_path(
        #     path, pos;
        #     use_envpath  = use_envpath,
        #     shell_escape = shell_escape)
    end

    if Base.Sys.isunix() && occursin(r"^~(?:/|$)", path)
        # if the path is just "~", don't consider the expanded username as a prefix
        if path == "~"
            dir, prefix = homedir(), ""
        else
            dir, prefix = splitdir(homedir() * path[2:end])
        end
    else
        dir, prefix = splitdir(path)
    end
    local files
    try
        if isempty(dir)
            files = readdir()
        elseif isdir(dir)
            files = readdir(dir)
        else
            return Completion[], 0:-1, false
        end
    catch
        return Completion[], 0:-1, false
    end

    matches = Set{String}()
    for file in files
        id = try isdir(joinpath(dir, file)) catch; false end
        # joinpath is not used because windows needs to complete with double-backslash
        push!(matches, id ? file * (@static Sys.iswindows() ? "\\\\" : "/") : file)
    end

    if use_envpath && length(dir) == 0
        # Look for files in PATH as well
        local pathdirs = split(ENV["PATH"], @static Sys.iswindows() ? ";" : ":")

        for pathdir in pathdirs
            local actualpath
            try
                actualpath = realpath(pathdir)
            catch
                # Bash doesn't expect every folder in PATH to exist, so neither shall we
                continue
            end

            if actualpath != pathdir && in(actualpath,pathdirs)
                # Remove paths which (after resolving links) are in the env path twice.
                # Many distros eg. point /bin to /usr/bin but have both in the env path.
                continue
            end

            local filesinpath
            try
                filesinpath = readdir(pathdir)
            catch e
                # Bash allows dirs in PATH that can't be read, so we should as well.
                if isa(e, SystemError)
                    continue
                else
                    # We only handle SystemErrors here
                    rethrow()
                end
            end

            for file in filesinpath
                # In a perfect world, we would filter on whether the file is executable
                # here, or even on whether the current user can execute the file in question.
                isfile(joinpath(pathdir, file)) && push!(matches, file)
            end
        end
    end

    matchList = Completion[PathCompletion(shell_escape ? replace(s, r"\s" => s"\\\0") : s, prefix) for s in matches]
    startpos = pos - lastindex(prefix) + 1 - count(isequal(' '), prefix)
    # The pos - lastindex(prefix) + 1 is correct due to `lastindex(prefix)-lastindex(prefix)==0`,
    # hence we need to add one to get the first index. This is also correct when considering
    # pos, because pos is the `lastindex` a larger string which `endswith(path)==true`.
    return sort_suggestions!(matchList), startpos:pos, !isempty(matchList)
end

function complete_expanduser(path::AbstractString, r)
    expanded = expanduser(path)
    return Completion[PathCompletion(expanded)], r, path != expanded
end

import REPL.REPLCompletions:
    should_method_complete,
    find_start_brace,
    get_value,
    get_type

# Method completion on function call expression that look like :(max(1))
function complete_methods(ex_org::Expr, context_module::Module=Main)
    func, found = get_value(ex_org.args[1], context_module)::Tuple{Any,Bool}
    !found && return Completion[]

    args_ex, kwargs_ex = complete_methods_args(ex_org.args[2:end], ex_org, context_module, true, true)

    out = Completion[]
    complete_methods!(out, func, args_ex, kwargs_ex)
    return out
end

function complete_any_methods(ex_org::Expr, callee_module::Module, context_module::Module, moreargs::Bool)
    out = Completion[]
    args_ex, kwargs_ex = try
        complete_methods_args(ex_org.args[2:end], ex_org, context_module, false, false)
    catch
        return out
    end

    for name in names(callee_module; all=true)
        if !Base.isdeprecated(callee_module, name) && isdefined(callee_module, name)
            func = getfield(callee_module, name)
            if !isa(func, Module)
                complete_methods!(out, func, args_ex, kwargs_ex, moreargs)
            elseif callee_module === Main::Module && isa(func, Module)
                callee_module2 = func
                for name in names(callee_module2)
                    if isdefined(callee_module2, name)
                        func = getfield(callee_module, name)
                        if !isa(func, Module)
                            complete_methods!(out, func, args_ex, kwargs_ex, moreargs)
                        end
                    end
                end
            end
        end
    end

    return out
end

function complete_methods_args(funargs::Vector{Any}, ex_org::Expr, context_module::Module, default_any::Bool, allow_broadcasting::Bool)
    args_ex = Any[]
    kwargs_ex = Pair{Symbol,Any}[]
    if allow_broadcasting && ex_org.head === :. && ex_org.args[2] isa Expr
        # handle broadcasting, but only handle number of arguments instead of
        # argument types
        for _ in (ex_org.args[2]::Expr).args
            push!(args_ex, Any)
        end
    else
        for ex in funargs
            if isexpr(ex, :parameters)
                for x in ex.args
                    n, v = isexpr(x, :kw) ? (x.args...,) : (x, x)
                    push!(kwargs_ex, n => get_type(get_type(v, context_module)..., default_any))
                end
            elseif isexpr(ex, :kw)
                n, v = (ex.args...,)
                push!(kwargs_ex, n => get_type(get_type(v, context_module)..., default_any))
            else
                push!(args_ex, get_type(get_type(ex, context_module)..., default_any))
            end
        end
    end
    return args_ex, kwargs_ex
end

function complete_methods!(out::Vector{Completion}, @nospecialize(func), args_ex::Vector{Any}, kwargs_ex::Vector{Pair{Symbol,Any}}, moreargs::Bool=true)
    ml = methods(func)
    # Input types and number of arguments
    if isempty(kwargs_ex)
        t_in = Tuple{Core.Typeof(func), args_ex...}
        na = length(t_in.parameters)::Int
        orig_ml = fill(nothing, length(ml))
    else
        isdefined(ml.mt, :kwsorter) || return out # kwargs given, but no kwarg method found
        kwfunc = ml.mt.kwsorter
        kwargt = NamedTuple{(first.(kwargs_ex)...,), Tuple{last.(kwargs_ex)...}}
        t_in = Tuple{Core.Typeof(kwfunc), kwargt, Core.Typeof(func), args_ex...}
        na = length(t_in.parameters)::Int
        orig_ml = ml
        ml = methods(kwfunc)
        func = kwfunc
    end
    if !moreargs
        na = typemax(Int)
    end

    for (method::Method, orig_method) in zip(ml, orig_ml)
        ms = method.sig

        # Check if the method's type signature intersects the input types
        if typeintersect(Base.rewrap_unionall(Tuple{(Base.unwrap_unionall(ms)::DataType).parameters[1 : min(na, end)]...}, ms), t_in) != Union{}
            push!(out, MethodCompletion(func, t_in, method, orig_method))
        end
    end
end

import REPL.REPLCompletions:
    latex_symbols,
    emoji_symbols,
    non_identifier_chars,
    whitespace_chars,
    bslash_separators,
    subscripts,
    subscript_regex,
    superscripts,
    superscript_regex,
    afterusing

function close_path_completion(str, startpos, r, paths, pos)
    length(paths) == 1 || return false  # Only close if there's a single choice...
    _path = str[startpos:prevind(str, first(r))] * (paths[1]::PathCompletion).path
    path = expanduser(replace(_path, r"\\ " => " "))
    # ...except if it's a directory...
    try
        isdir(path)
    catch e
        e isa Base.IOError || rethrow() # `path` cannot be determined to be a file
    end && return false
    # ...and except if there's already a " at the cursor.
    return lastindex(str) <= pos || str[nextind(str, pos)] != '"'
end

function bslash_completions(string, pos)
    slashpos = something(findprev(isequal('\\'), string, pos), 0)
    if (something(findprev(in(bslash_separators), string, pos), 0) < slashpos &&
        !(1 < slashpos && (string[prevind(string, slashpos)]=='\\')))
        # latex / emoji symbol substitution
        s = string[slashpos:pos]
        latex = get(latex_symbols, s, "")
        if !isempty(latex) # complete an exact match
            return (true, (Completion[BslashCompletion(latex)], slashpos:pos, true))
        elseif occursin(subscript_regex, s)
            sub = map(c -> subscripts[c], s[3:end])
            return (true, (Completion[BslashCompletion(sub)], slashpos:pos, true))
        elseif occursin(superscript_regex, s)
            sup = map(c -> superscripts[c], s[3:end])
            return (true, (Completion[BslashCompletion(sup)], slashpos:pos, true))
        end
        emoji = get(emoji_symbols, s, "")
        if !isempty(emoji)
            return (true, (Completion[BslashCompletion(emoji)], slashpos:pos, true))
        end
        # return possible matches; these cannot be mixed with regular
        # Julian completions as only latex / emoji symbols contain the leading \
        suggestions = Completion[BslashCompletion(k, s) for k in keys(startswith(s, "\\:") ? emoji_symbols : latex_symbols)]
        return (true, (sort_suggestions!(suggestions), slashpos:pos, true))
    end
    return (false, (Completion[], 0:-1, false))
end

function dict_identifier_key(str, tag, context_module = Main)
    if tag === :string
        str_close = str*"\""
    elseif tag === :cmd
        str_close = str*"`"
    else
        str_close = str
    end

    frange, end_of_identifier = find_start_brace(str_close, c_start='[', c_end=']')
    isempty(frange) && return (nothing, nothing, nothing)
    obj = context_module
    for name in split(str[frange[1]:end_of_identifier], '.')
        Base.isidentifier(name) || return (nothing, nothing, nothing)
        sym = Symbol(name)
        isdefined(obj, sym) || return (nothing, nothing, nothing)
        obj = getfield(obj, sym)
    end
    (isa(obj, AbstractDict) && length(obj)::Int < 1_000_000) || return (nothing, nothing, nothing)
    begin_of_key = something(findnext(!isspace, str, nextind(str, end_of_identifier) + 1), # +1 for [
                             lastindex(str)+1)
    return (obj::AbstractDict, str[begin_of_key:end], begin_of_key)
end

# This needs to be a separate non-inlined function, see #19441
@noinline find_dict_matches(identifier) = String[repr(key) for key in keys(identifier)]

@static if isdefined(Base, :parsed_toml)

function project_deps_get_completion_candidates(pkgstarts::String, project_file::String)
    loading_candidates = String[]
    d = Base.parsed_toml(project_file)
    pkg = get(d, "name", nothing)::Union{String, Nothing}
    if pkg !== nothing && startswith(pkg, pkgstarts)
        push!(loading_candidates, pkg)
    end
    deps = get(d, "deps", nothing)::Union{Dict{String, Any}, Nothing}
    if deps !== nothing
        for (pkg, _) in deps
            startswith(pkg, pkgstarts) && push!(loading_candidates, pkg)
        end
    end
    return Completion[PackageCompletion(name) for name in loading_candidates]
end

else # @static if isdefined(Base, :TOML)

function project_deps_get_completion_candidates(pkgstarts::String, project_file::String)
    loading_candidates = String[]
    open(project_file) do io
        state = :top
        for line in eachline(io)
            if occursin(Base.re_section, line)
                state = occursin(Base.re_section_deps, line) ? :deps : :other
            elseif state === :top
                if (m = match(Base.re_name_to_string, line)) !== nothing
                    root_name = String(m.captures[1])
                    push!(loading_candidates, root_name)
                end
            elseif state === :deps
                if (m = match(Base.re_key_to_string, line)) !== nothing
                    dep_name = m.captures[1]
                    push!(loading_candidates, dep_name)
                end
            end
        end
    end
    return Completion[PackageCompletion(name, pkgstarts) for name in loading_candidates]
end

end # @static if isdefined(Base, :TOML)

function completions(string, pos, context_module = Main)
    # First parse everything up to the current position
    partial = string[1:pos]
    inc_tag = Base.incomplete_tag(Meta.parse(partial, raise=false, depwarn=false))

    # ?(x, y)TAB lists methods you can call with these objects
    # ?(x, y TAB lists methods that take these objects as the first two arguments
    # MyModule.?(x, y)TAB restricts the search to names in MyModule
    rexm = match(r"(\w+\.|)\?\((.*)$", partial)
    if rexm !== nothing
        # Get the module scope
        if isempty(rexm.captures[1])
            callee_module = context_module
        else
            modname = Symbol(rexm.captures[1][1:end-1])
            if isdefined(context_module, modname)
                callee_module = getfield(context_module, modname)
                if !isa(callee_module, Module)
                    callee_module = context_module
                end
            else
                callee_module = context_module
            end
        end
        moreargs = !endswith(rexm.captures[2], ')')
        callstr = "_(" * rexm.captures[2]
        if moreargs
            callstr *= ')'
        end
        ex_org = Meta.parse(callstr, raise=false, depwarn=false)
        if isa(ex_org, Expr)
            return complete_any_methods(ex_org, callee_module::Module, context_module, moreargs), (0:length(rexm.captures[1])+1) .+ rexm.offset, false
        end
    end

    # if completing a key in a Dict
    identifier, partial_key, loc = dict_identifier_key(partial,inc_tag, context_module)
    if identifier !== nothing
        matches = find_dict_matches(identifier)::Vector{String}
        length(matches)==1 && (lastindex(string) <= pos || string[nextind(string,pos)] != ']') && (matches[1]*=']')
        if length(matches)>0
            suggestions = Completion[DictCompletion(identifier, match, partial_key) for match in matches]
            return sort_suggestions!(suggestions), loc:pos, false
        end
    end

    # otherwise...
    if inc_tag in [:cmd, :string]
        m = match(r"[\t\n\r\"`><=*?|]| (?!\\)", reverse(partial))
        startpos = nextind(partial, reverseind(partial, m.offset))
        r = startpos:pos

        expanded = complete_expanduser(replace(string[r], r"\\ " => " "), r)
        expanded[3] && return expanded  # If user expansion available, return it

        paths, r, success = complete_path(replace(string[r], r"\\ " => " "), pos)

        if inc_tag === :string && close_path_completion(string, startpos, r, paths, pos)
            paths[1] = PathCompletion((paths[1]::PathCompletion).path * "\"")
        end

        # within string scope, append latex completions as well
        if inc_tag === :string
            should_bs, bs_ret = bslash_completions(string, pos)
            if should_bs
                sort_suggestions!(append!(first(bs_ret), paths))
                return bs_ret
            end
        end
        #Latex symbols can be completed for strings
        (success || inc_tag === :cmd) && return paths, r, success
    end

    ok, ret = bslash_completions(string, pos)
    ok && return ret

    # Make sure that only bslash_completions is working on strings
    inc_tag==:string && return Completion[], 0:-1, false
    if inc_tag === :other && should_method_complete(partial)
        frange, method_name_end = find_start_brace(partial)
        # strip preceding ! operator
        s = replace(partial[frange], r"\!+([^=\(]+)" => s"\1")
        ex = Meta.parse(s * ")", raise=false, depwarn=false)

        if isa(ex, Expr)
            if ex.head === :call
                return complete_methods(ex, context_module), first(frange):method_name_end, false
            elseif ex.head === :. && begin
                    x = ex.args[2]
                    isa(x, Expr) && x.head === :tuple
                end
                return complete_methods(ex, context_module), first(frange):(method_name_end - 1), false
            end
        end
    elseif inc_tag === :comment
        return Completion[], 0:-1, false
    end

    dotpos = something(findprev(isequal('.'), string, pos), 0)
    startpos = nextind(string, something(findprev(in(non_identifier_chars), string, pos), 0))
    # strip preceding ! operator
    if (m = match(r"^\!+", string[startpos:pos])) !== nothing
        startpos += length(m.match)
    end

    ffunc = (mod,x)->true
    suggestions = Completion[]
    comp_keywords = true
    if afterusing(string, startpos)
        # We're right after using or import. Let's look only for packages
        # and modules we can reach from here

        # If there's no dot, we're in toplevel, so we should
        # also search for packages
        s = string[startpos:pos]
        if dotpos <= startpos
            for dir in Base.load_path()
                if basename(dir) in Base.project_names && isfile(dir)
                    append!(suggestions, project_deps_get_completion_candidates(s, dir))
                end
                isdir(dir) || continue
                for pname in readdir(dir)
                    if pname[1] != '.' && pname != "METADATA" && pname != "REQUIRE"
                        # Valid file paths are
                        #   <Mod>.jl
                        #   <Mod>/src/<Mod>.jl
                        #   <Mod>.jl/src/<Mod>.jl
                        if isfile(joinpath(dir, pname))
                            endswith(pname, ".jl") && push!(suggestions,
                                                            PackageCompletion(pname[1:prevind(pname, end-2)], s))
                        else
                            mod_name = if endswith(pname, ".jl")
                                pname[1:prevind(pname, end-2)]
                            else
                                pname
                            end
                            if isfile(joinpath(dir, pname, "src",
                                               "$mod_name.jl"))
                                push!(suggestions, PackageCompletion(mod_name, s))
                            end
                        end
                    end
                end
            end
        end
        ffunc = (mod,x)->(Base.isbindingresolved(mod, x) && isdefined(mod, x) && isa(getfield(mod, x), Module))
        comp_keywords = false
    end
    startpos == 0 && (pos = -1)
    dotpos < startpos && (dotpos = startpos - 1)
    s = string[startpos:pos]
    comp_keywords && append!(suggestions, complete_keyword(s))
    # The case where dot and start pos is equal could look like: "(""*"").d","". or  CompletionFoo.test_y_array[1].y
    # This case can be handled by finding the beginning of the expression. This is done below.
    if dotpos == startpos
        i = prevind(string, startpos)
        while 0 < i
            c = string[i]
            if c in [')', ']']
                if c==')'
                    c_start='('; c_end=')'
                elseif c==']'
                    c_start='['; c_end=']'
                end
                frange, end_of_identifier = find_start_brace(string[1:prevind(string, i)], c_start=c_start, c_end=c_end)
                startpos = first(frange)
                i = prevind(string, startpos)
            elseif c in ['\'', '\"', '\`']
                s = "$c$c"*string[startpos:pos]
                break
            else
                break
            end
            s = string[startpos:pos]
        end
    end
    append!(suggestions, complete_symbol(s, ffunc, context_module))
    return sort_suggestions!(suggestions), (dotpos+1):pos, false
end

@inline sort_suggestions!(suggestions) = sort!(suggestions, by=score, rev=true)

function shell_completions(string, pos)
    # First parse everything up to the current position
    scs = string[1:pos]
    local args, last_parse
    try
        args, last_parse = Base.shell_parse(scs, true)::Tuple{Expr,UnitRange{Int}}
    catch
        return Completion[], 0:-1, false
    end
    ex = args.args[end]::Expr
    # Now look at the last thing we parsed
    isempty(ex.args) && return Completion[], 0:-1, false
    arg = ex.args[end]
    if all(s -> isa(s, AbstractString), ex.args)
        arg = arg::AbstractString
        # Treat this as a path

        # As Base.shell_parse throws away trailing spaces (unless they are escaped),
        # we need to special case here.
        # If the last char was a space, but shell_parse ignored it search on "".
        ignore_last_word = arg != " " && scs[end] == ' '
        prefix = ignore_last_word ? "" : join(ex.args)

        # Also try looking into the env path if the user wants to complete the first argument
        use_envpath = !ignore_last_word && length(args.args) < 2

        return complete_path(prefix, pos, use_envpath=use_envpath, shell_escape=true)
    elseif isexpr(arg, :incomplete) || isexpr(arg, :error)
        partial = scs[last_parse]
        ret, range = completions(partial, lastindex(partial))
        range = range .+ (first(last_parse) - 1)
        return ret, range, true
    end
    return Completion[], 0:-1, false
end

end # module


text: # just to make package manager happy

ext: jl


text:
using MXNet
# main training function
function TrainNet(data, trainsize, noutputs, layerconfig, batchsize, epochs, savefile)
    # prepare data
    data, mX, sX = stnorm(data)
    Y = data[1:trainsize,1:noutputs]'
    YT = data[trainsize+1:end,1:noutputs]'
    X = data[1:trainsize,noutputs+1:end]'
    XT = data[trainsize+1:end,noutputs+1:end]'
    # sizes of layers
    L1size = layerconfig[1]
    L2size = layerconfig[2]
    L3size = layerconfig[3]
    L4size = layerconfig[4]
    Outputsize = size(Y,1)
    # set up MLP
    data = mx.Variable(:data)
    label = mx.Variable(:label)
    if L4size != 0
        net  = @mx.chain    mx.Variable(:data) =>
                            mx.FullyConnected(num_hidden=L1size) =>
                            mx.Activation(act_type=:relu) =>
                            mx.FullyConnected(num_hidden=L2size) =>
                            mx.Activation(act_type=:tanh) =>
                            mx.FullyConnected(num_hidden=L3size) =>
                            mx.Activation(act_type=:relu) =>
                            mx.FullyConnected(num_hidden=L4size) =>
                            mx.Activation(act_type=:tanh) =>
                            mx.FullyConnected(num_hidden=Outputsize) =>
                            mx.LinearRegressionOutput(label)
    elseif L3size != 0
        net  = @mx.chain    mx.Variable(:data) =>
                            mx.FullyConnected(num_hidden=L1size) =>
                            mx.Activation(act_type=:relu) =>
                            mx.FullyConnected(num_hidden=L2size) =>
                            mx.Activation(act_type=:tanh) =>
                            mx.FullyConnected(num_hidden=L3size) =>
                            mx.Activation(act_type=:relu) =>
                            mx.FullyConnected(num_hidden=Outputsize) =>
                            mx.LinearRegressionOutput(label)
    else
        net  = @mx.chain    mx.Variable(:data) =>
                            mx.FullyConnected(num_hidden=L1size) =>
                            mx.Activation(act_type=:relu) =>
                            mx.FullyConnected(num_hidden=L2size) =>
                            mx.Activation(act_type=:tanh) =>
                            mx.FullyConnected(name=:fc2, num_hidden=Outputsize) =>
                            mx.LinearRegressionOutput(label)
    end
    # final model definition, don't change, except if using gpu
    model = mx.FeedForward(net, context=mx.cpu())
    # set up the optimizer: select one, explore parameters, if desired
    #optimizer = mx.SGD(lr=0.01, momentum=0.9, weight_decay=0.00001)
    optimizer = mx.ADAM()
    trainprovider = mx.ArrayDataProvider(:data => X, batch_size=batchsize, shuffle=true, :label => Y)
    evalprovider = mx.ArrayDataProvider(:data => XT, batch_size=batchsize, shuffle=true, :label => YT)
    # train
    mx.fit(model, optimizer, eval_metric=mx.MSE(), initializer=mx.UniformInitializer(0.05), trainprovider, eval_data=evalprovider, n_epoch = epochs)
    # more training with larger batch size, saving the final fitted model
    batchsize = 10*batchsize
    mx.fit(model, optimizer, eval_metric=mx.MSE(), trainprovider, eval_data=evalprovider, n_epoch = epochs, callbacks=[mx.do_checkpoint(savefile, frequency=epochs, save_epoch_0=false)])
end

max_stars_repo_name: nilshg/Econometrics

max_issues_repo_name: nilshg/Econometrics

max_forks_repo_name: nilshg/Econometrics


text:
# this file used something like:

# julia-0.4 -e 'using Go; e = parse_log("evol-1.log", "evol-1.dump", 50, 0.5; limit=-1); plot_tramlines(e, "go.png"; min_axis=1900, tint=ancestors("05f051414412588b"))'

# where the ancestors tag is typically a recent, high scoring individual

const p_challenge = r"^\s*(?P<i>\d+)/(?P<n>\d+)\s+(?P<j>\d+)/(?P<m>\d+) (?P<surprise>(?:!| )) (?P<a>[a-f0-9]+):(?P<x>\d+)\s+(?:>|=|<)\s+(?P<b>[a-f0-9]+):(?P<y>\d+).*$"

# 34 random bytes d1862e709058196a => fddd8634d3f25690
# 30 random bits 755e948bc02367be => d68ad2353d0c378a
const p_random = r"^\s*\d+ random (?:bits|bytes|deltas) (?P<a>[a-f0-9]+) => (?P<c>[a-f0-9]+)\s*$"

# rotate 262 cb03e03f73d749c2 => 326e568542888c67 + merge 123 326e568542888c67, fddd8634d3f25690 => f2b69ec0fe5c1fd7
# rotate chunk 7 0043f8c791de1a53 => 0bf67111a514bd9d + merge 97 f7e67544a9a5ef3c, 0bf67111a514bd9d => da7fc58918b190c0
const p_rotate = r"^\s*rotate (?:chunk )?\d+ (?P<a>[a-f0-9]+) => (?P<x>[a-f0-9]+) \+ merge \d+ (?P<y>[a-f0-9]+), (?P<z>[a-f0-9]+) => (?P<c>[a-f0-9]+)\s*$"

# merge 264 0f9139cba0588c2a, d9a5f3b8d06028d9 => 9b6504a9fe224d56
const p_merge = r"^\s*merge \d+ (?P<a>[a-f0-9]+), (?P<b>[a-f0-9]+) => (?P<c>[a-f0-9]+)\s*$"

# temp 0.819
const p_temp = r"\s*temp \d+\.\d+\s*$"

const INK = C.RGB(0.2, 0.2, 0.2)
const PEN = 0.45

# this is not hashed - each instance is unique, and we can have
# multiple, distinct instances with the same tag.
type Individual
    tag::AbstractString
end

abstract Event

immutable VirginBirths <: Event
    ids::Vector{Individual}
end

type Birth <: Event
    id::Individual
    parents::Vector{Individual}
    deaths::Vector{Individual}
    Birth(id, parents) = new(id, parents, Individual[])
end

abstract Challenge <: Event

immutable FailedChallenge <: Challenge
    won::Individual
    lost::Individual
end

type SuccessfulChallenge <: Challenge
    won::Individual
    lost::Individual
    deaths::Vector{Individual}
    SuccessfulChallenge(won, lost) = new(won, lost, Individual[])
end

immutable Deaths <: Event
    ids::Vector{Individual}
    Deaths(ids::Vector{Individual}) = new(ids)
    Deaths(id::Individual) = new(Individual[id])
end

delete!(v::Vector, i::Integer) = splice!(v, i:i)

function prune!(events, id)

    death = Nullable{Event}()
    function clear_death()
        if ! isnull(death)
            delete!(get(death).deaths, findfirst(get(death).deaths, id))
        end
    end
    function set_death(e)
        clear_death()
        findfirst(e.deaths, id) == 0 && push!(e.deaths, id)
        death = Nullable{Event}(e)
    end

    deleted, positive, i, dies = 0, false, length(events), false
    for i = length(events):-1:1
        e = events[i]
        if isa(e, VirginBirths) && id in e.ids
            if ! positive && dies
                delete!(e.ids, findfirst(e.ids, id))
                deleted += 1
                clear_death()
            end
            return deleted
        elseif isa(e, Birth) && e.id == id
            if ! positive # && dies  so we clear out near start too
                events[i] = Deaths(e.deaths)
                deleted += 1
                clear_death()
            end
            return deleted
        elseif isa(e, Birth) && ! positive
            if dies
                set_death(e)
            end
            positive = id in e.parents
        elseif isa(e, SuccessfulChallenge) && e.won == id
            positive = true
        elseif isa(e, SuccessfulChallenge) && e.lost == id
            if ! positive && dies
                set_death(e)
            end
            positive = true
        elseif isa(e, SuccessfulChallenge) && ! positive && dies
            set_death(e)
        elseif isa(e, Deaths) && id in e.ids
            dies = true
        end
    end
    deleted
end

allids(e::VirginBirths) = e.ids
allids(e::Birth) = vcat(e.id, e.parents, e.deaths)
allids(e::Challenge) = vcat(e.won, e.lost, e.deaths)
allids(e::Deaths) = e.ids

function prune!(events)
    deleted, ids, i = 0, Set{Individual}(), length(events)
    while i > 0
        idcpy = collect(allids(events[i]))
        for id in idcpy
            if ! (id in ids)
                push!(ids, id)
                deleted += prune!(events, id)
            end
        end
        i = min(i-1, length(events))
    end
    println("deleted $(deleted) (from unsuccessful individuals)")
    deleted
end

function prune_type!(events, T)
    count = 0
    for i = length(events):-1:1
        if isa(events[i], T)
            delete!(events, i)
            count += 1
        end
    end
    println("deleted $(count) $(T)")
    events
end

parent(population, tag) = population[findfirst(id -> id.tag == tag, population)]

function parse_line(events, population, line, n, fraction)

#    println("$(length(population)) $(line)")

    if (m = match(p_challenge, line)) != nothing
        x, y = parse(Int, m[:x]), parse(Int, m[:y])
        a, b = population[x], population[y]
         @assert a.tag == m[:a]
        @assert b.tag == m[:b]
        # order so x/a is the higher rated
        if y < x
            x, y = y, x
            a, b = b, a
        end
        if m[:surprise] == "!"
            push!(events,  SuccessfulChallenge(b, a))
            population[x+1:y] = population[x:y-1]
            population[x] = b
        else
            push!(events, FailedChallenge(a, b))
        end
        return 1

    elseif (m = match(p_random, line)) != nothing
        parent1 = parent(population, m[:a])  # triggers error on initial junk
        push!(population, Individual(m[:c]))
        push!(events, Birth(population[end], [parent1]))
        return 1

    elseif (m = match(p_rotate, line)) != nothing
        parent1 = parent(population, m[:a])
        parent2 = parent(population, m[:x] == m[:y] ? m[:z] : m[:y])
        push!(population, Individual(m[:c]))
        push!(events, Birth(population[end], [parent1, parent2]))
        return 1

    elseif (m = match(p_merge, line)) != nothing
        parent1 = parent(population, m[:a])
        parent2 = parent(population, m[:b])
        push!(population, Individual(m[:c]))
        push!(events, Birth(population[end], [parent1, parent2]))
        return 1

    elseif match(p_temp, line) != nothing
        push!(events, Deaths(Individual[]))
        while length(population) > Int(fraction * n)
            id = pop!(population)
            push!(events[end].ids, id)
        end
        return 1

    else
        return 0
    end
end

function parse_log(log_path, dump_path, n, fraction; limit=-1)

    events, population = Vector{Event}(), Vector{Individual}()
    push!(events, VirginBirths(Individual[]))

    # take initial population from dump
    for expression in undump(dump_path)[2][1:n]
        push!(population, Individual(name(expression)))
        push!(events[1].ids, population[end])
    end

    open(log_path, "r") do io

        started = false

        for line in eachline(io)

            try
                if ! started
                    # discard any junk breeding
                    population = population[1:n]
                end
                k = parse_line(events, population, line[1:end-1], n, fraction)
                started = started || isa(events[end], Challenge)
                if started
                    limit -= k
                    limit == 0 && break
                end
            catch
                if started
                    rethrow()
                end
            end

        end
    end

    println("initial $(length(events)) events")
    prune_type!(events, FailedChallenge)
    while prune!(events) > 0 end
    prune_type!(events, Deaths)
    println("final $(length(events)) events")

    events
end

typealias ColourTable Tuple{Int, C.Colorant}

type PlotState
    max_pop::Int
    rows::Int
    cols::Int
    row::Int
    col::Int
    colour_index::Dict{Individual, Tuple{Int,Int,Int}}
    colour_table::Dict{Int, C.Colorant}
    PlotState(max_pop, rows, cols) = new(max_pop, rows, cols, 0, 0,
                                         Dict{Individual, Tuple{Int,Int,Int}}(),
                                         Dict{Int, C.Colorant}())
end

show(io::IO, s::PlotState) = show(io, "PlotState($(s.max_pop),$(s.col)/$(s.cols),$(s.row)/$(s.rows))")

function expand!(p::Vector{Vector{Individual}}, e::Birth)
    popn = copy(p[end])
    push!(popn, e.id)
    for id in e.deaths
        i = findfirst(popn, id)
        i > 0 && delete!(popn, i)
    end
    push!(p, popn)
end

function expand!(p::Vector{Vector{Individual}}, e::VirginBirths)
    popn = isempty(p) ? Individual[] : copy(p[end])
    for id in e.ids
        push!(popn, id)
    end
    push!(p, popn)
end

function expand!(p::Vector{Vector{Individual}}, e::SuccessfulChallenge)
    popn = copy(p[end])
    i = findfirst(popn, e.lost)
    j = findfirst(popn, e.won)
    @assert i < j (i, j)
    popn[i+1:j] = popn[i:j-1]
    popn[i] = e.won
    for id in e.deaths
        i = findfirst(popn, id)
        i > 0 && delete!(popn, i)
    end
    push!(p, popn)
end

function expand!(p::Vector{Vector{Individual}}, e::FailedChallenge)
    popn = copy(p[end])
    push!(p, popn)
end

function make_state(events, ny, ratio, min_scale, min_axis)
    nx = length(events)
    rows = max(1, floor(Int, sqrt(nx * ny / ratio) / ny))
    cols = ceil(Int, nx / rows) + 1
    s = PlotState(ny, rows, cols)
    n = max(min_scale * cols, min_axis / 1.2)
    scale = n / cols
    m = scale * ny * rows + 0.2 * n
    n = n + 0.2 * n
    println("($(nx),$(ny)) -> $s, ($n, $m), $(scale)")
    s, n, m
end

function break!(e, p, s::PlotState)
    for i in 1:(s.rows-1)
        j = i*s.cols
        splice!(e, j:j-1, [VirginBirths(Individual[])])
        splice!(p, j:j-1, Vector{Individual}[p[j]])
    end
end

function getc(s, id, k)
    (b, kk, a) = s.colour_index[id]
    if k <= kk
        s.colour_table[b]
    else
        s.colour_table[a]
    end
end

function grey!(s, popn, events)
    s.colour_table[0] = INK
    for p in popn
        for id in p
            if ! haskey(s.colour_index, id)
                s.colour_index[id] = (0, 0, 0)
            end
        end
    end
end

function multi(s, popn, events)
    for p in popn
        for id in p
            if ! haskey(s.colour_index, id)
                n = length(s.colour_table) + 1
                s.colour_table[n] = convert(C.RGB{N.UFixed8}, C.HSV(360 * rand(), 0.5 + 0.5 * rand(), 0.2 + 0.8 * rand()))
                s.colour_index[id] = (n, 0, n)
            end
        end
    end
end

function multig(s, popn, events)
    for p in popn
        for id in p
            if ! haskey(s.colour_index, id)
                n = length(s.colour_table) + 1
                s.colour_table[n] = convert(C.RGB{N.UFixed8}, C.HSV(0, 0, 0.8 * rand()))
                s.colour_index[id] = (n, 0, n)
            end
        end
    end
end

function find(tag, popn)
    for k in length(popn):-1:1
        for id in popn[k]
            if id.tag == tag
                return id
            end
        end
    end
    error("$(tag) does not exist")
end

function ancestors(tag...)

    function tint!(s, popn, events)

        parents = Dict{Individual, Tuple{Int, Vector{Individual}}}()
        for (k, e) in enumerate(events)
            if isa(e, Birth)
                parents[e.id] = (k, e.parents)
            end
        end

#        grey!(s, popn, events)
        multig(s, popn, events)
        zero = length(s.colour_table)

        for (h, t) in enumerate(tag)
            println("$h $t")
            print("tinting...")
            qcount, pcount = 0, 0
            id = find(t, popn)
            queue = Tuple{Int, Int, Individual}[(1, length(events), id)]
            while length(queue) > 0
                (depth, k, id) = splice!(queue, 1:1)[1]
                index = h * 1000 + depth + zero
                if ! haskey(s.colour_table, index)
                    z = 0.95^(depth-1)
                    s.colour_table[index] = convert(C.RGB{N.UFixed8}, C.HSV((h-1) * 60.0, z, 0.2+0.8*z))
                    print(" $depth,")
                end
                c = s.colour_table[index]
                if s.colour_index[id][2] == 0
                    s.colour_index[id] = (index, k, s.colour_index[id][3])
                    pcount += 1
                end
                (k, ps) = get(parents, id, (0, []))
                for p in ps
                    # tweak for efficiency - queue was huge
                    if s.colour_index[p][2] == 0
                        push!(queue, (depth+1, k, p))
                        qcount += 1
                    end
                end
            end
            println(" done")
            println("$(qcount) queued; $(pcount) altered")
        end
    end
end

function children(tag...)

    function tint!(s, popn, events)

        children = Dict{Individual, Vector{Individual}}()
        for (k, e) in enumerate(events)
            if isa(e, Birth)
                for p in e.parents
                    children[p] = push!(get(children, p, Individual[]), e.id)
                end
            end
        end

        multig(s, popn, events)
        zero = length(s.colour_table)

        for (h, t) in enumerate(tag)
            println("$h $t")
            print("tinting...")
            qcount, pcount = 0, 0
            id = find(t, popn)
            queue = Tuple{Int, Individual}[(1, id)]
            while length(queue) > 0
                (depth, id) = splice!(queue, 1:1)[1]
                index = h * 1000 + depth + zero
                if ! haskey(s.colour_table, index)
                    z = 0.95^(depth-1)
                    s.colour_table[index] = convert(C.RGB{N.UFixed8}, C.HSV((h-1) * 60.0, z, 0.2+0.8*z))
                    print(" $depth,")
                end
                c = s.colour_table[index]
                if s.colour_index[id][2] == 0
                    s.colour_index[id] = (index, 1, index)
                    pcount += 1
                end
                for c in get(children, id, [])
                    # tweak for efficiency - queue was huge
                    if s.colour_index[c][2] == 0
                        push!(queue, (depth+1, c))
                        qcount += 1
                    end
                end
            end
            println(" done")
            println("$(qcount) queued; $(pcount) altered")
        end
    end
end

function rshift(popn)
    popn2 = copy(popn)
    splice!(popn2, 1:0, popn2[1:1])
    popn2[1] = Individual[]
    popn2
end

function squish!(events, popn)
    e, p = Integer[], Integer[]
    for i in 2:(length(events)-1)
        if isa(events[i], Birth) && isa(events[i+1], SuccessfulChallenge) && events[i].id == events[i+1].won
            push!(e, i+1)
            push!(p, i)
        end
    end
    deleteat!(events, e)
    deleteat!(popn, p)
    println("squished $(length(e)) births")
end

function plot_tramlines(events, path;
                        ratio=1.4, min_scale=2, min_axis=1000, tint=grey!,
                        lo=1, hi=-1)

    popn = Vector{Vector{Individual}}()
    for e in events
        expand!(popn, e)
    end

    squish!(events, popn)  # only for evolve2

    if lo > 1
        popn = popn[lo:end]
        events = events[lo+1:end]
        splice!(events, 1:0, [VirginBirths(popn[1])])
    end
    if hi != -1
        popn = popn[1:hi-lo+1]
        events = events[1:hi-lo+1]
    end
    max_pop = maximum(map(length, popn))

    s, nx, ny = make_state(events, max_pop, ratio, min_scale, min_axis)
    break!(events, popn, s)

    dk = ceil(Int, length(events) / 10)
    tint(s, popn, events)
    D.with(
           D.PNG(path, ceil(Int, nx), ceil(Int, ny)),
#           D.PDF(path; size="letter"),
           D.Axes(scale=(s.cols, s.rows*s.max_pop)),
           D.Pen(PEN; cap="round", join="round")) do
        print("plotting $(length(events))...")
        for (k, (b, e, a)) in enumerate(zip(rshift(popn), events, popn))
            if k % dk == 0
                print(" $k,")
            end
            plot(s, b, e, a, k)
        end
        println(" done")
    end
end

gety(s, i) = (s.rows - (s.row-1)) * s.max_pop - i

function plot_step(x, yb, ya)
    D.move(x, yb)
    if yb-ya > 0.1
        D.line(x+0.5, yb-0.5)
        D.line(x+0.5, ya+0.5)
        D.line(x+1, ya)
    elseif yb-ya < -0.1
        D.line(x+0.5, yb+0.5)
        D.line(x+0.5, ya-0.5)
        D.line(x+1, ya)
    else
        D.line(x+1, ya)
    end
end

function plot_steps(s, b, a, k)
    for id in reverse(a)
        c, i, j = getc(s, id, k), findfirst(b, id), findfirst(a, id)
        if i > 0 && j > 0
            x = s.col
            yb, ya = gety(s, i), gety(s, j)
            D.draw(D.Ink(c)) do
                plot_step(x, yb, ya)
            end
        end
    end
    s.col += 1
end

plot(s, b, e::SuccessfulChallenge, a, k) = plot_steps(s, b, a, k)

fade(c) = D.WHITE - 0.5 * (D.WHITE - c)

function plot(s, b, e::Birth, a, k)
    x = s.col
    ya = gety(s, findfirst(a, e.id))
    for id in e.parents
        yb = gety(s, findfirst(b, id))
        D.draw(D.Ink(fade(getc(s, id, k)))) do
            plot_step(x, yb, ya)
        end
    end
    plot_steps(s, b, a, k)
end

function plot(s, b, e::VirginBirths, a, k)
    s.row += 1
    s.col = 1
end

# d = longest_lived("evol-1.log", "evol-1.dump", 50, 0.5; n2=100)
# evolve(d, 19, 1000, 25, 10000, "evol-2.dump")

function longest_lived(log_path, dump_path, n1, fraction; limit=-1, n2=-1)

    events = parse_log(log_path, dump_path, n1, fraction; limit=limit)
    popn = Vector{Vector{Individual}}()
    for e in events
        expand!(popn, e)
    end

    n2 = n2 < 1 ? n1 : n2
    score = Dict{AbstractString, Float64}()
    for p in popn
        for (rank, id) in enumerate(p)
            if ! haskey(score, id.tag)
                score[id.tag] = 0.0
            end
            score[id.tag] += 1.0 / rank
        end
    end

    pairs = [(score[id], id) for id in keys(score)]
    sort!(pairs; by = x -> x[1], rev=true)

    # read dump fike a second time - inefficient, but used only once
    named = undump(dump_path)[1]
    names = map(x -> x[2], pairs[1:n2])
    expressions = [named[name] for name in names]

    for (i, (name, expression)) in enumerate(zip(names, expressions))
        println("$(name):$i")
        dump_expression(expression)
    end

    return expressions
end

# d = final_gen("evol-2.log", "evol-2.dump", 100, 0.9)
# evolve(d, 19, 1000, 25, 10000, "evol-3.dump")

function final_gen(log_path, dump_path, n, fraction; limit=-1)

    events = parse_log(log_path, dump_path, n, fraction; limit=limit)
    popn = Vector{Vector{Individual}}()
    for e in events
        expand!(popn, e)
    end

    tags = AbstractString[]
    for id in popn[end]
        if !(id.tag in tags)
            push!(tags, id.tag)
        end
    end

    named = undump(dump_path)[1]
    expressions = [named[tag] for tag in tags]

    temp = 1.0
    ops = build_ops(n, temp)
    while length(expressions) < n
        expression = weighted_rand(ops)(expressions)
        if !(expression in expressions)
            println("adding $(name(expression))")
            push!(expressions, expression)
        end
    end
    while length(expressions) > n
        expression = pop!(expressions)
        println("dropping $(name(expression))")
    end

    for (i, expression) in enumerate(expressions)
        println("$(name(expression)):$i")
        dump_expression(expression)
    end

    return expressions
end


text:
"""
    generate_complex_potential_source(; pos = ComplexQuantity((0.0, 0.0)m)::ComplexQuantity, massflowout = 1.0m²/s)
    -> f::Z->R

We find the velocity potential function of a source or sink by realizing that
´´´
    2π·r·u_r = q
       where u_r is radial velocity
             q is flow
        ↕
    u_r = q / (2π·r)
´´´
The velocity potential, ϕ, is the imaginary valued function found by integrating
velocity along the velocity gradient. The sign and endpoint of the integration
can be defined as we want. With the limitation above:
´´´
    u_r  = ∇ϕ = δu / δr
    ⇑
    ϕ = ∫ q / (2πr) dr = q / (2πr) · ∫ 1 / r dr = q / (2πr) · ln(r)
´´´
"""
function generate_complex_potential_source(; pos = complex(0.0, 0.0)m, massflowout = 1.0m²/s)
    # Note: log means natural logarithm, not log10
    #
    # Note: natural logarithm of a quantity gives no meaning. We can drop the unit as long as we
    # use the scalar consistently.
    @assert massflowout != zero(massflowout)
    # Would call this ϕ_s(p)
    (p::ComplexQuantity) -> begin
        Δp = p - pos
        r = hypot(Δp)
        massflowout ∙ log(r / oneunit(r)) / (2π)
    end
end

"""
    generate_complex_potential_vortex(; pos = ComplexQuantity((0.0, 0.0)m)::ComplexQuantity, vorticity= 1.0m²/s)
    -> f::Z->R

We find the velocity potential function of a vortex by realizing that velocity
´´´
    u_θ = K / r
       where u_θ is tangential velocity, i.e. length per time in the tangential direction.
             K is vorticity
             Γ = K / 2π  is circulation
´´´
The velocity potential, ϕ, is the scalar valued function found by integrating
velocity along the velocity gradient. The sign and endpoint of the integration
can be defined as we want.
"""
function generate_complex_potential_vortex(; pos = complex(0.0, 0.0)m, vorticity = 1.0m²/s)
    # Would call this ϕ_v(p)
    (p::ComplexQuantity) -> begin
        Δp = p - pos
        Θ = -angle(Δp)
        vorticity ∙ Θ
    end
end

absolute_scale() = ColorSchemes.linear_grey_10_95_c0_n256
complex_arg0_scale() = ColorSchemes.linear_ternary_red_0_50_c52_n256

"Scale and move real values to fit in the closed interval 0.0 to 1.0, given minimum and maximum values in A"
normalize_datarange_real(mi, ma, A) = map(A) do x
    @assert ma != mi
    if isnan(x) || isinf(x)
        mi
    else
        (x - mi) / (ma - mi)
    end
end
"""
Scale and move the magnitude (modulus, absolute value) of complex values
to fit in the closed interval 0.0 to 1.0, without changing the argument (angle).
The argument of the minimum value(s) is lost, since the angle is undefined for zero.
"""
normalize_datarange_complex(mi, ma, A) = map(A) do x
    @assert ma != mi
    if isnan(x) || isinf(x)
        mi
    else
        θ = angle(x)
        r = (hypot(x) - mi) / (ma - mi)
        r * complex(cos(θ), sin(θ))
    end
end
"""
    normalize_datarange(A)

normalize the elements in matrix A from 0.0 to 1.0, type preserving type.
This includes complex values, in which case the magnitude is normalized.
"""
function normalize_datarange(A::AbstractArray{<:RealQuantity})
    # Minimum and maximum value
    mi, ma = lenient_min_max(A)
    normalize_datarange_real(mi, ma, A)
end
function normalize_datarange(A::AbstractArray{<:Real})
    # Minimum and maximum value
    mi, ma = lenient_min_max(A)
    normalize_datarange_real(mi, ma, A)
end
function normalize_datarange(A)
    # Minimum and maximum magnitude
    mi, ma = lenient_min_max(A)
    normalize_datarange_complex(mi, ma, A)
end


text:
using Test
import Weave: WeaveMarkdown
import Markdown

# Test markdown2html writer

html = WeaveMarkdown.html(Markdown.parse("""

# H1

## H2

## H3

Some **text** with different [^note] *formatting* and \$math\$ and text.

`some code` with [link](http://github.com)

[^note]: test note

---

\$more math\$

* List one
* List two

1. List one
2. List two

```julia
x = 3
```

!!! note "Something"

    Test admonition with ![Image](link/to/image.png)

> Some important quote

""", flavor = WeaveMarkdown.weavemd))

ref_html = """
H1

H2

H3

Some text with different [note] formatting and \$math\$ and text.

some code with link

note

test note

\\[
more math
\\]

List one

List two

List one

List two

x = 3

Something

Test admonition with

Some important quote

"""

@test html == ref_html

#Test Weave additions
md = Markdown.parse("""

Multiline equations

\$\$
x = 2
\$\$

And comments

""", flavor = WeaveMarkdown.weavemd);

@test md.content[2].formula  == "x = 2"
@test typeof(md.content[3].content[2]) == WeaveMarkdown.Comment
@test md.content[3].content[2].text == " inline "
@test md.content[4].text == "\nMultiple lines\n "

@test WeaveMarkdown.latex(md.content[2]) == "\\[\nx = 2\n\\]\n"
@test WeaveMarkdown.latex(md.content[4]) == "% \n% Multiple lines\n%  \n"

@test WeaveMarkdown.html(md.content[2]) == "
\\[\nx = 2\n\\]
"
@test WeaveMarkdown.html(md.content[4]) == "\n\n"


text:
using RecipesBase, Measures

@recipe function f(
    output::SisPosteriorMcmcOutput,
    S_true::InteractionSequence{Int}
    )

    S_sample = output.S_sample
    γ_sample = output.γ_sample
    layout := (2,1)
    legend --> false
    xguide --> "Index"
    yguide --> ["Distance from True Mode" "γ"]
    size --> (800, 600)
    margin --> 5mm
    y1 = map(x->output.posterior.dist(S_true,x), S_sample)
    y2 = γ_sample
    hcat(y1,y2)
end

@recipe function f(
    outputs::Vector{SisPosteriorMcmcOutput},
    S_true::InteractionSequence{Int}
    )
    d = outputs[1].posterior.dist
    S_samples = [x.S_sample for x in outputs]
    γ_samples = [x.γ_sample for x in outputs]
    layout := (2,1)
    k = length(outputs)
    subplot := [i ≤ k ? 1 : 2 for i in 1:(2*k)]
    legend --> [true false]
    label --> ["Chain $i" for j in 1:1, i in ncycle(1:k, 2)]
    xguide --> "Index"
    yguide --> ["Distance from Reference" "γ"]
    size --> (800, 600)
    margin --> 5mm
    alpha --> 0.5
    y1 = [map(y->d(S_true,y), x) for x in S_samples]
    y2 = γ_samples
    vcat(y1,y2)

end

@recipe function f(output::SisPosteriorModeConditionalMcmcOutput, S_true::Vector{Path{Int}})
    S_sample = output.S_sample
    d = output.posterior.dist
    xguide --> "Index"
    yguide --> "Distance from Truth"
    size --> (800, 300)
    label --> nothing
    margin --> 5mm
    map(x->d(S_true,x), S_sample)
end

@recipe function f(output::SisPosteriorDispersionConditionalMcmcOutput)
    xguide --> "Index"
    yguide --> "γ"
    size --> (800, 300)
    label --> nothing
    margin --> 5mm
    output.γ_sample
end

@recipe function f(
    output::SimPosteriorMcmcOutput,
    S_true::InteractionSequence{Int}
    )

    S_sample = output.S_sample
    γ_sample = output.γ_sample
    layout := (2,1)
    legend --> false
    xguide --> "Index"
    yguide --> ["Distance from True Mode" "γ"]
    size --> (800, 600)
    margin --> 5mm
    y1 = map(x->output.posterior.dist(S_true,x), S_sample)
    y2 = γ_sample
    hcat(y1,y2)
end

@recipe function f(
    outputs::Vector{SimPosteriorMcmcOutput},
    S_true::InteractionSequence{Int}
    )
    d = outputs[1].posterior.dist
    S_samples = [x.S_sample for x in outputs]
    γ_samples = [x.γ_sample for x in outputs]
    layout := (2,1)
    k = length(outputs)
    subplot := [i ≤ k ? 1 : 2 for i in 1:(2*k)]
    legend --> [true false]
    label --> ["Chain $i" for j in 1:1, i in ncycle(1:k, 2)]
    legend --> false
    xguide --> "Index"
    yguide --> ["Distance from True Reference" "γ"]
    size --> (800, 600)
    margin --> 5mm
    alpha --> [i ≤ k ? 1.0 : 0.5 for j in 1:1, i in 1:(2*k)]
    y1 = [map(y->d(S_true,y), x) for x in S_samples]
    y2 = γ_samples
    vcat(y1,y2)
end

@recipe function f(output::SimPosteriorModeConditionalMcmcOutput, S_true::Vector{Path{Int}})
    S_sample = output.S_sample
    d = output.posterior.dist
    xguide --> "Index"
    yguide --> "Distance from Truth"
    size --> (800, 300)
    label --> nothing
    margin --> 5mm
    map(x->d(S_true,x), S_sample)
end

@recipe function f(output::SimPosteriorDispersionConditionalMcmcOutput)
    xguide --> "Index"
    yguide --> "γ"
    size --> (800, 300)
    label --> nothing
    margin --> 5mm
    output.γ_sample
end

@recipe function f(output::T) where {T<:SpfPosteriorMcmcOutput}
    xguide --> "Index"
    yguide --> "γ"
    legend --> false
    size --> (800, 300)
    label := nothing
    output.γ_sample
end

@recipe function f(output::SpfPosteriorMcmcOutput{T}, I_true::Path{T}) where T <:Union{Int, String}
    I_sample = output.I_sample
    xguide --> "Index"
    yguide --> ["Distance from Truth" "γ"]
    legend --> false
    size --> (800, 600)
    label := [nothing nothing]
    layout := (2,1)
    y1 = map(x -> output.dist(I_true, x), I_sample)
    y2 = output.γ_sample
    hcat(y1, y2)
end

@recipe function f(output::SpfPosteriorModeConditionalMcmcOutput{T}, I_true::Path{T}) where T <:Union{Int, String}
    I_sample = output.I_sample
    xguide --> "Index"
    yguide --> "Distance from Truth"
    legend --> false
    size --> (800, 600)
    y1 = map(x -> output.posterior.dist(I_true, x), I_sample)
    y1
end

@recipe function f(output::SpfPosteriorDispersionConditionalMcmcOutput{T}) where T <:Union{Int, String}
    xguide --> "Index"
    yguide --> "γ"
    legend --> false
    size --> (800, 600)
    output.γ_sample
end


text:
export EuclidDistance, PackedEuclidDistance

"""
$(TYPEDEF)

Default linear offset between two scalar variables.
"""
struct EuclidDistance{T <: SamplableBelief} <: AbstractRelativeMinimize
  Z::T
end

EuclidDistance(::UniformScaling=LinearAlgebra.I) = EuclidDistance(Normal())

getDimension(::InstanceType{EuclidDistance{<:SamplableBelief}}) = 1

# before consolidation, see RoME.jl #244
getManifolds(::InstanceType{EuclidDistance{<:SamplableBelief}}) = (:Euclid,)
getDomain(::InstanceType{EuclidDistance{<:SamplableBelief}}) = ContinuousEuclid{1}

getSample(cf::CalcFactor{<:EuclidDistance}, N::Int=1) = (reshape(rand(cf.factor.Z,N),:,N), )

# new and simplified interface for both nonparametric and parametric
function (s::CalcFactor{<:EuclidDistance})(z, x1, x2)
  # v0.21+, should return residual and not have residual parameter
  return z .- norm(x2 .- x1)
end

"""
$(TYPEDEF)
Serialization type for `EuclidDistance` binary factor.
"""
mutable struct PackedEuclidDistance <: PackedInferenceType
  mimeType::String
  Z::String
  # PackedEuclidDistance() = new()
  # PackedEuclidDistance(z::AS) where {AS <: AbstractString} = new(z)
end

function convert(::Type{PackedEuclidDistance}, d::EuclidDistance)
  PackedEuclidDistance( "/application/JuliaLang/PackedSamplableBelief",
                        convert(PackedSamplableBelief, d.Z) )
end

function convert(::Type{<:EuclidDistance}, d::PackedEuclidDistance)
  EuclidDistance(convert(SamplableBelief, d.Z))
end

#


text:
"""
VarCorr

Information from the fitted random-effects variance-covariance matrices.

# Members
* `σρ`: a `NamedTuple` of `NamedTuple`s as returned from `σρs`
* `s`: the estimate of the scale parameter in the distribution of the conditional dist'n of Y

The main purpose of defining this type is to isolate the logic in the show method.
"""
struct VarCorr
    σρ::NamedTuple
    s
end
VarCorr(m::MixedModel) = VarCorr(σρs(m), dispersion_parameter(m) ? dispersion(m) : nothing)

function Base.show(io::IO, vc::VarCorr)
    σρ = vc.σρ
    nmvec = string.([keys(σρ)...])
    cnmvec = string.(foldl(vcat, [keys(sig)...] for sig in getproperty.(values(σρ), :σ)))
    σvec = vcat(collect.(values.(getproperty.(values(σρ), :σ)))...)
    if !isnothing(vc.s)
        push!(σvec, vc.s)
        push!(nmvec, "Residual")
    end
    nmwd = maximum(textwidth.(nmvec)) + 1
    cnmwd = maximum(textwidth.(cnmvec)) + 1
    nρ = maximum(length.(getproperty.(values(σρ), :ρ)))
    varvec = abs2.(σvec)
    showσvec = showoff(σvec, :plain)
    showvarvec = showoff(varvec, :plain)
    varwd = maximum(textwidth.(showvarvec)) + 1
    stdwd = maximum(textwidth.(showσvec)) + 1
    println(io, "Variance components:")
    write(io, " "^(nmwd))
    write(io, cpad("Column", cnmwd))
    write(io, cpad("Variance", varwd))
    write(io, cpad("Std.Dev.", stdwd))
    iszero(nρ) || write(io,"  Corr.")
    println(io)
    ind = 1
    for (i,v) in enumerate(values(vc.σρ))
        write(io, rpad(nmvec[i], nmwd))
        firstrow = true
        k = length(v.σ)   # number of columns in grp factor k
        ρ = v.ρ
        ρind = 0
        for j in 1:k
            !firstrow && write(io, " "^nmwd)
            write(io, rpad(cnmvec[ind], cnmwd))
            write(io, lpad(showvarvec[ind], varwd))
            write(io, lpad(showσvec[ind], stdwd))
            for l in 1:(j - 1)
                ρind += 1
                ρval = ρ[ρind]
                ρval === -0.0 ? write(io, "   .  ") : @printf(io, "%6.2f", ρval)
            end
            println(io)
            firstrow = false
            ind += 1
        end
    end
    if !isnothing(vc.s)
        write(io, rpad(last(nmvec), nmwd))
        write(io, " "^cnmwd)
        write(io, lpad(showvarvec[ind], varwd))
        write(io, lpad(showσvec[ind], stdwd))
    end
    println(io)
end


text:
function print_block(fun::Function,out=stdout)
    io = IOBuffer()
    fun(io)
    data = split(String(take!(io)), "\n")
    if length(data) == 1
        println(out, "[ $(data[1])")
    elseif length(data) > 1
        for (p,dl) in zip(vcat("⎡", repeat(["⎢"], length(data)-2), "⎣"),data)
            println(out, "$(p) $(dl)")
        end
    end
end

function print_bboxes(io::IO, tree::Tree, node::Node, maxdepth::N) where {Tree<:AbstractTree, Node<:AbstractNode, N<:Number}
    println(io, aabb(node))
    (maxdepth -= 1) == 0 && return
    is_inner(node) && print_block(io) do io
        print_bboxes(io, tree, tree[node.left], maxdepth)
        print_bboxes(io, tree, tree[node.right], maxdepth)
    end
end

print_bboxes(bvh::B, maxdepth=Inf) where {B<:BVH} =
    print_bboxes(stdout, bvh.tree, bvh.root,maxdepth)

shadebox(ray::Ray, node::BN, tree::Tree, factor::T) where {T,BN<:Union{BinaryNode{T},BinaryIntNode{T}},Tree<:AbstractTree{T}} =
    intersect(ray, node.bbox) ? factor*(shadebox(ray, tree[node.left], tree, factor) +
                                        shadebox(ray, tree[node.right], tree, factor)) : 1

shadebox(ray::Ray,node::N,tree::Tree, ::T) where {T,N<:Union{LeafNode,LeafIntNode},Tree<:AbstractTree} =
    0 # intersect(ray, node.bbox) ? 1 : 0

function debug_bvh_render(acc::Acc, cam::Camera,
                          factor::T=T(1.3),gain::T=T(1e-4);
                          show_model::Bool=true, kwargs...) where {T,O,B<:BVH{T,O},Acc<:BVHAccelerator{T,O,B}}
    e = one(T)
    z = zero(T)
    black = RGB(z,z,z)
    white = RGB(e,e,e)
    red = RGB(e,z,z)

    tree = acc.bvh.tree
    node = acc.bvh.root

    render(cam, kwargs...) do ray::Ray
        if intersect(ray, node.bbox)
            if show_model && intersect(ray, acc)
                red
            else
                gain*white*shadebox(ray,node,tree,factor)
            end
        else
            black
        end
    end
end


text:
"""
The `SIFTS` module of MIToS allows to obtain the
residue-level mapping between databases stored in the SIFTS XML files.
It makes easy to assign PDB residues to UniProt/Pfam positions.
Given the fact that pairwise alignments can lead to misleading association between residues in both sequences,
SIFTS offers  more reliable association between sequence and structure residue numbers.

**Features**

- Download and parse SIFTS XML files
- Store residue-level mapping in Julia
- Easy generation of `Dict`s between residues numbers

```julia

using MIToS.SIFTS
```
"""
module SIFTS

  using LightXML
  using AutoHashEquals
  using MIToS.Utils

  import Base: hash, ==, call, string, print, write, show, convert, isnull, parse
  import MIToS.Utils: isobject, findobjects, collectobjects, capture, collectcaptures, guess_type

  export DataBase, dbPDBe, dbInterPro, dbUniProt, dbPfam, dbNCBI, dbPDB, dbCATH, dbSCOP,
  SIFTSResidue, downloadsifts, siftsmapping, SIFTSXML,

  # Mitos.Utils
  capture, collectcaptures, isobject, findobjects, collectobjects, Is, Not, In

  include("XMLParser.jl")
  include("ResidueMapping.jl")

end


text:
#=

=#

using DataStructures
using IterTools
function value_iteration(S,A,P,R,Θ)
γ = 0.9
V₀ = DefaultDict{Tuple{Int,Int},Int}(0)
Vk::Array{Any} = [V₀]
k = 1
V₁::Dict{Tuple{Int,Int},Float64} = Dict{Tuple{Int,Int},Float64}()
# update value
while true
    k = k + 1
    push!(Vk,Dict{Tuple{Int,Int},Float64}())
    for s in S
        Vp = []
         for a in A
            #=@show s,a
            println(broadcast(s′ -> R(s,a,s′),[S...]))
            println(broadcast(s′ -> P(s′,s,a),[S...]))=#
            #push!(Vp,sum(broadcast(s′ -> P(s′,s,a)*R(s,a,s′),[S...])))
            push!(Vp, sum(broadcast(s′ -> P(s′,s,a)*(R(s,a,s′)+γ*Vk[k-1][s′]),[S...])))
        end
        #@show Vp
        Vk[k][s] = maximum(Vp)

    end
    abs(sum(abs,values(Vk[k])) - sum(abs,values(Vk[k-1]))) < Θ && break;
end
@show k
Vk[k]
end
ActionToVec = Dict(1 => (0,1),2 => (1,0), 3 => (0,-1), 4 => (-1,0))
function P(s′,s,a)

    Action = [0.1, 0.1, 0.1, 0.1]
    Action[a] = 0.7
    Vec = map(t -> ActionToVec[t],[1,2,3,4])
    Near = map(x -> s.+x,Vec)
    !any([s′] .== Near) && return 0

    Action[findin(Near,[s′])][1]

end
function R(s,a,s′)
    s == (0,0) && return 10
    s == (1,1) && return -10
    0
end

max_stars_repo_name: bollafa/MarkovDecissionProcess

max_issues_repo_name: bollafa/MarkovDecissionProcess

max_forks_repo_name: bollafa/MarkovDecissionProcess


text:
module SemiIsomorphism

using SimpleGraphs, SimpleGraphAlgorithms, ChooseOptimizer, SimpleTools
using Combinatorics, Permutations, JuMP, Gurobi, LinearAlgebra, ProgressMeter

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

export semi_iso, is_semi_iso, semi_auto, auto, semi_mates, semi_auto_1

"""
    semi_iso(G,H)

Find  permutation matrices `P` and `Q` such that `AP==QB`
where `A` and `B` are the
adjacency matrices of `G` and `H` respectively.
"""
function semi_iso(G::SimpleGraph, H::SimpleGraph)
    A = adjacency(G)
    B = adjacency(H)

    if !semi_iso_quick_check(G, H)
        error("These graphs are not semi-isomorphic")
    end

    return semi_iso(A, B)
end

function semi_iso_quick_check(G::SimpleGraph, H::SimpleGraph)
    NV(G) == NV(H) &&
        NE(G) == NE(H) &&
        deg(G) == deg(H) &&
        abs(detx(G)) == abs(detx(H)) &&
        permanent(G) == permanent(H)
end

function semi_iso(A::Matrix, B::Matrix)
    n, c = size(A)
    nn, cc = size(B)
    @assert n == c && nn == cc "Matrices must be square"
    @assert n == nn "Matrices must have the same size"

    m = Model(get_solver())

    @variable(m, P[1:n, 1:n], Bin)
    @variable(m, Q[1:n, 1:n], Bin)

    for i = 1:n
        for k = 1:n
            @constraint(
                m,
                sum(A[i, j] * P[j, k] for j = 1:n) == sum(Q[i, j] * B[j, k] for j = 1:n)
            )
        end
        @constraint(m, sum(P[i, j] for j = 1:n) == 1)
        @constraint(m, sum(P[j, i] for j = 1:n) == 1)

        @constraint(m, sum(Q[i, j] for j = 1:n) == 1)
        @constraint(m, sum(Q[j, i] for j = 1:n) == 1)
    end

    optimize!(m)
    status = Int(termination_status(m))

    if status != 1
        error("These graphs are not semi-isomorphic")
    end

    return Int.(value.(P)), Int.(value.(Q))
end

"""
    is_semi_iso(G,H)

returns `true` is the graphs are semi-isomorphic.
"""
function is_semi_iso(G::SimpleGraph, H::SimpleGraph)::Bool
    try
        semi_iso(G, H)
    catch
        return false
    end
    return true
end

"""
    semi_iso_1(G,H)

One sided semi-isomorphism.
"""
function semi_iso_1(A::Matrix, B::Matrix)
    n, c = size(A)
    nn, cc = size(B)
    @assert n == c && nn == cc "Matrices must be square"
    @assert n == nn "Matrices must have the same size"

    m = Model(get_solver())

    @variable(m, P[1:n, 1:n], Bin)

    for i = 1:n
        for k = 1:n
            @constraint(m, sum(A[i, j] * P[j, k] for j = 1:n) == B[i, k])
        end
        @constraint(m, sum(P[i, j] for j = 1:n) == 1)
        @constraint(m, sum(P[j, i] for j = 1:n) == 1)
    end

    optimize!(m)
    status = Int(termination_status(m))

    if status != 1
        error("These graphs are not simply semi-isomorphic")
    end

    return Int.(value.(P))
end

function semi_iso_1(G::SimpleGraph, H::SimpleGraph)
    A = adjacency(G)
    B = adjacency(H)
    return semi_iso_1(A, B)
end

function semi_auto_1(A::Matrix)
    n, c = size(A)
    @assert n == c "Matrix must be square"

    [Permutation(p) for p in permutations(1:n) if A[p, :] == A]
end

"""
    semi_auto_1(G)

Return all permutations s.t. `AP==A` where `A` is the
adjacency matrix of `G` and `P` is a permutation matrix.
"""
semi_auto_1(G::SimpleGraph) = semi_auto_1(adjacency(G))

function semi_auto(A::Matrix)
    n, c = size(A)
    @assert n == c "Matrix must be square"

    [
        (Permutation(p), Permutation(q)) for p in permutations(1:n) for
        q in permutations(1:n) if A[p, q] == A
    ]
end

"""
    is_adjacency(A::Matrix)
Check if the 0-1 matrix `A` is the adjacency matrix of a simple graph.
"""
function is_adjacency(A::Matrix)::Bool
    A == A' && tr(A) == 0
end

"""
    semi_auto(G)

List out the elements of the semi-automorphism group of G. These are pairs of permutations
such that `P' A Q == A`.
"""
semi_auto(G::SimpleGraph) = semi_auto(adjacency(G))

function auto(A::Matrix)
    n, c = size(A)
    @assert n == c "Matrix must be square"

    [Permutation(p) for p in permutations(1:n) if A[p, p] == A]
end

"""
    auto(G)

List out all the automorphisms of `G`.
"""
auto(G::SimpleGraph) = auto(adjacency(G))

"""
    make_allow(A::Matrix)

Given a matrix, create an "allow" dictionary where the
allowed values for key `k` are the indices of the zeros
in row `k`.
"""
function make_allow(A::Matrix)
    r, c = size(A)
    allow = Dict{Int,Vector{Int}}()

    for i = 1:r
        row = A[i, :]
        allow[i] = findall(t -> t == 0, row)
    end
    return allow
end

"""
    semi_mates(G::SimpleGraph)

Find other graphs that are semi-isomorphic to `G`. Optional named arguments:
* stopper (`false`): Stop at first match different from `G`
* uhash_check (`true`): Skip presumably isomorphic graphs found
* verbose (`true`): Show stuff going on
"""
function semi_mates(
    G::SimpleGraph;
    stopper::Bool = false,
    uhash_check::Bool = true,
    verbose::Bool = true,
)
    A = adjacency(G)
    n = NV(G)

    list = SimpleGraph{Int}[]
    hashes = UInt[]

    push!(list, relabel(G))

    if verbose
        @info "Generating allowable permutations"
    end
    allow = make_allow(A)

    push!(hashes, uhash(G))
    good_perms = Permutations.iter_perms(allow)

    count = 0
    ticker = 1000

    verbose && println("each dot is $ticker permutations")

    for p in good_perms
        count += 1
        if verbose && count % ticker == 0
            print(".")
        end
        B = A[:, p]
        if B != B'
            continue
        end
        if verbose
            println("\n", Permutation(p))
        end
        H = SimpleGraph(B)
        uH = uhash(H)

        if (!uhash_check || uH ∉ hashes) && H ∉ list
            push!(list, H)
            push!(hashes, uH)
            if stopper
                verbose && println("Quick stop")
                return H
            end
        end
    end
    if verbose
        println()
        @info "We considered $count permutations"
    end
    return list

end

function example_maker(A0::Matrix, emb::Bool = false)
    A = dcat(A0, A0)
    B = [0A0 A0; A0 0A0]

    G = SimpleGraph(A)
    H = SimpleGraph(B)

    if emb
        embed(G, :combined)
        embed(H, :combined)
    end
    return G, H
end

function example_maker(G::SimpleGraph, emb::Bool = false)
    A = adjacency(G)
    return example_maker(A, emb)
end

# include("frac_iso.jl")
include("one_side_group.jl")
include("invariants.jl")
include("make_graph.jl")

end # module


text:
export isLegal, lbfgsUpdate, lbfgsHvFunc2, ssbin, solveSubProblem, subHv, polyval, polyinterp

function isLegal(v)
    return sum(isnan.(v))==0 & sum(isinf.(v))==0
end

function lbfgsUpdate(y,s,corrections,debug,old_dirs,old_stps,Hdiag)
    ys = dot(y,s);
    if ys > 1e-10 || size(old_dirs,2)==0
        numCorrections = size(old_dirs,2);
        if numCorrections < corrections
            # Full Update
            old_dirs = [old_dirs s];
            old_stps = [old_stps y];
        else
            # Limited-Memory Update
            old_dirs = [old_dirs[:,2:corrections] s];
            old_stps = [old_stps[:,2:corrections] y];
        end

        # Update scale of initial Hessian approximation
        Hdiag = ys/dot(y,y);
    else
        if debug==3
            @printf("Skipping Update\n");
        end
    end
    return old_dirs, old_stps, Hdiag
end

function lbfgsHvFunc2(v,Hdiag,N,M)
    if cond(M)>(1/(eps(Float32)))
        pr =  Array{Float32}(ssbin(M,500))
        L = spdiagm((pr,),0)
        Hv = v/Hdiag - N*L*((L*M*L)\(L*(transpose(N)*v)))

    else
        Hv = v/Hdiag - N*(M\(transpose(N)*v))
    end

    return Hv
end

function ssbin(A,nmv)
    # Stochastic matrix-free binormalization for symmetric real A.
    # x = ssbin(A,nmv,n)
    #   A is a symmetric real matrix or function handle. If it is a
    #     function handle, then v = A(x) returns A*x.
    #   nmv is the number of matrix-vector products to perform.
    #   [n] is the size of the matrix. It is necessary to specify n
    #     only if A is a function handle.
    #   diag(x) A diag(x) is approximately binormalized.

    # Jan 2010. Algorithm and code by Andrew M. Bradley (ambrad@stanford.edu).
    # Aug 2010. Modified to record and use dp. New omega schedule after running
    #   some tests.
    # Jul 2011. New strategy to deal with reducible matrices: Use the old
    #   iteration in the early iterations; then switch to snbin-like behavior,
    #   which deals properly with oscillation.
    n = size(A,1)
    d = ones(Float32,n,1)
    dp = d
    for k = 1:nmv
      # Approximate matrix-vector product
      u = randn(Float32,n,1)
      s = u./sqrt.(dp)
      y = A*s
      # omega^k
      alpha = (k - 1)/nmv
      omega = (1 - alpha)*1/2 + alpha*1/nmv
      # Iteration
      d = (1-omega)*d/sum(d) + omega*y.^2/sum(y.^2)
      if (k < min(32,floor(nmv/2)))
        # First focus on making d a decent approximation
        dp = d
      else
        # This block makes ssbin behave like snbin except for omega
        tmp = dp
        dp = d
        d = tmp
      end
    end
    return 1f0./(d.*dp).^(1/4)
end

function solveSubProblem(x,g,H,funProj,optTol,progTol,maxIter,testOpt,feasibleInit,x_init)
# Uses SPG to solve for projected quasi-Newton direction
    options = spg_options(verbose=0,
                          optTol=optTol,
                          progTol=progTol,
                          maxIter = maxIter,
                          testOpt = testOpt,
                          feasibleInit = feasibleInit)
    funObj(p) = subHv(p,x,g,H);
    p, f, funEvals, subProjects, hist = minConf_SPG(funObj,x_init,funProj,options);
    return p, subProjects
end

function subHv(p,x,g,HvFunc)
    d = p-x;
    Hd = HvFunc(d);
    f = dot(g,d) + (1f0/2f0)*dot(d,Hd);
    g = g + Hd;
    return f, g
end

function polyval(p,x)
    value = 0
    order = length(p)
    for i=1:length(p)
        value = value + p[i]*x^(order-i)
    end
    return value
end

function polyinterp(points;xminBound=-Inf,xmaxBound=Inf)
# function [minPos] = polyinterp(points,doPlot,xminBound,xmaxBound)
#
#   Minimum of interpolating polynomial based on function and derivative
#   values
#
#   In can also be used for extrapolation if {xmin,xmax} are outside
#   the domain of the points.
#
#   Input:
#       points(pointNum,[x f g])
#       doPlot: set to 1 to plot, default: 0
#       xmin: min value that brackets minimum (default: min of points)
#       xmax: max value that brackets maximum (default: max of points)
#
#   set f or g to sqrt(-1) if they are not known
#   the order of the polynomial is the number of known f and g values minus 1

    nPoints = size(points,1);
    order = length(findin(imag(points[:,2:3]),0))-1;
    # Code for most common case:
    #   - cubic interpolation of 2 points
    #       w/ function and derivative values for both
    #   - no xminBound/xmaxBound

    if nPoints == 2 && order ==3 && isinf.(xmaxBound)
        # Solution in this case (where x2 is the farthest point):
        #    d1 = g1 + g2 - 3*(f1-f2)/(x1-x2);
        #    d2 = sqrt(d1^2 - g1*g2);
        #    minPos = x2 - (x2 - x1)*((g2 + d2 - d1)/(g2 - g1 + 2*d2));
        #    t_new = min(max(minPos,x1),x2);
        minVal, minPos = findmin(points[:,1]);
        notMinPos = -minPos+3;
        d1 = points[minPos,3] + points[notMinPos,3] - 3*(points[minPos,2]-points[notMinPos,2])/(points[minPos,1]-points[notMinPos,1]);
        d2 = sqrt.(d1^2 - points[minPos,3]*points[notMinPos,3]);
        if isreal(d2)
            t = points[notMinPos,1] - (points[notMinPos,1] - points[minPos,1])*((points[notMinPos,3] + d2 - d1)/(points[notMinPos,3] - points[minPos,3] + 2*d2));
            minPos = min(max(t,points[minPos,1]),points[notMinPos,1]);
        else
            minPos = mean(points[:,1]);
        end
        return minPos
    end

    xmin = minimum(points[:,1]);
    xmax = maximum(points[:,1]);

    # Compute Bounds of Interpolation Area
    if isinf.(xminBound)
        xminBound = xmin;
    end
    if isinf.(xmaxBound)
        xmaxBound = xmax;
    end

    # Constraints Based on available Function Values
    A = zeros(Float32,0,order+1);
    b = zeros(Float32,0);
    for i = 1:nPoints
        if imag(points[i,2])==0
            constraint = zeros(Float32,1,order+1);
            for j = order:-1:0
                constraint[order-j+1] = points[i,1]^j;
            end
            A = [A;constraint];
            b = [b;points[i,2]];
        end
    end

    # Constraints based on available Derivatives
    for i = 1:nPoints
        if isreal(points[i,3])
            constraint = zeros(Float32,1,order+1);
            for j = 1:order
                constraint[j] = (order-j+1)*points[i,1]^(order-j);
            end
            A = [A;constraint];
            b = [b;points[i,3]];
        end
    end

    # Find interpolating polynomial
    params = A\b;

    # Compute Critical Points
    dParams = zeros(Float32,order);
    for i = 1:length(params)-1
        dParams[i] = params[i]*(order-i+1);
    end

    if sum(isinf.(dParams)) >0
        cp = copy(transpose([xminBound;xmaxBound;points[:,1]]));
    else
        cp = copy(transpose([xminBound;xmaxBound;points[:,1];-roots(dParams)]));
    end

    # Test Critical Points
    fmin = inf;
    minPos = (xminBound+xmaxBound)/2; # Default to Bisection if no critical points valid
    for xCP = cp
        if imag(xCP)==0 && xCP >= xminBound && xCP <= xmaxBound
            fCP = polyval(params,xCP);
            if imag(fCP)==0 && fCP < fmin
                minPos = real(xCP);
                fmin = real(fCP);
            end
        end
    end

    return minPos, fmin
end


text:
struct Request
    connection::Connection
    method::AbstractString
    uri::AbstractString
    httpversion::AbstractString
    headers::Headers
    resource::AbstractString
    resource_base::AbstractString
    params::Dict{AbstractString,Any}
    query::Dict{AbstractString,Any}
    body::Dict{AbstractString,Any}
    cookies::Dict{AbstractString,Any}
    session::Dict{AbstractString,Any}
    files::Dict{AbstractString,Any}
    data::Dict{AbstractString,Any}
end # Request

function Request(connection::Connection)
    parseduri = URI(connection.uri)

    method = connection.method
    uri = connection.uri
    httpversion = connection.httpversion
    headers = connection.headers
    resource = parseduri.path
    resource_base = ""
    params = Dict{AbstractString, Any}()
    query = Dict{AbstractString, Any}()
    body = Dict{AbstractString, Any}()
    cookies = Dict{AbstractString, Any}()
    session = Dict{AbstractString, Any}()
    files = Dict{AbstractString, Any}()
    data = Dict{AbstractString, Any}()
    Request(connection, method, uri, httpversion, headers, resource,
        resource_base, params, query, body, cookies, session, files, data)
end

function Request(req::Request;
    method=req.method,
    uri=req.uri,
    httpversion=req.httpversion,
    headers=req.headers,
    resource=req.resource,
    resource_base=req.resource_base,
    params=req.params,
    query=req.query,
    body=req.body,
    cookies=req.cookies,
    session=req.session,
    files=req.files,
    data=req.data)

    Request(req.connection, method, uri, httpversion, headers, resource,
        resource_base, params, query, body, cookies, session, files, data)
end

"""
Retrieve a request's header.

The name is case-insensitive

```julia
getheader(req, "Content-Type") == getheader(req, "content-type") # true
```
"""
function getheader(req::Request, name::AbstractString, default=nothing)
    get(req.headers, name, default)
end

"""
Set a header for a request.

This function does not mutate the request or its headers, but instead
returns a new version of the request with the updated header.

If the value is `nothing`, the header will be deleted.

```julia
req2 = setheader(req, "X-Special-Header", "header value")
getheader(req, "X-Special-Header")  # nothing
getheader(req2, "X-Special-Header") # "header value"
```
"""
function setheader(req::Request, name::AbstractString, value::Union{AbstractString,Nothing})
    headers = set(req.headers, name, value)
    Request(req; headers=headers)
end

"""
Append a header for a request.

This function does not mutate the request or its headers, but instead
returns a new version of the request with the updated header.

If the header already exists, the header is appended to the existing header value
with a comma used as a separator.  Otherwise, the header is set to the value.

```julia
req = appendheader(req, "X-Special-Header", "header value 1")
getheader(req, "X-Special-Header") # "header value 1"

req = appendheader(req, "X-Special-Header", "header value 2")
getheader(req, "X-Special-Header") # "header value 1, header value 2"
```
"""
function appendheader(req::Request, name::AbstractString, value::AbstractString)
    headers = append(req.headers, name, value)
    Request(req; headers=headers)
end


text:
#=
Binary (or 0-1) knapsack problem

    Given a knapsack of some capacity C and n objects
       with object i having weight wi and profit pi,
       the goal is to choose some subset of the objects
       that can fit in the knapsack
       (i.e. the sum of their weights is no more than C)
       while maximizing profit.

    This can be formulated as a mixed-integer program as:
        maximize x' p subject to x∈{0,1}  and  w′x <= C
        x is a vector is size n, where x_i
        is 1 if we chose to keep the object in the knapsack,
           0 otherwise.

source: https://nbviewer.jupyter.org/github/JuliaOpt/Convex.jl/blob/master/
                                    examples/binary_knapsack.ipynb
=#

# Data taken from http://people.sc.fsu.edu/~jburkardt/datasets/knapsack_01/knapsack_01.html

w = [23; 31; 29; 44; 53; 38; 63; 85; 89; 82]
C = 165
p =  [92; 57; 49; 68; 60; 43; 67; 84; 87; 72];
n = length(w)

using Convex, GLPKMathProgInterface
x = Variable(n, :Bin)
problem = maximize(dot(p, x), dot(w, x) <= C)
solve!(problem, GLPKSolverMIP())


text:
module NarmaxLasso

export # Auxiliar functions
       generate_random_input,
       simulate,
       # Lasso regressors
       NarmaxRegressors,
       generate_all,
       # Lasso estimation
       narmax_lasso,
       LassoResult,
       # Basis
       Basis,
       Monomial

using RecipesBase
using Combinatorics
if isdir(Pkg.dir("GLMNet"))
    eval(Expr(:import, :GLMNet))
    const glmnet_loaded = true
else
    const glmnet_loaded = false
end

include("util.jl")
include("basis_functions.jl")
include("monomials.jl")
include("narmax_regressors.jl")
include("simulate.jl")
include("pathwise_coordinate_optimization.jl")
include("narmax_lasso.jl")
include("plot_recipes.jl")

end # module


text:
using Lyra
using Electron
using Test

Electron.prep_test_env()

@testset "Lya" begin

end


text:
# Author: Hiroharu Sugawara <hsugawa@gmail.com>
# Copyright: (C) 2020 Hiroharu Sugawara

# Author: Eric P. Hanson
# Copyright: (C) 2018? Eric P. Hanson

# Author: Martin Vuk <martin.vuk@fri.uni-lj.si>
# Copyright: (C) 2016 Martin Vuk

# License: BSD 3-clause

# PandocFIlters.jl

export walk!, toJSONFilter, AST_filter!

export Plain, Para, CodeBlock, RawBlock, BlockQuote
export OrderedList, BulletList, DefinitionList
export Header, HorizontalRule, Table, Div, Null

export Str, Emph, Strong, Strikeout, Superscript, Subscript, SmallCaps
export Quoted, Cite, Code, Space, LineBreak, Math, RawInline
export Link, Image, Note, SoftBreak, Span

"""
    walk!(a, action::Function, format, meta)

Function `walk!` will `walk!` `Pandoc` document abstract source tree (AST) and apply filter function on each element of the document AST.
Returns a modified tree.

  `action` must be a function which takes four arguments, `tag, content, format, meta`,
  and should return

  * `nothing` to leave the element unchanged
  * `[]` to delete the element
  * A Pandoc element to replace the element
  * or a list of Pandoc elements which will be spliced into the list the original object belongs to.
"""
function walk!(x :: Any, action :: Function, format, meta)
    return x
end

function walk!(x :: AbstractArray, action :: Function, format, meta)
  array = []
  w!(z) = walk!(z, action, format, meta)
  for item in x
    if (item isa AbstractDict) && haskey(item,"t")
      res = action(item["t"], get(item, "c", nothing), format, meta)
      if res === nothing
        push!(array, w!(item))
      elseif res isa AbstractArray
        for z in res
          push!(array, w!(z))
        end
      else
        push!(array, w!(res))
      end
    else
      push!(array, w!(item))
    end #if
  end #for
  return array
end

function walk!(dict :: AbstractDict, action :: Function, format, meta)
  # Python version (mutating):
  for k in keys(dict)
    dict[k] = walk!(dict[k], action, format, meta)
  end
  return dict
  # Dict(key => walk!(value,action, format, meta) for (key,value) in dict)
end

"""
    filter!(action::Function)
    filter!(actions::Vector{Function})

Converts an action or a list of actions into a filter that reads a JSON-formatted
pandoc document from stdin, transforms it by walk!ing the tree
with the actions, and returns a new JSON-formatted pandoc document
to stdout.  The argument is a list of functions action(key, value, format, meta),
where key is the type of the pandoc object (e.g. "Str", "Para"),
value is the contents of the object (e.g. a string for "Str",
a list of inline elements for "Para"), format is the target
output format (which will be taken for the first command line
argument if present), and meta is the document's metadata.
If the function returns None, the object to which it applies
will remain unchanged.  If it returns an object, the object will
be replaced.    If it returns a list, the list will be spliced in to
the list to which the target object belongs.    (So, returning an
empty list deletes the object.)
"""
filter(action::Function) = filter([action])

function filter(actions::Vector{T}) where {T<:Function}
  doc = JSON.parse(stdin)
  format = (length(ARGS) <= 0) ? "" : ARGS[1]
  if isa(doc,Dict) && haskey(doc,"meta")
    meta = doc["meta"]
  elseif doc isa AbstractArray  # old API
    meta = doc[1]["test"]
  else
    meta = Dict()
  end

  for action in actions
    doc = walk!(doc, action, format, meta)
  end
  JSON.print(stdout, doc)
end

"""
    AST_filter!(doc, action::Function, format="")
    AST_filter!(doc, actions::Vector{Function}, format="")

Converts an action or a list of actions into a filter that reads a JSON-formatted
pandoc document from stdin, transforms it by walk!ing the tree
with the actions, and returns a new JSON-formatted pandoc document
to stdout.  The argument is a list of functions action(key, value, format, meta),
where key is the type of the pandoc object (e.g. "Str", "Para"),
value is the contents of the object (e.g. a string for "Str",
a list of inline elements for "Para"), format is the target
output format (which will be taken for the first command line
argument if present), and meta is the document's metadata.
If the function returns None, the object to which it applies
will remain unchanged.  If it returns an object, the object will
be replaced.    If it returns a list, the list will be spliced in to
the list to which the target object belongs.    (So, returning an
empty list deletes the object.)
"""
AST_filter!(doc, action::Function; format = "") = AST_filter!(doc, [action], format)

function AST_filter!(doc, actions::Vector{T}; format = "") where {T<:Function}
  if isa(doc,Dict) && haskey(doc, "meta")
    meta = doc["meta"]
  elseif doc isa AbstractArray  # old API
    meta = doc[1]["test"]
  else
    meta = Dict()
  end

  for action in actions
    walk!(doc, action, format, meta)
  end
end

function elt(eltType, numargs)
    function fun(args...)
        lenargs = length(args)
        if lenargs != numargs
            throw(ArgumentError("$eltType expects $numargs arguments, but given $lenargs"))
        end
        if numargs == 0
            return Dict("t" => eltType)
        elseif numargs == 1
            xs = args[1]
        else
            xs = collect(args)
        end
        return Dict("t" => eltType, "c" => xs)
      end
    return fun
end
# Constructors for block elements

const Plain = elt("Plain", 1)
const Para = elt("Para", 1)
const CodeBlock = elt("CodeBlock", 2)
const RawBlock = elt("RawBlock", 2)
const BlockQuote = elt("BlockQuote", 1)
const OrderedList = elt("OrderedList", 2)
const BulletList = elt("BulletList", 1)
const DefinitionList = elt("DefinitionList", 1)
const Header = elt("Header", 3)
const HorizontalRule = elt("HorizontalRule", 0)
const Table = elt("Table", 5)
const Div = elt("Div", 2)
const Null = elt("Null", 0)

# Constructors for inline elements

const Str = elt("Str", 1)
const Emph = elt("Emph", 1)
const Strong = elt("Strong", 1)
const Strikeout = elt("Strikeout", 1)
const Superscript = elt("Superscript", 1)
const Subscript = elt("Subscript", 1)
const SmallCaps = elt("SmallCaps", 1)
const Quoted = elt("Quoted", 2)
const Cite = elt("Cite", 2)
const Code = elt("Code", 2)
const Space = elt("Space", 0)
const LineBreak = elt("LineBreak", 0)
const Math = elt("Math", 2)
const RawInline = elt("RawInline", 2)
const Link = elt("Link", 3)
const Image = elt("Image", 3)
const Note = elt("Note", 1)
const SoftBreak = elt("SoftBreak", 0)
const Span = elt("Span", 2)

# EOF


text:
using DifferentialEquations, Random, Distributions
using Plots, ColorSchemes
using DataFrames, CSV, LaTeXStrings
using Statistics, GLM, Measures
using LightGraphs, SimpleWeightedGraphs#, GraphPlot
using Compose, Cairo, Fontconfig
using LsqFit, Dates
import Gadfly
g = Gadfly # so we can use gadfly too!

pyplot(grid=false)
#  plotly()

fnt = Plots.font("Arial")
default(titlefont=fnt, guidefont=fnt, tickfont=fnt, legendfont=fnt)

epsilon = 1e-8

################################################################################
####  FIG S2B -- write data to CSV
################################################################################

function three_state_model(dc, c, k, t)
    dc[1] = k['a'] * c[2] - k['s'] * c[1]
    dc[2] = k['s'] * c[1] - k['a'] * c[2] - k['i'] * c[2] + k['p'] * c[3]
    dc[3] = k['i'] * c[2] - k['p'] * c[3]

    @assert abs(dc[1]+dc[2]+dc[3])<=epsilon
end

function isoutofdomain(c, k, t)
    for i in 1:3:length(c)
        if abs(1-(c[i] + c[i+1] + c[i+2])) >= epsilon
            println(c)
            println(k)
            println(t)
            return false
        end
    end
    return true
end

function run_recruitment(k_time_pairs, c0=Nothing)
    t = fill(Float64, 0)
    u = fill(Float64[], 0)

    if (c0 == Nothing)
        c0 = [1.0; 0.0; 0.0]
    end

    for pair in k_time_pairs
        k = pair[1]
        times = pair[2]
        prob = ODEProblem(three_state_model, c0, times, k)
        sol = solve(prob, #AutoTsit5(Rosenbrock23());
                    # isoutofdomain=isoutofdomain,
                    saveat=0.01,
                    d_discontinuities=[times[1], times[2]])
        t = [t; sol.t]
        u = [u; sol.u]
        c0 = sol.u[end]
    end
    as = [x[1] for x in u]
    rs = [x[2] for x in u]
    is = [x[3] for x in u]
    data = [as rs is]

    return [t, data]
end

ka = 0.75
ks = 2
ki = 0.5

k_recruit = Dict('a' => ka,
                 's' => ks,
                 'i' => ki,
                 'p' => 0)
k_release = Dict('a' => ka,
                 's' => 0,
                 'i' => 0,
                 'p' => 0)

function get_a(rec, rel)
    ktp = [[k_recruit, (0.0, rec)],
           [k_release, (rec, rec+rel)]]
    sol = run_recruitment(ktp)
    a = sol[2][end, 1]
    return a
end

function get_i(rec, rel)
    ktp = [[k_recruit, (0.0, rec)],
           [k_release, (rec, rec+rel)]]
    sol = run_recruitment(ktp)
    a = sol[2][end, 3]
    return a
end

recs = zeros(0)
rels = zeros(0)
vals = zeros(0)
labs = String[]

for rec in 0.01:0.1:5.01
    for rel in 0.01:0.1:5.01
        #  println(string(rec)*" "*string(rel))
        push!(recs, rec)
        push!(rels, rel)
        push!(labs, "Active")
        push!(vals, get_a(rec, rel))

        push!(recs, rec)
        push!(rels, rel)
        push!(labs, "Irreversibly Silent")
        push!(vals, get_i(rec, rel))
    end
end

data = DataFrame(Dict("Recruitment Duration" => recs,
                      "Release Duration" => rels,
                      "Fraction" => vals,
                      "Population" => labs))

CSV.write("recruitment_graded_response.csv", data)