text:
# AEdgeMap{T} defined in newtwork_sinterface.jl
valtype(m::AEdgeMap{T}) where {T} = T

"""
    mutable struct EdgeMap{G <: AGraphOrDiGraph, T, D} <: AEdgeMap{T}
        g::G
        vtype::Type{T}
        data::D
    end

Type implementing an edge map. The underlying container `data` can be a dictionary,
a matrix or a vector (for graphs with indexed edges).

    EdgeMap{T}(g, ::Type{T})

Returns a map that associates values of type `T`
to the vertices of  graph `g`. The underlying storage structures is chosen
accordingly.

    EdgeMap(g, data)

Construct a EdgeMap with `data` as underlying storage.
The storage type can be a matrix or an associative `edg => val` type or
a vector for graph with indexed edges.

    EdgeMap(g, f)

Construct an edge map with value `f(e)` for each `e` in `edges(g)`.
"""
mutable struct EdgeMap{G<:AGraphOrDiGraph, T, D} <: AEdgeMap{T}
    g::G
    vtype::Type{T}
    data::D
end

EdgeMap(g::AGraphOrDiGraph, d::AbstractMatrix{T}) where {T} = EdgeMap(g, T, d)
EdgeMap(g::AGraphOrDiGraph, d::AbstractVector{T}) where {T} = EdgeMap(g, T, d)
EdgeMap(g::AGraphOrDiGraph, d::Dict{Int, T}) where {T} = EdgeMap(g, T, d)
EdgeMap(g::AGraphOrDiGraph, d::Dict{E, T}) where {T,E<:AEdge} = EdgeMap(g, T, d)

function EdgeMap(g::AGraphOrDiGraph, ::Type{T}) where T
    E = edgetype(g)
    if E <: AIndexedEdge
        # data = Vector{T}(ne(g))
        data = Dict{Int,T}()
    else
        data = Dict{E,T}()
    end
    return EdgeMap(g, T, data)
end

function EdgeMap(g::AGraphOrDiGraph, f::Function)
    E = edgetype(g)
    if E <: AIndexedEdge
        data_ = [f(e) for e in edges(g)]
        # since edges(g) doesn't iterate according to index
        # have to permute after
        data = similar(data_)
        for (i, e) in enumerate(edges(g))
            data[idx(e)] = data_[i]
        end
        T = eltype(data)
    else
        data = Dict(e => f(e) for e in edges(g))
        T = valtype(data)
    end
    return EdgeMap(g, T, data)
end

length(m::EdgeMap) = length(m.data)

### ALL DATA
# matrix interface
getindex(m::EdgeMap, i::Integer, j::Integer) = getindex(m, edge(m.g, i, j))
setindex!(m::EdgeMap, x, i::Integer, j::Integer) = setindex!(m, x, edge(m.g, i, j))
haskey(m::EdgeMap, i::Integer, j::Integer) = haskey(m, edge(m.g, i, j))
size(m::EdgeMap) = (nv(m.g), nv(m.g))
size(m::EdgeMap, i::Integer)::Int = 1 <= i <= 2 ? nv(m.g) : error("wrong dimension")

### MATRIX DATA
_sort(i, j) = i <= j ? (i, j) : (j , i)
# Associative interface
getindex(m::EdgeMap{G,T,D}, e::AEdge) where {G<:AGraphOrDiGraph,T,D<:AbstractMatrix} =
    getindex(m, src(e), dst(e))
setindex!(m::EdgeMap{G,T,D}, x, e::AEdge) where {G<:AGraphOrDiGraph,T,D<:AbstractMatrix} =
    setindex!(m, x, src(e), dst(e))
get(m::EdgeMap{G,T,D}, e::AEdge, x) where {G<:AGraph,T,D<:AbstractMatrix} =
    get(m.data, _sort(Int(src(e)), Int(dst(e))), x)
get(m::EdgeMap{G,T,D}, e::AEdge, x) where {G<:ADiGraph,T,D<:AbstractMatrix} =
    get(m.data, (Int(src(e)), Int(dst(e))), x)
haskey(m::EdgeMap{G,T,D}, e::AEdge) where {G<:AGraphOrDiGraph,T,D<:AbstractMatrix} =
     haskey(m, src(e), dst(e))
haskey(m::EdgeMap{G,T,D}, i::Integer, j::Integer) where {G<:AGraphOrDiGraph,T,D<:AbstractMatrix} =
  (1 <= i <= size(m.data, 1)) && (1 <= j <= size(m.data, 1))

# matrix interface
getindex(m::EdgeMap{G,T,D}, i::Integer, j::Integer) where {G<:AGraph,T,D<:AbstractMatrix} =
    getindex(m.data, _sort(i, j)...)
getindex(m::EdgeMap{G,T,D}, i::Integer, j::Integer) where {G<:ADiGraph,T,D<:AbstractMatrix} =
    getindex(m.data, i, j)
setindex!(m::EdgeMap{G,T,D}, x, i::Integer, j::Integer) where {G<:AGraph,T,D<:AbstractMatrix} =
    setindex!(m.data, x, _sort(i, j)...)
setindex!(m::EdgeMap{G,T,D}, x, i::Integer, j::Integer) where {G<:ADiGraph,T,D<:AbstractMatrix} =
    setindex!(m.data, x, i, j)

### VECTOR DATA (only indexed edges)
# Associative interface
getindex(m::EdgeMap{G,T,D}, e::AIndexedEdge) where {G<:AGraphOrDiGraph,T,D<:AbstractVector} =
    getindex(m.data, idx(e))
setindex!(m::EdgeMap{G,T,D}, x, e::AIndexedEdge) where {G<:AGraphOrDiGraph,T,D<:AbstractVector} =
    setindex!(m.data, x, idx(e))
haskey(m::EdgeMap{G,T,D}, e::AIndexedEdge) where {G<:AGraphOrDiGraph,T,D<:AbstractVector} =
    1 <= idx(e) <= length(m.data)
get(m::EdgeMap{G,T,D}, e::AEdge, x) where {G<:AGraphOrDiGraph,T,D<:AbstractVector} =
    get(m.data, idx(e), x)

# TODO allow one dimensional indexing?
# it can be bugprone
# getindex{G<:AGraphOrDiGraph,T,D<:AbstractVector}(m::EdgeMap{G,T,D}, idx::Integer) =
#     getindex(m.data, idx)
# setindex!{G<:AGraphOrDiGraph,T,D<:AbstractVector}(m::EdgeMap{G,T,D}, x, idx::Integer) =
#     setindex!(m.data, x, idx)

### Dict{Int,T} DATA
# Associative interface
getindex(m::EdgeMap{G,T,Dict{Int,T}}, e::AIndexedEdge) where {G<:AGraphOrDiGraph,T} = getindex(m.data, idx(e))
setindex!(m::EdgeMap{G,T,Dict{Int,T}}, x, e::AIndexedEdge) where {G<:AGraphOrDiGraph,T} = setindex!(m.data, x, idx(e))
get(m::EdgeMap{G,T,Dict{Int,T}}, e::AIndexedEdge, x) where {G<:AGraphOrDiGraph,T} = get(m.data, idx(e), x)
haskey(m::EdgeMap{G,T,Dict{Int,T}}, e::AIndexedEdge) where {G<:AGraphOrDiGraph,T} = haskey(m.data, idx(e))

### Dict{E,T} DATA
# Associative interface
getindex(m::EdgeMap{G,T,Dict{E,T}}, e::E) where {G<:AGraphOrDiGraph,T,E<:AEdge} = getindex(m.data, e)
setindex!(m::EdgeMap{G,T,Dict{E,T}}, x, e::E) where {G<:AGraphOrDiGraph,T,E<:AEdge} = setindex!(m.data, x, e)
get(m::EdgeMap{G,T,Dict{E,T}}, e::E, x) where {G<:AGraphOrDiGraph,T,E<:AEdge} = get(m.data, e, x)
haskey(m::EdgeMap{G,T,Dict{E,T}}, e::E) where {G<:AGraphOrDiGraph,T,E<:AEdge} = haskey(m.data, e)

getindex(m::EdgeMap{G,T,Dict{E,T}}, e::AEdge) where {G<:AGraphOrDiGraph,T,E<:AEdge} = getindex(m.data, E(src(e),dst(e)))
setindex!(m::EdgeMap{G,T,Dict{E,T}}, x, e::AEdge) where {G<:AGraphOrDiGraph,T,E<:AEdge} = setindex!(m.data, x, E(src(e),dst(e)))
get(m::EdgeMap{G,T,Dict{E,T}}, e::AEdge, x) where {G<:AGraphOrDiGraph,T,E<:AEdge} = get(m.data, E(src(e),dst(e)), x)
haskey(m::EdgeMap{G,T,Dict{E,T}}, e::AEdge) where {G<:AGraphOrDiGraph,T,E<:AEdge} = haskey(m.data, E(src(e),dst(e)))
####
values(m::EdgeMap{G,T,D}) where {G,T,D<:Dict} = values(m.data)
values(m::EdgeMap{G,T,D}) where {G,T,D<:Array} = m.data
values(m::EdgeMap{G,T,D}) where {G,T,D<:AbstractSparseMatrix} = nonzeros(m.data)

==(m1::EdgeMap, m2::EdgeMap) = m1.data == m2.data

"""
    struct ConstEdgeMap{T} <: SimpleEdgeMap{T}
        val::T
    end

A type representing a constant vector map.
Any attempt to change the internal value, e.g. `emap[u,v] = 4`, will
fail silently.
"""
struct ConstEdgeMap{T} <: AEdgeMap{T}
    val::T
end

ConstEdgeMap(g::AGraphOrDiGraph, x) = ConstEdgeMap(x)

length(m::ConstEdgeMap) = typemax(Int)
getindex(m::ConstEdgeMap, e::AEdge) = m.val

setindex!(m::ConstEdgeMap, x, i::Integer, j::Integer) = nothing #TODO not throwin since it is used as
                                                                # a dummy map
# setindex!(m::ConstEdgeMap, x, i::Integer, j::Integer) = error("Cannot assign to ConstEdgeMap")
getindex(m::ConstEdgeMap, i::Integer, j::Integer) = m.val
setindex!(m::ConstEdgeMap, x, e::AEdge) = nothing
# setindex!(m::ConstEdgeMap, x, e::AEdge) = error("Cannot assign to ConstEdgeMap")
get(m::ConstEdgeMap, e::AEdge, x) = m.val

values(m::ConstEdgeMap) = [m.val]

"""
    edgemap2adjlist(emap)

Returns a vector of vectors containing the values of the edge map `emap` on graph `g`
following the same ordering of [`adjacency_list`](@ref)`(g)`.
"""
edgemap2adjlist(m::AEdgeMap) = [[m[e] for e in out_edges(m.g, i)] for i=1:nv(m.g)]

function Base.Matrix(emap::EdgeMap)
    g = emap.g
    M = zeros(valtype(emap), nv(g), nv(g))
    return fill_mat_from_map!(M, g, emap)
end

function sparse(emap::EdgeMap)
    g = emap.g
    M = spzeros(valtype(emap), nv(g), nv(g))
    return fill_mat_from_map!(M, g, emap)
end

function fill_mat_from_map!(M, g::ADiGraph, emap::AEdgeMap)
    for e in edges(g)
        u, v = src(e), dst(e)
        M[u,v] = emap[e]
    end
    return M
end

function fill_mat_from_map!(M, g::AGraph, emap::AEdgeMap)
    for e in edges(g)
        u, v = src(e), dst(e)
        M[u,v] = M[v,u] = emap[e]
    end
    return M
end

"""
    weights(g)

Returns an edge map containing the "weights" associated to edges.
For simple graphs, the return value is ConstEdgeMap(g, 1).
For networks, returns the "weights" edge property if defined,
otherwise the constant map.

Notice that the edge map returned by `weights` is the
default value for the edge weights used in many flow and
distance on graph algorithms.
"""
weights(g::AGraphOrDiGraph) = ConstEdgeMap(g, 1)

weights(g::ANetOrDiNet) = has_eprop(g, "weights") ?
        eprop(g, "weights") : ConstEdgeMap(g, 1)


text:
export whitespace_char, whitespace_maybe, whitespace
export horizontal_space_char, horizontal_space_maybe, horizontal_space
export vertical_space_char, vertical_space_maybe, vertical_space
export space, space_maybe, whitespace_newline
export trim, @trimmed
export newline, inline
export at_linestart, at_lineend

char_label = Dict(
    '\U0009' => "Horizontal tab (HT)",
    '\U0020' => "Space",
    '\U00A0' => "Non-break space",
    '\U1680' => "Ogham space mark",
    '\U180E' => "Mongolian vowel separator",
    '\U2000' => "En quad",
    '\U2001' => "Em quad",
    '\U2002' => "En space",
    '\U2003' => "Em space",
    '\U2004' => "Three-per-em space",
    '\U2005' => "Four-per-em space",
    '\U2006' => "Six-per-em space",
    '\U2007' => "Figure space",
    '\U2008' => "Punctuation space",
    '\U2009' => "Thin space",
    '\U200A' => "Hair space",
    '\U202F' => "Narrow no-break space",
    '\U205F' => "Medium mathematical space",
    '\U3000' => "Ideographic space",
    '\u2028' => "Line Separator",
    '\u200e' => "Left-to-right mark",
    '\v' => "Vertical tab",
    '\u2029' => "Paragraph separator",
    '\u85' => "Next line",
    '\u200f' => "Right-to-left mark",
    '\f' => "Form feed",
    '\U000A' => "Linefeed (LF)",
    '\U000B' => "Vertical tab (VT)",
    '\U000C' => "Form feed (FF)",
    '\U000D' => "Carriage return (CR)",
    '\U0085' => "Next line (NEL)",
    '\U2028' => "Line separator",
    '\U2029' => "Paragraph separator"
)

char_label_table(x::Union{ValueIn,ValueNotIn}) =
    char_label_table(x.sets)
function char_label_table(x)
    chars = sort(collect(x))
    rs = [ escape_string(repr(c)) for c in chars ]
    labels = [ get(char_label, c, "") for c in chars ]
    nr = maximum(length.(rs))
    n = maximum(length.(labels))
    println("| ",lpad("Char",nr)," | ", lpad("",n), " |")
    println("|-","-"^nr, "-|-", "-"^n,      "-|")
    for (c,l) in zip(rs, labels)
        print("| ")
        print(lpad(c,nr))
        println(" | ", lpad(l,n), " |")
    end
