text:
function maplight_get(auth::Auth, path, args; query = Dict(), options...)
    authentiate_query(query, auth)
    merge!(query, args)

    r = get(URI(MAPLIGHT_API; path = path); query = query, options...)

    handle_error(r)

    JSON.parse(r.data)
end

# -------

function bill_positions(auth::Auth, jurisdiction, session, prefix, number; options...)
    maplight_get(auth, "/services_open_api/map.bill_positions_v1.json", {
        "jurisdiction" => jurisdiction,
        "session" => session,
        "prefix" => prefix,
        "number" => number
        }; options...)
end

bill_positions(auth::String, jurisdiction, session, prefix, number; options...) = bill_positions(authenticate(auth), jurisdiction, session, prefix, number; options...)

bill_positions(jurisdiction, session, prefix, number; auth = authenticate(), options...) = bill_positions(auth, jurisdiction, session, prefix, number; options...)

# -------

function bill_search(auth::Auth, jurisdiction, search_str; options...)
    maplight_get(auth, "/services_open_api/map.bill_search_v1.json", { "jurisdiction" => jurisdiction, "search" => search_str }; options...)
end

bill_search(auth::String, jurisdiction, search_str; options...) = bill_search(authenticate(auth), jurisdiction, search_str; options...)

bill_search(jurisdiction, search_str; auth = authenticate(), options...) = bill_search(auth, jurisdiction, search_str; options...)

# -------

function bill_list(auth::Auth, jurisdiction, session; include_organizations = true, has_organizations = true, options...)
    maplight_get(auth, "/services_open_api/map.bill_list_v1.json", {
            "jurisdiction" => jurisdiction,
            "session" => session,
            "include_organizations" => include_organizations,
            "has_organizations" => has_organizations
            }; options...)
end

bill_list(auth::String, jurisdiction, session; include_organizations = true, has_organizations = true, options...) =
    bill_list(authenticate(auth), jurisdiction, session; include_organizations = include_organizations, has_organizations = has_organizations, options...)

bill_list(jurisdiction, session; auth = authenticate(), include_organizations = true, has_organizations = true, options...) =
    bill_list(auth, jurisdiction, session; include_organizations = include_organizations, has_organizations = has_organizations, options...)

# -------

function organization_positions(auth::Auth, jurisdiction, org_id; options...)
    maplight_get(auth, "/services_open_api/map.organization_positions_v1.json", {
        "jurisdiction" => jurisdiction,
        "organization_id" => org_id }; options...)
end

organization_positions(auth::String, jurisdiction, org_id; options...) = organization_positions(authenticate(auth), jurisdiction, org_id; options...)

organization_positions(jurisdiction, org_id; auth = authenticate(), options...) = organization_positions(auth, jurisdiction, org_id; options...)

# -------

function organization_search(auth::Auth, search_str, exact = false; options...)
    maplight_get(auth, "/services_open_api/map.organization_search_v1.json", {
        "search" => search_str,
        "exact" => int(exact) }; options...)
end

organization_search(auth::String, search_str, exact = false; options...) = organization_search(authenticate(auth), search_str, exact; options...)

organization_search(search_str, exact = false; auth = authenticate(), options...) = organization_search(auth, search_str, exact; options...)


text:
using BioSymbols
using BioSequences
using Base.Test

using Rifraf

import Rifraf.Proposal,
       Rifraf.Substitution,
       Rifraf.Insertion,
       Rifraf.Deletion,
       Rifraf.apply_proposals

srand(1)

@testset "Proposals" begin
    @testset "apply proposals" begin
        template = DNASeq("ACG")
        proposals = Proposal[Insertion(0, DNA_T),
                             Substitution(2, DNA_T),
                             Insertion(3, DNA_A),
                             Deletion(3)]
        expected = DNASeq("TATA")
        result = apply_proposals(template, proposals)
        @test result == expected
    end

    @testset "apply insertion" begin
        template = DNASeq("ATT")
        proposals = Proposal[Rifraf.Insertion(2,DNA_G)]
        expected = DNASeq("ATGT")
        result = apply_proposals(template, proposals)
        @test result == expected
    end

    @testset "apply ambiguous insertions" begin
        template = DNASeq("ACG")
        proposals = Proposal[Insertion(0, DNA_C),
                             Insertion(0, DNA_T)]
        @test_throws Rifraf.AmbiguousProposalsError apply_proposals(template, proposals)
    end

    @testset "apply ambiguous proposals" begin
        template = DNASeq("ACG")
        proposals = Proposal[Deletion(1),
                             Substitution(1, DNA_T)]
        @test_throws Rifraf.AmbiguousProposalsError apply_proposals(template, proposals)
    end
end

ext: jl


