text:
type LaserBounds{P<:LaserTagPOMDP}
    ubp::QMDPPolicy{P, Int}
    lbp::RandomPolicy
    max_depth::Int
    lb_not_same::Float64

    LaserBounds{P}() where P = new()
end

#=
function LaserBounds(p::LaserTagPOMDP)
    sol = QMDPSolver(max_iterations=1000)
    pol = solve(sol, p)
    lb_not_same = -p.step_cost/(1-discount(p))
    return LaserBounds{typeof(p)}(pol, lb)
end
=#

function bounds(l::LaserBounds, p::LaserTagPOMDP, b::ScenarioBelief)
    ps = particles(b)
    if !all(isterminal(p, s) for s in ps) && !isnull(previous_obs(b))
        if get(previous_obs(b)) == LaserTag.D_SAME_LOC
            @assert all(!isterminal(p, s) && s.robot==s.opponent for s in ps)
        else
            @assert any(isterminal(p, s) || s.robot!=s.opponent for s in ps)
        end
    end
    steps = l.max_depth-b.depth
    lb = ARDESPOT.branching_sim(p, l.lbp, b, steps)/length(b.scenarios)
    return bounds(l, p, ps, lb)
end

function bounds(l::LaserBounds, p::LaserTagPOMDP, particles, provided_lb=nothing)
    if all(isterminal(p, s) for s in particles)
        lb = 0.0
    elseif all(!isterminal(p, s) && s.robot==s.opponent for s in particles)
        lb = l.ubp.pomdp.tag_reward
    elseif !isa(provided_lb, Void)
        lb = provided_lb
    elseif length(particles) == 1
        lb = state_value(l.ubp, first(particles))
    else
        lb = l.lb_not_same
    end

    vsum = 0.0
    for s in particles
        vsum += state_value(l.ubp, s)
    end
    ub = vsum/length(particles)
    # if lb > ub
    #     warn("lb > ub")
    #     @show particles
    # end
    return lb, vsum/length(particles)
end

function init_bounds(l::LaserBounds, p::LaserTagPOMDP, max_depth::Int, rng=Base.GLOBAL_RNG)
    sol = QMDPSolver(max_iterations=1000)
    l.ubp = solve(sol, p)
    l.lbp = solve(RandomSolver(rng), p)
    l.max_depth = max_depth
    l.lb_not_same = -p.step_cost/(1-discount(p))
    return l
end

init_bounds(l::LaserBounds, p::LaserTagPOMDP, s::DESPOTSolver) = init_bounds(l, p, s.D, s.rng)
# DESPOT.init_bounds(l::LaserBounds, p::LaserTagPOMDP, c::DESPOT.DESPOTConfig) = init_bounds(l, p, c.search_depth)
#
# function DESPOT.bounds{S}(l::LaserBounds, p::LaserTagPOMDP, b::Vector{DESPOT.DESPOTParticle{S}}, ::DESPOT.DESPOTConfig)
#     return bounds(l, p, p.state for p in b)
# end
#
# DESPOT.default_action(l::LaserBounds, pomdp::POMDP, particles, c) = default_action(particles)

function default_action(particles)
    if all(s.robot==s.opponent for s in particles)
        return LaserTag.TAG_ACTION
    else
        warn("non-tag default")
        return rand(1:4)
    end
end

struct NoGapTag end

function ARDESPOT.default_action(ngt::NoGapTag, b, ex)
    if isa(ex, NoGap)
        return default_action(particles(b))
    else
        rethrow(ex)
    end
    return -1
end

function nogap_tag(b, ex)
    if isa(ex, NoGap)
        @assert all(!s.terminal && s.robot==s.opponent for s in iterator(b))
        # @assert ng.value == l.ubp.pomdp.tag_reward
        return LaserTag.TAG_ACTION
    else
        rethrow(ex)
    end
end

immutable InevitableInit end

POMCPOW.init_N(m::InevitableInit, pomdp::LaserTagPOMDP, h::BeliefNode, a::Int) = h.node == 1 ? 0 : 1

function POMCPOW.init_V(m::InevitableInit, pomdp::LaserTagPOMDP, h::BeliefNode, a::Int)
    # only works for POMCPOW now because of node access
    if a == LaserTag.TAG_ACTION
        if h.node == 1
            return 0.0
        elseif h.tree.o_labels[h.node] == LaserTag.C_SAME_LOC
            return pomdp.tag_reward
        else
            return -pomdp.tag_reward - pomdp.step_cost
        end
    else
        return -pomdp.step_cost
    end
end

struct RootToNextMLFirst
    rng::MersenneTwister
end

function MCTS.next_action(gen::RootToNextMLFirst, p::VDPTagPOMDP, b, node)
    if isroot(node) && n_children(node) < 1
        target_sum=MVector(0.0, 0.0)
        agent_sum=MVector(0.0, 0.0)
        for s in particles(b::ParticleCollection)
            target_sum += s.target
            agent_sum += s.agent
        end
        next = VDPTag.next_ml_target(mdp(p), target_sum/n_particles(b))
        diff = next-agent_sum/n_particles(b)
        return TagAction(false, atan2(diff[2], diff[1]))
    else
        return rand(gen.rng, actions(p))
    end
end

MCTS.next_action(gen::RootToNextMLFirst, p::Union{GenerativeBeliefMDP,MeanRewardBeliefMDP}, b, node) = next_action(gen, p.pomdp, b, node)

# XXXD
# function MCTS.next_action(gen::RootToNextMLFirst, p::VDPTag.DiscreteVDPTagProblem, b, node)
#     if isroot(node) && n_children(node) < 1
#         cpc = ParticleCollection([convert_s(state_type(cproblem(p)), s, p) for s in b.particles])
#         ca = next_action(gen, cproblem(p), cpc, node)
#         return convert_a(action_type(p), ca, p)
#     else
#         return rand(gen.rng, actions(p))
#     end
# end

struct ModeMDP <: Policy
    vi::ValueIterationPolicy
end

struct ModeMDPSolver <: Solver
    vi::ValueIterationSolver
end

ModeMDPSolver(;kwargs...) = ModeMDPSolver(ValueIterationSolver(;kwargs...))

POMDPs.solve(sol::ModeMDPSolver, pomdp::POMDP) = ModeMDP(solve(sol.vi, pomdp))

POMDPs.action(p::ModeMDP, b) = action(p.vi, mode(b))

struct VDPUpper end

function ARDESPOT.ubound(ub::VDPUpper, pomdp::POMDP, b::ScenarioBelief)
    if all(isterminal(pomdp, s) for s in particles(b))
        return 0.0
    else
        return mdp(cproblem(pomdp)).tag_reward
    end
end

struct SampleRollout{P<:Policy, RNG<:AbstractRNG}
    policy::P
    rng::RNG
end

function MCTS.estimate_value(ro::SampleRollout, mdp::Union{GenerativeBeliefMDP, MeanRewardBeliefMDP}, b, steps::Int)
    start_state = rand(ro.rng, b)
    sim = POMDPToolbox.RolloutSimulator(ro.rng,
                                        Nullable{Any}(start_state),
                                        Nullable{Float64}(),
                                        Nullable{Int}(steps))
    return POMDPToolbox.simulate(sim, mdp.pomdp, ro.policy, start_state)
end


text:
module Buffers

# TODO: Static buffers!!!

import ..reset!

export add!

abstract type AbstractBuffer end

function add! end
# function reset! end

struct Transition{O,A}
    observation::O
    action::A
    reward::Float64
end

function MonteCarloExperiences end
function TDExperiences end

include("DynamicBuffer.jl")
export DynamicBuffer, DynamicBinaryBuffer

include("SingleStepBuffer.jl")
export SingleStepBuffer, ready

include("MultiStepDynamicBuffer.jl")
export MultiStepDynamicBuffer

include("FixedSizeBuffer.jl")
export FixedSizeBuffer, take_last, sample!

end

ext: jl


text:
"""
    zi = interpola(z,ns,nt)

Interpolate field `z` to size `(ns,ns,nt)` using nearest neighbors.
"""
	function interpola(z,ns,nt)

        nss=size(z);
        if(length(nss)>=3)
            nas=nss[1]; nat=nss[3];
        else
            nas=nss[1]; nat=1;
        end

	sdim=div(ns,nas);
	tdim=div(nt,nat);
	ir=1:sdim; jr=ir; kr=1:tdim;
	rs=0:sdim:(ns-sdim);
	rt=0:tdim:(nt-tdim);
	zi=zeros(ns,ns,nt);
	for i=1:sdim
	 for j=1:sdim
	  for k=1:tdim
		zi[i.+rs,j.+rs,k.+rt]=z;
	  end
	 end
	end

        return zi
        end


text:
module Algorithms

using DifferentialForms
using StaticArrays

export barycentre
barycentre(xs::SVector{N,<:Form{D,1,T}}) where {N,D,T} = sum(xs) / length(xs)
barycentre(xs::SVector{N,SVector{D,T}}) where {N,D,T} = sum(xs) / length(xs)

# function circumcentre1(xs::SVector{R, <:Chain{V, 1, T}}) where {R, V, T}
#     # G. Westendorp, A formula for the N-circumsphere of an N-simplex,
#     # <https://westy31.home.xs4all.nl/Circumsphere/ncircumsphere.htm>,
#     # April 2013.
#     @assert iseuclidean(V)
#     D = ndims(V)
#     @assert R == D + 1
#
#     # Convert Euclidean to conformal basis
#     cxs = conformal.(xs)
#     # Circumsphere (this formula is why we are using conformal GA)
#     X = ∧(cxs)
#     # Hodge dual
#     sX = ⋆X
#     # Euclidean part is centre
#     cc = euclidean(sX)
#
#     # Calculate radius
#     # TODO: Move this into a test case
#     r2 = scalar(abs2(cc)).v - 2 * sX.v[1]
#     # Check radii
#     for i ∈ 1:R
#         ri2 = scalar(abs2(xs[i] - cc)).v
#         @assert abs(ri2 - r2) ≤ T(1.0e-12) * r2
#     end
#
#     cc::Chain{V, 1, T}
# end

export circumcentre
function circumcentre(xs::SVector{N,<:Form{D,1,T}}) where {N,D,T}
    D::Int
    @assert D ≥ 0
    N::Int
    @assert N ≥ 1
    @assert N ≤ D + 1
    # See arXiv:1103.3076v2 [cs.RA], section 10.1
    A = SMatrix{N + 1,N + 1}(i ≤ N && j ≤ N ? 2 * (xs[i] ⋅ xs[j])[] :
                             i == j ? zero(T) : one(T)
                             for i in 1:(N + 1), j in 1:(N + 1))
    b = SVector{N + 1}(i ≤ N ? (xs[i] ⋅ xs[i])[] : one(T) for i in 1:(N + 1))
    c = A \ b
    cc = sum(c[i] * xs[i] for i in 1:N)
    return cc::Form{D,1,T}
end

function circumcentre(xs::SVector{N,<:Form{D,1,T}},
                      ws::SVector{N,<:Form{D,0,T′}}) where {N,D,T,T′}
    # See arXiv:math/0508188, section 3.2:
    # See arXiv:1103.3076v2 [cs.RA], section 10.1:
    #
    # C = Ci xi + Cj xj + Ck xk
    # (C - xi) (C - xi) = r^2 + wi
    # (C - xj) (C - xj) = r^2 + wj
    # (C - xk) (C - xk) = r^2 + wk
    # Ci + Cj + Ck = 1
    #
    # (Ci xi + Cj xj + Ck xk - xi) (Ci xi + Cj xj + Ck xk - xi) = r^2 + wi
    # (Ci xi + Cj xj + Ck xk - xj) (Ci xi + Cj xj + Ck xk - xj) = r^2 + wj
    # (Ci xi + Cj xj + Ck xk - xk) (Ci xi + Cj xj + Ck xk - xk) = r^2 + wk
    # Ci + Cj + Ck = 1
    #
    # C^2 - 2 xi (Ci xi + Cj xj + Ck xk) + xi^2 = r^2 + wi
    # C^2 - 2 xj (Ci xi + Cj xj + Ck xk) + xj^2 = r^2 + wj
    # C^2 - 2 xk (Ci xi + Cj xj + Ck xk) + xk^2 = r^2 + wk
    # Ci + Cj + Ck = 1
    #
    # 2 xi (Ci xi + Cj xj + Ck xk) + (r^2 - C^2) = xi^2 - wi
    # 2 xj (Ci xi + Cj xj + Ck xk) + (r^2 - C^2) = xj^2 - wj
    # 2 xk (Ci xi + Cj xj + Ck xk) + (r^2 - C^2) = xk^2 - wk
    # Ci + Cj + Ck = 1

    D::Int
    @assert D ≥ 0
    N::Int
    @assert N ≥ 1
    @assert N ≤ D + 1
    # See arXiv:1103.3076v2 [cs.RA], section 10.1
    A = SMatrix{N + 1,N + 1}(i ≤ N && j ≤ N ? 2 * (xs[i] ⋅ xs[j])[] :
                             i == j ? zero(T) : one(T)
                             for i in 1:(N + 1), j in 1:(N + 1))
    b = SVector{N + 1}(i ≤ N ? (xs[i] ⋅ xs[i] - ws[i])[] : one(T)
                       for i in 1:(N + 1))
    c = A \ b
    cc = sum(c[i] * xs[i] for i in 1:N)
    return cc::Form{D,1}
end

export volume
"""
Unsigned volume
"""
function volume(xs::SVector{N,<:Form{D,1,T}}) where {N,D,T}
    D::Int
    @assert D ≥ 0
    N::Int
    @assert N ≥ 1
    @assert N ≤ D + 1
    ys = map(x -> x - xs[1], deleteat(xs, 1))
    vol0 = norm(∧(ys))
    if T <: Rational
        vol = rationalize(typeof(zero(T).den), vol0; tol=sqrt(eps(vol0)))
    else
        vol = vol0
    end
    vol::T
    vol /= factorial(N - 1)
    return vol::T
end
function volume(xs::SVector{N,SVector{D,T}}) where {N,D,T}
    return volume(SVector{N}(Form{D,1,T}(x) for x in xs))
end

export signed_volume
"""
Signed volume
"""
function signed_volume(xs::SVector{N,<:Form{D,1,T}},
                       signature::Int=1) where {N,D,T}
    D::Int
    @assert D ≥ 0
    N::Int
    R = N - 1
    @assert 0 ≤ R ≤ D
    # TODO: add metric signatures to DifferentialForms.jl
    # TODO: add outermorphisms to DifferentialForms.jl
    # TODO: move Plotting into its own module; just provide functions
    #       to obtain vertices, connectivity, etc.
    ys = map(x -> x - xs[1], deleteat(xs, 1))
    η = SVector{D,T}(a == D ? signature : 1 for a in 1:D)
    ys′ = map(y -> Form{D,1,T}(η .* y.elts), ys)
    y = ∧(ys)
    y′ = ∧(ys′)
    v2 = (y′ ⋅ y)[]::T
    sv2 = sign(v2)::T
    v = (sv2 * sqrt(abs(v2)))::T
    return v / factorial(R)
