text:
# Controls whether tuples and non-pointerfree immutables, which Julia
# stores as references, are stored inline in compound types when
# possible. Currently this is problematic because Julia fields of these
# types may be undefined.
const INLINE_TUPLE = false
const INLINE_POINTER_IMMUTABLE = false

const JLD_REF_TYPE = JldDatatype(HDF5Datatype(HDF5.H5T_STD_REF_OBJ, false), 0)
const BUILTIN_TYPES = Set([Symbol, Type, BigFloat, BigInt])
const JL_TYPENAME_TRANSLATE = Dict{String, String}()
const JLCONVERT_INFO = Dict{Any, Any}()
const H5CONVERT_INFO = Dict{Any, Any}()

const EMPTY_TUPLE_TYPE = Tuple{}
const TypesType = Core.SimpleVector
TupleType{T<:Tuple} = Type{T}
tupletypes(T::TupleType) = T.parameters
typetuple(types) = Tuple{types...}

## Helper functions

translate(oldname::String, newname::String) = JL_TYPENAME_TRANSLATE[oldname] = newname

# Holds information about the mapping between a Julia and HDF5 type
struct JldTypeInfo
    dtypes::Vector{JldDatatype}
    offsets::Vector{Int}
    size::Int
end

# Get information about the HDF5 types corresponding to Julia types
function JldTypeInfo(parent::JldFile, types::TypesType, commit::Bool)
    dtypes = Vector{JldDatatype}(undef, length(types))
    offsets = Vector{Int}(undef, length(types))
    offset = 0
    for i = 1:length(types)
        dtype = dtypes[i] = h5fieldtype(parent, types[i], commit)
        offsets[i] = offset
        offset += HDF5.h5t_get_size(dtype)
    end
    JldTypeInfo(dtypes, offsets, offset)
end
JldTypeInfo(parent::JldFile, @nospecialize(T), commit::Bool) =
    JldTypeInfo(parent, Base.unwrap_unionall(T).types, commit)

# Write an HDF5 datatype to the file
function commit_datatype(parent::JldFile, dtype::HDF5Datatype, @nospecialize(T))
    pparent = parent.plain
    if !exists(pparent, pathtypes)
        gtypes = g_create(pparent, pathtypes)
    else
        gtypes = pparent[pathtypes]
    end

    id = length(gtypes)+1
    try
        HDF5.t_commit(gtypes, @sprintf("%08d", id), dtype)
    finally
        close(gtypes)
    end
    a_write(dtype, name_type_attr, full_typename(parent, T))

    # Store in map
    parent.jlh5type[T] = JldDatatype(dtype, id)
end

# If parent is nothing, we are creating the datatype in memory for
# validation, so don't commit it
commit_datatype(parent::Nothing, dtype::HDF5Datatype, @nospecialize(T)) =
    JldDatatype(dtype, -1)

# The HDF5 library loses track of relationships among committed types
# after the file is saved. We mangle the names by appending a
# sequential identifier so that we can recover these relationships
# later.
mangle_name(jtype::JldDatatype, jlname) =
    jtype.index <= 0 ? string(jlname, "_") : string(jlname, "_", jtype.index)
Base.convert(::Type{HDF5.Hid}, x::JldDatatype) = x.dtype.id

## Serialization of datatypes to JLD
##
## h5fieldtype - gets the JldDatatype corresponding to a given
## Julia type, when the Julia type is stored as an element of an HDF5
## compound type or array. This is the only function that can operate
## on non-concrete types.
##
## h5type - gets the JldDatatype corresponding to an object of the
## given Julia type. For pointerfree types, this is usually the same as
## the h5fieldtype.
##
## h5convert! - converts data from Julia to HDF5 in a buffer. Most
## methods are dynamically generated by gen_h5convert, but methods for
## special built-in types are predefined.
##
## jlconvert - converts data from HDF5 to a Julia object.
##
## jlconvert! - converts data from HDF5 to Julia in a buffer. This is
## only applicable in cases where fields of that type may not be stored
## as references (e.g., not plain types).

## Special types
##
## To create a special serialization of a datatype, one should:
##
## - Define a method of h5fieldtype that dispatches to h5type
## - Define a method of h5type that constructs the type
## - Define h5convert! and jlconvert
## - If the type is an immutable, define jlconvert!
## - Add the type to BUILTIN_TYPES

## HDF5 bits kinds

# This construction prevents these methods from getting called on type unions
const BitsKindTypes = Union{map(x->Type{x}, Base.uniontypes(HDF5.HDF5BitsKind))...}

h5fieldtype(parent::JldFile, T::BitsKindTypes, ::Bool) =
    h5type(parent, T, false)

h5type(::JldFile, T::BitsKindTypes, ::Bool) =
    JldDatatype(HDF5Datatype(HDF5.hdf5_type_id(T), false), 0)

h5convert!(out::Ptr, ::JldFile, x::T, ::JldWriteSession) where {T<:HDF5.HDF5BitsKind} =
    unsafe_store!(convert(Ptr{T}, out), x)

_jlconvert_bits(::Type{T}, ptr::Ptr) where {T} = unsafe_load(convert(Ptr{T}, ptr))
_jlconvert_bits!(out::Ptr, ::Type{T}, ptr::Ptr) where {T} =
    (unsafe_store!(convert(Ptr{T}, out), unsafe_load(convert(Ptr{T}, ptr))); nothing)

jlconvert(T::BitsKindTypes, ::JldFile, ptr::Ptr) = _jlconvert_bits(T, ptr)
jlconvert!(out::Ptr, T::BitsKindTypes, ::JldFile, ptr::Ptr) = _jlconvert_bits!(out, T, ptr)

function jlconvert(T::Type{Char}, f::JldFile, ptr::Ptr)
    c = unsafe_load(convert(Ptr{Char}, ptr))
    if f.version < v"0.1.2"
        c = Char(reinterpret(UInt32, c))
    end
    return c
end
function jlconvert!(out::Ptr, T::Type{Char}, f::JldFile, ptr::Ptr)
    unsafe_store!(convert(Ptr{Char}, out), jlconvert(Char, f, ptr))
end

## Nothing

const NothingType = Type{Nothing}

jlconvert(T::NothingType, ::JldFile, ptr::Ptr) = nothing
jlconvert!(out::Ptr, T::NothingType, ::JldFile, ptr::Ptr) = (unsafe_store!(convert(Ptr{T}, out), nothing); nothing)
h5convert!(out::Ptr, ::JldFile, x::Nothing, ::JldWriteSession) = nothing

## Strings

h5fieldtype(parent::JldFile, ::Type{T}, ::Bool) where {T<:String} =
    h5type(parent, T, false)

# Stored as variable-length strings
function h5type(::JldFile, ::Type{T}, ::Bool) where T<:String
    type_id = HDF5.h5t_copy(HDF5.hdf5_type_id(T))
    HDF5.h5t_set_size(type_id, HDF5.H5T_VARIABLE)
    HDF5.h5t_set_cset(type_id, HDF5.cset(T))
    JldDatatype(HDF5Datatype(type_id, false), 0)
end

# no-inline needed to ensure gc-root for x
@noinline h5convert!(out::Ptr, ::JldFile, x::String, ::JldWriteSession) =
    unsafe_store!(convert(Ptr{Ptr{UInt8}}, out), pointer(x))

function jlconvert(T::Union{Type{String}}, ::JldFile, ptr::Ptr)
    strptr = unsafe_load(convert(Ptr{Ptr{UInt8}}, ptr))
    str = unsafe_string(strptr)
    Libc.free(strptr)
    str
end

## Symbols

h5fieldtype(parent::JldFile, ::Type{Symbol}, commit::Bool) =
    h5type(parent, Symbol, commit)

# Stored as a compound type that contains a variable length string
function h5type(parent::JldFile, ::Type{Symbol}, commit::Bool)
    haskey(parent.jlh5type, Symbol) && return parent.jlh5type[Symbol]
    id = HDF5.h5t_create(HDF5.H5T_COMPOUND, 8)
    HDF5.h5t_insert(id, "symbol_", 0, h5fieldtype(parent, String, commit))
    dtype = HDF5Datatype(id, parent.plain)
    commit ? commit_datatype(parent, dtype, Symbol) : JldDatatype(dtype, -1)
end

function h5convert!(out::Ptr, file::JldFile, x::Symbol, wsession::JldWriteSession)
    str = string(x)
    push!(wsession.persist, str)
    h5convert!(out, file, str, wsession)
end

jlconvert(::Type{Symbol}, file::JldFile, ptr::Ptr) = Symbol(jlconvert(String, file, ptr))

## BigInts and BigFloats

h5fieldtype(parent::JldFile, T::Union{Type{BigInt}, Type{BigFloat}}, commit::Bool) =
    h5type(parent, T, commit)

# Stored as a compound type that contains a variable length string
function h5type(parent::JldFile, T::Union{Type{BigInt}, Type{BigFloat}}, commit::Bool)
    haskey(parent.jlh5type, T) && return parent.jlh5type[T]
    id = HDF5.h5t_create(HDF5.H5T_COMPOUND, 8)
    HDF5.h5t_insert(id, "data_", 0, h5fieldtype(parent, String, commit))
    dtype = HDF5Datatype(id, parent.plain)
    commit ? commit_datatype(parent, dtype, T) : JldDatatype(dtype, -1)
end

function h5convert!(out::Ptr, file::JldFile, x::BigInt, wsession::JldWriteSession)
    str = string(x, base=62)
    push!(wsession.persist, str)
    h5convert!(out, file, str, wsession)
end
function h5convert!(out::Ptr, file::JldFile, x::BigFloat, wsession::JldWriteSession)
    str = string(x)
    push!(wsession.persist, str)
    h5convert!(out, file, str, wsession)
end

jlconvert(::Type{BigInt}, file::JldFile, ptr::Ptr) =
    parse(BigInt, jlconvert(String, file, ptr), base = 62)
jlconvert(::Type{BigFloat}, file::JldFile, ptr::Ptr) =
    parse(BigFloat, jlconvert(String, file, ptr))

## Types

h5fieldtype(parent::JldFile, ::Type{T}, commit::Bool) where {T<:Type} =
    h5type(parent, Type, commit)

# Stored as a compound type that contains a variable length string
function h5type(parent::JldFile, ::Type{T}, commit::Bool) where T<:Type
    haskey(parent.jlh5type, Type) && return parent.jlh5type[Type]
    id = HDF5.h5t_create(HDF5.H5T_COMPOUND, 8)
    HDF5.h5t_insert(id, "typename_", 0, h5fieldtype(parent, String, commit))
    dtype = HDF5Datatype(id, parent.plain)
    out = commit ? commit_datatype(parent, dtype, Type) : JldDatatype(dtype, -1)
end

function h5convert!(out::Ptr, file::JldFile, x::Type, wsession::JldWriteSession)
    str = full_typename(file, x)
    push!(wsession.persist, str)
    h5convert!(out, file, str, wsession)
end

jlconvert(::Type{T}, file::JldFile, ptr::Ptr) where {T<:Type} =
    julia_type(jlconvert(String, file, ptr))

## Union{}

h5fieldtype(parent::JldFile, ::Type{Union{}}, ::Bool) = JLD_REF_TYPE

## Arrays

# These show up as having T.size == 0, hence the need for
# specialization.
h5fieldtype(parent::JldFile, ::Type{Array{T,N}}, ::Bool) where {T,N} = JLD_REF_TYPE

## User-defined types
##
## Similar to special types, but h5convert!/jl_convert are dynamically
## generated.

## Tuples

if INLINE_TUPLE
    h5fieldtype(parent::JldFile, T::TupleType, commit::Bool) =
        isconcretetype(T) ? h5type(parent, T, commit) : JLD_REF_TYPE
