text:
#=
The code below was derived from the scipy implementation of the LOBPCG algorithm in https://github.com/scipy/scipy/blob/v1.1.0/scipy/sparse/linalg/eigen/lobpcg/lobpcg.py#L109-L568.

Since the link above mentions the license is BSD license, the notice for the BSD license 2.0 is hereby given below giving credit to the authors of the Python implementation.

Copyright (c) 2018, Robert Cimrman, Andrew Knyazev
All rights reserved.

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:
    * Redistributions of source code must retain the above copyright
      notice, this list of conditions and the following disclaimer.
    * Redistributions in binary form must reproduce the above copyright
      notice, this list of conditions and the following disclaimer in the
      documentation and/or other materials provided with the distribution.
    * The names of the contributors may not be used to endorse or promote products
      derived from this software without specific prior written permission.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDERS BE LIABLE FOR ANY
DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
=#

export lobpcg, lobpcg!, LOBPCGIterator

using LinearAlgebra
using Random

struct LOBPCGState{TR,TL}
    iteration::Int
    residual_norms::TR
    ritz_values::TL
end
function Base.show(io::IO, t::LOBPCGState)
    @printf io "%8d    %14e\n" t.iteration maximum(t.residual_norms)
    return
end

const LOBPCGTrace{TR,TL} = Vector{LOBPCGState{TR,TL}}
function Base.show(io::IO, tr::LOBPCGTrace)
    @printf io "Iteration    Maximum residual norm \n"
    @printf io "---------    ---------------------\n"
    for state in tr
        show(io, state)
    end
    return
end

struct LOBPCGResults{TL, TX, T, TR, TI, TM, TB, TTrace <: Union{LOBPCGTrace, AbstractVector{<:LOBPCGTrace}}}
    λ::TL
    X::TX
    tolerance::T
    residual_norms::TR
    iterations::TI
    maxiter::TM
    converged::TB
    trace::TTrace
end
function EmptyLOBPCGResults(X::TX, k::Integer, tolerance, maxiter) where {T, TX<:AbstractMatrix{T}}
    blocksize = size(X,2)
    λ = Vector{T}(undef, k)
    X = TX(undef, size(X, 1), k)

    iterations = zeros(Int, ceil(Int, k/blocksize))
    residual_norms = copy(λ)
    converged = falses(k)
    trace = fill(LOBPCGTrace{Vector{real(T)},Vector{T}}(), k÷blocksize+1)

    return LOBPCGResults(λ, X, tolerance, residual_norms, iterations, maxiter, converged, trace)
end

function Base.append!(r1::LOBPCGResults, r2::LOBPCGResults, n1, n2=length(r2.λ))
    n = n1 + n2
    r1.λ[n1+1:n] .= @view r2.λ[end-n2+1:end]
    r1.residual_norms[n1+1:n] .= @view r2.residual_norms[end-n2+1:end]
    r1.X[:, n1+1:n] .= @view r2.X[:, end-n2+1:end]

    ind = n1 ÷ length(r2.λ) + 1
    r1.iterations[ind] = r2.iterations
    r1.converged[n1+1:n] .= r2.converged
    r1.trace[ind] = r2.trace

    return r1
end
function Base.show(io::IO, r::LOBPCGResults)
    first_two(fr) = [x for (i, x) in enumerate(fr)][1:2]

    @printf io "Results of LOBPCG Algorithm\n"
    @printf io " * Algorithm: LOBPCG - CholQR\n"

    if length(join(r.λ, ",")) < 40 || length(r.λ) <= 2
        @printf io " * λ: [%s]\n" join(r.λ, ",")
    else
        @printf io " * λ: [%s, ...]\n" join(first_two(r.λ), ",")
    end

    if length(join(r.residual_norms, ",")) < 40 || length(r.residual_norms) <= 2
        @printf io " * Residual norm(s): [%s]\n" join(r.residual_norms, ",")
    else
        @printf io " * Residual norm(s): [%s, ...]\n" join(first_two(r.residual_norms), ",")
    end
    @printf io " * Convergence\n"
    @printf io "   * Iterations: %s\n" r.iterations
    @printf io "   * Converged: %s\n" all(r.converged)
    @printf io "   * Iterations limit: %s\n" r.maxiter

    return
end

struct Blocks{Generalized, T, TA<:AbstractArray{T}}
    block::TA # X, R or P
    A_block::TA # AX, AR or AP
    B_block::TA # BX, BR or BP
end
Blocks(X, AX) = Blocks{false, eltype(X), typeof(X)}(X, AX, X)
Blocks(X, AX, BX) = Blocks{true, eltype(X), typeof(X)}(X, AX, BX)
function A_mul_X!(b::Blocks, A)
    mul!(b.A_block, A, b.block)
    return
end
function A_mul_X!(b::Blocks, A, n)
    mul!(view(b.A_block, :, 1:n), A, view(b.block, :, 1:n))
    return
end
function B_mul_X!(b::Blocks{true}, B)
    mul!(b.B_block, B, b.block)
    return
end
function B_mul_X!(b::Blocks{true}, B, n)
    mul!(view(b.B_block, :, 1:n), B, view(b.block, :, 1:n))
    return
end
function B_mul_X!(b::Blocks{false}, B, n = 0)
    return
end

mutable struct Constraint{T, TVorM<:Union{AbstractVector{T}, AbstractMatrix{T}}, TM<:AbstractMatrix{T}, TC}
    Y::TVorM
    BY::TVorM
    gram_chol::TC
    gramYBV::TM # to be used in view
    tmp::TM # to be used in view
end
struct BWrapper end
struct NotBWrapper end

Constraint(::Nothing, B, X) = Constraint(nothing, B, X, BWrapper())
Constraint(::Nothing, B, X, ::NotBWrapper) = Constraint(nothing, B, X, BWrapper())
function Constraint(::Nothing, B, X, ::BWrapper)
    return Constraint{Nothing, Matrix{Nothing}, Matrix{Nothing}, Nothing}(Matrix{Nothing}(undef,0,0), Matrix{Nothing}(undef,0,0), nothing, Matrix{Nothing}(undef,0,0), Matrix{Nothing}(undef,0,0))
end

Constraint(Y, B, X) = Constraint(Y, B, X, BWrapper())
function Constraint(Y, B, X, ::BWrapper)
    T = eltype(X)
    if B isa Nothing
        BY = Y
    else
        BY = similar(Y)
        mul!(BY, B, Y)
    end
    return Constraint(Y, BY, X, NotBWrapper())
end
function Constraint(Y, BY, X, ::NotBWrapper)
    T = eltype(X)
    if Y isa SubArray
        gramYBY = @view Matrix{T}(I, size(Y.parent, 2), size(Y.parent, 2))[1:size(Y, 2), 1:size(Y, 2)]
        mul!(gramYBY, adjoint(Y), BY)
        gramYBV = @view zeros(T, size(Y.parent, 2), size(X, 2))[1:size(Y, 2), :]
    else
        gramYBY = adjoint(Y)*BY
        gramYBV = zeros(T, size(Y, 2), size(X, 2))
    end
    realdiag!(gramYBY)
    gramYBY_chol = cholesky!(Hermitian(gramYBY))
    tmp = deepcopy(gramYBV)

    return Constraint{eltype(Y), typeof(Y), typeof(gramYBV), typeof(gramYBY_chol)}(Y, BY, gramYBY_chol, gramYBV, tmp)
end

function update!(c::Constraint, X, BX)
    sizeY = size(c.Y, 2)
    sizeX = size(X, 2)
    c.Y.parent[:, sizeY+1:sizeY+sizeX] .= X
    if X !== BX
        c.BY.parent[:, sizeY+1:sizeY+sizeX] .= BX
    end
    sizeY += sizeX
    Y = @view c.Y.parent[:, 1:sizeY]
    BY = @view c.BY.parent[:, 1:sizeY]
    c.Y = Y
    c.BY = BY
    gram_chol = c.gram_chol
    new_factors = @view gram_chol.factors.parent[1:sizeY, 1:sizeY]
    c.gram_chol = typeof(gram_chol)(new_factors, gram_chol.uplo, convert(LinearAlgebra.BlasInt, 0))
    c.gramYBV = @view c.gramYBV.parent[1:sizeY, :]
    c.tmp = @view c.tmp.parent[1:sizeY, :]
    return c
end

function (constr!::Constraint{Nothing})(X, X_temp)
    nothing
end

function (constr!::Constraint)(X, X_temp)
    if size(constr!.Y, 2) > 0
        sizeX = size(X, 2)
        sizeY = size(constr!.Y, 2)
        gramYBV_view = view(constr!.gramYBV, 1:sizeY, 1:sizeX)
        mul!(gramYBV_view, adjoint(constr!.BY), X)
        tmp_view = view(constr!.tmp, 1:sizeY, 1:sizeX)
        ldiv!(tmp_view, constr!.gram_chol, gramYBV_view)
        mul!(X_temp, constr!.Y, tmp_view)
        @inbounds X .= X .- X_temp
    end
    nothing
end

struct RPreconditioner{TM, T, TA<:AbstractArray{T}}
    M::TM
    buffer::TA
    RPreconditioner{TM, T, TA}(M, X) where {TM, T, TA<:AbstractArray{T}} = new(M, similar(X))
end
RPreconditioner(M, X) = RPreconditioner{typeof(M), eltype(X), typeof(X)}(M, X)

function (precond!::RPreconditioner{Nothing})(X)
    nothing
end
function (precond!::RPreconditioner)(X)
    bs = size(X, 2)
    ldiv!(view(precond!.buffer, :, 1:bs), precond!.M, X)
    # Just returning buffer would be cheaper but struct at call site must be mutable
    @inbounds X .= @view precond!.buffer[:, 1:bs]
    nothing
end

struct BlockGram{Generalized, TA}
    XAX::TA
    XAR::TA
    XAP::TA
    RAR::TA
    RAP::TA
    PAP::TA
end
function BlockGram(XBlocks::Blocks{Generalized, T}) where {Generalized, T}
    sizeX = size(XBlocks.block, 2)
    XAX = zeros(T, sizeX, sizeX)
    XAP = zeros(T, sizeX, sizeX)
    XAR = zeros(T, sizeX, sizeX)
    RAR = zeros(T, sizeX, sizeX)
    RAP = zeros(T, sizeX, sizeX)
    PAP = zeros(T, sizeX, sizeX)
    return BlockGram{Generalized, Matrix{T}}(XAX, XAR, XAP, RAR, RAP, PAP)
end
XAX!(BlockGram, XBlocks) = mul!(BlockGram.XAX, adjoint(XBlocks.block), XBlocks.A_block)
XAP!(BlockGram, XBlocks, PBlocks, n) = mul!(view(BlockGram.XAP, :, 1:n), adjoint(XBlocks.block), view(PBlocks.A_block, :, 1:n))
XAR!(BlockGram, XBlocks, RBlocks, n) = mul!(view(BlockGram.XAR, :, 1:n), adjoint(XBlocks.block), view(RBlocks.A_block, :, 1:n))
RAR!(BlockGram, RBlocks, n) = mul!(view(BlockGram.RAR, 1:n, 1:n), adjoint(view(RBlocks.block, :, 1:n)), view(RBlocks.A_block, :, 1:n))
RAP!(BlockGram, RBlocks, PBlocks, n) = mul!(view(BlockGram.RAP, 1:n, 1:n), adjoint(view(RBlocks.A_block, :, 1:n)), view(PBlocks.block, :, 1:n))
PAP!(BlockGram, PBlocks, n) = mul!(view(BlockGram.PAP, 1:n, 1:n), adjoint(view(PBlocks.block, :, 1:n)), view(PBlocks.A_block, :, 1:n))
XBP!(BlockGram, XBlocks, PBlocks, n) = mul!(view(BlockGram.XAP, :, 1:n), adjoint(XBlocks.block), view(PBlocks.B_block, :, 1:n))
XBR!(BlockGram, XBlocks, RBlocks, n) = mul!(view(BlockGram.XAR, :, 1:n), adjoint(XBlocks.block), view(RBlocks.B_block, :, 1:n))
RBP!(BlockGram, RBlocks, PBlocks, n) = mul!(view(BlockGram.RAP, 1:n, 1:n), adjoint(view(RBlocks.B_block, :, 1:n)), view(PBlocks.block, :, 1:n))
#XBX!(BlockGram, XBlocks) = mul!(BlockGram.XAX, adjoint(XBlocks.block), XBlocks.B_block)
#RBR!(BlockGram, RBlocks, n) = mul!(view(BlockGram.RAR, 1:n, 1:n), adjoint(view(RBlocks.block, :, 1:n)), view(RBlocks.B_block, :, 1:n))
#PBP!(BlockGram, PBlocks, n) = mul!(view(BlockGram.PAP, 1:n, 1:n), adjoint(view(PBlocks.block, :, 1:n)), view(PBlocks.B_block, :, 1:n))

function I!(G, xr)
    @inbounds for j in xr, i in xr
        G[i, j] = ifelse(i==j, 1, 0)
    end
    return
end

function (g::BlockGram)(gram, lambda, n1::Int, n2::Int, n3::Int)
    xr = 1:n1
    rr = n1+1:n1+n2
    pr = n1+n2+1:n1+n2+n3
    @inbounds begin
        if n1 > 0
            #gram[xr, xr] .= view(g.XAX, 1:n1, 1:n1)
            gram[xr, xr] .= Diagonal(view(lambda, 1:n1))
        end
        if n2 > 0
            gram[rr, rr] .= view(g.RAR, 1:n2, 1:n2)
            gram[xr, rr] .= view(g.XAR, 1:n1, 1:n2)
            conj!(transpose!(view(gram, rr, xr), view(g.XAR, 1:n1, 1:n2)))
        end
        if n3 > 0
            gram[pr, pr] .= view(g.PAP, 1:n3, 1:n3)
            gram[rr, pr] .= view(g.RAP, 1:n2, 1:n3)
            gram[xr, pr] .= view(g.XAP, 1:n1, 1:n3)
            conj!(transpose!(view(gram, pr, rr), view(g.RAP, 1:n2, 1:n3)))
            conj!(transpose!(view(gram, pr, xr), view(g.XAP, 1:n1, 1:n3)))
        end
    end
    return
