text:
# =========================================================================
#    CORE FUNCTIONS AND MECHANISMS
# =========================================================================

const max_i = Ref{Any}(nothing)
rmax() = max_i[] #buat manggil nilai max di i

function callMax!()
    val = col(
        JuMP.Containers.DenseAxisArray{Float64}(undef, keys(b().V)), #u
        JuMP.Containers.DenseAxisArray{Float64}(undef, keys(b().V)), #v
        JuMP.Containers.DenseAxisArray{Float64}(undef, keys(b().V), keys(b().V)), #empty l
        JuMP.Containers.DenseAxisArray{Float64}(undef, keys(b().V)), #y
        JuMP.Containers.DenseAxisArray{Float64}(undef, keys(b().V)), #z
        JuMP.Containers.DenseAxisArray{Float64}(undef, keys(b().V), keys(b().V)) #empty x
    )

    for i in keys(b().V)
        z = [b().K[k].BP[i] for k in keys(b().K)
            if i in b().K[k].cover
        ]
        val.z[i] = findmax(z)[1] #dapet nilai max z tiap titik

        y = [sum(b().K[k].BP[p] for p in b().K[k].cover) for k in keys(b().K)
            if i in b().K[k].cover
        ]
        val.y[i] = findmax(y)[1] #dapet nilai max z tiap titik

        v = [b().K[k].BP[i] * b().K[k].Q for k in keys(b().K)
            if i in b().K[k].cover
        ]
        val.v[i] = findmax(v)[1]

        u = [sum(b().K[k].BP[p] for p in b().K[k].cover) * b().K[k].Q
            for k in keys(b().K) if i in b().K[k].cover
        ]
        val.u[i] = findmax(u)[1]
    end

    return max_i[] = val
end

function Q(key,R::Dict)
    if isa(key,β)
        q = Vector{NamedTuple}()

        for r in keys(R), k in keys(b().K), t in b().T
            if getproperty(R[r],key.q)[key.i,k,t] >= key.v
                push!(q,(r=r,k=k,t=t))
            end
        end

        return q
    elseif isa(key,Vector{β})
        if !isempty(key)
            q = Vector{Vector{NamedTuple}}()

            for b in key
                push!(q,Q(b,R))
            end

            if !isempty(q)
                #ngembaliin semua rkt
                return reduce(intersect,q)
            else
                #ngembaliin vector kosong
                return q
            end
        else
            q = Vector{NamedTuple}()

            for r in keys(R), k in keys(b().K), t in b().T
                push!(q,(r=r,k=k,t=t))
            end

            return q
        end
    end
end

function f(B,R,θ)
    if !isempty(Q(B,R))
        return sum(θ[q.r,q.k,q.t] - floor(θ[q.r,q.k,q.t]) for q in Q(B,R))
    else
        return 0
    end
end

function tot(B,R,θ)
    if !isempty(Q(B,R))
        return sum(θ[q.r,q.k,q.t] for q in Q(B,R))
    else
        return 0
    end
end

function master(n::node)
    mp = buildMaster(n)
    optimize!(mp)

    return mp
end

function buildMaster(n::node)
    mp = Model(get_optimizer())
    if silent()
        set_silent(mp)
    end

    R = Dict(1:length(n.columns) .=> n.columns)
    B = Dict(1:length(n.bounds) .=> n.bounds)
    ≲ = filter(b -> last(b).type == "≲",B)
    ≳ = filter(b -> last(b).type == "≳",B)

    # ================================
    #    MODEL CONSTRUCTION
    # ================================
    iter_k = collect(keys(b().K))
    iter_i = collect(keys(b().V))

    @variable(mp, θ[keys(R), iter_k, b().T] >= 0)
    @variable(mp, I[iter_i, vcat(first(b().T) - 1, b().T)])
    @variable(mp, 0 <= slack[i = iter_i, t = b().T] <= n.stab.slLim[i,t])
    @variable(mp, 0 <= surp[i = iter_i, t = b().T] <= n.stab.suLim[i,t])

    @objective(mp, Min,
        sum(
            θ[r,k,t] * (
                sum(
                    b().dist[i,j] * (
                        b().K[k].vx * R[r].x[i,j,k,t] +
                        b().K[k].vl * R[r].l[i,j,k,t]
                    )
                    for i in b().K[k].cover, j in b().K[k].cover
                ) +
                sum(
                    b().K[k].fd * R[r].u[i,k,t]
                    for i in b().K[k].cover
                ) +
                sum(
                    b().K[k].fp * R[r].z[i,k,t]
                    for i in b().K[k].cover
                )
            )
            for r in keys(R), k in keys(b().K), t in b().T
        ) + #column costs
        sum(
            b().V[i].h * I[i,t]
            for i in iter_i, t in b().T
        ) + #inventory costs
        sum(
            n.stab.slCoeff * slack[i,t]
            for i in iter_i, t in b().T
        ) - #stabilizer
        sum(
            n.stab.suCoeff * surp[i,t]
            for i in iter_i, t in b().T
        ) #stabilizer
    )

    @constraint(mp, λ[i = iter_i, t = b().T],
        I[i,t - 1] +
        sum(R[r].u[i,k,t] * θ[r,k,t] for r in keys(R), k in iter_k) +
        slack[i,t] - surp[i,t] ==
        sum(R[r].v[i,k,t] * θ[r,k,t] for r in keys(R), k in iter_k) +
        b().d[i,t] + I[i,t]
    )

    @constraint(mp, δ[k = iter_k, i = b().K[k].cover, t = b().T],
        sum(R[r].z[i,k,t] * θ[r,k,t] for r in keys(R)) <= b().K[k].BP[i]
    )

    @constraint(mp, ϵ[k = iter_k, t = b().T],
        sum(θ[r,k,t] for r in keys(R)) <= sum(b().K[k].BP[i] for i in b().K[k].cover)
    )

    @constraint(mp, [i = keys(b().V), t = b().T],
        b().V[i].MIN <= I[i,t] <= b().V[i].MAX
    )

    @constraint(mp, [i = keys(b().V)],
        I[i,first(b().T)-1] == b().V[i].START
    )

    # ================================
    #    BOUND GENERATOR
    # ================================
    @constraint(mp, ρ[b = keys(≲)],
        sum(θ[q.r,q.k,q.t] for q in Q(B[b].B,R)) <= B[b].κ
    )

    @constraint(mp, σ[b = keys(≳)],
        sum(θ[q.r,q.k,q.t] for q in Q(B[b].B,R)) >= B[b].κ
    )

    return mp
end

function getDuals(mp::Model)
    λ = dual.(mp.obj_dict[:λ])
    δ = dual.(mp.obj_dict[:δ])
    ϵ = dual.(mp.obj_dict[:ϵ])
    ρ = dual.(mp.obj_dict[:ρ])
    σ = dual.(mp.obj_dict[:σ])

    return dval(λ,δ,ϵ,ρ,σ)
end

function sub(n::node,duals::dval)
    sp = buildSub(n,duals)
    optimize!(sp)

    return sp
end

const subproblem_model = Ref{Any}(nothing)
colStructure() = subproblem_model[] #pemanggil model subproblem

function callSub!()
    sp = Model(get_optimizer())
    if solver_name(sp) == "Gurobi"
        set_optimizer_attribute(sp,"Cuts",2)
        set_optimizer_attribute(sp,"MIPFocus",2)
        set_optimizer_attribute(sp,"Presolve",2)
        set_optimizer_attribute(sp,"Cutoff",0)
        set_optimizer_attribute(sp,"NodefileStart",0.5)
    end

    # ================================
    #    MODEL CONSTRUCTION
    # ================================
    iter_k = collect(keys(b().K))
    iter_i = collect(keys(b().V))

    q = col(
        @variable(sp, u[iter_i, iter_k, b().T] >= 0, Int),
        @variable(sp, v[iter_i, iter_k, b().T] >= 0, Int),
        @variable(sp, l[iter_i, iter_i, iter_k, b().T] >= 0, Int),
        @variable(sp, y[iter_i, iter_k, b().T], Bin),
        @variable(sp, z[iter_i, iter_k, b().T], Bin),
        @variable(sp, x[iter_i, iter_i, iter_k, b().T], Bin)
    )

    @constraint(sp, [k = iter_k, t = b().T],
        sum(u[i,k,t] for i in b().K[k].cover) ==
        sum(v[i,k,t] for i in b().K[k].cover) #all pickup delivered
    )

    @constraint(sp, [k = iter_k, i = b().K[k].cover, t = b().T],
        sum(x[j,i,k,t] for j in b().K[k].cover) == y[i,k,t] + z[i,k,t] #traverse in
    )

    @constraint(sp, [k = iter_k, i = b().K[k].cover, t = b().T],
        sum(x[i,j,k,t] for j in b().K[k].cover) == y[i,k,t] + z[i,k,t] #traverse out
    )

    @constraint(sp, [k = iter_k, i = b().K[k].cover, t = b().T],
        sum(l[j,i,k,t] for j in b().K[k].cover) -
        sum(l[i,j,k,t] for j in b().K[k].cover) == u[i,k,t] - v[i,k,t] #load balance
    )

    @constraint(sp, [i = iter_i, k = iter_k, t = b().T],
        u[i,k,t] <= b().K[k].Q * y[i,k,t] #u-y corr
    )

    @constraint(sp, [i = iter_i, k = iter_k, t = b().T],
        v[i,k,t] <= b().K[k].Q * z[i,k,t] #u-y corr
    )

    @constraint(sp, [i = iter_i, j = iter_i, k = iter_k, t = b().T],
        l[i,j,k,t] <= b().K[k].Q * x[i,j,k,t] #l-x corr
    )

    @constraint(sp, [k = iter_k, t = b().T],
        sum(z[i,k,t] for i in b().K[k].cover) <= 1 #only one starting point
    )

    for k in keys(b().K), t in b().T
        n = [i for i in iter_i if !(i in b().K[k].cover)] #sets not in cover

        if !isempty(n)
            @constraint(sp, [forbidden = n], z[forbidden,k,t] == 0)
            @constraint(sp, [forbidden = n], y[forbidden,k,t] == 0)
        end
    end

    optimize!(sp) #first call biar model kebuild

    return subproblem_model[] = sp
end

function buildSub(n::node,duals::dval)
    sp = colStructure()

    if silent()
        set_silent(sp)
    end

    #BRANCHING PREP
    B = Dict(1:length(n.bounds) .=> n.bounds)
    ≲ = filter(b -> last(b).type == "≲",B)
    ≳ = filter(b -> last(b).type == "≳",B)

    g = @variable(sp, [keys(≲), keys(b().K), b().T], Bin) #anonymous variable addition
    h = @variable(sp, [keys(≳), keys(b().K), b().T], Bin) #anonymous variable addition

    #ADD OBJECTIVE
    @objective(sp, Min,
        sum(
            sum(
                b().dist[i,j] * (
                    b().K[k].vx * sp.obj_dict[:x][i,j,k,t] +
                    b().K[k].vl * sp.obj_dict[:l][i,j,k,t]
                )
                for i in b().K[k].cover, j in b().K[k].cover
            ) +
            sum(
                b().K[k].fd * sp.obj_dict[:u][i,k,t]
                for i in b().K[k].cover
            ) +
            sum(
                b().K[k].fp * sp.obj_dict[:z][i,k,t]
                for i in b().K[k].cover
            )
            for k in keys(b().K), t in b().T
        ) -
        sum(
            sum(
                (sp.obj_dict[:u][i,k,t] - sp.obj_dict[:v][i,k,t]) * duals.λ[i,t]
                for i in b().K[k].cover
            )
            for k in keys(b().K), t in b().T
        ) -
        sum(
            sum(
                sp.obj_dict[:z][i,k,t] * duals.δ[k,i,t]
                for i in b().K[k].cover
            )
            for k in keys(b().K), t in b().T
        ) -
        sum(
            duals.ϵ[k,t]
            for k in keys(b().K), t in b().T
        ) -
        sum(
            sum(
                g[j,k,t] * duals.ρ[j]
                for j in keys(≲)
            )
            for k in keys(b().K), t in b().T
        ) -
        sum(
            sum(
                h[j,k,t] * duals.σ[j]
                for j in keys(≳)
            )
            for k in keys(b().K), t in b().T
        )
    )

    #BRANCHING CONSTRAINTS FOR SUBPROBLEM
    for j in keys(≲), k in keys(b().K), t in b().T
        η = @variable(sp, [B[j].B], Bin)

        @constraint(sp, -g[j,k,t] >= 1 - sum(1 - η[e] for e in B[j].B))
        for e in B[j].B
            @constraint(sp, (getproperty(rmax(),e.q)[e.i] - e.v + 1) * η[e] >=
                getproperty(q,e.q)[e.i,k,t] - e.v + 1
            )
        end
    end

    for j in keys(≳), k in keys(b().K), t in b().T
        η = @variable(sp, [B[j].B], Bin)

        @constraint(sp, [e = B[j].B], h[j,k,t] <= η[e])
        for e in B[j].B
            @constraint(sp, e.v * η[e] <= getproperty(q,e.q)[e.i,k,t])
        end
    end

    return sp
end

function getCols(sp::Model)
    u = value.(sp.obj_dict[:u])
    v = value.(sp.obj_dict[:v])
    l = value.(sp.obj_dict[:l])
    y = value.(sp.obj_dict[:y])
    z = value.(sp.obj_dict[:z])
    x = value.(sp.obj_dict[:x])

    return col(u,v,l,y,z,x)
end

function updateStab!(stab::stabilizer,param::Float64)
    for i in first(stab.slLim.axes),t in last(stab.slLim.axes)
        stab.slLim[i,t] = param * stab.slLim[i,t]
        if stab.slLim[i,t] < 1
            stab.slLim[i,t] = 0
        end
    end

    for i in first(stab.suLim.axes),t in last(stab.suLim.axes)
        stab.suLim[i,t] = param * stab.suLim[i,t]
        if stab.suLim[i,t] < 1
            stab.suLim[i,t] = 0
        end
    end

    return stab
