text:
#=
Copyright 2021 BlackRock, Inc.

Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at

    http://www.apache.org/licenses/LICENSE-2.0

Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
=#

"""
    BoundedQuadratic(lb::Real, ub::Real, p::Real, q::Real, r::Real)

Represents bounded quadratic function
```math
f(x) = px^2 + qx + r, ∀ x ∈ [lb, ub]
```

The fields represent:
- `lb`: Lower bound of the function
- `ub`: Upper bound of the function
- `p`: Coefficient of the quadratic term
- `q`: Coefficient of the linear term
- `r`: Constant

# Example
```jldoctest
julia> bq = BoundedQuadratic(0., 5., 3., 2., 1.)
BoundedQuadratic: f(x) = 3.00000 x² + 2.00000 x + 1.00000, ∀x ∈ [0.00000, 5.00000]

```
"""
mutable struct BoundedQuadratic
    lb::Float64
    ub::Float64
    p::Float64
    q::Float64
    r::Float64

    function BoundedQuadratic(lb::Real, ub::Real, p::Real, q::Real, r::Real)
        @assert !isnan(lb)
        @assert !isnan(ub)
        @assert isfinite(p)
        @assert isfinite(q)
        @assert isfinite(r)
        return new(lb, ub, p, q, r)
    end
end

"""
    BoundedQuadratic()

Constructs an unbounded BoundedQuadratic `f(x) = 0, ∀ x`.
"""
BoundedQuadratic() = BoundedQuadratic(0.0, 0.0, 0.0)

"""
    BoundedQuadratic(p::Real, q::Real, r::Real)

Constructs a unbounded `BoundedQuadratic` function `f(x) = px^2 + qx + r, ∀ x ∈ ℝ`.
"""
BoundedQuadratic(p::Real, q::Real, r::Real) = BoundedQuadratic(Interval(), p, q, r)

"""
    BoundedQuadratic(dom::Interval, p::Real, q::Real, r::Real

Constructs a `BoundedQuadratic` function on an Interval `dom`.
"""
function BoundedQuadratic(dom::Interval, p::Real, q::Real, r::Real)
    return BoundedQuadratic(dom.lb, dom.ub, p, q, r)
end

function Base.show(io::IO, f::BoundedQuadratic)
    @printf(io, "BoundedQuadratic: f(x) = ")
    if f.p > 0
        @printf(io, "%.5f x² ", f.p)
    elseif f.p < 0
        @printf(io, "-%.5f x² ", -f.p)
    end
    if f.q > 0
        @printf(io, "+ %.5f x ", f.q)
    elseif f.q < 0
        @printf(io, "- %.5f x ", -f.q)
    end
    if f.r > 0
        @printf(io, "+ %.5f", f.r)
    elseif f.r < 0
        @printf(io, "- %.5f", -f.r)
    end
    if (f.p == 0) & (f.q == 0)
        @printf(io, "0")
    end
    @printf(io, ", ∀x ∈ ")
    print(io, domain(f))
    return
end

"""
    zero(::Type{BoundedQuadratic})

Construct an empty BoundedQuadratic with no constraints.

# Example
```jldoctest
julia> zero(BoundedQuadratic)
BoundedQuadratic: f(x) = 0, ∀x ∈ ℝ

```
"""
zero(::Type{BoundedQuadratic}) = BoundedQuadratic()
copy(f::BoundedQuadratic) = BoundedQuadratic(f.lb, f.ub, f.p, f.q, f.r)
function copy!(out::BoundedQuadratic, f::BoundedQuadratic)
    out.lb = f.lb
    out.ub = f.ub
    out.p = f.p
    out.q = f.q
    out.r = f.r
    return out
end

"""
    get_tangent(f::BoundedQuadratic, x::Real)

Construct a new unbounded `BoundedQuadratic` representing the (unbounded) tangent line to `f` at the point `x`.

# Example
```jldoctest
julia> bq = BoundedQuadratic(-1., 1., 1., 1., 1.)
BoundedQuadratic: f(x) = 1.00000 x² + 1.00000 x + 1.00000, ∀x ∈ [-1.00000, 1.00000]

julia> tangent = get_tangent(bq, 0.5)
BoundedQuadratic: f(x) = + 2.00000 x + 0.75000, ∀x ∈ ℝ

```
"""
function get_tangent(f::BoundedQuadratic, x::Real)
    q_tan = 2 * f.p * x + f.q
    r_tan = f.p * x^2 + f.q * x + f.r - q_tan * x
    return BoundedQuadratic(Interval(), 0.0, q_tan, r_tan)
end

"""
    get_line(x1::Real, y1::Real, x2::Real, y2::Real)

Construct a new unbounded `BoundedQuadratic` representing a line passing through `(x1, y1)` and `(x2, y2)`.

# Example
```jldoctest
julia> line = get_line(1., 2., 3., 4.)
BoundedQuadratic: f(x) = + 1.00000 x + 1.00000, ∀x ∈ ℝ

```
"""
function get_line(x1::Real, y1::Real, x2::Real, y2::Real)
    @assert x1 != x2
    q = (y2 - y1) / (x2 - x1)
    return BoundedQuadratic(0.0, q, y1 - q * x1)
end

"""
    domain(f::BoundedQuadratic)

Return the Interval on which the BoundedQuadratic is defined.

# Example
```jldoctest
julia> bq = BoundedQuadratic(0., 5., 3., 2., 1.)
BoundedQuadratic: f(x) = 3.00000 x² + 2.00000 x + 1.00000, ∀x ∈ [0.00000, 5.00000]

julia> dom = domain(bq)
[0.00000, 5.00000]

```
"""
domain(f::BoundedQuadratic) = Interval(f.lb, f.ub)

#####
##### Boolean
#####

"""
    isempty(f::BoundedQuadratic)

Return `true` if the domain of the BoundedQuadratic `f` is empty.
"""
isempty(f::BoundedQuadratic) = isempty(domain(f))

function ≤(g::BoundedQuadratic, f::BoundedQuadratic)
    @assert g.p == 0
    return (g.lb ≤ f.lb) & (g.ub ≥ f.lb) & (minimize(f - g)[2] ≥ 0)
end

≥(g::BoundedQuadratic, f::BoundedQuadratic) = f ≤ g

function ≲(g::BoundedQuadratic, f::BoundedQuadratic)
    @assert g.p == 0
    return (g.lb ≲ f.lb) & (g.ub ≳ f.ub) & (minimize(f - g)[2] ≳ 0.0)
end

≳(g::BoundedQuadratic, f::BoundedQuadratic) = f ≲ g

function ≈(f::BoundedQuadratic, g::BoundedQuadratic)
    return ((abs(f.lb - g.lb) ≤ ϵ) | (f.lb == g.lb)) &
           ((abs(f.ub - g.ub) ≤ ϵ) | (f.ub == g.ub)) &
           (abs(f.p - g.p) ≤ ϵ) &
           (abs(f.q - g.q) ≤ ϵ) &
           (abs(f.r - g.r) ≤ ϵ)
end

"""
    is_point(f::BoundedQuadratic)

Return `true` if the BoundedQuadratic is defined only on a single point (`lb == ub`).

See also: [`is_almost_point`](@ref)

# Example
```jldoctest
julia> bq = BoundedQuadratic(0., 5., 3., 2., 1.)
BoundedQuadratic: f(x) = 3.00000 x² + 2.00000 x + 1.00000, ∀x ∈ [0.00000, 5.00000]

julia> is_point(bq)
false

julia> bq = BoundedQuadratic(5., 5., 3., 2., 1.)
BoundedQuadratic: f(x) = 3.00000 x² + 2.00000 x + 1.00000, ∀x ∈ [5.00000, 5.00000]

julia> is_point(bq)
true

```
"""
is_point(f::BoundedQuadratic) = f.lb == f.ub

"""
    is_almost_point(f::BoundedQuadratic)

Return `true` if the BoundedQuadratic is approximately defined only on a single point (`lb ≈ ub`).

See also: [`is_point`](@ref)

# Example
```jldoctest
julia> bq = BoundedQuadratic(5., 5. + 1e-14, 3., 2., 1.)
BoundedQuadratic: f(x) = 3.00000 x² + 2.00000 x + 1.00000, ∀x ∈ [5.00000, 5.00000]

julia> is_point(bq)
false

julia> is_almost_point(bq)
true

```
"""
function is_almost_point(f::BoundedQuadratic)
    return (is_point(f)) || (f.lb ≈ f.ub)
end

"""
    continuous_and_overlapping(f::BoundedQuadratic, g::BoundedQuadratic)

Return `true` if `f`'s right endpoint corresponds with `g`'s left endpoint.
"""
function continuous_and_overlapping(f::BoundedQuadratic, g::BoundedQuadratic)
    return (f.ub ≈ g.lb) && (f(f.ub) ≈ g(g.lb))
end

is_convex(f::BoundedQuadratic) = f.p >= 0.0

#####
##### Operations
#####

"""
    (f::BoundedQuadratic)(x::Real)

Evaluate `f(x)` if `x` is in the domain of `f`, else return `Inf`.

# Example
```jldoctest
julia> bq = BoundedQuadratic(3., 2., 1.)
BoundedQuadratic: f(x) = 3.00000 x² + 2.00000 x + 1.00000, ∀x ∈ ℝ

julia> bq(1.)
6.0

```
"""
function (f::BoundedQuadratic)(x::Real)
    if x in domain(f)
        return f.p * x^2 + f.q * x + f.r
    else
        return Inf
    end
end

function +(f::BoundedQuadratic, a::Real)
    return BoundedQuadratic(f.lb, f.ub, f.p, f.q, f.r + a)
end

function -(f::BoundedQuadratic)
    @assert f.p == 0.0
    return BoundedQuadratic(f.lb, f.ub, -f.p, -f.q, -f.r)
end

function +(f::BoundedQuadratic, g::BoundedQuadratic)
    dom = domain(f) ∩ domain(g)
    return BoundedQuadratic(dom, f.p + g.p, f.q + g.q, f.r + g.r)
end

function -(f::BoundedQuadratic, g::BoundedQuadratic)
    return f + (-g)
end

*(f::BoundedQuadratic, α::Real) = BoundedQuadratic(f.lb, f.ub, α * f.p, α * f.q, α * f.r)
*(α::Real, f::BoundedQuadratic) = f * α

"""
    scale(f::BoundedQuadratic, α::Real)

Return a new `BoundedQuadratic` that has been scaled by `α`. That is, given `f(x)` and `α`, returns `f(αx)`.

Note: this operation requires scaling the domain.

See also: [`scale!`](@ref)

# Example
```jldoctest
julia> bq = BoundedQuadratic(-1., 1., 1., 1., 1.)
BoundedQuadratic: f(x) = 1.00000 x² + 1.00000 x + 1.00000, ∀x ∈ [-1.00000, 1.00000]

julia> scale(bq, 2.)
BoundedQuadratic: f(x) = 4.00000 x² + 2.00000 x + 1.00000, ∀x ∈ [-0.50000, 0.50000]

```
"""
function scale(f::BoundedQuadratic, α::Real)
    return BoundedQuadratic(f.lb / α, f.ub / α, α^2 * f.p, α * f.q, f.r)
end

"""
     scale!(f::BoundedQuadratic, α::Real)
     scale!(f::BoundedQuadratic, α::Real, out::BoundedQuadratic)

Scale `f` inplace.

See also: [`scale`](@ref)
"""
scale!(f::BoundedQuadratic, α::Real) = scale!(f, α, f)
scale!(f::BoundedQuadratic, α::Real, out::BoundedQuadratic) = copy!(out, scale(f, α))

"""
    perspective(f::BoundedQuadratic, α::Real)

Return the perspective function of `f`. That is, given `f(x)` and `α`, return `α * f(x / α)`.

Note: that this operation requires scaling of the domain.

See also: [`perspective!`](@ref)

# Example
```jldoctest
julia> bq = BoundedQuadratic(-1., 1., 1., 1., 1.)
BoundedQuadratic: f(x) = 1.00000 x² + 1.00000 x + 1.00000, ∀x ∈ [-1.00000, 1.00000]

julia> perspective(bq, 2.)
BoundedQuadratic: f(x) = 0.50000 x² + 1.00000 x + 2.00000, ∀x ∈ [-2.00000, 2.00000]

```
"""
function perspective(f::BoundedQuadratic, α::Real)
    return BoundedQuadratic(α * f.lb, α * f.ub, f.p / α, f.q, f.r * α)
end

"""
    perspective!(f::BoundedQuadratic, α::Real)
    perspective!(f::BoundedQuadratic, α::Real, out::BoundedQuadratic)

Shift perspective of `f` inplace.

See also: [`perspective`](@ref)
"""
perspective!(f::BoundedQuadratic, α::Real) = perspective!(f, α, f)
function perspective!(f::BoundedQuadratic, α::Real, out::BoundedQuadratic)
    return copy!(out, perspective(f, α))
