text:
const Initialized = Union{Type{Val{true}}, Type{Val{false}}}

struct OnDiskRepresentation{Offsets,JLTypes,H5Types} end

struct UnknownType{T}
    name::T
    parameters::Vector{Any}

    UnknownType{T}(name) where T = new(name)
    UnknownType{T}(name, parameters) where T = new(name, parameters)
end
UnknownType(name) = UnknownType{typeof(name)}(name)
UnknownType(name, parameters) = UnknownType{typeof(name)}(name, parameters)

struct Vlen{T}
    size::UInt32
    id::GlobalHeapID
end
sizeof{T<:Vlen}(::Type{T}) = 4 + sizeof(GlobalHeapID)

# Initial ODR for DataType
const DataTypeODR = OnDiskRepresentation{(0, sizeof(Vlen{String})),Tuple{String,Vector{Any}},Tuple{Vlen{String},Vlen{RelOffset}}}

const NULL_COMMITTED_DATATYPE = CommittedDatatype(RelOffset(0), 0)
# Look up the corresponding committed datatype for a given type
macro lookup_committed(f, T)
    quote
        cdt = get($(esc(f)).jlh5type, $(esc(T)), nothing)
        cdt !== nothing && return cdt::CommittedDatatype
    end
end

function track_weakref!(f::JLDFile, header_offset::RelOffset, v::ANY)
    header_offset !== NULL_REFERENCE && (f.jloffset[header_offset] = WeakRef(v))
    nothing
end

writeas(T::Type) = T

## Generic machinery

# Carries the type and on-disk representation of data to be read from
# the disk
ReadRepresentation(T, S) = ReadRepresentation{T,S}()
ReadRepresentation(T) = ReadRepresentation{T,T}()
sizeof{T,S}(::ReadRepresentation{T,S}) = sizeof(S)

# Determines whether a specific field type should be saved in the file
@noinline function hasfielddata(T::ANY)
    T === Union{} && return false
    !isleaftype(T) && return true
    T = T::DataType
    (T.mutable || T <: Type) && return true
    hasdata(T)
end

# Determines whether a specific type has fields that should be saved in the file
function hasdata(T::DataType)
    isempty(T.types) && T.size != 0 && return true
    for ty in T.types
        hasfielddata(writeas(ty)) && return true
    end
    false
end

# Gets the size of an on-disk representation
Base.@pure function sizeof{Offsets,JLTypes,H5Types}(::OnDiskRepresentation{Offsets,JLTypes,H5Types})
    Offsets[end]+sizeof(H5Types.parameters[end])
end

# Determines whether a type will have the same layout on disk as in memory
function samelayout(T::DataType)
    isempty(T.types) && return true
    offset = 0
    for i = 1:length(T.types)
        offset != fieldoffset(T, i) && return false
        ty = T.types[i]
        ty !== writeas(ty) && return false
        !samelayout(ty) && return false
        offset += ty.size
    end
    return offset == T.size
end

fieldnames{T<:Tuple}(x::Type{T}) = [Symbol(x) for x = 1:length(x.types)]
fieldnames(x::ANY) = Base.fieldnames(x)

# fieldodr gives the on-disk representation of a field of a given type,
# which is either always initialized (initialized=true) or potentially
# uninitialized (initialized=false)
@generated function fieldodr{T}(::Type{T}, initialized::Bool)
    if isleaftype(T)
        if !hasfielddata(T)
            # A ghost type, so no need to store at all
            return nothing
        elseif isa(T, DataType)
            if isbits(T)
                return :(odr(T))
            elseif !T.mutable
                return :(initialized ? odr(T) : RelOffset)
            end
        end
    end
    RelOffset
end

# h5fieldtype is fieldodr's HDF5 companion. It should give the HDF5
# datatype reflecting the on-disk representation.
@generated function h5fieldtype{T}(f::JLDFile, writeas::Type{T}, readas::Type,
                                   initialized::Initialized)
    if isleaftype(T)
        if !hasfielddata(T)
            return nothing
        elseif isbits(T) || (isa(initialized, Type{Type{Val{true}}}) && !T.mutable)
            return quote
                @lookup_committed f T
                $(if isempty(T.types)
                    # bitstype
                    :(return commit(f, OpaqueDatatype(sizeof(T)), T, readas))
                else
                    # Compound type
                    :(return commit_compound(f, fieldnames(T), T, readas))
                end)
            end
        end
    end
    ReferenceDatatype()
end

# Select an ODR, incorporating custom serialization only if the types do not
# match
CustomSerialization{WrittenAs}(::Type{WrittenAs}, ::Type{WrittenAs}, odr) = odr
CustomSerialization{WrittenAs,ReadAs}(::Type{WrittenAs}, ::Type{ReadAs}, odr) =
    CustomSerialization{WrittenAs,odr}

# objodr gives the on-disk representation of a given object. This is
# almost always the on-disk representation of the type. The only
# exception is strings, where the length is encoded in the datatype in
# HDF5, but in the object in Julia.
@inline function objodr(x)
    writtenas = writeas(typeof(x))
    _odr(writtenas, typeof(x), odr(writtenas))
end
_odr(writtenas::DataType, readas::DataType, odr) =
    CustomSerialization(writtenas, readas, odr)

# h5type is objodr's HDF5 companion. It should give the HDF5 datatype
# reflecting the on-disk representation
#
# Performance note: this should be inferrable.
function h5type(f::JLDFile, writtenas, x)
    check_writtenas_type(writtenas)
    T = typeof(x)
    @lookup_committed f T
    if !hasdata(writtenas)
        commit(f, OpaqueDatatype(1), writtenas, T, WrittenAttribute(f, :empty, UInt8(1)))
    elseif isempty(writtenas.types) # bitstype
        commit(f, OpaqueDatatype(writtenas.size), writtenas, T)
    else
        commit_compound(f, fieldnames(writtenas), writtenas, T)
    end
end
check_writtenas_type(::DataType) = nothing
check_writtenas_type(::Any) = throw(ArgumentError("writeas(leaftype) must return a leaf type"))
h5type(f::JLDFile, x) = h5type(f, writeas(typeof(x)), x)

# Make a compound datatype from a set of names and types
function commit_compound(f::JLDFile, names::AbstractVector{Symbol},
                         writtenas::DataType, readas::Type)
    types = writtenas.types
    offsets = Int[]
    h5names = Symbol[]
    members = H5Datatype[]
    fieldtypes = RelOffset[]
    hasfieldtype = false

    offset = 0
    for i = 1:length(types)
        fieldty = types[i]
        fieldwrittenas = writeas(fieldty)
        !hasfielddata(fieldwrittenas) && continue
        dtype = h5fieldtype(f, fieldwrittenas, fieldty, Val{i <= writtenas.ninitialized})
        dtype === nothing && continue
        push!(h5names, names[i])
        if isa(dtype, CommittedDatatype)
            # HDF5 cannot store relationships among committed
            # datatypes. We store these separately in an attribute.
            push!(fieldtypes, dtype.header_offset)
            dtype = f.datatypes[dtype.index]
            hasfieldtype = true
        else
            push!(fieldtypes, NULL_REFERENCE)
        end
        push!(members, dtype)
        push!(offsets, offset)
        offset += dtype.size::UInt32
    end

    @assert offset != 0
    compound = CompoundDatatype(offset, h5names, offsets, members)
    if hasfieldtype
        fieldtypeattr = WrittenAttribute(:field_datatypes,
                                         WriteDataspace(f, fieldtypes, DataType),
                                         ReferenceDatatype(),
                                         fieldtypes)
        commit(f, compound, writtenas, readas, fieldtypeattr)::CommittedDatatype
    else
        commit(f, compound, writtenas, readas)::CommittedDatatype
    end
end

# Write an HDF5 datatype to the file
function commit(f::JLDFile, dtype::H5Datatype, writeas::DataType, readas::DataType,
                attributes::WrittenAttribute...)
    io = f.io

    # This needs to be written this way or type inference gets unhappy...
    # Also needs to happen here so that we write the DataType type
    # before we try to find where this type will be written
    typeattr = WrittenAttribute(
        :julia_type, WriteDataspace(f, DataType, odr(DataType)), h5type(f, DataType, DataType), readas)

    offset = f.end_of_data

    seek(io, offset)
    id = length(f.datatypes)+1
    h5o = h5offset(f, offset)
    cdt = CommittedDatatype(h5o, id)
    f.datatype_locations[h5o] = cdt
    f.jlh5type[readas] = cdt
    push!(f.datatypes, dtype)

    if writeas !== readas
        wrtypeattr = WrittenAttribute(:written_type,
                                      WriteDataspace(f, DataType, odr(DataType)),
                                      h5type(f, DataType, DataType), writeas)
        f.h5jltype[cdt] = ReadRepresentation(readas, CustomSerialization{writeas, odr(writeas)})
        commit(f, dtype, tuple(typeattr, wrtypeattr, attributes...))
    else
        f.h5jltype[cdt] = ReadRepresentation(writeas, odr(writeas))
        commit(f, dtype, tuple(typeattr, attributes...))
    end

    cdt::CommittedDatatype
end

function read_field_datatypes(f::JLDFile, attrs::Vector{ReadAttribute})
    for attr in attrs
        if attr.name == :field_datatypes
            return read_attr_data(f, attr, ReferenceDatatype(),
                                  ReadRepresentation(RelOffset))
        end
    end
    RelOffset[]
end

function check_empty(attrs::Vector{ReadAttribute})
    for attr in attrs
        if attr.name == :empty
            return true
        end
    end
    false
end

# jltype is the inverse of h5type, providing a ReadRepresentation for an
# H5Datatype. We handle committed datatypes here, and other datatypes below.
function jltype(f::JLDFile, cdt::CommittedDatatype)
    haskey(f.h5jltype, cdt) && return f.h5jltype[cdt]::ReadRepresentation
    dt, attrs = read_committed_datatype(f, cdt)

    julia_type_attr = nothing
    written_type_attr = nothing
    for attr in attrs
        if attr.name == :julia_type
            julia_type_attr = attr
        elseif attr.name == :written_type
            written_type_attr = attr
        end
    end

    if isa(julia_type_attr, Void)
        throw(InvalidDataException())
    end
    julia_type_attr = julia_type_attr::ReadAttribute

    # If type of datatype is this datatype, then this is the committed
    # datatype that describes a datatype
    if h5offset(f, julia_type_attr.datatype_offset) == cdt.header_offset
        # Verify that the datatype matches our expectations
        if dt != H5TYPE_DATATYPE
            error("""The HDF5 datatype representing a Julia datatype does not match
                     the expectations of this version of JLD.

                     You may need to update JLD to read this file.""")
        end
        f.jlh5type[DataType] = cdt
        f.datatypes[cdt.index] = dt
        return (f.h5jltype[cdt] = ReadRepresentation(DataType, DataTypeODR()))
    end

    datatype = read_attr_data(f, julia_type_attr)
    if written_type_attr !== nothing
        # Custom serialization
        readas = datatype
        datatype = read_attr_data(f, written_type_attr)
        if isa(readas, UnknownType)
            warn("custom serialization of ", typestring(readas),
                 " encountered, but the type does not exist in the workspace; the data will be read unconverted")
            rr = (constructrr(f, datatype, dt, attrs)::Tuple{ReadRepresentation,Bool})[1]
            canonical = false
        else
            rr, canonical = constructrr(f, datatype, dt, attrs)::Tuple{ReadRepresentation,Bool}
            rrty = typeof(rr)
            rr = ReadRepresentation{readas, CustomSerialization{rrty.parameters[1], rrty.parameters[2]}}()
            canonical = canonical && writeas(readas) === datatype
        end
    else
        rr, canonical = constructrr(f, datatype, dt, attrs)::Tuple{ReadRepresentation,Bool}
    end

    canonical && (f.jlh5type[datatype] = cdt)
    f.datatypes[cdt.index] = dt
    f.h5jltype[cdt] = rr
end

# Constructs a ReadRepresentation for a given opaque (bitstype) type
function constructrr(::JLDFile, T::DataType, dt::BasicDatatype, attrs::Vector{ReadAttribute})
    dt.class == DT_OPAQUE || throw(UnsupportedFeatureException())
    if sizeof(T) == dt.size
        (ReadRepresentation(T), true)
    else
        empty = check_empty(attrs)
        if empty
            if !hasdata(T)
                (ReadRepresentation(T, nothing), true)
            else
                warn("$T has $(sizeof(T)*8) bytes, but written type was empty; reconstructing")
                reconstruct_bitstype(T.name.name, dt.size, empty)
            end
        else
            if isempty(T.types)
                warn("primitive type $T has $(sizeof(T)*8) bits, but written type has $(dt.size*8) bits; reconstructing")
            else
                warn("$T is a non-primitive type, but written type is a primitive type with $(dt.size*8) bits; reconstructing")
            end
            reconstruct_bitstype(T.name.name, dt.size, empty)
        end
    end
end

function constructrr(f::JLDFile, T::UnionAll, dt::CompoundDatatype,
                     attrs::Vector{ReadAttribute},
                     hard_failure::Bool=false)
    warn("read type $T is not a leaf type in workspace; reconstructing")
    return reconstruct_compound(f, string(T), dt, read_field_datatypes(f, attrs))
end

