text:
@doc raw"""
    AbstractGroupOperation

Abstract type for smooth binary operations $∘$ on elements of a Lie group $\mathcal{G}$:
```math
∘ : \mathcal{G} × \mathcal{G} → \mathcal{G}
```
An operation can be either defined for a specific [`AbstractGroupManifold`](@ref) over
number system `𝔽` or in general, by defining for an operation `Op` the following methods:

    identity!(::AbstractGroupManifold{𝔽,Op}, q, q)
    identity(::AbstractGroupManifold{𝔽,Op}, p)
    inv!(::AbstractGroupManifold{𝔽,Op}, q, p)
    inv(::AbstractGroupManifold{𝔽,Op}, p)
    compose(::AbstractGroupManifold{𝔽,Op}, p, q)
    compose!(::AbstractGroupManifold{𝔽,Op}, x, p, q)

Note that a manifold is connected with an operation by wrapping it with a decorator,
[`AbstractGroupManifold`](@ref). In typical cases the concrete wrapper
[`GroupManifold`](@ref) can be used.
"""
abstract type AbstractGroupOperation end

@doc raw"""
    AbstractGroupManifold{𝔽,O<:AbstractGroupOperation} <: AbstractDecoratorManifold{𝔽}

Abstract type for a Lie group, a group that is also a smooth manifold with an
[`AbstractGroupOperation`](@ref), a smooth binary operation. `AbstractGroupManifold`s must
implement at least [`inv`](@ref), [`identity`](@ref), [`compose`](@ref), and
[`translate_diff`](@ref).
"""
abstract type AbstractGroupManifold{𝔽,O<:AbstractGroupOperation,T<:AbstractEmbeddingType} <:
              AbstractEmbeddedManifold{𝔽,T} end

"""
    GroupManifold{𝔽,M<:Manifold{𝔽},O<:AbstractGroupOperation} <: AbstractGroupManifold{𝔽,O}

Decorator for a smooth manifold that equips the manifold with a group operation, thus making
it a Lie group. See [`AbstractGroupManifold`](@ref) for more details.

Group manifolds by default forward metric-related operations to the wrapped manifold.

# Constructor

    GroupManifold(manifold, op)
"""
struct GroupManifold{𝔽,M<:Manifold{𝔽},O<:AbstractGroupOperation} <:
       AbstractGroupManifold{𝔽,O,TransparentIsometricEmbedding}
    manifold::M
    op::O
end

Base.show(io::IO, G::GroupManifold) = print(io, "GroupManifold($(G.manifold), $(G.op))")

const GROUP_MANIFOLD_BASIS_DISAMBIGUATION =
    [AbstractDecoratorManifold, ValidationManifold, VectorBundle]

"""
    base_group(M::Manifold) -> AbstractGroupManifold

Un-decorate `M` until an `AbstractGroupManifold` is encountered.
Return an error if the [`base_manifold`](@ref) is reached without encountering a group.
"""
base_group(M::AbstractDecoratorManifold) = base_group(decorated_manifold(M))
function base_group(M::Manifold)
    return error("base_group: no base group found.")
end
base_group(G::AbstractGroupManifold) = G

base_manifold(G::GroupManifold) = G.manifold

decorator_group_dispatch(::Manifold) = Val(false)
function decorator_group_dispatch(M::AbstractDecoratorManifold)
    return decorator_group_dispatch(decorated_manifold(M))
end
decorator_group_dispatch(::AbstractGroupManifold) = Val(true)

function is_group_decorator(M::Manifold)
    return _extract_val(decorator_group_dispatch(M))
end

default_decorator_dispatch(::AbstractGroupManifold) = Val(false)

# piping syntax for decoration
if VERSION ≥ v"1.3"
    (op::AbstractGroupOperation)(M::Manifold) = GroupManifold(M, op)
    (::Type{T})(M::Manifold) where {T<:AbstractGroupOperation} = GroupManifold(M, T())
end

function decorator_transparent_dispatch(
    ::typeof(get_coordinates!),
    ::AbstractGroupManifold{<:Any,<:AbstractGroupOperation,TransparentIsometricEmbedding},
    args...,
)
    return Val(:transparent)
end
function decorator_transparent_dispatch(
    ::typeof(get_vector!),
    ::AbstractGroupManifold{<:Any,<:AbstractGroupOperation,TransparentIsometricEmbedding},
    args...,
)
    return Val(:transparent)
end

###################
# Action directions
###################

"""
    ActionDirection

Direction of action on a manifold, either [`LeftAction`](@ref) or [`RightAction`](@ref).
"""
abstract type ActionDirection end

"""
    LeftAction()

Left action of a group on a manifold.
"""
struct LeftAction <: ActionDirection end

"""
    RightAction()

Right action of a group on a manifold.
"""
struct RightAction <: ActionDirection end

"""
    switch_direction(::ActionDirection)

Returns a [`RightAction`](@ref) when given a [`LeftAction`](@ref) and vice versa.
"""
switch_direction(::ActionDirection)
switch_direction(::LeftAction) = RightAction()
switch_direction(::RightAction) = LeftAction()

##################################
# General Identity element methods
##################################

@doc raw"""
    Identity(G::AbstractGroupManifold, p)

The group identity element $e ∈ \mathcal{G}$ represented by point `p`.
"""
struct Identity{G<:AbstractGroupManifold,PT}
    group::G
    p::PT
end

Identity(M::AbstractDecoratorManifold, p) = Identity(decorated_manifold(M), p)
function Identity(M::Manifold, p)
    return error("Identity not implemented for manifold $(M) and point $(p).")
end

function Base.:(==)(e1::Identity, e2::Identity)
    return e1.p == e2.p && e1.group == e2.group
end

make_identity(M::Manifold, p) = Identity(M, identity(M, p))

Base.show(io::IO, e::Identity) = print(io, "Identity($(e.group), $(e.p))")

# To ensure allocate_result_type works
number_eltype(e::Identity) = Bool

Base.copyto!(e::TE, ::TE) where {TE<:Identity} = e
Base.copyto!(p, ::TE) where {TE<:Identity} = copyto!(p, e.p)
Base.copyto!(p::AbstractArray, e::TE) where {TE<:Identity} = copyto!(p, e.p)

Base.isapprox(p, e::Identity; kwargs...) = isapprox(e::Identity, p; kwargs...)
Base.isapprox(e::Identity, p; kwargs...) = isapprox(e.group, e, p; kwargs...)
Base.isapprox(e::E, ::E; kwargs...) where {E<:Identity} = true

function allocate_result(
    M::Manifold,
    f::typeof(get_coordinates),
    e::Identity,
    X,
    B::AbstractBasis,
)
    T = allocate_result_type(M, f, (e.p, X))
    return allocate(e.p, T, Size(number_of_coordinates(M, B)))
end

function decorator_transparent_dispatch(
    ::typeof(get_vector),
    ::AbstractGroupManifold,
    args...,
)
    return Val(:parent)
end
function allocate_result(M::Manifold, f::typeof(get_vector), e::Identity, Xⁱ)
    is_group_decorator(M) && return allocate_result(base_group(M), f, e, Xⁱ)
    return error("allocate_result not implemented for manifold $(M), function $(f), point $(e), and vector $(Xⁱ).")
end
function allocate_result(M::AbstractGroupManifold, f::typeof(get_vector), e::Identity, Xⁱ)
    return error("allocate_result not implemented for group manifold $(M), function $(f), $(e), and vector $(Xⁱ).")
end
function allocate_result(
    G::GT,
    ::typeof(get_vector),
    ::Identity{GT},
    Xⁱ,
) where {GT<:AbstractGroupManifold}
    B = VectorBundleFibers(TangentSpace, G)
    return allocate(Xⁱ, Size(representation_size(B)))
end

function decorator_transparent_dispatch(
    ::typeof(get_coordinates),
    ::AbstractGroupManifold,
    args...,
)
    return Val(:parent)
end
function allocate_result(
    M::AbstractDecoratorManifold,
    f::typeof(get_coordinates),
    e::Identity,
    X,
)
    is_group_decorator(M) && return allocate_result(base_group(M), f, e, X)
    return error("allocate_result not implemented for manifold $(M), function $(f), point $(e), and vector $(X).")
end
function allocate_result(
    M::AbstractGroupManifold,
    f::typeof(get_coordinates),
    e::Identity,
    X,
)
    return error("allocate_result not implemented for group manifold $(M), function $(f), $(e), and vector $(X).")
end
function allocate_result(
    G::GT,
    ::typeof(get_coordinates),
    ::Identity{GT},
    X,
) where {GT<:AbstractGroupManifold}
    return allocate(X, Size(manifold_dimension(G)))
end

function get_vector(M::AbstractGroupManifold, e::Identity, X, B::VeeOrthogonalBasis)
    M != e.group && error("On $(M) the identity $(e) does not match to perform get_vector.")
    return get_vector(decorated_manifold(M), e.p, X, B)
end
function get_vector(M::Manifold, e::Identity, X, B::VeeOrthogonalBasis)
    M != e.group.manifold &&
        error("On $(M) the identity $(e) does not match to perform get_vector.")
    return get_vector(M, e.p, X, B)
end
for MT in GROUP_MANIFOLD_BASIS_DISAMBIGUATION
    eval(
        quote
            @invoke_maker 1 Manifold get_vector(
                M::$MT,
                e::Identity,
                X,
                B::VeeOrthogonalBasis,
            )
        end,
    )
end
function get_vector!(M::AbstractGroupManifold, Y, e::Identity, X, B::VeeOrthogonalBasis)
    M != e.group && error("On $(M) the identity $(e) does not match to perform get_vector!")
    return get_vector!(decorated_manifold(M), Y, e.p, X, B)
end
function get_vector!(M::Manifold, Y, e::Identity, X, B::VeeOrthogonalBasis)
    M != e.group.manifold &&
        error("On $(M) the identity $(e) does not match to perform get_vector!")
    return get_vector!(M, Y, e.p, X, B)
end
for MT in GROUP_MANIFOLD_BASIS_DISAMBIGUATION
    eval(
        quote
            @invoke_maker 1 Manifold get_vector!(
                M::$MT,
                Y,
                e::Identity,
                X,
                B::VeeOrthogonalBasis,
            )
        end,
    )
end

function get_coordinates(M::AbstractGroupManifold, e::Identity, X, B::VeeOrthogonalBasis)
    M != e.group &&
        error("On $(M) the identity $(e) does not match to perform get_coordinates")
    return get_coordinates(decorated_manifold(M), e.p, X, B)
end
function get_coordinates(M::Manifold, e::Identity, X, B::VeeOrthogonalBasis)
    M != e.group.manifold &&
        error("On $(M) the identity $(e) does not match to perform get_coordinates")
    return get_coordinates(M, e.p, X, B)
end
for MT in GROUP_MANIFOLD_BASIS_DISAMBIGUATION
    eval(
        quote
            @invoke_maker 1 Manifold get_coordinates(
                M::$MT,
                e::Identity,
                X,
                B::VeeOrthogonalBasis,
            )
        end,
    )
end

function get_coordinates!(
    M::AbstractGroupManifold,
    Y,
    e::Identity,
    X,
    B::VeeOrthogonalBasis,
)
    M != e.group &&
        error("On $(M) the identity $(e) does not match to perform get_coordinates!")
    return get_coordinates!(decorated_manifold(M), Y, e.p, X, B)
end
function get_coordinates!(M::Manifold, Y, e::Identity, X, B::VeeOrthogonalBasis)
    M != e.group.manifold &&
        error("On $(M) the identity $(e) does not match to perform get_coordinates!")
    return get_coordinates!(M, Y, e.p, X, B)
end
for MT in GROUP_MANIFOLD_BASIS_DISAMBIGUATION
    eval(
        quote
            @invoke_maker 1 Manifold get_coordinates!(
                M::$MT,
                Y,
                e::Identity,
                X,
                B::VeeOrthogonalBasis,
            )
        end,
    )
end

manifold_dimension(G::GroupManifold) = manifold_dimension(G.manifold)

function check_manifold_point(G::AbstractGroupManifold, e::Identity; kwargs...)
    return DomainError(e, "The identity element $(e) does not belong to $(G).")
end
function check_manifold_point(
    G::GT,
    e::Identity{GT};
    kwargs...,
) where {GT<:AbstractGroupManifold}
    return nothing
end
##########################
# Group-specific functions
##########################

@doc raw"""
    inv(G::AbstractGroupManifold, p)

Inverse $p^{-1} ∈ \mathcal{G}$ of an element $p ∈ \mathcal{G}$, such that
$p \circ p^{-1} = p^{-1} \circ p = e ∈ \mathcal{G}$, where $e$ is the [`identity`](@ref)
element of $\mathcal{G}$.
"""
inv(::AbstractGroupManifold, ::Any...)
@decorator_transparent_function :intransparent function Base.inv(
    G::AbstractGroupManifold,
    p,
)
    q = allocate_result(G, inv, p)
    return inv!(G, q, p)
end

@decorator_transparent_function :intransparent function inv!(G::AbstractGroupManifold, q, p)
    return inv!(G.manifold, q, p)
end

@doc raw"""
    identity(G::AbstractGroupManifold, p)

Identity element $e ∈ \mathcal{G}$, such that for any element $p ∈ \mathcal{G}$,
$p \circ e = e \circ p = p$.
The returned element is of a similar type to `p`.
"""
identity(::AbstractGroupManifold, ::Any)
@decorator_transparent_signature Base.identity(G::AbstractDecoratorManifold, p)
function Base.identity(G::AbstractGroupManifold, p)
    y = allocate_result(G, identity, p)
    return identity!(G, y, p)
end

@decorator_transparent_signature identity!(G::AbstractDecoratorManifold, q, p)
function decorator_transparent_dispatch(::typeof(identity!), ::AbstractGroupManifold, q, p)
    return Val(:intransparent)
end

function Base.isapprox(
    G::GT,
    e::Identity{GT},
    p;
    kwargs...,
) where {GT<:AbstractGroupManifold}
    return isapprox(G, e.p, p; kwargs...)
end
function Base.isapprox(
    G::GT,
    p,
    e::Identity{GT};
    kwargs...,
) where {GT<:AbstractGroupManifold}
    return isapprox(G, e, p; kwargs...)
end
function Base.isapprox(
    ::GT,
    ::E,
    ::E;
    kwargs...,
) where {GT<:AbstractGroupManifold,E<:Identity{GT}}
    return true
end