end

"""
    shift(f::BoundedQuadratic, δ::Real)

Return `f` shifted along the `x`-axis by `δ`.

Note: for `δ > 0`, this is a right shift. For `δ < 0`, this is a left shift.

See also: [`shift!`](@ref)

# Example
```jldoctest
julia> bq = BoundedQuadratic(-1., 1., 1., 1., 1.)
BoundedQuadratic: f(x) = 1.00000 x² + 1.00000 x + 1.00000, ∀x ∈ [-1.00000, 1.00000]

julia> shift(bq, 2.)
BoundedQuadratic: f(x) = 1.00000 x² - 3.00000 x + 3.00000, ∀x ∈ [1.00000, 3.00000]

```
"""
function shift(f::BoundedQuadratic, δ::Real)
    return BoundedQuadratic(f.lb + δ, f.ub + δ, f.p, f.q - 2 * f.p * δ,
                            f.p * δ^2 - f.q * δ + f.r)
end

"""
    shift!(f::BoundedQuadratic, δ::Real)
    shift!(f::BoundedQuadratic, δ::Real, out::BoundedQuadratic)

Shift `f` inplace along the `x`-axis by `δ`.

See also: [`shift`](@ref)
"""
shift!(f::BoundedQuadratic, δ::Real) = shift!(f, δ, f)
function shift!(f::BoundedQuadratic, δ::Real, out::BoundedQuadratic)
    return copy!(out, shift(f, δ))
end

"""
    tilt(f::BoundedQuadratic, α::Real)

Return `f` tilted by `α`. This shifts linear coefficient `q` by `α`.

See also: [`tilt!`](@ref)

# Example
```jldoctest
julia> bq = BoundedQuadratic(-1., 1., 1., 1., 1.)
BoundedQuadratic: f(x) = 1.00000 x² + 1.00000 x + 1.00000, ∀x ∈ [-1.00000, 1.00000]

julia> tilt(bq, 2.)
BoundedQuadratic: f(x) = 1.00000 x² + 3.00000 x + 1.00000, ∀x ∈ [-1.00000, 1.00000]

```
"""
function tilt(f::BoundedQuadratic, α::Real)
    return BoundedQuadratic(f.lb, f.ub, f.p, f.q + α, f.r)
end

"""
    tilt!(f::BoundedQuadratic, α::Real)
    tilt!(f::BoundedQuadratic, α::Real, out::BoundedQuadratic)

Tilt `f` inplace.

See also: [`tilt`](@ref)
"""
tilt!(f::BoundedQuadratic, α::Real) = tilt!(f, α, f)
function tilt!(f::BoundedQuadratic, α::Real, out::BoundedQuadratic)
    return copy!(out, tilt(f, α))
end

"""
    restrict_dom(f::BoundedQuadratic, dom::Interval)
    restrict_dom(f::BoundedQuadratic, lb::Real, ub::Real)

Return a new BoundedQuadratic with domain restricted to the intersect of the passed domain.

See also: [`restrict_dom!`](@ref)

# Example
```jldoctest
julia> bq = BoundedQuadratic(-10., 10., 1., 1., 1.)
BoundedQuadratic: f(x) = 1.00000 x² + 1.00000 x + 1.00000, ∀x ∈ [-10.00000, 10.00000]

julia> restrict_dom(bq, 2., 3.)
BoundedQuadratic: f(x) = 1.00000 x² + 1.00000 x + 1.00000, ∀x ∈ [2.00000, 3.00000]

```
"""
function restrict_dom(f::BoundedQuadratic, dom::Interval)
    restricted_dom = domain(f) ∩ dom
    @assert restricted_dom.ub ≥ restricted_dom.lb
    return BoundedQuadratic(restricted_dom, f.p, f.q, f.r)
end
restrict_dom(f::BoundedQuadratic, lb::Real, ub::Real) = restrict_dom(f, Interval(lb, ub))

"""
    restrict_dom!(f::BoundedQuadratic, dom::Interval)
    restrict_dom!(f::BoundedQuadratic, lb::Real, ub::Real)
    restrict_dom!(f::BoundedQuadratic, lb::Real, ub::Real, out::BoundedQuadratic)

Restrict domain of `f` inplace.

See also: [`restrict_dom`](@ref)
"""
restrict_dom!(f::BoundedQuadratic, dom::Interval) = restrict_dom!(f, dom, f)
function restrict_dom!(f::BoundedQuadratic, dom::Interval, out::BoundedQuadratic)
    return copy!(out, restrict_dom(f, dom))
end
restrict_dom!(f::BoundedQuadratic, lb::Real, ub::Real) = restrict_dom!(f, lb, ub, f)
function restrict_dom!(f::BoundedQuadratic, lb::Real, ub::Real, out::BoundedQuadratic)
    return restrict_dom!(f, Interval(lb, ub), f)
end

"""
    extend_dom(f::BoundedQuadratic)

Return an unbounded version of `f`.

See also: [`extend_dom!`](@ref)

# Example
```jldoctest
julia> bq = BoundedQuadratic(-1., 1., 1., 1., 1.)
BoundedQuadratic: f(x) = 1.00000 x² + 1.00000 x + 1.00000, ∀x ∈ [-1.00000, 1.00000]

julia> extend_dom(bq)
BoundedQuadratic: f(x) = 1.00000 x² + 1.00000 x + 1.00000, ∀x ∈ ℝ

```
"""
extend_dom(f::BoundedQuadratic) = BoundedQuadratic(f.p, f.q, f.r)

"""
    extend_dom!(f::BoundedQuadratic)
    extend_dom!(f::BoundedQuadratic, out::BoundedQuadratic)

Remove the bounds of `f` inplace.

See also: [`extend_dom`](@ref)
"""
extend_dom!(f::BoundedQuadratic) = extend_dom!(f, f)
extend_dom!(f::BoundedQuadratic, out::BoundedQuadratic) = copy!(out, extend_dom(f))

"""
    reverse(f::BoundedQuadratic)

Return `f` reversed over the `y` axis. That is, given `f(x)`, return `f(-x)`.

See also: [`reverse!`](@ref)

# Example
```jldoctest
julia> bq = BoundedQuadratic(-1., 2., 1., 1., 1.)
BoundedQuadratic: f(x) = 1.00000 x² + 1.00000 x + 1.00000, ∀x ∈ [-1.00000, 2.00000]

julia> reverse(bq)
BoundedQuadratic: f(x) = 1.00000 x² - 1.00000 x + 1.00000, ∀x ∈ [-2.00000, 1.00000]

```
"""
reverse(f::BoundedQuadratic) = BoundedQuadratic(-f.ub, -f.lb, f.p, -f.q, f.r)

"""
    reverse!(f::BoundedQuadratic)
    reverse!(f::BoundedQuadratic, out::BoundedQuadratic)

Reverse `f` inplace.

See also: [`reverse`](@ref)
"""
reverse!(f::BoundedQuadratic) = reverse!(f, f)
function reverse!(f::BoundedQuadratic, out::BoundedQuadratic)
    return copy!(out, reverse(f))
end

"""
    derivative(f::BoundedQuadratic)

Return the derivative of `f`, `f'`

# Example
```jldoctest
julia> bq = BoundedQuadratic(3., 2., 1.)
BoundedQuadratic: f(x) = 3.00000 x² + 2.00000 x + 1.00000, ∀x ∈ ℝ

julia> derivative(bq)
BoundedQuadratic: f(x) = + 6.00000 x + 2.00000, ∀x ∈ ℝ

```
"""
function derivative(f::BoundedQuadratic)
    return BoundedQuadratic(f.lb, f.ub, 0.0, 2.0 * f.p, f.q)
end

"""
    derivative(f::BoundedQuadratic, x::Real)

Evaluate the derivative of `f` at `x`, `f'(x)`.

# Example
```jldoctest
julia> bq = BoundedQuadratic(3., 2., 1.)
BoundedQuadratic: f(x) = 3.00000 x² + 2.00000 x + 1.00000, ∀x ∈ ℝ

julia> derivative(bq, 2.)
14.0

```
"""
derivative(f::BoundedQuadratic, x::Real) = derivative(f)(x)

function _sum(f_list::Vector{BoundedQuadratic})
    dom, is_valid = _intersect(f_list)

    p = 0.0
    q = 0.0
    r = 0.0
    if is_valid
        for fi in f_list
            p += fi.p
            q += fi.q
            r += fi.r
        end
        return BoundedQuadratic(dom, p, q, r), is_valid
    else
        return BoundedQuadratic(dom, 0.0, 0.0, 0.0), is_valid
    end
end

"""
    sum(f_list::Vector{BoundedQuadratic})

Return the BoundedQuadratic sum of a list of BoundedQuadratics if it is valid, else `nothing`.

# Example
```jldoctest
julia> bq = BoundedQuadratic(-1., 5., 1., 2., 3.)
BoundedQuadratic: f(x) = 1.00000 x² + 2.00000 x + 3.00000, ∀x ∈ [-1.00000, 5.00000]

julia> sum([bq, bq])
BoundedQuadratic: f(x) = 2.00000 x² + 4.00000 x + 6.00000, ∀x ∈ [-1.00000, 5.00000]

```
"""
function sum(f_list::Vector{BoundedQuadratic})
    g, is_valid = _sum(f_list)
    if is_valid
        return g
    else
        return nothing
    end
end

function _intersect(f_list::Vector{BoundedQuadratic})
    k = length(f_list)
    doms = [Interval(), map(domain, f_list)...]
    dom = reduce(∩, doms)
    is_valid = !isempty(dom)
    return dom, is_valid
end

# TODO: the following example code fails as a doctest with Julia 1.0...
# it appears that the output of `Vector{T}` displays as `Array{T,1}` instead
"""
    intersect(f_list::Vector{BoundedQuadratic})

Intersect the domains of a list of BoundedQuadratics if possible.

# Example
```julia
julia> f_list = [BoundedQuadratic(-100, 100, 1, 1, 1),
          BoundedQuadratic(-50, 25, 1, 1, 1),
          BoundedQuadratic(0, 50, 1, 1, 1)]
3-element Vector{BoundedQuadratic}:
 BoundedQuadratic: f(x) = 1.00000 x² + 1.00000 x + 1.00000, ∀x ∈ [-100.00000, 100.00000]

 BoundedQuadratic: f(x) = 1.00000 x² + 1.00000 x + 1.00000, ∀x ∈ [-50.00000, 25.00000]

 BoundedQuadratic: f(x) = 1.00000 x² + 1.00000 x + 1.00000, ∀x ∈ [0.00000, 50.00000]

julia> out, valid = intersect(f_list);

julia> valid
true

julia> out
3-element Vector{BoundedQuadratic}:
 BoundedQuadratic: f(x) = 1.00000 x² + 1.00000 x + 1.00000, ∀x ∈ [0.00000, 25.00000]

 BoundedQuadratic: f(x) = 1.00000 x² + 1.00000 x + 1.00000, ∀x ∈ [0.00000, 25.00000]

 BoundedQuadratic: f(x) = 1.00000 x² + 1.00000 x + 1.00000, ∀x ∈ [0.00000, 25.00000]

julia> f_list = [BoundedQuadratic(-100, 0, 1,1,1)
            BoundedQuadratic(25, 50, 1,1,1)]  # non-overlapping
2-element Vector{BoundedQuadratic}:
 BoundedQuadratic: f(x) = 1.00000 x² + 1.00000 x + 1.00000, ∀x ∈ [-100.00000, 0.00000]

 BoundedQuadratic: f(x) = 1.00000 x² + 1.00000 x + 1.00000, ∀x ∈ [25.00000, 50.00000]

julia> out, valid = intersect(f_list);

julia> valid
false

julia> out
2-element Vector{BoundedQuadratic}:
 #undef
 #undef

```
"""
function intersect(f_list::Vector{BoundedQuadratic})
    k = length(f_list)
    dom, is_valid = _intersect(f_list)

    g_list = Vector{BoundedQuadratic}(undef, k)
    if is_valid
        for (i, f) in enumerate(f_list)
            g_list[i] = restrict_dom(f, dom)
        end
    end
    return g_list, is_valid
end

"""
    minimize(f::BoundedQuadratic)

Return the minimum `x` and `f(x)` of `f` over its domain.
"""
function minimize(f::BoundedQuadratic)
    if f.lb > f.ub  # Infeasible.
        xstar = NaN
        min_val = Inf
    else  # Feasible, but possibly unbounded.
        if f.p > 0
            xstar = -f.q / (2 * f.p)
            if xstar < f.lb
                xstar = f.lb
            elseif xstar > f.ub
                xstar = f.ub
            end
        elseif f.q > 0
            xstar = (f.lb == -Inf) ? NaN : f.lb
        elseif f.q < 0
            xstar = (f.ub == Inf) ? NaN : f.ub
        else  # Case of p = 0, q = 0.
            if f.lb > -Inf
                xstar = f.lb
            elseif f.ub < Inf
                xstar = f.ub
            else
                xstar = 0.0
            end
        end
        if isfinite(xstar)
            min_val = f.p * xstar^2 + f.q * xstar + f.r
        else
            min_val = -Inf
        end
    end
    return xstar, min_val
