text:
"""
    _parse_grr(gpa::SBML.GeneProductAssociation)::GeneAssociation

Parse `SBML.GeneProductAssociation` structure to the simpler GeneAssociation.
The input must be (implicitly) in a positive DNF.
"""
function _parse_grr(gpa::SBML.GeneProductAssociation)::GeneAssociation
    parse_ref(x) =
        typeof(x) == SBML.GPARef ? [x.gene_product] :
        begin
            @_models_log @warn "Could not parse a part of gene association, ignoring: $x"
            String[]
        end
    parse_and(x) =
        typeof(x) == SBML.GPAAnd ? vcat([parse_and(i) for i in x.terms]...) : parse_ref(x)
    parse_or(x) =
        typeof(x) == SBML.GPAOr ? vcat([parse_or(i) for i in x.terms]...) : [parse_and(x)]
    return parse_or(gpa)
end

"""
    _unparse_grr(
        ::Type{SBML.GeneProductAssociation},
        x::GeneAssociation,
    )::SBML.GeneAssociation

Convert a GeneAssociation to the corresponding `SBML.jl` structure.
"""
function _unparse_grr(
    ::Type{SBML.GeneProductAssociation},
    x::GeneAssociation,
)::SBML.GeneAssociation
    SBML.GPAOr([SBML.GPAAnd([SBML.GPARef(j) for j in i]) for i in x])
end

"""
    _parse_grr(s::String)::GeneAssociation

Parse a DNF gene association rule in format `(YIL010W and YLR043C) or (YIL010W
and YGR209C)` to `GeneAssociation. Also accepts `OR`, `|`, `||`, `AND`, `&`,
and `&&`.

# Example
```
julia> _parse_grr("(YIL010W and YLR043C) or (YIL010W and YGR209C)")
2-element Array{Array{String,1},1}:
 ["YIL010W", "YLR043C"]
 ["YIL010W", "YGR209C"]
```
"""
_parse_grr(s::String)::Maybe{GeneAssociation} = _maybemap(_parse_grr, _parse_grr_to_sbml(s))

"""
    _parse_grr_to_sbml(str::String)::Maybe{SBML.GeneProductAssociation}

Internal helper for parsing the string GRRs into SBML data structures. More
general than [`_parse_grr`](@ref).
"""
function _parse_grr_to_sbml(str::String)::Maybe{SBML.GeneProductAssociation}
    s = str
    toks = String[]
    m = Nothing
    while !isnothing(begin
        m = match(r"( +|[a-zA-Z0-9]+|[^ a-zA-Z0-9()]+|[(]|[)])(.*)", s)
    end)
        tok = strip(m.captures[1])
        !isempty(tok) && push!(toks, tok)
        s = m.captures[2]
    end

    fail() = throw(DomainError(str, "Could not parse GRR"))

    # shunting yard
    ops = Symbol[]
    vals = SBML.GeneProductAssociation[]
    fold(sym, op) =
        while !isempty(ops) && last(ops) == sym
            r = pop!(vals)
            l = pop!(vals)
            pop!(ops)
            push!(vals, op([l, r]))
        end
    for tok in toks
        if tok in ["and", "AND", "&", "&&"]
            push!(ops, :and)
        elseif tok in ["or", "OR", "|", "||"]
            fold(:and, SBML.GPAAnd)
            push!(ops, :or)
        elseif tok == "("
            push!(ops, :paren)
        elseif tok == ")"
            fold(:and, SBML.GPAAnd)
            fold(:or, SBML.GPAOr)
            if isempty(ops) || last(ops) != :paren
                fail()
            else
                pop!(ops)
            end
        else
            push!(vals, SBML.GPARef(tok))
        end
    end

    fold(:and, SBML.GPAAnd)
    fold(:or, SBML.GPAOr)

    if !isempty(ops) || length(vals) > 1
        fail()
    end

    if isempty(vals)
        nothing
    else
        first(vals)
    end
end

"""
    unparse_grr(grr::Vector{Vector{Gene}}

Converts a nested string gene reaction array  back into a gene reaction rule
string.

# Example
```
julia> _unparse_grr(String, [["YIL010W", "YLR043C"], ["YIL010W", "YGR209C"]])
"(YIL010W and YLR043C) or (YIL010W and YGR209C)"
```
"""
function _unparse_grr(::Type{String}, grr::GeneAssociation)::String
    grr_strings = String[]
    for gr in grr
        push!(grr_strings, "(" * join([g for g in gr], " and ") * ")")
    end
    grr_string = join(grr_strings, " or ")
    return grr_string
end


text:
export current_context,
    with_context,
    Limited,
    n_dropped,
    TAttrVal,
    DynamicAttrs,
    StaticAttrs,
    Resource,
    InstrumentationInfo

using UUIDs: uuid4

const CONTEXT_KEY = :OPEN_TELEMETRY_CONTEXT

# ??? performace issue with NamedTuple
struct Context{T<:NamedTuple}
    kv::T
end

Context() = Context(NamedTuple())

Base.merge(c1::Context, c2::Context) = Context(merge(c1.kv, c2.kv))

create_key(s) = Symbol(s, '-', uuid4())

for f in (:getindex, :haskey, :get)
    @eval Base.$f(ctx::Context, args...) = $f(ctx.kv, args...)
end

with_context(f; kv...) = with_context(f, current_context(); kv...)

"""
    with_context(f, [context]; kv...)

Run function `f` in the `context`. If extra `kv` pairs are provided, they will be merged with the `context` to form a new context. When `context` is not provided, the [`current_context`](@ref) will be used.
"""
function with_context(f, ctx::Context; kw...)
    task_local_storage(CONTEXT_KEY, merge(ctx, Context(values(kw)))) do
        f()
    end
end

"""
Return the `Context` associated with the caller's current execution unit.
"""
current_context() = get!(task_local_storage(), CONTEXT_KEY, Context())

# !!! intentional type piracy
function Base.schedule(t::Task)
    if isnothing(t.storage)
        t.storage = IdDict()
    end
    t.storage[CONTEXT_KEY] = current_context()
    Base.enq_work(t)
end

#####

"""
    Limited(container; limit=32)

Create a container wrapper with limited elements. It supports the following common containers:

  - `AbstractDict`. If the length of the dict contains more elements than `limit`. Then `pop!` will be called repeatedly until the length is equal to `limit`. Further new key-value pair insertions will be ignored. If the key already exists in the dict, then the corresponding value will always be updated. `n_dropped` will return the number of dropped insertions.
  - `AbstractVector`. Similar to `AbstractDict`.

The following methods from `Base` are defined on `Limited` which are then forwarded to the inner `container`. Feel free to create a PR if you find any method you need is missing:

  - `Base.getindex`
  - `Base.setindex!`
  - `Base.iterate`
  - `Base.length`
  - `Base.haskey`
  - `Base.push!`. Only defined on containers of `AbstractVector`.
"""
struct Limited{T}
    xs::T
    limit::Int
    n_dropped::Ref{Int}
end

function Limited(xs::Union{Dict,AbstractVector}; limit = 32)
    n_dropped = Ref(0)
    if length(xs) > limit
        n_dropped[] = length(xs) - limit
        for _ in 1:(length(xs)-limit)
            pop!(xs)
        end
    end
    Limited(xs, limit, n_dropped)
end

Base.getindex(x::Limited, args...) = getindex(x.xs, args...)
Base.haskey(x::Limited, k) = haskey(x.xs, k)
Base.length(x::Limited) = length(x.xs)
Base.iterate(x::Limited, args...) = iterate(x.xs, args...)
Base.pairs(A::Limited) = pairs(A.xs)

"""
Return the total number of dropped elements since creation.
"""
n_dropped(x::Limited) = x.n_dropped[]

function Base.setindex!(x::Limited{<:AbstractDict}, v, k)
    if haskey(x.xs, k)
        setindex!(x.xs, v, k)
    elseif length(x.xs) >= x.limit
        @warn "limit [$(x.limit)] exceeded, dropped."
        x.n_dropped[] += 1
    else
        setindex!(x.xs, v, k)
    end
end

Base.setindex!(x::Limited{<:AbstractVector}, v, k) = setindex!(x.xs, v, k)

function Base.push!(x::Limited{<:AbstractVector}, v)
    if length(x.xs) >= x.limit
        @warn "limit [$(x.limit)] exceeded, dropped."
        x.n_dropped[] += 1
    else
        push!(x.xs, v)
    end
end

#####

_truncate(x, limit) = x

function _truncate(s::String, limit)
    if length(s) <= limit
        s
    else
        @warn "value is truncated to length $limit"
        s[1:min(end, limit)]
    end
end

function _truncate(xs::Vector{String}, limit)
    is_truncated = false
    for (i, s) in enumerate(xs)
        if length(s) > limit
            xs[i] = s[1:limit]
            is_truncated = true
        end
    end
    is_truncated && @warn "some elements in the value are truncated."
    xs
end

const TAttrVal = Union{
    String,
    Bool,
    Int64,
    Int128,
    Float64,
    Vector{String},
    Vector{Bool},
    Vector{Int64},
    Vector{Float64},
}

is_valid_attr_val(nt::NamedTuple) = is_valid_attr_val(values(nt))
is_valid_attr_val(t::Tuple) = is_valid_attr_val(first(t)) && is_valid_attr_val(Base.tail(t))
is_valid_attr_val(x::TAttrVal) = true
is_valid_attr_val(t::Tuple{}) = true
is_valid_attr_val(x) = false

"""
    StaticAttrs(attrs::NamedTuple; value_length_limit=nothing)
    StaticAttrs(attrs::Pair{String, TAttrVal}...; value_length_limit=nothing)
    StaticAttrs(attrs::Pair{Symbol, TAttrVal}...; value_length_limit=nothing)

Here we use the `NamedTuple` internally to efficiently represent the immutable version of the [`Attributes`](https://github.com/open-telemetry/opentelemetry-specification/blob/main/specification/common/common.md#attributes) specification. If the `value_length_limit` is set to a positive int, the value of `String` or each element in a value of `Vector{String}` will be truncated to a maximum length of `value_length_limit`. By default we do not do the truncation.

See also [`DynamicAttrs`](@ref).
"""
struct StaticAttrs{T<:NamedTuple}
    attrs::T
    function StaticAttrs(attrs::T; value_length_limit = nothing) where {T<:NamedTuple}
        if is_valid_attr_val(attrs)
            if !isnothing(value_length_limit)
                attrs = map(x -> _truncate(x, value_length_limit), attrs)
            end
            new{T}(attrs)
        else
            throw(ArgumentError("Invalid value type found!"))
        end
    end
end

function Base.show(io::IO, s::StaticAttrs)
    print(io, "StaticAttrs")
    if isempty(s.attrs)
        print(io, "()")
    else
        print(io, s.attrs)
    end
end

StaticAttrs(; kw...) = StaticAttrs(NamedTuple(); kw...)

n_dropped(a::StaticAttrs) = 0

Base.keys(A::StaticAttrs) = keys(A.attrs)
Base.getindex(A::StaticAttrs, k::String) = getindex(A, Symbol(k))
Base.getindex(A::StaticAttrs, k::Symbol) = getindex(A.attrs, k)
Base.getindex(A::StaticAttrs, k::Tuple{Vararg{Symbol}}) =
    StaticAttrs(NamedTuple{k}(A.attrs))
Base.length(A::StaticAttrs) = length(A.attrs)
Base.iterate(A::StaticAttrs, args...) = iterate(A.attrs, args...)
Base.pairs(A::StaticAttrs) = pairs(A.attrs)

# ??? a more efficient approach?
Base.sort(A::StaticAttrs) = A[keys(A)|>collect|>sort|>Tuple]

"""
    DynamicAttrs(attrs::Dict{String, TAttrVal};count_limit=128, value_length_limit=nothing)
    DynamicAttrs(attrs::Pair{String, TAttrVal}...;count_limit=128, value_length_limit=nothing)

Here we use a `Dict` internally to represent the mutable version of the [`Attributes`](https://github.com/open-telemetry/opentelemetry-specification/blob/main/specification/common/common.md#attributes) specification. If the `value_length_limit` is set to a positive int, the value of `String` or each element in a value of `Vector{String}` will be truncated to a maximum length of `value_length_limit`. By default we do not do the truncation. When adding new pairs into it, if the number of attributes exceeds the `count_limit`, it will be dropped. You can get the number of dropped pairs via `n_dropped`.
"""
struct DynamicAttrs
    attrs::Limited{Dict{String,TAttrVal}}
    value_length_limit::Union{Int,Nothing}
    function DynamicAttrs(
        d::Dict{String,TAttrVal};
        count_limit = 128,
        value_length_limit = nothing,
    )
        if !isnothing(value_length_limit)
            for k in keys(d)
                d[k] = _truncate(d[k], value_length_limit)
            end
        end
        new(Limited(d; limit = count_limit), value_length_limit)
    end
end

DynamicAttrs(attrs::Pair{String}...; kw...) =
    DynamicAttrs(Dict{String,TAttrVal}(attrs...); kw...)

Base.getindex(A::DynamicAttrs, k::String) = getindex(A.attrs, k)
Base.haskey(A::DynamicAttrs, k::String) = haskey(A.attrs, k)
Base.length(A::DynamicAttrs) = length(A.attrs)
Base.iterate(A::DynamicAttrs, args...) = iterate(A.attrs, args...)
Base.pairs(A::DynamicAttrs) = pairs(A.attrs)

function Base.setindex!(d::DynamicAttrs, v::TAttrVal, k::String)
    if !isnothing(d.value_length_limit)
        v = _truncate(v, d.value_length_limit)
    end
    d.attrs[k] = v
end

n_dropped(a::DynamicAttrs) = n_dropped(a.attrs)

#####

Base.@kwdef struct Resource{A<:StaticAttrs}
    attributes::A = StaticAttrs()
    schema_url::String = ""
end

#####

Base.@kwdef struct InstrumentationInfo
    name::String = "Main"
    version::VersionNumber = v"0.0.1-dev"
end

ext: jl


text:
#--- 1D ---#
prim = zeros(3, 2)
prim[:, 1] .= [1.0, 3.0, 2.0]
prim[:, 2] .= [1 / 1864, -0.3, 1 / 1864]
# pure
Mu1, Mxi1, MuL1, MuR1 = KitBase.gauss_moments(prim[:, 1], 2)
Mu2, Mxi2, MuL2, MuR2 = KitBase.gauss_moments(prim[:, 2], 2)
# mixture
Mu, Mxi, MuL, MuR = KitBase.mixture_gauss_moments(prim, 2)

@test Mu[:, 1] ≈ Mu1 atol = 0.01
@test Mxi[:, 1] ≈ Mxi1 atol = 0.01
@test MuL[:, 1] ≈ MuL1 atol = 0.01
@test MuR[:, 1] ≈ MuR1 atol = 0.01
@test Mu[:, 2] ≈ Mu2 atol = 0.01
@test Mxi[:, 2] ≈ Mxi2 atol = 0.01
@test MuL[:, 2] ≈ MuL2 atol = 0.01
@test MuR[:, 2] ≈ MuR2 atol = 0.01

