text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

"""
Provide serialization of Julia objects via the functions
* [`serialize`](@ref)
* [`deserialize`](@ref)
"""
module Serialization

import Base: GMP, Bottom, unsafe_convert, uncompressed_ast
import Core: svec, SimpleVector
using Base: unaliascopy, unwrap_unionall, require_one_based_indexing, ntupleany
using Core.IR

export serialize, deserialize, AbstractSerializer, Serializer

abstract type AbstractSerializer end

mutable struct Serializer{I<:IO} <: AbstractSerializer
    io::I
    counter::Int
    table::IdDict{Any,Any}
    pending_refs::Vector{Int}
    known_object_data::Dict{UInt64,Any}
    version::Int
    Serializer{I}(io::I) where I<:IO = new(io, 0, IdDict(), Int[], Dict{UInt64,Any}(), ser_version)
end

Serializer(io::IO) = Serializer{typeof(io)}(io)

## serializing values ##

const n_int_literals = 33
const n_reserved_slots = 24
const n_reserved_tags = 8

const TAGS = Any[
    Symbol, Int8, UInt8, Int16, UInt16, Int32, UInt32, Int64, UInt64, Int128, UInt128,
    Float16, Float32, Float64, Char, DataType, Union, UnionAll, Core.TypeName, Tuple,
    Array, Expr, LineNumberNode, :__LabelNode__, GotoNode, QuoteNode, CodeInfo, TypeVar,
    Core.Box, Core.MethodInstance, Module, Task, String, SimpleVector, Method,
    GlobalRef, SlotNumber, TypedSlot, NewvarNode, SSAValue,

    # dummy entries for tags that don't correspond directly to types
    Symbol, # UNDEFREF_TAG
    Symbol, # BACKREF_TAG
    Symbol, # LONGBACKREF_TAG
    Symbol, # SHORTBACKREF_TAG
    Symbol, # LONGTUPLE_TAG
    Symbol, # LONGSYMBOL_TAG
    Symbol, # LONGEXPR_TAG
    Symbol, # LONGSTRING_TAG
    Symbol, # SHORTINT64_TAG
    Symbol, # FULL_DATATYPE_TAG
    Symbol, # WRAPPER_DATATYPE_TAG
    Symbol, # OBJECT_TAG
    Symbol, # REF_OBJECT_TAG
    Symbol, # FULL_GLOBALREF_TAG
    Symbol, # HEADER_TAG
    Symbol, # IDDICT_TAG
    Symbol, # SHARED_REF_TAG
    ReturnNode, GotoIfNot,
    fill(Symbol, n_reserved_tags)...,

    (), Bool, Any, Bottom, Core.TypeofBottom, Type, svec(), Tuple{}, false, true, nothing,
    :Any, :Array, :TypeVar, :Box, :Tuple, :Ptr, :return, :call, Symbol("::"), :Function,
    :(=), :(==), :(===), :gotoifnot, :A, :B, :C, :M, :N, :T, :S, :X, :Y, :a, :b, :c, :d, :e, :f,
    :g, :h, :i, :j, :k, :l, :m, :n, :o, :p, :q, :r, :s, :t, :u, :v, :w, :x, :y, :z, :add_int,
    :sub_int, :mul_int, :add_float, :sub_float, :new, :mul_float, :bitcast, :start, :done, :next,
    :indexed_iterate, :getfield, :meta, :eq_int, :slt_int, :sle_int, :ne_int, :push_loc, :pop_loc,
    :pop, :arrayset, :arrayref, :apply_type, :inbounds, :getindex, :setindex!, :Core, :!, :+,
    :Base, :static_parameter, :convert, :colon, Symbol("#self#"), Symbol("#temp#"), :tuple, Symbol(""),

    fill(:_reserved_, n_reserved_slots)...,

    (Int32(0):Int32(n_int_literals-1))...,
    (Int64(0):Int64(n_int_literals-1))...
]

@assert length(TAGS) == 255

const ser_version = 15 # do not make changes without bumping the version #!

format_version(::AbstractSerializer) = ser_version
format_version(s::Serializer) = s.version

const NTAGS = length(TAGS)

function sertag(@nospecialize(v))
    # NOTE: we use jl_value_ptr directly since we know at least one of the arguments
    # in the comparison below is a singleton.
    ptr = ccall(:jl_value_ptr, Ptr{Cvoid}, (Any,), v)
    ptags = convert(Ptr{Ptr{Cvoid}}, pointer(TAGS))
    # note: constant ints & reserved slots never returned here
    @inbounds for i in 1:(NTAGS-(n_reserved_slots+2*n_int_literals))
        ptr == unsafe_load(ptags,i) && return i%Int32
    end
    return Int32(-1)
end
desertag(i::Int32) = @inbounds(TAGS[i])

# tags >= this just represent themselves, their whole representation is 1 byte
const VALUE_TAGS = sertag(())
const ZERO32_TAG = Int32(NTAGS-(2*n_int_literals-1))
const ZERO64_TAG = Int64(NTAGS-(n_int_literals-1))
const TRUE_TAG = sertag(true)
const FALSE_TAG = sertag(false)
const EMPTYTUPLE_TAG = sertag(())
const TUPLE_TAG = sertag(Tuple)
const SIMPLEVECTOR_TAG = sertag(SimpleVector)
const SYMBOL_TAG = sertag(Symbol)
const INT8_TAG = sertag(Int8)
const ARRAY_TAG = findfirst(==(Array), TAGS)%Int32
const EXPR_TAG = sertag(Expr)
const MODULE_TAG = sertag(Module)
const METHODINSTANCE_TAG = sertag(Core.MethodInstance)
const METHOD_TAG = sertag(Method)
const TASK_TAG = sertag(Task)
const DATATYPE_TAG = sertag(DataType)
const TYPENAME_TAG = sertag(Core.TypeName)
const INT32_TAG = sertag(Int32)
const INT64_TAG = sertag(Int64)
const GLOBALREF_TAG = sertag(GlobalRef)
const BOTTOM_TAG = sertag(Bottom)
const UNIONALL_TAG = sertag(UnionAll)
const STRING_TAG = sertag(String)
const o0 = sertag(SSAValue)
const UNDEFREF_TAG         = Int32(o0+1)
const BACKREF_TAG          = Int32(o0+2)
const LONGBACKREF_TAG      = Int32(o0+3)
const SHORTBACKREF_TAG     = Int32(o0+4)
const LONGTUPLE_TAG        = Int32(o0+5)
const LONGSYMBOL_TAG       = Int32(o0+6)
const LONGEXPR_TAG         = Int32(o0+7)
const LONGSTRING_TAG       = Int32(o0+8)
const SHORTINT64_TAG       = Int32(o0+9)
const FULL_DATATYPE_TAG    = Int32(o0+10)
const WRAPPER_DATATYPE_TAG = Int32(o0+11)
const OBJECT_TAG           = Int32(o0+12)
const REF_OBJECT_TAG       = Int32(o0+13)
const FULL_GLOBALREF_TAG   = Int32(o0+14)
const HEADER_TAG           = Int32(o0+15)
const IDDICT_TAG           = Int32(o0+16)
const SHARED_REF_TAG       = Int32(o0+17)

writetag(s::IO, tag) = (write(s, UInt8(tag)); nothing)

function write_as_tag(s::IO, tag)
    tag < VALUE_TAGS && write(s, UInt8(0))
    write(s, UInt8(tag))
    nothing
end

# cycle handling
function serialize_cycle(s::AbstractSerializer, @nospecialize(x))
    offs = get(s.table, x, -1)::Int
    if offs != -1
        if offs <= typemax(UInt16)
            writetag(s.io, SHORTBACKREF_TAG)
            write(s.io, UInt16(offs))
        elseif offs <= typemax(Int32)
            writetag(s.io, BACKREF_TAG)
            write(s.io, Int32(offs))
        else
            writetag(s.io, LONGBACKREF_TAG)
            write(s.io, Int64(offs))
        end
        return true
    end
    s.table[x] = s.counter
    s.counter += 1
    return false
end

function serialize_cycle_header(s::AbstractSerializer, @nospecialize(x))
    serialize_cycle(s, x) && return true
    serialize_type(s, typeof(x), true)
    return false
end

function reset_state(s::AbstractSerializer)
    s.counter = 0
    empty!(s.table)
    empty!(s.pending_refs)
    s
end

serialize(s::AbstractSerializer, x::Bool) = x ? writetag(s.io, TRUE_TAG) :
                                                writetag(s.io, FALSE_TAG)

serialize(s::AbstractSerializer, p::Ptr) = serialize_any(s, oftype(p, C_NULL))

serialize(s::AbstractSerializer, ::Tuple{}) = writetag(s.io, EMPTYTUPLE_TAG)

function serialize(s::AbstractSerializer, t::Tuple)
    l = length(t)
    if l <= 255
        writetag(s.io, TUPLE_TAG)
        write(s.io, UInt8(l))
    else
        writetag(s.io, LONGTUPLE_TAG)
        write(s.io, Int32(l))
    end
    for x in t
        serialize(s, x)
    end
end

function serialize(s::AbstractSerializer, v::SimpleVector)
    writetag(s.io, SIMPLEVECTOR_TAG)
    write(s.io, Int32(length(v)))
    for x in v
        serialize(s, x)
    end
end

function serialize(s::AbstractSerializer, x::Symbol)
    tag = sertag(x)
    if tag > 0
        return write_as_tag(s.io, tag)
    end
    pname = unsafe_convert(Ptr{UInt8}, x)
    len = Int(ccall(:strlen, Csize_t, (Cstring,), pname))
    if len > 7
        serialize_cycle(s, x) && return
    end
    if len <= 255
        writetag(s.io, SYMBOL_TAG)
        write(s.io, UInt8(len))
    else
        writetag(s.io, LONGSYMBOL_TAG)
        write(s.io, Int32(len))
    end
    unsafe_write(s.io, pname, len)
    nothing
end

function serialize_array_data(s::IO, a)
    require_one_based_indexing(a)
    isempty(a) && return 0
    if eltype(a) === Bool
        last = a[1]::Bool
        count = 1
        for i = 2:length(a)
            if a[i]::Bool != last || count == 127
                write(s, UInt8((UInt8(last) << 7) | count))
                last = a[i]::Bool
                count = 1
            else
                count += 1
            end
        end
        write(s, UInt8((UInt8(last) << 7) | count))
    else
        write(s, a)
    end
end

function serialize(s::AbstractSerializer, a::Array)
    serialize_cycle(s, a) && return
    elty = eltype(a)
    writetag(s.io, ARRAY_TAG)
    if elty !== UInt8
        serialize(s, elty)
    end
    if ndims(a) != 1
        serialize(s, size(a))
    else
        serialize(s, length(a))
    end
    if isbitstype(elty)
        serialize_array_data(s.io, a)
    else
        sizehint!(s.table, div(length(a),4))  # prepare for lots of pointers
        @inbounds for i in eachindex(a)
            if isassigned(a, i)
                serialize(s, a[i])
            else
                writetag(s.io, UNDEFREF_TAG)
            end
        end
    end
end

function serialize(s::AbstractSerializer, a::SubArray{T,N,A}) where {T,N,A<:Array}
    # SubArray's copy only selects the relevant data (and reduces the size) but does not
    # preserve the type of the argument. This internal function does both:
    b = unaliascopy(a)
    serialize_any(s, b)
end

function serialize(s::AbstractSerializer, ss::String)
    len = sizeof(ss)
    if len > 7
        serialize_cycle(s, ss) && return
        writetag(s.io, SHARED_REF_TAG)
    end
    if len <= 255
        writetag(s.io, STRING_TAG)
        write(s.io, UInt8(len))
    else
        writetag(s.io, LONGSTRING_TAG)
        write(s.io, Int64(len))
    end
    write(s.io, ss)
    nothing
end

function serialize(s::AbstractSerializer, ss::SubString{String})
    # avoid saving a copy of the parent string, keeping the type of ss
    serialize_any(s, SubString(String(ss)))
end

# Don't serialize the pointers
function serialize(s::AbstractSerializer, r::Regex)
    serialize_type(s, typeof(r))
    serialize(s, r.pattern)
    serialize(s, r.compile_options)
    serialize(s, r.match_options)
end

function serialize(s::AbstractSerializer, n::BigInt)
    serialize_type(s, BigInt)
    serialize(s, string(n, base = 62))
end

function serialize(s::AbstractSerializer, ex::Expr)
    serialize_cycle(s, ex) && return
    l = length(ex.args)
    if l <= 255
        writetag(s.io, EXPR_TAG)
        write(s.io, UInt8(l))
    else
        writetag(s.io, LONGEXPR_TAG)
        write(s.io, Int32(l))
    end
    serialize(s, ex.head)
    for a in ex.args
        serialize(s, a)
    end
end

function serialize_dict_data(s::AbstractSerializer, d::AbstractDict)
    write(s.io, Int32(length(d)))
    for (k,v) in d
        serialize(s, k)
        serialize(s, v)
    end
end

function serialize(s::AbstractSerializer, d::Dict)
    serialize_cycle_header(s, d) && return
    serialize_dict_data(s, d)
end

function serialize(s::AbstractSerializer, d::IdDict)
    serialize_cycle(s, d) && return
    writetag(s.io, IDDICT_TAG)
    serialize_type_data(s, typeof(d))
    serialize_dict_data(s, d)
end

function serialize_mod_names(s::AbstractSerializer, m::Module)
    p = parentmodule(m)
    if p === m || m === Base
        key = Base.root_module_key(m)
        serialize(s, key.uuid === nothing ? nothing : key.uuid.value)
        serialize(s, Symbol(key.name))
    else
        serialize_mod_names(s, p)
        serialize(s, nameof(m))
    end
end

function serialize(s::AbstractSerializer, m::Module)
    writetag(s.io, MODULE_TAG)
    serialize_mod_names(s, m)
    writetag(s.io, EMPTYTUPLE_TAG)
end

# TODO: make this bidirectional, so objects can be sent back via the same key
const object_numbers = WeakKeyDict()
const obj_number_salt = Ref{UInt64}(0)
function object_number(s::AbstractSerializer, @nospecialize(l))
    global obj_number_salt, object_numbers
    if haskey(object_numbers, l)
        return object_numbers[l]
    end
    ln = obj_number_salt[]
    object_numbers[l] = ln
    obj_number_salt[] += 1
    return ln::UInt64
end

lookup_object_number(s::AbstractSerializer, n::UInt64) = nothing

remember_object(s::AbstractSerializer, @nospecialize(o), n::UInt64) = nothing

function lookup_object_number(s::Serializer, n::UInt64)
    return get(s.known_object_data, n, nothing)
end

function remember_object(s::Serializer, @nospecialize(o), n::UInt64)
    s.known_object_data[n] = o
    return nothing
