text:
import ClimaCore:
    Fields,
    Domains,
    Topologies,
    Meshes,
    DataLayouts,
    Operators,
    Geometry,
    Spaces

using OrdinaryDiffEq: ODEProblem, solve, SSPRK33

import Logging
import TerminalLoggers
Logging.global_logger(TerminalLoggers.TerminalLogger())
const FT = Float64

a = FT(0.0)
b = FT(4pi)
n = 128
α = FT(0.1)

domain = Domains.IntervalDomain(
    Geometry.ZPoint{FT}(a),
    Geometry.ZPoint{FT}(b),
    boundary_names = (:left, :right),
)
mesh = Meshes.IntervalMesh(domain, nelems = n)

cs = Spaces.CenterFiniteDifferenceSpace(mesh)
fs = Spaces.FaceFiniteDifferenceSpace(cs)

V = Geometry.WVector.(ones(FT, fs))
θ = sin.(Fields.coordinate_field(cs).z)

# Solve advection Equation: ∂θ/dt = -∂(vθ)

# upwinding
function tendency1!(dθ, θ, _, t)
    fcc = Operators.FluxCorrectionC2C(
        left = Operators.Extrapolate(),
        right = Operators.Extrapolate(),
    )
    fcf = Operators.FluxCorrectionF2F(
        left = Operators.Extrapolate(),
        right = Operators.Extrapolate(),
    )
    UB = Operators.UpwindBiasedProductC2F(
        left = Operators.SetValue(sin(a - t)),
        right = Operators.SetValue(sin(b - t)),
    )
    ∂ = Operators.DivergenceF2C()

    return @. dθ = -∂(UB(V, θ))
end
function tendency2!(dθ, θ, _, t)
    fcc = Operators.FluxCorrectionC2C(
        left = Operators.Extrapolate(),
        right = Operators.Extrapolate(),
    )
    fcf = Operators.FluxCorrectionF2F(
        left = Operators.Extrapolate(),
        right = Operators.Extrapolate(),
    )
    UB = Operators.UpwindBiasedProductC2F(
        left = Operators.SetValue(sin(a - t)),
        right = Operators.SetValue(sin(b - t)),
    )
    ∂ = Operators.DivergenceF2C()
    return @. dθ = -∂(UB(V, θ)) + fcc(V, θ)
end
# use the advection operator
function tendency3!(dθ, θ, _, t)

    fcc = Operators.FluxCorrectionC2C(
        left = Operators.Extrapolate(),
        right = Operators.Extrapolate(),
    )
    fcf = Operators.FluxCorrectionF2F(
        left = Operators.Extrapolate(),
        right = Operators.Extrapolate(),
    )
    A = Operators.AdvectionC2C(
        left = Operators.SetValue(sin(-t)),
        right = Operators.Extrapolate(),
    )
    return @. dθ = -A(V, θ)
end
# use the advection operator
function tendency4!(dθ, θ, _, t)

    fcc = Operators.FluxCorrectionC2C(
        left = Operators.Extrapolate(),
        right = Operators.Extrapolate(),
    )
    fcf = Operators.FluxCorrectionF2F(
        left = Operators.Extrapolate(),
        right = Operators.Extrapolate(),
    )
    A = Operators.AdvectionC2C(
        left = Operators.SetValue(sin(-t)),
        right = Operators.Extrapolate(),
    )
    return @. dθ = -A(V, θ) + fcc(V, θ)
end

# use the advection operator

@show tendency1!(similar(θ), θ, nothing, 0.0)
# Solve the ODE operator
Δt = 0.001
prob1 = ODEProblem(tendency1!, θ, (0.0, 10.0))
prob2 = ODEProblem(tendency2!, θ, (0.0, 10.0))
prob3 = ODEProblem(tendency3!, θ, (0.0, 10.0))
prob4 = ODEProblem(tendency4!, θ, (0.0, 10.0))
sol1 = solve(
    prob1,
    SSPRK33(),
    dt = Δt,
    saveat = 10 * Δt,
    progress = true,
    progress_message = (dt, u, p, t) -> t,
);
sol2 = solve(
    prob2,
    SSPRK33(),
    dt = Δt,
    saveat = 10 * Δt,
    progress = true,
    progress_message = (dt, u, p, t) -> t,
);
sol3 = solve(
    prob3,
    SSPRK33(),
    dt = Δt,
    saveat = 10 * Δt,
    progress = true,
    progress_message = (dt, u, p, t) -> t,
);
sol4 = solve(
    prob4,
    SSPRK33(),
    dt = Δt,
    saveat = 10 * Δt,
    progress = true,
    progress_message = (dt, u, p, t) -> t,
);

ENV["GKSwstype"] = "nul"
using ClimaCorePlots, Plots
Plots.GRBackend()

dirname = "advect"
path = joinpath(@__DIR__, "output", dirname)
mkpath(path)

anim = Plots.@animate for u in sol1.u
    Plots.plot(u, xlim = (-1, 1))
end
Plots.mp4(anim, joinpath(path, "UBP_advect.mp4"), fps = 10)
Plots.png(
    Plots.plot(sol1.u[end], xlim = (-1, 1)),
    joinpath(path, "sol1_advect_end.png"),
)

anim = Plots.@animate for u in sol2.u
    Plots.plot(u, xlim = (-1, 1))
end
Plots.mp4(anim, joinpath(path, "UBP_advect_fc.mp4"), fps = 10)
Plots.png(
    Plots.plot(sol2.u[end], xlim = (-1, 1)),
    joinpath(path, "sol2_advect_end.png"),
)

anim = Plots.@animate for u in sol3.u
    Plots.plot(u, xlim = (-1, 1))
end
Plots.mp4(anim, joinpath(path, "C2C_advect.mp4"), fps = 10)
Plots.png(
    Plots.plot(sol3.u[end], xlim = (-1, 1)),
    joinpath(path, "sol3_advect_end.png"),
)

anim = Plots.@animate for u in sol4.u
    Plots.plot(u, xlim = (-1, 1))
end
Plots.mp4(anim, joinpath(path, "C2C_advect_fc.mp4"), fps = 10)
Plots.png(
    Plots.plot(sol4.u[end], xlim = (-1, 1)),
    joinpath(path, "sol4_advect_end.png"),
)

p = Plots.plot(sol1.u[end], xlim = (-1, 1), ls = :dash, label = "UBP")
p = Plots.plot!(sol2.u[end], xlim = (-1, 1), ls = :dot, label = "UBP_FC")
p = Plots.plot!(sol3.u[end], xlim = (-1, 1), ls = :solid, label = "C2C")
p = Plots.plot!(sol4.u[end], xlim = (-1, 1), ls = :dashdot, label = "C2C_FC")
Plots.png(p, joinpath(path, "all_advect_end.png"))

function linkfig(figpath, alt = "")
    # buildkite-agent upload figpath
    # link figure in logs if we are running on CI
    if get(ENV, "BUILDKITE", "") == "true"
        artifact_url = "artifact://$figpath"
        print("\033]1338;url='$(artifact_url)';alt='$(alt)'\a\n")
    end
end

linkfig(
    relpath(joinpath(path, "advect_end.png"), joinpath(@__DIR__, "../..")),
    "Advect End Simulation",
)


text:
module Beancount

using PyCall
using Decimals
using Dates
using Parameters
using Combinatorics

import Base.==

const DEFAULT_CURRENCY = "EUR"
const aDict = Union{Nothing,Dict{Any,Any}}
const SN = Union{Nothing,String}
const DN = Union{Nothing,Decimal}

function __init__()
    global loader, load_file
    loader = pyimport("beancount.loader")
    load_file = loader.load_file
end

export Posting, Transaction, issimilar, balances, isTransaction, from_python, txn, pad

Base.homedir(path::AbstractString...) = joinpath(homedir(),path...)

struct Posting
    location::aDict
    account::String
    value::DN
    currency::SN
    comment::SN
    Posting(;account::String,value::Union{Nothing,Number}=nothing,currency::SN=nothing,comment::SN=nothing,location::aDict=nothing) = new(
        location,account,value == nothing ? value : convert(Decimal,value),currency,comment)
end

function ==(pA::Posting,pB::Posting)
    pA.account == pB.account && pA.value == pB.value && pA.currency == pB.currency
end

function ==(pA::Vector{Posting},pB::Vector{Posting})
    for a in permutations(pA)
        if all((a[i] == pB[i] for i in 1:length(pB)))
            return true
        end
        println(a)
    end
    return false
end

function Base.show(io::IO,p::Posting)
    @unpack account,value,currency,comment = p
    _comment = comment == nothing ? "" : " ; $comment"
    if value == nothing
        print(io,"$account$_comment")
    else
        _currency = currency == nothing ? DEFAULT_CURRENCY : currency
        print(io,"$account  $value $_currency$_comment")
    end
end

abstract type AbstractTransactionType end

struct txn <: AbstractTransactionType end
struct pad <: AbstractTransactionType end

const Collection1D{T} = Union{Vector{T},Set{T},NTuple{N,T}} where N

struct Transaction{typ<:AbstractTransactionType}
    location::aDict
    date::Date
    payee::SN
    description::String
    tags::Set{String}
    links::Set{String}
    postings::Vector{Posting}
    Transaction{typ}(;date::Date,payee::SN=nothing,description::String,
                     tags::Collection1D{String}=Set{String}(),links::Collection1D{String}=Set{String}(),postings::Vector{Posting}=Posting[],location::aDict=nothing) where typ <: AbstractTransactionType = new{typ}(
                         nothing,date,clean_text(payee),clean_text(description),
                         tags isa Set ? tags : Set(tags),
                         links isa Set ? links : Set(links),postings)
end

function ==(txA::Transaction{T},txB::Transaction{T}) where T <: AbstractTransactionType
    txA.date        == txA.date &&
    txA.payee       == txA.payee &&
    txA.description == txA.description &&
    txA.tags        == txA.tags &&
    txA.links       == txA.links &&
    txA.postings    == txA.postings
end

clean_text(::Nothing) = nothing

function clean_text(st::String)
    reduce(replace,[
        "\n" => " ",
        "\t" => " ",
        r" {2,}" => " ",
    ],init=st)
end

function from_python_posting(tup::NTuple{6,Any})
    Posting(
        location = tup[end],
        account=tup[1],
        value=convert(Decimal,tup[2][1]),
        currency=tup[2][2])
end

from_python(::Any) = nothing

function from_python(entry::NTuple{8,Any})
    if entry[3] in ("*","txn")
        typ = txn
    elseif entry[3] == "P"
        typ = pad
    else
        error("Unknown transaction type '$(entry[3])'.")
    end
    Transaction{typ}(location=entry[1],
                     date=entry[2],
                     payee=entry[4],
                     description=entry[5],
                     tags=Set{String}(string.(entry[6])),
                     links=Set{String}(string.(entry[7])),
                     postings=from_python_posting.(entry[8]))
end

function Base.show(io::IO,t::Transaction)
    @unpack date,payee,description,tags,links,postings = t
    _links = join("^".*links)
    _tags = join("#".*tags)
    _payee = payee == nothing ? "" : " \"$payee\""
    print(io,"$date *$_payee \"$description\" $_tags $_links")
    for p in postings
        print(io,"\n")
        print(io," "^4)
        show(io,p)
    end
end

function sum_positive(txA::Transaction)
    mapreduce(+,txA.postings) do p
        (p.value == nothing || p.value < 0) ? Decimal(0) : p.value
    end
end

function sum_negative(txA::Transaction)
    mapreduce(+,txA.postings) do p
        (p.value == nothing || p.value > 0) ? Decimal(0) : p.value
    end
end

function issimilar(txA::Transaction,txB::Transaction)
    (txA.date == txB.date) &&
    (sum_negative(txA) == sum_negative(txB) ||
    sum_positive(txA) == sum_positive(txB ))
end

function balances(tx::Transaction)
    sum_positive(tx) == sum_negative(tx) ||
    count([p.value == nothing for p in tx.postings]) == 1
end

isTransaction(::Transaction) = true
isTransaction(tup) = length(tup) == 8 && (tup[3] == "*" || tup[3] == "txn")

Base.push!(tr::Transaction,p::Posting) = push!(tr.postings,p)
Base.append!(tr::Transaction,p::Vector{Posting}) = append!(tr.postings,p)

end # module

ext: jl


text:
import Pkg.BinaryPlatforms: detect_libgfortran_version, detect_libstdcxx_version, detect_cxxstring_abi
using ObjectFile

"""
    detect_libgfortran_version(oh::ObjectHandle, platform::Platform)

Given an ObjectFile, examine its dynamic linkage to discover which (if any)
`libgfortran` it's linked against.  The major SOVERSION will determine which
GCC version we're restricted to.
"""
function detect_libgfortran_version(oh::ObjectHandle, platform::Platform)
    # We look for linkage to libgfortran
    libs = basename.(path.(DynamicLinks(oh)))
    fortran_libs = filter(l -> occursin("libgfortran", l), libs)
    if isempty(fortran_libs)
        return nothing
    end

    # If we find one, pass it off to Pkg.detect_libgfortran_version()
    return detect_libgfortran_version(first(fortran_libs), platform)
end

function check_libgfortran_version(oh::ObjectHandle, platform::Platform; verbose::Bool = false)
    libgfortran_version = nothing
    try
        libgfortran_version = detect_libgfortran_version(oh, platform)
    catch e
        if isa(e, InterruptException)
            rethrow(e)
        end
        @warn "$(path(oh)) could not be scanned for libgfortran dependency!" exception=(e, catch_backtrace())
        return true
    end

    if verbose && libgfortran_version != nothing
        @info("$(path(oh)) locks us to libgfortran v$(libgfortran_version)")
    end

    if compiler_abi(platform).libgfortran_version === nothing && libgfortran_version != nothing
        msg = strip(replace("""
        $(path(oh)) links to libgfortran!  This causes incompatibilities across
        major versions of GCC.  To remedy this, you must build a tarball for
        each major version of GCC.  To do this, immediately after your `platforms`
        definition in your `build_tarballs.jl` file, add the line:
        """, '\n' => ' '))
        msg *= "\n\n    platforms = expand_gfortran_versions(platforms)"
        @warn(msg)
        return false
    end
    return true
end

