text:
unpack(x) = unpack(x, Any)

"""
    unpack(bytes, T::Type = Any)

Return `unpack(IOBuffer(bytes), T)`.
"""
unpack(bytes, ::Type{T}) where {T} = unpack(IOBuffer(bytes), T)::T

"""
    unpack(msgpack_byte_stream::IO, T::Type = Any)

Return the Julia value of type `T` deserialized from `msgpack_byte_stream`.

`T` is assumed to have valid [`msgpack_type`](@ref) and [`from_msgpack`](@ref)
definitions.

If `msgpack_type(T) === AnyType()`, `unpack` will deserialize the next
MessagePack object from `msgpack_byte_stream` into the default Julia
representation corresponding to the object's MessagePack type. For details on
default Julia representations, see [`AbstractMsgPackType`](@ref).

See also: [`pack`](@ref)
"""
unpack(io::IO, ::Type{T}) where {T} = unpack_type(io, read(io, UInt8), msgpack_type(T), T)::T

#####
##### `AnyType`
#####

function unpack_type(io, byte, t::AnyType, U::Union)
    A, B = U.a, U.b # Unions are sorted, so `Nothing`/`Missing` would be first
    if A === Nothing || A === Missing
        byte === magic_byte(NilFormat) && return from_msgpack(A, nothing)
        return unpack_type(io, byte, msgpack_type(B), B)
    end
    return _unpack_any(io, byte, U)
end

@inline unpack_type(io, byte, ::AnyType, T::Type) = _unpack_any(io, byte, T)

function _unpack_any(io, byte, ::Type{T}) where {T}
    if byte <= magic_byte_max(IntFixPositiveFormat)
        return unpack_format(io, IntFixPositiveFormat(byte), T)
    elseif byte <= magic_byte_max(MapFixFormat)
        return unpack_format(io, MapFixFormat(byte), T)
    elseif byte <= magic_byte_max(ArrayFixFormat)
        return unpack_format(io, ArrayFixFormat(byte), T)
    elseif byte <= magic_byte_max(StrFixFormat)
        return unpack_format(io, StrFixFormat(byte), T)
    elseif byte === magic_byte(UInt8Format)
        return unpack_format(io, UInt8Format(), T)
    elseif byte === magic_byte(UInt16Format)
        return unpack_format(io, UInt16Format(), T)
    elseif byte === magic_byte(UInt32Format)
        return unpack_format(io, UInt32Format(), T)
    elseif byte === magic_byte(UInt64Format)
        return unpack_format(io, UInt64Format(), T)
    elseif byte === magic_byte(Int8Format)
        return unpack_format(io, Int8Format(), T)
    elseif byte === magic_byte(Int16Format)
        return unpack_format(io, Int16Format(), T)
    elseif byte === magic_byte(Int32Format)
        return unpack_format(io, Int32Format(), T)
    elseif byte === magic_byte(Int64Format)
        return unpack_format(io, Int64Format(), T)
    elseif byte === magic_byte(Float32Format)
        return unpack_format(io, Float32Format(), T)
    elseif byte === magic_byte(Float64Format)
        return unpack_format(io, Float64Format(), T)
    elseif byte === magic_byte(Str8Format)
        return unpack_format(io, Str8Format(), T)
    elseif byte === magic_byte(Str16Format)
        return unpack_format(io, Str16Format(), T)
    elseif byte === magic_byte(Str32Format)
        return unpack_format(io, Str32Format(), T)
    elseif byte === magic_byte(TrueFormat)
        return unpack_format(io, TrueFormat(), T)
    elseif byte === magic_byte(FalseFormat)
        return unpack_format(io, FalseFormat(), T)
    elseif byte === magic_byte(NilFormat)
        return unpack_format(io, NilFormat(), T)
    elseif byte === magic_byte(Array16Format)
        return unpack_format(io, Array16Format(), T)
    elseif byte === magic_byte(Array32Format)
        return unpack_format(io, Array32Format(), T)
    elseif byte === magic_byte(Map16Format)
        return unpack_format(io, Map16Format(), T)
    elseif byte === magic_byte(Map32Format)
        return unpack_format(io, Map32Format(), T)
    elseif byte === magic_byte(Ext8Format)
        return unpack_format(io, Ext8Format(), T)
    elseif byte === magic_byte(Ext16Format)
        return unpack_format(io, Ext16Format(), T)
    elseif byte === magic_byte(Ext32Format)
        return unpack_format(io, Ext32Format(), T)
    elseif byte === magic_byte(ExtFix1Format)
        return unpack_format(io, ExtFix1Format(), T)
    elseif byte === magic_byte(ExtFix2Format)
        return unpack_format(io, ExtFix2Format(), T)
    elseif byte === magic_byte(ExtFix4Format)
        return unpack_format(io, ExtFix4Format(), T)
    elseif byte === magic_byte(ExtFix8Format)
        return unpack_format(io, ExtFix8Format(), T)
    elseif byte === magic_byte(ExtFix16Format)
        return unpack_format(io, ExtFix16Format(), T)
    elseif byte === magic_byte(Bin8Format)
        return unpack_format(io, Bin8Format(), T)
    elseif byte === magic_byte(Bin16Format)
        return unpack_format(io, Bin16Format(), T)
    elseif byte === magic_byte(Bin32Format)
        return unpack_format(io, Bin32Format(), T)
    elseif byte >= magic_byte_min(IntFixNegativeFormat)
        return unpack_format(io, IntFixNegativeFormat(reinterpret(Int8, byte)), T)
    else
        invalid_unpack(io, byte, AnyType(), T)
    end
    # TODO Ext*Format
end

#####
##### `ImmutableStructType`
#####

struct Constructor{T} end

(::Constructor{T})(args...) where {T} = construct(T, args...)

construct(::Type{T}, args...) where {T} = T(args...)

function unpack_type(io, byte, ::ImmutableStructType, ::Type{T}) where {T}
    constructor = Constructor{T}()
    N = fieldcount(T)
    byte > magic_byte_max(MapFixFormat) && read(io, UInt8)
    Base.@nexprs 32 i -> begin
        F_i = fieldtype(T, i)
        unpack_type(io, read(io, UInt8), StringType(), Skip{Symbol})
        x_i = unpack_type(io, read(io, UInt8), msgpack_type(F_i), F_i)
        N == i && return Base.@ncall i constructor x
    end
    others = Any[]
    for i in 33:N
        F_i = fieldtype(T, i)
        push!(others, unpack_type(io, read(io, UInt8), msgpack_type(F_i), F_i))
    end
    return constructor(x1, x2, x3, x4, x5, x6, x7, x8, x9, x10, x11, x12, x13,
                       x14, x15, x16, x17, x18, x19, x20, x21, x22, x23, x24, x25,
                       x26, x27, x28, x29, x30, x31, x32, others...)
end

function unpack_type(io, byte, ::ImmutableStructType, ::Type{Skip{T}}) where {T}
    N = fieldcount(T)
    byte > magic_byte_max(MapFixFormat) && read(io, UInt8)
    Base.@nexprs 32 i -> begin
        F_i = fieldtype(T, i)
        unpack_type(io, read(io, UInt8), StringType(), Skip{Symbol})
        unpack_type(io, read(io, UInt8), msgpack_type(F_i), Skip{F_i})
        N == i && return Skip{T}()
    end
    for i in 33:N
        F_i = fieldtype(T, i)
        unpack_type(io, read(io, UInt8), msgpack_type(F_i), Skip{F_i})
    end
    return Skip{T}()
end

#####
##### `MutableStructType`
#####

function unpack_type(io, byte, t::MutableStructType, ::Type{T}) where {T}
    N = fieldcount(T)
    if byte <= magic_byte_max(MapFixFormat)
        pair_count = xor(byte, magic_byte_min(MapFixFormat))
    elseif byte === magic_byte(Map16Format)
        pair_count = ntoh(read(io, UInt16))
    elseif byte === magic_byte(Map32Format)
        pair_count = ntoh(read(io, UInt32))
    else
        invalid_unpack(io, byte, t, T)
    end
    object = T()
    for _ in 1:pair_count
        key = unpack(io, Symbol) # TODO validation check?
        Base.@nif(32,
                  i -> i <= N && fieldname(T, i) === key,
                  i -> setfield!(object, i, unpack(io, fieldtype(T, i))),
                  i -> begin
                      is_field_still_unread = true
                      for j in 33:N
                          fieldname(T, j) === key || continue
                          setfield!(object, j, unpack(io, fieldtype(T, j)))
                          is_field_still_unread = false
                          break
                      end
                      is_field_still_unread && unpack(io)
                  end)
    end
    return object
end

unpack_type(io, byte, ::MutableStructType, ::Type{<:Skip}) = unpack_type(io, byte, MapType(), Skip{Dict{Any,Any}})

#####
##### `IntegerType`
#####

function unpack_type(io, byte, t::IntegerType, ::Type{T}) where {T}
    if byte <= magic_byte_max(IntFixPositiveFormat)
        return unpack_format(io, IntFixPositiveFormat(byte), T)
    elseif byte === magic_byte(UInt8Format)
        return unpack_format(io, UInt8Format(), T)
    elseif byte === magic_byte(UInt16Format)
        return unpack_format(io, UInt16Format(), T)
    elseif byte === magic_byte(UInt32Format)
        return unpack_format(io, UInt32Format(), T)
    elseif byte === magic_byte(UInt64Format)
        return unpack_format(io, UInt64Format(), T)
    elseif byte === magic_byte(Int8Format)
        return unpack_format(io, Int8Format(), T)
    elseif byte === magic_byte(Int16Format)
        return unpack_format(io, Int16Format(), T)
    elseif byte === magic_byte(Int32Format)
        return unpack_format(io, Int32Format(), T)
    elseif byte === magic_byte(Int64Format)
        return unpack_format(io, Int64Format(), T)
    elseif byte >= magic_byte_min(IntFixNegativeFormat)
        return unpack_format(io, IntFixNegativeFormat(reinterpret(Int8, byte)), T)
    else
        invalid_unpack(io, byte, t, T)
    end
end

unpack_format(io, f::IntFixPositiveFormat, ::Type{T}) where {T} = from_msgpack(T, UInt8(f.byte))
unpack_format(io, f::IntFixNegativeFormat, ::Type{T}) where {T} = from_msgpack(T, Int8(f.byte))

unpack_format(io, f::IntFixPositiveFormat, ::Type{T}) where {T<:Skip} = T()
unpack_format(io, f::IntFixNegativeFormat, ::Type{T}) where {T<:Skip} = T()

unpack_format(io, ::UInt8Format, ::Type{T}) where {T} = from_msgpack(T, read(io, UInt8))
unpack_format(io, ::UInt16Format, ::Type{T}) where {T} = from_msgpack(T, ntoh(read(io, UInt16)))
unpack_format(io, ::UInt32Format, ::Type{T}) where {T} = from_msgpack(T, ntoh(read(io, UInt32)))
unpack_format(io, ::UInt64Format, ::Type{T}) where {T} = from_msgpack(T, ntoh(read(io, UInt64)))

unpack_format(io, ::UInt8Format, ::Type{T}) where {T<:Skip} = (skip(io, 1); T())
unpack_format(io, ::UInt16Format, ::Type{T}) where {T<:Skip} = (skip(io, 2); T())
unpack_format(io, ::UInt32Format, ::Type{T}) where {T<:Skip} = (skip(io, 4); T())
unpack_format(io, ::UInt64Format, ::Type{T}) where {T<:Skip} = (skip(io, 8); T())

unpack_format(io, ::Int8Format, ::Type{T}) where {T} = from_msgpack(T, read(io, Int8))
unpack_format(io, ::Int16Format, ::Type{T}) where {T} = from_msgpack(T, ntoh(read(io, Int16)))
unpack_format(io, ::Int32Format, ::Type{T}) where {T} = from_msgpack(T, ntoh(read(io, Int32)))
unpack_format(io, ::Int64Format, ::Type{T}) where {T} = from_msgpack(T, ntoh(read(io, Int64)))

unpack_format(io, ::Int8Format, ::Type{T}) where {T<:Skip} = (skip(io, 1); T())
unpack_format(io, ::Int16Format, ::Type{T}) where {T<:Skip} = (skip(io, 2); T())
unpack_format(io, ::Int32Format, ::Type{T}) where {T<:Skip} = (skip(io, 4); T())
unpack_format(io, ::Int64Format, ::Type{T}) where {T<:Skip} = (skip(io, 8); T())

#####
##### `NilType`
#####

function unpack_type(io, byte, ::NilType, ::Type{T}) where {T}
    byte === magic_byte(NilFormat) && return unpack_format(io, NilFormat(), T)
    invalid_unpack(io, f, T)
end

unpack_type(io, byte, ::NilType, ::Type{T}) where {T<:Skip} = T()

unpack_format(io, ::NilFormat, ::Type{T}) where {T} = from_msgpack(T, nothing)

#####
##### `BooleanType`
#####

function unpack_type(io, byte, t::BooleanType, ::Type{T}) where {T}
    byte === magic_byte(TrueFormat) && return unpack_format(io, TrueFormat(), T)
    byte === magic_byte(FalseFormat) && return unpack_format(io, FalseFormat(), T)
    invalid_unpack(io, byte, t, T)
end

unpack_type(io, byte, ::BooleanType, ::Type{T}) where {T<:Skip} = T()

