text:
@testset "prob_macro.jl" begin
    @testset "scalar" begin
        @model function demo(x)
            m ~ Normal()
            x ~ Normal(m, 1)
        end

        mval = 3
        xval = 2
        iters = 1000

        logprior = logpdf(Normal(), mval)
        loglike = logpdf(Normal(mval, 1), xval)
        logjoint = logprior + loglike

        model = demo(xval)
        @test logprob"m = mval | model = model" == logprior
        @test logprob"m = mval | x = xval, model = model" == logprior
        @test logprob"x = xval | m = mval, model = model" == loglike
        @test logprob"x = xval, m = mval | model = model" == logjoint

        varinfo = VarInfo(demo(missing))
        @test logprob"x = xval, m = mval | model = model, varinfo = varinfo" == logjoint

        varinfo = VarInfo(demo(xval))
        @test logprob"m = mval | model = model, varinfo = varinfo" == logprior
        @test logprob"m = mval | x = xval, model = model, varinfo = varinfo" == logprior
        @test logprob"x = xval | m = mval, model = model, varinfo = varinfo" == loglike
    end
    @testset "vector" begin
        n = 5
        @model function demo(x, n = n)
            m ~ MvNormal(n, 1.0)
            x ~ MvNormal(m, 1.0)
        end
        mval = rand(n)
        xval = rand(n)
        iters = 1000

        logprior = logpdf(MvNormal(n, 1.0), mval)
        loglike = logpdf(MvNormal(mval, 1.0), xval)
        logjoint = logprior + loglike

        model = demo(xval)
        @test logprob"m = mval | model = model" == logprior
        @test logprob"x = xval | m = mval, model = model" == loglike
        @test logprob"x = xval, m = mval | model = model" == logjoint

        varinfo = VarInfo(demo(xval))
        @test logprob"m = mval | model = model, varinfo = varinfo" == logprior
        @test logprob"x = xval | m = mval, model = model, varinfo = varinfo" == loglike
        # Currently, we cannot easily pre-allocate `VarInfo` for vector data
    end
    @testset "issue190" begin
        @model function gdemo(x, y)
            s ~ InverseGamma(2, 3)
            m ~ Normal(0, sqrt(s))
            x ~ filldist(Normal(m, sqrt(s)), length(y))
            for i in 1:length(y)
                y[i] ~ Normal(x[i], sqrt(s))
            end
        end
        c = Chains(rand(10, 2), [:m, :s])
        model_gdemo = gdemo([1.0, 0.0], [1.5, 0.0])
        r1 = prob"y = [1.5] | chain=c, model = model_gdemo, x = [1.0]"
        r2 = map(c[:s]) do s
            # exp(logpdf(..)) not pdf because this is exactly what the prob"" macro does, so we test r1 == r2
            exp(logpdf(Normal(1, sqrt(s)), 1.5))
        end
        @test r1 == r2
    end
end


text:
module TestUtils

using Test
using TMLEEpistasis
using MLJ
using TOML
using BGEN
using CSV
using DataFrames

include("helper_fns.jl")

@testset "Test retrieve coefficients" begin
    estimatorfile = joinpath("config", "tmle_config.toml")
    tmle_config = TOML.parsefile(estimatorfile)
    queries = ["query_1" => (t₁=[1, 0], t₂=[1, 0]),]
    tmle = TMLEEpistasis.estimators_from_toml(tmle_config, queries, TMLEEpistasis.PhenotypeTMLEEpistasis)

    n = 1000
    T = (t₁=categorical(rand([1, 0], n)), t₂=categorical(rand([1, 0], n)))
    W = (w₁=rand(n),)
    y = rand(n)
    mach = machine(tmle["continuous"], T, W, y)
    fit!(mach, verbosity=0)

    coefs = TMLEEpistasis.Qstack_coefs(mach)
    @test coefs isa Vector
    @test length(coefs) == 4

    repr = TMLEEpistasis.repr_Qstack_coefs(mach)
    @test repr isa String
    str_reprs = split(repr, ", ")
    modelstr = join([split(v, " => ")[1] for v in str_reprs], ",")
    @test modelstr == "HALRegressor,XGBoostRegressor,XGBoostRegressor,InteractionLMRegressor"

end

@testset "Test parse_queries" begin
    build_query_file()
    queries = TMLEEpistasis.parse_queries(queryfile)
    @test queries == [
        "QUERY_1" => (RSID_10 = ["AG", "GG"], RSID_100 = ["AG", "GG"]),
        "QUERY_2" => (RSID_10 = ["AG", "GG"], RSID_100 = ["AA", "GG"])
    ]
    rm(queryfile)

    @test TMLEEpistasis.querystring(queries[1][2]) == "RSID_10: AG -> GG & RSID_100: AG -> GG"
end

@testset "Test phenotypes parsing" begin
    allnames = TMLEEpistasis.phenotypesnames(phenotypefile)
    @test allnames == [:categorical_phenotype, :continuous_phenotype]

    @test TMLEEpistasis.phenotypes_list(nothing, [:categorical_phenotype], allnames) == [:continuous_phenotype]
    @test TMLEEpistasis.phenotypes_list(nothing, [], allnames) == allnames

    @test TMLEEpistasis.phenotypes_list("data/phen_list_1.csv", [], allnames) == [:continuous_phenotype]
    @test TMLEEpistasis.phenotypes_list("data/phen_list_2.csv", [:continuous_phenotype], allnames) == [:categorical_phenotype]

    # Test init_or_retrieve_results for PhenotypeTMLEEpistasis
    outfile = "temp_results.csv"
    @test TMLEEpistasis.init_or_retrieve_results(outfile, TMLEEpistasis.PhenotypeTMLEEpistasis) == Set()
    init_df = CSV.File(outfile) |> DataFrame
    @test names(init_df) == ["PHENOTYPE",
                            "QUERYNAME",
                            "QUERYSTRING",
                            "ESTIMATE",
                            "PVALUE",
                            "LOWER_BOUND",
                            "UPPER_BOUND",
                            "STD_ERROR",
                            "QSTACK_COEFS",
                            "NROWS",
                            "TCOUNTS"]

    CSV.write(outfile, DataFrame(PHENOTYPE=["toto", "tata", "toto"]))
    @test TMLEEpistasis.init_or_retrieve_results(outfile, TMLEEpistasis.PhenotypeTMLEEpistasis) == Set([:toto, :tata])
    rm(outfile)

    # Test init_or_retrieve_results for PhenotypeCrossValidation
    @test TMLEEpistasis.init_or_retrieve_results(outfile, TMLEEpistasis.PhenotypeCrossValidation) == Set()
    init_df = CSV.File(outfile) |> DataFrame
    @test names(init_df) == ["PHENOTYPE",
                            "VARIANTS",
                            "Q_RESULTSTRING",
                            "G_RESULTSTRING"]
    CSV.write(outfile, DataFrame(PHENOTYPE=["toto", "tata", "toto"]))
    @test TMLEEpistasis.init_or_retrieve_results(outfile, TMLEEpistasis.PhenotypeCrossValidation) == Set([:toto, :tata])
    rm(outfile)
end

@testset "Test TMLE REPORT PARSING" begin
    # Only one query, the reports is just a NamedTuple
    reports = (pvalue=0.05,)
    @test TMLEEpistasis.get_query_report(reports, 4) == reports
    reports = [(pvalue=0.05,), (pvalue=0.01,)]
    @test TMLEEpistasis.get_query_report(reports, 1) == reports[1]
    @test TMLEEpistasis.get_query_report(reports, 2) == reports[2]

end

@testset "Test encode" begin
    # Multiple treatments combined to joint treatment
    T = DataFrame(
        t₁=categorical(["AG", "GG", "AA", "AA", "AG"]),
        t₂=categorical(["CT", "CC", "TT", "CT", "CT"]),
        )

    t_target = TMLEEpistasis.encode(T)
    @test levels(t_target) == [1, 2, 9, 3, 5, 6, 8, 7, 4]
    @test t_target == categorical([5, 3, 7, 4, 5])

    # One treatment just returned
    T = DataFrame(t₁=categorical(["AG", "GG", "AA", "AA", "AG"]))

    t_target = TMLEEpistasis.encode(T)
    @test t_target == ["AG", "GG", "AA", "AA", "AG"]
end

@testset "Test TreatmentCountsRepr" begin
    T = DataFrame(
        RS_100=categorical(["AA", "AG", "GG", "AA", "AG"]),
        RS_10=categorical(["CC", "CG", "CG", "CC", "GG"])
        )

    @test TMLEEpistasis.TreatmentCountsRepr(T) ==
        "RS_100 & RS_10: AA CC 2 | AG CG 1 | AG GG 1 | GG CG 1"

    T = DataFrame(
        RS_100=categorical(["AA", "AG", "GG", "AA", "AG"]),
    )
    @test TMLEEpistasis.TreatmentCountsRepr(T) == "RS_100: AA 2 | AG 2 | GG 1"

end

end;

true

ext: jl


text:
function test_co_occurrence()
    println("-- Testing CoOccurrence recommender")

    data = DataAccessor(sparse([1 0 0; 4 5 0]))
    recommender = CoOccurrence(data, 1)
    build!(recommender)
    @test ranking(recommender, 1, 1) == 100.0
    @test ranking(recommender, 1, 2) == 50.0
    @test ranking(recommender, 1, 3) == 0.0
end

test_co_occurrence()


text:
println("Hello World!");

# https://github.com/docker-library/julia/pull/6
download("https://google.com");

max_stars_repo_name: benbc/docker-official-images

max_issues_repo_name: benbc/docker-official-images

max_forks_repo_name: benbc/docker-official-images


text:
# UFixed{T,f} maps UInts from 0 to 2^f-1 to the range [0.0, 1.0]
# For example, a UFixed8 maps 0x00 to 0.0 and 0xff to 1.0

immutable UFixed{T<:Unsigned,f} <: FixedPoint{T,f}
    i::T

    UFixed(i::Integer,_) = new(i%T)   # for setting by raw representation
    UFixed(x) = convert(UFixed{T,f}, x)
end

  rawtype{T,f}(::Type{UFixed{T,f}}) = T
  rawtype(x::Number) = rawtype(typeof(x))
nbitsfrac{T,f}(::Type{UFixed{T,f}}) = f
floattype{T<:UFixed}(::Type{T}) = floattype(supertype(T))

typealias UFixed8  UFixed{UInt8,8}
typealias UFixed10 UFixed{UInt16,10}
typealias UFixed12 UFixed{UInt16,12}
typealias UFixed14 UFixed{UInt16,14}
typealias UFixed16 UFixed{UInt16,16}

const UF = (UFixed8, UFixed10, UFixed12, UFixed14, UFixed16)

reinterpret{T<:Unsigned, f}(::Type{UFixed{T,f}}, x::T) = UFixed{T,f}(x, 0)

# The next lines mimic the floating-point literal syntax "3.2f0"
immutable UFixedConstructor{T,f} end
*{T,f}(n::Integer, ::UFixedConstructor{T,f}) = UFixed{T,f}(n,0)
const uf8  = UFixedConstructor{UInt8,8}()
const uf10 = UFixedConstructor{UInt16,10}()
const uf12 = UFixedConstructor{UInt16,12}()
const uf14 = UFixedConstructor{UInt16,14}()
const uf16 = UFixedConstructor{UInt16,16}()

zero{T,f}(::Type{UFixed{T,f}}) = UFixed{T,f}(zero(T),0)
function one{T<:UFixed}(::Type{T})
    T(typemax(rawtype(T)) >> (8*sizeof(T)-nbitsfrac(T)), 0)
end
zero(x::UFixed) = zero(typeof(x))
 one(x::UFixed) =  one(typeof(x))
rawone(v) = reinterpret(one(v))

# Conversions
convert{T<:UFixed}(::Type{T}, x::T) = x
convert{T1,T2,f}(::Type{UFixed{T1,f}}, x::UFixed{T2,f}) = UFixed{T1,f}(convert(T1, x.i), 0)
function convert{T,T2,f}(::Type{UFixed{T,f}}, x::UFixed{T2})
    U = UFixed{T,f}
    y = round((rawone(U)/rawone(x))*reinterpret(x))
    (0 <= y) & (y <= typemax(T)) || throw_converterror(U, x)
    reinterpret(U, _unsafe_trunc(T, y))
end
convert(::Type{UFixed16}, x::UFixed8) = reinterpret(UFixed16, convert(UInt16, 0x0101*reinterpret(x)))
convert{U<:UFixed}(::Type{U}, x::Real) = _convert(U, rawtype(U), x)
function _convert{U<:UFixed,T}(::Type{U}, ::Type{T}, x)
    y = round(widen1(rawone(U))*x)
    (0 <= y) & (y <= typemax(T)) || throw_converterror(U, x)
    U(_unsafe_trunc(T, y), 0)
end
function _convert{U<:UFixed}(::Type{U}, ::Type{UInt128}, x)
    y = round(rawone(U)*x)   # for UInt128, we can't widen
    (0 <= y) & (y <= typemax(UInt128)) & (x <= Float64(typemax(U))) || throw_converterror(U, x)
    U(_unsafe_trunc(UInt128, y), 0)
end

rem{T<:UFixed}(x::T, ::Type{T}) = x
rem{T<:UFixed}(x::UFixed, ::Type{T}) = reinterpret(T, _unsafe_trunc(rawtype(T), round((rawone(T)/rawone(x))*reinterpret(x))))
rem{T<:UFixed}(x::Real, ::Type{T}) = reinterpret(T, _unsafe_trunc(rawtype(T), round(rawone(T)*x)))

# convert(::Type{AbstractFloat}, x::UFixed) = convert(floattype(x), x)
float(x::UFixed) = convert(floattype(x), x)