end


text:
module AtomicData

using Unitful
@unit Ry "Ry" Rydberg 13.605_693_009u"eV" false;
@unit Ha "Ha" Hartree 27.211_386_02u"eV" false;

# This __init__ etc. is necessary for Unitful for the units
# to be fully available when this package is loaded.
# (http://ajkeller34.github.io/Unitful.jl/stable/extending/)
const _local_Unitful_basefactors = Unitful.basefactors
function __init__()
    # To use the atomic units outside this module
    merge!(Unitful.basefactors, _local_Unitful_basefactors)
    Unitful.register(AtomicData)
end
Unitful.register(AtomicData)

using DataFrames
using CSV
using HTTP

function parse_J(J::AbstractVector)
    function parse_number(j::AbstractString)
        if occursin("/", j)
            //(parse.(Ref(Int), split(j, "/"))...)
        else
            parse(Int, j)
        end
    end
    map(J) do j
        if ismissing(j) || j == "---"
            missing
        elseif occursin(",", j)
            parse_number.(split(j, ","))
        else
            parse_number(j)
        end
    end
end

function parse_eng(E::AbstractVector{<:Union{Missing,String}}, unit)
    map(E) do EE
        if !ismissing(EE)
            for c in " ()[]"
                EE = replace(EE, c=>"")
            end
            v = parse(Float64, EE)
            if unit == u"Ha"
                v*u"Ry" |> u"Ha"
            else
                v*unit
            end
        else
            EE
        end
    end |> Vector{Union{Missing,Quantity}}
end

function parse_eng(E::AbstractVector{<:Union{Missing,<:Real}}, unit)
    map(E) do EE
        if !ismissing(EE)
            if unit == u"Ha"
                EE*u"Ry" |> u"Ha"
            else
                EE*unit
            end
        else
            EE
        end
    end |> Vector{Union{Missing,Quantity}}
end

function get_nist_data(io::IO, unit)
    df = io |> f -> CSV.File(f, delim='\t') |> DataFrame

    [df[:, 1:2] DataFrame(J = parse_J(df[!, 3]), Level = parse_eng(df[!, 4], unit), Uncertainty = parse_eng(df[!, 5], unit)) df[:, 6:end]]
end

function get_nist_data(name::String, unit)
    units = Dict(u"cm"^(-1) => 0,
                 u"eV" => 1,
                 u"Ry" => 2,
                 u"Ha" => 2
                 )
    http_name = replace(name, " " => "+")
    url = "https://physics.nist.gov/cgi-bin/ASD/energy1.pl?encodedlist=XXT2&de=0&spectrum=$(http_name)&units=$(units[unit])&upper_limit=&parity_limit=both&conf_limit=All&conf_limit_begin=&conf_limit_end=&term_limit=All&term_limit_begin=&term_limit_end=&J_limit=&format=3&output=0&page_size=15&multiplet_ordered=0&conf_out=on&term_out=on&level_out=on&unc_out=on&j_out=on&lande_out=on&perc_out=on&biblio=on&temp=&submit=Retrieve+Data"

    req_body = HTTP.request("GET", url).body |> String
    occursin("Invalid element symbol", req_body) && error("Failed to retrieve atomic data for $(name)")
    req_body |> IOBuffer |> io -> get_nist_data(io, unit)
end

macro get_nist_data(varname, name, args...)
    varstr = string(varname)
    glob = Expr(:global, varname)
    quote
        if !isdefined(Main, Symbol($varstr))
            $glob
            $(esc(varname)) = get_nist_data($name, $(args...))
        end
    end
end

export get_nist_data, @get_nist_data

include("latex_tables.jl")

end # module

ext: jl


text:
struct Args{T,S}
    args::T
    kwargs::S
end

@inline function callwithargs(a::Args{T,S}, c::C, x::X) where {T,S,C,X}
    c(x, a.args...; a.kwargs...)
end
@inline function callwithargs(a::Args{T,Nothing}, c::C, x::X) where {T,C,X}
    c(x, a.args...)
end
@inline function callwithargs(a::Args{Nothing,S}, c::C, x::X) where {S,C,X}
    c(x; a.kwargs...)
end
@inline function callwithargs(a::Args{Nothing,Nothing}, c::C, x::X) where {C,X}
    c(x)
end

function tuplecast(f::F, x::NamedTuple{names,types}, a...) where {F,names,types}
    o = f.(values(x), a...)
    NamedTuple{names, typeof(o)}(o)
end
tuplecast(f::F, x::Tuple, a...) where F = f.(x, a...)
tuplecast(f::F, x::Nothing, a...) where F = nothing

function substituteargs(x::AbstractSelection, df)
    k = getselectionkeys(df, x)
    if length(k) > 1
        throw(error("Selection $(x) returned $(length(k)) columns $(k), instead of just one."))
    elseif length(k) == 0
        throw(error("Selection $(x) returned no columns."))
    else
        df[k[1]]
    end
end
substituteargs(x, df) = x
substituteargs(x::Nothing, df) = nothing

function selections_to_columns(t::Args{T,S}, df) where {T,S}
    Args(tuplecast(substituteargs, t.args, Ref(df)),
         tuplecast(substituteargs, t.kwargs, Ref(df)))
end

mutable struct StepStructor{S,F,G,A,B} <: AbstractStep
    selections::S
    params::Union{Nothing, Dict{Symbol,G}}
    t::Type{G}
    broadcast::Bool
    trained::Bool
    skip::Bool
    prehook::F
    args::Args{A,B}
end

function transform!(s::StepStructor, df)
    if s.broadcast
        for (k,f) in s.params
            df[k] = f.(df[k])
        end
    else
        for (k,f) in s.params
            df[k] = f(df[k])
        end
    end
    df
end

function transform(s::StepStructor, df)
    transform!(s, copy(df))
end

function fit!(s::StepStructor, df)
    args = selections_to_columns(s.args, df)
    s.params = Dict(col => callwithargs(args, s.t, s.prehook(df[col])) for col in getselectionkeys(df, s.selections))
    s.trained = true
end

function step_struct!(r::Recipe, s...; t=nothing, broadcast::Bool=true,
    skip::Bool=false, prehook=identity, args::Union{Tuple,Nothing}=nothing,
    kwargs::Union{NamedTuple,Nothing}=nothing, names::Union{Vector{Symbol},Nothing}=nothing,
    keep::Bool=true)
    t === nothing && throw(error("No struct supplied to `step_struct!`"))
    push!(r.steps,
          StepStructor([s...],
                       nothing,
                       t,
                       broadcast,
                       false,
                       skip,
                       prehook,
                       Args(args,kwargs)))
end


text:
using .ZeroOne

@testset "zeroone" begin
    @test zeroone("") == 0
    @test zeroone("01") == 0
    @test zeroone("10") == 1
    @test zeroone("1?0") == 4
    @test zeroone("0101?"^2) == 41
    @test zeroone("1?0?") == zeroone("1000") + zeroone("1001") + zeroone("1100") + zeroone("1101")
    @test zeroone.("?" .^ (1:10)) == [0, 1, 6, 24, 80, 240, 672, 1792, 4608, 11520]
    @test zeroone("?" ^ 23) == 530579456
    @test zeroone("?" ^ 24) == 157627897
    @test zeroone("?" ^ 500000) == 879720014
end


text:
using SSNVariability ; const S=SSNVariability
using Calculus, Distributions, LinearAlgebra, Statistics
using Random
using Test

Random.seed!(0)

@testset "input-output functions" begin
  g1=S.ReLu(1.2)
  @test g1(3.3) ≈ 3.3*1.2
  g2=S.ReQuad(4.4)
  @test g2(3.3) ≈ 3.3*3.3*4.4
  # gradients
  x=map( x->x+0.01*sign(x), rand(100) .- 0.2)
  for g in (g1,g2)
    @test g(-0.01) ≈ 0.0
    an = S.ioprime.(x,g)
    num = Calculus.gradient.(g,x)
    @test all( isapprox.(an,num;atol=0.01) )
  end
end

# Jacobian !
@testset "Jacobian matrix" begin
  ne,ni = 13,10
  ntot = ne+ni
  ntw = S.RecurrentNeuralNetwork(ne,ni)
  veli(u,i) = getindex(S.du_nonoise(u,ntw),i)
  @info "building the Jacobian numerically"
  utest = randn(ntot)
  Jnum = Matrix{Float64}(undef,ntot,ntot)
  for i in 1:ntot
    Jnum[i,:] =  Calculus.gradient( u -> veli(u,i), utest )
  end
  Jan=S.jacobian(utest,ntw)
  @test all(isapprox.(Jnum,Jan;rtol=1E-3))
end

@testset "2D Linear model Vs analytic results" begin

  ntw = let
    iofun=S.IOIdentity()
    weight_matrix = [0.5 -0.5 ; 2.5 -0.1]
    time_membrane = [ 0.8 , 0.4]
    input = [10. , 6.0]
    s1,s2,ρ = 1.0,0.3,0.5
    sigma_noise = Symmetric([ s1  s1*s2*ρ ; s1*s2*ρ s2 ])
    @assert isposdef(sigma_noise)
    S.RecurrentNeuralNetwork(iofun,weight_matrix,time_membrane,input,sigma_noise)
  end
  mu_an,sigma_an = S.mean_cov_linear_ntw(ntw)

  r_start = mu_an
  dt = 10E-3
  Ttot = 1000.0
  t,_,ei=S.run_network_noise(ntw,r_start,Ttot;stepsize=dt)
  sigma_num = cov(ei;dims=2)
  mu_num = mean(ei;dims=2)[:]

  @test all(isapprox.(mu_num,mu_an;atol=0.3))
  @test all(isapprox.(sigma_num,sigma_an;atol=0.4))

end

@testset "Find mean with `true` variance computed from simulation" begin
  ne,ni = 43,22
  ntot = ne+ni
  cov_noise_test=S.random_covariance_matrix(ntot,3.3)
  ntw = S.RecurrentNeuralNetwork(ne,ni;sigma_noise=cov_noise_test)

  # more weight variability
  wscal = rand(Uniform(0.5,1.5),S.n_neurons(ntw))
  broadcast!(*,ntw.weight_matrix,ntw.weight_matrix,wscal)
  mu_out,_ = S.mustart_sigmastart(ntw)
  r_sim_start = ntw.iofun.(mu_out)

  t,ev_v,ei=S.run_network_noise(ntw,r_sim_start,40.;stepsize=0.01)
  r_sim_out = median(ei;dims=2)[:]
  cov_diag_an = diag(cov(ev_v;dims=2))

  dmu_alloc = similar(mu_out)
  ε = 1E-3
  for i in 1:5000
    S.dmu!(dmu_alloc,mu_out,cov_diag_an,ntw)
    mu_out .+= ε*dmu_alloc
  end
  r_an_out = ntw.iofun.(mu_out)
  @test all(isapprox.(r_an_out,r_sim_out;atol=0.25))
end

@testset "Mean and covariance matrix with diagonal noise" begin

  Ttot = 300.0
  dt = 0.001
  ne,ni = 22,12
  ntot = ne+ni
  cov_noise_test = Symmetric(diagm(rand(Uniform(6.2,10.3),ntot)))
  ntw = S.RecurrentNeuralNetwork(ne,ni;sigma_noise=cov_noise_test)

  # more weight variability
  wscal = rand(Uniform(1.7,3.3),S.n_neurons(ntw))
  broadcast!(*,ntw.weight_matrix,ntw.weight_matrix,wscal)

  mu_an,sigma_an = S.mustart_sigmastart(ntw)
  r_start = ntw.iofun.(mu_an)
  t,ei_v,ei=S.run_network_noise(ntw,r_start,Ttot;stepsize=dt)
  sigma_num = cov(ei_v;dims=2)
  mu_num = mean(ei_v;dims=2)[:]

  dSigma_alloc = similar(sigma_num)
  dmu_alloc = similar(mu_num)

  ε = 1E-3
  for i in 1:5_000
    S.dmu!(dmu_alloc,mu_an,diag(sigma_an),ntw)
    mu_an .+= ε .* dmu_alloc
    S.dSigma!(dSigma_alloc,sigma_an,mu_an,ntw)
    sigma_an .+= ε .* dSigma_alloc
  end
  @test all(isapprox.(mu_an,mu_num;atol=0.05))
  @test all(isapprox.(sigma_an,sigma_num;atol=0.05))