unpack_format(io, ::TrueFormat, ::Type{T}) where {T} = from_msgpack(T, true)

unpack_format(io, ::FalseFormat, ::Type{T}) where {T} = from_msgpack(T, false)

#####
##### `FloatType`
#####

function unpack_type(io, byte, t::FloatType, ::Type{T}) where {T}
    if byte === magic_byte(Float32Format)
        return unpack_format(io, Float32Format(), T)
    elseif byte === magic_byte(Float64Format)
        return unpack_format(io, Float64Format(), T)
    else
        invalid_unpack(io, byte, t, T)
    end
end

unpack_format(io, ::Float32Format, ::Type{T}) where {T} = from_msgpack(T, ntoh(read(io, Float32)))
unpack_format(io, ::Float32Format, ::Type{T}) where {T<:Skip} = (skip(io, 4); T())
unpack_format(io, ::Float64Format, ::Type{T}) where {T} = from_msgpack(T, ntoh(read(io, Float64)))
unpack_format(io, ::Float64Format, ::Type{T}) where {T<:Skip} = (skip(io, 8); T())

#####
##### `StringType`
#####

function unpack_type(io, byte, t::StringType, ::Type{T}) where {T}
    if byte <= magic_byte_max(StrFixFormat)
        return unpack_format(io, StrFixFormat(byte), T)
    elseif byte === magic_byte(Str8Format)
        return unpack_format(io, Str8Format(), T)
    elseif byte === magic_byte(Str16Format)
        return unpack_format(io, Str16Format(), T)
    elseif byte === magic_byte(Str32Format)
        return unpack_format(io, Str32Format(), T)
    else
        invalid_unpack(io, byte, t, T)
    end
end

unpack_format(io, ::Str8Format, ::Type{T}) where {T} = _unpack_string(io, read(io, UInt8), T)
unpack_format(io, ::Str16Format, ::Type{T}) where {T} = _unpack_string(io, ntoh(read(io, UInt16)), T)
unpack_format(io, ::Str32Format, ::Type{T}) where {T} = _unpack_string(io, ntoh(read(io, UInt32)), T)
unpack_format(io, f::StrFixFormat, ::Type{T}) where {T} = _unpack_string(io, xor(f.byte, magic_byte_min(StrFixFormat)), T)

_unpack_string(io, n, ::Type{T}) where {T} = from_msgpack(T, String(read(io, n)))
_unpack_string(io, n, ::Type{T}) where {T<:Skip} = (skip(io, n); T())

# Below is a nice optimization we can apply when `io` wraps an indexable byte
# buffer; this is really nice for skipping an extra copy when we deserialize
# the field names of `T` where `msgpack_type(T)::MutableStruct`. Like much of
# this package, this trick can be traced back to Jacob Quinn's magnificent
# JSON3.jl!
function _unpack_string(io::Base.GenericIOBuffer, n, ::Type{T}) where {T}
    result = from_msgpack(T, PointerString(pointer(io.data, io.ptr), n))
    skip(io, n)
    return result
end

# necessary to resolve ambiguity
_unpack_string(io::Base.GenericIOBuffer, n, ::Type{T}) where {T<:Skip} = (skip(io, n); T())

#####
##### `BinaryType`
#####

function unpack_type(io, byte, t::BinaryType, ::Type{T}) where {T}
    if byte === magic_byte(Bin8Format)
        return unpack_format(io, Bin8Format(), T)
    elseif byte === magic_byte(Bin16Format)
        return unpack_format(io, Bin16Format(), T)
    elseif byte === magic_byte(Bin32Format)
        return unpack_format(io, Bin32Format(), T)
    else
        invalid_unpack(io, byte, t, T)
    end
end

unpack_format(io, ::Bin8Format, ::Type{T}) where {T} = _unpack_binary(io, read(io, UInt8), T)
unpack_format(io, ::Bin16Format, ::Type{T}) where {T} = _unpack_binary(io, ntoh(read(io, UInt16)), T)
unpack_format(io, ::Bin32Format, ::Type{T}) where {T} = _unpack_binary(io, ntoh(read(io, UInt32)), T)

_unpack_binary(io, n, ::Type{T}) where {T} = from_msgpack(T, read(io, n))
_unpack_binary(io, n, ::Type{T}) where {T<:Skip} = (skip(io, n); T())

#####
##### `ArrayType`
#####

function unpack_type(io, byte, t::ArrayType, ::Type{T}) where {T}
    if byte <= magic_byte_max(ArrayFixFormat)
        return unpack_format(io, ArrayFixFormat(byte), T)
    elseif byte === magic_byte(Array16Format)
        return unpack_format(io, Array16Format(), T)
    elseif byte === magic_byte(Array32Format)
        return unpack_format(io, Array32Format(), T)
    else
        invalid_unpack(io, byte, t, T)
    end
end

unpack_format(io, ::Array16Format, ::Type{T}) where {T} = _unpack_array(io, ntoh(read(io, UInt16)), T)
unpack_format(io, ::Array32Format, ::Type{T}) where {T} = _unpack_array(io, ntoh(read(io, UInt32)), T)
unpack_format(io, f::ArrayFixFormat, ::Type{T}) where {T} = _unpack_array(io, xor(f.byte, magic_byte_min(ArrayFixFormat)), T)

function _unpack_array(io, n, ::Type{T}) where {T}
    E = eltype(T)
    e = msgpack_type(E)
    result = Vector{E}(undef, n)
    for i in 1:n
        result[i] = unpack_type(io, read(io, UInt8), e, E)
    end
    return from_msgpack(T, result)
end

function _unpack_array(io, n, ::Type{Skip{T}}) where {T}
    E = eltype(T)
    e = msgpack_type(E)
    for _ in 1:n
        unpack_type(io, read(io, UInt8), e, Skip{E})
    end
    return Skip{T}()
end

function _unpack_array(io::Base.GenericIOBuffer, n, ::Type{T}) where {T<:ArrayView}
    E = eltype(T)
    e = msgpack_type(E)
    start = position(io)
    positions = Vector{UInt64}(undef, n)
    for i in 1:length(positions)
        positions[i] = (position(io) - start) + 1
        unpack_type(io, read(io, UInt8), e, Skip{E})
    end
    bytes = view(io.data, (start + 1):position(io))
    return ArrayView{E}(bytes, positions)
end

#####
##### `MapType`
#####

function unpack_type(io, byte, t::MapType, ::Type{T}) where {T}
    if byte <= magic_byte_max(MapFixFormat)
        return unpack_format(io, MapFixFormat(byte), T)
    elseif byte === magic_byte(Map16Format)
        return unpack_format(io, Map16Format(), T)
    elseif byte === magic_byte(Map32Format)
        return unpack_format(io, Map32Format(), T)
    else
        invalid_unpack(io, byte, t, T)
    end
end

unpack_format(io, ::Map16Format, ::Type{T}) where {T} = _unpack_map(io, ntoh(read(io, UInt16)), T)
unpack_format(io, ::Map32Format, ::Type{T}) where {T} = _unpack_map(io, ntoh(read(io, UInt32)), T)
unpack_format(io, f::MapFixFormat, ::Type{T}) where {T} = _unpack_map(io, xor(f.byte, magic_byte_min(MapFixFormat)), T)

_keytype(T) = Any
_keytype(::Type{T}) where {K,V,T<:AbstractDict{K,V}} = keytype(T)
_valtype(T) = Any
_valtype(::Type{T}) where {K,V,T<:AbstractDict{K,V}} = valtype(T)

function _unpack_map(io, n, ::Type{T}) where {T}
    K = _keytype(T)
    k = msgpack_type(K)
    V = _valtype(T)
    v = msgpack_type(V)
    dict = Dict{K,V}()
    for _ in 1:n
        key = unpack_type(io, read(io, UInt8), k, K)
        val = unpack_type(io, read(io, UInt8), v, V)
        dict[key] = val
    end
    return from_msgpack(T, dict)
end

function _unpack_map(io, n, ::Type{Skip{T}}) where {T}
    K = _keytype(T)
    k = msgpack_type(K)
    V = _valtype(T)
    v = msgpack_type(V)
    for i in 1:n
        unpack_type(io, read(io, UInt8), k, Skip{K})
        unpack_type(io, read(io, UInt8), v, Skip{V})
    end
    return Skip{T}()
end

function _unpack_map(io::Base.GenericIOBuffer, n, ::Type{T}) where {T<:MapView}
    K = _keytype(T)
    k = msgpack_type(K)
    V = _valtype(T)
    v = msgpack_type(V)
    start = position(io)
    positions = Dict{K,UnitRange{UInt64}}()
    for _ in 1:n
        key = unpack_type(io, read(io, UInt8), k, K)
        value_start = (position(io) - start) + 1
        unpack_type(io, read(io, UInt8), v, Skip{V})
        positions[key] = value_start:(position(io) - start)
    end
    bytes = view(io.data, (start + 1):position(io))
    return MapView{K,V}(bytes, positions)
end

#####
##### `ExtensionType`
#####

function unpack_type(io, byte, t::ExtensionType, ::Type{T}) where {T}
    byte === magic_byte(ExtFix1Format) && return unpack_format(io, ExtFix1Format(), T)
    byte === magic_byte(ExtFix2Format) && return unpack_format(io, ExtFix2Format(), T)
    byte === magic_byte(ExtFix4Format) && return unpack_format(io, ExtFix4Format(), T)
    byte === magic_byte(ExtFix8Format) && return unpack_format(io, ExtFix8Format(), T)
    byte === magic_byte(ExtFix16Format) && return unpack_format(io, ExtFix16Format(), T)
    byte === magic_byte(Ext8Format) && return unpack_format(io, Ext8Format(), T)
    byte === magic_byte(Ext16Format) && return unpack_format(io, Ext16Format(), T)
    byte === magic_byte(Ext32Format) && return unpack_format(io, Ext32Format(), T)
    invalid_unpack(io, byte, t, T)
end

unpack_format(io, ::ExtFix1Format, ::Type{T}) where {T} = _unpack_extension(io, 1, T)
unpack_format(io, ::ExtFix2Format, ::Type{T}) where {T} = _unpack_extension(io, 2, T)
unpack_format(io, ::ExtFix4Format, ::Type{T}) where {T} = _unpack_extension(io, 4, T)
unpack_format(io, ::ExtFix8Format, ::Type{T}) where {T} = _unpack_extension(io, 8, T)
unpack_format(io, ::ExtFix16Format, ::Type{T}) where {T} = _unpack_extension(io, 16, T)
unpack_format(io, ::Ext8Format, ::Type{T}) where {T} = _unpack_extension(io, read(io, UInt8), T)
unpack_format(io, ::Ext16Format, ::Type{T}) where {T} = _unpack_extension(io, ntoh(read(io, UInt16)), T)
unpack_format(io, ::Ext32Format, ::Type{T}) where {T} = _unpack_extension(io, ntoh(read(io, UInt32)), T)

function _unpack_extension(io, n, ::Type{T}) where {T}
    type = read(io, Int8)
    data = read(io, n)
    return from_msgpack(T, Extension(type, data))
end

#####
##### utilities
#####

@noinline function invalid_unpack(io, byte, m, T)
    error("invalid byte $(repr(byte)) encountered in $(io) attempting to read a MsgPack $(m) into a Julia $(T) at position $(position(io))")
end


text:
using Flux, PGFPlotsX, DelimitedFiles
using BSON: @load
using Statistics

"""
PlotParams holds information needed for plotting, e.g. type of plot, title, x and y labels, directory to save figure, etc.
"""
mutable struct PlotParams
    plot::Function
    title
    x_label
    y_label
    model_label
    x_min
    x_max
    x_range
    num_xticks
    y_min
    y_max
    y_range
    num_yticks
    err
    fig_dir
    colorbar_label
end

"""
PlotParams custom constuctor
"""
function new_params(plot;
                     title=nothing,
                     x_label=nothing,
                     y_label=nothing,
                     model_label=nothing,
                     x_min=nothing,
                     x_max = nothing,
                     x_range = nothing,
                     num_xticks = nothing,
                     y_min = nothing,
                     y_max = nothing,
                     y_range = nothing,
                     num_yticks=nothing,
                     err=nothing,
                     fig_dir = nothing,
                     colorbar_label=nothing)
    return PlotParams(plot, title, x_label, y_label, model_label, x_min, x_max, x_range, num_xticks, y_min, y_max, y_range, num_yticks, err, fig_dir, colorbar_label)
end

"""
Generate plots for each set of PlotParams given in params argument.
"""
function make_plots(model, data, plots; load_data=false, load_model=false, η=identity)
    load_model ? (@load model m) : local m = model;

    if load_data
        train_set, test_set = import_split(data_dir)
        augment!(test_set)
        test_x, gt = unpackage(test_set)
        gt = euler_to_index(gt)
    else
        (test_x, gt) = data
    end
    pred = η(map(argmax, eachcol(m(test_x))))
    for p ∈ plots
        if isnothing(p.err)
            p.plot((gt, pred), params=p)
        else
            p.plot((gt, pred), err=p.err, params=p)
        end
    end
end

