text:
export MonomialVector

# Invariant: Always sorted and no zero vector
struct MonomialVector{C} <: AbstractVector{Monomial{C}}
    vars::Vector{PolyVar{C}}
    Z::Vector{Vector{Int}}

    function MonomialVector{C}(vars::Vector{PolyVar{C}}, Z::Vector{Vector{Int}}) where {C}
        @assert !C || issorted(vars, rev=true)
        @assert all(z -> length(z) == length(vars), Z)
        @assert issorted(Z, rev=true, lt=grlex)
        new(vars, Z)
    end
end
MonomialVector(vars::Vector{PolyVar{C}}, Z::Vector{Vector{Int}}) where {C} = MonomialVector{C}(vars, Z)
MonomialVector{C}() where {C} = MonomialVector{C}(PolyVar{C}[], Vector{Int}[])

# Generate canonical reperesentation by removing variables that are not used
function canonical(m::MonomialVector)
    v = zeros(Bool, nvariables(m))
    for z in m.Z
        v = [v[i] || z[i] > 0 for i in eachindex(v)]
    end
    MonomialVector(_vars(m)[v], Vector{Int}[z[v] for z in m.Z])
end

function Base.hash(m::MonomialVector, u::UInt)
    cm = canonical(m)
    if length(cm.Z) == 0
        hash([], u)
    elseif length(cm.Z) == 1
        hash(Monomial(_vars(cm), cm.Z[1]), u)
    else
        hash(_vars(cm), hash(cm.Z, hash(u)))
    end
end

# vars copied as it may be modifed by `_add_variables!`
# `mutable_copy` recursively copies the vector or vector of integers.
MA.mutable_copy(m::MV) where {MV<:MonomialVector} = MV(copy(m.vars), MA.mutable_copy(m.Z))
Base.copy(m::MonomialVector) = MA.mutable_copy(m)
function Base.getindex(x::MonomialVector, I::AbstractVector{Bool})
    return typeof(x)(x.vars, x.Z[I])
end
function Base.getindex(x::MonomialVector, I::AbstractVector{<:Integer})
    return typeof(x)(x.vars, x.Z[sort(I)])
end
Base.getindex(x::MonomialVector, i::Integer) = Monomial(x.vars, x.Z[i])
Base.getindex(x::MonomialVector, i::CartesianIndex{1}) = x[i.I[1]]

function Base.deleteat!(x::MonomialVector, i)
    deleteat!(x.Z, i)
    return x
end

Base.firstindex(x::MonomialVector) = firstindex(x.Z)
Base.lastindex(x::MonomialVector) = lastindex(x.Z)
Base.size(x::MonomialVector) = (length(x),)
Base.length(x::MonomialVector) = length(x.Z)
Base.isempty(x::MonomialVector) = length(x) == 0
Base.iterate(x::MonomialVector) = isempty(x) ? nothing : (x[1], 1)
function Base.iterate(x::MonomialVector, state::Int)
    state < length(x) ? (x[state+1], state+1) : nothing
end

MultivariatePolynomials.extdegree(x::MonomialVector) = isempty(x) ? (0, 0) : extrema(sum.(x.Z))
MultivariatePolynomials.mindegree(x::MonomialVector) = isempty(x) ? 0 : minimum(sum.(x.Z))
MultivariatePolynomials.maxdegree(x::MonomialVector) = isempty(x) ? 0 : maximum(sum.(x.Z))

_vars(m::Union{Monomial, MonomialVector}) = m.vars

# Recognize arrays of monomials of this module
# [x, y] -> Vector{PolyVar}
# [x, x*y] -> Vector{Monomial}
# [1, x] -> Vector{Term{Int}}
const DMonoVecElemNonConstant{C} = Union{PolyVar{C}, Monomial{C}, Term{C}}
# [1] -> Vector{Int}
const DMonoVecElem{C} = Union{Int, DMonoVecElemNonConstant{C}}
const DMonoVec{C} = AbstractVector{<:DMonoVecElem{C}}

MP.emptymonovec(vars::AbstractVector{PolyVar{C}}) where {C} = MonomialVector{C}(vars, Vector{Int}[])
MP.emptymonovec(t::DMonoVecElemNonConstant) = emptymonovec(_vars(t))
MP.emptymonovec(::Type{<:DMonoVecElemNonConstant{C}}) where {C} = MonomialVector{C}()

function fillZfordeg!(Z, n, deg, ::Type{Val{true}}, filter::Function, ::Int)
    z = zeros(Int, n)
    z[1] = deg
    while true
        if filter(z)
            push!(Z, z)
            z = copy(z)
        end
        if z[end] == deg
            break
        end
        sum = 1
        for j in (n-1):-1:1
            if z[j] != 0
                z[j] -= 1
                z[j+1] += sum
                break
            else
                sum += z[j+1]
                z[j+1] = 0
            end
        end
    end
end
function fillZrec!(Z, z, i, n, deg, filter::Function)
    if deg == 0
        if filter(z)
            push!(Z, copy(z))
        end
    else
        for i in i:i+n-1
            z[i] += 1
            fillZrec!(Z, z, i, n, deg-1, filter)
            z[i] -= 1
        end
    end
end
function fillZfordeg!(Z, n, deg, ::Type{Val{false}}, filter::Function, maxdeg::Int)
    z = zeros(Int, maxdeg * n - maxdeg + 1)
    fillZrec!(Z, z, 1, n, deg, filter)
end
# List exponents in decreasing Graded Lexicographic Order
function getZfordegs(n, degs::AbstractVector{Int}, ::Type{Val{C}}, filter::Function) where C
    Z = Vector{Vector{Int}}()
    # For non-commutative, lower degree need to create a vector of exponent as large as for the highest degree
    maxdeg = isempty(degs) ? 0 : maximum(degs)
    for deg in sort(degs, rev=true)
        fillZfordeg!(Z, n, deg, Val{C}, filter, maxdeg)
    end
    @assert issorted(Z, rev=true, lt=grlex)
    Z
end

function MonomialVector(vars::Vector{PolyVar{true}}, degs::AbstractVector{Int}, filter::Function = x->true)
    MonomialVector{true}(vars, getZfordegs(length(vars), degs, Val{true}, z -> filter(Monomial(vars, z))))
end

function getvarsforlength(vars::Vector{PolyVar{false}}, len::Int)
    n = length(vars)
    map(i -> vars[((i-1) % n) + 1], 1:len)
end
function MonomialVector(vars::Vector{PolyVar{false}}, degs::AbstractVector{Int}, filter::Function = x->true)
    vars = unique!(sort(vars, rev=true))
    Z = getZfordegs(length(vars), degs, Val{false}, z -> filter(Monomial(getvarsforlength(vars, length(z)), z)))
    v = isempty(Z) ? vars : getvarsforlength(vars, length(first(Z)))
    MonomialVector{false}(v, Z)
end
MonomialVector(vars::Vector{<:PolyVar}, degs::Int, filter::Function = x->true) = MonomialVector(vars, [degs], filter)

MP.monomials(vars::AbstractVector{<:PolyVar}, args...) = MonomialVector(vars, args...)
MP.monomials(vars::Tuple{Vararg{PolyVar}}, args...) = monomials([vars...], args...)

#function MP.monomials(vars::TupOrVec{PolyVar{true}}, degs::AbstractVector{Int}, filter::Function = x->true)
#    Z = getZfordegs(length(vars), degs, true, z -> filter(Monomial(vars, z)))
#    [Monomial{true}(vars, z) for z in Z]
#end
#function MP.monomials(vars::TupOrVec{PolyVar{false}}, degs::AbstractVector{Int}, filter::Function = x->true)
#    Z = getZfordegs(length(vars), degs, false, z -> filter(Monomial(vars, z)))
#    v = isempty(Z) ? vars : getvarsforlength(vars, length(first(Z)))
#    [Monomial{false}(v, z) for z in Z]
#end
#MP.monomials(vars::TupOrVec{PV}, degs::Int, filter::Function = x->true) where {PV<:PolyVar} = monomials(vars, [degs], filter)

function buildZvarsvec(::Type{PV}, X::DMonoVec) where {PV<:PolyVar}
    varsvec = Vector{PV}[ (isa(x, DMonoVecElemNonConstant) ? _vars(x) : PolyVar[]) for x in X ]
    allvars, maps = mergevars(varsvec)
    nvars = length(allvars)
    Z = [zeros(Int, nvars) for i in 1:length(X)]
    offset = 0
    for (i, x) in enumerate(X)
        if isa(x, PolyVar)
            @assert length(maps[i]) == 1
            z = [1]
        elseif isa(x, Monomial)
            z = x.z
        elseif isa(x, Term)
            z = x.x.z
        else
            @assert isa(x, Int)
            z = Int[]
        end
        Z[i][maps[i]] = z
    end
    allvars, Z
end

MP.sortmonovec(X::MonomialVector) = (1:length(X), X)
function _sortmonovec(X::DMonoVec{C}) where {C}
    allvars, Z = buildZvarsvec(PolyVar{C}, X)
    σ = sortperm(Z, rev=true, lt=grlex)
    allvars, Z, σ
end
function _removedups!(Z::Vector{Vector{Int}}, σ::Vector{Int})
    dups = findall(i -> Z[σ[i]] == Z[σ[i-1]], 2:length(σ))
    deleteat!(σ, dups)
end
function MP.sortmonovec(X::DMonoVec{C}) where {C}
    if isempty(X)
        Int[], MonomialVector{C}()
    else
        allvars, Z, σ = _sortmonovec(X)
        _removedups!(Z, σ)
        σ, MonomialVector{C}(allvars, Z[σ])
    end
end

function MonomialVector{C}(X::DMonoVec{C}) where C
    allvars, Z = buildZvarsvec(PolyVar{C}, X)
    sort!(Z, rev=true, lt=grlex)
    dups = findall(i -> Z[i] == Z[i-1], 2:length(Z))
    deleteat!(Z, dups)
    MonomialVector{C}(allvars, Z)
end
function MonomialVector(X)
    monovectype(X)(X)
end

MP.monovectype(X::Union{DMonoVecElemNonConstant{C}, Type{<:DMonoVecElemNonConstant{C}}, DMonoVec{C}, Type{<:DMonoVec{C}}}) where {C} = MonomialVector{C}
function MP.monovec(X::DMonoVec)
    MonomialVector(X)
end
MP.monovec(a, mv::MonomialVector) = (a, mv)

function MP.mergemonovec(ms::Vector{MonomialVector{C}}) where {C}
    m = length(ms)
    I = ones(Int, length(ms))
    L = length.(ms)
    X = Vector{Monomial{C}}()
    while any(I .<= L)
        max = nothing
        for i in 1:m
            if I[i] <= L[i]
                x = ms[i][I[i]]
                if max === nothing || max < x
                    max = x
                end
            end
        end
        @assert max !== nothing
        # to ensure that max is no more a union
        max === nothing && return X
        push!(X, max)
        for i in 1:m
            if I[i] <= L[i] && max == ms[i][I[i]]
                I[i] += 1
            end
        end
    end
    # There is no duplicate by construction
    return MonomialVector{C}(buildZvarsvec(PolyVar{C}, X)...)
end

function _add_variables!(monos::MonomialVector{C}, allvars::Vector{PolyVar{C}}, map) where C
    Future.copy!(monos.vars, allvars)
    if !isempty(monos.Z)
        tmp = similar(first(monos.Z))
        for z in monos.Z
            Future.copy!(tmp, z)
            resize!(z, length(allvars))
            fill!(z, 0)
            z[map] = tmp
        end
    end
end

max_stars_repo_name: blegat/DynamicPolynomials

max_issues_repo_name: blegat/DynamicPolynomials

max_forks_repo_name: blegat/DynamicPolynomials


text:
export read_lpg_data

include("parsers.jl")

Base.@kwdef struct Molecule
    coords::Vector{NTuple{3,Float64}} = NTuple{3,Float64}[]
    types::Vector{Int16} = Int16[]
    charges::Vector{Float64} = Float64[]
    masses::Vector{Float64} = Float64[]
    bonds::Vector{Pair{NTuple{2,Int16},Int16}} = Pair{NTuple{2,Int16},Int16}[] #bond is (i, j) => btype
    angles::Vector{Pair{NTuple{3,Int16},Int16}} = Pair{NTuple{3,Int16},Int16}[] #angle is (i, j, k) => atype
    diheds::Vector{Pair{NTuple{4,Int16},Int16}} = Pair{NTuple{4,Int16},Int16}[] #dihedral is (i, j, k, l) => dtype
    improps::Vector{Pair{NTuple{4,Int16},Int16}} = Pair{NTuple{4,Int16},Int16}[] #improper is (i, j, k, l) => itype

    pair_coeffs::Vector{NTuple{2,Float64}} = NTuple{2,Float64}[]

    bond_coeffs::Vector{NTuple{2,Float64}} = NTuple{2,Float64}[]
    bond_map::Vector{Int16} = Int16[]

    angle_coeffs::Vector{NTuple{2,Float64}} = NTuple{2,Float64}[]
    angle_map::Vector{Int16} = Int16[]

    dihed_coeffs::Vector{NTuple{4,Float64}} = NTuple{4,Float64}[]
    dihed_map::Vector{Int16} = Int16[]

    improper_coeffs::Vector{Tuple{Float64,Int8,Int8}} = Tuple{Float64,Int8,Int8}[]
    improper_map::Vector{Int16} = Int16[]

    comment::String
end