end
function (g::BlockGram)(gram, n1::Int, n2::Int, n3::Int, normalized::Bool=true)
    xr = 1:n1
    rr = n1+1:n1+n2
    pr = n1+n2+1:n1+n2+n3
    if n1 > 0
        if normalized
            I!(gram, xr)
        #else
        #    @inbounds gram[xr, xr] .= view(g.XAX, 1:n1, 1:n1)
        end
    end
    if n2 > 0
        if normalized
            I!(gram, rr)
        #else
        #    @inbounds gram[rr, rr] .= view(g.RAR, 1:n2, 1:n2)
        end
        @inbounds gram[xr, rr] .= view(g.XAR, 1:n1, 1:n2)
        @inbounds conj!(transpose!(view(gram, rr, xr), view(g.XAR, 1:n1, 1:n2)))
    end
    if n3 > 0
        if normalized
            I!(gram, pr)
        #else
        #    @inbounds gram[pr, pr] .= view(g.PAP, 1:n3, 1:n3)
        end
        @inbounds gram[rr, pr] .= view(g.RAP, 1:n2, 1:n3)
        @inbounds gram[xr, pr] .= view(g.XAP, 1:n1, 1:n3)
        @inbounds conj!(transpose!(view(gram, pr, rr), view(g.RAP, 1:n2, 1:n3)))
        @inbounds conj!(transpose!(view(gram, pr, xr), view(g.XAP, 1:n1, 1:n3)))
    end
    return
end

abstract type AbstractOrtho end
struct CholQR{TA} <: AbstractOrtho
    gramVBV::TA # to be used in view
end

function rdiv!(A, B::UpperTriangular)
    s = size(A, 2)
    @inbounds A[:,1] .= view(A, :, 1) ./ B[1,1]
    @inbounds for i in 2:s
        for j in 1:i-1
            A[:,i] .= view(A, :, i) .- view(A, :, j) .* B[j,i]
        end
        A[:,i] .= view(A, :, i) ./ B[i,i]
    end
    return A
end

realdiag!(M) = M
function realdiag!(M::AbstractMatrix{TC}) where TC <: Complex
    @inbounds for i in 1:minimum(size(M))
        M[i,i] = real(M[i,i])
    end
    return M
end

function (ortho!::CholQR)(XBlocks::Blocks{Generalized}, sizeX = -1; update_AX=false, update_BX=false) where Generalized
    useview = sizeX != -1
    if sizeX == -1
        sizeX = size(XBlocks.block, 2)
    end
    X = XBlocks.block
    BX = XBlocks.B_block # Assumes it is premultiplied
    AX = XBlocks.A_block
    gram_view = view(ortho!.gramVBV, 1:sizeX, 1:sizeX)
    if useview
        mul!(gram_view, adjoint(view(X, :, 1:sizeX)), view(BX, :, 1:sizeX))
    else
        mul!(gram_view, adjoint(X), BX)
    end
    realdiag!(gram_view)
    cholf = cholesky!(Hermitian(gram_view))
    R = cholf.factors
    if useview
        rdiv!(view(X, :, 1:sizeX), UpperTriangular(R))
        update_AX && rdiv!(view(AX, :, 1:sizeX), UpperTriangular(R))
        Generalized && update_BX && rdiv!(view(BX, :, 1:sizeX), UpperTriangular(R))
    else
        rdiv!(X, UpperTriangular(R))
        update_AX && rdiv!(AX, UpperTriangular(R))
        Generalized && update_BX && rdiv!(BX, UpperTriangular(R))
    end

    return
end

struct LOBPCGIterator{Generalized, T, TA, TB, TL<:AbstractVector{T}, TR<:AbstractVector, TPerm<:AbstractVector{Int}, TV<:AbstractArray{T}, TBlocks<:Blocks{Generalized, T}, TO<:AbstractOrtho, TP, TC, TG, TM, TTrace}
    A::TA
    B::TB
    ritz_values::TL
    λperm::TPerm
    λ::TL
    V::TV
    residuals::TR
    largest::Bool
    XBlocks::TBlocks
    tempXBlocks::TBlocks
    PBlocks::TBlocks
    activePBlocks::TBlocks # to be used in view
    RBlocks::TBlocks
    activeRBlocks::TBlocks # to be used in view
    iteration::Base.RefValue{Int}
    currentBlockSize::Base.RefValue{Int}
    ortho!::TO
    precond!::TP
    constr!::TC
    gramABlock::TG
    gramBBlock::TG
    gramA::TV
    gramB::TV
    activeMask::TM
    trace::TTrace
end

"""
    LOBPCGIterator(A, largest::Bool, X, P=nothing, C=nothing) -> iterator

# Arguments

- `A`: linear operator;
- `X`: initial guess of the Ritz vectors; to be overwritten by the eigenvectors;
- `largest`: `true` if largest eigenvalues are desired and false if smallest;
- `P`: preconditioner of residual vectors, must overload `ldiv!`;
- `C`: constraint to deflate the residual and solution vectors orthogonal
    to a subspace; must overload `mul!`.
"""
LOBPCGIterator(A, largest::Bool, X, P=nothing, C=nothing) = LOBPCGIterator(A, nothing, largest, X, P, C)

"""
    LOBPCGIterator(A, B, largest::Bool, X, P=nothing, C=nothing) -> iterator

# Arguments

- `A`: linear operator;
- `B`: linear operator;
- `X`: initial guess of the Ritz vectors; to be overwritten by the eigenvectors;
- `largest`: `true` if largest eigenvalues are desired and false if smallest;
- `P`: preconditioner of residual vectors, must overload `ldiv!`;
- `C`: constraint to deflate the residual and solution vectors orthogonal
    to a subspace; must overload `mul!`;
"""
function LOBPCGIterator(A, B, largest::Bool, X, P=nothing, C=nothing)
    precond! = RPreconditioner(P, X)
    constr! = Constraint(C, B, X)
    return LOBPCGIterator(A, B, largest, X, precond!, constr!)
end
function LOBPCGIterator(A, B, largest::Bool, X, precond!::RPreconditioner, constr!::Constraint)
    T = eltype(X)
    nev = size(X, 2)
    if B isa Nothing
        XBlocks = Blocks(X, similar(X))
        tempXBlocks = Blocks(copy(X), similar(X))
        RBlocks = Blocks(similar(X), similar(X))
        activeRBlocks = Blocks(similar(X), similar(X))
        PBlocks = Blocks(similar(X), similar(X))
        activePBlocks = Blocks(similar(X), similar(X))
    else
        XBlocks = Blocks(X, similar(X), similar(X))
        tempXBlocks = Blocks(copy(X), similar(X), similar(X))
        RBlocks = Blocks(similar(X), similar(X), similar(X))
        activeRBlocks = Blocks(similar(X), similar(X), similar(X))
        PBlocks = Blocks(similar(X), similar(X), similar(X))
        activePBlocks = Blocks(similar(X), similar(X), similar(X))
    end
    ritz_values = zeros(T, nev*3)
    λ = zeros(T, nev)
    λperm = zeros(Int, nev*3)
    V = zeros(T, nev*3, nev*3)
    residuals = fill(real(T)(NaN), nev)
    iteration = Ref(1)
    currentBlockSize = Ref(nev)
    generalized = !(B isa Nothing)
    ortho! = CholQR(zeros(T, nev, nev))

    gramABlock = BlockGram(XBlocks)
    gramBBlock = BlockGram(XBlocks)

    gramA = zeros(T, 3*nev, 3*nev)
    gramB = zeros(T, 3*nev, 3*nev)

    activeMask = ones(Bool, nev)
    trace = LOBPCGTrace{Vector{real(T)},Vector{T}}()

    return LOBPCGIterator{generalized, T, typeof(A), typeof(B), typeof(λ), typeof(residuals), typeof(λperm), typeof(V), typeof(XBlocks), typeof(ortho!), typeof(precond!), typeof(constr!), typeof(gramABlock), typeof(activeMask), typeof(trace)}(A, B, ritz_values, λperm, λ, V, residuals, largest, XBlocks, tempXBlocks, PBlocks, activePBlocks, RBlocks, activeRBlocks, iteration, currentBlockSize, ortho!, precond!, constr!, gramABlock, gramBBlock, gramA, gramB, activeMask, trace)
end
function LOBPCGIterator(A, largest::Bool, X, nev::Int, P=nothing, C=nothing)
    LOBPCGIterator(A, nothing, largest, X, nev, P, C)
end
function LOBPCGIterator(A, B, largest::Bool, X, nev::Int, P=nothing, C=nothing)
    T = eltype(X)
    n = size(X, 1)
    sizeX = size(X, 2)
    if C isa Nothing
        sizeC = 0
        new_C = typeof(X)(undef, n, (nev÷sizeX)*sizeX)
    else
        sizeC = size(C,2)
        new_C = typeof(C)(undef, n, sizeC+(nev÷sizeX)*sizeX)
        new_C[:,1:sizeC] .= C
    end
    if B isa Nothing
        new_BC = new_C
    else
        new_BC = similar(new_C)
    end
    Y = @view new_C[:, 1:sizeC]
    BY = @view new_BC[:, 1:sizeC]
    if !(B isa Nothing)
        mul!(BY, B, Y)
    end
    constr! = Constraint(Y, BY, X, NotBWrapper())
    precond! = RPreconditioner(P, X)
    return LOBPCGIterator(A, B, largest, X, precond!, constr!)
end

function ortho_AB_mul_X!(blocks::Blocks, ortho!, A, B, bs=-1)
    # Finds BX
    bs == -1 ? B_mul_X!(blocks, B) : B_mul_X!(blocks, B, bs)
    # Orthonormalizes X and updates BX
    bs == -1 ? ortho!(blocks, update_BX=true) : ortho!(blocks, bs, update_BX=true)
    # Updates AX
    bs == -1 ? A_mul_X!(blocks, A) : A_mul_X!(blocks, A, bs)
    return
end
function residuals!(iterator)
    sizeX = size(iterator.XBlocks.block, 2)
    mul!(iterator.RBlocks.block, iterator.XBlocks.B_block, Diagonal(view(iterator.ritz_values, 1:sizeX)))
    @inbounds iterator.RBlocks.block .= iterator.XBlocks.A_block .- iterator.RBlocks.block
    # Finds residual norms
    @inbounds for j in 1:size(iterator.RBlocks.block, 2)
        iterator.residuals[j] = 0
        for i in 1:size(iterator.RBlocks.block, 1)
            x = iterator.RBlocks.block[i,j]
            iterator.residuals[j] += real(x*conj(x))
        end
        iterator.residuals[j] = sqrt(iterator.residuals[j])
    end
    return
end

function update_mask!(iterator, residualTolerance)
    sizeX = size(iterator.XBlocks.block, 2)
    # Update active vectors mask
    @inbounds iterator.activeMask .= view(iterator.residuals, 1:sizeX) .> residualTolerance
    iterator.currentBlockSize[] = sum(iterator.activeMask)
    return
end

function update_active!(mask, bs::Int, blockPairs...)
    @inbounds for (activeblock, block) in blockPairs
        activeblock[:, 1:bs] .= view(block, :, mask)
    end
    return
end

function precond_constr!(block, temp_block, bs, precond!, constr!)
    precond!(view(block, :, 1:bs))
    # Constrain the active residual vectors to be B-orthogonal to Y
    constr!(view(block, :, 1:bs), view(temp_block, :, 1:bs))
    return
end
function block_grams_1x1!(iterator)
    # Finds gram matrix X'AX
    XAX!(iterator.gramABlock, iterator.XBlocks)
    return
end
function block_grams_2x2!(iterator, bs)
    sizeX = size(iterator.XBlocks.block, 2)
    #XAX!(iterator.gramABlock, iterator.XBlocks)
    XAR!(iterator.gramABlock, iterator.XBlocks, iterator.activeRBlocks, bs)
    RAR!(iterator.gramABlock, iterator.activeRBlocks, bs)
    XBR!(iterator.gramBBlock, iterator.XBlocks, iterator.activeRBlocks, bs)
    iterator.gramABlock(iterator.gramA, view(iterator.ritz_values, 1:sizeX), sizeX, bs, 0)
    iterator.gramBBlock(iterator.gramB, sizeX, bs, 0, true)

    return
end
function block_grams_3x3!(iterator, bs)
    # Find R'AR, P'AP, X'AR, X'AP and R'AP
    sizeX = size(iterator.XBlocks.block, 2)
    #XAX!(iterator.gramABlock, iterator.XBlocks)
    XAR!(iterator.gramABlock, iterator.XBlocks, iterator.activeRBlocks, bs)
    XAP!(iterator.gramABlock, iterator.XBlocks, iterator.activePBlocks, bs)
    RAR!(iterator.gramABlock, iterator.activeRBlocks, bs)
    RAP!(iterator.gramABlock, iterator.activeRBlocks, iterator.activePBlocks, bs)
    PAP!(iterator.gramABlock, iterator.activePBlocks, bs)
    # Find X'BR, X'BP and P'BR
    XBR!(iterator.gramBBlock, iterator.XBlocks, iterator.activeRBlocks, bs)
    XBP!(iterator.gramBBlock, iterator.XBlocks, iterator.activePBlocks, bs)
    RBP!(iterator.gramBBlock, iterator.activeRBlocks, iterator.activePBlocks, bs)
    # Update the gram matrix [X R P]' A [X R P]
    iterator.gramABlock(iterator.gramA, view(iterator.ritz_values, 1:sizeX), sizeX, bs, bs)
    # Update the gram matrix [X R P]' B [X R P]
    iterator.gramBBlock(iterator.gramB, sizeX, bs, bs, true)

    return
end