end

whitespace_char = CharIn(
    "[:space:]",
    " \t\U0085\U200E\U200F\U2028\U2029"*"\U2029\U000C\U000B")

@with_names begin
    whitespace_maybe = !Atomic(Repeat(whitespace_char))
    whitespace = !Atomic(Repeat1(whitespace_char))
end

"""
    whitespace_char  = re"[[:space:]]"
    whitespace_maybe = re"(?>[[:space:]]*)"
    whitespace       = re"(?>[[:space:]]+)"

```jldoc
julia> CombinedParsers.char_label_table(whitespace_char)
|      Char |                     |
|-----------|---------------------|
|     '\\t' | Horizontal tab (HT) |
|     '\\v' |   Vertical tab (VT) |
|     '\\f' |      Form feed (FF) |
|       ' ' |               Space |
|   '\\u85' |     Next line (NEL) |
| '\\u200e' |  Left-to-right mark |
| '\\u200f' |  Right-to-left mark |
| '\\u2028' |      Line separator |
| '\\u2029' | Paragraph separator |
```
"""
whitespace_char, whitespace_maybe, whitespace

horizontal_space_char=CharIn("\\h",
    '\U0009', # "Horizontal tab (HT)"),
    '\U0020', # "Space"),
    '\U00A0', # "Non-break space"),
    '\U1680', # "Ogham space mark"),
    '\U180E', # "Mongolian vowel separator"),
    '\U2000', # "En quad"),
    '\U2001', # "Em quad"),
    '\U2002', # "En space"),
    '\U2003', # "Em space"),
    '\U2004', # "Three-per-em space"),
    '\U2005', # "Four-per-em space"),
    '\U2006', # "Six-per-em space"),
    '\U2007', # "Figure space"),
    '\U2008', # "Punctuation space"),
    '\U2009', # "Thin space"),
    '\U200A', # "Hair space"),
    '\U202F', # "Narrow no-break space"),
    '\U205F', # "Medium mathematical space"),
    '\U3000' # "Ideographic space"))
)

horizontal_space_maybe = with_name(:horizontal_space_maybe,
                                   !Atomic(Repeat(horizontal_space_char)),
                                   "\\h*")
horizontal_space = with_name(:horizontal_space,
                             !Atomic(Repeat1(horizontal_space_char)),
                             "\\h+")

"""
    horizontal_space_char  = re"[\\h]"
    horizontal_space_maybe = re"(?>[\\h]*)"
    horizontal_space       = re"(?>[\\h]+)"

```jldoc
julia> CombinedParsers.char_label_table(horizontal_space_char)
|      Char |                           |
|-----------|---------------------------|
|     '\\t' |       Horizontal tab (HT) |
|       ' ' |                     Space |
|       ' ' |           Non-break space |
|       ' ' |          Ogham space mark |
| '\\u180e' | Mongolian vowel separator |
|       ' ' |                   En quad |
|       ' ' |                   Em quad |
|       ' ' |                  En space |
|       ' ' |                  Em space |
|       ' ' |        Three-per-em space |
|       ' ' |         Four-per-em space |
|       ' ' |          Six-per-em space |
|       ' ' |              Figure space |
|       ' ' |         Punctuation space |
|       ' ' |                Thin space |
|       ' ' |                Hair space |
|       ' ' |     Narrow no-break space |
|       ' ' | Medium mathematical space |
|       '　' |         Ideographic space |
```
"""
horizontal_space_char, horizontal_space_maybe, horizontal_space

whitespace_horizontal = horizontal_space

"""
    trim(p...; whitespace=horizontal_space_maybe,
               left=whitespace, right=whitespace)

Ignore whitespace `left` and `right` of `sSequence(p...)`.
"""
trim(p...; whitespace=horizontal_space_maybe, left=whitespace, right=whitespace) =
    Sequence(left, sSequence(p...), right)[2]

trimmed(x) = x
function trimmed(node::Expr)
    if node.head == :(=) && length(node.args) == 2 && isa(node.args[1], Symbol)
        node.args[2] = Expr(:call, :trim, Expr(:call, :with_name, QuoteNode(node.args[1]), node.args[2]))
    end
    if node.head != :call
        node.args = map(trimmed, node.args)
    end
    node
end

"""
    @trimmed

Create parser within `whitespace_maybe` to match the variables they are asigned to.

See also [`trim`](@ref).

```@meta
DocTestFilters = r"map\\(.+\\)"
```

so, for example
```jldoctest
julia> @trimmed foo = AnyChar()
🗄 Sequence[2]
├─ (?>[\\h]*) ValueIn |> Repeat |> ! |> Atomic
├─ . AnyValue |> with_name(:foo)
└─ (?>[\\h]*) ValueIn |> Repeat |> ! |> Atomic
::Char

julia> parse(log_names(foo),"  ab  ")
   match foo@3-4:   ab
                    ^
'a': ASCII/Unicode U+0061 (category Ll: Letter, lowercase)
```

"""
macro trimmed(block)
    esc(trimmed(block))
end

vertical_space_char=CharIn(
    "\\v",
    '\U000A', # "Linefeed (LF)"),
    '\U000B', # "Vertical tab (VT)"),
    '\U000C', # "Form feed (FF)"),
    '\U000D', # "Carriage return (CR)"),
    '\U0085', # "Next line (NEL)"),
    '\U2028', # "Line separator"),
    '\U2029') # "Paragraph separator"))

vertical_space_maybe = with_name(:vertical_space_maybe,
                                 !Atomic(Repeat(vertical_space_char)),
                                 "\\v*")
vertical_space = with_name(:vertical_space,
                           !Atomic(Repeat1(vertical_space_char)),
                           "\\v+")

"""
    vertical_space_char  = re"[\\v]"
    vertical_space_maybe = re"(?>[\\v]*)"
    vertical_space       = re"(?>[\\v]+)"

```jldoc
julia> CombinedParsers.char_label_table(vertical_space_char)
|      Char |                      |
|-----------|----------------------|
|     '\\n' |        Linefeed (LF) |
|     '\\v' |    Vertical tab (VT) |
|     '\\f' |       Form feed (FF) |
|     '\\r' | Carriage return (CR) |
|   '\\u85' |      Next line (NEL) |
| '\\u2028' |       Line separator |
| '\\u2029' |  Paragraph separator |
```
"""
vertical_space_char, vertical_space_maybe, vertical_space

"Equivalent PRCE `\\h\\v`, [`horizontal_space_char`](@ref), [`vertical_space_char`](@ref)"
space_char  = CharIn("\\h\\v",horizontal_space_char,vertical_space_char)
@with_names begin
    space_maybe = Atomic(!Repeat(space_char))
    space = Atomic(!Repeat1(CharIn("\\h\\v",space_char)))
end

@deprecate whitespace_newline space

"""
    CombinedParsers.newline
    CombinedParsers.Regexp.bsr

newlines, PCRE `\\r` backslash R (BSR).

```jldoctest
julia> CombinedParsers.bsr
\\r with_name(:bsr)
::SubString{String}

julia> CombinedParsers.bsr.parser
(?>|🗄) Either |> Atomic |> !
├─ \\r\\n
└─ [\\n\\x0b\\f\\r\\x85] ValueIn
::SubString{String}
```

"""
bsr = with_name(
    :bsr, !Atomic(Either("\r\n",
                        CharIn(raw"\n\x0b\f\r\x85", '\n','\x0b','\f','\r','\U0085', '\U2028','\U2029'))),
    "\\r");

newline = bsr

"Equivalent PRCE `\\w`: Char with unicode class `L`, `N`, or `_`."
word_char=CharIn("\\w",UnicodeClass("L","N"),'_')

"SubString of at least 1 repeated [`CombinedParsers.word_char`](@ref)."
word = with_name(:word, !Repeat1(word_char), ## "[[:alpha:] ]+"
                 "\\w+")

"Vector of at least 1 repeated [`CombinedParsers.word`](@ref)s delimited by [`CombinedParsers.whitespace_horizontal`](@ref)."
words = join(word, whitespace_horizontal) ## "[[:alpha:] ]+"

non_word_char=CharNotIn("\\W",UnicodeClass("L","N"),'_')
non_word = with_name(:non_word, !Repeat1(non_word_char),
                     "\\W+")## "[[:alpha:] ]+"

"""
    beyond_word = Either(non_word_char,AtStart(),AtEnd())

Parser part of `word_boundary`.
"""
@with_names beyond_word = Either(non_word_char,AtStart(),AtEnd())

"""
    word_boundary = re"\b"
"""
@with_names word_boundary = Either(
    Sequence(PositiveLookbehind(word_char),PositiveLookahead(beyond_word)),
    Sequence(PositiveLookbehind(beyond_word),PositiveLookahead(word_char))
)

"""
    at_linestart

```jldoctest
julia> CombinedParsers.Regexp.at_linestart
|🗄 Either |> with_name(:at_linestart)
├─ ^ AtStart
└─ (?<=🗄)) Either |> Atomic |> with_name(:bsr) |> PositiveLookbehind
   ├─ \\r\\n
   └─ [\\n\\x0b\\f\\r\\x85] ValueIn |> !
::Union{AtStart, SubString}
```

!!! note
    used in `re"^"` if `Base.PCRE.MULTILINE` is set.
"""
@with_names at_linestart = Either(AtStart(),PositiveLookbehind(bsr))
# lineend   = Either(AtEnd(),bsr)

"""
    at_lineend

```jldoctest
julia> CombinedParsers.Regexp.at_lineend
|🗄 Either |> with_name(:at_lineend)
├─ \$ AtEnd
└─ (?=(?>|🗄)) Either |> Atomic |> with_name(:bsr) |> PositiveLookahead
   ├─ \\r\\n
   └─ [\\n\\x0b\\f\\r\\x85] ValueIn |> !
::Union{AtEnd, SubString{String}}
```

!!! note
    used in `re"\$"` if `Base.PCRE.MULTILINE` is set.
"""
@with_names at_lineend   = Either(AtEnd(),PositiveLookahead(bsr))

"""
    inline = !Atomic(Repeat(NegativeLookahead(at_lineend)*AnyChar()))

See [`at_lineend`](@ref).
"""
@with_names inline = !Atomic(Repeat(NegativeLookahead(at_lineend)*AnyChar()))

hex_digit = CharIn("[:xdigit:]",'A':'F','a':'f','0':'9')
export hex_digit, integer_base
"""
    integer_base(base,mind=0,maxd=Repeat_max)

Parser matching a integer format on base `base`.

!!! note
    Uses a second Base.parse call on match.

    A custom parser could aggregate result incrementally while matching.
"""
function integer_base(base=10,mind=0,maxd=Repeat_max)
    dig = if base == 16
        hex_digit
    elseif base <= 10
        ValueIn('0':('0'+(base-1)))
    else
        error("Base $base not supported")
    end
    map(!Repeat(mind:maxd,dig)) do v
        (isempty(v) ? 0 : parse(Int,convert(String,v),base=base))::Int
    end
end

ext: jl


text:
#--------------------------------------------------------------------
# DNSS.jl
# Soham 08-2019
# Test scalar field on Minkowski spacetime
#--------------------------------------------------------------------

using Test, NLsolve, LinearAlgebra, Random, Printf, ForwardDiff
# using LaTeXStrings, Printf, PyPlot

#-----------------------------------------
# Set up initial conditions on u0 and v0
#-----------------------------------------
function psibar(u::Real, v::Real)::Real
    return  exp(-u^2 / σ^2) * sin(pi*u) + exp(-v^2 / σ^2) * sin(pi*v)
end

function computeUboundary(PS::ProductSpace{S1, S2})::NTuple{1, Field{S2}} where {S1, S2}
    return (extractUboundary(Field(PS, psibar), :incoming), )
end

function computeVboundary(PS::ProductSpace{S1, S2})::NTuple{1, Field{S1}} where {S1, S2}
    return (extractVboundary(Field(PS, psibar), :incoming), )
end

#------------------------------------------------------------
# Source term or the change of coordinate functions
#------------------------------------------------------------
function twist(uv::Vector{T})::Vector{T} where {T<:Real}
    u, v = uv
    ω =  (pi/8) * cospi(u/2) * cospi(v/2)
    ũ =  u*cos(ω) + v*sin(ω)
    ṽ = -u*sin(ω) + v*cos(ω)
    return [ũ, ṽ]
end

function dxdX(u::Real, v::Real, component::Symbol)
    if component == :dUdŨ
        return inv(ForwardDiff.jacobian(twist, [u,v]))[1,1]
    elseif component == :dUdṼ
        return inv(ForwardDiff.jacobian(twist, [u,v]))[1,2]
    elseif component == :dVdŨ
        return inv(ForwardDiff.jacobian(twist, [u,v]))[2,1]
    elseif component == :dVdṼ
        return inv(ForwardDiff.jacobian(twist, [u,v]))[2,2]
    end
end

function source(u::Real, v::Real)::Real
    t = u + v
    r = v - u
    return 20 * sinpi(t) * sinpi(r) *  exp(-r^2 / 0.1)
end

