text:
"""
SnpArray
SnpArray(bednm, m)
SnpArray(plknm)
SnpArray(undef, m, n)

Raw .bed file as a shared, memory-mapped Matrix{UInt8}.  The number of rows, `m`
is stored separately because it is not uniquely determined by the size of the `data` field.
"""
struct SnpArray <: AbstractMatrix{UInt8}
    data::Matrix{UInt8}
    columncounts::Matrix{Int}
    rowcounts::Matrix{Int}
    m::Int
end

AbstractSnpArray = Union{SnpArray, SubArray{UInt8, 1, SnpArray}, SubArray{UInt8, 2, SnpArray}}

function SnpArray(bednm::AbstractString, m::Integer, args...; kwargs...)
    checkplinkfilename(bednm, "bed")
    data = makestream(bednm, args...; kwargs...) do io
        read(io, UInt16) == 0x1b6c || throw(ArgumentError("wrong magic number in file $bednm"))
        read(io, UInt8) == 0x01 || throw(ArgumentError(".bed file, $bednm, is not in correct orientation"))
        if endswith(bednm, ".bed")
            return Mmap.mmap(io)
        else
            return read(io)
        end
    end
    drows = (m + 3) >> 2   # the number of rows in the Matrix{UInt8}
    n, r = divrem(length(data), drows)
    iszero(r) || throw(ArgumentError("filesize of $bednm is not a multiple of $drows"))
    SnpArray(reshape(data, (drows, n)), zeros(Int, (4, n)), zeros(Int, (4, m)), m)
end

function SnpArray(
    bednm::AbstractString,
    args...;
    famnm::Union{AbstractString, Nothing}=nothing,
    kwargs...)
    checkplinkfilename(bednm, "bed")
    # check user supplied fam filename
    if famnm == nothing
        famnm = replace(bednm, ".bed" => ".fam")
        isfile(famnm) || throw(ArgumentError("fam file not found"))
    else # user supplied fam file name
        checkplinkfilename(famnm, "fam")
        isfile(famnm) || throw(ArgumentError("file $famnm not found"))
    end
    m = makestream(famnm) do stream
        countlines(stream)
    end
    SnpArray(bednm, m, args...; kwargs...)
end

function SnpArray(::UndefInitializer, m::Integer, n::Integer)
    SnpArray(Matrix{UInt8}(undef, ((m + 3) >> 2, n)), zeros(Int, (4, n)), zeros(Int, (4, m)), m)
end

function SnpArray(file::AbstractString, s::SnpArray)
    makestream(file, "w+") do io
        write(io, 0x1b6c)
        write(io, 0x01)
        write(io, s.data)
    end
    SnpArray(file, s.m, "r+")
end

function SnpArray(file::AbstractString, m::Integer, n::Integer)
    makestream(file, "w+") do io
        write(io, 0x1b6c)
        write(io, 0x01)
        write(io, fill(0x00, ((m + 3) >> 2, n)))
    end
    SnpArray(file, m, "r+")
end

StatsBase.counts(s::AbstractSnpArray; dims=:) = _counts(s, dims)

function _counts(s::AbstractSnpArray, dims::Integer)
    if isone(dims)
        cc = (typeof(s) == SnpArray) ? s.columncounts : zeros(Int, (4, size(s, 2)))
        if all(iszero, cc)
            m, n = size(s)
            @inbounds for j in 1:n
                for i in 1:m
                    cc[s[i, j] + 1, j] += 1
                end
            end
        end
        return cc
    elseif dims == 2
        rc = (typeof(s) == SnpArray) ? s.rowcounts : zeros(Int, (4, size(s, 1)))
        if all(iszero, rc)
            m, n = size(s)
            @inbounds for j in 1:n
                for i in 1:m
                    rc[s[i, j] + 1, i] += 1
                end
            end
        end
        return rc
    else
        throw(ArgumentError("counts(s::SnpArray, dims=k) only defined for k = 1 or 2"))
    end
end

_counts(s::AbstractSnpArray, ::Colon) = sum(_counts(s, 1), dims=2)

function Base.getindex(s::SnpArray, i::Int)  # Linear indexing
    d, r = divrem(i - 1, s.m)
    s[r + 1, d + 1]
end

@inline function Base.getindex(s::SnpArray, i::Integer, j::Integer)
    @boundscheck checkbounds(s, i, j)
    ip3 = i + 3
    (s.data[ip3 >> 2, j] >> ((ip3 & 0x03) << 1)) & 0x03
end

function Base.setindex!(s::SnpArray, x::UInt8, i::Int)  # Linear indexing
    d, r = divrem(i - 1, s.m)
    Base.setindex!(s, x, r + 1, d + 1)
end

@inline function Base.setindex!(s::SnpArray, x::UInt8, i::Integer, j::Integer)
    @boundscheck checkbounds(s, i, j)
    ip3 = i + 3
    shft = (ip3 & 0x03) << 1
    mask = ~(0x03 << shft)
    s.data[ip3 >> 2, j] = (s.data[ip3 >> 2, j] & mask) | (x << shft)
    x
end

Base.eltype(s::SnpArray) = UInt8

Base.length(s::SnpArray) = s.m * size(s.data, 2)

Statistics.mean(s::AbstractSnpArray; dims=:, model=ADDITIVE_MODEL) = _mean(s, dims, model)

function _mean(s::AbstractSnpArray,  dims::Integer, ::typeof(ADDITIVE_MODEL))
    m, n = size(s)
    if isone(dims)
        cc = _counts(s, 1)   # need to use extractor to force evaluation if needed
        means = Matrix{Float64}(undef, (1, n))
        @inbounds for j in 1:n
            means[j] = (cc[3, j] + 2cc[4, j]) / (cc[1, j] + cc[3, j] + cc[4, j])
        end
        return means
    elseif dims == 2
        rc = _counts(s, 2)
        means = Matrix{Float64}(undef, (m, 1))
        @inbounds for i in 1:m
            means[i] = (rc[3, i] + 2rc[4, i]) / (rc[1, i] + rc[3, i] + rc[4, i])
        end
        return means
    else
        throw(ArgumentError("mean(s::SnpArray, dims=k) only defined for k = 1 or 2"))
    end
end

function _mean(s::AbstractSnpArray,  dims::Integer, ::typeof(DOMINANT_MODEL))
    m, n = size(s)
    if isone(dims)
        cc = _counts(s, 1)   # need to use extractor to force evaluation if needed
        means = Matrix{Float64}(undef, (1, n))
        @inbounds for j in 1:n
            means[j] = (cc[3, j] + cc[4, j]) / (cc[1, j] + cc[3, j] + cc[4, j])
        end
        return means
    elseif dims == 2
        rc = _counts(s, 2)
        means = Matrix{Float64}(undef, (m, 1))
        @inbounds for i in 1:m
            means[i] = (rc[3, i] + rc[4, i]) / (rc[1, i] + rc[3, i] + rc[4, i])
        end
        return means
    else
        throw(ArgumentError("mean(s::SnpArray, dims=k) only defined for k = 1 or 2"))
    end
end

function _mean(s::AbstractSnpArray,  dims::Integer, ::typeof(RECESSIVE_MODEL))
    m, n = size(s)
    if isone(dims)
        cc = _counts(s, 1)   # need to use extractor to force evaluation if needed
        means = Matrix{Float64}(undef, (1, n))
        @inbounds for j in 1:n
            means[j] = cc[4, j] / (cc[1, j] + cc[3, j] + cc[4, j])
        end
        return means
    elseif dims == 2
        rc = _counts(s, 2)
        means = Matrix{Float64}(undef, (m, 1))
        @inbounds for i in 1:m
            means[i] = rc[4, i] / (rc[1, i] + rc[3, i] + rc[4, i])
        end
        return means
    else
        throw(ArgumentError("mean(s::SnpArray, dims=k) only defined for k = 1 or 2"))
    end
end

function _mean(s::AbstractSnpArray, ::Colon, ::typeof(ADDITIVE_MODEL))
    rc = _counts(s, 2)
    @views (sum(rc[3, :]) + 2sum(rc[4, :])) /  sum(rc[[1, 3, 4], :])
end

function _mean(s::AbstractSnpArray, ::Colon, ::typeof(DOMINANT_MODEL))
    rc = _counts(s, 2)
    @views (sum(rc[3, :]) + sum(rc[4, :])) /  sum(rc[[1, 3, 4], :])
end

function _mean(s::AbstractSnpArray, ::Colon, ::typeof(RECESSIVE_MODEL))
    rc = _counts(s, 2)
    @views sum(rc[4, :]) /  sum(rc[[1, 3, 4], :])
end

Base.size(s::SnpArray) = s.m, size(s.data, 2)

Base.size(s::SnpArray, k::Integer) =
k == 1 ? s.m : k == 2 ? size(s.data, 2) : k > 2 ? 1 : error("Dimension k out of range")

# TODO: need to implement `var` for different SNP models

Statistics.var(s::AbstractSnpArray; corrected::Bool=true, mean=nothing, dims=:) = _var(s, corrected, mean, dims)

function _var(s::AbstractSnpArray, corrected::Bool, mean, dims::Integer)
    m, n = size(s)
    means = something(mean, Statistics.mean(s, dims=dims))
    if isone(dims)
        cc = _counts(s, 1)
        vars = Matrix{Float64}(undef, (1, n))
        for j in 1:n
            mnj = means[j]
            vars[j] = (abs2(mnj) * cc[1, j] + abs2(1.0 - mnj) * cc[3, j] + abs2(2.0 - mnj) * cc[4,j]) /
            (cc[1, j] + cc[3, j] + cc[4, j] - (corrected ? 1 : 0))
        end
        return vars
    elseif dims == 2
        rc = _counts(s, 2)
        vars = Matrix{Float64}(undef, (m, 1))
        for i in 1:m
            mni = means[i]
            vars[i] = (abs2(mni) * rc[1, i] + abs2(1.0 - mni) * rc[3, i] + abs2(2.0 - mni) * rc[4, i]) /
            (rc[1, i] + rc[3, i] + rc[4, i] - (corrected ? 1 : 0))
        end
        return vars
    end
    throw(ArgumentError("var(s::SnpArray, dims=k) only defined for k = 1 or 2"))
end

"""
maf!(out, s)

Populate `out` with minor allele frequencies of SnpArray `s`.
"""
function maf!(out::AbstractVector{T}, s::AbstractSnpArray) where T <: AbstractFloat
    cc = _counts(s, 1)
    @inbounds for j in 1:size(s, 2)
        out[j] = (cc[3, j] + 2cc[4, j]) / 2(cc[1, j] + cc[3, j] + cc[4, j])
        (out[j] > 0.5) && (out[j] = 1 - out[j])
    end
    out
end
maf(s::AbstractSnpArray) = maf!(Vector{Float64}(undef, size(s, 2)), s)

"""
minorallele(out, s)

Populate `out` with minor allele indicators. `out[j] == true` means A2 is the minor
allele of `j`th column; `out[j] == false` means A1 is the minor allele.
"""
function minorallele!(out::AbstractVector{Bool}, s::AbstractSnpArray)
    cc = _counts(s, 1)
    @inbounds for j in 1:size(s, 2)
        out[j] = cc[1, j] > cc[4, j]
    end
    out
end

"""
minorallele(s)

Calculate minor allele indicators. `out[j] == true` means A2 is the minor
allele of `j`th column; `out[j] == false` means A1 is the minor allele.
"""
minorallele(s::AbstractSnpArray) = minorallele!(BitVector(undef, size(s, 2)), s)

@inline function convert(::Type{T}, x::UInt8, ::Val{1}) where T <: AbstractFloat
    iszero(x) ? zero(T) : isone(x) ? T(NaN) : T(x - 1)
end

@inline function convert(::Type{T}, x::UInt8, ::Val{2}) where T <: AbstractFloat
    iszero(x) ? zero(T) : isone(x) ? T(NaN) : one(T)
end

@inline function convert(::Type{T}, x::UInt8, ::Val{3}) where T <: AbstractFloat
    (iszero(x) || x == 2) ? zero(T) : isone(x) ? T(NaN) : one(T)
end

"""
Base.copyto!(v, s, model=ADDITIVE_MODEL, center=false, scale=false, impute=false)

Copy SnpArray `s` to numeric vector or matrix `v`.

# Arguments
- `model::Union{Val{1}, Val{2}, Val{3}}=ADDITIVE_MODEL`: `ADDITIVE_MODEL` (default), `DOMINANT_MODEL`, or `RECESSIVE_MODEL`.
- `center::Bool=false`: center column by mean.
- `scale::Bool=false`: scale column by theoretical variance.
- `impute::Bool=falase`: impute missing values by column mean.
"""
function Base.copyto!(
    v::AbstractVecOrMat{T},
    s::AbstractSnpArray;
    model::Union{Val{1}, Val{2}, Val{3}} = ADDITIVE_MODEL,
    center::Bool = false,
    scale::Bool = false,
    impute::Bool = false
    ) where T <: AbstractFloat
    m, n = size(s, 1), size(s, 2)
    # no center, scale, or impute
    if !center && !scale && !impute
        @inbounds for j in 1:n
            @simd for i in 1:m
                v[i, j] = SnpArrays.convert(T, s[i, j], model)
            end
        end
        return v
    end
    # center, scale, impute
    @inbounds for j in 1:n
        μj, mj = zero(T), 0
        @simd for i in 1:m
            vij = SnpArrays.convert(T, s[i, j], model)
            v[i, j] = vij
            μj += isnan(vij) ? zero(T) : vij
            mj += isnan(vij) ? 0 : 1
        end
        μj /= mj
        σj = model == ADDITIVE_MODEL ? sqrt(μj * (1 - μj / 2)) : sqrt(μj * (1 - μj))
        @simd for i in 1:m
            impute && isnan(v[i, j]) && (v[i, j] = μj)
            center && (v[i, j] -= μj)
            scale && σj > 0 && (v[i, j] /= σj)
        end
    end
    return v
end

"""
Base.convert(t, s, model=ADDITIVE_MODEL, center=false, scale=false, impute=false)

Convert a SnpArray `s` to a numeric vector or matrix of same shape as `s`.

# Arguments
- `t::Type{AbstractVecOrMat{T}}`: Vector or matrix type.
- `model::Union{Val{1}, Val{2}, Val{3}}=ADDITIVE_MODEL`: `ADDITIVE_MODEL` (default), `DOMINANT_MODEL`, or `RECESSIVE_MODEL`.
- `center::Bool=false`: center column by mean.
- `scale::Bool=false`: scale column by theoretical variance.
- `impute::Bool=falase`: impute missing values by column mean.
"""
function Base.convert(
    ::Type{T},
    s::AbstractSnpArray;
    kwargs...) where T <: Array
    T(s; kwargs...)