function sub_problem!(iterator, sizeX, bs1, bs2)
    subdim = sizeX+bs1+bs2
    if bs1 == 0
        gramAview = view(iterator.gramABlock.XAX, 1:subdim, 1:subdim)
        # Source of type instability
        realdiag!(gramAview)
        eigf = eigen!(Hermitian(gramAview))
    else
        gramAview = view(iterator.gramA, 1:subdim, 1:subdim)
        gramBview = view(iterator.gramB, 1:subdim, 1:subdim)
        # Source of type instability
        realdiag!(gramAview)
        realdiag!(gramBview)
        eigf = eigen!(Hermitian(gramAview), Hermitian(gramBview))
    end
    # Selects extremal eigenvalues and corresponding vectors
    partialsortperm!(view(iterator.λperm, 1:subdim), eigf.values, 1:subdim; rev=iterator.largest)
    @inbounds iterator.ritz_values[1:sizeX] .= view(eigf.values, view(iterator.λperm, 1:sizeX))
    @inbounds iterator.V[1:subdim, 1:sizeX] .= view(eigf.vectors, :, view(iterator.λperm, 1:sizeX))
    return
end

function update_X_P!(iterator::LOBPCGIterator{Generalized}, bs1, bs2) where Generalized
    sizeX = size(iterator.XBlocks.block, 2)
    x_eigview = view(iterator.V, 1:sizeX, 1:sizeX)
    r_eigview = view(iterator.V, sizeX+1:sizeX+bs1, 1:sizeX)
    p_eigview = view(iterator.V, sizeX+bs1+1:sizeX+bs1+bs2, 1:sizeX)
    r_blockview = view(iterator.activeRBlocks.block, :, 1:bs1)
    ra_blockview = view(iterator.activeRBlocks.A_block, :, 1:bs1)
    p_blockview = view(iterator.activePBlocks.block, :, 1:bs2)
    pa_blockview = view(iterator.activePBlocks.A_block, :, 1:bs2)
    if Generalized
        rb_blockview = view(iterator.activeRBlocks.B_block, :, 1:bs1)
        pb_blockview = view(iterator.activePBlocks.B_block, :, 1:bs2)
    end
    if bs1 > 0
        mul!(iterator.PBlocks.block, r_blockview, r_eigview)
        mul!(iterator.PBlocks.A_block, ra_blockview, r_eigview)
        if Generalized
            mul!(iterator.PBlocks.B_block, rb_blockview, r_eigview)
        end
    end
    if bs2 > 0
        mul!(iterator.tempXBlocks.block, p_blockview, p_eigview)
        mul!(iterator.tempXBlocks.A_block, pa_blockview, p_eigview)
        if Generalized
            mul!(iterator.tempXBlocks.B_block, pb_blockview, p_eigview)
        end
        @inbounds iterator.PBlocks.block .= iterator.PBlocks.block .+ iterator.tempXBlocks.block
        @inbounds iterator.PBlocks.A_block .= iterator.PBlocks.A_block .+ iterator.tempXBlocks.A_block
        if Generalized
            @inbounds iterator.PBlocks.B_block .= iterator.PBlocks.B_block .+ iterator.tempXBlocks.B_block
        end
    end
    block = iterator.XBlocks.block
    tempblock = iterator.tempXBlocks.block
    mul!(tempblock, block, x_eigview)
    block = iterator.XBlocks.A_block
    tempblock = iterator.tempXBlocks.A_block
    mul!(tempblock, block, x_eigview)
    if Generalized
        block = iterator.XBlocks.B_block
        tempblock = iterator.tempXBlocks.B_block
        mul!(tempblock, block, x_eigview)
    end
    @inbounds begin
        if bs1 > 0
            iterator.XBlocks.block .= iterator.tempXBlocks.block .+ iterator.PBlocks.block
            iterator.XBlocks.A_block .= iterator.tempXBlocks.A_block .+ iterator.PBlocks.A_block
            if Generalized
                iterator.XBlocks.B_block .= iterator.tempXBlocks.B_block .+ iterator.PBlocks.B_block
            end
        else
            iterator.XBlocks.block .= iterator.tempXBlocks.block
            iterator.XBlocks.A_block .= iterator.tempXBlocks.A_block
            if Generalized
                iterator.XBlocks.B_block .= iterator.tempXBlocks.B_block
            end
        end
    end
    return
end

function (iterator::LOBPCGIterator{Generalized})(residualTolerance, log) where {Generalized}
    sizeX = size(iterator.XBlocks.block, 2)
    iteration = iterator.iteration[]
    if iteration == 1
        ortho_AB_mul_X!(iterator.XBlocks, iterator.ortho!, iterator.A, iterator.B)
        # Finds gram matrix X'AX
        block_grams_1x1!(iterator)
        sub_problem!(iterator, sizeX, 0, 0)
        # Updates Ritz vectors X and updates AX and BX accordingly
        update_X_P!(iterator, 0, 0)
        residuals!(iterator)
        update_mask!(iterator, residualTolerance)
    elseif iteration == 2
        bs = iterator.currentBlockSize[]
        # Update active R blocks
        update_active!(iterator.activeMask, bs, (iterator.activeRBlocks.block, iterator.RBlocks.block))
        # Precondition and constrain the active residual vectors
        precond_constr!(iterator.activeRBlocks.block, iterator.tempXBlocks.block, bs, iterator.precond!, iterator.constr!)
        # Orthonormalizes R[:,1:bs] and finds AR[:,1:bs] and BR[:,1:bs]
        ortho_AB_mul_X!(iterator.activeRBlocks, iterator.ortho!, iterator.A, iterator.B, bs)
        # Find [X R] A [X R] and [X R]' B [X R]
        block_grams_2x2!(iterator, bs)
        # Solve the Rayleigh-Ritz sub-problem
        sub_problem!(iterator, sizeX, bs, 0)
        update_X_P!(iterator, bs, 0)
        residuals!(iterator)
        update_mask!(iterator, residualTolerance)
    else
        # Update active blocks
        bs = iterator.currentBlockSize[]
        # Update active R and P blocks
        update_active!(iterator.activeMask, bs,
            (iterator.activeRBlocks.block, iterator.RBlocks.block),
            (iterator.activePBlocks.block, iterator.PBlocks.block),
            (iterator.activePBlocks.A_block, iterator.PBlocks.A_block),
            (iterator.activePBlocks.B_block, iterator.PBlocks.B_block))
        # Precondition and constrain the active residual vectors
        precond_constr!(iterator.activeRBlocks.block, iterator.tempXBlocks.block, bs, iterator.precond!, iterator.constr!)
        # Orthonormalizes R[:,1:bs] and finds AR[:,1:bs] and BR[:,1:bs]
        ortho_AB_mul_X!(iterator.activeRBlocks, iterator.ortho!, iterator.A, iterator.B, bs)
        # Orthonormalizes P and updates AP
        iterator.ortho!(iterator.activePBlocks, bs, update_AX=true, update_BX=true)
        # Update the gram matrix [X R P]' A [X R P] and [X R P]' B [X R P]
        block_grams_3x3!(iterator, bs)
        # Solve the Rayleigh-Ritz sub-problem
        sub_problem!(iterator, sizeX, bs, bs)
        # Updates Ritz vectors X and updates AX and BX accordingly
        # And updates P, AP and BP
        update_X_P!(iterator, bs, bs)
        residuals!(iterator)
        update_mask!(iterator, residualTolerance)
    end
    if log
        return LOBPCGState(iteration, iterator.residuals[1:sizeX], iterator.ritz_values[1:sizeX])
    else
        return LOBPCGState(iteration, nothing, nothing)
    end
end

default_tolerance(::Type{T}) where {T<:Number} = eps(real(T))^(real(T)(3)/10)

"""
The Locally Optimal Block Preconditioned Conjugate Gradient Method (LOBPCG)

Finds the `nev` extremal eigenvalues and their corresponding eigenvectors satisfying `AX = λBX`.

`A` and `B` may be generic types but `Base.mul!(C, AorB, X)` must be defined for vectors and strided matrices `X` and `C`. `size(A, i::Int)` and `eltype(A)` must also be defined for `A`.

    lobpcg(A, [B,] largest, nev; kwargs...) -> results

# Arguments

- `A`: linear operator;
- `B`: linear operator;
- `largest`: `true` if largest eigenvalues are desired and false if smallest;
- `nev`: number of eigenvalues desired.

## Keywords

- `log::Bool`: default is `false`; if `true`, `results.trace` will store iterations
    states; if `false` only `results.trace` will be empty;

- `P`: preconditioner of residual vectors, must overload `ldiv!`;

- `C`: constraint to deflate the residual and solution vectors orthogonal
    to a subspace; must overload `mul!`;

- `maxiter`: maximum number of iterations; default is 200;

- `tol::Real`: tolerance to which residual vector norms must be under.

# Output

- `results`: a `LOBPCGResults` struct. `r.λ` and `r.X` store the eigenvalues and eigenvectors.
"""
function lobpcg(A, largest::Bool, nev::Int; kwargs...)
    lobpcg(A, nothing, largest, nev; kwargs...)
end
function lobpcg(A, B, largest::Bool, nev::Int; kwargs...)
    lobpcg(A, B, largest, rand(eltype(A), size(A, 1), nev); not_zeros=true, kwargs...)
end

"""
    lobpcg(A, [B,] largest, X0; kwargs...) -> results

# Arguments

- `A`: linear operator;
- `B`: linear operator;
- `largest`: `true` if largest eigenvalues are desired and false if smallest;
- `X0`: Initial guess, will not be modified. The number of columns is the number of eigenvectors desired.

## Keywords

- `not_zeros`: default is `false`. If `true`, `X0` will be assumed to not have any all-zeros column.

- `log::Bool`: default is `false`; if `true`, `results.trace` will store iterations
    states; if `false` only `results.trace` will be empty;

- `P`: preconditioner of residual vectors, must overload `ldiv!`;

- `C`: constraint to deflate the residual and solution vectors orthogonal
    to a subspace; must overload `mul!`;

- `maxiter`: maximum number of iterations; default is 200;

- `tol::Real`: tolerance to which residual vector norms must be under.

# Output

- `results`: a `LOBPCGResults` struct. `r.λ` and `r.X` store the eigenvalues and eigenvectors.
"""
function lobpcg(A, largest::Bool, X0; kwargs...)
    lobpcg(A, nothing, largest, X0; kwargs...)
end
function lobpcg(A, B, largest, X0;
                not_zeros=false, log=false, P=nothing, maxiter=200,
                C=nothing, tol::Real=default_tolerance(eltype(X0)))
    X = copy(X0)
    n = size(X, 1)
    sizeX = size(X, 2)
    sizeX > n && throw("X column dimension exceeds the row dimension")

    iterator = LOBPCGIterator(A, B, largest, X, P, C)

    return lobpcg!(iterator, log=log, tol=tol, maxiter=maxiter, not_zeros=not_zeros)
end

"""
    lobpcg!(iterator::LOBPCGIterator; kwargs...) -> results

# Arguments

- `iterator::LOBPCGIterator`: a struct having all the variables required
    for the LOBPCG algorithm.

## Keywords

- `not_zeros`: default is `false`. If `true`, the initial Ritz vectors will be assumed to not have any all-zeros column.

- `log::Bool`: default is `false`; if `true`, `results.trace` will store iterations
    states; if `false` only `results.trace` will be empty;

- `maxiter`: maximum number of iterations; default is 200;

- `tol::Real`: tolerance to which residual vector norms must be under.

# Output

- `results`: a `LOBPCGResults` struct. `r.λ` and `r.X` store the eigenvalues and eigenvectors.

"""
function lobpcg!(iterator::LOBPCGIterator; log=false, maxiter=200, not_zeros=false,
                 tol::Real=default_tolerance(eltype(iterator.XBlocks.block)))
    X = iterator.XBlocks.block
    iterator.constr!(iterator.XBlocks.block, iterator.tempXBlocks.block)
    if !not_zeros
        for j in 1:size(X,2)
            if all(x -> x==0, view(X, :, j))
                @inbounds X[:,j] .= rand.()
            end
        end
        iterator.constr!(iterator.XBlocks.block, iterator.tempXBlocks.block)
    end
    n = size(X, 1)
    sizeX = size(X, 2)
    iterator.iteration[] = 1
    while iterator.iteration[] <= maxiter
        state = iterator(tol, log)
        if log
            push!(iterator.trace, state)
        end
        iterator.currentBlockSize[] == 0 && break
        iterator.iteration[] += 1
    end
    @inbounds iterator.λ .= view(iterator.ritz_values, 1:sizeX)

    results = LOBPCGResults(iterator.λ, X, tol, iterator.residuals, iterator.iteration[], maxiter, all((x)->(norm(x)<=tol), view(iterator.residuals, 1:sizeX)), iterator.trace)

    return results
end

"""
lobpcg(A, [B,] largest, X0, nev; kwargs...) -> results

# Arguments

- `A`: linear operator;
- `B`: linear operator;
- `largest`: `true` if largest eigenvalues are desired and false if smallest;
- `X0`: block vectors such that the eigenvalues will be found size(X0, 2) at a time;
    the columns are also used to initialize the first batch of Ritz vectors;
- `nev`: number of eigenvalues desired.

## Keywords

- `log::Bool`: default is `false`; if `true`, `results.trace` will store iterations
    states; if `false` only `results.trace` will be empty;

- `P`: preconditioner of residual vectors, must overload `ldiv!`;

- `C`: constraint to deflate the residual and solution vectors orthogonal
    to a subspace; must overload `mul!`;

- `maxiter`: maximum number of iterations; default is 200;

- `tol::Real`: tolerance to which residual vector norms must be under.

# Output

- `results`: a `LOBPCGResults` struct. `r.λ` and `r.X` store the eigenvalues and eigenvectors.
"""
function lobpcg(A, largest::Bool, X0, nev::Int; kwargs...)
    lobpcg(A, nothing, largest, X0, nev; kwargs...)
