text:
# ISPC code generation.

@enum VARIABILITY unbound=1 uniform=2 variable=3
const qualifiers = Dict(unbound=>"", uniform=>"uniform", variable=>"variable")

immutable EmitContext
    header::Vector{Any}
    var_types::Dict
    modified::Set
    globals::Set
    cnames::Dict
    sizes::Dict
    declared::Set
    types::Dict
    variability::VARIABILITY
    loc_info::AbstractString
end

function nested_contex(ctx::EmitContext, variability::VARIABILITY)
    EmitContext(ctx.header, ctx.var_types, ctx.modified,
                ctx.globals, ctx.cnames, ctx.sizes,
                ctx.declared, ctx.types, variability, ctx.loc_info)
end

function with_loc_info(ctx::EmitContext, info::AbstractString)
    EmitContext(ctx.header, ctx.var_types, ctx.modified,
                ctx.globals, ctx.cnames, ctx.sizes,
                ctx.declared, ctx.types, ctx.variability, info)
end

const ispc_types = Dict(
    Void => "void",
    Bool => "bool",
    Float32 => "float",
    Float64 => "double",
    Int64 => "int64",
    UInt64 => "unsigned int64",
    Int32 => "int32",
    UInt32 => "unsigned int32",
    Int16 => "int16",
    UInt16 => "unsigned int16",
    Int8 => "int8",
    UInt8 => "unsigned int8",
)

type CType
    cname::ASCIIString
    is_array::Bool
    fields::Dict
end

function declare(ctype, cname, variability=unbound, prefix="")
    qualifier = qualifiers[variability]
    suffix = ctype.is_array? "[]" : ""
    return strip("$qualifier $(ctype.cname) $prefix $cname$suffix")
end

function declare_alias(ctype, cname, variability=unbound, prefix="")
    qualifier = qualifiers[variability]
    prefix = ctype.is_array? "* $qualifier" : ""
    return strip("$qualifier $(ctype.cname) $prefix $cname")
end

function get_ctype(T, ctx::EmitContext)
    ctype = get(ctx.types, T, nothing)
    if ctype == nothing
        ctype = new_ctype(T, ctx)
        ctx.types[T] = ctype
    end
    ctype
end

function new_ctype(T, ctx::EmitContext)
    if T in keys(ispc_types)
        # Some atomic types can be directly translated
        # to ISPC types:
        cname = ispc_types[T]
        return CType(cname, false, Dict())
    elseif length(T.types) > 1
        # Composite type
        return declare_struct(T, ctx)
    elseif issubtype(T, DenseArray)
        # check for DenseArray and not AbstractArray
        # because that would pull in all indexable types
        celtype = get_ctype(eltype(T), ctx)
        return CType(celtype.cname, true, Dict())
    else
        error("Unsupported type $T in $(ctx.loc_info)")
    end
end

function declare_struct(typ, ctx::EmitContext)
    # As a first-level identifier the type also needs a unique name:
    substitute_identifiers([typ.name], ctx.cnames, globals=ctx.globals)
    ctypename = ctx.cnames[typ.name]
    push!(ctx.globals, ctypename)

    field_decls = []
    cfnames = substitute_identifiers(fieldnames(typ))
    for (fname, ftype) in zip(fieldnames(typ), typ.types)
        cfname = cfnames[fname]
        cftype = get_ctype(ftype, ctx)
        fdecl = declare(cftype, cfname)
        push!(field_decls, "$fdecl;")
    end

    decl = """
    struct $(ctypename) {
    $(indent(field_decls))
    };"""
    push!(ctx.header, decl)

    return CType(ctypename, false, cfnames)
end

# A number of definitions to accomodate some quirks in the generated
# code.
const ispc_includes = """
// Use ISPC's multiple dispatch capabilities to deal with the fact
// that Julia uses the same function for bitwise and boolean NOT,
// whereas the ~ operator in ISPC does not work on booleans:
inline bool __not(bool val) {return !val;} // boolean NOT
inline int8 __not(int8 val) {return ~val;} // all others are bitwise
inline int16 __not(int16 val) {return ~val;}
inline int32 __not(int32 val) {return ~val;}
inline int64 __not(int64 val) {return ~val;}
inline unsigned int8 __not(unsigned int8 val) {return ~val;}
inline unsigned int16 __not(unsigned int16 val) {return ~val;}
inline unsigned int32 __not(unsigned int32 val) {return ~val;}
inline unsigned int64 __not(unsigned int64 val) {return ~val;}
"""

function indent(s::AbstractString)
    indent(split(s, '\n'))
end

function indent(lines::AbstractArray)
    join(["    " * string(line) for line in lines], "\n")
end

function emit_ispc(obj::Any, ctx::EmitContext)
    error("Unhandled object $obj ($(typeof(obj))) in $(ctx.loc_info)")
end

function emit_ispc(obj::Void, ctx::EmitContext)
    ""
end

function emit_ispc(num::Real, ctx::EmitContext)
    if num == Inf32
        return "floatbits(0x7F800000)"
    elseif num == -Inf32
        return "floatbits(0xFF800000)"
    end
    # Use hexadecimal float literals so that we
    # can represent the number exactly:
    return @sprintf("%a",num)
end

function emit_ispc(num::Integer, ctx::EmitContext)
    return "$num"
end

function emit_ispc(symbol::SymbolNode, ctx::EmitContext)
    return ctx.cnames[symbol.name]
end

function emit_ispc(symbol::Symbol, ctx::EmitContext)
    return ctx.cnames[symbol]
end

function emit_ispc(symbol::GenSym, ctx::EmitContext)
    return ctx.cnames[symbol]
end

function emit_ispc(node::QuoteNode, ctx::EmitContext)
    return node.value
end

function emit_ispc(top::TopNode, ctx::EmitContext)
    error("Unhandled: $(sexpr(top)) in $(ctx.loc_info)")
end

function emit_ispc(ref::GlobalRef, ctx::EmitContext)
    return emit_ispc(eval(ref), ctx)
end

function emit_ispc(typ::DataType, ctx::EmitContext)
    return get_ctype(typ, ctx).cname
end

function emit_ispc(expr::Expr, ctx::EmitContext)
    emit_ispc(Val{expr.head}, expr.args, with_loc_info(ctx, sexpr(expr)))
end

function emit_ispc(head::Any, args, ctx::EmitContext)
    error("Unhandled expression $head $(sexpr(args)) in $(ctx.loc_info)")
end

function emit_ispc(head::Type{Val{:new}}, args, ctx::EmitContext)
    objtype = args[1]
    values = args[2:end]
    cargs = join([emit_ispc(val, ctx) for val in values], ", ")
    return "{$cargs}"
end

function emit_ispc(head::Type{Val{:block}}, args, ctx::EmitContext)
    body = join([emit_ispc(arg, ctx) for arg in args], ";\n")
    return body
    # Remember: :block does not introduce a nested scope, so it's not
    # equivalent to curly brackets in C.
    # return """{
    # $body
    # }
    # """
end

function emit_ispc(head::Type{Val{:call}}, args, ctx::EmitContext)

    f_expr = args[1]
    f_args = args[2:end]
    f_argtypes = (map(typeof, f_args)...)

    # Catch a few corner cases:
    if f_expr == TopNode(:getfield) && f_argtypes == (GlobalRef, QuoteNode)
        # getfield(GlobalRef, QuoteNode) points to an object that we
        # should be able to resolve here:
        mod = eval(f_args[1])
        symbol = f_args[2].value
        return eval(GlobalRef(mod, symbol), ctx)
    end

    if f_expr == TopNode(:ccall)
        # We don't generally support ccall() in ISPC code,
        # but we can rewrite some to calls to ISPC's stdlib:
        func = eval(f_args[1])
        f_argtypes = eval(f_args[3])
        f_args = f_args[4:4+length(f_argtypes)-1]
        # there is sometimes one extra arg to ccall, so
        # make sure to take only as many as in arg_types.

    else
        # Just resolve the function:
        func = eval(f_expr)
    end

#    try
        return emit_func_call(func, f_args, ctx)
#    catch e
#        if isa(e, KeyError)
#            error("Error while translating $func in $(ctx.loc_info)")
#        end
#        rethrow(e)
#    end
end

function emit_func_call(f, args, ctx::EmitContext)
    if f == Base.arrayref
        return do_arrayref(ctx, args...)
    elseif f == Base.arrayset
        return do_arrayset(ctx, args...)
    elseif f == Base.getfield
        return do_getfield(ctx, args...)
    else
        fgen = basefuncs[f]
        fargs = [emit_ispc(arg, ctx) for arg in args]
        return fgen(fargs...)
    end
end

function do_getfield(ctx::EmitContext, obj, field::QuoteNode)
    cobj = emit_ispc(obj, ctx)
    ctype = get_ctype(ctx.var_types[obj], ctx)
    cfield = ctype.fields[field.value]
    return "$cobj.$cfield"
end

function col_major_index(sizes, indices)
    expr = indices[end]
    for i in length(sizes)-1:-1:1
        expr= :($(indices[i]) + $(sizes[i]) * $expr)
    end
    expr
end

function do_arrayref(ctx::EmitContext, array, I...)
    return indexed_array(ctx, array, I...)
end

function do_arrayset(ctx::EmitContext, array, value, I...)
    item = indexed_array(ctx, array, I...)
    val = emit_ispc(value, ctx)
    return "$item = $val;"
end

function indexed_array(ctx::EmitContext, array, I...)
    arr = emit_ispc(array, ctx)
    sizes = try
        ctx.sizes[array.name]
    catch e
        if isa(e, KeyError)
            error("Can't find array size for $(array.name) in $(ctx.sizes)")
        else
            rethrow(e)
        end
    end
    if length(I) == 1
        iexpr = :($(I[1]) - 1)
        idx = emit_ispc(iexpr, ctx)
    elseif length(I) == length(sizes)
        indices = [:($i-1) for i in I]
        idx = emit_ispc(col_major_index(sizes, indices), ctx)
    else
        error("There must be either one index or as many as dimensions (in $(ctx.loc_info))")
    end
    return "$arr[$idx]"
end

function emit_ispc(head::Type{Val{:return}}, args, ctx::EmitContext)
    if length(args) > 1
        error("Multiple return values are not supported (in $(ctx.loc_info)")
    end
    if args[1] == nothing
        return """return;"""
    else
        carg = emit_ispc(args[1], ctx)
        return """return $carg;"""
    end
end

function emit_ispc(head::Type{Val{:if}}, args, ctx::EmitContext)
    test = emit_ispc(args[1], ctx)
    then = args[2]
    coherent = false
    filter!(then.args) do stmt
        if stmt == Expr(:meta, :ispc, :coherent)
            coherent = true
            return false
        end
        return true
    end

    if_keyword = coherent ? "cif" : "if"
    thenblock = indent(emit_ispc(then, ctx))
    elsebranch = ""
    if length(args) >= 3
        elseblock = indent(emit_ispc(args[3], ctx))
        elsebranch = """
        else {
        $elseblock
        }"""
    end
    return """
    $(if_keyword) ($test) {
    $thenblock
    } $elsebranch"""
end

function emit_ispc(head::Type{Val{:dowhile}}, args, ctx::EmitContext)
    test = emit_ispc(args[1], ctx)
    block = args[2]
    coherent = false
    filter!(block.args) do stmt
        if stmt == Expr(:meta, :ispc, :coherent)
            coherent = true
            return false
        end
        return true
    end

    keyword = coherent ? "cdo" : "do"
    do_block = indent(emit_ispc(block, ctx))
    return """
    $(keyword) {
    $do_block
    } while ($test);"""
end

function emit_ispc(head::Type{Val{:while}}, args, ctx::EmitContext)
    test = emit_ispc(args[1], ctx)
    block = args[2]
    coherent = false
    filter!(block.args) do stmt
        if stmt == Expr(:meta, :ispc, :coherent)
            coherent = true
            return false
        end
        return true
    end

    keyword = coherent ? "cwhile" : "while"
    do_block = indent(emit_ispc(block, ctx))
    return """
    $(keyword)($test) {
    $do_block
    }"""
end

function emit_ispc(head::Type{Val{:continue}}, args, ctx::EmitContext)
    return """continue;"""
end

function emit_ispc(head::Type{Val{:break}}, args, ctx::EmitContext)
    return """break;"""
end

function emit_ispc(head::Type{Val{:foreach}}, args, ctx::EmitContext)
    targets, block = args
    body_ctx = nested_contex(ctx, unbound)
    body = indent(emit_ispc(block, body_ctx))

    active_exprs = (:active, QuoteNode(:active), Expr(:quote, :active))
    if length(targets) == 1 && targets[1].first in active_exprs
        idx = ctx.cnames[targets[1].second]
        return """
        foreach_active($idx) {
        $body
        }"""
    end

    iters = []
    for (object, index) in targets
        idx = ctx.cnames[index]
        if isa(object, SymbolNode)
            # It's an array symbol
            len = ctx.sizes[object.name][1]
            push!(iters, "$idx = 1 ... ($len+1)")
        elseif isa(object, Expr) && object.head == :(:)
            # It's a literal range
            start = emit_ispc(object.args[1], ctx)
            stop = emit_ispc(object.args[2], ctx)
            push!(iters, "$idx = $start ... ($stop+1)")
        else
            dump(object)
            error("unsupported iteration object: $object (in $(ctx.loc_info)")
        end
    end
    return """
    foreach($(iters...)) {
    $body
    }"""
