text:
@testset "Promotion" begin
    for T in [Mag, MagRef]
        @test promote_type(T, Arf) == promote_type(T, ArfRef) == Arf
        @test promote_type(T, Arb) == promote_type(T, ArbRef) == Arb
        @test promote_type(T, Acb) == promote_type(T, AcbRef) == Acb
        @test promote_type(T, ArbSeries) == ArbSeries
        @test promote_type(T, AcbSeries) == AcbSeries
    end
    @test promote_type(MagRef, Mag) == promote_type(Mag, MagRef) == Mag

    for T in [Arf, ArfRef]
        @test promote_type(T, Mag) == promote_type(T, MagRef) == Arf
        @test promote_type(T, Arb) == promote_type(T, ArbRef) == Arb
        @test promote_type(T, Acb) == promote_type(T, AcbRef) == Acb
        @test promote_type(T, ArbSeries) == ArbSeries
        @test promote_type(T, AcbSeries) == AcbSeries
    end
    @test promote_type(ArfRef, Arf) == promote_type(Arf, ArfRef) == Arf

    for T in [Arb, ArbRef]
        @test promote_type(T, Mag) == promote_type(T, MagRef) == Arb
        @test promote_type(T, Arf) == promote_type(T, ArfRef) == Arb
        @test promote_type(T, Acb) == promote_type(T, AcbRef) == Acb
        @test promote_type(T, ArbSeries) == ArbSeries
        @test promote_type(T, AcbSeries) == AcbSeries
    end
    @test promote_type(ArbRef, Arb) == promote_type(Arb, ArbRef) == Arb

    for T in [Acb, AcbRef]
        @test promote_type(T, Mag) == promote_type(T, MagRef) == Acb
        @test promote_type(T, Arf) == promote_type(T, ArfRef) == Acb
        @test promote_type(T, Arb) == promote_type(T, ArbRef) == Acb
        @test promote_type(T, ArbSeries) == AcbSeries
        @test promote_type(T, AcbSeries) == AcbSeries
    end
    @test promote_type(AcbRef, Acb) == promote_type(Acb, AcbRef) == Acb

    @test promote_type(ArbSeries, Mag) == promote_type(ArbSeries, MagRef) == ArbSeries
    @test promote_type(ArbSeries, Arf) == promote_type(ArbSeries, ArfRef) == ArbSeries
    @test promote_type(ArbSeries, Arb) == promote_type(ArbSeries, ArbRef) == ArbSeries
    @test promote_type(ArbSeries, Acb) == promote_type(ArbSeries, AcbRef) == AcbSeries
    @test promote_type(ArbSeries, AcbSeries) == AcbSeries

    @test promote_type(AcbSeries, Mag) == promote_type(AcbSeries, MagRef) == AcbSeries
    @test promote_type(AcbSeries, Arf) == promote_type(AcbSeries, ArfRef) == AcbSeries
    @test promote_type(AcbSeries, Arb) == promote_type(AcbSeries, ArbRef) == AcbSeries
    @test promote_type(AcbSeries, Acb) == promote_type(AcbSeries, AcbRef) == AcbSeries
    @test promote_type(AcbSeries, ArbSeries) == AcbSeries

    for T in [Mag, MagRef, Arf, ArfRef, Arb, ArbRef, Acb, AcbRef, ArbSeries, AcbSeries]
        @test promote_type(T, Float64) == Arblib._nonreftype(T)
    end

    for T in [Arf, ArfRef, Arb, ArbRef, Acb, AcbRef, ArbSeries, AcbSeries]
        @test promote_type(T, BigFloat) == Arblib._nonreftype(T)
    end

    @test promote_type(Arb, ComplexF64) == Acb
    @test promote_type(ArbRef, ComplexF64) == Acb
    @test promote_type(ArbSeries, ComplexF64) == AcbSeries
end


text:
using FrameFun.ApproximationProblems, FrameFun.ExtensionFrames, Test, BasisFunctions, FrameFun.Platforms, DomainSets

@testset "ApproximationProblems" begin
    ap1 = approximationproblem(Fourier(10))
    ap2 = approximationproblem(Fourier(10)^2)
    ap3 = approximationproblem(platform(Fourier(10)),10)
    ap4 = approximationproblem(Fourier(10), 0.0..0.5)

    @test length(elements(ap2))==2

    for ap in (ap1, ap2, ap3, ap4)
        @test samplingparam(ap) == nothing
    end

    for ap in (ap1, ap2, ap3, ap4)
        if ap isa ProductPlatformApproximation
            setsamplingparam!(ap, (10,10))
        else
            setsamplingparam!(ap, 10)
        end
    end

    @test SamplingStyle(ap1) == InterpolationStyle()
    @test SamplingStyle(ap2) == ProductSamplingStyle(InterpolationStyle(),InterpolationStyle())
    @test SamplingStyle(ap3) == InterpolationStyle()
    @test SamplingStyle(ap4) == OversamplingStyle()

    @test SolverStyle(SamplingStyle(ap1),ap1) == InverseStyle()
    @test SolverStyle(SamplingStyle(ap2),ap2) == ProductSolverStyle(InverseStyle(),InverseStyle())
    @test SolverStyle(SamplingStyle(ap3),ap3) == InverseStyle()
    @test SolverStyle(SamplingStyle(ap4),ap4) == DirectStyle()

    for ap in (ap1, ap2, ap3, ap4)
        @test samplingparam(ap) != nothing
    end

    @test dictionary(ap1) == Fourier(10)
    @test dictionary(ap2) == Fourier(10)^2
    @test dictionary(ap3) == Fourier(10)
    @test iscompatible(dictionary(ap4), extensionframe(Fourier(10), 0.0..0.5))

    ap1 = approximationproblem(platform(Fourier(10)))
    ap2 = approximationproblem(platform(Fourier(10)^2))
    ap3 = approximationproblem(platform(extensionframe(Fourier(10), 0.0..0.5)))

    for ap in (ap1,ap2,ap3)
        @test ap isa AdaptiveApproximation
    end

    @test dictionary(ap1,10)== Fourier(10)
    @test dictionary(ap2,10)== Fourier(10)^2
    @test iscompatible(dictionary(ap3,10),extensionframe(Fourier(10), 0.0..0.5))
end

ext: jl


text:
using Test
using StaticArrays, IntervalSets
import ClimaCore.DataLayouts: IJFH
import ClimaCore:
    Fields, Domains, Meshes, Topologies, Spaces, Operators, Geometry
using StaticArrays, IntervalSets, LinearAlgebra

FT = Float64

radius = FT(1)
domain = Domains.SphereDomain(radius)

sdiv = Operators.Divergence()
wdiv = Operators.WeakDivergence()

@testset "divergence" begin
    Ne = 7
    Nq = 6

    mesh = Meshes.EquiangularCubedSphere(domain, Ne)
    grid_topology = Topologies.Topology2D(mesh)

    quad = Spaces.Quadratures.GLL{Nq}()
    space = Spaces.SpectralElementSpace2D(grid_topology, quad)
    coords = Fields.coordinate_field(space)

    u_local = @. Geometry.UVVector(
        cosd(coords.long),
        -sind(coords.long) * sind(coords.lat),
    )

    u = Geometry.transform.(Ref(Geometry.Contravariant12Axis()), u_local)

    div_us = sdiv.(u)
    Spaces.weighted_dss!(div_us)

    div_uw = wdiv.(u)
    Spaces.weighted_dss!(div_uw)

    div_exact = @. -2 * sind(coords.long) * cosd(coords.lat) / radius

    @test div_us ≈ div_exact rtol = 1e-2
    @test div_uw ≈ div_exact rtol = 1e-2
end

"""
    convergence_rate(err, Δh)
Estimate convergence rate given vectors `err` and `Δh`
    err = C Δh^p+ H.O.T
    err_k ≈ C Δh_k^p
    err_k/err_m ≈ Δh_k^p/Δh_m^p
    log(err_k/err_m) ≈ log((Δh_k/Δh_m)^p)
    log(err_k/err_m) ≈ p*log(Δh_k/Δh_m)
    log(err_k/err_m)/log(Δh_k/Δh_m) ≈ p
"""
convergence_rate(err, Δh) =
    [log(err[i] / err[i - 1]) / log(Δh[i] / Δh[i - 1]) for i in 2:length(Δh)]

@testset "convergence tests for the divergence operator on the sphere" begin
    Nes = [3, 9, 27]
    Nqs = [4, 6]

    for (Iq, Nq) in enumerate(Nqs)
        err_sdiv = Array{FT}(undef, length(Nes))
        err_wdiv = Array{FT}(undef, length(Nes))
        Δh = Array{FT}(undef, length(Nes))

        for (Ie, Ne) in enumerate(Nes)
            mesh = Meshes.EquiangularCubedSphere(domain, Ne)
            grid_topology = Topologies.Topology2D(mesh)

            quad = Spaces.Quadratures.GLL{Nq}()
            space = Spaces.SpectralElementSpace2D(grid_topology, quad)
            coords = Fields.coordinate_field(space)

            u_local = @. Geometry.UVVector(
                cosd(coords.long),
                -sind(coords.long) * sind(coords.lat),
            )

            u =
                Geometry.transform.(
                    Ref(Geometry.Contravariant12Axis()),
                    u_local,
                )

            div_us = sdiv.(u)
            Spaces.weighted_dss!(div_us)

            div_uw = wdiv.(u)
            Spaces.weighted_dss!(div_uw)

            div_exact = @. -2 * sind(coords.long) * cosd(coords.lat) / radius

            err_sdiv[Ie] = norm(div_us .- div_exact)
            err_wdiv[Ie] = norm(div_uw .- div_exact)
            Δh[Ie] = 1 / Ne

        end

        convergence_rate_sdiv = convergence_rate(err_sdiv, Δh)
        convergence_rate_wdiv = convergence_rate(err_wdiv, Δh)

        for Ie in range(1, length = length(Nes) - 1)
            @test convergence_rate_sdiv[Ie] ≈ (Nq - 1) atol = 0.1
            @test convergence_rate_wdiv[Ie] ≈ (Nq - 1) atol = 0.1
        end
    end
end


text:
using Statistics, LinearAlgebra

# bonds of a petersen graph
const L1 = [(1, 6), (2, 7), (3, 8), (4, 9), (5, 10),
    (1, 2), (2, 3), (3, 4), (4, 5), (1, 5), (6, 8),
    (8, 10), (7, 10), (7, 9), (6, 9)]

# disconnected bonds of a petersen graph
const L2 = [(1, 3), (1, 4), (1, 7), (1, 8), (1, 9),
    (1, 10), (2, 4), (2, 5), (2, 6), (2, 8), (2, 9),
    (2, 10), (3, 5), (3, 6), (3, 7), (3, 9), (3, 10),
    (4, 6), (4, 7), (4, 8), (4, 10), (5, 6), (5, 7),
    (5, 8), (5, 9), (6, 7), (6, 10), (7, 8), (8, 9),
    (9, 10)]

function myvar_and_mean(v)
    mv  = Statistics.mean(v)
    sum((v .- mv).^2)./(length(v)-1), mv
end

function relu(x::T) where T
    x > 0 ? x : zero(T)
end

function embedding_loss(x)
    a = [LinearAlgebra.norm(x[:,i]-x[:,j]) for (i, j) in L1]
    b = [LinearAlgebra.norm(x[:,i]-x[:,j]) for (i, j) in L2]
    va, ma = myvar_and_mean(a)
    vb, mb = myvar_and_mean(b)
    va + vb + exp(relu(-mb + ma + 0.1)) - 1
end


text:
function check_point(M::Hyperbolic{N}, p::PoincareHalfSpacePoint; kwargs...) where {N}
    mpv = check_point(Euclidean(N), p.value; kwargs...)
    mpv === nothing || return mpv
    if !(last(p.value) > 0)
        return DomainError(
            norm(p.value),
            "The point $(p) does not lie on $(M) since its last entry is nonpositive.",
        )
    end