end

function checkStab(mp::Model)
    s = (sum(value.(mp.obj_dict[:slack])) + sum(value.(mp.obj_dict[:surp])))

    return s
end

function colGen(n::node;maxCG::Float64,track::Bool)
    terminate = false
    iter = 0
    mem = 0

    while !terminate
        if iter < maxCG
            mp = master(n)

            if has_values(mp) && has_duals(mp)
                if track #print master problem obj
                    println("obj: $(objective_value(mp))")
                end

                duals = getDuals(mp)
                sp = sub(n,duals)

                if track #print subproblem price
                    println("price: $(objective_value(sp))")
                end

                if isapprox(objective_value(sp),0,atol = 1e-8) || objective_value(sp) > 0
                    if isapprox(checkStab(mp),0,atol = 1e-8)
                        terminate = true #action
                        push!(n.status,"EVALUATED") #report
                        if track
                            println("EVALUATED")
                        end
                    else
                        updateStab!(n.stab,0.5) #action
                        push!(n.status,"STABILIZED") #report
                        if track
                            println("STABILIZED")
                        end
                    end
                else
                    push!(n.columns,getCols(sp)) #action
                    push!(n.status,"ADD_COLUMN") #report
                    if track
                        println("ADD_COLUMN")
                    end
                end

                iter += 1 #iteration update
            else
                terminate = true #action
                push!(n.status,"NO_SOLUTION")
                if track
                    println("NO_SOLUTION")
                end
            end
        else
            terminate = true #action
            push!(n.status,"EVALUATED") #report
            if track
                println("EVALUATED")
            end
        end
    end

    if n.status[end] == "NO_SOLUTION"
        println("NODE $(n.self) FAILED.")
    else
        println("NODE $(n.self) FINISHED.")
    end

    return n
end

function origin(n::node)
    z = JuMP.Containers.DenseAxisArray{Float64}(undef,keys(b().V),keys(b().K),b().T)
    y = JuMP.Containers.DenseAxisArray{Float64}(undef,keys(b().V),keys(b().K),b().T)
    u = JuMP.Containers.DenseAxisArray{Float64}(undef,keys(b().V),keys(b().K),b().T)
    v = JuMP.Containers.DenseAxisArray{Float64}(undef,keys(b().V),keys(b().K),b().T)
    x = JuMP.Containers.DenseAxisArray{Float64}(
        undef,collect(keys(b().V)),collect(keys(b().V)),collect(keys(b().K)),b().T
    )
    l = JuMP.Containers.DenseAxisArray{Float64}(
        undef,collect(keys(b().V)),collect(keys(b().V)),collect(keys(b().K)),b().T
    )

    R = Dict(1:length(n.columns) .=> n.columns)
    mp = master(n)
    θ = mp.obj_dict[:θ]

    for i in keys(b().V), k in keys(b().K), t in b().T
        z[i,k,t] = value(sum(R[r].z[i,k,t] * θ[r,k,t] for r in keys(R)))
        y[i,k,t] = value(sum(R[r].y[i,k,t] * θ[r,k,t] for r in keys(R)))
        u[i,k,t] = value(sum(R[r].u[i,k,t] * θ[r,k,t] for r in keys(R)))
        v[i,k,t] = value(sum(R[r].v[i,k,t] * θ[r,k,t] for r in keys(R)))

        for j in keys(b().V)
            x[i,j,k,t] = value(sum(R[r].x[i,j,k,t] * θ[r,k,t] for r in keys(R)))
            l[i,j,k,t] = value(sum(R[r].l[i,j,k,t] * θ[r,k,t] for r in keys(R)))
        end
    end

    return col(
        u,v,l,
        y,z,x
    )
end


text:
@testset "Variables" begin

    m = Model()

    @variable(m, x ≥ 0, u"m/s")
    @test x == UnitJuMP.UnitVariableRef(x.vref, u"m/s")
    @test unit(x) == u"m/s"
    @test owner_model(x) === m

    @variable(m, y[1:4], u"km/hr")
    @test y[1] == UnitJuMP.UnitVariableRef(y[1].vref, u"km/hr")

end

@testset "Constraints" begin

    m = Model()

    @variable(m, u[1:2] ≥ 0, u"m")
    @variable(m, v[1:4], u"km/hr")
    @variable(m, w, Bin)
    @variable(m, y)

    maxspeed = 4u"ft/s"
    speed = 2.3u"m/s"

    # Various combination of coefficient and variables with and without units
    @constraint(m, 2v[1] ≤ maxspeed)
    @test num_constraints(m, AffExpr, MOI.LessThan{Float64}) == 1

    @constraint(m, 2.3v[1] ≤ maxspeed)
    @test num_constraints(m, AffExpr, MOI.LessThan{Float64}) == 2

    @constraint(m, c1a, speed * w ≤ maxspeed)
    @test num_constraints(m, AffExpr, MOI.LessThan{Float64}) == 3
    @test unit(c1a) == u"m/s"

    @constraint(m, c1b, 40u"km/hr" * w ≤ 15u"m/s")
    @test num_constraints(m, AffExpr, MOI.LessThan{Float64}) == 4
    @test unit(c1b) == u"km/hr"

    @constraint(m, c1c, w * speed  ≤ maxspeed)
    @test num_constraints(m, AffExpr, MOI.LessThan{Float64}) == 5
    @test unit(c1c) == u"m/s"

    @constraint(m, c1d, sum(v[i] for i in 1:4)   ≤ maxspeed)
    @test num_constraints(m, AffExpr, MOI.LessThan{Float64}) == 6
    @test unit(c1d) == u"km/hr"

    @constraint(m, c1, 2v[1] + 4v[2] ≤ maxspeed)
    @test typeof(c1) <: UnitJuMP.UnitConstraintRef
    @test unit(c1) == u"km/hr"

    @constraint(m, c2, 2v[1] + 4v[2] ≤ maxspeed, u"m/s")
    @test unit(c2) == u"m/s"
    @test normalized_rhs(c2.cref) == convert(Float64, uconvert(u"m/s", maxspeed).val)

    @variable(m, z, Bin)

    maxlength = 1000u"yd"
    period = 1.5u"hr"
    @constraint(m, c3, u[2] + period * v[2] ≤ maxlength * z, u"cm")
    @test unit(c3) == u"cm"

    @constraint(m, c3b, u[2] + 1.5u"hr" * v[2] ≤ 1000u"yd" * z, u"cm")
    @test unit(c3b) == u"cm"

end

ext: jl


text:
using TikzPictures
using Test

# Pre-test cleanup (for repeated tests)
for file in ["testPic.pdf", "testPic.svg", "testDoc.pdf", "testDoc.tex"]
	if isfile(file)
		rm(file)
	end
end

# Run tests
tp = TikzPicture("\\draw (0,0) -- (10,10);\n\\draw (10,0) -- (0,10);\n\\node at (5,5) {tikz \$\\sqrt{\\pi}\$};", options="scale=0.25", preamble="")
td = TikzDocument()
push!(td, tp, caption="hello")

save(TEX("testPic"), tp)
@test isfile("testPic.tex")

if success(`lualatex -v`)
	save(PDF("testPic"), tp)
	@test isfile("testPic.pdf")

    save(SVG("testPic"), tp)
    @test isfile("testPic.svg")

    save(PDF("testDoc"), td)
    @test isfile("testDoc.pdf")
end


text:
using Ekztazy
using Distributed

client = Client()

ENV["JULIA_DEBUG"] = Ekztazy

TESTGUILD = ENV["TESTGUILD"]

on_message_create!(client) do (ctx)
    (!isme(client, ctx)) && reply(client, ctx, content="$(mention(ctx.message.author)), $(ctx.message.content) TEST")
end

command!(client, TESTGUILD, "boom", "Go boom!") do (ctx)
    reply(client, ctx, content="$(mention(ctx)) blew up!")
end

command!(client, TESTGUILD, "bam", "Go bam!") do (ctx)
    reply(client, ctx, content="$(mention(ctx)) slapped themselves!")
end

command!(client, TESTGUILD, "double", "Doubles a number!", legacy=false, options=[
    Option(Int, name="number", description="The number to double!")
]) do ctx, number
    reply(client, ctx, content="$(number*2)")
end

command!(client, TESTGUILD, "multiply", "Multiplies numbers!", legacy=false, options=Options(
    [Int, "a", "the first number"],
    [Int, "b", "the second number"]
)) do ctx, a::Int, b::Int
    reply(client, ctx, content="$(a*b)")
end

command!(client, TESTGUILD, "greet", "Greets a user", legacy=false, options=Options(
    [User, "u", "The user to greet"]
)) do ctx, u::Member
    reply(client, ctx, content="Hello, $(u)!")
end

command!(client, TESTGUILD, "water", "Water a plant", legacy=false, options=[
    Option(Int, "howmuch", "How long do you want to water the plant?")
]) do ctx, howmuch
    cm = component!(client, "magic"; auto_ack=false, type=2, style=1, label="Wow, a Button!?") do context
        update(client, context, content="You pressed the button!")
    end
    reply(client, ctx, components=[cm], content="$(mention(ctx)) watered their plant for $(howmuch) hours. So much that the plant grew taller than them!")
end

command!(client, TESTGUILD, "test", "Test something", legacy=false) do ctx
    cm = select!(
        client,
        "class",
        ("Rogue", "rogue"),
        ("Mage", "mage");
        max_values=1
    ) do context, choices
        reply(client, context, content="You chose $(choices[1])")
    end
    reply(client, ctx, components=[cm], content="What class you want noob?")
end

command!(client, TESTGUILD, "quit", "Ends the bot process!") do (ctx)
    reply(client, ctx, content="Shutting down the bot")
    close(client)
end

on_ready!(client) do (ctx)
    @info "Successfully logged in as $(ctx.user.username)" # obtain(client, User).username
end

start(client)


text:
function go_around(r, side)
    counter_horisontal = 0
    counter_vertical = 0
    while(isborder(r, side))
        move!(r, turn_side(side)) # turn_side(side)
        counter_vertical = counter_vertical + 1
    end
    move!(r, side)
    while(isborder(r, invert_side(turn_side(side)))) # invert_side(turn_side(side))
        move!(r, side)
        counter_horisontal = counter_horisontal + 1
    end
    move!(r, invert_side(turn_side(side)))
    c = counter_vertical
    while(c > 1)
        move!(r, invert_side(turn_side(side))) # invert_side(turn_side(side))
        c = c - 1
    end
    # counter_horisontal - в побочную сторону
    # counter_vertical - в главную сторону
    return (counter_horisontal, counter_vertical)
end
function turn_side(side)
    if side == Nord
        return West
    end
    if side == Sud
        return West
    end
    if side == Ost
        return Sud
    end
    if side == West
        return Sud
    end
end
function invert_side(side)
    if side == Nord
        return Sud
    end
    if side == Sud
        return Nord
    end
    if side == Ost
        return West
    end
    if side == West
        return Ost
    end
end
function gotobegin(r)
    a = 0
    b = 0
    while(!((isborder(r, West)) & (isborder(r, Sud))))
        a = a + gotoborder(r, West)
        b = b + gotoborder(r, Sud)
    end
    return(a, b)
end
function gotoborder(r, side)
    count = 0
    while(!isborder(r, side))
        move!(r, side)
        count = count + 1
    end
    return count
end
function measure_field(r)
    y = gotoborder(r, Nord)
    gotoborder(r, Sud)
    x = gotoborder(r, Ost)
    gotoborder(r, West)
    return (x,y)
end
function push_and_mark(r, side, counter)
    b = 0
    while(counter>0)
        putmarker!(r)
        if(isborder(r, side))
            b,count = go_around(r, side)
            counter = counter - count - 1
        else
            move!(r, side)
            counter = counter - 1
        end
    end
end
function push_hor(r, side, counter)
    b = 0
    while(counter>0)
        if(isborder(r, side))
            count,b = go_around(r, side)
            counter = counter - count - 1
        else
            move!(r, side)
            counter = counter - 1
        end
    end
end
function goto_coord(r, side, coord)
    while(coord > 0)
        coord = coord - 1
        move!(r, side)
    end
end
function main(r)
    c_hor,c_vert = gotobegin(r)
    rt_x,rt_y=c_hor,c_vert
    x_dim,y_dim = measure_field(r)
    side = Nord
    move!(r, Ost)
    d = 1
    obstacle_count = 0
    while(x_dim>0)
        d = d + 1
        x_dim = x_dim - 1
        c = y_dim
        while(c>0)
            if(isborder(r, side))
                count,b = go_around(r, side)
                c = c - count - 1
                if((b == 1) & (count*b != 0))
                    obstacle_count = obstacle_count + 1
                end
                # print(count, " ", b, "\n")
            else
                move!(r, side)
                c = c - 1
            end
        end
        gotobegin(r)
        if(x_dim != 0)
            goto_coord(r, Ost, d)
        end
    end
    obstacle_count = obstacle_count
    # print(obstacle_count)
    return obstacle_count
end


text:
# functions to help with creating submeshes

"""
  Copy the balues from the second array to the first, using the mask to map
  the last dimension.  For example:

  A[:, :, idx] = B[:, :, list[i]] if i != 0

  This function always copies the elements of the arrays. and preserves the
  order of `list`

  **Inputs**

   * src: source array
   * list: list of elements to copy (used for last dimension of array)

  **Inputs/Outputs**

   * dest: destination array
"""
function copy_masked(dest::AbstractVector, src::AbstractVector,
                     list::AbstractVector)

  @assert length(list) == size(dest, 1)

  pos = 1
  for i in list
    dest[pos] = copy(src[i])
    pos += 1
  end
end

function copy_masked(dest::AbstractMatrix, src::AbstractMatrix,
                     list::AbstractVector)

  @assert length(list) == size(dest, 2)
  @assert size(dest, 1) == size(src, 1)

  pos = 1
  for i in list
    for j=1:size(dest, 1)
      dest[j, pos] = copy(src[j, i])
    end
    pos += 1
  end
