text:
@testset "Primal Dual Traj" begin

    # Test PrimalDualTraj
    N = 10
    dt = 0.2
    p = 3
    model = UnicycleGame(p=p)
    probsize = ProblemSize(N, model)
    pdtraj = PrimalDualTraj(probsize, dt)
    @test size(pdtraj.pr)[1] == N
    @test size(state(pdtraj.pr.data[1]))[1] == model.n
    @test size(control(pdtraj.pr.data[1]))[1] == model.m
    @test pdtraj.pr.data[1].dt == dt
    @test size(pdtraj.du)[1] == p
    @test size(pdtraj.du[1])[1] == N-1
    @test size(pdtraj.du[1][1])[1] == model.n

    # Test init_traj!
    T = Float64
    x0 = rand(SVector{model.n,T})
    init_traj!(pdtraj, x0=x0, f=ones, amplitude=10.0)
    @test state(pdtraj.pr[1]) == x0
    @test state(pdtraj.pr[2]) == 10*ones(model.n)
    @test control(pdtraj.pr[1]) == 10*ones(model.m)
    @test pdtraj.du[1][1] == 10*ones(model.n)
    @test pdtraj.du[1][1] == 10*ones(model.n)

    # Test set_traj!
    n = model.n
    m = model.m
    p = model.p
    core = NewtonCore(probsize)
    x0 = rand(SVector{model.n,T})
    Δpdtraj = PrimalDualTraj(probsize, dt)
    Δtraj = ones(n*(N-1)+m*(N-1)+n*p*(N-1))
    init_traj!(Δpdtraj, x0=x0, f=ones, amplitude=10.0)
    set_traj!(core, Δpdtraj, Δtraj)

    @test state(Δpdtraj.pr[1]) == x0
    @test state(Δpdtraj.pr[2]) == ones(n)
    @test state(Δpdtraj.pr[end]) == ones(n)
    @test control(Δpdtraj.pr[1]) == ones(m)
    @test control(Δpdtraj.pr[end-1]) == ones(m)
    @test Δpdtraj.du[1][1] == ones(n)
    @test Δpdtraj.du[1][end] == ones(n)
    @test Δpdtraj.du[end][end] == ones(n)

    Δtraj = rand(n*(N-1)+m*(N-1)+n*p*(N-1))
    init_traj!(Δpdtraj, x0=x0, f=ones, amplitude=10.0)
    set_traj!(core, Δpdtraj, Δtraj)

    @test state(Δpdtraj.pr[1]) == x0
    @test state(Δpdtraj.pr[2]) == Δtraj[horizontal_idx(core, stampify(:x, 1, 2))]
    @test state(Δpdtraj.pr[end]) == Δtraj[horizontal_idx(core, stampify(:x, 1, N))]
    @test control(Δpdtraj.pr[1])[model.pu[1]] == Δtraj[horizontal_idx(core, stampify(:u, 1, 1))]
    @test control(Δpdtraj.pr[1])[model.pu[2]] == Δtraj[horizontal_idx(core, stampify(:u, 2, 1))]
    @test control(Δpdtraj.pr[2])[model.pu[1]] == Δtraj[horizontal_idx(core, stampify(:u, 1, 2))]
    @test control(Δpdtraj.pr[2])[model.pu[2]] == Δtraj[horizontal_idx(core, stampify(:u, 2, 2))]
    @test control(Δpdtraj.pr[2])[model.pu[3]] == Δtraj[horizontal_idx(core, stampify(:u, 3, 2))]
    @test control(Δpdtraj.pr[end-1])[model.pu[1]] == Δtraj[horizontal_idx(core, stampify(:u, 1, N-1))]
    @test Δpdtraj.du[1][1] == Δtraj[horizontal_idx(core, stampify(:λ, 1, 1))]
    @test Δpdtraj.du[1][end] == Δtraj[horizontal_idx(core, stampify(:λ, 1, N-1))]
    @test Δpdtraj.du[end][end] == Δtraj[horizontal_idx(core, stampify(:λ, 3, N-1))]

    # Test get_traj!
    N = 10
    dt = 0.1
    p = 3
    model = UnicycleGame(p=p)
    probsize = ProblemSize(n,model)
    core = NewtonCore(probsize)
    Δpdtraj = PrimalDualTraj(probsize, dt)
    Δtraj0 = rand(probsize.S)
    Δtraj1 = deepcopy(Δtraj0)
    Δtraj2 = deepcopy(Δtraj0)

    get_traj!(core, Δpdtraj, Δtraj0)
    Δtraj0
    Δtraj1

    @test !(Δtraj0 == Δtraj1)
    set_traj!(core, Δpdtraj, Δtraj1)
    get_traj!(core, Δpdtraj, Δtraj0)
    @test Δtraj0 == Δtraj1

    # Test update_traj!
    n = model.n
    m = model.m
    p = model.p
    x0 = rand(SVector{model.n,T})
    target = PrimalDualTraj(probsize, dt)
    source = PrimalDualTraj(probsize, dt)
    Δ = PrimalDualTraj(probsize, dt)
    init_traj!(target, x0=x0, f=ones, amplitude=0.0)
    init_traj!(source, x0=x0, f=ones, amplitude=10.0)
    init_traj!(Δ, x0=x0, f=ones, amplitude=100.0)
    α = 0.5
    update_traj!(target, source, α, Δ)

    @test state(target.pr[1]) == x0
    @test state(target.pr[2]) == 60*ones(n)
    @test state(target.pr[end]) == 60*ones(n)
    @test control(target.pr[1]) == 60*ones(m)
    @test control(target.pr[end-1]) == 60*ones(m)
    @test target.du[1][1] == 60*ones(n)
    @test target.du[1][end] == 60*ones(n)
    @test target.du[end][end] == 60*ones(n)

    # Test Δ_step
    T = Float64
    N = 10
    dt = 0.2
    p = 3
    model = UnicycleGame(p=p)
    probsize = ProblemSize(N, model)
    pdtraj = PrimalDualTraj(probsize, dt)
    n = model.n
    m = model.m
    x0 = 1e3*ones(SVector{model.n,T})
    Δ = PrimalDualTraj(probsize, dt)
    init_traj!(pdtraj, x0=x0, f=ones, amplitude=10.0)
    α = 0.5
    @test Δ_step(pdtraj, α) == 10.0*α

end


text: """rmp"""

ext: jl


text:
using FunctionalCollections

struct JuliaBackend <: ArrayBackend
end

const ObjectID = UInt64

function reify_value(val::Complex)
    if real(val) == 0
        :( $(imag(val)) * im ), Nullable{Symbol}()
    elseif imag(val) == 0
        real(val), Nullable{Symbol}()
    else
        :( $(real(val)) + $(imag(val)) * im ), Nullable{Symbol}()
    end
end

function reify_value(val::AbstractArray)
    if ndims(val) == 0
        return val[1], Nullable{Symbol}()
    else
        name = gensym()
        ast_node = Symbol(name)
        return ast_node, Nullable{Symbol}(ast_node)
    end
end

function reify_value(val::Number)
    val, Nullable{Symbol}()
end

function reify_value(val)
    name = gensym()
    ast_node = Symbol(name)
    return ast_node, Nullable{Symbol}(ast_node)
end

function reify_leaf(identifiers, val)

    val_id = object_id(val)
    if haskey(identifiers, val_id)
        ast_node, _ = identifiers[val_id]
    else
        ast_node, name = reify_value(val)
        if !isnull(name)
            identifiers = assoc(identifiers, val_id, get(name) => val)
        end
    end

    identifiers, ast_node
end

function reify(identifiers, aexpr::LazyArray)
    reify_leaf(identifiers, aexpr.value)
end

function reify(identifiers, aexpr::BackendArray)
    reify_leaf(identifiers, from_backend(aexpr))
end

function reify(identifiers, aexpr::Unknown)
    reify_leaf(identifiers, aexpr)
end

function unpack(x, i)
    x[i]
end

function repeat_outer(x, outer)
    repeat(x, outer=outer)
end

# Julia's reshape only accepts tuples
function reshape(x::AbstractArray, dims::AbstractArray)
    reshape(x, dims...)
end

# Zero-dimensional arrays are reified as scalars, and Julia doesn't have a reshape() method
# for scalars, so we have to provide our own.
function reshape(x::Number, dims...)
    reshape([x], dims...)
end

function random_uniform(state, dims...)
    state_copy = copy(state)
    res = rand(state_copy, Float64, dims...)
    res, state_copy
end

function random_normal(state, dims...)
    state_copy = copy(state)
    res = randn(state_copy, dims...)
    res, state_copy
end

function reify(identifiers, aexpr::CallNode)
    head = aexpr.func
    args = aexpr.args
    identifiers, arg_asts = reify(identifiers, args)
    identifiers, :( $head($(arg_asts...)) )
end

function reify{T}(identifiers, arr::Array{T, 1})
    result = Array{Any, 1}(length(arr))
    for (i, elem) in enumerate(arr)
        identifiers, elem_ast = reify(identifiers, elem)
        result[i] = elem_ast
    end
    identifiers, result
end

function reify{K, V}(identifiers, pair::Pair{K, V})
    fst, snd = pair
    identifiers, fst = reify(identifiers, fst)
    identifiers, snd = reify(identifiers, snd)
    identifiers, (fst => snd)
end

function reify{K, V}(identifiers, dict::Dict{K, V})
    result = Dict{Any, Any}()
    for (i, pair) in enumerate(dict)
        identifiers, pair_ast = reify(identifiers, pair)
        k, v = pair_ast
        result[k] = v
    end
    identifiers, result
end

struct JuliaArrayFunction <: BackendArrayFunction
    func :: Function
    known_arguments :: Array{AbstractArray, 1}
    mutated_parameters_indices :: Set{Int64}
end

function (cafunc::JuliaArrayFunction)(args...)

    # The parameters being mutated must be backend arrays.
    # Technically, this check should be done in a compiled array function in any backend
    if length(cafunc.mutated_parameters_indices) > 0
        for (i, arg) in enumerate(args)
            if in(i, cafunc.mutated_parameters_indices)
                @assert isa(arg, BackendValue) && array_backend(arg) == JuliaBackend()
            end
        end
    end

    unwrapped_args = map(unwrap, args)

    # Cannot call `cafunc.func` directly because of Julia's restrictions
    # on dynamic method definition.
    # See "applicable method may be too new" error, and
    # https://docs.julialang.org/en/latest/manual/methods.html#Redefining-Methods-1
    results = Base.invokelatest(cafunc.func, cafunc.known_arguments..., unwrapped_args...)

    wrapped_results = map(wrap, results)
    length(wrapped_results) == 1 ? wrapped_results[1] : wrapped_results
end

function compile_array_function(::JuliaBackend, afunc::ArrayFunction)

    serialized_aliases = serialize(afunc)

    # Replace known/unknown arrays in the expressions to be returned with temporary identifiers
    body_identifiers = PersistentHashMap{ObjectID, Any}()
    body_identifiers, reified_aliases = reify(body_identifiers, serialized_aliases)
    identifiers, reified_parameters = reify(body_identifiers, afunc.parameters)
    identifiers, reified_mutates = reify(identifiers, afunc.mutates)
    identifiers, reified_returns = reify(identifiers, afunc.returns)

    # At this point we do not need object IDs anymore,
    # a mapping from symbols to values will suffice
    body_identifiers = Dict(pair[2] for pair in body_identifiers)
    identifiers = Dict(pair[2] for pair in identifiers)

    # Since we cannot create a context for the compiled function,
    # known arrays will have to be passed explicitly
    known_parameter_map = filter((sym, val) -> !isa(val, Unknown), identifiers)
    # According to docs, `keys` and `values` will have the same order
    known_parameter_symbols = collect(keys(known_parameter_map))
    known_arguments = collect(values(known_parameter_map))

    # Identify unused parameters and replace them with placeholders
    used_in_expressions = Set(
        filter(unknown_sym -> haskey(body_identifiers, unknown_sym), reified_parameters))
    used_in_mutations = Set(values(reified_mutates))
    used_unknown_parameters = union(used_in_expressions, used_in_mutations)
    unknown_parameter_symbols = [
        in(unknown_sym, used_unknown_parameters) ? unknown_sym : :_
        for unknown_sym in reified_parameters]

    function_name = gensym("array_function")

    # Construct the function
    expr = quote
        function $function_name(
                $(known_parameter_symbols...),
                $(unknown_parameter_symbols...))

            # Assignment to aliases
            #$(((if length(targets) == 1
            #        :( $(targets[1]) = $expr )
            #    else
            #        :( ($(targets...),) = $expr )
            #    end)
            #    for (targets, expr) in reified_aliases)...)
            $((:( const $target = $expr )
                for (target, expr) in reified_aliases)...)

            # Copying to mutated parameters
            $((:( copy!($target, $alias) )
                for (alias, target) in reified_mutates)...)

            # Returns
            $(if length(reified_returns) > 0
                :( ($(reified_returns...),) )
            end)
        end
    end

    #println("Function to compile:")
    #println("-----------------------------------")
    #println(expr)
    #println("-----------------------------------")

    func = eval(expr)

    JuliaArrayFunction(func, map(unwrap, known_arguments), mutated_parameters_indices(afunc))
end

struct JuliaBackendRandomState <: BackendRandomState
    random_state :: MersenneTwister
end

struct JuliaBackendArray <: BackendArray
    array :: AbstractArray