end

@testset "Mean and covariance with non-diagonal noise" begin

  Ttot = 100.0
  dt = 0.05E-3
  ne,ni = 22,11
  ntot = ne+ni
  cov_noise_test = let c = S.random_covariance_matrix(ntot,1.0)
    rr = Diagonal(sqrt.(rand(Uniform(2.3,8.3),ntot)))
    Symmetric(rr*c*rr)
  end
  ntw = S.RecurrentNeuralNetwork(ne,ni;sigma_noise=cov_noise_test)

  wscal = rand(Uniform(1.7,2.0),S.n_neurons(ntw))
  broadcast!(*,ntw.weight_matrix,ntw.weight_matrix,wscal)

  mu_an,sigma_an = S.mustart_sigmastart(ntw)
  r_start = ntw.iofun.(mu_an)
  t,ei_v,ei=S.run_network_noise(ntw,r_start,Ttot;stepsize=dt)
  sigma_num = cov(ei_v[:,1:100:end];dims=2)
  mu_num = mean(ei_v;dims=2)[:]

  dSigma_alloc = similar(sigma_num)
  dmu_alloc = similar(mu_num)

  ε = 0.5E-3
  for i in 1:10_000
    S.dmu!(dmu_alloc,mu_an,diag(sigma_an),ntw)
    mu_an .+= ε .* dmu_alloc
    S.dSigma!(dSigma_alloc,sigma_an,mu_an,ntw)
    sigma_an .+= ε .* dSigma_alloc
  end
  @test all(isapprox.(mu_an,mu_num;atol=0.05))
  @test all(isapprox.(sigma_an,sigma_num;atol=0.05))
end


text:
mutable struct brookstone_supertype_t <: LCMType
  utime::Int64

  newvariable::generic_variable_t
  newfactor::generic_factor_t
  img::image_t
end

@lcmtypesetup(brookstone_supertype_t)


text:
#=
renormalize!

getTensorM
truncMatrixU
matU2tenU
getNewTensorT
retateTensorT90
normalize
=#

function renormalize!{T}(lattice::Classical2dSquareLattice{T}, dimM::Int)
	tensorT, impurityTensorTtilde = getTensorT(lattice)

	tenmatM = getTensorM(tensorT)
	matrixMMd = tenmatM.matrix * transpose(tenmatM.matrix)
	Ul, lambdaVector, Uld = svd(matrixMMd)
	trunUl = truncMatrixU(Ul,dimM)
	tensorUl = matU2tenU(trunUl)

	# contract tensors to get the bigger tensorT
	newTensorT = getNewTensorT(tensorT,tensorUl)

	# rotate tensorT 90 degrees to the clockwise
	newRotTensorT = rotateTensorT90(newTensorT)

	# normalize tensorT
	newRotNormTensorT, normalizationFactor = normalize(newRotTensorT)

	### similar thing to do to tensorTtilde
	newTensorTtilde = getNewTensorT(tensorT,tensorUl,impurityTensorTtilde)
	newRotTensorTtilde = rotateTensorT90(newTensorTtilde)
	newRotNormTensorTtilde = newRotTensorTtilde ./normalizationFactor

	setTensorT!(lattice, newRotNormTensorT)
	setTensorT!(lattice, newRotNormTensorTtilde; tilde = true)
	return normalizationFactor
end

function getTensorM{T}(tensorT::Array{T,4} ; returntenmat = true, mirror = false)
	if mirror == false
		@tensor begin
			tensorM[x1,x2,x1p,x2p,y,yp] := tensorT[x1,x1p,y,j] * tensorT[x2,x2p,j,yp]
		end
	else
		@tensor begin
			tensorM[x1,x2p,x1p,x2,y,yp] := tensorT[x1,x1p,y,j] * tensorT[x2,x2p,yp,j]
		end
	end
	returntenmat ? answer = tensor2tenmat(tensorM,[1,2],[3,4,5,6]) : answer = tensorM
	return answer
end

function truncMatrixU(matrixU::Matrix, dimM::Int)
	# truncate the column dimension of matrix U to the dimM
	if size(matrixU,2) > dimM
		matrixU = Array(view(matrixU,:,1:dimM))
	end
	return matrixU
end

function matU2tenU(matrixU::Matrix)
	size1 = Int(sqrt(size(matrixU,1)))
	tenmatU = Tenmat(matrixU,[1,2],[3],(size1,size1,size(matrixU,2)))
	tensorU = tenmat2tensor(tenmatU)
	return tensorU
end

function getNewTensorT{T}(tensorT::Array{T,4},tensorU::Array{T,3}, tensorTtilde::Array{T,4}= tensorT; mirror = false)
	# contract U,T,T,U tensors to make new tensorT
	# according to the Fig.1.(b)
	# the third input tensorTtilde is for impurity tensor
	# 	without the third input, it uses tensorT twice
	if mirror == false
		@tensor begin
			UT[x2,x,x1p,y,i] := tensorU[x1,x2,x] * tensorTtilde[x1,x1p,y,i]
			UTT[x,x1p,y,x2p,yp] := UT[x2,x,x1p,y,i] * tensorT[x2,x2p,i,yp]
			newTensorT[x,xp,y,yp] := UTT[x,x1p,y,x2p,yp] * tensorU[x1p,x2p,xp]
		end
	else
		@tensor begin
			UT[x2,x,x1p,y,i] := tensorU[x1,x2,x] * tensorTtilde[x1,x1p,y,i]
			UTT[x,x1p,y,x2p,yp] := UT[x2,x,x1p,y,i] * tensorT[x2,x2p,yp,i]
			newTensorT[x,xp,y,yp] := UTT[x,x1p,y,x2p,yp] * tensorU[x1p,x2p,xp]
		end
	end
	return newTensorT
end


text:
export ft_getbitmode
export ftbitmodedict

ftbitmodedict = Dict(
  FT_BITMODE_RESET => "FT_BITMODE_RESET",
  FT_BITMODE_ASYNC_BITBANG => "FT_BITMODE_ASYNC_BITBANG",
  FT_BITMODE_MPSSE => "FT_BITMODE_MPSSE",
  FT_BITMODE_SYNC_BITBANG => "FT_BITMODE_SYNC_BITBANG",
  FT_BITMODE_MCU_HOST => "FT_BITMODE_MCU_HOST",
  FT_BITMODE_FAST_SERIAL => "FT_BITMODE_FAST_SERIAL",
  FT_BITMODE_CBUS_BITBANG => "FT_BITMODE_CBUS_BITBANG",
  FT_BITMODE_SYNC_FIFO => "FT_BITMODE_SYNC_FIFO")

function ft_getbitmode(ft_handle::Culong)
  mode = Ref{Cuchar}()
  ft_status = ccall((:FT_GetBitMode, d2xx),
                     Cuint,
                     (Culong, Ref{Cuchar}),
                     ft_handle, mode)
  checkstatus(ft_status)
  return mode[]
end


text:
using OrdinaryDiffEq, Plots, DiffEqCallbacks
using TaylorIntegration, LinearAlgebra

T(p) = 1//2 * norm(p)^2
V(q) = 1//2 * (q[1]^2 + q[2]^2 + 2q[1]^2 * q[2]- 2//3 * q[2]^3)
H(p,q, params) = T(p) + V(q)

module InPlace
    using ParameterizedFunctions

    function q̇(dq,p,q,params,t)
        dq[1] = p[1]
        dq[2] = p[2]
    end

    function ṗ(dp,p,q,params,t)
        dp[1] = -q[1] * (1 + 2q[2])
        dp[2] = -q[2] - (q[1]^2 - q[2]^2)
    end

    const q0 = [0.1, 0.]
    const p0 = [0., 0.5]
    const u0 = vcat(p0,q0)

    henon = @ode_def HamiltonEqs begin
        dp1 = -q1 * (1 + 2q2)
        dp2 = -q2 - (q1^2 - q2^2)
        dq1 = p1
        dq2 = p2
    end

end

module OutOfPlace
    using StaticArrays

    function q̇(p, q, params, t)
        p
    end

    function ṗ(p, q, params, t)
        dp1 = -q[1] * (1 + 2q[2])
        dp2 = -q[2] - (q[1]^2 - q[2]^2)
        @SVector [dp1, dp2]
    end

    const q0 = @SVector [0.1, 0.]
    const p0 = @SVector [0., 0.5]
    const u0 = vcat(p0,q0)

    henon(z, p, t) = SVector(
        -z[3] * (1 + 2z[4]),
        -z[4] - (z[3]^2 - z[4]^2),
        z[1],
        z[2]
    )

end

function g(resid, u)
    resid[1] = H([u[1],u[2]], [u[3],u[4]], nothing) - E
    resid[2:4] .= 0
end

const cb = ManifoldProjection(g, nlopts=Dict(:ftol=>1e-13))

const E = H(InPlace.p0, InPlace.q0, nothing)

energy_err(sol) = map(i->H([sol[1,i], sol[2,i]], [sol[3,i], sol[4,i]], nothing)-E, 1:length(sol.u))
abs_energy_err(sol) = [abs.(H([sol[1,j], sol[2,j]], [sol[3,j], sol[4,j]], nothing) - E) for j=1:length(sol.u)]

function compare(mode=InPlace, all=true, plt=nothing; tmax=1e2)
    prob1 = DynamicalODEProblem(mode.ṗ, mode.q̇, mode.p0, mode.q0, (0., tmax))
    prob2 = ODEProblem(mode.henon, mode.u0, (0., tmax))

    GC.gc()
    (mode == InPlace  && all) && @time sol1 = solve(prob2, Vern9(), callback=cb, abstol=1e-14, reltol=1e-14)
    GC.gc()
    @time sol2 = solve(prob1, KahanLi8(), dt=1e-2, maxiters=1e10)
    GC.gc()
    @time sol3 = solve(prob1, SofSpa10(), dt=1e-2, maxiters=1e8)
    GC.gc()
    @time sol4 = solve(prob2, Vern9(), abstol=1e-14, reltol=1e-14)
    GC.gc()
    @time sol5 = solve(prob1, DPRKN12(), abstol=1e-14, reltol=1e-14)
    GC.gc()
    (mode == InPlace && all) && @time sol6 = solve(prob2, TaylorMethod(50), abstol=1e-20)

    (mode == InPlace && all) && println("Vern9 + ManifoldProjection max energy error:\t"*
        "$(maximum(abs_energy_err(sol1)))\tin\t$(length(sol1.u))\tsteps.")
    println("KahanLi8 max energy error:\t\t\t$(maximum(abs_energy_err(sol2)))\tin\t$(length(sol2.u))\tsteps.")
    println("SofSpa10 max energy error:\t\t\t$(maximum(abs_energy_err(sol3)))\tin\t$(length(sol3.u))\tsteps.")
    println("Vern9 max energy error:\t\t\t\t$(maximum(abs_energy_err(sol4)))\tin\t$(length(sol4.u))\tsteps.")
    println("DPRKN12 max energy error:\t\t\t$(maximum(abs_energy_err(sol5)))\tin\t$(length(sol5.u))\tsteps.")
    (mode == InPlace && all) && println("TaylorMethod max energy error:\t\t\t$(maximum(abs_energy_err(sol6)))"*
        "\tin\t$(length(sol6.u))\tsteps.")

    if plt == nothing
        plt = plot(xlabel="t", ylabel="Energy error")
    end
    (mode == InPlace && all) && plot!(sol1.t, energy_err(sol1), label="Vern9 + ManifoldProjection")
    plot!(sol2.t, energy_err(sol2), label="KahanLi8", ls=mode==InPlace ? :solid : :dash)
    plot!(sol3.t, energy_err(sol3), label="SofSpa10", ls=mode==InPlace ? :solid : :dash)
    plot!(sol4.t, energy_err(sol4), label="Vern9", ls=mode==InPlace ? :solid : :dash)
    plot!(sol5.t, energy_err(sol5), label="DPRKN12", ls=mode==InPlace ? :solid : :dash)
    (mode == InPlace && all) && plot!(sol6.t, energy_err(sol6), label="TaylorMethod")

    return plt
end

compare(tmax=1e2)

compare(tmax=1e3)

compare(tmax=1e4)

compare(tmax=5e4)

function in_vs_out(;all=false, tmax=1e2)
    println("In place versions:")
    plt = compare(InPlace, all, tmax=tmax)
    println("\nOut of place versions:")
    plt = compare(OutOfPlace, false, plt; tmax=tmax)
end

in_vs_out(all=true, tmax=1e3)

in_vs_out(tmax=1e2)

in_vs_out(tmax=1e3)

in_vs_out(tmax=1e4)

in_vs_out(tmax=5e4)

using DiffEqBenchmarks
DiffEqBenchmarks.bench_footer(WEAVE_ARGS[:folder],WEAVE_ARGS[:file])


