text:
function MCS(G)
    n = nv(G)
    w = zeros(UInt16, n)
    order = zeros(UInt16, n)
    for i = n:-1:1
        temp = w[[order[j]==0 for j=1:n]]
        v = findfirst(x->order[x]==0&&w[x]==maximum(temp), 1:n)
        w[neighbors(G, v)] .+= 1
        order[v] = i
    end
    return order
end

function MCSM!(G)
    n = nv(G)
    w = ones(UInt16, n)
    order = zeros(UInt16, n)
    for i = n:-1:1
        temp = w[[order[j]==0 for j=1:n]]
        maxw = maximum(temp)
        v = findfirst(x->order[x]==0&&w[x]==maxw, 1:n)
        order[v] = i
        unorder = findall(x->order[x]==0,1:n)
        unreach = copy(unorder)
        reach = [UInt16[] for j=1:maxw]
        S = UInt16[]
        for u in intersect(neighbors(G, v), unorder)
            push!(reach[w[u]], u)
            push!(S, u)
            deleteat!(unreach, bfind(unreach, length(unreach), u))
        end
        for k = 1:maxw-1
            while reach[k] != []
                z = pop!(reach[k])
                for u in intersect(neighbors(G, z), unreach)
                    deleteat!(unreach, bfind(unreach, length(unreach), u))
                    if w[u] > k
                        push!(reach[w[u]], u)
                        push!(S, u)
                    else
                        push!(reach[k], u)
                    end
                end
            end
        end
        for u in S
            w[u] += 1
            add_edge!(G, v, u)
        end
    end
    return G,order
end

function Candidate(G, u, v, i, order, n)
    flag = true
    for z in neighbors(G, u)
        if order[z] > i && has_edge(G, v, z) && !has_edge(G, findfirst(x->order[x]==i, 1:n), z)
            flag = false
        end
    end
    return flag
end

function MinimalChordal!(G, order, F)
    n = nv(G)
    for i = n:-1:1
        ind = UInt16[]
        inc = UInt16[]
        for j = 1:size(F[i], 2)
            if Candidate(G, F[i][1,j], F[i][2,j], i, order, n)
                push!(ind, j)
                push!(inc, F[i][1,j], F[i][2,j])
            end
        end
        if length(ind) > 0
            cand = F[i][:,ind]
            unique!(inc)
            sort!(inc)
            linc = length(inc)
            W = complete_graph(linc)
            for j = 1:size(cand, 2)
                rem_edge!(W, bfind(inc,linc,cand[1,j]), bfind(inc,linc,cand[2,j]))
            end
            keep,_ = MCSM!(W)
            for j = 1:size(cand,2)
                if !has_edge(keep, bfind(inc,linc,cand[1,j]), bfind(inc,linc,cand[2,j]))
                    rem_edge!(G, cand[1,j], cand[2,j])
                end
            end
        end
    end
    return G,MCS(G)
end

function bfind(A, l, a)
    if l == 0
        return 0
    end
    low = 1
    high = l
    while low <= high
        mid = Int(ceil(1/2*(low+high)))
        if A[mid] == a
           return mid
        elseif A[mid] < a
           low = mid + 1
        else
           high = mid - 1
        end
    end
    return 0
end

function GreedyOrder!(G; method="MF", minimize=false)
    n = nv(G)
    H = copy(G)
    order = zeros(UInt16, n)
    if minimize == true
        F = Vector{Array{UInt16,2}}(undef, n)
    end
    for i = 1:n
        if minimize == true
            F[i] = zeros(UInt16, 2, 1)
        end
        if method == "MF"
            num_fill = zeros(Int, n)
            for j = 1:n
                neib = neighbors(H, j)
                lneib = length(neib)
                if lneib > 0
                    sg = 0
                    for k = 1:lneib, l = k+1:lneib
                        if has_edge(H, neib[k], neib[l])
                            sg += 1
                        end
                    end
                    num_fill[j] = Int(lneib*(lneib-1)/2-sg)
                elseif order[j] != 0
                    num_fill[j] = typemax(Int)
                end
            end
            v = argmin(num_fill)
        else
            deg = LightGraphs.degree(H)
            for j = 1:n
                if order[j] != 0
                    deg[j] = typemax(Int)
                end
            end
            v = argmin(deg)
        end
        order[v] = i
        neib = copy(neighbors(H, v))
        for j = 1:length(neib)
            rem_edge!(H, v, neib[j])
            for k = j+1:length(neib)
                if add_edge!(H, neib[j], neib[k])
                    add_edge!(G, neib[j], neib[k])
                    if minimize == true
                        F[i] = [F[i] [neib[j]; neib[k]]]
                    end
                end
            end
        end
        if minimize == true
            F[i] = F[i][:, 2:end]
        end
    end
    if minimize == true
        G,order = MinimalChordal!(G, order, F)
    end
    return G,order
end

function chordal_cliques!(G; method="MF", minimize=false)
    G,order = GreedyOrder!(G, method=method, minimize=minimize)
    n = nv(G)
    candidate_cliques = Vector{Vector{UInt16}}(undef, n)
    for i = 1:n
        inter = intersect(neighbors(G, i), findall(x->order[x]>order[i], 1:n))
        inter = [inter;i]
        candidate_cliques[i] = sort!(inter)
    end
    sort!(candidate_cliques, by = x -> length(x))
    reverse!(candidate_cliques)
    maximal_cliques = [first(candidate_cliques)]
    for clique in Iterators.drop(candidate_cliques, 1)
        if all(other_clique -> !(clique ⊆ other_clique), maximal_cliques)
            push!(maximal_cliques, clique)
        end
    end
    cliquesize = length.(maximal_cliques)
    cql = length(cliquesize)
    return maximal_cliques,cql,cliquesize
end

function CheckChordal(G)
    n = nv(G)
    order = MCS(G)
    flag = true
    for i = 1:n
        v = findfirst(x->order[x]==i, 1:n)
        neib = intersect(neighbors(G, v), findall(x->order[x]>i, 1:n))
        for j = 1:length(neib)-1, k = j+1:length(neib)
            if !has_edge(G, neib[j], neib[k])
                flag = false
            end
        end
    end
    return flag,order
end

function add_clique!(G, nodes)
    for i in 1:length(nodes)-1, j in i+1:length(nodes)
        add_edge!(G, nodes[i], nodes[j])
    end
end

function max_cliques(G)
    cliques = convert(Vector{Vector{UInt16}}, maximal_cliques(G))
    sort!.(cliques)
    cliquesize = length.(cliques)
    cql = length(cliquesize)
    return cliques,cql,cliquesize
end

function fill!(G, order)
    n = nv(G)
    F = Vector{Array{UInt16,2}}(undef, n)
    for i = 1:n
        v = findfirst(x->order[x]==i, 1:n)
        neib = copy(intersect(neighbors(G, v), findall(x->order[x]>i, 1:n)))
        F[i] = zeros(UInt16, 2, 1)
        for j = 1:length(neib)-1, k = j+1:length(neib)
            if add_edge!(G, neib[j], neib[k])
                F[i] = [F[i] [neib[j]; neib[k]]]
            end
        end
        F[i] = F[i][:,2:end]
    end
    return G,F
end


text:
function remap(ln::AbstractString)
    return replace(
        ln,
        "Const " => "const ",
        "Private Function" => "function ",
        "Public Function " => "function ",
        "End Function" => "end",
        "If " => "if ",
        "End If" => "end",
        " Then" => "",
        "Else" => "else",
        "ElseIf " => "elseif ",
        "Do While" => "while",
        "Loop" => "end",
        "While " => "while ",
        "Wend" => "end",
        " Not " => " !",
        " And " => " && ",
        " Or " => " || ",
        "False" => "false",
        "True" => "true",
        "ByVal " => "",
        " As Boolean" => "::Bool",
        " As Double" => "::Float64",
        "Abs(" => "abs(",
        "Exp(" => "exp(",
        "Log(" => "log(",
        "Max(" => "max(",
        "Min(" => "min(",
        "Sqr(" => "abs2(",
        r"([0-9])#" => s"\1.0",
        "'" => "#",
    )
end

ext: jl


text:
abstract type AbstractEquationOfState{𝒯} end

struct BarotropicFluid{𝒯} <: AbstractEquationOfState{𝒯} end
struct IdealGas{𝒯} <: AbstractEquationOfState{𝒯} end

@inline function pressure(::BarotropicFluid{(:ρ, :ρu)}, state, aux, params)
    ρ = state.ρ
    cₛ = params.cₛ
    ρₒ = params.ρₒ

    return (cₛ * ρ)^2 / (2 * ρₒ)
end

@inline function pressure(::IdealGas{(:ρ, :ρu, :ρe)}, state, aux, params)
    ρ = state.ρ
    ρu = state.ρu
    ρe = state.ρe
    Φ = aux.Φ
    γ = params.γ

    return (γ - 1) * (ρe - dot(ρu, ρu) / 2ρ - ρ * Φ)
end

@inline function sound_speed(::BarotropicFluid{(:ρ, :ρu)}, state, aux, params)
    cₛ = params.cₛ
    ρₒ = params.ρₒ
    ρ = state.ρ

    return cₛ * sqrt(ρ / ρₒ)
end

@inline function sound_speed(eos::IdealGas{(:ρ, :ρu, :ρe)}, state, aux, params)
    ρ = state.ρ
    ρu = state.ρu
    ρe = state.ρe
    Φ = aux.Φ
    γ = params.γ

    p = calc_pressure(eos, state, aux, params)

    return sqrt(γ * p / ρ)
end

max_stars_repo_name: CliMA/CouplerMachine

max_issues_repo_name: CliMA/CouplerMachine


text:
module comparison

include("propwing.jl")
import CCBlade
using PyPlot
close("all")

Vinf = 50.0  # freestream speed
J = 0.85  # prop advance ratio
etaprop_center = [0.3, 0.6]  # spanwise location of prop (center)
cw = 1.0  # 1 if clockwise, -1 if counter-clockwise (when viewed from upstream)
xtowing = 0.2  # distance from prop to wing (for wake contraction) - assuming it is the same for all props
gamma = 2*pi/180  # flight path angle
alpha_wing_relative_to_fuselage = 2*pi/180 # angle of attack of wing relative to fuselage reference line
alpha = gamma + alpha_wing_relative_to_fuselage  # angle of attack
rho = 1.225  # air density  (irrelevant in normalization)

# --------- propeller --------
# geometry
Rhub = 0.035/2
Rtip = 0.236/2

# chord_prop = Rtip*[0.158788, 0.15083, 0.146158, 0.145793, 0.148978, 0.153175, 0.154064, 0.145548, 0.119747]
# twist_prop = pi/180*[54.9242, 49.6058, 44.6923, 40.1836, 36.0798, 32.3809, 29.0869, 26.1977, 23.7134]
# r = linspace(Rhub+0.05*Rtip, 0.95*Rtip, length(chord_prop))
# r = Rtip*[2.31e-01, 4.00e-01, 5.11e-01, 6.15e-01, 7.20e-01, 7.89e-01, 8.60e-01, 9.30e-01, 9.65e-01, 9.99e-01]
# chord_prop = Rtip*[1.56e-01, 1.46e-01, 1.47e-01, 1.52e-01, 1.54e-01, 1.52e-01, 1.38e-01, 1.10e-01, 9.13e-02, 6.87e-02]
# twist_prop = pi/180*[5.27e+01, 4.39e+01, 3.85e+01, 3.35e+01, 2.91e+01, 2.69e+01, 2.47e+01, 2.35e+01, 2.25e+01, 2.19e+01]
r = 1e-3*[30, 45, 60, 75, 90, 105]
chord_prop = 1e-3*[11.88, 15.59, 18.81, 19.55, 18.32, 3.96]
twist_prop = pi/180*[32.5, 26.5, 23.5, 19, 16.5, 14]

twist_75 = propwing.interp1(r, twist_prop, [0.75*Rtip])[1]
pitch = 25*pi/180 - twist_75
twist_prop += pitch + 6.6*pi/180

B = 4  # number of blades

aftype = CCBlade.af_from_aerodynfile("NACA64_A17.dat")
n = length(r)
af = Array{CCBlade.AirfoilData}(n)
for i = 1:n
    af[i] = aftype
end

Dprop = Rtip*2  # propeller diameter
Omega = Vinf/(J*Dprop)*2*pi  # rotation rate

uwake, vwake, T, Q, reff = propwing.propanalysis(Rhub, Rtip, r, chord_prop, twist_prop, B, af, Vinf, Omega, rho, xtowing)

nprop = length(etaprop_center)
println("T = ", T*nprop)
n = Omega/(2*pi)
CT = T / (rho * n^2 * Dprop^4)

# --------- wing geometry -------
AR = 5.33
sectionspan = [0.64]
twist_wing = [0, 0]*pi/180
b = sectionspan[1]*2
chord_wing = [b/AR, b/AR]
tc = [0.12, 0.12]  # irrelevant for the aerodynamics
sweep = [0]*pi/180
dihedral = [0]
N = [100]

alpha = 13.5*pi/180  # near stall
L, Di, CP, cl, cllocal, Vinfeff, alphaeff = propwing.winganalysis(sectionspan, chord_wing, twist_wing, tc, sweep, dihedral, N, alpha, rho, Vinf, reff, uwake, vwake, etaprop_center, cw)

semispan = sum(sectionspan)
figure()
plot(CP.y/semispan, cllocal)
plot([0, 1], [1.2, 1.2], "--")

figure()
plot(CP.y/semispan, alphaeff*180/pi)