end

function copy_masked(dest::Abstract3DArray, src::Abstract3DArray,
                     list::AbstractVector)

  @assert length(list) == size(dest, 3)
  @assert size(dest, 1) == size(src, 1)
  @assert size(dest, 2) == size(src, 2)

  pos = 1
  for i in list
    for k=1:size(dest, 2)
      for j=1:size(dest, 1)
        dest[j, k, pos] = copy(src[j, k, i])
      end
    end
    pos += 1
  end  # end if
end   # end function

function copy_masked(dest::Abstract4DArray, src::Abstract4DArray,
                     list::AbstractVector)

  @assert length(list) == size(dest, 4)
  @assert size(dest, 1) == size(src, 1)
  @assert size(dest, 2) == size(src, 2)
  @assert size(dest, 3) == size(src, 3)

  pos = 1
  for i in list
    for p=1:size(dest, 3)
      for k=1:size(dest, 2)
        for j=1:size(dest, 1)
          dest[j, k, p, pos] = copy(src[j, k, p, i])
        end
      end
    end
    pos += 1
  end  # end if
end   # end function

"""
  Copies the main data arrays from one mesh to a submesh

  TODO: list all arrays

  **Inputs**

   * oldmesh: the original mesh
   * el_list: list of elements to copy (order is preserved)

  **Inputs/Outputs**

   * mesh: the submesh
"""
function copy_data_arrays(mesh::PumiMesh, oldmesh::PumiMesh, el_list::AbstractVector)

  mesh.vert_coords = zeros(T1, mesh.dim, mesh.numVertPerElement, mesh.numEl)
  copy_masked(mesh.vert_coords, oldmesh.vert_coords, el_list)

  mesh.coords = zeros(T1, mesh.dim, mesh.numNodesPerElement, mesh.numEl)
  copy_masked(mesh.coords, oldmesh.coords, el_list)

  return nothing
end

"""
  original mesh -> submesh  for vectors (indexed using mesh.dofs)

  **Inputs**

   * oldmesh: the original mesh
   * oldq: the AbstractVector of values to be copied to the new mesh (for
           the elements that exist on the new mesh), length oldmesh.numdof
   * newmesh: the submesh

  **Inputs/Outputs**

   * newq: vector of length newmesh.numdof to put values into
"""
function injectionOperator(oldmesh::PumiMesh, oldq::AbstractVector,
                           newmesh::PumiMesh, newq::AbstractVector)

  @assert oldmesh.m_ptr == apf.getOldMesh(newmesh.subdata)

  parent_N = apf.getParentNumbering(newmesh.subdata)

  for i=1:newmesh.numEl
    el_i = newmesh.elements[i]
    oldel = apf.getNumberJ(parent_N, el_i, 0, 0) + 1

    for j=1:newmesh.numNodesPerElement
      for k=1:newmesh.numDofPerNode
        olddof = oldmesh.dofs[k, j, oldel]
        newdof = newmesh.dofs[k, j, i]
        newq[newdof] = oldq[olddof]
      end
    end
  end

  return nothing
end

"""
  Injection operator for 3D arrays, see other method for details.
"""
function injectionOperator(oldmesh::PumiMesh, oldq::Abstract3DArray,
                           newmesh::PumiMesh, newq::Abstract3DArray)

  @assert oldmesh.m_ptr == apf.getOldMesh(newmesh.subdata)

  parent_N = apf.getParentNumbering(newmesh.subdata)

  for i=1:newmesh.numEl
    el_i = newmesh.elements[i]
    oldel = apf.getNumberJ(parent_N, el_i, 0, 0) + 1

    for k=1:newmesh.numNodesPerElement
      for j=1:newmesh.numDofPerNode
        newq[j, k, i] = oldq[j, k, oldel]
      end
    end
  end

  return nothing
end

"""
  submesh -> original mesh for vectors, inverse function of
  [`injectionOperator`](@ref)

  **Inputs**

   * newmesh: the submesh
   * newq: vector of length newmesh.numdof
   * oldmesh: the original mesh

  **Inputs/Outputs**

   * oldq: vector of length oldmesh.numdof.  Entries that have corresponding
           dofs on the submesh are overwritten, other entries are unchanged.

"""
function rejectionOperator(newmesh::PumiMesh, newq::AbstractVector,
                           oldmesh::PumiMesh, oldq::AbstractVector)

  @assert oldmesh.m_ptr == apf.getOldMesh(newmesh.subdata)

  parent_N = apf.getParentNumbering(newmesh.subdata)

  for i=1:newmesh.numEl
    el_i = newmesh.elements[i]
    oldel = apf.getNumberJ(parent_N, el_i, 0, 0) + 1

    for j=1:newmesh.numNodesPerElement
      for k=1:newmesh.numDofPerNode
        olddof = oldmesh.dofs[k, j, oldel]
        newdof = newmesh.dofs[k, j, i]
        oldq[olddof] = newq[newdof]
      end
    end
  end

  return nothing
end

"""
  rejectionOperator() method for 3D arrays, see other method for details
"""
function rejectionOperator(newmesh::PumiMesh, newq::Abstract3DArray,
                           oldmesh::PumiMesh, oldq::Abstract3DArray)

  @assert oldmesh.m_ptr == apf.getOldMesh(newmesh.subdata)

  parent_N = apf.getParentNumbering(newmesh.subdata)

  for i=1:newmesh.numEl
    el_i = newmesh.elements[i]
    oldel = apf.getNumberJ(parent_N, el_i, 0, 0) + 1

    for k=1:newmesh.numNodesPerElement
      for j=1:newmesh.numDofPerNode
        oldq[j, k, oldel] = newq[j, k, i]
      end
    end
  end

  return nothing
end

"""
  This function returns an Array of Interface objects in the same order as
  newmesh.bndryfaces where the element that also exists on the submesh is
  elementL.  This can be useful for interpolating the solution for the new
  boundary condition needed for the submesh.
  This boundary condion must always be the last one.

  Note: 2D only (for now)

  **Inputs**

   * oldmesh: the original mesh
   * newmesh: the submesh

  **Outputs**

   * interface_arr: array of Interface objects
"""
function getBoundaryInterpArray(oldmesh::PumiMesh2D, newmesh::PumiMesh2D)

  rng = newmesh.bndry_offsets[end-1]:(newmesh.bndry_offsets[end]-1)
  interface_arr = Array{Interface}(length(rng))

  parent_N = apf.getParentNumbering(newmesh.subdata)
  edges = Array{Ptr{Void}}(12)  # apf::Downward
  el_arr = Array{Ptr{Void}}(2)  # number of elements that share an edge
  pos = 1  # position in interface_arr

  for i in rng
    el_i = newmesh.bndryfaces[i].element
    face_i = newmesh.bndryfaces[i].face
    el_ptr = newmesh.elements[el_i]

    parent_el = apf.getNumberJ(parent_N, el_ptr, 0, 0) + 1
    parent_elptr = oldmesh.elements[parent_el]

    apf.getDownward(oldmesh.m_ptr, parent_elptr, oldmesh.dim-1, edges)
    edge_ptr = edges[face_i]

    nup = apf.countAdjacent(oldmesh.m_ptr, edge_ptr, oldmesh.dim)
    @assert nup == 2
    apf.getAdjacent(el_arr)

    # figure out which is left and which is right
    elL_ptr = parent_elptr
    if el_arr[1] == elL_ptr
      elR_ptr = el_arr[2]
    else
      elR_ptr = el_arr[1]
    end

   # get local face numbers
   elnumL = apf.getNumberJ(oldmesh.el_Nptr, elL_ptr, 0, 0) + 1
   elnumR = apf.getNumberJ(oldmesh.el_Nptr, elR_ptr, 0, 0) + 1
   edgenum = apf.getNumberJ(oldmesh.face_Nptr, edge_ptr, 0, 0) + 1

   facelocalnumL = getFaceLocalNum(oldmesh, edgenum, elnumL)
   facelocalnumR = getFaceLocalNum(oldmesh, edgenum, elnumR)

   interface_arr[pos] = Interface(elnumL, elnumR, facelocalnumL, facelocalnumR, UInt8(1))  # orientation is always 1 for 2d
   pos += 1
 end

  return interface_arr
end


text:
using AbstractTrees
using AbstractTrees: ImplicitIndexStack, ImplicitNodeStack, stepstate, getnode

import AbstractTrees: children
import Base: &, <<, |, >>, ==, ~

struct ChildMask
    mask::BitVector
end
Base.copy(x::ChildMask) = ChildMask(copy(x.mask))
Base.ones(::Type{ChildMask}) = ChildMask(fill!(BitVector(256), true))
Base.broadcast(op::Union{typeof(|),typeof(&)},a::ChildMask, b::ChildMask) = ChildMask(broadcast(op,a.mask,b.mask))
Base.broadcast!(op::Union{typeof(|),typeof(&)},xs::ChildMask...) = ChildMask(broadcast!(op,map(x->x.mask, xs)...))
Base.any(m::ChildMask) = any(m.mask)
Base.broadcast(op::typeof(~),a::ChildMask) = ChildMask(broadcast(op,a.mask))
Base.:(==)(a::ChildMask, b::ChildMask) = a.mask == b.mask
Base.hash(c::ChildMask, u::UInt64) = Base.hash(c.mask, u)
Base.in(b::UInt8, c::ChildMask) = c.mask[b+1]

# Shoe the child mask using braille. One braille character per byte
function Base.show(io::IO, mask::ChildMask)
    for byte in reinterpret(UInt8, mask.mask.chunks)
        print(io, Char(0x2800+byte))
    end
end

struct IOTreeLeaf
    mask::ChildMask
    value::Any
end
function Base.showcompact(io::IO, node::IOTreeLeaf)
    Base.show(io, node.mask); print(io, ": ")
    Base.show(io, node.value)
end

struct IOTreeNode
    mask::ChildMask
    children::Vector{Union{IOTreeNode, IOTreeLeaf}}
end
Base.showcompact(io::IO, node::IOTreeNode) = Base.show(io, node.mask)
IOTreeNode(mask::ChildMask) = IOTreeNode(mask, Union{IOTreeNode, IOTreeLeaf}[])
IOTreeNode(mask::ChildMask, children::IOTreeNode) = IOTreeNode(mask, children.children)

struct IOTreeRoot
    children::Vector{Union{IOTreeNode, IOTreeLeaf}}
end
const IOTreeObject = Union{IOTreeNode, IOTreeLeaf, IOTreeRoot}
IOTreeRoot() = IOTreeRoot(IOTreeObject[])
children(tree::Union{IOTreeRoot,IOTreeNode}) = tree.children

function split!(tree, left_mask::ChildMask, right_mask::ChildMask)
    @assert !any(left_mask .& right_mask)
    total_mask = (left_mask .| right_mask)
    for (i,c) in enumerate(children(tree))
        isa(c, IOTreeLeaf) && continue
        # If the total mask is disjoint, we don't care
        any(total_mask .& c.mask) || continue
        # If the child is entirely contained within one of the masks,
        # we don't need to do anything.
        (!any(left_mask .& c.mask) ||
         !any(right_mask .& c.mask)) && continue
        # Ok, here we need to split. First, duplicate the children
        new_children = deepcopy(c.children)
        tree.children[i] = IOTreeNode(c.mask .& left_mask, c.children)
        push!(tree.children, IOTreeNode(c.mask .& right_mask, new_children))
    end
    tree
end
split!(tree, l, r) = split!(tree, ChildMask(l), ChildMask(r))

struct AbstractImmediate
  size::Int
end
Base.rem(x::AbstractImmediate, ::Type{T}) where {T} = T(0)

mutable struct GenerativeIOTree <: IO
    tree::IOTreeRoot
    state
    idx::Int
    hadImmediate::Bool
    nbytes_read::Int
end

struct AbstractBit
  # (-1%UInt8) indicates concrete
  byte_pos::UInt8
  bit_pos::UInt8
  # For concrete bits, is this 0 or 1
  set::Bool
end

struct AbstractBitBundle{N}
  tree::Union{GenerativeIOTree, Void}
  # 1 = low bit
  abstract_bits::NTuple{N, AbstractBit}
end
const AbstractByte = AbstractBitBundle{8}

function Base.:&(x::AbstractByte, mask::UInt8)
  AbstractByte(x.tree, tuple((((mask >> (i-1)) & 1) == 1 ? x.abstract_bits[i] : concrete_bit(0) for i = 1:8)...))
end

concrete_bit(x::Integer) = (@assert x == 0 || x == 1; AbstractBit(-1%UInt8,-1%UInt8,x!=0))
isconcrete(x::AbstractBit) = x.byte_pos == -1%UInt8 && x.bit_pos == -1%UInt8
function Base.convert(::Type{AbstractByte}, x::UInt8)
  AbstractByte(nothing, tuple((concrete_bit((x >> i) & 1) for i = 0:7)...))
end
Base.convert(::Type{AbstractByte}, x::Integer) = AbstractByte(convert(UInt8, x))
Base.showcompact(io::IO, tree::IOTreeRoot) = print(io, "<Root>")

(::Type{UInt16})(x::AbstractBitBundle{N}) where {N} =
  AbstractBitBundle{16}(x.tree, tuple((i > N ? concrete_bit(0) : x.abstract_bits[i] for i = 1:16)...))

function select_bit(x, y, i)
    i > length(x.abstract_bits) && return y.abstract_bits[i]
    i > length(y.abstract_bits) && return x.abstract_bits[i]
    xbit = x.abstract_bits[i]; ybit = y.abstract_bits[i]
    (isconcrete(xbit) && xbit.set) && return xbit
    (isconcrete(ybit) && ybit.set) && return ybit
    @assert !(!isconcrete(xbit) && !isconcrete(ybit))
    isconcrete(xbit) ? ybit : xbit