text:
import JuMP
@testset "Selection" begin

    # Helper methods
    nt(name) = named(name)
    tf(name) = t -> nt(name)(t)
    iv(size, name="in") = inputvertex(name, size)
    av(in, outsize, name) = absorbvertex(IndMem(MatMul(nout(in), outsize)), in; traitdecoration=tf(name))
    tv(in, name) = invariantvertex(IndMem(identity, in, nout(in)), in; traitdecoration=tf(name))

    concindmem(ins) = f ->  IndMem(f, ins, sum(nout, ins))
    cc(ins...; name) = conc(ins...; dims=1, traitdecoration=tf(name), outwrap=concindmem(ins))
    nc(name) = traitconf(nt(name))

    @testset "Split relaxed exact" begin
        import NaiveNASlib: split_exact_relaxed, Relaxed, Exact

        @test (:a => relaxed(3)) == (:a => 3 => Relaxed())

        @test split_exact_relaxed((:a=>2, :b=>3=>Relaxed(), :c=>4=>Exact()),) == (Dict(:a => 2, :c => 4), Dict(:b => 3))
        @test split_exact_relaxed((Dict(:a=>2, :b=>3=>Relaxed(), :c=>4=>Exact())),) == (Dict(:a => 2, :c => 4), Dict(:b => 3))
    end

    @testset "Parselect" begin
        import NaiveNASlib: parselect
        mat = reshape(collect(1:3*4), 3, 4)

        @test parselect(mat, 2 => [2,-1,-1,4]; newfun = (args...) -> 0) == [4 0 0 10;5 0 0 11;6 0 0 12]
        @test parselect(mat, 1 => [-1,1,3]; newfun = (args...) -> 0) == [0 0 0 0;1 4 7 10;3 6 9 12]

        @test parselect(mat, 1 => [2,-1,3,-1], 2 => [-1,1,-1,4]; newfun = (args...) -> 0) == [0 2 0 11;0 0 0 0;0 3 0 12;0 0 0 0]

        dfun = (T, d, s...) ->  d == 1 ? 10 : -10

        @test parselect(ones(2,2), 1 => [-1, 1], 2 => [1 , -1], newfun=dfun) == [10 -10; 1 -10]
        @test parselect(ones(2,2), 2 => [1 , -1], 1 => [-1, 1], newfun=dfun) == [10 10; 1 -10]
    end

    @testset "Make ΔNout" begin
        import NaiveNASlib: Exact, Relaxed, fallback
        @testset "All Exact" begin
            res = ΔNout(:a => 2, :b => 3 => Exact())
            @test res isa ΔNout{Exact}
            @test res.Δs == Dict(:a => 2, :b => 3)

            @test base(fallback(res)) isa ΔNout{Relaxed}
        end

        @testset "All Relaxed" begin
            import NaiveNASlib: Exact, Relaxed
            res = ΔNout(:a => 2 => Relaxed(), :b => relaxed(3))
            @test res isa ΔNout{Relaxed}
            @test res.Δs == Dict(:a => 2, :b => 3)

            @test base(fallback(res)) isa ΔSizeFailNoOp
        end

        @testset "Mix" begin
            import NaiveNASlib: ΔNoutMix
            res = ΔNout(:a => 2, :b => relaxed(3), :c => 4, :d => relaxed(5))
            @test res isa ΔNoutMix
            @test res.exact.Δs == Dict(:a => 2, :c => 4)
            @test res.relax.Δs == Dict(:b => 3, :d => 5)

            @test base(fallback(res)) isa ΔNout{Relaxed}
        end

        @testset "Fallback" begin
            @test fallback(ΔNout(:a => 2; fallback=ThrowΔSizeFailError())) isa ThrowΔSizeFailError
            @test fallback(ΔNout(:a => relaxed(2); fallback=ThrowΔSizeFailError())) isa ThrowΔSizeFailError
            @test fallback(ΔNout(:a => relaxed(2), :b => 3; fallback=ThrowΔSizeFailError())) isa ThrowΔSizeFailError
        end
    end

    @testset "Join extraced inds" begin
        import NaiveNASlib: join_extracted_inds

        @test join_extracted_inds([2,3,4,5], [0,0,0,0]) == [2,3,4,5]
        @test join_extracted_inds([2,3,4,5], [1,0,0,0,0,0,1]) == [2,-1,3,4,5,-1]
        @test join_extracted_inds([2,3,4,5], [1,2,3,4]) == [2,-1,3,-1,-1,4,-1,-1,-1,5,-1,-1,-1,-1,]
        @test join_extracted_inds([2], [2,4,0]) == [2,-1,-1,-1,-1,-1,-1]
        @test join_extracted_inds([], [3,0]) == [-1,-1,-1]
        @test join_extracted_inds([1,2,4], [0]) == [1,2,4]
        @test join_extracted_inds([1,2,3,4,8], [0,1,2]) == [1, 2, -1, 3, -1, -1, 4, 8]
    end

    @testset "Wrap and scale utility function" begin
        using NaiveNASlib: remaputilityfun
        inpt = iv(8)
        v1 = av(inpt, 9, "v1")
        v2 = av(v1, 7, "v2")
        v3 = av(v2, 8, "v3")

        @testset "Empty valid set" begin
            remapped = remaputilityfun(v -> 1e8, DefaultJuMPΔSizeStrategy(), [inpt], 1e-2 ,1e5)
            @test remapped(inpt) == 1e8
        end

        @testset "All zero" begin
            remapped = remaputilityfun(v -> 0, DefaultJuMPΔSizeStrategy(), [inpt, v1], 1e-2 ,1e5)
            @test remapped(inpt) == 0
            @test remapped(v1) == 0
        end

        @testset "No scaling" begin
            remapped = remaputilityfun(v -> v == v2 ? 0.0 : ones(nout(v)), DefaultJuMPΔSizeStrategy(), [inpt, v1, v2, v3], 1e-2, 1e5)

            @test remapped(inpt) == ones(nout(inpt))
            @test remapped(v1) == ones(nout(v1))
            @test remapped(v2) == 0.0
            @test remapped(v3) == ones(nout(v3))
        end

        @testset "Scale up" begin
            remapped = remaputilityfun(v -> 2.0 .^ -(1:nout(v)), DefaultJuMPΔSizeStrategy(), [inpt, v1, v2, v3], 1e-2, 1e5)

            @test remapped(inpt) ≈ 2.0 .^ -(1:nout(inpt))
            @test remapped(v1) ≈ [2.56, 1.28, 0.64, 0.32, 0.16, 0.08, 0.04, 0.02, 0.01]
            @test remapped(v2) ≈ [2.56, 1.28, 0.64, 0.32, 0.16, 0.08, 0.04]
            @test remapped(v3) ≈ [2.56, 1.28, 0.64, 0.32, 0.16, 0.08, 0.04, 0.02]
        end
        @testset "Scale down" begin
            remapped = remaputilityfun(v -> 10.0 .^ (1:nout(v)), DefaultJuMPΔSizeStrategy(), [inpt, v1, v2, v3], 1e-2, 1e5)

            @test remapped(inpt) ≈ 10.0 .^ (1:nout(inpt))
            @test remapped(v1) ≈ [0.01, 0.1, 1.0, 10.0, 100.0, 1000.0, 10000.0, 100000.0, 1.0e6]
            @test remapped(v2) ≈ [0.01, 0.1, 1.0, 10.0, 100.0, 1000.0, 10000.0]
            @test remapped(v3) ≈ [0.01, 0.1, 1.0, 10.0, 100.0, 1000.0, 10000.0, 100000.0]
        end

        @testset "Saturate upper" begin
            # We want to scale up alot, but we hit the upper bound
            remapped = remaputilityfun(DefaultJuMPΔSizeStrategy(), [inpt, v1, v2, v3], 1e-2, 1e5) do v
                v == inpt && return fill(100, nout(v))
                v == v1 && return 10.0 .^-(1:nout(v))
                v == v2 && return (1:nout(v)) .*  1e-3
                v == v3 && return fill(10, nout(v))
            end

            @test remapped(inpt) == fill(100, nout(inpt))
            @test remapped(v1) ≈ 10.0 .^ (3:-1:-5)
            @test remapped(v2) ≈ 10.0:10.0:70.0
            @test remapped(v3) ≈ fill(1e5, nout(v3))
        end

        @testset "Saturate lower" begin
            # We want to scale down alot, but we hit the upper bound
            remapped = remaputilityfun(DefaultJuMPΔSizeStrategy(), [inpt, v1, v2, v3], 1e-2, 1e5) do v
                v == inpt && return fill(100, nout(v))
                v == v1 && return 10.0 .^(1:nout(v))
                v == v2 && return 1:nout(v)
                v == v3 && return fill(10, nout(v))
            end

            @test remapped(inpt) ≈ fill(100, nout(inpt))
            @test remapped(v1) ≈ [0.1, 1.0, 10.0, 100.0, 1000.0, 10000.0, 100000.0, 1.0e6, 1.0e7]
            @test remapped(v2) ≈ 0.01:0.01:0.07
            @test remapped(v3) ≈ fill(0.1, nout(v3))
        end

        @testset "Saturate both" begin
            # We want to scale both up and down but we hit both the upper and lower bound
            remapped = remaputilityfun(DefaultJuMPΔSizeStrategy(), [inpt, v1, v2, v3], 1e-2, 1e5) do v
                v == inpt && return 13
                v == v1 && return 1e6
                v == v2 && return 1e-3
                v == v3 && return fill(10, nout(v))
            end

            @test remapped(inpt) == 13
            @test remapped(v1) ≈ 1e6
            @test remapped(v2) ≈ 1e-3
            @test remapped(v3) == fill(10, nout(v3))
        end

        @testset "Mixed signs" begin
            remapped = remaputilityfun(DefaultJuMPΔSizeStrategy(), [inpt, v1, v2, v3], 1e-2, 1e5) do v
                v == inpt && return 13
                v == v1 && return 1e3 .* repeat([1, -1], nout(v))
                v == v2 && return 1e-3 .* repeat([1, -1], nout(v))
                v == v3 && return fill(10, nout(v))
            end

            @test remapped(inpt) == 13
            @test remapped(v1) ≈ 1e4 .* repeat([1, -1], nout(v1))
            @test remapped(v2) ≈ 1e-2 .* repeat([1, -1], nout(v2))
            @test remapped(v3) == fill(100, nout(v3))
        end
    end

    @testset "sizeobjective" begin

        using NaiveNASlib: selectmodel, createselectvars!, sizeobjective!

        function solvesize(utilityfun, vs; case = NaiveNASlib.NeuronIndices(), strat= DefaultJuMPΔSizeStrategy())
            model = selectmodel(case, strat, vs, utilityfun)

            outselectvars, outinsertvars, noutdict = createselectvars!(case, strat, vs, (;model, utilityfun))
            JuMP.set_upper_bound.(Iterators.flatten(values(outinsertvars)), 20)

            objexpr = sizeobjective!(case, strat, vs, (;model, outselectvars, outinsertvars, noutdict, utilityfun))

            JuMP.@objective(model, Max, objexpr)
            JuMP.optimize!(model)
            return noutdict
        end

        @testset "all positive utility $label" for (label, utilfun) in (
            ("scalar", v -> 1.0),
            ("array", v -> fill(1.0, nout(v))),
        )
            inpt = iv(3)
            v1 = av(inpt, 5, "v1")
            v2 = av(v1, 6, "v2")

            vs = all_in_graph(inpt)
            noutdict = solvesize(utilfun, vs)

            @testset "Size of $(name(v))" for v in vs
                @test JuMP.value(noutdict[v]) == nout(v)
            end
        end

        @testset "all negative utility $label" for (label, utilfun) in (
            ("scalar", v -> -1.0),
            ("array", v -> fill(-1.0, nout(v))),
        )

            inpt = iv(3)
            v1 = av(inpt, 5, "v1")
            v2 = av(v1, 6, "v2")

            vs = all_in_graph(inpt)
            noutdict = solvesize(utilfun, vs)

            @testset "Size of $(name(v))" for v in vs
                @test JuMP.value(noutdict[v]) == 1
            end
        end

        @testset "Want half size" begin
            inpt = iv(4)
            v1 = av(inpt, 6, "v1")
            v2 = av(v1, 8, "v2")

            vs = all_in_graph(inpt)
            noutdict = solvesize(vs) do v
                repeat([1, -1], nout(v) ÷ 2)
            end

            @testset "Size of $(name(v))" for v in vs
                @test JuMP.value(noutdict[v]) == nout(v) ÷ 2
            end
        end
    end

    @testset "Absorb 2 Absorb" begin
        inpt = iv(3)
        v1 = av(inpt, 5, "v1")
        v2 = av(v1, 4, "v2")

        g = CompGraph(inpt, v2)

        @test Δnout!(v -> 1:nout(v), v1, -2)

        @test lastouts(v1) == lastins(v2) == [3,4,5]
        @test size(g(ones(3))) == (nout(v2),)

        @test Δnout!(v1, 3)

        @test lastouts(v1) == lastins(v2) == [1,2,3,-1,-1,-1]
        @test size(g(ones(3))) == (nout(v2),)

        @test Δnin!(v -> 1:nout(v), v2, -2)
        @test lastouts(v1) == lastins(v2) == [3,4,5,6]
        @test size(g(ones(3))) == (nout(v2),)
    end

    @testset "Argument plumbing" begin
        function graphgen()
            inpt = iv(3)
            v0 = av(inpt, 3, "v0")
            v1 = av(v0, 5, "v1")
            v2 = av(v1, 4, "v2")
            return v1, v2
        end

        @testset "Δnout!$(use_fun ? " with utility function" : "")" for use_fun in (false, true)
            f = use_fun ? (v -> (1:nout(v)),) : ()
            @testset "Single Δnout!$(wrap === identity ? "" : " $wrap")" for wrap in (identity, relaxed)
                v1,v2 = graphgen()
                @test Δnout!(f..., v1, wrap(3))
                @test nout.((v1, v2)) == (8, 4)
            end

            @testset "Single Δnout! pair$(wrap === identity ? "" : " $wrap")" for wrap in (identity, relaxed)
                v1,v2 = graphgen()
                @test Δnout!(f..., v1=>wrap(3))
                @test nout.((v1, v2)) == (8, 4)
            end

            @testset "Single Δnout! Dict$(wrap === identity ? "" : " $wrap")" for wrap in (identity, relaxed)
                v1,v2 = graphgen()
                @test Δnout!(f..., Dict(v1=>wrap(3)))
                @test nout.((v1, v2)) == (8, 4)
            end

            @testset "Multi Δnout! pair$(wrap === identity ? "" : " $wrap")" for wrap in (identity, relaxed)
                v1,v2 = graphgen()
                @test Δnout!(f..., v1=>wrap(3), v2=>wrap(2))
                @test nout.((v1, v2)) == (8, 6)
            end

            @testset "Multi Δnout! Dict$(wrap === identity ? "" : " $wrap")" for wrap in (identity, relaxed)
                v1,v2 = graphgen()
                @test Δnout!(f..., Dict(v1=>wrap(3), v2=>wrap(2)))
                @test nout.((v1, v2)) == (8, 6)
            end

            @testset "Multi Δnout! pair mix" begin
                v1,v2 = graphgen()
                @test Δnout!(f..., v1=>3, v2=>relaxed(2))
                @test nout.((v1, v2)) == (8, 6)
            end

            @testset "Multi Δnout! Dict mix" begin
                v1,v2 = graphgen()
                @test Δnout!(f..., Dict(v1=>3, v2=>relaxed(2)))
                @test nout.((v1, v2)) == (8, 6)
            end
        end

        @testset "Δnin!$(use_fun ? " with utility function" : "")" for use_fun in (false, true)
            f = use_fun ? (v -> (1:nout(v)),) : ()
            @testset "Single Δnin!$(wrap === identity ? "" : " $wrap")" for wrap in (identity, relaxed)
                v1,v2 = graphgen()
                @test Δnin!(f..., v1, wrap(3))
                @test nin.((v1, v2)) == ([6], [5])
            end

            @testset "Single Δnin! pair$(wrap === identity ? "" : " $(repr(wrap))")" for wrap in
            (identity, tuple, relaxed, relaxed ∘ tuple, tuple ∘ relaxed)
                v1,v2 = graphgen()
                @test Δnin!(f..., v1=>wrap(3))
                @test nin.((v1, v2)) == ([6], [5])
            end

            @testset "Single Δnin! Dict$(wrap === identity ? "" : " $(repr(wrap))")" for wrap in
            (identity, tuple, relaxed, relaxed ∘ tuple, tuple ∘ relaxed)
                v1,v2 = graphgen()
                @test Δnin!(f..., Dict(v1=>wrap(3)))
                @test nin.((v1, v2)) == ([6], [5])
            end

            @testset "Multi Δnin! pair mix $(repr(wrap1)) and $(repr(wrap2)) " for wrap1 in
            (identity, tuple, relaxed, relaxed ∘ tuple, tuple ∘ relaxed), wrap2 in
            (identity, tuple, relaxed, relaxed ∘ tuple, tuple ∘ relaxed)
                v1,v2 = graphgen()
                @test Δnin!(f..., v1=>wrap1(3), v2=>wrap2(2))
                @test nin.((v1, v2)) == ([6], [7])
            end

            @testset "Multi Δnin! Dict $(repr(wrap1)) and $(repr(wrap2)) " for wrap1 in
                (identity, tuple, relaxed, relaxed ∘ tuple, tuple ∘ relaxed), wrap2 in
                (identity, tuple, relaxed, relaxed ∘ tuple, tuple ∘ relaxed)
                v1,v2 = graphgen()
                @test Δnin!(f..., Dict(v1=>wrap1(3), v2=>wrap2(2)))
                @test nin.((v1, v2)) == ([6], [7])
            end
        end
    end

    @testset "add_participants" begin
        using NaiveNASlib: add_participants!
        using Logging: Info

        inpt = iv(3)
        v1 = av(inpt, 4, "v1")
        v2 = av(v1, 4, "v2")
        v3 = "v3" >> v1 + v2

        @testset "ΔNout $wraplabel" for (wraplabel, wrap) in (
            (LogΔSizeExec, s -> LogΔSizeExec("", Info, s)),
            (AlignNinToNout, s -> AlignNinToNout(vstrat=s)),
            (TruncateInIndsToValid, TruncateInIndsToValid),
            (WithUtilityFun, s -> WithUtilityFun(identity, s)),
            (TimeLimitΔSizeStrategy, s -> TimeLimitΔSizeStrategy(10, s)),
            (TimeOutAction,  s -> TimeOutAction(base=s)),
        )
            @test add_participants!(wrap(ΔNout(v2 => 3))) == [v2, v3 ,v1]
        end

        v4 = av(inpt, 5, "v4")
        @testset "ΔNoutMix" begin
            @test add_participants!(ΔNout(v1 => 2, v4 => relaxed(3), v2 => relaxed(2))) == [v1, v2, v3, v4]
        end
    end

    @testset "Absorb 2 Absorb fail error" begin
        import NaiveNASlib: ΔSizeFailError
        inpt = iv(3)
        v1 = av(inpt, 5, "v1")
        v2 = av(v1, 4, "v2")

        @test_throws ΔSizeFailError Δsize!(v -> 1:nout(v), ThrowΔSizeFailError("Success!?"), v1)
    end

    @testset "Absorb 2 Absorb fail NoOp" begin
        inpt = iv(3)
        v1 = av(inpt, 5, "v1")
        v2 = av(v1, 4, "v2")

        @test Δsize!(v -> 1:nout(v), ΔSizeFailNoOp(), v1) == false
    end

    @testset "Fail impossible" begin
        v0 = iv(3)
        v1 = tv(v0, "v1")
        v2 = av(v0, 3, "v2")
        v3 = "v3" >> v2 + v1
        @test Δnout!(v0, 1) == false
        @test @test_logs(
                (:warn, "Could not change nout of v3 by 1! Relaxing constraints..."),
                (:warn, "Could not change nout of v3 by 1! Attempt aborted!"),
                Δnout!(v3, 1)) == false
        @test @test_logs(
                (:warn, "Could not change nin of v3 by 1 and 1! Relaxing constraints..."),
                (:warn, "Could not change nin of v3 by 1 and 1! Attempt aborted!"),
                 Δnin!(v3, 1, 1)) == false
    end

    @testset "Δ = $wrap(0)" for wrap in (identity, relaxed)
        inpt = iv(3)
        v1 = av(inpt, 5, "v1")
        v2 = av(v1, 4, "v2")

        @test Δnout!(v1 => wrap(0)) == true
        @test [nout(v1)] == nin(v2) == [5]
    end

    @testset "SizeStack duplicate" begin
        inpt = iv(3)
        v1 = av(inpt, 7, "v1")
        v2 = av(inpt, 4, "v2")
        v3 = cc(v2, v1, name="v3")
        v4 = cc(v3, v2, name="v4")

        g = CompGraph(inpt, v4)
        @test size(g(ones(3,2))) == (nout(v4), 2)

        @test Δnout!(v4, -4)

        @test nout(v1) == 7
        @test nout(v2) == 2

        @test size(g(ones(3))) == (nout(v4),)

        @test Δnout!(v4, 6)

        @test nout(v1) == 11
        @test nout(v2) == 3

        @test size(g(ones(3, 3))) == (nout(v4), 3)

        @test Δnin!(v4, relaxed(-2), 3)

        @test nout(v1) == 6
        @test nout(v2) == 6

        @test size(g(ones(3))) == (nout(v4),)

        @test Δnin!(v4, -2, missing)

        @test nout(v1) == 5
        @test nout(v2) == 5

        @test size(g(ones(3))) == (nout(v4),)
    end

    @testset "SizeInvariant duplicate" begin
        inpt = iv(3)
        v1 = av(inpt, 7, "v1")
        v2 = av(inpt, 4, "v2")
        v3 = av(inpt, 3, "v3")

        v4 = cc(v1, v2, name="v4")
        v5 = cc(v2, v3, v2, name="v5")

        v6 = nc("v6") >> v4 + v5

        g = CompGraph(inpt, v6)
        @test size(g(ones(3))) == (nout(v6),)

        @test Δnout!(v->1:nout(v), v6, -4)

        @test nout(v1) == 5
        @test nout(v2) == 2
        @test nout(v3) == 3

        @test size(g(ones(3))) == (nout(v6),)

        @test Δnout!(v6, 6)

        @test nout(v1) == 9
        @test nout(v2) == 4
        @test nout(v3) == 5

        @test size(g(ones(3))) == (nout(v6),)

        @test Δnin!(v -> 1:nout(v), v6, -2, relaxed(-2))

        @test nout(v1) == 8
        @test nout(v2) == 3
        @test nout(v3) == 5

        @test size(g(ones(3))) == (nout(v6),)

        @test Δnin!(v -> 1:nout(v), v6, 3, missing)

        @test nout(v1) == 10
        @test nout(v2) == 4
        @test nout(v3) == 6

        @test size(g(ones(3))) == (nout(v6),)
    end

    @testset "SizeStack exact infeasible" begin
        inpt = iv(3)
        v1 = av(inpt, 4, "v1")
        v2a = tv(v1, "v2a")
        v2b = tv(v1, "v2b")
        v2c = tv(v1, "v2c")
        v3 = cc(v2a, v2b, v2c; name="v3")

        g = CompGraph(inpt, v3)
        @test size(g(ones(3))) == (nout(v3),)

        @test @test_logs (:warn, r"Could not change nout of") Δnout!(v3, 2)

        @test nout(v3) == 3nout(v1) == 15

        @test size(g(ones(3))) == (nout(v3),)
    end

    @testset "SizeInvariant exact Δnout! infeasible" begin
        inpt = iv(3)
        v1 = av(inpt, 4, "v1")
        v2 = av(v1, 4, "v2")
        v3 = "v3" >> v1 + v2

        @test @test_logs (:warn, r"Could not change nout of") Δnout!(v1 => -1, v2 =>0)
        @test nout(v3) == nout(v2) == nout(v1) == 3
    end

    @testset "SizeInvariant exact Δnin! infeasible" begin
        inpt = iv(3)
        v1 = av(inpt, 4, "v1")
        v2 = av(v1, 4, "v2")
        v3 = "v3" >> v1 + v2

        @test @test_logs (:warn, r"Could not change nin of") Δnin!(v3 => (-1, 0))
        @test nout(v3) == nout(v2) == nout(v1) == 3
    end

    @testset "SizeStack one immutable" begin
        inpt = iv(3)
        v1 = av(inpt, 5, "v1")
        v2 = cc(inpt, v1, name="v2")

        g = CompGraph(inpt, v2)
        @test size(g(ones(3))) == (nout(v2),)

        @test Δnout!(v1, -3) do v
            # "Tempt" optimizer to not select inputs from inpt
            v === v2 && return (-nout(inpt):nout(v1)-1)
            return 1:nout(v)
        end

        @test nin(v2) == [nout(inpt), nout(v1)] == [3, 2]
        @test nout(v2) == 5

        @test size(g(ones(3))) == (nout(v2),)
    end

    @testset "Decrease SizeStack and SizeInvariant" begin
        inpt = iv(3)
        v1 = av(inpt, 10, "v1")
        v2 = av(inpt, 5, "v2")
        v3 = av(inpt, 10, "v3")
        v4 = av(inpt, 5, "v4")

        v5 = cc(v1, v2, v3, name="v5")
        v6 = cc(v2, v1, v2, v4, name="v6")

        v7 = nc("v7") >> v5 + v6

        g = CompGraph(inpt, v7)
        @test size(g(ones(3))) == (nout(v7),)

        @test Δnout!(v7, -7)

        @test nout(v1) == 6
        @test nout(v2) == 3
        @test nout(v3) == 9
        @test nout(v4) == 6

        @test size(g(ones(3))) == (nout(v7),)
    end

    @testset "Decorating strategy" begin
        import NaiveNASlib: DecoratingJuMPΔSizeStrategy

        struct DummyDecorator{S} <: DecoratingJuMPΔSizeStrategy
            s::S
        end
        NaiveNASlib.base(s::DummyDecorator) = s.s
        NaiveNASlib.add_participants!(s::DummyDecorator, vs=AbstractVertex[]) = NaiveNASlib.add_participants!(base(s), vs)

        function graphgen()
            inpt = iv(3)
            v1 = av(inpt, 10, "v1")
            v2 = av(inpt, 5, "v2")
            v3 = av(inpt, 10, "v3")
            v4 = av(inpt, 5, "v4")

            v5 = cc(v1, v2, v3, name="v5")
            v6 = cc(v2, v1, v2, v4, name="v6")

            v7 = nc("v7") >> v5 + v6

            CompGraph(inpt, v7)
        end

        genstrat(::Type{ΔNout{Exact}}, g) = ΔNoutExact(g.outputs[1], -7)
        genstrat(::Type{ΔNout{Relaxed}}, g) = ΔNoutRelaxed(g.outputs[1], -7)
        genstrat(::typeof(ΔNinExact), g) = ΔNinExact(g.outputs[1], (-7, missing))

        @testset "$basestrat" for basestrat in (
            ΔNout{Exact},
            ΔNout{Relaxed},
            ΔNinExact,
        )
            g_exp, g_act = graphgen(), graphgen()
            s_exp = genstrat(basestrat, g_exp)
            s_act = genstrat(basestrat, g_act)

            @test Δsize!(v -> 1:nout(v), s_exp)
            @test Δsize!(v -> 1:nout(v), DummyDecorator(s_act))

            @test nout.(vertices(g_exp)) == nout.(vertices(g_act))
        end
    end

    @testset "Increase SizeStack and SizeInvariant" begin
        inpt = iv(3)
        v1 = av(inpt, 3, "v1")
        v2 = av(inpt, 4, "v2")
        v3 = av(inpt, 4, "v3")
        v4 = av(inpt, 5, "v4")

        v5 = cc(v1, v2, v4, name="v5")
        v6 = cc(v3, v4, v1, name="v6")

        v7 = nc("v7") >> v5 + v6

        g = CompGraph(inpt, v7)
        @test size(g(ones(3))) == (nout(v7),)

        oldg = deepcopy(g)

        Δnout!(v7, 6)

        @test nout(v1) == 5
        @test nout(v2) == 6
        @test nout(v3) == 6
        @test nout(v4) == 7

        @test size(g(ones(3))) == (nout(v7),)

        # Since we are just increasing we have the same mapping except a few zeros
        @test filter(>(0), g(1:3)) == oldg(1:3)
    end

    @testset "SizeStack increase decrease" begin
        inpt = iv(3)
        v1 = av(inpt, 3, "v1")
        v2 = av(inpt, 10, "v2")
        v3 = cc(v1,v2, name="v3")
        v4 = av(v3, 4, "v4")

        g = CompGraph(inpt, v3)
        @test size(g(ones(3))) == (nout(v3),)

        @test Δnout!(v -> 1:nout(v), v3 => relaxed(-5), v1 => 2, v2 => relaxed(0))

        @test nin(v3) == [nout(v1), nout(v2)] == [5, 7]
        @test nout(v3) == sum(nin(v3)) == 12

        @test lastins(v3) == [lastouts(v1), lastouts(v2)] == [[1,2,3,-1,-1], [4,5,6,7,8,9,10]]
        @test lastouts(v3) == lastins(v4) == [1,2,3,-1,-1,7,8,9,10,11,12,13]

        @test size(g(ones(3))) == (nout(v3),)
    end

    @testset "Constrained by remote subtree" begin
        inpt = iv(3)
        v0 = av(inpt, 10, "v0")

        # Main branch, the one we want to change
        v1 = cc(v0, v0, name="v1")

        # Path to subtree
        # First hop
        v2 = av(inpt, 10, "v2")
        v3 = nc("v3") >> v2 + v0

        # Second hop (v4a and v4b should not be touhed)
        v4a = av(inpt, 3, "v4a")
        v4b = av(inpt, 5, "v4b")
        v5 = cc(v4a, v3, v4b, name="v5")

        # Subtree
        v6 = av(inpt, 4, "v6")
        v7 = av(inpt, nout(v5) - nout(v6), "v7")
        v8 = cc(v6,v7,name="v8")

        # Aaaand connect it to the path
        v9 = nc("v9") >> v8 + v5
        v10 = av(v9, 5, "v10")

        g = CompGraph(inpt, [v1, v9])
        @test size.(g(ones(3))) == ((nout(v1),), (nout(v9),))

        Δnout!(v2, -1)

        @test nout(v6) == 4
        @test nout(v7) == 13
        @test nout(v0) == 9

        @test size.(g(ones(3))) == ((nout(v1),), (nout(v9),))

        Δnout!(v2, 5)

        @test nout(v6) == 4
        @test nout(v7) == 18
        @test nout(v0) == 14

        @test size.(g(ones(3))) == ((nout(v1),), (nout(v9),))
    end

    @testset "Low utility in not-connected" begin
        inpt = iv(3)
        v1 = av(inpt, 4, "v1")
        v2 = av(v1, 5, "v2")
        v3 = av(v2, 3, "v3")
        v4 = av(v2, 5, "v4")

        @test Δnout!(v1 => 2) do v
            v == v2 ? (-2:nout(v2)-3) : 1
        end
        @test nout(v1) == 6
        @test nout(v2) == 5
    end

    @testset "Increase at vertex removal" begin
        inpt = iv(3)
        v1 = av(inpt, 2, "v1")
        v2 = tv(v1, "v2")
        v3 = av(v2, 4, "v3")
        v4 = av(v3, 3, "v4")

        g = CompGraph(inpt, v4)
        @test size(g(ones(3))) == (nout(v4),)

        @test remove!(v3, RemoveStrategy())

        @test lastins(v4)== [1,2,3,4]
        @test lastouts(v1) == lastouts(v2) == [1, 2, -1, -1]

        @test size(g(ones(3))) == (nout(v4),)
    end

    @testset "Decrease at vertex removal" begin
        inpt = iv(3)
        v1 = av(inpt, 2, "v1")
        v2 = tv(v1, "v2")
        v3 = av(v2, 4, "v3")
        v4 = av(v3, 3, "v4")

        g = CompGraph(inpt, v4)
        @test size(g(ones(3))) == (nout(v4),)

        @test remove!(v3, RemoveStrategy(DecreaseBigger()))

        # What happened now is that nin(v4) got decreased from 4 to 2.
        # However, there was absolutely no need at all to select anything from v2 and before as they have not changed.

        @test lastins(v4) == lastouts(v3) == [3, 4]
        @test lastouts(v1) == lastouts(v2) == [1, 2]

        @test size(g(ones(3))) == (nout(v4),)
    end

    @testset "Increase at vertex removal SizeStack" begin
        inpt = iv(3)
        v1 = av(inpt, 3, "v1")
        v2 = av(inpt, 4, "v2")
        v3 = av(v1, 5, "v3")
        v4 = cc(v2, v3, name="v4")

        g = CompGraph(inpt, v4)
        @test size(g(ones(3))) == (nout(v4),)

        @test remove!(v3, RemoveStrategy())

        # v4 did not change
        @test lastins(v4) == [[1, 2, 3, 4], [1, 2, 3, 4, 5]]
        @test lastouts(v4) == 1:9

        @test [lastouts(v2), lastouts(v1)] == [[1,2,3,4], [1,2,3,-1,-1]]

        @test size(g(ones(3))) == (nout(v4),)
    end

    @testset "Decrease at vertex removal SizeStack" begin
        inpt = iv(3)
        v1 = av(inpt, 3, "v1")
        v2 = av(inpt, 4, "v2")
        v3 = av(v1, 5, "v3")
        v4 = cc(v2, v3, name="v4")

        g = CompGraph(inpt, v4)
        @test size(g(ones(Float32, 3))) == (nout(v4),)

        @test remove!(v3, RemoveStrategy(DecreaseBigger(mapstrat=WithUtilityFun(v -> 1:nout(v)))))

        @test lastins(v4) == [[1, 2, 3, 4], [3, 4, 5]]
        @test lastouts(v4) == [1, 2, 3, 4, 7, 8, 9]

        # v2 and v1 did not change
        @test [lastouts(v2), lastouts(v1)] == [[1,2,3,4], [1,2,3]]

        @test size(g(ones(3))) == (nout(v4),)
    end

    @testset "TruncateInIndsToValid" begin
        inpt = iv(3)
        v1 = av(inpt, 3, "v1")
        v2 = tv(v1, "v2")
        v3 = av(v2, 2, "v3")

        g = CompGraph(inpt, v2)
        @test size(g(ones(3))) == (nout(v2),)

        vnew = av(inpt, 7, "vnew")
        create_edge!(vnew, v2;strategy=PostAlign(TruncateInIndsToValid(AlignNinToNout())))
        @test nin(v2) == [nout(v1), nout(vnew)] == [7, 7]

        @test lastouts(vnew) == 1:7
        @test lastins(v2) == [[1,2,3,-1,-1,-1,-1], [1,2,3,-1,-1,-1,-1]]
        @test lastouts(v2) == [1,2,3,-1,-1,-1,-1]
    end

    @testset "WithUtilityFun" begin
        inpt = iv(3)
        v1 = av(inpt, 6, "v1")

        @test Δsize!(WithUtilityFun(v -> v === v1 ? [-1, 1, -1, 1, -1, 1] : 1:nout(v), DefaultJuMPΔSizeStrategy()), v1)

        @test nout(v1) == 3
        @test lastouts(v1) == [2,4,6]
    end

    @testset "CompConstraint" begin
        struct CompConstraint end
        function NaiveNASlib.compconstraint!(::NaiveNASlib.NeuronIndices, ::AbstractJuMPΔSizeStrategy, ::CompConstraint, data)
            var = data.outselectvars[data.vertex];
            JuMP.@constraint(data.model, var[[1, 3]] .== 0)
        end
        NaiveNASlib.nout(::CompConstraint) = 4
        NaiveNASlib.nin(::CompConstraint) = 3
        ccouts = nothing
        NaiveNASlib.Δsize!(::CompConstraint, ins::AbstractVector, outs::AbstractVector) = ccouts=outs

        inpt = iv(3)
        v1 = absorbvertex(CompConstraint(), inpt, traitdecoration = tf("v1"))
        v2 = av(v1, 5, "v2")

        Δnout!(v1, 3)
        @test ccouts == lastins(v2) == [2, 4, -1, -1, -1, -1, -1]
    end

    @testset "Time limit" begin
        import NaiveNASlib: TimeLimitΔSizeStrategy, selectmodel, NeuronIndices
        @test JuMP.time_limit_sec(selectmodel(NeuronIndices(), TimeLimitΔSizeStrategy(123), AbstractVertex[], () -> 1)) == 123
    end

    @testset "Timeout action" begin
        import NaiveNASlib: TimeLimitΔSizeStrategy, TimeOutAction, ΔSizeFailNoOp

        v0 = av(iv(1), 1000, "v0")
        v1 = cc(v0, v0, v0; name="v1")
        v2 = cc(v1, v0, v1; name="v2")
        v3 = "v3" >> v2 + v2

        @test @test_logs (:warn, "Solver timed out") Δsize!(TimeLimitΔSizeStrategy(0.001, TimeOutAction(); fallback=ΔSizeFailNoOp()), v3) == false

        modelres = nothing
        action = m -> (modelres = m; return false)

        @test Δsize!(TimeLimitΔSizeStrategy(0.001, TimeOutAction(;action); fallback=ΔSizeFailNoOp()), v3) == false

        @test modelres isa JuMP.Model
    end

    @testset "All negative utility $label" for (label, utilfun) in (
        ("scalar", v -> -1.0),
        ("array", v -> fill(-1.0, nout(v))),
    )
        inpt = iv(3)
        v1 = av(inpt, 4, "v1")
        v2 = av(v1, 5, "v2")

        @test Δsize!(utilfun, CompGraph(inpt, v2)) == true

        @test [nout(v1)] == nin(v2) == [nout(v2)] == [1]
    end

    @testset "WithKwargs" begin
        mutable struct WantKwargs
            insize::Int
            outsize::Int
        end
        NaiveNASlib.nout(w::WantKwargs) = w.outsize
        NaiveNASlib.nin(w::WantKwargs) = [w.insize]

        seenkwargs = Dict()
        function NaiveNASlib.Δsize!(w::WantKwargs, ins::AbstractVector, outs::AbstractVector; kwargs...)
            if ins[] !== missing
                w.insize = length(ins[])
            end
            w.outsize = length(outs)
            merge!(seenkwargs, kwargs)
        end

        inpt = iv(3)
        v1 = absorbvertex(WantKwargs(3, 5), inpt; traitdecoration=tf("v1"))

        @test Δsize!(WithKwargs(ΔNout(v1 => 2); a=3, b=4)) == true

        @test nout(v1) == 7
        @test seenkwargs == Dict(:a => 3, :b => 4)

        @test Δsize!(validateafterΔsize(base=WithKwargs(ΔNout(v1 => relaxed(-2)); c = 5))) == true
        @test nout(v1) == 5
        @test seenkwargs == Dict(:a => 3, :b => 4, :c=>5)
    end

    @testset "UpperInsertBound" begin
        import Logging
        inpt = iv(3)
        v1 = av(inpt, 4, "v1")
        v2 = av(v1, 5, "v2")

        strategy = ΔNoutExact(v1=>2)
        fallback = LogΔSizeExec("Success!?", Logging.Warn)
        @test @test_logs (:warn, "Success!?") Δsize!(UpperInsertBound(0;strategy, fallback)) == false

        @test @test_logs (:warn, "Success!?") Δsize!(TimeOutAction(;base=UpperInsertBound(0;strategy), fallback)) == false

        @test Δsize!(UpperInsertBound([v2], 0;strategy, fallback)) == true

    end
