text:
__precompile__()

module JLD
using HDF5, FileIO, Compat

import HDF5: close, dump, exists, file, getindex, setindex!, g_create, g_open, o_delete, name, names, read, write,
             HDF5ReferenceObj, HDF5BitsKind, ismmappable, readmmap
import Base: convert, length, endof, show, done, next, ndims, start, delete!, eltype,
             size, sizeof, unsafe_convert, datatype_pointerfree
import LegacyStrings: UTF16String

@noinline gcuse(x) = x # because of use of `pointer`, need to mark gc-use end explicitly

const magic_base = "Julia data file (HDF5), version "
const version_current = v"0.1.1"
const pathrefs = "/_refs"
const pathtypes = "/_types"
const pathrequire = "/_require"
const pathcreator = "/_creator"
const name_type_attr = "julia type"

const BitsKindOrString = Union{HDF5BitsKind, String}

function julia_type(s::AbstractString)
    s = replace(s, r"ASCIIString|UTF8String|ByteString", "String")
    s = replace(s, "Base.UTF16String", "LegacyStrings.UTF16String")
    _julia_type(s)
end

### Dummy types used for converting attribute strings to Julia types
mutable struct UnsupportedType; end
mutable struct UnconvertedType; end

struct JldDatatype
    dtype::HDF5Datatype
    index::Int
end
sizeof(T::JldDatatype) = sizeof(T.dtype)

struct JldWriteSession
    persist::Vector{Any} # To hold objects that should not be garbage-collected
    h5ref::ObjectIdDict  # To hold mapping from Object/Array -> HDF5ReferenceObject

    JldWriteSession() = new(Any[], ObjectIdDict())
end

# The Julia Data file type
# Purpose of the nrefs field:
# length(group) only returns the number of _completed_ items in a group. Since
# we'll write recursively, we need to keep track of the number of reference
# objects _started_.
mutable struct JldFile <: HDF5.DataFile
    plain::HDF5File
    version::VersionNumber
    toclose::Bool
    writeheader::Bool
    mmaparrays::Bool
    compatible::Bool
    compress::Bool
    h5jltype::Dict{Int,Type}
    jlh5type::Dict{Type,JldDatatype}
    jlref::Dict{HDF5ReferenceObj,WeakRef}
    truncatemodules::Vector{String}
    gref # Group references; can't annotate type here due to circularity
    nrefs::Int

    function JldFile(plain::HDF5File, version::VersionNumber=version_current, toclose::Bool=true,
                     writeheader::Bool=false, mmaparrays::Bool=false,
                     compatible::Bool=false, compress::Bool=false)
        f = new(plain, version, toclose, writeheader, mmaparrays, compatible,
                compress & !mmaparrays,
                Dict{HDF5Datatype,Type}(), Dict{Type,HDF5Datatype}(),
                Dict{HDF5ReferenceObj,WeakRef}(), String[])
        if toclose
            finalizer(f, close)
        end
        f
    end
end

struct JldGroup
    plain::HDF5Group
    file::JldFile
end

struct JldDataset
    plain::HDF5Dataset
    file::JldFile
end

iscompatible(f::JldFile) = f.compatible
iscompatible(g::JldGroup) = g.file.compatible

iscompressed(f::JldFile) = f.compress
iscompressed(g::JldGroup) = g.file.compress

ismmapped(f::JldFile) = f.mmaparrays
ismmapped(d::JldGroup) = d.file.mmaparrays

struct PointerException <: Exception; end
show(io::IO, ::PointerException) = print(io, "cannot write a pointer to JLD file")

struct TypeMismatchException <: Exception
    typename::String
end
show(io::IO, e::TypeMismatchException) =
    print(io, "stored type $(e.typename) does not match currently loaded type")

include("jld_types.jl")

file(x::JldFile) = x
file(x::Union{JldGroup, JldDataset}) = x.file

function close(f::JldFile)
    if f.toclose
        # Close types
        for x in values(f.jlh5type)
            close(x.dtype)
        end

        # Close reference group
        isdefined(f, :gref) && close(f.gref)

        # Ensure that all other datasets, groups, and datatypes are closed (ref #176)
        for obj_id in HDF5.h5f_get_obj_ids(f.plain.id, HDF5.H5F_OBJ_DATASET | HDF5.H5F_OBJ_GROUP | HDF5.H5F_OBJ_DATATYPE)
            HDF5.h5o_close(obj_id)
        end

        # Close file
        close(f.plain)
        if f.writeheader
            magic = zeros(UInt8, 512)
            tmp = string(magic_base, f.version)
            magic[1:length(tmp)] = Vector{UInt8}(tmp)
            rawfid = open(f.plain.filename, "r+")
            write(rawfid, magic)
            close(rawfid)
        end
        f.toclose = false
    end
    nothing
end
close(g::Union{JldGroup, JldDataset}) = close(g.plain)
show(io::IO, fid::JldFile) = isvalid(fid.plain) ? print(io, "Julia data file version ", fid.version, ": ", fid.plain.filename) : print(io, "Julia data file (closed): ", fid.plain.filename)

function jldopen(filename::AbstractString, rd::Bool, wr::Bool, cr::Bool, tr::Bool, ff::Bool; mmaparrays::Bool=false, compatible::Bool=false, compress::Bool=false)
    local fj
    if ff && !wr
        error("Cannot append to a write-only file")
    end
    if !cr && !isfile(filename)
        error("File ", filename, " cannot be found")
    end
    version = version_current
    pa = p_create(HDF5.H5P_FILE_ACCESS)
    # HDF5.h5p_set_libver_bounds(pa, HDF5.H5F_LIBVER_18, HDF5.H5F_LIBVER_18)
    try
        pa["fclose_degree"] = HDF5.H5F_CLOSE_STRONG
        if cr && (tr || !isfile(filename))
            # We're truncating, so we don't have to check the format of an existing file
            # Set the user block to 512 bytes, to save room for the header
            p = p_create(HDF5.H5P_FILE_CREATE)
            local f
            try
                p["userblock"] = 512
                f = HDF5.h5f_create(filename, HDF5.H5F_ACC_TRUNC, p.id, pa.id)
            finally
                close(p)
            end
            fj = JldFile(HDF5File(f, filename, false), version, true, true, mmaparrays, compatible, compress)
            # Record creator information. Don't use any fancy types here,
            # because we want to be able to read this even when formats change.
            write(fj, joinpath(pathcreator, "JULIA_MAJOR"), VERSION.major)
            write(fj, joinpath(pathcreator, "JULIA_MINOR"), VERSION.minor)
            write(fj, joinpath(pathcreator, "JULIA_PATCH"), VERSION.patch)
            if !isempty(VERSION.prerelease)
                write(fj, joinpath(pathcreator, "JULIA_PRERELEASE"), [VERSION.prerelease...])
            end
            if !isempty(VERSION.build)
                write(fj, joinpath(pathcreator, "JULIA_BUILD"), [VERSION.build...])
            end
            write(fj, joinpath(pathcreator, "WORD_SIZE"), Sys.WORD_SIZE)
            write(fj, joinpath(pathcreator, "ENDIAN_BOM"), ENDIAN_BOM)
        else
            # Test whether this is a jld file
            sz = filesize(filename)
            if sz < 512
                error("File size indicates $filename cannot be a Julia data file")
            end
            magic = Vector{UInt8}(512)
            rawfid = open(filename, "r")
            try
                magic = read!(rawfid, magic)
            finally
                close(rawfid)
            end
            if startswith(magic, Vector{UInt8}(magic_base))
                version = convert(VersionNumber, unsafe_string(pointer(magic) + length(magic_base)))
                gcuse(magic)
                if version < v"0.1.0"
                    fj = JLD00.jldopen(filename, rd, wr, cr, tr, ff; mmaparrays=mmaparrays)
                else
                    f = HDF5.h5f_open(filename, wr ? HDF5.H5F_ACC_RDWR : HDF5.H5F_ACC_RDONLY, pa.id)
                    fj = JldFile(HDF5File(f, filename, false), version, true, cr|wr, mmaparrays, compatible, compress)
                    # Load any required files/packages
                    if exists(fj, pathrequire)
                        r = read(fj, pathrequire)
                        for fn in r
                            mod = path2modsym(fn)
                            eval(Expr(:import, mod))
                        end
                    end
                end
            else
                if ishdf5(filename)
                    println("$filename is an HDF5 file, but it is not a recognized Julia data file. Opening anyway.")
                    fj = JldFile(h5open(filename, rd, wr, cr, tr, ff), version_current, true, false, mmaparrays, compatible, compress)
                else
                    error("$filename does not seem to be a Julia data or HDF5 file")
                end
            end
        end
    finally
        close(pa)
    end
    return fj
end

function jldopen(fname::AbstractString, mode::AbstractString="r"; mmaparrays::Bool=false, compatible::Bool=false, compress::Bool=false)
    mode == "r"  ? jldopen(fname, true , false, false, false, false, mmaparrays=mmaparrays, compatible=compatible, compress=compress) :
    mode == "r+" ? jldopen(fname, true , true , false, false, false, mmaparrays=mmaparrays, compatible=compatible, compress=compress) :
    mode == "w"  ? jldopen(fname, false, true , true , true , false, mmaparrays=mmaparrays, compatible=compatible, compress=compress) :
#     mode == "w+" ? jldopen(fname, true , true , true , true , false) :
#     mode == "a"  ? jldopen(fname, false, true , true , false, true ) :
#     mode == "a+" ? jldopen(fname, true , true , true , false, true ) :
    error("invalid open mode: ", mode)
end

function jldopen(f::Function, args...; kws...)
    jld = jldopen(args...; kws...)
    try
        f(jld)
    finally
        close(jld)
    end
end

function creator(file::JldFile, key::AbstractString)
    if file.version < v"0.1.1"
        return nothing
    end
    if key == "VERSION"
        path_prerelease = joinpath(pathcreator, "JULIA_PRERELEASE")
        prerelease = exists(file, path_prerelease) ? tuple(read(file, path_prerelease)...) : ()
        path_build = joinpath(pathcreator, "JULIA_BUILD")
        build = exists(file, path_build) ? tuple(read(file, path_build)...) : ()
        return VersionNumber(read(file, joinpath(pathcreator, "JULIA_MAJOR")),
                             read(file, joinpath(pathcreator, "JULIA_MINOR")),
                             read(file, joinpath(pathcreator, "JULIA_PATCH")),
                             prerelease,
                             build)
    elseif key in ("WORD_SIZE", "ENDIAN_BOM")
        return read(file, joinpath(pathcreator, key))
    else
        error("$key not recognized")
    end
end

function jldobject(obj_id::HDF5.Hid, parent)
    obj_type = HDF5.h5i_get_type(obj_id)
    obj_type == HDF5.H5I_GROUP ? JldGroup(HDF5Group(obj_id, file(parent.plain)), file(parent)) :
    obj_type == HDF5.H5I_DATATYPE ? HDF5Datatype(obj_id) :
    obj_type == HDF5.H5I_DATASET ? JldDataset(HDF5Dataset(obj_id, file(parent.plain)), file(parent)) :
    error("Invalid object type for path ", path)
end

getindex(parent::Union{JldFile, JldGroup}, path::String) =
    jldobject(HDF5.h5o_open(parent.plain.id, path), parent)

function getindex(parent::Union{JldFile, JldGroup, JldDataset}, r::HDF5ReferenceObj)
    if r == HDF5.HDF5ReferenceObj_NULL; error("Reference is null"); end
    obj_id = HDF5.h5r_dereference(parent.plain.id, HDF5.H5R_OBJECT, r)
    jldobject(obj_id, parent)
end

### "Inherited" behaviors
g_create(parent::Union{JldFile, JldGroup}, args...) = JldGroup(g_create(parent.plain, args...), file(parent))
function g_create(f::Function, parent::Union{JldFile, JldGroup}, args...)
    g = JldGroup(g_create(parent.plain, args...), file(parent))
    try
        f(g)
    finally
        close(g)
    end
end
g_open(parent::Union{JldFile, JldGroup}, args...) = JldGroup(g_open(parent.plain, args...), file(parent))
name(p::Union{JldFile, JldGroup, JldDataset}) = name(p.plain)
eltype(p::JldDataset) = eltype(p.plain)
exists(p::Union{JldFile, JldGroup, JldDataset}, path::String) = exists(p.plain, path)
root(p::Union{JldFile, JldGroup, JldDataset}) = g_open(file(p), "/")
o_delete(parent::Union{JldFile, JldGroup}, args...) = o_delete(parent.plain, args...)
function ensurepathsafe(path::String)
    if any([startswith(path, s) for s in (pathrefs,pathtypes,pathrequire)])
        error("$name is internal to the JLD format, use o_delete if you really want to delete it")
    end
end
function delete!(o::JldDataset)
    fullpath = name(o)
    ensurepathsafe(fullpath)
    o_delete(o.file, fullpath)
    refspath = joinpath(pathrefs, fullpath[2:end])
    exists(o.file, refspath) && o_delete(o.file, refspath)
end
function delete!(g::JldGroup)
    fullpath = name(g)
    ensurepathsafe(fullpath)
    for o in g typeof(o) == JldDataset && delete!(o) end
    o_delete(g.file,name(g))