end
function (|)(x::AbstractBitBundle{N}, y::AbstractBitBundle{M}) where {N,M}
  xtree = x.tree
  ytree = y.tree
  (xtree === nothing) && (xtree = ytree)
  (ytree === nothing) && (ytree = xtree)
  @assert xtree === ytree
  AbstractBitBundle{max(N,M)}(xtree, tuple((select_bit(x, y, i) for i = 1:max(N,M))...))
end
function (|)(x::AbstractBitBundle{N}, mask::Integer) where N
  AbstractBitBundle{N}(x.tree, tuple((((mask >> (i-1)) & 1) == 1 ? concrete_bit(1) : x.abstract_bits[i] for i = 1:N)...))
end
function (<<)(x::AbstractBitBundle{N}, shift::Integer) where N
  AbstractBitBundle{N}(x.tree, tuple(((i - shift >= 1) ? x.abstract_bits[i - shift] : concrete_bit(0) for i = 1:N)...))
end
function (>>)(x::AbstractBitBundle{N}, shift::Integer) where N
  AbstractBitBundle{N}(x.tree, tuple(((i + shift <= N) ? x.abstract_bits[i + shift] : concrete_bit(0) for i = 1:N)...))
end
~(x::AbstractBit) = AbstractBit(x.byte_pos, x.bit_pos, !x.set)
(~)(x::AbstractBitBundle{N}) where {N} = AbstractByte(x.tree, tuple((~x.abstract_bits[i] for i = 1:N)...))
function (&)(x::AbstractBitBundle{N}, mask::Integer) where N
  AbstractBitBundle{N}(x.tree, tuple((((mask >> (i-1)) & 1) == 1 ? x.abstract_bits[i] : concrete_bit(0) for i = 1:N)...))
end
(&)(mask::UInt8, x::AbstractByte) = x & mask

Base.rem(x::AbstractBitBundle, ::Type{UInt8}) =
  AbstractByte(x.tree, tuple((x.abstract_bits[i] for i = 1:8)...))

function descend(tree::IOTreeRoot, idxs)
    node = tree
    for idx in idxs
        node = node.children[idx]
    end
    node
end

function Base.read(tree::GenerativeIOTree, ::Type{UInt8})
    tree.idx += 1
    if tree.idx - 1 > length(tree.state.idx_stack)
        node = getnode(tree.tree, tree.state)
        isa(node, IOTreeLeaf) && @show node
        if isempty(node.children)
            new_node = IOTreeNode(ones(ChildMask), IOTreeObject[])
            push!(node.children, new_node)
        end
        # Make sure to maintain state invariants
        if !isempty(tree.state.idx_stack)
            push!(tree.state.node_stack, node)
        end
        push!(tree.state.idx_stack, 1)
    end
    tree.nbytes_read += 1
    bits = tuple((AbstractBit(tree.idx - 1, i, false) for i = 1:8)...)
    return AbstractByte(tree, bits)
end

function bitwise_mask(i)
    idxs = filter(x->(x & (UInt8(1)<<i)) != 0, 0:255)
    mask = fill!(BitVector(256), false)
    mask[idxs+1] = true
    ChildMask(mask)
end
const bitwise_masks = [bitwise_mask(i) for i = 0:7]

const ones_mask = ChildMask(fill!(BitVector(256), true))
const zeros_mask = ChildMask(fill!(BitVector(256), false))
function (==)(x::AbstractBitBundle, y::Integer)
    bits = sort(map(x->((i,b)=x; (b.byte_pos, b.bit_pos, (y & (1<<(i-1))) != 0)), enumerate(x.abstract_bits)))
    # Mask of values for which the condition is true
    mask = copy(ones_mask)
    bpos = bits[1][1]
    determined = true
    equal = true
    while true
        if isempty(bits) || (bit = shift!(bits); bit[1] != bpos)
            parent = bpos == 1 ? x.tree.tree : x.tree.state.node_stack[bpos-1]
            node = children(parent)[x.tree.state.idx_stack[bpos]]
            comb = node.mask .& mask
            if comb == node.mask
                # nothing to do
            elseif comb == zeros_mask
                equal = false
            else
                # We split this in such a way that the answer is always
                # false if we get here
                split!(parent, .~(mask), mask)
                determined = false
            end
            isempty(bits) && break
            mask = copy(ones_mask)
            bpos = bit[1]
        end
        (bpos == 0xff) && break
        mask .&= bit[3] ? bitwise_masks[bit[2]] : .~(bitwise_masks[bit[2]])
    end
    determined ? equal : false
end
(==)(x::Integer, y::AbstractBitBundle) = y == x

function Base.isless(a::AbstractBitBundle, b::Integer)
    error()
end
function Base.isless(a::Integer, b::AbstractBitBundle)
    error()
end
# TODO: Do this more efficienty
function Base.getindex(a::Dict, b::AbstractBitBundle)
    for (key, val) in a
        (key == b) && return val
    end
    error()
end
function Base.haskey(a::Dict, b::AbstractBitBundle)
    for key in keys(a)
        (key == b) && return true
    end
    return false
end
function Base.in(a::AbstractBitBundle, b::Vector)
    for x in b
        (a == x) && return true
    end
    return false
end

function next!(tree::GenerativeIOTree, v)
    parent = isempty(tree.state.node_stack) ? tree.tree : tree.state.node_stack[end]
    children(parent)[tree.state.idx_stack[end]] =
        IOTreeLeaf(children(parent)[tree.state.idx_stack[end]].mask, v)
    state = tree.state
    while true
        state = stepstate(Leaves(tree.tree), state)
        isnull(state) && return false
        state = get(state)
        !isa(getnode(tree.tree, state), IOTreeLeaf) && break
    end
    tree.state = state
    tree.idx = 1
    tree.hadImmediate = false
    tree.nbytes_read = 0
    return true
end


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule Wine_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("Wine")
JLLWrappers.@generate_main_file("Wine", UUID("9fae3aff-8997-5dd1-9b84-5d0cc5e0bffa"))
end  # module Wine_jll


text:
function _precompile_()
    ccall(:jl_generating_output, Cint, ()) == 1 || return nothing
    @assert precompile(Tuple{typeof(detect_bedgraph),IOStream})
    @assert precompile(Tuple{typeof(detect_noometiff),IOStream})
    @assert precompile(Tuple{typeof(detect_rdata),IOStream})
    @assert precompile(Tuple{typeof(detect_rdata_single),IOStream})
    @assert precompile(Tuple{typeof(detectwav),IOStream})

    @assert precompile(Tuple{typeof(load),File})
    @assert precompile(Tuple{typeof(load),Formatted})
    @assert precompile(Tuple{typeof(load),String})

    @assert precompile(Tuple{typeof(query),String})
    @assert precompile(Tuple{typeof(query),IOStream})
    @assert precompile(Tuple{typeof(query),IOStream,String})
    @assert precompile(Tuple{typeof(query),IOStream,Nothing})

    @assert precompile(Tuple{typeof(hasfunction),Function})
    @assert precompile(Tuple{typeof(hasmagic),Function})

    @assert precompile(Tuple{typeof(applicable_loaders),Type{<:DataFormat}})
    @assert precompile(Tuple{typeof(applicable_loaders),Formatted})
    @assert precompile(Tuple{typeof(applicable_savers),Type{<:DataFormat}})
    @assert precompile(Tuple{typeof(applicable_savers),Formatted})
    @assert precompile(Tuple{typeof(add_loader),Type{<:DataFormat},Symbol})
    @assert precompile(Tuple{typeof(add_saver),Type{<:DataFormat},Symbol})

    @assert precompile(Tuple{typeof(iter_eq),Array{UInt8,1},NTuple{10,UInt8}})
    @assert precompile(Tuple{typeof(iter_eq),Array{UInt8,1},NTuple{20,UInt8}})
    @assert precompile(Tuple{typeof(iter_eq),Array{UInt8,1},NTuple{30,UInt8}})
    @assert precompile(Tuple{typeof(iter_eq),Array{UInt8,1},NTuple{32,UInt8}})
    @assert precompile(Tuple{typeof(iter_eq),Array{UInt8,1},NTuple{35,UInt8}})
    @assert precompile(Tuple{typeof(iter_eq),Array{UInt8,1},NTuple{4,UInt8}})
    @assert precompile(Tuple{typeof(iter_eq),Array{UInt8,1},NTuple{6,UInt8}})
    @assert precompile(Tuple{typeof(iter_eq),Array{UInt8,1},NTuple{7,UInt8}})
    @assert precompile(Tuple{typeof(iter_eq),Array{UInt8,1},NTuple{8,UInt8}})
    @assert precompile(Tuple{typeof(iter_eq),Array{UInt8,1},Tuple{UInt8,UInt8,UInt8}})
    @assert precompile(Tuple{typeof(iter_eq),Array{UInt8,1},Tuple{UInt8,UInt8}})

    if isdefined(Base, :bodyfunction)
        m = which(query, (String,))
        f = Base.bodyfunction(m)
        @assert precompile(f, (Bool, typeof(query), String))
        m = which(load, (String,))
        f = Base.bodyfunction(m)
        @assert precompile(f, (Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}, typeof(load), String))
        m = which(load, (Formatted,))
        f = Base.bodyfunction(m)
        @assert precompile(f, (Any, typeof(load), Formatted))
        @assert precompile(f, (Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}, typeof(load), File))
    end

end


text:
struct LODF{Ax, L <: NTuple{2, Dict}} <: PowerNetworkMatrix{Float64}
    data::Array{Float64, 2}
    axes::Ax
    lookup::L
end

function _buildlodf(branches, nodes, dist_slack::Array{Float64} = [0.1])
    linecount = length(branches)
    ptdf, a = PowerSystems._buildptdf(branches, nodes, dist_slack)
    H = gemm('N', 'N', ptdf, a)
    ptdf_denominator = H
    for iline in 1:linecount
        if (1.0 - ptdf_denominator[iline, iline]) < 1.0E-06
            ptdf_denominator[iline, iline] = 0.0
        end
    end
    (Dem, dipiv, dinfo) = getrf!(
        Matrix{Float64}(LinearAlgebra.I, linecount, linecount) -
            Array(LinearAlgebra.Diagonal(ptdf_denominator)),
    )
    lodf = gemm('N', 'N', H, getri!(Dem, dipiv))
    lodf =
        lodf - Array(LinearAlgebra.Diagonal(lodf)) -
        Matrix{Float64}(LinearAlgebra.I, linecount, linecount)
    return lodf
end

function LODF(branches, nodes, dist_slack::Array{Float64} = [0.1])

    #Get axis names
    line_ax = [branch.name for branch in branches]
    lodf = _buildlodf(branches, nodes, dist_slack)

    axes = (line_ax, line_ax)
    look_up = (_make_ax_ref(line_ax), _make_ax_ref(line_ax))

    return LODF(lodf, axes, look_up)

end

function LODF(sys::System, dist_slack::Array{Float64} = [0.1])
    branches = get_components(ACBranch, sys)
    nodes = get_components(Bus, sys)

    return LODF(branches, nodes, dist_slack)

end


text:
@testset "Lattices" begin
  # Smoke test for genus symbol
  Qx, x = PolynomialRing(FlintQQ, "x")
  K, a = NumberField(x^2 - 2, "a")
  L = quadratic_lattice(K, identity_matrix(K, 10), gram_ambient_space = 35*identity_matrix(K, 10))
  P = prime_decomposition(maximal_order(K), 5)[1][1]
  @test Hecke._genus_symbol_kirschmer(L, P).data == [(10, 1, 1)]

  P = prime_decomposition(maximal_order(K), 2)[1][1]
  GS = Hecke._genus_symbol_kirschmer(L, P)
  @test GS.data[1] == [10]
  @test GS.data[2] == [0]
  @test GS.data[3] == [2]
  @test GS.data[4] == [35]
  @test GS.data[5] == []

  Qx, x = FlintQQ["x"];
  K, a = number_field(x^2 - 2, cached = true);
  Kt, t = K["t"];
  L, b= number_field(t^2 + 11, "b", check = true)
  p = prime_decomposition(maximal_order(K), 2)[1][1]
  P = prime_decomposition(maximal_order(L), p)[1][1]
  H = Hecke.lattice(hermitian_space(K, 2 * identity_matrix(K, 3)), pseudo_matrix(identity_matrix(K, 3), [p, p, p]))
  @test Hecke._genus_symbol_kirschmer(H, fmpz(2)) == Any[(3, 4, true, 4, -64)]
  @test islocally_isometric(H, H, fmpz(2))

  H = Hecke.lattice(hermitian_space(L, L(elem_in_nf(uniformizer(p))) * identity_matrix(L, 3)), pseudo_matrix(identity_matrix(L, 3), [P, P, P]))
  @test Hecke._genus_symbol_kirschmer(H, p) == Any[(3, 3, false)]
  @test islocally_isometric(H, H, p)
  @test Hecke.genus(H, p) == genus(HermLat, L, p, [(3, 3, -1)])

  V = hermitian_space(L, L(a) * identity_matrix(L, 3))
  M = Hecke.maximal_integral_lattice(V)
  @test Hecke.genus(M, p) == genus(HermLat, L, p, [(0, 2, 1), (1, 1, -1)])
  q = prime_decomposition(maximal_order(K), 11)[1][1]
  @test Hecke.genus(M, q) == genus(HermLat, L, q, [(-1, 2, 1), (0, 1, 1)])

  L, _ = number_field(t^2 - gen(K) * t + 1)
  V = hermitian_space(L, L(a) * identity_matrix(L, 3))
  M = Hecke.maximal_integral_lattice(V)
  @test Hecke.genus(M, p) == genus(HermLat, L, p, [(-2, 2, 1, 0), (0, 1, -1, 0)])

  V = hermitian_space(L, L(10) * identity_matrix(L, 3))
  M = Hecke.maximal_integral_lattice(V)
  @test Hecke.genus(M, p) == genus(HermLat, L, p, [(-2, 2, 1, 0), (0, 1, 1, 0)])

