text:
# V = G Q
struct BandedEigenvectors{T} <: AbstractMatrix{T}
    G::Vector{Givens{T}}
    Q::Matrix{T}
end

size(B::BandedEigenvectors) = size(B.Q)
getindex(B::BandedEigenvectors, i, j) = Matrix(B)[i, j]

# V = S⁻¹ Q W
struct BandedGeneralizedEigenvectors{T,M<:AbstractMatrix{T}} <: AbstractMatrix{T}
    S::SplitCholesky{T,M}
    Q::Vector{Givens{T}}
    W::BandedEigenvectors{T}
end

size(B::BandedGeneralizedEigenvectors) = size(B.W)
getindex(B::BandedGeneralizedEigenvectors, i, j) = Matrix(B)[i, j]

convert(::Type{Eigen{T, T, Matrix{T}, Vector{T}}}, F::Eigen{T, T, BandedEigenvectors{T}, Vector{T}}) where T = Eigen(F.values, Matrix(F.vectors))
convert(::Type{GeneralizedEigen{T, T, Matrix{T}, Vector{T}}}, F::GeneralizedEigen{T, T, BandedGeneralizedEigenvectors{T}, Vector{T}}) where T = GeneralizedEigen(F.values, Matrix(F.vectors))

compress(F::Eigen{T, T, BandedEigenvectors{T}, Vector{T}}) where T = convert(Eigen{T, T, Matrix{T}, Vector{T}}, F)
compress(F::GeneralizedEigen{T, T, BandedGeneralizedEigenvectors{T}, Vector{T}}) where T = convert(GeneralizedEigen{T, T, Matrix{T}, Vector{T}}, F)

eigen(A::Symmetric{T,<:BandedMatrix{T}}) where T <: Real = eigen!(copy(A))
eigen(A::Symmetric{T,<:BandedMatrix{T}}, B::Symmetric{T,<:BandedMatrix{T}}) where T <: Real = eigen!(copy(A), copy(B))

function eigen!(A::Symmetric{T,<:BandedMatrix{T}}) where T <: Real
    N = size(A, 1)
    KD = bandwidth(A)
    D = Vector{T}(undef, N)
    E = Vector{T}(undef, N-1)
    G = Vector{Givens{T}}(undef, 0)
    WORK = Vector{T}(undef, N)
    AB = symbandeddata(A)
    sbtrd!('V', A.uplo, N, KD, AB, D, E, G, WORK)
    Λ, Q = eigen(SymTridiagonal(D, E))
    Eigen(Λ, BandedEigenvectors(G, Q))
end

function eigen!(A::Symmetric{T,<:BandedMatrix{T}}, B::Symmetric{T,<:BandedMatrix{T}}) where T <: Real
    S = splitcholesky!(B)
    N = size(A, 1)
    KA = bandwidth(A)
    KB = bandwidth(B)
    Q = Vector{Givens{T}}(undef, 0)
    WORK = Vector{T}(undef, 2*N)
    AB = symbandeddata(A)
    BB = symbandeddata(B)
    sbgst!('V', A.uplo, N, KA, KB, AB, BB, Q, WORK)
    Λ, W = eigen!(A)
    GeneralizedEigen(Λ, BandedGeneralizedEigenvectors(S, Q, W))
end

function Matrix(B::BandedEigenvectors)
    Q = copy(B.Q)
    G = B.G
    for k in length(G):-1:1
        lmul!(G[k], Q)
    end
    return Q
end

function Matrix(B::BandedGeneralizedEigenvectors)
    V = Matrix(B.W)
    Q = B.Q
    for k in length(Q):-1:1
        lmul!(Q[k], V)
    end
    ldiv!(B.S, V)
    return V
end

function compress!(F::Eigen{T, T, BandedEigenvectors{T}, Vector{T}}) where T
    Q = F.vectors.Q
    G = F.vectors.G
    for k in length(G):-1:1
        lmul!(G[k], Q)
        pop!(G)
    end
    F
end

function mul!(y::Array{T,N}, B::BandedEigenvectors{T}, x::Array{T,N}) where {T,N}
    mul!(y, B.Q, x)
    G = B.G
    for k in length(G):-1:1
        lmul!(G[k], y)
    end
    y
end