end


text:
"""
    density(sp::AbstractSparseMatrix)

Compute the density of `sp`, the fraction of elements that
are non-zero. This assumes that there are no stored values equal to zero.
"""
function density(sp::AbstractSparseMatrix)
    return SparseArrays.nnz(sp) / length(sp)
end

"""
    summarystats(sp::SparseMatrixCSC)

Print some statistics for `sp`. The line
"number of non-zeros" counts structural non-zeros.
"""
function summarystats(sp::AbstractSparseMatrix)
    padding = 18
    for (description, statistic) in (
        ("size", size(sp)),
        ("num. elements", length(sp)),
        ("num. non-zeros", nnz(sp)),
        ("density", Printf.@sprintf("%.4e", density(sp))),
        ("num. stored zeros", count(iszero, SparseArrays.nonzeros(sp))))
        println(rpad(description, padding), ": ", statistic)
    end
    return nothing
end


text:
function get_map_bounds(routes_dict::OrderedDict{String,OrderedDict{Int64,Tuple{Float64,Float64}}})
    """
    Description:
        This function takes routes_dict and returns map_bounds for folium.
    Returns:
        (minimum(LAT),minimum(LOT)),(maximum(LAT),maximum(LOT))
    """
    lats = []
    lots = []

    for route in keys(routes_dict)
        [push!(lats,i[1]) for i in values(routes_dict[route])]
        [push!(lots,i[2]) for i in values(routes_dict[route])]
    end

    return (minimum(lats),minimum(lots)),(maximum(lats),maximum(lots))