end
Array{T,N}(s::AbstractSnpArray; kwargs...) where {T,N} = copyto!(Array{T,N}(undef, size(s)), s; kwargs...)

"""
missingpos(s::SnpArray)

Return a `SparseMatrixCSC{Bool,Int32}` of the same size as `s` indicating the positions with missing data.
"""
function missingpos(s::AbstractSnpArray)
    m, n = size(s)
    colptr = sizehint!(Int32[1], n + 1)
    rowval = Int32[]
    @inbounds for j in 1:n
        msngpos = Int32[]
        for i in 1:m
            isone(s[i, j]) && push!(msngpos, i)
        end
        append!(rowval, msngpos)
        push!(colptr, colptr[end] + length(msngpos))
    end
    SparseMatrixCSC(m, n, colptr, rowval, fill(true, length(rowval)))
end

function _missingrate!(out::AbstractVector{<:AbstractFloat}, s::AbstractSnpArray,  dims::Integer)
    m, n = size(s)
    if isone(dims)
        cc = _counts(s, 1)   # need to use extractor to force evaluation if needed
        @inbounds for j in 1:n
            out[j] = cc[2, j] / m
        end
    elseif dims == 2
        rc = _counts(s, 2)
        @inbounds for i in 1:m
            out[i] = rc[2, i] / n
        end
    else
        throw(ArgumentError("_missingrate(out, s::SnpArray, dims=k) only defined for k = 1 or 2"))
    end
    out
end

function missingrate(s::AbstractSnpArray, dims::Integer)
    if isone(dims)
        return _missingrate!(Vector{Float64}(undef, size(s, 2)), s, 1)
    elseif dims == 2
        return _missingrate!(Vector{Float64}(undef, s.m), s, 2)
    else
        throw(ArgumentError("missingrate(s::SnpArray, dims=k) only defined for k = 1 or 2"))
    end
end


text:
module Parsers

using Dates

include("utils.jl")

"""
    `Parsers.Options` is a structure for holding various parsing settings when calling `Parsers.parse`, `Parsers.tryparse`, and `Parsers.xparse`. They include:

  * `sentinel=nothing`: valid values include: `nothing` meaning don't check for sentinel values; `missing` meaning an "empty field" should be considered a sentinel value; or a `Vector{String}` of the various string values that should each be checked as a sentinel value. Note that sentinels will always be checked longest to shortest, with the longest valid match taking precedence.
  * `wh1=' '`: the first ascii character to be considered when ignoring leading/trailing whitespace in value parsing
  * `wh2='\t'`: the second ascii character to be considered when ignoring leading/trailing whitespace in value parsing
  * `openquotechar='"'`: the ascii character that signals a "quoted" field while parsing; subsequent characters will be treated as non-significant until a valid `closequotechar` is detected
  * `closequotechar='"'`: the ascii character that signals the end of a quoted field
  * `escapechar='"'`: an ascii character used to "escape" a `closequotechar` within a quoted field
  * `delim=nothing`: if `nothing`, no delimiter will be checked for; if a `Char` or `String`, a delimiter will be checked for directly after parsing a value or `closequotechar`; a newline (`\n`), return (`\r`), or CRLF (`"\r\n"`) are always considered "delimiters", in addition to EOF
  * `decimal='.'`: an ascii character to be used when parsing float values that separates a decimal value
  * `trues=nothing`: if `nothing`, `Bool` parsing will only check for the string `true` or an `Integer` value of `1` as valid values for `true`; as a `Vector{String}`, each string value will be checked to indicate a valid `true` value
  * `falses=nothing`: if `nothing`, `Bool` parsing will only check for the string `false` or an `Integer` value of `0` as valid values for `false`; as a `Vector{String}`, each string value will be checked to indicate a valid `false` value
  * `dateformat=nothing`: if `nothing`, `Date`, `DateTime`, and `Time` parsing will use a default `Dates.DateFormat` object while parsing; a `String` or `Dates.DateFormat` object can be provided for custom format parsing
  * `ignorerepeated=false`: if `true`, consecutive delimiter characters or strings will be consumed until a non-delimiter is encountered; if `false`, only a single delimiter character/string will be consumed. Useful for fixed-width delimited files where fields are padded with delimiters
  * `quoted=false`: whether parsing should check for `openquotechar` and `closequotechar` characters to signal quoted fields
  * `debug=false`: if `true`, various debug logging statements will be printed while parsing; useful when diagnosing why parsing returns certain `Parsers.ReturnCode` values
"""
struct Options{ignorerepeated, Q, debug, S, D, DF}
    sentinel::S # Union{Nothing, Missing, Vector{Tuple{Ptr{UInt8}, Int}}}
    wh1::UInt8
    wh2::UInt8
    oq::UInt8
    cq::UInt8
    e::UInt8
    delim::D # Union{Nothing, UInt8, Tuple{Ptr{UInt8}, Int}}
    decimal::UInt8
    trues::Union{Nothing, Vector{Tuple{Ptr{UInt8}, Int}}}
    falses::Union{Nothing, Vector{Tuple{Ptr{UInt8}, Int}}}
    dateformat::DF # Union{Nothing, Dates.DateFormat}
    strict::Bool
    silencewarnings::Bool
end

prepare(x::Vector{String}) = sort!(map(ptrlen, x), by=x->x[2], rev=true)
asciival(c::Char) = isascii(c)
asciival(b::UInt8) = b < 0x80

function Options(
            sentinel::Union{Nothing, Missing, Vector{String}},
            wh1::Union{UInt8, Char},
            wh2::Union{UInt8, Char},
            oq::Union{UInt8, Char},
            cq::Union{UInt8, Char},
            e::Union{UInt8, Char},
            delim::Union{Nothing, UInt8, Char, String},
            decimal::Union{UInt8, Char},
            trues::Union{Nothing, Vector{String}},
            falses::Union{Nothing, Vector{String}},
            dateformat::Union{Nothing, String, Dates.DateFormat},
            ignorerepeated, quoted, debug, strict=false, silencewarnings=false)
    asciival(wh1) && asciival(wh2) || throw(ArgumentError("whitespace characters must be ASCII"))
    asciival(oq) && asciival(cq) && asciival(e) || throw(ArgumentError("openquotechar, closequotechar, and escapechar must be ASCII characters"))
    (oq == delim) || (cq == delim) || (e == delim) && throw(ArgumentError("delim argument must be different than openquotechar, closequotechar, and escapechar arguments"))
    if sentinel isa Vector{String}
        for sent in sentinel
            if startswith(sent, string(Char(wh1))) || startswith(sent, string(Char(wh2)))
                throw(ArgumentError("sentinel value isn't allowed to start with wh1 or wh2 characters"))
            end
            if startswith(sent, string(Char(oq))) || startswith(sent, string(Char(cq)))
                throw(ArgumentError("sentinel value isn't allowed to start with openquotechar, closequotechar, or escapechar characters"))
            end
            if (delim isa UInt8 || delim isa Char) && startswith(sent, string(Char(delim)))
                throw(ArgumentError("sentinel value isn't allowed to start with a delimiter character"))
            elseif delim isa String && startswith(sent, delim)
                throw(ArgumentError("sentinel value isn't allowed to start with a delimiter string"))
            end
        end
    end
    sent = sentinel === nothing || sentinel === missing ? sentinel : prepare(sentinel)
    del = delim === nothing ? nothing : delim isa String ? ptrlen(delim) : delim % UInt8
    if del isa UInt8
        ((wh1 % UInt8) == del || (wh2 % UInt8) == del) && throw(ArgumentError("whitespace characters (`wh1=' '` and `wh2='\\t'` default keyword arguments) must be different than delim argument"))
    end
    trues = trues === nothing ? nothing : prepare(trues)
    falses = falses === nothing ? nothing : prepare(falses)
    df = dateformat === nothing ? nothing : dateformat isa String ? Dates.DateFormat(dateformat) : dateformat
    return Options{ignorerepeated, quoted, debug, typeof(sent), typeof(del), typeof(df)}(sent, wh1 % UInt8, wh2 % UInt8, oq % UInt8, cq % UInt8, e % UInt8, del, decimal % UInt8, trues, falses, df, strict, silencewarnings)
end

Options(;
    sentinel::Union{Nothing, Missing, Vector{String}}=nothing,
    wh1::Union{UInt8, Char}=UInt8(' '),
    wh2::Union{UInt8, Char}=UInt8('\t'),
    openquotechar::Union{UInt8, Char}=UInt8('"'),
    closequotechar::Union{UInt8, Char}=UInt8('"'),
    escapechar::Union{UInt8, Char}=UInt8('"'),
    delim::Union{Nothing, UInt8, Char, String}=nothing,
    decimal::Union{UInt8, Char}=UInt8('.'),
    trues::Union{Nothing, Vector{String}}=nothing,
    falses::Union{Nothing, Vector{String}}=nothing,
    dateformat::Union{Nothing, String, Dates.DateFormat}=nothing,
    ignorerepeated::Bool=false,
    quoted::Bool=false,
    debug::Bool=false,
) = Options(sentinel, wh1, wh2, openquotechar, closequotechar, escapechar, delim, decimal, trues, falses, dateformat, ignorerepeated, quoted, debug)

const OPTIONS = Options(nothing, UInt8(' '), UInt8('\t'), UInt8('"'), UInt8('"'), UInt8('"'), nothing, UInt8('.'), nothing, nothing, nothing, false, false, false)
const XOPTIONS = Options(missing, UInt8(' '), UInt8('\t'), UInt8('"'), UInt8('"'), UInt8('"'), UInt8(','), UInt8('.'), nothing, nothing, nothing, false, true, false)

# high-level convenience functions like in Base
"Attempt to parse a value of type `T` from string `buf`. Throws `Parsers.Error` on parser failures and invalid values."
function parse(::Type{T}, buf::Union{AbstractVector{UInt8}, AbstractString, IO}, options=OPTIONS; pos::Integer=1, len::Integer=buf isa IO ? 0 : sizeof(buf)) where {T}
    x, code, vpos, vlen, tlen = xparse(T, buf isa AbstractString ? codeunits(buf) : buf, pos, len, options)
    return ok(code) ? x : throw(Error(buf, T, code, pos, tlen))
end

"Attempt to parse a value of type `T` from `buf`. Returns `nothing` on parser failures and invalid values."
function tryparse(::Type{T}, buf::Union{AbstractVector{UInt8}, AbstractString, IO}, options=OPTIONS; pos::Integer=1, len::Integer=buf isa IO ? 0 : sizeof(buf)) where {T}
    x, code, vpos, vlen, tlen = xparse(T, buf isa AbstractString ? codeunits(buf) : buf, pos, len, options)
    return ok(code) ? x : nothing
end

default(::Type{T}) where {T <: Integer} = zero(T)
default(::Type{T}) where {T <: AbstractFloat} = T(0.0)
default(::Type{T}) where {T <: Dates.TimeType} = T(0)

# for testing purposes only, it's much too slow to dynamically create Options for every xparse call
"""
    Parsers.xparse(T, buf, pos, len, options) => (x, code, startpos, value_len, total_len)

    The core parsing function for any type `T`. Takes a `buf`, which can be a `Vector{UInt8}`, `Base.CodeUnits`,
    or an `IO`. `pos` is the byte position to begin parsing at. `len` is the total # of bytes in `buf` (signaling eof).
    `options` is an instance of `Parsers.Options`.

    `Parsers.xparse` returns a tuple of 5 values:
      * `x` is a value of type `T`, even if parsing does not succeed
      * `code` is a bitmask of parsing codes, use `Parsers.codes(code)` or `Parsers.text(code)` to see the various bit values set. See `?Parsers.ReturnCode` for additional details on the various parsing codes
      * `startpos`: the starting byte position of the value being parsed; will always equal the start `pos` passed in, except for quoted field where it will point instead to the first byte after the open quote character
      * `value_len`: the # of bytes consumed while parsing a value, will be equal to the total number of bytes consumed, except for quoted or delimited fields where the quote and delimiter characters will be subtracted out
      * `total_len`: the total # of bytes consumed while parsing a value, including any quote or delimiter characters; this can be added to the starting `pos` to allow calling `Parsers.xparse` again for a subsequent field/value
"""
function xparse end

function xparse(::Type{T}, buf::Union{AbstractVector{UInt8}, AbstractString, IO}; pos::Integer=1, len::Integer=buf isa IO ? 0 : sizeof(buf), sentinel=nothing, wh1::Union{UInt8, Char}=UInt8(' '), wh2::Union{UInt8, Char}=UInt8('\t'), quoted::Bool=true, openquotechar::Union{UInt8, Char}=UInt8('"'), closequotechar::Union{UInt8, Char}=UInt8('"'), escapechar::Union{UInt8, Char}=UInt8('"'), ignorerepeated::Bool=false, delim::Union{UInt8, Char, Tuple{Ptr{UInt8}, Int}, AbstractString, Nothing}=UInt8(','), decimal::Union{UInt8, Char}=UInt8('.'), trues=nothing, falses=nothing, dateformat::Union{Nothing, String, Dates.DateFormat}=nothing, debug::Bool=false) where {T}
    options = Options(sentinel, wh1, wh2, openquotechar, closequotechar, escapechar, delim, decimal, trues, falses, dateformat, ignorerepeated, quoted, debug)
    return xparse(T, buf isa AbstractString ? codeunits(buf) : buf, pos, len, options)
end

function xparse(::Type{T}, buf::AbstractString, pos, len, options::Options=XOPTIONS) where {T}
    return xparse(T, codeunits(buf), pos, len, options)
end