else
    h5fieldtype(parent::JldFile, T::TupleType, ::Bool) = JLD_REF_TYPE
end

function h5type(parent::JldFile, T::TupleType, commit::Bool)
    haskey(parent.jlh5type, T) && return parent.jlh5type[T]
    # Tuples should always be concretely typed, unless we're
    # reconstructing a tuple, in which case commit will be false
    !commit || isconcretetype(T) || error("unexpected non-concrete type $T")

    typeinfo = JldTypeInfo(parent, T, commit)
    if isopaque(T)
        id = HDF5.h5t_create(HDF5.H5T_OPAQUE, opaquesize(T))
    else
        id = HDF5.h5t_create(HDF5.H5T_COMPOUND, typeinfo.size)
    end
    for i = 1:length(typeinfo.offsets)
        fielddtype = typeinfo.dtypes[i]
        HDF5.h5t_insert(id, mangle_name(fielddtype, i), typeinfo.offsets[i], fielddtype)
    end

    dtype = HDF5Datatype(id, parent.plain)
    if commit
        jlddtype = commit_datatype(parent, dtype, T)
        if T == EMPTY_TUPLE_TYPE
            # to allow recovery of empty tuples, which HDF5 does not allow
            a_write(dtype, "empty", UInt8(1))
        end
        jlddtype
    else
        JldDatatype(dtype, -1)
    end
end

## All other objects

# For cases not defined above: If the type is mutable and non-empty,
# this is a reference. If the type is immutable, this is a type itself.
if INLINE_POINTER_IMMUTABLE
    h5fieldtype(parent::JldFile, @nospecialize(T), commit::Bool) =
        isconcretetype(T) && (!T.mutable || T.size == 0) ? h5type(parent, T, commit) : JLD_REF_TYPE
else
    h5fieldtype(parent::JldFile, @nospecialize(T), commit::Bool) =
        isconcretetype(T) && (!T.mutable || T.size == 0) && datatype_pointerfree(T) ? h5type(parent, T, commit) : JLD_REF_TYPE
end

function h5type(parent::JldFile, T::Type{Bool}, commit::Bool)
    if parent.version < v"0.1.3"
        h5type_default(parent, T, commit)
    else
        JldDatatype(HDF5Datatype(HDF5.hdf5_type_id(T), false), 0)
    end
end

h5type(parent::JldFile, @nospecialize(T), commit::Bool) = h5type_default(parent, T, commit)

function h5type_default(parent::JldFile, @nospecialize(T), commit::Bool)
    !isa(T, DataType) && unknown_type_err(T)
    T = T::DataType

    haskey(parent.jlh5type, T) && return parent.jlh5type[T]
    isconcretetype(T) || error("unexpected non-concrete type ", T)

    if isopaque(T)
        # Empty type or non-basic bitstype
        id = HDF5.h5t_create(HDF5.H5T_OPAQUE, opaquesize(T))
    else
        # Compound type
        typeinfo = JldTypeInfo(parent, T.types, commit)
        id = HDF5.h5t_create(HDF5.H5T_COMPOUND, typeinfo.size)
        for i = 1:length(typeinfo.offsets)
            fielddtype = typeinfo.dtypes[i]
            HDF5.h5t_insert(id, mangle_name(fielddtype, fieldnames(T)[i]), typeinfo.offsets[i], fielddtype)
        end
    end

    dtype = HDF5Datatype(id, parent.plain)
    if commit
        jlddtype = commit_datatype(parent, dtype, T)
        if T.size == 0
            # to allow recovery of empty types, which HDF5 does not allow
            a_write(dtype, "empty", UInt8(1))
        end
        jlddtype
    else
        JldDatatype(dtype, -1)
    end
end

# Normal objects
function _gen_jlconvert_type(typeinfo::JldTypeInfo, @nospecialize(T))
    ex = Expr(:block)
    args = ex.args
    push!(args, :(out = ccall(:jl_new_struct_uninit, Ref{T}, (Any,), T)))
    for i = 1:length(typeinfo.dtypes)
        h5offset = typeinfo.offsets[i]

        if HDF5.h5t_get_class(typeinfo.dtypes[i]) == HDF5.H5T_REFERENCE
            push!(args, quote
                ref = unsafe_load(convert(Ptr{HDF5ReferenceObj}, ptr)+$h5offset)
                if ref != HDF5.HDF5ReferenceObj_NULL
                    out.$(fieldnames(T)[i]) = convert($(T.types[i]), read_ref(file, ref))
                end
            end)
        else
            push!(args, :(out.$(fieldnames(T)[i]) = jlconvert($(T.types[i]), file, ptr+$h5offset)))
        end
    end
    push!(args, :(return out))
    return ex
end

function _gen_jlconvert_type!(typeinfo::JldTypeInfo, @nospecialize(T))
    error("unimplemented")
end

# Immutables
function _gen_jlconvert_immutable(typeinfo::JldTypeInfo, @nospecialize(T))
    ex = Expr(:block)
    args = ex.args
    jloffsets = map(idx->fieldoffset(T, idx), 1:fieldcount(T))
    if isbitstype(T)
        push!(args, :(out = Ref{T}()))
        push!(args, :(jlconvert!(unsafe_convert(Ptr{T}, out), T, file, ptr)))
        push!(args, :(return out[]))
    else
        nf = length(typeinfo.dtypes)
        push!(args, :(fieldvals = Vector{Any}(undef, $nf)))
        push!(args, :(ninit = 0))
        for i = 1:nf
            h5offset = typeinfo.offsets[i]
            jloffset = jloffsets[i]
            obj = gensym("obj")
            if isa(T.types[i], TupleType) && isbitstype(T.types[i])
                # We continue to store tuples as references for the sake of
                # backwards compatibility, but on 0.4 they are now stored
                # inline
                push!(args, quote
                    ref = unsafe_load(convert(Ptr{HDF5ReferenceObj}, ptr)+$h5offset)
                    if ref == HDF5.HDF5ReferenceObj_NULL
                        @warn("""A pointerfree tuple field was undefined.
                                 This is not supported in Julia 0.4 and the corresponding tuple will be uninitialized.""")
                    else
                        fieldvals[$i] = convert($(T.types[i]), read_ref(file, ref))
                        ninit += 1
                    end
                end)
            elseif HDF5.h5t_get_class(typeinfo.dtypes[i]) == HDF5.H5T_REFERENCE
                push!(args, quote
                    ref = unsafe_load(convert(Ptr{HDF5ReferenceObj}, ptr)+$h5offset)
                    if ref != HDF5.HDF5ReferenceObj_NULL
                        fieldvals[$i] = convert($(T.types[i]), read_ref(file, ref))
                        ninit += 1
                    end
                end)
            else
                push!(args, :(fieldvals[$i] = jlconvert($(T.types[i]), file, ptr+$h5offset)))
                push!(args, :(ninit += 1))
            end
        end
        push!(args, :(out = ccall(:jl_new_structv, Ref{T}, (Any,Ptr{Cvoid},UInt32), T, fieldvals, ninit)))
        push!(args, :(return out))
    end
    return ex
end

function _gen_jlconvert_immutable!(typeinfo::JldTypeInfo, @nospecialize(T))
    ex = Expr(:block)
    args = ex.args
    jloffsets = map(idx->fieldoffset(T, idx), 1:fieldcount(T))
    if isbitstype(T)
        for i = 1:length(typeinfo.dtypes)
            h5offset = typeinfo.offsets[i]
            jloffset = jloffsets[i]

            if isa(T.types[i], TupleType) && isbitstype(T.types[i])
                # We continue to store tuples as references for the sake of
                # backwards compatibility, but on 0.4 they are now stored
                # inline
                push!(args, quote
                    ref = unsafe_load(convert(Ptr{HDF5ReferenceObj}, ptr)+$h5offset)
                    if ref == HDF5.HDF5ReferenceObj_NULL
                        @warn("""A pointerfree tuple field was undefined.
                                 This is not supported in Julia 0.4 and the corresponding tuple will be uninitialized.""")
                    else
                        unsafe_store!(convert(Ptr{$(T.types[i])}, out)+$jloffset, read_ref(file, ref))
                    end
                end)
            elseif HDF5.h5t_get_class(typeinfo.dtypes[i]) == HDF5.H5T_REFERENCE
                error("reference encountered in pointerfree immutable; this is a bug")
            else
                push!(args, :(jlconvert!(out+$jloffset, $(T.types[i]), file, ptr+$h5offset)))
            end
        end
    else
        error("unimplemented")
    end
    push!(args, nothing)
    return ex
end

function _gen_jlconvert_tuple(typeinfo::JldTypeInfo, @nospecialize(T))
    # custom converter (rather than falling back to immutable datatype handling)
    # as we need to make sure to return a valid Tuple concrete type,
    # but some parameters of T may be Any
    ex = Expr(:block)
    args = ex.args
    tup = Expr(:tuple)
    tupargs = tup.args
    types = tupletypes(T)
    for i = 1:length(typeinfo.dtypes)
        h5offset = typeinfo.offsets[i]
        field = Symbol(string("field", i))

        if HDF5.h5t_get_class(typeinfo.dtypes[i]) == HDF5.H5T_REFERENCE
            push!(args, :($field = read_ref(file, unsafe_load(convert(Ptr{HDF5ReferenceObj}, ptr)+$h5offset))))
        else
            push!(args, :($field = jlconvert($(types[i]), file, ptr+$h5offset)))
        end
        push!(tupargs, field)
    end
    push!(args, :(return $tup))
    return ex
end

function gen_jlconvert(typeinfo::JldTypeInfo, @nospecialize(T))
    T === Nothing && return
    # TODO: this is probably invalid, so try to do this differently
    JLCONVERT_INFO[T] = typeinfo
    nothing
end

function gen_jlconvert(@nospecialize(T))
    typeinfo = JLCONVERT_INFO[T]::JldTypeInfo
    if isa(T, TupleType)
        return _gen_jlconvert_tuple(typeinfo, T)
    elseif isempty(fieldnames(T))
        if T.size == 0 && !T.mutable
            return T.instance
        else
           return :(_jlconvert_bits(T, ptr))
        end
    elseif T.size == 0
        return :(ccall(:jl_new_struct_uninit, Ref{T}, (Any,), T))
    elseif T.mutable
        return _gen_jlconvert_type(typeinfo, T)
    else
        return _gen_jlconvert_immutable(typeinfo, T)
    end
end

function gen_jlconvert!(@nospecialize(T))
    typeinfo = JLCONVERT_INFO[T]::JldTypeInfo
    if isa(T, TupleType)
        error("unimplemented")
    elseif isempty(fieldnames(T))
        if T.size == 0
            if !T.mutable
                return nothing
            else
                return :(unsafe_store!(convert(Ptr{Any}, out), $(T.instance)); nothing)
            end
        else
           return :(_jlconvert_bits!(out, T, ptr); nothing)
        end
    elseif T.size == 0
        return nothing
    elseif T.mutable
        return _gen_jlconvert_type!(typeinfo, T)
    else
        return _gen_jlconvert_immutable!(typeinfo, T)
    end
end

@generated function jlconvert!(out::Ptr, ::Type{T}, file::JldFile, ptr::Ptr) where T
    return gen_jlconvert!(T)
end

@generated function jlconvert(::Type{T}, file::JldFile, ptr::Ptr) where T
    return gen_jlconvert(T)
end

## Common functions for all non-special types (including gen_h5convert)

# Whether this datatype should be stored as opaque
isopaque(t::TupleType) = t == EMPTY_TUPLE_TYPE
# isopaque(t::DataType) = isempty(fieldnames(t))
isopaque(t::DataType) = isa(t, TupleType) ? t == EMPTY_TUPLE_TYPE : isempty(fieldnames(t))

# The size of this datatype in the HDF5 file (if opaque)
opaquesize(t::TupleType) = 1
opaquesize(t::DataType) = max(1, t.size)