"""
constructrr(f::JLDFile, T::DataType, dt::CompoundType, attrs::Vector{ReadAttribute},
            hard_failure::Bool=false)

Constructs a ReadRepresentation for a given type. This is the generic method for all
types not specially handled below.

If hard_failure is true, then throw a TypeMappingException instead of attempting
reconstruction. This helps in cases where we can't know if reconstructed parametric types
will have a matching memory layout without first inspecting the memory layout.
"""
struct TypeMappingException <: Exception; end
function constructrr(f::JLDFile, T::DataType, dt::CompoundDatatype,
                     attrs::Vector{ReadAttribute},
                     hard_failure::Bool=false)
    field_datatypes = read_field_datatypes(f, attrs)

    # If read type is not a leaf type, reconstruct
    if !isleaftype(T)
        warn("read type $T is not a leaf type in workspace; reconstructing")
        return reconstruct_compound(f, string(T), dt, field_datatypes)
    end

    # Map names in dt to their indices
    dtnames = Dict{Symbol,Int}()
    for i = 1:length(dt.names)
        dtnames[dt.names[i]] = i
    end
    mapped = falses(length(dt.names))

    offsets = Vector{Int}(length(T.types))
    types = Vector{Any}(length(T.types))
    odrs = Vector{Any}(length(T.types))
    fn = fieldnames(T)
    samelayout = isbits(T) && sizeof(T) == dt.size
    dtindex = 0
    for i = 1:length(T.types)
        wstype = T.types[i]
        writtenas = writeas(wstype)
        if !hasfielddata(writtenas)
            types[i] = wstype
            odrs[i] = wstype === writtenas ? nothing : CustomSerialization{writtenas,nothing}
            offsets[i] = dtindex == 0 ? 0 : (dt.offsets[dtindex] + dt.members[dtindex].size::UInt32)
        else
            if !haskey(dtnames, fn[i])
                hard_failure && throw(TypeMappingException())
                warn("saved type ", T, " is missing field ", fn[i], " in workspace type; reconstructing")
                return reconstruct_compound(f, string(T), dt, field_datatypes)
            end

            dtindex = dtnames[fn[i]]
            if !isempty(field_datatypes) && (ref = field_datatypes[dtindex]) != NULL_REFERENCE
                dtrr = jltype(f, f.datatype_locations[ref])
            else
                dtrr = jltype(f, dt.members[dtindex])
            end

            readtype, odrtype = typeof(dtrr).parameters

            if typeintersect(readtype, wstype) === Union{} &&
               !method_exists(convert, Tuple{Type{wstype}, readtype})
                # Saved type does not match type in workspace and no
                # convert method exists, so we definitely need to reconstruct.
                hard_failure && throw(TypeMappingException())
                warn("saved type ", T, " has field ", fn[i], "::", readtype,
                     ", but workspace type has field ", fn[i], "::", wstype,
                     ", and no applicable convert method exists; reconstructing")
                return reconstruct_compound(f, string(T), dt, field_datatypes)
            end

            types[i] = readtype
            odrs[i] = odrtype
            offsets[i] = dt.offsets[dtindex]
            samelayout = samelayout && offsets[i] == fieldoffset(T, i) && types[i] === wstype

            mapped[dtindex] = true
        end
    end

    if !all(mapped)
        warn("the following fields are present in type ", T,
             " saved in the file but not present in the type the workspace:\n\n",
             join(dt.names[.!mapped], "\n"),
             "\n\nData in these fields will not be accessible")
    end

    if samelayout
        (ReadRepresentation(T), true)
    else
        wodr = odr(T)
        offsets = (offsets...)
        if isa(wodr, OnDiskRepresentation)
            odr_offsets = typeof(wodr).parameters[1]
            odr_types = typeof(wodr).parameters[2].parameters
            odr_h5types = typeof(wodr).parameters[3].parameters
            tequal = length(odr_types) == length(types)
            if tequal
                for i = 1:length(types)
                    if !(odr_types[i] <: types[i]) || odr_h5types[i] != odrs[i]
                        tequal = false
                        break
                    end
                end
                if tequal && odr_types == offsets
                    return (ReadRepresentation(T, wodr), true)
                end
            end
        end
        return (ReadRepresentation(T, OnDiskRepresentation{offsets, Tuple{types...}, Tuple{odrs...}}()), false)
    end
end

# h5convert! stores the HDF5 representation of Julia data to a pointer. This
# method handles types with no padding or references where this is just a simple
# store
h5convert!{T}(out::Ptr, ::Type{T}, ::JLDFile, x, ::JLDWriteSession) =
    (unsafe_store!(convert(Ptr{T}, out), x); nothing)

# We pack types that have padding using a staged h5convert! method
@generated function h5convert!{Offsets,Types,H5Types}(out::Ptr, ::OnDiskRepresentation{Offsets,Types,H5Types}, file::JLDFile, x, wsession::JLDWriteSession)
    T = x
    types = Types.parameters
    members = H5Types.parameters

    getindex_fn = isa(T, Type{T} where T<:Tuple) ? (:getindex) : (:getfield)
    ex = Expr(:block)
    args = ex.args
    for i = 1:length(Offsets)
        member = members[i]
        isa(member, Void) && continue

        offset = Offsets[i]
        conv = :(h5convert!(out+$offset, $(member), file, convert($(types[i]), $getindex_fn(x, $i)), wsession))
        if i > T.ninitialized && (!isleaftype(x.types[i]) || !isbits(x.types[i]))
            push!(args, quote
                if !isdefined(x, $i)
                    h5convert_uninitialized!(out+$offset, $(member))
                else
                    $conv
                end
            end)
        else
            push!(args, conv)
        end
    end
    push!(args, nothing)
    ex
end

# jl_canbeuninitialized specifies whether a given ODR could be uninitialized,
# which determines whether or not we'll try to call jlconvert_isinitialized to
# determine whether or not it is actually defined when reconstructing types.
jlconvert_canbeuninitialized(::Any) = false

# jlconvert converts data from a pointer into a Julia object. This method
# handles types where this is just a simple load
@inline jlconvert{T}(::ReadRepresentation{T,T}, ::JLDFile, ptr::Ptr,
                     ::RelOffset) =
    unsafe_load(convert(Ptr{T}, ptr))

# When fields are undefined in the file but can't be in the workspace, we need
# to throw exceptions to prevent errors on null pointer loads
struct UndefinedFieldException
    ty::DataType
    fieldname::Symbol
end
Base.showerror(io::IO, x::UndefinedFieldException) =
    print(io, "field \"", x.fieldname, "\" of type ", x.ty,
          " must be defined in the current workspace, but was undefined in the file")

## Custom serialization

# wconvert and rconvert do type conversion before reading and writing,
# respectively. These fall back to convert.
wconvert(T, x) = convert(T, x)
rconvert(T, x) = convert(T, x)

sizeof{T,ODR}(::Type{CustomSerialization{T,ODR}}) = sizeof(ODR)

# Usually we want to convert the object and then write it.
@inline h5convert!{T,ODR}(out::Ptr, ::Type{CustomSerialization{T,ODR}}, f::JLDFile,
                          x, wsession::JLDWriteSession) =
    h5convert!(out, ODR, f, wconvert(T, x)::T, wsession)

# When writing as a reference, we don't want to convert the object first. That
# should happen automatically after write_dataset is called so that the written
# object gets the right written_type attribute.
@inline h5convert!{T}(out::Ptr, odr::Type{CustomSerialization{T,RelOffset}},
                      f::JLDFile, x, wsession::JLDWriteSession) =
    h5convert!(out, RelOffset, f, x, wsession)

# When writing as a reference to something that's being custom-serialized as an
# array, we have to convert the object first.
@inline h5convert!{T<:Array}(out::Ptr, odr::Type{CustomSerialization{T,RelOffset}},
                      f::JLDFile, x, wsession::JLDWriteSession) =
    h5convert!(out, RelOffset, f, wconvert(T, x)::T, wsession)

h5convert_uninitialized!{T,ODR}(out::Ptr, odr::Type{CustomSerialization{T,ODR}}) =
    h5convert_uninitialized!(out, ODR)

jlconvert_canbeuninitialized{T,S,ODR}(::ReadRepresentation{T,CustomSerialization{S,ODR}}) =
    jlconvert_canbeuninitialized(ODR)
jlconvert_isinitialized{T,S,ODR}(::ReadRepresentation{T,CustomSerialization{S,ODR}}, ptr::Ptr) =
    jlconvert_isinitialized(ReadRepresentation{S,ODR}(), ptr)
jlconvert{T,S,ODR}(::ReadRepresentation{T,CustomSerialization{S,ODR}},
                   f::JLDFile, ptr::Ptr, header_offset::RelOffset) =
    rconvert(T, jlconvert(ReadRepresentation{S,ODR}(), f, ptr, header_offset))::T

## Primitive datatypes
# These get special handling only in that they have different HDF5 type
# representations than ordinary opaque types

# This construction prevents these methods from getting called on type unions
const SignedTypes        = Union{Type{Int8}, Type{Int16}, Type{Int32}, Type{Int64}, Type{Int128}}
const UnsignedTypes      = Union{Type{UInt8}, Type{UInt16}, Type{UInt32}, Type{UInt64}, Type{UInt128}}
const FloatTypes         = Union{Type{Float16}, Type{Float32}, Type{Float64}}
const PrimitiveTypeTypes = Union{SignedTypes, UnsignedTypes, FloatTypes}
const PrimitiveTypes     = Union{Int8, Int16, Int32, Int64, Int128, UInt8, UInt16, UInt32,
                                 UInt64, UInt128, Float16, Float32, Float64}

for T in Base.uniontypes(SignedTypes)
    @eval h5fieldtype(::JLDFile, ::$T, ::$T, ::Initialized) =
        FixedPointDatatype(sizeof($(T.parameters[1])), true)
end
for T in Base.uniontypes(UnsignedTypes)
    @eval h5fieldtype(::JLDFile, ::$T, ::$T, ::Initialized) =
        FixedPointDatatype(sizeof($(T.parameters[1])), false)
end

function jltype(f::JLDFile, dt::FixedPointDatatype)
    signed = dt.bitfield1 == 0x08 ? true : dt.bitfield1 == 0x00 ? false : throw(UnsupportedFeatureException())
    ((dt.bitfield2 == 0x00) & (dt.bitfield3 == 0x00) & (dt.bitoffset == 0) & (dt.bitprecision == dt.size*8)) ||
        throw(UnsupportedFeatureException())
    if dt.size == 8
        return signed ? ReadRepresentation(Int64) : ReadRepresentation(UInt64)
    elseif dt.size == 1
        return signed ? ReadRepresentation(Int8) : ReadRepresentation(UInt8)
    elseif dt.size == 4
        return signed ? ReadRepresentation(Int32) : ReadRepresentation(UInt32)
    elseif dt.size == 2
        return signed ? ReadRepresentation(Int16) : ReadRepresentation(UInt16)
    elseif dt.size == 16
        return signed ? ReadRepresentation(Int128) : ReadRepresentation(UInt128)
    else
        throw(UnsupportedFeatureException())
    end
end

h5fieldtype(::JLDFile, ::Type{Float16}, ::Type{Float16}, ::Initialized) =
    FloatingPointDatatype(DT_FLOATING_POINT, 0x20, 0x0f, 0x00, 2, 0, 16, 10, 5, 0, 10, 0x0000000f)
h5fieldtype(::JLDFile, ::Type{Float32}, ::Type{Float32}, ::Initialized) =
    FloatingPointDatatype(DT_FLOATING_POINT, 0x20, 0x1f, 0x00, 4, 0, 32, 23, 8, 0, 23, 0x0000007f)
h5fieldtype(::JLDFile, ::Type{Float64}, ::Type{Float64}, ::Initialized) =
    FloatingPointDatatype(DT_FLOATING_POINT, 0x20, 0x3f, 0x00, 8, 0, 64, 52, 11, 0, 52, 0x000003ff)

function jltype(f::JLDFile, dt::FloatingPointDatatype)
    if dt == h5fieldtype(f, Float64, Float64, Val{true})
        return ReadRepresentation(Float64)
    elseif dt == h5fieldtype(f, Float32, Float32, Val{true})
        return ReadRepresentation(Float32)
    elseif dt == h5fieldtype(f, Float16, Float16, Val{true})
        return ReadRepresentation(Float16)
    else
        throw(UnsupportedFeatureException())
    end
end

function h5fieldtype(f::JLDFile, writeas::PrimitiveTypeTypes, readas::Type, init::Initialized)
    @lookup_committed f readas
    commit(f, h5fieldtype(f, writeas, writeas, init), writeas, readas)
end
h5type(f::JLDFile, writeas::PrimitiveTypeTypes, x) =
    h5fieldtype(f, writeas, typeof(x), Val{true})

# Used only for custom serialization
constructrr(f::JLDFile, T::PrimitiveTypeTypes, dt::Union{FixedPointDatatype,FloatingPointDatatype},
            ::Vector{ReadAttribute}) =
    dt == h5fieldtype(f, T, T, Val{true}) ? (ReadRepresentation(T), true) :
                                            throw(UnsupportedFeatureException())

## References

h5type(::JLDFile, ::Type{RelOffset}, ::RelOffset) = ReferenceDatatype()
odr(::Type{RelOffset}) = RelOffset

function h5convert!(out::Ptr, odr::Type{RelOffset}, f::JLDFile, x::Any, wsession::JLDWriteSession)
    unsafe_store!(convert(Ptr{RelOffset}, out), write_ref(f, x, wsession))
    nothing
end
h5convert_uninitialized!(out::Ptr, odr::Type{RelOffset}) =
    (unsafe_store!(convert(Ptr{RelOffset}, out), NULL_REFERENCE); nothing)

# Reading references as references
jlconvert(::ReadRepresentation{RelOffset,RelOffset}, f::JLDFile, ptr::Ptr,
          ::RelOffset) =
    unsafe_load(convert(Ptr{RelOffset}, ptr))
jlconvert_canbeuninitialized(::ReadRepresentation{RelOffset,RelOffset}) = false

# Reading references as other types
@inline function jlconvert{T}(::ReadRepresentation{T,RelOffset}, f::JLDFile, ptr::Ptr,
                              ::RelOffset)
    x = read_dataset(f, unsafe_load(convert(Ptr{RelOffset}, ptr)))
    (isa(x, T) ? x : convert(T, x))::T
end
jlconvert_canbeuninitialized{T}(::ReadRepresentation{T,RelOffset}) = true
jlconvert_isinitialized{T}(::ReadRepresentation{T,RelOffset}, ptr::Ptr) =
    unsafe_load(convert(Ptr{RelOffset}, ptr)) != NULL_REFERENCE

## Routines for variable-length datatypes