function decorator_transparent_dispatch(
    ::typeof(isapprox),
    M::AbstractDecoratorManifold,
    e::E,
    p,
) where {E<:Identity{<:AbstractGroupManifold}}
    return Val(:transparent)
end
function decorator_transparent_dispatch(
    ::typeof(isapprox),
    M::AbstractDecoratorManifold,
    p,
    e::E,
) where {E<:Identity{<:AbstractGroupManifold}}
    return Val(:transparent)
end
function decorator_transparent_dispatch(
    ::typeof(isapprox),
    M::AbstractDecoratorManifold,
    e1::E1,
    e2::E2,
) where {E1<:Identity{<:AbstractGroupManifold},E2<:Identity{<:AbstractGroupManifold}}
    return Val(:transparent)
end

@doc raw"""
    compose(G::AbstractGroupManifold, p, q)

Compose elements $p,q ∈ \mathcal{G}$ using the group operation $p \circ q$.
"""
compose(::AbstractGroupManifold, ::Any...)
@decorator_transparent_signature compose(M::AbstractDecoratorManifold, p, q)
function compose(G::AbstractGroupManifold, p, q)
    x = allocate_result(G, compose, p, q)
    return compose!(G, x, p, q)
end
function decorator_transparent_dispatch(
    ::typeof(compose),
    M::AbstractGroupManifold,
    x,
    p,
    q,
)
    return Val(:intransparent)
end

@decorator_transparent_signature compose!(M::AbstractDecoratorManifold, x, p, q)
function decorator_transparent_dispatch(
    ::typeof(compose!),
    M::AbstractGroupManifold,
    x,
    p,
    q,
)
    return Val(:intransparent)
end

_action_order(p, q, conv::LeftAction) = (p, q)
_action_order(p, q, conv::RightAction) = (q, p)

@doc raw"""
    translate(G::AbstractGroupManifold, p, q)
    translate(G::AbstractGroupManifold, p, q, conv::ActionDirection=LeftAction()])

Translate group element $q$ by $p$ with the translation $τ_p$ with the specified
`conv`ention, either left ($L_p$) or right ($R_p$), defined as
```math
\begin{aligned}
L_p &: q ↦ p \circ q\\
R_p &: q ↦ q \circ p.
\end{aligned}
```
"""
translate(::AbstractGroupManifold, ::Any...)
@decorator_transparent_signature translate(M::AbstractDecoratorManifold, p, q)
function translate(G::AbstractGroupManifold, p, q)
    return translate(G, p, q, LeftAction())
end
@decorator_transparent_signature translate(
    M::AbstractDecoratorManifold,
    p,
    q,
    conv::ActionDirection,
)
function translate(G::AbstractGroupManifold, p, q, conv::ActionDirection)
    return compose(G, _action_order(p, q, conv)...)
end

function decorator_transparent_dispatch(
    ::typeof(translate),
    ::AbstractGroupManifold,
    args...,
)
    return Val(:intransparent)
end

@decorator_transparent_signature translate!(M::AbstractDecoratorManifold, X, p, q)
function translate!(G::AbstractGroupManifold, X, p, q)
    return translate!(G, X, p, q, LeftAction())
end
@decorator_transparent_signature translate!(
    M::AbstractDecoratorManifold,
    X,
    p,
    q,
    conv::ActionDirection,
)
function translate!(G::AbstractGroupManifold, X, p, q, conv::ActionDirection)
    return compose!(G, X, _action_order(p, q, conv)...)
end

function decorator_transparent_dispatch(
    ::typeof(translate!),
    ::AbstractGroupManifold,
    args...,
)
    return Val(:intransparent)
end

@doc raw"""
    inverse_translate(G::AbstractGroupManifold, p, q)
    inverse_translate(G::AbstractGroupManifold, p, q, conv::ActionDirection=LeftAction())

Inverse translate group element $q$ by $p$ with the inverse translation $τ_p^{-1}$ with the
specified `conv`ention, either left ($L_p^{-1}$) or right ($R_p^{-1}$), defined as
```math
\begin{aligned}
L_p^{-1} &: q ↦ p^{-1} \circ q\\
R_p^{-1} &: q ↦ q \circ p^{-1}.
\end{aligned}
```
"""
inverse_translate(::AbstractGroupManifold, ::Any...)
@decorator_transparent_function :intransparent function inverse_translate(
    G::AbstractGroupManifold,
    p,
    q,
)
    return inverse_translate(G, p, q, LeftAction())
end
@decorator_transparent_function :intransparent function inverse_translate(
    G::AbstractGroupManifold,
    p,
    q,
    conv::ActionDirection,
)
    return translate(G, inv(G, p), q, conv)
end

@decorator_transparent_function :intransparent function inverse_translate!(
    G::AbstractGroupManifold,
    X,
    p,
    q,
)
    return inverse_translate!(G, X, p, q, LeftAction())
end
@decorator_transparent_function :intransparent function inverse_translate!(
    G::AbstractGroupManifold,
    X,
    p,
    q,
    conv::ActionDirection,
)
    return translate!(G, X, inv(G, p), q, conv)
end

@doc raw"""
    translate_diff(G::AbstractGroupManifold, p, q, X)
    translate_diff(G::AbstractGroupManifold, p, q, X, conv::ActionDirection=LeftAction())

For group elements $p, q ∈ \mathcal{G}$ and tangent vector $X ∈ T_q \mathcal{G}$, compute
the action of the differential of the translation $τ_p$ by $p$ on $X$, with the specified
left or right `conv`ention. The differential transports vectors:
```math
(\mathrm{d}τ_p)_q : T_q \mathcal{G} → T_{τ_p q} \mathcal{G}\\
```
"""
translate_diff(::AbstractGroupManifold, ::Any...)
@decorator_transparent_signature translate_diff(M::AbstractDecoratorManifold, p, q, X)
function translate_diff(G::AbstractGroupManifold, p, q, X)
    return translate_diff(G, p, q, X, LeftAction())
end
@decorator_transparent_signature translate_diff(
    M::AbstractDecoratorManifold,
    p,
    q,
    X,
    conv::ActionDirection,
)
function translate_diff(G::AbstractGroupManifold, p, q, X, conv::ActionDirection)
    Y = allocate_result(G, translate_diff, X, p, q)
    translate_diff!(G, Y, p, q, X, conv)
    return Y
end

@decorator_transparent_signature translate_diff!(M::AbstractDecoratorManifold, Y, p, q, X)
function translate_diff!(G::AbstractGroupManifold, Y, p, q, X)
    return translate_diff!(G, Y, p, q, X, LeftAction())
end
@decorator_transparent_signature translate_diff!(
    M::AbstractDecoratorManifold,
    Y,
    p,
    q,
    X,
    conv::ActionDirection,
)
function decorator_transparent_dispatch(
    ::typeof(translate_diff!),
    G::AbstractGroupManifold,
    Y,
    p,
    q,
    X,
    conv,
)
    return Val(:intransparent)
end

@doc raw"""
    inverse_translate_diff(G::AbstractGroupManifold, p, q, X)
    inverse_translate_diff(G::AbstractGroupManifold, p, q, X, conv::ActionDirection=LeftAction())

For group elements $p, q ∈ \mathcal{G}$ and tangent vector $X ∈ T_q \mathcal{G}$, compute
the action on $X$ of the differential of the inverse translation $τ_p$ by $p$, with the
specified left or right `conv`ention. The differential transports vectors:
```math
(\mathrm{d}τ_p^{-1})_q : T_q \mathcal{G} → T_{τ_p^{-1} q} \mathcal{G}\\
```
"""
inverse_translate_diff(::AbstractGroupManifold, ::Any...)
@decorator_transparent_function :intransparent function inverse_translate_diff(
    G::AbstractGroupManifold,
    p,
    q,
    X,
)
    return inverse_translate_diff(G, p, q, X, LeftAction())
end
@decorator_transparent_function :intransparent function inverse_translate_diff(
    G::AbstractGroupManifold,
    p,
    q,
    X,
    conv::ActionDirection,
)
    return translate_diff(G, inv(G, p), q, X, conv)
end

@decorator_transparent_function :intransparent function inverse_translate_diff!(
    G::AbstractGroupManifold,
    Y,
    p,
    q,
    X,
)
    return inverse_translate_diff!(G, Y, p, q, X, LeftAction())
end
@decorator_transparent_function :intransparent function inverse_translate_diff!(
    G::AbstractGroupManifold,
    Y,
    p,
    q,
    X,
    conv::ActionDirection,
)
    return translate_diff!(G, Y, inv(G, p), q, X, conv)
end

@doc raw"""
    group_exp(G::AbstractGroupManifold, X)

Compute the group exponential of the Lie algebra element `X`.

Given an element $X ∈ 𝔤 = T_e \mathcal{G}$, where $e$ is the [`identity`](@ref) element of
the group $\mathcal{G}$, and $𝔤$ is its Lie algebra, the group exponential is the map

````math
\exp : 𝔤 → \mathcal{G},
````
such that for $t,s ∈ ℝ$, $γ(t) = \exp (t X)$ defines a one-parameter subgroup with the
following properties:

````math
\begin{aligned}
γ(t) &= γ(-t)^{-1}\\
γ(t + s) &= γ(t) \circ γ(s) = γ(s) \circ γ(t)\\
γ(0) &= e\\
\lim_{t → 0} \frac{d}{dt} γ(t) &= X.
\end{aligned}
````

!!! note
    In general, the group exponential map is distinct from the Riemannian exponential map
    [`exp`](@ref).

```
group_exp(G::AbstractGroupManifold{𝔽,AdditionOperation}, X) where {𝔽}
```

Compute $q = X$.

    group_exp(G::AbstractGroupManifold{𝔽,MultiplicationOperation}, X) where {𝔽}

For `Number` and `AbstractMatrix` types of `X`, compute the usual numeric/matrix
exponential,

````math
\exp X = \operatorname{Exp} X = \sum_{n=0}^∞ \frac{1}{n!} X^n.
````
"""
group_exp(::AbstractGroupManifold, ::Any...)
@decorator_transparent_signature group_exp(G::AbstractDecoratorManifold, X)
function group_exp(G::AbstractGroupManifold, X)
    q = allocate_result(G, group_exp, X)
    return group_exp!(G, q, X)
end

@decorator_transparent_signature group_exp!(G::AbstractDecoratorManifold, q, X)

function decorator_transparent_dispatch(::typeof(group_exp), ::AbstractGroupManifold, q, X)
    return Val(:intransparent)
end
function decorator_transparent_dispatch(::typeof(group_exp!), ::AbstractGroupManifold, q, X)
    return Val(:intransparent)
end

@doc raw"""
    group_log(G::AbstractGroupManifold, q)

Compute the group logarithm of the group element `q`.

Given an element $q ∈ \mathcal{G}$, compute the right inverse of the group exponential map
[`group_exp`](@ref), that is, the element $\log q = X ∈ 𝔤 = T_e \mathcal{G}$, such that
$q = \exp X$

!!! note
    In general, the group logarithm map is distinct from the Riemannian logarithm map
    [`log`](@ref).

```
group_log(G::AbstractGroupManifold{𝔽,AdditionOperation}, q) where {𝔽}
```

Compute $X = q$.

    group_log(G::AbstractGroupManifold{𝔽,MultiplicationOperation}, q) where {𝔽}

For `Number` and `AbstractMatrix` types of `q`, compute the usual numeric/matrix logarithm:

````math
\log q = \operatorname{Log} q = \sum_{n=1}^∞ \frac{(-1)^{n+1}}{n} (q - e)^n,
````

where $e$ here is the [`identity`](@ref) element, that is, $1$ for numeric $q$ or the
identity matrix $I_m$ for matrix $q ∈ ℝ^{m × m}$.
"""
group_log(::AbstractGroupManifold, ::Any...)
@decorator_transparent_signature group_log(M::AbstractDecoratorManifold, q)
function group_log(G::AbstractGroupManifold, q)
    X = allocate_result(G, group_log, q)
    return group_log!(G, X, q)
end
function decorator_transparent_dispatch(::typeof(group_log), ::AbstractGroupManifold, X, q)
    return Val(:intransparent)
end

@decorator_transparent_signature group_log!(M::AbstractDecoratorManifold, X, q)
function decorator_transparent_dispatch(::typeof(group_log!), ::AbstractGroupManifold, X, q)
    return Val(:intransparent)
end
############################
# Group-specific Retractions
############################

"""
    GroupExponentialRetraction{D<:ActionDirection} <: AbstractRetractionMethod

Retraction using the group exponential [`group_exp`](@ref) "translated" to any point on the
manifold.

For more details, see
[`retract`](@ref retract(::GroupManifold, p, X, ::GroupExponentialRetraction)).

# Constructor

    GroupExponentialRetraction(conv::ActionDirection = LeftAction())
"""
struct GroupExponentialRetraction{D<:ActionDirection} <: AbstractRetractionMethod end

function GroupExponentialRetraction(conv::ActionDirection = LeftAction())
    return GroupExponentialRetraction{typeof(conv)}()
end

"""
    GroupLogarithmicInverseRetraction{D<:ActionDirection} <: AbstractInverseRetractionMethod

Retraction using the group logarithm [`group_log`](@ref) "translated" to any point on the
manifold.

For more details, see
[`inverse_retract`](@ref inverse_retract(::GroupManifold, p, q ::GroupLogarithmicInverseRetraction)).

# Constructor

    GroupLogarithmicInverseRetraction(conv::ActionDirection = LeftAction())
"""
struct GroupLogarithmicInverseRetraction{D<:ActionDirection} <:
       AbstractInverseRetractionMethod end

function GroupLogarithmicInverseRetraction(conv::ActionDirection = LeftAction())
    return GroupLogarithmicInverseRetraction{typeof(conv)}()
end

direction(::GroupExponentialRetraction{D}) where {D} = D()
direction(::GroupLogarithmicInverseRetraction{D}) where {D} = D()

@doc raw"""
    retract(
        G::AbstractGroupManifold,
        p,
        X,
        method::GroupExponentialRetraction{<:ActionDirection},
    )

Compute the retraction using the group exponential [`group_exp`](@ref) "translated" to any
point on the manifold.
With a group translation ([`translate`](@ref)) $τ_p$ in a specified direction, the
retraction is

````math
\operatorname{retr}_p = τ_p \circ \exp \circ (\mathrm{d}τ_p^{-1})_p,
````

where $\exp$ is the group exponential ([`group_exp`](@ref)), and $(\mathrm{d}τ_p^{-1})_p$ is
the action of the differential of inverse translation $τ_p^{-1}$ evaluated at $p$ (see
[`inverse_translate_diff`](@ref)).
"""
function retract(G::AbstractGroupManifold, p, X, method::GroupExponentialRetraction)
    conv = direction(method)
    Xₑ = inverse_translate_diff(G, p, p, X, conv)
    pinvq = group_exp(G, Xₑ)
    q = translate(G, p, pinvq, conv)
    return q