# Whether a type that is stored inline in HDF5 should be stored as a
# reference in Julia. This will only be called such that it returns
# true for some unions of special types defined above, unless either
# INLINE_TUPLE or INLINE_POINTER_IMMUTABLE is true.
uses_reference(T::DataType) = !datatype_pointerfree(T)
uses_reference(::TupleType) = true
uses_reference(::Union) = true

unknown_type_err(T) =
    error("""$T is not of a type supported by JLD
             Please report this error at https://github.com/JuliaIO/HDF5.jl""")

const BUILTIN_H5_types = Union{Nothing, Type, String, HDF5.HDF5BitsKind, Symbol, BigInt, BigFloat}
function gen_h5convert(parent::JldFile, @nospecialize(T))
    T <: BUILTIN_H5_types && return
    # TODO: this is probably invalid, so try to do this differently
    haskey(H5CONVERT_INFO, T) && return
    H5CONVERT_INFO[T] = parent

    dtype = parent.jlh5type[T].dtype
    istuple = isa(T, TupleType)

    if isopaque(T)
        return
    end

    if istuple
        types = tupletypes(T)
    else
        types = (T::DataType).types
    end

    n = HDF5.h5t_get_nmembers(dtype.id)
    for i = 1:n
        if HDF5.h5t_get_member_class(dtype.id, i-1) != HDF5.H5T_REFERENCE
            gen_h5convert(parent, types[i])
        end
    end
    nothing
end

# There is no point in specializing this
function _gen_h5convert!(@nospecialize(T))
    parent = H5CONVERT_INFO[T]::JldFile
    dtype = parent.jlh5type[T].dtype
    istuple = isa(T, TupleType)

    if isopaque(T)
        if T.size == 0
            return nothing
        else
            return :(unsafe_store!(convert(Ptr{T}, out), x))
        end
    end

    if istuple
        types = tupletypes(T)
    else
        types = (T::DataType).types
    end

    getindex_fn = istuple ? (:getindex) : (:getfield)
    ex = Expr(:block)
    args = ex.args
    n = HDF5.h5t_get_nmembers(dtype.id)
    for i = 1:n
        offset = HDF5.h5t_get_member_offset(dtype.id, i-1)
        if HDF5.h5t_get_member_class(dtype.id, i-1) == HDF5.H5T_REFERENCE
            if istuple
                push!(args, :(unsafe_store!(convert(Ptr{HDF5ReferenceObj}, out)+$offset,
                                            write_ref(file, $getindex_fn(x, $i), wsession))))
            else
                push!(args, quote
                    if isdefined(x, $i)
                        ref = write_ref(file, $getindex_fn(x, $i), wsession)
                    else
                        ref = HDF5.HDF5ReferenceObj_NULL
                    end
                    unsafe_store!(convert(Ptr{HDF5ReferenceObj}, out)+$offset, ref)
                end)
            end
        else
            push!(args, :(h5convert!(out+$offset, file, $getindex_fn(x, $i), wsession)))
        end
    end
    push!(args, nothing)
    return ex
end

@generated function h5convert!(out::Ptr, file::JldFile, x::T, wsession::JldWriteSession) where T
    return _gen_h5convert!(T)
end

## Find the corresponding Julia type for a given HDF5 type

# Type mapping function. Given an HDF5Datatype, find (or construct) the
# corresponding Julia type.
function jldatatype(parent::JldFile, dtype::HDF5Datatype)
    class_id = HDF5.h5t_get_class(dtype.id)
    if class_id == HDF5.H5T_STRING
        cset = HDF5.h5t_get_cset(dtype.id)
        if cset == HDF5.H5T_CSET_ASCII
            return String
        elseif cset == HDF5.H5T_CSET_UTF8
            return String
        else
            error("character set ", cset, " not recognized")
        end
    elseif class_id == HDF5.H5T_INTEGER || class_id == HDF5.H5T_FLOAT || class_id == HDF5.H5T_BITFIELD
        # This can be a performance hotspot
        HDF5.h5t_equal(dtype.id, HDF5.H5T_NATIVE_DOUBLE) > 0 && return Float64
        HDF5.h5t_equal(dtype.id, HDF5.H5T_NATIVE_INT64) > 0 && return Int64
        HDF5.h5t_equal(dtype.id, HDF5.H5T_NATIVE_FLOAT) > 0 && return Float32
        HDF5.h5t_equal(dtype.id, HDF5.H5T_NATIVE_INT32) > 0 && return Int32
        HDF5.h5t_equal(dtype.id, HDF5.H5T_NATIVE_UINT8) > 0 && return UInt8
        HDF5.h5t_equal(dtype.id, HDF5.H5T_NATIVE_UINT64) > 0 && return UInt64
        HDF5.h5t_equal(dtype.id, HDF5.H5T_NATIVE_UINT32) > 0 && return UInt32
        HDF5.h5t_equal(dtype.id, HDF5.H5T_NATIVE_INT8) > 0 && return Int8
        HDF5.h5t_equal(dtype.id, HDF5.H5T_NATIVE_INT16) > 0 && return Int16
        HDF5.h5t_equal(dtype.id, HDF5.H5T_NATIVE_UINT16) > 0 && return UInt16
        HDF5.h5t_equal(dtype.id, HDF5.H5T_NATIVE_B8) > 0 && return Bool
        error("unrecognized integer or float type")
    elseif class_id == HDF5.H5T_BITFIELD
        Bool
    elseif class_id == HDF5.H5T_COMPOUND || class_id == HDF5.H5T_OPAQUE
        addr = HDF5.objinfo(dtype).addr
        haskey(parent.h5jltype, addr) && return parent.h5jltype[addr]

        typename = a_read(dtype, name_type_attr)
        typename = get(JL_TYPENAME_TRANSLATE, typename, typename)
        T = julia_type(typename)
        if T == UnsupportedType
            @warn("type $typename not present in workspace; reconstructing")
            T = reconstruct_type(parent, dtype, typename)
        end

        T = readas(T)

        if !(T in BUILTIN_TYPES)
            # Call jldatatype on dependent types to validate them and
            # define jlconvert
            if class_id == HDF5.H5T_COMPOUND
                for i = 0:HDF5.h5t_get_nmembers(dtype.id)-1
                    member_name = HDF5.h5t_get_member_name(dtype.id, i)
                    idx = first(something(findlast("_", member_name), 0:-1))
                    if idx != sizeof(member_name)
                        member_dtype = HDF5.t_open(parent.plain, string(pathtypes, '/', lpad(member_name[idx+1:end], 8, '0')))
                        jldatatype(parent, member_dtype)
                    end
                end
            end

            gen_jlconvert(JldTypeInfo(parent, T, false), T)
        end

        # Verify that types match
        newtype = h5type(parent, T, false).dtype
        if T !== Expr  # in julia >= 0.7 Expr has 1 fewer fields; the trailing field can be ignored.
            dtype == newtype || throw(TypeMismatchException(typename))
        end

        # Store type in type index
        index = typeindex(parent, addr)
        parent.jlh5type[T] = JldDatatype(dtype, index)
        parent.h5jltype[addr] = T
        T
    else
        error("unrecognized HDF5 datatype class ", class_id)
    end
end

# Create a Julia type based on the HDF5Datatype from the file. Used
# when the type is no longer available.
function reconstruct_type(parent::JldFile, dtype::HDF5Datatype, savedname::AbstractString)
    name = gensym(savedname)
    class_id = HDF5.h5t_get_class(dtype.id)
    if class_id == HDF5.H5T_OPAQUE
        if exists(dtype, "empty")
            @eval (struct $name; end; $name)
        else
            sz = Int(HDF5.h5t_get_size(dtype.id))*8
            @eval (primitive type $name $sz end; $name)
        end
    else
        # Figure out field names and types
        nfields = HDF5.h5t_get_nmembers(dtype.id)
        fieldnames = Vector{Symbol}(undef, nfields)
        fieldtypes = Vector{Type}(undef, nfields)
        for i = 1:nfields
            membername = HDF5.h5t_get_member_name(dtype.id, i-1)
            idx = first(something(findlast("_", membername), 0:-1))
            fieldname = fieldnames[i] = Symbol(membername[1:idx-1])

            if idx != sizeof(membername)
                # There is something past the underscore in the HDF5 field
                # name, so the type is stored in file
                memberdtype = HDF5.t_open(parent.plain, string(pathtypes, '/', lpad(membername[idx+1:end], 8, '0')))
                fieldtypes[i] = jldatatype(parent, memberdtype)
            else
                memberclass = HDF5.h5t_get_member_class(dtype.id, i-1)
                if memberclass == HDF5.H5T_REFERENCE
                    # Field is a reference, so use Any
                    fieldtypes[i] = Any
                else
                    # Type is built-in
                    memberdtype = HDF5Datatype(HDF5.h5t_get_member_type(dtype.id, i-1), parent.plain)
                    fieldtypes[i] = jldatatype(parent, memberdtype)
                end
            end
        end

        if startswith(savedname, "(") || startswith(savedname, "Core.Tuple{")
            # We're reconstructing a tuple
            typetuple(fieldtypes)
        else
            # We're reconstructing some other type
            @eval begin
                struct $name
                    $([:($(fieldnames[i])::$(fieldtypes[i])) for i = 1:nfields]...)
                end
                $name
            end
        end
    end
end

# Get the index of a type in the types group. This could be cached, but
# it's already many times faster than calling H5Iget_name with a lot of
# data in the file, and it only needs to be called once per type.
# Revisit if this ever turns out to be a bottleneck.
function typeindex(parent::JldFile, addr::HDF5.Haddr)
    gtypes = parent.plain[pathtypes]
    i = 1
    for x in gtypes
        if HDF5.objinfo(x).addr == addr
            return i
        end
        i += 1
    end
end


text:
"""
    CSysModule

Module for management of coordinate systems.
"""
module CSysModule

using FinEtools.FTypesModule: FInt, FFlt, FCplxFlt, FFltVec, FIntVec, FFltMat, FIntMat, FMat, FVec, FDataDict
import LinearAlgebra: norm, cross

"""
Type for coordinate system transformations.
"""
struct CSys{F<:Function}
    isconstant::Bool
    isidentity::Bool
    updatebuffer!::F # function to update the coordinate system matrix
    csmat::Array{FFlt, 2} # the coordinate system matrix (buffer); see
end

"""
    CSys(sdim::FInt, mdim::FInt, computecsmat::F) where {F<:Function}

Construct ccoordinate system when the function to compute the
rotation matrix is given.

The function signature:
`update!(csmatout::FFltMat, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt)`
where
- `csmatout`= output matrix buffer,
- `XYZ`= location  in physical coordinates,
- `tangents`= tangent vector matrix, tangents to the parametric coordinate
  curves  in the element,
- `fe_label`= finite element label.
"""
function CSys(sdim::FInt, mdim::FInt, computecsmat::F) where {F<:Function}
    csmat = fill(zero(FFlt), sdim, mdim); # Allocate buffer, in preparation for the first call
    return CSys(false, false, computecsmat, csmat);
end

"""
    CSys(csmat::FFltMat)

Construct coordinate system when the rotation matrix is given.
"""
function CSys(csmat::FFltMat)
    function updatebuffer!(csmatout::FFltMat, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt)
        return csmatout # nothing to be done here, the matrix is already in the buffer
    end
    return CSys(true, false, updatebuffer!, deepcopy(csmat));# fill the buffer with the given matrix
end

"""
    CSys(dim::FInt)

Construct coordinate system when the rotation matrix is the identity.

`dim` = is the space dimension.
"""
function CSys(dim::FInt)
    return CSys([i==j ? one(FFlt) : zero(FFlt) for i=1:dim, j=1:dim]);
end