#--- 2D ---#
prim = zeros(4, 2)
prim[:, 1] .= [1.0, 3.0, 0.5, 1.3]
prim[:, 2] .= [1 / 1864, -1.3, 0.88, 2 / 1864]
# pure
Mu1, Mv1, Mxi1, MuL1, MuR1 = KitBase.gauss_moments(prim[:, 1], 1)
Mu2, Mv2, Mxi2, MuL2, MuR2 = KitBase.gauss_moments(prim[:, 2], 1)
# mixture
Mu, Mv, Mxi, MuL, MuR = KitBase.mixture_gauss_moments(prim, 1)

Muv1 = KitBase.moments_conserve(MuR1, Mv1, Mxi1, 3, 1, 1)
Muv2 = KitBase.moments_conserve(MuR2, Mv2, Mxi2, 3, 1, 1)
Muv = KitBase.mixture_moments_conserve(MuR, Mv, Mxi, 3, 1, 1)

@test Muv[:, 1] ≈ Muv1 atol = 0.01
@test Muv[:, 2] ≈ Muv2 atol = 0.01

#--- 3D ---#
prim = zeros(5, 2)
prim[:, 1] .= [1.0, 1.0, 0.5, 3.0, 2.0]
prim[:, 2] .= [1 / 1864, 0.1, -0.3, 2.0, 1 / 1864]
# pure
Mu1, Mv1, Mw1, MuL1, MuR1 = KitBase.gauss_moments(prim[:, 1], 0)
Mu2, Mv2, Mw2, MuL2, MuR2 = KitBase.gauss_moments(prim[:, 2], 0)
# mixture
Mu, Mv, Mw, MuL, MuR = KitBase.mixture_gauss_moments(prim, 0)

@test Mu[:, 1] ≈ Mu1 atol = 0.01
@test Mv[:, 1] ≈ Mv1 atol = 0.01
@test Mw[:, 1] ≈ Mw1 atol = 0.01
@test MuL[:, 1] ≈ MuL1 atol = 0.01
@test MuR[:, 1] ≈ MuR1 atol = 0.01

@test Mu[:, 2] ≈ Mu2 atol = 0.01
@test Mv[:, 2] ≈ Mv2 atol = 0.01
@test Mw[:, 2] ≈ Mw2 atol = 0.01
@test MuL[:, 2] ≈ MuL2 atol = 0.01
@test MuR[:, 2] ≈ MuR2 atol = 0.01

#--- generalized ---#
Mu1, Mv1, Mw1, MuL1, MuR1 = KitBase.gauss_moments(prim[:, 1], 0)
KitBase.moments_conserve_slope(zeros(5), Mu1, Mv1, Mw1, 2, 0, 0)

f1 = rand(16)
u1 = collect(-1:1/7.5:1)
w1 = ones(16)
KitBase.discrete_moments(f1, u1, w1, 1)
KitBase.moments_conserve(f1, u1, w1)
KitBase.moments_conserve(f1, f1, u1, w1)
KitBase.moments_conserve(f1, f1, f1, f1, u1, w1)

f2 = rand(16, 16)
u2 = rand(16, 16)
w2 = rand(16, 16)
KitBase.discrete_moments(f2, u2, w2, 1)
KitBase.moments_conserve(f2, u2, u2, w2)
KitBase.moments_conserve(f2, f2, u2, u2, w2)
KitBase.moments_conserve(f2, f2, f2, u2, u2, w2)

f3 = rand(16, 16, 16)
u3 = rand(16, 16, 16)
w3 = rand(16, 16, 16)
KitBase.discrete_moments(f3, u3, w3, 1)
KitBase.moments_conserve(f3, u3, u3, u3, w3)

KitBase.diatomic_moments_conserve(f1, f1, f1, u1, w1)
KitBase.diatomic_moments_conserve(f2, f2, f2, u2, u2, w2)

KitBase.pressure(f1, rand(3), u1, w1)
KitBase.pressure(f1, f1, rand(3), u1, w1, 2)
KitBase.pressure(f2, f2, rand(4), u2, u2, w2, 1)
KitBase.stress(f2, rand(4), u2, u2, w2)

KitBase.heat_flux(f1, rand(3), u1, w1)
KitBase.heat_flux(f1, f1, rand(3), u1, w1)
KitBase.heat_flux(f1, f1, f1, rand(3), u1, w1)
KitBase.heat_flux(f2, rand(4), u2, u2, w2)
KitBase.heat_flux(f2, f2, rand(4), u2, u2, w2)
KitBase.heat_flux(f2, f2, f2, rand(4), u2, u2, w2)
KitBase.heat_flux(f3, rand(5), u3, u3, u3, w3)


text:
using WAV
using ArgParse

push!(LOAD_PATH, pwd())

using Filter

function main(wav_file)

	s, fs = wavread(wav_file)

	sequence = decode(s, fs)

	print("The dialing code is: ")
	for s in sequence
		print(s)
	end
	println()
end

if PROGRAM_FILE == "SignalFiltering.jl"

	# Parse the argument
	s = ArgParseSettings("Audio filtering decoder")

	@add_arg_table s begin
		"--audio"
		help = "The audio to be decoded"
		required = true
	end

	parsed_args = parse_args(ARGS, s)

	main(parsed_args["audio"])
end

max_stars_repo_name: PatricioIribarneCatella/TP-Telefono

max_issues_repo_name: PatricioIribarneCatella/TP-Telefono

max_forks_repo_name: PatricioIribarneCatella/TP-Telefono


text:
# # Stratified Kelvin-Helmholtz instability
#
# ## Install dependencies
#
# First let's make sure we have all required packages installed.

# ```julia
# using Pkg
# pkg"add Oceananigans, JLD2, Plots"
# ```

# ## The physical domain
#
# We simulate Kelvin-Helmholtz instability in two-dimensions in ``x, z``
# and therefore assign `Flat` to the `y` direction,

using Oceananigans

grid = RectilinearGrid(size=(64, 64), x=(-5, 5), z=(-5, 5),
                       topology=(Periodic, Flat, Bounded))

# # The basic state
#
# We're simulating the instability of a sheared and stably-stratified basic state
# ``U(z)`` and ``B(z)``. Two parameters define our basic state: the Richardson number
# defined as,
#
# ```math
# Ri = \frac{∂_z B}{(∂_z U)^2}
# ```
#
# and the width of the stratification layer, ``h``.

shear_flow(x, y, z, t) = tanh(z)

stratification(x, y, z, t, p) = p.h * p.Ri * tanh(z / p.h)

U = BackgroundField(shear_flow)

B = BackgroundField(stratification, parameters=(Ri=0.1, h=1/4))

# Our basic state thus has a thin layer of stratification in the center of
# the channel, embedded within a thicker shear layer surrounded by unstratified fluid.

using Plots

zF = znodes(Face, grid)
zC = znodes(Center, grid)

Ri, h = B.parameters

kwargs = (ylabel="z", linewidth=3, label=nothing)

 U_plot = plot(shear_flow.(0, 0, zC, 0), zC; xlabel="U(z)", kwargs...)
 B_plot = plot([stratification(0, 0, z, 0, (Ri=Ri, h=h)) for z in zC], zC; xlabel="B(z)", color=:red, kwargs...)
Ri_plot = plot(@. Ri * sech(zF / h)^2 / sech(zF)^2, zF; xlabel="Ri(z)", color=:black, kwargs...) # Ri(z)= ∂_z B / (∂_z U)²; derivatives computed by hand

plot(U_plot, B_plot, Ri_plot, layout=(1, 3), size=(800, 400))

# In unstable flows it is often useful to determine the dominant spatial structure of the
# instability and the growth rate at which the instability grows.
# If the simulation idealizes a physical flow, this can be used to make
# predictions as to what should develop and how quickly.
# Since these instabilities are often attributed to a linear instability,
# we can determine information about the structure and the growth rate of the instability by analyzing the linear operator
# that governs small perturbations about a base state, or by solving for the linear dynamics.

# Here, we discuss first briefly linear instabilities and how one can obtain growth rates and structures of most unstable
# modes via eigenanalysis. Then we present an alternative method for approximating the eigenanalysis results when
# one does not have access to the linear dynamics or the linear operator about the base state.

# # Linear Instabilities
#
# The base state ``U(z)``, ``B(z)`` is a solution of the inviscid equations of motion. Whether the base state is
# stable or not is determined by whether small perturbations about this base state grow or decay. To formalize this,
# we study the linearized dynamics satisfied by perturbations about the base state:
# ```math
# \partial_t \Phi = L \Phi \, .
# ```
# where ``\Phi = (u, v, w, b)`` is a vector of the perturbation velocities ``u, v, w`` and perturbation buoyancy ``b``
# and ``L`` a linear operator that depends on the base state, ``L = L(U(z), B(z))`` (the `background_fields`).
# Eigenanalysis of the linear operator ``L`` determines the stability of the base state, such as the Kelvin-Helmholtz
# instability. That is, by using the ansantz
# ```math
# \Phi(x, y, z, t) = \phi(x, y, z) \, \exp(\lambda t) \, ,
# ```
# then ``\lambda`` and ``\phi`` are respectively eigenvalues and eigenmodes of ``L``, i.e., they obey
# ```math
# L \, \phi_j = \lambda_j \, \phi_j \quad j=1,2,\dots \, .
# ```
# From hereafter we'll use the convention that the eigenvalues are ordered according to their real part,
# ``\mathrm{Re}(\lambda_1) \ge \mathrm{Re}(\lambda_2) \ge \dotsb``.
#
# Remarks:
#
# As we touched upon briefly above, Oceananigans.jl, does not include a linearized version of the equations.
# Furthermore, Oceananigans.jl does not give us access to the linear operator ``L`` so that we can perform eigenanalysis.
# Below we discuss an alternative way of approximating the eigenanalysis results.
# The method boils down to solving the nonlinear equations while continually renormalize
# the magnitude of the perturbations to ensure that nonlinear terms
# (terms that are quadratic or higher in perturbations) remain negligibly small,
# i.e.,much smaller than the background flow.

# # The power method algorithm
#
# Successive application of ``L`` to a random initial state will eventually render it parallel
# with eigenmode ``\phi_1``:
# ```math
# \lim_{n \to \infty} L^n \Phi \propto \phi_1 \, .
# ```
# Of course, if ``\phi_1`` is an unstable mode (i.e., ``\sigma_1 = \mathrm{Re}(\lambda_1) > 0``), then successive application
# of ``L`` will lead to exponential amplification. (Similarly, if ``\sigma_1 < 0``, successive application of ``L`` will
# lead to exponential decay of ``\Phi`` down to machine precision.) Therefore, after each
# application of the linear operator ``L``, we rescale the output ``L \Phi`` back to a pre-selected amplitude.
#
# So, we initialize a `simulation` with random initial conditions with amplitude much less than those of
# the base state (which are ``O(1)``). Instead of "applying" ``L`` on our initial state, we evolve the
# (approximately) linear dynamics for interval ``\Delta \tau``. We measure how much the energy has grown
# during that interval, rescale the perturbations back to original energy amplitude and repeat.
# After some iterations the state will converge to the most unstable eigenmode.
#
# In summary, each iteration of the power method includes:
# - compute the perturbation energy, ``E_0``,
# - evolve the system for a time-interval ``\Delta \tau``,
# - compute the perturbation energy, ``E_1``,
# - determine the exponential growth of the most unstable mode during the interval
#   ``\Delta \tau`` as  ``\log(E_1 / E_0) / (2 \Delta \tau)``,
# - repeat the above until growth rate converges.
#
# By fiddling a bit with ``\Delta t`` we can get convergence after only a few iterations.
#
# Let's apply all these to our example.

# # The model

model = NonhydrostaticModel(timestepper = :RungeKutta3,
                              advection = UpwindBiasedFifthOrder(),
                                   grid = grid,
                               coriolis = nothing,
                      background_fields = (u=U, b=B),
                                closure = ScalarDiffusivity(ν=2e-4, κ=2e-4),
                               buoyancy = BuoyancyTracer(),
                                tracers = :b)

# We have included a "pinch" of viscosity and diffusivity in anticipation of what will follow furtherdown:
# viscosity and diffusivity will ensure numerical stability when we evolve the unstable mode to the point it becomes nonlinear.

# For this example, we take ``\Delta \tau = 15``.

simulation = Simulation(model, Δt=0.1, stop_iteration=150)

# Now some helper functions that will be used during for the power method algorithm.
#
# First a function that evolves the state for ``\Delta \tau`` and measure the energy growth
# over that period.

"""
    grow_instability!(simulation, e)

Grow an instability by running `simulation`.

Estimates the growth rate ``σ`` of the instability
using the fractional change in volume-mean kinetic energy,
over the course of the `simulation`

``
energy(t₀ + Δτ) / energy(t₀) ≈ exp(2 σ Δτ)
``

where ``t₀`` is the starting time of the simulation and ``t₀ + Δτ``
the ending time of the simulation. We thus find that the growth rate
is measured by

``
σ = log(energy(t₀ + Δτ) / energy(t₀)) / (2 Δτ) .
``
"""
function grow_instability!(simulation, energy)
    ## Initialize
    simulation.model.clock.iteration = 0
    t₀ = simulation.model.clock.time = 0
    compute!(energy)
    energy₀ = energy[1, 1, 1]

    ## Grow
    run!(simulation)

    ## Analyze
    compute!(energy)
    energy₁ = energy[1, 1, 1]
    Δτ = simulation.model.clock.time - t₀

    ## ½(u² + v²) ~ exp(2 σ Δτ)
    σ = growth_rate = log(energy₁ / energy₀) / 2Δτ

    return growth_rate
end
nothing # hide

# Finally, we write a function that rescales the state. The rescaling is done via computing the
# kinetic energy and then rescaling all flow fields so that the kinetic energy assumes a targetted value.
#
# (Measuring the perturbation growth via the kinetic energy works fine _unless_ an unstable mode _only_ has
# buoyancy structure. In that case, the total perturbation energy is more adequate.)

"""
    rescale!(model, energy; target_kinetic_energy=1e-3)

Rescales all model fields so that `energy = target_kinetic_energy`.
"""
function rescale!(model, energy; target_kinetic_energy=1e-6)
    compute!(energy)
    rescale_factor = √(target_kinetic_energy / energy[1, 1, 1])

    for f in merge(model.velocities, model.tracers)
        f .*= rescale_factor
    end

    return nothing
end

using Printf

# Another helper function for the power method,

"""
    convergence(σ)

Check if the growth rate has converged. If the array `σ` has at least 2 elements then returns the
relative difference between ``σ[end]`` and ``σ[end-1]``.
"""
convergence(σ) = length(σ) > 1 ? abs((σ[end] - σ[end-1]) / σ[end]) : 9.1e18 # pretty big (not Inf tho)
nothing # hide

# and the main function that performs the power method iteration.

