text:
__precompile__()
module CSTParser
global debug = true

using Tokenize
import Base: next, start, done, length, first, last, endof, getindex, setindex!
import Tokenize.Tokens
import Tokenize.Tokens: Token, iskeyword, isliteral, isoperator
import Tokenize.Lexers: Lexer, peekchar, iswhitespace

export ParseState, parse_expression

include("hints.jl")
import .Diagnostics: Diagnostic, LintCodes, FormatCodes

include("lexer.jl")
include("errors.jl")
include("spec.jl")
include("utils.jl")
include("components/array.jl")
include("components/curly.jl")
include("components/operators.jl")
include("components/do.jl")
include("components/functions.jl")
include("components/genericblocks.jl")
include("components/if.jl")
include("components/let.jl")
include("components/loops.jl")
include("components/macros.jl")
include("components/modules.jl")
include("components/prefixkw.jl")
include("components/quote.jl")
include("components/strings.jl")
include("components/try.jl")
include("components/types.jl")
include("components/tuples.jl")
include("conversion.jl")
include("display.jl")
include("scoping.jl")

"""
    parse_expression(ps)

Parses an expression until `closer(ps) == true`. Expects to enter the
`ParseState` the token before the the beginning of the expression and ends
on the last token.

Acceptable starting tokens are:
+ A keyword
+ An opening parentheses or brace.
+ An operator.
+ An instance (e.g. identifier, number, etc.)
+ An `@`.

"""
function parse_expression(ps::ParseState)
    startbyte = ps.nt.startbyte
    next(ps)
    if Tokens.begin_keywords < ps.t.kind < Tokens.end_keywords && ps.t.kind != Tokens.DO
        @catcherror ps startbyte ret = parse_kw(ps, Val{ps.t.kind})
    elseif ps.t.kind == Tokens.LPAREN
        @catcherror ps startbyte ret = parse_paren(ps)
    elseif ps.t.kind == Tokens.LSQUARE
        @catcherror ps startbyte ret = parse_array(ps)
    elseif ps.t.kind == Tokens.LBRACE
        @catcherror ps startbyte ret = parse_cell1d(ps)
    elseif isinstance(ps.t) || isoperator(ps.t)
        if ps.t.kind == Tokens.WHERE
            ret = IDENTIFIER(ps)
        else
            ret = INSTANCE(ps)
        end
        if (ret isa EXPR{OPERATOR{ColonOp,Tokens.COLON,false}}) && ps.nt.kind != Tokens.COMMA
            @catcherror ps startbyte ret = parse_unary(ps, ret)
        end
    elseif ps.t.kind == Tokens.AT_SIGN
        @catcherror ps startbyte ret = parse_macrocall(ps)
################################################################################
# Everything below here is an error
################################################################################
    elseif ps.t.kind == Tokens.ENDMARKER
        ps.errored = true
        return EXPR{ERROR}(EXPR[INSTANCE(ps)], 0, Variable[], "Unexpected end of file")
    elseif ps.t.kind == Tokens.COMMA
        ps.errored = true
        return EXPR{ERROR}(EXPR[INSTANCE(ps)], 0, Variable[], "Unexpected comma")
    elseif ps.t.kind == Tokens.RPAREN
        ps.errored = true
        return EXPR{ERROR}(EXPR[INSTANCE(ps)], 0, Variable[], "Unexpected )")
    elseif ps.t.kind == Tokens.RBRACE
        ps.errored = true
        return EXPR{ERROR}(EXPR[INSTANCE(ps)], 0, Variable[], "Unexpected }")
    elseif ps.t.kind == Tokens.RSQUARE
        ps.errored = true
        return EXPR{ERROR}(EXPR[INSTANCE(ps)], 0, Variable[], "Unexpected ]")
    else
        ps.errored = true
        return EXPR{ERROR}(EXPR[INSTANCE(ps)], 0, Variable[], "Unknown error")
    end

    while !closer(ps) && !(ps.closer.precedence == DotOp && ismacro(ret))
        @catcherror ps startbyte ret = parse_compound(ps, ret)
    end
    if ps.closer.precedence != DotOp && closer(ps) && ret isa EXPR{LITERAL{Tokens.MACRO}}
        ret = EXPR{MacroCall}(EXPR[ret], ret.span, Variable[], "")
    end

    return ret
end

"""
    parse_compound(ps, ret)

Handles cases where an expression - `ret` - is not followed by
`closer(ps) == true`. Possible juxtapositions are:
+ operators
+ `(`, calls
+ `[`, ref
+ `{`, curly
+ `,`, commas
+ `for`, generators
+ `do`
+ strings
+ an expression preceded by a unary operator
+ A number followed by an expression (with no seperating white space)
"""
function parse_compound(ps::ParseState, ret)
    startbyte = ps.nt.startbyte - ret.span
    if ps.nt.kind == Tokens.FOR
        ret = parse_generator(ps, ret)
    elseif ps.nt.kind == Tokens.DO
        ret = parse_do(ps, ret)
    elseif isajuxtaposition(ps, ret)
        op = EXPR{OPERATOR{TimesOp,Tokens.STAR,false}}(EXPR[], 0, Variable[], "")
        ret = parse_operator(ps, ret, op)
    elseif ps.nt.kind == Tokens.LPAREN && isemptyws(ps.ws)
        ret = @closer ps paren parse_call(ps, ret)
    elseif ps.nt.kind == Tokens.LBRACE && isemptyws(ps.ws)
        ret = parse_curly(ps, ret)
    elseif ps.nt.kind == Tokens.LSQUARE && isemptyws(ps.ws) && !(ret isa EXPR{OPERATOR{p,k,dot}} where {p, k, dot})
        ret = @nocloser ps block parse_ref(ps, ret)
    elseif ps.nt.kind == Tokens.COMMA
        ret = parse_tuple(ps, ret)
    elseif isunaryop(ret) && ps.nt.kind != Tokens.EQ
        ret = parse_unary(ps, ret)
    elseif isoperator(ps.nt)
        next(ps)
        op = INSTANCE(ps)
        format_op(ps, precedence(ps.t))
        ret = parse_operator(ps, ret, op)
    elseif (ret isa EXPR{IDENTIFIER} || (ret isa EXPR{BinarySyntaxOpCall} && ret.args[2] isa EXPR{OPERATOR{DotOp,Tokens.DOT,false}})) && (ps.nt.kind == Tokens.STRING || ps.nt.kind == Tokens.TRIPLE_STRING)
        next(ps)
        @catcherror ps startbyte arg = parse_string(ps, ret)
        ret = EXPR{x_Str}(EXPR[ret, arg], ret.span + arg.span, Variable[], "")
    # Suffix on x_str
    elseif ret isa EXPR{x_Str} && ps.nt.kind == Tokens.IDENTIFIER
        next(ps)
        arg = INSTANCE(ps)
        push!(ret.args, EXPR{LITERAL{Tokens.STRING}}(EXPR[], arg.span, Variable[], arg.val))
        ret.span += arg.span
    elseif (ret isa EXPR{IDENTIFIER} || (ret isa EXPR{BinarySyntaxOpCall} && ret.args[2] isa EXPR{OPERATOR{DotOp,Tokens.DOT,false}})) && ps.nt.kind == Tokens.CMD
        next(ps)
        @catcherror ps startbyte arg = parse_string(ps, ret)
        ret = EXPR{x_Cmd}(EXPR[ret, arg], ret.span + arg.span, Variable[], "")
    elseif ret isa EXPR{x_Cmd} && ps.nt.kind == Tokens.IDENTIFIER
        next(ps)
        arg = INSTANCE(ps)
        push!(ret.args, EXPR{LITERAL{Tokens.STRING}}(EXPR[], arg.span, Variable[], arg.val))
        ret.span += arg.span
    elseif ret isa EXPR{UnarySyntaxOpCall} && ret.args[2] isa EXPR{OPERATOR{16,Tokens.PRIME,d1}} where d1
        # prime operator followed by an identifier has an implicit multiplication
        @catcherror ps startbyte nextarg = @precedence ps 11 parse_expression(ps)
        ret = EXPR{BinaryOpCall}(EXPR[ret, EXPR{OPERATOR{TimesOp,Tokens.STAR,false}}(EXPR[], 0, Variable[], ""), nextarg], ret.span + nextarg.span, Variable[], "")
################################################################################
# Everything below here is an error
################################################################################
    elseif ps.nt.kind == Tokens.ENDMARKER
        ps.errored = true
        return EXPR{ERROR}(EXPR[INSTANCE(ps)], 0, Variable[], "Unexpected end of file")
    elseif ps.nt.kind == Tokens.LPAREN
        ps.errored = true
        return EXPR{ERROR}(EXPR[INSTANCE(ps)], 0, Variable[], "Unexpected (")
    elseif ps.nt.kind == Tokens.RPAREN
        ps.errored = true
        return EXPR{ERROR}(EXPR[INSTANCE(ps)], 0, Variable[], "Unexpected )")
    elseif ps.nt.kind == Tokens.LBRACE
        ps.errored = true
        return EXPR{ERROR}(EXPR[INSTANCE(ps)], 0, Variable[], "Unexpected {")
    elseif ps.nt.kind == Tokens.RBRACE
        ps.errored = true
        return EXPR{ERROR}(EXPR[INSTANCE(ps)], 0, Variable[], "Unexpected }")
    elseif ps.nt.kind == Tokens.LSQUARE
        ps.errored = true
        return EXPR{ERROR}(EXPR[INSTANCE(ps)], 0, Variable[], "Unexpected [")
    elseif ps.nt.kind == Tokens.RSQUARE
        ps.errored = true
        return EXPR{ERROR}(EXPR[INSTANCE(ps)], 0, Variable[], "Unexpected ]")
    elseif ret isa OPERATOR
        ps.errored = true
        return EXPR{ERROR}(EXPR[INSTANCE(ps)], 0, Variable[], "Unexpected operator")
    else
        ps.errored = true
        return EXPR{ERROR}(EXPR[INSTANCE(ps)], 0, Variable[], "Unknown error")
    end
    if ps.errored
        return EXPR{ERROR}(EXPR[INSTANCE(ps)], 0, Variable[], "Unknown error")
    end
    return ret
end

"""
    parse_list(ps)

Parses a list of comma seperated expressions finishing when the parent state
of `ps.closer` is met, newlines are ignored. Expects to start at the first item and ends on the last
item so surrounding punctuation must be handled externally.

**NOTE**
Should be replaced with the approach taken in `parse_call`
"""
function parse_list(ps::ParseState, puncs)
    startbyte = ps.nt.startbyte

    args = EXPR[]

    while !closer(ps)
        @catcherror ps startbyte a = @nocloser ps newline @closer ps comma parse_expression(ps)
        push!(args, a)
        if ps.nt.kind == Tokens.COMMA
            next(ps)
            push!(puncs, INSTANCE(ps))
            format_comma(ps)
        end
    end

    if ps.t.kind == Tokens.COMMA
        format_comma(ps)
    end
    return args
end

"""
    parse_paren(ps, ret)

Parses an expression starting with a `(`.
"""
function parse_paren(ps::ParseState)
    startbyte = ps.t.startbyte

    ret = EXPR{TupleH}(EXPR[INSTANCE(ps)], - startbyte, Variable[], "")
    format_lbracket(ps)

    @catcherror ps startbyte @default ps @nocloser ps inwhere @closer ps paren parse_comma_sep(ps, ret, false, true)

    if length(ret.args) == 2 && !(ret.args[2] isa EXPR{UnarySyntaxOpCall} && ret.args[2].args[2] isa EXPR{OPERATOR{DddotOp,Tokens.DDDOT,false}})

        if ps.ws.kind != SemiColonWS || (length(ret.args) == 2 && ret.args[2] isa EXPR{Block})
            ret = EXPR{InvisBrackets}(ret.args, ret.span, Variable[], "")
        end
    end

    # handle closing ')'
    next(ps)
    push!(ret.args, INSTANCE(ps))
    format_rbracket(ps)

    ret.span = ps.nt.startbyte - startbyte

    return ret
end

"""
    parse(str, cont = false)

Parses the passed string. If `cont` is true then will continue parsing until the end of the string returning the resulting expressions in a TOPLEVEL block.
"""
function parse(str::String, cont = false)
    ps = ParseState(str)
    x, ps = parse(ps, cont)
    if ps.errored
        x = EXPR{ERROR}(EXPR[], 0, Variable[], "Unknown error")
    end
    return x
end

function parse_doc(ps::ParseState)
    if ps.nt.kind == Tokens.STRING || ps.nt.kind == Tokens.TRIPLE_STRING
        next(ps)
        doc = INSTANCE(ps)
        if (ps.nt.kind == Tokens.ENDMARKER || ps.nt.kind == Tokens.END)
            return doc
        elseif isbinaryop(ps.nt) && !closer(ps)
            @catcherror ps startbyte ret = parse_compound(ps, doc)
            return ret
        end

        ret = parse_expression(ps)
        ret = EXPR{MacroCall}(EXPR[GlobalRefDOC, doc, ret], doc.span + ret.span, Variable[], "")
    elseif ps.nt.kind == Tokens.IDENTIFIER && ps.nt.val == "doc" && (ps.nnt.kind == Tokens.STRING || ps.nnt.kind == Tokens.TRIPLE_STRING)
        next(ps)
        doc = INSTANCE(ps)
        next(ps)
        @catcherror ps startbyte arg = parse_string(ps, doc)
        doc = EXPR{x_Str}(EXPR[doc, arg], doc.span + arg.span, Variable[], "")
        ret = parse_expression(ps)
        ret = EXPR{MacroCall}(EXPR[GlobalRefDOC, doc, ret], doc.span + ret.span, Variable[], "")
    else
        ret = parse_expression(ps)
    end
    return ret