@inline function xparse(::Type{T}, source::Union{AbstractVector{UInt8}, IO}, pos, len, options::Options{ignorerepeated, Q, debug, S, D, DF}=XOPTIONS) where {T, ignorerepeated, Q, debug, S, D, DF}
    startpos = vstartpos = vpos = pos
    sentinel = options.sentinel
    code = SUCCESS
    x = default(T)
    quoted = false
    sentinelpos = 0
    if debug
        println("parsing $T, pos=$pos, len=$len")
    end
    if eof(source, pos, len)
        code = (sentinel === missing ? SENTINEL : INVALID) | EOF
        @goto donedone
    end
    b = peekbyte(source, pos)
    if debug
        println("1) parsed: '$(escape_string(string(Char(b))))'")
    end
    # strip leading whitespace
    while b == options.wh1 || b == options.wh2
        if debug
            println("stripping leading whitespace")
        end
        pos += 1
        incr!(source)
        if eof(source, pos, len)
            code = INVALID | EOF
            @goto donedone
        end
        b = peekbyte(source, pos)
        if debug
            println("2) parsed: '$(escape_string(string(Char(b))))'")
        end
    end
    # check for start of quoted field
    if Q
        quoted = b == options.oq
        if quoted
            if debug
                println("detected open quote character")
            end
            code = QUOTED
            pos += 1
            vstartpos = pos
            incr!(source)
            if eof(source, pos, len)
                code |= INVALID_QUOTED_FIELD
                @goto donedone
            end
            b = peekbyte(source, pos)
            if debug
                println("3) parsed: '$(escape_string(string(Char(b))))'")
            end
            # ignore whitespace within quoted field
            while b == options.wh1 || b == options.wh2
                if debug
                    println("stripping whitespace within quoted field")
                end
                pos += 1
                incr!(source)
                if eof(source, pos, len)
                    code |= INVALID_QUOTED_FIELD | EOF
                    @goto donedone
                end
                b = peekbyte(source, pos)
                if debug
                    println("4) parsed: '$(escape_string(string(Char(b))))'")
                end
            end
        end
    end
    # check for sentinel values if applicable
    if sentinel !== nothing && sentinel !== missing
        if debug
            println("checking for sentinel value")
        end
        sentinelpos = checksentinel(source, pos, len, sentinel, debug)
    end
    x, code, pos = typeparser(T, source, pos, len, b, code, options)
    if sentinel !== nothing && sentinel !== missing && sentinelpos >= pos
        # if we matched a sentinel value that was as long or longer than our type value
        code &= ~(OK | INVALID | OVERFLOW)
        pos = sentinelpos
        fastseek!(source, pos - 1)
        code |= SENTINEL
        if eof(source, pos, len)
            code |= EOF
        end
    elseif sentinel === missing && pos == vstartpos
        code &= ~(OK | INVALID)
        code |= SENTINEL
    end
    vpos = pos
    if (code & EOF) == EOF
        if quoted
            # if we detected a quote character, it's an invalid quoted field due to eof in the middle
            code |= INVALID_QUOTED_FIELD
        end
        @goto donedone
    end

@label donevalue
    b = peekbyte(source, pos)
    if debug
        println("finished $T value parsing: pos=$pos, current character: '$(escape_string(string(Char(b))))'")
    end
    # donevalue means we finished parsing a value or sentinel, but didn't reach len, b is still the current byte
    # strip trailing whitespace
    while b == options.wh1 || b == options.wh2
        if debug
            println("stripping trailing whitespace")
        end
        pos += 1
        vpos += 1
        incr!(source)
        if eof(source, pos, len)
            code |= EOF
            if quoted
                code |= INVALID_QUOTED_FIELD
            end
            @goto donedone
        end
        b = peekbyte(source, pos)
        if debug
            println("8) parsed: '$(escape_string(string(Char(b))))'")
        end
    end
    if Q
        # for quoted fields, find the closing quote character
        # we should be positioned at the correct place to find the closing quote character if everything is as it should be
        # if we don't find the quote character immediately, something's wrong, so mark INVALID
        if quoted
            if debug
                println("looking for close quote character")
            end
            same = options.cq == options.e
            first = true
            while true
                vpos = pos
                pos += 1
                incr!(source)
                if same && b == options.e
                    if eof(source, pos, len)
                        code |= EOF
                        if !first
                            code |= INVALID
                        end
                        @goto donedone
                    elseif peekbyte(source, pos) != options.cq
                        if !first
                            code |= INVALID
                        end
                        break
                    end
                    code |= ESCAPED_STRING
                    pos += 1
                    incr!(source)
                elseif b == options.e
                    if eof(source, pos, len)
                        code |= INVALID_QUOTED_FIELD | EOF
                        @goto donedone
                    end
                    code |= ESCAPED_STRING
                    pos += 1
                    incr!(source)
                elseif b == options.cq
                    if !first
                        code |= INVALID
                    end
                    if eof(source, pos, len)
                        code |= EOF
                        @goto donedone
                    end
                    break
                end
                if eof(source, pos, len)
                    code |= INVALID_QUOTED_FIELD | EOF
                    @goto donedone
                end
                first = false
                b = peekbyte(source, pos)
                if debug
                    println("9) parsed: '$(escape_string(string(Char(b))))'")
                end
            end
            b = peekbyte(source, pos)
            if debug
                println("10) parsed: '$(escape_string(string(Char(b))))'")
            end
            # ignore whitespace after quoted field
            while b == options.wh1 || b == options.wh2
                if debug
                    println("stripping trailing whitespace after close quote character")
                end
                pos += 1
                incr!(source)
                if eof(source, pos, len)
                    code |= EOF
                    @goto donedone
                end
                b = peekbyte(source, pos)
                if debug
                    println("11) parsed: '$(escape_string(string(Char(b))))'")
                end
            end
        end
    end

    if options.delim !== nothing
        delim = options.delim
        # now we check for a delimiter; if we don't find it, keep parsing until we do
        if debug
            println("checking for delimiter: pos=$pos")
        end
        if !ignorerepeated
            # we're checking for a single appearance of a delimiter
            if delim isa UInt8
                if b == delim
                    pos += 1
                    incr!(source)
                    code |= DELIMITED
                    @goto donedone
                end
            else
                predelimpos = pos
                pos = checkdelim(source, pos, len, delim)
                if pos > predelimpos
                    # found the delimiter we were looking for
                    code |= DELIMITED
                    @goto donedone
                end
            end
        else
            # keep parsing as long as we keep matching delim
            if delim isa UInt8
                matched = false
                while b == delim
                    matched = true
                    pos += 1
                    incr!(source)
                    if eof(source, pos, len)
                        code |= DELIMITED
                        @goto donedone
                    end
                    b = peekbyte(source, pos)
                    if debug
                        println("12) parsed: '$(escape_string(string(Char(b))))'")
                    end
                end
                if matched
                    # if a newline is next, consume it as well
                    if b == UInt8('\n')
                        pos += 1
                        incr!(source)
                        code |= NEWLINE | ifelse(eof(source, pos, len), EOF, SUCCESS)
                    elseif b == UInt8('\r')
                        pos += 1
                        incr!(source)
                        if !eof(source, pos, len) && peekbyte(source, pos) == UInt8('\n')
                            pos += 1
                            incr!(source)
                        end
                        code |= NEWLINE | ifelse(eof(source, pos, len), EOF, SUCCESS)
                    end
                    code |= DELIMITED
                    @goto donedone
                end
            else
                matched = false
                predelimpos = pos
                pos = checkdelim(source, pos, len, delim)
                while pos > predelimpos
                    matched = true
                    if eof(source, pos, len)
                        code |= DELIMITED
                        @goto donedone
                    end
                    predelimpos = pos
                    pos = checkdelim(source, pos, len, delim)
                end
                if matched
                    # if a newline is next, consume it as well
                    b = peekbyte(source, pos)
                    if b == UInt8('\n')
                        pos += 1
                        incr!(source)
                        code |= NEWLINE | ifelse(eof(source, pos, len), EOF, SUCCESS)
                    elseif b == UInt8('\r')
                        pos += 1
                        incr!(source)
                        if !eof(source, pos, len) && peekbyte(source, pos) == UInt8('\n')
                            pos += 1
                            incr!(source)
                        end
                        code |= NEWLINE | ifelse(eof(source, pos, len), EOF, SUCCESS)
                    end
                    code |= DELIMITED
                    @goto donedone
                end
            end
        end
        # didn't find delimiter, but let's check for a newline character
        if b == UInt8('\n')
            pos += 1
            incr!(source)
            code |= NEWLINE | ifelse(eof(source, pos, len), EOF, SUCCESS)
            @goto donedone
        elseif b == UInt8('\r')
            pos += 1
            incr!(source)
            if !eof(source, pos, len) && peekbyte(source, pos) == UInt8('\n')
                pos += 1
                incr!(source)
            end
            code |= NEWLINE | ifelse(eof(source, pos, len), EOF, SUCCESS)
            @goto donedone
        end
        # didn't find delimiter or newline, so we're invalid, keep parsing until we find delimiter, newline, or len
        quo = Int(!quoted)
        code |= INVALID_DELIMITER
        while true
            pos += 1
            vpos += quo
            incr!(source)
            if eof(source, pos, len)
                code |= EOF
                @goto donedone
            end
            b = peekbyte(source, pos)
            if debug
                println("13) parsed: '$(escape_string(string(Char(b))))'")
            end
            if !ignorerepeated
                if delim isa UInt8
                    if b == delim
                        pos += 1
                        incr!(source)
                        code |= DELIMITED
                        @goto donedone
                    end
                else
                    predelimpos = pos
                    pos = checkdelim(source, pos, len, delim)
                    if pos > predelimpos
                        # found the delimiter we were looking for
                        code |= DELIMITED
                        @goto donedone
                    end
                end
            else
                if delim isa UInt8
                    matched = false
                    while b == delim
                        matched = true
                        pos += 1
                        incr!(source)
                        if eof(source, pos, len)
                            code |= DELIMITED
                            @goto donedone
                        end
                        b = peekbyte(source, pos)
                        if debug
                            println("12) parsed: '$(escape_string(string(Char(b))))'")
                        end
                    end
                    if matched
                        code |= DELIMITED
                        @goto donedone
                    end
                else
                    predelimpos = pos
                    pos = checkdelim(source, pos, len, delim)
                    while pos > predelimpos
                        matched = true
                        if eof(source, pos, len)
                            code |= DELIMITED
                            @goto donedone
                        end
                        predelimpos = pos
                        pos = checkdelim(source, pos, len, delim)
                    end
                    if matched
                        code |= DELIMITED
                        @goto donedone
                    end
                end
            end
            # didn't find delimiter, but let's check for a newline character
            if b == UInt8('\n')
                pos += 1
                incr!(source)
                code |= NEWLINE | ifelse(eof(source, pos, len), EOF, SUCCESS)
                @goto donedone
            elseif b == UInt8('\r')
                pos += 1
                incr!(source)
                if !eof(source, pos, len) && peekbyte(source, pos) == UInt8('\n')
                    pos += 1
                    incr!(source)
                end
                code |= NEWLINE | ifelse(eof(source, pos, len), EOF, SUCCESS)
                @goto donedone
            end
        end
    end

@label donedone
    if debug
        println("finished parsing: $(codes(code))")
    end
    return x, code, Int64(vstartpos), Int64(vpos - vstartpos), Int64(pos - startpos)
end

function checkdelim!(buf, pos, len, options::Options{ignorerepeated}) where {ignorerepeated}
    pos > len && return pos
    delim = options.delim
    @inbounds b = buf[pos]
    valuepos = pos
    if !ignorerepeated
        # we're checking for a single appearance of a delimiter
        if delim isa UInt8
            b == delim && return pos + 1
        else
            pos = checkdelim(buf, pos, len, delim)
            pos > valuepos && return pos
        end
    else
        # keep parsing as long as we keep matching delim
        if delim isa UInt8
            matched = false
            while b == delim
                matched = true
                pos += 1
                pos > len && return pos
                @inbounds b = buf[pos]
            end
            matched && return pos
        else
            matched = false
            predelimpos = pos
            pos = checkdelim(buf, pos, len, delim)
            while pos > predelimpos
                matched = true
                pos > len && return pos
                predelimpos = pos
                pos = checkdelim(buf, pos, len, delim)
            end
            matched && return pos
        end
    end
    return pos
end

include("ints.jl")
include("floats.jl")
include("strings.jl")
include("bools.jl")
include("dates.jl")
include("ryu.jl")

function __init__()
    # floats.jl globals
    Threads.resize_nthreads!(ONES)
    foreach(x->MPZ.init!(x), ONES)
    Threads.resize_nthreads!(NUMS)
    foreach(x->MPZ.init!(x), NUMS)
    Threads.resize_nthreads!(QUOS)
    foreach(x->MPZ.init!(x), QUOS)
    Threads.resize_nthreads!(REMS)
    foreach(x->MPZ.init!(x), REMS)
    Threads.resize_nthreads!(SCLS)
    foreach(x->MPZ.init!(x), SCLS)
    return
end

end # module

ext: jl


text:
using ITensors,
      LinearAlgebra,
      QuadGK

const βc = 0.5*log(sqrt(2.0)+1.0)

function ising_free_energy(β::Real,J::Real=1.0)
  k = β*J
  c = cosh(2.0*k)
  s = sinh(2.0*k)
  xmin = 0.0
  xmax = π
  integrand(x) = log(c^2+sqrt(s^4+1-2*s^2*cos(x)))
  integral,err = quadgk(integrand, xmin, xmax)::Tuple{Float64,Float64}
  return -(log(2.0)+integral/π)/(2.0*β)
end

function ising_magnetization(β::Real)
 β>βc && return (1.0-sinh(2.0*β)^(-4))^(1/8)
 return 0.0
end