q = 0.5*rho*Vinf^2
S = 2*sum(CP.chord.*CP.ds)
CLmax = L/(q*S)
println("Vinf = ", Vinf)
println("L = ", L)
println("Di = ", Di)

# -----  comparison case with 4 props --------
etaprop_center = [0.2, 0.4, 0.6, 0.8]  # spanwise location of prop (center)

# iterate on Vinf until I match lift.
Vinf = 48.7

Dprop = sqrt(2*Dprop^2*50^2/4/Vinf^2)  # I know what the diameter needs to be for constant thrust because J and solidity are constant.
# resize prop for new diameter
r *= (Dprop/2)/Rtip
Rhub *= (Dprop/2)/Rtip
chord_prop *= (Dprop/2)/Rtip
Rtip = Dprop/2
# compute new rotation rate (for same J)
Omega = Vinf/(J*Dprop)*2*pi

uwake, vwake, T, Q, reff = propwing.propanalysis(Rhub, Rtip, r, chord_prop, twist_prop, B, af, Vinf, Omega, rho, xtowing)

nprop = length(etaprop_center)
println("T = ", T*nprop)

alpha = 13.5*pi/180
L, Di, CP, cl, cllocal, Vinfeff, alphaeff = propwing.winganalysis(sectionspan, chord_wing, twist_wing, tc, sweep, dihedral, N, alpha, rho, Vinf, reff, uwake, vwake, etaprop_center, cw)

semispan = sum(sectionspan)
figure()
plot(CP.y/semispan, cllocal)
plot([0, 1], [1.2, 1.2], "--")

q = 0.5*rho*Vinf^2
S = 2*sum(CP.chord.*CP.ds)
CLmax = L/(q*S)
println("Vinf = ", Vinf)
println("L = ", L)
println("Di = ", Di)

end


text:
using Documenter, ControlSimTools

makedocs(;
    modules=[ControlSimTools],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
    ],
    repo="https://github.com/jonniedie/ControlSimTools.jl/blob/{commit}{path}#L{line}",
    sitename="ControlSimTools.jl",
    authors="Jonnie Diegelman",
    assets=String[],
)

deploydocs(;
    repo="github.com/jonniedie/ControlSimTools.jl",
)


text:
################################################################################
# This file is to defines commonly used constraints for power flow models
# This will hopefully make everything more compositional
################################################################################

"""
Checks that all cost models are present and of the same type
"""
function check_cost_models(pm::GenericPowerModel, nws)
    model = nothing

    for n in nws
        ref = pm.ref[:nw][n]
        for (i,gen) in ref[:gen]
            if haskey(gen, "cost")
                if model == nothing
                    model = gen["model"]
                else
                    if gen["model"] != model
                        error("cost models are inconsistent, the typical model is $(model) however model $(gen["model"]) is given on generator $(i)")
                    end
                end
            else
                error("no cost given for generator $(i)")
            end
        end
        for (i,dcline) in ref[:dcline]
            if haskey(dcline, "model")
                if model == nothing
                    model = dcline["model"]
                else
                    if dcline["model"] != model
                        error("cost models are inconsistent, the typical model is $(model) however model $(dcline["model"]) is given on dcline $(i)")
                    end
                end
            else
                error("no cost given for dcline $(i)")
            end
        end
    end

    return model
end

""
function objective_min_fuel_cost(pm::GenericPowerModel, nws=[pm.cnw])
    model = check_cost_models(pm, nws)

    if model == 1
        return objective_min_pwl_fuel_cost(pm, nws)
    elseif model == 2
        return objective_min_polynomial_fuel_cost(pm, nws)
    else
        error("Only cost models of types 1 and 2 are supported at this time, given cost model type of $(model)")
    end

end

"""
Checks that all cost models are polynomials, quadratic or less
"""
function check_polynomial_cost_models(pm::GenericPowerModel, nws)
    for n in nws
        ref = pm.ref[:nw][n]
        for (i,gen) in ref[:gen]
            @assert gen["model"] == 2
            if length(gen["cost"]) > 3
                error("only cost models of degree 3 or less are supported at this time, given cost model of degree $(length(gen["cost"])) on generator $(i)")
            end
        end
        for (i,dcline) in ref[:dcline]
            @assert dcline["model"] == 2
            if length(dcline["cost"]) > 3
                error("only cost models of degree 3 or less are supported at this time, given cost model of degree $(length(dcline["cost"])) on dcline $(i)")
            end
        end
    end
end

""
function objective_min_polynomial_fuel_cost(pm::GenericPowerModel, nws=[pm.cnw])
    check_polynomial_cost_models(pm, nws)

    pg = Dict(n => pm.var[:nw][n][:pg] for n in nws)
    dc_p = Dict(n => pm.var[:nw][n][:p_dc] for n in nws)

    from_idx = Dict()
    for n in nws
        ref = pm.ref[:nw][n]
        from_idx[n] = Dict(arc[1] => arc for arc in ref[:arcs_from_dc])
    end

    return @objective(pm.model, Min,
        sum(
            sum(gen["cost"][1]*pg[n][i]^2 + gen["cost"][2]*pg[n][i] + gen["cost"][3] for (i,gen) in pm.ref[:nw][n][:gen]) +
            sum(dcline["cost"][1]*dc_p[n][from_idx[n][i]]^2 + dcline["cost"][2]*dc_p[n][from_idx[n][i]] + dcline["cost"][3] for (i,dcline) in pm.ref[:nw][n][:dcline])
        for n in nws)
    )
end

""
function objective_min_polynomial_fuel_cost(pm::GenericPowerModel{T}, nws=[pm.cnw]) where T <: AbstractConicPowerFormulation
    check_polynomial_cost_models(pm, nws)

    pg = Dict(n => pm.var[:nw][n][:pg] for n in nws)
    dc_p = Dict(n => pm.var[:nw][n][:p_dc] for n in nws)

    from_idx = Dict()
    for n in nws
        ref = pm.ref[:nw][n]
        from_idx[n] = Dict(arc[1] => arc for arc in ref[:arcs_from_dc])
    end

    pg_sqr = Dict()
    dc_p_sqr = Dict()
    for n in nws
        pg_sqr[n] = pm.var[:nw][n][:pg_sqr] = @variable(pm.model,
            [i in keys(pm.ref[:nw][n][:gen])], basename="$(n)_pg_sqr",
            lowerbound = pm.ref[:nw][n][:gen][i]["pmin"]^2,
            upperbound = pm.ref[:nw][n][:gen][i]["pmax"]^2
        )
        for (i, gen) in pm.ref[:nw][n][:gen]
            @constraint(pm.model, norm([2*pg[n][i], pg_sqr[n][i]-1]) <= pg_sqr[n][i]+1)
        end

        dc_p_sqr[n] = pm.var[:nw][n][:dc_p_sqr] = @variable(pm.model,
            [i in keys(pm.ref[:nw][n][:dcline])], basename="$(n)_dc_p_sqr",
            lowerbound = pm.ref[:nw][n][:dcline][i]["pminf"]^2,
            upperbound = pm.ref[:nw][n][:dcline][i]["pmaxf"]^2
        )

        for (i, dcline) in pm.ref[:nw][n][:dcline]
            @constraint(pm.model, norm([2*dc_p[n][from_idx[n][i]], dc_p_sqr[n][i]-1]) <= dc_p_sqr[n][i]+1)
        end
    end

    return @objective(pm.model, Min,
        sum(
            sum( gen["cost"][1]*pg_sqr[n][i] + gen["cost"][2]*pg[n][i] + gen["cost"][3] for (i,gen) in pm.ref[:nw][n][:gen]) +
            sum(dcline["cost"][1]*dc_p_sqr[n][i]^2 + dcline["cost"][2]*dc_p[n][from_idx[n][i]] + dcline["cost"][3] for (i,dcline) in pm.ref[:nw][n][:dcline])
        for n in nws)
    )
end

"""
compute m and b from points pwl points
"""
function slope_intercepts(points::Array{T,1}) where T <: Real
    line_data = []

    for i in 3:2:length(points)
        x1 = points[i-2]
        y1 = points[i-1]
        x2 = points[i-0]
        y2 = points[i+1]

        m = (y2 - y1)/(x2 - x1)
        b = y1 - m * x1

        line = Dict{String,Any}(
            "slope" => m,
            "intercept" => b
        )
        push!(line_data, line)
    end

    return line_data
end

"""
compute lines in m and b from from pwl cost models
data is a list of components
"""
function get_lines(data)
    lines = Dict{Int,Any}()
    for (i,comp) in data
        @assert comp["model"] == 1
        line_data = slope_intercepts(comp["cost"])
        lines[i] = line_data
        for i in 2:length(line_data)
            if line_data[i-1]["slope"] > line_data[i]["slope"]
                error("non-convex pwl function found in points $(comp["cost"])\nlines: $(line_data)")
            end
        end
    end
    return lines
end

""
function objective_min_pwl_fuel_cost(pm::GenericPowerModel, nws=[pm.cnw])
    #check_polynomial_cost_models(pm, nws)

    pg = Dict(n => pm.var[:nw][n][:pg] for n in nws)
    dc_p = Dict(n => pm.var[:nw][n][:p_dc] for n in nws)

    pg_cost = Dict()
    dc_p_cost = Dict()
    for n in nws
        pg_cost[n] = pm.var[:nw][n][:pg_cost] = @variable(pm.model,
            [i in keys(pm.ref[:nw][n][:gen])], basename="$(n)_pg_cost"
        )

        # pwl cost
        gen_lines = get_lines(pm.ref[:nw][n][:gen])
        for (i, gen) in pm.ref[:nw][n][:gen]
            for line in gen_lines[i]
                @constraint(pm.model, pg_cost[n][i] >= line["slope"]*pg[n][i] + line["intercept"])
            end
        end

        dc_p_cost[n] = pm.var[:nw][n][:dc_p_cost] = @variable(pm.model,
            [i in keys(pm.ref[:nw][n][:dcline])], basename="$(n)_dc_p_cost",
        )

        # pwl cost
        dcline_lines = get_lines(pm.ref[:nw][n][:dcline])
        for (i, dcline) in pm.ref[:nw][n][:dcline]
            for line in dcline_lines[i]
                @constraint(pm.model, dc_p_cost[n][i] >= line["slope"]*dc_p[n][i] + line["intercept"])
            end
        end

        #for (i, dcline) in pm.ref[:nw][n][:dcline]
        #    @constraint(pm.model, norm([2*dc_p[n][from_idx[n][i]], dc_p_sqr[n][i]-1]) <= dc_p_sqr[n][i]+1)
        #end
    end

    #from_idx = Dict()
    #for n in nws
    #    ref = pm.ref[:nw][n]
    #    from_idx[n] = Dict(arc[1] => arc for arc in ref[:arcs_from_dc])
    #end

    return @objective(pm.model, Min,
        sum(
            sum( pg_cost[n][i] for (i,gen) in pm.ref[:nw][n][:gen]) +
            sum( dc_p_cost[n][i] for (i,dcline) in pm.ref[:nw][n][:dcline])
        for n in nws)
    )
end

"Cost of building branches"
function objective_tnep_cost(pm::GenericPowerModel, nws=[pm.cnw])
    branch_ne = Dict(n => pm.var[:nw][n][:branch_ne] for n in nws)
    branches = Dict(n => pm.ref[:nw][n][:ne_branch] for n in nws)

    return @objective(pm.model, Min,
        sum(
            sum( branch["construction_cost"]*branch_ne[n][i] for (i,branch) in branches[n])
        for n in nws)
    )
end


text:
module CthulhuTestSandbox

function testf_revise()
    T = rand() > 0.5 ? Int64 : Float64
    sum(rand(T, 100))
end

export testf_revise

end


text:
module Log

using Base.Test
using Bitmanip, Defs, FloatProps, Rounding

export vlog, vlog2, vlog10

# Note: The Coefficients were determined by a Mathematica script,
# available at <https://bitbucket.org/eschnett/vecmathlib>. They
# minimize the error in a least-square sense over the interval dx ∈
# [1/sqrt(2); sqrt(2)]. A better alternative would be to minimize the
# maximum of the error.

# One can use a smaller set of coefficients to increase speed at the
# extent of a slightly larger error.

# TODO: Use @evalpoly for this
@inline function vlog2_poly(y2::Float32)
    # Coefficients for Float32,
    # error=7.09807175879142775648452461821e-8
    r = mkFloat(Float32, 0.59723611417135718739797302426)
    r = muladd(r, y2, mkFloat(Float32, 0.961524413175528426101613434))
    r = muladd(r, y2, mkFloat(Float32, 2.88539097665498228703236701))
    return r
end

@inline function vlog2_poly(y2::Float64)
    # Coefficients for Float64,
    # error=1.48294180185938512675770096324e-16
    r = mkFloat(Float64, 0.243683403415639178527756320773)
    r = muladd(r, y2, mkFloat(Float64, 0.26136626803870009948502658))
    r = muladd(r, y2, mkFloat(Float64, 0.320619429891299265439389))
    r = muladd(r, y2, mkFloat(Float64, 0.4121983452028499242926))
    r = muladd(r, y2, mkFloat(Float64, 0.577078017761894161436))
    r = muladd(r, y2, mkFloat(Float64, 0.96179669392233355927))
    r = muladd(r, y2, mkFloat(Float64, 2.8853900817779295236))
    return r
end