end


text:
module TestMoments
    using Distributions
    using StatsBase
    using StreamStats
    using Base.Test

    # Moments of uniform draws
    for n in rand(1:1_000_000, 100)
        xs = rand(n)
        stat = StreamStats.Moments()
        for x in xs
            update!(stat, x)
        end
        online_ms, online_vs, online_ss, online_ks = state(stat)
        online_m = mean(xs)
        online_v = var(xs)
        online_s = skewness(xs)
        online_k = kurtosis(xs)
        online_n = nobs(stat)
        batch_m = mean(xs)
        batch_v = var(xs)
        batch_s = skewness(xs)
        batch_k = kurtosis(xs)

        @test_approx_eq(online_m, batch_m)
        @test_approx_eq(online_ms, batch_m)
        @test_approx_eq(online_v, batch_v)
        @test_approx_eq(online_vs, batch_v)
        @test abs(online_s - batch_s) < 1e-8 * abs(batch_s)
        @test abs(online_ss - batch_s) < 1e-8 * abs(batch_s)
        @test abs(online_k - batch_k) < 1e-8 * abs(batch_k)
        @test abs(online_ks - batch_k) < 1e-8 * abs(batch_k)
        @test online_n == n
    end
end


text:
@testset "216.combination-sum-iii.jl" begin
    @test combination_sum3(3, 7) == [[1, 2, 4]]
    @test combination_sum3(3, 9) == [[1,2,6], [1,3,5], [2,3,4]]
    @test combination_sum3(5, 20) == [[1,2,3,5,9],[1,2,3,6,8],[1,2,4,5,8],[1,2,4,6,7],[1,3,4,5,7],[2,3,4,5,6]]
    @test combination_sum3(5, 23) == [[1,2,3,8,9],[1,2,4,7,9],[1,2,5,6,9],[1,2,5,7,8],[1,3,4,6,9],[1,3,4,7,8],[1,3,5,6,8],[1,4,5,6,7],[2,3,4,5,9],[2,3,4,6,8],[2,3,5,6,7]]
end


text:
using SeisIO, SeisNoise, Dates, Glob, Base.Threads
export CC

function get_alltime_path_cc(FFTDIR::String, select_time::Bool, TIMETXT::String)
    if select_time == true
        alltime_path = readlines(TIMETXT)
    else
        alltime_path = sort(glob("*/*/*", FFTDIR))
    end
    return alltime_path
end

function get_SR_pairs_cc(select_sta::String, STATXT::String, SRTXT::String)
    All_NET_STA = Vector{Vector{String}}(undef, 0)
    SR_pairs = Vector{Vector{String}}(undef, 0)
    if select_sta == "no"
        nothing
    elseif select_sta == "sta"
        # a.read STATXT --> station_txt
        station_txt = readlines(STATXT)  # read sta.txt
        # b.make all stations --> rAll_NET_STA (for reading data[sac/mseed] from database and auto-corr)
        N_sta = size(station_txt, 1)
        rAll_NET_STA = Vector{Vector{String}}(undef, N_sta)
        for ll = 1:N_sta
            D_net, D_sta = split(station_txt[ll])
            rAll_NET_STA[ll] = [string(D_net), string(D_sta)]
        end
        # c.make s-r pairs --> rSR_pairs (for cross-corr)
        N_SR = Int(N_sta * (N_sta - 1) / 2)
        rSR_pairs = Vector{Vector{String}}(undef, N_SR)
        num = 0
        for ss = 1:N_sta-1
            for rr = ss+1:N_sta
                num += 1
                S_net, S_sta = split(station_txt[ss])
                R_net, R_sta = split(station_txt[rr])
                rSR_pairs[num] =
                    [string(S_net), string(S_sta), string(R_net), string(R_sta)]
            end
        end
        # d.return
        All_NET_STA = rAll_NET_STA
        SR_pairs = rSR_pairs
    elseif select_sta == "sr"
        # a.read SRTXT --> station_txt
        sr_txt = readlines(SRTXT)
        # b.generate pSR_pairs with repeated sta-pairs --> pSR_pairs
        source_mark = 1   # "1" mark the source of one read-line
        S_net = ""
        S_sta = ""
        pSR_pairs = Vector{Vector{String}}(undef, 0)
        for i = 1:size(sr_txt, 1)
            if sr_txt[i] != ""
                if source_mark == 1
                    s_net, s_sta = split(sr_txt[i])
                    S_net = s_net
                    S_sta = s_sta
                    source_mark = -1   # "-1" mark next line is not the source
                else
                    R_net, R_sta = split(sr_txt[i])
                    if !((S_net == R_net) && (S_sta == R_sta))
                        push!(
                            pSR_pairs,
                            [string(S_net), string(S_sta), string(R_net), string(R_sta)],
                        )
                    end
                end
            else
                source_mark = 1   # mark next line is the source
            end
        end
        # c.delete the repeated pSR_pairs --> SR_pairs
        rSR_pairs = Vector{Vector{String}}(undef, 0)
        repeat_mark = -1   # "-1" means not repeat
        for i = 1:size(pSR_pairs, 1)-1
            S_net_1, S_sta_1, R_net_1, R_sta_1 =
                pSR_pairs[i][1], pSR_pairs[i][2], pSR_pairs[i][3], pSR_pairs[i][4]
            for j = i+1:size(pSR_pairs, 1)
                S_net_2, S_sta_2, R_net_2, R_sta_2 =
                    pSR_pairs[j][1], pSR_pairs[j][2], pSR_pairs[j][3], pSR_pairs[j][4]
                condition_1 =
                    (S_net_1 == S_net_2) &&
                    (S_sta_1 == S_sta_2) &&
                    (R_net_1 == R_net_2) &&
                    (R_sta_1 == R_sta_2)   # Completely corresponding
                condition_2 =
                    (S_net_1 == R_net_2) &&
                    (S_sta_1 == R_sta_2) &&
                    (R_net_1 == S_net_2) &&
                    (R_sta_1 == S_sta_2)   # Reverse S and R, then corresponde
                if condition_1 || condition_2
                    repeat_mark = 1
                    break
                else
                    repeat_mark = -1
                end
            end
            if repeat_mark != 1
                push!(
                    rSR_pairs,
                    [string(S_net_1), string(S_sta_1), string(R_net_1), string(R_sta_1)],
                )
            end
        end
        # The last group will never repeat
        S_net_1, S_sta_1, R_net_1, R_sta_1 =
            pSR_pairs[end][1], pSR_pairs[end][2], pSR_pairs[end][3], pSR_pairs[end][4]
        push!(
            rSR_pairs,
            [string(S_net_1), string(S_sta_1), string(R_net_1), string(R_sta_1)],
        )

        # d.make all stations --> All_NET_STA (for read data[sac/mseed] from FFT-database and auto-corr)
        double_SR_pairs = Vector{Vector{String}}(undef, 0)
        for i = 1:size(rSR_pairs, 1)
            S_net, S_sta, R_net, R_sta =
                rSR_pairs[i][1], rSR_pairs[i][2], rSR_pairs[i][3], rSR_pairs[i][4]
            push!(double_SR_pairs, [S_net, S_sta])
            push!(double_SR_pairs, [R_net, R_sta])
        end
        rAll_NET_STA = Vector{Vector{String}}(undef, 0)
        NET_1 = ""
        STA_1 = ""
        repeat_mark = -1
        for i = 1:size(double_SR_pairs, 1)-1
            net_1, sta_1 = double_SR_pairs[i][1], double_SR_pairs[i][2]
            NET_1, STA_1 = net_1, sta_1
            for j = i+1:size(double_SR_pairs, 1)
                net_2, sta_2 = double_SR_pairs[j][1], double_SR_pairs[j][2]
                if NET_1 == net_2 && STA_1 == sta_2
                    repeat_mark = 1
                    break
                else
                    repeat_mark = -1
                end
            end
            if repeat_mark != 1
                push!(rAll_NET_STA, [string(NET_1), string(STA_1)])
            end
        end
        NET_1, STA_1 = double_SR_pairs[end][1], double_SR_pairs[end][2]
        push!(rAll_NET_STA, [string(NET_1), string(STA_1)])
        # e.conclusion --> All_NET_STA,SR_pairs
        All_NET_STA = rAll_NET_STA
        SR_pairs = rSR_pairs
    end

    ### modify sequence of SR_pairs
    net_1,sta_1,net_2,sta_2 = "","","",""
    for ii = 1:size(SR_pairs,1)
        net_1,sta_1,net_2,sta_2 = SR_pairs[ii][1],SR_pairs[ii][2],SR_pairs[ii][3],SR_pairs[ii][4]
        if net_1*sta_1 > net_2*sta_2
            SR_pairs[ii][1] = net_2
            SR_pairs[ii][2] = sta_2
            SR_pairs[ii][3] = net_1
            SR_pairs[ii][4] = sta_1
        end
    end
    return sort(All_NET_STA), sort(SR_pairs)
end

function get_allfile_list_cc(Onetime_path::String, channel_regular::String)
    files = Vector{String}(undef, 0)
    files_cha = glob(channel_regular, Onetime_path)
    append!(files, files_cha)
    sort!(files)
    return files
end

function get_FFT_dict_cc(Fs::Vector{FFTData})
    Fs_dict = Dict{String,FFTData}()
    for i = 1:size(Fs, 1)
        net, sta, kel, cha = split(Fs[i].name, '.')
        key = string(net, '.', sta)
        setindex!(Fs_dict, Fs[i], key)
    end
    return Fs_dict
end

function get_FFTData_cc(
    Onetime_path::String,
    channel_regular::String,
    comp::String,
    select_sta::String,
    STATXT::String,
    SRTXT::String,
)
    # 1.read one time-chunck paths
    files = get_allfile_list_cc(Onetime_path, channel_regular)
    # 2.read data
    N_files = size(files, 1)
    All_NET_STA = Vector{Vector{String}}(undef, 0)
    SR_pairs = Vector{Vector{String}}(undef, 0)
    Fs = Vector{FFTData}(undef, 0)
    if select_sta == "no"
        # a.read all data --> Fs
        for j = 1:N_files
            file = files[j]
            Fs_one = ""
            try
                fs_one = load_fft(file, comp)
                Fs_one = deepcopy(fs_one)
            catch e
                println("Errorpass in function 'get_FFTData':", e)
                continue
            end
            push!(Fs, Fs_one)
        end
        if size(Fs, 1) == 0   # one time-chunck include 0 files
            return Fs, 0, 0
        end
        # b.make all stations --> All_NET_STA (for read data[sac/mseed] from database and auto-corr)
        N_sta = size(Fs, 1)
        rAll_NET_STA = Vector{Vector{String}}(undef, 0)
        for j = 1:N_sta
            net, sta, kno, cha = split(Fs[j].name, '.')
            push!(rAll_NET_STA, [string(net), string(sta)])
        end
        All_NET_STA = deepcopy(rAll_NET_STA)
        # c.make s-r pairs --> SR_pairs (for cross-corr)
        N_SR = Int(N_sta * (N_sta - 1) / 2)
        rSR_pairs = Vector{Vector{String}}(undef, N_SR)
        num = 0
        for ss = 1:N_sta-1
            for rr = ss+1:N_sta
                num += 1
                S_net, S_sta = rAll_NET_STA[ss]
                R_net, R_sta = rAll_NET_STA[rr]
                rSR_pairs[num] = [S_net, S_sta, R_net, R_sta]
            end
        end
        SR_pairs = deepcopy(rSR_pairs)
    else
        # a.read sta-pairs --> rAll_NET_STA, rSR_pairs
        rAll_NET_STA, rSR_pairs = get_SR_pairs_cc(select_sta, STATXT, SRTXT)
        # b.read data from path with selected 'All_NET_STA' --> error_net_sta, Fs
        error_net_sta = Vector{Vector{String}}(undef, 0)
        read_mark = 1
        for ii = 1:size(rAll_NET_STA, 1)
            for j = 1:N_files
                file = files[j]
                if occursin(rAll_NET_STA[ii][1], basename(file)) &&
                   occursin(rAll_NET_STA[ii][2], basename(file))
                    read_mark = 1
                    try
                        Fs_one = load_fft(file, comp)
                        push!(Fs, Fs_one)
                        break
                    catch e
                        push!(error_net_sta, [rAll_NET_STA[ii][1], rAll_NET_STA[ii][2]])  # read error
                        println("Errorpass in function 'get_FFTData':", e)
                        continue
                    end
                else
                    read_mark = -1
                end
            end
            if read_mark == -1  # could not find file with net-sta provided by sta.txt
                push!(error_net_sta, [rAll_NET_STA[ii][1], rAll_NET_STA[ii][2]])
            end
        end
        if size(Fs, 1) == 0   # one time-chunck include 0 files
            return Fs, 0, 0
        end
        # c.modify rAll_NET_STA, rSR_pairs --> All_NET_STA, SR_pairs
        error_num = length(error_net_sta)
        if error_num == 0
            All_NET_STA = rAll_NET_STA
            SR_pairs = rSR_pairs
        else
            # modify All_NET_STA
            for i = 1:size(rAll_NET_STA, 1)
                net_1, sta_1 = rAll_NET_STA[i][1], rAll_NET_STA[i][2]
                mark_sta = -1
                for j = 1:error_num
                    net_2, sta_2 = error_net_sta[j][1], error_net_sta[j][2]
                    if net_1 == net_2 && sta_1 == sta_2
                        mark_sta = 1
                        break
                    end
                end
                if mark_sta == -1
                    push!(All_NET_STA, [net_1, sta_1])
                end
            end
            # modify SR_pairs
            for i = 1:size(rSR_pairs, 1)
                s_net_1, s_sta_1, r_net_1, r_sta_1 =
                    rSR_pairs[i][1], rSR_pairs[i][2], rSR_pairs[i][3], rSR_pairs[i][4]
                mark_sta = -1
                for j = 1:error_num
                    net_2, sta_2 = error_net_sta[j][1], error_net_sta[j][2]
                    if (s_net_1 == net_2 && s_sta_1 == sta_2) ||
                        (r_net_1 == net_2 && r_sta_1 == sta_2)
                        mark_sta = 1
                        break
                    end
                end
                if mark_sta == -1
                    push!(SR_pairs, [s_net_1, s_sta_1, r_net_1, r_sta_1])
                end
            end
        end
    end
    # 3.convert Fs --> Fs_dict
    Fs_dict = get_FFT_dict_cc(Fs)

    return Fs_dict, All_NET_STA, SR_pairs