function ising_mpo(sh::Tuple{Index,Index},sv::Tuple{Index,Index},
                   β::Real,J::Real=1.0;
                   sz::Bool=false,dual_lattice::Bool=true)
  d = dim(sh[1])
  T = ITensor(sh[1],sh[2],sv[1],sv[2])
  if dual_lattice
    for i = 1:d
      T[i,i,i,i] = 1.0
    end
    sz && (T[1,1,1,1] = -T[1,1,1,1])
    Q = [exp(β*J) exp(-β*J); exp(-β*J) exp(β*J)]
    D,U = eigen(Symmetric(Q))
    √Q = U*Diagonal(sqrt.(D))*U'
    Xh1 = ITensor(vec(√Q),sh[1],sh[1]')
    Xh2 = ITensor(vec(√Q),sh[2],sh[2]')
    Xv1 = ITensor(vec(√Q),sv[1],sv[1]')
    Xv2 = ITensor(vec(√Q),sv[2],sv[2]')
    T = replacetags(T*Xh1*Xh2*Xv1*Xv2,"1","0")
  else
    sig(s) = 1.0-2.0*(s-1)
    E0 = -4.0
    for s1 = 1:d, s2 = 1:d, s3 = 1:d, s4 = 1:d
      E = sig(s1)*sig(s2)+sig(s2)*sig(s3)+sig(s3)*sig(s4)+sig(s4)*sig(s1)
      val = exp(-β*(E-E0))
      T[sh[1](s1),sv[2](s2),sh[2](s3),sv[1](s4)] = val
    end
  end
  return T
end


text:
using Test
using OpenQuantumSystems
using Random, SparseArrays, LinearAlgebra, StableRNGs

@testset "rate_constant" begin

    Random.seed!(StableRNG(0), 1)

    D(op1::Array, op2::Array) = abs(norm(op1 - op2))
    D(x1::StateVector, x2::StateVector) = norm(x2 - x1)
    D(op1::AbstractOperator, op2::AbstractOperator) =
        abs(tracedistance_nh(dense(op1), dense(op2)))
    D(op1::AbstractSuperOperator, op2::AbstractSuperOperator) =
        abs(tracedistance_nh(dense(op1), dense(op2)))

    mode1 = Mode(0.2, 1.0)
    mode2 = Mode(0.3, 2.0)
    Energy = [0.0, 200.0]
    mol1 = Molecule([mode1], 3, [2.0, 200.0])
    mol2 = Molecule([mode2], 3, [3.0, 300.0])
    aggCore = AggregateCore([mol1, mol2])
    aggCore.coupling[2, 3] = 50
    aggCore.coupling[3, 2] = 50
    agg = setupAggregate(aggCore)
    aggTools = agg.tools
    aggOperators = agg.operators

    Ham_B = agg.operators.Ham_B
    Ham_I = agg.operators.Ham_I
    Ham_0 = agg.operators.Ham_0
    Ham = agg.operators.Ham

    Ham_0_lambda, Ham_0_S = eigen(Ham_0.data)
    Ham_0_Sinv = inv(Ham_0_S)
    Ham_0_lambda = diagm(Ham_0_lambda)

    basis = agg.tools.basis
    indicesLen = agg.tools.bSize
    indices = agg.tools.indices
    indicesMap = agg.tools.indicesMap
    FCFact = agg.tools.FCfactors
    FCProd = agg.tools.FCproduct

    tspan = get_tspan(0., 0.02, 100)
    W0, rho0, W0_bath = ultrafast_laser_excitation(10., [0.0, 0.3, 0.7], agg)

    tmp1 = copy(W0.data)
    tmp2 = copy(W0.data)
    p = (agg.core, agg.tools, agg.operators, W0, W0_bath, eltype(W0))

    t = 0.01

end

max_stars_repo_name: detrin/OpenQuantumSystems

max_issues_repo_name: detrin/OpenQuantumSystems

max_forks_repo_name: detrin/OpenQuantumSystems


text:
module JungleHelperTreasureChest

using ..Ahorn, Maple

@mapdef Entity "JungleHelper/TreasureChest" TreasureChest(x::Integer, y::Integer, sprite::String="")

const placements = Ahorn.PlacementDict(
    "Treasure Chest (Jungle Helper)" => Ahorn.EntityPlacement(
        TreasureChest
    )
)

function Ahorn.selection(entity::TreasureChest)
    x, y = Ahorn.position(entity)
    return Ahorn.Rectangle(x - 32, y - 24, 64, 48)
end

function Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::TreasureChest, room::Maple.Room)
    Ahorn.drawSprite(ctx, "JungleHelper/Treasure/TreasureIdle00", 0, 0)
end

end

max_stars_repo_name: Jackal-Celeste/JungleHelper

max_issues_repo_name: Jackal-Celeste/JungleHelper

max_forks_repo_name: Jackal-Celeste/JungleHelper


text:
using MRIReco
using Test
using HTTP

@testset "ISMRMRD" begin

filename = joinpath(datadir, "simple_gre.h5")

f = ISMRMRDFile(filename)
acq = RawAcquisitionData(f)

@test acq.profiles[1].head.version == Int16(0)
@test acq.profiles[1].head.measurement_uid == Int32(37)
@test acq.profiles[1].head.scan_counter == Int32(1)
@test acq.profiles[1].head.number_of_samples == Int16(256)
@test acq.profiles[1].head.available_channels == Int16(32)
@test acq.profiles[1].head.active_channels == Int16(32)

params = Dict{Symbol, Any}()
params[:reco] = "direct"
params[:reconSize] = (256,128) #this should be clear from context

Ireco = reconstruction(AcquisitionData(acq), params)
exportImage(joinpath(tmpdir,"recogre.png"), abs.(Ireco[:,:,1,1,1]))

filenameCopy = joinpath(tmpdir, "simple_gre_copy.h5")
fCopy = ISMRMRDFile(filenameCopy)
# store data in another ISMRMRD file
save(fCopy, acq)
acqCopy = RawAcquisitionData(fCopy)

@test convert(AcquisitionHeaderImmutable, acq.profiles[1].head) ==
       convert(AcquisitionHeaderImmutable, acqCopy.profiles[1].head)
@test acqCopy.profiles[1].data == acq.profiles[1].data

IrecoCopy = reconstruction(AcquisitionData(acqCopy), params)

@test IrecoCopy == Ireco

filename = joinpath(datadir, "simple_spiral.h5")

f = ISMRMRDFile(filename)

params = Dict{Symbol, Any}()
params[:reco] = "direct"
params[:reconSize] = (128,128) #this should be clear from context

Ireco = reconstruction(AcquisitionData(f), params)
exportImage(joinpath(tmpdir,"recospiral.png"), abs.(Ireco[:,:,1,1,1]))

filenameCopy = joinpath(tmpdir, "simple_spiral_copy.h5")
fCopy = ISMRMRDFile(filenameCopy)
# store data in another ISMRMRD file
acq = RawAcquisitionData(f)
save(fCopy, acq)
acqCopy = RawAcquisitionData(f)

@test convert(MRIReco.AcquisitionHeaderImmutable, acq.profiles[1].head) ==
       convert(MRIReco.AcquisitionHeaderImmutable, acqCopy.profiles[1].head)
@test acqCopy.profiles[1].traj == acq.profiles[1].traj
@test acqCopy.profiles[1].data == acq.profiles[1].data

# 3D_partial_fourier.h5

#=filename = "3D_partial_fourier.h5"
if !isfile(filename)
  HTTP.open("GET", "https://netix.dl.sourceforge.net/project/ismrmrd/data/3D_partial_fourier.h5") do http
    open(filename, "w") do file
      write(file, http)
    end
  end
end

f = ISMRMRDFile(filename)
=#

#@test size(f.data) == (909, 32, 160)

end

#=

For our own reference: Here is code that can be used with Python to load an
ISMRMRD file

import h5py
f = h5py.File('simple_gre.h5', 'r')
dset = f['/dataset/data']
dset[0]['head']

=#


text:
#-----------------------------------------------------------------------------------------------------
# Spin Matrices
#-----------------------------------------------------------------------------------------------------
# Dictionary
function spin_dict(D::Integer)
    J = (D-1)/2
    coeff = [sqrt(i*(D-i)) for i = 1:D-1]
    sp = sparse(1:D-1, 2:D, coeff, D, D)
    sm = sp'
    sx = (sp+sm)/2
    sy = (sp-sm)/2
    sz = sparse(1:D, 1:D, J:-1:-J)
    s0 = sparse(1.0I, D, D)
    Dict([('+', sp), ('-', sm), ('x', sx), ('y', sy), ('z', sz), ('1', s0)])
end
#-----------------------------------------------------------------------------------------------------
# Atomic spin matrix
function spin_atom(s::String, dic::Dict)
    n = length(s)
    ny = sum(si == 'y' for si in s)
    temp = n == 1 ? dic[s[1]] : kron([dic[si] for si in s]...)
    P = mod(ny, 2) == 0 ? (-1)^(ny÷2) : (1im)^ny
    P * temp
end
#-----------------------------------------------------------------------------------------------------
export spin
"""
    spin(spins; D=2)

Return matrix for spin operators.
The spins should be an iterable onject, each item is of the form (::Number, ::String).
"""
function spin(spins...; D::Integer=2)
    dic = spin_dict(D)
    res = sum(ci * spin_atom(si, dic) for (ci, si) in spins)
    Array(res)
end

spin(spins::AbstractVector{<:Tuple{<:Number, String}}; D::Integer=2) = spin(spins..., D=D)


text:
# ================ INITIALIZATION CODE DO NOT EDIT ====================================== %
# some global parameters -
BIG = 1e10
SMALL = 1e-8

# Load the experiment specifications -
tmp_value = JSON.parsefile("./experiments/Experiments.json")
experiment_array = tmp_value["experiment_array"]

# preload the data -
cached_data_dictionary = Dict()
for (experiment_index,experiment_dictionary) in enumerate(experiment_array)

  # Grab the data file -
  data_file_path = experiment_dictionary["data_file"]
  experiment_id = experiment_dictionary["experiment_id"]

  # Load the data -
  experimental_data_array = readdlm(data_file_path)

  # Cache the data w/experiment_id -
  cached_data_dictionary[experiment_id] = experimental_data_array
end

# preload the error function calls -
cached_error_function_dictionary = Dict()
for (experiment_index,experiment_dictionary) in enumerate(experiment_array)

  # Grab the experimental id -
  experiment_id = experiment_dictionary["experiment_id"]

  # Grab the error functions from the experiment_dictionary -
  error_function = eval(parse(experiment_dictionary["error_function"]))

  # cache these functions -
  cached_error_function_dictionary[experiment_id] = error_function
end
# ================ INITIALIZATION CODE DO NOT EDIT ====================================== %

function objective_function(parameter_array)

  # Script to solve the balance equations -
  time_start = 0.0
  time_stop = 120.0
  time_step_size = 0.01

  # Load the data dictionary -
  data_dictionary = DataDictionary(time_start,time_stop,time_step_size)

  # Update data dictionary to match new parameters before calculating obj
  parameter_mapping_array = data_dictionary["parameter_name_mapping_array"]
  for index = 1:length(parameter_mapping_array)
    if index <= data_dictionary["number_of_binding"]
      data_dictionary["binding_parameter_dictionary"][parameter_mapping_array[index]] = parameter_array[index]
    else
      data_dictionary["control_parameter_dictionary"][parameter_mapping_array[index]] = parameter_array[index]
    end
  end

  # how many objectives do we have?
  obj_array = BIG*ones(13,1)

  # Call simulation routine -
  # (run the model to SS, and then set the ICs to the SS for this parameter set)
  (time_array,simulation_state_array) = Simulation(time_start,time_stop,time_step_size,data_dictionary);

  # Call the error functions -
  # loop through the experimental dictionary, and call the appropriate error function -
  for (experiment_index,experiment_dictionary) in enumerate(experiment_array)

    # Get the error function pointer -
    error_function_pointer = experiment_dictionary["error_function"]
    output_index = parse(Int,experiment_dictionary["output_index"])
    experiment_id = experiment_dictionary["experiment_id"]
    species_symbol = experiment_dictionary["protein_symbol"];

    # Get the experimental data array -
    experimental_data_array = cached_data_dictionary[experiment_id]

    # Call the error function -
    error_function = cached_error_function_dictionary[experiment_id]
    error_value = error_function(experimental_data_array,time_array,simulation_state_array,output_index,species_symbol,data_dictionary)

    # Add the error to the objective array -
    obj_array[experiment_index] = error_value;
  end

  # @show obj_array

  # return the error array -
  return obj_array;
end

function neighbor_function(parameter_array)

  SIGMA = 0.01
  number_of_parameters = length(parameter_array)

  # calculate new parameters -
  perturbed_parameter_array = parameter_array.*(1+SIGMA*randn(number_of_parameters));

  # Parameters -
  # "n_gene_AP1_gene_AhR"	;	# 1
  # "K_gene_AP1_gene_AhR"	;	# 2
  # "n_gene_AP1_gene_PU1"	;	# 3
  # "K_gene_AP1_gene_PU1"	;	# 4
  # "n_gene_AP1_gene_PPARg"	;	# 5
  # "K_gene_AP1_gene_PPARg"	;	# 6
  # "n_gene_AhR_gene_Trigger"	;	# 7
  # "K_gene_AhR_gene_Trigger"	;	# 8
  # "n_gene_CD11b_gene_PU1_gene_cRAF"	;	# 9
  # "K_gene_CD11b_gene_PU1_gene_cRAF"	;	# 10
  # "n_gene_CD14_gene_PPARg_gene_CEBPa_gene_EGR1_gene_cRAF"	;	# 11
  # "K_gene_CD14_gene_PPARg_gene_CEBPa_gene_EGR1_gene_cRAF"	;	# 12
  # "n_gene_CD38_gene_IRF1_gene_PPARg_gene_Trigger_gene_cRAF"	;	# 13
  # "K_gene_CD38_gene_IRF1_gene_PPARg_gene_Trigger_gene_cRAF"	;	# 14
  # "n_gene_CEBPa_gene_Trigger"	;	# 15
  # "K_gene_CEBPa_gene_Trigger"	;	# 16
  # "n_gene_CEBPa_gene_PPARg"	;	# 17
  # "K_gene_CEBPa_gene_PPARg"	;	# 18
  # "n_gene_CEBPa_gene_CEBPa"	;	# 19
  # "K_gene_CEBPa_gene_CEBPa"	;	# 20
  # "n_gene_CEBPa_gene_GFI1"	;	# 21
  # "K_gene_CEBPa_gene_GFI1"	;	# 22
  # "n_gene_E2F_gene_E2F"	;	# 23
  # "K_gene_E2F_gene_E2F"	;	# 24
  # "n_gene_E2F_gene_PPARg"	;	# 25
  # "K_gene_E2F_gene_PPARg"	;	# 26
  # "n_gene_E2F_gene_CEBPa"	;	# 27
  # "K_gene_E2F_gene_CEBPa"	;	# 28
  # "n_gene_E2F_gene_GFI1"	;	# 29
  # "K_gene_E2F_gene_GFI1"	;	# 30
  # "n_gene_E2F_gene_cRAF"	;	# 31
  # "K_gene_E2F_gene_cRAF"	;	# 32
  # "n_gene_EGR1_gene_Trigger"	;	# 33
  # "K_gene_EGR1_gene_Trigger"	;	# 34
  # "n_gene_EGR1_gene_PU1"	;	# 35
  # "K_gene_EGR1_gene_PU1"	;	# 36
  # "n_gene_EGR1_gene_PPARg"	;	# 37
  # "K_gene_EGR1_gene_PPARg"	;	# 38
  # "n_gene_EGR1_gene_GFI1"	;	# 39
  # "K_gene_EGR1_gene_GFI1"	;	# 40
  # "n_gene_GFI1_gene_CEBPa"	;	# 41
  # "K_gene_GFI1_gene_CEBPa"	;	# 42
  # "n_gene_GFI1_gene_EGR1"	;	# 43
  # "K_gene_GFI1_gene_EGR1"	;	# 44
  # "n_gene_IRF1_gene_Trigger"	;	# 45
  # "K_gene_IRF1_gene_Trigger"	;	# 46
  # "n_gene_IRF1_gene_AhR"	;	# 47
  # "K_gene_IRF1_gene_AhR"	;	# 48
  # "n_gene_IRF1_gene_PPARg"	;	# 49
  # "K_gene_IRF1_gene_PPARg"	;	# 50
  # "n_gene_OCT1_gene_PPARg"	;	# 51
  # "K_gene_OCT1_gene_PPARg"	;	# 52
  # "n_gene_OCT4_gene_Trigger"	;	# 53
  # "K_gene_OCT4_gene_Trigger"	;	# 54
  # "n_gene_OCT4_gene_AhR"	;	# 55
  # "K_gene_OCT4_gene_AhR"	;	# 56
  # "n_gene_OCT4_gene_cRAF"	;	# 57
  # "K_gene_OCT4_gene_cRAF"	;	# 58
  # "n_gene_P21_gene_Trigger_gene_AP1_gene_PPARg_gene_PU1_gene_IRF1_gene_CEBPa_gene_cRAF"	;	# 59
  # "K_gene_P21_gene_Trigger_gene_AP1_gene_PPARg_gene_PU1_gene_IRF1_gene_CEBPa_gene_cRAF"	;	# 60
  # "n_gene_P21_gene_GFI1"	;	# 61
  # "K_gene_P21_gene_GFI1"	;	# 62
  # "n_gene_P47Phox_gene_PU1_gene_CEBPa_gene_cRAF"	;	# 63
  # "K_gene_P47Phox_gene_PU1_gene_CEBPa_gene_cRAF"	;	# 64
  # "n_gene_P47Phox_gene_PPARg"	;	# 65
  # "K_gene_P47Phox_gene_PPARg"	;	# 66
  # "n_gene_PPARg_gene_Trigger"	;	# 67
  # "K_gene_PPARg_gene_Trigger"	;	# 68
  # "n_gene_PPARg_gene_CEBPa"	;	# 69
  # "K_gene_PPARg_gene_CEBPa"	;	# 70
  # "n_gene_PPARg_gene_EGR1"	;	# 71
  # "K_gene_PPARg_gene_EGR1"	;	# 72
  # "n_gene_PPARg_gene_PU1"	;	# 73
  # "K_gene_PPARg_gene_PU1"	;	# 74
  # "n_gene_PPARg_gene_AP1"	;	# 75
  # "K_gene_PPARg_gene_AP1"	;	# 76
  # "n_gene_PU1_gene_Trigger"	;	# 77
  # "K_gene_PU1_gene_Trigger"	;	# 78
  # "n_gene_PU1_gene_CEBPa"	;	# 79
  # "K_gene_PU1_gene_CEBPa"	;	# 80
  # "n_gene_PU1_gene_PU1"	;	# 81
  # "K_gene_PU1_gene_PU1"	;	# 82
  # "n_gene_PU1_gene_AP1"	;	# 83
  # "K_gene_PU1_gene_AP1"	;	# 84
  # "n_gene_PU1_gene_OCT1"	;	# 85
  # "K_gene_PU1_gene_OCT1"	;	# 86
  # "n_gene_PU1_gene_AhR"	;	# 87
  # "K_gene_PU1_gene_AhR"	;	# 88
  # "n_gene_PU1_gene_GFI1"	;	# 89
  # "K_gene_PU1_gene_GFI1"	;	# 90
  # "W_gene_AP1_RNAP"	;	# 91
  # "W_gene_AP1_gene_AhR"	;	# 92
  # "W_gene_AP1_gene_PU1"	;	# 93
  # "W_gene_AP1_gene_PPARg"	;	# 94
  # "W_gene_AhR_RNAP"	;	# 95
  # "W_gene_AhR_gene_Trigger"	;	# 96
  # "W_gene_CD11b_RNAP"	;	# 97
  # "W_gene_CD11b_gene_PU1_gene_cRAF"	;	# 98
  # "W_gene_CD14_RNAP"	;	# 99
  # "W_gene_CD14_gene_PPARg_gene_CEBPa_gene_EGR1_gene_cRAF"	;	# 100
  # "W_gene_CD38_RNAP"	;	# 101
  # "W_gene_CD38_gene_IRF1_gene_PPARg_gene_Trigger_gene_cRAF"	;	# 102
  # "W_gene_CEBPa_RNAP"	;	# 103
  # "W_gene_CEBPa_gene_Trigger"	;	# 104
  # "W_gene_CEBPa_gene_PPARg"	;	# 105
  # "W_gene_CEBPa_gene_CEBPa"	;	# 106
  # "W_gene_CEBPa_gene_GFI1"	;	# 107
  # "W_gene_E2F_RNAP"	;	# 108
  # "W_gene_E2F_gene_E2F"	;	# 109
  # "W_gene_E2F_gene_PPARg"	;	# 110
  # "W_gene_E2F_gene_CEBPa"	;	# 111
  # "W_gene_E2F_gene_GFI1"	;	# 112
  # "W_gene_E2F_gene_cRAF"	;	# 113
  # "W_gene_EGR1_RNAP"	;	# 114
  # "W_gene_EGR1_gene_Trigger"	;	# 115
  # "W_gene_EGR1_gene_PU1"	;	# 116
  # "W_gene_EGR1_gene_PPARg"	;	# 117
  # "W_gene_EGR1_gene_GFI1"	;	# 118
  # "W_gene_GFI1_RNAP"	;	# 119
  # "W_gene_GFI1_gene_CEBPa"	;	# 120
  # "W_gene_GFI1_gene_EGR1"	;	# 121
  # "W_gene_IRF1_RNAP"	;	# 122
  # "W_gene_IRF1_gene_Trigger"	;	# 123
  # "W_gene_IRF1_gene_AhR"	;	# 124
  # "W_gene_IRF1_gene_PPARg"	;	# 125
  # "W_gene_OCT1_RNAP"	;	# 126
  # "W_gene_OCT1_gene_PPARg"	;	# 127
  # "W_gene_OCT4_RNAP"	;	# 128
  # "W_gene_OCT4_gene_Trigger"	;	# 129
  # "W_gene_OCT4_gene_AhR"	;	# 130
  # "W_gene_OCT4_gene_cRAF"	;	# 131
  # "W_gene_P21_RNAP"	;	# 132
  # "W_gene_P21_gene_Trigger_gene_AP1_gene_PPARg_gene_PU1_gene_IRF1_gene_CEBPa_gene_cRAF"	;	# 133
  # "W_gene_P21_gene_GFI1"	;	# 134
  # "W_gene_P47Phox_RNAP"	;	# 135
  # "W_gene_P47Phox_gene_PU1_gene_CEBPa_gene_cRAF"	;	# 136
  # "W_gene_P47Phox_gene_PPARg"	;	# 137
  # "W_gene_PPARg_RNAP"	;	# 138
  # "W_gene_PPARg_gene_Trigger"	;	# 139
  # "W_gene_PPARg_gene_CEBPa"	;	# 140
  # "W_gene_PPARg_gene_EGR1"	;	# 141
  # "W_gene_PPARg_gene_PU1"	;	# 142
  # "W_gene_PPARg_gene_AP1"	;	# 143
  # "W_gene_PU1_RNAP"	;	# 144
  # "W_gene_PU1_gene_Trigger"	;	# 145
  # "W_gene_PU1_gene_CEBPa"	;	# 146
  # "W_gene_PU1_gene_PU1"	;	# 147
  # "W_gene_PU1_gene_AP1"	;	# 148
  # "W_gene_PU1_gene_OCT1"	;	# 149
  # "W_gene_PU1_gene_AhR"	;	# 150
  # "W_gene_PU1_gene_GFI1"	;	# 151
  # "W_gene_Trigger_RNAP"	;	# 152
  # "W_gene_cRAF_RNAP"	;	# 153

  # Setup my upper bound, and lower bounds on parameters -
  lower_bound_array = zeros(number_of_parameters)
  upper_bound_array = [
    4.0   ;   # "n_gene_AP1_gene_AhR"	;	# 1
    500.0 ;   # "K_gene_AP1_gene_AhR"	;	# 2
    4.0   ;   # "n_gene_AP1_gene_PU1"	;	# 3
    500.0 ;   # "K_gene_AP1_gene_PU1"	;	# 4
    4.0   ;   # "n_gene_AP1_gene_PPARg"	;	# 5
    500.0 ;   # "K_gene_AP1_gene_PPARg"	;	# 6
    4.0   ;   # "n_gene_AhR_gene_Trigger"	;	# 7
    500.0 ;   # "K_gene_AhR_gene_Trigger"	;	# 8
    4.0   ;  # "n_gene_CD11b_gene_PU1_gene_cRAF"	;	# 9
    500.0 ;  # "K_gene_CD11b_gene_PU1_gene_cRAF"	;	# 10
    4.0   ;  # "n_gene_CD14_gene_PPARg_gene_CEBPa_gene_EGR1_gene_cRAF"	;	# 11
    500.0 ;  # "K_gene_CD14_gene_PPARg_gene_CEBPa_gene_EGR1_gene_cRAF"	;	# 12
    4.0   ;  # "n_gene_CD38_gene_IRF1_gene_PPARg_gene_Trigger_gene_cRAF"	;	# 13
    500.0 ;  # "K_gene_CD38_gene_IRF1_gene_PPARg_gene_Trigger_gene_cRAF"	;	# 14
    4.0   ;  # "n_gene_CEBPa_gene_Trigger"	;	# 15
    500.0 ;  # "K_gene_CEBPa_gene_Trigger"	;	# 16
    4.0   ;  # "n_gene_CEBPa_gene_PPARg"	;	# 17
    500.0 ;  # "K_gene_CEBPa_gene_PPARg"	;	# 18
    4.0   ;  # "n_gene_CEBPa_gene_CEBPa"	;	# 19
    500.0 ;  # "K_gene_CEBPa_gene_CEBPa"	;	# 20
    4.0   ;  # "n_gene_CEBPa_gene_GFI1"	;	# 21
    500.0 ;  # "K_gene_CEBPa_gene_GFI1"	;	# 22
    4.0   ;  # "n_gene_E2F_gene_E2F"	;	# 23
    500.0 ;  # "K_gene_E2F_gene_E2F"	;	# 24
    4.0   ;  # "n_gene_E2F_gene_PPARg"	;	# 25
    500.0 ;  # "K_gene_E2F_gene_PPARg"	;	# 26
    4.0   ;  # "n_gene_E2F_gene_CEBPa"	;	# 27
    500.0 ;  # "K_gene_E2F_gene_CEBPa"	;	# 28
    4.0   ;  # "n_gene_E2F_gene_GFI1"	;	# 29
    500.0 ;  # "K_gene_E2F_gene_GFI1"	;	# 30
    4.0   ;  # "n_gene_E2F_gene_cRAF"	;	# 31
    500.0 ;  # "K_gene_E2F_gene_cRAF"	;	# 32
    4.0   ;  # "n_gene_EGR1_gene_Trigger"	;	# 33
    500.0 ;  # "K_gene_EGR1_gene_Trigger"	;	# 34
    4.0   ;  # "n_gene_EGR1_gene_PU1"	;	# 35
    500.0 ;  # "K_gene_EGR1_gene_PU1"	;	# 36
    4.0   ;  # "n_gene_EGR1_gene_PPARg"	;	# 37
    500.0 ;  # "K_gene_EGR1_gene_PPARg"	;	# 38
    4.0   ;  # "n_gene_EGR1_gene_GFI1"	;	# 39
    500.0 ;  # "K_gene_EGR1_gene_GFI1"	;	# 40
    4.0   ;  # "n_gene_GFI1_gene_CEBPa"	;	# 41
    500.0 ;  # "K_gene_GFI1_gene_CEBPa"	;	# 42
    4.0   ;  # "n_gene_GFI1_gene_EGR1"	;	# 43
    500.0 ;  # "K_gene_GFI1_gene_EGR1"	;	# 44
    4.0   ;  # "n_gene_IRF1_gene_Trigger"	;	# 45
    500.0 ;  # "K_gene_IRF1_gene_Trigger"	;	# 46
    4.0   ;  # "n_gene_IRF1_gene_AhR"	;	# 47
    500.0 ;  # "K_gene_IRF1_gene_AhR"	;	# 48
    4.0   ;  # "n_gene_IRF1_gene_PPARg"	;	# 49
    500.0 ;  # "K_gene_IRF1_gene_PPARg"	;	# 50
    4.0   ;  # "n_gene_OCT1_gene_PPARg"	;	# 51
    500.0 ;  # "K_gene_OCT1_gene_PPARg"	;	# 52
    4.0   ;  # "n_gene_OCT4_gene_Trigger"	;	# 53
    500.0 ;  # "K_gene_OCT4_gene_Trigger"	;	# 54
    4.0   ;  # "n_gene_OCT4_gene_AhR"	;	# 55
    500.0 ;  # "K_gene_OCT4_gene_AhR"	;	# 56
    4.0   ;  # "n_gene_OCT4_gene_cRAF"	;	# 57
    500.0 ;  # "K_gene_OCT4_gene_cRAF"	;	# 58
    4.0   ;  # "n_gene_P21_gene_Trigger_gene_AP1_gene_PPARg_gene_PU1_gene_IRF1_gene_CEBPa_gene_cRAF"	;	# 59
    500.0 ;  # "K_gene_P21_gene_Trigger_gene_AP1_gene_PPARg_gene_PU1_gene_IRF1_gene_CEBPa_gene_cRAF"	;	# 60
    4.0   ;  # "n_gene_P21_gene_GFI1"	;	# 61
    500.0 ;  # "K_gene_P21_gene_GFI1"	;	# 62
    4.0   ;  # "n_gene_P47Phox_gene_PU1_gene_CEBPa_gene_cRAF"	;	# 63
    500.0 ;  # "K_gene_P47Phox_gene_PU1_gene_CEBPa_gene_cRAF"	;	# 64
    4.0   ;  # "n_gene_P47Phox_gene_PPARg"	;	# 65
    500.0 ;  # "K_gene_P47Phox_gene_PPARg"	;	# 66
    4.0   ;  # "n_gene_PPARg_gene_Trigger"	;	# 67
    500.0 ;  # "K_gene_PPARg_gene_Trigger"	;	# 68
    4.0   ;  # "n_gene_PPARg_gene_CEBPa"	;	# 69
    500.0 ;  # "K_gene_PPARg_gene_CEBPa"	;	# 70
    4.0   ;  # "n_gene_PPARg_gene_EGR1"	;	# 71
    500.0 ;  # "K_gene_PPARg_gene_EGR1"	;	# 72
    4.0   ;  # "n_gene_PPARg_gene_PU1"	;	# 73
    500.0 ;  # "K_gene_PPARg_gene_PU1"	;	# 74
    4.0   ;  # "n_gene_PPARg_gene_AP1"	;	# 75
    500.0 ;  # "K_gene_PPARg_gene_AP1"	;	# 76
    4.0   ;  # "n_gene_PU1_gene_Trigger"	;	# 77
    500.0 ;  # "K_gene_PU1_gene_Trigger"	;	# 78
    4.0   ;  # "n_gene_PU1_gene_CEBPa"	;	# 79
    500.0 ;  # "K_gene_PU1_gene_CEBPa"	;	# 80
    4.0   ;  # "n_gene_PU1_gene_PU1"	;	# 81
    500.0 ;  # "K_gene_PU1_gene_PU1"	;	# 82
    4.0   ;  # "n_gene_PU1_gene_AP1"	;	# 83
    500.0 ;  # "K_gene_PU1_gene_AP1"	;	# 84
    4.0   ;  # "n_gene_PU1_gene_OCT1"	;	# 85
    500.0 ;  # "K_gene_PU1_gene_OCT1"	;	# 86
    4.0   ;  # "n_gene_PU1_gene_AhR"	;	# 87
    500.0 ;  # "K_gene_PU1_gene_AhR"	;	# 88
    4.0   ;  # "n_gene_PU1_gene_GFI1"	;	# 89
    500.0 ;  # "K_gene_PU1_gene_GFI1"	;	# 90
    100.0 ; # "W_gene_AP1_RNAP"	;	# 91
    100.0 ;  # "W_gene_AP1_gene_AhR"	;	# 92
    100.0 ;  # "W_gene_AP1_gene_PU1"	;	# 93
    100.0 ;  # "W_gene_AP1_gene_PPARg"	;	# 94
    100.0 ;  # "W_gene_AhR_RNAP"	;	# 95
    100.0 ;  # "W_gene_AhR_gene_Trigger"	;	# 96
    100.0 ;  # "W_gene_CD11b_RNAP"	;	# 97
    100.0 ;  # "W_gene_CD11b_gene_PU1_gene_cRAF"	;	# 98
    100.0 ;  # "W_gene_CD14_RNAP"	;	# 99
    0.0   ;  # "W_gene_CD14_gene_PPARg_gene_CEBPa_gene_EGR1_gene_cRAF"	;	# 100
    100.0 ;  # "W_gene_CD38_RNAP"	;	# 101
    100.0 ;  # "W_gene_CD38_gene_IRF1_gene_PPARg_gene_Trigger_gene_cRAF"	;	# 102
    100.0 ;  # "W_gene_CEBPa_RNAP"	;	# 103
    100.0 ;  # "W_gene_CEBPa_gene_Trigger"	;	# 104
    100.0 ;  # "W_gene_CEBPa_gene_PPARg"	;	# 105
    100.0 ;  # "W_gene_CEBPa_gene_CEBPa"	;	# 106
    100.0 ;  # "W_gene_CEBPa_gene_GFI1"	;	# 107
    100.0 ;  # "W_gene_E2F_RNAP"	;	# 108
    100.0 ;  # "W_gene_E2F_gene_E2F"	;	# 109
    100.0 ;  # "W_gene_E2F_gene_PPARg"	;	# 110
    100.0 ;  # "W_gene_E2F_gene_CEBPa"	;	# 111
    100.0 ;  # "W_gene_E2F_gene_GFI1"	;	# 112
    100.0 ;  # "W_gene_E2F_gene_cRAF"	;	# 113
    100.0 ;  # "W_gene_EGR1_RNAP"	;	# 114
    100.0 ;  # "W_gene_EGR1_gene_Trigger"	;	# 115
    100.0 ;  # "W_gene_EGR1_gene_PU1"	;	# 116
    100.0 ;  # "W_gene_EGR1_gene_PPARg"	;	# 117
    100.0 ;  # "W_gene_EGR1_gene_GFI1"	;	# 118
    100.0 ;  # "W_gene_GFI1_RNAP"	;	# 119
    100.0 ;  # "W_gene_GFI1_gene_CEBPa"	;	# 120
    100.0 ;  # "W_gene_GFI1_gene_EGR1"	;	# 121
    100.0 ;  # "W_gene_IRF1_RNAP"	;	# 122
    100.0 ;  # "W_gene_IRF1_gene_Trigger"	;	# 123
    100.0 ;  # "W_gene_IRF1_gene_AhR"	;	# 124
    100.0 ;  # "W_gene_IRF1_gene_PPARg"	;	# 125
    100.0 ;  # "W_gene_OCT1_RNAP"	;	# 126
    100.0 ;  # "W_gene_OCT1_gene_PPARg"	;	# 127
    100.0 ;  # "W_gene_OCT4_RNAP"	;	# 128
    100.0 ;  # "W_gene_OCT4_gene_Trigger"	;	# 129
    100.0 ;  # "W_gene_OCT4_gene_AhR"	;	# 130
    100.0 ;  # "W_gene_OCT4_gene_cRAF"	;	# 131
    100.0 ;  # "W_gene_P21_RNAP"	;	# 132
    100.0 ;  # "W_gene_P21_gene_Trigger_gene_AP1_gene_PPARg_gene_PU1_gene_IRF1_gene_CEBPa_gene_cRAF"	;	# 133
    100.0 ;  # "W_gene_P21_gene_GFI1"	;	# 134
    100.0 ;  # "W_gene_P47Phox_RNAP"	;	# 135
    100.0 ;  # "W_gene_P47Phox_gene_PU1_gene_CEBPa_gene_cRAF"	;	# 136
    100.0 ;  # "W_gene_P47Phox_gene_PPARg"	;	# 137
    100.0 ;  # "W_gene_PPARg_RNAP"	;	# 138
    100.0 ;  # "W_gene_PPARg_gene_Trigger"	;	# 139
    100.0 ;  # "W_gene_PPARg_gene_CEBPa"	;	# 140
    100.0 ;  # "W_gene_PPARg_gene_EGR1"	;	# 141
    100.0 ;  # "W_gene_PPARg_gene_PU1"	;	# 142
    100.0 ;  # "W_gene_PPARg_gene_AP1"	;	# 143
    100.0 ;  # "W_gene_PU1_RNAP"	;	# 144
    100.0 ;  # "W_gene_PU1_gene_Trigger"	;	# 145
    100.0 ;  # "W_gene_PU1_gene_CEBPa"	;	# 146
    100.0 ;  # "W_gene_PU1_gene_PU1"	;	# 147
    100.0 ;  # "W_gene_PU1_gene_AP1"	;	# 148
    100.0 ;  # "W_gene_PU1_gene_OCT1"	;	# 149
    100.0 ;  # "W_gene_PU1_gene_AhR"	;	# 150
    100.0 ;  # "W_gene_PU1_gene_GFI1"	;	# 151
    100.0 ;  # "W_gene_Trigger_RNAP"	;	# 152
    100.0 ;  # "W_gene_cRAF_RNAP"	;	# 153
  ];

  # return the corrected value -
  return parameter_bounds_function(perturbed_parameter_array,lower_bound_array,upper_bound_array);
end

function acceptance_probability_function(rank_array,temperature)
  return (exp(-rank_array[end]/temperature))
end

function cooling_function(temperature)

  # define my new temperature -
  alpha = 0.90
  return alpha*temperature
end

function parameter_bounds_function(parameter_array,lower_bound_array,upper_bound_array)

  # iterate through and fix the parameters -
  new_parameter_array = copy(parameter_array)
  for (index,value) in enumerate(parameter_array)

    lower_bound = lower_bound_array[index]
    upper_bound = upper_bound_array[index]

    if (value<lower_bound)
      new_parameter_array[index] = lower_bound
    elseif (value>upper_bound)
      new_parameter_array[index] = upper_bound
    end
  end

  return new_parameter_array
end


text:
using Test

import MathOptInterface
const MOI = MathOptInterface

@testset "hs071" begin
    mock = MOI.Utilities.MockOptimizer(
        MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()),
        eval_objective_value=false
    )
    config = MOI.Test.TestConfig(optimal_status = MOI.LOCALLY_SOLVED)
    MOI.Utilities.set_mock_optimize!(
        mock,
        (mock) -> begin
            MOI.Utilities.mock_optimize!(
                mock, config.optimal_status,
                [1.0, 4.7429996418092970, 3.8211499817883077, 1.379408289755698]
            )
            MOI.set(mock, MOI.ObjectiveValue(), 17.014017145179164)
        end
    )
    MOI.Test.hs071_test(mock, config)
    MOI.Test.hs071_no_hessian_test(mock, config)

    d = MOI.Test.HS071(false)
    VI = MOI.VariableIndex
    @test MOI.objective_expr(d) == :(x[$(VI(1))] * x[$(VI(4))] * (x[$(VI(1))] +
                                     x[$(VI(2))] + x[$(VI(3))]) + x[$(VI(3))])
    @test MOI.constraint_expr(d, 1) ==
        :(x[$(VI(1))] * x[$(VI(2))] * x[$(VI(3))] * x[$(VI(4))] >= 25.0)
    @test MOI.constraint_expr(d, 2) ==
        :(x[$(VI(1))]^2 + x[$(VI(2))]^2 + x[$(VI(3))]^2 + x[$(VI(4))]^2 == 40.0)
    @test_throws ErrorException MOI.constraint_expr(d, 3)
end

@testset "mixed_complementarity" begin
    mock = MOI.Utilities.MockOptimizer(
        MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}())
    )
    config = MOI.Test.TestConfig(optimal_status = MOI.LOCALLY_SOLVED)
    MOI.Utilities.set_mock_optimize!(
        mock,
        (mock) -> MOI.Utilities.mock_optimize!(
            mock, config.optimal_status, [2.8, 0.0, 0.8, 1.2]
        )
    )
    MOI.Test.mixed_complementaritytest(mock, config)