end

import Base: size, ndims, eltype

eltype(arr::JuliaBackendArray) = eltype(arr.array)
size(arr::JuliaBackendArray) = size(arr.array)
size(arr::JuliaBackendArray, dim) = size(arr.array, dim)
ndims(arr::JuliaBackendArray) = ndims(arr.array)

array_backend(::JuliaBackendRandomState) = JuliaBackend()
array_backend(::JuliaBackendArray) = JuliaBackend()
array_backend(::JuliaArrayFunction) = JuliaBackend()

unwrap(value::Number) = value
unwrap(value::AbstractArray) = value
unwrap(value::JuliaBackendRandomState) = value.random_state
unwrap(value::JuliaBackendArray) = value.array
unwrap(value::LazyArray) = value.value

wrap(value::Number) = value
wrap(value::AbstractArray) = JuliaBackendArray(value)
wrap(value::MersenneTwister) = JuliaBackendRandomState(value)

function random_state(::JuliaBackend; seed=nothing)
    if seed === nothing
        JuliaBackendRandomState(MersenneTwister())
    else
        JuliaBackendRandomState(MersenneTwister(seed))
    end
end

empty_array(::JuliaBackend, dtype, shape...) = JuliaBackendArray(Array(dtype, shape...))

to_backend(::JuliaBackend, arr::JuliaBackendArray) = arr
to_backend(::JuliaBackend, arr::AbstractArray) = JuliaBackendArray(arr)
from_backend(arr::JuliaBackendArray) = arr.array


text:
using LiquidDSP
using Base.Test
import DSP