convert(::Type{BigFloat}, x::UFixed) = reinterpret(x)*(1/BigFloat(rawone(x)))
function convert{T<:AbstractFloat}(::Type{T}, x::UFixed)
    y = reinterpret(x)*(one(rawtype(x))/convert(T, rawone(x)))
    convert(T, y)  # needed for types like Float16 which promote arithmetic to Float32
end
convert(::Type{Bool}, x::UFixed) = x == zero(x) ? false : true
convert{T<:Integer}(::Type{T}, x::UFixed) = convert(T, x*(1/one(T)))
convert{Ti<:Integer}(::Type{Rational{Ti}}, x::UFixed) = convert(Ti, reinterpret(x))//convert(Ti, rawone(x))
convert(::Type{Rational}, x::UFixed) = reinterpret(x)//rawone(x)

# Traits
abs(x::UFixed) = x

(-){T<:UFixed}(x::T) = T(-reinterpret(x), 0)
(~){T<:UFixed}(x::T) = T(~reinterpret(x), 0)

+{T,f}(x::UFixed{T,f}, y::UFixed{T,f}) = UFixed{T,f}(convert(T, x.i+y.i),0)
-{T,f}(x::UFixed{T,f}, y::UFixed{T,f}) = UFixed{T,f}(convert(T, x.i-y.i),0)
*{T<:UFixed}(x::T, y::T) = convert(T,convert(floattype(T), x)*convert(floattype(T), y))
/{T<:UFixed}(x::T, y::T) = convert(T,convert(floattype(T), x)/convert(floattype(T), y))

# Comparisons
 <{T<:UFixed}(x::T, y::T) = reinterpret(x) < reinterpret(y)
<={T<:UFixed}(x::T, y::T) = reinterpret(x) <= reinterpret(y)

# Functions
trunc{T<:UFixed}(x::T) = T(div(reinterpret(x), rawone(T))*rawone(T),0)
floor{T<:UFixed}(x::T) = trunc(x)
function round{T,f}(x::UFixed{T,f})
    mask = convert(T, 1<<(f-1))
    y = trunc(x)
    return convert(T, reinterpret(x)-reinterpret(y)) & mask>0 ?
            UFixed{T,f}(y+one(UFixed{T,f})) : y
end
function ceil{T,f}(x::UFixed{T,f})
    k = 8*sizeof(T)-f
    mask = (typemax(T)<<k)>>k
    y = trunc(x)
    return convert(T, reinterpret(x)-reinterpret(y)) & (mask)>0 ?
            UFixed{T,f}(y+one(UFixed{T,f})) : y
end

trunc{T<:Integer}(::Type{T}, x::UFixed) = convert(T, div(reinterpret(x), rawone(x)))
round{T<:Integer}(::Type{T}, x::UFixed) = round(T, reinterpret(x)/rawone(x))
floor{T<:Integer}(::Type{T}, x::UFixed) = trunc(T, x)
 ceil{T<:Integer}(::Type{T}, x::UFixed) =  ceil(T, reinterpret(x)/rawone(x))

isfinite(x::UFixed) = true
isnan(x::UFixed) = false
isinf(x::UFixed) = false

bswap{f}(x::UFixed{UInt8,f}) = x
bswap(x::UFixed)  = typeof(x)(bswap(reinterpret(x)),0)

function minmax{T<:UFixed}(x::T, y::T)
    a, b = minmax(reinterpret(x), reinterpret(y))
    T(a,0), T(b,0)
end

# Iteration
# The main subtlety here is that iterating over 0x00uf8:0xffuf8 will wrap around
# unless we iterate using a wider type
@inline start{T<:UFixed}(r::StepRange{T}) = widen1(reinterpret(r.start))
@inline next{T<:UFixed}(r::StepRange{T}, i::Integer) = (T(i,0), i+reinterpret(r.step))
@inline function done{T<:UFixed}(r::StepRange{T}, i::Integer)
    i1, i2 = reinterpret(r.start), reinterpret(r.stop)
    isempty(r) | (i < min(i1, i2)) | (i > max(i1, i2))
end

function decompose(x::UFixed)
    g = gcd(reinterpret(x), rawone(x))
    div(reinterpret(x),g), 0, div(rawone(x),g)
end

# Promotions
promote_rule{T<:UFixed,Tf<:AbstractFloat}(::Type{T}, ::Type{Tf}) = promote_type(floattype(T), Tf)
promote_rule{T<:UFixed, R<:Rational}(::Type{T}, ::Type{R}) = R
function promote_rule{T<:UFixed, Ti<:Union{Signed,Unsigned}}(::Type{T}, ::Type{Ti})
    floattype(T)
end
@generated function promote_rule{T1,T2,f1,f2}(::Type{UFixed{T1,f1}}, ::Type{UFixed{T2,f2}})
    f = max(f1, f2)  # ensure we have enough precision
    T = promote_type(T1, T2)
    # make sure we have enough integer bits
    i1, i2 = 8*sizeof(T1)-f1, 8*sizeof(T2)-f2  # number of integer bits for each
    i = 8*sizeof(T)-f
    while i < max(i1, i2)
        T = widen1(T)
        i = 8*sizeof(T)-f
    end
    :(UFixed{$T,$f})
end

_unsafe_trunc{T}(::Type{T}, x::Integer) = x % T
_unsafe_trunc{T}(::Type{T}, x)          = unsafe_trunc(T, x)


text:
using LinearAlgebra

@testset "Test NCU algorithm" begin
    @testset "Test default NCU with 2CX" begin
        cct1 = QXZoo.Circuit.Circ(3)
        cct2 = QXZoo.Circuit.Circ(3)
        QXZoo.NCU.apply_ncu!(cct1, collect(1:2), [], 3, DefaultGates.GateSymbols.c_x)

        cct2 << GateCall2(GateSymbol(:c_x, 1), 3, 2, nothing)
        cct2 << GateCall2(DefaultGates.GateSymbols.c_x, 2, 1, DefaultGates.GateSymbols.x)
        cct2 << GateCall2(GateSymbol(:c_x, 1, true), 3, 2, nothing)
        cct2 << GateCall2(DefaultGates.GateSymbols.c_x, 2, 1, DefaultGates.GateSymbols.x)
        cct2 << GateCall2(GateSymbol(:c_x, 1), 3, 1, nothing)

        @test cct1.circ_ops.len == cct2.circ_ops.len
        @test cct1.gate_set == cct2.gate_set

        h1 = cct1.circ_ops.head
        h2 = cct2.circ_ops.head
        for i in 1:cct1.circ_ops.len
            @test h1.data === h2.data
            h1 = h1.next
            h2 = h2.next
        end
    end

    @testset "Test default NCU with 3CX" begin
        cct1 = QXZoo.Circuit.Circ(4)
        cct2 = QXZoo.Circuit.Circ(4)
        QXZoo.NCU.apply_ncu!(cct1, collect(1:3), [], 4, DefaultGates.GateSymbols.c_x)

        cct2 << GateCall2(GateSymbol(:c_x, 2), 4, 1, nothing)
        cct2 << DefaultGates.c_x(2,1)
        cct2 << GateCall2(GateSymbol(:c_x, 2, true), 4, 2, nothing)
        cct2 << DefaultGates.c_x(2,1)
        cct2 << GateCall2(GateSymbol(:c_x, 2), 4, 2, nothing)
        cct2 << DefaultGates.c_x(3,2)
        cct2 << GateCall2(GateSymbol(:c_x, 2, true), 4, 3, nothing)
        cct2 << DefaultGates.c_x(3,1)
        cct2 << GateCall2(GateSymbol(:c_x, 2), 4, 3, nothing)
        cct2 << DefaultGates.c_x(3,2)
        cct2 << GateCall2(GateSymbol(:c_x, 2, true), 4, 3, nothing)
        cct2 << DefaultGates.c_x(3,1)
        cct2 << GateCall2(GateSymbol(:c_x, 2), 4, 3, nothing)

        @test cct1.circ_ops.len == cct2.circ_ops.len
        @test cct1.gate_set == cct2.gate_set

        h1 = cct1.circ_ops.head
        h2 = cct2.circ_ops.head
        for i in 1:cct1.circ_ops.len
            @test h1.data === h2.data
            h1 = h1.next
            h2 = h2.next
        end
    end

    @testset "Test aux. assisted NCU with 5CZ" begin
        cct1 = QXZoo.Circuit.Circ(9)
        cct2 = QXZoo.Circuit.Circ(9)
        QXZoo.NCU.apply_ncu!(cct1, collect(1:5), collect(6:8), 9, DefaultGates.GateSymbols.c_z)

        cct2 << DefaultGates.h(9)
        QXZoo.NCU.apply_ncu!(cct2, [5,8], [], 9, DefaultGates.GateSymbols.c_x)
        QXZoo.NCU.apply_ncu!(cct2, [4,7], [], 8, DefaultGates.GateSymbols.c_x)
        QXZoo.NCU.apply_ncu!(cct2, [3,6], [], 7, DefaultGates.GateSymbols.c_x)
        QXZoo.NCU.apply_ncu!(cct2, [1,2], [], 6, DefaultGates.GateSymbols.c_x)
        QXZoo.NCU.apply_ncu!(cct2, [3,6], [], 7, DefaultGates.GateSymbols.c_x)
        QXZoo.NCU.apply_ncu!(cct2, [4,7], [], 8, DefaultGates.GateSymbols.c_x)
        QXZoo.NCU.apply_ncu!(cct2, [5,8], [], 9, DefaultGates.GateSymbols.c_x)
        QXZoo.NCU.apply_ncu!(cct2, [4,7], [], 8, DefaultGates.GateSymbols.c_x)
        QXZoo.NCU.apply_ncu!(cct2, [3,6], [], 7, DefaultGates.GateSymbols.c_x)
        QXZoo.NCU.apply_ncu!(cct2, [1,2], [], 6, DefaultGates.GateSymbols.c_x)
        QXZoo.NCU.apply_ncu!(cct2, [3,6], [], 7, DefaultGates.GateSymbols.c_x)
        QXZoo.NCU.apply_ncu!(cct2, [4,7], [], 8, DefaultGates.GateSymbols.c_x)
        cct2 << DefaultGates.h(9)

        @test cct1.circ_ops.len == cct2.circ_ops.len
        @test cct1.gate_set == cct2.gate_set

        h1 = cct1.circ_ops.head
        h2 = cct2.circ_ops.head
        for i in 1:cct1.circ_ops.len
            @test h1.data === h2.data
            h1 = h1.next
            h2 = h2.next
        end
    end

end


text:
using Test
using SymDiff

@test simplify(Val{:^}, :((2+0)^(1*a))) == :(2^a)


text:
function vec2scam(v::AbstractVector)
    Camera(SVector{3}(v[1:3]), SVector{3}(v[4:6]), v[7], SVector{2}(v[8:9]), SVector{2}(v[10:11]))
end

##################### objective #############################

function rodrigues_rotate_point(rot::AbstractVector{T}, X::AbstractVector{T}) where T
    sqtheta = sum(rot.*rot)
    if sqtheta > 1e-10
        theta = sqrt(sqtheta)
        costheta = cos(theta)
        sintheta = sin(theta)
        theta_inverse = 1. / theta

        w = theta_inverse * rot
        w_cross_X = cross(w,X)
        tmp = dot(w,X) * (1. - costheta)

        X*costheta + w_cross_X * sintheta + w * tmp
    else
        X + cross(rot,X)
    end
end

@inline function radial_distort(rad_params,proj)
    rsq = sum(proj.*proj)
    L = 1. + rad_params[1]*rsq + rad_params[2]*rsq*rsq
    proj*L
end

function project(cam, X)
    Xcam = rodrigues_rotate_point(cam.r, X - cam.c)
    distorted = radial_distort(cam.κ, SVector(Xcam[1], Xcam[2])/Xcam[3])
    distorted*cam.f + cam.x0
end

function compute_reproj_err(cam, X, w, feat)
    return w*(project(cam,X) - feat)
end

function pack(cam,X,w)
    vcat(cam.r, cam.c, cam.f, cam.x0, cam.κ, X, w)
end

function unpack(packed)
    N = length(packed)
    vec2scam(view(packed,1:N-4)), view(packed,N-3:N-1),packed[N]
end

function compute_w_err(w)
    1.0 - w*w
end

function compute_reproj_err_d(params, feat)
    cam, X, w = unpack(params)
    compute_reproj_err(cam,X,w,feat)
end

function compute_ba_J(::Val{:ForwardDiff}, cams::AbstractArray{T},X,w,obs,feats) where T
    p = size(obs,2)
    reproj_err_d = zeros(2*p, 11 + 3 + 1)
    for i in 1:p
        idx = (2*(i-1))+1
        @inbounds j = obs[1,i]
        @inbounds l = obs[2,i]
        @inbounds xi = SVector{15,T}(cams[:,j]...,X[:,l]...,w[i])
        @inbounds reproj_err_d[idx:idx+1,:] = ForwardDiff.jacobian(
            x -> compute_reproj_err_d(x, SVector{2,T}(feats[1,i], feats[2,i])),
              xi)
    end
    w_err_d = zeros(1,p)
    for i in 1:p
        @inbounds w_err_d[i] = ForwardDiff.derivative(compute_w_err, w[i])
    end
    (reproj_err_d, w_err_d)
end

function ba_objective(cams, X, w, obs, feats)
    reproj_err = similar(feats)
    for i in 1:size(feats,2)
        @inbounds l = obs[2,i]
        @inbounds reproj_err[:,i] .= compute_reproj_err(cams[obs[1,i]],SVector(X[1,l], X[2,l], X[3,l]),w[i],SVector(feats[1,i], feats[2,i]))
    end
    w_err = 1.0 .- w .* w
    (reproj_err, w_err)
end