const M_SQRT1_2 = 7.071067811865475244008443621048490392848359376884740365883398689953662392310596e-01
const M_SQRT2 = 1.414213562373095048801688724209698078569671875376948073176679737990732478462102e+00
@inline function vlog2{T}(x::T)
    x0 = x
    # Rescale

    e0 = asInt(T,x) & exponent_mask(T) >> mantissa_bits(T)
    e = e0 - exponent_offset(T)
    m = asInt(T,x) & mantissa_mask(T)
    # TODO: handle subnormals

    x = asFloat(T, m | (asInt(T, mkFloat(T,1)) & exponent_mask(T)))
    b = x > M_SQRT2
    x = vifthen(b, mkFloat(T,0.5)*x, x)
    e = vifthen(b, e+mkInt(T,1), e)
    @vassert mkFloat(T,M_SQRT1_2) <= x <= mkFloat(T,M_SQRT2)

    y = (x - mkFloat(T,1)) / (x + mkFloat(T,1))
    y2 = y*y
    r = vlog2_poly(y2)
    r *= y

    # Undo rescaling
    # TODO: convert to float by
    # - re-interpreting as subnormal, then multiplying with large number
    # - or'ing an exponent, then subtracting the hidden mantissa bit
    r += mkFloat(T, e)

    r = vifthen(x0 < mkFloat(T,0) || visnan(x0), mkFloat(T,NaN), r)
    r = vifthen(x0 == mkFloat(T,0), -mkFloat(T,Inf), r)
    r = vifthen(x0 == mkFloat(T,Inf), mkFloat(T,Inf), r)

    return r
end

const M_LN2 = 6.931471805599453094172321214581765680755001343602552541206800094933936219696955e-01
@inline function vexp{T}(x::T)
    # TODO: Re-calculate the coefficients for this function to
    # increase accuracy and avoid the multiplication
    return vlog2(x) * mkFloat(T,M_LN2)
end

const M_LN10 = 2.302585092994045684017991454684364207601101488628772976033327900967572609677367e+00
@inline function vlog10{T}(x::T)
    # TODO: Re-calculate the coefficients for this function to
    # increase accuracy and avoid the multiplication
    return vlog2(x) * mkFloat(T,M_LN10)
end

function selftest_type(T::Type)
    zero = convert(T,0)
    min = realmin(T)
    eps = Base.eps(T)
    max = realmax(T)
    vals = T[1.0,
             zero, #TODO nextfloat(zero),
             #TODO prevfloat(min), min, nextfloat(min),
             0.1*eps, eps, 10.0*eps,
             0.1*sqrt(eps), sqrt(eps), 10.0*sqrt(eps),
             1.0-eps, 1.0, 1.0+eps,
             2.0, 10.0, 100.0,
             0.1*sqrt(max), sqrt(max), 10.0*sqrt(max),
             0.1*max, prevfloat(max), max, Inf, NaN]
    for uval in vals
        for val in (uval, -uval)
            local log2_val::T
            try
                log2_val = log2(val)
            catch
                log2_val = convert(T,NaN)
            end
            vlog2_val = vlog2(val)
            @test (isequal(vlog2_val, log2_val) ||
                   isapprox(vlog2_val, log2_val, rtol=eps, atol=eps))
        end
    end
end

function selftest()
    selftest_type(Float32)
    selftest_type(Float64)
end

selftest()

end


text:
export BostonHousing

"""
Boston Housing Dataset.

Sources:
   (a) Origin:  This dataset was taken from the StatLib library which is
                maintained at Carnegie Mellon University.
   (b) Creator:  Harrison, D. and Rubinfeld, D.L. 'Hedonic prices and the
                 demand for clean air', J. Environ. Economics & Management,
                 vol.5, 81-102, 1978.
   (c) Date: July 7, 1993

Number of Instances: 506

Number of Attributes: 13 continuous attributes (including target
                            attribute "MEDV"), 1 binary-valued attribute.

Features:

    1. CRIM      per capita crime rate by town
    2. ZN        proportion of residential land zoned for lots over 25,000 sq.ft.
    3. INDUS     proportion of non-retail business acres per town
    4. CHAS      Charles River dummy variable (= 1 if tract bounds river; 0 otherwise)
    5. NOX       nitric oxides concentration (parts per 10 million)
    6. RM        average number of rooms per dwelling
    7. AGE       proportion of owner-occupied units built prior to 1940
    8. DIS       weighted distances to five Boston employment centres
    9. RAD       index of accessibility to radial highways
    10. TAX      full-value property-tax rate per 10,000 dollars
    11. PTRATIO  pupil-teacher ratio by town
    12. B        1000(Bk - 0.63)^2 where Bk is the proportion of blacks by town
    13. LSTAT    % lower status of the population

Target:

    14. MEDV     Median value of owner-occupied homes in 1000's of dollars

Note: Variable #14 seems to be censored at 50.00 (corresponding to a median price of \\\$50,000);
Censoring is suggested by the fact that the highest median price of exactly \\\$50,000 is reported in 16 cases,
while 15 cases have prices between \\\$40,000 and \\\$50,000, with prices rounded to the nearest hundred.
Harrison and Rubinfeld do not mention any censoring.

The data file stored in this repo is a copy of the This is a copy of UCI ML housing dataset.
https://archive.ics.uci.edu/ml/machine-learning-databases/housing/

# Interface

- [`BostonHousing.features`](@ref)
- [`BostonHousing.targets`](@ref)
- [`BostonHousing.feature_names`](@ref)
"""
module BostonHousing

using DataDeps
using DelimitedFiles

export features, targets, feature_names

const DATA = joinpath(@__DIR__, "boston_housing.csv")

"""
    targets(; dir = nothing)

Get the targets for the Boston housing dataset,
a 506 element array listing the targets for each example.

```jldoctest
julia> using MLDatasets: BostonHousing

julia> target = BostonHousing.targets();

julia> summary(target)
"1×506 Matrix{Float64}"

julia> target[1]
24.0
```
"""
function targets(; dir = nothing)
    housing = readdlm(DATA, ',')
    reshape(Vector{Float64}(housing[2:end,end]), (1, 506))
end

"""
    feature_names()

Return the  the names of the features provided in the dataset.
"""
function feature_names()
    ["crim","zn","indus","chas","nox","rm","age","dis","rad","tax","ptratio","b","lstat"]
end

"""
    features()

Return the features of the Boston Housing dataset. This is a 13x506 Matrix of Float64 datatypes.
The values are in the order ["crim","zn","indus","chas","nox","rm","age","dis","rad","tax","ptratio","b","lstat"].
It has 506 examples.

```jldoctest
julia> using MLDatasets: BostonHousing

julia> features = BostonHousing.features();

julia> summary(features)
"13×506 Matrix{Float64}"
```
"""
function features()
    housing = readdlm(DATA, ',')
    Matrix{Float64}(housing[2:end,1:13])' |> collect
end

end # module


text:
#!/usr/bin/env julia

using Memoize
using LRUCache

const DATA_DIR   = "../19_data"
const DQUOTE = "\""
const INPUT_FILE = "input.txt"

# ----------------------------------------
macro assert(bool_expr)
    message = string("Assertion: ", bool_expr, " failed")
    :($(esc(bool_expr)) || error($message))
end
# ----------------------------------------

abstract type Rule end

struct RuleRefs <: Rule
    branches::Array{Array{Int,1},1}
end

struct RuleSimple <: Rule
    str::String
end

# ----------------------------------------

function get_valid_count(rules::Dict{Int,Rule}, rid::Int, msgs::Array{String,1})::Int
    rid2msgs = rid->get_valid_msgs(rules, rid)

    @memoize LRU{Tuple{Any, Any}, Any}(maxsize=10000) #=
 =# function get_valid_msgs(rules::Dict{Int,Rule}, rid::Int)::Set{String}
        result = nothing
        rule = rules[rid]
        if typeof(rule) == RuleSimple
            return Set([string(rule.str)])
        else  # RuleRefs
            branch2msgs = branch->msg_prod(map(rid2msgs, branch))
            return msg_union(map(branch2msgs, rule.branches))
        end
    end

    valid_msgs = get_valid_msgs(rules::Dict{Int,Rule}, rid::Int)
    return length([1 for msg in msgs if msg in valid_msgs])
end

function map2(f, arr2)
    return [[f(x) for x in coll] for coll in arr2]
end

function msg_prod(sets::Array{Set{String},1})::Set{String}
    prod = (set1, set2) -> Set([a*b for a in set1 for b in set2])
    return foldl(prod, sets)
end

function msg_union(sets::Array{Set{String},1})::Set{String}
    return reduce(union, sets; init=Set())
end

function rule_parse(rule_str::String)::Rule
    # println("Parsing $(rule_str)")
    if occursin(DQUOTE, rule_str)
	m = match(r"\"(\w+)\"", rule_str)
        if isnothing(m)
            error("Unrecognized rule syntax (simple)")
        else
            str = string(m.captures[1])
            return RuleSimple(str)
        end
    else
        m = match(r"((?:\s|\d|\|)+)", rule_str)
        if isnothing(m)
            error("Unrecognized rule syntax (subrule)")
        else
            rule_refs = m.captures[1]
            strings2 = map(split_nonempty, split(rule_refs, "|"))
            ints2 = map2(s->parse(Int, s), strings2)
            return RuleRefs(ints2)
        end
    end
end

function rules_parse(rule_strs::Array{String,1})::Dict{Int,Rule}
    rules::Dict{Int,Rule} = Dict{Int,Rule}()
    for rule_str in rule_strs
        (rid_str, rule_body) = split(rule_str, ": ")
        rid::Int = parse(Int, rid_str)
        rules[rid] = rule_parse(string(rule_body))
    end
    return rules
end

function rules_print(rules::Dict{Int,Rule})
    println("===== Rules =====")
    for rid in keys(rules)
        println("$(rid): $(rules[rid])")
    end
end

split_nonempty = s->filter(z->!isempty(z), split(s, " "))

function split_on_empty(xs)
    return split_on_pred(xs, x->isempty(x))
end

function split_on_pred(xs, pred)
    result = []
    current = []
    for x in xs
       if pred(x) && !isempty(current)
           push!(result, current)
           current = []
       else
           push!(current, x)
       end
    end
    if !isempty(current)
        push!(result, current)
    end
    return result
end

# ----------------------------------------

function test_all()
    test_msg_prod()
end

function test_msg_prod()
    args = [Set(["a","b"]), Set(["x","y"])]
    result = msg_prod(args)
    @assert length(result) == 4
end

# ----------------------------------------

function main(filename::String)
    test_all()

    lines::Array{String,1} = readlines(filename)
    sections = split_on_empty(lines)
    if length(sections) != 2
        error("Input does not contain rules and messages sections")
    end
    rule_lines::Array{String,1}, msgs::Array{String,1} = tuple(sections...)
    println("# Rule count: $(length(rule_lines))")
    println("# Msgs count: $(length(msgs))")

    rules::Dict{Int,Rule} = rules_parse(rule_lines)
    rid = 0
    valid_count = get_valid_count(rules, rid, msgs)
    println("valid_count=$(valid_count)")
end

main("$DATA_DIR/$INPUT_FILE")

max_stars_repo_name: jaycoskey/AdventOfCode

max_issues_repo_name: jaycoskey/AdventOfCode

max_forks_repo_name: jaycoskey/AdventOfCode


text:
size(a::MultiValue) = size(a.array)

length(a::MultiValue) = length(a.array)

@propagate_inbounds function getindex(
    a::MultiValue{S,T,N}, I::Vararg{Integer,N}) where {S,T,N}
    a.array[I...]
end

@propagate_inbounds function getindex(a::MultiValue, i::Integer)
    a.array[i]
end

eltype(a::Type{MultiValue{S,T,N,L}}) where {S,T,N,L} = T

@inline iterate(a::MultiValue) = iterate(a.array)

@inline iterate(a::MultiValue, state) = iterate(a.array, state)

eachindex(a::MultiValue) = eachindex(a.array)

function CartesianIndices(a::MultiValue)
  CartesianIndices(a.array)
end

function LinearIndices(a::MultiValue)
  LinearIndices(a.array)
end


text:
module Qualification
export deprecate_qualifier_macro

@nospecialize
function deprecate_qualifier_macro(o, ln::LineNumberNode)
    s = string(o)
    trunc = min(length(s), 20)
    s = SubString(s, 1:trunc)
    @warn (
        "Deprecated use of qualifier $s at $(ln.file):$(ln.line):\n" *
        "Scoping specifiers such as 'internal', 'public' are deprecated. " *
        "Now the scope of a pattern is consistent with the visibility of the pattern object in current module."
    )
end
@specialize

end


text:
module TurbulenceClosures

export
    AbstractEddyViscosityClosure,
    VerticalScalarDiffusivity,
    HorizontalScalarDiffusivity,
    ScalarDiffusivity,
    VerticalScalarBiharmonicDiffusivity,
    HorizontalScalarBiharmonicDiffusivity,
    ScalarBiharmonicDiffusivity,
    TwoDimensionalLeith,
    SmagorinskyLilly,
    AnisotropicMinimumDissipation,
    ConvectiveAdjustmentVerticalDiffusivity,
    IsopycnalSkewSymmetricDiffusivity,

    ExplicitTimeDiscretization,
    VerticallyImplicitTimeDiscretization,

    DiffusivityFields,
    calculate_diffusivities!,

    ∇_dot_qᶜ,
    ∂ⱼ_τ₁ⱼ,
    ∂ⱼ_τ₂ⱼ,
    ∂ⱼ_τ₃ⱼ,

    cell_diffusion_timescale

using CUDA
using KernelAbstractions