end
function delete!(parent::Union{JldFile, JldGroup}, path::String)
    exists(parent, path) || error("$path does not exist in $parent")
    delete!(parent[path])
end
delete!(parent::Union{JldFile, JldGroup}, args::Tuple{Vararg{String}}) = for a in args delete!(parent,a) end
ismmappable(obj::JldDataset) = ismmappable(obj.plain)
readmmap(obj::JldDataset, args...) = readmmap(obj.plain, args...)
setindex!(parent::Union{JldFile, JldGroup}, val, path::String) = write(parent, path, val)

start(parent::Union{JldFile, JldGroup}) = (names(parent), 1)
done(parent::Union{JldFile, JldGroup}, state) = state[2] > length(state[1])
next(parent::Union{JldFile, JldGroup}, state) = parent[state[1][state[2]]], (state[1], state[2]+1)

### Julia data file format implementation ###

### Read ###

function read(parent::Union{JldFile, JldGroup}, name::String)
    local val
    obj = parent[name]
    try
        val = read(obj)
    finally
        close(obj)
    end
    readas(val)
end
read(parent::Union{JldFile,JldGroup}, name::Symbol) = read(parent, String(string(name)))

function read(obj::JldGroup)
    nms = names(obj)
    val = Dict{String, Any}()
    for nm in nms
        val[nm] = read(obj[nm])
    end
    return val
end

function read(obj::JldDataset)
    dtype = datatype(obj.plain)
    dspace_id = HDF5.h5d_get_space(obj.plain)
    extent_type = HDF5.h5s_get_simple_extent_type(dspace_id)
    try
        if extent_type == HDF5.H5S_SCALAR
            # Scalar value
            return read_scalar(obj, dtype, jldatatype(file(obj), dtype))
        elseif extent_type == HDF5.H5S_SIMPLE
            return read_array(obj, dtype, dspace_id, HDF5.H5S_ALL)
        elseif extent_type == HDF5.H5S_NULL
            # Empty array
            if HDF5.h5t_get_class(dtype) == HDF5.H5T_REFERENCE
                T = refarray_eltype(obj)
            else
                T = jldatatype(file(obj), dtype)
            end
            if exists(obj, "dims")
                dims = a_read(obj.plain, "dims")
                return Array{T}(dims...)
            else
                return T[]
            end
        end
    finally
        HDF5.h5s_close(dspace_id)
    end
end

## Scalars
read_scalar(obj::JldDataset, dtype::HDF5Datatype, ::Type{T}) where {T<:BitsKindOrString} =
    read(obj.plain, T)
function read_scalar(obj::JldDataset, dtype::HDF5Datatype, T::Type)
    buf = Vector{UInt8}(sizeof(dtype))
    HDF5.readarray(obj.plain, dtype.id, buf)
    sc = readas(jlconvert(T, file(obj), pointer(buf)))
    gcuse(buf)
    sc
end

## Arrays

# Read an array
function read_array(obj::JldDataset, dtype::HDF5Datatype, dspace_id::HDF5.Hid, dsel_id::HDF5.Hid,
                    dims::(Tuple{Vararg{Int}})=convert((Tuple{Vararg{Int}}), HDF5.h5s_get_simple_extent_dims(dspace_id)[1]))
    if HDF5.h5t_get_class(dtype) == HDF5.H5T_REFERENCE
        val = read_refs(obj, refarray_eltype(obj), dspace_id, dsel_id, dims)
    else
        val = read_vals(obj, dtype, jldatatype(file(obj), dtype), dspace_id, dsel_id, dims)
    end
    readas(val)
end

# Arrays of basic HDF5 kinds
function read_vals(obj::JldDataset, dtype::HDF5Datatype, T::Union{Type{S}, Type{Complex{S}}},
                   dspace_id::HDF5.Hid, dsel_id::HDF5.Hid, dims::Tuple{Vararg{Int}}) where {S<:HDF5BitsKind}
    if obj.file.mmaparrays && HDF5.iscontiguous(obj.plain) && dsel_id == HDF5.H5S_ALL
        readmmap(obj.plain, Array{T})
    else
        out = Array{T}(dims)
        HDF5.h5d_read(obj.plain.id, dtype.id, dspace_id, dsel_id, HDF5.H5P_DEFAULT, out)
        out
    end
end

# Arrays of immutables/bitstypes
function read_vals(obj::JldDataset, dtype::HDF5Datatype, T::Type, dspace_id::HDF5.Hid,
                   dsel_id::HDF5.Hid, dims::Tuple{Vararg{Int}})
    out = Array{T}(dims)
    # Empty objects don't need to be read at all
    T.size == 0 && !T.mutable && return out

    # Read from file
    n = prod(dims)
    h5sz = sizeof(dtype)
    buf = Vector{UInt8}(h5sz*n)
    HDF5.h5d_read(obj.plain.id, dtype.id, dspace_id, dsel_id, HDF5.H5P_DEFAULT, buf)

    f = file(obj)
    h5offset = pointer(buf)
    if datatype_pointerfree(T) && !T.mutable
        jloffset = pointer(out)
        jlsz = T.size

        # Perform conversion in buffer
        for i = 1:n
            jlconvert!(jloffset, T, f, h5offset)
            jloffset += jlsz
            h5offset += h5sz
        end
    else
        # Convert each item individually
        for i = 1:n
            out[i] = jlconvert(T, f, h5offset)
            h5offset += h5sz
        end
    end
    gcuse(buf)
    gcuse(out)
    out
end

# Arrays of references
function read_refs(obj::JldDataset, ::Type{T}, dspace_id::HDF5.Hid, dsel_id::HDF5.Hid,
                   dims::Tuple{Vararg{Int}}) where T
    refs = Array{HDF5ReferenceObj}(dims)
    HDF5.h5d_read(obj.plain.id, HDF5.H5T_STD_REF_OBJ, dspace_id, dsel_id, HDF5.H5P_DEFAULT, refs)

    out = Array{T}(dims)
    f = file(obj)
    for i = 1:length(refs)
        if refs[i] != HDF5.HDF5ReferenceObj_NULL
            out[i] = read_ref(f, refs[i])
        end
    end
    out
end

# Get element type of a reference array
function refarray_eltype(obj::JldDataset)
    typename = a_read(obj.plain, "julia eltype")
    T = julia_type(typename)
    if T == UnsupportedType
        warn("type $typename not present in workspace; interpreting array as Array{Any}")
        return Any
    end
    return T
end

## Reference
function read_ref(f::JldFile, ref::HDF5ReferenceObj)
    if haskey(f.jlref, ref)
        # Stored as WeakRefs and may no longer exist
        val = f.jlref[ref].value
        val !== nothing && return val
    end

    dset = f[ref]
    data = try
        read(dset)
    finally
        close(dset)
    end

    f.jlref[ref] = WeakRef(data)
    readas(data)
end

### Writing ###

write(parent::Union{JldFile, JldGroup}, name::String,
      data, wsession::JldWriteSession=JldWriteSession(); kargs...) =
    close(_write(parent, name, writeas(data), wsession; kargs...))

# Pick whether to use compact or default storage based on data size
function dset_create_properties(parent, sz::Int, obj, chunk=Int[]; mmap::Bool=false)
    if sz <= 8192 && !ismmapped(parent) && !mmap
        return compact_properties(), false
    end
    if iscompressed(parent) && !isempty(chunk)
        p = p_create(HDF5.H5P_DATASET_CREATE)
        p["chunk"] = chunk
        if iscompatible(parent)
            p["shuffle"] = ()
            p["compress"] = 5
        else
            p["blosc"] = 5
        end
        return p, true
    else
        return HDF5.DEFAULT_PROPERTIES, false
    end
end

# Write "basic" types
function _write(parent::Union{JldFile, JldGroup},
                    name::String,
                    data::Union{T, Array{T}},
                    wsession::JldWriteSession; kargs...) where T<:Union{HDF5BitsKind, String}
    chunk = T <: String ? Int[] : HDF5.heuristic_chunk(data)
    dprop, dprop_close = dset_create_properties(parent, sizeof(data), data, chunk; kargs...)
    dset, dtype = d_create(parent.plain, String(name), data, HDF5._link_properties(name), dprop)
    try
        # Write the attribute
        isa(data, Array) && isempty(data) && a_write(dset, "dims", [size(data)...])
        # Write the data
        HDF5.writearray(dset, dtype.id, data)
    finally
        close(dtype)
        dprop_close && close(dprop)
    end
    dset
end

# General array types
function _write(parent::Union{JldFile, JldGroup},
                path::String, data::Array{T},
                wsession::JldWriteSession; kargs...) where T
    f = file(parent)
    dtype = h5fieldtype(f, T, true)
    buf = h5convert_array(f, data, dtype, wsession)
    dims = convert(Vector{HDF5.Hsize}, [reverse(size(data))...])
    dspace = dataspace(data)
    chunk = HDF5.heuristic_chunk(dtype, size(data))
    dprop, dprop_close = dset_create_properties(parent, sizeof(buf),buf, chunk; kargs...)
    try
        dset = d_create(parent.plain, path, dtype.dtype, dspace, HDF5._link_properties(path), dprop)
        if dtype == JLD_REF_TYPE
            a_write(dset, "julia eltype", full_typename(f, T))
        end
        if isempty(data) && ndims(data) != 1
            a_write(dset, "dims", [size(data)...])
        else
            HDF5.writearray(dset, dtype.dtype.id, buf)
        end
        return dset
    finally
        dprop_close && close(dprop)
        close(dspace)
    end
end

# Dispatch correct method for Array{Union{}}
_write(parent::Union{JldFile, JldGroup}, path::String, data::Array{Union{}},
       wsession::JldWriteSession; kargs...) =
    # Keyword arguments do not currently work
    invoke(_write, Tuple{Union{JldFile, JldGroup},String,Array,JldWriteSession}, parent,
           path, data, wsession; kargs...)

# Convert an array to the format to be written to the HDF5 file, either
# references or values
function h5convert_array(f::JldFile, data::Array,
                         dtype::JldDatatype, wsession::JldWriteSession)
    if dtype == JLD_REF_TYPE
        # For type stability, return as Vector{UInt8}
        refs = VERSION < v"0.7.0-DEV.2083" ? reinterpret(UInt8,Vector{HDF5ReferenceObj}(length(data))) : Vector{UInt8}(length(data)*sizeof(HDF5ReferenceObj))
        arefs = reinterpret(HDF5ReferenceObj, refs)
        for i = 1:length(data)
            if isassigned(data, i)
                arefs[i] = write_ref(f, data[i], wsession)
            else
                arefs[i] = HDF5.HDF5ReferenceObj_NULL
            end
        end
        refs
    else
        gen_h5convert(f, eltype(data))
        h5convert_vals(f, data, dtype, wsession)
    end
end

# Hack to ensure that _h5convert_vals isn't compiled before h5convert!
function h5convert_vals(f::JldFile, @nospecialize(data), dtype::JldDatatype,
                        wsession::JldWriteSession)
    _h5convert_vals(f, data, dtype, wsession)
end

# Convert an array of immutables or bitstypes to a buffer representing
# HDF5 compound objects. A separate function so that it is specialized.
@noinline function _h5convert_vals(f::JldFile, data::Array,
                         dtype::JldDatatype, wsession::JldWriteSession)
    sz = HDF5.h5t_get_size(dtype)
    n = length(data)
    buf = Vector{UInt8}(sz*n)
    offset = pointer(buf)
    for i = 1:n
        h5convert!(offset, f, data[i], wsession)
        offset += sz
    end
    gcuse(buf)
    buf
end

# Get reference group, creating a new one if necessary
function get_gref(f::JldFile)
    isdefined(f, :gref) && return f.gref::JldGroup

    if !exists(f, pathrefs)
        gref = f.gref = g_create(f, pathrefs)
    else
        gref = f.gref = f[pathrefs]
    end
    f.nrefs = length(gref)
    gref
end

# Write a reference
function write_ref(parent::JldFile, data, wsession::JldWriteSession)
    # Check whether we have already written this object
    ref = get(wsession.h5ref, data, HDF5.HDF5ReferenceObj_NULL)
    ref != HDF5.HDF5ReferenceObj_NULL && return ref

    # Write an new reference
    gref = get_gref(parent)
    name = @sprintf "%08d" (parent.nrefs += 1)
    dset = _write(gref, name, writeas(data), wsession)

    # Add reference to reference list
    ref = HDF5ReferenceObj(HDF5.objinfo(dset).addr)
    close(dset)
    if !isa(data, Tuple) && typeof(data).mutable
        wsession.h5ref[data] = ref
    end
    ref
end
write_ref(parent::JldGroup, data, wsession::JldWriteSession) =
    write_ref(file(parent), data, wsession)

# Expressions, drop line numbers
function _write(parent::Union{JldFile, JldGroup},
                name::String, ex::Expr,
                wsession::JldWriteSession; kargs...)
    args = ex.args
    # Discard "line" expressions
    keep = trues(length(args))
    for i = 1:length(args)
        if (isa(args[i], Expr) && args[i].head == :line) || isa(args[i], LineNumberNode)
            keep[i] = false
        end
    end
    newex = Expr(ex.head)
    newex.args = args[keep]
    write_compound(parent, name, newex, wsession)
end