text:
if visualtests
  @testset "Plot recipes" begin
    Random.seed!(2018)
    data₁ = rand(LogNormal(2, 2), 5000)
    data₂ = rand(Pareto(.9, 10), 5000)

    Random.seed!(2018)
    xs = rand(LogNormal(0, 1), 5000)

    @testset "Return levels" begin
      @plottest returnplot(xs,label="log-normal") joinpath(datadir,"ReturnLevels.png") !istravis
    end

    @testset "Mean excess" begin
      @plottest begin
        plt1 = excessplot(data₁, label="Lognormal")
        plt2 = excessplot(data₂, label="Pareto")
        plot(plt1, plt2, size=(600,800), layout=(2,1))
      end joinpath(datadir,"MeanExcess.png") !istravis
    end

    @testset "Pareto quantile" begin
      @plottest begin
        plt1 = paretoplot(data₁, label="Lognormal")
        plt2 = paretoplot(data₂, label="Pareto")
        plot(plt1, plt2, size=(600,800), layout=(2,1))
      end joinpath(datadir,"ParetoQuantile.png") !istravis
    end
  end
end


text:
using EzXML

function Base.iterate(node::EzXML.Node)
    v = [node]
    t = EzXML.Node[]
    while !isempty(v)
        parentnode = pop!(v)
        for childnode in eachelement(parentnode)
            push!(v, childnode)
            push!(t, childnode)
        end
    end
    return t
end


text:
#Functions for building, initiating models and submodels.
function build_gwas()
    return Model(0,0, 0,0, 0, 0, 0, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, 0,0)
end

function initiate_model(model, interval, setsize, outputfilename, printchain, type)
    model.nlayer1 = 1
    model.nlayer2 = 0
    model.losoInterval = Float64(interval)
    model.layerlimit = 0.75 * model.losoInterval
    model.inlayer1 = true
    model.setsize = setsize
    model.notcomplete = true
    model.outputfilename = outputfilename
    model.printchain = printchain
    model.submodelcount = 0
    model.futures = Future[]
    #
    try
        model.m = size(model.snpsOld.binarygenotypes, 2)
        model.maxbp = parse(Float64, model.snpsOld.markerInfo[model.m,2])
    catch
        println("Model already initiated. Please build model again.")
    end
    #
    model.phenoOrder2, model.genoOrder = get_orders(model.phenotypes.IID, model.snpsOld.IID) # keep phenoOrder and phenoOrder2 separate as will have to use both when reading new adjusted phenotypes...
    get_phenotype_intersect!(model.phenotypes, model.phenoOrder2)
    model.weightOrder = model.phenoOrder2 # since pheno and weights have already been sorted/ordered. Weights won't need to be ordered again
    get_weight_intersect!(model.weights, model.weightOrder)
    model.snps = get_genotype_intersect_bin(model.snpsOld, model.genoOrder)
    model.snpsOld = nothing
    model.n = size(model.phenotypes.IID,1)
    model.k = size(model.phenotypes.trait,2)
    model.maxlayerlimit = get_maxlayerlimit(model)
    model.type = type
end

function inititate_submodel!(submodel)
    submodel.rng = MersenneTwister()
    submodel.K = get_1loci_k()
    submodel.output_stats = zeros(Float64,submodel.m, 10)
    if submodel.type == "additive"
        submodel.output_stats = zeros(Float64, submodel.m, 6)
        submodel.K = 0
    end
    submodel.yty = (submodel.phenotype.trait .* submodel.weights.rinverse)'submodel.phenotype.trait #(y .* rinv)'y' where rinv is nx1 vector - faster than y'Rinv*y, where Rinv is a diag matrix. Same allocations
end

# read directory of yloso phenotype files and find largest, calculate max layer limit
function get_maxlayerlimit(model)
    phenotypedir = dirname(model.phenotypesprefix)
    phenotypefiles = readdir(phenotypedir)
    maxn = 1
    maxlayer = 1
    phenotypefiles = phenotypefiles[occursin.(basename(model.phenotypesprefix), phenotypefiles)]
    for file in phenotypefiles
        n,layer = get_n_and_layer(file)
        if n >  maxn
            maxn = n
            maxlayer = layer
        elseif n == maxn
            if layer > maxlayer
                maxlayer = layer
                maxn = n
            end
        end
    end
    maxlayerlimit = maxn * 10000000
    if maxlayer == 2
        maxlayerlimit = maxlayerlimit + 5000000
    end
    return maxlayerlimit
end

function get_n_and_layer(x)
    sub = x[1:findlast(isequal('.'),x)-1]
    n = parse(Int64, split(sub[(findlast(isequal('.'),sub)+1):end], "n")[2])
    sub2 = sub[1:findlast(isequal('.'),sub)-1]
    layer = parse(Int64, split(sub2[(findlast(isequal('.'),sub2)+1):end], "layer")[2])
    return n, layer
end

max_stars_repo_name: egmreynolds/MANA

max_issues_repo_name: egmreynolds/MANA

max_forks_repo_name: egmreynolds/MANA


text:
function assert_gamma(γ)
    if !(-1 <= γ <= 0.0)
        throw(ArgumentError("unsupported γ: must be -1 <= γ <= 0)"))
    end
end

function assert_pade(pade)
    if pade ∉ 4:5
        throw(ArgumentError("4 or 5 pade approximation is supported"))
    end
end

function assert_fftlen(fftlen)
    if !ispow2(fftlen)
        throw(ArgumentError("fftlen must be power of 2"))
    end
end


text:
function conic_linear1_test()
    #=
        min -3x - 2y - 4z
    s.t.
        x +  y +  z == 3
        y +  z == 2
        x>=0 y>=0 z>=0
    =#
    model = TestModel{Float64}()

    v = MOI.add_variables(model, 3)

    vov = MOI.VectorOfVariables(v)

    vc = MOI.add_constraint(model, vov, MOI.Nonnegatives(3))

    c = MOI.add_constraint(
        model,
        MOI.VectorAffineFunction(
            MOI.VectorAffineTerm.(
                [1, 1, 1, 2, 2],
                MOI.ScalarAffineTerm.(1.0, [v; v[2]; v[3]]),
            ),
            [-3.0, -2.0],
        ),
        MOI.Zeros(2),
    )

    MOI.set(
        model,
        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
        MOI.ScalarAffineFunction(
            MOI.ScalarAffineTerm.([-3.0, -2.0, -4.0], v),
            0.0,
        ),
    )
    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)

    return model
end

function conic_linear2_test()
    #=
        min -3x - 2y - 4z
    s.t.
        x +  y +  z == 3
        y +  z == 2
        x>=0 y>=0 z>=0
    =#
    model = TestModel{Float64}()

    v = MOI.add_variables(model, 3)

    vov = MOI.VectorOfVariables(v)

    vc = MOI.add_constraint(
        model,
        MOI.VectorAffineFunction{Float64}(vov),
        MOI.Nonnegatives(3),
    )

    c = MOI.add_constraint(
        model,
        MOI.VectorAffineFunction(
            MOI.VectorAffineTerm.(
                [1, 1, 1, 2, 2],
                MOI.ScalarAffineTerm.(1.0, [v; v[2]; v[3]]),
            ),
            [-3.0, -2.0],
        ),
        MOI.Zeros(2),
    )

    MOI.set(
        model,
        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
        MOI.ScalarAffineFunction(
            MOI.ScalarAffineTerm.([-3.0, -2.0, -4.0], v),
            0.0,
        ),
    )
    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)

    return model
end

function conic_linear3_test()
    #=
        min  3x + 2y - 4z + 0s
    s.t.
           x           -  s  == -4    (i.e. x >= -4)
                y            == -3
           x      +  z       == 12
           x free
           y <= 0
           z >= 0
           s zero
    =#
    model = TestModel{Float64}()

    x, y, z, s = MOI.add_variables(model, 4)

    MOI.set(
        model,
        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
        MOI.ScalarAffineFunction(
            MOI.ScalarAffineTerm.([3.0, 2.0, -4.0], [x, y, z]),
            0.0,
        ),
    )
    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)

    c = MOI.add_constraint(
        model,
        MOI.VectorAffineFunction(
            MOI.VectorAffineTerm.(
                [1, 1, 2, 3, 3],
                MOI.ScalarAffineTerm.(
                    [1.0, -1.0, 1.0, 1.0, 1.0],
                    [x, s, y, x, z],
                ),
            ),
            [4.0, 3.0, -12.0],
        ),
        MOI.Zeros(3),
    )

    vov = MOI.VectorOfVariables([y])

    vc = MOI.add_constraint(model, vov, MOI.Nonpositives(1))

    # test fallback
    vz = MOI.add_constraint(model, [z], MOI.Nonnegatives(1))

    vov = MOI.VectorOfVariables([s])

    vs = MOI.add_constraint(model, vov, MOI.Zeros(1))

    return model
end

function conic_linear4_test()
    #=
        min  3x + 2y - 4z + 0s
    s.t.
           x           -  s  == -4    (i.e. x >= -4)
                y            == -3
           x      +  z       == 12
           x free
           y <= 0
           z >= 0
           s zero
    =#
    model = TestModel{Float64}()

    x, y, z, s = MOI.add_variables(model, 4)

    MOI.set(
        model,
        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
        MOI.ScalarAffineFunction(
            MOI.ScalarAffineTerm.([3.0, 2.0, -4.0], [x, y, z]),
            0.0,
        ),
    )
    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)

    c = MOI.add_constraint(
        model,
        MOI.VectorAffineFunction(
            MOI.VectorAffineTerm.(
                [1, 1, 2, 3, 3],
                MOI.ScalarAffineTerm.(
                    [1.0, -1.0, 1.0, 1.0, 1.0],
                    [x, s, y, x, z],
                ),
            ),
            [4.0, 3.0, -12.0],
        ),
        MOI.Zeros(3),
    )

    vov = MOI.VectorOfVariables([y])
    vc = MOI.add_constraint(
        model,
        MOI.VectorAffineFunction{Float64}(vov),
        MOI.Nonpositives(1),
    )

    vz = MOI.add_constraint(
        model,
        MOI.VectorAffineFunction(
            [MOI.VectorAffineTerm(1, MOI.ScalarAffineTerm(1.0, z))],
            [0.0],
        ),
        MOI.Nonnegatives(1),
    )

    vov = MOI.VectorOfVariables([s])

    vs = MOI.add_constraint(
        model,
        MOI.VectorAffineFunction{Float64}(vov),
        MOI.Zeros(1),
    )

    return model
end


text:
# Dual infeasible test problems
# Choose A, P and q in a way such that the problem becomes unbounded below
# here the last element of x appears in the cost function with a negative sign and is unconstrained

# This test is precision agnostic and will be run with T precision
if !@isdefined(UnitTestFloats)
    UnitTestFloats = [Float64] #if not run in full test setup, just do it for one float type
end

for T in UnitTestFloats
  local rng
  if precision(T) >= precision(Float64)

    rng = Random.MersenneTwister(555)
    nn = 1

    @testset "Dual infeasible QP problems - Testset 1" begin
      for iii = 1:nn

      # choose size of problem
      n = rand(rng,5:50)
      m = 2*n
      A = sprand(rng, T, m,n,0.7)*T(50)
      A[:,end] .= zero(T)

      P = spzeros(T, n,n)

      q = rand(rng, T, n, 1)*T(50)
      q[end] = -one(T)
      q = vec(q)
      xtrue = rand(rng, T, n, 1)*T(50)
      strue = rand(rng, T, m, 1)*T(50)
      b = A*xtrue+strue
      b = vec(b)

      constraint = COSMO.Constraint(-A,b,COSMO.Nonnegatives)

      settings = COSMO.Settings{T}(max_iter=10000,eps_abs = 1e-5,eps_rel=1e-5)
      model = COSMO.Model{T}()
      assemble!(model,P,q,[constraint], settings = settings)
      res = COSMO.optimize!(model);

      @test res.status == :Dual_infeasible
      end
    end
  end
end