# Write variable-length data and store the offset and length to out pointer
@inline function store_vlen!(out::Ptr, odr, f::JLDFile, x::AbstractVector, wsession::JLDWriteSession)
    unsafe_store!(convert(Ptr{UInt32}, out), length(x))
    unsafe_store!(convert(Ptr{GlobalHeapID}, out)+4, write_heap_object(f, odr, x, wsession))
    nothing
end
h5convert!{T}(out::Ptr, ::Type{Vlen{T}}, f::JLDFile, x, wsession::JLDWriteSession) =
    store_vlen!(out, T, f, x, wsession)

h5convert_uninitialized!{T<:Vlen}(out::Ptr, odr::Type{T}) =
    (unsafe_store!(convert(Ptr{Int128}, out), 0); nothing)

# Read variable-length data given offset and length in ptr
jlconvert{T,S}(::ReadRepresentation{T,Vlen{S}}, f::JLDFile, ptr::Ptr, ::RelOffset) =
    read_heap_object(f, unsafe_load(convert(Ptr{GlobalHeapID}, ptr+4)), ReadRepresentation{T, S}())
jlconvert_canbeuninitialized{T,S}(::ReadRepresentation{T,Vlen{S}}) = true
jlconvert_isinitialized{T,S}(::ReadRepresentation{T,Vlen{S}}, ptr::Ptr) =
    unsafe_load(convert(Ptr{GlobalHeapID}, ptr+4)) != GlobalHeapID(RelOffset(0), 0)

## Strings

const H5TYPE_VLEN_UTF8 = VariableLengthDatatype(DT_VARIABLE_LENGTH, 0x11, 0x01, 0x00,
                                               sizeof(Vlen{UInt8}),
                                               FixedPointDatatype(1, false))

h5fieldtype(::JLDFile, ::Type{String}, ::Type{String}, ::Initialized) =
    H5TYPE_VLEN_UTF8
function h5fieldtype(f::JLDFile, writeas::Type{String},
                     readas::Type, init::Initialized)
    @lookup_committed f readas
    commit(f, h5fieldtype(f, writeas, writeas, init), writeas, readas)
end

fieldodr(::Type{String}, ::Bool) = Vlen{String}

# Stored as variable-length strings
struct FixedLengthString{T<:AbstractString}
    length::Int
end
sizeof(x::FixedLengthString{String}) = x.length

h5type(f::JLDFile, writeas::Type{String}, x::String) =
    StringDatatype(typeof(x), sizeof(x))
h5type(f::JLDFile, writeas::Type{String}, x) =
    h5fieldtype(f, writeas, typeof(x), Val{true})
odr(::Type{String}) = fieldodr(String, true)
objodr(x::String) = FixedLengthString{String}(sizeof(x))

function jltype(f::JLDFile, dt::BasicDatatype)
    if dt.class == DT_STRING
        if (dt.bitfield1 == 0x01 || dt.bitfield1 == 0x11) && dt.bitfield2 == 0x00 && dt.bitfield3 == 0x00
            return FixedLengthString{String}(dt.size)
        else
            throw(UnsupportedFeatureException())
        end
    elseif dt.class == DT_OPAQUE
        error("attempted to read a bare (non-committed) opaque datatype")
    elseif dt.class == DT_REFERENCE
        return ReadRepresentation(Any, RelOffset)
    else
        throw(UnsupportedFeatureException())
    end
end

function jltype(f::JLDFile, dt::VariableLengthDatatype)
    if dt == H5TYPE_VLEN_UTF8
        return ReadRepresentation(String, Vlen{String})
    else
        throw(UnsupportedFeatureException())
    end
end

h5convert!(out::Ptr, ::FixedLengthString, f::JLDFile, x, ::JLDWriteSession) =
    (unsafe_copy!(convert(Ptr{UInt8}, out), pointer(x), sizeof(x)); nothing)
h5convert!(out::Ptr, ::Type{Vlen{String}}, f::JLDFile, x, wsession::JLDWriteSession) =
    store_vlen!(out, UInt8, f, Vector{UInt8}(x), wsession)

jlconvert(::ReadRepresentation{String,Vlen{String}}, f::JLDFile, ptr::Ptr, ::RelOffset) =
    String(jlconvert(ReadRepresentation(UInt8, Vlen{UInt8}), f, ptr, NULL_REFERENCE))
function jlconvert(rr::FixedLengthString{String}, ::JLDFile, ptr::Ptr, ::RelOffset)
    data = Vector{UInt8}(rr.length)
    unsafe_copy!(pointer(data), convert(Ptr{UInt8}, ptr), rr.length)
    String(data)
end

# Used only for custom serialization
constructrr(::JLDFile, ::Type{String}, dt::VariableLengthDatatype{FixedPointDatatype}, ::Vector{ReadAttribute}) =
    dt == H5TYPE_VLEN_UTF8 ?
        (ReadRepresentation(String, Vlen{String}), true) :
        throw(UnsupportedFeatureException())

## Symbols

function h5fieldtype(f::JLDFile, ::Type{Symbol}, readas::Type, ::Initialized)
    @lookup_committed f readas
    commit(f, H5TYPE_VLEN_UTF8, Symbol, readas)
end
fieldodr(::Type{Symbol}, ::Bool) = Vlen{String}

h5type(f::JLDFile, ::Type{Symbol}, x) = h5fieldtype(f, Symbol, typeof(x), Val{true})
odr(::Type{Symbol}) = Vlen{String}

h5convert!(out::Ptr, ::Type{Vlen{String}}, f::JLDFile, x::Symbol, ::JLDWriteSession) =
    store_vlen!(out, UInt8, f, Vector{UInt8}(String(x)), f.datatype_wsession)

constructrr(::JLDFile, ::Type{Symbol}, dt::VariableLengthDatatype{FixedPointDatatype}, ::Vector{ReadAttribute}) =
    dt == H5TYPE_VLEN_UTF8 ? (ReadRepresentation(Symbol, Vlen{String}), true) :
                             throw(UnsupportedFeatureException())

jlconvert(::ReadRepresentation{Symbol,Vlen{String}}, f::JLDFile, ptr::Ptr, ::RelOffset) =
    Symbol(jlconvert(ReadRepresentation(UInt8, Vlen{UInt8}), f, ptr, NULL_REFERENCE))

## BigInts and BigFloats

writeas(::Union{Type{BigInt},Type{BigFloat}}) = String
wconvert(::Type{String}, x::BigInt) = base(62, x)
wconvert(::Type{String}, x::BigFloat) = string(x)
rconvert(::Type{BigInt}, x::String) = parse(BigInt, x, 62)
rconvert(::Type{BigFloat}, x::String) = parse(BigFloat, x)

## DataTypes

const H5TYPE_DATATYPE = CompoundDatatype(
    sizeof(Vlen{String})+sizeof(Vlen{RelOffset}),
    ["name", "parameters"],
    [0, sizeof(Vlen{String})],
    [H5TYPE_VLEN_UTF8, VariableLengthDatatype(ReferenceDatatype())]
)

function h5fieldtype{T<:DataType}(f::JLDFile, ::Type{T}, readas::Type, ::Initialized)
    if !(readas <: DataType)
        @lookup_committed f readas
        return commit(f, H5TYPE_DATATYPE, DataType, readas)
    end

    @lookup_committed f DataType
    io = f.io
    offset = f.end_of_data

    seek(io, offset)
    id = length(f.datatypes)+1
    h5o = h5offset(f, offset)
    cdt = CommittedDatatype(h5o, id)
    f.datatype_locations[h5o] = cdt
    f.jlh5type[DataType] = cdt
    f.h5jltype[cdt] = ReadRepresentation(DataType, DataTypeODR())
    push!(f.datatypes, H5TYPE_DATATYPE)

    commit(f, H5TYPE_DATATYPE, (WrittenAttribute(:julia_type, WriteDataspace(f, DataType, odr(DataType)), cdt, DataType),))

    cdt
end
fieldodr{T<:DataType}(::Type{T}, ::Bool) = DataTypeODR()

h5type{T<:DataType}(f::JLDFile, ::Type{T}, x) =
    h5fieldtype(f, DataType, typeof(x), Val{true})
odr{T<:DataType}(::Type{T}) = DataTypeODR()

function typename(T::DataType)
    tn = Symbol[]
    m = T.name.module
    while m != module_parent(m)
        push!(tn, module_name(m))
        m = module_parent(m)
    end
    reverse!(tn)
    push!(tn, T.name.name)
    join(tn, ".")
end

function h5convert!(out::Ptr, ::DataTypeODR, f::JLDFile, T::DataType, wsession::JLDWriteSession)
    store_vlen!(out, UInt8, f, Vector{UInt8}(typename(T)), f.datatype_wsession)
    if !isempty(T.parameters)
        refs = RelOffset[begin
            if isa(x, DataType)
                # The heuristic here is that, if the field type is a committed
                # data type, then we commit the datatype and write it as a
                # reference to the committed datatype. Otherwise we write it
                # as a name. This ensures that type parameters that affect the
                # structure of a type are written to the file, so that we can
                # reconstruct the type when the layout depends on the
                # parameters.
                dt = h5fieldtype(f, writeas(x), x, Val{true})
                if isa(dt, CommittedDatatype)
                    (dt::CommittedDatatype).header_offset
                else
                    write_ref(f, x, wsession)
                end
            else
                write_ref(f, x, wsession)
            end
        end for x in T.parameters]
        store_vlen!(out+sizeof(Vlen{UInt8}), RelOffset, f, refs, f.datatype_wsession)
    end
    nothing
end

# Read a type. Returns an instance of UnknownType if the type or parameters
# could not be resolved.
function jlconvert{T}(::ReadRepresentation{T,DataTypeODR()}, f::JLDFile,
                      ptr::Ptr, header_offset::RelOffset)
    hasparams = unsafe_load(convert(Ptr{UInt32}, ptr+sizeof(Vlen{UInt8}))) != 0
    unknown_params = false
    if hasparams
        paramrefs = jlconvert(ReadRepresentation(RelOffset, Vlen{RelOffset}), f,
                              ptr+sizeof(Vlen{UInt8}), NULL_REFERENCE)
        params = Any[begin
            # If the reference is to a committed datatype, read the datatype
            nulldt = CommittedDatatype(UNDEFINED_ADDRESS, 0)
            cdt = get(f.datatype_locations, ref, nulldt)
            res = cdt !== nulldt ? (typeof(jltype(f, cdt)::ReadRepresentation)::DataType).parameters[1] : read_dataset(f, ref)
            unknown_params = unknown_params || isa(res, UnknownType)
            res
        end for ref in paramrefs]
    end

    path = String(jlconvert(ReadRepresentation(UInt8, Vlen{UInt8}), f, ptr, NULL_REFERENCE))
    parts = split(path, '.')
    m = Main
    for part in parts
        sym = Symbol(part)
        if !isdefined(m, sym)
            return hasparams ? UnknownType(path, params) : UnknownType(path)
        end
        m = getfield(m, sym)
    end
    if !isa(m, DataType) && !isa(m, UnionAll)
        return hasparams ? UnknownType(path, params) : UnknownType(path)
    end

    if hasparams
        unknown_params && return UnknownType(m, params)
        try
            m = m{params...}
        catch e
            return UnknownType(m, params)
        end
    elseif m === Tuple
        # Need to instantiate with no parameters, since Tuple is really
        # Tuple{Vararg{Any}}
        m = Tuple{}
    end
    track_weakref!(f, header_offset, m)
    m
end

constructrr{T<:DataType}(::JLDFile, ::Type{T}, dt::CompoundDatatype, ::Vector{ReadAttribute}) =
    dt == H5TYPE_DATATYPE ? (ReadRepresentation(DataType, DataTypeODR()), true) :
                            throw(UnsupportedFeatureException())

## Union Types

const H5TYPE_UNION = VariableLengthDatatype(H5TYPE_DATATYPE)

function h5fieldtype{T<:Union,S<:Union}(f::JLDFile, ::Type{T}, readas::Type{S}, ::Initialized)
    @lookup_committed f Union
    commit(f, H5TYPE_UNION, Union, Union)
end
function h5fieldtype{T<:Union}(f::JLDFile, ::Type{T}, readas::Type, ::Initialized)
    @lookup_committed f readas
    commit(f, H5TYPE_UNION, Union, readas)
end
fieldodr{T<:Union}(::Type{T}, ::Bool) = Vlen{DataTypeODR()}
h5fieldtype(f::JLDFile, ::Type{Union{}}, ::Initialized) = nothing
fieldodr(::Type{Union{}}, initialized::Bool) = nothing

h5type{T<:Union}(f::JLDFile, ::Type{T}, x::Any) =
    h5fieldtype(f, Union, typeof(x), Val{true})
odr(::Type{Union}) = fieldodr(Union, true)

h5convert!(out::Ptr, ::Type{Vlen{DataTypeODR()}}, f::JLDFile, x::Union, wsession::JLDWriteSession) =
    store_vlen!(out, DataTypeODR(), f, Vector{DataType}(Base.uniontypes(x)), wsession)

function jlconvert(::ReadRepresentation{Union, Vlen{DataTypeODR()}}, f::JLDFile,
                   ptr::Ptr, header_offset::RelOffset)
    v = Union{jlconvert(ReadRepresentation(DataType, Vlen{DataTypeODR()}), f, ptr, NULL_REFERENCE)...}
    track_weakref!(f, header_offset, v)
    v
end

constructrr{T<:Union}(::JLDFile, ::Type{T}, dt::VariableLengthDatatype, ::Vector{ReadAttribute}) =
    dt == H5TYPE_UNION ? (ReadRepresentation(Union, Vlen{DataTypeODR()}), true) :
                         throw(UnsupportedFeatureException())

## UnionAll

# This needs its own h5convert! method, since otherwise we will attempt to specialize the
# generic h5convert! method for the specific UnionAll type rather than for UnionAll
# more generally.
function h5convert!(out::Ptr, odr::OnDiskRepresentation{(0, 8),Tuple{TypeVar,Any},Tuple{JLD2.RelOffset,JLD2.RelOffset}},
                    f::JLDFile, x::UnionAll, wsession::JLDWriteSession)
    h5convert!(out, RelOffset, f, x.var, f.datatype_wsession)
    h5convert!(out+sizeof(RelOffset), RelOffset, f, x.body, f.datatype_wsession)