"""
    estimate_growth_rate(simulation, energy, ω; convergence_criterion=1e-3)

Estimates the growth rate iteratively until the relative change
in the estimated growth rate ``σ`` falls below `convergence_criterion`.

Returns ``σ``.
"""
function estimate_growth_rate(simulation, energy, ω, b; convergence_criterion=1e-3)
    σ = []
    power_method_data = []
    compute!(ω)
    push!(power_method_data, (ω=collect(interior(ω)[:, 1, :]), b=collect(interior(b)[:, 1, :]), σ=deepcopy(σ)))

    while convergence(σ) > convergence_criterion
        compute!(energy)

        @info @sprintf("About to start power method iteration %d; kinetic energy: %.2e", length(σ)+1, energy[1, 1, 1])
        push!(σ, grow_instability!(simulation, energy))
        compute!(energy)

        @info @sprintf("Power method iteration %d, kinetic energy: %.2e, σⁿ: %.2e, relative Δσ: %.2e",
                       length(σ), energy[1, 1, 1], σ[end], convergence(σ))

        compute!(ω)
        rescale!(simulation.model, energy)
        push!(power_method_data, (ω=collect(interior(ω)[:, 1, :]), b=collect(interior(b)[:, 1, :]), σ=deepcopy(σ)))
    end

    return σ, power_method_data
end
nothing # hide

# # Eigenplotting
#
# A good algorithm wouldn't be complete without a good visualization,

u, v, w = model.velocities
b = model.tracers.b

perturbation_vorticity = Field(∂z(u) - ∂x(w))

xF, yF, zF = nodes(perturbation_vorticity)
xC, yC, zC = nodes(b)

eigentitle(σ, t) = length(σ) > 0 ? @sprintf("Iteration #%i; growth rate %.2e", length(σ), σ[end]) : @sprintf("Initial perturbation fields")

function eigenplot(ω, b, σ, t; ω_lim=maximum(abs, ω)+1e-16, b_lim=maximum(abs, b)+1e-16)

    kwargs = (xlabel="x", ylabel="z", linewidth=0, label=nothing, color=:balance, aspectratio=1)
    ω_title(t) = t == nothing ? @sprintf("vorticity") : @sprintf("vorticity at t = %.2f", t)
    b_title(t) = t == nothing ? @sprintf("buoyancy")  : @sprintf("buoyancy at t = %.2f", t)

    plot_ω = contourf(xF, zF, clamp.(ω, -ω_lim, ω_lim)';
                      levels = range(-ω_lim, stop=ω_lim, length=20),
                       xlims = (xF[1], xF[grid.Nx]),
                       ylims = (zF[1], zF[grid.Nz]),
                       clims = (-ω_lim, ω_lim),
                       title = ω_title(t), kwargs...)

    plot_b = contourf(xC, zC, clamp.(b, -b_lim, b_lim)';
                    levels = range(-b_lim, stop=b_lim, length=20),
                     xlims = (xC[1], xC[grid.Nx]),
                     ylims = (zC[1], zC[grid.Nz]),
                     clims = (-b_lim, b_lim),
                     title = b_title(t), kwargs...)

    return plot(plot_ω, plot_b, layout=(1, 2), size=(800, 380))
end

function power_method_plot(ω, b, σ, t)
    plot_growthrates = scatter(σ, xlabel = "Power iteration",
                                  ylabel = "Growth rate",
                                   title = eigentitle(σ, nothing),
                                   label = nothing)

    plot_eigenmode = eigenplot(ω, b, σ, nothing)

    return plot(plot_growthrates, plot_eigenmode, layout=@layout([A{0.25h}; B]), size=(800, 600))
end
nothing # hide

# # Rev your engines...
#
# We initialize the power iteration with random noise and rescale to have a `target_kinetic_energy`

using Random, Statistics

mean_perturbation_kinetic_energy = Field(Average(1/2 * (u^2 + w^2)))
noise(x, y, z) = randn()
set!(model, u=noise, w=noise, b=noise)
rescale!(simulation.model, mean_perturbation_kinetic_energy, target_kinetic_energy=1e-6)
growth_rates, power_method_data = estimate_growth_rate(simulation, mean_perturbation_kinetic_energy, perturbation_vorticity, b)

@info "Power iterations converged! Estimated growth rate: $(growth_rates[end])"

# # Powerful convergence
#
# We animate the power method steps. A scatter plot illustrates how the growth rate converges
# as the power method iterates,

anim_powermethod = @animate for i in 1:length(power_method_data)
    power_method_plot(power_method_data[i].ω, power_method_data[i].b, power_method_data[i].σ, nothing)
end

mp4(anim_powermethod, "powermethod.mp4", fps = 1) # hide

# # Now for the fun part
#
# Now we simulate the nonlinear evolution of the eigenmode
# we've isolated for a few e-folding times ``1/\sigma``,

## Reset the clock
model.clock.iteration = 0
model.clock.time = 0

estimated_growth_rate = growth_rates[end]

simulation.stop_time = 5 / estimated_growth_rate
simulation.stop_iteration = 9.1e18 # pretty big (not Inf tho)

## Rescale the eigenmode
initial_eigenmode_energy = 5e-5
rescale!(simulation.model, mean_perturbation_kinetic_energy, target_kinetic_energy=initial_eigenmode_energy)

# Let's save and plot the perturbation vorticity and buoyancy and also the total vorticity and
# buoyancy (perturbation + basic state). It'll be also neat to plot the kinetic energy time-series
# and confirm it grows with the estimated growth rate.

total_vorticity = Field(∂z(u) + ∂z(model.background_fields.velocities.u) - ∂x(w))

total_b = Field(b + model.background_fields.tracers.b)

simulation.output_writers[:vorticity] =
    JLD2OutputWriter(model, (ω=perturbation_vorticity, Ω=total_vorticity, b=b, B=total_b, KE=mean_perturbation_kinetic_energy),
                     schedule = TimeInterval(0.10 / estimated_growth_rate),
                     filename = "kelvin_helmholtz_instability.jld2",
                     overwrite_existing = true)

# And now we...

@info "*** Running a simulation of Kelvin-Helmholtz instability..."
run!(simulation)

# ## Pretty things
#
# First we plot the nonlinear equilibration of the perturbation fields together
# with the evolution of the kinetic energy,

@info "Making a neat movie of stratified shear flow..."

function plot_energy_timeseries(t, measured_KE, σ, KEᵢ, stop_time)
    t_segment = [0, stop_time]
    predicted_KE = KEᵢ * exp.(2σ * [0, stop_time])

    energy_plot = plot(t_segment, predicted_KE,
             label = "~ exp(2 σ t)",
            legend = :topleft,
                lw = 2,
             color = :black,
             yaxis = :log,
             xlims = (0, simulation.stop_time),
             ylims = (initial_eigenmode_energy, 1e-1),
            xlabel = "time",
            ylabel = "kinetic energy")

    energy_plot = plot!(t, measured_KE,
              label = "perturbation kinetic energy",
              lw = 6,
              alpha = 0.5)

    return energy_plot
end

filepath = simulation.output_writers[:vorticity].filepath
ω_timeseries = FieldTimeSeries(filepath, "ω")
b_timeseries = FieldTimeSeries(filepath, "b")
KE_timeseries = FieldTimeSeries(filepath, "KE")
times = ω_timeseries.times
KE = []

anim_perturbations = @animate for (i, t) in enumerate(times)
    @info "Plotting frame $i from iteration $iteration..."

    ω_snapshot = interior(ω_timeseries, :, 1, :, i)
    b_snapshot = interior(b_timeseries, :, 1, :, i)
    push!(KE, KE_timeseries[i][1, 1, 1])

    energy_plot = plot_energy_timeseries(times[1:i], KE, estimated_growth_rate, initial_eigenmode_energy, simulation.stop_time)
    eigenmode_plot = eigenplot(ω_snapshot, b_snapshot, nothing, t; ω_lim=1, b_lim=0.05)

    layout = @layout [A{0.6h}; B]
    plot(eigenmode_plot, energy_plot; layout, size=(800, 600))
end

mp4(anim_perturbations, "kelvin_helmholtz_instability_perturbations.mp4", fps = 8) # hide

# And then the same for total vorticity & buoyancy of the fluid.

Ω_timeseries = FieldTimeSeries(filepath, "Ω")
B_timeseries = FieldTimeSeries(filepath, "B")
times = Ω_timeseries.times
KE = []

anim_total = @animate for (i, t) in enumerate(times)
    @info "Plotting frame $i from iteration $iteration..."

    Ω_snapshot = interior(Ω_timeseries, :, 1, :, i)
    B_snapshot = interior(B_timeseries, :, 1, :, i)
    push!(KE, KE_timeseries[i][1, 1, 1])

    energy_plot = plot_energy_timeseries(times[1:i], KE, estimated_growth_rate, initial_eigenmode_energy, simulation.stop_time)
    eigenmode_plot = eigenplot(Ω_snapshot, B_snapshot, nothing, t; ω_lim=1, b_lim=0.05)

    layout = @layout [A{0.6h}; B]
    plot(eigenmode_plot, energy_plot; layout, size=(800, 600))
end

mp4(anim_total, "kelvin_helmholtz_instability_total.mp4", fps = 8) # hide


text:
function measurement(z,c,μ,ξ,Ω,graph)
    n = size(graph,1)-1
    qnoise = diagm(0 => [5.0,0.02])
    nnew = maximum(c)-n

    if nnew > 0
        osize = length(ξ)
        nsize = osize+2*nnew
        Ω₁ = zeros(nsize,nsize)
        Ω₁[1:osize,1:osize] = Ω
        Ω = Ω₁
        append!(ξ,zeros(2*nnew))
        append!(μ,inverse_measurement(μ[1:3],z[:,c .> n])[:])
        padd = spzeros(Bool,nnew,n+1)
        padd1 = spzeros(Bool,n+nnew+1,nnew)
        graph = [[graph;padd] padd1]
    end

    jcb = jacobian_measurement(μ,c)
    ẑ = equation_measurement(μ,c)
    qnoiseinv = inv(qnoise)

    for (i,j) ∈ enumerate(c)
        h = jcb[:,:,i]
        dz = z[:,i]-ẑ[:,i]
        dz[2] = measure.(dz[2])
        ind = [1:3;2j+2;2j+3]

        ξ[ind] .+= h'*qnoiseinv*(dz+h*μ[ind])
        Ω[ind,ind] .+= h'*qnoiseinv*h
    end

    c = unique(c)
    rows = ones(length(c))
    cols = c .+ 1
    sg = size(graph,1)
    adjacency = sparse(rows,cols,true,sg,sg)
    graph .|= adjacency
    graph .|= adjacency'

    return Ω,graph

end


text:
"""
This package implements basic unit conversions between units commonly used in spectroscopy which are not of the same dimensions,
such as Hz and m. It exports the function `convert_unit` which follows the same convention as uconvert from Unitful.
"""
module SpectroscopicUnits

using Reexport
@reexport using Unitful
import PhysicalConstants.CODATA2018: c_0, ħ, h

export convert_unit

# Constants for the most commonly used units
export Hz, radHz, eV, m, cm, mm ,μm, nm

const Hz = u"Hz"
const radHz = u"rad*Hz"
const eV = u"eV"
const m = u"m"
const cm = u"cm"
const mm = u"mm"
const μm = u"μm"
const nm = u"nm"

function convert_unit(unit::Unitful.LengthUnits, quantity::Quantity{T,D,U}) where {T<:Real,D, U<:Unitful.FrequencyUnits}
    return uconvert(unit, c_0 / quantity)
end

function convert_unit(unit::Unitful.LengthUnits, quantity::Quantity{T,D,typeof(u"rad*Hz")}) where {T<:Real,D}
    return uconvert(unit, 2.0 * pi * c_0 / quantity)
end

function convert_unit(unit::Unitful.LengthUnits, quantity::Quantity{T,D,U}) where {U <: Unitful.EnergyUnits, T<:Real,D}
    energy = uconvert(u"J", quantity)
    freq = uconvert(Hz, energy / h)
    return convert_unit(unit, freq)
end

function convert_unit(unit::Unitful.LengthUnits, quantity::Quantity{T,D,typeof(cm^-1)}) where {T<:Real,D}
    return uconvert(unit, 1.0/quantity)
end

convert_unit(unit::U, quantity::Unitful.Length) where {U <: Unitful.FrequencyUnits} = uconvert(unit, c_0 / (quantity |> m))
convert_unit(unit::typeof(radHz), quantity::Unitful.Length) = uconvert(unit, 2.0 * pi * c_0 / (quantity |> m))

function convert_unit(unit::Unitful.EnergyUnits, quantity::Unitful.Length)
    freq = uconvert(Hz, c_0 / quantity )
    energy = uconvert(u"J", h * freq)
    return uconvert(unit, energy)
end

function convert_unit(unit::typeof(cm^-1), quantity::Unitful.Length)
    return uconvert(unit, 1.0/quantity)
end

function convert_unit(unit::U, quantity::Quantity{T,D,V}) where {T<:Real,D,U,V}
    temp = convert_unit(m, quantity)
    return convert_unit(unit, temp)
end

convert_unit(unit::U, quantity::Quantity{T,D,U}) where {T<:Real,D,U} = quantity

end


text:
function Base.map!(
    op::UnaryUnion, C::GBArray, A::GBArray;
    mask = nothing, accum = nothing, desc = nothing
)
    mask, accum = _handlenothings(mask, accum)
    desc = _handledescriptor(desc; in1=A)
    op = getoperator(op, eltype(A))
    accum = getaccum(accum, eltype(C))
    libgb.GrB_Matrix_apply(C, mask, accum, op, parent(A), desc)
    return C
end

Base.map!(op, C::GBArray, A::GBArray; mask = nothing, accum = nothing, desc = nothing) =
    map!(UnaryOp(op), C, A; mask, accum, desc)

function Base.map!(
    op::UnaryUnion, A::GBArray;
    mask = nothing, accum = nothing, desc = nothing
)
    return map!(op, A, A; mask, accum, desc)
end

Base.map!(op, A::GBArray; mask = nothing, accum = nothing, desc = nothing) =
    map!(UnaryOp(op), A; mask, accum, desc)

function Base.map(
    op::UnaryUnion, A::GBArray;
    mask = nothing, accum = nothing, desc = nothing
)
    t = inferoutputtype(A, op)
    return map!(op, similar(A, t), A; mask, accum, desc)
end

Base.map(op, A::GBArray; mask = nothing, accum = nothing, desc = nothing) =
    map(UnaryOp(op), A; mask, accum, desc)

function Base.map!(
    op::BinaryUnion, C::GBArray, x, A::GBArray;
    mask = nothing, accum = nothing, desc = nothing
)
    mask, accum = _handlenothings(mask, accum)
    desc = _handledescriptor(desc; in2=A)
    op = getoperator(op, optype(eltype(A), typeof(x)))
    accum = getaccum(accum, eltype(C))
    libgb.scalarmatapply1st[optype(typeof(x), eltype(A))](C, mask, accum, op, x, parent(A), desc)
    return C