end

function get_cc_substack_cc(
    Onetime_path::String,
    Fs_dict::Dict{String,FFTData},
    All_NET_STA::Vector{Vector{String}},
    SR_pairs::Vector{Vector{String}},
    auto_corr::Bool,
    cross_corr::Bool,
    cc_method::String,
    CCFDIR::String,
    maxlag::Float64,
    substack::Bool,
    substack_method::Vector{String},
)
    year = splitpath(Onetime_path)[end-2]
    yearday = splitpath(Onetime_path)[end-1]
    daytime = splitpath(Onetime_path)[end]

    println(All_NET_STA)
    println(SR_pairs)
    if auto_corr == true
        Threads.@threads for i = 1:size(All_NET_STA, 1)
            # debug:
            # println("thread num = ", Threads.threadid())
            #
            net, sta = All_NET_STA[i][1], All_NET_STA[i][2]
            key = string(net, '.', sta)
            C = correlate(Fs_dict[key], Fs_dict[key], maxlag, corr_type = cc_method)
            if substack == true
                for jj in substack_method
                    if jj == "mean"
                        C_mean = stack(C, allstack = true)
                        SPLITOUT_mean =
                            joinpath(CCFDIR, "substack_mean", year, yearday, daytime, key)
                        if !isdir(SPLITOUT_mean)
                            mkpath(SPLITOUT_mean)
                        end
                        save_corr(C_mean, SPLITOUT_mean)
                    end
                    if jj == "pws"
                        C_pws = pws(C)
                        SPLITOUT_pws =
                            joinpath(CCFDIR, "substack_pws", year, yearday, daytime, key)
                        if !isdir(SPLITOUT_pws)
                            mkpath(SPLITOUT_pws)
                        end
                        save_corr(C_pws, SPLITOUT_pws)
                    end
                    if jj == "robust"
                        C_robust = deepcopy(C)
                        C_robust.corr = robustpws(C.corr)
                        SPLITOUT_robust =
                            joinpath(CCFDIR, "substack_robust", year, yearday, daytime, key)
                        if !isdir(SPLITOUT_robust)
                            mkpath(SPLITOUT_robust)
                        end
                        save_corr(C_robust, SPLITOUT_robust)
                    end
                end
            else
                SPLITOUT_no = joinpath(CCFDIR, "substack_no", year, yearday, daytime, key)
                if !isdir(SPLITOUT_no)
                    mkpath(SPLITOUT_no)
                end
                save_corr(C, SPLITOUT_no)
            end
        end
    end

    if cross_corr == true
        Threads.@threads for i = 1:size(SR_pairs, 1)
            s_net, s_sta, r_net, r_sta =
                SR_pairs[i][1], SR_pairs[i][2], SR_pairs[i][3], SR_pairs[i][4]
            key_1 = string(s_net, '.', s_sta)
            key_2 = string(r_net, '.', r_sta)
            C = correlate(Fs_dict[key_1], Fs_dict[key_2], maxlag, corr_type = cc_method)
            if substack == true
                for jj in substack_method
                    if jj == "mean"
                        C_mean = stack(C, allstack = true)
                        SPLITOUT_mean =
                            joinpath(CCFDIR, "substack_mean", year, yearday, daytime, key_1)
                        if !isdir(SPLITOUT_mean)
                            mkpath(SPLITOUT_mean)
                        end
                        save_corr(C_mean, SPLITOUT_mean)
                    end
                    if jj == "pws"
                        C_pws = pws(C)
                        SPLITOUT_pws =
                            joinpath(CCFDIR, "substack_pws", year, yearday, daytime, key_1)
                        if !isdir(SPLITOUT_pws)
                            mkpath(SPLITOUT_pws)
                        end
                        save_corr(C_pws, SPLITOUT_pws)
                    end
                    if jj == "robust"
                        C_robust = deepcopy(C)
                        C_robust.corr = robustpws(C.corr)
                        SPLITOUT_robust = joinpath(
                            CCFDIR,
                            "substack_robust",
                            year,
                            yearday,
                            daytime,
                            key_1,
                        )
                        if !isdir(SPLITOUT_robust)
                            mkpath(SPLITOUT_robust)
                        end
                        save_corr(C_robust, SPLITOUT_robust)
                    end
                end
            else
                SPLITOUT_no = joinpath(CCFDIR, "substack_no", year, yearday, daytime, key_1)
                if !isdir(SPLITOUT_no)
                    mkpath(SPLITOUT_no)
                end
                save_corr(C, SPLITOUT_no)
            end
        end
    end

    return nothing
end

function main_cc(Onetime_path::String, Parameters::Dict{String,Any})
    flag = Parameters["flag"]
    channel_regular = Parameters["channel_regular"]
    comp = Parameters["comp"]
    select_sta = Parameters["select_sta"]
    STATXT = Parameters["STATXT"]
    SRTXT = Parameters["SRTXT"]
    auto_corr = Parameters["auto_corr"]
    cross_corr = Parameters["cross_corr"]
    cc_method = Parameters["cc_method"]
    CCFDIR = Parameters["CCFDIR"]
    maxlag = Parameters["maxlag"]
    substack = Parameters["substack"]
    substack_method = Parameters["substack_method"]
    println("Cross-correlating Onetime_day: $(Onetime_path) || now:", now())

    if flag == true
        ###### 1.read FFTData ######
        tALL = @elapsed begin
            Fs_dict, All_NET_STA, SR_pairs = get_FFTData_cc(
                Onetime_path,
                channel_regular,
                comp,
                select_sta,
                STATXT,
                SRTXT,
            )
            if All_NET_STA == 0 && SR_pairs == 0
                return nothing
            end
        end
        println("'Func:get_FFTData' took $tALL seconds")
        ###### 2.cc_substack ######
        tALL = @elapsed begin
            get_cc_substack_cc(
                Onetime_path,
                Fs_dict,
                All_NET_STA,
                SR_pairs,
                auto_corr,
                cross_corr,
                cc_method,
                CCFDIR,
                maxlag,
                substack,
                substack_method,
            )
        end
        println("'Func:get_cc_substack' took $tALL seconds")
    else
        ###### 1.read FFTData ######
        Fs_dict, All_NET_STA, SR_pairs =
            get_FFTData_cc(Onetime_path, channel_regular, comp, select_sta, STATXT, SRTXT)
        if All_NET_STA == 0 && SR_pairs == 0
            return nothing
        end
        ###### 2.cc_substack ######
        get_cc_substack_cc(
            Onetime_path,
            Fs_dict,
            All_NET_STA,
            SR_pairs,
            auto_corr,
            cross_corr,
            cc_method,
            CCFDIR,
            maxlag,
            substack,
            substack_method,
        )
    end
    return nothing
end

function CC(
    cores::Int64,
    threads::Int64,
    flag::Bool,
    precompile::Bool,
    FFTDIR::String,
    CCFDIR::String,
    select_time::Bool,
    TIMETXT::String,
    select_sta::String,
    STATXT::String,
    SRTXT::String,
    channel_regular::String,
    comp::String,
    cc_method::String,
    auto_corr::Bool,
    cross_corr::Bool,
    maxlag::Float64,
    substack::Bool,
    substack_method::Vector{String},
)
    ### a. parameters --> Dict-Parameters
    Parameters = Dict(
        "cores" => cores,
        "threads" => threads,
        "flag" => flag,
        "precompile" => precompile,
        "FFTDIR" => FFTDIR,
        "CCFDIR" => CCFDIR,
        "select_time" => select_time,
        "TIMETXT" => TIMETXT,
        "select_sta" => select_sta,
        "STATXT" => STATXT,
        "SRTXT" => SRTXT,
        "channel_regular" => channel_regular,
        "comp" => comp,
        "cc_method" => cc_method,
        "auto_corr" => auto_corr,
        "cross_corr" => cross_corr,
        "maxlag" => maxlag,
        "substack" => substack,
        "substack_method" => substack_method,
    )

    ### b. get_alltime_path
    alltime_path = get_alltime_path_cc(FFTDIR, select_time, TIMETXT)
    N_time_chunks = size(alltime_path, 1)

    ### c. precompile
    if precompile == true
        println("precompile begin...")
        tALL_pre = @elapsed begin
            pmap(main_cc, [alltime_path[1]], [Parameters])
            toremove = glob("*/*/*/*/*/*", CCFDIR)
            rm.(toremove)
        end
        println("      precompile took $tALL_pre seconds")
        println("precompile end...")
    end

    ### d. run
    println("\n\n******************* START ******************")
    println(
        "Beginning XCORR with $(cores) processors and each core with $(threads) threads... || N_time_chunks: $(N_time_chunks)",
    )
    tALL = @elapsed begin
        pmap(main_cc, alltime_path, fill(Parameters, N_time_chunks))
    end
    println("XCORR took $tALL seconds")
    println("******************* END ******************\n\n")
end

max_stars_repo_name: OUCyf/NoiseCC

max_issues_repo_name: OUCyf/NoiseCC

max_forks_repo_name: OUCyf/NoiseCC


text:
###################################################
# METODO DO GRADIENTE, COMO DESCRITO NO SLIDE 29 DE
#   https://leonardosecchin.github.io/files/otim1/4.1.Metodos_descida_gerais.pdf
# Autor: Leonardo D. Secchin
# Data : 31/01/2021
#
# Exemplos de uso:
#   gradiente(nlp);
#   gradiente(nlp, x0=[1;1]);
#   gradiente(nlp, eps=1e-8);
#   x, f, gradnorm, iter, status = gradiente(nlp, x0=[1;1], eps=1e-8, maxiter=2000, saidas=true);
# onde 'nlp' é a estrutura MathOptNLPModel do problema.
###################################################

# carrega pacotes necessários
using JuMP, NLPModels, NLPModelsJuMP, Printf, LinearAlgebra

# FUNÇÃO PRINCIPAL
function gradiente(nlp; x0=nothing, eps=1.0e-6, maxiter=1000, saidas=true)

    # DADOS DE ENTRADA
    # nlp     : estrutura MathOptNLPModel do problema
    # x0      : ponto inicial (opcional)
    # eps     : precisão para convergencia (opcional, padrão=1.0e-6)
    # maxiter : número máximo de iterações (opcional, padrão=1000)
    # saidas  : mensagens na tela?

    # DADOS DE SAÍDA
    # x       : último iterando
    # f       : f(x)
    # gradnorm: norma do infinito do gradiente de f
    # iter    : número de iterações
    # status  : 0=sucesso, 1=falha

    iter   = 0
    status = 1

    # parâmetro busca linear inexata (Armijo)
    eta = 0.5

    # define ponto inicial caso não fornecido
    if x0 == nothing
        # captura ponto inicial da estrutura nlp
        if nlp.meta.x0 != nothing
            x0 = nlp.meta.x0
        else
            # caso não exista, seta x0=origem
            x0 = zeros(Float64,nlp.meta.nvar)
        end
    end

    # f, gradiente e norma
    x        = float(x0)
    f        = obj(nlp, x)
    gradf    = grad(nlp, x)
    gradnorm = norm(gradf, Inf)

    # imprimi cabeçalho saídas
    if saidas
        @printf("\nit     \tf         |grad|\n==========================")
        @printf("\n%d\t%8.2e  %8.2e", iter, f, gradnorm)
    end

    while gradnorm > eps && iter < maxiter
        # direção descida
        d = -gradf

        # retorna novo iterando após busca linear (Armijo)
        x, f = armijo(nlp, x, f, gradf, d, eta)

        # atualiza dados do iterando
        gradf    = grad(nlp, x)
        gradnorm = norm(gradf, Inf)
        iter += 1

        # imprimi iteração corrente
        if saidas
            if mod(iter,20) == 0
                @printf("\n\nit     \tf         |grad|\n==========================")
            end
            @printf("\n%d\t%8.2e  %8.2e", iter, f, gradnorm)
        end
    end

    # status de sucesso
    if gradnorm <= eps
        status = 0
    end

    # saída
    if saidas
        println("\n\n*******************************")
        if status == 0
            println("PROBLEMA RESOLVIDO COM SUCESSO!")
        else
            println("FALHA NA RESOLUÇÃO.")
        end
        println("*******************************\n")
    end

    return x, f, gradnorm, iter, status
end

# BUSCA LINEAR
function armijo(nlp, x, f, gradf, d, eta)

    # calcula gradf' * d
    gtd = gradf' * d

    # passo inicial
    t = 1

    # Armijo com backtracking
    xnew = x + t*d
    fnew = obj(nlp, xnew)

    while t > 1e-20 && fnew > f + t*eta*gtd
        # divide t por 10
        t /= 10.0

        # nova tentativa
        xnew = x + t*d
        fnew = obj(nlp, xnew)
    end

    # retorna novo iterando
    return xnew, fnew
end

max_stars_repo_name: leonardosecchin/temp

max_issues_repo_name: leonardosecchin/temp

max_forks_repo_name: leonardosecchin/temp


text:
import RvSpectMLBase
using Test