end

function parse(ps::ParseState, cont = false)
    if ps.l.io.size == 0
        return (cont ? EXPR{FileH}(EXPR[], 0, Variable[], "") : nothing), ps
    end
    last_line = 0
    curr_line = 0

    if cont
        top = EXPR{FileH}(EXPR[], 0, Variable[], "")
        if ps.nt.kind == Tokens.WHITESPACE || ps.nt.kind == Tokens.COMMENT
            next(ps)
            push!(top.args, EXPR{LITERAL{nothing}}(EXPR[], ps.nt.startbyte, Variable[], "comments"))
        end

        while !ps.done && !ps.errored
            curr_line = ps.nt.startpos[1]
            ret = parse_doc(ps)

            # join semicolon sep items
            if curr_line == last_line && last(top.args) isa EXPR{TopLevel}
                push!(last(top.args).args, ret)
                last(top.args).span += ret.span
            elseif ps.ws.kind == SemiColonWS
                push!(top.args, EXPR{TopLevel}(EXPR[ret], ret.span, Variable[], ""))
            else
                push!(top.args, ret)
            end
            last_line = curr_line
        end
        top.span += ps.nt.startbyte
    else
        if ps.nt.kind == Tokens.WHITESPACE || ps.nt.kind == Tokens.COMMENT
            next(ps)
            top = EXPR{LITERAL{nothing}}(EXPR[], ps.nt.startbyte, Variable[], "comments")
        else
            top = parse_doc(ps)
            last_line = ps.nt.startpos[1]
            if ps.ws.kind == SemiColonWS
                top = EXPR{TopLevel}(EXPR[top], top.span, Variable[], "")
                while ps.ws.kind == SemiColonWS && ps.nt.startpos[1] == last_line && ps.nt.kind != Tokens.ENDMARKER
                    ret = parse_doc(ps)
                    push!(top.args, ret)
                    top.span += ret.span
                    last_line = ps.nt.startpos[1]
                end
            end
        end
    end

    return top, ps
end

function parse_file(path::String)
    x = parse(readstring(path), true)
    File([], [], path, x, [])
    # File([], (f -> (joinpath(dirname(path), f[1]), f[2])).(_get_includes(x)), path, x, [])
end

function parse_directory(path::String, proj = Project(path, []))
    for f in readdir(path)
        if isfile(joinpath(path, f)) && endswith(f, ".jl")
            try
                push!(proj.files, parse_file(joinpath(path, f)))
            catch
                println("$f failed to parse")
            end
        elseif isdir(joinpath(path, f))
            parse_directory(joinpath(path, f), proj)
        end
    end
    proj
end

ischainable(t::Token) = t.kind == Tokens.PLUS || t.kind == Tokens.STAR || t.kind == Tokens.APPROX

# include("precompile.jl")
# _precompile_()
end


text:
# add CSV      # pkg mode
using CSV, DataFrames
fname = "todos.csv"
df = CSV.read(fname, DataFrame)

ext: jl


text:
struct NNPDEHan{C1,C2,O} <: NeuralNetDiffEqAlgorithm
    u0::C1
    σᵀ∇u::C2
    opt::O
end
NNPDEHan(u0,σᵀ∇u;opt=Flux.ADAM(0.1)) = NNPDEHan(u0,σᵀ∇u,opt)

function DiffEqBase.solve(
    prob::TerminalPDEProblem,
    alg::NNPDEHan;
    abstol = 1f-6,
    verbose = false,
    maxiters = 300,
    save_everystep = false,
    dt,
    trajectories)

    X0 = prob.X0
    ts = prob.tspan[1]:dt:prob.tspan[2]
    d  = length(X0)
    g,f,μ,σ,p = prob.g,prob.f,prob.μ,prob.σ,prob.p

    data = Iterators.repeated((), maxiters)

    #hidden layer
    opt = alg.opt
    u0 = alg.u0
    σᵀ∇u = alg.σᵀ∇u
    ps = Flux.params(u0, σᵀ∇u...)

    function sol()
        map(1:trajectories) do j
            u = u0(X0)[1]
            X = X0
            for i in 1:length(ts)-1
                t = ts[i]
                _σᵀ∇u = σᵀ∇u[i](X)
                dW = sqrt(dt)*randn(d)
                u = u - f(X, u, _σᵀ∇u, p, t)*dt + _σᵀ∇u'*dW
                X  = X .+ μ(X,p,t)*dt .+ σ(X,p,t)*dW
            end
            X,u
        end
    end

    function loss()
        mean(sum(abs2,g(X) - u) for (X,u) in sol())
    end

    iters = eltype(X0)[]

    cb = function ()
        save_everystep && push!(iters, u0(X0)[1].data)
        l = loss()
        verbose && println("Current loss is: $l")
        l < abstol && Flux.stop()
    end

    Flux.train!(loss, ps, data, opt; cb = cb)

    save_everystep ? iters : u0(X0)[1].data
end #pde_solve


text:
using TerminalUserInterfaces
const TUI = TerminalUserInterfaces
using Test

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


text:
using Behavior
using Niancat.Gameface
using Niancat.Users
using Niancat.Games.NiancatGames

@when("Erik sends the command \"{String}\"") do context, commandstring
    svc = context[:service]
    user = context[:user]

    game, command = findcommand(svc, user, commandstring)
    context[:foundgame] = game
    context[:foundcommand] = command
end

@then("the service finds the game Niancat in the default instance") do context
    foundgame = context[:foundgame]
    expectedgame = context[:game]

    @expect gameinstanceid(foundgame) == gameinstanceid(expectedgame)
end

@then("the command is GetPuzzle()") do context
    command = context[:foundcommand]

    @expect command == NiancatGames.GetPuzzle()
end

@then("the command is Guess(\"{String}\")") do context, word
    command = context[:foundcommand]

    @expect command == NiancatGames.Guess(word)
end

@then("the command is SetPuzzle(\"{String}\")") do context, word
    command = context[:foundcommand]

    @expect command == NiancatGames.SetPuzzle(word)
end


text:
"""
$(TYPEDEF)
$(TYPEDFIELDS)

"""
struct CategoricalIndex <: AbstractDistDelayed
    probabilities::Vector{Float64}
end

function categoricalindex(
    probabilities::Vector{Float64}, sample_size::Int
)::IndexObjects.IndexVector
    iszero(sample_size) && return []

    sample = transpose(rand(Distributions.Multinomial(1, probabilities), sample_size))
    rval = sample * (1:length(probabilities))

    return IndexObjects.IndexVector(rval)
end


text:
using CairoMakie
using CSV

f = CSV.File("LinearSampling.csv")

random_sampling_trace = f.random_sampling_trace
optimal_sampling_trace_low = f.optimal_sampling_trace_low
optimal_sampling_trace_medium = f.optimal_sampling_trace_medium
samples = f.samples

noto_sans = "../resources/NotoSans-Regular.ttf"

tickfontsize    = 28
labelfontsize   = 32

f = Figure(font=noto_sans, resolution=(600,600))

ax = Axis(  f[1,1],
            xticklabelsize=tickfontsize,
            yticklabelsize=tickfontsize,
            yticklabelpad=2,
            xlabel="Sample", ylabel = "Tr(Σ)",
            xlabelsize=labelfontsize,
            ylabelsize=labelfontsize,
            yscale = log10,
            yminorticksvisible = true,
            yminorgridvisible = true,
            yminorticks = IntervalsBetween(10),
            xminorticksvisible = true,
            xminorgridvisible = true,
            xminorticks = IntervalsBetween(10))

l1 = lines!(ax, 1:max_steps, optimal_sampling_trace_low, color=:black)
l2 = lines!(ax, 1:max_steps, random_sampling_trace, color=:black, linestyle=:dash)
l3 = lines!(ax, 1:max_steps, random_sampling_trace ./ 2, color=:black, linestyle = :dot)

Legend(f[2,:],
       [l1,l2,l3],
       ["Dynamic Programming", "Random", "Random/Dimension"],
       orientation=:horizontal,
       labelsize=24)

ylims!(ax,(1e-3,1))
xlims!(ax,(0,1000))

save("LinearSampling.pdf",f)


text:
# ---
# title: 50. Pow(x, n)
# id: problem50
# author: zhwang
# date: 2022-03-02
# difficulty: Medium
# categories: Math, Binary Search
# link: <https://leetcode.com/problems/powx-n/description/>
# hidden: true
# ---
#
# Implement [pow( _x_ , _n_
# )](http://www.cplusplus.com/reference/valarray/pow/), which calculates  _x_
# raised to the power _n_ (i.e. xn).
#
#
#
# **Example 1:**
#
#
#
#     Input: x = 2.00000, n = 10
#     Output: 1024.00000
#
#
# **Example 2:**
#
#
#
#     Input: x = 2.10000, n = 3
#     Output: 9.26100
#
#
# **Example 3:**
#
#
#
#     Input: x = 2.00000, n = -2
#     Output: 0.25000
#     Explanation: 2-2 = 1/22 = 1/4 = 0.25
#
#
#
#
# **Constraints:**
#
#   * `-100.0 < x < 100.0`
#   * `-231 <= n <= 231-1`
#   * `-104 <= xn <= 104`
#
#
## @lc code=start
using LeetCode

function my_pow(x::Float64, n::Int)::Float64
    n == 0 && return 1
    if n < 0
        x, n = 1 / x, -n
    end
    res, extra = x, 1
    while n > 1
        n & 1 == 1 && (extra *= res)
        res *= res
        n >>= 1
    end
    return res * extra
end
## @lc code=end


text:
using PkgName, Base.Test

function ortho(V::AbstractMatrix{T})where{T}
    n = size(V,2)
    t = true
    for i in 1:n-1
        for j in (i+1):n
            if(vecdot(view(V, :, i),view(V, :, j))>(10*eps(T)))
                t = false
                break
            end
        end
    end
    return t
end

@test ortho(hessenberg(sprand(100,100,0.3),30)[2]) == true


text:
"""
```
AbstractVARModel{T} <: AbstractModel{T}
```

The AbstractVARModel is defined as a subtype of AbstractModel to implement
vector autoregression methods.
"""
abstract type AbstractVARModel{T} <: ModelConstructors.AbstractModel{T} end

function Base.show(io::IO, m::AbstractVARModel)
    @printf io "VAR Model\n"
    @printf io "observables:   %i\n" observables(m)
    @printf io "data_vintage:  %i\n" data_vintage(m)
    @printf io "description:\n %s\n" description(m)
end

get_parameters(m::AbstractVARModel) = m.parameters

"""
```
AbstractDSGEVARModel{T} <: AbstractVARModel{T}
```
The AbstractDSGEVARModel is defined as a subtype of AbstractVARModel to implement
DSGE-VAR methods.
"""
abstract type AbstractDSGEVARModel{T} <: AbstractVARModel{T} end

function Base.show(io::IO, m::AbstractDSGEVARModel)
    @printf io "DSGE-VAR Model\n"
    @printf io "observables:      %i\n" n_observables(m)
    @printf io "data_vintage:     %s\n" data_vintage(m)
    @printf io "DSGE model:       %s\n" spec(get_dsge(m))
    @printf io "DSGE description: %s\n" description(get_dsge(m))
end

"""
```
AbstractDSGEVECMModel{T} <: AbstractDSGEVARModel{T}
```
The AbstractDSGEVECMModel is defined as a subtype of AbstractDSGEVARModel to implement
DSGE-VECM methods.
"""
abstract type AbstractDSGEVECMModel{T} <: AbstractDSGEVARModel{T} end

# Number of anticipated policy shocks
n_anticipated_shocks(m::AbstractDSGEVARModel) = get_setting(m, :n_anticipated_shocks)
n_anticipated_shocks_padding(m::AbstractDSGEVARModel) = get_setting(m, :n_anticipated_shocks_padding)

# Dates, indices, number of periods for each regime
date_presample_start(m::AbstractDSGEVARModel) = get_setting(m, :date_presample_start)
date_mainsample_start(m::AbstractDSGEVARModel) = get_setting(m, :date_mainsample_start)
date_zlb_start(m::AbstractDSGEVARModel) = get_setting(m, :date_zlb_start)
# date_zlb_end(m::AbstractDSGEVARModel) = get_setting(m, :date_zlb_end)

date_presample_end(m::AbstractDSGEVARModel) = Dates.lastdayofquarter(get_setting(m, :date_mainsample_start) - Dates.Month(3))
date_prezlb_end(m::AbstractDSGEVARModel) = Dates.lastdayofquarter(get_setting(m, :date_zlb_start) - Dates.Month(3))
date_mainsample_end(m::AbstractDSGEVARModel) = Dates.lastdayofquarter(get_setting(m, :date_forecast_start) - Dates.Month(3))
date_conditional_end(m::AbstractDSGEVARModel) = get_setting(m, :date_conditional_end)

index_presample_start(m::AbstractDSGEVARModel) = 1
index_mainsample_start(m::AbstractDSGEVARModel) = subtract_quarters(date_mainsample_start(m), date_presample_start(m)) + 1
index_zlb_start(m::AbstractDSGEVARModel) = subtract_quarters(date_zlb_start(m), date_presample_start(m)) + 1
index_forecast_start(m::AbstractDSGEVARModel) = subtract_quarters(date_forecast_start(m), date_presample_start(m)) + 1

