text:
struct CenteredRBM{M,Ov,Oh}
    rbm::M
    offset_v::Ov
    offset_h::Oh
    """
        CenteredRBM(rbm, offset_v, offset_h)

    Creates a centered RBM, with offsets `offset_v` (visible) and `offset_h` (hidden).
    See <http://jmlr.org/papers/v17/14-237.html> for details.
    The resulting model is *not* equivalent to the original `rbm`
    (unless offset_v = 0 and offset_h = 0).
    """
    function CenteredRBM(rbm::RBM, offset_v::AbstractArray, offset_h::AbstractArray)
        @assert size(visible(rbm)) == size(offset_v)
        @assert size(hidden(rbm)) == size(offset_h)
        return new{typeof(rbm),typeof(offset_v),typeof(offset_h)}(rbm, offset_v, offset_h)
    end
end

function CenteredRBM(
    visible::AbstractLayer, hidden::AbstractLayer, w::AbstractArray,
    offset_v::AbstractArray, offset_h::AbstractArray
)
    return CenteredRBM(RBM(visible, hidden, w), offset_v, offset_h)
end

"""
    CenteredRBM(visible, hidden, w)

Creates a centered RBM, with offsets initialized to zero.
"""
function CenteredRBM(visible::AbstractLayer, hidden::AbstractLayer, w::AbstractArray)
    offset_v = zeros(eltype(w), size(visible))
    offset_h = zeros(eltype(w), size(hidden))
    return CenteredRBM(RBM(visible, hidden, w), offset_v, offset_h)
end

CenteredRBM(rbm::RBM) = CenteredRBM(visible(rbm), hidden(rbm), weights(rbm))

RBMs.visible(rbm::CenteredRBM) = visible(RBM(rbm))
RBMs.hidden(rbm::CenteredRBM) = hidden(RBM(rbm))
RBMs.weights(rbm::CenteredRBM) = weights(RBM(rbm))
visible_offset(rbm::CenteredRBM) = rbm.offset_v
hidden_offset(rbm::CenteredRBM) = rbm.offset_h

"""
    RBM(centered_rbm::CenteredRBM)

Returns an (uncentered) `RBM` which neglects the offsets of `centered_rbm`.
The resulting model is *not* equivalent to the original `centered_rbm`.
To construct an equivalent model, use the function
`uncenter(centered_rbm)` instead (see [`uncenter`](@ref)).
Shares parameters with `centered_rbm`.
"""
RBMs.RBM(rbm::CenteredRBM) = rbm.rbm

function RBMs.energy(rbm::CenteredRBM, v::AbstractArray, h::AbstractArray)
    Ev = RBMs.energy(visible(rbm), v)
    Eh = RBMs.energy(hidden(rbm), h)
    Ew = RBMs.interaction_energy(rbm, v, h)
    return Ev .+ Eh .+ Ew
end

function RBMs.interaction_energy(centered_rbm::CenteredRBM, v::AbstractArray, h::AbstractArray)
    centered_v = v .- visible_offset(centered_rbm)
    centered_h = h .- hidden_offset(centered_rbm)
    return RBMs.interaction_energy(RBM(centered_rbm), centered_v, centered_h)
end

function RBMs.inputs_v_to_h(centered_rbm::CenteredRBM, v::AbstractArray)
    centered_v = v .- visible_offset(centered_rbm)
    return RBMs.inputs_v_to_h(RBM(centered_rbm), centered_v)
end

function RBMs.inputs_h_to_v(centered_rbm::CenteredRBM, h::AbstractArray)
    centered_h = h .- hidden_offset(centered_rbm)
    return RBMs.inputs_h_to_v(RBM(centered_rbm), centered_h)
end

function RBMs.free_energy(rbm::CenteredRBM, v::AbstractArray; β::Real = true)
    E = RBMs.energy(visible(rbm), v)
    inputs = RBMs.inputs_v_to_h(rbm, v)
    F = RBMs.free_energy(hidden(rbm), inputs; β)
    ΔE = RBMs.energy(RBMs.Binary(hidden_offset(rbm)), inputs)
    return E + F - ΔE
end

function RBMs.mean_h_from_v(rbm::CenteredRBM, v::AbstractArray; β::Real=1)
    inputs = RBMs.inputs_v_to_h(rbm, v)
    return RBMs.transfer_mean(hidden(rbm), inputs; β)
end

function RBMs.mean_v_from_h(rbm::CenteredRBM, h::AbstractArray; β::Real = true)
    inputs = RBMs.inputs_h_to_v(rbm, h)
    return RBMs.transfer_mean(visible(rbm), inputs; β)
end

function RBMs.∂free_energy(
    rbm::CenteredRBM, v::AbstractArray;
    wts = nothing, stats = RBMs.suffstats(visible(rbm), v; wts)
)
    inputs = RBMs.inputs_v_to_h(rbm, v)
    h = RBMs.transfer_mean(hidden(rbm), inputs)
    ∂v = RBMs.∂energy(visible(rbm), stats)
    ∂h = RBMs.∂free_energy(hidden(rbm), inputs; wts)
    ∂w = RBMs.∂interaction_energy(rbm, v, h; wts)
    return (visible = ∂v, hidden = ∂h, w = ∂w)
end

function RBMs.∂interaction_energy(
    rbm::CenteredRBM, v::AbstractArray, h::AbstractArray; wts = nothing
)
    centered_v = v .- visible_offset(rbm)
    centered_h = h .- hidden_offset(rbm)
    ∂w = RBMs.∂interaction_energy(RBM(rbm), centered_v, centered_h; wts)
    return ∂w
end

function RBMs.log_pseudolikelihood(rbm::CenteredRBM, v::AbstractArray; β::Real=1)
    return RBMs.log_pseudolikelihood(uncenter(rbm), v; β)
end


text:
################################################################################
# PDict
#
# Implementation of a persistent dictionary type.
#
# @author Noah C. Benson
#
# MIT License
# Copyright (c) 2020-2021 Noah C. Benson

# #PLinearDict #################################################################
# We construct the types in this file in an unfortunaetly complex way in order
# to reduce code duplication. The three kinds of persistent sets are very
# similar, just using different equality tests, so we can keep most of the code
# the same. This method is used for PDict, PLinearDict, PSet, and PLinearSet.
macro _plindict_code(name::Symbol, eq, h)
    return quote
        struct $name{K,V} <: AbstractPDict{K,V}
            keys::Union{Nothing, Vector{K}}
            values::Union{Nothing, Vector{V}}
        end
        function $name{K,V}() where {K,V}
            return $name{K,V}(nothing, nothing)
        end
        function $name{K,V}(kv::Pair{J,U}) where {K,V,J,U}
            return $name{K,V}(K[kv[1]], V[kv[2]])
        end
        function $name{K,V}(kv::Tuple{J,U}) where {K,V,J,U}
            return $name{K,V}(K[kv[1]], V[kv[2]])
        end
        function $name{K,V}(d::$name{K,V}) where {K,V}
            return d
        end
        function $name{K,V}(d::AbstractDict) where {K,V}
            return reduce(push, d, init=$name{K,V}())
        end
        function $name{K,V}(kv::Pair, ps::Pair...) where {K,V}
            return reduce(push, ps, init=$name{K,V}(kv))
        end
        function $name{K,V}(kv::Tuple, ps::Tuple...) where {K,V}
            return reduce(push, ps, init=$name{K,V}(kv))
        end
        function $name{K,V}(itr) where {K,V}
            return reduce(push, itr, init=$name{K,V}())
        end
        # Document the equal/hash types.
        equalfn(::Type{T}) where {T <: $name} = $eq
        hashfn(::Type{T}) where {T <: $name} = $h
        equalfn(::$name) = $eq
        hashfn(::$name) = $h
        # Generic construction.
        $name() = $name{Any,Any}()
        $name(::Tuple{}) = $name{Any,Any}()
        $name(p::Pair{K,V}) where {K,V} = $name{K,V}(p)
        $name(d::AbstractDict{K,V}) where {K,V} = $name{K,V}(d)
        $name(ps::Pair...) = begin
            ps = _to_pairs(ps)
            T = Base.eltype(ps)
            K = T.parameters[1]
            V = T.parameters[2]
            return $name{K,V}(ps...)
        end
        $name(itr) = begin
            ps = _to_pairs(itr)
            T = Base.eltype(ps)
            K = T.parameters[1]
            V = T.parameters[2]
            return $name{K,V}(ps...)
        end
        # Base methods.
        Base.empty(s::$name{K,V}, ::Type{J}=K, ::Type{U}=V) where {K,V,J,U} = $name{J,U}()
        Base.length(s::$name) = begin
            ks = getfield(s, :keys)
            return (ks === nothing ? 0 : length(ks))
        end
        Base.IteratorEltype(::Type{$name{K,V}}) where {K,V} = Base.HasEltype()
        Base.eltype(::Type{$name{K,V}}) where {K,V} = Pair{K,V}
        Base.eltype(::$name{K,V}) where {K,V} = Pair{K,V}
        Base.IteratorSize(::Type{$name{K,V}}) where {K,V} = Base.HasLength()
        Base.iterate(u::$name{K,V}, ii::Int) where {K,V} = begin
            (ii > length(u)) && return nothing
            ks = getfield(u, :keys)::Vector{K}
            vs = getfield(u, :values)::Vector{V}
            return (Pair{K,V}(ks[ii], vs[ii]), #(@inbounds ks[ii], @inbounds vs[kk]),
                    ii + 1)
        end
        Base.iterate(u::$name{K,V}) where {K,V} = iterate(u, 1)
        Base.get(u::$name{K,V}, kk, df) where {K,V} = begin
            ks = getfield(u, :keys)
            (ks === nothing) && return df
            vs = getfield(u, :values)::Vector{V}
            for ii in 1:length(ks)
                k = ks[ii]
                $eq(k,kk) && return vs[ii]
            end
            return df
        end
        Base.in(kv::Pair, u::$name{K,V}, eqfn::F) where {K,V,F<:Function} = begin
            ks = getfield(u, :keys)
            (ks === nothing) && return false
            vs = getfield(u, :values)::Vector{V}
            for ii in 1:length(ks)
                k = ks[ii]
                $eq(k,kv[1]) && return eqfn(kv[2], vs[ii])
            end
            return false
        end
        Base.in(x::Pair, u::$name{K,V}) where {K,V} = in(x, u, (==))
        # And the Air methods.
        Air.push(u::$name{K,V}, kv::Pair) where {K,V} = begin
            ks = getfield(u, :keys)
            (ks === nothing) && return $name{K,V}(kv)
            vs = getfield(u, :values)::Vector{V}
            for ii in 1:length(ks)
                k = ks[ii]
                if $eq(k, kv[1])
                    if $eq(vs[ii], kv[2])
                        return u
                    else
                        return $name{K,V}(ks, setindex(vs, kv[2], ii))
                    end
                end
            end
            return $name{K,V}(push(ks, kv[1]), push(vs, kv[2]))
        end
        Air.setindex(u::$name{K,V}, v, k) where {K,V} = push(u, k => v)
        Air.delete(u::$name{K,V}, x::J) where {K,V,J} = begin
            ks = getfield(u, :keys)
            (ks === nothing) && return u
            vs = getfield(u, :values)::Vector{V}
            n = length(ks)
            if n == 1
                if $eq(ks[1], x)
                    return $name{K,V}(nothing, nothing)
                end
            else
                for ii in 1:n
                    if $eq(ks[ii], x)
                        return $name{K,V}(delete(ks, ii), delete(vs, ii))
                    end
                end
            end
            return u
        end
        Base.propertynames(u::$name) = (:count, :keys, :values)
        Base.getproperty(u::$name{K,V}, s::Symbol) where {K,V} = begin
            if s == :count
                return length(u)
            elseif s == :keys
                ks = getfield(u, :keys)
                return ks === nothing ? () : tuple(ks...)
            elseif s == :values
                vs = getfield(u, :values)
                return vs === nothing ? () : tuple(vs...)
            elseif s == :pairs
                ks = getfield(u, :keys)
                (ks === nothing) && return ()
                vs = getfield(u, :values)::Vector{V}
                return Tuple([k => v for (k,v) in zip(ks,vs)])
            else
                throw(ArgumentError("no such property $s of type $(typeof(u))"))
            end
        end
    end |> esc
end
# Declare the linear types.
@_plindict_code PLinearDict isequal hash
@_plindict_code PIdLinearDict (===) objectid
# And document them.
@doc """
    PLinearDict{K,V}

A persistent dict type that uses a simple array representation internally;
accordingly `PLinearDict` does not peform persistent operations such as
`setindex`, `push`, and `delete` efficiently. It is intended to be used by Air
internally as a sub-dict that stores items whose hash values are identical.

See also: [`PIdLinearDict`](@ref), [`PDict`](@ref).
""" PLinearDict
@doc """
    PIdLinearDict{K,V}

An identity-based persistent dict type that uses a simple array representation
internally; accordingly `PIdLinearDict` does not peform persistent operations
such as `setindex`, `push`, and `delete` efficiently. It is intended to be used
by Air internally as a sub-dict that stores items whose hash values are
identical.

See also: [`PLinearDict`](@ref), [`PIdDict`](@ref).
""" PIdLinearDict