# Generic (tuples, immutables, and compound types)
_write(parent::Union{JldFile, JldGroup}, name::String, s,
      wsession::JldWriteSession; kargs...) =
    write_compound(parent, name, s, wsession)
function write_compound(parent::Union{JldFile, JldGroup}, name::String,
                        s, wsession::JldWriteSession; kargs...)
    T = typeof(s)
    f = file(parent)
    dtype = h5type(f, T, true)
    gen_h5convert(f, T)

    buf = Vector{UInt8}(HDF5.h5t_get_size(dtype))
    h5convert!(pointer(buf), file(parent), s, wsession)
    gcuse(buf)

    dspace = HDF5Dataspace(HDF5.h5s_create(HDF5.H5S_SCALAR))
    dprop, dprop_close = dset_create_properties(parent, length(buf), buf; kargs...)
    try
        dset = HDF5.d_create(parent.plain, name, dtype.dtype, dspace, HDF5._link_properties(name), dprop)
        HDF5.writearray(dset, dtype.dtype.id, buf)
        return dset
    finally
        dprop_close && close(dprop)
        close(dspace)
    end
end

### Size, length, etc ###
size(dset::JldDataset) = size(dset.plain)
size(dset::JldDataset, d) = size(dset.plain, d)
length(dset::JldDataset) = prod(size(dset))
endof(dset::JldDataset) = length(dset)
ndims(dset::JldDataset) = ndims(dset.plain)

### Read/write via getindex/setindex! ###
function getindex(dset::JldDataset, indices::Union{AbstractRange{Int},Integer}...)
    sz = map(length, indices)
    dsel_id = HDF5.hyperslab(dset.plain, indices...)
    try
        dspace = HDF5._dataspace(sz)
        try
            return read_array(dset, datatype(dset.plain), dspace.id, dsel_id, sz)
        finally
            close(dspace)
        end
    finally
        HDF5.h5s_close(dsel_id)
    end
end

function setindex!(dset::JldDataset, X::AbstractArray{T,N}, indices::Union{AbstractRange{Int},Integer}...) where {T,N}
    f = file(dset)
    sz = map(length, indices)
    dsel_id = HDF5.hyperslab(dset.plain, indices...)
    try
        dtype = datatype(dset.plain)
        try
            # Convert array to writeable buffer
            if HDF5.h5t_get_class(dtype) == HDF5.H5T_REFERENCE
                written_eltype = refarray_eltype(dset)
                jldtype = JLD_REF_TYPE
            else
                written_eltype = jldatatype(f, dtype)
                jldtype = JldDatatype(dtype, -1)
            end

            buf = h5convert_array(f, convert(Array{written_eltype,N}, X), jldtype,
                                  JldWriteSession())

            dspace = HDF5._dataspace(sz)
            try
                HDF5.h5d_write(dset.plain.id, dtype, dspace, dsel_id, HDF5.H5P_DEFAULT, buf)
            finally
                close(dspace)
            end
        finally
            close(dtype)
        end
    finally
        HDF5.h5s_close(dsel_id)
    end
end
function setindex!(dset::JldDataset, x::Number, indices::Union{AbstractRange{Int},Integer}...)
    setindex!(dset, fill(x, map(length, indices)), indices...)
end

getindex(dset::JldDataset, I::Union{AbstractRange{Int},Integer,Colon}...) = getindex(dset, ntuple(i-> isa(I[i], Colon) ? (1:size(dset,i)) : I[i], length(I))...)
setindex!(dset::JldDataset, x, I::Union{AbstractRange{Int},Integer,Colon}...) = setindex!(dset, x, ntuple(i-> isa(I[i], Colon) ? (1:size(dset,i)) : I[i], length(I))...)

length(x::Union{JldFile, JldGroup}) = length(names(x))

### Custom serialization

readas(x) = x
writeas(x) = x

# Wrapper for associative keys
# We write this instead of the associative to avoid dependence on the
# Julia hash function
struct AssociativeWrapper{K,V,T<:Associative}
    keys::Vector{K}
    values::Vector{V}
end

readas(x::AssociativeWrapper{K,V,T}) where {K,V,T} = convert(T, x)
function writeas(x::T) where T<:Associative
    K, V = destructure(eltype(x))
    convert(AssociativeWrapper{K,V,T}, x)
end
destructure(::Type{Pair{K,V}}) where {K,V} = K, V  # not inferrable, julia#10880

# Special case for associative, to rehash keys
function convert(::Type{T}, x::AssociativeWrapper{K,V,T}) where {K,V,T<:Associative}
    ret = T()
    keys = x.keys
    values = x.values
    n = length(keys)
    if applicable(sizehint!, ret, n)
        sizehint!(ret, n)
    end
    for i = 1:n
        ret[keys[i]] = values[i]
    end
    ret
end

function convert(::Type{AssociativeWrapper{K,V,T}}, d::Associative) where {K,V,T}
    n = length(d)
    ks = Vector{K}(n)
    vs = Vector{V}(n)
    i = 0
    for (k,v) in d
        ks[i+=1] = k
        vs[i] = v
    end
    AssociativeWrapper{K,V,typeof(d)}(ks, vs)
end

# Special case for SimpleVector
# Wrapper for SimpleVector
struct SimpleVectorWrapper
    elements::Vector
end

# Special case for SimpleVector
readas(x::SimpleVectorWrapper) = Core.svec(x.elements...)
writeas(x::SimpleVector) = SimpleVectorWrapper([x...])

# function to convert string(mod::Module) back to mod::Module
function modname2mod(modname::AbstractString)
    parse(modname == "Main" ? modname : string("Main.", modname))
end

# Serializer for anonymous functions
# convert functions to lowered ast expressions
function func2expr(fun::Function)
    @assert !isgeneric(fun) "generic functions not supported"
    ast = Base.uncompressed_ast(fun.code)
    Expr(:function, Expr(:tuple, ast.args[1]...), Expr(:block, ast.args[3].args...))
end
struct AnonymousFunctionSerializer
    expr::Expr
    mod::AbstractString
    AnonymousFunctionSerializer(fun::Function) = new(func2expr(fun), string(fun.code.module))
end
readas(ast::AnonymousFunctionSerializer) = eval(modname2mod(ast.mod)).eval(ast.expr)
writeas(fun::Function) = AnonymousFunctionSerializer(fun)

# Serializer for GlobalRef
struct GlobalRefSerializer
    mod::AbstractString
    name::Symbol
    GlobalRefSerializer(g::GlobalRef) = new(string(g.mod), g.name)
end
readas(grs::GlobalRefSerializer) = GlobalRef(eval(modname2mod(grs.mod)), grs.name)
writeas(gr::GlobalRef) = GlobalRefSerializer(gr)

# StackFrame (Null the LambdaInfo in 0.5)
# or Core.MethodInstance in 0.6
JLD.writeas(data::StackFrame) =
    Base.StackFrame(data.func,
                    data.file,
                    data.line,
                    Nullable{Core.MethodInstance}(),
                    data.from_c,
                    data.inlined,
                    data.pointer)

### Converting attribute strings to Julia types

const _where_macrocall = Symbol("@where")
function expand_where_macro(e::Expr)
    e.head = :where
    shift!(e.args)
    Compat.macros_have_sourceloc && shift!(e.args)
    return true
end

is_valid_type_ex(s::Symbol) = true
is_valid_type_ex(s::QuoteNode) = true
is_valid_type_ex(s) = isbits(typeof(s))
function is_valid_type_ex(e::Expr)
    if e.head === :curly || e.head == :tuple || e.head == :.
        return all(is_valid_type_ex, e.args)
    elseif e.head === :where
        return is_valid_type_ex(e.args[1])
    elseif e.head === :let && length(e.args) == 2
        return is_valid_type_ex(e.args[1]) &&
               is_valid_type_ex(e.args[2].args[2])
    elseif e.head == :call
        f = e.args[1]
        if f isa Expr
            if f.head === :core
                f = f.args[1]
                return f === :Union || f === :TypeVar || f === :UnionAll
            end
        elseif f isa Symbol
            return f === :Union || f === :TypeVar || f === :symbol
        end
    end
    return false
end

const typemap_Core = Dict(
    :Uint8 => :UInt8,
    :Uint16 => :Uint16,
    :Uint32 => :UInt32,
    :Uint64 => :UInt64,
    :Nothing => :Void
)

const _typedict = Dict{String,Type}()

function fixtypes(typ)
    whereall = []
    typ = fixtypes(typ, whereall)
    while !isempty(whereall)
        var = pop!(whereall)
        typ = Expr(:let, Expr(:call, Expr(:core, :UnionAll), var.args[1], typ), var)
    end
    return typ
end
fixtypes(typ, whereall) = typ
function fixtypes(typ::Expr, whereall::Vector{Any})
    if typ.head === :macrocall && typ.args[1] === _where_macrocall
        expand_where_macro(typ) # @where => TypeVar format forwards compatibility
    end
    if typ.head === :.
        if length(typ.args) == 2 && typ.args[1] === :Core
            arg = typ.args[2].value
            return Expr(:., :Core, QuoteNode(get(typemap_Core, arg, arg)))
        else
            return typ
        end
    elseif typ == :(Core.Type{TypeVar(:T,Union(Core.Any,Core.Undef))}) || typ == :(Core.Type{TypeVar(:T)})
        # Work around https://github.com/JuliaLang/julia/issues/8226 and the removal of Top
        return :(Core.Type)
    end

    for i = 1:length(typ.args)
        typ.args[i] = fixtypes(typ.args[i], whereall)
    end

    if (typ.head === :call && !isempty(typ.args) &&
        typ.args[1] === :TypeVar) # TypeVar => where format backwards compatibility
        tv = gensym()
        push!(whereall, Expr(:(=), tv, typ))
        return tv
    end

    if (typ.head === :call && !isempty(typ.args) &&
        typ.args[1] === :Union)
        typ = Expr(:curly, typ.args...)
    end

    if typ.head === :tuple
        typ = Expr(:curly, :Tuple, typ.args...)
    end

    if typ.head === :curly
        # assume literal TypeVar should work like `T{<:S}`
        while !isempty(whereall)
            var = pop!(whereall)
            typ = Expr(:let, Expr(:call, Expr(:core, :UnionAll), var.args[1], typ), var)
        end
    end
    return typ
end

function _julia_type(s::AbstractString)
    typ = get(_typedict, s, UnconvertedType)
    if typ == UnconvertedType
        sp = parse(s, raise=false)
        if (isa(sp, Expr) && (sp.head == :error || sp.head == :continue || sp.head == :incomplete))
            println("error parsing type string ", s)
            eval(sp)
        end
        typ = julia_type(fixtypes(sp))
        if typ != UnsupportedType
            _typedict[s] = typ
        end
    end
    typ
end

function julia_type(e::Union{Symbol, Expr})
    if is_valid_type_ex(e)
        try # `try` needed to catch undefined symbols
            # `e` should be fully qualified, and thus reachable from Main
            typ = eval(Main, e)
            typ == Type && return Type
            isa(typ, Type) && return typ
        end
    end
    return UnsupportedType
end

### Converting Julia types to fully qualified names
function full_typename(io::IO, file::JldFile, jltype::Union)
    print(io, "Union(")
    types = Base.uniontypes(jltype)
    if !isempty(types)
        full_typename(io, file, types[1])
        for i = 2:length(types)
            print(io, ',')
            full_typename(io, file, types[i])
        end
    end
    print(io, ')')
end
function full_typename(io::IO, file::JldFile, x::typeof(Union{}))
    print(io, "Union()")
end
function full_typename(io::IO, file::JldFile, x::UnionAll)
    x == Type && return print(io, "Type")
    print(io, "@where(")
    full_typename(io, file, x.body)
    print(io, ',')
    tv = x.var
    if tv.lb === Union{} && tv.ub === Any
        print(io, tv.name)
    elseif tv.lb === Union{}
        print(io, tv.name)
        print(io, "<:")
        full_typename(io, file, tv.ub)
    else
        full_typename(io, file, tv.lb)
        print(io, "<:")
        print(io, tv.name)
        print(io, "<:")
        full_typename(io, file, tv.ub)
    end
    print(io, ')')
end
function full_typename(io::IO, file::JldFile, tv::TypeVar)
    print(io, tv.name)
end
function full_typename(io::IO, ::JldFile, x)
    # Only allow bitstypes that show as AST literals and make sure that they
    # read back exactly as they are saved. Use show here (instead of print) to
    # preserve as many Julian type distinctions as we can e.g., 1 vs 0x01
    # A different implementation will be required to support custom immutables
    # or things as simple as Int16(1).
    s = sprint(show, x)
    if isbits(x) && parse(s) === x && !isa(x, Tuple)
        print(io, s)
    else
        error("type parameters with objects of type ", typeof(x), " are currently unsupported")
    end
end
function full_typename(io::IO, ::JldFile, x::Symbol)
    s = string(x)
    if contains(s, " ")
        # escape spaces
        print_escaped(io, string("symbol(\"", string(x), "\")"), " ")
    else
        print(io, ":", x)
    end