#=
@testset "Spectra Timeseries Common Wavelengths" begin
    # TODO Write test for Spectra Timeseries Common Wavelengths type
    # Won't actually test it here, since code to make this from a timeseries is not in Base package
    using RvSpectMLBase.TheoreticalInstrument
    using DataFrames
    num_pixels = 1000
    num_orders = 10
    num_obs = 3
    lambda2d = collect(reshape(range(5000.0,stop=6000.0,length=num_pixels*num_orders),num_pixels,num_orders))
    inst = TheoreticalInstrument2D(lambda2d)
    times = sort(rand(num_obs))
    rvs = range(0.0, stop=10.0, length=length(times))
    df = DataFrame(:lambda=>[5501], :weight=>[0.5])
    @test_nowarn TheoreticalInstrument.generate_spectra_timeseries(times,df,inst, rvs)
    spectra = TheoreticalInstrument.generate_spectra_timeseries(times,df,inst, rvs)
    clt = spectra  # Is this right type?
    chunk_map = map(i->make_grid_for_chunk(clt,i,oversample_factor=oversample_factor), 1:num_chunks(clt))
    # need to figure out dimensions for arguments below
    @test_nowarn stscw = SpectralTimeSeriesCommonWavelengths(λ, flux, var, chunk_map, TheoreticalInstrument1D() )
end
=#

max_stars_repo_name: RvSpectML/RvSpectMLTemplate

max_issues_repo_name: RvSpectML/RvSpectMLTemplate

max_forks_repo_name: RvSpectML/RvSpectMLTemplate


text:
"""
ABM is a multi agent based simulation engine.

Motivations:
- Julia allow easily to parallelize the code also at simulation code.
- MASON does not provide double buffering for agents memories and positions.
- Julia provides lots of libraries for scientific computing.
- Julia data visualization.
"""
module ABM
using Distributed
using BenchmarkTools
using Compat
using Random
using Distributions
using UUIDs
using DataStructures
using LightGraphs
using DistributedArrays

export Agent
export Patch

export Position
export Real2D,
        toroidal, distance

export Location
export Field
export Field2D,
        setObjectLocation!,
        getObjectLocation, getObjectsAtLocation, getNeighborsWithinDistance, getAllObjects,
        numObjectsAtLocation,
        remove!, clean!
export Network
export Bounds,
        checkBoundCircle

export Priority
export Schedule,
        scheduleOnce!, scheduleRepeating!,
        stop!,step!

export SimState,
        addfield!, update!, random

export simulate! #, @simulate

include("priority.jl")
include("agent.jl")
include("patch.jl")
include("field/position.jl")
include("field/location.jl")

include("field/real2d.jl")
include("field/field.jl")
include("field/field2d.jl")
include("field/bounds.jl")
include("simstate.jl")
include("schedule.jl")

#TODO Add to schedule the interval for scheduling
#TODO Agent function is bad to pass in the method also the agent other methods?
#TODO More test for schedule and field2d

function simulate!(schedule::Schedule, nsteps::Int64)
    while schedule.steps < nsteps
        println("[",schedule.steps,"] time: ",schedule.time)
        step!(schedule)

    end
end

end # module


text:
const tol = 10*eps()
const days_per_year = 365.2422
const global_base_date = Date(2000,1,1)
const global_base_date_as_day = convert(Dates.Day, global_base_date)

"""
    years_between(a::Union{DateTime,Date}, b::Union{DateTime,Date})
The number of years between two dates. This is returned as a scalar and
assumes 365.2422 days per year.

### Inputs
* `a` - The end date
* `b` - The start date.
### Returns
* A `Real`.
"""
function years_between(a::Union{DateTime,Date}, b::Union{DateTime,Date})
    return (Dates.days(a) -Dates.days(b))/ days_per_year
end

"""
    years_from_global_base(a::Union{DateTime,Date})
The number of years (calculated by the `years_between` function) between a given date
and the 1st of January 2000.

### Inputs
* `date` - A date.
### Returns
* A `Real`.
"""
function years_from_global_base(date::Union{DateTime,Date})
    return years_between(date, global_base_date)
end

"""
    period_length(a::Dates.DatePeriod, base::Date = global_base_date)

Period length is designed to convert `TimePeriod` objects to a float in a consistent way to `years_from_global_base`.
So effectively the years_between method is calculated with start and end dates being those
at the start and end of a `Dates.DatePeriod`. This is slightly complicated because
a period like `Month(3)` might have slightly different numbers of total days depending
on when in the year it is. So a base date has to be input. The period is then
measured starting from this base date.

### Inputs
* `period` - A period.
* `base` - A date from which the period will be measured from.
### Returns
* A `Real`.
"""
function period_length(period::Dates.DatePeriod, base::Date = global_base_date)
    return years_between(base+period, base)
end


text:
include("../MLPClassifier/MLPClassifier.jl")
include("mutual-information.jl")
include("../datasets/dataset_reader.jl")
include("../datasets/Drug-Classification/dataset1.jl")
include("../datasets/Income-Classification/dataset2.jl")
include("../datasets/Wine-Classification/dataset3.jl")
include("../datasets/Drug-Consumism-Classification/dataset4.jl")
include("../datasets/Heart-Disease-Classification/dataset5.jl")
include("../datasets/Demo/datasetDemo.jl")
include("../datasets/Demo-overlap/datasetDemoOverlap.jl")
using Plots, Statistics, LaTeXStrings
using Flux: onehot

struct LinearizedLayer
    vector::Vector
    nbins::Int
    function LinearizedLayer(vector, nbins::Int)
        new(vector, nbins)
    end
    function LinearizedLayer(vector, nfeatures::Int, neuron_nbins::Int, output_nbins::Int)
        hists = Vector{Vector{Int}}(undef, nfeatures)
        columns = separate_columns(vector)
        for feature in 1:nfeatures
            feature_values = columns[feature]
            min = minimum(feature_values)
            max = maximum(feature_values)
            hists[feature] = create_bin_vector(feature_values, neuron_nbins, min, max)
        end
        bounds = repeat([neuron_nbins], nfeatures)
        indices = separate_columns(hists)
        combination_values = get_bin_from_indices.(indices, [bounds])
        max_value = foldr(*, bounds)
        linearized = create_bin_vector(combination_values, output_nbins, 1, max_value)
        new(linearized, output_nbins)
    end
    function LinearizedLayer(vector, neuron_nbins::Int, output_nbins::Int, lower_lim, upper_lim)
        columns = separate_columns(vector)
        hists = create_bin_histogram.(columns, [neuron_nbins], [lower_lim], [upper_lim])
        bounds = (get_weights_length.(hists)...,)
        indices = [getindex.(get_bin.(hists, v),1) for v in vector]
        indices = [correct_upper_bounds(ind, bounds) for ind in indices]
        combination_values = get_bin_from_indices.(indices, [bounds])
        #max_value = max(prod(bounds), typemax(Int64)) # foldr(*, bounds) == prod(bounds)
        max_value = map(n -> UInt128(n), bounds) |> prod
        max_value = max_value > typemax(Int64) ? typemax(Int64) : max_value
        linearized = create_bin_vector(combination_values, output_nbins, 1, max_value)
        new(linearized, output_nbins)
    end
end

# Check if the number of nodes explodes
# (map(n -> UInt128(n), nbins*ones(Int, num_nodes)) |> prod) > typemax(Int64)

# lims::::Tuple{T,T} where T <: Number
function layer_information_plane_position_MarkovChain_t(model::MLPClassifier, input, layer::Int, lims, X_flatten, Y_flatten, I_XY, nbins::Int=10)
    if length(model.net) >= layer && layer > 0
        subnet = model.net[1:layer]
        local subnet_output
        if length(model.net) == layer
            subnet_output = get_label_position.(predict.([model], input), [model.labels])
            nbins_subnet_output = length(model.labels)
            subnet_output = LinearizedLayer(subnet_output, nbins_subnet_output)
        else
            subnet_output = subnet.(input)
            nbins_subnet_output = nbins*length(subnet_output[1])#=nbins^length(subnet_output[1])=#
            subnet_output = LinearizedLayer(subnet_output, nbins, nbins_subnet_output, lims[1], lims[2])
        end
        I_XT = mutual_information_unit(X_flatten.vector, subnet_output.vector, (0, X_flatten.nbins), (0, subnet_output.nbins))
        I_TY = mutual_information(calculate_markov_chain_probability(Y_flatten.vector, X_flatten.vector, subnet_output.vector, (X_flatten.nbins, Y_flatten.nbins, subnet_output.nbins))...)
        (I_XT, I_TY/I_XY)
    else
        error("Invalid layer number.")
    end
end

get_label_position(value, label_vector) = findfirst(x -> x == value, label_vector)

#=
function calculate_joint_probability(X::Vector{<:Number}, Y::Vector{<:Number}, lim_X, lim_Y)
    min_X, max_X = lim_X
    min_Y, max_Y = lim_Y
    range_X = min_X:max_X
    range_Y = min_Y:max_Y
    p_x = fit(Histogram, (X), range_X, closed=:right)
    p_y = fit(Histogram, (Y), range_Y, closed=:right)
    p_xy = fit(Histogram, (X, Y), (range_X, range_Y), closed=:right)
    p_X = p_x.weights/sum(p_x.weights)
    p_Y = p_y.weights/sum(p_y.weights)
    p_XY = p_xy.weights/sum(p_xy.weights)
    return p_X, p_Y, p_XY
end
=#

function calculate_markov_chain_probability(Y::Vector{<:Number}, X::Vector{<:Number}, T::Vector{<:Number})
    p_t = fit(Histogram, (T))
    p_y = fit(Histogram, (Y))
    p_x = fit(Histogram, (X))
    p_xy = fit(Histogram, (X, Y))
    p_tx = fit(Histogram, (T, X))
    p_ty = zeros(Float64, (length(p_t.weights), length(p_y.weights)))
    for t = 1:length(p_t.weights)
        for y = 1:length(p_y.weights)
            for x = 1:length(p_x.weights)
                p1 = get_bin_probability(p_xy, (x,y)...)
                p2 = get_bin_probability(p_tx, (t,x)...)
                p3 = get_bin_probability(p_x, (x)...)
                if p2 == 0 && p3 == 0
                    p_ty[t,y] += 0
                else
                    p_ty[t,y] += p1*(p2/p3)
                end
            end
        end
    end
    p_T = p_t.weights/sum(p_t.weights)
    p_Y = p_y.weights/sum(p_y.weights)
    return p_T, p_Y, p_ty
end

function calculate_markov_chain_probability(Y::Vector{<:Number}, X::Vector{<:Number}, T::Vector{<:Number}, nbins)
    p_t = fit(Histogram, (T), nbins=nbins[3])
    p_y = fit(Histogram, (Y), nbins=nbins[2])
    p_x = fit(Histogram, (X), nbins=nbins[1])
    p_xy = fit(Histogram, (X, Y), nbins=(nbins[1:2]...,))
    p_tx = fit(Histogram, (T, X), nbins=(nbins[3], nbins[1]))
    p_ty = zeros(Float64, (length(p_t.weights), length(p_y.weights)))
    for t = 1:length(p_t.weights)
        for y = 1:length(p_y.weights)
            for x = 1:length(p_x.weights)
                p1 = get_bin_probability(p_xy, (x,y)...)
                p2 = get_bin_probability(p_tx, (t,x)...)
                p3 = get_bin_probability(p_x, (x)...)
                if p2 == 0 && p3 == 0
                    p_ty[t,y] += 0
                else
                    p_ty[t,y] += p1*(p2/p3)
                end
            end
        end
    end
    p_T = p_t.weights/sum(p_t.weights)
    p_Y = p_y.weights/sum(p_y.weights)
    return p_T, p_Y, p_ty
end

function get_corrected_bin(h::Histogram, val)
    bin = searchsortedfirst.(h.edges, val)
    bounds = size(h.weights)
    for i = 1:length(bin)
        if bin[i] >= bound[i]
            bin[i] = bound[i]
        end
    end
    bin
end

