text:
baremodule Semirings
    using ..Types
end

SemiringUnion = Union{AbstractSemiring, libgb.GrB_Semiring}

function semiringnames(name)
    simple = name[5:end]
    containername = Symbol(simple, "_RIG_T")
    exportedname = Symbol(simple)
    return containername, exportedname
end

function createsemirings()
    builtins = ["GxB_PLUS_FIRST", "GxB_TIMES_FIRST", "GxB_ANY_FIRST", "GxB_PLUS_SECOND", "GxB_TIMES_SECOND", "GxB_ANY_SECOND", "GxB_MIN_PAIR", "GxB_MAX_PAIR", "GxB_PLUS_PAIR", "GxB_TIMES_PAIR", "GxB_ANY_PAIR", "GxB_MIN_MIN", "GxB_TIMES_MIN", "GxB_ANY_MIN", "GxB_MAX_MAX", "GxB_PLUS_MAX", "GxB_TIMES_MAX", "GxB_ANY_MAX", "GxB_PLUS_PLUS", "GxB_TIMES_PLUS", "GxB_ANY_PLUS", "GxB_MIN_MINUS", "GxB_MAX_MINUS", "GxB_PLUS_MINUS", "GxB_TIMES_MINUS", "GxB_ANY_MINUS", "GxB_TIMES_TIMES", "GxB_ANY_TIMES", "GxB_MIN_DIV", "GxB_MAX_DIV", "GxB_PLUS_DIV", "GxB_TIMES_DIV", "GxB_ANY_DIV", "GxB_MIN_RDIV", "GxB_MAX_RDIV", "GxB_PLUS_RDIV", "GxB_TIMES_RDIV", "GxB_ANY_RDIV", "GxB_MIN_RMINUS", "GxB_MAX_RMINUS", "GxB_PLUS_RMINUS", "GxB_TIMES_RMINUS", "GxB_ANY_RMINUS", "GxB_MIN_ISEQ", "GxB_MAX_ISEQ", "GxB_PLUS_ISEQ", "GxB_TIMES_ISEQ", "GxB_ANY_ISEQ", "GxB_MIN_ISNE", "GxB_MAX_ISNE", "GxB_PLUS_ISNE", "GxB_TIMES_ISNE", "GxB_ANY_ISNE", "GxB_MIN_ISGT", "GxB_MAX_ISGT", "GxB_PLUS_ISGT", "GxB_TIMES_ISGT", "GxB_ANY_ISGT", "GxB_MIN_ISLT", "GxB_MAX_ISLT", "GxB_PLUS_ISLT", "GxB_TIMES_ISLT", "GxB_ANY_ISLT", "GxB_MIN_ISGE", "GxB_MAX_ISGE", "GxB_PLUS_ISGE", "GxB_TIMES_ISGE", "GxB_ANY_ISGE", "GxB_MIN_ISLE", "GxB_MAX_ISLE", "GxB_PLUS_ISLE", "GxB_TIMES_ISLE", "GxB_ANY_ISLE", "GxB_MIN_LOR", "GxB_MAX_LOR", "GxB_PLUS_LOR", "GxB_TIMES_LOR", "GxB_ANY_LOR", "GxB_MIN_LAND", "GxB_MAX_LAND", "GxB_PLUS_LAND", "GxB_TIMES_LAND", "GxB_ANY_LAND", "GxB_MIN_LXOR", "GxB_MAX_LXOR", "GxB_PLUS_LXOR", "GxB_TIMES_LXOR", "GxB_ANY_LXOR", "GxB_LOR_NE", "GxB_LOR_EQ", "GxB_LAND_EQ", "GxB_LXOR_EQ", "GxB_EQ_EQ", "GxB_ANY_EQ", "GxB_LAND_NE", "GxB_LXOR_NE", "GxB_EQ_NE", "GxB_ANY_NE", "GxB_LOR_GT", "GxB_LAND_GT", "GxB_LXOR_GT", "GxB_EQ_GT", "GxB_ANY_GT", "GxB_LOR_LT", "GxB_LAND_LT", "GxB_LXOR_LT", "GxB_EQ_LT", "GxB_ANY_LT", "GxB_LOR_GE", "GxB_LAND_GE", "GxB_LXOR_GE", "GxB_EQ_GE", "GxB_ANY_GE", "GxB_LOR_LE", "GxB_LAND_LE", "GxB_LXOR_LE", "GxB_EQ_LE", "GxB_ANY_LE", "GxB_LOR_FIRST", "GxB_LAND_FIRST", "GxB_LXOR_FIRST", "GxB_EQ_FIRST", "GxB_LOR_SECOND", "GxB_LAND_SECOND", "GxB_LXOR_SECOND", "GxB_EQ_SECOND", "GxB_LOR_PAIR", "GxB_LAND_PAIR", "GxB_LXOR_PAIR", "GxB_EQ_PAIR", "GxB_LOR_LOR", "GxB_LXOR_LOR", "GxB_EQ_LOR", "GxB_LAND_LAND", "GxB_EQ_LAND", "GxB_LOR_LXOR", "GxB_LAND_LXOR", "GxB_LXOR_LXOR", "GxB_EQ_LXOR", "GxB_BOR_BOR", "GxB_BOR_BAND", "GxB_BOR_BXOR", "GxB_BOR_BXNOR", "GxB_BAND_BOR", "GxB_BAND_BAND", "GxB_BAND_BXOR", "GxB_BAND_BXNOR", "GxB_BXOR_BOR", "GxB_BXOR_BAND", "GxB_BXOR_BXOR", "GxB_BXOR_BXNOR", "GxB_BXNOR_BOR", "GxB_BXNOR_BAND", "GxB_BXNOR_BXOR", "GxB_BXNOR_BXNOR", "GxB_MIN_FIRSTI", "GxB_MAX_FIRSTI", "GxB_ANY_FIRSTI", "GxB_PLUS_FIRSTI", "GxB_TIMES_FIRSTI", "GxB_MIN_FIRSTI1", "GxB_MAX_FIRSTI1", "GxB_ANY_FIRSTI1", "GxB_PLUS_FIRSTI1", "GxB_TIMES_FIRSTI1", "GxB_MIN_FIRSTJ", "GxB_MAX_FIRSTJ", "GxB_ANY_FIRSTJ", "GxB_PLUS_FIRSTJ", "GxB_TIMES_FIRSTJ", "GxB_MIN_FIRSTJ1", "GxB_MAX_FIRSTJ1", "GxB_ANY_FIRSTJ1", "GxB_PLUS_FIRSTJ1", "GxB_TIMES_FIRSTJ1", "GxB_MIN_SECONDI", "GxB_MAX_SECONDI", "GxB_ANY_SECONDI", "GxB_PLUS_SECONDI", "GxB_TIMES_SECONDI", "GxB_MIN_SECONDI1", "GxB_MAX_SECONDI1", "GxB_ANY_SECONDI1", "GxB_PLUS_SECONDI1", "GxB_TIMES_SECONDI1", "GxB_MIN_SECONDJ", "GxB_MAX_SECONDJ", "GxB_ANY_SECONDJ", "GxB_PLUS_SECONDJ", "GxB_TIMES_SECONDJ", "GxB_MIN_SECONDJ1", "GxB_MAX_SECONDJ1", "GxB_ANY_SECONDJ1", "GxB_PLUS_SECONDJ1", "GxB_TIMES_SECONDJ1", "GrB_PLUS_TIMES_SEMIRING", "GrB_PLUS_MIN_SEMIRING", "GrB_MIN_PLUS_SEMIRING", "GrB_MIN_TIMES_SEMIRING", "GrB_MIN_FIRST_SEMIRING", "GrB_MIN_SECOND_SEMIRING", "GrB_MIN_MAX_SEMIRING", "GrB_MAX_PLUS_SEMIRING", "GrB_MAX_TIMES_SEMIRING", "GrB_MAX_FIRST_SEMIRING", "GrB_MAX_SECOND_SEMIRING", "GrB_MAX_MIN_SEMIRING", "GrB_LOR_LAND_SEMIRING", "GrB_LAND_LOR_SEMIRING", "GrB_LXOR_LAND_SEMIRING", "GrB_LXNOR_LOR_SEMIRING"]
    for name ∈ builtins
        containername, exportedname = semiringnames(name)
        structquote = quote
            struct $containername <: AbstractSemiring
                pointers::Dict{DataType, libgb.GrB_Semiring}
                name::String
                $containername() = new(Dict{DataType, libgb.GrB_Semiring}(), $name)
            end
        end
        @eval(Types, $structquote)
        constquote = quote
            const $exportedname = Types.$containername()
            export $exportedname
        end
        @eval(Semirings,$constquote)
    end
end