################################################################################
# PDict
macro _pdict_code(name::Symbol, eq, h, lindict)
    return quote
        struct $name{K,V} <: AbstractPDict{K,V}
            count::Int
            root::PTree{$lindict{K,V}}
        end
        function $name{K,V}() where {K,V}
            return $name{K,V}(0, PTree{$lindict{K,V}}())
        end
        function $name{K,V}(d::$name{K,V}) where {K,V}
            return d
        end
        function $name{K,V}(d::AbstractDict) where {K,V}
            return reduce(push, d, init=$name{K,V}())
        end
        function $name{K,V}(ps::Pair...) where {K,V}
            return reduce(push, ps, init=$name{K,V}())
        end
        function $name{K,V}(ps::Tuple...) where {K,V}
            return reduce(push, ps, init=$name{K,V}())
        end
        function $name{K,V}(itr) where {K,V}
            return reduce(push, itr, init=$name{K,V}())
        end
        # Document the equal/hash types.
        equalfn(::Type{T}) where {T <: $name} = $eq
        hashfn(::Type{T}) where {T <: $name} = $h
        equalfn(::$name) = $eq
        hashfn(::$name) = $h
        # Generic construction.
        $name() = $name{Any,Any}()
        $name(::Tuple{}) = $name{Any,Any}()
        $name(p::Pair{K,V}) where {K,V} = $name{K,V}(p)
        $name(d::AbstractDict{K,V}) where {K,V} = $name{K,V}(d)
        $name(ps::Pair...) = begin
            ps = _to_pairs(ps)
            T = Base.eltype(ps)
            K = T.parameters[1]
            V = T.parameters[2]
            return $name{K,V}(ps...)
        end
        $name(itr) = begin
            ps = _to_pairs(itr)
            T = Base.eltype(ps)
            K = T.parameters[1]
            V = T.parameters[2]
            return $name{K,V}(ps...)
        end
        # Base methods.
        Base.empty(s::$name{K,V}, ::Type{J}=K, ::Type{U}=V) where {K,V,J,U} = $name{J,U}()
        Base.length(s::$name) = getfield(s, :count)
        Base.IteratorEltype(::Type{$name{K,V}}) where {K,V} = Base.HasEltype()
        Base.eltype(::Type{$name{K,V}}) where {K,V} = Pair{K,V}
        Base.eltype(::$name{K,V}) where {K,V} = Pair{K,V}
        Base.IteratorSize(::Type{$name{K,V}}) where {K,V} = Base.HasLength()
        Base.iterate(u::$name{K,V}) where {K,V} = begin
            x = itervals(getfield(u, :root))
            (x === nothing) && return nothing
            (rootel,rootii) = x
            (el,reliter) = iterate(rootel)
            return (el,(rootii,reliter,length(rootel)))
        end
        Base.iterate(u::$name{K,V}, tup::Tuple{HASH_T,Int,Int}) where {K,V} = begin
            (rootii, reliter, rellen) = tup
            root = getfield(u, :root)
            if reliter < rellen
                # There are more in that node:
                rootel = get(root, rootii, nothing)::$lindict{K,V}
                (el,reliter) = iterate(rootel, reliter)
                return (el,(rootii, reliter, rellen))
            end
            x = itervals(root, rootii)
            (x === nothing) && return nothing
            (rootel,rootii) = x
            (el,reliter) = iterate(rootel)
            return (el,(rootii,reliter,length(rootel)))
        end
        Base.get(u::$name{K,V}, k, df) where {K,V} = begin
            ld = get(getfield(u, :root), $h(k), nothing)
            return (ld === nothing ? df : get(ld, k, df))
        end
        Base.in(kv::Pair, u::$name{K,V}, eqfn::F) where {K,V,F<:Function} = begin
            ld = get(getfield(u, :root), $h(kv[1]), nothing)
            return (ld === nothing ? false : in(kv, ld, eqfn))
        end
        Base.in(x::Pair, u::$name{K,V}) where {K,V} = in(x, u, (==))
        # And the Air methods.
        Air.push(u::$name{K,V}, x::Pair) where {K,V} =
            setindex(u, x.second, x.first)
        Air.setindex(u::$name{K,V}, v, k) where {K,V} = begin
            hh = $h(k)
            root = getfield(u, :root)
            ld0 = get(root, hh, nothing)
            if ld0 === nothing
                return $name{K,V}(
                    getfield(u, :count) + 1,
                    setindex(root, $lindict{K,V}(Pair{K,V}(k,v)), hh))
            else
                ld1 = setindex(ld0, v, k)
                (ld1 === ld0) && return u
                n = getfield(u, :count) + length(ld1) - length(ld0)
                return $name{K,V}(n, setindex(root, ld1, hh))
            end
        end
        Air.delete(u::$name{K,V}, x::J) where {K,V,J} = begin
            hh = $h(x)
            root = getfield(u, :root)
            ld0 = get(root, hh, nothing)
            (ld0 === nothing) && return u
            ld1 = delete(ld0, x)
            (ld0 === ld1) && return u
            if length(ld1) == 0
                return $name{K,V}(getfield(u, :count) - 1, delete(root, hh))
            else
                return $name{K,V}(getfield(u, :count) - 1, setindex(root, ld1, hh))
            end
        end
    end |> esc
end
# Declare the types.
(@_pdict_code PIdDict (===) objectid PIdLinearDict)
(@_pdict_code PDict isequal hash PLinearDict)
# Document the types.
@doc """
    PDict{K,V}

A dictionay type roughly equivalent to `Dict{K,V}` but that stores data using a
persistent hash array mapped trie system that allows for efficient persistent
operations. These operations are represented by non-mutating versions of the
standard dictionary functions, such as `setindex` in place of `setindex!`,
`push` in place of `push!`, and `pop` in place of `pop!`. These operations are
performed in `O(log n)` time, and minimal data duplication is performed in
update operations.

See also: [`PIdDict`](@ref), `Dict`, [`push`](@ref), [`pop`](@ref),
[`setindex`](@ref), [`delete`](@ref), [`insert`](@ref), [`getpair`](@ref).

# Examples

```@meta
DocTestSetup = quote
    using Air
end
```

```jldoctest; filter=r"PDict{Any, ?Any}"
julia> PDict()
PDict{Any,Any}()
```

```jldoctest; filter=r"PDict{Symbol, ?Real} with 3 entries:"
julia> PDict(:a => 1, :b => 2, :c => 12.8)
PDict{Symbol,Real} with 3 entries:
  :c => 12.8
  :a => 1
  :b => 2
```

```jldoctest; filter=r"PDict{Symbol, ?Float64} with [34] entries:"
julia> d = PDict{Symbol,Float64}(:a => 1, :b => 2, :c => 12.8)
PDict{Symbol,Float64} with 3 entries:
  :c => 12.8
  :a => 1.0
  :b => 2.0

julia> :b in keys(d)
true

julia> d[:c]
12.8

julia> push(d, :d => 0.1)
PDict{Symbol,Float64} with 4 entries:
  :d => 0.1
  :c => 12.8
  :a => 1.0
  :b => 2.0
```
""" PDict
@doc """
    PIdDict{K,V}

A dictionay type roughly equivalent to `IdDict{K,V}` but that stores data using
a persistent hash array mapped trie system that allows for efficient persistent
operations. These operations are represented by non-mutating versions of the
standard dictionary functions, such as `setindex` in place of `setindex!`,
`push` in place of `push!`, and `pop` in place of `pop!`. These operations are
performed in `O(log n)` time, and minimal data duplication is performed in
update operations.

See also: [`PDict`](@ref), `Base.IdDict`, [`push`](@ref), [`pop`](@ref),
[`setindex`](@ref), [`delete`](@ref), [`insert`](@ref), [`getpair`](@ref).

# Examples

```@meta
DocTestSetup = quote
    using Air
end
```

```jldoctest; filter=r"PIdDict{Any, ?Any}\\(\\)"
julia> PIdDict()
PIdDict{Any,Any}()

```jldoctest; filter=r"PIdDict{Symbol, ?Real} with 3 entries:"
julia> PIdDict(:a => 1, :b => 2, :c => 12.8)
PIdDict{Symbol,Real} with 3 entries:
  :c => 12.8
  :a => 1
  :b => 2
```

```jldoctest; filter=r"PIdDict{Symbol, ?Float64} with [34] entries:"
julia> d = PIdDict{Symbol,Float64}(:a => 1, :b => 2, :c => 12.8)
PIdDict{Symbol,Float64} with 3 entries:
  :c => 12.8
  :a => 1.0
  :b => 2.0

julia> :b in keys(d)
true

julia> d[:c]
12.8

julia> push(d, :d => 0.1)
PIdDict{Symbol,Float64} with 4 entries:
  :d => 0.1
  :c => 12.8
  :a => 1.0
  :b => 2.0
""" PIdDict

# Export the relevant symbols.
export PDict, PIdDict

ext: jl


text:
using FactCheck
using SimpleHttpCommon
const S = SimpleHttpCommon
const N = S.N
const N_TYPE = S.N_TYPE
const TEST_DIR = dirname(@__FILE__)
const CRLF = S.CRLF
const STR_TYPE = S.STR_TYPE

function buffer(data = N; start=false)
    b = IOBuffer()
    if data != N
        write(b, data)
    end

    if start
        seekstart(b)
    end

    b
end

function mkres()
    HTTPResponse(IOSocket(buffer()))
end

type ExceptionNotRaised <: Exception
    info
end

type NOT_DEFINED end

function factcatch{T}(fn::Function, extype::Type{T}, einfo=N; kwargs...)
    try
        fn()
        throw(ExceptionNotRaised((extype, einfo)))
    catch e
        if !isa(e, extype)
            rethrow()
        end

        for (key, expected_value) in kwargs
            value = get(e.info, key, NOT_DEFINED)
            @fact value --> expected_value string("Key value mismatch $extype.info[$(repr(key))]. Info is: $(repr(e.info))")
        end
    end
end

fre(fn::Function, einfo=N; kwargs...) = factcatch(fn, S.ResponseException; kwargs...)
frqe(fn::Function, einfo=N; kwargs...) = factcatch(fn, S.RequestException; kwargs...)

const str = UTF8String

function writecrlf(io, data="")
    s = 0
    if data != ""
        s += Base.write(io, data)
    end
    s + Base.write(io, "\r\n")
end

const wclf = writecrlf


text:
using RCall
R"""
library("MPsychoR")
data("Rmotivation")
ind <- grep("hyb", colnames(Rmotivation))
HybMotivation <- na.omit(Rmotivation[,ind])
"""
@rget HybMotivation
α = cronalpha(HybMotivation)

@testset "reliabilitytests" begin
    @test round(α,digits=2) == 0.82
end


text:
export ARC_generic

# One dimensionnal ARC algorithm.
# We approximate our function h by a 2nd degree polynomial (Taylor) q enriched
# with a cubic term.
# If c is a good approximation of h then we move towards the minimizer of c
# Otherwise we reduce Δ in order to increase the weight of the cubic term.
# c(d) = h(t) + h'(t) * d + 0.5 * h''(t) * d² + (1/(3*Δ))*|d|³
#
# Specific key word arguments
# eps1 :: Float64. Below that treshold we don't move and reduce Δ.
#                  Treshold to determines if c is a bad approximation of h
# eps2 :: Float64. Over that treshold we augment the size of Δ
#                  Treshold to determines if c is a very good approximation of h
# red :: Float64. Rate at which we reduce Δ if we have a bad approximation
# aug :: Float64. Rate at which we augment Δ if we have a bad approximation
# Δ :: Float64. Size of the trust region at the begining of the algorithm.

function ARC_generic(h         :: AbstractNLPModel,
                     nlpstop   :: AbstractStopping;
                     verbose   :: Bool = true,
                     eps1      :: Float64 = 0.25,
                     eps2      :: Float64 = 0.75,
                     red       :: Float64 = 0.2,
                     aug       :: Float64 = 5.0,
                     Δ         :: Float64 = 0.5,
                     direction :: Symbol = :Nwt)

    Δ = [Δ]
    t = h.meta.x0; iter = 0;        # We establish our starting point t
    fₖ = obj(h, t);
    f₀ = copy(fₖ)
    gₖ = grad(h, t);
    g₀ = copy(gₖ)

    # H denotes the (approximation of the) second derivative
    if direction == :Nwt
      H = hess(h, t)
    elseif direction == :Sec || direction == :SecA
      H = [1.0]
    end

    OK = update_and_start!(nlpstop, x = t, fx = fₖ, gx = gₖ, g0 = g₀, Hx = H)

    verbose &&
        @printf(" iter  t         gₖ          Δ        pred         ared\n")
    verbose &&
        @printf(" %4d %7.2e  %7.2e  %7.2e \n", iter, t[1], gₖ[1], Δ[1])
    # We loop until we have a minimizer or we have reached the maximum number of
    # iterations.
    while !OK

        d = ARC_step_computation(H, gₖ, Δ) # We find the direction in which we
                                                # move
        # Numerical reduction computation
        ftestTR = obj(h, t + d)  # Value of h and h' at t + d
        gtestTR = grad(h, t + d)

        # We check to see if we have a good approximation of h using the ratio
        # of the actual reduction and the predicted reduction.
        (pred, ared, ratio) =
            pred_ared_computation(gₖ, fₖ, H, d, ftestTR, gtestTR)

        if (ratio .< eps1)
            # Bad approximation of h. We make the cubic term more prevalant
            Δ = red * Δ
            add_stop!(nlpstop.meta)
        else
            # Good approximation we move towards the minimizer of c
            (t, fₖ, gₖ, H) =
                step_computation(direction, h, t, d, gₖ, fₖ, ftestTR, gtestTR)

            if ratio > eps2
                # Very good approximation of h. We reduce the importance of the
                # cubic term.
                Δ = aug * Δ
            end
            OK = update_and_stop!(nlpstop, x = t, fx = fₖ, gx = gₖ, Hx = H)
        end

        iter += 1
        verbose && @printf(" %4d %7.2e  %7.2e  %7.2e %7.2e %7.2e\n",
                            iter, t[1], gₖ[1], Δ[1], pred[1], ared[1])
    end

    optimal = nlpstop.meta.optimal
    return optimal, nlpstop
end


text:
using NetworkInference
using Test
using DelimitedFiles

# These tests use a dataset generated from the 10-node Yeast1 network from http://gnw.sourceforge.net/
# GeneNetWeaver: In silico benchmark generation and performance profiling of network inference methods.
# Schaffter T, Marbach D, and Floreano D. Bioinformatics, 27(16):2263-70, 2011.

println("Getting nodes...")
data_path = joinpath(dirname(@__FILE__), "data")
data_file_path = joinpath(data_path, "yeast1_10_data.txt")
nodes = get_nodes(data_file_path)

println("Inferring networks...")
mi_network = InferredNetwork(MINetworkInference(), nodes)
clr_network = InferredNetwork(CLRNetworkInference(), nodes)
puc_network = InferredNetwork(PUCNetworkInference(), nodes)
pidc_network = InferredNetwork(PIDCNetworkInference(), nodes)

# The benchmarks were inferred prior to the implementation of this package. MI, PUC and PIDC were
# inferred using scripts based on InformationMeasures.jl and CLR was inferred using MINET:
# https://bioconductor.org/packages/release/bioc/html/minet.html