"""
Plot confusion matrix for data
"""
## Confusion matrix
function plot_cmatrix( data;
                       bound = CLASS_NUM - 1,
                       k=nothing,
                       params = defaultParams,
                       fig_dir=nothing)

    !isnothing(params.x_range) && (bound = params.x_range)
    if isnothing(k)
        (ontology_index == 2) && (k = 480)
        (ontology_index == 3) && (k = 480)
        (ontology_index == 4) && (k = 961)
    end
    (gt, pred) = data
    if ontology_index == 4
        a = euler_to_index(-bound)
        b = euler_to_index(bound)
        gt = euler_to_index(gt)
        pred = euler_to_index(pred)
    else
        a = 1
        b = bound
    end
    cm = confusion_matrix(k, gt, pred)
    cm = cm[a:b, a:b]
    x = collect(a:b)
    y = x;

    plot_heatmap(x, y, cm; params=params)
end

"""
Plot a PGFPlots heatmap.
"""
function plot_heatmap(x, y, z; params=defaultParams)

    if !isnothing(params.fig_dir)
        fig_dir = params.fig_dir
    else
        @error("No figure directory given.")
        return nothing
    end

    axis = @pgf Axis(
        {
            title=params.title,
            xlabel=params.x_label,
            ylabel=params.y_label,
            view = (0, 90),
            colorbar,
            "colormap/jet",
            "point meta max = 1000",
            colorbar_style = {
                                ylabel=params.colorbar_label
                             },
        },
        Plot3(
            {
                surf,
                shader = "flat",
            },
            Coordinates(x, y, z)
        )
    )
    !isnothing(fig_dir) && pgfsave(fig_dir, axis)
end

"""
Plot error on data. Error calculated by error function argument (defaults to relative error).
"""
## Error
function plot_error(data;
                   err = rel_error,
                   params = defaultParams,
                   x_min = -960,
                   x_max = 960,
                   x_range = nothing,
                   num_xticks = 9,
                   y_range = nothing,
                   y_min = nothing,
                   y_max = nothing,
                   num_yticks= 9,
                   fig_dir=nothing)
    (gt, pred) = data
    error, μₑ, σₑ = error_stats(gt, pred, err)
    max_error = maximum(error)
    min_error = minimum(error)
    x_min, x_max, y_min, y_max, num_xticks, num_yticks = _maybe_useparams(params, x_min, x_max, y_min, y_max, num_xticks, num_yticks)
    (isnothing(fig_dir) && !isnothing(params.fig_dir)) && (fig_dir = params.fig_dir)
    isnothing(y_min) && (y_min = min_error)
    isnothing(y_max) && (y_max = max_error)
    !isnothing(x_range) && (x_min = -abs(x_range)) & (x_max = abs(x_range))
    !isnothing(y_range) && (y_min = -abs(y_range)) & (y_max = abs(y_range))
    axis = @pgf Axis(
        {
            title=params.title,
            xlabel=params.x_label,
            ylabel=params.y_label,
            xmin = x_min,
            xmax = x_max,
            xtick = collect(range(x_min, x_max; length=num_xticks)),
            ytick = collect(range(y_min, y_max; length=num_yticks)),
            ymin = y_min,
            ymax = y_max,
            "minor y tick num = 9",
            "every axis/.append style" = {font = raw"\footnotesize" },
            "every node" = {font = raw"2em"},
            no_marks,
            legend_style={nodes={"scale=0.5"}, draw="none"},
        },
        Plot(
            {
                "only marks",
                mark_size = "0.1pt",
                solid,
                color => "black"
            },
            Table(;
                x=gt,
                y=error
            )
        ),
        LegendEntry(params.model_label),
        Plot({no_marks, mark_size="0.1pt", color="white"},Table(;x=x_min,y=y_min)),LegendEntry("\\mu = $μₑ"),
        Plot({no_marks, mark_size="0.1pt", color="white"},Table(;x=x_min,y=y_min)),LegendEntry("\\sigma = $σₑ"),
    )
    !isnothing(fig_dir) && pgfsave(fig_dir, axis)
end

"""
Plot error with colorbar. Color denotes data[2]
"""
## Error
function plot3_error(data;
                    err = rel_error,
                    params = defaultParams,
                    x_min = -960,
                    x_max = 960,
                    x_range = nothing,
                    num_xticks = 9,
                    y_max = nothing,
                    fig_name=nothing)
    colorbar_label = ""
    (gt,pred) = data
    if ontology_index == 4
        valid = index_to_euler(gt) .!= 0
        gt = index_to_euler(gt[valid])
        pred = index_to_euler(pred[valid])
        colorbar_label = "Predicted Euler"
    elseif ontology_index == 3
        colorbar_label = "Predicted \$h^{2,1}\$"
    elseif ontology_index == 2
        colorbar_label = "Predicted \$h^{1,1}\$"
    end
    error, μₑ, σₑ = error_stats(gt, pred, err)
    !isnothing(x_range) && (x_min = -abs(x_range)) & (x_max = abs(x_range))
    axis = @pgf Axis(
        {
            title=params.title,
            xlabel=params.x_label,
            ylabel=params.y_label,
            xmin = x_min,
            xmax = x_max,
            xtick = collect(range(x_min, x_max; length=num_xticks+1)),
            ymin = 0,
            ymax = y_max,
            "minor y tick num = 3",
            raw"every axis/.append style" = {font = raw"\footnotesize" },
            view = "{0}{90}",
            colorbar,
            "colormap/jet",
            colorbar_style = {
                                ylabel=colorbar_label
                             },
            legend_style={nodes={"scale=0.5"}, draw="none"},
        },
        Plot3(
            {
                "only marks",
                scatter,
                mark_size = "0.1pt",
                solid,
                color => "black"
            },
            Table(;
                x=gt,
                y=error,
                z=pred
            )
        ),
        LegendEntry(params.model_label),
        Plot({only_marks, mark_size="0.1pt", color="white"},Table(;x=0,y=0)),LegendEntry("\\mu = $μₑ"),
        Plot({only_marks, mark_size="0.1pt", color="white"},Table(;x=0,y=0)),LegendEntry("\\sigma = $σₑ"),
    )
    !isnothing(fig_name) && pgfsave("../../figs/$fig_name", axis)
end

function plot_error_vs_firstorder(data;
                   err = rel_error,
                   params = nothing,
                   x_min = -960,
                   x_max = 960,
                   x_range = nothing,
                   num_xticks = 9,
                   y_range = nothing,
                   y_min = nothing,
                   y_max = nothing,
                   num_yticks= 9,
                   fig_dir=nothing)
    (gt, X, pred) = data
    fo = calculate_leadingpiece(X)
    error, μₑ, σₑ = error_stats(gt, pred, err)
    error_fo, μ_fo, σ_fo = error_stats(fo, pred, err)
    max_error = maximum(error)
    min_error = minimum(error)
    x_min, x_max, y_min, y_max, num_xticks, num_yticks = _maybe_useparams(params, x_min, x_max, y_min, y_max, num_xticks, num_yticks)
    (isnothing(fig_dir) && !isnothing(params.fig_dir)) && (fig_dir = params.fig_dir)
    isnothing(y_min) && (y_min = min_error)
    isnothing(y_max) && (y_max = max_error)
    !isnothing(x_range) && (x_min = -abs(x_range)) & (x_max = abs(x_range))
    !isnothing(y_range) && (y_min = -abs(y_range)) & (y_max = abs(y_range))
    axis = @pgf Axis(
        {
            title=params.title,
            xlabel=params.x_label,
            ylabel=params.y_label,
            xmin = x_min,
            xmax = x_max,
            xtick = collect(range(x_min, x_max; length=num_xticks)),
            ytick = collect(range(y_min, y_max; length=num_yticks)),
            ymin = y_min,
            ymax = y_max,
            "minor y tick num = 9",
            "every axis/.append style" = {font = raw"\footnotesize" },
            "every node" = {font = raw"2em"},
            no_marks,
            legend_style={nodes={"scale=0.5"}, draw="none"},
        },
        Plot(
            {
                "only marks",
                mark_size = "0.1pt",
                solid,
                color => "blue",
                opacity = 0.01,
            },
            Table(;
                x=gt,
                y=error
            )
        ),
        Plot(
            {
                "only marks",
                mark_size = "0.1pt",
                solid,
                color => "red",
                opacity = 0.01,
            },
            Table(;
                x=fo,
                y=error_fo
            )
        ),
    )
    !isnothing(fig_dir) && pgfsave(fig_dir, axis, dpi=600)
end

"""
Plot data as a histogram
"""
## Distributions
function plot_datahist(data, ontology;
                   x_bound= CLASS_NUM - 1,
                   params = defaultParams,
                   fig_dir=nothing)
    (gt, p) = data
    if ontology == "euler"
        local o_range = -960:2:960
    elseif ontology == "h11" || ontology == "h21"
        local o_range = 1:480
        local x_ticks = [1,100,100,200,200,300,400,480]
    end
    number_density = [count(x->x==i, gt) for i ∈ o_range]
    bins = [(x,y) for (x,y) ∈ zip(collect(o_range), number_density)]
    (isnothing(fig_dir) && !isnothing(params.fig_dir)) && (fig_dir = params.fig_dir)
    axis = @pgf Axis(
        {
            title=params.title,
            xlabel=params.x_label,
            ylabel=params.y_label,
            xmin = params.x_min,
            xmax = params.x_max,
            ymin = params.y_min,
            xtick = x_ticks,
            grid = "none",
            "minor y tick num = 3",
        },
        Plot(
            {
                "mark=no",
                "ybar interval",
                fill = "black",
                color = "black",
            },
            Coordinates(bins)
        )
    )
    !isnothing(fig_dir) && pgfsave(fig_dir, axis)
end

"""
Plot data as a histogram
"""
## Distributions
function plot_datahist_with_prediction(data, ontology;
                   x_bound= CLASS_NUM - 1,
                   params = defaultParams,
                   fig_dir=nothing)
    (gt, p) = data
    if ontology == "euler"
        local o_range = -960:2:960
    elseif ontology == "h11" || ontology == "h21"
        local o_range = 1:480
    end
    number_density_gt = [count(x->x==i, gt) for i ∈ o_range]
    bins_gt = [(x,y) for (x,y) ∈ zip(collect(o_range), number_density_gt)]
    number_density_p = [count(x->x==i, p) for i ∈ o_range]
    bins_p = [(x,y) for (x,y) ∈ zip(collect(o_range), number_density_p)]
    (isnothing(fig_dir) && !isnothing(params.fig_dir)) && (fig_dir = params.fig_dir)

    axis = @pgf Axis(
        {
            title=params.title,
            xlabel=params.x_label,
            ylabel=params.y_label,
            xmin = params.x_min,
            xmax = params.x_max,
            ymin = params.y_min,
            xtick = collect(-900:25:900),
            grid = "none",
            "minor y tick num = 3",
        },
        Plot(
            {
                "mark=no",
                "ybar interval",
                fill = "blue",
                color = "blue",
                opacity = 0.5,
            },
            Coordinates(bins_gt)
        ),
        Plot(
            {
                "mark=no",
                "ybar interval",
                fill = "red",
                color = "red",
                opacity = 0.5,
            },
            Coordinates(bins_p)
        )
    )
    !isnothing(fig_dir) && pgfsave(fig_dir, axis)
end

"""
Plot mirror symmetry of data: h₁,₁ + h₂,₁ vs h₁,₁ - h₂,₁
"""
function plot_mirror(data;
                     params=defaultParams,
                     x_min=nothing,
                     x_max=nothing,
                     x_range = nothing,
                     num_xticks = 9,
                     y_range = nothing,
                     y_min = nothing,
                     y_max = nothing,
                     num_yticks= 9,
                     fig_dir=nothing)
    (ontology_index == 4) && (return @info("Ontology not set to hodge."))
    (h11,h21) = data
    xs = h11 .- h21
    ys = h11 .+ h21
    (isnothing(fig_dir) && !isnothing(params.fig_dir)) && (fig_dir = params.fig_dir)

    # x_min, x_max, y_min, y_max, num_xticks, num_yticks = _maybe_useparams(params,x_min, x_max, y_min, y_max, num_xticks, num_yticks)
    axis = @pgf Axis(
        {
            title=params.title,
            xlabel=params.x_label,
            ylabel=params.y_label,
            # xmin = x_min,
            # xmax = x_max,
            # xtick = collect(range(x_min, x_max; length=num_xticks)),
            # ytick = collect(range(y_min, y_max; length=num_yticks)),
            # ymin = y_min,
            # ymax = y_max,
            "minor y tick num = 9",
            "every axis/.append style" = {font = raw"\footnotesize" },
            "every node" = {font = raw"2em"},
            no_marks,
            legend_style={nodes={"scale=0.5"}, draw="none"},
        },
        Plot(
            {
                "only marks",
                mark_size = "0.05pt",
                solid,
                fill => "black",
                color => "black"
            },
            Table(;
                x=xs,
                y=ys
            )
        ),
        LegendEntry(params.model_label),
    )
    !isnothing(fig_dir) && pgfsave(fig_dir, axis)
end