end

function Base.map!(
    op, C::GBArray, x, A::GBArray;
    mask = nothing, accum = nothing, desc = nothing
)
    map!(BinaryOps.BinaryOp(op), C, x, A; mask, accum, desc)
end

function Base.map!(
    op::BinaryUnion, x, A::GBArray;
    mask = nothing, accum = nothing, desc = nothing
)
    return map!(op, A, x, A; mask, accum, desc)
end

function Base.map!(
    op, x, A::GBArray;
    mask = nothing, accum = nothing, desc = nothing
)
    map!(BinaryOps.BinaryOp(op), x, A; mask, accum, desc)
end

function Base.map(
    op::BinaryUnion, x, A::GBArray;
    mask = nothing, accum = nothing, desc = nothing
)
    t = inferoutputtype(A, op)
    return map!(op, similar(A, t), x, A; mask, accum, desc)
end

"""
    map(op::Union{Function, AbstractUnaryOp}, A::GBArray; kwargs...)::GBArray
    map(op::Union{Function, AbstractBinaryOp}, A::GBArray, x; kwargs...)::GBArray
    map(op::Union{Function, AbstractBinaryOp}, x, A::GBArray, kwargs...)::GBArray

Transform a GBArray by applying `op` to each element.

UnaryOps and single argument functions apply elementwise in the usual fashion.
BinaryOps and two argument functions require the additional argument `x` which is
    substituted as the first or second operand of `op` depending on its position.

# Arguments
- `op::Union{Function, AbstractUnaryOp, AbstractBinaryOp}`
- `A::GBArray`
- `x`: Position dependent argument to binary operators.

# Keywords
- `mask::Union{Nothing, GBVecOrMat} = nothing`: optional mask.
- `accum::Union{Nothing, AbstractBinaryOp} = nothing`: binary accumulator operation
    where `C[i,j] = accum(C[i,j], T[i,j])` where T is the result of this function before accum is applied.
- `desc::Union{Nothing, Descriptor} = nothing`
"""
function Base.map(
    op, x, A::GBArray;
    mask = nothing, accum = nothing, desc = nothing
)
    map(BinaryOps.BinaryOp(op), x, A; mask, accum, desc)
end

function Base.map!(
    op::BinaryUnion, C::GBArray, A::GBArray, x;
    mask = nothing, accum = nothing, desc = nothing
)
    mask, accum = _handlenothings(mask, accum)
    desc = _handledescriptor(desc; in1=A)
    op = getoperator(op, optype(eltype(A), typeof(x)))
    accum = getaccum(accum, eltype(C))
    libgb.scalarmatapply2nd[optype(typeof(x), eltype(A))](C, mask, accum, op, parent(A), x, desc)
    return C
end

function Base.map!(
    op, C::GBArray, A::GBArray, x;
    mask = nothing, accum = nothing, desc = nothing
)
    map!(BinaryOps.BinaryOp(op), C, A, x; mask, accum, desc)
end

function Base.map!(
    op::BinaryUnion, A::GBArray, x;
    mask = nothing, accum = nothing, desc = nothing
)
    return map!(op, A, A, x; mask, accum, desc)
end

function Base.map!(
    op, A::GBArray, x;
    mask = nothing, accum = nothing, desc = nothing
)
    map!(BinaryOps.BinaryOp(op), A, A, x; mask, accum, desc)
end

function Base.map(
    op::BinaryUnion, A::GBArray, x;
    mask = nothing, accum = nothing, desc = nothing
)
    t = inferoutputtype(A, op)
    return map!(op, similar(A, t), A, x; mask, accum, desc)
end

function Base.map(
    op, A::GBArray, x;
    mask = nothing, accum = nothing, desc = nothing
)
    map(BinaryOps.BinaryOp(op), A, x; mask, accum, desc)
end

Base.:+(x::valid_union, u::GBArray; mask = nothing, accum = nothing, desc = nothing) =
    map(BinaryOps.PLUS, x, u; mask, accum, desc)
Base.:+(u::GBArray, x::valid_union; mask = nothing, accum = nothing, desc = nothing) =
    map(BinaryOps.PLUS, u, x; mask, accum, desc)

Base.:*(x::valid_union, u::GBArray; mask = nothing, accum = nothing, desc = nothing) =
    map(BinaryOps.TIMES, x, u; mask, accum, desc)
Base.:*(u::GBArray, x::valid_union; mask = nothing, accum = nothing, desc = nothing) =
    map(BinaryOps.TIMES, u, x; mask, accum, desc)

Base.:-(x::valid_union, u::GBArray; mask = nothing, accum = nothing, desc = nothing) =
    map(BinaryOps.MINUS, x, u; mask, accum, desc)
Base.:-(u::GBArray, x::valid_union; mask = nothing, accum = nothing, desc = nothing) =
    map(BinaryOps.MINUS, u, x; mask, accum, desc)

Base.:-(u::GBArray) = map(UnaryOps.AINV, u)

"""
    mask!(C::GBArray, A::GBArray, mask::GBArray)

Apply a mask to matrix `A`, storing the results in C.

"""
function mask!(C::GBArray, A::GBArray, mask::GBArray; structural = false, complement = false)
    desc = Descriptor()
    structural && (desc.structural_mask=true)
    complement && (desc.complement_mask=true)
    map!(identity, C, A; mask, desc)
    return C
end

"""
    mask(A::GBArray, mask::GBArray)

Apply a mask to matrix `A`.
"""
function mask(A::GBArray, mask::GBArray; structural = false, complement = false)
    return mask!(similar(A), A, mask; structural, complement)
end


text:
using MultivariateTests
using Compat
using Compat.Test
using Compat.DelimitedFiles
using StatsBase
using HypothesisTests

# Columns are calcium, iron, protein, vitamin A, vitamin C
nutrient = readdlm(joinpath(@__DIR__, "data", "nutrient.txt"))[:,2:end]

# Columns are information, similarities, arithmetic, picture completion
wechsler = readdlm(joinpath(@__DIR__, "data", "wechsler.txt"))[:,2:end]

# Columns are type, length, left, right, bottom, top, diag
swiss = readdlm(joinpath(@__DIR__, "data", "swiss3.txt"))
genuine = convert(Matrix{Float64}, swiss[view(swiss, :, 1) .== "real", 2:end])
counterfeit = convert(Matrix{Float64}, swiss[view(swiss, :, 1) .== "fake", 2:end])

# Columns are survey answers: husbands' to questions 1-4, wives' to questions 1-4
spouse = readdlm(joinpath(@__DIR__, "data", "spouse.txt"))

@testset "Utility functions" begin
    MT = MultivariateTests

    @test MT.checkdims(rand(3, 2), rand(4, 2)) == (2, 3, 4)
    # Mismatched number of variables
    @test_throws DimensionMismatch MT.checkdims(rand(4, 3), rand(4, 6))
    # Empty observations
    @test_throws ArgumentError MT.checkdims(rand(0, 4), rand(2, 4))

    Sx = [ 3.0 -1.5 0.0
          -1.5  1.0 0.5
           0.0  0.5 1.0]
    Sy = fill(4.5, (3, 3))
    out = [3.5 0.5      1.5
           0.5 2.166667 1.833333
           1.5 1.833333 2.166667]
    @test MT._poolcov!(Sx, 2, Sy, 1) ≈ out atol=1e-6
    # Input is modified
    @test Sx ≈ out atol=1e-6

    # Positive semi-definite but not positive definite
    P = [ 1.0000  0.7426  0.1601 -0.7000 0.5500
          0.7426  1.0000 -0.2133 -0.5818 0.5000
          0.1601 -0.2133  1.0000 -0.1121 0.1000
         -0.7000 -0.5818 -0.1121  1.0000 0.4500
          0.5500  0.5000  0.1000  0.4500 1.0000]
    # Positive definite
    Q = [ 1.0 -0.5  0.0
         -0.5  1.0 -0.5
          0.0 -0.5  1.0]
    @test @inferred(MT.At_Binv_A(ones(5), P)) ≈ -0.8008792 atol=1e-6
    @test @inferred(MT.At_Binv_A(ones(3), Q)) ≈ 10.0
end

@testset "Generalizations of sample variance" begin
    @test @inferred(genvar(nutrient)) ≈ 2.8310418e19 atol=1e-6
    @test @inferred(totalvar(nutrient)) ≈ 2.83266877e6 atol=1e-6

    X = [1 2 5
         4 1 6
         4 0 4]
    @test @inferred(genvar(X)) ≈ 0.0
    @test @inferred(totalvar(X)) ≈ 5.0

    x = rand(Float32, 10)
    @test genvar(x) == totalvar(x) == var(x)

    it = (xᵢ for xᵢ in x)
    @test genvar(it) == totalvar(it) == var(it)
end

@testset "Partial correlation" begin
    @test @inferred(partialcor(wechsler[:,1], wechsler[:,2], wechsler[:,3:4])) ≈ 0.7118787 atol=1e-6

    w = PartialCorTest(wechsler[:,1], wechsler[:,2], wechsler[:,3:4])
    let out = sprint(show, w)
        @test contains(out, "reject h_0") && !contains(out, "fail to")
    end
    let ci = confint(w)
        @test first(ci) ≈ 0.4963917 atol=1e-6
        @test last(ci) ≈ 0.8447292 atol=1e-6
    end
    @test nobs(w) == 37
    @test dof(w) == 33
    @test pvalue(w) < 0.00001

    X = [ 2 1 0
          4 2 0
         15 3 1
         20 4 1]
    @test @inferred(partialcor(X[:,1], X[:,2], X[:,3])) ≈ 0.919145 atol=1e-6

    x = PartialCorTest(view(X,:,1), view(X,:,2), view(X,:,3))
    @test contains(sprint(show, x), "fail to reject")
    @test confint(x) == (-1.0, 1.0)
    @test nobs(x) == 4
    @test dof(x) == 1
    @test pvalue(x) ≈ 0.25776212 atol=1e-6
end

@testset "One sample Hotelling's T²" begin
    t = OneSampleHotellingT2(nutrient, [1000, 15, 60, 800, 75])
    @test nobs(t) == 737
    @test dof(t) == (5, 732)
    @test pvalue(t) ≈ 0.0 atol=eps()
    @test t.T² ≈ 1758.5413137 atol=1e-6
    @test t.F ≈ 349.7968048 atol=1e-6
    let out = sprint(show, t)
        @test contains(out, "reject h_0") && !contains(out, "fail to")
    end

    # Paired
    p = OneSampleHotellingT2(spouse[:,1:4], spouse[:,5:end])
    @test nobs(p) == 30
    @test dof(p) == (4, 26)
    @test pvalue(p) ≈ 0.039369144 atol=1e-6
    @test p.T² ≈ 13.127840261 atol=1e-6
    @test p.F ≈ 2.942446955 atol=1e-6
    let out = sprint(show, p)
        @test contains(out, "reject h_0") && !contains(out, "fail to")
    end
end

@testset "Two sample Hotelling's T²" begin
    eq = EqualCovHotellingT2(genuine, counterfeit)
    @test nobs(eq) == (100, 100)
    @test dof(eq) == (6, 193)
    @test pvalue(eq) ≈ 0.0 atol=eps()
    @test eq.T² ≈ 2412.4506855 atol=1e-6
    @test eq.F ≈ 391.9217023 atol=1e-6
    @test contains(sprint(show, eq), "reject h_0")

    un = UnequalCovHotellingT2(genuine, counterfeit)
    @test nobs(un) == (100, 100)
    @test dof(un) == (6, 193)
    @test pvalue(un) ≈ 0.0 atol=eps()
    @test un.T² ≈ 2412.4506855 atol=1e-6
    @test un.F ≈ 391.9217023 atol=1e-6
    @test contains(sprint(show, un), "reject h_0")
end

@testset "Bartlett's test" begin
    b = BartlettsTest(genuine, counterfeit)
    @test nobs(b) == (100, 100)
    @test dof(b) == 21
    @test pvalue(b) ≈ 0.0 atol=1e-10
    @test b.L′ ≈ 121.8991235 atol=1e-6
    @test contains(sprint(show, b), "reject h_0")
end


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

using Markdown

# ╔═╡ 5b2ee40e-a2b8-11ea-0fef-c35fe6918860
md"""
# The tower of Hanoi

The tower of hanoi is a famous puzzle.

![setup of the tower of a hanoi](https://upload.wikimedia.org/wikipedia/commons/0/07/Tower_of_Hanoi.jpeg)

The game consists of three rods with disks stacked on top of them. The puzzle will start with all the disk in a stack on one of the rods (like in the picture). The goal is to move all the discs to a single stack on  different rod.

To move the disks, you have to follow the following rules:

* You can move only one disk at a time.
* For each move, you have to take the upper disk from one of the stacks, and place it on top of another stack or empty rod.
* You cannot place a larger disk on top of a smaller disk.

This notebook will define a Julia implementation of the puzzle. It's up to you to write an algorithm that solves it.
"""

# ╔═╡ 95fbd0d2-a2b9-11ea-0682-fdf783251797
md"""
## Setting up the game pieces

What does a Julia implementation look like? We're not really interested in writing code that will manipulate physical disks. Our final goal is a function that will give us a _recipe_ to solve the tower of hanoi, which is just a list of moves to make. Because of that, we can use a lot of abstraction in our implementation, and keep the data structures as simple as possible.

To start, we have to define some representation of the disks and the stacks. The disks have one important property, which is that they are ordered. We can use integers to represent them.
"""

# ╔═╡ 620d6834-a2ba-11ea-150a-2132bb54e4b3
num_disks = 8

# ╔═╡ 35ada214-a32c-11ea-0da3-d5d494b28467
md"""(Side note: the number of disks is arbitrary. When testing your function, you may want to set it to 1 or 2 to start.)"""

# ╔═╡ 7243cc8e-a2ba-11ea-3f29-99356f0cdcf4
all_disks = 1:num_disks

# ╔═╡ 7e1ba2ac-a2ba-11ea-0134-2f61ed75be18
md"""
A single stack can be represented as an array with all the disks in it. We will list them from top to bottom.
"""

# ╔═╡ 43781a52-a339-11ea-3803-e56e2d08aa83
first_stack = collect(all_disks)

# ╔═╡ b648ab70-a2ba-11ea-2dcc-55b630e44325
md"""
Now we have to make three of those.
"""

# ╔═╡ 32f26f80-a2bb-11ea-0f2a-3fc631ada63d
starting_stacks = [first_stack, [], []]

# ╔═╡ e347f1de-a2bb-11ea-06e7-87cca6f2a240
md"""
## Defining the rules

Now that we have our "game board", we can implement the rules.

To start, we make two functions for states. A state of the game is just an array of stacks.

We will define a function that checks if a state is okay according to the rules. To be legal, all the stacks should be in the correct order, so no larger disks on top of smaller disks.

Another good thing to check: no disks should have appeared or disappeared since we started!
"""