text:
@defcomp SulphateForcing begin
    region = Index()

    se0_sulphateemissionsbase = Parameter(index=[region], unit="TgS/year")
    pse_sulphatevsbase = Parameter(index=[time, region], unit="%")
    se_sulphateemissions = Variable(index=[time, region], unit="TgS/year")
    area = Parameter(index=[region], unit="km2")
    area_e_eartharea = Parameter(unit="km2", default=5.1e8)

    sfx_sulphateflux = Variable(index=[time, region], unit="TgS/km2/yr")

    d_sulphateforcingbase = Parameter(unit="W/m2", default=-0.46666666666666673)
    ind_slopeSEforcing_indirect = Parameter(unit="W/m2", default=-0.2333333333333333)
    nf_naturalsfx = Parameter(index=[region], unit="TgS/km2/yr")

    fsd_directsf = Variable(index=[time, region], unit="W/m2")
    fsi_indirectsf = Variable(index=[time, region], unit="W/m2")
    fs_sulphateforcing = Variable(index=[time, region], unit="W/m2")

    function run_timestep(p, v, d, tt)
        bigSFX0 = p.se0_sulphateemissionsbase ./ p.area

        for rr in d.region
            v.se_sulphateemissions[tt, rr] = p.se0_sulphateemissionsbase[rr] * p.pse_sulphatevsbase[tt, rr] / 100

            # Eq.17 from Hope (2006) - sulfate flux
            v.sfx_sulphateflux[tt,rr] = v.se_sulphateemissions[tt,rr] / p.area[rr]
            # Update for Eq. 18 from Hope (2009) - sulfate radiative forcing effect
            bigSFD0 = p.d_sulphateforcingbase * bigSFX0[rr] / (sum(bigSFX0 .* p.area) / p.area_e_eartharea)
            v.fsd_directsf[tt, rr] = bigSFD0 * v.sfx_sulphateflux[tt,rr] / bigSFX0[rr]
            v.fsi_indirectsf[tt, rr] = p.ind_slopeSEforcing_indirect / log(2) * log((p.nf_naturalsfx[rr] + v.sfx_sulphateflux[tt, rr]) / p.nf_naturalsfx[rr])

            v.fs_sulphateforcing[tt, rr] = v.fsd_directsf[tt, rr] + v.fsi_indirectsf[tt, rr]
        end
    end
end


text:
# hydra-growing example
#

using ShapeFactory

compass = CartesianIndex.(((1, 0), (0, -1), (-1, 0), (0, 1)))

function scale_idx(idx, n)
    return idx + (idx - one(idx)) * (n-1)
end

function stencil(grid, index, hv)
    for di in CartesianIndices((0:1, 0:1))
        grid[index + di] = (grid[index + di] + hv) % 4
    end
end

function hydra_grid(height, width, maxgen=4)
    sub_grid = zeros(Int8, height * 2, width * 2)
    top_grid = Set(CartesianIndices((height, width)))
    while length(top_grid) > 0
        group = Dict{CartesianIndex, BitArray}()
        group[pop!(top_grid)] = [0, 0, 0, 0]

        gen = 0
        while length(group) > 0
            indexa, edgesa = rand(group)
            side = first(rand(filter(x -> last(x) == 0, collect(enumerate(edgesa)))))
            indexb = indexa + compass[side]
            edgesa[side] = 1
            if all(edgesa)
                pop!(group, indexa)
            end
            if indexb in top_grid
                di = compass[side]
                stencil(sub_grid, scale_idx(indexa, 2) + di, 2 + (side-1)%2)
                group[indexb] = [0, 0, 0, 0]
                delete!(top_grid, indexb)
                if (gen += 1) > maxgen
                    break
                end
            end
        end
    end
    return sub_grid
end

grid = hydra_grid(10, 10)

app = dispatch("CATIA.Application")
prt = app.ActiveDocument.Part

gst = prt.HybridBodies.Add()
gst.Name = "Output"

fac = prt.HybridShapeFactory
xyp = prt.FindObjectByName("xy plane")

pts = fac.AddNewPointCoord.((1000, 0, 1000, 0), (1000, 1000, 0, 0), 0)
ccs = fac.AddNewCircleCtrRadWithAngles.(pts, xyp, false, 500,
                                        (180, 270, 90, 0), (270, 360, 180, 90))

pts = fac.AddNewPointCoord.((500, 0, 500, 1000), (0, 500, 1000, 500), 0)
lns = fac.AddNewLinePtPt.(pts[1:2], pts[3:4])

gst.AppendHybridShape.(ccs)
gst.AppendHybridShape.(lns)

for y in 0:size(grid, 1)-1
    for x in 0:size(grid, 2)-1
        mx, my = x * 1000, y * 1000
        xv = fac.AddNewDirectionByCoord(mx, my, 1)
        xd = sqrt(mx^2 + my^2)
        st = grid[y + 1, x + 1]
        if st > 1
            tx = fac.AddNewTranslate(lns[st-1], xv, xd)
        else
            tx = fac.AddNewTranslate(ccs[((x+st)%2 + ((y+st)%2)*2) + 1], xv, xd)
        end
        tx.VectorType = 0
        gst.AppendHybridShape(tx)
    end
end

prt.Update()

max_stars_repo_name: zaika-generative/ShapeFactory

max_issues_repo_name: zaika-generative/ShapeFactory

max_forks_repo_name: zaika-generative/ShapeFactory


text:
mutable struct Way

    length::Float64
    pheromone::Float64

end


text:
module IPynbToJl

using JSON

export ipynb_to_jl

function ipynb_to_jl(infile, outfile=replace(infile, r"\.ipynb", ".jl"))
    f = JSON.parsefile(infile, use_mmap=false)
    f == nothing && error("Could not open [$infile]")
    outf = open(outfile, "w")
    write(outf, "#  ", abspath(outfile), "\n")
    write(outf, "#  AUTOGENERATED FROM ", abspath(infile), " on ", string(Dates.now()), "\n")

    for c in f["worksheets"][1]["cells"]
        cell_type = c["cell_type"]
        if cell_type == "code"
            write(outf, c["input"])
            write(outf, "\n")
        else
            write(outf, "#[$cell_type]\n")
            write(outf, map(l -> "# $l\n", c["source"]))
            write(outf, "\n")
        end
    end
    close(outf)
    print("Wrote $outfile\n") # If generating a source file, I like to ensure there's some output to stdout
end

end

max_stars_repo_name: dsyme/ADBench

max_issues_repo_name: dsyme/ADBench

max_forks_repo_name: dsyme/ADBench


text:
module GLParse
using GLUT

export createOpenGLFun

macro createOpenGLFun(cFun)
    cName = cFun.args[1].args[1]
    arguments = map(function (arg)
                        if isa(arg, Symbol)
                            arg = Expr(__ofTypeSym, arg)
                        end
                        return arg
                    end, cFun.args[1].args[2:end])

    # Get info out of arguments of `cFun`
    argumentNames = map(arg->arg.args[1], arguments)
    returnType    = cFun.args[2]
    inputTypes    = map(arg->arg.args[2], arguments)
    # Construct the result.
    functionPtr = glutGetProcAddress(cName)
    if functionPtr == C_NULL
        println("Warning: returned nullPointer -- Function:$(string(cFun)) not supported by OpenGL Driver")
        return ""
    end
    body = Expr(:ccall, functionPtr, returnType, Expr(:tuple,inputTypes...), argumentNames...)
    ret  = Expr(:function, Expr(:call, cName, argumentNames...), body)
    return esc(ret)
end

end #Module GLParse


text:
# ------------------------------------------------------------------
# Licensed under the ISC License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    iscompatible(model, task)

Check whether or not `model` can be used for
learning `task`.
"""
iscompatible(model::MLJBase.Model, task::AbstractLearningTask) = false
iscompatible(model::MLJBase.Model, task::RegressionTask) = issupervised(model) &&
  (MLJBase.target_scitype(model) == AbstractVector{MLJBase.Continuous})
iscompatible(model::MLJBase.Model, task::ClassificationTask) = issupervised(model) &&
  (MLJBase.target_scitype(model) == AbstractVector{<:MLJBase.Finite})
iscompatible(model::MLJBase.Model, task::ClusteringTask) = !issupervised(model)

"""
    isprobabilistic(model)

Check whether or not `model` is probabilistic.
"""
isprobabilistic(model::MLJBase.Model) = false
isprobabilistic(model::MLJBase.Probabilistic) = true

"""
    issupervised(model)

Check whether or not `model` is supervised.
"""
issupervised(model::MLJBase.Model) = false
issupervised(model::MLJBase.Supervised) = true

"""
    defaultloss(val)
    defaultloss(scitype)