# Interface for fields and size
n_lags(m::AbstractDSGEVARModel)          = m.lags
get_lags(m::AbstractDSGEVARModel)        = m.lags
get_shocks(m::AbstractDSGEVARModel)      = m.shocks
get_λ(m::AbstractDSGEVARModel)           = m.λ
get_dsge(m::AbstractDSGEVARModel)        = m.dsge
n_observables(m::AbstractDSGEVARModel)   = length(m.observables)
n_cointegrating(m::AbstractDSGEVECMModel) = length(m.cointegrating)
n_cointegrating_add(m::AbstractDSGEVECMModel) = length(m.cointegrating_add)
n_shocks(m::AbstractDSGEVARModel)        = length(m.shocks)

# Interface for accessing parameters
get_parameters(m::AbstractDSGEVARModel) = m.dsge.parameters
n_parameters(m::AbstractDSGEVARModel)   = n_parameters(m.dsge)
n_parameters_steady_state(m::AbstractDSGEVARModel) = n_parameters_steady_state(m.dsge)
n_parameters_free(m::AbstractDSGEVARModel) = n_parameters_free(m.dsge)

# Interface for accessing rng
get_rng(m::AbstractDSGEVARModel) = m.dsge.rng

# Interface for accessing settings dictionary
get_settings(m::AbstractDSGEVARModel) = m.dsge.settings

# Interface for accessing observables dictionary
get_observables(m::AbstractDSGEVARModel) = m.observables

# Interface for accessing cointegrating dictionaries
get_cointegrating(m::AbstractDSGEVECMModel) = m.cointegrating
get_cointegrating_add(m::AbstractDSGEVECMModel) = m.cointegrating_add

# Interface for data
cond_vintage(m::AbstractDSGEVARModel)    = get_setting(m, :cond_vintage)
cond_id(m::AbstractDSGEVARModel)         = get_setting(m, :cond_id)
cond_full_names(m::AbstractDSGEVARModel) = get_setting(m, :cond_full_names)
cond_semi_names(m::AbstractDSGEVARModel) = get_setting(m, :cond_semi_names)
use_population_forecast(m::AbstractDSGEVARModel) = get_setting(m, :use_population_forecast)
hpfilter_population(m::AbstractDSGEVARModel)     = get_setting(m, :hpfilter_population)
n_conditional_periods(m::AbstractDSGEVARModel)   = n_conditional_periods(get_dsge(m))

# Interface for general computation settings
use_parallel_workers(m::AbstractDSGEVARModel)    = get_setting(m, :use_parallel_workers)

# Interface for estimation settings
reoptimize(m::AbstractDSGEVARModel)          = get_setting(m, :reoptimize)
calculate_hessian(m::AbstractDSGEVARModel) = get_setting(m, :calculate_hessian)
hessian_path(m::AbstractDSGEVARModel)      = get_setting(m, :hessian_path)
n_hessian_test_params(m::AbstractDSGEVARModel) = get_setting(m, :n_hessian_test_params)

# Interface for Metropolis-Hastings settings
n_mh_blocks(m::AbstractDSGEVARModel)       =  get_setting(m, :n_mh_blocks)
n_mh_param_blocks(m::AbstractDSGEVARModel) =  get_setting(m, :n_mh_param_blocks)
n_mh_simulations(m::AbstractDSGEVARModel)  =  get_setting(m, :n_mh_simulations)
n_mh_burn(m::AbstractDSGEVARModel)         =  get_setting(m, :n_mh_burn)
mh_thin(m::AbstractDSGEVARModel)           =  get_setting(m, :mh_thin)

# Interface for forecast settings
date_forecast_start(m::AbstractDSGEVARModel)   = get_setting(m, :date_forecast_start)
forecast_block_size(m::AbstractDSGEVARModel)   = get_setting(m, :forecast_block_size)
forecast_start_block(m::AbstractDSGEVARModel)  = get_setting(m, :forecast_start_block)
forecast_input_file_overrides(m::AbstractDSGEVARModel) = get_setting(m, :forecast_input_file_overrides)
forecast_uncertainty_override(m::AbstractDSGEVARModel) = get_setting(m, :forecast_uncertainty_override)
forecast_smoother(m::AbstractDSGEVARModel)     = get_setting(m, :forecast_smoother)
forecast_tdist_df_val(m::AbstractDSGEVARModel) = get_setting(m, :forecast_tdist_df_val)
forecast_tdist_shocks(m::AbstractDSGEVARModel) = get_setting(m, :forecast_tdist_shocks)
forecast_zlb_value(m::AbstractDSGEVARModel)    = get_setting(m, :forecast_zlb_value)
impulse_response_horizons(m::AbstractDSGEVARModel) = get_setting(m, :impulse_response_horizons)
n_shockdec_periods(m::AbstractDSGEVARModel)    = index_shockdec_end(m) - index_shockdec_start(m) + 1

# Interface for alternative policy settings
alternative_policy(m::AbstractDSGEVARModel) = get_setting(m, :alternative_policy)

function date_forecast_end(m::AbstractDSGEVARModel)
    date = date_forecast_start(m) + Dates.Month(3 * (forecast_horizons(m)-1))
    return Dates.lastdayofquarter(date)
end

function forecast_horizons(m::AbstractDSGEVARModel; cond_type::Symbol = :none)
    horizons = get_setting(m, :forecast_horizons)
    if cond_type == :none
        return horizons
    else
        return horizons - n_conditional_periods(m)
    end
end

function date_shockdec_start(m::AbstractDSGEVARModel)
    startdate = get_setting(m, :shockdec_startdate)
    if (Nullables.isnull(startdate)) | (startdate===nothing)
        return DSGEVAR.date_mainsample_start(m)
    else
        return get(startdate)
    end
end

function date_shockdec_end(m::AbstractDSGEVARModel)
    enddate =  get_setting(m, :shockdec_enddate)
    if !isnull(enddate)
        return get(enddate)
    else
        return date_forecast_end(m)
    end
end

"""
```
transform_to_model_space!(m::AbstractDSGEVARModel, values::Vector{T}) where T<:AbstractFloat
```

Transforms `values` from the real line to the model space, and assigns `values[i]` to
`m.parameters[i].value` for non-steady-state parameters. Recomputes the steady-state
paramter values.

### Arguments
- `m`: the model object
- `values`: the new values to assign to non-steady-state parameters.
"""
function transform_to_model_space!(m::AbstractDSGEVARModel, values::Vector{T}) where {T<:AbstractFloat}
    new_values = transform_to_model_space(m.dsge.parameters, values)
    DSGE.update!(m, new_values)
    steadystate!(get_dsge(m))
end

# Updating parameters
function update!(m::AbstractDSGEVARModel{T}, values::Vector{T}) where {T<:Real}
    DSGE.update!(get_dsge(m), values)
    steadystate!(get_dsge(m))
end

function update!(m::AbstractDSGEVARModel{T}, values::ParameterVector{T}) where {T<:Real}
    DSGE.update!(get_dsge(m), values)
    steadystate!(get_dsge(m))
end

# Overload <= so that it applies to the underlying DSGE object
function (<=)(m::AbstractDSGEVARModel{T}, p::ModelConstructors.AbstractParameter{T}) where {T}
    dsge = get_dsge(m)
    dsge <= p
end

function (<=)(m::AbstractDSGEVARModel{T}, p::Union{ModelConstructors.SteadyStateParameter, ModelConstructors.SteadyStateParameterArray}) where {T}
    dsge = get_dsge(m)
    dsge <= p
end

function (<=)(m::AbstractDSGEVARModel{T}, p::ModelConstructors.SteadyStateParameterGrid) where {T}
    dsge = get_dsge(m)
    dsge <= p
end

function (<=)(m::AbstractDSGEVARModel, s::Setting)
    dsge = get_dsge(m)
    dsge <= s
end

# Setting access
function get_setting(m::AbstractDSGEVARModel, k::Symbol)
    return ModelConstructors.get_setting(get_dsge(m), k)
end

# Prior
function prior(m::AbstractDSGEVARModel)
    return ModelConstructors.prior(get_dsge(m))
end

# Not exposed to user. Actually create path and insert model string to file name.
function savepath(m::AbstractDSGEVARModel,
                  out_type::String,
                  sub_type::String,
                  file_name::String = "",
                  filestring_addl::Vector{String} = Vector{String}())
    # Containing directory
    dir = String(joinpath(saveroot(m), "output_data", spec(get_dsge(m)), subspec(get_dsge(m)),
                          spec(m), subspec(m), out_type, sub_type))

    if !isempty(file_name)
        base = filestring_base(m)
        return savepath(dir, file_name, base, filestring_addl)
    else
        return dir
    end
end

function savepath(dir::String,
                  file_name::String = "",
                  filestring_base::Vector{String} = Vector{String}(),
                  filestring_addl::Vector{String} = Vector{String}())
    if !isdir(dir)
        mkpath(dir)
    end

    if !isempty(file_name)
        (base, ext) = splitext(file_name)
        myfilestring = filestring(filestring_base, filestring_addl)
        file_name_detail = base * myfilestring * ext

        return joinpath(dir, file_name_detail)
    else
        return dir
    end
end

function filestring_base(m::AbstractDSGEVARModel)
    if !m.testing
        base = Vector{String}()
        for (skey, sval) in get_settings(m)
            if sval.print
                push!(base, ModelConstructors.to_filestring(sval))
            end
        end
        return base
    else
        return ["test"]
    end
end

# Helper to check for valid VAR systems
function check_valid_dsgevar(m::AbstractDSGEVARModel;
                             observables::Vector{Symbol} = collect(keys(m.observables)),
                             shocks::Vector{Symbol} = collect(keys(m.shocks)),
                             lags::Int = n_lags(m), λ::T = m.λ) where {T<:Real}
    # Check lags
    @assert lags >= 0 "Number of lags must be non-negative."

    # Check λ is non-negative
    @assert λ >= 0 "The weight on the DSGE prior λ must be non-negative."

    # Check observables
    missing_obs = Vector{Symbol}(undef, 0)
    for k in observables
        if !(k in keys(get_observables(get_dsge(m)))) &&
            !(k in keys(get_pseudo_observables(get_dsge(m))))
            push!(missing_obs, k)
        end
    end

    if !isempty(missing_obs)
        error("The following observables are not found in the underlying DSGE: \n" *
              join(string.(missing_obs), ", "))
    end

    # Check shocks
    missing_shocks = Vector{Symbol}(undef, 0)
    for k in shocks
        if !(k in keys(get_exogenous_shocks(get_dsge(m))))
            push!(missing_shocks, k)
        end
    end

    if !isempty(missing_shocks)
        error("The following shocks are not found in the underlying DSGE: \n" *
              join(string.(missing_shocks), ", "))
    end
end

# Helper to check for valid VECM systems
function check_valid_dsgevecm(m::AbstractDSGEVECMModel;
                              observables::Vector{Symbol} = collect(keys(m.observables)),
                              cointegrating::Vector{Symbol} = Vector{Symbol}(undef, 0),
                              cointegrating_add::Vector{Symbol} = Vector{Symbol}(undef, 0),
                              shocks::Vector{Symbol} = collect(keys(m.shocks)),
                              lags::Int = n_lags(m), λ::T = m.λ) where {T<:Real}

    # Check non-cointegrating fields
    check_valid_dsgevar(m; observables = observables,
                        shocks = shocks, lags = lags, λ = λ)

    # Check cointegrated observables
    missing_coint = Vector{Symbol}(undef, 0)
    for k in cointegrating
        if !(k in keys(get_observables(get_dsge(m)))) &&
            !(k in keys(get_pseudo_observables(get_dsge(m))))
            push!(missing_coint, k)
        end
    end

    if !isempty(missing_coint)
        error("The following cointegrating relationships are not found in the underlying DSGE: \n" *
              join(string.(missing_coint), ", "))
    end

    # Check additional cointegrated observables
    missing_coint = Vector{Symbol}(undef, 0)
    for k in cointegrating
        if !(k in keys(get_observables(get_dsge(m)))) &&
            !(k in keys(get_pseudo_observables(get_dsge(m))))
            push!(missing_coint, k)
        end
    end

    if !isempty(missing_coint)
        error("The following cointegrating relationships are not found in the underlying DSGE: \n" *
              join(string.(missing_coint), ", "))
    end
end

# Update VAR system, e.g. observables
function update!(m::AbstractDSGEVARModel; observables::Vector{Symbol} = Vector{Symbol}(undef, 0),
                 shocks::Vector{Symbol} = Vector{Symbol}(undef, 0),
                 lags::Int = 0, λ::T = m.λ, check_valid::Bool = true) where {T<:Real}

    if check_valid
        check_obs = isempty(observables) ? collect(keys(get_observables(m))) : observables
        check_sh  = isempty(shocks) ? collect(keys(get_shocks(m))) : shocks
        check_valid_dsgevar(m; observables = check_obs, shocks = check_sh,
                            lags = lags, λ = λ)
    end

    if m.λ != λ
        m.λ = λ
    end
    if !isempty(observables)
        for k in keys(m.observables)
            delete!(m.observables, k)
        end
        for (i,k) in enumerate(observables)
            m.observables[k] = i
        end
    end
    if !isempty(shocks)
        for k in keys(m.shocks)
            delete!(m.shocks, k)
        end
        for (i,k) in enumerate(shocks)
            m.shocks[k] = i
        end
    end
    if lags > 0
        m.lags = lags
    end

    return m
end