end
function lobpcg(A, B, largest::Bool, X0, nev::Int;
                not_zeros=false, log=false, P=nothing, maxiter=200,
                C=nothing, tol::Real=default_tolerance(eltype(X0)))
    n = size(X0, 1)
    sizeX = size(X0, 2)
    nev > n && throw("Number of eigenvectors desired exceeds the row dimension.")

    sizeX = min(nev, sizeX)
    X = X0[:, 1:sizeX]
    iterator = LOBPCGIterator(A, B, largest, X, nev, P, C)

    r = EmptyLOBPCGResults(X, nev, tol, maxiter)
    rnext = lobpcg!(iterator, log=log, tol=tol, maxiter=maxiter, not_zeros=not_zeros)
    append!(r, rnext, 0)
    converged_x = sizeX
    while converged_x < nev
        if nev-converged_x < sizeX
            cutoff = sizeX-(nev-converged_x)
            update!(iterator.constr!, view(iterator.XBlocks.block, :, 1:cutoff), view(iterator.XBlocks.B_block, :, 1:cutoff))
            X[:, 1:sizeX-cutoff] .= @view X[:, cutoff+1:sizeX]
            rand!(view(X, :, cutoff+1:sizeX))
            rnext = lobpcg!(iterator, log=log, tol=tol, maxiter=maxiter, not_zeros=true)
            append!(r, rnext, converged_x, sizeX-cutoff)
            converged_x += sizeX-cutoff
        else
            update!(iterator.constr!, iterator.XBlocks.block, iterator.XBlocks.B_block)
            rand!(X)
            rnext = lobpcg!(iterator, log=log, tol=tol, maxiter=maxiter, not_zeros=true)
            append!(r, rnext, converged_x)
            converged_x += sizeX
        end
    end
    return r
end


text:
module Scoring
    include("BaseScore.jl")
end  # module Scoring

ext: jl


text:
module S2Cell
using Markdown
using Printf

md"""Exception type for invalid cell IDs."""
struct InvalidCellID <: Exception
    msg:: String
end

md"""Exception type for invalid tokens."""
struct InvalidToken <: Exception
    msg:: String
end

md"""An S2 cell ID is an integer we package for type-checing"""
struct Cell
    id::UInt64
end

#
# S2 base constants needed for cell mapping
#

# The maximum level supported within an S2 cell ID. Each level is represented by two bits in the
# final cell ID
_S2_MAX_LEVEL = 30

# The maximum value within the I and J bits of an S2 cell ID
_S2_MAX_SIZE = 1 << _S2_MAX_LEVEL

# The number of bits in a S2 cell ID used for specifying the base face
_S2_FACE_BITS = 3

# The number of bits in a S2 cell ID used for specifying the position along the Hilbert curve
_S2_POS_BITS = 2 * _S2_MAX_LEVEL + 1

# The maximum value of the Si/Ti integers used when mapping from IJ to ST. This is twice the max
# value of I and J, since Si/Ti allow referencing both the center and edge of a leaf cell
_S2_MAX_SI_TI = 1 << (_S2_MAX_LEVEL + 1)

# Mask that specifies the swap orientation bit for the Hilbert curve
_S2_SWAP_MASK = 1

# Mask that specifies the invert orientation bit for the Hilbert curve
_S2_INVERT_MASK = 2

# The number of bits per I and J in the lookup tables
_S2_LOOKUP_BITS = 4

# Lookup table for mapping 10 bits of IJ + orientation to 10 bits of Hilbert curve position +
# orientation.
_S2_LOOKUP_POS = nothing

# Lookup table for mapping 10 bits of Hilbert curve position + orientation to 10 bits of IJ +
# orientation.
_S2_LOOKUP_IJ = nothing

# Lookup table of two bits of IJ from two bits of curve position, based also on the current curve
# orientation from the swap and invert bits
_S2_POS_TO_IJ = [
    [0, 1, 3, 2],  # 0: Normal order, no swap or invert
    [0, 2, 3, 1],  # 1: Swap bit set, swap I and J bits
    [3, 2, 0, 1],  # 2: Invert bit set, invert bits
    [3, 1, 0, 2],  # 3: Swap and invert bits set
]

# Lookup for the orientation update mask of one of the four sub-cells within a higher level cell.
# This mask is XOR'ed with the current orientation to get the sub-cell orientation.
_S2_POS_TO_ORIENTATION_MASK = [_S2_SWAP_MASK, 0, 0, _S2_SWAP_MASK | _S2_INVERT_MASK]

#
# S2 helper functions
#

md"""
    Convert S2 UV to ST.

    This is done using the quadratic projection that is used by default for S2. The C++ and Java S2
    libraries use a different definition of the ST cell-space, but the end result in IJ is the same.
    The below uses the C++ ST definition.

    See s2geometry/blob/c59d0ca01ae3976db7f8abdc83fcc871a3a95186/src/s2/s2coords.h#L317-L320

"""
function _s2_uv_to_st(component:: Real) :: Real
    if component >= 0.0
        return 0.5 * sqrt(1.0 + 3.0 * component)
    else
        return 1.0 - 0.5 * sqrt(1.0 - 3.0 * component)
    end
end

md"""
Convert S2 ST to UV.

This is done using the quadratic projection that is used by default for S2. The C++ and Java S2
libraries use a different definition of the ST cell-space, but the end result in IJ is the same.
The below uses the C++ ST definition.

See s2geometry/blob/c59d0ca01ae3976db7f8abdc83fcc871a3a95186/src/s2/s2coords.h#L312-L315

"""

function _s2_st_to_uv(component:: Real) :: Real
    if component >= 0.5
        return (1.0 / 3.0) * (4.0 * component^2 - 1.0)
    else
        return (1.0 / 3.0) * (1.0 - 4.0 * (1.0 - component)^2)
    end
end

md"""
Convert S2 ST to IJ.

The mapping here differs between C++ and Java versions, but the combination of
_st_to_ij(_uv_to_st(val)) is the same for both. The below uses the C++ ST definition.

See s2geometry/blob/2c02e21040e0b82aa5719e96033d02b8ce7c0eff/src/s2/s2coords.h#L333-L336

"""
function _s2_st_to_ij(component:: Real) :: Int
    # The reference implementation does round(_S2_MAX_SIZE * component - 0.5), which is equivalent
    # to math.floor(_S2_MAX_SIZE * component)
    return Int(max(0, min(_S2_MAX_SIZE - 1, floor(_S2_MAX_SIZE * component))))
end

md"""
Convert S2 Si/Ti to ST.

This converts an integer in range 0 to _S2_MAX_SI_TI into a float in range 0.0 to 1.0.

See s2geometry/blob/c59d0ca01ae3976db7f8abdc83fcc871a3a95186/src/s2/s2coords.h#L338-L341

"""
function _s2_si_ti_to_st(component:: Int) :: Real
    return (1.0 / _S2_MAX_SI_TI) * component
end

md"""
Convert face + UV to S2Point XYZ.

See s2geometry/blob/c59d0ca01ae3976db7f8abdc83fcc871a3a95186/src/s2/s2coords.h#L348-L357

Args:
    face: The S2 face for the input point.
    uv: The S2 face UV coordinates.

Returns:
    The unnormalised S2Point XYZ.

Raises:
    ValueError: If the face is not valid in range 0-5.

"""
function _s2_face_uv_to_xyz(face:: Integer, uv:: Tuple{<:Real,<:Real}):: Tuple{<:Real,<:Real,<:Real}
    # Face -> XYZ components -> indices with negation:
    # 0    -> ( 1,  u,  v)   -> ( /,  0,  1)
    # 1    -> (-u,  1,  v)   -> (-0,  /,  1)
    # 2    -> (-u, -v,  1)   -> (-0, -1,  /)
    # 3    -> (-1, -v, -u)   -> (-/, -1, -0) <- -1 here means -1 times the value in index 1,
    # 4    -> ( v, -1, -u)   -> ( 1, -/, -0)    not index -1
    # 5    -> ( v,  u, -1)   -> ( 1,  0, -/)
    if face == 0
        s2_point = (1, uv[1], uv[2])
    elseif face == 1
        s2_point = (-uv[1], 1, uv[2])
    elseif face == 2
        s2_point = (-uv[1], -uv[2], 1)
    elseif face == 3
        s2_point = (-1, -uv[2], -uv[1])
    elseif face == 4
        s2_point = (uv[2], -1, -uv[1])
    elseif face == 5
        s2_point = (uv[2], uv[1], -1)
    else
        throw(ArgumentError("Cannot convert UV to XYZ with invalid face: $face"))
    end

    return s2_point
end

#
# Cell ID <-> Token translation functions
#

md"""
Convert S2 cell ID to a S2 token.

Converts the S2 cell ID to hex and strips any trailing zeros. The 0 cell ID token is represented
as 'X' to prevent it being an empty string.

See s2geometry/blob/c59d0ca01ae3976db7f8abdc83fcc871a3a95186/src/s2/s2cell_id.cc#L204-L220

Args:
    cell_id: The S2 cell ID integer.

Returns:
    The S2 token string for the S2 cell ID.

Raises:
    TypeError: If the cell_id is not int.

"""
cell_id_to_token(cell:: Cell)::String = cell_id_to_token(cell.id)

function cell_id_to_token(cell_id:: UInt64)::String
    # The zero token is encoded as 'X' rather than as a zero-length string
    if cell_id == 0
        return "X"
    end

    # Convert cell ID to 16 character hex string and strip any implicit trailing zeros
    return rstrip((@sprintf "%016x" cell_id), '0')
end

md"""
Convert S2 token to S2 cell ID.

Restores the stripped 0 characters from the token and converts the hex string to integer.

See s2geometry/blob/c59d0ca01ae3976db7f8abdc83fcc871a3a95186/src/s2/s2cell_id.cc#L222-L239

Args:
    token: The S2 token string. Can be upper or lower case hex string.

Returns:
   The S2 cell ID for the S2 token.

Raises:
    TypeError: If the token is not str.
    InvalidToken: If the token length is over 16.

"""
function token_to_cell_id(token:: String)::Cell
    if length(token) == 0 || length(token) > 16
        throw(InvalidToken("Cannot convert S2 token with length > 16 characters"))
    end

    # Check for the zero cell ID represented by the character 'x' or 'X' rather than as the empty
    # string
    if lowercase(token) == "x"
        return Cell(0)
    end

    # Add stripped implicit zeros to create the full 16 character hex string
    token = token * ('0' ^ (16 - length(token)))

    # Convert to cell ID by converting hex to int
    return Cell(parse(UInt64, "0x" * token))
end

#
# Encode functions
#

md"""
Convert lat/lon to a S2 cell ID.

It is expected that the lat/lon provided are normalised, with latitude in the range -90 to 90.

Args:
    lat: The latitude to convert, in degrees.
    lon: The longitude to convert, in degrees.
    level: The level of the cell ID to generate, from 0 up to 30.

Returns:
    The S2 cell ID for the lat/lon location.

Raises:
    ValueError: When level is not an integer, is < 0 or is > 30.

"""
function lat_lon_to_cell_id(lat :: Real, lon :: Real, level :: Integer = 30)::Cell
    if level < 0 || level > _S2_MAX_LEVEL
        throw(ArgumentError("S2 level must be integer >= 0 and <= 30"))
    end

    # Populate _S2_LOOKUP_POS on first run.
    # See s2geometry/blob/c59d0ca01ae3976db7f8abdc83fcc871a3a95186/src/s2/s2cell_id.cc#L75-L109
    #
    # This table takes 10 bits of I and J and orientation and returns 10 bits of curve position and
    # new orientation

    # Reuse constant expressions
    lat_rad = deg2rad(lat)
    lon_rad = deg2rad(lon)
    sin_lat_rad = sin(lat_rad)
    cos_lat_rad = cos(lat_rad)
    sin_lon_rad = sin(lon_rad)
    cos_lon_rad = cos(lon_rad)

    # Convert to S2Point
    # This is effectively the unit non-geodetic ECEF vector
    s2_point = (
        cos_lat_rad * cos_lon_rad,  # X
        cos_lat_rad * sin_lon_rad,  # Y
        sin_lat_rad                 # Z
    )

    # Get cube face
    # See s2geometry/blob/2c02e21040e0b82aa5719e96033d02b8ce7c0eff/src/s2/s2coords.h#L380-L384
    #
    # The face is determined by the largest XYZ component of the S2Point vector. When the component
    # is negative, the second set of three faces is used.
    # Largest component -> face:
    # +x -> 0
    # +y -> 1
    # +z -> 2
    # -x -> 3
    # -y -> 4
    # -z -> 5
    face = findmax(abs.(s2_point))[2] - 1
    if s2_point[face + 1] < 0.0
        face += 3
    end

    # Convert face + XYZ to cube-space face + UV
    # See s2geometry/blob/2c02e21040e0b82aa5719e96033d02b8ce7c0eff/src/s2/s2coords.h#L366-L372
    #
    # The faces are oriented to ensure continuity of curve.
    # Face -> UV components -> indices with negation (without divisor, which is always the remaining
    # component (index: face % 3)):
    # 0 -> ( y,  z) -> ( 1,  2)
    # 1 -> (-x,  z) -> (-0,  2)
    # 2 -> (-x, -y) -> (-0, -1) <- -1 here means -1 times the value in index 1, not index -1
    # 3 -> ( z,  y) -> ( 2,  1)
    # 4 -> ( z, -x) -> ( 2, -0)
    # 5 -> (-y, -x) -> (-1, -0)
    #
    if face == 0
        uv =  (s2_point[2] / s2_point[1], s2_point[3] / s2_point[1])
    elseif face == 1
        uv = (-s2_point[1] / s2_point[2], s2_point[3] / s2_point[2])
    elseif face == 2
        uv = (-s2_point[1] / s2_point[3], -s2_point[2] / s2_point[3])
    elseif face == 3
        uv = (s2_point[3] / s2_point[1], s2_point[2] / s2_point[1])
    elseif face == 4
        uv = (s2_point[3] / s2_point[2], -s2_point[1] / s2_point[2])
    else
        uv = ( -s2_point[2] / s2_point[3], -s2_point[1] / s2_point[3])
    end

    # Project cube-space UV to cell-space ST
    # See s2geometry/blob/2c02e21040e0b82aa5719e96033d02b8ce7c0eff/src/s2/s2coords.h#L317-L320
    st = (_s2_uv_to_st(uv[1]), _s2_uv_to_st(uv[2]))

    # Convert ST to IJ integers
    # See s2geometry/blob/2c02e21040e0b82aa5719e96033d02b8ce7c0eff/src/s2/s2coords.h#L333-L336
    ij = (_s2_st_to_ij(st[1]), _s2_st_to_ij(st[2]))  # pylint: disable=invalid-name

    # Convert face + IJ to cell ID
    # See s2geometry/blob/c59d0ca01ae3976db7f8abdc83fcc871a3a95186/src/s2/s2cell_id.cc#L256-L298
    #
    # This is done by looking up 8 bits of I and J (4 each) at a time in the lookup table, along
    # with two bits of orientation (swap (1) and invert (2)). This gives back 8 bits of position
    # along the curve and two new orientation bits for the curve within the sub-cells in the next
    # step.
    #
    # The swap bit swaps I and J with each other
    # The invert bit inverts the bits of I and J, which means axes are negated
    #
    # Compared to the standard versions, we check the required number of steps we need to do for the
    # requested level and don't perform steps that will be completely overwritten in the truncation
    # below, rather than always doing every step. Each step does 4 bits each of I and J, which is 4
    # levels, so the required number of steps is ceil((level + 2) / 4), when level is > 0. The
    # additional 2 levels added are required to account for the top 3 bits (4 before right shift)
    # that are occupied by the face bits.
    bits = face & _S2_SWAP_MASK  # iiiijjjjoo. Initially set by by face
    cell_id = UInt64(face << (_S2_POS_BITS - 1))  # Insert face at second most signficant bits
    lookup_mask = (1 << _S2_LOOKUP_BITS) - 1  # Mask of 4 one bits: 0b1111
    if level > 0
        required_steps = Int(ceil((level + 2) / 4))
    else
        required_steps = 0
    end
    for k in 7:-1:(7 - required_steps)
        # Grab 4 bits of each of I and J
        offset = k * _S2_LOOKUP_BITS
        bits += ((ij[1] >> offset) & lookup_mask) << (_S2_LOOKUP_BITS + 2)
        bits += ((ij[2] >> offset) & lookup_mask) << 2

        # Map bits from iiiijjjjoo to ppppppppoo using lookup table
        bits = _S2_LOOKUP_POS[bits + 1]

        # Insert position bits into cell ID
        cell_id |= (bits >> 2) << (k * 2 * _S2_LOOKUP_BITS)

        # Remove position bits, leaving just new swap and invert bits for the next round
        bits &= _S2_SWAP_MASK | _S2_INVERT_MASK  # Mask: 0b11
    end

    # Left shift and add trailing bit
    # The trailing bit addition is disabled, as we are overwriting this below in the truncation
    # anyway. This line is kept as an example of the full method for S2 cell ID creation as is done
    # in the standard library versions.
    cell_id = cell_id << 1

    # Truncate to desired level
    # This is done by finding the mask of the trailing 1 bit for the specified level, then zeroing
    # out all bits less significant than this, then finally setting the trailing 1 bit. This is
    # still necessary to do even after a reduced number of steps `required_steps` above, since each
    # step contains multiple levels that may need partial overwrite. Additionally, we need to add
    # the trailing 1 bit, which is not yet set above.
    least_significant_bit_mask = 1 << (2 * (_S2_MAX_LEVEL - level))
    cell_id = (cell_id & -least_significant_bit_mask) | least_significant_bit_mask

    return Cell(cell_id)