end

@testset "math_program_complementarity_constraints" begin
    mock = MOI.Utilities.MockOptimizer(
        MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}())
    )
    config = MOI.Test.TestConfig(optimal_status = MOI.LOCALLY_SOLVED)
    MOI.Utilities.set_mock_optimize!(
        mock,
        (mock) -> MOI.Utilities.mock_optimize!(
            mock, config.optimal_status, [1.0, 0.0, 3.5, 0.0, 0.0, 0.0, 3.0, 6.0]
        )
    )
    MOI.Test.math_program_complementarity_constraintstest(mock, config)
end


text: Base.promote_rule(::Type{<:Combiner}, StorageT::Type{<:CuDense}) = StorageT


text:
function run_train(M,net_param,train_param,time_param,dt,utarg)

N        = length(M[1,:]);
w        = zeros(N,N);
extinput = zeros(N);
u        = zeros(N);

x, r, b = net_param[1:N], net_param[N+1:2*N], net_param[2*N+1];
stim, lambda, learn_every, nloop, target_type = train_param[1:N], train_param[N+1], train_param[N+2], train_param[N+3], train_param[N+4];
stim_on, stim_off, train_time = time_param[1], time_param[2], time_param[3];

# RHS of ODEs
dtheta(x_var,I_var) = 1 - cos(x_var) + I_var.*(1 + cos(x_var));
dr(r_var) = -b*r_var;