function update!(m::AbstractDSGEVARModel, dsge::AbstractDSGEModel;
                 observables::Vector{Symbol} = Vector{Symbol}(undef, 0),
                 shocks::Vector{Symbol} = Vector{Symbol}(undef, 0),
                 lags::Int = 0, λ::T = m.λ, check_valid::Bool = true) where {T<:Real}

    m.dsge = dsge
    update!(m; observables = observables, shocks = shocks, lags = lags,
            λ = λ, check_valid = check_valid)

    return m
end

# Update VECM system, e.g. observables
function update!(m::AbstractDSGEVECMModel; observables::Vector{Symbol} = Vector{Symbol}(undef, 0),
                 cointegrating::Vector{Symbol} = Vector{Symbol}(undef, 0),
                 cointegrating_add::Vector{Symbol} = Vector{Symbol}(undef, 0),
                 shocks::Vector{Symbol} = Vector{Symbol}(undef, 0),
                 lags::Int = 0, λ::T = m.λ, check_valid::Bool = true) where {T<:Real}

    if check_valid
        check_obs = isempty(observables) ? collect(keys(get_observables(m))) : observables
        check_sh  = isempty(shocks) ? collect(keys(get_shocks(m))) : shocks
        check_co  = isempty(cointegrating) ? collect(keys(get_cointegrating(m))) : cointegrating
        check_coa = isempty(cointegrating_add) ? collect(keys(get_cointegrating_add(m))) : cointegrating_add

        check_valid_dsgevecm(m; observables = check_obs,
                             cointegrating = check_co,
                             cointegrating_add = check_coa,
                             shocks = check_sh, lags = lags, λ = λ)
    end

    if m.λ != λ
        m.λ = λ
    end
    if !isempty(observables)
        for k in keys(get_observables(m))
            delete!(get_observables(m), k)
        end
        for (i,k) in enumerate(observables)
            m.observables[k] = i
        end
    end
    if !isempty(cointegrating)
        n_obs = length(m.observables)
        for k in keys(get_cointegrating(m))
            delete!(get_cointegrating(m), k)
        end
        for (i,k) in enumerate(cointegrating)
            m.cointegrating[k] = i + n_obs
        end
    end
    if !isempty(cointegrating_add)
        for k in keys(get_cointegrating_add(m))
            delete!(get_cointegrating_add(m), k)
        end
        for (i,k) in enumerate(cointegrating_add)
            m.cointegrating_add[k] = i
        end
    end
    if !isempty(shocks)
        for k in keys(m.shocks)
            delete!(m.shocks, k)
        end
        for (i,k) in enumerate(shocks)
            m.shocks[k] = i
        end
    end
    if lags > 0
        m.lags = lags
    end

    return m
end

function update!(m::AbstractDSGEVECMModel, dsge::AbstractDSGEModel;
                 observables::Vector{Symbol} = Vector{Symbol}(undef, 0),
                 cointegrating::Vector{Symbol} = Vector{Symbol}(undef, 0),
                 cointegrating_add::Vector{Symbol} = Vector{Symbol}(undef, 0),
                 shocks::Vector{Symbol} = Vector{Symbol}(undef, 0),
                 lags::Int = 0, λ::T = m.λ, check_valid::Bool = true) where {T<:Real}

    m.dsge = dsge
    update!(m; observables = observables, cointegrating = cointegrating,
            cointegrating_add = cointegrating_add,
            shocks = shocks, lags = lags,
            λ = λ, check_valid = check_valid)

    return m
end


text:
module ReuseDistance

export Node, Treap, search

# deps
import AbstractTrees
import ProgressMeter
import Setfield: @set

# implementation
include("reference.jl")
include("treap.jl")
include("redblack.jl")

function reuse(itr)
    #treap = Treap{Tuple{Int,eltype(itr)}}()
    treap = RedBlackTree{Tuple{Int,eltype(itr)}}()
    lastuse = Dict{eltype(itr),Int}()
    histogram = Dict{Int,Int}()

    ProgressMeter.@showprogress 1 for (t, i) in enumerate(itr)
        lasttime = get(lastuse, i, t)
        key = (lasttime, i)
        d = -1
        if lasttime != t
            d = nchildren(treap, key)
            delete!(treap, key)
        end
        histogram[d] = 1 + get(histogram, d, 0)
        push!(treap, (t, i))
        lastuse[i] = t
    end
    @assert istree(treap)
    # @assert isheap(treap)
    return histogram
end

end


text:
@safetestset "Body Init Tests" begin
    include("body_initialization_test.jl")
end

@safetestset "Fixed Joint Tests" begin
    include("fixed_joint_test.jl")
end

@safetestset "Prismatic Joint Tests" begin
    include("prismatic_joint_test.jl")
end

@safetestset "Revolute Joint Tests" begin
    include("revolute_joint_test.jl")
end

# include("initialization/cylindrical_joint_test.jl")
# include("initialization/spherical_joint_test.jl")
# include("initialization/cylindricalfree_joint_test.jl")
# include("initialization/planarfree_joint_test.jl")
# include("initialization/origincon_joint_test.jl")


text:
using ReTest # useful when setup.jl included separately

module Trace
export trace
using ReTest: @test

const RUN = []

function trace(x)
    push!(RUN, x)
    @test true
end
end

function check(x...; runtests=false, output::Union{Nothing,String}=nothing,
               verbose=true, stats=false, dry=false, strict::Bool=true,
               recursive=true, static=nothing, id=nothing, load=false)
    @assert !(runtests & (output !== nothing)) "unimplemented"
    args = x[1:end-1]
    expected = x[end]
    if expected isa AbstractString
        expected = split(x[end])
    end

    empty!(Trace.RUN)
    if runtests
        getfield(args[1], :runtests)(args[2:end]...; verbose=verbose, stats=stats, dry=dry,
                                     strict=strict, recursive=recursive, static=static,
                                     id=id, load=load)
    elseif output === nothing
        retest(args...; verbose=verbose, stats=stats, dry=dry, strict=strict,
               recursive=recursive, static=static, id=id, load=load)
    else
        mktemp() do path, io
            redirect_stdout(io) do
                retest(args...; verbose=verbose, stats=stats, dry=dry, strict=strict,
                       recursive=recursive, static=static, id=id, load=load)
            end
            seekstart(io)
            printed = join(map(rstrip, split(readchomp(io), '\n')), '\n')
            output = join(map(rstrip, split(chomp(output), '\n')), '\n')
            @test printed == output
        end
    end
    @test Trace.RUN == expected
end

macro chapter(title, x)
    title = string(title)

    if isempty(ARGS) || any(pat -> occursin(Regex(pat, "i"), title), ARGS)
        printstyled("\n\n", rpad("## $title #", 78, '#'), "\n\n", bold=true, color=:cyan)
        esc(x)
    end
end


text:
using ConcaveHull
using CurveProximityQueries
using LinearAlgebra
using Plots
import Random: seed!

pgfplots()
seed!(19)

@recipe function f(a::ConcaveHull.Hull)
    grid := false
    axis := false
    ticks := false
    legend := false
    aspect_ratio := :equal
    linewidth --> 1.5
    linecolor --> RGBA(220/255, 20/255, 60/255, 0.6) #Gray(0.7)
    linestyle --> :dash
    leg := false
    fill --> (0, RGBA(220/255, 20/255, 60/255, 0.1))
    @series begin
        pts = hcat(a.vertices...)
        x = vcat(pts[1, :], pts[1, 1])
        y = vcat(pts[2, :], pts[2, 1])
        x, y
    end
end

a = rand(Bernstein{2,12})
cp = a.f.control_points .- 0.5
a = Bernstein(cp..., 2*cp[1]-cp[2], cp[1])

b = Bernstein([[0.95, -0.2], [0.85, -0.2], [0.85, -0.1], [0.9125, -0.1], [0.975, -0.1], [1.0375, -0.1], [1.1, -0.1], [1.1, -0.1875], [1.1, -0.275], [1.1, -0.3625], [1.1, -0.45], [0.975, -0.45], [0.85, -0.45], [0.725, -0.45], [0.6, -0.45], [0.6, -0.3375], [0.6, -0.225], [0.6, -0.1125], [0.6, 0.0], [0.675, 0.0], [0.75, 0.0], [0.825, 0.0], [0.9, 0.0], [0.9, 0.1], [0.8, 0.1], [0.75, 0.1], [0.7, 0.1], [0.7, 0.05], [0.75, 0.0], [0.8, -0.05], [0.9, -0.05]])

points = []
da = CurveProximityQueries.differentiate(a.f)
for v = 0:0.002:1
    pt = a(v)
    mvec = [0. 1; -1 0]*normalize(da(v))
    push!(points, Array(pt+mvec*0.13))
end
ahull = concave_hull(points)

plot(ahull)
plot!(a,
      linecolor=:black,
      linewidth=1.0,
      fill=(0, Gray(0.7)),
     )
plot!([0.282932, 0.412878], [-0.195185, -0.198932], linecolor = Gray(0.2), markershape = :circle, markersize = 4, markercolor = :black, markerstrokecolor = :black)
plot!(b,
      linewidth=2.0,
      linecolor=:black,
     )
annotate!(0.07, -0.2, text("\$\\mathcal{B}\$", 25))
annotate!(0.8, -0.2, text("\$\\Psi_{\\mathcal{I}}\$", 25))
annotate!(0.348, -0.147, text("\$\\Delta\$", 22))
savefig("out/fig2b.pdf")


text:
#!/usr/bin/env julia
__precompile__(true)

module RHEOSafm

# installed from Julia package repository
using RHEOS
using NLopt

# Base and stdlib imports
import DelimitedFiles: readdlm, writedlm

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

# AFM_JPKimport.jl
export importJPK

# AFM_preprocessing.jl
export contact_point

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

include("AFM_JPKimport.jl")
include("AFM_preprocessing.jl")

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

end # module


text:
test_path = mktempdir()
@testset "ThermalGen data misspecification" begin
    # See https://discourse.julialang.org/t/how-to-use-test-warn/15557/5 about testing for warning throwing
    warn_message = "The data doesn't include devices of type ThermalStandard, consider changing the device models"
    model = DeviceModel(ThermalStandard, ThermalStandardUnitCommitment)
    c_sys5_re_only = PSB.build_system(PSITestSystems, "c_sys5_re_only")
    op_problem = OperationsProblem(MockOperationProblem, DCPPowerModel, c_sys5_re_only)
    @test_logs (:info,) (:warn, warn_message) match_mode = :any mock_construct_device!(
        op_problem,
        model,
    )
end

################################### Unit Commitment tests ##################################
@testset "Thermal UC With DC - PF" begin
    bin_variable_names = [
        PSI.make_variable_name(PSI.ON, PSY.ThermalStandard),
        PSI.make_variable_name(PSI.START, PSY.ThermalStandard),
        PSI.make_variable_name(PSI.STOP, PSY.ThermalStandard),
    ]
    uc_constraint_names = [
        PSI.make_constraint_name(PSI.RAMP_UP, PSY.ThermalStandard),
        PSI.make_constraint_name(PSI.RAMP_DOWN, PSY.ThermalStandard),
        PSI.make_constraint_name(PSI.DURATION_UP, PSY.ThermalStandard),
        PSI.make_constraint_name(PSI.DURATION_DOWN, PSY.ThermalStandard),
    ]
    model = DeviceModel(ThermalStandard, ThermalStandardUnitCommitment)

    c_sys5_uc = PSB.build_system(PSITestSystems, "c_sys5_uc")
    op_problem = OperationsProblem(MockOperationProblem, DCPPowerModel, c_sys5_uc)
    mock_construct_device!(op_problem, model)
    moi_tests(op_problem, false, 480, 0, 480, 120, 120, true)
    psi_constraint_test(op_problem, uc_constraint_names)
    psi_checkbinvar_test(op_problem, bin_variable_names)
    psi_checkobjfun_test(op_problem, GAEVF)

    op_problem = OperationsProblem(
        MockOperationProblem,
        DCPPowerModel,
        c_sys5_uc;
        use_parameters = true,
    )
    mock_construct_device!(op_problem, model)
    moi_tests(op_problem, true, 480, 0, 480, 120, 120, true)
    psi_constraint_test(op_problem, uc_constraint_names)
    psi_checkbinvar_test(op_problem, bin_variable_names)
    psi_checkobjfun_test(op_problem, GAEVF)

    c_sys14 = PSB.build_system(PSITestSystems, "c_sys14")
    for p in [true, false]
        op_problem = OperationsProblem(
            MockOperationProblem,
            DCPPowerModel,
            c_sys14;
            use_parameters = p,
        )
        mock_construct_device!(op_problem, model)
        moi_tests(op_problem, p, 480, 0, 240, 120, 120, true)
        psi_checkbinvar_test(op_problem, bin_variable_names)
        psi_checkobjfun_test(op_problem, GQEVF)
    end
end