end

md"""
Convert lat/lon to a S2 token.

Converts the S2 cell ID to hex and strips any trailing zeros. The 0 cell ID token is represented
as 'X' to prevent it being an empty string.

It is expected that the lat/lon provided are normalised, with latitude in the range -90 to 90.

See s2geometry/blob/c59d0ca01ae3976db7f8abdc83fcc871a3a95186/src/s2/s2cell_id.cc#L204-L220

Args:
    lat: The latitude to convert, in degrees.
    lon: The longitude to convert, in degrees.
    level: The level of the cell ID to generate, from 0 up to 30.

Returns:
    The S2 token string for the lat/lon location.

Raises:
    ValueError: When level is not an integer, is < 0 or is > 30.

"""
function lat_lon_to_token(lat:: Real, lon:: Real, level:: Integer = 30)::String
    # Generate cell ID and convert to token
    return cell_id_to_token(lat_lon_to_cell_id(lat, lon, level))
end

#
# Decode functions
#

md"""
Convert S2 cell ID to lat/lon.

Args:
    cell_id: The S2 cell ID integer.

Returns:
    The lat/lon (in degrees) tuple generated from the S2 cell ID.

Raises:
    TypeError: If the cell_id is not int.
    InvalidCellID: If the cell_id is invalid.

"""
cell_id_to_lat_lon(cell:: Cell):: Tuple{<:Real, <:Real} = cell_id_to_lat_lon(cell.id)

function cell_id_to_lat_lon(cell_id:: UInt64):: Tuple{<:Real, <:Real}
    # Check input
    if ! cell_id_is_valid(cell_id)
        throw(InvalidCellID("Cannot decode invalid S2 cell ID: $cell_id"))
    end

    # Populate _S2_LOOKUP_IJ on first run.
    # See s2geometry/blob/c59d0ca01ae3976db7f8abdc83fcc871a3a95186/src/s2/s2cell_id.cc#L75-L109
    # This table takes 10 bits of curve position and orientation and returns 10 bits of I and J and
    # new orientation

    # Extract face + IJ from cell ID
    # See s2geometry/blob/c59d0ca01ae3976db7f8abdc83fcc871a3a95186/src/s2/s2cell_id.cc#L312-L367
    #
    # This is done by looking up 8 bits of curve position at a time in the lookup table, along with
    # two bits of orientation (swap (1) and invert (2)). This gives back 8 bits of I and J (4 each)
    # and two new orientation bits for the curve within the sub-cells in the next step.
    #
    # The swap bit swaps I and J with each other
    # The invert bit inverts the bits of I and J, which means axes are negated
    #
    # In the first loop (most significant bits), the 3 bits occupied by the face need to be masked
    # out, since these are not set in the IJ to cell ID during encoding.
    #
    # The I and J returned here are of one of the two leaf (level 30) cells that are located
    # diagonally closest to the cell center. This happens because repeated ..00.. will select the
    # 'lower left' (for nominally oriented Hilbert curve segments) of the sub-cells. The ..10..
    # arising from the trailing bit, prior to the repeated ..00.. bits, ensures we first pick the
    # 'upper right' of the cell, then iterate in to lower left until we hit the leaf cell. This
    # means we pick the leaf cell to the north east of the parent cell center (again for nominal
    # orientation).
    # However, in the case of the swapped and inverted curve segment (4th sub-curve segment), the
    # ..10.. will select the 'lower left' and then iterate to the 'upper right' with each ..00..
    # following. In that case, we will be offset left and down by one leaf cell in each of I and J,
    # which needs to be fixed to have a consistent mapping. This is detectable by seeing that the
    # final bit of I or J is 1 (i.e we have picked an odd row/column, which will happen concurrently
    # in both I and J, so we only need to check one), except in case of level 29 where the logic is
    # inverted and the correction needs to be applied when we pick an even row/column (i.e I/J ends
    # in 0), since there are no trailing ..00..  available after the ``..10..`` when we are at level
    # 29+.
    #
    # This behaviour can be captured in the expression:
    # apply_correction = not leaf and (i ^ (is level 29)) & 1
    # apply_correction = not leaf and (i ^ (cell_id >> 2)) & 1
    #
    # We check for level 29 by looking for the trailing 1 in the third LSB, when we already know
    # that we are not a leaf cell (which could give false positive) by the initial check in the
    # expression.
    # See s2geometry/blob/c59d0ca01ae3976db7f8abdc83fcc871a3a95186/src/s2/s2cell_id.h#L503-L529
    #
    face = cell_id >> _S2_POS_BITS
    bits = face & _S2_SWAP_MASK  # ppppppppoo. Initially set by by face
    lookup_mask = (1 << _S2_LOOKUP_BITS) - 1  # Mask of 4 one bits: 0b1111
    i = 0
    j = 0
    for k in 7:-1:-1
        # Pull out 8 bits of cell ID, except in first loop where we pull out only 4
        if k == 7
            n_bits = (_S2_MAX_LEVEL - 7 * _S2_LOOKUP_BITS)
        else
            n_bits = _S2_LOOKUP_BITS
        end
        extract_mask = (1 << (2 * n_bits)) - 1  # 8 (or 4) one bits
        bits += ((cell_id >> (k * 2 * _S2_LOOKUP_BITS + 1)) & extract_mask) << 2

        # Map bits from ppppppppoo to iiiijjjjoo using lookup table
        bits = _S2_LOOKUP_IJ[bits + 1]

        # Extract I and J bits
        offset = k * _S2_LOOKUP_BITS
        i += (bits >> (_S2_LOOKUP_BITS + 2)) << offset  # Don't need lookup mask here
        j += ((bits >> 2) & lookup_mask) << offset

        # Remove I and J bits, leaving just new swap and invert bits for the next round
        bits &= _S2_SWAP_MASK | _S2_INVERT_MASK  # Mask: 0b11
    end

    # Resolve the center of the cell. For leaf cells, we add half the leaf cell size. For non-leaf
    # cells, we currently have one of either two cells diagonally around the cell center and want
    # to pick the leaf-cell edges that represent the parent cell center, as described above. The
    # center_correction_delta is 2x the offset, as we left shift I and J first.
    # This gives us the values Si and Ti, which are discrete representation of S and T in range 0 to
    # _S2_MAX_SI_TI. The extra power of 2 over IJ allows for identifying both the center and edge of
    # cells, whilst IJ is just the leaf cells.
    # See s2geometry/blob/c59d0ca01ae3976db7f8abdc83fcc871a3a95186/src/s2/s2coords.h#L57-L65
    is_leaf = Bool(cell_id & 1)  # Cell is leaf cell when trailing one bit is in LSB
    apply_correction = !is_leaf && Bool((i ⊻ (cell_id >> 2)) & 1)
    if is_leaf
        correction_delta = 1
    else
        correction_delta = apply_correction ? 2 : 0
    end
    si = (i << 1) + correction_delta
    ti = (j << 1) + correction_delta

    # Convert integer si/ti to double ST
    # See s2geometry/blob/c59d0ca01ae3976db7f8abdc83fcc871a3a95186/src/s2/s2coords.h#L338-L341
    st = (_s2_si_ti_to_st(si), _s2_si_ti_to_st(ti))

    # Project cell-space ST to cube-space UV
    # See s2geometry/blob/c59d0ca01ae3976db7f8abdc83fcc871a3a95186/src/s2/s2coords.h#L312-L315
    uv = (_s2_st_to_uv(st[1]), _s2_st_to_uv(st[2]))

    # Convert face + UV to S2Point XYZ
    # See s2geometry/blob/c59d0ca01ae3976db7f8abdc83fcc871a3a95186/src/s2/s2coords.h#L348-L357
    s2_point = _s2_face_uv_to_xyz(face, uv)

    # Normalise XYZ S2Point vector
    # This section is part of the reference implementation but is not necessary when mapping
    # straight into lat/lon, since the normalised and unnormalised triangles used to calculate the
    # angles are geometrically similar. If anything, the normalisation process loses precision when
    # tested against the reference implementation, albeit not at a level that is important either
    # way. The code below is left for demonstration of the normalisation process.
    # norm = math.sqrt(s2_point[0]^2 + s2_point[1]^2 + s2_point[2]^2)
    # s2_point = (s2_point[0] / norm, s2_point[1] / norm, s2_point[2] / norm)

    # Map into lat/lon
    # See s2geometry/blob/c59d0ca01ae3976db7f8abdc83fcc871a3a95186/src/s2/s2latlng.h#L196-L205
    lat_rad = atan(s2_point[3], sqrt(s2_point[1]^2 + s2_point[2]^2))
    lon_rad = atan(s2_point[2], s2_point[1])

    return (rad2deg(lat_rad), rad2deg(lon_rad))
end

md"""
Convert S2 token to lat/lon.

See s2geometry/blob/c59d0ca01ae3976db7f8abdc83fcc871a3a95186/src/s2/s2cell_id.cc#L222-L239

Args:
    token: The S2 token string. Can be upper or lower case hex string.

Returns:
    The lat/lon (in degrees) tuple generated from the S2 token.

Raises:
    TypeError: If the token is not str.
    InvalidToken: If the token length is over 16.
    InvalidToken: If the token is invalid.
    InvalidCellID: If the contained cell_id is invalid.

"""
function token_to_lat_lon(token:: String)::Tuple{<:Real,<:Real}
    # Check input
    if !token_is_valid(token)
        throw(InvalidToken("Cannot decode invalid S2 token: $token"))
    end

    # Convert to cell ID and decode to lat/lon
    return cell_id_to_lat_lon(token_to_cell_id(token))
end

#
# Token canonicalisation
#

md"""
Convert S2 token to a canonicalised S2 token.

This produces a token that matches the form generated by the reference C++ implementation:

- Lower case (except 'X' below)
- No whitespace
- Trailing '0' characters stripped
- Zero cell ID represented as 'X', not 'x' or ''

Args:
    token: The S2 token string to canonicalise.

Returns:
    The canonicalised S2 token.

"""
function token_to_canonical_token(token:: String):: String
    # Convert token to lower case.
    # Note that 'X' below will be returned upper case
    token = lowercase(token)

    # Strip any surrounding whitespace
    token = strip(token)

    # Strip any trailing zeros
    token = rstrip(token, '0')

    # If empty string or "x", return "X" token
    if token in ("", "x")
        token = "X"
    end

    return token
end

#
# Validation
#

md"""
Check that a S2 cell ID is valid.

Looks for valid face bits and a trailing 1 bit in one of the correct locations.

Args:
    cell_id: The S2 cell integer to validate.

Returns:
    True if the cell ID is valid, False otherwise.

Raises:
    TypeError: If the cell_id is not int.

"""
cell_id_is_valid(cell:: Cell):: Bool = cell_id_is_valid(cell.id)