# ╔═╡ 512fa6d2-a2bd-11ea-3dbe-b935b536967b
function islegal(stacks)
	#check order
	order_correct_per_stack = map(stacks) do stack
		stack == sort(stack)
	end
	order_correct = sum(order_correct_per_stack) == length(stacks)

	#check if we use the same disk set that we started with
	disks_in_state = sort([disk for stack in stacks for disk in stack])
	disks_complete = disks_in_state == all_disks

	order_correct && disks_complete
end

# ╔═╡ c56a5858-a2bd-11ea-1d96-77eaf5e74925
md"""
Another function for states: check if we are done! We can assume that we already checked if the state was legal. So we know that all the disks are there and they are ordered correctly.  To check if we are finished, we just need to check if the last stack contains all the disks.
"""

# ╔═╡ d5cc41e8-a2bd-11ea-39c7-b7df8de6ae3e
function iscomplete(stacks)
	last(stacks) == all_disks
end

# ╔═╡ 53374f0e-a2c0-11ea-0c91-97474780721e
md"""
Now the only rules left to implement are the rules for moving disks.

We could implement this as another check on states, but it's easier to just write a legal `move` function. Your solution will specify moves for the `move` function, so this will be the only way that the stacks are actually manipulated. That we, we are sure that nothing fishy is happening.

We will make our `move` function so that its input consists of a state of the game, and instructions for what to do. Its output will be the new state of the game.

So what should those instructions look like? It may seem intuitive to give a _disk_ that should be moved, but that's more than we need. After all, we are only allowed to take the top disk from one stack, and move it to the top of another. So we only have to say which _stacks_ we are moving between.

(Note that the `move` function is okay with moving a larger disk on top of a smaller disk. We already implemented that restriction in `islegal`.)
"""

# ╔═╡ e915394e-a2c0-11ea-0cd9-1df6fd3c7adf
function move(stacks, source::Int, target::Int)
	#check if the from stack if not empty
	if length(stacks[source]) == 0
		throw(ErrorException("Error: attempted to move disk from empty stack"))
	end

	#make a copy of the stacks
	new_stacks = [stack for stack in stacks]

	disk = first(new_stacks[source]) #take disk
	new_stacks[source] = new_stacks[source][2:end] #remove from old stack
	new_stacks[target] = vcat([disk], new_stacks[target]) #put on new stack

	new_stacks
end

# ╔═╡ 87b2d164-a2c4-11ea-3095-916628943879
md"""
## Solving the problem

We have implemented the game pieces and the rules, so you can start working on your solution.

To do this, you can fill in the `solve()` function. This function should give a solution for the `starting_stacks`, by moving all the disks from stack 1 to stack 3.

As output, `solve` should give a recipe, that tells us what to do. This recipe should be an array of moves. Each moves is a `(source, target)` tuple, specifying from which stack to which stack you should move.

For example, it might look like this:
"""

# ╔═╡ 29b410cc-a329-11ea-202a-795b31ce5ad5
wrong_solution = [(1,2), (2,3), (2,1), (1,3)]

# ╔═╡ ea24e778-a32e-11ea-3f11-dbe9d36b1011
md"""
Now you can work on building an actual solution. Some tips:
* `solve()` can keep track of the board if you want, but it doesn't have to.
* The section below will actually run your moves, which is very useful for checking them.
* If you want, you can change `num_disks` to 1 or 2. That can be a good starting point.
"""

# ╔═╡ 010dbdbc-a2c5-11ea-34c3-837eae17416f
function solve() ::Array{Tuple{Int, Int}}
	[]
end

# ╔═╡ 3eb3c4c0-a2c5-11ea-0bcc-c9b52094f660
md"""
## Checking solutions

This is where we can check a solution. We start with a function that takes our recipe and runs it.
"""

# ╔═╡ 4709db36-a327-11ea-13a3-bbfb18da84ce
function run_solution(moves, start = starting_stacks)
	all_states = map(x -> start, 1:length(moves) + 1)

	for i = 1:length(moves)
		try
			source, target = moves[i]
			all_states[i + 1] = move(all_states[i],source, target)
		catch

			all_states[i + 1] = []
		end
	end

	all_states
end

# ╔═╡ 372824b4-a330-11ea-2f26-7b9a1ad018f1
md"""
You can use this function to see what your solution does.

If `run_solution` encounters an error, it will give empty arrays from that point onwards. Look at what happens in the `wrong_solution` version and compare it to the moves in `wrong_solution`.
"""

# ╔═╡ 9173b174-a327-11ea-3a69-9f7525f2e7b4
run_solution(solve())

# ╔═╡ d2227b40-a329-11ea-105c-b585d5fcf970
run_solution(wrong_solution)

# ╔═╡ bb5088ec-a330-11ea-2c41-6b8b92724b3b
md"""
Now that we have way to run a recipe, we can check if its output is correct. We will check if all the intermediate states are legal and the final state is the finished puzzle.
"""

# ╔═╡ 10fb1c56-a2c5-11ea-2a06-0d8c36bfa138
function check_solution(moves, start = starting_stacks)
	try
		#run the solution
		all_states = run_solution(moves, start)

		#check if each state is legal
		legal_moves = map(islegal, all_states)
		all_legal = sum(legal_moves) == length(all_states)

		#check if the final state is is the completed puzzle
		complete = (iscomplete ∘ last)(all_states)

		all_legal && complete
	catch
		#return false if we encountered an error
		return false
	end
end

# ╔═╡ 8ea7f944-a329-11ea-22cc-4dbd11ec0610
check_solution(solve())

# ╔═╡ e54add0a-a330-11ea-2eeb-1d42f552ba38
if check_solution(solve())
	if num_disks >= 8
		md"""
		Congratulations, your solution works!
		"""
	else
		md"""
		Your solution works for $(num_disks) disks. Change `num_disks` to see if it works for 8 or more.
		"""
	end
else
	md"""
	The `solve()` function doesn't work yet. Keep working on it!
	"""
end

# ╔═╡ Cell order:
# ╟─5b2ee40e-a2b8-11ea-0fef-c35fe6918860
# ╟─95fbd0d2-a2b9-11ea-0682-fdf783251797
# ╠═620d6834-a2ba-11ea-150a-2132bb54e4b3
# ╟─35ada214-a32c-11ea-0da3-d5d494b28467
# ╠═7243cc8e-a2ba-11ea-3f29-99356f0cdcf4
# ╟─7e1ba2ac-a2ba-11ea-0134-2f61ed75be18
# ╠═43781a52-a339-11ea-3803-e56e2d08aa83
# ╟─b648ab70-a2ba-11ea-2dcc-55b630e44325
# ╠═32f26f80-a2bb-11ea-0f2a-3fc631ada63d
# ╟─e347f1de-a2bb-11ea-06e7-87cca6f2a240
# ╠═512fa6d2-a2bd-11ea-3dbe-b935b536967b
# ╟─c56a5858-a2bd-11ea-1d96-77eaf5e74925
# ╠═d5cc41e8-a2bd-11ea-39c7-b7df8de6ae3e
# ╟─53374f0e-a2c0-11ea-0c91-97474780721e
# ╠═e915394e-a2c0-11ea-0cd9-1df6fd3c7adf
# ╟─87b2d164-a2c4-11ea-3095-916628943879
# ╠═29b410cc-a329-11ea-202a-795b31ce5ad5
# ╟─ea24e778-a32e-11ea-3f11-dbe9d36b1011
# ╠═010dbdbc-a2c5-11ea-34c3-837eae17416f
# ╟─3eb3c4c0-a2c5-11ea-0bcc-c9b52094f660
# ╠═4709db36-a327-11ea-13a3-bbfb18da84ce
# ╟─372824b4-a330-11ea-2f26-7b9a1ad018f1
# ╠═9173b174-a327-11ea-3a69-9f7525f2e7b4
# ╠═d2227b40-a329-11ea-105c-b585d5fcf970
# ╟─bb5088ec-a330-11ea-2c41-6b8b92724b3b
# ╠═10fb1c56-a2c5-11ea-2a06-0d8c36bfa138
# ╠═8ea7f944-a329-11ea-22cc-4dbd11ec0610
# ╟─e54add0a-a330-11ea-2eeb-1d42f552ba38


text:
using JuMP, FactCheck, AmplNLWriter

## Solve test problem 1 (Synthesis of processing system) in
 #  M. Duran & I.E. Grossmann, "An outer approximation algorithm for
 #  a class of mixed integer nonlinear programs", Mathematical
 #  Programming 36, pp. 307-339, 1986.  The problem also appears as
 #  problem synthes1 in the MacMINLP test set.
 #
 #  min   5 x4 + 6 x5 + 8 x6 + 10 x1 - 7 x3 -18 math.log(x2 + 1)
 #       - 19.2 math.log(x1 - x2 + 1) + 10
 #  s.t.  0.8 math.log(x2 + 1) + 0.96 math.log(x1 - x2 + 1) - 0.8 x3 >= 0
 #        math.log(x2 + 1) + 1.2 math.log(x1 - x2 + 1) - x3 - 2 x6 >= -2
 #        x2 - x1 <= 0
 #        x2 - 2 x4 <= 0
 #        x1 - x2 - 2 x5 <= 0
 #        x4 + x5 <= 1
 #        0 <= x1 <= 2
 #        0 <= x2 <= 2
 #        0 <= x3 <= 1
 #        x1, x2, x3 continuous
 #        x4, x5, x6 binary
 #
 #
 #  The solution is (1.30098, 0, 1, 0, 1, 0).
 ##

if !isdefined(:solver); solver = BonminNLSolver(); end

m = Model(solver=solver)
x_U = [2,2,1]
@variable(m, x_U[i] >= x[i=1:3] >= 0)
@variable(m, y[4:6], Bin)

@NLobjective(m, Min, 10 + 10*x[1] - 7*x[3] + 5*y[4] + 6*y[5] + 8*y[6] - 18*log(x[2]+1) - 19.2*log(x[1]-x[2]+1))
@NLconstraint(m, 0.8*log(x[2] + 1) + 0.96*log(x[1] - x[2] + 1) - 0.8*x[3] >= 0)
@NLconstraint(m, log(x[2] + 1) + 1.2*log(x[1] - x[2] + 1) - x[3] - 2*y[6] >= -2)
@NLconstraint(m, x[2] - x[1] <= 0)
@NLconstraint(m, x[2] - 2*y[4] <= 0)
@NLconstraint(m, x[1] - x[2] - 2*y[5] <= 0)
@NLconstraint(m, y[4] + 2*y[5]*0.5 <= 1)

context("example: jump_minlp") do
    @fact solve(m) --> :Optimal
    @fact getvalue(x)[:] --> roughly([1.30098, 0.0, 1.0], 1e-5)
    @fact getvalue(y)[:] --> roughly([0.0, 1.0, 0.0], 1e-5)
    @fact getobjectivevalue(m) --> roughly(6.00975, 1e-5)
end


text:
isequaldlm(m1, m2, t) = isequal(m1, m2) && (eltype(m1) == eltype(m2) == t)

let dlm_data = readdlm(joinpath("perf", "kernel", "imdb-1.tsv"), '\t')
    @test size(dlm_data) == (31383,3)
    @test dlm_data[12345,2] == "Gladiator"
    @test dlm_data[31383,3] == 2005
    @test dlm_data[1,1] == "McClure, Marc (I)"
end

@test isequaldlm(readdlm(IOBuffer("1\t2\n3\t4\n5\t6\n")), [1. 2; 3 4; 5 6], Float64)
@test isequaldlm(readdlm(IOBuffer("1\t2\n3\t4\n5\t6\n"), Int), [1 2; 3 4; 5 6], Int)

@test size(readcsv(IOBuffer("1,2,3,4"))) == (1,4)
@test size(readcsv(IOBuffer("1,2,3,"))) == (1,4)
@test size(readcsv(IOBuffer("1,2,3,4\n"))) == (1,4)
@test size(readcsv(IOBuffer("1,2,3,\n"))) == (1,4)
@test size(readcsv(IOBuffer("1,2,3,4\n1,2,3,4"))) == (2,4)
@test size(readcsv(IOBuffer("1,2,3,4\n1,2,3,"))) == (2,4)
@test size(readcsv(IOBuffer("1,2,3,4\n1,2,3"))) == (2,4)

@test size(readcsv(IOBuffer("1,2,3,4\r\n"))) == (1,4)
@test size(readcsv(IOBuffer("1,2,3,4\r\n1,2,3\r\n"))) == (2,4)
@test size(readcsv(IOBuffer("1,2,3,4\r\n1,2,3,4\r\n"))) == (2,4)
@test size(readcsv(IOBuffer("1,2,3,\"4\"\r\n1,2,3,4\r\n"))) == (2,4)

@test size(readdlm(IOBuffer("1 2 3 4\n1 2 3"))) == (2,4)
@test size(readdlm(IOBuffer("1\t2 3 4\n1 2 3"))) == (2,4)
@test size(readdlm(IOBuffer("1\t 2 3 4\n1 2 3"))) == (2,4)
@test size(readdlm(IOBuffer("1\t 2 3 4\n1 2 3\n"))) == (2,4)
@test size(readdlm(IOBuffer("1,,2,3,4\n1,2,3\n"), ',')) == (2,5)

let result1 = reshape({"", "", "", "", "", "", 1.0, 1.0, "", "", "", "", "", 1.0, 2.0, "", 3.0, "", "", "", "", "", 4.0, "", "", ""}, 2, 13),
    result2 = reshape({1.0, 1.0, 2.0, 1.0, 3.0, "", 4.0, ""}, 2, 4)

    @test isequaldlm(readdlm(IOBuffer(",,,1,,,,2,3,,,4,\n,,,1,,,1\n"), ','), result1, Any)
    @test isequaldlm(readdlm(IOBuffer("   1    2 3   4 \n   1   1\n")), result2, Any)
    @test isequaldlm(readdlm(IOBuffer("   1    2 3   4 \n   1   1\n"), ' '), result1, Any)
    @test isequaldlm(readdlm(IOBuffer("1 2\n3 4 \n")), [[1.0, 3.0] [2.0, 4.0]], Float64)
end

let result1 = reshape({"", "", "", "", "", "", "भारत", 1.0, "", "", "", "", "", 1.0, 2.0, "", 3.0, "", "", "", "", "", 4.0, "", "", ""}, 2, 13)
    @test isequaldlm(readdlm(IOBuffer(",,,भारत,,,,2,3,,,4,\n,,,1,,,1\n"), ',') , result1, Any)
end