mi_benchmark = readdlm(joinpath(data_path, "mi.txt"))
clr_benchmark = readdlm(joinpath(data_path, "clr.txt"))
puc_benchmark = readdlm(joinpath(data_path, "puc.txt"))
pidc_benchmark = readdlm(joinpath(data_path, "pidc.txt"))

# Compare a few selected edges throughout the inferred network
for i in (1, 5, 10, 20, 40)
    @test mi_network.edges[i].weight ≈ mi_benchmark[2*i, 3] atol = 0.0001
    println("MI network inference $i passed")
    @test clr_network.edges[i].weight ≈ clr_benchmark[2*i, 3] atol = 0.0001
    println("CLR network inference $i passed")
    @test puc_network.edges[i].weight ≈ puc_benchmark[2*i, 3] atol = 0.0001
    println("PUC network inference $i passed")
    @test pidc_network.edges[i].weight ≈ pidc_benchmark[2*i, 3] atol = 0.0001
    println("PIDC network inference $i passed")
end

# Test infer_network
mi_network = infer_network(data_file_path, MINetworkInference())
@test mi_network.edges[40].weight ≈ mi_benchmark[2*40, 3] atol = 0.0001
println("infer_network passed")

# Test get_adjacency_matrix
present_node1, present_node2 = mi_network.edges[40].nodes
absent_node1, absent_node2 = mi_network.edges[41].nodes
adj_matrix_absolute, labels_to_ids, ids_to_labels = get_adjacency_matrix(mi_network, mi_network.edges[40].weight, absolute = true)
adj_matrix_percentage, labels_to_ids, ids_to_labels = get_adjacency_matrix(mi_network, 40 / length(mi_network.edges))
for adj_matrix in [adj_matrix_absolute, adj_matrix_percentage]
    @test adj_matrix[labels_to_ids[present_node1.label], labels_to_ids[present_node2.label]] === true
    @test adj_matrix[labels_to_ids[present_node2.label], labels_to_ids[present_node1.label]] === true
    @test adj_matrix[labels_to_ids[absent_node1.label], labels_to_ids[absent_node2.label]] === false
    @test adj_matrix[labels_to_ids[absent_node2.label], labels_to_ids[absent_node1.label]] === false
end
println("get_adjacency_matrix passed")

# Test read_network_file
test_network_filepath = joinpath(data_path, "mi.txt")
read_network = read_network_file(test_network_filepath)
@test length(read_network.edges) == 45
reference_nodes = ["G1", "G2", "G3", "G4", "G5", "G6", "G7", "G8", "G9", "10"]
@test Set([x.label for x in read_network.nodes]) == Set(reference_nodes)
# Compare a few edges
@test read_network.edges[1].weight ≈ 0.2972651133683284 atol = 0.00001
@test Set([x.label for x in read_network.edges[1].nodes]) == Set(["G1", "G4"])
@test read_network.edges[15].weight ≈ 0.07202931598519413 atol = 0.00001
@test Set([x.label for x in read_network.edges[15].nodes]) == Set(["G3", "G9"])
@test read_network.edges[34].weight ≈ 0.02896110920075312 atol = 0.00001
@test Set([x.label for x in read_network.edges[34].nodes]) == Set(["G2", "G8"])
@test read_network.edges[35].weight ≈ 0.02848989653696199 atol = 0.00001
@test Set([x.label for x in read_network.edges[35].nodes]) == Set(["G9", "10"])
println("read_network_file passed")

# These tests will only run if the EmpiricalBayes package exists:
include("empirical_bayes_glue_tests.jl")


text:
#
# Transforming Sampler
#
# applies a reversible function to the sample
#

type TransformSampler <: ModifyingSampler
  subsampler::Sampler
  fn::Function
  invfn::Function
  function TransformSampler(subsampler::Sampler, fn::Function, invfn::Function)
	  new(subsampler, fn, invfn)
  end
end

function sample(s::TransformSampler, support, cc::ChoiceContext)
	x, trace = sample(s.subsampler, map(s.invfn, support), cc)
	s.fn(x), Dict{Symbol, Any}(:sub=>trace)
end

amendtrace(s::TransformSampler, trace, x) = amendtrace(s.subsampler, trace[:sub], s.invfn(x))


text:
@testset "1020.number-of-enclaves.jl" begin
    @test num_of_enclaves_bfs([[0, 0, 0, 0], [1, 0, 1, 0], [0, 1, 1, 0], [0, 0, 0, 0]]) == 3
    @test num_of_enclaves_bfs([[0, 1, 1, 0], [0, 0, 1, 0], [0, 0, 1, 0], [0, 0, 0, 0]]) == 0
    @test num_of_enclaves_bfs([[0], [1], [1], [0], [0]]) == 0
    @test num_of_enclaves_bfs([
        [0, 0, 0, 1, 1, 1, 0, 1, 0, 0],
        [1, 1, 0, 0, 0, 1, 0, 1, 1, 1],
        [0, 0, 0, 1, 1, 1, 0, 1, 0, 0],
        [0, 1, 1, 0, 0, 0, 1, 0, 1, 0],
        [0, 1, 1, 1, 1, 1, 0, 0, 1, 0],
        [0, 0, 1, 0, 1, 1, 1, 1, 0, 1],
        [0, 1, 1, 0, 0, 0, 1, 1, 1, 1],
        [0, 0, 1, 0, 0, 1, 0, 1, 0, 1],
        [1, 0, 1, 0, 1, 1, 0, 0, 0, 0],
        [0, 0, 0, 0, 1, 1, 0, 0, 0, 1],
    ]) == 3
end


text:
"""
Callback which can be used to catch OpenGL errors.
"""
function openglerrorcallback(
        source::GLenum, typ::GLenum,
        id::GLuint, severity::GLenum,
        length::GLsizei, message::Ptr{GLchar},
        userParam::Ptr{Void}
    )
    errormessage = """
         ________________________________________________________________
        |
        | OpenGL Error!
        | source: $(GLENUM(source).name) :: type: $(GLENUM(typ).name)
        |  $(Compat.String(message, length))
        |________________________________________________________________
    """
    output = typ == GL_DEBUG_TYPE_ERROR ? error : info
    output(errormessage)
    nothing
end

global const _openglerrorcallback = cfunction(
    openglerrorcallback, Void,
    (GLenum, GLenum,GLuint, GLenum, GLsizei, Ptr{GLchar}, Ptr{Void})
)

"""
Screen constructor cnstructing a new screen from a parant screen.
"""
function Screen(
        parent::Screen;
        name = gensym(parent.name),
        area = parent.area,
        children::Vector{Screen} = Screen[],
        inputs::Dict{Symbol, Any} = copy(parent.inputs),
        renderlist::Tuple = (),
        hidden::Bool = parent.hidden,
        glcontext::GLContext = parent.glcontext,
        cameras = Dict{Symbol, Any}(),
        position = Vec3f0(2),
        lookat = Vec3f0(0),
        color = RGBA{Float32}(1,1,1,1)
    )
    screen = Screen(name,
        area, parent, children, inputs,
        renderlist, hidden, color,
        cameras, glcontext
    )
    pintersect = const_lift(x->intersect(zeroposition(value(parent.area)), x), area)
    relative_mousepos = const_lift(inputs[:mouseposition]) do mpos
        Point{2, Float64}(mpos[1]-value(pintersect).x, mpos[2]-value(pintersect).y)
    end
    #checks if mouse is inside screen and not inside any children
    insidescreen = droprepeats(const_lift(isinside, screen, relative_mousepos))
    merge!(screen.inputs, Dict(
        :mouseinside 	=> insidescreen,
        :mouseposition 	=> relative_mousepos,
        :window_area 	=> area
    ))
    # creates cameras for the sceen with the new inputs

    push!(parent.children, screen)
    screen
end

"""
On OSX retina screens, the window size is different from the
pixel size of the actual framebuffer. With this function we
can find out the scaling factor.
"""
function scaling_factor(window::Vec{2, Int}, fb::Vec{2, Int})
    (window[1] == 0 || window[2] == 0) && return Vec{2, Float64}(1.0)
    Vec{2, Float64}(fb) ./ Vec{2, Float64}(window)
end

"""
Correct OSX scaling issue and move the 0,0 coordinate to left bottom.
"""
function corrected_coordinates(
        window_size::Signal{Vec{2,Int}},
        framebuffer_width::Signal{Vec{2,Int}},
        mouse_position::Vec{2,Float64}
    )
    s = scaling_factor(window_size.value, framebuffer_width.value)
    Vec(mouse_position[1], window_size.value[2] - mouse_position[2]) .* s
end

"""
Standard set of callback functions
"""
function standard_callbacks()
    Function[
        window_open,
        window_size,
        window_position,
        keyboard_buttons,
        mouse_buttons,
        dropped_files,
        framebuffer_size,
        unicode_input,
        cursor_position,
        scroll,
        hasfocus,
        entered_window,
    ]
end

"""
Tries to create sensible context hints!
Taken from lessons learned at:
[GLFW](http://www.glfw.org/docs/latest/window.html)
"""
function standard_context_hints(major, minor)
    # this is spaar...Modern OpenGL !!!!
    major < 3 && error("OpenGL major needs to be at least 3.0. Given: $major")
    # core profile is only supported for OpenGL 3.2+ (and a must for OSX, so
    # for the sake of homogenity, we try to default to it for everyone!)
    if (major > 3 || (major == 3 && minor >= 2 ))
        profile = GLFW.OPENGL_CORE_PROFILE
    else
        profile = GLFW.OPENGL_ANY_PROFILE
    end
    [
        (GLFW.CONTEXT_VERSION_MAJOR, major),
        (GLFW.CONTEXT_VERSION_MINOR, minor),
        (GLFW.OPENGL_FORWARD_COMPAT, GL_TRUE),
        (GLFW.OPENGL_PROFILE, profile)
    ]
end

"""
Takes half the resolution of the primary monitor.
This should make for sensible defaults!
"""
function standard_screen_resolution()
    w, h = primarymonitorresolution()
    (div(w,2), div(h,2)) # half of total resolution seems like a good fit!
end

"""
Standard window hints for creating a plain context without any multisampling
or extra buffers beside the color buffer
"""
function standard_window_hints()
    [
        (GLFW.SAMPLES,      0),
        (GLFW.DEPTH_BITS,   0),

        (GLFW.ALPHA_BITS,   8),
        (GLFW.RED_BITS,     8),
        (GLFW.GREEN_BITS,   8),
        (GLFW.BLUE_BITS,    8),

        (GLFW.STENCIL_BITS, 0),
        (GLFW.AUX_BUFFERS,  0)
    ]
end
"""
Function to create a pure GLFW OpenGL window
"""
function create_glcontext(
        name = "GLWindow";
        resolution = standard_screen_resolution(),
        debugging = false,
        major = 3,
        minor = 3,# this is what GLVisualize needs to offer all features
        windowhints = standard_window_hints(),
        contexthints = standard_context_hints(major, minor)
    )

    for ch in contexthints
        GLFW.WindowHint(ch...)
    end
    for wh in windowhints
        GLFW.WindowHint(wh...)
    end

    @static if is_apple()
        if debugging
            warn("OpenGL debug message callback not available on osx")
            debugging = false
        end
    end
    GLFW.WindowHint(GLFW.OPENGL_DEBUG_CONTEXT, Cint(debugging))

    window = GLFW.CreateWindow(resolution..., Compat.String(name))
    GLFW.MakeContextCurrent(window)

    debugging && glDebugMessageCallbackARB(_openglerrorcallback, C_NULL)
    window
end

"""
Most basic Screen constructor, which is usually used to create a parent screen.
It creates an OpenGL context and registeres all the callbacks
from the kw_arg `callbacks`.
You can change the OpenGL version with `major` and `minor`.
Also `windowhints` and `contexthints` can be given.
You can query the standard context and window hints
with `GLWindow.standard_context_hints` and `GLWindow.standard_window_hints`.
Finally you have the kw_args color and resolution. The first sets the background
color of the window and the other the resolution of the window.
"""
function Screen(name = "GLWindow";
        resolution = standard_screen_resolution(),
        debugging = false,
        major = 3,
        minor = 3,# this is what GLVisualize needs to offer all features
        windowhints = standard_window_hints(),
        contexthints = standard_context_hints(major, minor),
        callbacks = standard_callbacks(),
        color = RGBA{Float32}(1,1,1,1)
    )
    # create glcontext
    window = create_glcontext(
        name,
        resolution=resolution, debugging=debugging,
        major=major, minor=minor,
        windowhints=windowhints, contexthints=contexthints
    )
    GLFW.ShowWindow(window)

    #create standard signals
    signal_dict = register_callbacks(window, callbacks)
    @materialize window_position, window_size, hasfocus = signal_dict
    @materialize framebuffer_size, cursor_position = signal_dict
    window_area = map(SimpleRectangle,
        Signal(Vec(0,0)),
        framebuffer_size
    )
    signal_dict[:window_area] = window_area
    # seems to be necessary to set this as early as possible
    fb_size = value(framebuffer_size)
    glViewport(0, 0, fb_size...)

    # GLFW uses different coordinates from OpenGL, and on osx, the coordinates
    # are not in pixel coordinates
    # we coorect them to always be in pixel coordinates with 0,0 in left down corner
    signal_dict[:mouseposition] = const_lift(corrected_coordinates,
        Signal(window_size), Signal(framebuffer_size), cursor_position
    )
    signal_dict[:mouse2id] = Signal(SelectionID{Int}(-1, -1))

    # TODO: free when context is freed. We don't have a good abstraction of a gl context yet, though
    # (It could be shared, so it does not map directly to one screen)
    preserve(map(signal_dict[:window_open]) do open
        if !open
            GLAbstraction.empty_shader_cache!()
        end
        nothing
    end)
    GLFW.SwapInterval(0) # deactivating vsync seems to make everything quite a bit smoother
    screen = Screen(Symbol(name),
        window_area, Screen[], signal_dict,
        (), false, color,
        Dict{Symbol, Any}(),
        GLContext(window, GLFramebuffer(framebuffer_size))
    )
    signal_dict[:mouseinside] = droprepeats(
        const_lift(isinside, screen, signal_dict[:mouseposition])
    )
    screen
end