"""
    CSys(sdim::FInt, mdim::FInt)

Construct coordinate system for isotropic-material used with isoparametric
finite elements.

- `sdim` = number of space dimensions,
- `mdim` = number of manifold dimensions of the finite element in which the
  coordinate system  is being evaluated.

# Warning
If  the coordinate system matrix  should be identity, better use the constructor
for this specific situation, `CSys(dim::FInt)`. That will be much more efficient.
# See also
`gen_iso_csmat`
"""
function CSys(sdim::FInt, mdim::FInt)
    csmat = fill(zero(FFlt), sdim, mdim); # Allocate buffer, prepare for the first call
    function updatebuffer!(csmatout::FFltMat, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt)
      gen_iso_csmat!(csmatout, XYZ, tangents, fe_label)
      return  csmatout
    end
    return CSys(false, false, updatebuffer!, csmat);
end

"""
    function updatecsmat!(self::CSys, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt)

Update the coordinate system orientation matrix.

The  coordinate system matrix is updated based upon the location `XYZ` of the
evaluation point, and possibly on the Jacobian matrix `tangents` within the
element in which the coordinate system matrix is evaluated,  or perhaps on the
label `fe_label` of the finite element.

After this function returns, the coordinate system matrix can be retrieved
from the buffer as `self.csmat`.
"""
function updatecsmat!(self::CSys, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt)
    self.updatebuffer!(self.csmat, XYZ, tangents, fe_label)
    return self.csmat
end

"""
    gen_iso_csmat(XYZ::FFltMat, tangents::FFltMat, fe_label::FInt)

Compute the coordinate system  for an isotropic material using information
available  by looking at the coordinate curves of isoparametric finite
elements.

- `XYZ`= location  in physical coordinates,
- `tangents`= tangent vector matrix, tangents to the parametric coordinate
  curves  in the element,
- `fe_label`= finite element label.

The basic assumption here is that the material is isotropic, and
therefore the choice of the material directions does not really matter as
long as they correspond to the dimensionality of the element. For
instance a one-dimensional element (L2 as an example) may be embedded
in a three-dimensional space.

This function assumes that it is being called for an mdim-dimensional manifold
element, which is embedded in a sdim-dimensional Euclidean space. If `mdim ==
sdim`, the coordinate system matrix is the identity; otherwise the local
coordinate directions are aligned with the linear subspace defined by the
tangent vectors.

# Warning

This *cannot* be reliably used to produce consistent stresses because each
quadrature point gets a local coordinate system which depends on the
orientation of the element.
"""
function gen_iso_csmat!(csmatout::FFltMat,
    XYZ::FFltMat, tangents::FFltMat, fe_label::FInt)
    sdim, mdim = size(tangents);
    if sdim == mdim # finite element embedded in space of the same dimension
        copyto!(csmatout, [i==j ? one(FFlt) : zero(FFlt) for i=1:sdim, j=1:sdim]);
    else # lower-dimensional finite element embedded in space of higher dimension
        @assert 0 < mdim < 3
        e1 = tangents[:,1]/norm(tangents[:,1]);
        if mdim == 1 # curve-like finite element
            copyto!(csmatout, e1);
        elseif mdim == 2 # surface-like finite element
            n = cross(e1, vec(tangents[:,2]/norm(tangents[:,2])));
            e2 = cross(n, e1);
            e2 = e2/norm(e2);
            csmatout[:,1] = e1
            csmatout[:,2] = e2
        end
    end
    return csmatout
end

end

ext: jl


text:
using Distributed; addprocs(1)
@everywhere using Outsource

# spawn worker
wc = outsource() do c
    state = (0, 1)
    while isopen(c)
        n = take!(c)
        # compute next n Fibonacci numbers
        for _ in 1:n
            state = (state[2], state[1] + state[2])
        end
        put!(c, state[1])
    end
end

# issue stateful work and retrieve result
put!(wc, 10)
take!(wc) # = 55
put!(wc, 10)
take!(wc) # = 6765


text:
using DataFrames, JLD2

# Here we define the parameters we want to try (define the grid)
# We will try every possible combination of these parameters
depths = collect(2:2:20)
maxfeatures = collect(2:2:10)
numestimators = [20, 40, 80, 160]

"""
	This function is going to create a DataFrame,
	where each row corresponds to the parameters for one experiment
"""
function getparameters(depths, maxfeatures, numestimators)
	# let's code this up live!

end

# Let's get these parameters into the right format

# and now let's save them to a JLD file

# Let's make a folder to save the results (if necessary)


text:
##############################################################################################
module IEEE39bus
##############################################################################################
# imports and exports
using Random, Distributions
using LinearAlgebra
export dx_dt

##############################################################################################
##############################################################################################
# Type union declarations for multiple dispatch and type aliases

# vectors and ensemble members of sample
VecA = Union{Vector{Float64}, SubArray{Float64, 1}}

# dictionary for model parameters
ParamDict = Union{Dict{String, Array{Float64}}, Dict{String, Vector{Float64}}}

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

function dx_dt(x::VecA, t::Float64, dx_params::ParamDict)

    """Time derivative of the phase and fequency of the effective-network swing eq model

    Input x is a 2 n_g vector of the phase and fequency at each of the n_g generator buses.
    The input dx_params is a ParamDict of all system parameters to be passed to the
    integration scheme.  The system is currenty defined autonomously to be run as an SDE,
    noise perturbed steady state."""

    # unpack the system parameters effective network of
    # Nishikawa, T., & Motter, A. E. (2015). Comparative analysis of existing
    # models for power-grid synchronization.
    A = dx_params["A"]::Array{Float64}
    D = dx_params["D"]::Array{Float64}
    H = dx_params["H"]::Array{Float64}
    K = dx_params["K"]::Array{Float64}
    γ = dx_params["γ"]::Array{Float64}
    ω = dx_params["ω"]::Array{Float64}

    # convert the effective bus coupling and passive injection to contain the change
    # of variable terms
    K = ω[1] * K / 2.0
    A = ω[1] * A / 2.0

    # unpack the phase and frequency at the n_g buses, with all phases listed first, then all
    # fequencies in the order of the bus index
    n_g = convert(Int, length(x) / 2)
    δ_1 = @view x[1:n_g]
    δ_2 = @view x[n_g+1:end]

    # define the vector of the derivatives
    dx = zeros(2 * n_g)

    # derivative of the phase equals frequency
    dx[1:n_g] .= δ_2

    # compute the derivative of the inertia normalized frequencies
    # entry j is defined as
    # A_j * ω/2 - D_j /2 * δ_2 - Σ_{i!=j} K * ω/2 * sin(δ_j - δ_i - γ_ij)
    for j in 1:n_g
        for i in 1:n_g
            if j != i
                # K is symmetric, we loop over the columns for faster memory access
                # with the same variable j as in the row index of the derivative
                dx[n_g + j] += -K[i, j] * sin(δ_1[j] - δ_1[i] - γ[i, j])
            end
        end
        # finally apply the remaining terms
        dx[n_g + j] += A[j] - δ_2[j] * D[j] / 2.0
    end
    # to compute the derivative of the frequencies, we finally
    # divide back out by the inertia
    dx[n_g + 1 : end] = dx[n_g + 1: end] ./ H
    return dx
end

##############################################################################################
# end module

end

max_stars_repo_name: sukhreens/DataAssimilationBenchmarks

max_issues_repo_name: sukhreens/DataAssimilationBenchmarks

max_forks_repo_name: sukhreens/DataAssimilationBenchmarks


text:
    # General string interface

    # Supertype for all strings in this package
    abstract type AbstractStaticString <: AbstractString end

    # A subtype for strings that are backed by a pointer and a length alone
    abstract type AbstractPointerString <: AbstractStaticString end

    # Lightweight type for taking a view into an existing string
    # Like MallocString, but not necessarily null-terminated
    struct StringView <: AbstractPointerString
        pointer::Ptr{UInt8}
        length::Int
    end
    @inline Base.length(s::StringView) = s.length

    # Custom replshow for interactive use (n.b. _NOT_ static-compilerable)
    function Base.show(io::IO, s::StringView)
        Base.print(io, "$(length(s))-byte StringView:\n \"")
        Base.escape_string(io, Base.unsafe_string(pointer(s), length(s)))
        Base.print(io, "\"")
    end

    # Custom printing
    @inline Base.print(s::AbstractStaticString) = printf(s)
    @inline Base.println(s::AbstractStaticString) = puts(s)
    @inline Base.print(fp::Ptr{FILE}, s::AbstractStaticString) = printf(fp, s)
    @inline Base.println(fp::Ptr{FILE}, s::AbstractStaticString) = puts(fp, s)

    # Fundamentals
    @inline Base.unsafe_convert(::Type{Ptr{T}}, s::AbstractPointerString) where {T} = Ptr{T}(s.pointer)
    @inline Base.pointer(s::AbstractPointerString) = s.pointer
    @inline Base.sizeof(s::AbstractPointerString) = s.length
    @inline function Base.:(==)(a::AbstractStaticString, b::AbstractStaticString)
        (N = length(a)) == length(b) || return false
        pa, pb = pointer(a), pointer(b)
        for n in 0:N-1
            unsafe_load(pa + n) == unsafe_load(pb + n) || return false
        end
        return true
    end
    @inline function Base.:(==)(a::AbstractStaticString, b::AbstractString)
        GC.@preserve a b begin
            (N = length(a)) == sizeof(b) || return false
            pa, pb = pointer(a), pointer(b)
            for n in 0:N-1
                unsafe_load(pa + n) == unsafe_load(pb + n) || return false
            end
            return true
        end
    end
    @inline function Base.:(==)(a::AbstractString, b::AbstractStaticString)
        GC.@preserve a b begin
            (N = sizeof(a)) == length(b) || return false
            pa, pb = pointer(a), pointer(b)
            for n in 0:N-1
                unsafe_load(pa + n) == unsafe_load(pb + n) || return false
            end
            return true
        end
    end

    # Some of the AbstractArray interface:
    @inline Base.firstindex(s::AbstractStaticString) = 1
    @inline Base.lastindex(s::AbstractPointerString) = s.length
    @inline Base.eachindex(s::AbstractPointerString) = 1:s.length
    @inline Base.getindex(s::AbstractStaticString, i::Int) = unsafe_load(pointer(s)+(i-1))
    @inline Base.setindex!(s::AbstractStaticString, x::UInt8, i::Integer) = unsafe_store!(pointer(s)+(i-1), x)
    @inline Base.setindex!(s::AbstractStaticString, x, i::Integer) = unsafe_store!(pointer(s)+(i-1), convert(UInt8,x))
    @inline Base.getindex(s::AbstractStaticString, ::Colon) = s
    @inline function Base.getindex(s::AbstractStaticString, r::UnitRange{<:Integer})
        i₀ = max(first(r), 1) - 1
        l = min(length(r), length(s))
        return StringView(pointer(s)+i₀, l)
    end
    @inline function Base.setindex!(s::AbstractStaticString, x, r::UnitRange{<:Integer})
        is₀ = first(r)-1
        ix₀ = firstindex(x)-1
        @inbounds for i = 1:length(r)
            s[i+is₀] = x[i+ix₀]
        end
    end
    @inline function Base.setindex!(s::AbstractStaticString, x, ::Colon)
        ix₀ = firstindex(x)-firstindex(s)
        @inbounds for i ∈ eachindex(s)
            s[i] = x[i+ix₀]
        end
    end

    # Some of the AbstractString interface
    @inline Base.ncodeunits(s::AbstractPointerString) = s.length
    @inline Base.codeunits(s::AbstractPointerString) = MallocVector{UInt8}(s.pointer, s.length)
    @inline Base.codeunit(s::AbstractStaticString) = UInt8
    @inline Base.codeunit(s::AbstractStaticString, i::Integer) = s[i]


