text:
colsupport(::TriangularLayout{'L'}, A, j) = isempty(j) ? (1:0) : colsupport(triangulardata(A), j) ∩ (minimum(j):size(A,1))
colsupport(::TriangularLayout{'U'}, A, j) = isempty(j) ? (1:0) : colsupport(triangulardata(A), j) ∩ OneTo(maximum(j))
rowsupport(::TriangularLayout{'U'}, A, j) = isempty(j) ? (1:0) : rowsupport(triangulardata(A), j) ∩ (minimum(j):size(A,2))
rowsupport(::TriangularLayout{'L'}, A, j) = isempty(j) ? (1:0) : rowsupport(triangulardata(A), j) ∩ OneTo(maximum(j))

###
# Lmul
###

mulreduce(M::Mul{<:TriangularLayout}) = Lmul(M)
mulreduce(M::Mul{<:TriangularLayout,<:TriangularLayout}) = Lmul(M)
mulreduce(M::Mul{<:Any,<:TriangularLayout}) = Rmul(M)
mulreduce(M::Mul{<:DiagonalLayout,<:TriangularLayout}) = Lmul(M)
mulreduce(M::Mul{<:TriangularLayout,<:DiagonalLayout}) = Rmul(M)

similar(M::Lmul{<:TriangularLayout{'U'},<:TriangularLayout{'U'}}) = UpperTriangular(Matrix{eltype(M)}(undef, size(M)))
similar(M::Lmul{<:TriangularLayout{'L'},<:TriangularLayout{'L'}}) = LowerTriangular(Matrix{eltype(M)}(undef, size(M)))
similar(M::Lmul{<:TriangularLayout{'U','U'},<:TriangularLayout{'U','U'}}) = UnitUpperTriangular(Matrix{eltype(M)}(undef, size(M)))
similar(M::Lmul{<:TriangularLayout{'L','U'},<:TriangularLayout{'L','U'}}) = UnitLowerTriangular(Matrix{eltype(M)}(undef, size(M)))

## Generic triangular multiplication
function materialize!(M::Lmul{<:TriangularLayout{'U','N'}})
    A,B = M.A,M.B
    m, n = size(B, 1), size(B, 2)
    if m != size(A, 1)
        throw(DimensionMismatch("right hand side B needs first dimension of size $(size(A,1)), has size $m"))
    end
    Adata = triangulardata(A)
    for j = rowsupport(B)
        cs = colsupport(B,j)
        for i = cs
            Bij = Adata[i,i]*B[i,j]
            for k = (i + 1:m) ∩ cs ∩ rowsupport(Adata,i)
                Bij += Adata[i,k]*B[k,j]
            end
            B[i,j] = Bij
        end
    end
    B
end

function materialize!(M::Lmul{<:TriangularLayout{'U','U'}})
    A,B = M.A,M.B
    m, n = size(B, 1), size(B, 2)
    if m != size(A, 1)
        throw(DimensionMismatch("right hand side B needs first dimension of size $(size(A,1)), has size $m"))
    end
    Adata = triangulardata(A)
    for j = rowsupport(B)
        cs = colsupport(B,j)
        for i = cs
            Bij = B[i,j]
            for k = (i + 1:m) ∩ cs ∩ rowsupport(Adata,i)
                Bij += Adata[i,k]*B[k,j]
            end
            B[i,j] = Bij
        end
    end
    B
end

function materialize!(M::Lmul{<:TriangularLayout{'L','N'}})
    A,B = M.A,M.B
    m, n = size(B, 1), size(B, 2)
    if m != size(A, 1)
        throw(DimensionMismatch("right hand side B needs first dimension of size $(size(A,1)), has size $m"))
    end
    Adata = triangulardata(A)
    for j = 1:n
        for i = reverse(colsupport(A,colsupport(B,j)))
            Bij = Adata[i,i]*B[i,j]
            for k = (1:i - 1) ∩ rowsupport(Adata,i)
                Bij += Adata[i,k]*B[k,j]
            end
            B[i,j] = Bij
        end
    end
    B
end
function materialize!(M::Lmul{<:TriangularLayout{'L','U'}})
    A,B = M.A,M.B
    m, n = size(B, 1), size(B, 2)
    if m != size(A, 1)
        throw(DimensionMismatch("right hand side B needs first dimension of size $(size(A,1)), has size $m"))
    end
    Adata = triangulardata(A)
    for j = 1:n
        for i = m:-1:1
            Bij = B[i,j]
            for k = (1:i - 1) ∩ rowsupport(Adata,i)
                Bij += Adata[i,k]*B[k,j]
            end
            B[i,j] = Bij
        end
    end
    B