let result1 = reshape({1.0, 1.0, 2.0, 2.0, 3.0, 3.0, 4.0, ""}, 2, 4)
    @test isequaldlm(readdlm(IOBuffer("1\t 2 3 4\n1 2 3")), result1, Any)
    @test isequaldlm(readdlm(IOBuffer("1\t 2 3 4\n1 2 3 ")), result1, Any)
    @test isequaldlm(readdlm(IOBuffer("1\t 2 3 4\n1 2 3\n")), result1, Any)
    @test isequaldlm(readdlm(IOBuffer("1,2,3,4\n1,2,3\n"), ','), result1, Any)
    @test isequaldlm(readdlm(IOBuffer("1,2,3,4\n1,2,3"), ','), result1, Any)
    @test isequaldlm(readdlm(IOBuffer("1,2,3,4\r\n1,2,3\r\n"), ','), result1, Any)
    @test isequaldlm(readdlm(IOBuffer("1,2,3,\"4\"\r\n1,2,3\r\n"), ','), result1, Any)
end

let result1 = reshape({"abc", "hello", "def,ghi", " \"quote\" ", "new\nline", "world"}, 2, 3),
    result2 = reshape({"abc", "line\"", "\"hello\"", "\"def", "", "\" \"\"quote\"\" \"", "ghi\"", "", "world", "\"new", "", ""}, 3, 4)

    @test isequaldlm(readdlm(IOBuffer("abc,\"def,ghi\",\"new\nline\"\n\"hello\",\" \"\"quote\"\" \",world"), ','), result1, Any)
    @test isequaldlm(readdlm(IOBuffer("abc,\"def,ghi\",\"new\nline\"\n\"hello\",\" \"\"quote\"\" \",world"), ',', quotes=false), result2, Any)
end

let result1 = reshape({"t", "c", "", "c"}, 2, 2),
    result2 = reshape({"t", "\"c", "t", "c"}, 2, 2)
    @test isequaldlm(readdlm(IOBuffer("t  \n\"c\" c")), result1, Any)
    @test isequaldlm(readdlm(IOBuffer("t t \n\"\"\"c\" c")), result2, Any)
end

@test isequaldlm(readcsv(IOBuffer("\n1,2,3\n4,5,6\n\n\n")), reshape({"",1.0,4.0,"","","",2.0,5.0,"","","",3.0,6.0,"",""}, 5, 3), Any)

let x = [1,2,3], y = [4,5,6], io = IOBuffer()
    writedlm(io, zip(x,y), ",  ")
    seek(io, 0)
    @test readcsv(io) == [x y]
end

let x = ["abc", "def\"ghi", "jk\nl"], y = [1, ",", "\"quoted\""], io = IOBuffer()
    writedlm(io, zip(x,y), ',')
    seek(io, 0)
    @test readcsv(io) == [x y]
end

let x = ["a" "b"; "d" ""], io = IOBuffer()
    writedlm(io, x)
    seek(io, 0)
    @test readdlm(io) == x
end

let x = ["\"hello\"", "world\""], io = IOBuffer()
    print(io, x)
    @assert takebuf_string(io) == "\"hello\"\nworld\"\n"

    writedlm(io, x)
    @assert takebuf_string(io) == "\"\"\"hello\"\"\"\n\"world\"\"\"\n"
end

# source: http://www.i18nguy.com/unicode/unicode-example-utf8.zip
let i18n_data = ["Origin (English)", "Name (English)", "Origin (Native)", "Name (Native)",
        "Australia", "Nicole Kidman", "Australia", "Nicole Kidman",
        "Austria", "Johann Strauss", "Österreich", "Johann Strauß",
        "Belgium (Flemish)", "Rene Magritte", "België", "René Magritte",
        "Belgium (French)", "Rene Magritte", "Belgique", "René Magritte",
        "Belgium (German)", "Rene Magritte", "Belgien", "René Magritte",
        "Bhutan", "Gonpo Dorji", "འབྲུག་ཡུལ།", "མགོན་པོ་རྡོ་རྗེ།",
        "Canada", "Celine Dion", "Canada", "Céline Dion",
        "Canada - Nunavut (Inuktitut)", "Susan Aglukark", "ᓄᓇᕗᒻᒥᐅᑦ", "ᓱᓴᓐ ᐊᒡᓗᒃᑲᖅ",
        "Democratic People's Rep. of Korea", "LEE Sol-Hee", "조선 민주주의 인민 공화국", "이설희",
        "Denmark", "Soren Hauch-Fausboll", "Danmark", "Søren Hauch-Fausbøll",
        "Denmark", "Soren Kierkegaard", "Danmark", "Søren Kierkegård",
        "Egypt", "Abdel Halim Hafez", "ﻣﺼﺮ", "ﻋﺑﺪﺍﻠﺣﻟﻳﻢ ﺤﺎﻓﻅ",
        "Egypt", "Om Kolthoum", "ﻣﺼﺮ", "ﺃﻡ ﻛﻟﺛﻭﻡ",
        "Eritrea", "Berhane Zeray", "ብርሃነ ዘርኣይ", "ኤርትራ",
        "Ethiopia", "Haile Gebreselassie", "ኃይሌ ገብረሥላሴ", "ኢትዮጵያ",
        "France", "Gerard Depardieu", "France", "Gérard Depardieu",
        "France", "Jean Reno", "France", "Jean Réno",
        "France", "Camille Saint-Saens", "France", "Camille Saint-Saëns",
        "France", "Mylene Demongeot", "France", "Mylène Demongeot",
        "France", "Francois Truffaut", "France", "François Truffaut",
        "France (Braille)", "Louis Braille", "⠋⠗⠁⠝⠉⠑", "⠇⠕⠥⠊⠎⠀<BR>⠃⠗⠁⠊⠇⠇⠑",
        "Georgia", "Eduard Shevardnadze", "საქართველო", "ედუარდ შევარდნაძე",
        "Germany", "Rudi Voeller", "Deutschland", "Rudi Völler",
        "Germany", "Walter Schultheiss", "Deutschland", "Walter Schultheiß",
        "Greece", "Giorgos Dalaras", "Ελλάς", "Γιώργος Νταλάρας",
        "Iceland", "Bjork Gudmundsdottir", "Ísland", "Björk Guðmundsdóttir",
        "India (Hindi)", "Madhuri Dixit", "भारत", "माधुरी दिछित",
        "Ireland", "Sinead O'Connor", "Éire", "Sinéad O'Connor",
        "Israel", "Yehoram Gaon", "ישראל", "יהורם גאון",
        "Italy", "Fabrizio DeAndre", "Italia", "Fabrizio De André",
        "Japan", "KUBOTA Toshinobu", "日本", "久保田    利伸",
        "Japan", "HAYASHIBARA Megumi", "日本", "林原 めぐみ",
        "Japan", "Mori Ogai", "日本", "森鷗外",
        "Japan", "Tex Texin", "日本", "テクス テクサン",
        "Norway", "Tor Age Bringsvaerd", "Noreg", "Tor Åge Bringsværd",
        "Pakistan (Urdu)", "Nusrat Fatah Ali Khan", "پاکستان", "نصرت فتح علی خان",
        "People's Rep. of China", "ZHANG Ziyi", "中国", "章子怡",
        "People's Rep. of China", "WONG Faye", "中国", "王菲",
        "Poland", "Lech Walesa", "Polska", "Lech Wałęsa",
        "Puerto Rico", "Olga Tanon", "Puerto Rico", "Olga Tañón",
        "Rep. of China", "Hsu Chi", "臺灣", "舒淇",
        "Rep. of China", "Ang Lee", "臺灣", "李安",
        "Rep. of Korea", "AHN Sung-Gi", "한민국", "안성기",
        "Rep. of Korea", "SHIM Eun-Ha", "한민국", "심은하",
        "Russia", "Mikhail Gorbachev", "Россия", "Михаил Горбачёв",
        "Russia", "Boris Grebenshchikov", "Россия", "Борис Гребенщиков",
        "Slovenia", "\"Frane \"\"Jezek\"\" Milcinski", "Slovenija", "Frane Milčinski - Ježek",
        "Syracuse (Sicily)", "Archimedes", "Συρακούσα", "Ἀρχιμήδης",
        "Thailand", "Thongchai McIntai", "ประเทศไทย", "ธงไชย แม็คอินไตย์",
        "U.S.A.", "Brad Pitt", "U.S.A.", "Brad Pitt",
        "Yugoslavia (Cyrillic)", "Djordje Balasevic", "Југославија", "Ђорђе Балашевић",
        "Yugoslavia (Latin)", "Djordje Balasevic", "Jugoslavija", "Đorđe Balašević"]

    i18n_arr = transpose(reshape(i18n_data, 4, int(floor(length(i18n_data)/4))))
    i18n_buff = PipeBuffer()
    writedlm(i18n_buff, i18n_arr, ',')
    @test i18n_arr == readcsv(i18n_buff)
end

max_stars_repo_name: alyst/julia

max_issues_repo_name: alyst/julia

max_forks_repo_name: alyst/julia


text:
using Mosek

function intProduction(productNames     :: Vector{String},
                       processNames     :: Vector{String},
                       timeRequirements :: Array{Float64,2},
                       timeResources    :: Vector{Float64},
                       profit           :: Vector{Float64})
    n = length(productNames)
    m = length(processNames)
    maketask() do t
        appendvars(t,n)
        appendcons(t,m)

        for j in 1:n
            putvarname(t,Int32(j),productNames[j])
            putvarbound(t,Int32(j),MSK_BK_LO,0.0,Inf)
            putvartype(t,Int32(j),MSK_VAR_TYPE_INT)
        end

        for i in 1:m
            putconname(t,Int32(i),processNames[i])
            putconbound(t,i, MSK_BK_UP, -Inf,timeResources[i])
            putarow(t,i,Int32[1:n...],timeRequirements[:,i])
        end

        putobjsense(t,MSK_OBJECTIVE_SENSE_MAXIMIZE)
        putclist(t,Int32[1:n...],profit)

        function callbackfunc(code::Callbackcode, dinf, iinf, liinf)
            if code == MSK_CALLBACK_NEW_INT_MIO && 0 != solutiondef(t,MSK_SOL_ITG)
                xx = getxx(t,MSK_SOL_ITG)
                println("New solution: $xx")
            end
            0
        end

        putcallbackfunc(t,callbackfunc)

        optimize(t,"mosek://solve.mosek.com:30080")

        xx = getxx(t,MSK_SOL_ITG)
        println("Final solution: $xx")
    end
end

timeRequirements = Float64[ 2.0  3.0  2.0  4.0
                            4.0  2.0  3.0  0.0
                            1.0  2.0  1.0  1.0
                            3.0  3.0  2.0  1.0 ]
timeResources    = Float64[100.0, 50.0, 30.0, 60.0]
profit           = Float64[1.7, 2.3, 1.0, 1.6]
processNames     = String["Assembly", "Polish","Varnish", "Packing" ]
productNames     = String["Chairs", "Tables", "Beds", "CuttingBoards"]

intProduction(productNames,processNames,timeRequirements,timeResources,profit)


text:
# NOTE: run with $ julia --depwarn=yes for these tests to pass
seed!(RNG, 1337)

@testset "histogram" begin
    x = randn(RNG, 10000)
    @test_logs (:warn, r"`symb`.+deprecated") histogram(x, symb = "#", closed = :right)
    @test_logs (:warn, r"`bins`.+deprecated") histogram(x, bins = 5, closed = :right)
    @test_logs (:warn, r"nbins::Int.+deprecated") histogram(x, 5, closed = :right)
end

@testset "barplot" begin
    @test_logs (:warn, r"`symb`.+deprecated") barplot(
        ["Paris", "Moskau"],
        [2.244, 11.92];
        symb = "#",
    )
end

@testset "plot" begin
    p = @inferred Plot([1, 2], [-1, 1])
    @test_logs (:warn, r"`annotate!`.+renamed") annotate!(p, :l, "text 1")
    @test_logs (:warn, r"`annotate!`.+renamed") annotate!(p, :l, "text 2", :red)
    @test_logs (:warn, r"`annotate!`.+renamed") annotate!(p, :l, "text 3", color = :blue)

    for i in 1:2
        @test_logs (:warn, r"`annotate!`.+renamed") annotate!(
            p,
            :l,
            i,
            "$i",
            color = :green,
        )
        @test_logs (:warn, r"`annotate!`.+renamed") annotate!(p, :l, i, "$i", :yellow)
    end

    @test_logs (:warn, r"`UnicodePlots.savefig`.+text") savefig(p, tempname() * ".png")
end

@testset "spy" begin
    @test_logs (:warn, r"`color = :automatic`.+deprecated") spy(
        _stable_sprand(RNG, 10, 10, 0.15),
        color = :automatic,
    )
end


text:
function Base.show(io::IO, tree::RTree)
    print(io, typeof(tree))
    print(io, "(variant="); print(io, tree.variant)
    print(io, ", tight_mbrs="); print(io, tree.tight_mbrs)
    print(io, ", nearmin_overlap="); print(io, tree.nearmin_overlap)
    print(io, ", fill_factor="); print(io, tree.fill_factor)
    print(io, ", split_factor="); print(io, tree.split_factor)
    print(io, ", reinsert_factor="); print(io, tree.reinsert_factor)

    print(io, ", leaf_capacity="); print(io, capacity(Leaf, tree))
    print(io, ", branch_capacity="); print(io, capacity(Branch, tree)); println(io, ")")
    print(io, length(tree)); print(io, " element(s) in "); print(io, height(tree));
    print(io, " level(s) ("); print(io, join(reverse(tree.nnodes_perlevel), ", ")); println(io, " node(s) per level):")
    _show(io, tree.root, recurse=true, indent=1)
end

function _show(io::IO, node::Node; recurse::Bool=false, indent::Int=0)
    for _ in 1:indent; print(io, ' '); end
    print(io, "level="); print(io, level(node));
    print(io, " nchildren="); print(io, length(node));
    print(io, " mbr=("); print(io, mbr(node).low); print(io, ", "); print(io, mbr(node).high); print(io, ")")
    if recurse && !isempty(node)
        print(io, ':')
        if node isa Leaf
            for child in children(node)
                println(io)
                for _ in 1:(indent+1); print(io, ' '); end
                show(io, child)
            end
        else
            for child in children(node)
                println(io)
                _show(io, child, recurse=true, indent=indent+1)
            end
        end
    end
end


text:
#A Small utility library for Quaternions and their functions
import Base.size, Base.getindex, Base.setindex!, Base.IndexStyle
import Base.*, Base.+, Base.-, Base./, Base.norm

mutable struct Quaternion <: AbstractArray{Float64,1}
    q1::Float64 #<- scalar part
    q2::Float64 #
    q3::Float64 # <- Vector part
    q4::Float64 #
    norm::Bool #Boolean to indicate whether quaternion is normalized
end

Base.size(::Quaternion) = (4,)
Base.IndexStyle(::Type{<:Quaternion}) = IndexLinear()
function Base.getindex(q::Quaternion, i::Int)
    if (i>0 && i<5)
        if (i == 1)
            return q.q1
        elseif i== 2
            return q.q2
        elseif i == 3
            return q.q3
        elseif i == 4
            return q.q4
        end
    else
        error("index out of range")
    end
end