"""
Function that creates a screenshot from `window` and saves it to `path`.
You can choose the channel of the framebuffer, which is usually:
`color`, `depth` and `objectid`
"""
screenshot(window; path="screenshot.png", channel=:color) =
   save(path, screenbuffer(window, channel), true)

"""
Returns the contents of the framebuffer of `window` as a Julia Array.
You can choose the channel of the framebuffer, which is usually:
`color`, `depth` and `objectid`
"""
function screenbuffer(window, channel=:color)
    fb = framebuffer(window)
    channels = fieldnames(fb)[2:end]
    if channel in channels
        img = gpu_data(getfield(fb, channel))[window.area.value]
        return rotl90(img)
    end
    error("Channel $channel does not exist. Only these channels are available: $channels")
end

widths(s::Screen) = widths(value(s.area))
ishidden(s::Screen) = s.hidden
framebuffer(s::Screen) = s.glcontext.framebuffer
nativewindow(s::Screen) = s.glcontext.window

"""
Check if a Screen is opened.
"""
function Base.isopen(window::Screen)
    isopen(nativewindow(window))
end
function Base.isopen(window::GLFW.Window)
    window.handle == C_NULL && return false
    !GLFW.WindowShouldClose(window)
end
"""
Swap the framebuffers on the Screen.
"""
function swapbuffers(window::Screen)
    swapbuffers(nativewindow(window))
end
function swapbuffers(window::GLFW.Window)
    window.handle == C_NULL && return
    GLFW.SwapBuffers(window)
    return
end
"""
Poll events on the screen which will propogate signals through react.
"""
function pollevents()
    GLFW.PollEvents()
end

function mouse2id(s::Screen)
    s.inputs[:mouse2id]
end
export mouse2id
function mouseposition(s::Screen)
    s.inputs[:mouseposition]
end
"""
Empties the content of the renderlist
"""
function Base.empty!(s::Screen)
    s.renderlist = ()
    for c in s.children
        empty!(c)
    end
    empty!(s.children)
    nothing
end

"""
returns a copy of the renderlist
"""
function GLAbstraction.renderlist(s::Screen)
    vcat(s.renderlist...)
end
function destroy!(screen::Screen)
    empty!(screen)
    nw = nativewindow(screen)
    if nw.handle != C_NULL
        GLFW.DestroyWindow(nw)
        nw.handle = C_NULL
    end
end

get_id(x::Integer) = x
get_id(x::RenderObject) = x.id
function delete_robj!(list, robj)
    for (i, id) in enumerate(list)
        if get_id(id) == robj.id
            splice!(list, i)
            return true, i
        end
    end
    false, 0
end

function Base.delete!(screen::Screen, robj::RenderObject)
    for renderlist in screen.renderlist
        deleted, i = delete_robj!(renderlist, robj)
        deleted && return true
    end
    false
end

function GLAbstraction.robj_from_camera(window, camera)
    cam = window.cameras[camera]
    return filter(renderlist(window)) do robj
        robj[:projection] == cam.projection
    end
end


text:
@inline function abs_fpdb(x::T) where {T<:AbstractFloat}
    return DoubleFloat(abs_fp_dd(x))
end

@inline function negabs_fp_db(x::T) where {T<:AbstractFloat}
    return DoubleFloat(negabs_fp_dd(x))
end

@inline function neg_fp_db(x::T) where {T<:AbstractFloat}
    return DoubleFloat(neg_fp_dd(x))
end

@inline function inv_fp_db(x::T) where {T<:AbstractFloat}
    return DoubleFloat(inv_fp_dd(x))
end

@inline function powr2_fp_db(x::T) where {T<:AbstractFloat}
    return DoubleFloat(powr2_fp_dd(x))
end

@inline function powr3_fp_db(x::T) where {T<:AbstractFloat}
    return DoubleFloat(powr3_fp_dd(x))
end

@inline function powr4_fp_db(x::T) where {T<:AbstractFloat}
    return DoubleFloat(powr4_fp_dd(x))
end

@inline function powr5_fp_db(x::T) where {T<:AbstractFloat}
    return DoubleFloat(powr5_fp_dd(x))
end

@inline function powr6_fp_db(x::T) where {T<:AbstractFloat}
    return DoubleFloat(powr6_fp_dd(x))
end

@inline function root2_fp_db(x::T) where {T<:AbstractFloat}
    return DoubleFloat(root2_fp_dd(x))
end

@inline function root3_fp_db(x::T) where {T<:AbstractFloat}
    return DoubleFloat(root3_fp_dd(x))
end

@inline function root4_fp_db(x::T) where {T<:AbstractFloat}
    return DoubleFloat(root4_fp_dd(x))
end

@inline function root5_fp_db(x::T) where {T<:AbstractFloat}
    throw(ErrorException("not implemented"))
end

@inline function root6_fp_db(x::T) where {T<:AbstractFloat}
    return DoubleFloat(root6_fp_dd(x))
end


text:
import Base: isempty, length, haskey, peek, getindex, get, get!, setindex!, delete!, empty!, empty, merge!, iterate
using Base.Order: Ordering, lt

# This file contains code that was formerly a part of Julia. License is MIT: http://julialang.org/license

# Binary heap indexing
heapleft(i::Integer) = 2i
heapright(i::Integer) = 2i + 1
heapparent(i::Integer) = div(i, 2)

# PriorityQueue
# -------------

"""
    PriorityQueue{K, V}([ord])

Construct a new [`PriorityQueue`](@ref), with keys of type
`K` and values/priorites of type `V`.
If an order is not given, the priority queue is min-ordered using
the default comparison for `V`.

A `PriorityQueue` acts like a `Dict`, mapping values to their
priorities, with the addition of a `dequeue!` function to remove the
lowest priority element.

```jldoctest
julia> PriorityQueue(Base.Order.Forward, "a" => 2, "b" => 3, "c" => 1)
PriorityQueue{String,Int64,Base.Order.ForwardOrdering} with 3 entries:
  "c" => 1
  "a" => 2
  "b" => 3
```
"""
struct PriorityQueue{K,V,O<:Ordering} <: AbstractDict{K,V}
    # Binary heap of (element, priority) pairs.
    xs::Array{Pair{K,V}, 1}
    o::O

    # Map elements to their index in xs
    index::Dict{K, Int}

    function PriorityQueue{K,V,O}(o::O) where {K,V,O<:Ordering}
        new{K,V,O}(Vector{Pair{K,V}}(), o, Dict{K, Int}())
    end

    PriorityQueue{K, V, O}(xs::Array{Pair{K,V}, 1}, o::O, index::Dict{K, Int}) where {K,V,O<:Ordering} = new(xs, o, index)

    function PriorityQueue{K,V,O}(o::O, itr) where {K,V,O<:Ordering}
        xs = Vector{Pair{K,V}}(undef, length(itr))
        index = Dict{K, Int}()
        for (i, (k, v)) in enumerate(itr)
            xs[i] = Pair{K,V}(k, v)
            if haskey(index, k)
                throw(ArgumentError("PriorityQueue keys must be unique"))
            end
            index[k] = i
        end
        pq = new{K,V,O}(xs, o, index)

        # heapify
        for i in heapparent(length(pq.xs)):-1:1
            percolate_down!(pq, i)
        end

        return pq
    end
end

# A copy constructor
PriorityQueue(xs::Array{Pair{K,V}, 1}, o::O, index::Dict{K, Int}) where {K,V,O<:Ordering} =
    PriorityQueue{K,V,O}(xs, o, index)

# Any-Any constructors
PriorityQueue(o::Ordering=Forward) = PriorityQueue{Any,Any,typeof(o)}(o)

# Construction from Pairs
PriorityQueue(ps::Pair...) = PriorityQueue(Forward, ps)
PriorityQueue(o::Ordering, ps::Pair...) = PriorityQueue(o, ps)
PriorityQueue{K,V}(ps::Pair...) where {K,V} = PriorityQueue{K,V,ForwardOrdering}(Forward, ps)
PriorityQueue{K,V}(o::Ord, ps::Pair...) where {K,V,Ord<:Ordering} = PriorityQueue{K,V,Ord}(o, ps)

# Construction specifying Key/Value types
# e.g., PriorityQueue{Int,Float64}([1=>1, 2=>2.0])
PriorityQueue{K,V}(kv) where {K,V} = PriorityQueue{K,V}(Forward, kv)
function PriorityQueue{K,V}(o::Ord, kv) where {K,V,Ord<:Ordering}
    try
        PriorityQueue{K,V,Ord}(o, kv)
    catch e
        if not_iterator_of_pairs(kv)
            throw(ArgumentError("PriorityQueue(kv): kv needs to be an iterator of tuples or pairs"))
        else
            rethrow(e)
        end
    end
end

# Construction inferring Key/Value types from input
# e.g. PriorityQueue{}

PriorityQueue(o1::Ordering, o2::Ordering) = throw(ArgumentError("PriorityQueue with two parameters must be called with an Ordering and an interable of pairs"))
PriorityQueue(kv, o::Ordering=Forward) = PriorityQueue(o, kv)
function PriorityQueue(o::Ordering, kv)
    try
        _priority_queue_with_eltype(o, kv, eltype(kv))
    catch e
        if not_iterator_of_pairs(kv)
            throw(ArgumentError("PriorityQueue(kv): kv needs to be an iterator of tuples or pairs"))
        else
            rethrow(e)
        end
    end
end

_priority_queue_with_eltype(o::Ord, ps, ::Type{Pair{K,V}} ) where {K,V,Ord} = PriorityQueue{  K,  V,Ord}(o, ps)
_priority_queue_with_eltype(o::Ord, kv, ::Type{Tuple{K,V}}) where {K,V,Ord} = PriorityQueue{  K,  V,Ord}(o, kv)
_priority_queue_with_eltype(o::Ord, ps, ::Type{Pair{K}}   ) where {K,  Ord} = PriorityQueue{  K,Any,Ord}(o, ps)
_priority_queue_with_eltype(o::Ord, kv, ::Type            ) where {    Ord} = PriorityQueue{Any,Any,Ord}(o, kv)

## TODO: It seems impossible (or at least very challenging) to create the eltype below.
##       If deemed possible, please create a test and uncomment this definition.
# _priority_queue_with_eltype{  D,Ord}(o::Ord, ps, ::Type{Pair{K,V} where K}) = PriorityQueue{Any,  D,Ord}(o, ps)

length(pq::PriorityQueue) = length(pq.xs)
isempty(pq::PriorityQueue) = isempty(pq.xs)
haskey(pq::PriorityQueue, key) = haskey(pq.index, key)

"""
    peek(pq)

Return the lowest priority key from a priority queue without removing that
key from the queue.
"""
peek(pq::PriorityQueue) = pq.xs[1]

function percolate_down!(pq::PriorityQueue, i::Integer)
    x = pq.xs[i]
    @inbounds while (l = heapleft(i)) <= length(pq)
        r = heapright(i)
        j = r > length(pq) || lt(pq.o, pq.xs[l].second, pq.xs[r].second) ? l : r
        if lt(pq.o, pq.xs[j].second, x.second)
            pq.index[pq.xs[j].first] = i
            pq.xs[i] = pq.xs[j]
            i = j
        else
            break
        end
    end
    pq.index[x.first] = i
    pq.xs[i] = x
end

function percolate_up!(pq::PriorityQueue, i::Integer)
    x = pq.xs[i]
    @inbounds while i > 1
        j = heapparent(i)
        if lt(pq.o, x.second, pq.xs[j].second)
            pq.index[pq.xs[j].first] = i
            pq.xs[i] = pq.xs[j]
            i = j
        else
            break
        end
    end
    pq.index[x.first] = i
    pq.xs[i] = x
end

# Equivalent to percolate_up! with an element having lower priority than any other
function force_up!(pq::PriorityQueue, i::Integer)
    x = pq.xs[i]
    @inbounds while i > 1
        j = heapparent(i)
        pq.index[pq.xs[j].first] = i
        pq.xs[i] = pq.xs[j]
        i = j
    end
    pq.index[x.first] = i
    pq.xs[i] = x
end

getindex(pq::PriorityQueue, key) = pq.xs[pq.index[key]].second

function get(pq::PriorityQueue, key, default)
    i = get(pq.index, key, 0)
    i == 0 ? default : pq.xs[i].second
end

function get!(pq::PriorityQueue, key, default)
    i = get(pq.index, key, 0)
    if i == 0
        enqueue!(pq, key, default)
        return default
    else
        return pq.xs[i].second
    end
end

# Change the priority of an existing element, or equeue it if it isn't present.
function setindex!(pq::PriorityQueue{K, V}, value, key) where {K,V}
    if haskey(pq, key)
        i = pq.index[key]
        oldvalue = pq.xs[i].second
        pq.xs[i] = Pair{K,V}(key, value)
        if lt(pq.o, oldvalue, value)
            percolate_down!(pq, i)
        else
            percolate_up!(pq, i)
        end
    else
        enqueue!(pq, key, value)
    end
    return value
end

"""
    enqueue!(pq, k=>v)

Insert the a key `k` into a priority queue `pq` with priority `v`.

```jldoctest
julia> a = PriorityQueue(PriorityQueue("a"=>1, "b"=>2, "c"=>3))
PriorityQueue{String,Int64,Base.Order.ForwardOrdering} with 3 entries:
  "a" => 1
  "b" => 2
  "c" => 3

julia> enqueue!(a, "d"=>4)
PriorityQueue{String,Int64,Base.Order.ForwardOrdering} with 4 entries:
  "a" => 1
  "b" => 2
  "c" => 3
  "d" => 4
```
"""
function enqueue!(pq::PriorityQueue{K,V}, pair::Pair{K,V}) where {K,V}
    key = pair.first
    if haskey(pq, key)
        throw(ArgumentError("PriorityQueue keys must be unique"))
    end
    push!(pq.xs, pair)
    pq.index[key] = length(pq)
    percolate_up!(pq, length(pq))

    return pq
end

"""
enqueue!(pq, k, v)

Insert the a key `k` into a priority queue `pq` with priority `v`.

"""
enqueue!(pq::PriorityQueue, key, value) = enqueue!(pq, key=>value)
enqueue!(pq::PriorityQueue{K,V}, kv) where {K,V} = enqueue!(pq, Pair{K,V}(kv.first, kv.second))