import Oceananigans.Utils: with_tracers

using Oceananigans
using Oceananigans.Architectures
using Oceananigans.Grids
using Oceananigans.Operators
using Oceananigans.BoundaryConditions
using Oceananigans.Fields
using Oceananigans.BuoyancyModels
using Oceananigans.Utils

using Oceananigans.Architectures: AbstractArchitecture, device

const VerticallyBoundedGrid{FT} = AbstractGrid{FT, <:Any, <:Any, <:Bounded}

#####
##### Abstract types
#####

"""
    AbstractTurbulenceClosure

Abstract supertype for turbulence closures.
"""
abstract type AbstractTurbulenceClosure{TimeDiscretization} end

@inline get_closure_i(i, closure::AbstractVector{<:AbstractTurbulenceClosure}) = @inbounds closure[i]
@inline get_closure_i(i, closure::AbstractTurbulenceClosure) = closure

@inline get_closure_ij(i, j, closure::AbstractMatrix{<:AbstractTurbulenceClosure}) = @inbounds closure[i, j]
@inline get_closure_ij(i, j, closure::AbstractTurbulenceClosure) = closure

# Fallbacks
validate_closure(closure) = closure
closure_summary(closure) = summary(closure)

const ClosureKinda = Union{Nothing, AbstractTurbulenceClosure, AbstractArray{<:AbstractTurbulenceClosure}}
add_closure_specific_boundary_conditions(closure::ClosureKinda, bcs, args...) = bcs

# To allow indexing a diffusivity with (Lx, Ly, Lz, i, j, k, grid)
struct DiscreteDiffusionFunction{F} <: Function
    func :: F
end

#####
##### Include module code
#####

include("implicit_explicit_time_discretization.jl")
include("turbulence_closure_utils.jl")
include("diffusion_operators.jl")
include("viscous_dissipation_operators.jl")
include("velocity_tracer_gradients.jl")
include("abstract_scalar_diffusivity_closure.jl")
include("abstract_eddy_viscosity_closure.jl")
include("abstract_scalar_biharmonic_diffusivity_closure.jl")
include("closure_tuples.jl")
include("isopycnal_rotation_tensor_components.jl")

# Implementations:
include("turbulence_closure_implementations/nothing_closure.jl")
include("turbulence_closure_implementations/scalar_diffusivity.jl")
include("turbulence_closure_implementations/scalar_biharmonic_diffusivity.jl")
include("turbulence_closure_implementations/leith_enstrophy_diffusivity.jl")
include("turbulence_closure_implementations/isopycnal_skew_symmetric_diffusivity.jl")
include("turbulence_closure_implementations/smagorinsky_lilly.jl")
include("turbulence_closure_implementations/anisotropic_minimum_dissipation.jl")
include("turbulence_closure_implementations/CATKEVerticalDiffusivities/CATKEVerticalDiffusivities.jl")
include("turbulence_closure_implementations/convective_adjustment_vertical_diffusivity.jl")

using .CATKEVerticalDiffusivities: CATKEVerticalDiffusivity

# Miscellaneous utilities
include("diffusivity_fields.jl")
include("turbulence_closure_diagnostics.jl")
include("vertically_implicit_diffusion_solver.jl")

#####
##### Some value judgements here
#####

end # module


text:
# Hello World! in Julia
# Navendu Pottekkat | 2020 | Apache-2.0 License
# https://www.apache.org/licenses/LICENSE-2.0

println("Hello World!")

max_stars_repo_path: Julia/Hello World/hello_world.jl

max_issues_repo_path: Julia/Hello World/hello_world.jl

max_forks_repo_path: Julia/Hello World/hello_world.jl


text:
using GeneralizedGenerated
using Random: GLOBAL_RNG
using TupleVectors: chainvec

export rand
EmptyNTtype = NamedTuple{(),Tuple{}} where T<:Tuple

function Base.rand(rng::AbstractRNG, d::ModelClosure, N::Int)
    r = chainvec(rand(rng, d), N)
    for j in 2:N
        @inbounds r[j] = rand(rng, d)
    end
    return r
end

Base.rand(d::ModelClosure, N::Int) = rand(GLOBAL_RNG, d, N)

@inline function Base.rand(m::ModelClosure; kwargs...)
    rand(GLOBAL_RNG, m; kwargs...)
end

@inline function Base.rand(rng::AbstractRNG, m::ModelClosure; cfg = NamedTuple(), ctx=NamedTuple())
    cfg′ = merge(cfg, (rng=rng,))
    gg_call(m, rand, cfg′, ctx, KeepReturn())
end

###############################################################################
# ctx::NamedTuple

@inline function tilde(::typeof(Base.rand), lens::typeof(identity), xname, x, d, cfg, ctx::NamedTuple)
    xnew = rand(cfg.rng, d)
    ctx′ = merge(ctx, NamedTuple{(dynamic(xname),)}((xnew,)))
    (xnew, ctx′, ctx′)
end

@inline function tilde(::typeof(Base.rand), lens, xname, x, d, cfg, ctx::NamedTuple)
    xnew = set(x, Lens!!(lens), rand(cfg.rng, d))
    ctx′ = merge(ctx, NamedTuple{(dynamic(xname),)}((xnew,)))
    (xnew, ctx′, ctx′)
end

###############################################################################
# ctx::Dict

@inline function tilde(::typeof(Base.rand), lens::typeof(identity), xname, x, d, cfg, ctx::Dict)
    x = rand(cfg.rng, d)
    ctx[dynamic(xname)] = x
    (x, ctx, ctx)
end

@inline function tilde(::typeof(Base.rand), lens, xname, x, m::AbstractConditionalModel, cfg, ctx::Dict)
    args = get(cfg.args, dynamic(xname), Dict())
    cfg = merge(cfg, (args = args,))
    tilde(rand, lens, xname, x, m(cfg.args), cfg, ctx)
end

###############################################################################
# ctx::Tuple{}

@inline function tilde(::typeof(Base.rand), lens::typeof(identity), xname, x, d, cfg, ctx::Tuple{})
    xnew = rand(cfg.rng, d)
    (xnew, (), xnew)
end

@inline function tilde(::typeof(Base.rand), lens, xname, x, d, cfg, ctx::Tuple{})
    xnew = set(x, Lens!!(lens), rand(cfg.rng, d))
    (xnew, (), xnew)
end

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


text:
export ft_read!

function ft_read!(ft_handle::Culong,
                 buffer::Array{UInt8,1},
                 bytestoread::Integer = length(buffer))
  bytesreturned = Ref{Cuint}()
  ft_status = ccall((:FT_Read, d2xx),
                      Cuint,
                      (Culong, Ptr{UInt8}, Cuint, Ref{Cuint}),
                      ft_handle, buffer, bytestoread, bytesreturned)
  checkstatus(ft_status)
  return (bytesreturned[])
end


text:
module TableFormat
    using JSON3
    export Format, Align, Group, Padding, Prefix, Scheme, Sign, DSymbol, Trim

    struct NamedValue{Name, T}
        value::T
        NamedValue{Name}(value::T) where {Name, Keys, T} = new{Name, T}(value)
    end

    struct TupleWithNamedValues{Name, Keys}
        values::NamedTuple
        TupleWithNamedValues{Name, Keys}(values) where {Name, Keys} = new{Name, Keys}(values)
    end

    Base.getproperty(tp::TupleWithNamedValues, prop::Symbol) = getfield(tp, :values)[prop]

    function tuple_with_named_values(name::Symbol, t::NamedTuple{Names}) where {Names}
        return TupleWithNamedValues{name, Names}((;zip(Names, NamedValue{name}.(values(t)))...))
    end

    function possible_values(::TupleWithNamedValues{Name, Keys}) where {Name, Keys}
        return join(string.(string(Name), ".", string.(Keys)), ", ")
    end

    const Align = tuple_with_named_values(:Align, (
        default = "",
        left = "<",
        right = ">",
        center = "^",
        right_sign = "="
    ))

    const Group = tuple_with_named_values(:Group, (
        no = "",
        yes = ","
    ))

    const Padding = tuple_with_named_values(:Padding, (no = "", yes =  0))

    const Prefix = tuple_with_named_values(:Prefix, (
        yocto = 10 ^ -24,
        zepto = 10 ^ -21,
        atto = 10 ^ -18,
        femto = 10 ^ -15,
        pico = 10 ^ -12,
        nano = 10 ^ -9,
        micro = 10 ^ -6,
        milli = 10 ^ -3,
        none = nothing,
        kilo = 10 ^ 3,
        mega = 10 ^ 6,
        giga = 10 ^ 9,
        tera = 10 ^ 12,
        peta = 10 ^ 15,
        exa = 10 ^ 18,
        zetta = 10 ^ 21,
        yotta = 10 ^ 24
    ))
    const Scheme = tuple_with_named_values(:Scheme, (
            default = "",
            decimal = "r",
            decimal_integer = "d",
            decimal_or_exponent = "g",
            decimal_si_prefix = "s",
            exponent = "e",
            fixed = "f",
            percentage = "%",
            percentage_rounded = "p",
            binary = "b",
            octal = "o",
            lower_case_hex = "x",
            upper_case_hex = "X",
            unicode = "c",
    ))
    const Sign = tuple_with_named_values(:Sign, (
        default = "",
        negative = "-",
        positive = "+",
        parantheses = "(",
        space = " "
    ))
    const DSymbol = tuple_with_named_values(:DSymbol, (
        no = "",
        yes = "\$",
        binary = "#b",
        octal = "#o",
        hex = "#x"
    ))
    const Trim = tuple_with_named_values(:Trim, (
        no = "",
        yes = "~"
    ))

    mutable struct Format
        locale
        nully
        prefix
        specifier
        function Format(;
                align = Align.default,
                fill = nothing,
                group = Group.no,
                padding = Padding.no,
                padding_width = nothing,
                precision = nothing,
                scheme = Scheme.default,
                sign = Sign.default,
                symbol = DSymbol.no,
                trim = Trim.no,

                symbol_prefix = nothing,
                symbol_suffix = nothing,
                decimal_delimiter = nothing,
                group_delimiter = nothing,
                groups = nothing,

                nully = "",

                si_prefix = Prefix.none
            )
            result = new(
                Dict(),
                "",
                Prefix.none.value,
                Dict{Symbol, Any}()
            )
            align!(result, align)
            fill!(result, fill)
            group!(result, group)
            padding!(result, padding)
            padding_width!(result, padding_width)
            precision!(result, precision)
            scheme!(result, scheme)
            sign!(result, sign)
            symbol!(result, symbol)
            trim!(result, trim)

            !isnothing(symbol_prefix) && symbol_prefix!(result, symbol_prefix)
            !isnothing(symbol_suffix) && symbol_suffix!(result, symbol_suffix)
            !isnothing(decimal_delimiter) && decimal_delimiter!(result, decimal_delimiter)
            !isnothing(group_delimiter) && group_delimiter!(result, group_delimiter)
            !isnothing(groups) && groups!(result, groups)
            !isnothing(nully) && nully!(result, nully)
            !isnothing(si_prefix) && si_prefix!(result, si_prefix)

            return result
        end
    end

    function check_named_value(t::TupleWithNamedValues{Name, Keys}, v::NamedValue{ValueName}) where {Name, ValueName, Keys}
        Name != ValueName && throw(ArgumentError("expected value to be one of $(possible_values(t))"))
        return true
    end
    function check_named_value(t::TupleWithNamedValues{Name, Keys}, v) where {Name, Keys}
        throw(ArgumentError("expected value to be one of $(possible_values(t))"))
    end

    check_char_value(v::Char) = string(v)
    function check_char_value(v::String)
        length(v) > 1 && throw(ArgumentError("expected char or string of length one"))
        return v
    end
    function check_char_value(v)
        throw(ArgumentError("expected char or string of length one"))
    end

    function check_int_value(v::Integer)
        v < 0 && throw(ArgumentError("expected value to be non-negative"))
    end

    function check_int_value(v)
        throw(ArgumentError("expected value to be Integer"))
    end

    function align!(f::Format, value)
        check_named_value(Align, value)
        f.specifier[:align] = value.value
    end

    function fill!(f::Format, value)
        if isnothing(value)
            f.specifier[:fill] = ""
            return
        end
        v = check_char_value(value)
        f.specifier[:fill] = v
    end

    function group!(f::Format, value)
        if value isa Bool
            value = value ? Group.yes : Group.no
        end
        check_named_value(Group, value)
        f.specifier[:group] = value.value
    end

    function padding!(f::Format, value)
        if value isa Bool
            value = value ? Padding.yes : Padding.no
        end
        check_named_value(Padding, value)
        f.specifier[:padding] = value.value
    end

    function padding_width!(f::Format, value)
        if isnothing(value)
            f.specifier[:width] = ""
        else
            check_int_value(value)
            f.specifier[:width] = value
        end
    end

    function precision!(f::Format, value)
        if isnothing(value)
            f.specifier[:precision] = ""
        else
            check_int_value(value)
            f.specifier[:precision] = ".$value"
        end
    end

    function scheme!(f::Format, value)
        check_named_value(Scheme, value)
        f.specifier[:type] = value.value
    end

    function sign!(f::Format, value)
        check_named_value(Sign, value)
        f.specifier[:sign] = value.value
    end

    function symbol!(f::Format, value)
        check_named_value(DSymbol, value)
        f.specifier[:symbol] = value.value
    end

    function trim!(f::Format, value)
        if value isa Bool
            value = value ? Trim.yes : Trim.no
        end
        check_named_value(Trim, value)
        f.specifier[:trim] = value.value
    end

    # Locale
    function symbol_prefix!(f::Format, value::AbstractString)
        if !haskey(f.locale, :symbol)
            f.locale[:symbol] = [value, ""]
        else
            f.locale[:symbol][1] = value
        end
    end

    function symbol_suffix!(f::Format, value::AbstractString)
        if !haskey(f.locale, :symbol)
            f.locale[:symbol] = ["", value]
        else
            f.locale[:symbol][2] = value
        end
    end

    function decimal_delimiter!(f::Format, value)
        v = check_char_value(value)
        f.locale[:decimal] = v
    end

    function group_delimiter!(f::Format, value)
        v = check_char_value(value)
        f.locale[:group] = v
    end

    function groups!(f::Format, value::Union{Vector{<:Integer}, <:Integer})
        groups = value isa Integer ? [value] : value
        isempty(groups) && throw(ArgumentError("groups cannot be empty"))

        for g in groups
            g < 0 && throw(ArgumentError("group entry must be non-negative integer"))
        end
        f.locale[:grouping] = groups
    end

    # Nully
    function nully!(f::Format, value)
        f.nully = value
    end

    # Prefix
    function si_prefix!(f::Format, value)
        check_named_value(Prefix, value)
        f.prefix = value.value
    end

    JSON3.StructTypes.StructType(::Type{Format}) = JSON3.RawType()

    function JSON3.rawbytes(f::Format)
        aligned = f.specifier[:align] != Align.default.value
        fill = aligned ? f.specifier[:fill] : ""
        spec_io = IOBuffer()
        print(spec_io,
            aligned ? f.specifier[:fill] : "",
            f.specifier[:align],
            f.specifier[:sign],
            f.specifier[:symbol],
            f.specifier[:padding],
            f.specifier[:width],
            f.specifier[:group],
            f.specifier[:precision],
            f.specifier[:trim],
            f.specifier[:type]
        )
        return JSON3.write(
            (
                locale = f.locale,
                nully = f.nully,
                prefix = f.prefix,
                specifier = String(take!(spec_io))
            )
        )
    end

    money(decimals, sign = Sign.default) = Format(
        group=Group.yes,
        precision=decimals,
        scheme=Scheme.fixed,
        sign=sign,
        symbol=Symbol.yes
    )

    function percentage(decimals, rounded::Bool=false)
        scheme = rounded ? Scheme.percentage_rounded : Scheme.percentage
        return Format(scheme = scheme, precision = decimals)
    end