text:
#######################################################################################################################
# CALCULATE MARGINAL DAMAGES FROM A CH₄ EMISSION PULSE WITH FUND.
#######################################################################################################################
# Description: This function uses the FUND model to calculate marginal climate damages that occur from a CH₄ emissions
#              pulse based on the SNEASY+CH4 climate projections.
#
# Function Arguments:
#
#       base_temp  = An array of baseline SNEASY+CH4 temperature projections (each row = climate projection using a new
#                    posterior parameter sample, each column = a new year).
#       base_co2    = An array of baseline SNEASY+CH4 CO₂ concentration projections (each row = climate projection using
#                     a new posterior parameter sample, each column = a new year).
#       pulse_temp = An array of SNEASY+CH4 temperature projections similar to `base_temp`, but with a pulse of CH₄
#                    emissions in a user-specified year.
#       pulse_co2  = An array of SNEASY+CH4 CO₂ concentration projections similar to `base_co2`, but with a pulse of CH₄
#                    emissions in a user-specified year.
#       end_year   = The final year to run DICE.
#----------------------------------------------------------------------------------------------------------------------

function fund_damages(base_temp::Array{Float64,2}, base_co2::Array{Float64,2}, pulse_temp::Array{Float64,2}, pulse_co2::Array{Float64,2}, end_year::Int)

    #Create vector of years FUND is run (note FUND initializes in 1950).
    fund_years = collect(1950:end_year)

    #Find indices in SNEASY+CH4 climate projections that correspond to FUND years (assuming SNEASY+CH4 initialized in 1765).
    sneasy2fund_indices = findall((in)(fund_years), 1765:end_year)

    #Extract climate output that corresponds to FUND years (each row is a new model run, each column a new year).
    base_temp  = base_temp[:, sneasy2fund_indices]
    base_co2   = base_co2[:, sneasy2fund_indices]
    pulse_temp = pulse_temp[:, sneasy2fund_indices]
    pulse_co2  = pulse_co2[:, sneasy2fund_indices]

    #Caluclate number of times to run FUND (equal to number of posterior parameter samples used in climate projections).
    number_samples = size(base_temp, 1)

    #Get a base and CH₄ emissions pulse version of FUND model.
    fund_base  = create_iam_fund(end_year=end_year)
    fund_pulse = create_iam_fund(end_year=end_year)

    # Pre-allocate arrays to store relevant FUND output (Note: FUND has 16 geographic regions).
    marginal_damages             = zeros(length(fund_years), 16, number_samples)
    regional_population_base     = zeros(length(fund_years), 16, number_samples)
    regional_consumption_base    = zeros(length(fund_years), 16, number_samples)

    # Caluclate marginal damages and per capita consumption values for each posterior parameter climate projection.
    for i in 1:number_samples

        # Wrap code in a try/catch statement in case extreme climate values cause a model error.
        try

            # FUND cannot accept negative temperature values (which may occur in a few early periods due to superimposed noise). In both IAMs, set negative anomolies to 1e-15 to maintain consistency.
            if minimum(base_temp[i, :]) < 0.0
                base_temp[i, findall(x -> x <= 0.0, base_temp[i,:])] .= 1e-15
                pulse_temp[i, findall(x -> x <= 0.0, pulse_temp[i,:])] .= 1e-15
            end

            # Assign SNEASY+CH4 projections for temperature and CO₂ concentrations to relevant FUND components.
            update_param!(fund_base, :temp, base_temp[i, :])
            update_param!(fund_base, :acco2, base_co2[i, :])

            update_param!(fund_pulse, :temp, pulse_temp[i, :])
            update_param!(fund_pulse, :acco2, pulse_co2[i, :])

            #Run the two versions of FUND.
            run(fund_base)
            run(fund_pulse)

            #Calculate annual marginal climate damages (note: 1950 damage = `missing`, so set to 0.0 for convenience).
            marginal_damages[2:end,:,i] = fund_pulse[:impactaggregation, :loss][2:end,:] .- fund_base[:impactaggregation, :loss][2:end,:]

            # Calculate population for baseline scenario (note: FUND population units = millions of people).
            regional_population_base[:,:,i] = fund_base[:socioeconomic, :population] .* 1e6

            # Calculate regional consumption.
            regional_consumption_base[:,:,i] = fund_base[:socioeconomic, :consumption]

        # If climate projections cause a model error, set FUND output to -99999.99.
        catch
            marginal_damages[:,:,i] .= -99999.99
            regional_population_base[:,:,i] .= -99999.99
            regional_consumption_base[:,:,i] .= -99999.99
        end
    end

    # Check if any samples caused a model error.
    error_indices = findall(x-> x == -99999.99, marginal_damages[1,1,:])
    good_indices  = findall(x-> x != -99999.99, marginal_damages[1,1,:])

    # Return results.
    return marginal_damages[:,:,good_indices], regional_population_base[:,:,good_indices], regional_consumption_base[:,:,good_indices], error_indices, good_indices
end

#######################################################################################################################
# CALCULATE THE SC-CH4 USING FUND.
#######################################################################################################################
# Description: This function uses the marginal climate damage estimates from DICE to calculate the SC-CH4.
#
# Function Arguments:
#
#       marginal_damages     = Annual marginal climate damages (row = year, column = FUND region, 3rd dimension = model run
#                              using a new posterior climate parameter sample).
#       regional_consumption = Regional consumption levels.
#       regional_population  = Regional population levels.
#       pulse_year           = The year that the CH₄ emission pulse is added.
#       end_year             = The final year to run FUND.
#       constant             = Boolean parameter for whether or not to use a constant discount rate (true = use constant discount rate).
#       η                    = Parameter representing the degree of intertemporal inequality aversion.
#       γ                    = Parameter representing the degree of regional inequality aversion.
#       ρ                    = Pure rate of time preference.
#       dollar_conversion    = Inflation conversion factor to scale dollar units to a different year.
#       equity_weighting     = Should the model produce equity-weighted SC-CH4 estimates? (true = produce equity-weighted estimates).
#----------------------------------------------------------------------------------------------------------------------

function fund_scch4(marginal_damages::Array{Float64,3}, regional_consumption::Array{Float64,3}, regional_population::Array{Float64,3}, pulse_year::Int, end_year::Int; constant::Bool=true, η::Float64=1.5, γ::Float64=1.5, ρ::Float64=0.01, dollar_conversion::Float64=1.35, equity_weighting=false)

    #Create vector of years FUND is run (note FUND initializes in 1950).
    fund_years = collect(1950:end_year)

    # Find timestep index when CH₄ emission pulse occurs.
    pulse_index = findall(x -> x == pulse_year, fund_years)[1]

    # Calculate number of SC-CH4 values to calculate.
    number_samples = size(marginal_damages, 3)

    # Calculate global and regional per capita consumption values (for discounting & equity-weighting).
    pc_consumption_regional = regional_consumption ./ regional_population
    pc_consumption_global   = sum(regional_consumption, dims=2) ./ sum(regional_population, dims=2)

    # Allocate an array to hold either constant or Ramsey discount factors.
    df = zeros(length(fund_years))

    # If using a constant discount rate, calculate constant discount factors for annual timesteps.
    if constant == true
        calculate_discount_factors!(df, ones(length(fund_years)), ρ, 0.0, fund_years, pulse_year, pulse_index)
    end

    #-----------------------------------------------------------------------
    # Loop through each marginal damage estimate and calculate the SC-CH4.
    #-----------------------------------------------------------------------

    # Calculate the SC-CH4 for each posterior parameter sample without equity-weighting.
    if equity_weighting == false

        # Pre-allocate array to store results and SC-CH4 estimates.
        discounted_global_marginal_damages = zeros(number_samples, length(fund_years))
        scch4 = zeros(number_samples)

        # Loop through damages resulting from each posterior parameter sample and calculate the SC-CH4.
        for i = 1:number_samples
            # If not using a constant discount rate, calculate a Ramsey discount rate for each posterior parameter sample based on global per capita consumption levels.
            if constant == false
                calculate_discount_factors!(df, pc_consumption_global[:,1,i], ρ, η, fund_years, pulse_year, pulse_index)
            end

            # Calculate discounted global annual marginal damages (summed across the FUND regions).
            discounted_global_marginal_damages[i,:] = @views sum(marginal_damages[:,:,i] .* df, dims=2)

            # Calculate the SC-CH4 as the sum of discounted annual marginal damages (scaled by inflation conversion factor for desired year).
            scch4[i] = @views sum(discounted_global_marginal_damages[i,:]) * dollar_conversion
        end

        # Return SC-CH4 and discounted annual marginal damage values.
        return scch4, discounted_global_marginal_damages

    # Calculate the SC-CH4 for each posterior parameter sample using equity-weighting.
    elseif equity_weighting == true

        # Pre-allocate array to store equity-weighted SC-CH4 estimates for each region as well as some temporary variables (for convenience).
        scch4           = zeros(number_samples, 16)
        tempvar_1       = zeros(16)
        annual_vals     = zeros(length(fund_years))
        cde_consumption = zeros(length(fund_years), number_samples)

        # Following the approach of Anthoff & Emmerling (2019), calculate equity-weighted SC-CH4 values for each normalization region.
        for i = 1:number_samples
            for r = 1:16
                for t = pulse_index:length(fund_years)
                    cde_consumption[t,i] = sum(regional_population[t,:,i] ./ sum(regional_population[t,:,i]) .* pc_consumption_regional[t,:,i] .^(1-γ)) ^ (1/(1-γ))
                    tempvar_1[:] .= @views (pc_consumption_regional[t,:,i] ./ pc_consumption_regional[pulse_index,r,i]) .^ -γ
                    tempvar_2    = (cde_consumption[t,i] ./ cde_consumption[pulse_index,i]) .^ (γ-η)
                    annual_vals[t] = @views sum(tempvar_1 .* tempvar_2 .* (1+ρ) ^ (-t+pulse_index) .* marginal_damages[t,:,i])
                end
                scch4[i,r] = sum(annual_vals) * dollar_conversion
            end
        end

        # Calculate indices for equity weighted estimates that produce an error (implausible parameter combinations can produce -Inf values).
        error_indices = unique(getindex.(findall(x -> !isfinite(x), scch4), 1))
        good_indices  = findall(!in(error_indices), collect(1:number_samples))

        # Return equity-weighted SC-CH4 values.
        return scch4[good_indices, :]
    end
end


text:
using BioSymbols
using BioSequences
using Base.Test

using Rifraf

include("test_utils.jl")

import Rifraf.forward,
       Rifraf.forward_moves,
       Rifraf.backward