end

function retract!(G::AbstractGroupManifold, q, p, X, method::GroupExponentialRetraction)
    conv = direction(method)
    Xₑ = inverse_translate_diff(G, p, p, X, conv)
    pinvq = group_exp(G, Xₑ)
    return translate!(G, q, p, pinvq, conv)
end

@doc raw"""
    inverse_retract(
        G::AbstractGroupManifold,
        p,
        X,
        method::GroupLogarithmicInverseRetraction{<:ActionDirection},
    )

Compute the inverse retraction using the group logarithm [`group_log`](@ref) "translated"
to any point on the manifold.
With a group translation ([`translate`](@ref)) $τ_p$ in a specified direction, the
retraction is

````math
\operatorname{retr}_p^{-1} = (\mathrm{d}τ_p)_e \circ \log \circ τ_p^{-1},
````

where $\log$ is the group logarithm ([`group_log`](@ref)), and $(\mathrm{d}τ_p)_e$ is the
action of the differential of translation $τ_p$ evaluated at the identity element $e$
(see [`translate_diff`](@ref)).
"""
function inverse_retract(G::GroupManifold, p, q, method::GroupLogarithmicInverseRetraction)
    conv = direction(method)
    pinvq = inverse_translate(G, p, q, conv)
    Xₑ = group_log(G, pinvq)
    return translate_diff(G, p, Identity(G, p), Xₑ, conv)
end

function inverse_retract!(
    G::AbstractGroupManifold,
    X,
    p,
    q,
    method::GroupLogarithmicInverseRetraction,
)
    conv = direction(method)
    pinvq = inverse_translate(G, p, q, conv)
    Xₑ = group_log(G, pinvq)
    return translate_diff!(G, X, p, Identity(G, p), Xₑ, conv)
end

#################################
# Overloads for AdditionOperation
#################################

"""
    AdditionOperation <: AbstractGroupOperation

Group operation that consists of simple addition.
"""
struct AdditionOperation <: AbstractGroupOperation end

const AdditionGroup = AbstractGroupManifold{𝔽,AdditionOperation} where {𝔽}

Base.:+(e::Identity{G}) where {G<:AdditionGroup} = e
Base.:+(p::Identity{G}, ::Identity{G}) where {G<:AdditionGroup} = p
Base.:+(::Identity{G}, p) where {G<:AdditionGroup} = p
Base.:+(p, ::Identity{G}) where {G<:AdditionGroup} = p
Base.:+(e::E, ::E) where {G<:AdditionGroup,E<:Identity{G}} = e

Base.:-(e::Identity{G}) where {G<:AdditionGroup} = e
Base.:-(e::Identity{G}, ::Identity{G}) where {G<:AdditionGroup} = e
Base.:-(::Identity{G}, p) where {G<:AdditionGroup} = -p
Base.:-(p, ::Identity{G}) where {G<:AdditionGroup} = p
Base.:-(e::E, ::E) where {G<:AdditionGroup,E<:Identity{G}} = e

Base.:*(e::Identity{G}, p) where {G<:AdditionGroup} = e
Base.:*(p, e::Identity{G}) where {G<:AdditionGroup} = e
Base.:*(e::E, ::E) where {G<:AdditionGroup,E<:Identity{G}} = e

Base.zero(e::Identity{G}) where {G<:AdditionGroup} = e

Base.identity(::AdditionGroup, p) = zero(p)

identity!(::AdditionGroup, q, p) = fill!(q, 0)

Base.inv(::AdditionGroup, p) = -p

inv!(::AdditionGroup, q, p) = copyto!(q, -p)

compose(::AdditionGroup, p, q) = p + q

function compose!(::GT, x, p, q) where {GT<:AdditionGroup}
    p isa Identity{GT} && return copyto!(x, q)
    q isa Identity{GT} && return copyto!(x, p)
    x .= p .+ q
    return x
end

translate_diff(::AdditionGroup, p, q, X, ::ActionDirection) = X

translate_diff!(::AdditionGroup, Y, p, q, X, ::ActionDirection) = copyto!(Y, X)

inverse_translate_diff(::AdditionGroup, p, q, X, ::ActionDirection) = X

function inverse_translate_diff!(::AdditionGroup, Y, p, q, X, ::ActionDirection)
    return copyto!(Y, X)
end

group_exp(::AdditionGroup, X) = X

group_exp!(::AdditionGroup, q, X) = copyto!(q, X)

group_log(::AdditionGroup, q) = q

group_log!(::AdditionGroup, X, q) = copyto!(X, q)

#######################################
# Overloads for MultiplicationOperation
#######################################

"""
    MultiplicationOperation <: AbstractGroupOperation

Group operation that consists of multiplication.
"""
struct MultiplicationOperation <: AbstractGroupOperation end

const MultiplicationGroup = AbstractGroupManifold{𝔽,MultiplicationOperation} where {𝔽}

Base.:*(e::Identity{G}) where {G<:MultiplicationGroup} = e
Base.:*(::Identity{G}, p) where {G<:MultiplicationGroup} = p
Base.:*(p, ::Identity{G}) where {G<:MultiplicationGroup} = p
Base.:*(e::E, ::E) where {G<:MultiplicationGroup,E<:Identity{G}} = e
Base.:*(::Identity{<:MultiplicationGroup}, e::Identity{<:AdditionGroup}) = e

Base.:/(p, ::Identity{G}) where {G<:MultiplicationGroup} = p
Base.:/(::Identity{G}, p) where {G<:MultiplicationGroup} = inv(p)
Base.:/(e::E, ::E) where {G<:MultiplicationGroup,E<:Identity{G}} = e

Base.:\(p, ::Identity{G}) where {G<:MultiplicationGroup} = inv(p)
Base.:\(::Identity{G}, p) where {G<:MultiplicationGroup} = p
Base.:\(e::E, ::E) where {G<:MultiplicationGroup,E<:Identity{G}} = e

Base.inv(e::Identity{G}) where {G<:MultiplicationGroup} = e

Base.one(e::Identity{G}) where {G<:MultiplicationGroup} = e

Base.transpose(e::Identity{G}) where {G<:MultiplicationGroup} = e

LinearAlgebra.det(::Identity{<:MultiplicationGroup}) = 1

LinearAlgebra.mul!(q, e::Identity{G}, p) where {G<:MultiplicationGroup} = copyto!(q, p)
LinearAlgebra.mul!(q, p, e::Identity{G}) where {G<:MultiplicationGroup} = copyto!(q, p)
function LinearAlgebra.mul!(q, e::E, ::E) where {G<:MultiplicationGroup,E<:Identity{G}}
    return identity!(e.group, q, e)
end

Base.identity(::MultiplicationGroup, p) = one(p)

function identity!(G::GT, q, p) where {GT<:MultiplicationGroup}
    isa(p, Identity{GT}) || return copyto!(q, one(p))
    return error("identity! not implemented on $(typeof(G)) for points $(typeof(q)) and $(typeof(p))")
end
identity!(::MultiplicationGroup, q::AbstractMatrix, p) = copyto!(q, I)

Base.inv(::MultiplicationGroup, p) = inv(p)

inv!(G::MultiplicationGroup, q, p) = copyto!(q, inv(G, p))

compose(::MultiplicationGroup, p, q) = p * q

# TODO: x might alias with p or q, we might be able to optimize it if it doesn't.
compose!(::MultiplicationGroup, x, p, q) = copyto!(x, p * q)

inverse_translate(::MultiplicationGroup, p, q, ::LeftAction) = p \ q
inverse_translate(::MultiplicationGroup, p, q, ::RightAction) = q / p

function inverse_translate!(G::MultiplicationGroup, x, p, q, conv::ActionDirection)
    return copyto!(x, inverse_translate(G, p, q, conv))
end

function group_exp!(G::MultiplicationGroup, q, X)
    X isa Union{Number,AbstractMatrix} && return copyto!(q, exp(X))
    return error("group_exp! not implemented on $(typeof(G)) for vector $(typeof(X)) and element $(typeof(q)).")
end


text:
using SimpleANOVA
using Test
using DataFrames

@testset "SimpleANOVA Tests" begin
include("test_anova.jl")
include("test_show.jl")
end

ext: jl


text:
# Test windowing.

## Regular Signals:
for loc in ([20,30,40], [2s,3s,4s]), win in ((-1s,1s), (-10,10), (-1.04s, 1.04s))
    t = 0.1s:0.1s:10.0s
    sig = Signal(t, collect(1:100))
    reps = window(sig, loc, win)
    @test reps.time == [2.0s,3.0s,4.0s]
    @test reps[1].time == -1s:0.1s:1s
    @test length(reps) == 1
    @test length(reps[1]) == 3
    @test reps[1][1] == collect(10:30)
    @test reps[1][2] == collect(20:40)
    @test reps[1][3] == collect(30:50)
end

## Irregular signals are tougher to test the edge-cases (with 4 different cases)
# Start with basic behavior like Regular Signals
# for loc in ([20,30,40], [2s,3s,4s]), win in ((-1s,1s), (-10,10), (-1.04s, 1.04s))
#     t = collect(0.1s:0.1s:10.0s)
#     sig = Signal(t, collect(1:100))
#     reps = window(sig, loc, win)
#     @test reps.time == [2.0s,3.0s,4.0s]
#     @test length(reps) == 1
#     @test length(reps[1]) == 3
#     # Unlike regular signals, reps[1] is an array, not a Signal.
#     @test reps[1][1].time == t[10:30] - t[20]
#     @test reps[1][2].time == t[20:40] - t[30]
#     @test reps[1][3].time == t[30:50] - t[40]
#     @test reps[1][1][1] == collect(10:30)
#     @test reps[1][2][1] == collect(20:40)
#     @test reps[1][3][1] == collect(30:50)
# end


text:
## plotting code using Plots

"""

Plotting of implicit functions.

An implicit function or equation is defined in terms of a logical
condition and a function of two variables. These produce `Predicate`
objects.

Predicates may be plotted over a specified region, the default begin [-5,5] x [-5,5].

The algorithm, breaks the region into blocks. The ultimate resolution
is given by `W=2^n` and `H=2^m`. The algorithm used, due to Tupper, colors
region if it knows the predicate is true or false, and otherwise
resolves the region into 4 equal-sized subregions and test each subregion
again to determine if it is true, false, or still maybe. This repeats until `W` and `H`
can't be subdivided.

The `plot`ting of a predicate simply plots each block that is knows satisfies the predicate "black" and each ambiguous block "red." By taking `m` and `n` larger the graphs look less "blocky" but take more time to render.

For text-based plots, the `asciigraph` function is available.

Examples:
```
using Plots, ImplicitEquations
pyplot()

a,b = -1,2
f(x,y) = y^4 - x^4 + a*y^2 + b*x^2
plot(f == 0)

## trident of Newton
c,d,e,h = 1,1,1,1
f(x,y) = x*y
g(x,y) =c*x^3 + d*x^2 + e*x + h
plot(eq(f,g), title="Trident of Newton") ## aka f ⩵ g (using Unicode\Equal[tab])

## inequality
f(x,y)= (y-5)*cos(4*sqrt((x-4)^2 + y^2))
g(x,y) = x*sin(2*sqrt(x^2 + y^2))
r = f < g
plot(r, (-10, 10), (-10, 10), N=9, M=9)  # (xmin, xmax), (ymin, ymax),
```
"""
plot_implicit = nothing

## Helpers to convert
linterp(A,B,a,b,W) = (a + A/W*(b-a),a + B/W*(b-a))

function xyrange(u, L, R, B, T, W, H; offset=0)
    tuple(linterp(u.x.val.lo, u.x.val.hi + offset, L, R, W)...,
          linterp(u.y.val.lo, u.y.val.hi + offset, B, T, H)...)
end

function get_xs_ys(rs, L, R, B, T, W, H)
    xs = Float64[]
    ys = Float64[]
    for u in rs
        x0,x1,y0,y1 = ImplicitEquations.xyrange(u, L,R,B,T,W,H)
        append!(xs, [x0,x1,x1,x0,x0]); push!(xs, NaN)
        append!(ys, [y0,y0,y1,y1,y0]); push!(ys, NaN)
    end
    if length(xs) > 0
        pop!(xs)
        pop!(ys)
    else
        xs = Float64[NaN]
        ys = Float64[NaN]
    end

    xs, ys
end

## A plot recipe
## x, y describe region to plot over
## N, M give no. of pixels 2^N by 2^M
## red and black are used for colors.
@recipe function f(p::Predicate, x=(-5,5), y=(-5,5);
                   N=8,
                   M=8,              # oddly m as keyword fails. 9/8 too slow
                   red=nothing,      # or :red ...
                   black=:black
                   )

    L, R = extrema(x)
    B, T = extrema(y)

    W = 2^N
    H = 2^M

    r, b, w = ImplicitEquations.GRAPH(p, L, R, B, T, W, H)

    ## add red as a series
    if length(r) > 0 && red != nothing
        @series begin
            xs, ys = get_xs_ys(r, L, R, B, T, W, H)

            seriestype := :shape
            fillcolor := red
            linewidth := 0
            x := xs
            y := ys

            ()
        end
    end

    ## do black
    seriestype --> :shape
    xlims --> [L, R]
    ylims --> [B, T]
    legend --> false

    fillcolor --> black
    linewidth --> 0

    xs, ys = get_xs_ys(b, L, R, B, T, W, H)
    x --> xs
    y --> ys

    ()

end


text:
# ---
# title: 169. Majority Element
# id: problem169
# author: zhwang
# date: 2022-03-10
# difficulty: Easy
# categories: Array, Divide and Conquer, Bit Manipulation
# link: <https://leetcode.com/problems/majority-element/description/>
# hidden: true
# ---
#
# Given an array of size _n_ , find the majority element. The majority element
# is the element that appears **more than** `⌊ n/2 ⌋` times.
#
# You may assume that the array is non-empty and the majority element always
# exist in the array.
#
# **Example 1:**
#
#
#
#     Input: [3,2,3]
#     Output: 3
#
# **Example 2:**
#
#
#
#     Input: [2,2,1,1,1,2,2]
#     Output: 2
#
#
#
## @lc code=start
using LeetCode