"""
_maybe_useparams decides which params to use
"""
function _maybe_useparams(params, xmin, xmax, ymin, ymax, numxticks, numyticks)
    x_min = xmin
    x_max = xmax
    y_min = ymin
    y_max = ymax
    num_xticks = numxticks
    num_yticks = numyticks
    !isnothing(params.x_range) && (x_min = -abs(params.x_range)) & (x_max = abs(params.x_range))
    !isnothing(params.y_range) && (y_min = -abs(params.y_range)) & (y_max = abs(params.y_range))
    !isnothing(params.x_min) && (x_min = params.x_min)
    !isnothing(params.x_max) && (x_max = params.x_max)
    !isnothing(params.y_min) && (y_min = params.y_min)
    !isnothing(params.y_max) && (y_max = params.y_max)
    !isnothing(params.num_xticks) && (num_xticks = params.num_xticks)
    !isnothing(params.num_yticks) && (num_yticks = params.num_yticks)
    !isnothing(params.num_xticks) && (num_xticks = params.num_xticks)
    return x_min, x_max, y_min, y_max, num_xticks, num_yticks
end

ext: jl


text:
__precompile__()

module Maple

include("attributetypes.jl")
include("macros.jl")
include("map_coder.jl")
include("metadata.jl")
include("side.jl")

@exportalltypes

# export functions manually
export encodeMap, loadMap, encodeSide, loadSide, entityMap, getRoomByName, getSideName, getPreferredType
export @pardef, @kwproxy, @mapdef, @valueequals

end

max_stars_repo_name: CelestialCartographers/Maple

max_issues_repo_name: CelestialCartographers/Maple

max_forks_repo_name: CelestialCartographers/Maple


text:
using CartesianFDM

n = (6, 6)
bc = ntuple(dirichlet, length(n))

ops = fdmoperators(bc, n)

### symbolic
ε = only(@variables(ε))

# metrics
V = replace!(mask(ops, scalar(:V, n))) do el
    iszero(el) ? ε : el
end
A = mask(ops, vector(:A, n))

# primary
T = scalar(:T, n)
D = scalar(:D, n)

# secondary
G = gradient(ops, A, V, T, D)
L = divergence(ops, A, G)

# remove ε
L = substitute.(L, Ref(Dict([ε => 0])))

### numerical
V̂ = map(mask(ops, ones(Float64, prod(n)))) do x
    iszero(x) ? eps(x) : x
end
Â = mask(ops, ntuple(length(n)) do i
                  ones(Bool, prod(n))
              end)
D̂ = zeros(Bool, prod(n))

Ĝ = gradient(ops, Â, V̂, T, D̂)
L̂ = divergence(ops, Â, Ĝ)

###
using SparseArrays

Ĵ = sparse(Symbolics.jacobian(L̂, T))
M = eval(first(build_function(Ĵ)))()

T̂ = rand(prod(n))

b̂₁ = M * T̂
b̂₂ = eval(first(build_function(L̂, T)))(T̂)

using LinearAlgebra

@show norm(b̂₁ - b̂₂)

#=
Unsteady heat conduction on a square
(homogeneous B.C. is non-periodic).

=#
#=
@variables t p
(rhs, rhs!) = build_function(L̂, T, p, t)

using DifferentialEquations

tspan = (0.0, 100.0)
T̂ = ones(prod(n))
prob = ODEProblem(eval(rhs!), T̂, tspan)
sol = solve(prob)

using Plots

heatmap(reshape(sol(0.1), n...))

=#


text:
@testset "VAE" begin

    using Distributions

    function vae_cost(We1, be1, We2, be2, We3, be3, We4, be4, Wd1, bd1, Wd2, bd2, Wd3, bd3, eps, x)
        dummy = 42.0
        # encoder
        he1 = tanh.(We1 * x .+ be1)
        he2 = tanh.(We2 * he1 .+ be2)
        mu = We3 * he2 .+ be3
        log_sigma2 = We4 * he2 .+ be4
        z = mu .+ sqrt.(exp.(log_sigma2)) .* eps
        # decoder
        hd1 = tanh.(Wd1 * z .+ bd1)
        hd2 = tanh.(Wd2 * hd1 .+ bd2)
        x_rec = logistic.(Wd3 * hd2 .+ bd3)
        # loss
        rec_loss = sum(x .* log.(1e-10 .+ x_rec) + (1 .- x) .* log.(1e-10 .+ 1 .- x_rec), dims=1)
        latent_loss = -0.5 * sum(1 .+ log_sigma2 .- mu .^ 2 .- exp.(log_sigma2), dims=1)
        cost = mean(rec_loss .+ latent_loss)
    end

    function xavier_init(dim_in, dim_out; c=1)
        low = -c * sqrt(6.0 / (dim_in + dim_out))
        high = c * sqrt(6.0 / (dim_in + dim_out))
        return rand(Uniform(low, high), dim_in, dim_out)
    end

    let
        We1 = xavier_init(500, 784); be1 = randn(500);
        We2 = xavier_init(500, 500); be2 = randn(500);
        We3 = xavier_init(20, 500); be3 = randn(20);
        We4 = xavier_init(20, 500); be4 = randn(20);
        Wd1 = xavier_init(500, 20); bd1 = randn(500);
        Wd2 = xavier_init(500, 500); bd2 = randn(500);
        Wd3 = xavier_init(784, 500); bd3 = randn(784);
        x = rand(784, 100); eps = rand(Normal(0, 1),  20)
        inputs = [:We1 => We1, :be1 => be1, :We2 => We2, :be2 => be2,
                  :We3 => We3, :be3 => be3, :We4 => We4, :be4 => be4,
                  :Wd1 => Wd1, :bd1 => bd1, :Wd2 => Wd2, :bd2 => bd2,
                  :Wd3 => Wd3, :bd3 => bd3, :eps => eps, :x => x]

        test_compare(vae_cost; inputs...)
    end
end


text:
# Copyright (c) 2021 California Institute of Technology (“Caltech”).
# U.S. Government sponsorship acknowledged.
#
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions are met:
#
# - Redistributions of source code must retain the above copyright notice, this
#   list of conditions and the following disclaimer.
#
# - Redistributions in binary form must reproduce the above copyright notice,
#   this list of conditions and the following disclaimer in the documentation
#   and/or other materials provided with the distribution.
#
# - Neither the name of Caltech nor its operating division, the Jet Propulsion
#   Laboratory, nor the names of its contributors may be used to endorse or
#   promote products derived from this software without specific prior written
#   permission.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS
# IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO,
# THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
# PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR
# CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
# EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
# PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS;
# OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY,
# WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR
# OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF
# ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

#################### 3D cryowindow #######################33

function find_acceptable_tolerance(Hrow, htol)
    Hsort = sort(abs.(Hrow))  # sort ascending
    Hrowtotal = sum(Hsort)
    s = 0.0
    i = length(Hsort)
    while 1 - s / Hrowtotal > htol && i > 1
        s += Hsort[i]
        #@show s,Hrowtotal,s/Hrowtotal,i
        i -= 1
    end
    tolerance = Hsort[i]
    return tolerance
end

# This version writes one large sparse matrix
function cryo_window(cfb; htol=0.01, do_calc=true, lmax=Inf)
    @time ℓ, H = angular_cryo_window(cfb; do_calc, lmax)  # ℓ contains full-sky (integer) ℓ
    @assert all(isfinite.(H))

    ell = Vector{Int}[]  # will contain continuous ℓ
    k = Vector{Float64}[] # will contain continuous k
    fulllen = 0
    @time for l in ℓ
        # find a closely matching L ≈ l, as that is most likely to give relevant k.
        iL = argmin(@. abs(l - cfb.angfab.ell))
        L = cfb.angfab.ell[iL]

        kj = get_kj(cfb, L)
        push!(ell, fill(l, length(kj)))
        push!(k, kj)
        fulllen += length(kj)
    end
    @show fulllen

    cryolen = 0
    I = Int[]
    J = Int[]
    V = Float64[]
    @time for i=1:size(H,1)
        #L = cfb.angfab.ell[i]  # cryo-ℓ
        #kN = get_kN(cfb, L)  # cryo-k
        kN_len = length(cfb.reff)
        cryolen += kN_len

        (!do_calc) && continue
        (cfb.angfab.ell[i] > lmax) && continue

        tolerance = find_acceptable_tolerance(H[i,:], htol)
        @show i,cfb.angfab.ell[i],sum(abs.(H[i,:])),tolerance
        flen = 0
        @time for j=1:size(H,2)
            flen += length(k[j])
            if abs(H[i,j]) >= tolerance
                ibase = cryolen - kN_len
                jbase = flen - length(k[j])
                Δkj = k[j][2] - k[j][1]
                b̃ll = radial_cryo_window2(cfb, i, ℓ[j], k[j])
                b̃ll2 = sparsify(abs2.(b̃ll), rtol=htol)
                bI, bJ, bV = findnz(b̃ll2)
                @assert all(isfinite.(bV))
                append!(I, ibase .+ bI)
                append!(J, jbase .+ bJ)
                append!(V, H[i,j] * Δkj .* bV)
            end
        end
    end
    W = sparse(I, J, V, cryolen, fulllen)
    k = [kj for kv in k for kj in kv]
    ell = [lj for lv in ell for lj in lv]
    return k, ell, W
end

function sparsify(m; rtol=0.01)
    for i=1:size(m,1)
        row = m[i,:]
        tolerance = find_acceptable_tolerance(row, rtol)
        s = @. abs(row) < tolerance
        m[i,s] .= 0
    end
    return sparse(m)
end

# vim: set sw=4 et sts=4 :


text:
# Mahalanobis distances

immutable Mahalanobis{T} <: Metric
    qmat::Matrix{T}
end

immutable SqMahalanobis{T} <: SemiMetric
    qmat::Matrix{T}
end

result_type{T}(::Mahalanobis{T}, ::AbstractArray, ::AbstractArray) = T
result_type{T}(::SqMahalanobis{T}, ::AbstractArray, ::AbstractArray) = T

# SqMahalanobis

function evaluate{T<:AbstractFloat}(dist::SqMahalanobis{T}, a::AbstractVector, b::AbstractVector)
    Q = dist.qmat
    z = a - b
    return dot(z, Q * z)
end

sqmahalanobis(a::AbstractVector, b::AbstractVector, Q::AbstractMatrix) = evaluate(SqMahalanobis(Q), a, b)

function colwise!{T<:AbstractFloat}(r::AbstractArray, dist::SqMahalanobis{T}, a::AbstractMatrix, b::AbstractMatrix)
    Q = dist.qmat
    m, n = get_colwise_dims(size(Q, 1), r, a, b)
    z = a - b
    dot_percol!(r, Q * z, z)
end

function colwise!{T<:AbstractFloat}(r::AbstractArray, dist::SqMahalanobis{T}, a::AbstractVector, b::AbstractMatrix)
    Q = dist.qmat
    m, n = get_colwise_dims(size(Q, 1), r, a, b)
    z = a .- b
    Qz = Q * z
    dot_percol!(r, Q * z, z)
end

function pairwise!{T<:AbstractFloat}(r::AbstractMatrix, dist::SqMahalanobis{T}, a::AbstractMatrix, b::AbstractMatrix)
    Q = dist.qmat
    m, na, nb = get_pairwise_dims(size(Q, 1), r, a, b)

    Qa = Q * a
    Qb = Q * b
    sa2 = dot_percol(a, Qa)
    sb2 = dot_percol(b, Qb)
    At_mul_B!(r, a, Qb)

    for j = 1 : nb
        @simd for i = 1 : na
            @inbounds r[i,j] = sa2[i] + sb2[j] - 2 * r[i,j]
        end
    end
    r
end

function pairwise!{T<:AbstractFloat}(r::AbstractMatrix, dist::SqMahalanobis{T}, a::AbstractMatrix)
    Q = dist.qmat
    m, n = get_pairwise_dims(size(Q, 1), r, a)

    Qa = Q * a
    sa2 = dot_percol(a, Qa)
    At_mul_B!(r, a, Qa)

    for j = 1 : n
        for i = 1 : j-1
            @inbounds r[i,j] = r[j,i]
        end
        r[j,j] = 0
        for i = j+1 : n
            @inbounds r[i,j] = sa2[i] + sa2[j] - 2 * r[i,j]
        end
    end
    r
end

# Mahalanobis

function evaluate{T<:AbstractFloat}(dist::Mahalanobis{T}, a::AbstractVector, b::AbstractVector)
    sqrt(evaluate(SqMahalanobis(dist.qmat), a, b))
end

mahalanobis(a::AbstractVector, b::AbstractVector, Q::AbstractMatrix) = evaluate(Mahalanobis(Q), a, b)

function colwise!{T<:AbstractFloat}(r::AbstractArray, dist::Mahalanobis{T}, a::AbstractMatrix, b::AbstractMatrix)
    sqrt!(colwise!(r, SqMahalanobis(dist.qmat), a, b))
end

function colwise!{T<:AbstractFloat}(r::AbstractArray, dist::Mahalanobis{T}, a::AbstractVector, b::AbstractMatrix)
    sqrt!(colwise!(r, SqMahalanobis(dist.qmat), a, b))
end

function pairwise!{T<:AbstractFloat}(r::AbstractMatrix, dist::Mahalanobis{T}, a::AbstractMatrix, b::AbstractMatrix)
    sqrt!(pairwise!(r, SqMahalanobis(dist.qmat), a, b))
end

function pairwise!{T<:AbstractFloat}(r::AbstractMatrix, dist::Mahalanobis{T}, a::AbstractMatrix)
    sqrt!(pairwise!(r, SqMahalanobis(dist.qmat), a))
end