# TODO: some codon tests
# TODO: test all update() variants
@testset "test forward and backward" begin
    const scores = Scores(-1.0, -1.0, -1.0, -Inf, -Inf)

    @testset "perfect forward" begin
        bandwidth = 1
        template = DNASeq("AA")
        seq = DNASeq("AA")
        lp = -3.0
        match = inv_log10(lp)
        log_p = fill(lp, length(seq))
        pseq = RifrafSequence(seq, log_p, bandwidth, scores)
        A = forward(template, pseq)
        # transpose because of column-major order
        expected = transpose(reshape([[0.0, lp + scores.deletion, 0.0];
                                      [lp + scores.insertion, match, match + lp + scores.deletion];
                                      [0.0, match + lp + scores.insertion, 2 * match]],
                                     (3, 3)))
        @test full(A) ≈ expected

        A2, _ = forward_moves(template, pseq)
        @test full(A2) ≈ full(A)
    end

    @testset "perfect backward" begin
        bandwidth = 1
        template = DNASeq("AA")
        seq = DNASeq("AA")
        lp = -3.0
        match = inv_log10(lp)
        log_p = fill(lp, length(seq))
        pseq = RifrafSequence(seq, log_p, bandwidth, scores)
        A = backward(template, pseq)
        # transpose because of column-major order
        expected = transpose(reshape([[2 * match, match + lp + scores.insertion, 0.0];
                                      [match + lp + scores.deletion, match, lp + scores.insertion];
                                      [0.0, lp + scores.deletion, 0.0]],
                                     (3, 3)))
        @test full(A) ≈ expected
    end

    @testset "imperfect forward" begin
        bandwidth = 1
        template = DNASeq("AA")
        seq = DNASeq("AT")
        lp = -3.0
        match = inv_log10(lp)
        log_p = fill(lp, length(seq))
        pseq = RifrafSequence(seq, log_p, bandwidth, scores)
        A1 = forward(template, pseq)
        B = backward(template, pseq)
        check_all_cols(A1, B, false)
        expected = transpose(reshape([[0.0, lp + scores.deletion, 0.0];
                                      [lp + scores.insertion, match, match + lp + scores.deletion];
                                      [0.0, match + lp + scores.insertion, match + lp + scores.mismatch]],
                                     (3, 3)))

        @test full(A1) ≈ expected atol=0.01

        A2, moves = forward_moves(template, pseq)
        @test full(A1) ≈ full(A2) atol=0.01
    end

    @testset "imperfect backward" begin
        bandwidth = 1
        template = DNASeq("AA")
        seq = DNASeq("AT")
        lp = -3.0
        match = inv_log10(lp)
        log_p = fill(lp, length(seq))
        pseq = RifrafSequence(seq, log_p, bandwidth, scores)
        B = backward(template, pseq)
        expected = transpose(reshape([[lp + scores.mismatch + match, lp + scores.insertion + match, 0.0];
                                      [2*lp + scores.deletion + scores.mismatch, lp + scores.mismatch, lp + scores.insertion];
                                      [0.0, lp + scores.deletion, 0.0]],
                                     (3, 3)))
        @test full(B) ≈ expected atol=0.01
    end

    @testset "forward/backward agreement" begin
        @testset "forward/backward agreement 1" begin
            template = DNASeq("TG")
            seq = DNASeq("GTCG")
            log_p = [-1.2, -0.8, -0.7, -1.0]
            bandwidth = 5
            local_scores = Scores(ErrorModel(2.0, 1.0, 1.0, 3.0, 3.0))
            pseq = RifrafSequence(seq, log_p, bandwidth, local_scores)
            A = forward(template, pseq)
            B = backward(template, pseq)
            check_all_cols(A, B, true)

            A2, _ = forward_moves(template, pseq)
            @test full(A) ≈ full(A2) atol=0.01
        end

        @testset "forward/backward agreement 2" begin
            template = DNASeq("GCACGGTC")
            seq = DNASeq("GACAC")
            log_p = [-1.1, -1.1, -0.4, -1.0, -0.7]
            bandwidth = 5
            local_scores = Scores(ErrorModel(2.0, 1.0, 1.0, 3.0, 3.0))
            pseq = RifrafSequence(seq, log_p, bandwidth, local_scores)
            A = forward(template, pseq)
            B = backward(template, pseq)
            check_all_cols(A, B, true)

            A2, _ = forward_moves(template, pseq)
            @test full(A) ≈ full(A2) atol=0.01
        end
    end

    @testset "insertion_agreement" begin
        template = DNASeq("AA")
        seq = DNASeq("ATA")
        bandwidth = 10
        log_p = [-5.0, -1.0, -6.0]
        pseq = RifrafSequence(seq, log_p, bandwidth, scores)
        A = forward(template, pseq)
        B = backward(template, pseq)
        score = (inv_log10(log_p[1]) +
                 log_p[2] + scores.insertion +
                 inv_log10(log_p[3]))
        @test A[end, end] ≈ score
        check_all_cols(A, B, false)

        A2, _ = forward_moves(template, pseq)
        @test full(A) == full(A2)
    end

    @testset "deletion agreement 1" begin
        template = DNASeq("GATAG")
        seq = DNASeq("GAAG")
        bandwidth = 10
        log_p = [-5.0, -2.0, -1.0, -6.0]
        pseq = RifrafSequence(seq, log_p, bandwidth, scores)
        A = forward(template, pseq)
        B = backward(template, pseq)
        score = (pseq.match_scores[1] +
                 pseq.match_scores[2] +
                 pseq.del_scores[3] +
                 pseq.match_scores[3] +
                 pseq.match_scores[4])
        @test A[end, end] ≈ score
        check_all_cols(A, B, false)

        A2, _ = forward_moves(template, pseq)
        @test full(A) == full(A2)
    end

    @testset "deletion agreement 2" begin
        template = DNASeq("ATA")
        seq = DNASeq("AA")
        bandwidth = 10
        log_p = [-2.0, -3.0]
        pseq = RifrafSequence(seq, log_p, bandwidth, scores)
        A = forward(template, pseq)
        B = backward(template, pseq)
        score = (pseq.match_scores[1] +
                 pseq.del_scores[2] +
                 pseq.match_scores[2])
        @test A[end, end] ≈ score
        check_all_cols(A, B, false)

        A2, _ = forward_moves(template, pseq)
        @test full(A) == full(A2)
    end
end

@testset "alignment" begin
    const errors = Rifraf.normalize(ErrorModel(1.0, 1.0, 1.0, 0.0, 0.0))
    const scores = Scores(errors)

    @testset "align 1" begin
        template = DNASeq("ATAA")
        seq = DNASeq("AAA")
        bandwidth = 10
        log_p = [-2.0, -3.0, -3.0]
        pseq = RifrafSequence(seq, log_p, bandwidth, scores)
        moves = Rifraf.align_moves(template, pseq)
        t, s = Rifraf.moves_to_aligned_seqs(moves, template, seq)
        @test t == DNASequence("ATAA")
        @test s == DNASequence("A-AA")
    end

    @testset "align 2" begin
        template = DNASeq("AACCTT")
        seq = DNASeq("AAACCCTT")
        bandwidth = 10
        log_p = fill(log10(0.1), length(seq))
        pseq = RifrafSequence(seq, log_p, bandwidth, scores)
        moves = Rifraf.align_moves(template, pseq)
        t, s = Rifraf.moves_to_aligned_seqs(moves, template, seq)
        @test t[end-1:end] == DNASequence("TT")
    end

    @testset "moves_to_indices 1" begin
        template = DNASeq("AAA")
        seq = DNASeq("AAA")
        bandwidth = 10
        log_p = fill(log10(0.1), length(seq))
        pseq = RifrafSequence(seq, log_p, bandwidth, scores)
        moves = Rifraf.align_moves(template, pseq)
        indices = Rifraf.moves_to_indices(moves, length(template), length(pseq))
        @test indices == collect(1:3)
    end

    @testset "moves_to_indices 2" begin
        template = DNASeq("AAA")
        seq = DNASeq("AAAT")
        bandwidth = 10
        log_p = fill(log10(0.1), length(seq))
        pseq = RifrafSequence(seq, log_p, bandwidth, scores)
        moves = Rifraf.align_moves(template, pseq)
        indices = Rifraf.moves_to_indices(moves, length(template), length(pseq))
        @test indices == collect(1:3)
    end

    @testset "moves_to_indices 3" begin
        template = DNASeq("AAAT")
        seq = DNASeq("AAA")
        bandwidth = 10
        log_p = fill(log10(0.1), length(seq))
        pseq = RifrafSequence(seq, log_p, bandwidth, scores)
        moves = Rifraf.align_moves(template, pseq)
        indices = Rifraf.moves_to_indices(moves, length(template), length(pseq))
        @test indices == [1, 2, 3, 3]
    end

    @testset "moves_to_indices 4" begin
        template = DNASeq("TAAA")
        seq = DNASeq("AAA")
        bandwidth = 10
        log_p = fill(log10(0.1), length(seq))
        pseq = RifrafSequence(seq, log_p, bandwidth, scores)
        moves = Rifraf.align_moves(template, pseq)
        indices = Rifraf.moves_to_indices(moves, length(template), length(pseq))
        @test indices == [0, 1, 2, 3]
    end

    @testset "align and skew" begin
        ref_errors = ErrorModel(10.0, 1e-10, 1e-10, 1.0, 1.0)
        ref_scores = Scores(ref_errors)

        ref = DNASeq("CGGCGATTT")
        consensus_errors = LogProb[-8., -8., -8., -1., -8., -10., -10.]
        consensus = RifrafSequence(DNASeq("CTGCCGA"), consensus_errors, 10, ref_scores)
        a, b = Rifraf.align(ref, consensus, skew_matches=true)
        expected_a = dna"CGG-CGATTT"
        expected_b = dna"CTGCCGA---"
        @test a == expected_a
        @test b == expected_b
    end
end

@testset "align with self" begin
    # make sequence with long insertion
    seqstring = DNASequence("AAAGGGTTTCCC")
    seq = DNASeq(seqstring)
    errors = fill(0.1, length(seq))
    errors[1:6] = 0.3
    errors[end-3:end] = 0.45
    err_log_p = log10.(errors)
    bandwidth = 3
    scores = Scores(ErrorModel(1.0, 10.0, 10.0, 0.0, 0.0))
    rseq = RifrafSequence(seq, err_log_p, bandwidth, scores)

    a, b = Rifraf.align(seq, rseq)
    @test a == b
    @test a == seqstring
end


text:
using Test, LinearAlgebra

@testset "test/integration/local_polytope_facets.jl" begin

using BellScenario: LocalPolytope, ConvexPolytope, Degeneracy, Behavior

dir = "./test/integration/files"

@testset "21->12 polytope" begin
    α = (2,1)
    β = (1,2)

    vertices = LocalPolytope.vertices(α, β, bits=1)

    # There should be (1*2)(1 + 1*(2 - 1)) vertices for fixed direction
    @test length(vertices) == 4
    @test size(vertices[1]) == (2,1)

    constraints = ConvexPolytope.facet_constraints(vertices, dir=dir)

    @test length(constraints[1]) == 4
    @test constraints[1] == [[0 -1 0],[0 0 -1],[-1 0 1],[-1 1 0]]
    @test constraints[2] == []
    @test constraints[3] == []

    canonical_groups = Degeneracy.canonical_facets(α,β,constraints[1],"normalized")

    @test issetequal(constraints[1], collect(Iterators.flatten(canonical_groups)))
    @test length(canonical_groups) == 1
    @test issetequal(canonical_groups[1], [
        [0 -1 0],[0 0 -1],[-1 0 1],[-1 1 0]
    ])
end

@testset "31->12 polytope" begin
    α = (3,1)
    β = (1,2)

    vertices = LocalPolytope.vertices(α, β, bits=1)

    # There should be (1*2)(1 + 1*(2 - 1)) vertices for fixed direction
    @test length(vertices) == 8
    @test size(vertices[1]) == (3,1)

    constraints = ConvexPolytope.facet_constraints(vertices, dir=dir)

    @test length(constraints[1]) == 6
    @test constraints[2] == []
    @test constraints[3] == []

    canonical_groups = Degeneracy.canonical_facets(α,β,constraints[1],"normalized")

    @test issetequal(constraints[1], collect(Iterators.flatten(canonical_groups)))
    @test length(canonical_groups) == 1
    @test issetequal(canonical_groups[1], [
        [0 -1 0 0],[0 0 -1 0],[0 0 0 -1],[-1 0 1 0],[-1 1 0 0],[-1 0 0 1]
    ])
end

