text:
#= concrete implementations of binary intersections between sets =#

export intersection

"""
    intersection(S::AbstractSingleton{N},
                 X::LazySet{N}
                ) where {N<:Real}

Return the intersection of a singleton with another set.

### Input

- `S` -- singleton
- `X` -- another set

### Output

If the sets intersect, the result is `S`.
Otherwise, the result is the empty set.
"""
function intersection(S::AbstractSingleton{N},
                      X::LazySet{N}
                     ) where {N<:Real}
    return element(S) ∈ X ? S : EmptySet{N}(dim(S))
end

# symmetric method
function intersection(X::LazySet{N},
                      S::AbstractSingleton{N}
                     ) where {N<:Real}
    return intersection(S, X)
end

# disambiguation
function intersection(S1::AbstractSingleton{N},
                      S2::AbstractSingleton{N}
                     ) where {N<:Real}
    return element(S1) == element(S2) ? S1 : EmptySet{N}(dim(S1))
end

"""
    intersection(L1::Line{N}, L2::Line{N}
                ) where {N<:Real}

Return the intersection of two 2D lines.

### Input

- `L1` -- first line
- `L2` -- second line

### Output

If the lines are identical, the result is the first line.
If the lines are parallel and not identical, the result is the empty set.
Otherwise the result is the only intersection point.

### Examples

The line ``y = -x + 1`` intersected with the line ``y = x``:

```jldoctest
julia> intersection(Line([-1., 1.], 0.), Line([1., 1.], 1.))
Singleton{Float64,Array{Float64,1}}([0.5, 0.5])

julia> intersection(Line([1., 1.], 1.), Line([1., 1.], 1.))
Line{Float64,Array{Float64,1}}([1.0, 1.0], 1.0)
```
"""
function intersection(L1::Line{N}, L2::Line{N}
                     ) where {N<:Real}
    b = [L1.b, L2.b]
    a = [transpose(L1.a); transpose(L2.a)]
    try
        # results in LAPACKException or SingularException if parallel
        return Singleton(a \ b)
    catch e
        @assert e isa LAPACKException || e isa SingularException "unexpected " *
            "$(typeof(e)) from LAPACK occurred while intersecting lines:\n$e"
        # lines are parallel
        if an_element(L1) ∈ L2
            # lines are identical
            return L1
        else
            # lines are parallel but not identical
            return EmptySet{N}(dim(L1))
        end
    end
end

"""
    intersection(H1::AbstractHyperrectangle{N},
                 H2::AbstractHyperrectangle{N}
                ) where {N<:Real}

Return the intersection of two hyperrectangles.

### Input

- `H1` -- first hyperrectangle
- `H2` -- second hyperrectangle

### Output

If the hyperrectangles do not intersect, the result is the empty set.
Otherwise the result is the hyperrectangle that describes the intersection.

### Algorithm

In each isolated direction `i` we compute the rightmost left border and the
leftmost right border of the hyperrectangles.
If these borders contradict, then the intersection is empty.
Otherwise the result uses these borders in each dimension.
"""
function intersection(H1::AbstractHyperrectangle{N},
                      H2::AbstractHyperrectangle{N}
                     ) where {N<:Real}
    n = dim(H1)
    c1 = center(H1)
    c2 = center(H2)
    r1 = radius_hyperrectangle(H1)
    r2 = radius_hyperrectangle(H2)
    high = Vector{N}(undef, n)
    low = Vector{N}(undef, n)
    for i in 1:n
        high1 = c1[i] + r1[i]
        low1 = c1[i] - r1[i]
        high2 = c2[i] + r2[i]
        low2 = c2[i] - r2[i]
        high[i] = min(high1, high2)
        low[i] = max(low1, low2)
        if high[i] < low[i]
            return EmptySet{N}(dim(H1))
        end
    end
    return Hyperrectangle(high=high, low=low)
end

# disambiguation
function intersection(S::AbstractSingleton{N},
                      H::AbstractHyperrectangle{N}
                     ) where {N<:Real}
    return invoke(intersection, Tuple{typeof(S), LazySet{N}}, S, H)
end
function intersection(H::AbstractHyperrectangle{N},
                      S::AbstractSingleton{N}
                     ) where {N<:Real}
    return invoke(intersection, Tuple{typeof(S), LazySet{N}}, S, H)
end

"""
    intersection(x::Interval{N}, y::Interval{N}
                ) where {N<:Real}

Return the intersection of two intervals.

### Input

- `x` -- first interval
- `y` -- second interval

### Output

If the intervals do not intersect, the result is the empty set.
Otherwise the result is the interval that describes the intersection.
"""
function intersection(x::Interval{N}, y::Interval{N}
                     ) where {N<:Real}
    if min(y) > max(x) || min(x) > max(y)
        return EmptySet{N}(1)
    else
        return Interval(max(min(x), min(y)), min(max(x), max(y)))
    end
end

"""
    intersection(X::Interval{N}, hs::HalfSpace{N}
                ) where {N<:Real}

Compute the intersection of an interval and a half-space.

### Input

- `X`  -- interval
- `hs` -- half-space

### Output

If the sets do not intersect, the result is the empty set.
If the interval is fully contained in the half-space, the result is the original
interval.
Otherwise the result is the interval that describes the intersection.

### Algorithm

We first handle the special case that the normal vector `a` of `hs` is close to
zero.
Then we distinguish the cases that `hs` is a lower or an upper bound.
"""
function intersection(X::Interval{N}, hs::HalfSpace{N}
                     ) where {N<:Real}
    @assert dim(hs) == 1 "cannot take the intersection between an interval " *
                         "and a $(dim(hs))-dimensional half-space"

    a = hs.a[1]
    b = hs.b
    if _isapprox(a, zero(N))
        if _geq(b, zero(N))
            # half-space is universal
            return X
        else
            # half-space is empty
            return EmptySet{N}(1)
        end
    end

    # idea:
    # ax ≤ b  <=>  x ≤ b/a  (if a > 0)
    # ax ≤ b  <=>  x ≥ b/a  (if a < 0)
    b_over_a = b / a
    lo = min(X)
    hi = max(X)

    if a > zero(N)
        # half-space is an upper bound
        # check whether ax ≤ b for x = lo
        if _leq(lo, b_over_a)
            # new upper bound: min(hi, b_over_a)
            if _leq(hi, b_over_a)
                return X
            else
                return Interval(lo, b_over_a)
            end
        else
            return EmptySet{N}(1)
        end
    else
        # half-space is a lower bound
        # check whether ax ≤ b for x = hi
        if _geq(hi, b_over_a)
            # new lower bound: max(lo, b_over_a)
            if _leq(b_over_a, lo)
                return X
            else
                return Interval(b_over_a, hi)
            end
        else
            return EmptySet{N}(1)
        end
    end
end

# symmetric method
function intersection(hs::HalfSpace{N}, X::Interval{N}
                     ) where {N<:Real}
    return intersection(X, hs)
end

"""
    intersection(X::Interval{N}, hp::Hyperplane{N}
                ) where {N<:Real}

Compute the intersection of an interval and a hyperplane.

### Input

- `X`  -- interval
- `hp` -- hyperplane

### Output

If the sets do not intersect, the result is the empty set.
Otherwise the result is the singleton that describes the intersection.
"""
function intersection(X::Interval{N}, hp::Hyperplane{N}
                     ) where {N<:Real}
    @assert dim(hp) == 1 "cannot take the intersection between an interval " *
                         "and a $(dim(hp))-dimensional hyperplane"

    # a one-dimensional hyperplane is just a point
    p = hp.b / hp.a[1]
    if _leq(min(X), p) && _leq(p, max(X))
        return Singleton([p])
    else
        return EmptySet{N}(1)
    end
end

# symmetric method
function intersection(hp::Hyperplane{N}, X::Interval{N}
                     ) where {N<:Real}
    return intersection(X, hp)
end

"""
    intersection(X::Interval{N}, Y::LazySet{N}
                ) where {N<:Real}

Compute the intersection of an interval and a convex set.

### Input

- `X` -- interval
- `Y` -- convex set

### Output

If the sets do not intersect, the result is the empty set.
Otherwise the result is the interval that describes the intersection, which may
be of type `Singleton` if the intersection is very small.
"""
function intersection(X::Interval{N}, Y::LazySet{N}
                     ) where {N<:Real}
    lower = max(min(X), -ρ(N[-1], Y))
    upper = min(max(X), ρ(N[1], Y))
    if _isapprox(lower, upper)
        return Singleton([lower])
    elseif lower < upper
        return Interval(lower, upper)
    else
        return EmptySet{N}(1)
    end