end

function plot_folium_map(routes_filtered::DataFrame, routes_dict::OrderedDict{String,OrderedDict{Int64,Tuple{Float64,Float64}}})
	"""
	Description:
		Plots TTC routes on the folium interactive Toronto map.
	Example:
		plot_folium_map(routes_filtered, routes_dict)
	"""
    m = flm.Map()
    for route in keys(routes_dict)
        #Get route info
        route_name = values(routes_filtered[routes_filtered[!, "route_id"] .== parse(Int64,split(route,"_")[1]),:][!, "route_long_name"])[1]
        route_type = values(routes_filtered[routes_filtered[!, "route_id"] .== parse(Int64,split(route,"_")[1]),:][!, "route_type"])[1]

        # Mark route type and assign differents colors
        if route_type == 1
            route_type = "SUBWAY"
            color = "red"
        elseif route_type == 0
            route_type = "STREETCAR"
            color = "green"
        elseif route_type == 3
            route_type = "BUS"
            color = "gray"
        end

        # Drawing a route
        flm.PolyLine(
        [values(sort(routes_dict[route]))],
        popup="Route id: $route Type: $route_type",
        tooltip="Route name: $route_name",
        color=color,
        weight = ifelse(route_type == "SUBWAY", 6, 4),
        opacity = 0.6
    ).add_to(m)

    end

	# Reduce vis to the map bounds and draw rectangle as a border.
    MAP_BOUNDS = get_map_bounds(routes_dict)
    flm.Rectangle(MAP_BOUNDS, color="black",weight=6).add_to(m)
    m.fit_bounds(MAP_BOUNDS)
    return m
end

println("Folium plotting functions were read.")

max_stars_repo_name: NykPol/EpidemicInUrbanNetworkToronto

max_issues_repo_name: NykPol/EpidemicInUrbanNetworkToronto

max_forks_repo_name: NykPol/EpidemicInUrbanNetworkToronto


text:
__precompile__()

module Twistronics

using LinearAlgebra
using DelimitedFiles
using Statistics
using Dates

using Random
using SparseArrays
using Arpack

using PyCall
const plt=PyNULL()
const colormap=PyNULL()
function __init__()
    copy!(plt,pyimport_conda("matplotlib.pyplot","matplotlib"))
    copy!(colormap,pyimport_conda("matplotlib.colors","matplotlib"))
end

include("Geometry_twisted.jl")
include("tight-binding_twisted.jl")
include("quantities_twisted.jl")
include("valley_operator_twisted.jl")

end


text:
module ExcelDataStreams

# package code goes here
using DataStreams
using Taro
using DataFrames
using DataArrays

export ExcelDataStream

initialized = false

type ExcelDataStream <: Data.Source
    data::DataFrame
end

"""
`ExcelDataStream( filename::AbstractString; kwargs... )`

Constructor for type which gives access to an Excel file as a DataStream

### Arguments

* `filename::AbstractString` : the name of the Excel file to stream

### Keyword Arguments

*   `columnindices::Vector{Int}` -- The indices of columns to use as fields.  Defaults to empty.
*   `names::Vector{Symbol}` -- The names to use for the fields.  Defaults to empty.
*   `eltypes::Vector{DataTypes}` -- The types expected for each field.  Conversion will be attemped.  Defaults to empty.
*   `sheet::String` -- Which sheet of the Excel file contains the data of interest.  Defaults to `sheet1`.
*   `skipstart::Int` -- How many initial rows of the Excel file to skip.  Defaults to `0`.
*   `nastrings::Vector{String}` -- Strings to consider to be missing data.  Defaults to empty.
"""
function ExcelDataStream(
    filename::AbstractString,
    ;
    columnindices::Vector{Int} = Int[],
    names::Vector{Symbol} = String[],
    eltypes::Vector{DataType} = DataType[],
    sheet::String = "Sheet1",
    skipstart::Int = 0,
    nastrings::Vector{String} = String[],
)
    global initialized
    if !initialized
        Taro.init()
        initialized = true
    end
    excel = Workbook( filename )
    sheetdata = getSheet( excel, sheet )

    columns = [DataArray( eltype, 0 ) for eltype in eltypes]
    done = false
    rownum = skipstart
    while !done
        row = getRow( sheetdata, rownum )
        done = isnull(row)
        if !done
            for i = 1:length(names)
                cell = getCell( row, columnindices[i] )
                cellvalue = getCellValue( cell )
                if cellvalue == nothing || cellvalue in nastrings
                    push!( columns[i], missing )
                else
                    push!( columns[i], cellvalue )
                    done = false
                end
            end
            if done
                # the last row is all missing
                [pop!( column ) for column in columns]
            end
            rownum += 1
        end
    end

    df = DataFrame()
    for i = 1:length(names)
        df[names[i]] = columns[i]
    end
    return ExcelDataStream( df )
end

Data.schema( ed::ExcelDataStream ) =
    Data.Schema( eltypes(ed.data), string.(names(ed.data)), size(ed.data,1) )

Data.isdone( ed::ExcelDataStream, row::Int, col::Int ) = any((row,col).>size(ed.data))

Data.streamtype( ::Type{ExcelDataStream}, ::Type{Data.Field} ) = true

Data.streamfrom(ed::ExcelDataStream, ::Type{Data.Field}, ::Type{T}, row::Int, col::Int) where {T} =
    T(ed.data[row,col])

Data.accesspattern( ::ExcelDataStream ) = Data.RandomAccess

end # module


text:
using IterationControl
using Test

const IC = IterationControl

include("_models_for_testing.jl")
const SquareRooter = IC.SquareRooter

@testset "square_rooter.jl" begin
    include("square_rooter.jl")
end

@testset "utilities" begin
    include("utilities.jl")
end

# this test must happen before test of controls.jl:
@testset "api" begin
     include("api.jl")
end

@testset "controls" begin
    include("controls.jl")
end

@testset "stopping_controls" begin
    include("stopping_controls.jl")
end

@testset "composite_controls" begin
    include("composite_controls.jl")
end

@testset "wrapped_controls" begin
    include("wrapped_controls.jl")
end

@testset "train!" begin
    include("train.jl")
end

@testset "deprecated" begin
    include("deprecated.jl")
end


text:
using OpenAIGymAPI
using Base.Test

# define helper functions
function get_remote_base(host = "http://127.0.0.1", port = "5000")
  return(host * ":" * port)
end

# setup background server
function setup_background_server()
  @spawn run(detach(`python /tmp/gym-http-api/gym_http_server.py`))
end

# teardown background server
function teardown_background_server()
  client = GymClient(get_remote_base())
  shutdown_server(client)
end

# define wrapper
function with_server(fn)
  setup_background_server()
  sleep(3) # otherwise API is called before server is fully setup
  result = fn()
  sleep(3) # otherwise API is shutdown before function is done
  teardown_background_server()
  return(result)
end

########## TESTS ##########

##### Valid use cases #####
function test_create_destroy()
  client = GymClient(get_remote_base())
  instance_id = env_create(client, "CartPole-v0")
  @test haskey(env_list_all(client), instance_id)
  env_close(client, instance_id)
  @test !haskey(env_list_all(client), instance_id)
end
with_server(test_create_destroy)

function test_action_space_discrete()
  client = GymClient(get_remote_base())
  instance_id = env_create(client, "CartPole-v0")
  action_info = env_action_space_info(client, instance_id)
  @test action_info["name"] == "Discrete"
  @test action_info["n"] == 2
end
with_server(test_action_space_discrete)

function test_action_space_sample()
  client = GymClient(get_remote_base())
  instance_id = env_create(client, "CartPole-v0")
  action = env_action_space_sample(client, instance_id)
  @test 0 <= action < 2
end
with_server(test_action_space_sample)

function test_action_space_contains()
  client = GymClient(get_remote_base())
  instance_id = env_create(client, "CartPole-v0")
  action_info = env_action_space_info(client, instance_id)
  @test action_info["n"] == 2
  @test env_action_space_contains(client, instance_id, 0) == true
  @test env_action_space_contains(client, instance_id, 1) == true
  @test env_action_space_contains(client, instance_id, 2) == false
end
with_server(test_action_space_contains)

function test_observation_space_box()
  client = GymClient(get_remote_base())
  instance_id = env_create(client, "CartPole-v0")
  obs_info = env_observation_space_info(client, instance_id)
  @test obs_info["name"] == "Box"
  @test length(obs_info["shape"]) == 1
  @test obs_info["shape"] == [4]
  @test length(obs_info["low"]) == 4
  @test length(obs_info["high"]) == 4
end
with_server(test_observation_space_box)

function test_reset()
  client = GymClient(get_remote_base())
  instance_id = env_create(client, "CartPole-v0")
  init_obs = env_reset(client, instance_id)
  @test length(init_obs) == 4
end
with_server(test_reset)

function test_step()
  client = GymClient(get_remote_base())
  instance_id = env_create(client, "CartPole-v0")
  env_reset(client, instance_id)
  result = env_step(client, instance_id, 0)
  observation, reward = result["observation"], result["reward"]
  done, info = result["done"], result["info"]
  @test length(observation) == 4
end
with_server(test_step)

function test_monitor_start_close_upload()
  # @test os.environ.get("OPENAI_GYM_API_KEY")
  client = GymClient(get_remote_base())
  instance_id = env_create(client, "CartPole-v0")
  env_monitor_start(client, instance_id, "tmp", force = true)
  env_reset(client, instance_id)
  env_step(client, instance_id, 1)
  env_monitor_close(client, instance_id)
  # upload(client, "tmp")
end
with_server(test_monitor_start_close_upload)


text:
using Abaco
using JSON3
using BenchmarkTools

function nop(s, ne, name, value, inputs)
end

function onresult(ts, ne, name, value, inputs)
    @info "age [$ts]: scope: [$ne] $name = $value"
end

interval = 5
ages = 4
abaco = Abaco.init(nop, interval, ages)
formula(abaco, "y = x1 + x2")

function main(sns)
    for ts in 0:5:80
        for i in 1:100
            ingest(abaco, ts, sns[i].first, :x1, sns[i].second)
            ingest(abaco, ts, sns[i].first, :x2, sns[i].second)
        end
    end
end

sns = Dict(i => "SN$i" => float(i) for i in 1:100)

main(sns)
ingest(abaco, 85, "SN", "x1", 1.0)
abaco.oncomplete = onresult

# expect (0 allocations: 0 bytes)
@btime ingest(abaco, 85, "SN", "x2", 1.0)


text:
# © 2016 Massachusetts Institute of Technology.  See LICENSE file for details.
using ArgParse
s = ArgParseSettings()
@add_arg_table s begin
    "--gpu"
        help = "Use GPU (cuda)"
        action = :store_true
    "--profile"
        help = "Use profiler"
        action = :store_true
end
const args = parse_args(ARGS, s)
const gpu = args["gpu"]
const profile = args["profile"]

if gpu
    using CUDArt
    CUDArt.init(0)
    device(0)
    make_node(x...) = SeepNode(CudaArray(x...))
    get_data(x) = to_host(x)
else
    make_node(x...) = SeepNode(x...)
    get_data(x) = copy(x)
end

using Seep

const X = linspace(-1, 3, 100)
const Y = sin(X.*X)