@testset "31->13 polytope" begin
    α = (3,1)
    β = (1,3)

    vertices = LocalPolytope.vertices(α, β, bits=1)

    # There should be (1*3)(1 + 3*(3 - 1)) vertices for fixed direction
    @test length(vertices) == 21
    @test size(vertices[1]) == (6,1)

    constraints = ConvexPolytope.facet_constraints(vertices, dir=dir)

    @test length(constraints[1]) == 15
    @test constraints[2] == []
    @test constraints[3] == []

    canonical_groups = Degeneracy.canonical_facets(α,β,constraints[1],"normalized")

    @test issetequal(constraints[1], collect(Iterators.flatten(canonical_groups)))
    @test length(canonical_groups) == 2
    @test issetequal(canonical_groups[1], [
        [0 -1 0 0 0 0 0], [0 0 -1 0 0 0 0], [0 0 0 -1 0 0 0],
        [0 0 0 0 -1 0 0], [0 0 0 0 0 -1 0], [0 0 0 0 0 0 -1],
        [-1 0 0 1 0 0 1], [-1 0 1 0 0 1 0], [-1 1 0 0 1 0 0]
    ])
    @test issetequal(canonical_groups[2], [
        [-1 -1 0 1 -1 1 0], [-1 -1 1 0 -1 0 1], [-1 0 -1 1 1 -1 0],
        [-1 0 1 -1 1 0 -1], [-1 1 -1 0 0 -1 1], [-1 1 0 -1 0 1 -1]
    ])
end

@testset "41->14 polytope" begin
    α = (4,1)
    β = (1,4)

    vertices = LocalPolytope.vertices(α, β, bits=1)

    # There should be (1*4)(1 + 7*(4 - 1)) vertices for fixed direction
    @test length(vertices) == 88
    @test size(vertices[1]) == (12,1)

    constraints = ConvexPolytope.facet_constraints(vertices, dir=dir)

    @test length(constraints[1]) == 664
    @test constraints[2] == []
    @test constraints[3] == []

    canonical_groups = Degeneracy.canonical_facets(α,β,constraints[1],"normalized")

    ggame_facet = [-2 1 0 0 0 0 1 0 0 0 0 1 0 0 0 0 1]*Behavior.norm_to_gen_proj(α,β)
    @test in(ggame_facet, canonical_groups[2])

    @test issetequal(constraints[1], collect(Iterators.flatten(canonical_groups)))
    @test length(canonical_groups) == 6
end

@testset "22->12 polytope" begin
    α = (2,2)
    β = (1,2)

    vertices = LocalPolytope.vertices(α, β, bits=1)

    # There should be (4*2)(1 + 1*(2 - 1)) vertices for fixed direction
    @test length(vertices) == 16
    @test size(vertices[1]) == (6,1)

    constraints = ConvexPolytope.facet_constraints(vertices, dir=dir)

    @test length(constraints[1]) == 8
    @test constraints[2] == []
    @test constraints[3] == []

    canonical_groups = Degeneracy.canonical_facets(α,β,constraints[1],"fixed-direction")

    @test issetequal(constraints[1], collect(Iterators.flatten(canonical_groups)))
    @test length(canonical_groups) == 1
end

@testset "22->22 polytope (fixed-direction)" begin
    α = (2,2)
    β = (2,2)

    vertices = LocalPolytope.vertices(α, β, bits=1, rep="fixed-direction")

    # There should be (4*4)(1 + 1*(4 - 1)) vertices for fixed direction
    @test length(vertices) == 64
    @test size(vertices[1]) == (10,1)

    constraints = ConvexPolytope.facet_constraints(vertices, dir=dir)
    @test constraints[2] == []
    @test constraints[3] == []

    @test length(constraints[1]) == 16

    canonical_groups = Degeneracy.canonical_facets(α,β,constraints[1],"fixed-direction")

    @test issetequal(constraints[1], collect(Iterators.flatten(canonical_groups)))
    @test length(canonical_groups) == 1
end

@testset "33->13 polytope" begin
    α = (3,3)
    β = (1,3)

    vertices = LocalPolytope.vertices(α, β, bits=1)

    # There should be (27*3)(1 + 3*(3 - 1)) vertices for fixed direction
    @test length(vertices) == 567
    @test size(vertices[1]) == (24,1)

    constraints = ConvexPolytope.facet_constraints(vertices, dir=dir)

    @test length(constraints[1]) == 33
    @test constraints[2] == []
    @test constraints[3] == []

    canonical_groups = Degeneracy.canonical_facets(α,β,constraints[1],"fixed-direction")

    @test issetequal(constraints[1], collect(Iterators.flatten(canonical_groups)))
    @test length(canonical_groups) == 2
end

@testset "32->22 polytope" begin
    α = (3,2)
    β = (2,2)

    vertices = LocalPolytope.vertices(α, β, bits=1)

    # There should be (27*3)(1 + 3*(3 - 1)) vertices for fixed direction
    @test length(vertices) == 320
    @test size(vertices[1]) == (15,1)

    constraints = ConvexPolytope.facet_constraints(vertices, dir=dir)

    @test length(constraints[1]) == 864
    @test constraints[2] == []
    @test constraints[3] == []

    canonical_groups = Degeneracy.canonical_facets(α,β,constraints[1],"fixed-direction")

    @test issetequal(constraints[1], collect(Iterators.flatten(canonical_groups)))
    @test length(canonical_groups) == 9
end

@testset "32->12 polytope" begin
    α = (3,2)
    β = (1,2)

    vertices = LocalPolytope.vertices(α, β, bits=1)

    @test length(vertices) == 64
    @test size(vertices[1]) == (9,1)

    constraints = ConvexPolytope.facet_constraints(vertices, dir=dir)

    @test length(constraints[1]) == 12
    @test constraints[2] == []
    @test constraints[3] == []

    canonical_groups = Degeneracy.canonical_facets(α,β,constraints[1],"fixed-direction")

    @test issetequal(constraints[1], collect(Iterators.flatten(canonical_groups)))
    @test length(canonical_groups) == 1
end

@testset "41->12 polytope" begin
    α = (4,1)
    β = (1,2)

    vertices = LocalPolytope.vertices(α, β, bits=1)

    @test length(vertices) == 16
    @test size(vertices[1]) == (4,1)

    constraints = ConvexPolytope.facet_constraints(vertices, dir=dir)

    @test length(constraints[1]) == 8
    @test constraints[2] == []
    @test constraints[3] == []

    canonical_groups = Degeneracy.canonical_facets(α,β,constraints[1],"normalized")

    @test issetequal(constraints[1], collect(Iterators.flatten(canonical_groups)))
    @test length(canonical_groups) == 1
end

@testset "33->12 polytope" begin
    α = (3,3)
    β = (1,2)

    vertices = LocalPolytope.vertices(α, β, bits=1)

    @test length(vertices) == 216
    @test size(vertices[1]) == (15,1)

    constraints = ConvexPolytope.facet_constraints(vertices, dir=dir)

    @test length(constraints[1]) == 18
    @test constraints[2] == []
    @test constraints[3] == []

    canonical_groups = Degeneracy.canonical_facets(α,β,constraints[1],"fixed-direction")

    @test issetequal(constraints[1], collect(Iterators.flatten(canonical_groups)))
    @test length(canonical_groups) == 1
end

@testset "32->13 polytope" begin
    α = (3,2)
    β = (1,3)

    vertices = LocalPolytope.vertices(α, β, bits=1)

    @test length(vertices) == 168
    @test size(vertices[1]) == (15,1)

    constraints = ConvexPolytope.facet_constraints(vertices, dir=dir)

    @test length(constraints[1]) == 24
    @test constraints[2] == []
    @test constraints[3] == []

    canonical_groups = Degeneracy.canonical_facets(α,β,constraints[1],"fixed-direction")

    @test issetequal(constraints[1], collect(Iterators.flatten(canonical_groups)))
    @test length(canonical_groups) == 2
end

@testset "41->13 polytope" begin
    α = (4,1)
    β = (1,3)

    vertices = LocalPolytope.vertices(α, β, bits=1)

    @test length(vertices) == 45
    @test size(vertices[1]) == (8,1)

    constraints = ConvexPolytope.facet_constraints(vertices, dir=dir)

    @test length(constraints[1]) == 36
    @test constraints[2] == []
    @test constraints[3] == []

    canonical_groups = Degeneracy.canonical_facets(α,β,constraints[1],"normalized")

    @test issetequal(constraints[1], collect(Iterators.flatten(canonical_groups)))
    @test length(canonical_groups) == 2
end

@testset "31->14 polytope" begin
    α = (3,1)
    β = (1,4)

    vertices = LocalPolytope.vertices(α, β, bits=1)

    @test length(vertices) == 40
    @test size(vertices[1]) == (9,1)

    constraints = ConvexPolytope.facet_constraints(vertices, dir=dir)

    @test length(constraints[1]) == 72
    @test constraints[2] == []
    @test constraints[3] == []

    canonical_groups = Degeneracy.canonical_facets(α,β,constraints[1],"normalized")

    @test issetequal(constraints[1], collect(Iterators.flatten(canonical_groups)))
    @test length(canonical_groups) == 3
end

end


text:
using Multipactor
using Test

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


text:
"""

    fibonacci(iterations, [num1, num2])

Builds a series of fibonacci numbers `n` iterations in a vector

    Args: iterations: length of fibonacci sequence to be calculated

    Optional: num1 and num2 default to 0 and 1 respectively, the starting numbers of the sequence.

# Examples
```julia-repl

# one argument -  starting values are 0 and 1
julia> fibonacci(10)
12-element Array{Int64,1}:
  0
  1
  1
  2
  3
  5
  8
 13
 21
 34
 55
 89

# two arguments - starting values are 4 and 1
julia> fibonacci(10, 4)
12-element Array{BigInt,1}:
   4
   1
   5
   6
  11
  17
  28
  45
  73
 118
 191
 309

# Three arguments - starting values are 3 and -1
julia> fibonacci(10, 3, -1)
12-element Array{Int64,1}:
  3
 -1
  2
  1
  3
  4
  7
 11
 18
 29
 47
 76
 ```
 """
function fibonacci(iterations, num1=0, num2=1)
    arr = (BigInt)[num1, num2]
    for i in 1:iterations
        push!(arr, +(arr[end], arr[end - 1]))
    end
    return arr
end

max_stars_repo_name: jawuku/julia

max_issues_repo_name: jawuku/julia

max_forks_repo_name: jawuku/julia


text:
module Abstract

"""
    UnionArray{T, N}
    UnionVector{T}
    UnionMatrix{T}

`UnionArray` stores heterogeneous elements without indirections (which
would happen in `Array{Any}`).

# Examples
```jldoctest
julia> using UnionArrays

julia> xs = UnionVector(Any[UInt8(1), 2.0, (a=1, b=2)]);

julia> xs[1]
0x01

julia> xs[2]
2.0

julia> xs[3]
(a = 1, b = 2)

julia> xs[1] = (a=3, b=4);

julia> xs[1]
(a = 3, b = 4)

julia> M = reshape(xs, (1, :)) :: UnionMatrix;

julia> M[1, 1]
(a = 3, b = 4)
```
"""
abstract type UnionArray{T, N} <: AbstractArray{T, N} end

const UnionVector{T} = UnionArray{T, 1}
const UnionMatrix{T} = UnionArray{T, 2}

end  # module


text:
using TensorNetworkAD
using Test
using Zygote, OMEinsum
using TensorNetworkAD: model_tensor, tensorfromclassical, trg_svd

@testset "unit" begin
    t = randn(10,10,10,10)
    u, v = trg_svd(t, 100, 0)
    @test t ≈ ein"ija,akl -> ijkl"(u,v)
end
@testset "real" begin
    χ, niter = 5, 5
    foo = β -> trg(model_tensor(Ising(),β), χ, niter)
    # the pytorch result with tensorgrad
    # https://github.com/wangleiphy/tensorgrad
    # clone this repo and type
    # $ python 1_ising_TRG/ising.py -chi 5 -Niter 5
    @test foo(0.4) ≈ 0.8919788686747141
    @test num_grad(foo, 0.4, δ=1e-6) ≈ Zygote.gradient(foo, 0.4)[1]