end

# symmetric method
function intersection(Y::LazySet{N}, X::Interval{N}
                     ) where {N<:Real}
    return intersection(X, Y)
end

# disambiguation
function intersection(X::Interval{N}, H::AbstractHyperrectangle{N}
                     ) where {N<:Real}
    return invoke(intersection, Tuple{Interval{N}, LazySet{N}}, X, H)
end
function intersection(H::AbstractHyperrectangle{N}, X::Interval{N}
                     ) where {N<:Real}
    return invoke(intersection, Tuple{Interval{N}, LazySet{N}}, X, H)
end
function intersection(X::Interval{N}, S::AbstractSingleton{N}
                     ) where {N<:Real}
    return invoke(intersection, Tuple{typeof(S), LazySet{N}}, S, X)
end
function intersection(S::AbstractSingleton{N}, X::Interval{N}
                     ) where {N<:Real}
    return invoke(intersection, Tuple{typeof(S), LazySet{N}}, S, X)
end

"""
    intersection(P1::AbstractHPolygon{N},
                 P2::AbstractHPolygon{N},
                 [prune]::Bool=true
                ) where {N<:Real}

Return the intersection of two polygons in constraint representation.

### Input

- `P1`    -- first polygon
- `P2`    -- second polygon
- `prune` -- (optional, default: `true`) flag for removing redundant constraints

### Output

If the polygons do not intersect, the result is the empty set.
Otherwise the result is the polygon that describes the intersection.

### Algorithm

We just combine the constraints of both polygons.
To obtain a linear-time algorithm, we interleave the constraints.
If there are two constraints with the same normal vector, we choose the tighter
one.

Redundancy of constraints is checked with
[`remove_redundant_constraints!(::AbstractHPolygon)`](@ref).
"""
function intersection(P1::AbstractHPolygon{N},
                      P2::AbstractHPolygon{N},
                      prune::Bool=true
                     ) where {N<:Real}
    # all constraints of one polygon are processed; now add the other polygon's
    # constraints
    @inline function add_remaining_constraints!(c, i, c1, i1, duplicates)
        c[i+1:length(c)-duplicates] = c1[i1:length(c1)]
        return true
    end

    # choose the constraint of c1; the directions are different
    @inline function choose_first_diff_dir!(c, i, i1, i2, c1, c2, duplicates)
        c[i] = c1[i1]
        if i1 == length(c1)
            return add_remaining_constraints!(c, i, c2, i2, duplicates)
        end
        return false
    end

    # choose the constraint of c1; the directions are equivalent (i.e., linearly
    # dependent)
    @inline function choose_first_same_dir!(c, i, i1, i2, c1, c2, duplicates)
        c[i] = c1[i1]
        if i1 == length(c1)
            if i2 < length(c2)
                return add_remaining_constraints!(c, i, c2, i2+1, duplicates)
            end
            return true
        elseif i2 == length(c2)
            return add_remaining_constraints!(c, i, c1, i1+1, duplicates)
        end
        return false
    end

    c1 = constraints_list(P1)
    c2 = constraints_list(P2)
    if length(c1) == 0
        return P2
    elseif length(c2) == 0
        return P1
    end
    # TODO: use common vector type of P1 and P2, see #2046
    c = Vector{LinearConstraint{N, Vector{N}}}(undef, length(c1) + length(c2))
    i1 = 1
    i2 = 1
    duplicates = 0
    for i in 1:length(c)
        if c1[i1].a <= c2[i2].a
            if c2[i2].a <= c1[i1].a
                duplicates += 1
                # constraints have the same normal vector: take the tighter one
                if is_tighter_same_dir_2D(c1[i1], c2[i2])
                    # first constraint is tighter
                    if choose_first_same_dir!(c, i, i1, i2, c1, c2, duplicates)
                        break
                    end
                else
                    # second constraint is tighter
                    if choose_first_same_dir!(c, i, i2, i1, c2, c1, duplicates)
                        break
                    end
                end
                i1 += 1
                i2 += 1
            else
                # first constraint comes first
                if choose_first_diff_dir!(c, i, i1, i2, c1, c2, duplicates)
                    break
                end
                i1 += 1
            end
        else
            # second constraint comes first
            if choose_first_diff_dir!(c, i, i2, i1, c2, c1, duplicates)
                break
            end
            i2 += 1
        end
    end
    if duplicates > 0
        deleteat!(c, length(c)-duplicates+1:length(c))
    end

    P = HPolygon(c, sort_constraints=false)
    if prune
        remove_redundant_constraints!(P)
        if isempty(P)
            return EmptySet{N}(2)
        end
    end
    return P
end

using MathProgBase.SolverInterface: AbstractMathProgSolver

"""
    intersection(P1::AbstractPolyhedron{N},
                 P2::AbstractPolyhedron{N};
                 backend=default_lp_solver(N)) where {N<:Real}

Compute the intersection of two polyhedra.

### Input

- `P1`      -- polyhedron
- `P2`      -- polyhedron
- `backend` -- (optional, default: `default_lp_solver(N)`) the LP solver used
               for the removal of redundant constraints; see the `Notes` section
               below for details

### Output

An `HPolyhedron` resulting from the intersection of `P1` and `P2`, with the
redundant constraints removed, or an empty set if the intersection is empty.
If one of the arguments is a polytope, the result is an `HPolytope` instead.

### Notes

The default value of the solver backend is `default_lp_solver(N)` and it is used
to run a feasiblity LP to remove the redundant constraints of the intersection.

If you want to use the `Polyhedra` library, pass an appropriate backend. For
example, to use the default Polyhedra library use
`default_polyhedra_backend(P, N)` or use `CDDLib.Library()` for the CDD library.

There are some shortcomings of the removal of constraints using the default
Polyhedra library; see e.g. #1038 and Polyhedra#146. It is safer to check for
emptiness of intersection before calling this function in those cases.

### Algorithm

This implementation unifies the constraints of the two sets obtained from the
`constraints_list` method.
"""
function intersection(P1::AbstractPolyhedron{N},
                      P2::AbstractPolyhedron{N};
                      backend=default_lp_solver(N)) where {N<:Real}

    # if one of P1 or P2 is bounded => the result is bounded
    HPOLY = (P1 isa AbstractPolytope || P2 isa AbstractPolytope) ?
        HPolytope : HPolyhedron

    # concatenate the linear constraints
    clist_P1 = _normal_Vector(P1) # TODO fix to similar type
    clist_P2 = _normal_Vector(P2)
    Q = HPOLY([clist_P1; clist_P2])

    # remove redundant constraints
    if backend isa AbstractMathProgSolver
        # if Q is empty => the feasiblity LP for the list of constraints of Q
        # is infeasible and remove_redundant_constraints! returns `false`
        if remove_redundant_constraints!(Q, backend=backend)
            return Q
        else
            return EmptySet{N}(dim(P1))
        end
    else
        # the correct way for this condition would be to check if `backend`
        # isa Polyhedra.Library; since that would require using Polyhedra: Library
        # and it is an optional dependency we opt to fallback without checking

        # convert to a Polyhedra's hrep
        Qph = polyhedron(Q; backend=backend)

        # remove the redundancies
        removehredundancy!(Qph)

        if isempty(Qph)
            return EmptySet{N}(dim(P1))
        else
            # convert back to HPOLY
            return convert(HPOLY, Qph)
        end
    end
end

# disambiguation
function intersection(S::AbstractSingleton{N},
                      P::AbstractPolyhedron{N}
                     ) where {N<:Real}
    return invoke(intersection, Tuple{typeof(S), LazySet{N}}, S, P)
end
function intersection(P::AbstractPolyhedron{N},
                      S::AbstractSingleton{N}
                     ) where {N<:Real}
    return invoke(intersection, Tuple{typeof(S), LazySet{N}}, S, P)
end
function intersection(X::Interval{N}, P::AbstractPolyhedron{N}
                     ) where {N<:Real}
    return invoke(intersection, Tuple{Interval{N}, LazySet{N}}, X, P)
end
function intersection(P::AbstractPolyhedron{N}, X::Interval{N}
                     ) where {N<:Real}
    return invoke(intersection, Tuple{Interval{N}, LazySet{N}}, X, P)
end