end

function serialize(s::AbstractSerializer, meth::Method)
    serialize_cycle(s, meth) && return
    writetag(s.io, METHOD_TAG)
    write(s.io, object_number(s, meth))
    serialize(s, meth.module)
    serialize(s, meth.name)
    serialize(s, meth.file)
    serialize(s, meth.line)
    serialize(s, meth.sig)
    serialize(s, meth.slot_syms)
    serialize(s, meth.nargs)
    serialize(s, meth.isva)
    serialize(s, meth.is_for_opaque_closure)
    serialize(s, meth.aggressive_constprop)
    if isdefined(meth, :source)
        serialize(s, Base._uncompressed_ast(meth, meth.source))
    else
        serialize(s, nothing)
    end
    if isdefined(meth, :generator)
        serialize(s, meth.generator)
    else
        serialize(s, nothing)
    end
    if isdefined(meth, :recursion_relation)
        serialize(s, method.recursion_relation)
    else
        serialize(s, nothing)
    end
    nothing
end

function serialize(s::AbstractSerializer, linfo::Core.MethodInstance)
    serialize_cycle(s, linfo) && return
    writetag(s.io, METHODINSTANCE_TAG)
    if isdefined(linfo, :uninferred)
        serialize(s, linfo.uninferred)
    else
        writetag(s.io, UNDEFREF_TAG)
    end
    serialize(s, nothing)  # for backwards compat
    serialize(s, linfo.sparam_vals)
    serialize(s, Any)  # for backwards compat
    serialize(s, linfo.specTypes)
    serialize(s, linfo.def)
    nothing
end

function serialize(s::AbstractSerializer, t::Task)
    serialize_cycle(s, t) && return
    if istaskstarted(t) && !istaskdone(t)
        error("cannot serialize a running Task")
    end
    state = [t.code, t.storage, t.state, t.result, t._isexception]
    writetag(s.io, TASK_TAG)
    for fld in state
        serialize(s, fld)
    end
end

function serialize(s::AbstractSerializer, g::GlobalRef)
    if (g.mod === __deserialized_types__ ) ||
        (g.mod === Main && isdefined(g.mod, g.name) && isconst(g.mod, g.name))

        v = getfield(g.mod, g.name)
        unw = unwrap_unionall(v)
        if isa(unw,DataType) && v === unw.name.wrapper && should_send_whole_type(s, unw)
            # handle references to types in Main by sending the whole type.
            # needed to be able to send nested functions (#15451).
            writetag(s.io, FULL_GLOBALREF_TAG)
            serialize(s, v)
            return
        end
    end
    writetag(s.io, GLOBALREF_TAG)
    serialize(s, g.mod)
    serialize(s, g.name)
end

function serialize(s::AbstractSerializer, t::Core.TypeName)
    serialize_cycle(s, t) && return
    writetag(s.io, TYPENAME_TAG)
    write(s.io, object_number(s, t))
    serialize_typename(s, t)
end

function serialize_typename(s::AbstractSerializer, t::Core.TypeName)
    serialize(s, t.name)
    serialize(s, t.names)
    primary = unwrap_unionall(t.wrapper)
    serialize(s, primary.super)
    serialize(s, primary.parameters)
    serialize(s, primary.types)
    serialize(s, isdefined(primary, :instance))
    serialize(s, primary.abstract)
    serialize(s, primary.mutable)
    serialize(s, primary.ninitialized)
    if isdefined(t, :mt) && t.mt !== Symbol.name.mt
        serialize(s, t.mt.name)
        serialize(s, collect(Base.MethodList(t.mt)))
        serialize(s, t.mt.max_args)
        if isdefined(t.mt, :kwsorter)
            serialize(s, t.mt.kwsorter)
        else
            writetag(s.io, UNDEFREF_TAG)
        end
    else
        writetag(s.io, UNDEFREF_TAG)
    end
    nothing
end

# decide whether to send all data for a type (instead of just its name)
function should_send_whole_type(s, t::DataType)
    tn = t.name
    if isdefined(tn, :mt)
        # TODO improve somehow
        # send whole type for anonymous functions in Main
        name = tn.mt.name
        mod = tn.module
        isanonfunction = mod === Main && # only Main
            t.super === Function && # only Functions
            unsafe_load(unsafe_convert(Ptr{UInt8}, tn.name)) == UInt8('#') && # hidden type
            (!isdefined(mod, name) || t != typeof(getfield(mod, name))) # XXX: 95% accurate test for this being an inner function
            # TODO: more accurate test? (tn.name !== "#" name)
        #TODO: iskw = startswith(tn.name, "#kw#") && ???
        #TODO: iskw && return send-as-kwftype
        return mod === __deserialized_types__ || isanonfunction
    end
    return false
end

function serialize_type_data(s, @nospecialize(t::DataType))
    whole = should_send_whole_type(s, t)
    iswrapper = (t === unwrap_unionall(t.name.wrapper))
    if whole && iswrapper
        writetag(s.io, WRAPPER_DATATYPE_TAG)
        serialize(s, t.name)
        return
    end
    serialize_cycle(s, t) && return
    if whole
        writetag(s.io, FULL_DATATYPE_TAG)
        serialize(s, t.name)
    else
        writetag(s.io, DATATYPE_TAG)
        tname = t.name.name
        serialize(s, tname)
        mod = t.name.module
        serialize(s, mod)
    end
    if !isempty(t.parameters)
        if iswrapper
            write(s.io, Int32(0))
        else
            write(s.io, Int32(length(t.parameters)))
            for p in t.parameters
                serialize(s, p)
            end
        end
    end
    nothing
end

function serialize(s::AbstractSerializer, t::DataType)
    tag = sertag(t)
    tag > 0 && return write_as_tag(s.io, tag)
    if t === Tuple
        # `sertag` is not able to find types === to `Tuple` because they
        # will not have been hash-consed. Plus `serialize_type_data` does not
        # handle this case correctly, since Tuple{} != Tuple. `Tuple` is the
        # only type with this property. issue #15849
        return write_as_tag(s.io, TUPLE_TAG)
    end
    serialize_type_data(s, t)
end

function serialize_type(s::AbstractSerializer, @nospecialize(t::DataType), ref::Bool = false)
    tag = sertag(t)
    tag > 0 && return writetag(s.io, tag)
    writetag(s.io, ref ? REF_OBJECT_TAG : OBJECT_TAG)
    serialize_type_data(s, t)
end

function serialize(s::AbstractSerializer, n::Int32)
    if 0 <= n <= (n_int_literals-1)
        write(s.io, UInt8(ZERO32_TAG+n))
    else
        writetag(s.io, INT32_TAG)
        write(s.io, n)
    end
    nothing
end

function serialize(s::AbstractSerializer, n::Int64)
    if 0 <= n <= (n_int_literals-1)
        write(s.io, UInt8(ZERO64_TAG+n))
    elseif typemin(Int32) <= n <= typemax(Int32)
        writetag(s.io, SHORTINT64_TAG)
        write(s.io, Int32(n))
    else
        writetag(s.io, INT64_TAG)
        write(s.io, n)
    end
    nothing
end

for i in 0:13
    tag = Int32(INT8_TAG + i)
    ty = TAGS[tag]
    (ty === Int32 || ty === Int64) && continue
    @eval serialize(s::AbstractSerializer, n::$ty) = (writetag(s.io, $tag); write(s.io, n); nothing)
end

serialize(s::AbstractSerializer, ::Type{Bottom}) = write_as_tag(s.io, BOTTOM_TAG)

function serialize(s::AbstractSerializer, u::UnionAll)
    writetag(s.io, UNIONALL_TAG)
    n = 0; t = u
    while isa(t, UnionAll)
        t = t.body
        n += 1
    end
    if isa(t, DataType) && t === unwrap_unionall(t.name.wrapper)
        write(s.io, UInt8(1))
        write(s.io, Int16(n))
        serialize(s, t)
    else
        write(s.io, UInt8(0))
        serialize(s, u.var)
        serialize(s, u.body)
    end
end

serialize(s::AbstractSerializer, @nospecialize(x)) = serialize_any(s, x)

function serialize_any(s::AbstractSerializer, @nospecialize(x))
    tag = sertag(x)
    if tag > 0
        return write_as_tag(s.io, tag)
    end
    t = typeof(x)::DataType
    nf = nfields(x)
    if nf == 0 && t.size > 0
        serialize_type(s, t)
        write(s.io, x)
    else
        if t.mutable
            serialize_cycle(s, x) && return
            serialize_type(s, t, true)
        else
            serialize_type(s, t, false)
        end
        for i in 1:nf
            if isdefined(x, i)
                serialize(s, getfield(x, i))
            else
                writetag(s.io, UNDEFREF_TAG)
            end
        end
    end
    nothing
end

"""
    Serialization.writeheader(s::AbstractSerializer)

Write an identifying header to the specified serializer. The header consists of
8 bytes as follows:

| Offset | Description                                     |
|:-------|:------------------------------------------------|
|   0    | tag byte (0x37)                                 |
|   1-2  | signature bytes "JL"                            |
|   3    | protocol version                                |
|   4    | bits 0-1: endianness: 0 = little, 1 = big       |
|   4    | bits 2-3: platform: 0 = 32-bit, 1 = 64-bit      |
|   5-7  | reserved                                        |
"""
function writeheader(s::AbstractSerializer)
    io = s.io
    writetag(io, HEADER_TAG)
    write(io, "JL")  # magic bytes
    write(io, UInt8(ser_version))
    endianness = (ENDIAN_BOM == 0x04030201 ? 0 :
                  ENDIAN_BOM == 0x01020304 ? 1 :
                  error("unsupported endianness in serializer"))
    machine = (sizeof(Int) == 4 ? 0 :
               sizeof(Int) == 8 ? 1 :
               error("unsupported word size in serializer"))
    write(io, UInt8(endianness) | (UInt8(machine) << 2))
    write(io, [0x00,0x00,0x00]) # 3 reserved bytes
    nothing
end

function readheader(s::AbstractSerializer)
    # Tag already read
    io = s.io
    m1 = read(io, UInt8)
    m2 = read(io, UInt8)
    if m1 != UInt8('J') || m2 != UInt8('L')
        error("Unsupported serialization format (got header magic bytes $m1 $m2)")
    end
    version    = read(io, UInt8)
    flags      = read(io, UInt8)
    reserved1  = read(io, UInt8)
    reserved2  = read(io, UInt8)
    reserved3  = read(io, UInt8)
    endianflag = flags & 0x3
    wordflag   = (flags >> 2) & 0x3
    wordsize = wordflag == 0 ? 4 :
               wordflag == 1 ? 8 :
               error("Unknown word size flag in header")
    endian_bom = endianflag == 0 ? 0x04030201 :
                 endianflag == 1 ? 0x01020304 :
                 error("Unknown endianness flag in header")
    # Check protocol compatibility.
    endian_bom == ENDIAN_BOM  || error("Serialized byte order mismatch ($(repr(endian_bom)))")
    # We don't check wordsize == sizeof(Int) here, as Int is encoded concretely
    # as Int32 or Int64, which should be enough to correctly deserialize a range
    # of data structures between Julia versions.
    if version > ser_version
        error("""Cannot read stream serialized with a newer version of Julia.
                 Got data version $version > current version $ser_version""")
    end
    s.version = version
    return
end

"""
    serialize(stream::IO, value)

Write an arbitrary value to a stream in an opaque format, such that it can be read back by
[`deserialize`](@ref). The read-back value will be as identical as possible to the original.
In general, this process will not work if the reading and writing are done by different
versions of Julia, or an instance of Julia with a different system image. `Ptr` values are
serialized as all-zero bit patterns (`NULL`).

An 8-byte identifying header is written to the stream first. To avoid writing the header,
construct a `Serializer` and use it as the first argument to `serialize` instead.
See also [`Serialization.writeheader`](@ref).
"""
function serialize(s::IO, x)
    ss = Serializer(s)
    writeheader(ss)
    serialize(ss, x)
end

"""
    serialize(filename::AbstractString, value)

Open a file and serialize the given value to it.

!!! compat "Julia 1.1"
    This method is available as of Julia 1.1.
"""
serialize(filename::AbstractString, x) = open(io->serialize(io, x), filename, "w")

## deserializing values ##

"""
    deserialize(stream)

Read a value written by [`serialize`](@ref). `deserialize` assumes the binary data read from
`stream` is correct and has been serialized by a compatible implementation of [`serialize`](@ref).
It has been designed with simplicity and performance as a goal and does not validate
the data read. Malformed data can result in process termination. The caller has to ensure
the integrity and correctness of data read from `stream`.
"""
deserialize(s::IO) = deserialize(Serializer(s))

"""
    deserialize(filename::AbstractString)

Open a file and deserialize its contents.

!!! compat "Julia 1.1"
    This method is available as of Julia 1.1.
"""
deserialize(filename::AbstractString) = open(deserialize, filename)

function deserialize(s::AbstractSerializer)
    handle_deserialize(s, Int32(read(s.io, UInt8)::UInt8))
end

function deserialize_cycle(s::AbstractSerializer, @nospecialize(x))
    slot = pop!(s.pending_refs)
    s.table[slot] = x
    nothing
end

# optimized version of:
#     slot = s.counter; s.counter += 1
#     push!(s.pending_refs, slot)
#     slot = pop!(s.pending_refs)
#     s.table[slot] = x
function resolve_ref_immediately(s::AbstractSerializer, @nospecialize(x))
    s.table[s.counter] = x
    s.counter += 1
    nothing
end

function gettable(s::AbstractSerializer, id::Int)
    get(s.table, id) do
        errmsg = """Inconsistent Serializer state when deserializing.
            Attempt to access internal table with key $id failed.

            This might occur if the Serializer contexts when serializing and deserializing are inconsistent.
            In particular, if multiple serialize calls use the same Serializer object then
            the corresponding deserialize calls should also use the same Serializer object.
        """
        error(errmsg)
    end
end