end

@inline function materialize!(M::BlasMatLmulVec{<:TriangularLayout{UPLO,UNIT,<:AbstractColumnMajor},
                                         <:AbstractStridedLayout}) where {UPLO,UNIT}
    A,x = M.A,M.B
    BLAS.trmv!(UPLO, 'N', UNIT, triangulardata(A), x)
end

@inline function materialize!(M::BlasMatLmulVec{<:TriangularLayout{'U',UNIT,<:AbstractRowMajor},
                                   <:AbstractStridedLayout}) where UNIT
    A,x = M.A,M.B
    BLAS.trmv!('L', 'T', UNIT, transpose(triangulardata(A)), x)
end

@inline function materialize!(M::BlasMatLmulVec{<:TriangularLayout{'L',UNIT,<:AbstractRowMajor},
                                   <:AbstractStridedLayout}) where UNIT
    A,x = M.A,M.B
    BLAS.trmv!('U', 'T', UNIT, transpose(triangulardata(A)), x)
end

@inline function materialize!(M::BlasMatLmulVec{<:TriangularLayout{'U',UNIT,<:ConjLayout{<:AbstractRowMajor}},
                                   <:AbstractStridedLayout,<:BlasComplex}) where UNIT
    A,x = M.A,M.B
    BLAS.trmv!('L', 'C', UNIT, triangulardata(A)', x)
end

@inline function materialize!(M::BlasMatLmulVec{<:TriangularLayout{'L',UNIT,<:ConjLayout{<:AbstractRowMajor}},
                                   <:AbstractStridedLayout,<:BlasComplex}) where UNIT
    A,x = M.A,M.B
    BLAS.trmv!('U', 'C', UNIT, triangulardata(A)', x)
end
# Triangular * Matrix

@inline function materialize!(M::BlasMatLmulMat{<:TriangularLayout{UPLO,UNIT,<:AbstractColumnMajor},
                                         <:AbstractStridedLayout, T}) where {UPLO,UNIT,T<:BlasFloat}
    A,x = M.A,M.B
    BLAS.trmm!('L', UPLO, 'N', UNIT, one(T), triangulardata(A), x)
end

@inline function materialize!(M::BlasMatLmulMat{<:TriangularLayout{'L',UNIT,<:AbstractRowMajor},
                                                <:AbstractStridedLayout, T}) where {UNIT,T<:BlasFloat}
    A,x = M.A,M.B
    BLAS.trmm!('L', 'U', 'T', UNIT, one(T), transpose(triangulardata(A)), x)
end

@inline function materialize!(M::BlasMatLmulMat{<:TriangularLayout{'U',UNIT,<:AbstractRowMajor},
                                                <:AbstractStridedLayout, T}) where {UNIT,T<:BlasFloat}
    A,x = M.A,M.B
    BLAS.trmm!('L', 'L', 'T', UNIT, one(T), transpose(triangulardata(A)), x)
end

@inline function materialize!(M::BlasMatLmulMat{<:TriangularLayout{'L',UNIT,<:ConjLayout{<:AbstractRowMajor}},
                                   <:AbstractStridedLayout, T}) where {UNIT,T<:BlasComplex}
    A,x = M.A,M.B
    BLAS.trmm!('L', 'U', 'C', UNIT, one(T), triangulardata(A)', x)
end

@inline function materialize!(M::BlasMatLmulMat{<:TriangularLayout{'U',UNIT,<:ConjLayout{<:AbstractRowMajor}},
                                                <:AbstractStridedLayout, T}) where {UNIT,T<:BlasComplex}
    A,x = M.A,M.B
    BLAS.trmm!('L', 'L', 'C', UNIT, one(T), triangulardata(A)', x)
end

###
# Rmul
###

@inline function materialize!(M::BlasMatRmulMat{<:AbstractStridedLayout,
                                                <:TriangularLayout{UPLO,UNIT,<:AbstractColumnMajor},T}) where {UPLO,UNIT,T<:BlasFloat}
    x,A = M.A,M.B
    BLAS.trmm!('R', UPLO, 'N', UNIT, one(T), triangulardata(A), x)
end

@inline function materialize!(M::BlasMatRmulMat{<:AbstractStridedLayout,
                                                <:TriangularLayout{'L',UNIT,<:AbstractRowMajor},T}) where {UNIT,T<:BlasFloat}
    x,A = M.A,M.B
    BLAS.trmm!('R', 'U', 'T', UNIT, one(T), transpose(triangulardata(A)), x)
end

@inline function materialize!(M::BlasMatRmulMat{<:AbstractStridedLayout,
                                                <:TriangularLayout{'U',UNIT,<:AbstractRowMajor},T}) where {UNIT,T<:BlasFloat}
x,A = M.A,M.B
BLAS.trmm!('R', 'L', 'T', UNIT, one(T), transpose(triangulardata(A)), x)
end

@inline function materialize!(M::BlasMatRmulMat{<:AbstractStridedLayout,
                                                <:TriangularLayout{'L',UNIT,<:ConjLayout{<:AbstractRowMajor}},T}) where {UPLO,UNIT,T<:BlasComplex}
    x,A = M.A,M.B
    BLAS.trmm!('R', 'U', 'C', UNIT, one(T), triangulardata(A)', x)
end

@inline function materialize!(M::BlasMatRmulMat{<:AbstractStridedLayout,
                                                <:TriangularLayout{'U',UNIT,<:ConjLayout{<:AbstractRowMajor}},T}) where {UPLO,UNIT,T<:BlasComplex}
x,A = M.A,M.B
BLAS.trmm!('R', 'L', 'C', UNIT, one(T), triangulardata(A)', x)
end

materialize!(M::MatRmulMat{<:AbstractStridedLayout,<:TriangularLayout}) = LinearAlgebra.rmul!(M.A, M.B)

########
# Ldiv
########

@inline function copyto!(dest::AbstractArray, M::Ldiv{<:Union{TriangularLayout,BidiagonalLayout,DiagonalLayout}})
    A, B = M.A, M.B
    dest ≡ B || copyto!(dest, B)
    ldiv!(A, dest)
end

for UNIT in ('U', 'N')
    for UPLO in ('L', 'U')
        @eval @inline materialize!(M::BlasMatLdivVec{<:TriangularLayout{$UPLO,$UNIT,<:AbstractColumnMajor},
                                            <:AbstractStridedLayout}) =
            BLAS.trsv!($UPLO, 'N', $UNIT, triangulardata(M.A), M.B)
    end

    @eval begin
        @inline materialize!(M::BlasMatLdivVec{<:TriangularLayout{'U',$UNIT,<:AbstractRowMajor},
                                                        <:AbstractStridedLayout}) =
            BLAS.trsv!('L', 'T', $UNIT, transpose(triangulardata(M.A)), M.B)

        @inline materialize!(M::BlasMatLdivVec{<:TriangularLayout{'L',$UNIT,<:AbstractRowMajor},
                                                        <:AbstractStridedLayout}) =
            BLAS.trsv!('U', 'T', $UNIT, transpose(triangulardata(M.A)), M.B)

        @inline materialize!(M::BlasMatLdivVec{<:TriangularLayout{'U',$UNIT,<:ConjLayout{<:AbstractRowMajor}},
                                                        <:AbstractStridedLayout}) =
            BLAS.trsv!('L', 'C', $UNIT, triangulardata(M.A)', M.B)

        @inline materialize!(M::BlasMatLdivVec{<:TriangularLayout{'L',$UNIT,<:ConjLayout{<:AbstractRowMajor}},
                                                        <:AbstractStridedLayout}) =
            BLAS.trsv!('U', 'C', $UNIT, triangulardata(M.A)', M.B)
    end
end

function materialize!(M::MatLdivMat{<:Union{TriangularLayout,BidiagonalLayout}})
    A,X = M.A,M.B
    size(A,2) == size(X,1) || throw(DimensionMismatch("Dimensions must match"))
    @views for j in axes(X,2)
        ldiv!(A, X[:,j])
    end
    X
end

function materialize!(M::MatLdivVec{<:TriangularLayout{'U','N'}})
    A,b = M.A,M.B
    require_one_based_indexing(A, b)
    n = size(A, 2)
    if !(n == length(b))
        throw(DimensionMismatch("second dimension of left hand side A, $n, and length of right hand side b, $(length(b)), must be equal"))
    end
    data = triangulardata(A)
    @inbounds for j in reverse(colsupport(b,1))
        iszero(data[j,j]) && throw(SingularException(j))
        bj = b[j] = data[j,j] \ b[j]
        for i in (1:j-1) ∩ colsupport(data,j)
            b[i] -= data[i,j] * bj
        end
    end
    b
end

function materialize!(M::MatLdivVec{<:TriangularLayout{'U','U'}})
    A,b = M.A,M.B
    require_one_based_indexing(A, b)
    n = size(A, 2)
    if !(n == length(b))
        throw(DimensionMismatch("second dimension of left hand side A, $n, and length of right hand side b, $(length(b)), must be equal"))
    end
    data = triangulardata(A)
    @inbounds for j in reverse(colsupport(b,1))
        iszero(data[j,j]) && throw(SingularException(j))
        bj = b[j]
        for i in (1:j-1) ∩ colsupport(data,j)
            b[i] -= data[i,j] * bj
        end
    end
    b
end

function materialize!(M::MatLdivVec{<:TriangularLayout{'L','N'}})
    A,b = M.A,M.B
    require_one_based_indexing(A, b)
    n = size(A, 2)
    if !(n == length(b))
        throw(DimensionMismatch("second dimension of left hand side A, $n, and length of right hand side b, $(length(b)), must be equal"))
    end
    data = triangulardata(A)
    @inbounds for j in 1:n
        iszero(data[j,j]) && throw(SingularException(j))
        bj = b[j] = data[j,j] \ b[j]
        for i in (j+1:n) ∩ colsupport(data,j)
            b[i] -= data[i,j] * bj
        end
    end
    b
end

function materialize!(M::MatLdivVec{<:TriangularLayout{'L','U'}})
    A,b = M.A,M.B
    require_one_based_indexing(A, b)
    n = size(A, 2)
    if !(n == length(b))
        throw(DimensionMismatch("second dimension of left hand side A, $n, and length of right hand side b, $(length(b)), must be equal"))
    end
    data = triangulardata(A)
    @inbounds for j in 1:n
        iszero(data[j,j]) && throw(SingularException(j))
        bj = b[j]
        for i in (j+1:n) ∩ colsupport(data,j)
            b[i] -= data[i,j] * bj
        end
    end
    b
end

function _bidiag_backsub!(M)
    A,b = M.A, M.B
    N = last(colsupport(b,1))
    dv = diagonaldata(A)
    ev = supdiagonaldata(A)
    b[N] = bj1 = dv[N]\b[N]

    @inbounds for j = (N - 1):-1:1
        bj  = b[j]
        bj -= ev[j] * bj1
        dvj = dv[j]
        if iszero(dvj)
            throw(SingularException(j))
        end
        bj   = dvj\bj
        b[j] = bj1 = bj
    end

    b
end

function _bidiag_forwardsub!(M)
    A, b = M.A, M.B
    dv = diagonaldata(A)
    ev = subdiagonaldata(A)
    N = length(b)
    b[1] = bj1 = dv[1]\b[1]
    @inbounds for j = 2:N
        bj  = b[j]
        bj -= ev[j - 1] * bj1
        dvj = dv[j]
        if iszero(dvj)
            throw(SingularException(j))
        end
        bj   = dvj\bj
        b[j] = bj1 = bj
    end
    b
end

#Generic solver using naive substitution, based on LinearAlgebra/src/bidiag.jl
function materialize!(M::MatLdivVec{<:BidiagonalLayout})
    A,b = M.A,M.B
    require_one_based_indexing(A, b)
    N = size(A, 2)
    if N != length(b)
        throw(DimensionMismatch("second dimension of A, $N, does not match one of the length of b, $(length(b))"))
    end

    if N == 0
        return b
    end

    if bidiagonaluplo(A) == 'L' #do forward substitution
        _bidiag_forwardsub!(M)
    else #do backward substitution
        _bidiag_backsub!(M)
    end
end


text:
"""
*Compute the bounds of classical Bell scenarios.*

The `LocalPolytope.jl` module provides tools for characterizing the convex polyhedral
structure of the correlations attainable using classical resources in a Bell scenario.

Each distinct strategy for a given Bell scenario corresponds to a point in vector
space.
The complete set of attainable strategies for a particular Bell scenario is denoted
``\\mathbf{P}``.
Given the constraints of normalization and non-negativity, the extreme points of
``\\mathbf{P}`` correspond to deterministic strategies (matrices with 0,1 elements).
Furthermore, when shared randomness is used in a Bell scenario, any two strategies
can be mixed together in a convex combination.
Hence, ``\\mathbf{P}`` is a convex polyhedron referred to as the *local polytope*.

A convex polytope has two equivalent descriptions:

1. **V-Description**: The polytope is the convex hull of a set of extreme points
    known as vertices ``\\mathbf{V}``,
    ``\\quad\\mathbf{P} = \\text{conv}(\\mathbf{V})``.
2. **H-Description**: The polytope is the intersection of a set of linear half-spaces
    known as facets ``\\mathbf{F}``,
    ``\\quad\\mathbf{P} = \\cap\\mathbf{F}``.

These two descriptions are equivalent,

```math
\\mathbf{P} = \\text{conv}(\\mathbf{V}) = \\cap\\mathbf{F},
```

and there exists a transformation between the set of vertices and facets
``\\mathbf{V} \\leftrightarrow \\mathbf{F}``.

For a given Bell scenario, the V-Description is typically easier to compute,
however, the H-Description describes the bounds of the local polytope as
as linear inequalities known as Bell inequalities.
Bell inequalities are important because they provide a simple test to verify that
a strategy is not contained by the local polyotpe.
That is, if a Bell inequality is violated, the classical resources considered for
the local polytope are not sufficient to reproduce the strategy.
Hence, a Bell violation witnesses the use of resources with greater operational value.
Bell violations are often used to characterize the advantages of quantum resources,
however, it is important to note that a Bell violation can simply mean that more classical
resources were used than anticipated.

### Module Exports:
* [`vertices`](@ref) - Compute the set of extreme-points for the local polytope.
* [`num_vertices`](@ref) - The number of vertices for the local polytope.
* [`vrep`](@ref) - Construct a [`Polyhedron`](https://juliapolyhedra.github.io/Polyhedra.jl/stable/polyhedron/)
    in the vertex representation.
* [`facets`](@ref) - Compute the linear inequalities which bound the local polytope.
* [`generator_vertex`](@ref) - Provide a canonical form for a vertex.
* [`generator_facet`](@ref) - Provide a canonical form for a facet.
* [`adjacency_decomposition`](@ref) - Efficiently compute the generator facets for the local
    polytope using the adjacency decomposition technique.
"""
module LocalPolytope

using LinearAlgebra, Combinatorics
using XPORTA, Polyhedra

import Polyhedra: vrep

using ..BellScenario

include("./vertices.jl")

"""
    vrep(scenario::Scenario; vertices_kwargs...) :: XPORTA.Polyhedron

Constructs a `Polyhedron` using the vertex representation.
See [Polyhedra.jl](https://github.com/JuliaPolyhedra/Polyhedra.jl)
for more details.
The `vertices_kwargs` keyword arguments are passed through to the `vertices`
function for each `Scenario`.

!!! note "Return Type"
    This function differs from the Polyhedra.jl implementation in that it returns
    a `Polyhedron` type rather than a `V-Representation`.
    This is done to reduce the number of steps required to construct a polyhedron.
"""
function vrep(scenario::Scenario; vertices_kwargs...) :: XPORTA.Polyhedron
    polyhedron(vrep(vertices(scenario; vertices_kwargs...)), XPORTA.Library())
end

include("./facets.jl")
include("./generators.jl")
include("./adjacency_decomposition.jl")
include("./utils.jl")

# legacy code
include("../Legacy/LocalPolytope.jl")

end

ext: jl


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/AsterReader.jl/blob/master/LICENSE

using AsterReader, Test
using AsterReader: RMEDFile, aster_read_data

datadir = first(splitext(basename(@__FILE__)))

@testset "read nodal field from result rmed file" begin
    rmedfile = joinpath(datadir, "rings.rmed")
    rmed = RMEDFile(rmedfile)
    temp = aster_read_data(rmed, "TEMP"; info_fields=true)
    @test isapprox(temp[15], 1.0)
    @test isapprox(temp[95], 2.0)
end

@testset "try to read gauss point field from result rmed file" begin
    rmedfile = joinpath(datadir, "rings.rmed")
    rmed = RMEDFile(rmedfile)
    @test_throws Exception aster_read_data(rmed, "FLUX_ELGA")
end


text:
"""
    Bisection <: UnivariateAlgorithm

Bisection search is a root-finding algorithm by maintainining a bracket of opposing
derivative signs. At each iteration the algorithm shrinks the bracket by half until a
critical point is located or the bracket is sufficiently small.
"""
struct Bisection <: UnivariateAlgorithm end

function _optimize(f, g, lower::T, upper::T, alg::Bisection;
        reltol, abstol, maxiter) where {T}
    y_lower, y_upper = g(lower), g(upper)
    sign(y_lower) == sign(y_upper) &&
        error("Bisection requires initial bracket with different derivative signs.")
    y_lower == 0 && (upper = lower)
    y_upper == 0 && (lower = upper)
    x = T(NaN)
    converged = false
    iter = 0
    while iter < maxiter
        iter += 1
        x = (lower + upper) / 2
        x_tol = reltol * abs(x) + abstol
        y = g(x)
        if upper - lower < 2x_tol || y == 0
            converged = true
            break
        end
        if sign(y) == sign(y_lower)
            lower = x
        else
            upper = x
        end
    end
    return Solution(converged, iter, x, f(x))
end

function _optimize(f, lower, upper, alg::Bisection; kwargs...)
    g(x) = Zygote.gradient(f, x)[1]
    _optimize(f, g, lower, upper, alg; kwargs...)
end