text:
module NEPTransformations

    # Transformations between types
    using SparseArrays

    include("nleigs_coefficients.jl")

    export taylor_expansion_pep
    export shift_and_scale
    export mobius_transform
    export CORKPencil
    export buildPencil
    export CorkLinearization
    export DefaultCorkLinearization
    export IarCorkLinearization
    export NleigsCorkLinearization
    export CORKPencilLR
    export lowRankCompress

    # We overload these
    import ..NEPCore.compute_Mder
    import ..NEPCore.compute_Mlincomb
    import ..NEPCore.compute_Mlincomb!
    import ..NEPCore.compute_MM
    import ..NEPCore.compute_resnorm
    import Base.size
    import SparseArrays.issparse

    """
        struct ShiftScaledNEP <: NEP
        ShiftScaleNEP(orgnep::NEP[,shift=0][,scale=1])

    Transforms a nep (orgnep) M(λ)v to a new nep T(λ)=M(scale*λ+shift). This can be used if the method does not have an easy implementation of shift and scaling. Usage of this transformation can slow down the algorithm.
    """
    struct ShiftScaledNEP <: NEP
        shift::Number
        scale::Number
        orgnep::NEP
        function ShiftScaledNEP(orgnep;shift=0,scale=1)
            return new(shift,scale,orgnep)
        end
    end
    """
        shift_and_scale(orgnep::NEP;shift=0,scale=1)

    Transforms the orgnep by defining a new NEP from the relation
    T(λ)=M(scale * λ+shift) where M is the orgnep. This function tries
     to preserve the NEP type, e.g., a shift_and_scale operation on
    an SPMF-object, return an SPMF object. If it cannot preserve
    the type, it will return a nep of the struct `ShiftScaledNEP`.

    # Example
    ```julia-repl
    julia> nep0=nep_gallery("pep0")
    julia> σ=3; α=10;
    julia> nep1=shift_and_scale(nep0,shift=σ,scale=α)
    julia> opnorm(compute_Mder(nep0,α*(4+4im)+σ)-compute_Mder(nep1,4+4im))
    8.875435870738592e-12
    ```
    """
    function shift_and_scale(orgnep::NEP;shift=0,scale=1)
        return ShiftScaledNEP(orgnep,shift=shift,scale=scale);
    end

    # Native implementation for PEPs by transformation of the coefficients. This is not necessarily fast.
    function shift_and_scale(orgnep::PEP;shift=0,scale=1)
        At=copy(shift*scale*orgnep.A); # Allocate correct type (not so fast)
        m=size(At,1)-1
        for j=0:m
            AA = zero(At[j+1])
            for i=j:m
                factor=((scale^j)*(shift^(i-j))*factorial(i)/factorial(i-j))/factorial(j)
                AA+=orgnep.A[i+1]*factor;

            end
            At[j+1]=AA;
        end
        return PEP(At)
    end

    # A shift and rescale of a DEP is DEP where the matrix coefficients and the delays are rescaled and another matrix coefficient with delay zero is added.
    function shift_and_scale(orgnep::DEP;shift=0,scale=1)
        J = Matrix{eltype(orgnep.A[1])}(I, size(orgnep.A[1]))
        return DEP([broadcast(*,orgnep.A,exp.(-orgnep.tauv*shift)/scale); [-shift/scale * J] ],[orgnep.tauv*scale; 0])
    end

    # Native implementation for SPMF. Only for generic SPMF (not AbstractSPMF)
    function shift_and_scale(orgnep::SPMF_NEP;shift=0,scale=1)
        orgfv=get_fv(orgnep);
        m=size(orgfv,1);
        fv = Vector{Function}(undef, m)
        # create anonymous functions corresponding to the
        # shift and scaled problem
        for i=1:m
            fv[i] = S -> orgfv[i](scale*S + shift*one(S))::((S isa Number) ? Number : Matrix)
        end
        # create a new SPMF with transformed functions
        return SPMF_NEP(get_Av(orgnep),fv;Schur_fact=orgnep.Schur_factorize_before);
    end

    # Size and issparse implemented by delegation
    function size(nep::ShiftScaledNEP,dim=-1)
        return size(nep.orgnep,dim)
    end
    function issparse(nep::ShiftScaledNEP)
        return issparse(nep.orgnep)
    end

    function compute_Mder(nep::ShiftScaledNEP,λ::Number,i::Integer=0)
        # Just scale the orgnep call. Chain rule for differentiation.
        return (nep.scale^i)*compute_Mder(nep.orgnep,nep.scale*λ+nep.shift,i);
    end

    function compute_MM(nep::ShiftScaledNEP,S,V)
        # Just call orgnep with a different S
        return compute_MM(nep.orgnep,S*nep.scale+eye(S)*nep.shift,V)
    end

    function compute_Mlincomb(nep::ShiftScaledNEP,λ::Number,V::AbstractVecOrMat,
                              a::Vector=ones(size(V,2)))
        # Multiply the coefficient matrix V with a diagonal matrix
        # corresponding to the chain rule
        p=size(V,2);
        z=nep.scale.^Array{eltype(V),1}(0:p-1)
        W=V * diagm(0 => z) # not very fast
        return compute_Mlincomb(nep.orgnep,nep.scale*λ+nep.shift,W,a);
    end

    """
        mobius_transform(orgnep::NEP,[,a=1][,b=0][,c=0][,d=1])

    Transforms a nep (orgnep) ``M(λ)v`` to a new nep ``T(λ)=M((aλ+b)/(cλ+d))``.
    This function tries to preserve the type such that `T`
    and `M` are of the same NEP-type (see `shift_and_scale()`).
    If it cannot be preserved it will return a `MobiusTransformedNEP`.
    It is in general advised to try to preserve the type, and
    the use of `MobiusTransformedNEP` can
    considerably slow down NEP-access.

    # Example

    ```julia-repl
    julia> nep0=nep_gallery("pep0")
    julia> a=1; b=3; c=4; d=5;
    julia> nep1=mobius_transform(nep0,a=a,b=b,c=c,d=d);
    julia> s=3;
    julia> opnorm(compute_Mder(nep0,(a*s+b)/(c*s+d))-compute_Mder(nep1,s))
    0.0
    ```
    """
    function mobius_transform(orgnep::NEP;a=1,b=0,c=0,d=1)
        return MobiusTransformedNEP(orgnep;a=a,b=b,c=c,d=d)
    end

    function mobius_transform(orgnep::SPMF_NEP;a=1,b=0,c=0,d=1)
        orgfv=get_fv(orgnep);
        m=size(orgfv,1);
        fv = Vector{Function}(undef, m)
        # create anonymous functions corresponding to the
        # möbius transformed problem
        for i=1:m
            fv[i] = S -> orgfv[i]((a*S + b*one(S)) / (c*S + d*one(S)))::((S isa Number) ? Number : Matrix)
        end
        # create a new SPMF with transformed functions
        return SPMF_NEP(get_Av(orgnep),fv;Schur_fact=orgnep.Schur_factorize_before);
    end

    """
    Hidden type representing a transformed NEP
    """
    struct MobiusTransformedNEP <: NEP
        a::Number
        b::Number
        c::Number
        d::Number
        orgnep::NEP
        function MobiusTransformedNEP(orgnep;a=1,b=0,c=0,d=1)
            return new(a,b,c,d,orgnep)
        end
    end

    # Size and issparse implemented by delegation
    function size(nep::MobiusTransformedNEP,dim=-1)
        return size(nep.orgnep,dim)
    end
    function issparse(nep::MobiusTransformedNEP)
        return issparse(nep.orgnep)
    end

    function compute_Mder(nep::MobiusTransformedNEP,λ::Number,i::Integer=0)
        # I did not found a better way
        return compute_Mder_from_MM(nep,λ,i);
    end

    function compute_MM(nep::MobiusTransformedNEP,S,V)
        # Just call orgnep with a different S
        return compute_MM(nep.orgnep, (nep.a*S + nep.b*I) / (nep.c*S + nep.d*I), V)
    end

    function compute_Mlincomb(nep::MobiusTransformedNEP,λ::Number,V::AbstractVecOrMat
                              ,a::Vector=ones(size(V,2)))
        # I did not found a better way
        return compute_Mlincomb_from_MM(nep,λ,V,a)
    end

    # consider renaming this function
    """
        taylor_expansion_pep(orgnep::NEP[,d=2])

    Compute the truncated (with `d` terms) Taylor series of the NEP.
    The output is a [`PEP`](@ref).
    """
    function taylor_expansion_pep(nep::NEP,d::Integer=2)
        A=Vector{Matrix}(undef,d+1)
        for i=0:d
            A[i+1]=compute_Mder(nep,0,i)/factorial(i);
        end
        return PEP(A)
    end

"""
    struct CORKPencil
    function CORKPencil(M,N,Av,Bv,Z);
    function CORKPencil(nep,is)

The struct `CORKPencil` represents a pencil with a particular structure,
as given in the reference. The data can either be constructed directly
via the first constructor, or from a NEP in the second constructor.
The second constructor takes a `NEP` and  `is` which specifies
a CORK-structure as well as an approximation method. This can be objects
of the type `IarCorkLinearization` or `NleigsCorkLinearization`,
which are the CORK-linearizations equivalent to (certain versions of)
[`iar`](@ref) and [`nleigs`](@ref).

See [`buildPencil`](@ref) how to build standard pencil.

# Example:

The following example constructs a  `CORKPencil` from
a general NEP and then computes approximations of NEPs
by the interpolation approach of `nleigs`.

```julia-repl
julia> using LinearAlgebra
julia> A=(1:4)*(1:4)'+I; B=diagm(1 => [1,2,3]); C=ones(4,4);
julia> f1= λ-> one(λ);
julia> f2= λ-> λ;
julia> f3= λ-> exp(sin(λ/2));
julia> nep=SPMF_NEP([A,B,C],[f1,f2,f3]);
julia> cp=CORKPencil(nep,NleigsCorkLinearization());
julia> (A,B)=buildPencil(cp) # Form the pencil
julia> λv=eigen(A,B).values;
julia> λ=λv[sortperm(abs.(λv))[1]]; # Smallest eigval
julia> minimum(svdvals(compute_Mder(nep,λ))) # It's a solution
2.4364382475487156e-11
```

# References:

*  Compact rational Krylov methods for nonlinear eigenvalue problems SIAM Journal on Matrix Analysis and Applications, 36 (2), 820-838, 2015.

    """
struct CORKPencil{T1<:AbstractMatrix,T2<:AbstractMatrix}
        M::T1
        N::T1
        Av::Vector{T2}   # Array of Array of matrices
        Bv::Vector{T2}   # Array of Array of matrices
    end

    abstract type CorkLinearization end;

    struct DefaultCorkLinearization <: CorkLinearization end;
    # TODO: define default

    struct IarCorkLinearization <: CorkLinearization
        d::Int
        function IarCorkLinearization(;d=10)
            return new(d);
        end
    end

    struct NleigsCorkLinearization <: CorkLinearization
        Σ::Vector
        Ξ::Vector
        maxdgr::Int
        tollin::Float64;
        function NleigsCorkLinearization(;Σ=[-1.0-1im,-1+1im,+1+1im,1-1im],Ξ = [Inf], maxdgr=100 , tollin=1e-6 )
            return new(Σ,Ξ,maxdgr,tollin);
        end
    end

    function CORKPencil(nep,is::IarCorkLinearization)
        M=diagm( 0 =>  ones(is.d) )[2:end,:]
        N=diagm( -1 =>  1 ./ (1:is.d-1) )[2:end,:]
        Av=Array{AbstractMatrix,1}(undef, is.d)
        Bv=Array{AbstractMatrix,1}(undef, is.d)
        Av[1]=-compute_Mder(nep,0,0)
        for j=2:is.d Av[j]=zero(Av[1])              end
        for j=1:is.d Bv[j]=compute_Mder(nep,0,j)/j  end
        return CORKPencil(M,N,Av,Bv)
    end

    function CORKPencil(nep,is::NleigsCorkLinearization)
        D,β,ξ,σ=nleigs_coefficients(nep,is.Σ,tollin=is.tollin,Ξ=is.Ξ)
        d=length(β)-1
        σ=σ[1:d+1]; β=β[1:d+1]; ξ=ξ[1:d+1]
        σ,β,ξ=promote(σ,β,ξ) # be sure that M and N have the same type
        M=diagm( -1 => σ[1:d], 0 =>  β[1:d] )[2:end-1,1:end-1]
        N=diagm( -1 => ones(d), 0 =>  β[1:d]./ξ[1:d] )[2:end-1,1:end-1]
        Av=Array{AbstractMatrix,1}(undef, d)
        Av[1:d-1]=D[1:d-1]; Av[d]=D[d]-σ[d]/β[d+1]*D[d+1]
        Bv=Array{AbstractMatrix,1}(undef, d)
        Bv[1:d-1]=D[1:d-1]/ξ[d+1]; Bv[d]=D[d]/ξ[d+1]-D[d+1]/β[d+1]
        return CORKPencil(M,N,Av,Bv)
    end

    # construct the linearization
"""
    (A,B)=buildPencil(cp)

Constructs a pencil from a `CORKPencil` or `CORKPencilLR`. The returned
matrices correspond to a generalized eigenvalue problem. See also [`CORKPencil`](@ref), [`CORKPencilLR`](@ref). See [`lowRankCompress`] for examples.

"""
    function buildPencil(cp::CORKPencil)
        n=size(cp.Av[1],1)
        if issparse(cp.Av[1])   II=sparse(I, n, n)
        else                    II=Matrix(I, n, n)    end
        return [hcat(cp.Av...); kron(cp.M,II)], [hcat(cp.Bv...); kron(cp.N,II)]
    end