end

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

# See also <https://github.com/JaneliaSciComp/Morton.jl>
export morton
function morton(xs::SVector{D,T}, xmin::SVector{D,T},
                xmax::SVector{D,T}) where {U<:Unsigned,D,T}
    return morton(UInt64, xs, xmin, xmax)
end
function morton(::Type{U}, xs::SVector{D,T}, xmin::SVector{D,T},
                xmax::SVector{D,T}) where {U<:Unsigned,D,T}
    D == 0 && return morton(SVector{D,U}())
    y = (xs - xmin) ./ (xmax - xmin)
    nbits = min(sizeof(U) - 1, 8 * sizeof(U) ÷ D)
    nbins = U(1) << nbits
    i = clamp.(floor.(U, nbins .* y), 0, nbins - 1)
    return morton(i)
end
function morton(i::SVector{D,U}) where {D,U<:Unsigned}
    D == 0 && return U(0)
    return reduce(|, SVector{D,U}(spread(Val(D), i[d]) << (d - 1) for d in 1:D))
end

spread(::Val{1}, i::Unsigned) = i
function spread(::Val{2}, i::UInt64)
    # can use 64 ÷ 2 = 32 bits
    @assert i ==
            i &
            0b0000000000000000000000000000000011111111111111111111111111111111
    i = i & 0b0000000000000000000000000000000011111111111111111111111111111111
    i = ((i << 16) | i) &
        0b0000000000000000111111111111111100000000000000001111111111111111
    i = ((i << 08) | i) &
        0b0000000011111111000000001111111100000000111111110000000011111111
    i = ((i << 04) | i) &
        0b0000111100001111000011110000111100001111000011110000111100001111
    i = ((i << 02) | i) &
        0b0011001100110011001100110011001100110011001100110011001100110011
    i = ((i << 01) | i) &
        0b0101010101010101010101010101010101010101010101010101010101010101
    return i
end
function spread(::Val{3}, i::UInt64)
    # can use 64 ÷ 3 = 21 bits
    @assert i ==
            i &
            0b000000000000000000000000000000000000000000111111111111111111111
    i = i & 0b000000000000000000000000000000000000000000111111111111111111111
    i = ((i << 32) | i) &
        0b000000000011111000000000000000000000000000000001111111111111111
    i = ((i << 16) | i) &
        0b000000000011111000000000000000011111111000000000000000011111111
    i = ((i << 08) | i) &
        0b001000000001111000000001111000000001111000000001111000000001111
    i = ((i << 04) | i) &
        0b001000011000011000011000011000011000011000011000011000011000011
    i = ((i << 02) | i) &
        0b001001001001001001001001001001001001001001001001001001001001001
    return i
end
function spread(::Val{4}, i::UInt64)
    # can use 64 ÷ 4 = 16 bits
    @assert i ==
            i &
            0b0000000000000000000000000000000000000000000000001111111111111111
    i = i & 0b0000000000000000000000000000000000000000000000001111111111111111
    i = ((i << 24) | i) &
        0b0000000000000000000000001111111100000000000000000000000011111111
    i = ((i << 12) | i) &
        0b0000000000001111000000000000111100000000000011110000000000001111
    i = ((i << 06) | i) &
        0b0000001100000011000000110000001100000011000000110000001100000011
    i = ((i << 03) | i) &
        0b0001000100010001000100010001000100010001000100010001000100010001
    return i
end
function spread(::Val{5}, i::UInt64)
    # can use 64 ÷ 5 = 12 bits
    @assert i ==
            i & 0b000000000000000000000000000000000000000000000000111111111111
    i = i & 0b000000000000000000000000000000000000000000000000111111111111
    i = ((i << 32) | i) &
        0b000000000000000011110000000000000000000000000000000011111111
    i = ((i << 16) | i) &
        0b000000000000000011110000000000000000111100000000000000001111
    i = ((i << 08) | i) &
        0b000000001100000000110000000011000000001100000000110000000011
    i = ((i << 04) | i) &
        0b000010000100001000010000100001000010000100001000010000100001
    return i
end
function spread(::Val{6}, i::UInt64)
    # can use 64 ÷ 6 = 10 bits
    @assert i ==
            i & 0b000000000000000000000000000000000000000000000000001111111111
    i = i & 0b000000000000000000000000000000000000000000000000001111111111
    i = ((i << 40) | i) &
        0b000000000011000000000000000000000000000000000000000011111111
    i = ((i << 20) | i) &
        0b000000000011000000000000000000001111000000000000000000001111
    i = ((i << 10) | i) &
        0b000000000011000000000011000000000011000000000011000000000011
    i = ((i << 05) | i) &
        0b000001000001000001000001000001000001000001000001000001000001
    return i
end

end


text:
using MathProgComplex, DataStructures, Memento, Printf, Pkg
MPC = MathProgComplex

using DataStructures, Memento, Test

@static if VERSION < v"0.7.0-DEV.2005"
    using Base.Test
else
    using Test
end

# suppress all messages less important than errors
setlevel!(getlogger(MathProgComplex), "error")

# @testset "MathProgComplex" begin
    include("basetypes.jl")
    include("cplx2real.jl")
    include("dat_importexport.jl")
    include("polynomial_types.jl")
    include("jump_export.jl")
# end


text:
import Base: convert, promote_rule, string, print, show
import Base.Dates: value, toms

# Convenience type for working with HH:MM:SS.
immutable Time <: TimePeriod
    seconds::Int
end
const ZERO = Time(0)

function Time(hour::Integer, minute::Integer, second::Integer)
    Time(hour * 3600 + minute * 60 + second)
end

function Time(s::AbstractString)
    # "-" represents 0:00 for some DST rules
    s == "-" && return ZERO
    parsed = map(n -> parse(Int, n), split(s, ':'))

    # Only can handle up to hour, minute, second.
    length(parsed) > 3 && throw(ArgumentError("Invalid Time string"))
    any(parsed[2:end] .< 0) && throw(ArgumentError("Invalid Time string"))

    # Handle variations where minutes and seconds may be excluded.
    values = [0,0,0]
    values[1:length(parsed)] = parsed

    if values[1] < 0
        for i in 2:length(values)
            values[i] = -values[i]
        end
    end

    return Time(values...)
end

# TimePeriod methods
value(t::Time) = t.seconds
toms(t::Time) = value(t) * 1000

hour(t::Time) = div(value(t), 3600)
minute(t::Time) = rem(div(value(t), 60), 60)
second(t::Time) = rem(value(t), 60)

function hourminutesecond(t::Time)
    h, r = divrem(value(t), 3600)
    m, s = divrem(r, 60)
    return h, m, s
end

convert(::Type{Second}, t::Time) = Second(value(t))
convert(::Type{Millisecond}, t::Time) = Millisecond(value(t) * 1000)
promote_rule{P<:Union{Week,Day,Hour,Minute,Second}}(::Type{P}, ::Type{Time}) = Second
promote_rule(::Type{Millisecond}, ::Type{Time}) = Millisecond

# Should be defined in Base.Dates
Base.isless(x::Period, y::Period) = isless(promote(x,y)...)

# https://en.wikipedia.org/wiki/ISO_8601#Times
function string(t::Time)
    neg = value(t) < 0 ? "-" : ""
    h, m, s = map(abs, hourminutesecond(t))
    @sprintf("%s%02d:%02d:%02d", neg, h, m, s)
end
print(io::IO, t::Time) = print(io, string(t))
show(io::IO, t::Time) = print(io, t)

# min/max offsets across all zones and all time.
# Note: A warning is given when we find an Olson zone that exceeds these values.
const MIN_GMT_OFFSET = Time("-15:56:00")  # Asia/Manilla
const MAX_GMT_OFFSET = Time("15:13:42")   # America/Metlakatla

# min/max save across all zones/rules and all time.
# Note: A warning is given when we find an Olson rule that exceeds these values.
const MIN_SAVE = Time("00:00")
const MAX_SAVE = Time("02:00")  # France, Germany, Port, Spain

const MIN_OFFSET = MIN_GMT_OFFSET + MIN_SAVE
const MAX_OFFSET = MAX_GMT_OFFSET + MAX_SAVE
const ABS_DIFF_OFFSET = abs(MAX_OFFSET - MIN_OFFSET)


text:
using Documenter, FASTX

# Build documentation.

makedocs(
    format = Documenter.HTML(),
    modules = [FASTX, FASTX.FASTQ, FASTX.FASTA],
    sitename = "FASTX.jl",
    doctest = false,
    strict = false,
    pages = [
        "Home" => "index.md",
        "Manual" => [
            "FASTA formatted files" => "manual/fasta.md",
            "FASTQ formatted files" => "manual/fastq.md"
        ],
        "Library" => [
            "Public" => "lib/public.md"
        ]
    ],
    authors = "Ben J. Ward, The BioJulia Organisation and other contributors."
)

deploydocs(
    repo = "github.com/BioJulia/FASTX.jl.git",
    deps = nothing,
    make = nothing
)


text:
a = readlines("mass_table.txt")
a = split.(a,"   ")
mass_dict = Dict()
for i in a
	mass_dict[i[1]] = parse(Float64,i[2])
end
protein = readlines("rosalind_prtm.txt")[1]
mass = sum(map(x->mass_dict[x],split(protein,"")))
round(mass,digits = 3)

max_stars_repo_name: Ved-Mahajan/JustForFun

max_issues_repo_name: Ved-Mahajan/JustForFun

max_forks_repo_name: Ved-Mahajan/JustForFun


text:
struct MBConv{E, D, X, P}
  expansion::E
  depthwise::D
  excitation::X
  projection::P

  do_expansion::Bool
  do_excitation::Bool
  do_skip::Bool
end
Flux.@functor MBConv

"""
    MBConv(
      in_channels, out_channels, kernel, stride;
      expansion_ratio, se_ratio)

Mobile Inverted Residual Bottleneck Block
([reference](https://arxiv.org/abs/1801.04381)).

# Arguments
- `in_channels`: Number of input channels.
- `out_channels`: Number of output channels.
- `expansion_ratio`:
  Expansion ratio defines the number of output channels.
  Set to `1` to disable expansion phase.
  `out_channels = input_channels * expansion_ratio`.
- `kernel`: Size of the kernel for the depthwise conv phase.
- `stride`: Size of the stride for the depthwise conv phase.
- `se_ratio`:
  Squeeze-Excitation ratio. Should be in `(0, 1]` range.
  Set to `-1` to disable.
"""
function MBConv(
  in_channels, out_channels, kernel, stride;
  expansion_ratio, se_ratio = 0.25,
)
  do_skip = stride == 1 && in_channels == out_channels
  do_expansion, do_excitation = expansion_ratio != 1, 0 < se_ratio ≤ 1
  pad, bias = SamePad(), false

  mid_channels = ceil(Int, in_channels * expansion_ratio)
  expansion = do_expansion ?
    Chain(
      Conv((1, 1), in_channels=>mid_channels; bias, pad),
      BatchNorm(mid_channels, swish)) :
    identity

  depthwise = Chain(
    Conv(kernel, mid_channels=>mid_channels; bias, stride, pad, groups=mid_channels),
    BatchNorm(mid_channels, swish))

  if do_excitation
    n_squeezed_channels = max(1, ceil(Int, in_channels * se_ratio))
    excitation = Chain(
      AdaptiveMeanPool((1, 1)),
      Conv((1, 1), mid_channels=>n_squeezed_channels, swish; pad),
      Conv((1, 1), n_squeezed_channels=>mid_channels; pad))
  else
    excitation = identity
  end

  projection = Chain(
    Conv((1, 1), mid_channels=>out_channels; pad, bias),
    BatchNorm(out_channels))
  MBConv(
    expansion, depthwise, excitation, projection, do_expansion,
    do_excitation, do_skip)
end

function (m::MBConv)(x)
  o = m.depthwise(m.expansion(x))

  if m.do_excitation
    o = σ.(m.excitation(o)) .* o
  end
  o = m.projection(o)
  if m.do_skip
    o = o + x
  end
  o
end


text:
using RNG.MersenneTwisters
using RNG.Random123
using RNG.PCG
using RNG.Xorshifts
import Base: rand

include("common.jl")

macro rngtest(rngs...)
    for rng in rngs
        if rng.head == :call
            r = @eval $rng
            name = rng.args[1]
        else
            r = @eval $(rng.args[1])
            name = rng.args[2]
        end
        @printf "%20s: %6.3f ns/64 bits\n" name speed_test(r)
    end
end

@rngtest(
Xorshift64(123),
Xorshift64Star(123),
Xorshift128Star(123),
Xorshift128Plus(123),
Xorshift1024Star(123),
Xorshift1024Plus(123),
Xoroshiro128Plus(123),
Xoroshiro128Star(123),

MersenneTwister(123),
MT19937(123),

(Threefry2x(UInt64, (123, 321)), Threefry2x64),
(Threefry4x(UInt64, (123, 321, 456, 654)), Threefry4x64),
(Threefry2x(UInt32, (123, 321)), Threefry2x32),
(Threefry4x(UInt32, (123, 321, 456, 654)), Threefry4x32),
(Philox2x(UInt64, 123), Philox2x64),
(Philox4x(UInt64, (123, 321)), Philox4x64),
(Philox2x(UInt32, 123), Philox2x32),
(Philox4x(UInt32, (123, 321)), Philox4x32),
(AESNI1x(123), AESNI1x128),
(ARS1x(123), ARS1x128),
(AESNI4x((123, 321, 456, 654)), AESNI4x32),
(ARS4x((123, 321, 456, 654)), ARS4x32),

(PCGStateOneseq(UInt64, PCG_XSH_RS, 123), PCG_XSH_RS_128),
(PCGStateOneseq(UInt64, PCG_XSH_RR, 123), PCG_XSH_RR_128),
(PCGStateOneseq(UInt64, PCG_RXS_M_XS, 123), PCG_RXS_M_XS_64),
(PCGStateOneseq(UInt64, PCG_XSL_RR, 123), PCG_XSL_RR_128),
(PCGStateOneseq(UInt64, PCG_XSL_RR_RR, 123), PCG_XSL_RR_RR_64),

(PCGStateOneseq(UInt32, PCG_XSH_RS, 123), PCG_XSH_RS_64),
(PCGStateOneseq(UInt32, PCG_XSH_RR, 123), PCG_XSH_RR_64),
(PCGStateOneseq(UInt32, PCG_RXS_M_XS, 123), PCG_RXS_M_XS_32),
(PCGStateOneseq(UInt32, PCG_XSL_RR, 123), PCG_XSL_RR_64),

(PCGStateOneseq(UInt128, PCG_RXS_M_XS, 123), PCG_RXS_M_XS_128),
(PCGStateOneseq(UInt128, PCG_XSL_RR_RR, 123), PCG_XSL_RR_RR_128),
)