function evaluate_information_plane(model::MLPClassifier, input, output, nbins::Int=10)
    layers_number = length(model.net)
    X_flatten = LinearizedLayer(input, model.nfeatures, nbins, nbins*model.nfeatures#=nbins^model.nfeatures=#)
    Y_flatten = LinearizedLayer(get_label_position.(getindex.(output,1), [model.labels]), length(model.labels))
    I_XY = mutual_information_unit(X_flatten.vector, Y_flatten.vector, (0, X_flatten.nbins), (0, Y_flatten.nbins))
    activation_lims = (-1, 1) # Limits of activation function tanh
    output_lims = (0, 1) # Limits of activation function softmax
    IP_points = layer_information_plane_position_MarkovChain_t.([model], [input], 1:(layers_number-2), [activation_lims], [X_flatten], [Y_flatten], [I_XY], [nbins])
    push!(IP_points, layer_information_plane_position_MarkovChain_t(model, input, layers_number, (minimum(Y_flatten.vector), maximum(Y_flatten.vector)), X_flatten, Y_flatten, I_XY, nbins))
    IP_points
end

function information_bottleneck(β::Float64, model::MLPClassifier, input, output;  nbins::Int=10)
    net_output = flat_vector(model(input), nbins = nbins)
    X_flatten = flat_vector(input, nbins = nbins)
    Y_flatten = flat_vector(output, nbins = nbins)
    I_XZ = mutual_information([X_flatten], [net_output])
    I_ZY = mutual_information([net_output], [Y_flatten])
    I_XZ - β*I_ZY
end

function flat_vector1(vector)
    columns = separate_columns(vector)
    h = fit(Histogram, (columns...,))
    bounds = size(h.weights)
    indexes = [searchsortedfirst.(h.edges, v) for v in vector]
    correct_upper_bounds!.(indexes, [bounds])
    [LinearIndices(h.weights)[ind...] for ind in indexes]
end

#=function flat_vector(vector)
    columns = separate_columns(vector)
    hists = fit.(Histogram, columns)
    bounds = (get_weights_length.(hists)...,)
    indices = [getindex.(get_bin.(hists, v),1) for v in vector]
    indices = [correct_upper_bounds(ind, bounds) for ind in indices]
    #linear_indexes = LinearIndices(bounds)
    get_bin_from_indices.(indices, [bounds])
end=#

include("information_plane_plot.jl")

#=
# Test
model = MLPClassifier(dataset1_nfeatures, dataset1_labels, [15,5,10,5])
x_train, y_train, x_test, y_test = get_dataset1("../datasets/Drug-Classification")
x_train, y_train, x_test, y_test = get_dataset2("../datasets/Income-Classification")
net = MLPClassifier(dataset2_nfeatures, dataset2_labels, [26,13])

plot_IP_average(dataset2_nfeatures, dataset2_labels, [26,13], 3, get_train_dataset2, "../datasets/Income-Classification", 100, 2, "IP-avg-10n-100ep-v2.png")
plot_IP_average(dataset2_nfeatures, dataset2_labels, [26,13], 3, data, 100, 2, "IP-avg-10n-100ep-v2.png")

#fit!(model, x_train, y_train ; max_epochs=20)
plot_IP(model, x_train, y_train)
println("-> Error: ", error_rate(zip(x_test, y_test), model))
fit!(model, x_train, y_train ; max_epochs=1000)
plot_IP(model, x_train, y_train, imagename="IP2.png")
println("-> Error: ", error_rate(zip(x_test, y_test), model))
fit!(model, x_train, y_train ; max_epochs=2000)
plot_IP(model, x_train, y_train, imagename="IP3.png")
println("-> Error: ", error_rate(zip(x_test, y_test), model))
=#

#=
x_train, y_train, x_test, y_test = get_datasetDemo("../datasets/Demo")
net = MLPClassifier(datasetDemo_nfeatures, datasetDemo_labels, [3,2])
fit!(net, x_train, y_train ; max_epochs=10)
println("-> Error: ", error_rate(zip(x_test, y_test), net))
data = zip(x_train, y_train)
plot_IP_average(datasetDemo_nfeatures, datasetDemo_labels, topology, number_nets, data, epochs, step_len, "IP-dtDemo-MLP-fixed.png")
=#

#=
topology = [3,2]
number_nets = 10
epochs = 100
step_len = 10
x_train, y_train, x_test, y_test = get_datasetDemoOverlap("../datasets/Demo-overlap")
data = zip(x_train, y_train)
net = MLPClassifier(datasetOverlap_nfeatures, datasetOverlap_labels, topology)
plot_IP_average(datasetOverlap_nfeatures, datasetOverlap_labels, topology, number_nets, data, epochs, step_len, "IP-dtOverlap-MLP-fixed.png")
plot_IP_average(datasetOverlap_nfeatures, datasetOverlap_labels, topology, number_nets, data, epochs, step_len, "IP-dtOverlap-MLP-fixed2.png")
=#

max_stars_repo_name: Lucas-Muniz/MNEAT

max_issues_repo_name: Lucas-Muniz/MNEAT

max_forks_repo_name: Lucas-Muniz/MNEAT


text:
"""
# Notes
- [1] is modified for hexa-x configuration; see [2].
# References
[1] T. Lee, M. Leok, and N. H. McClamroch,
“Geometric Tracking Control of a Quadrotor UAV on SE(3),”
in 49th IEEE Conference on Decision and Control (CDC), Atlanta, GA, Dec. 2010, pp. 5420–5425. doi: 10.1109/CDC.2010.5717652.
[2] PX4 Airframes Reference, https://docs.px4.io/master/en/airframes/airframe_reference.html.
[3] M. C. Achtelik, K. M. Doth, D. Gurdan, and J. Stumpf, “Design of a multi rotor MAV with regard to efficiency, dynamics and redundancy,” AIAA Guid. Navig. Control Conf. 2014, no. August, pp. 1–17, 2012, doi: 10.2514/6.2012-4779.
[4] https://github.com/fdcl-ftc/fault-tolerant-control/issues/62
# Variables
u ∈ R^6: rotor forces
"""
Base.@kwdef struct LeeHexacopter <: Hexacopter
    J = diagm([0.0820, 0.0845, 0.1377])  # kg m^2
    l = 0.315  # m
    kM = 8.004e-4  # m
    m = 4.34  # kg
    g = 9.81  # m / s^2
    B = [ 1   1             1              1             1              1;
         -l   l         0.5*l         -0.5*l        -0.5*l          0.5*l;
          0   0 0.5*sqrt(3)*l -0.5*sqrt(3)*l 0.5*sqrt(3)*l -0.5*sqrt(3)*l;
        -kM  kM           -kM             kM            kM            -kM]
    # actuator limit
    dim_input = 6
    u_min = zeros(dim_input)
    u_max = 0.6371 * (m*g) * ones(dim_input)  # [3] and [4]
end

function saturate(multicopter::LeeHexacopter, u)
    @unpack u_min, u_max = multicopter
    u_saturated = clamp.(u, u_min, u_max)
end

function input_to_force_moment(multicopter::LeeHexacopter, u)
    @unpack B = multicopter
    ν = B * u
end

function airframe_reference(multicopter::LeeHexacopter)
    :hexa_x
end


text:
@testset "Create new app" begin

  testdir = pwd()
  using Pkg

  @safetestset "Do not autostart app" begin
    using Genie

    workdir = Base.Filesystem.mktempdir()

    Genie.newapp(workdir, autostart = false, testmode = true)

    @test true === true
  end;

  # cd(testdir)
  # Pkg.activate(".")

  # @safetestset "Autostart app" begin
  #   using Genie

  #   workdir = Base.Filesystem.mktempdir()

  #   Genie.newapp(workdir, autostart = true, testmode = true)

  #   @test true === true
  # end;

  cd(testdir)
  Pkg.activate(".")

  @safetestset "Microstack file structure" begin
    using Genie

    workdir = Base.Filesystem.mktempdir()

    Genie.newapp(workdir, autostart = false, testmode = true)

    @test sort(readdir(workdir)) == sort([".gitattributes", ".gitignore", "Manifest.toml", "Project.toml", "bin",
                                          "bootstrap.jl", "config", "genie.jl", "public", "routes.jl", "src"])
    @test readdir(joinpath(workdir, Genie.config.path_initializers)) == ["converters.jl", "logging.jl"]

    # TODO: add test for files in /src /config /public and /bin
  end;

  cd(testdir)
  Pkg.activate(".")

  @safetestset "DB support file structure" begin
    using Genie

    workdir = Base.Filesystem.mktempdir()

    Genie.newapp(workdir, autostart = false, dbsupport = true, testmode = true)

    @test sort(readdir(workdir)) == sort([".gitattributes", ".gitignore", "Manifest.toml", "Project.toml", "bin",
                                          "bootstrap.jl", "config", "db", "genie.jl", "public", "routes.jl", "src"])
    @test sort(readdir(joinpath(workdir, Genie.config.path_db))) == sort(["connection.yml", "migrations", "seeds"])
    @test sort(readdir(joinpath(workdir, Genie.config.path_initializers))) == sort(["converters.jl", "logging.jl", "searchlight.jl"])
  end;

  cd(testdir)
  Pkg.activate(".")

  @safetestset "MVC support file structure" begin
    using Genie

    workdir = Base.Filesystem.mktempdir()

    Genie.newapp(workdir, autostart = false, mvcsupport = true, testmode = true)

    @test sort(readdir(workdir)) == sort([".gitattributes", ".gitignore", "Manifest.toml", "Project.toml", "app", "bin", "bootstrap.jl", "config", "genie.jl", "public", "routes.jl", "src"])
    @test sort(readdir(joinpath(workdir, Genie.config.path_app))) == sort(["helpers", "layouts", "resources"])
    @test sort(readdir(joinpath(workdir, Genie.config.path_initializers))) == sort(["converters.jl", "logging.jl"])
  end;

  cd(testdir)
  Pkg.activate(".")

  @safetestset "New controller" begin
    using Genie

    workdir = Base.Filesystem.mktempdir()
    cd(workdir)
    Genie.newcontroller("Yazoo")

    @test isdir(joinpath(workdir, "app", "resources", "yazoo")) == true
    @test isfile(joinpath(workdir, "app", "resources", "yazoo", "YazooController.jl")) == true
  end;

  cd(testdir)
  Pkg.activate(".")

  @safetestset "New resource" begin
    using Genie

    workdir = Base.Filesystem.mktempdir()
    cd(workdir)
    Genie.newresource("Kazoo")

    @test isdir(joinpath(workdir, "app", "resources", "kazoo")) == true
    @test isfile(joinpath(workdir, "app", "resources", "kazoo", "KazooController.jl")) == true
  end;

  cd(testdir)
  Pkg.activate(".")

  @safetestset "New task" begin
    using Genie, Genie.Exceptions

    workdir = Base.Filesystem.mktempdir()
    cd(workdir)
    Genie.newtask("Vavoom")

    @test isdir(joinpath(workdir, "tasks")) == true
    @test isfile(joinpath(workdir, "tasks", "VavoomTask.jl")) == true
    @test_throws FileExistsException Genie.newtask("Vavoom")
  end;

  cd(testdir)
  Pkg.activate(".")

end;


text:
# Run all tests.
module TestRunner
  using FactCheck
  using Orchestra.System

  include("test_util.jl")
  include(joinpath("orchestra", "test_transformers.jl"))
  include(joinpath("julia", "test_decisiontree.jl"))
  include(joinpath("julia", "test_mlbase.jl"))
  include(joinpath("julia", "test_dimensionalityreduction.jl"))
  include(joinpath("orchestra", "test_ensemble.jl"))
  if LIB_SKL_AVAILABLE
    include(joinpath("python", "test_scikit_learn.jl"))
  else
    info("Skipping scikit-learn tests.")
  end
  if LIB_CRT_AVAILABLE
    include(joinpath("r", "test_caret.jl"))
  else
    info("Skipping CARET tests.")
  end
  include("test_system.jl")

  exitstatus()
end # module


text:
using Kinetic, ProgressMeter, Plots, LinearAlgebra, JLD2

cd(@__DIR__)
D = Dict{Symbol,Any}()
begin
    D[:matter] = "gas"
    D[:case] = "sod"
    D[:space] = "1d2f1v"
    D[:flux] = "kfvs"
    D[:collision] = "bgk"
    D[:nSpecies] = 1
    D[:interpOrder] = 2
    D[:limiter] = "vanleer"
    D[:boundary] = "fix"
    D[:cfl] = 0.8
    D[:maxTime] = 0.2

    D[:x0] = 0.0
    D[:x1] = 1.0
    D[:nx] = 100
    D[:pMeshType] = "uniform"
    D[:nxg] = 1

    D[:umin] = -5.0
    D[:umax] = 5.0
    D[:nu] = 100
    D[:vMeshType] = "rectangle"
    D[:nug] = 0

    D[:knudsen] = 1e-4
    D[:mach] = 0.0
    D[:prandtl] = 1.0
    D[:inK] = 2.0
    D[:omega] = 0.81
    D[:alphaRef] = 1.0
    D[:omegaRef] = 0.5
end

ks = SolverSet(D)
ctr, face = init_fvm(ks, ks.ps, :dynamic_array; structarray = true)

t = 0.0
dt = timestep(ks, ctr, t)
nt = Int(floor(ks.set.maxTime / dt))
res = zeros(3)
for iter = 1:nt
    reconstruct!(ks, ctr)
    evolve!(ks, ctr, face, dt)
    Kinetic.update!(ks, ctr, face, dt, res)
end

plot_line(ks, ctr)

sol0 = zeros(ks.ps.nx, 3)
for i = 1:ks.ps.nx
    sol0[i, :] .= ctr.prim[i]
    sol0[i, 3] = 1 / sol0[i, 3]
end

max_stars_repo_name: CSMMLab/Flowmachine

max_issues_repo_name: CSMMLab/Flowmachine

max_forks_repo_name: CSMMLab/Flowmachine


text:
using Revise
using AdFem
using LinearAlgebra
using PyPlot
using DelimitedFiles
using MAT

function calc_residual_and_jacobian(θ, u)
    E = θ[1:get_ngauss(mmesh)]
    nu = θ[get_ngauss(mmesh)+1:end]
    μ = E / (2*(1 + nu))
    λ = E*nu/((1 + nu)*(1 - 2*nu))

    ψ, J = neo_hookean(u, μ, λ, mmesh)
    B = compute_fem_source_term(zeros(get_ngauss(mmesh)), -0.5*ones(get_ngauss(mmesh)), mmesh)
    ψ = ψ - B
    J, ψ = impose_Dirichlet_boundary_conditions(J, ψ, bdnode, zeros(length(bdnode)))
    ψ, J
end

mmesh = Mesh(joinpath(PDATA, "twoholes_large.stl"))
mmesh = Mesh(mmesh.nodes * 10, mmesh.elems)

θ = placeholder([10.0 * ones(get_ngauss(mmesh)); 0.3 * ones(get_ngauss(mmesh))])
bdnode = bcnode((x,y)->(x<1e-5) || (x>0.49-1e-5), mmesh)
bdnode = [bdnode; bdnode .+ mmesh.ndof]
nr = newton_raphson_with_grad(calc_residual_and_jacobian, zeros(2mmesh.ndof), θ)

ADCME.options.newton_raphson.verbose = true
sess = Session(); init(sess)
x = run(sess, nr)

matwrite("fenics/fwd1.mat", Dict("x"=>x))

close("all")
figure(figsize = (15, 5))
subplot(121)
u = x[1:mmesh.nnode]
v = x[mmesh.ndof+1:mmesh.ndof + mmesh.nnode]
visualize_scalar_on_fem_points(sqrt.(u.^2+v.^2), mmesh)
subplot(122)
visualize_displacement(x, mmesh)
savefig("fenics/fwd1.png")


text:
module P2P
# # current release: v0.18
# #     next planned release: v0.20 Aug 27, 2021.
# #     Please do NOT download source code until v0.20 release.

# #     !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
# #     !!! source code will available v0.20 Aug 27, 2022 !!!
# #     !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

end