"""
    struct CORKPencilLR
    function CORKPencilLR(M,N,Av,AvLR,Bv,BvLR,Z);

Represents / constructs a low-rank CORK-pencil. `AvLR`, `BvLR`
and `Z` correspond to the low-rank factorization of terms
`Av` and `Bv`. See [`CORKPencil`](@ref) and reference below.

# Example

The example illustrate a low-rank linearization of a Taylor expansion
of the NEP ``A0-λI+vv^Te^{-λ}``.

```julia-repl
julia> A0=[1.0 3.0; -1.0 2.0]/10;
julia> v=reshape([-1.0 ; 1]/sqrt(2),n,1);

julia> Av=[-A0-v*v']
julia> Bv=[-one(A0)-v*v']
julia> BvLR=[v/2, -v/3, v/4, -v/5, v/6, -v/7,  v/8, -v/9]
julia> AvLR=zero.(BvLR);
julia> Z=v;
julia> d=9;
julia> M=diagm( 0 =>  ones(d) )[2:end,:]
julia> N=diagm( -1 =>  1 ./ (1:d-1) )[2:end,:]
julia> cplr=CORKPencilLR(M,N,Av,AvLR,Bv,BvLR,Z);
julia> (AA,BB)=buildPencil(cplr);
julia> λ=eigen(AA,BB).values[end];
julia> minimum(svdvals(A0-λ*I+v*v'*exp(-λ)))
8.870165379112754e-13
```

# References:

*  Section 7 in Compact rational Krylov methods for nonlinear eigenvalue problems SIAM Journal on Matrix Analysis and Applications, 36 (2), 820-838, 2015.
"""
    struct CORKPencilLR
    	M::AbstractMatrix
    	N::AbstractMatrix
    	Av::Vector{AbstractMatrix}
    	AvLR::Vector{AbstractMatrix} # In CORK-paper Avtilde
    	Bv::Vector{AbstractMatrix}
    	BvLR::Vector{AbstractMatrix} # In CORK-paper Bvtilde
        Z::AbstractMatrix
    end

"""
    cplr=lowRankCompress(cp_org::CORKPencil,dtilde,rk)

Constructs a [`CORKPencilLR`](@ref) from a [`CORKPencil`](@ref). This is done by
assuming that terms higher than `dtilde` are of low rank, with rank `rk`.
More precisely, all `A[j]` and `B[j]` for `j>dtilde` are assumed
to be of the form ``C_jZ^T``.

# Example

This illustrates how to form a `CORKPencil` from a NEP, and
subsequently form a smaller pencil using `CORKPencilLR`.

```julia-repl
julia> A0=[1.0 3.0; -1.0 2.0]/10;
julia> v=[-1.0 ; 1]/sqrt(2);
julia> nep=DEP([A0,v*v']);
julia> cp_org=CORKPencil(nep,IarCorkLinearization(d=10));
julia> cplr=lowRankCompress(cp_org::CORKPencil,1,1);
julia> (AA,BB)=buildPencil(cplr);
julia> λ=eigen(AA,BB).values[end];
julia> minimum(svdvals(A0-λ*I+v*v'*exp(-λ))) # Check if it is an eigval
2.7621446071952216e-14
```

"""
function lowRankCompress(cp_org::CORKPencil,dtilde,rk)
    d=length(cp_org.Av);
    # Only use the first low-rank term as a basis (this
    # can cause if the first low-rank term has lower rank
    # than the other ones).
    Z=svd(cp_org.Bv[dtilde+1]).V[:,1:rk];

    # Check that the M and N-matrices have the block triangular structure
    if ((norm(cp_org.M[1:(dtilde-1),(dtilde+1):end]) > 0) ||
        (norm(cp_org.N[1:(dtilde-1),(dtilde+1):end]) > 0))
        error("The M-matrix does not have the required structure. Try increasing dtilde.");
    end

    # Manually extract the low-rank coefficients by
    # multiplying by Z. Works since Z is orthogonal.
    Bvtilde=map(i-> cp_org.Bv[i]*Z,  (dtilde+1):d);
    Avtilde=map(i-> cp_org.Av[i]*Z,  (dtilde+1):d);

    return CORKPencilLR(cp_org.M,cp_org.N,
                        cp_org.Av[1:dtilde],Avtilde,
                        cp_org.Bv[1:dtilde],Bvtilde,
                        Z);