"""
    detect_libstdcxx_version(oh::ObjectHandle, platform::Platform)

Given an ObjectFile, examine its dynamic linkage to discover which (if any)
`libgfortran` it's linked against.  The major SOVERSION will determine which
GCC version we're restricted to.
"""
function detect_libstdcxx_version(oh::ObjectHandle, platform::Platform)
    # We look for linkage to libstdc++
    libs = basename.(path.(DynamicLinks(oh)))
    libstdcxx_libs = filter(l -> occursin("libstdc++", l), libs)
    if isempty(libstdcxx_libs)
        return nothing
    end

    # Extract all pieces of `.gnu.version_d` from libstdc++.so, find the `GLIBCXX_*`
    # symbols, and use the maximum version of that to find the GLIBCXX ABI version number
    version_symbols = readmeta(first(libstdcxx_libs)) do oh
        unique(vcat((x -> x.names).(ELFVersionData(oh))...))
    end
    version_symbols = filter(x -> startswith(x, "GLIBCXX_"), version_symbols)
    if isempty(version_symbols)
        # This would be weird, but let's be prepared
        return nothing
    end
    return maximum([VersionNumber(split(v, "_")[2]) for v in version_symbols])
end

function check_libstdcxx_version(oh::ObjectHandle, platform::Platform; verbose::Bool = false)
    libstdcxx_version = nothing

    try
        libstdcxx_version = detect_libstdcxx_version(oh, platform)
    catch e
        if isa(e, InterruptException)
            rethrow(e)
        end
        @warn "$(path(oh)) could not be scanned for libstdcxx dependency!" exception=(e, catch_backtrace())
        return true
    end

    if verbose && libstdcxx_version != nothing
        @info("$(path(oh)) locks us to libstdc++ v$(libstdcxx_version)+")
    end

    # This actually isn't critical, so we don't complain.  Yet.
    # if compiler_abi(platform).libstdcxx_version === nothing && libstdcxx_version != nothing
    #     msg = strip(replace("""
    #     $(path(oh)) links to libstdc++!  This causes incompatibilities across
    #     major versions of GCC.  To remedy this, you must build a tarball for
    #     each major version of GCC.  To do this, immediately after your `products`
    #     definition in your `build_tarballs.jl` file, add the line:
    #     """, '\n' => ' '))
    #     msg *= "\n\n    products = expand_cxx_versions(products)"
    #     warn(io, msg)
    #     return false
    # end
    return true
end

function cppfilt(symbol_names::Vector, platform::Platform)
    input = IOBuffer()
    for name in symbol_names
        println(input, name)
    end
    seekstart(input)

    output = IOBuffer()
    mktemp() do t, io
        ur = preferred_runner()(dirname(t); cwd="/workspace/", platform=platform)
        run_interactive(ur, `/opt/bin/c++filt`; stdin=input, stdout=output)
    end

    return filter!(s -> !isempty(s), split(String(take!(output)), "\n"))
end

"""
    detect_cxxstring_abi(oh::ObjectHandle, platform::Platform)

Given an ObjectFile, examine its symbols to discover which (if any) C++11
std::string ABI it's using.  We do this by scanning the list of exported
symbols, triggering off of instances of `St7__cxx11` or `_ZNSs` to give
evidence toward a constraint on `cxx11`, `cxx03` or neither.
"""
function detect_cxxstring_abi(oh::ObjectHandle, platform::Platform)
    try
        # First, if this object doesn't link against `libstdc++`, it's a `:cxxany`
        if !any(occursin("libstdc++", l) for l in ObjectFile.path.(DynamicLinks(oh)))
            return nothing
        end

        symbol_names = cppfilt(symbol_name.(Symbols(oh)), platform)
        # Shove the symbol names through c++filt (since we don't want to have to
        # reimplement the parsing logic in Julia).  If anything has `cxx11` tags,
        # then mark it as such.
        if any(occursin("[abi:cxx11]", c) || occursin("std::__cxx11", c) for c in symbol_names)
            return :cxx11
        end
        # Otherwise, if we still have `std::string`'s or `std::list`'s in there, it's implicitly a
        # `cxx03` binary, even though we don't have a __cxx03 namespace or something.  Mark it.
        if any(occursin("std::string", c) || occursin("std::basic_string", c) ||
               occursin("std::list", c) for c in symbol_names)
            return :cxx03
        end
    catch e
        if isa(e, InterruptException)
            rethrow(e)
        end
        @warn "$(path(oh)) could not be scanned for cxx11 ABI!" exception=(e, catch_backtrace())
    end
    return nothing
end

function check_cxxstring_abi(oh::ObjectHandle, platform::Platform; io::IO = stdout, verbose::Bool = false)
    # First, check the stdlibc++ string ABI to see if it is a superset of `platform`.  If it's
    # not, then we have a problem!
    cxx_abi = detect_cxxstring_abi(oh, platform)

    # If no std::string symbols found, just exit out immediately
    if cxx_abi == nothing
        return true
    end

    if verbose && cxx_abi != nothing
        @info("$(path(oh)) locks us to $(cxx_abi)")
    end

    if compiler_abi(platform).cxxstring_abi == nothing && cxx_abi != nothing
        msg = strip(replace("""
        $(path(oh)) contains std::string values!  This causes incompatibilities across
        the GCC 4/5 version boundary.  To remedy this, you must build a tarball for
        both GCC 4 and GCC 5.  To do this, immediately after your `products`
        definition in your `build_tarballs.jl` file, add the line:
        """, '\n' => ' '))
        msg *= "\n\n    products = expand_cxx_versions(products)"
        @warn(msg)
        return false
    end

    if compiler_abi(platform).cxxstring_abi != cxx_abi
        msg = strip(replace("""
        $(path(oh)) contains $(cxx_abi) ABI std::string values within its public interface,
        but we are supposedly building for $(compiler_abi(platform).cxxstring_abi) ABI. This usually
        indicates that the build system is somehow ignoring our choice of compiler, as we manually
        insert the correct compiler flags for this ABI choice!
        """, '\n' => ' '))
        @warn(msg)
        return false
    end
    return true
end


text:
# use only activations without "kinks" for numerical gradient checking
# see e.g. https://stackoverflow.com/questions/40623512/how-to-check-relu-gradient
const ACTIVATIONS = [identity, σ, swish, softplus, logcosh, mish, tanhshrink, lisht]

@testset "testing simple matrix model" begin
    layerbuilder(k) = Flux.Dense(k, 2, NNlib.relu)
    x = ArrayNode(randn(Float32, 4, 5))
    m = reflectinmodel(x, layerbuilder)
    @test size(m(x).data) == (2, 5)
    @test m isa ArrayModel
    @test eltype(m(x).data) == Float32
    @inferred m(x)
end

@testset "testing simple aggregation model" begin
    layerbuilder(k) = Flux.Dense(k, 2, NNlib.relu)
    x = BagNode(ArrayNode(randn(Float32, 4, 4)), [1:2, 3:4])
    m = reflectinmodel(x, layerbuilder)
    @test size(m(x).data) == (2, 2)
    @test m isa BagModel
    @test eltype(m(x).data) == Float32
    @inferred m(x)
end

@testset "testing simple tuple models" begin
    layerbuilder(k) = Flux.Dense(k, 2, NNlib.relu)
    x = ProductNode((a=ArrayNode(randn(Float32, 3, 4)), b=ArrayNode(randn(Float32, 4, 4))))
    m = reflectinmodel(x, layerbuilder)
    @test eltype(m(x).data) == Float32
    @test size(m(x).data) == (2, 4)
    @test m isa ProductModel
    @test m.ms[:a] isa ArrayModel
    @test m.ms[:b] isa ArrayModel
    @inferred m(x)

    x = ProductNode((BagNode(ArrayNode(randn(Float32, 3, 4)), [1:2, 3:4]),
                     BagNode(ArrayNode(randn(Float32, 4, 4)), [1:1, 2:4])))
    m = reflectinmodel(x, layerbuilder)
    @test size(m(x).data) == (2, 2)
    @test m isa ProductModel
    @test m.ms[1] isa BagModel
    @test m.ms[1].im isa ArrayModel
    @test m.ms[1].bm isa ArrayModel
    @test m.ms[2] isa BagModel
    @test m.ms[2].im isa ArrayModel
    @test m.ms[2].bm isa ArrayModel
    @inferred m(x)
end

@testset "testing nested bag model" begin
    bn = BagNode(ArrayNode(randn(Float32, 2, 8)), [1:1, 2:2, 3:6, 7:8])
    x = BagNode(bn, [1:2, 3:4])
    m = reflectinmodel(x, d -> Flux.Dense(d, 2))
    @test size(m(x).data) == (2, 2)
    @test m isa BagModel
    @test m.im isa BagModel
    @test m.im.im isa ArrayModel
    @test m.im.bm isa ArrayModel
    @test m.bm isa ArrayModel
    @test eltype(m(x).data) == Float32
    @inferred m(x)

    a = BagNode(BagNode(ArrayNode(randn(2,2)),[1:2]),[1:1])
    b = BagNode(missing,[0:-1])
    c = BagNode(a.data[1:0], [0:-1])
    m = reflectinmodel(a)
    abc = catobs(a, b, c)
    bca = catobs(b, c, a)
    ma = m(a).data
    mb = m(b).data
    mc = m(c).data
    mabc = m(abc).data
    mbca = m(bca).data
    @test mb ≈ mc
    @test mabc[:,1] ≈ ma
    @test mabc[:,2] ≈ mb
    @test mabc[:,3] ≈ mc
    @test mbca[:,1] ≈ mb
    @test mbca[:,2] ≈ mc
    @test mbca[:,3] ≈ ma
    @inferred m(a)
    @inferred m(b)
    @inferred m(c)
end

@testset "reflectinmodel for missing + all_imputing" begin
    f1 = x -> reflectinmodel(x)
    f2 = x -> reflectinmodel(x, d -> Flux.Chain(Dense(d, 10), Dense(10, 10)))
    f3 = x -> reflectinmodel(x; all_imputing=true)
    _get_first(x::ArrayModel{<:Flux.Chain}) = x.m[1]
    _get_first(x) = x.m

    x1 = maybehot(2, 1:3) |> ArrayNode
    x2 = maybehot(missing, 1:3) |> ArrayNode
    for f in [f1, f2, f3], x in [x1, x2]
        @test _get_first(f(x)) isa PostImputingDense
        @inferred f(x)(x)
    end

    x1 = maybehotbatch([1, 2, 3], 1:3) |> ArrayNode
    x2 = maybehotbatch([1, 2, missing], 1:3) |> ArrayNode
    x3 = maybehotbatch(fill(missing, 3), 1:3) |> ArrayNode
    for f in [f1, f2, f3], x in [x1, x2, x3]
        @test _get_first(f(x)) isa PostImputingDense
        @inferred f(x)(x)
    end

    x1 = NGramMatrix(["a", "b", "c"]) |> ArrayNode
    x2 = NGramMatrix(["a", missing, "c"]) |> ArrayNode
    x3 = NGramMatrix(fill(missing, 3)) |> ArrayNode
    for f in [f1, f2]
        @test _get_first(f(x1)) isa Flux.Dense
        @test _get_first(f(x2)) isa PostImputingDense
        @test _get_first(f(x3)) isa PostImputingDense
        for x in [x1, x2, x3] @inferred f(x)(x) end
    end
    for x in [x1, x2, x3]
        @test _get_first(f3(x)) isa PostImputingDense
        @inferred f3(x)(x)
    end

    x1 = rand([1, 2], 3, 3) |> ArrayNode
    x2 = rand([1, 2, missing], 3, 3) |> ArrayNode
    x3 = fill(missing, 3, 3) |> ArrayNode
    for f in [f1, f2]
        @test _get_first(f(x1)) isa Flux.Dense
        @test _get_first(f(x2)) isa PreImputingDense
        @test _get_first(f(x3)) isa PreImputingDense
        for x in [x1, x2, x3] @inferred f(x)(x) end
    end
    for x in [x1, x2, x3]
        @test _get_first(f3(x)) isa PreImputingDense
        @inferred f3(x)(x)
    end
end

# pn.m should be identity for any product node pn with a single key
@testset "single key dictionary reflect in model" begin
    layerbuilder(k) = Flux.Dense(k, 2, NNlib.relu)
    fsm = Dict("" => layerbuilder)

    x1 = (ArrayNode(randn(Float32, 3, 4)),) |> ProductNode
    x2 = (a = ArrayNode(randn(Float32, 3, 4)),) |> ProductNode
    x3 = (a = ArrayNode(randn(Float32, 3, 4)), b = ArrayNode(randn(Float32, 3, 4))) |> ProductNode

    m1 = reflectinmodel(x1, layerbuilder; single_key_identity=false)
    m1_ski = reflectinmodel(x1, layerbuilder; single_key_identity=true)
    m1_ski_fsm = reflectinmodel(x1, layerbuilder; fsm=fsm, single_key_identity=true)
    m2 = reflectinmodel(x2, layerbuilder; single_key_identity=false)
    m2_ski = reflectinmodel(x2, layerbuilder; single_key_identity=true)
    m2_ski_fsm = reflectinmodel(x2, layerbuilder; fsm=fsm, single_key_identity=true)
    m3 = reflectinmodel(x3, layerbuilder; single_key_identity=false)
    m3_ski = reflectinmodel(x3, layerbuilder; single_key_identity=true)
    m3_ski_fsm = reflectinmodel(x3, layerbuilder; fsm=fsm, single_key_identity=true)

    for m in [m1, m1_ski, m1_ski_fsm]
        @test eltype(m(x1).data) == Float32
        @test size(m(x1).data) == (2, 4)
        @inferred m(x1)
        @test m isa ProductModel
        @test m.ms[1] isa ArrayModel
    end

    for m in [m2, m2_ski, m2_ski_fsm]
        @test eltype(m(x2).data) == Float32
        @test size(m(x2).data) == (2, 4)
        @inferred m(x2)
        @test m isa ProductModel
        @test m.ms[1] isa ArrayModel
    end

    for m in [m3, m3_ski, m3_ski_fsm]
        @test eltype(m(x3).data) == Float32
        @test size(m(x3).data) == (2, 4)
        @inferred m(x3)
        @test m isa ProductModel
        @test m.ms[1] isa ArrayModel
        @test m.ms[2] isa ArrayModel
    end

    @test m1.m isa ArrayModel{<:Dense}
    @test m2.m isa ArrayModel{<:Dense}
    @test m3.m isa ArrayModel{<:Dense}
    @test m1_ski.m isa IdentityModel
    @test m2_ski.m isa IdentityModel
    @test m3_ski.m isa ArrayModel{<:Dense}
    # fsm overrides ski
    @test m1_ski_fsm.m isa ArrayModel{<:Dense}
    @test m2_ski_fsm.m isa ArrayModel{<:Dense}
    @test m3_ski_fsm.m isa ArrayModel{<:Dense}
end