end

function emit_ispc(head::Type{Val{:unmasked}}, args, ctx::EmitContext)
    body = join([indent(emit_ispc(arg, ctx)) for arg in args], "\n")
    return """unmasked {
    $body
    }"""
end

function emit_ispc(head::Type{Val{:(=)}}, args, ctx::EmitContext)
    lhs, rhs = args
    rhs_code = emit_ispc(rhs, ctx)

    if isa(lhs, Symbol) && isa(rhs, SymbolNode)
        # Alias any array sizes:
        try
            ctx.sizes[lhs] = ctx.sizes[rhs.name]
        catch
        end
    end

    cname = ctx.cnames[lhs]
    declared = lhs in ctx.declared
    if !declared
        push!(ctx.declared, lhs)
        ctype = get_ctype(ctx.var_types[lhs], ctx)
        if ctype.is_array
            # all arrays are currently uniform
            # TODO get variability from type
            decl = declare_alias(ctype, cname, uniform)
        else
            decl = declare(ctype, cname, ctx.variability)
        end
        return "$decl = $rhs_code;"
    else
        return "$cname = $rhs_code;"
    end
end

function emit_ispc(newvar::NewvarNode, ctx::EmitContext)
    cname = get(ctx.cnames, newvar.name, nothing)
    if cname == nothing
        println(STDERR, "Warning: encountered bug https://github.com/JuliaLang/julia/issues/12620 -- correctness is not guaranteed.")
        prefix = string(newvar.name)
        for key in keys(ctx.cnames)
            if startswith(string(key), "##$(prefix)#")
                println(STDERR, "Guessing $newvar => $(NewvarNode(key))")
                return emit_ispc(NewvarNode(key), ctx)
            end
        end
        return ""
    end
    ctype = get_ctype(ctx.var_types[newvar.name], ctx)
    decl = declare(ctype, cname, ctx.variability)
    push!(ctx.declared, newvar.name)
    return "$(decl);"
end

import Graphs

function to_ispc(ast, ctx::EmitContext)

    # println(var_types)

    # Get rid of line number nodes:
    ast = strip_lineno(ast)

    # Simplify jump analysis by replacing labels with direct indexing:
    statements = labels_to_lines(ast.args)
    # print_statements(STDOUT, statements)

    graph = to_graph(statements)
    dot = Graphs.to_dot(graph)
    open("graph.dot", "w") do f
        write(f, dot)
    end

    # Recover a structured control flow:
    statements, _ = raise_ast(statements)
    # print_statements(STDOUT, statements)

    ast = Expr(:block, statements...)

    # Make trees from opening and closing :ispc tags:
    ast = meta_to_trees(ast)

    # Collect the indices to foreach statements:
    ast = collect_foreach_indices(ast)
    # println(ast)

    # Emit ISPC code:

    # code = try
        code = emit_ispc(ast, ctx)
    # catch e
    #    println(STDERR, ctx.var_types)
    #    print_statements(STDERR, statements)
    #    rethrow(e)
    # end
    return code
end

function gen_ispc_func(ctx::EmitContext, name, ret_type, body, arg_types, arg_symbols...)
    c_ret_type = get_ctype(ret_type, ctx).cname
    c_args = []
    for (argtype, argname) in zip(arg_types, arg_symbols)
        prefix = argname in ctx.modified ? "&" : ""
        carg = declare(get_ctype(argtype, ctx), ctx.cnames[argname], uniform, prefix)
        push!(c_args, carg)
    end
    c_args_decls = join(c_args, ", ")
    return """
    export $c_ret_type $name($c_args_decls) {
    $body
    }
    """
end

function ispc_codegen!(func::ISPCFunction, kernel_name)
    sizes = Dict()
    argtypes = []
    argnames = []
    for arg in func.arg_names
        if isa(arg, Symbol) || isa(arg, Expr)
            push!(argnames, arg)
            push!(argtypes, func.var_types[arg])
        else
            symbol = arg[1]
            sizes[symbol] = arg[2:end]
            for a in arg
                push!(argnames, a)
                push!(argtypes, func.var_types[a])
            end
        end
    end

    # Replace Julia identifier with C-compatible ones:
    cnames = substitute_identifiers(keys(func.var_types),
                                    globals=func.file.global_names)

    substitute_identifiers([kernel_name], cnames,
                                    globals=func.file.global_names)

    header = []
    declared = Set(argnames)
    ctx = EmitContext(header, func.var_types, func.modified,
                        func.file.global_names, cnames, sizes,
                        declared, Dict(), uniform, "")

    body = indent(to_ispc(func.ast, ctx))

    io = IOBuffer()
    for decl in ctx.header
        write(io, decl)
        write(io, "\n")
    end
    func.ispc_name = cnames[kernel_name]
    func_code = gen_ispc_func(ctx, func.ispc_name, Void, body,
                                   argtypes, argnames...)
    write(io, func_code)
    write(io, "\n")
    func.ispc_code = takebuf_string(io)
end

"""
Converts symbols to unique C-proof identifiers.
"""
function substitute_identifiers(symbols, subst=Dict(); globals=Set())
    idx = 1
    for s in unique(symbols)
        if isa(s, GenSym)
            base = "_gensym$(s.id)"
        else
            base = gen_valid_identifier(string(s))
        end
        name = base
        while name in values(subst) || name in globals
            name = "$(base)_$(idx)"
            idx += 1
        end
        subst[s] = name
    end
    return subst
end

# Converting Julia identifiers to C identifiers:
# We could do it the boring way and generate faceless numbered
# variables, but for clarity and aesthetics we'll be extra nice
# and even try to spell out the names of unicode characters.
# Sorry about the wasted microseconds. Maybe add a debug switch.

latex_symbols = Base.REPLCompletions.latex_symbols
const unicode_names = Dict([reverse(p) for p in latex_symbols])
const disallowed = r"[^a-zA-Z0-9_]"
# this does not catch all cases but is enough when the input is
# already a valid Julia identifier (eg. no initial digit etc.)

function gen_valid_identifier(s::ASCIIString)
    replace(s, disallowed, "_")
end

function gen_valid_identifier(s::UTF8String)
    # Replace all non-ascii characters with either the
    # corresponding latex command, or an underscore:
    s2 = join(map(collect(s)) do c
        isascii(c) ? c : get(unicode_names, "$c", "_")
    end, "")
    gen_valid_identifier(ascii(s2))
end


text:
using JSON

file = open("C:/Users/rmsms/OneDrive/바탕 화면/response.json")

ITEM = JSON.parse(file)

Active = []
for item in ITEM
    push!(Active, item["Active"])
end

plot(Active)

ext: jl

max_stars_repo_name: rmsmsgood/aespa

max_issues_repo_name: rmsmsgood/aespa

max_forks_repo_name: rmsmsgood/aespa


text:
module ChampionV3 # LOLTools

using ..LOLTools: AbstractDTO, lol_api_server, http_action, call_api, nothing_in_event
using HTTP

struct ChampionInfo <: AbstractDTO
    freeChampionIdsForNewPlayers
    freeChampionIds
    maxNewPlayerLevel
    ChampionInfo(freeChampionIdsForNewPlayers, freeChampionIds, maxNewPlayerLevel) = new(freeChampionIdsForNewPlayers, freeChampionIds, maxNewPlayerLevel)
end

"""
         champion_rotations(api_key::String,
                            platform::String ;
                            endpoint::HTTP.URI = lol_api_server(platform),
                            action::Function = http_action,
                            event::Function = nothing_in_event)::ChampionInfo
"""
function champion_rotations(api_key::String,
                            platform::String ;
                            endpoint::HTTP.URI = lol_api_server(platform),
                            action::Function = http_action,
                            event::Function = nothing_in_event)::ChampionInfo
    call_api(ChampionInfo, api_key, action, endpoint, "/lol/platform/v3/champion-rotations", event, champion_rotations)
end

end # module LOLTools.ChampionV3


text:
@doc raw"""
    guess_density(basis)

Build a superposition of atomic densities (SAD) guess density.

We take for the guess density a gaussian centered around the atom, of
length specified by `atom_decay_length`, normalized to get the right number of electrons
```math
\hat{ρ}(G) = Z \exp\left(-(2π \text{length} |G|)^2\right)
```
"""
guess_density(basis::PlaneWaveBasis) = guess_density(basis, basis.model.atoms)
@timing function guess_density(basis::PlaneWaveBasis{T}, atoms) where {T}
    model = basis.model
    ρ = zeros(complex(T), basis.fft_size)
    # If no atoms, start with a zero initial guess
    isempty(atoms) && return from_fourier(basis, ρ)
    # fill ρ with the (unnormalized) Fourier transform, ie ∫ e^{-iGx} ρ(x) dx
    for (spec, positions) in atoms
        n_el_val = n_elec_valence(spec)
        decay_length::T = atom_decay_length(spec)
        for (iG, G) in enumerate(G_vectors(basis))
            Gsq = sum(abs2, model.recip_lattice * G)
            form_factor::T = n_el_val * exp(-Gsq * decay_length^2)
            for r in positions
                ρ[iG] += form_factor * cis(-2T(π) * dot(G, r))
            end
        end
    end

    # projection in the normalized plane wave basis
    from_fourier(basis, ρ / sqrt(model.unit_cell_volume))
end

@doc raw"""
Get the lengthscale of the valence density for an atom with `n_elec_core` core
and `n_elec_valence` valence electrons.
```
"""
function atom_decay_length(n_elec_core, n_elec_valence)
    # Adapted from ABINIT/src/32_util/m_atomdata.F90,
    # from which also the data has been taken.

    n_elec_valence = round(Int, n_elec_valence)
    if n_elec_valence == 0
        return 0.0
    end

    data = if n_elec_core < 0.5
        # Bare ions: Adjusted on 1H and 2He only
        [0.6, 0.4, 0.3, 0.25, 0.2]
    elseif n_elec_core < 2.5
        # 1s2 core: Adjusted on 3Li, 6C, 7N, and 8O
        [1.8, 1.4, 1.0, 0.7, 0.6, 0.5, 0.4, 0.35, 0.3]
    elseif n_elec_core < 10.5
        # Ne core (1s2 2s2 2p6): Adjusted on 11na, 13al, 14si and 17cl
        [2.0, 1.6, 1.25, 1.1, 1.0, 0.9, 0.8, 0.7 , 0.7, 0.7, 0.6]
    elseif n_elec_core < 12.5
        # Mg core (1s2 2s2 2p6 3s2): Adjusted on 19k, and on n_elec_core==10
        [1.9, 1.5, 1.15, 1.0, 0.9, 0.8, 0.7, 0.6 , 0.6, 0.6, 0.5]
    elseif n_elec_core < 18.5
        # Ar core (Ne + 3s2 3p6): Adjusted on 20ca, 25mn and 30zn
        [2.0, 1.8, 1.5, 1.2, 1.0, 0.9, 0.85, 0.8, 0.75, 0.7, 0.65, 0.65, 0.6]
    elseif n_elec_core < 28.5
        # Full 3rd shell core (Ar + 3d10): Adjusted on 31ga, 34se and 38sr
        [1.5, 1.25, 1.15, 1.05, 1.00, 0.95, 0.95, 0.9, 0.9, 0.85, 0.85, 0.80,
         0.8 , 0.75, 0.7]
    elseif n_elec_core < 36.5
        # Krypton core (Ar + 3d10 4s2 4p6): Adjusted on 39y, 42mo and 48cd
        [2.0, 2.00, 1.60, 1.40, 1.25, 1.10, 1.00, 0.95, 0.90, 0.85, 0.80, 0.75, 0.7]
    else
        # For the remaining elements, consider a function of n_elec_valence only
        [2.0 , 2.00, 1.55, 1.25, 1.15, 1.10, 1.05, 1.0 , 0.95, 0.9, 0.85, 0.85, 0.8]
    end
    data[min(n_elec_valence, length(data))]
end
atom_decay_length(sp::Element) = atom_decay_length(n_elec_core(sp), n_elec_valence(sp))


text:
Base.@kwdef struct FluxAccumulator{B} <: AbstractPhysicsComponent
    bctype::B
end

function calc_component!(source, accumulator::FluxAccumulator, state, aux, physics)

    # save flux in atmos source for integration in time using atmos timestepping
    if boundary_mask(physics.parameters, aux.x, aux.y, aux.z)
        E, H = calc_ocean_sfc_fluxes(physics, accumulator.bctype, state, aux)
        source.F_ρe_accum = (E + H) # latent + sensible heat fluxes [W/m^2]
    end
end

"""
    calculate_land_sfc_fluxes(model::DryAtmosModel, state, aux)
- calculate furface fluxes using the bulk gradient diffusion theory
"""