end
function buildPencil(cp::CORKPencilLR)

    # Extract the variables
    n=size(cp.Av[1],1)
    dtilde=length(cp.Av); # Nof full rank terms
    d=length(cp.Av)+length(cp.AvLR);
    rk=size(cp.Z,2);

    # Type logic. The pencil will be of the largest type of the eltypes
    T=promote_type(eltype(cp.Z),
                   eltype(cp.Av[1]),
                   eltype(cp.AvLR[1]),
                   eltype(cp.Bv[1]),
                   eltype(cp.BvLR[1]),
                   eltype(cp.M),
                   eltype(cp.N));

    In=Matrix{T}(I,n,n);
    Idtilde=Matrix{T}(I,rk,rk);

    ## Note: this info is not in CORK paper:
    ##    M11: (dtilde-1) x (dtilde)
    ##    M21: (d-dtilde) x (dtilde)
    ##    M22: (d-dtilde) x (d-dtilde)
    ## same for N

    # Extract Mij and Nij:
    M11=cp.M[1:(dtilde-1),1:dtilde];
    M21=cp.M[(dtilde):end,1:dtilde]
    M22=cp.M[(dtilde):end,(dtilde+1):end]

    N11=cp.N[1:(dtilde-1),1:dtilde];
    N21=cp.N[(dtilde):end,1:dtilde]
    N22=cp.N[(dtilde):end,(dtilde+1):end]

    # Build the B-matrix
    Btilde1=[hcat(cp.Bv[1:dtilde]...) hcat(cp.BvLR...) ]
    Btilde2=[kron(N11,In) zeros((dtilde-1)*n,(d-dtilde)*rk)]
    Btilde3=[kron(N21,cp.Z') kron(N22,Idtilde)];
    Btilde=vcat(Btilde1,Btilde2,Btilde3);

    # Build the A-matrix
    Atilde1=[hcat(cp.Av[1:dtilde]...) hcat(cp.AvLR...) ]
    Atilde2=[kron(M11,In) zeros((dtilde-1)*n,(d-dtilde)*rk)]
    Atilde3=[kron(M21,cp.Z') kron(M22,Idtilde)];
    Atilde=vcat(Atilde1,Atilde2,Atilde3);

    return (Atilde,Btilde);
end

end  # End Module


text:
using LinearAlgebra, Statistics
using BenchmarkTools, Interpolations, Parameters, Plots, QuantEcon, Random, Roots

function coleman_egm(g, k_grid, β, u′, u′_inv, f, f′, shocks)

    # Allocate memory for value of consumption on endogenous grid points
    c = similar(k_grid)

    # Solve for updated consumption value
    for (i, k) in enumerate(k_grid)
        vals = u′.(g.(f(k) * shocks)) .* f′(k) .* shocks
        c[i] = u′_inv(β * mean(vals))
    end

    # Determine endogenous grid
    y = k_grid + c  # y_i = k_i + c_i

    # Update policy function and return
    Kg = LinearInterpolation(y,c, extrapolation_bc=Line())
    Kg_f(x) = Kg(x)
    return Kg_f
end

Model = @with_kw (α = 0.65, # productivity parameter
                  β = 0.95, # discount factor
                  γ = 1.0,  # risk aversion
                  μ = 0.0,  # lognorm(μ, σ)
                  s = 0.1,  # lognorm(μ, σ)
                  grid_min = 1e-6, # smallest grid point
                  grid_max = 4.0,  # largest grid point
                  grid_size = 200, # grid size
                  u = γ == 1 ? log : c->(c^(1-γ)-1)/(1-γ), # utility function
                  u′ = c-> c^(-γ), # u'
                  f = k-> k*α, # production function
                  f′ = k -> α, # f'
                  grid = range(grid_min, grid_max, length = grid_size)) # grid

mlog = Model(); # Log Linear model

Random.seed!(42); # For reproducible behavior.

shock_size = 250     # Number of shock draws in Monte Carlo integral
shocks = exp.(mlog.μ .+ mlog.s * randn(shock_size));

c_star(y) = (1 - mlog.α * mlog.β) * y

# some useful constants
ab = mlog.α * mlog.β
c1 = log(1 - ab) / (1 - mlog.β)
c2 = (mlog.μ + mlog.α * log(ab)) / (1 - mlog.α)
c3 = 1 / (1 - mlog.β)
c4 = 1 / (1 - ab)

v_star(y) = c1 + c2 * (c3 - c4) + c4 * log(y)

n = 15
function check_convergence(m, shocks, c_star, g_init, n_iter)
    k_grid = m.grid
    g = g_init
    plt = plot()
    plot!(plt, m.grid, g.(m.grid),
          color = RGBA(0,0,0,1), lw = 2, alpha = 0.6, label = "initial condition c(y) = y")
    for i in 1:n_iter
        new_g = coleman_egm(g, k_grid, m.β, m.u′, m.u′, m.f, m.f′, shocks)
        g = new_g
        plot!(plt, k_grid, new_g.(k_grid), alpha = 0.6, color = RGBA(0,0,(i / n_iter), 1),
              lw = 2, label = "")
    end

    plot!(plt, k_grid, c_star.(k_grid),
          color = :black, lw = 2, alpha = 0.8, label = "true policy function c*")
    plot!(plt, legend = :topleft)
end

check_convergence(mlog, shocks, c_star, identity, n)

max_stars_repo_name: uq-aibe/aim

max_issues_repo_name: uq-aibe/aim

max_forks_repo_name: uq-aibe/aim


text:
using CodecBGZF
using Test

# Make a buffer type which we can read from after closing it
struct Buffer <: IO
    x::Base.GenericIOBuffer{Vector{UInt8}}
end

Buffer() = Buffer(IOBuffer())
Buffer(x) = Buffer(IOBuffer(x))
Base.read(b::Buffer, ::Type{UInt8}) = read(b.x, UInt8)
Base.write(b::Buffer, x) = write(b.x, x)
Base.write(b::Buffer, x::UInt8) = write(b.x, x)
Base.write(b::Buffer, x::Vector{UInt8}) = write(b.x, x)
Base.take!(b::Buffer) = take!(b.x)
Base.close(b::Buffer) = nothing
Base.isopen(b::Buffer) = isopen(b.x)
Base.bytesavailable(b::Buffer) = bytesavailable(b.x)
Base.eof(b::Buffer) = eof(b.x)
Base.seek(b::Buffer, v::Integer) = seek(b.x, v)

@testset "VirtualOffset" begin
    voff = VirtualOffset(0, 0)
    @test voff == VirtualOffset(0, 0)
    @test voff != VirtualOffset(0, 1)
    @test voff <= VirtualOffset(0, 1)
    @test voff != VirtualOffset(1, 0)
    @test voff <= VirtualOffset(1, 0)
    @test convert(UInt64, VirtualOffset(0, 1)) == UInt(1)
    @test convert(UInt64, VirtualOffset(1, 0)) == UInt(1 << 16)
    @test convert(VirtualOffset, UInt64(0x10001)) == VirtualOffset(1, 1)
    @test isless(voff, VirtualOffset(1, 0))
    @test voff[1] == 0
    @test voff[2] == 0
    @test string(voff) == "VirtualOffset(0, 0)"

    voff += 1
    @test voff[1] == 0
    @test voff[2] == 1

    voff = VirtualOffset(1234, 555)
    @test voff[1] == 1234
    @test voff[2] == 555
    @test string(voff) == "VirtualOffset(1234, 555)"

    voff = VirtualOffset(1234, 555)
    buf = IOBuffer()
    @test write(buf, voff) == sizeof(UInt64)
    seekstart(buf)
    @test read(buf, VirtualOffset) === voff
    @test eof(buf)

    @test_throws ArgumentError VirtualOffset(1 << 48, 0)
    @test_throws ArgumentError VirtualOffset(0, 1 << 16)
end

@testset "Basics" begin
	stream = BGZFDecompressorStream(IOBuffer("foo"))
	@test stream isa BGZFDecompressorStream

	stream = BGZFCompressorStream(IOBuffer("foo"))
	@test stream isa BGZFCompressorStream

    # Bad gzip identifier
    stream = BGZFDecompressorStream(IOBuffer([0x1a, 0x1a]))
	@test_throws BGZFError read(stream, UInt8)

    # Too short
    stream = BGZFDecompressorStream(IOBuffer([0x1f, 0x8b]))
	@test_throws BGZFError read(stream, UInt8)

    # Bad compression flag
    stream = BGZFDecompressorStream(IOBuffer([0x1f, 0x8b, 0x00, 0x04]))
	@test_throws BGZFError read(stream, UInt8)

    # Too short
    stream = BGZFDecompressorStream(IOBuffer([0x1f, 0x8b, 0x08, 0x04]))
	@test_throws BGZFError read(stream, UInt8)

    # Bad flag
    stream = BGZFDecompressorStream(IOBuffer([0x1f, 0x8b, 0x08, 0xfa]))
	@test_throws BGZFError read(stream, UInt8)

	bad_subfield = UInt8[0x1f, 0x8b, 0x08, 0x04, 0x00, 0x00, 0x00, 0x00,
	                     0x00, 0xff, 0x06, 0x00, 0x42, 0x43, 0x03, 0x00,
	                     0x23, 0x00, 0x01, 0x05, 0x00, 0xfa, 0xff, 0x68,
	                     0x65, 0x6c, 0x6c, 0x6f, 0x86, 0xa6, 0x10, 0x36,
	                     0x05, 0x00, 0x00, 0x00]

    stream = BGZFDecompressorStream(IOBuffer(bad_subfield))
	@test_throws BGZFError read(stream, UInt8)

    no_bsize     = UInt8[0x1f, 0x8b, 0x08, 0x04, 0x00, 0x00, 0x00, 0x00,
	                     0x00, 0xff, 0x06, 0x00, 0x41, 0x43, 0x02, 0x00,
	                     0x23, 0x00, 0x01, 0x05, 0x00, 0xfa, 0xff, 0x68,
	                     0x65, 0x6c, 0x6c, 0x6f, 0x86, 0xa6, 0x10, 0x36,
	                     0x05, 0x00, 0x00, 0x00]

    stream = BGZFDecompressorStream(IOBuffer(no_bsize))
	@test_throws BGZFError read(stream, UInt8)

	bad_crc = copy(CodecBGZF.EOF_BLOCK)
	bad_crc[end-6] = 0x01
	@test_throws Exception read(BGZFDecompressorStream(Buffer(bad_crc)), UInt8)

    bad_isize = copy(CodecBGZF.EOF_BLOCK)
	bad_isize[end-2] = 0x01
	@test_throws Exception read(BGZFDecompressorStream(Buffer(bad_isize)), UInt8)
end

@testset "Simple I/O" begin
	# Reading
	filename = joinpath(dirname(@__FILE__), "bar.bgz")
    stream = BGZFDecompressorStream(open(filename))
    @test read(stream, UInt8) === UInt8('b')
    @test read(stream, UInt8) === UInt8('a')
    @test read(stream, UInt8) === UInt8('r')
    @test eof(stream)
    @test read(stream) == UInt8[]
    @test_throws EOFError read(stream, UInt8)
    @test flush(stream) === nothing
    @test close(stream) === nothing
    @test_throws ArgumentError read(stream, UInt8)

	# Writing
	dump = Buffer()
	stream = BGZFCompressorStream(dump)
	write(stream, "ba")
	write(stream, UInt8('r'))
	close(stream)
	data = take!(dump)
	@test data[1:2] == [0x1f, 0x8b]
	@test data[end-27:end] == CodecBGZF.EOF_BLOCK

	# Empty files
	stream = BGZFDecompressorStream(Buffer())
	data = read(stream)
	@test eof(stream)
	@test data == UInt8[]
    close(stream)

    buffer = Buffer()
	stream = BGZFCompressorStream(buffer)
	write(stream, "")
	close(stream)
	@test take!(buffer) == CodecBGZF.EOF_BLOCK

	# No EOF file
	buffer = Buffer()
	stream = BGZFCompressorStream(buffer)
	write(stream, "hello")
	close(stream)
	data = take!(buffer)
	noend = data[1:end-28]
	stream = BGZFDecompressorStream(Buffer(noend))
	@test_throws BGZFError read(stream)
end

@testset "Larger files" begin
    A = collect(reinterpret(UInt8, rand(1:2000, 100000)))
    buffer = Buffer()
    stream = BGZFCompressorStream(buffer)
    write(stream, A)
    close(stream)
    compressed = take!(buffer)

    @test length(compressed) < 350000
    @test compressed[end-27:end] == CodecBGZF.EOF_BLOCK

    stream = BGZFDecompressorStream(Buffer(compressed))
    @test read(stream) == A
end

@testset "More roundtrips" begin
    buffer = Buffer()
    stream = BGZFCompressorStream(buffer)
    write(stream, "")
    write(stream, "")
    write(stream, "Abra")
    write(stream, "")
    write(stream, "cadabra")
    write(stream, "!")
    close(stream)

    stream = BGZFDecompressorStream(Buffer(take!(buffer)))
    @test String(read(stream)) == "Abracadabra!"
end

@testset "Seeking" begin
	filename = joinpath(dirname(@__FILE__), "bar.bgz")
    stream = BGZFDecompressorStream(open(filename))
    read(stream, 1)
    @test read(stream, UInt8) == UInt8('a')
    seekstart(stream)
    @test read(stream, UInt8) == UInt8('b')
end

@testset "BGZF offsets" begin
    function test_virtualoffset(stream, data)
        v = VirtualOffset(stream)
        next = read(stream, 128)
        coff, uoff = offsets(v)
        @test data[coff+1:coff+2] == [0x1f, 0x8b]
        seekstart(stream)
        seek(stream, v)
        @test VirtualOffset(stream) == v
        @test read(stream, 128) == next
    end

    A = collect(reinterpret(UInt8, rand(1:2000, 100000)))
    buffer = Buffer()
    stream = BGZFCompressorStream(buffer)
    write(stream, A)
    close(stream)
    compressed = take!(buffer)

    stream = BGZFDecompressorStream(Buffer(compressed))
    read(stream, 25000)
    test_virtualoffset(stream, compressed)

    # Larger offset
    seekstart(stream)
    read(stream, 125000)
    test_virtualoffset(stream, compressed)
    close(stream)

    # Go back multiple blocks for offset, and also parse lots of empty blocks
    buffer = Buffer()
    for i in 1:20
        stream = BGZFCompressorStream(buffer)
        write(stream, rand(UInt8, 50))
        close(stream)
    end
    data = take!(buffer)
    buffer = Buffer(data)
    stream = BGZFDecompressorStream(buffer)
    read(stream, 10)
    test_virtualoffset(stream, data)
end

@testset "gzi" begin
    bytes = gzi(Buffer(CodecBGZF.EOF_BLOCK))
    @test length(bytes) == 24
    @test collect(reinterpret(Int, bytes)) == [1, 0, 0]

    filename = joinpath(dirname(@__FILE__), "bar.bgz")
    inbytes = open(read, filename)
    bytes = gzi(IOBuffer(inbytes))
    @test length(bytes) == 40
    reallen = [2, 0, 0, length(inbytes) - length(CodecBGZF.EOF_BLOCK), 3]
    @test collect(reinterpret(Int, bytes)) == reallen

    A = collect(reinterpret(UInt8, rand(1:2000, 100000)))
    buffer = Buffer()
    writer = BGZFCompressorStream(buffer)
    write(writer, A)
    close(writer)
    compressed = take!(buffer)
    bytes = gzi(Buffer(compressed))
    buffer = IOBuffer(compressed)
    nums = reinterpret(Int, bytes)
    @assert nums[1] * 2 + 1 == length(nums)
    reader = BGZFDecompressorStream(buffer)

    # Now seek every block offset. It will crash if its not a valid block
    for coffset in nums[2:2:end]
        seek(reader, VirtualOffset(coffset, 0))
        @test VirtualOffset(reader) == VirtualOffset(coffset, 0)
    end
    @test eof(reader)
end


text:
function read_suds(fname::String;
  full::Bool=false,
  v::Int64=KW.v,
  )

  sid = open(fname, "r")
  S = SeisData()

  # Tracking channels and indices
  xc = Array{UnitRange{Int64}, 1}(undef, 32)    # Array of arrays of channel indices
  xn = Array{Int64,1}(undef, 32)                # Array of lengths
  xt = Array{Int64,1}(undef, 32)                # Array of times
  xz = Array{Int64,1}(undef, 32)                # Array of segment sizes
  xnet = Array{UInt16,1}(undef, 32)             # Array of IDs mapped to UInt16
  xi = 1                                        # index into BUF.x
  xj = 0                                        # index into xc, xn
  nx = 0
  fs = 0.0
  fs_last = 0.0

  # File read ================================================================
  # Parse all structures
  cnt = 0
  while !eof(sid)
    cnt += 1
    read_struct_tag!(sid, v)
    # (struct_id, struct_len, nb) = read_struct_tag(sid, v)

    # Attempt to skip unsupported structures
    if SB.sid in unsupported
      nsk = SB.nbs + SB.nbx
      skip(sid, nsk)
      (v > 1) && @warn(string("SUDS struct #", SB.sid, " unsupported; skipped ", nsk, " B"))
      continue
    end

    v > 0 && println("reading structure code ", SB.sid)
    # 6: MUXDATA ------------------------------------------------------------
    if SB.sid == Int16(6)
      # (net, data_type, ts, loctime, Nc, fs, data_type, nz) = read_6!(S, sid, v, full)
      read_6!(S, sid, v, full)
      t0 = 1000000*round(Int64, SB.t_f64)
      (v > 2) && println("nz = ", SB.nz)
      if fs_last == 0.0
        for i in 1:SB.T.n_ch
          if S.fs[i] == 0.0
            S.fs[i] = SB.fs
          end
        end
        fs_last = SB.fs
      elseif fs_last != SB.fs
        @warn(string("fs changes in structure ", cnt))
      end
    else

      # 5: STATIONCOMP  ------------------------------------------------------
      if SB.sid == Int16(5)
        i = read_5!(S, sid, v, full)
        t0 = 1000000*Int64(SB.t_i32)

      # 7: DESCRIPTRACE  -----------------------------------------------------
      elseif SB.sid == Int16(7)
        # (id, data_type, nz, fs, ts) = read_7!(S, sid, v, full)
        read_7!(S, sid, v, full)
        t0 = 1000000*round(Int64, SB.t_f64)

      # 20: COMMENT  ---------------------------------------------------------
      elseif SB.sid == Int16(20)
        read_20!(S, sid, v, full)
        t0 = 0
        continue

      # 30: TIMECORRECTION ---------------------------------------------------
      elseif SB.sid == Int16(30)
        # (tc, rc, t_eff, id, net_flag) = read_30!(S, sid, v, full)
        read_30!(S, sid, v, full)
        SB.data_type = 0x00
        t0 = 0

      # 32: CHANSET ----------------------------------------------------------
      elseif SB.sid == Int16(32)
        read_32!(S, sid, v, full)
        SB.data_type = 0x00
        t0 = 0
        continue

      # 25-29 are skipped unless logged  -------------------------------------
    elseif SB.sid in 25:29 || SB.sid == Int16(31)
        if v > 1
          getfield(SUDS, Symbol(string("read_", SB.sid, "!")))(S, sid, v, full)
          SB.data_type = 0x00
          t0 = 0
        else
          nsk = SB.nbs + SB.nbx
          skip(sid, nsk)
          SB.data_type = 0x00
          t0 = 0
          continue
        end

      # ANYTHING ELSE  -------------------------------------------------------
      else
        getfield(SUDS, Symbol(string("read_", SB.sid, "!")))(S, sid, v, full)
        SB.data_type = 0x00
        t0 = 0
      end
      SB.nz = zero(Int32)
    end

    # Parse data
    if SB.nbx > 0
      # read and reinterpret data
      (v > 2) && println(stdout, "reading ", SB.nbx, " B")
      checkbuf_8!(BUF.buf, SB.nbx)
      readbytes!(sid, BUF.buf, SB.nbx)
      (y, sz) = suds_decode(BUF.buf, SB.data_type)
      nx = div(SB.nbx, sz)

      # Increment xj
      xj += 1
      checkbuf!(BUF.x, xi+nx)
      if xj > length(xn)
        L = length(xc)
        resize!(xc, L+32)
        resize!(xn, L+32)
        resize!(xt, L+32)
        resize!(xz, L+32)
        resize!(xnet, L+32)
      end
      copyto!(BUF.x, xi, y, 1, nx)

      # Store channel indices to xc
      if SB.sid == Int16(6)
        xc[xj] = 1:SB.T.n_ch
        xnet[xj] = SB.T.net
      elseif SB.sid == Int16(7)
        j =  findid(SB.id_str, S)
        xc[xj] = j:j
        xnet[xj] = 0x0000
      end

      # Store start indices to xn and per-segment sizes to xz. increment xi
      xn[xj] = xi
      xt[xj] = t0
      xz[xj] = SB.nz == zero(Int32) ? nx : SB.nz
      xi += nx
    end

    # if structure is a time correction, flush the suds buffers and write to S
    if SB.sid == Int16(30)

      # First identify which start times get corrected
      net = 0x0000
      tc = round(Int64, SB.tc*1.0e6)
      if SB.irig
        net = reinterpret(UInt16, SB.id[1:2])[1]
        for j = 1:xj
          if xnet[j] == net
            xt[j] += tc
          end
        end
      end
      flush_suds!(S, xc, xn, xt, xz, xj, v)

      # Adjust fs after flushing buffers, else fictitious time gaps appear
      if SB.irig
        net_str = String(copy(SB.id[1:2]))
        if SB.rc != 0.0f0
          for i = 1:S.n
            if (S.id[i][1:2] == net_str) && (S.fs[i] != SB.fs + SB.rc)
              (v > 2) && println("Adjusting S.fs[", i, "]")
              S.fs[i] += SB.rc
            end
          end
        end
      end
      xi = 1
      xj = 0
      SB.irig = false

    end
  end
  flush_suds!(S, xc, xn, xt, xz, xj, v)
  resize!(BUF.buf, 65535)
  resize!(BUF.x, 65535)

  return S
end

read_suds!(
  S::GphysData,
  fname::String;
  full::Bool=false,
  v::Int64=KW.v,
  ) = (U = read_suds(fname, full=full, v=v); append!(S,U))

function readsudsevt(fname::String;
  full::Bool=false,
  v::Int64=KW.v,
  )

  TD = read_suds(fname, full=full, v=v)
  src_auth = get(auth, Int32(SB.H.auth), "")
  m = SUDS.SB.H.shorts[6]
  mag_scale = try
    mag_scale[m]
  catch
    replace(join(Char.(reinterpret(UInt8, [m]))), "\0" => "")
  end

  # Generate empty event structs
  Mag = EQMag(val   = SB.H.mag[1],
              nst   = Int64(SB.H.shorts[7]),
              scale = mag_scale,
              src   = src_auth
              )

  # Create loc
  Loc = EQLoc(lat   = SB.H.lat,
              lon   = SB.H.lon,
              dep   = Float64(SB.H.floats[1]),
              dx    = Float64(SB.H.floats[2]),
              dy    = Float64(SB.H.floats[2]),
              dz    = Float64(SB.H.floats[3]),
              rms   = Float64(SB.H.floats[4]),
              gap   = Float64(SB.H.gap),
              dmin  = Float64(SB.H.d_min),
              nst   = Int64(SB.H.shorts[1]),
              flags = SB.H.chars[6] in (0x47, 0x4e, 0x53, 0x65, 0x66) ? 0x20 : 0x00,
              src   = get(loc_prog, Char(SB.H.chars[6]), "Unknown location program")
              )

  # Create header
  H = SeisHdr(id    = string(SB.H.evno),
              loc   = Loc,
              mag   = Mag,
              ot    = u2d(SB.H.ot),
              src   = fname
              )
  H.misc["auth"]  = src_auth
  H.misc["reg"]   = SB.H.reg
  H.misc["model"] = String(copy(SB.H.model))

  # Create event container
  Ev = SeisEvent(hdr = H, data = TD)

  return Ev
end


text:
export EmptyGridWorld

mutable struct EmptyGridWorld <: AbstractGridWorld
    world::GridWorldBase{Tuple{Empty,Wall,Goal}}
    agent_pos::CartesianIndex{2}
    agent::Agent
    goal_reward::Float64
    reward::Float64
end

function EmptyGridWorld(;n=8, agent_start_pos=CartesianIndex(2,2), agent_start_dir=RIGHT, goal_pos = CartesianIndex(n-1, n-1))
    objects = (EMPTY, WALL, GOAL)
    w = GridWorldBase(objects, n, n)

    w[WALL, [1,n], 1:n] .= true
    w[WALL, 1:n, [1,n]] .= true

    goal_reward = 1.0
    reward = 0.0

    env = EmptyGridWorld(w, agent_start_pos, Agent(dir=agent_start_dir), goal_reward, reward)

    reset!(env, agent_start_pos = agent_start_pos, agent_start_dir = agent_start_dir, goal_pos = goal_pos)

    return env
end

function (env::EmptyGridWorld)(::MoveForward)
    dir = get_dir(env.agent)
    dest = dir(env.agent_pos)
    env.reward = 0.0
    if !env.world[WALL, dest]
        env.agent_pos = dest
        if env.world[GOAL, env.agent_pos]
            env.reward = env.goal_reward
        end
    end
    env
end

RLBase.get_terminal(env::EmptyGridWorld) = env.world[GOAL, env.agent_pos]

function RLBase.reset!(env::EmptyGridWorld; agent_start_pos = CartesianIndex(2, 2), agent_start_dir = RIGHT, goal_pos = CartesianIndex(size(env.world)[end] - 1, size(env.world)[end] - 1))

    env.reward = 0.0
    env.agent_pos = agent_start_pos
    agent = get_agent(env)
    set_dir!(agent, agent_start_dir)

    env.world[EMPTY, :, :] .= .!env.world[WALL, :, :]
    env.world[GOAL, goal_pos] = true
    env.world[EMPTY, goal_pos] = false

    return env
end


text:
module AbstractIndices

using LinearAlgebra, Statistics, ArrayInterface, StaticRanges

using ArrayInterface: can_setindex

using Base: to_index, axes, broadcasted, AbstractCartesianIndex, @_propagate_inbounds_meta

import Base: getindex, setindex, iterate, checkindex, checkbounds

using Base.Broadcast:
    Broadcasted, BroadcastStyle, DefaultArrayStyle, AbstractArrayStyle, Unknown, combine_axes

using StaticRanges:
    can_set_length,
    can_set_first,
    can_set_last,
    ForwardOrdering,
    ReverseOrdering,
    ContinuousTrait,
    Continuous,
    DiscreteTrait,
    Discrete,
    similar_type,
    OneToRange,
    AbstractStepRange,
    AbstractStepRangeLen,
    AbstractLinRange,
    StaticUnitRange,
    Length,
    Size

import Base.Broadcast: combine_axes

# iterators
import Base: iterate, isdone, has_offset_axes
import Base.Iterators: Pairs

import Base: OneTo, tail, show, values, keys, @propagate_inbounds

export AbstractIndex,
       Index,
       SimpleIndex,
       IndicesArray,
       IndicesMatrix,
       IndicesVector,
       IArray,
       IMatrix,
       IVector,
       # methods
       dimnames,
       unname,
       # reexports
       mrange,
       srange,
       OneToSRange,
       OneToMRange,
       UnitMRange,
       UnitSRange,
       pop,
       popfirst

const TupOrVec{T} = Union{Tuple{Vararg{T}},AbstractVector{T}}

include("IndexCore/IndexCore.jl")
using .IndexCore

include("IndicesArrays/IndicesArrays.jl")
include("LinearAlgebra/LinearAlgebra.jl")

end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

name = "Qt5Base"
version = v"5.15.2"

# Collection of sources required to build qt5
sources = [
    ArchiveSource("https://download.qt.io/official_releases/qt/$(version.major).$(version.minor)/$version/submodules/qtbase-everywhere-src-$version.tar.xz",
                  "909fad2591ee367993a75d7e2ea50ad4db332f05e1c38dd7a5a274e156a4e0f8"),
    ArchiveSource("https://github.com/phracker/MacOSX-SDKs/releases/download/10.15/MacOSX10.14.sdk.tar.xz",
                  "0f03869f72df8705b832910517b47dd5b79eb4e160512602f593ed243b28715f")
]

script = raw"""
cd $WORKSPACE/srcdir

BIN_DIR="/opt/bin/${bb_full_target}"

mkdir build
cd build/

qtsrcdir=`ls -d ../qtbase-everywhere-src-*`

commonoptions=" \
-opensource -confirm-license \
-openssl-linked  -nomake examples -release \
"

export OPENSSL_LIBS="-L${libdir} -lssl -lcrypto"

cp -a ../qtbase-everywhere-src-*/mkspecs/linux-aarch64-gnu-g++ $qtsrcdir/mkspecs/musl
sed -i 's/aarch64-linux-gnu/x86_64-linux-musl/g' ../qtbase-everywhere-src-*/mkspecs/musl/qmake.conf

cp -a ../qtbase-everywhere-src-*/mkspecs/linux-aarch64-gnu-g++ $qtsrcdir/mkspecs/$target
sed -i "s/aarch64-linux-gnu/$target/g" ../qtbase-everywhere-src-*/mkspecs/$target/qmake.conf
echo "QMAKE_CFLAGS_ARCH_HASWELL =" >> ../qtbase-everywhere-src-*/mkspecs/$target/qmake.conf

case "$target" in

	*linux*)
        ../qtbase-everywhere-src-*/configure -L ${libdir} -I ${includedir} \
            -platform musl -xplatform $target -device-option CROSS_COMPILE=${BIN_DIR}/$target- \
            -prefix $prefix $commonoptions \
            -fontconfig
        ;;

	*apple-darwin*)
        echo "QMAKE_CFLAGS_ARCH_HASWELL =" >> ../qtbase-everywhere-src-*/mkspecs/macx-clang/qmake.conf
        cd $WORKSPACE/srcdir/MacOSX10.14.sdk
        rm -rf /opt/$target/$target/sys-root/System
        rsync -a usr/* /opt/$target/$target/sys-root/usr/
        cp -a System /opt/$target/$target/sys-root/

        cd $WORKSPACE/srcdir/qtbase-everywhere-src-*/

        cat <<EOT > mkspecs/features/mac/default_pre.prf
CONFIG = asset_catalogs rez \$\$CONFIG
load(default_pre)

QMAKE_ASSET_CATALOGS_APP_ICON = AppIcon

macx-xcode:qtConfig(static): \\
    QMAKE_XCODE_DEBUG_INFORMATION_FORMAT = dwarf

QMAKE_XCODE_LIBRARY_SUFFIX_SETTING = QT_LIBRARY_SUFFIX

xcode_copy_phase_strip_setting.name = COPY_PHASE_STRIP
xcode_copy_phase_strip_setting.value = NO
QMAKE_MAC_XCODE_SETTINGS += xcode_copy_phase_strip_setting
EOT

        sed -i '1s;^;QMAKE_MAC_SDK.macosx.Path = '"/opt/$target/$target/sys-root"'\
        QMAKE_MAC_SDK.macosx.SDKVersion = '"10.14"'\
        QMAKE_MAC_SDK.macosx.PlatformPath = '"/opt/$target"'\n;' 'mkspecs/features/mac/sdk.prf'
        echo "" >  mkspecs/features/mac/no_warn_empty_obj_files.prf

        sed -i "s?-fuse-ld=x86_64-apple-darwin14?-fuse-ld=${BIN_DIR}/x86_64-apple-darwin14-ld?g" ${BIN_DIR}/x86_64-apple-darwin14-clang++
        sed -i "s?-fuse-ld=x86_64-apple-darwin14?-fuse-ld=${BIN_DIR}/x86_64-apple-darwin14-ld?g" ${BIN_DIR}/x86_64-apple-darwin14-clang

        cd $WORKSPACE/srcdir/build

        export QT_MAC_SDK_NO_VERSION_CHECK=1
        ../qtbase-everywhere-src-*/configure \
            QMAKE_CXXFLAGS+=-F/opt/$target/$target/sys-root/System/Library/Frameworks \
            QMAKE_RANLIB=${BIN_DIR}/ranlib \
            QMAKE_MACOSX_DEPLOYMENT_TARGET=10.14 \
            -platform musl -xplatform macx-clang -device-option CROSS_COMPILE=${BIN_DIR}/$target- \
            -prefix ${prefix} $commonoptions \
            -skip qtwinextras
        ;;

    *mingw*)
        echo "QMAKE_CFLAGS_ARCH_HASWELL =" >> ../qtbase-everywhere-src-*/mkspecs/win32-g++/qmake.conf
        ../qtbase-everywhere-src-*/configure -I $WORKSPACE/srcdir/qtbase-everywhere-src-*/include/QtANGLE -platform musl -xplatform win32-g++ -device-option CROSS_COMPILE=${BIN_DIR}/$target- \
            -prefix $prefix $commonoptions \
            -opengl dynamic
		;;

    *x86_64-unknown-freebsd*)
        sed -i 's/load(qt_config)//' ../qtbase-everywhere-src-*/mkspecs/freebsd-g++/qmake.conf
        grep -A11 QMAKE_CC ../qtbase-everywhere-src-*/mkspecs/linux-aarch64-gnu-g++/qmake.conf | sed -e 's/aarch64-linux-gnu/x86_64-unknown-freebsd11.1/' >> ../qtbase-everywhere-src-*/mkspecs/freebsd-g++/qmake.conf
        echo "QMAKE_CFLAGS_ARCH_HASWELL =" >> ../qtbase-everywhere-src-*/mkspecs/freebsd-g++/qmake.conf

        ../qtbase-everywhere-src-*/configure -platform musl -xplatform freebsd-g++ -device-option CROSS_COMPILE=${BIN_DIR}/$target- \
            -extprefix $prefix $commonoptions \
            -skip qtwinextras -fontconfig -sysroot /opt/$target/bin/../$target/sys-root
		;;
esac

make -j${nproc}
make install
$BUILD_STRIP $bindir/qmake
install_license $WORKSPACE/srcdir/qtbase-everywhere-src-*/LICENSE.LGPLv3
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms_linux = [
    Platform("aarch64", "linux"; libc="glibc"),
    Platform("armv7l", "linux"; libc="glibc"),
    Platform("x86_64", "linux"; libc="glibc"),
    Platform("i686", "linux"; libc="glibc"),
    Platform("x86_64", "freebsd"),
    Platform("powerpc64le", "linux"; libc="glibc"),
]
platforms_linux = expand_cxxstring_abis(platforms_linux)
platforms_win = expand_cxxstring_abis([Platform("x86_64", "windows"), Platform("i686", "windows")])
platforms_macos = [ Platform("x86_64", "macos") ]

# The products that we will ensure are always built
products = [
    LibraryProduct(["Qt5Concurrent", "libQt5Concurrent", "QtConcurrent"], :libqt5concurrent),
    LibraryProduct(["Qt5Core", "libQt5Core", "QtCore"], :libqt5core),
    LibraryProduct(["Qt5DBus", "libQt5DBus", "QtDBus"], :libqt5dbus),
    LibraryProduct(["Qt5Gui", "libQt5Gui", "QtGui"], :libqt5gui),
    LibraryProduct(["Qt5Network", "libQt5Network", "QtNetwork"], :libqt5network),
    LibraryProduct(["Qt5OpenGL", "libQt5OpenGL", "QtOpenGL"], :libqt5opengl),
    LibraryProduct(["Qt5PrintSupport", "libQt5PrintSupport", "QtPrintSupport"], :libqt5printsupport),
    LibraryProduct(["Qt5Sql", "libQt5Sql", "QtSql"], :libqt5sql),
    LibraryProduct(["Qt5Test", "libQt5Test", "QtTest"], :libqt5test),
    LibraryProduct(["Qt5Widgets", "libQt5Widgets", "QtWidgets"], :libqt5widgets),
    LibraryProduct(["Qt5Xml", "libQt5Xml", "QtXml"], :libqt5xml),
]

products_macos = [
    FrameworkProduct("QtConcurrent", :libqt5concurrent),
    FrameworkProduct("QtCore", :libqt5core),
    FrameworkProduct("QtDBus", :libqt5dbus),
    FrameworkProduct("QtGui", :libqt5gui),
    FrameworkProduct("QtNetwork", :libqt5network),
    FrameworkProduct("QtOpenGL", :libqt5opengl),
    FrameworkProduct("QtPrintSupport", :libqt5printsupport),
    FrameworkProduct("QtSql", :libqt5sql),
    FrameworkProduct("QtTest", :libqt5test),
    FrameworkProduct("QtWidgets", :libqt5widgets),
    FrameworkProduct("QtXml", :libqt5xml),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    BuildDependency("Xorg_libX11_jll"),
    Dependency("Xorg_libXext_jll"),
    BuildDependency("Xorg_glproto_jll"),
    Dependency("Xorg_libxcb_jll"),
    Dependency("Xorg_xcb_util_wm_jll"),
    Dependency("Xorg_xcb_util_image_jll"),
    Dependency("Xorg_xcb_util_keysyms_jll"),
    Dependency("Xorg_xcb_util_renderutil_jll"),
    Dependency("xkbcommon_jll"),
    Dependency("Libglvnd_jll"),
    Dependency("Fontconfig_jll"),
    Dependency("Glib_jll"),
    Dependency("Zlib_jll"),
    Dependency("CompilerSupportLibraries_jll"),
    Dependency("OpenSSL_jll"),
]

include("../../fancy_toys.jl")

if any(should_build_platform.(triplet.(platforms_linux)))
    build_tarballs(ARGS, name, version, sources, script, platforms_linux, products, dependencies; preferred_gcc_version = v"7")
end
if any(should_build_platform.(triplet.(platforms_win)))
    build_tarballs(ARGS, name, version, sources, script, platforms_win, products, dependencies; preferred_gcc_version = v"8")
end
if any(should_build_platform.(triplet.(platforms_macos)))
    build_tarballs(ARGS, name, version, sources, script, platforms_macos, products_macos, dependencies; preferred_gcc_version = v"7")
end

max_stars_repo_name: waralex/Yggdrasil

max_issues_repo_name: waralex/Yggdrasil

max_forks_repo_name: waralex/Yggdrasil


text:
half(L) = div(L - 1, 2)

logZ(M::Matrix) = logdet(I + M)
Z(M::Matrix) = det(I + M)

function pairing(L, h)
    E(k) = sqrt(1 + h^2 + 2h * cospi((2k + 1) / L))
    ζ(k) = - (h + cospi((2k + 1) / L))
    u(k) = sqrt((1 + ζ(k) / E(k)) / 2)
    v(k) = sqrt((1 - ζ(k) / E(k)) / 2)
    g(r) = 2 / L * sum(sinpi((2k + 1) / L * r) * v(k) / u(k) for k in 0:div(L - 1, 2))

    M = zeros(L, L)
    for i in 1:L, j in 1:L
        M[i,j] = g(i -  j)
    end
    return M
end

function get_prob(x::BitVector, M::Matrix)
    amp = log(abs.(det(M[x .> 0, x .> 0])))
    return amp
end

function get_all_probs(M::Matrix)
    L = size(M, 1)
    logZ = log(abs(det(M)))
    x = i->bitarray(i, L)
    amps = [get_prob(x(i), M) for i in 0:2^L - 1]
    amps = exp.(amps)
    amps2 = normalize(amps, 1)
    return round.(amps, digits = 7), round.(amps2, digits = 7)
end

# for L in [2, 4,8,16,20, 32]
#     P = pairing(L, 1.0)
#     Z1 = logdet(I + P)
#     Z2 = log(sum(get_all_probs(P)[1]))
#     println("L is $L Z is $Z1 and $Z2")
# end


text:
using .Networks

using JuMP
using Dates
using DataStructures
using Printf
using Parameters

"""
REF_SCHEDULE_TYPE : Indicates wether to consider the preceding market or TSO schedule as a reference.
                    The reference schedule is used to get decided commitment and production levels if
                      tso actions are missing.
"""
@with_kw mutable struct TSOConfigs
    CONSIDER_N_1_CSTRS::Bool = false
    loss_of_load_penalty = 1e7
    limitation_penalty = 1e-03
    out_path = nothing
    problem_name = "TSO"
    REF_SCHEDULE_TYPE::Union{Market,TSO} = TSO();
end

@with_kw struct TSOLimitableModel <: AbstractLimitableModel
    #gen,ts,s
    delta_p = SortedDict{Tuple{String,DateTime,String},VariableRef}();
    #gen,ts,s
    p_injected = SortedDict{Tuple{String,DateTime,String},VariableRef}();
    #gen,ts,s
    p_limit = SortedDict{Tuple{String,DateTime,String},VariableRef}();
    #gen,ts,s
    b_is_limited = SortedDict{Tuple{String,DateTime,String},VariableRef}();
    #gen,ts,s
    p_limit_x_is_limited = SortedDict{Tuple{String,DateTime,String},VariableRef}();
end

@with_kw struct TSOPilotableModel <: AbstractPilotableModel
    #gen,ts,s
    delta_p = SortedDict{Tuple{String,DateTime,String},VariableRef}();
    #gen,ts,s
    p_injected = SortedDict{Tuple{String,DateTime,String},VariableRef}();
    #gen,ts,s
    b_start = SortedDict{Tuple{String,DateTime,String},VariableRef}();
    #gen,ts,s
    b_on = SortedDict{Tuple{String,DateTime,String},VariableRef}();
end

@with_kw struct TSOLoLModel <: AbstractLoLModel
    #bus,ts,s
    p_loss_of_load = SortedDict{Tuple{String,DateTime,String},VariableRef}();
end

@with_kw mutable struct TSOObjectiveModel <: AbstractObjectiveModel
    deltas = AffExpr(0)

    prop_cost = AffExpr(0)
    start_cost = AffExpr(0)

    penalty = AffExpr(0)

    full_obj_1 = AffExpr(0)
    full_obj_2 = AffExpr(0)
end

@with_kw mutable struct TSOModel <: AbstractModelContainer
    model::Model = Model()
    limitable_model::TSOLimitableModel = TSOLimitableModel()
    pilotable_model::TSOPilotableModel = TSOPilotableModel()
    lol_model::TSOLoLModel = TSOLoLModel()
    objective_model::TSOObjectiveModel = TSOObjectiveModel()
    #ts,s
    eod_constraint::SortedDict{Tuple{Dates.DateTime,String}, ConstraintRef} =
        SortedDict{Tuple{Dates.DateTime,String}, ConstraintRef}()
    #branch,ts,s,ptdf_case
    flows::SortedDict{Tuple{String,DateTime,String,String},AffExpr} =
        SortedDict{Tuple{String,DateTime,String,String},AffExpr}()
    rso_constraint::SortedDict{Tuple{String,DateTime,String,String},ConstraintRef} =
        SortedDict{Tuple{String,DateTime,String,String},ConstraintRef}()
end

function has_positive_slack(model_container::TSOModel)::Bool
    return has_positive_value(model_container.lol_model.p_loss_of_load)
end

function sum_capping(limitable_model::TSOLimitableModel, ts,s, network::Networks.Network)
    error("TODO : requires uncertainties cause capping=uncertainties-injection")
end

#TODO define a struct LocalisedLolModel to use it for TSOLoLModel and TSOBilevelTSOLoLModel
function sum_lol(lol_model::TSOLoLModel, ts, s, network::Networks.Network)
    sum_l = 0.
    for bus in Networks.get_buses(network)
        bus_id = Networks.get_id(bus)
        sum_l += lol_model.p_loss_of_load[bus_id,ts,s]
    end
    return sum_l
end

function add_tso_limitable_prop_cost!(obj_component::AffExpr,
                                uncertainties_at_ech::UncertaintiesAtEch,
                                p_injected::AbstractDict{T,V}, network)  where T <: Tuple where V <: VariableRef
    #NOTE: need to make sure uncertainty > injection

    for ((gen_id,ts,s), p_injected_var) in p_injected
        generator = Networks.get_generator(network, gen_id)
        gen_prop_cost = Networks.get_prop_cost(generator)
        uncertainty = get_uncertainties(uncertainties_at_ech, gen_id, ts, s)
        add_to_expression!(obj_component,
                            (uncertainty - p_injected_var) * gen_prop_cost)
    end

    return obj_component
end

function create_objectives!(model_container::TSOModel,
                            network, uncertainties_at_ech, gratis_starts, loss_of_load_cost, limitation_penalty)

    # cost for cutting load/consumption
    add_coeffxsum_cost!(model_container.objective_model.penalty,
                        model_container.lol_model.p_loss_of_load, loss_of_load_cost)

    # avoid limiting when not necessary
    for (_, var_is_limited) in model_container.limitable_model.b_is_limited
        model_container.objective_model.penalty += limitation_penalty * var_is_limited
    end

    ## Objective 1 :

    # cost for deviating from market schedule
    for (_, var_delta) in model_container.pilotable_model.delta_p
        model_container.objective_model.deltas += var_delta
    end
    for (_, var_delta) in model_container.limitable_model.delta_p
        model_container.objective_model.deltas += var_delta
    end

    ## Objective 2 :

    # cost for starting pilotables
    add_pilotable_start_cost!(model_container.objective_model.start_cost,
                            model_container.pilotable_model.b_start, network, gratis_starts)

    # cost for limitables : cost of capped limitable power
    add_tso_limitable_prop_cost!(model_container.objective_model.prop_cost,
                                uncertainties_at_ech,
                                model_container.limitable_model.p_injected, network)

    # cost for using pilotables
    add_prop_cost!(model_container.objective_model.prop_cost,
                            model_container.pilotable_model.p_injected, network)

    # Objective 1 :
    model_container.objective_model.full_obj_1 = ( model_container.objective_model.deltas +
                                                model_container.objective_model.penalty )
    # Objective 2 :
    model_container.objective_model.full_obj_2 = ( model_container.objective_model.start_cost +
                                                model_container.objective_model.prop_cost +
                                                model_container.objective_model.penalty )

    return model_container
end

function tso_out_fo(network::Networks.Network,
                    target_timepoints::Vector{Dates.DateTime},
                    generators_initial_state::SortedDict{String,GeneratorState},
                    scenarios::Vector{String},
                    uncertainties_at_ech::UncertaintiesAtEch,
                    firmness::Firmness,
                    preceding_market_schedule::Schedule,
                    preceding_tso_schedule::Schedule,
                    gratis_starts::Set{Tuple{String,Dates.DateTime}},
                    configs::TSOConfigs
                    )

    # TODO : check coherence between : preceding_reference_schedule and TSOActions.impositions cause we do not consider TSOActions
    if is_market(configs.REF_SCHEDULE_TYPE)
        reference_schedule = preceding_market_schedule
    elseif is_tso(configs.REF_SCHEDULE_TYPE)
        reference_schedule = preceding_tso_schedule
    else
        throw( error("Invalid REF_SCHEDULE_TYPE config.") )
    end

    pilotables_list_l = Networks.get_generators_of_type(network, Networks.PILOTABLE)
    limitables_list_l = Networks.get_generators_of_type(network, Networks.LIMITABLE)
    buses_list_l = Networks.get_buses(network)

    model_container_l = TSOModel()

    # Variables
    add_pilotables_vars!(model_container_l,
                            pilotables_list_l, target_timepoints, scenarios,
                            preceding_market_schedule,
                            injection_vars=true, commitment_vars=true, delta_vars=true)
    add_limitables_vars!(model_container_l, target_timepoints, scenarios,
                            limitables_list_l, preceding_market_schedule,
                            injection_vars=true, limit_vars=true, delta_vars=true
                            )
    add_lol_vars!(model_container_l, target_timepoints, scenarios,
                    buses_list_l,
                    local_lol_vars=true)

    # Constraints

    # Pilotables
    pilotable_power_constraints!(model_container_l.model,
                                model_container_l.pilotable_model, pilotables_list_l, target_timepoints, scenarios,
                                firmness, reference_schedule,
                                always_link_scenarios=false)
    unit_commitment_constraints!(model_container_l.model,
                                model_container_l.pilotable_model, pilotables_list_l,  target_timepoints, scenarios,
                                firmness, reference_schedule, generators_initial_state,
                                always_link_scenarios=false)

    # Limitables
    limitable_power_constraints!(model_container_l.model,
                                model_container_l.limitable_model, limitables_list_l, target_timepoints, scenarios,
                                firmness, uncertainties_at_ech,
                                always_link_scenarios=false)

    # LoL
    local_lol_constraints!(model_container_l.model,
                            model_container_l.lol_model, buses_list_l, target_timepoints, scenarios,
                            uncertainties_at_ech)

    # EOD
    eod_constraints!(model_container_l.model, model_container_l.eod_constraint,
                    model_container_l.pilotable_model,
                    model_container_l.limitable_model,
                    model_container_l.lol_model,
                    target_timepoints, scenarios,
                    uncertainties_at_ech, network
                    )

    # RSO
    add_rso_flows_exprs(model_container_l.flows,
                model_container_l.pilotable_model,
                model_container_l.limitable_model,
                model_container_l.lol_model,
                all_combinations(network, target_timepoints, scenarios),
                uncertainties_at_ech, network)
    combinations = (configs.CONSIDER_N_1_CSTRS) ? all_combinations(network, target_timepoints, scenarios) :
                                                    all_n_combinations(network, target_timepoints, scenarios)
    rso_constraints!(model_container_l.model,
                    model_container_l.rso_constraint,
                    model_container_l.flows,
                    combinations,
                    network,
                    prefix="tso")

    create_objectives!(model_container_l,
                        network, uncertainties_at_ech,
                        gratis_starts,
                        configs.loss_of_load_penalty, configs.limitation_penalty)

    solve_2steps_deltas!(model_container_l, configs)

    return model_container_l
end

max_stars_repo_name: klorel/pscopf

max_issues_repo_name: klorel/pscopf

max_forks_repo_name: klorel/pscopf


text:
@testset "gaussian.jl" begin
    model = Gaussian_Peak()

    @test parameter_size(model) == 3

    θ=Float64[1, 0, 1]

    @test eval_y(model,[0.0],θ) ≈ [1.0]

    X=rand(10)
    Y=alloc_y(model,X,θ)

    @test (@benchmark accumulate_y!($model,$Y,$X,$θ)).allocs == 0

end


text:
function readNewick(args...; filepath::String = "", format::String = "nwk")
    # read in the string of tree file
    instream = open(expanduser(filepath));
    instring = readstring(instream);
    close(instream);
    trees = split(instring,"");  # split the each element of the tree string
    trees = trees[trees.!="\r"];
    trees = trees[trees.!="\n"];

    # calculate the sum of "(" and "," in the tree in order to set up the edge dimension
      count1 = 0
      count2 = 0
      for i = 1:length(trees)
        if trees[i] == "("
            count1 = count1 + 1
        end
        if trees[i] == ","
            count2 = count2 + 1
        end
      end

    edge = Array{Float64}(count1 + count2, 2);
    node_label = Array{String}(count1, 1); # generate empty node_label vector
    tip_label = String[];  # generate empty tip_label vector

    #ei = Float64[]
    edge_length = Array{Float64}(count1 + count2, 1);
    ei = Array{Int64}(count1+1, 1);
   # edge_length = Array(Float64,1)

    currnode = 1;
    Nnode = currnode;
    i = j = k = 1;

#   while trees[i]!=";"
    for i = 1:(size(trees)[1] - 1)
        if trees[i] == "("
            edge[j, 1] = currnode
            i = i + 1

            # is the next element a label?
            if countnz(trees[i] .== ["("; ")"; ","; ":"; ";"]) == 0
                temp = getLabel(Tree = trees, Start = i, stop_char = [","; ":"; ")"; ";"])
                #btip_label[k] = temp["label"]
                push!(tip_label, temp["label"])
                i = temp["end"]
                edge[j,2] = -k
                k = k + 1
                # is there a branch length?
                if trees[i] == ":"
                    temp = getEdgeLength(Tree = trees, Start = i)
                    edge_length[j] = temp["edge_length"]
                    #push!(edge_length,temp["edge_length"])
                    i = temp["end"]
                end
               # push!(ei,NaN)

            elseif trees[i] == "("
                Nnode = Nnode + 1
                currnode = Nnode
                edge[j, 2] = currnode
                #push!(ei,currnode)
                ei[currnode] = j
            end
            j = j + 1

        elseif trees[i] == ")"
                i = i + 1
                # is the next element a label?
                if countnz(trees[i] .== ["("; ")"; ","; ":"; ";"]) == 0
                    temp = getLabel(Tree = trees, Start = i, stop_char = [","; ":"; ")"; ";"])
                    node_label[currnode] = temp["label"]
                    i = temp["end"]
                end
                ii = Int(ei[currnode])

                if trees[i] == ":"
                    temp = getEdgeLength(Tree = trees, Start = i)
                    if currnode >1
                         edge_length[ii] = temp["edge_length"]
                        #push!(edge_length,temp["edge_length"])
                    else
                        root_edge = temp["edge_length"]
                        i = temp["end"]
                    end
                end
                if currnode > 1
                    currnode = Int(edge[ii, 1])
                end

        elseif trees[i] == ","
            edge[j, 1] = currnode
            i = i + 1
            if countnz(trees[i] .== ["("; ")"; ","; ":"; ";"]) == 0
                temp = getLabel(Tree = trees, Start = i, stop_char = [","; ":"; ")"; ";"])
                push!(tip_label, temp["label"])
                i = temp["end"]
                edge[j, 2] = -k
                k = k + 1
                # is there a branch length?
                if trees[i] == ":"
                    temp = getEdgeLength(Tree = trees, Start = i)
                    edge_length[j] = temp["edge_length"]
                   # push!(edge_length,temp["edge_length"])

                    i = temp["end"]
                end

            elseif trees[i] == "("
                Nnode = Nnode + 1
                currnode = Nnode
                edge[j, 2] = currnode
                #push!(ei,j)
                 ei[currnode] = j
            end
            j = j + 1
        end

    end

    Ntip = k - 1
    temInd = edge[:,:] .> 0
    edge[temInd[:, 2], 2] = edge[temInd[:, 2], 2] .+ Ntip
    edge[temInd[:, 1], 1] = edge[temInd[:, 1], 1] .+ Ntip
    edge = abs.(edge)

    tree = ((edge), (Nnode), (tip_label), (edge_length), (node_label))
    tree  = Dict(zip(["edge", "Nnode", "tip_label", "edge_length", "node_label"], tree))
    return tree
end


text:
"""
    readvariablelength(f::IO)
Variable length numbers in MIDI files are represented as a sequence of bytes.
If the first bit is 0, we're looking at the last byte in the sequence. The remaining
7 bits indicate the number.
"""
function readvariablelength(f::IO)

    mask = 0b10000000
    notmask = ~mask
    # Read the first byte
    b = read(f, UInt8)
    bytes = UInt8[]
    if (b & mask) == 0
        # We're done here. The first bit isn't set, so the number is contained in the 7 remaining bits.
        convert(Int, b)
    else
        result = convert(UInt32, 0)
        while (b & mask) == mask
            result <<= 7
            result += (b & notmask)
            b = read(f, UInt8)
        end
        result = (result << 7) + b # No need to "& notmask", since the most significant bit is 0
        Int(reinterpret(Int32,result))
    end
end

"""
    writevariablelength(f::IO, number::Int)
Write on `f` the given `number`, firstly converting it to the "variable length" format.
See the documentation for more.
"""
function writevariablelength(f::IO, number::Int)
    if number < 0 || number > 0x0FFFFFFF
        error("Unable to write variable length value ",
              string(number),
              " as it does not fit in the range [0, 0x0FFFFFFF].")
    end
    if number < 128
        write(f, UInt8(number))
    else
        bytes = UInt8[]

        push!(bytes, UInt8(number & 0x7F)) # Get the bottom 7 bits
        number >>>= 7 # Is there a bug with Julia here? Testing in the REPL on negative numbers give >> and >>> the same result
        while number > 0
            push!(bytes, UInt8(((number & 0x7F) | 0x80)))
            number >>>= 7
            continuation = 0x80
        end
        reverse!(bytes)
        for b in bytes
            write(f, b)
        end
    end
end


text:
function OrekitEphemerisFile()
    return OrekitEphemerisFile(())
end

function add_satellite(obj::OrekitEphemerisFile, arg0::JString)
    return jcall(obj, "addSatellite", OrekitEphemerisFile_OrekitSatelliteEphemeris, (JString,), arg0)
end

function get_satellites(obj::OrekitEphemerisFile)
    return jcall(obj, "getSatellites", Map, ())
end


text:
@testset "Graphs" begin
    using LightGraphs, MetaGraphs

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

    cc(ins...; name) = conc(ins...; dims=2, traitdecoration=tf(name))
    nc(name) = traitconf(nt(name))

    @testset "SizeDiGraph" begin

        @testset "SizeAbsorb" begin
            v0 = iv(3)
            v1 = av(v0, 4, "v1")
            v2 = av(v1, 5, "v2")

            g = SizeDiGraph(v1)

            @test nv(g) == 2
            @test ne(g) == 1

            g = fullgraph(v0)

            @test nv(g) == 3
            @test ne(g) == 2

            @test g[1,:vertex] == v0
            @test g[2,:vertex] == v1

            for e in edges(g)
                @test nout(g[e.src,:vertex]) == get_prop(g, e, :size)
            end
        end

        @testset "SizeStack" begin
            v0 = iv(3, "in1")
            v1 = iv(4, "in2")
            v2 = cc(v0, v1, name="v2")

            g = fullgraph(v0)

            @test nv(g) == 3
            @test ne(g) == 2

            @test g[1,:vertex] == v0
            @test g[2,:vertex] == v2

            for e in edges(g)
                @test nout(g[e.src,:vertex]) == get_prop(g, e, :size)
            end
        end

        @testset "SizeInvariant" begin
            v0 = iv(3)
            v1 = tv(v0, "v1")
            v2 = nc("v2") >> v0 + v1

            g = fullgraph(v0)

            @test nv(g) == 3
            @test ne(g) == 3

            @test g[1,:vertex] == v0
            @test g[2,:vertex] == v1
            @test g[3,:vertex] == v2

            for e in edges(g)
                @test nout(g[e.src,:vertex]) == get_prop(g, e, :size)
            end
        end
    end

    @testset "ΔSizeGraph" begin
        changed(v) = filter(vi -> nout_org(vi) != nout(vi) || nin(vi) != nin_org(vi), all_in_graph(v))

        function testedge(::Input, src, dst, size)
            @test nout_org(src) == size
            @test nin_org(dst)[inputs(dst) .== src] == [size]
        end
        function testedge(::Output, src, dst, size)
            @test nout_org(dst) == size
            @test nin_org(src)[inputs(src) .== dst] == [size]
        end

        @testset "ΔnoutSizeGraph" begin
            v0 = iv(3)
            v1 = av(v0, 4, "v1")
            v2 = av(v0, 5, "v2")
            v3 = cc(v1, v2, name="v3")
            v4 = av(v1, 3, "v4")
            v5 = av(v3, 3, "v5")
            v6 = nc("v6") >> v4 + v5
            v7 = cc(v6, v2, name="v7")

            g = ΔnoutSizeGraph(v5)
            Δnout(v5, 2)

            clist = changed(v5)
            glist = NaiveNASlib.vertexproplist(g, :vertex)
            for v in clist
                @test v in glist
            end

            for e in edges(g)
                testedge(get_prop(g, e, :direction), g[e.src,:vertex], g[e.dst, :vertex], get_prop(g, e, :size))
            end

        end

        @testset "ΔninSizeGraph" begin
            v0 = iv(3)
            v1 = av(v0, 4, "v1")
            v2 = av(v0, 5, "v2")
            v3 = cc(v1, v2, name="v3")
            v4 = av(v0, 9, "v4")
            v5 = nc("v5") >> v4 + v3
            v6 = av(v5, 3, "v6")

            g = ΔninSizeGraph(v3)
            Δnin(v3, 2, 1)

            clist = changed(v3)
            @test Set(clist) == Set(NaiveNASlib.vertexproplist(g, :vertex))

            for e in edges(g)
                testedge(get_prop(g, e, :direction), g[e.src,:vertex], g[e.dst, :vertex], get_prop(g, e, :size))
            end
        end
    end
end


text:
using MOCNeutronTransport
include("../src/AngularQuadrature.jl")
@testset "AngularQuadrature" begin
    for type in [Float32, Float64]
        @testset "Chebyshev" begin
            # chebyshev_angular_quadrature
            angles, weights = generate_chebyshev_angular_quadrature(1, type)
            @test size(angles, 1) == size(weights, 1) == 1
            @test angles[1] == π/type(4)
            @test weights[1] == type(1)
            @test typeof(angles[1]) == typeof(weights[1]) == typeof(type(1))

            angles, weights = generate_chebyshev_angular_quadrature(2, type)
            @test size(angles, 1) == size(weights, 1) == 2
            @test angles  == type[3π/8, π/8]
            @test weights == type[1//2, 1//2]

            angles, weights = generate_chebyshev_angular_quadrature(3, type)
            @test size(angles, 1) == size(weights, 1) == 3
            @test angles  == type[5π/12, 3π/12, π/12]
            @test weights == type[1/3,  1/3,   1/3]
        end

        @testset "Chebyshev-Chebyshev" begin
            # Default data type is Float64
            quadrature = generate_angular_quadrature("Chebyshev-Chebyshev", 3, 2)
            @test quadrature isa ProductAngularQuadrature
            @test size(quadrature.γ, 1) == size(quadrature.w_γ, 1) == 6
            @test size(quadrature.θ, 1) == size(quadrature.w_θ, 1) == 2
            @test quadrature.γ   == Tuple([5π/12, 3π/12, π/12, 7π/12, 9π/12, 11π/12])
            @test quadrature.w_γ == Tuple([1/6, 1/6, 1/6, 1/6, 1/6, 1/6])
            @test quadrature.θ   == Tuple([3π/8, π/8])
            @test quadrature.w_θ == Tuple([1/2, 1/2])

            # Specified data type
            quadrature = generate_angular_quadrature("Chebyshev-Chebyshev", 3, 2, T=type)
            @test quadrature isa ProductAngularQuadrature
            @test size(quadrature.γ, 1) == size(quadrature.w_γ, 1) == 6
            @test size(quadrature.θ, 1) == size(quadrature.w_θ, 1) == 2
            γ_ref = Tuple(type[5π/12, 3π/12, π/12, 7π/12, 9π/12, 11π/12])
            for i = 1:6
                @test quadrature.γ[i] ≈ γ_ref[i]
            end
            @test quadrature.w_γ == Tuple(type[1/6, 1/6, 1/6, 1/6, 1/6, 1/6])
            @test quadrature.θ   == Tuple(type[3π/8, π/8])
            @test quadrature.w_θ == Tuple(type[1/2, 1/2])
        end
    end
end

max_stars_repo_name: khurrumsaleem/MOCNeutronTransport

max_issues_repo_name: khurrumsaleem/MOCNeutronTransport

max_forks_repo_name: khurrumsaleem/MOCNeutronTransport


text:
# automatically generated by the FlatBuffers compiler, do not modify

NamespaceA.NamespaceB.eval(quote

FlatBuffers.@with_kw mutable struct TableInNestedNS
    foo::Int32 = 0
end
FlatBuffers.@ALIGN(TableInNestedNS, 1)
FlatBuffers.slot_offsets(::Type{T}) where {T<:TableInNestedNS} = [
    0x00000004
]

TableInNestedNS(buf::AbstractVector{UInt8}) = FlatBuffers.read(TableInNestedNS, buf)
TableInNestedNS(io::IO) = FlatBuffers.deserialize(io, TableInNestedNS)

end)

max_stars_repo_name: rjkat/flatbuffers-julia

max_issues_repo_name: rjkat/flatbuffers-julia

max_forks_repo_name: rjkat/flatbuffers-julia