Default loss for value `val` or its scientific type `scitype`.
"""
defaultloss(val) = defaultloss(MLJBase.scitype(val))
defaultloss(::Type{<:MLJBase.Finite}) =
  MLJBase.misclassification_rate
defaultloss(::Type{<:MLJBase.Infinite}) =
  MLJBase.rms


text:
using GPlot, Colors, Random; PREVIEW = false; SAVEFIG = true; SAVEPATH = "../GPlotExamples.jl/examples/";

continuous_preview(false) # considerable speedup...

####
#### Simple line plot, no latex
#### NOTE: this is super quick (after precomp)
#### because there's no latex compilation pass
####

begin
# dummy data
    x1 = range(-2, stop=2, length=100)
    y1 = @. exp(-x1 * sin(x1))
    y2 = @. exp(-x1 * cos(x1))

    t = @elapsed begin
        f = Figure("simple_notex", reset=true)
        # plot things with different colors, markers etc
        plot(x1, y1, color="darkblue", lwidth=0.02)
        plot(x1, y2, color="darkred", lwidth=0.02) # will overwrite previous
        plot!(x1, y1, color="darkblue", lwidth=0.02)
        # -----
        PREVIEW && preview(f)
        SAVEFIG && savefig(f, format="pdf", path=SAVEPATH)
    end; println("$(f.id)...done in $(round(t, digits=2))s")
end

####
#### Simple line plot, with titles and latex
####

begin
    # dummy data
    x1 = range(-2, stop=2, length=100)
    y1 = @. exp(-x1 * sin(x1))
    y2 = @. exp(-x1 * cos(x1))
    x2 = range(0, stop=2, length=10)
    y3 = @. sqrt(x2)
    y4 = @. exp(sin(x1)-x1^2)*x1
    x3 = range(-1, stop=1, length=15)
    y5 = @. x3^2
    y6 = @. -y5+1

    t = @elapsed begin
        texpreamble = tex"""
            \usepackage[T1]{fontenc}
            \usepackage[default]{sourcesanspro}
        """
        f = Figure("simple_tex", texpreamble=texpreamble, reset=true)
        # plot things with different colors, markers etc
        plot!(x1, y1, color="darkblue", lwidth=0.02, lstyle=3)
        plot!(x1, y2, color="indianred")
        plot!(x2, y3, lstyle="none", marker="fcircle", msize=0.1, color="#0c88c2")
        plot!(x1, y4, color="#76116d", lwidth=0.1)
        plot!(x3, y5, ls="-", color="orange", lwidth=0.05, marker="o", label="yolo")
        plot!(x3, y6, ls="-", color="darkcyan", lwidth=0.05, marker=".")
        # tex strings for the titles
        xtitle!(tex"The $x$ axis $\int_0^\infty f(x)\mathrm{d}x$", color="blue",
            fontsize=12)
        x2title(tex"$x_2$")
        y2title(tex"axis $y_2$")
        ytitle(tex"Axis $y$", fontsize=12)
        title(tex"The new title $\mathcal N$")
        legend()
        # -----
        PREVIEW && preview(f)
        SAVEFIG && savefig(f; format="pdf", path=SAVEPATH)
    end; println("$(f.id)...done in $(round(t, digits=2))s")
end

####
#### Simple multiple line plots (in for loop) for square root with latex
#### [GLE EXAMPLE]
####

t = @elapsed begin
    f = Figure("sqroot_tex", size=(11, 8), latex=true, reset=true)
    x = range(0, stop=10, length=100);
    alphas = 1:10
    for alpha ∈ alphas
        plot!(x, sqrt.(alpha * x), col=RGB(alpha/10,0,0),
              label=tex"$\alpha=##alpha$")
    end
    xtitle!(t"$x$")
    ytitle!(t"$\sqrt{\alpha x}$")
    title!("Square Root Function")
    legend(pos="tl", fontsize=7)
    # -----
    PREVIEW && preview(gcf())
    SAVEFIG && savefig(f, format="pdf", path=SAVEPATH)
end; println("$(f.id)...done in $(round(t, digits=2))s")

####
#### Simple line plots with centered axis for sine with latex
#### [GLE EXAMPLE]
####

begin
    # dummy data
    x = range(-2pi, 2pi, length=100)
    y = sin.(x)

    t = @elapsed begin
        f = Figure("sine_tex", latex=true, fontsize=8, reset=true)
        plot!(x, y, col="red", smooth=true)
        title!(tex"$f(x)=\sin(x)$", dist=0.3)
        ax = gca()
        ax.math = true
        ax.xaxis.min = -2pi
        ax.xaxis.max = 2pi
        ax.xaxis.ticks.places = [-4, -3, -2, -1, 1, 2, 3, 4]/2*pi
        ax.xaxis.ticks.labels = GPlot.TicksLabels(names=[t"$-2\pi$", t"$-3\pi/2$", t"$-\pi$", t"$-\pi/2$", t"$\pi/2$", t"$\pi$", t"$3\pi/2$", t"$2\pi$"])
        ax.yaxis.ticks.places = [-4, -3, -2, -1, 1, 2, 3, 4]/4
        ax.yaxis.ticks.labels = GPlot.TicksLabels(names=["-1", "-3/4", "-1/2", "-1/4", "1/4", "1/2", "3/4", "1"])
        # -----
        PREVIEW && preview(gcf())
        SAVEFIG && savefig(f, format="pdf", path=SAVEPATH)
    end; println("$(f.id)...done in $(round(t, digits=2))s")
end

####
#### Simple histogram with pdf normalisation and pdf fit added
####

begin
    # dummy data
    x = randn(10_000)
    xx = range(-4, 4, length=100)
    y  = @. exp(-xx^2/2)/sqrt(2pi)

    t = @elapsed begin
        f = Figure("simple_hist_notex", reset=true)
        hist(x, fill="CornflowerBlue", color="white", scaling="pdf", nbins=50)
        plot!(xx, y)
        # -----
        PREVIEW && preview(f)
        SAVEFIG && savefig(f, format="pdf", path=SAVEPATH)
    end; println("$(f.id)...done in $(round(t, digits=2))s")
end

####
#### Fill between lines with transparency
#### NOTE: alpha=true must be there otherwise GLE errors!
####

begin
    # dummy data
    x = 0.1:0.1:5
    y1 = @. x^2 / exp(x)
    y2 = @. x^3 / exp(x)
    y3 = @. exp(-x^2)

    t = @elapsed begin
        f = Figure("simple_fill_transp_notex", alpha=true, reset=true)
        GPlot.add_axes2d!()
        fill_between!(x, y1, y2, alpha=0.5)
        fill_between!(x, 0, y3, color="red", alpha=0.5)
        # -----
        PREVIEW && preview(f)
        SAVEFIG && savefig(f, format="pdf", path=SAVEPATH)
    end; println("$(f.id)...done in $(round(t, digits=2))s")
end

####
#### Simple log plot, no latex
####

begin
    # dummy data
    x = range(1, stop=1000, length=100)
    y = @. log(x)

    t = @elapsed begin
        f = Figure("simple_logscale_notex", reset=true)
        plot(x, y, lstyle="--", lw=0.1)
        xscale!("log")
        # -----
        PREVIEW && preview(f)
        SAVEFIG && savefig(f, format="pdf", path=SAVEPATH)
    end; println("$(f.id)...done in $(round(t, digits=2))s")
end

####
#### Simple bar, no latex
####

begin
    # dummy data
    y = [50, 10, 20, 30, 5, 100]

    t = @elapsed begin
        f = Figure("simple_bar_notex", reset=true)
        bar(y, fill="indianred", color="orange", horiz=true)
        ylim!(0.5, 6.5)
        # -----
        PREVIEW && preview(f)
        SAVEFIG && savefig(f, format="pdf", path=SAVEPATH)
    end; println("$(f.id)...done in $(round(t, digits=2))s")
end

####
#### Simple grouped bar plot, no tex
####

begin
    # dummy data
    y = [10, 20, 15, 50]
    δ = [5, 3, 7, 10]
    y2 = y .+ δ
    x = 1:length(y)

    for stacked ∈ [true, false]
        t = @elapsed begin
            f = Figure("simple_gbar_stacked=$(stacked)_notex", reset=true)
            bar(x, y, y2; stacked=stacked,
                fills=["darkgreen", "cyan"],
                colors=["white", "white"])
            # -----
            PREVIEW && preview(f)
            SAVEFIG && savefig(f, format="pdf", path=SAVEPATH)
        end; println("$(f.id)...done in $(round(t, digits=2))s")
    end
end

####
#### Axis modifications: ticks, lims
####

begin
    x = range(0, stop=2, length=10)
    y = @. x^2 / exp(x)

    t = @elapsed begin
        f = Figure("simple_lim_ticks_notex", reset=true)

        plot(x, y, lw=0.02, smooth=false)

        xlim(min=-0.5, max=2.5)
        ylim(min=-0.5, max=1.0)

        xticks([0.0, 1.5], ["A", "BB"])
        yticks([-0.2, 0.2, 0.8], ["C", "dd", "ee"], color="blue")

        x2ticks([0.5, 1.0], tickscol="indianred")

        PREVIEW && preview(f)
        SAVEFIG && savefig(f, format="pdf", path=SAVEPATH)
    end; println("$(f.id)...done in $(round(t, digits=2))s")
end

####
#### Simple line plot with scatter and axis manips
####

begin
    x = range(-5, 5, length=100)
    xs = range(-5, 5, length=20)

    t = @elapsed begin
        f = Figure("simple_scatter_notex", reset=true)

        λ(x) = 1-exp(-sin(x) * cos(x)/x)

        plot(x, λ.(x), col="blue")
        scatter!(xs, λ.(xs), col="red")

        ylim(-.5, .75)

        y2axis("off")
        x2axis("off")

        grid()
        xticks(-5:2.5:5, tickscol="lightgray", angle=45)
        yticks(-.5:0.25:.75, tickscol="lightgray")

        PREVIEW && preview(f)
        SAVEFIG && savefig(f, format="pdf", path=SAVEPATH)
    end; println("$(f.id)...done in $(round(t, digits=2))s")
end

####
#### Missing values
####

begin
    x  = range(-1, 1, length=18)
    y  = @. exp(-x^2*cos(10x))
    Random.seed!(0)
    mask = [rand()<0.25 for i ∈ eachindex(y)]
    ym = convert(Vector{GPlot.CanMiss{Float64}}, y)
    ym[mask] .= missing

    t = @elapsed begin
        f = Figure("simple_withmissing_notex", reset=true)
        plot(x, ym)
        PREVIEW && preview(f)
        SAVEFIG && savefig(f, format="pdf", path=SAVEPATH)
    end; println("$(f.id)...done in $(round(t, digits=2))s")
end

####
#### Butterfly plot
####

begin
    t = range(0, 12pi, length=5000);
    r = @. exp(cos(t)) - 2*cos(4*t) + sin(t/12)^5;
    x = @. r * cos(t);
    y = @. r * sin(t);

    t = @elapsed begin
        f = Figure("butterfly_notex", reset=true)
        plot(x, y)
        math!()
        PREVIEW && preview(f)
        SAVEFIG && savefig(f, format="pdf", path=SAVEPATH)
    end; println("$(f.id)...done in $(round(t, digits=2))s")
end

####
#### Legend styling (and implicit)
####

begin
    t = @elapsed begin
        f = Figure("legend_styling_notex", reset=true)

        plot(tan, 0, pi/2)
        scatter!(sin, 0, pi/2)
        plot!(cos, 0, pi/2)
        ylim(0, 2)

        legend(font="tt", fontsize=12, offset=(1,1), nobox=false, bgcol="blue")
    end; println("$(f.id)...done in $(round(t, digits=2))s")
end


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

@inline nop_func(x...) = nothing

struct CombinedCallback{N,Fs<:NTuple{N,Function}} <: Function
    fs::Fs
end

function Base.show(io::IO, f::CombinedCallback)
    print(io, Base.typename(typeof(f)).name, "(")
    show(io, f.fs)
    print(io, ")")
end

function Base.show(io::IO, M::MIME"text/plain", f::CombinedCallback)
    print(io, Base.typename(typeof(f)).name, "(")
    show(io, M, f.fs)
    print(io, ")")
end

@inline _run_callbacks(fs::NTuple{0,Function}, args...; kwargs...) = nothing

@inline function _run_callbacks(fs::NTuple{1,Function}, args...; kwargs...)
    first(fs)(args...; kwargs...)
    nothing
end

@inline function _run_callbacks(fs::NTuple{N,Function}, args...; kwargs...) where N
    first(fs)(args...; kwargs...)
    _run_callbacks(Base.tail(fs), args...; kwargs...)
    nothing
end

@inline (f::CombinedCallback)(args...; kwargs...) = _run_callbacks(f.fs, args...; kwargs...)

combine_callbacks(::typeof(nop_func), f::Function) = f
combine_callbacks(f::Function, ::typeof(nop_func)) = f
combine_callbacks(::typeof(nop_func), ::typeof(nop_func)) = nop_func

combine_callbacks(f::Function, g::Function) = CombinedCallback((f, g))
combine_callbacks(f::CombinedCallback, g::Function) = CombinedCallback((f.fs..., g))
combine_callbacks(f::Function, g::CombinedCallback) = CombinedCallback((f, g.fs...))
combine_callbacks(f::CombinedCallback, g::CombinedCallback) = CombinedCallback((f.fs..., g.fs...))

near_neg_inf(::Type{T}) where T<:Real = T(-1E38) # Still fits into Float32


text:
function model_called(obj::ModelObserver, arg0::Vector{Orbit}, arg1::Map)
    return jcall(obj, "modelCalled", void, (Vector{Orbit}, Map), arg0, arg1)
end


text:
using SignatureGB
using Singular
P, (x1, x2, x3, x4, x5, x6, x7, x8)  = Singular.PolynomialRing(Fp(65521), ["x1", "x2", "x3", "x4", "x5", "x6", "x7", "x8"]);
gb = SignatureGB.f5([x1, x2, x3, x4, x5, x6, x7, x8]);F = [x1^2*x4^2-46*x1*x2^2*x4+166*x1*x3*x4*x5+52*x1*x4^2*x7+529*x2^4-3818*x2^2*x3*x5-1196*x2^2*x4*x7+961*x2^2*x5^2+434*x2^2*x5*x6+49*x2^2*x6^2-2604*x2*x3^2*x5-588*x2*x3^2*x6-4526*x2*x4^2*x5-1022*x2*x4^2*x6+4278*x2*x4*x5*x6+1488*x2*x4*x5*x8+966*x2*x4*x6^2+336*x2*x4*x6*x8+1764*x3^4+6132*x3^2*x4^2-5796*x3^2*x4*x6-2016*x3^2*x4*x8+6889*x3^2*x5^2+4316*x3*x4*x5*x7+5329*x4^4-10074*x4^3*x6-3504*x4^3*x8+4761*x4^2*x6^2+3312*x4^2*x6*x8+676*x4^2*x7^2+576*x4^2*x8^2-28*x1*x2*x4+110*x1*x4*x7+644*x2^3-2530*x2^2*x7-2324*x2*x3*x5-728*x2*x4*x7+9130*x3*x5*x7+2860*x4*x7^2+196*x2^2-1540*x2*x7+3025*x7^2, -92*x1*x2*x4+2116*x2^3-7636*x2*x3*x5-2392*x2*x4*x7+1922*x2*x5^2+868*x2*x5*x6+98*x2*x6^2-2604*x3^2*x5-588*x3^2*x6-4526*x4^2*x5-1022*x4^2*x6+4278*x4*x5*x6+1488*x4*x5*x8+966*x4*x6^2+336*x4*x6*x8-28*x1*x4+1932*x2^2-5060*x2*x7-2324*x3*x5-728*x4*x7+392*x2-1540*x7, 166*x1*x4*x5-3818*x2^2*x5-5208*x2*x3*x5-1176*x2*x3*x6+7056*x3^3+12264*x3*x4^2-11592*x3*x4*x6-4032*x3*x4*x8+13778*x3*x5^2+4316*x4*x5*x7-2324*x2*x5+9130*x5*x7, 2*x1^2*x4-46*x1*x2^2+166*x1*x3*x5+104*x1*x4*x7-1196*x2^2*x7-9052*x2*x4*x5-2044*x2*x4*x6+4278*x2*x5*x6+1488*x2*x5*x8+966*x2*x6^2+336*x2*x6*x8+12264*x3^2*x4-5796*x3^2*x6-2016*x3^2*x8+4316*x3*x5*x7+21316*x4^3-30222*x4^2*x6-10512*x4^2*x8+9522*x4*x6^2+6624*x4*x6*x8+1352*x4*x7^2+1152*x4*x8^2-28*x1*x2+110*x1*x7-728*x2*x7+2860*x7^2, 166*x1*x3*x4-3818*x2^2*x3+1922*x2^2*x5+434*x2^2*x6-2604*x2*x3^2-4526*x2*x4^2+4278*x2*x4*x6+1488*x2*x4*x8+13778*x3^2*x5+4316*x3*x4*x7-2324*x2*x3+9130*x3*x7, 434*x2^2*x5+98*x2^2*x6-588*x2*x3^2-1022*x2*x4^2+4278*x2*x4*x5+1932*x2*x4*x6+336*x2*x4*x8-5796*x3^2*x4-10074*x4^3+9522*x4^2*x6+3312*x4^2*x8, 52*x1*x4^2-1196*x2^2*x4+4316*x3*x4*x5+1352*x4^2*x7+110*x1*x4-2530*x2^2-728*x2*x4+9130*x3*x5+5720*x4*x7-1540*x2+6050*x7, 1488*x2*x4*x5+336*x2*x4*x6-2016*x3^2*x4-3504*x4^3+3312*x4^2*x6+1152*x4^2*x8];
sort!(F, by = p -> leading_monomial(p))
println(@elapsed SignatureGB.f5(F, verbose=true));exit()


text:
using Documenter, AsapMapper

makedocs(
    modules = [AsapMapper],
    sitename = "AsapMapper",
    format = Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true"),
    authors = "Mark Hildebrand, Arthur Hlaing",
    pages = [
        "Home" => "index.md",
        "tutorial.md"
    ]
)

deploydocs(
    repo = "github.com/hildebrandmw/AsapMapper.jl.git",
    target = "build",
)


text:
export Gaussian, GaussianMixture, sample

"""
Defines Gaussian PDF