function mul!(y::Array{T,N}, B::Adjoint{T,BandedEigenvectors{T}}, x::Array{T,N}) where {T,N}
    Q = B.parent.Q
    G = B.parent.G
    if length(G) > 0
        z = copy(x)
        for k in 1:length(G)
            lmul!(G[k]', z)
        end
        mul!(y, Q', z)
    else
        mul!(y, Q', x)
    end
    y
end

function mul!(y::Array{T,N}, B::Transpose{T,BandedEigenvectors{T}}, x::Array{T,N}) where {T,N}
    Q = B.parent.Q
    G = B.parent.G
    if length(G) > 0
        z = copy(x)
        for k in 1:length(G)
            lmul!(G[k]', z)
        end
        mul!(y, transpose(Q), z)
    else
        mul!(y, transpose(Q), x)
    end
    y
end

function mul!(y::Array{T,N}, B::BandedGeneralizedEigenvectors{T}, x::Array{T,N}) where {T,N}
    mul!(y, B.W, x)
    Q = B.Q
    for k in length(Q):-1:1
        lmul!(Q[k], y)
    end
    ldiv!(B.S, y)
end

function mul!(y::Array{T,N}, B::Adjoint{T,BandedGeneralizedEigenvectors{T,M}}, x::Array{T,N}) where {T,M,N}
    z = copy(x)
    ldiv!(B.parent.S', z)
    Q = B.parent.Q
    for k in 1:length(Q)
        lmul!(Q[k]', z)
    end
    mul!(y, B.parent.W', z)
end

function ldiv!(y::Array{T,N}, B::BandedGeneralizedEigenvectors{T}, x::Array{T,N}) where {T,N}
    z = copy(x)
    lmul!(B.S, z)
    Q = B.Q
    for k in 1:length(Q)
        lmul!(Q[k]', z)
    end
    mul!(y, B.W', z)
end

function mul!(y::Array{T,N}, x::Array{T,N}, B::BandedEigenvectors{T}) where {T,N}
    x .= x'
    mul!(y, B', x)
    x .= x'
    y .= y'
end

function mul!(y::Array{T,N}, x::Array{T,N}, B::BandedGeneralizedEigenvectors{T}) where {T,N}
    x .= x'
    mul!(y, B', x)
    x .= x'
    y .= y'
end

#
# The following is a Julia translation of *SBTRD.f in LAPACK that allows
# extraction of the Givens rotations.
#
#  Definition:
#  ===========
#
#       SUBROUTINE DSBTRD( VECT, UPLO, N, KD, AB, LDAB, D, E, Q, LDQ,
#                          WORK, INFO )
#
#       .. Scalar Arguments ..
#       CHARACTER          UPLO, VECT
#       INTEGER            INFO, KD, LDAB, LDQ, N
#       ..
#       .. Array Arguments ..
#       DOUBLE PRECISION   AB( LDAB, * ), D( * ), E( * ), Q( LDQ, * ),
#      $                   WORK( * )
#       ..
function sbtrd!(VECT::Char, UPLO::Char,
                N::Int, KD::Int, AB::AbstractMatrix{T},
                D::AbstractVector{T}, E::AbstractVector{T}, Q::Vector{Givens{T}},
                WORK::AbstractVector{T}) where T
    require_one_based_indexing(AB)
    chkstride1(AB)
    chkuplo(UPLO)
    chkvect(VECT)
    M    = size(AB,2)
    M ≠ N && throw(ArgumentError("Matrix must be square"))
    size(AB,1) < KD+1 && throw(ArgumentError("Not enough bands"))

    # In Julia, pointers to entries of a SubArray ignore stride.
    # For example, if `A = view(Matrix{Float64}(I, 10, 10), 1:2:9, 1:2:9)`,
    # then `unsafe_load(pointer(A, 3+5*2)) == 1`, even if `stride(A, 2) == 20.`
    # Beyond the use of `LDAB` in Julia pointers, there are two uses of the
    # actual (Fortran) `LDABF` for increment-setting in the LAPACK calls.
    LDAB = size(AB, 1)
    LDABF = max(1, stride(AB, 2))

    ZERO = zero(T)
    TEMP1 = Ref{T}()
    TEMP2 = Ref{T}()
    TEMP3 = Ref{T}()
    KD1 = KD + 1
    KDM1 = KD - 1
    INCX = LDABF - 1

    INCA = KD1*LDABF
    KDN = min(N-1, KD)

    UPPER = UPLO == 'U'
    WANTQ = true

    if UPPER
        if KD > 1
            # Reduce to tridiagonal form, working with upper triangle
            NR = 0
            J1 = KDN + 2
            J2 = 1
            for I = 1:N-2
                # Reduce i-th row of matrix to tridiagonal form
                for K = KDN+1:-1:2
                    J1 = J1 + KDN
                    J2 = J2 + KDN
                    if NR > 0
                        # generate plane rotations to annihilate nonzero
                        # elements which have been created outside the band
                        largv!(NR, pointer(AB, 1+LDAB*(J1-2)), INCA, pointer(WORK, J1), KD1, pointer(D, J1), KD1)
                        # apply rotations from the right
                        # Dependent on the the number of diagonals either
                        # DLARTV or DROT is used
                        if NR ≥ 2*KD-1
                            for L = 1:KD-1
                                lartv!(NR, pointer(AB, L+1+LDAB*(J1-2)), INCA, pointer(AB, L+LDAB*(J1-1)), INCA, pointer(D, J1), pointer(WORK, J1), KD1)
                            end
                        else
                            JEND = J1 + (NR-1)*KD1
                            for JINC = J1:KD1:JEND
                                rot!(KDM1, pointer(AB, 2+LDAB*(JINC-2)), 1, pointer(AB, 1+LDAB*(JINC-1)), 1, D[JINC], WORK[JINC])
                            end
                        end
                    end

                    if K > 2
                        if K ≤ N-I+1
                            # generate plane rotation to annihilate a(i,i+k-1)
                            # within the band
                            lartg!(AB[KD-K+3, I+K-2], AB[KD-K+2, I+K-1], TEMP1, TEMP2, TEMP3)
                            D[I+K-1] = TEMP1[]
                            WORK[I+K-1] = TEMP2[]
                            AB[KD-K+3, I+K-2] = TEMP3[]
                            # apply rotation from the right
                            rot!(K-3, pointer(AB, KD-K+4+LDAB*(I+K-3)), 1, pointer(AB, KD-K+3+LDAB*(I+K-2)), 1, D[I+K-1], WORK[I+K-1])
                        end
                        NR = NR + 1
                        J1 = J1 - KDN - 1
                    end

                    # apply plane rotations from both sides to diagonal
                    # blocks
                    if NR > 0
                        lar2v!(NR, pointer(AB, KD1+LDAB*(J1-2)), pointer(AB, KD1+LDAB*(J1-1)), pointer(AB, KD+LDAB*(J1-1)), INCA, pointer(D, J1), pointer(WORK, J1), KD1)
                    end
                    # apply plane rotations from the left
                    if NR > 0
                        if 2*KD-1 < NR
                            # Dependent on the the number of diagonals either
                            # DLARTV or DROT is used
                            for L = 1:KD-1
                                if J2+L > N
                                    NRT = NR - 1
                                else
                                    NRT = NR
                                end
                                if NRT > 0
                                    lartv!(NRT, pointer(AB, KD-L+LDAB*(J1+L-1)), INCA, pointer(AB, KD-L+1+LDAB*(J1+L-1)), INCA, pointer(D, J1), pointer(WORK, J1), KD1)
                                end
                            end
                        else
                            J1END = J1 + KD1*( NR-2 )
                            if J1END ≥ J1
                                for JIN = J1:KD1:J1END
                                    rot!(KDM1, pointer(AB, KDM1+LDAB*JIN), INCX, pointer(AB, KD+LDAB*JIN), INCX, D[JIN], WORK[JIN])
                                end
                            end
                            LEND = min( KDM1, N-J2 )
                            LAST = J1END + KD1
                            if LEND > 0
                                rot!(LEND, pointer(AB, KDM1+LDAB*LAST), INCX, pointer(AB, KD+LDAB*LAST), INCX, D[LAST], WORK[LAST])
                            end
                        end
                    end

                    if WANTQ
                        for J = J1:KD1:J2
                            push!(Q, Givens(J-1, J, D[J], -WORK[J]))
                        end
                    end
                    if J2+KDN > N
                        # adjust J2 to keep within the bounds of the matrix
                        NR = NR - 1
                        J2 = J2 - KDN - 1
                    end

                    for J = J1:KD1:J2
                        # create nonzero element a(j-1,j+kd) outside the band
                        # and store it in WORK
                        WORK[J+KD] = WORK[J]*AB[1, J+KD]
                        AB[1, J+KD] = D[J]*AB[1, J+KD]
                    end
                end
            end
        end # if

        if KD > 0
            # copy off-diagonal elements to E
            for I = 1:N-1
                E[I] = AB[KD, I+1]
            end
        else
            # set E to zero if original matrix was diagonal
            for I = 1:N-1
                E[I] = ZERO
            end
        end
        # copy diagonal elements to D
        for I = 1:N
            D[I] = AB[KD1, I]
        end
    else # if UPPER
        if KD > 1
            # Reduce to tridiagonal form, working with lower triangle
            NR = 0
            J1 = KDN + 2
            J2 = 1
            for I = 1:N-2
                # Reduce i-th column of matrix to tridiagonal form
                for K = KDN+1:-1:2
                    J1 = J1 + KDN
                    J2 = J2 + KDN
                    if NR > 0
                        # generate plane rotations to annihilate nonzero
                        # elements which have been created outside the band
                        largv!(NR, pointer(AB, KD1+LDAB*(J1-KD1-1)), INCA, pointer(WORK, J1), KD1, pointer(D, J1), KD1)
                        # apply plane rotations from one side
                        # Dependent on the the number of diagonals either
                        # DLARTV or DROT is used
                        if NR > 2*KD-1
                            for L = 1:KD-1
                                lartv!(NR, pointer(AB, KD1-L+LDAB*(J1-KD1+L-1)), INCA, pointer(AB, KD1-L+1+LDAB*(J1-KD1+L-1)), INCA, pointer(D, J1), pointer(WORK, J1), KD1)
                            end
                        else
                            JEND = J1 + KD1*(NR-1)
                            for JINC = J1:KD1:JEND
                                rot!(KDM1, pointer(AB, KD+LDAB*(JINC-KD1)), INCX, pointer(AB, KD1+LDAB*(JINC-KD1)), INCX, D[JINC], WORK[JINC])
                            end
                        end
                    end

                    if K > 2
                        if K ≤ N-I+1
                            # generate plane rotation to annihilate a(i+k-1,i)
                            # within the band
                            lartg!(AB[K-1, I], AB[K, I], TEMP1, TEMP2, TEMP3)
                            D[I+K-1] = TEMP1[]
                            WORK[I+K-1] = TEMP2[]
                            AB[K-1, I] = TEMP3[]
                            # apply rotation from the right
                            rot!(K-3, pointer(AB, K-2+LDAB*I), INCX, pointer(AB, K-1+LDAB*I), INCX, D[I+K-1], WORK[I+K-1])
                        end
                        NR = NR + 1
                        J1 = J1 - KDN - 1
                    end

                    # apply plane rotations from both sides to diagonal
                    # blocks
                    if NR > 0
                        lar2v!(NR, pointer(AB, 1+LDAB*(J1-2)), pointer(AB, 1+LDAB*(J1-1)), pointer(AB, 2+LDAB*(J1-2)), INCA, pointer(D, J1), pointer(WORK, J1), KD1)
                    end
                    # apply plane rotations from the right
                    # Dependent on the the number of diagonals either
                    # DLARTV or DROT is used

                    if NR > 0
                        if NR > 2*KD-1
                            for L = 1:KD-1
                                if J2+L > N
                                    NRT = NR - 1
                                else
                                    NRT = NR
                                end
                                if NRT > 0
                                    lartv!(NRT, pointer(AB, L+2+LDAB*(J1-2)), INCA, pointer(AB, L+1+LDAB*(J1-1)), INCA, pointer(D, J1), pointer(WORK, J1), KD1)
                                end
                            end
                        else
                            J1END = J1 + KD1*(NR-2)
                            if J1END ≥ J1
                                for J1INC = J1:KD1:J1END
                                    rot!(KDM1, pointer(AB, 3+LDAB*(J1INC-2)), 1, pointer(AB, 2+LDAB*(J1INC-1)), 1, D[J1INC], WORK[J1INC])
                                end
                            end
                            LEND = min(KDM1, N-J2)
                            LAST = J1END + KD1
                            if LEND > 0
                                rot!(LEND, pointer(AB, 3+LDAB*(LAST-2)), 1, pointer(AB, 2+LDAB*(LAST-1)), 1, D[LAST], WORK[LAST])
                            end
                        end
                    end

                    if WANTQ
                        for J = J1:KD1:J2
                            push!(Q, Givens(J-1, J, D[J], -WORK[J]))
                        end
                    end

                    if J2+KDN > N
                        # adjust J2 to keep within the bounds of the matrix
                        NR = NR - 1
                        J2 = J2 - KDN - 1
                    end

                    for J = J1:KD1:J2
                        # create nonzero element a(j+kd,j-1) outside the band
                        # and store it in WORK
                        WORK[J+KD] = WORK[J]*AB[KD1, J]
                        AB[KD1, J] = D[J]*AB[KD1, J]
                    end
                end
            end
        end # if
        if KD > 0
            # copy off-diagonal elements to E
            for I = 1:N-1
                E[I] = AB[2, I]
            end
        else
            # set E to zero if original matrix was diagonal
            for I = 1:N-1
                E[I] = ZERO
            end
        end
        # copy diagonal elements to D
        for I = 1:N
            D[I] = AB[1, I]
        end
    end # if UPPER
end

#
# The following is a Julia translation of *SBGST.f in LAPACK that allows
# extraction of the Givens rotations.
#
#  Definition:
#  ===========
#
#       SUBROUTINE DSBGST( VECT, UPLO, N, KA, KB, AB, LDAB, BB, LDBB, X,
#                          LDX, WORK, INFO )
#
#       .. Scalar Arguments ..
#       CHARACTER          UPLO, VECT
#       INTEGER            INFO, KA, KB, LDAB, LDBB, LDX, N
#       ..
#       .. Array Arguments ..
#       DOUBLE PRECISION   AB( LDAB, * ), BB( LDBB, * ), WORK( * ),
#      $                    X( LDX, * )
#       ..
function sbgst!(VECT::Char, UPLO::Char,
                N::Int, KA::Int, KB::Int, AB::AbstractMatrix{T},
                BB::AbstractMatrix{T}, X::Vector{Givens{T}},
                WORK::AbstractVector{T}) where T
    require_one_based_indexing(AB)
    require_one_based_indexing(BB)
    chkstride1(AB)
    chkstride1(BB)
    chkuplo(UPLO)
    chkvect(VECT)
    size(AB,2) == size(BB,2) == N || throw(ArgumentError("Matrices must be square"))
    size(AB,1) < KA+1 && throw(ArgumentError("Not enough bands in AB"))
    size(BB,1) < KB+1 && throw(ArgumentError("Not enough bands in BB"))

    # In Julia, pointers to entries of a SubArray ignore stride.
    # For example, if `A = view(Matrix{Float64}(I, 10, 10), 1:2:9, 1:2:9)`,
    # then `unsafe_load(pointer(A, 3+5*2)) == 1`, even if `stride(A, 2) == 20.`
    # Beyond the use of `LDAB` in Julia pointers, there are two uses of the
    # actual (Fortran) `LDABF` for increment-setting in the LAPACK calls.
    LDAB = size(AB, 1)
    LDABF = max(1, stride(AB, 2))

    ZERO = zero(T)
    TEMP1 = Ref{T}()
    TEMP2 = Ref{T}()
    TEMP3 = Ref{T}()
    KA1 = KA + 1
    KB1 = KB + 1
    INCA = LDABF*KA1
    M = ( N+KB ) ÷ 2

    UPPER = UPLO == 'U'
    WANTX = true

    UPDATE = true
    I = N + 1
    @label L10
    if UPDATE
        I = I - 1
        KBT = min( KB, I-1 )
        I0 = I - 1
        I1 = min( N, I+KA )
        I2 = I - KBT + KA1
        if I < M+1
            UPDATE = false
            I = I + 1
            I0 = M
            if KA == 0
                @goto L480
            end
            @goto L10
        end
    else
        I = I + KA
        if I > N-1
            @goto L480
        end
    end

    if UPPER
        # Transform A, working with the upper triangle
        if UPDATE
            # Form inv(S(i))**T * A * inv(S(i))
            BII = BB[KB1, I]
            for J = I:I1
                AB[I-J+KA1, J] = AB[I-J+KA1, J] / BII
            end
            for J = max(1, I-KA):I
                AB[J-I+KA1, I] = AB[J-I+KA1, I] / BII
            end
            for K = I-KBT:I-1
                for J = I-KBT:K
                    AB[J-K+KA1, K] = AB[J-K+KA1, K] - BB[J-I+KB1, I]*AB[K-I+KA1, I] - BB[K-I+KB1, I]*AB[J-I+KA1, I] + AB[KA1, I]*BB[J-I+KB1, I]*BB[K-I+KB1, I]
                end
                for J = max(1, I-KA):I-KBT-1
                    AB[J-K+KA1, K] = AB[J-K+KA1, K] - BB[K-I+KB1, I]*AB[J-I+KA1, I]
                end
            end
            for J = I:I1
                for K = max(J-KA, I-KBT):I-1
                    AB[K-J+KA1, J] = AB[K-J+KA1, J] - BB[K-I+KB1, I]*AB[I-J+KA1, J]
                end
            end
            # store a(i,i1) in RA1 for use in next loop over K
            RA1 = AB[I-I1+KA1, I1]
        end
        # Generate and apply vectors of rotations to chase all the
        # existing bulges KA positions down toward the bottom of the
        # band
        for K = 1:KB-1
            if UPDATE
                # Determine the rotations which would annihilate the bulge
                # which has in theory just been created
                if I-K+KA < N && I-K > 1
                    # generate rotation to annihilate a(i,i-k+ka+1)
                    lartg!(AB[K+1, I-K+KA], RA1, TEMP1, TEMP2, TEMP3)
                    WORK[N+I-K+KA-M] = TEMP1[]
                    WORK[I-K+KA-M] = TEMP2[]
                    RA = TEMP3[]
                    # create nonzero element a(i-k,i-k+ka+1) outside the
                    # band and store it in WORK(i-k)
                    t = -BB[KB1-K, I]*RA1
                    WORK[I-K] = WORK[N+I-K+KA-M]*t - WORK[I-K+KA-M]*AB[1, I-K+KA]
                    AB[1, I-K+KA] = WORK[I-K+KA-M]*t + WORK[N+I-K+KA-M]*AB[1, I-K+KA]
                    RA1 = RA
                end
            end
            J2 = I - K - 1 + max(1, K-I0+2)*KA1
            NR = ( N-J2+KA ) ÷ KA1
            J1 = J2 + ( NR-1 )*KA1
            if UPDATE
                J2T = max( J2, I+2*KA-K+1 )
            else
                J2T = J2
            end
            NRT = ( N-J2T+KA ) ÷ KA1
            for J = J2T:KA1:J1
                # create nonzero element a(j-ka,j+1) outside the band
                # and store it in WORK(j-m)
                WORK[J-M] = WORK[J-M]*AB[1, J+1]
                AB[1, J+1] = WORK[N+J-M]*AB[1, J+1]
            end
            # generate rotations in 1st set to annihilate elements which
            # have been created outside the band
            if NRT > 0
                largv!(NRT, pointer(AB, 1+LDAB*(J2T-1)), INCA, pointer(WORK, J2T-M), KA1, pointer(WORK, N+J2T-M), KA1)
            end
            if NR > 0
                # apply rotations in 1st set from the right
                for L = 1:KA-1
                    lartv!(NR, pointer(AB, KA1-L+LDAB*(J2-1)), INCA, pointer(AB, KA-L+LDAB*J2), INCA, pointer(WORK, N+J2-M), pointer(WORK, J2-M), KA1)
                end
                # apply rotations in 1st set from both sides to diagonal
                # blocks
                lar2v!(NR, pointer(AB, KA1+LDAB*(J2-1)), pointer(AB, KA1+LDAB*J2), pointer(AB, KA+LDAB*J2), INCA, pointer(WORK, N+J2-M), pointer(WORK, J2-M), KA1)
            end
            # start applying rotations in 1st set from the left
            for L = KA-1:-1:KB-K+1
                NRT = ( N-J2+L ) ÷ KA1
                if NRT > 0
                    lartv!(NRT, pointer(AB, L+LDAB*(J2+KA1-L-1)), INCA, pointer(AB, L+1+LDAB*(J2+KA1-L-1)), INCA, pointer(WORK, N+J2-M), pointer(WORK, J2-M), KA1)
                end
            end

            if WANTX
                # post-multiply X by product of rotations in 1st set
                for J = J2:KA1:J1
                    push!(X, Givens(J, J+1, WORK[N+J-M], -WORK[J-M]))
                end
            end
        end

        if UPDATE
            if I2 ≤ N && KBT > 0
                # create nonzero element a(i-kbt,i-kbt+ka+1) outside the
                # band and store it in WORK(i-kbt)
                WORK[I-KBT] = -BB[KB1-KBT, I]*RA1
            end
        end

        for K = KB:-1:1
            if UPDATE
                J2 = I - K - 1 + max( 2, K-I0+1 )*KA1
            else
                J2 = I - K - 1 + max( 1, K-I0+1 )*KA1
            end
            # finish applying rotations in 2nd set from the left
            for L = KB-K:-1:1
                NRT = ( N-J2+KA+L ) ÷ KA1
                if NRT > 0
                    lartv!(NRT, pointer(AB, L+LDAB*(J2-L)), INCA, pointer(AB, L+1+LDAB*(J2-L)), INCA, pointer(WORK, N+J2-KA), pointer(WORK, J2-KA), KA1)
                end
            end
            NR = ( N-J2+KA ) ÷ KA1
            J1 = J2 + ( NR-1 )*KA1
            for J = J1:-KA1:J2
                WORK[J] = WORK[J-KA]
                WORK[N+J] = WORK[N+J-KA]
            end
            for J = J2:KA1:J1
                # create nonzero element a(j-ka,j+1) outside the band
                # and store it in WORK(j)
                WORK[J] = WORK[J]*AB[1, J+1]
                AB[1, J+1] = WORK[N+J]*AB[1, J+1]
            end
            if UPDATE
                if I-K < N-KA && K ≤ KBT
                    WORK[I-K+KA] = WORK[I-K]
                end
            end
        end

        for K = KB:-1:1
            J2 = I - K - 1 + max( 1, K-I0+1 )*KA1
            NR = ( N-J2+KA ) ÷ KA1
            J1 = J2 + ( NR-1 )*KA1
            if NR > 0
                # generate rotations in 2nd set to annihilate elements
                # which have been created outside the band
                largv!(NR, pointer(AB, 1+LDAB*(J2-1)), INCA, pointer(WORK, J2), KA1, pointer(WORK, N+J2), KA1)
                # apply rotations in 2nd set from the right
                for L = 1:KA-1
                    lartv!(NR, pointer(AB, KA1-L+LDAB*(J2-1)), INCA, pointer(AB, KA-L+LDAB*J2), INCA, pointer(WORK, N+J2), pointer(WORK, J2), KA1)
                end
                # apply rotations in 2nd set from both sides to diagonal
                # blocks
                lar2v!(NR, pointer(AB, KA1+LDAB*(J2-1)), pointer(AB, KA1+LDAB*J2), pointer(AB, KA+LDAB*J2), INCA, pointer(WORK, N+J2), pointer(WORK, J2), KA1)
            end
            # start applying rotations in 2nd set from the left
            for L = KA-1:-1:KB-K+1
                NRT = ( N-J2+L ) ÷ KA1
                if NRT > 0
                    lartv!(NRT, pointer(AB, L+LDAB*(J2+KA1-L-1)), INCA, pointer(AB, L+1+LDAB*(J2+KA1-L-1)), INCA, pointer(WORK, N+J2), pointer(WORK, J2), KA1)
                end
            end

            if WANTX
                # post-multiply X by product of rotations in 2nd set
                for J = J2:KA1:J1
                    push!(X, Givens(J, J+1, WORK[N+J], -WORK[J]))
                end
            end
        end

        for K = 1:KB-1
            J2 = I - K - 1 + max( 1, K-I0+2 )*KA1
            for L = KB-K:-1:1
                NRT = ( N-J2+L ) ÷ KA1
                if NRT > 0
                    lartv!(NRT, pointer(AB, L+LDAB*(J2+KA1-L-1)), INCA, pointer(AB, L+1+LDAB*(J2+KA1-L-1)), INCA, pointer(WORK, N+J2-M), pointer(WORK, J2-M), KA1)
                end
            end
        end

        if KB > 1
            for J = N-1:-1:I-KB+2*KA+1
                WORK[N+J-M] = WORK[N+J-KA-M]
                WORK[J-M] = WORK[J-KA-M]
            end
        end
    else
        # Transform A, working with the lower triangle
        if UPDATE
            # Form inv(S(i))**T * A * inv(S(i))
            BII = BB[1, I]
            for J = I:I1
                AB[J-I+1, I] = AB[J-I+1, I] / BII
            end
            for J = max(1,I-KA):I
                AB[I-J+1, J] = AB[I-J+1, J] / BII
            end
            for K = I-KBT:I-1
                for J = I-KBT:K
                    AB[K-J+1,J] = AB[K-J+1, J] - BB[I-J+1, J]*AB[I-K+1, K] - BB[I-K+1, K]*AB[I-J+1, J] + AB[1, I]*BB[I-J+1, J]*BB[I-K+1, K]
                end
                for J = max(1, I-KA):I-KBT-1
                    AB[K-J+1, J] = AB[K-J+1, J] - BB[I-K+1, K]*AB[I-J+1, J]
                end
            end
            for J = I:I1
                for K = max(J-KA, I-KBT):I-1
                    AB[J-K+1, K] = AB[J-K+1, K] - BB[I-K+1, K]*AB[J-I+1, I]
                end
            end
            # store a(i1,i) in RA1 for use in next loop over K
            RA1 = AB[I1-I+1, I]
        end
        # Generate and apply vectors of rotations to chase all the
        # existing bulges KA positions down toward the bottom of the
        # band
        for K = 1:KB-1
            if UPDATE
                # Determine the rotations which would annihilate the bulge
                # which has in theory just been created
                if I-K+KA < N && I-K > 1
                    # generate rotation to annihilate a(i-k+ka+1,i)
                    lartg!(AB[KA1-K, I], RA1, TEMP1, TEMP2, TEMP3)
                    WORK[N+I-K+KA-M] = TEMP1[]
                    WORK[I-K+KA-M] = TEMP2[]
                    RA = TEMP3[]
                    t = -BB[K+1, I-K]*RA1
                    WORK[I-K] = WORK[N+I-K+KA-M]*t - WORK[I-K+KA-M]*AB[KA1, I-K]
                    AB[KA1, I-K] = WORK[I-K+KA-M]*t + WORK[N+I-K+KA-M]*AB[KA1, I-K]
                    RA1 = RA
                end
            end
            J2 = I - K - 1 + max( 1, K-I0+2 )*KA1
            NR = ( N-J2+KA ) ÷ KA1
            J1 = J2 + ( NR-1 )*KA1
            if UPDATE
                J2T = max( J2, I+2*KA-K+1 )
            else
                J2T = J2
            end
            NRT = ( N-J2T+KA ) ÷ KA1
            for J = J2T:KA1:J1
                # create nonzero element a(j+1,j-ka) outside the band
                # and store it in WORK(j-m)
                WORK[J-M] = WORK[J-M]*AB[KA1, J-KA+1]
                AB[KA1, J-KA+1] = WORK[N+J-M]*AB[KA1, J-KA+1]
            end
            # generate rotations in 1st set to annihilate elements which
            # have been created outside the band
            if NRT > 0
                largv!(NRT, pointer(AB, KA1+LDAB*(J2T-KA-1)), INCA, pointer(WORK, J2T-M), KA1, pointer(WORK, N+J2T-M), KA1)
            end
            if NR > 0
                # apply rotations in 1st set from the left
                for L = 1:KA-1
                    lartv!(NR, pointer(AB, L+1+LDAB*(J2-L-1)), INCA, pointer(AB, L+2+LDAB*(J2-L-1)), INCA, pointer(WORK, N+J2-M), pointer(WORK, J2-M), KA1)
                end
                # apply rotations in 1st set from both sides to diagonal
                # blocks
                lar2v!(NR, pointer(AB, 1+LDAB*(J2-1)), pointer(AB, 1+LDAB*J2), pointer(AB, 2+LDAB*(J2-1)), INCA, pointer(WORK, N+J2-M), pointer(WORK, J2-M), KA1)
            end
            # start applying rotations in 1st set from the right
            for L = KA-1:-1:KB-K+1
                NRT = ( N-J2+L ) ÷ KA1
                if NRT > 0
                    lartv!(NRT, pointer(AB, KA1-L+1+LDAB*(J2-1)), INCA, pointer(AB, KA1-L+LDAB*J2), INCA, pointer(WORK, N+J2-M), pointer(WORK, J2-M), KA1)
                end
            end

            if WANTX
                # post-multiply X by product of rotations in 1st set
                for J = J2:KA1:J1
                    push!(X, Givens(J, J+1, WORK[N+J-M], -WORK[J-M]))
                end
            end
        end

        if UPDATE
            if I2 ≤ N && KBT > 0
                # create nonzero element a(i-kbt+ka+1,i-kbt) outside the
                # band and store it in WORK(i-kbt)
                WORK[I-KBT] = -BB[KBT+1, I-KBT]*RA1
            end
        end

        for K = KB:-1:1
            if UPDATE
                J2 = I - K - 1 + max( 2, K-I0+1 )*KA1
            else
                J2 = I - K - 1 + max( 1, K-I0+1 )*KA1
            end
            # finish applying rotations in 2nd set from the right
            for L = KB-K:-1:1
                NRT = ( N-J2+KA+L ) ÷ KA1
                if NRT > 0
                    lartv!(NRT, pointer(AB, KA1-L+1+LDAB*(J2-KA-1)), INCA, pointer(AB, KA1-L+LDAB*(J2-KA)), INCA, pointer(WORK, N+J2-KA), pointer(WORK, J2-KA), KA1)
                end
            end
            NR = ( N-J2+KA ) ÷ KA1
            J1 = J2 + ( NR-1 )*KA1
            for J = J1:-KA1:J2
                WORK[J] = WORK[J-KA]
                WORK[N+J] = WORK[N+J-KA]
            end
            for J = J2:KA1:J1
                # create nonzero element a(j+1,j-ka) outside the band
                # and store it in WORK(j)
                WORK[J] = WORK[J]*AB[KA1, J-KA+1]
                AB[KA1, J-KA+1] = WORK[N+J]*AB[KA1, J-KA+1]
            end
            if UPDATE
                if I-K < N-KA && K ≤ KBT
                    WORK[I-K+KA] = WORK[I-K]
                end
            end
        end

        for K = KB:-1:1
            J2 = I - K - 1 + max( 1, K-I0+1 )*KA1
            NR = ( N-J2+KA ) ÷ KA1
            J1 = J2 + ( NR-1 )*KA1
            if NR > 0
                # generate rotations in 2nd set to annihilate elements
                # which have been created outside the band
                largv!(NR, pointer(AB, KA1+LDAB*(J2-KA-1)), INCA, pointer(WORK, J2), KA1, pointer(WORK, N+J2), KA1)
                # apply rotations in 2nd set from the left
                for L = 1:KA-1
                    lartv!(NR, pointer(AB, L+1+LDAB*(J2-L-1)), INCA, pointer(AB, L+2+LDAB*(J2-L-1)), INCA, pointer(WORK, N+J2), pointer(WORK, J2), KA1)
                end
                # apply rotations in 2nd set from both sides to diagonal
                # blocks
                lar2v!(NR, pointer(AB, 1+LDAB*(J2-1)), pointer(AB, 1+LDAB*J2), pointer(AB, 2+LDAB*(J2-1)), INCA, pointer(WORK, N+J2), pointer(WORK, J2), KA1)
            end
            # start applying rotations in 2nd set from the right
            for L = KA-1:-1:KB-K+1
                NRT = ( N-J2+L ) ÷ KA1
                if NRT > 0
                    lartv!(NRT, pointer(AB, KA1-L+1+LDAB*(J2-1)), INCA, pointer(AB, KA1-L+LDAB*J2), INCA, pointer(WORK, N+J2), pointer(WORK, J2), KA1)
                end
            end

            if WANTX
                # post-multiply X by product of rotations in 2nd set
                for J = J2:KA1:J1
                    push!(X, Givens(J, J+1, WORK[N+J], -WORK[J]))
                end
            end
        end

        for K = 1:KB-1
            J2 = I - K - 1 + max( 1, K-I0+2 )*KA1
            # finish applying rotations in 1st set from the right
            for L = KB-K:-1:1
                NRT = ( N-J2+L ) ÷ KA1
                if NRT > 0
                    lartv!(NRT, pointer(AB, KA1-L+1+LDAB*(J2-1)), INCA, pointer(AB, KA1-L+LDAB*J2), INCA, pointer(WORK, N+J2-M), pointer(WORK, J2-M), KA1)
                end
            end
        end

        if KB > 1
            for J = N-1:-1:I-KB+2*KA+1
                WORK[N+J-M] = WORK[N+J-KA-M]
                WORK[J-M] = WORK[J-KA-M]
            end
        end
    end
    @goto L10
    @label L480

    UPDATE = true
    I = 0
    @label L490
    if UPDATE
        I = I + 1
        KBT = min( KB, M-I )
        I0 = I + 1
        I1 = max( 1, I-KA )
        I2 = I + KBT - KA1
        if I > M
            UPDATE = false
            I = I - 1
            I0 = M + 1
            if KA == 0
                return
            end
            @goto L490
        end
    else
        I = I - KA
        if I < 2
            return
        end
    end

    if I < M - KBT
        NX = M
    else
        NX = N
    end

    if UPPER
        # Transform A, working with the upper triangle
        if UPDATE
            # Form inv(S(i))**T * A * inv(S(i))
            BII = BB[KB1, I]
            for J = I1:I
                AB[J-I+KA1, I] = AB[J-I+KA1, I] / BII
            end
            for J = I:min(N, I+KA)
                AB[I-J+KA1, J] = AB[I-J+KA1, J] / BII
            end
            for K = I+1:I+KBT
                for J = K:I+KBT
                    AB[K-J+KA1, J] = AB[K-J+KA1, J] - BB[I-J+KB1, J]*AB[I-K+KA1, K] - BB[I-K+KB1, K]*AB[I-J+KA1, J] + AB[KA1, I]*BB[I-J+KB1, J]*BB[I-K+KB1, K]
                end
                for J = I+KBT+1:min(N, I+KA)
                    AB[K-J+KA1, J] = AB[K-J+KA1, J] - BB[I-K+KB1, K]*AB[I-J+KA1, J]
                end
            end
            for J = I1:I
                for K = I+1:min(J+KA, I+KBT)
                    AB[J-K+KA1, K] = AB[J-K+KA1, K] - BB[I-K+KB1, K]*AB[J-I+KA1, I]
                end
            end
            # store a(i1,i) in RA1 for use in next loop over K
            RA1 = AB[I1-I+KA1, I]
        end

        for K = 1:KB-1
            if UPDATE
                # Determine the rotations which would annihilate the bulge
                # which has in theory just been created
                if I+K-KA1 > 0 && I+K < M
                    # generate rotation to annihilate a(i+k-ka-1,i)
                    lartg!(AB[K+1, I], RA1, TEMP1, TEMP2, TEMP3)
                    WORK[N+I+K-KA] = TEMP1[]
                    WORK[I+K-KA] = TEMP2[]
                    RA = TEMP3[]
                    # create nonzero element a(i+k-ka-1,i+k) outside the
                    # band and store it in WORK(m-kb+i+k)
                    t = -BB[KB1-K, I+K]*RA1
                    WORK[M-KB+I+K] = WORK[N+I+K-KA]*t - WORK[I+K-KA]*AB[1, I+K]
                    AB[1, I+K] = WORK[I+K-KA]*t + WORK[N+I+K-KA]*AB[1, I+K]
                    RA1 = RA
                end
            end
            J2 = I + K + 1 - max( 1, K+I0-M+1 )*KA1
            NR = ( J2+KA-1 ) ÷ KA1
            J1 = J2 - ( NR-1 )*KA1
            if UPDATE
                J2T = min(J2, I-2*KA+K-1)
            else
                J2T = J2
            end
            NRT = ( J2T+KA-1 ) ÷ KA1
            for J = J1:KA1:J2T
                # create nonzero element a(j-1,j+ka) outside the band
                # and store it in WORK(j)
                WORK[J] = WORK[J]*AB[1, J+KA-1]
                AB[1, J+KA-1] = WORK[N+J]*AB[1, J+KA-1]
            end
            # generate rotations in 1st set to annihilate elements which
            # have been created outside the band

            if NRT > 0
                largv!(NRT, pointer(AB, 1+LDAB*(J1+KA-1)), INCA, pointer(WORK, J1), KA1, pointer(WORK, N+J1), KA1)
            end
            if NR > 0
                # apply rotations in 1st set from the left
                for L = 1:KA-1
                    lartv!(NR, pointer(AB, KA1-L+LDAB*(J1+L-1)), INCA, pointer(AB, KA-L+LDAB*(J1+L-1)), INCA, pointer(WORK, N+J1), pointer(WORK, J1), KA1)
                end
                # apply rotations in 1st set from both sides to diagonal
                # blocks
                lar2v!(NR, pointer(AB, KA1+LDAB*(J1-1)), pointer(AB, KA1+LDAB*(J1-2)), pointer(AB, KA+LDAB*(J1-1)), INCA, pointer(WORK, N+J1), pointer(WORK, J1), KA1)
            end
            # start applying rotations in 1st set from the right
            for L = KA-1:-1:KB-K+1
                NRT = ( J2+L-1 ) ÷ KA1
                J1T = J2 - ( NRT-1 )*KA1
                if NRT > 0
                    lartv!(NRT, pointer(AB, L+LDAB*(J1T-1)), INCA, pointer(AB, L+1+LDAB*(J1T-2)), INCA, pointer(WORK, N+J1T), pointer(WORK, J1T), KA1)
                end
            end

            if WANTX
                # post-multiply X by product of rotations in 1st set
                for J = J1:KA1:J2
                    push!(X, Givens(J, J-1, WORK[N+J], -WORK[J]))
                end
            end
        end

        if UPDATE
            if I2 > 0 && KBT > 0
                # create nonzero element a(i+kbt-ka-1,i+kbt) outside the
                # band and store it in WORK(m-kb+i+kbt)
                WORK[M-KB+I+KBT] = -BB[KB1-KBT, I+KBT]*RA1
            end
        end

        for K = KB:-1:1
            if UPDATE
                J2 = I + K + 1 - max( 2, K+I0-M )*KA1
            else
                J2 = I + K + 1 - max( 1, K+I0-M )*KA1
            end
            # finish applying rotations in 2nd set from the right
            for L = KB-K:-1:1
                NRT = ( J2+KA+L-1 ) ÷ KA1
                J1T = J2 - ( NRT-1 )*KA1
                if NRT > 0
                    lartv!(NRT, pointer(AB, L+LDAB*(J1T+KA-1)), INCA, pointer(AB, L+1+LDAB*(J1T+KA-2)), INCA, pointer(WORK, N+M-KB+J1T+KA), pointer(WORK, M-KB+J1T+KA), KA1)
                end
            end
            NR = ( J2+KA-1 ) ÷ KA1
            J1 = J2 - ( NR-1 )*KA1
            for J = J1:KA1:J2
                WORK[M-KB+J] = WORK[M-KB+J+KA]
                WORK[N+M-KB+J] = WORK[N+M-KB+J+KA]
            end
            for J = J1:KA1:J2
                # create nonzero element a(j-1,j+ka) outside the band
                # and store it in WORK(m-kb+j)
                WORK[M-KB+J] = WORK[M-KB+J]*AB[1, J+KA-1]
                AB[1, J+KA-1] = WORK[N+M-KB+J]*AB[1, J+KA-1]
            end
            if UPDATE
                if I+K > KA1 && K ≤ KBT
                    WORK[M-KB+I+K-KA] = WORK[M-KB+I+K]
                end
            end
        end

        for K = KB:-1:1
            J2 = I + K + 1 - max( 1, K+I0-M )*KA1
            NR = ( J2+KA-1 ) ÷ KA1
            J1 = J2 - ( NR-1 )*KA1
            if NR > 0
                # generate rotations in 2nd set to annihilate elements
                # which have been created outside the band
                largv!(NR, pointer(AB, 1+LDAB*(J1+KA-1)), INCA, pointer(WORK, M-KB+J1), KA1, pointer(WORK, N+M-KB+J1), KA1)
                # apply rotations in 2nd set from the left
                for L = 1:KA-1
                    lartv!(NR, pointer(AB, KA1-L+LDAB*(J1+L-1)), INCA, pointer(AB, KA-L+LDAB*(J1+L-1)), INCA, pointer(WORK, N+M-KB+J1), pointer(WORK, M-KB+J1), KA1)
                end
                # apply rotations in 2nd set from both sides to diagonal
                # blocks
                lar2v!(NR, pointer(AB, KA1+LDAB*(J1-1)), pointer(AB, KA1+LDAB*(J1-2)), pointer(AB, KA+LDAB*(J1-1)), INCA, pointer(WORK, N+M-KB+J1), pointer(WORK, M-KB+J1), KA1)
            end
            # start applying rotations in 2nd set from the right
            for L = KA-1:-1:KB-K+1
                NRT = ( J2+L-1 ) ÷ KA1
                J1T = J2 - ( NRT-1 )*KA1
                if NRT > 0
                    lartv!(NRT, pointer(AB, L+LDAB*(J1T-1)), INCA, pointer(AB, L+1+LDAB*(J1T-2)), INCA, pointer(WORK, N+M-KB+J1T), pointer(WORK, M-KB+J1T), KA1)
                end
            end

            if WANTX
                # post-multiply X by product of rotations in 2nd set
                for J = J1:KA1:J2
                    push!(X, Givens(J, J-1, WORK[N+M-KB+J], -WORK[M-KB+J]))
                end
            end
        end

        for K = 1:KB-1
            J2 = I + K + 1 - max( 1, K+I0-M+1 )*KA1
            # finish applying rotations in 1st set from the right
            for L = KB-K:-1:1
                NRT = ( J2+L-1 ) ÷ KA1
                J1T = J2 - ( NRT-1 )*KA1
                if NRT > 0
                    lartv!(NRT, pointer(AB, L+LDAB*(J1T-1)), INCA, pointer(AB, L+1+LDAB*(J1T-2)), INCA, pointer(WORK, N+J1T), pointer(WORK, J1T), KA1)
                end
            end
        end

        if KB > 1
            for J = 2:min(I+KB, M)-2*KA-1
                WORK[N+J] = WORK[N+J+KA]
                WORK[J] = WORK[J+KA]
            end
        end
    else
        # Transform A, working with the lower triangle
        if UPDATE
            # Form inv(S(i))**T * A * inv(S(i))
            BII = BB[1, I]
            for J = I1:I
                AB[I-J+1, J] = AB[I-J+1, J] / BII
            end
            for J = I:min(N, I+KA)
                AB[J-I+1, I] = AB[J-I+1, I] / BII
            end
            for K = I+1:I+KBT
                for J = K:I+KBT
                    AB[J-K+1, K] = AB[J-K+1, K] - BB[J-I+1, I]*AB[K-I+1, I] - BB[K-I+1, I]*AB[J-I+1, I] + AB[1, I]*BB[J-I+1, I]*BB[K-I+1, I]
                end
                for J = I+KBT+1:min(N,I+KA)
                    AB[J-K+1, K] = AB[J-K+1, K] - BB[K-I+1, I]*AB[J-I+1, I]
                end
            end
            for J = I1:I
                for K = I+1:min(J+KA, I+KBT)
                    AB[K-J+1, J] = AB[K-J+1, J] - BB[K-I+1, I]*AB[I-J+1, J]
                end
            end
            # store a(i,i1) in RA1 for use in next loop over K
            RA1 = AB[I-I1+1, I1]
        end
        # Generate and apply vectors of rotations to chase all the
        # existing bulges KA positions up toward the top of the band
        for K = 1:KB-1
            if UPDATE
                # Determine the rotations which would annihilate the bulge
                # which has in theory just been created
                if I+K-KA1 > 0 && I+K < M
                    # generate rotation to annihilate a(i,i+k-ka-1)
                    lartg!(AB[KA1-K, I+K-KA], RA1, TEMP1, TEMP2, TEMP3)
                    WORK[N+I+K-KA] = TEMP1[]
                    WORK[I+K-KA] = TEMP2[]
                    RA = TEMP3[]
                    # create nonzero element a(i+k,i+k-ka-1) outside the
                    # band and store it in WORK(m-kb+i+k)
                    t = -BB[K+1, I]*RA1
                    WORK[M-KB+I+K] = WORK[N+I+K-KA]*t - WORK[I+K-KA]*AB[KA1, I+K-KA]
                    AB[KA1, I+K-KA] = WORK[I+K-KA]*t + WORK[N+I+K-KA]*AB[KA1, I+K-KA]
                    RA1 = RA
                end
            end
            J2 = I + K + 1 - max( 1, K+I0-M+1 )*KA1
            NR = ( J2+KA-1 ) ÷ KA1
            J1 = J2 - ( NR-1 )*KA1
            if UPDATE
                J2T = min( J2, I-2*KA+K-1 )
            else
                J2T = J2
            end
            NRT = ( J2T+KA-1 ) ÷ KA1
            for J = J1:KA1:J2T
                # create nonzero element a(j+ka,j-1) outside the band
                # and store it in WORK(j)
                WORK[J] = WORK[J]*AB[KA1, J-1]
                AB[KA1, J-1] = WORK[N+J]*AB[KA1, J-1]
            end
            # generate rotations in 1st set to annihilate elements which
            # have been created outside the band
            if NRT > 0
                largv!(NRT, pointer(AB, KA1+LDAB*(J1-1)), INCA, pointer(WORK, J1), KA1, pointer(WORK, N+J1), KA1)
            end
            if NR > 0
                # apply rotations in 1st set from the right
                for L = 1:KA-1
                    lartv!(NR, pointer(AB, L+1+LDAB*(J1-1)), INCA, pointer(AB, L+2+LDAB*(J1-2)), INCA, pointer(WORK, N+J1), pointer(WORK, J1), KA1)
                end
                # apply rotations in 1st set from both sides to diagonal
                # blocks
                lar2v!(NR, pointer(AB, 1+LDAB*(J1-1)), pointer(AB, 1+LDAB*(J1-2)), pointer(AB, 2+LDAB*(J1-2)), INCA, pointer(WORK, N+J1), pointer(WORK, J1), KA1)
            end
            # start applying rotations in 1st set from the left
            for L = KA-1:-1:KB-K+1
                NRT = ( J2+L-1 ) ÷ KA1
                J1T = J2 - ( NRT-1 )*KA1
                if NRT > 0
                    lartv!(NRT, pointer(AB, KA1-L+1+LDAB*(J1T-KA1+L-1)), INCA, pointer(AB, KA1-L+LDAB*(J1T-KA1+L-1)), INCA, pointer(WORK, N+J1T), pointer(WORK, J1T), KA1)
                end
            end

            if WANTX
                # post-multiply X by product of rotations in 1st set
                for J = J1:KA1:J2
                    push!(X, Givens(J, J-1, WORK[N+J], -WORK[J]))
                end
            end
        end

        if UPDATE
            if I2 > 0 && KBT > 0
                # create nonzero element a(i+kbt,i+kbt-ka-1) outside the
                # band and store it in WORK(m-kb+i+kbt)
                WORK[M-KB+I+KBT] = -BB[KBT+1, I]*RA1
            end
        end

        for K = KB:-1:1
            if UPDATE
                J2 = I + K + 1 - max( 2, K+I0-M )*KA1
            else
                J2 = I + K + 1 - max( 1, K+I0-M )*KA1
            end
            # finish applying rotations in 2nd set from the left
            for L = KB-K:-1:1
                NRT = ( J2+KA+L-1 ) ÷ KA1
                J1T = J2 - ( NRT-1 )*KA1
                if NRT > 0
                    lartv!(NRT, pointer(AB, KA1-L+1+LDAB*(J1T+L-2)), INCA, pointer(AB, KA1-L+LDAB*(J1T+L-2)), INCA, pointer(WORK, N+M-KB+J1T+KA), pointer(WORK, M-KB+J1T+KA), KA1)
                end
            end
            NR = ( J2+KA-1 ) ÷ KA1
            J1 = J2 - ( NR-1 )*KA1
            for J = J1:KA1:J2
                WORK[M-KB+J] = WORK[M-KB+J+KA]
                WORK[N+M-KB+J] = WORK[N+M-KB+J+KA]
            end
            for J = J1:KA1:J2
                # create nonzero element a(j+ka,j-1) outside the band
                # and store it in WORK(m-kb+j)
                WORK[M-KB+J] = WORK[M-KB+J]*AB[KA1, J-1]
                AB[KA1, J-1] = WORK[N+M-KB+J]*AB[KA1, J-1]
            end
            if UPDATE
                if I+K > KA1 && K ≤ KBT
                    WORK[M-KB+I+K-KA] = WORK[M-KB+I+K]
                end
            end
        end

        for K = KB:-1:1
            J2 = I + K + 1 - max( 1, K+I0-M )*KA1
            NR = ( J2+KA-1 ) ÷ KA1
            J1 = J2 - ( NR-1 )*KA1
            if NR > 0
                # generate rotations in 2nd set to annihilate elements
                # which have been created outside the band
                largv!(NR, pointer(AB, KA1+LDAB*(J1-1)), INCA, pointer(WORK, M-KB+J1), KA1, pointer(WORK, N+M-KB+J1), KA1)
                # apply rotations in 2nd set from the right
                for L = 1:KA-1
                    lartv!(NR, pointer(AB, L+1+LDAB*(J1-1)), INCA, pointer(AB, L+2+LDAB*(J1-2)), INCA, pointer(WORK, N+M-KB+J1), pointer(WORK, M-KB+J1), KA1)
                end
                # apply rotations in 2nd set from both sides to diagonal
                # blocks
                lar2v!(NR, pointer(AB, 1+LDAB*(J1-1)), pointer(AB, 1+LDAB*(J1-2)), pointer(AB, 2+LDAB*(J1-2)), INCA, pointer(WORK, N+M-KB+J1), pointer(WORK, M-KB+J1), KA1)
            end
            # start applying rotations in 2nd set from the left
            for L = KA-1:-1:KB-K+1
                NRT = ( J2+L-1 ) ÷ KA1
                J1T = J2 - ( NRT-1 )*KA1
                if NRT > 0
                    lartv!(NRT, pointer(AB, KA1-L+1+LDAB*(J1T-KA1+L-1)), INCA, pointer(AB, KA1-L+LDAB*(J1T-KA1+L-1)), INCA, pointer(WORK, N+M-KB+J1T), pointer(WORK, M-KB+J1T), KA1)
                end
            end

            if WANTX
                # post-multiply X by product of rotations in 2nd set
                for J = J1:KA1:J2
                    push!(X, Givens(J, J-1, WORK[N+M-KB+J], -WORK[M-KB+J]))
                end
            end
        end

        for K = 1:KB-1
            J2 = I + K + 1 - max( 1, K+I0-M+1 )*KA1
            # finish applying rotations in 1st set from the left
            for L = KB-K:-1:1
                NRT = ( J2+L-1 ) ÷ KA1
                J1T = J2 - ( NRT-1 )*KA1
                if NRT > 0
                    lartv!(NRT, pointer(AB, KA1-L+1+LDAB*(J1T-KA1+L-1)), INCA, pointer(AB, KA1-L+LDAB*(J1T-KA1+L-1)), INCA, pointer(WORK, N+J1T), pointer(WORK, J1T), KA1)
                end
            end
        end

        if KB > 1
            for J = 2:min(I+KB,M)-2*KA-1
                WORK[N+J] = WORK[N+J+KA]
                WORK[J] = WORK[J+KA]
            end
        end
    end
    @goto L490
end


text:
module FastDirectedRounding

import Base: (+), (-), (*), (/), (\), sqrt, hypot, inv

using AdjacentFloat

import ErrorfreeArithmetic: eftAdd, eftSub, eftMul, eftSqr,
    accInv, accDiv, accSqrt, accInvForSign, accDivForSign, accSqrtForSign

export sqr

include("eftRound.jl")
include("directRounding.jl")

end # module

ext: jl


text:
module Lib

using CxxWrap

const SRCDIR = @__DIR__
const PKGDIR = dirname(SRCDIR)
const DEPSDIR = joinpath(PKGDIR, "deps")

@wrapmodule(joinpath(DEPSDIR, "lib", "libpcm.so"))

function __init__()
    @initcxx
end

end


text:
# This file contains functions for loading events data

include("downloads.jl")
include("loader.jl")
include("competitions.jl")
include("games.jl")
include("teams.jl")
include("players.jl")
include("utils.jl")
include("events.jl")

"""
    function get_events_data()
Obtaining the events stream data.

# Optional Arguments
- `download=true`: specifies if the data need to be downloaded or not. If set to false, a `path` need to be provided.
- `path=""`; specifies the local path where the data are located if `download=false` and the path where the data will be stored if `download=true`.
"""
function get_events_data(source::Symbol; download::Bool = true, path::String = "")

    if download
        path = joinpath(path, "wyscout_data")
        mkdir(path)
    elseif !download && path != ""
        path = joinpath(path)
    else
        throw(ArgumentError("Please provide a local path"));
    end

    # At the moment we only support Wyscout data
    if source == :wyscout
        event_data_urls = PublicWyscoutURLs()
        if download || length(readdir(path)) == 0
            download_repo(event_data_urls, path)
        end
        loader = PublicWyscoutLoader(path = path)
        loader = create_match_index(loader)
    else
        throw(ArgumentError("Please provide a valid source name: `:wyscout`"));
    end

    return loader
end


text:
using TypeParams
using Test
using MacroTools: gensym_ids

@testset "@typeparams" begin
    @test @macroexpand(@typeparams struct T; end) == :(struct T{}; end)
    @test @macroexpand(@typeparams mutable struct T; end) == :(mutable struct T{}; end)
    @test @macroexpand(@typeparams struct T; a; end) == :(struct T{}; a; end)
    @test @macroexpand(@typeparams struct T; a::Int; end) == :(struct T{}; a::Int; end)
    @test @macroexpand(@typeparams struct T{A}; a::A; end) == :(struct T{A}; a::A; end)
    @test @macroexpand(@typeparams struct T; a; T(a) = new(a); end) == :(struct T{}; a; T(a) = new(a); end)

    @test @macroexpand(@typeparams struct T; a::{A}; end) == :(struct T{A}; a::A; end)
    @test @macroexpand(@typeparams struct T; a::{A<:Integer}; end) == :(struct T{A<:Integer}; a::A; end)
    @test @macroexpand(@typeparams struct T{A,B}; a::A; b::{B}; end) == :(struct T{A,B}; a::A; b::B; end)
    @test @macroexpand(@typeparams struct T; a1::{A}; a2::{A}; end) == :(struct T{A}; a1::A; a2::A; end)

    @test gensym_ids(@macroexpand(@typeparams struct T; a::{}; end)) == :(struct T{a_1}; a::a_1; end)
    @test gensym_ids(@macroexpand(@typeparams struct T; a::{<:Integer}; end)) == :(struct T{a_1<:Integer}; a::a_1; end)
    @test gensym_ids(@macroexpand(@typeparams struct T; a::{}; b::{}; end)) == :(struct T{a_1,b_2}; a::a_1; b::b_2; end)
end

Base.@kwdef @typeparams struct KWDef
    a::{} = 1
end
@testset "@kwdef interop" begin
    @test KWDef() == KWDef{Int}(1)
end


text:
############################################################################################
"""
$(TYPEDEF)
Default Configuration for NUTS sampler.

# Fields
$(TYPEDFIELDS)
"""
struct ConfigNUTS{K<:KineticEnergy} <: AbstractConfiguration
    ## Global targets
    "Target Acceptance Rate"
    δ::Float64
    "Default size for tuning iterations in each cycle."
    window::ConfigTuningWindow
    "Kinetic Energy used in Hamiltonian: GaussianKineticEnergy"
    energy::K
    function ConfigNUTS(
        δ::Float64,
        window::ConfigTuningWindow,
        energy::K,
    ) where {K<:KineticEnergy}
        @argcheck 0.0 < δ <= 1.0 "Acceptance rate not bounded between 0 and 1"
        return new{K}(
            δ,
            window,
            energy
        )
    end
end

"""
$(SIGNATURES)
Initialize NUTS custom configurations.

# Examples
```julia
```

"""
function init(
    ::Type{AbstractConfiguration},
    mcmc::Type{NUTS},
    objective::Objective,
    proposalconfig::ConfigProposal;
    ## Target Acceptance Rate
    δ=0.80,
    ## MCMC Phase tune variables based on standard HMC best practice target acceptance rate of 0.80
    window= ConfigTuningWindow(
        [1, 5, 1],
        [75, 50, 25],
        [Warmup(), Adaptionˢˡᵒʷ(), Adaptionᶠᵃˢᵗ(), Exploration()]
        ),
    ## Kinetic energy in Hamiltonian
    energy=GaussianKineticEnergy(
        init(
            objective.model.info.reconstruct.default.output, proposalconfig.metric, length(objective.tagged)
        )...,
    )
)
    return ConfigNUTS(
        δ,
        window,
        energy
    )
end

############################################################################################
#export
export ConfigNUTS


text:
# Various heap implementation

###########################################################
#
#   Heap interface specification
#
#   Each heap is associated with a handle type (H), and
#   a value type v.
#
#   Here, the value type must be comparable, and a handle
#   is an object through which one can refer to a specific
#   node of the heap and thus update its value.
#
#   Each heap type must implement all of the following
#   functions. Here, let h be a heap, i be a handle, and
#   v be a value.
#
#   - length(h)         returns the number of elements
#
#   - isempty(h)        returns whether the heap is
#                       empty
#
#   - push!(h, v)       add a value to the heap
#
#   - sizehint!(h)       set size hint to a heap
#
#   - top(h)            return the top value of a heap
#
#   - pop!(h)           removes the top value, and
#                       returns it
#
#  For mutable heaps, it should also support
#
#   - push!(h, v)       adds a value to the heap and
#                       returns a handle to v
#
#   - update!(h, i, v)  updates the value of an element
#                       (referred to by the handle i)
#
###########################################################

# HT: handle type
# VT: value type

@compat abstract type AbstractHeap{VT} end

@compat abstract type AbstractMutableHeap{VT,HT} <: AbstractHeap{VT} end

# comparer

immutable LessThan
end

immutable GreaterThan
end

compare(c::LessThan, x, y) = x < y
compare(c::GreaterThan, x, y) = x > y

# heap implementations

include("heaps/binary_heap.jl")
include("heaps/mutable_binary_heap.jl")
include("heaps/arrays_as_heaps.jl")

# generic functions

function extract_all!{VT}(h::AbstractHeap{VT})
    n = length(h)
    r = Vector{VT}(n)
    for i = 1 : n
        r[i] = pop!(h)
    end
    r
end

function extract_all_rev!{VT}(h::AbstractHeap{VT})
    n = length(h)
    r = Vector{VT}(n)
    for i = 1 : n
        r[n + 1 - i] = pop!(h)
    end
    r
end

# Array functions using heaps

function nextreme{T, Comp}(comp::Comp, n::Int, arr::AbstractVector{T})
    if n <= 0
        return T[] # sort(arr)[1:n] returns [] for n <= 0
    elseif n >= length(arr)
        return sort(arr, lt = (x, y) -> compare(comp, y, x))
    end

    buffer = BinaryHeap{T,Comp}(comp)

    for i = 1 : n
        @inbounds xi = arr[i]
        push!(buffer, xi)
    end

    for i = n + 1 : length(arr)
        @inbounds xi = arr[i]
        if compare(comp, top(buffer), xi)
            # This could use a pushpop method
            pop!(buffer)
            push!(buffer, xi)
        end
    end

    return extract_all_rev!(buffer)
end

@doc """
Returns the `n` largest elements of `arr`.

Equivalent to `sort(arr, lt = >)[1:min(n, end)]`
""" ->
function nlargest{T}(n::Int, arr::AbstractVector{T})
    return nextreme(LessThan(), n, arr)
end

@doc """
Returns the `n` smallest elements of `arr`.

Equivalent to `sort(arr, lt = <)[1:min(n, end)]`
""" ->
function nsmallest{T}(n::Int, arr::AbstractVector{T})
    return nextreme(GreaterThan(), n, arr)
end


text:
using GraphPlot, SimpleHypergraphs
import LightGraphs

using DataFrames
using CSV

#cd(raw"copy_the_path_here_if_your_Julia_cannot_find_files")

b = CSV.read("business.csv",delim="\t")

b.business_id = Symbol.(b.business_id);
for name in names(b)[2:15]
    b[!,name] = Symbol.(b[!, name])
end

revs = CSV.read("review2.csv",delim="\t")

revs.business_id = Symbol.(revs.business_id)

@assert sort(b.business_id) == sort(unique(revs.business_id))

using StatsBase
using Random

using SimpleHypergraphs

function report(revs::DataFrame, b::DataFrame, stars)
    revsg = revs[revs.stars .== stars, :]
    sort!(revsg,:business_id)

    business_ids_set = Set{Symbol}(revsg.business_id)
    #list of businesses filtered to the data
    bf = filter(row -> row[:business_id] in business_ids_set, b)
    sort!(bf,:business_id)
    bf.b_id = 1:nrow(bf)

    revsg.user_id = categorical(revsg.user_id,true)

    # Verices x hyperedges
    h = Hypergraph{Float64}(nrow(bf),length(levels(revsg.user_id)));
    for r in 1:nrow(revsg)
        bix = searchsorted(bf.business_id,revsg.business_id[r]).start
        h[bix,revsg.user_id[r].level] = 1
    end

    println("Constructed a hypergraph having size: " ,size(h))

    for gt in [:city, :state, :Alcohol,:NoiseLevel,:RestaurantsTakeOut,:category]
        partition = Vector{Set{Int}}()
        for bfg in groupby(bf,gt)
            push!(partition,Set{Int}(sort(bfg.b_id)))
        end
        println(gt," ",modularity(h,partition))
    end
end

for stars in 1:5
    println(stars)
    report(revs, b, stars)
end


text:
# Test linear operations on SAC traces
using Test
using SAC
if VERSION >= v"0.7-"
    import FFTW: rfft, irfft
end

function inplace_compare(s, f, f!, val)
    s = deepcopy(s)
    s′ = deepcopy(s)
    f!(s, val)
    f(s′, val) == s
end

function inplace_raw_compare(a::AbstractVector, b, delta, f!, val)
    s = SACtr(a, delta)
    f!(s, val)
    all(s.t .== b)
end

function inplace_raw_compare(a::AbstractVector, b, delta, f)
    s = SACtr(a, delta)
    modify!(s, f)
    all(s.t .== b)
end

@testset "Operations" begin
    ## In-place versions
    # Compare copying and inplace operations
    let s = SAC.sample(), val = rand()
        @test inplace_compare(s, add, add!, val)
        @test inplace_compare(s, subtract, subtract!, val)
        @test inplace_compare(s, divide, divide!, val)
        @test inplace_compare(s, multiply, multiply!, val)
    end
    # Compare inplace operations with raw array operations
    let a = rand(SAC.SACFloat, 100), val = rand(SAC.SACFloat)-SAC.SACFloat(2), delta = 0.01
        @test inplace_raw_compare(a, a.+val, delta, add!, val)
        @test inplace_raw_compare(a, a.-val, delta, add!, -val)
        @test inplace_raw_compare(a, a.*val, delta, multiply!, val)
        @test inplace_raw_compare(a, a./val, delta, divide!, val)
    end
    # General modification
    let a = rand(SAC.SACFloat, 100), delta = 2
        @test inplace_raw_compare(a, a.^2, delta, x->x^2)
        @test inplace_raw_compare(a, rem.(a, SAC.SACFloat(0.1)), delta, x->rem(x, SAC.SACFloat(0.1)))
        @test inplace_raw_compare(a, 2 .* (a .- 1)/20, delta, x->2*(x-1)/20)
    end

    ## Copying versions
    # Trace and constant
    let s = SAC.sample(), val = rand(SAC.SACFloat)
        # Addition/subtraction: raw values
        @test typeof(s + val) == SACtr
        @test all(s.t .+ val .== (s + val).t)
        @test all(s.t .+ val .== (val + s).t)
        @test all(val .+ s.t .== (val + s).t)
        @test all(val .+ s.t .== (s + val).t)
        @test s + 0 == s
        @test typeof(s - val) == SACtr
        @test all(s.t .- val .== (s - val).t)
        @test all(s.t .- val .== (-val + s).t)
        @test all(val .- s.t .== (val - s).t)
        @test all(val .- s.t .== (-s + val).t)
        @test s - 0 == s
        # Order
        @test s + val == val + s
        @test s - val == -val + s
        # Multiplication/divison: raw values
        val += SAC.SACFloat(1) # Ensure we're not dividing by zero
        @test all(s.t .* val .== (s*val).t)
        @test all(s.t .* val .== (val*s).t)
        @test s*1 == s
        @test all(s.t ./ val .== (s/val).t)
        @test all(val ./ s.t .== (val/s).t)
        @test s/1 == s
        # Order
        @test s*val == val*s
        @test s/val == s*1/val
        @test val/s == val*1/s
    end
    # General modification
    let s = SACtr(rand(100), 2.0)
        @test all(modify(s, x->x^2).t .== s.t.^2)
        @test isapprox(modify(s, x->rem(x, 0.1)).t, rem.(s.t, 0.1))
    end

    # Array of traces and constant
    let a = sample(:array), vals = rand(SAC.SACFloat, length(a))
        local N = length(a)
        # Array with constant
        val = vals[1]
        @test typeof(a .+ val) == Array{SACtr,1}
        n = rand(1:N)
        @test (a .+ val)[n] == a[n] + val
        @test (a .- val)[n] == a[n] - val
        @test (a .* val)[n] == a[n] * val
        @test (a ./ val)[n] == a[n] / val
        # Array with array
        @test typeof(a .+ vals) == Array{SACtr,1}
        @test a .+ vals == vals .+ a
        @test typeof(a .- vals) == Array{SACtr,1}
        @test a .- vals == -vals .+ a
        @test a .* vals == vals .* a
        @test typeof(a ./ vals) == Array{SACtr,1}
        @test a ./ vals == a .* 1 ./ vals
        @test all(a .+ zeros(N) .== a)
        @test all(a .- zeros(N) .== a)
        @test all(a .* ones(N) .== a)
        @test all(a ./ ones(N) .== a)
    end

    # Multiple trace versions
    let s1 = SAC.sample(), s2 = SAC.sample() + 1
        @test_throws ArgumentError SAC.check_traces_combine(s1, cut(s1, :a, 1, :e, -1))
        s1′ = deepcopy(s1)
        s1′.delta = 2s1.delta
        @test_throws ArgumentError SAC.check_traces_combine(s1, s1′)
        @test SAC.check_traces_combine(s1, s1)
        @test SAC.check_traces_combine(s1, s2)
        @test typeof(s1 + s2) == SACtr
        @test all((s1 + s2).t .== s1.t .+ s2.t)
        @test all((s1 - s2).t .== s1.t .- s2.t)
        @test all((s1 * s2).t .== s1.t .* s2.t)
        @test all((s1 / s2).t .== s1.t ./ s2.t)
        # Order
        @test s1 + s2 == s2 + s1
        @test s1 - s2 == -s2 + s1
        @test s1 * s2 == s2 * s1
        @test s1/s2 == s1*1/s2
    end

    ## Trend and mean
    # Single trace
    let s = SACtr([i%2 for i in 1:1000], 0.01)
        @test s[:depmen] ≈ 0.5
        @test rmean(s)[:depmen] ≈ 0.0
        rmean!(s)
        @test s[:depmen] ≈ 0.0
    end
    let s = SACtr(1:100, 0.1)
        @test all(isapprox.(rtrend(s)[:t], 0.0, atol=1e-4))
        rtrend!(s)
        @test all(isapprox.(s[:t], 0.0, atol=1e-4))
    end
    # Arrays
    let a = [SAC.sample() for _ in 1:5], s = SAC.sample()
        @test all(rmean(a) .== rmean(s))
        rmean!(a)
        rmean!(s)
        @test all(a .== s)
    end
    let a = [SAC.sample() for _ in 1:5], s = SAC.sample()
        @test all(rtrend(a) .== rtrend(s))
        rtrend!(a)
        rtrend!(s)
        @test all(s .== s)
    end

    ## Normalisation
    let s = SACtr(rand(SAC.SACFloat, 100), 0.1), s′ = deepcopy(s)
        @test all(normalise(s).t .== s′.t./maximum(abs, s.t))
        @test normalise!(s) == normalise(s′)
    end

    ## FFT
    let s = SACtr(rand(SAC.SACFloat, 10000), 0.01)
        f, S = fft(s)
        @test s.t ≈ irfft(S, s[:npts])
        @test length(f) == length(S) == s[:npts]÷2 + 1
        @test minimum(f) == zero(typeof(f[1]))
        @test maximum(f) ≈ 1/s[:delta]/2
        # Array version
        ss = [s for _ in 1:3]
        f, S = fft(ss)
        @test length(f) == length(S) == length(ss)
        @test all((f[1],) .== f)
        @test all((S[1],) .== S)
    end

    ## Interpolation
    let s = SACtr([0,1,2,3], 1), s′ = deepcopy(s)
        @test interpolate(s, delta=0.5) == interpolate(s, npts=7) == interpolate(s, n=2)
        @test interpolate(s, n=2).t ≈ range(0, stop=3, length=7)
        interpolate!(s, n=4)
        @test interpolate(s′, n=4) == s
    end

    ## Time shift
    let s = SACtr([0,1,2], 1), s′ = deepcopy(s)
        @test tshift(s, 1, false).t == [0, 0, 1]
        @test tshift(s, -1, false).t == [1, 2, 0]
        @test tshift(s, 1).t == [2, 0, 1]
        tshift!(s, 3)
        @test s == tshift(s′, 3)
    end
end


text:
ComponentArray{A}(::UndefInitializer, ax::Axes) where {A<:StaticArrays.StaticArray,Axes<:Tuple} =
    ComponentArray(similar(A), ax...)

# Sometimes you just need to get things done.
function Base.Broadcast.BroadcastStyle(::Type{s66}) where s66<:Transpose{T,s19} where s19<:ComponentArray{s62,s63,s64,s65} where s65<:Tuple{Vararg{AbstractAxis,N}} where s64<:StaticArrays.StaticArray{S,s62,N} where S<:Tuple where s63<:N where s62<:T where {T, N}
    return CAStyle(StaticArrays.StaticArrayStyle{2}(), getaxes(s19), 2)
end
function Base.Broadcast.BroadcastStyle(::Type{s66}) where s66<:Adjoint{T,s19} where s19<:ComponentArray{s62,s63,s64,s65} where s65<:Tuple{Vararg{AbstractAxis,N}} where s64<:StaticArrays.StaticArray{S,s62,N} where S<:Tuple where s63<:N where s62<:T where {T, N}
    return CAStyle(StaticArrays.StaticArrayStyle{2}(), getaxes(s19), 2)
end

# This is pretty unideal. It gets the answer right, but it allocates in the process. Need to
# figure this out because what's the point of StaticArrays if they are going to allocate?
function Base.copy(bc::Base.Broadcast.Broadcasted{Style,<:Ax,<:F,s19}) where {s19<:Tuple, Axes, Ax, F, Style<:CAStyle{s18,<:s17,<:s12}} where {s12, s17, s18<:StaticArrays.StaticArrayStyle}
    args = map(getdata, bc.args)
    # # style = BC.combine_styles(BC.BroadcastStyle.(typeof.(args))...)
    return ComponentArray{s17}(Base.copy(BC.broadcasted(bc.f, args...)))
end


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

using Test, Libdl

pkg_dir = dirname(Base.find_package("UMAT"))
if Sys.iswindows()
    lib_dir = joinpath(pkg_dir,"..","deps","usr","bin")
else
    lib_dir = joinpath(pkg_dir,"..","deps","usr","lib")
end

@test isfile(joinpath(lib_dir,"libelastic." * dlext))
@test isfile(joinpath(lib_dir,"libisotropic_plast_exp." * dlext))
@test isfile(joinpath(lib_dir,"libisotropic_plast_imp." * dlext))
@test isfile(joinpath(lib_dir,"libmises_umat." * dlext))


text:
# This file loads interface from OMEinsum.
# This file implements alternative driver for EinCode object with Tensordot.contract.
# OMEinsum outperforms Tensordot.contract under Zygote.jl but latter
#  is much more compatible (surviving Zygote.hessian).

# Disable OMEinsum backend by default.
use_einsum = false

"""
    einmm"ij,jk -> ik"(A,B)
String macro interface which understands `numpy.einsum`'s notation.
Translates strings into `EinCode`-structs that can be called to evaluate an `einsum`.

This specific block is copy-paste from under-Peter/OMEinsum.
"""
macro einmm_str(s::AbstractString)
    s = replace(s, " " => "")
    m = match(r"([\(\)a-z,α-ω]*)->([a-zα-ω]*)", s)
    m == nothing && throw(ArgumentError("invalid einsum specification $s"))
    sixs, siy = m.captures
    if '(' in sixs
        error("Nested Einsum not supported at the moment.")
    else
        iy  = Tuple(siy)
        ixs = Tuple(Tuple(ix) for ix in split(sixs,','))
        # TODO: Alias EinCode to switch to Tensordot.contract backend.
        return EinCode(ixs, iy)
    end
end

"""
    execute(s::EinCode{::Any, ::Any}, TL, TR)
Execute operation defined in an EinCode with OMEinsum dispatcher
or Tensordot.contract compatibility backend.
"""
function execute(s::EinCode, TL, TR)
    if use_einsum
        return s(TL, TR)
    else
        idL, idR, operm = ein_char2idx(s)
        return permutedims(contract(TL, TR, idL, idR), operm)
    end
end

"Converts Einstein-summation indices to dimension index of tensors."
function ein_char2idx(sxL::Tuple, sxR::Tuple, sy::Tuple)
    cidxL = []
    cidxR = []
    # Parse inputs.
    for (ixL, cxL)=enumerate(sxL)
        found = false
        for (ixR, cxR)=enumerate(sxR)
            if (cxL == cxR)
                cidxL = [cidxL..., ixL]
                cidxR = [cidxR..., ixR]
            end
        end
    end
    # Get tensordot out indices.
    lcL, lcR = ((lcT, cidxT) -> begin
        for ixT=reverse(sort(cidxT))
            deleteat!(lcT, ixT)
        end
        lcT
    end).([[sxL...], [sxR...]], [cidxL, cidxR])
    ostr = [lcL..., lcR...]
    # Get transpose operation for final output.
    operm = (cy -> begin
        for (iout, cout)=enumerate(ostr)
            if (cy == cout)
                return iout
            end
        end
    end).(sy)
    cidxL, cidxR, operm
end

"Converts Einstein-summation indices to dimension index of tensors."
function ein_char2idx(::EinCode{ixs, iy}) where {ixs, iy}
    # Currently only able to handle 2-tensor summations.
    sxL, sxR = ixs
    cidxL, cidxR, oidx = ein_char2idx(sxL, sxR, iy)
end

@adjoint ein_char2idx(sxL::Tuple, sxR::Tuple, sy::Tuple) = begin
    ein_char2idx(sxL, sxR, sy), _ -> nothing
end

@adjoint ein_char2idx(s::EinCode) = begin
    ein_char2idx(s), _ -> nothing
end


text:
module MaterialModels

using Reexport
@reexport using Tensors
using NLsolve
using Rotations
import ForwardDiff
import DiffResults
using StaticArrays

"""
    AbstractMaterial

Store material parameters here. It can also be used to store constant precomputed properties, e.g. the elastic stiffness tensor.
Ideally, the name should be chosen according to the first author of the initial publication of a model. For every `Material` there should be a keyword constructor
and a list of arguments in the docstrings. If possible, also include the reference to a publication.
"""
abstract type AbstractMaterial end

"""
    AbstractMaterialState

Store state variables here. For now, this should **not** be mutable, a new object should be constructed for every new state. (We can discuss if this is a good or a bad idea.)
"""
abstract type AbstractMaterialState end
abstract type AbstractResiduals end

"""
    StrainMeasure
Defines the type of strain measure the a material uses, i.e Deformation gradient, Green-Lagrange strain etc.
"""
abstract type StrainMeasure end

"""
    material_response(m::AbstractMaterial, Δε::SymmetricTensor{2,3}, state::AbstractMaterialState, Δt; cache, options)

Compute the stress, stress tangent and state variables for the given strain increment `Δε` and previous state `state`.

Instead of the strain increment, the total strain could be handed over. Good ideas on how to handle this in general are welcome.
For non-continuum kind of material models, the interface should be similar with stress-like and strain-like quantities.
(E.g. for cohesive laws traction instead of stress and displacement jump instead of strain.)
This function signature must be the same for all material models, even if they don't require all arguments.

"""
function material_response end

"""
    initial_material_state(::AbstractMaterial)

Return the `MaterialState` that belongs to the given `Material` and is initialized with zeros.
"""
function initial_material_state end

"""
    get_cache(m::AbstractMaterial)

Construct cache object for iteratively solving non-linear material models.

For material models which require an iterative solution procedure, it is recommended to allocate storage for the iterative solver only once and reuse it for all material points.
When multithreading is used, each threads needs its own cache.
"""
function get_cache end
"""
    update_cache!(cache::OnceDifferentiable, f)

Update the cache object with the residual function for the current time/load step.

As the residual functions depend i.a. on the strain increment, the function and its jacobian need to be updated for every load step.
"""
function update_cache! end

include("traits.jl")
include("LinearElastic.jl")
include("Plastic.jl")
include("CrystalViscoPlastic/slipsystems.jl")
include("CrystalViscoPlastic/CrystalViscoPlastic.jl")
include("CrystalViscoPlastic/CrystalViscoPlasticRed.jl")

include("FiniteStrain/neohook.jl")
include("FiniteStrain/yeoh.jl")
include("FiniteStrain/stvenant.jl")

include("nonlinear_solver.jl")
include("wrappers.jl")

export initial_material_state, get_cache
export material_response

export AbstractMaterial, AbstractMaterialState
export LinearElastic, Plastic, CrystalViscoPlastic
export LinearElasticState, PlasticState, CrystalViscoPlasticState
export AbstractDim, OneD, UniaxialStrain, UniaxialStress, PlaneStrain, PlaneStress
export ∂S∂C, ∂S∂E, ∂Pᵀ∂F

export NeoHook, Yeoh, StVenant
end


text:
const γ_exact = 7//5
const μ_exact = 1//100
ρ_g(t, x, y, z, ::Val{2}) = sin(pi*x)*cos(pi*t)*cos(pi*y) + 3
U_g(t, x, y, z, ::Val{2}) = (sin(pi*x)*cos(pi*t)*cos(pi*y) + 3)*sin(pi*x)*cos(pi*t)*cos(pi*y)
V_g(t, x, y, z, ::Val{2}) = (sin(pi*x)*cos(pi*t)*cos(pi*y) + 3)*sin(pi*x)*cos(pi*t)*cos(pi*y)
W_g(t, x, y, z, ::Val{2}) = 0
E_g(t, x, y, z, ::Val{2}) = sin(pi*x)*cos(pi*t)*cos(pi*y) + 100
Sρ_g(t, x, y, z, ::Val{2}) = pi*(2*sin(2*pi*x) - 2*sin(2*pi*y) - sin(pi*(2*t - 2*x)) + sin(pi*(2*t + 2*x)) + sin(pi*(2*t - 2*y)) - sin(pi*(2*t + 2*y)) + 2*sin(pi*(2*x + 2*y)) + sin(pi*(-2*t + 2*x + 2*y)) + sin(pi*(2*t + 2*x + 2*y)) + 10*cos(pi*(-t + x + y)) - 2*cos(pi*(t - x + y)) + 2*cos(pi*(t + x - y)) + 14*cos(pi*(t + x + y)))/8
SU_g(t, x, y, z, ::Val{2}) = pi*(-600*sin(pi*t)*sin(pi*x)^2*cos(pi*t)*cos(pi*y)^2 - 900*sin(pi*t)*sin(pi*x)*cos(pi*y) - 900*sin(pi*x)^3*sin(pi*y)*cos(pi*t)^3*cos(pi*y)^2 - 1800*sin(pi*x)^2*sin(pi*y)*cos(pi*t)^2*cos(pi*y) + 540*sin(pi*x)^2*cos(pi*t)^3*cos(pi*x)*cos(pi*y)^3 + 1080*sin(pi*x)*cos(pi*t)^2*cos(pi*x)*cos(pi*y)^2 + 7*pi*sin(pi*x)*cos(pi*t)*cos(pi*y) + pi*sin(pi*y)*cos(pi*t)*cos(pi*x) + 120*cos(pi*t)*cos(pi*x)*cos(pi*y))/300
SV_g(t, x, y, z, ::Val{2}) = pi*(-600*sin(pi*t)*sin(pi*x)^2*cos(pi*t)*cos(pi*y)^2 - 900*sin(pi*t)*sin(pi*x)*cos(pi*y) - 540*sin(pi*x)^3*sin(pi*y)*cos(pi*t)^3*cos(pi*y)^2 - 1080*sin(pi*x)^2*sin(pi*y)*cos(pi*t)^2*cos(pi*y) + 900*sin(pi*x)^2*cos(pi*t)^3*cos(pi*x)*cos(pi*y)^3 - 120*sin(pi*x)*sin(pi*y)*cos(pi*t) + 1800*sin(pi*x)*cos(pi*t)^2*cos(pi*x)*cos(pi*y)^2 + 7*pi*sin(pi*x)*cos(pi*t)*cos(pi*y) + pi*sin(pi*y)*cos(pi*t)*cos(pi*x))/300
SW_g(t, x, y, z, ::Val{2}) = 0
SE_g(t, x, y, z, ::Val{2}) = pi*(-300*sin(pi*t)*sin(pi*x)*cos(pi*y) + 480*sin(pi*x)^4*sin(pi*y)*cos(pi*t)^4*cos(pi*y)^3 + 1080*sin(pi*x)^3*sin(pi*y)*cos(pi*t)^3*cos(pi*y)^2 - 480*sin(pi*x)^3*cos(pi*t)^4*cos(pi*x)*cos(pi*y)^4 - 7*pi*sin(pi*x)^2*sin(pi*y)^2*cos(pi*t)^2 - 840*sin(pi*x)^2*sin(pi*y)*cos(pi*t)^2*cos(pi*y) - 1080*sin(pi*x)^2*cos(pi*t)^3*cos(pi*x)*cos(pi*y)^3 + 14*pi*sin(pi*x)^2*cos(pi*t)^2*cos(pi*y)^2 + 4*pi*sin(pi*x)*sin(pi*y)*cos(pi*t)^2*cos(pi*x)*cos(pi*y) - 42000*sin(pi*x)*sin(pi*y)*cos(pi*t) + 840*sin(pi*x)*cos(pi*t)^2*cos(pi*x)*cos(pi*y)^2 - 7*pi*cos(pi*t)^2*cos(pi*x)^2*cos(pi*y)^2 + 42000*cos(pi*t)*cos(pi*x)*cos(pi*y))/300
ρ_g(t, x, y, z, ::Val{3}) = sin(pi*x)*cos(pi*t)*cos(pi*y)*cos(pi*z) + 3
U_g(t, x, y, z, ::Val{3}) = (sin(pi*x)*cos(pi*t)*cos(pi*y)*cos(pi*z) + 3)*sin(pi*x)*cos(pi*t)*cos(pi*y)*cos(pi*z)
V_g(t, x, y, z, ::Val{3}) = (sin(pi*x)*cos(pi*t)*cos(pi*y)*cos(pi*z) + 3)*sin(pi*x)*cos(pi*t)*cos(pi*y)*cos(pi*z)
W_g(t, x, y, z, ::Val{3}) = (sin(pi*x)*cos(pi*t)*cos(pi*y)*cos(pi*z) + 3)*sin(pi*x)*sin(pi*z)*cos(pi*t)*cos(pi*y)
E_g(t, x, y, z, ::Val{3}) = sin(pi*x)*cos(pi*t)*cos(pi*y)*cos(pi*z) + 100
Sρ_g(t, x, y, z, ::Val{3}) = pi*(-sin(pi*t)*sin(pi*x)*cos(pi*y)*cos(pi*z) - 2*sin(pi*x)^2*sin(pi*y)*cos(pi*t)^2*cos(pi*y)*cos(pi*z)^2 - sin(pi*x)^2*sin(pi*z)^2*cos(pi*t)^2*cos(pi*y)^2 + sin(pi*x)^2*cos(pi*t)^2*cos(pi*y)^2*cos(pi*z)^2 - 3*sin(pi*x)*sin(pi*y)*cos(pi*t)*cos(pi*z) + 2*sin(pi*x)*cos(pi*t)^2*cos(pi*x)*cos(pi*y)^2*cos(pi*z)^2 + 3*sin(pi*x)*cos(pi*t)*cos(pi*y)*cos(pi*z) + 3*cos(pi*t)*cos(pi*x)*cos(pi*y)*cos(pi*z))
SU_g(t, x, y, z, ::Val{3}) = pi*(-600*sin(pi*t)*sin(pi*x)^2*cos(pi*t)*cos(pi*y)^2*cos(pi*z)^2 - 900*sin(pi*t)*sin(pi*x)*cos(pi*y)*cos(pi*z) - 900*sin(pi*x)^3*sin(pi*y)*cos(pi*t)^3*cos(pi*y)^2*cos(pi*z)^3 - 600*sin(pi*x)^3*sin(pi*z)^2*cos(pi*t)^3*cos(pi*y)^3*cos(pi*z) + 300*sin(pi*x)^3*cos(pi*t)^3*cos(pi*y)^3*cos(pi*z)^3 - 1800*sin(pi*x)^2*sin(pi*y)*cos(pi*t)^2*cos(pi*y)*cos(pi*z)^2 - 180*sin(pi*x)^2*sin(pi*z)^2*cos(pi*t)^3*cos(pi*x)*cos(pi*y)^3*cos(pi*z) - 900*sin(pi*x)^2*sin(pi*z)^2*cos(pi*t)^2*cos(pi*y)^2 + 540*sin(pi*x)^2*cos(pi*t)^3*cos(pi*x)*cos(pi*y)^3*cos(pi*z)^3 + 900*sin(pi*x)^2*cos(pi*t)^2*cos(pi*y)^2*cos(pi*z)^2 - 360*sin(pi*x)*sin(pi*z)^2*cos(pi*t)^2*cos(pi*x)*cos(pi*y)^2 + 1080*sin(pi*x)*cos(pi*t)^2*cos(pi*x)*cos(pi*y)^2*cos(pi*z)^2 + 10*pi*sin(pi*x)*cos(pi*t)*cos(pi*y)*cos(pi*z) + pi*sin(pi*y)*cos(pi*t)*cos(pi*x)*cos(pi*z) - pi*cos(pi*t)*cos(pi*x)*cos(pi*y)*cos(pi*z) + 120*cos(pi*t)*cos(pi*x)*cos(pi*y)*cos(pi*z))/300
SV_g(t, x, y, z, ::Val{3}) = pi*(-600*sin(pi*t)*sin(pi*x)^2*cos(pi*t)*cos(pi*y)^2*cos(pi*z)^2 - 900*sin(pi*t)*sin(pi*x)*cos(pi*y)*cos(pi*z) + 180*sin(pi*x)^3*sin(pi*y)*sin(pi*z)^2*cos(pi*t)^3*cos(pi*y)^2*cos(pi*z) - 540*sin(pi*x)^3*sin(pi*y)*cos(pi*t)^3*cos(pi*y)^2*cos(pi*z)^3 - 600*sin(pi*x)^3*sin(pi*z)^2*cos(pi*t)^3*cos(pi*y)^3*cos(pi*z) + 300*sin(pi*x)^3*cos(pi*t)^3*cos(pi*y)^3*cos(pi*z)^3 + 360*sin(pi*x)^2*sin(pi*y)*sin(pi*z)^2*cos(pi*t)^2*cos(pi*y) - 1080*sin(pi*x)^2*sin(pi*y)*cos(pi*t)^2*cos(pi*y)*cos(pi*z)^2 - 900*sin(pi*x)^2*sin(pi*z)^2*cos(pi*t)^2*cos(pi*y)^2 + 900*sin(pi*x)^2*cos(pi*t)^3*cos(pi*x)*cos(pi*y)^3*cos(pi*z)^3 + 900*sin(pi*x)^2*cos(pi*t)^2*cos(pi*y)^2*cos(pi*z)^2 - 120*sin(pi*x)*sin(pi*y)*cos(pi*t)*cos(pi*z) + pi*sin(pi*x)*sin(pi*y)*cos(pi*t)*cos(pi*z) + 1800*sin(pi*x)*cos(pi*t)^2*cos(pi*x)*cos(pi*y)^2*cos(pi*z)^2 + 10*pi*sin(pi*x)*cos(pi*t)*cos(pi*y)*cos(pi*z) + pi*sin(pi*y)*cos(pi*t)*cos(pi*x)*cos(pi*z))/300
SW_g(t, x, y, z, ::Val{3}) = pi*(-600*sin(pi*t)*sin(pi*x)^2*cos(pi*t)*cos(pi*y)^2*cos(pi*z) - 900*sin(pi*t)*sin(pi*x)*cos(pi*y) - 900*sin(pi*x)^3*sin(pi*y)*cos(pi*t)^3*cos(pi*y)^2*cos(pi*z)^2 - 240*sin(pi*x)^3*sin(pi*z)^2*cos(pi*t)^3*cos(pi*y)^3 + 840*sin(pi*x)^3*cos(pi*t)^3*cos(pi*y)^3*cos(pi*z)^2 - 1800*sin(pi*x)^2*sin(pi*y)*cos(pi*t)^2*cos(pi*y)*cos(pi*z) + 900*sin(pi*x)^2*cos(pi*t)^3*cos(pi*x)*cos(pi*y)^3*cos(pi*z)^2 + 2160*sin(pi*x)^2*cos(pi*t)^2*cos(pi*y)^2*cos(pi*z) - pi*sin(pi*x)*sin(pi*y)*cos(pi*t) + 1800*sin(pi*x)*cos(pi*t)^2*cos(pi*x)*cos(pi*y)^2*cos(pi*z) - 120*sin(pi*x)*cos(pi*t)*cos(pi*y) + 10*pi*sin(pi*x)*cos(pi*t)*cos(pi*y) + pi*cos(pi*t)*cos(pi*x)*cos(pi*y))*sin(pi*z)/300
SE_g(t, x, y, z, ::Val{3}) = pi*(45*(1 - cos(2*pi*x))^2*(1 - cos(4*pi*z))*(cos(2*pi*t) + 1)^2*(cos(2*pi*y) + 1)^2/128 - 300*sin(pi*t)*sin(pi*x)*cos(pi*y)*cos(pi*z) + 240*sin(pi*x)^4*sin(pi*y)*sin(pi*z)^2*cos(pi*t)^4*cos(pi*y)^3*cos(pi*z)^2 + 480*sin(pi*x)^4*sin(pi*y)*cos(pi*t)^4*cos(pi*y)^3*cos(pi*z)^4 + 60*sin(pi*x)^4*sin(pi*z)^4*cos(pi*t)^4*cos(pi*y)^4 - 120*sin(pi*x)^4*cos(pi*t)^4*cos(pi*y)^4*cos(pi*z)^4 + 540*sin(pi*x)^3*sin(pi*y)*sin(pi*z)^2*cos(pi*t)^3*cos(pi*y)^2*cos(pi*z) + 1080*sin(pi*x)^3*sin(pi*y)*cos(pi*t)^3*cos(pi*y)^2*cos(pi*z)^3 - 240*sin(pi*x)^3*sin(pi*z)^2*cos(pi*t)^4*cos(pi*x)*cos(pi*y)^4*cos(pi*z)^2 + 180*sin(pi*x)^3*sin(pi*z)^2*cos(pi*t)^3*cos(pi*y)^3*cos(pi*z) - 480*sin(pi*x)^3*cos(pi*t)^4*cos(pi*x)*cos(pi*y)^4*cos(pi*z)^4 - 360*sin(pi*x)^3*cos(pi*t)^3*cos(pi*y)^3*cos(pi*z)^3 - 3*pi*sin(pi*x)^2*sin(pi*y)^2*sin(pi*z)^2*cos(pi*t)^2 - 7*pi*sin(pi*x)^2*sin(pi*y)^2*cos(pi*t)^2*cos(pi*z)^2 - 7*pi*sin(pi*x)^2*sin(pi*y)*sin(pi*z)^2*cos(pi*t)^2*cos(pi*y) - 840*sin(pi*x)^2*sin(pi*y)*cos(pi*t)^2*cos(pi*y)*cos(pi*z)^2 - 3*pi*sin(pi*x)^2*sin(pi*y)*cos(pi*t)^2*cos(pi*y)*cos(pi*z)^2 - 540*sin(pi*x)^2*sin(pi*z)^2*cos(pi*t)^3*cos(pi*x)*cos(pi*y)^3*cos(pi*z) - 420*sin(pi*x)^2*sin(pi*z)^2*cos(pi*t)^2*cos(pi*y)^2 + 4*pi*sin(pi*x)^2*sin(pi*z)^2*cos(pi*t)^2*cos(pi*y)^2 - 1080*sin(pi*x)^2*cos(pi*t)^3*cos(pi*x)*cos(pi*y)^3*cos(pi*z)^3 + 16*pi*sin(pi*x)^2*cos(pi*t)^2*cos(pi*y)^2*cos(pi*z)^2 + 420*sin(pi*x)^2*cos(pi*t)^2*cos(pi*y)^2*cos(pi*z)^2 + 4*pi*sin(pi*x)*sin(pi*y)*cos(pi*t)^2*cos(pi*x)*cos(pi*y)*cos(pi*z)^2 - 42000*sin(pi*x)*sin(pi*y)*cos(pi*t)*cos(pi*z) + 7*pi*sin(pi*x)*sin(pi*z)^2*cos(pi*t)^2*cos(pi*x)*cos(pi*y)^2 + 3*pi*sin(pi*x)*cos(pi*t)^2*cos(pi*x)*cos(pi*y)^2*cos(pi*z)^2 + 840*sin(pi*x)*cos(pi*t)^2*cos(pi*x)*cos(pi*y)^2*cos(pi*z)^2 + 42000*sin(pi*x)*cos(pi*t)*cos(pi*y)*cos(pi*z) - 3*pi*sin(pi*z)^2*cos(pi*t)^2*cos(pi*x)^2*cos(pi*y)^2 - 7*pi*cos(pi*t)^2*cos(pi*x)^2*cos(pi*y)^2*cos(pi*z)^2 + 42000*cos(pi*t)*cos(pi*x)*cos(pi*y)*cos(pi*z))/300

max_stars_repo_name: christophernhill-dev/CLIMA

max_issues_repo_name: christophernhill-dev/CLIMA

max_forks_repo_name: christophernhill-dev/CLIMA


text:
#### Plotting the interactive effects of temperature and fertilisation

# In order to visualize our results as heatmaps (or a contour plots), we first split df by Z and turn into 2D arrays:

# subset df by body mass ratio
dfZ1 = filter(:Z => Z -> Z == 1, df)
dfZ10 = filter(:Z => Z -> Z == 10, df)
dfZ100 = filter(:Z => Z -> Z == 100, df)

# empty 2D arrays
TK_array_persistence1 = fill(NaN, length(allT), length(allI_K)) #preallocate a 2d array for persistence values Z = 1
TK_array_persistence10 = fill(NaN, length(allT), length(allI_K))
TK_array_persistence100 = fill(NaN, length(allT), length(allI_K))
TK_array_cc = fill(NaN, length(allT), length(allI_K))

###### Z = 1
# fill 2D array for Z = 1
for (j,k) in enumerate(allI_K)
	for (i,t) in enumerate(allT)
		tmp = dfZ1[(dfZ1.I_K .== k) .& (dfZ1.T .== t),:]
		TK_array_persistence1[i,j] = tmp.sp_survival[1]
	end
end

# plot for Z = 1
palette2 = [RGB(190/255, 242/255, 170/255), RGB(149/255, 188/255, 219/255), RGB(30/255, 80/255, 150/255)]
p1 = heatmap(string.(allI_K), string.(allT .- 273.15), TK_array_persistence1
    , title  = "3b, Z = 1", classes = 3, color = palette2
    , xlabel = "Fertilisation", ylabel ="Temperature (C)"
    )

###### Z = 10
# fill 2D array for Z = 10
for (j,k) in enumerate(allI_K)
	for (i,t) in enumerate(allT)
		tmp = dfZ10[(dfZ10.I_K .== k) .& (dfZ10.T .== t),:]
		TK_array_persistence10[i,j] = tmp.sp_survival[1]
	end
end

# plot for Z = 10
p10 = heatmap(string.(allI_K), string.(allT .- 273.15), TK_array_persistence10
    , title  = "3c, Z = 10", classes = 3, color = palette2
    , xlabel = "Fertilisation", ylabel ="Temperature (C)" )

###### Z = 100
# fill 2D array for Z = 100
for (j,k) in enumerate(allI_K)
	for (i,t) in enumerate(allT)
		tmp = dfZ100[(dfZ100.I_K .== k) .& (dfZ100.T .== t),:]
		TK_array_persistence100[i,j] = tmp.sp_survival[1]
	end
end

# plot for Z = 1
p100 = heatmap(string.(allI_K), string.(allT .- 273.15), TK_array_persistence100
    , title  = "3d, Z = 100", classes = 3, color = palette2
    , xlabel = "Fertilisation", ylabel ="Temperature (C)" )

###### Log 10 of the carrying capacity of the system
# fill 2D array for Z = 1 of carrying capacity
for (j,k) in enumerate(allI_K)
	for (i,t) in enumerate(allT)
		tmp = dfZ1[(dfZ1.I_K .== k) .& (dfZ1.T .== t),:]
		TK_array_cc[i,j] = log10(tmp.K[1])
	end
end

# plot for Z = 1
palette1 = reverse([RGB(255/255, 255/255, 255/255)
    , RGB(245/255, 242/255, 81/255)
    , RGB(255/255, 242/255, 0/255)
    , RGB(255/255, 196/255, 0/255)
    , RGB(255/255, 145/255, 0/255)
    , RGB(255/255, 111/255, 0/255)
    , RGB(255/255, 25/255, 0/255)])
pcc = heatmap(string.(allI_K), string.(allT .- 273.15), TK_array_cc
    , title  = "3a, log10(K)", classes = 7, color = palette1
    , xlabel = "Fertilisation", ylabel ="Temperature (C)" )

l3 = @layout [a b; c d]
plot(pcc, p1, p10, p100, layout = l3)


text:
# # Example: paired-end sequencing
#
# This is an example generated from this source
# file: [`pe-example.jl`](@__REPO_ROOT_URL__/examples/sequencing/pe-example.jl)
# You are seeing the
#md # online documentation version. The corresponding notebook can be found
#md # here: [`pe-example.ipynb`](@__NBVIEWER_ROOT_URL__/man/sequencing/examples/pe-example.ipynb)
#nb # jupyter notebook version. The corresponding online documentation page can
#nb # be found here: [`pe-example.html`](https://bioinfologics.github.io/Pseudoseq.jl/latest/examples/sequencing/pe-example),
# and the script can be found here: [`pe-example.jl`](./pe-example.jl)

# For the simulation we are going to:
#
# 1. Create a pool of 5000 copies of a reference genome.
# 2. Fragment the DNA molecules in the pool, to an average length of 700bp.
# 3. Subsample the molecules in the pool to achieve approximatly 50x coverage.
# 4. Create a set of 250bp paired-end reads.
# 5. Apply errors to the paired-end reads at a rate of 0.001 (.1%).
# 6. Generate an output FASTQ file.

using Pseudoseq.Sequencing

# Starting with a FASTA formatted file containing the genome we want to sequence,
# we create a pool with a copy of the genome.

m = Molecules("ecoli-ref.fasta")

# We want to first have 5000 full copies of the genome, so we will make an amplifier.

amp = amplify(5000)

# Next we want a fragmenter to make a pool of shorter DNA molecules, with an average
# length of 700bp.

frag = fragment(700bp)

# We will want a filter to make sure no extreme sized molecules make their way through.

size_filter = select(x -> 900 >= length(x) >= 450)

# Next we create a subsampler which will randomly sample molecules, to give us
# a desired expected coverage. We will want our reads to be 250bp in length.

ssmpl = subsample(50X, 2 * 250bp)

# We then want a read-maker that will give us paired end reads. We want our
# reads to be 250bp in length.

readmaker = makereads(2 * 250bp)

# Once we have reads, we will mark positions in the reads that are incorrectly
# detected by the sequencer: errors.
# We will construct a `FixedProbSubstitutions` function with a per base error
# probability of 0.001 and pass it to the `make_substitutions` method.
# This will make errors fall totally randomly over each read.

errmaker = make_substitutions(FixedProbSubstitutions(0.001))

# Now we can push our molecules through a pipeline of these processors, and out
# to a FASTQ file:

m |> amp |> frag |> size_filter |> ssmpl |> readmaker |> errmaker |> generate("pe-reads.fastq")

# You can also compose the processors together into one whole function.
# Typing \circ in the julia repl and then hitting tab gives you the circular
# composition symbol. Note how pipelining above progresses from left to right,
# but composition is right to left in order.

my_protocol = errmaker ∘ readmaker ∘ ssmpl ∘ size_filter ∘ frag ∘ amp

m |> my_protocol |> generate("pe-reads.fastq")


text:
module CommunalHelperCassetteFallingBlock

using ..Ahorn, Maple
using Ahorn.CommunalHelper
using Ahorn.CommunalHelperEntityPresets: CustomCassetteBlockData

@mapdefdata Entity "CommunalHelper/CassetteFallingBlock" CassetteFallingBlock CustomCassetteBlockData

const placements = Ahorn.PlacementDict(
    "Cassette Falling Block ($index - $color) (Communal Helper)" => Ahorn.EntityPlacement(
        CassetteFallingBlock,
        "rectangle",
        Dict{String,Any}(
            "index" => index,
        ),
    ) for (color, index) in cassetteColorNames
)

Ahorn.editingOptions(entity::CassetteFallingBlock) = Dict{String,Any}(
    "index" => cassetteColorNames
)

Ahorn.nodeLimits(entity::CassetteFallingBlock) = 1, 1

Ahorn.minimumSize(entity::CassetteFallingBlock) = 16, 16
Ahorn.resizable(entity::CassetteFallingBlock) = true, true

Ahorn.selection(entity::CassetteFallingBlock) = Ahorn.getEntityRectangle(entity)

const block = "objects/cassetteblock/solid"

function Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::CassetteFallingBlock)
    width = Int(get(entity.data, "width", 32))
    height = Int(get(entity.data, "height", 32))

    index = Int(get(entity.data, "index", 0))

    hexColor = String(get(entity.data, "customColor", ""))
    if hexColor != "" && length(hexColor) == 6
        renderCassetteBlock(ctx, 0, 0, width, height, index, hexToRGBA(hexColor))
    else
        renderCassetteBlock(ctx, 0, 0, width, height, index)
    end
end

end


text:
"""
	las2pointcloud(fnames::String...) -> PointCloud

Read more than one file `.las` and extrapolate the LAR model and the color of each point.
"""
# function las2pointcloud(fnames::String...)::PointCloud
# 	Vtot = Array{Float64,2}(undef, 3, 0)
# 	rgbtot = Array{LasIO.N0f16,2}(undef, 3, 0)
# 	for fname in fnames
# 		V = las2larpoints(fname)
# 		rgb = las2color(fname)
# 		Vtot = hcat(Vtot,V)
# 		rgbtot = hcat(rgbtot,rgb)
# 	end
# 	return PointCloud(Vtot,rgbtot)
# end

function las2pointcloud(fnames::String...)::PointCloud
	Vtot = Array{Float64,2}(undef, 3, 0)
	rgbtot = Array{LasIO.N0f16,2}(undef, 3, 0)

	for fname in fnames
		header, laspoints = read_LAS_LAZ(fname)
		npoints = length(laspoints)

		# COORDS
		x = [LasIO.xcoord(laspoints[k], header) for k in 1:npoints]
		y = [LasIO.ycoord(laspoints[k], header) for k in 1:npoints]
		z = [LasIO.zcoord(laspoints[k], header) for k in 1:npoints]
		V = vcat(x',y',z')

		# RGB
		type = LasIO.pointformat(header)
		rgb = Array{LasIO.N0f16,2}(undef, 3, 0)
		if type != LasPoint0 && type != LasPoint1
			r = LasIO.ColorTypes.red.(laspoints)
			g = LasIO.ColorTypes.green.(laspoints)
			b = LasIO.ColorTypes.blue.(laspoints)
			rgb = vcat(r',g',b')
		end

		Vtot = hcat(Vtot,V)
		rgbtot = hcat(rgbtot,rgb)
	end

	return PointCloud(Vtot,rgbtot)
end

"""
	las2larpoints(file::String) -> Points

Return coordinates of points in LAS file.
"""
function las2larpoints(file::String)::Points
	header, laspoints = read_LAS_LAZ(file)
	npoints = length(laspoints)
	x = [LasIO.xcoord(laspoints[k], header) for k in 1:npoints]
	y = [LasIO.ycoord(laspoints[k], header) for k in 1:npoints]
	z = [LasIO.zcoord(laspoints[k], header) for k in 1:npoints]
	return vcat(x',y',z')
end

"""
	las2aabb(file::String) -> AABB

Return LAS file's bounding box.
"""
function las2aabb(file::String)::AABB
	header = nothing
	open(file,"r") do s
	  LasIO.skiplasf(s)
	  header = read(s, LasHeader)
	end
	#header = LasIO.read(fname, LasIO.LasHeader)
	aabb = LasIO.boundingbox(header)
	return AABB(aabb.xmax, aabb.xmin, aabb.ymax, aabb.ymin, aabb.zmax, aabb.zmin)
end

function las2aabb(header::LasHeader)::AABB
	aabb = LasIO.boundingbox(header)
	return AABB(aabb.xmax, aabb.xmin, aabb.ymax, aabb.ymin, aabb.zmax, aabb.zmin)
end

"""
	las2color(file::String)::Points

Return color, rgb, associated to each point in LAS file.
"""
function las2color(file::String)::Points
	header, laspoints =  read_LAS_LAZ(file)
	npoints = length(laspoints)
	type = LasIO.pointformat(header)
	if type != LasPoint0 && type != LasPoint1
		r = LasIO.ColorTypes.red.(laspoints)
		g = LasIO.ColorTypes.green.(laspoints)
		b = LasIO.ColorTypes.blue.(laspoints)
		return vcat(r',g',b')
	end
	return rgbtot = Array{LasIO.N0f16,2}(undef, 3, 0)
end

"""
	color(p::LasPoint, header::LasHeader)

Return color of one point in LAS file.
"""
function color(p::LasPoint, header::LasHeader)
	type = LasIO.pointformat(header)
	if type != LasPoint0 && type != LasPoint1
		r = LasIO.ColorTypes.red(p)
		g = LasIO.ColorTypes.green(p)
		b = LasIO.ColorTypes.blue(p)
		return vcat(r',g',b')
	end
	return rgbtot = Array{LasIO.N0f16,2}(undef, 3, 0)
end

"""
	 xyz(p::LasPoint, h::LasHeader)

Return coords of this point `p`.
"""
function xyz(p::LasPoint, h::LasHeader)
	return [LasIO.xcoord(p, h); LasIO.ycoord(p, h); LasIO.zcoord(p, h)]
end

"""
	read_LAS_LAZ(fname::String)

Read point cloud files: LAS or LAZ.
"""
function read_LAS_LAZ(fname::String)
	if endswith(fname,".las")
		header, laspoints = LasIO.FileIO.load(fname)
	elseif endswith(fname,".laz")
		header, laspoints = LazIO.load(fname)
	end
	return header,laspoints
end