function majority_element(nums::Vector{Int})::Int
    res = vote = 0
    for num in nums
        if vote == 0
            vote += 1
            res = num
        else
            vote += (res == num) ? 1 : -1
        end
    end
    return res
end
## @lc code=end


text:
using AlphaZero

include("games.jl")
GAME = get(ENV, "GAME", "tictactoe")
SelectedGame = GAME_MODULE[GAME]
using .SelectedGame: Game

env = MCTS.Env{Game}(MCTS.RolloutOracle{Game}())
computer = MctsPlayer(env, niters=1, timeout=1.0, τ=ConstSchedule(0.5))

#interactive!(Game(), Human{Game}(), computer)
start_explorer(Explorer(computer))


text:
module FGb

using PolynomialRings
using PolynomialRings: construct_monomial, variablesymbols
using PolynomialRings.MonomialOrderings: MonomialOrder
using PolynomialRings.Monomials: enumeratenz, AbstractMonomial
using PolynomialRings.Polynomials: NamedPolynomial, PolynomialOver, terms
using PolynomialRings.Terms: coefficient, monomial

include("LibFGb.jl")

struct FGbPolynomial{T<:Polynomial,N}
    p::Ptr{Void}
end

in_FGb = false

function FGb_with(f::Function, ::Type{P}, num_vars::Integer) where P<:Polynomial
    global in_FGb
    assert(!in_FGb)

    enter_INT()
    init_urgent(4,2,"DRLDRL",100000,0)
    log_file_handle = Ref{UInt}(0)
    init(1,1,0,log_file_handle)

    reset_coeffs(1)
    varnames = ["xx$n" for n=1:num_vars]

    reset_expos(num_vars, 0, varnames)

    internal_threads(8)
    mod_internal_threads(8)

    in_FGb = true
    res = try
        f(p -> convert(FGbPolynomial{P, num_vars}, p))
    finally
        in_FGb = false
        reset_memory()
        exit_INT()
    end
    res
end

FGb_with(f::Function, ::Type{P}) where P<:NamedPolynomial = FGb_with(f, P, length(variablesymbols(P)))

import Base: convert, show

function show(io::IO, f::FGbPolynomial{T}) where T
    print(io, "FGb(")
    show(io, convert(T,f))
    print(io, ")")
end

function convert(::Type{FGbPolynomial{T,N}}, f::T) where T<:Polynomial where N
    p = creat_poly(length(terms(f)))
    for (i,t) in enumerate(terms(f))
        exp = zeros(UInt32, N)
        for (j,e) in enumeratenz(monomial(t))
            exp[j] = e
        end
        set_expos2(p,i-1,exp,length(exp))
        set_coeff_gmp(p,i-1,coefficient(t))
    end
    full_sort_poly2(p)
    FGbPolynomial{T,N}(p)
end

function convert(::Type{T}, f::FGbPolynomial{T,N}) where T<:Polynomial where N
    n_terms = nb_terms(f.p)
    exponents = Vector{UInt32}(n_terms * N)
    coeff_ptrs = Vector{Ptr{BigInt}}(n_terms)

    code = export_poly_INT_gmp2(N, n_terms, coeff_ptrs, exponents, f.p)

    coefficients = map(unsafe_load, coeff_ptrs)

    sum( c * construct_monomial(T, exponents[(1+k*N):((k+1)*N)]) for (k,c) in zip(0:(n_terms-1), coefficients) )
end

function groebner(F::Vector{T}) where T <: FGbPolynomial
    OUTPUTSIZE=100000
    output_basis = Vector{Ptr{Void}}(OUTPUTSIZE)
    env = SFGB_Comp_Desc()
    cputime = Ref{Float64}(0)
    n = groebner(map(f->f.p, F), length(F), output_basis, 1, 0, cputime,0, -1, 0, env)
    [T(output_basis[i]) for i=1:n]
end

import PolynomialRings
import PolynomialRings: gröbner_basis

struct FGbAlgorithm <: PolynomialRings.Backends.Gröbner.Backend end

const ApplicableBaserings = Union{BigInt, Rational{BigInt}}
const ApplicablePolynomial = PolynomialOver{<:ApplicableBaserings}

function gröbner_basis(::FGbAlgorithm, ::MonomialOrder{:degrevlex}, polynomials::AbstractArray{<:ApplicablePolynomial})
    length(polynomials) == 0 && return polynomials
    integral_polynomials = [p for (p, _) in integral_fraction.(polynomials)]
    P = eltype(polynomials)
    PP = eltype(integral_polynomials)
    num_vars = maximum(max_variable_index, integral_polynomials)
    FGb_with(PP, num_vars) do FGbPolynomial
        G = map(FGbPolynomial, integral_polynomials)
        gr = groebner(G)
        map(g -> convert(P, convert(PP,g)), gr)
    end
end

PolynomialRings.Backends.Gröbner.set_default(FGbAlgorithm())

end


text:
using Printf
using BenchmarkTools

# Get the matrix dimensions from the command line argument.
n, = size(ARGS)
if (n < 1)
    println("Usage: matrixMult.jl N")
    println("       ---> Please specify the dimensions.")
    exit()
end

N = parse(Int, ARGS[1])

println("----------------------------")
println(@sprintf "Matrix dimensions: %d" N)
println("----------------------------")
println(" ")

# Multiply matrices A and B.
A = randn(N, N)
B = randn(N, N)

@btime C = A*B


text:
"""
    abstract type AbstractFunctionConversionBridge{F, S} <: AbstractBridge end

Bridge a constraint `G`-in-`S` into a constraint `F`-in-`S` where `F` and `G`
are equivalent representations of the same function. By convention, the
transformed function is stored in the `constraint` field.
"""
abstract type AbstractFunctionConversionBridge{F,S} <: AbstractBridge end

function MOI.get(
    model::MOI.ModelLike,
    attr::MOI.AbstractConstraintAttribute,
    bridge::AbstractFunctionConversionBridge,
)
    if invariant_under_function_conversion(attr)
        return MOI.get(model, attr, bridge.constraint)
    else
        throw(
            ArgumentError(
                "Bridge of type `$(typeof(bridge))` does not support accessing the attribute `$attr` because `MOIB.Constraint.invariant_under_function_conversion($attr)` returns `false`.",
            ),
        )
    end
end

function MOI.supports(
    model::MOI.ModelLike,
    attr::MOI.AbstractConstraintAttribute,
    ::Type{<:AbstractFunctionConversionBridge{F,S}},
) where {F,S}
    return invariant_under_function_conversion(attr) &&
           MOI.supports(model, attr, MOI.ConstraintIndex{F,S})
end
function MOI.set(
    model::MOI.ModelLike,
    attr::MOI.AbstractConstraintAttribute,
    bridge::AbstractFunctionConversionBridge,
    value,
)
    if invariant_under_function_conversion(attr)
        return MOI.set(model, attr, bridge.constraint, value)
    else
        throw(
            ArgumentError(
                "Bridge of type `$(typeof(bridge))` does not support setting the attribute `$attr` because `MOIB.Constraint.invariant_under_function_conversion($attr)` returns `false`.",
            ),
        )
    end
end

"""
    invariant_under_function_conversion(attr::MOI.AbstractConstraintAttribute)

Returns whether the value of the attribute does not change if the constraint
`F`-in-`S` is transformed into a constraint `G`-in-`S` where `F` and `G` are
equivalent representations of the same function. If it returns true, then
subtypes of [`Constraint.AbstractFunctionConversionBridge`](@ref) such as
[`Constraint.ScalarFunctionizeBridge`](@ref) and
[`Constraint.VectorFunctionizeBridge`](@ref) will automatically support
[`MOI.get`](@ref) and [`MOI.set`](@ref) for `attr`.
"""
invariant_under_function_conversion(::MOI.AbstractConstraintAttribute) = false

function invariant_under_function_conversion(
    ::Union{
        MOI.ConstraintSet,
        MOI.ConstraintBasisStatus,
        MOI.ConstraintPrimal,
        MOI.ConstraintPrimalStart,
        MOI.ConstraintDual,
        MOI.ConstraintDualStart,
    },
)
    return true
end

include("functionize.jl")
const ScalarFunctionize{T,OT<:MOI.ModelLike} =
    SingleBridgeOptimizer{ScalarFunctionizeBridge{T},OT}
const VectorFunctionize{T,OT<:MOI.ModelLike} =
    SingleBridgeOptimizer{VectorFunctionizeBridge{T},OT}
# TODO add affine -> quadratic conversion bridge


text:
using PyCall
using LinearAlgebra
using TensorOperations
using Printf
using Test
include("hf_utils.jl")

# Setup
pyscf = pyimport("pyscf")

mol = pyscf.gto.M(atom = """O          0.00000        0.00000        0.11779
  H          0.00000        0.75545       -0.47116
  H          0.00000       -0.75545       -0.47116""", basis = "631g",)

mf = pyscf.scf.RHF(mol).run()

s = mol.intor_symmetric("int1e_ovlp")
h = mol.intor("int1e_kin") + mol.intor("int1e_nuc")
g = permutedims(mol.intor("int2e", aosym = 1), [4, 3, 2, 1])
nelectron = mol.nelectron
n_ao = convert(UInt8, mol.nao_nr())
e_nuc = mf.energy_nuc()
# End setup

@testset "HF Utils" begin
    @testset "vhf tools" begin
        dm = rand(n_ao, n_ao)
        dm = dm' + dm # Needs to be symmetric

        dm_row_major = permutedims(dm, [2, 1])
        pyscf_j, pyscf_k = pyscf.scf.hf.get_jk(mol, dm_row_major)
        pyscf_vhf = pyscf.scf.hf.get_veff(mol, dm_row_major)

        J = get_coloumb(g, dm)
        K = get_exchange(g, dm)
        @test norm(pyscf_j - J) < 1e-13
        @test norm(pyscf_k - K) < 1e-13

        vhf = get_vhf(g, dm)
        @test norm(pyscf_vhf - vhf) < 1e-13
    end

    @testset "rdm tools" begin
        pyscf_rdm1 = mf.make_rdm1()
        dm1 = make_rdm1(mf.mo_occ, mf.mo_coeff)
        @test norm(pyscf_rdm1 - dm1) < 1e-12

        mo_occ = zeros(n_ao)
        mo_occ[1:nelectron÷2] .= 2
        mo_coeff = rand(n_ao, n_ao)
        rand_error = make_rdm1(mo_occ, mo_coeff) - mf.make_rdm1(mo_coeff, mo_occ)
        @test norm(rand_error) < 1e-12
    end

    @testset "fock utils" begin
        dm = mf.make_rdm1()

        pyscf_fock = mf.get_fock()
        my_fock = get_fock(h, g, dm)
        @test norm(my_fock - pyscf_fock) < 1e-12

        mo_energy, mo_coeff = solve_fock(h, g, dm, s)
        @test mf.mo_energy ≈ mo_energy atol = 1e-6

        # Make sure MOs are orthogonal
        # mo_ovlp = mo_coeff * inv(mo_coeff)
        # mo_ovlp = mf.mo_coeff * inv(mf.mo_coeff)
        mo_ovlp = mo_coeff' * s * mo_coeff
        display(mo_coeff)
        println()
        display(mf.mo_coeff)
        println()
        res = svd(mo_ovlp)
        println(res.S)
        # Fₐₒ = get_fock(h, g, dm)
        # @tensor fock_mo[i, j] := mo_coeff[μ, i] * s[μ, ν] * Fₐₒ[ν, ρ] * s[ρ, σ] * mo_coeff[σ, j]
        # println(diag(fock_mo))
        # println(mo_energy)
        # println(mf.mo_energy)

        mo_occ = zeros(n_ao)
        mo_occ[1:nelectron÷2] .= 2

        dm1 = make_rdm1(mo_occ, mo_coeff)
        @test norm(mf.make_rdm1() - dm1) < 1e-6
        @tensor nelec_test = scalar(dm1[μ, ν] * s[μ, ν])
        println(nelec_test)
        @tensor nelec_test = scalar(mf.make_rdm1()[μ, ν] * s[μ, ν])
        println(nelec_test)
    end

    @testset "energy tools" begin
        my_energy = calc_energy_elec(h, g, mf.make_rdm1())[1] + e_nuc
        @test my_energy ≈ mf.e_tot atol = 1e-6
    end
end

# Clean up
using Glob
rm.(glob("tmp*"))


text:
using CompileBot

botconfig = BotConfig(
  "AcuteML";
  yml_path = "SnoopCompile.yml"        # parse `os` and `version` from `SnoopCompile.yml`
)

snoop_bot(
  botconfig,
  "$(@__DIR__)/example_script.jl",
)

max_stars_repo_name: aminya/AML


text:
#####
##### Main script for the JuliaCon 2020 demo
#####

ENV["JULIA_CUDA_MEMORY_POOL"] = "split" # "binned" / "split"

# Enables running the script on a distant machine without an X server
ENV["GKSwstype"]="nul"

using AlphaZero

const DUMMY_RUN = false
include("../scripts/lib/dummy_run.jl")

include("../games/connect-four/main.jl")
using .ConnectFour: Game, Training

params, benchmark = Training.params, Training.benchmark
if DUMMY_RUN
  params, benchmark = dummy_run_params(params, benchmark)
end

session = Session(
  Game,
  Training.Network{Game},
  params,
  Training.netparams,
  benchmark=benchmark,
  dir="sessions/connect-four",
  autosave=true,
  save_intermediate=false)

resume!(session)


text:
function km(ind, T, r)
    perm = sortperm(vec(T))
    lastT = T[perm[1]]
    n = 0
    d = length(perm)
    prod = 1.0
    Fr = ones(length(r))
    ri = 1
    for (i, p) in enumerate(perm)
        if T[p] != lastT
            prod *= (1 - n/d)
            while ri <= length(r) && r[ri] <= T[p]
                Fr[ri] = 1-prod
                ri += 1
            end
            n = 0
            d = length(perm) - i + 1
            lastT = T[p]
        end
        n += ind[p]
    end
    prod *= (1 - n/d)
    while ri <= length(r) && r[ri] <= lastT
        Fr[ri] = 1-prod
        ri += 1
    end
    Fr
end


text:
p = parse.(Int, readlines("Dec1_data.txt"))

# Question 1
fun(x) = div(x, 3) - 2

# test cases
fun(1969)    # should be 654
fun(100756)  # should be 33583

# result 1
sum(fun, p)

# Question 2
function fun2(rem)
    ret = fun(rem)
    ret <= 0 ? 0 : ret += fun2(ret)
end

# test cases
fun2(1969)    # should be 966
fun2(100756)  # should be 100756