text:
# Base for all Stationary Agents
"""
    StationaryAgentBase

`StationaryAgentBase` is the base type for all agents designed to handle stationary
arms. Even though it is not enforced, the assumption is the reward updated to
stationary agents are from a stationary arm model.
"""
abstract type StationaryAgentBase <: AgentBase end

"""
    get_arm_index( agent::StationaryAgentBase )

Gets the index of next arm to pull.

## Example
```julia-repl

```
"""
get_arm_index( agent::StationaryAgentBase ) = nothing

"""
    update_reward!( agent::StationaryAgentBase, r::Real )

Updates the reward to bandit algorithm `agent`. `r` must a real number within valid range.
"""
update_reward!( agent::StationaryAgentBase, r::Real ) = nothing

include( "EpsilonGreedy.jl" )


text:
# symbols or what
struct Measurement{N}

end

max_stars_repo_name: amiragha/mps

max_issues_repo_name: amiragha/mps

max_forks_repo_name: amiragha/mps


text:
variable_primal(model) = x -> MOI.get(model, MOI.VariablePrimal(), x)

function heuristicCreate(model::Model)
    m = heuristicModel()
    m.model = model
    return m
end

function removeConstrains!(m::heuristicModel)
    model = m.model
    functions = m.functions
    sets = m.sets
    violation = m.violation
    functionValue = m.functionValue

    for (F, S) in list_of_constraint_types(model)
        if F != VariableRef
            for con in all_constraints(model, F, S)
                f = JuMP.MOI.get(model, JuMP.MOI.ConstraintFunction(), con)
                s = JuMP.MOI.get(model, JuMP.MOI.ConstraintSet(), con)
                push!(functions,f)
                push!(sets,s)
                push!(violation,Inf)
                push!(functionValue,0.0)
                delete(model,con)
                unregister(model,:con)
            end
        end
    end
    nothing
end

function updateFunctionValue!(m::heuristicModel)
    model = m.model
    functions = m.functions
    functionValue = m.functionValue

    for i in 1:length(functions)
        v = functions[i].constant
        for t in functions[i].terms
            coef = t.coefficient
            indx = t.variable_index.value
            v += coef * value(all_variables(model)[indx])
        end
        functionValue[i] = v
    end
    nothing
end

function updateViolation!(m::heuristicModel)
    sets = m.sets
    functionValue = m.functionValue
    violation = m.violation
    functions = m.functions
    sets = m.sets

    for i in 1:length(functions)
        violation[i] = MOD.distance_to_set(MOD.DefaultDistance(), functionValue[i], sets[i])
    end
    nothing
end


text:
mutable struct CallbackData
    model_root::XpressProblem # should not use operations here
    data::Any # data for user
    model::XpressProblem # local model # ptr_model::Ptr{Nothing}
end

# must be mutable
mutable struct _CallbackUserData
    callback::Function
    model::XpressProblem
    data::Any
end
Base.cconvert(::Type{Ptr{Cvoid}}, x::_CallbackUserData) = x
function Base.unsafe_convert(::Type{Ptr{Cvoid}}, x::_CallbackUserData)
    return pointer_from_objref(x)::Ptr{Cvoid}
end

export CallbackData

function setcboptnode_wrapper(ptr_model::Lib.XPRSprob, my_object::Ptr{Cvoid}, feas::Ptr{Cint})
    usrdata = unsafe_pointer_to_objref(my_object)::_CallbackUserData
    callback, model, data = usrdata.callback, usrdata.model, usrdata.data
    model_inner = XpressProblem(ptr_model; finalize_env = false)
    callback(CallbackData(model, data, model_inner))
    return zero(Cint)
end

function set_callback_optnode!(model::XpressProblem, callback::Function, data::Any = nothing)
    callback_ptr = @cfunction(setcboptnode_wrapper, Cint, (Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cint}))
    usrdata = _CallbackUserData(callback, model, data)
    Lib.XPRSaddcboptnode(model.ptr, callback_ptr, usrdata, 0)
    # we need to keep a reference to the callback function
    # so that it isn't garbage collected
    return callback_ptr, usrdata
end

function setcbpreintsol_wrapper(ptr_model::Lib.XPRSprob, my_object::Ptr{Cvoid}, soltype::Ptr{Cint}, ifreject::Ptr{Cint}, cutoff::Ptr{Cdouble})
    usrdata = unsafe_pointer_to_objref(my_object)::_CallbackUserData
    callback, model, data = usrdata.callback, usrdata.model, usrdata.data
    model_inner = XpressProblem(ptr_model; finalize_env = false)
    callback(CallbackData(model, data, model_inner))
    return zero(Cint)
end

function set_callback_preintsol!(model::XpressProblem, callback::Function, data::Any = nothing)
    callback_ptr = @cfunction(setcbpreintsol_wrapper, Cint, (Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cint}, Ptr{Cint}, Ptr{Cdouble}))
    usrdata = _CallbackUserData(callback, model, data)
    Lib.XPRSaddcbpreintsol(model.ptr, callback_ptr, usrdata,0)
    # we need to keep a reference to the callback function
    # so that it isn't garbage collected
    return callback_ptr, usrdata
end

function addcboutput2screen_wrapper(ptr_model::Lib.XPRSprob, my_object::Ptr{Cvoid}, msg::Ptr{Cchar}, len::Cint, msgtype::Cint)

    usrdata = unsafe_pointer_to_objref(my_object)::_CallbackUserData
    show_warning = usrdata.data::Bool

    if msgtype == 4
        #= Error =#
        return zero(Cint)
    elseif msgtype == 3 && show_warning
        #= Warning =#
        msg_str = unsafe_string(msg)
        println(msg_str)
        return zero(Cint)
    elseif msgtype == 2
        #= Not used =#
        return zero(Cint)
    elseif msgtype == 1
         #= Information =#
         msg_str = unsafe_string(msg)
         println(msg_str)
         return zero(Cint)
    else
        #= Exiting - buffers need flushing  =#
        flush(stdout)
        return zero(Cint)
    end
end

function setoutputcb!(model::XpressProblem, show_warning::Bool)
    callback_ptr = @cfunction(addcboutput2screen_wrapper, Cint, (Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cchar}, Cint, Cint))
    usrdata = _CallbackUserData(_null, model, show_warning)
    Lib.XPRSaddcbmessage(model.ptr, callback_ptr, usrdata, 0)
    return callback_ptr, usrdata
end

function _null()
    return
end


text:
function _update_initial_conditions!(
    model::OperationModel,
    key::ICKey{T, U},
    source, # Store or State are used in simulations by default
) where {T <: InitialConditionType, U <: PSY.Component}
    if get_execution_count(model) < 1
        return
    end
    container = get_optimization_container(model)
    model_resolution = get_resolution(model.internal.store_parameters)
    ini_conditions_vector = get_initial_condition(container, key)
    timestamp = get_current_timestamp(model)
    previous_values = get_condition.(ini_conditions_vector)
    # The implementation of specific update_initial_conditions! is located in the files
    # update_initial_conditions_in_memory_store.jl and update_initial_conditions_simulation.jl
    update_initial_conditions!(ini_conditions_vector, source, model_resolution)
    for (i, initial_condition) in enumerate(ini_conditions_vector)
        IS.@record :execution InitialConditionUpdateEvent(
            timestamp,
            initial_condition,
            previous_values[i],
            get_name(model),
        )
    end
    return
end

# Note to devs: Implemented this way to avoid ambiguities and future proof custom ic updating
function update_initial_conditions!(
    model::DecisionModel,
    key::ICKey{T, U},
    source, # Store or State are used in simulations by default
) where {T <: InitialConditionType, U <: PSY.Component}
    _update_initial_conditions!(model, key, source)
    return
end
function update_initial_conditions!(
    model::EmulationModel,
    key::ICKey{T, U},
    source, # Store or State are used in simulations by default
) where {T <: InitialConditionType, U <: PSY.Component}
    _update_initial_conditions!(model, key, source)
    return
end


text:
"""
Module containing routines to solve and perturb nonlinear eigenvalue problems.
"""
module NLEVP
include("NLEVP_exports.jl")
#header
using SparseArrays

include("./NLEVP/algebra.jl")
include("./NLEVP/polys_pade.jl")
include("./NLEVP/LinOpFam.jl")
include("./NLEVP/perturbation.jl")
include("./NLEVP/Householder.jl")
include("./NLEVP/nicoud.jl")
include("./NLEVP/picard.jl")
include("./NLEVP/iterative_solvers.jl")
#include("./NLEVP/mehrmann.jl")
include("./NLEVP/beyn.jl")
include("./NLEVP/solver.jl")
include("./NLEVP/toml.jl")
include("./NLEVP/save.jl")
include("./NLEVP/gallery.jl")
#using .TOML
end


text:
rx = r"[\e\[[0-9;]*m]*"

t = "┌──────────\e[31m────\e[39m──┬\e[1m\e[34m────────────────┬\e[22m\e[39m\e[39m──────────────┬────────────────┬──────────────end"

collect(eachmatch(rx, t))


text:
module PHonon

using AbInitioSoftwareBase.Inputs: Namelist, groupname
using Compat: only
using Crystallography: ReciprocalPoint
using PyFortran90Namelists: fparse
using QuantumESPRESSOBase.Inputs.PHonon:
    ReciprocalPoint,
    QPointsCard,
    PhInput,
    Q2rInput,
    DynmatInput,
    MatdynInput,
    PhNamelist,
    Q2rNamelist,
    MatdynNamelist,
    DynmatNamelist

const Q_POINTS_SPECIAL_BLOCK_REGEX = r"""
^ [ \t]* qPointsSpecs [ \t]*$\R
^ [ \t]* \S+ [ \t]* $\R  # nqs
(?P<block>
 (?:
  ^ [ \t]* \S+ [ \t]+ \S+ [ \t]+ \S+ [ \t]+ \S+ [ \t]* $\R?
 )+
)
"""imx
const Q_POINTS_SPECIAL_ITEM_REGEX = r"""
^ [ \t]*
([\-|\+]?(?:\d*[\.]\d+ | \d+[\.]?\d*)
    (?:[E|e|d|D][+|-]?\d+)?)\h*
([\-|\+]?(?:\d*[\.]\d+ | \d+[\.]?\d*)
    (?:[E|e|d|D][+|-]?\d+)?)\h*
([\-|\+]?(?:\d*[\.]\d+ | \d+[\.]?\d*)
    (?:[E|e|d|D][+|-]?\d+)?)\h*
([\-|\+]?(?:\d*[\.]\d+ | \d+[\.]?\d*)
    (?:[E|e|d|D][+|-]?\d+)?)\h*
"""imx

function Base.tryparse(::Type{QPointsCard}, str::AbstractString)
    m = match(Q_POINTS_SPECIAL_BLOCK_REGEX, str)
    if m === nothing
        @info "no q-points are provided."
        return
    else
        captured = only(m.captures)
        data = map(eachmatch(Q_POINTS_SPECIAL_ITEM_REGEX, captured)) do matched
            # TODO: Match `nqs`
            ReciprocalPoint(map(x -> fparse(Float64, x), matched.captures))
        end
        return QPointsCard(data)
    end
end # function Base.tryparse
function Base.tryparse(::Type{PhInput}, str::AbstractString)
    title_line = ""
    for line in split(str, r"\r\n|\r|\n")
        if !isempty(line)
            title_line = line
            break
        end
    end
    args = []
    for T in (PhNamelist, QPointsCard)
        push!(args, tryparse(T, str))
    end
    if all(x === nothing for x in args)
        return
    else
        return PhInput(title_line, args[1], args[2])
    end
end # function Base.tryparse
function Base.tryparse(::Type{Q2rInput}, str::AbstractString)
    x = parse(Q2rNamelist, str)
    if x !== nothing
        return Q2rInput(parse(Q2rNamelist, str))
    else
        return
    end
end # function Base.tryparse
function Base.parse(::Type{MatdynInput}, str::AbstractString)
    args = []
    for T in (MatdynNamelist, QPointsCard)
        push!(args, tryparse(T, str))
    end
    if all(x === nothing for x in args)
        return
    else
        return MatdynInput(args[1], args[2])
    end
end # function Base.tryparse
function Base.tryparse(::Type{DynmatInput}, str::AbstractString)
    if parse(DynmatNamelist, str)
        return DynmatInput(parse(DynmatNamelist, str))
    else
        return
    end
end # function Base.tryparse

function Base.parse(
    ::Type{T},
    str::AbstractString,
) where {T<:Union{QPointsCard,PhInput,Q2rInput,DynmatInput,MatdynInput}}
    x = tryparse(T, str)
    if x === nothing
        throw(Meta.ParseError("cannot find card `$(groupname(T))`!"))
    else
        return x
    end
end # function Base.parse

end


text:
module athing
include("another_thing.jl")
# Write your package code here.

end

max_stars_repo_name: Sarah-W/athing

max_issues_repo_name: Sarah-W/athing

max_forks_repo_name: Sarah-W/athing


text:
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.11.2
#   kernelspec:
#     display_name: Julia 1.7.0-DEV
#     language: julia
#     name: julia-1.7
# ---

# %%
using DataFrames, CSV
csvstring(df) = String(take!(CSV.write(IOBuffer(), df)))
csv2df(str) = CSV.read(IOBuffer(str), DataFrame)

# %%
using Random
Random.seed!(4649373)

# %%
data = round.(10randn(7, 5), digits=2)
names = string.('a':'e')
df = DataFrame(data, names)