# deserialize_ is an internal function to dispatch on the tag
# describing the serialized representation. the number of
# representations is fixed, so deserialize_ does not get extended.
function handle_deserialize(s::AbstractSerializer, b::Int32)
    if b == 0
        return desertag(Int32(read(s.io, UInt8)::UInt8))
    end
    if b >= VALUE_TAGS
        return desertag(b)
    elseif b == TUPLE_TAG
        return deserialize_tuple(s, Int(read(s.io, UInt8)::UInt8))
    elseif b == SHORTBACKREF_TAG
        id = read(s.io, UInt16)::UInt16
        return gettable(s, Int(id))
    elseif b == BACKREF_TAG
        id = read(s.io, Int32)::Int32
        return gettable(s, Int(id))
    elseif b == ARRAY_TAG
        return deserialize_array(s)
    elseif b == DATATYPE_TAG
        return deserialize_datatype(s, false)
    elseif b == FULL_DATATYPE_TAG
        return deserialize_datatype(s, true)
    elseif b == WRAPPER_DATATYPE_TAG
        tname = deserialize(s)::Core.TypeName
        return unwrap_unionall(tname.wrapper)
    elseif b == OBJECT_TAG
        t = deserialize(s)
        if t === Missing
            return missing
        end
        return deserialize(s, t)
    elseif b == REF_OBJECT_TAG
        slot = s.counter; s.counter += 1
        push!(s.pending_refs, slot)
        t = deserialize(s)
        return deserialize(s, t)
    elseif b == SHARED_REF_TAG
        slot = s.counter; s.counter += 1
        obj = deserialize(s)
        s.table[slot] = obj
        return obj
    elseif b == SYMBOL_TAG
        return deserialize_symbol(s, Int(read(s.io, UInt8)::UInt8))
    elseif b == SHORTINT64_TAG
        return Int64(read(s.io, Int32)::Int32)
    elseif b == EXPR_TAG
        return deserialize_expr(s, Int(read(s.io, UInt8)::UInt8))
    elseif b == MODULE_TAG
        return deserialize_module(s)
    elseif b == STRING_TAG
        return deserialize_string(s, Int(read(s.io, UInt8)::UInt8))
    elseif b == LONGSTRING_TAG
        return deserialize_string(s, Int(read(s.io, Int64)::Int64))
    elseif b == SIMPLEVECTOR_TAG
        return deserialize_svec(s)
    elseif b == GLOBALREF_TAG
        return GlobalRef(deserialize(s)::Module, deserialize(s)::Symbol)
    elseif b == FULL_GLOBALREF_TAG
        ty = deserialize(s)
        tn = unwrap_unionall(ty).name
        return GlobalRef(tn.module, tn.name)
    elseif b == LONGTUPLE_TAG
        return deserialize_tuple(s, Int(read(s.io, Int32)::Int32))
    elseif b == LONGEXPR_TAG
        return deserialize_expr(s, Int(read(s.io, Int32)::Int32))
    elseif b == LONGBACKREF_TAG
        id = read(s.io, Int64)::Int64
        return gettable(s, Int(id))
    elseif b == LONGSYMBOL_TAG
        return deserialize_symbol(s, Int(read(s.io, Int32)::Int32))
    elseif b == HEADER_TAG
        readheader(s)
        return deserialize(s)
    elseif b == INT8_TAG
        return read(s.io, Int8)
    elseif b == INT8_TAG+1
        return read(s.io, UInt8)
    elseif b == INT8_TAG+2
        return read(s.io, Int16)
    elseif b == INT8_TAG+3
        return read(s.io, UInt16)
    elseif b == INT32_TAG
        return read(s.io, Int32)
    elseif b == INT8_TAG+5
        return read(s.io, UInt32)
    elseif b == INT64_TAG
        return read(s.io, Int64)
    elseif b == INT8_TAG+7
        return read(s.io, UInt64)
    elseif b == INT8_TAG+8
        return read(s.io, Int128)
    elseif b == INT8_TAG+9
        return read(s.io, UInt128)
    elseif b == INT8_TAG+10
        return read(s.io, Float16)
    elseif b == INT8_TAG+11
        return read(s.io, Float32)
    elseif b == INT8_TAG+12
        return read(s.io, Float64)
    elseif b == INT8_TAG+13
        return read(s.io, Char)
    elseif b == IDDICT_TAG
        slot = s.counter; s.counter += 1
        push!(s.pending_refs, slot)
        t = deserialize(s)
        return deserialize_dict(s, t)
    end
    t = desertag(b)::DataType
    if t.mutable && length(t.types) > 0  # manual specialization of fieldcount
        slot = s.counter; s.counter += 1
        push!(s.pending_refs, slot)
    end
    return deserialize(s, t)
end

function deserialize_symbol(s::AbstractSerializer, len::Int)
    str = Base._string_n(len)
    unsafe_read(s.io, pointer(str), len)
    sym = Symbol(str)
    if len > 7
        resolve_ref_immediately(s, sym)
    end
    return sym
end

deserialize_tuple(s::AbstractSerializer, len) = ntupleany(i->deserialize(s), len)

function deserialize_svec(s::AbstractSerializer)
    n = read(s.io, Int32)
    svec(Any[ deserialize(s) for i=1:n ]...)
end

function deserialize_module(s::AbstractSerializer)
    mkey = deserialize(s)
    if isa(mkey, Tuple)
        # old version, TODO: remove
        if mkey === ()
            return Main
        end
        m = Base.root_module(mkey[1])
        for i = 2:length(mkey)
            m = getfield(m, mkey[i])::Module
        end
    else
        name = String(deserialize(s)::Symbol)
        pkg = (mkey === nothing) ? Base.PkgId(name) : Base.PkgId(Base.UUID(mkey), name)
        m = Base.root_module(pkg)
        mname = deserialize(s)
        while mname !== ()
            m = getfield(m, mname)::Module
            mname = deserialize(s)
        end
    end
    return m
end

function deserialize(s::AbstractSerializer, ::Type{Method})
    lnumber = read(s.io, UInt64)
    meth = lookup_object_number(s, lnumber)
    if meth !== nothing
        meth = meth::Method
        makenew = false
    else
        meth = ccall(:jl_new_method_uninit, Ref{Method}, (Any,), Main)
        makenew = true
    end
    deserialize_cycle(s, meth)
    mod = deserialize(s)::Module
    name = deserialize(s)::Symbol
    file = deserialize(s)::Symbol
    line = deserialize(s)::Int32
    sig = deserialize(s)::Type
    syms = deserialize(s)
    if syms isa SimpleVector
        # < v1.2
        _ambig = deserialize(s)
    else
        slot_syms = syms::String
    end
    nargs = deserialize(s)::Int32
    isva = deserialize(s)::Bool
    is_for_opaque_closure = false
    aggressive_constprop = false
    template_or_is_opaque = deserialize(s)
    if isa(template_or_is_opaque, Bool)
        is_for_opaque_closure = template_or_is_opaque
        if format_version(s) >= 14
            aggressive_constprop = deserialize(s)::Bool
        end
        template = deserialize(s)
    else
        template = template_or_is_opaque
    end
    generator = deserialize(s)
    recursion_relation = nothing
    if format_version(s) >= 15
        recursion_relation = deserialize(s)
    end
    if makenew
        meth.module = mod
        meth.name = name
        meth.file = file
        meth.line = line
        meth.sig = sig
        meth.nargs = nargs
        meth.isva = isva
        meth.is_for_opaque_closure = is_for_opaque_closure
        meth.aggressive_constprop = aggressive_constprop
        if template !== nothing
            # TODO: compress template
            meth.source = template::CodeInfo
            meth.pure = template.pure
            if !@isdefined(slot_syms)
                slot_syms = ccall(:jl_compress_argnames, Ref{String}, (Any,), meth.source.slotnames)
            end
        end
        meth.slot_syms = slot_syms
        if generator !== nothing
            meth.generator = generator
        end
        if recursion_relation !== nothing
            meth.recursion_relation = recursion_relation
        end
        if !is_for_opaque_closure
            mt = ccall(:jl_method_table_for, Any, (Any,), sig)
            if mt !== nothing && nothing === ccall(:jl_methtable_lookup, Any, (Any, Any, UInt), mt, sig, typemax(UInt))
                ccall(:jl_method_table_insert, Cvoid, (Any, Any, Ptr{Cvoid}), mt, meth, C_NULL)
            end
        end
        remember_object(s, meth, lnumber)
    end
    return meth
end

function deserialize(s::AbstractSerializer, ::Type{Core.MethodInstance})
    linfo = ccall(:jl_new_method_instance_uninit, Ref{Core.MethodInstance}, (Ptr{Cvoid},), C_NULL)
    deserialize_cycle(s, linfo)
    tag = Int32(read(s.io, UInt8)::UInt8)
    if tag != UNDEFREF_TAG
        linfo.uninferred = handle_deserialize(s, tag)::CodeInfo
    end
    tag = Int32(read(s.io, UInt8)::UInt8)
    if tag != UNDEFREF_TAG
        # for reading files prior to v1.2
        handle_deserialize(s, tag)
    end
    linfo.sparam_vals = deserialize(s)::SimpleVector
    _rettype = deserialize(s)  # for backwards compat
    linfo.specTypes = deserialize(s)
    linfo.def = deserialize(s)
    return linfo
end

function deserialize(s::AbstractSerializer, ::Type{Core.LineInfoNode})
    mod = deserialize(s)
    if mod isa Module
        method = deserialize(s)
    else
        # files post v1.2 and pre v1.6 are broken
        method = mod
        mod = Main
    end
    return Core.LineInfoNode(mod, method, deserialize(s)::Symbol, deserialize(s)::Int, deserialize(s)::Int)
end

function deserialize(s::AbstractSerializer, ::Type{PhiNode})
    edges = deserialize(s)
    if edges isa Vector{Any}
        edges = Vector{Int32}(edges)
    end
    values = deserialize(s)::Vector{Any}
    return PhiNode(edges, values)
end

function deserialize(s::AbstractSerializer, ::Type{CodeInfo})
    ci = ccall(:jl_new_code_info_uninit, Ref{CodeInfo}, ())
    deserialize_cycle(s, ci)
    code = deserialize(s)::Vector{Any}
    ci.code = code
    # allow older-style IR with return and gotoifnot Exprs
    for i in 1:length(code)
        stmt = code[i]
        if isa(stmt, Expr)
            ex = stmt::Expr
            if ex.head === :return
                code[i] = ReturnNode(isempty(ex.args) ? nothing : ex.args[1])
            elseif ex.head === :gotoifnot
                code[i] = GotoIfNot(ex.args[1], ex.args[2])
            end
        end
    end
    ci.codelocs = deserialize(s)::Vector{Int32}
    _x = deserialize(s)
    if _x isa Array || _x isa Int
        pre_12 = false
        ci.ssavaluetypes = _x
    else
        pre_12 = true
        # < v1.2
        ci.method_for_inference_limit_heuristics = _x
        ci.ssavaluetypes = deserialize(s)
        ci.linetable = deserialize(s)
    end
    ci.ssaflags = deserialize(s)
    if pre_12
        ci.slotflags = deserialize(s)
    else
        ci.method_for_inference_limit_heuristics = deserialize(s)
        ci.linetable = deserialize(s)
    end
    ci.slotnames = deserialize(s)
    if !pre_12
        ci.slotflags = deserialize(s)
        ci.slottypes = deserialize(s)
        ci.rettype = deserialize(s)
        ci.parent = deserialize(s)
        world_or_edges = deserialize(s)
        pre_13 = isa(world_or_edges, Integer)
        if pre_13
            ci.min_world = world_or_edges
        else
            ci.edges = world_or_edges
            ci.min_world = reinterpret(UInt, deserialize(s))
            ci.max_world = reinterpret(UInt, deserialize(s))
        end
    end
    ci.inferred = deserialize(s)
    ci.inlineable = deserialize(s)
    ci.propagate_inbounds = deserialize(s)
    ci.pure = deserialize(s)
    if format_version(s) >= 14
        ci.aggressive_constprop = deserialize(s)::Bool
    end
    return ci
end

if Int === Int64
const OtherInt = Int32
else
const OtherInt = Int64
end

function deserialize_array(s::AbstractSerializer)
    slot = s.counter; s.counter += 1
    d1 = deserialize(s)
    if isa(d1, Type)
        elty = d1
        d1 = deserialize(s)
    else
        elty = UInt8
    end
    if isa(d1, Int32) || isa(d1, Int64)
        if elty !== Bool && isbitstype(elty)
            a = Vector{elty}(undef, d1)
            s.table[slot] = a
            return read!(s.io, a)
        end
        dims = (Int(d1),)
    elseif d1 isa Dims
        dims = d1::Dims
    else
        dims = convert(Dims, d1::Tuple{Vararg{OtherInt}})::Dims
    end
    if isbitstype(elty)
        n = prod(dims)::Int
        if elty === Bool && n > 0
            A = Array{Bool, length(dims)}(undef, dims)
            i = 1
            while i <= n
                b = read(s.io, UInt8)::UInt8
                v = (b >> 7) != 0
                count = b & 0x7f
                nxt = i + count
                while i < nxt
                    A[i] = v
                    i += 1
                end
            end
        else
            A = read!(s.io, Array{elty}(undef, dims))
        end
        s.table[slot] = A
        return A
    end
    A = Array{elty, length(dims)}(undef, dims)
    s.table[slot] = A
    sizehint!(s.table, s.counter + div(length(A)::Int,4))
    deserialize_fillarray!(A, s)
    return A
end

function deserialize_fillarray!(A::Array{T}, s::AbstractSerializer) where {T}
    for i = eachindex(A)
        tag = Int32(read(s.io, UInt8)::UInt8)
        if tag != UNDEFREF_TAG
            @inbounds A[i] = handle_deserialize(s, tag)
        end
    end
    return A
end

function deserialize_expr(s::AbstractSerializer, len)
    e = Expr(:temp)
    resolve_ref_immediately(s, e)
    e.head = deserialize(s)::Symbol
    e.args = Any[ deserialize(s) for i = 1:len ]
    e
end

module __deserialized_types__ end

function deserialize(s::AbstractSerializer, ::Type{Core.TypeName})
    number = read(s.io, UInt64)
    return deserialize_typename(s, number)
end