function calc_ocean_sfc_fluxes(physics, bctype::BulkFormulaTemperature, state⁻, aux⁻)
    # Apply bulks laws using the tangential velocity as energy flux
    ρ = state⁻.ρ
    ρu = state⁻.ρu
    ρq = state⁻.ρq
    eos = physics.eos
    parameters = physics.parameters

    # vertical unit vector
    n̂ = aux⁻.∇Φ / parameters.g

    # obtain surface fields from bcs
    ϕ = lat(aux⁻.x, aux⁻.y, aux⁻.z)
    Cₕ = bctype.drag_coef_temperature(parameters, ϕ)
    Cₑ = bctype.drag_coef_moisture(parameters, ϕ)
    T_sfc = bctype.surface_temperature(parameters, ϕ)
    LH_v0 = parameters.LH_v0

    # magnitude of tangential velocity (usually called speed)
    u = ρu / ρ
    speed_tangential = norm((I - n̂ ⊗ n̂) * u)

    # sensible heat flux
    cp = calc_heat_capacity_at_constant_pressure(eos, state⁻, parameters)
    T = calc_air_temperature(eos, state⁻, aux⁻, parameters)
    H = ρ * Cₕ * speed_tangential * cp * (T_sfc - T)

    # latent heat flux
    q = ρq / ρ
    q_tot_sfc = calc_saturation_specific_humidity(ρ, T_sfc, parameters)
    E = ρ * Cₑ * speed_tangential * LH_v0 * (q_tot_sfc - q)

    return E, H
end

function calc_ocean_sfc_fluxes(physics, bctype::CoupledPrimarySlabOceanBC, state⁻, aux⁻; MO_params = nothing) # should pass in the coupler state (also move to coupler), so can access states of both models derectly -e.g. callback?

    # Apply bulks laws using the tangential velocity as energy flux
    ρ = state⁻.ρ
    ρu = state⁻.ρu
    ρq = state⁻.ρq
    eos = physics.eos
    parameters = physics.parameters

    # vertical unit vector
    n̂ = aux⁻.∇Φ / parameters.g

    # obtain surface fields from bcs
    Cₕ = parameters.Cₗ
    Cₑ = parameters.Cₑ
    LH_v0 = parameters.LH_v0
    T_sfc = aux⁻.T_sfc

    # magnitude of tangential velocity (usually called speed)
    u = ρu / ρ
    speed_tangential = norm((I - n̂ ⊗ n̂) * u)

    # sensible heat flux
    cp = calc_heat_capacity_at_constant_pressure(eos, state⁻, parameters)
    T = calc_air_temperature(eos, state⁻, aux⁻, parameters)
    H = ρ * Cₕ * speed_tangential * cp * (T_sfc - T)

    # latent heat flux
    q = ρq / ρ
    q_tot_sfc = calc_saturation_specific_humidity(ρ, T_sfc, parameters)
    E = ρ * Cₑ * speed_tangential * LH_v0 * (q_tot_sfc - q)

    E = isnan(E) ? zero(E) : E # set the fluxes of the non-sfc layers to 0 (otherwise NaNs if moisture) TODO: if calc_component imports BL, then could apply the FluxAccumulator on the boundary only and this line won't be needed

    return E, H
end

max_stars_repo_name: CliMA/CouplerMachine

max_issues_repo_name: CliMA/CouplerMachine


text:
using Documenter, RRRMC

CIbuild = get(ENV, "CI", nothing) == "true"

makedocs(
    modules  = [RRRMC],
    format   = Documenter.HTML(prettyurls = CIbuild),
    sitename = "RRRMC.jl",
    pages    = Any[
        "Home" => "index.md",
        "Sampling algorithms" => "algorithms.md",
        "Graph types" => "graph-types.md",
        "Built-in graphs" => "graphs-builtin.md",
        "Graphs interface" => "interface.md"
       ]
    )

deploydocs(
    repo   = "github.com/carlobaldassi/RRRMC.jl.git",
    target = "build",
    deps = nothing,
    make = nothing,
)


text:
# Update W
function update_W!(s::StateFS, c::ConstantsFS, d::DataFS)
  wr = s.W_star .* s.r
  s.theta.W .= wr ./ sum(wr, dims=2)
end

# Update W star
function update_W_star!(s::StateFS, c::ConstantsFS, d::DataFS, tuners::TunersFS)
  for i in 1:d.data.I
    for k in 1:c.constants.K
      update_W_star!(i, k, s, c, d, tuners)
      # NOTE: Make sure to always update W immediately after updating r or
      #       W_star!
      update_W!(s, c, d)
    end
  end

  return
end

function update_W_star!(i::Int, k::Int, s::StateFS, c::ConstantsFS, d::DataFS,
                        tuners::TunersFS)
  if s.r[i, k] == 0
    # Sample from prior
    s.W_star[i, k] = rand(c.W_star_prior)
  else  # s.r[i, k] == 1
    logprior(w::Float64) = logpdf(c.W_star_prior, w)

    function loglike(w::Float64)::Float64
      log_numer = sum(s.theta.lam[i] .== k) * log(w)
      w_other = s.W_star[i, 1:end .!= k]
      r_other = s.r[i, 1:end .!= k]
      log_denom = d.data.N[i] * log(w + sum(w_other .* r_other))
      return log_numer - log_denom
    end

    log_prob(w::Float64)::Float64 = logprior(w) + loglike(w)

    # Metropolis step to update W*_{ik}
    s.W_star[i, k] = MCMC.metLogAdaptive(s.W_star[i, k], log_prob, tuners.W_star[i, k])
  end
end


text:
function car_dynamics(x::Array{T, 1} where {T <: Real},
	                  u::Array{T, 1} where {T <: Real})
  lr = 1.5
  lf = 1.8

  beta = atan(lr / (lr + lf) * tan(u[2]))
  dx1 = x[4] * cos(x[3] + beta)
  dx2 = x[4] * sin(x[3] + beta)
  dx3 = x[4] / lr * sin(beta)
  dx4 = u[1]
  return [dx1, dx2, dx3, dx4]
end

function car_dynamics_overt(range_dict::Dict{Symbol, Array{T, 1}} where {T <: Real},
                            N_OVERT::Int,
							t_idx::Union{Int, Nothing}=nothing)
	lr = 1.5
    lf = 1.8

	if isnothing(t_idx)
		v1 = :(atan($(lr / (lr + lf)) * tan(u2)))
	    v1_oA = overapprox(v1, range_dict; N=N_OVERT)
	    range_dict = merge(range_dict, Dict(v1_oA.output => v1_oA.output_range))

	    v2 = :($(v1_oA.output) + x3)
	    v2_oA = overapprox(v2, range_dict; N=N_OVERT)
	    range_dict = merge(range_dict, Dict(v2_oA.output => v2_oA.output_range))

	    v3 = :(x4 * cos($(v2_oA.output)))
	    v3_oA = overapprox(v3, range_dict; N=N_OVERT)
	    range_dict = merge(range_dict, Dict(v3_oA.output=> v3_oA.output_range))

	    v4 = :(x4 * sin($(v2_oA.output)))
	    v4_oA = overapprox(v4, range_dict; N=N_OVERT)
	    range_dict = merge(range_dict, Dict(v4_oA.output=> v4_oA.output_range))

	    v5 = :($(1 / lr) * x4 * sin($(v1_oA.output)))
	    v5_oA = overapprox(v5, range_dict; N=N_OVERT)
	    range_dict = merge(range_dict, Dict(v5_oA.output=> v5_oA.output_range))
	else
    	v1 = "atan($(lr / (lr + lf)) * tan(u2_$t_idx))"
		v1 = Meta.parse(v1)
	    v1_oA = overapprox(v1, range_dict; N=N_OVERT)
	    range_dict = merge(range_dict, Dict(v1_oA.output => v1_oA.output_range))

	    v2 = "$(v1_oA.output) + x3_$t_idx"
		v2 = Meta.parse(v2)
	    v2_oA = overapprox(v2, range_dict; N=N_OVERT)
	    range_dict = merge(range_dict, Dict(v2_oA.output => v2_oA.output_range))

	    v3 = "x4_$t_idx * cos($(v2_oA.output))"
		v3 = Meta.parse(v3)
	    v3_oA = overapprox(v3, range_dict; N=N_OVERT)
	    range_dict = merge(range_dict, Dict(v3_oA.output=> v3_oA.output_range))

	    v4 = "x4_$t_idx * sin($(v2_oA.output))"
		v4 = Meta.parse(v4)
	    v4_oA = overapprox(v4, range_dict; N=N_OVERT)
	    range_dict = merge(range_dict, Dict(v4_oA.output=> v4_oA.output_range))

	    v5 = "$(1 / lr) * x4_$t_idx * sin($(v1_oA.output))"
		v5 = Meta.parse(v5)
	    v5_oA = overapprox(v5, range_dict; N=N_OVERT)
	    range_dict = merge(range_dict, Dict(v5_oA.output=> v5_oA.output_range))
	end
    oA_out = add_overapproximate([v1_oA, v2_oA, v3_oA, v4_oA, v5_oA])
    return oA_out, [v3_oA.output, v4_oA.output, v5_oA.output]
end

function car_update_rule(input_vars::Array{Symbol, 1},
                         control_vars::Array{Symbol, 1},
						 overt_output_vars::Array{Symbol, 1})
    integration_map = Dict(input_vars[1] => overt_output_vars[1],
                           input_vars[2] => overt_output_vars[2],
                           input_vars[3] => overt_output_vars[3],
                           input_vars[4] => control_vars[1])
    return integration_map
end

car_input_vars = [:x1, :x2, :x3, :x4]
car_control_vars = [:u1, :u2]

Car = OvertProblem(
	car_dynamics,
	car_dynamics_overt,
	car_update_rule,
	car_input_vars,
	car_control_vars
)


text:
using ExaPF
using ExaPF.Parse
using ExaPF.PowerSystem

# file locations
# datafile = "GO-Data/datasets/Trial_3_Real-Time/Network_30R-025/scenario_1/case.raw"
# datafile = "GO-Data/datasets/Trial_3_Real-Time/Network_13R-015/scenario_11/case.raw"
# datafile = "test/case14.raw"
# npartition: Number of partitions for the additive Schwarz preconditioner
function pf(datafile, npartition=2, solver="default")
    # read data
    println(solver)

    pf = PowerSystem.PowerNetwork(datafile)
    x = ExaPF.PowerSystem.get_x(pf)
    u = ExaPF.PowerSystem.get_u(pf)
    p = ExaPF.PowerSystem.get_p(pf)
    return solve(pf, x, u, p, npartition, solver);
end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

"""
    GitAnnotated(repo::GitRepo, commit_id::GitHash)
    GitAnnotated(repo::GitRepo, ref::GitReference)
    GitAnnotated(repo::GitRepo, fh::FetchHead)
    GitAnnotated(repo::GitRepo, comittish::AbstractString)

An annotated git commit carries with it information about how it was looked up and
why, so that rebase or merge operations have more information about the context of
the commit. Conflict files contain information about the source/target branches in
the merge which are conflicting, for instance. An annotated commit can refer to the
tip of a remote branch, for instance when a [`FetchHead`](@ref) is passed, or to a
branch head described using `GitReference`.
"""
function GitAnnotated(repo::GitRepo, commit_id::GitHash)
    ann_ptr_ptr = Ref{Ptr{Void}}(C_NULL)
    @check ccall((:git_annotated_commit_lookup, :libgit2), Cint,
                  (Ptr{Ptr{Void}}, Ptr{Void}, Ptr{GitHash}),
                   ann_ptr_ptr, repo.ptr, Ref(commit_id))
    return GitAnnotated(repo, ann_ptr_ptr[])
end

function GitAnnotated(repo::GitRepo, ref::GitReference)
    ann_ref_ref = Ref{Ptr{Void}}(C_NULL)
    @check ccall((:git_annotated_commit_from_ref, :libgit2), Cint,
                  (Ptr{Ptr{Void}}, Ptr{Void}, Ptr{Void}),
                   ann_ref_ref, repo.ptr, ref.ptr)
    return GitAnnotated(repo, ann_ref_ref[])
end

function GitAnnotated(repo::GitRepo, fh::FetchHead)
    ann_ref_ref = Ref{Ptr{Void}}(C_NULL)
    @check ccall((:git_annotated_commit_from_fetchhead, :libgit2), Cint,
                  (Ptr{Ptr{Void}}, Ptr{Void}, Cstring, Cstring, Ptr{GitHash}),
                   ann_ref_ref, repo.ptr, fh.name, fh.url, Ref(fh.oid))
    return GitAnnotated(repo, ann_ref_ref[])
end

function GitAnnotated(repo::GitRepo, comittish::AbstractString)
    obj = GitObject(repo, comittish)
    cmt = peel(GitCommit, obj)
    return GitAnnotated(repo, GitHash(cmt))
end

function GitHash(ann::GitAnnotated)
    Base.@gc_preserve ann begin
        oid = unsafe_load(ccall((:git_annotated_commit_id, :libgit2), Ptr{GitHash}, (Ptr{Void},), ann.ptr))
    end
    return oid
end