end

## Pointers

struct PointerException <: Exception; end
Base.showerror(io::IO, ::PointerException) = print(io, "cannot write a pointer to JLD file")
h5fieldtype{T<:Ptr}(::JLDFile, ::Type{T}, ::Type, ::Initialized) = throw(PointerException())
h5type{T<:Ptr}(::JLDFile, ::Type{T}, ::ANY) = throw(PointerException())

## Arrays

h5fieldtype{T<:Array}(::JLDFile, ::Type{T}, ::Type{T}, ::Initialized) =
    ReferenceDatatype()
fieldodr{T<:Array}(::Type{T}, ::Bool) = RelOffset

@inline function odr{T,N}(::Type{Array{T,N}})
    writtenas = writeas(T)
    CustomSerialization(writtenas, T, fieldodr(writtenas, false))
end

# This is all so that when you define a writeas method to write something as an
# array, it writes a reference to the actual array where the datatype is
# committed and has a written_type attribute.
function h5fieldtype{T<:Array}(f::JLDFile, ::Type{T}, readas::DataType,
                               ::Initialized)
    @lookup_committed f readas
    commit(f, ReferenceDatatype(), T, readas)
end
h5type{T<:Array}(f::JLDFile, ::Type{T}, x) =
    h5fieldtype(f, T, typeof(x), Val{true})
_odr{T<:Array}(writtenas::Type{T}, readas::Type{T}, odr) = odr
_odr{T<:Array}(writtenas::Type{T}, readas::DataType, odr) =
    CustomSerialization{writtenas,RelOffset}

function constructrr{T<:Array}(::JLDFile, ::Type{T}, dt::BasicDatatype,
                               attrs::Vector{ReadAttribute})
    dt.class == DT_REFERENCE || throw(UnsupportedFeatureException())
    (ReadRepresentation{Array, RelOffset}(), true)
end

## SimpleVectors

writeas(::Type{SimpleVector}) = Vector{Any}
wconvert(::Type{Vector{Any}}, x::SimpleVector) = Any[x for x in x]
rconvert(::Type{SimpleVector}, x::Vector{Any}) = Core.svec(x...)

## Dicts

writeas{K,V}(::Type{Dict{K,V}}) = Vector{Pair{K,V}}
writeas(::Type{ObjectIdDict}) = Vector{Pair{Any,Any}}
wconvert{K,V}(::Type{Vector{Pair{K,V}}}, x::Associative{K,V}) = collect(x)
function rconvert{T<:Associative,K,V}(::Type{T}, x::Vector{Pair{K,V}})
    d = T()
    isa(d, Dict) && sizehint!(d::Dict, length(x))
    for (k,v) in x
        d[k] = v
    end
    d
end

## Type reconstruction

module ReconstructedTypes end

function reconstruct_bitstype(name::Union{Symbol,String}, size::Integer, empty::Bool)
    sym = gensym(name)
    eval(ReconstructedTypes, empty ? :(struct $(sym) end) : :(primitive type $(sym) $(size*8) end))
    T = getfield(ReconstructedTypes, sym)
    (ReadRepresentation(T, empty ? nothing : T), false)
end

function constructrr(f::JLDFile, unk::UnknownType, dt::BasicDatatype,
                     attrs::Vector{ReadAttribute})
    warn("type ", typestring(unk), " does not exist in workspace; reconstructing")
    reconstruct_bitstype(typestring(unk), dt.size, check_empty(attrs))
end

# Convert an ordinary type or an UnknownType to a corresponding string. This is
# only used to create gensymmed names for reconstructed types.
function typestring(T::UnknownType)
    tn = IOBuffer()
    if isa(T.name, DataType)
        write(tn, typename(T.name))
    else
        print(tn, T.name)
    end
    if isdefined(T, :parameters)
        write(tn, '{')
        for i = 1:length(T.parameters)
            x = T.parameters[i]
            i != 1 && write(tn, ',')
            if isa(x, UnknownType)
                write(tn, typestring(x))
            else
                print(tn, x)
            end
        end
        write(tn, '}')
    end
    String(take!(tn))
end

# If we managed to resolve the name, but not the parameters, then try calling the
# constructor.
function constructrr(f::JLDFile, unk::UnknownType{T} where T<:Union{UnionAll,DataType},
                     dt::CompoundDatatype, attrs::Vector{ReadAttribute})
    if unk.name == Tuple
        field_datatypes = read_field_datatypes(f, attrs)
        # For a tuple with unknown fields, we should reconstruct the fields
        rodr = reconstruct_odr(f, dt, field_datatypes)
        # This is a "pseudo-RR" since the tuple is not fully parametrized, but
        # the parameters must depend on the types actually encoded in the file
        (ReadRepresentation(Tuple, rodr), false)
    else
    end
end

"""
    behead(T)

Given a UnionAll type, recursively eliminates the `where` clauses
"""
behead(T::UnionAll) = behead(T.body)
behead(T::ANY) = T

function getfieldtypes(f::JLDFile, dt::CompoundDatatype, attrs::Vector{ReadAttribute})
    datatypes = read_field_datatypes(f, attrs)
    jltypes = Vector{Any}(length(datatypes))
    for i = 1:length(datatypes)
        ref = datatypes[i]
        jltypes[i] = jltype(f, ref == NULL_REFERENCE ? dt.members[i] : f.datatype_locations[ref])
    end
    jltypes
end

function constructrr(f::JLDFile, unk::UnknownType{UnionAll}, dt::CompoundDatatype,
                     attrs::Vector{ReadAttribute})
    field_datatypes = read_field_datatypes(f, attrs)
    body = behead(unk.name)
    if length(body.parameters) != length(unk.parameters)
        warn("read type ", typestring(unk), " has a different number of parameters from type ",
             unk.name, " in workspace; reconstructing")
        reconstruct_compound(f, typestring(unk), dt, field_datatypes)
    else
        params = copy(unk.parameters)
        for i = 1:length(params)
            if isa(params[i], UnknownType)
                param = body.parameters[i]::TypeVar
                params[i] = param.ub
            end
        end

        # Try to construct the rr for the relaxed type. On failure, fall back to
        # reconstruct_compound.
        local T
        try
            T = unk.name{params...}
        catch err
            warn("type parameters for ", typestring(unk), " do not match type ", unk.name, " in workspace; reconstructing")
            return reconstruct_compound(f, typestring(unk), dt, field_datatypes)
        end

        try
            (rr,) = constructrr(f, T, dt, attrs, true)
            warn("some parameters could not be resolved for type ", typestring(unk), "; reading as ", T)
            return (rr, false)
        catch err
            !isa(err, TypeMappingException) && rethrow(err)
            warn("some parameters could not be resolved for type ", typestring(unk), "; reconstructing")
            return reconstruct_compound(f, typestring(unk), dt, field_datatypes)
        end
    end
end

# An UnknownType for which we were not able to resolve the name.
function constructrr(f::JLDFile, unk::UnknownType{String}, dt::CompoundDatatype,
                     attrs::Vector{ReadAttribute})
    ts = typestring(unk)
    warn("type ", ts, " does not exist in workspace; reconstructing")
    reconstruct_compound(f, ts, dt, read_field_datatypes(f, attrs))
end

# Reconstruct the ODR of a type from the CompoundDatatype and field_datatypes
# attribute
function reconstruct_odr(f::JLDFile, dt::CompoundDatatype,
                         field_datatypes::Vector{RelOffset})
    # Get the type and ODR information for each field
    types = Vector{Any}(length(dt.names))
    h5types = Vector{Any}(length(dt.names))
    for i = 1:length(dt.names)
        if !isempty(field_datatypes) && (ref = field_datatypes[i]) != NULL_REFERENCE
            dtrr = jltype(f, f.datatype_locations[ref])
        else
            dtrr = jltype(f, dt.members[i])
        end
        types[i], h5types[i] = typeof(dtrr).parameters
    end
    return OnDiskRepresentation{(dt.offsets...), Tuple{types...}, Tuple{h5types...}}()
end

# Reconstruct type that is a "lost cause": either we were not able to resolve
# the name, or the workspace type has additional fields, or cannot convert
# fields to workspace types
function reconstruct_compound(f::JLDFile, T::String, dt::H5Datatype,
                              field_datatypes::Union{Vector{RelOffset},Void})
    rodr = reconstruct_odr(f, dt, field_datatypes)
    types = typeof(rodr).parameters[2].parameters

    # Now reconstruct the type
    reconname = gensym(T)
    eval(ReconstructedTypes, Expr(:composite_type, reconname, Core.svec(),
                                  Core.svec(dt.names...), Any,
                                  types, false, 0))
    T = getfield(ReconstructedTypes, reconname)

    (ReadRepresentation(T, rodr), false)
end

# These need to go at the bottom. Also, JLD2 doesn't support custom serialization because
# these methods are not guaranteed to work if you add methods to `writeas`.

@generated function h5type{T<:Array}(f::JLDFile, ::Type{T}, ::T)
    if T <: Array{Union{}}
        return :(ReferenceDatatype())
    end
    ty = T.parameters[1]
    writtenas = writeas(ty)
    !hasfielddata(writtenas) ? :(h5type(f, $writtenas, $(Expr(:new, ty)))) : :(h5fieldtype(f, $writtenas, $ty, Val{false}))
end

# jlconvert for empty objects
@generated function jlconvert{T}(::ReadRepresentation{T,nothing}, f::JLDFile, ptr::Ptr,
                                 header_offset::RelOffset)
   T.size == 0 && return Expr(:new, T)

   # In this case, T is a non-empty object, but the written data was empty
   # because the custom serializers for the fields all resulted in empty
   # objects
   return Expr(:new, T, [begin
       writtenas = writeas(ty)
       @assert writtenas.size == 0
       if writtenas === ty
           Expr(:new, ty)
       else
           :(rconvert($ty, $(Expr(:new, writtenas))))
       end
   end for ty in T.types]...)
end

# This jlconvert method handles compound types with padding or references
@generated function jlconvert{T,S}(::ReadRepresentation{T,S}, f::JLDFile, ptr::Ptr,
                                   header_offset::RelOffset)
    isa(S, DataType) && return :(convert(T, unsafe_load(convert(Ptr{S}, ptr))))
    @assert isa(S, OnDiskRepresentation)

    offsets = typeof(S).parameters[1]
    types = typeof(S).parameters[2].parameters
    odrs = typeof(S).parameters[3].parameters

    blk = Expr(:block)
    args = blk.args
    if isbits(T)
        # For bits types, we should always inline, because otherwise we'll just
        # pass a lot of crap around in registers
        push!(args, Expr(:meta, :inline))
    elseif T.mutable
        push!(args, quote
            obj = $(Expr(:new, T))
            track_weakref!(f, header_offset, obj)
        end)
    end
    fsyms = []
    fn = T === Tuple ? [Symbol(i) for i = 1:length(types)] : fieldnames(T)
    for i = 1:length(types)
        offset = offsets[i]
        rtype = types[i]
        odr = odrs[i]

        fsym = T.mutable ? Expr(:., :obj, QuoteNode(fn[i])) : Symbol("field_", fn[i])
        push!(fsyms, fsym)

        rr = ReadRepresentation{rtype,odr}()

        if odr === nothing
            # Type is not stored or single instance
            if T.types[i] == Union{}
                # This cannot be defined
                push!(args, Expr(:return, Expr(:new, T, fsyms[1:i-1]...)))
                return blk
            else
                push!(args, :($fsym = $(Expr(:new, T.types[i]))))
            end
        else
            if jlconvert_canbeuninitialized(rr)
                push!(args, quote
                    if !jlconvert_isinitialized($rr, ptr+$offset)
                        $(if T <: Tuple || i <= T.ninitialized
                            # Reference must always be initialized
                            :(throw(UndefinedFieldException(T,$(QuoteNode(fn[i])))))
                        elseif T.mutable
                            # Reference could be uninitialized
                            :(return obj)
                        else
                            Expr(:return, Expr(:new, T, fsyms[1:i-1]...))
                        end)
                    end
                end)
            end

            if T === Tuple
                # Special case for reconstructed tuples, where we don't know the
                # field types in advance
                push!(args, :($fsym = jlconvert($rr, f, ptr+$offset, NULL_REFERENCE)::$rtype))
            else
                ttype = T.types[i]
                push!(args, :($fsym = convert($ttype, jlconvert($rr, f, ptr+$offset, NULL_REFERENCE)::$rtype)::$ttype))
            end
        end
    end

    push!(args, T.mutable ? (:obj) : T <: Tuple ? Expr(:tuple, fsyms...) : Expr(:new, T, fsyms...))

    blk
end

# odr gives the on-disk representation of a given type, similar to
# fieldodr, but actually encoding the data for things that odr stores
# as references
@generated function odr{T}(::Type{T})
    if !hasdata(T)
        # A pointer singleton or ghost. We need to write something, but we'll
        # just write a single byte.
        return nothing
    elseif isbits(T) && samelayout(T)
        # Has a specialized convert method or is an unpadded type
        return T
    end

    offsets = zeros(Int, length(T.types))
    odrs = Vector{Any}(length(T.types))
    offset = 0
    for i = 1:length(T.types)
        ty = T.types[i]
        writtenas = writeas(ty)
        fodr = fieldodr(writtenas, i <= T.ninitialized)
        if writtenas !== ty && fodr !== nothing
            odrs[i] = CustomSerialization{writtenas,fodr}
        else
            odrs[i] = fodr
        end
        offsets[i] = offset
        offset += sizeof(fodr)
    end

    OnDiskRepresentation{(offsets...), Tuple{T.types...}, Tuple{odrs...}}()
end


text:
module GuillotineIfAnd

include("./../../constraints/geometric/constants.jl")
using JuMP
using .Constants

export all_dimensions_and