"""
    dequeue!(pq)

Remove and return the lowest priority key from a priority queue.

```jldoctest
julia> a = PriorityQueue(Base.Order.Forward, ["a" => 2, "b" => 3, "c" => 1])
PriorityQueue{String,Int64,Base.Order.ForwardOrdering} with 3 entries:
  "c" => 1
  "a" => 2
  "b" => 3

julia> dequeue!(a)
"c"

julia> a
PriorityQueue{String,Int64,Base.Order.ForwardOrdering} with 2 entries:
  "a" => 2
  "b" => 3
```
"""
function dequeue!(pq::PriorityQueue)
    x = pq.xs[1]
    y = pop!(pq.xs)
    if !isempty(pq)
        pq.xs[1] = y
        pq.index[y.first] = 1
        percolate_down!(pq, 1)
    end
    delete!(pq.index, x.first)
    return x.first
end

function dequeue!(pq::PriorityQueue, key)
    idx = pq.index[key]
    force_up!(pq, idx)
    dequeue!(pq)
    return key
end

"""
    dequeue_pair!(pq)

Remove and return a the lowest priority key and value from a priority queue as a pair.

```jldoctest
julia> a = PriorityQueue(Base.Order.Forward, "a" => 2, "b" => 3, "c" => 1)
PriorityQueue{String,Int64,Base.Order.ForwardOrdering} with 3 entries:
  "c" => 1
  "a" => 2
  "b" => 3

julia> dequeue_pair!(a)
"c" => 1

julia> a
PriorityQueue{String,Int64,Base.Order.ForwardOrdering} with 2 entries:
  "a" => 2
  "b" => 3
```
"""
function dequeue_pair!(pq::PriorityQueue)
    x = pq.xs[1]
    y = pop!(pq.xs)
    if !isempty(pq)
        pq.xs[1] = y
        pq.index[y.first] = 1
        percolate_down!(pq, 1)
    end
    delete!(pq.index, x.first)
    return x
end

function dequeue_pair!(pq::PriorityQueue, key)
    idx = pq.index[key]
    force_up!(pq, idx)
    dequeue_pair!(pq)
end

"""
    delete!(pq, key)
Delete the mapping for the given key in a priority queue, and return the priority queue.
# Examples
```jldoctest
julia> q = PriorityQueue(Base.Order.Forward, "a"=>2, "b"=>3, "c"=>1)
PriorityQueue{String,Int64,Base.Order.ForwardOrdering} with 3 entries:
  "c" => 1
  "a" => 2
  "b" => 3
julia> delete!(q, "b")
PriorityQueue{String,Int64,Base.Order.ForwardOrdering} with 2 entries:
  "c" => 1
  "a" => 2
```
"""
function delete!(pq::PriorityQueue, key)
    dequeue_pair!(pq, key)
    return pq
end

function empty!(pq::PriorityQueue)
    empty!(pq.xs)
    empty!(pq.index)
    return pq
end

empty(pq::PriorityQueue) = PriorityQueue(empty(pq.xs), pq.o, empty(pq.index))

#merge!(d::SortedDict, other::PriorityQueue) = invoke(merge!, Tuple{AbstractDict, PriorityQueue}, d, other)

function merge!(d::AbstractDict, other::PriorityQueue)
    next = iterate(other, false)
    while next !== nothing
        (k, v), state = next
        d[k] = v
        next = iterate(other, state)
    end
    return d
end

function merge!(combine::Function, d::AbstractDict, other::PriorityQueue)
    next = iterate(other, false)
    while next !== nothing
        (k, v), state = next
        d[k] = haskey(d, k) ? combine(d[k], v) : v
        next = iterate(other, state)
    end
    return d
end

# Opaque not to be exported.
mutable struct _PQIteratorState{K, V, O <: Ordering}
    pq::PriorityQueue{K, V, O}
    _PQIteratorState{K, V, O}(pq::PriorityQueue{K, V, O}) where {K, V, O <: Ordering} = new(pq)
end

_PQIteratorState(pq::PriorityQueue{K, V, O}) where {K, V, O <: Ordering} = _PQIteratorState{K, V, O}(pq)

# Unordered iteration through key value pairs in a PriorityQueue
# O(n) iteration.
function _iterate(pq::PriorityQueue, state)
    (k, idx), i = state
    return (pq.xs[idx], i)
end
_iterate(pq::PriorityQueue, ::Nothing) = nothing

iterate(pq::PriorityQueue, ::Nothing) = nothing

function iterate(pq::PriorityQueue, ordered::Bool=true)
    if ordered
        isempty(pq) && return nothing
        state = _PQIteratorState(PriorityQueue(copy(pq.xs), pq.o, copy(pq.index)))
        return dequeue_pair!(state.pq), state
    else
        _iterate(pq, iterate(pq.index))
    end
end

function iterate(pq::PriorityQueue, state::_PQIteratorState)
    isempty(state.pq) && return nothing
    return dequeue_pair!(state.pq), state
end

iterate(pq::PriorityQueue, i) = _iterate(pq, iterate(pq.index, i))


text:
# ---
# title: 488. Zuma Game
# id: problem488
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Depth-first Search
# link: <https://leetcode.com/problems/zuma-game/description/>
# hidden: true
# ---
#
# Think about Zuma Game. You have a row of balls on the table, colored red(R),
# yellow(Y), blue(B), green(G), and white(W). You also have several balls in
# your hand.
#
# Each time, you may choose a ball in your hand, and insert it into the row
# (including the leftmost place and rightmost place). Then, if there is a group
# of 3 or more balls in the same color touching, remove these balls. Keep doing
# this until no more balls can be removed.
#
# Find the minimal balls you have to insert to remove all the balls on the
# table. If you cannot remove all the balls, output -1.
#
#
#
# **Example 1:**
#
#
#
#     Input: board = "WRRBBW", hand = "RB"
#     Output: -1
#     Explanation: WRRBBW -> WRR[R]BBW -> WBBW -> WBB[B]W -> WW
#
#
# **Example 2:**
#
#
#
#     Input: board = "WWRRBBWW", hand = "WRBRW"
#     Output: 2
#     Explanation: WWRRBBWW -> WWRR[R]BBWW -> WWBBWW -> WWBB[B]WW -> WWWW -> empty
#
#
# **Example 3:**
#
#
#
#     Input: board = "G", hand = "GGGGG"
#     Output: 2
#     Explanation: G -> G[G] -> GG[G] -> empty
#
#
# **Example 4:**
#
#
#
#     Input: board = "RBYYBBRRB", hand = "YRBGB"
#     Output: 3
#     Explanation: RBYYBBRRB -> RBYY[Y]BBRRB -> RBBBRRB -> RRRB -> B -> B[B] -> BB[B] -> empty
#
#
#
#
# **Constraints:**
#
#   * You may assume that the initial row of balls on the table won't have any 3 or more consecutive balls with the same color.
#   * `1 <= board.length <= 16`
#   * `1 <= hand.length <= 5`
#   * Both input strings will be non-empty and only contain characters 'R','Y','B','G','W'.
#
#
## @lc code=start
using LeetCode

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


text:
abstract type Overview end

Base.show(io::IO, mime::MIME"text/plain", ov::Overview) = _show(io, mime, ov)
Base.show(io::IO, mime::MIME"text/html", ov::Overview) = _show(io, mime, ov)

_show(io, mime, ov) = show(io, mime, Markdown.MD(ov))

function Markdown.MD(ov::Overview)
    io = IOBuffer()
    show(io, "text/markdown", ov)
    seek(io, 0)
    return Markdown.parse(io)
end

preprocess(m::Module; include = hasdoc, options...) = preprocess(
    DocumentationOverview.find(m; include = include);
    sortby = fullname,
    options...,
)

function preprocess(apis; sortby = nothing, signature = nothing, options...)
    apis = mapfoldl(
        signaturesetter(signature) ∘ DocumentationOverview.API,
        push!,
        apis;
        init = API[],
    )
    if sortby !== nothing
        sort!(apis; by = sortby)
    end
    return apis, options
end

function preprocess(
    fullnames::Expr;
    namespace = Main,
    sortby = nothing,
    signature = nothing,
    options...,
)
    if !Meta.isexpr(fullnames, :vect)
        error("expected an expression of form `:[a.b.c, d.e.f, ...]`; got: ", fullnames)
    end
    setter = signaturesetter(signature)
    apis = mapfoldl(push!, fullnames.args; init = API[]) do ex
        vars = fullname_symbols(ex)
        vars === nothing && error("unsupported: ", ex)
        vars = vars::Vector{Symbol}
        setter(resolvein(vars, namespace))
    end
    if sortby !== nothing
        sort!(apis; by = sortby)
    end
    return apis, options
end


text:
@everywhere using AccelerationBenchmark.OACCEL2017

tests = [(:A,100), (:A,200),
         (:B,100), (:B,200),
         (:C,100), (:C,200),
         (:D,500), (:D,1000),
         (:E,100), (:E,200),
         (:F,200), (:F,500),
         (:G,100), (:G,200),
         (:D, 50000), (:D, 100000),
         (:E, 50000), (:E, 100000)]

seeds = 0:999
OACCEL2017.runmany(tests,seeds)


text:
@testset "480.sliding-window-median.jl" begin
    @test median_sliding_window([1,3,-1,-3,5,3,6,7], 3) == [1,-1,-1,3,5,6.]
end


text:
module RandomExtensionsTests

using Test
include("tests.jl")

end


text:
import MaximumLikelihoodProblems

import GLM
import Statistics
import Test

β_true = [1.0, 2.0]

# @show β_true
# @show β_hat

Test.@test typeof(β_hat) == typeof(β_true)
Test.@test ndims(β_hat) == ndims(β_true)
Test.@test size(β_hat) == size(β_true)
absolute_error = abs.(β_hat - β_true)
square_error = abs2.(β_hat - β_true)
absolute_error_proportional = absolute_error ./ abs.(β_true)
square_error_proportional = square_error ./ abs.(β_true)
Test.@test sum(absolute_error) < 0.150
Test.@test sum(square_error) < 0.010
Test.@test sum(absolute_error_proportional) < 0.100
Test.@test sum(square_error_proportional) < 0.010
Test.@test maximum(absolute_error) < 0.100
Test.@test maximum(square_error) < 0.010
Test.@test maximum(absolute_error_proportional) < 0.100
Test.@test maximum(square_error_proportional) < 0.005
Test.@test Statistics.mean(absolute_error) < 0.070
Test.@test Statistics.mean(square_error) < 0.007
Test.@test Statistics.mean(absolute_error_proportional) < 0.050
Test.@test Statistics.mean(square_error_proportional) < 0.010

gradient_vector_at_θ_hat = MaximumLikelihoodProblems.Internal.gradient_vector(transformed_gradient_problem,
                                                                              θ_hat)
hessian_matrix_at_θ_hat = MaximumLikelihoodProblems.Internal.hessian_matrix(transformed_gradient_problem,
                                                                            θ_hat)
Test.@test MaximumLikelihoodProblems.Internal._is_approximately_zero(gradient_vector_at_θ_hat)
Test.@test MaximumLikelihoodProblems.Internal._is_approximately_hermitian(hessian_matrix_at_θ_hat)
Test.@test MaximumLikelihoodProblems.Internal._is_approximately_negative_definite(hessian_matrix_at_θ_hat;
                                                                                fuzz_factor = 10)
Test.@test MaximumLikelihoodProblems.Internal._is_approximately_positive_definite(-hessian_matrix_at_θ_hat;
                                                                                fuzz_factor = 10)

external_model_glm = GLM.glm(X, y, GLM.Binomial(), GLM.LogitLink())
beta_hat_external_model_glm = GLM.coef(external_model_glm)
Test.@test isapprox(β_hat, beta_hat_external_model_glm)


text:
export whitespace_char
whitespace_string = " \t\U0085\U200E\U200F\U2028\U2029"*"\U2029\U000C\U000B"
whitespace_char = CharIn("[:space:]",whitespace_string)

export horizontal_space_char
# The horizontal space characters are:
horizontal_space_char=CharIn(
    '\U0009', # "Horizontal tab (HT)"),
    '\U0020', # "Space"),
    '\U00A0', # "Non-break space"),
    '\U1680', # "Ogham space mark"),
    '\U180E', # "Mongolian vowel separator"),
    '\U2000', # "En quad"),
    '\U2001', # "Em quad"),
    '\U2002', # "En space"),
    '\U2003', # "Em space"),
    '\U2004', # "Three-per-em space"),
    '\U2005', # "Four-per-em space"),
    '\U2006', # "Six-per-em space"),
    '\U2007', # "Figure space"),
    '\U2008', # "Punctuation space"),
    '\U2009', # "Thin space"),
    '\U200A', # "Hair space"),
    '\U202F', # "Narrow no-break space"),
    '\U205F', # "Medium mathematical space"),
    '\U3000' # "Ideographic space"))
)

export vertical_space_char
# The vertical space characters are:
vertical_space_char=CharIn(
    '\U000A', # "Linefeed (LF)"),
    '\U000B', # "Vertical tab (VT)"),
    '\U000C', # "Form feed (FF)"),
    '\U000D', # "Carriage return (CR)"),
    '\U0085', # "Next line (NEL)"),
    '\U2028', # "Line separator"),
    '\U2029') # "Paragraph separator"))

export bsr
"""
```jldoctest
julia> CombinedParsers.Regexp.bsr
(?>|🗄...) Either |> Atomic
├─ \r\n
└─ [\n\x0b\f\r\x85] CharIn
::Union{Char, SubString}
```

PCRE backslash R (BSR), for newlines.
"""
bsr = Atomic(Either("\r\n",
                    !CharIn(raw"\n\x0b\f\r\x85", '\n','\x0b','\f','\r','\U0085', '\U2028','\U2029')));

export at_linestart, at_lineend
"""
```jldoctest
julia> CombinedParsers.Regexp.at_linestart
|🗄... Either
├─ ^ AtStart
└─ (?<=🗄...)) PositiveLookbehind
   ├─ \n\r
   └─ [\n\x0b\f\r\x85] CharIn
::Union{AtStart, Char, SubString}
```

Note: in PCRE
```julia
Either(
    on_options(Base.PCRE.MULTILINE,
           '^' => at_linestart),
    parser('^' => AtStart())
)
```
"""
@with_names at_linestart = Either(AtStart(),PositiveLookbehind(bsr))
@with_names at_lineend   = Either(AtEnd(),PositiveLookahead(bsr))