function load(rig::AbstractSemiring)
    booleans = ["GxB_LOR_FIRST", "GxB_LAND_FIRST", "GxB_LXOR_FIRST", "GxB_EQ_FIRST", "GxB_ANY_FIRST", "GxB_LOR_SECOND", "GxB_LAND_SECOND", "GxB_LXOR_SECOND", "GxB_EQ_SECOND", "GxB_ANY_SECOND", "GxB_LOR_PAIR", "GxB_LAND_PAIR", "GxB_LXOR_PAIR", "GxB_EQ_PAIR", "GxB_ANY_PAIR", "GxB_LOR_LOR", "GxB_LXOR_LOR", "GxB_EQ_LOR", "GxB_ANY_LOR", "GxB_LAND_LAND", "GxB_EQ_LAND", "GxB_ANY_LAND", "GxB_LOR_LXOR", "GxB_LAND_LXOR", "GxB_LXOR_LXOR", "GxB_EQ_LXOR", "GxB_ANY_LXOR", "GxB_LOR_EQ", "GxB_LAND_EQ", "GxB_LXOR_EQ", "GxB_EQ_EQ", "GxB_ANY_EQ", "GxB_LOR_GT", "GxB_LAND_GT", "GxB_LXOR_GT", "GxB_EQ_GT", "GxB_ANY_GT", "GxB_LOR_LT", "GxB_LAND_LT", "GxB_LXOR_LT", "GxB_EQ_LT", "GxB_ANY_LT", "GxB_LOR_GE", "GxB_LAND_GE", "GxB_LXOR_GE", "GxB_EQ_GE", "GxB_ANY_GE", "GxB_LOR_LE", "GxB_LAND_LE", "GxB_LXOR_LE", "GxB_EQ_LE", "GxB_ANY_LE", "GrB_LOR_LAND_SEMIRING", "GrB_LAND_LOR_SEMIRING", "GrB_LXOR_LAND_SEMIRING", "GrB_LXNOR_LOR_SEMIRING"]

    integers = ["GxB_PLUS_FIRST", "GxB_TIMES_FIRST", "GxB_ANY_FIRST", "GxB_PLUS_SECOND", "GxB_TIMES_SECOND", "GxB_ANY_SECOND", "GxB_MIN_PAIR", "GxB_MAX_PAIR", "GxB_PLUS_PAIR", "GxB_TIMES_PAIR", "GxB_ANY_PAIR", "GxB_MIN_MIN", "GxB_TIMES_MIN", "GxB_ANY_MIN", "GxB_MAX_MAX", "GxB_PLUS_MAX", "GxB_TIMES_MAX", "GxB_ANY_MAX", "GxB_PLUS_PLUS", "GxB_TIMES_PLUS", "GxB_ANY_PLUS", "GxB_MIN_MINUS", "GxB_MAX_MINUS", "GxB_PLUS_MINUS", "GxB_TIMES_MINUS", "GxB_ANY_MINUS", "GxB_TIMES_TIMES", "GxB_ANY_TIMES", "GxB_MIN_DIV", "GxB_MAX_DIV", "GxB_PLUS_DIV", "GxB_TIMES_DIV", "GxB_ANY_DIV", "GxB_MIN_RDIV", "GxB_MAX_RDIV", "GxB_PLUS_RDIV", "GxB_TIMES_RDIV", "GxB_ANY_RDIV", "GxB_MIN_RMINUS", "GxB_MAX_RMINUS", "GxB_PLUS_RMINUS", "GxB_TIMES_RMINUS", "GxB_ANY_RMINUS", "GxB_MIN_ISEQ", "GxB_MAX_ISEQ", "GxB_PLUS_ISEQ", "GxB_TIMES_ISEQ", "GxB_ANY_ISEQ", "GxB_MIN_ISNE", "GxB_MAX_ISNE", "GxB_PLUS_ISNE", "GxB_TIMES_ISNE", "GxB_ANY_ISNE", "GxB_MIN_ISGT", "GxB_MAX_ISGT", "GxB_PLUS_ISGT", "GxB_TIMES_ISGT", "GxB_ANY_ISGT", "GxB_MIN_ISLT", "GxB_MAX_ISLT", "GxB_PLUS_ISLT", "GxB_TIMES_ISLT", "GxB_ANY_ISLT", "GxB_MIN_ISGE", "GxB_MAX_ISGE", "GxB_PLUS_ISGE", "GxB_TIMES_ISGE", "GxB_ANY_ISGE", "GxB_MIN_ISLE", "GxB_MAX_ISLE", "GxB_PLUS_ISLE", "GxB_TIMES_ISLE", "GxB_ANY_ISLE", "GxB_MIN_LOR", "GxB_MAX_LOR", "GxB_PLUS_LOR", "GxB_TIMES_LOR", "GxB_ANY_LOR", "GxB_MIN_LAND", "GxB_MAX_LAND", "GxB_PLUS_LAND", "GxB_TIMES_LAND", "GxB_ANY_LAND", "GxB_MIN_LXOR", "GxB_MAX_LXOR", "GxB_PLUS_LXOR", "GxB_TIMES_LXOR", "GxB_ANY_LXOR", "GxB_LOR_NE", "GxB_LOR_EQ", "GxB_LAND_EQ", "GxB_LXOR_EQ", "GxB_EQ_EQ", "GxB_ANY_EQ", "GxB_LAND_NE", "GxB_LXOR_NE", "GxB_EQ_NE", "GxB_ANY_NE", "GxB_LOR_GT", "GxB_LAND_GT", "GxB_LXOR_GT", "GxB_EQ_GT", "GxB_ANY_GT", "GxB_LOR_LT", "GxB_LAND_LT", "GxB_LXOR_LT", "GxB_EQ_LT", "GxB_ANY_LT", "GxB_LOR_GE", "GxB_LAND_GE", "GxB_LXOR_GE", "GxB_EQ_GE", "GxB_ANY_GE", "GxB_LOR_LE", "GxB_LAND_LE", "GxB_LXOR_LE", "GxB_EQ_LE", "GxB_ANY_LE", "GrB_PLUS_TIMES_SEMIRING", "GrB_PLUS_MIN_SEMIRING", "GrB_MIN_PLUS_SEMIRING", "GrB_MIN_TIMES_SEMIRING", "GrB_MIN_FIRST_SEMIRING", "GrB_MIN_SECOND_SEMIRING", "GrB_MIN_MAX_SEMIRING", "GrB_MAX_PLUS_SEMIRING", "GrB_MAX_TIMES_SEMIRING", "GrB_MAX_FIRST_SEMIRING", "GrB_MAX_SECOND_SEMIRING", "GrB_MAX_MIN_SEMIRING"]

    unsignedintegers = ["GxB_PLUS_FIRST", "GxB_TIMES_FIRST", "GxB_ANY_FIRST", "GxB_PLUS_SECOND", "GxB_TIMES_SECOND", "GxB_ANY_SECOND", "GxB_MIN_PAIR", "GxB_MAX_PAIR", "GxB_PLUS_PAIR", "GxB_TIMES_PAIR", "GxB_ANY_PAIR", "GxB_MIN_MIN", "GxB_TIMES_MIN", "GxB_ANY_MIN", "GxB_MAX_MAX", "GxB_PLUS_MAX", "GxB_TIMES_MAX", "GxB_ANY_MAX", "GxB_PLUS_PLUS", "GxB_TIMES_PLUS", "GxB_ANY_PLUS", "GxB_MIN_MINUS", "GxB_MAX_MINUS", "GxB_PLUS_MINUS", "GxB_TIMES_MINUS", "GxB_ANY_MINUS", "GxB_TIMES_TIMES", "GxB_ANY_TIMES", "GxB_MIN_DIV", "GxB_MAX_DIV", "GxB_PLUS_DIV", "GxB_TIMES_DIV", "GxB_ANY_DIV", "GxB_MIN_RDIV", "GxB_MAX_RDIV", "GxB_PLUS_RDIV", "GxB_TIMES_RDIV", "GxB_ANY_RDIV", "GxB_MIN_RMINUS", "GxB_MAX_RMINUS", "GxB_PLUS_RMINUS", "GxB_TIMES_RMINUS", "GxB_ANY_RMINUS", "GxB_MIN_ISEQ", "GxB_MAX_ISEQ", "GxB_PLUS_ISEQ", "GxB_TIMES_ISEQ", "GxB_ANY_ISEQ", "GxB_MIN_ISNE", "GxB_MAX_ISNE", "GxB_PLUS_ISNE", "GxB_TIMES_ISNE", "GxB_ANY_ISNE", "GxB_MIN_ISGT", "GxB_MAX_ISGT", "GxB_PLUS_ISGT", "GxB_TIMES_ISGT", "GxB_ANY_ISGT", "GxB_MIN_ISLT", "GxB_MAX_ISLT", "GxB_PLUS_ISLT", "GxB_TIMES_ISLT", "GxB_ANY_ISLT", "GxB_MIN_ISGE", "GxB_MAX_ISGE", "GxB_PLUS_ISGE", "GxB_TIMES_ISGE", "GxB_ANY_ISGE", "GxB_MIN_ISLE", "GxB_MAX_ISLE", "GxB_PLUS_ISLE", "GxB_TIMES_ISLE", "GxB_ANY_ISLE", "GxB_MIN_LOR", "GxB_MAX_LOR", "GxB_PLUS_LOR", "GxB_TIMES_LOR", "GxB_ANY_LOR", "GxB_MIN_LAND", "GxB_MAX_LAND", "GxB_PLUS_LAND", "GxB_TIMES_LAND", "GxB_ANY_LAND", "GxB_MIN_LXOR", "GxB_MAX_LXOR", "GxB_PLUS_LXOR", "GxB_TIMES_LXOR", "GxB_ANY_LXOR", "GxB_LOR_NE", "GxB_LOR_EQ", "GxB_LAND_EQ", "GxB_LXOR_EQ", "GxB_EQ_EQ", "GxB_ANY_EQ", "GxB_LAND_NE", "GxB_LXOR_NE", "GxB_EQ_NE", "GxB_ANY_NE", "GxB_LOR_GT", "GxB_LAND_GT", "GxB_LXOR_GT", "GxB_EQ_GT", "GxB_ANY_GT", "GxB_LOR_LT", "GxB_LAND_LT", "GxB_LXOR_LT", "GxB_EQ_LT", "GxB_ANY_LT", "GxB_LOR_GE", "GxB_LAND_GE", "GxB_LXOR_GE", "GxB_EQ_GE", "GxB_ANY_GE", "GxB_LOR_LE", "GxB_LAND_LE", "GxB_LXOR_LE", "GxB_EQ_LE", "GxB_ANY_LE", "GxB_BOR_BOR", "GxB_BOR_BAND", "GxB_BOR_BXOR", "GxB_BOR_BXNOR", "GxB_BAND_BOR", "GxB_BAND_BAND", "GxB_BAND_BXOR", "GxB_BAND_BXNOR", "GxB_BXOR_BOR", "GxB_BXOR_BAND", "GxB_BXOR_BXOR", "GxB_BXOR_BXNOR", "GxB_BXNOR_BOR", "GxB_BXNOR_BAND", "GxB_BXNOR_BXOR", "GxB_BXNOR_BXNOR", "GrB_PLUS_TIMES_SEMIRING", "GrB_PLUS_MIN_SEMIRING", "GrB_MIN_PLUS_SEMIRING", "GrB_MIN_TIMES_SEMIRING", "GrB_MIN_FIRST_SEMIRING", "GrB_MIN_SECOND_SEMIRING", "GrB_MIN_MAX_SEMIRING", "GrB_MAX_PLUS_SEMIRING", "GrB_MAX_TIMES_SEMIRING", "GrB_MAX_FIRST_SEMIRING", "GrB_MAX_SECOND_SEMIRING", "GrB_MAX_MIN_SEMIRING"]

    floats = ["GxB_PLUS_FIRST", "GxB_TIMES_FIRST", "GxB_ANY_FIRST", "GxB_PLUS_SECOND", "GxB_TIMES_SECOND", "GxB_ANY_SECOND", "GxB_MIN_PAIR", "GxB_MAX_PAIR", "GxB_PLUS_PAIR", "GxB_TIMES_PAIR", "GxB_ANY_PAIR", "GxB_MIN_MIN", "GxB_TIMES_MIN", "GxB_ANY_MIN", "GxB_MAX_MAX", "GxB_PLUS_MAX", "GxB_TIMES_MAX", "GxB_ANY_MAX", "GxB_PLUS_PLUS", "GxB_TIMES_PLUS", "GxB_ANY_PLUS", "GxB_MIN_MINUS", "GxB_MAX_MINUS", "GxB_PLUS_MINUS", "GxB_TIMES_MINUS", "GxB_ANY_MINUS", "GxB_TIMES_TIMES", "GxB_ANY_TIMES", "GxB_MIN_DIV", "GxB_MAX_DIV", "GxB_PLUS_DIV", "GxB_TIMES_DIV", "GxB_ANY_DIV", "GxB_MIN_RDIV", "GxB_MAX_RDIV", "GxB_PLUS_RDIV", "GxB_TIMES_RDIV", "GxB_ANY_RDIV", "GxB_MIN_RMINUS", "GxB_MAX_RMINUS", "GxB_PLUS_RMINUS", "GxB_TIMES_RMINUS", "GxB_ANY_RMINUS", "GxB_MIN_ISEQ", "GxB_MAX_ISEQ", "GxB_PLUS_ISEQ", "GxB_TIMES_ISEQ", "GxB_ANY_ISEQ", "GxB_MIN_ISNE", "GxB_MAX_ISNE", "GxB_PLUS_ISNE", "GxB_TIMES_ISNE", "GxB_ANY_ISNE", "GxB_MIN_ISGT", "GxB_MAX_ISGT", "GxB_PLUS_ISGT", "GxB_TIMES_ISGT", "GxB_ANY_ISGT", "GxB_MIN_ISLT", "GxB_MAX_ISLT", "GxB_PLUS_ISLT", "GxB_TIMES_ISLT", "GxB_ANY_ISLT", "GxB_MIN_ISGE", "GxB_MAX_ISGE", "GxB_PLUS_ISGE", "GxB_TIMES_ISGE", "GxB_ANY_ISGE", "GxB_MIN_ISLE", "GxB_MAX_ISLE", "GxB_PLUS_ISLE", "GxB_TIMES_ISLE", "GxB_ANY_ISLE", "GxB_MIN_LOR", "GxB_MAX_LOR", "GxB_PLUS_LOR", "GxB_TIMES_LOR", "GxB_ANY_LOR", "GxB_MIN_LAND", "GxB_MAX_LAND", "GxB_PLUS_LAND", "GxB_TIMES_LAND", "GxB_ANY_LAND", "GxB_MIN_LXOR", "GxB_MAX_LXOR", "GxB_PLUS_LXOR", "GxB_TIMES_LXOR", "GxB_ANY_LXOR", "GxB_LOR_NE", "GxB_LOR_EQ", "GxB_LAND_EQ", "GxB_LXOR_EQ", "GxB_EQ_EQ", "GxB_ANY_EQ", "GxB_LAND_NE", "GxB_LXOR_NE", "GxB_EQ_NE", "GxB_ANY_NE", "GxB_LOR_GT", "GxB_LAND_GT", "GxB_LXOR_GT", "GxB_EQ_GT", "GxB_ANY_GT", "GxB_LOR_LT", "GxB_LAND_LT", "GxB_LXOR_LT", "GxB_EQ_LT", "GxB_ANY_LT", "GxB_LOR_GE", "GxB_LAND_GE", "GxB_LXOR_GE", "GxB_EQ_GE", "GxB_ANY_GE", "GxB_LOR_LE", "GxB_LAND_LE", "GxB_LXOR_LE", "GxB_EQ_LE", "GxB_ANY_LE", "GrB_PLUS_TIMES_SEMIRING", "GrB_PLUS_MIN_SEMIRING", "GrB_MIN_PLUS_SEMIRING", "GrB_MIN_TIMES_SEMIRING", "GrB_MIN_FIRST_SEMIRING", "GrB_MIN_SECOND_SEMIRING", "GrB_MIN_MAX_SEMIRING", "GrB_MAX_PLUS_SEMIRING", "GrB_MAX_TIMES_SEMIRING", "GrB_MAX_FIRST_SEMIRING", "GrB_MAX_SECOND_SEMIRING", "GrB_MAX_MIN_SEMIRING"]

    positionals = ["GxB_MIN_FIRSTI", "GxB_MAX_FIRSTI", "GxB_ANY_FIRSTI", "GxB_PLUS_FIRSTI", "GxB_TIMES_FIRSTI", "GxB_MIN_FIRSTI1", "GxB_MAX_FIRSTI1", "GxB_ANY_FIRSTI1", "GxB_PLUS_FIRSTI1", "GxB_TIMES_FIRSTI1", "GxB_MIN_FIRSTJ", "GxB_MAX_FIRSTJ", "GxB_ANY_FIRSTJ", "GxB_PLUS_FIRSTJ", "GxB_TIMES_FIRSTJ", "GxB_MIN_FIRSTJ1", "GxB_MAX_FIRSTJ1", "GxB_ANY_FIRSTJ1", "GxB_PLUS_FIRSTJ1", "GxB_TIMES_FIRSTJ1", "GxB_MIN_SECONDI", "GxB_MAX_SECONDI", "GxB_ANY_SECONDI", "GxB_PLUS_SECONDI", "GxB_TIMES_SECONDI", "GxB_MIN_SECONDI1", "GxB_MAX_SECONDI1", "GxB_ANY_SECONDI1", "GxB_PLUS_SECONDI1", "GxB_TIMES_SECONDI1", "GxB_MIN_SECONDJ", "GxB_MAX_SECONDJ", "GxB_ANY_SECONDJ", "GxB_PLUS_SECONDJ", "GxB_TIMES_SECONDJ", "GxB_MIN_SECONDJ1", "GxB_MAX_SECONDJ1", "GxB_ANY_SECONDJ1", "GxB_PLUS_SECONDJ1", "GxB_TIMES_SECONDJ1"]
    name = rig.name
    if name ∈ booleans
        rig.pointers[Bool] = load_global(name * "_BOOL")
    end

    if name ∈ integers
        rig.pointers[Int8] =load_global(name * "_INT8")
        rig.pointers[Int16] = load_global(name * "_INT16")
        rig.pointers[Int32] = load_global(name * "_INT32")
        rig.pointers[Int64] = load_global(name * "_INT64")
    end

    if name ∈ unsignedintegers
        rig.pointers[UInt8] =load_global(name * "_UINT8")
        rig.pointers[UInt16] = load_global(name * "_UINT16")
        rig.pointers[UInt32] = load_global(name * "_UINT32")
        rig.pointers[UInt64] = load_global(name * "_UINT64")
    end

    if name ∈ floats
        rig.pointers[Float32] = load_global(name * "_FP32")
        rig.pointers[Float64] = load_global(name * "_FP64")
    end
end

multiplyop(rig::libgb.GrB_Semiring) = libgb.GxB_Semiring_multiply(rig)
addop(rig::libgb.GrB_Semiring) = libgb.GxB_Semiring_add(rig)

xtype(rig::SemiringUnion) = xtype(addop(rig))
ytype(rig::SemiringUnion) = ytype(addop(rig))
ztype(rig::SemiringUnion) = ztype(addop(rig))


text:
# PREAMBLE

# PKG_SETUP

# ## Setup

using IterativeLQR
using LinearAlgebra
using Plots

# ## horizon
T = 101

# ## acrobot
num_state = 4
num_action = 1
num_parameter = 0