# array model for matrices with one row should implement identity
@testset "single scalar as identity" begin
    layerbuilder(k) = Flux.Dense(k, 2, NNlib.relu)
    x1 = ArrayNode(randn(Float32, 1, 3))
    x2 = BagNode(ArrayNode(randn(Float32, 1, 5)), [1:2, 3:5])
    x3 = (a = ArrayNode(randn(Float32, 1, 4)), b = ArrayNode(randn(Float32, 2, 4))) |> ProductNode

    m1 = reflectinmodel(x1, layerbuilder; single_scalar_identity=false)
    m1_sci = reflectinmodel(x1, layerbuilder; single_scalar_identity=true)
    m2 = reflectinmodel(x2, layerbuilder; single_scalar_identity=false)
    m2_sci = reflectinmodel(x2, layerbuilder; single_scalar_identity=true)
    m3 = reflectinmodel(x3, layerbuilder; single_scalar_identity=false)
    m3_sci = reflectinmodel(x3, layerbuilder; single_scalar_identity=true)

    @test size(m1(x1).data) == (2, 3)
    @test eltype(m1(x1).data) == Float32
    @test m1 isa ArrayModel
    @test size(m1_sci(x1).data) == (1, 3)
    @test eltype(m1_sci(x1).data) == Float32
    @test m1_sci isa ArrayModel
    @inferred m1(x1)

    @test size(m2(x2).data) == (2, 2)
    @test eltype(m2(x2).data) == Float32
    @test m2 isa BagModel
    @test m2.im isa ArrayModel{<:Dense}
    @test size(m2_sci(x2).data) == (2, 2)
    @test eltype(m2_sci(x2).data) == Float32
    @test m2_sci isa BagModel
    @test m2_sci.im isa IdentityModel
    @inferred m2(x2)

    @test size(m3(x3).data) == (2, 4)
    @test eltype(m3(x3).data) == Float32
    @test m3 isa ProductModel
    @test m3.ms[1] isa ArrayModel{<:Dense}
    @test m3.ms[2] isa ArrayModel{<:Dense}
    @test size(m3_sci(x3).data) == (2, 4)
    @test eltype(m3_sci(x3).data) == Float32
    @test m3_sci isa ProductModel
    @test m3_sci.ms[1] isa IdentityModel
    @test m3_sci.ms[2] isa ArrayModel{<:Dense}
    @inferred m3(x3)
end

@testset "model aggregation grad check w.r.t. inputs" begin
    for (bags1, bags2, bags3) in BAGS3
        layerbuilder(k) = Dense(k, 2, rand(ACTIVATIONS)) |> f64
        abuilder(d) = BagCount(all_aggregations(Float64, d))
        x = randn(4, 4)
        y = randn(3, 4)
        z = randn(2, 8)

        for ds in [x -> ArrayNode(x),
                   x -> BagNode(ArrayNode(x), bags1)]
            m = reflectinmodel(ds(x), layerbuilder, abuilder)
            @inferred m(ds(x))
            @test gradtest(x -> m(ds(x)).data, x)
        end

        for ds in [(x, y) -> ProductNode((ArrayNode(x), ArrayNode(y))),
                   (x, y) -> ProductNode((a=BagNode(ArrayNode(x), bags1), b=BagNode(ArrayNode(y), bags2)))]
            m = reflectinmodel(ds(x, y), layerbuilder, abuilder)
            @inferred m(ds(x, y))
            @test gradtest((x, y) -> m(ds(x, y)).data, x, y)
        end

        ds = z -> BagNode(BagNode(ArrayNode(z), bags3), bags1)
        m = reflectinmodel(ds(z), layerbuilder, abuilder)
        @inferred m(ds(z))
        @test gradtest(z -> m(ds(z)).data, z)
    end
end

@testset "model aggregation grad check w.r.t. inputs weighted" begin
    for (bags1, bags2, bags3) in BAGS3
        layerbuilder(k) = Dense(k, 2, rand(ACTIVATIONS)) |> f64
        abuilder(d) = BagCount(all_aggregations(Float64, d))
        x = randn(4, 4)
        y = randn(3, 4)
        z = randn(2, 8)
        w1 = abs.(randn(4)) .+ 0.1
        w2 = abs.(randn(4)) .+ 0.1
        w3 = abs.(randn(8)) .+ 0.1

        ds = x -> WeightedBagNode(ArrayNode(x), bags1, w1)
        m = reflectinmodel(ds(x), layerbuilder, abuilder)
        @inferred m(ds(x))
        @test gradtest(x -> m(ds(x)).data, x)

        ds = (x, y) -> ProductNode((WeightedBagNode(ArrayNode(x), bags1, w1),
                                    WeightedBagNode(ArrayNode(y), bags2, w2)))
        m = reflectinmodel(ds(x, y), layerbuilder, abuilder)
        @inferred m(ds(x, y))
        @test gradtest((x, y) -> m(ds(x, y)).data, x, y)

        ds = z -> WeightedBagNode(WeightedBagNode(ArrayNode(z), bags3, w3), bags1, w1)
        m = reflectinmodel(ds(z), layerbuilder, abuilder)
        @inferred m(ds(z))
        @test gradtest(z -> m(ds(z)).data, z)
    end
end

@testset "model aggregation grad check w.r.t. params" begin
    for (bags1, bags2, bags3) in BAGS3
        layerbuilder(k) = Dense(k, 2, rand(ACTIVATIONS)) |> f64
        abuilder(d) = BagCount(all_aggregations(Float64, d))
        x = randn(4, 4)
        y = randn(3, 4)
        z = randn(2, 8)

        for ds in [
                   ArrayNode(x),
                   BagNode(ArrayNode(x), bags1),
                   ProductNode((ArrayNode(x), ArrayNode(y))),
                   ProductNode((a=BagNode(ArrayNode(y), bags1), b=BagNode(ArrayNode(x), bags2))),
                   BagNode(BagNode(ArrayNode(z), bags3), bags1)
                  ]
            m = reflectinmodel(ds, layerbuilder, abuilder)
            @inferred m(ds)
            @test gradtest(() -> m(ds).data, Flux.params(m))
        end
    end
end

@testset "model aggregation grad check w.r.t. params weighted" begin
    for (bags1, bags2, bags3) in BAGS3
        layerbuilder(k) = Dense(k, 2) |> f64
        abuilder(d) = BagCount(all_aggregations(Float64, d))
        x = randn(4, 4)
        y = randn(3, 4)
        z = randn(2, 8)
        w1 = abs.(randn(4)) .+ 0.1
        w2 = abs.(randn(4)) .+ 0.1
        w3 = abs.(randn(8)) .+ 0.1

        for ds in [
                   WeightedBagNode(ArrayNode(x), bags1, w1),
                   ProductNode((WeightedBagNode(ArrayNode(x), bags1, w1),
                                WeightedBagNode(ArrayNode(y), bags2, w2))) ,
                   WeightedBagNode(WeightedBagNode(ArrayNode(z), bags3, w3), bags1, w1)
                  ]
            m = reflectinmodel(ds, layerbuilder, abuilder)
            @inferred m(ds)
            @test gradtest(() -> m(ds).data, Flux.params(m))
        end
    end
end


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

using BAT
using Test

using Random, StatsBase, Distributions, ArraysOfArrays

import UltraNest, PyCall

# Disable UltraNest progress output
# ToDo: Find a cleaner way to do this.
PyCall.py"""
import os
import sys
sys.stdout = open(os.devnull, 'w')
"""

@testset "test_ultranest" begin
    dist = product_distribution([
        MixtureModel([truncated(Normal(-1, 0.1), -2, 0), truncated(Normal(1, 0.1), 0, 2)], [0.5, 0.5]),
        MixtureModel([truncated(Normal(-2, 0.25), -3, -1), truncated(Normal(2, 0.25), 1, 3)], [0.3, 0.7]),
        MixtureModel([truncated(Normal(-5, 0.25), -6, -4), truncated(Normal(5, 0.25), 4, 6)], [0.2, 0.8]),
    ])

    prior = product_distribution(Uniform.(minimum.(dist.v), maximum.(dist.v)))

    likelihood = let dist = dist
        function (v::AbstractVector{<:Real})
            ll = logpdf(dist, v)
            # lofpdf on MixtureModel returns NaN in gaps between distributions, and UltraNest
            # doesn't like -Inf, so return -1E10
            T = promote_type(Float32, typeof(ll))
            # isnan(ll) here only required for Distributions < v0.25
            (log = isnan(ll) || isinf(ll) && ll < 0 ? T(-1E10) : T(ll),)
        end
    end

    posterior = PosteriorDensity(likelihood, prior)
    algorithm = ReactiveNestedSampling()
    r = BAT.sample_and_verify(posterior, algorithm, dist)
    @test r.verified

    smpls = r.result
    @test logdensityof(posterior).(smpls.v) ≈ smpls.logd

    uwsmpls = r.uwresult
    @test logdensityof(posterior).(uwsmpls.v) ≈ uwsmpls.logd
    @test all(isequal(1), uwsmpls.weight)

    logz_expected = -log(prod(maximum.(prior.v) .- minimum.(prior.v)))
    @test isapprox(r.logintegral.val, logz_expected, atol = 10 * r.logintegral.err)

    @test r.ess > 50
end


text:
module ClusterOfGeneMutation

const data_dir = "../data"

# package code goes here
include("preprocess.jl")

end # module


text:
function main()
    N, K = parse.(Int, split(readline()))
    as = parse.(Int, split(readline()))
    bs = parse.(Int, split(readline()))

    sort!(bs)

    l, r = 0, maximum(as) * bs[end]
    while r - l > 1
        m = (l + r) ÷ 2
        cnt = 0
        for a in as
            cnt += searchsortedlast(bs, m ÷ a)
        end

        if cnt ≥ K
            r = m
        else
            l = m
        end
    end

    println(r)
end

main()


text:
#### Bibliography parser

# import fileinput
# from pprint import pprint

export Parser
export clear_comments

using JSON

"""Main struct for Bibtex parsing"""
struct Parser
    data
    token
    token_type
    _next_token
    hashtable
    mode
    records
    line
    white
    nl
    token_re
end

"""Return the bibtex content without comments"""
function clear_comments(data)
    res = replace(data, r"(%.*\n)" => "")
    res = replace(res, r"(comment [^\n]*\n)" => "")
    return res
end

"""Returns a token iterator"""
function tokenize(self)
    for item in self.token_re.finditer(self.data)
        i = item.group(0)
        if self.white.match(i)
            if self.nl.match(i)
                self.line += 1
            end
            continue
        else
            yield i
        end
    end
end

"""Parses self.data and stores the parsed bibtex to self.rec"""
function parse(self)
    while true
        try
            self.next_token()
            while self.database()
            end
        catch
            error("StopIteration")
        end
    end
end

"""Outer constructor"""
function Parser(self, data)
    token = nothing
    token_type = nothing
    token_re = r"([^\s\"#%'(){}@,=]+|\n|@|\"|{|}|=|,)"
    _next_token = self.tokenize().next
    hashtable = Dict()
    mode = nothing
    records = Dict()
    line = 1

    # compile some regexes
    white = r"[\n|\s]+"
    nl = r"[\n]"
    return Parser(data,
        token,
        token_type,
        _next_token,
        hashtable,
        mode,
        records,
        line,
        white,
        nl,
        token_re)
end

"""Returns next token"""
function next_token(self)
    self.token = self._next_token()
    #print self.line, self.token
end

"""Database"""
function database(self)
    if self.token == "@"
        self.next_token()
        self.entry()
    end
end

"""Entry"""
function entry(self)
    if lowercase(self.token) == "string"
        self.mode = "string"
        self.string()
    else
        self.mode = "record"
        self.record()
    end
    self.mode = nothing
end

"""String"""
function string(self)
    if lowercase(self.token) == "string"
        self.next_token()
        if self.token == "{"
            self.next_token()
            self.field()

            self.token == "}" || error("} missing")
        end
    end
end

"""Field"""
function field(self)
    name = self.name()
    if self.token == '='
        self.next_token()
        value = self.value()
        if self.mode == 'string'
            self.hashtable[name] = value
        end
        return (name, value)
    end
end

"""Value"""
function value(self)
    value = ""
    val = []

    while true
        if self.token == '"'
            while true
                self.next_token()
                if self.token == '"'
                    break
                else
                    push!(val, self.token)
                end
            end
            if self.token == '"'
                self.next_token()
            else
                error("\" missing")
            end
        elseif self.token == '{'
            brac_counter = 0
            while true
                self.next_token()
                if self.token == '{'
                    brac_counter += 1
                end
                if self.token == '}'
                    brac_counter -= 1
                end
                if brac_counter < 0
                    break
                else
                    push!(val, self.token)
                end
            end
            if self.token == '}'
                self.next_token()
            else
                error("} missing")
            end
        elseif self.token != "=" && re.match(r"\w|#|,", self.token)
            value = self.query_hashtable(self.token)
            push!(val, value)
            while true
                self.next_token()
                # if token is in hashtable then replace
                value = self.query_hashtable(self.token)
                if re.match(r"[^\w#]|,|}|{", self.token) #self.token == '' :
                    break
                else
                    push!(val, value)
                end
            end

        elseif self.token.isdigit()
            value = self.token
            self.next_token()
        else
            if self.token in self.hashtable
                value = self.hashtable[ self.token ]
            else
                value = self.token
            end
            self.next_token()
        end

        if re.match(r"}|,", self.token)
            break
        end

    value = join(val, " ")
    return value
end

query_hashtable( self, s ) =
  s in self.hashtable ? self.hashtable[ self.token ] : s

"""Returns parsed Name"""
function name(self)
    name = self.token
    self.next_token()
    return name
end

"""Returns parsed Key"""
function key(self)
    key = self.token
    self.next_token()
    return key
end

"""Record"""
function record(self)
    if !(self.token in ["comment", "string", "preamble"])
        record_type = self.token
        self.next_token()
        if self.token == '{'
            self.next_token()
            key = self.key()
            self.records[ key ] = Dict()
            self.records[ key ]["type"] = record_type
            self.records[ key ]["id"] = key
            if self.token == ","
                while true
                    self.next_token()
                    field = self.field()
                    if field
                        k = field[0]
                        val = field[1]

                        if k == "author"
                            val = self.parse_authors(val)
                        end

                        if k == "year"
                            val = Dict("literal" => val)
                            k = "issued"
                        end

                        if k == "pages"
                            val = replace(val, "--" => "-")
                            k = "page"
                        end

                        if k == "title"
                            #   Preserve capitalization, as described in http://tex.stackexchange.com/questions/7288/preserving-capitalization-in-bibtex-titles
                            #   This will likely choke on nested curly-brackets, but that doesn't seem like an ordinary practice.
                            function capitalize(s)
                                return s.group(1)*uppercase(s.group(2))
                            end
                            while val.find('{') > -1
                                caps = (val.find('{'), val.find('}'))
                                val = val.replace(val[caps[0]:caps[1]+1], re.sub("(^|\s)(\S)", capitalize, val[caps[0]+1:caps[1]]).strip())
                            end

                        self.records[ key ][k] = val
                    end
                    self.token == "," || break
                end
                if self.token != "}" && self.token != "@"
                    error("@ missing")
                end
            end
        end
    end