"""
    intersection(P1::Union{VPolytope{N}, VPolygon{N}},
                 P2::Union{VPolytope{N}, VPolygon{N}};
                 [backend]=default_polyhedra_backend(P1, N),
                 [prunefunc]=removevredundancy!) where {N<:Real}

Compute the intersection of two polytopes in vertex representation.

### Input

- `P1`        -- polytope in vertex representation
- `P2`        -- polytope in vertex representation
- `backend`   -- (optional, default: `default_polyhedra_backend(P1, N)`) the
                 backend for polyhedral computations
- `prunefunc` -- (optional, default: `removevredundancy!`) function to prune
                 the vertices of the result

### Output

A `VPolygon` if both arguments are `VPolygon`s, and a `VPolytope` otherwise.
"""
function intersection(P1::Union{VPolytope{N}, VPolygon{N}},
                      P2::Union{VPolytope{N}, VPolygon{N}};
                      backend=default_polyhedra_backend(P1, N),
                      prunefunc=removevredundancy!) where {N<:Real}
    Q1 = polyhedron(convert(VPolytope, P1); backend=backend)
    Q2 = polyhedron(convert(VPolytope, P2); backend=backend)
    Pint = Polyhedra.intersect(Q1, Q2)
    prunefunc(Pint)
    res = VPolytope(Pint)
    if P1 isa VPolygon && P2 isa VPolygon
        return convert(VPolygon, res)
    end
    return res
end

"""
    intersection(cup::UnionSet{N}, X::LazySet{N}) where {N<:Real}

Return the intersection of a union of two convex sets and another convex set.

### Input

- `cup` -- union of two convex sets
- `X`   -- convex set

### Output

The union of the pairwise intersections, expressed as a `UnionSet`.
If one of those sets is empty, only the other set is returned.
"""
function intersection(cup::UnionSet{N}, X::LazySet{N}) where {N<:Real}
    return intersection(cup.X, X) ∪ intersection(cup.Y, X)
end

# symmetric method
function intersection(X::LazySet{N}, cup::UnionSet{N}) where {N<:Real}
    return intersection(cup, X)
end

# disambiguation
function intersection(cup::UnionSet{N}, S::AbstractSingleton{N}) where {N<:Real}
    return element(S) ∈ cup ? S : EmptySet{N}(dim(S))
end
function intersection(S::AbstractSingleton{N}, cup::UnionSet{N}) where {N<:Real}
    return invoke(intersection, Tuple{UnionSet{N}, typeof(S)}, cup, S)
end

"""
    intersection(cup::UnionSetArray{N}, X::LazySet{N}) where {N<:Real}

Return the intersection of a union of a finite number of convex sets and another
convex set.

### Input

- `cup` -- union of a finite number of convex sets
- `X`   -- convex set

### Output

The union of the pairwise intersections, expressed as a `UnionSetArray`.
"""
function intersection(cup::UnionSetArray{N}, X::LazySet{N}) where {N<:Real}
    return UnionSetArray([intersection(Y, X) for Y in array(cup)])
end

# symmetric method
function intersection(X::LazySet{N}, cup::UnionSetArray{N}) where {N<:Real}
    return intersection(cup, X)
end

# disambiguation
function intersection(cup::UnionSetArray{N},
                      S::AbstractSingleton{N}) where {N<:Real}
    return element(S) ∈ cup ? S : EmptySet{N}(dim(S))
end
function intersection(S::AbstractSingleton{N},
                      cup::UnionSetArray{N}) where {N<:Real}
    return invoke(intersection, Tuple{UnionSetArray{N}, typeof(S)}, cup, S)
end

"""
    intersection(L::LinearMap{N}, S::LazySet{N}) where {N<:Real}

Return the intersection of a lazy linear map and a convex set.

### Input

- `L` -- linear map
- `S` -- convex set

### Output

The polytope obtained by the intersection of `l.M * L.X` and `S`.
"""
function intersection(L::LinearMap{N}, S::LazySet{N}) where {N<:Real}
    return intersection(linear_map(L.M, L.X), S)
end

# symmetric method
function intersection(S::LazySet{N}, L::LinearMap{N}) where {N}
    return intersection(L, S)
end

# disambiguation
function intersection(L1::LinearMap{N}, L2::LinearMap{N}) where {N}
    return intersection(linear_map(L1.M, L1.X), linear_map(L2.M, L2.X))
end

"""
    intersection(U::Universe{N}, X::LazySet{N}) where {N<:Real}

Return the intersection of a universe and a convex set.

### Input

- `U` -- universe
- `X` -- convex set

### Output

The set `X`.
"""
function intersection(U::Universe{N}, X::LazySet{N}) where {N<:Real}
    return X
end

# symmetric method
function intersection(X::LazySet{N}, U::Universe{N}) where {N<:Real}
    return X
end

# disambiguation
function intersection(U::Universe{N}, ::Universe{N}) where {N<:Real}
    return U
end
function intersection(U::Universe{N}, P::AbstractPolyhedron{N}) where {N<:Real}
    return P
end
function intersection(P::AbstractPolyhedron{N}, U::Universe{N}) where {N<:Real}
    return P
end
function intersection(U::Universe{N}, S::AbstractSingleton{N}) where {N<:Real}
    return S
end
function intersection(S::AbstractSingleton{N}, U::Universe{N}) where {N<:Real}
    return S
end
function intersection(X::Interval{N}, U::Universe{N}) where {N<:Real}
    return X
end
function intersection(U::Universe{N}, X::Interval{N}) where {N<:Real}
    return X
end

"""
    intersection(P::AbstractPolyhedron{N}, rm::ResetMap{N}) where {N<:Real}

Return the intersection of a polyhedron and a polyhedral reset map.

### Input

- `P`  -- polyhedron
- `rm` -- polyhedral reset map

### Output

A polyhedron.

### Notes

We assume that `rm` is polyhedral, i.e., has a `constraints_list` method
defined.
"""
function intersection(P::AbstractPolyhedron{N}, rm::ResetMap{N}) where {N<:Real}
    return intersection(P, HPolyhedron(constraints_list(rm)))
end

# symmetric method
function intersection(rm::ResetMap{N}, P::AbstractPolyhedron{N}) where {N<:Real}
    return intersection(P, rm)
end

# more efficient version for polytopic
function intersection(P::AbstractPolyhedron{N},
                      rm::ResetMap{N, <:AbstractPolytope}) where {N<:Real}
    return intersection(P, HPolytope(constraints_list(rm)))
end

# symmetric method
function intersection(rm::ResetMap{N, <:AbstractPolytope},
                      P::AbstractPolyhedron{N}) where {N<:Real}
    return intersection(P, rm)
end

function intersection(U::Universe{N}, X::CartesianProductArray{N}) where {N<:Real}
    return X
end

# symmetric method
function intersection(X::CartesianProductArray{N}, U::Universe{N}) where {N<:Real}
    return intersection(U, X)
end

"""
        intersection(X::CartesianProductArray{N}, Y::CartesianProductArray{N})

Return the intersection between cartesian products of a finite number of convex sets.

### Input

 - `X` -- cartesian product of a finite number of convex sets
 - `Y` -- cartesian product of a finite number of convex sets

### Output

The decomposed set which represents concrete intersection between `X` and `Y`

### Algorithm

This algorithm intersect corresponding blocks between sets.
"""
function intersection(X::CartesianProductArray{N}, Y::CartesianProductArray{N}) where {N<:Real}
    @assert same_block_structure(array(X), array(Y)) "block structure has to be the same"

    return CartesianProductArray([intersection(array(X)[i], array(Y)[i]) for i in eachindex(array(X))])
end