end

@testset "complex" begin
    β, χ, niter = 0.4, 12, 3
    @test trg(model_tensor(Ising(),β), χ, niter) ≈
        real(trg(model_tensor(Ising(),β) .+ 0im, χ, niter))
    @test Zygote.gradient(β -> trg(model_tensor(Ising(),β), χ, niter), 0.4)[1] ≈
        real(Zygote.gradient(β -> real(trg(model_tensor(Ising(),β) .+ 0im, χ, niter)), 0.4)[1])
end


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    Hexahedron(p1, p2, ..., p8)

A hexahedron with points `p1`, `p2`, ..., `p8`.
"""
struct Hexahedron{Dim,T,V<:AbstractVector{Point{Dim,T}}} <: Polyhedron{Dim,T}
  vertices::V
end

nvertices(::Type{<:Hexahedron}) = 8
nvertices(h::Hexahedron) = nvertices(typeof(h))

function (h::Hexahedron)(u, v, w)
  A1, A2, A4, A3,
  A5, A6, A8, A7 = coordinates.(h.vertices)
  Point((1-u)*(1-v)*(1-w)*A1 +
            u*(1-v)*(1-w)*A2 +
            (1-u)*v*(1-w)*A3 +
                u*v*(1-w)*A4 +
            (1-u)*(1-v)*w*A5 +
                u*(1-v)*w*A6 +
                (1-u)*v*w*A7 +
                    u*v*w*A8)
end


text:
using WeavePynb

using CwJWeaveTpl

fnames = [
          "derivatives", ## more questions
          "numeric_derivatives",

          "mean_value_theorem",
          "optimization",
          "curve_sketching",

          "linearization",
          "newtons_method",
          "lhopitals_rule",  ## Okay - -but could beef up questions..

          "implicit_differentiation", ## add more questions?
          "related_rates",
          "taylor_series_polynomials"
]

process_file(nm; cache=:off) = CwJWeaveTpl.mmd(nm * ".jmd", cache=cache)

function process_files(;cache=:user)
    for f in fnames
        @show f
        process_file(f, cache=cache)
    end
end

"""
## TODO derivatives

tangent lines intersect at avearge for a parabola

Should we have derivative results: inverse functions, logarithmic differentiation...
"""


text: @show filter(iseven, 1:10)

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#   Tests related to the transformations of local reference frames.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

# File: ./src/transformations/local_frame.jl
# ==========================================

# Functions: ecef_to_ned and ned_to_ecef
# --------------------------------------

@testset "Functions ecef_to_ned and ned_to_ecef" begin
    r_ecef = SVector(R0, 0, -R0)
    lat = 0.5
    lon = 0.44
    h = 130e3

    r_ned = ecef_to_ned(r_ecef, lat, lon, h; translate = true)

    @test r_ned[1] ≈  -8.345950969454717e6
    @test r_ned[2] ≈  -2.7167002618976594e6
    @test r_ned[3] ≈   4.496865561149651e6

    r_ecef_conv = ned_to_ecef(r_ned, lat, lon, h; translate = true)

    @test r_ecef_conv[1] ≈ r_ecef[1]
    @test r_ecef_conv[2] ≈ r_ecef[2] atol = 1e-9 # r_ecef[2] is 0.
    @test r_ecef_conv[3] ≈ r_ecef[3]

    r_ned = ecef_to_ned(r_ecef, lat, lon, h)
    @test norm(r_ned) == norm(r_ecef)

    B_l = SVector(24178.985570422887, -39579.98354881559, -14897.692106044478)

    # This is a validated code that rotates NED to ECEF without translations.
    D_ecef_l = angle_to_dcm(lat, -lon, 0, :YZX)
    B_ecef_exp = D_ecef_l*SVector(-B_l[3], +B_l[2], +B_l[1])
    B_ecef = ned_to_ecef(B_l, lat, lon, h)

    @test B_ecef[1] ≈ B_ecef_exp[1]
    @test B_ecef[2] ≈ B_ecef_exp[2]
    @test B_ecef[3] ≈ B_ecef_exp[3]
end


text:
#include("treesMain.jl")
include("../utilities.jl")

function testIntervalTrees()
    #test add nodes
    rootNode = IntervalTreeNode(5,10,1)
    addNodeToIntervalTree!(rootNode, IntervalTreeNode(15,25))
    addNodeToIntervalTree!(rootNode, IntervalTreeNode(1,12,2))
    addNodeToIntervalTree!(rootNode, IntervalTreeNode(8,16))
    addNodeToIntervalTree!(rootNode, IntervalTreeNode(14,20))
    addNodeToIntervalTree!(rootNode, IntervalTreeNode(18,21))
    addNodeToIntervalTree!(rootNode, IntervalTreeNode(2,8))
    addNodeToIntervalTree!(rootNode, IntervalTreeNode(2.1,8))

    print("Ordered Intervals From Tree \n")
    printTree(rootNode)
    print("\n\n")

    #testTreeCopy, display at end
    treeCopy = deepCopyIntervalTreeNodes(rootNode, false)
    treeCopyMerged = deepCopyIntervalTreeNodes(rootNode, true)

    #test query interval
    print("Query Interval (8,10,1) \n")
    query1 = Interval(8,10)
    results = IntervalTreeNode[]
    findIntersectingIntervals!(query1, rootNode, results)
    for result in results
        print(result,"\n")
    end

    #test query interval again
    query2 = Interval(20,22)
    results = IntervalTreeNode[]
    print("\n\n")
    print("Query Interval (20,22,1) \n")
    findIntersectingIntervals!(query2, rootNode, results)
    for result in results
        print(result,"\n")
    end
    print("\n\n")

    #test delete interval
    #=
    query3 = Interval(5,5,1)
    results = IntervalTreeNode[]
    findIntersectingIntervals!(query3, rootNode, results)
    print("Tree Without Intervals Intersecting '[]' \n")
    for result in results
        print("Deleting: ", result, "\n")
        rootNode = deleteNodeFromIntervalTree!(result)
    end
    printTree(rootNode)
    =#

    #test split on point
    print("Tree Without 8 in it \n")
    rootNode = removePointFromIntervalTree!(rootNode, 8)
    print("Ordered Intervals From Tree \n")
    printTree(rootNode)
    print("\n\n")

    #test split on point
    print("Tree Without [5,10] in it \n")
    rootNode = removeIntervalFromIntervalTree!(rootNode, Interval(5,10))
    print("Ordered Intervals From Tree \n")
    printTree(rootNode)
    print("\n\n")

    print("Impose max of 15 \n")
    rootNode = imposeMaximumOnIntervalTree!(rootNode, 15)
    print("Ordered Intervals From Tree \n")
    printTree(rootNode)
    print("\n\n")

    print("Impose min of 2 \n")
    #rootNode = imposeMinimumOnIntervalTree!(rootNode, 2)
    print("Ordered Intervals From Tree \n")
    printTree(rootNode)
    print("\n\n")

    print("Impose open max of 14 \n")
    rootNode = imposeOpenMaximumOnIntervalTree!(rootNode, 14)
    print("Ordered Intervals From Tree \n")
    printTree(rootNode)
    print("\n\n")

    println("Impose open min of 3")
    rootNode = imposeOpenMinimumOnIntervalTree!(rootNode, 3)
    println("Ordered Intervals From Tree")
    printTree(rootNode)
    print("\n\n")

    println("Add the interval [8,15] (1) with merging")
    rootNode = addIntervalToIntervalTreeWithMerging!(rootNode, Interval(8,15))
    println("Ordered Intervals From Tree")
    printTree(rootNode)
    print("\n\n")

    println("Add the interval [26,30] (1) with merging")
    rootNode = addIntervalToIntervalTreeWithMerging!(rootNode, Interval(26,30))
    println("This is the 'Post Manipulation' Tree")
    printTree(rootNode)
    print("\n\n")

    println("Original Tree Copy")
    println("Ordered Intervals From Tree")
    printTree(treeCopy)
    print("\n\n")

    println("Original Tree Merged")
    println("Ordered Intervals From Tree")
    printTree(treeCopyMerged)
    print("\n\n")

    println("Union of Original and Post-Manipulation")
    println("Ordered Intervals From Tree")
    union = unionIntervalTrees(treeCopy, rootNode, false)
    printTree(union)
    print("\n\n")

    println("Union of Original and Post-Manipulation With Merging")
    println("Ordered Intervals From Tree")
    union2 = unionIntervalTrees(treeCopy, rootNode, true)
    printTree(union2)
    print("\n\n")

    ex1 = IntervalTreeNode(26,30,1)
    addNodeToIntervalTree!(ex1, IntervalTreeNode(15,25))
    addNodeToIntervalTree!(ex1, IntervalTreeNode(100,200,3))
    addNodeToIntervalTree!(ex1, IntervalTreeNode(300,400,6))

    ex2 = IntervalTreeNode(5,10,1)
    addNodeToIntervalTree!(ex2, IntervalTreeNode(15,25))
    addNodeToIntervalTree!(ex2, IntervalTreeNode(1,12,2))
    addNodeToIntervalTree!(ex2, IntervalTreeNode(8,16))
    addNodeToIntervalTree!(ex2, IntervalTreeNode(14,20))
    addNodeToIntervalTree!(ex2, IntervalTreeNode(18,21))
    addNodeToIntervalTree!(ex2, IntervalTreeNode(2,8))
    addNodeToIntervalTree!(ex2, IntervalTreeNode(90.4,150,1.6))
    addNodeToIntervalTree!(ex2, IntervalTreeNode(169,200,0.6))
    addNodeToIntervalTree!(ex2, IntervalTreeNode(310,400,10))

    println("Ordered Intervals From New Tree 1")
    printTree(ex1)
    print("\n\n")
    println("Ordered Intervals From New Tree 2")
    printTree(ex2)
    print("\n\n")

    println("Intersection of Previous 2 Trees")
    println("Ordered Intervals From Tree")
    intersection = intersectIntervalTrees(ex1, ex2)
    printTree(intersection)
    print("\n\n")

    println("Test of creating a tree from an array")
    a = [1,2,3,4,5,7,8,9,10]
    println(a)
    treeFromSet = setToIntervalTree(a)
    printTree(treeFromSet)
    print("\n\n")

    println("Test of creating a tree from an array")
    a = [1,2,3,4,5,7,8,9,11]
    println(a)
    treeFromSet = setToIntervalTree(a)
    printTree(treeFromSet)
    print("\n\n")

    println("Test of creating a tree from an array")
    a = [1,2,5,7,4,11,11.5,13.4,13.2]
    println(a)
    treeFromSet = setToIntervalTree(a)
    printTree(treeFromSet)
    print("\n\n")

    println("Test of creating a tree from an array")
    a = [1,1,1,1]
    println(a)
    treeFromSet = setToIntervalTree(a)
    printTree(treeFromSet)
    print("\n\n")

    println("Test of creating a tree from an array")
    a = [1,10]
    println(a)
    treeFromSet = setToIntervalTree(a)
    printTree(treeFromSet)
    print("\n\n")

    println("Test of creating a tree from an array")
    a = [1]
    println(a)
    treeFromSet = setToIntervalTree(a)
    printTree(treeFromSet)
    print("\n\n")

    println("Test of creating a tree from a set")
    a = Set((1,2,3,4,5,7,8,9,10))
    println(a)
    treeFromSet = setToIntervalTree(a)
    printTree(treeFromSet)
    print("\n\n")

    println("Test of disjoint union")
    b = Set((2,4,8,10,12))
    println(a)
    println(b)
    println("Expected Answer: {12} {1,3,5,7,9}")
    treeFromSet2 = setToIntervalTree(b)
    setA, setB = disjunctiveUnion(treeFromSet, treeFromSet2)
    println(setA)
    println(setB)
end