# %%
str = csvstring(df)

# %%
df1 = csv2df(str)

# %%
df1 == df

# %%
print(str)

# %%
# copy and paste the above

df2 = csv2df("""
a,b,c,d,e
-5.91,19.26,19.44,1.03,1.85
-3.67,-10.33,1.57,-17.67,-0.47
6.72,8.81,-3.87,-11.0,-0.56
-6.2,10.66,11.52,6.19,-7.78
-8.64,-1.6,9.28,-1.8,3.72
-1.01,-17.56,7.02,1.08,18.56
-20.48,-9.66,5.31,-17.65,11.51
""")

# %%
df2 == df

# %%

max_stars_repo_path: 0003/dataframe to string.jl

max_stars_repo_name: genkuroki/public

max_issues_repo_path: 0003/dataframe to string.jl

max_issues_repo_name: genkuroki/public

max_forks_repo_path: 0003/dataframe to string.jl

max_forks_repo_name: genkuroki/public


text:
using Documenter
using Pkg
using Test
using SymbolicUtils
import IfElse: ifelse

DocMeta.setdocmeta!(
    SymbolicUtils,
    :DocTestSetup,
    :(using SymbolicUtils);
    recursive=true
)

# Only test one Julia version to avoid differences due to changes in printing.
if v"1.6" ≤ VERSION < v"1.7-beta3.0"
    doctest(SymbolicUtils)
else
    @warn "Skipping doctests"
end

# == / != syntax is nice, let's use it in tests
macro eqtest(expr)
    @assert expr.head == :call && expr.args[1] in [:(==), :(!=)]
    if expr.args[1] == :(==)
        :(@test isequal($(expr.args[2]), $(expr.args[3])))
    else
        :(@test !isequal($(expr.args[2]), $(expr.args[3])))
    end |> esc
end
SymbolicUtils.show_simplified[] = false

include("basics.jl")
include("order.jl")
include("rewrite.jl")
include("rulesets.jl")
include("code.jl")
include("nf.jl")
include("interface.jl")
include("fuzz.jl")
include("adjoints.jl")


text:
using DelimitedFiles

filename = "./log-conv-euler-2D-quad-beta8.5-wopost-CFL.9"
file = readdlm(filename)

num_scheme = 1
Narr = [1;2;3;4]
Karr = [2^2*2;4^2*2;8^2*2;16^2*2;32^2*2]
block_size = 4
num_lines = size(file,1)
num_block = div(num_lines,block_size)
L1_err_dict = Dict()
L2_err_dict = Dict()
for b in 1:num_block
    l = 1+block_size*(b-1)
    N = parse(Int,file[l,3][1])
    K = file[l,6]
    L1_err = file[l+1,4]
    L2_err = file[l+2,4]
    if !haskey(L1_err_dict,(N,K))
        L1_err_dict[(N,K)] = [L1_err]
    else
        push!(L1_err_dict[(N,K)],L1_err)
    end

    if !haskey(L2_err_dict,(N,K))
        L2_err_dict[(N,K)] = [L2_err]
    else
        push!(L2_err_dict[(N,K)],L2_err)
    end
end

num_row = length(Karr)-1
num_col = length(Narr)
conv_table_L1 = zeros(num_row,num_col)
conv_table_L2 = zeros(num_row,num_col)
for i = 1:num_row
    for j = 1:num_col
        K = Karr[i]
        N = Narr[j]
        conv_table_L1[i,j] = log2(L1_err_dict[(N,K)][1]/L1_err_dict[(N,4*K)][1])
        conv_table_L2[i,j] = log2(L2_err_dict[(N,K)][1]/L2_err_dict[(N,4*K)][1])
    end
end

err_table = zeros(length(Karr),2*length(Narr))
for i = 1:length(Karr)
    for j = 1:2*length(Narr)
        K = Karr[i]
        N = Narr[div(j-1,2)+1]
        is_L1 = mod1(j,2) == 1
        if is_L1
            err_table[i,j] = L1_err_dict[(N,K)][1]
        else
            err_table[i,j] = L2_err_dict[(N,K)][1]
        end
    end
end

conv_table = [conv_table_L1[:,1] conv_table_L2[:,1] conv_table_L1[:,2] conv_table_L2[:,2] conv_table_L1[:,3] conv_table_L2[:,3] conv_table_L1[:,4] conv_table_L2[:,4]]

max_stars_repo_name: yiminllin/ESDG-PosLimit

max_issues_repo_name: yiminllin/ESDG-PosLimit

max_forks_repo_name: yiminllin/ESDG-PosLimit


text:
"""
Code declaring abstract types for instruments' trait system.

Author: Eric Ford
Created: August 2020
"""

""" Abstract Base type for Instruments """
abstract type AbstractInstrument end

""" Abstract Base type for 1-D extracted spectra
    Should be specialized for specific instruments. """
abstract type AbstractInstrument1D <: AbstractInstrument end

""" Abstract Base type for 2-D extracted spectra.
    Should be specialized for each instrument """
abstract type AbstractInstrument2D <: AbstractInstrument end

#=
# TODO: Decide if these are really a good idea or if we shouldn't allow them.
""" Trait for genric 1D Extracted spectra """
struct Generic1D <: AbstractInstrument1D end

""" Trait for generic 2D Extracted spectra  """
struct Generic2D <: AbstractInstrument2D end
=#

max_stars_repo_name: RvSpectML/RvSpectMLTemplate

max_issues_repo_name: RvSpectML/RvSpectMLTemplate

max_forks_repo_name: RvSpectML/RvSpectMLTemplate


text:
facts("VectorQFunction with in-place bellman!") do
    qfunc = VectorQFunction(fixture_num_states())
    trans = TransitionProbability(fixture_transition())
    reward = Reward(fixture_reward())
    bellman!(qfunc, fixture_initial_value(), trans, reward, fixture_discount())
    @fact value(qfunc) => roughly(fixture_expected_value())
    @fact policy(qfunc) == fixture_expected_policy() => true
end

facts("VectorQFunction with bellman") do
    qfunc = VectorQFunction(fixture_num_states())
    trans = TransitionProbability(fixture_transition())
    reward = Reward(fixture_reward())
    newq = bellman(qfunc, fixture_initial_value(), trans, reward, fixture_discount())
    @fact value(newq) => roughly(fixture_expected_value())
    @fact policy(newq) == fixture_expected_policy() => true
    @fact all(map((x, y) -> isapprox(x, y), value(qfunc), fixture_expected_value())) => false
    @fact policy(qfunc) == fixture_expected_policy() => false
end

facts("ArrayQFunction with bellman!") do
    qfunc = ArrayQFunction(fixture_num_states(), fixture_num_actions())
    trans = TransitionProbability(fixture_transition())
    reward = Reward(fixture_reward())
    bellman!(qfunc, fixture_initial_value(), trans, reward, fixture_discount())
    @fact value(qfunc) => roughly(fixture_expected_value())
    @fact policy(qfunc) == fixture_expected_policy() => true
end

facts("SparseArrayTransitionProbability with bellman!") do
    qfunc = ArrayQFunction(fixture_num_states(), fixture_num_actions())
    T = SparseMatrixCSC{Float64,Int}
    array = fixture_transition()
    array = T[ sparse(array[:, :, a]) for a=1:fixture_num_actions() ]
    trans = TransitionProbability(array)
    reward = Reward(fixture_reward())
    bellman!(qfunc, fixture_initial_value(), trans, reward, fixture_discount())
    @fact value(qfunc) => roughly(fixture_expected_value())
    @fact policy(qfunc) == fixture_expected_policy() => true
end

facts("FunctionTransitionProbability with bellman!") do
    qfunc = ArrayQFunction(fixture_num_states(), fixture_num_actions())
    array = fixture_transition()
    trans = TransitionProbability((s, t, a) -> array[s, t, a], fixture_num_states(), fixture_num_actions())
    reward = Reward(fixture_reward())
    bellman!(qfunc, fixture_initial_value(), trans, reward, fixture_discount())
    @fact value(qfunc) => roughly(fixture_expected_value())
    @fact policy(qfunc) == fixture_expected_policy() => true
end

facts("ArrayReward with bellman!") do
    qfunc = VectorQFunction(fixture_num_states())
    trans = TransitionProbability(fixture_transition())
    reward = Reward(fixture_reward(:array))
    bellman!(qfunc, fixture_initial_value(), trans, reward, fixture_discount())
    @fact value(qfunc) => roughly(fixture_expected_value(:array))
    @fact policy(qfunc) == fixture_expected_policy(:array) => true
end

facts("bellman can be called with arrays") do
    Q = zeros(5, 3)
    newQ = bellman(Q, fixture_initial_value(), fixture_transition(), fixture_reward(), fixture_discount())
    @fact typeof(newQ) => Array{Float64,2}
    @fact Q == newQ => false
    @fact Float64[ maximum(newQ[s, :]) for s=1:5 ] => roughly(fixture_expected_value())
    @fact Int[ indmax(newQ[s, :]) for s=1:5 ] == fixture_expected_policy() => true
end

facts("bellman can be called with vectors") do
    value = zeros(5)
    policy = zeros(Int, 5)
    newVal, newPol = bellman(value, policy, fixture_initial_value(), fixture_transition(), fixture_reward(), fixture_discount())
    @fact typeof(newVal) => Vector{Float64}
    @fact typeof(newPol) => Vector{Int}
    @fact value == newVal => false
    @fact policy == newPol => false
    @fact newVal => roughly(fixture_expected_value())
    @fact newPol == fixture_expected_policy() => true
end

facts("bellman! can be called with arrays") do
    Q = zeros(5, 3)
    newQ = bellman!(Q, fixture_initial_value(), fixture_transition(), fixture_reward(), fixture_discount())
    @fact Q === newQ => true
    @fact Float64[ maximum(newQ[s, :]) for s=1:5 ] => roughly(fixture_expected_value())
    @fact Int[ indmax(newQ[s, :]) for s=1:5 ] == fixture_expected_policy() => true
end

facts("bellman! can be called with vectors") do
    value = zeros(5)
    policy = zeros(Int, 5)
    newVal, newPol = bellman!(value, policy, fixture_initial_value(), fixture_transition(), fixture_reward(), fixture_discount())
    @fact value === newVal => true
    @fact policy === newPol => true
    @fact value => roughly(fixture_expected_value())
    @fact policy == fixture_expected_policy() => true
end


text:
using RealAlg, CSV, DataFrames

# Cell_Df = CSV.read("/Users/bradyplanden/Documents/Git/RealisationAlgorithms/RealAlg/Data/Sample-Cell.csv", DataFrame)
# Constants_Df = tuple.(Tuple.(eachrow(dropmissing(Cell_Df[:,1:2])))...)
# Geometry_Df = tuple.(Tuple.(eachrow(dropmissing(Cell_Df[:,3:4])))...)
# Negative_Df = tuple.(Tuple.(eachrow(dropmissing(Cell_Df[:,5:6])))...)
# Positive_Df = tuple.(Tuple.(eachrow(dropmissing(Cell_Df[:,7:8])))...)
# Seperator_Df = tuple.(Tuple.(eachrow(dropmissing(Cell_Df[:,9:10])))...)

# CellData = Cell(Constants(1,1.2,1.3,1.4,1.5,1.6),Geometry(1,1.2,1.3,1.4,1.5),Negative(1,1.2,1.3,1.4,1.5,1.6,1.7,1.8,1.9,1.10,11.11,1.12,1.13,1.14),Positive(1,1.2,1.3,1.4,1.5,1.6,1.7,1.8,1.9,1.10,11.11,1.12,1.13,1.14),Seperator(1,1.2))
CellData = Cell(Constants(),Geometry(),Negative(),Positive(),Seperator())

M=4
Fs = 2
Tlen = 32768
Nfft = 2^(ceil(log2(Fs*Tlen)))
f = 0:Nfft-1
s = (2im.*Fs)*tan.(pi.*f./Nfft)

# loc = Any[0,1]
loc = Any[0, 128e-6, 204e-6, 394e-6]
Def = "Neg"

test_ce = C_e(CellData,s,M)
println("Ce:",length(test_ce))
test_j = j(CellData,s,loc,Def)
println("j:",typeof(test_j))
test_cse = C_se(CellData,s,loc,Def)
println("Cse:",typeof(test_cse))
test_ϕ = Phi_s(CellData,s,loc,Def)
println("ϕ_s:",typeof(test_ϕ))
test_ϕ_e = Phi_e(CellData,s,loc,Def)
println("ϕ_e:", length(test_ϕ_e))

# c = Cell_Df[:,2]
# c = dropmissing(Cell_Df)
# x = copy(Cell_Df)
# y = Tuple(x)

# c = dropmissing(Cell_Df[:,1:2])
# I1 = [x[1] for x in j]

# q1 = Tuple.(eachrow(dropmissing(Cell_Df[:,1:2])))
# q1 = tuple.(q1...)

max_stars_repo_name: BradyPlanden/RealisationAlgorithms

max_issues_repo_name: BradyPlanden/RealisationAlgorithms

max_forks_repo_name: BradyPlanden/RealisationAlgorithms


text:
module Waterflow

using LightGraphs, MetaGraphs, Distributed, DataFrames
using Infiltrator