"""
    intersection(cpa::CartesianProductArray{N}, P::AbstractPolyhedron{N}
                ) where {N<:Real}

Compute the intersection of a Cartesian product of a finite number of polyhedral
sets with a polyhedron.

### Input

- `cpa` -- Cartesian product of a finite number of polyhedral sets
- `P`   -- polyhedron

### Output

A Cartesian product of a finite number of polyhedral sets.
See the *Algorithm* section below for details about the structure.

### Notes

The restriction to polyhedral sets in `cpa` only applies to the blocks that are
actually intersected with `P` (see the *Algorithm* section below for details).
All other blocks are not considered by the intersection and remain identical.

### Algorithm

The underlying idea of the algorithm is to exploit the unconstrained dimensions
of `P`.
Without loss of generality, assume that `cpa` has the structure ``X × Y × Z`` such that only the
dimensions of ``Y`` are constrained in ``P``, and denoting a suitable projection
of ``P`` to the dimensions of ``Y`` with ``P|_Y``, we have the following
equivalence:

```math
    (X × Y × Z) ∩ P = X × (Y ∩ P|_Y) × Z
```

Note that ``Y`` may still consist of many blocks.
However, due to the structural restriction of a Cartesian product, we cannot
break down this set further even if ``P|_Y`` is still unconstrained in some
dimensions of blocks in ``Y``.
This would require a restructuring of the dimensions.
Consider this example:

```math
    Y := [0, 1] × [1, 2] × [2, 3]
    P|_Y := x₁ + x₃ ≤ 2
    Y ∩ P|_Y = 0 ≤ x₁ ∧ 1 ≤ x₂ ≤ 2 ∧ 2 ≤ x₃ ∧ x₁ + x₃ ≤ 2
```

Even though the constraints of dimension ``x₂`` are decoupled from the rest,
due to the last constraint the Cartesian product cannot be broken down further.
In particular, the result ``Y ∩ P|_Y`` is a polyhedron in this implementation.

Now we explain the implementation of the above idea.
We first identify the dimensions in which `P` is constrained.
Then we identify the block dimensions of ``X × Y × Z`` such that ``Y`` has
minimal dimension.
Finally, we convert ``Y`` to a polyhedron and intersect it with a suitable
projection of `P`.
"""
function intersection(cpa::CartesianProductArray{N}, P::AbstractPolyhedron{N}
                     ) where {N<:Real}
    # search for the indices of the block trisection into
    # "unconstrained | constrained | unconstrained" (the first and third section
    # may be empty)
    constrained_dims = constrained_dimensions(P)
    if isempty(constrained_dims)
        return cpa
    end
    minimal = minimum(constrained_dims)
    maximal = maximum(constrained_dims)
    lower_bound = 1
    blocks = array(cpa)
    cb_start = 0
    cb_end = length(blocks)
    dim_start = 0
    dim_end = dim(P)
    for (i, block) in enumerate(blocks)
        dim_block = dim(block)
        upper_bound = lower_bound + dim_block - 1
        interval = lower_bound:upper_bound
        if cb_start == 0 && minimal ∈ interval
            dim_start = lower_bound
            cb_start = i
        end
        if maximal ∈ interval
            dim_end = upper_bound
            cb_end = i
            break
        end
        lower_bound = upper_bound + 1
    end
    # compute intersection with constrained blocks
    cap = _intersection_cpa_polyhedron(
        CartesianProductArray(blocks[cb_start:cb_end]), P, dim_start:dim_end)

    # construct result
    result_array = vcat(blocks[1:(cb_start - 1)],  # unconstrained
                        [cap],  # constrained, intersected
                        blocks[(cb_end + 1):length(blocks)])  # unconstrained
    return CartesianProductArray(result_array)
end

function _intersection_cpa_polyhedron(cpa::CartesianProductArray{N},
                                      P::AbstractPolyhedron{N},
                                      constrained_dims) where {N<:Real}
    T = isbounded(cpa) ? HPolytope : HPolyhedron
    hpoly_low_dim = T(constraints_list(cpa))
    cap_low_dim = intersection(hpoly_low_dim, project(P, constrained_dims))
end

# symmetric method
function intersection(P::AbstractPolyhedron{N}, cpa::CartesianProductArray{N}
                     ) where {N<:Real}
    return intersection(cpa, P)
end


text:
"""
    @htl string-expression

Create a `Result` object with string interpolation (`\$`) that uses
context-sensitive hypertext escaping. Before Julia 1.6, interpolated
string literals, e.g. `\$("Strunk & White")`, are treated as errors
since they cannot be reliably detected (see Julia issue #38501).
"""
macro htl(expr)
    if typeof(expr) == String
        return interpolate([expr])
    end
    if !Meta.isexpr(expr, :string)
        throw(DomainError(expr, "a string literal is required"))
    end
    args = expr.args
    for part in expr.args
        if Meta.isexpr(part, :(=))
            throw(DomainError(part,
             "assignments are not permitted in an interpolation"))
        end
    end
    if VERSION < v"1.6.0-DEV"
        # Find cases where we may have an interpolated string literal and
        # raise an exception (till Julia issue #38501 is addressed)
        if length(args) == 1 && args[1] isa String
            throw("interpolated string literals are not supported")
        end
        for idx in 2:length(args)
            if args[idx] isa String && args[idx-1] isa String
                throw("interpolated string literals are not supported")
            end
        end
    end
    return interpolate(expr.args)
end

"""
    @htl_str -> Result

Create a `Result` object with string interpolation (`\$`) that uses
context-sensitive hypertext escaping. Unlike the `@htl` macro, this
string literal does not include escaping feature [1]. To include `\$`
within user content one must write `$`. Observe that `"` and
any other HTML ampersand escape sequence can be used as appropriate.

In this syntax, interpolation is extended beyond regular Julia strings
to handle three additional cases: tuples, named tuples (for attributes),
and generators. See Julia #38734 for the feature request so that this
could also work within the `@htl` macro syntax.

[1] There are also a few edge cases, see `@raw_str` documentation and
Julia #22926 for more detail.
"""
macro htl_str(expr::String)
    # Essentially this is an ad-hoc scanner of the string, splitting
    # it by `$` to find interpolated parts and delegating the hard work
    # to `Meta.parse`, treating everything else as a literal string.
    args = Any[]
    start = idx = 1
    strlen = lastindex(expr)
    while true
        idx = findnext(isequal('$'), expr, start)
        if idx == nothing
           chunk = expr[start:strlen]
           push!(args, expr[start:strlen])
           break
        end
        push!(args, expr[start:prevind(expr, idx)])
        start = nextind(expr, idx)
        (nest, tail) = Meta.parse(expr, start; greedy=false)
        @assert nest != nothing
        if !(expr[start] == '(' || nest isa Symbol)
            throw(DomainError(nest,
             "interpolations must be symbols or parenthesized"))
        end
        if Meta.isexpr(nest, :(=))
            throw(DomainError(nest,
             "assignments are not permitted in an interpolation"))
        end
        if nest isa String
            # this is an interpolated string literal
            nest = Expr(:string, nest)
        end
        push!(args, nest)
        start = tail
    end
    return interpolate(args)
end

"""
    Result(fn)

This object wraps a function produced by the `@htl` macro. This function
prints a the evaluated to the given `io`. This object is also showable
via `"text/html"` so it may be used in an HTML display context.
"""
struct Result
    content::Function

    Result(fn::Function) = new(fn)
end

Result(ob) = Result(io::IO -> print(io, ob))

function Result(xs...)
    Result() do io::IO
        for x in xs
            print(io, x)
        end
    end
end

Base.show(io::IO, h::Result) = h.content(EscapeProxy(io))
Base.show(io::IO, m::MIME"text/html", h::Result) = h.content(EscapeProxy(io))
Base.show(io::EscapeProxy, h::Result) = h.content(io)
Base.show(io::EscapeProxy, m::MIME"text/html", h::Result) = h.content(io)

ext: jl


text:
## Do inference in parallel

using Distributed

## add number of processors desired

addprocs(4)

## Load relevant code into scope for each worker

@everywhere begin
    using Pkg;Pkg.activate(".")
end
@everywhere using JLD2,LinearAlgebra,TransformVariables,Parameters,LogDensityProblems,MCMCChains,Dates,DynamicHMC,BlackBoxOptim,OrdinaryDiffEq,DiffEqCallbacks,Random
@everywhere using Distributions

@everywhere import KenyaCoVSD
@everywhere include("inference_variable_transmissibility.jl");
## Load data

@load("data/linelist_data_with_pos_neg_20feb_to_27apr_c_age.jld2")
@load("data/serological_data_with_20feb_to_31dec_age.jld2")
@load("data/deaths_data20210521.jld2")
@load("data/N_kenya.jld2")

## Parameter transformations and initial KE guess for HMC

trans = as((ct_min1 = as(Real, 0.4, 0.999),
        R₀ = as(Real, 0.0, 7.0),
        ϵ = as(Real, 0.7, 0.999),
        χ₁=as(Real, 0.1, 30.0),χ₂=as(Real, 0.1, 30.0),
        p_test₁=as(Real, 5e-6, 5e-3),p_test₂ = as(Real, 5e-6, 5e-3),
        p_choose1=as(Real, 0.1, 0.999),
        P_eff=as(Real, 0.5, 0.999),
        schooleffect = as(Real, 0.5, 0.999),
        extra_transmissibility = as(Real, 0.5, 2.5),
        influx_exposed_new_variant = as(Real, 0.0, 10e3),
        p_test_boost = as(Real, 1.0, 3.0),
        χ_boost= as(Real, 1.0, 3.0),
        E₀=as(Real, 0.0, 10e3)))

## Transformation of variables
D = Diagonal(0.1*ones(TransformVariables.dimension(trans)))