end

@doc raw"""
    convert(::Type{PoincareHalfSpacePoint}, p::PoincareBallPoint)

convert a point [`PoincareBallPoint`](@ref) `p` (from $ℝ^n$) from the
Poincaré ball model of the [`Hyperbolic`](@ref) manifold $\mathcal H^n$ to a [`PoincareHalfSpacePoint`](@ref) $π(p) ∈ ℝ^n$.
Denote by $\tilde p = (p_1,\ldots,p_{n-1})$. Then the isometry is defined by

````math
π(p) = \frac{1}{\lVert \tilde p \rVert^2 - (p_n-1)^2}
\begin{pmatrix}2p_1\\⋮\\2p_{n-1}\\1-\lVert p\rVert^2\end{pmatrix}.
````
"""
function convert(::Type{PoincareHalfSpacePoint}, p::PoincareBallPoint)
    return PoincareHalfSpacePoint(
        1 / (norm(p.value[1:(end - 1)])^2 + (last(p.value) - 1)^2) .*
        vcat(2 .* p.value[1:(end - 1)], 1 - norm(p.value)^2),
    )
end

@doc raw"""
    convert(::Type{PoincareHalfSpacePoint}, p::Hyperboloid)
    convert(::Type{PoincareHalfSpacePoint}, p)

convert a [`HyperboloidPoint`](@ref) or `Vector``p` (from $ℝ^{n+1}$) from the
Hyperboloid model of the [`Hyperbolic`](@ref) manifold $\mathcal H^n$ to a [`PoincareHalfSpacePoint`](@ref) $π(x) ∈ ℝ^{n}$.

This is done in two steps, namely transforming it to a Poincare ball point and from there further on to a PoincareHalfSpacePoint point.
"""
convert(::Type{PoincareHalfSpacePoint}, ::Any)
function convert(t::Type{PoincareHalfSpacePoint}, p::HyperboloidPoint)
    return convert(t, convert(PoincareBallPoint, p))
end
function convert(t::Type{PoincareHalfSpacePoint}, p::T) where {T<:AbstractVector}
    return convert(t, convert(PoincareBallPoint, p))
end

@doc raw"""
    convert(::Type{PoincareHalfSpaceTVector}, p::PoincareBallPoint, X::PoincareBallTVector)

convert a [`PoincareBallTVector`](@ref) `X` at `p` to a [`PoincareHalfSpacePoint`](@ref)
on the [`Hyperbolic`](@ref) manifold $\mathcal H^n$ by computing the push forward $π_*(p)[X]$ of
the isometry $π$ that maps from the Poincaré ball to the Poincaré half space,
cf. [`convert(::Type{PoincareHalfSpacePoint}, ::PoincareBallPoint)`](@ref).

The formula reads

````math
π_*(p)[X] =
\frac{1}{\lVert \tilde p\rVert^2 + (1-p_n)^2}
\begin{pmatrix}
2X_1\\
⋮\\
2X_{n-1}\\
-2⟨X,p⟩
\end{pmatrix}
-
\frac{2}{(\lVert \tilde p\rVert^2 + (1-p_n)^2)^2}
\begin{pmatrix}
2p_1(⟨X,p⟩-X_n)\\
⋮\\
2p_{n-1}(⟨X,p⟩-X_n)\\
(\lVert p \rVert^2-1)(⟨X,p⟩-X_n)
\end{pmatrix}
````
where $\tilde p = \begin{pmatrix}p_1\\⋮\\p_{n-1}\end{pmatrix}$.
"""
function convert(
    ::Type{PoincareHalfSpaceTVector},
    p::PoincareBallPoint,
    X::PoincareBallTVector,
)
    den = norm(p.value[1:(end - 1)])^2 + (last(p.value) - 1)^2
    scp = dot(p.value, X.value)
    c1 =
        (2 / den .* X.value[1:(end - 1)]) .-
        (4 * (scp - last(X.value)) / (den^2)) .* p.value[1:(end - 1)]
    c2 = -2 * scp / den - 2 * (1 - norm(p.value)^2) * (scp - last(X.value)) / (den^2)
    return PoincareHalfSpaceTVector(vcat(c1, c2))
end

@doc raw"""
    convert(::Type{PoincareHalfSpaceTVector}, p::HyperboloidPoint, ::HyperboloidTVector)
    convert(::Type{PoincareHalfSpaceTVector}, p::P, X::T) where {P<:AbstractVector, T<:AbstractVector}

convert a [`HyperboloidTVector`](@ref) `X` at `p` to a [`PoincareHalfSpaceTVector`](@ref)
on the [`Hyperbolic`](@ref) manifold $\mathcal H^n$ by computing the push forward $π_*(p)[X]$ of
the isometry $π$ that maps from the Hyperboloid to the Poincaré half space,
cf. [`convert(::Type{PoincareHalfSpacePoint}, ::HyperboloidPoint)`](@ref).

This is done similarly to the approach there, i.e. by using the Poincaré ball model as
an intermediate step.
"""
convert(::Type{PoincareHalfSpaceTVector}, ::Any)
function convert(
    t::Type{PoincareHalfSpaceTVector},
    p::HyperboloidPoint,
    X::HyperboloidTVector,
)
    return convert(t, convert(AbstractVector, p), convert(AbstractVector, X))
end
function convert(
    ::Type{PoincareHalfSpaceTVector},
    p::P,
    X::T,
) where {P<:AbstractVector,T<:AbstractVector}
    return convert(
        PoincareHalfSpaceTVector,
        convert(Tuple{PoincareBallPoint,PoincareBallTVector}, (p, X))...,
    )
end

@doc raw"""
    convert(
        ::Type{Tuple{PoincareHalfSpacePoint,PoincareHalfSpaceTVector}},
        (p,X)::Tuple{PoincareBallPoint,PoincareBallTVector}
    )

Convert a [`PoincareBallPoint`](@ref) `p` and a [`PoincareBallTVector`](@ref) `X`
to a [`PoincareHalfSpacePoint`](@ref) and a [`PoincareHalfSpaceTVector`](@ref) simultaneously,
see [`convert(::Type{PoincareHalfSpacePoint}, ::PoincareBallPoint)`](@ref) and
[`convert(::Type{PoincareHalfSpaceTVector}, ::PoincareBallPoint,::PoincareBallTVector)`](@ref)
for the formulae.
"""
function convert(
    ::Type{Tuple{PoincareHalfSpacePoint,PoincareHalfSpaceTVector}},
    (p, X)::Tuple{PoincareBallPoint,PoincareBallTVector},
)
    return (convert(PoincareHalfSpacePoint, p), convert(PoincareHalfSpaceTVector, p, X))
end

@doc raw"""
    convert(
        ::Type{Tuple{PoincareHalfSpacePoint,PoincareHalfSpaceTVector}},
        (p,X)::Tuple{HyperboloidPoint,HyperboloidTVector}
    )
    convert(
        ::Type{Tuple{PoincareHalfSpacePoint,PoincareHalfSpaceTVector}},
        (p, X)::Tuple{P,T},
    ) where {P<:AbstractVector, T <: AbstractVector}

Convert a [`HyperboloidPoint`](@ref) `p` and a [`HyperboloidTVector`](@ref) `X`
to a [`PoincareHalfSpacePoint`](@ref) and a [`PoincareHalfSpaceTVector`](@ref) simultaneously,
see [`convert(::Type{PoincareHalfSpacePoint}, ::HyperboloidPoint)`](@ref) and
[`convert(::Type{PoincareHalfSpaceTVector}, ::Tuple{HyperboloidPoint,HyperboloidTVector})`](@ref)
for the formulae.
"""
function convert(
    ::Type{Tuple{PoincareHalfSpacePoint,PoincareHalfSpaceTVector}},
    (p, X)::Tuple{HyperboloidPoint,HyperboloidTVector},
)
    return (convert(PoincareHalfSpacePoint, p), convert(PoincareHalfSpaceTVector, p, X))
end
function convert(
    ::Type{Tuple{PoincareHalfSpacePoint,PoincareHalfSpaceTVector}},
    (p, X)::Tuple{P,T},
) where {P<:AbstractVector,T<:AbstractVector}
    return (convert(PoincareHalfSpacePoint, p), convert(PoincareHalfSpaceTVector, p, X))
end

@doc raw"""
    distance(::Hyperbolic, p::PoincareHalfSpacePoint, q::PoincareHalfSpacePoint)

Compute the distance on the [`Hyperbolic`](@ref) manifold $\mathcal H^n$ represented in the
Poincaré half space model. The formula reads

````math
d_{\mathcal H^n}(p,q) = \operatorname{acosh}\Bigl( 1 + \frac{\lVert p - q \rVert^2}{2 p_n q_n} \Bigr)
````
"""
function distance(::Hyperbolic, p::PoincareHalfSpacePoint, q::PoincareHalfSpacePoint)
    return acosh(1 + norm(p.value .- q.value)^2 / (2 * p.value[end] * q.value[end]))
end

@doc raw"""
    inner(
        ::Hyperbolic{n},
        p::PoincareHalfSpacePoint,
        X::PoincareHalfSpaceTVector,
        Y::PoincareHalfSpaceTVector
    )

Compute the inner product in the Poincaré half space model. The formula reads
````math
g_p(X,Y) = \frac{⟨X,Y⟩}{p_n^2}.
````
"""
function inner(
    ::Hyperbolic,
    p::PoincareHalfSpacePoint,
    X::PoincareHalfSpaceTVector,
    Y::PoincareHalfSpaceTVector,
)
    return dot(X.value, Y.value) / last(p.value)^2
end

@doc raw"""
    project(::Hyperbolic, ::PoincareHalfSpacePoint ::PoincareHalfSpaceTVector)

projction of tangent vectors in the Poincaré half space model is just the identity, since
the tangent space consists of all $ℝ^n$.
"""
project(::Hyperbolic, ::PoincareHalfSpacePoint::PoincareHalfSpaceTVector)

function allocate_result(
    ::Hyperbolic,
    ::typeof(project),
    X::PoincareHalfSpaceTVector,
    ::PoincareHalfSpacePoint,
)
    return PoincareHalfSpaceTVector(allocate(X.value))
end

function project!(
    ::Hyperbolic,
    Y::PoincareHalfSpaceTVector,
    ::PoincareHalfSpacePoint,
    X::PoincareHalfSpaceTVector,
)
    return (Y.value .= X.value)
end


text:
using MonkeyLang
using Test
using IOCapture

const m = MonkeyLang

@testset "MonkeyLang.jl" begin
    include("test_helpers.jl")

    include("lexer_test.jl")
    include("ast_test.jl")
    include("parser_test.jl")
    include("object_test.jl")

    include("quote_test.jl")
    include("ast_modify_test.jl")

    include("code_test.jl")
    include("symbol_table_test.jl")
    include("analyzer_test.jl")
    include("compiler_test.jl")

    include("backend_test.jl")
    include("evaluator_test.jl")
    include("vm_test.jl")
    include("transpilers/transpiler_test.jl")

    include("repl_test.jl")
    include("cli_test.jl")
end