# set up correlation matrix
P = Dict{Int64,Array{Float64,2}}();
Px = Dict{Int64,BitArray{2}}();
for ni=1:Int(N)
    ni_prenum = sum(M[ni,:].!=0);
    ni_preind = M[ni,:].!=0;
    P[ni] = (1.0/lambda)*eye(ni_prenum);
    Px[ni] = ni_preind;
end

t = 0.;
for iloop = 1:Int(nloop)
    print("Loop $(iloop)...",'\n')
    for i=1:Int(train_time/dt)
          # External stimulus
          if i > Int(stim_on/dt) && i <= Int(stim_off/dt)
                extinput[:] = stim[:];
          else
                extinput[:] = zeros(N)[:];
          end

          # Update neuron phase
          k1 = dt*dtheta(x,      u+extinput);
          k2 = dt*dtheta(x+k1/2, u+extinput);
          k3 = dt*dtheta(x+k2/2, u+extinput);
          k4 = dt*dtheta(x+k3,   u+extinput);
          xnext  = x + (k1+2*k2+2*k3+k4)/6;

          # Update filtered spikes
          l1 = dt*dr(r);
          l2 = dt*dr(r+l1/2);
          l3 = dt*dr(r+l2/2);
          l4 = dt*dr(r+l3);

          # Spike detection.
          idx1 = xnext - x .> 0.;
          idx2 = xnext - x .> mod(pi - mod(x,2*pi),2*pi); # distance to next pi.
          idx  = idx1.*idx2;
          ind  = collect(1:N)[idx[:]]';

          # Update each neuron's spike
          if !isempty(ind)
              r[idx[:]] = r[idx[:]] + b;
          end

          # Update x, r, z, t.
          x = x + (k1+2*k2+2*k3+k4)/6;
          r = r + (l1+2*l2+2*l3+l4)/6;
          u = M*r;
          t = t + dt;

          # Train w
          if i > Int(stim_off/dt)+1 && i <= Int(train_time/dt) && mod(i, learn_every/dt) == 0
              for ni=1:Int(N)
                	# update inverse correlation matrix
                	k     = P[ni]*r[Px[ni][:]];
                	vPv   = r[Px[ni][:]].'*k;
                	den   = 1.0/(1.0 + vPv[1]);
                	P[ni] = P[ni] - k*(k.'*den);

                	# update recurrent weights
                  e  = M[ni,Px[ni][:]]*r[Px[ni][:]] - utarg[i,ni];
                	dw = -e[1]*k*den;
                	M[ni,Px[ni][:]] = M[ni,Px[ni][:]] + dw';
              end
          end
    end
end

return M

end


text:
using Test, JLD2, LinearAlgebra
#
# gradient computation by parial automatic differentiation and adjoint merthods
# of
#           E(v) = ∫ |U(v)|² dx
# where
# A(v) U(v) = F(v) the discrete version of - Δ u = f
# and v are the vertices of the voronoi cells
#

# load polygonal mesh
filename, nc = "squarepolmesh_coarse", 100
mesh_filename = "$(@__DIR__)/data/$(filename)_$(nc).jld2"
println("read file $(mesh_filename)")
JLD2.@load(mesh_filename, pv, cellsb, cellsbt, t, pb, tb)
 # boundary points / rhs
 rhs, boundary_condition = if occursin("square", mesh_filename)
    PolygonalFem.rhs_sqr, PolygonalFem.boundary_condition_sqr
else
    PolygonalFem.rhs_L, PolygonalFem.boundary_condition_L
end

np, ε = size(pv, 1), 1e-4
x, dp = pv, 2 * rand(size(pv)...) .- 1
Ib = unique(PolygonalFem.btri(t)[:])
xp, xm = x .+ ε * dp, x .- ε * dp
# finite differences
f = x -> PolygonalFem.costnormU(x, cellsb, t, rhs, boundary_condition, Ib)
fp, fm = f(xp), f(xm)
@time g = PolygonalFem.∇costnormU(x, cellsb, t, rhs, boundary_condition, Ib)
@time g = PolygonalFem.∇costnormU(x, cellsb, t, rhs, boundary_condition, Ib)
println("")
println(fp)
println(fm)
println(" ")
println((fp - fm) / (2 * ε))
println(dot(g, dp))
println("")
@testset "test eval gradient dot function for polygonal costnormU" begin
    @test abs(dot(g, dp) - (fp - fm) / (2 * ε)) < 1e-4
end


text:
using Stan

## testing accessors to CMDSTAN_HOME
let oldpath = Stan.CMDSTAN_HOME
    newpath = Stan.CMDSTAN_HOME * "##test##"
    set_CMDSTAN_HOME!(newpath)
    @test Stan.CMDSTAN_HOME == newpath
    set_CMDSTAN_HOME!(oldpath)
    @test Stan.CMDSTAN_HOME == oldpath
end


text:
using ArgParse
using DrWatson
@quickactivate
using BSON
using Random
using FileIO
using DataFrames
using Base.Threads: @threads
using GenerativeAD

const CHUNK_SIZE = 1000

s = ArgParseSettings()
@add_arg_table! s begin
    "source_prefix"
		arg_type = String
		default = "evaluation/images"
		help = "Data prefix of generated files."
	"filename"
		arg_type = String
		default = "evaluation/images_eval.bson"
		help = "Where to store the cached DataFrame."
	"-f", "--force"
    	action = :store_true
		help = "Overwrite all generated files."
end

"""
	collect_stats(source_prefix::String)

Collects evaluation DataFrames from a given data folder prefix `source_prefix`.
"""
function collect_stats(source_prefix::String)
	source = datadir(source_prefix)
	@info "Collecting files from $source folder."
	files = GenerativeAD.Evaluation.collect_files_th(source)
	@info "Collected $(length(files)) files from $source folder."

	parts = collect(Iterators.partition(files, CHUNK_SIZE))
	frames = Vector{DataFrame}(undef, length(parts))
	@info "Loading into $(length(parts)) partitions."
	@threads for i in 1:length(parts)
		df = reduce(vcat, [load(f)[:df] for f in parts[i]])
		frames[i] = df
		print("-")
	end
	@info "Loaded $(length(parts)) partitions."
	df = reduce(vcat, frames)
end

function main(args)
	df = collect_stats(args["source_prefix"])
	f = datadir(args["filename"])
	if (isfile(f) && args["force"]) || ~isfile(f)
		@info "Saving $(nrow(df)) rows to $f."
		wsave(f, Dict(:df => df))
	end
	@info "---------------- DONE -----------------"
end

main(parse_args(ARGS, s))


text:
"""
    AuroraKNN(;kKNN=1000, loocv=true, tree=KDTree,
               pca=false, pca_dimension=false)

Aurora with `k`-Nearest neighbors.  If `looc=false`, then `k` is chosen equal to `kKNN`,
while if `loocv=true`, then `k` is selected for each held-out replicate by
Leave-One-Out Cross-validation among the choices 1,...,`kKNN`.

`tree` describes the nearest neighbor computation strategy. The following options are
available: `:auto`, as well as ,
`:kdtree`, `:balltree` and `:brutetree` from the `NearestNeighbors.jl` package.

If `pca=true`, a dimension reduction strategy is employed to find nearest neighbors using
PCA (principal component analysis). For each held-out replicate, the order statistics
are projected into the principal component subspace of dimension `pca_dimension`. Note that
in this case, the resulting nearest neighbors may only be interpreted as approximate
nearest neighbors.
"""
Base.@kwdef struct AuroraKNN <: AbstractAurora
    kKNN::Int = 1000
    loocv::Bool = true
    tree::Symbol = :auto
    pca::Bool = false
    pca_dimension::Int = 10
end

Base.@kwdef struct FittedAuroraKNN{Ms, Ss} <: AbstractFittedAurora
    μs_mat::Ms
    μs::Ss
    kskNN::Vector{Int}
    aurora
end

function StatsBase.fit(aurora::AuroraKNN, Zs::AbstractVector{<:ReplicatedSample})
    K = nobs(Zs[1])
    any(nobs.(Zs) != K) || throw("All Zs should have the same number of replicates.")

    n = length(Zs)
    kKNN  = min(aurora.kKNN, n-2)
    loocv = aurora.loocv
    tree_symbol = aurora.tree

    if tree_symbol == :auto
        tree_symbol = K <= 12 ? :kdtree : :brutetree
    end

    if tree_symbol == :kdtree
        tree = KDTree
    elseif tree_symbol == :balltree
        tree = BallTree
    elseif tree_symbol == :brutetree
        tree = BruteTree
    else
        throw("Nearest neighbor strategy chosen is not available.")
    end

    pca = aurora.pca && aurora.pca_dimension < K-1

    μs_mat = zeros(n,K)
    cache_mat = zeros(K-1, n)
    cache_Y = zeros(n)

    if loocv
        cache_knn = zeros(kKNN, n)
        cache_lnn_loo = zeros(kKNN)
        kskNN = zeros(Int, K)
    else
        kskNN = fill(kKNN, K)
    end

    for j in Base.OneTo(K)
        design_matrix!(cache_mat, Zs, j)
        cache_Y[:] = StatsBase.response.(Zs, j)

        if !pca
            kdd = tree(cache_mat)
            _idxs,_dists = knn(kdd, cache_mat, kKNN, loocv);
        else
            pca_fit = fit(PCA, cache_mat; maxoutdim=aurora.pca_dimension, pratio=1)
            pca_projected_X = transform(pca_fit, cache_mat)
            kdd = tree(pca_projected_X)
            _idxs,_dists = knn(kdd, pca_projected_X, kKNN, loocv);
        end

        if loocv
            for i in Base.OneTo(n)
                cumsum!(view(cache_knn, :, i), view(cache_Y, _idxs[i]))
            end
            cache_knn .-= cache_knn[1,:]'
            cache_knn ./= [1; 1:(kKNN-1)]
            cache_knn .-= cache_Y'
            cache_lnn_loo .= vec(mean(abs2, cache_knn; dims=2))
            k_opt =  argmin(cache_lnn_loo)
            kskNN[j] = k_opt
            μs_mat[:,j] .= ((cache_knn[k_opt, :] .+ cache_Y).*(k_opt-1) .+ cache_Y)./ k_opt
        else
            for i in Base.OneTo(n)
                μs_mat[i,j] = mean(view(cache_Y, _idxs[i]))
            end
        end
    end
    μs = mean(μs_mat; dims=2)
    FittedAuroraKNN(μs = μs, μs_mat = μs_mat, aurora=aurora, kskNN = kskNN)
end


text:
# TODO: finish this!
struct RepeatedProductIterator{N, I}
    iterator::I
end

Iterators.product(N::Int, it) = RepeatedProductIterator(N, it)

Base.IteratorSize(::Type{RepeatedProductIterator{N, I}}) where {N, I} = Base.IteratorSize(I)

max_stars_repo_name: Roger-luo/QuHamiltonian

max_issues_repo_name: Roger-luo/QuHamiltonian

max_forks_repo_name: Roger-luo/QuHamiltonian


text:
module RiemannComplexNumbers

import Base.inv, Base.Complex, Base.show
import Base.+, Base.-, Base.*, Base./, Base.==, Base.hash
import Base: isinf, isnan, iszero, Complex, isequal

export RC, ComplexInf, ComplexNaN, IM

struct RC{T<:Complex} <: Number
    val::T
    nan_flag::Bool
    inf_flag::Bool
end

const ComplexNaN = RC((0+0im)/0,true,false)
const ComplexInf = RC(0+0im,false,true)
const IM = RC(im,false,false)

function RC(z::Complex)
    if isnan(z)
        return ComplexNaN
    end
    if isinf(z)
        return ComplexInf
    end
    return RC(z,false,false)
end

function Complex(a::RC)::Complex
    if isinf(a)
        return Inf + Inf*im
    end
    if isnan(a)
        return NaN
    end
    return a.val
end

RC(z::Real) = RC(z+0im)

isinf(z::RC) = z.inf_flag
isnan(z::RC) = z.nan_flag

function iszero(z::RC)
    if isnan(z) || isinf(z)
        return false
    end
    return iszero(z.val)
end

function show(io::IO, z::RC)
    if isinf(z)
        print(io,"ComplexInf")
    elseif isnan(z)
        print(io,"ComplexNaN")
    else
        sz = string(z.val)[1:end-2] * "IM"
        print(io,sz)
    end
end

function (==)(a::RC, b::RC)::Bool
    if isnan(a) || isnan(b)
        return false
    end
    if isinf(a) && isinf(b)
        return true
    end
    if isinf(a) || isinf(b)
        return false
    end
    return a.val == b.val
end

isequal(a::RC,b::Number) = isequal(promote(a,b)...)
isequal(a::Number,b::RC) = isequal(promote(a,b)...)

function isequal(a::RC, b::RC)::Bool
    # for isequal, nan's compare true
    if isnan(a) && isnan(b)
        return true
    end
    # but if only one is nan, then it's false
    if isnan(a) || isnan(b)
        return false
    end
    ## complex infinites are equal
    if isinf(a) && isinf(b)
        return true
    end
    if isinf(a) || isinf(b)
        return false
    end
    # finally, fall back on isequal for Complex
    return isequal(a.val, b.val)
end

function hash(a::RC, h::UInt=UInt(0))
    if isinf(a)
        return hash(Inf,h)
    end
    if isnan(a)
        return hash(NaN,h)
    end
    return hash(a.val,h)
end

import Base.promote_rule

promote_rule(::Type{RC{T}}, ::Type{S}) where {T,S} = RC

include("arithmetic.jl")
include("functions.jl")
end  # end of module


text:
abstract SGDUpdater

immutable MomentumUpdater <: SGDUpdater
    learningRate::Float64
    decayL1::Float64
    decayL2::Float64
    momentum::Float64
    lastUpdate::Vector{Float64}
end

function updateGradients(m::MomentumUpdater, state::GradientState, batchSize::Int)
    rateL1 = m.decayL1 * state.multiplierL1
    rateL2 = m.decayL2 * state.multiplierL2
    for i in 1:length(state.vol.w)
        gradL1 = rateL1 * (state.vol.w[i] > 0 ? 1 : -1)
        gradL2 = rateL2 * (state.vol.w[i] * state.vol.w[i] / 2.0)

        rawGradient = (gradL1 + gradL2 + state.vol.dw[i]) / batchSize
        momentumGradient = m.momentum * m.lastUpdate[i] - m.learningRate * rawGradient

        state.vol.w[i] += momentumGradient
        m.lastUpdate[i] = momentumGradient
    end
    fill!(state.vol.dw, 0.0)
end


text:
using GLVisualize.ComposeBackend, Gadfly
using Colors, GLVisualize

if !isdefined(:runtests)
    window = glscreen()
    composebackend = ComposeBackend.GLVisualizeBackend(window)
end

const static_example = true

using Distributions
x1 = rand(40)
y1 = 4.*x1 .+ 2 .+randn(40)
x2 = rand(40)
y2 = -6.*x2 .+ 3 .+ randn(40)
x  = [x1;x2]
y  = [y1;y2]
col = [fill("Slope 4",40); fill("Slope -6",40)]
p = plot(x=x,y=y,colour=col, Geom.point, Geom.smooth(method=:lm))

draw(composebackend, p)

if !isdefined(:runtests)
	renderloop(window)
end


text:
#!/usr/bin/env julia
################################################################################
#
#  Copyright (c) 2015 Wojciech Migda
#  All rights reserved
#  Distributed under the terms of the MIT license
#
################################################################################
#
#  Filename: rf_clf.jl
#
#  Decription:
#       Classify data using RandomForest.
#
#  Authors:
#       Wojciech Migda
#
################################################################################
#
#  History:
#  --------
#  Date         Who  Ticket     Description
#  ----------   ---  ---------  ------------------------------------------------
#  2015-11-22   wm              Initial version
#
################################################################################

using HDF5
using DecisionTree

function pp(M)
    const MIN = minimum(M, 2)
    const NORM = (M .- MIN) ./ (maximum(M, 2) - MIN)
    const SHIFT = NORM + 0.5 .- mean(NORM, 2)

    return map(x -> max(x, 0.0), map(x -> min(x, 1.0), SHIFT))
end

function main(DATA_PATH::String)

    const X_TRAIN = pp(h5read("$(DATA_PATH)/trainX.h5", "gsv/train/X/data"))
    const Y_TRAIN = h5read("$(DATA_PATH)/trainY.h5", "gsv/train/Y/data")

    println("Training the classifier")
    # Train random forest with
    # 20 for number of features chosen at each random split,
    # 50 for number of trees,
    # and 1.0 for ratio of subsampling.
    model = build_forest(Y_TRAIN, X_TRAIN, 20, 50, 1.0)

    const X_TEST = pp(h5read("$(DATA_PATH)/testX.h5", "gsv/test/X/data"))
    const Y_TEST = char(apply_forest(model, X_TEST))
    println("Prediction done.")

    const TEST_LABELS = h5read("$(DATA_PATH)/testX.h5", "gsv/test/X/labels")

    # Add header option to writedlm(), https://github.com/JuliaLang/julia/issues/10284
    fid = open("$(DATA_PATH)/mySubmission.csv", "w")
    write(fid, "ID,Class\n")
    writecsv(fid, zip(TEST_LABELS, Y_TEST))
    close(fid)

    const ACCURACY = nfoldCV_forest(Y_TRAIN, X_TRAIN, 20, 50, 4, 1.0);
    println("4 fold accuracy: $(mean(ACCURACY))")
end

if ~isinteractive()
    main("/repo/kaggle/street-view-getting-started-with-julia/KAGGLE-first-steps-with-julia/data")
end


text:
using Cat

# "Real" numbers
const R = Float64

# The only objects in Smooth are Nothing, R, and products of these
@category Smooth #where objects = obj -> obj in (Nothing, R, Tuple)

# Basic algebra
@arrow Smooth Plus :: (R, R) --> R
@arrow Smooth Mult :: (R, R) --> R
@arrow Smooth Neg :: R --> R
# Calc 1
@arrow Smooth Exp :: R --> R
@arrow Smooth Sin :: R --> R
@arrow Smooth Cos :: R --> R

# Just an unbound input is equivalent to the identity morphism
Placeholder = Smooth.Identity{R}

@alias Smooth Plus Base.:+ 2
@alias Smooth Neg Base.:- 1
@alias Smooth Mult Base.:* 2

Base.exp(a::Smooth.Arrow) = Exp(a)
Base.sin(a::Smooth.Arrow) = Sin(a)
Base.cos(a::Smooth.Arrow) = Cos(a)

# Evaluate an expression tree
@interpretation Eval (=>) Smooth

@interpret function (e::Eval)(m::Plus, a, b)
    a + b
end

@interpret function (e::Eval)(m::Mult, a, b)
    a * b
end

@interpret function (e::Eval)(m::Neg, a)
    -a
end

@interpret function (e::Eval)(m::Exp, a)
    exp(a)
end

@interpret function (e::Eval)(m::Sin, a)
    sin(a)
end

@interpret function (e::Eval)(m::Cos, a)
    cos(a)
end

#=
Reverse mode diff:

In reverse mode AD, we take a function a -> b
and convert it to a function (a, db) -> (b, da)
where "da" is the space of "differentials of a";
for our basic setup here, where we work in R^N, da == a,
but this could be something like e.g. a tangent/cotangent bundle.
Basically, the RDiff function takes an input and an output differential
and converts them to an output and an input differential.

Mathematically (from a categorical point of view), RDiff is a functor
that takes the Smooth category to a new category where:
(1) the objects are the same as Smooth
and
(2) a morphism (arrow) from A ~> B corresponds to an arrow (A x dB) -> (B x dA) in Smooth
so, we have the same objects/morphisms as smooth, but we are basically changing what composition
means: A ~> B  and B ~> C compose to A ~> C, but the composition rules aren't the same
as Smooth.  This is kind of similar to how Kleisli composition works.

We handle this in _code_ by simply re-using Smooth and changing out RDiff acts on composition.
The downside to this (compared to, say, ForwardDiff) is that we lose some of the automagicness
of letting Cat.jl define how the functor acts on all of the built-in arrows (product, terminal, constant, etc...)
=#
@functor RDiff :: Smooth => Smooth (T -> T)

#= RDiff(m::Smooth.Identity{T}) where {T} = Smooth.Identity{RDiff(T)}() =#
# Cartesian
#= RDiff(m::Smooth.Product) = Smooth.Product([RDiff(x) for x in m.factors]...) =#
# A ~> (B, C), e.g., needs to go to (A, (dB, dC)) -> ((B, C), dA)
RDiff(m::Smooth.Product) = begin
    dum = Smooth.Identity{Tuple{source(m), target(m)}}()
    # Each factor maps to an (A, dB) -> (B, dA)
    factors = []
    for k in 1:length(m.factors)
        x = m.factors[k]
        push!(factors, compose(RDiff(x), Smooth.Product(dum[1], dum[2][k])))
    end
    dA = sum([f[2] for f in factors])
    b = Smooth.Product([f[1] for f in factors]...)
    out = Smooth.Product(b, dA)
    out
end

RDiff(m::Smooth.Proj{A, B}) where {A,B} = begin
    @assert false # this is actually not impelmented properly but we don't need it for this example
end

# Terminal is T ~> Nothing, so (T, Nothing) -> (Nothing, T)
RDiff(m::Smooth.Terminal{T}) where {T <: AbstractFloat} = compose(Smooth.Constant( (nothing, 0.0) ),
                                                                  Smooth.Terminal{Tuple{T, Nothing}}())
# Functions compose differently under reverse mode
function RDiff(m::Smooth.Composed)
    # f::a -> b
    # g::b -> c
    inps = Smooth.Identity{Tuple{source(m.f), target(m.g)}}()
    df = RDiff(m.f) # df :: (a, db) -> (b, da)
    dg = RDiff(m.g) # dg :: (b, dc) -> (c, db)
    # We need to "intertwine" these to get:
    # dg o df :: (a, dc) -> (c, da)
    #=
    z = dg(f(inps[1]), inps[2])
    w = df(inps[1], z[2])
    (w[2], z[1])
    =#
    f = compose(m.f, inps[1]) # A better implementation would memoize f
    z = compose(dg, Smooth.Product(f, inps[2])) # z :: (~) ~> (c, db)
    w = compose(df, Smooth.Product(inps[1], z[2]))
    Smooth.Product(z[1], w[2])
end

function split_inputs(m::Smooth.Arrow)
    a = Smooth.Identity{Tuple{source(m), target(m)}}()
    return (a[1], a[2])
end

function RDiff(m::Plus)
    a, db = split_inputs(m)
    Smooth.Product(a[1] + a[2],
                   Smooth.Product(db, db))
end

function RDiff(m::Mult)
    a, db = split_inputs(m)
    Smooth.Product(a[1]*a[2],
                   Smooth.Product(a[1]*db, a[2]*db))
end

function RDiff(m::Smooth.Constant)
    # map inputs into Nothing, then return a constant
    compose(Smooth.Product(Smooth.Constant(m.val),
                           Smooth.Constant(nothing)),
            Smooth.Terminal{Tuple{Nothing, target(m)}}())
end

function RDiff(m::Exp)
    a, db = split_inputs(m)
    u = exp(a)
    Smooth.Product(u, db*u)
end

function RDiff(m::Sin)
    a, db = split_inputs(m)
    Smooth.Product(sin(a), db*cos(a))
end

function RDiff(m::Cos)
    a, db = split_inputs(m)
    Smooth.Product(cos(a), -db*sin(a))
end

# Here is the actual "user code"
x = Placeholder()
y = cos(sin(x))
y = exp(y*y) + 2.0
dy = RDiff(y)

a = 1.0
ya = Eval()(y, a)
ya2 = Eval()(y, a+1e-3)
println("y($a)=$(ya)")
println("dy($a)=$(Eval()(dy, a, 1.0))")
println("Numerical deriv: $(1e3*(ya2-ya))")


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct VirtualMachineScaleSetVMReimageParameters <: SwaggerModel
    tempDisk::Any # spec type: Union{ Nothing, Bool } # spec name: tempDisk

    function VirtualMachineScaleSetVMReimageParameters(;tempDisk=nothing)
        o = new()
        validate_property(VirtualMachineScaleSetVMReimageParameters, Symbol("tempDisk"), tempDisk)
        setfield!(o, Symbol("tempDisk"), tempDisk)
        o
    end
end # type VirtualMachineScaleSetVMReimageParameters

const _property_map_VirtualMachineScaleSetVMReimageParameters = Dict{Symbol,Symbol}(Symbol("tempDisk")=>Symbol("tempDisk"))
const _property_types_VirtualMachineScaleSetVMReimageParameters = Dict{Symbol,String}(Symbol("tempDisk")=>"Bool")
Base.propertynames(::Type{ VirtualMachineScaleSetVMReimageParameters }) = collect(keys(_property_map_VirtualMachineScaleSetVMReimageParameters))
Swagger.property_type(::Type{ VirtualMachineScaleSetVMReimageParameters }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_VirtualMachineScaleSetVMReimageParameters[name]))}
Swagger.field_name(::Type{ VirtualMachineScaleSetVMReimageParameters }, property_name::Symbol) =  _property_map_VirtualMachineScaleSetVMReimageParameters[property_name]