const x = make_node(Float64, 1, 1)
const y = make_node(Float64, 1, 1)
ai = x

const v = Set{SeepNode}()

n1 = 100
const w1 = make_node(randn(n1, 1))
const b1 = make_node(zeros(n1, 1))
ai = tanh(w1*ai+b1)
link!(objective(1e-3*(w1.*w1)), v)

n2 = 100
const w2 = make_node(randn(n2, n1)/sqrt(n1))
const b2 = make_node(zeros(n2, 1))
ai = tanh(w2*ai+b2)
link!(objective(1e-3*(w2.*w2)), v)

const w3 = make_node(randn(1, n2)/sqrt(n2))
const b3 = make_node(zeros(1, 1))
ai = w3*ai+b3
link!(objective(1e-3*(w3.*w3)), v)

const e = ai-y
link!(objective(e.*e), v)

const a = ai

const E = Float64[]
;

fill!(w1.backward, 0)
fill!(b1.backward, 0)
fill!(w2.backward, 0)
fill!(b2.backward, 0)
fill!(w3.backward, 0)
fill!(b3.backward, 0)

update = momentum(1e-3, 0.8)

if gpu
  for n in v
    Seep.streamify!(n)
  end
end

open(f->graphviz(f, first(v)), "example.gv", "w")

n = 0

function step()
    global n

    for i in randperm(length(X))
        ni = n::Int + 1
        n = ni
        x[ni] = collect(X[i])
        y[ni] = Y[i:i]

        backward!(w1, ni)
        backward!(b1, ni)
        backward!(w2, ni)
        backward!(b2, ni)
        backward!(w3, ni)
        backward!(b3, ni)
    end

    push!(E, get_data(forward!(e, n))[1]^2)
    weight_update(update, w1, b1, w2, b2, w3, b3)
end

step()

for i in 1:10
    @time step()
end

Profile.clear_malloc_data()

for i in 1:50
    step()
end

kill!(e)

if gpu
CUDArt.close(0)
end

println("finished cleaning up")


text:
module Constraint

using OrderedCollections # for OrderedDict in Map

using MathOptInterface
const MOI = MathOptInterface
const MOIU = MOI.Utilities
const MOIB = MOI.Bridges

const CI = MOI.ConstraintIndex

# Definition of a constraint bridge
include("bridge.jl")

# Mapping between constraint indices and bridges
include("map.jl")

# Bridge optimizer bridging a specific constraint bridge
include("single_bridge_optimizer.jl")

# TODO(odow): the compiler in Julia <= 1.2 (and in later versions unless
# fixed) gets stuck compiling add_constrained_variable for some inputs. This
# method seemed necessary to fix it.
# See https://github.com/JuliaLang/julia/issues/32167 for more.
function MOI.Bridges.Constraint.bridge_constraint(BridgeType, b, f, s)
    return throw(MOI.UnsupportedConstraint{typeof(f),typeof(s)}())
end

# Constraint bridges
# Function conversion bridges
include("function_conversion.jl")
# Transformation between a set S1 and a set S2 such that A*S1 = S2 for some linear map A
include("set_map.jl")

include("vectorize.jl")
const Vectorize{T,OT<:MOI.ModelLike} =
    SingleBridgeOptimizer{VectorizeBridge{T},OT}
include("scalarize.jl")
const Scalarize{T,OT<:MOI.ModelLike} =
    SingleBridgeOptimizer{ScalarizeBridge{T},OT}
include("slack.jl")
const ScalarSlack{T,OT<:MOI.ModelLike} =
    SingleBridgeOptimizer{ScalarSlackBridge{T},OT}
const VectorSlack{T,OT<:MOI.ModelLike} =
    SingleBridgeOptimizer{VectorSlackBridge{T},OT}
include("interval.jl")
const SplitInterval{T,OT<:MOI.ModelLike} =
    SingleBridgeOptimizer{SplitIntervalBridge{T},OT}
include("quad_to_soc.jl")
const QuadtoSOC{T,OT<:MOI.ModelLike} =
    SingleBridgeOptimizer{QuadtoSOCBridge{T},OT}
include("soc_to_nonconvex_quad.jl") # do not add these bridges by default
const SOCtoNonConvexQuad{T,OT<:MOI.ModelLike} =
    SingleBridgeOptimizer{SOCtoNonConvexQuadBridge{T},OT}
const RSOCtoNonConvexQuad{T,OT<:MOI.ModelLike} =
    SingleBridgeOptimizer{RSOCtoNonConvexQuadBridge{T},OT}
include("norm_to_lp.jl")
const NormInfinity{T,OT<:MOI.ModelLike} =
    SingleBridgeOptimizer{NormInfinityBridge{T},OT}
const NormOne{T,OT<:MOI.ModelLike} = SingleBridgeOptimizer{NormOneBridge{T},OT}
include("geomean_to_relentr.jl")
const GeoMeantoRelEntr{T,OT<:MOI.ModelLike} =
    SingleBridgeOptimizer{GeoMeantoRelEntrBridge{T},OT}
include("geomean.jl")
const GeoMean{T,OT<:MOI.ModelLike} = SingleBridgeOptimizer{GeoMeanBridge{T},OT}
include("relentr_to_exp.jl")
const RelativeEntropy{T,OT<:MOI.ModelLike} =
    SingleBridgeOptimizer{RelativeEntropyBridge{T},OT}
include("norm_spec_nuc_to_psd.jl")
const NormSpectral{T,OT<:MOI.ModelLike} =
    SingleBridgeOptimizer{NormSpectralBridge{T},OT}
const NormNuclear{T,OT<:MOI.ModelLike} =
    SingleBridgeOptimizer{NormNuclearBridge{T},OT}
include("square.jl")
const Square{T,OT<:MOI.ModelLike} = SingleBridgeOptimizer{SquareBridge{T},OT}
include("det.jl")
const LogDet{T,OT<:MOI.ModelLike} = SingleBridgeOptimizer{LogDetBridge{T},OT}
const RootDet{T,OT<:MOI.ModelLike} = SingleBridgeOptimizer{RootDetBridge{T},OT}
include("soc_to_psd.jl")
const SOCtoPSD{T,OT<:MOI.ModelLike} =
    SingleBridgeOptimizer{SOCtoPSDBridge{T},OT}
const RSOCtoPSD{T,OT<:MOI.ModelLike} =
    SingleBridgeOptimizer{RSOCtoPSDBridge{T},OT}
include("indicator_activate_on_zero.jl")
include("indicator_sos.jl")
const IndicatortoSOS1{T,BC<:MOI.AbstractScalarSet,MaybeBC} =
    SingleBridgeOptimizer{IndicatorSOS1Bridge{T,BC,MaybeBC}}
include("semi_to_binary.jl")
const SemiToBinary{T,OT<:MOI.ModelLike} =
    SingleBridgeOptimizer{SemiToBinaryBridge{T},OT}
include("zero_one.jl")
const ZeroOne{T,OT<:MOI.ModelLike} = SingleBridgeOptimizer{ZeroOneBridge{T},OT}

"""
    add_all_bridges(bridged_model, ::Type{T})

Add all bridges defined in the `Bridges.Constraint` submodule to
`bridged_model`. The coefficient type used is `T`.
"""
function add_all_bridges(bridged_model, ::Type{T}) where {T}
    if T <: AbstractFloat
        MOIB.add_bridge(bridged_model, GreaterToIntervalBridge{T})
        MOIB.add_bridge(bridged_model, LessToIntervalBridge{T})
    end

    MOIB.add_bridge(bridged_model, GreaterToLessBridge{T})
    MOIB.add_bridge(bridged_model, LessToGreaterBridge{T})
    MOIB.add_bridge(bridged_model, NonnegToNonposBridge{T})
    MOIB.add_bridge(bridged_model, NonposToNonnegBridge{T})
    MOIB.add_bridge(bridged_model, ScalarizeBridge{T})
    MOIB.add_bridge(bridged_model, VectorizeBridge{T})
    MOIB.add_bridge(bridged_model, ScalarSlackBridge{T})
    MOIB.add_bridge(bridged_model, VectorSlackBridge{T})
    MOIB.add_bridge(bridged_model, ScalarFunctionizeBridge{T})
    MOIB.add_bridge(bridged_model, VectorFunctionizeBridge{T})
    MOIB.add_bridge(bridged_model, SplitIntervalBridge{T})
    MOIB.add_bridge(bridged_model, QuadtoSOCBridge{T})
    # We do not add `(R)SOCtoNonConvexQuad` because it starts with a convex
    # conic constraint and generate a non-convex constraint (in the QCP
    # interpretation).
    MOIB.add_bridge(bridged_model, NormInfinityBridge{T})
    MOIB.add_bridge(bridged_model, NormOneBridge{T})
    MOIB.add_bridge(bridged_model, GeoMeantoRelEntrBridge{T})
    MOIB.add_bridge(bridged_model, GeoMeanBridge{T})
    MOIB.add_bridge(bridged_model, RelativeEntropyBridge{T})
    MOIB.add_bridge(bridged_model, NormSpectralBridge{T})
    MOIB.add_bridge(bridged_model, NormNuclearBridge{T})
    MOIB.add_bridge(bridged_model, SquareBridge{T})
    MOIB.add_bridge(bridged_model, LogDetBridge{T})
    MOIB.add_bridge(bridged_model, RootDetBridge{T})
    MOIB.add_bridge(bridged_model, RSOCBridge{T})
    MOIB.add_bridge(bridged_model, SOCRBridge{T})
    # We do not add `SOCtoPSDBridge` as transforming the `SOC` to `RSOC` and
    # then to `PSD` produces a smaller SDP constraint.
    MOIB.add_bridge(bridged_model, RSOCtoPSDBridge{T})
    MOIB.add_bridge(bridged_model, IndicatorActiveOnFalseBridge{T})
    MOIB.add_bridge(bridged_model, IndicatorSOS1Bridge{T})
    MOIB.add_bridge(bridged_model, SemiToBinaryBridge{T})
    MOIB.add_bridge(bridged_model, ZeroOneBridge{T})
    return
end

end


text:
"""
    setulb(n, m, x, l, u, nbd, f, g, factr, pgtol, wa, iwa, task, iprint, csave, lsave, isave, dsave)
This function wraps L-BFGS-B's fortran subroutine `setulb` which partitions the working arrays `wa` and `iwa`,
and then uses the limited memory BFGS method to solve the bound constrained optimization problem by calling `mainlb`.

Note, `nmax` is the dimension of the largest problem to be solved, `mmax` is the maximum number of limited memory corrections.

# Arguments
- `n::Ref{Cint}`: the dimension of the problem.
- `m::Ref{Cint}`: the maximum number of variable metric corrections used to define the limited memory matrix.
- `x::Vector{Cdouble}`: an approximation to the solution, of length n.
- `l::Vector{Cdouble}`: the lower bound on x, of length n.
- `u::Vector{Cdouble}`: the upper bound on x, of length n.
- `nbd::Vector{Cint}`: represents the type of bounds imposed on the variables, and must be specified as follows:
    1. nbd[i]=0 if x[i] is unbounded;
    2. nbd[i]=1 if x(i) has only a lower bound;
    3. nbd[i]=2 if x[i] has both lower and upper bounds;
    4. nbd[i]=3 if x[i] has only an upper bound.
- `f::Ref{Cdouble}`: the value of the function at x(on first entry f is unspecified).
- `g::Vector{Cdouble}`: the value of the gradient at x(on first entry g is unspecified).
- `factr::Ref{Cdouble}`: the iteration will stop when (f^k - f^{k+1})/max{|f^k|,|f^{k+1}|,1} <= factr*epsmch,
    where epsmch is the machine precision, which is automatically generated by the code. Typical values for factr:
    * 1.e12 for low accuracy;
    * 1.e7 for moderate accuracy;
    * 1.e1 for extremely high accuracy.
- `pgtol::Ref{Cdouble}`: the iteration will stop when max{|proj g_i | i = 1, ..., n} <= pgtol where pg_i is the ith component of the projected gradient.
- `wa::Vector{Cdouble}`: the working array of length (2mmax + 5)nmax + 12mmax^2 + 12mmax.
- `iwa::Vector{Cint}`: the integer working array of length 3nmax.
- `task::Vector{Cuchar}`: a working string of characters of length 60 indicating the current job when entering and quitting this function.
- `iprint::Ref{Cint}`: it controls the frequency and type of output generated:
    1. iprint<0 no output is generated;
    2. iprint=0 print only one line at the last iteration;
    3. 0<iprint<99 print also f and |proj g| every iprint iterations;
    4. iprint=99 print details of every iteration except n-vectors;
    5. iprint=100 print also the changes of active set and final x;
    6. iprint>100 print details of every iteration including x and g.
    7. when iprint > 0, the file iterate.dat will be created to summarize the iteration.
- `csave::Vector{Cuchar}`: a working string of characters of length 60.
- `lsave::Vector{Bool}`: a logical working array of dimension 4. if `task` = NEW_X, the following information is available:
    1. if lsave[1] == true then the initial X has been replaced by its projection in the feasible set;
    2. if lsave[2] == true then the problem is constrained;
    3. if lsave[3] == true then each variable has upper and lower bounds.
- `isave::Vector{Cint}`: an integer working array of dimension 44. if `task` = NEW_X, the following information is available:
    1. isave[22] = the total number of intervals explored in the search of Cauchy points;
    2. isave[26] = the total number of skipped BFGS updates before the current iteration;
    3. isave[30] = the number of current iteration;
    4. isave[31] = the total number of BFGS updates prior the current iteration;
    5. isave[33] = the number of intervals explored in the search of Cauchy point in the current iteration;
    6. isave[34] = the total number of function and gradient evaluations;
    7. isave[36] = the number of function value or gradient evaluations in the current iteration;
    8. if isave[37] = 0  then the subspace argmin is within the box;
    9. if isave[37] = 1  then the subspace argmin is beyond the box;
    10. isave[38] = the number of free variables in the current iteration;
    11. isave[39] = the number of active constraints in the current iteration;
    12. n + 1 - isave[40] = the number of variables leaving the set of active constraints in the current iteration;
    13. isave[41] = the number of variables entering the set of active constraints in the current iteration.
- `dsave::Vector{Cdouble}`: a double precision working array of dimension 29. if `task` = NEW_X, the following information is available:
    1. dsave[1] = current 'theta' in the BFGS matrix;
    2. dsave[2] = f(x) in the previous iteration;
    3. dsave[3] = factr*epsmch;
    4. dsave[4] = 2-norm of the line search direction vector;
    5. dsave[5] = the machine precision epsmch generated by the code;
    6. dsave[7] = the accumulated time spent on searching for Cauchy points;
    7. dsave[8] = the accumulated time spent on subspace minimization;
    8. dsave[9] = the accumulated time spent on line search;
    9. dsave[11] = the slope of the line search function at the current point of line search;
    10. dsave[12] = the maximum relative step length imposed in line search;
    11. dsave[13] = the infinity norm of the projected gradient;
    12. dsave[14] = the relative step length in the line search;
    13. dsave[15] = the slope of the line search function at the starting point of the line search;
    14. dsave[16] = the square of the 2-norm of the line search direction vector.
"""
function setulb(n, m, x, l, u, nbd, f, g, factr, pgtol, wa, iwa, task, iprint, csave, lsave, isave, dsave)
    global liblbfgsb
    hdl = Libdl.dlopen_e(liblbfgsb)
    @assert hdl != C_NULL "Could not open $liblbfgsb"
    setulb_ = Libdl.dlsym_e(hdl, :setulb_)
    @assert setulb_ != C_NULL "Could not find `setulb` within $liblbfgsb"
    ccall(setulb_, Cvoid, (Ref{Cint}, Ref{Cint}, Ref{Cdouble}, Ref{Cdouble}, Ref{Cdouble},
          Ref{Cint}, Ref{Cdouble}, Ref{Cdouble}, Ref{Cdouble}, Ref{Cdouble}, Ref{Cdouble},
          Ref{Cint}, Ref{Cuchar}, Ref{Cint}, Ref{Cuchar}, Ref{Cint}, Ref{Cint}, Ref{Cdouble}, Csize_t, Csize_t),
          n, m, x, l, u, nbd, f, g, factr, pgtol, wa, iwa, task, iprint, csave, lsave, isave, dsave, 60, 60)