text:
# Attention mechanisms:
#
"""
    abstract type AttentionMechanism

Attention mechanisms follow the same interface and common signatures.

If possible, the algorithm allows precomputing of the projections
of the context vector
generated by the encoder in a encoder-decoder-architecture
(i.e. in case of an RNN encoder the accumulated encoder hidden states).

By default attention scores are scaled according to
Vaswani et al., 2017 *(Vaswani et al., Attention Is All You Need,
CoRR, 2017)*.

All algorithms use soft attention.

## Constructors:
    Attn*Mechanism*(dec_units, enc_units; scale=true)
    Attn*Mechanism*(units; scale=true)

The one-argument version can be used, if encoder dimensions and decoder
dimensions are the same.

## Common Signatures:
    function (attn::AttentionMechanism)(h_t, h_enc; reset=false)
    function (attn::AttentionMechanism)(; reset=false)

### Arguments:
+ `h_t`:    decoder hidden state. If ``h_t`` is a vector, its length
            equals the number of decoder units. If it is a matrix,
            ``h_t`` includes the states for a minibatch of samples and has
            the size [units, mb].
+ `h_enc`:  encoder hidden states, 2d or 3d. If ``h_{enc}`` is a
            matrix [units, steps] with the hidden states of all encoder steps.
            If 3d: [units, mb, steps] encoder states for all minibatches.
+ `reset=false`: If the keyword argument is set to `true`, projections of
            the encoder states are computed. By default projections are
            stored in the object and reused until the object is resetted.
            For attention mechanisms that don't allow precomputation
            the argument is ignored.

The short form `(::AttentionMechanism)(reset=true)` can be used to reset
the precomputed projections.

### Return values
All functions return `c` and `α` where α is a matrix of size [mb,steps]
with the attention factors for each step and minibatch.
`c` is a matrix of size [units, mb] with the context vector for
each sample of the minibatch, calculated as the α-weighted sum of
all encoder hidden states ``h_{enc}`` for each minibatch.

## Attention Mechanisms:

All attention mechanisms calculate attention factors α from scores
derived from projections of the encoder hidden states:
```math
\\alpha = \\mathrm{softmax}(\\mathrm{score}(h_{enc},h_{t}) \\cdot 1/\\sqrt{n}))
```

Attention mechanisms implemented:
"""
abstract type AttentionMechanism
end

function (attn::AttentionMechanism)(;reset=false)

    if reset && hasfield(typeof(attn), :projections)
        attn.projections = nothing
    end
end

"""
    mutable struct AttnBahdanau <: AttentionMechanism

Bahdanau-style (additive, concat) attention mechanism according
to the paper:

*D. Bahdanau, KH. Co, Y. Bengio,
Neural Machine Translation by jointlylearning to align and translate,
ICLR, 2015*.

```math
\\mathrm{score}(h_{t},h_{enc}) = v_{a}^{\\top}\\cdot\\tanh(W[h_{t},h_{enc}])
```

### Constructors:
    AttnBahdanau(dec_units, enc_units; scale=true)
    AttnBahdanau(units; scale=true)
"""
mutable struct AttnBahdanau <: AttentionMechanism
    enc
    dec
    combine
    scale
    projections
    AttnBahdanau(dec_units, enc_units; scale=true) = new(Linear(enc_units, dec_units, bias=false),
                                Linear(dec_units, dec_units, bias=false),
                                Linear(dec_units, 1, bias=false),
                                scale ? 1/sqrt(enc_units) : 1.0,
                                nothing)
    AttnBahdanau(units; scale=true) = AttnBahdanau(units, units,
                                scale=scale)
end

function (attn::AttnBahdanau)(h_t, h_enc; reset=false)
                                    # h_t is a (n_units, <n_mb>),
                                    # h_enc is (n_units, <n_mb>, n_steps)
    # make all 3d:
    #
    h_encR = reshape(h_enc, size(h_enc)[1], :, size(h_enc)[ndims(h_enc)])
    units, mb, steps = size(h_encR)
    h_tR = reshape(h_t, size(h_t)[1], :)

    # this is possible, because the Linear-layers are used:
    #
    if reset || attn.projections == nothing
        attn.projections = attn.enc(h_encR)
    end
    score = attn.combine(tanh.(attn.projections .+ attn.dec(h_tR)))
    score *= attn.scale
    α = softmax(score, dims=3)

    # calc. context from encoder states:
    #
    c = sum(α .* h_encR, dims=3)

    # remove unneeded dims:
    #
    c = reshape(c, units, mb)
    α = reshape(α, mb, steps)
    return c, α
end

"""
    mutable struct AttnLuong <: AttentionMechanism

Luong-style (multiplicative) attention mechanism according to
the paper (referred as *General*-type attention):
*M.-T. Luong, H. Pham, C.D. Manning,
Effective Approaches to Attention-based Neural Machine Translation,
CoRR, 2015*.
```math
\\mathrm{score}(h_{t},h_{enc}) = h_{t}^{\\top} W h_{enc}
```

### Constructors:
    AttnLuong(dec_units, enc_units; scale=true)
    AttnLuong(units; scale=true)
"""
mutable struct AttnLuong <: AttentionMechanism
    enc
    scale
    projections
    AttnLuong(dec_units, enc_units; scale=true) = new(Linear(enc_units,
                dec_units, bias=false),
                scale ? 1/sqrt(enc_units) : 1.0,
                nothing)
    AttnLuong(units; scale=true) = AttnLuong(units, units, scale=scale)
end

function (attn::AttnLuong)(h_t, h_enc; reset=false)
    # make all 3d:
    #
    h_encR = reshape(h_enc, size(h_enc)[1], :, size(h_enc)[ndims(h_enc)])
    units, mb, steps = size(h_encR)
    h_tR = reshape(h_t, size(h_t)[1], :)

    if reset || attn.projections == nothing
        attn.projections = attn.enc(h_encR)
    end

    score = sum(attn.projections .* h_tR, dims=1)
    score *= attn.scale
    α = softmax(score, dims=3)

    # calc. context from encoder states:
    #
    c = sum(α .* h_encR, dims=3)

    # remove unneeded dims:
    #
    c = reshape(c, units, mb)
    α = reshape(α, mb, steps)
    return c, α
end

"""
    mutable struct AttnDot <: AttentionMechanism

Dot-product attention (without trainable parameters)
according to the Luong, et al. (2015) paper.

``\\mathrm{score}(h_{t},h_{enc}) = h_{t}^{\\top} h_{enc}``

### Constructors:
    AttnDot(; scale=true)
"""
mutable struct AttnDot <: AttentionMechanism
    scale
    AttnDot(;scale=true) = new(scale)
end

function (attn::AttnDot)(h_t, h_enc; reset=false)
    # make all 3d:
    #
    h_encR = reshape(h_enc, size(h_enc)[1], :, size(h_enc)[ndims(h_enc)])
    units, mb, steps = size(h_encR)
    h_tR = reshape(h_t, size(h_t)[1], :)

    score = sum(h_encR .* h_tR, dims=1)
    if attn.scale
        score = score ./ Float32(sqrt(units))
    end
    α = softmax(score, dims=3)

    # calc. context from encoder states:
    #
    c = sum(α .* h_encR, dims=3)

    # remove unneeded dims:
    #
    c = reshape(c, units, mb)
    α = reshape(α, mb, steps)
    return c, α
end

"""
    mutable struct AttnLocation <: AttentionMechanism

Location-based attention that only depends on the current
decoder state ``h_t`` and not on the encoder states,
according to the Luong, et al. (2015) paper.

``\\mathrm{score}(h_{t}) = W h_{t}``

### Constructors:
    AttnLocation(len, dec_units; scale=true)

+ `len`: maximum sequence length of the encoder to be considered
        for attention. If the actual length of ``h_{enc}`` is bigger as the
        length of α, attention factors for the remaining states are set to
        0.0. If the actual length of h_enc is smaller than α, only the matching
        attention factors are applied.
+ `dec_units`: number of decoder units.
"""
mutable struct AttnLocation <: AttentionMechanism
    dec
    len
    scale
    AttnLocation(len, dec_units; scale=true) = new(Linear(
                                                    dec_units, len, bias=false),
                                                len, scale)
end

function (attn::AttnLocation)(h_t, h_enc)
    # make all 3d:
    #
    h_encR = reshape(h_enc, size(h_enc)[1], :, size(h_enc)[ndims(h_enc)])
    units, mb, steps = size(h_encR)
    h_tR = reshape(h_t, size(h_t)[1], :)

    score = attn.dec(h_tR)
    if attn.scale
        score = score ./ sqrt(units)
    end
    α = softmax(score, dims=1)

    if attn.len > steps
        α = α[1:steps,:]
    elseif attn.len < steps
        α = vcat(init0(steps-attn.len, mb), α)
    end

    α = permutedims(α, (2,1))
    α = reshape(α, :, mb, steps)

    # calc. context from encoder states:
    #
    c = sum(α .* h_encR, dims=3)

    # remove unneeded dims:
    #
    c = reshape(c, units, mb)
    α = reshape(α, mb, steps)
    return c, α
end

"""
    mutable struct AttnInFeed <: AttentionMechanism

Input-feeding attention that depends on the current
decoder state ``h_t`` and the next input to the decoder ``i_{t+1}``,
according to the Luong, et al. (2015) paper.

Infeed attention provides a semantic attention that depends on the
next input token.

``\\mathrm{score}(h_{t}, i_{t+1}) = W_h h_{t} + W_i i_{t+1} = W [h_t, i_{t+1}]``

### Constructors:
    AttnInFeed(len, dec_units, fan_in; scale=true)

+ `len`: maximum sequence length of the encoder to be considered
        for attention. If the actual length of ``h_{enc}`` is bigger as the
        length of α, attention factors for the remaining states are set to
        0.0. If the actual length of h_enc is smaller than α, only the matching
        attention factors are applied.
+ `dec_units`: number of decoder units.
+ `fan_in`: size of the decoder input.

### Signature:
    function (attn::AttnInFeed)(h_t, inp, h_enc)

+ `h_t`:    decoder hidden state. If ``h_t`` is a vector, its length
            equals the number of decoder units. If it is a matrix,
            ``h_t`` includes the states for a minibatch of samples and has
            the size [units, mb].
+ `inp`: next decoder input ``i_{t+1}``
            (e.g. next embedded token of sequence)
+ `h_enc`:  encoder hidden states, 2d or 3d. If ``h_{enc}`` is a
            matrix [units, steps] with the hidden states of all encoder steps.
            If 3d: [units, mb, steps] encoder states for all minibatches.
"""
mutable struct AttnInFeed <: AttentionMechanism
    dec
    emb
    len
    scale
    AttnInFeed(len, dec_units, fan_in; scale=true) =
                    new(Linear(dec_units, len, bias=false),
                        Linear(fan_in, len, bias=false),
                        len, scale)
end

function (attn::AttnInFeed)(h_t, inp, h_enc)
    # make all 3d:
    #
    h_encR = reshape(h_enc, size(h_enc)[1], :, size(h_enc)[ndims(h_enc)])
    units, mb, steps = size(h_encR)
    h_tR = reshape(h_t, size(h_t)[1], :)

    score = attn.dec(h_tR) .+ attn.emb(inp)
    if attn.scale
        score = score ./ sqrt(units)
    end
    α = softmax(score, dims=1)

    if attn.len > steps
        α = α[1:steps,:]
    elseif attn.len < steps
        α = vcat(init0(steps-attn.len, mb), α)
    end

    α = permutedims(α, (2,1))
    α = reshape(α, :, mb, steps)

    # calc. context from encoder states:
    #
    c = sum(α .* h_encR, dims=3)

    # remove unneeded dims:
    #
    c = reshape(c, units, mb)
    α = reshape(α, mb, steps)
    return c, α
end


text:
# Autogenerated wrapper script for Shoco_jll for armv6l-linux-musleabihf
export libshoco

JLLWrappers.@generate_wrapper_header("Shoco")
JLLWrappers.@declare_library_product(libshoco, "libshoco.so")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libshoco,
        "lib/libshoco.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