"""
    merge_analysis(repo::GitRepo, anns::Vector{GitAnnotated}) -> analysis, preference

Run analysis on the branches pointed to by the annotated branch tips `anns` and
determine under what circumstances they can be merged. For instance, if `anns[1]`
is simply an ancestor of `ann[2]`, then `merge_analysis` will report that a
fast-forward merge is possible.

Return two outputs, `analysis` and `preference`. `analysis` has several possible values:
    * `MERGE_ANALYSIS_NONE`: it is not possible to merge the elements of `anns`.
    * `MERGE_ANALYSIS_NORMAL`: a regular merge, when HEAD and the commits that the
      user wishes to merge have all diverged from a common ancestor. In this case the
      changes have to be resolved and conflicts may occur.
    * `MERGE_ANALYSIS_UP_TO_DATE`: all the input commits the user wishes to merge can
      be reached from HEAD, so no merge needs to be performed.
    * `MERGE_ANALYSIS_FASTFORWARD`: the input commit is a descendant of HEAD and so no
      merge needs to be performed - instead, the user can simply checkout the
      input commit(s).
    * `MERGE_ANALYSIS_UNBORN`: the HEAD of the repository refers to a commit which does not
      exist. It is not possible to merge, but it may be possible to checkout the input
      commits.
`preference` also has several possible values:
    * `MERGE_PREFERENCE_NONE`: the user has no preference.
    * `MERGE_PREFERENCE_NO_FASTFORWARD`: do not allow any fast-forward merges.
    * `MERGE_PREFERENCE_FASTFORWARD_ONLY`: allow only fast-forward merges and no
      other type (which may introduce conflicts).
`preference` can be controlled through the repository or global git configuration.
"""
function merge_analysis(repo::GitRepo, anns::Vector{GitAnnotated})
    analysis = Ref{Cint}(0)
    preference = Ref{Cint}(0)
    anns_ref = Ref(map(a->a.ptr, anns))
    anns_size = Csize_t(length(anns))
    @check ccall((:git_merge_analysis, :libgit2), Cint,
                  (Ptr{Cint}, Ptr{Cint}, Ptr{Void}, Ptr{Ptr{Void}}, Csize_t),
                   analysis, preference, repo.ptr, anns_ref, anns_size)
    return analysis[], preference[]
end

"""
    ffmerge!(repo::GitRepo, ann::GitAnnotated)

Fastforward merge changes into current HEAD. This is only possible if the commit
referred to by `ann` is descended from the current HEAD (e.g. if pulling changes
from a remote branch which is simply ahead of the local branch tip).
"""
function ffmerge!(repo::GitRepo, ann::GitAnnotated)
    cmt = GitCommit(repo, GitHash(ann))

    checkout_tree(repo, cmt)
    with(head(repo)) do head_ref
        cmt_oid = GitHash(cmt)
        msg = "libgit2.merge: fastforward $(string(cmt_oid)) into $(name(head_ref))"
        new_head_ref = if reftype(head_ref) == Consts.REF_OID
            target!(head_ref, cmt_oid, msg=msg)
        else
            GitReference(repo, cmt_oid, fullname(head_ref), msg=msg)
        end
        close(new_head_ref)
    end
    return true
end

# Merge changes into current head
"""
    merge!(repo::GitRepo, anns::Vector{GitAnnotated}; kwargs...) -> Bool

Merge changes from the annotated commits (captured as [`GitAnnotated`](@ref) objects)
`anns` into the HEAD of the repository `repo`. The keyword arguments are:
  * `merge_opts::MergeOptions = MergeOptions()`: options for how to perform the merge,
    including whether fastforwarding is allowed. See [`MergeOptions`](@ref) for more
    information.
  * `checkout_opts::CheckoutOptions = CheckoutOptions()`: options for how to perform
    the checkout. See [`CheckoutOptions`](@ref) for more information.

`anns` may refer to remote or local branch heads. Return `true` if the merge is
successful, otherwise return `false` (for instance, if no merge is possible
because the branches have no common ancestor).

# Examples
```julia
upst_ann = LibGit2.GitAnnotated(repo, "branch/a")

# merge the branch in
LibGit2.merge!(repo, [upst_ann])
```
"""
function merge!(repo::GitRepo, anns::Vector{GitAnnotated};
                merge_opts::MergeOptions = MergeOptions(),
                checkout_opts::CheckoutOptions = CheckoutOptions())
    anns_size = Csize_t(length(anns))
    @check ccall((:git_merge, :libgit2), Cint,
                  (Ptr{Void}, Ptr{Ptr{Void}}, Csize_t,
                   Ptr{MergeOptions}, Ptr{CheckoutOptions}),
                   repo.ptr, map(x->x.ptr, anns), anns_size,
                   Ref(merge_opts), Ref(checkout_opts))
    info("Review and commit merged changes.")
    return true
end

# Internal implementation of merge.
# Returns `true` if merge was successful, otherwise `false`
"""
    merge!(repo::GitRepo, anns::Vector{GitAnnotated}, fastforward::Bool; kwargs...) -> Bool

Merge changes from the annotated commits (captured as [`GitAnnotated`](@ref) objects)
`anns` into the HEAD of the repository `repo`. If `fastforward` is `true`, *only* a
fastforward merge is allowed. In this case, if conflicts occur, the merge will fail.
Otherwise, if `fastforward` is `false`, the merge may produce a conflict file which
the user will need to resolve.

The keyword arguments are:
  * `merge_opts::MergeOptions = MergeOptions()`: options for how to perform the merge,
    including whether fastforwarding is allowed. See [`MergeOptions`](@ref) for more
    information.
  * `checkout_opts::CheckoutOptions = CheckoutOptions()`: options for how to perform
    the checkout. See [`CheckoutOptions`](@ref) for more information.

`anns` may refer to remote or local branch heads. Return `true` if the merge is
successful, otherwise return `false` (for instance, if no merge is possible
because the branches have no common ancestor).

# Examples
```julia
upst_ann_1 = LibGit2.GitAnnotated(repo, "branch/a")

# merge the branch in, fastforward
LibGit2.merge!(repo, [upst_ann_1], true)

# merge conflicts!
upst_ann_2 = LibGit2.GitAnnotated(repo, "branch/b")
# merge the branch in, try to fastforward
LibGit2.merge!(repo, [upst_ann_2], true) # will return false
LibGit2.merge!(repo, [upst_ann_2], false) # will return true
```
"""
function merge!(repo::GitRepo, anns::Vector{GitAnnotated}, fastforward::Bool;
                merge_opts::MergeOptions = MergeOptions(),
                checkout_opts::CheckoutOptions = CheckoutOptions())
    ma, mp = merge_analysis(repo, anns)
    if isset(ma, Cint(Consts.MERGE_ANALYSIS_UP_TO_DATE))
        return true # no merge - everything is up to date
    end

    ffpref = if fastforward
        Consts.MERGE_PREFERENCE_FASTFORWARD_ONLY
    elseif isset(mp, Cint(Consts.MERGE_PREFERENCE_NONE))
        Consts.MERGE_PREFERENCE_NONE
    elseif isset(mp, Cint(Consts.MERGE_PREFERENCE_NO_FASTFORWARD))
        Consts.MERGE_PREFERENCE_NO_FASTFORWARD
    elseif isset(mp, Cint(Consts.MERGE_PREFERENCE_FASTFORWARD_ONLY))
        Consts.MERGE_PREFERENCE_FASTFORWARD_ONLY
    else
        throw(ArgumentError("unknown merge preference: $(mp)."))
    end

    merge_result = if ffpref == Consts.MERGE_PREFERENCE_NONE
        if isset(ma, Cint(Consts.MERGE_ANALYSIS_FASTFORWARD))
            if length(anns) > 1
                warn("Unable to perform Fast-Forward merge with mith multiple merge heads.")
                false
            else
                ffmerge!(repo, anns[1])
            end
        elseif isset(ma, Cint(Consts.MERGE_ANALYSIS_NORMAL))
            merge!(repo, anns,
                   merge_opts=merge_opts,
                   checkout_opts=checkout_opts)
        end
    elseif ffpref == Consts.MERGE_PREFERENCE_FASTFORWARD_ONLY
        if isset(ma, Cint(Consts.MERGE_ANALYSIS_FASTFORWARD))
            if length(anns) > 1
                warn("Unable to perform Fast-Forward merge with mith multiple merge heads.")
                false
            else
                ffmerge!(repo, anns[1])
            end
        else
            warn("Cannot perform fast-forward merge.")
            false
        end
    elseif ffpref == Consts.MERGE_PREFERENCE_NO_FASTFORWARD
        if isset(ma, Cint(Consts.MERGE_ANALYSIS_NORMAL))
            merge!(repo, anns,
                   merge_opts=merge_opts,
                   checkout_opts=checkout_opts)
        end
    else
        throw(ArgumentError("unknown merge analysis result: $(ma)"))
    end
    return merge_result
end

"""
    merge_base(repo::GitRepo, one::AbstractString, two::AbstractString) -> GitHash

Find a merge base (a common ancestor) between the commits `one` and `two`.
`one` and `two` may both be in string form. Return the `GitHash` of the merge base.
"""
function merge_base(repo::GitRepo, one::AbstractString, two::AbstractString)
    oid1_ptr = Ref(GitHash(one))
    oid2_ptr = Ref(GitHash(two))
    moid_ptr = Ref(GitHash())
    moid = try
        @check ccall((:git_merge_base, :libgit2), Cint,
                (Ptr{GitHash}, Ptr{Void}, Ptr{GitHash}, Ptr{GitHash}),
                moid_ptr, repo.ptr, oid1_ptr, oid2_ptr)
        moid_ptr[]
    catch e
        #warn("Pkg:",path(repo),"=>",e.msg)
        GitHash()
    end
    return moid
end

max_stars_repo_name: ihnorton/julia

max_issues_repo_name: JuliaLangEs/julieta

max_forks_repo_name: JuliaLangEs/julieta


text:
using LinearOptimizationSolver
using Plots
using Polyhedra

A = [2 1; 1 2]
b = [4, 4]
c = [4, 3]
input = create(A, b, c, tol=1E-10, solver = 1,verbose=false,standard=false)
output = solve(input)

x = output.X[1,:]
y = output.X[2,:]
n = length(x)
u = zeros(n)
v = zeros(n)
for i = 1:n-1
    u[i] = x[i+1] - x[i]
    v[i] = y[i+1] - y[i]
end

h = hrep([A; -1 0; 0 -1], [b; 0; 0])
p = polyhedron(h)

image = plot(p)
quiver!(x, y, quiver = (u, v))
quiver!([x[end]], [y[end]], quiver = ([c[1]]/4, [c[2]]/4))
plot!(x,y; seriestype = :scatter, title = "Problema da produção",legend = false)
savefig(image, "examples\\prod_path")

gap = output.gap

image = plot(log10.(gap),title = "Problema da produção",legend = false)
savefig(image, "examples\\prod_conv")


text:
export LexicographicPathOrder, LPO

abstract type AbstractOrder end

"""
Generates lexicographic path order based on given signature order.

# Examples
```jldoctest
julia> >ₗₚₒ = LPO((:i, 1), (:f, 2), (:e, 0));
```
"""
struct LexicographicPathOrder <: AbstractOrder
    order::Vector{Tuple{Symbol,Int}}
end
LexicographicPathOrder(order::Tuple...) = LexicographicPathOrder(collect(order))
Base.in(t, lpo::LexicographicPathOrder) = t ∈ lpo.order
const LPO = LexicographicPathOrder

(>ₗₚₒ::LPO)(s, t) = false
(>ₗₚₒ::LPO)(s::Term, t::Variable) = occursin(t, s)

function (>ₗₚₒ::LPO)(s::Fn, t::Fn)
    F, G = (s.name, length(s)), (t.name, length(t))
    F ∈ >ₗₚₒ || throw(ArgumentError("$F is not contained in order"))
    G ∈ >ₗₚₒ || throw(ArgumentError("$G is not contained in order"))

    s == t && return false

    any(s) do sᵢ
        sᵢ >ₗₚₒ t || sᵢ == t
    end && return true

    all(t) do tⱼ
        s >ₗₚₒ tⱼ
    end || return false

    findfirst(==(F), >ₗₚₒ.order) < findfirst(==(G), >ₗₚₒ.order) && return true

    if F == G
        for (sᵢ, tᵢ) ∈ zip(s, t)
            sᵢ == tᵢ && continue
            sᵢ >ₗₚₒ tᵢ && return true
        end
    end

    return false
end


text:
using ForwardDiff
using DiffRules

import ForwardDiff: Dual, Partials, single_seed, partials, Chunk, Tag, seed!, value,
              gradient, valtype, extract_gradient!, derivative,extract_derivative,
              checktag, vector_mode_gradient, vector_mode_jacobian,
              chunk_mode_gradient

const AMBIGUOUS_TYPES = (AbstractFloat, Irrational, Integer, Rational,
                         Real, RoundingMode, ComplexF64)

include("dual.jl")
include("config.jl")
include("derivative.jl")
include("gradient.jl")
include("jacobian.jl")
include("apiutils.jl")
include("diffrules.jl")


text:
using AbstractGPs
using BlockDiagonals
using ChainRulesCore
using FillArrays
using FiniteDifferences
using LinearAlgebra
using KernelFunctions
using Random
using StaticArrays
using StructArrays
using TemporalGPs
using Test
using Zygote

using FiniteDifferences: rand_tangent
using AbstractGPs: var
using TemporalGPs: AbstractLGSSM, _filter, NoContext
using Zygote: Context, _pullback

include("test_util.jl")