end

"""
    timer(x)
The double precision cpu timing subroutine in the L-BFGS-B code.
"""
function timer(x)
    global liblbfgsb
    hdl = Libdl.dlopen_e(liblbfgsb)
    @assert hdl != C_NULL "Could not open $liblbfgsb"
    timer_ = Libdl.dlsym_e(hdl, :timer_)
    @assert timer_ != C_NULL "Could not find `timer` within $liblbfgsb"
    ccall(timer_, Cvoid, (Ref{Cdouble},), x)
end


text:
"""
$(SIGNATURES)

CarrierReplicaCPU that holds possible carrier representations.
The type is unknown, when TrackingState is initialized.
It is determined based on the signal, which is given in the track
function.
"""
struct CarrierReplicaCPU{
        CF32 <: StructArray{ComplexF32},
        CF64 <: StructArray{ComplexF64}
    }
    carrier_f32::CF32
    carrier_f64::CF64
end

function CarrierReplicaCPU()
    CarrierReplicaCPU(
        StructArray{Complex{Float32}}(undef, 0),
        StructArray{Complex{Float64}}(undef, 0)
    )
end
"""
$(SIGNATURES)

DownconvertedSignalCPU that holds possible downconverted signal representations.
The type is unknown, when TrackingState is initialized.
It is determined based on the signal, which is given in the track
function.
"""
struct DownconvertedSignalCPU{
        DSF32 <: StructArray{ComplexF32},
        DSF64 <: StructArray{ComplexF64}
    }
    downconverted_signal_f32::DSF32
    downconverted_signal_f64::DSF64
end

function DownconvertedSignalCPU(num_ants::NumAnts{1})
    DownconvertedSignalCPU(
        StructArray{Complex{Float32}}(undef, 0),
        StructArray{Complex{Float64}}(undef, 0)
    )
end

function DownconvertedSignalCPU(num_ants::NumAnts{N}) where N
    DownconvertedSignalCPU(
        StructArray{Complex{Float32}}(undef, 0, N),
        StructArray{Complex{Float64}}(undef, 0, N)
    )
end

"""
$(SIGNATURES)

TrackingState that holds the tracking state.
"""
struct TrackingState{
        S <: AbstractGNSS,
        C <: AbstractCorrelator,
        CALF <: AbstractLoopFilter,
        COLF <: AbstractLoopFilter,
        CN <: AbstractCN0Estimator,
        DS <: Union{DownconvertedSignalCPU, Nothing},
        CAR <: Union{CarrierReplicaCPU, Nothing},
        COR <: Union{Vector{Int8}, Nothing}
    }
    prn::Int
    system::S
    init_carrier_doppler::typeof(1.0Hz)
    init_code_doppler::typeof(1.0Hz)
    carrier_doppler::typeof(1.0Hz)
    code_doppler::typeof(1.0Hz)
    carrier_phase::Float64
    code_phase::Float64
    correlator::C
    carrier_loop_filter::CALF
    code_loop_filter::COLF
    sc_bit_detector::SecondaryCodeOrBitDetector
    integrated_samples::Int
    prompt_accumulator::ComplexF64
    cn0_estimator::CN
    downconverted_signal::DS
    carrier::CAR
    code::COR
end

"""
$(SIGNATURES)

Convenient TrackingState constructor. Mandatory parameters are the GNSS system, the
carrier doppler `carrier_doppler` and the code phase `code_phase`. Optional parameters are
- Code doppler `code_doppler`, that defaults to carrier doppler multiplied with the code /
  center frequency ratio
- Carrier phase `carrier_phase`, that defaults to 0
- Carrier loop filter `carrier_loop_filter`, that defaults to the `ThirdOrderBilinarLF()`
- Code loop filter `code_loop_filter`, that defaults to the `SecondOrderBilinearLF()`
- Number of antenna elemants `num_ants`, that defaults to `NumAnts(1)`
- Correlator `correlator`, that defaults to `get_default_correlator(...)`
- Integrated samples `integrated_samples`, that defaults to 0
- Prompt accumulator for bit detection `prompt_accumulator`, that defaults to 0
- CN0 estimator `cn0_estimator`, that defaults to `MomentsCN0Estimator(20)`
"""
function TrackingState(
    prn::Integer,
    system::S,
    carrier_doppler,
    code_phase;
    code_doppler = carrier_doppler * get_code_center_frequency_ratio(system),
    carrier_phase = 0.0,
    carrier_loop_filter::CALF = ThirdOrderBilinearLF(),
    code_loop_filter::COLF = SecondOrderBilinearLF(),
    sc_bit_detector = SecondaryCodeOrBitDetector(),
    num_ants = NumAnts(1),
    correlator::C = get_default_correlator(system, num_ants),
    integrated_samples = 0,
    prompt_accumulator = zero(ComplexF64),
    cn0_estimator::CN = MomentsCN0Estimator(20),
    num_samples = 0
) where {
    T <: AbstractMatrix,
    S <: AbstractGNSS{T},
    C <: AbstractCorrelator,
    CALF <: AbstractLoopFilter,
    COLF <: AbstractLoopFilter,
    CN <: AbstractCN0Estimator
}
    if found(sc_bit_detector)
        code_phase = mod(code_phase, get_code_length(system) *
            get_secondary_code_length(system))
    else
        code_phase = mod(code_phase, get_code_length(system))
    end
    downconverted_signal = DownconvertedSignalCPU(num_ants)
    carrier = CarrierReplicaCPU()
    code = Vector{Int8}(undef, 0)

    TrackingState{S, C, CALF, COLF, CN, typeof(downconverted_signal), typeof(carrier), typeof(code)}(
        prn,
        system,
        carrier_doppler,
        code_doppler,
        carrier_doppler,
        code_doppler,
        carrier_phase / 2π,
        code_phase,
        correlator,
        carrier_loop_filter,
        code_loop_filter,
        sc_bit_detector,
        integrated_samples,
        prompt_accumulator,
        cn0_estimator,
        downconverted_signal,
        carrier,
        code
    )
end

# CUDA dispatch
function TrackingState(
    prn::Integer,
    system::S,
    carrier_doppler,
    code_phase;
    num_samples,
    code_doppler = carrier_doppler * get_code_center_frequency_ratio(system),
    carrier_phase = 0.0,
    carrier_loop_filter::CALF = ThirdOrderBilinearLF(),
    code_loop_filter::COLF = SecondOrderBilinearLF(),
    sc_bit_detector = SecondaryCodeOrBitDetector(),
    num_ants = NumAnts(1),
    correlator::C = get_default_correlator(system, num_ants),
    integrated_samples = 0,
    prompt_accumulator = zero(ComplexF64),
    cn0_estimator::CN = MomentsCN0Estimator(20)
) where {
    T <: CuMatrix,
    S <: AbstractGNSS{T},
    C <: AbstractCorrelator,
    CALF <: AbstractLoopFilter,
    COLF <: AbstractLoopFilter,
    CN <: AbstractCN0Estimator
}
    if found(sc_bit_detector)
        code_phase = mod(code_phase, get_code_length(system) *
            get_secondary_code_length(system))
    else
        code_phase = mod(code_phase, get_code_length(system))
    end
    downconverted_signal = nothing
    carrier = nothing
    code = nothing

    TrackingState{S, C, CALF, COLF, CN, Nothing, Nothing, Nothing}(
        prn,
        system,
        carrier_doppler,
        code_doppler,
        carrier_doppler,
        code_doppler,
        carrier_phase / 2π,
        code_phase,
        correlator,
        carrier_loop_filter,
        code_loop_filter,
        sc_bit_detector,
        integrated_samples,
        prompt_accumulator,
        cn0_estimator,
        downconverted_signal,
        carrier,
        code
    )
end

@inline get_system(state::TrackingState) = state.system
@inline get_code_phase(state::TrackingState) = state.code_phase
@inline get_carrier_phase(state::TrackingState) = state.carrier_phase
@inline get_init_code_doppler(state::TrackingState) = state.init_code_doppler
@inline get_init_carrier_doppler(state::TrackingState) = state.init_carrier_doppler
@inline get_code_doppler(state::TrackingState) = state.code_doppler
@inline get_carrier_doppler(state::TrackingState) = state.carrier_doppler
@inline get_correlator(state::TrackingState) = state.correlator
@inline get_sc_bit_detector(state::TrackingState) = state.sc_bit_detector
@inline get_carrier_loop_filter(state::TrackingState) = state.carrier_loop_filter
@inline get_code_loop_filter(state::TrackingState) = state.code_loop_filter
@inline get_prompt_accumulator(state::TrackingState) = state.prompt_accumulator
@inline get_integrated_samples(state::TrackingState) = state.integrated_samples
@inline get_cn0_estimator(state::TrackingState) = state.cn0_estimator
@inline get_downconverted_signal(state::TrackingState) = state.downconverted_signal
@inline get_carrier(state::TrackingState) = state.carrier
@inline get_code(state::TrackingState) = state.code
@inline get_prn(state::TrackingState) = state.prn


text:
# MIT License
#
# Copyright (c) 2018 Martin Biel
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

@with_kw mutable struct DryFrictionParameters{T}
    γ::T = 0.9
    r::T = 1e-3
end

"""
    DryFrictionProximal

Functor object for using dry-friction acceleration in the prox step of a quasigradient algorithm. Create by supplying a [`DryFriction`](@ref) object through `prox ` to `QuasiGradient.Optimizer` or by setting the [`Prox`](@ref) attribute.

...
# Parameters
- `prox::AbstractProx = Polyhedron`: Inner prox step
- `γ::AbstractFloat = 0.9`: Heavy-ball parameter
- `r::AbstractFloat = 1e-3`: Dry-friction parameter
...
"""
struct DryFrictionProximal{T <: AbstractFloat, P <: AbstractProximal} <: AbstractProximal
    parameters::DryFrictionParameters{T}

    prox::P
    xprev::Vector{T}

    function DryFrictionProximal(proximal::AbstractProximal, ::Type{T}; kw...) where T <: AbstractFloat
        P = typeof(proximal)
        return new{T,P}(DryFrictionParameters{T}(; kw...),
                        proximal,
                        Vector{T}())
    end
end

function initialize_prox!(quasigradient::AbstractQuasiGradient, dryfriction::DryFrictionProximal)
    # Initialize inner
    initialize_prox!(quasigradient, dryfriction.prox)
    return nothing
end

function restore_proximal_master!(quasigradient::AbstractQuasiGradient, dryfriction::DryFrictionProximal)
    # Restore inner
    restore_proximal_master!(quasigradient, dryfriction.prox)
    return nothing
end

function prox!(quasigradient::AbstractQuasiGradient, dryfriction::DryFrictionProximal, x::AbstractVector, ∇f::AbstractVector, h::AbstractFloat)
    if length(dryfriction.xprev) == 0
        # First iteration standard prox
        resize!(dryfriction.xprev, length(x))
        dryfriction.xprev .= x
        prox!(quasigradient, dryfriction.prox, x, ∇f, h)
        return nothing
    end
    @unpack γ, r = dryfriction.parameters
    λ = h / (1 + h * γ)
    z = (1 / (1 + h * γ)) * (x - dryfriction.xprev) / h - λ * ∇f
    proj_z = (1 - λ*r / max(λ*r, norm(z))) * z
    # Check for restart
    if norm(proj_z) <= sqrt(eps())
        dryfriction.parameters.r = 0.1 * r
        proj_z .= -∇f
    end
    dryfriction.xprev .= x
    # Regular prox step
    prox!(quasigradient, dryfriction.prox, x, -proj_z, h)
    return nothing
end

# API
# ------------------------------------------------------------
"""
    DryFriction

Factory object for [`DryFrictionProximal`](@ref). Pass to `prox` in `Quasigradient.Optimizer` or set the [`Prox`](@ref) attribute. See ?DryFrictionProximal for parameter descriptions.

"""
mutable struct DryFriction <: AbstractProx
    prox::AbstractProx
    parameters::DryFrictionParameters
end
DryFriction(; prox::AbstractProx = Polyhedron(), kw...) = DryFriction(prox, DryFrictionParameters(; kw...))

function (dryfriction::DryFriction)(structure::StageDecompositionStructure, x₀::AbstractVector, ::Type{T}) where T <: AbstractFloat
    proximal = dryfriction.prox(structure, x₀, T)
    return DryFrictionProximal(proximal, T; type2dict(dryfriction.parameters)...)