function all_dimensions_and(m)

    #The "left" side of a box
    @variable(m, guillotine_flbx[i=1:length(boxes)-1, j=1:length(boxes)-1; i!=j], Bin)
    @constraint(m, guillotine_x_cons_flbx[i=1:length(boxes)-1, j=1:length(boxes)-1; i!=j],
        m[:if1_eq_flbx][i, j] + m[:y_for_x_guillo][i, j] >=
        2 * m[:guillotine_flbx][i, j]
    )

    @variable(m, guillotine_flbx_j_plus[i=1:length(boxes)-1, j=1:length(boxes)-1; i!=j], Bin)
    @constraint(m, guillotine_x_cons_j_plus[i=1:length(boxes)-1, j=1:length(boxes)-1; i!=j],
        m[:if1_eq_flbx_j_plus][i, j] + m[:y_for_x_guillo][i, j] >=
        2 * m[:guillotine_flbx_j_plus][i, j]
    )

    #The "right" side of a box
    @variable(m, guillotine_flbx_ij_plus[i=1:length(boxes)-1, j=1:length(boxes)-1; i!=j], Bin)
    @constraint(m, guillotine_x_cons_ij_plus[i=1:length(boxes)-1, j=1:length(boxes)-1; i!=j],
        m[:if1_eq_flbx_ij_plus][i, j] + m[:y_for_x_guillo][i, j] >=
        2 * m[:guillotine_flbx_ij_plus][i, j]
    )

    @variable(m, guillotine_flbx_i_plus[i=1:length(boxes)-1, j=1:length(boxes)-1; i!=j], Bin)
    @constraint(m, guillotine_x_cons_i_plus[i=1:length(boxes)-1, j=1:length(boxes)-1; i!=j],
        m[:if1_eq_flbx_i_plus][i, j] + m[:y_for_x_guillo][i, j] >=
        2 * m[:guillotine_flbx_i_plus][i, j]
    )

    # ===== Y AXIS =====

    #The "left" side of a box
    @variable(m, guillotine_flby[i=1:length(boxes)-1, j=1:length(boxes)-1; i!=j], Bin)
    @constraint(m, guillotine_y_cons_flby[i=1:length(boxes)-1, j=1:length(boxes)-1; i!=j],
        m[:if1_eq_flby][i, j] + m[:x_for_y_guillo][i, j] >=
        2 * m[:guillotine_flby][i, j]
    )

    @variable(m, guillotine_flby_j_plus[i=1:length(boxes)-1, j=1:length(boxes)-1; i!=j], Bin)
    @constraint(m, guillotine_y_cons_j_plus[i=1:length(boxes)-1, j=1:length(boxes)-1; i!=j],
        m[:if1_eq_flby_j_plus][i, j] + m[:x_for_y_guillo][i, j] >=
        2 * m[:guillotine_flby_j_plus][i, j]
    )

    #The "right" side of a box
    @variable(m, guillotine_flby_ij_plus[i=1:length(boxes)-1, j=1:length(boxes)-1; i!=j], Bin)
    @constraint(m, guillotine_y_cons_ij_plus[i=1:length(boxes)-1, j=1:length(boxes)-1; i!=j],
        m[:if1_eq_flby_ij_plus][i, j] + m[:x_for_y_guillo][i, j] >=
        2 * m[:guillotine_flby_ij_plus][i, j]
    )

    @variable(m, guillotine_flby_i_plus[i=1:length(boxes)-1, j=1:length(boxes)-1; i!=j], Bin)
    @constraint(m, guillotine_y_cons_i_plus[i=1:length(boxes)-1, j=1:length(boxes)-1; i!=j],
        m[:if1_eq_flby_i_plus][i, j] + m[:x_for_y_guillo][i, j] >=
        2 * m[:guillotine_flby_i_plus][i, j]
    )

    @expression(m, guillotine_and_objective[i=1:length(boxes)-1, j=1:length(boxes)-1; i!=j],
        #x
        (1 - m[:guillotine_flbx][i, j]) +
        (1 - m[:guillotine_flbx_i_plus][i, j]) +
        (1 - m[:guillotine_flbx_j_plus][i, j]) +
        (1 - m[:guillotine_flbx_ij_plus][i, j]) +

        #y
        (1 - m[:guillotine_flby][i, j]) +
        (1 - m[:guillotine_flby_i_plus][i, j]) +
        (1 - m[:guillotine_flby_j_plus][i, j]) +
        (1 - m[:guillotine_flby_ij_plus][i, j])
    )

    return m
end

end

ext: jl


text:
using CSV, JLD2, Distributions, Dates
using InteractionNetworkModels

parse(Int, "10")
CSV.File()
Distributions.Uniform(0.0,10.0)


text:
type HistPlot <: Graph
	dataGroups::Vector{DataGroup}
end

HistPlot(dg::DataGroup) = HistPlot([dg])

function plotData(g::HistPlot, showLegend::Bool, fid)
	global dgCount

	for x in 1:length(g.dataGroups)
		if x == 1
			println(fid, "histogram")
		end

		dg = g.dataGroups[x]
		if dg.lineColor == "blank"
			dg.lineColor = goodColors[mod(x, length(goodColors))]
		end

		if dg.markerColor == "blank"
			dg.markerColor = goodColors[mod(x, length(goodColors))]
		end

		if dg.markerType == "blank"
			dg.markerType = goodSymbols[mod(x, length(goodSymbols))]
		end

		println(fid, "color "*dg.lineColor)
		println(fid, "line-style "*dg.lineStyle)
		println(fid, "marker-color "*dg.markerColor)
		println(fid, "marker "*dg.markerType)
		println(fid, "marker-line-width $(dg.markerLineWidth)")
		println(fid, "marker-min-scale $(dg.markerScale[1])")
		println(fid, "marker-scale $(dg.markerScale[2])")
		println(fid, "error-bar-line-width $(dg.errorLineWidth)")
		if showLegend
			dgCount += 1
			if dg.legend.label == "\'Data Group \'"
				println(fid, "legend "*dg.legend.label*string(dgCount))
			else
				println(fid, "legend "*dg.legend.label)
			end
		end

		push!(dataFs, sparrowD*"data.sparrow"*randstring(5))
		if isempty(find(DataFrames.names(dg) .== :yerr))
			writedlm(dataFs[end], [dg.data[:x] dg.data[:y]])
			println(fid, "plot $(dataFs[end])@\$1:\$2")
		else
			writedlm(dataFs[end], [dg.data[:x] dg.data[:y] dg.data[:yerr]])
			println(fid, "plot $(dataFs[end])@\$1:\$2:yerr=\$3")
		end
	end
end

function histplot(dgs::Vector{DataGroup})
	plot(PlotFrame(HistPlot(dgs)))
end

function histplot(dg::DataGroup)
	plot(PlotFrame(HistPlot(dg)))
end

function histplot(dfs::Vector{DataFrame})
	dgs = DataGroup[]
	for df in dfs
		push!(dgs, DataGroup(df))
	end
	histplot(dgs)
end

function histplot(df::DataFrame)
	histplot(DataGroup(df))
end

function histplot(ys::Array{Float64, 2})
	dfs = DataFrame[]
	for r = 1:size(y, 1)
		push!(dfs, DataFrame(x = collect(1:length(y)), y = y[r,:]))
	end
	histplot(dfs)
end

function histplot(xs::Array{Float64, 2}, ys::Array{Float64, 2})
	dfs = DataFrame[]
	for r = 1:size(y, 1)
		push!(dfs, DataFrame(x = x[r,:], y = y[r,:]))
	end
	histplot(dfs)
end

function histplot(xs::Array{Float64, 2}, ys::Array{Float64, 2}, yerrs::Array{Float64, 2})
	dfs = DataFrame[]
	for r = 1:size(ys, 1)
		push!(dfs, DataFrame(x = vec(xs[r,:]), y = vec(ys[r,:]), yerr = vec(yerrs[r,:])))
	end
	histplot(dfs)
end


text:
## Figure 3 - classifier

using JuMP, Ipopt, Distributions, LaTeXStrings, Plots
pyplot()

# Generating random multivariate data.
data11 = rand(Normal(10, 1.5), 100)
data12 = rand(Normal(8, 1.5), 100)
data1 = [data11 data12]

data01 = rand(Normal(2, 1.5), 100)
data02 = rand(Normal(2, 1.5), 100)
data0 = [data01 data02]

scatter(
    data1[:, 1],
    data1[:, 2],
    label = "Pos. obs.",
    xaxis = (L"x_1", (-5,15)),
    yaxis = (L"x_2", (-5,15))
)

scatter!(data0[:, 1], data0[:, 2], label = "Neg. obs.")

savefig("classes_no_classifier.pdf")

function classify(data1, data0, δ)
    m = Model(with_optimizer(Ipopt.Optimizer))

    @variable(m, u[1:100] >= 0)
    @variable(m, v[1:100] >= 0)
    @variable(m, a[1:2])
    @variable(m, b)

    @constraint(m, error1[i = 1:100],
        sum(a[j] * data1[i, j] for j = 1:2) + b - u[i] <= -1
    )
    @constraint(m, error0[i = 1:100],
        sum(a[j] * data0[i, j] for j = 1:2) + b + v[i] >= 1
    )

    # This parameter trades off accuracy and robustness for the classifier
    @objective(m, Min,
        sum(u[i] for i = 1:100) + sum(v[i] for i = 1:100) +
        δ * (sum(a[i]^2 for i = 1:2))
    )

    optimize!(m)

    return value.(a), value(b)
end

# Classify without robustness
a, b = classify(data1, data0, 0.0)

# Calculating the points on plane defined by a and b
x1 = range(-10, stop = 20, length = 100)
x2 = -(b .+ a[1] .* x1) / a[2]

plot!(
    x1,
    x2,
    xaxis = (L"x_1", (-5,15)),
    yaxis = (L"x_2", (-5,15)),
    label = "Classifier",
    color = 3,
)

savefig("classes_with_classifier.pdf")

# Classify with robustness
δ = 0.1
a, b = classify(data1, data0, δ)

# Boundaries of the slabs
x2u = -(b .+ 1 .+ a[1] .* x1) / a[2]
x2d = -(b .- 1 .+ a[1] .* x1) / a[2]
x2 = -(b .+ a[1] .* x1) / a[2]

plot!(
    x1,
    [x2u x2 x2d],
    xaxis = (L"x_1", (-5,15)),
    yaxis = (L"x_2"),
    label = [" " L"Classifier $\gamma = 0.1$" " "],
    color = [4 4 4],
    linestyle = [:dot :solid :dot]
)

savefig("classes_with_robust_classifier.pdf")


text:
#############################################################################
# diagm.jl
# Converts a vector of size n into an n x n diagonal
# All expressions and atoms are subtpyes of AbstractExpr.
# Please read expressions.jl first.
#############################################################################

import LinearAlgebra.diagm, LinearAlgebra.Diagonal
export diagm, Diagonal

struct DiagMatrixAtom <: AbstractExpr
    head::Symbol
    id_hash::UInt64
    children::Tuple{AbstractExpr}
    size::Tuple{Int, Int}

    function DiagMatrixAtom(x::AbstractExpr)
        (num_rows, num_cols) = x.size

        if num_rows == 1
            sz = num_cols
        elseif num_cols == 1
            sz = num_rows
        else
            throw(ArgumentError("Only vectors are allowed for diagm/Diagonal. Did you mean to use diag?"))
        end

        children = (x, )
        return new(:diagm, hash(children), children, (sz, sz))
    end
end

function sign(x::DiagMatrixAtom)
    return sign(x.children[1])
end

# The monotonicity
function monotonicity(x::DiagMatrixAtom)
    return (Nondecreasing(),)
end

# If we have h(x) = f o g(x), the chain rule says h''(x) = g'(x)^T f''(g(x))g'(x) + f'(g(x))g''(x);
# this represents the first term
function curvature(x::DiagMatrixAtom)
    return ConstVexity()
end

function evaluate(x::DiagMatrixAtom)
    return Diagonal(vec(evaluate(x.children[1])))
end

function diagm((d, x)::Pair{<:Integer, <:AbstractExpr})
    d == 0 || throw(ArgumentError("only the main diagonal is supported"))
    return DiagMatrixAtom(x)
end
Diagonal(x::AbstractExpr) = DiagMatrixAtom(x)

function conic_form!(x::DiagMatrixAtom, unique_conic_forms::UniqueConicForms=UniqueConicForms())
    if !has_conic_form(unique_conic_forms, x)
        sz = x.size[1]

        I = 1:sz+1:sz*sz
        J = 1:sz
        coeff = sparse(I, J, 1.0, sz * sz, sz)
        objective = conic_form!(x.children[1], unique_conic_forms)
        new_obj = coeff * objective
        cache_conic_form!(unique_conic_forms, x, new_obj)
    end
    return get_conic_form(unique_conic_forms, x)
end


text:
mutable struct AccumulatorObsGrad{A,B,C,D,E,F} <: AbstractObservableAccumulator
    res::A
    ∇res::B

    ψ_σ::C       # Buffers alls the <σ|ψ> of the denominator
    ∇ψ_σ::D
    mel_buf::E      # ⟨σ|Ô|σ'⟩ in the buffer

    ∇logψ_acc::F
end

function AccumulatorObsGrad(net::NeuralNetwork, hilb, v, batch_sz)
    bnet      = cached(net, batch_sz)
    CT        = Complex{real(out_type(net))}
    v         = state(hilb, bnet)

    acc       = AccumulatorLogGradPsi(bnet, v, batch_sz)

    w         = trainable_first(net)
    RT        = real(eltype(w))

    res       = zero(eltype(out_similar(bnet)))
    ∇res      = grad_cache(CT, net)

    ∇0_buf    = vec_data(grad_cache(net))[1]
    mel_buf   = similar(w, CT, 1, batch_sz)

    return AccumulatorObsGrad(
        res, ∇res,
        res, ∇0_buf,
        mel_buf, acc)
end

@forward AccumulatorObsGrad.∇logψ_acc Base.length, Base.count,
    isfull

accum(a::AccumulatorObsGrad) = a.∇logψ_acc