# result 2
sum(fun2, p)


text:
function forward(backend::GPUBackend, state::CropLayerState, inputs::Vector{Blob})
  crop_size = state.layer.crop_size
  for i = 1:length(inputs)
    input = inputs[i]
    output = state.blobs[i]
    if state.layer.random_crop
      w_off = abs(rand(Int)) % (size(input,1)-crop_size[1]+1)
      h_off = abs(rand(Int)) % (size(input,2)-crop_size[2]+1)
    else
      w_off = div(size(input,1)-crop_size[1], 2)
      h_off = div(size(input,2)-crop_size[2], 2)
    end
    w_off2 = size(input,1) - crop_size[1] - w_off
    h_off2 = size(input,2) - crop_size[2] - h_off

    if state.layer.random_mirror && rand(UInt)%2 == 0
      mirror = true;
    else
      mirror = false;
    end
    padded2dense!(backend, output, input, (w_off, h_off), (w_off2, h_off2), mirror)
  end
end


text:
# functions to decompose a proper polygon into triangles
export triangulate

"""
`alt_mod(k,n)` is `mod(k,n)` unless the result is zero,
in which case we return `n`.
"""
function alt_mod(k::Int, n::Int)
    r = mod(k,n)
    return r==0 ? n : r
end

"""
`check_diagonal(X::HPolygon, k::Int)` checks to see if the diagonal
between vertices `k-1` and `k+1` (a) does not intersect the boundary of `X`
and (b) goes through the interior of `X`.
"""
function check_diagonal(X::HPolygon, k::Int)::Bool
    n = npoints(X)

    if n <= 3
        return true
    end

    kprev = alt_mod(k-1,n)
    knext = alt_mod(k+1,n)

    a = X.plist[kprev]
    b = X.plist[k]      # this is the current vertex (index k)
    c = X.plist[knext]

    dg = a+c  # this is the diagonal we're testing

    # check if ac contains a vertex (other than a or c)
    for i=1:n
        if i != kprev && i != knext
            v = X.plist[i]
            if in(v,dg)
                return false
            end
        end
    end

    # now check that ac doesn't intersect any of the sides of the polygon
    f(x) = alt_mod(x,n)
    exclude = [ f(k+j) for j=-2:1 ]
    for i=1:n
        if in(i,exclude)
            continue
        end
        S = X.plist[f(i)] +X.plist[f(i+1)]
        if meet_check(dg,S)
            return false
        end

    end

    # check that the diagonal is interior to X
    m = midpoint(dg)
    if !in(m,X)
        return false
    end

    return true
end

"""
`find_ear_diagonal(X::HPolygon)` returns an index `k` such that the
diagonal from `k-1` to `k+1` passes `check_diagonal`.
Returns `0` if no ear diagonal is found.
"""
function find_ear_diagonal(X::HPolygon)::Int
    for k=1:npoints(X)
        if check_diagonal(X,k)
            return k
        end
    end
    return 0
end

"""
`triangulate(X::HPolygon)` returns a list of triangles that
triangulate `X`. The polygon should have at least three sides, and
have no bad angles, and not self-intersect.
"""
function triangulate(X::HPolygon)::Array{HTriangle,1}
    n = npoints(X)

    if n < 3
        return HTriangle[]

    end

    if n == 3
        T = HTriangle(X)
        return [T]
    end

    k = find_ear_diagonal(X)

    if k==0
        @error "The polygon cannot be triangulated. Is it simple and proper?"
    end

    kprev = alt_mod(k-1,n)
    knext = alt_mod(k+1,n)
    a = X.plist[kprev]
    b = X.plist[k]
    c = X.plist[knext]
    T = a+b+c

    qlist = [ HPoint(X.plist[j]) for j=1:n if j!=k ]
    Y = HPolygon(qlist)
    result = triangulate(Y)
    pushfirst!(result, T)
    return result
end

"""
`area(X::HPolygon)` returns the area of the polygon, but is only
reliable if `X` does not self-intersect. See `is_simple`.
"""
function area(X::HPolygon)::Float64
    if npoints(X) < 3
        return 0.0
    end
    TT = triangulate(X)
    return sum(area(T) for T in TT)
end

# DEBUG STUFF

function draw_diagonal(X::HPolygon, k::Int)
    n = npoints(X)
    kprev = alt_mod(k-1,n)
    knext = alt_mod(k+1,n)
    a = X.plist[kprev]
    b = X.plist[knext]
    S = a+b
    set_line_style(S,:dot)
    draw(S)
end

function visualize(P::HPolygon)
    n = npoints(P)
    set_thickness(P,2)
    plot()
    draw(P)
    f(x) = alt_mod(x,n)
    for i=1:n
        a = P.plist[f(i-1)]
        ii = alt_mod(i+1,n)
        b = P.plist[ii]
        S = a+b
        if check_diagonal(P,i)
            set_color(S,:green)
        else
            continue
        end
        if in(midpoint(S),P)
            draw(S)
        end
    end

    for pt in P.plist
        draw(pt)
    end
    finish()
end


text:
using StreamReader
using FactCheck

facts("All") do
    @fact DEFAULT_PART_SIZE --> 1024
    @fact calculate(10) --> (1, 10, 10)
    @fact calculate(10, 3) --> (4, 3, 1)
    @fact calculate(1024) --> (1, 1024, 1024)
    @fact calculate(2048) --> (2, 1024, 1024)
    @fact calculate(2024) --> (2, 1024, 1000)

    pit = PartsIterator(10)
    @fact pit.size --> 10
    @fact pit.part_size --> 10
    @fact pit.left --> 0
    @fact pit.length --> 1
    @fact pit.part --> 0
    @fact pit.current_size --> 0
    @fact pit.loaded --> 0

    pit = PartsIterator(10, 3)

    @fact pit.size --> 10
    @fact pit.part_size --> 3
    @fact pit.left --> 0
    @fact pit.length --> 4
    @fact pit.part --> 0
    @fact pit.current_size --> 0
    @fact pit.loaded --> 0

    s = start(pit)
    @fact s --> 3
    @fact pit.part_size --> 3
    @fact pit.left --> 10
    @fact pit.length --> 4
    @fact pit.part --> 0
    @fact pit.current_size --> 0
    @fact pit.loaded --> 0

    @fact done(pit, s) --> false
    @fact pit.part_size --> 3
    @fact pit.left --> 10
    @fact pit.length --> 4
    @fact pit.part --> 0
    @fact pit.current_size --> 0
    @fact pit.loaded --> 0

    d, s = next(pit, s)
    @fact d --> 3
    @fact s --> 3
    @fact pit.part_size --> 3
    @fact pit.left --> 7
    @fact pit.length --> 4
    @fact pit.part --> 1
    @fact pit.current_size --> 3
    @fact pit.loaded --> 3

    @fact done(pit, s) --> false

    d, s = next(pit, s)
    @fact d --> 3
    @fact s --> 3
    @fact pit.part_size --> 3
    @fact pit.left --> 4
    @fact pit.length --> 4
    @fact pit.part --> 2
    @fact pit.current_size --> 3
    @fact pit.loaded --> 6

    @fact done(pit, s) --> false

    d, s = next(pit, s)
    @fact d --> 3
    @fact s --> 1
    @fact pit.part_size --> 3
    @fact pit.left --> 1
    @fact pit.length --> 4
    @fact pit.part --> 3
    @fact pit.current_size --> 3
    @fact pit.loaded --> 9

    @fact done(pit, s) --> false

    d, s = next(pit, s)
    @fact d --> 1
    @fact s --> 0
    @fact pit.part_size --> 3
    @fact pit.left --> 0
    @fact pit.length --> 4
    @fact pit.part --> 4
    @fact pit.current_size --> 1
    @fact pit.loaded --> 10

    @fact done(pit, s) --> true
    @fact [(1, 3, 3, 7, 30.0), (2, 3, 6, 4, 60.0), (3, 3, 9, 1, 90.0), (4, 1, 10, 0, 100.0)] --> [(pit.part, d, pit.loaded, pit.left, ceil(loadedpct(pit), 2)) for d in pit]

    #### READER ITERATOR TEST #####

    io = IOBuffer()
    s = write(io, "God is Love!!")
    seekstart(io)
    ir = IOStringReaderIterator(io, s, 4)
    edata = String["God ", "is L", "ove!", "!"]
    data = collect(ir)
    @fact edata --> data

    #### READER LIST ITERATOR TEST #####

    data = "God is Love!!"
    io = IOBuffer()
    s = write(io, data)
    seekstart(io)
    num_parts, parts_size, last_part_size = calculate(s, 4)
    parts_size_arr = Integer[parts_size for i in 1:num_parts-1]
    push!(parts_size_arr, last_part_size)

    edata = {{"Go", "d "}, {"is", " L"}, {"ov", "e!"}, {"!"}}

    rli = ReaderListIterator((mri) -> IOStringReaderIterator(io, parts_size_arr[mri.num], 2), num_parts)
    for (num, sri) in rli
        for d in sri
            @fact d --> edata[num][sri.num] "Mismatch edata[$num][$(sri.num)]"
        end
    end
end


text:
module Simplices

using LinearAlgebra

include("Vertices.jl")

function convergenceconfig(dim::Int, T::Type; kwargs...)
  kwargs = Dict(kwargs)
  xtol_abs = get(kwargs, :xtol_abs, zeros(real(T))) .* ones(Bool, dim)
  xtol_rel = get(kwargs, :xtol_rel, eps(real(T))) .* ones(Bool, dim)
  ftol_abs = get(kwargs, :ftol_abs, zero(real(T)))
  ftol_rel = get(kwargs, :ftol_rel, eps(real(T)))
  stopval = get(kwargs, :stopval, eps(real(T)))

  return (xtol_abs=xtol_abs, xtol_rel=xtol_rel,
          ftol_abs=ftol_abs, ftol_rel=ftol_rel,
          stopval=stopval)
end

struct Simplex{T<:Number, U<:Number}
  vertices::Vector{Vertex{T,U}}
  permabs::Vector{Int}
  function Simplex(vertices::Vector{Vertex{T,U}}) where {T<:Number, U<:Number}
    output = new{T,U}(vertices, zeros(Int64, length(vertices)))
    sort!(output)
    return output
  end
end

function Simplex(f::T, ic::AbstractVector{U}, step) where {T, U<:Number}
  return Simplex(f, vertexpositions(ic, step .* ones(Bool, length(ic))))
end