end

function str(::DryFriction)
    return ""
end


text:
"""
STOCHASTIC GRADIENT HAMILTONIAN MONTE CARLO

The following codes explore the use of SGHMC for Bayesian Inference
"""
immutable SGHMC
  dU      ::Function
  dK      ::Function
  dKΣ     ::Array{Float64}
  C       ::Array{Float64}
  V       ::Array{Float64}
  init_est::Array{Float64}
  d       ::Int64
end

function mcmc(parameters::SGHMC;
  leapfrog_params::Dict{Symbol, Real} = Dict([:ɛ => .05, :τ => 20]),
  set_seed::Int64 = 123,
  r::Int64 = 1000)

  dU, dK, dKΣ, C, V, w, d = parameters.dU, parameters.dK, parameters.dKΣ, parameters.C, parameters.V, parameters.init_est, parameters.d
  ɛ, τ = leapfrog_params[:ɛ], leapfrog_params[:τ]

  if typeof(set_seed) == Int64
    srand(set_seed)
  end

  chain = zeros(r, d);
  chain[1, :] = w
  B = .5 * V * ɛ
  D = sqrt(2 * (C - B) * ɛ)
  if size(B) != size(C)
    error("C and V should have the same dimension.")
  else
    if sum(size(B)) > 1
      if det(B) > det(C)
        error("ɛ is too big. Consider decreasing it.")
      end
    else
      if det(B[1]) > det(C[1])
        error("ɛ is too big. Consider decreasing it.")
      end
    end
  end

  for i in 1:(r - 1)
    p = randn(d, 1)

    for j in 1:τ
      p = p - dU(w) * ɛ - C * inv(dKΣ) * p + D * randn(d, 1);
      w = w + dK(p) * ɛ;
    end

    chain[i + 1, :] = w
  end

  return chain
end


text:
@testset "polyagamma" begin
    @test mean(PolyaGamma(1, 0)) == 1 / 4
    @test mean(PolyaGamma(1, 2.0)) == tanh(1.0) / 4

    for (b, c) in ((1, 0), (1, 2.0), (3, 0), (3, 2.5), (3, 3.2), (1.2, 3.2))
        p = PolyaGamma(b, c)
        @test logpdf(p, rand(p)) isa Real
        @test mean(rand(p, 10000)) ≈ mean(p) atol = 1e-2
    end

    p = PolyaGamma(1, 0)
    @test Distributions.insupport(p, 0)
    @test !Distributions.insupport(p, -1)
    @test minimum(p) === 0
    @test maximum(p) == Inf
    @test Distributions.params(p) == (1, 0)
end


text:
using GalacticGCMAES, GalacticOptim, ForwardDiff
using Test

@testset "GalacticGCMAES.jl" begin
    rosenbrock(x, p) = (p[1] - x[1])^2 + p[2] * (x[2] - x[1]^2)^2
    x0 = zeros(2)
    _p = [1.0, 100.0]
    l1 = rosenbrock(x0, _p)
    f_ad = OptimizationFunction(rosenbrock, GalacticOptim.AutoForwardDiff())
    f_noad = OptimizationFunction(rosenbrock)

    prob = GalacticOptim.OptimizationProblem(f_ad, x0, _p, lb=[-1.0, -1.0], ub=[1.0, 1.0])
    sol = solve(prob, GCMAESOpt(), maxiters=1000)
    @test 10 * sol.minimum < l1

    prob = GalacticOptim.OptimizationProblem(f_noad, x0, _p, lb=[-1.0, -1.0], ub=[1.0, 1.0])
    sol = solve(prob, GCMAESOpt(), maxiters=1000)
    @test 10 * sol.minimum < l1
end


text:
using Revise
using SnpArrays
using DelimitedFiles
using Random
using MendelIHT
using Knockoffs
using LinearAlgebra
using GLMNet
using Distributions
using DataFrames
using CSV
using Printf

fdr = 0.6
seed = 12
cd("/scratch/users/bbchu/ukb/prs/fdr$fdr/sim$seed")

function R2(X::AbstractMatrix, y::AbstractVector, β̂::AbstractVector)
    μ = y - X * β̂
    tss = y .- mean(y)
    return 1 - dot(μ, μ) / dot(tss, tss)
end

chr = 10
plinkname = "/scratch/users/bbchu/ukb/subset/ukb.10k.chr10"
knockoffname = "/scratch/users/bbchu/ukb/subset/ukb.10k.merged.chr10"
x = SnpArray(plinkname * ".bed")
xko = SnpArray(knockoffname * ".bed")
# xla = SnpLinAlg{Float64}(x, center=true, scale=true, impute=true)
xla = convert(Matrix{Float64}, x, center=true, scale=true, impute=true)
xko_la = convert(Matrix{Float64}, xko, center=true, scale=true, impute=true)
original = vec(readdlm("/scratch/users/bbchu/ukb/subset/ukb.chr$chr.original.snp.index", Int))
knockoff = vec(readdlm("/scratch/users/bbchu/ukb/subset/ukb.chr$chr.knockoff.snp.index", Int))

#
# simulate phenotypes using UKB chr10 subset
#
n, p = size(x)
# simulate β
Random.seed!(seed)
k = 100 # number of causal SNPs
h2 = 0.5 # heritability
d = Normal(0, sqrt(h2 / (2k))) # from paper: Efficient Implementation of Penalized Regression for Genetic Risk Prediction
β = zeros(p)
β[1:k] .= rand(d, k)
shuffle!(β)
# simulate y
ϵ = Normal(0, 1 - h2)
y = xla * β + rand(ϵ, n)

#
# compare R2 across populations, save result in a dataframe
#
β_iht = vec(readdlm("iht.beta"))
β_iht_knockoff = extract_beta(vec(readdlm("iht.knockoff.beta")), fdr, original, knockoff)
β_iht_knockoff_cv = extract_beta(vec(readdlm("iht.knockoff.cv.beta")), fdr, original, knockoff)
β_iht_knockoff_combined_cv = extract_combine_beta(vec(readdlm("iht.knockoff.cv.beta")), original, knockoff)
β_lasso = vec(readdlm("lasso.beta"))
β_lasso_knockoff = extract_beta(vec(readdlm("lasso.knockoff.beta")), fdr, original, knockoff)
β_lasso_knockoff_combined = extract_combine_beta(vec(readdlm("lasso.knockoff.beta")), original, knockoff)

pop = "african"

xpop = SnpArray("/scratch/users/bbchu/ukb/populations/chr10/ukb.chr$chr.$pop.bed")
Xpop = SnpLinAlg{Float64}(xpop, center=true, scale=true, impute=true)
# simulate "true" phenotypes for these populations
ytrue = Xpop * β + rand(ϵ, size(Xpop, 1))

# IHT
iht_r2 = R2(Xpop, ytrue, β_iht)
iht_ko_r2 = R2(Xpop, ytrue, β_iht_knockoff) # knockoff IHT
iht_ko_cv_r2 = R2(Xpop, ytrue, β_iht_knockoff_cv) # better cv routine
iht_ko_cv_combined_r2 = R2(Xpop, ytrue, β_iht_knockoff_combined_cv) # original β combined with its knockoff

# LASSO
lasso_r2 = R2(Xpop, ytrue, β_lasso)
lasso_ko_r2 = R2(Xpop, ytrue, β_lasso_knockoff) # knockoff IHT
lasso_ko_combined_r2 = R2(Xpop, ytrue, β_lasso_knockoff_combined) # original β combined with its knockoff


text:
using Juno
using Profile
using BenchmarkTools
import TrajectoryOptimization: reset_evals, dynamics_jacobian, outer_loop_update

# Set up
N = 101
model, obj = Dynamics.dubinscar
model, obj = Dynamics.quadrotor
n,m = model.n, model.m

x0 = [0.0;0.0;0.]
xf = [0.0;1.0;0.]
tf =  3.
Qf = 100.0*Diagonal(I,n)
Q = (1e-3)*Diagonal(I,n)
R = (1e-2)*Diagonal(I,m)

obj = LQRObjective(Q, R, Qf, tf, x0, xf)

x_min = [-0.25; -0.001; -Inf]
x_max = [0.25; 1.001; Inf]

obj_con = ConstrainedObjective(obj,x_min=x_min,x_max=x_max)

solver = Solver(model,obj,N=N)
solver.opts.verbose = true
solver.opts.square_root = true
solver.state.infeasible
n,m = get_sizes(solver)

solver2 = Solver(model,obj_con,N=N)
solver2.opts.al_type = :algencan

X0 = Array{Float64,2}(undef,0,0)
U0 = ones(m,N)

logger = default_logger(solver)
global_logger(logger)

#****************************#
#       INITIALIZATION       #
#****************************#
mbar,mm = get_num_controls(solver)
get_num_constraints(solver)

results = init_results(solver, X0, U0)
results2 = init_results(solver2, X0, U0)
results.ρ[1] = 0
results2.ρ[1] = 0
X,U = results.X, results.U
X_,U_ = results.X_, results.U_

#****************************#
#           SOLVER           #
#****************************#
## Initial rollout
if !solver.state.infeasible
    X[1] = solver.obj.x0
    flag = rollout!(results,solver) # rollout new state trajectoy
    # rollout!(results2,solver2)

    if !flag
        @info "Bad initial control sequence, setting initial control to zero"
        results.U .= zeros(mm,N)
        rollout!(results,solver)
    end
end

if solver.state.infeasible
    if solver.control_integration == :foh
        calculate_derivatives!(results, solver, results.X, results.U)
        calculate_midpoints!(results, solver, results.X, results.U)
    end
    update_constraints!(results,solver,results.X,results.U)
    update_constraints!(results_new,solver_new,results.X,results.U)
end
# p = plot()
# plot_trajectory!(results)
# plot(to_array(results.X)')

# results.X == results_new.X.x
# results.μ == results_new.μ.x

reset_evals(solver)
Fc = dynamics_jacobian(solver)
Fd = discrete_dynamics_jacobian(solver)
x = rand(n)
u = rand(m)
xdot = zero(x)
A,B = zeros(n,n),zeros(n,m)
solver.fd(xdot,x,u)
Fc(x,u)
Fd(A,B,x,u)
evals(solver,:f)

results.S
J_prev = cost(solver, results)
TrajectoryOptimization.update_jacobians!(results, solver)
reset_evals(solver)
results.ρ[1]
Δv = _backwardpass_sqrt!(results, solver)
S1 = copy(results.S[1])
Sfull = S1'S1
S2 = qr(S1).R
S2'S2 ≈ Sfull

A = rand(4,4)
B = rand(6,4)
P = zeros(10,4)
@btime chol_plus($A,$B)
@btime chol_plus!(P,A,B)

K1 = copy(results.K[1])
_backwardpass_sqrt_noqr!(results, solver)

@btime _backwardpass_sqrt_noqr!(results, solver)
@btime _backwardpass_sqrt!(results, solver)

J = forwardpass!(results, solver, Δv, J_prev)
c_max = max_violation(results)
results.X .= deepcopy(results.X_)
results.U .= deepcopy(results.U_)

print_header(logger, InnerLoop)
J_prev = cost(solver2, results2)
update_jacobians!(results2, solver2)
Δv = _backwardpass_new!(results2, solver2)
J = forwardpass!(results2, solver2, Δv, J_prev)
c_max = max_violation(results2)
print_row(logger, InnerLoop)

copyto!(results2.X,results2.X_)
copyto!(results2.U,results2.U_)

outer_loop_update(results2, solver2)

#
# b2 = @benchmark copyto!(results_new.X,results_new.X_)
# b1 = @benchmark results.X .= deepcopy(results.X_)
# judge(median(b2),median(b1))

outer_loop_update(results,solver)
cost(solver,results)
max_violation(results)
update_constraints!(results, solver)

outer_loop_update(results_new,solver_new)
update_constraints!(results_new, solver_new)

max_violation(results)
max_violation(results_new)

b1 = @benchmark outer_loop_update($results,$solver)
b2 = @benchmark outer_loop_update($results_new,$solver_new)
judge(median(b2),median(b1))

b1 = @benchmark update_constraints!($results,$solver)
b2 = @benchmark update_constraints!($results_new,$solver_new)
judge(median(b2),median(b1))

b1 = @benchmark max_violation($results)
b2 = @benchmark max_violation($results_new)
judge(median(b2),median(b1))

# Warm start
λ = deepcopy(results.λ)
push!(λ, results.λN)
U0_warm = to_array(results.U)
res_warm = init_results(solver, X0, U0_warm, λ=λ)
rollout!(res_warm,solver)

J_prev = cost(solver, res_warm)
update_jacobians!(res_warm, solver)
Δv = backwardpass!(res_warm, solver)
J = forwardpass!(res_warm, solver, Δv)


text:
module PlutoStaticHTML

if isdefined(Base, :Experimental) && isdefined(Base.Experimental, Symbol("@max_methods"))
    @eval Base.Experimental.@max_methods 1
end

import Base:
    show,
    string
import Pluto:
    PlutoRunner,
    WorkspaceManager

using Base64: base64encode
using Dates
using Pkg:
    Types.Context,
    Types.UUID,
    Operations
using Pluto:
    BondValue,
    Cell,
    CellOutput,
    Configuration.CompilerOptions,
    Notebook,
    PkgCompat.dependencies,
    Pluto,
    PlutoRunner,
    ServerSession,
    SessionActions,
    WorkspaceManager,
    generate_html,
    load_notebook_nobackup,
    update_dependency_cache!,
    update_run!,
    update_save_run!
using PrecompileSignatures: @precompile_signatures
using SHA: sha256
using TOML: parse as parsetoml

include("module_doc.jl")
include("context.jl")
include("cache.jl")
include("mimeoverride.jl")
include("with_terminal.jl")
include("html.jl")
include("style.jl")
include("build.jl")
include("documenter.jl")

export HTMLOptions
export documenter_output, franklin_output, html_output
export BuildOptions, build_notebooks

@precompile_signatures(PlutoStaticHTML)

end # module


text:
struct K80rel <: K80
  κ::Float64
  function K80rel(κ::Float64;
                  safe::Bool=true)
    if safe
      if κ <= 0.
        error("K80 parameter κ must be positive")
      end
    end
    new(κ)
  end
end