function init!(c::AccumulatorObsGrad, σ, ψ_σ, ∇ψ_σ)
    # Store temporaries
    c.ψ_σ      = ψ_σ
    c.∇ψ_σ     = ∇ψ_σ

    # Set accumulators to 0
    c.res      = 0
    for ∇res=vec_data(c.∇res)
        ∇res   .= 0
    end

    # initialize the lower accumulator.
    reset!(c, σ)
    return nothing
end

# This resets the value at the end of a process_accumulatr
function reset!(c::AccumulatorObsGrad, σ)
    c.mel_buf .= 0

    init!(accum(c), σ)
    return nothing
end

function (c::AccumulatorObsGrad)(mel::Number, cngs_l, cngs_r, v)
    isfull(c) && (process_accumulator!(c); reset!(c, v))

    # If the matrix element is zero, don't do anything
    mel == 0.0 && return c
    n_cngs_l = isnothing(cngs_l) ? 0 : length(cngs_l)
    n_cngs_r = isnothing(cngs_r) ? 0 : length(cngs_r)

    # If there are no changes, just sum it
    if n_cngs_l == 0 && n_cngs_r == 0 && false
        c.res += mel
    else
        accum(c)(cngs_l, cngs_r)
        @inbounds c.mel_buf[count(c)] = mel
    end

    return c
end

function (c::AccumulatorObsGrad)(mel::Number, cngs, v)
    isfull(c) && (process_accumulator!(c); reset!(c, v))

    # If the matrix element is zero, don't do anything
    mel == 0.0 && return acc

    n_cngs = isnothing(cngs) ? 0 : length(cngs)

    # If there are no changes, just sum it
    if n_cngs == 0 && false
        c.res += mel
    else
        accum(c)(cngs)
        @inbounds c.mel_buf[count(c)] = mel
    end

    return c
end

function process_accumulator!(c::AccumulatorObsGrad)
    count(c) == 0 && return nothing

    count_tmp = count(c)
    ∇ψ_σp, ψ_σp = process_accumulator!(accum(c))

    ∇res  = vec_data(c.∇res)[1]
    ∇ψ_σp = vec_data(∇ψ_σp)[1]
    ∇ψ_σ  = c.∇ψ_σ

    c.mel_buf .*= exp.(ψ_σp .- c.ψ_σ)
    ∇ψ_σp     .= c.mel_buf .* (∇ψ_σp)# .- ∇ψ_σ)
    if count_tmp != length(c)
        uview(∇ψ_σp, :, (count_tmp+1):length(c)) .= 0.0
    end

    sum!(∇res, ∇ψ_σp, init=false)
    c.res += sum(c.mel_buf)

    return nothing
end

function finalize!(c::AccumulatorObsGrad)
    process_accumulator!(c)
    return c.res, c.∇res
end


text:
using CairoMakie
CairoMakie.activate!()

# Um gráfico simples de linha pode ser alcançado com uma síntaxe igualmente simples:
fig = scatterlines(1:10, 1:10)

# É importante salientar que todas as funções como scatterlines criam e retornam
# novos objetos Figure, Axis e plot em uma coleção chamada FigureAxisPlot. Estes
# são os métodos não-mutantes. Em oposição, temos os métodos mutantes (e.g.
# scatterlines!, note o !) que retornam um objeto plot que pode ser adicionado
# em qualquer Axis ou aplicado com current_figure().

# Dessa forma, é possível desempacotar o retorno em 3 objetos e, por exemplo,
# adicionar mais gráficosà figura e eixo retornados pelo método:

fig, ax, pltobj = scatterlines(1:10)
lines!(ax, x, (1:10).-1; label="line")
fig

# Ainda usando as funções não mutantes, é possível costumizar os gráficos através
# de atributos. A lista de atributos disponíveis pode ser obtida com:
pltobj.attributes
# onde pltobj é um objeto plot retornado pelos métodos mutantes ou não-mutantes
# Não apenas os objetos de plot possuem atributos, claro! Axis e Figuras também
# podem ser costumizados.
# Use, por exemplo, ?Axis no terminal para obter a lista de atributos do Axis.

# Usando as funções não-mutantes, é possível definir as configurações do Axis, bem
# como da Figure.

lines(
    1:10,
    (1:10).^2;
    label="x²",
    linewidth=2,
    linestyle=nothing,
    figure=(;
        figure_padding=5,
        resolution=(600, 400),
        ont="sans",
        backgroundcolor=:grey90,
        fontsize=16
    ),
    axis=(;
        xlabel="x",
        title="title",
        xgridstyle=:dash,
        ygridstyle=:dash
    )
)
scatterlines!(1:10, (10:-1:1).^2; label="Reverse(x)²")
axislegend("legend"; position=:ct)
current_figure()

# Contudo, minha forma predileta, até agora, é usar os métodos mutantes e controlar
# separadamente o Axis e a Figure para ter mais controle e dar mais clareza ao código.

x = 1:10
y = (10:-1:1).^2
z = rand(10)

fig = Figure(resolution=(600, 400), backgroundcolor=:grey90)
ax = Axis(
    fig[1, 1],
    backgroundcolor=:white,
    xgridstyle=:dash,
    ygridstyle=:dash
)
pltobj = scatter!(ax, x, y; color=z, label="scatters")
lines!(ax, x, 1.1y; label="line")
Legend(fig[2, 1:2], ax, "labels", orientation=:horizontal)
Colorbar(fig[1, 2], pltobj, label="colorbar")
fig

max_stars_repo_name: ericmiguel/julia-notes

max_issues_repo_name: ericmiguel/julia-notes

max_forks_repo_name: ericmiguel/julia-notes


text:
# Starting with 1 and spiralling anticlockwise in the following way, a square
# spiral with side length 7 is formed.
#
# 37 36 35 34 33 32 31
# 38 17 16 15 14 13 30
# 39 18  5  4  3 12 29
# 40 19  6  1  2 11 28
# 41 20  7  8  9 10 27
# 42 21 22 23 24 25 26
# 43 44 45 46 47 48 49
#
# It is interesting to note that the odd squares lie along the bottom right
# diagonal, but what is more interesting is that 8 out of the 13 numbers lying
# along both diagonals are prime; that is, a ratio of 8/13 ≈ 62%.
#
# If one complete new layer is wrapped around the spiral above, a square spiral
# with side length 9 will be formed. If this process is continued, what is the
# side length of the square spiral for which the ratio of primes along both
# diagonals first falls below 10%?

using ProjectEulerSolutions

# Create the sequence of corner values and check each for primeness.  Keep a
# running total to verify that the desired threshold has been met.
function p058solution(thresh::Float64=0.44)::Integer
    nprimes = 0
    total_tested = 1
    startindex = 3
    i = 2
    while i == 2 || nprimes / total_tested > thresh
        nprimes += sum(isprime.(collect(startindex .+ (0:i:i*3))))
        startindex = startindex + 4*i + 2
        total_tested += 4
        i += 2
    end
    return i - 1
end

p058 = Problems.Problem(p058solution)

Problems.benchmark(p058, 0.1)

max_stars_repo_name: gnujosh/julia-euler-project

max_issues_repo_name: gnujosh/julia-euler-project

max_forks_repo_name: gnujosh/julia-euler-project


text:
datadim(X::AbstractMatrix) = size(X, 1)
datadim(X::AbstractArray) = Base.front(size(X))

_getproperty(d::AbstractDataset, ::Val{T}) where {T} = getfield(d, T)
Base.getproperty(d::AbstractDataset, k::Symbol) = _getproperty(d, Val(k))

_getproperty(d::AbstractDataset, ::Val{:dim})       = datadim(d.X)
_getproperty(d::AbstractDataset, ::Val{:n_data})    = last(size(d.X))
_getproperty(d::AbstractDataset, ::Val{:n_test})    = last(size(d.Xt))
_getproperty(d::AbstractDataset, ::Val{:n_display}) = n_display(d)
_getproperty(d::AbstractDataset, ::Val{:x_display}) = d.X[:,1:n_display(d)]
_getproperty(d::AbstractDataset, ::Val{:vis})       = (args...) -> vis(d, args...)

const DATASET_NAMES = (
    "gaussian", "2dring", "3dring", # simple
    "mnist", "cifar10",             # image
    "griffiths2011", "xu2019",      # feature
)

function Dataset(name::String, n_data::Int, args...; is_preview=true, kwargs...)
    if name == "gaussian"
        dataset = GaussianDataset(n_data, args...; kwargs...)
    end
    if name == "2dring"
        dataset = TwoDimRingDataset(n_data, args...; kwargs...)
    end
    if name == "3dring"
        dataset = ThreeDimRingDataset(n_data, args...; kwargs...)
    end
    if name == "mnist"
        dataset = MNISTDataset(n_data, args...; kwargs...)
    end
    if name == "cifar10"
        dataset = CIFAR10Dataset(n_data, args...; kwargs...)
    end
    if name == "griffiths2011"
        dataset = FeatureDataset(n_data, get_features_griffiths2011(), args...; kwargs...)
    end
    if name == "xu2019"
        dataset = FeatureDataset(n_data, get_features_xu2019(), args...; kwargs...)
    end

    if is_preview
        nd_half = div(n_display(dataset), 2)
        fig = vis(
            dataset,
            (train=dataset.X[:,1:nd_half], test=dataset.Xt[:,1:nd_half])
        )
        display(fig)
        println()
    end

    return dataset
end


text:
using ApproxFun
using EffectiveWaves

T = Float64
θin = 0.13
k=1.; ho = 2
medium = Medium(1.0,1.0+0.0im)
ω = real(k*medium.c)
specie = Specie(ρ=0.1,r=0.1, c=0.5, volfrac=0.1)

# Guess a course mesh
k_eff0 = wavenumber_low_volfrac(ω, medium, [specie]; tol = 1e-12)
max_x = 10 .* k/imag(k_eff0)
x = 0.0:0.1:max_x

wave0 = EffectiveWave(ω, k_eff0, medium, [specie]; θin = θin, hankel_order = ho, tol=1e-8)
wave_avg0 = AverageWave(x, wave0)

R_eff = reflection_coefficient(ω, wave0, medium, [specie]; θin = θin, hankel_order = ho)
R = reflection_coefficient(ω, wave_avg0, medium, specie; θin = θin)

amps = wave_avg0;
# k = ω/medium.c
M = amps.hankel_order
σ = trap_scheme(amps.x);

R_ms = T(2)*specie.num_density/(cos(θin)*k^2)*
[
    sum(
        im^T(m)*exp(-im*θin*m)*amps.amplitudes[j,m+M+1,1]*exp(im*amps.x[j]*cos(θin))*σ[j]
    for j in eachindex(amps.x))
for m = -M:M];
sum(R_ms) - R

x = 0.:0.02:10

# function space
S = Chebyshev(minimum(x)..maximum(x));
n1 = length(x); m1 = Int(round(n*0.4));
V = Array{Float64}(n1,m1); # Vandermonde matrix
for k = 1:m1
    V[:,k] = Fun(S,[zeros(k-1);1]).(x)
end

v = exp(x.*(-1.0 + 1.0im))
σ = trap_scheme(x)
sum(σ[i]*v[i] for i in eachindex(x))

f = Fun(S,V\v)
sum(Fun(S,V\v))
1/(1.0+1.0im)

R_funs = map(-M:M) do m
    v = [
        im^T(m)*exp(-im*θin*m)*amps.amplitudes[j,m+M+1,1]*exp(im*amps.x[j]*cos(θin))
    for j in eachindex(amps.x)]
    sum(Fun(S,V\v))
end

R_fun = sum(R_funs)*T(2)*specie.num_density/(cos(θin)*k^2)

    R_m = sum(Fun(S,V\v));
sum(R_m)


text:
using CoverRanges
using CoverRanges: LeftToCover, PartialCover
using Test

@testset "Impossible task" begin
    r = 1:5
    K = 5
    C = [1:4]
    @test solve_min_covering(r, K, C) === nothing
end

@testset "Simple example with type $T" for T in (Int64, Int32)
    r = T(1):T(6)
    K = T(6)
    C = [T(1):T(1), T(6):T(6), T(2):T(5), T(1):T(3), T(4):T(6)]
    @test Set(C[solve_min_covering(r, K, C)]) == Set([T(1):T(3), T(4):T(6)])
end

@testset "Simple example with negative numbers" begin
    r = -6:-1
    K = 6
    C = [-1:-1, -6:-6, -5:-2, -3:-1, -6:-4]
    @test Set(C[solve_min_covering(r, K, C)]) == Set([-3:-1, -6:-4])
end

function check_coverage(range_to_cover, ranges)
    num_covered = 0
    for node in range_to_cover
        num_covered += any( r -> node in r, ranges)
    end
    return num_covered
end

function make_random_subrange(r::UnitRange{T}) where {T}
    a, b = first(r), last(r)
    k = rand(a:b)
    l = k + abs(round(T, (length(r)/100)*randn()))
    k:l
end

@testset "Random examples of size $n" for n ∈ (50, 200, 500)
    r = -n:n
    C = [ make_random_subrange(r) for _ = 1:3n ]
    for K in round.(Int, [.2, .4, .8, 1] .* (2n))
        possible = check_coverage(r, C) >= K
        result = solve_min_covering(r, K, C)
        if !possible
            @test result === nothing
        else
            @test check_coverage(r, C[result]) >= K
        end
    end

end

@testset "AccessDict{$T}" for T in (Int16, Int32, Int64)
    n = 1000
    r = T(1):T(n)
    C = [ make_random_subrange(r) for _ = 1:2n ]
    K = round(T, .9*n)
    cache = AccessDict{LeftToCover{T}, PartialCover{T}}()
    solve_min_covering(r, K, C; cache=cache)
    @test total_sets(cache) == length(cache)
    @test total_gets(cache) <= length(cache)
end