@testset "Thermal UC With AC - PF" begin
    bin_variable_names = [
        PSI.make_variable_name(PSI.ON, PSY.ThermalStandard),
        PSI.make_variable_name(PSI.START, PSY.ThermalStandard),
        PSI.make_variable_name(PSI.STOP, PSY.ThermalStandard),
    ]
    uc_constraint_names = [
        PSI.make_constraint_name(PSI.RAMP_UP, PSY.ThermalStandard),
        PSI.make_constraint_name(PSI.RAMP_DOWN, PSY.ThermalStandard),
        PSI.make_constraint_name(PSI.DURATION_UP, PSY.ThermalStandard),
        PSI.make_constraint_name(PSI.DURATION_DOWN, PSY.ThermalStandard),
    ]
    model = DeviceModel(ThermalStandard, ThermalStandardUnitCommitment)

    c_sys5_uc = PSB.build_system(PSITestSystems, "c_sys5_uc")
    op_problem = OperationsProblem(MockOperationProblem, ACPPowerModel, c_sys5_uc)
    mock_construct_device!(op_problem, model)
    moi_tests(op_problem, false, 600, 0, 600, 240, 120, true)
    psi_constraint_test(op_problem, uc_constraint_names)
    psi_checkbinvar_test(op_problem, bin_variable_names)
    psi_checkobjfun_test(op_problem, GAEVF)

    op_problem = OperationsProblem(
        MockOperationProblem,
        ACPPowerModel,
        c_sys5_uc;
        use_parameters = true,
    )
    mock_construct_device!(op_problem, model)
    moi_tests(op_problem, true, 600, 0, 600, 240, 120, true)
    psi_constraint_test(op_problem, uc_constraint_names)
    psi_checkbinvar_test(op_problem, bin_variable_names)
    psi_checkobjfun_test(op_problem, GAEVF)

    c_sys14 = PSB.build_system(PSITestSystems, "c_sys14")
    for p in [true, false]
        op_problem = OperationsProblem(
            MockOperationProblem,
            ACPPowerModel,
            c_sys14;
            use_parameters = p,
        )
        mock_construct_device!(op_problem, model)
        moi_tests(op_problem, p, 600, 0, 360, 240, 120, true)
        psi_checkbinvar_test(op_problem, bin_variable_names)
        psi_checkobjfun_test(op_problem, GQEVF)
    end
end

@testset "Thermal MultiStart UC With DC - PF" begin
    bin_variable_names = [
        PSI.make_variable_name(PSI.ON, PSY.ThermalMultiStart),
        PSI.make_variable_name(PSI.START, PSY.ThermalMultiStart),
        PSI.make_variable_name(PSI.STOP, PSY.ThermalMultiStart),
    ]
    uc_constraint_names = [
        PSI.make_constraint_name(PSI.RAMP_UP, PSY.ThermalMultiStart),
        PSI.make_constraint_name(PSI.RAMP_DOWN, PSY.ThermalMultiStart),
        PSI.make_constraint_name(PSI.DURATION_UP, PSY.ThermalMultiStart),
        PSI.make_constraint_name(PSI.DURATION_DOWN, PSY.ThermalMultiStart),
    ]
    model = DeviceModel(ThermalMultiStart, ThermalStandardUnitCommitment)

    c_sys5_uc = PSB.build_system(PSITestSystems, "c_sys5_pglib")
    for p in [true, false]
        op_problem = OperationsProblem(
            MockOperationProblem,
            DCPPowerModel,
            c_sys5_uc;
            use_parameters = p,
        )
        mock_construct_device!(op_problem, model)
        moi_tests(op_problem, p, 384, 0, 240, 48, 96, true)
        psi_constraint_test(op_problem, uc_constraint_names)
        psi_checkbinvar_test(op_problem, bin_variable_names)
        psi_checkobjfun_test(op_problem, GAEVF)
    end
end

@testset "Thermal MultiStart UC With AC - PF" begin
    bin_variable_names = [
        PSI.make_variable_name(PSI.ON, PSY.ThermalMultiStart),
        PSI.make_variable_name(PSI.START, PSY.ThermalMultiStart),
        PSI.make_variable_name(PSI.STOP, PSY.ThermalMultiStart),
    ]
    uc_constraint_names = [
        PSI.make_constraint_name(PSI.RAMP_UP, PSY.ThermalMultiStart),
        PSI.make_constraint_name(PSI.RAMP_DOWN, PSY.ThermalMultiStart),
        PSI.make_constraint_name(PSI.DURATION_UP, PSY.ThermalMultiStart),
        PSI.make_constraint_name(PSI.DURATION_DOWN, PSY.ThermalMultiStart),
    ]
    model = DeviceModel(ThermalMultiStart, ThermalStandardUnitCommitment)

    c_sys5_uc = PSB.build_system(PSITestSystems, "c_sys5_pglib")
    for p in [true, false]
        op_problem = OperationsProblem(
            MockOperationProblem,
            ACPPowerModel,
            c_sys5_uc;
            use_parameters = p,
        )
        mock_construct_device!(op_problem, model)
        moi_tests(op_problem, p, 432, 0, 288, 96, 96, true)
        psi_constraint_test(op_problem, uc_constraint_names)
        psi_checkbinvar_test(op_problem, bin_variable_names)
        psi_checkobjfun_test(op_problem, GAEVF)
    end
end

################################### Basic Unit Commitment tests ############################
@testset "Thermal Basic UC With DC - PF" begin
    bin_variable_names = [
        PSI.make_variable_name(PSI.ON, PSY.ThermalStandard),
        PSI.make_variable_name(PSI.START, PSY.ThermalStandard),
        PSI.make_variable_name(PSI.STOP, PSY.ThermalStandard),
    ]
    model = DeviceModel(ThermalStandard, ThermalBasicUnitCommitment)

    c_sys5_uc = PSB.build_system(PSITestSystems, "c_sys5_uc")
    op_problem = OperationsProblem(MockOperationProblem, DCPPowerModel, c_sys5_uc)
    mock_construct_device!(op_problem, model)
    moi_tests(op_problem, false, 480, 0, 240, 120, 120, true)
    psi_checkbinvar_test(op_problem, bin_variable_names)
    psi_checkobjfun_test(op_problem, GAEVF)

    op_problem = OperationsProblem(
        MockOperationProblem,
        DCPPowerModel,
        c_sys5_uc;
        use_parameters = true,
    )
    mock_construct_device!(op_problem, model)
    moi_tests(op_problem, true, 480, 0, 240, 120, 120, true)
    psi_checkbinvar_test(op_problem, bin_variable_names)
    psi_checkobjfun_test(op_problem, GAEVF)

    c_sys14 = PSB.build_system(PSITestSystems, "c_sys14")
    for p in [true, false]
        op_problem = OperationsProblem(
            MockOperationProblem,
            DCPPowerModel,
            c_sys14;
            use_parameters = p,
        )
        mock_construct_device!(op_problem, model)
        moi_tests(op_problem, p, 480, 0, 240, 120, 120, true)
        psi_checkbinvar_test(op_problem, bin_variable_names)
        psi_checkobjfun_test(op_problem, GQEVF)
    end
end

@testset "Thermal Basic UC With AC - PF" begin
    bin_variable_names = [
        PSI.make_variable_name(PSI.ON, PSY.ThermalStandard),
        PSI.make_variable_name(PSI.START, PSY.ThermalStandard),
        PSI.make_variable_name(PSI.STOP, PSY.ThermalStandard),
    ]
    model = DeviceModel(ThermalStandard, ThermalBasicUnitCommitment)

    c_sys5_uc = PSB.build_system(PSITestSystems, "c_sys5_uc")
    op_problem = OperationsProblem(MockOperationProblem, ACPPowerModel, c_sys5_uc)
    mock_construct_device!(op_problem, model)
    moi_tests(op_problem, false, 600, 0, 360, 240, 120, true)
    psi_checkbinvar_test(op_problem, bin_variable_names)
    psi_checkobjfun_test(op_problem, GAEVF)

    op_problem = OperationsProblem(
        MockOperationProblem,
        ACPPowerModel,
        c_sys5_uc;
        use_parameters = true,
    )
    mock_construct_device!(op_problem, model)
    moi_tests(op_problem, true, 600, 0, 360, 240, 120, true)
    psi_checkbinvar_test(op_problem, bin_variable_names)
    psi_checkobjfun_test(op_problem, GAEVF)

    c_sys14 = PSB.build_system(PSITestSystems, "c_sys14")
    for p in [true, false]
        op_problem = OperationsProblem(
            MockOperationProblem,
            ACPPowerModel,
            c_sys14;
            use_parameters = p,
        )
        mock_construct_device!(op_problem, model)
        moi_tests(op_problem, p, 600, 0, 360, 240, 120, true)
        psi_checkbinvar_test(op_problem, bin_variable_names)
        psi_checkobjfun_test(op_problem, GQEVF)
    end
end

@testset "Thermal MultiStart Basic UC With DC - PF" begin
    bin_variable_names = [
        PSI.make_variable_name(PSI.ON, PSY.ThermalMultiStart),
        PSI.make_variable_name(PSI.START, PSY.ThermalMultiStart),
        PSI.make_variable_name(PSI.STOP, PSY.ThermalMultiStart),
    ]
    model = DeviceModel(ThermalMultiStart, ThermalBasicUnitCommitment)

    c_sys5_uc = PSB.build_system(PSITestSystems, "c_sys5_pglib")
    for p in [true, false]
        op_problem = OperationsProblem(
            MockOperationProblem,
            DCPPowerModel,
            c_sys5_uc;
            use_parameters = p,
        )
        mock_construct_device!(op_problem, model)
        moi_tests(op_problem, p, 384, 0, 96, 48, 96, true)
        psi_checkbinvar_test(op_problem, bin_variable_names)
        psi_checkobjfun_test(op_problem, GAEVF)
    end
end

@testset "Thermal MultiStart Basic UC With AC - PF" begin
    bin_variable_names = [
        PSI.make_variable_name(PSI.ON, PSY.ThermalMultiStart),
        PSI.make_variable_name(PSI.START, PSY.ThermalMultiStart),
        PSI.make_variable_name(PSI.STOP, PSY.ThermalMultiStart),
    ]
    model = DeviceModel(ThermalMultiStart, ThermalBasicUnitCommitment)

    c_sys5_uc = PSB.build_system(PSITestSystems, "c_sys5_pglib")
    for p in [true, false]
        op_problem = OperationsProblem(
            MockOperationProblem,
            ACPPowerModel,
            c_sys5_uc;
            use_parameters = p,
        )
        mock_construct_device!(op_problem, model)
        moi_tests(op_problem, p, 432, 0, 144, 96, 96, true)
        psi_checkbinvar_test(op_problem, bin_variable_names)
        psi_checkobjfun_test(op_problem, GAEVF)
    end
end

################################### Basic Dispatch tests ###################################
@testset "Thermal Dispatch With DC - PF" begin
    model = DeviceModel(ThermalStandard, ThermalDispatch)
    c_sys5 = PSB.build_system(PSITestSystems, "c_sys5")
    for p in [true, false]
        op_problem = OperationsProblem(
            MockOperationProblem,
            DCPPowerModel,
            c_sys5;
            use_parameters = p,
        )
        mock_construct_device!(op_problem, model)
        moi_tests(op_problem, p, 120, 0, 120, 120, 0, false)
        psi_checkobjfun_test(op_problem, GAEVF)
    end

    c_sys14 = PSB.build_system(PSITestSystems, "c_sys14")
    for p in [true, false]
        op_problem = OperationsProblem(
            MockOperationProblem,
            DCPPowerModel,
            c_sys14;
            use_parameters = p,
        )
        mock_construct_device!(op_problem, model)
        moi_tests(op_problem, p, 120, 0, 120, 120, 0, false)
        psi_checkobjfun_test(op_problem, GQEVF)
    end
end

@testset "Thermal Dispatch With AC - PF" begin
    model = DeviceModel(ThermalStandard, ThermalDispatch)
    c_sys5 = PSB.build_system(PSITestSystems, "c_sys5")
    for p in [true, false]
        op_problem = OperationsProblem(
            MockOperationProblem,
            ACPPowerModel,
            c_sys5;
            use_parameters = p,
        )
        mock_construct_device!(op_problem, model)
        moi_tests(op_problem, p, 240, 0, 240, 240, 0, false)
        psi_checkobjfun_test(op_problem, GAEVF)
    end

    c_sys14 = PSB.build_system(PSITestSystems, "c_sys14")
    for p in [true, false]
        op_problem = OperationsProblem(
            MockOperationProblem,
            ACPPowerModel,
            c_sys14;
            use_parameters = p,
        )
        mock_construct_device!(op_problem, model)
        moi_tests(op_problem, p, 240, 0, 240, 240, 0, false)
        psi_checkobjfun_test(op_problem, GQEVF)
    end
end

# This Formulation is currently broken
@testset "ThermalMultiStart Dispatch With DC - PF" begin
    model = DeviceModel(ThermalMultiStart, ThermalDispatch)
    c_sys5 = PSB.build_system(PSITestSystems, "c_sys5_pglib")
    for p in [true, false]
        op_problem = OperationsProblem(
            MockOperationProblem,
            DCPPowerModel,
            c_sys5;
            use_parameters = p,
        )
        mock_construct_device!(op_problem, model)
        moi_tests(op_problem, p, 240, 0, 48, 48, 48, false)
        psi_checkobjfun_test(op_problem, GAEVF)
    end
end

@testset "ThermalMultiStart Dispatch With AC - PF" begin
    model = DeviceModel(ThermalMultiStart, ThermalDispatch)
    c_sys5 = PSB.build_system(PSITestSystems, "c_sys5_pglib")
    for p in [true, false]
        op_problem = OperationsProblem(
            MockOperationProblem,
            ACPPowerModel,
            c_sys5;
            use_parameters = p,
        )
        mock_construct_device!(op_problem, model)
        moi_tests(op_problem, p, 288, 0, 96, 96, 48, false)
        psi_checkobjfun_test(op_problem, GAEVF)
    end
end

################################### No Minimum Dispatch tests ##############################