## Define scenarios
scenarios = [(σ = 0.16,ω = 1/180,ι = 0.5, scenario_name = "reduced_inf"),
              (σ = 0.5,ω = 1/180,ι = 0.5,scenario_name = "_med_susceptibility_reduced_inf"),
              (σ = 1.0,ω = 1/180,ι = 0.5,scenario_name = "_full_susceptibility_reduced_inf"),
              (σ = 1.0,ω = 1/90,ι = 0.5,scenario_name = "_full_susceptibility_fast_wane_reduced_inf")]

##Fit each scenario
@distributed for i in [2]#1:length(scenarios)
    scenario = scenarios[i]
    model = KenyaCoVSD.CoVAreaModel("Nairobi",ll_twogroup_newvariant_infboost_var_trans,KenyaCoVSD.priors_twogroup_newvariant_cities;
                case_data = linelist_data_with_pos_neg,
                sero_data = serological_data,
                death_data = deaths_data,
                pop_data = N_kenya)
    model.prob = KenyaCoVSD.build_prob(Val(:lower_trans))
    inferparameters!(model,2000,trans,0.05,D;serowaningrate = 0.0,σ = scenario.σ,ω = scenario.ω,ι = scenario.ι)

    @save("sensitivity_analysis/Nairobi_model$(scenario.scenario_name).jld2",model)
end


text:
using LinearAlgebra

# k_local(Elastic Modulus,thickness,height,length,number of nodes,row/column one, row/column two)

#creates local stiffness matrix and maps to global
function k_local(k,l,n,rc_1,rc_2)
    k_global = zeros(n,n)
    k_global[rc_1,rc_1] = k
    k_global[rc_1,rc_2] = -k
    k_global[rc_2,rc_1] = -k
    k_global[rc_2,rc_2] = k
    return k_global
end

# constants from problem statement
k1 = 40
k2 = 30
k3 = 60
k4 = 50
# combining all globally mapped matrices
kG = k_local(k1,l,5,1,2)+
     k_local(k2,l,5,2,4)+
     k_local(k3,l,5,3,4)+
     k_local(k4,l,5,4,5)

# Adjusts Global Matrix for fixed-end reactions
function adj_fix(kG,x,n)
    default = zeros(Float64,1,n)
    default[x] = 1
    kG[x,:] = default
    kG[:,x] = transpose(default)
end
adj_fix(kG,1,5)
adj_fix(kG,3,5)

# Force Vector
f = Float64[0 0 0 100 60]
F = transpose(f)

# Displacement vector
u = kG\F
println("Nodal displacement vector (mm):")
display(u)


text:
function watch_files(dir)
    for entry in readdir(dir)
        path = joinpath(dir, entry)
        if isfile(path) && endswith(entry, ".jl")
            println("Watching $path")
            watch_file(callback, path)
        elseif isdir(path)
            watch_files(path)
        else
            #println("Ignoring $path")
        end
    end
end

function callback(fn, args...)
    println("$fn changed, running...\n")
    flush(STDOUT)
    try
        run(`julia test/runtests.jl`)
    catch
    end
end

watch_files(pwd())
wait()


text:
struct FlippedMesh{M}
        base::M
end

Base.:-(m::FlippedMesh) = m.base

struct FlippedMeshCellIterator{CA}
    cellarray::CA
end

cells(m::FlippedMesh) = FlippedMeshCellIterator(m.base.faces)

function Base.iterate(iter::FlippedMeshCellIterator)
    it = iterate(iter.cellarray)
    (it == nothing) ? it : (flip(it[1]), it[2],)
end

function Base.iterate(iter::FlippedMeshCellIterator, state)
    it = iterate(iter.cellarray, state)
    (it == nothing) ? it : (flip(it[1]), it[2],)
end

# Base.start(iter::FlippedMeshCellIterator) = start(iter.cellarray)
# Base.done(iter::FlippedMeshCellIterator, state) = done(iter.cellarray, state)
# function Base.next(iter::FlippedMeshCellIterator, state)
#     item, newstate = next(iter.cellarray, state)
#     flip(item), newstate
# end