text:
# one layer of a multilayer perceptron (neural net)
# ninᵢ := noutᵢ₋₁ + 1  (we add the bias term automatically, so there's one extra input)
# forward value is f(wx + b), where f is the activation function
# Σ := wx + b
# note: w is a parameter for the case of tied weights (it can be a TransposeView!)
type Layer{A <: Mapping, MATF <: AbstractMatrix{Float64}, GSTATE <: ParameterUpdaterState} <: NeuralNetLayer
  nin::Int
  nout::Int
  activation::A
  ploss::ParameterLoss
  p::Float64  # dropout retention probability

  # the state of the layer

  dwState::Matrix{GSTATE}
  dbState::Matrix{GSTATE}

  x::VecF  # nin x 1 -- input

  w::MATF  # nout x nin -- weights connecting previous layer to this layer
  b::VecF  # nout x 1 -- bias terms
  δ::VecF  # nout x 1 -- sensitivities (calculated during backward pass)
  Σ::VecF  # nout x 1 -- inner products (calculated during forward pass)
  a::VecF  # nout x 1 -- f(Σ) (calculated during forward pass)

  r::VecF  # nin x 1 -- vector of dropout retention... 0 if we drop this incoming weight, 1 if we keep it
  nextr::VecF  # nout x 1 -- retention of the nodes of this layer (as opposed to r which applies to the incoming weights)
end

function Layer(nin::Integer, nout::Integer, activation::Mapping,
               updater::ParameterUpdater, p::Float64 = 1.0;
               ploss::ParameterLoss = NoParameterLoss(),
               wgt = nothing, weightInit::Function = _initialWeights)
  w = weightInit(nin, nout, activation)
  Layer(nin, nout, activation, ploss, p,
            ParameterUpdaterState(updater, nout, nin),
            ParameterUpdaterState(updater, nout, 1),
            zeros(nin),
            w,
            [zeros(nout) for i in 1:4]...,
            ones(nin),
            ones(nout))
end

function Base.print{A,M,G}(io::IO, l::Layer{A,M,G})
  print(io, "Layer{$(l.nin)=>$(l.nout) $(l.activation) p=$(l.p) ‖δ‖₁=$(sumabs(l.δ)) $(M<:TransposeView ? "T" : "")}")
end

# gemv! :: Σ += p * w * x  (note: 'T' would imply p * w' * x)
@generated function dosigmamult!(layer::Layer, α::Float64)
  if layer.parameters[2] <: TransposeView
    return quote
      copy!(layer.Σ, layer.b)
      BLAS.gemv!('T', α, layer.w.mat, layer.x, 1.0, layer.Σ)
    end
  else
    return quote
      copy!(layer.Σ, layer.b)
      BLAS.gemv!('N', α, layer.w, layer.x, 1.0, layer.Σ)
    end
  end
end

# takes input vector, and computes Σⱼ = wⱼ'x + bⱼ  and  Oⱼ = A(Σⱼ)
function forward!(layer::Layer, x::AVecF, istraining::Bool)

  if istraining
    # train... randomly drop out incoming nodes
    # note: I said incoming, since this layers weights are the weights connecting the previous layer to this one
    #       So on dropout, we are actually dropping out thr previous layer's nodes...
    for i in 1:layer.nin
      ri = rand() <= layer.p ? 1.0 : 0.0
      layer.r[i] = ri
      layer.x[i] = ri * x[i]
    end

    dosigmamult!(layer, 1.0)
  else
    # test... need to multiply weights by dropout prob p
    fill!(layer.r, 1.0)
    copy!(layer.x, x)

    dosigmamult!(layer, layer.p)
  end

  value!(layer.a, layer.activation, layer.Σ)     # activate
end

# # backward step for the final (output) layer
# # note: costmult is the amount to multiply against f'(Σ)... L2 case should be: (yhat-y)
# function updateSensitivities!(layer::Layer, costmult::AVecF, multiplyDerivative::Bool)
#   copy!(layer.δ, costmult)
#   if multiplyDerivative
#     for i in 1:layer.nout
#       layer.δ[i] *= backward(layer.activation, layer.Σ[i])
#     end
#   end
# end

function updateSensitivities!(layer::Layer, mloss::ModelLoss, output::AVec, target::AVec)
    sensitivity!(layer.δ, layer.activation, mloss, layer.x, output, target)
end

# this is the backward step for a hidden layer
# notes: we are figuring out the effect of each node's activation value on the next sensitivities
# TODO: replace with LearnBase version when implemented
function updateSensitivities!(layer::Layer, nextlayer::Layer)
  for i in 1:layer.nout
    δi = 0.0
    for j in 1:nextlayer.nout
      δi += nextlayer.w[j,i] * nextlayer.nextr[j] * nextlayer.δ[j]
    end
    layer.δ[i] = δi * deriv(layer.activation, layer.Σ[i])
  end
end

# function param_change!(state::ParameterUpdaterState, updater::ParameterUpdater, ploss::ParameterLoss, gradient::Real, param::Real)

# update weights/bias one column at a time... skipping over the dropped out nodes
function updateWeights!(layer::Layer, updater::ParameterUpdater)

  # note: i refers to the output, j refers to the input

  for i in 1:layer.nout

    if layer.nextr[i] > 0.0

      # if this node is retained, we can update incoming bias
      bGrad = layer.δ[i]  # δi is the gradient
      dbi = param_change!(layer.dbState[i], updater, layer.ploss, bGrad, layer.b[i])
      layer.b[i] += dbi

      for j in 1:layer.nin

        # if this input node is retained, then we can also update the weight
        if layer.r[j] > 0.0
          wGrad = bGrad * layer.x[j]
          dwij = param_change!(layer.dwState[i,j], updater, layer.ploss, wGrad, layer.w[i,j])
          layer.w[i,j] += dwij
        end

      end
    end
  end
end


text:
#from https://github.com/jlustigy/ExoJulia/blob/master/ExoJulia/Orbit/keplerSolver.jl
function kepler_solver(M, ecc)
    Mred = mod(M,2pi)

    #Initial Guess:
    E = Mred + sign(sin(Mred))*0.85*ecc

    #Setting Initial Parameter
    di3 = 1.0

    #Setting the tolerance
    tol = 1e-12

    #Beginning the loop for quartic kepler_solver
    niter = 0
    while (abs(di3) > tol) & (niter < 30)
        SE = ecc*sin(E); CE = ecc*cos(E)
        f_of_E = E-SE-Mred; df_of_E=1-CE; d2f_of_E=SE; d3f_of_E=CE
        di1 = -f_of_E/df_of_E
        di2 = -f_of_E/(df_of_E+0.5*di1*d2f_of_E)
        di3 = -f_of_E/(df_of_E+0.5*di2*d2f_of_E+di2^2/6.0*d3f_of_E)

        E+=di3
        niter = niter + 1
    end

    E += (M-Mred)

    if(niter == 30)
        println("Error: Reached niter = 30")
    end

    E
end

#modified version of https://github.com/jlustigy/ExoJulia/blob/master/spring16/hw2/Agol_Agol/rv_model_one.jl
function analytic_rv(t, period, tp, ecc, m1, m2, a, w, inc)

    M = (t - tp)*2.0*pi/period

    eanom = kepler_solver(M, ecc)
    f = 2.0*atan(sqrt((1.0+ecc)/(1.0-ecc))*tan(0.5*eanom))

    G = NbodyGradient.GNEWT

    rvmod = ((sqrt(G/((m1+m2)*a*(1-ecc^2)))*m2*sin(inc))*(cos(w+f)+ecc*cos(w)))

end

function test_rvs()

    #defining the times of radial velocity (equivalent to how many data points)
    tmax = 40.0
    t = collect(0.0:1.0:tmax);

    #Parameters of Analytic Radial Velocity
    #period
    P = 17.5
    #time of periastron
    t_p = 0.0
    #eccentricity
    ecc = 0.001
    #star mass
    m_s = 1.034
    #planet mass
    m_p = 5.91e-5
    #semi-major axis
    a = 0.13340608724009007
    #argument of periastron
    ω = pi/2
    #inclination
    inc = pi/2

    #get radial velocities from analytic method
    rvs_analytic = analytic_rv.(t,P,t_p,ecc,m_s,m_p,a,ω,inc);

    #move to NbodyGradient method
    #Setting Initial Conditions for NbodyGradient Radial Velocity
    star = Elements(m = m_s)

    t0 = 0.0

    planet = Elements(
        m = m_p,
        P = P,
        t0 = t0,
        esinϖ = ecc*sin(pi+ω),
        ecosϖ = ecc*cos(pi+ω),
        Ω = float(pi/2),
        I = inc)

    N = 2 #number of bodies

    ic = ElementsIC(0.0,N,star,planet)
    s = State(ic)

    h = planet.P/30

    rvs = RadialVelocities(t,ic);

    #run Integrator
    intr = Integrator(h,tmax)
    intr(s,rvs)

    #get radial velocities from NbodyGradient method
    rvs_nbody = rvs.rvs;

    #test the accuracy between the analytic model and the NbodyGradient rv_model
    @test isapprox(rvs_analytic, rvs_nbody)
end


text:
module Preimpute

include("Regress.jl")

end


text:
using BumpStdlibs
using Test

import Dates
import Random
import TimeZones

@testset "BumpStdlibs.jl" begin
    include("unit-tests.jl")
    include("integration-tests.jl")
end


text:
module Vulnerability
include("BallisticLimit.jl")
using .BallisticLimit
end


text:
# ---
# title: 116. Populating Next Right Pointers in Each Node
# id: problem116
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Tree, Depth-first Search, Breadth-first Search
# link: <https://leetcode.com/problems/populating-next-right-pointers-in-each-node/description/>
# hidden: true
# ---
#
# You are given a **perfect binary tree**  where all leaves are on the same
# level, and every parent has two children. The binary tree has the following
# definition:
#
#
#
#     struct Node {
#       int val;
#       Node *left;
#       Node *right;
#       Node *next;
#     }
#
#
# Populate each next pointer to point to its next right node. If there is no
# next right node, the next pointer should be set to `NULL`.
#
# Initially, all next pointers are set to `NULL`.
#
#
#
# **Follow up:**
#
#   * You may only use constant extra space.
#   * Recursive approach is fine, you may assume implicit stack space does not count as extra space for this problem.
#
#
#
# **Example 1:**
#
# ![](https://assets.leetcode.com/uploads/2019/02/14/116_sample.png)
#
#
#
#     Input: root = [1,2,3,4,5,6,7]
#     Output: [1,#,2,3,#,4,5,6,7,#]
#     Explanation: Given the above perfect binary tree (Figure A), your function should populate each next pointer to point to its next right node, just like in Figure B. The serialized output is in level order as connected by the next pointers, with '#' signifying the end of each level.
#
#
#
#
# **Constraints:**
#
#   * The number of nodes in the given tree is less than `4096`.
#   * `-1000 <= node.val <= 1000`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
using Test

import LinearAlgebra.I
import LinearAlgebra.det
import LinearAlgebra.pinv
import LinearAlgebra.lu

include("../../../lib/KTBC.jl")
import .KTBC.CreateKTBC

# 14
K₃, = CreateKTBC(3)

# Calculate value
function makeSingularReduceDiagonal(M, step=.1, detV=.1)
    n, m = size(M)
    if (n != m)
        throw(DimensionMismatch("M must be square matrix!"))
    end
    while true
        M = M - I(n) * step;
        det(M) <= detV && break
    end
    return M
end

@test_throws DimensionMismatch makeSingularReduceDiagonal([
  1 2 3
  1 3 4
])

singularK₃ = makeSingularReduceDiagonal(K₃, .0000001, .0000001);
singularK₃

# 1.4142134996579785 <- very close!
singularK₃[1,1]
det(singularK₃)

# Solve the equation and give the exact answer √2
M = [
  √2 -1 0
  -1 √2 -1
  0 -1 √2
]

M[1,1]

det([
  √2 -1 0
  -1 √2 -1
  0 -1 √2
])

u = [
  1
  √2
  1
]

@test round.(Int, M * u) == [
  0
  0
  0
]

# 17
# Cat't solve in case of dependent cols. Only linear regression possible!
A = [
  2 3 1
  2 3 1
  3 2 1
]
invA = pinv(A)
b = [
  1
  2
  0
]

x = invA * b

A * x

@test A * x != b

# 18

A = [
  2 3
  4 5
]
B = [
  1 2
  2 4
]

A * B

# 21
A = [
  2 0
  3 1
]

B = [
  1 2
  4 5
]

@test (A * B) * (A * B) == (A * B)^2


text:
module DirectDepWithoutCompatEntry

import Example
import MLJModelInterface

f(x) = x + x

end # module


text:
# convert from Julia expression into NodeData form

function expr_to_nodedata(ex::Expr)
    nd = NodeData[]
    values = Float64[]
    expr_to_nodedata(ex,nd,values,-1)
    return nd,values
end

function expr_to_nodedata(ex::Expr,nd::Vector{NodeData},values::Vector{Float64},parentid)

    myid = length(nd) + 1
    if isexpr(ex,:call)
        op = ex.args[1]
        if length(ex.args) == 2
            push!(nd,NodeData(CALLUNIVAR, univariate_operator_to_id[op], parentid))
        elseif op in comparison_operators
            push!(nd,NodeData(COMPARISON, comparison_operator_to_id[op], parentid))
        else
            push!(nd,NodeData(CALL, operator_to_id[op], parentid))
        end
        for k in 2:length(ex.args)
            expr_to_nodedata(ex.args[k],nd,values,myid)
        end
    elseif isexpr(ex, :ref)
        @assert ex.args[1] == :x
        push!(nd,NodeData(VARIABLE, ex.args[2], parentid))
    elseif isexpr(ex, :comparison)
        op = ex.args[2]
        opid = comparison_operator_to_id[op]
        for k in 2:2:length(ex.args)-1
            @assert ex.args[k] == op
        end
        push!(nd, NodeData(COMPARISON, opid, parentid))
        for k in 1:2:length(ex.args)
            expr_to_nodedata(ex.args[k],nd,values,myid)
        end
    elseif isexpr(ex,:&&) || isexpr(ex,:||)
        @assert length(ex.args) == 2
        op = ex.head
        opid = logic_operator_to_id[op]
        push!(nd, NodeData(LOGIC, opid, parentid))
        expr_to_nodedata(ex.args[1],nd,values,myid)
        expr_to_nodedata(ex.args[2],nd,values,myid)
    else
        error("Unrecognized expression $ex: $(ex.head)")
    end
    nothing