function Base.show(io::IO, object::K80rel)
  print(io, "\r\e[0m\e[1mK\e[0mimura 19\e[1m80\e[0m model (relative rate form)
κ = $(object.κ)")
end

@inline function Q(mod::K80rel)
  κ = mod.κ

  Q₁ = 0.25 * κ
  Q₂ = 0.25
  Q₃ = -(2 * Q₂ + Q₁)

  return Qmatrix(Q₃, Q₂, Q₁, Q₂,
                 Q₂, Q₃, Q₂, Q₁,
                 Q₁, Q₂, Q₃, Q₂,
                 Q₂, Q₁, Q₂, Q₃)
end

@inline function P(mod::K80rel, t::Float64)
  if t < 0.0
    error("t must be positive")
  end
  κ = mod.κ

  e₁ = 0.25 * exp(-t)
  e₂ = 0.5 * exp(-0.5 * t * (κ + 1))

  P₁ = 0.25 + e₁ + e₂
  P₂ = 0.25 + e₁ - e₂
  P₃ = 0.25 - e₁

  return Pmatrix(P₁, P₃, P₂, P₃,
                 P₃, P₁, P₃, P₂,
                 P₂, P₃, P₁, P₃,
                 P₃, P₂, P₃, P₁)
end


text:
include("activation_functions.jl")

"""
Make a linear forward calc.
"""
function linear_forward(A, W, b)
    # Make a linear forward and return the inputs as cache
    Z = BigFloat.((W * A) .+ b)
    cache = (A, W, b)

    @assert size(Z) == (size(W, 1), size(A, 2))

    return (Z = Z, cache = cache)
end # function

"""
Make a forward activation
"""
function linear_forward_activate(A_prev, W, b; activation_function = "tanh")
    @assert activation_function ∈ ("sigmoid", "relu", "softmax", "tanh", "swish")
    Z, linear_cache = linear_forward(A_prev, W, b)

    if activation_function == "sigmoid"
        A, activation_cache = sigmoid(Z)
    elseif activation_function == "relu"
        A, activation_cache = relu(Z)
    elseif activation_function == "softmax"
        A, activation_cache = softmax(Z)
    elseif activation_function == "tanh"
        A, activation_cache = tanhact(Z)
    elseif activation_function == "swish"
        A, activation_cache = swish(Z)
    end # if

    cache = (linear_step_cache = linear_cache, activation_step_cache = activation_cache)

    @assert size(A)[1] == size(W, 1)[1]

    return A, cache
end # function

max_stars_repo_name: Squalm/NeuralNetworks

max_issues_repo_name: Squalm/NeuralNetworks

max_forks_repo_name: Squalm/NeuralNetworks


text:
include("r3frac.jl")
include("r3mv.jl")
include("LatticeVars.jl")
include("AtomicVars.jl")
include("AtomicSpeciesVars.jl")
include("MuffinTinRadialVars.jl")
include("APWLOVars.jl")
include("readspecies.jl")
include("mtdmin.jl")
include("checkmt.jl")
include("genrmesh.jl")
include("radnucl.jl")
include("potnucl.jl")

include("XC_funcs/XC_x_slater.jl")
include("XC_funcs/XC_c_vwn.jl")

include("polynm.jl")
include("splint.jl")
include("splintwp.jl")
include("wsplint.jl")
include("wsplintp.jl")

include("solve_atom.jl")
include("rdirac.jl")
include("rdiracint.jl")


text:
# Na razie tutaj krotki "doc"

# Albo wywolujemy z tablicami liczb, wtedy trzeba podac N
# Np.
# K=[3 4 5]
# M=[4 3 2 4 5]
# C = Spritz_Encrypt(8,K,M);
# Spritz_Decrypt(8,K,C)
# Spritz_Hash(8,[3 4 5 4 2],32)
#
# Albo z "napisami", wtedy przyjmowane jest N=8, kazdy znak 'x' zamieniany jest na int('x'), np.
# Spritz_Hash("arcfour",32)
#
#

type Spritz
    i::Int64
    j::Int64
    k::Int64
    z::Int64
    a::Int64
    w::Int64
    N::Int64
    S

   #Spritz(n::Int64)=new(0,0,0,0,0,0,n,0)
   function Spritz(n::Int64)

    S=[]
     for v=1:n
	  S=[S,v-1]
      end

    return new(0,0,0,0,0,1,n,S)
   end

end

function Spritz_Show(sp::Spritz)
  println("i=",sp.i," j=",sp.j, " k=",sp.k, " z=",sp.z, " a=",sp.a, " w=",sp.w, " N=",sp.N)
  println("S=",sp.S )
end

function Spritz_Swap(sp::Spritz, i::Int64, j::Int64)
# println("Spritz_Swap(sp,",i,",",j,")")
  tmp = sp.S[i+1]
  sp.S[i+1]=sp.S[j+1]
  sp.S[j+1]=tmp
end

function Spritz_Update(sp::Spritz)
  sp.i=(sp.i+sp.w) % sp.N
  #sp.j=(sp.k+sp.S[sp.j+1 + sp.S[sp.i+1 % sp.N +1 ] % sp.N +1 ] ) % sp.N
 # println("i=",sp.i)
 # println("sp.j+ sp.S[sp.i] +1= ", sp.S[(sp.j+ sp.S[(sp.i%sp.N)+1])%sp.N+1]  )

 # sp.j= (sp.k+sp.S[(sp.j+ sp.S[((sp.i+1)%sp.N)+1])%sp.N+1])%sp.N
   sp.j= (sp.k+sp.S[ (sp.j+ sp.S[sp.i+1]) %sp.N+1])%sp.N

  sp.k=(sp.i+sp.k+sp.S[sp.j+1])%sp.N
  #sp.k=(sp.i + sp.k + sp.S[sp.j+1 % sp.N + 1] ) % sp.N

  Spritz_Swap(sp,sp.i,sp.j)
end

function Spritz_Whip(sp::Spritz, r::Int64)
  for v in 0:(r-1)
    Spritz_Update(sp)
  end
  #tutaj powinno byc: do w=w+1 until GCD(w,N)=1
  sp.w=sp.w+1
  while gcd(sp.w,sp.N)>1
     sp.w=sp.w+1
  end
end

function Spritz_Crush(sp::Spritz)
#println("convert(Int64,(floor(sp.N/2)-1)) = ", convert(Int64,(floor(sp.N/2)-1)))

for v in 0:convert(Int64,(floor(sp.N/2)-1))
    #println("v=",v)
    if sp.S[v+1]>sp.S[sp.N-1-v+1]
      #println("zamiana ",v+1, " z " , sp.N-1-v+1, ", czyli  Spritz_Swap(sp,",v,",",sp.N-1-v," <--)" )
      #Spritz_Show(sp)
    Spritz_Swap(sp,v,sp.N-1-v) end
end
end

function Spritz_Shuffle(sp::Spritz)
  Spritz_Whip(sp,2*sp.N)
  Spritz_Crush(sp)
  Spritz_Whip(sp,2*sp.N)
  Spritz_Crush(sp)
  Spritz_Whip(sp,2*sp.N)
  sp.a=0
end

function Spritz_AbsorbStop(sp::Spritz)
  if sp.a==convert(Int64,floor(sp.N/2)) Spritz_Shuffle(sp) end
  sp.a = sp.a+1
end

function Spritz_AbsorbNibble(sp::Spritz,x::Int64)
  if sp.a==convert(Int64,floor(sp.N/2)) Spritz_Shuffle(sp) end
  Spritz_Swap(sp,sp.a, (convert(Int64,floor(sp.N/2))+x )%sp.N)

  sp.a = sp.a + 1
end

function Spritz_AbsorbByte(sp::Spritz, b::Int64)
  D=convert(Int64,ceil(sqrt(sp.N)))
  low= b% D
  high  = convert(Int64, floor(b/D))
  # println("!!!!!!!!!! low=", low, ", high=",high, ", D=",D)
  Spritz_AbsorbNibble(sp,low)
  Spritz_AbsorbNibble(sp,high)
end

function Spritz_Absorb(sp::Spritz,I::Array{Int64,2})
#println("size(I)=",size(I))
len=size(I)[2]
#println("len=",len)
for v in 0:(len-1) Spritz_AbsorbByte(sp,I[v+1]) end
end

function Spritz_Output(sp::Spritz)
szk= sp.S[(sp.z+sp.k)%sp.N+1]
si=sp.S[(sp.i+szk)%sp.N+1]
sp.z=sp.S[(sp.j+si)%sp.N+1]
return sp.z
end

function Spritz_Drip(sp::Spritz)
  if sp.a>0 Spritz_Shuffle(sp) end
  Spritz_Update(sp)
  return Spritz_Output(sp)
end

function Spritz_Squeeze(sp::Spritz, r::Int64)
  if sp.a>0 Spritz_Shuffle(sp) end
  P=Int64[]
  for v in 0:(r-1)
     P=[P,Spritz_Drip(sp)]
  end
  return P
end

function Spritz_Initialize(sp,n::Int64)
    sp.i=0
    sp.j=0
    sp.k=0
    sp.z=0
    sp.a=0
    sp.w=1
    sp.N=n
    sp.S=[]
     for v=1:n s
	  sp.S=[sp.S,v-1]
      end
end

function Spritz_KeySetup(sp::Spritz,K::Array{Int64,2})
  Spritz_Absorb(sp,K)
end

function Spritz_Encrypt(n::Int64, K::Array{Int64,2}, M::Array{Int64,2}, debug=0)
#debug = 0, return C
#debug = 1, return permutation S
  sp=Spritz(n)
  Spritz_KeySetup(sp,K)
  #Spritz_Show(sp)
  SQUEEZE=Spritz_Squeeze(sp,size(M)[2])
  C=[]
  for v in 1:size(M)[2]
    C=[C, (M[v]+ SQUEEZE[v])%n]
  end
 # println("C = ", C);
 # print("debug = ", debug,"\t")
  if(debug==0)
    return C'
  else
    return sp.S'
    end

end

function Spritz_Encrypt(K_str::ASCIIString, M_str::ASCIIString)
  n=256
  sp=Spritz(n)

  K=[]

  for v in 1:length(K_str)
        K=[K,int(K_str[v])]
  end

 M=[]
  for v in 1:length(M_str)
        M=[M,int(M_str[v])]
  end

 K=K'
 M=M'
 #println("K=",K," , M=",M, " , size(M)=")
 #println(size(M))

  Spritz_KeySetup(sp,K')
  Spritz_Show(sp)
  SQUEEZE=Spritz_Squeeze(sp,size(M)[2])
  C=[]
  for v in 1:size(M)[2]
    C=[C, char((M[v]+ SQUEEZE[v])%n)]
  end
  #println("C = ", C);
  return C'
end

function Spritz_Decrypt(n::Int64, K::Array{Int64,2}, C::Array{Int64,2})
  println("Hello")
  sp=Spritz(n)
  Spritz_KeySetup(sp,K)
  SQUEEZE=Spritz_Squeeze(sp,size(C)[2])

  M=[]
  for v in 1:size(C)[2]

    if (C[v]-SQUEEZE[v]>=0)
      tmp = (C[v]-SQUEEZE[v])%n
    else
      tmp = (C[v]-SQUEEZE[v]+n)%n
    end
    M = [M, tmp]
  end
  println("M = ",M);
  return M'
end

function Spritz_Hash(n::Int64, M::Array{Int64,2}, r::Int64)
 sp=Spritz(n)
 Spritz_Absorb(sp,M)
 Spritz_AbsorbStop(sp)
 Spritz_AbsorbByte(sp,r)
 R = Spritz_Squeeze(sp,r)
 R = R'
 return R
end

function Spritz_Hash(M_str::ASCIIString, r::Int64)
 n=256
 sp=Spritz(n)

 M=[]
  for v in 1:length(M_str)
        M=[M,int(M_str[v])]
  end

 M=M'

 Spritz_Absorb(sp,M)
 Spritz_AbsorbStop(sp)
 Spritz_AbsorbByte(sp,r)
 SQUEEZED=Spritz_Squeeze(sp,r)
 SQUEEZED=SQUEEZED'

 R=[]
 for v in 1:size(SQUEEZED)[2]
     R = [R, hex(SQUEEZED[v])]
    # print(char(SQUEEZED[v])," ")
  end

     R=R'

 return R
end

max_stars_repo_name: PawelPamula/MarkovTools

max_issues_repo_name: PawelPamula/MarkovTools

max_forks_repo_name: PawelPamula/MarkovTools


text:
#############
# This demonstrates the Chebyshev transform and inverse transform,
# explaining precisely the normalization and points
#############

using FastTransforms

# first kind points -> first kind polynomials
n = 20
p_1 = chebyshevpoints(Float64, n; kind=1)
f = exp.(p_1)
f̌ = chebyshevtransform(f; kind=1)

f̃ = x -> [cos(k*acos(x)) for k=0:n-1]' * f̌
f̃(0.1) ≈ exp(0.1)

# first kind polynomials -> first kind points
ichebyshevtransform(f̌; kind=1) ≈ exp.(p_1)

# second kind points -> first kind polynomials
p_2 = chebyshevpoints(Float64, n; kind=2)
f = exp.(p_2)
f̌ = chebyshevtransform(f; kind=2)

f̃ = x -> [cos(k*acos(x)) for k=0:n-1]' * f̌
f̃(0.1) ≈ exp(0.1)

# first kind polynomials -> second kind points
ichebyshevtransform(f̌; kind=2) ≈ exp.(p_2)

# first kind points -> second kind polynomials
n = 20
p_1 = chebyshevpoints(Float64, n; kind=1)
f = exp.(p_1)
f̌ = chebyshevutransform(f; kind=1)
f̃ = x -> [sin((k+1)*acos(x))/sin(acos(x)) for k=0:n-1]' * f̌
f̃(0.1) ≈ exp(0.1)

# second kind polynomials -> first kind points
ichebyshevutransform(f̌; kind=1) ≈ exp.(p_1)

# second kind points -> second kind polynomials
p_2 = chebyshevpoints(Float64, n; kind=2)[2:n-1]
f = exp.(p_2)
f̌ = chebyshevutransform(f; kind=2)
f̃ = x -> [sin((k+1)*acos(x))/sin(acos(x)) for k=0:n-3]' * f̌
f̃(0.1) ≈ exp(0.1)

# second kind polynomials -> second kind points
ichebyshevutransform(f̌; kind=2) ≈ exp.(p_2)