"""
    read_lpg_data(f)

Read the datafile produced by LigParGen from `f`. `f` may
be an I/O stream or a file name.
"""
function read_lpg_data(io::IO)
    comment = readline(io)
    molstruct = Molecule(comment = comment[1] == '#' ? comment : "#" * comment)
    natom, nbond, nangle, ndihedral, nimproper,
    ntypes, nbtypes, natypes, ndtypes, nitypes = parse_header(io)

    resize!(molstruct.coords, natom)
    resize!(molstruct.types, natom)
    resize!(molstruct.masses, natom)
    resize!(molstruct.charges, natom)

    resize!(molstruct.bond_map, nbtypes)
    resize!(molstruct.angle_map, natypes)
    resize!(molstruct.dihed_map, ndtypes)
    resize!(molstruct.improper_map, nitypes)

    @assert startswith(readline(io), "Masses")

    read_masses!(molstruct.masses, io, natom)

    readline(io)

    @assert startswith(readline(io), "Pair")

    read_pcoeff!(molstruct.pair_coeffs, molstruct.types, io, ntypes)

    readline(io)

    @assert startswith(readline(io), "Bond")

    read_bcoeff!(molstruct.bond_coeffs, molstruct.bond_map, io, nbtypes)

    readline(io)

    @assert startswith(readline(io), "Angle")

    read_acoeff!(molstruct.angle_coeffs, molstruct.angle_map, io, natypes)

    readline(io)

    @assert startswith(readline(io), "Dihedral")

    read_dcoeff!(molstruct.dihed_coeffs, molstruct.dihed_map, io, ndtypes)

    readline(io)

    @assert startswith(readline(io), "Improper")

    read_icoeff!(molstruct.improper_coeffs, molstruct.improper_map, io, nitypes)

    readline(io)

    @assert startswith(readline(io), "Atoms")

    read_atoms!(molstruct.coords, molstruct.charges, io, natom)

    readline(io)

    @assert startswith(readline(io), "Bonds")

    read_bonds!(molstruct.bonds, io, nbond, molstruct.bond_map)

    readline(io)

    @assert startswith(readline(io), "Angles")

    read_angles!(molstruct.angles, io, nangle, molstruct.angle_map)

    readline(io)

    @assert startswith(readline(io), "Dihedrals")

    read_dihed!(molstruct.diheds, io, ndihedral, molstruct.dihed_map)

    readline(io)

    @assert startswith(readline(io), "Impropers")

    read_impropers!(molstruct.improps, io, nimproper, molstruct.improper_map)

    return molstruct
end

read_lpg_data(fname::AbstractString) = open(fname) do io read_lpg_data(io) end

ext: jl


text:
@inline function eatwhitespaces(str::Union{VectorBackedUTF8String, String}, i=1, len=lastindex(str))
    while i<=len
        @inbounds b = codeunit(str, i)

        if b==0x20 # This is ' '
            i += 1
        else
            break
        end
    end
    return i
end

@inline function eatnewlines(str::Union{VectorBackedUTF8String, String}, i=1, len=lastindex(str))
    count = 0
    while i<=len
        @inbounds b = codeunit(str, i)
        if b == 0xd # '\r'
            i += 1
            if i<=len
                @inbounds b = codeunit(str, i)
                if b == 0xa # '\n'
                    i += 1
                end
            end
            count += 1
        elseif b == 0xa
            i += 1
            if i<=len
                @inbounds b = codeunit(str, i)
                if b == 0xd
                    i += 1
                end
            end
            count += 1
        else
            break
        end
    end

    return i, count
end

@inline function tryparsenext_base10_digit(T,str::Union{VectorBackedUTF8String, String},i, len)
    i > len && @goto error
    @inbounds b = codeunit(str,i)
    diff = b-0x30
    diff >= UInt8(10) && @goto error
    return convert(T, diff), i+1

    @label error
    return nothing
end

@inline _isdigit(b::UInt8) = ( (0x30 ≤ b) & (b ≤ 0x39) )

@inline function parse_uint_and_stop(str::Union{VectorBackedUTF8String, String}, i, len, n::T) where {T <: Integer}
    ten = T(10)
    # specialize handling of the first digit so we can return an error
    max_without_overflow = div(typemax(T)-9,10) # the larg
    i <= len || return n, false, i
    @inbounds b = codeunit(str, i)
    diff = b-0x30
    if diff < UInt8(10) && n <= max_without_overflow
        n *= ten
        n += T(diff)
    else
        return n, false, i
    end
    i += 1

    while i <= len && n <= max_without_overflow
        @inbounds b = codeunit(str, i)
        diff = b-0x30
        if diff < UInt8(10)
            n *= ten
            n += T(diff)
        else
            return n, true, i
        end
        i += 1
    end
    return n, true, i
end

@inline function read_digits(str::Union{VectorBackedUTF8String, String}, i, len)
    # slurp up extra digits
    while i <= len
        @inbounds b = codeunit(str, i)
        if !_isdigit(b) # do nothing
            return i
        end
        i += 1
    end
    return i
end

@inline function _is_e(str::Union{VectorBackedUTF8String, String}, i)
    @inbounds b = codeunit(str, i)
    return  (b==0x65) | (b==0x45)
end

@inline function _is_negative(str::Union{VectorBackedUTF8String, String}, i)
    @inbounds b = codeunit(str, i)
    return b==0x2d
end

@inline function _is_positive(str::Union{VectorBackedUTF8String, String}, i)
    @inbounds b = codeunit(str, i)
    return b==0x2b
end

const pre_comp_exp = Float64[10.0^i for i=0:22]

@inline function tryparsenext(::Numeric{F}, str::Union{VectorBackedUTF8String, String}, i, len) where {F<:AbstractFloat}
    R = Nullable{F}

    i>len && @goto error

    negate = false
    @inbounds b = codeunit(str, i)
    if b==0x2d # '-'
        negate = true
        i += 1
    elseif b==0x2b # '+'
        i +=1
    end

    f1::Int64 = 0

    # read an integer up to the decimal point
    f1, rval1, idecpt = parse_uint_and_stop(str, i, len, f1)
    idecpt = read_digits(str, idecpt, len) # get any trailing digits
    i = idecpt

    ie = i
    frac_digits = 0

    # next thing must be dec pt.
    if i <= len && @inbounds(codeunit(str, i)) == 0x2e # Check for '.'
        i += 1
        f1, rval2, ie = parse_uint_and_stop(str, i, len, f1)
        frac_digits = ie - i

        ie = read_digits(str, ie, len) # get any trailing digits
    elseif !rval1 # no first number, and now no deciaml point => invalid
        @goto error
    end

    # Next thing must be exponent
    i = ie
    eval::Int32 = 0

    if i <= len && _is_e(str, i)
        i += 1

        enegate = false
        if i<=len
            if _is_negative(str, i)
                enegate = true
                i += 1
            elseif _is_positive(str, i)
                i += 1
            end
        end
        eval, rval3, i = parse_uint_and_stop(str, i, len, eval)
        if enegate
            eval *= Int32(-1)
        end
    end

    exp = eval - frac_digits

    maxexp = 308
    minexp = -307

    if frac_digits <= 15 && -22 <= exp <= 22
        if exp >= 0
            f = F(f1)*pre_comp_exp[exp+1]
        else
            f = F(f1)/pre_comp_exp[-exp+1]
        end
    else
          f = convert_to_double(f1, exp)
    end

    if negate
        f = -f
    end

    @label done
    return R(convert(F, f)), i

    @label error
    return R(), i
end

function tryparsenext(f::Field{T}, str::Union{VectorBackedUTF8String, String}, i, len, opts::LocalOpts{T_ENDCHAR}) where {T, T_ENDCHAR<:UInt8}
    R = Nullable{T}
    i > len && @goto error
    if f.ignore_init_whitespace
        i = eatwhitespaces(str, i, len)
    end
    @chk2 res, i = tryparsenext(f.inner, str, i, len, opts)

    if f.ignore_end_whitespace
        i0 = i

        while i<=len
            @inbounds b = codeunit(str, i)

            !opts.spacedelim && opts.endchar == 0x09 && b == 0x09 && (i = i+1; @goto done) # 0x09 is \t

            b!=0x20 && b!=0x09 && break
            i=i+1
        end

        opts.spacedelim && i > i0 && @goto done
    end
    # todo don't ignore whitespace AND spacedelim

    if i > len
        if f.eoldelim
            @goto done
        else
            @goto error
        end
    end

    i>len && error("Internal error.")
    @inbounds b = codeunit(str, i)
    opts.spacedelim && (b!=0x20 || b!=0x09) && (i+=1; @goto done)
    !opts.spacedelim && opts.endchar == b && (i+=1; @goto done)

    if f.eoldelim
        if b == 0x0d # '\r'
            i+=1
            if i<=len
                @inbounds b = codeunit(str, i)
                if b == 0x0a # '\n'
                    i+=1
                end
            end
            @goto done
        elseif b == 0x0a # '\n'
            i+=1
            if i<=len
                @inbounds b = codeunit(str, i)
                if b == 0x0d # '\r'
                    i+=1
                end
            end
            @goto done
        end
    end

    @label error
    return R(), i

    @label done
    return R(convert(T, res)), i
end

function tryparsenext(q::Quoted{T,S,<:UInt8,<:UInt8}, str::Union{VectorBackedUTF8String, String}, i, len, opts::LocalOpts{<:UInt8,<:UInt8,<:UInt8}) where {T,S}
    if i>len
        q.required && @goto error
        # check to see if inner thing is ok with an empty field
        @chk2 x, i = tryparsenext(q.inner, str, i, len, opts) error
        @goto done
    end
    @inbounds b = codeunit(str, i)
    ii = i+1
    quotestarted = false
    if q.quotechar == b
        quotestarted = true
        if !q.includequotes
            i = ii
        end

        if q.stripwhitespaces
            i = eatwhitespaces(str, i, len)
        end
    else
        q.required && @goto error
    end

    if quotestarted
        qopts = LocalOpts(q.quotechar, false, q.quotechar, q.escapechar,
                         q.includequotes, q.includenewlines)
        @chk2 x, i = tryparsenext(q.inner, str, i, len, qopts)
    else
        @chk2 x, i = tryparsenext(q.inner, str, i, len, opts)
    end

    if i > len
        if quotestarted && !q.includequotes
            @goto error
        end
        @goto done
    end

    if q.stripwhitespaces
        i = eatwhitespaces(str, i, len)
    end
    i>len && error("Internal error.")
    @inbounds b = codeunit(str, i)
    ii = i + 1

    if quotestarted && !q.includequotes
        b != q.quotechar && @goto error
        i = ii
    end

    @label done
    return Nullable{T}(x), i

    @label error
    return Nullable{T}(), i
end


text:
# ---
# title: 818. Race Car
# id: problem818
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Dynamic Programming, Heap
# link: <https://leetcode.com/problems/race-car/description/>
# hidden: true
# ---
#
# Your car starts at position 0 and speed +1 on an infinite number line.  (Your
# car can go into negative positions.)
#
# Your car drives automatically according to a sequence of instructions A
# (accelerate) and R (reverse).
#
# When you get an instruction "A", your car does the following: `position +=
# speed, speed *= 2`.
#
# When you get an instruction "R", your car does the following: if your speed is
# positive then `speed = -1` , otherwise `speed = 1`.  (Your position stays the
# same.)
#
# For example, after commands "AAR", your car goes to positions 0->1->3->3, and
# your speed goes to 1->2->4->-1.
#
# Now for some target position, say the **length** of the shortest sequence of
# instructions to get there.
#
#
#
#     **Example 1:**
#     Input:
#     target = 3
#     Output: 2
#     Explanation:
#     The shortest instruction sequence is "AA".
#     Your position goes from 0->1->3.
#
#
#
#     **Example 2:**
#     Input:
#     target = 6
#     Output: 5
#     Explanation:
#     The shortest instruction sequence is "AAARA".
#     Your position goes from 0->1->3->7->7->6.
#
#
#
#
# **Note:**
#
#   * `1 <= target <= 10000`.
#
#
## @lc code=start
using LeetCode

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


text:
function test()
    a = 0
    for i in 1:10000
        a = a + 1
    end
    println(a)
end

@time test()


text:
using StatsBase: mean
using NQCDynamics: TimeCorrelationFunctions, DynamicsOutputs
using .TimeCorrelationFunctions: TimeCorrelationFunction

"""
Sum the outputs from each trajectory.
"""
struct SumReduction end
function (reduction::SumReduction)(u, batch, I)
    sum_outputs!(u, batch)
    return (u, false)
end

"""
Average the outputs over all trajectories.
"""
mutable struct MeanReduction
    counter::Int
end

function (reduction::MeanReduction)(u, batch, I)
    u .*= reduction.counter
    reduction.counter += length(batch)
    sum_outputs!(u, batch)
    u ./= reduction.counter
    return (u, false)
end

function sum_outputs!(u, batch)
    for i=1:length(batch)
        if size(batch[i]) != size(u)
            @warn "DimensionMismatch encountered when reducing outputs. \
                Discarding trajectory. Check `u_init` matches the output size."
        else
            u .+= batch[i]
        end
    end
end

struct AppendReduction end
(::AppendReduction)(u,data,I) = (append!(u,data), false)

"""
    Reduction(reduction::Symbol)

Converts reduction keyword into function that performs reduction.

* `:discard` is used internally when using callbacks to save outputs instead.
"""
function Reduction(reduction::Symbol)
    if reduction === :mean
        return MeanReduction(0)
    elseif reduction === :sum
        return SumReduction()
    elseif reduction === :append
        return AppendReduction()
    else
        throw(ArgumentError("`reduction` $reduction not recognised."))
    end
end

function get_u_init(::Union{MeanReduction, SumReduction}, saveat, stripped_kwargs, tspan, u0, output::TimeCorrelationFunction)
    savepoints = get_savepoints(saveat, stripped_kwargs, tspan)
    u_init = [TimeCorrelationFunctions.correlation_template(output) for _ ∈ savepoints]
    return u_init
end

function get_u_init(::Union{MeanReduction, SumReduction}, saveat, stripped_kwargs, tspan, u0, output::DynamicsOutputs.EnsembleSaver)
    savepoints = get_savepoints(saveat, stripped_kwargs, tspan)
    u_init = DynamicsOutputs.output_template(output, savepoints, u0)
    return u_init
end

function get_u_init(::Union{MeanReduction, SumReduction}, saveat, stripped_kwargs, tspan, u0, output)
    return output_template(output, u0)
end

function get_u_init(::AppendReduction, saveat, stripped_kwargs, tspan, u0, output)
    return nothing
end