function acrobot(x, u, w)
    mass1 = 1.0
    inertia1 = 0.33
    length1 = 1.0
    lengthcom1 = 0.5

    mass2 = 1.0
    inertia2 = 0.33
    length2 = 1.0
    lengthcom2 = 0.5

    gravity = 9.81
    friction1 = 0.1
    friction2 = 0.1

    function M(x, w)
        a = (inertia1 + inertia2 + mass2 * length1 * length1
            + 2.0 * mass2 * length1 * lengthcom2 * cos(x[2]))

        b = inertia2 + mass2 * length1 * lengthcom2 * cos(x[2])

        c = inertia2

       return [a b; b c]
    end

    function Minv(x, w)
        m = M(x, w)
        a = m[1, 1]
        b = m[1, 2]
        c = m[2, 1]
        d = m[2, 2]
        1.0 / (a * d - b * c) * [d -b;-c a]
    end

    function τ(x, w)
        a = (-1.0 * mass1 * gravity * lengthcom1 * sin(x[1])
            - mass2 * gravity * (length1 * sin(x[1])
            + lengthcom2 * sin(x[1] + x[2])))

        b = -1.0 * mass2 * gravity * lengthcom2 * sin(x[1] + x[2])

        return [a; b]
    end

    function C(x, w)
        a = -2.0 * mass2 * length1 * lengthcom2 * sin(x[2]) * x[4]
        b = -1.0 * mass2 * length1 * lengthcom2 * sin(x[2]) * x[4]
        c = mass2 * length1 * lengthcom2 * sin(x[2]) * x[3]
        d = 0.0

        return [a b; c d]
    end

    function B(x, w)
        [0.0; 1.0]
    end

    q = view(x, 1:2)
    v = view(x, 3:4)

    qdd = Minv(q, w) * (-1.0 * C(x, w) * v
            + τ(q, w) + B(q, w) * u[1] - [friction1; friction2] .* v)

    return [x[3]; x[4]; qdd[1]; qdd[2]]
end

function midpoint_explicit(x, u, w)
    h = 0.1 # timestep
    x + h * acrobot(x + 0.5 * h * acrobot(x, u, w), u, w)
end

# ## model
dynamics = Dynamics(midpoint_explicit, num_state, num_action, num_parameter)
model = [dynamics for t = 1:T-1]

# ## initialization
x1 = [0.0; 0.0; 0.0; 0.0]
xT = [0.0; π; 0.0; 0.0]
ū = [1.0 * randn(num_action) for t = 1:T-1]
w = [zeros(num_parameter) for t = 1:T]
x̄ = rollout(model, x1, ū, w)

# ## objective
ot = (x, u, w) -> 0.1 * dot(x[3:4], x[3:4]) + 0.1 * dot(u, u)
oT = (x, u, w) -> 0.1 * dot(x[3:4], x[3:4])
ct = Cost(ot, num_state, num_action, num_parameter)
cT = Cost(oT, num_state, 0, num_parameter)
objective = [[ct for t = 1:T-1]..., cT]

# ## constraints
goal(x, u, w) = x - xT

cont = Constraint()
conT = Constraint(goal, num_state, 0)
constraints = [[cont for t = 1:T-1]..., conT]

# ## problem
prob = Solver(model, objective, constraints)
initialize_controls!(prob, ū)
initialize_states!(prob, x̄)

# ## solve
solve!(prob)

# ## solution
x_sol, u_sol = get_trajectory(prob)