end
function full_typename(io::IO, file::JldFile, jltype::DataType)
    mod = jltype.name.module
    if mod != Main
        mname = string(mod)
        for x in file.truncatemodules
            if startswith(mname, x)
                mname = length(x) == length(mname) ? "" : mname[sizeof(x)+1:end]
                if startswith(mname, ".")
                    mname = mname[2:end]
                end
                break
            end
        end

        if !isempty(mname)
            print(io, mname)
            print(io, '.')
        end
    end

    print(io, jltype.name.name)
    if !isempty(jltype.parameters)
        print(io, '{')
        full_typename(io, file, jltype.parameters[1])
        for i = 2:length(jltype.parameters)
            print(io, ',')
            full_typename(io, file, jltype.parameters[i])
        end
        print(io, '}')
    elseif jltype <: Tuple
        print(io, "{}")
    end
end
function full_typename(file::JldFile, x)
    io = IOBuffer(Vector{UInt8}(64), true, true)
    truncate(io, 0)
    full_typename(io, file, x)
    String(take!(io))
end

function truncate_module_path(file::JldFile, mod::Module)
    push!(file.truncatemodules, string(mod))
end

function names(parent::Union{JldFile, JldGroup})
    n = names(parent.plain)
    keep = trues(length(n))
    reserved = [pathrefs[2:end], pathtypes[2:end], pathrequire[2:end], pathcreator[2:end]]
    for i = 1:length(n)
        if in(n[i], reserved)
            keep[i] = false
        end
    end
    n[keep]
end

function save_write(f, s, vname, wsession::JldWriteSession)
    if !isa(vname, Function)
        try
            write(f, s, vname)
        catch e
            if isa(e, PointerException)
                warn("Skipping $vname because it contains a pointer")
            end
        end
    end
end

macro save(filename, vars...)
    if isempty(vars)
        # Save all variables in the current module
        writeexprs = Vector{Expr}(0)
        m = current_module()
        for vname in names(m)
            s = string(vname)
            if !ismatch(r"^_+[0-9]*$", s) # skip IJulia history vars
                v = eval(m, vname)
                if !isa(v, Module)
                    push!(writeexprs, :(save_write(f, $s, $(esc(vname)), wsession)))
                end
            end
        end
    else
        writeexprs = Vector{Expr}(length(vars))
        for i = 1:length(vars)
            writeexprs[i] = :(write(f, $(string(vars[i])), $(esc(vars[i])), wsession))
        end
    end

    quote
        local f = jldopen($(esc(filename)), "w")
        wsession = JldWriteSession()
        try
            $(Expr(:block, writeexprs...))
        finally
            close(f)
        end
    end
end

macro load(filename, vars...)
    if isempty(vars)
        if isa(filename, Expr)
            warn("""@load-ing a file without specifying the variables to be loaded may produce
                    unexpected behavior unless the file is specified as a string literal. Future
                    versions of JLD will require that the file is specified as a string literal
                    in this case.""")
            filename = eval(current_module(), filename)
        end
        # Load all variables in the top level of the file
        readexprs = Expr[]
        vars = Symbol[]
        f = jldopen(filename)
        try
            for n in names(f)
                obj = f[n]
                try
                    if isa(obj, JldDataset)
                        push!(vars, Symbol(n))
                    end
                finally
                    close(obj)
                end
            end
        finally
            close(f)
        end
    end
    return quote
        f = jldopen($(esc(filename)))
        ($([esc(x) for x in vars]...),) = try
            ($([:(read(f, $(string(x)))) for x in vars]...),)
        finally
            close(f)
        end
        $(Symbol[v for v in vars]) # convert to Array
    end
end

# Save all the key-value pairs in the dict as top-level variables of the JLD
function FileIO.save(f::File{format"JLD"}, dict::Associative; compatible::Bool=false, compress::Bool=false)
    jldopen(FileIO.filename(f), "w"; compatible=compatible, compress=compress) do file
        wsession = JldWriteSession()
        for (k,v) in dict
            if !isa(k, AbstractString)
                error("Keys must be strings (the names of variables), got $k")
            end
            write(file, String(k), v, wsession)
        end
    end
end
# Or the names and values may be specified as alternating pairs
function FileIO.save(f::File{format"JLD"}, name::AbstractString, value, pairs...; compatible::Bool=false, compress::Bool=false)
    if isodd(length(pairs)) || !isa(pairs[1:2:end], Tuple{Vararg{AbstractString}})
        throw(ArgumentError("arguments must be in name-value pairs"))
    end
    jldopen(FileIO.filename(f), "w"; compatible=compatible, compress=compress) do file
        wsession = JldWriteSession()
        write(file, String(name), value, wsession)
        for i=1:2:length(pairs)
            write(file, String(pairs[i]), pairs[i+1], wsession)
        end
    end
end
FileIO.save(f::File{format"JLD"}, value...; kwargs...) = error("Must supply a name for each variable")

# load with just a filename returns a dictionary containing all the variables
function FileIO.load(f::File{format"JLD"})
    jldopen(FileIO.filename(f), "r") do file
        Dict{String,Any}([(var, read(file, var)) for var in names(file)])
    end
end
# When called with explicitly requested variable names, return each one
function FileIO.load(f::File{format"JLD"}, varname::AbstractString)
    jldopen(FileIO.filename(f), "r") do file
        read(file, varname)
    end
end
FileIO.load(f::File{format"JLD"}, varnames::AbstractString...) = load(f, varnames)
function FileIO.load(f::File{format"JLD"}, varnames::Tuple{Vararg{AbstractString}})
    jldopen(FileIO.filename(f), "r") do file
        map((var)->read(file, var), varnames)
    end
end

# As of this version, packages aren't loaded into Main by default, so the root
# module check verifies that packages are still identified as being top level
# even if a binding to them is not present in Main.
if VERSION >= v"0.7.0-DEV.1877"
    _istoplevel(m::Module) = module_parent(m) == Main || Base.is_root_module(m)
else
    _istoplevel(m::Module) = module_parent(m) == Main
end

function addrequire(file::JldFile, mod::Module)
    _istoplevel(mod) || error("must be a toplevel module")
    addrequire(file, module_name(mod))
end

function addrequire(file::JldFile, modsym::Symbol)
    has_require = exists(file, pathrequire)
    modules = has_require ? map(path2modsym, read(file, pathrequire)) : Symbol[]
    push!(modules, modsym)
    has_require && o_delete(file, pathrequire)
    write(file, pathrequire, modules)
end

function addrequire(file::JldFile, filename::AbstractString)
    warn("\"addrequire(file, filename)\" is deprecated, please use \"addrequire(file, module)\"")
    addrequire(file, path2modsym(filename))
end

# Cope with JLD0.1 format
path2modsym(s::Symbol) = s
function path2modsym(filename::AbstractString)
    bname = basename(filename)
    if endswith(bname, ".jl")
        bname = bname[1:end-3]
    end
    Symbol(bname)
end

export
    addrequire,
    creator,
    ismmappable,
    jldopen,
    o_delete,
    readmmap,
    @load,
    @save,
    load,
    save,
    translate,
    truncate_module_path

const _runtime_properties = Ref{HDF5.HDF5Properties}()
compact_properties() = _runtime_properties[]

function __init__()
    COMPACT_PROPERTIES = p_create(HDF5.H5P_DATASET_CREATE)
    HDF5.h5p_set_layout(COMPACT_PROPERTIES.id, HDF5.H5D_COMPACT)

    global _runtime_properties
    _runtime_properties[] = COMPACT_PROPERTIES

    HDF5.rehash!(_typedict, length(_typedict.keys))
    HDF5.rehash!(BUILTIN_TYPES.dict, length(BUILTIN_TYPES.dict.keys))

    nothing
end

include("JLD00.jl")
end


text:
using MyFirstPackage
using Test

@testset "MyFirstPackage.jl" begin
    @test my_f(2,1) == 7
    @test my_f(2,3) == 13
end

ext: jl


text:
#! /usr/bin/julia

# Rosetta Code, Probabilistic choice

p = [1/i for i in 5:11]
plen = length(p)
q = [0.0, [sum(p[1:i]) for i = 1:plen]]
plab = [char(i) for i in 0x05d0:(0x05d0+plen)]
hi = 10^6
push!(p, 1.0 - sum(p))
plen += 1

accum = zeros(Int, plen)

for i in 1:hi
    accum[sum(rand() .>= q)] += 1
end

r = accum/hi

println("Rates at which items are selected (", hi, " trials).")
println(" Item  Expected   Actual")
for i in 1:plen
    println(@sprintf("   \u2067%s   %8.6f  %8.6f", plab[i], p[i], r[i]))
end

println()
println("Rates at which items are selected (", hi, " trials).")
println(" Item         Count   Expected   Actual")
for i in 1:plen
    println(@sprintf("   %s yields  %6d   %8.6f  %8.6f",
                     plab[i], accum[i], p[i], r[i]))
end

max_stars_repo_name: MichaeLeroy/rosetta-code

max_issues_repo_name: MichaeLeroy/rosetta-code

max_forks_repo_name: MichaeLeroy/rosetta-code


text:
function beta(obj::GNSSFieldAttitudeContext, arg0::FieldAbsoluteDate)
    return jcall(obj, "beta", RealFieldElement, (FieldAbsoluteDate,), arg0)
end

function beta_ds(obj::GNSSFieldAttitudeContext)
    return jcall(obj, "betaDS", FieldDerivativeStructure, ())
end

function compute_phi(obj::GNSSFieldAttitudeContext, arg0::RealFieldElement, arg1::RealFieldElement)
    return jcall(obj, "computePhi", RealFieldElement, (RealFieldElement, RealFieldElement), arg0, arg1)
end

function get_date(obj::GNSSFieldAttitudeContext)
    return jcall(obj, "getDate", FieldAbsoluteDate, ())
end

function get_delta_ds(obj::GNSSFieldAttitudeContext)
    return jcall(obj, "getDeltaDS", FieldDerivativeStructure, ())
end

function get_mu_rate(obj::GNSSFieldAttitudeContext)
    return jcall(obj, "getMuRate", RealFieldElement, ())
end

function get_orbit_angle_since_midnight(obj::GNSSFieldAttitudeContext)
    return jcall(obj, "getOrbitAngleSinceMidnight", RealFieldElement, ())
end

function get_secured_beta(obj::GNSSFieldAttitudeContext)
    return jcall(obj, "getSecuredBeta", RealFieldElement, ())
end

function get_sv_bcos(obj::GNSSFieldAttitudeContext)
    return jcall(obj, "getSVBcos", RealFieldElement, ())
end

function get_turn_span(obj::GNSSFieldAttitudeContext)
    return jcall(obj, "getTurnSpan", FieldTurnSpan, ())
end

function get_yaw_end(obj::GNSSFieldAttitudeContext, arg0::RealFieldElement)
    return jcall(obj, "getYawEnd", RealFieldElement, (RealFieldElement,), arg0)
end

function get_yaw_start(obj::GNSSFieldAttitudeContext, arg0::RealFieldElement)
    return jcall(obj, "getYawStart", RealFieldElement, (RealFieldElement,), arg0)
end

function in_orbit_plane_absolute_angle(obj::GNSSFieldAttitudeContext, arg0::RealFieldElement)
    return jcall(obj, "inOrbitPlaneAbsoluteAngle", RealFieldElement, (RealFieldElement,), arg0)
end

function in_sun_side(obj::GNSSFieldAttitudeContext)
    return jcall(obj, "inSunSide", jboolean, ())
end

function in_turn_time_range(obj::GNSSFieldAttitudeContext)
    return jcall(obj, "inTurnTimeRange", jboolean, ())
end

function linear_model_still_active(obj::GNSSFieldAttitudeContext, arg0::RealFieldElement, arg1::RealFieldElement)
    return jcall(obj, "linearModelStillActive", jboolean, (RealFieldElement, RealFieldElement), arg0, arg1)
end

function nominal_yaw(obj::GNSSFieldAttitudeContext, arg0::FieldAbsoluteDate)
    return jcall(obj, "nominalYaw", TimeStampedFieldAngularCoordinates, (FieldAbsoluteDate,), arg0)
end

function orbit_normal_yaw(obj::GNSSFieldAttitudeContext)
    return jcall(obj, "orbitNormalYaw", TimeStampedFieldAngularCoordinates, ())
end

function set_half_span(obj::GNSSFieldAttitudeContext, arg0::RealFieldElement, arg1::jdouble)
    return jcall(obj, "setHalfSpan", void, (RealFieldElement, jdouble), arg0, arg1)
end

function set_up_turn_region(obj::GNSSFieldAttitudeContext, arg0::jdouble, arg1::jdouble)
    return jcall(obj, "setUpTurnRegion", jboolean, (jdouble, jdouble), arg0, arg1)
end

function time_since_turn_start(obj::GNSSFieldAttitudeContext)
    return jcall(obj, "timeSinceTurnStart", RealFieldElement, ())
end

function turn_corrected_attitude(obj::GNSSFieldAttitudeContext, arg0::FieldDerivativeStructure)
    return jcall(obj, "turnCorrectedAttitude", TimeStampedFieldAngularCoordinates, (FieldDerivativeStructure,), arg0)
end

function turn_corrected_attitude(obj::GNSSFieldAttitudeContext, arg0::RealFieldElement, arg1::RealFieldElement)
    return jcall(obj, "turnCorrectedAttitude", TimeStampedFieldAngularCoordinates, (RealFieldElement, RealFieldElement), arg0, arg1)
end

function yaw_rate(obj::GNSSFieldAttitudeContext, arg0::RealFieldElement)
    return jcall(obj, "yawRate", RealFieldElement, (RealFieldElement,), arg0)