function get_savepoints(saveat, stripped_kwargs, tspan)
    if saveat != []
        if saveat isa Number
            savepoints = tspan[1]:saveat:tspan[2]
        else
            savepoints = saveat
        end
    elseif haskey(stripped_kwargs, :dt)
        dt = stripped_kwargs[:dt]
        savepoints = tspan[1]:dt:tspan[2]
    else
        throw(ArgumentError("Make sure to pass either `saveat` or `dt` as keyword arguments
            to compute average/sum over many trajectories."))
    end
end


text:
@testset "Triangle" begin

p1 = Point3(0)
p2 = Point3(1, 0, 0)
p3 = Point3(0, 1, 0)

Tri = GRAYTR.Triangle(p1, p2, p3)

@test GRAYTR.can_intersect(Tri)
@test GRAYTR.area(Tri) ≈ 0.5

@test GRAYTR.obj_bounds(Tri) == GRAYTR.BoundingBox(Point3(0), Point3(1, 1, 0))
@test GRAYTR.world_bounds(Tri) == GRAYTR.BoundingBox(Point3(0), Point3(1, 1, 0))

T = GRAYTR.translation(-0.5, -0.5, -1)
@test GRAYTR.world_bounds(T(Tri)) == GRAYTR.BoundingBox(Point3(-0.5, -0.5, -1), Point3(0.5, 0.5, -1))

Tri2 = GRAYTR.Triangle(T)

R = GRAYTR.Ray(Point3(0.2, 0.2, 3), Vector3(0, 0, -1))
@test GRAYTR.intersectP(R, Tri)

R = GRAYTR.Ray(Point3(0.6, 0.6, 3), Vector3(0, 0, -1))
@test !GRAYTR.intersectP(R, Tri)

R = GRAYTR.Ray(Point3(2, 1, 3), Vector3(0, 0, -1))
@test !GRAYTR.intersectP(R, Tri)

R = GRAYTR.Ray(Point3(2, 1, 3), Vector3(0, 0, 1))
@test !GRAYTR.intersectP(R, Tri)

R = GRAYTR.Ray(Point3(0.2, 0.2, 0), Vector3(1, 0, 0))
@test !GRAYTR.intersectP(R, Tri)

R = GRAYTR.Ray(Point3(0.2, 0.2, 2), Vector3(0, 0, -1), 5.0, Inf, 1)
@test !GRAYTR.intersectP(R, Tri)

R = GRAYTR.Ray(Point3(0.2, 0.2, 3), Vector3(0, 0, -1))
dg, t, e = GRAYTR.shape_intersect(R, Tri)
@test dg.n == Normal3(0, 0, 1)
@test t ≈ 3.0

R = GRAYTR.Ray(Point3(0.6, 0.6, 3), Vector3(0, 0, -1))
dg, t, e = GRAYTR.shape_intersect(R, Tri)
@test dg == nothing

R = GRAYTR.Ray(Point3(2, 1, 3), Vector3(0, 0, -1))
dg, t, e = GRAYTR.shape_intersect(R, Tri)
@test dg == nothing

R = GRAYTR.Ray(Point3(2, 1, 3), Vector3(0, 0, 1))
dg, t, e = GRAYTR.shape_intersect(R, Tri)
@test dg == nothing

R = GRAYTR.Ray(Point3(0.2, 0.2, 0), Vector3(1, 0, 0))
dg, t, e = GRAYTR.shape_intersect(R, Tri)
@test dg == nothing

R = GRAYTR.Ray(Point3(0.2, 0.2, 2), Vector3(0, 0, -1), 5.0, Inf, 1)
dg, t, e = GRAYTR.shape_intersect(R, Tri)
@test dg == nothing

end

max_stars_repo_name: dronir/GRAYTR

max_issues_repo_name: dronir/GRAYTR

max_forks_repo_name: dronir/GRAYTR


text:
module Day14
export day14

function parseMask(mask)
    fixed = 0
    variable = 0
    for (power, index) in enumerate(reverse(1:length(mask)))
        char = mask[index]
        if char == 'X'
            variable += 1 << (power-1)
        elseif char == '1'
            fixed += 1 << (power-1)
        end
    end
    :mask, fixed, variable
end

function parseLine(line)
    mask = match(r"^mask = ([X01]{36})$", line)
    if mask !== nothing
        return parseMask(mask[1])
    end

    mem = match(r"^mem\[(\d*)\] = (\d*)$", line)
    if mem !== nothing
        return (:mem, parse(Int, mem[1]), parse(Int, mem[2]))
    end
    error("Strange line: ", line)
end

"Parse input to a list of lines.
Possible types:
 - :mem, adress, value
 - :mask, fixedbits, variablebits"
function readInput(inputfile="inputs/input14.txt")
    [parseLine(line) for line in eachline(inputfile)]
end

function day14(inputfile="inputs/input14.txt")
    values = readInput(inputfile)
    (day14Part1(values), day14Part2(values))
end

function day14Part1(input)
    currentFixed = 0
    currentVariable = 0
    memory = Dict{Int, Int}()
    for (type, p1, p2) in input
        if type === :mask
            currentFixed = p1
            currentVariable = p2
        elseif type === :mem
            memory[p1] = (currentFixed + (currentVariable & p2))
        end
    end
    sum(values(memory))
end

function computeFluctuations(value)
    power = 0
    non_zeros = []
    while value > 0
        value, r = divrem(value, 2)
        if r == 1
            push!(non_zeros, 0:2^power:2^power)
        end
        power += 1
    end
    vec(map(sum, Iterators.product(non_zeros...)))
end

function day14Part2(input)
    # use fixed Mask with OR (|) now -> 0 no change -> 1 overwritten
    currentFixed = 0
    currentVariable = [0]
    memory = Dict{Int, Int}()
    for (type, p1, p2) in input
        if type === :mask
            currentFixed = p1
            currentVariable = computeFluctuations(p2)
        elseif type === :mem
            for addr in currentVariable .+ (p1 | currentFixed)
                memory[addr] = p2
            end
        end
    end
    sum(values(memory))
end

end #module


text:
# Test suite for DualClockFifo
@testset "Testing Dual Clock Fifo" begin
    # Instantiate a small fifo for easier testing. Everything should be
    # parameterized to scale to larger fifos and fifos with different elements.
    buffersize = 3
    fifo = AsapSim.DualClockFifo(Int64, buffersize)

    @test AsapSim.buffersize(fifo) == buffersize
    @test AsapSim.isreadempty(fifo) == true

    # Do some testing of the "isfull" function, which is independent of the
    # actual fifo implementation.
    @test AsapSim.isfull(1, 1, 3) == false
    @test AsapSim.isfull(2, 1, 3) == false
    @test AsapSim.isfull(3, 1, 3) == true
    @test AsapSim.isfull(3, 2, 3) == false
    @test AsapSim.isfull(1, 2, 3) == true
    @test AsapSim.isfull(2, 3, 3) == true

    # Verify that we can write to the fifo. Should be able to write two times
    # before receiving an invalid write signal.
    @test AsapSim.iswriteready(fifo) == true
    AsapSim.write(fifo, 10)
    AsapSim.writeupdate!(fifo)

    @test AsapSim.iswriteready(fifo) == true
    AsapSim.write(fifo, 20)
    AsapSim.writeupdate!(fifo)

    @test AsapSim.iswriteready(fifo) == false

    # ----------------- #
    # Test some reading #
    # ----------------- #
    # No read should be ready because we haven't requested a read yet.
    @test AsapSim.isreadready(fifo) == false
    AsapSim.read(fifo)
    AsapSim.readupdate!(fifo)

    @test AsapSim.isreadready(fifo) == true
    @test AsapSim.read(fifo) == 10

    # Set another read request before the next update.
    AsapSim.x!(fifo)
    AsapSim.readupdate!(fifo)
    @test AsapSim.isreadready(fifo) == true
    @test AsapSim.read(fifo) == 20
    @test AsapSim.isreadempty(fifo) == true

    # Ensure we can keep reading from the fifo until we ask it to increment again.
    AsapSim.readupdate!(fifo)
    @test AsapSim.isreadready(fifo) == true
    AsapSim.increment!(fifo)
    AsapSim.readupdate!(fifo)
    @test AsapSim.isreadready(fifo) == false

    # Do another write
    AsapSim.write(fifo, 30)
    AsapSim.writeupdate!(fifo)
    AsapSim.readupdate!(fifo)
    @test AsapSim.isreadready(fifo) == true
    @test AsapSim.read(fifo) == 30
end

@testset "Running Long Fifo Test" begin

    # The goal here is to simulate two different clocks writing to the fifo over
    # a long period of time.
    #
    # We will sweep through different "frequencies" for the read and write side,
    # writing random numbers and reading. A circular dequeue from DataStrucutres
    # will be used to ensure that data is being maintained correctly.
    #
    # 5 different cases will be tried:
    #
    # - write side much faster than read side
    # - write side slightly faster than read side
    # - both sides equal
    # - write side slightly slower than read side
    # - write side much slower than read side.

    write_periods = (1, 9, 10, 11, 20)
    read_periods  = (10, 10, 10, 10, 10)

    # Set up some parameters to allow variable testing if desired.
    buffersize = 32
    datatype = Int

    # Use the CircularDeque as a golden reference.
    que = DataStructures.CircularDeque{datatype}(buffersize)
    fifo = AsapSim.DualClockFifo(datatype, buffersize)

    # Number of "clock cycles" to run each test for.
    test_cycles = 10_000

    # Initialize random number generator seed for reproducibility
    srand(0)
    for (wp, rp) in zip(write_periods, read_periods)
        for i in 1:test_cycles
            # Do a read if possible and if we're on the correct cycle.
            if mod(i, rp) == 0
                if AsapSim.isreadready(fifo)
                    AsapSim.increment!(fifo)
                    readdata = AsapSim.read(fifo)
                    # Check that it matches the data from the que
                    @test readdata == shift!(que)
                end
                # clock the fifo
                AsapSim.readupdate!(fifo)
            end

            # Write if possible and on the correct cycle
            if mod(i, wp) == 0
                if AsapSim.iswriteready(fifo)
                    writedata = rand(datatype)
                    push!(que, writedata)
                    AsapSim.write(fifo, writedata)
                end
                AsapSim.writeupdate!(fifo)
            end
        end
    end
end


text:
# load a BioNetGen defined BCR network, solve the ODEs
# plot a specific observable and compare to the BioNetGen solution

using DiffEqBase, Catalyst, Plots, OrdinaryDiffEq, Sundials, DataFrames, CSVFiles, LinearAlgebra, SparseArrays
using ReactionNetworkImporters
using TimerOutputs

# parameters
doplot = false
tf = 10000.
build_jac = false
sparse_jac = false
zeroout_jac = false

# BNG simulation data
datadir  = joinpath(@__DIR__,"../data/bcr")
fname    = joinpath(datadir, "bcr.net")
gdatfile = joinpath(datadir, "bcr.gdat")
print("getting gdat file...")
gdatdf = DataFrame(load(File(format"CSV", gdatfile), header_exists=true, spacedelim=true) )
println("done")

# we'll time the DiffEq solvers
const to = TimerOutput()
reset_timer!(to)

# BioNetGen network
@timeit to "bionetgen" prnbng = loadrxnetwork(BNGNetwork(),fname);
show(to)
rnbng = prnbng.rn; u0 = prnbng.u₀; p = prnbng.p;
@timeit to "bODESystem" bosys = convert(ODESystem, rnbng)
show(to)
@timeit to "bODEProb" boprob = ODEProblem(bosys, Pair.(species(rnbng),u0), (0.,tf), Pair.(parameters(rnbng),p))
show(to)
u = copy(u0);
du = similar(u);
@timeit to "f1" boprob.f(du,u,p,0.)
@timeit to "f2" boprob.f(du,u,p,0.)
if build_jac
    J = zeros(length(u),length(u))
    #J = similar(rnbng.odefun.jac_prototype)
    @timeit to "J1" rnbng.jac(J,u,p,0.)
    @timeit to "J2" rnbng.jac(J,u,p,0.)
end
show(to)
println()

# BNG simulation results for Activated Syk
asykgroups = prnbng.groupstoids[:Activated_Syk]
#asyksyms = findall(x -> x ∈ asykgroups, rnbng.syms_to_ints)
# asynbng = zeros(length(gdatdf[:time]))
# for sym in asyksyms
#     global asynbng
#     asynbng += gdatdf[sym]
# end

# DiffEq solver
#reset_timer!(to);
#@timeit to "BNG_CVODE_BDF-LU-1" begin bsol = solve(boprob, CVODE_BDF(), saveat=1., abstol=1e-8, reltol=1e-8); end; show(to)
#@timeit to "BNG_CVODE_BDF-LU-2" begin bsol = solve(boprob, CVODE_BDF(), saveat=1., abstol=1e-8, reltol=1e-8); end; show(to)
#@timeit to "BNG_CVODE_BDF-GMRES-1" begin bsol = solve(boprob, CVODE_BDF(linear_solver=:GMRES), saveat=1., abstol=1e-8, reltol=1e-8); end; show(to)
#@timeit to "BNG_CVODE_BDF-GMRES-2" begin bsol = solve(boprob, CVODE_BDF(linear_solver=:GMRES), saveat=1., abstol=1e-8, reltol=1e-8); end; show(to)
# #reset_timer!(to); @timeit to "BNG_RODAS5_BDF" begin bsol2 = solve(boprob, rodas5(autodiff=false), saveat=1., abstol=1e-8, reltol=1e-8); end; show(to)
#@timeit to "KenCarp4-1" begin sol = solve(boprob,KenCarp4(autodiff=false,linsolve=LinSolveFactorize(lu)), saveat=1., abstol=1e-8, reltol=1e-8); end; show(to)
#@timeit to "KenCarp4-2" begin bsol = solve(boprob,KenCarp4(autodiff=false,linsolve=LinSolveFactorize(lu)), saveat=1., abstol=1e-8, reltol=1e-8); end; show(to)
@timeit to "KenCarp4-1" begin sol = solve(boprob,KenCarp4(autodiff=false), abstol=1e-8, reltol=1e-8, saveat=1.); end; show(to)
@timeit to "KenCarp4-2" begin bsol = solve(boprob,KenCarp4(autodiff=false), abstol=1e-8, reltol=1e-8, saveat=1.); end; show(to)

# Activated Syk from DiffEq
basyk = sum(bsol[asykgroups,:], dims=1)

if doplot
    plotlyjs()
    plot(gdatdf[!,:time][2:end], gdatdf[!,:Activated_Syk][2:end], xscale=:log10, label="AsykGroup", linestyle=:dot)
#     # plot!(cdatdf[:time][2:end], asynbng[2:end], xscale=:log10, label=:AsykSum)
    plot!(bsol.t[2:end], basyk'[2:end], label="AsykDEBio", xscale=:log10)
end

# test the error, note may be large in abs value though small relatively
# #norm(gdatdf[:Activated_Syk] - asynbng, Inf)
# #norm(asynbng - basyk', Inf)
norm(gdatdf[!,:Activated_Syk] - basyk', Inf)

#@assert all(abs.(gdatdf[!,:Activated_Syk] - basyk') .< 1e-6 * abs.(gdatdf[:Activated_Syk]))


text: include("ObservationsTests.jl")


text:
# Internal routine: integrate f over the union of the open intervals
# (s[1],s[2]), (s[2],s[3]), ..., (s[end-1],s[end]), using h-adaptive
# integration with the order-n Kronrod rule and weights of type Tw,
# with absolute tolerance atol and relative tolerance rtol,
# with maxevals an approximate maximum number of f evaluations.
function do_quadgk(f::F, s::NTuple{N,T}, n, atol, rtol, maxevals, nrm) where {T,N,F}
    x,w,gw = cachedrule(eltype(s),n)

    @assert N ≥ 2
    segs = ntuple(i -> evalrule(f, s[i],s[i+1], x,w,gw, nrm), Val{N-1}())
    I = sum(s -> s.I, segs)
    E = sum(s -> s.E, segs)
    numevals = (2n+1) * (N-1)

    # logic here is mainly to handle dimensionful quantities: we
    # don't know the correct type of atol, in particular, until
    # this point where we have the type of E from f.  Also, follow
    # Base.isapprox in that if atol≠0 is supplied by the user, rtol
    # defaults to zero.
    atol_ = something(atol, zero(E))
    rtol_ = something(rtol, iszero(atol_) ? sqrt(eps(one(eltype(x)))) : zero(eltype(x)))

    # optimize common case of no subdivision
    if E ≤ atol_ || E ≤ rtol_ * nrm(I) || numevals ≥ maxevals
        return (I, E) # fast return when no subdivisions required
    end

    return adapt(f, heapify!(collect(segs), Reverse), I, E, numevals, x,w,gw,n, atol_, rtol_, maxevals, nrm)
end

# internal routine to perform the h-adaptive refinement of the integration segments (segs)
function adapt(f, segs::Vector{T}, I, E, numevals, x,w,gw,n, atol, rtol, maxevals, nrm) where {T}
    # Pop the biggest-error segment and subdivide (h-adaptation)
    # until convergence is achieved or maxevals is exceeded.
    while E > atol && E > rtol * nrm(I) && numevals < maxevals
        s = heappop!(segs, Reverse)
        mid = (s.a + s.b) / 2
        s1 = evalrule(f, s.a, mid, x,w,gw, nrm)
        s2 = evalrule(f, mid, s.b, x,w,gw, nrm)
        I = (I - s.I) + s1.I + s2.I
        E = (E - s.E) + s1.E + s2.E
        numevals += 4n+2

        # handle type-unstable functions by converting to a wider type if needed
        Tj = promote_type(typeof(s1), promote_type(typeof(s2), T))
        if Tj !== T
            return adapt(f, heappush!(heappush!(Vector{Tj}(segs), s1, Reverse), s2, Reverse),
                         I, E, numevals, x,w,gw,n, atol, rtol, maxevals, nrm)
        end

        heappush!(segs, s1, Reverse)
        heappush!(segs, s2, Reverse)
    end

    # re-sum (paranoia about accumulated roundoff)
    I = segs[1].I
    E = segs[1].E
    for i in 2:length(segs)
        I += segs[i].I
        E += segs[i].E
    end
    return (I, E)
end

# transform f and the endpoints s to handle infinite intervals, if any,
# and pass transformed data to workfunc(f, s, tfunc)
function handle_infinities(workfunc, f, s)
    s1, s2 = s[1], s[end]
    if s1 isa Real && s2 isa Real # check for infinite or semi-infinite intervals
        inf1, inf2 = isinf(s1), isinf(s2)
        if inf1 || inf2
            if inf1 && inf2 # x = t/(1-t^2) coordinate transformation
                return workfunc(t -> begin t2 = t*t; den = 1 / (1 - t2);
                                            f(t*den) * (1+t2)*den*den; end,
                                map(x -> isinf(x) ? copysign(one(x), x) : 2x / (1+hypot(1,2x)), s),
                                t -> t / (1 - t^2))
            end
            let (s0,si) = inf1 ? (s2,s1) : (s1,s2) # let is needed for JuliaLang/julia#15276
                if si < 0 # x = s0 - t/(1-t)
                    return workfunc(t -> begin den = 1 / (1 - t);
                                                f(s0 - t*den) * den*den; end,
                                    reverse(map(x -> 1 / (1 + 1 / (s0 - x)), s)),
                                    t -> s0 - t/(1-t))
                else # x = s0 + t/(1-t)
                    return workfunc(t -> begin den = 1 / (1 - t);
                                                f(s0 + t*den) * den*den; end,
                                    map(x -> 1 / (1 + 1 / (x - s0)), s),
                                    t -> s0 + t/(1-t))
                end
            end
        end
    end
    return workfunc(f, s, identity)
end

# Gauss-Kronrod quadrature of f from a to b to c...

"""
    quadgk(f, a,b,c...; rtol=sqrt(eps), atol=0, maxevals=10^7, order=7, norm=norm)

Numerically integrate the function `f(x)` from `a` to `b`, and optionally over additional
intervals `b` to `c` and so on. Keyword options include a relative error tolerance `rtol`
(if `atol==0`, defaults to `sqrt(eps)` in the precision of the endpoints), an absolute error tolerance
`atol` (defaults to 0), a maximum number of function evaluations `maxevals` (defaults to
`10^7`), and the `order` of the integration rule (defaults to 7).

Returns a pair `(I,E)` of the estimated integral `I` and an estimated upper bound on the
absolute error `E`. If `maxevals` is not exceeded then `E <= max(atol, rtol*norm(I))`
will hold. (Note that it is useful to specify a positive `atol` in cases where `norm(I)`
may be zero.)

The endpoints `a` et cetera can also be complex (in which case the integral is performed over
straight-line segments in the complex plane). If the endpoints are `BigFloat`, then the
integration will be performed in `BigFloat` precision as well.

!!! note
    It is advisable to increase the integration `order` in rough proportion to the
    precision, for smooth integrands.

More generally, the precision is set by the precision of the integration
endpoints (promoted to floating-point types).

The integrand `f(x)` can return any numeric scalar, vector, or matrix type, or in fact any
type supporting `+`, `-`, multiplication by real values, and a `norm` (i.e., any normed
vector space). Alternatively, a different norm can be specified by passing a `norm`-like
function as the `norm` keyword argument (which defaults to `norm`).

!!! note
    Only one-dimensional integrals are provided by this function. For multi-dimensional
    integration (cubature), there are many different algorithms (often much better than simple
    nested 1d integrals) and the optimal choice tends to be very problem-dependent. See the
    Julia external-package listing for available algorithms for multidimensional integration or
    other specialized tasks (such as integrals of highly oscillatory or singular functions).

The algorithm is an adaptive Gauss-Kronrod integration technique: the integral in each
interval is estimated using a Kronrod rule (`2*order+1` points) and the error is estimated
using an embedded Gauss rule (`order` points). The interval with the largest error is then
subdivided into two intervals and the process is repeated until the desired error tolerance
is achieved.

These quadrature rules work best for smooth functions within each interval, so if your
function has a known discontinuity or other singularity, it is best to subdivide your
interval to put the singularity at an endpoint. For example, if `f` has a discontinuity at
`x=0.7` and you want to integrate from 0 to 1, you should use `quadgk(f, 0,0.7,1)` to
subdivide the interval at the point of discontinuity. The integrand is never evaluated
exactly at the endpoints of the intervals, so it is possible to integrate functions that
diverge at the endpoints as long as the singularity is integrable (for example, a `log(x)`
or `1/sqrt(x)` singularity).

For real-valued endpoints, the starting and/or ending points may be infinite. (A coordinate
transformation is performed internally to map the infinite interval to a finite one.)
"""
quadgk(f, a, b, c...; kws...) =
    quadgk(f, promote(a, b, c...)...; kws...)

quadgk(f, a::T,b::T,c::T...;
       atol=nothing, rtol=nothing, maxevals=10^7, order=7, norm=norm) where {T} =
    handle_infinities(f, (a, b, c...)) do f, s, _
        do_quadgk(f, s, order, atol, rtol, maxevals, norm)
    end


text:
function copper_plate(canonical::Canonical, expression::Symbol, bus_count::Int64)

    time_steps = model_time_steps(canonical)
    devices_netinjection = _remove_undef!(canonical.expressions[expression])

    canonical.constraints[:CopperPlateBalance] = JuMPConstraintArray(undef, time_steps)

    for t in time_steps
        canonical.constraints[:CopperPlateBalance][t] = JuMP.@constraint(canonical.JuMPmodel, sum(canonical.expressions[expression].data[1:bus_count, t]) == 0)
    end

    return

end


text:
using Oceananigans
using GLMakie
using JLD2

filename = "eddying_channel_catke"

zonal_file = jldopen(filename * "_zonal_average.jld2")

B = FieldTimeSeries(filename * "_zonal_average.jld2", "b")
U = FieldTimeSeries(filename * "_zonal_average.jld2", "u")
V = FieldTimeSeries(filename * "_zonal_average.jld2", "v")
bb = FieldTimeSeries(filename * "_zonal_average.jld2", "bb")
wb = FieldTimeSeries(filename * "_zonal_average.jld2", "wb")
vb = FieldTimeSeries(filename * "_zonal_average.jld2", "vb")

fig = Figure(resolution=(1200, 1600))

axbb = Axis(fig[1, 1])
axwb = Axis(fig[2, 1])
axvb = Axis(fig[3, 1])
axN = Axis(fig[4, 1])
axS = Axis(fig[5, 1])

Nt = length(bb.times)

slider = Slider(fig[6, :], range=1:Nt, startvalue=1)
n = slider.value

x, y, z = nodes(bb)

bbn = @lift interior(bb[$n])[1, :, :]
wbn = @lift interior(wb[$n])[1, :, :]
vbn = @lift interior(vb[$n])[1, :, :]
Bn = @lift interior(B[$n])[1, :, :]
Nsq = @lift begin
    Bn = B[$n]
    Nsq = Field(∂z(Bn))
    compute!(Nsq)
    interior(Nsq)[1, :, :]
end

Sh = @lift begin
    Un = U[$n]
    # Sh = Field(sqrt(∂z(Un)^2))
    Vn = V[$n]
    Sh = Field(sqrt(∂z(Un)^2 + ∂z(Vn)^2))
    compute!(Sh)
    interior(Sh)[1, :, :]
end

buoyancycontours!(ax) = contour!(ax, y * 1e-3, z, Bn; levels=25, color = (:black, 0.6), linewidth = 3)

heatmap!(axbb, y * 1e-3, z, bbn, colormap=:speed)
buoyancycontours!(axbb)

heatmap!(axwb, y * 1e-3, z, wbn, colormap=:curl)
buoyancycontours!(axwb)

heatmap!(axvb, y * 1e-3, z, vbn, colormap=:curl)
buoyancycontours!(axvb)

heatmap!(axN, y * 1e-3, z, Nsq, colormap=:thermal, colorrange=(1e-6, 1e-5))
buoyancycontours!(axN)

heatmap!(axS, y * 1e-3, z, Sh, colormap=:solar, colorrange=(0, 1e-4))
buoyancycontours!(axS)

display(fig)

#=
record(fig, "eddy_fluxes.mp4", 1:Nt; framerate=12) do nn
    @info "Drawing frame $nn of $Nt..."
    n[] = nn
end
=#

max_stars_repo_name: navidcy/mesoscale-parametrization

max_issues_repo_name: navidcy/mesoscale-parametrization

max_forks_repo_name: navidcy/mesoscale-parametrization


text:
import AmplNLWriter
import Ipopt
import MINLPTests

using Test

const OPTIMIZER = () -> AmplNLWriter.Optimizer(Ipopt.amplexe, ["print_level=0"])

@testset "MINLPTests" begin
    ###
    ### src/nlp tests.
    ###

    MINLPTests.test_nlp(
        OPTIMIZER,
        exclude = [
            "005_011",  # Uses the function `\`
            "006_010",  # User-defined function
            "007_010",  # Infeasible model
        ],
        objective_tol = 1e-5,
        primal_tol = 1e-5,
        dual_tol = NaN,
    )

    @testset "nlp_007_010" begin
        MINLPTests.nlp_007_010(
            OPTIMIZER,
            1e-5,
            NaN,
            NaN,
            Dict(MINLPTests.INFEASIBLE_PROBLEM => AmplNLWriter.MOI.INFEASIBLE),
            Dict(MINLPTests.INFEASIBLE_PROBLEM => AmplNLWriter.MOI.NO_SOLUTION),
        )
    end

    ###
    ### src/nlp-cvx tests.
    ###

    MINLPTests.test_nlp_cvx(
        OPTIMIZER,
        exclude = [
            "109_010"  # Ipopt fails to converge
        ]
    )
end


text:
"""
This submodule contains useful implementations of the various parts of a probabilistic tuning problem.
"""
module PBTLibrary

"""
Quick callback function to be used in the sciml_train optimization process. Displays current loss.
"""
function basic_pbt_callback(p, loss)
    display(loss)
    return false
end

function benchmark_callback(p, loss, tempt, templ, initial_time)
    if Base.Libc.time() - initial_time  >= 600
        return true
    else
        push!(templ, loss)
        display(loss)
        push!(tempt,Base.Libc.time())
        return false
    end
end

export confidence_interval
function confidence_interval(losses, delta)
    N = length(losses)+4
    d = (sum(losses[:].<delta)+2)/N
    c_high=d+(d*(1-d)/N)^0.5
    c_low=d-(d*(1-d)/N)^0.5
    return d, c_low, c_high
end

end


text:
using SpecialFunctions

abstract type HelmholtzOperator2D <: IntegralOperator end
scalartype(::HelmholtzOperator2D) = ComplexF64

struct SingleLayer{T} <: HelmholtzOperator2D
    wavenumber::T
end

struct HyperSingular{T} <: HelmholtzOperator2D
    wavenumber::T
end

struct DoubleLayer{T} <: HelmholtzOperator2D
    wavenumber::T
end

struct DoubleLayerTransposed{T} <: HelmholtzOperator2D
    wavenumber::T
end

function cellcellinteractions!(biop::HelmholtzOperator2D, tshs, bshs, tcell, bcell, z)

    regularcellcellinteractions!(biop, tshs, bshs, tcell, bcell, z)

end

function testfunc1()
    print("test function!")
end

defaultquadstrat(op::HelmholtzOperator2D, tfs, bfs) = DoubleNumWiltonSauterQStrat(4,3,4,3,4,4,4,4)

function quaddata(op::HelmholtzOperator2D, g::LagrangeRefSpace, f::LagrangeRefSpace, tels, bels,
        qs::DoubleNumWiltonSauterQStrat)

  tqd = quadpoints(g, tels, (qs.outer_rule_far,))
  bqd = quadpoints(f, bels, (qs.inner_rule_far,))

  return (tpoints=tqd, bpoints=bqd)
end

function quadrule(op::HelmholtzOperator2D, g::LagrangeRefSpace, f::LagrangeRefSpace, i, τ, j, σ, qd,
        qs::DoubleNumWiltonSauterQStrat)

    DoubleQuadRule(
        qd.tpoints[1,i],
        qd.bpoints[1,j]
    )

end

mutable struct KernelValsHelmholtz2D
    wavenumber
    vect
    dist
    green
    gradgreen
    txty
end

function kernelvals(biop::HelmholtzOperator2D, tgeo, bgeo)

    k = biop.wavenumber
    r = tgeo.cart - bgeo.cart
    R = norm(r)

    kr = k * R
    hankels = hankelh2.([0 1], kr)
    green = -im / 4 * hankels[1]
    gradgreen = k * im / 4 * hankels[2] * r / R

    txty = dot(normal(tgeo), normal(bgeo))

    KernelValsHelmholtz2D(k, r, R, green, gradgreen, txty)
end

shapevals(op::HelmholtzOperator2D, ϕ, ts) = shapevals(ValDiff(), ϕ, ts)

function integrand(biop::SingleLayer, kerneldata, tvals,
    tgeo, bvals, bgeo)

    gx = tvals[1]
    fy = bvals[1]

    gx * kerneldata.green * fy
end

function integrand(biop::HyperSingular, kernel, tvals, tgeo,
    bvals, bgeo)

    gx = tvals[1]
    fy = bvals[1]

    dgx = tvals[2]
    dfy = bvals[2]

    k    = kernel.wavenumber
    G    = kernel.green
    txty = kernel.txty

    (dgx * dfy - k*k * txty * gx * fy) * G
end

function integrand(biop::DoubleLayer, kernel, fp, mp, fq, mq)
    nq = normal(mq)
    fp[1] * dot(nq, -kernel.gradgreen) * fq[1]
end

function integrand(biop::DoubleLayerTransposed, kernel, fp, mp, fq, mq)
    np = normal(mp)
    fp[1] * dot(np, kernel.gradgreen) * fq[1]
end

mutable struct PlaneWaveDirichlet{T,P} <: Functional
    wavenumber::T
    direction::P
end

mutable struct PlaneWaveNeumann{T,P} <: Functional
    wavenumber::T
    direction::P
end

mutable struct ScalarTrace{F} <: Functional
    f::F
end

strace(f, mesh::Mesh) = ScalarTrace(f)

(s::ScalarTrace)(x) = s.f(cartesian(x))
integrand(s::ScalarTrace, tx, fx) = dot(tx.value, fx)

shapevals(f::Functional, ϕ, ts) = shapevals(ValOnly(), ϕ, ts)

function (field::PlaneWaveDirichlet)(mp)

    wavenumber = field.wavenumber
    direction  = field.direction

    cart = cartesian(mp)
    exp(-im * wavenumber * dot(direction, cart))
end

function(field::PlaneWaveNeumann)(mp)

    wavenumber = field.wavenumber
    direction  = field.direction

    cart = cartesian(mp)
    norm = normal(mp)

    wave = exp(-im * wavenumber * dot(direction, cart))
    grad = -im * wavenumber * direction * wave

    d = norm[1] * grad[1]
    for i in 2:length(norm)  d += norm[i] * grad[i]  end
    return d
end

function integrand(pw::PlaneWaveDirichlet, sv, fx)
    tx = sv[1]
    return dot(tx, fx)
end

function integrand(pw::PlaneWaveNeumann, sv, fx)
    tx = sv[1]
    d = tx[1] * fx[1]
    for i in 2:length(tx) d += tx[i]*fx[i] end
    return d
end


text:
function epitope_feature_matrix(
    data::AbstractHLAData; rank_threshold::Real = 2, allele_depth::Int = 1,
)
    df = epitope_prediction(data, rank_threshold = rank_threshold,
        allele_depth = allele_depth)

    hla_types = Escape.hla_types(data)
    alleles = sort(unique_alleles(hla_types, allele_depth = allele_depth))

    positions = sequence_length(records(data))

    m = zeros(positions, length(alleles))

    for row in eachrow(df)
        allele_idx = findfirst(x -> x == parse_allele(row[:best_allele]), alleles)
        start_idx = row[:pos] + 1
        stop_idx = start_idx + length(row[:peptide]) - 1

        m[start_idx:stop_idx, allele_idx] .= 1
    end

    return m
end

function epitope_prediction(data::AbstractHLAData;
    rank_threshold::Real = 2, allele_depth::Int = 1
)

    records = Escape.records(data)
    consensus = replace(consensus_sequence(records), '-' => 'X')
    alleles = unique_alleles(hla_types(data), allele_depth = allele_depth)

    df = epitope_prediction(consensus, alleles, rank_threshold = rank_threshold)

    return df
end

function epitope_prediction(
    query::String, alleles::Vector{HLAAllele};
    rank_threshold::Real = 100
)
    @assert all(hla_accuracy.(alleles) .== hla_accuracy(alleles[1]))
    accuracy = min(2, hla_accuracy(alleles[1]))

    valid_4_digits = valid_alleles()
    filter!(x -> parse(Int, x.field_2) <= 1, valid_4_digits)

    valid_matching_digits = limit_hla_accuracy.(
        valid_4_digits,
        allele_depth = accuracy
    )

    alleles = valid_4_digits[map(x -> x in alleles, valid_matching_digits)]
    allele_strings = String[]
    output_file = tempname()

    for allele in string.(alleles)
        push!(allele_strings, allele[1:5] * allele[7:end])
    end

    @suppress run(
        `$(Conda.BINDIR)/mhcflurry-predict-scan
            --alleles $allele_strings
            --sequences $query
            --out $output_file
            --results-all`
    )

    df = DataFrame(CSV.File(output_file))
    filter!(x -> (x[:affinity_percentile] <= rank_threshold) |
        (x[:presentation_percentile] <= rank_threshold), df)

    df[!, :best_allele] = Escape.limit_hla_accuracy.(
        parse_allele.(df[!, :best_allele]),
        allele_depth = accuracy
    ) .|> string

    return df
end

function valid_alleles()
    mhcflurry_output = @suppress read(
        `$(Conda.BINDIR)/mhcflurry-predict --list-supported-alleles`, String
    )
    valid_allele_strings = split(mhcflurry_output, '\n')
    filter!(x -> startswith(x, "HLA"), valid_allele_strings)

    valid_alleles = parse_allele.(valid_allele_strings)

    return unique(skipmissing(valid_alleles))
end

function consensus_sequence(records::Vector{FASTX.FASTA.Record})
    consensus = Vector{String}(undef, sequence_length(records))

    for i in 1:sequence_length(records)
        sequence = [string(FASTA.sequence(record)[i]) for record in records]
        counts = countmap(sequence)
        delete!(counts, "-")
        delete!(counts, "X")

        if isempty(counts)
            consensus[i] = "-"
        else
            consensus[i] = findfirst(x -> x == maximum(values(counts)), counts)
        end
    end

    return string.(consensus...)
end

function sequence_length(records::Vector{FASTX.FASTA.Record})
    sequence_length = maximum([length(FASTA.sequence(r)) for r in records])

    return sequence_length
end


text:
function smoothness(lambda::Array{<:Real}, laplacian::Array{<:Real})::Float64
    return lambda' * laplacian * lambda
end

function objective(mu::Array{<:Real}, omega::Array{<:Real}, lambda::Array{<:Real})::Float64
    return omega' * (mu - lambda).^2 / 2
end

max_stars_repo_name: tomaskocak/BAI

max_issues_repo_name: tomaskocak/BAI

max_forks_repo_name: tomaskocak/BAI


text:
#------Distance Matrix Functions--------
"""
    dist_matrix(distr1, distr2; dist_met = kmer_seeded_edit_dist)

distances[i, j] is distance from distr1[i] to distr2[j].
"""
function dist_matrix(distr1, distr2; dist_met = kmer_seeded_edit_dist)
    if distr1 == distr2
        return symmetric_dist_matrix(distr1, dist_met = dist_met)
    end

    distances = zeros(length(distr1), length(distr2))
    for i in 1:length(distr1)
        for j in 1:length(distr2)
            distances[i, j] = dist_met(distr1[i], distr2[j])
        end
    end
    return distances
end

"""
    symmetric_dist_matrix(distr; dist_met = kmer_seeded_edit_dist)

Similar to dist_matrix but where distr1 == distr2.
Automatically called by dist_matrix if the two distributions are identical.
"""
function symmetric_dist_matrix(distr; dist_met = kmer_seeded_edit_dist)
    distances = zeros(length(distr), length(distr))
    for i in 1:length(distr)
        for j in i+1:length(distr)
            distances[i, j] = distances[j, i]  = dist_met(distr[i], distr[j])
        end
    end
    return distances
end

#-------OTHER UTILITY FUNCTIONS----------
"""
    print_fasta(seqs, names)

Prints fasta format to the terminal, for copypasting into alignment/blast etc.
"""
function print_fasta(seqs, names)
    for i in 1:length(seqs)
        println(">", names[i])
        println(seqs[i])
    end
end

"""
    degap(s::String)

Returns given string without '-' gap symbols.
"""
function degap(s::String)
    return replace(s, "-", "")
end

"""
    degap(s::DNASequence)

Returns given string without '-' gap symbols.
"""
function degap(s::DNASequence)
    return DNASequence(degap(String(s)))
end

"""
    dash_count(inStr::String)

Counts number of gap symbols '-' in given string.
"""
function dash_count(inStr::String)
    sum([char == '-' for char in inStr])
end

"""
    single_gap(str::String, ind::Int)

True if `str` has a single gap '-' at index `ind`, else false.
"""
function single_gap(str::String, ind::Int)
    return str[ind] == '-' && (ind == 1 || str[ind-1] != '-') && (ind == length(str) || str[ind+1] != '-')
end

"""
    single_mod_three_gap(str::String, ind::Int)

True if `str` has a gap length of 1 mod 3 at given index.
"""
function single_mod_three_gap(str::String, ind::Int)
    hi = ind
    while hi <= length(str) && str[hi] == '-'
        hi += 1
    end
    lo = ind
    while lo >= 1 && str[lo] == '-'
        lo -= 1
    end
    return (hi - lo - 1)%3 == 1
end

"""
    triple_gap(str::String, ind::Int)

True if index is within a gap of length a multiple of three, else false.
"""
function triple_gap(str::String, ind::Int)
    hi = ind
    while hi <= length(str) && str[hi] == '-'
        hi += 1
    end
    lo = ind
    while lo >= 1 && str[lo] == '-'
        lo -= 1
    end
    return (hi - lo - 1)%3 == 0
end

"""
    seq_details(fasta_path)

Gives names, sequences, error rates, and lengths from given filepath, which may end in '.fasta' or '.fastq'.
"""
function seq_details(fasta_path)
	if last(fasta_path) == 'a'
		names, seqs = read_fasta_with_names(fasta_path)
    	error_rates = [parse(Float64, split(i,['|', '='])[3]) for i in names]
	else
		seqs, errs, names = read_fastq(fasta_path)
		error_rates = [mean(phred_to_p(e)) for e in errs]
    end
    string_lengths = [length(i) for i in seqs]
    return names, seqs, error_rates, string_lengths
end

"""
    print_rgb(r, g, b, t)

Prints in colors `r`,`g`,`b` to terminal.
"""
function print_rgb(r, g, b, t)
    print("\e[1m\e[38;2;$r;$g;$b;249m", t)
end

"""
    reverse_complement(dna_string::String)

Returns the complement of the reverse of given nucleotide sequence.
"""
function reverse_complement(dna_string::String)
    return String(reverse_complement(DNASequence(dna_string)))
end

"""
    print_diffs(s1, s2; width=5, prefix="")

Prints two already aligned sequences with differences in color to terminal.
"""
function print_diffs(s1, s2; width=5, prefix="")
    if length(s1) != length(s2)
        error("Aligned strings are meant to be the same length.")
    end
    for i in 1:length(s1)
        if s1[i] != s2[i]
            println(prefix*"Position: ", i)
            # println(s1[max(i-width, 1):min(i+width, length(s1))])
            print_rgb(0, 0, 0, prefix*s1[max(i-width, 1):max(i-1, 1)])
            print_rgb(117, 41, 48, s1[i])
            print_rgb(0, 0, 0, s1[min(i+1, length(s1)):min(i+width, length(s1))])
            println("")

            print_rgb(0, 0, 0, prefix*s2[max(i-width, 1):max(i-1, 1)])
            print_rgb(117, 41, 48, s2[i])
            print_rgb(0, 0, 0, s2[min(i+1, length(s1)):min(i+width, length(s1))])
            println("")
        end
    end
end

"""
    trim_ends_indices(seq, ref; edge_reduction=0.1)

Align `seq` to `ref` with default low penalties for gaps on ends,
and trim insertions on the ends of `seq`.
Returns (start, stop) indices.
"""
function trim_ends_indices(seq, ref; edge_reduction=0.1)
    saln, raln = nw_align(seq, ref; edge_reduction=edge_reduction)
    x = 1
    while raln[x] == '-'
        x += 1
    end
    y = length(raln)
    while raln[y] == '-'
        y -= 1
    end
    n_end = length(raln) - y
    y = length(seq) - n_end
    return x, y
end

#----Amino Acids---

"""
    translate_to_aa(s::String)

Return amino acid string translation of nucleotide sequence using BioSequences conversion.
"""
function translate_to_aa(s::String)
    rna = convert(RNASequence, DNASequence(s))
    return string(translate(rna))
end

"""
    generate_aa_seqs(str::String)

Return sequence translated to amino acids in each reading frame
(returns three amino acid sequences)."""
function generate_aa_seqs(str::String)
    aa1 = translate_to_aa(str[1:length(str) - length(str)%3])
    aa2 = translate_to_aa(str[2:(length(str) - (length(str) - 1) % 3)])
    aa3 = translate_to_aa(str[3:(length(str) - (length(str) - 2) % 3)])
    return aa1, aa2, aa3
end

#----Length Filtering----

"""
    length_filter_inds(seqs::Array{String, 1}, minlength::Int, maxlength::Int)

Return indices of sequences in given array with length within given range (inclusive).
"""
function length_filter_inds(seqs::Array{String, 1}, minlength::Int, maxlength::Int)
    inds = map(s->(minlength <= length(s) <= maxlength), seqs)
    return inds
end

"""
    length_filter(seqs::Array{String, 1}, phreds::Union{Array{Vector{Phred},1},Void}, names::Union{Array{String,1},Void},
                  minlength::Int, maxlength::Int)

Filter sequences and corresponding names and phreds (which may be `nothing`) by length.
"""
function length_filter(seqs::Array{String, 1}, phreds::Union{Array{Vector{Phred},1},Void}, names::Union{Array{String,1},Void},
                           minlength::Int, maxlength::Int)
    inds = length_filter_inds(seqs, minlength, maxlength)
    seqs = seqs[inds]
    phreds = (phreds != nothing && length(phreds) > 0) ? phreds[inds] : phreds
    names = (names != nothing && length(names) > 0) ? names[inds] : names
    return seqs, phreds, names
end

"""
    length_filter(seqs::Array{String, 1}, minlength::Int64, maxlength::Int64)

Filter sequences by length.
"""
function length_filter(seqs::Array{String, 1}, minlength::Int64, maxlength::Int64)
    inds = length_filter_inds(seqs, minlength, maxlength)
    return seqs[inds]
end

"""
    filter_by_length(args...)

Deprecated. See `length_filter`.
"""
function filter_by_length(args...)
    Base.depwarn("Use length_filter instead.", :filter_by_length)
    return length_filter(args...)
end

"""
    concat_fastas(filepaths::Array{String, 1}, outfile::String)

Write contents of all given files to a single .fasta file.
"""
function concat_fastas(filepaths::Array{String, 1}, outfile::String)
    if length(outfile) < 7 || outfile[end-5:end] != ".fasta"
        error("Outfile must end in .fasta")
    end
    all_seqs, all_names = String[], String[]
    for (j, infile) in enumerate(infiles_list)
        fl = infile
        sqs, nms = read_fasta_with_names_in_other_order(fl, seqtype=String)
        for i in 1:length(sqs)
            push!(all_seqs, sqs[i])
            push!(all_names, nms[i])
        end
    end
    write_fasta(outfile, all_seqs, names=all_names)
end

#------Frequency Fxns------

"""
    maxfreq(vec)

Return the frequency of the most common element in `vec`.
"""
function maxfreq(vec)
    mo = mode(vec)
    proportionmap(vec)[mo]
end

"""
    freq(vec, elem)

Return the frequency of given element in given array;
if the element is not present, return 0.0.
"""
function freq(vec, elem)
    map = proportionmap(vec)
    if haskey(map, elem)
        return map[elem]
        else return 0.0
    end
end

"""
    sorted_freqs(vec)

Return tuples of (freq, elem) of unique elements of `vec`
in order of decreasing frequency.
"""
function sorted_freqs(vec)
    propDict = proportionmap(vec)
    seqkeys = keys(propDict)
    return reverse(sort([(propDict[k],k) for k in seqkeys]))
end

"""
    freq_dict_print(dictin; thresh=0)

Prints frequency:element of elements of `dictin` above given threshold,
where `dictin` is a proportionmap of elements (see `proportionmap` in StatsBase).
"""
function freq_dict_print(dictin; thresh=0)
    for i in keys(dictin)
        if dictin[i]>thresh
            println(dictin[i], ":", i)
        end
    end
end

"""
    logsum(lga, lgb)

Compute numerically stable logsum. Returns `-Inf` if either input is `-Inf`.
"""
function logsum(lga, lgb)
    if lga == lgb == -Inf
        return -Inf
    end
    # log(x) + log(1 + exp(log(y) - log(x))) for x > y
    return lga > lgb ? (lga + log(1 + exp(lgb - lga))) : (lgb + log(1 + exp(lga - lgb)))
end

"""
	get_group_inds(arr, func)

Apply a function to a given set of objects, group by equal outputs, and return indices of elements in each group
"""
function get_group_inds(arr,func)
   groupings = [func(i) for i in arr]
   groups = union(groupings)
   groupinds = [find(groupings.==i) for i in groups]
   return groupinds
end

const nl43env = "GAGCAGAAGACAGTGGCAATGAGAGTGAAGGAGAAGTATCAGCACTTGTGGAGATGGGGGTGGAAA"*
"TGGGGCACCATGCTCCTTGGGATATTGATGATCTGTAGTGCTACAGAAAAATTGTGGGTCACAGTCTATTATGGGGT"*
"ACCTGTGTGGAAGGAAGCAACCACCACTCTATTTTGTGCATCAGATGCTAAAGCATATGATACAGAGGTACATAATG"*
"TTTGGGCCACACATGCCTGTGTACCCACAGACCCCAACCCACAAGAAGTAGTATTGGTAAATGTGACAGAAAATTTT"*
"AACATGTGGAAAAATGACATGGTAGAACAGATGCATGAGGATATAATCAGTTTATGGGATCAAAGCCTAAAGCCATG"*
"TGTAAAATTAACCCCACTCTGTGTTAGTTTAAAGTGCACTGATTTGAAGAATGATACTAATACCAATAGTAGTAGCG"*
"GGAGAATGATAATGGAGAAAGGAGAGATAAAAAACTGCTCTTTCAATATCAGCACAAGCATAAGAGATAAGGTGCAG"*
"AAAGAATATGCATTCTTTTATAAACTTGATATAGTACCAATAGATAATACCAGCTATAGGTTGATAAGTTGTAACAC"*
"CTCAGTCATTACACAGGCCTGTCCAAAGGTATCCTTTGAGCCAATTCCCATACATTATTGTGCCCCGGCTGGTTTTG"*
"CGATTCTAAAATGTAATAATAAGACGTTCAATGGAACAGGACCATGTACAAATGTCAGCACAGTACAATGTACACAT"*
"GGAATCAGGCCAGTAGTATCAACTCAACTGCTGTTAAATGGCAGTCTAGCAGAAGAAGATGTAGTAATTAGATCTGC"*
"CAATTTCACAGACAATGCTAAAACCATAATAGTACAGCTGAACACATCTGTAGAAATTAATTGTACAAGACCCAACA"*
"ACAATACAAGAAAAAGTATCCGTATCCAGAGGGGACCAGGGAGAGCATTTGTTACAATAGGAAAAATAGGAAATATG"*
"AGACAAGCACATTGTAACATTAGTAGAGCAAAATGGAATGCCACTTTAAAACAGATAGCTAGCAAATTAAGAGAACA"*
"ATTTGGAAATAATAAAACAATAATCTTTAAGCAATCCTCAGGAGGGGACCCAGAAATTGTAACGCACAGTTTTAATT"*
"GTGGAGGGGAATTTTTCTACTGTAATTCAACACAACTGTTTAATAGTACTTGGTTTAATAGTACTTGGAGTACTGAA"*
"GGGTCAAATAACACTGAAGGAAGTGACACAATCACACTCCCATGCAGAATAAAACAATTTATAAACATGTGGCAGGA"*
"AGTAGGAAAAGCAATGTATGCCCCTCCCATCAGTGGACAAATTAGATGTTCATCAAATATTACTGGGCTGCTATTAA"*
"CAAGAGATGGTGGTAATAACAACAATGGGTCCGAGATCTTCAGACCTGGAGGAGGCGATATGAGGGACAATTGGAGA"*
"AGTGAATTATATAAATATAAAGTAGTAAAAATTGAACCATTAGGAGTAGCACCCACCAAGGCAAAGAGAAGAGTGGT"*
"GCAGAGAGAAAAAAGAGCAGTGGGAATAGGAGCTTTGTTCCTTGGGTTCTTGGGAGCAGCAGGAAGCACTATGGGCG"*
"CAGCGTCAATGACGCTGACGGTACAGGCCAGACAATTATTGTCTGATATAGTGCAGCAGCAGAACAATTTGCTGAGG"*
"GCTATTGAGGCGCAACAGCATCTGTTGCAACTCACAGTCTGGGGCATCAAACAGCTCCAGGCAAGAATCCTGGCTGT"*
"GGAAAGATACCTAAAGGATCAACAGCTCCTGGGGATTTGGGGTTGCTCTGGAAAACTCATTTGCACCACTGCTGTGC"*
"CTTGGAATGCTAGTTGGAGTAATAAATCTCTGGAACAGATTTGGAATAACATGACCTGGATGGAGTGGGACAGAGAA"*
"ATTAACAATTACACAAGCTTAATACACTCCTTAATTGAAGAATCGCAAAACCAGCAAGAAAAGAATGAACAAGAATT"*
"ATTGGAATTAGATAAATGGGCAAGTTTGTGGAATTGGTTTAACATAACAAATTGGCTGTGGTATATAAAATTATTCA"*
"TAATGATAGTAGGAGGCTTGGTAGGTTTAAGAATAGTTTTTGCTGTACTTTCTATAGTGAATAGAGTTAGGCAGGGA"*
"TATTCACCATTATCGTTTCAGACCCACCTCCCAATCCCGAGGGGACCCGACAGGCCCGAAGGAATAGAAGAAGAAGG"*
"TGGAGAGAGAGACAGAGACAGATCCATTCGATTAGTGAACGGATCCTTAGCACTTATCTGGGACGATCTGCGGAGCC"*
"TGTGCCTCTTCAGCTACCACCGCTTGAGAGACTTACTCTTGATTGTAACGAGGATTGTGGAACTTCTGGGACGCAGG"*
"GGGTGGGAAGCCCTCAAATATTGGTGGAATCTCCTACAGTATTGGAGTCAGGAACTAAAGAATAGTGCTGTTAACTT"*
"GCTCAATGCCACAGCCATAGCAGTAGCTGAGGGGACAGATAGGGTTATAGAAGTATTACAAGCAGCTTATAGAGCTA"*
"TTCGCCACATACCTAGAAGAATAAGACAGGGCTTGGAAAGGATTTTGCTATAAGATGGGTGGCAAGTGG";


text:
# Licensed to the Apache Software Foundation (ASF) under one
# or more contributor license agreements.  See the NOTICE file
# distributed with this work for additional information
# regarding copyright ownership.  The ASF licenses this file
# to you under the Apache License, Version 2.0 (the
# "License"); you may not use this file except in compliance
# with the License.  You may obtain a copy of the License at
#
#   http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing,
# software distributed under the License is distributed on an
# "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
# KIND, either express or implied.  See the License for the
# specific language governing permissions and limitations
# under the License.

# activation functions

@doc doc"""
    σ.(x::NDArray)
    sigmoid.(x::NDArray)

Computes sigmoid of x element-wise.

```math
σ(x) = \frac{1}{(1 + exp(-x))}
```

The storage type of `sigmoid` output is always dense.
"""
function σ end
const sigmoid = σ
_nddoc[:σ] = false
@_remap broadcasted(::typeof(σ), x::NDArray) sigmoid(x)

@doc doc"""
    relu.(x::NDArray)

Computes rectified linear.

```math
\max(x, 0)
```
"""
function relu end
_nddoc[:relu] = false
@_remap broadcasted(::typeof(relu), x::NDArray) relu(x)

@doc doc"""
    softmax.(x::NDArray, [dim = ndims(x)])

Applies the softmax function.

The resulting array contains elements in the range `(0, 1)`
and the elements along the given axis sum up to 1.

```math
softmax(\mathbf{z})_j = \frac{e^{z_j}}{\sum_{k=1}^K e^{z_k}}
```
"""
function softmax end
_nddoc[:softmax] = false
@_remap broadcasted(::typeof(softmax), x::NDArray)           softmax(x; axis = -ndims(x))
@_remap broadcasted(::typeof(softmax), x::NDArray, dim::Int) softmax(x; axis = -dim)

"""
    log_softmax.(x::NDArray, [dim = ndims(x)])

Computes the log softmax of the input.
This is equivalent to computing softmax followed by log.

julia> x
2×3 mx.NDArray{Float64,2} @ CPU0:
 1.0  2.0  0.1
 0.1  2.0  1.0

julia> mx.log_softmax.(x)
2×3 mx.NDArray{Float64,2} @ CPU0:
 -1.41703  -0.41703  -2.31703
 -2.31703  -0.41703  -1.41703
"""
function log_softmax end
_nddoc[:log_softmax] = false
@_remap broadcasted(::typeof(log_softmax), x::NDArray)           log_softmax(x; axis = -ndims(x))
@_remap broadcasted(::typeof(log_softmax), x::NDArray, dim::Int) log_softmax(x; axis = -dim)


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule CUTENSOR_CUDA110_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("CUTENSOR_CUDA110")
JLLWrappers.@generate_main_file("CUTENSOR_CUDA110", UUID("004fe424-f3b5-51a0-a814-026e7c39908e"))
end  # module CUTENSOR_CUDA110_jll


text:
using Levenshtein
using Base.Test
using Compat

function levenshtein_base(a::AbstractString, b::AbstractString, deletion_cost::Real, insertion_cost::Real, substitution_cost::Real)
    local costs::Array{Real, 2} = zeros(Real, length(a) + 1, length(b) + 1)
    # Initialize the deletion costs
    for i in 1:length(a)
        costs[i + 1,1] = i * deletion_cost
    end
    # Initialize the insertion costs
    for i in 1:length(b)
        costs[1, i + 1] = i * insertion_cost
    end

    local row = 2
    for r in a
        local col = 2
        for c in b
            local deletion::Real = costs[row - 1, col] + deletion_cost
            local insertion::Real = costs[row, col - 1] + insertion_cost
            local substitution::Real = costs[row - 1, col - 1] + (r == c ? 0 : substitution_cost)
            costs[row, col] = min(deletion, insertion, substitution)
            col += 1
        end
        row += 1
    end

    return costs[length(a) + 1, length(b) + 1]
end

# Testing base implementation
# Tests checked with http://www.kurzhals.info/static/samples/levenshtein_distance/
@test 4 == levenshtein_base("Hi, my name is", "my name is", 1, 1, 1)
@test 21 == levenshtein_base("%^@!^@#^@#!! Snoooooooop", "Dro!p it!!!! like it's hot", 1, 1, 1)
@test 7 == levenshtein_base("Alborgów", "amoniak", 1, 1, 1)

@test 3 == levenshtein_base("sitting", "kitten", 1, 1, 1)
@test 1 == levenshtein_base("sitting", "kitten", 1, 0, 0)
@test 0 == levenshtein_base("sitting", "kitten", 0, 1, 0)
@test 0 == levenshtein_base("sitting", "kitten", 0, 0, 1)
@test 5 == levenshtein_base("sitting", "kitten", 1, 1, Inf)
@test 12 == levenshtein_base("sitting", "kitten", 10, 1, 1)

# write your own tests here
@test 4 == levenshtein("Hi, my name is", "my name is")
@test 21 == levenshtein("%^@!^@#^@#!! Snoooooooop", "Dro!p it!!!! like it's hot")
@test 7 == levenshtein("Alborgów", "amoniak")

@test 3 == levenshtein("sitting", "kitten")
@test 1 == levenshtein("sitting", "kitten", 1, 0, 0)
@test 0 == levenshtein("sitting", "kitten", 0, 1, 0)
@test 0 == levenshtein("sitting", "kitten", 0, 0, 1)
@test 5 == levenshtein("sitting", "kitten", 1, 1, Inf)
@test 12 == levenshtein("sitting", "kitten", 10, 1, 1)

@test levenshtein("a", "ab", 1, 2, 0) == levenshtein("ab", "a", 2, 1, 0)

@test levenshtein("a", "", 10, 0, 0) == 10
@test levenshtein("", "a", 0, 10, 0) == 10

for i in 1:100
    s = randstring(floor(Int,rand() * 100))
    t = randstring(floor(Int,rand() * 100))
    ins, del, sub = rand() * 1000, rand() * 1000, rand() * 1000

    calculated, truth = levenshtein(s, t, del, ins, sub), levenshtein_base(s, t, del, ins, sub)

    if calculated != truth
        println("s = $s, t = $t, ins = $ins, del = $del, sub = $sub, calculated = $calculated, truth = $truth")
    end

    @test abs(calculated - truth) < 1e-9
end


text:
using Photon, Printf
using Knet:value
using Photon.Layers: Dense
using Photon.Losses: MSELoss

# An example that demonstrates how to use Photon for linear regression problem

"""
The most simple model possible, just a single layer with 1 output
"""
model = Dense(1)

"""
This is the function we want to approximate
"""
fn(x,y,z) = 2x - 3y + 4z .- 5

"""
Now we create the training data
"""
X = [randn(Float32, 3, 8) for _ in 1:50]
Y = [fn(x[1,:], x[2,:], x[3,:]) for x in X]

"""
And finally we create the workout and start the training.
We use Mean Square Error loss and SGD as the optimizer (default).
"""
workout = Workout(model, MSELoss())
train!(workout, zip(X,Y), epochs=5)

# let's see how well we approximated the function
@printf "%.1fx + %.1fy  + %.1fz + %.1f\n" Array(value(model.params.w))... Array(value(model.params.b))...


text:
export FixedRange, FixedEnd, @fixed_range

""" `FixedRange{A, B}()` is like `UnitRange{Int64}(A, B)`, but the bounds are encoded
in the type. """
struct FixedRange{A, B} end
""" `FixedRange{2, FixedEnd()}()` behaves like a type-stable 2:end """
struct FixedEnd{N} end
FixedEnd() = FixedEnd{0}()

Unrolled.type_length(::Type{FixedRange{A, B}}) where {A, B} = B - A + 1
Base.length(fr::FixedRange{A, B}) where {A, B} = B - A + 1
Base.maximum(fr::FixedRange{A, B}) where {A, B} = B::Int
Base.minimum(fr::FixedRange{A, B}) where {A, B} = A::Int
Base.getindex(fr::FixedRange, i::Int) = minimum(fr) + i - 1

Base.iterate(fr::FixedRange, state=minimum(fr)) =
    state > maximum(fr) ? nothing : (state, state+1)
Base.lastindex(fr::FixedRange) = length(fr)

replace_end(::FixedEnd{N}, ::Type{SEQ}) where {N, SEQ} = type_length(SEQ) - N
replace_end(n::Int, ::Type) = n

@generated Base.getindex(seq, ::FixedRange{A, B}) where {A, B} =
    :(tuple($((:(seq[$i]) for i in replace_end(A, seq):replace_end(B, seq))...)))

""" `@fixed_range 3:10` behaves like the standard range `3:10`, but is stored within
the type system, so that `some_tuple[@fixed_range 3:10]` is type-stable. Also supports
`some_tuple[@fixed_range 3:end-5]` """
macro fixed_range(r::Expr)
    process(x::Int) = x
    process(x::Symbol) = x === :end ? :($Unrolled.FixedEnd()) : x
    function process(x::Expr)
        @assert @capture(x, en_-m_) "`@fixed_range` macro cannot handle $x"
        @assert en === :end
        :(FixedEnd{$m}())
    end
    expand(a, b) =
        :($Unrolled.FixedRange{$(process(a)), $(process(b))}())
    @match r begin
        a_:b_ => esc(expand(a, b))
        s_[a_:b_] => esc(:($s[$(expand(a, b))]))
        any_ => error("Bad @fixed_range")
    end
end


text:
module TensorFunctions

using LinearAlgebra,TensorOperations

export @tensorfunc

#= global setting =#
function order(ex::Expr,ord::Expr)
    if ord.head == :tuple
        QuoteNode[ord.args...]
    else
        error("not implemented")
    end
end
function order(ex::Expr)
    [nothing,]
end

tracefunc=tensortrace
contractfunc=tensorcontract

include("tensorfunc.jl")

end # module


text:
# No scattering in either the added layer or the composite layer.
function interaction_helper!(::ScatteringInterface_00, SFI,
                                composite_layer::CompositeLayer{FT},
                                added_layer::AddedLayer{FT},
                                I_static::AbstractArray{FT2}) where {FT<:Union{AbstractFloat, ForwardDiff.Dual},FT2}

    # If SFI, interact source function in no scattering
    if SFI
        J₀⁺, J₀⁻ = similar(composite_layer.J₀⁺), similar(composite_layer.J₀⁺)
        J₀⁺ = added_layer.J₀⁺ .+ added_layer.t⁺⁺ ⊠ composite_layer.J₀⁺
        J₀⁻ = composite_layer.J₀⁻ .+ composite_layer.T⁻⁻ ⊠ added_layer.J₀⁻
        composite_layer.J₀⁺ = J₀⁺
        composite_layer.J₀⁻ = J₀⁻
    end

    # Batched multiplication between added and composite
    composite_layer.T⁻⁻[:] = added_layer.t⁻⁻ ⊠ composite_layer.T⁻⁻
    composite_layer.T⁺⁺[:] = added_layer.t⁺⁺ ⊠ composite_layer.T⁺⁺
end

# No scattering in inhomogeneous composite layer.
# Scattering in homogeneous layer, added to bottom of the composite layer.
# Produces a new, scattering composite layer.
function interaction_helper!(::ScatteringInterface_01, SFI,
                                composite_layer::CompositeLayer{FT},
                                added_layer::AddedLayer{FT},
                                I_static::AbstractArray{FT2}) where {FT<:Union{AbstractFloat, ForwardDiff.Dual},FT2}

    if SFI
        J₀⁺, J₀⁻ = similar(composite_layer.J₀⁺), similar(composite_layer.J₀⁺)
        J₀⁻ = composite_layer.J₀⁻ .+ composite_layer.T⁻⁻ ⊠ (added_layer.r⁻⁺ ⊠ composite_layer.J₀⁺ .+ added_layer.J₀⁻)
        J₀⁺ = added_layer.J₀⁺ .+ added_layer.t⁺⁺ ⊠ composite_layer.J₀⁺
        composite_layer.J₀⁺ = J₀⁺
        composite_layer.J₀⁻ = J₀⁻
    end

    # Batched multiplication between added and composite
    composite_layer.R⁻⁺[:] = composite_layer.T⁻⁻ ⊠ added_layer.r⁻⁺ ⊠ composite_layer.T⁺⁺
    composite_layer.R⁺⁻[:] = added_layer.r⁺⁻
    composite_layer.T⁺⁺[:] = added_layer.t⁺⁺ ⊠ composite_layer.T⁺⁺
    composite_layer.T⁻⁻[:] = composite_layer.T⁻⁻ ⊠ added_layer.t⁻⁻
end

# Scattering in inhomogeneous composite layer.
# no scattering in homogeneous layer which is
# added to the bottom of the composite layer.
# Produces a new, scattering composite layer.
function interaction_helper!(::ScatteringInterface_10, SFI,
                                composite_layer::CompositeLayer{FT},
                                added_layer::AddedLayer{FT},
                                I_static::AbstractArray{FT2}) where {FT<:Union{AbstractFloat, ForwardDiff.Dual},FT2}

    if SFI
        J₀⁺, J₀⁻ = similar(composite_layer.J₀⁺), similar(composite_layer.J₀⁺)
        J₀⁺ = added_layer.J₀⁺ .+ added_layer.t⁺⁺ ⊠ (composite_layer.J₀⁺ .+ composite_layer.R⁺⁻ ⊠ added_layer.J₀⁻)
        J₀⁻ = composite_layer.J₀⁻ .+ composite_layer.T⁻⁻ ⊠ added_layer.J₀⁻
        composite_layer.J₀⁺ = J₀⁺
        composite_layer.J₀⁻ = J₀⁻
    end

    # Batched multiplication between added and composite
    composite_layer.T⁺⁺[:] = added_layer.t⁺⁺ ⊠ composite_layer.T⁺⁺
    composite_layer.T⁻⁻[:] = composite_layer.T⁻⁻ ⊠ added_layer.t⁻⁻
    composite_layer.R⁺⁻[:] = added_layer.t⁺⁺ ⊠ composite_layer.R⁺⁻ ⊠ added_layer.t⁻⁻
end

# Scattering in inhomogeneous composite layer.
# Scattering in homogeneous layer which is added to the bottom of the composite layer.
# Produces a new, scattering composite layer.
function interaction_helper!(::ScatteringInterface_11, SFI,
                                composite_layer::CompositeLayer{FT},
                                added_layer::AddedLayer{FT},
                                I_static::AbstractArray{FT2}) where {FT<:Union{AbstractFloat, ForwardDiff.Dual},FT2}

    @unpack r⁺⁻, r⁻⁺, t⁻⁻, t⁺⁺ = added_layer
    @unpack R⁻⁺, R⁺⁻, T⁺⁺, T⁻⁻, J₀⁺, J₀⁻ = composite_layer

    # Used to store `(I - R⁺⁻ * r⁻⁺)⁻¹`
    tmp_inv = similar(t⁺⁺)

    # Compute and store `(I - R⁺⁻ * r⁻⁺)⁻¹`
    @timeit "interaction inv1" batch_inv!(tmp_inv, I_static .- r⁻⁺ ⊠ R⁺⁻) #Suniti
    # Temporary arrays:
    # T₁₂(I-R₀₁R₂₁)⁻¹
    T01_inv = T⁻⁻ ⊠ tmp_inv;

    # R₂₀ = R₁₀ + T₀₁(I-R₂₁R₀₁)⁻¹ R₂₁T₁₀
    composite_layer.R⁻⁺[:] = R⁻⁺ .+ T01_inv ⊠ r⁻⁺ ⊠ T⁺⁺ #Suniti
    # T₀₂ = T₀₁(1-R₂₁R₀₁)⁻¹T₁₂
    composite_layer.T⁻⁻[:] = T01_inv ⊠ t⁻⁻ #Suniti

    if SFI
        #J₀₂⁻ = J₀₁⁻ + T₀₁(1-R₂₁R₀₁)⁻¹(R₂₁J₁₀⁺+J₁₂⁻)
        composite_layer.J₀⁻[:] = J₀⁻ .+ T01_inv ⊠ (r⁻⁺ ⊠ J₀⁺ .+ added_layer.J₀⁻)
    end

    # Repeating for mirror-reflected directions

    # Compute and store `(I - r⁻⁺ * R⁺⁻)⁻¹`
    @timeit "interaction inv2" batch_inv!(tmp_inv, I_static .- R⁺⁻ ⊠ r⁻⁺) #Suniti
    # T₂₁(I-R₀₁R₂₁)⁻¹
    T21_inv = t⁺⁺ ⊠ tmp_inv

    # T₂₀ = T₂₁(I-R₀₁R₂₁)⁻¹T₁₀
    composite_layer.T⁺⁺[:] = T21_inv  ⊠ T⁺⁺ #Suniti
    # R₀₂ = R₁₂ + T₂₁(1-R₀₁R₂₁)⁻¹R₀₁T₁₂
    composite_layer.R⁺⁻[:] = r⁺⁻ .+ T21_inv ⊠ R⁺⁻ ⊠ t⁻⁻ #Suniti

    if SFI
        # J₂₀⁺ = J₂₁⁺ + T₂₁(I-R₀₁R₂₁)⁻¹(J₁₀ + R₀₁J₁₂⁻ )
        composite_layer.J₀⁺[:] = added_layer.J₀⁺ .+ T21_inv ⊠ (J₀⁺ + R⁺⁻ ⊠ added_layer.J₀⁻)
    end
end

"Compute interaction between composite and added layers"
function interaction!(scattering_interface::AbstractScatteringInterface, SFI,
                        composite_layer::CompositeLayer{FT},
                        added_layer::AddedLayer{FT},
                        I_static::AbstractArray{FT2}) where {FT<:Union{AbstractFloat, ForwardDiff.Dual},FT2}

    interaction_helper!(scattering_interface, SFI, composite_layer, added_layer, I_static)
    synchronize_if_gpu()

end


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

mutable struct IoK8sApiAutoscalingV2beta2ExternalMetricSource <: SwaggerModel
    metric::Any # spec type: Union{ Nothing, IoK8sApiAutoscalingV2beta2MetricIdentifier } # spec name: metric
    target::Any # spec type: Union{ Nothing, IoK8sApiAutoscalingV2beta2MetricTarget } # spec name: target

    function IoK8sApiAutoscalingV2beta2ExternalMetricSource(;metric=nothing, target=nothing)
        o = new()
        validate_property(IoK8sApiAutoscalingV2beta2ExternalMetricSource, Symbol("metric"), metric)
        setfield!(o, Symbol("metric"), metric)
        validate_property(IoK8sApiAutoscalingV2beta2ExternalMetricSource, Symbol("target"), target)
        setfield!(o, Symbol("target"), target)
        o
    end
end # type IoK8sApiAutoscalingV2beta2ExternalMetricSource

const _property_map_IoK8sApiAutoscalingV2beta2ExternalMetricSource = Dict{Symbol,Symbol}(Symbol("metric")=>Symbol("metric"), Symbol("target")=>Symbol("target"))
const _property_types_IoK8sApiAutoscalingV2beta2ExternalMetricSource = Dict{Symbol,String}(Symbol("metric")=>"IoK8sApiAutoscalingV2beta2MetricIdentifier", Symbol("target")=>"IoK8sApiAutoscalingV2beta2MetricTarget")
Base.propertynames(::Type{ IoK8sApiAutoscalingV2beta2ExternalMetricSource }) = collect(keys(_property_map_IoK8sApiAutoscalingV2beta2ExternalMetricSource))
Swagger.property_type(::Type{ IoK8sApiAutoscalingV2beta2ExternalMetricSource }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiAutoscalingV2beta2ExternalMetricSource[name]))}
Swagger.field_name(::Type{ IoK8sApiAutoscalingV2beta2ExternalMetricSource }, property_name::Symbol) =  _property_map_IoK8sApiAutoscalingV2beta2ExternalMetricSource[property_name]

function check_required(o::IoK8sApiAutoscalingV2beta2ExternalMetricSource)
    (getproperty(o, Symbol("metric")) === nothing) && (return false)
    (getproperty(o, Symbol("target")) === nothing) && (return false)
    true
end

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


text:
using Yao, CuYao, CUDA
using BenchmarkTools

reg = rand_state(12; nbatch=1000)
creg = reg |> cu
@benchmark CUDA.@sync creg |> put(12, 3=>Z)
@benchmark CUDA.@sync creg |> put(12, 3=>X)
@benchmark reg |> put(12, 3=>Z)
@benchmark CUDA.@sync creg |> control(12, 6, 3=>X)
@benchmark reg |> control(12, 6, 3=>X)
@benchmark CUDA.@sync creg |> put(12, 3=>rot(X, 0.3))
@benchmark reg |> put(12, 3=>rot(X, 0.3))

reg = rand_state(20)
creg = reg |> cu
g = swap(12, 7, 2)
@benchmark reg |> g
@benchmark CUDA.@sync creg |> g


text:
struct InfeasibleSearchProblemError <: Exception
    msg::String
end

function generic_graph_search(p::SearchProblem, fringe_priority::Function)
    # the closed set, states that we don't need to expand anymore
    closed_set = Set{state_type(p)}()
    # the fringe is a priority queue of SearchNode that are left to be expanded
    fringe = PriorityQueue{search_node_type(p), Float64}()

    n0 = root_node(p)
    enqueue!(fringe, n0, fringe_priority(n0))
    while true
        if isempty(fringe)
            throw(InfeasibleSearchProblemError("Fringe was empty, but no solution found"))
        end
        current_search_node = dequeue!(fringe)
        # We have found a path to the goal state
        if is_goal_state(p, leaf_state(current_search_node))
            return action_sequence(current_search_node), state_sequence(current_search_node), cost(current_search_node)
        elseif !(leaf_state(current_search_node) in closed_set)
            # make sure we don't explore this state again
            push!(closed_set, leaf_state(current_search_node))
            # expand the node
            for child_search_node in expand(current_search_node, p)
                enqueue!(fringe, child_search_node, fringe_priority(child_search_node))
            end
        else
            continue
        end
    end
end

astar_search(p::SearchProblem, heuristic::Function) = weighted_astar_search(p, heuristic, 0.0)

function weighted_astar_search(p::SearchProblem, heuristic::Function, eps::Float64)
    @assert eps >= 0.0
    astar_priority = (n::SearchNode) -> cost(n) + (1+eps)*heuristic(leaf_state(n))
    return generic_graph_search(p, astar_priority)
end


text:
using Revise
using SerialCOBS
using ProtoBuf
using Printf
using Test
using BenchmarkTools

include("$(@__DIR__)/msgs/imu_msg_pb.jl")
imu = IMU(acc_x=0., acc_y=0., acc_z=0.,
          gyr_x=0., gyr_y=0., gyr_z=0.,
          time=time())
in_vicon = VICON(pos_x=0., pos_y=0., pos_z=0.,
                 quat_w=0., quat_x=0., quat_y=0., quat_z=0.,
                 time=time())
imu_vicon = IMU_VICON(imu=imu, vicon=in_vicon)

ard = Arduino("/dev/tty.usbmodem92225501", 57600);

out_vicon = VICON(pos_x=0., pos_y=0., pos_z=0.,
                  quat_w=0., quat_x=0., quat_y=0., quat_z=0.,
                  time=time())
holybro = Arduino("/dev/tty.usbmodem14201", 57600);

# %% Benchmark recieve
open(ard)
results = @benchmark recieve!($ard, $imu_vicon)
close(ard)
display(results)

# %% Benchmark message
open(holybro)

# %%
message(holybro, out_vicon)

# %%
try
    for i in 1:1000
        message(holybro, out_vicon)
        sleep(0.0001)
    end
catch
    close(holybro)
end

# %%

open(holybro) do sp
    message(ard, out_vicon)
end
# close(holybro)

# %%
decoded_msg = 0
encoded_msg = 0
cnt = 0

open(ard)
    start_time = time()
# results = @benchmark begin
    for i in 1:100000
        encoded_msg = SerialCOBS.retrieve_encoded_msg(ard)
        if encoded_msg !== nothing
            decoded_msg = decode(ard, encoded_msg)
            readproto(IOBuffer(decoded_msg), imu_vicon)

            cnt += 1
            # @printf("IMU accel: \t[%1.3f, %1.3f, %1.3f]\n",
            #         imu_vicon.imu.acc_x, imu_vicon.imu.acc_y, imu_vicon.imu.acc_z)
            # @printf("Vicon pos: \t[%1.3f, %1.3f, %1.3f]\n",
            #         imu_vicon.vicon.pos_x, imu_vicon.vicon.pos_x, imu_vicon.vicon.pos_x)
        end
    end
end_time = time()

# end
close(ard)

println(cnt/(end_time - start_time))

# display(results)

# %% Begin Testing here
@testset "Encode/Decode Are Inverse Functions" begin
    test_message1 = rand(UInt8, 25)
    ard = Main.SerialCOBS.Arduino("/dev/tty.usbmodem14201", 57600)

    encoded_msg1 = Main.SerialCOBS.encode(ard, test_message1)
    decoded_msg1 = Main.SerialCOBS.decode(ard, encoded_msg1)
    @test decoded_msg1 == test_message1

    # Large message size
    test_message2 = rand(UInt8, 254)
    ard = Main.SerialCOBS.Arduino("/dev/tty.usbmodem14201", 57600)

    encoded_msg2 = Main.SerialCOBS.encode(ard, test_message2)
    decoded_msg2 = Main.SerialCOBS.decode(ard, encoded_msg2)
    @test decoded_msg2 == test_message2
end

# %%
ard = Arduino(portname, holybro_baudrate)

# %%
@testset "Opening/Closing Arduino Port" begin
    @test !isopen(ard)

    open(ard)
    @test isopen(ard)

    close(ard)
    @test !isopen(ard)
end

# %% Check for Strings
@testset "Arduino Echo" begin
    # Send arduino the byte array encoding the string "Test"
    input = Vector{UInt8}("Test")
    output = 0

    open(ard) do sp
        message(ard, input)
        for i in 1:1000
            # Wait for the arudino to echo back with the reversed message
            if bytesavailable(ard) > 0
                output = recieve(ard)
                break
            end
            sleep(0.01)
            i == 1000 ? error("Is arduino running PacketSerialEcho.ino connected to port $(ard.portname)?") : nothing
        end

    end
    @test isequal(input, output)

    # Do same test but with different Numeric Types
    input = Vector{Float64}([3., 5., 5.])
    output = 0

    open(ard) do sp
        message(ard, Vector(reinterpret(UInt8, input)))

        while true
            if bytesavailable(ard) > 0
                output = reinterpret(Float64, recieve(ard))
                break
            end
        end
    end
    @test input == output
end


text:
macro warnpcfail(ex::Expr)
    modl = __module__
    file = __source__.file === nothing ? "?" : String(__source__.file)
    line = __source__.line
    quote
        $(esc(ex)) || @warn """precompile directive
     $($(Expr(:quote, ex)))
 failed. Please report an issue in $($modl) (after checking for duplicates) or remove this directive.""" _file=$file _line=$line
    end
end

function _precompile_()
    for ST in (N0f8, Float32, Float64)
        for T in (Gray{ST}, RGB{ST})
            @warnpcfail precompile(restrict, (Vector{T},))
            @warnpcfail precompile(restrict, (Matrix{T},))
        end
    end
end # _precompile_


text:
# Utils for modifying plot style variables

import .LibCImPlot: GetStyle, PopStyleVar, PopStyleColor

function PushStyleColor(idx::ImPlotCol_, col::UInt32)
    LibCImPlot.PushStyleColorU32(idx, col)
end

function PushStyleColor(idx::ImPlotCol_, col::ImVec4)
    LibCImPlot.PushStyleColorVec4(idx, col)
end

function PushStyleVar(idx::ImPlotStyleVar_, val::AbstractFloat)
    LibCImPlot.PushStyleVarFloat(idx, Float32(val))
end

function PushStyleVar(idx::ImPlotStyleVar_, val::Integer)
    LibCImPlot.PushStyleVarInt(idx, Cint(val))
end

function PushStyleVar(idx::ImPlotStyleVar_, val1::AbstractFloat, val2::AbstractFloat)
    LibCImPlot.PushStyleVarVec2(idx, ImVec2(Float32(val1), Float32(val2)))
end


text:
using Test
using Spandex

@testset "push! 1" begin
    local pq = PriorityQueue{Int64}(min)

    push!(pq, 3)
    @test 3 == top(pq)
    @test 1 == height(pq)

    push!(pq, 4)
    @test 3 == top(pq)
    @test 2 == height(pq)

    push!(pq, 5)
    @test 3 == top(pq)
    @test 2 == height(pq)

    push!(pq, 2)
    @test 2 == top(pq)
    @test 3 == height(pq)

    push!(pq, 1)
    @test 1 == top(pq)
    @test 3 == height(pq)
end

@testset "push! 2" begin
    local pq = PriorityQueue{Int64}(min)

    for i = 9:-1:1
        push!(pq, i)
    end

    @test 1 == top(pq)
    @test 4 == height(pq)

    push!(pq, 9)
    @test 1 == top(pq)
    @test 4 == height(pq)
end

@testset "push! 3" begin
    local pq = PriorityQueue{Int64}(min, 5)

    for i = 9:-1:1
        push!(pq, i)
    end

    @test 1 == top(pq)
    @test 4 == height(pq)

    push!(pq, 9)
    @test 1 == top(pq)
    @test 4 == height(pq)
end

@testset "pop! 1" begin
    local pq = PriorityQueue{Int64}(min)

    push!(pq, 18)
    push!(pq, 19)
    push!(pq, 20)
    @test 18 == top(pq)

    pop!(pq)
    @test 19 == top(pq)

    pop!(pq)
    @test 20 == top(pq)

    pop!(pq)
    @test 0 == pq.size
end

@testset "pop! 2" begin
    local pq = PriorityQueue{Int64}(min)

    for i = 9:-1:1
        push!(pq, i)
    end

    @test 1 == top(pq)
    @test 4 == height(pq)

    pop!(pq)
    @test 2 == top(pq)
    @test 4 == height(pq)

    pop!(pq)
    @test 3 == top(pq)
    @test 3 == height(pq)

    pop!(pq)
    @test 4 == top(pq)
    @test 3 == height(pq)

    pop!(pq)
    @test 5 == top(pq)
    @test 3 == height(pq)
end

@testset "pop! 3" begin
    local N = 20
    local pq = PriorityQueue{Int64}(min)

    for i = (N+1):-1:1
        push!(pq, i)
    end

    @test 1 == top(pq)

    for i = 1:Int64(N / 2)
        pop!(pq)
    end
    @test (1 + N / 2) == top(pq)

    for i = 1:Int64((N / 2) - 1)
        pop!(pq)
    end

    @test N == top(pq)
end


text:
using Revise, ScatteringTransform, Plots, LinearAlgebra
pyplot()

k = 2^9
t = (-1+1/k):1/k:1; size(t) # interval
n = 20       # degree of polynomial
testSignal = hcat([t.^i for i=0:n]...)* randn(n+1)
plot(t, testSignal)

layers = stParallel(2, length(testSignal))
# 11 is *way* too far out
layers = stParallel(2, length(testSignal), CWTType = WT.Morlet(11), nScales=[8 for i=1:3])
layers.shears[1].averagingLength

# let's look at the actual filters that are used
n = ScatteringTransform.getNs(size(testSignal), layers)
i=1

layers.shears[i]
daughters,ξ = computeWavelets(n[i], layers.shears[i])
layers.shears[1].frameBound
ω = [0:ceil(Int, n[1]); -floor(Int,n[1])+1:-1]*2π;
plot(daughters[:, :], legend=false)
plot(daughters[:, end], legend=false)

plt = Array{Any,1}(undef,4)
for (k,x) in enumerate([(:none,"Unnormalized (may be using this one)"), (:sqrtScaling,
                                                 "2-norm Scaling"),
                        (:absScaling, "1-norm scaling"),
                        (:maxOne, "∞-norm scaling (they use)")])
    tmpDaughters = zeros(33,n[1]+1)
    for i=0:32
        tmpDaughters[i+1, :] = ScatteringTransform.Daughter(layers.shears[1],
                                                            2.0^(2+i/8), ω,
                                                            x[1])[1:(n[1]+1)]
    end
    plt[k] = plot(tmpDaughters[:,:]', legend=false, title="$(x[2])")
end
plot(plt...,layout=(4,1))
savefig("possibleScalings.pdf")
norm(ScatteringTransform.Daughter(layers.shears[i], 2.0^(3/8), ω)[1:(n[1]+1)])
layers.shears[i]

stResult = st(testSignal, layers, absType())


text:
@testset "Attachments" begin
  println(testname("Attachment tests"))
  print("  [  ] Upload attachment (blob mode) ")
  data = createdoc(db, Dict("item" => "screenshot"))
  result = put_attachment(db, data["id"], data["rev"], "test.png", "image/png", "../data/test.png")
  @test result["ok"] == true
  println("\r  [OK] Upload attachment (blob mode)")

  print("  [  ] Retrieve attachment (blob mode) ")
  att = get_attachment(db, result["id"], "test.png"; rev=result["rev"])
  open("../data/fetched.png", "w") do f
    write(f, att)
  end

  rm("../data/fetched.png")
  println("\r  [OK] Retrieve attachment (blob mode)")

  print("  [  ] Delete attachment (blob mode) ")
  result = delete_attachment(db, result["id"], result["rev"], "test.png")
  @test result["ok"] == true
  println("\r  [OK] Delete attachment (blob mode)")
end

max_stars_repo_name: xpqz/couchzilla

max_issues_repo_name: xpqz/couchzilla

max_forks_repo_name: xpqz/couchzilla