text:
"""
    polyeval(z, p, N)

Evaluate the polynomial ``\\sum_k c[k] z^{k-1}`` for the coefficients `c[1]`, `c[2]`, ..., `c[N]`;
that is, the coefficients are given in ascending order by power of `z`.  This macro expands
to efficient inline code that uses either Horner's method.

```jldoctest

julia> p = randn(5)
5-element Array{Float64,1}:
  0.106455
  0.0351716
 -0.204764
 -0.87642
  0.260407

julia> polyeval(0.5, p, 5)
-0.020427265391245605

```
"""
function polyeval(x, p::Vector, N)
    return evalpoly(x,@view p[1:N])
end

function polyeval(x, p::Tuple, N)
    return evalpoly(x,p[1:N])
end
"""
    ConvergenceError([msg, val, niter, err])

Exception that is to be thrown when an iterative procedure fails
to converge.

"""
mutable struct ConvergenceError <: Exception
    msg::String
    val::Float64
    niter::Int
    err::Float64
    ConvergenceError() = new("Calculations failed to converge", NaN, -1, NaN)
    ConvergenceError(msg) = new(msg, NaN, -1, NaN)
    ConvergenceError(msg, val) = new(msg, val, -1, NaN)
    ConvergenceError(msg, val, niter) = new(msg, val, niter, NaN)
    ConvergenceError(msg, val, niter, err) = new(msg, val, niter, err)
end

"""
    calcz(vm, b0, c0, [EPS, [MAXITER, [relax]]])

Calculates the compressibility factor using a virial equation:

`z = 1 + b₀/z + c₀/vₘ²`

The algorithm transforms the values to a cubic polynomial and solves the roots using
cardan's method.

 * `b0` b₀ parameter of the virial equation
 * `c0` c₀ parameter of the virial equation
"""
function calcz(b0, c0)
    #transformation:
    #z3 = z2 + b₀z + c0
    #z3-z2- b₀z- c0
    poly =(-c0,-b0,-1.0,1.0)
    res =  cardan(poly)

    posibleres = filter(i->abs(imag(i)) < sqrt(eps(b0)),res)
    realres = map(real,posibleres)
    if length(realres)>=1
        return maximum(realres)
    else
    return error("not found")
    end
end

function cardan(poly::NTuple{4,T}) where {T<:AbstractFloat}
    # Cubic equation solver for complex polynomial (degree=3)
    # http://en.wikipedia.org/wiki/Cubic_function   Lagrange's method
    third = 0.3333333333333333 #speeds up the calculation
    a1  =  T(1.0) #in the original cardan, there is a divition here (1/a[4]). in this particular case its not necessary
    E1  = -complex(poly[3])*a1
    E2  =  complex(poly[2])*a1
    E3  = -complex(poly[1])*a1
    s0  =  E1
    E12 =  E1*E1
    A   =  2*E1*E12 - 9*E1*E2 + 27*E3 # = s1^3 + s2^3
    B   =  E12 - 3*E2                 # = s1 s2
    # quadratic equation: z^2 - Az + B^3=0  where roots are equal to s1^3 and s2^3
    Δ = sqrt(A*A - 4*B*B*B)
    if real(conj(A)*Δ)>=0 # scalar product to decide the sign yielding bigger magnitude
        s1 = exp(log(0.5 * (A + Δ)) * third)
    else
        s1 = exp(log(0.5 * (A - Δ)) * third)
    end
    if s1 == 0
        s2 = s1
    else
        s2 = B / s1
    end
    zeta1 = complex(-0.5, sqrt(T(3.0))*0.5)
    zeta2 = conj(zeta1)
    return third*(s0 + s1 + s2), third*(s0 + s1*zeta2 + s2*zeta1), third*(s0 + s1*zeta1 + s2*zeta2)
end


text:
# ## Coupled Van der Pol oscillator

# ## Model

module coupled_vanderpol #jl

using ReachabilityAnalysis

@taylorize function coupled_vanderpol!(du, u, p, t)
    du[1] = u[2]
    du[2] = (1.0 - u[1]^2)*u[2] - u[1] + (u[3]-u[1])
    du[3] = u[4]
    du[4] = (1.0 - u[3]^2)*u[4] - u[3] + (u[1]-u[3])
    return du
end