md"""
Check that a S2 cell ID is valid.

Looks for valid face bits and a trailing 1 bit in one of the correct locations.

Args:
    cell_id: The S2 cell integer to validate.

Returns:
    True if the cell ID is valid, False otherwise.

Raises:
    TypeError: If the cell_id is not int.

"""
function cell_id_is_valid(cell_id:: Integer):: Bool
    # Check for zero ID
    # This avoids overflow warnings below when 1 gets added to max uint64
    if cell_id == 0
        return false
    end

    # Check face bits
    if (cell_id >> _S2_POS_BITS) > 5
        return false
    end

    # Check trailing 1 bit is in one of the even bit positions allowed for the 30 levels, using the
    # mask: 0b0001010101010101010101010101010101010101010101010101010101010101 = 0x1555555555555555
    lowest_set_bit = cell_id & (~cell_id + 1)
    if (lowest_set_bit & 0x1555555555555555) == 0
        return false
    end

    return true  # Checks have passed, cell ID must be valid
end

md"""
Check that a S2 token is valid.

Looks for valid characters, then checks that the contained S2 cell ID is also valid. Note that
the '', 'x' and 'X' tokens are considered invalid, since the cell IDs they represent are
invalid.

Args:
    token: The S2 token string to validate.

Returns:
    True if the token is valid, False otherwise.

Raises:
    TypeError: If the token is not str.

"""
function token_is_valid(token:: String):: Bool
    # First check string with regex
    if match(r"^[0-9a-fA-F]{1,16}$", token) === nothing
        return false
    end

    # Check the contained cell ID is valid
    return cell_id_is_valid(token_to_cell_id(token))
end

#
# Level extraction functions
#

md"""
Get the level for a S2 cell ID.

See s2geometry/blob/c59d0ca01ae3976db7f8abdc83fcc871a3a95186/src/s2/s2cell_id.h#L543-L551

Args:
    cell_id: The S2 cell ID integer.

Returns:
    The level of the S2 cell ID.

Raises:
    TypeError: If the cell_id is not int.
    InvalidCellID: If the cell_id is invalid.

"""
cell_id_to_level(cell:: Cell):: Integer = cell_id_to_level(cell.id)

md"""
Get the level for a S2 cell ID.

See s2geometry/blob/c59d0ca01ae3976db7f8abdc83fcc871a3a95186/src/s2/s2cell_id.h#L543-L551

Args:
    cell_id: The S2 cell ID integer.

Returns:
    The level of the S2 cell ID.

Raises:
    TypeError: If the cell_id is not int.
    InvalidCellID: If the cell_id is invalid.

"""
function cell_id_to_level(cell_id:: Integer):: Integer
    # Check input
    if ! cell_id_is_valid(cell_id)
        throw(InvalidCellID("Cannot decode invalid S2 cell ID: $cell_id"))
    end

    # Find the position of the lowest set one bit, which will be the trailing one bit. The level is
    # given by the max level (30) minus the floored division by two of the position of the lowest
    # set bit.
    #
    # The position of the lowest set bit is found using 'count trailing zeros', which would be
    # equivalent to the C++20 function std::countr_zero() or the ctz instruction.
    lsb_pos = 0
    while cell_id != 0
        if Bool(cell_id & 1)
            break
        end
        lsb_pos += 1
        cell_id >>= 1
    end

    return _S2_MAX_LEVEL - (lsb_pos >> 1)
end

md"""
Get the level for a S2 token.

See s2geometry/blob/c59d0ca01ae3976db7f8abdc83fcc871a3a95186/src/s2/s2cell_id.h#L543-L551

Args:
    token: The S2 token string. Can be upper or lower case hex string.

Returns:
    The level of the S2 token.

Raises:
    TypeError: If the token is not str.
    InvalidToken: If the token length is over 16.
    InvalidToken: If the token is invalid.
    InvalidCellID: If the contained cell_id is invalid.

"""
function token_to_level(token:: String):: Integer
    # Check input
    if !token_is_valid(token)
        throw(InvalidToken("Cannot decode invalid S2 token: $token"))
    end

    # Convert to cell ID and get the level for that
    return cell_id_to_level(token_to_cell_id(token))
end

#
# Parent functions
#

md"""
Get the parent cell ID of a S2 cell ID.

Args:
    cell_id: The S2 cell ID integer.
    level: The parent level to get the cell ID for. Must be less than or equal to the current
        level of the provided cell ID. If unspecified, or None, the direct parent cell ID will
        be returned.

Returns:
    The parent cell ID at the specified level.

Raises:
    TypeError: If the cell_id is not int.
    InvalidCellID: If the cell_id is invalid.
    ValueError: If cell ID is already level 0 and level is None.
    ValueError: When level is not an integer, is < 0 or is > 30.
    ValueError: If level is greater than the provided cell ID level.

"""
cell_id_to_parent_cell_id(cell:: Cell, level = nothing)::Integer = cell_id_to_parent_cell_id(cell.id, level)

md"""
Get the parent cell ID of a S2 cell ID.

Args:
    cell_id: The S2 cell ID integer.
    level: The parent level to get the cell ID for. Must be less than or equal to the current
        level of the provided cell ID. If unspecified, or None, the direct parent cell ID will
        be returned.

Returns:
    The parent cell ID at the specified level.

Raises:
    TypeError: If the cell_id is not int.
    InvalidCellID: If the cell_id is invalid.
    ValueError: If cell ID is already level 0 and level is None.
    ValueError: When level is not an integer, is < 0 or is > 30.
    ValueError: If level is greater than the provided cell ID level.

"""
function cell_id_to_parent_cell_id(cell_id:: Integer, level = nothing)::Integer
    # Check input
    if ! cell_id_is_valid(cell_id)
        throw(InvalidCellID("Cannot decode invalid S2 cell ID: $cell_id"))
    end

    # Get current level of the cell ID and check it is suitable with the requested level
    current_level = cell_id_to_level(cell_id)
    if level == -1 && current_level == 0
        throw(ArgumentError("Cannot get parent cell ID of a level 0 cell ID"))
    end
    if level === nothing
        level = current_level - 1
    end

    if !(level isa Integer) || level < 0 || level > _S2_MAX_LEVEL
        throw(ArgumentError("S2 level must be integer >= 0 and <= 30 but was $level"))
    end

    if level > current_level
        throw(ArgumentError("Cannot get level $level parent cell ID of cell ID with level $current_level"))
    end
    if level == current_level
        # Requested parent level is current level, return cell ID itself
        return cell_id
    end

    # Truncate to desired level
    # This is done by finding the mask of the trailing 1 bit for the specified level, then zeroing
    # out all bits less significant than this, then finally setting the trailing 1 bit. This is
    # still necessary to do even after a reduced number of steps `required_steps` above, since each
    # step contains multiple levels that may need partial overwrite. Additionally, we need to add
    # the trailing 1 bit, which is not yet set above.
    least_significant_bit_mask = 1 << (2 * (_S2_MAX_LEVEL - level))
    cell_id = (cell_id & -least_significant_bit_mask) | least_significant_bit_mask

    return cell_id
end

md"""
Get the parent token of a S2 token.

Args:
    token: The S2 token string. Can be upper or lower case hex string.
    level: The parent level to get the token for. Must be less than or equal to the current
        level of the provided toke. If unspecified, or None, the direct parent token will be
        returned.

Returns:
    The parent token at the specified level.

Raises:
    TypeError: If the token is not str.
    InvalidToken: If the token length is over 16.
    InvalidToken: If the token is invalid.
    InvalidCellID: If the contained cell_id is invalid.
    ValueError: If token is already level 0 and level is None.
    ValueError: When level is not an integer, is < 0 or is > 30.
    ValueError: If level is greater than the provided token level.

"""
function token_to_parent_token(token:: String, level = nothing)::String
    # Check input
    if ! token_is_valid(token)
        throw(InvalidToken("Cannot decode invalid S2 token: $token"))
    end

    # Convert to cell ID and get parent and convert back to token
    return cell_id_to_token(cell_id_to_parent_cell_id(token_to_cell_id(token), level))
end

md"""
Initialise the S2 lookups in global vars _S2_LOOKUP_POS and _S2_LOOKUP_IJ.

This generates 4 variations of a 4 level deep Hilbert curve, one for each swap/invert bit
combination. This allows mapping between 8 bits (+2 orientation) of Hilbert curve position and 8
bits (+2 orientation) of I and J, and vice versa. The new orientation bits read from the mapping
tell us the base orientation of the curve segments within the next deeper level of sub-cells.

This implementation differs in structure from the reference implementation, since it is
iterative rather than recursive. The end result is the same lookup table.

See s2geometry/blob/c59d0ca01ae3976db7f8abdc83fcc871a3a95186/src/s2/s2cell_id.cc#L75-L109

"""
function __init__()
    global _S2_LOOKUP_IJ, _S2_LOOKUP_POS
    # Initialise empty lookup tables
    lookup_length = 1 << (2 * _S2_LOOKUP_BITS + 2)  # = 1024
    _S2_LOOKUP_POS = [0 for i in 1:lookup_length]
    _S2_LOOKUP_IJ = [0 for i in 1:lookup_length]

    # Generate lookups for each of the base orientations given by the swap and invert bits
    # 0-3 effectively
    for base_orientation in [0, _S2_SWAP_MASK, _S2_INVERT_MASK, _S2_SWAP_MASK | _S2_INVERT_MASK]
        # Walk the 256 possible positions within a level 4 curve. This implementation is not
        # the fastest since it does not reuse the common ancestor of neighbouring positions, but
        # is simpler to read
        for pos in 0:(4^4 - 1)  # 4 levels of sub-divisions
            ij = 0  # Has pattern iiiijjjj, not ijijijij
            orientation = base_orientation

            # Walk the pairs of bits of pos, from most significant to least, getting IJ and
            # orientation as we go
            for bit_pair_offset in 0:3
                # Bit pair is effectively the sub-cell index
                bit_pair = (pos >> ((3 - bit_pair_offset) * 2)) & 0b11

                # Get the I and J for the sub-cell index. These need to be spread into iiiijjjj
                # by inserting as bit positions 4 and 0
                ij_bits = _S2_POS_TO_IJ[orientation + 1][bit_pair + 1]
                ij = (
                    (ij << 1)  # Free up position 4 and 0 from old IJ
                    | ((ij_bits & 2) << 3)  # I bit in position 4
                    | (ij_bits & 1)  # J bit in position 0
                )

                # Update the orientation with the new sub-cell orientation
                orientation = orientation ⊻ _S2_POS_TO_ORIENTATION_MASK[bit_pair + 1]
            end

            # Shift IJ and position to allow orientation bits in LSBs of lookup
            ij <<= 2
            pos <<= 2

            # Write lookups
            _S2_LOOKUP_POS[(ij | base_orientation) + 1] = pos | orientation
            _S2_LOOKUP_IJ[(pos | base_orientation) + 1] = ij | orientation
        end
    end
end

end # module


text:
using DataArrays, Documenter

makedocs(
    modules = [DataArrays],
    clean = false,
    format = :html,
    sitename = "DataArrays.jl",
    authors = "Simon Kornblith, John Myles White, and other contributors",
    pages = [
        "Home" => "index.md",
        "Missing Data and Arrays" => "da.md",
        "Utilities" => "util.md",
    ],
)

deploydocs(
    repo = "github.com/JuliaStats/DataArrays.jl.git",
    target = "build",
    deps = nothing,
    make = nothing,
)


text:
module Nastran

import Base: show, +, -, zero

using StaticArrays

export NastranDeck, GenericNastranDeck, NastranModel
export CoordSet,get_coord,get_global_xyz
export MassCG, XYZ, Mat3x3

const XYZ = SVector{3,Float64}
const Mat3x3 = SMatrix{3,3,Float64}

include("cards.jl")
using .Cards

include("deck.jl")

include("coords.jl")
include("model.jl")
include("weights.jl")

end # module


text:
if basename(pwd()) == "aoc"
    cd("2016/19")
end

max_stars_repo_name: stellartux/advent-of-code

max_issues_repo_name: stellartux/advent-of-code

max_forks_repo_name: stellartux/advent-of-code


text:
using Suppressor

function main()
    programs = ["hello", "slow", "long", "DS"] .* ".jl"

    for program in programs
        output = @capture_err run(`time julia $program`)
        @show output
        break
    end
end

main()

max_stars_repo_name: dmolina/DaemonMode

max_issues_repo_name: dmolina/DaemonMode

max_forks_repo_name: dmolina/DaemonMode


text:
using GPLikelihoods
using AbstractGPs
using Test
using Random
using Functors
using Distributions
using StatsFuns

@testset "GPLikelihoods.jl" begin

    include("test_utils.jl")
    include("links.jl")
    @testset "likelihoods" begin
        include(joinpath("likelihoods", "bernoulli.jl"))
        include(joinpath("likelihoods", "categorical.jl"))
        include(joinpath("likelihoods", "gaussian.jl"))
        include(joinpath("likelihoods", "poisson.jl"))
        include(joinpath("likelihoods", "gamma.jl"))
        include(joinpath("likelihoods", "exponential.jl"))
    end

end


text:
using VortexDistributions
using Test

@testset "Point vortex " begin

    # unitary
    v0 = [.2 .4 1]
    w0 = PointVortex(v0)
    @test vortex_array(w0) == v0

    v1 = [.2 .4 1;0.7 1.5 -1;-.3 1.2 1]
    w1 = PointVortex(v1)
    @test vortex_array(w1) == v1

    w3 = rand_pointvortex(1000)
    v3 = vortex_array(w3)
    @test v3 == vortex_array(w3)

end