end

function expr_to_nodedata(ex::Number,nd::Vector{NodeData},values::Vector{Float64},parentid)
    valueidx = length(values)+1
    push!(values,ex)
    push!(nd, NodeData(VALUE, valueidx, parentid))
    nothing
end

export expr_to_nodedata

# (i,j) nonzero means there's an edge *from* j to *i*
# since we get a column-oriented matrix, this gives us a fast way to look up the
# edges leaving any node (i.e., the children)
function adjmat(nd::Vector{NodeData})
    len = length(nd)
    I = Array(Int,len)
    J = Array(Int,len)
    realnz = 0
    for nz in 1:len
        par = nd[nz].parent
        par > 0 || continue
        realnz += 1
        I[realnz] = nz
        J[realnz] = par
    end
    resize!(I,realnz)
    resize!(J,realnz)
    return sparse(I,J,ones(Bool,realnz),len,len)
end

export adjmat


text:
using Test
using Altro
using BenchmarkTools
using TrajectoryOptimization
using RobotDynamics
using RobotZoo
using StaticArrays, LinearAlgebra
using ForwardDiff
using FiniteDiff
using Logging
const RD = RobotDynamics
const TO = TrajectoryOptimization

## Create a relatively large problem
N = 101
tf = 5.0
dt = tf / (N - 1)
model = RobotZoo.DoubleIntegrator(50)
n,m = size(model)
Q = Diagonal(@SVector fill(1.0, n))
R = Diagonal(@SVector fill(1e-2, m))
Qf = Q * (N - 1)
xf = @SVector fill(10.0, n)
x0 = @SVector zeros(n)
cons = ConstraintList(n,m,N)
obs = CircleConstraint(n, fill(4.0,50), rand(50)*10, fill(0.1,50))
goal = GoalConstraint(xf)
# unorm = NormConstraint(n, m, 1000, TO.SecondOrderCone(), :control)
# add_constraint!(cons, unorm, 1:N-1)
add_constraint!(cons, obs, 2:N-1)
add_constraint!(cons, goal, N)

obj = LQRObjective(Q, R, Qf, xf, N)
# prob = Problem(model, obj, x0, tf, xf = xf)
prob = Problem(model, obj, x0, tf, xf = xf, constraints=cons)

# Initialize the solver
t = @elapsed altrosolver = ALTROSolver(prob, dynamics_funsig=RD.InPlace())
@test t < 60  # it should finish initializing the solver in under a minute (usually about 8 seconds on a desktop)
t

##
altrosolver.opts.verbose = 2
altrosolver.opts.static_bp = false
altrosolver.opts.projected_newton = false
altrosolver.opts.trim_stats = false

solver = altrosolver.solver_al
conset = get_constraints(solver)
J_prev = cost(solver)

ilqr = Altro.get_ilqr(altrosolver)
# t_step = @elapsed Altro.step!(ilqr, J_prev)
t_step = @elapsed let solver = ilqr
    J_prev = cost(solver)
    Altro.errstate_jacobians!(solver.model, solver.G, solver.Z)
    Altro.dynamics_expansion!(solver)
    Altro.error_expansion!(solver.model, solver.D, solver.G)
    Altro.cost_expansion!(solver.obj, solver.Efull, solver.Z)
    Altro.error_expansion!(solver.model, solver.Eerr, solver.Efull, solver.G, solver.Z)

    # Get next iterate
    Altro.backwardpass!(solver)
    Jnew = Altro.forwardpass!(solver, J_prev)
end
@test t_step < 45*10  # Should be about 45 seconds.

##
t_solve = @elapsed solve!(altrosolver)
@test t_solve < 120  # should be close to 12 seconds
t_solve

## Projected Newton
pn = altrosolver.solver_pn
t_pn = @elapsed solve!(pn)
@test t_pn < 10*17  # should be about 17 seconds


text:
abstract type ClusterModel end

"""
    TotalClusterSS( Clustered::ClusterModel )

Returns a scalar of the total sum of squares for a ClusterModel object.
"""
function TotalClusterSS( Clustered::ClusterModel )
    GrandMean = StatsBase.mean( Clustered.X , dims = 1 )
    return sum((Clustered.X .- GrandMean).^ 2)
end

"""
    WithinClusterSS( Clustered::ClusterModel )

Returns a scalar of the within cluter sum of squares for a ClusterModel object.
"""
function WithinClusterSS( Clustered::ClusterModel )
    Clusters = unique( Clustered.Assignments )
    WithinSS = zeros( size(Clustered.X) )
    for cluster in Clusters
        ElementsInCluster = vec(Clustered.Assignments .== cluster)
        samples = Clustered.X[ ElementsInCluster, : ]
        clusterMean = StatsBase.mean( samples , dims = 1 )
        WithinSS[ElementsInCluster,:] = (samples .- clusterMean) .^ 2
    end
    return sum( WithinSS )
end

"""
    BetweenClusterSS( Clustered::ClusterModel )

Returns a scalar of the between cluster sum of squares for a ClusterModel object.
"""
function BetweenClusterSS( Clustered::ClusterModel )
    Clusters = unique( Clustered.Assignments )
    GrandMean = StatsBase.mean( Clustered.X , dims = 1 )
    BetweenSS = 0.0
    for cluster in Clusters
        ElementsInCluster = vec(Clustered.Assignments .== cluster)
        samples = Clustered.X[ ElementsInCluster, : ]
        clusterMean = StatsBase.mean( samples , dims = 1 )
        BetweenSS += size(samples)[1] * sum( ( clusterMean .- GrandMean ) .^ 2)
    end
    return BetweenSS
end

struct KMeansClustering <: ClusterModel
    X
    Centroids
    Assignments
end

"""
    KMeans( X, Clusters; tolerance = 1e-8, maxiters = 200 )

Returns a ClusterModel object after finding clusterings for data in `X` via MacQueens K-Means algorithm. `Clusters` is the K parameter, or the # of clusters.

MacQueen, J. B. (1967). Some Methods for classification and Analysis of Multivariate Observations. Proceedings of 5th Berkeley Symposium on Mathematical Statistics and Probability. 1. University of California Press. pp. 281–297.
"""
function KMeans( X, Clusters; tolerance = 1e-8, maxiters = 200 )
    (Xrows, Xcols) = size( X )
    ResultVector = zeros( Xrows )
    #1) Initialize K-Means centroids: Randomly select C samples
    rndsmpl = collect( StatsBase.sample( collect(1 : Xrows), Clusters, replace = false ) )
    Centroids = X[ rndsmpl, :]
    newCentroid = zeros( Clusters, Xcols )
    #2) finding which many points are nearest to our centroid
    #3) Calculate the mean value of each dimension to make a new centroid and repeat.
    for iter in 1 : maxiters
        DistsToCentroids = SquareEuclideanDistance( X, Centroids )
        ResultVector = last.(Tuple.(argmin( DistsToCentroids, dims = 2 )))

        for cluster in 1:Clusters
            Assigned = (ResultVector .== cluster)[:,1]
            newCentroid[ cluster, : ] = StatsBase.mean( X[ Assigned, : ], dims = 1 )
        end
        ToleranceCheck = ( 1.0 / Clusters ) * sum( ( Centroids .- newCentroid ) .^ 2)
        if ToleranceCheck <= tolerance ; break; end
        Centroids = newCentroid
    end
    return KMeansClustering( X, Centroids, ResultVector )
end


text:
for N in [Float32, Float64, Rational{Int}]
    # Example 6 from [1].
    #
    # [1] Tommaso Dreossi, Thao Dang, and Carla Piazza. *Reachability computation for polynomial dynamical systems.*
    #      Formal Methods in System Design 50.1 (2017): 1-38.
    D = N[-1 0 0; -1 -1 0; 0 0 -1]
    c = N[-0.80, -0.95, 0.0, 0.85, 1.0, 0.0]
    P = HParallelotope(D, c)

    # test getter functions
    @test directions(P) == D
    @test offset(P) == c
    @test dim(P) == 3

    # check constructor's size check
    @test_throws AssertionError HParallelotope(D, vcat(c, c))

    # test computation of the list of constraints
    @test constraints_list(P) == [HalfSpace(D[1, :], c[1]),
                                  HalfSpace(D[2, :], c[2]),
                                  HalfSpace(D[3, :], c[3]),
                                  HalfSpace(-D[1, :], c[4]),
                                  HalfSpace(-D[2, :], c[5]),
                                  HalfSpace(-D[3, :], c[6])]

    P = HParallelotope(N[1 0; 0 1], N[1, 1, 1, 1])

    # test center
    @test center(P) == N[0, 0]

    # test vertices functions
    @test base_vertex(P) == N[-1, -1]
    @test extremal_vertices(P) == [N[1, -1], N[-1, 1]]

    # test generators getters
    @test genmat(P) == N[1 0; 0 1]

    # random parallelotope
    rand(HParallelotope)
end

for N in [Float32, Float64]
    # conversion from zonotope
    Z = Zonotope(N[0, 0], N[1 0; 0 1])
    @test convert(HParallelotope, Z) ==
        HParallelotope(N[0 -1; 1 0], N[1, 1, 1, 1])
    Z = Zonotope(N[0, 0], N[1 0 1; 0 1 1])
    @test_throws AssertionError convert(HParallelotope, Z)
end

for N in [Float64]
    P = HParallelotope(N[1 0; 0 1], N[1, 1, 1, 1])
    @test collect(generators(P)) == [N[1, 0], N[0, 1]]
end


text:
#Test code
#-------------------------------------------------------------------------------

#No real test code yet... just run demos:

using Quantities
Q=Quantities

const AbstractLength{T<:Number, QT<:Q.AbstractLengthType, U<:Q.Unit} = Q.Quantity{T,QT,U}
const Length_m{T} = Q.Quantity{T,Q.LengthType,Q.Meter}

#Sample data:
#-------------------------------------------------------------------------------
l1m=Length(1.5, :m)
l2m=Length(1.2, :m)
d1m=Dist(5, :m)
d2m=Dist(8.01, :m)
l1in=Dist(55, :in)
d1in=Dist(83.01, :in)
inch=Length(1, :in)

#g(x::AbstractLength) = "CATCHALL: AbstractLength"
#==Problems:
 -Although above expression works, ::AbstractLength generates signatures that
  conflict with almost all declarations below
  NOTE: Not sure what the reason is for this, because you CAN do a proper
  catchall (see below).
==#

#Catchall that works:
g{T,QT<:Q.AbstractLengthType,U}(x::Q.Quantity{T,QT,U}) = "AbstractLength{ANY}"

g{T,QT<:Q.AbstractLengthType}(x::Q.Quantity{T,QT,Q.Meter}) = "AbstractLength{:m}"
g{T,QT<:Q.AbstractLengthType}(x::Q.Quantity{T,QT,Q.Inch}) = "AbstractLength{:in}"
g(x::Length_m) = "Length{:m}"
g{T}(x::_Length{T, Q.Unit{:in}}) = "_Length{:in}"
#g{T}(x::Q.Quantity{T, Q.LengthType, Q.Inch}) = "Length{:in}" #Equiv to previous

symbollist = [:l1m, :l2m, :d1m, :d2m, :l1in, :d1in]
for symb in symbollist
	val = eval(symb)
	print("$symb => ", )
		show(val)
		println()
end

valuelist = [eval(symb) for symb in symbollist]
for val in valuelist
	print("g($val) => ", g(val))
		println()
end

@show g(Length(3.3, :in))
@show g(Dist(3.8, :MYTYPE))

println("\nA few mathematical operations...")
@show d1m, d2m
@show d1m+d2m
@show d1m+d1m
@show 3*d1m
@show d1m*4.0
@show 5inch
@show 11.3inch
@show value(l1m)

:Test_Complete


text:
using Gadfly

set_default_plot_size(6inch, 8inch)

x = rand(40)
y = 4 .* x .+ 2 .+ randn(40)
plot(x=x,y=y, Geom.point, Geom.smooth(method=:lm))


text:
using NativeSVG

dr = Drawing(width = "300px", height = "100px") do
    desc() do
        str("Example InitialCoords - SVG's initial coordinate system")
    end
    g(fill = "none", stroke = "black", stroke_width = "3") do
        line(x1 = "0", y1 = "1.5", x2 = "300", y2 = "1.5")
        line(x1 = "1.5", y1 = "0", x2 = "1.5", y2 = "100")
    end
    g(fill = "red", stroke = "none") do
        rect(x = "0", y = "0", width = "3", height = "3")
        rect(x = "297", y = "0", width = "3", height = "3")
        rect(x = "0", y = "97", width = "3", height = "3")
    end
    g(font_size = "14", font_family = "Verdana") do
        text(x = "10", y = "20") do
            str("(0,0)")
        end
        text(x = "240", y = "20") do
            str("(300,0)")
        end
        text(x = "10", y = "90") do
            str("(0,100)")
        end
    end
end
display(dr)

dr =
    Drawing(
        width = "300px",
        height = "200px",
        viewBox = "0 0 1500 1000",
        preserveAspectRatio = "none"
    ) do
        desc() do
            str("""Example ViewBox - uses the viewBox
                   attribute to automatically create an initial user coordinate
                   system which causes the graphic to scale to fit into the
                   SVG viewport no matter what size the SVG viewport is.""")
        end
    # This rectangle goes from (0,0) to (1500,1000) in user coordinate system.
    # Because of the viewBox attribute above,
    # the rectangle will end up filling the entire area
    # reserved for the SVG content.
        rect(
            x = "0",
            y = "0",
            width = "1500",
            height = "1000",
            fill = "yellow",
            stroke = "blue",
            stroke_width = "12"
        )
    # A large, red triangle
        path(fill = "red", d = "M 750,100 L 250,900 L 1250,900 z")
    # A text string that spans most of the SVG viewport
        text(x = "100", y = "600", font_size = "200", font_family = "Verdana") do
            str("Stretch to fit")
        end
    end
display(dr)

max_stars_repo_name: BenLauwens/NativeSVG

max_issues_repo_name: BenLauwens/NativeSVG

max_forks_repo_name: BenLauwens/NativeSVG