export newline, inline, whitespace_maybe, whitespace
newline = bsr
inline = !Repeat(CharNotIn(vertical_space_char))
whitespace_maybe = !Repeat(CharIn("\\h",horizontal_space_char))
whitespace_horizontal = !Repeat1(CharIn("\\h",horizontal_space_char))
whitespace = whitespace_horizontal

export whitespace_newline
whitespace_newline = Repeat1(CharIn("\\h\\v",horizontal_space_char,vertical_space_char))


text:
"""
    _setdiff(x::Interval{T}, y::Interval{T})

Computes the set difference x\\y and always returns a tuple of two intervals.
If the set difference is only one interval or is empty, then the returned tuple contains 1
or 2 empty intervals.
"""
function _setdiff(x::Interval{T}, y::Interval{T}) where T
    intersection = x ∩ y

    isempty(intersection) && return (x, emptyinterval(T))
    intersection == x && return (emptyinterval(T), emptyinterval(T))  # x is subset of y; setdiff is empty

    x.lo == intersection.lo && return (Interval(intersection.hi, x.hi), emptyinterval(T))
    x.hi == intersection.hi && return (Interval(x.lo, intersection.lo), emptyinterval(T))

    return (Interval(x.lo, y.lo), Interval(y.hi, x.hi))

end

"""
    setdiff(A::IntervalBox{N,T}, B::IntervalBox{N,T})

Returns a vector of `IntervalBox`es that are in the set difference `A ∖ B`,
i.e. the set of `x` that are in `A` but not in `B`.

Algorithm: Start from the total overlap (in all directions);
expand each direction in turn.
"""
function setdiff(A::IntervalBox{N,T}, B::IntervalBox{N,T}) where {N,T}

    intersection = A ∩ B
    isempty(intersection) && return [A]

    result_list = fill(IntervalBox(emptyinterval(T), N), 2*N)
    offset = 0
    x = A.v
    @inbounds for i = 1:N
        tmp = _setdiff(A[i], B[i])
        @inbounds for j = 1:2
            x = setindex(x, tmp[j], i)
            result_list[offset+j] = IntervalBox{N, T}(x)
        end
        offset += 2
        x = setindex(x, intersection[i], i)
    end
    filter!(!isempty, result_list)
end


text:
# assert file to test polynomial implementation
using Compat
using Base.Test
using Polynomials