end


text:
macro forward(ex, fs)
    T, prop = @match ex begin
        :($T.$prop) => (T, prop)
        _ => error("Invalid expression $ex, expected <Type>.<prop>")
    end

    fs = @match fs begin
        :(($(fs...),)) => fs
        _ => error("Expected a tuple of functions, got $fs")
    end

    defs = map(fs) do f
        esc(:($f(x::$T, args...; kwargs...) = $f(x.$prop, args...; kwargs...)))
    end

    Expr(:block, defs...)
end

macro inst(ex)
    result_id, inst = @match ex begin
        :($result_id = $inst) => (result_id, inst)
        _ => (nothing, ex)
    end

    type_id, call = @match inst begin
        :($call::$T) => (T, call)
        _ => (nothing, inst)
    end

    opcode, args = @match call begin
        :($opcode($(args...))) => (esc(opcode), esc.(args))
        _ => error("Invalid call $call")
    end

    :(Instruction($opcode, $(esc(type_id)), $(esc(result_id)), $(args...)))
end

isline(x) = false
isline(x::LineNumberNode) = true

function rmlines(ex)
    @match ex begin
        Expr(:macrocall, m, _...) => Expr(:macrocall, m, nothing, filter(x -> !isline(x), ex.args[3:end])...)
        ::Expr                    => Expr(ex.head, filter(!isline, ex.args)...)
        a                         => a
    end
end

macro broadcastref(ex)
    T = @match ex begin
        :(struct $T; $(fields...); end) => T
        :(abstract type $T end) => T
    end

    T = @match T begin
        :($T <: $AT) => T
        _ => T
    end

    T = @match T begin
        ::Symbol => T
        :($T{$(_...)}) => T
    end

    quote
        Base.@__doc__ $(esc(ex))
        Base.broadcastable(x::$(esc(T))) = Ref(x)
    end
end

function source_version(language::SourceLanguage, version::Integer)
    @match language begin
        &SourceLanguageGLSL => begin
            major = version ÷ 100
            minor = (version - 100 * major) ÷ 10
            VersionNumber(major, minor)
        end
    end
end

function source_version(language::SourceLanguage, version::VersionNumber)::UInt32
    @match language begin
        &SourceLanguageGLSL => begin
            @unpack major, minor = version
            10 * minor + 100 * major
        end
    end
end


text:
root_dir = "../."
using Pkg; Pkg.activate(root_dir)
using Strathreads, Temporal, Indicators, Dates
using DotEnv

pathCoins = "../data/test/"
outDir = "results"

_assets = Base.Filesystem.readdir("data/test/")
d = load("../data/opts.jld")
for c in _assets
    opts = d["macd"][c]
    assets = [c]

    universe = Universe(assets)
    function datasource(asset::String)::TS
        savedata_path = joinpath(pathCoins, asset)
        return Temporal.tsread(savedata_path, indextype=UInt64, format="yyyy-mm-dd HH:MM:SS")
    end
    universe = Universe(assets)
    universe = gather(assets, source=datasource)
    function fun(x::TS; args...)::TS
        close_prices = x[:close]
        macd_ma = macd(close_prices; args...)
        output = [close_prices macd_ma]
        #output.fields = [:close, :MACD] -> shit renames fields
        return output
    end
    arg_names = [:nfast, :nslow]
    arg_defaults = [trunc(Int, opts[1]),
                    trunc(Int, opts[2])]
    arg_ranges = [5:1:12, 15:1:35]
    paramset = ParameterSet(arg_names, arg_defaults, arg_ranges)
    indicator = Indicator(fun, paramset)

    siglong = @signal  MACD ↑ Signal
    sigexit = @signal MACD == Signal
    sigshort = @signal MACD ↓ Signal

    # define the trading rules
    longrule = @rule siglong → long 100
    shortrule = @rule sigshort → short 100
    exitrule = @rule sigexit → liquidate 1.0
    rules = (longrule, shortrule, exitrule)

    # run strategy
    strat = Strategy(universe, indicator, rules)
    bt = backtest(strat, px_trade=:open, px_close=:close, verbose=true)
    tswrite(bt[c], joinpath(outDir, c))
end

max_stars_repo_name: lithiferous/Strathreads

max_issues_repo_name: lithiferous/Strathreads

max_forks_repo_name: lithiferous/Strathreads


text:
import Interpolations
"""
	interpdf(data,timecol,timvec)

Re-sample dataframe applying linear interpolation to all columns

**Input**
* data: DataFrame where at least one column contains DateTime
* timevec: output time vector (=Vector{DateTime,1})
* timecol: column containing DateTime (default value = :datetime)

**Output**
* dataframe containing all re-sampled columns + timevec

**Example**
```
data = DataFrame(Temp=[10,11,12,14],
			datetime=[DateTime(2010,1,1,0),DateTime(2010,1,1,1),
			DateTime(2010,1,1,2),DateTime(2010,1,1,4)]);
datai = interpdf(data,[DateTime(2010,1,1,12,0,0)]);

```
"""
function interpdf(data::DataFrame,timevec::Vector{DateTime};timecol=:datetime)
	# prepare input time vector for interpolation
	x,xi = preptime(data[timecol],timevec);
	# run for all input columns except for DateTime (x vector) and Types not
	# suitable for interpolation
	dfi = DataFrame(datetime=timevec);
	for i in names(data)
		if (i != timecol) && (eltype(data[i]) <: Real)
			dfi[i] = interp1(x,data[i],xi);
		end
	end
    return dfi
end

"""
	interp1(x,y,xi)
Interpolate 1D data

**Input**
* x: x coordinate vector
* y: y values (vector)
* xi: x coordinate to be interpolated

**Output**
* yi: interpolated value

**Example**
```
yi = interp1([1,2,3,4],[10,20,30,40],1.5);
```
"""
function interp1(x,y,xi)
	# Declare output vector
	out = zeros(length(xi));
	# Create interpolation object
    itp = Interpolations.interpolate((x,),y,
            Interpolations.Gridded(Interpolations.Linear()));
	# exclude (set to NA) extrapolated values
	maxx = maximum(x);minx = minimum(x);
	for i = 1:length(out)
        if (xi[i] <= maxx) && (xi[i] >= minx)
            out[i] = itp(xi[i]);
		else
			out[i] = NaN;
        end
    end
    return out;
end

function interp1(x::Vector{DateTime},y::Vector{Float64},xi::Vector{DateTime})
	t1,t2 = preptime(x,xi);
	return interp1(t1,y,t2)
end
"""
	preptime(dft,timevec)

Auxiliary function to prepare time in DateTime format for interpolation
Will convert DateTime to Number/Int and declare/create output (length) dataframe

"""
function preptime(dft::Vector{DateTime},timevec::Vector{DateTime})
	x = Dates.value.(dft);
	xi = Dates.value.(timevec);
	return x, xi
end

"""
	meshgrid(x,y)

Function to create meshgrid, i.e. returns the same as 2D Matlab meshgrid [function](https://de.mathworks.com/help/matlab/ref/meshgrid.html).

**Input**
* x: vector or range (Int or Float)
* y: vector or range (Int or Float)

**Output**
* xi: matrix with x coordinates  (Int or Float)
* yi: matrix with y coordinates  (Int or Float)

**Example**

```
xi,yi = meshgrid(0:1:10,100:1:110);

```

"""
function meshgrid(x,y)
    xi = [j for i in y, j in x];
    yi = [i for i in y, j in x];
    return xi, yi;
end

"""
	mesh2vec(xi,yi)

Convert meshgrid back to vector

**Input**
* xi: matrix of x coordinates (in mesghrid format)
* yi: matrix of y coordinates (in mesghrid format)

**Output**
* x: vector of x coordinates
* y: vector of x coordinates

**Example**

```
xi = [0 1 2;0 1 2;0 1 2];
yi = [10 10 10;11 11 11;12 12 12];
x,y = mesh2vec(xi,yi);

```
"""
function mesh2vec(x,y)
    return x[1,:], y[:,1]
end

"""
	interp2(x,y,z,xi,yi)

Re-sample 2D data using linear interpolation

**Input**
* x: vector with x coordinates
* y: vector with y coordinates
* z: matrix with (x,y)-dependent values (meshgrid-like)
* xi: new x coordinates
* yi: new y coordinates

**Output**
* zi: interpolated values

**Example**
```
# Prepare input
x = [1.,2.,3.,4.,5.,];
y = [2.,3.,4.,5.,6.,7.,8.];
z = ones(Float64,(length(y),length(x)));
z[2,3] = 2.;
xi = [4.1,3.0,3.0,3.5,5.0,5.01]
yi = [2.1,3.0,3.5,3.5,2.0,3.30]
# compute
zi = interp2(x,y,z,xi,yi)
```
"""
function interp2(x::Vector{Float64},y::Vector{Float64},z::Matrix{Float64},
				 xi::Vector{Float64},yi::Vector{Float64})
	minx,miny,maxx,maxy,itp = prep2interp(x,y,z);
	# declare output
	zi = zeros(Float64,size(xi));
	for i in 1:length(xi[:])
		# do not extrapolate!
		if (xi[i] >= minx && yi[i] >= miny) && (xi[i] <= maxx && yi[i] <= maxy)
		   zi[i] = itp(yi[i],xi[i]) # use yi,xi order as matlab sorts masgrid in different order
		else
		   zi[i] = NaN;
		end
	end
	return zi
end
function interp2(x::Vector{Float64},y::Vector{Float64},z::Matrix{Float64},
				 xi::Matrix{Float64},yi::Matrix{Float64})
	zi = interp2(x,y,z,xi[:],yi[:]);
	return reshape(zi,size(xi));
end
function interp2(x::Matrix{Float64},y::Matrix{Float64},z::Matrix{Float64},
				 xi::Matrix{Float64},yi::Matrix{Float64})
	xm,ym = mesh2vec(x,y);
	return interp2(xm,ym,z,xi,yi);
end
function interp2(x::Matrix{Float64},y::Matrix{Float64},z::Matrix{Float64},
				 xi::Float64,yi::Float64)
	xm,ym = mesh2vec(x,y);
	return interp2(xm,ym,z,xi,yi);
end
function interp2(x::Vector{Float64},y::Vector{Float64},z::Matrix{Float64},
				 xi::Float64,yi::Float64)
	return interp2(x,y,z,[xi],[yi])[1]
end

"""
	Auxilliary function for determination of min and max values of input vectors
"""
function prep2interp(x::Vector{Float64},y::Vector{Float64},z::Matrix{Float64})
    # Need to switch x <--> y as matlab meshgrid uses transposed matrix
    itp = Interpolations.interpolate((y,x),z,
            Interpolations.Gridded(Interpolations.Linear()));
    return minimum(x), minimum(y), maximum(x), maximum(y), itp;
end