#------------------------------------------------------------
# Set up computation on each patch using a non-linear solver
#------------------------------------------------------------
function computePatch(PS::ProductSpace{S1, S2}, ubnd::NTuple{1,Field{S2}}, vbnd::NTuple{1,Field{S1}})::NTuple{1, Field{ProductSpace{S1, S2}}} where {S1, S2}
    if solver == :nonlin
        # Define all the operators we use local to the patch.
        B = incomingboundary(PS)
        I = identity(PS)
        DU, DV = derivative(PS)
        psibnd = combineUVboundary(ubnd, vbnd, :incoming)

        dUdŨ   = Field(PS, (u,v)-> dxdX(u,v, :dUdŨ))
        dVdŨ   = Field(PS, (u,v)-> dxdX(u,v, :dVdŨ))
        dUdṼ   = Field(PS, (u,v)-> dxdX(u,v, :dUdṼ))
        dVdṼ   = Field(PS, (u,v)-> dxdX(u,v, :dVdṼ))

        # Compute the residual. This is slow but works. Sets the residual in the
        # interior as well as the boundary.
        function residual(psi::NTuple{1, Field{S}})::NTuple{1, Field{S}} where {S}
            psi   = first(psi)
            # FIXME: The twisted coordinates do not work
            dpsidŨ    = dUdŨ * (DU * psi) + dVdŨ * (DV * psi)
            ddpsidŨdṼ = dUdṼ * (DU * dpsidŨ) + dVdṼ * (DV * dpsidŨ)
            F = (I - B) * ddpsidŨdṼ  + B * (psi - first(psibnd))
            return (F,)
        end

        # Residual to pass to the non-linear solver in NLsolve
        function f!(fvec::Array{T,1}, x::Array{T,1}) where {T}
            fvec[:]  = reshapeFromTuple(residual(reshapeToTuple(PS, 1, x)))
        end

        # Specify an initial guess
        function initialguess(PS::ProductSpace{S1, S2})::NTuple{1, Field{ProductSpace{S1, S2}}} where {S1, S2}
            return (Field(PS, (u,v)->1), )
        end

        # Compute the initial guess to pass to the solver
        U0 = reshapeFromTuple(initialguess(PS))

        # Call the non-linear solver
        U = nlsolve(f!, U0; method=:trust_region, autodiff=:forward, show_trace=true, ftol=1e-8, iterations=100)
        return reshapeToTuple(PS, 1, U.zero)

    elseif solver == :lin
        # Define all the operators we use local to the patch.
        B = incomingboundary(PS)
        DU, DV = derivative(PS)
        bnd = combineUVboundary(ubnd, vbnd, :incoming)
        L   = 2*DU*DV
        s   = 0 * Field(PS, source)
        b   = first(bnd)

        # Using the linearity of the system to solve the system
        # L u  = s # Wave equation with a source
        # B u =  b # Boundary conditions
        # (L + B) u = b
        U = linsolve(L + B,  s + b)
        return (U, )
    end
end

#-----------------------------------------
# Compare with analytic solution
#-----------------------------------------
function deltapsi(U::NTuple{1, Field})::NTuple{1, Field} where {S}
    # Compute the error at twice the number of points
    psi = project(first(U), prolongate(first(U).space))
    psiexact = Field(psi.space, psibar)
    return (psi - psiexact, )
end

#-----------------------------------------
# h and p convergence functions
#-----------------------------------------
function pconv(range::StepRange, nh::Int=1)
    println("Testing p convergence")
    n_ = collect(range)
    l_ = zeros(size(n_))
    for index in CartesianIndices(n_)
        np = n_[index]
        params = Parameters((np, np), (nh,nh), urange, vrange, 1)
        q = distribute(params, computePatch, computeUboundary, computeVboundary)
        l_[index] = rmse(extract(map(deltapsi, q), 1))
        @printf("  n = %i, rmse = %e\n", n_[index], l_[index])
    end
    return (n_, l_)
end

function hconv(range::StepRange, np::Int=4)
    println("Testing h convergence")
    n_ = collect(range)
    l_ = zeros(size(n_))
    for index in CartesianIndices(n_)
        nh = n_[index] = 2^n_[index]
        params = Parameters((np, np), (nh, nh), urange, vrange, 1)
        l_[index] = rmse(extract(map(deltapsi, distribute(params, computePatch, computeUboundary, computeVboundary)), 1))
        l0 = index[1] > 1 ? l_[index[1] - 1] : 1
        @printf("  n = %3i, rmse[2h] / rmse[h] = %e\n", n_[index], l0 / l_[index])
    end
    return (n_, l_)
end

#-----------------------------------------
# Setup simulation grid parameters
#-----------------------------------------
# TODO: Check if we find convergence for different sizes of domains. What
# happens if the domain size increases? Is the solution unstable for larger
# domain sizes? What we find is that one needs more spacetime elements for
# larger domain sizes to see p-convergence. For a single or fewer elements,
# the p-convergence is either non-existent, or slow.
# FIXME: In any case, we saturate around 1e-10 with all the errors. Is this due
# to the floor set by the condition number of the second derivative matrices?
# What happens if we choose to work with BigFloat?
npoints  = (20, 20)
npatches = (1, 1)
urange   = (-1.0, 1.0)
vrange   = (-1.0, 1.0)
nfields  = 1
Grid     = Parameters(npoints, npatches, urange, vrange, nfields)
solver   = :lin # Use a linear solver instead of the non-linear solver
σ        = 0.5

#-----------------------------------------
# Compute the scalar field over the grid
#-----------------------------------------
# ϕ_ = distribute(Grid, computePatch, computeUboundary, computeVboundary)
# save("output/data/minkowski_2D/sigma_0p5_np_30_nh_4_source/minkowski_2D_sigma_0p5_psi", extract(ϕ_, 1))
# save("output/data/minkowski_2D/sigma_0p5_np_30_nh_4_source/minkowski_2D_sigma_0p5_deltapsi", extract(map(deltapsi,ϕ_), 1))

#-----------------------------------------
# Now test h-p convergence of the error
# in the solution
#-----------------------------------------
np, ep = pconv(2:2:40, 12)
nh, eh = hconv(0:1:9, 4)

#-----------------------------------------
# Save all the data for solution and convergence
#-----------------------------------------
# save("output/data/minkowski_2D_sigma_0p9_psi", extract(ϕ_, 1))
# save("output/data/minkowski_2D_sigma_0p9_deltapsi_pconv", np, ep)
# save("output/data/minkowski_2D_sigma_0p9_deltapsi_hconv", nh, eh)


text:
""" Vertical circulation associated with horizontal velocities u, v. """
@inline Γᶠᶠᵃ(i, j, k, grid, u, v) = δxᶠᵃᵃ(i, j, k, grid, Δy_vᶜᶠᵃ, v) - δyᵃᶠᵃ(i, j, k, grid, Δx_uᶠᶜᵃ, u)

""" Vertical vorticity associated with horizontal velocities u, v. """
@inline ζ₃ᶠᶠᵃ(i, j, k, grid, u, v) = Γᶠᶠᵃ(i, j, k, grid, u, v) / Azᶠᶠᵃ(i, j, k, grid)


text:
struct Direction{T} end

const Left = Direction{:Left}()
const Right = Direction{:Right}()

const Beginning = Left
const Ending = Right

struct Endpoint{T, D}
    endpoint::T
    included::Bool

    function Endpoint{T, D}(ep::T, included::Bool) where {T, D}
        @assert D isa Direction
        new{T, D}(ep, included)
    end
end

const LeftEndpoint{T} = Endpoint{T, Left}
const RightEndpoint{T} = Endpoint{T, Right}

LeftEndpoint(ep::T, included::Bool) where T = LeftEndpoint{T}(ep, included)
RightEndpoint(ep::T, included::Bool) where T = RightEndpoint{T}(ep, included)

LeftEndpoint(i::AbstractInterval{T}) where T = LeftEndpoint{T}(first(i), first(inclusivity(i)))
RightEndpoint(i::AbstractInterval{T}) where T = RightEndpoint{T}(last(i), last(inclusivity(i)))

function Base.hash(x::Endpoint{T, D}, h::UInt) where {T, D}
    # Note: we shouldn't need to hash `T` as this is covered by the endpoint field.
    h = hash(:Endpoint, h)
    h = hash(D, h)
    h = hash(x.endpoint, h)
    h = hash(x.included, h)
    return h
end

# https://github.com/JuliaLang/julia/pull/26601
if VERSION >= v"0.7.0-DEV.4743"
    Base.broadcastable(e::Endpoint) = Ref(e)
end

"""
    ==(a::Endpoint, b::Endpoint) -> Bool

Determine if two endpoints are equal. When both endpoints are left or right then the points
and inclusiveness must be the same.

Checking the equality of left-endpoint and a right-endpoint is slightly more difficult. A
left-endpoint and a right-endpoint are only equal when they use the same point and are
both included. Note that left/right endpoints which are both not included are not equal
as the left-endpoint contains values below that point while the right-endpoint only contains
values that are above that point.

Visualizing two contiguous intervals can assist in understanding this logic:

    [x..y][y..z] -> RightEndpoint == LeftEndpoint
    [x..y)[y..z] -> RightEndpoint != LeftEndpoint
    [x..y](y..z] -> RightEndpoint != LeftEndpoint
    [x..y)(y..z] -> RightEndpoint != LeftEndpoint
"""
function Base.:(==)(a::Endpoint, b::Endpoint)
    a.endpoint == b.endpoint && a.included == b.included
end

function Base.:(==)(a::LeftEndpoint, b::RightEndpoint)
    a.endpoint == b.endpoint && a.included && b.included
end

function Base.:(==)(a::RightEndpoint, b::LeftEndpoint)
    b == a
end

function Base.isequal(a::Endpoint, b::Endpoint)
    isequal(a.endpoint, b.endpoint) && isequal(a.included, b.included)
end

function Base.isequal(a::LeftEndpoint, b::RightEndpoint)
    isequal(a.endpoint, b.endpoint) && a.included && b.included
end

function Base.isequal(a::RightEndpoint, b::LeftEndpoint)
    isequal(b, a)
end

function Base.isless(a::LeftEndpoint, b::LeftEndpoint)
    a.endpoint < b.endpoint || (a.endpoint == b.endpoint && a.included && !b.included)
end

function Base.isless(a::RightEndpoint, b::RightEndpoint)
    a.endpoint < b.endpoint || (a.endpoint == b.endpoint && !a.included && b.included)
end

function Base.isless(a::LeftEndpoint, b::RightEndpoint)
    a.endpoint < b.endpoint
end

function Base.isless(a::RightEndpoint, b::LeftEndpoint)
    a.endpoint < b.endpoint || (a.endpoint == b.endpoint && !(a.included && b.included))
end

# Comparisons between Scalars and Endpoints
Base.:(==)(a, b::Endpoint) = a == b.endpoint && b.included
Base.:(==)(a::Endpoint, b) = b == a

Base.isless(a, b::LeftEndpoint)  = a < b.endpoint || (a == b.endpoint && !b.included)
Base.isless(a, b::RightEndpoint) = a < b.endpoint
Base.isless(a::LeftEndpoint, b)  = a.endpoint < b
Base.isless(a::RightEndpoint, b) = a.endpoint < b || (a.endpoint == b && !a.included)


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "HiGHS"

version = v"0.3.1"