"""
    create_system(xx,yy,Tc,Th,coupling,
                  p;potential_type=:tilted_channel) -> (u0,params)

Initializes a system to be solved by the finite volume solvers in
BrownianDynamics. `u0` is an initial condition for the system and `params`
is the parameters of the system of type `FVMParameters`.
"""
function create_system(xx,yy,Tc,Th,coupling,p;potential_type=:tilted_channel)
    if potential_type == :tilted_channel
        potential = (x,y) -> -p[1]*x+p[2]cos(2π*x)+p[3]*y^2
    elseif potential_type == :curved_channel
        potential = (x,y) -> (-p[1]*x-p[2]*exp(-0.5*(-sin(2π*x)-sin(2π*(0.5y+0.5))+2cos(π*(0.5y+0.5)))^2)
                              +p[3]*(y)^4-p[4]*cos(2π*x)*cos(2π*(0.5y+0.5)))
    elseif potential_type == :double_channel
        potential = (x,y) -> ( exp(-10*(y+0.5)^2)*(-p[4]+p[2]*cos(2π*(x+π/4)))
                              +exp(-10*(y-0.5)^2)*(-p[5]+p[3]*cos(2π*(x-π/4)))
                              -p[1]*x+p[6]*y^2)
    end

    A = 1.0
    density_init = (x,y) -> exp(-potential(x,y))
    temperature_init = (x,y) -> ((Tc-Th)/(yy[end]-yy[1]))*y+Th-((Tc-Th)/(yy[end]-yy[1]))*yy[1]

    params = create_params((xx,yy),potential,A,coupling,temperature_init)
    u0 = create_initial_conditions((xx,yy),density_init,temperature_init)

    u0,params
end

function create_params(mesh,potential,A,coupling,temperature_init)
    xx,yy = mesh
    nx = length(xx)
    ny = length(yy)
    dx,dy = xx[2]-xx[1],yy[2]-yy[1]
    # Discrete potential. We will need the potential at the center of the volumes, as well
    # as the edges. The gradient of the potential only needs to be evaluated at the edges
    # of the cells.
    V = [potential(x,y) for x in xx, y in yy]
    V_x = [ForwardDiff.derivative(xp->potential(xp,y),x-0.5dx) for x in [xx;xx[end]+dx], y in [yy;yy[end]+dy]]
    V_y = [ForwardDiff.derivative(yp->potential(x,yp),y-0.5dy) for x in [xx;xx[end]+dx], y in [yy;yy[end]+dy]]
    Vxshift = [potential(x-0.5dx,y) for x in [xx;xx[end]+dx], y in [yy;yy[end]+dy]]
    Vyshift = [potential(x,y-0.5dy) for x in [xx;xx[end]+dx], y in [yy;yy[end]+dy]]
    Pmat = OffsetArray{eltype(xx)}(undef,0:nx+1,0:ny+1)
    fill!(Pmat,zero(eltype(xx)))
    Pmat[0:nx+1,0] *= 0
    Pmat[0:nx+1,ny+1] *= 0
    Tmat = OffsetArray{eltype(xx)}(undef,0:nx+1,0:ny+1)
    fill!(Tmat,zero(eltype(xx)))
    Tmat[0:nx+1,0:ny+1] = [temperature_init(x,y) for x in [xx[1]-dx;xx;xx[end]+dx],
                                                     y in [yy[1]-dy;yy;yy[end]+dy]]
    Jx = Array{eltype(xx)}(undef,nx+1,ny+1)
    Jy = Array{eltype(xx)}(undef,nx+1,ny+1)

    # jac = spzeros(eltype(xx),2nx*ny,2nx*ny)
    # jac_prototype = spzeros(eltype(xx),nx*ny,nx*ny)
    # jac_prototype[diagind(jac_prototype,0)] .= one(eltype(xx))*eps()
    # jac_prototype[diagind(jac_prototype,1)] .= one(eltype(xx))*eps()
    # jac_prototype[diagind(jac_prototype,-1)] .= one(eltype(xx))*eps()
    # jac_prototype[diagind(jac_prototype,nx*ny-1)] .= one(eltype(xx))*eps()
    # jac_prototype[diagind(jac_prototype,-nx*ny+1)] .= one(eltype(xx))*eps()
    # jac_prototype[diagind(jac_prototype,nx*ny-nx)] .= one(eltype(xx))*eps()
    # jac_prototype[diagind(jac_prototype,-nx*ny+nx)] .= one(eltype(xx))*eps()
    # jac_prototype[diagind(jac_prototype,nx*ny-1)] .= one(eltype(xx))*eps()
    # jac_prototype[diagind(jac_prototype,-nx*ny+1)] .= one(eltype(xx))*eps()
    # jac_prototype[diagind(jac_prototype,nx*ny-nx)] .= one(eltype(xx))*eps()
    # jac_prototype[diagind(jac_prototype,-nx*ny+nx)] .= one(eltype(xx))*eps()
    jac_prototype = spdiagm(0 => ones(eltype(xx),nx*ny),
                            1 => ones(eltype(xx),nx*ny-1),
                            -1 => ones(eltype(xx),nx*ny-1),
                            nx => ones(eltype(xx),nx*ny-nx),
                            -nx => ones(eltype(xx),nx*ny-nx),
                            nx*ny-nx => ones(eltype(xx),nx),
                            -nx*ny+nx => ones(eltype(xx),nx),
                            nx*ny-1 => ones(eltype(xx),1),
                            -nx*ny+1 => ones(eltype(xx),1))
    jac = vcat(hcat(jac_prototype,jac_prototype),hcat(jac_prototype,jac_prototype))

    FVMParameters(Pmat,Tmat,jac,Jx,Jy,V,Vxshift,Vyshift,V_x,V_y,dx,dy,A,coupling)
end

function create_initial_conditions(mesh,density_init,temperature_init)
    xx,yy = mesh
    nx = length(xx)
    ny = length(yy)
    dx,dy = xx[2]-xx[1],yy[2]-yy[1]
    # Discrete versions of the temperature and the probability.
    P = [density_init(x,y) for x in xx, y in yy]
    P /= sum(P)*dx*dy  # Normalize probability.
    T = [temperature_init(x,y) for x in xx, y in yy]

    u0 = [P[:];T[:]]
end

function evolve_to_steady_state!(integrator,params;steadytol=1e-3,maxiters=500,
                                                   convergence_warning=true)
    du = Array{eltype(integrator.u)}(undef,size(integrator.u))
    iters = 0
    # TODO right now du will be calculated somewhere in the step! function.
    # I then recalculate du, which is a bit of a waste.
    while ((norm(integrator.f(du,integrator.u,params,0))/size(integrator.u)[1] > steadytol)
            && (iters < maxiters))
        step!(integrator)
        iters += 1
    end
    (iters > maxiters && convergence_warning) && println("Maximum number of iterations reached, exiting.")

    integrator.u
end