end

function parse_authors(self, authors)
    res = []
    authors = split(authors, "and")
    for author in authors
        _author = split(author, ",")
        family = rstrip(strip(_author[0]))
        rec = Dict("family" => family)
        try
            given = rstrip(strip(_author[1]))
            rec["given"] = given
        catch
        end
        push!(res, rec)
    end
    return res
end

"""Returns json formated records"""
json(self) = JSON.dumps(Dict("items"=>self.records.values()))


text:
println( "\n" )

# ===
# TEST: test.accessor.js
println( "===\nTEST: test.accessor.js\n\n" )
include( "test.accessor.jl" )
println( "\n" )

# ===
# TEST: test.array.js
println( "===\nTEST: test.array.js\n\n" )
include( "test.array.jl" )
println( "\n" )

# ===
# TEST: test.deepset.js
println( "===\nTEST: test.deepset.js\n\n" )
include( "test.deepset.jl" )
println( "\n" )

# ===
# TEST: test.matrix.js
println( "===\nTEST: test.matrix.js\n\n" )
include( "test.matrix.jl" )
println( "\n" )

# ===
# TEST: test.number.js
println( "===\nTEST: test.number.js\n\n" )
include( "test.number.jl" )
println( "\n" )

# ===
# TEST: test.partial.js
println( "===\nTEST: test.partial.js\n\n" )
include( "test.partial.jl" )
println( "\n" )

# ===
# TEST: test.typedarray.js
println( "===\nTEST: test.typedarray.js\n\n" )
include( "test.typedarray.jl" )
println( "\n" )

max_stars_repo_name: distributions-io/normal-mgf

max_issues_repo_name: distributions-io/geometric-mgf

max_forks_repo_name: distributions-io/geometric-mgf


text:
# environments
# ------------
# - images
# - headers
# - links
# - code inline
# - code blocks
# - maths inline
# - maths blocks
# - div blocks
# - lx commands
# - hfun commands
# - tables
# - admonition (is it a thing?)
# - blockquotes
# - escaped blocks
# - def blocks
# - indented blocks
# - lists
# - raw html inclusion
# - raw html block with separating line
# - footnotes
# - citations

F.def_GLOBAL_VARS!()

@testset "close:headers" begin
    h = """
    A
    # H
    B
    """ |> fd2html
    @test h // """

A

H

B
"""
    h = """
    A
    # H1
    ## H2
    B
    ## H2
    C
    """ |> fd2html
    @test h // """

A

H1

H2

B

H2

C
"""
end

@testset "close:div" begin
    h = """
    A
    @@b,c D @@
    E
    """ |> fd2html
    @test h // """

A

D

E
"""
end

# NOTE: will fail upon use of CommonMark (list)
@testset "close:list+i" begin
    h = """
    A
    * B
    * `C`
    E
    """ |> fd2html
    @test h // """

A

B

C

E
"""
end

@testset "ending p" begin
    # Franklin.LOGGING[] = true
    h = raw"""
    @def date_format = "e, d u Y"
    ```julia:ex1
    #hideall
    using Dates
    println(Franklin.fd_date(DateTime("1996-01-01T12:30:00")))
    ```
    A \textoutput{ex1} B
    """ |> fd2html
    @test h // "
A Mon, 1 Jan 1996 B
"
end


text:
"""
    DualDict{K1, K2, V}

A convenient dictionary type that efficiently implements a dictionary with
two different key types `K1` and `K2` for one common value type `V`. General
usage is experimental.

**Example**
```julia-repl
julia> DualDict{String, Int, Float64}() # Initialize empty DualDict
DualDict{String,Int64,Float64} with 0 entries:

julia> dd = DualDict{String, Int, Float64}("a" => 2.1, 3 => 5.7, "v" => 3.1)
DualDict{String,Int64,Float64} with 3 entries:
  "v" => 3.1
  "a" => 2.1
  3 => 5.7

julia> dd[3] # Index normally
5.7

julia> dd["v"] = 42 # Set values normally
42

julia> dd[2] = 8.3 # Add new pairs
8.3

julia> dd
DualDict{String,Int64,Float64} with 4 entries:
  "v" => 42.0
  "a" => 2.1
  2 => 8.3
  3 => 5.7

julia> delete!(dd, "a") # Delete pairs as usual
DualDict{String,Int64,Float64} with 3 entries:
  "v" => 42.0
  2 => 8.3
  3 => 5.7

julia> haskey(dd, "a") # Use haskey
false

julia> [k => v for (k, v) in dd] # Iterate as normal Dict
3-element Array{Pair{A,Float64} where A,1}:
 "v" => 42.0
   2 => 8.3
   3 => 5.7

julia> empty!(dd) # Empty it out
DualDict{String,Int64,Float64} with 0 entries:
```
"""
struct DualDict{K1, K2, V}
    dict1::Dict{K1, V}
    dict2::Dict{K2, V}
end

# Default constructor
function DualDict{K1, K2, V}() where {K1, K2, V}
    return DualDict{K1, K2, V}(Dict{K1, V}(), Dict{K2, V}())
end

# Tuple constructor
function DualDict{K1, K2, V}(tup::NTuple{N, Pair{<:Union{K1, K2}, <:V}}
                             ) where {N, K1, K2, V}
    pairs1 = Iterators.filter(x -> x[1] isa K1, tup)
    pairs2 = Iterators.filter(x -> x[1] isa K2, tup)
    return DualDict{K1, K2, V}(Dict{K1, V}(pairs1), Dict{K2, V}(pairs2))
end

# Pair constructor
function DualDict{K1, K2, V}(pairs::Pair{<:Union{K1, K2}, <:V}...
                             ) where {K1, K2, V}
    return DualDict{K1, K2, V}(pairs)
end

# Extend Base.:(==)
function Base.:(==)(dd1::DualDict, dd2::DualDict)::Bool
    return dd1.dict1 == dd2.dict1 && dd1.dict2 == dd2.dict2
end

# Extend Base.copy
function Base.copy(dd::DualDict{K1, K2, V})::DualDict{K1, K2, V} where {K1, K2, V}
    return DualDict(copy(dd.dict1), copy(dd.dict2))
end

# Extend Base.in when checking the first dict
function Base.in(pair::Pair{K1, V}, dd::DualDict{K1, K2, V})::Bool where {K1, K2, V}
    return pair in dd.dict1
end

# Extend Base.in when checking the second dict
function Base.in(pair::Pair{K2, V}, dd::DualDict{K1, K2, V})::Bool where {K1, K2, V}
    return pair in dd.dict2
end

# Extend Base.haskey when the index is that of dict1
function Base.haskey(dd::DualDict{K1, K2, V}, index::K1)::Bool where {K1, K2, V}
    return haskey(dd.dict1, index)
end

# Extend Base.haskey when the index is that of dict1
function Base.haskey(dd::DualDict{K1, K2, V}, index::K2)::Bool where {K1, K2, V}
    return haskey(dd.dict2, index)
end

# Extend Base.length
Base.length(dd::DualDict)::Int = length(dd.dict1) + length(dd.dict2)

# Extend Base.isempty
Base.isempty(dd::DualDict)::Bool = isempty(dd.dict1) && isempty(dd.dict2)

# Extend Base.getindex when the index is that of dict1
function Base.getindex(dd::DualDict{K1, K2, V}, index::K1)::V where {K1, K2, V}
    haskey(dd, index) || throw(BoundsError("Index $index not found in $dd."))
    return dd.dict1[index]
end

# Extend Base.getindex when the index is that of dict2
function Base.getindex(dd::DualDict{K1, K2, V}, index::K2)::V where {K1, K2, V}
    haskey(dd, index) || throw(BoundsError("Index $index not found in $dd."))
    return dd.dict2[index]
end

# Extend Base.setindex! when the index is that of dict1
function Base.setindex!(dd::DualDict{K1, K2, V}, value,
                        index::K1)::V where {K1, K2, V}
    return dd.dict1[index] = value
end

# Extend Base.setindex! when the index is that of dict2
function Base.setindex!(dd::DualDict{K1, K2, V}, value,
                        index::K2)::V where {K1, K2, V}
    return dd.dict2[index] = value
end

# Extend Base.empty!
function Base.empty!(dd::DualDict{K1, K2, V})::DualDict{K1, K2, V} where {K1, K2, V}
    empty!(dd.dict1)
    empty!(dd.dict2)
    return dd
end

# Extend Base.delete! when the index is that of dict1
function Base.delete!(dd::DualDict{K1, K2, V},
                      index::K1)::DualDict{K1, K2, V} where {K1, K2, V}
    delete!(dd.dict1, index)
    return dd
end

# Extend Base.delete! when the index is that of dict2
function Base.delete!(dd::DualDict{K1, K2, V},
                      index::K2)::DualDict{K1, K2, V} where {K1, K2, V}
    delete!(dd.dict2, index)
    return dd
end

## Helper functions for iterate
# iterate returned nothing
_process_iterate(id::Int, result::Nothing)::Nothing = result

# iterate returned a Tuple
function _process_iterate(id::Int, state::Tuple{Pair{K, V}, Int}
                          )::Tuple{Pair{K, V}, Tuple{Int, Int}} where {K, V}
    return state[1], (id, state[2])
end

# Extend Base.iterate (initial)
function Base.iterate(dd::DualDict{K1, K2, V}
                      )::Union{Nothing, Tuple{Pair{K1,V},Tuple{Int, Int}},
                               Tuple{Pair{K2,V},Tuple{Int, Int}}} where {K1, K2, V}
    if !isempty(dd.dict1)
        return _process_iterate(1, iterate(dd.dict1))
    else
        return _process_iterate(2, iterate(dd.dict2))
    end
end

# Extend Base.iterate (later)
function Base.iterate(dd::DualDict{K1, K2, V},
                      state::Tuple{Int, Int}
                      )::Union{Nothing, Tuple{Pair{K1,V},Tuple{Int, Int}},
                               Tuple{Pair{K2,V},Tuple{Int, Int}}} where {K1, K2, V}
    if state[1] == 1
        result = _process_iterate(1, iterate(dd.dict1, state[2]))
        return result !== nothing ? result : _process_iterate(2, iterate(dd.dict2))
    else
        return _process_iterate(2, iterate(dd.dict2, state[2]))
    end
end

# Extend Base.keys
function Base.keys(dd::DualDict)::Vector
    return [collect(keys(dd.dict1))..., collect(keys(dd.dict2))...]
end

# Extend Base.values
function Base.values(dd::DualDict{K1, K2, V})::Vector{V} where {K1, K2, V}
    return [values(dd.dict1)..., values(dd.dict2)...]
end

# Make the show string
function show_string(dd::DualDict)::String
    num_pairs = length(dd)
    if num_pairs == 1
        str = string(typeof(dd), " with 1 entry:")
    else
        str = string(typeof(dd), " with ", num_pairs, " entries:")
    end
    for pair in dd
        str *= string("\n  ", pair)
    end
    return str
end

# Show DualDicts in REPLMode
function Base.show(io::IO, dd::DualDict)
    print(io, show_string(dd))
    return
end

# Show DualDicts in IJuliaMode
function Base.show(io::IO, ::MIME"text/latex", dd::DualDict)
    print(io, show_string(dd))
end


text:
# Example on how to load types and methods from C++ into an existing module
module ExtendedTypes

include(joinpath(@__DIR__, "testcommon.jl"))

@readmodule CxxWrap.CxxWrapCore.libextended()
@wraptypes
@wrapfunctions

function __init__()
  @initcxx
end

export ExtendedWorld, greet

end

using .ExtendedTypes
using Test

@testset "$(basename(@__FILE__)[1:end-3])" begin

w = ExtendedWorld()
@test greet(w) == "default hello"

end


text:
using StructuralEquationModels, ProximalSEM, ProximalOperators, CSV, DataFrames

# load data
dat = DataFrame(CSV.File("data_dem.csv"))

############################################################################
### define models
############################################################################

graph_1 = """
    ind60 =∼ 1*x1 + x2 + x3
    dem60 =∼ 1*y1 + y2 + y3 + y4
    dem65 =∼ 1*y5 + y6 + y7 + y8
    dem60 ∼ ind60
    dem65 ∼ dem60
    dem65 ∼ ind60
    ind60 ∼∼ ind60
    dem60 ∼∼ dem60
    dem65 ∼∼ dem65
    x1 ∼∼ x1
    x2 ∼∼ x2
    x3 ∼∼ x3
    y1 ∼∼ y1
    y2 ∼∼ y2
    y3 ∼∼ y3
    y4 ∼∼ y4
    y5 ∼∼ y5
    y6 ∼∼ y6
    y7 ∼∼ y7
    y8 ∼∼ y8
    y1 ∼∼ y5
    y2 ∼∼ y4 + y6
    y3 ∼∼ y7
    y4 ∼∼ y8
    y6 ∼∼ y8