@testset "Thermal Dispatch NoMin With DC - PF" begin
    model = DeviceModel(ThermalStandard, ThermalDispatchNoMin)
    c_sys5 = PSB.build_system(PSITestSystems, "c_sys5")
    for p in [true, false]
        op_problem = OperationsProblem(
            MockOperationProblem,
            DCPPowerModel,
            c_sys5;
            use_parameters = p,
        )
        mock_construct_device!(op_problem, model)
        moi_tests(op_problem, p, 120, 0, 120, 120, 0, false)
        moi_lbvalue_test(op_problem, :P_lb__ThermalStandard__RangeConstraint, 0.0)
        psi_checkobjfun_test(op_problem, GAEVF)
    end

    c_sys14 = PSB.build_system(PSITestSystems, "c_sys14")
    for p in [true, false]
        op_problem = OperationsProblem(
            MockOperationProblem,
            DCPPowerModel,
            c_sys14;
            use_parameters = p,
        )
        mock_construct_device!(op_problem, model)
        moi_tests(op_problem, p, 120, 0, 120, 120, 0, false)
        moi_lbvalue_test(op_problem, :P_lb__ThermalStandard__RangeConstraint, 0.0)
        psi_checkobjfun_test(op_problem, GQEVF)
    end
end

@testset "Thermal Dispatch NoMin With AC - PF" begin
    model = DeviceModel(ThermalStandard, ThermalDispatchNoMin)
    c_sys5 = PSB.build_system(PSITestSystems, "c_sys5")
    for p in [true, false]
        op_problem = OperationsProblem(
            MockOperationProblem,
            ACPPowerModel,
            c_sys5;
            use_parameters = p,
        )
        mock_construct_device!(op_problem, model)
        moi_tests(op_problem, p, 240, 0, 240, 240, 0, false)
        moi_lbvalue_test(op_problem, :P_lb__ThermalStandard__RangeConstraint, 0.0)
        psi_checkobjfun_test(op_problem, GAEVF)
    end

    c_sys14 = PSB.build_system(PSITestSystems, "c_sys14")
    for p in [true, false]
        op_problem = OperationsProblem(
            MockOperationProblem,
            ACPPowerModel,
            c_sys14;
            use_parameters = p,
        )
        mock_construct_device!(op_problem, model)
        moi_tests(op_problem, p, 240, 0, 240, 240, 0, false)
        moi_lbvalue_test(op_problem, :P_lb__ThermalStandard__RangeConstraint, 0.0)
        psi_checkobjfun_test(op_problem, GQEVF)
    end
end

# This Formulation is currently broken
#=
@testset "Thermal Dispatch NoMin With DC - PF" begin
    model = DeviceModel(ThermalMultiStart, ThermalDispatchNoMin)
    c_sys5 = PSB.build_system(PSITestSystems, "c_sys5_pglib")
    for p in [true, false]
        op_problem = OperationsProblem(
            MockOperationProblem,
            DCPPowerModel,
            c_sys5;
            use_parameters = p,
        )
        mock_construct_device!(op_problem, model)
        moi_tests(op_problem, p, 240, 0, 48, 48, 48, false)
        moi_lbvalue_test(op_problem, :P_lb__ThermalMultiStart__RangeConstraint, 0.0)
        psi_checkobjfun_test(op_problem, GAEVF)
    end
end

@testset "ThermalMultiStart Dispatch NoMin With AC - PF" begin
    model = DeviceModel(ThermalMultiStart, ThermalDispatchNoMin)
    c_sys5 = PSB.build_system(PSITestSystems, "c_sys5_pglib")
    for p in [true, false]
        op_problem = OperationsProblem(
            MockOperationProblem,
            ACPPowerModel,
            c_sys5;
            use_parameters = p,
        )
        mock_construct_device!(op_problem, model)
        moi_tests(op_problem, p, 288, 0, 96, 96, 48, false)
        moi_lbvalue_test(op_problem, :P_lb__ThermalMultiStart__RangeConstraint, 0.0)
        psi_checkobjfun_test(op_problem, GAEVF)
    end
end
=#
################################### Ramp Limited Testing ##################################
@testset "Thermal Ramp Limited Dispatch With DC - PF" begin
    constraint_names = [
        PSI.make_constraint_name(PSI.RAMP_UP, PSY.ThermalStandard),
        PSI.make_constraint_name(PSI.RAMP_DOWN, PSY.ThermalStandard),
    ]
    model = DeviceModel(ThermalStandard, ThermalRampLimited)
    c_sys5_uc = PSB.build_system(PSITestSystems, "c_sys5_uc")
    for p in [true, false]
        op_problem = OperationsProblem(
            MockOperationProblem,
            DCPPowerModel,
            c_sys5_uc;
            use_parameters = p,
        )
        mock_construct_device!(op_problem, model)
        moi_tests(op_problem, p, 120, 0, 216, 120, 0, false)
        psi_constraint_test(op_problem, constraint_names)
        psi_checkobjfun_test(op_problem, GAEVF)
    end

    c_sys14 = PSB.build_system(PSITestSystems, "c_sys14")
    for p in [true, false]
        op_problem = OperationsProblem(
            MockOperationProblem,
            DCPPowerModel,
            c_sys14;
            use_parameters = p,
        )
        mock_construct_device!(op_problem, model)
        moi_tests(op_problem, p, 120, 0, 120, 120, 0, false)
        psi_checkobjfun_test(op_problem, GQEVF)
    end
end

@testset "Thermal Ramp Limited Dispatch With AC - PF" begin
    constraint_names = [
        PSI.make_constraint_name(PSI.RAMP_UP, PSY.ThermalStandard),
        PSI.make_constraint_name(PSI.RAMP_DOWN, PSY.ThermalStandard),
    ]
    model = DeviceModel(ThermalStandard, ThermalRampLimited)
    c_sys5_uc = PSB.build_system(PSITestSystems, "c_sys5_uc")
    for p in [true, false]
        op_problem = OperationsProblem(
            MockOperationProblem,
            ACPPowerModel,
            c_sys5_uc;
            use_parameters = p,
        )
        mock_construct_device!(op_problem, model)
        moi_tests(op_problem, p, 240, 0, 336, 240, 0, false)
        psi_constraint_test(op_problem, constraint_names)
        psi_checkobjfun_test(op_problem, GAEVF)
    end

    c_sys14 = PSB.build_system(PSITestSystems, "c_sys14")
    for p in [true, false]
        op_problem = OperationsProblem(
            MockOperationProblem,
            ACPPowerModel,
            c_sys14;
            use_parameters = p,
        )
        mock_construct_device!(op_problem, model)
        moi_tests(op_problem, p, 240, 0, 240, 240, 0, false)
        psi_checkobjfun_test(op_problem, GQEVF)
    end
end

@testset "Thermal Ramp Limited Dispatch With DC - PF" begin
    constraint_names = [
        PSI.make_constraint_name(PSI.RAMP_UP, PSY.ThermalMultiStart),
        PSI.make_constraint_name(PSI.RAMP_DOWN, PSY.ThermalMultiStart),
    ]
    model = DeviceModel(ThermalMultiStart, ThermalRampLimited)
    c_sys5_uc = PSB.build_system(PSITestSystems, "c_sys5_pglib")
    for p in [true, false]
        op_problem = OperationsProblem(
            MockOperationProblem,
            DCPPowerModel,
            c_sys5_uc;
            use_parameters = p,
        )
        mock_construct_device!(op_problem, model)
        moi_tests(op_problem, p, 240, 0, 144, 48, 48, false)
        psi_constraint_test(op_problem, constraint_names)
        psi_checkobjfun_test(op_problem, GAEVF)
    end
end

@testset "Thermal Ramp Limited Dispatch With AC - PF" begin
    constraint_names = [
        PSI.make_constraint_name(PSI.RAMP_UP, PSY.ThermalMultiStart),
        PSI.make_constraint_name(PSI.RAMP_DOWN, PSY.ThermalMultiStart),
    ]
    model = DeviceModel(ThermalMultiStart, ThermalRampLimited)
    c_sys5_uc = PSB.build_system(PSITestSystems, "c_sys5_pglib")
    for p in [true, false]
        op_problem = OperationsProblem(
            MockOperationProblem,
            ACPPowerModel,
            c_sys5_uc;
            use_parameters = p,
        )
        mock_construct_device!(op_problem, model)
        moi_tests(op_problem, p, 288, 0, 192, 96, 48, false)
        psi_constraint_test(op_problem, constraint_names)
        psi_checkobjfun_test(op_problem, GAEVF)
    end
end

################################### ThermalMultiStart Testing ##################################

@testset "Thermal MultiStart with MultiStart UC and DC - PF" begin
    constraint_names = [
        PSI.make_constraint_name(PSI.ACTIVE_RANGE_IC, PSY.ThermalMultiStart),
        PSI.make_constraint_name(PSI.START_TYPE, PSY.ThermalMultiStart),
        PSI.make_constraint_name(PSI.MUST_RUN_LB, PSY.ThermalMultiStart),
        PSI.make_constraint_name(PSI.STARTUP_TIMELIMIT_WARM, PSY.ThermalMultiStart),
        PSI.make_constraint_name(PSI.STARTUP_TIMELIMIT_HOT, PSY.ThermalMultiStart),
        PSI.make_constraint_name(PSI.STARTUP_INITIAL_CONDITION_LB, PSY.ThermalMultiStart),
        PSI.make_constraint_name(PSI.STARTUP_INITIAL_CONDITION_UB, PSY.ThermalMultiStart),
    ]
    model = DeviceModel(PSY.ThermalMultiStart, PSI.ThermalMultiStartUnitCommitment)
    no_less_than = Dict(true => 334, false => 330)
    c_sys5_pglib = PSB.build_system(PSITestSystems, "c_sys5_pglib")
    for p in [true, false]
        op_problem = OperationsProblem(
            MockOperationProblem,
            DCPPowerModel,
            c_sys5_pglib;
            use_parameters = p,
        )
        mock_construct_device!(op_problem, model)
        moi_tests(op_problem, p, 528, 0, no_less_than[p], 60, 144, true)
        psi_constraint_test(op_problem, constraint_names)
        psi_checkobjfun_test(op_problem, GAEVF)
    end
end

@testset "Thermal MultiStart with MultiStart UC and AC - PF" begin
    constraint_names = [
        PSI.make_constraint_name(PSI.ACTIVE_RANGE_IC, PSY.ThermalMultiStart),
        PSI.make_constraint_name(PSI.START_TYPE, PSY.ThermalMultiStart),
        PSI.make_constraint_name(PSI.MUST_RUN_LB, PSY.ThermalMultiStart),
        PSI.make_constraint_name(PSI.STARTUP_TIMELIMIT_WARM, PSY.ThermalMultiStart),
        PSI.make_constraint_name(PSI.STARTUP_TIMELIMIT_HOT, PSY.ThermalMultiStart),
        PSI.make_constraint_name(PSI.STARTUP_INITIAL_CONDITION_LB, PSY.ThermalMultiStart),
        PSI.make_constraint_name(PSI.STARTUP_INITIAL_CONDITION_UB, PSY.ThermalMultiStart),
    ]
    model = DeviceModel(PSY.ThermalMultiStart, PSI.ThermalMultiStartUnitCommitment)
    no_less_than = Dict(true => 382, false => 378)
    c_sys5_pglib = PSB.build_system(PSITestSystems, "c_sys5_pglib")
    for p in [true, false]
        op_problem = OperationsProblem(
            MockOperationProblem,
            ACPPowerModel,
            c_sys5_pglib;
            use_parameters = p,
        )
        mock_construct_device!(op_problem, model)
        moi_tests(op_problem, p, 576, 0, no_less_than[p], 108, 144, true)
        psi_constraint_test(op_problem, constraint_names)
        psi_checkobjfun_test(op_problem, GAEVF)
    end
end

################################### Thermal Compact UC Testing ##################################
@testset "Thermal Standard with Compact UC and DC - PF" begin
    model = DeviceModel(PSY.ThermalStandard, PSI.ThermalCompactUnitCommitment)
    c_sys5 = PSB.build_system(PSITestSystems, "c_sys5")
    for p in [true, false]
        op_problem = OperationsProblem(
            MockOperationProblem,
            DCPPowerModel,
            c_sys5;
            use_parameters = p,
        )
        mock_construct_device!(op_problem, model)
        moi_tests(op_problem, p, 480, 0, 595, 0, 120, true)
        psi_checkobjfun_test(op_problem, GAEVF)
    end
end

@testset "Thermal MultiStart with Compact UC and DC - PF" begin
    model = DeviceModel(PSY.ThermalMultiStart, PSI.ThermalCompactUnitCommitment)
    c_sys5_pglib = PSB.build_system(PSITestSystems, "c_sys5_pglib")
    for p in [true, false]
        op_problem = OperationsProblem(
            MockOperationProblem,
            DCPPowerModel,
            c_sys5_pglib;
            use_parameters = p,
        )
        mock_construct_device!(op_problem, model)
        moi_tests(op_problem, p, 384, 0, 286, 0, 96, true)
        psi_checkobjfun_test(op_problem, GAEVF)
    end
end

@testset "Thermal Standard with Compact UC and AC - PF" begin
    model = DeviceModel(PSY.ThermalStandard, PSI.ThermalCompactUnitCommitment)
    c_sys5 = PSB.build_system(PSITestSystems, "c_sys5")
    for p in [true, false]
        op_problem = OperationsProblem(
            MockOperationProblem,
            ACPPowerModel,
            c_sys5;
            use_parameters = p,
        )
        mock_construct_device!(op_problem, model)
        moi_tests(op_problem, p, 600, 0, 715, 120, 120, true)
        psi_checkobjfun_test(op_problem, GAEVF)
    end