@testset "TemporalGPs.jl" begin

    println("util:")
    @testset "util" begin
        include(joinpath("util", "harmonise.jl"))
        include(joinpath("util", "scan.jl"))
        include(joinpath("util", "zygote_friendly_map.jl"))
        include(joinpath("util", "zygote_rules.jl"))
        include(joinpath("util", "gaussian.jl"))
        include(joinpath("util", "mul.jl"))
        include(joinpath("util", "regular_data.jl"))
    end

    println("models:")
    include(joinpath("models", "model_test_utils.jl"))
    include(joinpath("models", "test_model_test_utils.jl"))
    @testset "models" begin
        include(joinpath("models", "linear_gaussian_conditionals.jl"))
        include(joinpath("models", "gauss_markov_model.jl"))
        include(joinpath("models", "lgssm.jl"))
        include(joinpath("models", "missings.jl"))
    end

    println("gp:")
    @testset "gp" begin
        include(joinpath("gp", "lti_sde.jl"))
        include(joinpath("gp", "posterior_lti_sde.jl"))
    end

    println("space_time:")
    @testset "space_time" begin
        include(joinpath("space_time", "rectilinear_grid.jl"))
        include(joinpath("space_time", "regular_in_time.jl"))
        include(joinpath("space_time", "separable_kernel.jl"))
        include(joinpath("space_time", "to_gauss_markov.jl"))
        include(joinpath("space_time", "pseudo_point.jl"))
    end
end


text:
asset(url...) = @compat readstring(joinpath(dirname(@__FILE__), "..", "assets", "bower_components", url...))

#Vega Scaffold: https://github.com/vega/vega/wiki/Runtime
function writehtml(io::IO, v::VegaLiteVis; title="Vega.jl Visualization")
  d3          = asset("d3","d3.min.js")
  vega        = asset("vega", "vega.js")
  vegalite    = asset("vega-lite", "vega-lite.js")
  vegaembed   = asset("vega-embed", "vega-embed.js")

  divid = "vg" * randstring(3)

  println(io,
  """

      <title>$title</title>
      <meta charset="UTF-8">

  """)
end

function show(io::IO, v::VegaLiteVis)

    if displayable("text/html")
        v
    else
        # create a temporary file
        tmppath = string(tempname(), ".vegalite.html")
        io = open(tmppath, "w")
        writehtml(io, v)
        close(io)

        # Open the browser
        @osx_only     run(`open $tmppath`)
        @windows_only run(`cmd /c start $tmppath`)
        @linux_only   run(`xdg-open $tmppath`)
    end

    return
end

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

# asset(url...) = readall(Pkg.dir("Vega", "assets", "bower_components", url...))

# function writemime(io::IO, ::MIME"text/html", v::VegaVisualization)
#     divid = "vg" * randstring(3)
#     script_contents = scriptstr(v, divid)
#     display("text/html", """
#
#

#
#
#     """)
# end
#
# import Patchwork: Elem

# function patchwork_repr(v::VegaLiteVis)
#     divid = "vg" * randstring(3)
#     script_contents = scriptstr(v, divid)
#     Elem(:div, [
#         Elem(:div, "") & Dict(:id=>divid, :style=>Dict("min-height"=>"$(v.height + 110)px")),
#         Elem(:script, script_contents) & Dict(:type=>"text/javascript")
#     ])
# end

# function writemime(io::IO, m::MIME"text/html", v::VegaLiteVis)
#
#
#   fh = """
#

#
#
#
#
#
#     <dom-module id="vega-lite-plot">
#       <template>
#

#             <content></content>
#

#       </template>
#
#     </dom-module>
#
#     <vega-lite-plot json="$(v.vis)" svg ></vega-lite-plot>
#
#

#   """
#
#   # writemime(io, m, fh)
#   write(io, fh)
# end

# function jslibpath(url...)
#   libpath = joinpath(dirname(@__FILE__), "..", "assets", "bower_components", url...)
#   replace(libpath, "\\", "/")  # for windows...
# end
#
#
# function writemime(io::IO, m::MIME"text/html", v::VegaLiteVis)
#   divid = "vl" * randstring(3) # generated id for this plot

#   script = """
#     require.config({
#       paths: {
#         d3:        "https://d3js.org/d3.v3.min",
#         vega:      "https://vega.github.io/vega/vega.min",
#         vegalite:  "https://vega.github.io/vega-lite/vega-lite",
#         vegaembed: "https://vega.github.io/vega-editor/vendor/vega-embed"
#       }
#     });
#
#     require(["d3"], function(d3){
#
#       window.d3 = d3;
#
#       require(["vega", "vegalite", "vegaembed"],
#              function(vg, vgl, vge){
#
#         // window.vg = vg;
#
#         vg.parse.spec({
#   "description": "A simple bar chart with embedded data.",
#   "data": {
#     "values": [
#       {"a": "A","b": 28}, {"a": "B","b": 55}, {"a": "C","b": 43},
#       {"a": "D","b": 91}, {"a": "E","b": 81}, {"a": "F","b": 53},
#       {"a": "G","b": 19}, {"a": "H","b": 87}, {"a": "I","b": 52}
#     ]
#   },
#   "mark": "bar",
#   "encoding": {
#     "x": {"field": "a", "type": "ordinal"},
#     "y": {"field": "b", "type": "quantitative"}
#   }
# }, function(chart) { chart({el:"#$divid"}).update(); });
#
#         document.getElementById("$divid").insertAdjacentHTML('beforeend', '
Save yourgh ! as PNG')
#
#
#       });
#     });
#   """
#
#   spec = """
#   {
#     "width": 400,
#     "height": 200,
#     "padding": {"top": 10, "left": 30, "bottom": 30, "right": 10},
#     "data": [
#       {
#         "name": "table",
#         "values": [
#           {"x": 1,  "y": 28}, {"x": 2,  "y": 55},
#           {"x": 3,  "y": 43}, {"x": 4,  "y": 91},
#           {"x": 5,  "y": 81}, {"x": 6,  "y": 53},
#           {"x": 7,  "y": 19}, {"x": 8,  "y": 87},
#           {"x": 9,  "y": 52}, {"x": 10, "y": 48},
#           {"x": 11, "y": 24}, {"x": 12, "y": 49},
#           {"x": 13, "y": 87}, {"x": 14, "y": 66},
#           {"x": 15, "y": 17}, {"x": 16, "y": 27},
#           {"x": 17, "y": 68}, {"x": 18, "y": 16},
#           {"x": 19, "y": 49}, {"x": 20, "y": 15}
#         ]
#       }
#     ],
#     "scales": [
#       {
#         "name": "x",
#         "type": "ordinal",
#         "range": "width",
#         "domain": {"data": "table", "field": "x"}
#       },
#       {
#         "name": "y",
#         "type": "linear",
#         "range": "height",
#         "domain": {"data": "table", "field": "y"},
#         "nice": true
#       }
#     ],
#     "axes": [
#       {"type": "x", "scale": "x"},
#       {"type": "y", "scale": "y"}
#     ],
#     "marks": [
#       {
#         "type": "rect",
#         "from": {"data": "table"},
#         "properties": {
#           "enter": {
#             "x": {"scale": "x", "field": "x"},
#             "width": {"scale": "x", "band": true, "offset": -1},
#             "y": {"scale": "y", "field": "y"},
#             "y2": {"scale": "y", "value": 0}
#           },
#           "update": {
#             "fill": {"value": "steelblue"}
#           },
#           "hover": {
#             "fill": {"value": "red"}
#           }
#         }
#       }
#     ]
#   }
#   """
#
#   spec2 = """
#   {
#     "description": "A simple bar chart with embedded data.",
#     "data": {
#       "values": [
#         {"a": "A","b": 28}, {"a": "B","b": 55}, {"a": "C","b": 43},
#         {"a": "D","b": 91}, {"a": "E","b": 81}, {"a": "F","b": 53},
#         {"a": "G","b": 19}, {"a": "H","b": 87}, {"a": "I","b": 52}
#       ]
#     },
#     "mark": "bar",
#     "encoding": {
#       "x": {"field": "a", "type": "ordinal"},
#       "y": {"field": "b", "type": "quantitative"}
#     }
#   }
#   """
#
#   script = """
#         require.config({
#           paths: {
#             d3: "https://vega.github.io/vega-editor/vendor/d3.min",
#             vega: "https://vega.github.io/vega/vega.min",
#             cloud: "https://vega.github.io/vega-editor/vendor/d3.layout.cloud",
#             topojson: "https://vega.github.io/vega-editor/vendor/topojson",
#             vegalite: "https://vega.github.io/vega-lite/vega-lite",
#             embed: "https://vega.github.io/vega-editor/vendor/vega-embed"
#           }
#         });
#
#         require(["d3", "topojson", "cloud"], function(d3, topojson, cloud){
#
#           window.d3 = d3;
#           window.topojson = topojson;
#           window.d3.layout.cloud = cloud;
#           console.log("d3")
#
#
#             require(["vega"], function(vg) {
#
#               window.vg = vg
#               console.log("vg")
#
#               require(["vegalite"], function(vgl) {
#
#                 console.log("vegalite")
#
#                 var embedSpec = {
#                   renderer: "svg",
#                   actions: false,
#                   spec: $spec2
#                 };
#
#                 var vgSpec = vgl.compile(embedSpec.spec).spec;
#                 vg.parse.spec(vgSpec, function(chart) { chart({el:\"#$divid\"}).update(); });
#
#               }); //vegaembed require end
#             }); //vega require end
#
#         }); //d3 require end
#     """
    # vg.embed("#$divid", embedSpec, function(error, result) { });
    # vg.parse.spec($spec, function(chart) { chart({el:\"#$divid\"}).update(); });

    # window.setTimeout(function() {
    # var pnglink = document.getElementById(\"$divid\").getElementsByTagName(\"canvas\")[0].toDataURL(\"image/png\")
    # document.getElementById(\"$divid\").insertAdjacentHTML('beforeend', '
Save as PNG')
    #
    # }, 20);

  # var vgSpec = vgl.compile(embedSpec.spec).spec;
  # var embedSpec = {
  # mode: "vega-lite",
  # renderer: "svg",
  # actions: false,
  # spec: $(JSON.json(v.vis))
  # };
  #
  # vg.embed("#$divid", embedSpec, function(error, result) { result.view });

  # patch = Elem(:div, [
  #             Elem(:div, "") & Dict(:id=>divid, :style=>Dict("min-height"=>"300px")),
  #             # Elem(:div, "") & Dict(:id=>divid),
  #             Elem(:script, script) & Dict(:type=>"text/javascript")
  #         ])
  #
  # return writemime(io, m, patch)

  # fh = """
  #

  #
  #
  #
  #

  # // d3: "$(jslibpath("d3","d3.min.js"))",
  # // vega: "$(jslibpath("vega", "vega.js"))",
  # // vegalite: "$(jslibpath("vega-lite", "vega-lite.js"))",
  # // vegaembed: "$(jslibpath("vega-embed", "vega-embed.js"))"
#
#   fh = """
#

#

#
#     <meta charset="utf-8">
#
#
#
#

#   """
#   # FIXME : understand why vega-embed can't be loaded
#   # vg.embed("#$divid", embedSpec, function(error, result) {
#
#   println(fh)
#   # writemime(io, m, fh)
#   write(io, fh)
# end
#
#
# export writemime


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "Gloo"
version = v"0.0.20200317"