function model(X0)
    S = @system(x' = coupled_vanderpol!(x), dim:4)
    return IVP(S, X0)
end

end #jl

# ## References

max_stars_repo_name: mforets/ReachabilityLibrary

max_issues_repo_name: mforets/ReachabilityLibrary

max_forks_repo_name: mforets/ReachabilityLibrary


text:
import Entropies: TransferOperator, invariantmeasure, InvariantMeasure, Probabilities
using Entropies.GroupSlices
export TransferOperator

"""
	marginal_indices(x)

Returns a column vector `v` with the same number of elements as there are unique
elements in `x`. `v[i]` is the indices of elements in `x` matching `v[i]`.

For example, if the third unique element in `x`, and the element `u₃ = unique(x)[3]`
appears four times in `x`, then `v[3]` is a vector of four integers indicating the
position of the elements matching `u₃`.
"""
function marginal_indices(visited_bins, selected_axes)
    marginal_pts = [x[selected_axes] for x in visited_bins]
    groupinds(groupslices(marginal_pts))
end

"""
    marginal_probs_from_μ(seleced_axes, visited_bins, iv::InvariantMeasure, inds_μpositive)

Estimate marginal probabilities from a pre-computed invariant measure, given a set
of visited bins, an invariant measure and the indices of the positive-measure bins.
The indices in `selected_axes` determines which marginals are selected.
"""
function marginal_probs_from_μ(seleced_axes, visited_bins, iv::InvariantMeasure, inds_μpositive)

    marginal_inds::Vector{Vector{Int}} =
        marginal_indices(visited_bins, seleced_axes)

    # When the invariant measure over the joint space is already known, we don't
    # need to estimate histograms. We simply sum over the nonzero entries of the
    # (already estimated) invariant distribution `iv` in the marginal space
    # (whose indices are given by `seleced_axes`).
    μpos = iv.ρ[inds_μpositive]
    marginal = zeros(Float64, length(marginal_inds))
    @inbounds for i = 1:size(marginal_inds, 1)
        marginal[i] = sum(μpos[marginal_inds[i]])
    end
    return marginal
end

function transferentropy(s, t, est::TransferOperator{<:RectangularBinning}; base = 2, q = 1,
        τT = -1, τS = -1, η𝒯 = 1, dT = 1, dS = 1, d𝒯 = 1)

    emb = EmbeddingTE(τT = τT, τS = τS, η𝒯 = η𝒯, dT = dT, dS = dS, d𝒯 = d𝒯)

    joint_pts, vars, τs, js = te_embed(s, t, emb)
    iv = invariantmeasure(joint_pts, est.ϵ)

    # The bins visited by the orbit
    unique_visited_bins = unique(iv.to.bins)

    # # The subset of visited bins with nonzero measure
    inds_μpositive = findall(iv.ρ .> 0)
    positive_measure_bins = unique_visited_bins[inds_μpositive]

    # Estimate marginal probability distributions from joint measure
    cols_ST = [vars.S; vars.T]
    cols_T𝒯 = [vars.𝒯; vars.T]
    cols_T = vars.T
    p_T  = marginal_probs_from_μ(cols_T, positive_measure_bins, iv, inds_μpositive)
    p_ST = marginal_probs_from_μ(cols_ST, positive_measure_bins, iv, inds_μpositive)
    p_T𝒯 = marginal_probs_from_μ(cols_T𝒯, positive_measure_bins, iv, inds_μpositive)
    p_joint = iv.ρ[inds_μpositive]

    te = genentropy(Probabilities(p_ST), base = base, q = q) +
        genentropy(Probabilities(p_T𝒯), base = base, q = q) -
        genentropy(Probabilities(p_T), base = base, q = q) -
        genentropy(Probabilities(p_joint), base = base, q = q)
end

function transferentropy(s, t, c, est::TransferOperator{<:RectangularBinning}; base = 2, q = 1,
        τT = -1, τS = -1, τC = -1, η𝒯 = 1, dT = 1, dS = 1, dC = 1, d𝒯 = 1)

    emb = EmbeddingTE(τT = τT, τS = τS, τC = τC, η𝒯 = η𝒯, dT = dT, dS = dS, d𝒯 = d𝒯, dC = dC)

    joint_pts, vars, τs, js = te_embed(s, t, c, emb)
    iv = invariantmeasure(joint_pts, est.ϵ)

    # The bins visited by the orbit
    unique_visited_bins = unique(iv.to.bins)

    # # The subset of visited bins with nonzero measure
    inds_μpositive = findall(iv.ρ .> 0)
    positive_measure_bins = unique_visited_bins[inds_μpositive]

    # Estimate marginal probability distributions from joint measure
    cols_ST = [vars.S; vars.T; vars.C]
    cols_T𝒯 = [vars.𝒯; vars.T; vars.C]
    cols_T = [vars.T; vars.C]
    p_T  = marginal_probs_from_μ(cols_T, positive_measure_bins, iv, inds_μpositive)
    p_ST = marginal_probs_from_μ(cols_ST, positive_measure_bins, iv, inds_μpositive)
    p_T𝒯 = marginal_probs_from_μ(cols_T𝒯, positive_measure_bins, iv, inds_μpositive)
    p_joint = iv.ρ[inds_μpositive]

    te = genentropy(Probabilities(p_ST), base = base, q = q) +
        genentropy(Probabilities(p_T𝒯), base = base, q = q) -
        genentropy(Probabilities(p_T), base = base, q = q) -
        genentropy(Probabilities(p_joint), base = base, q = q)
end


text:
# functions related to chi-square distribution

import .RFunctions:
    chisqpdf,
    chisqlogpdf,
    chisqcdf,
    chisqccdf,
    chisqlogcdf,
    chisqlogccdf,
    chisqinvcdf,
    chisqinvccdf,
    chisqinvlogcdf,
    chisqinvlogccdf

# pdf for numbers with generic types
function chisqpdf(k::Real, x::Number)
  hk = k / 2  # half k
  1 / (2^(hk) * gamma(hk)) * x^(hk - 1) * exp(-x / 2)
end

# logpdf for numbers with generic types
function chisqlogpdf(k::Real, x::Number)
  hk = k / 2  # half k
  -hk * log(oftype(hk, 2)) - loggamma(hk) + (hk - 1) * log(x) - x / 2
end


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

# BAT-internal:
const MCMCOutputWithChains = Tuple{DensitySampleVector, MCMCBasicStats, AbstractVector{<:MCMCIterator}}

# BAT-internal:
function MCMCOutputWithChains(chainspec::MCMCSpec)
    dummy_chain = chainspec(zero(Int64))

    (
        DensitySampleVector(dummy_chain),
        MCMCBasicStats(dummy_chain),
        Vector{typeof(dummy_chain)}()
    )
end

# BAT-internal:
const MCMCOutput = Tuple{DensitySampleVector, MCMCBasicStats}

# BAT-internal:
function MCMCOutput(chainspec::MCMCSpec)
    samples, stats = MCMCOutputWithChains(chainspec::MCMCSpec)
    (samples, stats)
end

# BAT-internal:
function mcmc_sample(
    chainspec::MCMCSpec,
    nsamples::Integer,
    nchains::Integer;
    max_nsteps::Int64 = Int64(10 * nsamples),
    max_time::Float64 = Inf,
    tuner_config::AbstractMCMCTuningStrategy = AbstractMCMCTuningStrategy(chainspec.algorithm),
    convergence_test::MCMCConvergenceTest = BrooksGelmanConvergence(),
    init_strategy::MCMCInitStrategy = MCMCInitStrategy(tuner_config),
    burnin_strategy::MCMCBurninStrategy = MCMCBurninStrategy(chainspec.algorithm, nsamples, max_nsteps, tuner_config),
    granularity::Int = 1,
    strict_mode::Bool = false
)
    result = MCMCOutputWithChains(chainspec)

    result_samples, result_stats, result_chains = result

    (chains, tuners) = mcmc_init(
        chainspec,
        nchains,
        tuner_config,
        init_strategy
    )

    mcmc_tune_burnin!(
        (),
        tuners,
        chains,
        convergence_test,
        burnin_strategy;
        strict_mode = strict_mode
    )

    append!(result_chains, chains)

    mcmc_sample!(
        (result_samples, result_stats),
        result_chains,
        nsamples;
        max_nsteps = max_nsteps,
        max_time = max_time,
        granularity = granularity
    )

    result
end

# BAT-internal:
function mcmc_sample!(
    result::MCMCOutput,
    chains::AbstractVector{<:MCMCIterator},
    nsamples::Integer;
    max_nsteps::Int64 = Int64(100 * nsamples),
    max_time::Float64 = Inf,
    granularity::Int = 1
)
    result_samples, result_stats = result

    samples = DensitySampleVector.(chains)
    stats = MCMCBasicStats.(chains)

    nonzero_weights = granularity <= 1
    callbacks = [
        MCMCMultiCallback(
            MCMCAppendCallback(samples[i], nonzero_weights),
            MCMCAppendCallback(stats[i], nonzero_weights)
        ) for i in eachindex(chains)
    ]

    mcmc_iterate!(
        callbacks,
        chains;
        max_nsamples = Int64(nsamples),
        max_nsteps = max_nsteps,
        max_time = max_time
    )

    for x in samples
        merge!(result_samples, x)
    end

    for x in stats
        merge!(result_stats, x)
    end

    result
end

default_sampling_algorithm(posterior::AbstractPosteriorDensity) = MetropolisHastings()

"""
    function bat_sample(
        rng::AbstractRNG,
        posterior::AbstractPosteriorDensity,
        n::Union{Integer,Tuple{Integer,Integer}},
        algorithm::MCMCAlgorithm;
        max_nsteps::Integer,
        max_time::Real,
        tuning::AbstractMCMCTuningStrategy,
        init::MCMCInitStrategy,
        burnin::MCMCBurninStrategy,
        convergence::MCMCConvergenceTest,
        strict::Bool = false,
        filter::Bool = true
    )

Sample `posterior` via Markov chain Monte Carlo (MCMC).

`n` must be either a tuple `(nsteps, nchains)` or an integer. `nchains`
specifies the (approximate) number of MCMC steps per chain, `nchains` the
number of MCMC chains. If n is an integer, it is interpreted as
`nsteps * nchains`, and the number of steps and chains are chosen
automatically.
"""
function bat_sample(
    rng::AbstractRNG,
    posterior::AbstractPosteriorDensity,
    n::Tuple{Integer,Integer},
    algorithm::MCMCAlgorithm;
    max_nsteps::Integer = 10 * n[1],
    max_time::Real = Inf,
    tuning::AbstractMCMCTuningStrategy = AbstractMCMCTuningStrategy(algorithm),
    init::MCMCInitStrategy = MCMCInitStrategy(tuning),
    burnin::MCMCBurninStrategy = MCMCBurninStrategy(algorithm, n[1], max_nsteps, tuning),
    convergence::MCMCConvergenceTest = BrooksGelmanConvergence(),
    strict::Bool = false,
    filter::Bool = true
)
    #!!!! Temporary: ignore rng value, to be fixed:
    chainspec = MCMCSpec(algorithm, posterior)

    nsamples_per_chain, nchains = n

    samples, mcmc_stats, chains = mcmc_sample(
        chainspec,
        nsamples_per_chain,
        nchains;
        tuner_config = tuning,
        convergence_test = convergence,
        init_strategy = init,
        burnin_strategy = burnin,
        max_nsteps = Int64(max_nsteps),
        max_time = Float64(max_time),
        granularity = filter ? 1 : 2,
        strict_mode = strict
    )

    stats = bat_stats(mcmc_stats)

    (samples = samples, stats = stats, chains = chains)
end

function bat_sample(
    rng::AbstractRNG,
    posterior::AbstractPosteriorDensity,
    n::Integer,
    algorithm::MCMCAlgorithm;
    kwargs...
)
    nchains = 4
    nsamples = div(n, nchains)
    bat_sample(rng, posterior, (nsamples, nchains), algorithm; kwargs...)
end


text:
mutable struct RobustRestorer
    obj_val_R::Float64
    f_R::Vector{Float64}
    x_ref::Vector{Float64}

    theta_ref::Float64
    D_R::Vector{Float64}
    obj_val_R_trial::Float64

    pp::Vector{Float64}
    nn::Vector{Float64}
    zp::Vector{Float64}
    zn::Vector{Float64}

    dpp::Vector{Float64}
    dnn::Vector{Float64}
    dzp::Vector{Float64}
    dzn::Vector{Float64}

    pp_trial::Vector{Float64}
    nn_trial::Vector{Float64}

    inf_pr_R::Float64
    inf_du_R::Float64
    inf_compl_R::Float64

    mu_R::Float64
    tau_R::Float64
    zeta::Float64

    filter::Vector{Tuple{Float64,Float64}}
end

function RobustRestorer(ips::AbstractInteriorPointSolver)

    nn = Vector{Float64}(undef,ips.m)
    zp = Vector{Float64}(undef,ips.m)
    zn = Vector{Float64}(undef,ips.m)
    dpp= Vector{Float64}(undef,ips.m)
    dnn= Vector{Float64}(undef,ips.m)
    dzp= Vector{Float64}(undef,ips.m)
    dzn= Vector{Float64}(undef,ips.m)
    pp_trial = Vector{Float64}(undef,ips.m)
    nn_trial = Vector{Float64}(undef,ips.m)

    return RobustRestorer(
        0.,
        primal(ips._w2),
        primal(ips._w1),
        0.,
        primal(ips._w3),
        0.,
        dual(ips._w3),
        dual(ips._w4),
        zp, zn,
        dpp, dnn, dzp, dzn,
        dual(ips._w2),
        dual(ips._w1),
        0.,0.,0.,0.,0.,0.,
        Tuple{Float64,Float64}[],
    )
end

function initialize_robust_restorer!(ips::AbstractInteriorPointSolver)
    @trace(ips.logger,"Initializing restoration phase variables.")
    ips.RR == nothing && (ips.RR = RobustRestorer(ips))
    RR = ips.RR

    RR.x_ref .= ips.x
    RR.theta_ref = get_theta(ips.c)
    RR.D_R   .= min.(1,1 ./abs.(RR.x_ref))

    RR.mu_R = max(ips.mu,norm(ips.c,Inf))
    RR.tau_R= max(ips.opt.tau_min,1-RR.mu_R)
    RR.zeta = sqrt(RR.mu_R)

    RR.nn .= (RR.mu_R.-ips.opt.rho.*ips.c)./2 ./ips.opt.rho .+
        sqrt.(((RR.mu_R.-ips.opt.rho.*ips.c)./2 ./ips.opt.rho).^2 .+ RR.mu_R.*ips.c./2 ./ips.opt.rho)
    RR.pp .= ips.c .+ RR.nn
    RR.zp .= RR.mu_R./RR.pp
    RR.zn .= RR.mu_R./RR.nn

    RR.obj_val_R = get_obj_val_R(RR.pp,RR.nn,RR.D_R,ips.x,RR.x_ref,ips.opt.rho,RR.zeta)
    RR.f_R.=0
    empty!(RR.filter)
    push!(RR.filter,(ips.theta_max,-Inf))

    ips.l .= 0.
    ips.zl_r .= min.(ips.opt.rho,ips.zl_r)
    ips.zu_r .= min.(ips.opt.rho,ips.zu_r)
    ips.cnt.t = 0

    # misc
    ips.del_w = 0
end


text:
using SOtoZulip
using SOtoZulip: prune_questions
using Underscores
using LoggingFacilities, Logging

include("../configuration.jl")

logger = TimestampTransformerLogger(current_logger(), BeginningMessageLocation();
                                              format = "yyyy-mm-dd HH:MM:SS")
global_logger(logger)

global_zulip!(email = EMAIL, apikey = API_KEY, ep = ZULIP_EP)

const db = getdb(SODB)

qs = @_ get_questions(db)

for qss in Iterators.partition(qs, 50)
    prune_questions(db, qss)
end


text:
# Copyright (c) 2021 Idiap Research Institute, http://www.idiap.ch/
#  Niccolò Antonello <nantonel@idiap.ch>

using ArgParse

s = ArgParseSettings()
@add_arg_table! s begin
    "--conf"
        help = "configuration setup"
        arg_type = String
        default = "1a"
end
parsed_args = parse_args(ARGS, s)
setup = parsed_args["conf"]

using HMMGradients, Flux, Zygote
using Random, Statistics, LinearAlgebra
using FiniteStateTransducers
using DSP, MFCC
using BSON, JLD2, LibSndFile, FileIO, UUIDs

include("WFSTs.jl")
include("Models.jl")
include("Utils.jl")
include("conf/$(setup)/feat_conf.jl")
include("conf/$(setup)/model_conf.jl")

# get transition matrix
lexicon, ilexicon = get_lexicon()
H, L = get_HL(lexicon)
a, A = get_aA(H)
Ns = size(A,1)

# init model
modely = get_convnet(Nf,Ns;
                     Nks=Nks,
                     Nhs=Nhs,
                     strides=strides,
                     dilations=dilations,
                     dropout=dropout,
                     fout=fout)

# load training data
feat_dir = get_feat_dir(setup)
data = load(joinpath(feat_dir,"train.jld2"))
uttID2feats, uttID2tr, uttID2text, uttID2phones =
data["uttID2feats"], data["uttID2tr"], data["uttID2text"], data["uttID2phones"]
# curriculum data (isolated words)
uttID2text_cur = filter(x->length(split(x.second))==1,uttID2text)

# load dev data
data = load(joinpath(feat_dir,"dev.jld2"))
uttID2feats_dev, uttID2tr_dev, uttID2text_dev, uttID2phones_dev =
data["uttID2feats"], data["uttID2tr"], data["uttID2text"], data["uttID2phones"]
# curriculum data (isolated words)
uttID2text_dev_cur = filter(x->length(split(x.second))==1,uttID2text_dev)

# dataloaders
Xs_cur  = [uttID2feats[uttID]     for uttID in keys(uttID2text_cur) ]
Ys_cur  = [uttID2tr[uttID]        for uttID in keys(uttID2text_cur) ]
Xs      = [uttID2feats[uttID]     for uttID in keys(uttID2feats)    ]
Ys      = [uttID2tr[uttID]        for uttID in keys(uttID2tr)       ]

Xs_test_cur = [uttID2feats_dev[uttID] for uttID in keys(uttID2text_dev_cur)]
Ys_test_cur = [uttID2phones_dev[uttID]  for uttID in keys(uttID2text_dev_cur)]
Xs_test = [uttID2feats_dev[uttID] for uttID in keys(uttID2feats_dev)]
Ys_test = [uttID2phones_dev[uttID]  for uttID in keys(uttID2feats_dev)]

N_cur  = length(Xs_cur)
N      = length(Xs)
N_test = length(Xs_test)

cur_data   = Flux.Data.DataLoader((Xs_cur ,Ys_cur ), batchsize=Nb, shuffle=false)
train_data = Flux.Data.DataLoader((Xs     ,Ys     ), batchsize=Nb, shuffle=true)
test_data_cur  = Flux.Data.DataLoader((Xs_test_cur,Ys_test_cur), batchsize=Nb)
test_data  = Flux.Data.DataLoader((Xs_test,Ys_test), batchsize=Nb)

# define maximum likelihood function
function loss(Nt,t2tr,A,x,λ1)
  y  = modely(x)
  yp = exp.(y)
  f  = nlogMLlog(Nt,t2tr,A,y) +  λ1 * norm(yp,1)
  return f
end

function test(modely,a,A,ippsym,test_data)
  Flux.testmode!(modely)
  Nw = 0
  err = 0
  for (x,ps) in test_data
    x = feats_post.(x)
    Nts = ceil.(Int,size.(x,1)/3)
    xb = zeropad(x)
    y  = modely(xb)
    for i in eachindex(Nts)
      gamma = logposterior(Nts[i],a,A,view(y,:,:,i))
      ps_dec  = posterior2phones(ippsym,gamma)
      Nw  += length(ps[i])
      err += levenshtein(ps[i],ps_dec)
    end
  end
  per = err / Nw
  Flux.trainmode!(modely)
  return per
end

function train!(modely,a,A,H,opt,λ1,epochs,train_data,test_data)
  Flux.trainmode!(modely)
  N = length(train_data.data[1])
  ps = Flux.params(modely)
  best_per = Inf
  best_modely = deepcopy(modely)
  ippsym = get_iisym(H)
  for e in 1:epochs
    cost = 0
    for (x,t2trs) in train_data
      x = feats_post.(x)
      Nts = length.(t2trs) .+ 1
      xb = zeropad(x)
      train_loss, back =
      Zygote.pullback(() -> loss(Nts,t2trs,A,xb,λ1), ps)
      if isnan(train_loss) | isinf(train_loss)
        error("Nan/Inf cost function!!")
      end
      cost += train_loss
      gs = back(one(Float32))
      Flux.update!(opt, ps, gs)
    end
    per = test(modely,a,A,ippsym,test_data)
    save_best = per <= best_per
    if save_best
      best_modely = deepcopy(modely)
      best_per = per
      BSON.@save "models/$setup/current_modely.bson" best_modely
    end
    println("epoch: $e cost: $(round(cost/N,digits=4)) PER: $(round(per*100,digits=3))" * (save_best ? " ⋆ " : ""))
  end
  Flux.testmode!(best_modely)
  Flux.testmode!(modely)
  return best_modely, modely
end

model_folder = joinpath("models","$setup")
mkpath(model_folder)
println("Using setup: $setup")
println(read("conf/$setup/model_conf.jl",String))
if curriculum_training
  println("Curriculum training with $N_cur isolated words")
  best_modely, modely = train!(modely,a,A,H,opt,λ1,epochs_cur,cur_data,test_data_cur)
  BSON.@save joinpath(model_folder,"best_modely_curriculum.bson") best_modely
  modely = deepcopy(best_modely)
end
println("Training with $N utterances")
best_modely, modely = train!(modely,a,A,H,opt,λ1,epochs,train_data,test_data)
BSON.@save joinpath(model_folder,"best_modely_final.bson") best_modely
BSON.@save joinpath(model_folder,"modely.bson") modely


text:
"""
	landmark_isomap(coords, search=:knn, neigh=16, anchors=1500)

This function runs Landmark Isomap at 3-D points and returns a coordinate matrix
with the unfolded points.

## Parameters:

* `coords`  - coordinate matrix of the reference points.
* `search`  - search type to build neighbors graph for Isomap (:knn for
  k-nearest neighbor or :radius for radius search)
* `neigh`   - number of neighbors (for `search`=:knn) or radius
  distance (for `search`=:radius) to build neighbors graph for Isomap.
* `anchors` - number of anchors/landmark points for the dimensionality reduction.
"""
function landmark_isomap(coords::AbstractMatrix, search=:knn, neigh=16, anchors=1500)
	n, dim = size(coords, 2), 2
	n < anchors && (anchors = n)

	g, ianchors = graph_and_anchors(coords, search, neigh, anchors)
	ADM = Array{Float64}(undef, (anchors, anchors))
	dissmatrix!(ADM, g, ianchors)

	if n > anchors
		iothers = setdiff(1:n, ianchors)
		atcoords, M1, M3 = anchors_mds(ADM, dim)

		otcoords = Array{Float64}(undef, (3, length(iothers)))
		Threads.@threads for i in 1:length(iothers)
			otcoords[:,i] .= triangulation(g, iothers[i], ianchors, M1, M3)
		end
		tcoords = Array{Float64}(undef, (3,n))
		tcoords[:,ianchors] .= atcoords
		tcoords[:,iothers] .= otcoords
	else
		M = fit(MDS, ADM, maxoutdim=dim, distances=true)
		tcoords = vcat(transform(M),zeros(1, n))
	end
 	tcoords
end

function graph_and_anchors(ref_coords, nhood, neigh_val, nanchors)
	n = size(ref_coords,2)
	idxs, dists = get_neighbors(ref_coords, nhood, neigh_val, true)
	src, dest, dwgt = Int64[], Int64[], Float64[]

	for i in 1:n
		for j in 1:length(idxs[i])
			idxs[i][j] == i && continue
			push!(src, i)
			push!(dest, idxs[i][j])
			push!(dwgt, dists[i][j])
		end
	end

	g = SimpleWeightedGraph(Int.(src), Int.(dest), dwgt)
	comps = length(connected_components(g))
	@assert comps==1 "$comps subgroups of isolated vertices, need to increase number of neighbors"

	wgt = nhood==:radius ? [1/length(x) for x in idxs] : [mean(x) for x in dists]
	ianchors = n > nanchors ? sort!(sample(1:n, Weights(wgt), nanchors, replace=false)) : collect(1:n)
	g, ianchors
end

function dissmatrix!(ADM, g, iax::Vector{Int})
	n = size(ADM,1)
	for (i,ia) in enumerate(iax)
		dcols = dijkstra_shortest_paths(g, ia).dists[iax[i:n]]
		for (d,j) in zip(dcols,i:n)
			ADM[i,j] = ADM[j,i] = d
		end
	end
end

function anchors_mds(ADM, dims)
	nx = size(ADM,1)
	G = dmat2gram(ADM)
    F = eigen(Symmetric(G))
	λ = F.values
	sorti = sortperm(F.values,rev=true)
	sortλ = λ[sorti]
    EM = (F.vectors[:,sorti])[:,1:dims]
    sq_eigenvals = sortλ[1:dims].^0.5
    AM = Diagonal(sq_eigenvals)
    atcoords = vcat(permutedims(EM*AM),zeros(1,nx))

	# to use later for another points allocations
	M1 = permutedims(EM)
	M1 ./= reshape(sq_eigenvals,dims,1)
	M3 = Array{Float64}(undef,nx)
	mean!(M3,ADM.^2)

	atcoords, M1, M3
end

function triangulation(g::SimpleWeightedGraph, i, j, M1, M3)
	M2 = dijkstra_shortest_paths(g, i).dists[j] .^ 2
	vcat(-0.5*M1*(M2-M3),0)
end


text:
module TriggerWrapper

using JavaCall

const AbsoluteDate = @jimport org.orekit.time.AbsoluteDate
const AbstractDetector = @jimport org.orekit.propagation.events.AbstractDetector
const Action = @jimport org.hipparchus.ode.events.Action
const Class = @jimport java.lang.Class
const DateBasedManeuverTriggers = @jimport org.orekit.forces.maneuvers.trigger.DateBasedManeuverTriggers
const EventBasedManeuverTriggers = @jimport org.orekit.forces.maneuvers.trigger.EventBasedManeuverTriggers
const EventDetector = @jimport org.orekit.propagation.events.EventDetector
const EventHandler = @jimport org.orekit.propagation.events.handlers.EventHandler
const Field = @jimport org.hipparchus.Field
const FieldAbsoluteDate = @jimport org.orekit.time.FieldAbsoluteDate
const JString = @jimport java.lang.String
const ManeuverTriggers = @jimport org.orekit.forces.maneuvers.trigger.ManeuverTriggers
const Object = @jimport java.lang.Object
const ParameterDriver = @jimport org.orekit.utils.ParameterDriver
const RealFieldElement = @jimport org.hipparchus.RealFieldElement
const SpacecraftState = @jimport org.orekit.propagation.SpacecraftState
const Stream = @jimport java.util.stream.Stream

export DateBasedManeuverTriggers
export EventBasedManeuverTriggers
export ManeuverTriggers
export equals
export event_occurred
export get_class
export get_duration
export get_end_date
export get_events_detectors
export get_field_events_detectors
export get_name
export get_parameters_drivers
export get_start_date
export get_start_firing_detector
export get_stop_firing_detector
export get_triggered_end
export get_triggered_start
export hash_code
export init
export is_firing
export notify
export notify_all
export reset_state
export set_firing
export to_string
export wait

include("date_based_maneuver_triggers.jl")
include("event_based_maneuver_triggers.jl")
include("maneuver_triggers.jl")

end


text:
const term_c = (Tokens.RPAREN, Tokens.RSQUARE, Tokens.RBRACE, Tokens.END, Tokens.ELSE, Tokens.ELSEIF, Tokens.CATCH, Tokens.FINALLY, Tokens.ENDMARKER)

"""
Continue parsing statements until an element of `closers` is hit (usually
`end`). Statements are grouped in a `Block` EXPR.
"""
function parse_block(ps::ParseState, ret::Vector{EXPR}=EXPR[], closers=(Tokens.END,), docable=false; allow_const_field = false)
    prevpos = position(ps)
    while kindof(ps.nt) ∉ closers # loop until an expected closer is hit
        if kindof(ps.nt) ∈ term_c # error handling if an unexpected closer is hit
            if kindof(ps.nt) === Tokens.ENDMARKER
                break
            else
                push!(ret, mErrorToken(ps, EXPR(next(ps)), UnexpectedToken))
            end
        else
            if docable
                a = parse_doc(ps)
            else
                a = parse_expression(ps; allow_const_field = allow_const_field)
            end
            push!(ret, a)
        end
        prevpos = loop_check(ps, prevpos)
    end
    return ret
end

"""
Parses an iterator, allowing for the preceding keyword `outer`. Returns an
error expression if an invalid expression is parsed (anything other than
`=`, `in`, `∈`).
"""
function parse_iterator(ps::ParseState, outer = parse_outer(ps))
    arg = @closer ps :range @closer ps :ws @nocloser ps :wsop parse_expression(ps)
    if !is_range(arg)
        arg = mErrorToken(ps, arg, InvalidIterator)
    else
        arg = adjust_iter(arg)
    end
    if outer !== nothing
        arg.args[1] = setparent!(EXPR(:outer, EXPR[arg.args[1]], EXPR[outer]), arg)
        arg.fullspan += outer.fullspan
        arg.span = outer.fullspan + arg.span
    end
    return arg
end

function adjust_iter(x::EXPR)
    # Assumes x is a valid iterator
    if x.head === :call # isoperator(x.args[1]) && x.args[1].val in ("in", "∈")
        EXPR(EXPR(:OPERATOR, 0, 0, "="), EXPR[x.args[2], x.args[3]], EXPR[x.args[1]])
    elseif isassignment(x) && length(x.args) == 2 && x.args[1].head == :call && x.args[2].head == :block && length(x.args[2].args) == 1
        x.args[2] = setparent!(x.args[2].args[1], x)
        x
    else
        x
    end
end

"""
    is_range(x::EXPR)

Is `x` a valid iterator for use in `for` loops or generators?
"""
is_range(x::EXPR) = isassignment(x) || (x.head === :call && (is_in(x.args[1]) || is_elof(x.args[1])))

function parse_outer(ps)
    if kindof(ps.nt) === Tokens.OUTER && kindof(ps.nws) !== EmptyWS && !(Tokens.isoperator(kindof(ps.nnt)) && isbinaryop(ps.nnt) && kindof(ps.nnws) !== EmptyWS)
        EXPR(next(ps))
    end
end

"""
    parse_iterators(ps::ParseState, allowfilter = false)

Parses a group of iterators e.g. used in a `for` loop or generator. Can allow
for a succeeding `Filter` expression.
"""
function parse_iterators(ps::ParseState, iters, trivia)
    push!(iters, parse_iterator(ps))
    if iscomma(ps.nt)
        push!(trivia, accept_comma(ps))
        parse_iterators(ps, iters, trivia)
    end
end

"""
parse_filter(ps::ParseState, arg)

Parse a conditional filter following a generator.
"""
function parse_filter(ps::ParseState, arg)
    if kindof(ps.nt) === Tokens.IF # assumes we're inside a generator
        trivia = EXPR[EXPR(next(ps))]
        cond = @closer ps :range parse_expression(ps)
        if headof(arg) === :block
            arg = EXPR(:filter, EXPR[cond; arg.args], trivia)
        else
            arg = EXPR(:filter, EXPR[cond, arg], trivia)
        end
    end
    return arg
end

"""
    parse_call(ps, ret)

Parses a function call. Expects to start before the opening parentheses and is passed the expression declaring the function name, `ret`.
"""
function parse_call(ps::ParseState, ret::EXPR, ismacro=false)
    if is_minus(ret) || is_not(ret) || is_approx(ret)
        arg = @closer ps :unary @closer ps :inwhere @precedence ps PowerOp parse_expression(ps)
        if istuple(arg)
            pushfirst!(arg.args, ret)
            ret = EXPR(:call, arg.args, arg.trivia)
        else
            ret = EXPR(:call, EXPR[ret, arg], nothing)
        end
    elseif is_and(ret) || is_decl(ret) || is_exor(ret)
        arg = @precedence ps 20 parse_expression(ps)
        if is_exor(ret) && istuple(arg) && length(arg) == 3 && issplat(arg.args[1])
            arg = EXPR(:brackets, arg.args)
        end
        ret = EXPR(ret, EXPR[arg], nothing)
    else
        !ismacro && ismacroname(ret) && (ismacro = true)
        syntaxcall = is_issubt(ret) || is_issupt(ret)
        args = if syntaxcall
            EXPR[]
        elseif ismacro
            EXPR[ret, EXPR(:NOTHING, 0, 0)]
        else
            EXPR[ret]
        end
        # args = ismacro ? EXPR[ret, EXPR(:NOTHING, 0, 0)] : EXPR[ret]
        trivia = EXPR[EXPR(next(ps))]
        @closeparen ps @default ps parse_comma_sep(ps, args, trivia, !ismacro, insert_params_at = ismacro ? 3 : 2)
        accept_rparen(ps, trivia)
        ret = EXPR(ismacro ? :macrocall : syntaxcall ? ret : :call, args, trivia)
    end
    return ret
end

"""
Parses a comma separated list, optionally allowing for conversion of
assignment (`=`) expressions to `Kw`.
"""
function parse_comma_sep(ps::ParseState, args::Vector{EXPR}, trivia::Vector{EXPR}, kw = true, block = false, istuple = false; insert_params_at = 2)
    prevpos = position(ps)
    @nocloser ps :inwhere @nocloser ps :newline @closer ps :comma while !closer(ps)
        a = parse_expression(ps)
        if kw && _do_kw_convert(ps, a)
            a = _kw_convert(a)
        end
        push!(args, a)
        if iscomma(ps.nt)
            if istuple
                # As soon as we hit a comma in a list we can no longer convert to a block expr
                block = false
            end
            accept_comma(ps, trivia)
        elseif kindof(ps.ws) == SemiColonWS
            break
        elseif !closer(ps)
            # We've not hit a closing token, nor separating punctuation so let's insert an error-wrapped comma
            push!(trivia, EXPR(:errortoken, EXPR[EXPR(:COMMA, 0, 0)], nothing))
        end
        prevpos = loop_check(ps, prevpos)
    end
    if istuple && length(args) > 1
        block = false
    end

    if kindof(ps.ws) == SemiColonWS
        if @nocloser ps :inwhere @nocloser ps :newline @closer ps :comma @nocloser ps :semicolon closer(ps)
            if block && !(length(args) == 0 && ispunctuation(trivia[1])) && !(isunarycall(last(args)) && is_dddot(last(args).args[2]))
                push!(args, EXPR(:block, EXPR[pop!(args)]))
            else
                insert!(args, insert_params_at, EXPR(:parameters, EXPR[], nothing, 0, 0))
            end
        else
            a = @nocloser ps :newline @closer ps :comma @nocloser ps :inwhere parse_expression(ps)
            if block && !(length(args) == 0 && ispunctuation(trivia[1])) && !issplat(last(args)) && !(istuple && iscomma(ps.nt))
                args1 = EXPR[pop!(args), a]
                prevpos = position(ps)
                @nocloser ps :inwhere @nocloser ps :newline @closer ps :comma while @nocloser ps :semicolon !closer(ps)
                    a = parse_expression(ps)
                    push!(args1, a)
                    prevpos = loop_check(ps, prevpos)
                end
                body = EXPR(:block, args1)
                push!(args, body)
                args = body
            else
                parse_parameters(ps, args, EXPR[a], insert_params_at)
            end
        end
    end
    return
end

"""
    parse_parameters(ps::ParseState, args::Vector{EXPR}, args1::Vector{EXPR} = EXPR[]; usekw = true)

Parses parameter arguments for a function call (e.g. following a semicolon).
"""
function parse_parameters(ps::ParseState, args::Vector{EXPR}, args1::Vector{EXPR} = EXPR[], insert_params_at = 2; usekw = true)
    trivia = EXPR[]
    isfirst = isempty(args1)
    prevpos = position(ps)
    @nocloser ps :inwhere @nocloser ps :newline  @closer ps :comma while !isfirst || (@nocloser ps :semicolon !closer(ps))
        a = isfirst ? parse_expression(ps) : first(args1)
        if usekw && _do_kw_convert(ps, a)
            a = _kw_convert(a)
        end
        if isfirst
            push!(args1, a)
        else
            pop!(args1)
            push!(args1, a)
        end
        if iscomma(ps.nt)
            accept_comma(ps, trivia)
        elseif kindof(ps.ws) !== SemiColonWS && !closer(ps)
            push!(trivia, EXPR(:errortoken, EXPR[EXPR(:COMMA, 0, 0)], nothing))
        end
        if kindof(ps.ws) == SemiColonWS
            parse_parameters(ps, args1, EXPR[], 1; usekw=usekw)
        end
        prevpos = isfirst ? loop_check(ps, prevpos) : position(ps)
        isfirst = true
    end
    if !isempty(args1)
        insert!(args, insert_params_at, EXPR(:parameters, args1, trivia))
    end
    return
end

function parse_macroname(ps)
    at = EXPR(ps)
    if isidentifier(ps.nt) || isoperator(ps.nt) || iskeyword(ps.nt)
        if !isemptyws(ps.ws)
            ws = ps.ws.endbyte - ps.ws.startbyte + 1
            mname = INSTANCE(next(ps))
            mname.val = valof(mname) isa String ? string("@", " "^ws, valof(mname)) : string("@", " "^ws)
            mname.span += ws
            mname.fullspan += ws
            mname = mErrorToken(ps, mname, UnexpectedWhiteSpace)
        else
            next(ps)
            if VERSION > v"1.3.0-" && val(ps.t, ps) == "var" && (kindof(ps.nt) === Tokens.STRING || kindof(ps.nt) === Tokens.TRIPLE_STRING)
                mname = EXPR(:IDENTIFIER, max(1, ps.nt.startbyte - ps.t.startbyte + 1), max(1, ps.t.endbyte - ps.t.startbyte + 2), string("@", val(ps.t, ps)))
                arg = parse_string_or_cmd(next(ps), mname)

                EXPR(:NONSTDIDENTIFIER, EXPR[mname, arg], nothing)

            else
                # set span/fullspan min length at 1 to account for the case of a lonely '@'
                EXPR(:IDENTIFIER, max(1, ps.nt.startbyte - ps.t.startbyte + 1), max(1, ps.t.endbyte - ps.t.startbyte + 2), string("@", val(ps.t, ps)))
            end
        end
    else
        mErrorToken(ps, at, MalformedMacroName)
    end
end

"""
    parse_macrocall(ps)

Parses a macro call. Expects to start on the `@`.
"""
function parse_macrocall(ps::ParseState)
    mname = parse_macroname(ps)
    # Handle cases with @ at start of dotted expressions
    if kindof(ps.nt) === Tokens.DOT && isemptyws(ps.ws)
        prevpos = position(ps)
        while kindof(ps.nt) === Tokens.DOT
            op = EXPR(:OPERATOR, next(ps))
            nextarg = EXPR(:IDENTIFIER, next(ps))
            mname = EXPR(op, EXPR[mname, EXPR(:quotenode, EXPR[nextarg], nothing)], nothing)
            prevpos = loop_check(ps, prevpos)
        end
    end

    if iscomma(ps.nt) || isnewlinews(ps.ws)
        return EXPR(:macrocall, EXPR[mname, EXPR(:NOTHING, 0, 0)], nothing, mname.fullspan, mname.span)
    elseif isemptyws(ps.ws) && kindof(ps.nt) === Tokens.LPAREN
        return parse_call(ps, mname, true)
    elseif isemptyws(ps.ws) && kindof(ps.nt) === Tokens.LBRACE
        next(ps)
        return EXPR(:macrocall, EXPR[mname, EXPR(:NOTHING, 0, 0), @default ps @closebrace ps parse_braces(ps)], nothing)
    elseif isemptyws(ps.ws) && kindof(ps.nt) === Tokens.LSQUARE
        next(ps)
        return EXPR(:macrocall, EXPR[mname, EXPR(:NOTHING, 0, 0), @default ps parse_array(ps)], nothing)
    else
        args = EXPR[mname, EXPR(:NOTHING, 0, 0)]
        insquare = ps.closer.insquare
        prevpos = position(ps)
        @default ps while !closer(ps)
            if insquare
                a = @closer ps :insquare @closer ps :inmacro @closer ps :ws @closer ps :wsop parse_expression(ps)
            else
                a = @closer ps :inmacro @closer ps :ws @closer ps :wsop parse_expression(ps)
            end
            push!(args, a)

            if valof(mname) == "@doc" && ps.t.endpos[1] + 1 == ps.nt.startpos[1] && length(args) == 3 && kindof(ps.nt) !== Tokens.ENDMARKER
                a = parse_expression(ps)
                push!(args, a)
                break
            end
            if ps.closer.for_generator && kindof(ps.nt) === Tokens.FOR
                break
            end
            prevpos = loop_check(ps, prevpos)
        end
        return EXPR(:macrocall, args, nothing)
    end
end

"""
parse_generator(ps)

Having hit `for` not at the beginning of an expression return a generator.
Comprehensions are parsed as SQUAREs containing a generator.
"""
function parse_generator(ps::ParseState, first::EXPR)
    kw = EXPR(next(ps))
    iters, trivia = EXPR[], EXPR[]
    @closesquare ps parse_iterators(ps, iters, trivia)
    if kindof(ps.nt) === Tokens.IF # filter
        push!(trivia, EXPR(next(ps)))
        cond = @closer ps :range parse_expression(ps)
        pushfirst!(iters, cond)
        iters = EXPR(:filter, iters, trivia)
        trivia = EXPR[]
    end
    if kindof(ps.nt) === Tokens.FOR
        return EXPR(:flatten, EXPR[EXPR(:generator, EXPR[parse_generator(ps, first); iters], [kw; trivia])], nothing)
    else
        return EXPR(:generator, EXPR[first; iters], EXPR[kw; trivia])
    end
end

function get_appropriate_child_to_expand(x)
    if headof(x) === :generator && !(headof(x.args[1]) in (:generator, :flatten))
        return x, x.args[1]
    elseif headof(x) === :flatten &&  headof(x.args[1]) === :generator && headof(x.args[1].args[1]) === :generator
        x.args[1], x.args[1].args[1]
    else
        get_appropriate_child_to_expand(x.args[1])
    end
end

is_nonstd_identifier(ps) = VERSION > v"1.3.0-" && isidentifier(ps.nt) && isemptyws(ps.nws) && (kindof(ps.nnt) === Tokens.STRING || kindof(ps.nnt) === Tokens.TRIPLE_STRING)

function parse_nonstd_identifier(ps)
    id = INSTANCE(next(ps))
    if valof(id) == "var"
        EXPR(:NONSTDIDENTIFIER, EXPR[id, INSTANCE(next(ps))])
    else
        mErrorToken(ps, EXPR(:NONSTDIDENTIFIER, EXPR[id, INSTANCE(next(ps))]), UnexpectedToken)
    end
end

function parse_importexport_item(ps, is_colon = false)
    if kindof(ps.nt) === Tokens.AT_SIGN
        parse_macroname(next(ps))
    elseif kindof(ps.nt) === Tokens.LPAREN
        a = EXPR(:brackets, EXPR[], EXPR[EXPR(next(ps))])
        push!(a, @closeparen ps parse_expression(ps))
        pushtotrivia!(a, accept_rparen(ps))
        a
    elseif kindof(ps.nt) === Tokens.EX_OR
        parse_unary(ps, INSTANCE(next(ps)))
    elseif !is_colon && isoperator(ps.nt)
        next(ps)
        EXPR(:OPERATOR, ps.nt.startbyte - ps.t.startbyte,  1 + ps.t.endbyte - ps.t.startbyte, val(ps.t, ps))
    elseif is_nonstd_identifier(ps)
        parse_nonstd_identifier(ps)
    else
        INSTANCE(next(ps))
    end
end
"""
Helper function for parsing import/using statements.
"""
function parse_dot_mod(ps::ParseState, is_colon=false, allow_as=false)
    ret = EXPR(EXPR(:OPERATOR, 0, 0, "."), EXPR[], EXPR[])

    prevpos = position(ps)
    while kindof(ps.nt) === Tokens.DOT || kindof(ps.nt) === Tokens.DDOT || kindof(ps.nt) === Tokens.DDDOT
        d = EXPR(:OPERATOR, next(ps))
        trailing_ws = d.fullspan - d.span
        if is_dot(d)
            push!(ret, EXPR(:OPERATOR, 1 + trailing_ws, 1, "."))
        elseif is_ddot(d)
            push!(ret, EXPR(:OPERATOR, 1, 1, "."))
            push!(ret, EXPR(:OPERATOR, 1 + trailing_ws, 1, "."))
        elseif is_dddot(d)
            push!(ret, EXPR(:OPERATOR, 1, 1, "."))
            push!(ret, EXPR(:OPERATOR, 1, 1, "."))
            push!(ret, EXPR(:OPERATOR, 1 + trailing_ws, 1, "."))
        end
        prevpos = loop_check(ps, prevpos)
    end

    prevpos = position(ps)
    while true
        push!(ret, parse_importexport_item(ps, is_colon))

        if kindof(ps.nt) === Tokens.DOT
            pushtotrivia!(ret, EXPR(next(ps)))
            if kindof(ps.nt) === Tokens.COLON
                u = parse_unary_colon(ps, EXPR(:IDENTIFIER, next(ps)))
                push!(ret, u)
                return ret
            end
        elseif isoperator(ps.nt) && (ps.nt.dotop || kindof(ps.nt) === Tokens.DOT)
            pushtotrivia!(ret, EXPR(:DOT, 1, 1))
            ps.nt = RawToken(kindof(ps.nt), ps.nt.startpos, ps.nt.endpos, ps.nt.startbyte + 1, ps.nt.endbyte, ps.nt.token_error, false, ps.nt.suffix)
        else
            @static if VERSION > v"1.6-"
                if allow_as && !isnewlinews(ps.ws) && isidentifier(ps.nt) && val(ps.nt, ps) == "as"
                    as = EXPR(:AS, next(ps))
                    as_val = parse_importexport_item(ps, is_colon)
                    ret = EXPR(:as, EXPR[ret, as_val], EXPR[as])
                end
            end
            break
        end
        prevpos = loop_check(ps, prevpos)
    end
    ret
end


text:
using Minibatch
using Revtok
using DataStructures

struct Vocab
    itos::Vector
    stoi::Associative
end
Base.length(vocab::Vocab) = length(vocab.itos)

function process(corpus, vocabsize=200, batchsize=4)
    corpus = tokenize.(split(corpus, '\n'))
    segments = buildvocab(counter(Iterators.flatten(corpus)), vocabsize)
    corpus = segment.(corpus, segments)
    itos = keys(segments)
    stoi = Dict((reverse(p) for p in enumerate(itos)))
    vocab = Vocab(collect(itos), stoi)
    corpus = [[vocab.stoi[x] for x in ex] for ex in corpus]
    corpus = Iterators.partition(corpus, batchsize)
    corpus = map(b -> MaskedBatch(b, (true,)), corpus)
    return corpus, vocab
end

english = """This is an account of how magical thinking made us modern.
When people talk about magical thinking, it is usually as a cognitive feature of children, uneducated people, the mushy-minded, or the mentally ill.
If we notice magical thinking in ourselves, it is with a pang of shame: literate adults are supposed to be more sophisticated than that.
At the same time, magical thinking is obviously rampant in the world.
It’s hard not to be fascinated, even if it’s a horrified fascination.
I think that the reason that it has been so difficult to precisely define “magical thinking” is that what we call “magical thinking” is a collection of stigmatized examples of a more general, and generally useful, cognitive capacity.
This is the ability to think in “as if” mode: “as if” inanimate objects had minds, “as if” thoughts could affect reality, “as if” symbols had power over their referents."""

spanish = """Este es un relato de cómo el pensamiento mágico nos hizo modernos.
Cuando la gente habla de pensamiento mágico, generalmente es como un rasgo cognitivo de los niños, de las personas sin educación, de los que tienen la mente blanda o de los enfermos mentales.
Si nos damos cuenta del pensamiento mágico en nosotros mismos, es con una espiral de vergüenza: se supone que los adultos alfabetizados son más sofisticados que eso.
Al mismo tiempo, el pensamiento mágico es obviamente desenfrenado en el mundo.
Es difícil no estar fascinado, aunque sea una fascinación horrorizada.
Creo que la razón por la que ha sido tan difícil definir con precisión el “pensamiento mágico” es que lo que llamamos "pensamiento mágico" es una colección de ejemplos estigmatizados de una capacidad cognitiva más general, y en general útil.
Esta es la capacidad de pensar en modo “como si”: “como si” los objetos inanimados tuvieran mente, “como si” los pensamientos pudieran afectar a la realidad, “como si” los símbolos tuvieran poder sobre sus referentes."""

#corpus = zip(english, spanish)
en, vocab_en = process(english)
es, vocab_es = process(spanish)


text:
@testset "Solve LP" begin
    cp = test_simpleLP()
    optimizer = GLPK.Optimizer
    lp = optimize_model(cp, optimizer)
    @test termination_status(lp) === MOI.OPTIMAL
    sol = JuMP.value.(lp[:x])
    @test sol ≈ [1.0, 2.0]

    optimizer = Clp.Optimizer
    lp = optimize_model(cp, optimizer)
    @test termination_status(lp) === MOI.OPTIMAL
    sol = JuMP.value.(lp[:x])
    @test sol ≈ [1.0, 2.0]
end


text:
function part1(input)
    containers = sort(parse.(Int, eachline(input)), rev = true)
    return combinations(containers, 150)
end

function part2(input)
    containers = sort(parse.(Int, eachline(input)), rev = true)
    return last(combinations2(containers, 150))
end

function combinations(containers, n, i = 1,
                      cache = Dict{Tuple{Int, Int}, Int }())
    n < 0 && return 0
    n == 0 && return 1
    i > length(containers) && return 0
    if !haskey(cache, (n, i))
        c = sum(combinations(containers, n - containers[j], j + 1, cache)
                for j = i:length(containers))
        cache[(n, i)] = c
    end
    return cache[(n, i)]
end

function combinations2(containers, n, i = 1,
                       cache = Dict{Tuple{Int, Int}, Tuple{Int, Int}}())
    n < 0 && return (0, 0)
    n == 0 && return (0, 1)
    i > length(containers) && return (0, 0)
    if !haskey(cache, (n, i))
        min_cups, num_combinations = typemax(Int), 0
        for j = i:length(containers)
            mc, nc = combinations2(containers, n - containers[j], j + 1, cache)
            nc == 0 && continue
            if mc < min_cups
                min_cups = mc
                num_combinations = nc
            elseif mc == min_cups
                num_combinations += nc
            end
        end
        if num_combinations == 0
            cache[(n, i)] = (0, 0)
        else
            cache[(n, i)] = (min_cups + 1, num_combinations)
        end
    end
    return cache[(n, i)]
end


text:
# Set up Problem
model = Dynamics.car
n,m = model.n, model.m
N = 51
Q = Diagonal(I,n)*0.01
R = Diagonal(I,m)*0.01
Qf = Diagonal(I,n)*0.01
xf = [0,1,0]
obj = LQRObjective(Q,R,Qf,xf,N)
n,m = model.n, model.m
bnd = BoundConstraint(n,m, x_min=[-0.5, -0.01, -Inf], x_max=[0.5, 1.01, Inf], u_min=[0.1,-2], u_max=1.2)
bnd1 = BoundConstraint(n,m, u_min=bnd.u_min, u_max=bnd.u_max)
bnd_x = BoundConstraint(n,m, x_min=[-0.5, -0.01, -Inf], x_max=[0.5, 1.01, Inf])
goal = goal_constraint(xf)
obs = (([0.2, 0.6], 0.25),
       ([-0.5, 0.5], 0.4),
       ([-0.15,0.85], 0.1))
obs1 = TO.planar_obstacle_constraint(n,m, obs[1]..., :obstacle1)
obs2 = TO.planar_obstacle_constraint(n,m, obs[2]..., :obstacle2)
obs3 = TO.planar_obstacle_constraint(n,m, obs[3]..., :obstacle3)
con = Constraints(N)
con[1] += bnd1
for k = 2:N-1
    con[k] += bnd  + obs1 + obs2 + obs3
end
con[N] += goal
prob = Problem(rk4(model), obj, constraints=con, tf=3)

# Solve with ALTRO
initial_controls!(prob, ones(m,N-1))
ilqr = iLQRSolverOptions()
al = AugmentedLagrangianSolverOptions(opts_uncon=ilqr)
al.penalty_initial = 0.001
al.constraint_tolerance = 1e-2
al.constraint_tolerance_intermediate = 1e-2
al.verbose = true
res, = solve(prob, al)
max_violation(res)

plot()
TO.plot_circle!(obs[1]...)
TO.plot_circle!(obs[2]...)
TO.plot_circle!(obs[3]...)
plot_trajectory!(res.X, markershape=:circle)
plot(res.U)

# Projected Newton
opts = ProjectedNewtonSolverOptions{Float64}(verbose=true)
solver = ProjectedNewtonSolver(res,opts)
TO.update!(res, solver)
V0 = copy(solver.V.V)

Z = TO.primals(solver.V)
a = solver.a.duals
Y,y = TO.active_constraints(prob, solver)
viol0 = norm(y,Inf)
HinvY = Diagonal(solver.H)\Y'
S = cholesky(Symmetric(Y*HinvY))

solver.opts.feasibility_tolerance = 1e-12
solver.opts.verbose = true
solver.opts.active_set_tolerance = 1e-6
begin
    copyto!(solver.V.V, V0)
    TO.update!(res, solver)
    TO.projection_solve!(prob, solver)
end
copyto!(solver.V.V, V0)
TO.update!(res, solver)
TO._projection_solve!(prob, solver, solver.V, true)
TO._projection_linesearch!(prob, solver, solver.V, S, Hinv*Y')

dZ = -Hinv*Y'*(S\y)
Z .+= dZ

TO.dynamics_constraints!(prob, solver)
TO.update_constraints!(prob, solver)
y = solver.y[a]
viol = norm(y,Inf)

TO.projection!(res, solver)

V_ = TO.newton_step!(res, solver)

plot!(solver.V.U)

A = rand(10,10)
B = A + 1e-4*I
b = rand(10)

x = reg_solve(A,b,1e-2)
norm(A*x-b)

max_issues_repo_name: MehrdadZareian/GMARA

max_forks_repo_name: MehrdadZareian/GMARA


text:
import RvSpectMLBase
using Test

@testset "Theoretical Instrument" begin
    using RvSpectMLBase.TheoreticalInstrument
    num_pixels = 100
    num_orders = 10
    #resolution = 1e5

    @test_nowarn TheoreticalInstrument1D()
    n1 = TheoreticalInstrument1D()

    @testset "1D Extracted Spectra Traits" begin
            @test typeof(n1) <: AbstractInstrument
            @test typeof(n1) <: AnyTheoreticalInstrument
            @test min_pixel(n1) == 1
            @test max_pixel(n1) == 128*8192
            @test first(pixels_all(n1)) == min_pixel(n1)
            @test last(pixels_all(n1)) == max_pixel(n1)
            @test default_ccf_mask_v_width(n1) > 0
            @test get_inst_module(n1) == TheoreticalInstrument
    end

    lambda2d = collect(reshape(range(5000.0,stop=6000.0,length=num_pixels*num_orders),num_pixels,num_orders))
    @test_nowarn TheoreticalInstrument2D(lambda2d)
    n2 = TheoreticalInstrument2D(lambda2d)

    @testset "2D Extracted Spectra Traits" begin
            @test typeof(n2) <: AbstractInstrument
            @test typeof(n2) <: AnyTheoreticalInstrument
            @test min_order(n2) == 1
            @test max_order(n2) == num_orders
            @test first(orders_all(n2)) == min_order(n2)
            @test last(orders_all(n2)) == max_order(n2)
            @test orders_to_use_default(n2) == min_order(n2):max_order(n2)
            @test min_pixel_in_order(n2) == 1
            @test max_pixel_in_order(n2) == num_pixels
            @test first(pixels_all(n2)) == min_pixel_in_order(n2)
            @test last(pixels_all(n2)) == max_pixel_in_order(n2)*(max_order(n2)-min_order(n2)+1)
            @test default_ccf_mask_v_width(n2) > 0
            @test get_inst_module(n2) == TheoreticalInstrument
    end
end

max_stars_repo_name: RvSpectML/RvSpectMLTemplate

max_issues_repo_name: RvSpectML/RvSpectMLTemplate

max_forks_repo_name: RvSpectML/RvSpectMLTemplate


text:
# adapted from Nemo
function (>){P}(x::ArbFloat{P}, y::ArbFloat{P})
    return Bool(ccall(@libarb(arb_gt), Cint, (Ptr{ArbFloat{P}}, Ptr{ArbFloat{P}}), &x, &y))
end
function (>=){P}(x::ArbFloat{P}, y::ArbFloat{P})
    return Bool(ccall(@libarb(arb_ge), Cint, (Ptr{ArbFloat{P}}, Ptr{ArbFloat{P}}), &x, &y))
end
function (<){P}(x::ArbFloat{P}, y::ArbFloat{P})
    return Bool(ccall(@libarb(arb_lt), Cint, (Ptr{ArbFloat{P}}, Ptr{ArbFloat{P}}), &x, &y))
end
function (<=){P}(x::ArbFloat{P}, y::ArbFloat{P})
    return Bool(ccall(@libarb(libarb), Cint, (Ptr{ArbFloat{P}}, Ptr{ArbFloat{P}}), &x, &y))
end

(==){R<:Real,P}(x::ArbFloat{P}, y::R) = x == ArbFloat{P}(y)
(!=){R<:Real,P}(x::ArbFloat{P}, y::R) = x != ArbFloat{P}(y)
(<=){R<:Real,P}(x::ArbFloat{P}, y::R) = x <= ArbFloat{P}(y)
(>=){R<:Real,P}(x::ArbFloat{P}, y::R) = x >= ArbFloat{P}(y)
(<){R<:Real,P}(x::ArbFloat{P}, y::R) = x < ArbFloat{P}(y)
(>){R<:Real,P}(x::ArbFloat{P}, y::R) = x > ArbFloat{P}(y)

(==){R<:Real,P}(x::R, y::ArbFloat{P}) = ArbFloat{P}(x) == y
(!=){R<:Real,P}(x::R, y::ArbFloat{P}) = ArbFloat{P}(x) != y
(<=){R<:Real,P}(x::R, y::ArbFloat{P}) = ArbFloat{P}(x) <= y
(>=){R<:Real,P}(x::R, y::ArbFloat{P}) = ArbFloat{P}(x) >= y
(<){R<:Real,P}(x::R, y::ArbFloat{P}) = ArbFloat{P}(x) < y
(>){R<:Real,P}(x::R, y::ArbFloat{P}) = ArbFloat{P}(x) > y

# see predicates.jl for isequal{P}(x::ArbFloat{P}, y::ArbFloat{P})
isless{P}(x::ArbFloat{P}, y::ArbFloat{P}) = (x<y)


text:
function getTblMatch(;production_::Symbol=:S,token_::Union{Symbol,Int}=:return)
	tkn = -1
	if typeof(token_) == Int
		tkn = Symbol(uppercase(tks_names[token_]))

	else
		tkn = token_
	end

	tkn = Symbol(lowercase(string(tkn)))
	return table[tkn][getProd_idx(production_)]
end
function prodToString(prd::Production)
	output = ""

	output = prd.lexem*"-> "

	for i=1:length(prd.subprods)

		for elem in prd.subprods[i]
			if typeof(elem) == Symbol
				output = output*string(elem)
			else
				output = output*lowercase(string("'",elem,"'"))
			end
			output = output*" "
		end

		if i != length(prd.subprods)
			output = output*" | "
		end
	end
	return output*"\n"
end

function stringToProd(production::String)
	try
		prd = uppercase(production)
		splited = split(prd,"=")
		name = splited[1]
		enum = getProd_idx(Symbol(name))
		prd = splited[2]
		prd = split(prd," ")
		prd = map(Symbol,prd)

		output = Production(prd)
		output.enum = enum
		output.lexem = name
		return output
	catch
		return []
	end
end

function reverseRule(production_::Union{String,Production})
	production = -1
	if typeof(production_) == String
		production = stringToProd(production_)
	end
	output = Symbol[]
	sprods = production.subprods
	for p=length(sprods):-1:1
		output = vcat(output,sprods[p])
	end
	return output
end

function isProduction(p::Symbol)
	try
		getProd(p)
		return true
	catch
		return false
	end
end

function print_analyzer()
	println("Pilha: ",top(stack))
	println("Fita: ",head.lexem)
	println("Regra: ",rule)
	println()
end

function makeStepLine(lexem::String)
	len = length(lexem)
	df = 5 - len
	x = lexem

	for i=1:df
		x=x*" "
	end
	x=x*";"
	return x
end

stack = Stack(Union{Symbol,Int})

push!(stack,:S)

head = false
head = nextToken();
move_head = false
rule = ""
steps = "entrada;Pilha;Regra\n"

tabs=0
	while(head == false || head.categ_num != Int(EOF)) #head == false ---> end of line
		if head != false #Se não tiver um eol na fita

			newline = string(makeStepLine(head.lexem),string(collect(stack)))
			newline = @sprintf("%s",newline)
			if isProduction(top(stack)) #Se tiver uma produção no topo da pilha

				rule = getTblMatch(production_=top(stack),token_=head.categ_num)
				steps = steps*newline*";"*rule*"\n"

				tree_buff=""
				#=
				for i=1:tabs
					print("\t")
					#@printf("%s","\t")
					tree_buff=vcat(tree_buff,@sprintf("%s","\t")) #yeah, I know calling @printf and @sprintf twice is lame, But its bugging.
				end
				=#
				tree_buff=vcat(tree_buff,@sprintf("%10s"," "))
				println("          ",rule)
				tree_buff=vcat(tree_buff,@sprintf("%s%s","",rule*"\n"))
				write(tree_file,tree_buff)

				tabs=tabs+1
				r_rule = reverseRule(rule)
				pop!(stack)

				for i in r_rule
					push!(stack,i)
				end
				#sleep(1)
				#@show "loop"
			else

				if Int(eval(top(stack))) == Int(EPSILON)
					pop!(stack)
				else
					pop!(stack)
					tkn = printToken(head)
					write(tree_file,tkn*"\n")
					println(tkn)

					println()
					head = nextToken();
					tabs=0
					if head == false #Se for um eol, pula a fita de novo
						head = nextToken();
					end
				end
			end

			if typeof(head)!=Bool
				if head.categ_num == Int(EOF)
					if top(stack) == :S

						newline = string(makeStepLine(head.lexem),string(collect(stack)))
						newline = @sprintf("%s",newline)
						pop!(stack)
						steps = steps*newline*";S=epsilon"*"\n"

						newline = string(makeStepLine(head.lexem),string(collect(stack)))
						newline = @sprintf("%s",newline)
						steps = steps*newline*";---"*"\n"
					end
				end
			else
				head = nextToken();
			end
		end
		if manual
			readline(STDIN)
		end
	end

if (head.categ_num == Int(EOF))
	if (length(stack) == 0)
		info("Entrada aceita")
	elseif (length(stack)==1 && top(stack) == :S) #S -> epsilon

		pop!(stack)
		info("Entrada aceita")
	else
		info("Entrada rejeitada")
	end
else
		info("Entrada rejeitada")
end

z = open("../outputs/steps/$input_name.csv","w+")

#This bunch of replaces is to force Github to render the CSV. It's not elegant. There might be a better way to solve it.
steps=replace(steps,"\"","") #Para forçar o github a exibir o arquivo como CSV. Por algum motivo o CSV não é renderizado se tiver aspas
steps=replace(steps,"Any","") #Para forçar o github a exibir o arquivo como CSV. Por algum motivo o CSV não é renderizado se tiver aspas
steps=replace(steps,", ","،") #Para forçar o github a exibir o arquivo como CSV. Por algum motivo o CSV não é renderizado se tiver aspas
steps=replace(steps,";",",") #Para forçar o github a exibir o arquivo como CSV. Por algum motivo o CSV não é renderizado se tiver aspas

write(z,String(steps))
close(z)
steps = CSV.read("../outputs/steps/$input_name.csv";delim=",")

max_stars_repo_name: NaelsonDouglas/namelesscompiler

max_issues_repo_name: NaelsonDouglas/namelesscompiler

max_forks_repo_name: NaelsonDouglas/namelesscompiler


text:
# Writer the content of a into w.
function write{T}(w::WritableFile, a::Array{T})
	# If this is not provided, Base.IO write methods will write
	# arrays one element at a time.
	if isbits(T)
		write(w, pointer(a), length(a)*sizeof(T))
	else
		invoke(write, (IO, Array), w, a)
	end
end

# Writer the content of a into w.
if isdefined(:CartesianRange)
	function write{T,N,A<:Array}(w::WritableFile, a::SubArray{T,N,A})
		# This function is copied from Julia base/io.jl
		if !isbits(T) || stride(a,1)!=1
			return invoke(write, (Any, AbstractArray), s, a)
		end
		colsz = size(a,1)*sizeof(T)
		if N<=1
			return write(w, pointer(a, 1), colsz)
		else
			# cartesianmap was deprecated by v0.4.0-RC1.
			for idx in CartesianRange(tuple(1, size(a)[2:end]...))
				write(w, pointer(a, idx.I), colsz)
			end
			return colsz*Base.trailingsize(a,2)
		end
	end
else
	function write{T,N,A<:Array}(w::WritableFile, a::SubArray{T,N,A})
		# This function is copied from Julia base/io.jl
		if !isbits(T) || stride(a,1)!=1
			return invoke(write, (Any, AbstractArray), s, a)
		end
		colsz = size(a,1)*sizeof(T)
		if N<=1
			return write(w, pointer(a, 1), colsz)
		else
			cartesianmap((idxs...)->write(w, pointer(a, idxs), colsz),
				tuple(1, size(a)[2:end]...))
			return colsz*Base.trailingsize(a,2)
		end
	end
end

# Writer the byte b in w.
function write(w::WritableFile, b::UInt8)
	write(w, UInt8[b])
end

# Read and return a byte from f. Throws EOFError if there is no more byte to read.
function read(f::ReadableFile, ::Type{UInt8})
	# This function needs to be fast because readbytes, readall, etc.
	# uses it. Avoid function calls when possible.
	b = Array(UInt8, 1)
	read(f, b)
	b[1]
end


text:
function OilOutflowSen(Depth,Thp,Bsw,OilPVT,GasPVT,WaterPVT, Qmax;ϵ=0.0006, Ts=80, ∇T=1, GOR=false, di=2.99, Qn=10, Dn=50, tol=0.05, sg=0.6)

    # Create the range of Depth and Rates
    DepthRange=range(0, stop=Depth, length=Dn)
    if size(Qmax,1)==1
        Qrange=range(100, stop=Qmax, length=Qn)
    else
        Qrange=Qmax
        Qn=size(Qmax,1)
    end

    #Size of sensibilities
    Thpn=size(Thp,1)
    gorn=size(GOR,1)
    din=size(di,1)
    bswn=size(Bsw,1)

    #Create the interpolations function for PVT
    ## Oil
    Rs=LinearInterpolation(OilPVT.P.*1, OilPVT.Rs.*1)
    ρo=LinearInterpolation(OilPVT.P.*1, OilPVT.Rho.*1)
    βo=LinearInterpolation(OilPVT.P.*1, OilPVT.Bo.*1)
    μo=LinearInterpolation(OilPVT.P.*1, OilPVT.Muo.*1)
    Co=LinearInterpolation(OilPVT.P.*1, OilPVT.Co.*1)
    σo=LinearInterpolation(OilPVT.P.*1, OilPVT.Ten.*1)

    ## Gas
    Z=LinearInterpolation(GasPVT.P.*1, GasPVT.Z.*1)
    ρg=LinearInterpolation(GasPVT.P.*1, GasPVT.Rho.*1)
    βg=LinearInterpolation(GasPVT.P.*1, GasPVT.Bg.*1)
    μg=LinearInterpolation(GasPVT.P.*1, GasPVT.Mug.*1)
    Cg=LinearInterpolation(GasPVT.P.*1, GasPVT.Cg.*1)

    ## Water
    ρw=LinearInterpolation(WaterPVT.P.*1, WaterPVT.Rho.*1)
    βw=LinearInterpolation(WaterPVT.P.*1, WaterPVT.Bw.*1)
    μw=LinearInterpolation(WaterPVT.P.*1, WaterPVT.Muw.*1)
    Cw=LinearInterpolation(WaterPVT.P.*1, WaterPVT.Cw.*1)
    σw=LinearInterpolation(WaterPVT.P.*1, WaterPVT.Ten.*1)

    #Temperature
    T=LinearInterpolation([0, Depth], [Ts, Ts+Depth*(∇T/100)])

    #Calculate variables
    GOR==false ? Rsi=maximum(OilPVT.Rs) : Rsi=GOR

    leg=Array{String,1}(undef,Thpn*gorn*din*bswn)
    Pwf=zeros(Thpn*gorn*din*bswn,Qn)
    Pprofile=zeros(Dn,Qn,Thpn*gorn*din*bswn)

    count=0
    for k=1:Thpn  #Loop for Thp
        for m=1:gorn #Loop for Gor
            for n=1:din #Loop for Diameter
                for o=1:bswn
            count=count+1
        #Create Variables
    ql=zeros(Qn)     #Vector of Liquid rate
    qo=zeros(Qn)     #Vector of Oil Rate
    usl=zeros(Qn)    #Vector of liquid Velocities
    um=zeros(Dn)    #Vector of liquid Velocities
    P=zeros(Dn,Qn)   #Matrix of Pressures
    ∇P=zeros(Dn,Qn)  #Matriz of perssure gradient
    e=tol+0.01       #Reset tolerance value
    ∇Pnew=0.0        #Reset Gradient Value
    TransArea=π*(di[n]/(12*2))^2  #Transversal Area [Ft2]
    P[1,:].=Thp[k]       #Initial pressure is Thp
    ∇P[1,:].=(ρo(Thp[k])*(1-Bsw[o])+ρw(Thp[k])*Bsw[o])*(0.433/62.4)         #Initial PRessure gradient is the average density

    #General Loop for each Rate
    for j=1:Qn

        #Calculate Liquid Velocities
        ql[j]=Qrange[j]*(5.615/86400)           #Liquid Rate [ft3/s]
        qo[j]=(Qrange[j]*(1-Bsw[o]))*(5.615/86400) #Oil Rate [ft3/s]
        usl[j]=ql[j]/TransArea                  #Liquid Velocity [ft/s]

        #General Loop for each depth
        for i=2:Dn

         ∇Pguess=∇P[i-1,j]
          e=tol+0.01

          #While Iteration
           while e>=tol

             Pguess=∇Pguess*(DepthRange[i]-DepthRange[i-1])+P[i-1,j]                             #Estimate Pressure Guess [psi]
             FreeGas=(Rsi[m]-Rs[Pguess])*qo[j]                                            #Estimate FreeGas Rate [scf/d]
             usg= (4*FreeGas*Z(Pguess)*(460+T(DepthRange[i]))*14.7)/(Pguess*520*π*(di[n]/12)^2)  #Estimate Gas Velocity
             um[i]=usg+usl[j]  #Total Velocity
             σl= (σo(Pguess)*(1-Bsw[o])+σw(Pguess)*Bsw[o])
             ρl= (ρo(Pguess)*(1-Bsw[o])+ρw(Pguess)*Bsw[o])                                   #Liquid Density Lbm/ft3
             μl= (μo(Pguess)*(1-Bsw[o])+μw(Pguess)*Bsw[o])                                   #Liquid Viscosity [Cp]
             λg=usg/um[i]
             LBi=1.071-0.2218*(um[i]^2 /(di[n]/12))
             LBi>=0.13 ? LB=LBi : LB=0.13

                #Evaluate if Bubble flow Exist
                if LB>λg

                  yl=1-0.5*(1+(um[i]/0.8)-sqrt((1+(um[i]/0.8))^2-4*(usg/0.8)))

                else

                 Nvl= 1.938*usl[j]*(ρl/σl)^0.25                                    #Liquid Velocity Number
                 Nvg=1.938*usg*(ρl/σl)^0.25                                        #Gas Velocity Number
                 Nd=120.872*(di[n]/12)*(ρl/σl)^0.5                                    #Pipe Diameter Number
                 Nl=0.15726*μl*(1/(ρl*σl^3))^0.25                                  #Liquid Viscosity Number
                 #CNl=(0.0019+0.0322*Nl-0.6642*Nl^2+4.9951*Nl^3)/(1+10.0147*Nl-33.8696*Nl^2+277.2817*Nl^3) # original
                 CNl=(0.0019+0.0505*Nl-0.0929*Nl^2+0.061*Nl^3)   #pengtools
                 H=(CNl*Nvl*Pguess^0.1)/(Nvg^0.575 * 14.7^0.1 * Nd)
                 ψy=((0.0047+1123.32*H-729489.64*H^2)/(1+1097.1566*H-722153.97*H^2))^0.5
                 B=(Nvg*Nl^0.38)/(Nd^2.14)
                # ψ=(1.0886+69.9473*B-2334.3497*B^2+12896.683*B^3)/(1+53.4401*B-1517.9369*B^2+8419.8115*B^3)

                     if B>0.055
                         ψ=2.5714*B+1.5962
                    elseif B>0.025
                         ψ=-533.33*B^2+58.524*B+0.1171
                    else
                         ψ=27170*B^3-317.52*B^2+0.5472*B+0.9999
                    end

                     yl=ψy*ψ                                                                   #Liquid HoldUp
                end

                 ρg=(28.97*Pguess*sg)/(Z(Pguess)*10.73*T(DepthRange[i]))
                 ma= TransArea*(usl[j]*ρl+usg*ρg)*86400
                 Nre=(0.022*ma)/((di[n]/12)*μl^yl*μg(Pguess)^(1-yl))
                 f=(1/(-4*log((ϵ/3.7065)-(5.0452/Nre)*log((ϵ^1.1098/2.8257)+(7.149/Nre)^0.8981))))^2   #Friction factor
                 ρa=yl*ρl+(1-yl)*ρg
                 ∇Pnew=(1/144)*(ρa+((f*ma^2)/(7.413e10*(di[n]/12)^5*ρa))+((ρa*(um[i]-um[i-1]))/(32.2*2*(DepthRange[i]-DepthRange[i-1]))))
                  e=  abs(∇Pguess-∇Pnew)/∇Pnew
                  ∇Pguess=∇Pnew

            end #End While

            ∇P[i,j]=∇Pnew
            P[i,j]=∇P[i,j]*(DepthRange[i]-DepthRange[i-1])+P[i-1,j]

        end #End for Dn

    end #End for Qn
      Pwf[count,:].=P[end,:]
      Pprofile[:,:,count]=P
      leg[count]="Thp=$(Thp[k]) psi ; Bsw=$(Bsw[o]); Gor=$(Rsi[m]) scf/bbl; di=$(di[n]) in"
                end
            end
        end
    end
    return DepthRange, Qrange, Pwf, leg, Pprofile

    end #End Function


text:
using Compat
using Base.Test

v = 1

d = Dict{Int,Int}()
d[1] = 1
@test Compat.@Dict(1 => 1) == d
@test Compat.@Dict(1 => v) == d

@test typeof(@compat(Dict(1 => 1))) == Dict{Int,Int}
@test @compat(Dict(1 => 1)) == d
@test @compat(Dict(1 => v)) == d

ad = Dict{Any,Any}()
ad[1] = 1
@test Compat.@AnyDict(1 => 1) == ad
@test Compat.@AnyDict(1 => v) == ad

@test typeof(@compat(Dict{Any,Any}(1 => 1))) == Dict{Any,Any}
@test @compat(Dict{Any,Any}(1 => 1)) == ad
@test @compat(Dict{Any,Any}(1 => v)) == ad

td = Dict{Int,Float64}()
td[1] = 1.0

@test typeof(@compat(Dict{Int,Float64}(1 => 1))) == Dict{Int,Float64}
@test @compat(Dict{Int,Float64}(1 => 1)) == td
@test @compat(Dict{Int,Float64}(1 => v)) == td

@test @compat(Dict()) == Dict()
@test @compat(Dict{Any,Any}()) == Dict{Any,Any}()
if VERSION >= v"0.3.0-"
	@test @compat(Dict([(1, 1)])) == d
end

d2 = Dict{Symbol,Dict{Symbol,Int}}()
d2[:a] = Dict{Symbol,Int}()
d2[:a][:b] = 1
@test @compat(Dict(:a => Dict(:b => 1))) == d2

d = Dict(zip([1, 2], [3, 4]))
@test d == @compat Dict(1=>3, 2=>4)

@compat function f()
	a = :a
	b = Dict(:b => 1)
	Dict(a => b)
end
@test f() == d2

ns = length(d.slots)
@test length(sizehint!(d, ns + 1).slots) > ns

@test @compat split("a,b,,c", ',', limit=2) == ["a", "b,,c"]
@test @compat split("a,b,,c", ',', limit=2,keep=true) == ["a", "b,,c"]
@test @compat split("a,b,,c", ',', keep=false) == ["a", "b", "c"]
@test @compat split("a,b,,c", ',', keep=true) == ["a", "b", "", "c"]

@test @compat rsplit("a,b,,c", ',', limit=2) == ["a,b,", "c"]
@test @compat rsplit("a,b,,c", ',', limit=2,keep=true) == ["a,b,", "c"]
@test @compat rsplit("a,b,,c", ',', keep=false) == ["a", "b", "c"]
@test @compat rsplit("a,b,,c", ',', keep=true) == ["a", "b", "", "c"]

if VERSION < v"0.4.0-dev+1387"
    @test isdefined(Main, :AbstractString)
end


text:
struct ShaderCompilationError <: Exception
    msg
end

Base.showerror(io::IO, err::ShaderCompilationError) = print(io, "ShaderCompilationError:\n\n$(err.msg)")

# WARNING: this is type piracy.
function ShaderModule(device::Device, code::Vector{UInt8})
    size = cld(length(code), 4)
    ShaderModule(device, size, reinterpret(UInt32, resize!(copy(code), size * 4)))
end

spirv_code(shader::ShaderFile{<:FormatSPIRV}) = read(shader.file)
spirv_code(shader::ShaderFile) = spirv_code(compile(shader))

"""
    compile(shader)

Compile a shader file in text format to SPIR-V.
"""
function compile(shader::ShaderFile{FormatGLSL}; extra_flags=[], validate_spirv=true)::ShaderFile{FormatSPIRV}
    if !isfile(shader.file)
        throw(ArgumentError("File $(shader.file) does not exist"))
    end

    flags = ["-V"]
    validate_spirv && "--SPIRV-val" ∉ extra_flags ? push!(flags, "--SPIRV-val") : nothing
    dst = tempname()
    err = IOBuffer()
    try
        run(pipeline(`$glslangValidator $flags -o $dst $(shader.file)`, stdout=err))
    catch e
        e isa ProcessFailedException
        err_str = String(take!(err))
        throw(ShaderCompilationError(err_str))
    end
    ShaderFile(dst, FormatSPIRV(), shader.stage)
end

compile(shader::ShaderFile{FormatHLSL}) = compile(convert(ShaderFile{FormatGLSL}, shader), extra_flags=["-D"])


text:
###
### DynamicHMC backend - https://github.com/tpapp/DynamicHMC.jl
###

struct DynamicNUTS{AD, space} <: Hamiltonian{AD}
    n_iters   ::  Int   # number of samples
end
DynamicNUTS{AD}(n_iters::Int, space::Tuple) where {AD} = DynamicNUTS{AD, space}(n_iters)

"""
    DynamicNUTS(n_iters::Integer)

Dynamic No U-Turn Sampling algorithm provided by the DynamicHMC package.
To use it, make sure you have the DynamicHMC package installed.

"""
DynamicNUTS(args...) = DynamicNUTS{ADBackend()}(args...)
function DynamicNUTS{AD}(n_iters::Integer, space::Symbol...) where AD
    DynamicNUTS{AD}(n_iters, space)
end

getspace(::Type{<:DynamicNUTS{<:Any, space}}) where {space} = space
getspace(alg::DynamicNUTS{<:Any, space}) where {space} = space

function Sampler(alg::DynamicNUTS{T}, s::Selector=Selector()) where T <: Hamiltonian
  return Sampler(alg, Dict{Symbol,Any}(), s)
end

function sample(model::Model,
                alg::DynamicNUTS{AD},
                ) where AD

    spl = Sampler(alg)

    n = alg.n_iters
    samples = Array{Sample}(undef, n)
    weight = 1 / n
    for i = 1:n
        samples[i] = Sample(weight, Dict{Symbol, Any}())
    end

    vi = VarInfo(model)
    runmodel!(model, vi, SampleFromUniform())

    if spl.selector.tag == :default
        link!(vi, spl)
        runmodel!(model, vi, spl)
    end

    function _lp(x)
        value, deriv = gradient_logp(x, vi, model, spl)
        return ValueGradient(value, deriv)
    end

    chn_dynamic, _ = NUTS_init_tune_mcmc(FunctionLogDensity(length(vi[spl]), _lp), alg.n_iters)

    for i = 1:alg.n_iters
        vi[spl] = chn_dynamic[i].q
        samples[i].value = Sample(vi, spl).value
    end

    return Chain(0.0, samples)
end


text:
"""
    get_movers(idx_symbol::T,
               dire::T="up",
               chg::T="percent") where T<:AbstractString

Get top 10 (up or down) movers of an idx:
https://developer.tdameritrade.com/movers/apis/get/marketdata/{index}/movers

# Arguments
dire = "up" or "down"
chg = "value" or "percent"

# Examples
```julia
get_movers(raw"\$DJI")
get_movers("\$DJI", "up", "value")
```
"""
function get_movers(idx_symbol::T,
                    dire::T="up",
                    chg::T="percent") where T<:AbstractString

    uri = construct_api("marketdata/$idx_symbol/movers",
                        (direction=dire, change=chg)
                       )

    # Authorization is a header
    head = ["Authorization" => "Bearer "*AUTH_KEY.ACCESS_TOKEN]
    return @pipe HTTP.get(uri, head).body |>
                 JSON3.read
end