# Can't use string, DLL location has to be a const
# (which makes sense but still, many hours wasted!)
# https://github.com/JuliaLang/julia/issues/29602
const ihacres = "../../ihacres_nim/lib/ihacres.dll"

"""@def macro

Inline code to avoid repetitious declarations.
"""
macro def(name, definition)
    return quote
        macro $(esc(name))()
            esc($(Expr(:quote, definition)))
        end
    end
end

@def add_preprefix begin
    if !isnothing(id_prefix)
        prefix = id_prefix * prefix
    end
end

include("Node.jl")
include("StreamNode.jl")
include("DamNode.jl")
include("Climate.jl")

"""Run a model attached to a node.

water_order: Volume of water to be extracted (in ML/timestep)
exchange: Volume of flux (in ML/timestep), where negative values are losses to the groundwater system.
"""
function run_node!(mg::MetaGraph, g::AbstractGraph, node_id::Int, climate::Climate, timestep::Int;
                   water_order=nothing, exchange=nothing)

    curr_node = get_prop(mg, node_id, :node)
    if checkbounds(Bool, curr_node.outflow, timestep)
        # already ran for this time step so no need to recurse further
        return curr_node.outflow[timestep], curr_node.level[timestep]
    end

    outflow = 0.0
    inflow = 0.0

    ins = inneighbors(g, node_id)
    if length(ins) == 0
        inflow = 0.0
        src_name = "Out-of-Catchment"
    else
        inflow = 0.0
        for i in ins
            src_name = get_prop(mg, i, :name)
            # Get inflow from previous node
            upstream_flow, upstream_level = run_node!(mg, g, i, climate, timestep)
            inflow += upstream_flow
        end
    end

    gauge_id = get_prop(mg, node_id, :name)
    rain, et = climate_values(curr_node, climate, timestep)
    if ismissing(rain) | ismissing(et)
        @warn "No climate data found for node $(node_id) ($(gauge_id)) using column markers $(climate.rainfall_id) and $(climate.et_id)"
        return 0.0
    end

    wo = 0.0
    if !isnothing(water_order)
        release_col = filter(x -> occursin(gauge_id, string(x))
                                  & occursin("releases", string(x)),
                                  names(water_order))

        wo = checkbounds(Bool, water_order.Date, timestep) ? water_order[timestep, release_col][1] : 0.0
    end

    ex = 0.0
    if !isnothing(exchange)
        exchange_col = filter(x -> occursin(gauge_id, x)
                                  & occursin("exchange", x),
                                  names(exchange))
        ex = checkbounds(Bool, exchange.Date, timestep) ? exchange[timestep, exchange_col][1] : 0.0
    end

    # Calculate outflow for this node
    func = get_prop(mg, node_id, :nfunc)
    if curr_node isa StreamNode
        outflow, level = func(curr_node, rain, et, inflow, 0.0)
    elseif curr_node isa DamNode
        outflow, level = func(curr_node, rain, et, inflow, wo, ex)
    else
        throw(ArgumentError("Unknown node type!"))
    end

    return outflow, level
end

include("Network.jl")

export @def
export ihacres, StreamNode, DamNode, Climate
export find_inlets_and_outlets, create_network, create_node
export update_params!, run_node!, reset!, sim_length
export climate_values

end  # end module

max_stars_repo_name: ConnectedSystems/camgraph

max_issues_repo_name: ConnectedSystems/camgraph

max_forks_repo_name: ConnectedSystems/camgraph


text:
for _dummy_ in 1:1 # avoid global variable warnings
    # reseeding with random seed
    rng = LazySets.GLOBAL_RNG
    seed = rand(1:10000)
    LazySets.reseed(rng, seed)
    n1 = rand(Int)
    LazySets.reseed(rng, seed)
    n2 = rand(Int)
    @test n1 == n2

    # StrictlyIncreasingIndices
    vectors = Vector{AbstractVector{Int}}()
    for v in LazySets.StrictlyIncreasingIndices(5, 4)
        push!(vectors, copy(v))
    end
    @test vectors == [[1, 2, 3, 4], [1, 2, 3, 5], [1, 2, 4, 5], [1, 3, 4, 5], [2, 3, 4, 5]]

    # square matrix
    @test LazySets.issquare([2 3; 0 0])
    @test LazySets.issquare(sparse([1], [1], [1], 3, 3))
    @test !LazySets.issquare(hcat([2 3]))
    @test !LazySets.issquare(sparse([1], [1], [1], 3, 4))

    # invertible matrix: (1) invertible, (2) non-invertible, (3) non-square
    # dense matrix
    @test LazySets.isinvertible([2 3; 1 2])
    @test !LazySets.isinvertible([2 3; 0 0])
    @test !LazySets.isinvertible(hcat([2 3]))
    # sparse matrix
    @test LazySets.isinvertible(sparse([2 3; 1 2]))
    @test !LazySets.isinvertible(sparse([2 3; 0 0]))
    @test !LazySets.isinvertible(sparse(hcat([2 3])))
    # diagonal matrix
    @test LazySets.isinvertible(Diagonal([2 0; 0 2]))
    @test !LazySets.isinvertible(Diagonal([2 0; 0 0]))

    # matrix rank
    A = sprandn(2, 10, 0.4)
    @test rank(Matrix(A)) == rank(A) == rank(view(A, :, :))

    for N in [Float64, Rational{Int}, Float32]
        # substitution
        x = N[1, 2, 3]
        substitution = Dict(1 => N(4), 3 => N(0))
        @test LazySets.substitute(substitution, x) == N[4, 2, 0]
        LazySets.substitute!(substitution, x)
        @test x == N[4, 2, 0]

        A = N[1 4; 2 5; 3 6]
        x1 = N[0, 2, 0]
        y1 = N[3, 0]
        x2 = SingleEntryVector(2, 3, N(2))
        y2 = SingleEntryVector(1, 2, N(3))
        @test -x2 == SingleEntryVector(2, 3, N(-2))
        @test inner(x1, A, y1) == dot(x1, A * y1) == inner(x2, A, y2) ==
              dot(x2, A * y2) == N(12)

        x = N[0, 1, -1]
        @test LazySets.rectify(x) == N[0, 1, 0]

        # approximate permutation check
        v1 = [N[1, 2], N[2, 3], N[3, 4]]
        tol = Base.rtoldefault(N)/2
        v2 = [N[3 + tol, 4 - tol], N[2 - tol, 3 + tol], N[1 + tol, 2 + tol]]
        v3 = [N[2, 3], N[1, 2], N[4, 4]]
        @test ispermutation(v1, v2) && !ispermutation(v2, v3)

        # to_negative_vector
        for v in [N[-1, 0, 1], sparsevec([1, 3], N[-1, 1], 3), N(-1):N(1):N(1)]
            u = LazySets.Arrays.to_negative_vector(v)
            @test u isa Vector{N} && u == N[1, 0, -1]
        end
    end

    for N in [Float64, Float32]
        # modified dot product
        @test isnan(dot(N[1, 0], N[Inf, -Inf]))
        @test LazySets.dot_zero(N[1, 0], N[Inf, -Inf]) == N(Inf)

        # remove zero columns
        m = 3
        n = 10
        A = rand(N, m, n)
        A[:, 1] = A[:, 5] = A[:, n] = zeros(N, m)
        B = LazySets.delete_zero_columns!(A)
        @test size(B) == (m, n-3)
    end
end


text:
#==
This example shows multiobjective Paretofronts with 3 objectives.
The front is convex. The NBI method does not catch the full Pareto front,
e.g. as seen by the missing arcs between f1 and f2, whilst f3 is 0.
To capture this area, we need some of the elements in β to be negative.

Taken from
Il Yong Kim and OL De Weck.
Adaptive weighted sum method for multiobjective optimization:
a new method for pareto front generation.
Structural and Multidisciplinary Optimization, 31(2):105–116, 2006.
==#

using MultiJuMP, JuMP
using Ipopt

m = multi_model(Ipopt.Optimizer)

@variable(m, x[i=1:3] >= 0)
@NLexpression(m, f1, x[1])
@NLexpression(m, f2, x[2])
@NLexpression(m, f3, x[3])
@NLconstraint(m, x[1]^4+2x[2]^3+5x[3]^2<=1)

obj1 = SingleObjective(f1, sense = MOI.MAX_SENSE)#,
                       #iv = Dict{String, Any}("x[$i]" => j for (i,j) in enumerate([0, 0., 1.2])))

obj2 = SingleObjective(f2, sense = MOI.MAX_SENSE)#,
                       #iv = Dict{String, Any}("x[$i]" => j for (i,j) in enumerate([0, 0., 1.2])))

obj3 = SingleObjective(f3, sense = MOI.MAX_SENSE)#,
                       #iv = Dict{String, Any}("x[$i]" => j for (i,j) in enumerate([0, 0., 2.])))

md = get_multidata(m)
md.objectives = [obj1, obj2, obj3]
#md.objectives = [obj1, obj2]
md.pointsperdim = 10
optimize!(m, method = NBI())

using Plots
pyplot()
plot(md)


text:
## Julia program for Random walk Metropolis-Hastings for Normal distribution
## author: weiya <szcfweiya@gmail.com>
## date: 2018-08-21

# random walk Metropolis-Hastings
function rmh(T, delta, f::Function, initval = 5)
    x = ones(T+1)
    x[1] = initval
    for t = 1:T
        # generate Yt
        epsi = rand() * 2 * delta - delta
        Yt = epsi + x[t]
        # accept or not
        u = rand()
        r = f(Yt)/f(x[t])
        if r >= 1
            x[t+1] = Yt
        else
            if u <= r
                x[t+1] = Yt
            else
                x[t+1] = x[t]
            end
        end
    end
    return(x)
end

# density function of N(0, 1) without normalization
function dnorm(x)
    return(exp(-0.5 * x^2))
end

# example
deltalist = [0.1, 0.5, 1]
N = 15000
# results of mean and variance
muvar = zeros(3, 2)
using Statistics
ps = []
for i = 1:3
    x = rmh(N, deltalist[i], dnorm)
    push!(ps, plot(x, legend=false))
    muvar[i, 1] = mean(x)
    muvar[i, 2] = var(x)
end
# plot
using Plots
plot(ps[1], ps[2], ps[3], layout=(3,1))
savefig("res_rmh.png")

max_stars_repo_name: szcf-weiya/MCMC

max_issues_repo_name: szcf-weiya/MCMC

max_forks_repo_name: szcf-weiya/MCMC


text:
using SymbolicUtils
using Test

using SymbolicUtils: showraw, Symbolic

function rand_input(T)
    if T == Bool
        return rand(Bool)
    elseif T <: Integer
        x = rand(-100:100)
        while iszero(x)
            x = rand(-100:100)
        end
        return x
    elseif T == Rational
        return Rational(rand_input(Int), rand(1:50)) # no 0 denominator tests yet!
    elseif T == Real
        return rand_input(rand([Int, Float64, Rational]))
    elseif T == Complex
        return rand_input(Real) + rand_input(Real) * im
    elseif T == Number
        # more real than complex
        return rand_input(rand([Real, Real, Real, Complex{Float64}]))
    elseif !isabstracttype(T)
        return rand(T)
    else
        error("Don't know how to generate input of type $T")
    end
end

rand_input(i::Symbolic{T}) where {T} = rand_input(T)

const num_spec = let
    @syms a b::Real c::Integer d::Float64 e::Rational f

    leaf_funcs = [()->rand_input(Real),
                  ()->rand_input(Complex),
                  ()->rand([a b c d e f]),
                  ()->rand([a b c d e f])]

    binops = SymbolicUtils.diadic
    nopow  = filter(x->x!==(^), binops)
    twoargfns = vcat(nopow, (x,y)->x isa Union{Int, Rational, Complex{<:Rational}} ? x * y : x^y)
    fns = vcat(1 .=> vcat(SymbolicUtils.monadic, [one, zero]),
               2 .=> vcat(twoargfns, fill(+, 5), [-,-], fill(*, 5)),
    3 .=> [+, *])

    (leaves=leaf_funcs, funcs=fns, input=rand_input)
end

const bool_spec = let
    @syms b::Bool x::Real y::Real z::Complex

    bool_leaf_funcs = [()->rand(Bool),
                       ()->rand([b, (x, b) => ((x > 0) | b), (x,)=>(x < 0), (y,z) => y==z, (y, z) => y!=z])]

    fns = vcat(1 .=> [(!), (~)],
               2 .=> [(|), (&), xor],
               3 .=> [cond]) # cond will still stay in bool by condtruction

    (leaves=bool_leaf_funcs,
     funcs=fns,
     input=rand_input
    )
end

function gen_rand_expr(inputs;
                       spec=num_spec,
                       leaf_prob=0.92,
                       depth=0,
                       min_depth=1,
                       max_depth=5)

    if depth > max_depth  || (min_depth <= depth && rand() < leaf_prob)
        leaf = rand(spec.leaves)()
        if leaf isa SymbolicUtils.Sym
            push!(inputs, leaf)
        elseif leaf isa Pair
            foreach(i->push!(inputs, i), leaf[1])
            return leaf[2]
        end
        return leaf
    end
    arity, f = rand(spec.funcs)
    args = [gen_rand_expr(inputs,
                          spec=spec,
                          leaf_prob=leaf_prob,
                          depth=depth+1,
                          min_depth=min_depth,
                          max_depth=max_depth) for i in 1:arity]
    try
        return f(args...)
    catch err
        if err isa DomainError || err isa DivideError || err isa MethodError
            return gen_rand_expr(inputs,
                                 spec=spec,
                                 leaf_prob=leaf_prob,
                                 depth=depth,
                                 min_depth=min_depth,
                                 max_depth=max_depth)
        else
            rethrow(err)
        end
    end