end

@testset "Thermal MultiStart with Compact UC and AC - PF" begin
    model = DeviceModel(PSY.ThermalMultiStart, PSI.ThermalCompactUnitCommitment)
    c_sys5_pglib = PSB.build_system(PSITestSystems, "c_sys5_pglib")
    for p in [true, false]
        op_problem = OperationsProblem(
            MockOperationProblem,
            ACPPowerModel,
            c_sys5_pglib;
            use_parameters = p,
        )
        mock_construct_device!(op_problem, model)
        moi_tests(op_problem, p, 432, 0, 334, 48, 96, true)
        psi_checkobjfun_test(op_problem, GAEVF)
    end
end

################################### Thermal Compact Dispatch Testing ##################################

@testset "Thermal Standard with Compact Dispatch and DC - PF" begin
    model = DeviceModel(PSY.ThermalStandard, PSI.ThermalCompactDispatch)
    c_sys5 = PSB.build_system(PSITestSystems, "c_sys5")
    for p in [true, false]
        op_problem = OperationsProblem(
            MockOperationProblem,
            DCPPowerModel,
            c_sys5;
            use_parameters = p,
        )
        mock_construct_device!(op_problem, model)
        moi_tests(op_problem, p, 120, 0, 168, 120, 0, false)
        psi_checkobjfun_test(op_problem, GAEVF)
    end
end

@testset "Thermal MultiStart with Compact Dispatch and DC - PF" begin
    model = DeviceModel(PSY.ThermalMultiStart, PSI.ThermalCompactDispatch)
    c_sys5_pglib = PSB.build_system(PSITestSystems, "c_sys5_pglib")
    for p in [true, false]
        op_problem = OperationsProblem(
            MockOperationProblem,
            DCPPowerModel,
            c_sys5_pglib;
            use_parameters = p,
        )
        mock_construct_device!(op_problem, model)
        moi_tests(op_problem, p, 240, 0, 144, 48, 48, false)
        psi_checkobjfun_test(op_problem, GAEVF)
    end
end

@testset "Thermal Standard with Compact Dispatch and AC - PF" begin
    model = DeviceModel(PSY.ThermalStandard, PSI.ThermalCompactDispatch)
    c_sys5 = PSB.build_system(PSITestSystems, "c_sys5")
    for p in [true, false]
        op_problem = OperationsProblem(
            MockOperationProblem,
            ACPPowerModel,
            c_sys5;
            use_parameters = p,
        )
        mock_construct_device!(op_problem, model)
        moi_tests(op_problem, p, 240, 0, 288, 240, 0, false)
        psi_checkobjfun_test(op_problem, GAEVF)
    end
end

@testset "Thermal MultiStart with Compact Dispatch and AC - PF" begin
    model = DeviceModel(PSY.ThermalMultiStart, PSI.ThermalCompactDispatch)
    no_less_than = Dict(true => 382, false => 378)
    c_sys5_pglib = PSB.build_system(PSITestSystems, "c_sys5_pglib")
    for p in [true, false]
        op_problem = OperationsProblem(
            MockOperationProblem,
            ACPPowerModel,
            c_sys5_pglib;
            use_parameters = p,
        )
        mock_construct_device!(op_problem, model)
        moi_tests(op_problem, p, 288, 0, 192, 96, 48, false)
        psi_checkobjfun_test(op_problem, GAEVF)
    end
end

############################# Model validation tests #######################################
@testset "Solving ED with CopperPlate for testing Ramping Constraints" begin
    ramp_test_sys = PSB.build_system(PSITestSystems, "c_ramp_test")
    template = OperationsProblemTemplate(CopperPlatePowerModel)
    set_device_model!(template, ThermalStandard, ThermalRampLimited)
    set_device_model!(template, PowerLoad, StaticPowerLoad)
    ED = OperationsProblem(
        EconomicDispatchProblem,
        template,
        ramp_test_sys;
        optimizer = Cbc_optimizer,
    )
    @test build!(ED; output_dir = mktempdir(cleanup = true)) == PSI.BuildStatus.BUILT
    moi_tests(ED, false, 10, 0, 20, 10, 5, false)
    res = solve!(ED)
    psi_checksolve_test(ED, [MOI.OPTIMAL], 11191.00)
end

# Testing Duration Constraints
@testset "Solving UC with CopperPlate for testing Duration Constraints" begin
    template = get_thermal_standard_uc_template()
    UC = OperationsProblem(
        UnitCommitmentProblem,
        template,
        PSB.build_system(PSITestSystems, "c_duration_test");
        optimizer = Cbc_optimizer,
        use_parameters = true,
    )
    @test build!(UC; output_dir = mktempdir(cleanup = true)) == PSI.BuildStatus.BUILT
    moi_tests(UC, true, 56, 0, 56, 14, 21, true)
    psi_checksolve_test(UC, [MOI.OPTIMAL], 8223.50)
end

## PWL linear Cost implementation test
@testset "Solving UC with CopperPlate testing Convex PWL" begin
    template = get_thermal_standard_uc_template()
    UC = OperationsProblem(
        UnitCommitmentProblem,
        template,
        PSB.build_system(PSITestSystems, "c_linear_pwl_test");
        optimizer = Cbc_optimizer,
        use_parameters = true,
    )
    @test build!(UC; output_dir = mktempdir(cleanup = true)) == PSI.BuildStatus.BUILT
    moi_tests(UC, true, 32, 0, 8, 4, 10, true)
    psi_checksolve_test(UC, [MOI.OPTIMAL], 9336.736919354838)
end

@testset "Solving UC with CopperPlate testing PWL-SOS2 implementation" begin
    template = get_thermal_standard_uc_template()
    UC = OperationsProblem(
        UnitCommitmentProblem,
        template,
        PSB.build_system(PSITestSystems, "c_sos_pwl_test");
        optimizer = Cbc_optimizer,
        use_parameters = true,
    )
    @test build!(UC; output_dir = mktempdir(cleanup = true)) == PSI.BuildStatus.BUILT
    moi_tests(UC, true, 32, 0, 8, 4, 14, true)
    psi_checksolve_test(UC, [MOI.OPTIMAL], 8500.89716, 10.0)
end

@testset "UC with MarketBid Cost in ThermalGenerators" begin
    template = get_thermal_standard_uc_template()
    set_device_model!(
        template,
        DeviceModel(ThermalMultiStart, ThermalMultiStartUnitCommitment),
    )
    UC = OperationsProblem(
        UnitCommitmentProblem,
        template,
        PSB.build_system(PSITestSystems, "c_market_bid_cost");
        optimizer = Cbc_optimizer,
        use_parameters = true,
    )
    @test build!(UC; output_dir = mktempdir(cleanup = true)) == PSI.BuildStatus.BUILT
    moi_tests(UC, true, 38, 0, 18, 8, 13, true)
end

@testset "Operation ModelThermalDispatchNoMin - and PWL Non Convex" begin
    c_sys5_pwl_ed_nonconvex = PSB.build_system(PSITestSystems, "c_sys5_pwl_ed_nonconvex")
    template = get_thermal_dispatch_template_network()
    set_device_model!(template, DeviceModel(ThermalStandard, ThermalDispatchNoMin))
    op_problem = OperationsProblem(
        MockOperationProblem,
        CopperPlatePowerModel,
        c_sys5_pwl_ed_nonconvex;
        use_parameters = true,
        export_pwl_vars = true,
    )
    @test_throws IS.InvalidValue mock_construct_device!(
        op_problem,
        DeviceModel(ThermalStandard, ThermalDispatchNoMin),
    )
end

#TODO: Add test for newer UC models
@testset "Solving UC Models with Linear Networks" begin
    c_sys5 = PSB.build_system(PSITestSystems, "c_sys5")
    c_sys5_dc = PSB.build_system(PSITestSystems, "c_sys5_dc")
    parameters_value = [true, false]
    systems = [c_sys5, c_sys5_dc]
    networks = [DCPPowerModel, NFAPowerModel, StandardPTDFModel, CopperPlatePowerModel]
    PTDF_ref = IdDict{System, PTDF}(c_sys5 => PTDF(c_sys5), c_sys5_dc => PTDF(c_sys5_dc))

    for net in networks, p in parameters_value, sys in systems
        template = get_thermal_dispatch_template_network(net)
        set_device_model!(template, ThermalStandard, ThermalStandardUnitCommitment)
        UC = OperationsProblem(
            template,
            sys;
            optimizer = GLPK_optimizer,
            use_parameters = p,
            PTDF = PTDF_ref[sys],
        )
        @test build!(UC; output_dir = mktempdir(cleanup = true)) == PSI.BuildStatus.BUILT
        psi_checksolve_test(UC, [MOI.OPTIMAL, MOI.LOCALLY_SOLVED], 340000, 100000)
    end
end


text:
#! /usr/bin/julia

# Rosetta Code, Heronian triangles

type IntegerTriangle{T<:Integer}
    a::T
    b::T
    c::T
    p::T
    σ::T
end

function IntegerTriangle{T<:Integer}(a::T, b::T, c::T)
    p = a + b + c
    s = div(p, 2)
    σ = isqrt(s*(s-a)*(s-b)*(s-c))
    (x, y, z) = sort([a, b, c])
    IntegerTriangle(x, y, z, p, σ)
end

function isprimheronian{T<:Integer}(a::T, b::T, c::T)
    p = a + b + c
    iseven(p) || return false
    gcd(a, b, c) == 1 || return false
    s = div(p, 2)
    t = s*(s-a)*(s-b)*(s-c)
    0 < t || return false
    σ = isqrt(t)
    σ^2 == t
end

slim = 200

ht = IntegerTriangle[]

for a in 1:slim, b in a:slim, c in b:slim
    isprimheronian(a, b, c) || continue
    push!(ht, IntegerTriangle(a, b, c))
end

sort!(ht, by=x->(x.σ, x.p, x.c))

print("The number of primitive Hernonian triangles having sides ≤ ")
println(slim, " is ", length(ht))

tlim = 10
tlim = min(tlim, length(ht))

println()
println("Tabulating the first (by σ, p, c) ", tlim, " of these:")
println("    a   b   c    σ    p")
for t in ht[1:tlim]
    println(@sprintf "%6d %3d %3d %4d %4d" t.a t.b t.c t.σ t.p)
end

tlim = 210
println()
println("Tabulating those having σ = ", tlim, ":")
println("    a   b   c    σ    p")
for t in ht
    t.σ == tlim || continue
    t.σ == tlim || break
    println(@sprintf "%6d %3d %3d %4d %4d" t.a t.b t.c t.σ t.p)
end

max_stars_repo_name: MichaeLeroy/rosetta-code

max_issues_repo_name: MichaeLeroy/rosetta-code

max_forks_repo_name: MichaeLeroy/rosetta-code


text:
"""
    reliabilityplot(probabilities, frequencies; deviation=true, kwargs...)
    reliabilityplot(probabilities, frequencies, low_high; deviation=true, kwargs...)
    reliabilityplot(probabilities, frequencies, low, high; deviation=true, kwargs...)

Plot a reliability diagram of the observed `frequencies` versus the predicted
`probabilities`, optionally with consistency bars ranging from `low` to `high` for
each probability.

If `deviation` is `true` (default), the difference `frequencies - probabilities` is
plotted versus `probabilities`. This can be helpful for inspecting how closely the
observed frequencies match the predicted probabilities.

# Examples

```jldoctest
julia> probabilities = sort(rand(10));

julia> frequencies = rand(10);

julia> reliabilityplot(probabilities, frequencies);

julia> # plot frequencies instead of deviations
       reliabilityplot(probabilities, frequencies; deviation=false);

julia> # use different colors for markers
       reliabilityplot(probabilities, frequencies; markercolor=:red);

julia> # random consistency bars
       low = max.(0, probabilities .- rand.());

julia> high = min.(1, probabilities .+ rand.());

julia> # plot pre-computed consistency bars
       reliabilityplot(probabilities, frequencies, low, high);

julia> # alternatively consistency bars can be specified as vector of tuples
       reliabilityplot(probabilities, frequencies, map(tuple, low, high));
```

# References

Bröcker, J. and Smith, L.A. (2007). *Increasing the reliability of reliability diagrams*.
Weather and forecasting, 22(3), pp. 651-661.
"""
RecipesBase.@userplot struct ReliabilityPlot{X<:AbstractVector,Y<:AbstractVector,T}
    x::X
    y::Y
    low_high::T
end

"""
    reliabilityplot(
        probabilities::AbstractVector{<:Real},
        outcomes::AbstractVector{Bool};
        binning=EqualMass(),
        consistencybars=ConsistencyBars(),
        kwargs...,
    )

Plot a reliability diagram of the observed frequencies versus the mean probabilities
in different clusters together with a set of consistency bars.

Optionally, the binning algorithm that is used to form the clusters and the consistency
bars can be configured with the `binning` and `consistencybars` keyword arguments. If
`consistencybars === nothing`, then no consistency bars are computed.

# Examples

```jldoctest
julia> probabilities = rand(10);

julia> outcomes = rand(Bool, 10);

julia> # default binning and consistency bars
       reliabilityplot(probabilities, outcomes);

julia> # custom options: without consistency bars
       reliabilityplot(probabilities, outcomes; consistencybars=nothing);
```

See also: [`EqualMass`](@ref), [`EqualSize`](@ref), [`ConsistencyBars`](@ref)
"""
reliabilityplot(::AbstractVector{<:Real}, ::AbstractVector{Bool}; kwargs...)