"""

observed_vars = ["x1", "x2", "x3", "y1", "y2", "y3", "y4", "y5", "y6", "y7", "y8"]
latent_vars = ["ind60", "dem60", "dem65"]

partable = ParameterTable(latent_vars, observed_vars, graph_1)

# use ridge from StructuralEquationModels
model_ridge = Sem(
    specification = partable,
    data = dat,
    loss = (SemML, SemRidge,),
    α_ridge = .001,
    which_ridge = 16:20
)

solution_ridge = sem_fit(model_ridge)

# use ridge from ProximalSEM
λ = zeros(31); λ[16:20] .= 0.002

model_prox = Sem(
    specification = partable,
    data = dat,
    loss = (SemML,),
    diff = SemDiffProximal,
    operator_g = SqrNormL2(λ)
)

solution_prox = sem_fit(model_prox)

@testset "ridge_solution" begin
    @test isapprox(solution_prox.solution, solution_ridge.solution; rtol = 1e-4)
end


text:
using LinearAlgebra: cross, dot, norm, normalize

@testset "p06e" begin
    eps0, psia, oma, bpa, bqa, pia, bpia, epsa, chia, za, zetaa, thetaa, pa, gam, phi, psi = p06e(2400000.5, 52541.0)
    @test eps0 ≈ 0.4090926006005828715 atol=1e-14
    @test psia ≈ 0.6664369630191613431e-3 atol=1e-14
    @test oma ≈ 0.4090925973783255982 atol=1e-14
    @test bpa ≈ 0.5561149371265209445e-6 atol=1e-14
    @test bqa ≈ -0.6191517193290621270e-5 atol=1e-14
    @test pia ≈ 0.6216441751884382923e-5 atol=1e-14
    @test bpia ≈ 3.052014180023779882 atol=1e-14
    @test epsa ≈ 0.4090864054922431688 atol=1e-14
    @test chia ≈ 0.1387703379530915364e-5 atol=1e-14
    @test za ≈ 0.2921789846651790546e-3 atol=1e-14
    @test zetaa ≈ 0.3178773290332009310e-3 atol=1e-14
    @test thetaa ≈ 0.2650932701657497181e-3 atol=1e-14
    @test pa ≈ 0.6651637681381016344e-3 atol=1e-14
    @test gam ≈ 0.1398077115963754987e-5 atol=1e-14
    @test phi ≈ 0.4090864090837462602 atol=1e-14
    @test psi ≈ 0.6664464807480920325e-3 atol=1e-14
end

@testset "p2pv" begin
    p = [0.25,1.2,3.0]
    pv = ERFA._p2pv(p)
    @test pv[1][1] ≈ 0.25 atol=1e-12
    @test pv[1][2] ≈ 1.2 atol=1e-12
    @test pv[1][3] ≈ 3.0 atol=1e-12
    @test pv[2][1] ≈ 0.0 atol=1e-12
    @test pv[2][2] ≈ 0.0 atol=1e-12
    @test pv[2][3] ≈ 0.0 atol=1e-12
    @test_throws ArgumentError ERFA._p2pv([0.25,1.2])
    @test pv == [p, zeros(3)]
end

@testset "p2s" begin
    theta, phi, r = p2s([100.,-50.,25.])
    @test theta ≈ -0.4636476090008061162 atol=1e-12
    @test phi ≈ 0.2199879773954594463 atol=1e-12
    @test r ≈ 114.5643923738960002 atol=1e-12
    @test_throws ArgumentError p2s([100.,-50.])
end

@testset "pap" begin
    a = [1.,0.1,0.2]
    b = [-3.,1e-3,0.2]
    theta = pap(a, b)
    @test theta ≈ 0.3671514267841113674 atol=1e-12
    @test_throws ArgumentError pap(a[1:2], b)
    @test_throws ArgumentError pap(a, b[1:2])
end

@testset "pas" begin
    p = pas(1.0, 0.1, 0.2, -1.0)
    @test p ≈ -2.724544922932270424 atol=1e-12
end

@testset "pb06" begin
    bzeta, bz, btheta = pb06(2400000.5, 50123.9999)
    @test bzeta ≈ -0.5092634016326478238e-3 atol=1e-12
    @test bz ≈ -0.3602772060566044413e-3 atol=1e-12
    @test btheta ≈ -0.3779735537167811177e-3 atol=1e-12
end

@testset "pfw06" begin
    gamb, phib, psib, epsa = pfw06(2400000.5, 50123.9999)
    @test gamb ≈ -0.2243387670997995690e-5 atol=1e-16
    @test phib ≈ 0.4091014602391312808 atol=1e-12
    @test psib ≈ -0.9501954178013031895e-3 atol=1e-14
    @test epsa ≈ 0.4091014316587367491 atol=1e-12
end

@testset "pdp" begin
    a = [2.,2.,3.]
    b = [1.,3.,4.]
    ab = ERFA._pdp(a, b)
    @test ab ≈ 20 atol=1e-12
    @test ab == dot(a, b)
    @test_throws ArgumentError ERFA._pdp([2.,2.,3.], [1.,3.])
    @test_throws ArgumentError ERFA._pdp([2.,3.], [1.,3.,4.])
end

@testset "plan94" begin
    @test_throws ERFAException plan94(2400000.5, -320000., 10)
    p, v = @test_logs (:warn,) plan94(2400000.5, -320000., 3)
    @test p[1] ≈ 0.9308038666832975759 atol=1e-11
    @test p[2] ≈ 0.3258319040261346000 atol=1e-11
    @test p[3] ≈ 0.1422794544481140560 atol=1e-11
    @test v[1] ≈ -0.6429458958255170006e-2 atol=1e-11
    @test v[2] ≈ 0.1468570657704237764e-1 atol=1e-11
    @test v[3] ≈ 0.6406996426270981189e-2 atol=1e-11

    p, v = plan94(2400000.5, 43999.9, 1)
    @test p[1] ≈ 0.2945293959257430832 atol=1e-11
    @test p[2] ≈ -0.2452204176601049596 atol=1e-11
    @test p[3] ≈ -0.1615427700571978153 atol=1e-11
    @test v[1] ≈ 0.1413867871404614441e-1 atol=1e-11
    @test v[2] ≈ 0.1946548301104706582e-1 atol=1e-11
    @test v[3] ≈ 0.8929809783898904786e-2 atol=1e-11
end

@testset "pm" begin
    p = [0.3,1.2,-2.5]
    m = ERFA._pm(p)
    @test m ≈ 2.789265136196270604 atol=1e-14
    @test_throws ArgumentError ERFA._pm([0.3,1.2])
    @test m ≈ norm(p) atol=1e-14
end

@testset "pmp" begin
    a = [2.0,2.0,3.0]
    b = [1.0,3.0,4.0]
    amb = ERFA._pmp(a, b)
    @test amb[1] ≈ 1.0 atol=1e-12
    @test amb[2] ≈ -1.0 atol=1e-12
    @test amb[3] ≈ -1.0 atol=1e-12
    @test_throws ArgumentError ERFA._pmp(a[1:2], b)
    @test_throws ArgumentError ERFA._pmp(a, b[1:2])
    @test amb == a .- b
end

@testset "pmpx" begin
    rc = 1.234
    dc = 0.789
    pr = 1e-5
    pd = -2e-5
    px = 1e-2
    rv = 10.0
    pmt = 8.75
    vob = [0.9, 0.4, 0.1]
    pco = pmpx(rc, dc, pr, pd, px, rv, pmt, vob)
    @test pco[1] ≈ 0.2328137623960308440 atol=1e-12
    @test pco[2] ≈ 0.6651097085397855317 atol=1e-12
    @test pco[3] ≈ 0.7095257765896359847 atol=1e-12
    @test_throws ArgumentError pmpx(rc, dc, pr, pd, px, rv, pmt, vob[1:2])
end

@testset "pmsafe" begin
    ra1 = 1.234
    dec1 = 0.789
    pmr1 = 1e-5
    pmd1 = -2e-5
    px1 = 1e-2
    rv1 = 10.0
    ep1a = 2400000.5
    ep1b = 48348.5625
    ep2a = 2400000.5
    ep2b = 51544.5
    ra2, dec2, pmr2, pmd2, px2, rv2 = pmsafe(ra1, dec1, pmr1, pmd1, px1, rv1, ep1a, ep1b, ep2a, ep2b)
    @test ra2 ≈ 1.234087484501017061 atol=1e-12
    @test dec2 ≈ 0.7888249982450468567 atol=1e-12
    @test pmr2 ≈ 0.9996457663586073988e-5 atol=1e-12
    @test pmd2 ≈ -0.2000040085106754565e-4 atol=1e-16
    @test px2 ≈ 0.9999997295356830666e-2 atol=1e-12
    @test rv2 ≈ 10.38468380293920069 atol=1e-10
end

@testset "pmat00" begin
    rbp = pmat00(2400000.5, 50123.9999)
    @test rbp[1,1] ≈ 0.9999995505175087260 atol=1e-12
    @test rbp[1,2] ≈ 0.8695405883617884705e-3 atol=1e-14
    @test rbp[1,3] ≈ 0.3779734722239007105e-3 atol=1e-14
    @test rbp[2,1] ≈ -0.8695405990410863719e-3 atol=1e-14
    @test rbp[2,2] ≈ 0.9999996219494925900 atol=1e-12
    @test rbp[2,3] ≈ -0.1360775820404982209e-6 atol=1e-14
    @test rbp[3,1] ≈ -0.3779734476558184991e-3 atol=1e-14
    @test rbp[3,2] ≈ -0.1925857585832024058e-6 atol=1e-14
    @test rbp[3,3] ≈ 0.9999999285680153377 atol=1e-12
end

@testset "pmat06" begin
    rbp = pmat06(2400000.5, 50123.9999)
    @test rbp[1,1] ≈ 0.9999995505176007047 atol=1e-12
    @test rbp[1,2] ≈ 0.8695404617348208406e-3 atol=1e-14
    @test rbp[1,3] ≈ 0.3779735201865589104e-3 atol=1e-14
    @test rbp[2,1] ≈ -0.8695404723772031414e-3 atol=1e-14
    @test rbp[2,2] ≈ 0.9999996219496027161 atol=1e-12
    @test rbp[2,3] ≈ -0.1361752497080270143e-6 atol=1e-14
    @test rbp[3,1] ≈ -0.3779734957034089490e-3 atol=1e-14
    @test rbp[3,2] ≈ -0.1924880847894457113e-6 atol=1e-14
    @test rbp[3,3] ≈ 0.9999999285679971958 atol=1e-12
end

@testset "pmat76" begin
    rmatp = pmat76(2400000.5, 50123.9999)
    @test rmatp[1,1] ≈ 0.9999995504328350733 atol=1e-12
    @test rmatp[1,2] ≈ 0.8696632209480960785e-3 atol=1e-14
    @test rmatp[1,3] ≈ 0.3779153474959888345e-3 atol=1e-14
    @test rmatp[2,1] ≈ -0.8696632209485112192e-3 atol=1e-14
    @test rmatp[2,2] ≈ 0.9999996218428560614 atol=1e-12
    @test rmatp[2,3] ≈ -0.1643284776111886407e-6 atol=1e-14
    @test rmatp[3,1] ≈ -0.3779153474950335077e-3 atol=1e-14
    @test rmatp[3,2] ≈ -0.1643306746147366896e-6 atol=1e-14
    @test rmatp[3,3] ≈ 0.9999999285899790119 atol=1e-12
end

@testset "pn" begin
    p = [0.3,1.2,-2.5]
    r, u = ERFA._pn(p)
    @test r ≈ 2.789265136196270604 atol=1e-12
    @test u[1] ≈ 0.1075552109073112058 atol=1e-12
    @test u[2] ≈ 0.4302208436292448232 atol=1e-12
    @test u[3] ≈ -0.8962934242275933816 atol=1e-12
    @test_throws ArgumentError ERFA._pn([0.3,1.2])
    r1, u1 = norm(p), normalize(p)
    @test r ≈ r1 atol=1e-12
    @test u ≈ u1 atol=1e-12
end

@testset "pn00" begin
    dpsi = -0.9632552291149335877e-5
    deps =  0.4063197106621141414e-4
    epsa, rb, rp, rbp, rn, rbpn = pn00(2400000.5, 53736.0, dpsi, deps)
    @test epsa ≈ 0.4090791789404229916 atol=1e-12
    @test rb[1,1] ≈ 0.9999999999999942498 atol=1e-12
    @test rb[1,2] ≈ -0.7078279744199196626e-7 atol=1e-18
    @test rb[1,3] ≈ 0.8056217146976134152e-7 atol=1e-18
    @test rb[2,1] ≈ 0.7078279477857337206e-7 atol=1e-18
    @test rb[2,2] ≈ 0.9999999999999969484 atol=1e-12
    @test rb[2,3] ≈ 0.3306041454222136517e-7 atol=1e-18
    @test rb[3,1] ≈ -0.8056217380986972157e-7 atol=1e-18
    @test rb[3,2] ≈ -0.3306040883980552500e-7 atol=1e-18
    @test rb[3,3] ≈ 0.9999999999999962084 atol=1e-12
    @test rp[1,1] ≈ 0.9999989300532289018 atol=1e-12
    @test rp[1,2] ≈ -0.1341647226791824349e-2 atol=1e-14
    @test rp[1,3] ≈ -0.5829880927190296547e-3 atol=1e-14
    @test rp[2,1] ≈ 0.1341647231069759008e-2 atol=1e-14
    @test rp[2,2] ≈ 0.9999990999908750433 atol=1e-12
    @test rp[2,3] ≈ -0.3837444441583715468e-6 atol=1e-14
    @test rp[3,1] ≈ 0.5829880828740957684e-3 atol=1e-14
    @test rp[3,2] ≈ -0.3984203267708834759e-6 atol=1e-14
    @test rp[3,3] ≈ 0.9999998300623538046 atol=1e-12
    @test rbp[1,1] ≈ 0.9999989300052243993 atol=1e-12
    @test rbp[1,2] ≈ -0.1341717990239703727e-2 atol=1e-14
    @test rbp[1,3] ≈ -0.5829075749891684053e-3 atol=1e-14
    @test rbp[2,1] ≈ 0.1341718013831739992e-2 atol=1e-14
    @test rbp[2,2] ≈ 0.9999990998959191343 atol=1e-12
    @test rbp[2,3] ≈ -0.3505759733565421170e-6 atol=1e-14
    @test rbp[3,1] ≈ 0.5829075206857717883e-3 atol=1e-14
    @test rbp[3,2] ≈ -0.4315219955198608970e-6 atol=1e-14
    @test rbp[3,3] ≈ 0.9999998301093036269 atol=1e-12
    @test rn[1,1] ≈ 0.9999999999536069682 atol=1e-12
    @test rn[1,2] ≈ 0.8837746144872140812e-5 atol=1e-16
    @test rn[1,3] ≈ 0.3831488838252590008e-5 atol=1e-16
    @test rn[2,1] ≈ -0.8837590456633197506e-5 atol=1e-16
    @test rn[2,2] ≈ 0.9999999991354692733 atol=1e-12
    @test rn[2,3] ≈ -0.4063198798559573702e-4 atol=1e-16
    @test rn[3,1] ≈ -0.3831847930135328368e-5 atol=1e-16
    @test rn[3,2] ≈ 0.4063195412258150427e-4 atol=1e-16
    @test rn[3,3] ≈ 0.9999999991671806225 atol=1e-12
    @test rbpn[1,1] ≈ 0.9999989440499982806 atol=1e-12
    @test rbpn[1,2] ≈ -0.1332880253640848301e-2 atol=1e-14
    @test rbpn[1,3] ≈ -0.5790760898731087295e-3 atol=1e-14
    @test rbpn[2,1] ≈ 0.1332856746979948745e-2 atol=1e-14
    @test rbpn[2,2] ≈ 0.9999991109064768883 atol=1e-12
    @test rbpn[2,3] ≈ -0.4097740555723063806e-4 atol=1e-14
    @test rbpn[3,1] ≈ 0.5791301929950205000e-3 atol=1e-14
    @test rbpn[3,2] ≈ 0.4020553681373702931e-4 atol=1e-14
    @test rbpn[3,3] ≈ 0.9999998314958529887 atol=1e-12
end

@testset "pn00a" begin
    dpsi, deps, epsa, rb, rp, rbp, rn, rbpn = pn00a(2400000.5, 53736.0)
    @test dpsi ≈ -0.9630909107115518431e-5 atol=1e-12
    @test deps ≈ 0.4063239174001678710e-4 atol=1e-12
    @test epsa ≈ 0.4090791789404229916 atol=1e-12
    @test rb[1,1] ≈ 0.9999999999999942498 atol=1e-12
    @test rb[1,2] ≈ -0.7078279744199196626e-7 atol=1e-16
    @test rb[1,3] ≈ 0.8056217146976134152e-7 atol=1e-16
    @test rb[2,1] ≈ 0.7078279477857337206e-7 atol=1e-16
    @test rb[2,2] ≈ 0.9999999999999969484 atol=1e-12
    @test rb[2,3] ≈ 0.3306041454222136517e-7 atol=1e-16
    @test rb[3,1] ≈ -0.8056217380986972157e-7 atol=1e-16
    @test rb[3,2] ≈ -0.3306040883980552500e-7 atol=1e-16
    @test rb[3,3] ≈ 0.9999999999999962084 atol=1e-12
    @test rp[1,1] ≈ 0.9999989300532289018 atol=1e-12
    @test rp[1,2] ≈ -0.1341647226791824349e-2 atol=1e-14
    @test rp[1,3] ≈ -0.5829880927190296547e-3 atol=1e-14
    @test rp[2,1] ≈ 0.1341647231069759008e-2 atol=1e-14
    @test rp[2,2] ≈ 0.9999990999908750433 atol=1e-12
    @test rp[2,3] ≈ -0.3837444441583715468e-6 atol=1e-14
    @test rp[3,1] ≈ 0.5829880828740957684e-3 atol=1e-14
    @test rp[3,2] ≈ -0.3984203267708834759e-6 atol=1e-14
    @test rp[3,3] ≈ 0.9999998300623538046 atol=1e-12
    @test rbp[1,1] ≈ 0.9999989300052243993 atol=1e-12
    @test rbp[1,2] ≈ -0.1341717990239703727e-2 atol=1e-14
    @test rbp[1,3] ≈ -0.5829075749891684053e-3 atol=1e-14
    @test rbp[2,1] ≈ 0.1341718013831739992e-2 atol=1e-14
    @test rbp[2,2] ≈ 0.9999990998959191343 atol=1e-12
    @test rbp[2,3] ≈ -0.3505759733565421170e-6 atol=1e-14
    @test rbp[3,1] ≈ 0.5829075206857717883e-3 atol=1e-14
    @test rbp[3,2] ≈ -0.4315219955198608970e-6 atol=1e-14
    @test rbp[3,3] ≈ 0.9999998301093036269 atol=1e-12
    @test rn[1,1] ≈ 0.9999999999536227949 atol=1e-12
    @test rn[1,2] ≈ 0.8836238544090873336e-5 atol=1e-14
    @test rn[1,3] ≈ 0.3830835237722400669e-5 atol=1e-14
    @test rn[2,1] ≈ -0.8836082880798569274e-5 atol=1e-14
    @test rn[2,2] ≈ 0.9999999991354655028 atol=1e-12
    @test rn[2,3] ≈ -0.4063240865362499850e-4 atol=1e-14
    @test rn[3,1] ≈ -0.3831194272065995866e-5 atol=1e-14
    @test rn[3,2] ≈ 0.4063237480216291775e-4 atol=1e-14
    @test rn[3,3] ≈ 0.9999999991671660338 atol=1e-12
    @test rbpn[1,1] ≈ 0.9999989440476103435 atol=1e-12
    @test rbpn[1,2] ≈ -0.1332881761240011763e-2 atol=1e-14
    @test rbpn[1,3] ≈ -0.5790767434730085751e-3 atol=1e-14
    @test rbpn[2,1] ≈ 0.1332858254308954658e-2 atol=1e-14
    @test rbpn[2,2] ≈ 0.9999991109044505577 atol=1e-12
    @test rbpn[2,3] ≈ -0.4097782710396580452e-4 atol=1e-14
    @test rbpn[3,1] ≈ 0.5791308472168152904e-3 atol=1e-14
    @test rbpn[3,2] ≈ 0.4020595661591500259e-4 atol=1e-14
    @test rbpn[3,3] ≈ 0.9999998314954572304 atol=1e-12
end

@testset "pn06" begin
    dpsi = -0.9632552291149335877e-5
    deps =  0.4063197106621141414e-4
    epsa, rb, rp, rbp, rn, rbpn = pn06(2400000.5, 53736.0, dpsi, deps)
    @test epsa ≈ 0.4090789763356509926 atol=1e-12
    @test rb[1,1] ≈ 0.9999999999999942497 atol=1e-12
    @test rb[1,2] ≈ -0.7078368960971557145e-7 atol=1e-14
    @test rb[1,3] ≈ 0.8056213977613185606e-7 atol=1e-14
    @test rb[2,1] ≈ 0.7078368694637674333e-7 atol=1e-14
    @test rb[2,2] ≈ 0.9999999999999969484 atol=1e-12
    @test rb[2,3] ≈ 0.3305943742989134124e-7 atol=1e-14
    @test rb[3,1] ≈ -0.8056214211620056792e-7 atol=1e-14
    @test rb[3,2] ≈ -0.3305943172740586950e-7 atol=1e-14
    @test rb[3,3] ≈ 0.9999999999999962084 atol=1e-12
    @test rp[1,1] ≈ 0.9999989300536854831 atol=1e-12
    @test rp[1,2] ≈ -0.1341646886204443795e-2 atol=1e-14
    @test rp[1,3] ≈ -0.5829880933488627759e-3 atol=1e-14
    @test rp[2,1] ≈ 0.1341646890569782183e-2 atol=1e-14
    @test rp[2,2] ≈ 0.9999990999913319321 atol=1e-12
    @test rp[2,3] ≈ -0.3835944216374477457e-6 atol=1e-14
    @test rp[3,1] ≈ 0.5829880833027867368e-3 atol=1e-14
    @test rp[3,2] ≈ -0.3985701514686976112e-6 atol=1e-14
    @test rp[3,3] ≈ 0.9999998300623534950 atol=1e-12
    @test rbp[1,1] ≈ 0.9999989300056797893 atol=1e-12
    @test rbp[1,2] ≈ -0.1341717650545059598e-2 atol=1e-14
    @test rbp[1,3] ≈ -0.5829075756493728856e-3 atol=1e-14
    @test rbp[2,1] ≈ 0.1341717674223918101e-2 atol=1e-14
    @test rbp[2,2] ≈ 0.9999990998963748448 atol=1e-12
    @test rbp[2,3] ≈ -0.3504269280170069029e-6 atol=1e-14
    @test rbp[3,1] ≈ 0.5829075211461454599e-3 atol=1e-14
    @test rbp[3,2] ≈ -0.4316708436255949093e-6 atol=1e-14
    @test rbp[3,3] ≈ 0.9999998301093032943 atol=1e-12
    @test rn[1,1] ≈ 0.9999999999536069682 atol=1e-12
    @test rn[1,2] ≈ 0.8837746921149881914e-5 atol=1e-14
    @test rn[1,3] ≈ 0.3831487047682968703e-5 atol=1e-14
    @test rn[2,1] ≈ -0.8837591232983692340e-5 atol=1e-14
    @test rn[2,2] ≈ 0.9999999991354692664 atol=1e-12
    @test rn[2,3] ≈ -0.4063198798558931215e-4 atol=1e-14
    @test rn[3,1] ≈ -0.3831846139597250235e-5 atol=1e-14
    @test rn[3,2] ≈ 0.4063195412258792914e-4 atol=1e-14
    @test rn[3,3] ≈ 0.9999999991671806293 atol=1e-12
    @test rbpn[1,1] ≈ 0.9999989440504506688 atol=1e-12
    @test rbpn[1,2] ≈ -0.1332879913170492655e-2 atol=1e-14
    @test rbpn[1,3] ≈ -0.5790760923225655753e-3 atol=1e-14
    @test rbpn[2,1] ≈ 0.1332856406595754748e-2 atol=1e-14
    @test rbpn[2,2] ≈ 0.9999991109069366795 atol=1e-12
    @test rbpn[2,3] ≈ -0.4097725651142641812e-4 atol=1e-14
    @test rbpn[3,1] ≈ 0.5791301952321296716e-3 atol=1e-14
    @test rbpn[3,2] ≈ 0.4020538796195230577e-4 atol=1e-14
    @test rbpn[3,3] ≈ 0.9999998314958576778 atol=1e-12
end

@testset "pn06a" begin
    dpsi, deps, epsa, rb, rp, rbp, rn, rbpn = pn06a(2400000.5, 53736.0)
    @test dpsi ≈ -0.9630912025820308797e-5 atol=1e-12
    @test deps ≈ 0.4063238496887249798e-4 atol=1e-12
    @test epsa ≈ 0.4090789763356509926 atol=1e-12
    @test rb[1,1] ≈ 0.9999999999999942497 atol=1e-12
    @test rb[1,2] ≈ -0.7078368960971557145e-7 atol=1e-14
    @test rb[1,3] ≈ 0.8056213977613185606e-7 atol=1e-14
    @test rb[2,1] ≈ 0.7078368694637674333e-7 atol=1e-14
    @test rb[2,2] ≈ 0.9999999999999969484 atol=1e-12
    @test rb[2,3] ≈ 0.3305943742989134124e-7 atol=1e-14
    @test rb[3,1] ≈ -0.8056214211620056792e-7 atol=1e-14
    @test rb[3,2] ≈ -0.3305943172740586950e-7 atol=1e-14
    @test rb[3,3] ≈ 0.9999999999999962084 atol=1e-12
    @test rp[1,1] ≈ 0.9999989300536854831 atol=1e-12
    @test rp[1,2] ≈ -0.1341646886204443795e-2 atol=1e-14
    @test rp[1,3] ≈ -0.5829880933488627759e-3 atol=1e-14
    @test rp[2,1] ≈ 0.1341646890569782183e-2 atol=1e-14
    @test rp[2,2] ≈ 0.9999990999913319321 atol=1e-12
    @test rp[2,3] ≈ -0.3835944216374477457e-6 atol=1e-14
    @test rp[3,1] ≈ 0.5829880833027867368e-3 atol=1e-14
    @test rp[3,2] ≈ -0.3985701514686976112e-6 atol=1e-14
    @test rp[3,3] ≈ 0.9999998300623534950 atol=1e-12
    @test rbp[1,1] ≈ 0.9999989300056797893 atol=1e-12
    @test rbp[1,2] ≈ -0.1341717650545059598e-2 atol=1e-14
    @test rbp[1,3] ≈ -0.5829075756493728856e-3 atol=1e-14
    @test rbp[2,1] ≈ 0.1341717674223918101e-2 atol=1e-14
    @test rbp[2,2] ≈ 0.9999990998963748448 atol=1e-12
    @test rbp[2,3] ≈ -0.3504269280170069029e-6 atol=1e-14
    @test rbp[3,1] ≈ 0.5829075211461454599e-3 atol=1e-14
    @test rbp[3,2] ≈ -0.4316708436255949093e-6 atol=1e-14
    @test rbp[3,3] ≈ 0.9999998301093032943 atol=1e-12
    @test rn[1,1] ≈ 0.9999999999536227668 atol=1e-12
    @test rn[1,2] ≈ 0.8836241998111535233e-5 atol=1e-14
    @test rn[1,3] ≈ 0.3830834608415287707e-5 atol=1e-14
    @test rn[2,1] ≈ -0.8836086334870740138e-5 atol=1e-14
    @test rn[2,2] ≈ 0.9999999991354657474 atol=1e-12
    @test rn[2,3] ≈ -0.4063240188248455065e-4 atol=1e-14
    @test rn[3,1] ≈ -0.3831193642839398128e-5 atol=1e-14
    @test rn[3,2] ≈ 0.4063236803101479770e-4 atol=1e-14
    @test rn[3,3] ≈ 0.9999999991671663114 atol=1e-12
    @test rbpn[1,1] ≈ 0.9999989440480669738 atol=1e-12
    @test rbpn[1,2] ≈ -0.1332881418091915973e-2 atol=1e-14
    @test rbpn[1,3] ≈ -0.5790767447612042565e-3 atol=1e-14
    @test rbpn[2,1] ≈ 0.1332857911250989133e-2 atol=1e-14
    @test rbpn[2,2] ≈ 0.9999991109049141908 atol=1e-12
    @test rbpn[2,3] ≈ -0.4097767128546784878e-4 atol=1e-14
    @test rbpn[3,1] ≈ 0.5791308482835292617e-3 atol=1e-14
    @test rbpn[3,2] ≈ 0.4020580099454020310e-4 atol=1e-14
    @test rbpn[3,3] ≈ 0.9999998314954628695 atol=1e-12
end

@testset "pnm00a" begin
    rbpn = pnm00a(2400000.5, 50123.9999)
    @test rbpn[1,1] ≈ 0.9999995832793134257 atol=1e-12
    @test rbpn[1,2] ≈ 0.8372384254137809439e-3 atol=1e-14
    @test rbpn[1,3] ≈ 0.3639684306407150645e-3 atol=1e-14
    @test rbpn[2,1] ≈ -0.8372535226570394543e-3 atol=1e-14
    @test rbpn[2,2] ≈ 0.9999996486491582471 atol=1e-12
    @test rbpn[2,3] ≈ 0.4132915262664072381e-4 atol=1e-14
    @test rbpn[3,1] ≈ -0.3639337004054317729e-3 atol=1e-14
    @test rbpn[3,2] ≈ -0.4163386925461775873e-4 atol=1e-14
    @test rbpn[3,3] ≈ 0.9999999329094390695 atol=1e-12
end

@testset "pn00b" begin
    dpsi, deps, epsa, rb, rp, rbp, rn, rbpn = pn00b(2400000.5, 53736.0)
    @test dpsi ≈ -0.9632552291148362783e-5 atol=1e-12
    @test deps ≈ 0.4063197106621159367e-4 atol=1e-12
    @test epsa ≈ 0.4090791789404229916 atol=1e-12
    @test rb[1,1] ≈ 0.9999999999999942498 atol=1e-12
    @test rb[1,2] ≈ -0.7078279744199196626e-7 atol=1e-16
    @test rb[1,3] ≈ 0.8056217146976134152e-7 atol=1e-16
    @test rb[2,1] ≈ 0.7078279477857337206e-7 atol=1e-16
    @test rb[2,2] ≈ 0.9999999999999969484 atol=1e-12
    @test rb[2,3] ≈ 0.3306041454222136517e-7 atol=1e-16
    @test rb[3,1] ≈ -0.8056217380986972157e-7 atol=1e-16
    @test rb[3,2] ≈ -0.3306040883980552500e-7 atol=1e-16
    @test rb[3,3] ≈ 0.9999999999999962084 atol=1e-12
    @test rp[1,1] ≈ 0.9999989300532289018 atol=1e-12
    @test rp[1,2] ≈ -0.1341647226791824349e-2 atol=1e-14
    @test rp[1,3] ≈ -0.5829880927190296547e-3 atol=1e-14
    @test rp[2,1] ≈ 0.1341647231069759008e-2 atol=1e-14
    @test rp[2,2] ≈ 0.9999990999908750433 atol=1e-12
    @test rp[2,3] ≈ -0.3837444441583715468e-6 atol=1e-14
    @test rp[3,1] ≈ 0.5829880828740957684e-3 atol=1e-14
    @test rp[3,2] ≈ -0.3984203267708834759e-6 atol=1e-14
    @test rp[3,3] ≈ 0.9999998300623538046 atol=1e-12
    @test rbp[1,1] ≈ 0.9999989300052243993 atol=1e-12
    @test rbp[1,2] ≈ -0.1341717990239703727e-2 atol=1e-14
    @test rbp[1,3] ≈ -0.5829075749891684053e-3 atol=1e-14
    @test rbp[2,1] ≈ 0.1341718013831739992e-2 atol=1e-14
    @test rbp[2,2] ≈ 0.9999990998959191343 atol=1e-12
    @test rbp[2,3] ≈ -0.3505759733565421170e-6 atol=1e-14
    @test rbp[3,1] ≈ 0.5829075206857717883e-3 atol=1e-14
    @test rbp[3,2] ≈ -0.4315219955198608970e-6 atol=1e-14
    @test rbp[3,3] ≈ 0.9999998301093036269 atol=1e-12
    @test rn[1,1] ≈ 0.9999999999536069682 atol=1e-12
    @test rn[1,2] ≈ 0.8837746144871248011e-5 atol=1e-14
    @test rn[1,3] ≈ 0.3831488838252202945e-5 atol=1e-14
    @test rn[2,1] ≈ -0.8837590456632304720e-5 atol=1e-14
    @test rn[2,2] ≈ 0.9999999991354692733 atol=1e-12
    @test rn[2,3] ≈ -0.4063198798559591654e-4 atol=1e-14
    @test rn[3,1] ≈ -0.3831847930134941271e-5 atol=1e-14
    @test rn[3,2] ≈ 0.4063195412258168380e-4 atol=1e-14
    @test rn[3,3] ≈ 0.9999999991671806225 atol=1e-12
    @test rbpn[1,1] ≈ 0.9999989440499982806 atol=1e-12
    @test rbpn[1,2] ≈ -0.1332880253640849194e-2 atol=1e-14
    @test rbpn[1,3] ≈ -0.5790760898731091166e-3 atol=1e-14
    @test rbpn[2,1] ≈ 0.1332856746979949638e-2 atol=1e-14
    @test rbpn[2,2] ≈ 0.9999991109064768883 atol=1e-12
    @test rbpn[2,3] ≈ -0.4097740555723081811e-4 atol=1e-14
    @test rbpn[3,1] ≈ 0.5791301929950208873e-3 atol=1e-14
    @test rbpn[3,2] ≈ 0.4020553681373720832e-4 atol=1e-14
    @test rbpn[3,3] ≈ 0.9999998314958529887 atol=1e-12
end

@testset "pnm00b" begin
    rbpn = pnm00b(2400000.5, 50123.9999)
    @test rbpn[1,1] ≈ 0.9999995832776208280 atol=1e-12
    @test rbpn[1,2] ≈ 0.8372401264429654837e-3 atol=1e-14
    @test rbpn[1,3] ≈ 0.3639691681450271771e-3 atol=1e-14
    @test rbpn[2,1] ≈ -0.8372552234147137424e-3 atol=1e-14
    @test rbpn[2,2] ≈ 0.9999996486477686123 atol=1e-12
    @test rbpn[2,3] ≈ 0.4132832190946052890e-4 atol=1e-14
    @test rbpn[3,1] ≈ -0.3639344385341866407e-3 atol=1e-14
    @test rbpn[3,2] ≈ -0.4163303977421522785e-4 atol=1e-14
    @test rbpn[3,3] ≈ 0.9999999329092049734 atol=1e-12
end

@testset "pnm06a" begin
    rbpn = pnm06a(2400000.5, 50123.9999)
    @test rbpn[1,1] ≈ 0.9999995832794205484 atol=1e-12
    @test rbpn[1,2] ≈ 0.8372382772630962111e-3 atol=1e-14
    @test rbpn[1,3] ≈ 0.3639684771140623099e-3 atol=1e-14
    @test rbpn[2,1] ≈ -0.8372533744743683605e-3 atol=1e-14
    @test rbpn[2,2] ≈ 0.9999996486492861646 atol=1e-12
    @test rbpn[2,3] ≈ 0.4132905944611019498e-4 atol=1e-14
    @test rbpn[3,1] ≈ -0.3639337469629464969e-3 atol=1e-14
    @test rbpn[3,2] ≈ -0.4163377605910663999e-4 atol=1e-14
    @test rbpn[3,3] ≈ 0.9999999329094260057 atol=1e-12
end

@testset "pnm80" begin
    rmatpn = pnm80(2400000.5, 50123.9999)
    @test rmatpn[1,1] ≈ 0.9999995831934611169 atol=1e-12
    @test rmatpn[1,2] ≈ 0.8373654045728124011e-3 atol=1e-14
    @test rmatpn[1,3] ≈ 0.3639121916933106191e-3 atol=1e-14
    @test rmatpn[2,1] ≈ -0.8373804896118301316e-3 atol=1e-14
    @test rmatpn[2,2] ≈ 0.9999996485439674092 atol=1e-12
    @test rmatpn[2,3] ≈ 0.4130202510421549752e-4 atol=1e-14
    @test rmatpn[3,1] ≈ -0.3638774789072144473e-3 atol=1e-14
    @test rmatpn[3,2] ≈ -0.4160674085851722359e-4 atol=1e-14
    @test rmatpn[3,3] ≈ 0.9999999329310274805 atol=1e-12
end

@testset "pom00" begin
    xp =  2.55060238e-7
    yp =  1.860359247e-6
    sp = -0.1367174580728891460e-10
    rpom = pom00(xp, yp, sp)
    @test rpom[1,1] ≈ 0.9999999999999674721 atol=1e-12
    @test rpom[1,2] ≈ -0.1367174580728846989e-10 atol=1e-16
    @test rpom[1,3] ≈ 0.2550602379999972345e-6 atol=1e-16
    @test rpom[2,1] ≈ 0.1414624947957029801e-10 atol=1e-16
    @test rpom[2,2] ≈ 0.9999999999982695317 atol=1e-12
    @test rpom[2,3] ≈ -0.1860359246998866389e-5 atol=1e-16
    @test rpom[3,1] ≈ -0.2550602379741215021e-6 atol=1e-16
    @test rpom[3,2] ≈ 0.1860359247002414021e-5 atol=1e-16
    @test rpom[3,3] ≈ 0.9999999999982370039 atol=1e-12
end

@testset "ppp" begin
    a = [2.0,2.0,3.0]
    b = [1.0,3.0,4.0]
    apb = ERFA._ppp(a, b)
    @test apb[1] ≈ 3.0 atol=1e-12
    @test apb[2] ≈ 5.0 atol=1e-12
    @test apb[3] ≈ 7.0 atol=1e-12
    @test_throws ArgumentError ERFA._ppp([2.0,3.0], [1.0,3.0,4.0])
    @test_throws ArgumentError ERFA._ppp([2.0,2.0,3.0], [3.0,4.0])
    @test apb == a .+ b
end

@testset "ppsp" begin
    a = [2.0,2.0,3.0]
    s = 5.0
    b = [1.0,3.0,4.0]
    apsb = ERFA._ppsp(a, s, b)
    @test apsb[1] ≈ 7.0 atol=1e-12
    @test apsb[2] ≈ 17.0 atol=1e-12
    @test apsb[3] ≈ 23.0 atol=1e-12
    @test_throws ArgumentError ERFA._ppsp([2.0,2.0], 5.0, [1.0,3.0,4.0])
    @test_throws ArgumentError ERFA._ppsp([2.0,2.0,0.0], 5.0, [3.0,4.0])
    @test apsb == a .+ s .* b
end

@testset "pr00" begin
    dpsipr, depspr = pr00(2400000.5, 53736.)
    @test dpsipr ≈ -0.8716465172668347629e-7 atol=1e-22
    @test depspr ≈ -0.7342018386722813087e-8 atol=1e-22
end

@testset "prec76" begin
    ep01 = 2400000.5
    ep02 = 33282.0
    ep11 = 2400000.5
    ep12 = 51544.0
    zeta, z, theta = prec76(ep01, ep02, ep11, ep12)
    @test zeta ≈ 0.5588961642000161243e-2 atol=1e-12
    @test z ≈ 0.5589922365870680624e-2 atol=1e-12
    @test theta ≈ 0.4858945471687296760e-2 atol=1e-12
end

@testset "pv2p" begin
    pv = [[0.3,1.2,-2.5],[-0.5,3.1,0.9]]
    p = ERFA._pv2p(pv)
    @test p[1] ≈ 0.3 atol=0.0
    @test p[2] ≈ 1.2 atol=0.0
    @test p[3] ≈ -2.5 atol=0.0
    @test_throws ArgumentError ERFA._pv2p([[0.3,1.2],[-0.5,3.1,0.9]])
    @test p == first(pv)
end

@testset "pv2s" begin
    pv = [[-0.4514964673880165,0.03093394277342585,0.05594668105108779],
          [1.292270850663260e-5,2.652814182060692e-6,2.568431853930293e-6]]
    theta, phi, r, td, pd, rd = pv2s(pv)
    @test theta ≈ 3.073185307179586515 atol=1e-12
    @test phi ≈ 0.1229999999999999992 atol=1e-12
    @test r ≈ 0.4559999999999999757 atol=1e-12
    @test td ≈ -0.7800000000000000364e-5 atol=1e-16
    @test pd ≈ 0.9010000000000001639e-5 atol=1e-16
    @test rd ≈ -0.1229999999999999832e-4 atol=1e-16
    pve = [[0.03093394277342585,0.05594668105108779],
          [1.292270850663260e-5,2.652814182060692e-6,2.568431853930293e-6]]
    @test_throws ArgumentError pv2s(pve)
end

@testset "pvdpv" begin
    a = [[2.,2.,3.],[6.,0.,4.]]
    b = [[1.,3.,4.],[0.,2.,8.]]
    adb = pvdpv(a, b)
    @test adb[1] ≈ 20.0 atol=1e-12
    @test adb[2] ≈ 50.0 atol=1e-12
    ae = [[2.,3.],[6.,0.,4.]]
    be = [[3.,4.],[0.,2.,8.]]
    @test_throws ArgumentError pvdpv(ae, b)
    @test_throws ArgumentError pvdpv(a, be)
end

@testset "pvm" begin
    pv = [[0.3,1.2,-2.5],[0.45,-0.25,1.1]]
    r, s = ERFA._pvm(pv)
    @test r ≈ 2.789265136196270604 atol=1e-12
    @test s ≈ 1.214495780149111922 atol=1e-12
    pve = [[1.2,-2.5],[0.45,-0.25,1.1]]
    @test_throws ArgumentError ERFA._pvm(pve)
    @test [r, s] == norm.(pv)
end

@testset "pvmpv" begin
    a = [[2.0,2.0,3.0],[5.0,6.0,3.0]]
    b = [[1.0,3.0,4.0],[3.0,2.0,1.0]]
    amb = ERFA._pvmpv(a, b)
    @test amb[1][1] ≈ 1.0 atol=1e-12
    @test amb[1][2] ≈ -1.0 atol=1e-12
    @test amb[1][3] ≈ -1.0 atol=1e-12
    @test amb[2][1] ≈ 2.0 atol=1e-12
    @test amb[2][2] ≈ 4.0 atol=1e-12
    @test amb[2][3] ≈ 2.0 atol=1e-12
    ae = [[2.0,3.0],[5.0,6.0,3.0]]
    @test_throws ArgumentError ERFA._pvmpv(ae, b)
    @test_throws ArgumentError ERFA._pvmpv(a, ae)
    @test amb == a .- b
end

@testset "pvppv" begin
    a = [[2.0,2.0,3.0],[5.0,6.0,3.0]]
    b = [[1.0,3.0,4.0],[3.0,2.0,1.0]]
    apb = ERFA._pvppv(a, b)
    @test apb[1][1] ≈ 3.0 atol=1e-12
    @test apb[1][2] ≈ 5.0 atol=1e-12
    @test apb[1][3] ≈ 7.0 atol=1e-12
    @test apb[2][1] ≈ 8.0 atol=1e-12
    @test apb[2][2] ≈ 8.0 atol=1e-12
    @test apb[2][3] ≈ 4.0 atol=1e-12
    ae = [[2.0,3.0],[5.0,6.0,3.0]]
    @test_throws ArgumentError ERFA._pvppv(ae, b)
    @test_throws ArgumentError ERFA._pvppv(a, ae)
    @test apb == a .+ b
end

@testset "pvxpv" begin
    a = [[2.0,2.0,3.0],[6.0,0.0,4.0]]
    b = [[1.0,3.0,4.0],[0.0,2.0,8.0]]
    axb = pvxpv(a, b)
    @test axb[1][1] ≈ -1.0 atol=1e-12
    @test axb[1][2] ≈ -5.0 atol=1e-12
    @test axb[1][3] ≈ 4.0 atol=1e-12
    @test axb[2][1] ≈ -2.0 atol=1e-12
    @test axb[2][2] ≈ -36.0 atol=1e-12
    @test axb[2][3] ≈ 22.0 atol=1e-12
    ae = [[2.0,3.0],[5.0,6.0,3.0]]
    @test_throws ArgumentError pvxpv(ae, b)
    @test_throws ArgumentError pvxpv(a, ae)
end

@testset "pvstar" begin
    pv = [[126668.5912743160601,2136.792716839935195,-245251.2339876830091],
          [-0.4051854035740712739e-2,-0.6253919754866173866e-2,0.1189353719774107189e-1]]
    ra, dec, pmr, pmd, px, rv = pvstar(pv)
    @test ra ≈ 0.1686756e-1 atol=1e-12
    @test dec ≈ -1.093989828 atol=1e-12
    @test pmr ≈ -0.1783235160000472788e-4 atol=1e-16
    @test pmd ≈ 0.2336024047000619347e-5 atol=1e-16
    @test px ≈ 0.74723 atol=1e-12
    @test rv ≈ -21.60000010107306010 atol=1e-11
    pve = [[2137.792716839935195,-245251.2339876830091],
          [-0.4051854035740712739e-2,-0.6253919754866173866e-2,0.1189353719774107189e-1]]
    @test_throws ArgumentError pvstar(pve)
    @test_throws ERFAException pvstar([pv[1], [1e12, 0.0, 0.0]])
    @test_logs (:warn,) pvstar([zeros(3), zeros(3)])
end

@testset "pvtob" begin
    elong = 2.0
    phi = 0.5
    hm = 3000.0
    xp = 1e-6
    yp = -0.5e-6
    sp = 1e-8
    theta = 5.0
    pv = pvtob(elong, phi, hm, xp, yp, sp, theta)
    @test pv[1][1] ≈ 4225081.367071159207 atol=1e-5
    @test pv[1][2] ≈ 3681943.215856198144 atol=1e-5
    @test pv[1][3] ≈ 3041149.399241260785 atol=1e-5
    @test pv[2][1] ≈ -268.4915389365998787 atol=1e-9
    @test pv[2][2] ≈ 308.0977983288903123 atol=1e-9
    @test pv[2][3] ≈ 0 atol=1e-0
end

@testset "pvu" begin
    dt = 2920.0
    pv = [[126668.5912743160734,2136.792716839935565,-245251.2339876830229],
          [-0.4051854035740713039e-2,-0.6253919754866175788e-2,0.1189353719774107615e-1]]
    upv = pvu(dt, pv)
    @test upv[1][1] ≈ 126656.7598605317105 atol=1e-6
    @test upv[1][2] ≈ 2118.531271155726332 atol=1e-8
    @test upv[1][3] ≈ -245216.5048590656190 atol=1e-6
    @test upv[2][1] ≈ -0.4051854035740713039e-2 atol=1e-12
    @test upv[2][2] ≈ -0.6253919754866175788e-2 atol=1e-12
    @test upv[2][3] ≈ 0.1189353719774107615e-1 atol=1e-12
    pve = [[2136.792716839935565,-245251.2339876830229],
          [-0.4051854035740713039e-2,-0.6253919754866175788e-2,0.1189353719774107615e-1]]
    @test_throws ArgumentError pvu(dt, pve)
end

@testset "pvup" begin
    dt = 2920.0
    pv = [[126668.5912743160734,2136.792716839935565,-245251.2339876830229],
          [-0.4051854035740713039e-2,-0.6253919754866175788e-2,0.1189353719774107615e-1]]
    p = pvup(dt, pv)
    @test p[1] ≈ 126656.7598605317105 atol=1e-6
    @test p[2] ≈ 2118.531271155726332 atol=1e-8
    @test p[3] ≈ -245216.5048590656190 atol=1e-6
    pve = [[2136.792716839935565,-245251.2339876830229],
          [-0.4051854035740713039e-2,-0.6253919754866175788e-2,0.1189353719774107615e-1]]
    @test_throws ArgumentError pvup(dt, pve)
end

@testset "pxp" begin
    a = [2.0,2.0,3.0]
    b = [1.0,3.0,4.0]
    axb = ERFA._pxp(a, b)
    @test axb[1] ≈ -1.0 atol=1e-12
    @test axb[2] ≈ -5.0 atol=1e-12
    @test axb[3] ≈ 4.0 atol=1e-12
    @test_throws ArgumentError ERFA._pxp([2.0,2.0], [1.0,3.0,4.0])
    @test_throws ArgumentError ERFA._pxp([2.0,2.0,3.0], [3.0,4.0])
    @test axb ≈ cross(a, b)
end


text:
module YaDistances

using LinearAlgebra

export cosine_similarity

const AVT = AbstractArray{T, N} where {T <: Any, N}

function cosine_similarity(v₁::AVT, v₂::AVT)::Float64 # where {T <: Any}
	dot(v₁, v₂) / (norm(v₁) * norm(v₂))
end

end  ## Module

max_stars_repo_path: Data Science From Scratch/src/YaDistances.jl

max_stars_repo_name: pascal-p/julia-notebooks

max_issues_repo_name: pascal-p/julia-notebooks

max_forks_repo_name: pascal-p/julia-notebooks


text:
"Single World Intervention Programs!"
module SWIP

import Distributions
import Cassette

# export ~, ⟂, ud

include("miniomega.jl")
include("intervene.jl")
# input("graph.jl")

end # module


text:
# using Genie, Genie.Router
# using Genie.Renderer, Genie.Renderer.Html, Genie.Renderer.Json

using JSON
using SwaggerMarkdown

@swagger """
/doge:
  get:
    description: doge
    responses:
      '200':
        description: Returns a doge.