text:
# Skipping MacroDefinition: LLVM_FOR_EACH_VALUE_SUBCLASS ( macro ) macro ( Argument ) macro ( BasicBlock ) macro ( InlineAsm ) macro ( MDNode ) macro ( MDString ) macro ( User ) macro ( Constant ) macro ( BlockAddress ) macro ( ConstantAggregateZero ) macro ( ConstantArray ) macro ( ConstantDataSequential ) macro ( ConstantDataArray ) macro ( ConstantDataVector ) macro ( ConstantExpr ) macro ( ConstantFP ) macro ( ConstantInt ) macro ( ConstantPointerNull ) macro ( ConstantStruct ) macro ( ConstantVector ) macro ( GlobalValue ) macro ( GlobalAlias ) macro ( GlobalObject ) macro ( Function ) macro ( GlobalVariable ) macro ( UndefValue ) macro ( Instruction ) macro ( BinaryOperator ) macro ( CallInst ) macro ( IntrinsicInst ) macro ( DbgInfoIntrinsic ) macro ( DbgDeclareInst ) macro ( MemIntrinsic ) macro ( MemCpyInst ) macro ( MemMoveInst ) macro ( MemSetInst ) macro ( CmpInst ) macro ( FCmpInst ) macro ( ICmpInst ) macro ( ExtractElementInst ) macro ( GetElementPtrInst ) macro ( InsertElementInst ) macro ( InsertValueInst ) macro ( LandingPadInst ) macro ( PHINode ) macro ( SelectInst ) macro ( ShuffleVectorInst ) macro ( StoreInst ) macro ( TerminatorInst ) macro ( BranchInst ) macro ( IndirectBrInst ) macro ( InvokeInst ) macro ( ReturnInst ) macro ( SwitchInst ) macro ( UnreachableInst ) macro ( ResumeInst ) macro ( UnaryInstruction ) macro ( AllocaInst ) macro ( CastInst ) macro ( AddrSpaceCastInst ) macro ( BitCastInst ) macro ( FPExtInst ) macro ( FPToSIInst ) macro ( FPToUIInst ) macro ( FPTruncInst ) macro ( IntToPtrInst ) macro ( PtrToIntInst ) macro ( SExtInst ) macro ( SIToFPInst ) macro ( TruncInst ) macro ( UIToFPInst ) macro ( ZExtInst ) macro ( ExtractValueInst ) macro ( LoadInst ) macro ( VAArgInst )
# Skipping MacroDefinition: LLVM_DECLARE_VALUE_CAST ( name ) LLVMValueRef LLVMIsA ## name ( LLVMValueRef Val ) ;
typealias LLVMBool Cint
type LLVMOpaqueMemoryBuffer
end
typealias LLVMMemoryBufferRef Ptr{LLVMOpaqueMemoryBuffer}
type LLVMOpaqueContext
end
typealias LLVMContextRef Ptr{LLVMOpaqueContext}
type LLVMOpaqueModule
end
typealias LLVMModuleRef Ptr{LLVMOpaqueModule}
type LLVMOpaqueType
end
typealias LLVMTypeRef Ptr{LLVMOpaqueType}
type LLVMOpaqueValue
end
typealias LLVMValueRef Ptr{LLVMOpaqueValue}
type LLVMOpaqueBasicBlock
end
typealias LLVMBasicBlockRef Ptr{LLVMOpaqueBasicBlock}
type LLVMOpaqueBuilder
end
typealias LLVMBuilderRef Ptr{LLVMOpaqueBuilder}
type LLVMOpaqueModuleProvider
end
typealias LLVMModuleProviderRef Ptr{LLVMOpaqueModuleProvider}
type LLVMOpaquePassManager
end
typealias LLVMPassManagerRef Ptr{LLVMOpaquePassManager}
type LLVMOpaquePassRegistry
end
typealias LLVMPassRegistryRef Ptr{LLVMOpaquePassRegistry}
type LLVMOpaqueUse
end
typealias LLVMUseRef Ptr{LLVMOpaqueUse}
type LLVMOpaqueDiagnosticInfo
end
typealias LLVMDiagnosticInfoRef Ptr{LLVMOpaqueDiagnosticInfo}
# begin enum ANONYMOUS_1
typealias ANONYMOUS_1 Cint
const LLVMZExtAttribute = int32(1)
const LLVMSExtAttribute = int32(2)
const LLVMNoReturnAttribute = int32(4)
const LLVMInRegAttribute = int32(8)
const LLVMStructRetAttribute = int32(16)
const LLVMNoUnwindAttribute = int32(32)
const LLVMNoAliasAttribute = int32(64)
const LLVMByValAttribute = int32(128)
const LLVMNestAttribute = int32(256)
const LLVMReadNoneAttribute = int32(512)
const LLVMReadOnlyAttribute = int32(1024)
const LLVMNoInlineAttribute = int32(2048)
const LLVMAlwaysInlineAttribute = int32(4096)
const LLVMOptimizeForSizeAttribute = int32(8192)
const LLVMStackProtectAttribute = int32(16384)
const LLVMStackProtectReqAttribute = int32(32768)
const LLVMAlignment = int32(2031616)
const LLVMNoCaptureAttribute = int32(2097152)
const LLVMNoRedZoneAttribute = int32(4194304)
const LLVMNoImplicitFloatAttribute = int32(8388608)
const LLVMNakedAttribute = int32(16777216)
const LLVMInlineHintAttribute = int32(33554432)
const LLVMStackAlignment = int32(469762048)
const LLVMReturnsTwice = int32(536870912)
const LLVMUWTable = int32(1073741824)
const LLVMNonLazyBind = int32(-2147483648)
# end enum ANONYMOUS_1
# begin enum LLVMAttribute
typealias LLVMAttribute Cint
const LLVMZExtAttribute = int32(1)
const LLVMSExtAttribute = int32(2)
const LLVMNoReturnAttribute = int32(4)
const LLVMInRegAttribute = int32(8)
const LLVMStructRetAttribute = int32(16)
const LLVMNoUnwindAttribute = int32(32)
const LLVMNoAliasAttribute = int32(64)
const LLVMByValAttribute = int32(128)
const LLVMNestAttribute = int32(256)
const LLVMReadNoneAttribute = int32(512)
const LLVMReadOnlyAttribute = int32(1024)
const LLVMNoInlineAttribute = int32(2048)
const LLVMAlwaysInlineAttribute = int32(4096)
const LLVMOptimizeForSizeAttribute = int32(8192)
const LLVMStackProtectAttribute = int32(16384)
const LLVMStackProtectReqAttribute = int32(32768)
const LLVMAlignment = int32(2031616)
const LLVMNoCaptureAttribute = int32(2097152)
const LLVMNoRedZoneAttribute = int32(4194304)
const LLVMNoImplicitFloatAttribute = int32(8388608)
const LLVMNakedAttribute = int32(16777216)
const LLVMInlineHintAttribute = int32(33554432)
const LLVMStackAlignment = int32(469762048)
const LLVMReturnsTwice = int32(536870912)
const LLVMUWTable = int32(1073741824)
const LLVMNonLazyBind = int32(-2147483648)
# end enum LLVMAttribute
# begin enum LLVMOpcode
typealias LLVMOpcode Uint32
const LLVMRet = uint32(1)
const LLVMBr = uint32(2)
const LLVMSwitch = uint32(3)
const LLVMIndirectBr = uint32(4)
const LLVMInvoke = uint32(5)
const LLVMUnreachable = uint32(7)
const LLVMAdd = uint32(8)
const LLVMFAdd = uint32(9)
const LLVMSub = uint32(10)
const LLVMFSub = uint32(11)
const LLVMMul = uint32(12)
const LLVMFMul = uint32(13)
const LLVMUDiv = uint32(14)
const LLVMSDiv = uint32(15)
const LLVMFDiv = uint32(16)
const LLVMURem = uint32(17)
const LLVMSRem = uint32(18)
const LLVMFRem = uint32(19)
const LLVMShl = uint32(20)
const LLVMLShr = uint32(21)
const LLVMAShr = uint32(22)
const LLVMAnd = uint32(23)
const LLVMOr = uint32(24)
const LLVMXor = uint32(25)
const LLVMAlloca = uint32(26)
const LLVMLoad = uint32(27)
const LLVMStore = uint32(28)
const LLVMGetElementPtr = uint32(29)
const LLVMTrunc = uint32(30)
const LLVMZExt = uint32(31)
const LLVMSExt = uint32(32)
const LLVMFPToUI = uint32(33)
const LLVMFPToSI = uint32(34)
const LLVMUIToFP = uint32(35)
const LLVMSIToFP = uint32(36)
const LLVMFPTrunc = uint32(37)
const LLVMFPExt = uint32(38)
const LLVMPtrToInt = uint32(39)
const LLVMIntToPtr = uint32(40)
const LLVMBitCast = uint32(41)
const LLVMAddrSpaceCast = uint32(60)
const LLVMICmp = uint32(42)
const LLVMFCmp = uint32(43)
const LLVMPHI = uint32(44)
const LLVMCall = uint32(45)
const LLVMSelect = uint32(46)
const LLVMUserOp1 = uint32(47)
const LLVMUserOp2 = uint32(48)
const LLVMVAArg = uint32(49)
const LLVMExtractElement = uint32(50)
const LLVMInsertElement = uint32(51)
const LLVMShuffleVector = uint32(52)
const LLVMExtractValue = uint32(53)
const LLVMInsertValue = uint32(54)
const LLVMFence = uint32(55)
const LLVMAtomicCmpXchg = uint32(56)
const LLVMAtomicRMW = uint32(57)
const LLVMResume = uint32(58)
const LLVMLandingPad = uint32(59)
# end enum LLVMOpcode
# begin enum LLVMTypeKind
typealias LLVMTypeKind Uint32
const LLVMVoidTypeKind = uint32(0)
const LLVMHalfTypeKind = uint32(1)
const LLVMFloatTypeKind = uint32(2)
const LLVMDoubleTypeKind = uint32(3)
const LLVMX86_FP80TypeKind = uint32(4)
const LLVMFP128TypeKind = uint32(5)
const LLVMPPC_FP128TypeKind = uint32(6)
const LLVMLabelTypeKind = uint32(7)
const LLVMIntegerTypeKind = uint32(8)
const LLVMFunctionTypeKind = uint32(9)
const LLVMStructTypeKind = uint32(10)
const LLVMArrayTypeKind = uint32(11)
const LLVMPointerTypeKind = uint32(12)
const LLVMVectorTypeKind = uint32(13)
const LLVMMetadataTypeKind = uint32(14)
const LLVMX86_MMXTypeKind = uint32(15)
# end enum LLVMTypeKind
# begin enum LLVMLinkage
typealias LLVMLinkage Uint32
const LLVMExternalLinkage = uint32(0)
const LLVMAvailableExternallyLinkage = uint32(1)
const LLVMLinkOnceAnyLinkage = uint32(2)
const LLVMLinkOnceODRLinkage = uint32(3)
const LLVMLinkOnceODRAutoHideLinkage = uint32(4)
const LLVMWeakAnyLinkage = uint32(5)
const LLVMWeakODRLinkage = uint32(6)
const LLVMAppendingLinkage = uint32(7)
const LLVMInternalLinkage = uint32(8)
const LLVMPrivateLinkage = uint32(9)
const LLVMDLLImportLinkage = uint32(10)
const LLVMDLLExportLinkage = uint32(11)
const LLVMExternalWeakLinkage = uint32(12)
const LLVMGhostLinkage = uint32(13)
const LLVMCommonLinkage = uint32(14)
const LLVMLinkerPrivateLinkage = uint32(15)
const LLVMLinkerPrivateWeakLinkage = uint32(16)
# end enum LLVMLinkage
# begin enum LLVMVisibility
typealias LLVMVisibility Uint32
const LLVMDefaultVisibility = uint32(0)
const LLVMHiddenVisibility = uint32(1)
const LLVMProtectedVisibility = uint32(2)
# end enum LLVMVisibility
# begin enum LLVMDLLStorageClass
typealias LLVMDLLStorageClass Uint32
const LLVMDefaultStorageClass = uint32(0)
const LLVMDLLImportStorageClass = uint32(1)
const LLVMDLLExportStorageClass = uint32(2)
# end enum LLVMDLLStorageClass
# begin enum LLVMCallConv
typealias LLVMCallConv Uint32
const LLVMCCallConv = uint32(0)
const LLVMFastCallConv = uint32(8)
const LLVMColdCallConv = uint32(9)
const LLVMWebKitJSCallConv = uint32(12)
const LLVMAnyRegCallConv = uint32(13)
const LLVMX86StdcallCallConv = uint32(64)
const LLVMX86FastcallCallConv = uint32(65)
# end enum LLVMCallConv
# begin enum LLVMIntPredicate
typealias LLVMIntPredicate Uint32
const LLVMIntEQ = uint32(32)
const LLVMIntNE = uint32(33)
const LLVMIntUGT = uint32(34)
const LLVMIntUGE = uint32(35)
const LLVMIntULT = uint32(36)
const LLVMIntULE = uint32(37)
const LLVMIntSGT = uint32(38)
const LLVMIntSGE = uint32(39)
const LLVMIntSLT = uint32(40)
const LLVMIntSLE = uint32(41)
# end enum LLVMIntPredicate
# begin enum LLVMRealPredicate
typealias LLVMRealPredicate Uint32
const LLVMRealPredicateFalse = uint32(0)
const LLVMRealOEQ = uint32(1)
const LLVMRealOGT = uint32(2)
const LLVMRealOGE = uint32(3)
const LLVMRealOLT = uint32(4)
const LLVMRealOLE = uint32(5)
const LLVMRealONE = uint32(6)
const LLVMRealORD = uint32(7)
const LLVMRealUNO = uint32(8)
const LLVMRealUEQ = uint32(9)
const LLVMRealUGT = uint32(10)
const LLVMRealUGE = uint32(11)
const LLVMRealULT = uint32(12)
const LLVMRealULE = uint32(13)
const LLVMRealUNE = uint32(14)
const LLVMRealPredicateTrue = uint32(15)
# end enum LLVMRealPredicate
# begin enum LLVMLandingPadClauseTy
typealias LLVMLandingPadClauseTy Uint32
const LLVMLandingPadCatch = uint32(0)
const LLVMLandingPadFilter = uint32(1)
# end enum LLVMLandingPadClauseTy
# begin enum LLVMThreadLocalMode
typealias LLVMThreadLocalMode Uint32
const LLVMNotThreadLocal = uint32(0)
const LLVMGeneralDynamicTLSModel = uint32(1)
const LLVMLocalDynamicTLSModel = uint32(2)
const LLVMInitialExecTLSModel = uint32(3)
const LLVMLocalExecTLSModel = uint32(4)
# end enum LLVMThreadLocalMode
# begin enum LLVMAtomicOrdering
typealias LLVMAtomicOrdering Uint32
const LLVMAtomicOrderingNotAtomic = uint32(0)
const LLVMAtomicOrderingUnordered = uint32(1)
const LLVMAtomicOrderingMonotonic = uint32(2)
const LLVMAtomicOrderingAcquire = uint32(4)
const LLVMAtomicOrderingRelease = uint32(5)
const LLVMAtomicOrderingAcquireRelease = uint32(6)
const LLVMAtomicOrderingSequentiallyConsistent = uint32(7)
# end enum LLVMAtomicOrdering
# begin enum LLVMAtomicRMWBinOp
typealias LLVMAtomicRMWBinOp Uint32
const LLVMAtomicRMWBinOpXchg = uint32(0)
const LLVMAtomicRMWBinOpAdd = uint32(1)
const LLVMAtomicRMWBinOpSub = uint32(2)
const LLVMAtomicRMWBinOpAnd = uint32(3)
const LLVMAtomicRMWBinOpNand = uint32(4)
const LLVMAtomicRMWBinOpOr = uint32(5)
const LLVMAtomicRMWBinOpXor = uint32(6)
const LLVMAtomicRMWBinOpMax = uint32(7)
const LLVMAtomicRMWBinOpMin = uint32(8)
const LLVMAtomicRMWBinOpUMax = uint32(9)
const LLVMAtomicRMWBinOpUMin = uint32(10)
# end enum LLVMAtomicRMWBinOp
# begin enum LLVMDiagnosticSeverity
typealias LLVMDiagnosticSeverity Uint32
const LLVMDSError = uint32(0)
const LLVMDSWarning = uint32(1)
const LLVMDSRemark = uint32(2)
const LLVMDSNote = uint32(3)
# end enum LLVMDiagnosticSeverity
typealias LLVMFatalErrorHandler Ptr{Void}
typealias LLVMDiagnosticHandler Ptr{Void}
typealias LLVMYieldCallback Ptr{Void}
# Skipping MacroDefinition: LLVM_TARGET ( TargetName ) void LLVMInitialize ## TargetName ## TargetInfo ( void ) ;
# Skipping MacroDefinition: LLVM_ASM_PRINTER ( TargetName ) void LLVMInitialize ## TargetName ## AsmPrinter ( void ) ;
# Skipping MacroDefinition: LLVM_ASM_PARSER ( TargetName ) void LLVMInitialize ## TargetName ## AsmParser ( void ) ;
# Skipping MacroDefinition: LLVM_DISASSEMBLER ( TargetName ) void LLVMInitialize ## TargetName ## Disassembler ( void ) ;
# begin enum LLVMByteOrdering
typealias LLVMByteOrdering Uint32
const LLVMBigEndian = uint32(0)
const LLVMLittleEndian = uint32(1)
# end enum LLVMByteOrdering
type LLVMOpaqueTargetData
end
typealias LLVMTargetDataRef Ptr{LLVMOpaqueTargetData}
type LLVMOpaqueTargetLibraryInfotData
end
typealias LLVMTargetLibraryInfoRef Ptr{LLVMOpaqueTargetLibraryInfotData}
type LLVMOpaqueTargetMachine
end
typealias LLVMTargetMachineRef Ptr{LLVMOpaqueTargetMachine}
type LLVMTarget
end
typealias LLVMTargetRef Ptr{LLVMTarget}
# begin enum LLVMCodeGenOptLevel
typealias LLVMCodeGenOptLevel Uint32
const LLVMCodeGenLevelNone = uint32(0)
const LLVMCodeGenLevelLess = uint32(1)
const LLVMCodeGenLevelDefault = uint32(2)
const LLVMCodeGenLevelAggressive = uint32(3)
# end enum LLVMCodeGenOptLevel
# begin enum LLVMRelocMode
typealias LLVMRelocMode Uint32
const LLVMRelocDefault = uint32(0)
const LLVMRelocStatic = uint32(1)
const LLVMRelocPIC = uint32(2)
const LLVMRelocDynamicNoPic = uint32(3)
# end enum LLVMRelocMode
# begin enum LLVMCodeModel
typealias LLVMCodeModel Uint32
const LLVMCodeModelDefault = uint32(0)
const LLVMCodeModelJITDefault = uint32(1)
const LLVMCodeModelSmall = uint32(2)
const LLVMCodeModelKernel = uint32(3)
const LLVMCodeModelMedium = uint32(4)
const LLVMCodeModelLarge = uint32(5)
# end enum LLVMCodeModel
# begin enum LLVMCodeGenFileType
typealias LLVMCodeGenFileType Uint32
const LLVMAssemblyFile = uint32(0)
const LLVMObjectFile = uint32(1)
# end enum LLVMCodeGenFileType
# begin enum LLVMVerifierFailureAction
typealias LLVMVerifierFailureAction Uint32
const LLVMAbortProcessAction = uint32(0)
const LLVMPrintMessageAction = uint32(1)
const LLVMReturnStatusAction = uint32(2)
# end enum LLVMVerifierFailureAction
# begin enum LLVMLinkerMode
typealias LLVMLinkerMode Uint32
const LLVMLinkerDestroySource = uint32(0)
const LLVMLinkerPreserveSource = uint32(1)
# end enum LLVMLinkerMode
typealias llvm_lto_t Ptr{Void}
# begin enum llvm_lto_status
typealias llvm_lto_status Uint32
const LLVM_LTO_UNKNOWN = uint32(0)
const LLVM_LTO_OPT_SUCCESS = uint32(1)
const LLVM_LTO_READ_SUCCESS = uint32(2)
const LLVM_LTO_READ_FAILURE = uint32(3)
const LLVM_LTO_WRITE_FAILURE = uint32(4)
const LLVM_LTO_NO_TARGET = uint32(5)
const LLVM_LTO_NO_WORK = uint32(6)
const LLVM_LTO_MODULE_MERGE_FAILURE = uint32(7)
const LLVM_LTO_ASM_FAILURE = uint32(8)
const LLVM_LTO_NULL_OBJECT = uint32(9)
# end enum llvm_lto_status
# begin enum llvm_lto_status_t
typealias llvm_lto_status_t Uint32
const LLVM_LTO_UNKNOWN = uint32(0)
const LLVM_LTO_OPT_SUCCESS = uint32(1)
const LLVM_LTO_READ_SUCCESS = uint32(2)
const LLVM_LTO_READ_FAILURE = uint32(3)
const LLVM_LTO_WRITE_FAILURE = uint32(4)
const LLVM_LTO_NO_TARGET = uint32(5)
const LLVM_LTO_NO_WORK = uint32(6)
const LLVM_LTO_MODULE_MERGE_FAILURE = uint32(7)
const LLVM_LTO_ASM_FAILURE = uint32(8)
const LLVM_LTO_NULL_OBJECT = uint32(9)
# end enum llvm_lto_status_t
type LLVMOpaqueObjectFile
end
typealias LLVMObjectFileRef Ptr{LLVMOpaqueObjectFile}
type LLVMOpaqueSectionIterator
end
typealias LLVMSectionIteratorRef Ptr{LLVMOpaqueSectionIterator}
type LLVMOpaqueSymbolIterator
end
typealias LLVMSymbolIteratorRef Ptr{LLVMOpaqueSymbolIterator}
type LLVMOpaqueRelocationIterator
end
typealias LLVMRelocationIteratorRef Ptr{LLVMOpaqueRelocationIterator}
type LLVMOpaqueGenericValue
end
typealias LLVMGenericValueRef Ptr{LLVMOpaqueGenericValue}
type LLVMOpaqueExecutionEngine
end
typealias LLVMExecutionEngineRef Ptr{LLVMOpaqueExecutionEngine}
type LLVMOpaqueMCJITMemoryManager
end
typealias LLVMMCJITMemoryManagerRef Ptr{LLVMOpaqueMCJITMemoryManager}
type LLVMMCJITCompilerOptions
    OptLevel::Uint32
    CodeModel::LLVMCodeModel
    NoFramePointerElim::LLVMBool
    EnableFastISel::LLVMBool
    MCJMM::LLVMMCJITMemoryManagerRef