text:
"""
    Duality(optimizer)

Duality uses Lagrangian relaxation to compute over-approximated bounds for a network

# Problem requirement
1. Network: any depth, any activation function that is monotone
2. Input: hyperrectangle
3. Output: halfspace

# Return
`BasicResult`

# Method
Lagrangian relaxation.
Default `optimizer` is `GLPKSolverMIP()`.

# Property
Sound but not complete.

# Reference
[K. Dvijotham, R. Stanforth, S. Gowal, T. Mann, and P. Kohli,
"A Dual Approach to Scalable Verification of Deep Networks,"
*ArXiv Preprint ArXiv:1803.06567*, 2018.](https://arxiv.org/abs/1803.06567)
"""
@with_kw struct Duality <: Solver
    optimizer = GLPK.Optimizer
end

function solve(solver::Duality, problem::Problem)

    @assert is_halfspace_equivalent(problem.output) "Duality only accepts HalfSpace output sets. Got a $(length(d))-d polytope."

    model = Model(solver)
    c, d = tosimplehrep(problem.output)
    d = d[1]

    λ = init_vars(model, problem.network, :λ)
    μ = init_vars(model, problem.network, :μ)
    o = dual_value(solver, problem, model, λ, μ)
    @constraint(model, last(λ) .== -c)
    optimize!(model)

    # Interpret result:
    if termination_status(model) != OPTIMAL
        return BasicResult(:unknown)
    elseif value(o) - d <= 0.0
        return BasicResult(:holds)
    else
        return BasicResult(:violated)
    end
end