"""

@swagger """
/to:
  put:
    description: to
    responses:
      '200':
        description: Returns a doge.
"""

@swagger """
/the:
  post:
    description: the
    responses:
      '200':
        description: Returns a doge.
"""

@swagger """
/moon:
  get:
    description: moon
    responses:
      '200':
        description: Returns a moon.
"""

function run_example()
    info = Dict{String, Any}([
        "title"   => "Doge to the moon"
        "version" => "1.0.5"
    ])

    optional_fields = Dict{String, Any}([
        "host"     => "localhost"
        "basePath" => "/doge"
        # "bad"      => "bad"
    ])

    openApi = OpenAPI("2.0", info, optional_fields=optional_fields)

    spec = build(openApi)
    println(JSON.json(spec))
end

run_example()


text:
using GitHub, JSON, PkgBenchmark

commit = benchmarkpkg("Krylov")  # current state of repository
# master = benchmarkpkg("Krylov", "master")
# judgement = judge(commit, master)
# export_markdown("benchmark.md", judgement)
export_markdown("benchmark.md", commit)

gist_json = JSON.parse(
			"""
            {
            "description": "A benchmark for Krylov repository",
            "public": true,
            "files": {
                "benchmark.md": {
                "content": "$(escape_string(sprint(export_markdown, commit)))"
                }
            }
            }
            """
        )

# Need to add GITHUB_AUTH to your .bashrc
myauth = GitHub.authenticate(ENV["GITHUB_AUTH"])
posted_gist = create_gist(params = gist_json, auth = myauth)
println(posted_gist.html_url)


text:
using DelimitedFiles

function read_data(filename)
    open(filename, "r") do f
        data = split(replace(read(f, String), "| " => "", "\n" => " "), " ", keepempty = false)
    end
end

function count_nums(num_list, num)
    if num == 1
        numlen = 2
    elseif num == 7
        numlen = 3
    elseif num == 4
        numlen = 4
    elseif num == 8
        numlen = 7
    else
        numlen = 0
    end

    num_count = 0
    for n in num_list
        if length(n) == numlen
            num_count += 1
        end
    end
    return num_count
end

ones = 0
sevens = 0
fours = 0
eights = 0

j = 0
filename = "data_8.dat"
data = read_data(filename)
for ii in 1:Int(length(data)/14)
    all_nums = data[1 + j*14:10 + j*14]
    decode_nums = data[11 + j*14:14 + j*14]
    global ones += count_nums(decode_nums, 1)
    global sevens += count_nums(decode_nums, 7)
    global fours += count_nums(decode_nums, 4)
    global eights += count_nums(decode_nums, 8)

    global j += 1
end

println(ones + sevens + fours + eights)


text:
# Indexing and dimensions
export workitemIdx, workgroupIdx, workgroupDim, gridDim, gridDimWG
export threadIdx, blockIdx, blockDim

@generated function _index(::Val{fname}, ::Val{name}, ::Val{range}) where {fname, name, range}
    Context() do ctx
        T_int32 = LLVM.Int32Type(ctx)

        # create function
        llvm_f, _ = create_function(T_int32)
        mod = LLVM.parent(llvm_f)

        # generate IR
        Builder(ctx) do builder
            entry = BasicBlock(llvm_f, "entry"; ctx)
            position!(builder, entry)

            # call the indexing intrinsic
            intr_typ = LLVM.FunctionType(T_int32)
            intr = LLVM.Function(mod, "llvm.amdgcn.$fname.id.$name", intr_typ)
            idx = call!(builder, intr)

            # attach range metadata
            range_metadata = MDNode([ConstantInt(UInt32(range.start); ctx),
                                     ConstantInt(UInt32(range.stop); ctx)];
                                    ctx)
            metadata(idx)[LLVM.MD_range] = range_metadata
            ret!(builder, idx)
        end

        call_function(llvm_f, UInt32)
    end
end

@generated function _dim(::Val{base}, ::Val{off}, ::Val{range}, ::Type{T}) where {base, off, range, T}
    Context() do ctx
        T_int8 = LLVM.Int8Type(ctx)
        T_int32 = LLVM.Int32Type(ctx)
        _as = convert(Int, AS.Constant)
        T_ptr_i8 = LLVM.PointerType(T_int8, _as)
        T_ptr_i32 = LLVM.PointerType(T_int32, _as)
        T_ptr_T = LLVM.PointerType(convert(LLVMType, T; ctx), _as)

        # create function
        llvm_f, _ = create_function(T_int32)
        mod = LLVM.parent(llvm_f)

        # generate IR
        Builder(ctx) do builder
            entry = BasicBlock(llvm_f, "entry"; ctx)
            position!(builder, entry)

            # get the kernel dispatch pointer
            intr_typ = LLVM.FunctionType(T_ptr_i8)
            intr = LLVM.Function(mod, "llvm.amdgcn.dispatch.ptr", intr_typ)
            ptr = call!(builder, intr)

            # load the index
            offset = base+((off-1)*sizeof(T))
            idx_ptr_i8 = inbounds_gep!(builder, ptr, [ConstantInt(offset; ctx)])
            idx_ptr_T = bitcast!(builder, idx_ptr_i8, T_ptr_T)
            idx_T = load!(builder, idx_ptr_T)
            idx = zext!(builder, idx_T, T_int32)

            # attach range metadata
            range_metadata = MDNode([ConstantInt(T(range.start); ctx),
                                     ConstantInt(T(range.stop); ctx)];
                                    ctx)
            metadata(idx_T)[LLVM.MD_range] = range_metadata
            ret!(builder, idx)
        end

        call_function(llvm_f, UInt32)
    end
end

# TODO: look these up for the current agent/queue
# TODO: grids can be up to typemax(UInt64)
const _max_group_size = (x=1023, y=1023, z=1023)
const _max_groups     = (x=typemax(UInt32)-1, y=typemax(UInt32)-1, z=typemax(UInt32)-1)
const _max_grid_size  = (x=typemax(UInt32)-1, y=typemax(UInt32)-1, z=typemax(UInt32)-1)

for dim in (:x, :y, :z)
    # Workitem index
    fname = Symbol("workitem")
    fn = Symbol("workitemIdx_$dim")
    intr = Symbol("$dim")
    @eval @inline $fn() = Int(_index($(Val(fname)), $(Val(intr)), $(Val(0:_max_group_size[dim])))) + 1
    cufn = Symbol("threadIdx_$dim")
    @eval @inline $cufn() = $fn()

    # Workgroup index
    fname = Symbol("workgroup")
    fn = Symbol("workgroupIdx_$dim")
    intr = Symbol("$dim")
    @eval @inline $fn() = Int(_index($(Val(fname)), $(Val(intr)), $(Val(0:_max_groups[dim])))) + 1
    cufn = Symbol("blockIdx_$dim")
    @eval @inline $cufn() = $fn()
end
for (dim,off) in ((:x,1), (:y,2), (:z,3))
    # Workgroup dimension (in workitems)
    fn = Symbol("workgroupDim_$dim")
    base = _packet_offsets[findfirst(x->x==:workgroup_size_x,_packet_names)]
    @eval @inline $fn() = Int(_dim($(Val(base)), $(Val(off)), $(Val(0:_max_group_size[dim])), UInt16))
    cufn = Symbol("blockDim_$dim")
    @eval @inline $cufn() = $fn()

    # Grid dimension (in workitems)
    fn = Symbol("gridDim_$dim")
    base = _packet_offsets[findfirst(x->x==:grid_size_x,_packet_names)]
    @eval @inline $fn() = Int(_dim($(Val(base)), $(Val(off)), $(Val(0:_max_grid_size[dim])), UInt32))
    # Grid dimension (in workgroups)
    fn_wg = Symbol("gridDimWG_$dim")
    fn_wg_dim = Symbol("workgroupDim_$dim")
    @eval @inline $fn_wg() = div($fn(), $fn_wg_dim())
end

"""
    workitemIdx()::ROCDim3