end
typealias LLVMMemoryManagerAllocateCodeSectionCallback Ptr{Void}
typealias LLVMMemoryManagerAllocateDataSectionCallback Ptr{Void}
typealias LLVMMemoryManagerFinalizeMemoryCallback Ptr{Void}
typealias LLVMMemoryManagerDestroyCallback Ptr{Void}
const LLVMDisassembler_VariantKind_None = 0
const LLVMDisassembler_VariantKind_ARM_HI16 = 1
const LLVMDisassembler_VariantKind_ARM_LO16 = 2
const LLVMDisassembler_VariantKind_ARM64_PAGE = 1
const LLVMDisassembler_VariantKind_ARM64_PAGEOFF = 2
const LLVMDisassembler_VariantKind_ARM64_GOTPAGE = 3
const LLVMDisassembler_VariantKind_ARM64_GOTPAGEOFF = 4
const LLVMDisassembler_VariantKind_ARM64_TLVP = 5
const LLVMDisassembler_VariantKind_ARM64_TLVOFF = 6
const LLVMDisassembler_ReferenceType_InOut_None = 0
const LLVMDisassembler_ReferenceType_In_Branch = 1
const LLVMDisassembler_ReferenceType_In_PCrel_Load = 2
const LLVMDisassembler_ReferenceType_In_ARM64_ADRP = 0x0000000100000001
const LLVMDisassembler_ReferenceType_In_ARM64_ADDXri = 0x0000000100000002
const LLVMDisassembler_ReferenceType_In_ARM64_LDRXui = 0x0000000100000003
const LLVMDisassembler_ReferenceType_In_ARM64_LDRXl = 0x0000000100000004
const LLVMDisassembler_ReferenceType_In_ARM64_ADR = 0x0000000100000005
const LLVMDisassembler_ReferenceType_Out_SymbolStub = 1
const LLVMDisassembler_ReferenceType_Out_LitPool_SymAddr = 2
const LLVMDisassembler_ReferenceType_Out_LitPool_CstrAddr = 3
const LLVMDisassembler_ReferenceType_Out_Objc_CFString_Ref = 4
const LLVMDisassembler_ReferenceType_Out_Objc_Message = 5
const LLVMDisassembler_ReferenceType_Out_Objc_Message_Ref = 6
const LLVMDisassembler_ReferenceType_Out_Objc_Selector_Ref = 7
const LLVMDisassembler_ReferenceType_Out_Objc_Class_Ref = 8
const LLVMDisassembler_ReferenceType_DeMangled_Name = 9
const LLVMDisassembler_Option_UseMarkup = 1
const LLVMDisassembler_Option_PrintImmHex = 2
const LLVMDisassembler_Option_AsmPrinterVariant = 4
const LLVMDisassembler_Option_SetInstrComments = 8
const LLVMDisassembler_Option_PrintLatency = 16
typealias LLVMDisasmContextRef Ptr{Void}
typealias LLVMOpInfoCallback Ptr{Void}
type LLVMOpInfoSymbol1
    Present::Cint
    Name::Ptr{Uint8}
    Value::Cint