sources = [
    GitSource("https://github.com/ERGO-Code/HiGHS.git", "40fa96a10554d76b1413bf33cc79a2153aa1db13"),
    DirectorySource("./bundled"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
if [[ "${target}" == *86*-linux-musl* ]]; then
    pushd /opt/${target}/lib/gcc/${target}/*/include
    # Fix bug in Musl C library, see
    # https://github.com/JuliaPackaging/BinaryBuilder.jl/issues/387
    atomic_patch -p0 $WORKSPACE/srcdir/patches/mm_malloc.patch
    popd
fi
mkdir -p HiGHS/build
cd HiGHS/build
apk add --upgrade cmake
cmake -DCMAKE_INSTALL_PREFIX=$prefix \
    -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} \
    -DCMAKE_BUILD_TYPE=Release \
    -DBUILD_SHARED_LIBS=ON \
    -DFAST_BUILD=ON \
    -DJULIA=ON \
    -DIPX=ON ..
if [[ "${target}" == *-linux-* ]]; then
        make -j $nproc
else
    if [[ "${target}" == *-mingw* ]]; then
        cmake --build . --config Release
    else
        cmake --build . --config Release --parallel
    fi
fi
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = expand_cxxstring_abis(supported_platforms())

# Useful for testing to add the MIP and QP solvers.
# filter!(Sys.iswindows, platforms)

# The products that we will ensure are always built
products = [
    LibraryProduct("libhighs", :libhighs),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency(PackageSpec(name="CompilerSupportLibraries_jll", uuid="e66e0078-7015-5450-92f7-15fbd957f2ae"))
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; preferred_gcc_version = v"4.9")


text:
using Test

@testset "Core" begin
  include("Directed.jl")
  include("CPortGraphs.jl")
  include("Undirected.jl")
end

@testset "Monoidal" begin
  include("MonoidalDirected.jl")
  include("MonoidalUndirected.jl")
end

@testset "Algebras" begin
  include("Algebras.jl")
end

@testset "Algorithms" begin
  include("Algorithms.jl")
end

@testset "Expressions" begin
  include("Expressions.jl")
end

@testset "Scheduling" begin
  include("ScheduleUndirected.jl")
end

@testset "Serialization" begin
  include("GraphML.jl")
  include("JSON.jl")
end


text:
###########
# This gives an example of Gauss–Seidel
#
# It is _much_ more efficient to use BlockBandedMatrices.jl
###########

using BandedMatrices, LazyArrays, LinearAlgebra

function finitedifference_2d(n)
    h = 1/n
    D² = BandedMatrix(0 => Fill(-2,n), 1 => Fill(1,n-1), -1 => Fill(1,n-1))/h^2
    D_xx = kron(D², Eye(n))
    D_yy = kron(Eye(n), D²)
    D_xx + D_yy
end

function gaussseidel(L, U, b, x=copy(b), M=5)
    for _=1:M
        @view(x[1:end-1]) .= Mul(U , @view(x[2:end]))
        x[end] = 0
        x .= b .- x
        x .= Ldiv(L, x)
    end
    x
end

n = 400
Δt = (1/n^2)/4; Δ = finitedifference_2d(n); A = I - Δt*Δ  # 160k x 160k discretization

L = LowerTriangular(A)
U = UpperTriangular(@view A[1:end-1,2:end])

b = randn(size(A,1));
@time u = A\b; # 6s
x = copy(b)
@time gaussseidel(L,U, b, x, 20) # 1.6s
norm(x - u) # 6*10^(-10)


text:
using FEniCS
using PyPlot

L = 1; W = 0.2
mu = 1
rho = 1
delta = W/L
gamma = 0.4*delta^2
beta = 1.25
lambda_ = beta
g = gamma

# Create mesh and define function space
mesh = BoxMesh(Point([0, 0, 0]), Point([L, W, W]), 10, 3, 3)
V = VectorFunctionSpace(mesh, "P", 1)
c = Constant((0,0,0))
bc = DirichletBC(V, c, "on_boundary && x[1]<1E-14")
tol = 1E-14

#bc1 = fenics.DirichletBC(V.pyobject, c.pyobject, "on_boundary and x[1]<tol",tol=1E-14)
function epsilon(u)
     0.5*(nabla_grad(u)+Transpose(nabla_grad(u)))
end

function sigma(u)
     lambda_*nabla_div(u)*Identity(d) + 2*mu*epsilon(u)
end

u = TrialFunction(V)
d = geometric_dimension(u)# space dimension
v = TestFunction(V)
f = Constant((0, 0, -rho*g))
T = Constant((0, 0, 0))
a = inner(sigma(u), epsilon(v))*dx
L = dot(f, v)*dx + dot(T, v)*ds
u=FeFunction(V)
lvsolve(a,L,u,bc)
s = sigma(u) - (1.0/3)*tr(sigma(u))*Identity(d)
#fenics.plot(U.pyobject,mode="displacement")
von_Mises = sqrt(3.0/2*inner(s, s))
V = FunctionSpace(mesh, "P", 1)

von_Mises = project(von_Mises, V)
u_magnitude = sqrt(dot(u, u))
u_magnitude = project(u_magnitude, V)

# How the heck do we plot these things!!!
# FEniCS.Plot(von_Mises, title="Stress intensity")
# FEniCS.Plot(u_magnitude)


text:
AlwaysAssertionError() = AlwaysAssertionError("")

# The documentation for the `Base.@assert` macro says: "Warning: An assert might be
# disabled at various optimization levels."
# Therefore, we have the `always_assert` function. `always_assert` is like
# `Base.@assert`, except that `always_assert` will always run and will never be
# disabled.
function always_assert(cond::Bool)
    cond || throw(AlwaysAssertionError())
    return nothing
end


text:
#   Copyright (c) 2020-2022 Velexi Corporation
#
#   Licensed under the Apache License, Version 2.0 (the "License");
#   you may not use this file except in compliance with the License.
#   You may obtain a copy of the License at
#
#       http://www.apache.org/licenses/LICENSE-2.0
#
#   Unless required by applicable law or agreed to in writing, software
#   distributed under the License is distributed on an "AS IS" BASIS,
#   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
#   See the License for the specific language governing permissions and
#   limitations under the License.

"""
Unit tests for the Zero type.
"""

# --- Imports

# Standard library
import InteractiveUtils.subtypes
using Test

# GeometricAlgebra.jl
using GeometricAlgebra

# --- Constructor tests

@testset "Zero: constructor tests" begin
    # --- Inner constructor

    for precision_type in subtypes(AbstractFloat)
        B = Zero{precision_type}()
        @test B isa Zero{precision_type}
        @test B === Zero{precision_type}()
    end

    # --- Outer constructor

    B = Zero()
    @test B isa Zero{Float64}
end

@testset "Zero: zero()" begin
    # --- zero(M::AbstractMultivector)

    for precision_type in subtypes(AbstractFloat)
        # --- M::AbstractScalar

        B = zero(Zero{precision_type}())
        @test B isa Zero{precision_type}

        B = zero(One{precision_type}())
        @test B isa Zero{precision_type}

        B = zero(Scalar{precision_type}(5))
        @test B isa Zero{precision_type}

        # --- M::AbstractBlade

        B = zero(Blade{precision_type}([1 2 3]))
        @test B isa Zero{precision_type}

        B = zero(Pseudoscalar{precision_type}(10, 5))
        @test B isa Zero{precision_type}

        # --- M::AbstractMultivector

        B = zero(Multivector{precision_type}([Scalar(3)]))
        @test B isa Zero{precision_type}
    end

    # --- zero(::Type{<:AbstractMultivector})

    B = zero(AbstractScalar)
    @test B isa Zero{Float64}

    B = zero(Zero)
    @test B isa Zero{Float64}

    B = zero(One)
    @test B isa Zero{Float64}

    B = zero(Scalar)
    @test B isa Zero{Float64}

    B = zero(AbstractBlade)
    @test B isa Zero{Float64}

    B = zero(Blade)
    @test B isa Zero{Float64}

    B = zero(Pseudoscalar)
    @test B isa Zero{Float64}

    B = zero(AbstractMultivector)
    @test B isa Zero{Float64}

    B = zero(Multivector)
    @test B isa Zero{Float64}

    # --- zero(::Type{<:AbstractMultivector{T}})

    for precision_type in subtypes(AbstractFloat)
        # --- Type{<:AbstractScalar{T}}

        B = zero(AbstractScalar{precision_type})
        @test B isa Zero{precision_type}

        B = zero(Zero{precision_type})
        @test B isa Zero{precision_type}

        B = zero(One{precision_type})
        @test B isa Zero{precision_type}

        B = zero(Scalar{precision_type})
        @test B isa Zero{precision_type}

        # --- Type{<:AbstractBlade{T}}

        B = zero(AbstractBlade{precision_type})
        @test B isa Zero{precision_type}

        B = zero(Blade{precision_type})
        @test B isa Zero{precision_type}

        B = zero(Pseudoscalar{precision_type})
        @test B isa Zero{precision_type}

        # --- Type{<:AbstractMultivector{T}}

        B = zero(AbstractMultivector{precision_type})
        @test B isa Zero{precision_type}

        B = zero(Multivector{precision_type})
        @test B isa Zero{precision_type}
    end
end

# --- Test comparison operations

@testset "Zero: iszero(B)" begin
    # Basic functions
    for precision_type in subtypes(AbstractFloat)
        # Blade
        vectors = [3 3; 4 4; 0 1; 0 1]
        B = Blade{precision_type}(vectors)
        @test !iszero(B)

        vectors = [1 2; 1 2; 0 0]
        B = Blade{precision_type}(vectors)
        @test iszero(B)

        # Pseudoscalar
        B = Pseudoscalar{precision_type}(10, 3)
        @test !iszero(B)

        # Scalar
        B = Scalar{precision_type}(0)
        @test iszero(B)

        B = Scalar{precision_type}(3)
        @test !iszero(B)

        # Zero
        B = Zero{precision_type}()
        @test iszero(B)

        # One
        B = One{precision_type}()
        @test !iszero(B)
    end
end

# --- Test attribute methods

@testset "Zero: AbstractMultivector attribute functions" begin
    # Basic functions
    for precision_type in subtypes(AbstractFloat)
        B = Zero{precision_type}()
        @test dim(B) == 0
        @test grades(B) == [0]
        @test blades(B) == [B]
        @test norm(B) isa precision_type
        @test norm(B) == 0

        @test B[0] == [B]
        @test B[1] == []
    end
end

@testset "Zero: AbstractBlade attribute functions" begin
    # Basic functions
    for precision_type in subtypes(AbstractFloat)
        B = Zero{precision_type}()
        @test grade(B) == 0
        @test basis(B) == 1
        @test volume(B) isa precision_type
        @test volume(B) == 0
        @test sign(B) == 0
    end
end

@testset "Zero: AbstractScalar attribute functions" begin
    for precision_type in subtypes(AbstractFloat)
        B = Zero{precision_type}()
        @test value(B) isa precision_type
        @test value(B) == 0
    end
end

# --- Tests for AbstractMultivector interface functions

@testset "Zero: -(B)" begin
    for precision_type in subtypes(AbstractFloat)
        B = Zero{precision_type}()
        @test -B === B
    end
end

@testset "Zero: reverse(B)" begin
    for precision_type in subtypes(AbstractFloat)
        B = Zero{precision_type}()
        @test reverse(B) === B
    end
end

@testset "Zero: dual(B)" begin
    for precision_type in subtypes(AbstractFloat)
        B = Zero{precision_type}()
        expected_message = "The dual of Zero is not well-defined"

        test_dim = 5
        @test_throws ErrorException(expected_message) dual(B, test_dim)
        @test_throws ErrorException(expected_message) dual(B)
    end
end

@testset "Zero: convert(B)" begin
    for precision_type_converted in subtypes(AbstractFloat)
        for precision_type_src in subtypes(AbstractFloat)
            B = Zero{precision_type_src}()
            B_converted = convert(precision_type_converted, B)
            @test B_converted isa Zero{precision_type_converted}
        end
    end
end

# --- Tests for AbstractBlade interface functions

@testset "Zero: inverse(B)" begin
    for precision_type in subtypes(AbstractFloat)
        B = Zero{precision_type}()

        inverse_B = inv(B)
        @test inverse_B isa Scalar{precision_type}
        @test inverse_B == Scalar{precision_type}(Inf)
    end
end


text:
export buildD, findSEA

struct DistanceMatrix
    atoms::Vector{Any}
    distances::AbstractArray
end

struct SEA
    sets::Tuple
end

function distance(A::Atom, B::Atom)
    a = A.xyz
    b = B.xyz
    return (((a[1]-b[1])^2) + ((a[2]-b[2])^2) + ((a[3]-b[3])^2))^0.5
end

function buildD(V)
    len = size(V)[1]
    atoms = []
    arr = zeros(Float64, (len,len))
    for i = 1:len
        push!(atoms, V[i].Z)
    end
    for i = 1:len, j = i+1:len
        arr[i,j] = distance(V[i],V[j])
        arr[j,i] = arr[i,j]
    end
    return DistanceMatrix(atoms, arr)
end

function checkSEA(A::Vector{Float64}, B::Vector{Float64}, δ)::Bool
    a = sort(A)
    b = sort(B)
    z = a - b
    for i in z
        if i < 10.0^(-δ)
            continue
        else
            return false
        end
    end
    return true
end

function findSEA(D::DistanceMatrix, δ::Int)
    len = size(D.distances)[1]
    out = Tuple[]
    for i = 1:len, j = i+1:len
        if cmp(D.atoms[i], D.atoms[j]) == 0
            c = checkSEA(D.distances[:,i], D.distances[:,j], δ)
            if c
                push!(out, (i,j))
            end
        end
    end
    nonos = Int[]
    biggerun = Vector{Int}[]
    for i = 1:len
        if i in nonos
            continue
        else
            biggun = Int[]
            push!(biggun,i)
        end

        for k in out
            if i in k
                if i == k[1]
                    push!(biggun, k[2])
                    push!(nonos, k[2])
                else
                    push!(biggun, k[1])
                    push!(nonos, k[1])
                end
            end
        end
    push!(biggerun, biggun)
    end

    return biggerun
end


text:
# ******************************************************************************************
# Notices:
#
# Copyright © 2022 United States Government as represented by the Administrator of the
# National Aeronautics and Space Administration.  All Rights Reserved.
#
# Disclaimers
#
# No Warranty: THE SUBJECT SOFTWARE IS PROVIDED "AS IS" WITHOUT ANY WARRANTY OF ANY KIND,
# EITHER EXPRESSED, IMPLIED, OR STATUTORY, INCLUDING, BUT NOT LIMITED TO, ANY WARRANTY
# THAT THE SUBJECT SOFTWARE WILL CONFORM TO SPECIFICATIONS, ANY IMPLIED WARRANTIES OF
# MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE, OR FREEDOM FROM INFRINGEMENT, ANY
# WARRANTY THAT THE SUBJECT SOFTWARE WILL BE ERROR FREE, OR ANY WARRANTY THAT DOCUMENTATION,
# IF PROVIDED, WILL CONFORM TO THE SUBJECT SOFTWARE. THIS AGREEMENT DOES NOT, IN ANY MANNER,
# CONSTITUTE AN ENDORSEMENT BY GOVERNMENT AGENCY OR ANY PRIOR RECIPIENT OF ANY RESULTS,
# RESULTING DESIGNS, HARDWARE, SOFTWARE PRODUCTS OR ANY OTHER APPLICATIONS RESULTING FROM
# USE OF THE SUBJECT SOFTWARE.  FURTHER, GOVERNMENT AGENCY DISCLAIMS ALL WARRANTIES AND
# LIABILITIES REGARDING THIRD-PARTY SOFTWARE, IF PRESENT IN THE ORIGINAL SOFTWARE, AND
# DISTRIBUTES IT "AS IS."
#
# Waiver and Indemnity:  RECIPIENT AGREES TO WAIVE ANY AND ALL CLAIMS AGAINST THE UNITED
# STATES GOVERNMENT, ITS CONTRACTORS AND SUBCONTRACTORS, AS WELL AS ANY PRIOR RECIPIENT.
# IF RECIPIENT'S USE OF THE SUBJECT SOFTWARE RESULTS IN ANY LIABILITIES, DEMANDS, DAMAGES,
# EXPENSES OR LOSSES ARISING FROM SUCH USE, INCLUDING ANY DAMAGES FROM PRODUCTS BASED ON,
# OR RESULTING FROM, RECIPIENT'S USE OF THE SUBJECT SOFTWARE, RECIPIENT SHALL INDEMNIFY AND
# HOLD HARMLESS THE UNITED STATES GOVERNMENT, ITS CONTRACTORS AND SUBCONTRACTORS, AS WELL
# AS ANY PRIOR RECIPIENT, TO THE EXTENT PERMITTED BY LAW.  RECIPIENT'S SOLE REMEDY FOR ANY
# SUCH MATTER SHALL BE THE IMMEDIATE, UNILATERAL TERMINATION OF THIS AGREEMENT.
# ******************************************************************************************

"""
Replay buffer.
"""
mutable struct ReplayBuffer
    obs_buf::Matrix{Float32}        # stores observations
    act_buf::Matrix{Float32}        # stores actions
    rew_buf::Vector{Float32}        # stores rewards
    obs2_buf::Matrix{Float32}       # stores post-action observations
    done_buf::Vector{Float32}       # stores termination signal
    ptr::Int                        # current free index
    size::Int                       # current size
    max_size::Int                   # maximum size
end

function ReplayBuffer(obs_dim::Int, act_dim::Int, max_size::Int)
    obs_buf = zeros(Float32, obs_dim, max_size)
    act_buf = zeros(Float32, act_dim, max_size)
    rew_buf = zeros(Float32, max_size)
    obs2_buf = zeros(Float32, obs_dim, max_size)
    done_buf = zeros(Float32, max_size)
    ptr = 1
    size = 0
    return ReplayBuffer(obs_buf, act_buf, rew_buf, obs2_buf, done_buf, ptr, size, max_size)
end

"""
Store data from one step of simulation.
"""
function store!(
    buf::ReplayBuffer,
    obs::Vector{<:Real},
    act::Vector{<:Real},
    rew::Real,
    next_obs::Vector{<:Real},
    done::Bool
)
    buf.obs_buf[:,buf.ptr] = obs
    buf.act_buf[:,buf.ptr] = act
    buf.rew_buf[buf.ptr] = rew
    buf.obs2_buf[:,buf.ptr] = next_obs
    buf.done_buf[buf.ptr] = done
    buf.ptr = mod1(buf.ptr + 1, buf.max_size)
    buf.size = min(buf.size + 1, buf.max_size)
    return buf.size
end

"""
Randomly sample from replay buffer.
"""
function sample_batch(buf::ReplayBuffer, batch_size::Int=32)
    idxs = rand(1:buf.size, batch_size)
    batch = (
        obs=buf.obs_buf[:,idxs],
        act=buf.act_buf[:,idxs],
        rew=buf.rew_buf[idxs],
        obs2=buf.obs2_buf[:,idxs],
        done=buf.done_buf[idxs]
    )
    return batch
end

"""
Lazily clear replay buffer.
"""
function clear!(buf::ReplayBuffer)
    buf.ptr = 1
    buf.size = 0
    return nothing
end


text:
function reset!(alg::TrustRegion)
    alg.Δ[] = alg.Δ₀
    return nothing
end

function step!(x, alg::TrustRegion, f::Function, state)
    Δ = alg.Δ[]
    Δₘₐₓ = alg.Δₘₐₓ
    η = alg.η
    ηₛ = alg.ηₛ
    ηₑ = alg.ηₑ
    fx = state.f
    ∇fx = state.∇f
    B = approx_hessian(alg.hessian, state)
    model(p) = fx + ∇fx⋅p + (p'*B*p)/2
    p = alg.method(state, B, Δ)
    ρ = (fx - f(x + p)) / (model(zero(p)) - model(p))

    # Update radius
    if ρ ≥ ηₑ && norm(p) ≈ Δ # highly successful
        alg.Δ[] = min(alg.σₑ*Δ, Δₘₐₓ)
    elseif ρ < ηₛ
        alg.Δ[] = alg.σₛ*Δ
    end

    # Apply step
    if ρ ≥ η
        x .+= p
    end

    return nothing
end


text:
using StaticArrays

"""
    Diff = Union{Nothing,NTuple}

Type alias for Union{Nothing,NTuple}. This is used as type parameters for the fields of struct Step.
Mainly for convenient construction of any combination of ::Nothing, ::TTuple.
"""
const Diff_svec = Union{Nothing,SVector}
Diff_svec(::Nothing) = nothing
Diff_svec(t::Tuple) = SVector(t)

const Diff_tuple = Union{Nothing,Tuple}
Diff_tuple(::Nothing) = nothing

struct StepA{T<:Diff_svec, S<:Diff_svec, R<:Diff_svec, Q<:Diff_svec}
    drop_orbs :: T
    drop_kinks :: S
    add_orbs :: R
    add_kinks :: Q
end

StepA() = StepA(nothing, nothing, nothing, nothing)
"""
    StepA(a,b,c,d)

Construct a Step for any combination of ::Nothing and ::SVector.
"""
# Without using Diff_svec it would be needed to define 4! constructors for all combinations
StepA(a,b,c,d) = StepA(Diff_svec(a), Diff_svec(b), Diff_svec(c), Diff_svec(d))

struct StepB{T<:Diff_tuple, S<:Diff_tuple, R<:Diff_tuple, Q<:Diff_tuple}
    drop_orbs :: T
    drop_kinks :: S
    add_orbs :: R
    add_kinks :: Q
end

StepB() = StepB(nothing, nothing, nothing, nothing)
StepB(a,b,c,d) = StepB(Diff_tuple(a), Diff_tuple(b), Diff_tuple(c), Diff_tuple(d))

struct StepC
    drop_orbs :: Union{SVector, Nothing}
    drop_kinks :: Union{SVector, Nothing}
    add_orbs :: Union{SVector, Nothing}
    add_kinks :: Union{SVector, Nothing}
end

StepC() = StepC(nothing, nothing, nothing, nothing)

struct StepD{T<:SVector, S<:SVector, R<:SVector, Q<:SVector}
    drop_orbs :: T
    drop_kinks :: S
    add_orbs :: R
    add_kinks :: Q
end

emptyStepD(T::DataType) = StepD(SVector{0,T}(),SVector{0,T}(),SVector{0,T}(),SVector{0,T}())
StepD(a,b,c,d) = StepD(Diff_svec(a), Diff_svec(b), Diff_svec(c), Diff_svec(d))

struct StepE
    drop_orbs :: SVector
    drop_kinks :: SVector
    add_orbs :: SVector
    add_kinks :: SVector
end

emptyStepE(T::DataType) = StepE(SVector{0,T}(),SVector{0,T}(),SVector{0,T}(),SVector{0,T}())

struct StepF{T<:Tuple, S<:Tuple, R<:Tuple, Q<:Tuple}
    drop_orbs :: T
    drop_kinks :: S
    add_orbs :: R
    add_kinks :: Q
end

StepF() = StepF((),(),(),())
StepF(a,b,c,d) = StepF(Diff_tuple(a), Diff_tuple(b), Diff_tuple(c), Diff_tuple(d))

struct StepF{T<:Tuple, S<:Tuple, R<:Tuple, Q<:Tuple}
    drop_orbs :: T
    drop_kinks :: S
    add_orbs :: R
    add_kinks :: Q
end

StepF() = StepF((),(),(),())
StepF(a,b,c,d) = StepF(Diff_tuple(a), Diff_tuple(b), Diff_tuple(c), Diff_tuple(d))

struct StepG
    drop_orbs :: Tuple
    drop_kinks :: Tuple
    add_orbs :: Tuple
    add_kinks :: Tuple
end

StepG() = StepG((),(),(),())

drop_orbs!(occ, ::Nothing) = nothing
drop_orbs!(occ, s::SVector{0,T}) where T = nothing
drop_orbs!(occ, s::Tuple{}) where T = nothing
drop_orbs!(occ, s) = setdiff!(occ, s)

add_orbs!(occ, ::Nothing) = nothing
add_orbs!(occ, s::SVector{0,T}) where T = nothing
add_orbs!(occ, s::Tuple{}) where T = nothing
add_orbs!(occ, s) = union!(occ, s)

drop_kinks!(ck, ::Nothing) = nothing
drop_kinks!(ck, s::SVector{0,T}) where T = nothing
drop_kinks!(ck, s::Tuple{}) where T = nothing
drop_kinks!(ck, s) = setdiff!(ck , s)

add_kinks!(ck, ::Nothing) = nothing
add_orbs!(ck, s::SVector{0,T}) where T = nothing
add_orbs!(ck, s::Tuple{}) where T = nothing
function add_kinks!(ck, s)
    for k in s
        insert!(ck, searchsortedfirst(ck, by=first, first(k)), first(k) => last(k))
    end
end

function apply_step!(c, s)
    drop_orbs!(c.occupations, s.drop_orbs)
    drop_kinks!(c.kinks, s.drop_kinks)
    add_orbs!(c.occupations, s.add_orbs)
    add_kinks!(c.kinks, s.add_kinks)
end

import CPIMC: T2, T4, Kink, ImgTime, Configuration, Kinks, Orbital

struct Step{S<:Union{Nothing,<:Orbital,Configuration},T<:Union{Nothing,<:Orbital,Configuration}}
    drop :: S
    add :: T
end

Step() = Step(nothing,nothing)
Step(drop::Set{T}, add) where {T <: Orbital} = Step(Configuration(drop), add)
Step(drop, add::Set{T}) where {T <: Orbital} = Step(drop, Configuration(add))
Step(drop::Set{T}, add::Set{T}) where {T <: Orbital} = Step(Configuration(drop), Configuration(add))

function apply_step!(c::Configuration, step::Step)
    drop!(c, step.drop)
    add!(c, step.add)
    nothing
end

drop!(c::Configuration, n::Nothing) = nothing

drop!(c::Configuration{T}, o::T) where {T <: Orbital} = delete!(c.occupations, o)
function drop!(c::Configuration, oc)
    for o in oc
        drop!(c, o)
    end
end
drop!(c::Configuration, τ::ImgTime) = setdiff!(c.kinks, [τ => c.kinks[τ]])
drop!(c::Configuration, τ, k) = setdiff!(c.kinks, τ => k)
function drop!(c::Configuration{T}, ck::Kinks) where {T <: Orbital}
    for (τ, k) in ck
        drop!(c, τ, k)
    end
end
function drop!(c1::Configuration{T}, c2::Configuration{T}) where {T <: Orbital}
    drop!(c1, c2.occupations)
    drop!(c1, c2.kinks)
end
function add!(c::Configuration, n::Nothing)
    nothing
end

add!(c::Configuration{T}, o::T) where {T <: Orbital} = push!(c.occupations, o)
add!(c::Configuration, oc) = union!(c.occupations, oc)
add!(ck::Kinks, τ::ImgTime, k::Kink) = insert!(ck, searchsortedfirst(ck, by=first, τ), τ => k)
add!(c::Configuration{T}, p::Pair ) where {T <: Orbital} = add!(c.kinks, first(p), last(p))
function add!(ck::Kinks, ps::Pair{ImgTime, <:Kink{T}}...) where {T <: Orbital}
    for k in ps
        add!(ck, first(k), last(k))
    end
    return ck
end
add!(c::Configuration{T}, ck::Kinks) where {T <: Orbital} = add!(c.kinks, ck...)
function add!(c1::Configuration{T}, c2::Configuration{T}) where {T <: Orbital}
    add!(c1, c2.occupations)
    add!(c1, c2.kinks)
end

using CPIMC.PlaneWaves

S = sphere_with_same_spin(PlaneWave((0,0,0)),dk=1)
a = PlaneWave((-2,0,0))
b = PlaneWave((3,0,0))
c = PlaneWave((0,0,0))
d = PlaneWave((1,0,0))
e = PlaneWave((5,9,9))
f = PlaneWave((1,1,1))

g = PlaneWave(a.vec + e.vec - f.vec, Up)
h = PlaneWave(c.vec + d.vec - g.vec, Up)

i = PlaneWave((39,100,-5))
j = PlaneWave((100,0,6))

sd = Kinks( ImgTime(0.2) => T4(a,b,c,d),
       ImgTime(0.5) => T4(c,d,a,b),
       ImgTime(0.6) => T4(b,a,d,c),
       ImgTime(0.8) => T4(d,c,b,a) )

conf = Configuration(sphere(PlaneWave((0,0,0),Up),dk=4),sd)

droporbs = (c,d,f)
addorbs = (e,i,j)
dropkinks = (ImgTime(0.2) => T4(a,b,c,d), ImgTime(0.5) => T4(c,d,a,b))
addkinks = (ImgTime(0.9) => T2(a,b), ImgTime(0.1) => T4(a,a,a,a))

droporbs_set = Set((c,d,f))
addorbs_set = Set((e,i,j))
dropkinks_set = Set( (ImgTime(0.2) => T4(a,b,c,d), ImgTime(0.5) => T4(c,d,a,b)) )
addkinks_set = Set( (ImgTime(0.9) => T2(a,b), ImgTime(0.1) => T4(a,a,a,a)) )

using BenchmarkTools

println("\t Benchmark ")
println(repeat("=", 80))

println("\t construct empty Step:")
println(repeat(" -", 40), "\n")

print(typeof(StepA()), "\t:\t")
@btime StepA();
println("\n")
print(typeof(StepB()), "\t:\t")
@btime StepB();
println("\n")
print(typeof(StepC()), "\t:\t")
@btime StepC();
println("\n")
print(typeof(emptyStepD(Nothing)), "\t:\t")
@btime emptyStepD(Nothing);
println("\n")
print(typeof(emptyStepE(Nothing)), "\t:\t")
@btime emptyStepE(Nothing);
println("\n")
print(typeof(StepF()), "\t:\t")
@btime StepF();
println("\n")
print(typeof(StepG()), "\t:\t")
@btime StepG();
println("\n")
print("old version \t:\t")
@btime Step();
println("\n")

println("\t construct nonempty Step:\t")
print("drop ", length(droporbs), " orbs and ", length(dropkinks), " kinks \t")
println("add ", length(addorbs), " orbs and ", length(addkinks), " kinks")
println(repeat(" -", 40), "\n")
println("\n")

print(StepA, "\t:\t")
@btime StepA(droporbs,dropkinks,addorbs,addkinks);
println("\n")
print(StepB, "\t:\t")
@btime StepB(droporbs,dropkinks,addorbs,addkinks);
println("\n")
print(StepC, "\t:\t")
@btime StepC(droporbs,dropkinks,addorbs,addkinks);
println("\n")
print(StepD, "\t:\t")
@btime StepD(droporbs,dropkinks,addorbs,addkinks);
println("\n")
print(StepE, "\t:\t")
@btime StepE(droporbs,dropkinks,addorbs,addkinks);
println("\n")
print(StepF, "\t:\t")
@btime StepF(droporbs,dropkinks,addorbs,addkinks);
println("\n")
print(StepG, "\t:\t")
@btime StepG(droporbs,dropkinks,addorbs,addkinks);
println("\n")
print("old version \t:\t")
@btime Step(Configuration(Set(droporbs), dropkinks...), Configuration(Set(addorbs), addkinks...));
println("\n")

print("\t apply_step!:\t")
print("drop ", length(droporbs), " orbs and ", length(dropkinks), " kinks \t")
println("add ", length(addorbs), " orbs and ", length(addkinks), " kinks")
println(repeat(" -", 40), "\n")

x = StepA(droporbs,dropkinks,addorbs,addkinks);
print(typeof(x), "\t:\t")
@btime apply_step!(co,$x) setup=(co=$conf);
println("\n")
x = StepB(droporbs,dropkinks,addorbs,addkinks);
print(typeof(x), "\t:\t")
@btime apply_step!(co,$x) setup=(co=$conf);
println("\n")
x = StepC(droporbs,dropkinks,addorbs,addkinks);
print(typeof(x), "\t:\t")
@btime apply_step!(co,$x) setup=(co=$conf);
println("\n")
x = StepD(droporbs,dropkinks,addorbs,addkinks);
print(typeof(x), "\t:\t")
@btime apply_step!(co,$x) setup=(co=$conf);
println("\n")
x = StepE(droporbs,dropkinks,addorbs,addkinks);
print(typeof(x), "\t:\t")
@btime apply_step!(co,$x) setup=(co=$conf);
println("\n")
x = StepF(droporbs,dropkinks,addorbs,addkinks);
print(typeof(x), "\t:\t")
@btime apply_step!(co,$x) setup=(co=$conf);
println("\n")
x = StepG(droporbs,dropkinks,addorbs,addkinks);
print(typeof(x), "\t:\t")
@btime apply_step!(co,$x) setup=(co=$conf);
println("\n")
x = Step(Configuration(Set(droporbs), dropkinks...), Configuration(Set(addorbs), addkinks...));
print("old version \t:\t")
@btime apply_step!(co,$x) setup=(co=$conf);
println("\n")


text:
struct KernelRegression{R,T}
    x::R
    y::Vector{T}
    k::Function
    h::Float64
end

# Want a version for x as a vector (this) and X as a matrix (doesn't exist yet)
# function kernelregression{R <: Any,
#                           T <: Real}(x::Vector{R},
#                                      y::Vector{T};
#                                      kernel::Symbol = :epanechnikov,
#                                      bandwidth::Real = NaN,
#                                      getbandwidth::Function = kNN.bandwidth)
#     if isnan(bandwidth)
#         h_x = getbandwidth(x)
#         h_y = getbandwidth(y)
#         h = sqrt(h_x * h_y)
#     else
#         h = bandwidth
#     end
#     k = SmoothingKernels.kernels[kernel]
#     return KernelRegression(x, y, k, h)
# end

function StatsBase.predict(model::KernelRegression, x::Real)
    y, n = 0.0, length(model.x)
    h, k = model.h, model.k
    normalizer = 0.0
    for i in 1:n
        # TODO: Get definitions and usage of bandwidth right
        # w_i = (1 / h) * k((x - model.x[i]) / h)
        w_i = h * k(h * (x - model.x[i]))
        y += model.y[i] * w_i
        normalizer += w_i
    end
    return y / normalizer
end

function StatsBase.predict!(ys::Vector, model::KernelRegression, xs::AbstractVector{T}) where {T}
    n = length(xs)
    # @assert length(ys) == n
    for i in 1:n
        ys[i] = predict(model, xs[i])
    end
    return
end

function StatsBase.predict(model::KernelRegression, xs::AbstractVector{T}) where {T}
    ys = Array{T}(undef,length(xs))
    predict!(ys, model, xs)
    return ys
end


text:
using LazySets: block_to_dimension_indices, substitute_blocks, get_constrained_lowdimset
"""
    overapproximate(X::S, ::Type{S}, args...) where {S<:LazySet}

Overapproximating a set of type `S` with type `S` is a no-op.

### Input

- `X`       -- set
- `Type{S}` -- target set type
- `args`    -- further arguments (ignored)

### Output

The input set.
"""
function overapproximate(X::S, ::Type{S}, args...) where {S<:LazySet}
    return X
end

"""
    overapproximate(S::LazySet{N},
                    ::Type{<:HPolygon},
                    [ε]::Real=Inf)::HPolygon where {N<:Real}

Return an approximation of a given 2D convex set.
If no error tolerance is given, or is `Inf`, the result is a box-shaped polygon.
Otherwise the result is an ε-close approximation as a polygon.

### Input

- `S`        -- convex set, assumed to be two-dimensional
- `HPolygon` -- type for dispatch
- `ε`        -- (optional, default: `Inf`) error bound

### Output

A polygon in constraint representation.
"""
function overapproximate(S::LazySet{N},
                         ::Type{<:HPolygon},
                         ε::Real=Inf
                        )::HPolygon where {N<:Real}
    @assert dim(S) == 2
    if ε == Inf
        constraints = Vector{LinearConstraint{N}}(undef, 4)
        constraints[1] = LinearConstraint(DIR_EAST(N), ρ(DIR_EAST(N), S))
        constraints[2] = LinearConstraint(DIR_NORTH(N), ρ(DIR_NORTH(N), S))
        constraints[3] = LinearConstraint(DIR_WEST(N), ρ(DIR_WEST(N), S))
        constraints[4] = LinearConstraint(DIR_SOUTH(N), ρ(DIR_SOUTH(N), S))
        return HPolygon(constraints, sort_constraints=false)
    else
        return tohrep(approximate(S, ε))
    end
end

"""
    overapproximate(S::LazySet, ε::Real)::HPolygon

Alias for `overapproximate(S, HPolygon, ε)`.
"""
function overapproximate(S::LazySet,
                         ε::Real;
                        )::HPolygon
    return overapproximate(S, HPolygon, ε)
end

"""
    overapproximate(S::LazySet,
                    Type{<:Hyperrectangle})::Union{Hyperrectangle, EmptySet}

Return an approximation of a given set as a hyperrectangle.

### Input

- `S`              -- set
- `Hyperrectangle` -- type for dispatch

### Output

A hyperrectangle.
"""
function overapproximate(S::LazySet,
                         ::Type{<:Hyperrectangle};
                        )::Union{Hyperrectangle, EmptySet}
    return box_approximation(S)
end

"""
    overapproximate(S::CartesianProductArray{N, <:AbstractHyperrectangle{N}},
                    ::Type{<:Hyperrectangle}) where {N<:Real}

Return a tight overapproximation of the cartesian product array of a finite
number of convex sets with and hyperrectangle.

### Input

- `S`              -- cartesian product array of a finite number of convex set
- `Hyperrectangle` -- type for dispatch

### Output

A hyperrectangle.

### Algorithm

This method falls back to the corresponding `convert` method. Since the sets wrapped
by the cartesian product array are hyperrectangles, it can be done efficiently
without overapproximation.
"""
function overapproximate(S::CartesianProductArray{N, <:AbstractHyperrectangle{N}},
                          ::Type{<:Hyperrectangle}) where {N<:Real}
    return convert(Hyperrectangle, S)
end

"""
    overapproximate(S::CartesianProduct{N, <:AbstractHyperrectangle{N}, <:AbstractHyperrectangle{N}},
                    ::Type{<:Hyperrectangle}) where {N<:Real}

Return a tight overapproximation of the cartesian product of two
hyperrectangles by a new hyperrectangle.

### Input

- `S`              -- cartesian product of two hyperrectangular sets
- `Hyperrectangle` -- type for dispatch

### Output

A hyperrectangle.

### Algorithm

This method falls back to the corresponding `convert` method. Since the sets wrapped
by the cartesian product are hyperrectangles, it can be done efficiently
without overapproximation.
"""
function overapproximate(S::CartesianProduct{N, <:AbstractHyperrectangle{N}, <:AbstractHyperrectangle{N}},
                          ::Type{<:Hyperrectangle}) where {N<:Real}
    return convert(Hyperrectangle, S)
end

"""
    overapproximate(lm::LinearMap{N, <:AbstractHyperrectangle{N}},
                    ::Type{Hyperrectangle}) where {N}

Return a tight overapproximation of the linear map of a hyperrectangular set
using a hyperrectangle.

### Input

- `S`              -- linear map of a hyperrectangular set
- `Hyperrectangle` -- type for dispatch

### Output

A hyperrectangle.

### Algorithm

If `c` and `r` denote the center and vector radius of a hyperrectangle `H`,
a tight hyperrectangular overapproximation of `M * H` is obtained by transforming
`c ↦ M*c` and `r ↦ abs.(M) * c`, where `abs.(⋅)` denotes the element-wise absolute
value operator.
"""
function overapproximate(lm::LinearMap{N, <:AbstractHyperrectangle{N}},
                         ::Type{Hyperrectangle}) where {N<:Real}
    M, X = lm.M, lm.X
    center_MX = M * center(X)
    radius_MX = abs.(M) * radius_hyperrectangle(X)
    return Hyperrectangle(center_MX, radius_MX)
end

"""
    overapproximate(S::LazySet)::Union{Hyperrectangle, EmptySet}

Alias for `overapproximate(S, Hyperrectangle)`.
"""
overapproximate(S::LazySet;
               )::Union{Hyperrectangle, EmptySet} =
    overapproximate(S, Hyperrectangle)

"""
    overapproximate(S::ConvexHull{N, Zonotope{N}, Zonotope{N}},
                    ::Type{<:Zonotope})::Zonotope where {N<:Real}

Overapproximate the convex hull of two zonotopes.

### Input

- `S`        -- convex hull of two zonotopes
- `Zonotope` -- type for dispatch

### Algorithm

This function implements the method proposed in [1].
The convex hull of two zonotopes ``Z₁`` and ``Z₂`` of the same order,
that we write

```math
Z_j = ⟨c^{(j)}, g^{(j)}_1, …, g^{(j)}_p⟩
```
for ``j = 1, 2``, can be overapproximated as follows:

```math
CH(Z_1, Z_2) ⊆ \\frac{1}{2}⟨c^{(1)}+c^{(2)}, g^{(1)}_1+g^{(2)}_1, …, g^{(1)}_p+g^{(2)}_p, c^{(1)}-c^{(2)}, g^{(1)}_1-g^{(2)}_1, …, g^{(1)}_p-g^{(2)}_p⟩.
```

If the zonotope order is not the same, this algorithm calls
`reduce_order` to reduce the order to the minimum of the arguments.

It should be noted that the output zonotope is not necessarily the minimal
enclosing zonotope, which is in general expensive in high dimensions. This is
further investigated in [2].

[1] Reachability of Uncertain Linear Systems Using Zonotopes, A. Girard.
    HSCC 2005.

[2] Zonotopes as bounding volumes, L. J. Guibas et al, Proc. of Symposium on
    Discrete Algorithms, pp. 803-812.
"""
function overapproximate(S::ConvexHull{N, Zonotope{N}, Zonotope{N}},
                         ::Type{<:Zonotope})::Zonotope where {N<:Real}
    Z1, Z2 = S.X, S.Y

    # reduce to the same order if possible
    if order(Z1) != order(Z2)
        min_order = min(order(Z1), order(Z2))
        Z1 = reduce_order(Z1, min_order)
        Z2 = reduce_order(Z2, min_order)
    end

    if order(Z1) >= order(Z2)
        c, G = _overapproximate_convex_hull_zonotope(Z1, Z2)
    else
        c, G = _overapproximate_convex_hull_zonotope(Z2, Z1)
    end
    return Zonotope(c, G)
end

# assumes that dim(Z1) == dim(Z2) and order(Z1) >= order(Z2)
function _overapproximate_convex_hull_zonotope(Z1::Zonotope{N}, Z2::Zonotope{N}) where {N}
    c = (Z1.center + Z2.center)/N(2)

    # the case of equal order is treated separately to avoid a slicing (this creates a copy)
    if order(Z1) == order(Z2)
        G = hcat(Z1.generators .+ Z2.generators,
                 Z1.center - Z2.center,
                 Z1.generators .- Z2.generators)/N(2)
    else
        G = hcat((Z1.generators[:, 1:ngens(Z2)] .+ Z2.generators)/N(2),
                 (Z1.center - Z2.center)/N(2),
                 (Z1.generators[:, 1:ngens(Z2)] .- Z2.generators)/N(2),
                 Z1.generators[:, ngens(Z2)+1:end])
    end
    return c, G
end

"""
    overapproximate(Z::Zonotope, ::Type{<:Hyperrectangle})::Hyperrectangle

Return a tight overapproximation of a zonotope with an axis-aligned box.

### Input

- `Z`              -- zonotope
- `Hyperrectangle` -- type for dispatch

### Output

A hyperrectangle.

### Algorithm

This function implements the method in [Section 5.1.2, 1]. A zonotope
``Z = ⟨c, G⟩`` can be overapproximated tightly by an axis-aligned box
(i.e. a `Hyperrectangle`) such that its center is ``c`` and the radius along
dimension ``i`` is the column-sum of the absolute values of the ``i``-th row
of ``G`` for ``i = 1,…, p``, where ``p`` is the number of generators of ``Z``.

[1] *Althoff, M., Stursberg, O., & Buss, M. (2010). Computing reachable sets of
hybrid systems using a combination of zonotopes and polytopes. Nonlinear analysis:
hybrid systems, 4(2), 233-249.*
"""
function overapproximate(Z::Zonotope, ::Type{<:Hyperrectangle})::Hyperrectangle
    r = sum(abs.(Z.generators), dims=2)[:]
    return Hyperrectangle(Z.center, r)
end

"""
    overapproximate(X::LazySet{N}, dir::AbstractDirections{N})::HPolytope{N}
        where {N}

Overapproximating a set with template directions.

### Input

- `X`   -- set
- `dir` -- (concrete) direction representation

### Output

An `HPolytope` overapproximating the set `X` with the directions from `dir`.
"""
function overapproximate(X::LazySet{N},
                         dir::AbstractDirections{N}
                        )::HPolytope{N} where {N}
    halfspaces = Vector{LinearConstraint{N}}()
    sizehint!(halfspaces, length(dir))
    H = HPolytope(halfspaces)
    for d in dir
        addconstraint!(H, LinearConstraint(d, ρ(d, X)))
    end
    return H
end

"""
    overapproximate(X::LazySet{N},
                    dir::Type{<:AbstractDirections})::HPolytope{N} where {N}

Overapproximating a set with template directions.

### Input

- `X`   -- set
- `dir` -- type of direction representation

### Output

A `HPolytope` overapproximating the set `X` with the directions from `dir`.
"""
function overapproximate(X::LazySet{N},
                         dir::Type{<:AbstractDirections}
                        )::HPolytope{N} where {N}
    return overapproximate(X, dir{N}(dim(X)))
end

"""
    overapproximate(S::LazySet{N}, ::Type{<:Interval}) where {N<:Real}

Return the overapproximation of a real unidimensional set with an interval.

### Input

- `S`        -- one-dimensional set
- `Interval` -- type for dispatch

### Output

An interval.

### Algorithm

The method relies on the exact conversion to `Interval`. Two support
function evaluations are needed in general.
"""
function overapproximate(S::LazySet{N}, ::Type{<:Interval}) where {N<:Real}
    @assert dim(S) == 1 "cannot overapproximate a $(dim(S))-dimensional set with an `Interval`"
    return convert(Interval, S)
end
function overapproximate_cap_helper(X::LazySet{N},             # compact set
                                    P::AbstractPolyhedron{N},  # polyhedron
                                    dir::AbstractDirections{N};
                                    kwargs...
                                   ) where {N<:Real}
    Hi = constraints_list(P)
    m = length(Hi)
    Q = HPolytope{N}()

    for di in dir
        ρ_X_Hi_min = ρ(di, X ∩ Hi[1], kwargs...)
        for i in 2:m
            ρ_X_Hi = ρ(di, X ∩ Hi[i], kwargs...)
            if ρ_X_Hi < ρ_X_Hi_min
                ρ_X_Hi_min = ρ_X_Hi
            end
        end
        addconstraint!(Q, HalfSpace(di, ρ_X_Hi_min))
    end
    return Q
end

"""
    overapproximate(cap::Intersection{N, <:LazySet, <:AbstractPolyhedron{N}},
                    dir::AbstractDirections{N};
                    kwargs...
                   ) where {N<:Real}

Return the overapproximation of the intersection between a compact set and a
polytope given a set of template directions.

### Input

- `cap`         -- intersection of a compact set and a polytope
- `dir`         -- template directions
- `kwargs`      -- additional arguments that are passed to the support function
                   algorithm

### Output

A polytope in H-representation such that the normal direction of each half-space
is given by an element of `dir`.

### Algorithm

Let `di` be a direction drawn from the set of template directions `dir`.
Let `X` be the compact set and let `P` be the polytope. We overapproximate the
set `X ∩ H` with a polytope in constraint representation using a given set of
template directions `dir`.

The idea is to solve the univariate optimization problem `ρ(di, X ∩ Hi)` for
each half-space in the set `P` and then take the minimum.
This gives an overapproximation of the exact support function.

This algorithm is inspired from [G. Frehse, R. Ray. Flowpipe-Guard Intersection
for Reachability Computations with Support
Functions](https://www.sciencedirect.com/science/article/pii/S1474667015371809).

### Notes

This method relies on having available the `constraints_list` of the polytope
`P`.

This method of overapproximations can return a non-empty set even if the original
intersection is empty.
"""
function overapproximate(cap::Intersection{N,
                                           <:LazySet,
                                           <:AbstractPolyhedron{N}},
                         dir::AbstractDirections{N};
                         kwargs...
                        ) where {N<:Real}
    return overapproximate_cap_helper(cap.X, cap.Y, dir; kwargs...)
end

# symmetric method
function overapproximate(cap::Intersection{N,
                                           <:AbstractPolyhedron{N},
                                           <:LazySet},
                         dir::AbstractDirections{N};
                         kwargs...
                        ) where {N<:Real}
    return overapproximate_cap_helper(cap.Y, cap.X, dir; kwargs...)
end

# disambiguation
function overapproximate(cap::Intersection{N,
                                           <:AbstractPolytope{N},
                                           <:AbstractPolyhedron{N}},
                         dir::AbstractDirections{N};
                         kwargs...
                        ) where {N<:Real}
    return overapproximate_cap_helper(cap.X, cap.Y, dir; kwargs...)
end

# symmetric method
function overapproximate(cap::Intersection{N,
                                           <:AbstractPolyhedron{N},
                                           <:AbstractPolytope{N}},
                         dir::AbstractDirections{N};
                         kwargs...
                        ) where {N<:Real}
    return overapproximate_cap_helper(cap.Y, cap.X, dir; kwargs...)
end

# disambiguation
function overapproximate(cap::Intersection{N,
                                           <:AbstractPolytope{N},
                                           <:AbstractPolytope{N}},
                         dir::AbstractDirections{N};
                         kwargs...
                        ) where {N<:Real}
    return overapproximate_cap_helper(cap.X, cap.Y, dir; kwargs...)
end

"""
    overapproximate(cap::Intersection{N, <:HalfSpace{N}, <:AbstractPolytope{N}},
                    dir::AbstractDirections{N};
                    [kwargs]...
                   ) where {N<:Real}

Return the overapproximation of the intersection between a half-space and a
polytope given a set of template directions.

### Input

- `cap`         -- intersection of a half-space and a polytope
- `dir`         -- template directions
- `kwargs`      -- additional arguments that are passed to the support function
                   algorithm

### Output

A polytope in H-representation such that the normal direction of each half-space
is given by an element of `dir`.
"""
function overapproximate(cap::Intersection{N,
                                           <:HalfSpace{N},
                                           <:AbstractPolytope{N}},
                         dir::AbstractDirections{N};
                         kwargs...
                        ) where {N<:Real}
    H = HPolytope{N}()
    c = constraints_list(H)
    append!(c, constraints_list(cap.X))
    append!(c, constraints_list(cap.Y))
    return overapproximate(H, dir; kwargs...)
end

# symmetric method
function overapproximate(cap::Intersection{N,
                                           <:AbstractPolytope{N},
                                           <:HalfSpace{N}},
                         dir::AbstractDirections{N};
                         kwargs...
                        ) where {N<:Real}
    return overapproximate(swap(cap), dir; kwargs...)
end

# ==========================================
# Functionality that requires TaylorModels
# ==========================================

# function to be loaded by Requires
function load_taylormodels_overapproximation()

return quote

using .TaylorModels: Taylor1, TaylorN, TaylorModelN, TaylorModel1,
                     polynomial, remainder, domain,
                     normalize_taylor, linear_polynomial,
                     constant_term, evaluate, mid

# helper functions
@inline get_linear_coeffs(p::Taylor1) = linear_polynomial(p).coeffs[2:end]
@inline get_linear_coeffs(p::TaylorN) = linear_polynomial(p).coeffs[2].coeffs

"""
    overapproximate(vTM::Vector{TaylorModel1{T, S}},
                    ::Type{Zonotope}) where {T, S}

Overapproximate a taylor model in one variable with a zonotope.

### Input

- `vTM`      -- `TaylorModel1`
- `Zonotope` --  type for dispatch

### Output

A zonotope that overapproximates the range of the given taylor model.

### Examples

If the polynomials are linear, this functions exactly transforms to a zonotope.
However, the nonlinear case necessarily introduces overapproximation error.
Consider the linear case first:

```julia
julia> using LazySets, TaylorModels

julia> const IA = IntervalArithmetic;

julia> I = IA.Interval(-0.5, 0.5) # interval remainder
[-0.5, 0.5]

julia> x₀ = IA.Interval(0.0) # expansion point
[0, 0]

julia> D = IA.Interval(-3.0, 1.0)
[-3, 1]

julia> p1 = Taylor1([2.0, 1.0], 2) # define a linear polynomial
 2.0 + 1.0 t + 𝒪(t³)

julia> p2 = Taylor1([0.9, 3.0], 2) # define another linear polynomial
 0.9 + 3.0 t + 𝒪(t³)

julia> vTM = [TaylorModel1(pi, I, x₀, D) for pi in [p1, p2]]
2-element Array{TaylorModel1{Float64,Float64},1}:
 2.0 + 1.0 t + [-0.5, 0.5]
 0.9 + 3.0 t + [-0.5, 0.5]
```

Here, `vTM` is a taylor model vector, since each component is a taylor model in
one variable (`TaylorModel1`). Using `overapproximate(vTM, Zonotope)` we can
compute its associated zonotope in generator representation:

```julia
julia> using LazySets.Approximations

julia> Z = overapproximate(vTM, Zonotope);

julia> center(Z)
2-element Array{Float64,1}:
  1.0
 -2.1

julia> Matrix(genmat(Z))
2×3 Array{Float64,2}:
 2.0  0.5  0.0
 6.0  0.0  0.5
```

Note how the generators of this zonotope mainly consist of two pieces: one comes
from the linear part of the polynomials, and another one that corresponds to the
interval remainder. This conversion gives the same upper and lower bounds as the
range evaluation using interval arithmetic:

```julia
julia> X = box_approximation(Z)
Hyperrectangle{Float64}([1.0, -2.1], [2.5, 6.5])

julia> Y = evaluate(vTM[1], vTM[1].dom) × evaluate(vTM[2], vTM[2].dom)
[-1.5, 3.5] × [-8.60001, 4.40001]

julia> H = convert(Hyperrectangle, Y) # this IntevalBox is the same as X
Hyperrectangle{Float64}([1.0, -2.1], [2.5, 6.5])
```
However, the zonotope returns better results if we want to approximate the `TM`,
since it is not axis-aligned:

```julia
julia> d = [-0.35, 0.93];

julia> ρ(d, Z) < ρ(d, X)
true
```

This function also works if the polynomials are non-linear; for example suppose
that we add a third polynomial with a quadratic term:

```julia
julia> p3 = Taylor1([0.9, 3.0, 1.0], 3);

julia> vTM = [TaylorModel1(pi, I, x₀, D) for pi in [p1, p2, p3]]
3-element Array{TaylorModel1{Float64,Float64},1}:
           2.0 + 1.0 t + [-0.5, 0.5]
           0.9 + 3.0 t + [-0.5, 0.5]
  0.9 + 3.0 t + 1.0 t² + [-0.5, 0.5]

julia> Z = overapproximate(vTM, Zonotope);

julia> center(Z)
3-element Array{Float64,1}:
  1.0
 -2.1
  0.8999999999999999

julia> Matrix(genmat(Z))
3×4 Array{Float64,2}:
 2.0  0.5  0.0  0.0
 6.0  0.0  0.5  0.0
 6.0  0.0  0.0  6.5
```
The fourth and last generator corresponds to the addition between the interval
remainder and the box overapproximation of the nonlinear part of `p3` over the
domain.

### Algorithm

Let ``\\text{vTM} = (p, I)`` be a vector of ``m`` taylor models, where ``I``
is the interval remainder in ``\\mathbb{R}^m``. Let ``p_{lin}``
(resp. ``p_{nonlin}``) correspond to the linear (resp. nonlinear) part of each
scalar polynomial.

The range of ``\\text{vTM}`` can be enclosed by a zonotope with center ``c``
and matrix of generators ``G``, ``Z = ⟨c, G⟩``, by performing a conservative
linearization of ``\\text{vTM}``:

```math
    vTM' = (p', I') := (p_{lin} − p_{nonlin} , I + \\text{Int}(p_{nonlin})).
```

This algorithm proceeds in two steps:

1- Conservatively linearize ``\\text{vTM}`` as above and compute a box
   overapproximation of the nonlinear part.
2- Transform the linear taylor model to a zonotope exactly through variable
   normalization onto the symmetric intervals ``[-1, 1]``.
"""
function overapproximate(vTM::Vector{TaylorModel1{T, S}},
                            ::Type{Zonotope}) where {T, S}
    m = length(vTM)

    # preallocations
    c = Vector{T}(undef, m) # center of the zonotope
    gen_lin = Matrix{T}(undef, m, 1) # generator of the linear part
    gen_rem = Vector{T}(undef, m) # generators of the remainder

    # compute overapproximation
    return _overapproximate_vTM_zonotope!(vTM, c, gen_lin, gen_rem)
end

"""
    overapproximate(vTM::Vector{TaylorModelN{N, T, S}},
                    ::Type{Zonotope}) where {N,T, S}

Overapproximate a multivariate taylor model with a zonotope.

### Input

- `vTM`      -- `TaylorModelN`
- `Zonotope` -- type for dispatch

### Output

A zonotope that overapproximates the range of the given taylor model.

### Examples

Consider a vector of two 2-dimensional taylor models of order 2 and 4
respectively.

```julia
julia> using LazySets, LazySets.Approximations, TaylorModels

julia> const IA = IntervalArithmetic;

julia> x₁, x₂ = set_variables(Float64, ["x₁", "x₂"], order=8)
2-element Array{TaylorN{Float64},1}:
  1.0 x₁ + 𝒪(‖x‖⁹)
  1.0 x₂ + 𝒪(‖x‖⁹)

julia> x₀ = IntervalBox(0..0, 2) # expansion point
[0, 0] × [0, 0]

julia> Dx₁ = IA.Interval(0.0, 3.0) # domain for x₁
[0, 3]

julia> Dx₂ = IA.Interval(-1.0, 1.0) # domain for x₂
[-1, 1]

julia> D = Dx₁ × Dx₂ # take the cartesian product of the domain on each variable
[0, 3] × [-1, 1]

julia> r = IA.Interval(-0.5, 0.5) # interval remainder
[-0.5, 0.5]

julia> p1 = 1 + x₁^2 - x₂
 1.0 - 1.0 x₂ + 1.0 x₁² + 𝒪(‖x‖⁹)

julia> p2 = x₂^3 + 3x₁^4 + x₁ + 1
 1.0 + 1.0 x₁ + 1.0 x₂³ + 3.0 x₁⁴ + 𝒪(‖x‖⁹)

julia> vTM = [TaylorModelN(pi, r, x₀, D) for pi in [p1, p2]]
2-element Array{TaylorModelN{2,Float64,Float64},1}:
             1.0 - 1.0 x₂ + 1.0 x₁² + [-0.5, 0.5]
   1.0 + 1.0 x₁ + 1.0 x₂³ + 3.0 x₁⁴ + [-0.5, 0.5]

julia> Z = overapproximate(vTM, Zonotope);

julia> center(Z)
2-element Array{Float64,1}:
   5.5
 124.0

julia> Matrix(genmat(Z))
2×4 Array{Float64,2}:
 0.0  -1.0  5.0    0.0
 1.5   0.0  0.0  123.0
```

### Algorithm

We refer to the algorithm description for the univariate case.
"""
function overapproximate(vTM::Vector{TaylorModelN{N, T, S}},
                         ::Type{Zonotope}) where {N,T, S}
    m = length(vTM)
    n = N # number of variables is get_numvars() in TaylorSeries

    # preallocations
    c = Vector{T}(undef, m) # center of the zonotope
    gen_lin = Matrix{T}(undef, m, n) # generator of the linear part
    gen_rem = Vector{T}(undef, m) # generators for the remainder

    # compute overapproximation
    return _overapproximate_vTM_zonotope!(vTM, c, gen_lin, gen_rem)
end

function _overapproximate_vTM_zonotope!(vTM, c, gen_lin, gen_rem)
    @inbounds for (i, x) in enumerate(vTM)
        xpol, xdom = polynomial(x), domain(x)

        # linearize the TM
        pol_lin = constant_term(xpol) + linear_polynomial(xpol)
        rem_nonlin = remainder(x)

        # build an overapproximation of the nonlinear terms
        pol_nonlin = xpol - pol_lin
        rem_nonlin += evaluate(pol_nonlin, xdom)

        # normalize the linear polynomial to the symmetric interval [-1, 1]
        Q = normalize_taylor(pol_lin, xdom, true)

        # build the generators
        α = mid(rem_nonlin)
        c[i] = constant_term(Q) + α  # constant terms
        gen_lin[i, :] = get_linear_coeffs(Q) # linear terms
        gen_rem[i] = abs(rem_nonlin.hi - α)
    end
    return Zonotope(c, hcat(gen_lin, Diagonal(gen_rem)))
end

end # quote
end # load_taylormodels_overapproximation

# ==========================================
# Lazy linear maps of cartesian products
# ==========================================

"""
    overapproximate(lm::LinearMap{N, <:CartesianProductArray{N}},
                    ::Type{CartesianProductArray{N, S}}
                   ) where {N, S<:LazySet{N}}

Decompose a lazy linear map of a cartesian product array while keeping the
original block structure.

### Input

- `lm`                    -- lazy linear map of cartesian product array
- `CartesianProductArray` -- type for dispatch

### Output

A `CartesianProductArray` representing the decomposed linear map.
"""
function overapproximate(lm::LinearMap{N, <:CartesianProductArray{N}},
                         ::Type{CartesianProductArray{N, S}}
                        ) where {N, S<:LazySet{N}}
    cpa = array(lm.X)
    arr = Vector{S}(undef, length(cpa))
    return _overapproximate_lm_cpa!(arr, lm.M, cpa, S)
end

"""
    overapproximate(lm::LinearMap{N, <:CartesianProductArray{N}},
                    ::Type{<:CartesianProductArray},
                    dir::Type{<:AbstractDirections}) where {N}

Decompose a lazy linear map of a cartesian product array with template
directions while keeping the original block structure.

### Input

- `lm`                    -- lazy linear map of a cartesian product array
- `CartesianProductArray` -- type for dispatch
- `dir`                   -- template directions for overapproximation

### Output

A `CartesianProductArray` representing the decomposed linear map.
"""
function overapproximate(lm::LinearMap{N, <:CartesianProductArray{N}},
                         ::Type{<:CartesianProductArray},
                         dir::Type{<:AbstractDirections}) where {N}
    cpa = array(lm.X)
    arr = Vector{HPolytope{N}}(undef, length(cpa))
    return _overapproximate_lm_cpa!(arr, lm.M, cpa, dir)
end

"""
    overapproximate(lm::LinearMap{N, <:CartesianProductArray{N}},
                    ::Type{<:CartesianProductArray},
                    set_type::Type{<:LazySet}) where {N}

Decompose a lazy linear map of a cartesian product array with a given set type
while keeping the original block structure.

### Input

- `lm`                    -- lazy linear map of a cartesian product array
- `CartesianProductArray` -- type for dispatch
- `set_type`              -- set type for overapproximation

### Output

A `CartesianProductArray` representing the decomposed linear map.
"""
function overapproximate(lm::LinearMap{N, <:CartesianProductArray{N}},
                         ::Type{<:CartesianProductArray},
                         set_type::Type{<:LazySet}) where {N}
    cpa = array(lm.X)
    arr = Vector{set_type{N}}(undef, length(cpa))
    return _overapproximate_lm_cpa!(arr, lm.M, cpa, set_type)
end

function _overapproximate_lm_cpa!(arr, M, cpa, overapprox_option)
    # construct Minkowski sum for block row
    function _block_row(cpa::Vector{S}, M::AbstractMatrix{N},
                        row_range::UnitRange{Int}) where {N, S<:LazySet{N}}
        arr_inner = Vector{LinearMap{N, <:S}}(undef, length(cpa))
        col_start_ind, col_end_ind = 1, 0
        @inbounds for (j, bj) in enumerate(cpa)
            col_end_ind += dim(bj)
            arr_inner[j] =
                LinearMap(M[row_range, col_start_ind:col_end_ind], bj)
            col_start_ind = col_end_ind + 1
        end
        return MinkowskiSumArray(arr_inner)
    end

    row_start_ind, row_end_ind = 1, 0
    @inbounds for (i, bi) in enumerate(cpa)
        row_end_ind += dim(bi)
        ms = _block_row(cpa, M, row_start_ind:row_end_ind)
        arr[i] = overapproximate(ms, overapprox_option)
        row_start_ind = row_end_ind + 1
    end

    return CartesianProductArray(arr)
end

"""
    overapproximate(cap::Intersection{N,
                            <:CartesianProductArray{N},
                            <:AbstractPolyhedron{N}},
                       ::Type{CartesianProductArray}, oa) where {N}

Return the intersection of the cartesian product of a finite number of convex sets and a polyhedron.

### Input

 - `cap` -- Lazy intersection of cartesian product array and polyhedron
 - `CartesianProductArray` -- type for dispatch
 - `oa`  -- template directions

### Output

The intersection between `cpa` and `P` with overapproximation in each constrained block.

### Algorithm

The intersection is only needed to be taken in the elements of the cartesian product array (subsets of
variables, or "blocks") which are constrained in `P`.
Hence we first search for constrained blocks and then take the intersection of
a lower-dimensional Cartesian product of these blocks with the projection of `Y`
onto the variables of these blocks. (This projection is syntactic and exact.)
The result is a `CartesianProductArray` with the same block structure as in `X`.
However, when we decompose back our set we overapproximate during projection operation,
therefore we need to specify overapproximation strategy (it is Hyperrectangle by default)
"""
function overapproximate(cap::Intersection{N,
                                            <:CartesianProductArray{N},
                                            <:AbstractPolyhedron{N}},
                            ::Type{CartesianProductArray}, oa) where {N}

    cpa, P = cap.X, cap.Y

    cpa_low_dim, vars, block_structure, blocks = get_constrained_lowdimset(cpa, P)

    low_intersection = intersection(cpa_low_dim, project(P, vars))

    if isempty(low_intersection)
        return EmptySet{N}()
    end

    decomposed_low_set = decompose(low_intersection, block_structure, oa)

    return substitute_blocks(decomposed_low_set, cpa, blocks)
end

# symmetric method
function overapproximate(cap::Intersection{N,
                                            <:AbstractPolyhedron{N},
                                            <:CartesianProductArray{N}},
                            ::Type{CartesianProductArray}, oa) where {N}
    overapproximate(Intersection(cap.Y, cap.X), oa)
end


text:
using LinearAlgebra

export Generic2DLattice

abstract type Abstract2DLattice end

"""

Structure containing the real and the reciprical basis vectors of a surface
"""
struct Generic2DLattice <: Abstract2DLattice
    a1
    a2
    b1
    b2
end

"""
# Arguments
* a1 first lattice vector
* a2 second lattice vector

# Example
` Generic2DLattice([1,0],[1,0]) `
"""
function Generic2DLattice(a1, a2)
    area = a2[2]*a1[1]-a2[1]*a1[2]
    b1 = 2π * [ a2[2],-a2[1]] / area
    b2 = 2π * [-a1[2], a1[1]] / area
    Generic2DLattice(a1,a2,b1,b2)
end

struct HexGrid <: Abstract2DLattice
    a
    l::Generic2DLattice
end

HexGrid(a) = HexGrid(a,Generic2DLattice([a,0],[cos(pi/180*120),sin(pi/180*120)].*a))

function rotate_reciprocal(lattice::Generic2DLattice,direction)
	unitx = [lattice.b1 lattice.b2] * direction
	ang = -atan(unitx[2],unitx[1])
	R = [cos(ang) -sin(ang)
	     sin(ang) cos(ang)]
	R*lattice.b1,R*lattice.b2
end


text:
const PACKAGE_ROOT = dirname(dirname(@__FILE__))
const TESTROOT = joinpath(PACKAGE_ROOT, "test")
const DOCROOT = joinpath(PACKAGE_ROOT, "docs")
const DOCSOURCE = joinpath(DOCROOT, "src")
const EXAMPLESROOT = joinpath(PACKAGE_ROOT, "examples")
const EXAMPLES = [
    # ("Linear regression", "linear-regression"),
    ]


text:
"""
Archimedes' Approximation of Pi

## References

- https://itech.fgcu.edu/faculty/clindsey/mhf4404/archimedes/archimedes.html
- https://betterexplained.com/articles/prehistoric-calculus-discovering-pi/
"""
module ArchimedesPi

#=

"If an angle of a triangle be bisected and the straight line cutting the angle cut the base also,
the segments of the base will have the same ratio as the remaining sides of the triangle."

-- Euclid's Elements

Let $r = 1$, $OSide_{n} = a$, $ISide_{n} = b$, $OSide_{2n} = x$, $ISide_{2n} = y$, we have

$$
(x / 2) / (a / 2 - x / 2) = 1 / (a / b)  =>
x = a * b / (a + b)
$$

and

$$
(y / 2) / (x / 2) = (b / 2) / y  =>
2 (y ^ 2) = x * b
$$

So the perimeter of circumscribing polygon of 2n sides is (the harmonic mean)

$$
O_{2n} = 2n * OSide_{2n} = 2n / (1 / a + 1 / b) = \frac{2}{1 / I_{n} + 1 / O_{n}}
$$

and the perimeter of inscribing polygon of 2n sides is (the geometric mean)

$$
I_{2n} = 2n * ISide_{2n}
       = 2n * \sqrt{x * b} / \sqrt{2}
       = n * \sqrt{2} * \sqrt{x * b}
       = \sqrt{2} * \sqrt{(OSide_{2n} * n) * (ISide_{n} * n)}
       = \sqrt{2} * \sqrt{(O_{2n} / 2) * I_{n}}
       = \sqrt{I_{n} * O_{2n}}
$$
=#

function calculate_pi(n_digits::Int=100)::BigFloat
    prec = Int(ceil(log2(10) * (n_digits + 1)))  # +1 because pi = 3.(...)
    setprecision(prec * 2) do
        # Starts from a hexagon
        outer = big"2.0" * sqrt(big"3.0")
        inner = big"3.0"
        df = outer - inner

        while true
            outer = big"2.0" / (one(BigFloat) / inner + one(BigFloat) / outer)
            inner = sqrt(inner * outer)

            df = begin
                new_df = outer - inner
                if new_df >= df
                    break
                end
                new_df
            end
        end

        pi_archimedes = (outer + inner) / big"2.0"
        # Still uses the default rounding, intentionally
        setprecision(prec)
        pi_archimedes + big(0)
    end
end

#=
Verifies the answer:

    pi_archimedes = calculate_pi()
    setprecision(precision(pi_archimedes))  # Warning: Set the global BigFloat precision!
    println(pi_archimedes - big(pi))
=#

end # module

max_stars_repo_name: gyk/TrivialSolutions

max_issues_repo_name: gyk/TrivialSolutions

max_forks_repo_name: gyk/TrivialSolutions


text:
export expandsampledata

"""
    expandsampledata(x, experiment[, colnames])

Return a DataFrame containing the sample data for all or some of the column names in the chosen `experiment`.
Columns are the same as those in `sampledata(x)`.

If `colnames` is supplied, each row of the returned `DataFrame` corresponds to an entry of `colnames` and contains the data for the sample matching that column in the specified experiment.

If `colnames` is not supplied, each row of the returned `DataFrame` corresponds to a column of the specified experiment.

An error is raised if the requested columns do not have a matching sample in `samplemap(x)`.
Use [`dropunused`](@ref) to remove unused columns from each experiment prior to calling this function.

A warning is raised if `sampledata(x)` contains duplicate sample names.
In such cases, data is taken from the first entry for each sample.

A warning is raised if `samplemap(x)` contains multiple occurrences of the same experiment/colname combination with a different sample.
In such cases, the first occurrence of the combination is used.

# Examples
```jldoctest
julia> using MultiAssayExperiments;

julia> x = exampleobject();

julia> expandsampledata(x, "foo")
10×2 DataFrame
 Row │ name      disease
     │ String    String
─────┼───────────────────
   1 │ Patient1  good
   2 │ Patient1  good
   3 │ Patient1  good
   4 │ Patient2  bad
   5 │ Patient2  bad
   6 │ Patient2  bad
   7 │ Patient3  good
   8 │ Patient3  good
   9 │ Patient3  good
  10 │ Patient4  bad

julia> expandsampledata(x, "foo", ["foo2", "foo1"])
2×2 DataFrame
 Row │ name      disease
     │ String    String
─────┼───────────────────
   1 │ Patient1  good
   2 │ Patient1  good
```
"""
function expandsampledata(x::MultiAssayExperiment, exp::String)
    names = SummarizedExperiments.coldata(experiment(x, exp))[!,"name"]
    return expandsampledata(x, exp, names)
end

function expandsampledata(x::MultiAssayExperiment, exp::String, colnames::Vector{String})
    sm = samplemap(x)
    collated = Dict{String,String}()
    for r in eachrow(sm)
        if r.experiment == exp
            if haskey(collated, r.colname)
                if collated[r.colname] != r.sample
                    @warn "multiple rows in 'samplemap(x)' for '" * sample * "' in '" * exp * "', using the first only"
                end
            else
                collated[r.colname] = r.sample
            end
        end
    end

    samples = Dict{String,Int}()
    sd = sampledata(x)
    sn = sd[!,"name"]
    for s in 1:length(sn)
        sample = sn[s]
        if haskey(samples, sample)
            @warn "multiple rows in 'sampledata(x)' for '" * sample * "', using the first only"
            continue
        end
        samples[sample] = s
    end

    indices = Vector{Int}(undef, length(colnames))
    for c in 1:length(colnames)
        colname = colnames[c]
        if !haskey(collated, colname)
            throw(ErrorException("failed to find column name '" * colname * "' in the sample mapping for experiment '" * exp * "'"))
        end

        x = collated[colname];
        if !haskey(samples, x)
            throw(ErrorException("failed to find sample name '" * x * "' in the sample data"))
        end

        indices[c] = samples[x]
    end

    return sd[indices,:]
end


text:
using EquationsSolver
using Documenter

DocMeta.setdocmeta!(EquationsSolver, :DocTestSetup, :(using EquationsSolver); recursive=true)

makedocs(;
    modules=[EquationsSolver],
    authors="yjy <522432938@qq.com> and contributors",
    #repo="https://github.com/jake484/EquationsSolver.jl/blob/{commit}{path}#{line}",
    sitename="EquationsSolver.jl",
    format=Documenter.HTML(;
        # prettyurls=get(ENV, "CI", "false") == "true",
        # canonical="https://jake484.github.io/EquationsSolver.jl",
        # assets=["assets/logo.png"],
    ),
    pages=[
        "Home" => "index.md",
        "Tutorials" => Any[
            "tutorials/LinearProblem.md",
            "tutorials/NLProblem.md",
        ]
    ],
)

deploydocs(;
    repo="github.com/jake484/EquationsSolver.jl",
    devbranch="main",
)