end

struct Errored
    err
end

function fuzz_test(ntrials, spec, simplify=simplify;kwargs...)
    inputs = Set()
    expr = gen_rand_expr(inputs; spec=spec, kwargs...)
    inputs = collect(inputs)
    unsimplifiedstr = """
    function $(tuple(inputs...))
        $(sprint(io->showraw(io, expr)))
    end
    """

    simplifiedstr = """
    function $(tuple(inputs...))
        $(sprint(io->showraw(io, simplify(expr))))
    end
    """
    f = include_string(Main, unsimplifiedstr)
    g = include_string(Main, simplifiedstr)

    for i=1:ntrials
        args = [spec.input(i) for i in inputs]
        unsimplified = try
            Base.invokelatest(f, args...)
        catch err
            Errored(err)
        end

        simplified  = try
            Base.invokelatest(g, args...)
        catch err
            Errored(err)
        end
        if unsimplified isa Errored
            if !(simplified isa Errored)
                @test_skip false
                @goto print_err
            end
            @test true
        elseif isnan(unsimplified)
            if !isnan(simplified)
                @test_skip false
                @goto print_err
            end
            @test true
        else
            if !(unsimplified ≈ simplified)
                @test_skip false
                @goto print_err
            end
            @test true
        end
        continue

        @label print_err
        println("""Test failed for expression
                    $(sprint(io->showraw(io, expr))) = $unsimplified
                Simplified:
                    $(sprint(io->showraw(io, simplify(expr)))) = $simplified
                Inputs:
                    $inputs = $args
                """)
    end
end


text:
layout_data(layout::Component) = layout
layout_data(layout::Function) = layout()
function process_layout(request::HTTP.Request, state::HandlerState)
    return HTTP.Response(
        200,
        ["Content-Type" => "application/json"],
        body = JSON3.write(layout_data(state.app.layout))
    )
end


text:
export inf, tropicalpolynomial

import Base: +, *, ==, ^, -

"""
    Tropical(val [, isinf=false])

Create a number in the Tropical semi ring where ``a ⨁ b = max(a,b)`` and ``a ⨀ b = a + b``.
"""
Tropical(val::T, isinf=false) where {T<:Real} = Tropical{T}(val, isinf)
inf(val::T) where {T<:Real} = Tropical(val, true)
inf(::Type{T}) where {T<:Real} = Tropical(zero(T), true)
inf(t::Tropical{T}) where {T<:Real} = Tropical(t.val, true)

Base.isinf(t::Tropical) = t.isinf

function (+)(a::Tropical{T}, b::Tropical{T}) where {T<:Real}
    if isinf(a)
        b
    elseif isinf(b)
        a
    else
        Tropical(max(a.val, b.val))
    end
end

# ⨁(a::Tropical{T}, b::Tropical{T}) where {T<:Real} = a + b
# ⨁(a::Tropical, b::Tropical) = +(promote(a,b)...)
# ⨁(a::T, b::T) where {T<:Real} = Tropical(max(a, b))
# ⨁(a::Real, b::Real) = ⨁(promote(a, b)...)

# this doesn't make sense, but we define it nonetheless to make MultivariatePolynomials
# happy
(-)(a::Tropical) = a

function (*)(a::Tropical{T}, b::Tropical{T}) where {T<:Real}
    if isinf(a)
        a
    elseif isinf(b)
        b
    else
        Tropical(a.val + b.val)
    end
end

# ⨀(a::Tropical{T}, b::Tropical{T}) where {T<:Real} = a * b
# ⨀(a::Tropical, b::Tropical) = *(promote(a,b)...)
# ⨀(a::T, b::T) where {T<:Real} = Tropical(a + b)
# ⨀(a::Real, b::Real) = ⨀(promote(a, b)...)

(^)(a::Tropical, p::Integer) = isinf(a) ? a : Tropical(p * a.val)

function Base.isequal(a::Tropical{T}, b::Tropical{T}) where T
    if isinf(a) || isinf(b)
        isinf(a) == isinf(b)
    else
        a.val == b.val
    end
end
(==)(a::Tropical{T}, b::Tropical{T}) where T = isequal(a,b)

Base.one(::Tropical{T}) where T = Tropical(zero(T))
Base.one(::Type{Tropical{T}}) where T = Tropical(zero(T))

Base.zero(::Tropical{T}) where {T} = Tropical(one(T), true)
Base.zero(::Type{Tropical{T}}) where {T} = Tropical(one(T), true)

Base.show(io::IO, t::Tropical{T}) where {T} = isinf(t) ? print(io, "-∞") : print(io, t.val)

function Base.promote_rule(::Type{Tropical{T}}, ::Type{Tropical{S}}) where {S,T}
    Tropical{promote_type(S,T)}
end
function Base.promote_rule(::Type{T}, ::Type{Tropical{S}}) where {S<:Real,T<:Real}
    Tropical{promote_type(S,T)}
end

Base.convert(::Type{Tropical{T}}, t::Tropical{T}) where {T<:Real} = t
function Base.convert(::Type{Tropical{T}}, t::Tropical{<:Real}) where {T<:Real}
    Tropical{T}(convert(T, t.val), t.isinf)
end
function Base.convert(::Type{Tropical{T}}, x::Real) where {T<:Real}
    if _isone(x)
        one(Tropical{T})
    elseif iszero(x)
        zero(Tropical{T})
    else
        Tropical(convert(T, x))
    end
end

# function Base.convert(::Type{T}, t::Tropical{<:Real}) where {T<:Number}
#     if _isone(t)
#         one(T)
#     elseif iszero(t)
#         zero(T)
#     else
#         convert(T, t.val)
#     end
# end

"""
    tropicalpolynomial(p::MP.AbstractPolynomial)

Interprets a normal polynomials a tropical one. This will especially convert all coefficents which are
1 (neutral element w.r.t. multiplication) to 0 (neutral element wr.t. tropical multiplication).
"""
function tropicalpolynomial(p::MP.AbstractPolynomial)
    mapcoefficients(c -> convert(Tropical{Float64}, c), p)
end

Base.float(t::Tropical) = isinf(t) ? -Inf : float(t.val)


text:
"""
    struct CodeGenOptions <: Any
Hold a pointer to a `clang::CodeGenOptions` object.
"""
struct CodeGenOptions
    ptr::CXCodeGenOptions
end

Base.unsafe_convert(::Type{CXCodeGenOptions}, x::CodeGenOptions) = x.ptr
Base.cconvert(::Type{CXCodeGenOptions}, x::CodeGenOptions) = x


text:
# ==========================================
#    Cell - a single piece of table data
# ==========================================

"""
A `Cell` is a single piece of data annotated by a column name
"""
immutable Cell{Name, ElType}
    data::ElType
    function Cell{T}(x::T)
        check_Cell(Val{Name},ElType)
        new(convert(ElType,x))
    end
end
@compat @inline (::Type{Cell{Name}}){Name, ElType}(x::ElType) = Cell{Name,ElType}(x)

Base.convert{Name, T1, T2}(::Type{Cell{Name, T2}}, x::Cell{Name,T1}) = Cell{Name,T2}(x.data)

@generated function check_Cell{Name, ElType}(::Type{Val{Name}}, ::Type{ElType})
    if !isa(Name, Symbol)
        return :(error("Field name $F should be a symbol"))
    elseif Name == :Row
        return :( error("Field name cannot be :Row") )
    elseif !isa(ElType, DataType)
        return :(error("ElType $ElType should be a data type"))
    else
        return nothing
    end
end

@compat Base.:(==){Name}(cell1::Cell{Name}, cell2::Cell{Name}) = (cell1.data == cell2.data)

@inline rename{Name1, Name2, ElType}(x::Cell{Name1, ElType}, ::Type{Val{Name2}}) = Cell{Name2, ElType}(x.data)

@inline name{Name,ElType}(::Cell{Name,ElType}) = Name
@inline name{Name}(::Type{Cell{Name}}) = Name
@inline name{Name,ElType}(::Type{Cell{Name,ElType}}) = Name
@inline Base.eltype{Name,ElType}(::Cell{Name,ElType}) = ElType
@inline Base.eltype{Name,ElType}(::Type{Cell{Name,ElType}}) = ElType
@inline Base.length{Name,ElType}(::Cell{Name,ElType}) = 1
@inline Base.length{Name,ElType}(::Type{Cell{Name,ElType}}) = 1

@inline nrow(::Cell) = 1
@inline ncol(::Cell) = 1
@inline ncol{C <: Cell}(::Type{C}) = 1

Base.getindex{Name}(c::Cell{Name}, ::Type{Val{Name}}) = c.data
Base.getindex{Name1, Name2}(c::Cell{Name1}, ::Type{Val{Name2}}) = error("Tried to index cell of field name :$Name1 with field name :$Name2")

Base.start(c::Cell) = false # Similar iterators as Julia scalars
Base.next(c::Cell, i::Bool) = (c.data, true)
Base.done(c::Cell, i::Bool) = i
Base.endof(c::Cell) = 1

Base.getindex(c::Cell) = c.data
Base.getindex(c::Cell, i::Integer) = ((i == 1) ? c.data : throw(BoundsError())) # This matches the behaviour of other scalars in Julia
Base.getindex(c::Cell, ::Colon) = c

Base.copy{F,ElType}(cell::Cell{F,ElType}) = Cell{F,ElType}(copy(cell.data))

# @Column and @Cell are very similar
macro Cell(expr)
    if expr.head != :(=) && expr.head != :(kw) # strange Julia bug, see issue 7669
        error("A Expecting expression like @Cell(name::Type = value) or @Cell(name = value)")
    end
    local field
    value = expr.args[2]
    if isa(expr.args[1], Symbol)
        name = expr.args[1]
        return :( TypedTables.Cell{$(QuoteNode(name))}($(esc(value))) )
    elseif isa(expr.args[1],Expr)
        if expr.args[1].head != :(::) || length(expr.args[1].args) != 2 || !isa(expr.args[1].args[1], Symbol)
            error("B Expecting expression like @Cell(name::Type = value) or @Cell(name = value)")
        end
        name = expr.args[1].args[1]
        eltype = expr.args[1].args[2]
        field = :( TypedTables.Cell{$(QuoteNode(name)), $(esc(eltype))}($(esc(value))) )
    else
        error("C Expecting expression like @Cell(name::Type = value) or @Cell(name = value)")
    end
end


text:
## connect do Optim.jl as backend

function sem_fit(model::Sem{O, I, L, D}) where
    {O <: SemObs, L <: Loss, I <: Imply, D <: SemForwardDiff}
    result = Optim.optimize(
                par -> model(par),
                model.imply.start_val,
                model.diff.algorithm,
                autodiff = :forward,
                model.diff.options)
    return result
end

function sem_fit(model::Sem{O, I, L, D}) where
    {O <: SemObs, L <: Loss, I <: Imply, D <: SemFiniteDiff}
    result = Optim.optimize(
                par -> model(par),
                model.imply.start_val,
                model.diff.algorithm,
                model.diff.options)
    return result
end

function sem_fit(model::Sem{O, I, L, D}) where
    {O <: SemObs, L <: Loss, I <: Imply, D <: SemAnalyticDiff}
    if isnothing(model.diff.hessian_functions)
        result = Optim.optimize(
                model,
                (grad, par) -> model(par, grad),
                model.imply.start_val,
                model.diff.algorithm,
                model.diff.options)
    else
        result = Optim.optimize(
                model,
                (grad, par) -> model(par, grad),
                (H, par) -> model(par, H),
                model.imply.start_val,
                model.diff.algorithm,
                model.diff.options)
    end
    return result
end

function sem_fit(model::A, g!) where
    {A <: AbstractSem}
    result = Optim.optimize(
                par -> model(par),
                g!,
                model.imply.start_val,
                model.diff.algorithm,
                model.diff.options)#;
                #inplace = false)
    return result
end

function sem_fit(model::A, start_val::B) where
    {A <: AbstractSem, B <: AbstractArray}
    result = Optim.optimize(
                par -> model(par),
                start_val,
                model.sem_vec[1].diff.algorithm,
                model.sem_vec[1].diff.options)
    return result
end

function sem_fit(model::A, g!, h!) where
    {A <: AbstractSem}
    result = Optim.optimize(
                par -> model(par),
                g!,
                h!,
                model.imply.start_val,
                Newton(),
                model.diff.options)#;
                #inplace = false)
    return result
end

#function sem_fit(model::Sem{O, I, L, D}) where
#    {O <: SemObs, L <: Loss, I <: Imply, D <: SemReverseDiff}
#    result = optimize(
#                par -> model(par),
#                par -> Zygote.gradient(model, par)[1],
#                model.imply.start_val,
#                model.diff.algorithm,
#                model.diff.options;
#                inplace = false)
#    return result
#end

#function sem_fit(model::A, start_val) where
#    {A <: AbstractSem}
#    result = optimize(
#                par -> model(par),
#                start_val,
#                model.sem_vec[1].diff.algorithm,
#                model.sem_vec[1].diff.options)
#    return result
#end