function setindex!(q::Quaternion, v::T, i::Int) where {T <: Real}
    if (i>0 && i<5)
        if (i == 1)
            q.q1 = Float64(v)
        elseif i== 2
            q.q2 = Float64(v)
        elseif i == 3
            q.q3 = Float64(v)
        elseif i == 4
            q.q4 = Float64(v)
        end
        q.norm = false
    else
        error("index out of range")
    end
    return q
end

Quaternion(q1::S, qvec::Vector{T}, norm=false) where {T,S <: Real} = Quaternion(Float64(q1), Float64(qvec[1]), Float64(qvec[2]), Float64(qvec[3]), norm)
Quaternion(qvec::Vector{T}) where {T <: Real} = Quaternion(0., qvec)
Quaternion(q1::T) where {T <: Real} = Quaternion(Float64(q1), Float64(zeros(3)))
vec(q::Quaternion) = [[q[2] q[3] q[4]]...]
scalar(q::Quaternion) = q[1]
conjugate(q::Quaternion) = Quaternion(scalar(q), [-vec(q)...], isnormalized(q))
function inverse(q::Quaternion)
    if isnormalized(q)
        return conjugate(q)
    else
        return conjugate(q)/(norm(q))^2
    end
end

norm(q::Quaternion) = sqrt(q[1]^2 + q[2]^2 + q[3]^2 + q[4]^2)
function isnormalized(q::Quaternion)
    q.norm ? true : false
end

function normalize(q::Quaternion)
    if !q.norm
        q = q/norm(q)
        q.norm = true
    end
    return q
end

(-)(q::Quaternion) = Quaternion(-q[1], -q[2], -q[3], -q[4], false)
(+)(q1::Quaternion, q2::Quaternion) = Quaternion(q1[1] + q2[1], q1[2] + q2[2], q1[3] + q2[3], q1[4] + q2[4],false)
(-)(q1::Quaternion, q2::Quaternion) = Quaternion(q1[1] - q2[1], q1[2] - q2[2], q1[3] - q2[3], q1[4] - q2[4],false)
(*)(q1::Quaternion, q2::Quaternion) = Quaternion(q1[1]*q2[1] - dot(vec(q1),vec(q2)),
    q2[1].*vec(q1) + q1[1].*vec(q2) + cross(vec(q1),vec(q2)),false)
(/)(q::Quaternion, s::T) where {T <: Real} = Quaternion(q[1]/s, q[2:4]./s, false)

angleaxis2quat(angle::T, axis::Vector{B}) where {T <: Real, B <: Real} = Quaternion(cos(angle/2), axis.*sin(angle/2),true)

"""
```julia
DCM = quat2DCM(q::Quaternion)
```
Computes the DCM corresponding the input Quaternion.
"""
function quat2DCM(q::Quaternion)

  if norm(q) < eps()
    q = Quaternion(1., zeros(3))
  end

  q12 = q[1]^2
  q22 = q[2]^2
  q32 = q[3]^2
  q42 = q[4]^2
  q1q2 = q[1]*q[2]
  q3q4 = q[3]*q[4]
  q1q3 = q[1]*q[3]
  q2q4 = q[2]*q[4]
  q2q3 = q[2]*q[3]
  q1q4 = q[1]*q[4]

  return transpose([[(q12 + q22 - q32 - q42) 2*(q2q3 - q1q4) 2*(q1q3 + q2q4)];
          [2*(q2q3 + q1q4) (q12 - q22 + q32 - q42) 2*(q3q4 - q1q2)];
          [2*(-q1q3 + q2q4) 2*(q1q2 + q3q4) (q12 - q22 - q32 + q42)]])
end

#Helper function for rotation a vector with the Quaternion
function RotateVector(q::Quaternion, v::Vector{T}) where {T <: Real}

  if norm(q) < eps()
    q = Quaternion(1., zeros(3))
  end

  q12 = q[1]^2
  q22 = q[2]^2
  q32 = q[3]^2
  q42 = q[4]^2
  q1q2 = q[1]*q[2]
  q3q4 = q[3]*q[4]
  q1q3 = q[1]*q[3]
  q2q4 = q[2]*q[4]
  q2q3 = q[2]*q[3]
  q1q4 = q[1]*q[4]

  return [(q12 + q22 - q32 - q42)*v[1] + 2*(q2q3 + q1q4)*v[2] + 2*(-q1q3 + q2q4)v[3],
                2*(q2q3 - q1q4)*v[1] + (q12 - q22 + q32 - q42)*v[2] +  2*(q1q2 + q3q4)*v[3],
                2*(q1q3 + q2q4)*v[1] + 2*(q3q4 - q1q2)*v[2] + (q12 - q22 - q32 + q42)*v[3]]
end


text:
module LiveHTTPSample

using HTTP
using JSON3
using Revise

include("LiveHTTP.jl")

# Configuration
ip_addr = "0.0.0.0"
http_port = 8080
websocket_port = 8081

# Main entry point to start the program
function main()
    println("URL: http://host:$http_port  websocket_port: $websocket_port")
    start_http_and_websocket_servers(ip_addr, http_port, http_router, websocket_port, websocket_router)
end

# Map HTTP URLs to functions
function http_router(request)
    route_http_request(request;
        public=Dict(
            "/"      => url_root,
            "/somepage" => url_somepage,
        ),
        public_path=Dict(
            "/static/" => url_static,
        )
    )
end

# Map WebSocket messages to functions
function websocket_router(websocket, request)
    route_websocket_request(websocket, request;
        public=Dict(
            "button_clicked" => ws_button_clicked,
            "error" => ws_error,
        ),
        private=Dict(
            "private" => ws_private,
        ),
        default=ws_default
    )
end

function ws_button_clicked(ws, request)
    println("Received button_clicked message from browser: ", request[:msg])
    write(ws, JSON3.write(Dict(:type=>"update_span", :msg=>request[:msg])))
end

function ws_private(ws, request)
    println("Received private message from browser: ", request[:msg])
end

function ws_error(ws, request)
    println("Received error message from browser: ", request[:msg])
end

function ws_default(ws, request)
    println("Server received unexpected message from browser: ", JSON3.write(request))
    write(ws, JSON3.write(Dict(:type=>"log", :msg=>"Unexpected message was sent to server")))
end

# Home Page
function url_root(request, sessionid)
    response(sessionid, 200, [], body="""

                <title>Home Page</title>

Home Page

                Static Page

                Some Page

        """
    )
end

# Some Page
function url_somepage(request, sessionid)
    css = """
        body { background:lightgray; }
        """

    js = """
        var socket;

        window.addEventListener('load', function(event) {

            var elem = document.getElementById("clicked")
            elem.focus();

            console.log("Connecting websocket to " + location.hostname + "...");

            // Create WebSocket connection
            socket = new WebSocket('ws://' + location.hostname + ':$websocket_port');

            // Connection opened
            socket.addEventListener('open', function (event) {
                console.log("Websocket connected");
            });

            // Listen for messages
            socket.addEventListener('message', function (event) {
                var request = JSON.parse(event.data);
                switch(request.type) {
                    case 'update_span':
                        ws_update_span(request);
                        break;
                    case 'log':
                        console.log("Msg: ", request.msg);
                        break;
                    case 'alert':
                        alert("Msg: " + request.msg);
                        break;
                    case 'error':
                        console.log("Error: ", request.msg);
                        break;
                    default:
                        console.log("Server sent an unexpected message to browser: ", request);
                }
            });
        });

        window.addEventListener('unload', function(event) {
            socket.close();
        });

        function button_onclick() {
            console.log("Sending button_clicked message to server");
            socket.send(JSON.stringify({'type': 'button_clicked', 'msg': event.target.innerHTML}));
        }

        function ws_update_span(request) {
            console.log("Received update_span message from server");
            var clicked = document.getElementById('clicked');
            clicked.value = 'You clicked ' + request.msg;
        }
        """

    html = """

                <title>Some Page</title>

Some Page

            Home

            Static Page

            <input id="clicked"></input>

            <button id="button1" onclick="button_onclick();">Button 1</button>
            <button id="button2" onclick="button_onclick();">Button 2</button>

        """

    response(sessionid, 200, [], body=html)
end

function url_static(request, path, sessionid)
    pkg_path = joinpath((splitpath(@__DIR__)[1:end-1])...)
    path = url_to_path(path, joinpath(pkg_path, "static"))
    println("Path: ", path)
    if !ismissing(path)
        if isdir(path)
            response(sessionid, 500, [], body="Sorry, no folders yet")
        elseif endswith(path, ".css")
            response(sessionid, 200, ["Content-Type"=>"text/css"], body=read(path))
        else
            response(sessionid, 200, [], body=read(path))
        end
    else
        response(sessionid, 404, [], body="Sorry, bad request")
    end
end

# Return a filesystem path under some directory or missing if it does not exist
function url_to_path(req_path::AbstractString, dir::AbstractString)
    r_parts = HTTP.URIs.unescapeuri.(split(HTTP.URI(req_path).path[2:end], "/"))
    fs_path = normpath(joinpath(r_parts...))
    println(fs_path)
    fs_path = joinpath(dir, fs_path)
    println(fs_path, " ", isfile(fs_path), " " , isdir(fs_path))
    isfile(fs_path) || isdir(fs_path) ? fs_path : missing
end

end # module


text:
include("BMA.jl")
include("GLM.jl")
include("RFM.jl")

max_stars_repo_name: jungangl/forecast-machine-learning

max_issues_repo_name: jungangl/forecast-machine-learning

max_forks_repo_name: jungangl/forecast-machine-learning


text:
module RandomV06

    import Random

    include("random_v06.jl")

    const V06 = Val{6}
    const V07 = Val{7}
    const Vcur = Val{7}

    include("rand_ver.jl")

    export rand_ver,
        seed_ver!,
        rand_ver!,
        randn_ver,
        randn_ver!,
        randexp_ver,
        randexp_ver!,
        bitrand_ver,
        randstring_ver,
        randsubseq_ver,
        randsubseq_ver!,
        shuffle_ver,
        shuffle_ver!,
        randperm_ver,
        randcycle_ver

end


text:
module FastDMTransform

const KDM = 4.148808e3 # MHz^2 s pc^-1 cm^3

include("cpu.jl")

end


text: @testset "posterior" begin end


text:
using MathOptInterface
const MOI = MathOptInterface
const MOIT = MOI.Test

const OPTIMIZER = SCIP.Optimizer(display_verblevel=0)
const CONFIG_DIRECT = MOIT.Config(
    atol=5e-3, rtol=1e-4,
    exclude=Any[
        MOI.ConstraintDual, MOI.ConstraintName, MOI.DualObjectiveValue, MOI.VariableBasisStatus, MOI.ConstraintBasisStatus,
    ],
)

@testset "MOI unit tests" begin
    excluded = copy(MOI_BASE_EXCLUDED)
    append!(excluded, [
        "test_linear_integration", # Can not delete variable while model contains constraints
        "test_basic_VectorOfVariables_SecondOrderCone",
        "test_conic_SecondOrderCone_nonnegative_post_bound",
        "test_variable_delete_SecondOrderCone",
        "test_modification_func_scalaraffine_",
        "test_modification_func_vectoraffine_",
    ])
    MOIT.runtests(
        OPTIMIZER,
        CONFIG_DIRECT,
        warn_unsupported=false,
        exclude = excluded,
    )

end


text:
using Zygote:@show,hook
function dem(x, msg)
    # @show (norm(x), typeof(x), msg)
    return x
end

# todo version that doesn't have an fft built in
function (shears::ConvFFT{D,OT,A,B,C,PD,P})(x) where {D,OT,A,B,C,PD,P <: Tuple}
    if typeof(shears.weight) <: CuArray && !(typeof(x) <: Flux.CuArray ||
                                            typeof(x) <: CuArray)
        error("don't try to apply a CuArray to a non-CuArray")
    end
    xbc, usedInds = applyBC(x, shears.bc, ndims(shears.weight) - 1)

    Forward = shears.fftPlan[1]
    if size(xbc) != size(Forward)
        xbc = reshape(xbc, size(Forward))
    end
    x̂ = Forward * xbc

    nextLayer = hook(x -> dem(x, "should be after internalConvFFT"), internalConvFFT(x̂,
                                shears.weight, usedInds, shears.fftPlan[2],
                                shears.bias, shears.analytic))
    tmp = hook(x -> dem(x, "nonlin"), shears.σ.(nextLayer))
    return hook(x -> dem(x, "entering ConvFFT"), tmp)
end

function (shears::ConvFFT)(x)
    if typeof(shears.weight) <: CuArray && !(typeof(x) <: CuArray)
        error("don't try to apply a CuArray to a non-CuArray")
    end
    xbc, usedInds = applyBC(x, shears.bc, ndims(shears.weight) - 1)

    F = shears.fftPlan
    if F isa Nothing
        F = makePlan(eltype(x), outType(shears), shears.w, size(x), shears.bc)
        Forward = F[1]
    end

    if F isa Tuple
        Forward = F[1]
    else
        Forward = F
    end
    if size(xbc) != size(Forward)
        xbc = reshape(xbc, size(Forward))
    end
    x̂ = Forward * xbc

    nextLayer = internalConvFFT(x̂, shears.weight, usedInds, F,
                                shears.bias, shears.analytic)
    return shears.σ.(nextLayer)
end

abstract type TransformTypes end # If the input/wavelets are either real or
# analytic, there are efficiency gains to be had

struct AnalyticWavelet <: TransformTypes end
struct RealWaveletRealSignal <: TransformTypes end
struct RealWaveletComplexSignal <: TransformTypes end
struct NonAnalyticMatching <: TransformTypes end

function internalConvFFT(x̂, shears::AbstractArray{<:Number,N}, usedInds,
                         fftPlan, bias, isAnalytic) where N
    axShear = axes(shears)
    axx = axes(x̂)[N:end - 1]
    @views łλ(ii, bias) = applyWeight(x̂, shears[axShear[1:end - 1]..., ii], usedInds, fftPlan, bias[ii], isAnalytic[ii])
    @views łλ(ii, bias::Nothing) = applyWeight(x̂, shears[axShear[1:end - 1]..., ii], usedInds, fftPlan, bias, isAnalytic[ii])
    @views mapped = map(ii -> łλ(ii, bias), 1:size(shears)[end])
    return permutedims(cat(mapped..., dims=1), ((2:N)..., 1, (N + 1):ndims(mapped[1])...))
end

# no bias, not analytic and both match (either both real or both complex)
function applyWeight(x̂, shear, usedInds, fftPlan, bias::Nothing, An::NonAnalyticMatching)
    tmp = fftPlan \ (shear .* x̂) # filter
    @views tmp = tmp[usedInds..., axes(tmp)[length(usedInds) + 1:end]...] # get rid of the padding
    tmp = reshape(tmp, (1, size(tmp)...)) # add a dummy dimension to join over
    return tmp
end