p(x) = α/(sqrt.((2*π)^d * prod(σ) )) *  exp.(-sum(((x .- μ)./sqrt.(σ)).^2,dims=1)./2)

"""
struct Gaussian{R<:Real, A <: AbstractVector{R}}
    d::Int
    σ::A
    μ::A
    α::R
end
Gaussian(R,d::Int) = Gaussian(d,ones(R,d),zeros(R,d),one(R))
Gaussian(d::Int) = Gaussian(d,ones(d),zeros(d),1.0)
Gaussian(d::Int,σ::AbstractVector{R},μ::AbstractVector{R},α=one(R)) where R<: Real = Gaussian(d,σ,μ,α)

mean(G::Gaussian) = G.μ
std(G::Gaussian) = G.μ

function (G::Gaussian)(X::AbstractArray{R,2}) where R <: Real
    t1 = G.α./( (R(2*pi))^(G.d/2) * sqrt(prod(G.σ)))
    t2 = exp.(-sum(((X .- G.μ)./sqrt.(G.σ)).^2,dims=1)./R(2))
    return vec(t1 .* t2)
end

function sample(G::Gaussian{R,Vector{R}}, n) where R <: Real
    X = randn(R,G.d,n)
    X = sqrt.(G.σ) .* X .+ G.μ
    return X
end

struct GaussianMixture{R<:Real, A <: AbstractVector{R}}
    Gs::Array{Gaussian{R,A}}
end

function (GM::GaussianMixture{R})(X::AbstractArray{R,2}) where R <: Real
    p = GM.Gs[1](X)
    for k=2:length(GM.Gs)
        p += GM.Gs[k](X)
    end
    return p
end

function sample(GM::GaussianMixture{R,Vector{R}},n::Int) where R <: Real
    # determine proportion of samples from each Gaussian
    totalMass = zero(R)
    for k=1:length(GM.Gs)
        totalMass += GM.Gs[k].α
    end
    X = zeros(R,GM.Gs[1].d,0)
    for k=1:length(GM.Gs)-1
        X = [X sample(GM.Gs[k],Int(round(n*GM.Gs[k].α/totalMass)))]
    end
    X = [X sample(GM.Gs[end],n-size(X,2))]
    return X
end

max_stars_repo_name: EmoryMLIP/MFGnet

max_issues_repo_name: EmoryMLIP/MFGnet

max_forks_repo_name: EmoryMLIP/MFGnet


text:
# Point conversion routines
using Test
import CryptoGroups.Specs: bits2octet, octet2bits, int2octet, octet2int, ECP, octet, point, EC2N, PB
import CryptoGroups: @bin_str, bitlength, @hex_str

let
    x = UInt8[0, 129]
    @test bits2octet(octet2bits(x)) == x

    a = octet2bits(x)
    # First leftmostbits are padded to zero to make a full octet
    @test bits2octet(a[4:end]) == bits2octet(a)

    @test octet2bits(x, 13) == octet2bits(x)[4:end]
end

@test int2octet(123456789, 4) == hex"075BCD15"

@test octet2int(hex"0003ABF1CD") == 61600205

@test int2octet(94311, bitlength(104729)) == hex"017067"

@test bits2octet(bin"11011011011101111001101111110110111110001") == hex"01B6EF37EDF1"

@test octet2int(hex"01E74E") == 124750

@test octet2bits(hex"0117B2939ACC", 41) == bin"10001011110110010100100111001101011001100"

# Field element to integer conversion

@test octet2int(bits2octet(bin"11111111001000010011110000110011110101110")) == 2191548508078

let

    prime_curve = ECP(
        p = 6277101735386680763835789423207666416083908700390324961279,
        a = 6277101735386680763835789423207666416083908700390324961276,
        b = 2455155546008943817740293915197451784769108058161191238065
    )

    xp = 602046282375688656758213480587526111916698976636884684818
    yp = 174050332293622031404857552280219410364023488927386650641

    @test octet(xp, yp, bitlength(prime_curve); mode = :compressed) == hex"03 188DA80E B03090F6 7CBF20EB 43A18800 F4FF0AFD 82FF1012"
    @test point(hex"03 188DA80E B03090F6 7CBF20EB 43A18800 F4FF0AFD 82FF1012", prime_curve) == (xp, yp)

    @test octet(xp, yp, bitlength(prime_curve); mode = :uncompressed) == hex"04 188DA80E B03090F6 7CBF20EB 43A18800 F4FF0AFD 82FF1012 07192B95 FFC8DA78 631011ED 6B24CDD5 73F977A1 1E794811"
    @test point(hex"04 188DA80E B03090F6 7CBF20EB 43A18800 F4FF0AFD 82FF1012 07192B95 FFC8DA78 631011ED 6B24CDD5 73F977A1 1E794811", prime_curve) == (xp, yp)

    @test octet(xp, yp, bitlength(prime_curve); mode = :hybrid) == hex"07 188DA80E B03090F6 7CBF20EB 43A18800 F4FF0AFD 82FF1012 07192B95 FFC8DA78 631011ED 6B24CDD5 73F977A1 1E794811"
    @test point(hex"07 188DA80E B03090F6 7CBF20EB 43A18800 F4FF0AFD 82FF1012 07192B95 FFC8DA78 631011ED 6B24CDD5 73F977A1 1E794811", prime_curve) == (xp, yp)

end

### Now the binary curve with polynomial basis
let

    f = octet2bits(hex"800000000000000000000000000000000000000000000201", 192)

    a = octet2bits(hex"2866537B 67675263 6A68F565 54E12640 276B649E F7526267", 191)
    b = octet2bits(hex"2866537B 67675263 6A68F565 54E12640 276B649E F7526267", 191)

    binary_curve = EC2N{PB}(
        f = f,
        a = a,
        b = b,
        n = 0,
    )

    xp = bin"01101101011001111011010111110001010001000110010000001101111100111000100111100101001100111010111101100100001101010011100001101101001000100110111111100101100100001001010111000011010101000001101"
    yp = bin"11101100101101111100111001101000011001110110011111110010101111000110011001010010011001011100111000011101010001001000101110010100010010000011000111010100000111011111001100000000001100011111011"

    @test octet(xp, yp, PB(f); mode = :uncompressed) == hex"04 36B3DAF8 A23206F9 C4F299D7 B21A9C36 9137F2C8 4AE1AA0D 765BE734 33B3F95E 332932E7 0EA245CA 2418EA0E F98018FB"
    @test point(hex"04 36B3DAF8 A23206F9 C4F299D7 B21A9C36 9137F2C8 4AE1AA0D 765BE734 33B3F95E 332932E7 0EA245CA 2418EA0E F98018FB", binary_curve) == (xp, yp)

    @test octet(xp, yp, PB(f); mode = :compressed) == hex"02 36B3DAF8 A23206F9 C4F299D7 B21A9C36 9137F2C8 4AE1AA0D"

    @test octet(xp, yp, PB(f); mode = :hybrid) == hex"06 36B3DAF8 A23206F9 C4F299D7 B21A9C36 9137F2C8 4AE1AA0D 765BE734 33B3F95E 332932E7 0EA245CA 2418EA0E F98018FB"
    @test point(hex"06 36B3DAF8 A23206F9 C4F299D7 B21A9C36 9137F2C8 4AE1AA0D 765BE734 33B3F95E 332932E7 0EA245CA 2418EA0E F98018FB", binary_curve) == (xp, yp)

end

### Now reading of the point
let

    prime_curve = ECP(
        p = 6277101735386680763835789423207666416083908700390324961279,
        a = 6277101735386680763835789423207666416083908700390324961276,
        b = 5005402392289390203552069470771117084861899307801456990547
    )

    @test point(hex"03 EEA2BAE7 E1497842 F2DE7769 CFE9C989 C072AD69 6F48034A", prime_curve) == (5851329466723574623122023978072381191095567081251774399306, 2487701625881228691269808880535093938601070911264778280469)

end

### Last part of converting binary curve to a compressed point

let

    f = octet2bits(hex"800000000000000000000000000000000000000000000201", 192)

    a = octet2bits(hex"40102877 4D7777C7 B7666D13 66EA4320 71274F89 FF01E718", 191)
    b = octet2bits(hex"0620048D 28BCBD03 B6249C99 182B7C8C D19700C3 62C46A01", 191)

    binary_curve = EC2N{PB}(
        f = f,
        a = a,
        b = b,
    )

    PO = hex"02 3809B2B7 CC1B28CC 5A87926A AD83FD28 789E81E2 C9E3BF10"

    xp = bin"01110000000100110110010101101111100110000011011001010001100110001011010100001111001001001101010101011011000001111111101001010000111100010011110100000011110001011001001111000111011111100010000"
    yp = bin"00101110100001101000011100001100110001001101101000101001111001111011011111100000001011101100000110110010010000100111010001111100001110011110011011110101110110001000011011111010110011010001010"

    # point(_hex2bytes(PO), binary_curve) == (xp, yp)
end

using CryptoGroups: spec, specialize, ECPoint, AffinePoint, Weierstrass, FP, generator, <|, F2PB, BinaryCurve, F2GNB

import CryptoGroups

let

    _spec = spec(:P_192)

    P = specialize(ECPoint{AffinePoint{Weierstrass, FP}}, _spec)
    g = P <| generator(_spec)
    @test P <| octet(g; mode = :uncompressed) == g
    @test P <| octet(g; mode = :compressed) == g
    @test P <| octet(g; mode = :hybrid) == g

end

let

    _spec = CryptoGroups.Specs.Curve_B_163_PB

    P = specialize(ECPoint{AffinePoint{BinaryCurve, F2PB}}, _spec)
    g = P <| generator(_spec)
    @test P <| octet(g; mode = :uncompressed) == g
    #@test P <| octet(g; mode = :compressed) == g
    @test P <| octet(g; mode = :hybrid) == g
end

let
    _spec = CryptoGroups.Specs.Curve_B_163_GNB

    P = specialize(ECPoint{AffinePoint{BinaryCurve, F2GNB}}, _spec)
    g = P <| generator(_spec)
    @test P <| octet(g; mode = :uncompressed) == g
    #@test P <| octet(g; mode = :compressed) == g
    @test P <| octet(g; mode = :hybrid) == g
end

let

    _spec = CryptoGroups.Specs.Curve_K_163_PB

    P = specialize(ECPoint{AffinePoint{BinaryCurve, F2PB}}, _spec)
    g = P <| generator(_spec)
    @test P <| octet(g; mode = :uncompressed) == g
    #@test P <| octet(g; mode = :compressed) == g
    @test P <| octet(g; mode = :hybrid) == g
end


text:
module Scripts

  using ..AlphaZero
  using ..AlphaZero.UserInterface: Session, resume!

  include("dummy_run.jl")
  export dummy_run

  dummy_run(s::String; args...) = dummy_run(Examples.experiments[s]; args...)

  include("test_game.jl")
  export test_game

  test_game(e::Experiment; args...) = test_game(e.gspec; args...)

  test_game(s::String; args...) = test_game(Examples.experiments[s]; args...)

  """
      train(experiment; [dir, autosave, save_intermediate])

  Start or resume a training session.

  The optional keyword arguments are passed
  directly to the [`Session`](@ref Session(::Experiment)) constructor.
  """
  train(e::Experiment; args...) = UserInterface.resume!(Session(e; args...))

  train(s::String; args...) = train(Examples.experiments[s]; args...)

  """
      explore(experiment; [dir])

  Use the interactive explorer to visualize the current agent.
  """
  explore(e::Experiment; args...) = UserInterface.explore(Session(e; args...))

  explore(s::String; args...) = explore(Examples.experiments[s]; args...)

  function play(e::Experiment; args...)
    session = Session(e; args...)
    if GI.two_players(e.gspec)
      interactive!(session.env.gspec, AlphaZeroPlayer(session), Human())
    else
      interactive!(session.env.gspec, Human())
    end
  end

  """
      play(experiment; [dir])

  Play an interactive game against the current agent.
  """
  play(s::String; args...) = play(Examples.experiments[s]; args...)

  include("test_grad_updates.jl")

  test_grad_updates(s::String; args...) =
    test_grad_updates(Examples.experiments[s]; args...)

end