function deserialize_typename(s::AbstractSerializer, number)
    name = deserialize(s)::Symbol
    tn = lookup_object_number(s, number)
    if tn !== nothing
        makenew = false
    else
        # reuse the same name for the type, if possible, for nicer debugging
        tn_name = isdefined(__deserialized_types__, name) ? gensym() : name
        tn = ccall(:jl_new_typename_in, Ref{Core.TypeName}, (Any, Any),
                   tn_name, __deserialized_types__)
        makenew = true
    end
    remember_object(s, tn, number)
    deserialize_cycle(s, tn)

    names = deserialize(s)::SimpleVector
    super = deserialize(s)::Type
    parameters = deserialize(s)::SimpleVector
    types = deserialize(s)::SimpleVector
    has_instance = deserialize(s)::Bool
    abstr = deserialize(s)::Bool
    mutabl = deserialize(s)::Bool
    ninitialized = deserialize(s)::Int32

    if makenew
        tn.names = names
        # TODO: there's an unhanded cycle in the dependency graph at this point:
        # while deserializing super and/or types, we may have encountered
        # tn.wrapper and throw UndefRefException before we get to this point
        ndt = ccall(:jl_new_datatype, Any, (Any, Any, Any, Any, Any, Any, Cint, Cint, Cint),
                    tn, tn.module, super, parameters, names, types,
                    abstr, mutabl, ninitialized)
        tn.wrapper = ndt.name.wrapper
        ccall(:jl_set_const, Cvoid, (Any, Any, Any), tn.module, tn.name, tn.wrapper)
        ty = tn.wrapper
        if has_instance && !isdefined(ty, :instance)
            # use setfield! directly to avoid `fieldtype` lowering expecting to see a Singleton object already on ty
            Core.setfield!(ty, :instance, ccall(:jl_new_struct, Any, (Any, Any...), ty))
        end
    end

    tag = Int32(read(s.io, UInt8)::UInt8)
    if tag != UNDEFREF_TAG
        mtname = handle_deserialize(s, tag)
        defs = deserialize(s)
        maxa = deserialize(s)::Int
        if makenew
            tn.mt = ccall(:jl_new_method_table, Any, (Any, Any), name, tn.module)
            if !isempty(parameters)
                tn.mt.offs = 0
            end
            tn.mt.name = mtname
            tn.mt.max_args = maxa
            for def in defs
                if isdefined(def, :sig)
                    ccall(:jl_method_table_insert, Cvoid, (Any, Any, Ptr{Cvoid}), tn.mt, def, C_NULL)
                end
            end
        end
        tag = Int32(read(s.io, UInt8)::UInt8)
        if tag != UNDEFREF_TAG
            kws = handle_deserialize(s, tag)
            if makenew
                tn.mt.kwsorter = kws
            end
        end
    elseif makenew
        tn.mt = Symbol.name.mt
    end
    return tn::Core.TypeName
end

function deserialize_datatype(s::AbstractSerializer, full::Bool)
    slot = s.counter; s.counter += 1
    if full
        tname = deserialize(s)::Core.TypeName
        ty = tname.wrapper
    else
        name = deserialize(s)::Symbol
        mod = deserialize(s)::Module
        ty = getfield(mod,name)
    end
    if isa(ty,DataType) && isempty(ty.parameters)
        t = ty
    else
        np = Int(read(s.io, Int32)::Int32)
        if np == 0
            t = unwrap_unionall(ty)
        elseif ty === Tuple
            # note np==0 has its own tag
            if np == 1
                t = Tuple{deserialize(s)}
            elseif np == 2
                t = Tuple{deserialize(s), deserialize(s)}
            elseif np == 3
                t = Tuple{deserialize(s), deserialize(s), deserialize(s)}
            elseif np == 4
                t = Tuple{deserialize(s), deserialize(s), deserialize(s), deserialize(s)}
            else
                t = Tuple{Any[ deserialize(s) for i=1:np ]...}
            end
        else
            t = ty
            for i = 1:np
                t = t{deserialize(s)}
            end
        end
    end
    s.table[slot] = t
    return t
end

function deserialize(s::AbstractSerializer, ::Type{UnionAll})
    form = read(s.io, UInt8)
    if form == 0
        var = deserialize(s)
        body = deserialize(s)
        return UnionAll(var, body)
    else
        n = read(s.io, Int16)
        t = deserialize(s)::DataType
        w = t.name.wrapper
        k = 0
        while isa(w, UnionAll)
            w = w.body
            k += 1
        end
        w = t.name.wrapper
        k -= n
        while k > 0
            w = w.body
            k -= 1
        end
        return w
    end
end

function deserialize(s::AbstractSerializer, ::Type{Task})
    t = Task(()->nothing)
    deserialize_cycle(s, t)
    t.code = deserialize(s)
    t.storage = deserialize(s)
    state = deserialize(s)
    if state === :runnable
        t._state = Base.task_state_runnable
    elseif state === :done
        t._state = Base.task_state_done
    elseif state === :failed
        t._state = Base.task_state_failed
    else
        @assert false
    end
    t.result = deserialize(s)
    exc = deserialize(s)
    if exc === nothing
        t._isexception = false
    elseif exc isa Bool
        t._isexception = exc
    else
        t._isexception = true
        t.result = exc
    end
    t
end

function deserialize_string(s::AbstractSerializer, len::Int)
    out = ccall(:jl_alloc_string, Ref{String}, (Csize_t,), len)
    unsafe_read(s.io, pointer(out), len)
    return out
end

# default DataType deserializer
function deserialize(s::AbstractSerializer, t::DataType)
    nf = length(t.types)
    if nf == 0 && t.size > 0
        # bits type
        return read(s.io, t)
    elseif t.mutable
        x = ccall(:jl_new_struct_uninit, Any, (Any,), t)
        deserialize_cycle(s, x)
        for i in 1:nf
            tag = Int32(read(s.io, UInt8)::UInt8)
            if tag != UNDEFREF_TAG
                ccall(:jl_set_nth_field, Cvoid, (Any, Csize_t, Any), x, i-1, handle_deserialize(s, tag))
            end
        end
        return x
    elseif nf == 0
        return ccall(:jl_new_struct_uninit, Any, (Any,), t)
    else
        na = nf
        vflds = Vector{Any}(undef, nf)
        for i in 1:nf
            tag = Int32(read(s.io, UInt8)::UInt8)
            if tag != UNDEFREF_TAG
                f = handle_deserialize(s, tag)
                na >= i && (vflds[i] = f)
            else
                na >= i && (na = i - 1) # rest of tail must be undefined values
            end
        end
        return ccall(:jl_new_structv, Any, (Any, Ptr{Any}, UInt32), t, vflds, na)
    end
end

function deserialize_dict(s::AbstractSerializer, T::Type{<:AbstractDict})
    n = read(s.io, Int32)
    t = T(); sizehint!(t, n)
    deserialize_cycle(s, t)
    for i = 1:n
        k = deserialize(s)
        v = deserialize(s)
        t[k] = v
    end
    return t
end

function deserialize(s::AbstractSerializer, T::Type{Dict{K,V}}) where {K,V}
    return deserialize_dict(s, T)
end

deserialize(s::AbstractSerializer, ::Type{BigInt}) = parse(BigInt, deserialize(s), base = 62)

function deserialize(s::AbstractSerializer, t::Type{Regex})
    pattern = deserialize(s)
    compile_options = deserialize(s)
    match_options = deserialize(s)
    return Regex(pattern, compile_options, match_options)
end

## StackTraces

# provide a custom serializer that skips attempting to serialize the `outer_linfo`
# which is likely to contain complex references, types, and module references
# that may not exist on the receiver end
function serialize(s::AbstractSerializer, frame::Base.StackTraces.StackFrame)
    serialize_type(s, typeof(frame))
    serialize(s, frame.func)
    serialize(s, frame.file)
    write(s.io, frame.line)
    write(s.io, frame.from_c)
    write(s.io, frame.inlined)
    write(s.io, frame.pointer)
    nothing
end

function deserialize(s::AbstractSerializer, ::Type{Base.StackTraces.StackFrame})
    func = deserialize(s)
    file = deserialize(s)
    line = read(s.io, Int)
    from_c = read(s.io, Bool)
    inlined = read(s.io, Bool)
    pointer = read(s.io, UInt64)
    return Base.StackTraces.StackFrame(func, file, line, nothing, from_c, inlined, pointer)
end

end

max_stars_repo_name: steveno/julia

max_issues_repo_name: steveno/julia

max_forks_repo_name: steveno/julia


text:
module RocketRefCountOperatorTest

using Test
using Rocket

include("../test_helpers.jl")

@testset "operator: ref_count()" begin

    println("Testing: operator ref_count()")

    run_proxyshowcheck("RefCount", ref_count(), args = (custom_source = Rocket.connectable(Subject(Any), never(Any)), check_subscription = true, ))

    @testset begin
        source  = from(1:5) |> multicast(Subject(Int)) |> ref_count()

        actor1 = keep(Int)
        actor2 = keep(Int)

        subscription1 = subscribe!(source, actor1)
        subscription2 = subscribe!(source, actor2)

        unsubscribe!(subscription1)
        unsubscribe!(subscription2)

        @test actor1.values == [ 1, 2, 3, 4, 5 ]
        @test actor2.values == [ ]
    end

    @testset begin
        subject = Subject(Int)
        source  = subject |> multicast(Subject(Int)) |> ref_count()

        actor1 = keep(Int)
        actor2 = keep(Int)
        actor3 = keep(Int)

        subscription3 = subscribe!(subject, actor3)

        next!(subject, 1)

        @test actor1.values == [ ]
        @test actor2.values == [ ]
        @test actor3.values == [ 1 ]

        subscription1 = subscribe!(source, actor1)
        subscription2 = subscribe!(source, actor2)

        @test actor1.values == [ ]
        @test actor2.values == [ ]
        @test actor3.values == [ 1 ]

        next!(subject, 2)

        @test actor1.values == [ 2 ]
        @test actor2.values == [ 2 ]
        @test actor3.values == [ 1, 2 ]

        unsubscribe!(subscription1)

        next!(subject, 3)

        @test actor1.values == [ 2 ]
        @test actor2.values == [ 2, 3 ]
        @test actor3.values == [ 1, 2, 3 ]

        unsubscribe!(subscription2)
        subscription1 = subscribe!(source, actor1)

        next!(subject, 4)

        @test actor1.values == [ 2, 4 ]
        @test actor2.values == [ 2, 3 ]
        @test actor3.values == [ 1, 2, 3, 4 ]

        unsubscribe!(subscription1)
        unsubscribe!(subscription2)

        next!(subject, 5)

        @test actor1.values == [ 2, 4 ]
        @test actor2.values == [ 2, 3 ]
        @test actor3.values == [ 1, 2, 3, 4, 5]
    end

end

end

ext: jl


text:
function [D,J,JInv,X]=JacobiSphere(ksi,F,Grid)
Rad=Grid.Rad;
X1=Grid.Nodes(F.N(1)).P(1)...
  +(Grid.Nodes(F.N(2)).P(1)-Grid.Nodes(F.N(1)).P(1))*ksi(1)...
  +(Grid.Nodes(F.N(4)).P(1)-Grid.Nodes(F.N(1)).P(1))*ksi(2)...
  +(Grid.Nodes(F.N(3)).P(1)-Grid.Nodes(F.N(4)).P(1)-Grid.Nodes(F.N(2)).P(1)+Grid.Nodes(F.N(1)).P(1))*ksi(1)*ksi(2);
X2=Grid.Nodes(F.N(1)).P(2)...
  +(Grid.Nodes(F.N(2)).P(2)-Grid.Nodes(F.N(1)).P(2))*ksi(1)...
  +(Grid.Nodes(F.N(4)).P(2)-Grid.Nodes(F.N(1)).P(2))*ksi(2)...
  +(Grid.Nodes(F.N(3)).P(2)-Grid.Nodes(F.N(4)).P(2)-Grid.Nodes(F.N(2)).P(2)+Grid.Nodes(F.N(1)).P(2))*ksi(1)*ksi(2);
X3=Grid.Nodes(F.N(1)).P(3)...
  +(Grid.Nodes(F.N(2)).P(3)-Grid.Nodes(F.N(1)).P(3))*ksi(1)...
  +(Grid.Nodes(F.N(4)).P(3)-Grid.Nodes(F.N(1)).P(3))*ksi(2)...
  +(Grid.Nodes(F.N(3)).P(3)-Grid.Nodes(F.N(4)).P(3)-Grid.Nodes(F.N(2)).P(3)+Grid.Nodes(F.N(1)).P(3))*ksi(1)*ksi(2);
f=Rad*(X1^2+X2^2+X3^2)^(-3/2);

dx1dX1=f*(X2^2+X3^2);
dx1dX2=-f*X1*X2;
dx1dX3=-f*X1*X3;
dx2dX1=dx1dX2;
dx2dX2=f*(X1^2+X3^2);
dx2dX3=-f*X2*X3;
dx3dX1=dx1dX3;
dx3dX2=dx2dX3;
dx3dX3=f*(X1^2+X2^2);

% dx1dX1=1;
% dx1dX2=-0;
% dx1dX3=0;
% dx2dX1=0;
% dx2dX2=1;
% dx2dX3=0;
% dx3dX1=0;
% dx3dX2=0;
% dx3dX3=1;

dX1dksi1=(Grid.Nodes(F.N(2)).P(1)-Grid.Nodes(F.N(1)).P(1))...
  +(Grid.Nodes(F.N(3)).P(1)-Grid.Nodes(F.N(4)).P(1)...
  -Grid.Nodes(F.N(2)).P(1)+Grid.Nodes(F.N(1)).P(1))*ksi(2);
dX1dksi2=(Grid.Nodes(F.N(4)).P(1)-Grid.Nodes(F.N(1)).P(1))...
  +(Grid.Nodes(F.N(3)).P(1)-Grid.Nodes(F.N(4)).P(1)...
  -Grid.Nodes(F.N(2)).P(1)+Grid.Nodes(F.N(1)).P(1))*ksi(1);
dX2dksi1=(Grid.Nodes(F.N(2)).P(2)-Grid.Nodes(F.N(1)).P(2))...
  +(Grid.Nodes(F.N(3)).P(2)-Grid.Nodes(F.N(4)).P(2)...
  -Grid.Nodes(F.N(2)).P(2)+Grid.Nodes(F.N(1)).P(2))*ksi(2);
dX2dksi2=(Grid.Nodes(F.N(4)).P(2)-Grid.Nodes(F.N(1)).P(2))...
  +(Grid.Nodes(F.N(3)).P(2)-Grid.Nodes(F.N(4)).P(2)...
  -Grid.Nodes(F.N(2)).P(2)+Grid.Nodes(F.N(1)).P(2))*ksi(1);
dX3dksi1=(Grid.Nodes(F.N(2)).P(3)-Grid.Nodes(F.N(1)).P(3))...
  +(Grid.Nodes(F.N(3)).P(3)-Grid.Nodes(F.N(4)).P(3)...
  -Grid.Nodes(F.N(2)).P(3)+Grid.Nodes(F.N(1)).P(3))*ksi(2);
dX3dksi2=(Grid.Nodes(F.N(4)).P(3)-Grid.Nodes(F.N(1)).P(3))...
  +(Grid.Nodes(F.N(3)).P(3)-Grid.Nodes(F.N(4)).P(3)...
  -Grid.Nodes(F.N(2)).P(3)+Grid.Nodes(F.N(1)).P(3))*ksi(1);

J=zeros(3,2);

J(1,1)=dx1dX1*dX1dksi1+dx1dX2*dX2dksi1+dx1dX3*dX3dksi1;
J(2,1)=dx2dX1*dX1dksi1+dx2dX2*dX2dksi1+dx2dX3*dX3dksi1;
J(3,1)=dx3dX1*dX1dksi1+dx3dX2*dX2dksi1+dx3dX3*dX3dksi1;
J(1,2)=dx1dX1*dX1dksi2+dx1dX2*dX2dksi2+dx1dX3*dX3dksi2;
J(2,2)=dx2dX1*dX1dksi2+dx2dX2*dX2dksi2+dx2dX3*dX3dksi2;
J(3,2)=dx3dX1*dX1dksi2+dx3dX2*dX2dksi2+dx3dX3*dX3dksi2;