@testset "Single vortex accuracy" begin

    # simple test field
    Nx = 400; Ny = 400
    Lx = 200; Ly = Lx
    x = LinRange(-Lx/2,Lx/2, Nx+1)[1:end-1]; y = LinRange(-Ly/2,Ly/2, Ny+1)[1:end-1];
    psi0 = one.(x*y') |> complex
    psi = Torus(copy(psi0),x,y)
    xp = 1.133
    yp = 2.787
    vp = PointVortex(xp,yp,1)
    sp = ScalarVortex(vp)
    vortex!(psi,sp)
    vfound = findvortices(psi)
    @test isapprox(vfound[1].xv,xp,atol=0.01)
    @test isapprox(vfound[1].yv,yp,atol=0.01)
end

@testset "Multivortex creation and detection" begin
    @test found_near(30)
end

@testset "Periodic dipole" begin
    vp = PointVortex(.1,.3,1)
    vn = PointVortex(-.2,.7,-1)
    dip = Dipole(vp,vn)
    @test typeof(dip) <: VortexGroup
end


text:
function sample_random(shape::Tuple{Int64},redFac;calsize::Int64=0,kargs...)
  return sample_random((shape[1],1),redFac;calsize=calsize,kargs...)
end

function sample_random(shape::Tuple{Int64,Int64,Int64},redFac;calsize::Int64=0,kargs...)
  error("Not implemented")
end

"""
    sample_random(shape::Tuple{Int64,Int64},redFac::Float64;calsize::Int64=0,kargs...)

generates a random sampling pattern for an Array of size `shape` with a subsampling factor `redFac`.

# Arguments
* `shape::NTuple{N,Int64}` - size of the Array to be sampled
* `redFac::Float64`        - subsampling factor
* (`calsize::Int64=0`)     - size of the fully sampled calibration area
"""
function sample_random(shape::Tuple{Int64,Int64},redFac::Float64;calsize::Int64=0,kargs...)
  M,N = shape
  A = zeros(Int64,M,N)
  ChosenSamples = 0
  NumberSamples = floor(Int64,M*N/redFac)
  A, ChosenSamples = setcalibrationarea(A,M,N,redFac,calsize)

  chosen = (LinearIndices(A))[findall(x->x!=0, A)]
  B = abs.(A .- 1)
  selection = (LinearIndices(B))[findall(x->x!=0, B)]
  #srand(1234)
  while ChosenSamples < NumberSamples
    randnum = rand(1:length(selection))
    push!(chosen,selection[randnum])
    deleteat!(selection,randnum)
    ChosenSamples += 1
  end
  A = reshape(A,M*N,1)
  A[chosen] .= 1

  A = reshape(A,M,N)

  return chosen
end


text:
using TensorKit
using TensorOperations
using TensorKitAD
using Zygote
using KrylovKit
using LinearAlgebra
using Plots

using Revise
using statmech_tm_solver

"""
function effective_spect(A::TensorMap{ComplexSpace, 2, 1})
    _, AL = left_canonical_QR(A)
    _, AR = right_canonical_QR(A)
    map_AC, map_C = tangent_map_tn(O, AL, AR)
    W, _ = eig(map_C)
    W = diag(W.data)

    return W
end
function effective_spect(AL::TensorMap{ComplexSpace, 2, 1}, AR::TensorMap{ComplexSpace, 2, 1})
    map_AC, map_C = tangent_map_tn(O, AL, AR)
    WAC, _ = eig(map_AC)
    WAC = diag(WAC.data)

    WC, _ = eig(map_C)
    WC = diag(WC.data)

    return WAC, WC
end
"""
T = mpo_triangular_AF_ising()
T_adj = reshape(T.data, (2, 2, 2, 2))
T_adj = permutedims(conj.(T_adj), (1, 3, 2, 4))
T_adj = TensorMap(T_adj, ℂ^2*ℂ^2, ℂ^2*ℂ^2)

chi = 4
A = TensorMap(rand, ComplexF64, ℂ^chi*ℂ^2, ℂ^chi)
A = act(T, A); chi=chi*2
_, AL = left_canonical_QR(A)
_, AR = right_canonical_QR(A)
BL, BR = AL, AR

for ix in 1:3
    A = act(T, AL)
    B = act(T_adj, BL)
    chi=chi*2
    _, AL = left_canonical_QR(A)
    _, AR = right_canonical_QR(A)
    _, BL = left_canonical_QR(B)
    _, BR = right_canonical_QR(B)

    flag = 1

    for ix in 1:20
        fmap_AC, fmap_C = tangent_map(T, AL, AR)

        AC = TensorMap(rand, ComplexF64, ℂ^chi*ℂ^2, ℂ^chi)
        C = TensorMap(rand, ComplexF64, ℂ^chi, ℂ^chi)
        _, AC = eigsolve(fmap_AC, AC, 1, :LR)
        AC = AC[1]
        _, C = eigsolve(fmap_C, C, 1, :LR)
        C = C[1]

        AL1, AR1, ϵL, ϵR = calculate_ALR(AC, C)

        #WAC, WC = effective_spect(AL, AR)
        #Wmax = findmax(abs.(WAC))[1] + 0.1
        #plot(real.(WAC), imag.(WAC), seriestype=:scatter, xlims=(-Wmax, Wmax), ylims=(-Wmax, Wmax), size=(500,500))
        #plot!(real.(WC), imag.(WC), seriestype=:scatter)

        flag1 = nonherm_cost_func(T, AL1)
        if flag1 > flag
            flag = 1
            break
        else
            flag = flag1
            AL, AR = AL1, AR1
        end

        @show nonherm_cost_func(T, AL), nonherm_cost_func(T, AR)
        @show ϵL, ϵR
        @show free_energy(T, AL), free_energy(T, AR)
    end

    for ix in 1:20
        fmap_BC, fmap_C = tangent_map(T_adj, BL, BR)

        BC = TensorMap(rand, ComplexF64, ℂ^chi*ℂ^2, ℂ^chi)
        C = TensorMap(rand, ComplexF64, ℂ^chi, ℂ^chi)
        λBC, BC = eigsolve(fmap_BC, BC, 1, :LR)
        BC = BC[1]
        λC, C = eigsolve(fmap_C, C, 1, :LR)
        C = C[1]

        BL1, BR1, ϵL, ϵR = calculate_ALR(BC, C)

        #WAC, WC = effective_spect(AL, AR)
        #Wmax = findmax(abs.(WAC))[1] + 0.1
        #plot(real.(WAC), imag.(WAC), seriestype=:scatter, xlims=(-Wmax, Wmax), ylims=(-Wmax, Wmax), size=(500,500))
        #plot!(real.(WC), imag.(WC), seriestype=:scatter)

        flag1 = nonherm_cost_func(T_adj, BL1)
        if flag1 > flag
            flag = 1
            break
        else
            flag = flag1
            BL, BR = BL1, BR1
        end

        @show nonherm_cost_func(T_adj, BL), nonherm_cost_func(T_adj, BR)
        @show free_energy(T_adj, BL), free_energy(T_adj, BR)
    end

    @show ovlp(BL, act(T, AL)) / ovlp(BL, AL) |> log

    quicksave("AL_chi$(chi)", AL)
    quicksave("AR_chi$(chi)", AR)
    quicksave("BL_chi$(chi)", BL)
    quicksave("BR_chi$(chi)", BR)
end


text:
struct CompleteGraph{T<:Integer} <: LG.AbstractGraph{T}
    nv::Int
end

CompleteGraph(nv::Integer) = CompleteGraph{Int}(nv)

LG.edgetype(::CompleteGraph) = LG.Edge{Int}
LG.is_directed(::Type{<:CompleteGraph}) = false
LG.nv(g::CompleteGraph) = g.nv
LG.ne(g::CompleteGraph) = div(nv(g) * (nv(g)-1), 2)
LG.vertices(g::CompleteGraph) = 1:LG.nv(g)

LG.edges(g::CompleteGraph) = [Edge(i, j) for i in 1:nv(g)-1 for j in i+1:nv(g)]

function LG.outneighbors(g::CompleteGraph, v)
    return [v1 for v1 in vertices(g) if v1 != v]
end

LG.inneighbors(g::CompleteGraph, v) = outneighbors(g, v)

LG.has_vertex(g::CompleteGraph, v) = 1 <= v <= LG.nv(g)

function LG.has_edge(g::CompleteGraph, v1, v2)
    if !has_vertex(g, v1) || !has_vertex(g, v2)
        return false
    end
    return v1 != v2
end


text:
facts("Preprocessing") do

    fname = joinpath(dirname(@__FILE__),  "..", "data", "test_Hz19.5-testing.bdf")
    s = read_SSR(fname)

    context("Triggers") do

        context("Validation") do

            validate_triggers(s.triggers)

            s1 = deepcopy(s)
            delete!(s1.triggers, "Index")
            @fact_throws KeyError validate_triggers(s1.triggers)

            s1 = deepcopy(s)
            delete!(s1.triggers, "Code")
            @test_throws KeyError validate_triggers(s1.triggers)

            s1 = deepcopy(s)
            delete!(s1.triggers, "Duration")
            @test_throws KeyError validate_triggers(s1.triggers)

            s1 = deepcopy(s)
            s1.triggers["test"] = 1
            validate_triggers(s1.triggers)

            s1 = deepcopy(s)
            s1.triggers["Duration"] = s1.triggers["Duration"][1:4]
            @test_throws KeyError validate_triggers(s1.triggers)

            s1 = deepcopy(s)
            s1.triggers["Code"] = s1.triggers["Code"][1:4]
            @test_throws KeyError validate_triggers(s1.triggers)

        end
    end
end

facts("Referencing") do

    signals = [0 1 2] .* ones(5, 3)

    context("Remove template") do

	signals = [0 1 2] .* ones(5, 3)
	template = vec(2 * ones(5))
	@fact remove_template(signals, template) -->  [-2 -1 0] .* ones(5, 3)
	@fact_throws remove_template(signals, [1, 2, 3])

    end
    context("Reference to channel") do

	@fact rereference(signals, 3) --> [-2 -1 0] .* ones(5, 3)
	@fact rereference(signals, "C2", ["C1", "C2", "C3"]) --> [-1 0 1] .* ones(5, 3)

    end

    context("Reference to group of channels") do

	@fact rereference(signals, [1, 2, 3]) --> [-1 0 1] .* ones(5, 3)
	@fact rereference(signals, ["C2", "C1", "C3"], ["C1", "C2", "C3"]) --> [-1 0 1] .* ones(5, 3)
	@fact rereference(signals, "car", ["C1", "C2", "C3"]) --> [-1 0 1] .* ones(5, 3)
	@fact rereference(signals, "average", ["C1", "C2", "C3"]) --> [-1 0 1] .* ones(5, 3)

    end
end


text:
# Skipping MacroDefinition: imcdCanvasPutImage(_canvas,_image,_x,_y,_w,_h,_xmin,_xmax,_ymin,_ymax){if(_image->color_space==IM_RGB){if(_image->has_alpha)cdCanvasPutImageRectRGBA(_canvas,_image->width,_image->height,(unsignedchar*)_image->data[0],(unsignedchar*)_image->data[1],(unsignedchar*)_image->data[2],(unsignedchar*)_image->data[3],_x,_y,_w,_h,_xmin,_xmax,_ymin,_ymax);elsecdCanvasPutImageRectRGB(_canvas,_image->width,_image->height,(unsignedchar*)_image->data[0],(unsignedchar*)_image->data[1],(unsignedchar*)_image->data[2],_x,_y,_w,_h,_xmin,_xmax,_ymin,_ymax);}elsecdCanvasPutImageRectMap(_canvas,_image->width,_image->height,(unsignedchar*)_image->data[0],_image->palette,_x,_y,_w,_h,_xmin,_xmax,_ymin,_ymax);}
immutable imImage
    width::Cint
    height::Cint
    color_space::Cint
    data_type::Cint
    has_alpha::Cint
    depth::Cint
    line_size::Cint
    plane_size::Cint
    size::Cint
    count::Cint
    data::Ptr{Ptr{Void}}
    palette::Ptr{Clong}
    palette_count::Cint
    attrib_table::Ptr{Void}
end


text:
module TestColfuncs
	using Base.Test
	using DataArrays
	using DataFrames

	m = [1 2 3; 3 4 6]
	df = DataFrame(m)
	@assert isequal(DataFrame([2.0 3.0 4.5;], names(df)), colmeans(df))
	@assert isequal(DataFrame([4 6 9;], names(df)), colsums(df))
end


text:
module AcuteML

import EzXML.Node

# aml macro
include("types.jl")
include("xmlutils.jl")
include("@aml/@aml_parse.jl")
include("@aml/@aml_create.jl")

# io
include("io.jl")

# templating
include("templating.jl")

# main macro
export @aml
# XML types
export Node, Document

################################################################
"""
  @aml

# Type Definition
Use `@aml` macro to define a Julia type, and then the package automatically creates a xml or html associated with the defined type.

### Document Definition
* Use doc literal before the root name to define a HTML or XML document. For HTML documents root should always be "html".
```julia
@aml mutable struct Doc doc"xml_root"
# add fields (elements) here
end
```
```julia
@aml mutable struct Web doc"html"
# add fields (elements) here
end
```

### Nodes (Elements) Definition
* Specify the html/xml struct name as a string after the struct name after a space
```julia
@aml mutable struct Person "person"
# add fields (elements) here
end
```
* If the html/xml name is the same as struct name, you can use `"~"` instead
```julia
@aml mutable struct person "~"
# add fields (elements) here
end
```

### Fields Names
* Sepecify the html/xml field name as a string in front of the field after `,`
```julia
field, "study-field"
```
* If the html/xml name is the same as variable name, you can use `"~"` instead
```julia
age::UInt, "~"
```

### Attributes
* If the value is going to be an attribute put `att` before its name
```julia
id::Int64, att"~"
```

### Default Value
* You can specify the default value for an argument by using `= defVal` syntax
```julia
GPA::Float64 = 4.5, "~"
```

### Value Types
You can use Julia types or defined types for values. see and [Supported Value Types](https://aminya.github.io/AcuteML.jl/dev/supportedValueTypes/)  [Custom Value Types](https://aminya.github.io/AcuteML.jl/dev/customValueTypes/) for more information.

* If you don't specify the type of a variable, it is considered to be string for aml manipulations:
```julia
field, "study-field"
```
However, for a high performance code specify String type (`field::String, "study-field"`)

* For already `@aml` defined types, name should be the same as the defined type root name
```julia
university::University, "university"
```

* Table types are supported through PrettyTables.jl.

### Value Checking
You can define any restriction for values using functions.

* To define any restrictions for the values of one field, define a function that checks a criteria for the field value and returns Bool, and put its name after a `,` after the field name:
```julia
GPA::Float64, "~", GPAcheck
```

* To define any restrictions for multiple values of a struct, define a function that gets all the variables and checks a criteria and returns Bool, and put its name after a `,` after the struct name:
```julia
@aml mutable struct Person "person", check_course
# ...
end
```

Refer to https://aminya.github.io/AcuteML.jl/dev/valueChecking/ for some of these functions examples.

### Optional Fields
* If a field is optional, don't forget to define its type as `UN{}` (Union with Nothing), and set the default value as `nothing`.
```julia
residence::UN{String}=nothing, "residence-stay" # optional with nothing as default value
```
```julia
funds::UN{String}, "financial-funds"   # optional, but you should pass nothing manually in construction
```

### Text Nodes
If the value is going to be in a Text node:
- use `txt"index"` for non-vector field type, which `index` is an Integer that shows the positon of text node. If you give `txt""` it considers it like `txt"1"`.

```julia
textnode_single:String, txt"2"
```

- use `txt"indices"` for vector field type, which `indices` is an array index that shows the positons of the text nodes. If you give `txt""` it considers it like `txt"[:]"`

```julia
textnode_vector::Vector{String}, txt"[2:3]"
```

Note that the vector Text nodes should only be used as the last field of a struct (because possible positons for text node should be known). Alternatively, you can make non-vector separate fields with correct position in the struct definition.

### Empty Elements (Self-Closing) Definition
* Use `sc"name"` to define a self-closing (empty) element (e.g. `<rest />`)
```julia
@aml struct rest empty"~"
end
```
-------------------------------------------------------

# Example - Simple
```julia
using AcuteML

@aml mutable struct body "~"
    h1, "~"
    p::Vector{String}, "~"
end

@aml mutable struct html doc"html"
    body::body, "~"
end

b = body(h1 = "My heading", p = ["Paragraph1", "Paragraph2"])
d = html(body = b)
```

```html
julia> pprint(d)
<?xml version="1.0" encoding="UTF-8" standalone="yes"?>
<!DOCTYPE html PUBLIC "-//W3C//DTD HTML 4.0 Transitional//EN" "http://www.w3.org/TR/REC-html40/loose.dtd">
<html>
  <body>
    <h1>My heading</h1>
    <p>Paragraph1</p>
    <p>Paragraph2</p>
  </body>
 </html>
```
-------------------------------------------------------

# Example - Struct Definition

First, we define the structs using `@aml` to store the data in:

```julia
using AcuteML

# Types definition

# Person Type
@aml mutable struct Person "person", check_course
    age::UInt64, "~"
    field, "study-field"
    GPA::Float64 = 4.5, "~", GPAcheck
    courses::Vector{String}, "taken-courses"
    professors::UN{DataFrame} = nothing, "table"
    id::Int64, att"~"
    comment::UN{String} = nothing, txt"end"
end

@aml mutable struct University doc"university"
    name, att"university-name"
    people::Vector{Person}, "person"
end

```

```julia
# Value Checking Functions
GPAcheck(x) = x <= 4.5 && x >= 0

function check_course(age, field, GPA, courses, professors, id, comment)

    if field == "Mechanical Engineering"
        relevant = ["Artificial Intelligence", "Robotics", "Machine Design"]
    elseif field == "Computer Engineering"
        relevant = ["Julia", "Algorithms"]
    else
        error("study field is not known")
    end

    return any(in.(courses, Ref(relevant)))
end
```
-------------------------------------------------------

# Example - Creator

After we defined the structs, we can create instances of them by passing our data to the fields:

```julia

P1 = Person(age=24, field="Mechanical Engineering", courses = ["Artificial Intelligence", "Robotics"], id = 1, comment = "He is a genius")
P2 = Person(age=18, field="Computer Engineering", GPA=4, courses=["Julia"], id = 2)

U = University(name="Julia University", people=[P1, P2])

U.people[2].GPA=4.2 # mutability support after Doc creation

```

```julia
# An example that doesn't meet the criteria function for GPA because GPA is more than 4.5
P3 = Person(age=99, field="Macro Wizard", GPA=10, courses=["Julia Magic"], id = 3)
julia>
GPA doesn't meet criteria function
```

```html
julia> pprint(P1) # or print(P1.aml)
<person id="1">
  <age>24</age>
  <study-field>Mechanical Engineering</study-field>
  <GPA>4.5</GPA>
  <taken-courses>Artificial Intelligence</taken-courses>
  <taken-courses>Robotics</taken-courses>
  He is a genius
</person>

julia> pprint(U) # or print(U.aml)
<?xml version="1.0" encoding="UTF-8"?>
<university university-name="Julia University">
  <person id="1">
    <age>24</age>
    <study-field>Mechanical Engineering</study-field>
    <GPA>4.5</GPA>
    <taken-courses>Artificial Intelligence</taken-courses>
    <taken-courses>Robotics</taken-courses>
    He is a genius
  </person>
  <person id="2">
    <age>18</age>
    <study-field>Computer Engineering</study-field>
    <GPA>4.2</GPA>
    <taken-courses>Julia</taken-courses>
  </person>
</university>
```

P3 with Tables.jl type:
```julia
Profs1 = DataFrame(course = ["Artificial Intelligence", "Robotics"], professor = ["Prof. A", "Prof. B"] )

P3 = Person(age=24, field="Mechanical Engineering", courses = ["Artificial Intelligence", "Robotics"], professors= Profs1, id = 1)
```
```html
julia> pprint(P3)

<person id="1">
<age>24</age>
<study-field>Mechanical Engineering</study-field>
<GPA>4.5</GPA>
<taken-courses>Artificial Intelligence</taken-courses>
<taken-courses>Robotics</taken-courses>
<table>
<tr class="header">
<th style="text-align: right; ">course</th>
<th style="text-align: right; ">professor</th>
</tr>
<tr class="subheader headerLastRow">
<th style="text-align: right; ">String</th>
<th style="text-align: right; ">String</th>
</tr>
<tr>
<td style="text-align: right; ">Artificial Intelligence</td>
<td style="text-align: right; ">Prof. A</td>
</tr>
<tr>
<td style="text-align: right; ">Robotics</td>
<td style="text-align: right; ">Prof. B</td>
</tr>
</table>
</person>
```
-------------------------------------------------------

# Example - Extractor

After we defined the structs, we can automatically extract and store the data in their fields:

```julia
using AcuteML

xml = parsexml(\"\"\"
<?xml version="1.0" encoding="UTF-8"?>
<university university-name="Julia University">
  <person id="1">
    <age>24</age>
    <study-field>Mechanical Engineering</study-field>
    <GPA>4.5</GPA>
    <taken-courses>Artificial Intelligence</taken-courses>
    <taken-courses>Robotics</taken-courses>
    He is a genius
  </person>
  <person id="2">
    <age>18</age>
    <study-field>Computer Engineering</study-field>
    <GPA>4.2</GPA>
    <taken-courses>Julia</taken-courses>
  </person>
</university>
\"\"\")

# extract University
U = University(xml) # StructName(xml) extracts the data and stores them in proper format

# Now you can access all of the data by calling the fieldnames

julia>U.name
"Julia University"

# extract Person
P1 = U.people[1]

julia>P1.age
24

julia>P1.field
Mechanical Engineering

julia>P1.GPA
4.5

julia>P1.courses
["Artificial Intelligence", "Robotics"]

julia>P1.id
1

julia> P1.comment
"He is a genius"
```
"""
macro aml(expr)
    expr = macroexpand(__module__, expr) # to expand literal macros and @static
    # expr = macroexpand(@__MODULE__, expr) # for functions debuging.

    #  check if aml is used before struct
    if expr isa Expr && expr.head == :struct

        # expr.args[3] # arguments
        # args_param.args # empty
        expr.args[3], args_param, args_defaultvalue, args_type, args_var, args_name, args_function, args_literaltype, struct_name, struct_nodetype, struct_function, is_struct_mutable, args_custom_creator, args_custom_extractor, args_custom_updater, T = aml_parse(expr)

        out = aml_create(expr, args_param, args_defaultvalue, args_type, args_var, args_name, args_function, args_literaltype, struct_name, struct_nodetype, struct_function, is_struct_mutable, args_custom_creator, args_custom_extractor, args_custom_updater, T)

    # elseif expr isa Expr && expr.head == :tuple
    #     amlTypesSupport(expr)
    else
        error("Invalid usage of @aml")
    end

    return out
end
################################################################

include("../deps/SnoopCompile/precompile/precompile_AcuteML.jl")
_precompile_()

end

max_stars_repo_name: JunoLab/atom-indent-detective


text:
precompile(Ipopt.eval_jac_g_wrapper, (Int32, Ptr{Float64}, Int32, Int32, Int32, Ptr{Int32}, Ptr{Int32}, Ptr{Float64}, Ptr{Void},))
precompile(Ipopt.eval_h_wrapper, (Int32, Ptr{Float64}, Int32, Float64, Int32, Ptr{Float64}, Int32, Int32, Ptr{Int32}, Ptr{Int32}, Ptr{Float64}, Ptr{Void},))
precompile(Ipopt.createProblem, (Int64, Array{Float64, 1}, Array{Float64, 1}, Int64, Array{Float64, 1}, Array{Float64, 1}, Int64, Int64, Function, Function, Function, Function, Function,))
precompile(Ipopt.eval_f_wrapper, (Int32, Ptr{Float64}, Int32, Ptr{Float64}, Ptr{Void},))
precompile(Ipopt.eval_g_wrapper, (Int32, Ptr{Float64}, Int32, Int32, Ptr{Float64}, Ptr{Void},))
precompile(Ipopt.eval_grad_f_wrapper, (Int32, Ptr{Float64}, Int32, Ptr{Float64}, Ptr{Void},))
precompile(Ipopt.call, (Type{Ipopt.IpoptProblem}, Ptr{Void}, Int64, Int64, Function, Function, Function, Function, Function,))
precompile(Ipopt.call, (Array{Any, 1}, Type{Ipopt.IpoptMathProgModel},))
precompile(Ipopt.optimize!, (Ipopt.IpoptMathProgModel,))
precompile(Ipopt.call, (Array{Any, 1}, Type{Ipopt.IpoptSolver},))
precompile(Ipopt.addOption, (Ipopt.IpoptProblem, ASCIIString, ASCIIString,))
precompile(Ipopt.NonlinearModel, (Ipopt.IpoptSolver,))
precompile(Ipopt.status, (Ipopt.IpoptMathProgModel,))
precompile(Ipopt.solveProblem, (Ipopt.IpoptProblem,))
precompile(Ipopt.call, (Type{Ipopt.IpoptMathProgModel},))
precompile(Ipopt.addOption, (Ipopt.IpoptProblem, ASCIIString, Int64,))
precompile(Ipopt.getobjval, (Ipopt.IpoptMathProgModel,))
precompile(Ipopt.setwarmstart!, (Ipopt.IpoptMathProgModel, Array{Float64, 1},))
precompile(Ipopt.getsolution, (Ipopt.IpoptMathProgModel,))


text:
using YaoBlocksQASM, DocThemeIndigo
using Documenter
indigo = DocThemeIndigo.install(YaoBlocksQASM)

makedocs(;
    modules = [YaoBlocksQASM],
    authors = "Arsh Sharma",
    repo = "https://github.com/QuantumBFS/YaoBlocksQASM.jl/blob/{commit}{path}#{line}",
    sitename = "YaoBlocksQASM.jl",
    format = Documenter.HTML(;
        prettyurls = get(ENV, "CI", "false") == "true",
        canonical = "https://QuantumBFS.github.io/YaoBlocksQASM.jl",
        assets = String[],
    ),
    pages = ["Quickstart" => "index.md"],
)

deploydocs(; repo = "github.com/QuantumBFS/YaoBlocksQASM.jl", devbranch = "main")

max_stars_repo_name: Sov-trotter/YaoBlocksQASM

max_issues_repo_name: Sov-trotter/YaoBlocksQASM

max_forks_repo_name: Sov-trotter/YaoBlocksQASM


text:
# *********************************************************************************
# REopt, Copyright (c) 2019-2020, Alliance for Sustainable Energy, LLC.
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without modification,
# are permitted provided that the following conditions are met:
#
# Redistributions of source code must retain the above copyright notice, this list
# of conditions and the following disclaimer.
#
# Redistributions in binary form must reproduce the above copyright notice, this
# list of conditions and the following disclaimer in the documentation and/or other
# materials provided with the distribution.
#
# Neither the name of the copyright holder nor the names of its contributors may be
# used to endorse or promote products derived from this software without specific
# prior written permission.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
# ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
# WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
# IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT,
# INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
# BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
# DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
# LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE
# OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED
# OF THE POSSIBILITY OF SUCH DAMAGE.
# *********************************************************************************
"""
    add_financial_results(m::JuMP.AbstractModel, p::REoptInputs, d::Dict; _n="")

Adds the Financial results to the dictionary passed back from `run_reopt` using the solved model `m` and the `REoptInputs` for node `_n`.
Note: the node number is an empty string if evaluating a single `Site`.

Financial results:
- `lcc_us_dollars` Optimal lifecycle cost
- `net_capital_costs_plus_om_us_dollars` Capital cost for all technologies plus present value of operations and maintenance over anlaysis period
- `net_capital_costs` Net capital costs for all technologies, in present value, including replacement costs and incentives.
"""
function add_financial_results(m::JuMP.AbstractModel, p::REoptInputs, d::Dict; _n="")
    r = Dict{String, Any}()
    r["lcc_us_dollars"] = round(value(m[Symbol("Costs"*_n)]) + 0.0001 * value(m[Symbol("MinChargeAdder"*_n)]))
    r["net_capital_costs_plus_om_us_dollars"] = round(
        value(m[Symbol("TotalTechCapCosts"*_n)] + m[Symbol("TotalStorageCapCosts"*_n)]) +
        value(m[Symbol("TotalPerUnitSizeOMCosts"*_n)]) * (1 - p.owner_tax_pct), digits=0
    )
    r["net_capital_costs"] = round(value(m[Symbol("TotalTechCapCosts"*_n)] + m[Symbol("TotalStorageCapCosts"*_n)]),
                                   digits=2)
    d["Financial"] = r
    nothing
end

max_stars_repo_name: NLaws/REoptLite

max_issues_repo_name: NLaws/REoptLite

max_forks_repo_name: NLaws/REoptLite