function check_required(o::VirtualMachineScaleSetVMReimageParameters)
    true
end

function validate_property(::Type{ VirtualMachineScaleSetVMReimageParameters }, name::Symbol, val)
end


text:
module TestInference

using Baselet
using Test

valof(::Val{x}) where {x} = x

macro test_inferred(ex)
    quote
        f() = $(esc(ex))
        Test.@test Test.@inferred(f()) isa Any
    end
end

@testset "getindex" begin
    @test_inferred Baselet.getindex((Val(1), Val(2), Val(3), Val(4)), 2:3)
end

@testset "zip" begin
    @test_inferred Baselet.zip((Val(1), Val(2), Val(3), Val(4)))
    @test_inferred Baselet.zip(
        (Val(1), Val(2), Val(3), Val(4)),
        (Val(5), Val(6), Val(7), Val(8)),
    )
end

@testset "flatten" begin
    @test_inferred Baselet.flatten(((Val(1), Val(2)), (Val(3), Val(4), Val(5))))
end

@testset "enumerate" begin
    if VERSION >= v"1.6-"
        @info "Skip inference test for `enumerate` on Julia $VERSION"
    else
        @test_inferred Val(Baselet.enumerate((Val(1), Val(2), Val(3), Val(4))))
    end
end

@testset "any" begin
    @test_inferred Val(Baselet.any(valof, (Val(false), Val(false), Val(false))))
    @test_inferred Val(Baselet.any(valof, (Val(false), Val(false), Val(true))))