Returns the work item index within the work group.
See also: [`threadIdx`](@ref)
"""
@inline workitemIdx() = (x=workitemIdx_x(), y=workitemIdx_y(), z=workitemIdx_z())

"""
    workgroupIdx()::ROCDim3

Returns the work group index.
See also: [`blockIdx`](@ref)
"""
@inline workgroupIdx() = (x=workgroupIdx_x(), y=workgroupIdx_y(), z=workgroupIdx_z())

"""
    workgroupDim()::ROCDim3

Returns the size of each workgroup in workitems.
See also: [`blockDim`](@ref)
"""
@inline workgroupDim() = (x=workgroupDim_x(), y=workgroupDim_y(), z=workgroupDim_z())

"""
    gridDim()::ROCDim3

Returns the size of the grid in workitems.
This behaviour is different from CUDA where `gridDim` gives the size of the grid in blocks.
"""
@inline gridDim() = (x=gridDim_x(), y=gridDim_y(), z=gridDim_z())

"""
    gridDimWG()::ROCDim3

Returns the size of the grid in workgroups.
This is equivalent to CUDA's `gridDim`.
"""
@inline gridDimWG() = (x=gridDimWG_x(), y=gridDimWG_y(), z=gridDimWG_z())

# For compat with CUDAnative et. al

"""
    threadIdx()::ROCDim3

Returns the thread index within the block.
See also: [`workitemIdx`](@ref)
"""
@inline threadIdx() = (x=threadIdx_x(), y=threadIdx_y(), z=threadIdx_z())

"""
    blockIdx()::ROCDim3

Returns the block index within the grid.
See also: [`workgroupIdx`](@ref)
"""
@inline blockIdx() = (x=blockIdx_x(), y=blockIdx_y(), z=blockIdx_z())

"""
    blockDim()::ROCDim3

Returns the dimensions of the block.
See also: [`workgroupDim`](@ref)
"""
@inline blockDim() = (x=blockDim_x(), y=blockDim_y(), z=blockDim_z())