# no bias, analytic wavelet (so complex valued, but only the positive half of x̂
# matters)
function applyWeight(x̂, shear, usedInds, fftPlan, bias::Nothing, An::AnalyticWavelet)
    outer = axes(x̂)[ndims(shear) + 1:end]
    isSourceOdd = mod(size(fftPlan, 1) + 1, 2)
    accessedAxes = axes(shear)
    @views tmp = shear .* x̂[accessedAxes..., outer...] # filter
    wave = cat(tmp, adapt(tmp, zeros(eltype(tmp), size(shear, 1) - 1 - isSourceOdd,
                                     size(tmp)[2:end]...)), dims=1) # symmetrize
    tmp = fftPlan \ wave       # back to time domain
    @views tmp = tmp[usedInds..., axes(tmp)[length(usedInds) + 1:end]...] # get rid of the padding
    tmp = reshape(tmp, (1, size(tmp)...)) # add a dummy dimension to join over
    return tmp
end

# no bias, not analytic, still symmetric (i.e. real, this is for the averaging
# wavelet where the other wavelets are analytic/complex)
function applyWeight(x̂, shear, usedInds, fftPlan, bias::Nothing, An::RealWaveletRealSignal)
    outer = axes(x̂)[ndims(shear) + 1:end]
    isSourceOdd = mod(size(fftPlan, 1) + 1, 2)
    tmp = shear .* x̂ # filter
    wave = cat(tmp, reverse(conj.(tmp[2:end - isSourceOdd, outer...]), dims=1), dims=1) # symmetrize
    tmp = fftPlan \ wave        # back to time domain
    tmp = tmp[usedInds..., axes(tmp)[length(usedInds) + 1:end]...] # get rid of the padding
    tmp = reshape(tmp, (1, size(tmp)...)) # add a dummy dimension to join over
    return tmp
end

# no bias, signal asymmetric/complex, but the wavelet is real (averaging function)
function applyWeight(x̂, shear, usedInds, fftPlan, bias::Nothing, An::RealWaveletComplexSignal)
    outer = axes(x̂)[ndims(shear) + 1:end]
    isSourceOdd = mod(size(fftPlan, 1) + 1, 2)
    tmp = cat(shear, reverse(conj.(shear[2:end - isSourceOdd]),
                             dims=1), dims=1) # construct the full wavelet
    tmp = tmp .* x̂ # filter
    tmp = fftPlan \ tmp        # back to time domain
    tmp = tmp[usedInds..., axes(tmp)[length(usedInds) + 1:end]...] # get rid of the padding
    tmp = reshape(tmp, (1, size(tmp)...)) # add a dummy dimension to join over
    return tmp
end
# biased (and one of the others, doesn't matter which)
function applyWeight(x̂, shear, usedInds, fftPlan, bias, An)
    return applyWeight(x̂, shear, usedInds, fftPlan, nothing, An) .+
        bias
end


text:
@testset "Spine2D" begin
    @polyvar x y
    f = x^2*y + y^2 + 3x^2*y^3 + y^4 + x^4*y^4

    S = spine(f)
    @test length(components_complement(S)) == 5
    @test sprint(show, S) == "Spine2D{DynamicPolynomials.Polynomial{true, Tropical{Float64}}} with:\n * 4 vertices\n * 4 half-rays\n * 1 bounded cells"
    S = spine(f, minimal_component_size=0.05)
    @test length(components_complement(S)) == 5

    S = spine(f, domain=(-5, 5, -5, 5), minimal_component_size=0.01)
    @test length(components_complement(S)) == 5
end


text:
"""
`gamut_vis(type, title, lum, whitepoint, colmatch,name::AbstractString,altitude::Real,slice=collect(5:5:95))`

creates 3D flyaround of gamut given by illuminant

Example: `gamut_vis(Colors.DIN99o,"DIN99o at 6500K", D_series_luminance(6500),D_series_whitepoint(6500,CMF2012_10),CMF2012_10,"DIN99o_6500K_",35,collect(5:5:95))`
"""
function gamut_vis(type, title::AbstractString, lum::LSpec, whitepoint, colmatch, name::AbstractString,altitude::Real,slice=collect(5:5:95))
    println("Preparing gamut... this may take about 30s.")
    a,b,c=sl(type,lum,whitepoint,colmatch,390.0,1.0,830.0)
    println("Creating animation frames... this may take a while.")
    n=name
    if type == xyY
        xl = (0,1)
        yl = (0,1)
        zl = (0,1)
        xg = "x"
        yg = "y"
        zg = "Y"
    elseif type == XYZ
        xl = (0,1.5)
        yl = (0,1.5)
        zl = (0,1.5)
        xg = "X"
        yg = "Z"
        zg = "Y"
    elseif type == Lab
        xl = (-130,130)
        yl = (-130,130)
        zl = (0,100)
        xg = "a"
        yg = "b"
        zg = "L"
    elseif type == Luv
        xl = (-130,130)
        yl = (-130,130)
        zl = (0,100)
        xg = "u"
        yg = "v"
        zg = "L"
    elseif type in (DIN99,DIN99d,DIN99o)
        xl = (-60,60)
        yl = (-60,60)
        zl = (0,100)
        xg = "a99"
        yg = "b99"
        zg = "L99"
    end

    @inbounds for azimuth in 0:359
        gamutslice(type, title, a, b, c,(xl,yl,zl,xg,yg,zg), azimuth, altitude, slice)
        png("$(name)_$(string(azimuth, pad=3))")
        perc=round((azimuth+1)/360*100,digits=2)
        println("Frame $azimuth finished. $perc%")

    end
end

"""
`gamut_vis(type,T1::Real,dT::Real,T2::Real,colmatch,name::AbstractString,altitude::Real,azimuth::Real,slice=collect(5:5:95))`

Creates gamut given by Temperature range
"""
function gamut_vis(type,T1::Real,dT::Real,T2::Real,colmatch,name::AbstractString,altitude::Real,azimuth::Real,slice=collect(5:5:95))
    println("Preparing gamut... this may take about 30s.")
    for T in T1:dT:T2
        lum=normalize_spec(blackbody_illuminant(SPECENV, T))
        wp=blackbody_whitepoint(SPECENV, T)
        a,b,c=sl(type,lum,wp,colmatch,390.0,1.0,830.0)
        println("Creating animation frames... this may take a while.")
        n=name
        if type == xyY
            xl = (0,1)
            yl = (0,1)
            zl = (0,1)
            xg = "x"
            yg = "y"
            zg = "Y"
        elseif type == XYZ
            xl = (0,1.5)
            yl = (0,1.5)
            zl = (0,1.5)
            xg = "X"
            yg = "Z"
            zg = "Y"
        elseif type == Lab
            xl = (-130,130)
            yl = (-130,130)
            zl = (0,100)
            xg = "a"
            yg = "b"
            zg = "L"
        elseif type == Luv
            xl = (-130,130)
            yl = (-130,130)
            zl = (0,100)
            xg = "u"
            yg = "v"
            zg = "L"
        elseif type in (DIN99,DIN99d,DIN99o)
            xl = (-60,60)
            yl = (-60,60)
            zl = (0,100)
            xg = "a99"
            yg = "b99"
            zg = "L99"
        end
        title=T
        gamutslice(type, title, a, b, c,(xl,yl,zl,xg,yg,zg), azimuth, altitude, slice)
        png("$(name)_$(string(tstr, pad=5))")
        println("Frame for CCT $tstr K finished.")
    end
end

function gamutslice(type, title, a, b, c, (xl,yl,zl,xg,yg,zg), azimuth, altitude, slice)
    scatter()
    @inbounds for s in slice
        h,i,j=getslice(a,b,c,s,type)
        plot!(h,i,j,labels=false,color=:black,linewidth=0.5,markersize=0.35,xlims=xl,xguide=xg,ylims=yl,yguide=yg,zlims=zl,zguide=zg,camera=(azimuth,altitude),title="$(title)")
    end
end


text:
################################################################################
# Toy Example
################################################################################
using Plots
T = Float64

# Define the dynamics of the system
p = 2 # Number of players
model = DoubleIntegratorGame(p=p) # game with 3 players with unicycle dynamics
model = UnicycleGame(p=p) # game with 3 players with unicycle dynamics
model = BicycleGame(p=p) # game with 3 players with unicycle dynamics
n = model.n
m = model.m

# Define the horizon of the problem
N = 20 # N time steps
dt = 0.1 # each step lasts 0.1 second
probsize = ProblemSize(N,model) # Structure holding the relevant sizes of the problem

# Define the objective of each player
# We use a LQR cost
Q = [Diagonal(10*ones(SVector{model.ni[i],T})) for i=1:p] # Quadratic state cost
# Q = [Diagonal(10*ones(SVector{model.ni[1],T})), Diagonal(1*ones(SVector{model.ni[2],T}))] # Quadratic state cost
# Q = [Diagonal(1*ones(SVector{model.ni[1],T})), Diagonal(10*ones(SVector{model.ni[2],T}))] # Quadratic state cost
R = [Diagonal(0.1*ones(SVector{model.mi[i],T})) for i=1:p] # Quadratic control cost
# R = [Diagonal(0.1*ones(SVector{model.mi[1],T})), Diagonal(0.01*ones(SVector{model.mi[2],T}))] # Quadratic control cost
# R = [Diagonal(0.01*ones(SVector{model.mi[1],T})), Diagonal(0.1*ones(SVector{model.mi[2],T}))] # Quadratic control cost
# Desrired state
xf = [SVector{model.ni[1],T}([2, 0.0,0,0]),
      SVector{model.ni[2],T}([2, 0.0,0,0]),
      # SVector{model.ni[3],T}([3,-0.4,0,0]),
      # SVector{model.ni[4],T}([3,+0.8,0,0]),
      ]
# Desired control
uf = [zeros(SVector{model.mi[i],T}) for i=1:p]

# # Objectives of the game
game_obj = GameObjective(Q,R,xf,uf,N,model)
# radius = 0.5*ones(p)
# μ = 4.0*ones(p)
# μ = 4.0*[1.0, 0.1]
# μ = 4.0*[0.1, 1.0]
# add_collision_cost!(game_obj, radius, μ)

# Define the constraints that each player must respect
game_con = GameConstraintValues(probsize)
# Add collision avoidance
radius = 0.19
add_collision_avoidance!(game_con, radius)

# # Add control bounds
# u_max =  5*ones(SVector{m,T})
# u_min = -5*ones(SVector{m,T})
# add_control_bound!(game_con, u_max, u_min)

# # Add wall constraint
# walls = [Wall([-1.0,-1.75], [1.0,0.25], [1.,-1.]/sqrt(2))]
# add_wall_constraint!(game_con, walls)

# # Add circle constraint
# xc = [1., 2., 3.]
# yc = [1., 2., 3.]
# radius = [0.1, 0.2, 0.3]
# add_circle_constraint!(game_con, xc, yc, radius)

# Define the initial state of the system
x0 = SVector{model.n,T}([
    0.0, 0.0,# 0.5, #0.0,
   -0.2, 0.2,# 0.7, #0.6,
    0.0, 0.0,# 0.0, #0.0,
    0.0, 0.0,# 0.0, #0.0,
    ])

# Define the Options of the solver
opts = Options()
opts.ls_iter = 15
opts.outer_iter = 20
opts.inner_iter = 20
opts.ρ_0 = 1e0
opts.reg_0 = 1e-7
opts.α_dual = 1.0
opts.λ_max = 1.0*1e7
opts.ϵ_dyn = 1e-6
opts.ϵ_sta = 1e-6
opts.ϵ_con = 1e-6
opts.ϵ_opt = 1e-6
opts.regularize = true
# opts.αx_dual = ones(p)
# opts.αx_dual = [sqrt(10), 1/sqrt(10)]
# opts.αx_dual = [1/sqrt(10), sqrt(10)]

s = 0.5
opts.αx_dual = [sqrt(s), 1/sqrt(s)]
# Define the game problem
prob = GameProblem(N,dt,x0,model,opts,game_obj,game_con)

prob.game_con.state_conval[1][1]
# Solve the problem
@time newton_solve!(prob)
# @profiler newton_solve!(prob)

plot!(prob.model, prob.pdtraj.pr)
plot!(prob.stats)

function subspace_animation(prob::GameProblem)
    p = prob.probsize.p
    S = exp.(Vector(-1.0:0.05:1.0) .* log(10))
    plt = plot(legend=false)

    anim = @animate for s in S
        @show s
        prob.opts.αx_dual = [sqrt(s), 1/sqrt(s)]
        set_constraint_params!(prob.game_con, prob.opts)
        newton_solve!(prob)
        plot!(prob.model, prob.pdtraj.pr, plt=plt)
        text(plt,4, 0.5, "hello")
        # plot!(prob.stats)
    end
    path = joinpath("~/Documents/p2_subspace_din.gif")
    gif(anim, path, fps = 10)
    return nothing
end

subspace_animation(prob)


text:
"""
Code for creating, loading and maniuplating line lists and masks.

Author: Eric Ford
Created: August 2020
"""

"""
    Read line list in ESPRESSO csv format.
ESPRESSO format: lambda and weight.
convert_air_to_vacuum determines whether to convert to vacuum wavelengths.
Warning: ESPRESSO masks don't provide line depth and sometimes include one entry for a blend of lines.
"""
function read_linelist_espresso(fn::String; convert_air_to_vacuum::Bool = true)
    local df = CSV.read(fn,DataFrame,threaded=false,header=["lambda","weight"],delim=' ',ignorerepeated=true)
    @assert hasproperty(df, :lambda)
    @assert hasproperty(df, :weight)
    if convert_air_to_vacuum
        df[!,:lambda] .= λ_air_to_vac.(df[!,:lambda])
    end
    return df
end

""" Read line list in VALD csv format.
   VALD format: lambda_lo, lambdaa_hi and depth.
   convert_air_to_vacuum determines whether to convert to vacuum wavelengths.
"""
function read_linelist_vald(fn::String; convert_air_to_vacuum::Bool = true)
    local df = CSV.read(fn,DataFrame,threaded=false,header=["lambda","depth"])
    @assert hasproperty(df, :lambda)
    @assert hasproperty(df, :depth)
    if convert_air_to_vacuum
        df[!,:lambda] .= λ_air_to_vac.(df[!,:lambda])
    end
    df[!,:weight] = df[!,:depth] # TODO: Decide out what we want to do about tracking depths and weights sepoarately
    return df
end

""" Read line list in csv format.
   format: lambda, weight, lambda_lo, lambdaa_hi.
   Assumes air to vacuumb wavelength conversion has already been applied.
"""
function read_linelist_rvspectml(fn::String)
    local df = CSV.read(fn,DataFrame,threaded=false)
    @assert hasproperty(df, :lambda)
    @assert hasproperty(df, :weight) || hasproperty(df, :depth)
    if hasproperty(df, :depth) && !hasproperty(df, :weight)
        df[!,:weight] = df[!,:depth] # TODO: Decide out what we want to do about tracking depths and weights sepoarately
    end
    return df
end

max_stars_repo_name: RvSpectML/EchelleCCFs

max_issues_repo_name: RvSpectML/EchelleCCFs

max_forks_repo_name: RvSpectML/EchelleCCFs