end

@testset "all" begin
    @test_inferred Val(Baselet.all(valof, (Val(true), Val(true), Val(true))))
    @test_inferred Val(Baselet.all(valof, (Val(true), Val(true), Val(false))))
end

@testset "sort" begin
    @test_inferred Baselet.sort((Val(3), Val(1), Val(2), Val(0)), by = valof)
end

end  # module


text:
import DSP.resample
export resample, resample!

function mkresample(rate::T) where T<:Real
  r = rationalize(rate)
  h = T.(resample_filter(r))
  ff = FIRFilter(h, r)
  τ = timedelay(ff)
  setphase!(ff, τ)
  return ff
end

@doc """
    resample!(S::SeisData [, chans=CC, fs=FS])
    resample(S::SeisData [, chans=CC, fs=FS])

Resample data in S to `FS`. If keyword `fs` is not specified, data are resampled
to the lowest non-zero value in `S.fs[CC]`.Note that a poor choice of `FS` can
lead to upsampling and other undesirable behavior.

Use keyword `chans=CC` to only resample channel numbers `CC`. By default, all
channels `i` with `S.fs[i] > 0.0` are resampled.

    resample!(C::SeisChannel, fs::Float64)
    resample(C::SeisChannel, fs::Float64)

Resample `C.x` to `fs`.
""" resample!
function resample!(S::GphysData;
  chans::Union{Integer, UnitRange, Array{Int64,1}}=Int64[],
  fs::Float64=0.0)

  if chans == Int64[]
    chans = 1:S.n
  end

  f0 = fs == 0.0 ? minimum(S.fs[chans[S.fs[chans] .> 0.0]]) : fs

  # This setup is very similar to filtfilt!
  N = nx_max(S)
  T = unique([eltype(i) for i in S.x[chans]])
  nT = length(T)

  sz = 0
  yy = Any
  for i = 1:nT
    zz = sizeof(T[i])
    if zz > sz
      yy = T[i]
      sz = zz
    end
  end
  Y = Array{yy,1}(undef, N+1)
  Z = similar(Y)

  # Get groups
  GRPS = get_unique(S, ["fs", "eltype"], chans)
  i = GRPS[1][1]
  ty = eltype(S.x[i])
  fs = ty(S.fs[i])
  j = 1
  while isapprox(fs, f0)
    j += 1
    j > length(GRPS) && return nothing
    i = GRPS[j][1]
    fs = ty(S.fs[i])
  end
  ff = mkresample(ty(f0/fs))
  rate_old = f0/fs
  si::Int = 0
  ei::Int = 0

  for grp in GRPS
    # get fs, eltype
    c = grp[1]
    ty = eltype(S.x[c])
    fs = ty(S.fs[c])
    fs == f0 && continue
    rate = f0/fs

    # update rate, filter
    if rate != rate_old || ty != eltype(Y)
      ff = mkresample(ty(f0/fs))
      rate_old = rate
    end

    # reinterpret Y if needed
    if ty != eltype(Y)
      Y = reinterpret(ty, isa(Y, Base.ReinterpretArray) ? Y.parent : Y)
      Z = reinterpret(ty, isa(Z, Base.ReinterpretArray) ? Z.parent : Z)
    end

    # Here, things change. We loop in reverse over each segment in each
    # element of grp
    for i in grp
      n_seg = size(S.t[i],1)-1
      gap_inds = Array{Int64,1}(undef, n_seg+1)
      gap_inds[1] = 1
      for k = n_seg:-1:1
        si            = S.t[i][k,1]
        ei            = S.t[i][k+1,1] - (k == n_seg ? 0 : 1)

        # determine boundaries
        nx_in         = ei-si+1
        nx_out        = ceil(Int, nx_in*rate)
        nz_out        = inputlength(ff, nx_out)

        # create padded version of this segment of S.x[i]
        if nz_out > length(Z)
          append!(Z, zeros(ty, nz_out-length(Z)))
        end
        copyto!(Z, 1, S.x[i], si, nx_in)
        if nz_out > nx_in
          Z[nx_in+1:nz_out] = zeros(ty, nz_out-nx_in)
        end
        n_out          = outputlength(ff, nz_out)

        # basically the filt() call in DSP.jl/stream_filt.jl
        ybuf           = view(Y, 1 : n_out)
        ny             = filt!(ybuf, ff, Z[1:nz_out])
        copyto!(S.x[i], si, Y, 1, ny)
        deleteat!(S.x[i], si+ny:ei)
        gap_inds[k+1]  = ceil(Int, ei*rate)
      end
      copyto!(S.t[i], 1, gap_inds, 1, n_seg+1)
      setindex!(S.fs, f0, i)
      note!(S, i, string("resample!, fs=", repr("text/plain", f0, context=:compact=>true)))
    end
  end
  return nothing
end

function resample!(C::GphysChannel, f0::Float64)
  C.fs > 0.0 || error("Can't resample non-timeseries data!")
  (C.fs == f0) && return nothing
  @assert f0 > 0.0

  # This setup is very similar to filtfilt!
  N = nx_max(C)
  ty = eltype(C.x)
  Y = Array{ty,1}(undef, N+1)
  Z = similar(Y)
  rate = f0/C.fs
  ff = mkresample(rate)

  if size(C.t,1) == 2
    nx_in         = length(C.x)
    nx_out        = ceil(Int, nx_in*rate)
    nz_out        = inputlength(ff, nx_out)
    if nz_out > length(Z)
      append!(Z, zeros(ty, nz_out-length(Z)))
    end
    copyto!(Z, 1, C.x, 1, nx_in)
    if nz_out > nx_in
      Z[nx_in+1:nz_out] = zeros(ty, nz_out-nx_in)
    end
    n_out          = outputlength(ff, nz_out)
    ybuf           = view(Y, 1 : n_out)
    ny             = filt!(ybuf, ff, Z[1:nz_out])
    copyto!(C.x, 1, Y, 1, ny)
    deleteat!(C.x, ny+1:nx_in)
  else
    n_seg           = size(C.t, 1)-1
    gap_inds        = ones(Int64, n_seg+1)
    for k = n_seg:-1:1
      si            = C.t[k,1]
      ei            = C.t[k+1,1] - (k == n_seg ? 0 : 1)

      # determine boundaries
      nx_in         = ei-si+1
      nx_out        = ceil(Int, nx_in*rate)
      nz_out        = inputlength(ff, nx_out)

      # create padded version of this segment of C.x
      if nz_out > length(Z)
        append!(Z, zeros(ty, nz_out-length(Z)))
      end
      copyto!(Z, 1, C.x, si, nx_in)
      if nz_out > nx_in
        Z[nx_in+1:nz_out] = zeros(ty, nz_out-nx_in)
      end
      n_out          = outputlength(ff, nz_out)

      # basically the filt() call in DSP.jl/stream_filt.jl
      ybuf           = view(Y, 1 : n_out)
      ny             = filt!(ybuf, ff, Z[1:nz_out])
      copyto!(C.x, si, Y, 1, ny)
      deleteat!(C.x, si+ny:ei)
      gap_inds[k+1]  = ceil(Int, ei*rate)
    end
    copyto!(C.t, 1, gap_inds, 1, n_seg+1)
  end
  setfield!(C, :fs, f0)
  note!(C, string("resample!, fs=", repr("text/plain", f0, context=:compact=>true)))
  return nothing
end

@doc (@doc resample!)
function resample(S::GphysData;
  chans::Union{Integer, UnitRange, Array{Int64,1}}=Int64[],
  fs::Float64=0.0)

  U = deepcopy(S)
  resample!(U, chans=chans, fs=fs)
  return U
end

function resample(C::GphysChannel, f0::Float64)
  U = deepcopy(C)
  resample!(U, f0)
  return U
end