text:
using LinearAlgebra
using DiffEqBase
using LsqFit

abstract type AbstractInterface end
export AbstractInterface

abstract type AbstractBoundaryInterface <: AbstractInterface end
export AbstractBoundaryInterface

abstract type AbstractInternalInterface <: AbstractInterface end
export AbstractInternalInterface

abstract type AbstractReactiveInternalInterface <: AbstractInternalInterface end
export AbstractReactiveInternalInterface

struct EmptyInterface <: AbstractInterface end
export EmptyInterface

struct ReactiveInternalInterface{T,B,C,C2,N,Q<:AbstractReaction,X} <: AbstractReactiveInternalInterface
    domain1::T
    domain2::N
    reactions::Array{Q,1}
    veckinetics::X
    veckineticsinds::Array{Int64,1}
    rxnarray::B
    stoichmatrix::C
    Nrp1::C2
    Nrp2::C2
    A::Float64
    parameterindexes::Array{Int64,1}
    domaininds::Array{Int64,1}
    p::Array{Float64,1}
    reversibililty::Array{Bool,1}
end
function ReactiveInternalInterface(domain1,domain2,reactions,A)
    vectuple,vecinds,otherrxns,otherrxninds,posinds = getveckinetics(reactions)
    rxns = vcat(reactions[vecinds],reactions[otherrxninds])
    rxns = [ElementaryReaction(index=i,reactants=rxn.reactants,reactantinds=rxn.reactantinds,products=rxn.products,
        productinds=rxn.productinds,kinetics=rxn.kinetics,radicalchange=rxn.radicalchange,reversible=rxn.reversible,pairs=rxn.pairs) for (i,rxn) in enumerate(rxns)]
    rxnarray = getinterfacereactioninds(domain1,domain2,rxns)
    M,Nrp1,Nrp2 = getstoichmatrix(domain1,domain2,reactions)
    reversibility = Array{Bool,1}(getfield.(rxns,:reversible))
    return ReactiveInternalInterface(domain1,domain2,
            rxns,vectuple,posinds,rxnarray,M,Nrp1,Nrp2,A,[1,length(reactions)],
            [0,1],ones(length(rxns)),reversibility),ones(length(rxns))
end
export ReactiveInternalInterface

function getkfskrevs(ri::ReactiveInternalInterface,T1,T2,phi1,phi2,Gs1,Gs2,cstot::Array{Q,1}) where {Q}
    kfs = getkfs(ri,T1,0.0,0.0,Array{Q,1}(),ri.A,phi1)
    Kc = getKc.(ri.reactions,ri.domain1.phase,ri.domain2.phase,Ref(Gs1),Ref(Gs2),T1,phi1)
    krevs = kfs./Kc
    return kfs,krevs
end

function evaluate(ri::ReactiveInternalInterface,dydt,domains,T1,T2,phi1,phi2,Gs1,Gs2,cstot,p::W) where {W<:DiffEqBase.NullParameters}
    kfs,krevs = getkfskrevs(ri,T1,T2,phi1,phi2,Gs1,Gs2,cstot)
    addreactionratecontributions!(dydt,ri.rxnarray,cstot,kfs,krevs,ri.A)
end

function evaluate(ri::ReactiveInternalInterface,dydt,domains,T1,T2,phi1,phi2,Gs1,Gs2,cstot,p)
    kfs,krevs = getkfskrevs(ri,T1,T2,phi1,phi2,Gs1,Gs2,cstot)
    addreactionratecontributions!(dydt,ri.rxnarray,cstot,kfs.*p[ri.parameterindexes[1]:ri.parameterindexes[2]],krevs.*p[ri.parameterindexes[1]:ri.parameterindexes[2]],ri.A)
end
export evaluate

struct ReactiveInternalInterfaceConstantTPhi{J,N,B,B2,B3,C,C2,Q<:AbstractReaction} <: AbstractReactiveInternalInterface
    domain1::J
    domain2::N
    reactions::Array{Q,1}
    rxnarray::B
    stoichmatrix::C
    Nrp1::C2
    Nrp2::C2
    kfs::B2
    krevs::B3
    T::Float64
    A::Float64
    parameterindexes::Array{Int64,1}
    domaininds::Array{Int64,1}
    p::Array{Float64,1}
    reversibility::Array{Bool,1}
end
function ReactiveInternalInterfaceConstantTPhi(domain1,domain2,reactions,T,A,phi=0.0)
    @assert domain1.T == domain2.T
    reactions = upgradekinetics(reactions,domain1,domain2)
    rxnarray = getinterfacereactioninds(domain1,domain2,reactions)
    kfs = getkf.(reactions,nothing,T,0.0,0.0,Ref([]),A,phi)
    Kc = getKc.(reactions,domain1.phase,domain2.phase,Ref(domain1.Gs),Ref(domain2.Gs),T,phi)
    krevs = kfs./Kc
    M,Nrp1,Nrp2 = getstoichmatrix(domain1,domain2,reactions)
    reversibility = Array{Bool,1}(getfield.(reactions,:reversible))
    if isa(reactions,Vector{Any})
        reactions = convert(Vector{ElementaryReaction},reactions)
    end
    if isa(kfs,Vector{Any})
        kfs = convert(Vector{Float64},kfs)
    end
    return ReactiveInternalInterfaceConstantTPhi(domain1,domain2,reactions,
            rxnarray,M,Nrp1,Nrp2,kfs,krevs,T,A,[1,length(reactions)],
            [0,1],kfs[1:end],reversibility),kfs[1:end]
end
export ReactiveInternalInterfaceConstantTPhi

function getkfskrevs(ri::ReactiveInternalInterfaceConstantTPhi,T1,T2,phi1,phi2,Gs1,Gs2,cstot)
    return ri.kfs,ri.krevs
end

function evaluate(ri::ReactiveInternalInterfaceConstantTPhi,dydt,domains,T1,T2,phi1,phi2,Gs1,Gs2,cstot,p::W) where {W<:DiffEqBase.NullParameters}
    addreactionratecontributions!(dydt,ri.rxnarray,cstot,ri.kfs,ri.krevs,ri.A)
end

function evaluate(ri::ReactiveInternalInterfaceConstantTPhi,dydt,domains,T1,T2,phi1,phi2,Gs1,Gs2,cstot,p)
    if p[ri.parameterindexes[1]:ri.parameterindexes[2]] == ri.kfs
        kfs = ri.kfs
    else
        kfs = p[ri.parameterindexes[1]:ri.parameterindexes[2]]
    end
    if length(Gs1) == 0 || length(Gs2) == 0 || (all(Gs1 .== ri.domain1.Gs) && all(Gs2 .== ri.domain2.Gs))
        krevs = ri.krevs
    else
        Kc = getKcs(ri,T1,Gs1,Gs2)
        krevs = kfs./Kc
    end
    addreactionratecontributions!(dydt,ri.rxnarray,cstot,kfs,krevs,ri.A)
end
export evaluate

"""
construct the stochiometric matrix for the reactions crossing both domains and the reaction molecule # change
"""
function getstoichmatrix(domain1,domain2,rxns)
    M = spzeros(length(rxns),length(domain1.phase.species)+length(domain2.phase.species))
    Nrp1 = zeros(length(rxns))
    Nrp2 = zeros(length(rxns))
    N1 = length(domain1.phase.species)
    spcs1 = domain1.phase.species
    spcs2 = domain2.phase.species
    for (i,rxn) in enumerate(rxns)
        Nrp1[i] = Float64(length([x for x in rxn.products if x in spcs1]) - length([x for x in rxn.reactants if x in spcs1]))
        Nrp2[i] = Float64(length([x for x in rxn.products if x in spcs2]) - length([x for x in rxn.reactants if x in spcs2]))
        for (j,r) in enumerate(rxn.reactants)
            isfirst = true
            ind = findfirst(isequal(r),domain1.phase.species)
            if ind === nothing
                isfirst = false
                ind = findfirst(isequal(r),domain2.phase.species)
            end
            M[i,isfirst ? ind : ind+N1] += 1
        end
        for (j,r) in enumerate(rxn.products)
            isfirst = true
            ind = findfirst(isequal(r),domain1.phase.species)
            if ind === nothing
                isfirst = false
                ind = findfirst(isequal(r),domain2.phase.species)
            end
            M[i,isfirst ? ind : ind+N1] -= 1
        end
    end
    return M,Nrp1,Nrp2
end

function getinterfacereactioninds(domain1,domain2,reactions)
    indices = zeros(Int64,(6,length(reactions)))
    N1 = length(domain1.phase.species)
    for (i,rxn) in enumerate(reactions)
        for (j,r) in enumerate(rxn.reactants)
            isfirst = true
            ind = findfirst(isequal(r),domain1.phase.species)
            if ind === nothing
                isfirst = false
                ind = findfirst(isequal(r),domain2.phase.species)
            end
            indices[j,i] = isfirst ? ind : ind+N1
        end
        for (j,r) in enumerate(rxn.products)
            isfirst = true
            ind = findfirst(isequal(r),domain1.phase.species)
            if ind === nothing
                isfirst = false
                ind = findfirst(isequal(r),domain2.phase.species)
            end
            indices[j+3,i] = isfirst ? ind : ind+N1
        end
    end
    return indices
end

function upgradekinetics(rxns,domain1,domain2)
    domain1surf = hasproperty(domain1.phase,:sitedensity)
    domain2surf = hasproperty(domain2.phase,:sitedensity)
    @assert !(domain1surf && domain2surf)
    if domain1surf
        surfdomain = domain1
    elseif domain2surf
        surfdomain = domain2
    end
    newrxns = Array{ElementaryReaction,1}(undef,length(rxns))
    for (i,rxn) in enumerate(rxns)
        if isa(rxn.kinetics,StickingCoefficient)
            spc = [spc for spc in rxn.reactants if !(spc in surfdomain.phase.species)]
            @assert length(spc) == 1
            kin = stickingcoefficient2arrhenius(rxn.kinetics,surfdomain.phase.sitedensity,length(rxn.reactants)-1,spc[1].molecularweight)
            newrxns[i] = ElementaryReaction(index=rxn.index,reactants=rxn.reactants,reactantinds=rxn.reactantinds,products=rxn.products,
                productinds=rxn.productinds,kinetics=kin,radicalchange=rxn.radicalchange,reversible=rxn.reversible,pairs=rxn.pairs)
        else
            newrxns[i] = rxn
        end
    end
    return [rxn for rxn in newrxns]
end

function stickingcoefficient2arrhenius(sc,sitedensity,N,mw;Tmin=300.0,Tmax=2000.0)
    mass = mw/Na
    ksc(T) = sc(T)/sitedensity^N*sqrt(kB*T/(2.0*pi*mass))
    Ts = Array{Float64,1}(Tmin:10:Tmax);
    kscvals = ksc.(Ts)
    k(T,p) = log.(abs(p[1]).*T.^p[2].*exp.(-p[3]./(R.*T)))
    p0 = [sc.A/sitedensity*sqrt(kB*1000.0/(2.0*pi*mass)),0.5,sc.Ea]
    fit = curve_fit(k,Ts,log.(kscvals),p0;x_tol=1e-18)
    @assert fit.converged
    p = fit.param
    p[1] = abs(p[1])
    return Arrhenius(;A=p[1],n=p[2],Ea=p[3])
end

struct Inlet{Q<:Real,S,V<:AbstractArray,U<:Real,X<:Real,FF<:Function} <: AbstractBoundaryInterface
    domain::S
    y::V
    F::FF
    T::U
    P::X
    H::Q
end

function Inlet(domain::V,conddict::Dict{String,X},F::FF) where {V,X<:Real,B<:Real,FF<:Function}
    y = makespcsvector(domain.phase,conddict)
    T = conddict["T"]
    P = conddict["P"]
    yout = y./sum(y)
    H = dot(getEnthalpy.(getfield.(domain.phase.species,:thermo),T),yout)
    return Inlet(domain,yout,F,T,P,H)
end

export Inlet

struct Outlet{V,FF<:Function} <: AbstractBoundaryInterface
    domain::V
    F::FF
end
export Outlet


text:
function check_for_new_closed(client, config, bot_path)
    new_closed = glob(joinpath(bot_path[2:end], "temp/new_closed*.json"))
    if length(new_closed) > 0
        yr, m = split(string(now()), "-")[1:2]
        current_closed_orders_path = joinpath(bot_path, "history/" * yr * "/" * m * ".json")
        !isdir(dirname(current_closed_orders_path)) ? mkpath(dirname(current_closed_orders_path)) : ""

        closed_orders = isfile(current_closed_orders_path) ? JSON.parsefile(current_closed_orders_path) : Dict()
        !haskey(closed_orders, string(today())) ? closed_orders[string(today())] = [] : ""
        #
        ## If testing is active
        test_bot = config["testing"] == "1" ? true : false
        if test_bot
            #
            ## Profit is estimated using the assumption that BNB is used for fee payments
            ### the associated BNB value is subtracted from the raw profit
            for new in new_closed
                tmp = JSON.parsefile(new)
                if tmp["type"] == "normal"
                    paid = tmp["buy_order"]["cummulativeQuoteQty"]
                    got = tmp["order"]["cummulativeQuoteQty"]
                elseif tmp["type"] == "initial"
                    price = typeof(config["init_state"]["price"]) != Float64 ? parse(Float64, config["init_state"]["price"]) : config["init_state"]["price"]
                    got = typeof(tmp["order"]["cummulativeQuoteQty"]) != Float64 ? parse(Float64, tmp["order"]["cummulativeQuoteQty"]) : tmp["order"]["cummulativeQuoteQty"]
                    paid = got * price/parse(Float64, tmp["order"]["price"])
                end
                #
                ## Fees paid using BNB assumed (0.075% per trade)
                fees = (paid + got) * 0.00075
                profit = got - paid - fees
                tmp["profit"] = profit
                # Calculate thes profit and profit percent
                if !haskey(tmp, "current_volume")
                    current_volume = typeof(config["total_value"]) == Float64 ? config["total_value"] : parse(Float64, config["total_value"])
                    if length(config["extra_volume"]) > 0
                        current_volume += sum(config["extra_volume"])
                    end
                    current_volume += config["thes_volume"]
                else
                    current_volume = parse(Float64, config["total_value"])
                end
                tmp["profit_percent"] = tmp["profit"] * 100 / current_volume
                tmp["profit_thes"] = tmp["profit"] * config["thes_factor"]

                push!(closed_orders[string(today())], tmp)

                open(current_closed_orders_path, "w") do f
                    JSON.print(f, closed_orders, 4)
                end
                open(new, "w") do f
                    JSON.print(f, tmp, 4)
                end

                if tmp["type"] == "normal"
                    mv(new, joinpath(bot_path, "closed_orders/" * basename(new)))
                elseif tmp["type"] == "initial"
                    mv(new, joinpath(bot_path, "closed_init_orders/" * basename(new)))
                end
            end
        else
            #
            ## Real trades
            Pinguu.anaylse_closed_orders(client, config, bot_path; limit=2)
        end
        return 0
    else
        return 1
    end