# Collection of sources required to complete build
sources = [
    GitSource("https://github.com/facebookincubator/gloo.git", "113bde13035594cafdca247be953610b53026553"),
    DirectorySource("./bundled"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/gloo
atomic_patch -p1 ../patches/mingw32.patch
mkdir build
cd build
if [[ $target != *w64-mingw32* ]]; then
    cmake_extra_args="-DUSE_LIBUV=ON"
fi
cmake \
    -DCMAKE_INSTALL_PREFIX=$prefix \
    -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} \
    -DCMAKE_BUILD_TYPE=Release \
    -DBUILD_SHARED_LIBS=ON \
    $cmake_extra_args \
    ..
cmake --build . -- -j $nproc
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()
filter!(p -> nbits(p) == 64, platforms) # Gloo can only be built on 64-bit systems
platforms = expand_cxxstring_abis(platforms)

# The products that we will ensure are always built
products = [
    LibraryProduct("libgloo", :libgloo),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    BuildDependency("LibUV_jll"),
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies;
    preferred_gcc_version=v"5",
    julia_compat="1.6")

max_stars_repo_name: cmhamel/Yggdrasil

max_issues_repo_name: cmhamel/Yggdrasil

max_forks_repo_name: cmhamel/Yggdrasil


text:
# #########################################################################
# ##############################   GaN  ###################################
# #########################################################################

# function force(r::Vector{<:Real}, p::GaN, Type1::Symbol, Type2::Symbol)

#     if (Type1 == :Ga) && (Type2 == :Ga) # Ga-Ga interaction
#         return force(r, p.c) + force(r, p.lj_Ga_Ga)
#     elseif (Type1 == :N) && (Type2 == :N) # N-N interaction
#         return force(r, p.c) + force(r, p.lj_N_N)
#     else
#         return force(r, p.c) + force(r, p.bm_Ga_N)
#     end
# end

# function force(r::Vector{Atom}, p::MixedPotential)
#     n = length(r)
#     # f = Array{Float64}(undef, 3, n)
#     f = [zeros(3) for j = 1:n]
#     for i = 1:(n-1)
#         ri = r[i].Position
#         for j = (i+1):n
#             rj = r[j].Position
#             rtemp = ri - rj
#             f[i] +=  force(rtemp, p, r[i].Type, r[j].Type)
#             f[j] -= f[i]
#         end
#     end
#     return f
# end

# function force(c::Configuration, p::MixedPotential)
#     return force(c.Atoms, p)
# end

# function force(r::Vector{Configuration}, p::MixedPotential)
#     n = length(r)

#     f = [[zeros(3) for j = 1:length(r[i].Atoms)] for i = 1:n]

#     for i = 1:n
#         f[i] = force(r[i], p)
#     end
#     return f
# end


text:
function GradientFunction(arg0::MultivariateDifferentiableFunction)
    return GradientFunction((MultivariateDifferentiableFunction,), arg0)
end

function value(obj::GradientFunction, arg0::Vector{jdouble})
    return jcall(obj, "value", Vector{jdouble}, (Vector{jdouble},), arg0)
end


text:
InputFolder = './Images/CellsCyto/';
OutputFolder = './Results/Images/CellsCyto/';
Lbl = 1;

@iA = '*_C0000*.tif';		% Image filter
@iS = '*_C0001*.tif';		% Image seeds filter

@fxgs_lSeededPropagate [iA, iS] > [L, P];
params.Thr = 0;
params.Power = 0.5;
params.AnalyzeCC = 0;
/endf

%show iA > iS;
/show iA > L;
/keep L > tif;

max_stars_repo_name: SebastienTs/LOBSTER

max_issues_repo_name: SebastienTs/LOBSTER

max_forks_repo_name: SebastienTs/LOBSTER


text:
######################################################
#
# Example testing all features of multiscene
#=
 Features are  up to now
 - Arrangement of several subscenes with  in a given scene grid
 - ',' key switches between gallery view and focus view showing one  subscene
 - space key can block/unblock main task
 - handle interactive variables associated with  hotkeys via keyboard interaction,
   separately for each subscene. Rationale: sliders would eat up too much screen real estate
 - Store context dict in the attributes of LScene
=#

import Pkg
Pkg.activate(mktempdir())
Pkg.add(["ExtendableGrids","GLMakie"]

using ExtendableGrids.MultiScene
using GLMakie

setmakie!(GLMakie)

# Scene with two internal variables, shown only once
function subscene1!(ax)
    N=10
    s=1.0
    makedata(s,N)=s*rand(N)
    makestatus(s,N)="s=$(round(s,digits=2)) N=$(N)"
    data=Node(makedata(s,N))
    status=Node(makestatus(s,N))

    scene=lines(lift(a->a,data), color = :red, linewidth = 4)

    scene_interaction(scene,[:n,:s]) do delta,key
        # key: which key  is "switched on" (of those passed in the array)
        # delta: increment/decrement value
        if key==:n
            N=max(2,N+delta)
        elseif key==:s
            s=max(0.1,s+0.01*Float64(delta))
        end
        data[]=makedata(s,N)
        status[]=makestatus(s,N)
        update!(scene)
    end
    add_scene!(ax,scene,title="lines",status=status)
end

# Scene with two internal variables, shown multiple times
# as a consequence, data and nodes are stored in a context dict
# which is stored in ax.attributes. Not sure if this is a good idea...
function subscene2!(ax,k,l)
    ctx=subscenecontext(ax)

    N=100
    makedata(k,l)=[sinpi(2*k*i/N)*sinpi(2*l*j/N) for i=1:N, j=1:N]
    makestatus(k,l)="k=$(round(k,digits=2)) l=$(round(k,digits=2))"

    if isempty(ctx)
        data=Node(makedata(k,l))
        status=Node(makestatus(k,l))
        scene=heatmap(lift(a->a,data))
        ctx[:data]=data
        ctx[:status]=status
        ctx[:k]=k
        ctx[:l]=l

        scene_interaction(scene,[:k,:l]) do delta,key
            ctx[key]=max(0.1,ctx[key]+0.1*delta)
            status[]=makestatus(ctx[:k], ctx[:l])
            data[]=makedata(ctx[:k], ctx[:l])
        end
        add_scene!(ax,scene,title="heatmap2d",status=status)
        ctx
    else
        ctx[:k]=k
        ctx[:l]=l
        ctx[:status][]=makestatus(ctx[:k], ctx[:l])
        ctx[:data][]=makedata(ctx[:k], ctx[:l])
        yieldwait()
        ctx
    end
end

# simple subscene with one interior variable
# so no need to care about switching
function subscene3!(ax)
    N=50
    data=Node(rand(N,3))
    status=Node("N=$(N)")
    scene=scatter(lift(a->a,data), color = rand(RGBf0))
    scene_interaction(scene) do delta,key
        N=max(4,N+delta)
        status[]="N=$(N)"
        data[]=rand(N,3)
    end
    add_scene!(ax,scene,title="scatter3d",status=status)
end

# simple subscene with one interior variable
function subscene4!(ax)
    iso=1.7

    r = LinRange(-1, 1, 100)
    cube = [(x.^2 + y.^2 + z.^2) for x = r, y = r, z = r]
    cubedata=cube .* (cube .> 1.4)
    data=Node(iso)
    status=Node("iso=$(round(iso,digits=2))")
    scene=volume(cubedata, algorithm = :iso, isorange = 0.05, isovalue = lift(a->a,data))

    scene_interaction(scene) do delta,key
        iso=min(2.5,max(1.5,iso+0.05*delta))
        data[]=iso
        status[]="iso=$(round(iso,digits=2))"
    end
    add_scene!(ax,scene,title="cube3d",status=status)
end

function test_multiscene()

    parent,subscenes=multiscene(layout=(2,2))

    subscene1!(subscenes[1,1])
    subscene2!(subscenes[1,2],1.9,3.5)
    subscene3!(subscenes[2,1])
    subscene4!(subscenes[2,2])

    display(parent)

    # this loop runs "forever" and can be temporarily
    # stopped by the space key
    k=2.0
    l=2.0
    dir=1.0
    while true
        for i=1:1000
            k+=dir*0.01
            l+=dir*0.01
            subscene2!(subscenes[1,2],k,l)
        end
        dir=-dir
    end
end


text:
# Test the wigner 3-j symbol calculations
@testset "wigner3j" begin

    # Meta-parameters
    j_max = 300     # Range of j-values to test
    N = 1000        # Number of randomized tests to perform

    # Compute Wigner matrices
    wigner_A, wigner_B = compute_wigner_values((2j_max + 1), j_max + 1, 2j_max + 1)

    println("Verifying Wigner Symbols...")

    # Outputs
    phase_function_results_A = Array{Float64,1}()
    phase_function_results_B = Array{Float64,1}()
    wigner_symbols_results_A = Array{Float64,1}()
    wigner_symbols_results_B = Array{Float64,1}()

    # Set counter to 1
    count = 0

    # Loop until we see N non-zeros Wigner values
    while count < N

        # Random inputs
        m = rand(1:j_max)
        n = rand(1:j_max)
        l = rand(1:j_max)

        # Result from PhaseFunction module
        push!(phase_function_results_A, wigner_A[m, n, l])
        push!(phase_function_results_B, wigner_B[m, n, l])

        # Result from external Wigner Symbols package
        # If there's a domain error, replace with 0.0.
        # (The 3j symbol *should* be zero outside the domain)
        try
            push!(wigner_symbols_results_A, Float64(wigner3j(m, n, l-1, -1, 1, 0)))
        catch
            push!(wigner_symbols_results_A, 0.0)
        end

        try
            push!(wigner_symbols_results_B, Float64(wigner3j(m, n, l-1, -1, -1, 2)))
        catch
            push!(wigner_symbols_results_B, 0.0)
        end

        # If a discrepancy ever pops up, print the discrepancy so it can be reproduced:
        if (phase_function_results_A[end] ≉ wigner_symbols_results_A[end] ||
            phase_function_results_B[end] ≉ wigner_symbols_results_B[end])
            println("Error with: ", (m, n, l))
            println("PhaseFunction output: ", (phase_function_results_A[end], phase_function_results_B[end]))
            println("WignerSymbols output: ", (wigner_symbols_results_A[end], wigner_symbols_results_B[end]))
        end

        # Only increment the counter if non-zero
        (phase_function_results_A[end] > 0) && (count = count + 1)

    end

    # Compare the result with WignerSymbols package values
    @test phase_function_results_A ≈ wigner_symbols_results_A
    @test phase_function_results_B ≈ wigner_symbols_results_B
end

# Test the Aerosol Optics calculations (both NAI2 and Siewert)
@testset "aerosol_optics" begin

    println("Testing NAI2 and PCW equivalence...")

    # STEP 1: Create the Aerosol

    # Aerosol particle distribution and properties
    μ  = 0.3                # Log mean radius
    σ  = 6.82               # Log stddev of radius
    r_max = 30.0            # Maximum radius
    nquad_radius = 2500     # Number of quadrature points for integrating of size dist.
    nᵣ = 1.3                # Real part of refractive index
    nᵢ = 0.001              # Imag part of refractive index
    size_distribution = LogNormal(log(μ), log(σ))

    # Create the aerosol
    aero = make_univariate_aerosol(size_distribution, r_max, nquad_radius, nᵣ, nᵢ)

    # STEP 2: Create the Mie Calculations model

    λ = 0.55   # Incident wavelength
    polarization_type = Stokes_IQUV()
    truncation_type = δBGE(10, 10)
    model_NAI2 = make_mie_model(NAI2(), aero, λ, polarization_type, truncation_type)

    ###
    ### NOTE: Temporarily removing PCW tests because they are too heavy to run on Travis
    ###

    # Get saved wigner matrices
    # ftp = FTP("ftp://fluo.gps.caltech.edu/XYZT_hitran/")
    # println("Downloading full Wigner values...")
    # download(ftp, "wigner_values.jld", "/tmp/wigner_values.jld");

    # println("Loading full Wigner values...")
    # wigner_A, wigner_B = load_wigner_values("/home/rjeyaram/RadiativeTransfer/src/Scattering/Mie/wigner_values.jld")
    # model_PCW = make_mie_model(PCW(), aero, λ, polarization_type, truncation_type, wigner_A, wigner_B)

    # STEP 3: Perform the Mie Calculations and compare the results

    aerosol_optics_NAI2 = compute_aerosol_optical_properties(model_NAI2);
    # aerosol_optics_PCW = compute_aerosol_optical_properties(model_PCW);

    # Load truth values computed from PCW
    @load "helper/PCW_AerosolOptics.jld" aerosol_optics_PCW

    @test aerosol_optics_NAI2 ≈ aerosol_optics_PCW

    println("Testing aerosol_optical autodiff...")

    # Test whether autodiff works
    compute_aerosol_optical_properties(model_NAI2, autodiff=true);

end


text:
using StackOverflow

"""
    searcherror(intitle::String)

Adding the ability to search for an error through StackOverflow.jl while in the Julia Terminal.
"""
function searcherror(intitle::String)
    data = replace(intitle, " " => "%20")
    r = HTTP.request("GET", "https://api.stackexchange.com/2.2/search?order=desc&sort=activity&intitle=$(data)&site=stackoverflow")
    json_obj = StackOverflow.convert_HTTP_Response_To_JSON(r)
    for (k, v) in json_obj
        if occursin("items", k)
            for item in v
                for (key, value) in item
                    if key == "title" || key == "link"
                        println("$(key) $(value)")
                    end
                end
                println("\n")
            end
        end
    end
end


text:
struct HookeStiffness
    stiffness::Any
    lamecoeffs
    function HookeStiffness(stiffness,lamecoeffs)
        @assert length(stiffness) == length(lamecoeffs) == 2
        new(stiffness,lamecoeffs)
    end
end

function HookeStiffness(lambda1, mu1, lambda2, mu2)
    stiffness1 = plane_strain_voigt_hooke_matrix(lambda1, mu1)
    stiffness2 = plane_strain_voigt_hooke_matrix(lambda2, mu2)
    lamecoeffs = [(lambda1,mu1),(lambda2,mu2)]
    return HookeStiffness([stiffness1, stiffness2],lamecoeffs)
end

function Base.getindex(hs::HookeStiffness, s)
    row = cell_sign_to_row(s)
    return hs.stiffness[row]
end

function Base.show(io::IO,stiffness::HookeStiffness)
    lambda1,mu1 = lame_coefficients(stiffness,+1)
    lambda2,mu2 = lame_coefficients(stiffness,-1)
    str = "HookeStiffness\n\tlambda1, mu1 = $lambda1, $mu1\n\tlambda2, mu2 = $lambda2, $mu2"
    print(io,str)
end

function lame_coefficients(hs::HookeStiffness,s)
    row = cell_sign_to_row(s)
    return hs.lamecoeffs[row]
end

function plane_strain_voigt_hooke_matrix(lambda, mu)
    l2mu = lambda + 2mu
    return [
        l2mu lambda 0.0
        lambda l2mu 0.0
        0.0 0.0 mu
    ]
end


text:
"""
`ParametersCGVI`: Parameters for right-hand side function of continuous Galerkin variational Integrator.

### Parameters

* `Θ`: function of the noncanonical one-form (∂L/∂v)
* `f`: function of the force (∂L/∂q)
* `Δt`: time step
* `b`: weights of the quadrature rule
* `c`: nodes of the quadrature rule
* `x`: nodes of the basis
* `m`: mass matrix
* `a`: derivative matrix
* `r₀`: reconstruction coefficients at the beginning of the interval
* `r₁`: reconstruction coefficients at the end of the interval
* `t`: initial time
* `q`: solution of q at time t
* `p`: solution of p at time t
"""
mutable struct ParametersCGVI{DT,TT,D,S,R,ΘT,FT} <: Parameters{DT,TT}
    Θ::ΘT
    f::FT

    Δt::TT

    b::Vector{TT}
    c::Vector{TT}

    x::Vector{TT}
    m::Matrix{TT}
    a::Matrix{TT}
    r₀::Vector{TT}
    r₁::Vector{TT}

    t::TT
    q::Vector{DT}
    p::Vector{DT}

    function ParametersCGVI{DT,TT,D,S,R,ΘT,FT}(Θ::ΘT, f::FT, Δt::TT, b, c, x, m, a, r₀, r₁) where {DT,TT,D,S,R,ΘT,FT}
        new(Θ, f, Δt, b, c, x, m, a, r₀, r₁, zero(TT), zeros(DT,D), zeros(DT,D))
    end
end

function ParametersCGVI(equ::IODE{DT,TT,ΘT,FT}, Δt::TT, b, c, x, m, a, r₀, r₁) where {DT,TT,ΘT,FT}
    ParametersCGVI{DT,TT,ndims(equ),length(x),length(c),ΘT,FT}(equ.ϑ, equ.f, Δt, b, c, x, m, a, r₀, r₁)
end

struct IntegratorCacheCGVI{ST,D,S,R} <: IODEIntegratorCache{ST,D}
    q̃::Vector{ST}
    p̃::Vector{ST}
    ṽ::Vector{ST}
    f̃::Vector{ST}
    s̃::Vector{ST}

    X::Vector{Vector{ST}}
    Q::Vector{Vector{ST}}
    P::Vector{Vector{ST}}
    V::Vector{Vector{ST}}
    F::Vector{Vector{ST}}

    function IntegratorCacheCGVI{ST,D,S,R}() where {ST,D,S,R}
        # create temporary vectors
        q̃ = zeros(ST,D)
        p̃ = zeros(ST,D)
        ṽ = zeros(ST,D)
        f̃ = zeros(ST,D)
        s̃ = zeros(ST,D)

        # create internal stage vectors
        X = create_internal_stage_vector(ST,D,S)
        Q = create_internal_stage_vector(ST,D,R)
        P = create_internal_stage_vector(ST,D,R)
        V = create_internal_stage_vector(ST,D,R)
        F = create_internal_stage_vector(ST,D,R)

        new(q̃, p̃, ṽ, f̃, s̃, X, Q, P, V, F)
    end
end

function update_params!(params::ParametersCGVI, sol::AtomicSolutionPODE)
    # set time for nonlinear solver and copy previous solution
    params.t  = sol.t
    params.q .= sol.q
    params.p .= sol.p
end

"Compute stages of variational partitioned Runge-Kutta methods."
@generated function function_stages!(x::Vector{ST}, b::Vector{ST}, params::ParametersCGVI{DT,TT,D,S,R}) where {ST,DT,TT,D,S,R}
    cache = IntegratorCacheCGVI{ST,D,S,R}()

    quote
        @assert length(x) == length(b)

        compute_stages!(x, $cache, params)

        compute_rhs!(b, $cache.X, $cache.Q, $cache.P, $cache.F, $cache.p̃, params)
    end
end

function compute_stages!(x, cache::IntegratorCacheCGVI, params::ParametersCGVI)
    compute_stages!(x, cache.X, cache.Q, cache.V, cache.P, cache.F, cache.q̃, cache.p̃, params)
end

function compute_stages!(x, X, Q, V, P, F, q, p, params::ParametersCGVI{DT,TT,D,S,R}) where {DT,TT,D,S,R}

    # copy x to X
    for i in eachindex(X)
        for k in eachindex(X[i])
            X[i][k] = x[D*(i-1)+k]
        end
    end

    # copy x to p
    for k in eachindex(p)
        p[k] = x[D*S+k]
    end

    # compute Q
    compute_stages_q!(X, Q, q, params)

    # compute V
    compute_stages_v!(X, V, params)

    # compute P and F
    compute_stages_p!(Q, V, P, F, params)
end

function compute_stages_q!(X::Vector{Vector{ST}}, Q::Vector{Vector{ST}}, q::Vector{ST}, params::ParametersCGVI{DT,TT,D,S,R}) where {ST,DT,TT,D,S,R}
    @assert R == length(Q)
    @assert S == length(X)

    local y::ST

    # compute Q
    for i in eachindex(Q)
        @assert D == length(Q[i])
        for k in eachindex(Q[i])
            y = 0
            for j in eachindex(X)
                @assert D == length(X[j])
                y += params.m[i,j] * X[j][k]
            end
            Q[i][k] = y
        end
    end

    # compute q
    for k in eachindex(q)
        y = 0
        for i in eachindex(X)
            y += params.r₁[i] * X[i][k]
        end
        q[k] = y
    end
end

function compute_stages_v!(X::Vector{Vector{ST}}, V::Vector{Vector{ST}}, params::ParametersCGVI{DT,TT,D,S,R}) where {ST,DT,TT,D,S,R}
    @assert R == length(V)
    @assert S == length(X)

    local y::ST

    # compute V
    for i in eachindex(V)
        @assert D == length(V[i])
        for k in eachindex(V[i])
            y = 0
            for j in eachindex(X)
                @assert D == length(X[j])
                y += params.a[i,j] * X[j][k]
            end
            V[i][k] = y / params.Δt
        end
    end
end

function compute_stages_p!(Q::Vector{Vector{ST}}, V::Vector{Vector{ST}},
                           P::Vector{Vector{ST}}, F::Vector{Vector{ST}},
                           params::ParametersCGVI{DT,TT,D,S,R}) where {ST,DT,TT,D,S,R}

    @assert R == length(Q) == length(V) == length(P) == length(F)

    local tᵢ::TT

    # compute P=ϑ(Q,V) and F=f(Q,V)
    for i in eachindex(Q,V,P,F)
        @assert D == length(Q[i]) == length(V[i]) == length(P[i]) == length(F[i])
        tᵢ = params.t + params.Δt * params.c[i]
        params.Θ(tᵢ, Q[i], V[i], P[i])
        params.f(tᵢ, Q[i], V[i], F[i])
    end
end

function compute_rhs!(b::Vector{ST}, X::Vector{Vector{ST}}, Q::Vector{Vector{ST}},
                                     P::Vector{Vector{ST}}, F::Vector{Vector{ST}},
                                     p::Vector{ST}, params::ParametersCGVI{DT,TT,D,S,R}) where {ST,DT,TT,D,S,R}
    local y::ST
    local z::ST

    # compute b = - [(P-AF)]
    for i in 1:S
        for k in 1:D
            z = 0
            for j in eachindex(P,F)
                z += params.b[j] * params.m[j,i] * F[j][k]
                z += params.b[j] * params.a[j,i] * P[j][k] / params.Δt
            end
            b[D*(i-1)+k] = z - (params.r₁[i] * p[k] - params.r₀[i] * params.p[k]) / params.Δt
        end
    end

    # compute b = - [(q-r₀Q)]
    for k in eachindex(params.q)
        y = 0
        for j in eachindex(X)
            y += params.r₀[j] * X[j][k]
        end
        b[D*S+k] = y - params.q[k]
    end
end

"Continuous Galerkin Variational Integrator."
struct IntegratorCGVI{DT,TT,D,S,R,ΘT,FT,GT,HT,VT,PT,ST,IT,BT<:Basis} <: DeterministicIntegrator{DT,TT}
    equation::IODE{DT,TT,ΘT,FT,GT,HT,VT}

    basis::BT
    quadrature::Quadrature{TT,R}

    Δt::TT

    params::PT
    solver::ST
    iguess::InitialGuessPODE{DT,TT,VT,FT,IT}
    cache::IntegratorCacheCGVI{DT,D,S,R}
end

function IntegratorCGVI(equation::IODE{DT,TT,ΘT,FT,GT,HT,VT}, basis::Basis{TT,P}, quadrature::Quadrature{TT,R}, Δt::TT;
    interpolation=HermiteInterpolation{DT}) where {DT,TT,ΘT,FT,GT,HT,VT,P,R}
    D = equation.d
    S = nbasis(basis)

    N = D*(S+1)

    # create solution vector for nonlinear solver
    x = zeros(DT,N)

    # compute coefficients
    r₀ = zeros(TT, nbasis(basis))
    r₁ = zeros(TT, nbasis(basis))
    m  = zeros(TT, nnodes(quadrature), S)
    a  = zeros(TT, nnodes(quadrature), S)

    for i in 1:nbasis(basis)
        r₀[i] = evaluate(basis, i, zero(TT))
        r₁[i] = evaluate(basis, i, one(TT))
        for j in 1:nnodes(quadrature)
            m[j,i] = evaluate(basis, i, nodes(quadrature)[j])
            a[j,i] = derivative(basis, i, nodes(quadrature)[j])
        end
    end

    if get_config(:verbosity) > 1
        println()
        println("  Continuous Galerkin Variational Integrator")
        println("  ==========================================")
        println()
        println("    c  = ", nodes(quadrature))
        println("    b  = ", weights(quadrature))
        println("    x  = ", nodes(basis))
        println("    m  = ", m)
        println("    a  = ", a)
        println("    r₀ = ", r₀)
        println("    r₁ = ", r₁)
        println()
    end

    # create cache for internal stage vectors and update vectors
    # cache = NonlinearFunctionCacheCGVI{DT}(D,S,R)

    # create params
    params = ParametersCGVI(equation, Δt, weights(quadrature), nodes(quadrature), nodes(basis), m, a, r₀, r₁)

    # create rhs function for nonlinear solver
    function_stages = (x,b) -> function_stages!(x, b, params)

    # create nonlinear solver
    solver = get_config(:nls_solver)(x, function_stages)

    # create initial guess
    iguess = InitialGuessPODE(interpolation, equation, Δt)

    # create cache
    cache = IntegratorCacheCGVI{DT,D,S,R}()

    # create integrator
    IntegratorCGVI{DT, TT, D, S, R, ΘT, FT, GT, HT, VT, typeof(params), typeof(solver), typeof(iguess.int), typeof(basis)}(
                equation, basis, quadrature, Δt, params, solver, iguess, cache)
end

equation(integrator::IntegratorCGVI) = integrator.equation
timestep(integrator::IntegratorCGVI) = integrator.Δt

function create_integrator_cache(int::IntegratorCGVI{DT,TT}) where {DT,TT}
    IntegratorCacheCGVI{DT, TT, ndims(int), nbasis(int.basis), nnodes(int.quadrature)}()
end

function initialize!(int::IntegratorCGVI, sol::AtomicSolutionPODE)
    sol.t̅ = sol.t - timestep(int)

    equation(int).v(sol.t, sol.q, sol.p, sol.v)
    equation(int).f(sol.t, sol.q, sol.p, sol.f)

    initialize!(int.iguess, sol.t, sol.q, sol.p, sol.v, sol.f,
                            sol.t̅, sol.q̅, sol.p̅, sol.v̅, sol.f̅)
end

function initial_guess!(int::IntegratorCGVI{DT,TT,D,S,R}, sol::AtomicSolutionPODE{DT,TT}) where {DT,TT,D,S,R}
    if nnodes(int.basis) > 0
        for i in 1:S
            evaluate!(int.iguess, sol.q, sol.p, sol.v, sol.f,
                                  sol.q̅, sol.p̅, sol.v̅, sol.f̅,
                                  int.cache.q̃, int.cache.p̃,
                                  int.cache.ṽ, int.cache.f̃,
                                  nodes(int.basis)[i], nodes(int.basis)[i])

            for k in 1:D
                int.solver.x[D*(i-1)+k] = int.cache.q̃[k]
            end
        end
    else
        for i in 1:S
            for k in 1:D
                int.solver.x[D*(i-1)+k] = 0
            end
        end
    end

    evaluate!(int.iguess, sol.q, sol.p, sol.v, sol.f,
                          sol.q̅, sol.p̅, sol.v̅, sol.f̅,
                          int.cache.q̃, int.cache.p̃,
                          one(TT), one(TT))

    for k in 1:D
        int.solver.x[D*S+k] = int.cache.p̃[k]
    end
end

function update_solution!(sol::AtomicSolutionPODE, cache::IntegratorCacheCGVI)
    sol.q .= cache.q̃
    sol.p .= cache.p̃
end

"Integrate ODE with variational partitioned Runge-Kutta integrator."
function integrate_step!(int::IntegratorCGVI{DT,TT}, sol::AtomicSolutionPODE{DT,TT}) where {DT,TT}
    # update nonlinear solver parameters from cache
    update_params!(int.params, sol)

    # compute initial guess
    initial_guess!(int, sol)

    # reset cache
    reset!(sol, timestep(int))

    # call nonlinear solver
    solve!(int.solver)

    # print solver status
    print_solver_status(int.solver.status, int.solver.params)

    # check if solution contains NaNs or error bounds are violated
    check_solver_status(int.solver.status, int.solver.params)

    # compute vector fields at internal stages
    compute_stages!(int.solver.x, int.cache, int.params)

    # compute final update
    update_solution!(sol, int.cache)

    # copy solution to initial guess
    update!(int.iguess, sol.t, sol.q, sol.p, sol.v, sol.f)
end


text:
using HopfieldNets
using Pkg

include(Pkg.dir("HopfieldNets", "demo", "letters.jl"))

patterns = hcat(X, O)

n = size(patterns, 1)

net = DiscreteHopfieldNet(n)

train!(net, patterns)

settle!(net, 10, true)

Xcorrupt = copy(X)
for i = 2:7
     Xcorrupt[i] = 1
end

Xrestored = associate!(net, Xcorrupt[:,3])
all(Xcorrupt .== X)
all(Xrestored .== X)


text:
using BenchmarkTools
using LightGraphs
using Ripserer
using Eirene
using Random

include(joinpath(@__DIR__, "utils.jl"))

function gcycle_scaling(Ns=[[10]; 100:100:1000]; seconds=1000, samples=10)
    map(Ns) do N
        dist = floyd_warshall_shortest_paths(cycle_graph(N)).dists

        println("N=$N")
        println("running cohomology")
        bc = @benchmark ripserer(
            $dist
        ) seconds=seconds samples=samples evals=1 gcsample=true
        tc = minimum(bc).time

        println("  running involuted ")
        bi = @benchmark ripserer(
            $dist; alg=:involuted
        ) seconds=seconds samples=samples evals=1 gcsample=true
        ti = minimum(bi).time

        println("  running eirene    ")
        be = @benchmark eirene(
            $dist; model="vr"
        ) seconds=seconds samples=samples evals=1 gcsample=true
        te = minimum(be).time
        println("  done.             ")

        println((;N, tc, ti, te))
        return (;N, tc, ti, te)
    end
end

function random16_scaling(dims=1:8)
    dist = load_data(joinpath(@__DIR__, "..", "datasets", "random16.dist"))

    map(dims) do dim
        println("dim=$dim")

        println("running cohomology")
        bc = @benchmark ripserer(
            $dist; dim_max=$dim
        ) seconds=1000 samples=10 evals=1 gcsample=true
        tc = minimum(bc).time

        println("  running involuted ")
        bi = @benchmark ripserer(
            $dist; dim_max=$dim, alg=:involuted
        ) seconds=1000 samples=10 evals=1 gcsample=true
        ti = minimum(bi).time

        println("  running eirene    ")
        be = @benchmark eirene(
            $dist; model=$("vr"), maxdim=$dim
        ) seconds=1000 samples=10 evals=1 gcsample=true
        te = minimum(be).time
        println("  done.             ")

        (;dim, tc, ti, te)
    end
end

function run_gcycle_and_random16_scaling()
    CSV.write("random16.csv", random16_scaling())
    CSV.write("gcycle.csv", gcycle_scaling())
end

"""
    subset_scaling(data, subsets; eirene_max=subsets[end], samples=10, dim_max=1)

Benchmark Ripserer and Eirene at various subsets of the `data`. `data` is shuffled first,
then increasingly large subsets are taken.
"""
function subset_scaling(
    data, subsets; eirene_max=subsets[end], samples=10, seconds=1000, dim_max=1, threshold=nothing
)
    data = perturb(data)
    if !isnothing(threshold)
        thr_r = (;threshold)
        thr_e = (;maxrad=threshold)
    else
        thr_r = NamedTuple()
        thr_e = NamedTuple()
    end
    map(subsets) do N
        ripserer_data = subset(data, N)
        if ripserer_data isa Vector
            eirene_data = Ripserer.to_matrix(ripserer_data)
            model="pc"
        else
            eirene_data = ripserer_data
            model="vr"
        end

        println("  running cohomology on $N points")
        bc = @benchmark ripserer(
            $ripserer_data, dim_max=$dim_max, $thr_r...
        ) seconds=seconds samples=samples evals=1 gcsample=true
        tc = minimum(bc).time

        println("  running involuted on $N points")
        bi = @benchmark ripserer(
            $ripserer_data; alg=:involuted, dim_max=$dim_max, $thr_r...
        ) seconds=seconds samples=samples evals=1 gcsample=true
        ti = minimum(bi).time

        if N ≤ eirene_max
            println("  running eirene on $N points")
            be = @benchmark eirene(
                $eirene_data; model=$model, maxdim=$dim_max, $thr_e...
            ) seconds=seconds samples=samples evals=1 gcsample=true
            te = minimum(be).time
        else
            te = nothing
        end

        (;N, tc, ti, te)
     end
end

"""
    dim_scaling(data, dims; eirene_max, samples=10)

Run Ripserer and Eirene on the same data in various dimensions.
"""
function dim_scaling(data, dims; eirene_max, samples=10, seconds=1000, threshold=nothing)
    if !isnothing(threshold)
        thr_r = (;threshold)
        thr_e = (;maxrad=threshold)
    else
        thr_r = NamedTuple()
        thr_e = NamedTuple()
    end
    map(dims) do dim_max
        model = data isa AbstractMatrix ? "vr" : "pc"
        bc = @benchmark ripserer(
            $data; dim_max=$dim_max, $thr_r...
        ) seconds=seconds samples=samples evals=1 gcsample=true
        tc = minimum(bc).time

        println("  running involuted ")
        bi = @benchmark ripserer(
            $data; alg=:involuted, dim_max=$dim_max, $thr_r...
        ) seconds=seconds samples=samples evals=1 gcsample=true
        ti = minimum(bi).time

        println("  running eirene    ")
        be = @benchmark eirene(
            $data; model=$model, maxdim=$dim_max, $thr_e...
        ) seconds=seconds samples=samples evals=1 gcsample=true
        te = minimum(be).time

        (;N, tc, ti, te)
     end
end

perturb(v::AbstractVector) = shuffle(v)
function perturb(m::AbstractMatrix)
    ord = shuffle(1:size(m, 1))
    return m[ord, ord]
end
subset(v::AbstractVector, N) = v[1:N]
subset(m::AbstractMatrix, N) = m[1:N, 1:N]

function run_subset_scaling()
    for (file,              name,         dim_max, threshold, subsets) in [
        ("dragon1000.dist", "dragon1000", 1,       nothing,   100:100:1000),
        ("o3_1024.pts",     "o3_1024",    3,       1.8,       124:100:1024),
        ("o3_4096.pts",     "o3_4096",    3,       1.4,       196:100:4096),
        ("hiv.dist",        "hiv",        1,       nothing,   188:100:1088),
    ]
        println(" ++ $name ++")
        data = load_data(joinpath(@__DIR__, "..", "datasets", file))
        res = subset_scaling(
            data, subsets; eirene_max=subsets[end], samples=10, dim_max=dim_max
        )
        CSV.write("$name.csv", res)
    end
end


text:
module SpinGlassDynamics

using MKL
using LinearAlgebra
using Distributions
using SpinGlassNetworks
using DifferentialEquations
using DocStringExtensions

include("CIM.jl")

end #module


text:
module Mise

# Base
using Base.Filesystem
using Base.Iterators: partition, zip
using LinearAlgebra
using Printf
using Random
using Statistics

# Math
using NumericalIntegration
using Dierckx

# Statistic Tools

# Statistics
using StatsBase: mean, std, mean_and_std, zscore, crosscor
import StatsBase: zscore!
using StatsBase
using Distributions
using KernelDensity
using CurveFit

# Statistical Models
using TimeSeries
using HypothesisTests
using Loess
using StateSpaceModels

# DiffEq
using DifferentialEquations

# ML
using Flux
using Zygote

# GPU
using CUDA

# MCMC
#using Mamba

# Tables
using DataFrames, DataFramesMeta, Missings, Query
using DataValues
using Dates, TimeZones
using JuliaDB

# IO
using BSON, CSV
using DelimitedFiles
using FileIO
using Formatting
using Glob

# utils
using ArgParse
using MicroLogging
using ProgressMeter

# Plots
using Plots
using Plots.PlotMeasures
using Colors, ColorTypes
using StatsPlots: @df, StatsPlots

if isa(Sys.which("python3"), String)
    using ExcelReaders
end

ENV["GKSwstype"] = "100"

# to use Plots in headless system
# https://github.com/JuliaPlots/Plots.jl/issues/1076#issuecomment-327509819
include("utils.jl")
include("data.jl")
include("input.jl")
include("output.jl")
include("seasonality.jl")
include("loss.jl")
include("activation.jl")
include("evaluation.jl")
include("plots.jl")

# Statistical Methods
include("analysis/analysis.jl")
include("analysis/plot_analysis.jl")
include("ARIMA/ARIMA.jl")
include("ARIMA/plot_ARIMA.jl")
include("OU/model.jl")
include("OU_season/model.jl")

# Machine Learning
include("DNN/preprocess.jl")
include("DNN/model.jl")
include("DNN_season/model.jl")
include("RNN/model.jl")
#include("RNN_Attention/model.jl")

include("postprocess/post_DNN.jl")

# input
export join_data, filter_raw_data, filter_station, filter_jongro, read_jongro,
     parse_aerosols,
# utils
    extract_col_statvals, zscore!, unzscore, unzscore!,
    minmax_scaling, minmax_scaling!,
    unminmax_scaling, unminmax_scaling!,
    exclude_elem, findrow,
    WHO_PM10, WHO_PM25,
# data
    get_date_range, validate_dates, window_df,
    split_sizes3, split_sizes2,
    remove_sparse_input!, is_sparse_Y,
    getX, getY, make_pair_DNN, make_pair_date_DNN, make_batch_DNN,
    serializeBatch, whcn2cnh, matrix2arrays,
    make_pair_RNN, make_pair_date_RNN, make_batch_RNN, make_batch_date_RNN,
    zero2Missing!,
    padded_push!, construct_annual_table, compose_seasonality,
# activation
# evaluation
    evaluations1, evaluations2, evaluation1, evaluation2,
    RMSE, MAE, MSPE, MAPE,
    RSR, NSE, PBIAS, IOA, RefinedIOA, R2, AdjR2, classification,
# loss
    huber_loss, huber_loss_mean, mse_rnn,
# smoothing
    periodic_mean, populate_periodic_mean,
    season_adj_lee, compute_annual_mean,
    decompose_seasonality,
# analysis & plot
    compute_inttscale, mean_aucotor, pdf,
    plot_anal_lineplot,
    plot_anal_pdf, plot_anal_autocor,
    plot_anal_correlogram, plot_anal_violin,
    plot_anal_periodic_mean, plot_anal_periodic_fluc,
# ARIMA
    smoothing_mean,
    predict_ARIMA, plot_ARIMA_acf,
    plot_seasonality,
    plot_ARIMA_fluc, plot_ARIMA1_mean, plot_ARIMA_mean_smoothing,
# preprocess
    load_data_DNN, filter_station_DNN, process_raw_data_DNN!, read_station,
# DNN
    train_DNN, train_season_DNN, corr_input,
# RNN
    train_RNN, train_Attention,
# OU
    evolve_OU, evolve_OU_season,
# ARIMA
    predict_ARIMA,
# post processing
    compute_corr, test_features, test_station, test_classification,
# output
    predict_DNN_model_zscore, predict_DNN_model_zscore_season,
    predict_DNN_model_minmax, predict_DNN_model_minmax_season,
    predict_DNN_model_logzscore, predict_DNN_model_logzscore_season,
    predict_DNN_model_logminmax, predict_DNN_model_logminmax_season,
    predict_DNN_model_invzscore, predict_DNN_model_invzscore_season,
    predict_RNN_model_zscore, predict_RNN_model_minmax,
    export_CSV,
# plot
    plot_histogram,
    plot_lineplot_total,
    plot_corr_input,
    plot_pcorr,
    plot_corr,
    plot_ARIMA_corr,
    plot_ARIMA_scatter,
    plot_ARIMA_lineplot,
    plot_DNN_scatter,
    plot_DNN_histogram,
    plot_DNN_lineplot,
    plot_RNN_lineplot
end


text:
h₁(i::Int) = max(6 - abs(i-7), 0)
h₂(i::Int) = h₁(i - 8)
h₃(i::Int) = h₁(i - 4)

"""
  triangular_test_function(c1::Int, c2::Int, c3::Int, L::Int=32)

Generates a set of triangluar test functions with 3 classes.
"""
function triangular_test_functions(c1::Int, c2::Int, c3::Int; L::Int=32, shuffle::Bool=false)
  @assert c1 >= 0
  @assert c2 >= 0
  @assert c3 >= 0

  u = rand(Uniform(0,1),1)[1]
  ϵ = rand(Normal(0,1),(L,c1+c2+c3))

  y = string.(vcat(ones(c1), ones(c2) .+ 1, ones(c3) .+ 2))

  H₁ = Array{Float64,2}(undef,L,c1)
  H₂ = Array{Float64,2}(undef,L,c2)
  H₃ = Array{Float64,2}(undef,L,c3)
  for i in 1:L
    H₁[i,:] .= u * h₁(i) + (1 - u) * h₂(i)
    H₂[i,:] .= u * h₁(i) + (1 - u) * h₃(i)
    H₃[i,:] .= u * h₂(i) + (1 - u) * h₃(i)
  end

  H = hcat(H₁, H₂, H₃) + ϵ

  if shuffle
    idx = [1:(c1+c2+c3)...]
    shuffle!(idx)
    return H[:,idx], y[idx]
  end

  return H, y
end

max_stars_repo_name: BoundaryValueProblems/LDBExperiments

max_issues_repo_name: BoundaryValueProblems/LDBExperiments

max_forks_repo_name: BoundaryValueProblems/LDBExperiments