text:
function background_evolve(Phi0,Pi0,pvals,funcs;verbose=false,
    options::BackgroundOptions = BackgroundOptions())
    Eh = funcs["Eh"]
    accel = funcs["Pi_eom"]
    V = funcs["V"]
    Hub = funcs["H"]
    d = length(Phi0)
    function eom(du,u,p,t)
        # name fields
        Phi  = u[1:d]
        Pi   = u[d+1:2d]
        par  = [Phi;Pi;pvals]
        #e    = Eh(par...)
        # eom
        du[1:d] = Pi
        du[d+1:2*d] = accel(par...)
    end
    function end_inflation(vars,t,integrator)
        if t > options.min_efolds
            e = Eh(vars...,pvals...)
            return options.end_epsilon-e # trigger when this hits zero from above
        else
            return 1.0
        end
    end
    e0 = Eh([Phi0;Pi0;pvals]...)
    if e0 >= 1.0 && iszero(options.min_efolds)
        #println("Initial ϵH = $e0")
        error("Inflation cannot start, initial eH = $e0 > 1")
    end
    H0 = Hub([Phi0;Pi0;pvals]...)
    if verbose
        vars = [Phi0;Pi0;pvals]
        ev0 = funcs["Ev"](vars...)
        et0 = funcs["Et"](vars...)
        o0 = funcs["Om"](vars...)
        println("Initial H: $H0")
        println("Initial ϵV: $ev0")
        println("Initial ϵH: $e0")
        println("Initial η: $et0")
        println("Initial Ω: $o0")
    end

    func = ODEFunction(eom)#;jac=eom_jac)
    SolT = eltype(Phi0)
    prob = ODEProblem(func,SolT[Phi0;Pi0],(zero(SolT),SolT(options.max_efolds)))
    cb = ContinuousCallback(end_inflation,terminate!)
    if verbose
        println("Solving background eom")
    end
    sol = DifferentialEquations.solve(prob,options.background_solver,callback=cb,
    reltol=options.background_rtol,abstol=options.background_atol,
    maxiters=options.max_steps,
    save_everystep=options.save_everystep,saveat=options.saveat)# automatically choose algorithm?
    if verbose
        println("Done. Nend = $(sol.t[end])")
    end
    Nend = sol.t[end]
    if Nend > options.Nexit_to_end && haskey(funcs,"ns_hca") && funcs["ns_hca"] != Nothing
        vars = [sol(Nend-options.Nexit_to_end);pvals]
        println("Estimated Pζ (HCA):\n $(funcs["Pz_hca"](vars...))")
        println("Estimated ns (HCA):\n $(funcs["ns_hca"](vars...))")
    end

    return sol
end

function geodesic(Phi0,Pi0,pvals,funcs;verbose=false,max_efolds=max_efolds)
    accel = funcs["Geo_eom"]
    function eom(du,u,p,t)
        # name fields
        Phi  = u[1:d]
        Pi   = u[d+1:2d]
        par  = [Phi;Pi;pvals]
        # eom
        du[1:d] = Pi
        du[d+1:2*d] = accel(par...)
    end

    func = ODEFunction(eom)
    prob = ODEProblem(func,Float64[Phi0;Pi0],(0.0,max_efolds))
    if verbose
        println("Solving background eom")
    end
    sol = DifferentialEquations.solve(prob,background_solver,reltol=background_rtol,abstol=background_atol,maxiters=max_steps,save_everystep=false,saveat=0.001)# automatically choose algorithm?
    if verbose
        println("Done. Nend = $(sol.t[end])")
    end

    return sol

end

function get_ks(sol,funcs,pvals,options)
    d = length(sol.u[1])÷2
    Hub(N) = funcs["H"]([sol(N);pvals]...)
    Nend = sol.t[end]
    Npivot = Nend - options.Nexit_to_end
    loga_init = log(options.kpivot) + log(Mpc_to_Lpl) - log(Hub(Npivot)) - Npivot

    Nkrange = options.Nkrange
    Nstep = 2*Nkrange / options.desired_ksteps
    Ns = Npivot-Nkrange:Nstep:Npivot+Nkrange
    Hs = [Hub(n) for n in Ns]
    logks = Ns .+ loga_init .+ log.(Hs)
    return loga_init,logks,Ns
end

