text:
mutable struct VirtualGlobalVariable
    # actual profiled type
    t::Any
    # keeps `id` of `JETInterpreter` that defined this
    id::Symbol
    # `Symbol` of a dummy generic function that generates dummy backedge (i.e. `li`)
    edge_sym::Symbol
    # dummy backedge, which will be invalidated on update of `t`
    li::MethodInstance
    # whether this virtual global variable is declarared as constant or not
    iscd::Bool

    function VirtualGlobalVariable(@nospecialize(t),
                                   id::Symbol,
                                   edge_sym::Symbol,
                                   li::MethodInstance,
                                   iscd::Bool,
                                   )
        return new(t, id, edge_sym, li, iscd)
    end
end

"""
    function overload_abstract_call_gf_by_type!()
        ...
    end
    push_inithook!(overload_abstract_call_gf_by_type!)

the aims of this overload are:
1. report `NoMethodErrorReport` on empty method signature matching
2. keep inference on non-concrete call sites in toplevel frame created by [`virtual_process!`](@ref)
3. don't bail out even after the current return type grows up to `Any` and collect as much
   error points as possible; of course it slows down inference performance, but hopefully it
   stays to be "practical" speed (because the number of matching methods is limited beforehand)
4. force constant prop' even if the inference result can't be improved anymore when `rettype`
   is already `Const`; this is because constant prop' can still produce more "correct"
   analysis by throwing away the error reports in the callee frames
5. always add backedges (even if a new method can't refine the return type grew up to`Any`),
   because a new method always may change the JET analysis result
"""
function overload_abstract_call_gf_by_type!()

# %% for easier interactive update of abstract_call_gf_by_type
Core.eval(CC, quote

# TODO:
# - report "too many method matched"
# - maybe "cound not identify method table for call" won't happen since we eagerly propagate bottom for e.g. undef var case, etc.

function abstract_call_gf_by_type(interp::$(JETInterpreter), @nospecialize(f), argtypes::Vector{Any}, @nospecialize(atype), sv::InferenceState,
                                  max_methods::Int = InferenceParams(interp).MAX_METHODS)
    if sv.params.unoptimize_throw_blocks && sv.currpc in sv.throw_blocks
        return CallMeta(Any, false)
    end
    valid_worlds = WorldRange()
    atype_params = unwrap_unionall(atype).parameters
    splitunions = 1 < unionsplitcost(atype_params) <= InferenceParams(interp).MAX_UNION_SPLITTING
    mts = Core.MethodTable[]
    fullmatch = Bool[]
    if splitunions
        splitsigs = switchtupleunion(atype)
        applicable = Any[]
        infos = MethodMatchInfo[]
        for sig_n in splitsigs
            mt = ccall(:jl_method_table_for, Any, (Any,), sig_n)
            if mt === nothing
                add_remark!(interp, sv, "Could not identify method table for call")
                return CallMeta(Any, false)
            end
            mt = mt::Core.MethodTable
            matches = findall(sig_n, method_table(interp); limit=max_methods)
            if matches === missing
                add_remark!(interp, sv, "For one of the union split cases, too many methods matched")
                return CallMeta(Any, false)
            end
            #=== abstract_call_gf_by_type patch point 1-1 start ===#
            info = MethodMatchInfo(matches)
            if $(is_empty_match)(info)
                # report `NoMethodErrorReport` for union-split signatures
                $(report!)(interp, $(NoMethodErrorReport)(interp, sv, true, atype))
            end
            push!(infos, info)
            #=== abstract_call_gf_by_type patch point 1-1 end ===#
            append!(applicable, matches)
            valid_worlds = intersect(valid_worlds, matches.valid_worlds)
            thisfullmatch = _any(match->(match::MethodMatch).fully_covers, matches)
            found = false
            for (i, mt′) in enumerate(mts)
                if mt′ === mt
                    fullmatch[i] &= thisfullmatch
                    found = true
                    break
                end
            end
            if !found
                push!(mts, mt)
                push!(fullmatch, thisfullmatch)
            end
        end
        info = UnionSplitInfo(infos)
    else
        mt = ccall(:jl_method_table_for, Any, (Any,), atype)
        if mt === nothing
            add_remark!(interp, sv, "Could not identify method table for call")
            return CallMeta(Any, false)
        end
        mt = mt::Core.MethodTable
        matches = findall(atype, method_table(interp, sv); limit=max_methods)
        if matches === missing
            # this means too many methods matched
            # (assume this will always be true, so we don't compute / update valid age in this case)
            add_remark!(interp, sv, "Too many methods matched")
            return CallMeta(Any, false)
        end
        push!(mts, mt)
        push!(fullmatch, _any(match->(match::MethodMatch).fully_covers, matches))
        info = MethodMatchInfo(matches)
        #=== abstract_call_gf_by_type patch point 1-2 start ===#
        if $(is_empty_match)(info)
            # report `NoMethodErrorReport` for this call signature
            $(report!)(interp, $(NoMethodErrorReport)(interp, sv, false, atype))
        end
        #=== abstract_call_gf_by_type patch point 1-2 end ===#
        applicable = matches.matches
        valid_worlds = matches.valid_worlds
    end
    update_valid_age!(sv, valid_worlds)
    applicable = applicable::Array{Any,1}
    napplicable = length(applicable)
    rettype = Bottom
    edgecycle = false
    edges = Any[]
    nonbot = 0  # the index of the only non-Bottom inference result if > 0
    seen = 0    # number of signatures actually inferred
    istoplevel = sv.linfo.def isa Module
    multiple_matches = napplicable > 1

    if f !== nothing && napplicable == 1 && is_method_pure(applicable[1]::MethodMatch)
        val = pure_eval_call(f, argtypes)
        if val !== false
            # TODO: add some sort of edge(s)
            return CallMeta(val, MethodResultPure())
        end
    end

    #=== abstract_call_gf_by_type patch point 4-1 start ===#
    nreports = length(interp.reports)
    #=== abstract_call_gf_by_type patch point 4-1 end ===#

    for i in 1:napplicable
        match = applicable[i]::MethodMatch
        method = match.method
        sig = match.spec_types
        #=== abstract_call_gf_by_type patch point 2 start ===#
        if istoplevel && !isdispatchtuple(sig) && !$(istoplevel)(sv) # keep going for "our" toplevel frame
        #=== abstract_call_gf_by_type patch point 2 end ===#
            # only infer concrete call sites in top-level expressions
            add_remark!(interp, sv, "Refusing to infer non-concrete call site in top-level expression")
            rettype = Any
            break
        end
        sigtuple = unwrap_unionall(sig)::DataType
        splitunions = false
        this_rt = Bottom
        # TODO: splitunions = 1 < unionsplitcost(sigtuple.parameters) * napplicable <= InferenceParams(interp).MAX_UNION_SPLITTING
        # currently this triggers a bug in inference recursion detection
        if splitunions
            splitsigs = switchtupleunion(sig)
            for sig_n in splitsigs
                rt, edgecycle1, edge = abstract_call_method(interp, method, sig_n, svec(), multiple_matches, sv)
                if edge !== nothing
                    push!(edges, edge)
                end
                edgecycle |= edgecycle1::Bool
                this_rt = tmerge(this_rt, rt)
                #=== abstract_call_gf_by_type patch point 3-1 start ===#
                # this_rt === Any && break # keep going and collect as much error reports as possible
                #=== abstract_call_gf_by_type patch point 3-1 end ===#
            end
        else
            this_rt, edgecycle1, edge = abstract_call_method(interp, method, sig, match.sparams, multiple_matches, sv)
            edgecycle |= edgecycle1::Bool
            if edge !== nothing
                push!(edges, edge)
            end
        end
        if this_rt !== Bottom
            if nonbot === 0
                nonbot = i
            else
                nonbot = -1
            end
        end
        seen += 1
        rettype = tmerge(rettype, this_rt)
        #=== abstract_call_gf_by_type patch point 3-2 start ===#
        # rettype === Any && break # keep going and collect as much error reports as possible
        #=== abstract_call_gf_by_type patch point 3-2 end ===#
    end

    #=== abstract_call_gf_by_type patch point 4-2 start ===#
    # check if constant propagation can improve analysis by throwing away possibly false positive reports
    has_been_reported = (length(interp.reports) - nreports) > 0
    #=== abstract_call_gf_by_type patch point 4-2 end ===#

    # try constant propagation if only 1 method is inferred to non-Bottom
    # this is in preparation for inlining, or improving the return result
    is_unused = call_result_unused(sv)
    #=== abstract_call_gf_by_type patch point 4-3 start ===#
    if nonbot > 0 && seen == napplicable && (!edgecycle || !is_unused) &&
            (is_improvable(rettype) || has_been_reported) && InferenceParams(interp).ipo_constant_propagation
    #=== abstract_call_gf_by_type patch point 4-3 end ===#
        # if there's a possibility we could constant-propagate a better result
        # (hopefully without doing too much work), try to do that now
        # TODO: it feels like this could be better integrated into abstract_call_method / typeinf_edge
        const_rettype = abstract_call_method_with_const_args(interp, rettype, f, argtypes, applicable[nonbot]::MethodMatch, sv, edgecycle)
        if const_rettype ⊑ rettype
            # use the better result, if it's a refinement of rettype
            rettype = const_rettype
        end
    end
    if is_unused && !(rettype === Bottom)
        add_remark!(interp, sv, "Call result type was widened because the return value is unused")
        # We're mainly only here because the optimizer might want this code,
        # but we ourselves locally don't typically care about it locally
        # (beyond checking if it always throws).
        # So avoid adding an edge, since we don't want to bother attempting
        # to improve our result even if it does change (to always throw),
        # and avoid keeping track of a more complex result type.
        rettype = Any
    end
    #=== abstract_call_gf_by_type patch point 5 start ===#
    # a new method may refine analysis, so we always add backedges
    if true # !(rettype === Any) # adding a new method couldn't refine (widen) this type
    #=== abstract_call_gf_by_type patch point 5 end ===#
        for edge in edges
            add_backedge!(edge::MethodInstance, sv)
        end
        for (thisfullmatch, mt) in zip(fullmatch, mts)
            if !thisfullmatch
                # also need an edge to the method table in case something gets
                # added that did not intersect with any existing method
                add_mt_backedge!(mt, atype, sv)
            end
        end
    end
    #print("=> ", rettype, "\n")
    return CallMeta(rettype, info)
end

end) # Core.eval(CC, quote
# %% for easier interactive update of abstract_call_gf_by_type

end # function overload_abstract_call_gf_by_type!()
push_inithook!(overload_abstract_call_gf_by_type!)

function is_empty_match(info)
    res = info.results
    isa(res, MethodLookupResult) || return false # when does this happen ?
    return isempty(res.matches)
end

"""
    function overload_abstract_call_method_with_const_args!()
        ...
    end
    push_inithook!(overload_abstract_call_method_with_const_args!)

the aim of this overloads is:
1. force constant prop' even if the inference result can't be improved anymore when `rettype`
   is already `Const`; this is because constant prop' can still produce more "correct"
   analysis by throwing away the error reports in the callee frames
2. when local cache for constant analysis is hit, reconstruct the "after-constant-prop' analysis"
   state using `JETInterpreter`'s local report caches
"""
function overload_abstract_call_method_with_const_args!()