function dual_value(solver::Duality,
                    problem::Problem,
                    model::Model,
                    λ::Vector{Vector{VariableRef}},
                    μ::Vector{Vector{VariableRef}})
    bounds = get_bounds(problem)
    layers = problem.network.layers

    # prepend an array of 0s to λ so the λᵢ
    # term cancels in the input layer
    λ = [zeros(dim(bounds[1])), λ...]

    o = 0
    for i in 1:length(layers)
        Lᵢ, Bᵢ = layers[i], bounds[i]
        W, b = Lᵢ.weights, Lᵢ.bias
        c, r = Bᵢ.center, Bᵢ.radius
        B̂ᵢ₊₁ = approximate_affine_map(Lᵢ, Bᵢ)

        # layer value
        o += λ[i]'*c - μ[i]'*(W*c + b) + sum(symbolic_abs.(λ[i] .- W'*μ[i]) .* r)

        # activation value
        o += activation_value(Lᵢ.activation, μ[i], λ[i+1], low(B̂ᵢ₊₁), high(B̂ᵢ₊₁))
    end

    @objective(model, Min, o)
    return o
end

function activation_value(σ::ReLU, μᵢ, λᵢ, l̂ᵢ, ûᵢ)

    gᵢl̂ᵢ = @. μᵢ*l̂ᵢ - λᵢ*σ(l̂ᵢ)
    gᵢûᵢ = @. μᵢ*ûᵢ - λᵢ*σ(ûᵢ)

    max = symbolic_max

    # NOTE: ifelse evaluates all of its arguments in advance.
    # This may seem like a bug since symbolic_max mutates the model
    # and introduces new variables, but it is fine. Only those
    # variables that are part of the objective will affect the problem,
    # and those will only the the ones summed *after* taking the ifelse
    return sum(@. ifelse(l̂ᵢ < 0 < ûᵢ,
                         max(gᵢl̂ᵢ, gᵢûᵢ, 0),
                         max(gᵢl̂ᵢ, gᵢûᵢ)))
end

function activation_value(σ::Id, μᵢ, λᵢ, l̂ᵢ, ûᵢ)
    max = symbolic_max
    sum(@. max(μᵢ*l̂ᵢ - λᵢ*l̂ᵢ, μᵢ*ûᵢ - λᵢ*ûᵢ))
end

function activation_value(σ::Any, μᵢ, λᵢ, l̂ᵢ, ûᵢ)
    max = symbolic_max
    sum(@. max(μᵢ*l̂ᵢ, μᵢ*ûᵢ) + max(-λᵢ*σ(l̂ᵢ), -λᵢ*σ(ûᵢ)))
end


text:
module UnROOT

export ROOTFile, array

import Base: keys, get, getindex, show, length, iterate, position

using CodecZlib, CodecLz4, CodecXz
using Mixers
using Parameters
using StaticArrays

include("constants.jl")
include("utils.jl")
include("io.jl")
include("types.jl")
include("streamers.jl")
include("bootstrap.jl")
include("root.jl")
include("custom.jl")

end # module


text:
using JLD2, FileIO, Test, Base64

# In #206 the encoding of Union types was changed
# but old files should still be readable
@testset "Legacy Union Type" begin
    # Test this by loading a file created with the most recent
    # JLD2 version before that change v0.1.14
    # Code:
    #   jldopen(fn, "w") do f
    #       f["U"] = Union{Int, String}
    #       f["V"] = Union{Int, Float64}[1, 2.0, 3]
    #   end
    # Here encoded in base64.
    oldfileb64 = """[encoded data, 7500 characters]"""

    fn = joinpath(mktempdir(), "test.jld2")
    write(fn, base64decode(oldfileb64))

    u, v = jldopen(fn, "r") do f
        f["U"], f["V"]
    end

    @test u == Union{Int64, String}
    @test v == Union{Int64, Float64}[Int64(1), 2.0, Int64(3)]
end


text:
@testset "op_dd_dd" begin

    @test @df(neg_dd_dd, HILO(gamma_df64)) == HILO(-gamma_df64)
    @test @df(abs_dd_dd, HILO(-gamma_df64)) == HILO(gamma_df64)
    @test sum(@df(negabs_dd_dd, HILO(-(gamma_df32))) .- HILO(-(abs(-gamma_df32)))) < 1.0e-30
    @test sum(abs.(HILO(gamma_df64) .- @df(inv_dd_dd, @df(inv_dd_dd, HILO(gamma_df64))))) < 1.0e-30

end


text:
import COBRA

export reactionsIrrev, addFluxCost, remove_reaction, remove_reactions

"""
Creates an extended model, where reactions are split in two,
forward and backward with non-negative fluxes. For each
reaction in the original model, a new reaction is appended
at the end which represents the backward flux.
"""
function reactionsIrrev(cobra_model::COBRA.LPproblem)
    S = [(cobra_model.S) (-cobra_model.S)]
    lb = [max.(0, cobra_model.lb); max.(0, -cobra_model.ub)]
    ub = [max.(0, cobra_model.ub); max.(0, -cobra_model.lb)]
    b = copy(cobra_model.b)
    c = [cobra_model.c; -cobra_model.c]
    rxns = [cobra_model.rxns; cobra_model.rxns .* "_bk"]
    COBRA.LPproblem(S, b, c, lb, ub, cobra_model.osense,
                    copy(cobra_model.csense), rxns,
                    copy(cobra_model.mets))
end

"""
Constrains the sum of all fluxes (weighted by costs) to be less than 1.
Adds a pseudometabolite that is produced by all reactions in the original
model, and a pseudoreaction that is the only consumer of this pseudometabolite.
Then limits the flux of the pseudoreaction to be less than 1.
"""
function addFluxCost(cobra_model::COBRA.LPproblem, costs::Vector{Float64})
    (m,n) = size(cobra_model.S)
    @assert length(costs) == n
    S = [(cobra_model.S) zeros(m); costs' -1]
    lb = [cobra_model.lb; 0]
    ub = [cobra_model.ub; 1]
    b = [cobra_model.b; 0]
    c = [cobra_model.c; 0]
    rxns = [cobra_model.rxns; "PseudoRxnCost"]
    mets = [cobra_model.mets; "PseudoMetCost"]
    csense = [cobra_model.csense; '=']
    COBRA.LPproblem(S, b, c, lb, ub,
                    cobra_model.osense, csense,
                    rxns, mets)
end

"""
Index of reaction rxn.
"""
function reaction_index(cobra_model::COBRA.LPproblem, rxn::String)
    k = findfirst(cobra_model.rxns, rxn)
    k == 0 && throw(KeyError(rxn))
    return k
end

"""
Index of metabolite met.
"""
function metabolite_index(cobra_model::COBRA.LPproblem, met::String)
    i = findfirst(cobra_model.mets, met)
    i == 0 && throw(KeyError(met))
    return i
end

"""
Returns list of metabolite indices that participate in the k'th reaction.
"""
function reaction_metabolites(cobra_model::COBRA.LPproblem, k::Int)
    @assert 1 ≤ k ≤ length(cobra_model.rxns)
    find(cobra_model.S[1:end, k])
end

"""
Returns list of metabolite indices that are consumed in the k'th reaction.
"""
function reaction_substrates(cobra_model::COBRA.LPproblem, k::Int)
    @assert 1 ≤ k ≤ length(cobra_model.rxns)
    find(s -> s < 0, cobra_model.S[1:end, k])
end

"""
Returns list of metabolite indices that are produced in the k'th reaction.
"""
function reaction_products(cobra_model::COBRA.LPproblem, k::Int)
    @assert 1 ≤ k ≤ length(cobra_model.rxns)
    find(s -> s > 0, cobra_model.S[1:end, k])
end

reaction_metabolites(cobra_model::COBRA.LPproblem, rxn::String) = reaction_metabolites(cobra_model, reaction_index(cobra_model, rxn))
reaction_substrates(cobra_model::COBRA.LPproblem, rxn::String)  = reaction_substrates(cobra_model,  reaction_index(cobra_model, rxn))
reaction_products(cobra_model::COBRA.LPproblem, rxn::String)    = reaction_products(cobra_model,    reaction_index(cobra_model, rxn))

"""
List of reaction indices that involve the i'th metabolite.
"""
function metabolite_reactions end

function metabolite_reactions(cobra_model::COBRA.LPproblem, i::Integer)
    @assert 1 ≤ i ≤ length(cobra_model.mets)
    find(cobra_model.S[i,:])
end

metabolite_reactions(cobra_model::COBRA.LPproblem, met::String) = metabolite_reactions(cobra_model, metabolite_index(cobra_model, met))

"""
Returns list of reactions that exchange the i'th metabolite
"""
function exchange_reactions(cobra_model::COBRA.LPproblem, i::Integer)
    indexes = Int[]
    for k in metabolite_reactions(cobra_model, i)
        if isempty(reaction_substrates(cobra_model, k)) || isempty(reaction_products(cobra_model, k))
            push!(indexes, k)
        end
    end
    return indexes
end

exchange_reactions(cobra_model::COBRA.LPproblem, met::String) = exchange_reactions(cobra_model, metabolite_index(cobra_model, met))

"""
Returns a copy of cobra_model where reaction i has been eliminated
"""
remove_reaction(cobra_model::COBRA.LPproblem, i::Integer) = remove_reactions(cobra_model, [i])

"""
Returns a copy of cobra_model where the reactions in list have been eliminated
"""
function remove_reactions(cobra_model::COBRA.LPproblem, list::AbstractVector{Int})
    keep = [k for k = 1 : length(cobra_model.rxns) if k ∉ list]
    keep_reactions(cobra_model, keep)
end

"""
Returns a copy of cobra_model where all reactions have been eliminated except those in list.
"""
function keep_reactions(cobra_model::COBRA.LPproblem, list::AbstractVector{Int})
    @assert all(1 .≤ list .≤ length(cobra_model.rxns))
    list = unique(list)

    S = cobra_model.S[:, list]
    c = cobra_model.c[list]
    lb = cobra_model.lb[list]
    ub = cobra_model.ub[list]
    rxns = cobra_model.rxns[list]

    COBRA.LPproblem(S, cobra_model.b, c, lb, ub, cobra_model.osense,
                    cobra_model.csense, rxns, cobra_model.mets)
end

"""
Returns a copy of model where metabolites in list have been eliminated.
"""
function remove_metabolites(model::COBRA.LPproblem, list::AbstractVector{Int})
    keep = [i for i = 1 : length(model.mets) if i ∉ list]
    keep_metabolites(model, keep)
end

"""
Returns a copy of model where metabolites in list have been eliminated.
"""
remove_metabolite(model::COBRA.LPproblem, i::Integer) = remove_metabolites(model, [i])

"""
Returns a copy of model where all metabolites have been eliminated except those in list.
"""
function keep_metabolites(model::COBRA.LPproblem, list::AbstractVector{Int})
    @assert all(1 .≤ list .≤ length(model.mets))
    list = unique(list)

    S = model.S[list, :]
    b = model.b[list]
    csense = model.csense[list]
    mets = model.mets[list]

    COBRA.LPproblem(S, b, model.c, model.lb, model.ub,
                    model.osense, csense, model.rxns, mets)
end

"""
Returns a copy of model where reactions that can't carry flux have been removed.
Null reactions are those for which max(|fmin|, |fmax|) < tol, where tol = 1e-10
by default. This requires fmin, fmax from fva. If you have computed fmin, fmax
already, you can pass them as arguments to save computing time.
If fixbounds = true (default), sets lb and ub of reactions to fmin, fmax.
"""
function reduce_model end

function reduce_model(model::COBRA.LPproblem; tol::Real = 1e-10, fixbounds::Bool = true)
    @assert 0 ≤ tol < Inf
    fmin, fmax = fva(model)
    reduce_model(model, fmin, fmax; tol=tol, fixbounds=fixbounds)
end

function reduce_model(model::COBRA.LPproblem, fmin::Vector{Float64}, fmax::Vector{Float64};
                      tol::Real = 1e-6, fixbounds::Bool = true)
    @assert 0 ≤ tol < Inf
    @assert length(fmin) == length(fmax) == length(model.rxns)
    model = deepcopy(model)
    if fixbounds
        model.lb .= fmin
        model.ub .= fmax
    end
    rxns = [k for k = 1 : length(model.rxns) if abs(fmin[k]) < tol && abs(fmax[k]) < tol]
    model = remove_reactions(model, rxns)
    mets = [i for i = 1 : length(model.mets) if length(metabolite_reactions(model, i)) == 0]
    model = remove_metabolites(model, mets)
    return model
end

"""
Sets the objective reaction as the i'th reaction.
If sense == 1 (default), flux is maximized. If
sense == -1, flux is minimized.
"""
function set_objective!(model::COBRA.LPproblem, i::Integer, sense::Integer = 1)
    model.c .= 0
    model.c[i] = 1
    model.osense = sign(sense)
end

lower_bound(model::COBRA.LPproblem, rxn::String) = model.lb[COBRAUtils.reaction_index(model, rxn)]
upper_bound(model::COBRA.LPproblem, rxn::String) = model.ub[COBRAUtils.reaction_index(model, rxn)]


text:
using Documenter, MixedModels, StatsBase

makedocs(
    sitename = "MixedModels",
    pages = ["index.md",
             "constructors.md",
             "optimization.md",
             "GaussHermite.md",
             "bootstrap.md",
             "SimpleLMM.md",
             "MultipleTerms.md",
             "SingularCovariance.md",
             "SubjectItem.md",
             "benchmarks.md"]
)

deploydocs(
    repo    = "github.com/dmbates/MixedModels.jl.git"
)


text:
using Test
using Pkg
using Comonicon

Pkg.activate(pkgdir(Comonicon, "test", "projects", "Hello"))
Pkg.develop(path = pkgdir(Comonicon))
Pkg.test()

Pkg.activate(pkgdir(Comonicon, "test", "projects", "FakePkg"))
Pkg.develop(path = pkgdir(Comonicon))
Pkg.test()


text:
"""
Clouds months tends to have lower radiance due to attenuation. The bias_correction function uses the number of cloud-free observations to adjust the radiance accordingly.
```julia
bias_correction(radiance, clouds)
```
"""
function bias_correction(radiance, clouds::Array{T, 1}, smoothing_parameter=10.0) where T <: Any
    missings = findall(ismissing, radiance)
    y = filter!(!ismissing, copy(radiance))
    x = Array{Union{Float64, Missing}}(copy(clouds))
    x[missings] .= missing
    x = Array{Float64}(filter!(!ismissing, x))
    y_ = Array{Float64}(detrend_ts(y))
    m = predict(fit(SmoothingSpline, x, y_, smoothing_parameter), x)
    topval = maximum(m)
    fixed1 = y_ -m .+ topval
    fixed2 = fixed1 + y - y_ # adding the time trend
    for i in 1:length(y)
        if fixed2[i] > maximum(y)
            fixed2[i] = y[i]
        end
    end
    fixed2 = Array{Union{Float64, Missing}}(fixed2)
    for i in missings #putting back the missing values
        fixed2 = insert!(fixed2, i, missing)
    end
    return fixed2
end

"""
The bias correction function can use the datacubes of radiance and the number of cloud-free observations to correct for attenuation in radiance due to low number of cloud-free observations.
```julia
bias_correction(radiance, clouds)
```
"""
function bias_correction(radiance_datacube, clouds_datacube::Array{T, 3}, mask=ones(Int8, (size(radiance_datacube)[1],size(radiance_datacube)[2]))) where T <: Any
    rad_corrected_datacube = radiance_datacube
    @showprogress for i in 1:size(radiance_datacube)[1]
        for j in 1:size(radiance_datacube)[2]
            if count(i->(ismissing(i)),radiance_datacube[i, j, :])/length(radiance_datacube[i, j, :]) > 0.50
                continue
            end
            if mask[i, j]==0
                continue
            end
            missings = findall(ismissing, radiance_datacube[i, j, :])
            radiance_arr = filter!(!ismissing, copy(radiance_datacube[i, j, :]))
            clouds_arr = Array{Union{Float64, Missing}}(copy(clouds_datacube[i , j, :]))
            clouds_arr[missings] .= missing
            clouds_arr = filter!(!ismissing, clouds_arr)
            if rank_correlation_test(radiance_arr, clouds_arr) <0.05
                rad_corrected_datacube[i, j, :]= bias_correction(copy(radiance_datacube[i, j, :]), copy(clouds_datacube[i , j, :]))
            else
                rad_corrected_datacube[i, j, :]= radiance_datacube[i, j, :]
            end
        end
    end
    return rad_corrected_datacube
end


text:
using PayGroups
using Test

@testset "PayGroups.jl" begin
    # Write your tests here.
end


text:
#julia -p 10 -L rollingHarness.jl test_Rolling2.jl outPathStub d
#passed arguments
	#ARGS[1] is partial path for output.
	#ARGS[2] is d must be one of 20, 50, 1000
#Use this grid of N
#N_grid = [10, 20, 40]

using Distributed

const spath = ARGS[1]
const d = parse(Int, ARGS[2])
@assert d in [20, 50, 1000] "Only bins of size 20, 50, 1000 currently supported"

const param_path = "AdjSales_NoWeekends_ShuffleWithinCol"
const s = .95
K_grid = vcat(round.(Int, 2 .^(3:.25:10)), 1115)
outPath = "$(spath)_Ross_$(s)__$(d)_$(param_path)"

#First read in the data and parse it appropriately
ps_full = readdlm("../RossmanKaggleData/CleanedData/ps_full$(d).csv", ',')
supp_full = readdlm("../RossmanKaggleData/CleanedData/support$(d).csv", ',')

#Load shuffled data
tdata = readdlm("../RossmanKaggleData/CleanedData/$(param_path)_Binned$(d).csv", ',')

binned_data = tdata[2:end, 2:end]  #drop column header = stores, drop row header = dates
dates= tdata[2:end, 1] #keep track of dates for fun

#Separate into 10 processors.
#Divy them
half_pts = floor(Int, size(binned_data, 1)/2)
start_time = time_ns()

#For N=10, 20 split K_grid into 2 and data into 2 (sacrifice 1 point)
file10_aa = @spawn rollingTest(K_grid[1:2:end], supp_full, ps_full,
					 binned_data[1:half_pts, :], dates[1:half_pts],
					 "$(outPath)_thread_10_a_a",
					 [10], s)
file10_ab = @spawn rollingTest(K_grid[1:2:end], supp_full, ps_full,
					 binned_data[(half_pts + 1):end, :], dates[(half_pts + 1):end],
					 "$(outPath)_thread_10_a_b",
					 [10], s)
file10_ba = @spawn rollingTest(K_grid[2:2:end], supp_full, ps_full,
					 binned_data[1:half_pts, :], dates[1:half_pts],
					 "$(outPath)_thread_10_b_a",
					 [10], s)
file10_bb = @spawn rollingTest(K_grid[2:2:end], supp_full, ps_full,
					 binned_data[(half_pts + 1):end, :], dates[(half_pts + 1):end],
					 "$(outPath)_thread_10_b_b",
					 [10], s)

file20_aa = @spawn rollingTest(K_grid[1:2:end], supp_full, ps_full,
					 binned_data[1:half_pts, :], dates[1:half_pts],
					 "$(outPath)_thread_20_a_a",
					 [20], s)
file20_ab = @spawn rollingTest(K_grid[1:2:end], supp_full, ps_full,
					 binned_data[(half_pts + 1):end, :], dates[(half_pts + 1):end],
					 "$(outPath)_thread_20_a_b",
					 [20], s)
file20_ba = @spawn rollingTest(K_grid[2:2:end], supp_full, ps_full,
					 binned_data[1:half_pts, :], dates[1:half_pts],
					 "$(outPath)_thread_20_b_a",
					 [20], s)
file20_bb = @spawn rollingTest(K_grid[2:2:end], supp_full, ps_full,
					 binned_data[(half_pts + 1):end, :], dates[(half_pts + 1):end],
					 "$(outPath)_thread_20_b_b",
					[20], s)

#N = 40, just split K = 2 processors
file40_a = @spawn rollingTest(K_grid[1:2:end], supp_full, ps_full,
					 binned_data, dates,
					 "$(outPath)_thread_40_a",
					[40], s)
file40_b = @spawn rollingTest(K_grid[2:2:end], supp_full, ps_full,
					 binned_data, dates,
					 "$(outPath)_thread_40_b",
					[40], s)

# ######
file10_aa = fetch(file10_aa)
file10_ab = fetch(file10_ab)
file10_ba = fetch(file10_ba)
file10_bb = fetch(file10_bb)

file20_aa = fetch(file20_aa)
file20_ab = fetch(file20_ab)
file20_ba = fetch(file20_ba)
file20_bb = fetch(file20_bb)

file40_a = fetch(file40_a)
file40_b = fetch(file40_b)

time_stamp = (time_ns() - start_time) * 1e-9

##read everyone in, throw away a line
data_10_aa, header = readdlm(file10_aa, ',', header=true)
data_10_ab, header = readdlm(file10_ab, ',', header=true)
data_10_ba, header = readdlm(file10_ba, ',', header=true)
data_10_bb, header = readdlm(file10_bb, ',', header=true)

data_20_aa, header = readdlm(file20_aa, ',', header=true)
data_20_ab, header = readdlm(file20_ab, ',', header=true)
data_20_ba, header = readdlm(file20_ba, ',', header=true)
data_20_bb, header = readdlm(file20_bb, ',', header=true)

data_40_a, header = readdlm(file40_a, ',', header=true)
data_40_b, header = readdlm(file40_b, ',', header=true)

data = reduce(vcat, (data_10_aa, data_10_ab, data_10_ba, data_10_bb,
					 data_20_aa, data_20_ab, data_20_ba, data_20_bb,
					 data_40_a, data_40_b) )

f = open("$(outPath).csv", "w")
writedlm(f, header, ',')
writedlm(f, data, ',')
close(f)

println("Time:", time_stamp )
println("$(outPath).csv")


text:
module TestProjectExtras

include("preamble.jl")
using Aqua: is_julia12_or_later, ispass, ⊜
using Base: PkgId, UUID

@testset "is_julia12_or_later" begin
    @test is_julia12_or_later("1.2")
    @test is_julia12_or_later("1.3")
    @test is_julia12_or_later("1.3, 1.4")
    @test is_julia12_or_later("1.3 - 1.4, 1.6")
    @test !is_julia12_or_later("1")
    @test !is_julia12_or_later("1.1")
    @test !is_julia12_or_later("1.0 - 1.1")
    @test !is_julia12_or_later("1.0 - 1.3")
end

with_sample_pkgs() do

    results = Dict(zip(
        [p.name for p in AquaTesting.SAMPLE_PKGIDS],
        Aqua.analyze_project_extras(collect(AquaTesting.SAMPLE_PKGIDS)),
    ))
    pkgids = Dict([p.name => p for p in AquaTesting.SAMPLE_PKGIDS])

    @testset "PkgWithIncompatibleTestProject" begin
        r = results["PkgWithIncompatibleTestProject"]
        @test !ispass(r)
        @test r ⊜ false
        msg = sprint(show, "text/plain", r)
        @test occursin("Root and test projects are inconsistent.", msg)
        @test occursin("Test dependencies not in root project", msg)
        @test occursin("Dependencies not in test project", msg)
        @test occursin("Random =", msg)
        @test occursin("REPL =", msg)
        @test !occursin("Test =", msg)
    end

    @testset "PkgWithCompatibleTestProject" begin
        r = results["PkgWithCompatibleTestProject"]
        @test ispass(r)
        Aqua.test_project_extras(pkgids["PkgWithCompatibleTestProject"])
        msg = sprint(show, "text/plain", r)
        @test occursin("Root and test projects are consistent.", msg)
        @test occursin("Root project:", msg)
        @test occursin("Test project:", msg)
    end

    @testset "PkgWithPostJulia12Support" begin
        r = results["PkgWithPostJulia12Support"]
        @test ispass(r)
        Aqua.test_project_extras(pkgids["PkgWithPostJulia12Support"])
        msg = sprint(show, "text/plain", r)
        @test occursin("Supporting only post-1.2 `julia`", msg)
        @test occursin("ignoring root project", msg)
    end

    @testset "PkgWithoutTestProject" begin
        r = results["PkgWithoutTestProject"]
        @test ispass(r)
        Aqua.test_project_extras(pkgids["PkgWithoutTestProject"])
        msg = sprint(show, "text/plain", r)
        @test occursin("test/Project.toml file does not exist", msg)
    end

    @testset "PkgWithoutProject" begin
        r = results["PkgWithoutProject"]
        @test !ispass(r)
        @test r ⊜ false
        msg = sprint(show, "text/plain", r)
        @test occursin("Project.toml file at project directory does not exist", msg)
    end

end

end  # module


text:
fi = open("../fv.bin", "r")
const pc_on_sq_old = read!( fi, Array{Int16}(tuple(pos_n, NumSQ)))
const kkp_old = read!( fi, Array{Int16}(tuple(kkp_end, NumSQ, NumSQ)))
#println("type of pc_on_sq = ", typeof(pc_on_sq))
#println("type of kkp = ", typeof(kkp))
close(fi)

const pc_on_sq = pc_on_sq_old
kkp = kkp_old
# kkp = zeros(Int16,NumSQ,NumSQ,kkp_end)
# for i = 1:kkp_end
#     for j = 1:NumSQ
#         for k = 1:NumSQ
#             kkp[k,j,i] = kkp_old[i,j,k]
#         end
#     end
# end

# println("pc_on_sq[0,0] = ", pc_on_sq[1,1])
# println("pc_on_sq[0,1] = ", pc_on_sq[1,2])
# println("pc_on_sq[1,0] = ", pc_on_sq[2,1])
# println("pc_on_sq[1,1] = ", pc_on_sq[2,2])
# println("kkp[0,0,0] = ", kkp[1,1,1])
# println("kkp[0,0,1] = ", kkp[1,1,2])
# println("kkp[0,1,0] = ", kkp[1,2,1])
# println("kkp[0,1,1] = ", kkp[1,2,2])

max_stars_repo_name: kimrin/MechajyoForWindows

max_issues_repo_name: kimrin/MechajyoForWindows

max_forks_repo_name: kimrin/MechajyoForWindows


text:
# This file is part of RLinearAlgebra.jl

module ProceduralTestLSVROneRandCyclic

using Test, RLinearAlgebra, Random

Random.seed!(1010)

@testset "LSVR One Hop Random Cyclic -- Procedural" begin

    # Verify appropriate super type
    @test supertype(LinSysVecRowOneRandCyclic) == LinSysVecRowSampler

    # Test whether row ordering remains fixed
    @test let
        A = rand(10,3)
        b = rand(10)
        x = rand(3)

        cyc = LinSysVecRowOneRandCyclic()

        α, β = RLinearAlgebra.sample(cyc, A, b, x, 1)

        order = copy(cyc.order)

        flag = true
        for j = 2:100
            α, β = RLinearAlgebra.sample(cyc, A, b, x, j)

            #Ordering should not change
            flag = flag & (cyc.order == order)
        end

        flag
    end

end

end


text:
using TimeSeries, MarketData, Base.Dates
FactCheck.setstyle(:compact)
FactCheck.onlystats(true)

facts("deprecated methods") do

    context("deprecated bydate methods correctly subset") do
        @fact by(mdata, 1, period=dayofweek).meta --> "Apple"

        @fact by(cl,4, period=day).timestamp[1]                 --> Date(2000,1,4)
        @fact by(cl, "Friday", period=dayname).timestamp[1]     --> Date(2000,1,7)
        @fact by(cl,5, period=week).timestamp[1]                --> Date(2000,1,31)
        @fact by(cl,5, period=month).timestamp[1]               --> Date(2000,5,1)
        @fact by(cl, "June", period=monthname).timestamp[1]     --> Date(2000,6,1)
        @fact by(cl,2001, period=year).timestamp[1]             --> Date(2001,1,2)
        @fact by(cl,1, period=dayofweek).timestamp[1]           --> Date(2000,1,3)
        # all the days in the nth week of each month
        @fact by(cl,5, period=dayofweekofmonth).timestamp[1]    --> Date(2000,1,31)
        @fact by(cl,365, period=dayofyear).timestamp[1]         --> Date(2001,12,31)
        @fact by(cl,4, period=quarterofyear).timestamp[1]       --> Date(2000,10,2)
        @fact by(cl,1, period=dayofquarter).timestamp[1]        --> Date(2001,10,1)
    end

    context("deprecated to / from methods select correctly") do
        @fact to(cl, 2000, 01, 01).values    --> to(cl, Date(2000, 01, 01)).values
        @fact to(cl, 2001, 01, 01).values    --> to(cl, Date(2001, 01, 01)).values
        @fact to(cl, 2002, 01, 01).values    --> to(cl, Date(2002, 01, 01)).values
        @fact from(cl, 2000, 01, 01).values  --> from(cl, Date(2000, 01, 01)).values
        @fact from(cl, 2001, 01, 01).values  --> from(cl, Date(2001, 01, 01)).values
        @fact from(cl, 2002, 01, 01).values  --> from(cl, Date(2002, 01, 01)).values
    end

    context("deprecated percentchange methods compute correct values") do
        @fact percentchange(op, method="simple").values  --> percentchange(op, :simple).values
        @fact percentchange(op, method="log").values     --> percentchange(op, :log).values
    end

    context("deprecated findall returns correct indices") do
        @fact findall(cl .> op) --> find(cl .> op)
    end

    context("deprecated collapse squishes correctly") do
        @fact collapse(cl, last).values[1]                  --> 99.50
        @fact collapse(cl, last).timestamp[1]               --> Date(2000,1,7)
        @fact collapse(cl, last, period=month).values[1]    --> 103.75
        @fact collapse(cl, last, period=month).timestamp[1] --> Date(2000,1,31)
    end
end


text:
const makiegallery_dir = dirname(dirname(@__DIR__))

const current_ref_version = Ref{String}("v0.4.2")

"""
    ref_image_dir(version = string(current_ref_version[]))

Returns the directory in which reference images are stored.
"""
function ref_image_dir(version = string(current_ref_version[]))
    return download_reference(version)
end

"Purges the reference image cache!"
purge_refimage_cache() = rm(joinpath(makiegallery_dir, "testimages"); recursive = true)

"""
    download_reference(version = string(current_ref_version[]))

Downloads the reference images from ReferenceImages for a specific version
"""
function download_reference(version = string(current_ref_version[]))
    refpath_version = version[1] == 'v' && version[2] in '0':'9' ? version[2:end] : version

    if isabspath(version)
        return version
    end

    download_dir = joinpath(makiegallery_dir, "testimages")
    isdir(download_dir) || mkpath(download_dir)
    tarfile = joinpath(download_dir, "gallery.tar.gz")
    url = "https://github.com/JuliaPlots/MakieReferenceImages/archive/$(version).tar.gz"
    refpath = joinpath(download_dir, "MakieReferenceImages-$(refpath_version)", "gallery")
    if !isdir(refpath) # if not yet downloaded
        download_images() = download(url, tarfile)
        try
            @info "downloading reference images for version $version"
            download_images()
        catch e
            if isa(e, ErrorException) && occursin("Hash Mismatch", e.msg)
                rm(tarfile, force = true)
                download_images()
            else
                rethrow(e)
            end
        end
        BinaryProvider.unpack(tarfile, download_dir)
        # check again after download
        if !isdir(refpath)
            error("Something went wrong while downloading reference images. Plots can't be compared to references")
        end
    end
    refpath
end

is_image_file(path) = lowercase(splitext(path)[2]) in (".png", ".jpg", ".jpeg")

function extract_frames(video, frame_folder)
    path = joinpath(frame_folder, "frames%04d.jpg")
    FFMPEG.ffmpeg_exe(`-loglevel quiet -i $video -y $path`)
end

function compare_media(a, b; sigma = [1,1], eps = 0.02)
    file, ext = splitext(a)
    if ext in (".png", ".jpg", ".jpeg", ".JPEG", ".JPG")
        imga = load(a)
        imgb = load(b)
        if size(imga) != size(imgb)
            @warn "images don't have the same size, difference will be Inf"
            return Inf
        end
        return approx_difference(imga, imgb, sigma, eps)
    elseif ext in (".mp4", ".gif")
        mktempdir() do folder
            afolder = joinpath(folder, "a")
            bfolder = joinpath(folder, "b")
            mkpath(afolder); mkpath(bfolder)
            extract_frames(a, afolder)
            extract_frames(b, bfolder)
            aframes = joinpath.(afolder, readdir(afolder))
            bframes = joinpath.(bfolder, readdir(bfolder))
            if length(aframes) > 10
                # we don't want to compare too many frames since it's time costly
                # so we just compare 10 random frames if more than 10
                samples = rand(1:length(aframes), 10)
                aframes = aframes[samples]
                bframes = bframes[samples]
            end
            # test by maximum diff
            return mean(compare_media.(aframes, bframes; sigma = sigma, eps = eps))
        end
    else
        error("Unknown media extension: $ext")
    end
end

"""
Compares all media recursively in two recorded folders!
"""
function run_comparison(
        test_record_path, test_diff_path,
        reference = MakieGallery.download_reference();
        maxdiff = 10.032
    )
    @testset "Reference Image Tests" begin
        folders = joinpath.(test_record_path, readdir(test_record_path))
        l = length(folders)
        count = 1
        for folder in folders
            @debug "Running index $count" progress=count/l
            count += 1
            if isdir(folder) && isdir(joinpath(reference, basename(folder)))
                media = joinpath(folder, "media")
                ref_folder = joinpath(reference, basename(folder), "media")
                test_folder = joinpath(test_record_path, media)
                ref_media = filter(x-> x != "thumb.jpg", sort(readdir(ref_folder)))
                test_media = filter(x-> x != "thumb.jpg", sort(readdir(test_folder)))
                @testset "$(basename(folder))" begin
                    if isempty(test_media)
                        @warn("recodings are missing for $folder")
                    else
                        for test in test_media
                            ref = joinpath(ref_folder, test)
                            test = joinpath(test_folder, test)
                            diff = compare_media(ref, test)
                            if diff >= maxdiff
                                refdiff = joinpath(test_diff_path, basename(folder), string("ref_", basename(ref)))
                                testdiff = joinpath(test_diff_path, basename(folder), string("test_", basename(test)))
                                mkpath(dirname(refdiff)); mkpath(dirname(testdiff))
                                cp(ref, refdiff, force = true)
                                cp(test, testdiff, force = true)
                                @test diff <= maxdiff
                            else
                                @test diff <= maxdiff
                            end
                        end
                    end
                end
            end
        end
    end
end

function load_test_database()
    database = load_database()
    # which one are the slowest and kicked those out!
    slow_examples = Set([
        "Animated time series",
        "Animation",
        "Lots of Heatmaps",
        "Chess Game",
        "Line changing colour",
        "Colormap collection",
        "Record Video",
        "Animated surface and wireframe",
        "Moire",
        "Line GIF",
        "Interaction with mouse",
        "Electrostatic repulsion",
        "pong",
        "pulsing marker",
        "Travelling wave",
        "Axis theming",
        "Legend",
        "Color Legend",
        "DifferentialEquations path animation",
        "Interactive Differential Equation",
        "Spacecraft from a galaxy far, far away",
        "WorldClim visualization",
        "Image on Geometry (Moon)",
        "Image on Geometry (Earth)",
        "Interaction with mouse",
        "Air Particulates",
        "Window resizing",
        "Aspect ratios stretching circles"
    ])
    # DIffeq errors with stackoverflow
    # The others look fine on the CI, but the measured difference is too high -.-
    # Maybe related to the axis changes, will investigate later
    filter!(database) do entry
        !("diffeq" in entry.tags) &&
        !(entry.unique_name in (:cobweb_plot, :analysis, :colormap_collection, :lots_of_heatmaps, :interaction_with_mouse, :normals_of_a_cat)) &&
        !(entry.title in slow_examples) &&
        !("download" in entry.tags)
     end
     return database
end


text:
using SparseMatrixAssemblers
using SparseArrays
using Test

@testset "Constructors" begin
    # Mainly testing that expressions compile
    @test typeof(SparseMatrixAssembler()) == SparseMatrixAssembler{Float64,Int}
    @test typeof(SparseMatrixAssembler{Int32}()) == SparseMatrixAssembler{Int32,Int}
    @test typeof(SparseMatrixAssembler{Int32,Int32}()) == SparseMatrixAssembler{Int32,Int32}
end

@testset "Conversion" begin
    # Mainly testing that expressions compile
    @test typeof(convert(SparseMatrixAssembler, SparseMatrixAssembler())) == SparseMatrixAssembler{Float64,Int}
    @test typeof(convert(SparseMatrixAssembler{Int32}, SparseMatrixAssembler())) == SparseMatrixAssembler{Int32,Int}
    @test typeof(convert(SparseMatrixAssembler{Int32,Int32}, SparseMatrixAssembler())) == SparseMatrixAssembler{Int32,Int32}
end

@testset "Assignment" begin
    @testset "Scalar" begin
        A = SparseMatrixAssembler{Int}()
        A[1,2] = 1
        @test sparse(A) == [0 1]
    end

    @testset "Col vector" begin
        A = SparseMatrixAssembler{Int}()
        A[1:2,2] = [1,2]
        @test sparse(A) == [0 1; 0 2]
    end

    @testset "Row vector" begin
        A = SparseMatrixAssembler{Int}()
        A[2,1:2] = [1,2]
        @test sparse(A) == [0 0; 1 2]
    end

    @testset "Matrix" begin
        A = SparseMatrixAssembler{Int}()
        A[2:3,1:2] = [1 2; 3 4]
        @test sparse(A) == [0 0; 1 2; 3 4]
    end
end

@testset "Conversion to other types" begin
    # Mainly testing that expressions compile
    i = [1,1]
    j = [1,1]
    v = [1,2]
    @test sparse(SparseMatrixAssembler(i,j,v)) == sparse(i,j,v)
    @test sparse(SparseMatrixAssembler(i,j,v),2,2) == sparse(i,j,v,2,2)
    @test sparse(SparseMatrixAssembler(i,j,v),2,2,*) == sparse(i,j,v,2,2, *)

    @test SparseMatrixCSC(SparseMatrixAssembler(i,j,v)) == sparse(i,j,v)
    @test SparseMatrixCSC(SparseMatrixAssembler(i,j,v),2,2) == sparse(i,j,v,2,2)
    @test SparseMatrixCSC(SparseMatrixAssembler(i,j,v),2,2,*) == sparse(i,j,v,2,2, *)

    @test Matrix(SparseMatrixAssembler(i,j,v)) == sparse(i,j,v)
    @test Matrix(SparseMatrixAssembler(i,j,v),2,2) == sparse(i,j,v,2,2)
    @test Matrix(SparseMatrixAssembler(i,j,v),2,2,*) == sparse(i,j,v,2,2, *)
end


text:
using KohnLeiboviciSzkup3
using Test

@testset "KohnLeiboviciSzkup3.jl" begin
    # Write your tests here.
end