end


text:
"""
    apply_kernel_to_field(k,f...) -> Field

Returns a field obtained by applying the kernel `k` to the
values of the fields in `f`. That is, the returned field evaluated at
a vector of points `x` provides the value obtained by applying kernel `k` to the
values of the fields `f` at the vector of points `x`. Formally, the resulting field at a
vector of points
 `x` is defined as

    fx = evaluate_fields(f,x)
    apply_kernel(k,fx...)

In order to be able to call the [`field_gradient`](@ref) function of the
resulting field, one needs to define the gradient operator
associated with the underlying kernel.
This is done by adding a new method to [`apply_kernel_gradient(k,f...)`](@ref) for each kernel type.
"""
@inline function apply_kernel_to_field(k,f...)
  AppliedField(k,f...)
end

"""
    apply_kernel_gradient(k,f...)

Returns a field representing the gradient of the field obtained with

    apply_kernel_to_field(k,f...)
"""
function apply_kernel_gradient(k,f...)
  @abstractmethod
end

# Result of applying a kernel to the value of some fields

struct AppliedField{K,F} <: Field
  k::K
  f::F
  @inline function AppliedField(k,f...)
    new{typeof(k),typeof(f)}(k,f)
  end
end

function field_return_type(f::AppliedField,x)
  Ts = field_return_types(f.f,x)
  kernel_return_type(f.k, testvalues(Ts...)...)
end

function field_cache(f::AppliedField,x)
  cf = field_caches(f.f,x)
  fx = evaluate_fields!(cf,f.f,x)
  ck = kernel_cache(f.k,fx...)
  (ck,cf)
end

@inline function evaluate_field!(cache,f::AppliedField,x)
  ck, cf = cache
  fx = evaluate_fields!(cf,f.f,x)
  apply_kernel!(ck,f.k,fx...)
end

@inline function field_gradient(f::AppliedField)
  apply_kernel_gradient(f.k,f.f...)
end


text:
include("./qol.jl")
include("./timer.jl")


text:
# SPDX-License-Identifier: X11
# 2020-11-23
# Determinant (100pt)

a,b = parse.(Int, (split ∘ readline)())
c,d = parse.(Int, (split ∘ readline)())

println(a*d - b*c)

max_stars_repo_name: jthvai/competitive-programming

max_issues_repo_name: jthvai/competitive-programming

max_forks_repo_name: jthvai/competitive-programming


text:
using JLD   # 数据存取
using Gadfly    # 绘图
using DataFrames # 数据分类
using Cairo # PDF保存
using Statistics    # 计算置信区间

# ]add JLD, Gadfly, DataFrames, Cairo, Statistics

# 需要计算置信区间的函数
function error_95(data)
    1.96*std(data)/sqrt(length(data))
end

function isequ(a, b)
    if a == b
        return true
    else
        return false
    end
end

r = load("MBL_ED_test.jld")
rDF = r["result"]
# rDF = filter(:E => E -> !(ismissing(E) || isnothing(E) || isnan(E)), rDF)
rDF = filter(( [:W, :Class] => (W, Class) -> (!( isequ(W, 8.5) || isequ(W, 9.0)|| isequ(W, 10.0)) || !( isequ(Class, "QP"))) ), rDF)

#########################################################
# 读取数据并封装至Dataframe
# r = load("fin_MBL_CPU.jld")
# d_c = Vector{String}(undef, length(r["result"]))
# # d_n = Vector{Int}(undef, length(r["result"]))
# d_n = Vector{String}(undef, length(r["result"]))
# d_w = Vector{Float64}(undef, length(r["result"]))
# d_r = Vector{Float64}(undef, length(r["result"]))
# d_e = Vector{Float64}(undef, length(r["result"]))
# for i in 1:length(r["result"])
#     d_c[i] = r["result"][i][1]
#     # d_n[i] = r["result"][i][2]
#     d_n[i] = string(r["result"][i][2])
#     d_w[i] = r["result"][i][3]
#     d_r[i] = r["result"][i][4]
#     d_e[i] = r["result"][i][5]
# end
# rDF = DataFrame(Class = d_c, N_tol = d_n, W = d_w, r = d_r, E = d_e)
# d_c, d_n, d_w, d_r, d_e = [nothing for i in 1:5]
# rDF = filter(:E => E -> !(ismissing(E) || isnothing(E) || isnan(E)), rDF)

# r = load("test_MBL_3.jld")
# d_c = Vector{String}(undef, length(r["result"]))
# # d_n = Vector{Int}(undef, length(r["result"]))
# d_n = Vector{String}(undef, length(r["result"]))
# d_w = Vector{Float64}(undef, length(r["result"]))
# d_r = Vector{Float64}(undef, length(r["result"]))
# d_e = Vector{Float64}(undef, length(r["result"]))
# for i in 1:length(r["result"])
#     d_c[i] = r["result"][i][1]
#     # d_n[i] = r["result"][i][2]
#     d_n[i] = string(r["result"][i][2])
#     d_w[i] = r["result"][i][3]
#     d_r[i] = r["result"][i][4]
#     d_e[i] = r["result"][i][5]
# end
# rDF_2 = DataFrame(Class = d_c, N_tol = d_n, W = d_w, r = d_r, E = d_e)
# d_c, d_n, d_w, d_r, d_e = [nothing for i in 1:5]
# rDF_2 = filter(:E => E -> !(ismissing(E) || isnothing(E) || isnan(E)), rDF_2)
# append!(rDF, rDF_2)
# rDF = filter(( [:W, :Class] => (W, Class) -> (!( isequ(W, 8.5) || isequ(W, 9.0)|| isequ(W, 10.0)) || !( isequ(Class, "QP"))) ), rDF)

#########################################################

# 数据分类，并计算平均值与置信区间
gdf = groupby(rDF, [:Class, :W, :N_tol])
fin_DF = combine(gdf, :E => error_95 , :E => mean, :r => error_95 , :r => mean)
fin_DF = groupby(fin_DF, :Class)

# 绘图
# 准备刻度
r_ticks = [0.38, 0.43, 0.48, 0.53]
e_ticks = [0.0, 0.2, 0.4, 0.6, 0.8, 1.0]
mysmoothing = 0.65 # 平滑度
# R Model plot
my_buffer = fin_DF[(Class = "R",)]
myp_Rr = plot(my_buffer, x=:W, y=:r_mean, color=:N_tol, shape=:N_tol
    , ymin = my_buffer.r_mean - my_buffer.r_error_95
    , ymax = my_buffer.r_mean + my_buffer.r_error_95
    , Geom.errorbar
    , Geom.point
    # , Geom.line
    , Guide.xticks(ticks=0:2:10)
    , Guide.yticks(ticks=r_ticks, label=true)
    , yintercept = [0.39, 0.53]
    , Geom.hline(style = [:dash, :dash], color=["gray", "gray"])
    , Guide.xlabel("W")
    , Guide.ylabel("Mean of r")
    , Guide.title("R Model")
    , layer( Stat.smooth(method=:loess, levels=[0.95], smoothing = mysmoothing), Geom.line )
    )

myp_RE = plot(my_buffer, x=:W, y=:E_mean, color=:N_tol, shape=:N_tol
    , ymin = my_buffer.E_mean - my_buffer.E_error_95
    , ymax = my_buffer.E_mean + my_buffer.E_error_95
    , Geom.errorbar
    , Geom.point
    # , Geom.line
    , Guide.xticks(ticks=0:2:10)
    , Guide.yticks(ticks=e_ticks)
    , Guide.xlabel("W")
    , Guide.ylabel("Entroy")
    , Guide.title("R Model")
    , layer( Stat.smooth(method=:loess, levels=[0.95], smoothing = mysmoothing), Geom.line )
    )

# img = SVG("r_w_test.svg", 6inch, 4inch)
# draw(img, myp)

# QP Model plot
my_buffer = fin_DF[(Class = "QP",)]

myp_QPr = plot(my_buffer, x=:W, y=:r_mean, color=:N_tol, shape=:N_tol
    , ymin = my_buffer.r_mean - my_buffer.r_error_95
    , ymax = my_buffer.r_mean + my_buffer.r_error_95
    , Geom.errorbar
    , Geom.point
    # , Geom.line
    , yintercept = [0.39, 0.53]
    , Geom.hline(style = [:dash, :dash], color=["gray", "gray"])
    , Guide.xticks(ticks=0:2:8)
    , Guide.yticks(ticks=r_ticks)
    , Guide.xlabel("W")
    , Guide.ylabel("Mean of r")
    , Guide.title("QP Model")
    , layer( Stat.smooth(method=:loess, levels=[0.95], smoothing = mysmoothing), Geom.line )
    )

myp_QPE = plot(my_buffer, x=:W, y=:E_mean, color=:N_tol, shape=:N_tol
    , ymin = my_buffer.E_mean - my_buffer.E_error_95
    , ymax = my_buffer.E_mean + my_buffer.E_error_95
    , Geom.errorbar
    , Geom.point
    # , Geom.line
    , Guide.xticks(ticks=0:2:8)
    , Guide.yticks(ticks=e_ticks)
    , Guide.xlabel("W")
    , Guide.ylabel("Entroy")
    , Guide.title("QP Model")
    , layer( Stat.smooth(method=:loess, levels=[0.95], smoothing = mysmoothing), Geom.line )
    )

my_buffer = nothing

# img = SVG("E_w_test.svg", 6inch, 4inch)
# draw(img, myp)

fin_pl = gridstack([myp_QPE myp_RE; myp_QPr myp_Rr])

# img = SVG("5_result_CUDA_test.svg", 12inch, 8inch)
img = PDF("fin_MBL_CPU.pdf", 12inch, 8inch)
# img_PNG = PNG("result_test.svg", 12inch, 8inch)
draw(img, fin_pl)


text:
function compute_lagrangian(data::Dict)
  return data["function"](data);
end