function solve_steady_state(integrator,params::SpectralParameters{T};
                            steadytol=1e-3,maxiters=50) where T <: Number
    u = copy(integrator.u)
    nn = round(Int,length(u)/2)
    nd2 = round(Int,nn/2)
    du = similar(u)
    # We will use the boundary jacobian to assert that `u[nd2+1] == 1`.
    boundary_jac = zeros(eltype(u),2nn)
    boundary_jac[nd2+1] = 1
    jac = Array{eltype(u)}(undef,length(u),length(u))
    du[:] .= integrator.f(du,u,params,0)
    jac[:,:] .= integrator.f(Val{:jac},jac,u,params,0)
    # Newton's method.
    iters = 0
    residual_list = Float64[]
    while (norm([du;u[nd2+1]-1])/nn>steadytol) && (iters<maxiters)
        du[:] .= integrator.f(du,u,params,0)
        push!(residual_list,norm(du)/nn)
        jac[:,:] .= integrator.f(Val{:jac},jac,u,params,0)

        u[:] .+= -(([jac boundary_jac;boundary_jac' 0.0])\[du;u[nd2+1]-1])[1:2nn]
        iters += 1
    end
    iters > maxiters && println("Maximum number of iterations reached, exiting.")

    u
end

function solve_steady_state_uncoupled(integrator,params;
                                      steadytol=1e-3,maxiters=10,print_residual=true)
    Pmat,Tmat,Jx,Jy,V_x,V_y,dx,dy = params
    u = copy(integrator.u)
    nn = round(Int,length(u))

    du = similar(u)
    # We will use the boundary jacobian to assert that `sum(u[1:nx*ny])*dx*dy == 1`.
    boundary_jac = ones(nn)
    jac = spzeros(nn,nn)
    du[:] .= integrator.f(du,u,params,0)
    jac[:,:] .= integrator.f(Val{:jac},jac,u,params,0)
    # Newton's method.
    iters = 0
    while (norm([du;sum(u[1:nn])*dx*dy-1])/nn>steadytol) && (iters<maxiters)
        jac[:,:] .= integrator.f(Val{:jac},jac,u,params,0)

        du[:] .= integrator.f(du,u,params,0)
        u[:] .+= -(([jac boundary_jac;boundary_jac' 0.0])\[du;sum(u)*dx*dy-1])[1:nx*ny]
        iters += 1
    end
    iters > maxiters && println("Maximum number of iterations reached, exiting.")
    print_residual && @show norm([du;sum(u[1:nn])*dx*dy-1])/nn

    u
end

function solve_steady_state(u_init,params::FVMParameters;steadytol::Float64=1e-10,
                            maxiters::Int=5,print_residual::Bool=true,autodiff=false)
    dx,dy = params.dx,params.dy
    jac = params.jac_
    u = copy(u_init)
    nn = round(Int,length(u)/2)
    nx = round(Int,sqrt(nn))

    du = similar(u)
    # We will use the boundary jacobian to assert that `sum(u[1:nx*ny])*dx*dy == 1`.
    boundary_jac = zeros(eltype(u),2nn)
    boundary_jac[1:nn] .= one(eltype(u))
    du[:] .= flux!(du,u,params,0)
    # Autodiff code.
    if autodiff
        @unpack Pmat,Tmat,jac_,Jx,Jy,V,Vxshift,Vyshift,V_x,V_y,dx,dy,A,coupling = params
        params_autodiff = (V,Vxshift,Vyshift,V_x,V_y,dx,dy,A,coupling)
        f_autodiff = (du,u) -> flux_autodiff!(du,u,params_autodiff,0)
        tmp = ForwardDiff.jacobian!(jac,f_autodiff,du,u)
    end
    # Newton's method.
    du[:] .= flux!(du,u,params,0)
    iters = zero(Int64)
    while (norm([du;sum(u[1:nn])*dx*dy-1])/nn>steadytol) && (iters<maxiters)
        if autodiff
            jac[:,:] .= ForwardDiff.jacobian!(jac,f_autodiff,du,u)
            flux_autodiff!(du,u,params_autodiff,0)
        else
            flux!(Val{:jac},jac,u,params,0)
            flux!(du,u,params,0)
        end

        u[:] .+= -(([jac boundary_jac;boundary_jac' zero(eltype(u))])\[du;sum(u[1:nn])*dx*dy-one(eltype(u))])[1:2nn]
        iters += one(Int64)
    end
    iters > maxiters && println("Maximum number of iterations reached, exiting.")
    print_residual && @show norm([du;sum(u[1:nn])*dx*dy-1])/nn

    u
end

function stencil_indices!(inds,dims,row)
    n1,n2 = dims
    inds[:] .= SVector(row,row-1,row-n1,row+1,row+n2).-one(eltype(inds))
    i2s = CartesianIndices((n1,n2))
    i_j = i2s[row]
    i,j = i_j[1],i_j[2]
    if i==1
        inds[2] += n1-one(eltype(inds))
    elseif i==n1
        inds[4] -= n1-one(eltype(inds))
    end
    for i in 1:length(inds)
        inds[i] = mod(inds[i],n1*n2)+one(eltype(inds))
    end
    inds
end


text:
# world
abstract type World end
abstract type R2 <: World end
abstract type R3 <: World end

# friction cone
abstract type FrictionCone end
abstract type LinearizedCone <: FrictionCone end
abstract type NonlinearCone <: FrictionCone end

# environment
struct Environment{W <: World, F <: FrictionCone}
    surf::Any
	surf_grad::Any
end

function environment_2D_flat(; cone = LinearizedCone)
	Environment{R2,cone}(x -> 0.0, x-> zero(x))
end

function environment_3D_flat(; cone = LinearizedCone)
	Environment{R3,cone}(x -> 0.0, x-> zero(x))
end

function environment_2D(surf; cone = LinearizedCone)
	# Generate two functions: they both take as input a vector and return a vector.
	@variables q[1:1]
	s = surf(q)
	s = Symbolics.simplify.(s)
	ds = Symbolics.jacobian([s], q, simplify = true)
	ds = reshape(ds, 1)

	surf_fast = eval(Symbolics.build_function(s, q, checkbounds=true))
	surf_grad_fast = eval(Symbolics.build_function(ds, q, checkbounds=true)[1])

	Environment{R2,cone}(surf_fast, surf_grad_fast)
end

function environment_3D(surf; cone = LinearizedCone)
	# Generate two functions: they both take as input a vector and return a vector.
	@variables q[1:2]
	s = surf(q)
	s = Symbolics.simplify.(s)
	ds = Symbolics.jacobian([s], q, simplify = true)
	ds = reshape(ds, 2)

	surf_fast = eval(Symbolics.build_function(s, q, checkbounds=true))
	surf_grad_fast = eval(Symbolics.build_function(ds, q, checkbounds=true)[1])

	Environment{R3,cone}(surf_fast, surf_grad_fast)
end

function skew(x)
	SMatrix{3,3}([0.0 -x[3] x[2];
	               x[3] 0.0 -x[1];
				   -x[2] x[1] 0.0])
end

# rotation matrix rotating unit vector a onto unit vector b
function rot(a, b)
	v = cross(a, b)
	s = sqrt(transpose(v) * v)
	c = transpose(a) * b

	R = Diagonal(@SVector ones(3)) + skew(v) + 1.0 / (1.0 + c) * skew(v) * skew(v)
end

function rotation(env::Environment{R3,<:FrictionCone}, q)
	# unit surface normal (3D)
	n = [-1.0 * env.surf_grad(q[1:2]); 1.0]
	ns = n ./ sqrt(transpose(n) * n)

	# world-frame normal
	nw = @SVector [0.0, 0.0, 1.0]

	rot(ns, nw)
end

function rotation(env::Environment{R2,<:FrictionCone}, q)
	# unit surface normal (3D)
	sg = env.surf_grad(q[1:1])[1]
	n = [-1.0 * sg; 1.0]
	ns = n ./ sqrt(transpose(n) * n)

	# world-frame normal
	nw = @SVector [0.0, 1.0]

	ang = atan(nw[2], nw[1]) - atan(ns[2], ns[1])

	# Rw->s
	SMatrix{2,2}([cos(ang) -sin(ang); sin(ang) cos(ang)]) # this the rotation from world frame to surface frame
end

function rotation_3d(env::Environment{R3,<:FrictionCone}, q)
	r = rotation(env, q)
end

function rotation_3d(env::Environment{R2,<:FrictionCone}, q)
	r = rotation(env, q)
	r3d = [r[1,1] 0 r[1,2];
		     0    1   0   ;
		   r[2,1] 0 r[2,2]]
end

function friction_mapping(env::Environment{R2,LinearizedCone})
    SMatrix{1, 2}([1.0 -1.0])
end

function friction_mapping(env::Environment{R3,LinearizedCone})
    SMatrix{2, 4}([1.0 0.0 -1.0 0.0;
                   0.0 1.0 0.0 -1.0])
end

function friction_mapping(env::Environment{R2,NonlinearCone})
    SMatrix{1, 1}([1.0])
end

function friction_mapping(env::Environment{R3,NonlinearCone})
    SMatrix{2, 2}([1.0 0.0;
                   0.0 1.0])
end

dim(env::Environment{R2,<:FrictionCone}) = 2
dim(env::Environment{R3,<:FrictionCone}) = 3

friction_dim(env::Environment{R2,LinearizedCone}) = 2
friction_dim(env::Environment{R3,LinearizedCone}) = 4

friction_dim(env::Environment{R2,NonlinearCone}) = 1
friction_dim(env::Environment{R3,NonlinearCone}) = 2

function verify_2D_surface(env::Environment{R2,<:FrictionCone}, x;
	x_range = range(0.0, stop = 5.0, length = 1000))

	nw = [0.0; 1.0]
	tw = [1.0; 0.0]

	s = [env.surf([xx])[1] for xx in x_range]
	plt = plot(x_range, s, label = "", aspect_ratio = :equal)

	for t_cand in x
		s_cand = env.surf([t_cand])[1]
		p_cand = [t_cand; s_cand]
		tan_norm = rotation(env, [t_cand])' * tw
		norm_norm = rotation(env, [t_cand])' * nw

		plt = plot!([p_cand[1], p_cand[1] + tan_norm[1]],
			[p_cand[2], p_cand[2] + tan_norm[2]], color = :green, width = 2.0, label = "")
		plt = plot!([p_cand[1], p_cand[1] + norm_norm[1]],
			[p_cand[2], p_cand[2] + norm_norm[2]], color = :cyan, width = 2.0, label = "")
	end

	return plt
end


text:
module HTM

using Hyperscript
export @htm_str

const UENDTAG = "<//>"

include("util.jl")
include("parse.jl")

"""
    create_element(tag, attrs[, children...])

Create a DOM element.

This is an alternative syntax and (currently) serves as a rather trivial
abstraction layer inspired by
[`React.createElement`](https://pt-br.reactjs.org/docs/react-api.html#createelement).

```jldoctest
julia> HTM.create_element("div", [:class => "fruit"], "🍍")
<div class="fruit">🍍</div>
```
"""
@inline create_element(tag, attrs, children...) = Hyperscript.Node(Hyperscript.DEFAULT_HTMLSVG_CONTEXT, tag, children, attrs)

"""
    Node(tag[, attrs, children])

Compile time internal representation of a node.

```jldoctest
julia> HTM.Node(["div"], [:class => ["fruit"]], ["🍍"])
HTM.Node(["div"], [:class => ["fruit"]], Union{String, HTM.Node}["🍍"])
```
"""
struct Node
    tag::Vector{String}
    attrs::Vector{Pair{Symbol,Vector{String}}}
    children::Vector{Union{String,HTM.Node}}
    Node(tag, attrs=Pair{Symbol,Vector{String}}[], children=Union{String,HTM.Node}[]) = new(tag, attrs, children)
end
Base.:(==)(🍍::Node, 🍌::Node) = 🍍.tag == 🍌.tag && 🍍.attrs == 🍌.attrs && 🍍.children == 🍌.children

@doc raw"""
    render(x)

Generic function that defines how a Julia object is rendered.

This should normally return using [`@htm_str`](@ref).

This is an alternative to `Base.show(io::IO, m::MIME"text/html", x)` inspired
by WebIO and should only be redefined when Julia's display system is not
powerful enough for your needs.

```jldoctest
julia> struct Fruit
           name::String
           emoji::Char
       end

julia> HTM.render(🍎::Fruit) = htm"$(🍎.name): <div class=fruit>$(🍎.emoji)</div>"

julia> htm"<p>$(Fruit(\\"pineapple\\", '🍍'))</p>"
<p>pineapple: <div class="fruit">🍍</div></p>
```
"""
@inline render(🍎) = 🍎
@inline render(::Bool) = nothing

@inline create_tag(🍎) = 🍎
@inline create_tag(v::AbstractVector) = string(create_tag.(v)...)
@inline create_tag(v::AbstractVector{S}) where {S<:AbstractString} = *(create_tag.(v)...)

@inline create_value(🍎) = 🍎
@inline create_value(b::Bool) = b ? nothing : error("should have been disabled")
@inline create_value(v::AbstractVector{S}) where {S<:AbstractString} = *(v...)

@inline isenabled(🍎) = true
@inline isenabled(b::Bool) = b
@inline isenabled(::Nothing) = false

@inline create_attrs(d::AbstractDict) = create_attrs(collect(d))
@inline create_attrs(v::AbstractVector) = (attrs = Pair{Symbol,Any}[]; foreach(p -> isenabled(last(p)) && pushattr!(attrs, create_attr(p)), v); attrs)
@inline pushattr!(v::AbstractVector, p::AbstractVector) = append!(v, p)
@inline pushattr!(v::AbstractVector, p::Pair) = push!(v, p)

@inline create_attr(p::Pair) = create_attr(first(p), last(p))
@inline create_attr(🔑::Symbol, v) = create_attr(Val(🔑), v)
@inline create_attr(🔑, v) = create_attr(Symbol(🔑), v)

# attribute fallback
@inline create_attr(::Val{K}, x) where {K} = K => create_value(x)  # no interps in keys: use spread attributes

# spread attributes
@inline create_attr(::Val{Symbol()}, d) = create_attrs(d)

# style attribute
@inline create_attr(🔑::Val{:style}, d::AbstractDict) = create_attr(🔑, collect(d))
@inline create_attr(::Val{:style}, v::AbstractVector{Pair{K,V}}) where {K,V} = :style => string(cssprop.(v)...)  # TODO: remove last character?
@inline cssprop(p::Pair) = string(first(p), ':', last(p), ';')  # TODO: add space between key/value
@inline cssprop(p::Pair{Symbol,V}) where {V} = cssprop(string(first(p)) => last(p))

# class attribute
@inline create_attr(::Val{:class}, s::AbstractString) = :class => create_value(s)
@inline create_attr(::Val{:class}, m::AbstractSet) = :class => string((c -> *(c, ' ')).(m)...)  # TODO: adjust spaces around classes?
@inline create_attr(🔑::Val{:class}, v) = create_attr(🔑, Set(v))

"""
    @htm_str

Create a DOM object from a literal string.

Parsing is done via [`HTM.parse`](@ref).
"""
macro htm_str(s)
    htm = parse(s)
    esc(macrotoexpr(htm))
end

# expressions for generic contexts
@inline toexpr(s::AbstractString) = (length(s) > 1 && startswith(s, '$')) ? Meta.parse(s[nextind(s, begin):end]) : s
@inline toexpr(v::AbstractVector) = length(v) === 1 ? :($(toexpr(first(v)))) : vectoexpr(v)
@inline toexpr(p::Pair) = (v = last(p); :($(:(first($(p)))) => $(valtoexpr(last(p)))))  # no interps in keys
@inline toexpr(🍍::Node) = :($(create_element)($(tagtoexpr(🍍.tag)), $(attrstoexpr(🍍.attrs)), $(childrentoexpr(🍍.children))))

# expressions for specific contexts
@inline tagtoexpr(x) = isempty(x) ? "" : :($(create_tag)($(toexpr(x))))
@inline attrstoexpr(x) = isempty(x) ? Pair{Symbol,Any}[] : :($(create_attrs)($(vectoexpr(x))))
@inline childrentoexpr(x) = isempty(x) ? Any[] : :($(render)($(toexpr(x))))

@inline vectoexpr(v::AbstractVector) = :([$(toexpr.(v)...)])
@inline valtoexpr(v::AbstractVector) = length(v) > 1 ? toexpr(v) : (isempty(v) ? "" : toexpr(first(v)))
@inline function macrotoexpr(v::AbstractVector)
    length(v) > 1 && return toexpr(v)
    isempty(v) && return nothing
    return toexpr(first(v))
end

end


text:
"Restrict Hartree-Fock"
function hartree_fock(mol::Molecule;
        max_iter::Int = 100, δe::Float64 = 1e-10,
        large_system::Bool = false,
        return_density::Bool = false
    )

    if large_system == false
        small_system_hartree_fock(mol; max_iter = max_iter, δe = δe, return_density = return_density)
    else
        large_system_hartree_fock(mol; max_iter = max_iter, δe = δe, return_density = return_density)
    end
end

function small_system_hartree_fock(mol::Molecule;
        max_iter = 100, δe = 1e-10,
        return_density::Bool = false
    )
    basis = Vector{STO3G}()
    Ne::Int = -mol.charge
    for atom::Atom in mol.atoms
        Ne += atom.Ne
        for orbital in atom.orbitals
            push!(basis, orbital)
        end
    end
    @assert iseven(Ne) "odd electron number is not suppport yet"
    Nb::Int = length(basis)
    Nn::Int = length(mol.atoms)

    # overlap matrix
    S = zeros(Float64, Nb, Nb) + I
    for i = 1:Nb
        for j = (i+1):Nb
            S[i, j] = S[j,i] = overlap_integral(basis[i], basis[j])
        end
    end

    # kinetic matrix
    T = zeros(Float64, Nb, Nb)
    for i = 1:Nb
        for j = i:Nb
            T[i,j] = T[j,i] = kinetic_integral(basis[i], basis[j])
        end
    end

    # nuclear elergy
    nuclear_energy = 0.0
    for I = 1:Nn
        for J = (I+1):Nn
            nuclear_energy += mol.atoms[I].Z * mol.atoms[J].Z / abs(mol.atoms[I].R - mol.atoms[J].R)
        end
    end

    # nuclear attraction matrix
    Vn = zeros(Float64, Nb, Nb)
    for atom in mol.atoms
        for i = 1:Nb
            for j = i:Nb
                _v = nuclear_acctaction_integral(atom.Z, atom.R, basis[i], basis[j])
                Vn[i,j] += _v
                if i != j
                    Vn[j, i] += _v
                end
            end
        end
    end

    # one body matrix
    H0 = T + Vn

    # use libint2 for electron repulsion integral instead
    Ve = zeros(Float64, Nb, Nb, Nb, Nb)
    Ve!(Ve, basis)

    d, U = eigen(Symmetric(S))
    X = U * diagm(d.^(-1/2)) * U'

    ρ = zeros(Float64, Nb, Nb)  # density matrix
    G = zeros(Float64, Nb, Nb)  # effective interaction matrix
    F = zeros(Float64, Nb, Nb)  # Fock matrix
    H = zeros(Float64, Nb, Nb)  # orthogonal Fock matrix
    C = zeros(Float64, Nb, Nb)  # orthogonal coefficients
    Cx = zeros(Float64, Nb, Nb) # real coefficients

    electron_energy = 0.0
    total_energy = 0.0
    old_energy = Inf64

    # start hartree-fock iteration
    iter_times = 0
    for hf_iter = 1:max_iter
        @. F = H0 + G
        H .= X * F * X
        C .= eigvecs(H)
        Cx .= real(X * C)

        # update density matrix
        ρ .= 0
        for i = 1:div(Ne, 2)
            for α = 1:Nb
                for a = 1:Nb
                    @inbounds ρ[α, a] += 2*Cx[α, i] * Cx[a, i]
                end
            end
        end

        # update energy
        electron_energy = 0.0
        for α = 1:Nb
            for a = 1:Nb
                @inbounds electron_energy += ρ[a, α] * (H0[α, a] + F[α, a])
            end
        end
        electron_energy *= 0.5

        if abs(old_energy - electron_energy) < δe
            break
        end

        old_energy = electron_energy

        # update effective interaction matrix
        G .= 0
        for α = 1:Nb
            for β = 1:Nb
                for b = 1:Nb
                    for a = 1:Nb
                        G[α, a] += (Ve[α, β, b, a] - 0.5Ve[α, β, a, b]) * ρ[b, β]
                    end
                end
            end
        end
        iter_times = hf_iter
    end

    total_energy = electron_energy + nuclear_energy
    if iter_times == max_iter
        println("Warning: hartree_fock iterator time reaches `max_iter`, but may not reach the convergence condition")
    end
    if return_density
        return total_energy, ρ
    else
        return total_energy
    end
end

function large_system_hartree_fock(mol::Molecule;
        max_iter = 100, δe = 1e-10,
        return_density = false
    )
    basis = Vector{STO3G}()
    Ne::Int = -mol.charge
    for atom::Atom in mol.atoms
        Ne += atom.Ne
        for orbital in atom.orbitals
            push!(basis, orbital)
        end
    end
    @assert iseven(Ne) "odd electron number is not suppport yet"
    Nb::Int = length(basis)
    Nn::Int = length(mol.atoms)

    # overlap matrix
    S = zeros(Float64, Nb, Nb) + I
    for i = 1:Nb
        for j = (i+1):Nb
            S[i, j] = S[j,i] = overlap_integral(basis[i], basis[j])
        end
    end

    # kinetic matrix
    T = zeros(Float64, Nb, Nb)
    for i = 1:Nb
        for j = i:Nb
            T[i,j] = T[j,i] = kinetic_integral(basis[i], basis[j])
        end
    end

    # nuclear elergy
    nuclear_energy = 0.0
    for I = 1:Nn
        for J = (I+1):Nn
            nuclear_energy += mol.atoms[I].Z * mol.atoms[J].Z / abs(mol.atoms[I].R - mol.atoms[J].R)
        end
    end

    # nuclear attraction matrix
    Vn = zeros(Float64, Nb, Nb)
    for atom in mol.atoms
        for i = 1:Nb
            for j = i:Nb
                _v = nuclear_acctaction_integral(atom.Z, atom.R, basis[i], basis[j])
                Vn[i,j] += _v
                if i != j
                    Vn[j, i] += _v
                end
            end
        end
    end

    # one body matrix
    H0 = T + Vn

    d, U = eigen(Symmetric(S))
    X = U * diagm(d.^(-1/2)) * U'

    ρ = zeros(Float64, Nb, Nb)  # density matrix
    G = zeros(Float64, Nb, Nb)  # effective interaction matrix
    F = zeros(Float64, Nb, Nb)  # Fock matrix
    H = zeros(Float64, Nb, Nb)  # orthogonal Fock matrix
    C = zeros(Float64, Nb, Nb)  # orthogonal coefficients
    Cx = zeros(Float64, Nb, Nb) # real coefficients

    # prepare data for Veff! function, which links to libint2
    Nshell, l, alpha, coeff, R = make_shells(basis)

    electron_energy = 0.0
    total_energy = 0.0
    old_energy = Inf64

    # start hartree-fock iteration
    iter_times = 0
    for hf_iter = 1:max_iter
        Hx = H0 + G
        H = X * Hx * X
        ε, C = eigen(H)
        Cx = real(X * C)

        # update density matrix
        ρ .= 0
        for i = 1:div(Ne, 2)
            for α = 1:Nb
                for a = 1:Nb
                    ρ[α, a] += 2*Cx[α, i] * Cx[a, i]
                end
            end
        end

        # update energy
        electron_energy = 0.0
        for α = 1:Nb
            for a = 1:Nb
                electron_energy += ρ[a, α] * (H0[α, a] + Hx[α, a])
            end
        end
        electron_energy *= 0.5

        if abs(old_energy - electron_energy) < δe
            break
        end

        old_energy = electron_energy

        # update effective interaction matrix
        Veff!(G, ρ, Nshell, l, alpha, coeff, R)
        iter_times = hf_iter
    end

    total_energy = electron_energy + nuclear_energy
    if iter_times == max_iter
        println("Warning: hartree_fock iterator time reaches `max_iter`, but may not reach the convergence condition")
    end
    if return_density
        return total_energy, ρ
    else
        return total_energy
    end
end


text:
@testset "Scheme Types" begin
    @testset "Parse" begin
        der = hex2bytes("304402201f62993ee03fca342fcb45929993fa6ee885e00ddad8de154f268d98f083991402201e1ca12ad140c04e0e022c38f7ce31da426b8009d02832f0b44f39a6b178b7a1")
        sig = Signature{:ECDSA}(big"0x1f62993ee03fca342fcb45929993fa6ee885e00ddad8de154f268d98f0839914",
                                big"0x1e1ca12ad140c04e0e022c38f7ce31da426b8009d02832f0b44f39a6b178b7a1")
        @test Secp256k1.Signature(der) == sig
    end
    @testset "Serialize" begin
        testcases = (
            (1, 2),
            (rand(big.(0:big(2)^255)), rand(big.(0:big(2)^255))),
            (rand(big.(0:big(2)^255)), rand(big.(0:big(2)^255))),
            (rand(big.(0:big(2)^255)), rand(big.(0:big(2)^223))))
        for x in testcases
            sig  = Signature{:ECDSA}(x[1], x[2])
            der  = Secp256k1.serialize(sig)
            sig2 = Secp256k1.Signature(der)
            @test sig2 == sig
        end
    end
end


text:
#
# Regressions for take-up of disability benefits.
# using just what's in the model datasets
#
# include("intro.jl")
# using STBRegressions

using ScottishTaxBenefitModel
using .Definitions
using .RunSettings: Settings
using CSV,DataFrames,GLM,RegressionTables

const settings = Settings()

#
# UK wide with Scottish dummy
#
frshh = CSV.File("$(MODEL_DATA_DIR)/model_households.tab" ) |> DataFrame
frspeople = CSV.File("$(MODEL_DATA_DIR)/model_people.tab") |> DataFrame

fm = innerjoin( frshh, frspeople, on=[:data_year, :hid ], makeunique=true )

fm = innerjoin( frshh, frspeople, on=[:data_year, :hid ], makeunique=true )
fm.age_sq = fm.age.^2
fm.deaf_blind=fm.registered_blind .| fm.registered_deaf .| fm.registered_partially_sighted
fm.yr = fm.data_year .- 2014
fm.any_dis = (
    fm.disability_vision .|
    fm.disability_hearing .|
    fm.disability_mobility .|
    fm.disability_dexterity .|
    fm.disability_learning .|
    fm.disability_memory .|
    fm.disability_other_difficulty .|
    fm.disability_mental_health .|
    fm.disability_stamina .|
    fm.disability_socially )
fm.adls_bad=fm.adls_are_reduced.==1
fm.adls_mid=fm.adls_are_reduced.==2
fm.rec_carers = fm.income_carers_allowance .> 0
fm.rec_dla = ( fm.income_dlamobility.>0.0) .| ( fm.income_dlaself_care .>0.0 )
fm.rec_dla_child = ( fm.rec_dla ) .& (fm.age .<= 16 )
fm.rec_dla_adult = ( fm.rec_dla ) .& (fm.age .> 16 )

fm.rec_dla_care = ( fm.income_dlaself_care .>0.0 )
fm.rec_dla_mob = ( fm.income_dlamobility.>0.0 )
fm.rec_pip = ( fm.income_personal_independence_payment_mobility.>0.0) .| ( fm.income_personal_independence_payment_daily_living .>0.0 )
fm.rec_pip_care = ( fm.income_personal_independence_payment_daily_living .>0.0 )
fm.rec_pip_mob = ( fm.income_personal_independence_payment_mobility.>0.0)
fm.rec_esa = ( fm.income_employment_and_support_allowance.>0.0)
fm.rec_aa = ( fm.income_attendance_allowance.>0.0)
fm.rec_carers = ( fm.income_carers_allowance.>0.0)
fm_rec_aa = ( fm.income_attendance_allowance.>0.0)
fm.scotland = fm.region .== 299999999
fm.male = fm.sex .== 1

# subsets
fm_working_age = fm[((fm.age .< 65).&(fm.age .> 16 )),:]
fm_pension_age = fm[(fm.age .>= 65),:]
fm_children = fm[(fm.age .<= 16),:]

#
# Regress on everything.
#
main_out = []

bens = [
  # :rec_dla_child, # in children - too many regressors
  :rec_carers,
  :rec_pip,
  :rec_pip_mob,
  :rec_pip_care,
  :rec_esa,
  :rec_dla_adult, # in adults
  :rec_dla_mob,
  :rec_dla_care,
  :rec_aa ]

  datasets = Dict(
    :rec_dla_child=>fm_children,
    :rec_carers=>fm,
    :rec_pip=>fm_working_age,
    :rec_pip_mob=>fm_working_age,
    :rec_pip_care=>fm_working_age,
    :rec_esa=>fm_working_age,
    :rec_dla_adult=>fm_working_age,
    :rec_dla_mob=>fm_working_age,
    :rec_dla_care=>fm_working_age,
    :rec_aa=>fm_pension_age )

    main_out=[]

for target in bens
      v = glm(
        Term( target ) ~
          Term( :scotland ) +
          Term( :data_year )+
          Term( :male )+
          Term( :age )+
          Term( :age_sq ) +
          Term( :has_long_standing_illness )+
          Term( :adls_bad )+
          Term( :adls_mid )+
          Term( :registered_blind )+
          Term( :registered_partially_sighted )+
          Term( :registered_deaf )+

          Term( :disability_dexterity )+
          Term( :disability_learning )+
          Term( :disability_memory )+
          Term( :disability_mental_health )+
          Term( :disability_mobility )+
          Term( :disability_socially )+
          Term( :disability_stamina )+
          Term( :disability_other_difficulty ),
        datasets[target ],
        Binomial(),
        ProbitLink(), contrasts=Dict( :data_year=>DummyCoding()))
    push!( main_out, v )
end

RegressionTables.regtable( main_out ...;  renderSettings = latexOutput("disability_regressions.tex"))
RegressionTables.regtable( main_out ... )

##
# significant ones only
#
# just over
rec_carers = glm( @formula( rec_carers ~ male+age+age_sq+adls_bad+adls_mid+registered_deaf+disability_mental_health), fm, Binomial(), ProbitLink() )
rec_pip =  glm( @formula( rec_pip ~ scotland+data_year+has_long_standing_illness+adls_bad+adls_mid+disability_memory+disability_mental_health+disability_mobility ), fm_working_age, Binomial(), ProbitLink(), contrasts=Dict( :data_year=>DummyCoding()) )
rec_pip_care =  glm( @formula( rec_pip_care ~ scotland+data_year+has_long_standing_illness+adls_bad+adls_mid+disability_memory+disability_mental_health+disability_mobility ), fm_working_age, Binomial(), ProbitLink(), contrasts=Dict( :data_year=>DummyCoding()) )
rec_pip_mob =  glm( @formula( rec_pip_mob ~ scotland+data_year+registered_blind+has_long_standing_illness+adls_bad+adls_mid+disability_memory+disability_mental_health+disability_mobility ), fm_working_age, Binomial(), ProbitLink(), contrasts=Dict( :data_year=>DummyCoding()) )
rec_aa = glm( @formula( rec_aa ~ male+adls_bad+adls_mid+registered_blind+disability_dexterity+disability_mental_health+disability_mobility ), fm_pension_age, Binomial(), ProbitLink())

#
# v. simple child one that combines mob and care
#
rec_dla_child = glm( @formula( rec_dla_child ~ male+adls_bad+disability_mental_health ), fm_children, Binomial(), ProbitLink() )

RegressionTables.regtable( rec_carers, rec_pip, rec_pip_care, rec_pip_mob, rec_aa, rec_dla_child )

fm_working_age.pip_care_prob = predict( rec_pip_care )
fm_working_age.pip_mob_prob = predict( rec_pip_mob )
fm_pension_age.aa_prob = predict( rec_aa )
fm_children.dla_prob = predict( rec_dla_child )

#
# Everybody *Not* receiving pip_care or DLA, ranked by probability of receiving it
#
positive_candidates_pip_care = fm_working_age[(((fm_working_age.rec_pip_care .== false) .& (fm_working_age.rec_dla_care .== false)) .& (fm_working_age.scotland .== true)), [:data_year,:hid,:pid,:weight, :pip_care_prob,:rec_pip_care,:rec_dla_care,:scotland]]
sort!(positive_candidates_pip_care, :pip_care_prob, rev=true )
CSV.write( "data/disability/positive_candidates_pip_care.csv", positive_candidates_pip_care )

#
# *least* likely people to be on PIP Care who are on it
#
negative_candidates_pip_care = fm_working_age[((fm_working_age.rec_pip_care .== true) .& (fm_working_age.scotland .== true)), [:data_year,:hid,:pid,:weight, :pip_care_prob,:rec_pip_care,:rec_dla_care,:scotland]]
sort!(negative_candidates_pip_care, :pip_care_prob  )
CSV.write( "data/disability/negative_candidates_pip_care.csv", negative_candidates_pip_care )

#
# Everybody *Not* receiving pip_mobility or DLA, ranked by probability of receiving it. For a more generous test, we go down this list.
#
positive_candidates_pip_mob = fm_working_age[((fm_working_age.rec_pip_mob .== false) .& (fm_working_age.rec_dla_mob .== false) .& (fm_working_age.scotland .== true)), [:data_year,:hid,:pid,:weight, :pip_mob_prob,:rec_pip_mob,:rec_dla_care,:scotland]]
sort!(positive_candidates_pip_mob, :pip_mob_prob, rev=true )
CSV.write( "data/disability/positive_candidates_pip_mob.csv", positive_candidates_pip_mob )

#
# *least* likely people to be on PIP Mobility who are actually on it. If we want a less generous test, we go down this list.
#
negative_candidates_pip_mob = fm_working_age[((fm_working_age.rec_pip_mob .== true) .& (fm_working_age.scotland .== true)), [:data_year,:hid,:pid,:weight, :pip_mob_prob,:rec_pip_mob,:rec_dla_care,:scotland]]
sort!(negative_candidates_pip_mob, :pip_mob_prob )
CSV.write( "data/disability/negative_candidates_pip_mob.csv", negative_candidates_pip_mob )

#
# ditto for AA amongst pension age
#
positive_candidates_aa = fm_pension_age[((fm_pension_age.rec_aa .== false) .& (fm_pension_age.scotland .== true)), [:data_year,:hid,:pid,:weight, :aa_prob,:rec_aa,:scotland]]
sort!(positive_candidates_aa, :aa_prob, rev=true )
CSV.write( "data/disability/positive_candidates_aa.csv", positive_candidates_aa )
#
# *least* likely people to be on AA who are actually on it. If we want a less generous test, we go down this list.
#
negative_candidates_aa = fm_pension_age[((fm_pension_age.rec_aa .== true) .& (fm_pension_age.scotland .== true)), [:data_year,:hid,:pid,:weight, :aa_prob,:rec_aa,:scotland]]
sort!(negative_candidates_aa, :aa_prob )

CSV.write( "data/disability/negative_candidates_aa.csv", negative_candidates_aa )

#
# ditto for DLA amongst children
#
positive_candidates_aa = fm_pension_age[((fm_pension_age.rec_aa .== false) .& (fm_pension_age.scotland .== true)), [:data_year,:hid,:pid,:weight, :aa_prob,:rec_aa,:scotland]]
sort!(positive_candidates_aa, :aa_prob, rev=true )
CSV.write( "data/disability/positive_candidates_aa.csv", positive_candidates_aa )
#
# *least* likely people to be on AA who are actually on it. If we want a less generous test, we go down this list.
#
negative_candidates_aa = fm_pension_age[((fm_pension_age.rec_aa .== true) .& (fm_pension_age.scotland .== true)), [:data_year,:hid,:pid,:weight, :aa_prob,:rec_aa,:scotland]]
sort!(negative_candidates_aa, :aa_prob )

CSV.write( "data/disability/negative_candidates_aa.csv", negative_candidates_aa )

positive_candidates_dla_children = fm_children[((fm_children.rec_dla .== false) .& (fm_children.scotland .== true)), [:data_year,:hid,:pid,:weight, :dla_prob,:rec_dla,:scotland]]
sort!(positive_candidates_dla_children, :dla_prob, rev=true )
CSV.write( "data/disability/positive_candidates_dla_children.csv", positive_candidates_dla_children )
#
# *least* likely people to be on AA who are actually on it. If we want a less generous test, we go down this list.
#
negative_candidates_dla_children = fm_children[((fm_children.rec_dla .== true) .& (fm_children.scotland .== true)), [:data_year,:hid,:pid,:weight, :dla_prob,:rec_dla,:scotland]]
sort!(negative_candidates_dla_children, :dla_prob )
CSV.write( "data/disability/negative_candidates_dla_children.csv", negative_candidates_dla_children )

max_stars_repo_name: grahamstark/ScottishTaxBenefitModel

max_issues_repo_name: grahamstark/ScottishTaxBenefitModel

max_forks_repo_name: grahamstark/ScottishTaxBenefitModel


text:
# p25.jl - stability regions for ODE formulas

# Adams-Bashforth:
clf(); subplot(221);
zplot(z) = plot(real(z),imag(z));
plot([-8,8],[0,0],"k-"); plot([0,0],[-8,8],"k-");
z = exp.(1im*pi*(0:200)/100); r = z-1;
s = 1; zplot(r./s);                                  # order 1
s = (3-1./z)/2; zplot(r./s);                         # order 2
s = (23-16./z+5./z.^2)/12; zplot(r./s);             # order 3
axis("square"); axis([-2.5,.5,-1.5,1.5]); grid(true)
title("Adams—Bashforth");

# Adams-Moulton:
subplot(222);
plot([-8,8],[0,0],"k-"); plot([0,0],[-8,8],"k-");
s = (5*z+8-1./z)/12; zplot(r./s);                    # order 3
s = (9*z+19-5./z+1./z.^2)/24; zplot(r./s);           # order 4
s = (251*z+646-264./z+106./z.^2-19./z.^3)/720; zplot(r./s) ;   # 5
d = 1-1./z;
s = 1-d/2-d.^2/12-d.^3/24-19*d.^4/720-3*d.^5/160; zplot(d./s); # 6
axis("square"); axis([-7,1,-4,4]); grid(true);
title("Adams—Moulton");

# Backward differentiation:
subplot(223);
plot([-40,40],[0,0],"k-"); plot([0,0],[-40,40],"k-");
r = 0;
for i = 1:6  # orders 1-6
    r += (d.^i)/i;
    zplot(r);
end
axis("square"); axis([-15,35,-25,25]); grid(true);
title("backward differentiation")

# Runge-Kutta:
subplot(224);
plot([-8,8],[0,0],"k-"); plot([0,0],[-8,8],"k-");
w = zeros(4)+0im; W = w.';
for i = 2:length(z)
    # orders 1-4
    w[1] -= (1+w[1]-z[i]);
    w[2] -= (1+w[2]+.5*w[2]^2-z[i]^2)/(1+w[2]);
    w[3] -= (1+w[3]+.5*w[3]^2+w[3]^3/6-z[i]^3)/(1+w[3]+w[3]^2/2);
    w[4] -= (1+w[4]+.5*w[4]^2+w[4]^3/6+w[4].^4/24-z[i]^4)/(1+w[4]+w[4]^2/2+w[4]^3/6);
    W = [W; w.'];
end
zplot(W); axis("square"); axis([-5,2,-3.5,3.5]); grid(true);
title("Runge—Kutta");


text:
@enum Conversion::UInt8 CONV_UNSPECIFIED CONV_REPR CONV_STRING

# replacement field
struct Field
    argument::Argument
    conv::Conversion
    spec::Spec
end


text:
"""
    getinitialvalue(::Type{GeneralizedExtremeValue},y::Vector{<:Real})::Vector{<:Real}

Compute the initial values of the GEV parameters given the data `y`.

If the probability weighted moment estimations are valid, then those values are
returned. Otherwise, the probability weighted moment estimations of the Gumbel
distribution are returned.

# Example
```julia-repl
julia> Extremes.getinitialvalue(GeneralizedExtremeValue, y)
```
"""
function getinitialvalue(::Type{GeneralizedExtremeValue},y::Vector{<:Real})::Vector{<:Real}

    # Fit the model with by the probability weigthed moments
    fm = gevfitpwm(y)

    # Convert to fitted model in a Distribution object
    fd = getdistribution(fm.model, fm.θ̂)[]

    # check if initial values are in the domain of the GEV
    valid_initialvalues = all(insupport(fd,y))

    #= If one at least one value does not lie in the support, then the initial
     values are replaced by the Gumbel initial values. =#
    if valid_initialvalues
        μ₀ = location(fd)
        ϕ₀ = log(scale(fd))
        ξ₀ = Distributions.shape(fd)
    else
        fm = gumbelfitpwm(y)
        μ₀ = fm.θ̂[1]
        ϕ₀ = fm.θ̂[2]
        ξ₀ = 0.0
    end

    initialvalues = [μ₀, ϕ₀, ξ₀]

    return initialvalues

end

"""
    getinitialvalue(::Type{GeneralizedPareto},y::Vector{<:Real})::Vector{<:Real}

Compute the initial values of the GP distribution parameters given the data `y`.

If the probability weighted moment estimations are valid, then those values are
returned. Otherwise, the moment estimation of the exponential distribution is returned.

# Example
```julia-repl
julia> Extremes.getinitialvalue(GeneralizedPareto, y)
```
"""
function getinitialvalue(::Type{GeneralizedPareto},y::Vector{<:Real})::Vector{<:Real}

    # Fit the model with by the probability weigthed moments
    fm = gpfitpwm(y)

    # Convert to fitted model in a Distribution object
    fd = getdistribution(fm.model, fm.θ̂)[]

    if all(insupport(fd,y))
        σ₀ = scale(fd)
        ξ₀ = Distributions.shape(fd)
    else
        σ₀ = mean(y)
        ξ₀ = 0.0
    end

    initialvalues = [log(σ₀), ξ₀]

    return initialvalues

end

"""
    slicematrix(A::AbstractMatrix{T}; dims::Int=1)::Array{Array{T,1},1} where T

Convert a Matrix in a Vector of Vectors.

The slicing dimension is defined with the optional `dims` keyword.

# Examples
```julia-repl
julia> A = rand(1:10,5,4)
julia> V₁ = Extremes.slicematrix(A)
julia> V₂ = Extremes.slicematrix(A, dims=2)
```
"""
function slicematrix(A::AbstractMatrix{T}; dims::Int=1)::Array{Array{T,1},1} where T

    @assert dims∈(1,2) "dims should be either 1 or 2."

    n, m = size(A)

    if dims==1
        B = Vector{T}[Vector{T}(undef, n) for _ in 1:m]
        for i in 1:m
            B[i] .= A[:, i]
        end
    else
        B = Vector{T}[Vector{T}(undef, m) for _ in 1:n]
        for i in 1:n
            B[i] .= A[i, :]
        end
    end

    return B

end

"""
    unslicematrix(B::Array{Array{T,1},1}; dims::Int=1)::AbstractMatrix{T} where T

Convert a Vector of Vectors in a Matrix if all vectors share the same length.

The concatenation dimension is defined with the optional `dims` keyword.

# Examples
```julia-repl
julia> V = Extremes.slicematrix(rand(1:10, 5, 4))
julia> A₁ = Extremes.unslicematrix(V)
julia> A₂ = Extremes.unslicematrix(V, dims=2)
```
"""
function unslicematrix(B::Array{Array{T,1},1}; dims::Int=1)::AbstractMatrix{T} where T

    @assert dims∈(1,2) "dims should be either 1 or 2."

    m = length.(B)

    @assert all( ==(m[1]), m) "all vectors must have the same length."

    n = length(B)

    m = m[1]

    if dims == 1

        A = Matrix{T}(undef, m, n)

        for i in 1:n
            A[:,i] = B[i]
        end

    else

        A = Matrix{T}(undef, n, m)

        for i in 1:n
            A[i,:] = B[i]
        end

    end

    return A

end

"""
    buildVariables(df::DataFrame, ids::Vector{Symbol})::Vector{Variable}

Build the `Variable` type from the columns `ids` of the DataFrame `df`.

# Example
```julia-repl
julia> df = Extremes.dataset("fremantle")
julia> Extremes.buildVariables(df, [:Year, :SOI])
```
"""
function buildVariables(df::DataFrame, ids::Vector{Symbol})::Vector{Variable}

    variables = Vector{Variable}()

    for id in ids
        push!(variables, Variable(string(id), df[:,id]))
    end

    return variables

end


text:
using H3, H3Geometry

"""
Finding the history for a point involves computing the fine-grained
H3 reference (typically at level 15) as well as the H3 tile (typically
level 3). The tile is used to form the file name from which to
read the data and the reference is used to filter the data read from
the file. There may be multiple files to read, each with the same
name, but in different directories.
"""
function queryData(latitude, longitude, dirs; refLevel=15, tileLevel=3)
    ref = H3.API.geoToH3(latitude, longitude, refLevel)
    queryData(ref, dirs; tileLevel=tileLevel)
end

function queryData(refs::Vector, dirs::Vector{String}; tileLevel=3)
    # the refs will almost always be on a single tile
    tiles = unique(H3.API.h3ToParent.(refs, tileLevel))

    df = Vector()
    for tile in tiles
        fname = "x-$tileLevel-$tile"
        for d in dirs
            fx = Parquet.Table(joinpath(d, fname))
            rx = DataFrame(fx) |> @filter(_.h3 in Int64.(refs)) |> DataFrame
            push!(df, rx)
        end
    end
    r = reduce(vcat, df, init=DataFrame())
    r[!, :time] = Dates.unix2datetime.(r.t ./ 1000)
#    r[!, :lat] .= latitude
#    r[!, :lon] .= longitude
    @info "result" size(r)
    r
end

"""
Reads the index that lets us find the nearest data point from a file.

The structure of the index is that of a dictionary keyed by h3 indexes
at resolution 5 and 6. The value is a list of h3 indexes at level
15 which have been observed in the data files and for which, presumably,
we can find data.
"""
function readIndex(fname)
    Dict(eachrow(DataFrame(Arrow.Table(fname))))
end

"""
Finds the closest grid points in an index to a specified lat/long pair

All grid points which are an equal grid distance from the original point
will be returned. Since we are using grid distance, there will likely
be more than one result and each should be used in querying.
"""
function getGridPoints(latitude, longitude, index)
    ref = H3.API.geoToH3(latitude, longitude, 15)
    candidates = index[H3.API.h3ToParent(ref, 6)]
    distances = H3.API.h3Distance.(ref, candidates)
    candidates[distances .== minimum(distances)]
end

"""
Finds up-leveled hexes that overlap with a ring around the specified point.

This is slightly tricky because the parent of a hex doesn't necessarily contain
the center of the hex.
"""
function uplevel(ref::UInt64; level, radius=2)
    ref |> x->kRing(x, radius) .|> h3ToGeo .|> x -> geoToH3(x, level)
end

"""
Finds up-leveled hexes that overlap with a ring around the specified point.
"""
function uplevel(ref::Int64; level, radius=2)
    uplevel(UInt64(ref), level; radius=radius)
end

"""
Finds up-leveled hexes that overlap with a ring around the specified points.
"""
function uplevel(refs::Vector; level::Int, radius=2)
    unique(reduce(vcat, uplevel.(refs, level, radius=radius)))
end

"""
Finds the tiles that will hold the nearest grid points. Since we
know that the nearest grid point will be well less than 1000 grid steps
away, we could just form the set of h3 references within that distance
(called a k-ring) and figure out the set of tiles associated with any
of those points. That is a lot of computation, however. Instead, we
find the equivalent neighborhood at a coarser level of resolution and
then find the tiles for that neighborhood.

As a point of reference, the radius 1000 k-ring for a point has about
2-3 million elements while the level 10 k-ring with radius 7 covers
very nearly the same area with fewer than 200 hexes and these hexes
can be used to find the interesting tiles far faster (84us versus 2.6
seconds) than the 2-3 million with almost precisely the same result.
"""
function findTiles(ref; radius=7, level=10, tileLevel=3)
    locality = uplevel(ref, 10, radius=2)
    uplevel(locality, 3, radius=10)
end

"""
Finds the closest grid points in an index to a specified lat/long pair

All grid points which are an equal grid distance from the original point
will be returned. Since we are using grid distance, there will likely
be more than one result and each should be used in querying.
"""
function getGridPoints(ref; level=6, basedir)
    upref = uplevel(ref, 8, radius=0)[1]
    r = []
    tiles = findTiles(ref)
    for tile in tiles
        fname = joinpath(basedir, "x-3-$tile")
        df = DataFrame(Parquet.Table(fname)) |>
            @mutate(dx = h3Distance(upref, uplevel(UInt64(get(_.h3)), 8, radius=0)[1])) |>
            @filter(10_000 > _.dx >= 0) |>
            @take(10) |>
            @mutate(tile=findTiles(get(_.h3))) |>
            DataFrame
        push!(r, df)
    end
    df = vcat(r...)
    df = df |> @orderby(_.dx) |> DataFrame
    df[!, :origin] .= ref
    df[!, :rank] = 1:size(df)[1]
    df
end

function nearRegion(ref::UInt64, radius::Int, level::Int)
    reduce(vcat, kRing.(reduce(vcat, h3ToParent.(kRing(ref, 2), level)), radius))
end

max_stars_repo_name: agstack/weather-server

max_issues_repo_name: agstack/weather-server

max_forks_repo_name: agstack/weather-server


text:
iscapitalletter(b) = UInt(0x0041) <= b <= UInt(0x05A)
islowercaseletter(b) = UInt(0x0061) <= b <= UInt(0x07A)
isnamestart(b) = iscapitalletter(b) || islowercaseletter(b) || b == UInt('_')
isnamecontinue(b) = isnamestart(b) || isdigit(Char(b))

isunicodechar(b) = iscapitalletter(b) || islowercaseletter(b) || isdigit(Char(b))

isstringcontinue(b) = issourcechar(b) && !islineterminator(b)

iscommentstart(b) = b == UInt('#')
iscommentcontinue(b) = issourcechar(b) && !islineterminator(b)
issourcechar(b) = iswhitespace(b) || islineterminator(b) || UInt(0x0020) <= b <= UInt(0xFFFF)

iswhitespace(b) = b == UInt8('\t') || b == UInt8(' ')
islineterminator(b) = b == UInt8('\n') || b == UInt8('\r')
isunicodeBOM(b) = b == UInt(0xFEFF)
isnoncommentignored(b) = iswhitespace(b) || islineterminator(b) || isunicodeBOM(b) || b == UInt(',')

isexponent(b) = b == UInt('E') || b == UInt('e')

# Approach from JSON3.jl
macro eof()
    esc(quote
        if pos > len
            invalid("Unexpected End of File", buf, pos)
        end
    end)
end

macro skip_ignored()
    esc(quote
        b = getbyte(buf, pos)
        while isnoncommentignored(b) || iscommentstart(b)
            if iscommentstart(b)
                pos += 1
                line += 1
                @eof
                b = getbyte(buf, pos)
                while iscommentcontinue(b)
                    pos += 1
                    line += 1
                    @eof
                    b = getbyte(buf, pos)
                end
            else
                if islineterminator(b)
                    line += 1
                    column = 1
                else
                    column += 1
                end
                pos += 1
                @eof
                b = getbyte(buf, pos)
            end
        end
    end)
end

macro eof_skip_ignored()
    esc(quote
        @eof
        @skip_ignored
    end)
end

function only_ignored_left(buf, pos, len)
    while pos <= len
        b = getbyte(buf, pos)
        if isnoncommentignored(b)
            pos += 1
            continue
        end
        if iscommentstart(b)
            while pos <= len && iscommentcontinue(b)
                b = getbyte(buf, pos)
                pos += 1
            end
            continue
        end
        return false
    end
    return true
end

@noinline function invalid(error_text, buf, pos)
    n_chars = 25
    len = length(buf)
    start_char = pos - n_chars
    end_char = pos + n_chars
    indent = n_chars + 1
    if len - pos < n_chars
        end_char = length(buf)
    end
    if pos <= n_chars
        start_char = 1
        indent = pos
    end

    escaped_str = escape_string(String(@view(buf[start_char:end_char])))
    display_text_up_to_indent = @view(buf[start_char:min(end, indent + start_char - 1)])
    indent += length(escape_string(String(display_text_up_to_indent))) - length(display_text_up_to_indent)

    throw(
        ArgumentError(
            """
            invalid GraphQL string at byte position $pos while parsing
                $error_text
                $escaped_str
                $(lpad("^", indent))
            """
        )
    )
end

function getbyte(buf::AbstractVector{UInt8}, pos)
    @inbounds b = buf[pos]
    return b
end