pNULL = Poly(Float32[])
p0 = Poly([0])
p1 = Poly([1,0,0,0,0,0,0,0,0,0,0,0,0,0])
p2 = Poly([1,1,0,0])
p3 = Poly([1,2,1,0,0,0,0])
p4 = Poly([1,3,3,1,0,0])
p5 = Poly([1,4,6,4,1,0,0,0,0,0,0,0,0,0,0,0,0,0])
pN = Poly([276,3,87,15,24,0])
pR = Poly([3//4, -2//1, 1//1])
X = Poly([0.0, 1.0])
p1000 = Poly(randn(1000))

@test length(pNULL) == 1
@test length(p1000-p1000) == 1
@test length(p1000^0) == 1
@test length(0*p1000) == 1
@test length(p1000) == 1000
sprint(show, p1000)
sprint(show, pNULL)

@test p3 == Poly([1,2,1])
@test pN*10 == Poly([2760, 30, 870, 150, 240])
@test pN/10 == Poly([27.6, 0.3, 8.7, 1.5, 2.4])
@test 10*pNULL + pN == pN
@test 10*p0 + pN == pN
@test p5 + 2*p1 == Poly([3,4,6,4,1])
@test 10*pNULL - pN == -pN
@test p0 - pN == -pN
@test p5 - 2*p1 == Poly([-1,4,6,4,1])
@test p2*p2*p2 == p4
@test p2^4 == p5
@test pNULL^3 == pNULL
@test pNULL*pNULL == pNULL

@test map(degree, [pNULL,p0,p1,p2,p3,p4,p5,pN,pR,p1000]) == [0,0,0,1,2,3,4,4,2,999]

@test polyval(pN, -.125) == 276.9609375
@test polyval(pNULL, 10) == 0
@test polyval(p0, -10) == 0
@test isa(polyval(p0, X), Poly)
@test isa(polyval(p1, X), Poly)
@test polyval(poly([1//2, 3//2]), 1//2) == 0//1
@test polyder(polyint(pN)) == pN
@test polyder(pR) == Poly([-2//1,2//1])
@test polyint(pNULL,1) == p1
@test polyint(Poly(Rational[1,2,3])) == Poly(Rational[0, 1, 1, 1])
@test polyder(p3) == Poly([2,2])
@test polyder(p1) == polyder(p0) == polyder(pNULL) == pNULL

if VERSION >= v"0.4"
    @test pN(-.125) == 276.9609375
    @test pN([0.1, 0.2, 0.3]) == polyval(pN, [0.1, 0.2, 0.3])
end

@test poly([-1,-1]) == p3
@test roots(p0)==roots(p1)==roots(pNULL)==[]
@test roots(p2) == [-1]
a_roots = copy(pN.a)
@test all(abs(sort(roots(poly(a_roots))) - sort(a_roots)) .< 1e6)
@test length(roots(p5)) == 4
@test roots(pNULL) == []
@test roots(pR) == [1//2, 3//2]

@test pNULL + 2 == p0 + 2 == 2 + p0 == Poly([2])
@test p2 - 2 == -2 + p2 == Poly([-1,1])
@test 2 - p2 == Poly([1,-1])

p0 = Poly([0])
p1 = Poly([1])
p2 = Poly([5, 6, -3, 2 ,4])
p3 = Poly([7, -3, 2, 6])
p4 = p2 * p3
@test divrem(p4, p2) == (p3, zero(p3))
@test p3%p2 == p3
@test all((abs((p2 ÷ p3 - Poly([1/9,2/3])).a)) .< eps())
@test divrem(p0,p1) == (p0,p0)
@test divrem(p1,p1) == (p1,p0)
@test divrem(p2,p2) == (p1,p0)
@test divrem(pR, pR) == (one(pR), zero(pR))
@test_throws DivideError p1 ÷ p0
@test_throws DivideError divrem(p0,p0)

#Tests for multivariable support
pX = Poly([1, 2, 3, 4, 5])
pS1 = Poly([1, 2, 3, 4, 5], "s")
pS2 = Poly([1, 2, 3, 4, 5], 's')
pS3 = Poly([1, 2, 3, 4, 5], :s)
@test pX != pS1
@test pS1 == pS2
@test pS1 == pS3
@test_throws ErrorException pS1 + pX
@test_throws ErrorException pS1 - pX
@test_throws ErrorException pS1 * pX
@test_throws ErrorException pS1 ÷ pX
@test_throws ErrorException pS1 % pX

#Testing copying.
pcpy1 = Poly([1,2,3,4,5], :y)
pcpy2 = copy(pcpy1)
@test pcpy1 == pcpy2

#Tests for Pade approximants

println("Test for the exponential function.")
a = Poly(1.//convert(Vector{BigInt},gamma(BigFloat(1):BigFloat(17))),"x")
PQexp = Pade(a,8,8)
@test isapprox(convert(Float64, padeval(PQexp,1.0)), exp(1.0))
@test isapprox(convert(Float64, padeval(PQexp,-1.0)), exp(-1.0))

println("Test for the sine function.")
b = Poly(convert(Vector{BigInt},sinpi((0:16)/2)).//convert(Vector{BigInt},gamma(BigFloat(1):BigFloat(17))),"x")
PQsin = Pade(b,8,7)
@test isapprox(convert(Float64, padeval(PQsin,1.0)), sin(1.0))
@test isapprox(convert(Float64, padeval(PQsin,-1.0)),sin(-1.0))

println("Test for the cosine function.")
c = Poly(convert(Vector{BigInt},sinpi((1:17)/2)).//convert(Vector{BigInt},gamma(BigFloat(1):BigFloat(17))),"x")
PQcos = Pade(c,8,8)
@test isapprox(convert(Float64, padeval(PQcos,1.0)), cos(1.0))
@test isapprox(convert(Float64, padeval(PQcos,-1.0)), cos(-1.0))

println("Test for the summation of a factorially divergent series.")
d = Poly(convert(Vector{BigInt},(-1).^(0:60).*gamma(BigFloat(1):BigFloat(61.0))).//1,"x")
PQexpint = Pade(d,30,30)
@compat println("The approximate sum of the divergent series is:  ", Float64(padeval(PQexpint,1.0)))
println("The approximate sum of the convergent series is: ",exp(1)*(-γ-sum([(-1).^k/k./gamma(k+1) for k=1:20])))
@test isapprox(convert(Float64, padeval(PQexpint,1.0)),
               exp(1)*(-γ-sum([(-1).^k/k./gamma(k+1) for k=1:20])))

## polyfit
xs = linspace(0, pi, 10)
ys = sin(xs)
p = polyfit(xs, ys)
p = polyfit(xs, ys, :t)
p = polyfit(xs, ys, 2)
@test maximum(abs(map(x->polyval(p, x), xs) - ys)) <= 0.03

## truncation
p1 = Poly([1,1]/10)
p2 = Poly([1,2]/10)
p3 = Poly([1,3]/10)
psum = p1 + p2 - p3
@test degree(psum) == 1         # will have wrong degree
@test degree(truncate(psum)) == 0 # the degree should be correct after truncation

@test truncate(Poly([2,1]),reltol=1/2,abstol=0) == Poly([2])
@test truncate(Poly([2,1]),reltol=1,abstol=0)   == Poly([0])
@test truncate(Poly([2,1]),reltol=0,abstol=1)   == Poly([2])

## setindex!
println("Test for setindex!()")
p1    = Poly([1,2,1])
p1[5] = 1
@test p1[5] == 1
@test p1 == Poly([1,2,1,0,0,1])

## elementwise operations #52
println("Test for element-wise operations")
p1  = Poly([1, 2])
p2  = Poly([3, 1.])
p   = [p1, p2]
q   = [3, p1]
@test isa(q,Vector{Poly{Int64}})
psum  = p+3
pprod = p*3
pmin  = p-3
@test isa(psum, Vector{Poly{Float64}})
@test isa(pprod,Vector{Poly{Float64}})
@test isa(pmin, Vector{Poly{Float64}})

## getindex with ranges #43
p1 = Poly([4,5,6])
@test all(p1[0:1] .== [4,5])
@test all(p1[0:end] .== [4,5,6])
p1[0:1] = [7,8]
@test all(p1[0:end] .== [7,8,6])


text:
# Automatically generated using Clang.jl

const GEOS_VERSION_MAJOR = 3
const GEOS_VERSION_MINOR = 7
const GEOS_VERSION_PATCH = 2
const GEOS_VERSION = "3.7.2"
const GEOS_JTS_PORT = "1.13.0"
const GEOS_CAPI_VERSION_MAJOR = 1
const GEOS_CAPI_VERSION_MINOR = 11
const GEOS_CAPI_VERSION_PATCH = 2
const GEOS_CAPI_VERSION = "3.7.2-CAPI-1.11.2"
const GEOS_CAPI_FIRST_INTERFACE = GEOS_CAPI_VERSION_MAJOR
const GEOS_CAPI_LAST_INTERFACE = GEOS_CAPI_VERSION_MAJOR + GEOS_CAPI_VERSION_MINOR
const GEOS_PREC_NO_TOPO = 1 << 0
const GEOS_PREC_KEEP_COLLAPSED = 1 << 1
const GEOSContextHandle_HS = Cvoid
const GEOSContextHandle_t = Ptr{GEOSContextHandle_HS}
const GEOSMessageHandler = Ptr{Cvoid}
const GEOSMessageHandler_r = Ptr{Cvoid}
const GEOSGeom_t = Cvoid
const GEOSGeometry = GEOSGeom_t
const GEOSPrepGeom_t = Cvoid
const GEOSPreparedGeometry = GEOSPrepGeom_t
const GEOSCoordSeq_t = Cvoid
const GEOSCoordSequence = GEOSCoordSeq_t
const GEOSSTRtree_t = Cvoid
const GEOSSTRtree = GEOSSTRtree_t
const GEOSBufParams_t = Cvoid
const GEOSBufferParams = GEOSBufParams_t
const GEOSGeom = Ptr{GEOSGeometry}
const GEOSCoordSeq = Ptr{GEOSCoordSequence}

@cenum GEOSGeomTypes::UInt32 begin
    GEOS_POINT = 0
    GEOS_LINESTRING = 1
    GEOS_LINEARRING = 2
    GEOS_POLYGON = 3
    GEOS_MULTIPOINT = 4
    GEOS_MULTILINESTRING = 5
    GEOS_MULTIPOLYGON = 6
    GEOS_GEOMETRYCOLLECTION = 7
end

@cenum GEOSByteOrders::UInt32 begin
    GEOS_WKB_XDR = 0
    GEOS_WKB_NDR = 1
end

const GEOSQueryCallback = Ptr{Cvoid}
const GEOSDistanceCallback = Ptr{Cvoid}
const GEOSInterruptCallback = Cvoid

@cenum GEOSBufCapStyles::UInt32 begin
    GEOSBUF_CAP_ROUND = 1
    GEOSBUF_CAP_FLAT = 2
    GEOSBUF_CAP_SQUARE = 3
end

@cenum GEOSBufJoinStyles::UInt32 begin
    GEOSBUF_JOIN_ROUND = 1
    GEOSBUF_JOIN_MITRE = 2
    GEOSBUF_JOIN_BEVEL = 3
end

@cenum GEOSRelateBoundaryNodeRules::UInt32 begin
    GEOSRELATE_BNR_MOD2 = 1
    GEOSRELATE_BNR_OGC = 1
    GEOSRELATE_BNR_ENDPOINT = 2
    GEOSRELATE_BNR_MULTIVALENT_ENDPOINT = 3
    GEOSRELATE_BNR_MONOVALENT_ENDPOINT = 4
end

@cenum GEOSValidFlags::UInt32 begin
    GEOSVALID_ALLOW_SELFTOUCHING_RING_FORMING_HOLE = 1
end

const GEOSWKTReader_t = Cvoid
const GEOSWKTReader = GEOSWKTReader_t
const GEOSWKTWriter_t = Cvoid
const GEOSWKTWriter = GEOSWKTWriter_t
const GEOSWKBReader_t = Cvoid
const GEOSWKBReader = GEOSWKBReader_t
const GEOSWKBWriter_t = Cvoid
const GEOSWKBWriter = GEOSWKBWriter_t


text:
### A Pluto.jl notebook ###
# v0.14.5

using Markdown
using InteractiveUtils

# ╔═╡ 68fd4890-b438-11eb-3d17-a3cc7fb859ba
md"""
# [Day 1: No Time for a Taxicab](https://adventofcode.com/2016/day/1)

Santa's sleigh uses a very high-precision clock to guide its movements, and the clock's oscillator is regulated by stars. Unfortunately, the stars have been stolen... by the Easter Bunny. To save Christmas, Santa needs you to retrieve all fifty stars by December 25th.

Collect stars by solving puzzles. Two puzzles will be made available on each day in the Advent calendar; the second puzzle is unlocked when you complete the first. Each puzzle grants one star. Good luck!

You're airdropped near Easter Bunny Headquarters in a city somewhere. "Near", unfortunately, is as close as you can get - the instructions on the Easter Bunny Recruiting Document the Elves intercepted start here, and nobody had time to work them out further.

The Document indicates that you should start at the given coordinates (where you just landed) and face North. Then, follow the provided sequence: either turn left (`L`) or right (`R`) 90 degrees, then walk forward the given number of blocks, ending at a new intersection.

There's no time to follow such ridiculous instructions on foot, though, so you take a moment and work out the destination. Given that you can only walk on the [street grid of the city](https://en.wikipedia.org/wiki/Taxicab_geometry), how far is the shortest path to the destination?

For example:

-   Following `R2, L3` leaves you `2` blocks East and `3` blocks North, or `5` blocks away.
-   `R2, R2, R2` leaves you `2` blocks due South of your starting position, which is `2` blocks away.
-   `R5, L5, R5, R3` leaves you `12` blocks away.

**How many blocks away is Easter Bunny HQ?**
"""

# ╔═╡ 10b9f984-a80f-4bdc-8e91-182224860a6a
puzzle_input = "L4, R2, R4, L5, L3, L1, R4, R5, R1, R3, L3, L2, L2, R5, R1, L1, L2, R2, R2, L5, R5, R5, L2, R1, R2, L2, L4, L1, R5, R2, R1, R1, L2, L3, R2, L5, L186, L5, L3, R3, L5, R4, R2, L5, R1, R4, L1, L3, R3, R1, L1, R4, R2, L1, L4, R5, L1, R50, L4, R3, R78, R4, R2, L4, R3, L4, R4, L1, R5, L4, R1, L2, R3, L2, R5, R5, L4, L1, L2, R185, L5, R2, R1, L3, R4, L5, R2, R4, L3, R4, L2, L5, R1, R2, L2, L1, L2, R2, L2, R1, L5, L3, L4, L3, L4, L2, L5, L5, R2, L3, L4, R4, R4, R5, L4, L2, R4, L5, R3, R1, L1, R3, L2, R2, R1, R5, L4, R5, L3, R2, R3, R1, R4, L4, R1, R3, L5, L1, L3, R2, R1, R4, L4, R3, L3, R3, R2, L3, L3, R4, L2, R4, L3, L4, R5, R1, L1, R5, R3, R1, R3, R4, L1, R4, R3, R1, L5, L5, L4, R4, R3, L2, R1, R5, L3, R4, R5, L4, L5, R2"

# ╔═╡ a5bc4a35-a5cb-4411-805f-49edb1ac09c4

# ╔═╡ 67881f6b-641c-4a8b-9536-05cad11fffe7
function exec_moves(input)
	moves = split(input, ", ")
	pos = [0,0]
	direction = 1 # north

	for move in moves
		distance = parse(Int, move[2])
		if move[1] == "L"
			direction -= 1
			if direction < 1
				direction = 4
			end
		elseif move[1] == "R"
			direction += 1
			if direction > 4
				direction = 1
			end
		end

		if direction == 1 # north
			pos[1] += distance
		elseif direction == 2 # east
			pos[2] += distance
		elseif direction == 3 # south
			pos[1] -= distance
		elseif direction == 4 # west
			pos[2] -= distance
		end
	end

	pos
end

# ╔═╡ ddf66bd4-e128-4c19-a296-1aa80210df61
@assert exec_moves("R2, L3") == (3,2)

# ╔═╡ 414b8bb2-d348-4f3f-aea2-dc09a0106018
exec_moves("R2, L3")

# ╔═╡ Cell order:
# ╟─68fd4890-b438-11eb-3d17-a3cc7fb859ba
# ╟─10b9f984-a80f-4bdc-8e91-182224860a6a
# ╠═a5bc4a35-a5cb-4411-805f-49edb1ac09c4
# ╠═67881f6b-641c-4a8b-9536-05cad11fffe7
# ╠═ddf66bd4-e128-4c19-a296-1aa80210df61
# ╠═414b8bb2-d348-4f3f-aea2-dc09a0106018


text:
function gp_formula(formula::AbstractString)
    spl1 = findfirst(==(':'), formula)
    spl2 = last(findfirst(==('~'), formula))
    isnothing(spl2) && throw(ArgumentError("Invalid formula specification"))
    barind = nextind(formula, spl2)
    formula[barind] == '|' || throw(ArgumentError("Invalid formula specification"))

    if isnothing(spl1) || spl1 > spl2
        lik = GaPLAC.Gaussian()
        spl1 = prevind(formula, spl2)
    else
        spl1 < spl2 || throw(ArgumentError("Invalid formula specification"))
        lik = strip(formula[nextind(formula, spl1):prevind(formula, spl2)])
        lik = isempty(lik) ? GaPLAC.Gaussian() : GaPLAC.eval(Meta.parse(lik))
    end

    resp = strip(formula[1:prevind(formula, spl1)])

    gp = strip(formula[nextind(formula, barind):end])
    gp = GaPLAC.eval(Meta.parse(gp))

    return GPSpec(Symbol(resp), lik, gp)
end

max_stars_repo_name: biobakery/GaPLAC

max_issues_repo_name: biobakery/GaPLAC

max_forks_repo_name: biobakery/GaPLAC


text:
using PkgTemplateXKDR
using Test

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


text:
module TestUtilities

using Test
using CategoricalDistributions
using CategoricalArrays
using Random
using StableRNGs

rng = StableRNGs.StableRNG(123)

import CategoricalDistributions: classes, transform, decoder, int

@testset "classes" begin
    v = categorical(collect("asqfasqffqsaaaa"), ordered=true)
    @test classes(v[1]) == levels(v)
    @test classes(v) == levels(v)
    levels!(v, reverse(levels(v)))
    @test classes(v[1]) == levels(v)
    @test classes(v) == levels(v)
    vsub = view(v, 1:2)
    @test classes(vsub) == classes(v)
end

@testset "int, classes, decoder" begin
    N = 10
    mix = shuffle(rng, 0:N - 1)

    Xraw = broadcast(x->mod(x,N), rand(rng, Int, 2N, 3N))
    Yraw = string.(Xraw)

    # to turn a categ matrix into a ordinary array with categorical
    # elements. Needed because broacasting the identity gives a
    # categorical array in CategoricalArrays >0.5.2
    function matrix_(X)
        ret = Array{Any}(undef, size(X))
        for i in eachindex(X)
            ret[i] = X[i]
        end
        return ret
    end

    X = categorical(Xraw)
    x = X[1]
    Y = categorical(Yraw)
    y = Y[1]
    V = matrix_(X)
    W = matrix_(Y)

    # broadcasted encoding:
    @test int(X) == int(V)
    @test int(Y) == int(W)

    # encoding is right-inverse to decoding:
    d = decoder(x)
    @test d(int(V)) == V # ie have the same elements
    e = decoder(y)
    @test e(int(W)) == W

    @test int(classes(y)) == 1:length(classes(x))

    # int is based on ordering not index
    v = categorical(['a', 'b', 'c'], ordered=true)
    @test int(v) == 1:3
    levels!(v, ['c', 'a', 'b'])
    @test int(v) == [2, 3, 1]
    @test int(missing) |> ismissing
    v_int = int(categorical(['a', missing]))
    @test length(v_int) == 2
    @test v_int[1] == 1
    @test ismissing(v_int[2])

    # Errors
    @test_throws DomainError int("g")
end

@testset "transforming from raw values and categorical values" begin
    values = vcat([missing, ], collect("asdfjklqwerpoi"))
    Xraw = rand(rng, values, 15, 10)
    X = categorical(Xraw)
    element = skipmissing(X) |> first

    @test transform(element, missing) |> ismissing

    # errors:
    @test_throws(CategoricalDistributions.err_missing_class("non-existent"),
                 transform(element, "non-existent"))

    raw = first(skipmissing(Xraw))
    c = transform(element, raw)
    @test Set(classes(c)) == Set(classes(X))
    @test c == first(skipmissing(X))

    RAW = Xraw[2:end-1,2:end-1]
    C = transform(element, RAW)
    @test Set(classes(C)) == Set(classes(X))
    @test identity.(skipmissing(C)) ==
        identity.(skipmissing(X[2:end-1,2:end-1]))

    raw = first(skipmissing(Xraw))
    c = transform(X, raw)
    @test Set(classes(c)) == Set(classes(X))
    @test c == first(skipmissing(X))

    RAW = Xraw[2:end-1,2:end-1]
    C = transform(X, RAW)
    @test Set(classes(C)) == Set(classes(X))
    @test identity.(skipmissing(C)) ==
        identity.(skipmissing(X[2:end-1,2:end-1]))

    pool = CategoricalArrays.pool(element)
    a = transform(element, 'a')
    @test a == a
    @test a isa CategoricalValue
end

end

true


text:
using Base.Test
using toJuMP
using JuMP
using Ipopt

testdir = Pkg.dir("toJuMP")

include("$(testdir)/test/gms.jl")
include("$(testdir)/test/mod.jl")


text:
using Oceananigans: AbstractModel, AbstractOutputWriter, AbstractDiagnostic

using Oceananigans.Architectures: AbstractArchitecture, CPU
using Oceananigans.AbstractOperations: @at
using Oceananigans.Distributed
using Oceananigans.Advection: CenteredSecondOrder
using Oceananigans.BoundaryConditions: regularize_field_boundary_conditions
using Oceananigans.Fields: Field, tracernames, TracerFields, XFaceField, YFaceField, CenterField, compute!
using Oceananigans.Forcings: model_forcing
using Oceananigans.Grids: with_halo, topology, inflate_halo_size, halo_size, Flat, architecture, RectilinearGrid, Face, Center
using Oceananigans.ImmersedBoundaries: ImmersedBoundaryGrid
using Oceananigans.TimeSteppers: Clock, TimeStepper, update_state!
using Oceananigans.TurbulenceClosures: with_tracers, DiffusivityFields
using Oceananigans.Utils: tupleit
using Oceananigans.Models.HydrostaticFreeSurfaceModels: validate_tracer_advection
import Oceananigans.Architectures: architecture

const RectilinearGrids =  Union{RectilinearGrid, ImmersedBoundaryGrid{<:Any, <:Any, <:Any, <:Any, <:RectilinearGrid}}

function ShallowWaterTendencyFields(grid, tracer_names, prognostic_names)
    u =  XFaceField(grid)
    v =  YFaceField(grid)
    h = CenterField(grid)
    tracers = TracerFields(tracer_names, grid)

    return NamedTuple{prognostic_names}((u, v, h, Tuple(tracers)...))
end

function ShallowWaterSolutionFields(grid, bcs, prognostic_names)
    u =  XFaceField(grid, boundary_conditions = getproperty(bcs, prognostic_names[1]))
    v =  YFaceField(grid, boundary_conditions = getproperty(bcs, prognostic_names[2]))
    h = CenterField(grid, boundary_conditions = getproperty(bcs, prognostic_names[3]))

    return NamedTuple{prognostic_names[1:3]}((u, v, h))
end

mutable struct ShallowWaterModel{G, A<:AbstractArchitecture, T, V, R, F, E, B, Q, C, K, TS, FR} <: AbstractModel{TS}

                          grid :: G         # Grid of physical points on which `Model` is solved
                  architecture :: A         # Computer `Architecture` on which `Model` is run
                         clock :: Clock{T}  # Tracks iteration number and simulation time of `Model`
    gravitational_acceleration :: T         # Gravitational acceleration, full, or reduced
                     advection :: V         # Advection scheme for velocities, mass and tracers
                      coriolis :: R         # Set of parameters for the background rotation rate of `Model`
                       forcing :: F         # Container for forcing functions defined by the user
                       closure :: E         # Diffusive 'turbulence closure' for all model fields
                    bathymetry :: B         # Bathymetry/Topography for the model
                      solution :: Q         # Container for transports `uh`, `vh`, and height `h`
                       tracers :: C         # Container for tracer fields
            diffusivity_fields :: K         # Container for turbulent diffusivities
                   timestepper :: TS        # Object containing timestepper fields and parameters
                   formulation :: FR        # Either conservative or vector-invariant
end

struct ConservativeFormulation end

struct VectorInvariantFormulation end

"""
    ShallowWaterModel(; grid,
                        gravitational_acceleration,
                              clock = Clock{eltype(grid)}(0, 0, 1),
                 momentum_advection = UpwindBiasedFifthOrder(),
                   tracer_advection = WENO5(),
                     mass_advection = WENO5(),
                           coriolis = nothing,
                forcing::NamedTuple = NamedTuple(),
                            closure = nothing,
                         bathymetry = nothing,
                            tracers = (),
                 diffusivity_fields = nothing,
    boundary_conditions::NamedTuple = NamedTuple(),
                timestepper::Symbol = :RungeKutta3,
                        formulation = ConservativeFormulation())

Construct a shallow water model on `grid` with `gravitational_acceleration` constant.

Keyword arguments
=================

  - `grid`: (required) The resolution and discrete geometry on which `model` is solved. The
            architecture (CPU/GPU) that the model is solve is inferred from the architecture
            of the grid.
  - `gravitational_acceleration`: (required) The gravitational acceleration constant.
  - `clock`: The `clock` for the model.
  - `momentum_advection`: The scheme that advects velocities. See `Oceananigans.Advection`.
    Default: `UpwindBiasedFifthOrder()`.
  - `tracer_advection`: The scheme that advects tracers. See `Oceananigans.Advection`. Default: `WENO5()`.
  - `mass_advection`: The scheme that advects the mass equation. See `Oceananigans.Advection`. Default:
    `WENO5()`.
  - `coriolis`: Parameters for the background rotation rate of the model.
  - `forcing`: `NamedTuple` of user-defined forcing functions that contribute to solution tendencies.
  - `closure`: The turbulence closure for `model`. See `Oceananigans.TurbulenceClosures`.
  - `bathymetry`: The bottom bathymetry.
  - `tracers`: A tuple of symbols defining the names of the modeled tracers, or a `NamedTuple` of
               preallocated `CenterField`s.
  - `diffusivity_fields`: Stores diffusivity fields when the closures require a diffusivity to be
                          calculated at each timestep.
  - `boundary_conditions`: `NamedTuple` containing field boundary conditions.
  - `timestepper`: A symbol that specifies the time-stepping method. Either `:QuasiAdamsBashforth2` or
                   `:RungeKutta3` (default).
  - `formulation`: Whether the dynamics are expressed in conservative form (`ConservativeFormulation()`;
                   default) or in non-conservative form with a vector-invariant formulation for the
                   non-linear terms (`VectorInvariantFormulation()`).

!!! warning "Formulation-grid compatibility requirements"
    The `ConservativeFormulation()` requires `RectilinearGrid`.
    Use `VectorInvariantFormulation()` with `LatitudeLongitudeGrid`.
"""
function ShallowWaterModel(;
                           grid,
                           gravitational_acceleration,
                               clock = Clock{eltype(grid)}(0, 0, 1),
                  momentum_advection = UpwindBiasedFifthOrder(),
                    tracer_advection = WENO5(),
                      mass_advection = WENO5(),
                            coriolis = nothing,
                 forcing::NamedTuple = NamedTuple(),
                             closure = nothing,
                          bathymetry = nothing,
                             tracers = (),
                  diffusivity_fields = nothing,
     boundary_conditions::NamedTuple = NamedTuple(),
                 timestepper::Symbol = :RungeKutta3,
                         formulation = ConservativeFormulation())

    arch = architecture(grid)

    tracers = tupleit(tracers) # supports tracers=:c keyword argument (for example)

    topology(grid, 3) === Flat ||
        throw(ArgumentError("ShallowWaterModel requires `topology(grid, 3) === Flat`. " *
                            "Use `topology = ($(topology(grid, 1)), $(topology(grid, 2)), Flat)` " *
                            "when constructing `grid`."))

    (typeof(grid) <: RectilinearGrids || formulation == VectorInvariantFormulation()) ||
        throw(ArgumentError("`ConservativeFormulation()` requires a rectilinear `grid`. \n" *
                            "Use `VectorInvariantFormulation()` or change your grid to a rectilinear one."))

    Hx, Hy, Hz = inflate_halo_size(grid.Hx, grid.Hy, 0, topology(grid), momentum_advection, tracer_advection, mass_advection, closure)
    any((grid.Hx, grid.Hy, grid.Hz) .< (Hx, Hy, 0)) && # halos are too small, remake grid
        (grid = with_halo((Hx, Hy, 0), grid))

    prognostic_field_names = formulation isa ConservativeFormulation ? (:uh, :vh, :h, tracers...) :  (:u, :v, :h, tracers...)
    default_boundary_conditions = NamedTuple{prognostic_field_names}(Tuple(FieldBoundaryConditions()
                                                                           for name in prognostic_field_names))

    momentum_advection = validate_momentum_advection(momentum_advection, formulation)

    if isnothing(tracer_advection)
        tracer_advection_tuple = NamedTuple{tracernames(tracers)}(nothing for tracer in 1:length(tracers))
    else
        default_tracer_advection, tracer_advection = validate_tracer_advection(tracer_advection, grid)

        # Advection schemes
        tracer_advection_tuple = with_tracers(tracernames(tracers),
                                            tracer_advection,
                                            (name, tracer_advection) -> default_tracer_advection,
                                            with_velocities=false)
    end

    advection = merge((momentum=momentum_advection, mass=mass_advection), tracer_advection_tuple)

    bathymetry_field = CenterField(grid)
    if !isnothing(bathymetry)
        set!(bathymetry_field, bathymetry)
        fill_halo_regions!(bathymetry_field)
    else
        fill!(bathymetry_field, 0.0)
    end

    boundary_conditions = merge(default_boundary_conditions, boundary_conditions)
    boundary_conditions = regularize_field_boundary_conditions(boundary_conditions, grid, prognostic_field_names)

    solution           = ShallowWaterSolutionFields(grid, boundary_conditions, prognostic_field_names)
    tracers            = TracerFields(tracers, grid, boundary_conditions)
    diffusivity_fields = DiffusivityFields(diffusivity_fields, grid, tracernames(tracers), boundary_conditions, closure)

    # Instantiate timestepper if not already instantiated
    timestepper = TimeStepper(timestepper, grid, tracernames(tracers);
                              Gⁿ = ShallowWaterTendencyFields(grid, tracernames(tracers), prognostic_field_names),
                              G⁻ = ShallowWaterTendencyFields(grid, tracernames(tracers), prognostic_field_names))

    # Regularize forcing and closure for model tracer and velocity fields.
    model_fields = merge(solution, tracers)
    forcing = model_forcing(model_fields; forcing...)
    closure = with_tracers(tracernames(tracers), closure)

    model = ShallowWaterModel(grid,
                              arch,
                              clock,
                              eltype(grid)(gravitational_acceleration),
                              advection,
                              coriolis,
                              forcing,
                              closure,
                              bathymetry_field,
                              solution,
                              tracers,
                              diffusivity_fields,
                              timestepper,
                              formulation)

    update_state!(model)

    return model
end

using Oceananigans.Advection: VectorInvariantSchemes

validate_momentum_advection(momentum_advection, formulation) = momentum_advection
validate_momentum_advection(momentum_advection, ::VectorInvariantFormulation) =
    throw(ArgumentError("VectorInvariantFormulation requires a vector invariant momentum advection scheme. \n"*
                        "Use `momentum_advection = VectorInvariant()`."))
validate_momentum_advection(momentum_advection::VectorInvariantSchemes, ::VectorInvariantFormulation) = momentum_advection

formulation(model::ShallowWaterModel)  = model.formulation
architecture(model::ShallowWaterModel) = model.architecture

function shallow_water_velocities(model::ShallowWaterModel)
    if formulation(model) isa VectorInvariantFormulation
        return (model.solution.u, model.solution.v)
    else
        u = Field(@at (Face, Center, Center) model.solution.uh / model.solution.h)
        v = Field(@at (Center, Face, Center) model.solution.vh / model.solution.h)

        compute!(u)
        compute!(v)

        return (u, v)
    end
end


text:
"""
    reverse_complement(s::AbstractString)
Given: A DNA string s of length at most 1000 bp.

Return: The reverse complement s^c of s.
"""

function reverse_complement(s::AbstractString)
    rules = Dict('A' => 'T', 'T' => 'A', 'G' => 'C', 'C' => 'G')
    return map(x -> rules[x], reverse(s))
end

max_issues_repo_name: AugustoCL/Julia

max_forks_repo_name: AugustoCL/Julia


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#   Compute the fundamental arguments related to the IAU-2006 theory.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# References
# ==============================================================================
#
#   [1] Vallado, D. A (2013). Fundamentals of Astrodynamics and Applications.
#       Microcosm Press, Hawthorn, CA, USA.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

export luni_solar_args_iau2006, planetary_args_iau2006

"""
    luni_solar_args_iau2006(jd_tt::Number)

Compute the fundamental arguments related to the luni-solar effect for the
IAU-2006 theory **[1]**(p. 211).

The returned values are in [rad].

# References

- **[1]**: Vallado, D. A (2013). Fundamentals of Astrodynamics and Applications.
    Microcosm Press, Hawthorn, CA, USA.
"""
function luni_solar_args_iau2006(jd_tt::Number)
    # Compute the Julian Centuries from `jd_tt`.
    t_tt = (jd_tt - JD_J2000) / 36525

    # Auxiliary variables
    # ===================

    a2d = 1 / 3600
    d2r = π / 180

    # Delaunay arguments of the Sun and Moon
    # ======================================
    #
    # Evaluate the Delaunay arguments associated with the Moon and the Sun in
    # [arcsec] [1, p. 210].

    M_m  = @evalpoly(
        t_tt,
        +485868.249036,
        +1717915923.2178,
        +31.8792,
        +0.051635,
        -0.00024470
    )

    M_s  = @evalpoly(
        t_tt,
        +1287104.79305,
        +129596581.0481,
        -0.5532,
        +0.000136,
        -0.00001149
    )

    u_Mm = @evalpoly(
        t_tt,
        +335779.526232,
        +1739527262.8478,
        -12.7512,
        -0.001037,
        +0.00000417
    )

    D_s  = @evalpoly(
        t_tt,
        +1072260.70369,
        +1602961601.2090,
        -6.3706,
        +0.006593,
        -0.00003169
    )

    Ω_m  = @evalpoly(
        t_tt,
        +450160.398036,
        -6962890.5431,
        +7.4722,
        +0.007702,
        -0.00005939
    )

    # Convert to the interval [0,2π].
    M_s  = mod( M_s * a2d, 360) * d2r
    M_m  = mod( M_m * a2d, 360) * d2r
    u_Mm = mod(u_Mm * a2d, 360) * d2r
    D_s  = mod( D_s * a2d, 360) * d2r
    Ω_m  = mod( Ω_m * a2d, 360) * d2r

    return M_s, M_m, u_Mm, D_s, Ω_m
end

"""
    planetary_args_iau2006(jd_tt::Number)

Compute the fundamental arguments related to the planetary effects for the
IAU-2006 theory **[1]**(p. 211).

The returned values are in [rad].

# References

- **[1]**: Vallado, D. A (2013). Fundamentals of Astrodynamics and Applications.
    Microcosm Press, Hawthorn, CA, USA.
"""
function planetary_args_iau2006(jd_tt::Number)
    # Compute the Julian Centuries from `jd_tt`.
    t_tt = (jd_tt - JD_J2000) / 36525

    # Mean Heliocentric longitudes of the planets.
    #
    # TODO: In the example in [1, p. 221], the value related to Uranus is slight
    # different. The value used here is shown in [1](p. 211).

    λ_M☿ = @evalpoly(t_tt, 4.402_608_842, 2608.790_314_157_4)
    λ_M♀ = @evalpoly(t_tt, 3.176_146_697, 1021.328_554_621_1)
    λ_Me = @evalpoly(t_tt, 1.753_470_314,  628.307_584_999_1)
    λ_M♂ = @evalpoly(t_tt, 6.203_480_913,  334.061_242_670_0)
    λ_M♃ = @evalpoly(t_tt, 0.599_546_497,   52.969_096_264_1)
    λ_M♄ = @evalpoly(t_tt, 0.874_016_757,   21.329_910_496_0)
    λ_M⛢ = @evalpoly(t_tt, 5.481_293_872,    7.478_159_856_7)
    λ_M♆ = @evalpoly(t_tt, 5.311_886_287,    3.813_303_563_8)

    # General precession in longitude.
    p_λ = @evalpoly(t_tt, 0, 0.024_381_75, 0.000_005_386_91)

    # Convert to the interval [0,2π].
    λ_M☿ = mod(λ_M☿, 2π)
    λ_M♀ = mod(λ_M♀, 2π)
    λ_Me = mod(λ_Me, 2π)
    λ_M♂ = mod(λ_M♂, 2π)
    λ_M♃ = mod(λ_M♃, 2π)
    λ_M♄ = mod(λ_M♄, 2π)
    λ_M⛢ = mod(λ_M⛢, 2π)
    λ_M♆ = mod(λ_M♆, 2π)
    p_λ  = mod(p_λ,  2π)

    return λ_M☿, λ_M♀, λ_Me, λ_M♂, λ_M♃, λ_M♄, λ_M⛢, λ_M♆, p_λ
end


text:
# Tests for invertible neural network module
# Author: Philipp Witte, pwitte3@gatech.edu
# Date: January 2020
using InvertibleNetworks, Test

const test_suite = get(ENV, "test_suite", "all") # "all", "basics", "layers" or "networks"

basics = ["test_utils/test_objectives.jl",
          "test_utils/test_sequential.jl",
          "test_utils/test_nnlib_convolution.jl",
          "test_utils/test_activations.jl",
          "test_utils/test_squeeze.jl",
          "test_utils/test_jacobian.jl",
          "test_utils/test_chainrules.jl"]

          # Layers
layers = ["test_layers/test_residual_block.jl",
          "test_layers/test_flux_block.jl",
          "test_layers/test_resnet.jl",
          "test_layers/test_layer_conv1x1.jl",
          "test_layers/test_coupling_layer_basic.jl",
          "test_layers/test_coupling_layer_irim.jl",
          "test_layers/test_coupling_layer_glow.jl",
          "test_layers/test_coupling_layer_hint.jl",
          "test_layers/test_conditional_layer_hint.jl",
          "test_layers/test_conditional_res_block.jl",
          "test_layers/test_hyperbolic_layer.jl",
          "test_layers/test_actnorm.jl",
          "test_layers/test_layer_affine.jl"]

networks = ["test_networks/test_unrolled_loop.jl",
            "test_networks/test_generator.jl",
            "test_networks/test_glow.jl",
            "test_networks/test_hyperbolic_network.jl",
            "test_networks/test_multiscale_hint_network.jl",
            "test_networks/test_conditional_hint_network.jl"]

if test_suite == "all" || test_suite == "basics"
    @testset "Basics" begin
        for t=basics
            @testset "Test $t" begin
                @time include(t)
            end
        end
    end
end

if test_suite == "all" || test_suite == "layers"
    @testset "Layers" begin
        for t=layers
            @testset  "Test $t" begin
                @time include(t)
            end
        end
    end
end

# Networks
if test_suite == "all" || test_suite == "networks"
    @testset "Networks" begin
        for t=networks
            @testset  "Test $t" begin
                @time include(t)
            end
        end
    end
end