# %% for easier interactive update of abstract_call_method_with_const_args
Core.eval(CC, quote

function abstract_call_method_with_const_args(interp::$(JETInterpreter), @nospecialize(rettype), @nospecialize(f), argtypes::Vector{Any}, match::MethodMatch, sv::InferenceState, edgecycle::Bool)
    method = match.method
    nargs::Int = method.nargs
    method.isva && (nargs -= 1)
    length(argtypes) >= nargs || return Any
    haveconst = false
    allconst = true
    # see if any or all of the arguments are constant and propagating constants may be worthwhile
    for a in argtypes
        a = widenconditional(a)
        if allconst && !isa(a, Const) && !isconstType(a) && !isa(a, PartialStruct)
            allconst = false
        end
        if !haveconst && has_nontrivial_const_info(a) && const_prop_profitable(a)
            haveconst = true
        end
        if haveconst && !allconst
            break
        end
    end
    #=== abstract_call_method_with_const_args patch point 1 start ===#
    # force constant propagation even if it doesn't improve return type;
    # constant prop' may improve report accuracy
    haveconst || #= improvable_via_constant_propagation(rettype) || =# return Any
    #=== abstract_call_method_with_const_args patch point 1 end ===#
    if nargs > 1
        if istopfunction(f, :getindex) || istopfunction(f, :setindex!)
            arrty = argtypes[2]
            # don't propagate constant index into indexing of non-constant array
            if arrty isa Type && arrty <: AbstractArray && !issingletontype(arrty)
                return Any
            elseif arrty ⊑ Array
                return Any
            end
        elseif istopfunction(f, :iterate)
            itrty = argtypes[2]
            if itrty ⊑ Array
                return Any
            end
        end
    end
    if !allconst && (istopfunction(f, :+) || istopfunction(f, :-) || istopfunction(f, :*) ||
                     istopfunction(f, :(==)) || istopfunction(f, :!=) ||
                     istopfunction(f, :<=) || istopfunction(f, :>=) || istopfunction(f, :<) || istopfunction(f, :>) ||
                     istopfunction(f, :<<) || istopfunction(f, :>>))
        return Any
    end
    force_inference = allconst || InferenceParams(interp).aggressive_constant_propagation
    if istopfunction(f, :getproperty) || istopfunction(f, :setproperty!)
        force_inference = true
    end
    mi = specialize_method(match, !force_inference)
    mi === nothing && return Any
    mi = mi::MethodInstance
    # decide if it's likely to be worthwhile
    if !force_inference && !const_prop_heuristic(interp, method, mi)
        return Any
    end
    inf_cache = get_inference_cache(interp)
    inf_result = cache_lookup(mi, argtypes, inf_cache)
    if inf_result === nothing
        if edgecycle
            # if there might be a cycle, check to make sure we don't end up
            # calling ourselves here.
            infstate = sv
            cyclei = 0
            while !(infstate === nothing)
                if method === infstate.linfo.def && any(infstate.result.overridden_by_const)
                    return Any
                end
                if cyclei < length(infstate.callers_in_cycle)
                    cyclei += 1
                    infstate = infstate.callers_in_cycle[cyclei]
                else
                    cyclei = 0
                    infstate = infstate.parent
                end
            end
        end
        inf_result = InferenceResult(mi, argtypes)
        frame = InferenceState(inf_result, #=cache=#false, interp)
        frame === nothing && return Any # this is probably a bad generated function (unsound), but just ignore it
        frame.limited = true
        frame.parent = sv
        push!(inf_cache, inf_result)
        typeinf(interp, frame) || return Any
    #=== abstract_call_method_with_const_args patch point 2 start ===#
    else
        # local cache for this constant analysis is hit
        inf_result = inf_result::InferenceResult
        if !isa(inf_result.result, InferenceState)
            # corresponds to report throw away logic in `_typeinf(interp::JETInterpreter, frame::InferenceState)`
            $(filter!)(r->$(!is_lineage)(r.lineage, sv, inf_result.linfo), interp.reports)

            local_cache = $(get)(interp.cache, argtypes, nothing)
            if isa(local_cache, $(Vector{InferenceErrorReportCache}))
                $(foreach)(local_cache) do cached
                    $(restore_cached_report!)(cached, interp, sv)
                end
            end
        end
    #=== abstract_call_method_with_const_args patch point 2 end ===#
    end
    result = inf_result.result
    # if constant inference hits a cycle, just bail out
    isa(result, InferenceState) && return Any
    add_backedge!(inf_result.linfo, sv)
    return result
end

end) # Core.eval(CC, quote
# %% for easier interactive update of abstract_call_method_with_const_args

end # function overload_abstract_call_method_with_const_args!()
push_inithook!(overload_abstract_call_method_with_const_args!)

function CC.abstract_eval_special_value(interp::JETInterpreter, @nospecialize(e), vtypes::VarTable, sv::InferenceState)
    ret = @invoke abstract_eval_special_value(interp::AbstractInterpreter, e, vtypes::VarTable, sv::InferenceState)

    if isa(ret, Const)
        # unwrap virtual global variable to actual type
        val = ret.val
        if isa(val, VirtualGlobalVariable)
            # add dummy backedge, which will be invalidated on update of this vitual global variable
            add_backedge!(val.li, sv)

            ret = val.t
        end
    elseif isa(e, GlobalRef)
        mod, name = e.mod, e.name
        if isdefined(mod, name)
            # we don't track types of global variables except when we're in toplevel frame,
            # and here we just annotate this as `Any`; NOTE: this is becasue:
            # - we can't track side-effects of assignments of global variables that happen in
            #   (possibly deeply nested) callees, and it might be possible if we just ignore
            #   assignments happens in callees that aren't reached by type inference by
            #   the widening heuristics
            # - consistency with Julia's native type inference
            # - it's hard to track side effects for cached frames

            # TODO: add report pass here (for performance linting)

            # special case and propagate `Main` module as constant adn this is somewhat
            # critical for performance when self-profiling
            if name === :Main
                ret = Const(Main)
            end
        else
            # report access to undefined global variable
            report!(interp, GlobalUndefVarErrorReport(interp, sv, mod, name))

            # `ret` at this point should be annotated as `Any` by `NativeInterpreter`, and
            # we just pass it as is to collect as much error points as possible within this
            # frame
            # IDEA: we can change it to `Bottom` to suppress any further abstract interpretation
            # with this variable, and the later analysis after the update on this (currently)
            # undefined variable just works because we will invalidate the cache for this frame
            # anyway
        end
    end

    return ret
end

function CC.abstract_eval_value(interp::JETInterpreter, @nospecialize(e), vtypes::VarTable, sv::InferenceState)
    ret = @invoke abstract_eval_value(interp::AbstractInterpreter, e, vtypes::VarTable, sv::InferenceState)

    # report non-boolean condition error
    stmt = get_cur_stmt(sv)
    if isa(stmt, GotoIfNot)
        t = widenconst(ret)
        if t !== Bottom && !⊑(Bool, t)
            report!(interp, NonBooleanCondErrorReport(interp, sv, t))
            ret = Bottom
        end
    end

    return ret
end

function CC.abstract_eval_statement(interp::JETInterpreter, @nospecialize(e), vtypes::VarTable, sv::InferenceState)
    if istoplevel(sv)
        if interp.concretized[get_cur_pc(sv)]
            return Any # bail out if it has been interpreted by `ConcreteInterpreter`
        end
    end

    ret = @invoke abstract_eval_statement(interp::AbstractInterpreter, e, vtypes::VarTable, sv::InferenceState)

    # assign virtual global variable
    # for toplevel frames, we do virtual global variable assignments, whose types are
    # propagated even if they're non-constant
    if istoplevel(sv)
        stmt = get_cur_stmt(sv)
        if isexpr(stmt, :(=))
            lhs = first(stmt.args)

            if isa(lhs, GlobalRef)
                set_virtual_globalvar!(interp, lhs.mod, lhs.name, ret, sv)
            end
        end
    end

    return ret
end

function set_virtual_globalvar!(interp, mod, name, @nospecialize(t), sv)
    local update::Bool = false
    id = get_id(interp)

    iscd = is_constant_declared(mod, name, sv.src.code)

    t′, id′, (edge_sym, li) = if isdefined(mod, name)
        val = getfield(mod, name)
        if isa(val, VirtualGlobalVariable)
            t′ = val.t
            if val.iscd && widenconst(t′) !== widenconst(t)
                report!(interp, InvalidConstantRedefinition(interp, sv, mod, name, widenconst(t′), widenconst(t)))
                return
            end

            # update previously-defined virtual global variable
            update = true
            t′, val.id, (val.edge_sym, val.li)
        else
            if isconst(mod, name)
                t′ = typeof(val)
                if t′ !== widenconst(t)
                    report!(interp, InvalidConstantRedefinition(interp, sv, mod, name, t′, widenconst(t)))
                    return
                end
            end

            # this pass hopefully won't happen within the current design
            @warn "JET.jl can't trace updates of global variable that already have values" mod name val
            return
        end
    else
        # define new virtual global variable
        Bottom, id, gen_dummy_backedge(mod)
    end

    # if this is constant declared and it's value is known to be constant, let's concretize
    # it for good reasons; this will help us analyse on code with global type aliases, etc.
    if !(t ⊑ VirtualGlobalVariable) && isa(t, Const) && iscd
        return Core.eval(mod, :(const $(name) = $(t.val)))
    end

    # at this point undefined slots may still be annotated as `NOT_FOUND` (e.g. when there is
    # undefined slot, etc.), which can't be used as abstract value for later profiling,
    # so replace it with Bottom
    if t === NOT_FOUND
        t = Bottom
    end

    if begin
            # if the previous virtual global variable assignment happened in the same inference process,
            # JET needs to perform type merge
            id === id′ ||
            # if this assignment happens in an non-deterministic way, we still need to perform type merge
            # NOTE: this may happen multiple times for the same statement (within an iteration for
            # maximum fixed point computation), so pre-computing basic blocks before entering a toplevel
            # inference frame might be better
            is_nondeterministic(sv)
        end
        t = tmerge(t′, t)
    end

    if id !== id′
        # invalidate the dummy backedge that is bound to this virtual global variable,
        # so that depending `MethodInstance` will run fresh type inference on the next hit
        li = force_invalidate!(mod, edge_sym)
    end

    ex = if update
        quote
            local name = $(name)
            name.t = $(t)
            name.id = $(QuoteNode(id))
            name.edge_sym = $(QuoteNode(edge_sym))
            name.li = $(li)
            name
        end
    else
        vgv = VirtualGlobalVariable(t, id, edge_sym, li, iscd)
        :(const $(name) = $(vgv))
    end
    return Core.eval(mod, ex)::VirtualGlobalVariable
end

function is_constant_declared(mod, name, stmts)
    # `fix_global_symbols!` replaces all the symbols in a toplevel frame with `GlobalRef`
    gr = GlobalRef(mod, name)

    return any(stmts) do x
        if isexpr(x, :const)
            arg = first((x::Expr).args)
            isa(arg, GlobalRef) && return arg == gr
        end
        return false
    end
end

is_nondeterministic(sv) = is_nondeterministic(get_cur_pc(sv), compute_basic_blocks(sv.src.code))

# XXX: does this approach really cover all the control flow ?
function is_nondeterministic(pc, bbs)
    isnd = false

    for (idx, bb) in enumerate(bbs.blocks)
        if pc in bb.stmts
            for bb′ in bbs.blocks
                if idx in bb′.succs
                    isnd |= length(bb′.succs) > 1
                end
            end
        end
    end

    return isnd
end

function gen_dummy_backedge(mod)
    @gensym edge_sym
    return edge_sym, force_invalidate!(mod, edge_sym) # just generate dummy `MethodInstance` to be invalidated
end

# TODO: find a more fine-grained way to do this ? re-evaluating an entire function seems to be over-kill for this
function force_invalidate!(mod, edge_sym)
    λ = Core.eval(mod, :($(edge_sym)() = return))::Function
    m = first(methods(λ))
    return specialize_method(m, Tuple{typeof(λ)}, svec())::MethodInstance
end


text:
@testset "Conversion from and to MATLAB model" begin
    filename = model_paths["iJO1366.mat"]

    mm = load_mat_model(filename)
    sm = convert(StandardModel, mm)
    mm2 = convert(MATModel, sm)

    @test Set(reactions(mm)) == Set(reactions(sm))
    @test Set(reactions(mm)) == Set(reactions(mm2))
end

@testset "MATModel generic interface" begin
    model = load_model(model_paths["e_coli_core.mat"])

    @test reaction_stoichiometry(model, "EX_ac_e") == Dict("ac_e" => -1)
    @test reaction_stoichiometry(model, 44) == Dict("ac_e" => -1)
end

ext: jl


text:
print("$(ARGS[1]) + $(ARGS[2]) = ")
println(parse(Int, ARGS[1]) + parse(Int, ARGS[2]))


text:
module Weno
export weno!,WenoData
using LinearAlgebra
#using Devectorize
#===================================================
Weno 5  method in 1d.
===================================================#
struct WenoData
    c::Array{Float64,2}
    dright::Array{Float64,1}
    dleft::Array{Float64,1}
    b0::Float64
    b1::Float64
    epsilon::Float64
    size::Int32
    InC::Array{Float64,1}
    right::Array{Float64,1}
    left::Array{Float64,1}
    reconstructed::Array{Float64,1}
    work::Array{Float64,1}
    alpharight::Array{Float64,1}
    alphaleft::Array{Float64,1}
    beta::Array{Float64,1}
    WenoData(size)=new(transpose(reshape(
        [11.0/6.0,-7.0/6.0,1.0/3.0,1.0/3.0,5.0/6.0,-1.0/6.0,
         -1.0/6.0,5.0/6.0,1.0/3.0,1.0/3.0,-7/6.0,11.0/6.0],3,4)),
               Array([3.0/10.0,3.0/5.0,1.0/10.0]),
               Array([1.0/10.0,3.0/5.0,3.0/10]),
               13.0/12.0,
               1.0/4.0,
                   1.0e-6,
                   size,
                   Array{Float64}(undef,size+4), #InC
                   Array{Float64}(undef,3),Array{Float64}(undef,3),#right,left
                   Array{Float64}(undef,2*size+8),#reconstructed
                   Array{Float64}(undef,size+2),#work
                   Array{Float64}(undef,3),Array{Float64}(undef,3),#alpha{right,left}
                   Array{Float64}(undef,3)#beta
               )

end

function weno!(W,F,L,In::Array{Float64},Out::Array{Float64})

    size=W.size

    h1 = -1.0/(L/size)

    # build an extended array with phantom cells to deal with periodicity:
    W.InC[1]=In[size-1]
    W.InC[2]=In[size]

    # same computing time (on my machine) with A and B:

    # A:
    @views W.InC[3:2+size]=In[:] #slower!
    # B:
    # @simd for i=1:size
    #    W.InC[2+i]=In[i]
    # end
    W.InC[size+3]=In[1]
    W.InC[size+4]=In[2]
    # precompute for regularity coefficients:
    @simd for vol=1:size+2
        W.work[vol]= (W.InC[vol]-2.0*W.InC[vol+1]+W.InC[vol+2])^2
    end
    # lets's start computation:
    @views for vol= 3:2+size
        # @simd for r= 0:2
        #     W.left[r+1] =dot(W.c[r+1,:],W.InC[vol-r:vol-r+2])
        #     W.right[r+1]=dot(W.c[r+2,:],W.InC[vol-r:vol-r+2])
        # end

        # faster than above, as it does not allocate memory:----
        for r=0:2
            W.left[r+1]=0.0
            W.right[r+1]=0.0
            k=vol-r-1
            @simd for i=1:3
                 W.left[r+1]+=W.c[r+1,i]*W.InC[k+i]
            end
            @simd for i=1:3
                W.right[r+1]+=W.c[r+2,i]*W.InC[k+i]
            end
        end
        #-------------------------------------------
        # regularity coefficients
        W.beta[1]=W.b0* W.work[vol]+
	W.b1*(3.0*W.InC[vol]-4.0*W.InC[vol+1]+W.InC[vol+2])^2

        W.beta[2]=W.b0* W.work[vol-1]+
        W.b1*(W.InC[vol-1]-W.InC[vol+1])^2

        W.beta[3]=W.b0*W.work[vol-2]+
	W.b1*(W.InC[vol-2]-4.0*W.InC[vol-1]+3*W.InC[vol])^2

        sleft=0.0
        sright=0.0
        @simd for r=1:3
            s=1.0/(W.epsilon+W.beta[r])^2
            W.alphaleft[r]=W.dleft[r]*s
            W.alpharight[r]=W.dright[r]*s
            sleft+=W.alphaleft[r]
            sright+=W.alpharight[r]
        end

        recleft=dot(W.alphaleft,W.left)
        recright=dot(W.alpharight,W.right)
        W.reconstructed[2*vol-1]  = recleft/sleft
        W.reconstructed[2*vol]= recright/sright
    end

    W.reconstructed[2*size+5:2*size+8]=W.reconstructed[5:8]
    #Numerical flux at boundaries of volumes:
    @simd for vol in 1:size
        W.work[vol+1]=F(W.reconstructed[2*vol+4], W.reconstructed[2*vol+5])
    end

    W.work[1]=W.work[size+1]
    # result:
    @simd for vol in 1:size
        Out[vol]=h1*(W.work[vol+1]-W.work[vol])
    end
    # slower:
    #@. Out[1:size] = h1*(W.work[2:size+1]-W.work[1:size])
    nothing
end

end

max_stars_repo_name: danielcort/benchmarks-python-julia-c

max_issues_repo_name: danielcort/benchmarks-python-julia-c

max_forks_repo_name: danielcort/benchmarks-python-julia-c


text:
@testset "MSA Stats" begin

    msa_types = (
        Matrix{Residue},
        NamedResidueMatrix{Array{Residue,2}},
        MultipleSequenceAlignment,
        AnnotatedMultipleSequenceAlignment
        )

    pf09645_sto = joinpath(pwd(), "data", "PF09645_full.stockholm")
    gaoetal2011 = joinpath(pwd(), "data", "Gaoetal2011.fasta")

    gaoetal_msas = [ read(gaoetal2011, FASTA, T) for T in msa_types ]
    pfam_msas    = [ read(pf09645_sto, Stockholm, T) for T in msa_types ]

    @testset "Gaps" begin

        for msa in gaoetal_msas
            @test gapfraction(msa) == 0.0
            @test gapfraction(getsequence(msa,1)) == 0.0
            @test gapfraction(msa[1,:]) == 0.0

            @test residuefraction(msa) == 1.0

            @test sum(coverage(msa)) == 6.0
        end

        for msa in pfam_msas
            @test gapfraction(getsequence(msa,4)) == 0.0
            @test gapfraction(msa[:,1]) == 0.75
            @test gapfraction(msa[:,1]) == 0.75

            @test residuefraction(getsequence(msa, 4)) == 1.0
            @test residuefraction(msa[:,1]) == 0.25

            @test coverage(msa)[4] == 1.0
        end

    end
end


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 = "oxigraph_server"
version = v"0.2.5"

# Collection of sources required to complete build
sources = [
    GitSource("https://github.com/oxigraph/oxigraph.git", "a21dcbb4f7355d7a00a86fbc5ad2c350a53629c4"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/oxigraph/server

cargo build --release --no-default-features --features=sled

install_license $WORKSPACE/srcdir/oxigraph/LICENSE-MIT

cp ../target/${rust_target}/release/oxigraph_server${exeext} ${bindir}/
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms(; experimental=true)

# Rust toolchain for i686 Windows is unusable
filter!(p -> !Sys.iswindows(p) || arch(p) != "i686", platforms)

# The products that we will ensure are always built
products = Product[
    ExecutableProduct("oxigraph_server", :oxigraph_server),
]

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

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies;
                compilers=[:c, :rust], preferred_gcc_version=v"7", lock_microarchitecture=false, julia_compat="1.6")

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
"""
    asfunc(r2)

Evolve `alpha_S(mu_R^2)`. Does not use `mu^2` grid or weight tables.

# Returns
- `alphas::FLoat64`: `alpha_S` value.
- `nf::Integer`: number of flavours at scale r2.
- `ierr::Integer`: error code
"""
function asfunc(r2::Float64)

    r2 = Ref{Float64}(r2)
    nf = Ref{Int32}()
    ierr = Ref{Int32}()

    alphas = @ccall asfunc_(r2::Ref{Float64}, nf::Ref{Int32},
                          ierr::Ref{Int32})::Float64

    alphas, nf[], ierr[]
end

"""
    evolfg(itype, func, def, iq0)

Evolve the flavour pdf set.

# Arguments
- `itype::Integer`: select un-polarised (1), polarised (2) or
time-like (3) evolution.
- `func`: User-defined function that returns input `x * f_j(x)`
at `iq0`. `j` is from `0` to `2 * nf`.
- `def::Array{Float64}`: input array containing the contribution of
quark species `i` to the input distribution `j`.
- `iq0::Integer`: grid index of the starting scale `mu_0^2`.

# Returns
- `epsi::Float64`: max deviation of the quadratic spline interpolation
from linear interpolation mid-between grid points.
"""
function evolfg(itype::Integer, func, def::Array{Float64}, iq0::Integer)

    itype = Ref{Int32}(itype)
    iq0 = Ref{Int32}(iq0)
    epsi = Ref{Float64}()

    @ccall evolfg_(itype::Ref{Int32}, func::Ptr{Cvoid}, def::Ref{Float64},
                   iq0::Ref{Int32}, epsi::Ref{Float64})::Nothing

    epsi[]
end

"""
    evsgns(itype, func, isns, n, iq0)

Evolve an arbitrary set of single/non-singlet pdfs. The
evolution can only run in FFNS or MFNS mode, as it is not
possible to correctly match at the thresholds as in evolfg.

# Arguments
The arguments are as for evolfg, expect def::Array{Float64}
is replaced with:
- `isns::Array{Int32,1}`: Input int array specifing the
evolution type. Entries can be (+1, -1, +-2) corresponding to
singlet, valence non-singlet and +/- q_ns singlets respectively.
- `n::Integer`: Number of singlet/non-singlet pdfs to evolve

# Returns
- `epsi::Float64`: Maximum deviation of the quadratic spline from
linear interpolation mid-between the grid points.
"""
function evsgns(itype::Integer, func, isns::Array{Int32,1}, n::Integer, iq0::Integer)

    itype = Ref{Int32}(itype)
    n = Ref{Int32}(n)
    iq0 = Ref{Int32}(iq0)
    epsi = Ref{Float64}()

    @ccall evsgns_(itype::Ref{Int32}, func::Ptr{Cvoid}, isns::Ref{Int32},
                   n::Ref{Int32}, iq0::Ref{Int32}, epsi::Ref{Float64})::Nothing

    epsi[]
end

"""
    evsgnsp(itype, func, isns, n, iq0)

Prototype parallel version on evsgns.
"""
function evsgnsp(itype::Integer, func, isns::Array{Int32,1}, n::Integer, iq0::Integer, jrun::Integer)

    itype = Ref{Int32}(itype)
    n = Ref{Int32}(n)
    iq0 = Ref{Int32}(iq0)
    jrun = Ref{Int32}(jrun)
    epsi = Ref{Float64}()

    @ccall evsgnsp_(itype::Ref{Int32}, func::Ptr{Cvoid}, isns::Ref{Int32},
                   n::Ref{Int32}, iq0::Ref{Int32}, epsi::Ref{Float64}, jrun::Ref{Int32})::Nothing

    epsi[]
end

"""
    usrpdf(fun, iset, n, offset)

Create a user-defined type-5 pdfset (same type as the output
of evsgns).

# Arguments
- `fun`: User-defined function with the signature
fun(ipdf::Integer, x::Float64, qq::Float64, first::UInt8)::Float64
specifying the values at x and qq of pdfset ipdf.
- `iset::Integer`: Pdfset identifier, between 1 and 24.
- `n::Integer`: Number of pdf tables in addition to gluon tables.
- `offset::Float64`: Relative offset at the thresholds mu_h^2, used
to catch matching discontinuities.

# Returns
- `epsi::Float64`: Maximum deviation of the quadratic spline from
linear interpolation mid-between the grid points.
"""
function usrpdf(fun, iset::Integer, n::Integer, offset::Float64)

    iset = Ref{Int32}(iset)
    n = Ref{Int32}(n)
    offset = Ref{Float64}(offset)
    epsi = Ref{Float64}()

    @ccall usrpdf_(fun::Ptr{Cvoid}, iset::Ref{Int32}, n::Ref{Int32},
        offset::Ref{Float64}, epsi::Ref{Float64})::Nothing

    epsi[]
end

"""
    nptabs(iset)

Get the number of pdf tables in set `iset::Integer`.
"""
function nptabs(iset::Integer)

    iset = Ref{Int32}(iset)

    ntabs = @ccall nptabs_(iset::Ref{Int32})::Int32

    ntabs[]
end

"""
    ievtype(iset)

Get the pdf evolution type for set `iset::Integer`.
"""
function ievtyp(iset::Integer)

    iset = Ref{Int32}(iset)

    ityp = @ccall ievtyp_(iset::Ref{Int32})::Int32

    ityp[]
end


text:
module VarianceComponentsHDFE

# using DataDeps
using CSV
using DataFrames, DataFramesMeta, Parameters
using LinearAlgebra, SparseArrays, Random, Statistics
using SparseArrays, LightGraphs, VectorizedRoutines, CSVFiles, DataFramesMeta
using Distributions, Arpack
using LinearOperators, FastClosures, Krylov
using ArgParse

include("init.jl")
include("leave_out_correction.jl")
include("parameters_settings.jl")
include("laplacians/Laplacians.jl")

using .Laplacians
include("solvers.jl")

export find_connected_set,prunning_connected_set,drop_single_obs, index_constr
export compute_movers, check_clustering, eff_res
export do_Pii, lincom_KSS, compute_matchid, leave_out_estimation, compute_whole
export Settings, JLAAlgorithm, ExactAlgorithm, AbstractLeverageAlgorithm

# Exporting these functions for ease of benchmarking/testing
export computeLDLinv, approxcholOperator, approxcholSolver

function parse_commandline()
    argparsesettings_obj = ArgParseSettings()

    @add_arg_table! argparsesettings_obj begin
        "path"
            help = "path to CSV file containing data"
            required = true
        "--first_id"
            help = "column index in CSV file for the first ID (e.g. Person).  Use the most granular type."
            arg_type = Int
            default = 1
        "--second_id"
            help = "column index in CSV file for the second ID (e.g. Firm).  Use the less granular type."
            arg_type = Int
            default = 2
        "--observation_id"
            help = "column index in CSV file for observation (e.g. Wage)."
            arg_type = Int
            default = 4
        "--algorithm"
            help = "type of algorithm: exact or JLA"
            arg_type = String
            default = "Exact"
        "--simulations"
            help = "number of simulations in the JLA algorithm"
            arg_type = Int
            default = 100
        "--header"
            help = "CSV file contains header"
            action = :store_true
        "--first_id_display"
            help = "The display text associated with first_id (e.g. Person)."
            arg_type = String
            default = "Person"
        "--second_id_display"
            help = "The display text associated with second_id (e.g. Firm)"
            arg_type = String
            default = "Firm"
        #=
        "--person_effects"
            help = "compute person effects"
            action = :store_true
        "--cov_effects"
            help = "compute cov effects"
            action = :store_true
        =#
        "--write_CSV"
            help = "write output to a CSV"
            action = :store_true
        "--output_path"
            help = "path to output CSV"
            arg_type = String
            default = "VarianceComponents.csv"
    end

    return parse_args(argparsesettings_obj)
end

function julia_main()::Cint
    try
        real_main()
    catch
        Base.invokelatest(Base.display_error, Base.catch_stack())
        return 1
    end
    return 0
end

function real_main()

    println("Number of threads: $(Threads.nthreads())")

    parsed_args = parse_commandline()

    path = parsed_args["path"]
    header = parsed_args["header"]
    first_idx = parsed_args["first_id"]
    second_idx = parsed_args["second_id"]
    observation_idx = parsed_args["observation_id"]
    algorithm = parsed_args["algorithm"]
    # person_effects = parsed_args["person_effects"]
    # cov_effects = parsed_args["cov_effects"]
    simulations = parsed_args["simulations"]
    first_id_display = parsed_args["first_id_display"]
    second_id_display = parsed_args["second_id_display"]

    data  = DataFrame!(CSV.File(path; header=header))
    id = data[:,first_idx]
    firmid = data[:,second_idx]
    y = data[:,observation_idx]

    controls = nothing

    if algorithm == "Exact"
        settings = Settings(leverage_algorithm = ExactAlgorithm(), person_effects=true, cov_effects=true)
    else
        settings = Settings(leverage_algorithm = JLAAlgorithm(num_simulations=simulations), person_effects=true, cov_effects=true)
    end

    θFE, θPE, θCOV, obs, β, Dalpha, Fpsi, Pii, Bii_pe, Bii_fe, Bii_cov = compute_whole(y,id,firmid,controls,settings;verbose=true)

    println("Bias-Corrected Variance Components:")
    println("Bias-Corrected Variance of $first_id_display: $θPE")
    println("Bias-Corrected Variance of $second_id_display: $θFE")
    println("Bias-Corrected Covariance of $first_id_display-$second_id_display Effects: $θCOV")

    if parsed_args["write_CSV"]

        y_output = y[obs]
        id_output = id[obs]
        firmid_output = firmid[obs]

        max_length = length(obs)

        output = DataFrame(observation = obs,
                           worker_id = id_output,
                           firm_id = firmid_output,
                           y = y_output,
                           beta = vcat(β,missings(max(max_length-length(β),0))),
                           D_alpha = Dalpha,
                           F_psi = Fpsi,
                           Pii = Pii,
                           Bii_pe = Bii_pe,
                           Bii_fe = Bii_fe,
                           Bii_cov = Bii_cov,
                           variance_comp_firm_effects = [θFE; missings(max_length-1)],
                           variance_comp_person_effects = [θPE; missings(max_length-1)],
                           covariance_comp_effects = [θCOV; missings(max_length-1)])
        output_path = parsed_args["output_path"]
        CSV.write(output_path,output)
    end
    return
end

end


text:
### A Pluto.jl notebook ###
# v0.12.21

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ 3f8349ba-8be8-11eb-32b9-55388c7242aa
using Plots, Colors, PlutoUI, StatsBase

# ╔═╡ 5d3d4988-8be8-11eb-1de8-3b114233e526
function simulate(N,p)
	v = fill(0,N,N)
	time = 0
	while any( v .== 0 ) && time<100
		time += 1
		for i= 1:N, j=1:N
		     if rand()<p && v[i,j]==0
				   v[i,j] = time
			 end
		end
	end
	v
end

# ╔═╡ 5d3d7b56-8be8-11eb-1bb1-ddecbefefc49
md"""
N= $(@bind N Slider(2:20, show_value=true, default=8))
"""

# ╔═╡ 5d4409e4-8be8-11eb-2d06-03e4aa311fc0
md"""
p= $(@bind prob Slider(0.01:.01:1, show_value=true, default=.1))

t = $(@bind tt Slider(1:100, show_value=true, default=1))
"""

# ╔═╡ 5d414452-8be8-11eb-233e-2d81eaacbfb6
simulation = simulate(N, prob)

# ╔═╡ 82299e50-8bec-11eb-3591-6b02b1b2a7de
begin

	plot()
	cdf= [ count(simulation .≤ i) for i=0:100]
	bar(cdf, c=:purple, legend=false, xlim=(0,tt),alpha=0.8)
end

# ╔═╡ 075669c8-8bef-11eb-288e-791816cc0d5b
begin
	newcdf = [ count(simulation .> i) for i=0:100]
	bar!( newcdf, c=RGB(0,1,0), legend=false, xlim=(0,tt),alpha=0.8)
end

# ╔═╡ 482ebeb0-8bec-11eb-2aaf-b5522fa606d4
bar(countmap(simulation[:]), c=:red, legend=false, xlim=(0,tt+.5))

# ╔═╡ c1e720f4-8bf5-11eb-386b-b32d313a2996
@bind bernoulliwidth Slider(10:10:500, show_value=true)

# ╔═╡ 3b02051e-8bf4-11eb-011b-3b7131b245a6
HTML("""<img src="https://upload.wikimedia.org/wikipedia/commons/thumb/b/b3/ETH-BIB-Bernoulli%2C_Daniel_%281700-1782%29-Portrait-Portr_10971.tif_%28cropped%29.jpg/440px-ETH-BIB-Bernoulli%2C_Daniel_%281700-1782%29-Portrait-Portr_10971.tif_%28cropped%29.jpg" alt="" width=$(bernoulliwidth))>""")

# ╔═╡ 5f5a0caa-8bf6-11eb-1242-a91551de2922
begin
	b = "hello"
	HTML(b)
end

# ╔═╡ 5d460802-8be8-11eb-164e-71074e4e4b66
begin
	rectangle(w, h, x, y) = Shape(x .+ [0,w,w,0], y .+ [0,0,h,h])
	circle(r,x,y) = (θ = LinRange(0,2π,500); (x.+r.*cos.(θ), y.+r.*sin.(θ)))
end

# ╔═╡ 5d44663c-8be8-11eb-0986-bfc7546ee2ab
begin
	w = .9
	h = .9
	c = [RGB(0,1,0),RGB(1,0,0),:purple][1 .+ (simulation .< tt) .+ (simulation .<  (tt.-1))]

	plot()

	for i=1:N, j=1:N
		plot!( rectangle(w,h, i, j), c=:black)
		plot!( circle(.3,i+.45,j+.45), c=c[i,j], fill=true)
	end

	for i=1:N, j=1:N
		if simulation[i,j] < tt
	       annotate!(i+.45,j+.5,text("$(simulation[i,j])",font(7),:white))
		end
	end

	plot!(legend=false,ratio=1, axis=false, lims=(0.5,N+1.1), ticks=false,title="time = $(tt-1)  red count: $(sum(simulation.<tt))")

end

# ╔═╡ Cell order:
# ╠═3f8349ba-8be8-11eb-32b9-55388c7242aa
# ╠═5d3d4988-8be8-11eb-1de8-3b114233e526
# ╠═5d3d7b56-8be8-11eb-1bb1-ddecbefefc49
# ╠═5d414452-8be8-11eb-233e-2d81eaacbfb6
# ╟─5d4409e4-8be8-11eb-2d06-03e4aa311fc0
# ╠═5d44663c-8be8-11eb-0986-bfc7546ee2ab
# ╠═82299e50-8bec-11eb-3591-6b02b1b2a7de
# ╠═075669c8-8bef-11eb-288e-791816cc0d5b
# ╠═482ebeb0-8bec-11eb-2aaf-b5522fa606d4
# ╠═c1e720f4-8bf5-11eb-386b-b32d313a2996
# ╠═3b02051e-8bf4-11eb-011b-3b7131b245a6
# ╠═5f5a0caa-8bf6-11eb-1242-a91551de2922
# ╠═5d460802-8be8-11eb-164e-71074e4e4b66


text:
@defcomp co2forcing begin
    c_CO2concentration = Parameter(index=[time], unit="ppbv")
    f0_CO2baseforcing = Parameter(unit="W/m2", default=1.68)
    fslope_CO2forcingslope = Parameter(unit="W/m2", default=5.5)
    c0_baseCO2conc = Parameter(unit="ppbv", default=400859.5833)
    f_CO2forcing = Variable(index=[time], unit="W/m2")

    function run_timestep(p, v, d, t)

        # eq.13 in Hope 2006
        # the max() condition was added to prevent numerical issues due to  non-positive CO2 concentration levels which occurs in 2300 for low-emission scenarios
        # if the growth effects are paired with GDP-emission feedback
        if p.c_CO2concentration[t] > 0
            v.f_CO2forcing[t] = max(0, p.f0_CO2baseforcing + p.fslope_CO2forcingslope * log(p.c_CO2concentration[t] / p.c0_baseCO2conc) )
        else
            v.f_CO2forcing[t] = 0
        end
    end
end


text:
using ExtremeStats
println("Reading data")
@time x,lons,lats,nlon,nlat,NpY,N_years = load_X("/Net/Groups/BGI/people/mjung/FLUXCOM/_ENSEMBLES/8daily/TERall/","TERall_","TERall",2001:2012,[20,40],[40, 20]);
println("Calculating anomalies")
gc()
@time x=get_anomalies(x,NpY,nlon,nlat)
println("Labelling extremes")
gc()
@time el=label_Extremes(x,0.05);
println("Calculating time bounds")
gc()
@time getTbounds(el);
println("Calculating features")
@time f1=getFeatures(el,Features.Mean,Features.Max_z,Features.Min_z,Features.Duration,Features.Size,Features.NumPixel,Features.Quantile);
gc()
@time f2=getFeatures(el,Features.TS_ZValue,Features.TS_Area)

max_stars_repo_name: CAB-LAB/ExtremeStats

max_issues_repo_name: CAB-LAB/ExtremeStats

max_forks_repo_name: CAB-LAB/ExtremeStats


text:
using Makie, CairoMakie, Plots
using NotebookToLaTeX
import NotebookToLaTeX:dispatch_output
import NotebookToLaTeX:createfolders
import NotebookToLaTeX:createtemplate
import NotebookToLaTeX:extractnotebook
import NotebookToLaTeX:collectoutputs
using Test

@testset "NotebookToLatex.jl" begin
    # Write your tests here.
    path = "./project/"

    @testset "Creating Folders and Files" begin
        createfolders(path)
        @test isdir(path)
        @test isdir(path * "/notebooks")
        @test isdir(path * "/figures")
        @test isdir(path * "/frontmatter")

        createtemplate(path, :book)
        @test isfile(path * "/main.tex")
    end

    #= @testset "Parsing Pluto Notebook" begin =#
    #=     nb = extractnotebook("./notebooktest.jl") =#
    #=     outputs  = collectoutputs(nb, path) =#
    #=     println(outputs) =#
    #=     @test outputs[1] == (:nothing, "") =#
    #=     @test outputs[2] == (:plot, "notebooktest_figure1.png") =#
    #=     @test isfile(path * "/figures/notebooktest_figure1.png") =#
    #= end =#

    @testset "Pluto to Latex" begin
        plutopath = "./build_latex/"
        fontpath = "/home/davibarreira/.local/share/fonts/Unknown Vendor/TrueType/JuliaMono/"
        notebooktolatex("./pluto/notebooktest.jl", template=:book, fontpath=fontpath)
        @test isfile(plutopath * "main.tex")
        @test isfile(plutopath * "preface.tex")
        @test isfile(plutopath * "frontmatter/copyright.tex")
        @test isfile(plutopath * "frontmatter/titlepage.tex")
        @test isfile(plutopath * "julia_font.tex")
        @test isfile(plutopath * "julia_listings.tex")
        @test isfile(plutopath * "julia_listings_unicode.tex")
        @test isfile(plutopath * "/notebooks/notebooktest.tex")
        @test !isfile(plutopath * "/fonts/JuliaMono_Bold.ttf")
        @test !isfile(plutopath * "/fonts/JuliaMono_Medium.ttf")
        @test !isfile(plutopath * "/fonts/JuliaMono_Regular.ttf")
        @test isfile(plutopath * "/figures/figure.pdf")
        @test isfile(plutopath * "/figures/notebooktest_figure1.png")
        @test isfile(plutopath * "/figures/notebooktest_figure2.pdf")
        @test isfile(plutopath * "/figures/plotexample.png")

        notebooktolatex("./pluto/foo.jl", template=:book, fontpath=fontpath)
        notebooktolatex("./pluto/foo2.jl", template=:book, fontpath=fontpath)
        @test isfile(plutopath * "/notebooks/foo.tex")
        @test isfile(plutopath * "/notebooks/foo2.tex")
    end

    @testset "Jupyter to Latex" begin
        jupyterpath = "./jupyter/build_latex/"
        notebooktolatex("./jupyter/jupyternotebook.ipynb", jupyterpath, template=:mathbook)
        @test isfile(jupyterpath * "main.tex")
        @test isfile(jupyterpath * "preface.tex")
        @test isfile(jupyterpath * "frontmatter/copyright.tex")
        @test isfile(jupyterpath * "frontmatter/titlepage.tex")
        @test isfile(jupyterpath * "julia_font.tex")
        @test isfile(jupyterpath * "julia_listings.tex")
        @test isfile(jupyterpath * "julia_listings_unicode.tex")
        @test isfile(jupyterpath * "/notebooks/jupyternotebook.tex")
        @test isfile(jupyterpath * "/fonts/JuliaMono_Bold.ttf")
        @test isfile(jupyterpath * "/fonts/JuliaMono_Medium.ttf")
        @test isfile(jupyterpath * "/fonts/JuliaMono_Regular.ttf")
        @test isfile(jupyterpath * "/figures/figure.pdf")
        @test isfile(jupyterpath * "/figures/jupyternotebook_figure1.pdf")
        @test isfile(jupyterpath * "/figures/jupyternotebook_figure1.svg")
        @test isfile(jupyterpath * "/figures/jupyternotebook_figure2.png")
        @test isfile(jupyterpath * "/figures/plotexample.png")

        rm(jupyterpath, recursive=true)
    end

    @testset "Template :article and :matharticle" begin
        jupyterpath = "./jupyter/build_latex/"

        for template in [:article, :matharticle]
            notebooktolatex("./jupyter/jupyternotebook.ipynb", jupyterpath, template=template)
            @test isfile(jupyterpath * "main.tex")
            @test !isfile(jupyterpath * "preface.tex")
            @test !isfile(jupyterpath * "frontmatter/copyright.tex")
            @test !isfile(jupyterpath * "frontmatter/titlepage.tex")
            @test isfile(jupyterpath * "julia_font.tex")
            @test isfile(jupyterpath * "julia_listings.tex")
            @test isfile(jupyterpath * "julia_listings_unicode.tex")
            @test isfile(jupyterpath * "/notebooks/jupyternotebook.tex")
            @test isfile(jupyterpath * "/fonts/JuliaMono_Bold.ttf")
            @test isfile(jupyterpath * "/fonts/JuliaMono_Medium.ttf")
            @test isfile(jupyterpath * "/fonts/JuliaMono_Regular.ttf")
            @test isfile(jupyterpath * "/figures/figure.pdf")
            @test isfile(jupyterpath * "/figures/jupyternotebook_figure1.pdf")
            @test isfile(jupyterpath * "/figures/jupyternotebook_figure1.svg")
            @test isfile(jupyterpath * "/figures/jupyternotebook_figure2.png")
            @test isfile(jupyterpath * "/figures/plotexample.png")
            @test contains(read(jupyterpath * "/notebooks/jupyternotebook.tex", String), "\\section{My Notebook}")
            @test contains(read(jupyterpath * "/notebooks/jupyternotebook.tex", String), "\\subsection{Start importing}")
            @test contains(read(jupyterpath * "/notebooks/jupyternotebook.tex", String), "\\subsubsection{Some theory}")

            @test contains(read(jupyterpath * "/notebooks/jupyternotebook.tex", String), "\\subsubsection{Some theory}")

            if template == :matharticle
                mathspecific = ["\\usepackage{amsfonts, amsthm,amsmath,amssymb,mathtools}",
                    "\\usepackage{bbm}",
                    "\\usepackage{bm}",
                    "\\usepackage{mathtools}",
                    "\\usepackage{thmtools}"]
                for mathpackage in mathspecific
                    @test contains(read(jupyterpath * "main.tex", String), mathpackage)
                end
            end

            rm(jupyterpath, recursive=true)
        end
    end

        rm(path, recursive=true)
    rm("./build_latex/", recursive=true)
end


text:
# Types for managing event selection
abstract type BitField end

fields(x::BitField)     = fields(typeof(x))
name(x::BitField)       = name(typeof(x))
defaults(x::BitField)   = defaults(typeof(x))
value(x::BitField)      = x.val
Base.convert(::Type{T}, x::T) where {T <: BitField} = x

Base.:(==)(x::BitField, y::BitField) = value(x) == value(y)
Base.hash(x::BitField, u::UInt = UInt(0x12038497)) = hash(value(x), u)

(::Type{T})(;kw...) where {T <: BitField} = construct(T, (;kw...))

# Main entry point
function construct(::Type{T}, kw::NamedTuple) where {T <: BitField}
    # Apply defaults to the provided keywords.
    kw = merge(defaults(T), kw)
    val = zero(UInt64)
    fieldpairs = fields(T)
    for (k, v) in pairs(kw)
        # Get the start and stop points of this bit field.
        start, stop = fieldpairs[k]

        # Mask the provided value to the requested number of bits.
        v = v & mask(0, stop - start)

        # Set the corresponding bits in `x`.
        # Since `x` is initialized to zero, we do not need to explicitly clear these bits.
        val |= (v << start)
    end
    return T(val)
end

function Base.show(io::IO, reg::BitField)
    print(io, name(reg), ": ")
    for k in keys(fields(reg))
        print(io, " $k=$(string(reg[k]; base = 16))")
    end
    println(io)
    return nothing
end
hex(i::BitField) = hex(value(i))

# Allow this to be written to an IO by forwarding to the wrapped value
Base.write(io::IO, reg::BitField) = write(io, value(reg))
Base.convert(::Type{T}, reg::BitField) where {T <: Integer} = convert(T, value(reg))

# For convenience purposes
function Base.getindex(reg::BitField, name::Symbol)
    start, stop = fields(reg)[name]
    return (value(reg) & mask(start, stop)) >> start
end

### Core Event Select
"""
    CoreSelectRegister(; kw...)

Construct a bitmask for programming `Core` level counters.

Keywords
========
- `event::UInt`: Select the event to be counted. Default: `0x00`
- `umask::Uint`: Select the subevent to be counted within the selected event. Default: `0x00`
- `usr::Bool`: Specifies the counter should be active when the processor is operating at
    privilege modes 1, 2, and 3. Default: `true`.
- `os::Bool`: Specifies the counter should be active when the processor is operating at
    privilege mode 0. Default: `true`.
- `e::Bool`: Edge detect. Default: `false`.
- `en::Bool`: Enable the counter. Default: `true`.
- `inv::Bool`: When set, inverts the counter-mask (CMASK) comparison, so that both greater
    than or equal to and less than comparisons can be made (0: greater than or equal; 1:
    less than). Note if counter-mask is programmed to zero, INV flag is ignored. Default: `false`.
- `cmask::Bool`: When this field is not zero, a logical processor compares this mask to the
    events count of the detected microarchitectural condition during a single cycle. If the
    event count is greater than or equal to this mask, the counter is incremented by one.
    Otherwise the counter is not incremented.

    This mask is intended for software to characterize microarchitectural conditions that
    can count multiple occurrences per cycle (for example, two or more instructions retired
    per clock; or bus queue occupations). If the counter-mask field is 0, then the counter
    is incremented each cycle by the event count associated with multiple occurrences.
"""
struct CoreSelectRegister <: BitField
    val::UInt64
end

function fields(::Type{CoreSelectRegister})
    nt = (
        event = (0, 7),
        umask = (8, 15),
        usr   = (16, 16),
        os    = (17, 17),
        e     = (18, 18),
        #pc    = (19, 19), # not implemented on Skylake+
        #int   = (20, 20), # not worrying about this
        en    = (22, 22),
        inv   = (23, 23),
        cmask = (24, 31),
    )
    return nt
end
name(::Type{CoreSelectRegister}) = "Core Event Select Register"
defaults(::Type{CoreSelectRegister}) = (usr = true, os = true, en = true)

### Uncore Select
"""
    UncoreSelectRegister(; kw...)

Construct a bitmask for programming `Uncore` level counters.

Keywords
========

- `event::UInt`: Select event to be counted. Default: `0x00`
- `umask::UInt`: Select subevents to be counted within the selected event. Default: `0x00`
- `reset::Bool`: When set to 1, the corresponding counter will be cleared to 0. Default: `false`
- `edge_detact::Bool`: When set to 1, rather than measuring the event in each cycle it is
    active, the corresponding counter will increment when a 0 to 1 transition (i.e. rising edge)
    is detected.

    When 0, the counter will increment in each cycle that the event is asserted.

    NOTE: `edge_detect` is in series following `thresh`, Due to this, the `thresh` field
    must be set to a non-0 value. For events that increment by no more than 1 per cycle,
    set `thresh` to 0x1. Default: `false`.
- `overflow_enable::Bool`: When this bit is set to 1 and the corresponding counter overflows,
    an overflow message is sent to the UBox’s global logic. The message identifies the unit
    that sent it. Default: `false`.
- `en::Bool`: Local Counter Enable. Default: `true`.
- `invert::Bool`: Invert comparison against Threshold.

    0 - comparison will be ‘is event increment >= threshold?’.

    1 - comparison is inverted - ‘is event increment < threshold?’

    e.g. for a 64-entry queue, if SW wanted to know how many cycles the queue had fewer
    than 4 entries, SW should set the threshold to 4 and set the invert bit to 1.
    Default: `false`.
- `thresh::UInt`: Threshold is used, along with the invert bit, to compare against the
    counter’s incoming increment value. i.e. the value that will be added to the counter.

    For events that increment by more than 1 per cycle, if the threshold is set to a value
    greater than 1, the data register will accumulate instances in which the event
    increment is >= threshold. Default: `0x00`.
"""
struct UncoreSelectRegister <: BitField
    val::UInt64
end

function fields(::Type{UncoreSelectRegister})
    nt = (
        event           = (0, 7),
        umask           = (8, 15),
        reset           = (17, 17),
        edge_detect     = (18, 18),
        tid_enable      = (19, 19),      # Only applicable for CHAS
        overflow_enable = (20, 20),
        en              = (22, 22),
        invert          = (23, 23),
        thresh          = (24, 31),
    )
    return nt
end
name(::Type{UncoreSelectRegister}) = "Uncore Event Select Register"
defaults(::Type{UncoreSelectRegister}) = (en = true,)

### CHA Filter 0
"""
    CHAFilter0(; kw...)
"""
struct CHAFilter0 <: BitField
    val::UInt64
end
fields(::Type{CHAFilter0}) = (
    thread_id   = (0, 2),
    core_id     = (3, 8),
    LLC_I_state = (17, 17),
    SF_S_state  = (18, 18),
    SF_E_state  = (19, 19),
    SF_H_state  = (20, 20),
    LLC_S_state = (21, 21),
    LLC_E_state = (22, 22),
    LLC_M_state = (23, 23),
    LLC_F_state = (24, 24),
)
name(::Type{CHAFilter0}) = "CHA Filter 0"
defaults(::Type{CHAFilter0}) = (
    LLC_I_state = true,
    LLC_S_state = true,
    LLC_E_state = true,
    LLC_M_state = true,
    LLC_F_state = true,
)

# ### CHA Filter 1
"""
    CHAFilter1(; kw...)
"""
struct CHAFilter1 <: BitField
    val::UInt64
end

fields(::Type{CHAFilter1}) = (
    remote          = (0, 0),
    loc             = (1, 1),
    all_opc         = (3, 3),
    near_memory     = (4, 4),
    not_near_memory = (5, 5),
    opc0            = (9, 18),
    opc1            = (19, 28),
    non_coherent    = (30, 30),
    isoc            = (31, 31),
)
name(::Type{CHAFilter1}) = "CHA Filter 1"
defaults(::Type{CHAFilter1}) = (
    remote          = true,
    loc             = true,
    all_opc         = true,
    near_memory     = true,
    not_near_memory = true
)


text:
########################################
# Incoming messages processing
########################################

struct Message
    stream::String
    topic::String
    type::String
    sender_id::Int
    content::String
end
StructTypes.StructType(::Type{Message}) = StructTypes.OrderedStruct()

struct TimedMessage
    id::Int
    createts::Int
    exects::Int
    msg::Message
end
StructTypes.StructType(::Type{TimedMessage}) = StructTypes.OrderedStruct()
TimedMessage(createts, exects, msg) = TimedMessage(-1, createts, exects, msg)
TimedMessage(createts::DateTime, exects::DateTime, msg) = TimedMessage(-1, toepoch(createts), toepoch(exects), msg)
TimedMessage(id, createts::DateTime, exects::DateTime, msg) = TimedMessage(id, toepoch(createts), toepoch(exects), msg)

tablename(::Type{TimedMessage}) = "messages"
idproperty(::Type{TimedMessage}) = :id
autoincrement(::Type{TimedMessage}) = :id

########################################
# TimeZone information
########################################

struct Sender
    id::Int
    tz::String
end
StructTypes.StructType(::Type{Sender}) = StructTypes.OrderedStruct()
tablename(::Type{Sender}) = "senders"
idproperty(::Type{Sender}) = :id


text:
```

ДАНО: Робот - в произвольной клетке поля

РЕЗУЛЬТАТ: Робот - в исходном положении, и все клетки по периметру внешней рамки промакированы

```

```

1. Предлагаю сделать нового робота у которого функция move! = try_move!. Т.е. каждым шагом он будет пытаться обойти что-то.
2. Я не помню как объявлять функцию заранее для try_move_var2!
UPD. 1. оказалость не нужным
UPD. Для кошмарных внутренних перегородок можно юзать робота Beta и его функцию get_back!, а для
обхода периметра robota.beta = originRobot

```
using HorizonSideRobots
originRobot = Robot(animate = true,"15example.sit")
include("../AbstractType.jl")
include("../functions.jl")

function master(robot::Robot)
    back_path = moving_in_angle!(robot)
    putmarker!(robot)
    for side in [Nord, Ost, Sud, West]
        flag = true
        while (flag)
            if !(isborder(robot, side))
                move!(robot, side)
                putmarker!(robot)
            else
                if (try_move!(robot, side))
                    putmarker!(robot)
                else
                    flag = false
                end
            end
        end
    end
    moving_back_to_start!(robot, back_path)
end

master!(originRobot)

max_stars_repo_name: Igetsuga/AlgorithmicsHomework

max_issues_repo_name: Igetsuga/AlgorithmicsHomework

max_forks_repo_name: Igetsuga/AlgorithmicsHomework


text:
add_together(x,y) = x + y

add_together(1)

max_stars_repo_name: icankeep/jupyterlab-lsp

max_issues_repo_name: icankeep/jupyterlab-lsp

max_forks_repo_name: icankeep/jupyterlab-lsp


text:
"radon transform of LoG"
# $$e_{\sigma}(t)=\left(\frac{1}{\sqrt{2 \pi} \sigma^{5}}\right)\left(t^{2}-\sigma^{2}\right) \exp ^{-t^{2} / 2 \sigma^{2}}$$
# note that we consider 3D

"Compute Radon transform of normalized Laplacian of Gaussian in 3D"
function radon_log(q, sigma; z0=0)
    _radon_log(t, ang; z0=0) = (t.^2 .- sigma^2 .+ z0^2) .* exp.(-(t.^2 + z0^2 ) / ( 2*sigma^2 ))

    nangles, detcount = size(q)
    p_log = zeros(nangles, detcount)
    tt = (-detcount/2+0.5):1.0:(detcount/2-0.5)

    for ang in 1:nangles
        for (i, t) in enumerate(tt)
            dotp = q[ang, :] .* _radon_log.(t .- tt, sigma, z0=z0) ./ ( sqrt(2*pi) * sigma^3 )

            p_log[ang, i] = sum(dotp)
        end
    end

    return p_log
end

"""
    radon_filter(p, angles, fun_filter, filt_sz)

Radon transform of a filter slice by slice specified by a 1D fucntion `fun_filter(theta,angle)`
"""
function radon_filter(p::Array{T, 3}, angles::Array{T, 1}, fun_filter, szfilter::Int) where {T<:AbstractFloat}
    if !isodd(szfilter)
        println("! error: only support odd number of filter size, szfilter is increased by 1")
        szfilter += 1
    end

    nangles, nslice, detcount = size(p)

    halfsz = div(szfilter-1, 2)
    tt = collect(-halfsz:1:halfsz)
    out = zeros(size(p))
    filt_ = zeros(2*halfsz+1)

    for ang=1:nangles
        filt_ .= fun_filter.(tt, angles[ang])
        for slice=1:nslice
            for i = 1:2*halfsz+1
                let i=i # needed for simd
                    fi = filt_[i]
                    @simd for j = halfsz+1:(detcount-halfsz)
                        @inbounds out[ang, slice, i+j-(halfsz+1)] += p[ang, slice, j] * fi
                    end
                end
            end
        end
    end
    return out
end

max_stars_repo_name: JuliaTomo/DynamicTomographyWithFlow

max_issues_repo_name: JuliaTomo/DynamicTomographyWithFlow

max_forks_repo_name: JuliaTomo/DynamicTomographyWithFlow


text:
using DiffEqFlux, DelayDiffEq, OrdinaryDiffEq, StochasticDiffEq, Test, Random

mp = Float32[0.1,0.1]
x = Float32[2.; 0.]
xs = Float32.(hcat([0.; 0.], [1.; 0.], [2.; 0.]))
tspan = (0.0f0,1.0f0)
dudt = Flux.Chain(Flux.Dense(2,50,tanh),Flux.Dense(50,2))
fastdudt = FastChain(FastDense(2,50,tanh),FastDense(50,2))
fastcdudt = FastChain(FastDense(2,50,tanh,precache=true,numcols=size(xs)[2]),FastDense(50,2,precache=true,numcols=size(xs)[2]))
staticdudt = FastChain(StaticDense(2,50,tanh),StaticDense(50,2))

NeuralODE(dudt,tspan,Tsit5(),save_everystep=false,save_start=false)(x)
NeuralODE(dudt,tspan,Tsit5(),saveat=0.1)(x)
NeuralODE(dudt,tspan,Tsit5(),saveat=0.1,sensealg=TrackerAdjoint())(x)

NeuralODE(dudt,tspan,Tsit5(),save_everystep=false,save_start=false)(xs)
NeuralODE(dudt,tspan,Tsit5(),saveat=0.1)(xs)
NeuralODE(dudt,tspan,Tsit5(),saveat=0.1,sensealg=TrackerAdjoint())(xs)

@info "Test some gradients"

node = NeuralODE(dudt,tspan,Tsit5(),save_everystep=false,save_start=false)
grads = Zygote.gradient(()->sum(node(x)),Flux.params(x,node))
@test ! iszero(grads[x])
@test ! iszero(grads[node.p])

grads = Zygote.gradient(()->sum(node(xs)),Flux.params(xs,node))
@test ! iszero(grads[xs])
@test ! iszero(grads[node.p])

node = NeuralODE(dudt,tspan,Tsit5(),save_everystep=false,save_start=false,sensealg=TrackerAdjoint())
grads = Zygote.gradient(()->sum(node(x)),Flux.params(x,node))
@test ! iszero(grads[x])
@test ! iszero(grads[node.p])

grads = Zygote.gradient(()->sum(node(xs)),Flux.params(xs,node))
@test ! iszero(grads[xs])
@test ! iszero(grads[node.p])

node = NeuralODE(dudt,tspan,Tsit5(),save_everystep=false,save_start=false,sensealg=BacksolveAdjoint())
grads = Zygote.gradient(()->sum(node(x)),Flux.params(x,node))
@test ! iszero(grads[x])
@test ! iszero(grads[node.p])

grads = Zygote.gradient(()->sum(node(xs)),Flux.params(xs,node))
@test ! iszero(grads[xs])
@test ! iszero(grads[node.p])

node = NeuralODE(dudt,tspan,Tsit5(),save_everystep=false,save_start=false,sensealg=BacksolveAdjoint(autojacvec=ZygoteVJP()))
grads = Zygote.gradient(()->sum(node(x)),Flux.params(x,node))
@test ! iszero(grads[x])
@test ! iszero(grads[node.p])

grads = Zygote.gradient(()->sum(node(xs)),Flux.params(xs,node))
@test ! iszero(grads[xs])
@test ! iszero(grads[node.p])

## Fast

@info "Test some fast layers"

node = NeuralODE(fastdudt,tspan,Tsit5(),save_everystep=false,save_start=false)
pd = Flux.params(node)[1]
gradsnc = Zygote.gradient(()->sum(node(x)),Flux.params(x,node)) # no cache
@test ! iszero(gradsnc[x])
@test ! iszero(gradsnc[node.p])

gradsnc2 = Zygote.gradient(()->sum(node(xs)),Flux.params(xs,node))
@test ! iszero(gradsnc2[xs])
@test ! iszero(gradsnc2[node.p])

nodec = NeuralODE(fastcdudt,tspan,Tsit5(),save_everystep=false,save_start=false,p=pd)
gradsc = Zygote.gradient(()->sum(nodec(x)),Flux.params(x,nodec)) #with cache
@test ! iszero(gradsc[x])
@test ! iszero(gradsc[nodec.p])
@test gradsnc[x] ≈ gradsc[x] rtol=1e-6
@test gradsnc[node.p] ≈ gradsc[nodec.p] rtol=1e-6

gradsc2 = Zygote.gradient(()->sum(nodec(xs)),Flux.params(xs,nodec))
@test ! iszero(gradsc2[xs])
@test ! iszero(gradsc2[nodec.p])
@test gradsnc2[xs] ≈ gradsc2[xs] rtol=1e-6
@test gradsnc2[node.p] ≈ gradsc2[nodec.p] rtol=1e-6
#test with low tolerance ode solver
node = NeuralODE(fastdudt, tspan, Tsit5(), abstol=1e-12, reltol=1e-12, save_everystep=false, save_start=false)
pd = Flux.params(node)[1]
gradsnc = Zygote.gradient(()->sum(node(x)),Flux.params(x,node))
@test ! iszero(gradsnc[x])
@test ! iszero(gradsnc[node.p])

nodec = NeuralODE(fastcdudt, tspan, Tsit5(), abstol=1e-12, reltol=1e-12, save_everystep=false, save_start=false,p=pd)
gradsc = Zygote.gradient(()->sum(nodec(x)),Flux.params(x,nodec))
@test gradsnc[x] ≈ gradsc[x] rtol=1e-3
@test gradsnc[node.p] ≈ gradsc[nodec.p] rtol=1e-3

node = NeuralODE(fastdudt,tspan,Tsit5(),save_everystep=false,save_start=false,sensealg=TrackerAdjoint())
grads = Zygote.gradient(()->sum(node(x)),Flux.params(x,node))
@test ! iszero(grads[x])
@test ! iszero(grads[node.p])

grads = Zygote.gradient(()->sum(node(xs)),Flux.params(xs,node))
@test ! iszero(grads[xs])
@test ! iszero(grads[node.p])

goodgrad = grads[node.p]
p = node.p

node = NeuralODE(fastcdudt,tspan,Tsit5(),save_everystep=false,save_start=false,sensealg=TrackerAdjoint())
grads = Zygote.gradient(()->sum(node(x)),Flux.params(x,node))
@test ! iszero(grads[x])
@test ! iszero(grads[node.p])

grads = Zygote.gradient(()->sum(node(xs)),Flux.params(xs,node))
@test ! iszero(grads[xs])
@test ! iszero(grads[node.p])

goodgradc = grads[node.p]
pc = node.p

node = NeuralODE(fastdudt,tspan,Tsit5(),save_everystep=false,save_start=false,sensealg=BacksolveAdjoint(),p=p)
grads = Zygote.gradient(()->sum(node(x)),Flux.params(x,node))
@test ! iszero(grads[x])
@test ! iszero(grads[node.p])

grads = Zygote.gradient(()->sum(node(xs)),Flux.params(xs,node))
@test ! iszero(grads[xs])
@test ! iszero(grads[node.p])
goodgrad2 = grads[node.p]
@test goodgrad ≈ goodgrad2 # Make sure adjoint overloads are correct

node = NeuralODE(fastcdudt,tspan,Tsit5(),save_everystep=false,save_start=false, sensealg=BacksolveAdjoint(),p=pc)
grads = Zygote.gradient(()->sum(node(x)),Flux.params(x,node))
@test ! iszero(grads[x])
@test ! iszero(grads[node.p])

grads = Zygote.gradient(()->sum(node(xs)),Flux.params(xs,node))
@test !iszero(grads[xs])
@test ! iszero(grads[node.p])
goodgrad2 = grads[node.p]
@test goodgradc ≈ goodgrad2 rtol=1e-6

node = NeuralODE(staticdudt,tspan,Tsit5(),save_everystep=false,save_start=false,sensealg=BacksolveAdjoint(autojacvec=false),p=p)
grads = Zygote.gradient(()->sum(node(x)),Flux.params(x,node))
@test ! iszero(grads[x])
@test ! iszero(grads[node.p])

grads = Zygote.gradient(()->sum(node(xs)),Flux.params(xs,node))
goodgrad2 = grads[node.p]
@test goodgrad ≈ goodgrad2 rtol = 1e-6

node = NeuralODE(staticdudt,tspan,Tsit5(),save_everystep=false,save_start=false,sensealg=BacksolveAdjoint(autojacvec=ZygoteVJP()),p=p)
grads = Zygote.gradient(()->sum(node(x)),Flux.params(x,node))
@test !iszero(grads[x])
@test !iszero(grads[node.p])

@test_throws ErrorException grads = Zygote.gradient(()->sum(node(xs)),Flux.params(xs,node))

@info "Test some adjoints"

# Adjoint
@testset "adjoint mode" begin
    node = NeuralODE(dudt,tspan,Tsit5(),save_everystep=false,save_start=false)
    grads = Zygote.gradient(()->sum(node(x)),Flux.params(x,node))
    @test ! iszero(grads[x])
    @test ! iszero(grads[node.p])

    grads = Zygote.gradient(()->sum(node(xs)),Flux.params(xs,node))
    @test ! iszero(grads[xs])
    @test ! iszero(grads[node.p])

    node = NeuralODE(dudt,tspan,Tsit5(),saveat=0.0:0.1:1.0)
    grads = Zygote.gradient(()->sum(node(x)),Flux.params(x,node))
    @test ! iszero(grads[x])
    @test ! iszero(grads[node.p])

    grads = Zygote.gradient(()->sum(node(xs)),Flux.params(xs,node))
    @test ! iszero(grads[xs])
    @test ! iszero(grads[node.p])

    node = NeuralODE(dudt,tspan,Tsit5(),saveat=0.1)
    grads = Zygote.gradient(()->sum(node(x)),Flux.params(x,node))
    @test ! iszero(grads[x])
    @test ! iszero(grads[node.p])

    grads = Zygote.gradient(()->sum(node(xs)),Flux.params(xs,node))
    @test ! iszero(grads[xs])
    @test ! iszero(grads[node.p])

    node = NeuralODE(fastdudt,tspan,Tsit5(),save_everystep=false,save_start=false)
    grads = Zygote.gradient(()->sum(node(x)),Flux.params(x,node))
    @test ! iszero(grads[x])
    @test ! iszero(grads[node.p])

    @test_broken grads = Zygote.gradient(()->sum(node(xs)),Flux.params(xs,node)) isa Tuple
    @test_broken ! iszero(grads[xs])
    @test_broken ! iszero(grads[node.p])

    node = NeuralODE(fastcdudt,tspan,Tsit5(),save_everystep=false,save_start=false)
    grads = Zygote.gradient(()->sum(node(x)),Flux.params(x,node))
    @test ! iszero(grads[x])
    @test ! iszero(grads[node.p])

    @test_broken grads = Zygote.gradient(()->sum(node(xs)),Flux.params(xs,node)) isa Tuple
    @test_broken ! iszero(grads[xs])
    @test_broken ! iszero(grads[node.p])

    node = NeuralODE(fastdudt,tspan,Tsit5(),saveat=0.0:0.1:1.0)
    grads = Zygote.gradient(()->sum(node(x)),Flux.params(x,node))
    @test ! iszero(grads[x])
    @test ! iszero(grads[node.p])

    @test_broken grads = Zygote.gradient(()->sum(node(xs)),Flux.params(xs,node)) isa Tuple
    @test_broken ! iszero(grads[xs])
    @test_broken ! iszero(grads[node.p])

    node = NeuralODE(fastcdudt,tspan,Tsit5(),saveat=0.0:0.1:1.0)
    grads = Zygote.gradient(()->sum(node(x)),Flux.params(x,node))
    @test ! iszero(grads[x])
    @test ! iszero(grads[node.p])

    @test_broken grads = Zygote.gradient(()->sum(node(xs)),Flux.params(xs,node)) isa Tuple
    @test_broken ! iszero(grads[xs])
    @test_broken ! iszero(grads[node.p])

    node = NeuralODE(fastdudt,tspan,Tsit5(),saveat=0.1)
    grads = Zygote.gradient(()->sum(node(x)),Flux.params(x,node))
    @test ! iszero(grads[x])
    @test ! iszero(grads[node.p])

    @test_broken grads = Zygote.gradient(()->sum(node(xs)),Flux.params(xs,node)) isa Tuple
    @test_broken ! iszero(grads[xs])
    @test_broken ! iszero(grads[node.p])

    node = NeuralODE(fastcdudt,tspan,Tsit5(),saveat=0.1)
    grads = Zygote.gradient(()->sum(node(x)),Flux.params(x,node))
    @test ! iszero(grads[x])
    @test ! iszero(grads[node.p])

    @test_broken grads = Zygote.gradient(()->sum(node(xs)),Flux.params(xs,node)) isa Tuple
    @test_broken ! iszero(grads[xs])
    @test_broken ! iszero(grads[node.p])
end

@info "Test Tracker"

# RD
@testset "Tracker mode" begin
    node = NeuralODE(dudt,tspan,Tsit5(),save_everystep=false,save_start=false,sensealg=TrackerAdjoint())
    grads = Zygote.gradient(()->sum(node(x)),Flux.params(x,node))
    @test ! iszero(grads[x])
    @test ! iszero(grads[node.p])

    grads = Zygote.gradient(()->sum(node(xs)),Flux.params(xs,node))
    @test ! iszero(grads[xs])
    @test ! iszero(grads[node.p])

    node = NeuralODE(dudt,tspan,Tsit5(),saveat=0.0:0.1:1.0,sensealg=TrackerAdjoint())
    grads = Zygote.gradient(()->sum(node(x)),Flux.params(x,node))
    @test ! iszero(grads[x])
    @test ! iszero(grads[node.p])

    grads = Zygote.gradient(()->sum(node(xs)),Flux.params(xs,node))
    @test ! iszero(grads[xs])
    @test ! iszero(grads[node.p])

    node = NeuralODE(dudt,tspan,Tsit5(),saveat=0.1,sensealg=TrackerAdjoint())
    grads = Zygote.gradient(()->sum(node(x)),Flux.params(x,node))
    @test ! iszero(grads[x])
    @test ! iszero(grads[node.p])

    grads = Zygote.gradient(()->sum(node(xs)),Flux.params(xs,node))
    @test ! iszero(grads[xs])
    @test ! iszero(grads[node.p])

    node = NeuralODE(fastdudt,tspan,Tsit5(),save_everystep=false,save_start=false,sensealg=TrackerAdjoint())
    grads = Zygote.gradient(()->sum(node(x)),Flux.params(x,node))
    @test ! iszero(grads[x])
    @test ! iszero(grads[node.p])

    @test_broken grads = Zygote.gradient(()->sum(node(xs)),Flux.params(xs,node)) isa Tuple
    @test_broken ! iszero(grads[xs])
    @test_broken ! iszero(grads[node.p])

    node = NeuralODE(fastcdudt,tspan,Tsit5(),save_everystep=false,save_start=false,sensealg=TrackerAdjoint())
    grads = Zygote.gradient(()->sum(node(x)),Flux.params(x,node))
    @test ! iszero(grads[x])
    @test ! iszero(grads[node.p])

    @test_broken grads = Zygote.gradient(()->sum(node(xs)),Flux.params(xs,node)) isa Tuple
    @test_broken ! iszero(grads[xs])
    @test_broken ! iszero(grads[node.p])

    node = NeuralODE(fastdudt,tspan,Tsit5(),saveat=0.0:0.1:1.0,sensealg=TrackerAdjoint())
    grads = Zygote.gradient(()->sum(node(x)),Flux.params(x,node))
    @test ! iszero(grads[x])
    @test ! iszero(grads[node.p])

    @test_broken grads = Zygote.gradient(()->sum(node(xs)),Flux.params(xs,node)) isa Tuple
    @test_broken ! iszero(grads[xs])
    @test_broken ! iszero(grads[node.p])

    node = NeuralODE(fastcdudt,tspan,Tsit5(),saveat=0.0:0.1:1.0,sensealg=TrackerAdjoint())
    grads = Zygote.gradient(()->sum(node(x)),Flux.params(x,node))
    @test ! iszero(grads[x])
    @test ! iszero(grads[node.p])

    @test_broken grads = Zygote.gradient(()->sum(node(xs)),Flux.params(xs,node)) isa Tuple
    @test_broken ! iszero(grads[xs])
    @test_broken ! iszero(grads[node.p])

    node = NeuralODE(fastdudt,tspan,Tsit5(),saveat=0.1,sensealg=TrackerAdjoint())
    grads = Zygote.gradient(()->sum(node(x)),Flux.params(x,node))
    @test ! iszero(grads[x])
    @test ! iszero(grads[node.p])

    @test_broken grads = Zygote.gradient(()->sum(node(xs)),Flux.params(xs,node)) isa Tuple
    @test_broken ! iszero(grads[xs])
    @test_broken ! iszero(grads[node.p])

    node = NeuralODE(fastcdudt,tspan,Tsit5(),saveat=0.1,sensealg=TrackerAdjoint())
    grads = Zygote.gradient(()->sum(node(x)),Flux.params(x,node))
    @test ! iszero(grads[x])
    @test ! iszero(grads[node.p])

    @test_broken grads = Zygote.gradient(()->sum(node(xs)),Flux.params(xs,node)) isa Tuple
    @test_broken ! iszero(grads[xs])
    @test_broken ! iszero(grads[node.p])
end

@info "Test non-ODEs"

dudt2 = Flux.Chain(Flux.Dense(2,50,tanh),Flux.Dense(50,2))
fastdudt2 = FastChain(FastDense(2,50,tanh),FastDense(50,2))
fastcdudt2 = FastChain(FastDense(2,50,tanh,numcols=size(xs)[2],precache=true),FastDense(50,2,numcols=size(xs)[2],precache=true))
NeuralDSDE(dudt,dudt2,(0.0f0,.1f0),SOSRI(),saveat=0.1)(x)
sode = NeuralDSDE(dudt,dudt2,(0.0f0,.1f0),SOSRI(),saveat=0.0:0.01:0.1)

grads = Zygote.gradient(()->sum(sode(x)),Flux.params(x,sode))
@test ! iszero(grads[x])
@test ! iszero(grads[sode.p])
@test ! iszero(grads[sode.p][end])

grads = Zygote.gradient(()->sum(sode(xs)),Flux.params(xs,sode))
@test ! iszero(grads[xs])
@test ! iszero(grads[sode.p])
@test ! iszero(grads[sode.p][end])

sode = NeuralDSDE(fastdudt,fastdudt2,(0.0f0,.1f0),SOSRI(),saveat=0.0:0.01:0.1)
pd = Flux.params(sode)[1]
Random.seed!(1234)
gradsnc = Zygote.gradient(()->sum(sode(x)),Flux.params(x,sode)) #no cacahe
@test ! iszero(gradsnc[x])
@test ! iszero(gradsnc[sode.p])
@test ! iszero(gradsnc[sode.p][end])

gradsnc2 = Zygote.gradient(()->sum(sode(xs)),Flux.params(xs,sode))
@test_broken gradsnc2 isa Tuple
@test ! iszero(gradsnc2[xs])
@test ! iszero(gradsnc2[sode.p])
@test ! iszero(gradsnc2[sode.p][end])

sodec = NeuralDSDE(fastcdudt,fastcdudt2,(0.0f0,.1f0),SOSRI(),saveat=0.0:0.01:0.1,p=pd)
Random.seed!(1234)
gradsc = Zygote.gradient(()->sum(sodec(x)),Flux.params(x,sodec))
@test ! iszero(gradsc[x])
@test ! iszero(gradsc[sodec.p])
@test ! iszero(gradsc[sodec.p][end])
@test gradsnc[x] ≈ gradsc[x] rtol=1e-6
@test gradsnc[sode.p] ≈ gradsc[sodec.p] rtol=1e-6

gradsc2 = Zygote.gradient(()->sum(sodec(xs)),Flux.params(xs,sodec))
@test_broken gradsc2 isa Tuple
@test ! iszero(gradsc2[xs])
@test ! iszero(gradsc2[sodec.p])
@test ! iszero(gradsc2[sodec.p][end])

dudt22 = Flux.Chain(Flux.Dense(2,50,tanh),Flux.Dense(50,4),x->reshape(x,2,2))
fastdudt22 = FastChain(FastDense(2,50,tanh),FastDense(50,4),(x,p)->reshape(x,2,2))
fastcdudt22 = FastChain(FastDense(2,50,tanh,numcols=size(xs)[2],precache=true),FastDense(50,4,numcols=size(xs)[2],precache=true),(x,p)->reshape(x,2,2))
NeuralSDE(dudt,dudt22,(0.0f0,.1f0),2,LambaEM(),saveat=0.01)(x)

sode = NeuralSDE(dudt,dudt22,(0.0f0,0.1f0),2,LambaEM(),saveat=0.0:0.01:0.1)

grads = Zygote.gradient(()->sum(sode(x)),Flux.params(x,sode))
@test ! iszero(grads[x])
@test ! iszero(grads[sode.p])
@test ! iszero(grads[sode.p][end])

@test_broken grads = Zygote.gradient(()->sum(sode(xs)),Flux.params(xs,sode))
@test_broken ! iszero(grads[xs])
@test ! iszero(grads[sode.p])
@test ! iszero(grads[sode.p][end])

sode = NeuralSDE(fastdudt,fastdudt22,(0.0f0,0.1f0),2,LambaEM(),saveat=0.0:0.01:0.1)
pd = Flux.params(sode)[1]
Random.seed!(1234)
gradsnc = Zygote.gradient(()->sum(sode(x)),Flux.params(x,sode))
@test ! iszero(gradsnc[x])
@test ! iszero(gradsnc[sode.p])
@test ! iszero(gradsnc[sode.p][end])

@test_broken gradsnc = Zygote.gradient(()->sum(sode(xs)),Flux.params(xs,sode))
@test_broken ! iszero(gradsnc[xs])
@test ! iszero(gradsnc[sode.p])
@test ! iszero(gradsnc[sode.p][end])

sodec = NeuralSDE(fastcdudt,fastcdudt22,(0.0f0,0.1f0),2,LambaEM(),saveat=0.0:0.01:0.1,p=pd)
Random.seed!(1234)
gradsc = Zygote.gradient(()->sum(sodec(x)),Flux.params(x,sodec))
@test ! iszero(gradsc[x])
@test ! iszero(gradsc[sodec.p])
@test ! iszero(gradsc[sodec.p][end])
@test gradsnc[x] ≈ gradsc[x] rtol=1e-6
@test gradsnc[sode.p] ≈ gradsc[sodec.p] rtol=1e-6

@test_broken gradsc = Zygote.gradient(()->sum(sodec(xs)),Flux.params(xs,sodec))
@test_broken ! iszero(gradsc[xs])
@test ! iszero(gradsc[sodec.p])
@test ! iszero(gradsc[sodec.p][end])

ddudt = Flux.Chain(Flux.Dense(6,50,tanh),Flux.Dense(50,2))
NeuralCDDE(ddudt,(0.0f0,2.0f0),(p,t)->zero(x),(1f-1,2f-1),MethodOfSteps(Tsit5()),saveat=0.1)(x)
dode = NeuralCDDE(ddudt,(0.0f0,2.0f0),(p,t)->zero(x),(1f-1,2f-1),MethodOfSteps(Tsit5()),saveat=0.0:0.1:2.0)

grads = Zygote.gradient(()->sum(dode(x)),Flux.params(x,dode))
@test ! iszero(grads[x])
@test ! iszero(grads[dode.p])

@test_broken grads = Zygote.gradient(()->sum(dode(xs)),Flux.params(xs,dode)) isa Tuple
@test_broken ! iszero(grads[xs])
@test ! iszero(grads[dode.p])

fastddudt = FastChain(FastDense(6,50,tanh),FastDense(50,2))
NeuralCDDE(fastddudt,(0.0f0,2.0f0),(p,t)->zero(x),(1f-1,2f-1),MethodOfSteps(Tsit5()),saveat=0.1)(x)
dode = NeuralCDDE(fastddudt,(0.0f0,2.0f0),(p,t)->zero(x),(1f-1,2f-1),MethodOfSteps(Tsit5()),saveat=0.0:0.1:2.0)
pd = Flux.params(dode)[1]
gradsnc = Zygote.gradient(()->sum(dode(x)),Flux.params(x,dode))
@test ! iszero(gradsnc[x])
@test ! iszero(gradsnc[dode.p])

@test_broken gradsnc = Zygote.gradient(()->sum(dode(xs)),Flux.params(xs,dode)) isa Tuple
@test_broken ! iszero(gradsnc[xs])
@test ! iszero(gradsnc[dode.p])

fastcddudt = FastChain(FastDense(6,50,tanh,numcols=size(xs)[2],precache=true),FastDense(50,2,numcols=size(xs)[2],precache=true))
NeuralCDDE(fastcddudt,(0.0f0,2.0f0),(p,t)->zero(x),(1f-1,2f-1),MethodOfSteps(Tsit5()),saveat=0.1)(x)
dodec = NeuralCDDE(fastcddudt,(0.0f0,2.0f0),(p,t)->zero(x),(1f-1,2f-1),MethodOfSteps(Tsit5()),saveat=0.0:0.1:2.0,p=pd)

gradsc = Zygote.gradient(()->sum(dodec(x)),Flux.params(x,dodec))
@test ! iszero(gradsc[x])
@test ! iszero(gradsc[dodec.p])
@test gradsnc[x] ≈ gradsc[x] rtol=1e-6
@test gradsnc[dode.p] ≈ gradsc[dodec.p] rtol=1e-6

@test_broken gradsc = Zygote.gradient(()->sum(dodec(xs)),Flux.params(xs,dodec)) isa Tuple
@test_broken ! iszero(gradsc[xs])
@test ! iszero(gradsc[dodec.p])


text:
######################################################################
# Functions that are less imperative and more functional. These
# functions do a bit of smart magic in order to accomplish a given task.
# -----
# Licensed under MIT License

export negate
negate(x) = (args...; kwargs...) -> !x(args...; kwargs...)

export isathing
const isathing = negate(isnothing)

export truthy, falsy
truthy(::Nothing) = false
truthy(b::Bool)   = b
truthy(n::Number) = n != 0
truthy(_) = true
const falsy = negate(truthy)

function Base.split(cond, arr::AbstractArray)
    atrue  = []
    afalse = []
    for item ∈ arr
        if cond(item)
            push!(atrue, item)
        else
            push!(afalse, item)
        end
    end
    atrue, afalse
end
Base.split(cond, tpl::Tuple) = filter(item->cond(item), tpl), filter(item->!cond(item), tpl)

"""`insert!(arr::Vector{T}, elem::T; before = nothing, after = nothing, strict::Bool = false)`

 Insert `elem` into `arr` either immediately before `before` or immediately after `after`. Specifying both keyword
 arguments `before` and `after` as well as neither is illegal - exactly one must be specified. The specified relative
 element must exist in the array.

 If `strict` is false, simple equality (==) is used, otherwise strict equality (===)."""
function Base.insert!(arr::Vector{T}, elem::T; before = nothing, after = nothing, by = identity, strict::Bool = false) where T
    if before === nothing && after === nothing
        throw(ArgumentError("Neither `before` nor `after` specified"))
    elseif before !== nothing && after !== nothing
        throw(ArgumentError("Both `before` and `after` specified"))
    end

    idx = if before !== nothing
        indexof(arr, before; by=by, strict=strict)
    else
        idx = indexof(arr, after; by=by, strict=strict)
        idx === nothing ? nothing : idx+1
    end

    if idx === nothing
        ArgumentError("Relative element not found")
    end
    insert!(arr, idx, elem)
end

export indexed
abstract type Indexability end
struct Indexed <: Indexability end
struct Collectible <: Indexability end
struct Singular <: Indexability end
@generated function indexability(x)
    if hassignature(getindex, x, Integer) && hassignature(length, x) && hassignature(firstindex, x) && hassignature(lastindex, x)
        :(Indexed())
    elseif hassignature(iterate, x)
        :(Collectible())
    else
        :(Singular())
    end
end
indexed(v) = indexed(indexability(v), v)
indexed(::Indexed, v) = v
indexed(::Collectible, v) = collect(v)
indexed(::Singular, x) = (x,)

export iterable
abstract type Iterability end
struct Iterable <: Iterability end
struct NonIterable <: Iterability end
@generated function iterability(x)
    if hassignature(iterate, x)
        :(Iterable())
    else
        :(NonIterable())
    end
end
"""`iterable(x)`
If `x` is iterable (i.e. `iterate(x)` exists), return `x`. Otherwise, return an iterable container around `x`, e.g.
`tuple(x)`. The result of this function will always be iterable."""
iterable(x) = iterable(iterability(x), x)
iterable(::Iterable, x) = x
iterable(::NonIterable, x) = (x,)

export bytes
bytes(b::Union{UInt8, Int8}) = [convert(UInt8, b)]
bytes(t::Tuple) = bytes(collect(t))
function bytes(x)
    buf = IOBuffer()
    write(buf, x)
    take!(buf)
end


text:
"""
TransformedSet

Set obtained by some transformation from a feasible set (`generatingset`).
Should support `forward!` and `backward!` methods.
"""
abstract type TransformedSet <: FeasibleSet end

function generatingset(ts::TransformedSet)
error("Cannot find the generating set for $ts.")
end

"""
forward!(ts::TransformedSet)

In-place forward transform assosiated with the set `ts`: `x ∈ ts ⇔ f(x) ∈ ts.set`
Use `forward!(ts)(q,p)` to obtain `q = f(p)`.
"""
function forward!(ts::TransformedSet)
error("The forward transform is not defined for $ts.")
end

function backward!(ts::TransformedSet)
error("The backward transform is not defined for $ts.")
end

# it's easy to get an element of a transformed set
function getelement(ts::TransformedSet)
    el = getelement(generatingset(ts))
    backward!(ts)(getelement(generatingset(ts)))
end

"""
LinearTransformedSet

Subtype of TransformedSet where `forward` and `backward` transformations are given by multiplication
by forward and backward "plans" (i.e. -- precomputed matrices).
"""
abstract type LinearTransformedSet <: TransformedSet end

generatingset(s::LinearTransformedSet) = s.set
forward!(s::LinearTransformedSet) = ((xf, x) -> mul!(xf, s.fplan, x))
backward!(s::LinearTransformedSet) = (x -> mul!(x, s.bplan, x))

project!(xp, x, s::LinearTransformedSet) = (forward!(s)(xp, x); project!(xp, s.set); backward!(s)(xp))

struct FourierTransformedSet{TS,PF,PB} <: LinearTransformedSet where {TS <: FeasibleSet,PF <: AbstractFFTs.Plan,PB <: AbstractFFTs.Plan}
set::TS
fplan::PF
bplan::PB
end


text:
using VCFTools
using Test
using VariantCallFormat
using CodecZlib
using DelimitedFiles

# packages needed only for testing
using Random
using CSV
using DataFrames
using StatsBase
using LinearAlgebra

include("gtstats_test.jl")
include("conformgt_test.jl")
include("convert_test.jl")
include("filter_test.jl")
include("aim_test.jl")
include("grm_test.jl")


text:
import Term: make_logo

@testset "logo" begin
    @test_nothrow make_logo()

    logo = make_logo()
    @test logo.measure.w == 76
    @test logo.measure.h == 31

    testpanel(logo, 76, 31)
end


text:
@doc raw"""
    ScalarProduct()

The scalar product is an inner product of the form:

```math
f(\mathbf{x}, \mathbf{y}) = \mathbf{x}^{\intercal}\mathbf{y}
```
"""
struct ScalarProduct <: InnerProduct end
@inline base_aggregate(::ScalarProduct, s::T, x::T, y::T) where {T} = s + x*y


text:
using Reexport
using DataFrames
using CSV
using HTTP
using TOML
using Dates
using Statistics
using Distributions
using Optim
using Convex
using SCS
using JSON
using MathOptInterface
using StatsBase

# load my codes -
include("Types.jl")
include("Base.jl")
include("Network.jl")
include("Files.jl")
include("Compute.jl")
include("Transform.jl")
include("Factory.jl")
include("Simulations.jl")
include("Options.jl")

max_stars_repo_name: Pooksoft/Serenity

max_issues_repo_name: Alba-Intelligence/Serenity

max_forks_repo_name: Alba-Intelligence/Serenity