end
export check_for_new_closed
# checkFrNewClosed(config_paths, users) = @task check_for_new_closed(client, config, status, bot_path)

ext: jl


text:
mutable struct Irrep
    name
    characters
end

function Cn_irr(n)
    irreps = [Irrep("A",ones(n))]
    if n % 2 == 0
        b = ones(n)
        for i = 1:n
            if i % 2 == 0
                b[i] *= -1
            end
        end
        push!(irreps, Irrep("B", b))
    end
    if 2 < n < 5
        θ = 2*pi / n
        v = zeros(n)
        for j = 0:n-1
            v[j+1] += 2*cos(j*θ)
        end
        push!(irreps, Irrep("E",v))

    elseif n > 2
        θ = 2*pi / n
        l = round(Int, (n - length(irreps)) / 2)
        for i = 1:l
            v = zeros(n)
            for j = 0:n-1
                v[j+1] += 2*cos(i*j*θ)
            end
            push!(irreps, Irrep("E$i",v))
        end
    end
    return irreps
end

function Cn_irrmat(n)
    names = ["A"]
    classes = collect(1:n)
    chars = ones(n)'
    if n % 2 == 0
        push!(names,"B")
        bi = ones(n)
        for i=1:n
            if i % 2 == 0
                bi[i] *= -1
            end
        end
        chars = vcat(chars, bi')
    end
    if 2 < n < 5
        push!(names, "E")
        θ = 2*pi / n
        v = zeros(n)
        for j = 0:n-1
            v[j+1] += 2*cos(j*θ)
        end
        chars = vcat(chars, v')
    elseif n ≥ 5
        θ = 2*π / n
        l = round(Int, ((n-length(names))/2))
        for i = 1:l
            push!(names, "E$i")
            v = zeros(n)
            for j = 0:n-1
                v[j+1] += 2*cos(i*j*θ)
            end
            chars = vcat(chars, v')
        end
    end
    return names, classes, chars
end

function Cnv_irr(n)
    names, classes, chars = Cn_irrmat(n)
    names[1] = "A1"
    insert!(names, 2, "A2")
    chars = vcat(chars[1,:]', chars[1,:]', chars[2:end,:])
    for i = 2:n
        if i == n-i+2 || i > n-i+2
            break
        end
        deleteat!(classes, n-i+2)
        chars = chars[:,1:end.!=n-i+2]
    end
    if n % 2 == 0
        nirr = round(Int,(n/2)+3)
        names[3] = "B1"
        insert!(names, 4, "B2")
        chars = vcat(chars[1:3,:], chars[3,:]', chars[4:end,:])
        σv = zeros(nirr)
        σd = zeros(nirr)
        σv[1:4] = [1,-1,1,-1]
        σd[1:4] = [1,-1,-1,1]
        chars = hcat(chars, σv, σd)
    else
        nirr = round(Int, (n-1)/2+2)
        σv = zeros(nirr)
        σv[1:2] = [1, -1]
        chars = hcat(chars, σv)
    end
    return names, classes, chars
end

function Cnh_irr(n)
    names, classes, cnchars = Cn_irrmat(n)
    classes = collect(1:2*n)
    if n % 2 == 0
        newnames = []
        for i = 1:length(names)
            push!(newnames, names[i]*"u")
            names[i] = names[i]*"g"
        end
        append!(names, newnames)
        cncharsi = -1 * cnchars
        top = hcat(cnchars, cnchars)
        bot = hcat(cnchars, cncharsi)
        chars = vcat(top, bot)
    else
        newnames = []
        for i = 1:length(names)
            push!(newnames, names[i]*"″")
            names[i] = names[i]*"′"
        end
        append!(names, newnames)
        cncharsi = -1 * cnchars
        top = hcat(cnchars, cnchars)
        bot = hcat(cnchars, cncharsi)
        chars = vcat(top, bot)
    end
    return names, classes, chars
end

function Sn_irr(n)
    if n % 4 == 0
        names, classes, chars = Cn_irrmat(n)
    elseif n % 2 == 0
        ni = round(Int, n/2)
        names, classes, cnchars = Cn_irrmat(ni)
        classes = collect(1:2*n)
        newnames = []
        for i = 1:length(names)
            push!(newnames, names[i]*"u")
            names[i] = names[i]*"g"
        end
        append!(names, newnames)
        cncharsi = -1 * cnchars
        top = hcat(cnchars, cnchars)
        bot = hcat(cnchars, cncharsi)
        chars = vcat(top, bot)
    else
        throw(ArgumentError("Odd number n for S group"))
    end
    return names, classes, chars
end

function Dn_irr(n)
    if n == 2
        names = ["A", "B1", "B2", "B3"]
        classes = collect(1:4)
        chars = [1 1 1 1; 1 1 -1 -1; 1 -1 1 -1; 1 -1 -1 1]
        return names, classes, chars
    end
    names, classes, chars = Cn_irrmat(n)
    names[1] = "A1"
    insert!(names, 2, "A2")
    chars = vcat(chars[1,:]', chars[1,:]', chars[2:end,:])
    for i = 2:n
        if i == n-i+2 || i > n-i+2
            break
        end
        deleteat!(classes, n-i+2)
        chars = chars[:,1:end.!=n-i+2]
    end
    if n % 2 == 0
        nirr = round(Int,(n/2)+3)
        names[3] = "B1"
        insert!(names, 4, "B2")
        chars = vcat(chars[1:3,:], chars[3,:]', chars[4:end,:])
        C2p = zeros(nirr)
        C2pp = zeros(nirr)
        C2p[1:4] = [1,-1,1,-1]
        C2pp[1:4] = [1,-1,-1,1]
        chars = hcat(chars, C2p, C2pp)
    else
        nirr = round(Int, (n-1)/2+2)
        C2p = zeros(nirr)
        C2p[1:2] = [1, -1]
        chars = hcat(chars, C2p)
    end
    return names, classes, chars
end

function Dnh_irr(n)
    names, classes, dnchars = Dn_irr(n)
    classes = collect(1:2*n)
    if n % 2 == 0
        newnames = []
        for i = 1:length(names)
            push!(newnames, names[i]*"u")
            names[i] = names[i]*"g"
        end
        append!(names, newnames)
        dncharsi = -1 * dnchars
        top = hcat(dnchars, dnchars)
        bot = hcat(dnchars, dncharsi)
        chars = vcat(top, bot)
    else
        newnames = []
        for i = 1:length(names)
            push!(newnames, names[i]*"″")
            names[i] = names[i]*"′"
        end
        append!(names, newnames)
        dncharsi = -1 * dnchars
        top = hcat(dnchars, dnchars)
        bot = hcat(dnchars, dncharsi)
        chars = vcat(top, bot)
    end
    return names, classes, chars
end

function Dnd_irr(n)
    if n % 2 == 0
        n2 = 2*n
        names, classes, chars = Sn_irr(n2)
        classes = collect(1:2*n2)
        names[1] = "A1"
        insert!(names, 2, "A2")
        chars = vcat(chars[1,:]', chars[1,:]', chars[2:end,:])
        for i = 2:n2
            if i == n2-i+2 || i > n2-i+2
                break
            end
            deleteat!(classes, n2-i+2)
            chars = chars[:,1:end.!=n2-i+2]
        end
        nirr = n+3
        names[3] = "B1"
        insert!(names, 4, "B2")
        chars = vcat(chars[1:3,:], chars[3,:]', chars[4:end,:])
        C2p = zeros(nirr)
        σd = zeros(nirr)
        C2p[1:4] = [1,-1,1,-1]
        σd[1:4] = [1,-1,-1,1]
        chars = hcat(chars, C2p, σd)
    else
        names, classes, dnchars = Dn_irr(n)
        classes = collect(1:2*n)
        newnames = []
        for i = 1:length(names)
            push!(newnames, names[i]*"u")
            names[i] = names[i]*"g"
        end
        append!(names, newnames)
        dncharsi = -1 * dnchars
        top = hcat(dnchars, dnchars)
        bot = hcat(dnchars, dncharsi)
        chars = vcat(top, bot)
    end
    return names, classes, chars
end


text:
# Model
include_model("box")

# Horizon
T = 26

# Time step
tf = 2.5
h = tf / (T-1)

# Bounds

_uu = Inf * ones(model.m)
_uu[model.idx_u] .= 0.0 # set control to zero
_ul = zeros(model.m)
ul, uu = control_bounds(model, T, _ul, _uu)

# Initial and final states
mrp = MRP(UnitQuaternion(RotY(pi / 10.0) * RotX(pi / 15.0)))

q1 = [0.0, 0.0, 2.5, mrp.x, mrp.y, mrp.z]
v1 = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
v2 = v1 - G_func(model,q1) * h
q2 = q1 + 0.5 * h * (v1 + v2)

x1 = [q1; q2]

xl, xu = state_bounds(model, T, x1 = x1)

# Objective
obj = PenaltyObjective(1.0e5, model.m)

# Constraints
con_contact = contact_constraints(model, T)

# Problem
prob = trajectory_optimization_problem(model,
               obj,
               T,
               h = h,
               xl = xl,
               xu = xu,
               ul = ul,
               uu = uu,
               con = con_contact)

# Trajectory initialization
x0 = [0.01 * rand(model.n) for t = 1:T] #linear_interpolation(x1, x1, T) # linear interpolation on state
u0 = [0.001 * rand(model.m) for t = 1:T-1] # random controls

# Pack trajectories into vector
z0 = pack(x0, u0, prob)

#NOTE: may need to run examples multiple times to get good trajectories
# Solve nominal problem

@time z̄ , info = solve(prob, copy(z0), tol = 1.0e-5, c_tol = 1.0e-5, mapl = 5)

check_slack(z̄, prob)
x̄, ū = unpack(z̄, prob)

include(joinpath(pwd(), "models/visualize.jl"))
vis = Visualizer()
open(vis)
visualize!(vis, model, state_to_configuration(x̄), Δt = h)
settransform!(vis["/Cameras/default"], compose(Translation(1.0, 1.0, 3.0), LinearMap(RotZ(0.0))))


text:
using EngEconomics, Roots, Plots

# Given
aString = "Buy New Forklift"
aInit = -100000
aAnnual = 0
aLife = 5
aScrap = 20000

bString = "Rent Forklift"
bInit = 0
bAnnual = -20000
bLife = 5

cString = "Buy Used Forklift"
cInit = -50000
cSalvage = 0
cLife = 5
cAnnual = -5000
cAnnualG = -1000

allRevenue = 30000

MARR = 0.08

# Find
# - BC
IRR_A_eq(x) = aInit + allRevenue * seriesPresentAmountFactor(x, aLife) + aScrap * presentWorthFactor(x, aLife)
find_zero(IRR_A_eq, MARR)

IRR_B_eq(x) = (bAnnual + allRevenue) * seriesPresentAmountFactor(x, bLife)
find_zero(IRR_B_eq, MARR)

IRR_C_eq(x) = cInit + (allRevenue + cAnnual + cAnnualG * sinkingFundPaymentArithmetic(x, cLife)) * seriesPresentAmountFactor(x, cLife)
find_zero(IRR_C_eq, MARR)

# CB
ΔCB_Init = cInit - bInit
ΔCB_Scrap = cScrap - bScrap
ΔCB_Annual_eq(x) = cAnnual + cAnnualG * sinkingFundPaymentArithmetic(x, cLife) - bAnnual
ΔCB_IRR_eq(x) = ΔCB_Init + ΔCB_Scrap * presentWorthFactor(x, aLife) + ΔCB_Annual_eq(x) * seriesPresentAmountFactor(x, aLife)
xVec = collect(0:0.01:1)
yVec = zeros(size(xVec))
plot(xVec, ΔCB_IRR_eq.(xVec))
plot!(xVec, yVec)

find_zero(ΔCB_IRR_eq, MARR)

# AC
ΔAC_Init = aInit - cInit
ΔAC_Scrap = aScrap - cScrap
ΔAC_Annual_eq(x) = aAnnual - (cAnnual + cAnnualG * sinkingFundPaymentArithmetic(x, cLife))
ΔAC_IRR_eq(x) = ΔAC_Init + ΔAC_Scrap * presentWorthFactor(x, aLife) + ΔAC_Annual_eq(x) * seriesPresentAmountFactor(x, aLife)

plot(xVec, ΔAC_IRR_eq.(xVec))
plot!(xVec, yVec)

find_zero(ΔCB_IRR_eq, MARR)


text:
module QuasiHamiltonianRTN
using Expokit # expmv
using ExpmV
using SparseArrays
using LinearAlgebra

include("SpecialUnitary.jl")
include("Utilities.jl")

include("Lattice1D.jl")
include("StarGraph.jl")

export quasiHamiltonian, evolution, bloch_vector, density_operator, localized_state

"""
    Vnoise(n, γ)

Returns the generator ``V`` for the transition probabilities
between the various noise configurations for `n` fluctuators,
with switching rate ``\\gamma``.

It is defined as

``\\mathbf{V} = \\sum_{i=1}^N \\mathbf{V}_\\gamma^{(i)}``

with ``\\mathbf{V}_\\gamma^{(i)} = \\mathbf{I}_2^{\\otimes i-1} \\otimes \\mathbf{V}_\\gamma \\otimes \\mathbf{I}_2^{\\otimes N-i}``.
"""
function Vnoise(n, γ)
    V = sparse([-γ γ; γ -γ])

    function Vi(i)
        foldl(kron, [j == i ? V : sparse(1.0I, 2, 2) for j = 1:n])
    end

    sum([Vi(i) for i = 1:n])
end

"""
    quasiHamiltonian(Hamiltonian, γ)

Returns the quasi-Hamiltonian ``H_q`` for the Hamiltonian specified by the
array of matrices `H`, representing all the possible configurations of the
Hamiltonian, with switching rate γ.

The quasi-Hamiltonian is defined as in Eq. (20) of Joynt et al.
This function returns ``-i H_q``, in order to deal with real matrices instead of
imaginary ones.

## Arguments

* `Hamiltonian` is an array of matrices, each of them representing the Hamiltonian
with a different configuration of the noise.
* `γ` is the switching rate of the noise.
"""
function quasiHamiltonian(Hamiltonian, γ)::SparseMatrixCSC{Float64,Int64}
    n = size(Hamiltonian[1], 1)
    Nc = length(Hamiltonian)
    RTN_number = Int(log2(Nc))
    λ::Array{SparseArrays.SparseMatrixCSC{Complex{Float64},Int64},1} = sun_generators(n)
    f::Array{SparseArrays.SparseMatrixCSC{Float64,Int64},3} = structure_constants(λ)

    res::Array{SparseArrays.SparseMatrixCSC{Float64,Int64},1} = map(h -> igen(f,λ,h), Hamiltonian)
    Hq = blockdiag(res...) + kron(-Vnoise(RTN_number,γ), sparse(1.0I, n^2-1, n^2-1))

    return Hq
end

"""
    evolution(Hq, r0, t; [expmpkg=:ExpmV])

Evaluates the state of the system with quasi-Hamiltonian matrix `Hq`, starting
from an initial Bloch vector `r0`, at time `t`.

If `t` is a `StepRangeLen` object of length ```N_t``, return a ``N_r \\times N_t``
array containing the Bloch vector of the system at each time step as columns.
"""
function evolution(Hq::SparseMatrixCSC, r0::Vector, t; expmpkg=:ExpmV)
    @assert expmpkg ∈ (:ExpmV, :Expokit) "expmpkg must be either ExpmV or Expokit"

    n = size(Hq, 1)
    # n = Nc * Nq
    Nq = length(r0)
    Nc = Int(n / Nq)

    rtnstate = ones(Nc)/sqrt(Nc)

    v0 = kron(rtnstate, r0)
    y = kron(rtnstate, sparse(1.0I, Nq, Nq))'

    if expmpkg == :ExpmV
        res = ExpmV.expmv(-t, Hq, v0)

    elseif expmpkg == :Expokit
        res = Array{Float64}(undef, n, length(t))
        for (i, ti) in enumerate(t)
            res[:, i] = Expokit.expmv(-ti, Hq, v0)
        end
    end

    return y * res
end

function evolution(Hamiltonian::Array{SparseMatrixCSC{Float64,Int64},1},
            r0::Vector, t; γ=1, expmpkg=:ExpmV)
    n = size(Hamiltonian[1], 1)
    Nc = length(Hamiltonian)
    Nq = length(r0)

    if n^2-1 != Nq
        throw(ArgumentError("Hamiltonian and r0 have incompatible dimensions"))
    end

    Hq = quasiHamiltonian(Hamiltonian, γ)

    return evolution(Hq, r0, t; expmpkg=expmpkg)
end

end # module


text:
"""
    h = evalTF(tf, z)

Evaluates the transfer function described by `tf` at the point(s) given in the
`z` vector. `tf` must be a TransferFunction object.
"""
function evalTF(tf, z)
    if size(z) == ()
        return tf(z)[1, 1]
    else
        return tf(z, false)[:, 1, 1]
    end
end

max_stars_repo_name: Mahmoud-Kharsa/DeltaSigma

max_issues_repo_name: Mahmoud-Kharsa/DeltaSigma

max_forks_repo_name: Mahmoud-Kharsa/DeltaSigma


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

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local iv = try Base.loaded_modules[Base.PkgId(Base.UUID("6e696c72-6542-2067-7265-42206c756150"), "AbstractPlutoDingetjes")].Bonds.initial_value catch; b -> missing; end
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : iv(el)
        el
    end
end

# ╔═╡ 9a8e1c63-f27c-441b-8eb6-aca30f4ca2d8
# ╠═╡ show_logs = false
begin
	using Logging
	global_logger(NullLogger())
	display("")
end

# ╔═╡ f0545c67-5cfd-438f-a9ef-92c35ebaefa4
# ╠═╡ show_logs = false
#Set-up packages
begin

	using DataFrames, HTTP, CSV, Dates, Plots, PlutoUI, Printf, LaTeXStrings, HypertextLiteral, ShortCodes, XLSX

	gr();
	Plots.GRBackend()

	#Define html elements
	nbsp = html"&nbsp" #non-breaking space
	vspace = html"""

"""
	br = html"
"

	#Sets the width of cells, caps the cell width by 90% of screen width
	#(setting overwritten by cell below)
	# @bind screenWidth @htl("""
	#

	#
	#

	# """)

	# cellWidth= min(1000, screenWidth*0.9)
	# @htl("""
	#
	# """)

	#Sets the width of the cells
	#begin
	#	html"""


text:
module ConstrainedDynamicsVis

using LinearAlgebra
using StaticArrays
using ConstrainedDynamics
using ConstrainedDynamics: AbstractMechanism, Shape, vrotate, getid

using CoordinateTransformations
using GeometryBasics: GeometryBasics, Point, Vec
using Blink
using MeshCat
using MeshCat: js_scaling, js_position, js_quaternion

export visualize,
    convertvideo

include("convertshape.jl")
include("visualize.jl")

include("video.jl")

end


text:
#######################################
## Contains osm2network, to transform an osm file into a network
#######################################

const ROAD_CLASSES = Dict(
    "motorway" => 1,
    "trunk" => 2,
    "primary" => 3,
    "secondary" => 4,
    "tertiary" => 5,
    "unclassified" => 6,
    "residential" => 6,
    "service" => 7,
    "motorway_link" => 1,
    "trunk_link" => 2,
    "primary_link" => 3,
    "secondary_link" => 4,
    "tertiary_link" => 5,
    "living_street" => 8,
    "pedestrian" => 8,
    "road" => 6)

#TOOLS from OpenStreetMapParser
highway(w::OSMWay) = haskey(w.tags, "highway")
roadway(w::OSMWay) = get(ROAD_CLASSES,w.tags["highway"],0.)
function oneway(w::OSMWay)
    v = get(w.tags,"oneway", "")

    if v == "false" || v == "no" || v == "0"
        return false
    elseif v == "-1" || v == "true" || v == "yes" || v == "1"
        return true
    end

    highway = get(w.tags,"highway", "")
    junction = get(w.tags,"junction", "")

    return (highway == "motorway" ||
            highway == "motorway_link" ||
            junction == "roundabout")
end

visible{T <: OSMElement}(obj::T) = (get(obj.tags, "visible", "") != "false")
services(w::OSMWay) = (get(w.tags,"highway", "") == "services")
reverse(w::OSMWay) = (get(w.tags,"oneway", "") == "-1")
toradians(degree::AbstractFloat) = degree * π / 180.0

"""
Opens a OSM file,
and returns a network object from the OSM data
"""
function osm2network(filename::AbstractString)
  osm = parseOSM(filename);

  # fetch all nodes ids
  osm_id = Int[n.id for n in osm.nodes]
  # construct inverse dictionary
  node_id = Dict{Int,Int}()
  sizehint!(node_id, length(osm_id))
  for i in 1:length(osm_id)
      node_id[osm_id[i]] = i
  end
  #construct the "Node" objects
  lonlat = Tuple{Float64,Float64}[c.lonlat for c in osm.nodes]
  bounds = boundingBox(lonlat)
  center = ((bounds[2]+bounds[1])/2, (bounds[4]+bounds[3])/2)
  nodes = Array{Node}(length(lonlat))
  for (i,(lon,lat)) in enumerate(lonlat)
      x,y = toENU(lon,lat,center)
      nodes[i] = Node(x,y,lon,lat)
  end

  g = DiGraph(length(nodes))

  #We transform ways into edges
  roads = Dict{Tuple{Int,Int},Road}()
  for w in osm.ways
    if highway(w) && visible(w) && !(services(w)) && roadway(w) != 0
      for i in 2:length(w.nodes)
        if reverse(w)
          edge = Edge(node_id[w.nodes[i]],node_id[w.nodes[i-1]])
        else
          edge = Edge(node_id[w.nodes[i]],node_id[w.nodes[i-1]])
        end
        if !has_edge(g,edge)
          add_edge!(g,edge)
          roads[src(edge),dst(edge)] = Road(src(edge),dst(edge),distanceGeo(nodes[src(edge)],nodes[dst(edge)]), roadway(w))
        end
        if !oneway(w)
          edge = Edge(dst(edge),src(edge))
          if !has_edge(g,edge)
            add_edge!(g,edge)
            roads[src(edge),dst(edge)] = Road(src(edge),dst(edge),distanceGeo(nodes[src(edge)],nodes[dst(edge)]), roadway(w))
          end
        end
      end
    end
  end
  n = Network(g, nodes, roads)
  n.projcenter = center
  return n
end

max_stars_repo_name: andrew-zhm/RoutingNetworks

max_issues_repo_name: andrew-zhm/RoutingNetworks

max_forks_repo_name: andrew-zhm/RoutingNetworks


text:
# ---
# title: 1331. Rank Transform of an Array
# id: problem1331
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: Array
# link: <https://leetcode.com/problems/rank-transform-of-an-array/description/>
# hidden: true
# ---
#
# Given an array of integers `arr`, replace each element with its rank.
#
# The rank represents how large the element is. The rank has the following
# rules:
#
#   * Rank is an integer starting from 1.
#   * The larger the element, the larger the rank. If two elements are equal, their rank must be the same.
#   * Rank should be as small as possible.
#
#
#
# **Example 1:**
#
#
#
#     Input: arr = [40,10,20,30]
#     Output: [4,1,2,3]
#     **Explanation** : 40 is the largest element. 10 is the smallest. 20 is the second smallest. 30 is the third smallest.
#
# **Example 2:**
#
#
#
#     Input: arr = [100,100,100]
#     Output: [1,1,1]
#     **Explanation** : Same elements share the same rank.
#
#
# **Example 3:**
#
#
#
#     Input: arr = [37,12,28,9,100,56,80,5,12]
#     Output: [5,3,4,2,8,6,7,1,3]
#
#
#
#
# **Constraints:**
#
#   * `0 <= arr.length <= 105`
#   * `-109 <= arr[i] <= 109`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
# Note that since Base Julia already has a `Symbol` object, we use the term
# `SymtabEntry` to refer to the concrete instantiation of a symbol within an
# object file, and the term `SymbolRef` as a reference to that object.

# Export Symbols API
export Symbols,
       getindex, length, iterate, lastindex, eltype,
       handle, header

# Export SymtabEntry API
export SymtabEntry,
       deref, symbol_name, symbol_value, isundef, isglobal, islocal, isweak

# Export SymbolRef API
export SymbolRef,
       symbol_number

# Import iteration protocol
import Base: length, iterate, lastindex

"""
    Symbols

An abstraction over the concept of a collection of symbol (`SymtabEntry`) types
within an object file.  One can think of the `Symbols` object containing the
table of symbols within the object file, whereas the `SymtabEntry`/`SymbolRef`
objects contain the actual symbol data itself.  The list of available API
operations is given below, with methods that subclasses must implement marked
in emphasis:

### Creation
  - *Symbols()*

### Iteration
  - getindex()
  - *lastindex()*
  - length()
  - iterate()
  - eltype()

### Misc.
  - *handle()*
"""
abstract type Symbols{H<:ObjectHandle} end

# Fairly simple iteration interface specification
@mustimplement lastindex(syms::Symbols)
length(syms::Symbols) = lastindex(syms)
iterate(syms::Symbols, idx=1) = idx > length(syms) ? nothing : (syms[idx], idx+1)
eltype(::Type{S}) where {S <: Symbols} = SymbolRef

function getindex(syms::Symbols{H}, idx) where {H <: ObjectHandle}
    # Punt off to `getindex_ref`
    oh = handle(syms)
    return getindex_ref(
        syms,
        symtab_entry_offset(oh),
        symtab_entry_size(oh),
        symtab_entry_type(oh),
        SymbolRef,
        idx
    )
end

"""
    SymtabEntry

An abstraction over the concept of a symbol within an object file.  This type
does not use the `Symbol` name as this would conflict with the builtin Julia
`Symbol` type, so the name `SymtabEntry` is used instead.  As a user, the
`SymbolRef` type should be the primary method of interacting with symbols, as a
developer adding new object file formats, some methods must support
`SymtabEntry`s, others must support only `SymbolRef`s.  Note that any method
that works on a `SymtabEntry` must also work with a `SymbolRef`, see the
`@derefmethod` macro for a convenient helper macro to generate `SymbolRef` ->
`SymtabEntry` wrapper methods. The list of available API operations is given
below, with methods that subclasses must implement marked in emphasis:

### Creation:
  - *SymtabEntry()*

### Util:
  - deref()

### Properties:
  - *symbol_name()*
  - *symbol_value()*
  - *isundef()*
  - *isglobal()*
  - *islocal()*
  - *isweak()*
"""
abstract type SymtabEntry{H<:ObjectHandle} end

# Dummy `deref()`
deref(sym::SymtabEntry) = sym

"""
    symbol_name(sym::SymtabEntry)

Return the name of the given section as a string.  In order to return a true
name, it is necessary to perform a lookup within the object's string table,
which cannot be done using just a `SymtabEntry` object; use a `SymbolRef`
object instead if you need that.  For sanity sake, this method will return a
string, but the contents of the string may be something like the offset within
the string table pointing to this `SymtabEntry`'s name, e.g. "@strtab.123"
"""
@mustimplement symbol_name(sym::SymtabEntry)

"""
    symbol_value(sym::SymtabEntry)

Return the value of the given symbol
"""
@mustimplement symbol_value(sym::SymtabEntry)

"""
    isundef(sym::SymtabEntry)

Return `true` if the given symbol is undefined
"""
@mustimplement isundef(sym::SymtabEntry)

"""
    isglobal(sym::SymtabEntry)

Return `true` if the given symbol is global
"""
@mustimplement isglobal(sym::SymtabEntry)

"""
    islocal(sym::SymtabEntry)

Return `true` if the given symbol is local
"""
@mustimplement islocal(sym::SymtabEntry)

"""
    isweak(sym::SymtabEntry)

Return `true` if the given symbol is weak
"""
@mustimplement isweak(sym::SymtabEntry)

"""
    SymbolRef

Provides a reference to a `SymtabEntry`, along with a reference to the
`ObjectHandle` this `SymtabEntry` comes from.  This should be the primary
method by which users interact with symbols inside object files.  The list of
available API operations is given below, with methods that subclasses must
implement marked in emphasis.  Note that this overlaps heavily with the
`SymtabEntry` object API, this is by design as many of the methods are simply
passthroughs to the underlying `SymtabEntry` API calls for ease of use.

### Creation:
  - *SymbolRef()*

### Util:
  - *deref()*
  - *Symbols()*
  - handle()

### Properties:
  - *symbol_number()*
  - *symbol_name()*
  - symbol_value()
  - isundef()
"""
abstract type SymbolRef{H<:ObjectHandle} end

"""
SymbolRef(symbols::Symbols, sym::SymtabEntry, idx::UInt32)

Construct a `SymbolRef` object pointing to the given `SymtabEntry`, which
itself represents the `idx`'th symbol within the given `Symbols` collection.
"""
@mustimplement SymbolRef(symbols::Symbols, sym::SymtabEntry, idx)
@mustimplement deref(sym::SymbolRef)

"""
    Symbols(sym::SymbolRef)

Return the `Symbols` object that this `SymbolRef` belongs to.
"""
@mustimplement Symbols(sym::SymbolRef)

"""
    handle(sym::SymbolRef)

Return the `ObjectHandle` that this `SymbolRef` belongs to.
"""
handle(sym::SymbolRef) = handle(Symbols(sym))

"""
    symbol_number(sym::SymbolRef)

Return the number (index) of the given symbol.
"""
@mustimplement symbol_number(sym::SymbolRef)

"""
    symbol_name(sym::SymbolRef)

Return the name of the given symbol as a string.  This method often performs
some kind of lookup within the string table of the object to get the full name
of the symbol.
"""
@mustimplement symbol_name(sym::SymbolRef)

@derefmethod symbol_value(sym::SymbolRef)
@derefmethod isundef(sym::SymbolRef)
@derefmethod isglobal(sym::SymbolRef)
@derefmethod islocal(sym::SymbolRef)
@derefmethod isweak(sym::SymbolRef)

## Printing
function show(io::IO, sym::Union{SymtabEntry{H},SymbolRef{H}}) where {H <: ObjectHandle}
    print(io, "$(format_string(H)) Symbol")

    if !get(io, :compact, false)
        println(io)
        println(io, "       Name: $(symbol_name(sym))")
        println(io, "      Value: $(symbol_value(sym))")
        println(io, "    Defined: $(isundef(sym) ? "No" : "Yes")")
        println(io, "     Strong: $(isweak(sym) ? "No" : "Yes")")
        print(io,   "   Locality: $(isglobal(sym) ? "Global" : "Local")")
    else
        print(io, " \"$(symbol_name(sym))\"")
    end
end

show(io::IO, syms::Symbols{H}) where {H <: ObjectHandle} = show_collection(io, syms, H)


text:
# A spec file generator for AWS APIs.
#
# Since AWS does not publish WSDL files anymore,
# this crawls AWS API specification web pages and
# extracts relevant information as a Julia specification file.
#
# Parsing (of copmplex types particularly) is not very accurate.
# Since scraping web pages is inherently flaky, it seems better
# to keep the logic simple and easy to change. Scraping logic
# may need tweaks if changes to the HTML pages break it. The
# generated file needs manual verification and minor edits.

using Cascadia
using Gumbo
using Requests

function attr_type(elem)
    txt = nodeText(elem)
    is_array = contains(txt, "Array") || contains(txt, "array")
    parts = split(txt)
    elem_type = is_array ? parts[4] : parts[2]
    if lowercase(elem_type) == "boolean"
        elem_type = "Bool"
    elseif lowercase(elem_type) in ("string", "strings")
        elem_type = "String"
    elseif lowercase(elem_type) in ("integer", "long", "integers", "longs")
        elem_type = "Int"
    elseif lowercase(elem_type) in ("float", "floats", "double", "doubles")
        elem_type = "Float64"
    elseif lowercase(elem_type) == "base64-encoded"
        elem_type = "Vector{UInt8}"
    elseif lowercase(elem_type) == "timestamp"
        elem_type = "Base.Dates.DateTime"
    else
        elem_type = String(elem_type)
    end
    is_array ? "Vector{$(elem_type)}" : elem_type
end

function attr_name(elem)
    elem_name = nodeText(elem)
    if endswith(elem_name, ".N")
        elem_name = rsplit(elem_name, '.'; limit=2)[1]
    elseif contains(elem_name, "(request)") && contains(elem_name, "(response)") && contains(elem_name, ",")
        for name in split(elem_name, ','; limit=2)
            if contains(name, "(request)")
                elem_name = split(name)[1]
                break
            end
        end
    end
    String(elem_name)
end

function parse_type_attrs(elem)
    qs = matchall(Selector(".term"), elem)
    attr_names = map(attr_name, qs)

    qs = matchall(Selector("p:contains(\"Type:\")"), elem)
    attr_types = map(attr_type, qs)

    attributes = Vector{Pair{String,String}}()
    for (a,b) in zip(attr_names, attr_types)
        push!(attributes, a=>b)
    end

    attributes
end

function parse_type_name(elem)
    qs = matchall(Selector(".topictitle"), elem)
    nodeText(qs[1])
end

function parse_type_id(elem)
    qs = matchall(Selector(".topictitle"), elem)
    qs[1].attributes["id"]
end

function retry_get(url)
    while true
        try
            return get(url)
        catch ex
            warn("will retry ", url, " - error: ", ex)
            sleep(5)
        end
    end
end

"""
Generate type specification for an individual type.
"""
function gen_type(url)
    r = retry_get(url)
    h = parsehtml(bytestring(r.data))

    typename = parse_type_name(h.root)
    attributes = parse_type_attrs(h.root)

    info(typename, " => ", attributes)
    typename => attributes
end

"""
Parse the type collection page and
iterate through individual type pages.
"""
function gen_types(url)
    all_types = Vector{Pair{String,Vector{Pair{String,String}}}}()

    info("type list: ", url)
    r = retry_get(url)
    h = parsehtml(bytestring(r.data))

    qs = matchall(Selector(".listitem > p > a"), h.root)
    baseurl = rsplit(url, '/'; limit=2)[1]
    info("got ", length(qs), " types")
    for elem in qs
        typeuri = elem.attributes["href"]
        typeurl = baseurl * "/" * typeuri
        info("type spec ", length(all_types)+1, " of ", length(qs), ": ", typeurl)
        push!(all_types, gen_type(typeurl))
    end

    all_types
end

"""
Generate API specification for an individual operation.
"""
function gen_api(url)
    r = retry_get(url)
    h = parsehtml(bytestring(r.data))

    apiname = parse_type_name(h.root)
    apiid = parse_type_id(h.root)

    q_reqp = matchall(Selector("#"*apiid*"_RequestParameters"), h.root)
    q_resp = matchall(Selector("#"*apiid*"_ResponseElements"), h.root)
    params_list = String[]
    isempty(q_reqp) || push!(params_list, "Request")
    isempty(q_resp) || push!(params_list, "Result")
    var_list = matchall(Selector(".variablelist"), h.root)

    req_resp = Dict()
    for idx in 1:length(params_list)
        param_type = params_list[idx]
        vars = var_list[idx]
        attributes = parse_type_attrs(vars)
        req_resp[param_type] = attributes
    end

    info(apiname, " => ", req_resp)
    apiname => req_resp
end

"""
Parse the operations collection page and
iterate through individual operation pages.
"""
function gen_apis(url)
    all_apis = Dict()

    info("api list: ", url)
    r = retry_get(url)
    h = parsehtml(bytestring(r.data))

    qs = matchall(Selector(".listitem > p > a"), h.root)
    baseurl = rsplit(url, '/'; limit=2)[1]
    info("got ", length(qs), " apis")
    for elem in qs
        apiuri = elem.attributes["href"]
        apiurl = baseurl * "/" * apiuri
        info("api spec ", length(all_apis)+1, " of ", length(qs), ": ", apiurl)
        push!(all_apis, gen_api(apiurl))
    end

    all_apis
end

"""
Generate API specification for `apiname` from the provided
URLs pointing to types and operations specifications.
Also create a Julia specification file with appropriate name
that can be used subsequently with `gen.jl`.
"""
function write_spec(apiname, typesurl, apisurl)
    all_types = gen_types(typesurl)
    all_apis = gen_apis(apisurl)

    specfile = lowercase(apiname)*"_spec.jl"
    open(specfile, "w") do outfile
        println(outfile, "const ", apiname, "Types = [")
        isfirst = true
        for (type_name,type_spec) in all_types
            println(outfile, isfirst ? "    :" : "    ,:", type_name, " => [")
            isfirstattr = true
            for (attr_name,attr_type) in type_spec
                println(outfile, isfirstattr ? "        :" : "        ,:", attr_name, " => ", attr_type)
                isfirstattr = false
            end
            println(outfile, "    ]")
            isfirst = false
        end
        println(outfile, "]")
        println(outfile, "")
        println(outfile, "const ", apiname, "Api = Dict(")
            isfirst = true
            for (api_name,api_spec) in all_apis
                println(outfile, isfirst ? "    :" : "    ,:", api_name, " => Dict(")
                isfirstparam = true
                for (param_name,param_spec) in api_spec
                    println(outfile, isfirstparam ? "        :" : "        ,:", param_name, " => [")
                    isfirstattr = true
                    for (attr_name,attr_type) in param_spec
                        println(outfile, isfirstattr ? "            :" : "            ,:", attr_name, " => ", attr_type)
                        isfirstattr = false
                    end
                    println(outfile, "        ]")
                    isfirstparam = false
                end
                println(outfile, "    )")
                isfirst = false
            end
        println(outfile, ")")
    end

    typefile = lowercase(apiname)*"_gen_types.jl"
    open(typefile, "w") do outfile
        for (type_name,type_spec) in all_types
            println(outfile, "type $type_name end")
        end
    end
end

# generate and write out specification for each API class
write_spec("AutoScaling", "http://docs.aws.amazon.com/AutoScaling/latest/APIReference/API_Types.html", "http://docs.aws.amazon.com/AutoScaling/latest/APIReference/API_Operations.html")
write_spec("SQS", "http://docs.aws.amazon.com/AWSSimpleQueueService/latest/APIReference/API_Types.html", "http://docs.aws.amazon.com/AWSSimpleQueueService/latest/APIReference/API_Operations.html")
write_spec("EC2", "http://docs.aws.amazon.com/AWSEC2/latest/APIReference/API_Types.html", "http://docs.aws.amazon.com/AWSEC2/latest/APIReference/API_Operations.html")


text:
using FiniteVolumes, StaticArrays, OrdinaryDiffEq, Plots

mesh = PeriodicCartesianMesh(50, 50)
u0 = [SVector(1.0 + exp(-100*((x - 0.4)^2 + (y - 0.3)^2)), 0.0, 0.0) for (x, y) in cell_centers(mesh)]

function f(v, n)
    h, ux, uy = v[1], v[2]/v[1], v[3]/v[1]
    un = ux*n[1] + uy*n[2]
    return SVector(h*un, h*ux*un + h^2*n[1]*9.81/2, h*uy*un + h^2*n[2]*9.81/2)
end

dudt(u, t, p) = -FiniteVolumes.div(f, mesh, u, Upwind())

problem = ODEProblem(dudt, u0, (0.0, 0.2))
sol = solve(problem, Euler(), dt=0.001, saveat=0.01)

plot(
    plot(mesh, sol(0.0), 1),
    plot(mesh, sol(0.06), 1),
    plot(mesh, sol(0.13), 1),
    plot(mesh, sol(0.20), 1),
)


text:
function sphere_vertex(x, y, z, scale)
    length = sqrt(x^2 + y^2 + z^2)
    return [(i * scale) / length for i in (x, y, z)]
end
function sphere_vertex(a::Vector{T}, scale) where T
    @assert length(a) == 3
    x, y, z = a[1], a[2], a[3]
    return sphere_vertex(x, y, z, scale)
end

function middle_point!(point_1, point_2, middle_point_cache, verts, scale)
    # We check if we have already cut this edge first
    # to avoid duplicated verts
    smaller_index = min(point_1, point_2)
    greater_index = max(point_1, point_2)
    key = "$smaller_index, $greater_index"
    if haskey(middle_point_cache, key)
        return middle_point_cache[key]
    end
    vert_1 = verts[point_1]
    vert_2 = verts[point_2]
    middle = [sum(i) / 2 for i in zip(vert_1, vert_2)]
    append!(verts, [sphere_vertex(middle, scale)])
    index = length(verts)
    middle_point_cache[key] = index
    return index
end # function

"""
Create icosphere
"""
function icosphere(subdiv = 3, scale = 1.0)
    PHI = (1. + sqrt(5.)) / 2.
    verts = [(sphere_vertex(-1., PHI, 0., scale)),
        (sphere_vertex(1., PHI, 0., scale)),
        (sphere_vertex(-1., -PHI, 0., scale)),
        (sphere_vertex(1., -PHI, 0., scale)),
        (sphere_vertex(0., -1., PHI, scale)),
        (sphere_vertex(0., 1., PHI, scale)),
        (sphere_vertex(0., -1., -PHI, scale)),
        (sphere_vertex(0., 1., -PHI, scale)),
        (sphere_vertex(PHI, 0., -1., scale)),
        (sphere_vertex(PHI, 0., 1., scale)),
        (sphere_vertex(-PHI, 0., -1., scale)),
        (sphere_vertex(-PHI, 0., 1., scale))]
    faces = [[1, 12, 6], # 5 faces around point 0
            [1, 6, 2],
            [1, 2, 8],
            [1, 8, 11],
            [1, 11, 12],
            [2, 6, 10], # Adjacent faces
            [6, 12, 5],
            [12, 11, 3],
            [11, 8, 7],
            [8, 2, 9],
            [4, 10, 5],  # 5 faces around 3
            [4, 5, 3],
            [4, 3, 7],
            [4, 7, 9],
            [4, 9, 10],
            [5, 10, 6],# Adjacent faces
            [3, 5, 12],
            [7, 3, 11],
            [9, 7, 8],
            [10, 9, 2]]

    middle_point_cache = Dict{String,Int}()
    for i in 1:subdiv
        faces_subdiv = Array{Int64,1}[]
        for tri in faces
            v1 = middle_point!(tri[1], tri[2], middle_point_cache, verts, scale)
            v2 = middle_point!(tri[2], tri[3], middle_point_cache, verts, scale)
            v3 = middle_point!(tri[3], tri[1], middle_point_cache, verts, scale)
            append!(faces_subdiv, [[tri[1], v1, v3]])
            append!(faces_subdiv, [[tri[2], v2, v1]])
            append!(faces_subdiv, [[tri[3], v3, v2]])
            append!(faces_subdiv, [[v1, v2, v3]])
            faces = faces_subdiv
        end
    end
    v = Vector{Point3f0}(verts)
    f = Vector{TriangleFace{Int}}(faces)
    sphere = normal_mesh(Mesh(v, f))
    return sphere
end


text:
##
using AeroMDAO

##
foilpath = joinpath((dirname ∘ dirname ∘ pathof)(AeroMDAO), "test/CRM.dat")

## Foil processing
coords = read_foil(foilpath)

## Cosine spacing
cos_foil = cosine_foil(coords, 51)

## Camber-thickness transformations
xcamthick = foil_camthick(cos_foil)
foiler = camthick_foil(xcamthick[:,1], xcamthick[:,2], xcamthick[:,3])

## Fitting to CST
num_dv = 4

# Coordinates fitting
up, low          = split_foil(cos_foil)
alpha_u, alpha_l = coords_to_CST(up, num_dv), coords_to_CST(low, num_dv)
cst_foil         = kulfan_CST(alpha_u, alpha_l, (1e-4, -1e-4), 0.0)

## Camber-thickness fitting
alphas   = camthick_to_CST(cos_foil, num_dv)
cam_foil = camber_CST(alphas..., (0., 2e-4), 0)

## Kulfan CST
alpha_u = [0.1, 0.3, 0.2, 0.15, 0.2]
alpha_l = [-0.1, -0.1, -0.1, -0.001, -0.02]
dzs     = (1e-4, 1e-4)
foil    = kulfan_CST(alpha_u, alpha_l, dzs, 0.2)

## NACA 4-digit airfoils
naca = naca4((2,4,1,2))

## Plotting library
using Plots
plotlyjs();

## Cosine and camber-thickness
plot(getindex.(cos_foil, 1), getindex.(cos_foil, 2), label = "Cosine")
plot!(xcamthick[:,1], xcamthick[:,2], label = "Camber")
plot!(xcamthick[:,1], xcamthick[:,3], label = "Thickness", aspectratio = 1)
plot!(getindex.(cst_foil, 1), getindex.(cst_foil, 2), label = "CST Coordinates Fit")
plot!(getindex.(cam_foil, 1), getindex.(cam_foil, 2), label = "CST Camber-Thickness Fit", aspectratio = 1)
plot!(getindex.(foil, 1), getindex.(foil, 2), label = "CST", aspectratio = 1)
plot!(getindex.(naca, 1), getindex.(naca, 2), label = "NACA", aspectratio = 1)

max_forks_repo_name: GodotMisogi/AeroMDAO


text:
using Juriba.Particles
const PT = Particles
using QuadGK

@Base.kwdef mutable struct Powerlaw <: AbstractParticles
    pspec::Float64 = 0
    plnorm::Float64 = 0
    pmin::Float64 = 0
    pmax::Float64 = 0

    size::Int
    mass_gr::Float64   # particle mass in grams
    mass_kev::Float64  # same as above but in keV, using electrons as "reference"
    p::Vector{Float64} # array of particle momenta
    ndens::Vector{Float64}			# array of number density per unit volume, per unit momentum
    gamma::Vector{Float64}			# array of particle kinetic energies for each momentum
    gdens::Vector{Float64}		    # array of number density per unit volume, per unit gamma
    gdens_diff::Vector{Float64}     # array with differential of number density for radiation calculation

end

function powerlaw(s)
    size = s

    plnorm = 1.0

    mass_gr = emgm
    mass_kev = emgm*gr_to_kev

    p = zeros(Float64, size)
    ndens = similar(p)
    gamma = similar(p)
    gdens = similar(p)
    gdens_diff = similar(p)

    return Powerlaw(plnorm=plnorm, size = size, mass_gr = mass_gr, mass_kev = mass_kev, p=p, ndens=ndens,
                    gamma=gamma, gdens=gdens, gdens_diff=gdens_diff)
end

function set_p(pl::Powerlaw, min, ucom, bfield, betaeff, r, fsc)

    pl.pmin = min
    pl.pmax = max_p(pl, ucom, bfield, betaeff, r, fsc)

    pinc = (log10(pmax) - log10(pmin))/(pl.size - 1)

    for i in eachindex(pl.p)
        pl.p[i] = 10^(log10(pl.pmin)+i*pinc)
        pl.gamma[i] = sqrt((pl.p[i]/(pl.mass_gr*cee))^2 + 1)
    end

end

function set_p(pl::Powerlaw, min, gmax)
    pl.pmin = min
    pl.pmax = sqrt(gmax^2 - 1)*pl.mass_gr*cee

    pinc = (log10(pmax) - log10(pmin))/(size - 1)

    for i in eachindex(pl.p)
        pl.p[i] = 10^(log10(pl.pmin) + i*pinc)
        pl.gamma[i] = sqrt((pl.p[i]/(pl.mass_gr*cee))^2 + 1)
    end

end

function set_ndens(pl::Powerlaw)

    for i = 1:pl.size
        pl.ndens[i] = pl.plnorm*pl.p[i]^(-pl.psec)*exp(-pl.p[i]/pl.pmax)
    end

    PT.initialize_gdens(pl)
    PT.gdens_differentiate(pl)
end

function set_pspec(pl::Powerlaw, s1)
    pl.spec = s1
end

function set_norm(pl::Powerlaw, n)
    pl.plnorm = n*(1.0 - pl.pspec)/(pl.p[end]^(1 - pl.spec) - pl.p[1]^(1.0 - pl.pspec))
end

""" Injection function to be integrated in cooling """
function injection_pl_int(x, p::PT.injection_pl_params)
    params = p

    s = params.s
    n = params.n
    m = params.m
    max = params.max

    mom_int = sqrt((x^2) - 1)*m*cee

    return n*mom_int^(-s)*exp(-mom_int/max)
end

function cooling_steadystate(pl::Powerlaw, ucom, n0, bfield, r, betaeff)

    Urad = bfield^2/8π + ucom
    pdot_ad = betaeff*cee/r
    pdot_rad = (4.0*sigtom*cee*Urad)/(3.0*pl.mass_gr*cee^2)
    tinj = r/cee

    params = PT.injection_pl_params(pl.psec, pl.plnorm, pl.mass_gr, pl.pmax)

    func(x) = PT.injection_pl_int(x, params)

    for i = 1:pl.size-1
        integral, error = quadgk(func, pl.gamma[i], pl.gamma[i+1], rtol=1e1, atol=1e1, maxevals=100)

        pl.ndens[i] = (integral/tinj)/(pdot_ad*pl.p[i]/(pl.mass_gr*cee)+pl.pdot_rad*(pl.gamma[i]*pl.p[i]/(pl.mass_gr*cee)))
    end

    pl.ndens[end] = pl.ndens[end-1]*(pl.p[end]/pl-p[end-2])^(-pl.pspec - 1)*exp(-1)

    # the last bin is set by arbitrarily assuming cooled distribution; this is necessary because the integral
    # above is undefined for the last bin

    # The last step requires a renormalization. The reason is that the result of gsl_integration_qag strongly
    # depends on the value of "size". Without doing anything fancy, this can be fixed simply by ensuring that
    # the total integrated number of density equals n0 (which we know), and rescaling the array ndens[i] by
    # the appropriate constant.

    renorm = PT.count_particles(pl)/n0

    ndens[i] ./= renorm

    PT.initialize_gdens(pl)
    PT.gdens_differentiate(pl)

end

"""
Method to calculate maximum momentum of non thermal particles based on acceleration and cooling timescales
The estimate is identical to the old agnjet but in momentum space; see Lucchini et al. 2019 for the math of
the old version
"""
function max_p(pl::Powerlaw, ucom, bfield, betaeff, r, fsc)

    Urad = bfield^2/8π + ucom
    escom = betaeff*cee/r
    syncon = (4.0*sigtom*Urad)/(3.0*pl.mass_gr*cee)

    accon = (3.0*fsc*pl.charg*bfield)/(4.0*pl.mass_gr*cee)

    b = escom/syncon
    c = accon/syncon

    gmax = (-b + sqrt(b^2 + 4.0*c))/2.0
    return sqrt(gmax^2 - 1)*pl.mass_gr*cee

end

function test(pl::Powerlaw)

    println("Power-law distribution")
    println("pspec: $(pl.pspec)")
    println("Array size: $(pl.size)")
    println("Default normalization: $(pl.plnorm)")
    println("Particle mass in grams: $(pl.mass_gr)")

end


text:
using OptionalData
using Test

struct TestType
    a::Float64
    b::Float64
    c::Float64
end

@OptionalData opt1 Symbol
@OptionalData opt2 TestType "Try again."

@testset "OptionalData" begin
    @test string(opt1) == "OptData{Symbol}()"
    @test_throws NoDataError get(opt1)
    @test sprint(showerror, NoDataError(opt1)) == "`opt1` is not available."
    push!(opt1, :Test)
    @test string(opt1) == "OptData{Symbol}(Test)"
    @test get(opt1) == :Test
    push!(opt1, "Test")
    @test get(opt1) == :Test

    @test string(opt2) == "OptData{TestType}()"
    @test_throws NoDataError get(opt2)
    @test sprint(showerror, NoDataError(opt2)) == "`opt2` is not available. Try again."
    push!(opt2, TestType(1, 2, 3))
    @test string(opt2) == "OptData{TestType}(TestType(1.0, 2.0, 3.0))"
    @test get(opt2) == TestType(1, 2, 3)
    push!(opt2, 1, 2, 3)
    @test get(opt2) == TestType(1, 2, 3)
end

max_stars_repo_name: helgee/OptionalData

max_issues_repo_name: helgee/OptionalData

max_forks_repo_name: helgee/OptionalData


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

using Markdown
using InteractiveUtils

# ╔═╡ e077ab16-22e5-11eb-15ed-5f0654bf7328
using Pkg, DrWatson

# ╔═╡ e75a9a06-22e5-11eb-2e36-8d4f62b830ed
begin
	@quickactivate "StatisticsWithJuliaPlutoNotebooks"
	using Random, Distributions, DataFrames, Plots
	Random.seed!(0)
end;

# ╔═╡ ca5fbdc8-22e5-11eb-0a60-cb0a127a3ca5
md"## Listing6.14"

# ╔═╡ 1aa41898-3afa-11eb-2f49-570853ce6435
begin
	lambda = 0.1
	dist = Exponential(1/lambda)
	actualMedian = median(dist)

	M = 10^3
	N = 10^4
	nRange = 2:2:10
	alpha = 0.05
end

# ╔═╡ fba0be84-3b2f-11eb-1e0f-5f8a1f264dba
begin
	df = DataFrame()
	for n in nRange
		coverageCount = 0
		for _ in 1:M
			sampleData = rand(dist, n)
			bootstrapSampleMeans = [median(rand(sampleData, n)) for _ in 1:N]
			L = quantile(bootstrapSampleMeans, alpha/2)
			U = quantile(bootstrapSampleMeans, 1-alpha/2)
			coverageCount += L < actualMedian && actualMedian < U
		end
		append!(df, DataFrame(:n => n, :coverage => coverageCount/M))
	end
	df
end

# ╔═╡ 475ff888-22e6-11eb-2354-09f8f40a8e12
md"## End of listing6.14"

# ╔═╡ Cell order:
# ╟─ca5fbdc8-22e5-11eb-0a60-cb0a127a3ca5
# ╠═e077ab16-22e5-11eb-15ed-5f0654bf7328
# ╠═e75a9a06-22e5-11eb-2e36-8d4f62b830ed
# ╠═1aa41898-3afa-11eb-2f49-570853ce6435
# ╠═fba0be84-3b2f-11eb-1e0f-5f8a1f264dba
# ╟─475ff888-22e6-11eb-2354-09f8f40a8e12


text:
using Documenter, FLOWFarm

DocMeta.setdocmeta!(FLOWFarm, :DocTestSetup, :(using FLOWFarm); recursive=true)

makedocs(;
    modules=[FLOWFarm],
    format=Documenter.HTML(),
    pages=[
        "Intro" => "index.md",
        "Quick Start" => "Tutorial.md",
        "How-to Guide" => "How_to.md",
        "Theory" => "Explanation.md",
        "Reference" => "Reference.md"
    ],
    repo="https://github.com/byuflowlab/FLOWFarm.jl/blob/{commit}{path}#L{line}",
    sitename="FLOWFarm.jl",
    authors="Jared J. Thomas <jaredthomas68@gmail.com>, PJ Stanley <stanley_andrewpj@yahoo.com>",
    doctest=true
)

deploydocs(
    repo = "github.com/byuflowlab/FLOWFarm.jl.git",
    versions = nothing,
    devbranch = "develop"
)


text:
module TestLevel

using FactCheck
using Log4jl
import Log4jl.Level

facts("Levels") do
    context("are comparable") do
        @fact Level.OFF --> less_than(Level.ALL)
        @fact Level.TRACE --> less_than(Level.ALL)
        @fact Level.OFF --> less_than(Level.TRACE)
    end
    context("suport custom values") do
        lvl = :DIAG
        vlvl = Int32(550)
        Level.add(lvl, vlvl)
        @fact_throws AssertionError Level.add(lvl, vlvl)
        @fact Level.OFF --> less_than(Level.DIAG)
        @fact Level.DIAG --> less_than(Level.TRACE)
    end
end

end


text:
using StaticArrays, Distributions, Setfield

struct MCTSNode
    A::UInt8
    B::UInt8
    state::DynamicState
    move::SVector{3, UInt8}
    index::UInt16
    parent::UInt16
    dec_children::SMatrix{8, 8, UInt16, 64}
    chance_children::SVector{2, UInt16}
    σ₁::SVector{8, Float64}
    x₁::SVector{8, Float64}
    n₁::SVector{8, UInt32}
    σ₂::SVector{8, Float64}
    x₂::SVector{8, Float64}
    n₂::SVector{8, UInt32}
end

get_σ(decisions::UInt8, decision::UInt8) =
    is_possible(decisions, decision) ? (1 / Base.ctpop_int(decisions)) : 0.0

function MCTSNode(s::DynamicState, static_s::StaticState, parent::UInt16,
    index::UInt16, move::SVector{3, UInt8})
    A, B = get_possible_decisions(s, static_s)
    return MCTSNode(
        A,
        B,
        s,
        move,
        index,
        parent,
        (@SMatrix zeros(UInt16, 8, 8)),
        (@SVector zeros(UInt16, 2)),
        (@SVector [get_σ(A, i) for i = 0x01:0x08]),
        (@SVector zeros(8)),
        (@SVector zeros(UInt16, 8)),
        (@SVector [get_σ(B, i) for i = 0x01:0x08]),
        (@SVector zeros(8)),
        (@SVector zeros(UInt16, 8)),
    )
end

function chance_node_in_tree(
    tree::SizedVector{20000, MCTSNode, Vector{MCTSNode}}, s::MCTSNode,
    i::UInt16, static_s::StaticState)
    state_1, state_2 = s.state, s.state
    odds = Int8(100)

    chance = get_chance(s.state)
    if chance == 0x05
        odds = Int8(50)
        state_1 = get_chance_state_1(s.state, static_s, chance)
        state_2 = get_chance_state_2(s.state, static_s, chance)
    elseif chance != 0x00
        active1, active2 = get_active(s.state)
        agent = chance < 0x03 ? 0x01 : 0x02
        odds = isodd(chance) ?
            static_s[agent][active1].charged_move_1.buffChance :
            static_s[agent][active2].charged_move_2.buffChance
        state_1 = get_chance_state_1(s.state, static_s, chance)
        state_2 = get_chance_state_2(s.state, static_s, chance)
    end

    new_i = i
    while tree[new_i].index != 0x0000
        new_i += 0x0001
    end
    chance_index = rand() < odds / 100 ? 1 : 2
    if s.chance_children[chance_index] == 0x0000
        tree[new_i] = MCTSNode(chance_index == 1 ? state_1 : state_2, static_s,
            s.index, new_i, @SVector [0x00, 0x00, UInt8(chance_index)])
        s = @set s.chance_children[chance_index] = new_i
    end
    tree[s.index] = s
    u₁ = MCTS!(tree, tree[s.chance_children[chance_index]], new_i, static_s)
    parent = tree[s.parent]
    parent = @set parent.n₁[Int64(s.move[1])] += 1
    parent = @set parent.n₂[Int64(s.move[2])] += 1
    parent = @set parent.x₁[Int64(s.move[1])] += u₁
    parent = @set parent.x₂[Int64(s.move[2])] -= u₁
    tree[parent.index] = parent
    return u₁
end

function fill_missing_children(
    tree::SizedVector{20000, MCTSNode, Vector{MCTSNode}}, s::MCTSNode,
    i::UInt16, static_s::StaticState)
    a = get_decision(s.A, s.B,
        UInt8((length(findall(!iszero, s.dec_children))) %
            Base.ctpop_int(s.A) + 1),
        UInt8((length(findall(!iszero, s.dec_children))) ÷
            Base.ctpop_int(s.A) + 1))
    state = play_turn(s.state, static_s, (a[1], a[2]))
    new_i = i
    while tree[new_i].index != 0x0000
        new_i += 0x0001
    end
    tree[new_i] = MCTSNode(state, static_s, s.index, new_i,
        @SVector [a[1], a[2], 0x00])
    u₁ = Base.cmp(0.5, play_battle(tree[new_i].state, static_s))
    tree[s.index] = @set s.dec_children[Int64(a[1]), Int64(a[2])] = new_i
    update_mcts!(tree, tree[new_i], a[1], a[2], u₁)
    return u₁
end

function select_mcts(tree::SizedVector{20000, MCTSNode, Vector{MCTSNode}},
    s::MCTSNode, γ::Float64 = 0.1)
    η = γ / Base.ctpop_int(s.A)
    max_x1, max_x2 = -Inf, -Inf
    for i = 1:8
        if !iszero(s.A & 0x02^(i-0x01))
            if s.x₁[i] > max_x1
                max_x1 = max(max_x1, s.x₁[i])
            end
        end
        if !iszero(s.B & 0x02^(i-0x01))
            if s.x₂[i] > max_x2
                max_x2 = max(max_x2, s.x₂[i])
            end
        end
    end
    w₁(i) = s.x₁[i] - max_x1
    divisor = mapreduce(x ->
        (iszero(s.A & 2^(x-1)) ? 0.0 : exp(η * w₁(x))), +, 0x01:0x08)
    tree[s.index] = @set s.σ₁ = @SVector [
        (iszero(s.A & 0x02^(i-0x01)) ? 0.0 :
        (((1 - γ) * exp(η * w₁(i))) / divisor) +
        γ / Base.ctpop_int(s.A)) for i = 0x01:0x08]
    η = γ / Base.ctpop_int(s.B)
    w₂(i) = s.x₂[i] - max_x2
    divisor = mapreduce(x -> iszero(s.B & 2^(x-1)) ? 0.0 :
        exp(η * w₂(x)), +, 1:0x08)
    tree[s.index] = @set s.σ₂ = @SVector [(iszero(s.B & 2^(i-1)) ? 0.0 :
        ((1 - γ) * exp(η * w₂(i))) / divisor +
        γ / Base.ctpop_int(s.B)) for i = 0x01:0x08]

    return UInt8(rand(Categorical(s.σ₁))), UInt8(rand(Categorical(s.σ₂)))
end

function update_mcts!(tree::SizedVector{20000, MCTSNode, Vector{MCTSNode}},
    s::MCTSNode, a₁::UInt8, a₂::UInt8, u₁::Int64)
    iszero(s.parent) && return
    parent = tree[s.parent]
    if a₁ != 0x0000 && a₂ != 0x0000
        parent = @set parent.n₁[Int64(a₁)] += 1
        parent = @set parent.n₂[Int64(a₂)] += 1
        parent = @set parent.x₁[Int64(a₁)] += u₁ / parent.σ₁[Int64(a₁)]
        parent = @set parent.x₂[Int64(a₂)] -= u₁ / parent.σ₂[Int64(a₂)]
        tree[parent.index] = parent
    end
    update_mcts!(tree, tree[s.parent], tree[s.parent].move[1],
        tree[s.parent].move[2], u₁)
end

function MCTS!(tree::SizedVector{20000, MCTSNode, Vector{MCTSNode}},
    s::MCTSNode, i::UInt16, static_s::StaticState)
    if Base.ctpop_int(s.A) == 0x00 || Base.ctpop_int(s.B) == 0x00
        score = battle_score(s.state, static_s)
        return score > 0.5 ? 1 : score < 0.5 ? -1 : 0
    elseif tree[s.index] == s && get_chance(s.state) != 0x0000
        return chance_node_in_tree(tree, s, i, static_s)
    elseif tree[s.index] == s &&
        length(findall(!iszero, s.dec_children)) !=
        Base.ctpop_int(s.A) * Base.ctpop_int(s.B)
        return fill_missing_children(tree, s, i, static_s)
    else
        a = select_mcts(tree, s)
        s₁ = tree[s.dec_children[a[1], a[2]]]
        u₁ = MCTS!(tree, s₁, i, static_s)
        update_mcts!(tree, s₁, a[1], a[2], u₁)
        return u₁
    end
end

const default_node = MCTSNode(0x00, 0x00,
        DynamicState(StaticState(["mew" "mew" "mew" "mew" "mew" "mew"])),
        (@SVector zeros(UInt8, 3)), 0x0000, 0x0000,
        (@SMatrix [0x0000 for i = 1:8, j = 1:8]), (@SVector zeros(UInt16, 2)),
        (@SVector zeros(8)), (@SVector zeros(8)), (@SVector zeros(UInt16, 8)),
        (@SVector zeros(8)), (@SVector zeros(8)), (@SVector zeros(UInt16, 8)))

function delete_nodes!(tree::SizedVector{20000, MCTSNode, Vector{MCTSNode}},
    index::UInt16, to_keep::UInt16)
    for i = 1:8, j = 1:8
        if tree[index].dec_children[i, j] != 0x0000 &&
            tree[index].dec_children[i, j] != to_keep
            delete_nodes!(tree, tree[index].dec_children[i, j], to_keep)
        end
    end
    for i = 1:2
        if tree[index].chance_children[i] != 0x0000 &&
            tree[index].chance_children[i] != to_keep
            delete_nodes!(tree, tree[index].chance_children[i], to_keep)
        end
    end
    tree[index] = default_node
end

function update_tree_MCTS!(
    tree::SizedVector{20000, MCTSNode, Vector{MCTSNode}}, curr_index::UInt16,
    static_s::StaticState)

    for i = 0x0001:0x2710
        if tree[i].index == 0x0000
            MCTS!(tree, tree[curr_index], i, static_s)
        end
    end
end

function select_decisions_MCTS(dynamic_state::DynamicState,
    static_s::StaticState)
    tree = SizedVector{20_000}([default_node for i = 0x0001:0x4e20])
    tree[0x0001] = MCTSNode(dynamic_state, static_s, 0x0000, 0x0001,
        @SVector [0x00, 0x00, 0x00])
    curr_index = 0x0001
    strat = Strategy([], [], [], [])

    while Base.ctpop_int(tree[curr_index].A) != 0x00 &&
        Base.ctpop_int(tree[curr_index].B) != 0x00
        update_tree_MCTS!(tree, curr_index, static_s)
        d1, d2 = select_mcts(tree, tree[curr_index])
        prev_index = curr_index
        curr_index = tree[curr_index].dec_children[d1, d2]

        active = get_active(tree[curr_index].state)
        chance = get_chance(tree[curr_index].state)
        if chance == 0x05
            curr_index = rand() < 0.5 ? tree[curr_index].chance_children[1] :
                tree[curr_index].chance_children[2]
        elseif chance != 0x00
            agent = chance < 0x03 ? 0x01 : 0x02
            move = isodd(chance) ?
                static_s[agent][active[agent]].charged_move_1 :
                static_s[agent][active[agent]].charged_move_2
            curr_index = rand() < move.buffChance / 100 ?
                tree[curr_index].chance_children[1] :
                tree[curr_index].chance_children[2]
        end

        push!(strat.decisions, (d1, d2))
        push!(strat.scores, battle_score(tree[curr_index].state, static_s))
        push!(strat.activeMons, active)
        push!(strat.hps,
            ((get_hp(tree[curr_index].state[0x01][0x01]),
            get_hp(tree[curr_index].state[0x01][0x02]),
            get_hp(tree[curr_index].state[0x01][0x03])),
            (get_hp(tree[curr_index].state[0x02][0x01]),
            get_hp(tree[curr_index].state[0x02][0x02]),
            get_hp(tree[curr_index].state[0x02][0x03]))))

        delete_nodes!(tree, prev_index, curr_index)
        new_s = tree[curr_index]
        tree[curr_index] = @set new_s.parent = 0x0000
    end

    return strat
end


text:
module IDLCall

using Compat

# Find IDL library directory if on Linux
if is_unix()
    idl_dir = dirname(chomp(readstring(`which idl`)))
    idl_lib_dir = chomp(readstring(`bash -c "ls -d $(idl_dir)/bin.*"`))
    const idlcall = idl_lib_dir*"/libidl"
    const idlrpc = idl_lib_dir*"/libidl_rpc"
else
    const idlcall = "libidl"
    const idlrpc = "libidl_rpc"
end

export init, get_var, put_var, execute, @get_var, @put_var, idl_repl

jl_idl_type = get(ENV, "JL_IDL_TYPE", is_windows() ? "CALLABLE" : "RPC")

jl_idl_type == "RPC" ? include("IDLRPC.jl") :
jl_idl_type == "CALLABLE" ? include("IDLCallable.jl") :
is_windows() ? error("JL_IDL_TYPE must be CALLABLE on windows") :
error("JL_IDL_TYPE must be RPC or CALLABLE")

include("IDLREPL.jl")

init()
idl_repl()
repl = idl_repl

end


text:
# space.jl

"""
A domain is always a subset of a geometric space. A geometric space is completely
characterized by the type `T` of its elements. There are no restrictions on
concrete values: each instance of type `T` belongs to `GeometricSpace{T}`.
Conversely, each numeric type `T` can give rise to only one geometric space.

The goal of `GeometricSpace` is to provide a lightweight framework for converting
elements between spaces. Two features are embeddings and isomorphisms. It is
lightweight because the space is characterized by the type and has no state, hence
it can be inferred from any variable based only on its type. A restriction that
results from this choice is that there can only be one type of embedding from one
space to another. More general embeddings could be implemented using domains or
specific maps.

Embeddings:
A space `A` is embedded in a space `B` if each element of `A` corresponds to an
element of `B`. In that case an element of `A` can be promoted (using `promote_space`)
to an element of type `B`. This is not necessarily the same as the promotion of
the type of `A` to the type of `B`. For example, a scalar can be embedded into a
two-dimensional space (say with the second component equal to zero). It would be
undesirable to implement that using the standard Julia promotion system, since
it would apply to all Julia code used concurrently.

Isomorphisms:
If `A` is embedded in `B` and `B` is embedded in `A` than they `A` and `B` are
isomorphic. In that case, an element of `A` can be converted to an element of
`B` and vice-versa. The conversion is invertible. One example is the isomorphism
between ℝ^2 and ℂ.

New spaces can be defined by defining a new numeric type. Embeddings and
isomorphisms are defined by defining new `embedding_reduction` and
`isomorphism_reduction` rules, along with conversiong using `convert_space`
and `restrict_space`. See the documentation of these functions for
information on how to use them.
"""
struct GeometricSpace{T}
end

# An internal shorthand in order to avoid excessively long lines when there
# are several arguments of type GeometricSpace :-)
# We export a longer name in order to avoid name clashes with other Space types
const GSpace = GeometricSpace

"AnySpace is the superset of all possible geometric spaces."
const AnySpace = GeometricSpace{Any}

eltype(::GeometricSpace{T}) where {T} = T
eltype(::Type{GeometricSpace{T}}) where {T} = T
eltype(::Type{S}) where {S <: GeometricSpace} = eltype(supertype(S))

"""
Some geometric spaces have an eltype that is composed in terms of an underlying
eltype. For example, `Complex{Float64}` and `SVector{2,Float64}` are based on
`Float64`, which is their `subeltype`.
"""
subeltype(::Type{GeometricSpace{T}}) where {T} = T

"""
Create a space that is similar to the given space, but with a different `subeltype`.

For example, `similar_space(ComplexSpace{Float64}, BigFloat)` yields a
`ComplexSpace{BigFloat}`.
"""
similar_space(::Type{GeometricSpace{T}}, ::Type{S}) where {T,S} = GeometricSpace{S}

isreal(space::GeometricSpace{T}) where {T} = isreal(T)
isreal(::Type{GeometricSpace{T}}) where {T} = isreal(T)

# Return the zero element
zero(space::GeometricSpace{T}) where {T} = zero(T)
zero(::Type{GeometricSpace{T}}) where {T} = zero(T)

"The origin of a space is its zero element."
origin(space::GeometricSpace) = zero(space)

# Definition of element membership is strictly based on type:
in(x, ::GeometricSpace) = false
in(x, ::Type{GeometricSpace{T}}) where {T} = false
in(x::T, ::Type{GeometricSpace{T}}) where {T} = true

# spaces cannot be empty
isempty(::GeometricSpace) = false

# Make the space bigger by widening the numeric type
widen(A::Type{GeometricSpace{T}}) where {T} = _widen(A, eltype(A), subeltype(A))
# eltype and subeltype are the same: widen those
_widen(A, ::Type{T}, ::Type{T}) where {T} = GeometricSpace{widen(T)}
# eltype and subeltype are different: widen the subeltype
_widen(A, ::Type{T}, ::Type{S}) where {T,S} = similar_space(A, widen(S))

"Return the geometric space type with eltype `T`."
spacetype(::Type{T}) where {T} = GeometricSpace{T}

"Return the geometric space of all elements with the same type as `x`."
spaceof(x::T) where {T} = spacetype(T)
# For convenience
spaceof(x::GeometricSpace{T}) where {T} = typeof(x)

"""
Return the superspace of the given geometric space `A`. The superspace of `A`
should be larger than `A`, and `A` should be embedded in it.

Superspaces are used to automatically discover embeddings and promotion rules.
Its role is analogous to that of `supertype` in the Julia type system. Indeed,
`AnySpace` is a superspace of all spaces, much like `Any` is a supertype of
all types.
"""
superspace(::Type{GeometricSpace{T}}) where {T} = AnySpace

superspaceof(x) = superspace(spaceof(x))

"""
A space `A` is a subspace of space `B` if `B` is a supertype of `A`.
"""  # Any space is a subspace of itself
issubspace(A::Type{GeometricSpace{T}}, B::Type{GeometricSpace{T}}) where {T} = true
# So is AnySpace. This is covered by the line above, but we make it explicit for clarity
issubspace(A::Type{AnySpace}, B::Type{AnySpace}) = true
# AnySpace is not a subspace of any other space
issubspace(A::Type{AnySpace}, B::Type{GeometricSpace{T}}) where {T} = false
# In all other cases, we replace A by its superspace and recurse
issubspace(A::Type{GeometricSpace{T}}, B::Type{GeometricSpace{S}}) where {T,S} = issubspace(superspace(A), B)
# TODO add issubspace to specific spaces


text:
struct Node
    X::Float64
    Y::Float64
end

struct Support
    Node::Int
    XRelease::Int
    YRelease::Int
end

struct Material
    E::Float64
end

struct Section
    I::Float64
end

struct Element
    NodeB::Int
    NodeE::Int
    Material
    Section
end

struct CLoad
    value::Float64
    distance::Float64
end

struct MLoad
    value::Float64
    distance::Float64
end

struct NLoad{T}
    loads::T
end

struct DLoad
    value1::Float64
    value2::Float64
    distance1::Float64
    distance2::Float64
end

struct ALoad
    value::Float64
    distance::Float64
end

struct SurfaceShearLoad
    value::Float64
    distance1::Float64
    distance2::Float64
end

function gK(Xb, Yb, Xe, Ye, E, A)
  δx = Xe - Xb
  δy = Ye - Yb
  L = sqrt(δx^2 + δy^2)
  _cos = δx/L
  _sin = δy/L
  _cos² = _cos * _cos
  _sin² = _sin * _sin
  _cos_sin = _cos * _sin
  (A*E/L)*[_cos² _cos_sin -_cos² -_cos_sin
   _cos_sin _sin² -_cos_sin -_sin²
   -_cos² -_cos_sin _cos² _cos_sin
   -_cos_sin -_sin² _cos_sin _sin²]
end

function T(Xb, Yb, Xe, Ye)
    δx = Xe - Xb
    δy = Ye - Yb
    L = sqrt(δx^2 + δy^2)
    _cos = δx/L
    _sin = δy/L

    [_cos _sin 0 0
     -_sin _cos 0 0
     0 0 _cos _sin
     0 0 -_sin _cos]
end

function k(Xb, Yb, Xe, Ye, E, I)
    δx = Xe - Xb
    δy = Ye - Yb
    L = sqrt(δx^2 + δy^2)

    (E*I/L^3)*[12. 6L -12 6L
               6L 4L^2 -6L 2L^2
               -12 -6L 12 -6L
               6L 2L^2 -6L 4L^2]
end

function Qf(Xb, Yb, Xe, Ye, cLoad::CLoad)
    δx = Xe - Xb
    δy = Ye - Yb
    L = sqrt(δx^2 + δy^2)
    W = cLoad.value
    l₁ = cLoad.distance
    l₂ = L - cLoad.distance

    FSb = (W*l₂^2/L^3)*(3l₁ + l₂)
    FMb = (W*l₁*l₂^2/L^2)
    FSe = (W*l₁^2/L^3)*(l₁ + 3l₂)
    FMe = (-W*l₁^2*l₂/L^2)

    [FSb; FMb; FSe; FMe]
end

function Qf(Xb, Yb, Xe, Ye, dLoad::DLoad)
    δx = Xe - Xb
    δy = Ye - Yb
    L = sqrt(δx^2 + δy^2)
    w₁ = dLoad.value1
    w₂ = dLoad.value2
    l₁ = dLoad.distance1
    l₂ = L - dLoad.distance2

    isUniform = w₁ == w₂

    FSb, FMb, FSe, FMe, = (0., 0., 0., 0.)

    if isUniform
        w = w₁

        FSb = (w*L/2)*(1 - (l₁/L^4)*(2L^3 - 2l₁^2*L + l₁^3) - (l₂^3/L^4)*(2L - l₂))
        FMb = (w*L^2/12)*(1 - (l₁^2/L^4)*(6L^2 - 8l₁*L + 3l₁^2) - (l₂^3/L^4)*(4L - 3l₂))
        FSe = (w*L/2)*(1 - (l₁^3/L^4)*(2L - l₁) - (l₂/L^4)*(2L^3 - 2l₂^2 * L + l₂^3))
        FMe = (-w*L^2/12)*(1 - (l₁^3/L^4)*(4L - 3l₁) - (l₂^2/L^4)*(6L^2 - 8l₂*L + 3l₂^2))
    else
        FSb = w₁*(L-l₁)^3/(20*L^3)*( (7L + 8l₁) - (l₂*(3L + 2l₁)/(L-l₁))
              *(1 + (l₂/(L-l₁)) + (l₂^2/(L-l₁)^2)) + (2l₂^4/(L-l₁)^3))
              +(w₂*(L-l₁)^3/20L^3)*((3L + 2l₁)*(1 + l₂/(L - l₁)
                + l₂^2/(L-l₁)^2) - (l₂^3/(L-l₁)^2)*(2 + (15L - 8l₂)/(L-l₁)))

        FMb = (w₁*(L-l₁)^3/60L^2)*(3(L+4l₁) - (l₂*(2L + 3l₁)/(L-l₁))*
               (1 + (l₂/(L-l₁)) + (l₂^2/(L-l₁)^2)) + (3l₂^4/(L-l₁)^3))
               + (w₂*(L-l₁)^3/60L^2)*((2L + 3l₁)*(1 + (l₂/(L-l₁))
               + (l₂^2/(L-l₁)^2)) - (3l₂^3/(L-l₁)^2)*(1 + (5L-4l₂)/(L-l₁)))

        FSe = ((w₁ + w₂)/2)*(L - l₁ - l₂) - FSb

        FMe = ((L- l₁ - l₂)/60)*(w₁*(-2L + 2l₁ - l₂) - w₂*(L - l₁ + 2l₂))
               + FSb*L - FMb

    end
    [FSb; FMb; FSe; FMe]
end

function Qf(Xb, Yb, Xe, Ye, loads::Array{CLoad,1})
    sum(map(load -> Qf(Xb, Yb, Xe, Ye, load), loads))
end

function Qf(Xb, Yb, Xe, Ye, loads::Array{DLoad,1})
    sum(map(load -> Qf(Xb, Yb, Xe, Ye, load), loads))
end

function Qf(Xb, Yb, Xe, Ye, loads::Array{Any,1})
    if size(loads, 1) == 0
        return [0.0; 0.0; 0.0; 0.0]
    else
        sum(map(load -> Qf(Xb, Yb, Xe, Ye, load), loads))
    end
end

function nsc(COORD, MSUP, NCJT)
    #Determine number of joints (NJ)
    NJ = size(COORD, 1)

    #Determine Number of Restraints (NR)
    NR = sum(map(s -> s.XRelease + s.YRelease, MSUP))

    numDOF = NJ * NCJT - NR

    NSC = Array{Int}(2 * NJ)                                  # Create the NSC with uninitialized values

    freeCoord = 1
    restraintCoord = numDOF + 1
    for nodeIndex in 1:size(COORD, 1)                           # Loop through all of the COORD
        supportIndexes = [s.Node for s in MSUP]
        if in(nodeIndex, supportIndexes)                        #   if this node is a support then...
            supportIndex = findfirst(supportIndexes, nodeIndex) #     get the support Index for this node
            if MSUP[supportIndex].XRelease == 1             #     if node is restrained in the X direction...
                 NSC[2*nodeIndex - 1] = restraintCoord           #       NSC X value for this node is the next restraint coordinate
                restraintCoord += 1
            else                                                #     otherwise, this node is free in the X direction
                 NSC[2*nodeIndex - 1] = freeCoord         #       NSC X value for this node is the next free coordinate
                freeCoord += 1
            end
            if MSUP[supportIndex].YRelease == 1             #     if node is restrained in the Y direction...
                NSC[2*nodeIndex] = restraintCoord         #       NSC Y value for this node is the next restraint coordinate
                restraintCoord += 1
            else                                                #     otherwise, this node is free in the Y direction
                 NSC[2*nodeIndex] = freeCoord             #       NSC Y value for this node is the next free coordinate
                freeCoord += 1
            end
        else                                                    #   otherwise, this node is free
             NSC[2*nodeIndex - 1] = freeCoord             #     NSC X value for this node is the next free coordinate
             NSC[2*nodeIndex] = freeCoord + 1             #     NSC Y value for this node is the next free coordinate
            freeCoord += 2
        end
    end
    NSC
end

function gs(COORD, MSUP, MPRP, EM, CP, NSC, NCJT)
    #Determine number of joints (NJ)
    NJ = size(COORD, 1)

    #Determine Number of Restraints (NR)
    NR = sum(map(s -> s.XRelease + s.YRelease, MSUP))

    numDOF = NJ * NCJT - NR

    S = zeros(numDOF, numDOF)

    for memberIndex in 1:size(MPRP, 1)
        memberInfo = MPRP[memberIndex]
        nodeBegin, nodeEnd = memberInfo.NodeB, memberInfo.NodeE
        Xb = COORD[nodeBegin].X
        Yb = COORD[nodeBegin].Y
        Xe = COORD[nodeEnd].X
        Ye = COORD[nodeEnd].Y
        E = EM[memberInfo.Material].E
        I = CP[memberInfo.Section].I
        GK = k(Xb, Yb, Xe, Ye, E, I)

        NSC_Indecies =
        vcat([NSC[(nodeBegin - 1)*NCJT + x] for x in 1:NCJT],
             [NSC[(nodeEnd - 1)*NCJT + x] for x in 1:NCJT])

        iCount = 0
        for i in NSC_Indecies
            if i == NSC_Indecies[1]
                iCount = 0
            end
            iCount += 1
            if i <= numDOF
                jCount = 0
                for j in NSC_Indecies
                    jCount +=1
                    if j <= numDOF
                         S[j,i] = S[j,i] + GK[jCount, iCount]
                    end
                end
            end
        end
    end
    S
end

function p_pf(COORD, MSUP, NSC, memberLoads, nodeLoads, NCJT)
    NJ = size(COORD, 1)
    NR = sum(map(s -> s.XRelease + s.YRelease, MSUP))
    numDOF = NJ * NCJT - NR
    NML = length(memberLoads)
    NJL = length(nodeLoads)

    P = zeros(numDOF)::Array{Float64,1}
    Pf = zeros(numDOF)::Array{Float64,1}

    dofIndecies = find(x -> x <= numDOF, NSC)

    for i in 1:size(dofIndecies,1)

        dofIndex = dofIndecies[i]
        remainder = dofIndex%NCJT
        dofOfNode = (remainder == 0 ? NCJT : remainder)
        node = Int(1 + (dofIndex - dofOfNode)/NCJT)

        _nodeLoads =  map(load -> load[2], filter(load -> load[1] == node, nodeLoads))

         for load in _nodeLoads
            P[i] = P[i] + load.loads[dofOfNode]
         end

        for j in 1:size(MPRP,1)
            member = MPRP[j]
            nodeBegin, nodeEnd = member.NodeB, member.NodeE
            Xb = COORD[nodeBegin].X
            Yb = COORD[nodeBegin].Y
            Xe = COORD[nodeEnd].X
            Ye = COORD[nodeEnd].Y

            _memberLoads = map(load -> load[2], filter(load -> load[1] == j, memberLoads))
            _Qf = Qf(Xb, Yb, Xe, Ye, _memberLoads)

            if nodeBegin == node
                Pf[i] = Pf[i] + _Qf[dofOfNode]
            end
            if nodeEnd == node
                Pf[i] = Pf[i] + _Qf[2 + dofOfNode]
            end
        end
    end
    P, Pf
end

function get_d(S, P, Pf)
    inv(S)*(P - Pf)
end

function forces_reactions(COORD, NSC, MSUP, EM, CP, MPRP, d, NCJT)
    NJ = size(COORD, 1)
    NR = sum(map(s -> s.XRelease + s.YRelease, MSUP))
    numDOF = NJ * NCJT - NR

    R = zeros(NR)
    forces = Array{Array{Float64,1}}(size(MPRP, 1))

    for member in 1:size(MPRP,1)
        memberInfo = MPRP[member]
        nodeBegin, nodeEnd = memberInfo.NodeB, memberInfo.NodeE
        Xb = COORD[nodeBegin].X
        Yb = COORD[nodeBegin].Y
        Xe = COORD[nodeEnd].X
        Ye = COORD[nodeEnd].Y
        E = EM[memberInfo.Material].E
        I = CP[memberInfo.Section].I

        v = zeros(NCJT*2)

        NSC_Indecies =
        vcat([NSC[(nodeBegin - 1)*NCJT + x] for x in 1:NCJT],
             [NSC[(nodeEnd - 1)*NCJT + x] for x in 1:NCJT])

        for i in 1:(NCJT*2)
            NSC_Index = NSC_Indecies[i]
            if NSC_Index <= numDOF
                v[i] = d[NSC_Index]
            end
        end

        _T = T(Xb, Yb, Xe, Ye)

        U = _T * v

        BK = k(Xb, Yb, Xe, Ye, I, E)

        _memberLoads = map(load -> load[2], filter(load -> load[1] == member, memberLoads))
        _Qf = Qf(Xb, Yb, Xe, Ye, _memberLoads)

        Q = BK * U + _Qf
        forces[member] = Q

        F = _T'Q

        for i in 1:(NCJT*2)
            nscIndex =
                if i<=NCJT
                    (nodeBegin -1)*NCJT + i
                else
                    (nodeEnd - 1)*NCJT + (i - NCJT)
                end
            N = NSC[nscIndex]
            if N > numDOF
                R[N - numDOF] = R[N - numDOF] + F[i]
            end
        end
    end
    forces, R
end

function run_analysis(COORD, MSUP, EM, CP, MPRP, memberLoads, nodeLoads)
    NCJT = 2
    NSC = nsc(COORD, MSUP, NCJT)
    GS = gs(COORD, MSUP, MPRP, EM, CP, NSC, NCJT)
    P, Pf = p_pf(COORD, MSUP, NSC, memberLoads, nodeLoads, NCJT)
    d = get_d(GS, P, Pf)
    forces, reactions = forces_reactions(COORD, NSC, MSUP, EM, CP, MPRP, d, NCJT)
end

max_stars_repo_name: Julia-SE/matrix-analysis

max_issues_repo_name: Julia-SE/matrix-analysis

max_forks_repo_name: Julia-SE/matrix-analysis


text:
function run_ed(hours::UnitRange{Int64}, args...)

end

max_issues_repo_name: bobbyxiong/redispatch-ptg

max_forks_repo_name: bobbyxiong/redispatch-ptg


text:
###############################################################################
#
# Optimize flow rates in the sunlit and shaded layers
#
###############################################################################
"""
    optimize_flows!(
                node::SPACSimple{FT}
    ) where {FT<:AbstractFloat}

Optimize the flow rates in sunlit and shaded layers, given
- `node` [`SPACSimple`] type struct
"""
function optimize_flows!(
            node::SPACSimple{FT}
) where {FT<:AbstractFloat}
    # unpack required parameters
    @unpack frac_sh, frac_sl = node.container2L;

    # calculate the ecrit
    node.ec = critical_flow(node.hs, node.ec);
    ec_sl   = node.ec * frac_sl;
    ec_sh   = node.ec * frac_sh;

    # optimize the sunlit and shade layers
    f_sl = min(ec_sl / FT(1.01), node.opt_f_sl);
    f_sh = min(ec_sh / FT(1.01), node.opt_f_sh);

    ms = ReduceStepMethodND{FT}(x_mins = FT[0,0], x_maxs = [ec_sl,ec_sh], x_inis = [f_sl, f_sh], Δ_inis = FT[0.1,0.1]);
    st = SolutionToleranceND{FT}(FT[9e-4, 9e-4], 50);

    @inline f(x) = (
        leaf_gas_exchange_nonopt!(node, x[1], x[2]);
        return node.containerOP
    );

    fs = find_peak(f, ms, st);

    # update the optimal flow rates
    node.opt_f_sl = fs[1];
    node.opt_f_sh = fs[2];

    return nothing
end


text:
using BinaryBuilder

name = "ITK"
version = v"5.0.1"

# Collection of sources required to build zlib
sources = [
    "https://github.com/InsightSoftwareConsortium/ITK/releases/download/v$(version)/InsightToolkit-$(version).tar.gz" =>
    "613b125cbf58481e8d1e36bdeacf7e21aba4b129b4e524b112f70c4d4e6d15a6",
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/InsightToolkit-*
mkdir build && cd build

if [[ "${target}" == i686-linux-* ]]; then
    export EMU=
elif [[ "${target}" == x86_64-linux-* ]]; then
    export EMU=
elif [[ "${target}" == arm-linux-* ]]; then
    export EMU=
elif [[ "${target}" == powerpc64le-linux-* ]]; then
    export EMU=
elif [[ "${target}" == x86_64-apple-* ]]; then
    export EMU=
elif [[ "${target}" == i686-w64-mingw32 ]]; then
    export EMU=
elif [[ "${target}" == x86_64-w64-mingw32 ]]; then
    export EMU=
elif [[ "${target}" == *freebsd* ]]; then
    export EMU=
elif [[ "${target}" == aarch64-linux-* ]]; then
    export EMU=
else
    export EMU=
fi

cmake ../ -DCMAKE_CROSSCOMPILING_EMULATOR=${EMU} -DBUILD_EXAMPLES=OFF -DBUILD_TESTING=OFF -DCMAKE_BUILD_TYPE=Release -DBUILD_SHARED_LIBS=ON -DCMAKE_INSTALL_PREFIX=$prefix -DCMAKE_TOOLCHAIN_FILE=/opt/$target/$target.toolchain

cmake -C TryRunResults.cmake ../
make install -j${nproc}
"""

# Build for ALL THE PLATFORMS!
platforms = supported_platforms()

# The products that we will ensure are always built
products = prefix -> [
    LibraryProduct(prefix, "libitk", :libitk),
]

# Dependencies that must be installed before this package can be built
dependencies = [
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: IanButterworth/ITKBuilder

max_issues_repo_name: IanButterworth/ITKBuilder

max_forks_repo_name: IanButterworth/ITKBuilder


text:
using PyPlot, AcousticWave
# =============generate synthetic data set================================
nz = 201; nx = 201; ext = 20; iflag = 1; dz = 5.; dx = 5.; dt = 1e-3;
v  = 3000. * ones(nz, nx);
# initialize source term
pos_src = [101 75; 101 125]; f0 = 30.0; ot = [0.0, 0.01];
src = InitSources(pos_src, f0, ot, dt);
path = "/Users/wenlei/Desktop/srcSpt.bin"
src2srcSpts(path, src, nz, nx)
p = readSrcSpt(path, 47); p = reshape(p, nz, nx);

vmax = maximum(v); vmin = minimum(v);
fidMtx = InitFidMtx(nz, nx, ext, iflag, dz, dx, dt, vmax, vmin, f0, v);
iz = collect(5:5:201); ix = 20 * ones(Int64, length(iz));
pos_rec = hcat(iz, ix);
shot = MultiStepForward(pos_rec, src, fidMtx, tmax=0.5);

max_stars_repo_name: WneleiGao/AcousticWaveLab

max_issues_repo_name: WneleiGao/AcousticWaveLab

max_forks_repo_name: WneleiGao/AcousticWaveLab


text:
# Array of structs

using BenchmarkTools, Statistics, CSV

struct TripPayment
    vendor_id::Int
    tpep_pickup_datetime::String
    tpep_dropoff_datetime::String
    passenger_count::Int
    trip_distance::Float64
    fare_amount::Float64
    extra::Float64
    mta_tax::Float64
    tip_amount::Float64
    tolls_amount::Float64
    improvement_surcharge::Float64
    total_amount::Float64
end

# Use CVS.jl to parse the file into a vecot or TripPayment objects
function read_trip_payment_file(file)
    f = CSV.File(file, datarow = 3)
    records = Vector{TripPayment}(undef, length(f))
    for (i, row) in enumerate(f)
        records[i] = TripPayment(row.VendorID,
                                 row.tpep_pickup_datetime,
                                 row.tpep_dropoff_datetime,
                                 row.passenger_count,
                                 row.trip_distance,
                                 row.fare_amount,
                                 row.extra,
                                 row.mta_tax,
                                 row.tip_amount,
                                 row.tolls_amount,
                                 row.improvement_surcharge,
                                 row.total_amount)
    end
    return records
end

records = read_trip_payment_file("yellow_tripdata_2018-12_100k.csv")

#=
julia> records
100000-element Array{TripPayment,1}:
 TripPayment(1, "2018-12-01 00:28:22", "2018-12-01 00:44:07", 2, 2.5, 12.0, 0.5, 0.5, 3.95, 0.0, 0.3, 17.25)
 TripPayment(1, "2018-12-01 00:52:29", "2018-12-01 01:11:37", 3, 2.3, 13.0, 0.5, 0.5, 2.85, 0.0, 0.3, 17.15)
 TripPayment(2, "2018-12-01 00:12:52", "2018-12-01 00:36:23", 1, 0.0, 2.5, 0.5, 0.5, 0.0, 0.0, 0.3, 3.8)
 ⋮
=#

# find mean of fare_amount field

#=
julia> mean(r.fare_amount for r in records)
11.985990599999994
=#

# how much time does it take?

#=
julia> @btime mean(r.fare_amount for r in $records);
  639.603 μs (1 allocation: 16 bytes)
=#

# what if we have the fare_amount data in an array?
fare_amounts = [r.fare_amount for r in records];

#=
julia> @btime mean(fare_amounts)
  27.090 μs (1 allocation: 16 bytes)
=#

# how much faster?
#=
julia> 640 / 27
23.703703703703702
=#


text:
mutable struct Problem
  F::Firm

  N::Int
  S::Int
  p::Matrix{Float64}

  dtol::Float64
  maxiter::Int
end

function updatefirmτ!(Φ::Problem, τ::Vector{Float64})::Problem
  Φ.F.τ .= τ
  refreshR!(Φ.F, Φ.p)
  findpolicy!(Φ::Problem)
  return Φ
end

#gets the maximum by column of a square matrix. Storage space must be provided
function colmaxsquare!(A::AbstractMatrix{Float64},
  V::AbstractVector{Float64},
  Π::AbstractVector{Int})::Nothing

  N::Int = size(A,2)
  @assert length(V) == N
  @assert length(Π) == N
  @assert size(A,1) == N

  local maxv::Float64
  local maxi::Int

  @fastmath @inbounds for c ∈ 1:N

    maxv = A[1,c]
    maxi = 1

    for r ∈ 2:N
      (maxv < A[r,c]) && (maxv = A[r,c]; maxi=r)
    end

    V[c] = maxv
    Π[c] = maxi
  end

  return nothing
end

function findpolicy!(Φ::Problem)::Problem
  local iter::Int = 0
  local derr::Float64 = 1.0
  #local vrow::Matrix{Float64} = Matrix{Float64}(undef,1,Φ.N) #captures updated values
  #local πrow::Matrix{CartesianIndex} = Matrix{CartesianIndex}(undef,1,Φ.N) #captures updated policy
  local vold::Array{Float64,3} = deepcopy(Φ.F.V) #holds the previous value for comparison purposes

  #each view here maps all current states to a single future state
  local rslices = Matrix{SubArray{Float64}}(undef, Φ.S, Φ.S)
  local mslices = Array{SubArray{Float64},3}(undef, Φ.S, Φ.S, Φ.S)
  local vslices = Matrix{SubArray{Float64}}(undef, Φ.S, Φ.S) #will a view for each lagged and current state
  local πslices = Matrix{SubArray{Int}}(undef, Φ.S, Φ.S)

  #make our convenience views
  for sₗ ∈ 1:Φ.S, s ∈ 1:Φ.S
    rslices[sₗ, s] = view(Φ.F.R, :, :, sₗ, s)
    vslices[sₗ, s] = view(Φ.F.V, :, sₗ, s)
    πslices[sₗ, s] = view(Φ.F.Π, :, sₗ, s)
    mslices[sₗ, s,:] .= ((sₜ₁::Int)->view(Φ.F.M, :, :, sₗ, s, sₜ₁)).(1:Φ.S)
  end

  while (derr > Φ.dtol) && (iter < Φ.maxiter)
    for sₗ ∈ 1:Φ.S, s ∈ 1:Φ.S #need to use vec since findmax returns column vectors

      colmaxsquare!(rslices[sₗ, s] .+
        (Φ.p[1, s] .* mslices[sₗ, s,1] .* vslices[s, 1]) .+
        (Φ.p[2, s] .* mslices[sₗ, s,2] .* vslices[s, 2]), vslices[sₗ, s], πslices[sₗ, s])
      #vrow, πrow = findmax(rslices[sₗ, s] .+
      #  (Φ.p[1, s] .* mslices[sₗ, s,1] .* vslices[s, 1]) .+
      #  (Φ.p[2, s] .* mslices[sₗ, s,2] .* vslices[s, 2]), dims=1)

      #need to do this to fix the dimensionality of the returned values from findmax
      #vslices[sₗ, s] .= vec(vrow)
      #πslices[sₗ, s] .= vec(πrow)

    end

    derr = maximum((abs).(vold .- Φ.F.V))
    vold .= Φ.F.V

    iter += 1

    (iter % 500 == 0) && println("iter: $iter, derr: $derr")

  end

  (iter == Φ.maxiter) && @warn("Convergence failure iter: $iter, derr: $derr")

  return Φ
end

function Problem(F::Firm{NTuple{2,Float64}}, p::Matrix{Float64};
    dtol::Float64=error("dtol is required"),
    maxiter::Int = error("maxiter is required"),
    verbose::Bool = length(F.ψ[1]) < 1000, findpolicy::Bool = true)::Problem

  N::Int = length(F.ψ[1])
  S::Int = size(p,1)

  Φ::Problem = Problem(F, N, S, p, dtol, maxiter)
  findpolicy && findpolicy!(Φ)

  #display(Φ.F.V)
  if verbose && findpolicy
    highpolicy::Matrix{Float64} = Matrix{Int}(undef, Φ.N, Φ.S)
    lowpolicy::Matrix{Float64} = Matrix{Int}(undef, Φ.N, Φ.S)

    highvalue::Matrix{Float64} = Matrix{Float64}(undef, Φ.N, Φ.S)
    lowvalue::Matrix{Float64} = Matrix{Float64}(undef, Φ.N, Φ.S)
    for s ∈ 1:Φ.S
      highpolicy[:,s] .= ((n)->(maximum(Φ.F.R[:,n,1,s]) > 0.0) ? (Φ.F.Π[n,1,s]) : -1.0).(1:Φ.N)
      lowpolicy[:,s] .= ((n)->(maximum(Φ.F.R[:,n,2,s]) > 0.0) ? (Φ.F.Π[n,2,s]) : -1.0).(1:Φ.N)

      highvalue[:,s] = ((n)->(maximum(Φ.F.R[:,n,1,s]) > 0.0) ? (Φ.F.V[n,1,s]) : -1.0).(1:Φ.N)
      lowvalue[:,s]  = ((n)->(maximum(Φ.F.R[:,n,2,s]) > 0.0) ? (Φ.F.V[n,2,s]) : -1.0).(1:Φ.N)

      highvalue[:,s] .= (v::Float64->round(v,digits=1)).(highvalue[:,s])
      lowvalue[:,s] .= (v::Float64->round(v,digits=1)).(lowvalue[:,s])
    end

    for n ∈ 1:Φ.N
      k::Float64 = round(Φ.F.ψ[1][n], digits=2)
      yh::Float64 = round(Φ.F.ψ[2][n], digits=2)

      ylgh::Float64 = getyl(Φ.F.ψ[2][n], Φ.F.ψ[1][n], Φ.p[1,1], Φ.F)
      ylgl::Float64 = getyl(Φ.F.ψ[2][n], Φ.F.ψ[1][n], Φ.p[1,2], Φ.F)

      ylgh = round(ylgh, digits=2)
      ylgl = round(ylgl, digits=2)

      println("n: $n; k,yh: ($k, $yh); yl, v,p|sₗ=h: $ylgh, $(highvalue[n,:]), $(highpolicy[n,:])")
      #println("n: $n; k,yh: ($k, $yh); yl, p|sₗ=h: $ylgh, $(highpolicy[n,:]); p|sₗ=l: $ylgl, $(lowpolicy[n,:] )")
      #println(Φ.F.Π[n,1,1])
    end
  end

  return Φ
end

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
using PowerSystems
using NLsolve
const PSY = PowerSystems

include(joinpath(dirname(@__FILE__), "dynamic_test_data.jl"))
include(joinpath(dirname(@__FILE__), "data_utils.jl"))
############### Data Network ########################

raw_file_dir = joinpath(dirname(@__FILE__), "../benchmarks/psse/GENROU/ThreeBusMulti.raw")
sys = System(raw_file_dir);

#Replace Gen101 by Source
remove_component!(ThermalStandard, sys, "generator-101-1");
add_source_to_ref(sys)

function dyn_gen_genrou(generator)
    return PSY.DynamicGenerator(
        name = get_name(generator),
        ω_ref = 1.0, #ω_ref
        machine = machine_genrou(), #machine
        shaft = shaft_genrou(), #shaft
        avr = avr_type1(), #avr
        prime_mover = tg_none(), #tg
        pss = pss_none(),
    ) #pss
end

for l in get_components(PSY.PowerLoad, sys)
    PSY.set_model!(l, PSY.LoadModels.ConstantImpedance)
end

#Add GENROU to System
g = get_component(ThermalStandard, sys, "generator-102-1")
dyn_gen = dyn_gen_genrou(g)
add_component!(sys, dyn_gen, g);