D=norm(cross(J(:,1),J(:,2)),2);
if nargout > 2
  X=[X1 X2 X3]*(Rad/sqrt((X1^2+X2^2+X3^2)));   X=[X1 X2 X3];
end
JInv=inv(J'*J)*J'*sqrt(det(J'*J));
end


text:
#!/usr/bin/env julia

#######################################################################
# Generation of a list of the top N most starred Julia packages
#######################################################################

#--------------------------------------------------
# Imports
#--------------------------------------------------

include("src/JuliaPkgsList.jl")
using Main.JuliaPkgsList

using ArgParse

#--------------------------------------------------
# Command line arguments
#--------------------------------------------------

# → Dict (arguments)
function parse_command_line_args()
    argsStr = ArgParseSettings(
        description = "Generates the list of top <pkgnum> most starred Julia packages using data from Julia Hub. Use -h or --help for usage details."
    )
    @add_arg_table! argsStr begin
        "pkgnum"
            help = "number of top most starred packages of interest (0 means all packages)"
            arg_type = Int
            required = true

        "--reload", "-r"
            help = "flag specifying if packages information must be reloaded"
            action = :store_true
        "--pkginfo", "-p"
            help = "JSON file with packages information"
            arg_type = String
            default = JuliaPkgsList.PKGS_INFO_FILE

        "--name", "-n"
            help = "flag specifying if the output should be packages' names instead of repositories"
            action = :store_true
        "--includeversion"
            help = "flag specifying if the output should include the latest version of packages"
            action = :store_true

        "--excluded", "-e"
            help = "file with a list of repositories to exclude"
            arg_type = String
            default = JuliaPkgsList.EXCLUDED_REPOS_FILE
        "--noexclude"
            help = "flag specifying if all packages should be included (i.e. ignore the file with excluded)"
            action = :store_true

        "--show", "-s"
            help = "flag specifying if the output should be printed to console instead of a file"
            action = :store_true

        "--out", "-o"
            help = "output file with the top packages list"
            arg_type = String
            default = JuliaPkgsList.PKGS_LIST_FILE
        "--nopkgnum"
            help = "flag specifying that the output file should not be prepended with the number of packages"
            action = :store_true
    end
    parse_args(argsStr)
end

# All script parameters
const PARAMS = parse_command_line_args()

#--------------------------------------------------
# Main
#--------------------------------------------------

genTopPkgsList(
    PARAMS["pkgnum"], PARAMS["reload"]
;
    pkgsListOutputFile = PARAMS["out"],
    addPkgNum = !PARAMS["nopkgnum"],
    showName = PARAMS["name"],
    includeVersion = PARAMS["includeversion"],
    pkgsInfoFile = PARAMS["pkginfo"],
    excludedReposFile = PARAMS["excluded"],
    useStdOut = PARAMS["show"],
    useExcluded = !PARAMS["noexclude"]
)


text:
module Statsbomb

using DataFrames
using JSON: Parser, parsefile
using HTTP

struct SBLoader
    mode::String
    dir_or_url::String
    SBLoader(mode="remote", dir_or_url="https://raw.githubusercontent.com/statsbomb/open-data/master/data/") = new(mode, dir_or_url)
end

"""
load_competitions

Function to get the names of all available competitions in the free data

Arguments

loader

Returns

DataFrame
"""
function load_competitions(loader::SBLoader)

    if loader.mode == "local"

        comp_file = joinpath(loader.dir_or_url, "competitions.json")
        comp = parsefile(comp_file)
        return vcat(DataFrame.(comp)..., cols=:union)

    elseif loader.mode == "remote"

        resp = HTTP.get(loader.dir_or_url*"competitions.json")
        str = String(resp.body)
        comp = Parser.parse(str)
        return vcat(DataFrame.(comp)..., cols=:union)

    end
end

"""
load_matches

Function to get matches available from a particular season and competition

Arguments

loader
comp_id
season_id

Returns

DataFrame
"""
function load_matches(loader::SBLoader, comp_id, season_id)

    if loader.mode == "local"
        matches = parsefile(joinpath(joinpath(loader.dir_or_url*"\\matches", string(comp_id)), string(season_id)*".json"));

    elseif loader.mode == "remote"
        url = loader.dir_or_url*"/matches/"*string(comp_id)*"/"*string(season_id)*".json"
        resp = HTTP.get(url)
        str = String(resp.body)
        matches = Parser.parse(str)
    end

    df = vcat(DataFrame.(matches)..., cols=:union)
    df_home = vcat(DataFrame.(df."home_team")..., cols=:union)
    df_away = vcat(DataFrame.(df."away_team")..., cols=:union)

    df = hcat([df, df_home, df_away]..., makeunique=true)
    select!(df, Not([:home_team, :away_team]))

    df
end

"""
load_events

Function to get all events of a match

Arguments

loader
match_id

Returns

DataFrame
"""
function load_events(loader::SBLoader, match_id)

    if loader.mode == "local"

        data = parsefile(joinpath(loader.dir_or_url, "events\\"*string(match_id)*".json"))

    elseif loader.mode == "remote"

        url = loader.dir_or_url*"/events/"*string(match_id)*".json"
        resp = HTTP.get(url)
        str = String(resp.body)
        data = Parser.parse(str)
    end

    data = [Dict(isa(v, Array) ? k=>string(v) : k=> v for (k,v) in d) for d in data]     ##converting arrays to string to stop duplicating rows
    events = vcat(DataFrame.(data)..., cols=:union)
end

end ##end of module


text:
module HierarchicalAD

using DrWatson
@quickactivate
using Luxor
using Augmentor
using ImageIO, FileIO
using StatsBase
using Images
using BSON
using IPMeasures
using Flux
using CSV
using DataFrames
using MLDatasets
using Random
using Dates
using ValueHistories
using DataDeps
using PyCall
using UCI
using NearestNeighbors
using LinearAlgebra
using Statistics
using EvalMetrics
using Suppressor

include("digits.jl")
include("model_utils.jl")
include("disentanglement.jl")
include("training_utils.jl")
include("experiments.jl")
include("shapes2d.jl")
include("data.jl")
include("evaluation.jl")

include("knn.jl")
include("abstract_vlae.jl")
include("vlae.jl")
include("fvlae.jl")
include("had.jl")
include("conv_classifier.jl")

export VLAE, FVLAE, AbstractVLAE

end

max_stars_repo_name: vitskvara/HierarchicalAD

max_issues_repo_name: vitskvara/HierarchicalAD

max_forks_repo_name: vitskvara/HierarchicalAD


text:
export spectrum, eigenspace, jordan_normal_form, rational_canonical_form, companion_matrix

@doc Markdown.doc"""
    spectrum(M::MatElem{T}) where T <: FieldElem -> Dict{T, Int}

Returns the spectrum of a matrix, i.e. the set of eigenvalues of $M$ with multiplicities.
"""
function spectrum(M::MatElem{T}) where T <: FieldElem
  @assert issquare_matrix(M)
  K = base_ring(M)
  f = charpoly(M)
  fac = factor(f) #Should use roots! But needs to take into account
                  #multiplicities
  D = Dict{elem_type(K), Int}()
  for (g, v) in fac
    if degree(g) > 1
      continue
    end
    lambda = -divexact(coeff(g, 0), leading_coefficient(g))
    D[lambda] = v
  end
  return D
end

@doc Markdown.doc"""
    spectrum(M::MatElem{T}, L) where T <: FieldElem -> Dict{T, Int}

Returns the spectrum of a matrix over the field $L$, i.e. the set of eigenvalues of $M$ with multiplicities.
"""
function spectrum(M::MatElem{T}, L) where T <: FieldElem
  @assert issquare_matrix(M)
  M1 = change_base_ring(M, L)
  return spectrum(M1)
end

function issquare_matrix(M::MatElem)
  return ncols(M) == nrows(M)
end

eigvals(M::MatElem{T}) where T <: FieldElem = spectrum(M)
eigvals(M::MatElem{T}, L) where T <: FieldElem = spectrum(M, L)

@doc Markdown.doc"""
    eigenspace(M::MatElem{T}, lambda::T) where T <: FieldElem -> Vector{MatElem{T}}

Returns a basis of the eigenspace of $M$ with respect to the eigenvalues $\lambda$.
"""
function eigenspace(M::MatElem{T}, lambda::T) where T <: FieldElem
  @assert issquare_matrix(M)
  N = deepcopy(M)
  for i = 1:ncols(N)
    N[i, i] -= lambda
  end
  res = Hecke.left_kernel_basis(N)
  resvect = [matrix(base_ring(M), 1, ncols(M), x) for x in res]
  return resvect
end

function eigenspaces(M::MatElem{T}) where T<:Hecke.FieldElem

  S = spectrum(M)
  L = Dict{elem_type(base_ring(M)), typeof(M)}()
  for k in keys(S)
    push!(L, k => Hecke.vcat(Hecke.eigenspace(M, k)))
  end
  return L
end

function closure_with_pol(v::MatElem{T}, M::MatElem{T}) where T <: FieldElem
  i = 1
  E = rref(v)[2]
  w = v*M
  res = Hecke.cleanvect(E, w)
  while !iszero(res)
    v = vcat(v, w)
    E = vcat(E, res)
    rref!(E)
    i += 1
    w = w*M
    res = Hecke.cleanvect(E, w)
  end
  fl, c = Hecke.can_solve_with_solution(v, w, side = :left)
  @assert fl
  return v, c
end

function decompose_primary(M::MatElem{T}) where T <: FieldElem
  dimsum = 0
  i = 1
  K = base_ring(M)
  Kt, t = PolynomialRing(K, "t", cached = false)
  L = Dict{typeof(t), Vector{Tuple{typeof(M), typeof(M)}}}()
  S = zero_matrix(K, 0, ncols(M))
  while isempty(L) || sum(nrows(y[1]) for x in values(L) for y in x) < ncols(M)
    ei = zero_matrix(K, 1, ncols(M))
    ei[1, i] = 1
    res = Hecke.cleanvect(S, ei)
    while iszero(res)
      ei[1, i] = 0
      i += 1
      ei[1, i] = 1
      res = Hecke.cleanvect(S, ei)
    end
    C, coeffs = closure_with_pol(ei, M)
    S = vcat(S, C)
    rref!(S)
    cccc = [-coeffs[1, l] for l = 1:ncols(coeffs)]
    push!(cccc, K(1))
    m = Kt(cccc)
    p = one(Kt)
    for (g, s) in L
      f, fu = ppio(m, g)
      if isone(f)
        continue
      end
      w = ei * Hecke._subst(fu, M)
      p = p*f
      W = closure_with_pol(w, M)[1]
      W1 = W
      for (sub, gensub) in s
        W1 = vcat(W1, sub)
      end
      rk = rank(W1)
      if rk == nrows(W1) - nrows(W)
        continue
      end
      if rk == nrows(W1)
        #The space is independent from the one I had before
        push!(s, (W, w))
      else
        independent_from = Int[]
        intersect_not_contained = Int[]
        is_contained = []
        for j = 1:length(s)
          Ws = vcat(W, s[j][1])
          rk = rank(Ws)
          if rk == nrows(Ws)
            push!(independent_from, j)
          elseif rk == nrows(W)
            push!(is_contained, j)
          else
            push!(intersect_not_contained, j)
          end
        end
        if isempty(intersect_not_contained)
          #W containes some subspaces but it is independent from the others
          snew = Vector{Tuple{typeof(M), typeof(M)}}(undef, length(s) - length(is_contained) +1)
          indn = 1
          for z = 1:length(s)
            if !(z in is_contained)
              snew[indn] = s[z]
              indn += 1
            end
          end
          snew[end] = (W, w)
          L[g] = snew
        else
          snew = Vector{Tuple{typeof(M), typeof(M)}}(undef, length(s) - length(is_contained) -length(intersect_not_contained) +1)
          indn = 1
          newsub = W
          for z = 1:length(s)
            if !(z in is_contained) && !(z in intersect_not_contained)
              snew[indn] = s[z]
              indn += 1
            elseif (z in intersect_not_contained)
              newsub = vcat(s[z][1], newsub)
            end
          end
          rk = rref!(newsub)
          newsub1 = zero_matrix(K, rk, ncols(M))
          for r = 1:rk
            for cc = r:ncols(M)
              newsub1[r, cc] = newsub[r, cc]
            end
          end
          snew[end] = (newsub1, zero_matrix(K, 1, ncols(M)))
          L[g] = snew
        end
      end
    end
    ei = ei*Hecke._subst(p, M)
    m = divexact(m, p)
    if !isone(m)
      fac = factor(m)
      for (q, e) in fac
        qq = divexact(m, q^e)
        w = ei*Hecke._subst(qq, M)
        W = closure_with_pol(w, M)[1]
        L[q] = [(W, w)]
      end
    end
  end
  return L
end

function _copy_matrix_into_matrix(A::MatElem, i::Int, j::Int, B::MatElem)
  for k in 0:nrows(B) - 1
    for l in 0:ncols(B) - 1
      setindex!(A, B[k+1, l+1], i+k, j+l)
    end
  end
  return nothing
end

@doc Markdown.doc"""
    companion_matrix(p::PolyElem) -> MatElem

Returns the companion matrix of $p = \sum_{i=0}^n a_ix^i$, i.e. the matrix

    0 1  0  $\dots$  0
    0  0  1  $\dots$  0
    $\vdots$  $\vdots$  $\ddots$  $\vdots$  $\vdots$
    0  0  $\vdots$  0  1
    $-a_0$  $-a_1$  $-a_2$  $\dots$  $-a_{n-1}$

"""
function companion_matrix(p::PolyElem)
  K = base_ring(p)
  p1 = divexact(p, leading_coefficient(p))
  M = zero_matrix(K, degree(p), degree(p))
  for i = 1:degree(p)-1
    setindex!(M, one(K), i, i+1)
  end
  for i = 1:degree(p)
    setindex!(M, -coeff(p1, i-1), degree(p), i)
  end
  return M
end

function jordan_block(p::PolyElem, e::Int)
  K = base_ring(p)
  M = zero_matrix(K, degree(p)*e, degree(p)*e)
  C = companion_matrix(p)
  for i = 1:e
    _copy_matrix_into_matrix(M, 1+degree(p)*(i-1), 1+degree(p)*(i-1), C)
  end
  for i = 1:e-1
    M[i*degree(p), i*degree(p)+1] = K(1)
  end
  return M
end

@doc Markdown.doc"""
    issimilar(M::MatElem{T}, N::MatElem{T}) where T <: FieldElem -> Bool

Returns true if the matrices are similar (conjugated) and false otherwise.
"""
function issimilar(M::MatElem{T}, N::MatElem{T}) where T <: FieldElem
  CM = rational_canonical_form(M)[1]
  CN = rational_canonical_form(N)[1]
  return CM == CN
end

@doc Markdown.doc"""
    conjugating_matrix(M::MatElem{T}, N::MatElem{T}) where T <: FieldElem -> MatElem{T}

Returns a matrix $S$ such that $S\times N \times S^{-1} = M$.
"""
function conjugating_matrix(M::MatElem{T}, N::MatElem{T}) where T <: FieldElem
  CM, SM = rational_canonical_form(M)
  CN, SN = rational_canonical_form(N)
  if CM != CN
    error("Matrices are not similar!")
  end
  return inv(SM)*SN
end

################################################################################
#
#  Rational canonical form - Geck's algorithm
#
################################################################################

function maximal_vector(M::MatElem{T}, Kt, max_deg::Int = -1) where T <: FieldElem
  K = base_ring(M)
  v = zero_matrix(K, 1, ncols(M))
  bound = ncols(M)
  if max_deg != -1
    bound = min(max_deg, bound)
  end
  v[1, 1] = 1
  to_reduce, coeffs_vect = closure_with_pol(v, M)
  coeffs = T[-coeffs_vect[1, i] for i = 1:length(coeffs_vect)]
  push!(coeffs, K(1))
  mp = Kt(coeffs)
  if nrows(to_reduce) == bound
    return v, mp
  end
  rref!(to_reduce)
  for i = 2:nrows(M)
    v1 = zero_matrix(K, 1, ncols(M))
    v1[1, i] = 1
    res = Hecke.cleanvect(to_reduce, v1)
    if iszero(res)
      continue
    end
    C1, coeffs_vect2 = closure_with_pol(v1, M)
    coeffs2 = T[-coeffs_vect2[1, i] for i = 1:length(coeffs_vect2)]
    to_reduce = vcat(to_reduce, C1)
    rref!(to_reduce)
    push!(coeffs2, K(1))
    mp2 = Kt(coeffs2)
    if divides(mp, mp2)[1]
      continue
    end
    if nrows(C1) == bound
      return v1, mp2
    end
    f1, f2 = coprime_fact(mp, mp2)
    if f1 == mp && f2 == mp2
      v = v + v1
    else
      v = v*Hecke._subst(divexact(mp, f1), M) + v1*Hecke._subst(divexact(mp2, f2), M)
    end
    C, coeffs_vect = closure_with_pol(v, M)
    coeffs = T[-coeffs_vect[1, i] for i = 1:length(coeffs_vect)]
    push!(coeffs, K(1))
    mp = Kt(coeffs)
    if nrows(C) == bound
      return v, mp
    end
  end
  return v, mp
end

#Finds a1, b1 such that a1*b1 = a*b/gcd(a, b) and a1, b1 are coprime
function coprime_fact(a::PolyElem{T}, b::PolyElem{T}) where T <: FieldElem
  d = gcd(a, b)
  if isone(d)
    return a, b
  end
  atilde = divexact(a, d)
  btilde = divexact(b, d)
  b1 = gcd(btilde^(degree(b)), b)
  bprime = divexact(b, b1)
  a1 = bprime*atilde
  return a1, b1
end

function find_pivots(M::MatElem{T}) where T <: FieldElem
  pivots = Vector{Int}(undef, nrows(M))
  j = 1
  for i = 1:nrows(M)
    while iszero(M[i, j])
      j += 1
    end
    pivots[i] = j
  end
  return pivots
end

@doc Markdown.doc"""
    complete_to_basis(C::MatElem{T}) where T <: FieldElem -> MatElem{T}

Returns a matrix representing a basis of $K^n$ whose first elements are given by the columns of $C$.
"""
function complete_to_basis(C::MatElem{T}) where T <: FieldElem
  CEF = rref(C)[2]
  pivots = find_pivots(CEF)
  S = zero_matrix(base_ring(C), ncols(C), ncols(C))
  _copy_matrix_into_matrix(S, 1, 1, C)
  j = nrows(C)+1
  for i = 1:ncols(C)
    if i in pivots
      continue
    end
    S[j, i] = 1
    j += 1
  end
  return S
end

function find_invariant_complement(M::MatElem{T}, C::MatElem{T}) where T <: FieldElem
  #I follow Geck's approach. It would be nice if we were able to do the computation without
  # a base change
  #First, I have to complete to a basis C
  S = complete_to_basis(C)
  Sinv = inv(S)
  TB = S*M*Sinv
  ed = zero_matrix(base_ring(M), ncols(C), 1)
  ed[nrows(C), 1] = 1
  N = zero_matrix(base_ring(M), ncols(C), ncols(C))
  N[nrows(C), 1] = 1
  for i = 2:nrows(C)
    ed = TB*ed
    for j = 1:ncols(C)
      N[j, i] = ed[j, 1]
    end
  end
  k, K = kernel(N, side = :left)
  return view(K, 1:k, 1:ncols(K))*S
end

function _closure(M::MatElem{T}, v::MatElem{T}, d::Int) where T <: FieldElem
  res = zero_matrix(base_ring(M), d, ncols(M))
  for i = 1:ncols(M)
    res[1, i] = v[1, i]
  end
  w = v
  for i = 2:d
    w = w*M
    for j = 1:ncols(M)
      res[i, j] = w[1, j]
    end
  end
  return res
end

#M represents a linear map on a vector space and S gives a basis for an invariant subspace.
#The function returns a matrix representing the restriction of the linear map to the subspace
function restriction(M::MatElem{T}, S::MatElem{T}) where T <: FieldElem
  TR = S*M
  fl, R = Hecke.can_solve_with_solution(S, TR, side = :left)
  if !fl
    error("The subspace is not invariant!")
  end
  return view(R, 1:nrows(S), 1:nrows(S))
end

@doc Markdown.doc"""
    _rational_canonical_form_setup(M)

Returns minpolys, the basis transformation and the vectors generating the blocks of the
rational canonical form of $M$.
"""
function _rational_canonical_form_setup(M::MatElem{T}) where T <: FieldElem
  K = base_ring(M)
  Kt, t = PolynomialRing(K, "t", cached = false)
  v, mp = maximal_vector(M, Kt)
  pols = typeof(mp)[mp]
  basis_transf = typeof(M)[]
  gens = typeof(v)[v]
  C = _closure(M, v, degree(mp))
  push!(basis_transf, C)
  K = find_invariant_complement(M, C)
  N = restriction(M, K)
  bound = degree(mp)
  while sum(nrows(x) for x in basis_transf) < nrows(M)
    w, mp1 = maximal_vector(N, Kt, bound)
    push!(pols, mp1)
    push!(gens, w*K)
    subclos = _closure(N, w, degree(mp1))
    push!(basis_transf, subclos*K)
    if nrows(subclos) == nrows(N)
      break
    end
    K1 = find_invariant_complement(N, subclos)
    K = K1*K
    N = restriction(M, K)
    bound = degree(mp1)
  end
  return pols, basis_transf, gens
end

@doc Markdown.doc"""
    rational_canonical_form(M::MatElem{T}) where T <: FieldElem -> MatElem{T}, MatElem{T}

Returns matrices $C$ and $S$ such that $C = SMS^{-1}$ and $C$ is in rational canonical form.
"""
function rational_canonical_form(M::MatElem{T}) where T <: FieldElem
  @assert issquare_matrix(M)
  pols, basis_transf, gens = _rational_canonical_form_setup(M)
  N = zero(M)
  S = zero(M)
  ind = nrows(M)+1
  for i = 1:length(pols)
    C = companion_matrix(pols[i])
    ind -= nrows(C)
    _copy_matrix_into_matrix(N, ind, ind, C)
    _copy_matrix_into_matrix(S, ind, 1, basis_transf[i])
  end
  return N, S
end

function pre_factorization(pols::Vector)
  coprime_factors = Vector{typeof(pols[1])}(undef, length(pols))
  coprime_factors[length(pols)] = pols[length(pols)]
  for i = length(pols)-1:-1:1
    coprime_factors[i] = ppio(pols[i], pols[i+1])[2]
  end
  factors = Vector{typeof(coprime_factors[1])}()
  for i = 1:length(coprime_factors)
    fac = factor(coprime_factors[i])
    for (p, v) in fac
      push!(factors, divexact(p, leading_coefficient(p)))
    end
  end
  return factors
end

function factor_over(f::PolyElem{T}, l::Vector) where T <: FieldElem
  exps = Vector{Int}(undef, length(l))
  for i = 1:length(l)
    exps[i] = Int(valuation(f, l[i]))
  end
  return exps
end

function refine_for_jordan(pols::Vector, gens::Vector, M::MatElem)
  factors = pre_factorization(pols)
  gens_polys_mults = Vector{Tuple{typeof(gens[1]), typeof(pols[1]), Int}}()
  for i = length(pols):-1:1
    lef = factor_over(pols[i], factors)
    for j = 1:length(lef)
      if iszero(lef[j])
        continue
      end
      vj = gens[i]*Hecke._subst(divexact(pols[i], factors[j]^lef[j]), M)
      push!(gens_polys_mults, (vj, factors[j], lef[j]))
    end
  end
  return factors, gens_polys_mults
end

@doc Markdown.doc"""
    jordan_normal_form(M::MatElem{T}) where T <: FieldElem -> MatElem{T}, MatElem{T}

Returns matrices $J$ and $S$ such that $J = SMS^{-1}$ and $J$ is in Jordan normal form.
"""
function jordan_normal_form(M::MatElem{T}) where T <: FieldElem
  @assert issquare_matrix(M)
  pols, basis_transf, gens = _rational_canonical_form_setup(M)
  factors, gens_polys_mults = refine_for_jordan(pols, gens, M)
  J = zero(M)
  S = zero(M)
  ind = 1
  for i = 1:length(factors)
    for j = 1:length(gens_polys_mults)
      el = gens_polys_mults[j]
      if el[2] != factors[i]
        continue
      end
      JB = jordan_block(factors[i], el[3])
      _copy_matrix_into_matrix(J, ind, ind, JB)
      #Now, the transformation
      N = Hecke._subst(factors[i], M)
      w = el[1]
      for k = 1:el[3]
        aux = w
        for t = 1:degree(el[2])
          for s = 1:ncols(M)
            S[ind+(k-1)*degree(el[2])+t-1, s] = aux[1, s]
          end
          aux = aux*M
        end
        w = w*N
      end
      ind += nrows(JB)
    end
  end
  return J, S
end

@doc Markdown.doc"""
    jordan_decomposition(M::MatElem{T}) where T <:FieldElem -> MatElem{T}, MatElem{T}

Returns matrices $S$ and $N$ such that $N$ is nilpotent, $S$ is semisimple and $M = S+N$.
"""
function jordan_decomposition(M::MatElem{T}) where T <: FieldElem
  @assert issquare_matrix(M)
  K = base_ring(M)
  pols, basis_transf, gens = _rational_canonical_form_setup(M)
  factors, gens_polys_mults = refine_for_jordan(pols, gens, M)
  J = similar(M)
  B = similar(M)
  N = similar(M)
  ind = 1
  for i = 1:length(factors)
    for j = 1:length(gens_polys_mults)
      el = gens_polys_mults[j]
      if el[2] != factors[i]
        continue
      end
      #Now, the transformation
      N1 = Hecke._subst(factors[i], M)
      w = el[1]
      for k = 1:el[3]
        aux = w
        for t = 1:degree(el[2])
          for s = 1:ncols(M)
            B[ind+(k-1)*degree(el[2])+t-1, s] = aux[1, s]
          end
          aux = aux*M
        end
        w = w*N1
      end
      for c = 1:el[3]-1
        N[ind+c*degree(el[2])-1, ind+c*degree(el[2])] = one(K)
      end
      JB = jordan_block(factors[i], 1)
      for c = 1:el[3]
        _copy_matrix_into_matrix(J, ind, ind, JB)
        ind += degree(el[2])
      end
    end
  end
  Binv = inv(B)
  return Binv*J*B, Binv*N*B
end

@doc Markdown.doc"""
    multiplicative_jordan_decomposition(M::MatElem{T}) where T <:FieldElem -> MatElem{T}, MatElem{T}

Returns matrices $S$ and $U$ such that $U$ is unipotent, $S$ is semisimple and $M = SU$.
"""
function multiplicative_jordan_decomposition(M::MatElem{T}) where T <:FieldElem
  S, N = jordan_decomposition(M)
  U = inv(S)*N + identity_matrix(base_ring(M), ncols(M))
  return S, U
end

################################################################################
#
#  Algorithm from Steel - Requires factorization!
#
################################################################################

function split_primary(L::Dict, M::MatElem{T}) where T <: FieldElem
  for (g, S) in L
    newl = Vector{Tuple{MatElem{T}, MatElem{T}}}()
    for (W, w) in S
      if !iszero(w)
        push!(newl, (W, w))
        continue
      end
      #I need the sequence of the kernels of g^i(M) \cap W
      #So I restrict to W the endomorphism
      MW = restriction(M, W)
      #Now, MW is the restriction of M to W.
      #I compute g(MW) and then the kernels of the powers
      gMW = Hecke._subst(g, MW)
      kernels = Vector{typeof(M)}()
      push!(kernels, zero_matrix(base_ring(MW), 0, ncols(gMW)))
      d, K = kernel(gMW, side = :left)
      rref!(K)
      push!(kernels, K)
      M1 = gMW
      while true
        M1 *= gMW
        if iszero(M1)
          push!(kernels, identity_matrix(base_ring(MW), ncols(M1)))
          break
        end
        d1, K = kernel(M1, side = :left)
        rref!(K)
        push!(kernels, K)
      end
      #Now, I start taking vectors from the kernels till I generate the entire space
      subgen = zero_matrix(base_ring(M), 0, ncols(K))
      for i = length(kernels)-1:-1:1
        j = 1
        to_reduce = vcat(kernels[i], subgen)
        rref!(to_reduce)
        while j <= nrows(kernels[i+1])
          ej = sub(kernels[i+1], j:j, 1:ncols(K))
          res = Hecke.cleanvect(to_reduce, ej)
          while iszero(res) && j < nrows(kernels[i+1])
            j += 1
            ej = sub(kernels[i+1], j:j, 1:ncols(K))
            res = Hecke.cleanvect(to_reduce, ej)
          end
          if iszero(res)
            break
          end
          C = closure_with_pol(ej, MW)[1]
          push!(newl, (C*W, ej*W))
          subgen = vcat(subgen, C)
          rref!(subgen)
          if nrows(subgen) == nrows(W)
            i = 1
            break
          end
          to_reduce = vcat(kernels[i], subgen)
          rref!(to_reduce)
          j += 1
        end
      end
    end
   L[g] = newl
  end
  return L
end

function jordan_normal_form1(M::MatElem{T}) where T <: FieldElem
  K = base_ring(M)
  L = decompose_primary(M)
  L = split_primary(L, M)
  J = similar(M)
  B = similar(M)
  ind = 1
  for (g, v) in L
    d = degree(g)
    for i = 1:length(v)
      e = divexact(nrows(v[i][1]), d)
      J1 = jordan_block(g, e)
      w = v[i][2]
      for k = 1:e
        aux = w
        for j = 1:degree(g)
          aux = aux*M
          for s = 1:ncols(M)
            B[ind+(k-1)*d+j-1, s] = aux[1, s]
          end
        end
        w = w*Hecke._subst(g, M)
      end
      _copy_matrix_into_matrix(J, ind, ind, J1)
      ind += d*e
    end
  end
  return J, B
end

function rational_canonical_form1(M::MatElem{T}) where T <: FieldElem
  K = base_ring(M)
  L = decompose_primary(M)
  L = split_primary(L, M)
  #order the subspace by dimension
  max_length = 0
  for (q, v) in L
    v1 = sort(v, by = x -> nrows(x[1]))
    if length(v1) > max_length
      max_length = length(v1)
    end
    L[q] = v1
  end
  vectors_and_minpolys = Vector{Tuple{typeof(M), typeof(first(keys(L)))}}(undef, max_length)
  Kt = parent(first(keys(L)))
  for i = 1:max_length
    mp = one(Kt)
    w = zero_matrix(K, 1, ncols(M))
    for (q, v) in L
      if length(v) < i
        continue
      end
      W = v[length(v)-i+1][1]
      gw = v[length(v)-i+1][2]
      mult = divexact(nrows(W), degree(q))
      w += gw
      mp *= q^mult
    end
    vectors_and_minpolys[i] = (w, mp)
  end
  #Now, the matrix and the transformation matrix
  CF = similar(M)
  TM = similar(M)
  ind = 1
  for i = max_length:-1:1
    w, mp = vectors_and_minpolys[i]
    C = companion_matrix(mp)
    _copy_matrix_into_matrix(CF, ind, ind, C)
    S = closure_with_pol(w, M)[1]
    _copy_matrix_into_matrix(TM, ind, 1, S)
    ind += nrows(S)
  end
  return CF, TM
end

################################################################################
#
#  Simultaneous diagonalization
#
################################################################################

function commute_pairwise(L::Vector{S}) where S <: Hecke.MatElem{T} where T <:Hecke.FieldElem

  n = length(L)

  # Check pairwise commutativity
  for i in 1:n
    for j in i+1:n
      if L[i]*L[j] != L[j]*L[i]
        return false
      end
    end
  end
  return true
end

function issimultaneous_diagonalizable(L::Vector{S}) where S <: Hecke.MatElem{T} where T <:Hecke.FieldElem

  if !commute_pairwise(L)
    return false
  end

  for i = 1:length(L)
    f = minpoly(L[i])
    if !issquarefree(f)
      return false
    end
    lr = roots(f)
    if length(lr) != degree(f)
      return false
    end
  end
  return true
end

function _intersect(A::Hecke.MatElem{T}, B::Hecke.MatElem{T}) where T <: Hecke.FieldElem
  n = Hecke.nrows(A)
  M = Hecke.vcat(A, B)
  a, N = Hecke.kernel(M, side=:left)
  return N[1:a, 1:n]*A
end

function simultaneous_diagonalization(L...; check::Bool = true)
  return simultaneous_diagonalization(collect(L), check = check)
end

@doc Markdown.doc"""
    simultaneous_diagonalization(L::Vector{S}; check::Bool=false)

Returns a tuple whose first entry is the transformation matrix and whose
second entry is an array of matrices containing the diagonal forms of
the elements of $L$. If "check" is set to `true`, the algorithm checks whether
the matrices in $L$ are simultaneous diagonalizable before computing the transformation matrix. Default value is "true".
"""
function simultaneous_diagonalization(L::Vector{S}; check::Bool = true) where S <: MatElem{T} where T <: FieldElem

  if check
    if !issimultaneous_diagonalizable(L)
      error("Only one matrix as input!")
    end
  end

  n = length(L)
  if n == 1
    error("Only one matrix as input!")
  end
  s = Hecke.nrows(L[1])

  # Compute transformation marix
  egs = [eigenspaces(L[i]) for i = 1:length(L)]
  _egs = [Dict([x] => v for (x, v) in y) for y in egs]
  CE = common_eigenspaces(_egs)
  A =  Hecke.vcat(collect(values(CE)))

  # Compute diagonal forms
  D = [similar(L[1]) for i in 1:length(L)]
  m = 0
  for (v, k) in CE
    nr = Hecke.nrows(k)
    for j in 1:nr
      for i in 1:n
        D[i][m + j, m + j] = v[i]
      end
    end
    m = m + nr
  end
  return A, D
end

@doc Markdown.doc"""
    simultaneous_diagonalization(L::Vector{MatElem}, K::Field; check::Bool=false)

Returns a tuple whose first entry is the transformation matrix and whose
second entry is an array of matrices containing the diagonal forms of
the elements of $L$ computed over the field $K$. If "check" is set to `true`, the algorithm checks whether
the matrices in $L$ are simultaneous diagonalizable before computing the transformation matrix. Default value is "true".
"""
function simultaneous_diagonalization(L::Vector{S}, K::W; check::Bool = true) where S <: MatElem{T} where T <: FieldElem where W<:Field
  L1 = [change_base_ring(x, K) for x in L]
  return simultaneous_diagonalization(L1, check = check)
end

function common_eigenspaces(L::Vector{Dict{Vector{T}, S}}) where S<:Hecke.MatElem{T} where T<:Hecke.FieldElem

  n = length(L)
  if n==1
    return L[1]
  end
  k = BigInt(floor(n/2))
  return intersect_eigenspaces(common_eigenspaces(L[1:k]), common_eigenspaces(L[k+1:n]))
end

function intersect_eigenspaces(L1::Dict{Vector{T}, S}, L2::Dict{Vector{T}, S}) where S<:Hecke.MatElem{T} where T <: Hecke.FieldElem

  L = Dict{keytype(L1), valtype(L1)}()
  for (k1, v1) in L1
    for (k2, v2) in L2
      I = intersect_spaces(v1, v2)
      if !iszero(I)
        push!(L, vcat(k1, k2)  => I)
      end
    end
  end
  return L
end


text:
# Unit tests for functions in gates.jl

@testset "Tests: elementary universal gates in gates.jl" begin
    @test isapprox(QCO.U2Gate(0,-π/4), QCO.U3Gate(π/2,0,-π/4), atol=tol_0)
    @test isapprox(QCO.U1Gate(-π/4), QCO.U3Gate(0,0,-π/4), atol=tol_0)

    ctrl_qbit_1 = Array{Complex{Float64},2}([1 0; 0 0])
    ctrl_qbit_2 = Array{Complex{Float64},2}([0 0; 0 1])

    @test isapprox(kron(ctrl_qbit_1, QCO.IGate(1)) + kron(ctrl_qbit_2, QCO.XGate()), QCO.CXGate(), atol = tol_0)
    @test isapprox(kron(ctrl_qbit_1, QCO.IGate(1)) + kron(ctrl_qbit_2, QCO.YGate()), QCO.CYGate(), atol = tol_0)
    @test isapprox(kron(ctrl_qbit_1, QCO.IGate(1)) + kron(ctrl_qbit_2, QCO.ZGate()), QCO.CZGate(), atol = tol_0)
    @test isapprox(kron(ctrl_qbit_1, QCO.IGate(1)) + kron(ctrl_qbit_2, QCO.HGate()), QCO.CHGate(), atol = tol_0)
    @test isapprox(QCO.CXGate(), QCO.CNotGate(), atol = tol_0)

    ϴ1 = π/3
    ϴ2 = -2*π/3
    @test isapprox(kron(ctrl_qbit_1, QCO.IGate(1)) + kron(ctrl_qbit_2, QCO.RXGate(ϴ1)), QCO.CRXGate(ϴ1), atol = tol_0)
    @test isapprox(kron(ctrl_qbit_1, QCO.IGate(1)) + kron(ctrl_qbit_2, QCO.RXGate(ϴ2)), QCO.CRXGate(ϴ2), atol = tol_0)
    @test isapprox(kron(ctrl_qbit_1, QCO.IGate(1)) + kron(ctrl_qbit_2, QCO.RYGate(ϴ1)), QCO.CRYGate(ϴ1), atol = tol_0)
    @test isapprox(kron(ctrl_qbit_1, QCO.IGate(1)) + kron(ctrl_qbit_2, QCO.RYGate(ϴ2)), QCO.CRYGate(ϴ2), atol = tol_0)
    @test isapprox(kron(ctrl_qbit_1, QCO.IGate(1)) + kron(ctrl_qbit_2, QCO.RZGate(ϴ1)), QCO.CRZGate(ϴ1), atol = tol_0)
    @test isapprox(kron(ctrl_qbit_1, QCO.IGate(1)) + kron(ctrl_qbit_2, QCO.RZGate(ϴ2)), QCO.CRZGate(ϴ2), atol = tol_0)

    θ = π/3
    ϕ = -2*π/3
    λ = π/6
    @test isapprox(kron(ctrl_qbit_1, QCO.IGate(1)) + kron(ctrl_qbit_2, QCO.U3Gate(θ,ϕ,λ)), QCO.CU3Gate(θ,ϕ,λ), atol = tol_0)

    @test isapprox(QCO.SwapGate(), QCO.CNotGate() * QCO.CNotRevGate() * QCO.CNotGate(), atol = tol_0)

    @test isapprox(QCO.iSwapGate(), QCO.CNotRevGate() * QCO.CNotGate() * kron(QCO.SGate(), QCO.IGate(1)) * QCO.CNotRevGate(), atol = tol_0)

    @test isapprox(QCO.PhaseGate(λ), QCO.U3Gate(0,0,λ), atol = tol_0)

    @test isapprox(QCO.DCXGate(), QCO.CNotGate() * QCO.CNotRevGate(), atol = tol_0)

    S1_S2 = kron(QCO.SGate(), QCO.SGate())
    H2 = kron(QCO.IGate(1), QCO.HGate())
    @test isapprox(QCO.MGate(), QCO.CNotRevGate() * H2 * S1_S2)

    Z1 = QCO.get_full_sized_gate("Z_1", 2);
    Z2 = QCO.get_full_sized_gate("Z_2", 2);
    T2 = QCO.get_full_sized_gate("T_2", 2);
    Y_2 = QCO.get_full_sized_gate("Y_2", 2);
    cnot_12 = QCO.get_full_sized_gate("CNot_1_2", 2);
    cnot_21 = QCO.get_full_sized_gate("CNot_2_1", 2);
    Sdagger1 = QCO.get_full_sized_gate("Sdagger_1", 2);
    Tdagger1 = QCO.get_full_sized_gate("Tdagger_1", 2);
    SX1 = QCO.get_full_sized_gate("SX_1", 2);
    SXdagger2 = QCO.get_full_sized_gate("SXdagger_2", 2);
    @test isapprox(-QCO.HCoinGate(), Z2 * cnot_21 * SXdagger2 * Y_2 * cnot_21 * Tdagger1 * Z1 * T2 * Sdagger1 * cnot_12 * Sdagger1 * SX1 * cnot_21 * cnot_12, atol = tol_0)

    CV_23 = QCO.get_full_sized_gate("CV_2_3", 3);
    cnot_12 = QCO.get_full_sized_gate("CNot_1_2", 3);
    CVdagger_23 = QCO.get_full_sized_gate("CVdagger_2_3", 3);
    CV_13 = QCO.get_full_sized_gate("CV_1_3", 3);
    @test isapprox(QCO.ToffoliGate(), CV_23 * cnot_12 * CVdagger_23 * cnot_12 * CV_13, atol = tol_0)

    CVdagger_12 = QCO.get_full_sized_gate("CVdagger_1_2", 3);
    CVdagger_31 = QCO.get_full_sized_gate("CVdagger_3_1", 3);
    @test isapprox(QCO.IGate(3), CVdagger_12 * CVdagger_12 * CVdagger_12 * CVdagger_12 * CVdagger_31 * CVdagger_31 * CVdagger_31 * CVdagger_31, atol = tol_0)

    H_2 = QCO.get_full_sized_gate("H_2", 2);
    T_1 = QCO.get_full_sized_gate("T_1", 2);
    Tdagger2 = QCO.get_full_sized_gate("Tdagger_2", 2);
    @test isapprox(QCO.CVdaggerGate(), H_2 * Tdagger1 * cnot_21 * T_1 * Tdagger2 * cnot_21 * H_2, atol = tol_0)

    cnot_12 = QCO.get_full_sized_gate("CNot_1_2", 2);
    @test isapprox(QCO.CVGate(), cnot_12 * QCO.CVdaggerGate(), atol = tol_0)

    # Next 3 tests from: https://american-cse.org/csci2015/data/9795a059.pdf

    CV_31 = QCO.get_full_sized_gate("CV_3_1", 3);
    CNot_23 = QCO.get_full_sized_gate("CNot_2_3", 3);
    CVdagger_31 = QCO.get_full_sized_gate("CVdagger_3_1", 3);
    CV_21 = QCO.get_full_sized_gate("CV_2_1", 3);

    CVdagger_21 = QCO.get_full_sized_gate("CVdagger_2_1", 3);

    @test isapprox(CV_31 * CNot_23 * CVdagger_31 * CNot_23 * CV_21, CVdagger_21 * CNot_23 * CV_31 * CNot_23 * CVdagger_31, atol = tol_0)

    CV_12 = QCO.get_full_sized_gate("CV_1_2", 3);
    CV_32 = QCO.get_full_sized_gate("CV_3_2", 3);
    CVdagger_32 = QCO.get_full_sized_gate("CVdagger_3_2", 3);
    CVdagger_12 = QCO.get_full_sized_gate("CVdagger_1_2", 3);
    CNot_13 = QCO.get_full_sized_gate("CNot_1_3", 3);

    @test isapprox(CV_32 * CNot_13 * CVdagger_32 * CNot_13 * CV_12, CVdagger_32 * CNot_13 * CV_32 * CNot_13 * CVdagger_12, atol = tol_0)

    CVdagger_13 = QCO.get_full_sized_gate("CVdagger_1_3", 3);
    cnot_21 = QCO.get_full_sized_gate("CNot_2_1", 3);
    @test isapprox(CV_13 * cnot_21 * CVdagger_13 * cnot_21 * CV_23, CVdagger_13 * cnot_21 * CV_13 * cnot_21 * CVdagger_23, atol = tol_0)

    # 2-qubit Grover's diffusion operator (Ref: https://arxiv.org/pdf/1804.03719.pdf)
    H_1 = QCO.get_full_sized_gate("H_1", 2);
    X_2 = QCO.get_full_sized_gate("X_2", 2);
    CNot_1_2 = QCO.get_full_sized_gate("CNot_1_2", 2);

    @test isapprox(QCO.GroverDiffusionGate(), Y_2 * H_1 * X_2 * CNot_1_2 * H_1 * Y_2, atol=tol_0)

    # 2 Qubit QFT (Ref: https://www.cs.bham.ac.uk/internal/courses/intro-mqc/current/lecture06_handout.pdf)

    SWAP = QCO.get_full_sized_gate("Swap_1_2", 2);
    CU = QCO.get_full_sized_gate("CU3_2_1", 2, angle = [0, π/4, π/4]);

    @test isapprox(QCO.QFT2Gate(), H_1 * CU * H_2 * SWAP)

    # CRY Decomp (Ref: https://quantumcomputing.stackexchange.com/questions/2143/how-can-a-controlled-ry-be-made-from-cnots-and-rotations)

    CNOT12 = QCO.get_full_sized_gate("CNot_1_2", 2);
    U3_negpi = QCO.get_full_sized_gate("U3_2", 2, angle=[-pi/2, 0, 0]);
    U3_pi = QCO.get_full_sized_gate("U3_2", 2, angle=[pi/2, 0, 0]);

    @test isapprox(QCO.CRYGate(pi), CNOT12 * U3_negpi * CNOT12 * U3_pi, atol = tol_0)

    # Identity tests

    Id = QCO.IGate(2);
    CRXRev = QCO.get_full_sized_gate("CRX_2_1", 2, angle=pi);

    @test isapprox(Id, CRXRev * CRXRev * CRXRev * CRXRev, atol = tol_0)

    CRYRev = QCO.get_full_sized_gate("CRY_2_1", 2, angle=pi);
    @test isapprox(Id, CRYRev * CRYRev * CRYRev * CRYRev, atol = tol_0)

    CRZRev = QCO.get_full_sized_gate("CRZ_2_1", 2, angle=pi);
    @test isapprox(Id, CRZRev * CRZRev * CRZRev * CRZRev, atol = tol_0)

    # Rev gate tests for involution
    @test isapprox(QCO.CXRevGate(), QCO.CNotRevGate(), atol = tol_0)
    @test isapprox(QCO.CXGate() * QCO.SwapGate() * QCO.CXRevGate() * QCO.SwapGate(), QCO.IGate(2), atol = tol_0)
    @test isapprox(QCO.CYGate() * QCO.SwapGate() * QCO.CYRevGate() * QCO.SwapGate(), QCO.IGate(2), atol = tol_0)
    @test isapprox(QCO.CZGate() * QCO.SwapGate() * QCO.CZRevGate() * QCO.SwapGate(), QCO.IGate(2), atol = tol_0)
    @test isapprox(QCO.CHGate() * QCO.SwapGate() * QCO.CHRevGate() * QCO.SwapGate(), QCO.IGate(2), atol = tol_0)

    # Rev gate tests for non-involution
    @test isapprox(QCO.CVGate() * QCO.SwapGate() * QCO.CVRevGate() * QCO.SwapGate(), QCO.CVGate()^2, atol = tol_0)
    @test isapprox(QCO.CSXGate() * QCO.SwapGate() * QCO.CSXRevGate() * QCO.SwapGate(), QCO.CSXGate()^2, atol = tol_0)

    @test isapprox(QCO.WGate(), QCO.HCoinGate(), atol = tol_0)

    # Fredkin test = CSwapGate
    CNot_32 = QCO.get_full_sized_gate("CNot_3_2", 3)
    CV_23 = QCO.get_full_sized_gate("CV_2_3", 3)
    CV_13 = QCO.get_full_sized_gate("CV_1_3", 3)
    CNot_12 = QCO.get_full_sized_gate("CNot_1_2", 3)
    CVdagger_23 = QCO.get_full_sized_gate("CVdagger_2_3", 3)
    @test isapprox(QCO.CSwapGate(), CNot_32 * CV_23 * CV_13 * CNot_12 * CVdagger_23 * CNot_12 * CNot_32, atol = tol_0)

    # CCZGate test: CCZ is equivalent to Toffoli when the target qubit is conjugated by Hadamard gates
    H_3 = QCO.get_full_sized_gate("H_3", 3)
    @test isapprox(QCO.ToffoliGate(), H_3 * QCO.CCZGate() * H_3, atol = tol_0)

    # Peres test
    CVdagger_13 = QCO.get_full_sized_gate("CVdagger_1_3", 3)
    @test isapprox(QCO.PeresGate(), QCO.ToffoliGate() * CNot_12, atol = tol_0)
    @test isapprox(QCO.PeresGate(), CVdagger_13 * CVdagger_23 * CNot_12 * CV_23, atol = tol_0)

    # iSwap test
    @test isapprox(QCO.get_full_sized_gate("iSwap_2_1", 3), QCO.get_full_sized_gate("iSwap_1_2", 3), atol = tol_0)
    @test isapprox(QCO.get_full_sized_gate("iSwap_2_3", 4), QCO.get_full_sized_gate("iSwap_3_2", 4), atol = tol_0)
    @test isapprox(QCO.get_full_sized_gate("iSwap_3_5", 5), QCO.get_full_sized_gate("iSwap_5_3", 5), atol = tol_0)

    # Sycamore test
    @test isapprox(QCO.SycamoreGate()^12, QCO.IGate(2), atol = tol_0)

    # RCCX test
    T_3 = QCO.get_full_sized_gate("T_3", 3)
    Tdagger_3 = QCO.get_full_sized_gate("Tdagger_3", 3)
    @test isapprox(QCO.RCCXGate(),  H_3 * Tdagger_3 * CNot_23 * T_3 * CNot_13 * Tdagger_3 * CNot_23 * T_3 * H_3, atol = tol_0)

    # Margolus test
    RY_a = QCO.get_full_sized_gate("RY_3", 3, angle = -π/4)
    RY_b = QCO.get_full_sized_gate("RY_3", 3, angle = π/4)
    @test isapprox(QCO.MargolusGate(),  RY_a * CNot_23 * RY_a * CNot_13 * RY_b * CNot_23 * RY_b, atol = tol_0)

    # CiSwap test
    @test isapprox(QCO.CiSwapGate(),  CNot_32 * H_3 * CNot_23 * T_3 * CNot_13 * Tdagger_3 * CNot_23 * T_3 * CNot_13 * Tdagger_3 * H_3 * CNot_32, atol = tol_0)

end


text:
#=
assume two populations and assign observations to the populations and simultaneously fit
polynomial functions to these populations; on a given domain containing the observations,
the functions must satisfy certain shape constraints

given domain D ⊆ ℝⁿ, observations (Xᵢ ∈ D, Yᵢ ∈ ℝ₊), solve for p₁(x), p₂(x):
min  ‖z‖ :                              residual loss
p₂(x) ≥ 0, ∀x ∈ D                       p2 nonnegative
p₁(x) ≥ p₂(x), ∀x ∈ D                   p1 above p2
1 ≥ p₁(x), ∀x ∈ D                       p1 below 1
bᵢ ∈ {0,1}, ∀i                          assignment of points
0.4m <= Σᵢbᵢ <= 0.6m                    approximately balance assignments
|zᵢ - (Yᵢ - p₁(Xᵢ))| ≤ Mᵢ(1 - bᵢ), ∀i   residuals for p1
|zᵢ - (Yᵢ - p₂(Xᵢ))| ≤ Mᵢbᵢ, ∀i         residuals for p2
=#

struct PolyRegression <: ExampleInstance
    n::Int # dimension of independent variables / polynomial domain
    halfdeg::Int # polynomial regressor half-degree
    m::Int # number of observations
    signal_ratio::Float64 # signal-to-noise ratio (zero if no noise)
    use_nat::Bool # use WSOS cone formulation, else SDP formulation
end

function build(inst::PolyRegression)
    (n, m, signal_ratio) = (inst.n, inst.m, inst.signal_ratio)
    @assert m > n >= 1

    # setup interpolation
    D = PolyUtils.BoxDomain{Float64}(zeros(n), ones(n))
    (U, _, Ps, V) = PolyUtils.interpolate(D, inst.halfdeg, calc_V = true)
    @assert m >= 2U
    F = qr!(Array(V'), ColumnNorm())

    # generate noisy data in D = [0, 1]ⁿ from two underlying functions
    X = rand(m, n)
    V_X = PolyUtils.make_chebyshev_vandermonde(X, 2 * inst.halfdeg)
    p_X = F \ V_X'

    fs = make_nonneg_polys(2, Ps)
    f1 = fs[1]
    f2 = fs[2]
    @assert minimum(p_X' * f1) >= -1e-4
    @assert minimum(p_X' * f2) >= -1e-4

    f1 .+= f2
    f1_max = -find_poly_min(-f1, Ps)
    @assert f1_max >= -1e-4
    f1 ./= f1_max
    f2 ./= f1_max

    m1 = div(m, 2)
    Y1 = p_X' * f1
    Y2 = p_X' * f2
    @assert minimum(Y1) >= -1e-4
    @assert maximum(Y1) <= 1 + 1e-4
    @assert minimum(Y2) >= -1e-4
    @assert maximum(Y2) <= 1 + 1e-4
    Y = vcat(Y1[1:m1], Y2[(m1 + 1):m])

    if !iszero(signal_ratio)
        noise = randn(m)
        noise .*= norm(Y) / sqrt(signal_ratio) / norm(noise)
        Y .+= noise
        Y = max.(Y, 0.0)
        Y = min.(Y, 1.0)
    end
    @assert all(0 .<= Y .<= 1)

    # build model
    model = JuMP.Model()
    JuMP.@variable(model, p1[1:U])
    JuMP.@variable(model, p2[1:U])
    JuMP.@variable(model, r1[1:m])
    JuMP.@variable(model, r2[1:m])
    JuMP.@variable(model, z[1:m])
    JuMP.@variable(model, b[1:m], Bin)

    JuMP.@variable(model, epi)
    JuMP.@objective(model, Min, epi)

    # TODO tighten big-Ms, inefficient sparsity
    JuMP.@constraints(model, begin
        vcat(epi, z) in Hypatia.EpiNormEuclCone{Float64}(1 + m)
        r1 .== z - Y + p_X' * p1
        r2 .== z - Y + p_X' * p2
        ceil(0.4 * m) <= sum(b)
        sum(b) <= floor(0.6 * m)
        -(1 .- b) .<= r1
        r1 .<= 1 .- b
        -b .<= r2
        r2 .<= b
    end)

    # shape constraints
    add_wsos(aff) = (inst.use_nat ? add_wsos_nat : add_wsos_ext)(Ps, aff, model)
    add_wsos(p2)
    add_wsos(p1 - p2)
    add_wsos(1 .- p1)

    # save for use in tests
    model.ext[:b] = b
    model.ext[:p1] = p1
    model.ext[:p2] = p2
    model.ext[:Ps] = Ps

    return model
end

function test_extra(inst::PolyRegression, model::JuMP.Model)
    stat = JuMP.termination_status(model)
    @test stat in OPT_OR_LIMIT
    (stat == MOI.OPTIMAL) || return

    tol = eps()^0.2
    b = JuMP.value.(model.ext[:b])
    @test b ≈ round.(Int, b) atol = tol rtol = tol
    @test all(-tol .<= b .<= 1 + tol)
    # check WSOS feasibility
    Ps = model.ext[:Ps]
    p1 = JuMP.value.(model.ext[:p1])
    p2 = JuMP.value.(model.ext[:p2])
    @test check_nonneg(p2 .+ tol, Ps)
    @test check_nonneg(p1 - p2 .+ tol, Ps)
    @test check_nonneg((1 + tol) .- p1, Ps)
    # objective is 0 if no noise
    if iszero(inst.signal_ratio)
        @test JuMP.objective_value(model) ≈ 0 atol = tol rtol = tol
    end
    return
end


text:
module Hello

using UnicodePlots

Base.@ccallable function julia_main(ARGS::Vector{String})::Cint
    println("hello, world")
    @show sin(0.0)
    println(lineplot(1:100, sin.(range(0, stop=2π, length=100))))
    return 0
end

end


text:
using ITensors, Test

@testset "basic DMRG" begin
  N = 100
  sites = spinOneSites(N)

  ampo = AutoMPO(sites)
  for j=1:N-1
    add!(ampo,"Sz",j,"Sz",j+1)
    add!(ampo,0.5,"S+",j,"S-",j+1)
    add!(ampo,0.5,"S-",j,"S+",j+1)
  end
  H = toMPO(ampo)

  psi = randomMPS(sites)

  sweeps = Sweeps(5)
  @test length(sweeps) == 5
  maxdim!(sweeps,10,20,100,100)
  mindim!(sweeps,1,10,20,20)
  cutoff!(sweeps,1E-11)
  energy,psi = dmrg(H,psi,sweeps,maxiter=2,quiet=true)
  @test energy ≈ -138.94 rtol=1e-3
  # test with SVD too!
  psi = randomMPS(sites)
  energy,psi = dmrg(H,psi,sweeps,maxiter=2,quiet=true,which_factorization="svd")
  @test energy ≈ -138.94 rtol=1e-3
end


text:
# waiting for official Timestamp support in Julia. Calendar.jl is too slow to load.
const date_pattern = anchored_regex("(\\d{4})-(\\d{2})-(\\d{2})T(\\d{2}):(\\d{2}):(\\d{2})Z")

immutable DateTime
    year::Int
    month::Int
    date::Int
    hour::Int
    minute::Int
    second::Int
end

function ymd_hms(args...)
    DateTime(args[1:6]...)
end


text:
module ReinforcementLearningDatasets

const RLDatasets = ReinforcementLearningDatasets
export RLDatasets

using DataDeps

include("deep_ope/d4rl/d4rl_policies.jl")

include("d4rl/d4rl/register.jl")
include("d4rl/d4rl_pybullet/register.jl")
include("atari/register.jl")
include("rl_unplugged/atari/register.jl")
include("rl_unplugged/bsuite/register.jl")
include("rl_unplugged/dm/register.jl")
include("deep_ope/d4rl/register.jl")

include("common.jl")
include("init.jl")

include("d4rl/d4rl_dataset.jl")
include("atari/atari_dataset.jl")
include("rl_unplugged/util.jl")
include("rl_unplugged/atari/rl_unplugged_atari.jl")
include("rl_unplugged/bsuite/bsuite.jl")
include("rl_unplugged/dm/rl_unplugged_dm.jl")
include("deep_ope/d4rl/d4rl_policy.jl")

include("deep_ope/d4rl/evaluate.jl")

end