# ## visuals
plot(hcat(x_sol...)')
plot(hcat(u_sol...)', linetype=:steppost)

# ## benchmark allocations + timing
info = @benchmark solve!($prob, x̄, ū) setup=(x̄=deepcopy(x̄), ū=deepcopy(ū))

ext: jl


text:
export precompile_activator, precompile_deactivator, precompile_pather, BotConfig, @snoopi_bot, @snoopi_bench

const UStrings = Union{AbstractString,Regex,AbstractChar}
################################################################
"""
    BotConfig

Config object that holds the options and configuration for the SnoopCompile bot. This object is fed to the `@snoopi_bot`.

# Arguments:
- `packageName::String`
- `subst::Vector{Pair{UStrings, UStrings}}` : to replace a packages precompile setences with another's package like `["ImageTest" => "Images"]`
- `blacklist::Vector{UStrings}` : to remove some precompile sentences

`const UStrings ==  Union{AbstractString,Regex,AbstractChar}` # every string like type that `replace()` has a method for.
"""
struct BotConfig
    packageName::String
    subst::Vector{Pair{T1, T2}} where {T1<:UStrings, T2 <: UStrings}
    blacklist::Vector{T3} where {T3<:UStrings}
end

function BotConfig(packageName::String; subst::Vector{Pair{T1, T2}} where {T1<:UStrings, T2 <: UStrings} = Vector{Pair{String, String}}(), blacklist::Vector{T3} where {T3<:UStrings}= String[])
    return BotConfig(packageName, subst, blacklist)
end

include("bot/botutils.jl")
include("bot/precompileInclude.jl")
include("bot/snoopi_bot.jl")
include("bot/snoopi_bench.jl")

# deprecation and backward compatiblity
macro snoopiBot(args...)
     f, l = __source__.file, __source__.line
     Base.depwarn("`@snoopiBot` at $f:$l is deprecated, rename the macro to `@snoopi_bot`.", Symbol("@snoopiBot"))
     return esc(:(@snoopi_bot($(args...))))
end
macro snoopiBench(args...)
    f, l = __source__.file, __source__.line
    Base.depwarn("`@snoopiBench` at $f:$l is deprecated, rename the macro to `@snoopi_bench`.", Symbol("@snoopiBench"))
    return esc(:(@snoopi_bench($(args...))))
end

@eval @deprecate $(Symbol("@snoopiBot")) $(Symbol("@snoopi_bot"))
@eval @deprecate $(Symbol("@snoopiBench")) $(Symbol("@snoopi_bench"))


text:
using Pkg.Artifacts
using polymake_jll

# adapted from https://github.com/JuliaLang/julia/pull/38797
function full_artifact_dir(m::Module)
   artifacts_toml = joinpath(dirname(dirname(Base.pathof(m))), "StdlibArtifacts.toml")

   # If this file exists, it's a stdlib JLL and we must download the artifact ourselves
   if isfile(artifacts_toml)
       # we need to remove the _jll for the artifact name
       meta = artifact_meta(string(m)[1:end-4], artifacts_toml)
       hash = Base.SHA1(meta["git-tree-sha1"])
       if !artifact_exists(hash)
           dl_info = first(meta["download"])
           download_artifact(hash, dl_info["url"], dl_info["sha256"])
       end
       return artifact_path(hash)
   else
      # Otherwise, we can just use the artifact directory given to us by the module
      return m.artifact_dir
   end
end

function prepare_deps_tree()
   mutable_artifacts_toml = joinpath(dirname(pathof(polymake_jll)), "..", "MutableArtifacts.toml")
   polymake_tree = "polymake_tree"
   polymake_tree_hash = artifact_hash(polymake_tree, mutable_artifacts_toml)

   # create a directory tree for polymake with links to dependencies
   # looking similiar to the tree in the build environment
   # for compiling wrappers at run-time
   polymake_tree_hash = create_artifact() do art_dir
      mkpath(joinpath(art_dir,"deps"))
      deps = [ polymake_jll.FLINT_jll,
               polymake_jll.GMP_jll,
               polymake_jll.MPFR_jll,
               polymake_jll.PPL_jll,
               polymake_jll.Perl_jll,
               polymake_jll.bliss_jll,
               polymake_jll.boost_jll,
               polymake_jll.cddlib_jll,
               polymake_jll.lrslib_jll,
               polymake_jll.normaliz_jll ]

      # dependencies
      for dep in deps
         symlink(full_artifact_dir(dep), joinpath(art_dir,"deps","$dep"))
      end

      # polymake prefix directories
      for dir in filter(d -> d!="bin", readdir(polymake_jll.artifact_dir))
         symlink(joinpath(polymake_jll.artifact_dir,dir), joinpath(art_dir,dir))
      end

      # polymake perl scripts
      bindir(name) = joinpath(art_dir,"bin", name)
      mkpath(bindir(""))
      for file in ("polymake", "polymake-config")
         symlink(joinpath(polymake_jll.artifact_dir, "bin", file), bindir(file))
      end

      # wrappers for 4ti2 utilities which set the correct LIBPATH
      bindir4ti2 = joinpath(polymake_jll.lib4ti2_jll.artifact_dir,"bin")
      bin4ti2 = joinpath(bindir4ti2,"\$(basename \$0)")

      wrapper = bindir("4ti2_wrapper")
      write(wrapper, """
         #!/bin/sh
         # Since we cannot run these binaries through the usual julia commands we need
         # this wrapper that sets up the correct library paths.
         export $(polymake_jll.JLLWrappers.LIBPATH_env)="$(polymake_jll.lib4ti2_jll.LIBPATH[])"
         exec $(bin4ti2) "\$@"
         """)

      scriptwrapper = bindir("4ti2_script_wrapper")
      chmod(wrapper, 0o755)
      write(scriptwrapper, """
         #!/bin/sh
         # We cannot run the 4ti2 scripts directly as macOS will remove
         # DYLD_FALLBACK_LIBRARY_PATH for any subshells.
         # Instead we source the 4ti2 scripts, this means \$0 will still point to the
         # symlink in our wrapper directory. Since these scripts use \$0 to determine the
         # path of the real executable, they will execute our wrapper for the binary
         # (which will set up the library paths) instead of the original binary.
         # So we do not need to adjust any library paths here.
         source $(bin4ti2) "\$@"
         """)
      chmod(scriptwrapper, 0o755)

      for bin in readdir(bindir4ti2)
         if isfile(joinpath(bindir4ti2, bin))
            # a shell script starts with '#!'
            if read(joinpath(bindir4ti2, bin), 2) == UInt8['#', '!']
               symlink("4ti2_script_wrapper", bindir(bin))
            else
               symlink("4ti2_wrapper", bindir(bin))
            end
         end
      end
   end

   bind_artifact!(mutable_artifacts_toml,
      polymake_tree,
      polymake_tree_hash;
      force=true
   )

   # Point polymake to our custom tree
   return artifact_path(polymake_tree_hash)
end

max_stars_repo_name: plut/Yggdrasil

max_issues_repo_name: plut/Yggdrasil

max_forks_repo_name: plut/Yggdrasil


text:
abstract type Projector end
struct ThresholdProjector <: Projector
	β
	η
end

max_stars_repo_name: MKAbdElrahman/TopologyOptimizationTools

max_issues_repo_name: MKAbdElrahman/TopologyOptimizationTools

max_forks_repo_name: MKAbdElrahman/TopologyOptimizationTools


text:
@testset "Batchnorm" begin
    v = CUDA.rand(Float32, 2)
    m = CUDA.rand(Float32, 2, 5)

    @testset for training in (true, false), track_stats in (true, false)
        kws = (training=training, track_stats=track_stats)

        # Normal
        NNlibCUDA.batchnorm(v, v, m, v, v, 1.0; kws...)
        NNlibCUDA.∇batchnorm(v, v, m, m, v, v, 1.0; kws...)

        # No affine
        NNlibCUDA.batchnorm(nothing, nothing, m, v, v, 1.0; kws...)
        NNlibCUDA.∇batchnorm(nothing, nothing, m, m, v, v, 1.0; kws...)

        # No tracking
        NNlibCUDA.batchnorm(v, v, m, nothing, nothing, 1.0; kws...)
        NNlibCUDA.∇batchnorm(v, v, m, m, nothing, nothing, 1.0; kws...)

        # Both or neither tracked or affine params must be set
        for (α, β) in ((v, nothing), (nothing, v))
            @test_throws MethodError NNlibCUDA.batchnorm(α, β, m, v, v, 1.0; kws...)
            @test_throws MethodError NNlibCUDA.∇batchnorm(α, β, m, m, v, v, 1.0; kws...)
            @test_throws ArgumentError NNlibCUDA.batchnorm(v, v, m, α, β, 1.0; kws...)
        end
    end
end


text:
import HDF5

const HDF5_TS_ROOT_PATH = "time_series"
const TIME_SERIES_DATA_FORMAT_VERSION = "1.0.1"
const TIME_SERIES_VERSION_KEY = "data_format_version"
const COMPONENT_REFERENCES_KEY = "component_references"

"""
Stores all time series data in an HDF5 file.

The file used is assumed to be temporary and will be automatically deleted when there are
no more references to the storage object.
"""
mutable struct Hdf5TimeSeriesStorage <: TimeSeriesStorage
    file_path::String
    read_only::Bool
    compression::CompressionSettings
end

"""
Constructs Hdf5TimeSeriesStorage by creating a temp file.
"""
function Hdf5TimeSeriesStorage()
    return Hdf5TimeSeriesStorage(true)
end

"""
Constructs Hdf5TimeSeriesStorage.

# Arguments
- `create_file::Bool`: create new file
- `filename=nothing`: if nothing, create a temp file, else use this name.
- `directory=nothing`: if set and filename is nothing, create a temp file in this
   directory. Use tempdir() if not set. This should be set if the time series data is larger
   than the tmp filesystem can hold.
- `read_only = false`: If true, don't allow changes to the file. Allows simultaneous read
   access.
"""
function Hdf5TimeSeriesStorage(
    create_file::Bool;
    filename = nothing,
    directory = nothing,
    read_only = false,
    compression = CompressionSettings(),
)
    if create_file
        if isnothing(filename)
            if isnothing(directory)
                directory = tempdir()
            end
            filename, io = mktemp(directory)
            close(io)
        end

        storage = Hdf5TimeSeriesStorage(filename, read_only, compression)
        _make_file(storage)
    else
        storage = Hdf5TimeSeriesStorage(filename, read_only, compression)
    end

    @debug "Constructed new Hdf5TimeSeriesStorage" _group = LOG_GROUP_TIME_SERIES storage.file_path read_only compression

    return storage
end

"""
Constructs Hdf5TimeSeriesStorage from an existing file.
"""
function from_file(
    ::Type{Hdf5TimeSeriesStorage},
    filename::AbstractString;
    read_only = false,
    directory = nothing,
)
    if !isfile(filename)
        error("time series storage $filename does not exist")
    end

    if read_only
        file_path = abspath(filename)
    else
        parent = isnothing(directory) ? tempdir() : directory
        file_path, io = mktemp(parent)
        close(io)
        copy_file(filename, file_path)
    end

    storage = Hdf5TimeSeriesStorage(false; filename = file_path, read_only = read_only)
    if !read_only
        version = read_data_format_version(storage)
        if version == "1.0.0"
            _convert_from_1_0_0!(storage)
        end
        _deserialize_compression_settings!(storage)
    end

    @info "Loaded time series from storage file existing=$filename new=$(storage.file_path) compression=$(storage.compression)"
    return storage
end

"""
Copy the time series data to a new file. This should get called when the system is
undergoing a deepcopy.

# Arguments
- `storage::Hdf5TimeSeriesStorage`: storage instance
- `directory::String`: If nothing, use tempdir
"""
function copy_to_new_file!(storage::Hdf5TimeSeriesStorage, directory = nothing)
    if directory === nothing
        directory = tempdir()
    end

    # If we ever choose to keep the HDF5 file open then this will break.
    # Any open buffers will need to be flushed.
    filename, io = mktemp(directory)
    close(io)
    copy_file(get_file_path(storage), filename)
    storage.file_path = filename
end

get_compression_settings(storage::Hdf5TimeSeriesStorage) = storage.compression

get_file_path(storage::Hdf5TimeSeriesStorage) = storage.file_path

function read_data_format_version(storage::Hdf5TimeSeriesStorage)
    HDF5.h5open(storage.file_path, "r") do file
        root = _get_root(storage, file)
        if !haskey(HDF5.attributes(root), TIME_SERIES_VERSION_KEY)
            return "1.0.0"
        end
        return HDF5.read(HDF5.attributes(root)[TIME_SERIES_VERSION_KEY])
    end
end

function serialize_time_series!(
    storage::Hdf5TimeSeriesStorage,
    component_uuid::UUIDs.UUID,
    name::AbstractString,
    ts::TimeSeriesData,
)
    check_read_only(storage)
    uuid = string(get_uuid(ts))
    component_name = make_component_name(component_uuid, name)

    HDF5.h5open(storage.file_path, "r+") do file
        root = _get_root(storage, file)
        if !haskey(root, uuid)
            HDF5.create_group(root, uuid)
            path = root[uuid]
            # Create a group to store component references as attributes.
            # Use this instead of this time series' group or the dataset so that
            # the only attributes are component references.
            component_refs = HDF5.create_group(path, COMPONENT_REFERENCES_KEY)
            data = get_array_for_hdf(ts)
            settings = storage.compression
            if settings.enabled
                if settings.type == CompressionTypes.BLOSC
                    path["data", blosc = settings.level] = data
                elseif settings.type == CompressionTypes.DEFLATE
                    if settings.shuffle
                        path["data", shuffle = (), deflate = settings.level] = data
                    else
                        path["data", deflate = settings.level] = data
                    end
                else
                    error("not implemented for type=$(settings.type)")
                end
            else
                path["data"] = data
            end
            _write_time_series_attributes!(storage, ts, path)
            @debug "Create new time series entry." _group = LOG_GROUP_TIME_SERIES uuid component_uuid name
        else
            component_refs = root[uuid][COMPONENT_REFERENCES_KEY]
            @debug "Add reference to existing time series entry." _group =
                LOG_GROUP_TIME_SERIES uuid component_uuid name
        end
        HDF5.attributes(component_refs)[component_name] = true
    end

    return
end

"""
Return a String for the data type of the forecast data, this implementation avoids the use of `eval` on arbitrary code stored in HDF dataset.
"""
function get_data_type(ts::TimeSeriesData)
    data_type = eltype_data(ts)
    if data_type <: CONSTANT
        return "CONSTANT"
    elseif data_type == POLYNOMIAL
        return "POLYNOMIAL"
    elseif data_type == PWL
        return "PWL"
    elseif data_type <: Integer
        # We currently don't convert integers stored in TimeSeries.TimeArrays to floats.
        # This is a workaround.
        return "CONSTANT"
    else
        error("$data_type is not supported in forecast data")
    end
end

function _write_time_series_attributes!(
    storage::Hdf5TimeSeriesStorage,
    ts::T,
    path,
) where {T <: StaticTimeSeries}
    _write_time_series_attributes_common!(storage, ts, path)
end

function _write_time_series_attributes!(
    storage::Hdf5TimeSeriesStorage,
    ts::T,
    path,
) where {T <: Forecast}
    _write_time_series_attributes_common!(storage, ts, path)
    interval = get_interval(ts)
    HDF5.attributes(path)["interval"] = time_period_conversion(interval).value
end

function _write_time_series_attributes_common!(storage::Hdf5TimeSeriesStorage, ts, path)
    initial_timestamp = Dates.datetime2epochms(get_initial_timestamp(ts))
    resolution = get_resolution(ts)
    data_type = get_data_type(ts)
    HDF5.attributes(path)["module"] = string(parentmodule(typeof(ts)))
    HDF5.attributes(path)["type"] = string(nameof(typeof(ts)))
    HDF5.attributes(path)["initial_timestamp"] = initial_timestamp
    HDF5.attributes(path)["resolution"] = time_period_conversion(resolution).value
    HDF5.attributes(path)["data_type"] = data_type
end

function _read_time_series_attributes(
    storage::Hdf5TimeSeriesStorage,
    path,
    rows,
    ::Type{T},
) where {T <: StaticTimeSeries}
    return _read_time_series_attributes_common(storage, path, rows)
end

function _read_time_series_attributes(
    storage::Hdf5TimeSeriesStorage,
    path,
    rows,
    ::Type{T},
) where {T <: Forecast}
    data = _read_time_series_attributes_common(storage, path, rows)
    data["interval"] = Dates.Millisecond(HDF5.read(HDF5.attributes(path)["interval"]))
    return data
end

const _TYPE_DICT = Dict("CONSTANT" => CONSTANT, "POLYNOMIAL" => POLYNOMIAL, "PWL" => PWL)

function _read_time_series_attributes_common(storage::Hdf5TimeSeriesStorage, path, rows)
    initial_timestamp =
        Dates.epochms2datetime(HDF5.read(HDF5.attributes(path)["initial_timestamp"]))
    resolution = Dates.Millisecond(HDF5.read(HDF5.attributes(path)["resolution"]))
    data_type = _TYPE_DICT[HDF5.read(HDF5.attributes(path)["data_type"])]
    return Dict(
        "type" => _read_time_series_type(path),
        "initial_timestamp" => initial_timestamp,
        "resolution" => resolution,
        "dataset_size" => size(path["data"]),
        "start_time" => initial_timestamp + resolution * (rows.start - 1),
        "data_type" => data_type,
    )
end

function _read_time_series_type(path)
    module_str = HDF5.read(HDF5.attributes(path)["module"])
    type_str = HDF5.read(HDF5.attributes(path)["type"])
    return get_type_from_strings(module_str, type_str)
end

function add_time_series_reference!(
    storage::Hdf5TimeSeriesStorage,
    component_uuid::UUIDs.UUID,
    name::AbstractString,
    ts_uuid::UUIDs.UUID,
)
    check_read_only(storage)
    uuid = string(ts_uuid)
    component_name = make_component_name(component_uuid, name)
    HDF5.h5open(storage.file_path, "r+") do file
        root = _get_root(storage, file)
        path = root[uuid][COMPONENT_REFERENCES_KEY]

        # It's possible that this is overly restrictive, but as of now there is not a good
        # reason for a caller to add a reference multiple times. This should be a bug.
        @assert !haskey(HDF5.attributes(path), component_name) "There is already a reference to $component_name for time series $ts_uuid"

        HDF5.attributes(path)[component_name] = true
        @debug "Add reference to existing time series entry." _group = LOG_GROUP_TIME_SERIES uuid component_uuid name
    end
end

# TODO: This needs to change if we want to directly convert Hdf5TimeSeriesStorage to
# InMemoryTimeSeriesStorage, which is currently not supported at System deserialization.
function iterate_time_series(storage::Hdf5TimeSeriesStorage)
    Channel() do channel
        HDF5.h5open(storage.file_path, "r") do file
            root = _get_root(storage, file)
            for uuid_group in root
                uuid_path = HDF5.name(uuid_group)
                uuid_str = uuid_path[(findlast("/", uuid_path).start + 1):end]
                uuid = UUIDs.UUID(uuid_str)

                data = uuid_group["data"][:]
                attributes = Dict()
                for name in keys(HDF5.attributes(uuid_group))
                    attributes[name] = HDF5.read(HDF5.attributes(uuid_group)[name])
                end
                refs = uuid_group[COMPONENT_REFERENCES_KEY]
                for ref in keys(HDF5.attributes(refs))
                    component, name = deserialize_component_name(ref)
                    put!(channel, (component, name, data, attributes))
                end
            end
        end
    end
end

#=
# This could be used if we deserialize the type directly from HDF.
function _make_rows_columns(dataset, ::Type{T}) where T <: StaticTimeSeries
    rows = UnitRange(1, size(dataset)[1])
    columns = UnitRange(1, 1)
    return (rows, columns)
end

function _make_rows_columns(dataset, ::Type{T}) where T <: Forecast
    rows = UnitRange(1, size(dataset)[1])
    columns = UnitRange(1, size(dataset)[2])
    return (rows, columns)
end
=#

function remove_time_series!(
    storage::Hdf5TimeSeriesStorage,
    uuid::UUIDs.UUID,
    component_uuid::UUIDs.UUID,
    name::AbstractString,
)
    check_read_only(storage)
    HDF5.h5open(storage.file_path, "r+") do file
        root = _get_root(storage, file)
        path = _get_time_series_path(root, uuid)
        components = path[COMPONENT_REFERENCES_KEY]
        HDF5.delete_attribute(components, make_component_name(component_uuid, name))
        if isempty(keys(HDF5.attributes(components)))
            @debug "$path has no more references; delete it." _group = LOG_GROUP_TIME_SERIES
            HDF5.delete_object(path)
        end
    end
end

function deserialize_time_series(
    ::Type{T},
    storage::Hdf5TimeSeriesStorage,
    ts_metadata::TimeSeriesMetadata,
    rows::UnitRange,
    columns::UnitRange,
) where {T <: StaticTimeSeries}
    # Note that all range checks must occur at a higher level.
    return HDF5.h5open(storage.file_path, "r") do file
        root = _get_root(storage, file)
        uuid = get_time_series_uuid(ts_metadata)
        path = _get_time_series_path(root, uuid)
        attributes = _read_time_series_attributes(storage, path, rows, T)
        @assert_op attributes["type"] == T
        @debug "deserializing a StaticTimeSeries" _group = LOG_GROUP_TIME_SERIES T
        data_type = attributes["data_type"]
        data = get_hdf_array(path["data"], data_type, rows)
        return T(
            ts_metadata,
            TimeSeries.TimeArray(
                range(
                    attributes["start_time"];
                    length = length(rows),
                    step = attributes["resolution"],
                ),
                data,
            ),
        )
    end
end

function deserialize_time_series(
    ::Type{T},
    storage::Hdf5TimeSeriesStorage,
    ts_metadata::TimeSeriesMetadata,
    rows::UnitRange,
    columns::UnitRange,
) where {T <: AbstractDeterministic}
    # Note that all range checks must occur at a higher level.
    return HDF5.h5open(storage.file_path, "r") do file
        root = _get_root(storage, file)
        uuid = get_time_series_uuid(ts_metadata)
        path = _get_time_series_path(root, uuid)
        actual_type = _read_time_series_type(path)
        if actual_type == SingleTimeSeries
            last_index = size(path["data"])[1]
            return deserialize_deterministic_from_single_time_series(
                storage,
                ts_metadata,
                rows,
                columns,
                last_index,
            )
        end

        attributes = _read_time_series_attributes(storage, path, rows, T)
        @assert actual_type <: T "actual_type = $actual_type T = $T"
        @debug "deserializing a Forecast" _group = LOG_GROUP_TIME_SERIES T
        data_type = attributes["data_type"]
        data = get_hdf_array(path["data"], data_type, attributes, rows, columns)
        new_ts = actual_type(ts_metadata, data)
    end
end

function get_hdf_array(
    dataset,
    ::Type{<:CONSTANT},
    attributes::Dict{String, Any},
    rows::UnitRange{Int},
    columns::UnitRange{Int},
)
    data = SortedDict{Dates.DateTime, Vector{Float64}}()
    initial_timestamp = attributes["start_time"]
    interval = attributes["interval"]
    start_time = initial_timestamp + interval * (columns.start - 1)
    if length(columns) == 1
        data[start_time] = dataset[rows, columns.start]
    else
        data_read = dataset[rows, columns]
        for (i, it) in
            enumerate(range(start_time; length = length(columns), step = interval))
            data[it] = @view data_read[1:length(rows), i]
        end
    end
    return data
end

function get_hdf_array(
    dataset,
    type::Type{POLYNOMIAL},
    attributes::Dict{String, Any},
    rows::UnitRange{Int},
    columns::UnitRange{Int},
)
    data = SortedDict{Dates.DateTime, Vector{POLYNOMIAL}}()
    initial_timestamp = attributes["start_time"]
    interval = attributes["interval"]
    start_time = initial_timestamp + interval * (columns.start - 1)
    if length(columns) == 1
        data[start_time] = retransform_hdf_array(dataset[rows, columns.start, :], type)
    else
        data_read = retransform_hdf_array(dataset[rows, columns, :], type)
        for (i, it) in
            enumerate(range(start_time; length = length(columns), step = interval))
            data[it] = @view data_read[1:length(rows), i]
        end
    end
    return data
end

function get_hdf_array(
    dataset,
    type::Type{PWL},
    attributes::Dict{String, Any},
    rows::UnitRange{Int},
    columns::UnitRange{Int},
)
    data = SortedDict{Dates.DateTime, Vector{PWL}}()
    initial_timestamp = attributes["start_time"]
    interval = attributes["interval"]
    start_time = initial_timestamp + interval * (columns.start - 1)
    if length(columns) == 1
        data[start_time] = retransform_hdf_array(dataset[rows, columns.start, :, :], type)
    else
        data_read = retransform_hdf_array(dataset[rows, columns, :, :], type)
        for (i, it) in
            enumerate(range(start_time; length = length(columns), step = interval))
            data[it] = @view data_read[1:length(rows), i]
        end
    end
    return data
end

function get_hdf_array(dataset, type::Type{<:CONSTANT}, rows::UnitRange{Int})
    data = retransform_hdf_array(dataset[rows], type)
    return data
end

function get_hdf_array(dataset, type::Type{POLYNOMIAL}, rows::UnitRange{Int})
    data = retransform_hdf_array(dataset[rows, :, :], type)
    return data
end

function get_hdf_array(dataset, type::Type{PWL}, rows::UnitRange{Int})
    data = retransform_hdf_array(dataset[rows, :, :, :], type)
    return data
end

function retransform_hdf_array(data::Array, ::Type{<:CONSTANT})
    return data
end

function retransform_hdf_array(data::Array, T::Type{POLYNOMIAL})
    row, column, tuple_length = get_data_dims(data, T)
    if isnothing(column)
        t_data = Array{POLYNOMIAL}(undef, row)
        for r in 1:row
            t_data[r] = tuple(data[r, 1:tuple_length]...)
        end
    else
        t_data = Array{POLYNOMIAL}(undef, row, column)
        for r in 1:row, c in 1:column
            t_data[r, c] = tuple(data[r, c, 1:tuple_length]...)
        end
    end
    return t_data
end

function retransform_hdf_array(data::Array, T::Type{PWL})
    row, column, tuple_length, array_length = get_data_dims(data, T)
    if isnothing(column)
        t_data = Array{PWL}(undef, row)
        for r in 1:row
            tuple_array = Array{POLYNOMIAL}(undef, array_length)
            for l in 1:array_length
                tuple_array[l] = tuple(data[r, 1:tuple_length, l]...)
            end
            t_data[r] = tuple_array
        end
    else
        t_data = Array{PWL}(undef, row, column)
        for r in 1:row, c in 1:column
            tuple_array = Array{POLYNOMIAL}(undef, array_length)
            for l in 1:array_length
                tuple_array[l] = tuple(data[r, c, 1:tuple_length, l]...)
            end
            t_data[r, c] = tuple_array
        end
    end
    return t_data
end

function get_data_dims(data::Array, ::Type{POLYNOMIAL})
    if length(size(data)) == 2
        row, tuple_length = size(data)
        return (row, nothing, tuple_length)
    elseif length(size(data)) == 3
        return size(data)
    else
        error("Hdf data array is $(length(size(data)))-D array, expected 2-D or 3-D array.")
    end
end

function get_data_dims(data::Array, ::Type{PWL})
    if length(size(data)) == 3
        row, tuple_length, array_length = size(data)
        return (row, nothing, tuple_length, array_length)
    elseif length(size(data)) == 4
        return size(data)
    else
        error("Hdf data array is $(length(size(data)))-D array, expected 3-D or 4-D array.")
    end
end

function deserialize_time_series(
    ::Type{T},
    storage::Hdf5TimeSeriesStorage,
    ts_metadata::TimeSeriesMetadata,
    rows::UnitRange,
    columns::UnitRange,
) where {T <: Probabilistic}
    # Note that all range checks must occur at a higher level.
    total_percentiles = length(get_percentiles(ts_metadata))

    return HDF5.h5open(storage.file_path, "r") do file
        root = _get_root(storage, file)
        uuid = get_time_series_uuid(ts_metadata)
        path = _get_time_series_path(root, uuid)
        attributes = _read_time_series_attributes(storage, path, rows, T)
        @assert_op attributes["type"] == T
        @assert_op length(attributes["dataset_size"]) == 3
        @debug "deserializing a Forecast" _group = LOG_GROUP_TIME_SERIES T
        data = SortedDict{Dates.DateTime, Matrix{attributes["data_type"]}}()
        initial_timestamp = attributes["start_time"]
        interval = attributes["interval"]
        start_time = initial_timestamp + interval * (first(columns) - 1)
        if length(columns) == 1
            data[start_time] =
                transpose(path["data"][1:total_percentiles, rows, first(columns)])
        else
            data_read = PermutedDimsArray(
                path["data"][1:total_percentiles, rows, columns],
                [3, 2, 1],
            )
            for (i, it) in enumerate(
                range(start_time; length = length(columns), step = attributes["interval"]),
            )
                data[it] = @view data_read[i, 1:length(rows), 1:total_percentiles]
            end
        end

        new_ts = T(ts_metadata, data)
    end
end

function deserialize_time_series(
    ::Type{T},
    storage::Hdf5TimeSeriesStorage,
    ts_metadata::TimeSeriesMetadata,
    rows::UnitRange,
    columns::UnitRange,
) where {T <: Scenarios}
    # Note that all range checks must occur at a higher level.
    total_scenarios = get_scenario_count(ts_metadata)

    return HDF5.h5open(storage.file_path, "r") do file
        root = _get_root(storage, file)
        uuid = get_time_series_uuid(ts_metadata)
        path = _get_time_series_path(root, uuid)
        attributes = _read_time_series_attributes(storage, path, rows, T)
        @assert_op attributes["type"] == T
        @assert_op length(attributes["dataset_size"]) == 3
        @debug "deserializing a Forecast" _group = LOG_GROUP_TIME_SERIES T
        data = SortedDict{Dates.DateTime, Matrix{attributes["data_type"]}}()
        initial_timestamp = attributes["start_time"]
        interval = attributes["interval"]
        start_time = initial_timestamp + interval * (first(columns) - 1)
        if length(columns) == 1
            data[start_time] =
                transpose(path["data"][1:total_scenarios, rows, first(columns)])
        else
            data_read =
                PermutedDimsArray(path["data"][1:total_scenarios, rows, columns], [3, 2, 1])
            for (i, it) in enumerate(
                range(start_time; length = length(columns), step = attributes["interval"]),
            )
                data[it] = @view data_read[i, 1:length(rows), 1:total_scenarios]
            end
        end

        new_ts = T(ts_metadata, data)
    end
end

function clear_time_series!(storage::Hdf5TimeSeriesStorage)
    check_read_only(storage)
    # Re-create the file. HDF5 will not actually free up the deleted space until h5repack
    # is run on the file.
    _make_file(storage)
    @info "Cleared all time series."
end

function get_num_time_series(storage::Hdf5TimeSeriesStorage)
    num = 0

    HDF5.h5open(storage.file_path, "r") do file
        root = _get_root(storage, file)
        for component in root
            num += 1
        end
    end

    return num
end

function _make_file(storage::Hdf5TimeSeriesStorage)
    HDF5.h5open(storage.file_path, "w") do file
        root = HDF5.create_group(file, HDF5_TS_ROOT_PATH)
        HDF5.attributes(root)[TIME_SERIES_VERSION_KEY] = TIME_SERIES_DATA_FORMAT_VERSION
        _serialize_compression_settings(storage, root)
    end
end

function _serialize_compression_settings(storage::Hdf5TimeSeriesStorage, root)
    HDF5.attributes(root)["compression_enabled"] = storage.compression.enabled
    HDF5.attributes(root)["compression_type"] = string(storage.compression.type)
    HDF5.attributes(root)["compression_level"] = storage.compression.level
    HDF5.attributes(root)["compression_shuffle"] = storage.compression.shuffle
end

function _deserialize_compression_settings!(storage::Hdf5TimeSeriesStorage)
    HDF5.h5open(storage.file_path, "r+") do file
        root = _get_root(storage, file)
        storage.compression = CompressionSettings(
            enabled = HDF5.read(HDF5.attributes(root)["compression_enabled"]),
            type = CompressionTypes(HDF5.read(HDF5.attributes(root)["compression_type"])),
            level = HDF5.read(HDF5.attributes(root)["compression_level"]),
            shuffle = HDF5.read(HDF5.attributes(root)["compression_shuffle"]),
        )
    end
    return
end

_get_root(storage::Hdf5TimeSeriesStorage, file) = file[HDF5_TS_ROOT_PATH]

function _get_time_series_path(root::HDF5.Group, uuid::UUIDs.UUID)
    uuid_str = string(uuid)
    if !haskey(root, uuid_str)
        throw(ArgumentError("UUID $uuid_str does not exist"))
    end

    return root[uuid_str]
end

function check_read_only(storage::Hdf5TimeSeriesStorage)
    if storage.read_only
        error("Operation not permitted; this time series file is read-only")
    end
end

is_read_only(storage::Hdf5TimeSeriesStorage) = storage.read_only

function compare_values(
    x::Hdf5TimeSeriesStorage,
    y::Hdf5TimeSeriesStorage;
    compare_uuids = false,
)
    item_x = sort!(collect(iterate_time_series(x)), by = z -> z[1])
    item_y = sort!(collect(iterate_time_series(y)), by = z -> z[1])
    if length(item_x) != length(item_y)
        @error "lengths don't match" length(item_x) length(item_y)
        return false
    end

    for ((uuid_x, name_x, data_x, attrs_x), (uuid_y, name_y, data_y, attrs_y)) in
        zip(item_x, item_y)
        if uuid_x != uuid_y
            @error "component UUIDs don't match" uuid_x uuid_y
            return false
        end
        if name_x != name_y
            @error "names don't match" name_x name_y
            return false
        end
        if data_x != data_y
            @error "data doesn't match" data_x data_y
            return false
        end
        if sort!(collect(keys(attrs_x))) != sort!(collect(keys(attrs_y)))
            @error "attr keys don't match" attrs_x attrs_y
        end
        if collect(values(attrs_x)) != collect(values(attrs_y))
            @error "attr values don't match" attrs_x attrs_y
        end
    end

    return true
end

function _convert_from_1_0_0!(storage::Hdf5TimeSeriesStorage)
    # 1.0.0 version did not support compression.
    # 1.0.0 stored component name/UUID pairs in a dataset.
    # That wasn't efficient if a user added many shared references.
    HDF5.h5open(storage.file_path, "r+") do file
        root = _get_root(storage, file)
        for uuid_group in root
            components = HDF5.create_group(uuid_group, COMPONENT_REFERENCES_KEY)
            component_names = uuid_group["components"][:]
            for name in component_names
                HDF5.attributes(components)[name] = true
            end
            HDF5.delete_object(uuid_group["components"])
        end

        HDF5.attributes(root)[TIME_SERIES_VERSION_KEY] = TIME_SERIES_DATA_FORMAT_VERSION
        compression = CompressionSettings()
        _serialize_compression_settings(storage, root)
    end

    @debug "Converted file from 1.0.0 format" _group = LOG_GROUP_TIME_SERIES
end


text:
module DVeTest
using Base.Test
using Quantum

xb = basis(:X,[1:5])
d = DiracVector([1:5], xb)
@test d==d''
@test d[:]==d.coeffs
@test d'*d == 55
@test d*d' == DiracMatrix(kron(d.coeffs, d.coeffs'), d.basis, d.basis')
@test isdual(d,d')
@test filterstates(s->label(s)%2==0, d) == DiracVector([2,4], basis(xb[2], xb[4]))
@test filtercoeffs(c->c%2==0, d) == DiracVector([2,4], basis(xb[2], xb[4]))
@test map(x->x*2, d)==d+d
@test xb[1]'*d == 1
@test d'*xb[1] == 1
@test norm(normalize(d)) < 1+1e-8 && norm(normalize(d)) > 1-1e-8
@test reduce(+,[i*xb[i]' for i=1:length(xb)]) == d'
@test (1/sqrt(2)*xb[1]) + (1/sqrt(2)*xb[2]) == dvec([1/sqrt(2), 1/sqrt(2)], basis(xb[1], xb[2]))
@test ket(:X,1)+ket(:X,1)+ket(:X,1) = dvec([3], basis(xb[1]))
@test ket(:X,1)+(bra(:S,"1")*ket(:X, 1)) * ket(:X,1)==dvec([1+1*bra(:S,"1")*ket(:X,1)], basis(ket(:X,1)))

end


text:
function value_of(::Type{PositionAngle}, arg0::JString)
    return jcall(PositionAngle, "valueOf", PositionAngle, (JString,), arg0)
end

function values(::Type{PositionAngle})
    return jcall(PositionAngle, "values", Vector{PositionAngle}, ())
end


text:
import CombinedParsers: ParserTypes
using CombinedParsers.Regexp

import CombinedParsers.Regexp: char, integer_base, escape_sequence, escaped_character
@testset "char" begin
    ##@test match(parse(char,with_options(Base.PCRE.CASELESS,"A")) =='a'
    @test parse(char,"A") == CharIn('A')
    @test parse(char,"\\^") == CharIn('^')
    @test tryparse(char,"^") === nothing
    @test parse(escape_sequence(),raw"\Q[].\E")=="[]."
    ##@btime _iterate(pattern,".")
    ## @btime _iterate(pattern,"\\N")
    @test parse(integer_base(8),"765")==501
    @test parse(integer_base(10),"765")==765
    @test parse(integer_base(16),"765")==1893
    @test parse(integer_base(16),"")==0
    @test parse(Sequence(v->Char(v[2]),"\\x{",integer_base(16),"}"),"\\x{10}") == '\x10'
    @test parse(Repeat(escaped_character),raw"\a\t\r\n") == collect("\a\t\r\n")
    @test parse(parser(parse(escaped_character,"\\o{100}")),"@")=='@'
    @test parse(parser(parse(escaped_character,"\\x10")),"\U0010")=='\x10'
    @test parse(parser(parse(escaped_character,"\\x{010}")),"\U0010")=='\x10'
    @test parse(escaped_character,"\\t")==('\t')
    @test parse(escaped_character,"\\x{0065}") == ('e')
    @test parse(integer_base(8,3,3),"100")==64
    match(r"\100","@")
    match(r"\x100","@")
    match(r"\o{100}","@")
    match(r"\x1","\U0010")
end

import CombinedParsers.Regexp: @test_pcre
@testset "char groups" begin
    @test_pcre "\\d+" "1123"
    @test_pcre "\\D+" "abcd"
    @test_pcre "\\R" "\n"
    @test_pcre "\\W" " "
end

import CombinedParsers.Regexp: bracket_char, bracket, pcre_options
@testset "brackets" begin
    @test parse(bracket_char,"a")=='a'
    @test parse(bracket,"[a-z]").sets==CharIn('a':'z').sets
    @test parse(bracket,with_options(Base.PCRE.EXTENDED_MORE,"[a-z ]")).sets==CharIn('a':'z').sets
    @test parse(bracket,"[\x3f-\x5F]").sets==CharIn('?':'_').sets
    @test parse(parser(']'),"]")==']'
    @test parse(parse(bracket,"[]abc]"),"]")==']'
    @test parse(parse(bracket,"[]-_abc-]"),"]")==']'
    @test match(re"[-]abc]","a") === nothing
    @test parse(re"[]abc-]","-")=='-'
    @test parse(re" [a- z]"xxx,"b")=='b'
    @test tryparse(re" [a- z]"xx,"b")==nothing
    @test tryparse(re"[a- z]","b")==nothing
    @test parse(re"[a- z]"xx,"b")=='b'
    ##@test match(r" [a-z]"xxx,"b")==match(re" [a-z]"xxx,"b")
    @test result_type(bracket) == CombinedParser
    @test parse(bracket,"[a\\E]") == CharIn('a')
    ##parse(bracket,with_options("i","[^]a]"))
    @test_throws ArgumentError Regcomb("[")
end

parse(re"\v(?#ss=)","\n")

import CombinedParsers.Regexp: option_sequences,skip_whitespace_on, ParserWithCaptures
@testset "internal options" begin
    @test match(r"(?xx)[a- x]","b")==match(re"(?xx)[a- x]","b")
    @test parse(pcre_options,"i")==Base.PCRE.CASELESS
    @test parse(parse(option_sequences,"a|b(?i)a"),"bA")==('b','A')
    @test re" a"x == CharIn('a')

@testset "comments" begin
    ##@test parse(comment_par,"(?# comment)") == with_log("comment",Always())
    @test match(re"(?# comment)a","a").match == "a"
    ## match(re"(?# comment)a","a").match == "a"

    ## parse(seq(parse(skip_whitespace_and_comments,with_options("x"," (?#xxx) (?#yyy) "))...),"a")
## push!(pattern,lazy(opt(comment_par)));

## parse_all(lazy(opt(comment_par)),"(?#a)")

end

end

import CombinedParsers: Repeat_max
import CombinedParsers.Regexp: quantified, repetition, sequence
@testset "repetitions and optional" begin
    @test parse(Optional('a'),"b")===missing
    @test parse(Optional('a'),"a")==='a'
    @test parse(repetition,"{1,}") == 1:Repeat_max
    @test parse(repetition,"*") == 0:Repeat_max
    @test parse(repetition,"+") == 1:Repeat_max
    @test parse(repetition,"{3}") == 3:3
    @test parse(repetition,"?") == 0:1
    @test match(parse(quantified,"a*"),"aaab").match=="aaa"
    @test parse(parse(quantified,"a*?"),"aa")==Char[]
    ## lazy support
    @test_pcre "a*?(abc)?" "abc"
    ## lazy rep
    @test_pcre "a(?:b|(c|e){1,2}?|d)+?(.)" "ace"
    @test parse(parse(quantified,"a*"),"aa")==Char['a','a']
    @test parse(quantified,"a?")|>regex_string == "a?"
    parse(quantified,"a{3}")
    parse(quantified,with_options(Base.PCRE.EXTENDED,"a {3}"))
    parse(parse(sequence,"ab*"), "abbb")
    pp=parse(sequence,"ab*")
    ##@btime _iterate(pp, "abbb")
    ##@btime match(r"ab*","abbb")
    @test parse(quantified,with_options(Base.PCRE.CASELESS,"a*?")) == Lazy(Repeat(CharIn('a','A')))
end

import CombinedParsers.Regexp: alternation
@testset "alternations" begin
    parse(alternation,"[-]abc\\]")
    parse(alternation,"a|b")
    parse(alternation,
          with_options(Base.PCRE.EXTENDED|Base.PCRE.EXTENDED_MORE,
                       "a {3}bc | d | [a - e]# comment?"))
    @test parse(Sequence(2,AtStart(),alternation,AtEnd()),"")==Sequence()
end

import CombinedParsers.Regexp: captured, subpattern, atomic_group, backreference, subroutine
@testset "sequences, captures" begin
    match(parse(captured,"(ab)"),"ab")
    parse(captured,"()")
    pp = parse(captured,"(?ab)")
    parse(parse(captured,"(ab)"), "ab")
    parse(parse(subpattern,"(?:ab)"), "ab")
    @test parse(atomic_group,"(?>ab*)")|>regex_string == "(?>ab*)"
    @test parse(backreference,"(?P=ab)").name == :ab
    @test parse(backreference,"\\1").index == 1
    @test parse(backreference,"\\g-1").index == -1
    @test parse(subroutine,"(?-1)").delta == Symbol("-")
    @test_pcre "((((a)(b))))(c)" "abc" true
    @test_pcre "(a)|(d)" "a" true
    @test_pcre "(a)|(d)" "d" true
    @test_pcre "(1a(2b?)*)*0" "1a1a21a2b22b0" true
    @test_pcre "(1a(2b?)*)*0" "1a1a21a2b22b0" true
    @test_pcre "(ab)*c" "ababc" true
    @test_pcre "^(ab)*c\$" "ababc" true
    @test_pcre "(?<ab>foo)\\1" "foofoo"
    re"(?<ab>foo)(?P=ab)"
    @test_pcre "(?<ab>foo)(?P=ab)" "foofoo" true
    re"(?<ab>foo|bar)(?1)"
    @test_pcre "(?<ab>foo|bar)(?1)" "foobar" true
    @test_pcre "(ab)(?<ab>foo)\\g-1" "abfoofoo" true
    ## recursive back references
    @test_pcre "(a|b\\1)+" "aaa" true
    @test_pcre "(a|b\\1)+" "aba" true
    @test_pcre "(a|b\\1)+" "ababbaa" true
    @test_pcre "(ab)(?<ab>foo)\\g-1\\g-2" "abfoofooab" true
    @test_pcre "(se|respo)nse and ((?1)nse)" "sense and response" true
    re"(?ab)(?#comment)"
    @test_pcre "^ab*(?<ab>c)" "ac"
    @test_pcre raw"\Q \Ea" " a" true "x"
    @test_pcre "a*abc?xyz+pqr{3}ab{2,}xy{4,5}pq{0,6}AB{0,}zz" "abcxyzpqrrrabbxyyyypqAzz" true
end

@test_pcre raw"\\\\\"" "\\\"" true
@test_pcre raw"\"" "\"" true

@testset "look around" begin
    @test_pcre "ab*(?<=ab)c" "abc" true
end

@test_pcre "(ab|a|b)+c" "abbabc" true

@test regex_string(AtEnd())=="\$"

@test_pcre "^(a|b|abc)+c\$" "abcbabc"

@test_pcre "^aaa(?<!c)b" "aaab"

m = match(PositiveLookahead("a"), "aaab")
@test m.start == 1
@test m.stop == 1
@test m.state == CombinedParsers.MatchState()
@test get(m) == "a"


text:
function getRobject(file::String, name::String)
    R"""
    library(nodiv)
    load($file)
    """
    getRobject(name)
end

function getRobject(name::String)
    R"""
    obj = get($name)
    if(!inherits(obj, "distrib_data"))
        stop(paste(name, "is not a distrib_data object!"))

    obj$shape <- NULL
    obj$sitestat = obj$coords@data
    obj$coords = sp::coordinates(obj$coords)
    """
    rcopy(R"obj")
end

function SpatialEcology.Assemblage(rdict::Dict{T, S}) where {T <: Symbol, S}

    cd_type = rdict[:type] == "grid" ? griddata : pointdata #this code only works as long as there are only those two types
    Assemblage(rdict[:comm], rdict[:coords] ,
        cdtype = cd_type, traits = rdict[:species_stats],
        sitestats = rdict[:sitestat])
end


text:
#Define forces
include("Frames.jl")
include("Transformables.jl")
abstract type AbstractForceMoment <: AbstractTransformable  end

mutable struct Force <: AbstractForceMoment
    Fx::Float64
    Fy::Float64
    Fz::Float64
    frame::Type{<:AbstractFrame}
end

function (force::Force)(newval::Vector{Float64})
    force.Fx = newval[1]
    force.Fy = newval[2]
    force.Fz = newval[3]
    return force
end

function (force::Force)()
    return [force.Fx force.Fy force.Fz]
end

mutable struct Moment <: AbstractForceMoment
    Mx::Float64
    My::Float64
    Mz::Float64
    frame::Type{<:AbstractFrame}
end

function (moment::Moment)(newval::Vector{Float64})
    moment.Mx = newval[1]
    moment.My = newval[2]
    moment.Mz = newval[3]
    return moment
end

function (moment::Moment)()
    return [moment.Mx moment.My moment.Mz]
end

FrameOf(fm::Type{<:AbstractForceMoment}) = fm.frame


text:
export iauAticqn
"""
Quick CIRS to ICRS astrometric place transformation, given the star-
independent astrometry parameters plus a list of light-deflecting
bodies.

Use of this function is appropriate when efficiency is important and
where many star positions are all to be transformed for one date.
The star-independent astrometry parameters can be obtained by
calling one of the functions iauApci[13], iauApcg[13], iauApco[13]
or iauApcs[13].

If the only light-deflecting body to be taken into account is the
Sun, the iauAticq function can be used instead.

This function is part of the International Astronomical Union's
SOFA (Standards of Fundamental Astronomy) software collection.

Status:  support function.

Given:
   ri,di  double      CIRS RA,Dec (radians)
   astrom iauASTROM*  star-independent astrometry parameters:
      pmt    double       PM time interval (SSB, Julian years)
      eb     double[3]    SSB to observer (vector, au)
      eh     double[3]    Sun to observer (unit vector)
      em     double       distance from Sun to observer (au)
      v      double[3]    barycentric observer velocity (vector, c)
      bm1    double       sqrt(1-|v|^2): reciprocal of Lorenz factor
      bpn    double[3][3] bias-precession-nutation matrix
      along  double       longitude + s' (radians)
      xpl    double       polar motion xp wrt local meridian (radians)
      ypl    double       polar motion yp wrt local meridian (radians)
      sphi   double       sine of geodetic latitude
      cphi   double       cosine of geodetic latitude
      diurab double       magnitude of diurnal aberration vector
      eral   double       "local" Earth rotation angle (radians)
      refa   double       refraction constant A (radians)
      refb   double       refraction constant B (radians)
   n     int           number of bodies (Note 3)
   b     iauLDBODY[n] data for each of the n bodies (Notes 3,4):
      bm    double       mass of the body (solar masses, Note 5)
      dl    double       deflection limiter (Note 6)
      pv    [2][3]       barycentric PV of the body (au, au/day)

Returned:
   rc,dc  double     ICRS astrometric RA,Dec (radians)

Notes:

   1. Iterative techniques are used for the aberration and light
   deflection corrections so that the functions iauAticqn and
   iauAtciqn are accurate inverses; even at the edge of the Sun's
   disk the discrepancy is only about 1 nanoarcsecond.

   2. If the only light-deflecting body to be taken into account is the
   Sun, the iauAticq function can be used instead.

   3. The struct b contains n entries, one for each body to be
   considered.  If n = 0, no gravitational light deflection will be
   applied, not even for the Sun.

   4. The struct b should include an entry for the Sun as well as for
   any planet or other body to be taken into account.  The entries
   should be in the order in which the light passes the body.

   5. In the entry in the b struct for body i, the mass parameter
   b[i].bm can, as required, be adjusted in order to allow for such
   effects as quadrupole field.

   6. The deflection limiter parameter b[i].dl is phi^2/2, where phi is
   the angular separation (in radians) between star and body at
   which limiting is applied.  As phi shrinks below the chosen
   threshold, the deflection is artificially reduced, reaching zero
   for phi = 0.   Example values suitable for a terrestrial
   observer, together with masses, are as follows:

      body i     b[i].bm        b[i].dl

      Sun        1.0            6e-6
      Jupiter    0.00095435     3e-9
      Saturn     0.00028574     3e-10

   7. For efficiency, validation of the contents of the b array is
   omitted.  The supplied masses must be greater than zero, the
   position and velocity vectors must be right, and the deflection
   limiter greater than zero.

Called:
   iauS2c       spherical coordinates to unit vector
   iauTrxp      product of transpose of r-matrix and p-vector
   iauZp        zero p-vector
   iauAb        stellar aberration
   iauLdn       light deflection by n bodies
   iauC2s       p-vector to spherical
   iauAnp       normalize angle into range +/- pi

This revision:   2013 October 9

SOFA release 2018-01-30

Copyright (C) 2018 IAU SOFA Board.  See notes at end.
"""
# void iauAticqn(double ri, double di, iauASTROM *astrom,
#                int n, iauLDBODY b[], double *rc, double *dc)

function iauAticqn(ri::Real, di::Real, astrom::iauASTROM,
                   n::Int, b::Array{iauLDBODY, 1})

   # Allocate return value
   ref_astrom = Ref{iauASTROM}(astrom)
   ref_rc     = Ref{Float64}(0.0)
   ref_dc     = Ref{Float64}(0.0)
   ref_b      = Ref{iauLDBODY}(b[1])

   # println(b)

   # Copy b, and transpose pv field
   # for ld in b
   #    ld.pv = [ld.pv[1, 1] ld.pv[1, 3] ld.pv[2, 2];
   #             ld.pv[1, 2] ld.pv[2, 1] ld.pv[2, 3]]
   # end
   for ld in b
      ld.pv = ld.pv'
   end

   # println(b)

   ccall((:iauAticqn, libsofa_c), Cvoid,
            (Cdouble, Cdouble, Ref{iauASTROM}, Cint, Ref{iauLDBODY},
            Ref{Cdouble}, Ref{Cdouble}),
            convert(Float64, ri), convert(Float64, di),
            ref_astrom, convert(Int32, n), ref_b, ref_rc, ref_dc)

   # Reverse transpose operation since Julia is pass by reference
   for ld in b
      ld.pv = ld.pv'
   end

   return ref_rc[], ref_dc[]
end


text:
# idea: use a type parameterized on functions as the stack element
using DataStructures
import Base: convert, promote_rule, +, -, *, /, >, <, log, exp, <=, >=

type RDiff <: Real
    val::Float64
    adj::Float64
    RDiff(val, adj) = new(val, adj)
end
typealias Converts Union{Float64, Int, Irrational, Rational}
RDiff(x::Converts) = RDiff(x, 0.)

convert(::Type{RDiff}, x::Converts) = RDiff(x)
promote_rule(::Type{Bool}, ::Type{RDiff}) = RDiff
promote_rule(::Type{Float64}, ::Type{RDiff}) = RDiff
promote_rule{T<:Integer}(::Type{T}, ::Type{RDiff}) = RDiff

value(x::RDiff) = x.val
adjoint(x::RDiff) = x.adj
value(x::Array{RDiff}) = map(value, x)
adjoint(x::Array{RDiff}) = map(adjoint, x)
function incadjoint!{N}(x::Array{RDiff, N}, y::Array{Float64, N})
    for i in eachindex(x)
        x[i].adj += y[i]
    end
end
>(x::RDiff, y::RDiff) = x.val > y.val
<(x::RDiff, y::RDiff) = x.val < y.val
>=(x::RDiff, y::RDiff) = x.val >= y.val
<=(x::RDiff, y::RDiff) = x.val <= y.val

type RNode{F<:Function, T1<:Tuple, T2}
    op::F
    inputs::T1
    outputs::T2
end

const global varstack = Stack(RNode)

function grad(f, N)
    y = Array{RDiff}(N)

    function ∇f(x, storage)
        copy!(y, x)
        f(y)
        for node in varstack
            backprop(node)
        end
        clear!(varstack)
        for i in eachindex(storage)
            storage[i] = y[i].val
        end
    end

    return ∇f
end

function clear!(s::Stack)
    while !isempty(s)
        pop!(s)
    end
    nothing
end

backprop(n::RNode) = backprop(n.op, n.inputs..., n.outputs)

macro diff_rule(sig, code)
    fn = _fn_of_vals(sig)
    op = sig.args[1]
    op_args = sig.args[2:end]
    argsyms = Symbol[s.args[1] for s in sig.args[2:end]]
    ex = quote
        function $op($(op_args...))
            res = map(RDiff, $fn)
            n = RNode($op, tuple($(argsyms...)), res)
            push!(varstack, n)
            res
        end

        function backprop(::typeof($op), $(op_args...), res)
            $code
            nothing
        end
    end
    esc(ex)
end

_fn_of_vals(x) = x
function _fn_of_vals(ex::Expr)
    if ex.head == :(::)
        out = :(value($(ex.args[1])))
    else
        out = copy(ex)
        for i in eachindex(out.args)
            out.args[i] = _fn_of_vals(out.args[i])
        end
    end
    out
end

@diff_rule log(x::RDiff) begin
    x.adj += res.adj / x.val
end

@diff_rule exp(x::RDiff) begin
    x.adj += res.adj * x.val
end

@diff_rule +(x::RDiff, y::RDiff) begin
    x.adj += res.adj
    y.adj += res.adj
end

@diff_rule -(x::RDiff, y::RDiff) begin
    x.adj += res.adj
    y.adj += -res.adj
end

@diff_rule *(x::RDiff, y::RDiff) begin
    x.adj += res.adj * y.val
    y.adj += res.adj * x.val
end

@diff_rule /(x::RDiff, y::RDiff) begin
    x.adj += res.adj / y.val
    y.adj += -res.adj * x.val / y.val^2
end

@diff_rule -(x::RDiff) begin
    x.adj += -res.adj
end

@diff_rule *(x::Array{RDiff}, y::Array{RDiff}) begin
    incadjoint!(x, adjoint(res) * value(y)')
    incadjoint!(y, value(x)' * adjoint(res))
end

@diff_rule +(x::Array{RDiff}, y::Array{RDiff}) begin
    incadjoint!(x, adjoint(res))
    incadjoint!(y, adjoint(res))
end

#= TODO
@diff_rule map(f::Function, x::Array{RDiff}) begin
    incadjoint!(x, adjoint(res) .* grad_expr(f)(x))
end
=#


text:
function test_add_remove(storage::IS.TimeSeriesStorage)
    name = "component1"
    label = "val"
    component = IS.TestComponent(name, 5)
    ts = create_time_series_data()
    IS.add_time_series!(storage, IS.get_uuid(component), label, ts)

    ts_data = IS.get_time_series(storage, IS.get_uuid(ts))

    @test TimeSeries.timestamp(ts_data) == TimeSeries.timestamp(ts.data)
    @test TimeSeries.values(ts_data) == TimeSeries.values(ts.data)

    component2 = IS.TestComponent("component2", 6)
    IS.add_time_series!(storage, IS.get_uuid(component2), label, ts)

    IS.get_num_time_series(storage) == 2

    IS.remove_time_series!(storage, IS.get_uuid(ts), IS.get_uuid(component2), label)

    # There should still be one reference to the data.
    ts_data2 = IS.get_time_series(storage, IS.get_uuid(ts))
    @test ts_data2 isa TimeSeries.TimeArray

    IS.remove_time_series!(storage, IS.get_uuid(ts), IS.get_uuid(component), label)
    @test_throws ArgumentError IS.get_time_series(storage, IS.get_uuid(ts))
    IS.get_num_time_series(storage) == 0
end

function test_get_subset(storage::IS.TimeSeriesStorage)
    name = "component1"
    label = "val"
    component = IS.TestComponent(name, 1)
    ts = create_time_series_data()
    IS.add_time_series!(storage, IS.get_uuid(component), label, ts)
    ts_data = IS.get_time_series(storage, IS.get_uuid(ts))

    @test TimeSeries.timestamp(ts_data) == TimeSeries.timestamp(ts.data)
    index = 3
    len = 5
    ts_subset = IS.get_time_series(storage, IS.get_uuid(ts); index=index, len=len)
    @test ts_subset[1] == ts_data[index]
    @test length(ts_subset) == len
end

function test_clear(storage::IS.TimeSeriesStorage)
    name = "component1"
    label = "val"
    component = IS.TestComponent(name, 5)
    ts = create_time_series_data()
    IS.add_time_series!(storage, IS.get_uuid(component), label, ts)

    ts_data = IS.get_time_series(storage, IS.get_uuid(ts))

    @test TimeSeries.timestamp(ts_data) == TimeSeries.timestamp(ts.data)
    @test TimeSeries.values(ts_data) == TimeSeries.values(ts.data)

    IS.clear_time_series!(storage)
    @test_throws ArgumentError IS.get_time_series(storage, IS.get_uuid(ts))
end

@testset "Test time series storage implementations" begin
    for in_memory in (true, false)
        test_add_remove(IS.make_time_series_storage(; in_memory = in_memory))
        test_get_subset(IS.make_time_series_storage(; in_memory = in_memory))
        test_clear(IS.make_time_series_storage(; in_memory = in_memory))
    end
end


text:
using gcSolver
using Plots
using OrdinaryDiffEq
using DiffEqDevTools

rxntfR = ones(gcSolver.Nparams) * 0.2

prob = gcSolver.runCkineSetup([500.0], rxntfR)

setups = [
    Dict(:alg => AutoTsit5(Rodas5())),
    Dict(:alg => AutoTsit5(Kvaerno5())),
    Dict(:alg => AutoVern9(Rodas5())),
    Dict(:alg => AutoVern9(KenCarp5())),
    Dict(:alg => AutoVern9(Kvaerno5())),
]

test_sol = solve(prob, AutoTsit5(Rodas5()), reltol = 1.0e-16, abstol = 1.0e-16)

tols = 1.0 ./ 10.0 .^ LinRange(4, 12, 30)

wp = WorkPrecisionSet(prob, tols, tols, setups; save_everystep = false, appxsol = test_sol, maxiters = Int(1e6))

plot(wp)


text:
"""
    install_dash(url, tag)

Clone python dash into `dash` folder and install (reinstall) it to current python enviroment
"""
function install_dash(url, tag)
    Conda.pip_interop(true)
    Conda.pip("uninstall -y", "dash")
    rm("dash", force = true, recursive = true)
    mkdir("dash")
    repo = clone(url, "dash")

    obj = try
        LibGit2.GitObject(repo, tag)
    catch
        LibGit2.GitObject(repo, "origin/$tag")
    end
    source_hash = LibGit2.string(LibGit2.GitHash(obj))
    LibGit2.checkout!(repo, source_hash)

    cd("dash") do
        Conda.pip("install", "./")
    end
    dashmodule = pyimport("dash")
    println(dashmodule)
    return dashmodule.__version__
end

dash_version() = VersionNumber(pyimport("dash").__version__)
dash_module_dir() = dirname(pyimport("dash").__file__)

function convert_deps_part(relative, external)
    result = Dict{Symbol, String}[]
    for i in eachindex(relative)
        d = OrderedDict{Symbol, String}()
        d[:relative_package_path] = relative[i]
        if !isnothing(external)
            d[:external_url] = external[i]
        end
        push!(result, d)
    end
    return result
end
function convert_deps(pyresources)
    relative_paths = pyresources["relative_package_path"]
    external_urls = get(pyresources, "external_url", nothing)
    result = OrderedDict{Symbol, Any}()
    result[:namespace] = pyresources["namespace"]
    if haskey(relative_paths, "prod")
        result[:prod] = convert_deps_part(
            relative_paths["prod"],
            isnothing(external_urls) ? nothing : external_urls["prod"]
            )
    end
    if haskey(relative_paths, "dev")
        result[:dev] = convert_deps_part(
            relative_paths["dev"],
            isnothing(external_urls) ? nothing : external_urls["dev"]
            )
    end
    return result
end

function _process_dist_part!(dict, resource, type)
    ns_symbol = Symbol(resource["namespace"])
    if !haskey(dict, ns_symbol)
        dict[ns_symbol] = OrderedDict(
            :namespace => resource["namespace"],
            :resources => Dict[]
        )
    end
    data = filter(v->v.first != "namespace", resource)
    data[:type] = type
    push!(
        dict[ns_symbol][:resources], data
        )
end
function convert_resources(js_dist, css_dist)
    result = OrderedDict{Symbol, Any}()
    _process_dist_part!.(Ref(result), js_dist, :js)
    _process_dist_part!.(Ref(result), css_dist, :css)
    return collect(values(result))
end

function deps_files(pyresources)
    relative_paths = pyresources["relative_package_path"]
    result = String[]
    if haskey(relative_paths, "prod")
        append!(result, relative_paths["prod"])
    end
    if haskey(relative_paths, "dev")
        append!(result, relative_paths["dev"])
    end
    return result
end

function resources_files(pyresources)
    result = String[]
    for res in pyresources
        if haskey(res, "relative_package_path")
            push!(result, res["relative_package_path"])
        end
        if haskey(res, "dev_package_path")
            push!(result, res["dev_package_path"])
        end
    end
    return result
end

function copy_files(resource_dir, files, dest_dir)
    for f in files
        full_path = joinpath(resource_dir, f)
        path_dir = dirname(f)
        mkpath(joinpath(dest_dir, path_dir))
        if isfile(full_path)
            cp(full_path, joinpath(dest_dir, f); force = true)
        else
            @warn "File $f don't exists in $resource_dir, skipped"
        end
    end
end

function renderer_resources(module_name)
    m = pyimport(module_name)

    meta =  OrderedDict(
        :version => m.__version__,
        :js_dist_dependencies =>  convert_deps(m._js_dist_dependencies[1]),
        :deps => convert_resources(m._js_dist, hasproperty(m, :_css_dist) ? m._css_dist : [])
    )
    files = vcat(deps_files(m._js_dist_dependencies[1]), resources_files(m._js_dist))
    return (dirname(m.__file__), meta, files)
end

function fill_components_deps!(deps::OrderedDict, files::Vector, module_name)
    m = pyimport(module_name)
    js_dist = hasproperty(m, :_js_dist) ? m._js_dist : []
    css_dist = hasproperty(m, :_css_dist) ? m._css_dist : []
    _process_dist_part!.(Ref(deps), js_dist, :js)
    _process_dist_part!.(Ref(deps), css_dist, :css)
    append!(files, vcat(resources_files(js_dist), resources_files(css_dist)))
end

function components_module_resources(module_name; name, prefix, metadata_file)
    m = pyimport(module_name)
    moduledir = dirname(m.__file__)
    metafile = joinpath(moduledir, metadata_file)
    if !isfile(metafile)
        error("meta file $(metadata_file) don't exists for module $(name)")
    end
    components = process_components_meta(metafile)
    meta = OrderedDict(
        :version => m.__version__,
        :name => name,
        :prefix => prefix,
        :components => components
    )

    return meta
end


text:
"""
    inout_map_reach(partition::AbstractVector{<:AbstractVector{Int}},
                    blocks::AbstractVector{Int},
                    n::Int)

Compute a mapping from input dimension to output dimension.
The map returns `0` for all input dimensions not present in the output.

### Input

- `partition` -- block partition; elements are start and end indices of a block
- `blocks`    -- list of all output block indices in the partition
- `n`         -- total number of input dimensions

### Output

A vector mapping from input dimension (index) to output dimension (entry).
"""
function inout_map_reach(partition::AbstractVector{<:AbstractVector{Int}},
                         blocks::AbstractVector{Int},
                         n::Int)
    inout = Vector{Int}(undef, n)
    blocks_idx = 1
    next_block_idx = blocks[blocks_idx]
    out_idx = 1
    for (i, bi) in enumerate(partition)
        if next_block_idx == i
            # block/variables present in output
            inout[bi] = out_idx:(out_idx + length(bi) - 1)
            out_idx += length(bi)
            if blocks_idx < length(blocks)
                blocks_idx += 1
                next_block_idx = blocks[blocks_idx]
            else
                next_block_idx = 0
            end
        else
            @assert next_block_idx > i || next_block_idx == 0
            # block/variables not present in output
            inout[bi] .= 0
        end
    end
    return inout
end


text:
"""
A `GridBasis` is a discrete basis that is associated with a grid.

The domain of the grid basis is the index set of the grid.
"""
struct GridBasis{T,G <: AbstractGrid} <: DiscreteDictionary{LinearIndex,T}
    grid    ::  G
end

const GridBasis1d{T,G <: AbstractGrid1d} = GridBasis{T,G}
const ProductGridBasis{T,G <: ProductGrid} = GridBasis{T,G}

coefficienttype(gb::GridBasis{T,G}) where {T,G} = T

# If not specified, we guess the codomaintype of the grid from the element
# type of its  elements.
GridBasis(grid::AbstractGrid) = GridBasis{subeltype(grid)}(grid)

GridBasis{T}(grid::AbstractGrid) where {T} = GridBasis{T,typeof(grid)}(grid)

# gridbasis(grid::AbstractGrid) = GridBasis(grid)
# gridbasis(grid::AbstractGrid, T) = GridBasis{T}(grid)
GridBasis(d::Dictionary, g::AbstractGrid = interpolation_grid(d)) = GridBasis{coefficienttype(d)}(g)

grid(b::GridBasis) = b.grid

@forward GridBasis.grid size, eachindex

dimension(b::GridBasis) = GridArrays.dimension(grid(b))

show(io::IO, mime::MIME"text/plain", d::GridBasis) = composite_show(io, mime, d)
Display.displaystencil(d::GridBasis{T}) where T = ["GridBasis{$(T)}(", grid(d), ")"]

tensorproduct(dicts::GridBasis...) = GridBasis{promote_type(map(coefficienttype,dicts)...)}(productgrid(map(grid, dicts)...))

support(b::GridBasis) = coverdomain(grid(b))

# Convenience function: add grid as extra parameter to hastransform
hastransform(s1::Dictionary, s2::GridBasis) =
	hasgrid_transform(s1, s2, grid(s2))
hastransform(s1::GridBasis, s2::Dictionary) =
	hasgrid_transform(s2, s1, grid(s1))
# and provide a default
hasgrid_transform(s1::Dictionary, s2, grid) = false

components(s::ProductGridBasis) = map(GridBasis{coefficienttype(s)}, components(grid(s)))
component(s::ProductGridBasis, i) = GridBasis{coefficienttype(s)}(component(grid(s), i))

apply_map(s::GridBasis, map) = GridBasis{coefficienttype(s)}(apply_map(grid(s), map))

sample(s::GridBasis, f) = sample(grid(s), f, coefficienttype(s))

native_index(d::ProductGridBasis, idx) = product_native_index(size(d), idx)
ordering(d::ProductGridBasis) = ProductIndexList{dimension(grid(d))}(size(d))

function gridextension(T, src::GridBasis, dest::GridBasis, src_grid::IndexSubGrid, dest_grid::AbstractGrid; options...)
    @assert supergrid(src_grid) == dest_grid
    IndexExtension{T}(src, dest, subindices(src_grid))
end

function gridrestriction(T, src::GridBasis, dest::GridBasis, src_grid::AbstractGrid, dest_grid::IndexSubGrid; options...)
    @assert supergrid(dest_grid) == src_grid
    IndexRestriction{T}(src, dest, subindices(dest_grid))
end


text:
using ThreeBodyProblem
using LinearAlgebra
using Test

μ₁ = 398600     # Earth gravitational parameter
μ₂ = 4902       # Moon gravitational paramter
d = 384400      # Avg distance between Earth and Moon
p = [μ₁; μ₂; d]
μ = μ₂/(μ₁ + μ₂)
@test computeR1R2(μ) == (-μ, 1-μ)
@test computeL1(μ) == [0.8369247061700067;0;0]
@test computeL1(p) == [321713.8570517506;0;0]
@test computeL2(μ) == [1.1556746768583537;0;0]
@test computeL2(p) == [444241.3457843512;0;0]
@test computeL3(μ) == [-1.005061834632099;0;0]
@test computeL3(p) == [-386345.76923257887;0;0]
@test computeL4(μ) == [0.48785136133154233;0.8660254037844386;0]
@test computeL4(p) == [187530.06329584488;332900.16521473817;0]
@test computeL5(μ) == [0.48785136133154233;-0.8660254037844386;0]
@test computeL5(p) == [187530.06329584488;-332900.16521473817;0]

@test computeLpts(μ,tol=1e-3) != computeLpts(μ)
@test computeLpts(p,tol=1e-3) != computeLpts(p)

rv = [1,0,0,0,1,0]
@test computeUeff(rv,p) < 0
@test computeUeff([-μ;0;0],μ) == -Inf

@test_broken computeC(rv,p) > 0
@test_broken computeC([-μ;0;0],μ) == Inf

Φ = ones(6,6)
@test_broken stability_index(ϕ) > 1

sys = earth_moon()
Lpts = computeLpts(sys)
rv = vcat(Lpts[1],zeros(3))
C₁ = computeC(rv, sys)

A = [-1. 0.  0.;
      0. 1.  0.;
      0. 0. -1.]
@test det(roty(π) - A) == 0
@test det(rotyd(180) - A) == 0

mytime = [2000, 11, 17.625065277778]
mytime2 = [2000,11,17,15,0,5.63999999]
mytime3 = [1000, 1, 1]
mytime4 = [1582, 11, 1]
mytime5 = [1582, 10, 1]
mytime6 = [1582, 10, 11]
mytime7 = [1582, 10, 5]
mytime8 = [1993, 4, 19]
@test date2mjd(mytime2) == date2mjd(mytime)
@test date2mjd(mytime3) == -314577.0
@test date2mjd(mytime4) == -100823.0
@test date2mjd(mytime5) == -101728.0
@test date2mjd(mytime6) == -100844.0
@test isnan(date2mjd(mytime7))
@test date2mjd(mytime8) == 49096.0

@test_throws ErrorException("ut1_date should be [Y,M,D] or [Y,M,D,h,m,s]") date2mjd(1)

@test wrapto360(314) == 314
@test wrapto360(372) == 12
@test wrapto360(721) == 1
@test wrapto360(-1) == 359
@test wrapto360(-361) == 359

@test wrapto180(181) == -179
@test wrapto180(179) == 179
@test wrapto180(-181) == 179

@test wrapto2pi(2π) == 2π
@test wrapto2pi(3π) - π == 0.
@test wrapto2pi(-π) - π == 0.

@test wraptopi(π) == π
@test wraptopi(2π) == 0.
@test wraptopi(-2π) == 0.

@test rotlatlon(π/2,-π/2,ang_unit=:rad) == I
@test_throws ErrorException("ang_unit should be :deg or :rad") rotlatlon(0,0,ang_unit=:foo)

@test mjd2gmst(51544.5,ang_unit=:deg) == 280.4606
@test_throws ErrorException("ang_unit should be :deg or :rad") mjd2gmst(0,ang_unit=:foo)

xyz_sphere,N = deserno_hemisphere(100,[1,0,0])
@test abs(minimum([norm(xyz_sphere[:,i]) for i = 1:N]) - 1) < 1e-10
@test abs(maximum([norm(xyz_sphere[:,i]) for i = 1:N]) - 1) < 1e-10


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

# ToDo: bat_default(::typeof(bat_integrate), ::Val{:algorithm}, ::AnySampleable) = BridgeSampling()
bat_default(::typeof(bat_integrate), ::Val{:algorithm}, ::SampledDensity) = BridgeSampling()


text:
import PyCall
PyCall.pyimport_conda("sympy.physics.wigner",       "sympy")
PyCall.pyimport_conda("sympy.physics.optics",       "sympy")
PyCall.pyimport_conda("sympy.physics.quantum.spin", "sympy")

## to use
## import_from(sympy.physics.wigner)
## wigner3j(sympy.S(6),0,4,0,2,0)
## Wigner3j(sympy.S(6),0,4,0,2,0).doit()


text:
"""
    asymptotic_monopole_wavenumbers(ω, medium::Acoustic, species; num_wavenumbers = 2)

Calculates the asymptotic effective wavenumbers for monopole scatterers by assuming a large number of wavenumbers. The 2D results are deduced in [Section 5](https://arxiv.org/pdf/1905.06996.pdf).
"""
function asymptotic_monopole_wavenumbers(ω::T, medium::Acoustic{T,2}, species::Species{T,2};
        num_wavenumbers = 2
    ) where T

    P = Int(round(num_wavenumbers / 2) + 1)
    k = ω / medium.c

    b = mean(s1.exclusion_distance * outer_radius(s1) + s2.exclusion_distance * outer_radius(s2) for s1 in species, s2 in species)
    num_density = mean(number_density(s) for s in species)

    T0 = mean(t_matrix(s.particle, medium, ω, 0) for s in species)[1]

    c = sqrt(2pi) * num_density * b^2 * T0 * hankelh1(0,k * b) * exp(-im * pi / 4.0)

    rc = abs(c)
    θc = angle(c)

    kps = map(-ceil(θc / (2pi)):P) do p
        σ = θc + 2pi * p
        kp = (σ + im * log(abs(σ)^(3/2) / rc)) / b
    end

    kns = map(ceil(θc / (2pi) - 3 / 4):P) do p
        σ = θc - 3pi / 2 - 2pi * p
        kn = (σ + im * log(abs(σ)^(3/2) / rc)) / b
    end

   return [kns; kps]
end

function asymptotic_monopole_wavenumbers(ω::T, medium::Acoustic{T,3}, species::Species{T,3};
        num_wavenumbers = 2
    ) where T

   @warn "asymptotic_monopole_wavenumbers has not yet been implemented for 3D acoustics, will use 2D acoustics instead."

   medium2D = Acoustic(2; ρ=medium.ρ, c=medium.c)

   species2D = [
       Specie(Acoustic(2; ρ = s.particle.medium.ρ, c = s.particle.medium.c), Circle(outer_radius(s));
           volume_fraction = volume_fraction(s))
   for s in species]

   return asymptotic_monopole_wavenumbers(ω, medium2D, species2D; num_wavenumbers = num_wavenumbers)
end

max_stars_repo_name: arturgower/EffectiveWaves


text:
module LCAlgebra

# +       NUM     LC      QEQ
# NUM     NUM     LC      QEQ
# LC      LC      LC      QEQ
# QEQ     QEQ     QEQ     ERR

# *       NUM     LC      QEQ
# NUM     NUM     LC      QEQ
# LC      LC      QEQ     ERR
# QEQ     QEQ     ERR     ERR

using FF
using SparseArrays, BitIntegers
import SparseArrays: sparsevec, nonzeroinds
import Base: +, -, *, /

export Signal, SignalState, LinearCombination, QEQ, islc, ids, iszero, substitute

struct LinearCombination
  v::SparseVector{ff,Int64}
end

struct Signal
  id::Int64
end

struct QEQ
  a::LinearCombination
  b::LinearCombination
  c::LinearCombination
end

const maxConstraintSize=1<<30

LCA = Union{Integer, Signal, LinearCombination}
LCB = Union{Signal, LinearCombination}
LCP = Union{LinearCombination, QEQ}
LCU = Union{LCA, LCP}

Signal(x::Signal) = x
LinearCombination(x::Signal) = LinearCombination(dropzeros!(sparsevec([x.id], [ff(1)], maxConstraintSize)))
LinearCombination(x::t) where t<:Integer = LinearCombination(dropzeros!(sparsevec([1], [ff(x)], maxConstraintSize)))
LinearCombination(x::LinearCombination) = x
QEQ(x::LinearCombination) = QEQ(LinearCombination(SparseVector(maxConstraintSize, Int64[], ff[])), LinearCombination(SparseVector(maxConstraintSize, Int64[], ff[])), x)
QEQ(x::T) where T<: LCA = QEQ(LinearCombination(x))
QEQ(x::QEQ) = x

(-)(x::Signal) = -LinearCombination(x)
(-)(x::LinearCombination) = LinearCombination(-x.v)
(-)(x::QEQ) = QEQ(-x.a, x.b, -x.c)

(+)(x::LinearCombination, y::LinearCombination) = LinearCombination(x.v+y.v)
(+)(x::QEQ, y::QEQ) = error("additive operations for (::QEQ, ::QEQ) are not determined")

(+)(x::T, y::QEQ) where T<:LCA = QEQ(y.a, y.b, y.c+LinearCombination(x))
(+)(x::QEQ, y::T) where T<:LCA = QEQ(x.a, x.b, x.c+LinearCombination(y))
(+)(x::A, y::B) where {A<:LCA, B<:LCA}  = LinearCombination(x)+LinearCombination(y)
(-)(x::A, y::B) where {A<:LCU, B<:LCU} = x + (-y)

(*)(x::P, y::T) where {P<:LCB, T<:LCB}= QEQ(LinearCombination(x), LinearCombination(y), LinearCombination(SparseVector(maxConstraintSize, Int64[], ff[])))
(*)(x::QEQ, y::QEQ) = error("multiplicative operations for (::QEQ, ::QEQ) are not determined")

(*)(x::T, y::QEQ) where T<: LCB = error("multiplicative operations for (::LinearCombination, ::QEQ) are not determined")
(*)(x::QEQ, y::T) where T<: LCB = error("multiplicative operations for (::QEQ, ::LinearCombination) are not determined")
(*)(x::QEQ, y::T) where T<:Integer = QEQ(x.a*ff(y), x.b, x.c*ff(y))
(/)(x::QEQ, y::T) where T<:Integer = QEQ(x.a/ff(y), x.b, x.c/ff(y))
(*)(x::T, y::QEQ) where T<:Integer = QEQ(y.a*ff(y), y.b, y.c*ff(y))

(*)(x::P, y::T) where {P<:LCB, T<:Integer} = LinearCombination(LinearCombination(x).v*ff(y))
(/)(x::P, y::T) where {P<:LCB, T<:Integer} = LinearCombination(LinearCombination(x).v/ff(y))
(*)(x::P, y::T) where {P<:Integer, T<:LCB} = LinearCombination(LinearCombination(y).v*ff(x))

cmp(x::LinearCombination, y::LinearCombination) = sign(cmp(nnz(x.v), nnz(y.v))*4 + cmp(nonzeros(x.v), nonzeros(y.v))*2 + cmp(nonzeroinds(x.v), nonzeroinds(y.v)))

# function isZero(a) {
#   if (a.type == "NUMBER") {
#       return a.value.isZero();
#   } else if (a.type == "LINEARCOMBINATION") {
#       for (let k in a.values) {
#           if (!a.values[k].isZero()) return false;
#       }
#       return true;
#   } else if (a.type == "QEQ") {
#       return (isZero(a.a) || isZero(a.b)) && isZero(a.c);
#   } else if (a.type == "ERROR") {
#       return false;
#   } else {
#       return false;
#   }
# }

iszero(::Signal) = false
iszero(x::LinearCombination) = nnz(x.v)==0
iszero(x::QEQ) = iszero(x.c) && (iszero(x.a) || iszero(x.b))
islc(x::QEQ) = iszero(x.a) || iszero(x.b)
ids(x::LinearCombination) = nonzeroinds(x.v)
ids(x::QEQ) = unique(vcat(nonzeroinds(x.a.v), nonzeroinds(x.b.v), nonzeroinds(x.c.v)))

function substitute(x::LinearCombination, s::Signal, zero::LinearCombination)
  t = zero.v[s.id]
  if t == 0
    return x
  end
  return  x - x.v[s.id]/t*zero
end

substitute(x::QEQ, s::Signal, zero::LinearCombination)=canonize!(QEQ(substitute(x.a, s, zero), substitute(x.b, s, zero), substitute(x.c, s, zero)))

evaluate(x::Signal, values::T) where T<: AbstractVector = values[x.id]
evaluate(x::LinearCombination, values::T) where T<: AbstractVector = sum(x .* values)
evaluate(x::QEQ, values::T) where T<: AbstractVector = evaluate(x.a, values) * evaluate(x.b, values) + evaluate(x.c, values)

function canonize!(x::QEQ)
  a = x.a.v[1]
  b = x.b.v[1]
  if iszero(x.a) || iszero(x.b)
    copyto!(x.a.v, sparsevec(ff[], maxConstraintSize))
    copyto!(x.b.v, sparsevec(ff[], maxConstraintSize))
  elseif a!=0 || b!=0
    x.a.v[1] = 0
    x.b.v[1] = 0
    copyto!(x.c.v, (x.c + x.a * b + x.b * a + a*b).v)
  end
  return x
end

end # module