function curved_nullspace(g,v)
    d = length(v)
    proj(u,v) = u * (u'*g*v) / (u'*g*u)
    norm(v) = sqrt(v'*g*v)
    ginv = inv(g)
    v ./= norm(v)
    perpuu = ginv .- v'.*v
    ws = [ perpuu[:,i] ./ norm(perpuu[:,i]) for i in 2:d ]
    ortho_ws = Array{eltype(v),2}(undef,d,d)
    ortho_ws[:,1] = v
    for i in 2:d
        w = ws[i-1]
        for j in 1:i-1
            w .-= proj(ortho_ws[:,j],w)
        end
        w ./= norm(w)
        ortho_ws[:,i] = w
    end
    return ortho_ws[:,2:end]
end

function transport_perturbations(sol,pvals,funcs;verbose=false,options::PerturbationOptions =PerturbationOptions())
    d = length(sol.u[1])÷2
    Nend = sol.t[end]
    Nsubhorizon = options.Nsubhorizon
    Nexit_to_end = options.Nexit_to_end
    Nkrange = options.Nkrange
    Nearly = Nend - Nexit_to_end - Nsubhorizon - Nkrange
    if Nearly < 0
        error("Not enough inflation! Attempted to evolve perturbations, but would've started at Ne = $Nearly")
    end
    logainit,logks,Nexits = get_ks(sol,funcs,pvals,options)
    if verbose
        println("Solving perturbations, with $(length(logks)) k-values")
        println("Nexit_to_end: $Nexit_to_end")
    end
    H(N) = funcs["H"]([sol(N);pvals]...)
    Mab(N) = funcs["Mab"]([sol(N)[1:2*d];pvals]...)
    Eh(N) = funcs["Eh"]([sol(N)[1:2*d];pvals]...)
    Ev(N) = funcs["Ev"](sol(N)...,pvals...)
    Vab(N) = reshape(funcs["Vab"](sol(N)...,pvals...),d,d)
    Om(N) = funcs["Om"]([sol(N)[1:2*d];pvals]...)
    Om2(N) = funcs["Om2"]([sol(N)[1:2*d];pvals]...)
    ginv(N) = funcs["ginv"]([sol(N)[1:2*d];pvals]...)
    g(N) = funcs["g"]([sol(N)[1:2*d];pvals]...)
    #eom = funcs["Pert_eom"]
    uab(t,loga,logk) = reshape(funcs["Pert_uab"](t,loga,logk,[sol(t)[1:2*d];pvals]...),2*d,2*d)
    hp(t) = reshape(funcs["Pert_curv"]([sol(t)[1:2*d];pvals]...),d,d)
    #jac = funcs["Pert_jac"]
    #function Peom_jac(GJ,G,logk,t)
    #    GJ[:,:] = reshape(jac(t,logainit,logk,G...,[sol(t)[1:2*d];pvals]...),(2*d)^2,(2*d)^2)
    #    nothing
    #end
    function Gprime(dG,G,logk,t)
        @assert(1.0+1e-8 >= Eh(t) >= 0.0, "ϵH weird value: N=$t, ϵH=$(Eh(t))")
        Gπ = hp(t)
        G11 = @view(G[1:d,:])
        G12 = @view(G[d+1:end,:])
        #G21 = @view(G[:,1:d])
        #G22 = @view(G[:,d+1:end])
        dG .= uab(t,logainit,logk)*G .- vcat(Gπ*G11, Gπ*G12)# .+ hcat(G21*Gπ, G22*Gπ)
        #dG[:,:] = reshape(eom(t,logainit,logk,G...,[sol(t)[1:2*d];pvals]...),2*d,2*d)
        nothing
    end
    Nend = sol.t[end]
    Nhc = sol.t[end] - Nexit_to_end
    Gs = Function[]
    odefunc = ODEFunction(Gprime)#,jac=Peom_jac)
    SolT = eltype(sol(0.0))
    for (logk,Nexit) in zip(logks,Nexits)
        if verbose
            println("Solving logk $logk")
        end
        Nstart = Nexit - Nsubhorizon
        prob = ODEProblem(odefunc,Matrix{SolT}(I,2*d,2*d),(SolT(Nstart),SolT(Nend)),logk)
        psol = DifferentialEquations.solve(prob,options.scalars_solver,
        reltol=options.scalars_rtol,abstol=options.scalars_atol,
        maxiters=options.scalars_max_steps,save_everystep=options.save_everystep
        ,saveat=options.saveat)# automatically choose algorithm?
        G(N) = psol(N)
        push!(Gs,G)
    end
    if verbose
        println("Tracing out powerspectra")
    end
    Pi(N) = sol(N)[d+1:2*d]
    gNe(N) = vcat(-reshape(g(N),d,d)*Pi(N)/(2*Eh(N)),zeros((d)))
    estar = div(length(Nexits),2)+1
    Nstar = Nexits[estar]
    τ = -1/(exp(Nstar+logainit)*H(Nstar))
    Σ_BD = Array{SolT}(undef,length(logks),2*d,2*d)
    for (ek,logk) in enumerate(logks)
        Nstart = Nexits[ek] - Nsubhorizon
        kτ = -exp(logk - Nstart - logainit)/H(Nstart)
        gi = reshape(ginv(Nstart),d,d)
        S = 1/2 * H(Nstart)^2 * abs(kτ)^2 *
        [gi -gi;
         -gi gi.*abs(kτ)^2]
        Σ_BD[ek,:,:] = S
    end
    function Σ(ek,N)
        G = convert(Array{BigFloat},Gs[ek](N))
        S = convert(Array{BigFloat},Σ_BD[ek,:,:])
        @assert !any(isnan,S)
        @assert !any(isinf,S)
        @assert !any(isnan,G)
        @assert !any(isinf,G)
        sig = G*S*G'
        #sig = Array{Float64}(undef,2*d,2*d)
        #for a in 1:2*d
        #    for b in 1:2*d
        #        sig[a,b] = 0.0
        #        for c in 1:2*d
        #            for d in 1:2*d
        #                sig[a,b] += G[a,c]*G[b,d]*S[c,d]
        #                if isnan(sig[a,b])
        #                    println("$(G[a,c]), $(G[b,d]), $(S[c,d])")
        #                end
        #            end
        #        end
        #    end
        #end
        @assert !any(isnan,sig)
        return sig
    end
    function Pζ(ek,N)
        gN = convert(Array{BigFloat},gNe(N))
        @assert !any(isnan,gN)
        #if any([x<0 for x in gN])
        #    println(ek,N)
        #end
        #@assert(all([x>0 for x in gN]),"negative in gN")
        S = convert(Array{BigFloat},Σ(ek,N))
        @assert !any(isnan,S)
        #@assert(all([x>0 for x in S]),"negative in Σ")
        #if ek == estar
        #    println(Gs[ek](N))
        #    println(gN)
        #    println(S)
        #end
        P = gN'*S*gN
        if isnan(P)
            println(gN)
            println(S)
        end
        @assert !isnan(P)
        #P = 0.0
        #for a in 1:2*d
        #    for b in 1:2*d
        #        P += gN[a]*gN[b]*S[a,b]
        #    end
        #end
        P /= (2*pi^2)
        #@assert(P>0,"negative P")
        return P
    end

    # adapted for curved space -- TODO need a unit test
    function Piso(ek,N)
        gN = gNe(N)
        gi = reshape(ginv(N),d,d)
        viα = LinearAlgebra.nullspace(convert(Array{SolT,2},(gi*gN[1:d])'))
        S = Σ(ek,N)
        P = 0.0
        for a in 1:d
            for b in 1:d
                for i in 1:(d-1)
                    P += viα[a,i]*S[a,b]*viα[b,i]
                end
            end
        end
        P *= 1/(2*Eh(N))
        P /= (2*pi^2)
        #if verbose
            #@assert(P>0,"negative P")
        #end
        return P
    end

    Pζs = BigFloat[]
    #Pζs_hc = Float64[]
    Pisos = BigFloat[]
    for (ek,logk) in enumerate(logks)
        append!(Pζs,Pζ(ek,Nend))
        if d>1
            append!(Pisos,Piso(ek,Nend))
        end
        #append!(Pζs_hc,Pζ(ek,Nhc+3.0))
    end

    # fit for ns, as
    # logPζ = As + (ns-1) log(k/k⋆) + 1/2 αs log(k/k⋆)^2 + ...
    if verbose
        println("Fitting powerspectra")
    end
    p = poly_fit(logks.-logks[estar],log.(Pζs),2)
    As = exp(p[1])
    ns = 1+p[2]
    as = p[3]*2
    if verbose
        println("As = $As")
        println("ns = $ns")
        println("αs = $as")
    end
    #p = poly_fit(logks.-logks[estar],log.(Pζs_hc),2)
    #Ashc = exp(p[1])
    #nshc = 1+p[2]
    #ashc = p[3]*2
    #if verbose
    #    println("Ashc = $Ashc")
    #    println("nshc = $nshc")
    #    println("αshc = $ashc")
    #end
    #println(errors)
    if options.compute_isocurvature && d > 1
        p = poly_fit(logks.-logks[estar],log.(Pisos),2)
        Asiso = exp(p[1])
        nsiso = 1+p[2]
        asiso = p[3]*2
        if verbose
            println("Asiso = $Asiso")
            println("nsiso = $nsiso")
            println("αsiso = $asiso")
        end
    end

    Ns = (Nend-Nexit_to_end-Nsubhorizon):0.1:Nend
    rs = Array{BigFloat}(undef,length(Ns))
    if options.compute_tensors
        global rs
        if verbose
            println("Computing tensor perturbations at kstar")
        end
        global r
        r = NaN
        tensor_eom = funcs["Tensor_eom"]
        #tensor_jac = funcs["Tensor_jac"]

        #function Gamjac(J,Gam,logk,t)
        #    J[:,:] = reshape(tensor_jac(t,logainit,logk,Gam[1,1],Gam[1,2],Gam[2,2],[sol(t)[1:2*d];pvals]...),4,4)
        #end
        function Gamprime(dGam,Gam,logk,t)
            dGam[:,:] = reshape(tensor_eom(t,logainit,logk,Gam[1,1],Gam[1,2],Gam[2,2],[sol(t)[1:2*d];pvals]...),2,2)
        end
        Nexit = Nexits[estar]
        logk = logks[estar]
        if verbose
            println("Solving logk $logk")
        end
        Nstart = Nexit - Nsubhorizon
        kτH = -exp(logk - Nstart - logainit)
        kτ = kτH/H(Nstart)
        GamIC = [ 1 -1 ; -1  abs(kτ)^2] * abs(kτH)^2
        odefunc = ODEFunction(Gamprime)#,jac=Gamjac)
        tprob = ODEProblem(odefunc,GamIC,(Nstart,Nend),logk)
        tsol = DifferentialEquations.solve(tprob,options.tensors_solver,
        reltol=options.tensors_rtol,abstol=options.tensors_atol,
        maxiters=options.tensors_max_steps,save_everystep=options.save_everystep,
        saveat=options.saveat)# automatically choose algorithm?
        Gam(N) = tsol(N)
        #push!(Gs,G)
        r = 4*Gam(Nend)[1,1]/As / (2*pi^2)
        if verbose
            println("Tensor to scalar ratio: $r")
        end
        rf(N) = 4*Gam(N)[1,1]/As / (2*pi^2)
        rs = rf.(Ns)
    else
        global rs
        global r
        r = NaN
        rs = []
    end

    PζsN = [Pζ(estar,N) for N in Ns]
    PisosN = [Piso(estar,N) for N in Ns]

    tsol = Dict()
    tsol["As"] = As
    tsol["ns"] = ns
    tsol["as"] = as
    tsol["r"] = r
    tsol["PzN"]=PζsN
    tsol["PisoN"]=PisosN
    tsol["rN"] = rs
    tsol["Pzk"]=Pζs
    tsol["Pisok"]=Pisos
    tsol["lnk"] = logks.-logks[estar]
    tsol["pert_Ne"]=Ns

    return tsol
end


text:
using OhMyJulia
using DataFrames
using PyCall
using NaturalSort
using Fire
using JSON
@pyimport gurobipy as grb

include("path.jl")
include("weight.jl")

#const data_list = ["abilene", "gscale", "Cernet", "Cesnet201006", "Chinanet", "Esnet", "Garr200112",
#                   "Grnet", "Ibm", "jelly16", "Nordu1997", "Sprint", "Tinet"]
const data_list = ["Cernet" , "Globalcenter" , "Goodnet" , "Gridnet" , "gscale" , "Janetbackbone" ,
"Rediris" , "Renater2010" , "Sinet" , "Sprint" , "SprintINET" , "Surfnet" ,
"SwitchL3" , "TataNld" , "Tinet" , "Uninett2011" , "Uunet" , "VtlWavenet2011" , "Xeex" , "Xspedius"]

const algo_list = ["Ecmp", "Edksp", "Ksp", "Custom", "Raeke", "Vlb"]

function read_topo(file)
    nodes, edges = Set{String}(), String[]
    for line in eachline(file) @when (m = match(r"(s\d+) -> (s\d+)", line)) != nothing
        push!(nodes, m[1], m[2])
        push!(edges, "$(m[1]) $(m[2])")
    end
    nodes, edges
end

function read_demand(hosts, demands)
    results = []
    hosts = readlines(hosts)
    n = length(hosts)
    for line in eachline(demands)
        line = map(x->parse(f64, x), split(line, ' '))
        demand = Dict{String, Float64}()
        for i in 1:n, j in 1:n @when i != j
            demand["$(hosts[i]) $(hosts[j])"] = line[(i-1) * n + j] / 2^20 # to Megabyte
        end
        push!(results, demand)
    end
    results
end

function read_budget(hosts, budgets)
    Dict(zip(map(x->replace(x, 'h', 's'), readlines(hosts)),
             map(x->parse(Int, x),        readlines(budgets))))
end

function write_scheme(io, scheme, x=Dict())
    for (pair, paths) in scheme
        println(io, pair, ':')
        x[pair] = xx = Dict()
        for (path, w) in paths
            p = map(x->split(x, ' '), path)
            println(io, join(map(car, p) ++ cadr(p[end]), " -> "), " @ ", w)
            xx[join(map(car, p) ++ cadr(p[end]), " -> ")] = w
        end
        println(io)
    end
end

function optimal_mcf(nodes, edges, demand)
    tpair(pair) = split(replace(pair, 'h', 's'), ' ')

    le = Dict(e => Dict() for e in edges)

    conn = Dict(n => (
        [s for s in nodes if "$s $n" in edges],
        [t for t in nodes if "$n $t" in edges]
    ) for n in nodes)

    m = grb.Model("optimal_mcf")
    m[:setParam]("OutputFlag", false)
    Z = m[:addVar](name="Z")

    for (pair, d) in demand
        u, v = tpair(pair)
        luv = Dict(e => (le[e][pair] = m[:addVar](name="l_$(e)_$(u)_$(v)")) for e in edges)
        m[:update]()

        # constraint 1: influx == outflux for core nodes
        for n in nodes @when n ∉ (u, v)
            ti = [luv["$u $n"] for u in car(conn[n])]
            to = [luv["$n $v"] for v in cadr(conn[n])]

            m[:addConstr](py"sum($ti) == sum($to)")
        end

        # constraint 2: demands are fulfilled
        ti = [luv["$x $v"] for x in car(conn[v])]
        to = [luv["$v $x"] for x in cadr(conn[v])]
        m[:addConstr](py"sum($ti) - sum($to) == $d")

        ti = [luv["$x $u"] for x in car(conn[u])]
        to = [luv["$u $x"] for x in cadr(conn[u])]
        m[:addConstr](py"sum($to) - sum($ti) == $d")
    end

    # constraint 3: minimize maximum link utilization
    for (edge, flows) in le
        flows = collect(values(flows))
        m[:addConstr](py"sum($flows) <= 1000 * $Z")
    end

    m[:setObjective](Z)
    m[:optimize]()

    json = Dict("Z" => m[:getObjective]()[:getValue]())

    for (edge, flows) in le, (pair, flow) in flows
        json["$pair $edge"] = flow[:X]
    end

    json
end

function single_failure(nodes, edges, demand, scheme)
    rev(x) = join(reverse(split(x, ' ')), ' ')

    stat(edges, scheme) = begin
        m1_scheme, Z = minimize_maximum_link_utilization(nodes, edges, demand, scheme)
        m2_scheme, Zs = minimize_maximum_link_utilization_then_maximize_throuput(nodes, edges, demand, scheme)
        Z, sum(d / Zs[pair] for (pair, d) in demand)
    end

    json = Dict()

    for e in edges @when (<)(split(e, ' ')...)
        edges_f = filter(x->x ∉ (e, rev(e)), edges)
        scheme_f = Dict()
        for (pair, paths) in scheme
            paths_f = [(path, weight) for (path, weight) in paths if e ∉ path && rev(e) ∉ path]
            # no need to renormalize since we will overwrite them anyway
            isempty(paths_f) && @goto disconnected
            scheme_f[pair] = paths_f
        end

        json[e] = stat(edges_f, scheme_f)
        continue

        @label disconnected
        json[e] = "some pairs have no path"
    end

    json
end

@main function main(data; optimal::Bool=false)
    @assert data in data_list
    io, process = open(`yates -budget 16 $data`)
    path_sets = parse_schemes(io)
    nodes, edges = read_topo("data/$data.dot")
    demands = read_demand("data/$data.hosts", "data/$data.demands")
    budgets = read_budget("data/$data.hosts", "data/$data.budgets")

    for (i, demand) in enumerate(demands)
        if optimal
            open("results/$data-$i-optimal-mcf.json", "w") do fout
                JSON.print(fout, optimal_mcf(nodes, edges, demand), 2)
            end
            continue
        end

        for algo in algo_list, (name, select) in (("program", select_program),
                                                  ("greedy", select_greedy),
                                                  ("hardnop", select_hard_nop)) @when length(nodes) < 16 || name != "greedy"
            json = Dict()

            open("results/$data-$i-$algo-$name.result", "w") do fout
                tic()
                raw_scheme = select(nodes, edges, path_sets[algo], budgets)
                time1 = toq()

                any(isempty, values(raw_scheme)) && return println(fout, "cannot solve")

                tic()
                m1_scheme, Z = minimize_maximum_link_utilization(nodes, edges, demand, raw_scheme)
                time2 = toq()

                tic()
                m2_scheme, Zs = minimize_maximum_link_utilization_then_maximize_throuput(nodes, edges, demand, raw_scheme)
                time3 = toq()

                println(fout, "number of nodes (switch): ", length(nodes))
                json["number of nodes (switch)"] = length(nodes)

                println(fout, "number of edges (core): ", length(edges))
                json["number of edges (core)"] = length(edges)

                println(fout, "number of demand pairs: ", length(demand))
                json["number of demand pairs"] = length(demand)

                println(fout, "average budgets: ", mean(values(budgets)))
                json["average budgets"] = mean(values(budgets))

                println(fout)

                println(fout, "time elapsed during path selection: ", time1, 's')
                json["time elapsed during path selection"] = time1

                println(fout, "time elapsed during first step: ", time2, 's')
                json["time elapsed during first step"] = time2

                println(fout, "time elapsed during first + second step: ", time3, 's')
                json["time elapsed during first + second step"] = time3

                println(fout)

                println(fout, "Z: ", Z)
                json["Z"] = Z

                println(fout, "total throuput before the second step:", sum(values(demand)) / max(Z, 1))
                json["total throuput before the second step"] = sum(values(demand)) / max(Z, 1)

                total_throuput = sum(d / Zs[pair] for (pair, d) in demand)
                println(fout, "total throuput after the second step:", total_throuput)
                json["total throuput after the second step"] = total_throuput

                println(fout)

                println(fout, "number of path through each node:")
                json["number of path through each node"] = x = Dict()

                node_dict = Dict(n => 0 for n in nodes)
                for (pair, paths) in m2_scheme, (path, weight) in paths, edge in path
                    n = split(edge, ' ') |> car
                    if car(n) == 's'
                        node_dict[n] += 1
                    end
                end
                for (node, number) in sort(collect(node_dict), lt=natural, by=car)
                    println(fout, "  ", node, ": ", number)
                    x[node] = number
                end

                println(fout, "  average: ", sum(values(node_dict)) / length(node_dict))
                x["average"] = sum(values(node_dict)) / length(node_dict)

                println(fout)

                println(fout, "number of path, total flow, and link utiliaztion of each edge:")
                json["number of path, total flow, and link utiliaztion of each edge"] = x = Dict()

                edge_dict = Dict(e => f64[] for e in edges)
                for (pair, paths) in m2_scheme, (path, weight) in paths, edge in path @when 'h' ∉ edge
                    push!(edge_dict[edge], demand[pair] * weight / Zs[pair])
                end
                for (edge, flows) in sort(collect(edge_dict), lt=natural, by=car)
                    println(fout, "  ", edge, ": ", length(flows), ", ", sum(flows), ", ", sum(flows) / 10, '%')
                    x[edge] = [length(flows), sum(flows), sum(flows) / 10]
                end

                println(fout, "  average: ", sum(length(v) for v in values(edge_dict)) / length(edge_dict),
                                       ", ", sum(sum(v) for v in values(edge_dict)) / length(edge_dict),
                                       ", ", sum(sum(v) for v in values(edge_dict)) / length(edge_dict) / 10, '%')
                x["average"] = [sum(length(v) for v in values(edge_dict)) / length(edge_dict),
                                sum(sum(v) for v in values(edge_dict)) / length(edge_dict),
                                sum(sum(v) for v in values(edge_dict)) / length(edge_dict) / 10]

                println(fout)

                println(fout, "number of path and demand ratio of each st pair:")
                json["number of path and demand ratio of each st pair"] = x = Dict()

                pair_dict = Dict(pair => length(paths) for (pair, paths) in m2_scheme)
                for (pair, nflows) in sort(collect(pair_dict), lt=natural, by=car)
                    println(fout, "  ", pair, ": ", nflows, ", ", 1 / Zs[pair] * 100, '%')
                    x[pair] = [nflows, 1 / Zs[pair] * 100]
                end
                println(fout, "  average: ", sum(values(pair_dict)) / length(pair_dict),
                                       ", ", sum(d / Zs[pair] for (pair, d) in demand) / sum(values(demand)) * 100, '%')
                x["average"] = [sum(values(pair_dict)) / length(pair_dict),
                                sum(d / Zs[pair] for (pair, d) in demand) / sum(values(demand)) * 100]

                println(fout)

                println(fout, "Z, total throuput, and throuput ratio at single edge failure:")
                json["Z, total throuput, and throuput ratio at single edge failure"] = x = Dict()

                edge_dict = single_failure(nodes, edges, demand, raw_scheme)
                for (edge, res) in sort(collect(edge_dict), lt=natural, by=car)
                    if res isa String
                        println(fout, "  ", edge, ": ", res)
                        x[edge] = res
                    else
                        Z, tt = res
                        println(fout, "  ", edge, ": ", Z, ", ", tt, ", ", tt / total_throuput * 100, '%')
                        x[edge] = Z, tt, tt / total_throuput * 100
                    end
                end
                for r in (90, 95, 99, 99.9, 99.99)
                    println(fout, "  P(T > $r%): ", mean(p isa Tuple && p[3] >= r for (e, p) in x if !startswith(e, 'P')))
                    x["P(T > $r%)"] = mean(p isa Tuple && p[3] >= r for (e, p) in x if !startswith(e, 'P'))
                end
                println(fout)

                println(fout, "path set:")
                json["path set"] = x = Dict()
                write_scheme(fout, path_sets[algo], x)
                println(fout)

                println(fout, "paths selected:")
                json["paths selected"] = x = Dict()
                write_scheme(fout, raw_scheme, x)
                println(fout)

                println(fout, "paths after the first step (minimize maximum link utilization):")
                json["paths after the first step (minimize maximum link utilization)"] = x = Dict()
                write_scheme(fout, m1_scheme, x)
                println(fout)

                println(fout, "paths after the second step (maximize total throuput):")
                json["paths after the second step (maximize total throuput)"] = x = Dict()
                write_scheme(fout, m2_scheme, x)
            end

            open("results/$data-$i-$algo-$name.json", "w") do fout
                JSON.print(fout, json, 2)
            end
        end
    end
end

ext: jl

max_stars_repo_name: ylxdzsw/TEexp

max_issues_repo_name: ylxdzsw/TEexp

max_forks_repo_name: ylxdzsw/TEexp


text:
module StdMsgs
include("std_msgs.jl")
export std_msgs
end


text:
module ReproducePlotUtils

using Reproduce, Plots, RollingFunctions, Statistics, FileIO

include("stat_utils.jl")

function get_agg(agg, ddict, key)
    agg(ddict["results"][key])
end

get_MEAN(ddict, key) = get_agg(mean, ddict, key)
get_AUC(ddict, key) = get_agg(sum, ddict, key)
get_AUE(ddict, key, perc=0.1) = get_agg(ddict, key) do x
    sum(x[end-max(1, Int(floor(length(x)*perc))):end])
end

get_MUE(ddict, key, perc=0.1) = get_agg(ddict, key) do x
    mean(x[end-max(1, Int(floor(length(x)*perc))):end])
end

function get_rolling_mean_line(ddict, key, n)
    if n > length(ddict["results"][key])
        n = length(ddict["results"][key])
    end
    rollmean(ddict["results"][key], n)
end

function get_extended_line(ddict, key1, key2; n=0)
    ret = zeros(eltype(ddict["results"][key1]), sum(ddict["results"][key2]))
    cur_idx = 1
    for i in 1:length(ddict["results"][key1])
        ret[cur_idx:(cur_idx + ddict["results"][key2][i] - 1)] .= ddict["results"][key1][i]
        cur_idx += ddict["results"][key2][i]
    end

    if n == 0
        ret
    else
        rollmean(ret, n)[1:n:end]
    end
end

function load_data(loc)
    ic = ItemCollection(loc)
    ic, diff(ic)
end

include("data.jl")
include("recipes.jl")

end


text:
# Generate scales.

import JSON

# This ID is incremented for every problem generated.  Barring changes
# to the generation, this should generate the same ids on subsequent
# runs.
cur_id = 1

function emit(question, answer)
    global cur_id
    println("INSERT INTO \"probs\" VALUES($cur_id, '$question', '$answer');")
    cur_id += 1
end

key_names = [ "C", "G", "D", "A", "E", "B", "F♯", "G♭", "D♭", "A♭", "E♭", "B♭", "F" ]

function gen_scale(ival, hands, name, hand; style="updown", octaves=2)
    for key in key_names
        answer = Dict("type" => "scale",
            "base" => key,
            "intervals" => ival,
            "hands" => hands,
            "octaves" => octaves,
            "style" => style)
        emit("$hand-scale $key $name", JSON.json(answer))
    end
end

function gen_scale(ival, name; style="updown", octaves=2)
    gen_scale(ival, 1, name, "RH", style=style, octaves=octaves)
    gen_scale(ival, 1, name, "LH", style=style, octaves=octaves)
    gen_scale(ival, 2, name, "2H", style=style, octaves=octaves)
end

println("DELETE FROM probs;")
# Learn major scales with both hands.
gen_scale("WWHWWWH", "major")

# The rest of the scales, just right hand.
gen_scale("WHWWWHW", 2, "minor (dorian)", "2H")
gen_scale("WWHWWHW", 2, "dominant (mixolydian)", "2H")
gen_scale("HWWHWWW", 2, "half diminished (locrian)", "2H")
#gen_scale("WHWHWHWH", 2, "diminished (whole-half)", "2H")
gen_scale("WHWHWHWH", "diminished (whole-half)")
gen_scale("HWHWHWHW", 2, "sym-dom (half-whole)", "2H")

# Simple patterns on major.
gen_scale("WWHWWWH", 1, "major 3rds", "RH", style="3up",
    octaves=1)
gen_scale("WWHWWWH", 1, "major 3rds", "LH", style="3up",
    octaves=1)
gen_scale("WWHWWWH", 2, "major 3rds", "2H", style="3up",
    octaves=1)
gen_scale("WWHWWWH", 1, "major 3rds rev", "RH", style="3upr",
    octaves=1)
gen_scale("WWHWWWH", 1, "major 3rds rev", "LH", style="3upr",
    octaves=1)
gen_scale("WWHWWWH", 2, "major 3rds rev", "2H", style="3upr",
    octaves=1)


text:
using Banana
using Test

# ll = LinkedList(["apple", "pear", "banana"])
#
# state = start(ll)
# item, state = next(ll, state)
# insert!(ll, state, "kiwi")
# for item in ll println(item) end

@test true

# @test first(ll) == "kiwi"

# end


text:
module test_bukdu_not_found

using Test
using Bukdu

result = Router.call(get, "/Lorem ipsum dolor sit amet consectetur")
@test result.route.action === Bukdu.System.not_found

result = Router.call(get, "/Lorem ipsum                           ")
@test result.route.action === Bukdu.System.not_found

result = Router.call(get, "/오우놀줄아는놈인가")
@test result.route.action === Bukdu.System.not_found

end # module test_bukdu_not_found


text:
function Generator()
    return Generator(())
end

function add_propagator(obj::Generator, arg0::Propagator)
    return jcall(obj, "addPropagator", ObservableSatellite, (Propagator,), arg0)
end

function add_scheduler(obj::Generator, arg0::Scheduler)
    return jcall(obj, "addScheduler", void, (Scheduler,), arg0)
end

function generate(obj::Generator, arg0::AbsoluteDate, arg1::AbsoluteDate)
    return jcall(obj, "generate", SortedSet, (AbsoluteDate, AbsoluteDate), arg0, arg1)
end

function get_propagator(obj::Generator, arg0::ObservableSatellite)
    return jcall(obj, "getPropagator", Propagator, (ObservableSatellite,), arg0)
end


text:
__precompile__()

module PyProj

import PyCall

const pyproj = PyCall.PyNULL()
const Proj = PyCall.PyNULL()
const Geod = PyCall.PyNULL()
const _transform = PyCall.PyNULL()

function __init__()
    copy!(pyproj, PyCall.pyimport_conda("pyproj", "pyproj"))

    copy!(Proj, pyproj[:Proj])
    copy!(Geod, pyproj[:Geod])
    copy!(_transform, pyproj[:transform])
end

function transform(p1, p2, x, y, z, radians=false)
    _transform(p1, p2, x, y, z, radians)
end

end # module


text:
using GLVisualize, GeometryTypes, FileIO
using GLAbstraction, Colors, Reactive
if !isdefined(:runtests)
	window = glscreen()
	timesignal = bounce(0f0:0.1:1f0)
end

primitive = SimpleRectangle(0f0,-0.5f0,1f0,1f0)
positions = rand(10f0:0.01f0:200f0, 10)

function interpolate(a, positions, t)
    [ae+((be-ae)*t) for (ae, be) in zip(a,positions)]
end
t = const_lift(*, timesignal, 10f0)
interpolated = foldp((positions,positions,positions), t) do v0_v1_ip, td
    v0,v1,ip = v0_v1_ip
    pol = td%1
    if isapprox(pol, 0.0)
        v0 = v1
        v1 = map(x-> rand(linspace(-50f0, 60f0, 100)), v0)
    end
    v0, v1, interpolate(v0, v1, pol)
end
b_sig = map(last, interpolated)
bars = visualize(
    (RECTANGLE, b_sig),
    intensity=b_sig,
    ranges=linspace(0,600, 10),
    color_norm=Vec2f0(-40,200),
    color_map=GLVisualize.default(Vector{RGBA})
)
_view(bars, window, camera=:orthographic_pixel)

if !isdefined(:runtests)
	renderloop(window)
end


text:
"""
    partition(rows::AbstractVector{Int}, fractions...; shuffle=false)

Splits the vector `rows` into a tuple of vectors whose lengths are
given by the corresponding `fractions` of `length(rows)`. The last
fraction is not provided, as it is inferred from the preceding
ones. So, for example,

    julia> partition(1:1000, 0.2, 0.7)
    (1:200, 201:900, 901:1000)

"""
function partition(rows::AbstractVector{Int}, fractions...; shuffle::Bool=false)
    rows = collect(rows)
    shuffle && shuffle!(rows)
    rowss = []
    if sum(fractions) >= 1
        throw(DomainError)
    end
    n_patterns = length(rows)
    first = 1
    for p in fractions
        n = round(Int, p*n_patterns)
        n == 0 ? (@warn "A split has only one element"; n = 1) : nothing
        push!(rowss, rows[first:(first + n - 1)])
        first = first + n
    end
    if first > n_patterns
        @warn "Last vector in the split has only one element."
        first = n_patterns
    end
    push!(rowss, rows[first:n_patterns])
    return tuple(rowss...)
end

function finaltypes(T::Type)
    s = InteractiveUtils.subtypes(T)
    if isempty(s)
        return [T, ]
    else
        return reduce(vcat, [finaltypes(S) for S in s])
    end
end

import Random
import Base.length

function unique_inverse(A::AbstractArray)
    out = Array{eltype(A)}(undef, 0)
    out_idx = Array{Vector{Int}}(undef, 0)
    seen = Dict{eltype(A), Int}()
    for (idx, x) in enumerate(A)
        if !in(x, keys(seen))
            seen[x] = length(seen) + 1
            push!(out, x)
            push!(out_idx, Int[])
        end
        push!(out_idx[seen[x]], idx)
    end
    out, out_idx
end

"""
    StratifiedKFold(strata,k)
Struct for StratifiedKFold
provide strata's and number of partitions(k) and simply collect the object for the indices.
Taken from MLBase (https://github.com/JuliaStats/MLBase.jl).
"""
struct StratifiedKFold
    n::Int                         #Total number of observations
    permseqs::Vector{Vector{Int}}  #Vectors of vectors of indexes for each stratum
    k::Int                         #Number of splits
    coeffs::Vector{Float64}        #About how many observations per strata are in a val set
    function StratifiedKFold(strata, k)
        2 <= k <= length(strata) || error("The value of k must be in [2, length(strata)].")
        strata_labels, permseqs = unique_inverse(strata)
        map(Random.shuffle!, permseqs)
        coeffs = Float64[]
        for (stratum, permseq) in zip(strata_labels, permseqs)
            k <= length(permseq) || error("k is greater than the length of stratum $stratum")
            push!(coeffs, length(permseq) / k)
        end
        new(length(strata), permseqs, k, coeffs)
    end
end
length(c::StratifiedKFold) = c.k
function Base.iterate(c::StratifiedKFold, s::Int=1)
    (s > c.k) && return nothing
    r = Int[]
    for (permseq, coeff) in zip(c.permseqs, c.coeffs)
        a, b = round.(Integer, [s-1, s] .* coeff)
        append!(r, view(permseq, a+1:b))
    end
    return setdiff(1:c.n, r), s+1
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct EncryptionSettings <: SwaggerModel
    enabled::Any # spec type: Union{ Nothing, Bool } # spec name: enabled
    diskEncryptionKey::Any # spec type: Union{ Nothing, KeyVaultAndSecretReference } # spec name: diskEncryptionKey
    keyEncryptionKey::Any # spec type: Union{ Nothing, KeyVaultAndKeyReference } # spec name: keyEncryptionKey

    function EncryptionSettings(;enabled=nothing, diskEncryptionKey=nothing, keyEncryptionKey=nothing)
        o = new()
        validate_property(EncryptionSettings, Symbol("enabled"), enabled)
        setfield!(o, Symbol("enabled"), enabled)
        validate_property(EncryptionSettings, Symbol("diskEncryptionKey"), diskEncryptionKey)
        setfield!(o, Symbol("diskEncryptionKey"), diskEncryptionKey)
        validate_property(EncryptionSettings, Symbol("keyEncryptionKey"), keyEncryptionKey)
        setfield!(o, Symbol("keyEncryptionKey"), keyEncryptionKey)
        o
    end
end # type EncryptionSettings

const _property_map_EncryptionSettings = Dict{Symbol,Symbol}(Symbol("enabled")=>Symbol("enabled"), Symbol("diskEncryptionKey")=>Symbol("diskEncryptionKey"), Symbol("keyEncryptionKey")=>Symbol("keyEncryptionKey"))
const _property_types_EncryptionSettings = Dict{Symbol,String}(Symbol("enabled")=>"Bool", Symbol("diskEncryptionKey")=>"KeyVaultAndSecretReference", Symbol("keyEncryptionKey")=>"KeyVaultAndKeyReference")
Base.propertynames(::Type{ EncryptionSettings }) = collect(keys(_property_map_EncryptionSettings))
Swagger.property_type(::Type{ EncryptionSettings }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_EncryptionSettings[name]))}
Swagger.field_name(::Type{ EncryptionSettings }, property_name::Symbol) =  _property_map_EncryptionSettings[property_name]

function check_required(o::EncryptionSettings)
    true
end

function validate_property(::Type{ EncryptionSettings }, name::Symbol, val)
end


text:
function bob(stimulus)
    is_empty = all([isspace(c) for c in stimulus])

    if !is_empty
        is_question = last(
            [c for c in stimulus if !isspace(c)]) == '?'
    else
        is_question = false
    end

    is_allcaps = (
        (uppercase(stimulus) == stimulus) &
        any([isletter(c) for c in stimulus]))
    # All caps mean 1. there are letters
    # 2. letters are already uppercase

    if is_empty
        return "Fine. Be that way!"
    elseif is_question
        if is_allcaps
            return "Calm down, I know what I'm doing!"
        else
            return "Sure."
        end
    elseif is_allcaps
        return "Whoa, chill out!"
    else
        return "Whatever."
    end
end

max_stars_repo_name: jfcliang/jl-exercism-solutions

max_issues_repo_name: jfcliang/jl-exercism-solutions

max_forks_repo_name: jfcliang/jl-exercism-solutions


text:
export getNumSubscribers, getTopic, publish, shutdown

"""
    getNumSubscribers(publisher)

Returns the number of subscribers that are currently connected to this Publisher.
"""
getNumSubscribers(publisher::Union{Cxx.CxxCore.CppPtr{Cxx.CxxCore.CppValue{Cxx.CxxCore.CxxQualType{Cxx.CxxCore.CppBaseType{Symbol("ros::Publisher")},(false, false, false)},N} where N,(false, false, false)}, Cxx.CxxCore.CppValue{Cxx.CxxCore.CxxQualType{Cxx.CxxCore.CppBaseType{Symbol("ros::Publisher")},(false, false, false)},N} where N}) = @cxx publisher->getNumSubscribers()

"""
    getTopic(publisher)

Returns the topic that this Publisher will publish on.
"""
getTopic(publisher::Union{Cxx.CxxCore.CppPtr{Cxx.CxxCore.CppValue{Cxx.CxxCore.CxxQualType{Cxx.CxxCore.CppBaseType{Symbol("ros::Publisher")},(false, false, false)},N} where N,(false, false, false)}, Cxx.CxxCore.CppValue{Cxx.CxxCore.CxxQualType{Cxx.CxxCore.CppBaseType{Symbol("ros::Publisher")},(false, false, false)},N} where N}) = unsafe_string(@cxx publisher->getTopic())

"""
    isLatched(publisher)

Returns whether or not this topic is latched.
"""
isLatched(publisher::Union{Cxx.CxxCore.CppPtr{Cxx.CxxCore.CppValue{Cxx.CxxCore.CxxQualType{Cxx.CxxCore.CppBaseType{Symbol("ros::Publisher")},(false, false, false)},N} where N,(false, false, false)}, Cxx.CxxCore.CppValue{Cxx.CxxCore.CxxQualType{Cxx.CxxCore.CppBaseType{Symbol("ros::Publisher")},(false, false, false)},N} where N}) = @cxx publisher->isLatched()

"""
    publish(publisher, msg)

Publish a message on the topic associated with this Publisher.
"""
publish(publisher::Union{Cxx.CxxCore.CppPtr{Cxx.CxxCore.CppValue{Cxx.CxxCore.CxxQualType{Cxx.CxxCore.CppBaseType{Symbol("ros::Publisher")},(false, false, false)},N} where N,(false, false, false)}, Cxx.CxxCore.CppValue{Cxx.CxxCore.CxxQualType{Cxx.CxxCore.CppBaseType{Symbol("ros::Publisher")},(false, false, false)},N} where N}, msg::Cxx.CxxCore.CppValue) = @cxx publisher->publish(msg)

publish(publisher::Union{Cxx.CxxCore.CppPtr{Cxx.CxxCore.CppValue{Cxx.CxxCore.CxxQualType{Cxx.CxxCore.CppBaseType{Symbol("ros::Publisher")},(false, false, false)},N} where N,(false, false, false)}, Cxx.CxxCore.CppValue{Cxx.CxxCore.CxxQualType{Cxx.CxxCore.CppBaseType{Symbol("ros::Publisher")},(false, false, false)},N} where N}, msg::Cxx.CxxCore.CppPtr) = @cxx publisher->publish(*(msg))

"""
    shutdown(publisher)

Shutdown the advertisement associated with this Publisher.
"""
shutdown(p::Union{Cxx.CxxCore.CppPtr{Cxx.CxxCore.CppValue{Cxx.CxxCore.CxxQualType{Cxx.CxxCore.CppBaseType{Symbol("ros::Publisher")},(false, false, false)},N} where N,(false, false, false)}, Cxx.CxxCore.CppValue{Cxx.CxxCore.CxxQualType{Cxx.CxxCore.CppBaseType{Symbol("ros::Publisher")},(false, false, false)},N} where N}) = @cxx p->shutdown()


text:
import Libdl
import Printf
using MATLAB

# include elements from the julia-config.jl script, needed for getting flags for
# building the mex file

isDebug() = ccall(:jl_is_debugbuild, Cint, ()) != 0

threadingOn() = ccall(:jl_threading_enabled, Cint, ()) != 0

function shell_escape(str)
    str = replace(str, "'" => "'\''")
    return "\"$str\""
end

function libDir()
    return if isDebug()
        dirname(abspath(Libdl.dlpath("libjulia-debug")))
    else
        dirname(abspath(Libdl.dlpath("libjulia")))
    end
end

private_libDir() = abspath(Sys.BINDIR, Base.PRIVATE_LIBDIR)

function includeDir()
    return abspath(Sys.BINDIR, Base.INCLUDEDIR, "julia")
end

function cflags()
    flags = IOBuffer()
    include = shell_escape(includeDir())
    if isDebug()
        print(flags, "-g")
    else
        print(flags, "-O")
    end
    print(flags, " -I", include)
    if threadingOn()
        print(flags, " -DJULIA_ENABLE_THREADING")
    end
    return String(take!(flags))
end

function ldflags()
    flags = IOBuffer()
    print(flags, "-L$(shell_escape(libDir()))")
    if Sys.isunix()
        print(flags, " -Wl,-rpath $(shell_escape(libDir()))")
    end
    return String(take!(flags))
end

function ldlibs()
    libname = if isDebug()
        "julia-debug"
    else
        "julia"
    end
    if Sys.isunix()
        return "-l$libname -ldl"
    else
        return realpath(joinpath(libDir(), "..", "lib", "lib$libname.dll.a"))
    end
end

# get build parameters
is_debug = isDebug()
threading_enabled = threadingOn()
julia_bin = is_debug ?
    joinpath(unsafe_string(Base.JLOptions().julia_bindir), "julia-debug") :
    joinpath(unsafe_string(Base.JLOptions().julia_bindir), "julia")
julia_home = unsafe_string(Base.JLOptions().julia_bindir)
sys_image = unsafe_string(Base.JLOptions().image_file)
lib_base = is_debug ? "julia-debug" : "julia"
lib_path = Libdl.dlpath("lib$lib_base")
lib_dir = Sys.iswindows() ? joinpath(julia_home, "..", "lib") : splitdir(lib_path)[1]
inc_dir = includeDir()
build_cflags = cflags()
build_ldflags = ldflags()
build_ldlibs = ldlibs()
build_src = abspath("mexjulia.cpp")
outdir = joinpath(pwd(),"..","mexjulia")
mex_cmd = "mex LDFLAGS=\'$(build_ldflags) \$LDFLAGS\' -v -largeArrayDims -outdir \"$outdir\" $(build_cflags) \"$(build_src)\" $(build_ldlibs)"

# Save build parameters to .mat file
mat"""
% save build parameters to .mat file
is_debug = $is_debug;
threading_enabled = $threading_enabled;
julia_bin = $julia_bin;
julia_home = $julia_home;
sys_image = $sys_image;
lib_base = $lib_base;
lib_path = $lib_path;
lib_dir = $lib_dir;
inc_dir = $inc_dir;
build_cflags = $build_cflags;
build_ldflags = $build_ldflags;
build_ldlibs = $build_ldlibs;
build_src = $build_src;
mex_cmd = $mex_cmd;
save("jldict.mat", "is_debug", "threading_enabled", "julia_bin", "julia_home",...
    "sys_image", "lib_base", "lib_path", "lib_dir", "inc_dir", "build_cflags",...
    "build_ldflags", "build_ldlibs", "build_src", "mex_cmd");
"""

# build mex file
s1 = MSession()
eval_string(s1, mex_cmd)
close(s1)

# move jldict.mat to mexjulia folder
mv("jldict.mat", joinpath(outdir, "jldict.mat"), force=true)

# add mexjulia directory to MATLAB path
mat"""
% check if current directory is on MATLAB path
path_dirs = regexp(path, pathsep, 'split');
if ispc
    on_path = any(strcmpi($outdir, path_dirs));
else
    on_path = any(strcmp($outdir, path_dirs));
end

% if not, add it and save
if ~on_path
    fprintf('%s is not on the MATLAB path. Adding it and saving...\\n\', $(pwd()));
    path($outdir, path);
    savepath;
end
"""


text:
"""
    ischordal(g)
return true if `g` is chordal, false otherwise.
"""
function ischordal(g)
    !is_connected(g) && return false
    ne(g) == nv(g) * (nv(g) - 1) / 2 && return true # if g is complete
    ordering = lbfs(g)
    for i in 3:length(ordering)
        index_neighbor = get_closest_neighbor(g, ordering, i)
        if index_neighbor == -1
            continue
        else
            prior_current = get_prior_neighbors(g, ordering, i, [index_neighbor])
            prior_closest = get_prior_neighbors(g, ordering, index_neighbor)
            if !issubset(prior_current, prior_closest)
                return false
            end
        end
    end
    return true
end

max_stars_repo_name: charlyalizadeh/PowerFlowNetworks

max_issues_repo_name: charlyalizadeh/PowerFlowNetworks

max_forks_repo_name: charlyalizadeh/PowerFlowNetworks


text:
using PyPlot
using Statistics
using PolynomialBasis
using ImplicitDomainQuadrature
using Revise
using CutCellDG
include("../../../../test/useful_routines.jl")
include("../transformation-elasticity-solver.jl")
TES = TransformationElasticitySolver

function perturbation(x, frequency, amplitude)
    return amplitude * sin.(2 * pi * frequency * x)
end

perturbed_distancefunction(x, initialposition, frequency, amplitude) =
    plane_distance_function(x, [1.0, 0.0], [initialposition, 0.0]) +
    perturbation(x[2,:], frequency, amplitude)

function plot_potential_components(
    ycoords,
    component,
    initialposition,
    frequency,
    amplitude,
    interfacescale,
    pdscale;
    filepath = "",
    ylabel = "",
)
    interfaceposition =
        initialposition .+ perturbation.(ycoords, frequency, amplitude)
    planeinterfacepd = mean(component)

    interfaceylim = interfacescale * amplitude
    pdylim = pdscale

    fig, ax = PyPlot.subplots(2, 1, sharex = true)
    ax[1].plot(ycoords, interfaceposition, color = "black")
    ax[1].set_ylim(
        initialposition - interfaceylim,
        initialposition + interfaceylim,
    )
    ax[1].grid()
    ax[1].set_ylabel("Interface position")

    ax[2].plot(ycoords, component, color = "black")
    ax[2].set_ylim(planeinterfacepd - pdylim, planeinterfacepd + pdylim)
    ax[2].grid()
    ax[2].set_ylabel(ylabel)
    ax[2].set_xlabel("y")
    fig.tight_layout()

    if length(filepath) > 0
        fig.savefig(filepath)
        return fig
    else
        return fig
    end
end

function plot_dilatation_and_normal_stress(
    ycoords,
    parentdilatation,
    productdilatation,
    parentsrr,
    productsrr;
    filepath = "",
)

    fig, ax = PyPlot.subplots(2, 1, sharex = true)
    ax[1].plot(ycoords, 1.0 .+ parentdilatation, label = "parent")
    ax[1].plot(ycoords, 1.0 .+ productdilatation, label = "product")
    ax[1].grid()
    ax[1].legend()
    ax[1].set_ylabel(L"1 + \epsilon_{kk}")

    ax[2].plot(ycoords, parentsrr, label = "parent")
    ax[2].plot(ycoords, productsrr, label = "product")
    ax[2].grid()
    ax[2].set_ylabel(L"\sigma_{nn}")
    ax[2].set_xlabel("y")

    if length(filepath) > 0
        fig.savefig(filepath)
        return fig
    else
        return fig
    end
end

function plot_normal_stress(
    ycoords,
    parentsrr,
    productsrr;
    filepath = "",
)

    fig, ax = PyPlot.subplots()
    ax.plot(ycoords, parentsrr, label = "parent")
    ax.plot(ycoords, productsrr, label = "product")
    ax.grid()
    ax.set_ylabel(L"\sigma_{nn}")
    ax.set_xlabel("y")
    ax.legend()

    if length(filepath) > 0
        fig.savefig(filepath)
        return fig
    else
        return fig
    end
end
distancefunction(x) =
    perturbed_distancefunction(x, initialposition, frequency, amplitude)

initialposition = 0.5
frequency = 2.5
amplitude = 1e-2
# amplitude = 1e-10
polyorder = 2
nelmts = 17
penaltyfactor = 1e3
meshwidth = [1.0, 1.0]
numqp = required_quadrature_order(polyorder) + 2

K1, K2 = 247.0, 247.0    # GPa
mu1, mu2 = 126.0, 126.0   # GPa

lambda1 = lame_lambda(K1, mu1)
lambda2 = lame_lambda(K2, mu2)
stiffness = CutCellDG.HookeStiffness(lambda1, mu1, lambda2, mu2)
rho1 = 3.93e3           # Kg/m^3
rho2 = 3.68e3           # Kg/m^3
V01 = 1.0 / rho1
V02 = 1.0 / rho2
ΔG0Jmol = -14351.0
molarmass = 0.147
ΔG0 = ΔG0Jmol / molarmass
theta0 = -0.067

numquerypoints = 1000
querypoints = vcat(
    repeat([initialposition], numquerypoints)',
    range(0, 1, length = numquerypoints)',
)
ycoords = querypoints[2, :]

transfstress =
    CutCellDG.plane_strain_transformation_stress(lambda1, mu1, theta0)

elasticitybasis = LagrangeTensorProductBasis(2, polyorder)
levelsetbasis = LagrangeTensorProductBasis(2,3)
dim, numpts = size(interpolation_points(levelsetbasis))
basispts = interpolation_points(elasticitybasis)

cgmesh = CutCellDG.CGMesh([0.0, 0.0], meshwidth, [nelmts, nelmts], numpts)
dgmesh = CutCellDG.DGMesh([0.0, 0.0], meshwidth, [nelmts, nelmts], basispts)
# CutCellDG.make_vertical_periodic!(dgmesh)

levelset = CutCellDG.LevelSet(distancefunction, cgmesh, levelsetbasis)

elementsize = CutCellDG.element_size(cgmesh)
minelmtsize = minimum(elementsize)
maxelmtsize = maximum(elementsize)
penalty = penaltyfactor / minelmtsize * 0.5 * (lambda1 + mu1 + lambda2 + mu2)
tol = minelmtsize^(polyorder + 1)
boundingradius = 1.5 * maxelmtsize

cutmesh = CutCellDG.CutMesh(dgmesh, levelset)
refseedpoints, seedcellids = CutCellDG.seed_zero_levelset(2, levelset, cutmesh)
spatialseedpoints =
    CutCellDG.map_to_spatial(refseedpoints, seedcellids, cutmesh)

mesh, cellquads, facequads, interfacequads =
    TES.construct_merged_mesh_and_quadratures(cutmesh, levelset, numqp)

nodaldisplacement = TES.nodal_displacement(
    mesh,
    elasticitybasis,
    cellquads,
    facequads,
    interfacequads,
    stiffness,
    theta0,
    penalty,
)

closestpoints, closestcellids = CutCellDG.closest_points_on_zero_levelset(
    querypoints,
    spatialseedpoints,
    seedcellids,
    levelset,
    tol,
    boundingradius,
)

normals = CutCellDG.collect_normals_at_spatial_points(
    closestpoints,
    closestcellids,
    levelset,
)

parentclosestrefpoints = CutCellDG.map_to_reference_on_merged_mesh(
    closestpoints,
    closestcellids,
    -1,
    mesh,
)
productclosestrefpoints = CutCellDG.map_to_reference_on_merged_mesh(
    closestpoints,
    closestcellids,
    +1,
    mesh,
)

################################################################################
parentstrain = CutCellDG.parent_strain(
    nodaldisplacement,
    elasticitybasis,
    parentclosestrefpoints,
    closestcellids,
    mesh,
)
parentstress = CutCellDG.parent_stress(parentstrain, stiffness)
parentstrainenergy = V02 * CutCellDG.strain_energy(parentstress, parentstrain)

parentradialtraction = CutCellDG.traction_force_at_points(parentstress, normals)
parentsrr = CutCellDG.traction_component(parentradialtraction, normals)
parentdilatation = CutCellDG.dilatation(parentstrain)

parentcompwork = V02 * (1.0 .+ parentdilatation) .* parentsrr
################################################################################

################################################################################
productstrain = CutCellDG.product_elastic_strain(
    nodaldisplacement,
    elasticitybasis,
    theta0,
    productclosestrefpoints,
    closestcellids,
    mesh,
)
productstress = CutCellDG.product_stress(productstrain, stiffness, theta0)
productstrainenergy =
    V01 * CutCellDG.strain_energy(productstress, productstrain)
productradialtraction =
    CutCellDG.traction_force_at_points(productstress, normals)
productsrr = CutCellDG.traction_component(productradialtraction, normals)
productdilatation = CutCellDG.dilatation(productstrain)
productcompwork = V01 * (1.0 .+ productdilatation) .* productsrr
################################################################################

srrmean = 0.5 * (parentsrr + productsrr)

jse = productstrainenergy - parentstrainenergy
jsediff = (maximum(jse) - minimum(jse)) / 2

# jcw = productcompwork - parentcompwork
jcw =
    (V01 * (1.0 .+ productdilatation) - V02 * (1.0 .+ parentdilatation)) .*
    srrmean
jcwdiff = (maximum(jcw) - minimum(jcw)) / 2

pd = jse - jcw
pddiff = (maximum(pd) - minimum(pd)) / 2

foldername = "examples\\time-step\\perturbed-plane-interface\\"*
             "hermite-levelset\\potential-components\\direct-srr-results\\"
interfacescale = 5
################################################################################

pdscale = 1.5pddiff
plot_potential_components(
    ycoords,
    pd,
    initialposition,
    frequency,
    amplitude,
    interfacescale,
    pdscale,
)

# plot_normal_stress(ycoords,parentsrr,productsrr)

# Error in [σnn]
# srrerr = productsrr - parentsrr
# srrerrdiff = (maximum(srrerr) - minimum(srrerr))/2
# pdscale = 0.03
# plot_potential_components(
#     ycoords,
#     srrerr,
#     initialposition,
#     frequency,
#     amplitude,
#     interfacescale,
#     pdscale
# )

# pdscale = 8e-6
# plot_potential_components(
#     ycoords,
#     jse,
#     initialposition,
#     frequency,
#     amplitude,
#     interfacescale,
#     pdscale,
#     ylabel = "Jump in strain energy",
#     # filepath = foldername * "\\strain-energy-jump.png",
# )

# pdscale = 8e-6
# plot_potential_components(
#     ycoords,
#     jcw,
#     initialposition,
#     frequency,
#     amplitude,
#     interfacescale,
#     pdscale,
#     ylabel = "Jump in compression work",
#     # filepath = foldername * "\\compression-work-jump.png",
# )
#
# plot_dilatation_and_normal_stress(
#     ycoords,
#     parentdilatation,
#     productdilatation,
#     parentsrr,
#     productsrr,
#     # filepath = foldername * "\\dilatation-and-normal-stress.png",
# )

# parentsrrdiff = (maximum(parentsrr) - minimum(parentsrr)) / 2
# yscale = 0.3
# plot_potential_components(
#     ycoords,
#     parentsrr,
#     initialposition,
#     frequency,
#     amplitude,
#     interfacescale,
#     yscale,
#     ylabel = L"\mathrm{Parent} \ \sigma_{nn}",
# )

# productsrrdiff = (maximum(productsrr) - minimum(productsrr)) / 2
# yscale = 0.3
# plot_potential_components(
#     ycoords,
#     productsrr,
#     initialposition,
#     frequency,
#     amplitude,
#     interfacescale,
#     yscale,
#     ylabel = L"\mathrm{Product} \ \sigma_{nn}",
# )


text:
# included in module Wut

struct BitPat
    b::BitArray{1}
    BitPat(b) = new(BitArray(b))
end
export BitPat

function Base.show(io::IO, m::BitPat)
    print(io,"[")
    for (i,v) in enumerate(m.b)
        v ? print(io," 1") : print(io," 0")
    end
    print(io," ]")
end

#end # module


text:
#containing functions trying to applicate the "Soft Actor Critic" from https://arxiv.org/abs/1801.01290
#First we need a policy with a variable gaussian distribution:

@everywhere function grad_all_networks(traj::Trajectory, Vf::Any, Vf_tar::Any, Qf1::Any, Qf2::Any, Pl::Any, params::Params)
    """
    This function determines the gradients for all networks over a Trajectory
    """
    #Take from struct
    len = traj.len
    a_dim = length(params.u_range)
    s_dim = length(params.s_init)
    γ = params.γ_disc

    #Preallocate target data:
    Vf_st_arr      = zeros(len)
    Vf_tar_st2_arr = zeros(len)
    pre_log_pi_arr = zeros(len)
    log_pi_arr     = zeros(len)
    Qf_st_new_arr  = zeros(len)
    Qf1_st_new_arr = zeros(len)
    rand_arr       = zeros(len,a_dim)
    in_q_arr       = zeros(len,s_dim+a_dim)

    """
    Precalculation to generate targets for the loss functions:
    This has to be done with the latest Versions of networks
    """
    for i_t = 1: traj.len
        #Take from struct
        s_t   = traj.s_t[i_t,:] #state
        s_t2  = traj.s_t2[i_t,:] #new state
        a_t   = traj.a_t[i_t,:] #old actions
        p_a_t = traj.p_a_t[i_t,:] #old probabilities
        in_q  = cat(s_t,a_t; dims=1)

        # Pre-evaluation
        Vf_st           = mlp_nlayer_value_norm(Vf,     s_t, params)[1]
        Vf_tar_st2      = mlp_nlayer_value_norm(Vf_tar, s_t2,params)[1]

        #-----------------------------------------------------------------------
        # Generate new acions: MAXIMUM ENTROPY
        #=
        Pl_st, std1     = mlp_nlayer_policy_SAC_norm(Pl,s_t)
        rand_vec        = randn(a_dim)
        a_t_new, sample = sample_action_SAC(Pl_st, std1, rand_vec, params)
        log_pi          = sum(logpdf_SAC(sample,Pl_st,std1))
        =#

        #Generate new actions: CLASSIC
        Pl_st           = mlp_nlayer_policy_norm(Pl,s_t,params)
        rand_vec        = randn(a_dim)
        a_t_new         = sample_action(Pl_st, params.σ, rand_vec)
        p_a_t_new       = get_π_prob(Pl_st,params.σ,a_t_new)
        log_pi          = sum(log.(p_a_t_new))
        #-----------------------------------------------------------------------

        in_q_new        = cat(s_t,a_t_new; dims = 1)
        Qf1_st_new      = mlp_nlayer_value_norm(Qf1,in_q_new,params)[1]
        Qf2_st_new      = mlp_nlayer_value_norm(Qf2,in_q_new,params)[1]

        #Save target data in Arrays:
        Vf_st_arr[i_t]      = Vf_st
        Vf_tar_st2_arr[i_t] = Vf_tar_st2
        pre_log_pi_arr[i_t] = sum(log.(p_a_t))
        log_pi_arr[i_t]     = log_pi
        Qf_st_new_arr[i_t]  = max(Qf1_st_new,Qf2_st_new) #using best of q-learning
        Qf1_st_new_arr[i_t] = Qf1_st_new #using best of q-learning
        rand_arr[i_t,:]     = rand_vec
        in_q_arr[i_t,:]     = in_q
    end

    """
    Gradient of Value function:
    """
    ĝ_vf = ∇_Vf(Vf, traj, Qf_st_new_arr, log_pi_arr, pre_log_pi_arr, params)

    """
    Gradient of Q functions:
    """
    ĝ_qf1 = ∇_Qf(Qf1, traj, Vf_tar_st2_arr, in_q_arr, params)
    ĝ_qf2 = ∇_Qf(Qf2, traj, Vf_tar_st2_arr, in_q_arr, params)

    """
    Gradient of Policy function: Select classic of maximum entropy
    """
    #ĝ_pl = ∇_Pl_SAC(Pl, traj, Qf1_st_new_arr, Vf_st_arr, rand_arr, log_pi_arr, pre_log_pi_arr, params)
    ĝ_pl = ∇_Pl(Pl, traj, Qf1_st_new_arr, Vf_st_arr, rand_arr, log_pi_arr, pre_log_pi_arr, params)

    return ĝ_vf, ĝ_qf1, ĝ_qf2, ĝ_pl
end

#Functions to Calculate the Value function loss + gradient
@everywhere function V_loss(Vf::Any, traj::Trajectory, Qf_st_new_arr::AbstractVector,
    log_pi_arr::AbstractVector, pre_log_pi_arr::AbstractVector, params::Params)

    loss = 0
    for i_t = 1: traj.len
        s_t = traj.s_t[i_t,:]
        V_t = mlp_nlayer_value_norm(Vf, s_t,params)[1]

        #Gather the V-function loss (zeros for classic approach, ones for maximum entropy)
        #loss += (V_t - Qf_st_new_arr[i_t] + log_pi_arr[i_t] * 1 - pre_log_pi_arr[i_t] * 1)^2
        loss += (V_t - Qf_st_new_arr[i_t])^2
    end

    return loss
end
@everywhere ∇_Vf = grad(V_loss,1)

#Functions to Calculate the Q function loss + gradient
@everywhere function Q_loss(Qf::Any, traj::Trajectory, Vf_tar_st2_arr::AbstractVector,
    in_q_arr::Array{Float64,2}, params::Params)

    loss = 0
    γ = params.γ_disc
    for i_t = 1: traj.len
        in_q = in_q_arr[i_t,:]
        r_t = traj.r_t[i_t]
        Q_t = mlp_nlayer_value_norm(Qf,in_q,params)[1]

        #Gather the Q-function loss:
        loss += (Q_t - r_t - γ*Vf_tar_st2_arr[i_t])^2
    end
    return loss
end
@everywhere ∇_Qf = grad(Q_loss,1)

#Functions to Calculate the Policy function loss + gradient
@everywhere function P_loss_SAC(Pl::Any, traj::Trajectory, Qf1_st_new_arr::AbstractVector,
    Vf_st_arr::AbstractVector, rand_arr::Array{Float64,2}, log_pi_arr::Array{Float64,1},
    pre_log_pi_arr::AbstractVector, params::Params)

    fac_mu  = 0
    fac_std = 1e-3
    fac_tan = 0
    u_range = params.u_range

    loss = 0

    for i_t = 1: traj.len
        #Take from struct
        s_t      = traj.s_t[i_t,:]
        rand_vec = rand_arr[i_t,:]

        #Regenerate to get gradient
        out_t, std      = mlp_nlayer_policy_SAC(Pl,s_t)
        sample          = out_t .+ rand_vec .* std
        log_pi          = sum(logpdf_SAC(sample,out_t,std))

        #The main Policy loss:
        policy_loss =  log_pi * (log_pi_arr[i_t]*1 - Qf1_st_new_arr[i_t] +  Vf_st_arr[i_t] - pre_log_pi_arr[i_t]*1)

        #Regression losses:
        mean_reg_loss = mean(fac_mu  * (out_t  .* out_t ))
        std_reg_loss  = mean(fac_std * (std    .* std   ))
        tan_reg_loss  = mean(fac_tan * (sample .* sample))

        #Policy Regression loss:
        policy_reg_loss = mean_reg_loss + std_reg_loss + tan_reg_loss

        #The final policy loss:
        loss += policy_loss + policy_reg_loss
    end

    return loss
end
@everywhere ∇_Pl_SAC = grad(P_loss_SAC,1)

#Functions to Calculate the Policy function loss + gradient
@everywhere function P_loss(Pl::Any, traj::Trajectory, Qf1_st_new_arr::AbstractVector,
    Vf_st_arr::AbstractVector, rand_arr::Array{Float64,2}, log_pi_arr::Array{Float64,1},
    pre_log_pi_arr::AbstractVector, params::Params)

    σ = params.σ

    #Initialize loss
    loss = 0

    for i_t = 1: traj.len
        #Take from struct
        s_t      = traj.s_t[i_t,:]
        rand_vec = rand_arr[i_t,:]

        #Regenerate to get gradient
        Pl_st           = mlp_nlayer_policy_norm(Pl,s_t,params)
        a_t_new         = Pl_st .+ (rand_vec .* σ)
        #p_a_t           = get_π_prob(Pl_st, σ, a_t_new)
        p_a_t           = (1 ./ (σ.*sqrt(2*pi))) .* exp.(((Pl_st .- a_t_new) .* (Pl_st .- a_t_new)) ./ (-2 .* σ .* σ))
        log_pi          = sum(log.(p_a_t))

        #The main Policy loss:
        policy_loss =  log_pi * (Qf1_st_new_arr[i_t] -  Vf_st_arr[i_t])

        #The final policy loss:
        loss += policy_loss
    end

    return loss
end
@everywhere ∇_Pl = grad(P_loss,1)

max_stars_repo_name: NikonPic/REINFORCEMENT-JULIA

max_issues_repo_name: NikonPic/REINFORCEMENT-JULIA

max_forks_repo_name: NikonPic/REINFORCEMENT-JULIA


text:
abstract type FormatSize end

struct SizeExact <: FormatSize end
struct SizeUpperbound <: FormatSize end
struct SizeUnknown <: FormatSize end

FormatSize(::Type{<:Integer}) = SizeExact()


text:
module MFCC

using DSP, FFTW
include("Mel.jl")
include("Filterbanks.jl")
using MFCC.Mel
using MFCC.Filterbanks

"""
    mfcc(signal::AbstractVector, window::Function = ones, windowsize = 400, overlap = 160, sr = 44100, numcep = 26)

Compute the Mel-Frequency Cepstral Coefficients of the given signal and return
an Array where each column contains the coefficients of a given frame obtained
from the stft.

**Parameters:**

`signal`: The signal on which the coefficients will be computed.

`window`: The window function to be applied on the signal by the STFT.

`windowsize`: The size of the STFT window.

`overlap`: The overlap of two adjacent windows of the STFT in samples.

`sr`: Sample rate of the signal.

`numcep`: The amount of cepstral coefficients to be extracted from the signal.

See also: [`Filterbank`](@ref)
"""
function mfcc(signal::AbstractVector,
	window::Function = ones,
	windowsize::Int = 400,
	overlap::Int = 160,
	sr::Int = 44100,
	numcep::Int = 26
)

	#Each column is a frame of stft
	sigstft = stft(signal, windowsize, overlap, window=window)

	#If the signal is real then the absolute value of it's spectrum is symmetric
	eltype(signal)<:Real && (windowsize = windowsize ÷ 2 + 1; sr ÷ 2)

	#Evaluate the periodogram of the signal
	periodgrm = map(frame->(abs2.(frame)/windowsize), sigstft)

	#Create the filterbank
	windows, windowstarts = mffilters(numcep, windowsize, sr)
	melfilterbank = Filterbank(windowsize, windows[1:end], windowstarts[1:end])

	#Apply the filterbank and sum the result of each filter
	filteredsignal = mapslices(frame->melfilterbank*frame, periodgrm, dims = 1)

	#Take the Type-III DCT of the log for each coeficient
	return mapslices(frame->idct(log.(frame)), filteredsignal, dims=1)

end

export mfcc

end # module


text:
# this file is the Julia code for the system of two approaches in Christensenet al.(2021) when all individuals have all omics data
# the R code can be found in Christensenet al.(2021)
using JWAS,DataFrames,CSV,Statistics,JWAS.Datasets,DelimitedFiles,Random, LinearAlgebra, StatsBase
Random.seed!(1)

res_all=zeros(20)

for data in 1:20  #data1,...,data20
    n      = 1055
    nTrain = 955
    nTest  = n-nTrain
    nfeat  = 1200  #nOmics features

    σ2_α  = 0.01   #variance for omics effect α
    σ2_mi = 2      #variance for omics each feature mi
    h2    = 0.47   #total heritability in equation (1)
    h2_m  = 0.61   #omics heritability in equation (2)
    h2_ar = 0.2    #residual bv heritability, part of h2.

    nOmics = 1200

    # calculate variance component for data simulation
    σ2_omics_gi  = h2_m * σ2_mi
    σ2_omics_ei  = (1-h2_m)*σ2_mi
    σ2_am        = nOmics * σ2_omics_gi * σ2_α #total omics genetic variance (=var(Gα))
    σ2_y         = σ2_am/(h2-h2_ar)
    σ2_ar        = h2_ar * σ2_y
    σ2_ϵ         = σ2_y - nOmics*σ2_mi*σ2_α - σ2_ar

    zeta   = σ2_omics_ei/σ2_omics_gi
    eta1   = σ2_ϵ/σ2_α

    #read GRM
    GRM=readdlm("/Users/tianjing/Box/Omicsdata/andrew1055/data$data/GRM_geno.txt",'\t')
    size(GRM)
    GRM=GRM+I*0.00001
    Ginv=inv(GRM)

    #read omics data
    M = readdlm("/Users/tianjing/Box/Omicsdata/andrew1055/data$data/omics.txt",'\t')

    #read phenotype
    yy = vec(readdlm("/Users/tianjing/Box/Omicsdata/andrew1055/data$data/y_no_ar.txt",','))

    ### computing G.hat:  (5) in paper Christensenet al.(2021)
    LHS11 = ones(n)'ones(n)
    LHS12 = ones(n)'I
    LHS22 = I'I + zeta*Ginv

    LHS = [LHS11  LHS12
           LHS12' LHS22]
    RHS=[ones(n)'M
         I'M]
    sol_all=inv(LHS)*RHS
    G_hat = sol_all[2:end,:]
    size(G_hat)

    ### (4) in paper Christensenet al.(2021)
    LHS11 = ones(nTrain)'ones(nTrain)
    LHS12 = ones(nTrain)'M[1:nTrain,:]
    LHS22 = M[1:nTrain,:]' * M[1:nTrain,:] + eta1*I

    LHS = [LHS11  LHS12
           LHS12' LHS22]
    RHS=[ones(nTrain)'yy[1:nTrain]
         M[1:nTrain,:]'yy[1:nTrain]]

    SOL = inv(LHS)*RHS
    alpha=SOL[2:end]

    hat_a_m = G_hat*alpha

    ### prediction accuracy
    tbv=vec(readdlm("/Users/tianjing/Box/Omicsdata/andrew1055/data$data/tbv_m.txt"))

    @show cor(hat_a_m,tbv)
    @show cor(hat_a_m[956:1055],tbv[956:1055])

    println("----data$data")
    @show round(cor(hat_a_m,tbv),digits=3)
    @show round(cor(hat_a_m[956:1055],tbv[956:1055]),digits=3)
    res_all[data]=cor(hat_a_m[956:1055],tbv[956:1055])
end

res_all

writedlm("/Users/tianjing/Box/omics/legarra_fullomics_res_20data.txt",res_all)

max_stars_repo_name: zhaotianjing/NN-LMM

max_issues_repo_name: zhaotianjing/NN-LMM

max_forks_repo_name: zhaotianjing/NN-LMM


text:
# TODO: The scaling of pseudo-matrices with an element scales the ideals and not the matrix ...

function _translate_pseudo_hnf(H::PMat, f)
  OL = maximal_order(codomain(f))
  return _translate_pseudo_hnf(H, f, OL)
end

function _translate_pseudo_hnf(H::PMat, f, target_order)
  OE = target_order
  E = nf(OE)
  coeff_ideals = fractional_ideal_type(OE)[]
  for i in 1:length(H.coeffs)
    push!(coeff_ideals, fractional_ideal(OE, elem_type(E)[f(x) for x in gens(H.coeffs[i])]))
  end
  m = zero_matrix(E, nrows(H), ncols(H))
  MH = matrix(H)
  for i in 1:nrows(H)
    for j in 1:ncols(H)
      m[i, j] = f(MH[i, j])
    end
  end
  pm = pseudo_matrix(m, coeff_ideals)
  return pm
end

function pseudo_hnf_via_absolute(H::PMat; full_rank::Bool = true, shape::Symbol = :lowerleft, nonzero::Bool = false)
  O = base_ring(H)
  E = nf(O)
  Eabs, EabsToE, KtoE = absolute_field(E)
  return pseudo_hnf_via_absolute(H, EabsToE, full_rank = full_rank,
                                             shape = shape,
                                             nonzero = nonzero)
end

function pseudo_hnf_via_absolute(H::PMat, EabsToE; full_rank::Bool = true, shape::Symbol = :lowerleft, nonzero::Bool = false)
  O = base_ring(H)
  E = nf(O)
  pm = _translate_pseudo_hnf(H, pseudo_inv(EabsToE))
  if full_rank
    HH = pseudo_hnf_full_rank(pm, shape)
  else
    HH = pseudo_hnf(pm, shape)
  end

  if nonzero
    r = 0
    if shape === :lowerleft
      for i in 1:nrows(HH)
        if !iszero_row(matrix(HH), i)
          r = i
          break
        end
      end
      HH = sub(HH, r:nrows(pm), 1:ncols(pm))
    elseif shape === :upperright
      for i in nrows(HH):-1:1
        if !iszero_row(matrix(HH), i)
          r = i
          break
        end
      end
      H = sub(HH, 1:r, 1:ncols(HH))
    end
  end

  return _translate_pseudo_hnf(HH, EabsToE, O)
end

function _sum_modules_with_map(a::PMat{<: NfRelElem, <: NfRelOrdFracIdl}, b::PMat, f, full_rank = true)
  H = vcat(a, b)
  return pseudo_hnf_via_absolute(H, f, shape = :lowerleft, nonzero = true)
end

_sum_modules(L::QuadLat, a::PMat, b::PMat, full_rank = true) = _sum_modules(a, b, full_rank)

function _sum_modules(L::HermLat, a::PMat, b::PMat, full_rank = true)
  _,f,_ = absolute_field(ambient_space(L))
  return _sum_modules_with_map(a, b, f, full_rank)
end

function _sum_modules(a::PMat{<: NfRelElem, <: NfRelOrdFracIdl}, b::PMat, full_rank = true)
  H = vcat(a, b)
  return pseudo_hnf_via_absolute(H, shape = :lowerleft, nonzero = true)
end

function _sum_modules(a::PMat, b::PMat, full_rank = true)
  H = vcat(a, b)
  H = pseudo_hnf(H, :lowerleft)
  r = 0
  for i in 1:nrows(H)
    if !iszero_row(matrix(H), i)
      r = i
      break
    end
  end
  @assert r != 0
  return sub(H, r:nrows(H), 1:ncols(H))
end

function _intersect_modules(a::PMat{<: NfRelElem, <: NfRelOrdFracIdl}, b::PMat, full_rank = true)
  _, f, _ = absolute_field(nf(base_ring(a)))
  return _intersect_modules_with_map(a, b, f, full_rank)
end

function _intersect_modules_with_map(a::PMat{<: NfRelElem, <: NfRelOrdFracIdl}, b::PMat, f, full_rank = true)
  OE = maximal_order(domain(f))
  aE = _translate_pseudo_hnf(a, pseudo_inv(f), OE)
  bE = _translate_pseudo_hnf(b, pseudo_inv(f), OE)
  c = _intersect_modules(aE, bE, full_rank)
  return _translate_pseudo_hnf(c, f, base_ring(a))
end

function _intersect_modules(a::PMat, b::PMat, full_rank = true)
  M1 = hcat(a, deepcopy(a))
  d = ncols(b)
  z = zero_matrix(base_ring(matrix(a)), d, d)
  M2 = hcat(pseudo_matrix(z, b.coeffs), b)
  M = vcat(M1, M2)
  if full_rank
    H = sub(pseudo_hnf(M, :lowerleft), 1:d, 1:d)
    return H
  else
    H = pseudo_hnf_kb(M, :lowerleft)
    i = 1
    while iszero_row(H.matrix, i)
      i += 1
    end
    return sub(H, i:d, 1:d)
  end
end

function _modules_equality(a::PMat, b::PMat)
  _spans_subset_of_pseudohnf(a, b, :lowerleft) && _spans_subset_of_pseudohnf(b, a, :lowerleft)
end

function _module_scale_ideal(a::NfAbsOrdIdl, b::PMat)
  return pseudo_matrix(matrix(b), [ a * c for c in coefficient_ideals(b)])
end

_module_scale_ideal(a::PMat, b::NfAbsOrdIdl) = _module_scale_ideal(b, a)

function _module_scale_ideal(a::NfOrdFracIdl, b::PMat)
  return pseudo_matrix(matrix(b), Ref(a) .* coefficient_ideals(b))
end

_module_scale_ideal(a::PMat, b::NfOrdFracIdl) = _module_scale_ideal(b, a)

function _module_scale_ideal(a::NfRelOrdIdl, b::PMat)
  return pseudo_matrix(matrix(b), Ref(a) .* coefficient_ideals(b))
end

_module_scale_ideal(a::PMat, b::NfRelOrdIdl) = _module_scale_ideal(b, a)

function _module_scale_ideal(a::NfRelOrdFracIdl, b::PMat)
  return pseudo_matrix(matrix(b), Ref(a) .* coefficient_ideals(b))
end

_module_scale_ideal(a::PMat, b::NfRelOrdFracIdl) = _module_scale_ideal(b, a)

*(a::NfRelOrdFracIdl, b::PMat) = _module_scale_ideal(a, b)

################################################################################
#
#  Local basis matrices
#
################################################################################

# Given a pseudo matrix over K with row span M and p a prime ideal of K, find a
# basis matrix of M \otimes OK_p.
function _local_basis_matrix(a::PMat, p::NfOrdIdl)
  @assert base_ring(a) == order(p)
  uni = uniformizer(p)
  z = zero_matrix(base_ring(matrix(a)), nrows(a), ncols(a))
  for i in 1:nrows(a)
    c = coefficient_ideals(a)[i]
    x = uni^valuation(c, p)
    for j in 1:ncols(a)
      z[i, j] = x * matrix(a)[i, j]
    end
  end
  return z
end

# Given a pseudo matrix over E/K with row span M and a p a prime ideal of K,
# find a basis matrix of M \otimes OE_p.
function _local_basis_matrix_prime_below(a::PMat, p::T) where T
  @assert base_ring(base_ring(a)) == order(p)
  R = base_ring(a)
  D = prime_decomposition(R, p)
  unis = elem_type(R)[p_uniformizer(q[1]) for q in D]
  @assert all(valuation(unis[i], D[i][1]) == 1 for i in 1:length(D))
  @assert all(sum(valuation(unis[i], D[j][1]) for j in 1:length(D)) == 1 for i in 1:length(D))
  z = zero_matrix(base_ring(matrix(a)), nrows(a), ncols(a))
  _c = coefficient_ideals(a)
  for i in 1:nrows(a)
    c = _c[i]
    x = unis[1]^(valuation(c, D[1][1]))
    for k in 2:length(D)
      x = x * unis[k]^valuation(c, D[k][1])
    end
    for j in 1:ncols(a)
      z[i, j] = x * matrix(a)[i, j]
    end
  end
  return z
end

# Given a pseudo matrix over E/K with row span M and a p a prime ideal of K,
# find a basis matrix N of M \otimes OE_p. The row span of N is a submodule of
# M.
function _local_basis_matrix_prime_below_submodule(a::PMat, p)
  @assert base_ring(base_ring(a)) == order(p)
  R = base_ring(a)
  D = prime_decomposition(R, p)
  z = zero_matrix(base_ring(matrix(a)), nrows(a), ncols(a))
  _c = coefficient_ideals(a)
  for i in 1:nrows(a)
    c = _c[i]
    f = factor(c)
    if length(f) == 0
      x = one(nf(R))
    else
      for Q in D
        if !(haskey(f, Q[1]))
          f[Q[1]] = 0
        end
      end
      x = approximate(Int[e for (_, e) in f], ideal_type(base_ring(a))[p for (p, _) in f])
      #@assert all(valuation(x, p) == e for (p, e) in f)
    end
    #@assert valuation(x, D[1][1]) == valuation(c, D[1][1])
    #x = unis[1]^valuation(c, D[1][1])
    #for k in 2:length(D)
    #  x = x * unis[k]^valuation(c, D[k][1])
    #end
    for j in 1:ncols(a)
      z[i, j] = x * matrix(a)[i, j]
    end
  end
  if false
    _z = pseudo_matrix(z, [one(nf(R)) * R for i in 1:nrows(z)])
    @assert _spans_subset_of_pseudohnf(_z, a, :lowerleft)
    @assert valuation(det(_z), D[1][1]) == valuation(det(a), D[1][1])
  end
  return z
end

# Given a pseudo matrix over K with row span M and a p a prime ideal of K,
# find a basis matrix N of M \otimes OK_p. The row span of N is a submodule of
# M.
function _local_basis_submodule_matrix(a::PMat, p)
  #@assert base_ring(base_ring(a)) == order(p)
  z = zero_matrix(base_ring(matrix(a)), nrows(a), ncols(a))
  for i in 1:nrows(a)
    c = coefficient_ideals(a)[i]
    vpc = valuation(c, p)
    found = false
    local x
    for b in absolute_basis(c) # could use generators here
      if valuation(b, p) == vpc
        found = true
        x = b
        break
      end
    end
    @assert found
    for j in 1:ncols(a)
      z[i, j] = x * matrix(a)[i, j]
    end
  end
  return z
end

function _local_basis_supermodule_matrix(a::PMat, p::NfOrdIdl)
  throw(NotImplemented())
end


text:
@testset "test_simplification" begin
## SMINREAL ##
G = ss([-5 0 0 0; 0 -1 -2.5 0; 0 4 0 0; 0 0 0 -6], [2 0; 0 1; 0 0; 0 2],
       [0 3 0 0; -2 0 0 1], [0 0; 1 0])
@test sminreal(G) == G
@test sminreal(G[1, 1]) == ss(0)
@test sminreal(G[1, 2]) == ss([-1 -2.5; 4 0], [1; 0], [3 0], [0])
@test sminreal(G[2, 1]) == ss([-5], [2], [-2], [1])
@test sminreal(G[2, 2]) == ss([-6], [2], [1], [0])

## MINREAL ##

s = tf("s")
P = [1/(s+1) 2/(s+3); 1/(s+1) 1/(s+1)]
sys = ss(P)
sysmin = minreal(sys)

@test size(sysmin.A,1) == 3 # Test that the reduction of sys worked

@test norminf(sys - sysmin)[1] < 1e-15 # And that the answer is correct

@test_broken balreal(sys-sysmin)

@test all(sigma(sys-sysmin, [0.0, 1.0, 2.0])[1] .< 1e-15)  # Previously crashed because of zero dimensions in tzero

t = 0:0.1:10
y1,x1 = step(sys,t)[[1,3]]
y2,x2 = step(sysmin,t)[[1,3]]
@test sum(abs2,y1.-y2) < 1e-6 # Test that the output from the two systems are the same

end


text:
@testset "Vertex" begin
    v = Vertex(1, -1)
    vr = repr(v)
    @test occursin("ID: $(v.id)\n", vr)
    @test occursin("Coords: ($(v.x), $(v.y))\n", vr)

    @test occursin("Edge: $(v.edge === nothing ? "nothing" : v.edge.id)\n", vr)
end

@testset "Edge" begin
    e = Edge()
    er = repr(e)
    @test occursin("ID: $(e.id)\n", er)
    @test occursin("Next: $(e.next === nothing ? "nothing" : e.next.id)\n", er)
    @test occursin("Prev: $(e.prev === nothing ? "nothing" : e.prev.id)\n", er)
    @test occursin("Twin: $(e.twin === nothing ? "nothing" : e.twin.id)\n", er)
    @test occursin("Dest: $(e.dest === nothing ? "nothing" : e.dest.id)\n", er)
    @test occursin("Face: $(e.face === nothing ? "nothing" : e.face.id)\n", er)
end

@testset "Face" begin
    f = Face()
    fr = repr(f)
    @test occursin("ID: $(f.id)\n", fr)
    @test occursin("Edge: $(f.edge === nothing ? "nothing" : f.edge.id)\n", fr)
end

@testset "Tessellation" begin
    tess = Tessellation()
    tr = repr(tess)
    @test occursin("$(length(tess.vertices)) vertices\n", tr)
    @test occursin("$(length(tess.edges)) edges\n", tr)
    @test occursin("$(length(tess.faces)) faces\n", tr)
end


text:
module Convex1d
export minimize

let
    s = read(joinpath(dirname(@__DIR__), "README.md"), String)
    s = replace(s, "```julia" => "```jldoctest")
    @doc s Convex1d
end

function middle(x,y)
    # lets not depend on Statistics for this
    mid = (1//2 * x) + (1//2*y)
    @assert x < mid <y
    mid
end

function find_initial_domain(f, (x_left, x_right)=(-1, 1))

    function grow_once(;x_from,y_from,x_to,y_to)
        if y_to > y_from
            return (x_from=x_from, y_from=y_from, x_to=x_to, y_to=y_to, converged=true)
        end
        step = x_to - x_from
        if step == 0
            msg = """
            Zero step in domain growing:
            x_from = $x_from
            y_from = $y_from
            x_to   = $x_to
            y_to   = $y_to
            step   = $step
            """
            error(msg)
        end
        x_to_new = step + x_to
        y_to_new = f(x_to_new)
        return (x_from=x_from, y_from=y_from, x_to=x_to_new, y_to=y_to_new, converged=false)
    end

    function grow(;x_from, y_from, x_to, y_to)
        item = grow_once(x_from=x_from, y_from=y_from, x_to=x_to, y_to=y_to)
        while !(item.converged)
            item = grow_once(x_from=item.x_from, y_from=item.y_from, x_to=item.x_to, y_to=item.y_to)
        end
        return item
    end

    x_left = float(x_left)
    x_right = float(x_right)
    x_center = middle(x_left, x_right)
    @assert x_left < x_center < x_right
    y_left = f(x_left)
    y_center = f(x_center)
    y_right = f(x_right)
    left = grow(x_from=x_center, y_from=y_center, x_to=x_left, y_to=y_left)
    right = grow(x_from=x_center, y_from=y_center, x_to=x_right, y_to=y_right)
    @assert left.x_to  < x_center < right.x_to
    @assert left.y_to  >= y_center
    @assert right.y_to >= y_center
    (left.x_to, right.x_to)
end

"""
    minimize(f, (x_left, x_right); atol)

Minimize a convex function `f` on the interval `[x_left, x_right]`.

# Arguments
* f: Function that should be minimized, does not have to be differentiable.
* atol: Search tolereance, such that `isapprox(true_minimizer, sol.minimizer, atol=atol)`
will hold.
"""
function minimize(f, (x_left, x_right)=find_initial_domain(f);
        atol=nothing
    )
    if !(x_left < x_right)
        msg = """
        x_left < x_right must hold. Got
        x_left  = $x_left
        x_right = $x_right
        """
        throw(ArgumentError(msg))
    end
    if atol === nothing
        atol = max(sqrt(eps(x_left)), sqrt(eps(x_right)))
    end
    # idea:
    # subdivide interval at five points
    # evaluate f at each point
    # recurse using x_left and x_right the neighbors of the minimizing point
    #
    # also we carefully reuse evals of f

    x1 = x_left
    x5 = x_right
    x3 = middle(x1, x5)
    x2 = middle(x1, x3)
    x4 = middle(x3, x5)
    xs = promote(x1,x2,x3,x4,x5)
    fs = map(f, xs)
    minimize_interval_subdivision_kernel(f, xs, fs, atol)
end

@noinline function minimize_interval_subdivision_kernel(f::F, xs, fs, atol) where {F}
    nevals = 5
    while true
        i = argmin(fs)
        at_left_bdry  = i == 1
        at_right_bdry = i == 5
        at_bdry = at_left_bdry | at_right_bdry
        if at_left_bdry
            ilo = 1; ihi = 2
        elseif at_right_bdry
            ilo = 4; ihi = 5
        else # interiour
            ilo = i-1
            ihi = i+1
        end
        @assert ilo < ihi
        x1 = xs[ilo]
        x5 = xs[ihi]
        f1 = fs[ilo]
        f5 = fs[ihi]
        if at_bdry
            x3 = middle(x1, x5)
            f3 = f(x3) ; nevals += 1
        else
            x3 = xs[i]
            f3 = fs[i]
        end

        converged = (x5 - x1) < atol
        if converged
            return minimize_interval_subdivision_finish((x1,x3,x5), (f1,f3,f5), nevals)
        end

        x2 = middle(x1, x3)
        f2 = f(x2); nevals += 1
        x4 = middle(x3, x5)
        f4 = f(x4); nevals += 1
        xs = (x1,x2,x3,x4,x5)
        fs = (f1,f2,f3,f4,f5)
    end
end

function minimize_interval_subdivision_finish(xs::NTuple{3}, fs::NTuple{3}, nevals)
    x1,x2,x3 = xs
    f1,f2,f3 = fs
    i_opt    = argmin(fs)
    # calculations for ybounds are too unstable
    # m_left   = (f2 - f1) / (x2 - x1)
    # m_right  = (f3 - f2) / (x3 - x2)
    # y_left   = f2 - m_right * (x2 - x1)
    # y_right  = f2 + m_left  * (x3 - x2)
    # @assert y_left <= f1
    # @assert y_right <= f3
    if i_opt == 1
        xbounds = (x1, x2)
        # ybounds = (y_left, f1)
    elseif i_opt == 2
        xbounds = (x1, x3)
        # ybounds = minmax(y_left, y_right)
    elseif i_opt == 3
        xbounds = (x2, x3)
        # ybounds = (y_right, f3)
    else
        error("Unreachable")
    end
    @assert issorted(xbounds)
    # @assert issorted(ybounds)
    return (
        minimizer = xs[i_opt],
        minimum   = fs[i_opt],
        nevals    = nevals,
        xbounds   = xbounds,
        # ybounds   = ybounds,
    )

end

end # module


text:
module Ex5_4c
using FaultDetectionTools
using DescriptorSystems
using LinearAlgebra
using Polynomials
using Test

# Example 5.4c - Solution of an EFDP using EFDSYN
println("Example 5.4c")

# define s as an improper transfer function
s = rtf('s');
# define Gu(s) and Gd(s)
Gu = [(s+1)/(s-2); (s+2)/(s-3)];     # enter Gu(s)
Gd = [(s-1)/(s+2); 0];               # enter Gd(s)
Gf = [(s+1)/(s-2) 0; (s+2)/(s-3) 1]; # enter Gf(s)
p = 2; mu = 1; md = 1; mf = 2;       # set dimensions

# setup the synthesis model with faults with Gf(s) = [ Gu(s) [0;1]]
sysf = fdimodset(dss([Gu Gd]),c =1,d = 2,f = 1,fs = 2);

# call of EFDSYN with the options for stability degree -3 and the synthesis
# of a scalar output filter
Qt, Rft, info = efdsyn(sysf, sdeg = -3, rdim = 1);

# normalize Q and Rf to match example
scale = evalfr(Rft.sys[1,1],Inf)[1,1];
Q = FDFilter(Qt.sys/scale,Qt.outputs,Qt.controls);
Rf = FDFilterIF(Rft.sys/scale,faults = Rft.faults);
println("Q = $(dss2rm(Q.sys, atol = 1.e-7))")
println("Rf = $(dss2rm(Rf.sys, atol = 1.e-7))")

# check synthesis results
R = fdIFeval(Q, sysf; minimal = true, atol = 1.e-7)
@test iszero(R.sys[:,[R.controls; R.disturbances]], atol = 1.e-7) &&
      iszero(R.sys[:,R.faults]-Rf.sys) &&
      fdisspec(Rf) == Bool[1 1] &&
      gpole(Q.sys) ≈ [-3] && gpole(Rf.sys) ≈ [-3]

end # module

max_stars_repo_name: andreasvarga/FaultDetectionTools

max_issues_repo_name: andreasvarga/FaultDetectionTools

max_forks_repo_name: andreasvarga/FaultDetectionTools


text:
using DataEnvelopmentAnalysis
using LinearAlgebra
using SparseArrays
using Test

@testset "DataEnvelopmentAnalysis" begin

    include("dea.jl")
    include("deaadd.jl")
    include("deaddf.jl")
    include("deagdf.jl")
    include("dearussell.jl")
    include("deaerg.jl")
    include("deamddf.jl")
    include("deaholder.jl")
    include("dearddf.jl")

    include("econoptim.jl")
    include("deacost.jl")
    include("dearevenue.jl")
    include("deaprofit.jl")
    include("deaprofitability.jl")

    include("malmquist.jl")

    include("peers.jl")

end


text:
using ADCME

db = Database("simulation.db")

execute(db, """
CREATE TABLE IF NOT EXISTS acoustic
(
    level text,
    size integer,
    forward real,
    backward real,
    primary key (level, size)
)
""")

function insert_record(level, size, t1, t2)
    execute(db, """
    INSERT OR REPLACE into acoustic VALUES (\"$level\", $size, $t1, $t2)
    """)
    close(db)
end


text:
using PackageCompilerConvenienceTools
using Test

import Pkg

@testset "PackageCompilerConvenienceTools.jl" begin
    @testset "utils.jl" begin
        PackageCompilerConvenienceTools._with_different_project(PackageCompilerConvenienceTools.default_project()) do activated_project
        end
    end
    @testset "public.jl" begin
        @testset "fast_and_simple" begin
            PackageCompilerConvenienceTools._with_temp_project() do my_temp_project
                PackageCompilerConvenienceTools._with_temp_dir() do my_temp_dir
                    Pkg.add("Crayons")
                    sysimage_path = joinpath(my_temp_dir, "ExampleSysimage.so")
                    @test !isfile(sysimage_path)
                    @test !ispath(sysimage_path)
                    @test_throws ArgumentError PackageCompilerConvenienceTools.fast_and_simple(my_temp_project;
                                                                                shared_project = false,
                                                                                replace_default = true)
                    @test !isfile(sysimage_path)
                    @test !ispath(sysimage_path)
                    PackageCompilerConvenienceTools.fast_and_simple(my_temp_project;
                                                          shared_project = false,
                                                          sysimage_path = sysimage_path)
                    @test isfile(sysimage_path)
                    @test ispath(sysimage_path)
                    p_1 = run(`$(Base.julia_cmd()) -J$(sysimage_path) --project=$(my_temp_project) -e 'import Crayons'`)
                    wait(p_1)
                    @test success(p_1)
                    p_2 = run(`$(Base.julia_cmd()) -J$(sysimage_path) --project=$(my_temp_project) -e 'import Pkg; Pkg.test("Crayons")'`)
                    wait(p_2)
                    @test success(p_2)
                end
            end
        end
    end

    if lowercase(strip(get(ENV, "CI", ""))) == "true"
        @testset "CI-only tests" begin
            @testset "CI public.jl" begin
                @testset "CI fast_and_simple" begin
                    PackageCompilerConvenienceTools._with_temp_project() do my_temp_project
                        PackageCompilerConvenienceTools._with_temp_dir() do my_temp__dir
                            Pkg.add("Crayons")
                            @test_throws ArgumentError PackageCompilerConvenienceTools.fast_and_simple(my_temp_project;
                                                                                shared_project = false,
                                                                                replace_default = true)
                            if !(Sys.iswindows())
                                PackageCompilerConvenienceTools.fast_and_simple(my_temp_project;
                                                                      shared_project = false,
                                                                      replace_default = true,
                                                                      force_replace_default = true)
                                p_1 = run(`$(Base.julia_cmd()) --project=$(my_temp_project) -e 'import Crayons'`)
                                wait(p_1)
                                @test success(p_1)
                                p_2 = run(`$(Base.julia_cmd()) --project=$(my_temp_project) -e 'import Pkg; Pkg.test("Crayons")'`)
                                wait(p_2)
                                @test success(p_2)
                            end
                        end
                    end
                end
            end
        end
    end
end


text:
function assign_lines_to_orders(line_list::DataFrame, order_info::DataFrame; v_center::Real=0, Δv_to_avoid_tellurics::Real = RvSpectMLBase.max_bc )
  @assert hasproperty(line_list, :lambda)
  @assert hasproperty(line_list, :lambda_lo)
  @assert hasproperty(line_list, :lambda_hi)
  ll1 = copy(line_list)
  ll2 = copy(line_list)
  ll3 = copy(line_list)
  ll1[!,:order] .= 0
  ll2[!,:order] .= 0
  ll3[!,:order] .= 0
  if !(hasproperty(line_list,:lambda_lo) && hasproperty(line_list,:lambda_hi))
    add_line_boundaries_to_line_list(line_list; Δv_to_avoid_tellurics=Δv_to_avoid_tellurics, v_center_to_avoid_tellurics=v_center )
  end
  for (i,line) in enumerate(eachrow(line_list))
    order_indices = findall( map(order-> order.λ_min <= line.lambda_lo <= order.λ_max && order.λ_min <= line.lambda_hi <= order.λ_max, eachrow(order_info) ) )
    if length(order_indices) == 1
      ll1[i,:order] = order_info.order[order_indices[1]]
    elseif length(order_indices) == 2
      ll1[i,:order] = order_info.order[order_indices[1]]
      ll2[i,:order] = order_info.order[order_indices[2]]
    elseif length(order_indices) == 3
      ll1[i,:order] = order_info.order[order_indices[1]]
      ll2[i,:order] = order_info.order[order_indices[2]]
      ll3[i,:order] = order_info.order[order_indices[3]]
    elseif length(order_indices) > 3
      @warn "One line appears in >3 orders.  Not implemented yet."
      println("line = ", line, "  order_indices = ", order_indices, "  line.order = ",order_info.order[order_indices])
      #=  df_tmp = DataFrame()
      for key in keys(eachcol(line_list))
        df_tmp[!,key] = fill(line[key], length(order_indices)-1 )
      end
      df_tmp[:order] .= order_indices[2:end]
      append!(line_list,df_tmp)
      =#
      #break
    end
  end
  append!(ll1,ll2)
  append!(ll1,ll3)
  ll1 |> @filter( _.order >= 1 ) |> @orderby(_.order) |> @thenby(_.lambda_lo) |> DataFrame
  #ll1 |> @filter( _.order >= 1 ) |> @orderby(_.lambda_lo) |> @thenby(_.order) |> DataFrame
end

function add_line_boundaries_to_line_list(line_list::DataFrame; Δv_to_avoid_tellurics::Real = 0.0, v_center_to_avoid_tellurics::Real = 0.0 )
    @assert 0.5*RvSpectMLBase.max_bc_earth_rotation <= Δv_to_avoid_tellurics <= 4*RvSpectMLBase.max_bc
    @assert hasproperty(line_list,:lambda)
    @assert !hasproperty(line_list,:lambda_lo)
    @assert !hasproperty(line_list,:lambda_hi)
    line_list_to_search_for_tellurics = copy(line_list)
    line_list_to_search_for_tellurics.lambda    = line_list_to_search_for_tellurics.lambda.*(calc_doppler_factor(v_center_to_avoid_tellurics))
    line_list_to_search_for_tellurics.lambda_lo = line_list_to_search_for_tellurics.lambda.*(calc_doppler_factor(v_center_to_avoid_tellurics)/calc_doppler_factor(Δv_to_avoid_tellurics))
    line_list_to_search_for_tellurics.lambda_hi = line_list_to_search_for_tellurics.lambda.*(calc_doppler_factor(v_center_to_avoid_tellurics)*calc_doppler_factor(Δv_to_avoid_tellurics))
    return line_list_to_search_for_tellurics
end

function expand_line_boundaries_in_line_list(line_list::DataFrame; Δv_to_avoid_tellurics::Real = 0.0, v_center_to_avoid_tellurics::Real = 0.0 )
    @assert 0.5*RvSpectMLBase.max_bc_earth_rotation <= Δv_to_avoid_tellurics <= 4*RvSpectMLBase.max_bc
    @assert hasproperty(line_list,:lambda)
    @assert hasproperty(line_list,:lambda_lo)
    @assert hasproperty(line_list,:lambda_hi)
    line_list_to_search_for_tellurics = copy(line_list)
    line_list_to_search_for_tellurics.lambda    = line_list_to_search_for_tellurics.lambda.*(calc_doppler_factor(v_center_to_avoid_tellurics))
    line_list_to_search_for_tellurics.lambda_lo = line_list_to_search_for_tellurics.lambda_lo.*(calc_doppler_factor(v_center_to_avoid_tellurics)/calc_doppler_factor(Δv_to_avoid_tellurics))
    line_list_to_search_for_tellurics.lambda_hi = line_list_to_search_for_tellurics.lambda_hi.*(calc_doppler_factor(v_center_to_avoid_tellurics)*calc_doppler_factor(Δv_to_avoid_tellurics))
    return line_list_to_search_for_tellurics
end

max_stars_repo_name: RvSpectML/RvSpectMLTemplate

max_issues_repo_name: RvSpectML/RvSpectMLTemplate

max_forks_repo_name: RvSpectML/RvSpectMLTemplate