end
type LLVMOpInfo1
    AddSymbol::LLVMOpInfoSymbol1
    SubtractSymbol::LLVMOpInfoSymbol1
    Value::Cint
    VariantKind::Cint
end
typealias LLVMSymbolLookupCallback Ptr{Void}
type LLVMOpaquePassManagerBuilder
end
typealias LLVMPassManagerBuilderRef Ptr{LLVMOpaquePassManagerBuilder}


text:
# The SimpleLayer abstract type and methods for it.
# To be `included` in MERAKit.jl.

"""
    SimpleLayer <: Layer

A `SimpleLayer` is a MERA layer that consists of a set of isometries and/or unitaries, and
nothing else. This allows writing convenient generic versions of many methods, reducing code
duplication for the concrete `Layer` types. Every subtype of `SimpleLayer` should implement
the iteration and indexing interfaces to return the various tensors of the layer in the same
order in which the constructor takes them in.
"""
abstract type SimpleLayer <: Layer end

Base.convert(::Type{T}, t::Tuple) where T <: SimpleLayer = T(t...)
Base.copy(layer::SimpleLayer) = typeof(layer)(map(deepcopy, _tuple(layer))...)

@inline function Base.iterate(layer::SimpleLayer, ::Val{i} = Val(1)) where {i}
    t = _tuple(layer)
    if i > length(t)
        return nothing
    else
        return t[i], Val(i+1)
    end
end

Base.length(layer::SimpleLayer) = length(_tuple(layer))

function remove_symmetry(layer::SimpleLayer)
    return baselayertype(layer)(map(remove_symmetry, _tuple(layer))...)
end

function TensorKitManifolds.projectisometric(layer::SimpleLayer)
    return typeof(layer)(map(projectisometric, _tuple(layer))...)
end

function TensorKitManifolds.projectisometric!(layer::SimpleLayer)
    foreach(projectisometric!, _tuple(layer))
    return layer
end

function pseudoserialize(layer::SimpleLayer)
    return repr(typeof(layer)), map(pseudoserialize, _tuple(layer))
end

function depseudoserialize(::Type{T}, data) where T <: SimpleLayer
    return T(map(d->depseudoserialize(d...), data)...)
end

function tensorwise_scale(layer::SimpleLayer, alpha::Number)
    return baselayertype(layer)((_tuple(layer) .* alpha)...)
end

function tensorwise_sum(l1::SimpleLayer, l2::SimpleLayer)
    @assert baselayertype(l1) == baselayertype(l2)
    return baselayertype(l1)((_tuple(l1) .+ _tuple(l2))...)
end

# # # Manifold functions

function TensorKitManifolds.inner(
    l::SimpleLayer, l1::SimpleLayer, l2::SimpleLayer; metric = :euclidean
)
    function custominner(t, t1, t2)
        mtc = isa(t1, Stiefel.StiefelTangent) ? metric : :euclidean
        return inner(t, t1, t2; metric = mtc)
    end
    return sum(custominner.(_tuple(l), _tuple(l1), _tuple(l2)))
end

function TensorKitManifolds.retract(
    l::SimpleLayer, ltan::SimpleLayer, alpha::Real; alg = :exp
)
    ts_and_ttans = retract.(_tuple(l), _tuple(ltan), alpha; alg = alg)
    ts = first.(ts_and_ttans)
    ttans = last.(ts_and_ttans)
    return baselayertype(l)(ts...), baselayertype(l)(ttans...)
end

function TensorKitManifolds.transport!(
    lvec::SimpleLayer, l::SimpleLayer, ltan::SimpleLayer, alpha::Real, lend::SimpleLayer;
    alg = :exp
)
    ttans = transport!.(
        _tuple(lvec), _tuple(l), _tuple(ltan), alpha, _tuple(lend);
        alg = alg
    )
    return baselayertype(l)(ttans...)
end

"""
    gradient_normsq(layer::Layer, env::Layer; metric = :euclidean)

Compute the norm of the gradient, given the enviroment layer `env` and the base point
`layer`.

See also: [`gradient`](@ref)
"""
function gradient_normsq(layer::SimpleLayer, env::SimpleLayer; metric = :euclidean)
    grad = gradient(layer, env; metric = metric)
    return sum(inner(x, z, z; metric = metric) for (x, z) in zip(layer, grad))
end