function compute_lagrangian_dictonly(data::Dict)
  x        = data["x"];
  d        = data["d"];
  s        = data["s"];
  lambda_s = data["lambda_s"];
  mu_s     = data["mu_s"];
  z        = data["z"];
  t        = data["t"];
  lambda_t = data["lambda_t"];
  mu_t     = data["mu_t"];
  beta     = data["beta"];

  (convH, convW, dictSize) = size(d);
  convLen = convH*convW;
  numImgs = size(x, 3);
  rerr = 0.0;
  spar = 0.0;
  cons = 0.0;
  d_fft = fft2(d);
  for j=1:numImgs;
    rerr += 0.5*vecnorm( x[:, :, j]-real(ifft2(sum(d_fft.*fft2(reshape(z[:, j], convH, convW, dictSize)), 3))) )^2;
    spar += beta*norm(t[:, j], 1);
    cons += mu_t*0.5*norm(z[:, j]-t[:, j])^2 + [lambda_t[:, j]'*(z[:, j]-t[:, j])][1];
  end
  obj = rerr + spar + cons + mu_s*0.5*norm(d[:]-s[:])^2 + [lambda_s[:]'*(d[:]-s[:])][1];

  return (obj, rerr, spar, 0.0);
end

max_stars_repo_name: bkong/FCSC


text:
using WMesa
using Test

@testset "WMesa.jl" begin
    # Write your tests here.
end


text:
using NewsAPI
using Test
using Dates

@testset "NewsAPI.jl" begin
    if ENV["NEWS_API_KEY"] != "0"
        @test length(NewsAPI.sources()) > 0
        @test length(NewsAPI.everything(q="bitcoin", from=today() - Day(1))) > 0
        @test length(NewsAPI.topheadlines(country="us")) > 0
    end
end


text:
# This file includes code that was formerly a part of Julia. License is MIT: http://julialang.org/license

__precompile__(true)

module LegacyStrings

export
    DirectIndexString,
    ByteString,
    ASCIIString,
    RepString,
    RevString,
    UTF8String,
    UTF16String,
    UTF32String,
    WString,
    ascii, # distinct from Base.ascii
    utf8,
    utf16,
    utf32,
    wstring

import Base:
    containsnul,
    convert,
    getindex,
    isvalid,
    length,
    lowercase,
    map,
    nextind,
    pointer,
    prevind,
    reverse,
    reverseind,
    show,
    sizeof,
    string,
    unsafe_convert,
    uppercase,
    write

using Compat
using Compat: IOBuffer
import Compat:
    lastindex,
    codeunit,
    ncodeunits

if VERSION < v"0.7-"
    import Base: lcfirst
    import Base: next
    import Base: rsearch
    import Base: search
    import Base: ucfirst
    import Base: UnicodeError
    using Base: DirectIndexString
else
    import Base: iterate
    include("unicodeerror.jl")
    include("directindex.jl")
end

struct ASCIIString <: DirectIndexString
    data::Vector{UInt8}
    ASCIIString(data::String) = new(Vector{UInt8}(codeunits(data)))
    ASCIIString(data) = new(data)
end

struct UTF8String <: AbstractString
    data::Vector{UInt8}
    UTF8String(data::String) = new(Vector{UInt8}(codeunits(data)))
    UTF8String(data) = new(data)
end

struct UTF16String <: AbstractString
    data::Vector{UInt16} # includes 16-bit NULL termination after string chars
    function UTF16String(data::Vector{UInt16})
        if length(data) < 1 || data[end] != 0
            throw(UnicodeError(UTF_ERR_NULL_16_TERMINATE, 0, 0))
        end
        new(data)
    end
end

struct UTF32String <: DirectIndexString
    data::Vector{UInt32} # includes 32-bit NULL termination after string chars
    function UTF32String(data::Vector{UInt32})
        if length(data) < 1 || data[end] != 0
            throw(UnicodeError(UTF_ERR_NULL_32_TERMINATE, 0, 0))
        end
        new(data)
    end
end

const ByteString = Union{ASCIIString,UTF8String}

include("support.jl")
include("ascii.jl")
include("utf8.jl")
include("utf16.jl")
include("utf32.jl")
include("rep.jl")

if isdefined(Base, :RevString)
    using Base: RevString
else
    include("rev.jl")
end

const AllLegacyStringTypes = Union{ASCIIString,UTF8String,UTF16String,UTF32String,RepString,RevString}

codeunit(s::SubString{<:AllLegacyStringTypes}) = codeunit(s.string)
ncodeunits(s::SubString{<:AllLegacyStringTypes}) = isdefined(s, :ncodeunits) ? s.ncodeunits : s.endof

if !isdefined(Base, :iterate)
    iterate(s::Union{String,SubString,AllLegacyStringTypes}, i::Int) = next(s, i)
end

end # module


text:
#Demo 9: Waveforms with NaN
#-------------------------------------------------------------------------------
using InspectDR
using Colors

#==Input
===============================================================================#

#Constants
#-------------------------------------------------------------------------------
black = RGB24(0, 0, 0)
white = RGB24(1, 1, 1)
red = RGB24(1, 0, 0)
green = RGB24(0, 1, 0)
blue = RGB24(0, 0, 1)

#Input data
#-------------------------------------------------------------------------------
npts = 1000
T=.1e-6 #Clock period
x = collect(range(0, stop=1e-6, length=npts))
y = sin.(2pi*(x./T))

#Invalidate some data points:
#NOTE: Cannot use F1-acceleration when x has NaN values
x[200] = y[200] = NaN #Causes break
x[460] = NaN
x[505] = NaN
y[580] = NaN

for i in 800:900
	y[i] = NaN
end

#==Generate plot
===============================================================================#
mplot = InspectDR.Multiplot(title="Waveforms with NaN")
mplot.layout[:ncolumns] = 1

plot = add(mplot, InspectDR.transientplot(:lin, title="Transient Data",
	xlabel="Time (s)", ylabels=["Voltage (V)"])
)
	plot.displayNaN = true #Enable hilighting of NaN values (slower)
	wfrm = add(plot, x, y)
	wfrm.line = line(color=blue, width=3)

gplot = display(InspectDR.GtkDisplay(), mplot)

:DONE


text:
## Copyright (c) 2013 Miguel Bazdresch
##
## This file is distributed under the 2-clause BSD License.

# Gaston test framework
#
# There are two kinds of tests:
#   - Tests that should error -- things that Gaston should not allow, such as
#     building a figure with an invalid plotstyle.
#   - Tests that should succeed -- things that should work without producing
#     an error.

using Gaston
using Base.Test

test_success_expected(r::Test.Success) = println("Test passed.")
test_success_expected(r) = error("Success expected; test failed: $(r.expr).")
test_failure_expected(r::Test.Success) = error("Error expected; test failed: $(r.expr)")
test_failure_expected(r) = println("Test passed.")

function run_tests_success_expected()
	set_filename(tempname())
	closeall()
	Test.with_handler(test_success_expected) do
		@test figure() == 1
		@test figure() == 2
		@test figure(4) == 4
		@test closefigure(4) == 4
		@test closefigure() == 2
		@test closeall() == 1
		@test closeall() == 0
		@test plot(1:10) == 1
		@test plot(2,1:10) == 2
		@test plot(4,1:10) == 4
		@test closeall() == 3
		@test begin
			plot(sin(-3:0.01:3),
			"legend", "sine",
			"plotstyle", "lines",
			"color","blue",
			"marker","ecircle",
			"linewidth",2,
			"pointsize",1.1,
			"title","test plot 1",
			"xlabel","x",
			"ylabel","y",
			"box","inside horizontal left top",
			"axis","loglog") == 1
		end
		@test histogram(rand(1000)) == 1
		@test begin
			histogram(randn(1000),
            "bins",100,
            "norm",1,
            "color","blue",
            "linewidth",2,
            "title","test histogram",
            "xlabel","x",
            "ylabel","y",
            "box","inside horizontal left top") == 1
		end
		z = rand(5,6)
		@test imagesc(z,"title","test imagesc 1","xlabel","xx","ylabel","yy") == 1
		@test imagesc(1:5,z,"title","test imagesc 2","xlabel","xx","ylabel","yy") == 1
		@test imagesc(1:5,1:6,z,"title","test imagesc 3","xlabel","xx","ylabel","yy") == 1
		@test surf(rand(10,10)) == 1
		@test surf(2:11,rand(10,10)) == 1
		@test surf(0:9,2:11,rand(10,10)) == 1
	    @test surf(0:9,2:11,(x,y)->x*y) == 1
		@test surf(0:9,2:11,(x,y)->x*y,"title","test","plotstyle","pm3d") == 1
		@test printfigure() == 1
		@test begin
			plot(2,1:10)
			printfigure(2,"png")
			closefigure() == 2
		end
		@test printfigure("eps") == 1
		@test printfigure("pdf") == 1
		@test begin
			set_print_size("640,480")
			printfigure("svg") == 1
		end
		@test printfigure("gif") == 1
		@test plot(1:10,"xrange","[2:3]") == 1
		@test plot(1:10,"xrange","[-1.1:3.4]") == 1
		@test plot(1:10,"xrange","[:3.4]") == 1
		@test plot(1:10,"xrange","[3.4:]") == 1
		@test plot(1:10,"xrange","[3.4:*]") == 1
		@test plot(1:10,"xrange","[*:3.4]") == 1
		closeall()
	end
end

function run_tests_failure_expected()
	closeall()
	Test.with_handler(test_failure_expected) do
		# high-level functions: figure
		@test figure("invalid")
    	@test figure(1.0)
	    @test figure(1:2)
		# high-level functions: closefigure
		@test closefigure(-1)
		@test closefigure("invalid")
		@test closefigure(1.0)
		@test closefigure(1:2)
		# high-level functions: plot and histogram
		for op = (:plot, :histogram)
			@test op("linewidth")
			@test op(0:10,0:11)
			@test op(0:10+im*0:10)
			@test op(0:10,"legend",0)
			@test op(0:10,"plotstyle","invalid")
			@test op(0:10,"marker","invalid")
			@test op(0:10,"marker",0)
			@test op(0:10,"linewidth","b")
			@test op(0:10,"linewidth",im)
			@test op(0:10,"pointsize","b")
			@test op(0:10,"pointsize",im)
			@test op(0:10,"title",0)
			@test op(0:10,"xlabel",0)
			@test op(0:10,"ylabel",0)
			@test op(0:10,"zlabel","z")
			@test op(0:10,"axis","invalid")
		end
		@test plot(1:10,"xrange","2:3")
		@test plot(1:10,"xrange","ab")
	end
end

println("Running Gaston tests.")
run_tests_success_expected()
run_tests_failure_expected()


text:
#
# lgemv.jl -
#
# Lazily Generalized Matrix-Vector mutiplication.
#
#-------------------------------------------------------------------------------
#
# This file is part of LazyAlgebra (https://github.com/emmt/LazyAlgebra.jl)
# released under the MIT "Expat" license.
#
# Copyright (c) 2017-2020 Éric Thiébaut.
#

"""
# Lazily Generalized Matrix-Vector mutiplication

```julia
lgemv([α=1,] [tr='N',] A, x) -> y
```

yields `y = α*op(A)*x` where `op(A)` is `A` if `tr` is `'N'`, `transpose(A)` if
`tr` is `'T'` and `adjoint(A)` if `tr` is `'C'`.  The expression `op(A)*x` is a
matrix-vector multiplication but interpreted by grouping consecutive dimensions
of `A` and `x` as follows:

- if `tr` is `'N'`, the trailing dimensions of `x` must match those of `A` and
  the leading dimensions of `A` are those of the result `y`;

- if `tr` is `'T'` or `'C'`, the leading dimensions of `x` must match those of
  `A` and the trailing dimensions of `A` are those of the result `y`.

The in-place version is called by:

```julia
lgemv!([α=1,] [tr='N',] A, x, [β=0,] y) -> y
```

and overwrites the contents of `y` with `α*op(A)*x + β*y`.  Note that `x` and
`y` must not be aliased.

The multipliers `α` and `β` must be both specified or omitted, they can be any
scalar numbers but are respectively converted to `promote_eltype(A,x)` and
`eltype(y)` which may throw an `InexactError` exception.

See also: [`lgemm`](@ref), [`LinearAlgebra.BLAS.gemv`](@ref),
[`LinearAlgebra.BLAS.gemv!`](@ref).

"""
function lgemv(α::Number,
               trans::Char,
               A::AbstractArray{<:Floats},
               x::AbstractArray{<:Floats})
    return _lgemv(Implementation(Val(:lgemv), α, trans, A, x),
                  α, trans, A, x)
end

function lgemv!(α::Number,
                trans::Char,
                A::AbstractArray{<:Floats},
                x::AbstractArray{<:Floats},
                β::Number,
                y::AbstractArray{<:Floats})
    return _lgemv!(Implementation(Val(:lgemv), α, trans, A, x, β, y),
                   α, trans, A, x, β, y)
end

@doc @doc(lgemv) lgemv!

lgemv(A::AbstractArray, x::AbstractArray) = lgemv(1, 'N', A, x)

lgemv!(A::AbstractArray, x::AbstractArray, y::AbstractArray) =
    lgemv!(1, 'N', A, x, 0, y)

lgemv(trans::Char, A::AbstractArray, x::AbstractArray) =
    lgemv(1, trans, A, x)

lgemv!(trans::Char, A::AbstractArray, x::AbstractArray,
       y::AbstractArray) = lgemv!(1, trans, A, x, 0, y)

lgemv(α::Number, A::AbstractArray, x::AbstractArray) = lgemv(α, 'N', A, x)

lgemv!(α::Number, A::AbstractArray, x::AbstractArray, β::Number,
       y::AbstractArray) = lgemv!(α, 'N', A, x, β, y)

@static if isdefined(LinearAlgebra, :Transpose)
    lgemv(A::LinearAlgebra.Transpose, x::AbstractArray) =
        lgemv(1, 'T', A.parent, x)
    lgemv!(A::LinearAlgebra.Transpose, x::AbstractArray, y::AbstractArray) =
        lgemv!(1, 'T', A.parent, x, 0, y)
    lgemv(α::Number, A::LinearAlgebra.Transpose, x::AbstractArray) =
        lgemv(α, 'T', A.parent, x)
    lgemv!(α::Number, A::LinearAlgebra.Transpose, x::AbstractArray,
           β::Number, y::AbstractArray) = lgemv!(α, 'T', A.parent, x, β, y)
end

@static if isdefined(LinearAlgebra, :Adjoint)
    lgemv(A::LinearAlgebra.Adjoint, x::AbstractArray) =
        lgemv(1, 'C', A.parent, x)
    lgemv!(A::LinearAlgebra.Adjoint, x::AbstractArray, y::AbstractArray) =
        lgemv!(1, 'C', A.parent, x, 0, y)
    lgemv(α::Number, A::LinearAlgebra.Adjoint, x::AbstractArray) =
        lgemv(α, 'C', A.parent, x)
    lgemv!(α::Number, A::LinearAlgebra.Adjoint, x::AbstractArray, β::Number,
           y::AbstractArray) = lgemv!(α, 'C', A.parent, x, β, y)
end

# Best implementations for lgemv and lgemv!

for (atyp, eltyp) in ((Real,   BlasReal),
                      (Number, BlasComplex))
    @eval begin
        function Implementation(::Val{:lgemv},
                                α::$atyp,
                                trans::Char,
                                A::Matrix{T},
                                x::Vector{T}) where {T<:$eltyp}
            return Blas()
        end

        function Implementation(::Val{:lgemv},
                                α::$atyp,
                                trans::Char,
                                A::Matrix{T},
                                x::Vector{T},
                                β::$atyp,
                                y::Vector{T}) where {T<:$eltyp}
            return Blas()
        end

        function Implementation(::Val{:lgemv},
                                α::$atyp,
                                trans::Char,
                                A::AbstractMatrix{T},
                                x::AbstractVector{T}) where {T<:$eltyp}
            return (is_flat_array(A, x) ? Blas() : Basic())
        end

        function Implementation(::Val{:lgemv},
                                α::$atyp,
                                trans::Char,
                                A::AbstractMatrix{T},
                                x::AbstractVector{T},
                                β::$atyp,
                                y::AbstractVector{T}) where {T<:$eltyp}
            return (is_flat_array(A, x, y) ? Blas() : Basic())
        end

        function Implementation(::Val{:lgemv},
                                α::$atyp,
                                trans::Char,
                                A::Array{T},
                                x::Array{T}) where {T<:$eltyp}
            return Blas()
        end

        function Implementation(::Val{:lgemv},
                                α::$atyp,
                                trans::Char,
                                A::Array{T},
                                x::Array{T},
                                β::$atyp,
                                y::Array{T}) where {T<:$eltyp}
            return Blas()
        end

        function Implementation(::Val{:lgemv},
                                α::$atyp,
                                trans::Char,
                                A::AbstractArray{T},
                                x::AbstractArray{T}) where {T<:$eltyp}
            return (is_flat_array(A, x) ? Blas() : Basic())
        end

        function Implementation(::Val{:lgemv},
                                α::$atyp,
                                trans::Char,
                                A::AbstractArray{T},
                                x::AbstractArray{T},
                                β::$atyp,
                                y::AbstractArray{T}) where {T<:$eltyp}
            return (is_flat_array(A, x, y) ? Blas() : Basic())
        end
    end
end

function Implementation(::Val{:lgemv},
                        α::Number,
                        trans::Char,
                        A::AbstractMatrix,
                        x::AbstractVector)
    return Basic()
end

function Implementation(::Val{:lgemv},
                        α::Number,
                        trans::Char,
                        A::AbstractMatrix,
                        x::AbstractVector,
                        β::Number,
                        y::AbstractVector)
    return Basic()
end

function Implementation(::Val{:lgemv},
                        α::Number,
                        trans::Char,
                        A::AbstractArray,
                        x::AbstractArray)
    return (is_flat_array(A, x) ? Linear() : Generic())
end

function Implementation(::Val{:lgemv},
                        α::Number,
                        trans::Char,
                        A::AbstractArray,
                        x::AbstractArray,
                        β::Number,
                        y::AbstractArray)
    return (is_flat_array(A, x, y) ? Linear() : Generic())
end

# BLAS implementations for (generalized) matrices and vectors.
# Linear indexing is assumed (this must have been checked before).

function _lgemv(::Blas,
                α::Number,
                trans::Char,
                A::AbstractArray{T},
                x::AbstractArray{T}) where {T<:BlasFloat}
    nrows, ncols, shape = _lgemv_dims(trans, A, x)
    return _blas_lgemv!(nrows, ncols, convert(T, α), trans, A, x,
                        zero(T), Array{T}(undef, shape))
end

function _lgemv!(::Blas,
                 α::Number,
                 trans::Char,
                 A::AbstractArray{T},
                 x::AbstractArray{T},
                 β::Number,
                 y::AbstractArray{T}) where {T<:BlasFloat}
    nrows, ncols = _lgemv_dims(trans, A, x, y)
    return _blas_lgemv!(nrows, ncols, convert(T, α), trans, A, x,
                        convert(T, β), y)
end

# Julia implementations for (generalized) matrices and vectors.
# Linear indexing is assumed (this must have been checked before).

function _lgemv(::Linear,
                α::Number,
                trans::Char,
                A::AbstractArray{<:Floats},
                x::AbstractArray{<:Floats})
    nrows, ncols, shape = _lgemv_dims(trans, A, x)
    T = _lgemv_type(α, A, x)
    return _linear_lgemv!(nrows, ncols,
                          promote_multiplier(α, A, x), trans, A, x,
                          promote_multiplier(0, T), Array{T}(undef, shape))
end

function _lgemv!(::Linear,
                 α::Number,
                 trans::Char,
                 A::AbstractArray{<:Floats},
                 x::AbstractArray{<:Floats},
                 β::Number,
                 y::AbstractArray{<:Floats})
    nrows, ncols = _lgemv_dims(trans, A, x, y)
    return _linear_lgemv!(nrows, ncols,
                          promote_multiplier(α, A, x), trans, A, x,
                          promote_multiplier(β, y), y)
end

# Basic Julia implementations when vectors and matrices are, respectively, 1D
# and 2D arrays.

function _lgemv(::Basic,
                α::Number,
                trans::Char,
                A::AbstractMatrix{<:Floats},
                x::AbstractVector{<:Floats})
    rows, cols = _lgemv_indices(trans, A, x)
    T = _lgemv_type(α, A, x)
    return _generic_lgemv!(rows, cols,
                           promote_multiplier(α, A, x), trans, A, x,
                           promote_multiplier(0, T),
                           similar(Array{T}, trans == 'N' ? rows : cols))
end

function _lgemv!(::Basic,
                 α::Number,
                 trans::Char,
                 A::AbstractMatrix{<:Floats},
                 x::AbstractVector{<:Floats},
                 β::Number,
                 y::AbstractVector{<:Floats})
    rows, cols = _lgemv_indices(trans, A, x, y)
    return _generic_lgemv!(rows, cols,
                           promote_multiplier(α, A, x), trans, A, x,
                           promote_multiplier(β, y), y)
end

# Generic implementations for any other cases.

function _lgemv(::Generic,
                α::Number,
                trans::Char,
                A::AbstractArray{<:Floats},
                x::AbstractArray{<:Floats})
    rows, cols = _lgemv_indices(trans, A, x)
    T = _lgemv_type(α, A, x)
    return _generic_lgemv!(cartesian_indices(rows), cartesian_indices(cols),
                           promote_multiplier(α, A, x), trans, A, x,
                           promote_multiplier(0, T),
                           similar(Array{T}, trans == 'N' ? rows : cols))
end

function _lgemv!(::Generic,
                 α::Number,
                 trans::Char,
                 A::AbstractArray{<:Floats},
                 x::AbstractArray{<:Floats},
                 β::Number,
                 y::AbstractArray{<:Floats})
    rows, cols = _lgemv_indices(trans, A, x, y)
    return _generic_lgemv!(cartesian_indices(rows), cartesian_indices(cols),
                           promote_multiplier(α, A, x), trans, A, x,
                           promote_multiplier(β, y), y)
end

#
# Call low-level BLAS version.  The differences with LinearAlgebra.BLAS.gemv!
# are that inputs are assumed to be flat arrays (see is_flat_array) and that
# multipliers are automatically converted.
#
for (f, T) in ((:dgemv_, Float64),
               (:sgemv_, Float32),
               (:zgemv_, ComplexF64),
               (:cgemv_, ComplexF32))
    @eval begin
        #
        # FORTRAN prototype:
        #     SUBROUTINE ${pfx}GEMV(TRANS,M,N,ALPHA,A,LDA,X,INCX,BETA,Y,INCY)
        #         ${T} ALPHA,BETA
        #         INTEGER M,N,LDA,INCX,INCY
        #         CHARACTER TRANS
        #         ${T} A(LDA,*),X(*),Y(*)
        #
        # Scalar arguments, α and β, can just be `Number` and integer arguments
        # can just be `Integer` but we want to keep the signature strict
        # because it is a low-level private method.
        #
        function _blas_lgemv!(nrows::Int,
                              ncols::Int,
                              α::($T),
                              trans::Char,
                              A::AbstractArray{$T},
                              x::AbstractArray{$T},
                              β::($T),
                              y::AbstractArray{$T})
            #@static if DEBUG
            #    (length(x) == (trans == 'N' ? ncols : nrows) &&
            #     length(y) == (trans == 'N' ? nrows : ncols) &&
            #     length(A) == nrows*ncols) ||
            #     bad_size("incompatible sizes")
            #end
            ccall((@blasfunc($f), libblas), Cvoid,
                  (Ref{UInt8}, Ref{BlasInt}, Ref{BlasInt}, Ref{$T},
                   Ptr{$T}, Ref{BlasInt}, Ptr{$T}, Ref{BlasInt},
                   Ref{$T}, Ptr{$T}, Ref{BlasInt}),
                  trans, nrows, ncols, α, A, nrows, x, 1, β, y, 1)
            return y
        end
    end
end
#
# Reference (non-BLAS) version for *flat* arrays which can be linearly indexed.
# Loops are ordered assuming the coefficients of A have column-major storage
# order.
#
function _linear_lgemv!(nrows::Int, ncols::Int,
                        α::Floats,
                        trans::Char,
                        A::AbstractArray{Ta},
                        x::AbstractArray{Tx},
                        β::Floats,
                        y::AbstractArray{Ty}) where {Ta<:Floats,
                                                     Tx<:Floats,
                                                     Ty<:Floats}
    if α == 0 || trans == 'N'
        # Form: y := β⋅y
        if β == 0
            @inbounds @simd for k in eachindex(y)
                y[k] = zero(Ty)
            end
        elseif β != 1
            @inbounds @simd for k in eachindex(y)
                y[k] *= β
            end
        end
    end
    if α != 0
        if trans == 'N'
            @inbounds for j in 1:ncols
                temp = α*x[j]
                if temp != zero(temp)
                    off = (j - 1)*nrows
                    @simd for i in 1:nrows
                        y[i] += A[off + i]*temp
                    end
                end
            end
        else
            T = promote_type(Ta, Tx)
            if Ta <: Reals || trans == 'T'
                @inbounds for j in 1:ncols
                    off = (j - 1)*nrows
                    temp = zero(T)
                    @simd for i in 1:nrows
                        temp += A[off + i]*x[i]
                    end
                    y[j] = (β == 0 ? α*temp : α*temp + β*y[j])
                end
            else
                @inbounds for j in 1:ncols
                    off = (j - 1)*nrows
                    temp = zero(T)
                    @simd for i in 1:nrows
                        temp += conj(A[off + i])*x[i]
                    end
                    y[j] = (β == 0 ? α*temp : α*temp + β*y[j])
                end
            end
        end
    end
    return y
end
#
# At the lowest level, the same code can serve for the very general case
# (multi-dimensional Cartesian indices) and the basic case (A is a 2D array
# while x and y are both vectors).
#
# The elements of A are accessed sequentially with one pass through A assuming
# they are stored in colum-major order.
#
function _generic_lgemv!(I, J,
                         α::Floats,
                         trans::Char,
                         A::AbstractArray{Ta},
                         x::AbstractArray{Tx},
                         β::Floats,
                         y::AbstractArray{Ty}) where {Ta<:Floats,
                                                      Tx<:Floats,
                                                      Ty<:Floats}
    if α == 0 || trans == 'N'
        # Form: y := β⋅y
        if β == 0
            @inbounds @simd for k in eachindex(y)
                y[k] = zero(Ty)
            end
        elseif β != 1
            @inbounds @simd for k in eachindex(y)
                y[k] *= β
            end
        end
    end
    if α != 0
        if trans == 'N'
            #
            # Form  y := α*A*x + y.
            #
            @inbounds for j in J
                temp = α*x[j]
                if temp != zero(temp)
                    @simd for i in I
                        y[i] += A[i,j]*temp
                    end
                end
            end
        else
            T = promote_type(Ta,Tx)
            if Ta <: Real || trans == 'T'
                #
                # Form  y := α*A^T*x + y
                #
                @inbounds for j in J
                    temp = zero(T)
                    @simd for i in I
                        temp += A[i,j]*x[i]
                    end
                    y[j] = (β == 0 ? α*temp : α*temp + β*y[j])
                end
            else
                #
                # Form  y := α*A^H*x + y.
                #
                @inbounds for j in J
                    temp = zero(T)
                    @simd for i in I
                        temp += conj(A[i,j])*x[i]
                    end
                    y[j] = (β == 0 ? α*temp : α*temp + β*y[j])
                end
            end
        end
    end
    return y
end
#
# This method yields the type of the elements for the result of lgemv.
#
_lgemv_type(α::Real, A::AbstractArray, x::AbstractArray) =
    promote_eltype(A, x)
_lgemv_type(α::Complex, A::AbstractArray, x::AbstractArray) =
    complex(promote_eltype(A, x))
#
# This method yields the number of rows and columns for lgemv assuming linear
# indexing and check arguments.
#
@inline function _lgemv_dims(trans::Char,
                             A::AbstractMatrix,
                             x::AbstractVector)
    nrows, ncols = size(A, 1), size(A, 2)
    if trans == 'N'
        length(x) == ncols || incompatible_dimensions()
    elseif trans == 'T' || trans == 'C'
        length(x) == nrows || incompatible_dimensions()
    else
        invalid_transpose_character()
    end
    return nrows, ncols, (trans == 'N' ? nrows : ncols)
end
#
# Idem for general "flat" arrays.
#
function _lgemv_dims(trans::Char,
                     A::AbstractArray{<:Any,Na},
                     x::AbstractArray{<:Any,Nx}) where {Na,Nx}
    1 ≤ Nx < Na || incompatible_dimensions()
    @inbounds begin
        Ny = Na - Nx
        if trans == 'N'
            ncols = 1
            for d in 1:Nx
                dim = size(x, d)
                size(A, Ny + d) == dim || incompatible_dimensions()
                ncols *= dim
            end
            shape = ntuple(d -> size(A, d), Ny)
            nrows = prod(shape)
        elseif trans == 'T' || trans == 'C'
            nrows = 1
            for d in 1:Nx
                dim = size(x, d)
                size(A, d) == dim || incompatible_dimensions()
                nrows *= dim
            end
            shape = ntuple(d -> size(A, Nx + d), Ny)
            ncols = prod(shape)
        else
            invalid_transpose_character()
        end
        return nrows, ncols, shape
    end
end
#
# This method yields the number of rows and columns for lgemv! assuming linear
# indexing and check arguments.
#
@inline function _lgemv_dims(trans::Char,
                             A::AbstractMatrix,
                             x::AbstractVector,
                             y::AbstractVector)
    nrows, ncols = size(A)
    if trans == 'N'
        (length(x) == ncols && length(y) == nrows) || incompatible_dimensions()
    elseif trans == 'T' || trans == 'C'
        (length(x) == nrows && length(y) == ncols) || incompatible_dimensions()
    else
        invalid_transpose_character()
    end
    return nrows, ncols
end
#
# Idem for general arrays.
#
function _lgemv_dims(trans::Char,
                     A::AbstractArray{<:Any,Na},
                     x::AbstractArray{<:Any,Nx},
                     y::AbstractArray{<:Any,Ny}) where {Na,Nx,Ny}
    (Na == Nx + Ny && Nx ≥ 1 && Ny ≥ 1) || incompatible_dimensions()
    nrows = ncols = 1
    @inbounds begin
        if trans == 'N'
            for d in 1:Ny
                dim = size(y, d)
                size(A, d) == dim || incompatible_dimensions()
                nrows *= dim
            end
            for d in 1:Nx
                dim = size(x, d)
                size(A, d + Ny) == dim || incompatible_dimensions()
                ncols *= dim
            end
        elseif trans == 'T' || trans == 'C'
            for d in 1:Nx
                dim = size(x, d)
                size(A, d) == dim || incompatible_dimensions()
                nrows *= dim
            end
            for d in 1:Ny
                dim = size(y, d)
                size(A, d + Nx) == dim || incompatible_dimensions()
                ncols *= dim
            end
        else
            invalid_transpose_character()
        end
    end
    return (nrows, ncols)
end
#
# Build tuples rows and cols of index intervals to access A[i,j] in lgemv and
# check arguments.
#
@inline function _lgemv_indices(trans::Char,
                                A::AbstractMatrix,
                                x::AbstractVector)
    rows, cols = axes(A)
    if trans == 'N'
        axes(x) == (cols,) || incompatible_dimensions()
    elseif trans == 'T' || trans == 'C'
        axes(x) == (rows,) || incompatible_dimensions()
    else
        invalid_transpose_character()
    end
    return rows, cols
end
#
# Idem for general arrays.
#
@inline function _lgemv_indices(trans::Char,
                                A::AbstractArray{<:Any,Na},
                                x::AbstractArray{<:Any,Nx}) where {Na,Nx}
    1 ≤ Nx < Na || incompatible_dimensions()
    @inbounds begin
        Ny = Na - Nx
        if trans == 'N'
            rows = ntuple(d -> axes(A, d), Ny)
            cols = ntuple(d -> axes(A, d + Ny), Nx)
            for d in 1:Nx
                axes(x, d) == cols[d] || incompatible_dimensions()
            end
        elseif trans == 'T' || trans == 'C'
            rows = ntuple(d -> axes(A, d), Nx)
            cols = ntuple(d -> axes(A, d + Nx), Ny)
            for d in 1:Nx
                axes(x, d) == rows[d] || incompatible_dimensions()
            end
        else
            invalid_transpose_character()
        end
        return rows, cols
    end
end
#
# Build tuples rows and cols of index intervals to access A[i,j] in lgemv! and
# check arguments.
#
@inline function _lgemv_indices(trans::Char,
                                A::AbstractMatrix,
                                x::AbstractVector,
                                y::AbstractVector)
    rows, cols = axes(A)
    if trans == 'N'
        (axes(x) == (cols,) && axes(y) == (rows,)) || incompatible_dimensions()
    elseif trans == 'T' || trans == 'C'
        (axes(x) == (rows,) && axes(y) == (cols,)) || incompatible_dimensions()
    else
        invalid_transpose_character()
    end
    return rows, cols
end
#
# Idem for general arrays.
#
@inline function _lgemv_indices(trans::Char,
                                A::AbstractArray{<:Any,Na},
                                x::AbstractArray{<:Any,Nx},
                                y::AbstractArray{<:Any,Ny}) where {Na,Nx,Ny}
    (Na == Nx + Ny && Nx ≥ 1 && Ny ≥ 1) || incompatible_dimensions()
    @inbounds begin
        if trans == 'N'
            rows = ntuple(d -> axes(A, d), Ny)
            cols = ntuple(d -> axes(A, d + Ny), Nx)
            for d in 1:Nx
                axes(x, d) == cols[d] || incompatible_dimensions()
            end
            for d in 1:Ny
                axes(y, d) == rows[d] || incompatible_dimensions()
            end
        elseif trans == 'T' || trans == 'C'
            rows = ntuple(d -> axes(A, d), Nx)
            cols = ntuple(d -> axes(A, d + Nx), Ny)
            for d in 1:Nx
                axes(x, d) == rows[d] || incompatible_dimensions()
            end
            for d in 1:Ny
                axes(y, d) == cols[d] || incompatible_dimensions()
            end
        else
            invalid_transpose_character()
        end
        return rows, cols
    end
end


text:
# Variables #
# ========================== #
function JuMP.build_variable(_error::Function, variable::JuMP.ScalarVariable, set::DecisionSet)
    return VariableConstrainedOnCreation(variable, decision_set(variable, set))
end

function JuMP.build_variable(_error::Function, variables::Vector{<:JuMP.ScalarVariable}, set::DecisionSet)
    return VariablesConstrainedOnCreation(variables, decision_set(variables, set))
end

function JuMP.build_variable(_error::Function, variable::JuMP.ScalarVariable, set::KnownSet)
    return VariableConstrainedOnCreation(variable, decision_set(variable, set))
end

function JuMP.build_variable(_error::Function, variables::Vector{<:JuMP.ScalarVariable}, set::KnownSet)
    return VariablesConstrainedOnCreation(variables, decision_set(variables, set))
end

function JuMP.add_variable(model::Model, variable::VariableConstrainedOnCreation{<:SingleDecisionSet}, name::String)
    decisions = get_decisions(model, variable.set.stage)
    if decisions isa IgnoreDecisions
        # Create a regular JuMP variable if decisions are not handled
        return JuMP.add_variable(model, variable.scalar_variable, name)
    end
    var_index, con_index = MOI.add_constrained_variable(backend(model), variable.set)
    model.ext[:stage_map][var_index] = variable.set.stage
    # Map to model decisions after indices are known
    if !has_decision(decisions, var_index)
        # Store decision if is seen for the first time
        set_decision!(decisions, var_index, variable.set.decision)
    else
        # Reuse if decision has been created already
        MOI.set(backend(model), MOI.ConstraintSet(), con_index, reuse(variable.set, decision(decisions, var_index)))
    end
    # Add any given decision constraints
    if state(variable.set.decision) != Known
        _moi_constrain_decision(backend(model), var_index, variable.scalar_variable.info, variable.set)
    end
    # Finally, set any given name
    if !isempty(name)
        MOI.set(backend(model), MOI.VariableName(), var_index, name)
    end
    # Return created decision as DecisionRef
    return DecisionRef(model, var_index)
end

function JuMP.add_variable(model::Model, variable::VariablesConstrainedOnCreation{<:MultipleDecisionSet}, names)
    decisions = get_decisions(model, variable.set.stage)
    if decisions isa IgnoreDecisions
        # Create regular JuMP variables if decisions are not handled
        var_refs = map(zip(variable.scalar_variables, names)) do (scalar_variable, name)
            JuMP.add_variable(model, scalar_variable, name)
        end
        return reshape_vector(var_refs, variable.shape)
    end
    if variable.set.constraint isa NoSpecifiedConstraint
        # Sufficient to use single decision bridges if decisions
        # are not initially constrained.
        refs = map(zip(variable.scalar_variables, variable.set.decisions, names)) do (scalar_variable, decision, name)
            set = SingleDecisionSet(variable.set.stage, decision, variable.set.constraint, variable.set.is_recourse)
            add_variable(model, VariableConstrainedOnCreation(scalar_variable, set), name)
        end
        return reshape_vector(refs, variable.shape)
    else
        var_indices, con_index = MOI.add_constrained_variables(backend(model), variable.set)
        # Map to model decisions after indices are known
        seen_decisions = Vector{Decision{Float64}}()
        for (i, var_index) in enumerate(var_indices)
            model.ext[:stage_map][var_index] = variable.set.stage
            if !has_decision(decisions, var_index)
                # Store decision if is seen for the first time
                set_decision!(decisions, var_index, variable.set.decisions[i])
            else
                # Reuse if decision has been created already
                push!(seen_decisions, decision(decisions, var_index))
            end
        end
        if !isempty(seen_decisions)
            # Sanity check
            length(seen_decisions) == length(variable.scalar_variables) || error("Inconsistency in number of seen decisions and created variables.")
            # Update decision set for reuse
            MOI.set(backend(model), MOI.ConstraintSet(), con_index, reuse(variable.set, seen_decisions))
        end
        # Add any given decision constraints
        for (index, scalar_variable, decision) in zip(var_indices, variable.scalar_variables, variable.set.decisions)
            if decision.state != Known
                _moi_constrain_decision(backend(model), index, scalar_variable.info, variable.set)
            end
        end
        # Finally, set any given names
        for (var_index, name) in zip(var_indices, JuMP.vectorize(names, variable.shape))
            if !isempty(name)
                MOI.set(backend(model), MOI.VariableName(), var_index, name)
            end
        end
        # Return created decisions as DecisionRef
        refs = map(var_indices) do index
            DecisionRef(model, index)
        end
        return reshape_vector(refs, variable.shape)
    end
end

# Containers #
# ========================== #
struct DecisionDenseAxisArray{A <: DenseAxisArray, S <: Union{DecisionSet, KnownSet}}
    array::A
    set::S

    function DecisionDenseAxisArray(array::DenseAxisArray, set::Union{DecisionSet, KnownSet})
        A = typeof(array)
        S = typeof(set)
        return new{A, S}(array, set)
    end
end

struct DecisionSparseAxisArray{A <: SparseAxisArray, S <: Union{DecisionSet, KnownSet}}
    array::A
    set::S

    function DecisionSparseAxisArray(array::SparseAxisArray, set::Union{DecisionSet, KnownSet})
        A = typeof(array)
        S = typeof(set)
        return new{A, S}(array, set)
    end
end

function JuMP.build_variable(_error::Function, variables::DenseAxisArray{<:JuMP.ScalarVariable}, set::Union{DecisionSet,KnownSet})
    return DecisionDenseAxisArray(variables, set)
end

function JuMP.build_variable(_error::Function, variables::SparseAxisArray{<:JuMP.ScalarVariable}, set::Union{DecisionSet,KnownSet})
    return DecisionSparseAxisArray(variables, set)
end

function JuMP.add_variable(model::Model, variable::DecisionDenseAxisArray, names::DenseAxisArray{String})
    array = variable.array
    refs = DenseAxisArray{DecisionRef}(undef, array.axes...)
    for idx in eachindex(array)
        var = array[idx]
        refs[idx] = JuMP.add_variable(model,
                                      VariableConstrainedOnCreation(var, decision_set(var, variable.set)),
                                      names[idx])
    end
    return refs
end

function JuMP.add_variable(model::Model,
                           variable::DecisionSparseAxisArray{SparseAxisArray{T,N,K}},
                           names::SparseAxisArray{String}) where {T, N, K}
    refs = SparseAxisArray(Dict{K,DecisionRef}())
    for (idx, var) in variable.array.data
        refs[idx] = JuMP.add_variable(model,
                                      VariableConstrainedOnCreation(var, decision_set(var, variable.set)),
                                      names[idx])
    end
    return refs
end

# Constraints #
# ========================== #
function JuMP.build_constraint(_error::Function, aff::Union{DecisionAffExpr, DecisionQuadExpr}, set::MOI.AbstractScalarSet)
    offset = constant(aff.variables)
    JuMP.add_to_expression!(aff.variables, -offset)
    shifted_set = MOIU.shift_constant(set, -offset)
    return JuMP.ScalarConstraint(aff, shifted_set)
end

function JuMP.build_constraint(_error::Function, aff::Union{DecisionAffExpr, DecisionQuadExpr}, lb, ub)
    JuMP.build_constraint(_error, aff, MOI.Interval(lb, ub))
end

# Helper function #
# ========================== #
function _moi_constrain_decision(backend::MOI.ModelLike, index, info, set::Union{SingleDecisionSet, MultipleDecisionSet})
    # We don't call the _moi* versions (e.g., _moi_set_lower_bound) because they
    # have extra checks that are not necessary for newly created variables.
    nothing_added = true
    if info.has_lb
        MOI.add_constraint(backend, SingleDecision(index),
                           MOI.GreaterThan{Float64}(info.lower_bound))
        nothing_added &= false
    end
    if info.has_ub
        MOI.add_constraint(backend, SingleDecision(index),
                           MOI.LessThan{Float64}(info.upper_bound))
        nothing_added &= false
    end
    if info.has_fix
        MOI.add_constraint(backend, SingleDecision(index),
                           MOI.EqualTo{Float64}(info.fixed_value))
        nothing_added &= false
    end
    if info.binary
        MOI.add_constraint(backend, SingleDecision(index),
                           MOI.ZeroOne())
        nothing_added &= false
    end
    if info.integer
        MOI.add_constraint(backend, SingleDecision(index), MOI.Integer())
        nothing_added &= false
    end
    if info.has_start
        MOI.set(backend, MOI.VariablePrimalStart(), index,
                Float64(info.start))
        nothing_added &= false
    end
end