function Simplex(f::T, positions::U
    ) where {T, W<:Number, V<:AbstractVector{W}, U<:AbstractVector{V}}
  if length(unique(length.(positions))) != 1
    throw(ArgumentError("All entries in positions $positions must be the same
                        length"))
  end
  vertices = [Vertex(p, f(p)) for p in positions]
  return Simplex(vertices)
end

function partitionhypercube(f::T, a, b) where {T}
  dims = length(a)
  @assert dims == length(b)
  @assert all(a .< b)
  l = b .- a
  positions = partitionunitypositions(dims, Bool)
  return (Simplex(f, [p .* l .+ a for p in ps]) for ps in positions)
end

function partitionunitypositions(dims::Int, t::Type{T}=Bool) where {T}
  output = Vector{Vector{Vector{T}}}()
  @inbounds for i ∈ CartesianIndices(Tuple(ones(Int64, dims) .* dims))
    v = Tuple(i)
    any(j->any(v[j] .== v[1:j-1]), 1:dims) && continue
    inner = Vector{Vector{T}}()
    push!(inner, zeros(T, dims))
    for j ∈ 2:dims
      x = zeros(T, dims)
      for k ∈ 1:(j-1)
        x[v[k]] = one(T)
      end
      push!(inner, x)
    end
    push!(inner, ones(T, dims))
    @assert length(inner) == dims + 1
    push!(output, inner)
  end
  @assert length(output) == factorial(dims)
  return output
end
function sharesvertex(v::Vertex, s::Simplex)
  for w in s
    all(v.position .== w.position) && return true
  end
  return false
end
function Base.:∈(x, s::Simplex)
  return ∈(Vertex(x, one(value(first(s)))), s)
end
function Base.:∈(v::Vertex, s::Simplex)
  sharesvertex(v, s) && return true
  h = hash(s)
  function inner(i)
    swap!(s, i, v)
    x = hypervolume(s)
    swap!(s, v, i)
    return x
  end
  volume = mapreduce(inner, +, s.vertices)
  @assert h == hash(s) "The simplex has changed and it shouldn't have."
  return isapprox(volume, hypervolume(s), rtol=sqrt(eps()), atol=0)
end

function Base.:∈(a::Simplex, b::Simplex)
  isequal(a, b) && return true
  for v in a
    sharesvertex(v, b) && continue # doesn't count, could just be neighbours
    v ∈ b && return true
  end
  return false
end
import Base: length, iterate, push!, iterate, getindex
import Base: eachindex, sort!, hash, extrema
Base.length(s::Simplex) = length(s.vertices)
function Base.push!(s::Simplex, v::Vertex)
  push!(s.vertices, v)
  l = length(s.vertices)
  return nothing
end
Base.iterate(s::Simplex) = iterate(s.vertices)
Base.iterate(s::Simplex, counter) = iterate(s.vertices, counter)
Base.getindex(s::Simplex, index) = s.vertices[index]

function Base.sort!(s::Simplex)
  sort!(s.vertices, lt=vertexsorter)
  sortperm!(s.permabs, s.vertices, by=x->abs(value(x)))
  return nothing
end

function Base.extrema(s::Simplex)
  return [extrema(position(v)[i] for v in s) for i in 1:dimensionality(s)]
end

dimensionality(s::Simplex) = length(s) - 1

remove!(s::Simplex, v::Vertex) = filter!(x -> !isequal(x, v), s.vertices)
remove!(s::Simplex, x::Vector) = deleteat!(s.vertices, x)

issortedbyangle(s::Simplex) = issorted(s, by=v->angle(value(v)))

selectabs(s, index) = s.vertices[s.permabs[index]]

bestvertex(s::Simplex) = selectabs(s, 1)
worstvertex(s::Simplex) = selectabs(s, length(s))
secondworstvertex(s::Simplex) = selectabs(s, length(s) - 1)

function centroidposition(s::Simplex, ignoredvertex=worstvertex(s))
  g(v) = isequal(v, ignoredvertex) ? zero(position(v)) : position(v)
  return mapreduce(g, +, s) / (length(s) - 1)
end

centre(s::Simplex) = mapreduce(v->position(v), +, s) / length(s)

hypervolume(s::Simplex) = hypervolume(s.vertices)

function hypervolume(vs::AbstractVector{V}) where {V<:Vertex}
  m = hcat((vcat(position(v), 1) for v in vs)...)
  d = length(position(vs[1]))
  return abs(det(m)) / factorial(d)
end

function swap!(s::Simplex, this::Vertex, forthat::Vertex)
  before = deepcopy(s)
  @assert this ∈ s.vertices
  lengthbefore = length(s)
  remove!(s, this)
  @assert length(s) == lengthbefore - 1 "$(length(s)), $lengthbefore"
  push!(s, forthat)
  @assert forthat ∈ s.vertices
  sort!(s)
  @assert length(s) == lengthbefore
  return nothing
end

swapworst!(s::Simplex, forthis::Vertex) = swap!(s, worstvertex(s), forthis)

function closestomiddlevertex(s::Simplex)
  mid = mapreduce(position, +, s) ./ length(s)
  _, index = findmin(map(v->sum((position(v) - mid).^2), s))
  return s[index]
end

function assessconvergence(simplex, config)

  if abs(value(bestvertex(simplex))) <= config[:stopval]
    return :STOPVAL_REACHED
  end

  toprocess = Set{Int}(1)
  processed = Set{Int}()
  while !isempty(toprocess)
    vi = pop!(toprocess)
    v = simplex.vertices[vi]
    connectedto = Set{Int}()
    for (qi, q) ∈ enumerate(simplex)
      thisxtol = true
      for (i, (pv, pq)) ∈ enumerate(zip(position(v), position(q)))
        thisxtol &= isapprox(pv, pq, rtol=config[:xtol_rel][i],
                                     atol=config[:xtol_abs][i])
      end
      thisxtol && push!(connectedto, qi)
      thisxtol && for i in connectedto if i ∉ processed push!(toprocess, i) end end
    end
    push!(processed, vi)
  end
  allxtol = all(i ∈ processed for i ∈ 1:length(simplex))
  allxtol && return :XTOL_REACHED

  allftol = true
  for (vi, v) ∈ enumerate(simplex)
    for qi ∈ vi+1:length(simplex)
      q = simplex.vertices[qi]
      allftol &= all(isapprox(value(v), value(q),
                              rtol=config[:ftol_rel], atol=config[:ftol_abs]))
      all(position(v) .== position(q)) && return :XTOL_DEGENERATE_SIMPLEX
    end
  end
  allftol && return :FTOL_REACHED

  return :CONTINUE
end

function _πtoπ(ϕ)
  ϕ < -π && return _πtoπ(ϕ + 2π)
  ϕ > π && return _πtoπ(ϕ - 2π)
  return ϕ
end

windingangle(s::Simplex{T, <:Real}) where {T} = zero(T)
function windingangle(s::Simplex{T, <:Complex}) where {T}
  return sum(_πtoπ.(angle.(value.(circshift(s.vertices, -1))) .-
                    angle.(value.(s.vertices))))
end

windingnumber(s::Simplex{T, <:Real}) where {T} = 0
function windingnumber(s::Simplex{T, <:Complex}) where {T}
  radians = windingangle(s)
  return isfinite(radians) ? Int64(round(radians / 2π)) : Int64(0)
end

end


text:
@testset "SHIP Multi-Species" begin

##

@info("-------- TEST 🚢  Multi-Species-Basis ---------")
using SHIPs, JuLIP, BenchmarkTools, LinearAlgebra, Test, Random, StaticArrays
using SHIPs: eval_basis!, eval_basis, PolyCutoff1s, PolyCutoff2s
using JuLIP.MLIPs: IPSuperBasis
using JuLIP.Testing: print_tf
using JuLIP.Potentials: evaluate, evaluate_d

function randR()
   R = rand(JVecF) .- 0.5
   return (0.9 + 2 * rand()) * R/norm(R)
end

randR(N) = [ randR() for n=1:N ], zeros(Int16, N)

randR(N, syms::NTuple{M, Symbol}) where {M} =
      randR(N, Int16.(atomic_number.(syms)))

randR(N, Zs::NTuple{M, <:Integer}) where {M} =
      randR(N, Int16.(Zs))

randR(N, Zs::NTuple{M, Int16}) where {M}  =
      randR(N)[1], rand(Zs, N), rand(Zs)

function randiso()
   K = @SMatrix rand(3,3)
   K = K - K'
   Q = rand([-1,1]) * exp(K)
end

function randiso(Rs, Zs)
   Q = randiso()
   p = shuffle(1:length(Rs))
   return [ Q * R for  R in Rs[p]], Zs[p]
end

##
@info("Testing creation and (de-)dictionisation of a few BasisSpecs")

spec = SparseSHIP(:X, 3, 10)
println(@test spec == SparseSHIP(3, 10))
println(@test decode_dict(Dict(spec)) == spec)

spec = SparseSHIP([:H, :He], 3, 6)
println(@test spec == SparseSHIP([1,2], 3, 6))
println(@test decode_dict(Dict(spec)) == spec)

spec = SparseSHIP([:H, :He, :Li], 2, 6)
println(@test spec == SparseSHIP([1,2,3], 2, 6))
println(@test decode_dict(Dict(spec)) == spec)

##

trans = PolyTransform(2, 1.0)
cutf = PolyCutoff2s(2, 0.5, 3.0)

ship2 = SHIPBasis(SparseSHIP([1,2], 2, 10, wL=2.0), trans, cutf)
ship3 = SHIPBasis(SparseSHIP([1,2], 3,  8, wL=2.0), trans, cutf)
ship4 = SHIPBasis(SparseSHIP([1,2], 4,  6, wL=1.5), trans, cutf)
ship5 = SHIPBasis(SparseSHIP([1,2], 5,  5, wL=1.5), trans, cutf)
ships = [ship2, ship3, ship4, ship5]

##

@info("Test (de-)dictionisation of basis sets")
for ship in ships
   println(@test (decode_dict(Dict(ship)) == ship))
end

Rs, Zs, iz = randR(20, (1,2))
🚢 = ship3
eval_basis(🚢, Rs, Zs, iz)

@info("Test isometry invariance for 3B-6B 🚢 s")
for ntest = 1:20
   Rs, Zs, iz = randR(20, (1,2))
   BB = [ eval_basis(🚢, Rs, Zs, iz) for 🚢 in ships ]
   RsX, ZsX = randiso(Rs, Zs)
   BBX = [ eval_basis(🚢, RsX, ZsX, iz) for 🚢 in ships ]
   for (B, BX) in zip(BB, BBX)
      print_tf(@test B ≈ BX)
   end
end
println()

##

@info("Test gradients for 3-6B 🚢-basis")
for 🚢 in ships
   @info("  body-order = $(SHIPs.bodyorder(🚢)):")
   Rs, Zs, z = randR(20, (1,2))
   tmp = SHIPs.alloc_temp_d(🚢, Rs)
   # SHIPs.precompute_dA!(tmp, 🚢, Rs, Zs)
   B = eval_basis(🚢, Rs, Zs, z)
   dB = SHIPs.alloc_dB(🚢, Rs)
   SHIPs.eval_basis_d!(dB, tmp, 🚢, Rs, Zs, z)
   @info("      finite-difference test into random directions")
   for ndirections = 1:20
      Us, _ = randR(length(Rs))
      errs = Float64[]
      for p = 2:10
         h = 0.1^p
         Bh = eval_basis(🚢, Rs+h*Us, Zs, z)
         dBh = (Bh - B) / h
         dBxU = sum( dot.(Ref(Us[n]), dB[n,:])  for n = 1:length(Rs) )
         push!(errs, norm(dBh - dBxU, Inf))
      end
      success = (/(extrema(errs)...) < 1e-3) || (minimum(errs) < 1e-10)
      print_tf(@test success)
   end
   println()
end

##

@info("Check Correctness of Multi-species SHIPBasis calculators")

randcoeffs(B) = 2 * (rand(length(B)) .- 0.5) .* (1:length(B)).^(-2)

m_naive_energy(basis::SHIPBasis, at) =
      sum( eval_basis(basis, R, at.Z[j], at.Z[i])
            for (i, j, R) in sites(at, cutoff(basis)) )

for basis in ships
   @info("   body-order = $(SHIPs.bodyorder(basis))")
   at = bulk(:Si) * 3
   at.Z[:] .= 1
   at.Z[2:2:end] .= 2
   rattle!(at, 0.1)
   print("     energy: ")
   println(@test energy(basis, at) ≈ m_naive_energy(basis, at) )
   print("site-energy: ")
   println(@test energy(basis, at) ≈ sum( site_energy(basis, at, n)
                                         for n = 1:length(at) ) )
   # we can test consistency of forces, site energy etc by taking
   # random inner products with coefficients
   @info("     a few random combinations")
   for n = 1:10
      c = randcoeffs(basis)
      sh = JuLIP.MLIPs.combine(basis, c)
      print_tf(@test energy(sh, at) ≈ dot(c, energy(basis, at)))
      print_tf(@test forces(sh, at) ≈ sum(c*f for (c, f) in zip(c, forces(basis, at))) )
      print_tf(@test site_energy(sh, at, 5) ≈ dot(c, site_energy(basis, at, 5)))
      print_tf(@test site_energy_d(sh, at, 5) ≈ sum(c*f for (c, f) in zip(c, site_energy_d(basis, at, 5))) )
   end
   println()
end

# ----------------------------------------------------------------------
#  Calculator Tests (everything above is a Basis test)
# ----------------------------------------------------------------------

shipsB = ships

@info("--------------- Fast Multi-🚢 Implementation ---------------")

@info("Testing correctness of `SHIP` against `SHIPBasis`")
for B in shipsB
   @info("   bodyorder = $(SHIPs.bodyorder(B))")
   coeffs = randcoeffs(B)
   ship = SHIP(B, coeffs)
   @info("   test (de-)dictionisation")
   println(@test decode_dict(Dict(ship)) == ship)
   @show length(B), length(ship)
   tmp = SHIPs.alloc_temp(ship, 0)
   @info("      check that SHIPBasis ≈ SHIP")
   for ntest = 1:30
      Rs, Zs, z0 = randR(10, (1,2))
      Es = SHIPs.evaluate!(tmp, ship, Rs, Zs, z0)
      Bs = dot(coeffs, SHIPs.eval_basis(B, Rs, Zs, z0))
      print_tf(@test Es ≈ Bs)
   end
   println()
   # ------------------------------------------------------------
   # @info("      Quick timing test")
   # Nr = 30
   # Rs, Zs, z0 = randR(Nr, (1,2))
   # b = SHIPs.alloc_B(B)
   # tmp = SHIPs.alloc_temp(ship, Nr)
   # tmpB = SHIPs.alloc_temp(B, Nr)
   # print("       SHIPBasis : "); @btime SHIPs.eval_basis!($b, $tmpB, $B, $Rs, $Zs, $z0)
   # print("            SHIP : "); @btime SHIPs.evaluate!($tmp, $ship, $Rs, $Zs, $z0)
   # println()
end

##
@info("Check Correctness of SHIP gradients")
for B in shipsB
   @info("   body-order = $(SHIPs.bodyorder(B))")
   coeffs = randcoeffs(B)
   ship = SHIP(B, coeffs)
   Rs, Zs, z0 = randR(10, (1,2))
   tmp = SHIPs.alloc_temp_d(ship, length(Rs))
   dEs = zeros(JVecF, length(Rs))
   SHIPs.evaluate_d!(dEs, tmp, ship, Rs, Zs, z0)
   Es = SHIPs.evaluate!(tmp, ship, Rs, Zs, z0)
   println(@test Es ≈ evaluate(ship, Rs, Zs, z0))
   println(@test dEs ≈ evaluate_d(ship, Rs, Zs, z0))
   @info("      Correctness of directional derivatives")
   for ndir = 1:20
      U = [rand(JVecF) .- 0.5 for _=1:length(Rs)]
      errs = Float64[]
      for p = 2:10
         h = 0.1^p
         dEs_U = dot(dEs, U)
         dEs_h = (SHIPs.evaluate!(tmp, ship, Rs + h * U, Zs, z0) - Es) / h
         push!(errs, abs(dEs_h - dEs_U))
      end
      success = (/(extrema(errs)...) < 1e-3) || (minimum(errs) < 1e-10)
      print_tf(@test success)
   end
   println()
end

##
@info("Check Correctness of SHIP calculators")

m_naive_energy(ship::SHIP, at) =
      sum( evaluate(ship, R, at.Z[j], at.Z[i])
            for (i, j, R) in sites(at, cutoff(ship)) )

for B in shipsB
   @info("   body-order = $(SHIPs.bodyorder(B))")
   coeffs = randcoeffs(B)
   ship = SHIP(B, coeffs)
   at = bulk(:Si) * 3
   at.Z[1:2:end] .= 1
   at.Z[2:2:end] .= 2
   rattle!(at, 0.1)
   print("     energy: ")
   println(@test energy(ship, at) ≈ m_naive_energy(ship, at) )
   print("site-energy: ")
   println(@test energy(ship, at) ≈ sum( site_energy(ship, at, n)
                                         for n = 1:length(at) ) )
   println("forces: ")
   println(@test JuLIP.Testing.fdtest(ship, at))
   println("site-forces: ")
   println(@test JuLIP.Testing.fdtest( x -> site_energy(ship, set_dofs!(at, x), 3),
                                       x -> mat(site_energy_d(ship, set_dofs!(at, x), 3))[:],
                                       dofs(at) ) )
end

##

@info("Multi-species filtering test")
trans = PolyTransform(2, 1.0)
cutf = PolyCutoff2s(2, 0.5, 3.0)
ship = SHIPBasis(SparseSHIP((1, 2), 5,  6; wL = 1.0), trans, cutf, filter=false)
@show maxgrp = maximum(SHIPs.alllen_bgrp(ship, 1))
@time fship = SHIPs.alg_filter_rpi_basis(ship)
@show length(fship), length(ship)
println(@test length(fship) < length(ship))

end

max_stars_repo_name: cortner/SHIPs

max_issues_repo_name: cortner/SHIPs

max_forks_repo_name: cortner/SHIPs


text:
push!(LOAD_PATH, "./src/")

using DelimitedFiles
using JLD2

using QuantumModels
using MatrixProductStateTools
using ExactDiagonalizationTools
using SymTensors
using GaussianFermions
using GutzwillerMPS
using TensorNetAlgs

ly = 6
lx = 30
L = lx * ly
t1, t2, t3, mu = 1.0, 1.0, 1.0, 0.0
j, k = 2.0, 0.6
boundary = (:PBC, :OBC)
b_str = "PBCY"

hamil = generatebdg(
    triangularhopping((ly, lx), t1, t2, t3, mu, boundary=boundary)
)
cm = correlationmatrix(hamil, div(L, 2))
fgs = corrmat2gmps(cm)

sz, sp, sm = spinoperators(1/2, symmetry=U1)

infos = "$b_str,ly=$ly,lx=$lx,t1=$t1,t2=$t2,t3=$t3,mu=$mu"
for m in [2000]

    print("Fermionic mps m=$m ")
    filename = "mps,triangular,hop,$infos,m=$m.jld2"
    if isfile(filename)
        println("already exists! loading...")
        @load filename mps
    else
        println("making ...")
        mps = gmps2mps(fgs, m, symmetry=U1)
        #@save filename mps

        # zzdata = measure(mps, sz, sz)
        # pmdata = measure(mps, sp, sm)
        vnedata = entanglemententropy(mps)

        println("Measurements for hop m = $m ...")
        # writedlm( "zz,triangular,hop,$infos,m=$m.dat", zzdata)
        # writedlm( "pm,triangular,hop,$infos,m=$m.dat", pmdata)
        writedlm( "vne,triangular,hop,$infos,m=$m.dat", vnedata)
    end

    for M in [12000]
        println("ZipandGutzwiller for M = $M ...")

        #filename = "mps,triangular,gutzF23,$infos,m=$m,M=$M.jld2"
        if isfile(filename)
            println("Already exists! dismiss.")
        else
            println("Making ...")
            mpsgutz = zipandgutzwiller!(mps, mps, mode=:F23, maxdim=M)
            normalize!(mpsgutz)
            #@save  filename mpsgutz
            println("Measurements for M = $M ...")
            # zzdata = measure(mpsgutz, sz, sz)
            # pmdata = measure(mpsgutz, sp, sm)
            vnedata = entanglemententropy(mpsgutz)
            # writedlm("zz,triangular,gutzF23,$infos,m=$m,M=$M.dat", zzdata)
            # writedlm("pm,triangular,gutzF23,$infos,m=$m,M=$M.dat", pmdata)
            writedlm("vne,triangular,gutzF23,$infos,m=$m,M=$M.dat", vnedata)

            mpo = generatempo(triangularspinmodel((ly, lx), j, j, j, k, k, k,
                                                  boundary=boundary, symmetry=U1))
            reducempo!(mpo)
            center_at!(mpsgutz, 1)
            env = initialenv(mps, mpo)
            dmrg2sitesweep!(mps, mpo, env, maxdim=M, verbose=true)

            dmrginfo = "j=$j,k=$k,M=$M"
            # zzdata = measure(mpsgutz, sz, sz)
            # pmdata = measure(mpsgutz, sp, sm)
            vnedata = entanglemententropy(mpsgutz)
            # writedlm("zz,triangular,1dmrg,$dmrginfo,gutzF23,$infos,m=$m,M=$M.dat", zzdata)
            # writedlm("pm,triangular,1dmrg,$dmrginfo,gutzF23,$infos,m=$m,M=$M.dat", pmdata)
            writedlm("vne,triangular,1dmrg,$dmrginfo,gutzF23,$infos,m=$m,M=$M.dat", vnedata)

        end
    end
end

max_stars_repo_name: amiragha/mps

max_issues_repo_name: amiragha/mps

max_forks_repo_name: amiragha/mps


text:
"""
    read_from_gfa(filename::String; weight_file::String)

Read graph from GFA file and optional weight_file (contains two column `node` and `weight`) and return `GFAResult` struct which
has `g` - the graph, `w` - weight array, `l` - node label array, `e` - edge label array and `p` - path array.
"""
function read_from_gfa(filename::String; weight_file::String="")
    u = readlines(filename);
    l_dict = Dict()
    w_dict = Dict()
    l_w_dict = Dict()
    l_s_dict = Dict()
    l_i_dict = Dict()
    links = Tuple[]
    e_o_dict = Dict()
    e_i_dict = Dict()
    p_dict = Dict()
    p::Array{Path} = Path[]

    if weight_file != ""
        w_f_dict = Dict()
        f = CSV.File(weight_file)
        for row in f
            w_f_dict[string(row.node)] = row.weight
        end
    end

    i = 1
    for line in u
        if line[1] == 'S'
            s_dummy = split(line, "\t")
            l_dict[i] = s_dummy[2] * "\t+"
            if weight_file != ""
                w_dict[i] = w_f_dict[s_dummy[2]]
            else
                w_dict[i] = length(s_dummy[3])
            end
            l_w_dict[l_dict[i]] = w_dict[i]
            l_w_dict[s_dummy[2] * "\t-"] = 0
            l_s_dict[l_dict[i]] = s_dummy[3]
            l_i_dict[l_dict[i]] = SubString(line, length(s_dummy[1] * "\t" * s_dummy[2] * "\t" * s_dummy[3] * "\t") + 1)
            i += 1
        elseif line[1] == 'L'
            l_dummy = split(line, "\t")
            k_start = l_dummy[2] * "\t" * l_dummy[3]
            if l_dummy[3] == "-"
                if l_w_dict[l_dummy[2] * "\t" * l_dummy[3]] == 0
                    l_dict[i] = l_dummy[2] * "\t-"
                    w_dict[i] = l_w_dict[l_dummy[2] * "\t+"]
                    l_w_dict[l_dict[i]] = w_dict[i]
                    l_s_dict[l_dict[i]] = l_s_dict[l_dummy[2] * "\t+"]
                    l_i_dict[l_dict[i]] = l_i_dict[l_dummy[2] * "\t+"]
                    i += 1
                end
            end
            k_end = l_dummy[4] * "\t" * l_dummy[5]
            if l_dummy[5] == "-"
                if l_w_dict[l_dummy[4] * "\t" * l_dummy[5]] == 0
                    l_dict[i] = l_dummy[4] * "\t-"
                    w_dict[i] = l_w_dict[l_dummy[4] * "\t+"]
                    l_w_dict[l_dict[i]] = w_dict[i]
                    l_s_dict[l_dict[i]] = l_s_dict[l_dummy[4] * "\t+"]
                    l_i_dict[l_dict[i]] = l_i_dict[l_dummy[4] * "\t+"]
                    i += 1
                end
            end
            e_name = k_start * "::" * k_end
            e_o_dict[e_name] = l_dummy[6]
            e_i_dict[e_name] = SubString(line, length(l_dummy[1] * "\t" * l_dummy[2] * "\t" * l_dummy[3] * "\t" * l_dummy[4] * "\t" * l_dummy[5] * "\t" * l_dummy[6] * "\t") + 1)
            push!(links, (k_start, k_end))
        elseif line[1] == 'P'
            p_dummy = split(line, "\t")
            path_name = p_dummy[2]
            p_extract = split(p_dummy[3], ",")
            p_arr = []
            for i = 1:(length(p_extract))
                last_char = p_extract[i][length(p_extract[i])]
                extracted = chop(p_extract[i])
                extracted = extracted * "\t" * last_char
                push!(p_arr, extracted)
            end
            p_dict[path_name] = p_arr
        end
    end

    w = _weight_from_dict(w_dict)
    l = _label_from_dict(l_dict, l_s_dict, l_i_dict)
    e = _edge_label_from_dict(e_o_dict, e_i_dict)

    g = SimpleDiGraph(length(w))

    reverse_l_dict = Dict(value => key for (key, value) in l_dict)

    for link in links
        add_edge!(g, reverse_l_dict[link[1]], reverse_l_dict[link[2]])
    end

    for (name, p_arr) in p_dict
        path = []
        for i = 1:length(p_arr)
            push!(path, reverse_l_dict[p_arr[i]])
        end
        push!(p, Path(name, path))
    end
    if weight_file != ""
        return GFAResult(g, w, l, e, p, false)
    else
        return GFAResult(g, w, l, e, p, false)
    end
end


text:
#######################################################################
function infertemperature(randpcaline, Spca, Sgal, maxiter)
#######################################################################

    #----------------------------------------------------
    # Sample PCA parameters that describe total flux
    #----------------------------------------------------

    param = [randpcaline() for s=1:Spca]

    W = [p[1] for p in param]

    μ = [p[2] for p in param]

    #----------------------------------------------------
    # Define function for VI
    #----------------------------------------------------

    function logp(T, scale)

        @assert(T >= mintemperature())

        @assert(T <= maxtemperature())

        local ℓ = zeros(eltype(T), Spca)

        local σ = 0.2

        local y = vec(scale * observedunitgalaxyvector(T))

        for s in 1:Spca

            local M = W[s]'*W[s] # + Matrix(Diagonal(σ.*σ))

            local proj(x) = vec(M \ (W[s]'*(x - μ[s])))[1]

            local fwd(z)  = vec(W[s]*z + μ[s])

            ℓ[s]    = logpdf(MvNormal(fwd(proj(y)), σ), y)

        end

        return logsumexp(ℓ)

    end

    function auxiliary(param)

        local T = converttemperature(param[1])

        local scale = convertscale(param[2])

        logp(T, scale)

    end

    #----------------------------------------------------
    # Initial mean for VI
    #----------------------------------------------------

    T0     = convertunctemperature(mintemperature() + rand()*(maxtemperature()-mintemperature()))

    scale0 = log(rand()*100.0)

    #--------------------------------------------------------
    # Call optimiser
    #--------------------------------------------------------

    opt    = Optim.Options(iterations = 10000, show_trace = true, show_every=2)
    result = Optim.optimize(x->-auxiliary(x), [T0;scale0], NelderMead(), opt)

    paramOpt = result.minimizer
    TOpt     = paramOpt[1]
    scaleOpt = paramOpt[2]

    @show converttemperature(TOpt)

    # #----------------------------------------------------
    # # Run VI
    # #----------------------------------------------------
    #
    # mu, Sigma, neglogev = VI(auxiliary, [TOpt; scaleOpt], S = Sgal, maxiter = maxiter, optimiser=NelderMead())
    #
    #
    # #----------------------------------------------------
    # # Output sampling function
    # #----------------------------------------------------
    #
    # jointPDFunc = MvNormal(mu, Sigma)
    #
    # function randobservedgalaxyvector()
    #
    #     local p = rand(jointPDFunc)
    #
    #     converttemperature(p[1]), convertscale(p[2])
    #
    # end
    #
    # return randobservedgalaxyvector

end


text:
using Plots
#komvos1 komvos2 komvos3
xyz = 0.1*[0 1 -1;
           0 2 -1;
           0 0  0
]; #orismos komvon san xyz
n = size(xyz)[2]; #arithmos komvon

#komvos1_start komvos1_end
numel = [1 2;
         1 3
]; #orismos ton stixion os grammes me toys akraious komvous

#orismos elastikotitas stixion
E = 220e9.*ones(N,1); A = 4e-4.*ones(N,1);

N = size(numel)[1]; #arithmos stoixion

#orismos ton desmeumenon komvon
BC = zeros(3*n,1)
#for 2d uncomment below line
dimconstr = 1:3:3*n; BC[dimconstr] .= 1;
BC[4:end] .= 1; #0==free
fixedofs = findall(x->x==1,BC[:])
freedofs = findall(x->x==0,BC[:])
#boundary conditions
U = zeros(3*n,1); U[4:end].=0;
U[dimconstr].=0; #uncomment for 2d
F = zeros(3*n,1); F[[1 2]]=[-100 100] ;
#orismos olikou mitroou stivarotitas
global K = zeros(3*n,3*n);
#start main loop to calculate stiffness matrix
for i=1:N
    i1 = numel[i,1]
    i2 = numel[i,2]
    #calclulate directions
    x1 = xyz[1,i1]; x2 = xyz[1,i2]; dx = x2-x1;
    y1 = xyz[2,i1]; y2 = xyz[2,i2]; dy = y2-y1;
    z1 = xyz[3,i1]; z2 = xyz[3,i2]; dz = z2-z1;
    L = sqrt(dx^2+dy^2+dz^2);
    nx = dx/L; ny = dy/L; nz = dz/L;
    #direction matrix
    I = [nx^2  nx*ny nx*nz;
         ny*nx ny^2  ny*nz;
         nz*nx nz*ny nz^2];
    k = E[i]*A[i]/L
    Ke = k*[I -I; -I I]
    #println(size(Ke),"\n")
    edofs = [3*i1-2 3*i1-1 3*i1  3*i2-2 3*i2-1 3*i2]
    K[edofs,edofs] = K[edofs,edofs][1,:,1,:] + Ke;
    #println(K,"\n")
end
Kfinal = K[freedofs,freedofs]
rhs = F
rhs = rhs[freedofs]
rhs = rhs - K[freedofs,fixedofs]*U[fixedofs]

Ufree = Kfinal\rhs
println(Ufree)

max_stars_repo_path: axnai julia/diktioma-efelkismos-3d.jl

max_stars_repo_name: billyon/first

max_issues_repo_path: axnai julia/diktioma-efelkismos-3d.jl

max_issues_repo_name: billyon/first

max_forks_repo_path: axnai julia/diktioma-efelkismos-3d.jl

max_forks_repo_name: billyon/first


text:
"""
    TransformedOutput(f, init; tspan::AbstractRange, kw...)

An output that stores the result of some function `f` of the grid/s

# Arguments

- `f`: a function or functor that accepts an `AbstractArray` or `NamedTuple` of
    `AbstractArray` with names matching `init`. The `AbstractArray` will be a view into
    the grid the same size as the init grids, removing any padding that has been added.
- `init`: initialisation `Array` or `NamedTuple` of `Array`

# Keywords

- `tspan`: `AbstractRange` timespan for the simulation
- `aux`: NamedTuple of arbitrary input data. Use `get(data, Aux(:key), I...)`
    to access from a `Rule` in a type-stable way.
- `mask`: `BitArray` for defining cells that will/will not be run.
- `padval`: padding value for grids with neighborhood rules. The default is `zero(eltype(init))`.

$EXPERIMENTAL
"""
mutable struct TransformedOutput{T,A<:AbstractVector{T},E,F,B} <: Output{T,A}
    frames::A
    running::Bool
    extent::E
    f::F
    buffer::B
end
function TransformedOutput(f::Function, init::Union{NamedTuple,AbstractMatrix}; extent=nothing, kw...)
    # We have to handle some things manually as we are changing the standard output frames
    extent = extent isa Nothing ? Extent(; init=init, kw...) : extent
    # Define buffers to copy to before applying `f`
    buffer = init isa NamedTuple ? map(zero, init) : zero(init)
    zeroframe = f(buffer)
    # Build simulation frames from the output of `f` for empty frames
    frames = [deepcopy(zeroframe) for f in eachindex(tspan(extent))]
    # Set the first frame to the output of `f` for `init`
    frames[1] = f(init)

    return TransformedOutput(frames, false, extent, f, buffer)
end
function TransformedOutput(init; kw...)
    throw(ArgumentError("TransformedOutput must be passed a function and the init grid(s) as arguments"))
end

function storeframe!(o::TransformedOutput, data::AbstractSimData)
    transformed = _transform_grids(o, grids(data))
    i = frameindex(o, data)
    # Copy the transformed grid/s to the output frames,
    # instead of just assigning (see issue #169)
    o[i] = _copytransformed!(o[i], transformed)
end

# Copy arrays manually as reducing functions can return the original object without copy.
_copytransformed!(dest::NamedTuple, src::NamedTuple) = map(_copytransformed!, dest, src)
_copytransformed!(dest::AbstractArray, src::AbstractArray) = dest .= src
# Non-array output is just assigned
_copytransformed!(dest, src) = src

# Multi/named grid simulation, f is passed a NamedTuple
function _transform_grids(o::TransformedOutput, grids::NamedTuple)
    # Make a new named tuple of raw arrays without wrappers, copying
    # to the buffer where an OffsetArray was used for padding.
    # Often it's faster to copy than use a view when f is sum/mean etc.
    nt = map(grids, o.buffer) do g, b
        source(g) isa OffsetArray ? copy!(b, sourceview(g)) : source(g)
    end
    o.f(nt)
end
# Single unnamed grid simulation, f is passed an AbstractArray
function _transform_grids(o::TransformedOutput, grids::NamedTuple{(DEFAULT_KEY,)})
    g = first(grids)
    A = source(g) isa OffsetArray ? copy!(o.buffer, sourceview(g)) : source(g)
    o.f(A)
end

init_output_grids!(o::TransformedOutput, init) = nothing
_initdata!(o::TransformedOutput, init) = nothing


text:
module MCMC

Base.warn("MCMC.jl has moved to Lora.jl. Development in Lora.jl will continue. MCMC.jl is a placeholder for the future merge of various independent MCMC implementations in Julia, including Lora.jl.")

end


text:
# By default, Julia/LLVM does not use fused multiply-add operations (FMAs).
# Since these FMAs can increase the performance of many numerical algorithms,
# we need to opt-in explicitly.
# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details.
@muladd begin

"""
    UnstructuredQuadSortedBoundaryTypes{N, BCs<:NTuple{N, Any}}

General container to sort the boundary conditions by type for the unstructured quadrilateral solver.
It stores a set of global indices for each boundary condition type to expedite computation
during the call to `calc_boundary_flux!`. The original dictionary form of the boundary conditions
set by the user in the elixir file is also stored for printing.

!!! warning "Experimental code"
    This boundary condition container is experimental and can change any time.
"""
mutable struct UnstructuredQuadSortedBoundaryTypes{N, BCs<:NTuple{N, Any}}
  boundary_condition_types::BCs # specific boundary condition type(s), e.g. BoundaryConditionWall
  boundary_indices::NTuple{N, Vector{Int}} # integer vectors containing global boundary indices
  boundary_dictionary::Dict{Symbol, Any} # boundary conditions as set by the user in the elixir file
end

# constructor that "eats" the original boundary condition dictionary and sorts the information
# from the `UnstructuredBoundaryContainer2D` in cache.boundaries according to the boundary types
# and stores the associated global boundary indexing in NTuple
function UnstructuredQuadSortedBoundaryTypes(boundary_conditions::Dict, cache)
  # extract the unique boundary function routines from the dictionary
  boundary_condition_types = Tuple(unique(collect(values(boundary_conditions))))
  n_boundary_types = length(boundary_condition_types)
  boundary_indices = ntuple(_ -> [], n_boundary_types)

  container = UnstructuredQuadSortedBoundaryTypes{n_boundary_types, typeof(boundary_condition_types)}(
    boundary_condition_types, boundary_indices, boundary_conditions)

  initialize!(container, cache)
end

function initialize!(boundary_types_container::UnstructuredQuadSortedBoundaryTypes{N}, cache) where N
  @unpack boundary_dictionary, boundary_condition_types = boundary_types_container

  unique_names = unique(cache.boundaries.name)

  # Verify that each Dict key is a valid boundary name
  for key in keys(boundary_dictionary)
    if !(key in unique_names)
      error("Key $(repr(key)) is not a valid boundary name")
    end
  end

  # Verify that each boundary has a boundary condition
  for name in unique_names
    if name !== Symbol("---") && !haskey(boundary_dictionary, name)
      error("No boundary condition specified for boundary $(repr(name))")
    end
  end

  # pull and sort the indexing for each boundary type
  _boundary_indices = Vector{Any}(nothing, N)
  for j in 1:N
    indices_for_current_type = Int[]
    for (test_name, test_condition) in boundary_dictionary
      temp_indices = findall(x->x===test_name, cache.boundaries.name)
      if test_condition === boundary_condition_types[j]
        indices_for_current_type = vcat(indices_for_current_type, temp_indices)
      end
    end
    _boundary_indices[j] = sort!(indices_for_current_type)
  end

  # convert the work array with the boundary indices into a tuple
  boundary_types_container.boundary_indices = Tuple(_boundary_indices)

  return boundary_types_container
end

end # @muladd


text:
# Rohan McLure, Australian National University (2019)
# Distributed implementation of the Reduce Parallel Research Kernel
@everywhere using ArrayChannels

function main()
    local iterations, payload
    if length(ARGS) != 2 || nprocs() == 1
        println("Usage: reduce -p n <# iterations> <vector_length>")
        exit(1)
    end

    argv = map(x -> parse(Int, x), ARGS)
    iterations, payload = argv

    ranks = nworkers()

    A = ArrayChannel(Float64, workers(), payload)

    futures = Vector{Future}(undef, ranks)
    @sync for proc in workers()
        futures[proc-1] = @spawnat proc work(A, 2, iterations, payload)
    end

    # Fetch results
    results = map(fetch, futures)
    value = results[1][1]
    time = maximum(map(x -> x[2], results))

    ground_truth = iterations+2.0+(iterations*iterations+5.0*iterations+4.0)*(ranks-1)/2;

    if abs(value - ground_truth) <= 1.e-8
        println("Solution validates")
    else
        println("Value provided: $value, expected: $ground_truth")
    end
    avgtime = time / iterations
    throughput = 1.e-6 * (2.0*ranks-1.0)*payload/avgtime
    println("Rate (MFlops/s): $throughput Avg time (s): $avgtime")
end

@everywhere mean(l) = sum(l) / length(l)

@everywhere function work(A, root, iterations, payload)
    local t0, t1
    fill!(A, 1.0)
    A_data = A.buffer.src
    constant_vector = ones(payload)
    for k in 0 : iterations
        if k == 1
            t0 = time_ns()
        end

        A_data .+= constant_vector
        reduce!(+, A, root)
    end

    t1 = time_ns()
    time = (t1-t0)*1.e-9

    if myid() == root
        return (A[1], time)
    else
        return (nothing, time)
    end
end

main()


text:
import RegressionAndOtherStories: model_summary

function model_summary(model::SampleModel, params::Vector{Symbol};
    round_estimates=true)

    sdf = read_summary(model)

    stats=names(sdf)
    stats[8] = "n_eff"
    items = [2,4,5,6,7,8,10]

    # Remove unwanted rows and create parameter pairs
    count = 1
    df = DataFrame()
    parameters = Pair{Symbol, Int}[]
    for (index, par) in enumerate(sdf[:, :parameters])
        if par in params
            push!(df, sdf[index, items])
            append!(parameters, [par => count])
            count += 1
        end
    end

    statistics = Pair{String, Int}[]
    for (index, stat) in enumerate(stats[items])
        append!(statistics, [stat => index])
    end

    if round_estimates
        estimates = round.(Array(df); digits=3)
    else
        estimates = Array(df)
    end

    return NamedArray(
        estimates,
        (OrderedDict(parameters...), OrderedDict(statistics...)),
        ("Par", "Stat")
    )
end

function model_summary(model::SampleModel, params::Vector{String};
    round_estimates=true)

    sdf = read_summary(model)

    stats=names(sdf)
    stats[8] = "n_eff"
    items = [2,4,5,6,7,8,10]

    # Remove unwanted rows and create parameter pairs
    count = 1
    df = DataFrame()
    parameters = Pair{String, Int}[]
    for (index, par) in enumerate(String.(sdf[:, :parameters]))
        if par in params
            push!(df, sdf[index, items])
            append!(parameters, [par => count])
            count += 1
        end
    end

    statistics = Pair{String, Int}[]
    for (index, stat) in enumerate(stats[items])
        append!(statistics, [stat => index])
    end

    if round_estimates
        estimates = round.(Array(df); digits=3)
    else
        estimates = Array(df)
    end

    return NamedArray(
        estimates,
        (OrderedDict(parameters...), OrderedDict(statistics...)),
        ("Par", "Stat")
    )

end


text:
using DrWatson
@quickactivate "PhaseBases"
using PhaseBases
using Test
using LinearAlgebra
using RecursiveArrayTools

@testset "Zernike generation" begin
    @test zernike(.5, .5, 3) == (z = [1.0, 0.5, 0.5, 0.5, 0.0, 0.0, 0.25, -0.25, -0.25, -0.25], zx = [0.0, 0.0, 1.0, 1.0, 2.0, 1.0, 1.5, 1.5, 1.0, 0.0], zy = [0.0, 1.0, 0.0, 1.0, 2.0, -1.0, 0.0, 1.0, 1.5, -1.5])
    z = ZernikeBW(5, 4);
    @test z.elements[:,:,15] ==
   [-4.0     -0.4375  1.0     -0.4375  -4.0
    -0.4375  -0.25    0.0625  -0.25    -0.4375
     1.0      0.0625  0.0      0.0625   1.0
    -0.4375  -0.25    0.0625  -0.25    -0.4375
    -4.0     -0.4375  1.0     -0.4375  -4.0]
    @test PhaseBases.elements(z, [5,10])[1] ==
        [3.0 1.5 1.0 1.5 3.0; 1.5 0.0 -0.5 0.0 1.5; 1.0 -0.5 -1.0 -0.5 1.0; 1.5 0.0 -0.5 0.0 1.5; 3.0 1.5 1.0 1.5 3.0]
    @test  PhaseBases.aperturedelements(z, 1) == PhaseBases.aperturedelements(z, [1,2])[:,:,1] == [
     0.0  0.0  1.0  0.0  0.0
     0.0  1.0  1.0  1.0  0.0
     1.0  1.0  1.0  1.0  1.0
     0.0  1.0  1.0  1.0  0.0
     0.0  0.0  1.0  0.0  0.0]

     coef = zeros(length(z))
     coef[10:11] = [-1,1]
     ph = compose(z, coef)
     @test ph == [
        -2.0   -2.875  0.0   2.875   2.0
        1.75  -0.25   0.0   0.25   -1.75
        1.0    0.125  0.0  -0.125  -1.0
       -1.25  -0.25   0.0   0.25    1.25
       -2.0    0.125  0.0  -0.125   2.0
     ]

end

@testset "inner product" begin
    z = ZernikeBW(5, 4);
    el = z.elements;
    @test PhaseBases.inner(el[1],el[1]) == 25.0
    @test PhaseBases.inner(el[1], el[2]) == 0

    a = reshape(1:27,(3,3,3))
    b= VectorOfArray([a[:,:,i] for i =1:last( size(a))])
    @test PhaseBases.inner(b,[1,100,10000]) == PhaseBases.inner([1,100,10000],b) == [
        191001  221304  251607
        201102  231405  261708
        211203  241506  271809   ]

end

@testset "aperture" begin
dom1= PhaseBases.CartesianDomain2D(-1:.2:1, -.9:.2:.9)
ap, mask = PhaseBases.aperture(dom1, 1);
ap, mask = PhaseBases.aperture(dom1, 1);
@test ap == [
    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  1  0  0  0  0  0
    0  0  0  1  1  1  1  1  0  0  0
    0  0  0  1  1  1  1  1  0  0  0
    0  0  0  1  1  1  1  1  0  0  0
    0  0  0  1  1  1  1  1  0  0  0
    0  0  0  0  0  1  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

# The sampled basis is not orthogonal by default!
zbig = ZernikeBW(512, 10);
m = PhaseBases.innermatrix(zbig);
norms = sqrt.(diag(m));
[zbig.elements[i] *= 1 / norms[i] for i in eachindex(zbig.elements)];
m = PhaseBases.innermatrix(zbig);
using ZChop
m = zchop(m);
r = m - I(length(zbig.elements));
maximum(abs.(r))


text:
using StanModels, Distributions

# ### snippet 8.25

N = 100                                          # individuals
height  = rand(Normal(10,2), N) ; # sim total height of each
leg_prop = rand(Uniform(0.4,0.5), N); # leg as proportion of height

# sim left leg as proportion + error
leg_left = leg_prop .* height .+  rand(Normal( 0 , 0.02 ), N);
# sim right leg as proportion + error
leg_right = leg_prop .* height .+  rand(Normal( 0 , 0.02 ), N);

# combine into data frame

df =  DataFrame(height=height, leg_left = leg_left, leg_right = leg_right);

# Show first 5 rows

first(df, 5)

# Define the Stan language model

m8_8s = "
data{
    int N;
    real height[N];
    real leg_right[N];
    real leg_left[N];
}
parameters{
    real a;
    real bl;
    real br;
    real sigma;
}
model{
    vector[N] mu;
    sigma ~ cauchy( 0 , 1 );
    br ~ normal( 2 , 10 );
    bl ~ normal( 2 , 10 );
    a ~ normal( 10 , 100 );
    for ( i in 1:100 ) {
        mu[i] = a + bl * leg_left[i] + br * leg_right[i];
    }
    height ~ normal( mu , sigma );
}
";

# Define the Stanmodel and set the output format to :mcmcchains.

sm = SampleModel("m8.8s", m8_8s);

# Input data for cmdstan

m8_8_data = Dict("N" => size(df, 1), "height" => df[!, :height],
    "leg_left" => df[!, :leg_left], "leg_right" => df[!, :leg_right]);

# Sample using cmdstan

(sample_file, log_file) = stan_sample(sm, data=m8_8_data);

# Describe the draws

if !(sample_file == nothing)
  chn = read_samples(sm)
  describe(chn)
end