for (Th, Tx) in [(Float32, Float32), (Float32, Complex64), (Complex64, Complex64)]
    # Create taps & signal
    h    = DSP.digitalfilter(DSP.Lowpass(0.5), DSP.FIRWindow(DSP.hamming(21)))
    x    = rand(Tx, 100)

    # Convert to complex if necessary
    h += Tx<:Complex ? im : 0

    # Create filter objects
    ff_d = DSP.FIRFilter(h, 1//2)
    ff_l = FIRDecim(Tx, 2, h)

    # Test gettaps
    h_l = gettaps(ff_l)
    @test_approx_eq h_l h

    # Test case where Tx assumed to be the same as Th
    if Tx == Th
        ff_l = FIRDecim(2, h)
    end

    # Test filtering
    y_d  = DSP.filt(ff_d, x)
    y_l  = execute(ff_l, x)
    @test_approx_eq y_l y_d

    # Test reset
    reset!(ff_l)
    y_l  = execute(ff_l, x)
    @test_approx_eq y_l y_d

    # TODO: This doesn't appear to work, not testing for now
    # # Make sure create with automatic taps doesn't throw an error
    # delay = 11
    # atten = 60
    # ff_l  = FIRDecim(Tx, 2, atten, delay)
    # # Check for execption when delay is 0
    # @test_throws ArgumentError FIRDecim(Tx, 2, 60, 0)

    # Make sure print doesn't throw an error
    print(ff_l)

    # Test destroy
    destroy(ff_l)
    @test_throws ErrorException execute(ff_l, x)
end

max_stars_repo_name: JayKickliter/LiquidDSP

max_issues_repo_name: JayKickliter/LiquidDSP

max_forks_repo_name: JayKickliter/LiquidDSP


text:
getptr(ptr::AbstractVector, i::Integer) = ptr[i]:(ptr[i+1]-1)

"""
    csr_to_csc(nrow, ncol, nz, rowptr, colval)

Convert the sparse matrix CSR format to CSC format and it returns `new_nz`, `colptr`, `rowval`.

# Arguments

- `nrow::Integer`: Number of rows
- `ncol::Integer`: Number of columns.
- `nz::AbstractVector`: Nonzero elements.
- `rowptr::AbstractVector`: Row pointers.
- `colval::AbstractVector`: Column indices.

Ref: https://github.com/scipy/scipy/blob/3b36a574dc657d1ca116f6e230be694f3de31afc/scipy/sparse/sparsetools/csr.h#L380
"""
function csr_to_csc(nrow::Integer, ncol::Integer, nz::AbstractVector, rowptr::AbstractVector{T}, colval::AbstractVector{T}) where {T<:Integer}
    nnz = rowptr[end] - 1
    new_nz = similar(nz)
    rowval = similar(colval)

    # calculate number of nonzero elements in each column, O(nnz)
    colptr = zeros(T, ncol+1)
    for i in 1:nnz
        colptr[colval[i]] += one(T)
    end

    # calculate offsets, O(ncol)
    cumsum = one(T)
    for j in 1:ncol
        temp = colptr[j]
        colptr[j] = cumsum
        cumsum += temp
    end
    colptr[end] = nnz

    # calculate index, O(nnz)
    for i in 1:nrow
        for jj in getptr(rowptr, i)
            j = colval[jj]
            dest = colptr[j]

            rowval[dest] = i
            new_nz[dest] = nz[jj]

            colptr[j] += 1
        end
    end

    # calculate offsets, O(ncol)
    last = one(T)
    for j in 1:(ncol+1)
        colptr[j], last = last, colptr[j]
    end
    return new_nz, colptr, rowval
end


text:
function contains(string,token)
    return occursin(token,string)
end

function is_file_path_ok(path_to_file::String)

    if (isfile(path_to_file) == false)
        throw(ArgumentError("$(path_to_file) does not exist."))
    end
end

function is_dir_path_ok(path_to_file::String)

    # get the dir of the path -
    dir_name = dirname(path_to_file)

    # check, is this a legit dir?
    if (isdir(dir_name) == false)
        throw(ArgumentError("$(dir_name) does not exist."))
    end
end

function write_file_to_path(path_to_file::String,buffer::Array{String,1})

    # check, is this a legit dir?
    is_dir_path_ok(path_to_file)

    # write -
    open("$(path_to_file)", "w") do f

        for line_item in buffer
            write(f,"$(line_item)\n")
        end
    end
end

function read_file_from_path(path_to_file::String)::Array{String,1}

    # is this mapping file path legit?
    is_file_path_ok(path_to_file)

    # initialize -
    buffer = String[]

    # Read in the file -
    open("$(path_to_file)", "r") do file
        for line in eachline(file)
            +(buffer,line)
        end
    end

    # return -
    return buffer
end

function extract_section(file_buffer_array::Array{String,1},start_section::String,end_section::String)

    # initialize -
    section_buffer = String[]

    # find the SECTION START AND END -
    section_line_start = 1
    section_line_end = 1
    for (index,line) in enumerate(file_buffer_array)
        if (occursin(start_section,line) == true)
            section_line_start = index
        elseif (occursin(end_section,line) == true)
            section_line_end = index
        end
    end

    for line_index = (section_line_start+1):(section_line_end-1)
        line_item = file_buffer_array[line_index]
        push!(section_buffer,line_item)
    end

    # return -
    return section_buffer
end

function build_mapping_dictionary(path_to_mapping_file::String;record_delim::String="=",field_delim::String=",")::Dict{String,Mapping}

    # check the file -
    is_file_path_ok(path_to_mapping_file)

    # initalize -
    mapping_dictionary = Dict{String,Mapping}()

    # ok, read -
    file_record_array = read_file_from_path(path_to_mapping_file)

    # process the record -
    for record in file_record_array

        # ignore comments -
        if (occursin("//",record) == false)

            # create map entry -
            mapping = Mapping()

            # create set -
            token_set = Set{String}()

            # record -
            token_array = split(record,record_delim)

            # the key is the first entry -
            key = token_array[1]
            mapping.key = key

            # push the tokens (genes, ecnumbers etc) into set -
            token_array = split(token_array[2],field_delim)
            for token in token_array

                # grab -
                push!(token_set,token)
            end

            # set the value -
            mapping.value = token_set

            # add entry to dictionary -
            mapping_dictionary[key] = mapping
        end
    end

    # return -
    return mapping_dictionary
end

max_stars_repo_name: varnerlab/CBModelTools

max_issues_repo_name: varnerlab/CBModelTools

max_forks_repo_name: varnerlab/CBModelTools


text:
module PokerHands

export chances

using Cards
using Combinatorics
using Random
using StaticArrays

Base.isless(c1::Card, c2::Card) = Cards.rank(c1) < Cards.rank(c2)

const VALUES = Dict(
    'A' => 14,
    'K' => 13,
    'Q' => 12,
    'J' => 11,
    'T' => 10,
    '9' => 9,
    '8' => 8,
    '7' => 7,
    '6' => 6,
    '5' => 5,
    '4' => 4,
    '3' => 3,
    '2' => 2
)

const SUITS = Dict('c'=>0, 'd'=>1, 'h'=>2, 's'=>3)

const POKER_HANDS =
    (:highcard, :pair, :twopair, :set, :straight, :flush, :full, :poker, :straightflush)

"""
    parse_hand(s)

Parse a four-element string `s` as two-element `StaticVector`.
"""
parse_hand(s) = SVector(Card(VALUES[s[1]], SUITS[s[2]]), Card(VALUES[s[3]], SUITS[s[4]]))

"""
    type_value(s)

Get the poker hand value of `Symbol` `s`, where `s` can be in `$POKER_HANDS`, `:highcard`
returns `1` and `:straightflush` returns `9`.
"""
type_value(s) = findfirst(h -> h === s, POKER_HANDS)

"""
    PokerHand

A `struct` to store a poker hand of five cards including its type, which can be in
`$POKER_HANDS`.
"""
struct PokerHand
    type::Symbol
    cards::NTuple{5, Card}
end

type_value(p::PokerHand) = type_value(p.type)

"""
    PokerHand(cards)

Determine the type of poker hand for 5 cards `cards` and return a `PokerHand` of
respective type.
"""
function PokerHand(cards)
    sorted_cards = sort(cards, rev=true)
    sorted_ranks = Cards.rank.(sorted_cards)
    isflush = length(unique(Cards.suit.(cards))) == 1
    is_ace_straight = sorted_ranks == SVector{5}(14,5,4,3,2)
    isstraight = all(==(-1), diff(sorted_ranks)) || is_ace_straight
    if is_ace_straight
        sorted_cards = SVector{5}(circshift(sorted_cards, -1))
    end

    if isflush && isstraight
        return PokerHand(:straightflush, Tuple(sorted_cards))
    elseif isflush
        return PokerHand(:flush, Tuple(sorted_cards))
    elseif isstraight
        return PokerHand(:straight, Tuple(sorted_cards))
    end

    ucards = length(unique(sorted_ranks))
    if ucards == 2
        c1 = count_equal(sorted_ranks, 1)
        if c1 == 1
            return PokerHand(:poker, Tuple(circshift(sorted_cards, -1)))
        elseif c1 == 2
            return PokerHand(:full, Tuple(circshift(sorted_cards, -2)))
        elseif c1 == 3
            return PokerHand(:full, Tuple(sorted_cards))
        elseif c1 == 4
            return PokerHand(:poker, Tuple(sorted_cards))
        end
    elseif ucards == 3
        c1 = count_equal(sorted_ranks, 1)
        c2 = count_equal(sorted_ranks, c1 + 1)
        c3 = 5 - c1 - c2
        if c1 == 3
            return PokerHand(:set, Tuple(sorted_cards))
        elseif c2 == 3
            return PokerHand(
                :set,
                (
                    sorted_cards[2],
                    sorted_cards[1],
                    sorted_cards[3],
                    sorted_cards[4],
                    sorted_cards[5],
                ),
            )
        elseif c3 == 3
            return PokerHand(:set, Tuple(circshift(sorted_cards, -2)))
        elseif c1 == 1
            return PokerHand(:twopair, Tuple(circshift(sorted_cards, -1)))
        elseif c2 == 1
            return PokerHand(
                :twopair,
                (
                    sorted_cards[1],
                    sorted_cards[2],
                    sorted_cards[4],
                    sorted_cards[5],
                    sorted_cards[3],
                ),
            )
        else
            return PokerHand(:twopair, Tuple(sorted_cards))
        end
    elseif ucards == 4
        if sorted_ranks[1] == sorted_ranks[2]
            return PokerHand(:pair, Tuple(sorted_cards))
        elseif sorted_ranks[2] == sorted_ranks[3]
            return PokerHand(
                :pair,
                (
                    sorted_cards[2],
                    sorted_cards[3],
                    sorted_cards[1],
                    sorted_cards[4],
                    sorted_cards[5],
                ),
            )
        elseif sorted_ranks[3] == sorted_ranks[4]
            return PokerHand(
                :pair,
                (
                    sorted_cards[3],
                    sorted_cards[4],
                    sorted_cards[1],
                    sorted_cards[2],
                    sorted_cards[5],
                ),
            )
        else
            return PokerHand(:pair, Tuple(circshift(sorted_cards, -3)))
        end
    else
        return PokerHand(:highcard, Tuple(sorted_cards))
    end
end

"""
    count_equal(x, i)

Count how many elements in `x` are equal to its `i`th element.
`x` is assumed to be sorted and elements before `x[i]` are not considered.
"""
function count_equal(x, i)
    j = i + 1
    while x[j] == x[i]
        j += 1
    end
    return j - i
end

"""
    play(p, q)

Compares two `PokerHand`s `p` and `q` and returns `1` if `p` wins, `2` if `q` wins and `0`
for a split pot.
"""
function play(p, q)
    pv = type_value(p)
    qv = type_value(q)
    if pv < qv
        return 2
    elseif qv < pv
        return 1
    end
    for i in 1:5
        pr = Cards.rank(p.cards[i])
        qr = Cards.rank(q.cards[i])
        if pr < qr
            return 2
        elseif qr < pr
            return 1
        end
    end
    return 0
end

"""
    get_best_hand(cards)

Select the five cards of `cards` that provide the highest ranked `PokerHand` and return it.
"""
function get_best_hand(cards)
    current_hand = PokerHand(
        :highcard, (Card(2, 0), Card(3, 1), Card(4, 2), Card(5, 3), Card(7, 0))
    )
    for card in combinations(cards, 5)
        next_hand = PokerHand(card)
        if play(next_hand, current_hand) == 1
            current_hand = next_hand
        end
    end
    return current_hand
end

"""
    play(hand, others, flop)

Compares the `hand` to all hands in `others` considering the best possible poker hand
that can be built using the cards in `flop`.
Returns `1` if hand wins, `0` if `hand` splits the pot with any other hand and else `2`.
"""
function play(hand, others, flop)
    my_poker_hand = get_best_hand(SVector{7}(hand..., flop...))
    split = false
    for cards in others
        other_poker_hand = get_best_hand(SVector{7}(cards..., flop...))
        r = play(my_poker_hand, other_poker_hand)
        if r == 2
            return 2
        elseif r == 0
            split = true
        end
    end
    return split ? 0 : 1
end

"""
    function chances(hand, p=8; n=1000) -> (win, split)

Estimates the probabilites of `hand` winning against `p` other players, by simulating `n`
random games.
"""
function chances(hand, p=8; n=10000)
    all_cards = collect(Cards.deck)
    setdiff!(all_cards, hand)
    w, s = 0, 0
    for i in 1:n
        inds = randperm(50)
        others = Tuple(SVector{2}(all_cards[inds[2i-1:2i]]) for i in 1:p)
        flop = SVector{5}(all_cards[inds[2p+1:2p+5]])
        r = play(hand, others, flop)
        w += r==1
        s += r==0
    end
    return round((w / n) * 100, digits = 2), round((s / n) * 100, digits = 2)
end
chances(s::AbstractString, p=8; n=10000) = chances(parse_hand(s), p; n=n)

end


text:
"""
    AbstractProjectiveSpace{𝔽} <: AbstractDecoratorManifold{𝔽}

An abstract type to represent a projective space over `𝔽` that is represented isometrically
in the embedding.
"""
abstract type AbstractProjectiveSpace{𝔽} <: AbstractDecoratorManifold{𝔽} end

@doc raw"""
    ProjectiveSpace{n,𝔽} <: AbstractProjectiveSpace{𝔽}

The projective space $𝔽ℙ^n$ is the manifold of all lines in $𝔽^{n+1}$.
The default representation is in the embedding, i.e. as unit norm vectors in
$𝔽^{n+1}$:
````math
𝔽ℙ^n := \bigl\{ [p] ⊂ 𝔽^{n+1} \ \big|\ \lVert p \rVert = 1, λ ∈ 𝔽, |λ| = 1, p ∼ p λ \bigr\},
````
where $[p]$ is an equivalence class of points $p$, and $∼$ indicates equivalence.
For example, the real projective space $ℝℙ^n$ is represented as the unit sphere $𝕊^n$, where
antipodal points are considered equivalent.

The tangent space at point $p$ is given by

````math
T_p 𝔽ℙ^{n} := \bigl\{ X ∈ 𝔽^{n+1}\ \big|\ ⟨p,X⟩ = 0 \bigr \},
````
where $⟨⋅,⋅⟩$ denotes the inner product in the embedding $𝔽^{n+1}$.

When $𝔽 = ℍ$, this implementation of $ℍℙ^n$ is the right-quaternionic projective
space.

# Constructor

    ProjectiveSpace(n[, field=ℝ])

Generate the projective space $𝔽ℙ^{n} ⊂ 𝔽^{n+1}$, defaulting to the real projective space
$ℝℙ^n$, where `field` can also be used to generate the complex- and right-quaternionic
projective spaces.
"""
struct ProjectiveSpace{N,𝔽} <: AbstractProjectiveSpace{𝔽} end
ProjectiveSpace(n::Int, field::AbstractNumbers=ℝ) = ProjectiveSpace{n,field}()

function active_traits(f, ::AbstractProjectiveSpace, args...)
    return merge_traits(IsIsometricEmbeddedManifold())
end

@doc raw"""
    ArrayProjectiveSpace{T<:Tuple,𝔽} <: AbstractProjectiveSpace{𝔽}

The projective space $𝔽ℙ^{n₁,n₂,…,nᵢ}$ is the manifold of all lines in $𝔽^{n₁,n₂,…,nᵢ}$.
The default representation is in the embedding, i.e. as unit (Frobenius) norm matrices in
$𝔽^{n₁,n₂,…,nᵢ}$:

````math
𝔽ℙ^{n_1, n_2, …, n_i} := \bigl\{ [p] ⊂ 𝔽^{n_1, n_2, …, n_i} \ \big|\ \lVert p \rVert_{\mathrm{F}} = 1, λ ∈ 𝔽, |λ| = 1, p ∼ p λ \bigr\}.
````
where $[p]$ is an equivalence class of points $p$, $\sim$ indicates equivalence, and
$\lVert ⋅ \rVert_{\mathrm{F}}$ is the Frobenius norm.
Note that unlike [`ProjectiveSpace`](@ref), the argument for `ArrayProjectiveSpace`
is given by the size of the embedding.
This means that [`ProjectiveSpace(2)`](@ref) and `ArrayProjectiveSpace(3)` are the same
manifold.
Additionally, `ArrayProjectiveSpace(n,1;field=𝔽)` and [`Grassmann(n,1;field=𝔽)`](@ref) are
the same.

The tangent space at point $p$ is given by

````math
T_p 𝔽ℙ^{n_1, n_2, …, n_i} := \bigl\{ X ∈ 𝔽^{n_1, n_2, …, n_i}\ |\ ⟨p,X⟩_{\mathrm{F}} = 0 \bigr \},
````

where $⟨⋅,⋅⟩_{\mathrm{F}}$ denotes the (Frobenius) inner product in the embedding
$𝔽^{n_1, n_2, …, n_i}$.

# Constructor

    ArrayProjectiveSpace(n₁,n₂,...,nᵢ; field=ℝ)

Generate the projective space $𝔽ℙ^{n_1, n_2, …, n_i}$, defaulting to the real projective
space, where `field` can also be used to generate the complex- and right-quaternionic
projective spaces.
"""
struct ArrayProjectiveSpace{N,𝔽} <: AbstractProjectiveSpace{𝔽} where {N<:Tuple} end
function ArrayProjectiveSpace(n::Vararg{Int,I}; field::AbstractNumbers=ℝ) where {I}
    return ArrayProjectiveSpace{Tuple{n...},field}()
end

function allocation_promotion_function(::AbstractProjectiveSpace{ℂ}, f, args::Tuple)
    return complex
end

@doc raw"""
    check_point(M::AbstractProjectiveSpace, p; kwargs...)

Check whether `p` is a valid point on the [`AbstractProjectiveSpace`](@ref) `M`, i.e.
that it has the same size as elements of the embedding and has unit Frobenius norm.
The tolerance for the norm check can be set using the `kwargs...`.
"""
function check_point(M::AbstractProjectiveSpace, p; kwargs...)
    if !isapprox(norm(p), 1; kwargs...)
        return DomainError(
            norm(p),
            "The point $(p) does not lie on the $(M) since its norm is not 1.",
        )
    end
    return nothing
end

@doc raw"""
    check_vector(M::AbstractProjectiveSpace, p, X; kwargs... )

Check whether `X` is a tangent vector in the tangent space of `p` on the
[`AbstractProjectiveSpace`](@ref) `M`, i.e. that `X` has the same size as elements of the
tangent space of the embedding and that the Frobenius inner product
$⟨p, X⟩_{\mathrm{F}} = 0$.
"""
function check_vector(M::AbstractProjectiveSpace, p, X; kwargs...)
    if !isapprox(dot(p, X), 0; kwargs...)
        return DomainError(
            dot(p, X),
            "The vector $(X) is not a tangent vector to $(p) on $(M), since it is not" *
            " orthogonal in the embedding.",
        )
    end
    return nothing
end

function decorated_manifold(M::AbstractProjectiveSpace{𝔽}) where {𝔽}
    return Euclidean(representation_size(M)...; field=𝔽)
end

get_embedding(M::AbstractProjectiveSpace) = decorated_manifold(M)

embed(::AbstractProjectiveSpace, p) = p
embed(::AbstractProjectiveSpace, p, X) = X

@doc raw"""
    distance(M::AbstractProjectiveSpace, p, q)

Compute the Riemannian distance on [`AbstractProjectiveSpace`](@ref) `M`$=𝔽ℙ^n$ between
points `p` and `q`, i.e.
````math
d_{𝔽ℙ^n}(p, q) = \arccos\bigl| ⟨p, q⟩_{\mathrm{F}} \bigr|.
````

Note that this definition is similar to that of the [`AbstractSphere`](@ref).
However, the absolute value ensures that all equivalent `p` and `q` have the same pairwise
distance.
"""
distance(::AbstractProjectiveSpace, p, q) = acos(min(abs(dot(p, q)), 1))

function exp!(M::AbstractProjectiveSpace, q, p, X)
    θ = norm(M, p, X)
    q .= cos(θ) .* p .+ usinc(θ) .* X
    return q
end

function get_basis(::ProjectiveSpace{n,ℝ}, p, B::DiagonalizingOrthonormalBasis{ℝ}) where {n}
    return get_basis(Sphere{n,ℝ}(), p, B)
end

@doc raw"""
    get_coordinates(M::AbstractProjectiveSpace, p, X, B::DefaultOrthonormalBasis{ℝ})

Represent the tangent vector $X$ at point $p$ from the [`AbstractProjectiveSpace`](@ref)
$M = 𝔽ℙ^n$ in an orthonormal basis by unitarily transforming the hyperplane containing $X$,
whose normal is $p$, to the hyperplane whose normal is the $x$-axis.

Given $q = p \overline{λ} + x$, where
$λ = \frac{⟨x, p⟩_{\mathrm{F}}}{|⟨x, p⟩_{\mathrm{F}}|}$, $⟨⋅, ⋅⟩_{\mathrm{F}}$ denotes the
Frobenius inner product, and $\overline{⋅}$ denotes complex or quaternionic conjugation, the
formula for $Y$ is
````math
\begin{pmatrix}0 \\ Y\end{pmatrix} = \left(X - q\frac{2 ⟨q, X⟩_{\mathrm{F}}}{⟨q, q⟩_{\mathrm{F}}}\right)\overline{λ}.
````
"""
get_coordinates(::AbstractProjectiveSpace{ℝ}, p, X, ::DefaultOrthonormalBasis)

function get_coordinates_orthonormal!(
    M::AbstractProjectiveSpace{𝔽},
    Y,
    p,
    X,
    ::RealNumbers,
) where {𝔽}
    n = div(manifold_dimension(M), real_dimension(𝔽))
    z = p[1]
    cosθ = abs(z)
    λ = nzsign(z, cosθ)
    pend, Xend = view(p, 2:(n + 1)), view(X, 2:(n + 1))
    factor = λ' * X[1] / (1 + cosθ)
    Y .= (Xend .- pend .* factor) .* λ'
    return Y
end

@doc raw"""
    get_vector(M::AbstractProjectiveSpace, p, X, B::DefaultOrthonormalBasis{ℝ})

Convert a one-dimensional vector of coefficients $X$ in the basis `B` of the tangent space
at $p$ on the [`AbstractProjectiveSpace`](@ref) $M=𝔽ℙ^n$ to a tangent vector $Y$ at $p$ by
unitarily transforming the hyperplane containing $X$, whose normal is the $x$-axis, to the
hyperplane whose normal is $p$.

Given $q = p \overline{λ} + x$, where
$λ = \frac{⟨x, p⟩_{\mathrm{F}}}{|⟨x, p⟩_{\mathrm{F}}|}$, $⟨⋅, ⋅⟩_{\mathrm{F}}$ denotes the
Frobenius inner product, and $\overline{⋅}$ denotes complex or quaternionic conjugation, the
formula for $Y$ is
````math
Y = \left(X - q\frac{2 \left\langle q, \begin{pmatrix}0 \\ X\end{pmatrix}\right\rangle_{\mathrm{F}}}{⟨q, q⟩_{\mathrm{F}}}\right) λ.
````
"""
get_vector(::AbstractProjectiveSpace, p, X, ::DefaultOrthonormalBasis{ℝ})

function get_vector_orthonormal!(
    M::AbstractProjectiveSpace{𝔽},
    Y,
    p,
    X,
    ::RealNumbers,
) where {𝔽}
    n = div(manifold_dimension(M), real_dimension(𝔽))
    z = p[1]
    cosθ = abs(z)
    λ = nzsign(z, cosθ)
    pend = view(p, 2:(n + 1))
    pX = dot(pend, X)
    Y[1] = -λ * pX * λ
    Y[2:(n + 1)] .= (X .- pend .* (pX / (1 + cosθ))) .* λ
    return Y
end

injectivity_radius(::AbstractProjectiveSpace) = π / 2
injectivity_radius(::AbstractProjectiveSpace, p) = π / 2
injectivity_radius(::AbstractProjectiveSpace, ::AbstractRetractionMethod) = π / 2
injectivity_radius(::AbstractProjectiveSpace, p, ::AbstractRetractionMethod) = π / 2

@doc raw"""
    inverse_retract(M::AbstractProjectiveSpace, p, q, method::ProjectionInverseRetraction)
    inverse_retract(M::AbstractProjectiveSpace, p, q, method::PolarInverseRetraction)
    inverse_retract(M::AbstractProjectiveSpace, p, q, method::QRInverseRetraction)

Compute the equivalent inverse retraction [`ProjectionInverseRetraction`](https://juliamanifolds.github.io/ManifoldsBase.jl/stable/retractions.html#ManifoldsBase.ProjectionInverseRetraction),
[`PolarInverseRetraction`](https://juliamanifolds.github.io/ManifoldsBase.jl/stable/retractions.html#ManifoldsBase.PolarInverseRetraction), and [`QRInverseRetraction`](https://juliamanifolds.github.io/ManifoldsBase.jl/stable/retractions.html#ManifoldsBase.QRInverseRetraction) on the
[`AbstractProjectiveSpace`](@ref) manifold `M`$=𝔽ℙ^n$, i.e.
````math
\operatorname{retr}_p^{-1} q = q \frac{1}{⟨p, q⟩_{\mathrm{F}}} - p,
````
where $⟨⋅, ⋅⟩_{\mathrm{F}}$ is the Frobenius inner product.

Note that this inverse retraction is equivalent to the three corresponding inverse
retractions on [`Grassmann(n+1,1,𝔽)`](@ref), where the three inverse retractions in this
case coincide.
For $ℝℙ^n$, it is the same as the `ProjectionInverseRetraction` on the real
[`Sphere`](@ref).
"""
inverse_retract(
    ::AbstractProjectiveSpace,
    p,
    q,
    ::Union{ProjectionInverseRetraction,PolarInverseRetraction,QRInverseRetraction},
)

function inverse_retract_qr!(::AbstractProjectiveSpace, X, p, q)
    X .= q ./ dot(p, q) .- p
    return X
end
function inverse_retract_polar!(::AbstractProjectiveSpace, X, p, q)
    X .= q ./ dot(p, q) .- p
    return X
end
function inverse_retract_project!(::AbstractProjectiveSpace, X, p, q)
    X .= q ./ dot(p, q) .- p
    return X
end

@doc raw"""
    isapprox(M::AbstractProjectiveSpace, p, q; kwargs...)

Check that points `p` and `q` on the [`AbstractProjectiveSpace`](@ref) `M`$=𝔽ℙ^n$ are
members of the same equivalence class, i.e. that $p = q λ$ for some element $λ ∈ 𝔽$ with
unit absolute value, that is, $|λ| = 1$.
This is equivalent to the Riemannian
[`distance`](@ref distance(::AbstractProjectiveSpace, p, q)) being 0.
"""
function Base.isapprox(::AbstractProjectiveSpace, p, q; kwargs...)
    return isapprox(abs(dot(p, q)), 1; kwargs...)
end

@doc raw"""
    log(M::AbstractProjectiveSpace, p, q)

Compute the logarithmic map on [`AbstractProjectiveSpace`](@ref) `M`$ = 𝔽ℙ^n$,
i.e. the tangent vector whose corresponding [`geodesic`](https://juliamanifolds.github.io/ManifoldsBase.jl/stable/functions.html#ManifoldsBase.geodesic-Tuple{AbstractManifold,%20Any,%20Any}) starting from `p`
reaches `q` after time 1 on `M`. The formula reads

````math
\log_p q = (q λ - \cos θ p) \frac{θ}{\sin θ},
````
where $θ = \arccos|⟨q, p⟩_{\mathrm{F}}|$ is the
[`distance`](@ref distance(::AbstractProjectiveSpace, p, q)) between $p$ and $q$,
$⟨⋅, ⋅⟩_{\mathrm{F}}$ is the Frobenius inner product, and
$λ = \frac{⟨q, p⟩_{\mathrm{F}}}{|⟨q, p⟩_{\mathrm{F}}|} ∈ 𝔽$ is the unit scalar that
minimizes $d_{𝔽^{n+1}}(p - q λ)$.
That is, $q λ$ is the member of the equivalence class $[q]$ that is closest to $p$ in the
embedding.
As a result, $\exp_p \circ \log_p \colon q ↦ q λ$.

The logarithmic maps for the real [`AbstractSphere`](@ref) $𝕊^n$ and the real projective
space $ℝℙ^n$ are identical when $p$ and $q$ are in the same hemisphere.
"""
log(::AbstractProjectiveSpace, p, q)

function log!(M::AbstractProjectiveSpace, X, p, q)
    z = dot(q, p)
    cosθ = abs(z)
    λ = nzsign(z, cosθ)
    X .= (q .* λ .- cosθ .* p) ./ usinc_from_cos(cosθ)
    return project!(M, X, p, X)
end

@doc raw"""
    manifold_dimension(M::AbstractProjectiveSpace{𝔽}) where {𝔽}

Return the real dimension of the [`AbstractProjectiveSpace`](@ref) `M`, respectively i.e.
the real dimension of the embedding minus the real dimension of the field `𝔽`.
"""
function manifold_dimension(M::AbstractProjectiveSpace{𝔽}) where {𝔽}
    return manifold_dimension(get_embedding(M)) - real_dimension(𝔽)
end

"""
    mean(
        M::AbstractProjectiveSpace,
        x::AbstractVector,
        [w::AbstractWeights,]
        method = GeodesicInterpolationWithinRadius(π/4);
        kwargs...,
    )

Compute the Riemannian [`mean`](@ref mean(M::AbstractManifold, args...)) of points in vector `x`
using [`GeodesicInterpolationWithinRadius`](@ref).
"""
mean(::AbstractProjectiveSpace, ::Any...)

function default_estimation_method(::AbstractProjectiveSpace, ::typeof(mean))
    return GeodesicInterpolationWithinRadius(π / 4)
end

function mid_point!(M::ProjectiveSpace, q, p1, p2)
    z = dot(p2, p1)
    λ = nzsign(z)
    q .= p1 .+ p2 .* λ
    project!(M, q, q)
    return q
end

@doc raw"""
    project(M::AbstractProjectiveSpace, p)

Orthogonally project the point `p` from the embedding onto the
[`AbstractProjectiveSpace`](@ref) `M`:
````math
\operatorname{proj}(p) = \frac{p}{\lVert p \rVert}_{\mathrm{F}},
````
where $\lVert ⋅ \rVert_{\mathrm{F}}$ denotes the Frobenius norm.
This is identical to projection onto the [`AbstractSphere`](@ref).
"""
project(::AbstractProjectiveSpace, ::Any)

project!(::AbstractProjectiveSpace, q, p) = (q .= p ./ norm(p))

@doc raw"""
    project(M::AbstractProjectiveSpace, p, X)

Orthogonally project the point `X` onto the tangent space at `p` on the
[`AbstractProjectiveSpace`](@ref) `M`:

````math
\operatorname{proj}_p (X) = X - p⟨p, X⟩_{\mathrm{F}},
````
where $⟨⋅, ⋅⟩_{\mathrm{F}}$ denotes the Frobenius inner product.
For the real [`AbstractSphere`](@ref) and `AbstractProjectiveSpace`, this projection is the
same.
"""
project(::AbstractProjectiveSpace, ::Any, ::Any)

project!(::AbstractProjectiveSpace, Y, p, X) = (Y .= X .- p .* dot(p, X))

@doc raw"""
    representation_size(M::AbstractProjectiveSpace)

Return the size points on the [`AbstractProjectiveSpace`](@ref) `M` are represented as,
i.e., the representation size of the embedding.
"""
@generated representation_size(::ArrayProjectiveSpace{N}) where {N} = size_to_tuple(N)
@generated representation_size(::ProjectiveSpace{N}) where {N} = (N + 1,)

@doc raw"""
    retract(M::AbstractProjectiveSpace, p, X, method::ProjectionRetraction)
    retract(M::AbstractProjectiveSpace, p, X, method::PolarRetraction)
    retract(M::AbstractProjectiveSpace, p, X, method::QRRetraction)

Compute the equivalent retraction [`ProjectionRetraction`](https://juliamanifolds.github.io/ManifoldsBase.jl/stable/retractions.html#ManifoldsBase.ProjectionRetraction), [`PolarRetraction`](https://juliamanifolds.github.io/ManifoldsBase.jl/stable/retractions.html#ManifoldsBase.PolarRetraction),
and [`QRRetraction`](https://juliamanifolds.github.io/ManifoldsBase.jl/stable/retractions.html#ManifoldsBase.QRRetraction) on the [`AbstractProjectiveSpace`](@ref) manifold `M`$=𝔽ℙ^n$,
i.e.
````math
\operatorname{retr}_p X = \operatorname{proj}_p(p + X).
````

Note that this retraction is equivalent to the three corresponding retractions on
[`Grassmann(n+1,1,𝔽)`](@ref), where in this case they coincide.
For $ℝℙ^n$, it is the same as the `ProjectionRetraction` on the real [`Sphere`](@ref).
"""
retract(
    ::AbstractProjectiveSpace,
    p,
    X,
    ::Union{ProjectionRetraction,PolarRetraction,QRRetraction},
)

function retract_polar!(M::AbstractProjectiveSpace, q, p, X)
    q .= p .+ X
    return project!(M, q, q)
end
function retract_project!(M::AbstractProjectiveSpace, q, p, X)
    q .= p .+ X
    return project!(M, q, q)
end
function retract_qr!(M::AbstractProjectiveSpace, q, p, X)
    q .= p .+ X
    return project!(M, q, q)
end

function Base.show(io::IO, ::ProjectiveSpace{n,𝔽}) where {n,𝔽}
    return print(io, "ProjectiveSpace($(n), $(𝔽))")
end
function Base.show(io::IO, ::ArrayProjectiveSpace{N,𝔽}) where {N,𝔽}
    return print(io, "ArrayProjectiveSpace($(join(N.parameters, ", ")); field = $(𝔽))")
end

"""
    uniform_distribution(M::ProjectiveSpace{n,ℝ}, p) where {n}

Uniform distribution on given [`ProjectiveSpace`](@ref) `M`. Generated points will be of
similar type as `p`.
"""
function uniform_distribution(M::ProjectiveSpace{n,ℝ}, p) where {n}
    d = Distributions.MvNormal(zero(p), 1.0)
    return ProjectedPointDistribution(M, d, project!, p)
end

@doc raw"""
    parallel_transport_to(M::AbstractProjectiveSpace, p, X, q)

Parallel transport a vector `X` from the tangent space at a point `p` on the
[`AbstractProjectiveSpace`](@ref) `M`$=𝔽ℙ^n$ to the tangent space at another point `q`.

This implementation proceeds by transporting $X$ to $T_{q λ} M$ using the same approach as
[`parallel_transport_direction`](@ref parallel_transport_direction(::AbstractProjectiveSpace, p, X, d)),
where $λ = \frac{⟨q, p⟩_{\mathrm{F}}}{|⟨q, p⟩_{\mathrm{F}}|} ∈ 𝔽$ is the unit scalar that
takes $q$ to the member $q λ$ of its equivalence class $[q]$ closest to $p$ in the
embedding.
It then maps the transported vector from $T_{q λ} M$ to $T_{q} M$.
The resulting transport to $T_{q} M$ is
````math
\mathcal{P}_{q ← p}(X) = \left(X - \left(p \frac{\sin θ}{θ} + d \frac{1 - \cos θ}{θ^2}\right) ⟨d, X⟩_p\right) \overline{λ},
````
where $d = \log_p q$ is the direction of the transport, $θ = \lVert d \rVert_p$ is the
[`distance`](@ref distance(::AbstractProjectiveSpace, p, q)) between $p$ and $q$, and
$\overline{⋅}$ denotes complex or quaternionic conjugation.
"""
parallel_transport_to(::AbstractProjectiveSpace, ::Any, ::Any, ::Any)

function parallel_transport_to!(::AbstractProjectiveSpace, Y, p, X, q)
    z = dot(q, p)
    λ = nzsign(z)
    m = p .+ q .* λ # un-normalized midpoint
    mnorm2 = real(dot(m, m))
    factor = λ' * dot(q, X) * (2 / mnorm2) # λ' * dot(q, X) ≡ dot(q * λ, X)
    # multiply by λ' to bring from T_{\exp_p(\log_p q)} M to T_q M
    # this ensures that subsequent functions like `exp(M, q, Y)` do the right thing
    Y .= (X .- m .* factor) .* λ'
    return Y
end
function vector_transport_to_project!(M::AbstractProjectiveSpace, Y, p, X, q)
    project!(M, Y, q, X)
    return Y
end

@doc raw"""
    parallel_transport_direction(M::AbstractProjectiveSpace, p, X, d)

Parallel transport a vector `X` from the tangent space at a point `p` on the
[`AbstractProjectiveSpace`](@ref) `M` along the [`geodesic`](https://juliamanifolds.github.io/ManifoldsBase.jl/stable/functions.html#ManifoldsBase.geodesic-Tuple{AbstractManifold,%20Any,%20Any}) in the direction
indicated by the tangent vector `d`, i.e.
````math
\mathcal{P}_{\exp_p (d) ← p}(X) = X - \left(p \frac{\sin θ}{θ} + d \frac{1 - \cos θ}{θ^2}\right) ⟨d, X⟩_p,
````
where $θ = \lVert d \rVert$, and $⟨⋅, ⋅⟩_p$ is the [`inner`](@ref) product at the point $p$.
For the real projective space, this is equivalent to the same vector transport on the real
[`AbstractSphere`](@ref).
"""
parallel_transport_direction(::AbstractProjectiveSpace, ::Any, ::Any, ::Any)

function parallel_transport_direction!(M::AbstractProjectiveSpace, Y, p, X, d)
    θ = norm(M, p, d)
    cosθ = cos(θ)
    dX = inner(M, p, d, X)
    α = usinc(θ) * dX
    β = ifelse(iszero(θ), zero(cosθ), (1 - cosθ) / θ^2) * dX
    Y .= X .- p .* α .- d .* β
    return Y
end


text:
# Common utilities

##### common types

"""
Base type for the output of clustering algorithm.
"""
abstract type ClusteringResult end

# generic functions

"""
    nclusters(R::ClusteringResult)

Get the number of clusters.
"""
nclusters(R::ClusteringResult) = length(R.counts)

"""
    counts(R::ClusteringResult)

Get the vector of cluster sizes.

`counts(R)[k]` is the number of points assigned to the ``k``-th cluster.
"""
counts(R::ClusteringResult) = R.counts

"""
    assignments(R::ClusteringResult)

Get the vector of cluster indices for each point.

`assignments(R)[i]` is the index of the cluster to which the ``i``-th point
is assigned.
"""
assignments(R::ClusteringResult) = R.assignments

##### convert display symbol to disp level

display_level(s::Symbol) =
    s == :none ? 0 :
    s == :final ? 1 :
    s == :iter ? 2 :
    error("Invalid value for the option 'display'.")

##### update minimum value

function updatemin!(r::AbstractArray, x::AbstractArray)
    n = length(r)
    length(x) == n || throw(DimensionMismatch("Inconsistent array lengths."))
    @inbounds for i = 1:n
        xi = x[i]
        if xi < r[i]
            r[i] = xi
        end
    end
    return r
end


text:
module QuantumNumbers

using Base.Iterators: Reverse,flatten,product,reverse
using Printf: @printf,@sprintf
using DataStructures: OrderedDict
using Random: MersenneTwister,seed!,shuffle!
using LinearAlgebra: norm
using ..Combinatorics: Combinations
using ...Prerequisites: Float
using ...Prerequisites.NamedVectors: HomoNamedVector

import ...Interfaces: ⊕,⊗,dimension,expand,permute,decompose,regularize!,regularize

export AbelianNumber,@abeliannumber,SQN,PQN,SPQN,Z2QN,periods
export qnscounts,qnsindptr,qnscompression,qnsexpansion,qnscontents,qnsindices,qnsbruteforce,qnsmontecarlo
export AbelianNumbers,SQNS,PQNS,SzPQNS,SPQNS,Z2QNS
export ukron,toordereddict

"Abstract type for all concrete quantum numbers for a single basis."
abstract type AbelianNumber <: HomoNamedVector{Float} end

"""
    dimension(::Type{<:AbelianNumber}) -> Int
    dimension(::AbelianNumber) -> Int

The dimension of the Hilbert space an `AbelianNumber` represents. Apparently, this is always 1.
"""
dimension(::Type{<:AbelianNumber})=1
dimension(::AbelianNumber)=1

"""
    regularize!(::Type{QN},array::AbstractVector{<:Real}) where QN<:AbelianNumber -> typeof(array)
    regularize!(::Type{QN},array::AbstractMatrix{<:Real}) where QN<:AbelianNumber -> typeof(array)

Regularize the elements of an array in place so that it can represent quantum numbers.
"""
function regularize!(::Type{QN},array::AbstractVector{<:Real}) where QN<:AbelianNumber
    @assert array|>length==QN|>length "regularize! error: not consistent shape of input array and $QN."
    for (i,period) in enumerate(QN|>periods)
        period==Inf || @inbounds begin
            remainder=array[i]%period
            array[i]=remainder<0 ? remainder+period : remainder
        end
    end
    return array
end
function regularize!(::Type{QN},array::AbstractMatrix{<:Real}) where QN<:AbelianNumber
    @assert size(array,1)==QN|>length "regularize! error: not consistent shape of input array and $QN."
    for (i,period) in enumerate(QN|>periods)
        period==Inf || @inbounds for j in 1:size(array,2)
            remainder=array[i,j]%period
            array[i,j]=remainder<0 ? remainder+period : remainder
        end
    end
    return array
end

"""
    regularize(::Type{QN},array::Union{AbstractVector{<:Real},AbstractMatrix{<:Real}}) where {QN<:AbelianNumber} -> typeof(array)

Regularize the elements of an array and return a copy that can represent quantum numbers.
"""
function regularize(::Type{QN},array::Union{AbstractVector{<:Real},AbstractMatrix{<:Real}}) where {QN<:AbelianNumber}
    result=copy(array)
    regularize!(QN,result)
    return result
end

"""
    @abeliannumber typename fieldnames fieldperiods

Construct a concrete `AbelianNumber` with the type name being `typename`, fieldnames specified by `fieldnames` and periods specified by `fieldperiods`.
"""
macro abeliannumber(typename,fieldnames,fieldperiods)
    typename=Symbol(typename)
    fieldnames=tuple(eval(fieldnames)...)
    fieldperiods=tuple(eval(fieldperiods)...)
    arguments=ntuple(i->Symbol(:v,i),length(fieldnames))
    @assert length(fieldnames)==length(fieldperiods) "abeliannumber error: number of fieldnames($(length(fieldnames))) and fieldperiods($(length(fieldperiods))) not equal."
    @assert all(isa(name,Symbol) for name in fieldnames) "abeliannumber error: all field names should be Symbol."
    @assert all(fieldperiods.>0) "abeliannumber error: all field periods should be positive."
    if all(fieldperiods.==Inf)
        title=Expr(:call,:($(esc(typename))),(:($arg::Float) for arg in arguments)...)
        body=Expr(:call,:new,arguments...)
    else
        title=Expr(:call,:($(esc(typename))),(:($arg::Float) for arg in arguments)...,Expr(:kw,:(regularize::Bool),:true))
        body=Expr(:call,:new,(p==Inf ? arg : :(regularize ? ($arg)%($p)+(($arg)%($p)<0 ? $p : 0) : $arg) for (p,arg) in zip(fieldperiods,arguments))...)
    end
    newtype=Expr(:struct,false,:($(esc(typename))<:AbelianNumber),Expr(:block,(:($field::Float) for field in fieldnames)...,Expr(:(=),title,body)))
    functions=Expr(:block,:(Base.fieldnames(::Type{<:$(esc(typename))})=$fieldnames),:(periods(::Type{<:$(esc(typename))})=$fieldperiods))
    return Expr(:block,:(global periods),:(Base.@__doc__($newtype)),functions)
end

"""
    SQN(Sz::Real)

The concrete `AbelianNumber` of a quantum system with spin z-component `Sz` conserved.
"""
@abeliannumber "SQN" (:Sz,) (Inf,)

"""
    PQN(N::Real)

The concrete `AbelianNumber` of a quantum system with particle number `N` conserved.
"""
@abeliannumber "PQN" (:N,) (Inf,)

"""
    SPQN(N::Real,Sz::Real)

The concrete `AbelianNumber` of a quantum system with both particle number `N` and spin z-component `Sz` conserved.
"""
@abeliannumber "SPQN" (:N,:Sz) (Inf,Inf)

"""
    Z2QN(N::Real)

The concrete `AbelianNumber` of a quantum system with a Z₂-like conserved quantity.
"""
@abeliannumber "Z2QN" (:N,) (2,)

abstract type QnsProtocol end
struct QnsIndptr <: QnsProtocol end
struct QnsCounts <: QnsProtocol end
"""
    qnsindptr

Indicate that methods with `AbelianNumbers` use the index pointer of the compressed contents.
"""
const qnsindptr=QnsIndptr()
"""
    qnscounts

Indicate that methods with `AbelianNumbers` use the count number of the compressed contents.
"""
const qnscounts=QnsCounts()

"""
    AbelianNumbers(form::Char,contents::Vector{<:AbelianNumber},counts::Vector{Int},::QnsCounts)
    AbelianNumbers(form::Char,contents::Vector{<:AbelianNumber},indptr::Vector{Int},::QnsIndptr)

The whole quantum numbers of the total bases of a Hilbert space.

The default constructors construct an `AbelianNumbers` from a vector of concrete quantum numbers and an vector containing their counts or indptr.
"""
struct AbelianNumbers{QN<:AbelianNumber}
    form::Char
    contents::Vector{QN}
    indptr::Vector{Int}
    function AbelianNumbers(form::Char,contents::Vector{<:AbelianNumber},counts::Vector{Int},::QnsCounts)
        @assert form|>uppercase ∈ ('G','U','C') "AbelianNumbers error: 'form'($form) is not 'G','U' or 'C'."
        @assert length(contents)==length(counts) "AbelianNumbers error: dismatch lengths of contents and counts ($(length(contents))!=$length(counts))."
        return new{contents|>eltype}(form|>uppercase,contents,[0,cumsum(counts)...])
    end
    function AbelianNumbers(form::Char,contents::Vector{<:AbelianNumber},indptr::Vector{Int},::QnsIndptr)
        @assert form|>uppercase ∈ ('G','U','C') "AbelianNumbers error: 'form'($form) is not 'G','U' or 'C'."
        @assert length(contents)+1==length(indptr) "AbelianNumbers error: dismatch shapes of contents and indptr ($(length(contents))+1!=$length(indptr))."
        return new{contents|>eltype}(form|>uppercase,contents,indptr)
    end
end

"""
    AbelianNumbers(qn::AbelianNumber,count::Int=1)

Construct an `AbelianNumbers` with one unique quantum number which occurs `count` times.
"""
AbelianNumbers(qn::AbelianNumber,count::Int=1)=AbelianNumbers('C',[qn],[0,count],qnsindptr)

"""
    AbelianNumbers(od::OrderedDict{<:AbelianNumber,Int})

Construct an `AbelianNumbers` from an ordered dict containing concrete quantum numbers and their counts.
"""
function AbelianNumbers(od::OrderedDict{<:AbelianNumber,Int})
    contents=Vector{od|>keytype}(undef,length(od))
    indptr=zeros(Int,length(od)+1)
    for (i,(qn,count)) in enumerate(od)
        @inbounds contents[i]=qn
        @inbounds indptr[i+1]=indptr[i]+count
    end
    AbelianNumbers('U',contents,indptr,qnsindptr)
end

"""
    AbelianNumbers(od::OrderedDict{<:AbelianNumber,UnitRange{Int}})

Construct an `AbelianNumbers` from an ordered dict containing concrete quantum numbers and their slices.
"""
function AbelianNumbers(od::OrderedDict{<:AbelianNumber,UnitRange{Int}})
    contents=Vector{od|>keytype}(undef,length(od))
    indptr=zeros(Int,length(od)+1)
    for (i,(qn,slice)) in enumerate(od)
        @inbounds contents[i]=qn
        @inbounds @assert indptr[i]+1==slice.start "AbelianNumbers error: slice not consistent."
        @inbounds indptr[i+1]=slice.stop
    end
    AbelianNumbers('U',contents,indptr,qnsindptr)
end

"""
    ==(qns1::AbelianNumbers,qns2::AbelianNumbers) -> Bool
    isequal(qns1::AbelianNumbers,qns2::AbelianNumbers) -> Bool

Overloaded equivalent operator. Two `AbelianNumbers`es are equal to each other if and only if both their `contents`es and `indptr`s are elementwise equal to each other.
!!! note
    It is not necessary for two `AbelianNumbers`es to have the same eltype nor the same form to be equal to each other.
"""
Base.:(==)(qns1::AbelianNumbers,qns2::AbelianNumbers)=qns1.contents==qns2.contents && qns1.indptr==qns2.indptr
Base.isequal(qns1::AbelianNumbers,qns2::AbelianNumbers)=isequal(qns1.contents,qns2.contents) && isequal(qns1.indptr,qns2.indptr)

"""
    show(io::IO,qns::AbelianNumbers)

Show an `AbelianNumbers`.
"""
Base.show(io::IO,qns::AbelianNumbers)=@printf io "QNS(%s)" join((@sprintf("%s=>%s",qn,slice) for (qn,slice) in pairs(qns,qnsindptr)),',')

"""
    string(qns::AbelianNumbers) -> String

Convert an `AbelianNumbers` to string.
"""
Base.string(qns::AbelianNumbers)=@sprintf "QNS(%s,%s)" qns|>length qns|>dimension

"""
    length(qns::AbelianNumbers) -> Int

Get the number of unduplicate qunatum numbers in the `AbelianNumbers`.
"""
Base.length(qns::AbelianNumbers)=length(qns.contents)

"""
    eltype(::Type{<:AbelianNumbers{QN}}) where {QN<:AbelianNumber}
    eltype(qns::AbelianNumbers)

Get the type of the concrete `AbelianNumber` contained in an `AbelianNumbers`.
"""
Base.eltype(::Type{<:AbelianNumbers{QN}}) where {QN<:AbelianNumber}=QN
Base.eltype(qns::AbelianNumbers)=qns|>typeof|>eltype

"""
    getindex(qns::AbelianNumbers,index::Int) -> eltype(qns)
    getindex(qns::AbelianNumbers,slice::UnitRange{Int}) -> AbelianNumbers
    getindex(qns::AbelianNumbers,indices::Vector{Int}) -> AbelianNumbers

Overloaded `[]` operator.
!!! note
    1. For an `AbelianNumbers`, all these `getindex` functions act on its `contents`, i.e. its compressed data, but not on its expansion, i.e. the uncompressed data. This definition is consistent with the [`length`](@ref) function.
    2. When the index is an integer, the result is an `AbelianNumber`, while when the index is a unit range or a vector of intgers, the result is an `AbelianNumbers`. The logic is quite reasonable because such behaviors are much alike to those of a vector container.
"""
Base.getindex(qns::AbelianNumbers,index::Int)=qns.contents[index]
function Base.getindex(qns::AbelianNumbers,slice::UnitRange{Int})
    contents=qns.contents[slice]
    indptr=qns.indptr[slice.start:slice.stop+1]
    indptr.=indptr.-qns.indptr[slice.start]
    return AbelianNumbers(qns.form,contents,indptr,qnsindptr)
end
function Base.getindex(qns::AbelianNumbers,indices::Vector{Int})
    contents=Vector{qns|>eltype}(undef,length(indices))
    indptr=zeros(Int,length(indices)+1)
    for (i,index) in enumerate(indices)
        contents[i]=qns.contents[index]
        indptr[i+1]=indptr[i]+qns.indptr[index+1]-qns.indptr[index]
    end
    return AbelianNumbers('U',contents,indptr,qnsindptr)
end

"""
    iterate(qns::AbelianNumbers,state::Int=1)
    iterate(rv::Iterators.Reverse{<:AbelianNumbers},state::Int=length(rv.itr,false))

Iterate or reversely iterate over the concrete `AbelianNumber`s contained in an `AbelianNumbers`.
"""
Base.iterate(qns::AbelianNumbers,state::Int=1)=state>length(qns) ? nothing : (@inbounds(qns.contents[state]),state+1)
Base.iterate(rv::Reverse{<:AbelianNumbers},state::Int=length(rv.itr))=state<1 ? nothing : (@inbounds(rv.itr.contents[state]),state-1)

"""
    keys(qns::AbelianNumbers) -> Vector{qns|>eltype}

Iterate over the concrete `AbelianNumber`s contained in an `AbelianNumbers`.
"""
Base.keys(qns::AbelianNumbers)=qns.contents

"""
    values(qns::AbelianNumbers,::QnsIndptr)
    values(qns::AbelianNumbers,::QnsCounts)

Iterate over the slices/counts of the `AbelianNumbers`.
"""
@views Base.values(qns::AbelianNumbers,::QnsIndptr)=((start+1):stop for (start,stop) in zip(qns.indptr[1:end-1],qns.indptr[2:end]))
@views Base.values(qns::AbelianNumbers,::QnsCounts)=(stop-start for (start,stop) in zip(qns.indptr[1:end-1],qns.indptr[2:end]))

"""
    pairs(qns::AbelianNumbers,choice::Union{QnsIndptr,QnsCounts})

Iterate over the `AbelianNumber=>slice` or `AbelianNumber=>count` pairs.
"""
Base.pairs(qns::AbelianNumbers,choice::Union{QnsIndptr,QnsCounts})=Base.Generator(=>,keys(qns),values(qns,choice))

"""
    toordereddict(qns::AbelianNumbers,::QnsIndptr) -> OrderedDict{qns|>eltype,UnitRange{Int}}
    toordereddict(qns::AbelianNumbers,::QnsCounts) -> OrderedDict{qns|>eltype,Int}

Convert an `AbelianNumbers` to an ordered dict.
"""
function toordereddict(qns::AbelianNumbers,::QnsIndptr)
    @assert qns.form != 'G' "toordereddict error: input `AbelianNumbers` cannot be `G` formed."
    result=OrderedDict{qns|>eltype,UnitRange{Int}}()
    for i=1:length(qns)
        @inbounds result[qns.contents[i]]=qns.indptr[i]+1:qns.indptr[i+1]
    end
    return result
end
function toordereddict(qns::AbelianNumbers,::QnsCounts)
    @assert qns.form != 'G' "toordereddict error: input `AbelianNumbers` cannot be `G` formed."
    result=OrderedDict{qns|>eltype,Int}()
    for i=1:length(qns)
        @inbounds result[qns.contents[i]]=qns.indptr[i+1]-qns.indptr[i]
    end
    return result
end

"""
    dimension(qns::AbelianNumbers) -> Int

The dimension of the Hilbert space an `AbelianNumbers` represents.
"""
dimension(qns::AbelianNumbers)=@inbounds(qns.indptr[end])

"""
    sort(qns::AbelianNumbers) -> AbelianNumbers,Vector{Int}

Sort the quantum numbers of an `AbelianNumber`, return the sorted `AbelianNumber` and the permutation array that sorts the expansion of the original `AbelianNumbers`.
"""
function Base.sort(qns::AbelianNumbers)
    ctpts=sortperm(qns.contents,alg=Base.Sort.QuickSort)
    masks=Vector{Bool}(undef,length(qns.contents));masks[1]=true
    unduplicate=1
    for i=2:length(ctpts)
        @inbounds masks[i]=qns.contents[ctpts[i]]≠qns.contents[ctpts[i-1]]
        @inbounds masks[i] && (unduplicate+=1)
    end
    contents=Vector{qns|>eltype}(undef,unduplicate)
    indptr=zeros(Int,unduplicate+1)
    permutation=Vector{Int}(undef,dimension(qns))
    qncount,ptcount=0,0
    for (mask,index) in zip(masks,ctpts)
        @inbounds mask && (qncount+=1;contents[qncount]=qns.contents[index];indptr[qncount+1]=indptr[qncount])
        @inbounds indptr[qncount+1]+=qns.indptr[index+1]-qns.indptr[index]
        for (i,p) in enumerate(@inbounds(qns.indptr[index]+1:qns.indptr[index+1]))
            @inbounds permutation[ptcount+i]=p
        end
        @inbounds ptcount+=qns.indptr[index+1]-qns.indptr[index]
    end
    return AbelianNumbers('C',contents,indptr,qnsindptr),permutation
end

struct QnsCompression   <: QnsProtocol end
struct QnsExpansion     <: QnsProtocol end
"""
    qnscompression

Indicate that [`findall`](@ref) and [`permute`](@ref) use the compressed contents.
"""
const qnscompression=QnsCompression()
"""
    qnsexpansion

Indicate that [`findall`](@ref) and [`permute`](@ref) use the expanded contents.
"""
const qnsexpansion=QnsExpansion()

"""
    findall(target::QN,qns::AbelianNumbers{QN},::QnsCompression) where QN<:AbelianNumber -> Vector{Int}
    findall(target::QN,qns::AbelianNumbers{QN},::QnsExpansion) where QN<:AbelianNumber -> Vector{Int}
    findall(targets::NTuple{N,QN},qns::AbelianNumbers{QN},::QnsCompression) where {N,QN<:AbelianNumber} -> Vector{Int}
    findall(targets::NTuple{N,QN},qns::AbelianNumbers{QN},::QnsExpansion) where {N,QN<:AbelianNumber} -> Vector{Int}

Find all the indices of the target quantum numbers in the contents ([`qnscompression`](@ref) case) or the expansion ([`qnsexpansion`](@ref) case) of an `AbelianNumbers`.
"""
Base.findall(target::QN,qns::AbelianNumbers{QN},::QnsCompression) where QN<:AbelianNumber=findall((target,),qns,QnsCompression())
Base.findall(target::QN,qns::AbelianNumbers{QN},::QnsExpansion) where QN<:AbelianNumber=findall((target,),qns,QnsExpansion())
function Base.findall(targets::NTuple{N,QN},qns::AbelianNumbers{QN},::QnsCompression) where {N,QN<:AbelianNumber}
    result=Int[]
    if qns.form=='C'
        for qn in targets
            range=searchsorted(qns.contents,qn)
            range.start<=range.stop && push!(result,range.start)
        end
    else
        for (i,qn) in enumerate(qns)
            qn ∈ targets && push!(result,i)
        end
    end
    result
end
function Base.findall(targets::NTuple{N,QN},qns::AbelianNumbers{QN},::QnsExpansion) where {N,QN<:AbelianNumber}
    result=Int[]
    for index in findall(targets,qns,qnscompression)
        @inbounds append!(result,(qns.indptr[index]+1):qns.indptr[index+1])
    end
    result
end

"""
    filter(target::QN,qns::AbelianNumbers{QN}) where QN<:AbelianNumber -> AbelianNumbers{QN}
    filter(targets::NTuple{N,QN},qns::AbelianNumbers{QN}) where {N,QN<:AbelianNumber} -> AbelianNumbers{QN}

Find a subset of an `AbelianNumbers` by picking out the quantum numbers in targets.
"""
Base.filter(target::QN,qns::AbelianNumbers{QN}) where QN<:AbelianNumber=qns[findall(target,qns,qnscompression)]
Base.filter(targets::NTuple{N,QN},qns::AbelianNumbers{QN}) where {N,QN<:AbelianNumber}=qns[findall(targets,qns,qnscompression)]

"""
    permute(qns::AbelianNumbers,permutation::Vector{Int},::QnsCompression) -> AbelianNumbers
    permute(qns::AbelianNumbers,permutation::Vector{Int},::QnsExpansion) -> AbelianNumbers

Reorder the quantum numbers contained in an `AbelianNumbers` with a permutation and return the new one.

For [`qnscompression`](@ref) case, the permutation is for the compressed contents of the original `AbelianNumbers` while for [`qnsexpansion`](@ref) case, the permutation is for the expanded contents of the original `AbelianNumbers`.
"""
permute(qns::AbelianNumbers,permutation::Vector{Int},::QnsCompression)=qns[permutation]
function permute(qns::AbelianNumbers,permutation::Vector{Int},::QnsExpansion)
    contents=Vector{qns|>eltype}(undef,length(permutation))
    indptr=zeros(Int,length(permutation)+1)
    expansion=expand(qns,qnsindices)
    for (i,p) in enumerate(permutation)
        contents[i]=qns.contents[expansion[p]]
        indptr[i+1]=indptr[i]+1
    end
    return AbelianNumbers('G',contents,indptr,qnsindptr)
end

struct QnsContents  <: QnsProtocol end
struct QnsIndices   <: QnsProtocol end
"""
    qnscontents

Indicate that [`expand`](@ref) uses the compressed/expanded contents.
"""
const qnscontents=QnsContents()
"""
    qnsindices

Indicate that [`expand`](@ref) uses the indices of the compressed/expanded contents.
"""
const qnsindices=QnsIndices()

"""
    expand(qns::AbelianNumbers,::QnsContents) -> Vector{qns|>eltype}
    expand(qns::AbelianNumbers,::QnsIndices) -> Vector{Int}

Expand the contents ([`qnscontents`](@ref) case) or indices ([`qnsindices`](@ref) case) of an `AbelianNumbers` to the uncompressed form.
"""
function expand(qns::AbelianNumbers,::QnsContents)
    result=Vector{qns|>eltype}(undef,dimension(qns))
    for i=1:length(qns)
        for j=@inbounds((qns.indptr[i]+1):qns.indptr[i+1])
            @inbounds result[j]=qns.contents[i]
        end
    end
    return result
end
function expand(qns::AbelianNumbers,::QnsIndices)
    result=Vector{Int}(undef,dimension(qns))
    for i=1:length(qns)
        for j=@inbounds((qns.indptr[i]+1):qns.indptr[i+1])
            @inbounds result[j]=i
        end
    end
    return result
end

"""
    +(qn::AbelianNumber) -> typeof(qn)
    +(qn::QN,qns::QN...) where QN<:AbelianNumber -> QN
    +(qns::AbelianNumbers) -> AbelianNumbers
    +(qn::QN,qns::AbelianNumbers{QN}) where QN<:AbelianNumber -> AbelianNumbers{QN}
    +(qns::AbelianNumbers{QN},qn::QN) where QN<:AbelianNumber -> AbelianNumbers{QN}

Overloaded `+` operator for `AbelianNumber` and `AbelianNumbers`.
!!! note
    1. The addition between an `AbelianNumbers` and an `AbelianNumber` is just a global shift of the contents of the `AbelianNumbers` by the `AbelianNumber`, therefore, the result is an `AbelianNumbers`.
    2. `+` cannot be used between two `AbelianNumbers` because the result is ambiguous. Instead, use `⊕` for direct sum and `⊗` for direct product.
    3. To ensure type stability, two `AbelianNumber` can be added together if and only if they are of the same type.
    4. Similarly, an `AbelianNumber` and an `AbelianNumbers` can be added together if and only if the former's type is the same with the latter's eltype.
"""
Base.:+(qn::AbelianNumber)=qn
Base.:+(qn::QN,qns::QN...) where QN<:AbelianNumber=map(+,qn,qns...)
Base.:+(qns::AbelianNumbers)=qns
Base.:+(qn::QN,qns::AbelianNumbers{QN}) where QN<:AbelianNumber=qns+qn
Base.:+(qns::AbelianNumbers{QN},qn::QN) where QN<:AbelianNumber=AbelianNumbers(qns.form=='G' ? 'G' : 'U',[iqn+qn for iqn in qns],qns.indptr,qnsindptr)

"""
    -(qn::AbelianNumber) -> typeof(qn)
    -(qn1::QN,qn2::QN) where QN<:AbelianNumber -> QN
    -(qns::AbelianNumbers) -> AbelianNumbers
    -(qn::QN,qns::AbelianNumbers{QN}) where QN<:AbelianNumber -> AbelianNumbers{QN}
    -(qns::AbelianNumbers{QN},qn::QN) where QN<:AbelianNumber -> AbelianNumbers{QN}

Overloaded `-` operator for `AbelianNumber` and `AbelianNumbers`.
!!! note
    1. The subtraction between an `AbelianNumbers` and an `AbelianNumber` is just a global shift of the contents of the `AbelianNumbers` by the `AbelianNumber`, therefore, the result is an `AbelianNumbers`.
    2. `-` cannot be used between two `AbelianNumbers` because the result is ambiguous. Instead, use `⊕` with signs for direct sum and `⊗` with signs for direct product.
    3. To ensure type stability, an `AbelianNumber` can be subtracted by another `AbelianNumber` if and only if they are of the same type.
    4. Similarly, an `AbelianNumber` can be subtracted by an `AbelianNumbers` or vice versa if and only if the former's type is the same with the latter's eltype.
"""
Base.:-(qn::AbelianNumber)=map(-,qn)
Base.:-(qn1::QN,qn2::QN) where QN<:AbelianNumber=map(-,qn1,qn2)
Base.:-(qns::AbelianNumbers)=AbelianNumbers(qns.form=='G' ? 'G' : 'U',-qns.contents,qns.indptr,qnsindptr)
Base.:-(qn::QN,qns::AbelianNumbers{QN}) where QN<:AbelianNumber=AbelianNumbers(qns.form=='G' ? 'G' : 'U',[qn-iqn for iqn in qns],qns.indptr,qnsindptr)
Base.:-(qns::AbelianNumbers{QN},qn::QN) where QN<:AbelianNumber=AbelianNumbers(qns.form=='G' ? 'G' : 'U',[iqn-qn for iqn in qns],qns.indptr,qnsindptr)

"""
    *(qn::AbelianNumber,factor::Integer) -> typeof(qn)
    *(factor::Integer,qn::AbelianNumber) -> typeof(qn)
    *(qns::AbelianNumbers,factor::Integer) -> AbelianNumbers
    *(factor::Integer,qns::AbelianNumbers) -> AbelianNumbers

Overloaded `*` operator for the multiplication between an integer and an `AbelianNumber` or an `AbelianNumbers`.
"""
@generated function Base.:*(qn::AbelianNumber,factor::Integer)
    exprs=Expr[:(getfield(qn,$i)*factor) for i=1:(qn|>fieldnames|>length)]
    return :(typeof(qn)($(exprs...)))
end
Base.:*(factor::Integer,qn::AbelianNumber)=qn*factor
Base.:*(qns::AbelianNumbers,factor::Integer)=AbelianNumbers(qns.form=='G' ? 'G' : 'U',[qn*factor for qn in qns.contents],qns.indptr,qnsindptr)
Base.:*(factor::Integer,qns::AbelianNumbers)=qns*factor

"""
    ^(qn::AbelianNumber,factor::Integer) -> typeof(qn)
    ^(qns::AbelianNumbers,factor::Integer) -> AbelianNumbers

Overloaded `^` operator for `AbelianNumber` and `AbelianNumbers`. This operation translates into the direct product of `factor` copies of `qn` or `qns`.
"""
Base.:^(qn::AbelianNumber,factor::Integer)=kron(NTuple{factor,qn|>typeof}(qn for i=1:factor)...)
Base.:^(qns::AbelianNumbers,factor::Integer)=kron(NTuple{factor,qns|>typeof}(qns for i=1:factor)...)

"""
    ⊕(qns::AbelianNumber...) -> AbelianNumbers{qns|>eltype}
    ⊕(qnses::AbelianNumbers...) -> qnses|>eltype

Get the direct sum of some `AbelianNumber`s or `AbelianNumbers`es.
"""
⊕(qns::AbelianNumber...)=union(qns...)
⊕(qnses::AbelianNumbers...)=union(qnses...)

"""
    ⊗(qns::AbelianNumber...) -> eltype(qns)
    ⊗(qnses::AbelianNumbers...) -> eltype(qnses)

Get the direct product of some `AbelianNumber`s or `AbelianNumbers`es.
"""
⊗(qns::AbelianNumber...)=kron(qns...)
⊗(qnses::AbelianNumbers...)=kron(qnses...)

"""
    union(qns::Vararg{<:AbelianNumber,N};signs::NTuple{N,Int}=ntuple(i->1,N)) where N -> AbelianNumbers
    union(qnses::Vararg{AbelianNumbers{QN},N};signs::NTuple{N,Int}=ntuple(i->1,N)) where {N,QN<:AbelianNumber} -> AbelianNumbers{QN}

Get the direct sum of some `AbelianNumber`s or `AbelianNumbers`es.
!!! note
    1. Physically, the direct sum of a couple of `AbelianNumber`s or `AbelianNumbers`es is defined by the direct sum of the bases of the Hilbert spaces they represent. Therefore, the input `AbelianNumber`s or `AbelianNumbers`es must be homogenous. Inhomogenous 'AbelianNumber's must be direct producted first to ensure homogenity before the direct sum.
    2. Apparently, the dimension of the result equals the summation of those of the inputs, which means, even for `AbelianNumber`s, the result will be naturally an `AbelianNumbers` because the dimension of the result is larger than 1.
    3. Signs of `AbelianNumber`s or `AbelianNumbers`es can be provided when getting their direct sums.
"""
function Base.union(qns::Vararg{<:AbelianNumber,N};signs::NTuple{N,Int}=ntuple(i->1,N)) where N
    AbelianNumbers('G',[sign==1 ? qn : -qn for (sign,qn) in zip(signs,qns)],fill(1,qns|>length),qnscounts)
end
function Base.union(qnses::Vararg{AbelianNumbers{QN},N};signs::NTuple{N,Int}=ntuple(i->1,N)) where {N,QN<:AbelianNumber}
    lengths=NTuple{N,Int}(length(qns) for qns in qnses)
    contents=Vector{QN}(undef,sum(lengths))
    indptr=zeros(Int,sum(lengths)+1)
    count=0
    for (sign,qns,length) in zip(signs,qnses,lengths)
        for i=1:length
            @inbounds contents[count+i]=sign==1 ? qns.contents[i] : -qns.contents[i]
            @inbounds indptr[count+i+1]=indptr[count+i]+qns.indptr[i+1]-qns.indptr[i]
        end
        count+=length
    end
    AbelianNumbers('G',contents,indptr,qnsindptr)
end

"""
    kron(::Type{QN},qn1::AbelianNumber,qn2::AbelianNumber) where QN<:AbelianNumber -> QN
    kron(qns::Vararg{<:AbelianNumber,N};signs::NTuple{N,Int}=ntuple(i->1,N)) where N -> eltype(qns)
    kron(qnses::Vararg{AbelianNumbers{QN},N};signs::NTuple{N,Int}=ntuple(i->1,N)) where {N,QN<:AbelianNumber} -> AbelianNumbers{QN}

Get the direct product of some `AbelianNumber`s or `AbelianNumbers`es.
!!! note
    1. Physically, the direct product of a couple of `AbelianNumber`s or `AbelianNumbers`es are defined by the direct product of the bases of the Hilbert spaces they represent. Therefore, `AbelianNumbers` with differenct types or `AbelianNumbers`es with differenct eltypes are allowed to be direct producted in principle. However, for simplicity, we only implement a method which handle the situation of two `AbelianNumber`s with differenct types. The type of the result should be provided as the first parameter. Note that in this situation, the `fieldnames` and `periods` of the result type must be exactly equal to the flattened fieldnames and periods of the two input `AbelianNumber`s, which means, even the order of the input `AbelianNumber`s matters.
    2. Apparently, the dimension of the result equals the product of those of the inputs. Therefore, the direct product of `AbelianNumber`s is also an `AbelianNumber` since its dimension is still one.
    3. For other situations except the one mentioned in Note.1, the input `AbelianNumber`s or `AbelianNumbers`es must be homogenous. Meanwhile, signs can also be provided for these situations. Note that each quantum number in the contents of the result is obtained by a summation of the corresponding quanum numbers out of the inputs with the correct signs. This is a direct observation of the Abelian nature of our quantum numbers.
"""
function Base.kron(::Type{QN},qn1::AbelianNumber,qn2::AbelianNumber) where QN<:AbelianNumber
    qnnames,qn1names,qn2names=QN|>fieldnames,qn1|>typeof|>fieldnames,qn2|>typeof|>fieldnames
    qnperiods,qn1periods,qn2periods=QN|>periods,qn1|>typeof|>periods,qn2|>typeof|>periods
    @assert qnnames==NTuple{QN|>length,Symbol}(flatten((qn1names,qn2names))) "kron error: fieldnames not match ($(qnnames),$(qn1names),$(qn2names))."
    @assert qnperiods==NTuple{QN|>length,Float}(flatten((qn1periods,qn2periods))) "kron error: periods not match ($(qnperiods),$(qn1periods),$(qn2periods))."
    QN(qn1...,qn2...)
end
Base.kron(qns::Vararg{<:AbelianNumber,N};signs::NTuple{N,Int}=ntuple(i->1,N)) where N=sum(sign==1 ? qn : -qn for (sign,qn) in zip(signs,qns))
function Base.kron(qnses::Vararg{AbelianNumbers{QN},N};signs::NTuple{N,Int}=ntuple(i->1,N)) where {N,QN<:AbelianNumber}
    lengths=NTuple{N,Int}(i<N ? dimension(qns) : length(qns) for (i,qns) in enumerate(qnses))
    contents=Vector{QN}(undef,prod(lengths))
    indptr=zeros(Int,prod(lengths)+1)
    cache=Vector{QN}(undef,N)
    @inbounds expansions=NTuple{N-1,Vector{Int}}(expand(qnses[i],qnsindices) for i=1:(N-1))
    for (i,indices) in enumerate(product(NTuple{N,UnitRange{Int64}}(1:length for length in reverse(lengths))...))
        for (j,(sign,qns,index)) in enumerate(zip(signs,qnses,reverse(indices)))
            @inbounds pos=j<N ? expansions[j][index] : index
            @inbounds cache[j]=sign==1 ? qns.contents[pos] : -qns.contents[pos]
        end
        @inbounds indptr[i+1]=indptr[i]+qnses[end].indptr[indices[1]+1]-qnses[end].indptr[indices[1]]
        @inbounds contents[i]=sum(cache)
    end
    AbelianNumbers('G',contents,indptr,qnsindptr)
end

"""
    ukron(qnses::Vararg{AbelianNumbers{QN},N};signs::NTuple{N,Int}=ntuple(i->1,N)) where {N,QN<:AbelianNumber} -> AbelianNumbers{QN},Dict{QN,Dict{NTuple{N,QN},UnitRange{Int}}}

Unitary Kronecker product of several `AbelianNumbers`es. The product result as well as the records of the product will be returned.
!!! note
    1. All input `AbelianNumbers` must be 'U' formed or 'C' formed.
    2. Since duplicate quantum number are not allowed in 'U' formed and 'C' formed `AbelianNumbers`es, in general, there exists a merge process of duplicate quantum numbers in the product result. Therefore, records are needed to keep track of this process, which will be returned along with the product result. The records are stored in a `Dict{QN,Dict{NTuple{N,QN},UnitRange{Int}}}` typed dict, in which, for each unduplicate quantum number `qn` in the product result, there exist a record `Dict((qn₁,qn₂,...)=>start:stop,...)` telling what quantum numbers `(qn₁,qn₂,...)` a mereged duplicate `qn` comes from and what slice `start:stop` this merged duplicate corresponds.
"""
function ukron(qnses::Vararg{AbelianNumbers{QN},N};signs::NTuple{N,Int}=ntuple(i->1,N)) where {N,QN<:AbelianNumber}
    @assert all(qns.form=='U' || qns.form=='C' for qns in qnses) "ukron error: all input qnses should be 'U' formed or 'C' formed."
    lengths=NTuple{N,Int}(length(qns) for qns in qnses)
    cache=Vector{QN}(undef,N)
    container=OrderedDict{QN,Int}()
    records=Dict{QN,Dict{NTuple{N,QN},UnitRange{Int}}}()
    for indices in product(NTuple{N,UnitRange{Int64}}(1:length for length in reverse(lengths))...)
        qn,count=QN|>zero,1
        for (j,(sign,qns,index)) in enumerate(zip(signs,qnses,reverse(indices)))
            @inbounds cache[j]=qns.contents[index]
            @inbounds qn=(sign==1 ? qn+cache[j] : qn-cache[j])
            @inbounds count=count*(qns.indptr[index+1]-qns.indptr[index])
        end
        container[qn]=get(container,qn,0)+count
        !haskey(records,qn) && (records[qn]=Dict{NTuple{N,QN},UnitRange{Int}}())
        records[qn][NTuple{N,QN}(cache)]=(container[qn]-count+1):container[qn]
    end
    contents=QN[qn for qn in keys(container)]
    sort!(contents)
    indptr=zeros(Int,length(container)+1)
    for (i,qn) in enumerate(contents)
        @inbounds indptr[i+1]=indptr[i]+container[qn]
    end
    AbelianNumbers('C',contents,indptr,qnsindptr),records
end

struct QnsBruteForce <: QnsProtocol end
struct QnsMonteCarlo <: QnsProtocol end
"""
    qnsbruteforce

Indicate that [`decompose`](@ref) uses the brute force method.
"""
const qnsbruteforce=QnsBruteForce()

"""
    qnsmontecarlo

Indicate that [`decompose`](@ref) uses the Monte Carlo method.
"""
const qnsmontecarlo=QnsMonteCarlo()

"""
    decompose(qnses::NTuple{N,AbelianNumbers{QN}},target::QN,signs::NTuple{N,Int},::QnsBruteForce;nmax::Int=20) where {N,QN<:AbelianNumber} -> Vector{NTuple{N,Int}}
    decompose(qnses::NTuple{N,AbelianNumbers{QN}},target::QN,signs::NTuple{N,Int},::QnsMonteCarlo;nmax::Int=20) where {N,QN<:AbelianNumber} -> Vector{NTuple{N,Int}}

Find a couple of decompositions of `target` with respect to `qnses`.
!!! note
    A tuple of integers `(i₁,i₂,...)` is called a decomposition of a given `target` with respect to the given `qnses` if and only if they satisfy the "decomposition rule":
    ```math
    \\sum_\\text{j} \\text{signs}[\\text{j}]\\times\\text{qnses}[\\text{j}][\\text{i}_{\\text{j}}]==\\text{target}
    ```
    This equation is in fact a kind of a set of restricted [linear Diophantine equations](https://en.wikipedia.org/wiki/Diophantine_equation#Linear_Diophantine_equations). Indeed, our quantum numbers are always discrete Abelian ones and all instances of a concrete `AbelianNumber` forms a [module](https://en.wikipedia.org/wiki/Module_(mathematics)) over the [ring](https://en.wikipedia.org/wiki/Ring_(mathematics)) of integers. Therefore, each quantum number can be represented as a integral multiple of the unit element of the Abelian module, which results in the final reduction of the above equation to a set of linear Diophantine equations. Then finding a decomposition is equivalent to find a solution of the reduced linear Diophantine equations, with the restriction that the quantum numbers constructed from the solution should be in the corresponding `qnses`. Here we provide two methods to find such decompositions, one is by brute force ([`qnsbruteforce`](@ref) case), and the other is by Monte Carlo simultatioins ([`qnsmontecarlo`](@ref) case).
"""
function decompose(qnses::NTuple{N,AbelianNumbers{QN}},target::QN,signs::NTuple{N,Int},::QnsBruteForce;nmax::Int=20) where {N,QN<:AbelianNumber}
    result=Set{NTuple{N,Int}}()
    cache=Vector{Int}(undef,N)
    dimensions=NTuple{N,Int}(dimension(qns) for qns in reverse(qnses))
    indices=findall(target,kron(qnses...,signs=signs),qnsexpansion)
    nmax<length(indices) && (shuffle!(MersenneTwister(),indices);indices=@views indices[1:nmax])
    for index in indices
        for (i,dimension) in enumerate(dimensions)
            @inbounds cache[end+1-i]=(index-1)%dimension+1
            @inbounds index=(index-1)÷dimension+1
        end
        push!(result,NTuple{N,Int}(cache))
    end
    return collect(NTuple{N,Int},result)
end
function decompose(qnses::NTuple{N,AbelianNumbers{QN}},target::QN,signs::NTuple{N,Int},::QnsMonteCarlo;nmax::Int=20) where {N,QN<:AbelianNumber}
    seed!()
    result=Set{NTuple{N,Int}}()
    expansions=Vector{QN}[expand(qns,qnscontents) for qns in qnses]
    dimensions=NTuple{N,Int}(dimension(qns) for qns in qnses)
    diff=indices->norm(sum(sign==1 ? expansion[index] : -expansion[index] for (sign,expansion,index) in zip(signs,expansions,indices))-target)
    count=1
    while true
        oldindices,newindices=Int[rand(1:dimension) for dimension in dimensions],ones(Int,N)
        olddiff,newdiff=diff(oldindices),diff(newindices)
        while newdiff>0
            pos=rand(1:N)
            index=rand(1:dimensions[pos]-1)
            @inbounds newindices[pos]=index<newindices[pos] ? index : index+1
            newdiff=diff(newindices)
            if newdiff<=olddiff || exp(olddiff-newdiff)>rand()
                oldindices[:]=newindices
                olddiff=newdiff
            end
            newindices[:]=oldindices
        end
        count=count+1
        push!(result,NTuple{N,Int}(newindices))
        (length(result)>=nmax || count>nmax*5) && break
    end
    return collect(NTuple{N,Int},result)
end

"""
    SQNS(S::Real) -> AbelianNumbers{SQN}

Construct the `AbelianNumbers` of the Hilbert space of a signle spin `S`.
"""
SQNS(S::Real)=AbelianNumbers('C',[SQN(sz) for sz=-S:S],collect(0:Int(2*S+1)),qnsindptr)

"""
    PQNS(N::Real) -> AbelianNumbers{PQN}

Construct the `AbelianNumbers` of the Hilbert space of a single-particle state with at most `N` identical particles.
"""
PQNS(N::Real)=AbelianNumbers('C',[PQN(np) for np=0:N],collect(0:Int(N)+1),qnsindptr)

"""
    SzPQNS(Sz::Real) -> AbelianNumbers{SPQN}

Construct the `AbelianNumbers` of the Hilbert space of a single-paritcle state with at most one particle whose spin-z component is `Sz`.
"""
SzPQNS(Sz::Real)=AbelianNumbers('C',[SPQN(0.0,0.0),SPQN(1.0,Sz)],[0,1,2],qnsindptr)

"""
    SPQNS(S::Real) -> AbelianNumbers{SPQN}

Construct the `AbelianNumbers` of the Hilbert space of a single site with internal degrees of freedom that can be ascribed to a spin `S`.
"""
function SPQNS(S::Real)
    sqns=[SQN(sz) for sz=-S:S]
    contents=[SPQN(0.0,0.0)]
    for n=1:length(sqns)
        pn=PQN(n*1.0)
        for spins in Combinations{n}(sqns)
            push!(contents,kron(SPQN,pn,sum(spins)))
        end
    end
    return sort(AbelianNumbers('G',contents,collect(0:length(contents)),qnsindptr))[1]
end

"""
    Z2QNS() -> AbelianNumbers{Z2QN}

Construct the `AbelianNumbers` of a ``Z_2`` Hilbert space.
"""
Z2QNS()=AbelianNumbers('C',[Z2QN(0.0),Z2QN(1.0)],[0,1,2],qnsindptr)

end #module


text:
#!/usr/bin/env julia

using Pkg
haskey(Pkg.installed(), "NarrativeTest") || Pkg.clone("https://github.com/rbt-lang/NarrativeTest.jl")

using DataKnots
using NarrativeTest

# Ignore the difference in the output of `print(Int)` between 32-bit and 64-bit platforms.
push!(NarrativeTest.EXPECTMAP, r"Int64" => s"Int(32|64)")

# Normalize printing of `Vector{Bool}`.
Base.show(io::IO, b::Bool) = print(io, get(io, :typeinfo, Any) === Bool ? (b ? "1" : "0") : (b ? "true" : "false"))

# Normalize `LoadError` output under 1.2.
if VERSION >= v"1.2.0-DEV"
    function Base.showerror(io::IO, ex::LoadError, bt; backtrace=true)
        print(io, "LoadError: ")
        showerror(io, ex.error, bt, backtrace=backtrace)
        print(io, "\nin expression starting at $(ex.file):$(ex.line)")
    end
end

# Normalize printing of Enum values.
if VERSION < v"1.2.0-DEV"
    Base.show(io::IO, c::DataKnots.Cardinality) =
        print(io, Symbol(c))
end

# Set the width to 72 so that MD->PDF via pandoc
ENV["COLUMNS"] = "72"

package_path(x) = relpath(joinpath(dirname(abspath(PROGRAM_FILE)), "..", x))
args = !isempty(ARGS) ? ARGS : package_path.(["doc/src", "README.md"])
exit(!runtests(args))


text:
using PlyIO

ply = Ply()
push!(ply, PlyComment("An example ply file"))

nverts = 1000

# Random vertices with position and color
vertex = PlyElement("vertex",
                    ArrayProperty("x", randn(nverts)),
                    ArrayProperty("y", randn(nverts)),
                    ArrayProperty("z", randn(nverts)),
                    ArrayProperty("r", rand(nverts)),
                    ArrayProperty("g", rand(nverts)),
                    ArrayProperty("b", rand(nverts)))
push!(ply, vertex)

# Some triangular faces
vertex_index = ListProperty("vertex_index", Int32, Int32)
for i=1:nverts
   push!(vertex_index, rand(0:nverts-1,3))
end
push!(ply, PlyElement("face", vertex_index))

# Some edges
vertex_index = ListProperty("vertex_index", Int32, Int32)
for i=1:nverts
   push!(vertex_index, rand(0:nverts-1,2))
end
push!(ply, PlyElement("edge", vertex_index))

# For the sake of the example, ascii format is used, the default binary mode is faster.
save_ply(ply, "example1.ply", ascii=true)


text:
"""
    FRAPMesh{T} <: Mesh{T}

A QBD-RAP discretisation scheme for a DiscretisedFluidQueue.

# Arguments
- `nodes`: the cell edges
- `me`: the MatrixExponential used to approximate model the fluid queue on each cell.
"""
struct FRAPMesh{T} <: Mesh{T}
    nodes::T
    me::AbstractMatrixExponential
    function FRAPMesh{T}(nodes::T,me::AbstractMatrixExponential) where T
        (nodes[1]!=0.0)&&throw(DomainError("first node must be 0"))
        return new{T}(nodes,me)
    end
end
FRAPMesh(nodes::T,me::AbstractMatrixExponential) where T =
    FRAPMesh{T}(nodes::T,me::AbstractMatrixExponential)
FRAPMesh{T}(nodes::T,n_bases::Int) where T = FRAPMesh{T}(nodes,build_me(cme_params[n_bases]))
FRAPMesh(nodes::T,n_bases::Int) where T = FRAPMesh{T}(nodes,build_me(cme_params[n_bases]))

"""

    n_bases_per_cell(mesh::FRAPMesh)

Number of bases in a cell
"""
n_bases_per_cell(mesh::FRAPMesh) = _order(mesh.me)

"""

    cell_nodes(mesh::FRAPMesh)

The cell centre
"""
cell_nodes(mesh::FRAPMesh) = Array(((mesh.nodes[1:end-1] + mesh.nodes[2:end]) / 2 )')

function build_lazy_generator(
    dq::DiscretisedFluidQueue{FRAPMesh{T}};
    v::Bool=false,
) where T
    me = dq.mesh.me
    blocks = (me.s*me.a, me.S, me.s*me.a)
    boundary_flux = OneBoundaryFlux(me.s[:],me.a[:])
    D = me.D
    out = LazyGenerator(dq,blocks,boundary_flux,D)
    v && println("UPDATE: LazyGenerator object created with keys ", keys(out))
    return out
end


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#   Types and structures.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

export PrettyTablesConf

"""
    T_BACKENDS

Types that define the supported backends.
"""
const T_BACKENDS = Union{Val{:auto}, Val{:text}, Val{:html}, Val{:latex}}

"""
    struct ColumnTable

This structure helps to access elements that comply with the column access
specification of Tables.jl.
"""
struct ColumnTable
    data::Any                    # .............................. Original table
    table::Any                   # ....... Table converted using `Tables.column`
    column_names::Vector{Symbol} # ................................ Column names
    size::Tuple{Int, Int}        # ........................... Size of the table
end

"""
    struct RowTable

This structure helps to access elements that comply with the row access
specification of Tables.jl.
"""
struct RowTable
    data::Any                    # .............................. Original table
    table::Any                   # ......... Table converted using `Tables.rows`
    column_names::Vector{Symbol} # ................................ Column names
    size::Tuple{Int, Int}        # ........................... Size of the table
end

"""
    struct PrintInfo

This structure stores the information required so that the backends can print
the tables.
"""
struct PrintInfo
    data::Any
    header::Any
    id_cols::AbstractVector{Int}
    id_rows::AbstractVector{Int}
    num_rows::Int
    num_cols::Int
    num_printed_cols::Int
    num_printed_rows::Int
    header_num_rows::Int
    header_num_cols::Int
    show_row_number::Bool
    row_number_column_title::String
    show_row_names::Bool
    row_names::Any
    row_name_alignment::Symbol
    row_name_column_title::String
    alignment::Vector{Symbol}
    cell_alignment::Ref{Any}
    formatters::Ref{Any}
    compact_printing::Bool
    title::String
    title_alignment::Symbol
    header_alignment::Vector{Symbol}
    header_cell_alignment::Ref{Any}
    cell_first_line_only::Bool
    renderer::Union{Val{:print}, Val{:show}}
    limit_printing::Bool
end

"""
    struct PrettyTablesConf

Type of the object that holds a pre-defined set of configurations for
PrettyTables.jl.
"""
struct PrettyTablesConf
    confs::Dict{Symbol, Any}
end

PrettyTablesConf() = PrettyTablesConf(Dict{Symbol, Any}())