RecipesBase.@recipe function f(plot::ReliabilityPlot)
    # extract x and y values
    x = plot.x
    y = plot.y

    # compute mean probabilities and frequencies and consistency bars
    probs, freqs_nodeviation, low_high = if y isa AbstractVector{Bool}
        binning = get(plotattributes, :binning, EqualMass())
        consistencybars = get(plotattributes, :consistencybars, ConsistencyBars())
        means_and_bars(x, y; binning=binning, consistencybars=consistencybars)
    else
        x, y, plot.low_high
    end
    deviation::Bool = get(plotattributes, :deviation, true)
    freqs = deviation ? freqs_nodeviation - probs : freqs_nodeviation - zero(probs)

    # default attributes
    seriestype --> :line
    markershape --> :circle

    # add consistency bars, taking into account the `deviation` keyword argument
    if low_high !== nothing
        low = map(low_high, freqs_nodeviation) do (low, _), freq
            return freq - low
        end
        high = map(low_high, freqs_nodeviation) do (_, high), freq
            return high - freq
        end
        yerror := (low, high)
    end

    return probs, freqs
end


text:
module Example

export
    Foo

src = [
    "foo",
]

for s = src
    include("$(s).jl")
end

end # module


text:
################################################################################
# Types of Empirical Potentials
################################################################################

include("../empirical/lj.jl")
include("../empirical/bm.jl")
include("../empirical/coulomb.jl")
include("../empirical/zbl.jl")
include("../empirical/morse.jl")

export LennardJones, BornMayer, Coulomb, ZBL, Morse

################################################################################
# InteratomicPotentials API implementations for empirical potentials
################################################################################

force(r::SVector{3}, p::EmpiricalPotential) = force(norm(r), r, p)

function energy_and_force(s::AbstractSystem, p::EmpiricalPotential)
    nnlist = neighborlist(s, get_rcutoff(p))

    e = 0.0
    f = fill(SVector{3}(zeros(3)), length(s))
    for ii in 1:length(s)
        for (jj, r, R) in zip(nnlist.j[ii], nnlist.r[ii], nnlist.R[ii])
            if (ismissing(get_species(p)) || (atomic_symbol(s, ii), atomic_symbol(s, jj)) ⊆ get_species(p))
                e += potential_energy(R, p)
                fo = force(R, r, p)
                f[ii] = f[ii] + fo
                f[jj] = f[jj] - fo
            end
        end
    end
    (; e=e * ENERGY_UNIT, f=f * FORCE_UNIT)
end

function virial_stress(s::AbstractSystem, p::EmpiricalPotential)
    nnlist = neighborlist(s, get_rcutoff(p))

    v = @SVector zeros(6)
    for ii in 1:length(s)
        for (r, R) in zip(nnlist.r[ii], nnlist.R[ii])
            vi = r * force(R, r, p)'
            v += @SVector [vi[1, 1], vi[2, 2], vi[3, 3], vi[3, 2], vi[3, 1], vi[2, 1]]
        end
    end
    v * ENERGY_UNIT
end


text:
using GenomicAnnotations
using BioSequences
using Test

@testset "GenomicAnnotations" begin
    @testset "GenBank parsing" begin
        s = "     gene            1..1"
        @test GenBank.parseposition(s) == (:gene, Locus(1:1, '+'))
        s = "     gene            complement(order(3300..4037,4047..4052))"
        @test GenBank.parseposition(s) == (:gene, Locus(3300:4052, '-', true, true, UnitRange{Int}[3300:4037, 4047:4052], false))
        chrs = collect(open(GenBank.Reader, "example.gbk"))
        @test length(chrs) == 2
        @test chrs[2].name == "plasmid1"
    end

    chr = collect(open(GenBank.Reader, "example.gbk"))[1]

    @testset "GFF parsing" begin
        open(GFF.Writer, "example.gff") do w
            write(w, chr)
        end
        gff = collect(open(GFF.Reader, "example.gff"))[1]
        @test begin
            gbkbuf = IOBuffer()
            gffbuf = IOBuffer()
            print(gbkbuf, chr.genes[2:4])
            print(gffbuf, gff.genes[2:4])
            String(take!(gbkbuf)) == String(take!(gffbuf))
        end
    end

    @testset "Extended methods" begin
        @test length(chr.genes[1]) == length(sequence(chr.genes[1]))
        @test length(chr.sequence) == length(sequence(chr))
    end

    @testset "Gene properties" begin
        @test length(propertynames(chr.genes[1])) == 12
        @test chr.genes[2].locus_tag == "tag01"
        @test (chr.genes[2].locus_tag = "tag01") == "tag01"
        @test begin
            chr.genes[2].test = "New column"
            chr.genes[2].test == "New column" && all(ismissing, chr.genes[[1,3,4,5,6]].test)
        end
        @test get(chr.genes[1], :locus_tag, "") == ""
        @test get(chr.genes[2], :locus_tag, "") == "tag01"
        @test get(chr.genes[2], :locustag, "") == ""
        @test begin
            GenomicAnnotations.pushproperty!(chr.genes[2], :db_xref, "GI:123")
            chr.genes[2].db_xref == ["GI:1293614", "GI:123"]
        end
        @test GenomicAnnotations.vectorise(Union{Missing,Int}[1,1,1]) == [[1],[1],[1]]
    end

    @testset "Iteration" begin
        @test length([g.locus_tag for g in chr.genes]) == 7
    end

    @testset "Adding/removing genes" begin
        addgene!(chr, :CDS, Locus(300:390, '+'), locus_tag = "tag04")
        @test chr.genes[end].locus_tag == "tag04"
        delete!(chr.genes[end])
        @test chr.genes[end].locus_tag == "reg01"
    end

    @testset "@genes" begin
        @test length(union(@genes(chr, CDS), @genes(chr, !CDS))) == length(chr.genes)
        @test length(intersect(@genes(chr, CDS), @genes(chr, !CDS))) == 0
        @test length(union(@genes(chr, gene), @genes(chr, !gene))) == length(chr.genes)
        @test length(intersect(@genes(chr, gene), @genes(chr, !gene))) == 0
        @test length(@genes(chr)) == length(chr.genes)
        @test @genes(chr, feature(gene) == $:CDS) == chr.genes[[2,4,6]]
        @test @genes(chr, feature(gene) == $:CDS) == @genes(chr, CDS)
        @test @genes(chr, iscomplement(gene)) == chr.genes[[5,6,7]]
        @test @genes(chr, feature(gene) == $:CDS, !iscomplement(gene)) == chr.genes[[2,4]]
        @test @genes(chr, length(gene) < 300)[1] == chr.genes[2]
        @test length(@genes(chr, get(gene, :locus_tag, "") == "")) == 3
        gene = chr.genes[3]
        @test @genes(chr, gene == $gene)[1] == chr.genes[3]
        d = Dict(:a => "tag01")
        @test @genes(chr, :locus_tag == d[$:a]) == @genes(chr, :locus_tag == "tag01")
    end

    @testset "Broadcast" begin
        # Broadcasted assignment on existing property
        chr.genes[3:4].gene .= "AXL2P"
        @test all(chr.genes[3:4].gene .== "AXL2P")
        # Broadcasted assignment on previously missing property
        chr.genes[3:4].newproperty .= true
        @test all(chr.genes[3:4].newproperty .== true)
        # Broadcasted assignment with @genes
        @genes(chr, gene).newproperty .= false
        @test all(chr.genes[[3,5]].newproperty .== false)
    end

    @testset "Locus" begin
        loc = Locus(1:1, '.', true, true, UnitRange{Int}[], false)
        @test Locus() == loc
        @test locus(chr.genes[2]) < locus(chr.genes[4])
        @test locus(chr.genes[2]) == locus(chr.genes[2])
        @test iscomplement(chr.genes[2]) == false
        @test iscomplement(chr.genes[5]) == true
    end

    seq = dna"atgtccatatacaacggtatctccacctcaggtttagatctcaacaacggaaccattgccgacatgagacagttaggtatcgtcgagagttacaagctaaaacgagcagtagtcagctctgcatctgaagccgctgaagttctactaagggtggataacatcatccgtgcaagaccaagaaccgccaatagacaacatatgtaa"
    @test sequence(chr.genes[2]) == seq
    @test length(chr.genes[2]) == length(seq)

    @testset "Empty Record" begin
        chr = GenBank.Record()
        @test chr.name == ""
        @test chr.sequence == dna""
        @test chr.header == ""
        @test chr.genes == Gene[]
        @test names(chr.genedata) == []
    end
end


text:
function HansenThirdBodyLinear(arg0::jint, arg1::jint)
    return HansenThirdBodyLinear((jint, jint), arg0, arg1)
end

function compute_init_values(obj::HansenThirdBodyLinear, arg0::jdouble, arg1::jdouble, arg2::jdouble)
    return jcall(obj, "computeInitValues", void, (jdouble, jdouble, jdouble), arg0, arg1, arg2)
end

function get_derivative(obj::HansenThirdBodyLinear, arg0::jint, arg1::jdouble)
    return jcall(obj, "getDerivative", jdouble, (jint, jdouble), arg0, arg1)
end

function get_value(obj::HansenThirdBodyLinear, arg0::jint, arg1::jdouble)
    return jcall(obj, "getValue", jdouble, (jint, jdouble), arg0, arg1)
end


text:
import AbstractAlgebra

include("interpolation.jl")

# just for convenience
AA = AbstractAlgebra

mutable struct polynomial_interpolation
    groebner
    polyring
    xs
    ys
end

function eat_point(P::polynomial_interpolation, p)

    x, y = p

    push!(P.xs, x)
    push!(P.ys, y)

    if length(P.xs) == 40
        f = interpolate_polynomial(P.polyring, P.xs, P.ys)
        return f
    end

    return P.polyring(0)
end

function update(P::polynomial_interpolation, f)
    P.xs = []
    P.ys = []
end

mutable struct rational_interpolation
    polynomial
    predicted_degrees
    init
    power
    current
end

function generate_point(R::rational_interpolation)

    base = R.polynomial.groebner.base
    R.current = R.current * R.init
    R.power += 1

    dA, dB = R.predicted_degrees

    return R.current .* [base(2), base(3)] .* (base(5))^(dA - dB)
end

function interpolate_coeff(R::rational_interpolation)

    # power -> interpolated coeff
    polyring = R.polynomial.polyring

    A, B = polyring(0), polyring(0)

    while iszero(A) && iszero(B)
        nevals = sum(R.predicted_degrees) + 2
        xs = [ generate_point(R) for i in 1:nevals ]
        ys = [ evaluate(R.polynomial.groebner, xs[i]) for i in 1:nevals ]
        xs = [ x[1] for x in xs ]
        Q = interpolate_rational_function(R.polynomial.polyring, xs, ys)

        Ai, Bi = map(leading_coefficient, [numerator(Q), denominator(Q)])
        println(Q)

        A = eat_point(R.polynomial, Ai)
        B = eat_point(R.polynomial, Bi)
    end

    return A // B
end

function update(R::rational_interpolation, f)
    R.power = 0
    R.current = R.init^R.power
end

mutable struct groebner_evaluation
    f
    index
    base
end

function update(G::groebner_evaluation, f)
    G.index += 1
end

function evaluate(G::groebner_evaluation, p)
    return G.f(p...)[G.index]
end

FF = AA.GF(2^31 - 1)

# emulates groebner basis coeffs computation
f = (x, y) -> FF.([ (x^2 + y) // (x + 11), (x*y + 1) // (8x^3 + 13) ])

function predict_degrees(f)
    #return [(2, 1), (2, 3)]
    return (2, 1)
end

function interpolate(f)
    # predict shape (i.e positions of coeffs) of f
    predicted_degrees = predict_degrees(f)
    ncoeffs = 2

    base = FF

    interpolation_ring, = AA.PolynomialRing(base, "x")

    groebner = groebner_evaluation(f, 1, base)

    polynomial = polynomial_interpolation(groebner, interpolation_ring, [], [])

    rational = rational_interpolation(polynomial, predicted_degrees, base(2), 0, base(2)^0)

    i = 0
    while i < ncoeffs
        c = interpolate_coeff(rational)
        @info "interpolated " c
        i += 1
    end

end

interpolate(f)


text:
facts("Lemma") do
    context("constructor") do
        lemma = Lemma("fruit n 3 3 @ ~ + 3 3 13134947 04612722 07294550", 123)

        @fact lemma.id --> 123
        @fact lemma.word --> "fruit"
        @fact lemma.pos --> 'n'
        @fact lemma.pointer_syms --> ["@", "~", "+"]
        @fact lemma.tagsense_count --> 3
        @fact lemma.synset_offsets --> [13134947, 4612722, 7294550]
    end

    context("show") do
        lemma = Lemma("fruit n 3 3 @ ~ + 3 3 13134947 04612722 07294550", 123)

        @fact io_to_string(io -> show(io, lemma)) --> "fruit.n"
    end
end


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

struct VersionApi <: SwaggerApi
    client::Swagger.Client
end

function _swaggerinternal_getCodeVersion(_api::VersionApi; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", IoK8sApimachineryPkgVersionInfo, "/version/", ["BearerToken"])
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

"""

get the code version
Return: IoK8sApimachineryPkgVersionInfo
"""
function getCodeVersion(_api::VersionApi; _mediaType=nothing)
    _ctx = _swaggerinternal_getCodeVersion(_api; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function getCodeVersion(_api::VersionApi, response_stream::Channel; _mediaType=nothing)
    _ctx = _swaggerinternal_getCodeVersion(_api; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

export getCodeVersion