text:
# class constructors
# convenience type constructor:
# enable `R[:x]` as short for `UnivariatePolynomial{R,:x}`
getindex(R::Type{<:Ring}, X::Symbol) = UnivariatePolynomial{R,X}

### Constructors
category_trait(P::Type{<:Polynomial}) = category_trait_poly(P, category_trait(basetype(P)))
category_trait_poly(::Type{<:UnivariatePolynomial}, Z::Type{<:FieldTrait}) = EuclidianDomainTrait
category_trait_poly(::Type{<:Polynomial}, Z::Type{<:UniqueFactorizationDomainTrait}) = UniqueFactorizationDomainTrait
category_trait_poly(::Type, Z::Type{<:IntegralDomainTrait}) = IntegralDomainTrait
category_trait_poly(::Type, ::Type) = CommutativeRingTrait
category_trait_fraction(::Type{<:IntegralDomainTrait}) = FieldTrait
category_trait_fraction(::Type) = CommutativeRingTrait
isprimemod(Z::Type{<:Quotient{<:UnivariatePolynomial}}) = isirreducible(modulus(Z))

basetype(::Type{<:Polynomial{T}}) where T = T
depth(::Type{T}) where T<:Polynomial = depth(basetype(T)) + 1
function lcunit(a::Polynomial)
    lco = LC(a)
    isunit(lco) ? lco : one(lco)
end

### access to polynomial coefficients
function getindex(u::UnivariatePolynomial, i::Integer)
    0 <= i <= deg(u) ? u.coeff[i+1] : zero(basetype(u))
end

"""
    varname(P)

Return variable name of univariate polynomial or polynomial type `P` as a symbol.
"""
varname(::Type{T}) where {R,X,T<:UnivariatePolynomial{R,X}} = X
varname(::T) where T<:Polynomial = varname(T)

"""
    iscoprime(a, b)

Return if there is a common non-unit divisor of `a` and `b`.
"""
function iscoprime(a::T, b::T) where T<:UnivariatePolynomial
    g, u, v, f = pgcdx(a, b)
    isunit(g) || !isunit(f)
end
"""
    isinvertible(a, b)

Return if there is an inverse of `a` mod ` b`.
"""
function isinvertible(a::T, b::T) where T<:UnivariatePolynomial
    g, u, v, f = pgcdx(a, b)
    isunit(g) && isunit(f)
end

function issimpler(a::T, b::T) where T<:Polynomial
    da, db = deg(a), deg(b)
    da < db || da == db && issimpler(LC(a), LC(b))
end

# promotion and conversion
_promote_rule(::Type{UnivariatePolynomial{R,X}}, ::Type{UnivariatePolynomial{S,Y}}) where {X,Y,R,S} = Base.Bottom # throw(DomainError((X,Y), "cannot promote univariate polynomials with differnet variables"))
_promote_rule(::Type{UnivariatePolynomial{R,X}}, ::Type{UnivariatePolynomial{S,X}}) where {X,R,S} = UnivariatePolynomial{promote_type(R,S),X}
_promote_rule(::Type{UnivariatePolynomial{R,X}}, ::Type{S}) where {X,R,S<:Ring} = UnivariatePolynomial{promote_type(R,S),X}
promote_rule(::Type{UnivariatePolynomial{R,X}}, ::Type{S}) where {X,R,S<:Integer} = UnivariatePolynomial{promote_type(R,S),X}
promote_rule(::Type{UnivariatePolynomial{R,X}}, ::Type{S}) where {X,R,S<:Rational} = UnivariatePolynomial{promote_type(R,S),X}

(P::Type{<:UnivariatePolynomial{S}})(a::S) where {S} = P([a])
(P::Type{<:UnivariatePolynomial{S}})(a::T) where {S,T} = P([S(a)])

# convert coefficient vector to polynomial
function UnivariatePolynomial{T,X}(v::Vector{S}) where {X,T<:Ring,S<:T}
    x = Symbol(X)
    n = length(v)
    while n > 0 && iszero(v[n])
        n -= 1
    end
    if n < length(v)
        v = copyto!(similar(v, n), 1, v, 1, n)
    end
    UnivariatePolynomial{S,x}(v, NOCHECK)
end
# convert coefficient vector to polynomial, element type of vector determines class
UnivariatePolynomial(x::Symbol, v::Vector{S}) where {S} = UnivariatePolynomial{S,x}(v)
UnivariatePolynomial(x::Symbol, v::S) where {S} = UnivariatePolynomial{S,x}([v])
# convert polynomial with same symbol and type aliasing original
UnivariatePolynomial{S,X}(p::UnivariatePolynomial{S,X}) where {X,S<:Ring} = p
# convert polynomial with different coefficient type
function UnivariatePolynomial{S,X}(p::UnivariatePolynomial{T,Y}) where {X,Y,S,T}
    if X == Y
        co = [S(c) for c in p.coeff]
        UnivariatePolynomial{S,X}(co)
    elseif S <: UnivariatePolynomial
        co = [S(p)]
        UnivariatePolynomial{S,X}(co)
    else
        throw(DomainError((X,Y), "cannot convert polynomials with differing variable names"))
    end
end

function monom(P::Type{<:UnivariatePolynomial}, xv::Vector{<:Integer})
    length(xv) > 1 && throw(ArgumentError("univariate monom has only one variable"))
    monom(P, xv...)
end

function monom(P::Type{<:UnivariatePolynomial{S}}, k::Integer=1) where S
    k = max(k, -1)
    v = zeros(S,k+1)
    v[k+1] = one(S)
    P(v)
end

"""
    UnivariatePolynomials{X,S}(vector)

Construct a new polynomial Ring-element.
Allow all coefficient classes, which can be mapped to S, that means
the canonical homomorphism is used.
"""
function UnivariatePolynomial{S,X}(v::AbstractVector) where {X,S}
    isempty(v) ? UnivariatePolynomial{S,X}(S[]) : UnivariatePolynomial{S,X}([S(x) for x = v])
end

# canonical embedding homomorphism from base ring
UnivariatePolynomial(r::R) where {R<:Ring} = UnivariatePolynomial{R,:x}([r])

# make new copy
copy(p::UnivariatePolynomial) = typeof(p)(copy(p.coeff))

### Arithmetic
function +(p::T, q::T) where T<:UnivariatePolynomial
    vp = p.coeff
    vq = q.coeff
    np = length(vp)
    nq = length(vq)
    np >= nq || return +(q, p)
    v = copyto!(similar(vp, np), 1, vp, 1, np)
    for i = 1:nq
        v[i] += vq[i]
    end
    while np > 0 && iszero(v[np])
        np -= 1
    end
    if np < length(v)
        resize!(v, np)
    end
    T(v, NOCHECK)
end
+(p::T, q::Integer) where {S,T<:UnivariatePolynomial{S}} = p + T([S(q)])
+(q::Integer, p::T) where {S,T<:UnivariatePolynomial{S}} = +(p, q)
+(q::S, p::T) where {S<:Ring,T<:UnivariatePolynomial{S}} = +(p, q)

function -(p::T) where T<:UnivariatePolynomial
    vp = p.coeff
    np = length(vp)
    v = similar(vp)
    for i = 1:np
        v[i] = -vp[i]
    end
    T(v, NOCHECK)
end
-(p::T, q::T) where T<:UnivariatePolynomial = +(p, -q)
-(p::T, q::Ring) where T<:UnivariatePolynomial = +(p, -q)
-(p::T, q::Integer) where T<:UnivariatePolynomial = +(p, -q)
-(q::Integer, p::T) where T<:UnivariatePolynomial = +(-p, q)
-(q::S, p::T) where {S<:Ring,T<:UnivariatePolynomial{S}} = +(-p, q)

function *(p::T, q::T) where T<:UnivariatePolynomial
    vp = p.coeff
    vq = q.coeff
    np = length(p.coeff)
    nq = length(q.coeff)
    np <= nq || return *(q, p)
    nv = np + nq - 1
    if np == 0 || nq == 0
        return zero(T)
    elseif ismonom(p)
        v = multmono(p, np, vp, q, nq, vq)
    elseif ismonom(q)
        v = multmono(q, nq, vq, p, np, vp)
    else
        v = similar(vp, nv)
        for k = 1:nv
            i1 = max(k+1-np,1)
            i2 = min(k,nq)
            vk = vp[k-i1+1] * vq[i1]
            for j = i1+1:i2
                vk += vp[k-j+1] * vq[j]
            end
            v[k] = vk
        end
    end
    v === vp ? p : v === vq ? q : T(v, NOCHECK)
end

function *(p::UnivariatePolynomial{R}, q::R) where R<:Ring
    if iszero(q)
        zero(p)
    else
        T = typeof(p)
        S = basetype(T)
        # make broadcast recognize q as scalar
        T(p.coeff .* Ref(q), NOCHECK)
    end
end
*(p::UnivariatePolynomial{S}, q::Integer) where {S} = *(p, S(q))
*(q::Integer, p::UnivariatePolynomial) = *(p, q)
*(q::R, p::UnivariatePolynomial{R}) where R<:Ring = *(p, q)

function /(p::T, q::T) where {S,T<:UnivariatePolynomial{S}}
    d, r = divrem(p, q)
    iszero(r) ? d : throw(DomainError((p, q), "cannot divide a / b"))
end
function /(p::T, q::Ring) where {S,T<:UnivariatePolynomial{S}}
    T(p.coeff ./ S(q), NOCHECK)
end
/(p::UnivariatePolynomial{S}, q::Integer) where S = /(p, S(q))

function ^(p::P, k::Integer) where P<:Polynomial
    n = deg(p)
    if n == 0
        P(LC(p)^k)
    elseif k < 0
        throw(DomainError((p,k), "polynom power negative exponent"))
    elseif k == 0
        one(p)
    elseif k == 1 || n < 0
        p
    elseif n > 0 && ismonom(p)
        indv = multideg(p)
        lco = LC(p)
        mon = monom(P, k * indv)
        if !isone(lco) && !isempty(indv)
            mon.coeff[end] = lco^k
        end
        mon
    else
        Base.power_by_squaring(p, k)
    end
end

"""
    multideg(p::Polynomial)

Return vector of variable exponents of the leading monomial of `p`.
"""
multideg(p::UnivariatePolynomial) = [deg(p)]

function smul!(v::Vector{S}, r, m::S) where S
    for i in r
        v[i] *= m
    end
    v
end

# division and remainder algorithm. d, r = divrem(p, q) => d * q + r == p.
# if the third argument is true, perform `pseudo-division` by multiplying
# the divident by Ring element in order to allow division.
# In the other case, divide by leading term of q. If remainder is not zero
# the degree of d is not reduced (lower than that of q).
# Attempt to divide by null polynomial throws DomainError.
function divrem(vp::Vector{S}, vq::Vector{S}, ::Val{F}) where {S<:Ring,F}
    np = length(vp)
    nq = length(vq)
    nq > 0 || throw(DomainError(vq, "Cannot divide by zero polynomial."))
    f = one(S)
    if np < nq
        return S[], vp, f
    end
    divi = vq[nq]
    fac = F && !isunit(divi)
    n = fac ? nq - 1 : np
    if nq == 1
        isone(divi) && return vp, S[], f
        vf = copy(vp)
        if fac
            f = divi / gcd(divi, gcd(vp))
            !isone(f) && smul!(vf, 1:np, f)
        end
        vr = similar(vf)
        for i = 1:np
            vf[i], vr[i] = divrem(vf[i], divi)
        end
    else
        vf = similar(vp, np - nq + 1)
        vr = copy(vp)
        for i = np:-1:nq
            vri = vr[i]
            multi, r = divrem(vri, divi)
            if fac && !iszero(r)
                g = divi / gcd(divi, r)
                f *= g
                vri *= g
                smul!(vr, 1:i-1, g)
                smul!(vf, i-nq+2:np-nq+1, g)
                multi = vri / divi
            end
            vr[i] = r
            for j = 1:nq-1
                vr[j+i-nq] -= vq[j] * multi
            end
            vf[i-nq+1] = multi
        end
    end
    while n > 0 && iszero(vr[n])
        n -= 1
    end
    resize!(vr, n)
    n = length(vf)
    while n > 0 && iszero(vf[n])
        n -= 1
    end
    resize!(vf, n)
    vf, vr, f
end

function rem(p::T, q::T) where T<:UnivariatePolynomial
    m = deg(p)
    n = deg(q)
    m < n && return p
    uc = LC(q)
    isunit(uc) || return divrem(p, q)[2]
    n > 0 || return zero(typeof(q))
    c = p.coeff[m-n+2:m+1]
    cp = p.coeff
    cq = q.coeff
    for i = m-n+1:-1:1
        fc = c[n] / uc
        for k = n:-1:2
            c[k] = c[k-1] - cq[k] * fc
        end
        c[1] = cp[i] - cq[1] * fc
    end
    typeof(q)(c)
end

function divrem(p::T, q::T) where T<:UnivariatePolynomial
    cp = p.coeff; cq = q.coeff
    d, r = divrem(cp, cq, Val(false))
    tweak(d, cp, p), tweak(r, cp, p)
end

function tweak(d, cp, p::T) where T<:UnivariatePolynomial
    d === cp ? p : T(d, NOCHECK)
end

function div(p::T, q::T) where T<:UnivariatePolynomial
    cp = p.coeff; cq = q.coeff
    d = divrem(cp, cq, Val(false))[1]
    tweak(d, cp, p)
end

function divrem(f::T, g::AbstractVector{T}) where T<:UnivariatePolynomial
    n = length(g)
    a = zeros(T, n)
    r = copy(f.coeff)
    nf = length(r)
    while nf > 0
        modi = false
        for i = 1:n
            gi = g[i].coeff
            ni = length(gi)
            if ni <= nf
                lcgi = gi[ni]
                d = div(r[nf], lcgi)
                if !iszero(d)
                    modi = true
                    a[i] += d * monom(T, nf-ni)
                    nj = nf - ni
                    for j = nf-ni+1:nf
                        rj = r[j] - gi[j-nf+ni] * d
                        r[j] = rj
                        if !iszero(rj)
                            nj = j
                        end
                    end
                    nf = nj
                end
            end
        end
        if !modi
            nf -= 1
        end
        while nf > 0 && iszero(r[nf])
            nf -= 1
        end
    end
    a, T(r)
end

"""
    q, r, f = pdivrem(a::T, b::T) where T<:UnivariatePolynomial{R}

A variant of divrem, if leading term of divisor is not a unit.
The polynomial `a` is multiplied by a minimal factor `f ∈ R` with
`f * a = q * b + r`.
"""
function pdivrem(p::T, q::T) where {S,T<:UnivariatePolynomial{S}}
    cp = p.coeff; cq = q.coeff
    vd, vr, f = divrem(cp, cq, Val(true))
    tweak(vd, cp, p), tweak(vr, cp, p), f
end

"""
    content(p::Polynomial)

Return the content of the polynomial p, i.e. the `gcd` of its coefficients,
negated if leading coefficient is negative.
"""
function content(p::Polynomial)
    g = gcd(p.coeff)
    isnegative(LC(p)) ? -g : g
end
function content(q::Polynomial{Q}) where Q<:Union{QQ,Quotient}
    c = lcm(getfield.(q.coeff, :den))
    g = gcd(getfield.(q.coeff, :num) .* ( div.(c, getfield.(q.coeff, :den))))
    g = isnegative(LC(q)) ? -g : g
    Q(g, c)
end
"""
    isnegative(a::Ring)

Check if a < 0.
"""
function isnegative(a::RingInt)
    a < zero(a)
end

"""
    primpart(p::Polynomial)

The primitive part of the polynomial p, that means the `gcd` of its coefficients is `1`,
"""
primpart(p::Polynomial) = p / content(p)

#=
Return the degree of the polynomial p, i.e. the highest exponent in the polynomial that has a
nonzero coefficient. The degree of the zero polynomial is defined to be -1.
=#
deg(p::UnivariatePolynomial) = length(p.coeff) - 1

function inv(p::T) where T<:Polynomial
    if isunit(p)
        return T(inv(p[0]))
    else
        throw(DomainError(p, "Only unit polynomials can be inverted"))
    end
end

isunit(p::Polynomial) = deg(p) == 0 && isunit(LC(p))
isone(p::Polynomial) = deg(p) == 0 && isone(LC(p))
iszero(p::Polynomial) = deg(p) < 0
zero(::Type{T}) where {S,T<:UnivariatePolynomial{S}} = T(S[])
one(::Type{T}) where {S,T<:UnivariatePolynomial{S}} = T([one(S)])
==(p::T, q::T) where {T<:UnivariatePolynomial} = p.coeff == q.coeff
function ==(p::S, q::T) where {S<:UnivariatePolynomial,T<:UnivariatePolynomial}
    (varname(S) == varname(T) || deg(p) == 0) && p.coeff == q.coeff
end
==(p::Polynomial, q::Polynomial) = false
function hash(p::UnivariatePolynomial{S,X}, h::UInt) where {X,S}
    n = length(p.coeff)
    n == 0 ? hash(0, h) : n == 1 ? hash(p[0]) : hash(X, hash(p.coeff, h))
end

"""
    ismonom(p::Polynomial)

Return iff polynomial `p` consists is identical to its leading term.
"""
ismonom(p::UnivariatePolynomial) = all(iszero.(view(p.coeff, 1:deg(p))))
"""
    ismonic(p::Polynomial)

Return iff leading coefficient of polynomial `p` is one.
"""
ismonic(p::Polynomial) = isone(LC(p))

# induced homomorphism
function (h::Hom{F,R,S})(p::UnivariatePolynomial{<:R,X}) where {X,F,R,S}
    UnivariatePolynomial{S,X}((h.f).(p.coeff))
end

# auxiliary functions

"""
    LC(p::Polynomial)

Return the leading coefficient of a non-zero polynomial. This coefficient
cannot be zero. Return zero for zero polynomial.
"""
function LC(p::Polynomial)
    c = p.coeff
    n = length(c)
    n == 0 ? zero(basetype(p)) : c[n]
end

function LM(p::P) where {S,P<:UnivariatePolynomial{S}}
    n = length(p.coeff)
    n == 0 && return zero(P)
    coeff = zeros(S, n)
    coeff[n] = one(S)
    P(coeff)
end

"""
    LT(p::Polynomial)

Return leading term of polynomial `p`. Coefficient is taken from `p`.
"""
function LT(p::P) where {S,P<:UnivariatePolynomial{S}}
    n = length(p.coeff)
    n == 0 && return zero(P)
    coeff = zeros(S, n)
    coeff[n] = p.coeff[n]
    P(coeff)
end

"""
    pgcd(a, b)

Pseudo gcd of univariate polynomials gcd`a` and `b`.
Uses subresultant sequence to accomplish non-field coeffient types.
"""
function pgcd(a::T, b::T) where {S,T<:UnivariatePolynomial{S}}
    iszero(b) && return a
    iszero(a) && return b
    a, cc = presultant_seq(a, b, Val(false))
    a = primpart(a)
    a / lcunit(a) * cc
end

"""
    resultant(a, b)

Calculate resultant of two univariate polynomials of general coeffient types.
"""
function resultant(a::T, b::T) where {S,T<:UnivariatePolynomial{S}}
    _, _, r = presultant_seq(a, b, Val(true))
    r(0)
end
resultant(a::T, b::T) where T = iszero(a) || iszero(b) ? zero(T) : oneunit(T)
resultant(a, b) = resultant(promote(a,b)...)

"""
    discriminant(a)

Calculate discriminant of a univariate polynomial `a`.
"""
function discriminant(a::UnivariatePolynomial)
    la = LC(a)
    s = iseven(deg(a) >> 1) ? la : -la
    resultant(a, derive(a)) / s
end

"""
    presultant_seq(a, b)

Modification of Euclid's algorithm to produce `subresultant sequence of pseudo-remainders`.
The next to last calculated remainder is a scalar multiple of the gcd.
Another interpretation of this remainder yields the resultant of `a` and `b`.
See: `https://en.wikipedia.org/wiki/Polynomial_greatest_common_divisor#Subresultant_pseudo-remainder_sequence`
and TAoCP 2nd Ed. 4.6.1.
"""
function presultant_seq(a::T, b::T, ::Val{Usedet})  where {Usedet,S,T<:UnivariatePolynomial{S}}
    E = one(S)
    da = deg(a)
    db = deg(b)
    d = da - db
    s = E
    if d < 0
        a, b, da, db, d = b, a, db, da, -d
        if isodd(da) && isodd(db)
            s = -s
        end
    end
    cc, a, b = normgcd(a, b)
    iszero(b) && return b, cc, b
    s *= cc^(da+db)
    ψ = -E
    β = iseven(d) ? -E : E
    det = isodd(da) && isodd(db) ? -E : E
    while true
        γ = LC(b)
        δ = γ^(d+1)
        a = a * δ
        c = rem(a, b)
        iszero(c) && break
        dc = deg(c)
        c /= β
        # prepare for next turn
        if Usedet
            det = det * δ ^ db / γ ^ (da - dc) / β ^ db

        end
        if isodd(db) && isodd(dc)
            det = -det
        end
        ψ = iszero(d) ? ψ : (-γ)^d / ψ^(d-1)
        a, b = b, c
        da, db = db, dc
        d = da - db
        β = -γ * ψ^d
    end
    r = db > 0 ? zero(b) : d < 1 ? one(b) : d == 1 ? b : LC(b)^(d-1) * b
    b, cc, r * s / det
end

"""
    signed_subresultant_polynomials(P::T, Q::T) where {S,T<:UnivariatePolynomial{S}}

This code is taken from "Algorithm 8.76 (Signed Subresultant Polynomials)"
"Algorithms in Real Algebraic Geometry" by Basu, Pollak, Roy - 2016.
Its use is restricted to the case of `S` is an integral domain (there is no non-trivial divisor of zero).
"""
function signed_subresultant_polynomials(P::T, Q::T)  where {S,T<:UnivariatePolynomial{S}}
    # epsi(n) = (-1) ^ (n*(n-1)÷2)
    epsi(n::Int) = iseven(n >> 1) ? 1 : -1
    p, q = deg(P), deg(Q)
    p > q || throw(ArgumentError("degree of polynomials: $p is not > $q"))
    sresp = zeros(T, p+1)
    s = zeros(S, p+1)
    t = zeros(S, p+1)
    sresp[p+1] = P
    s[p+1] = t[p+1] = 1
    sresp[p] = Q
    bq = LC(Q)
    t[p] = bq
    if p > q - 1
        bqp = bq ^ (p - q - 1)
        sresp[q+1] = (epsi(p-q) * bqp) * Q
    end
    s[q+1] = LC(sresp[q+1])
    i = p + 1
    j = p
    while j > 0 && !iszero(sresp[j])
        k = deg(sresp[j])
        if k == j - 1
            s[j] = t[j]
            if k > 0
                sresp[k] = -rem(s[j]^2 * sresp[i], sresp[j]) / (s[j+1] * t[i])
            end
            elseif k < j - 1
            s[j] = 0
            sig = -1
            for d = 1:j-k-1
                t[j-d] = (t[j]*t[j-d+1]) / s[j+1] * sig
                sig = -sig
            end
            s[k+1] = t[k+1]
            sresp[k+1] = s[k+1] * sresp[j] / t[j]
            if k > 0
                sresp[k] = -rem((t[j] * s[k+1]) * sresp[i], sresp[j]) / (s[j+1] * t[i])
            end
        end
        if k > 0
            t[k] = LC(sresp[k])
        end
        i, j = j, k
    end
    sresp
end

"""
    g, u, v, f = pgcdx(a, b)

Extended pseudo GCD algorithm.
Return `g == pgcd(a, b)` and `u, v, f` with `a * u + b * v == g * f`.
"""
function pgcdx(a::T, b::T) where {X,S,T<:UnivariatePolynomial{S}}
    E = one(S)
    iszero(b) && return a, E, zero(S), a
    iszero(a) && return b, zero(S), E, b
    da = deg(a)
    db = deg(b)
    d = da - db
    if d < 0
        g, u, v, f = pgcdx(b, a)
        return g, v, u, f
    end
    cc, a, b = normgcd(a, b)
    ψ = -E
    β = iseven(d) ? E : -E
    EE = one(T)
    ZZ = zero(T)
    s1, s2, t1, t2 = EE, ZZ, ZZ, EE
    while true
        γ = LC(b)
        γd = γ^(d+1)
        a = a * γd
        q, c = divrem(a, b)
        c /= β
        a, b = b, c
        iszero(b) && break
        s1, s2 = s2, (s1 * γd - s2 * q) / β
        t1, t2 = t2, (t1 * γd - t2 * q) / β
        # prepare for next turn
        da = db
        db = deg(c)
        ψ = d == 0 ? ψ : (-γ)^d / ψ^(d-1)
        d = da - db
        β = -γ * ψ^d
    end
    cs = gcd(content(s2), content(t2))
    a = a / cs
    f = content(a)
    a / (f / cc), s2 / cs, t2 / cs, f
end

"""
    sylvester(u::P, v::P) where P<:UnivariatePolynomial

Return sylvester matrix of polynomials `u` and `v`.
"""
function LinearAlgebra.sylvester(v::P, u::P) where {Z,P<:UnivariatePolynomial{Z}}
    nu = deg(u)
    nv = deg(v)
    n = max(nu, 0) + max(nv, 0)
    S = zeros(Z, n, n)
    if nv >= 0
        for k = 1:nu
            S[k:k+nv,k] .= reverse(v.coeff)
        end
    end
    if nu >= 0
        for k = 1:nv
            S[k:k+nu,k+nu] .= reverse(u.coeff)
        end
    end
    S
end

"""
    resultant_naive(u, v)

Reference implementation of resultant (determinant of sylvester matrix)
"""
function resultant_naive(u::P, v::P) where {Z,P<:UnivariatePolynomial{Z}}
    S = sylvester(u, v)
    det(S)
end

"""
    g, ag, bg = normgcd(a, b)

Divided `a` and `b` by the gcd of their contents.
"""
function normgcd(a, b)
    ca = content(a)
    cb = content(b)
    g = gcd(ca, cb)
    isunit(g) ? (one(g), a, b) : (g, a / g, b / g)
end

"""
    invert(p, q)

Inverse of `p` modulo `q`, where both are polynomials of the same type.
If the base type is not a field, typically no inverse exists.
"""
function invert(p::T, q::T) where T <: UnivariatePolynomial
    g, u, v, f = pgcdx(p, q)
    if isunit(g) && isunit(f)
        u / g / f
    else
        throw(DomainError((p, q), "p is not invertible modulo q"))
    end
end

# multiply p * q with a monom p
function multmono(p, np, vp, q, nq, vq)
    fact = vp[np]
    if isone(fact) && np == 1
        return vq
    end
    v = similar(vq, np + nq - 1)
    z = zero(fact)
    for i = 1:np-1
        v[i] = z
    end
    for i = 1:nq
        v[i+np-1] = vq[i] * fact
    end
    v
end

"""
    evaluate(p, y)

Evaluate polynomial by replacing variable `:x` by `y`. `y` may be an object which
can be converted to `basetype(p)` or another polynomial.
Convenient method ot evaluate is is `p(y)`.
"""
function evaluate(p::UnivariatePolynomial{S}, x::T) where {S,T}
    _evaluate(p, x)
end
function evaluate(p::U, x::V) where {U<:UnivariatePolynomial,Y,T,V<:UnivariatePolynomial{T,Y}}
    if ismonic(x) && ismonom(x)
        spread(p, deg(x), Y, T)
     else
        _evaluate(p, x)
    end
end

function _evaluate(p::UnivariatePolynomial{S}, x::T) where {S,T}
    c = p.coeff
    n = length(c)
    R = promote_type(S,T)
    n == 0 && return zero(R)
    n == 1 && return R(c[1])
    a = convert(R, c[n])
    for k = n-1:-1:1
        a *= x
        a += c[k]
    end
    a
end
(p::UnivariatePolynomial)(a, b...) = evaluate(p, a, b...)

"""
    derive(p::UnivariatePolynomial)

Return formal derivative of polynomial `p`.

For `k in 1:deg(p)` we have `derive(p)[k-1] = k * p[k]`.
If `deg(p) * LC(p) == 0` degree: `deg(derive(p)) < deg(p) - 1`.
"""
function derive(p::P) where P<:UnivariatePolynomial
    n = deg(p)
    n < 0 && return p
    c = similar(p.coeff, n)
    for k = 1:n
        c[k] = p[k] * k
    end
    while n > 0 && iszero(c[n])
        n -= 1
    end
    resize!(c, n)
    P(c, NOCHECK)
end

# efficient implementation of `p(x^m)`.
function spread(p::P, m::Integer) where {T,P<:UnivariatePolynomial{T}}
    P(_spread(p.coeff, m, T))
end
function spread(p::P, m::Integer, Y, S) where {T,P<:UnivariatePolynomial{T}}
    c = _spread(p.coeff, m, S)
    UnivariatePolynomial(Y, c)
end

function _spread(c::Vector{T}, m::Integer, ::Type{S}) where {T,S}
    n = length(c)
    R = promote_type(S, T)
    v = zeros(R, max(0, (n - 1) * m + 1))
    for k = 1:n
        v[(k-1)*m+1] = c[k]
    end
    v
end

function Base.isless(p::T, q::T) where T<:UnivariatePolynomial
    cp = p.coeff
    cq = q.coeff
    lp = length(cp)
    lq = length(cq)
    lp < lq && return true
    lq < lp && return false
    for k = lp:-1:1
        isless(cp[k], cq[k]) && return true
        isless(cq[k], cp[k]) && return false
    end
    false
end

### Display functions

import Base: show

issimple(::Union{ZZ,ZZmod,QQ,Number}) = true
issimple(::Quotient{<:UnivariatePolynomial{S,:γ}}) where S = true
issimple(p::Polynomial) = ismonom(p)
issimple(::Any) = false

function showvar(io::IO, ::UnivariatePolynomial{S,X}, n::Integer) where {X,S}
    if n == 2
        print(io, X)
    elseif n != 1
        print(io, string(X, '^', n-1))
    end
end

isconstterm(::UnivariatePolynomial, n::Integer) = n == 1

function show(io::IO, p::Polynomial{T}) where T
    N = length(p.coeff)
    N <= 0 && return show(io, zero(T))
    start = true
    for n = N:-1:1
        el = p.coeff[n]
        iszero(el) && !start && continue
        !start && print(io, ' ')
        if isconstterm(p, n)
            showelem(io, el, start)
        else
            if isone(el)
                !start && print(io, "+ ")
            elseif isone(-el)
                print(io, '-')
                !start && print(io, ' ')
            elseif !isconstterm(p, n)
                if !issimple(el)
                    !start && print(io, "+ ")
                    print(io, '(')
                    show(io, el)
                    print(io, ')')
                else
                    showelem(io, el, start)
                end
                print(io, '*')
            end
        end
        showvar(io, p, n)
        start = false
    end
end

function showelem(io::IO, el, start::Bool)
    v = sprint(show, el, context=io)
    v1 = v[1]
    if v1 == '-'
        print(io, '-')
        !start && print(io, ' ')
        print(io, SubString(v, nextind(v, 1)))
    else
        start || print(io, "+ ")
        print(io, v1 == '+' ? SubString(v, nextind(v, 1)) : v)
    end
end

LinearAlgebra.det(a::AbstractMatrix{D}) where D<:Ring = det!(copy(a))
det!(a::AbstractMatrix{D}) where D<:Union{Ring,Number}  = _det(a, category_trait(D))
_det(a::AbstractMatrix, ::Type{<:IntegralDomainTrait}) = det_DJB(a)
_det(a::AbstractMatrix{D}, ::Type{<:IntegralDomainTrait}) where D<:QuotientRing = det_DJB(a)
_det(a::AbstractMatrix{D}, ::Type{<:CommutativeRingTrait}) where D<:QuotientRing = det_QR(a)
_det(a::AbstractMatrix, ::Type{<:CommutativeRingTrait}) = det_Bird(a)

"""
    det_DJB(a::Matrix{D}) where D<:Union{Ring,Number}

This code is taken from "Algorithm 8.36 (Dogdson-Jordan-Bareiss)" of
"Algorithms in Real Algebraic Geometry" by Basu, Pollak, Roy - 2016.

Its use is restricted to the case of `D` an integral domain.
"""
function det_DJB(a::AbstractMatrix{D}) where D<:Union{Ring,Number}
    n = checksquare(a)
    n == 0 && return one(D)
    b = copy(a)
    b00 = one(D)
    s = 1
    for k = 1:n-1
        j0 = 0
        for j = k:n
            if !iszero(b[k,j])
                j0 = j
                break
            end
        end
        if j0 == 0
            return zero(D)
        end
        if j0 != k
            for i = k:n
                b[i,j0], b[i,k] = b[i,k], b[i,j0]
            end
            s = -s
        end
        bkk = b[k,k]
        for i = k+1:n
            bik = b[i,k]
            b[i,k] = 0
            for j = k+1:n
                bkj = b[k,j]
                bij = b[i,j]
                b[i,j] = div(bkk * bij - bik * bkj, b00)
            end
        end
        b00 = bkk
    end
    return b[n,n] * s
end

"""
    det_QR(a::Matrix{D}) where D<:QuotientRing

This code is an extension of det_DJB to certain quotient rings which are domains.
I.e. `ZZ/p` where `p` is not prime and `P[:x]/q` where `q` is the product of polynomials.
"""
det_QR(a::AbstractMatrix{D}) where {Z,D<:QuotientRing{Z}} = det_QR!(copy(a))

function det_QR!(b::AbstractMatrix{D}) where {Z,D<:QuotientRing{Z}}
    println("det_QR!($D)")
    n = checksquare(b)
    n == 0 && return one(D)
    s = one(D)
    for k = 1:n-1
        a = view(b, k:n, k:n)
        ij = findfirst(isunit, a)
        # println("k = $k, ij = $ij")
        # display(a)
        if ij !== nothing
            ij += CartesianIndex(k-1, k-1)
            s = flipall!(b, ij, k, s)
        else
            vmin = zero(Z)
            for i = k:n
                v, ij, s = rowdivgcd!(b, i, k, ij, s)
                vmin = max(vmin, v)
            end
            if ij !== nothing
                s = flipall!(b, ij, k, s)
            else
                m = modulus(D)
                v, w = splitmod(vmin, m)
                # println("splitmod($vmin, $m) = $v, $w")
                ZV = Z/v
                dv = det!(ZV.(a))
                isone(w) && return dv
                ZW = Z/w
                dw = det!(ZW.(a))
                return crt(D(value(dv)), D(value(dw)), v, w) * s
            end
        end
        bkk = one(D)
        for i = k+1:n
            bik = b[i,k]
            for j = k+1:n
                bkj = b[k,j]
                bij = b[i,j]
                b[i,j] = bkk * bij - bik * bkj
            end
        end
    end
    return b[n,n] * s
end

function flipall!(a::AbstractMatrix, ij::CartesianIndex, k, s)
    # println("flipall($ij, $k)")
    # display(a)
    i, j = ij[1], ij[2]
    if i != k
        s = -s
        for u = k:size(a, 2)
            a[i,u], a[k,u] = a[k,u], a[i,u]
        end
    end
    if j != k
        s = -s
        for u = k:size(a, 1)
            a[u,j], a[u,k] = a[u,k], a[u,j]
        end
    end
    akk = a[k,k]
    if isunit(akk)
        for u = k:size(a, 2)
            # println("a[$k,$u] / akk == $(a[k,u]) / $akk = $(a[k,u] / akk)")
            a[k,u] /= akk
        end
        s *= akk
    end
    # println("k = $k:")
    # display(a)
    s
end

function rowdivgcd!(b::AbstractMatrix{D}, i, k, ij, s) where {Z,D<:QuotientRing{Z}}
    n = size(b, 2)
    g0 = zero(Z)
    g1 = one(Z)
    for j = k:n
        g1 = gcd(g0, value(b[i,j]))
        isone(g1) && break
        g0 = g1
    end
    if !isone(g1)
        s *= g1
        for j = k:n
            bij = D(value(b[i,j]) ÷ g1)
            b[i,j] = bij
            if ij === nothing && isunit(bij)
                ij = CartesianIndex(i, j)
            end
        end
    end
    g0, ij, s
end

function crt(x, y, p, q)
    g, u, v = gcdx(p, q)
    y * u * p + x * v * q
end
"""
    v = splitmod(a, m)

Assume `gcd(a, m) != 1`.
Find `v` in a way that `gcd(v, m/v) == 1`.
"""
function splitmod(v, m)
    a = gcd(v, m)
    b = m ÷ a
    x = 1

    while !isone(b)
        g = gcd(a, b)
        if !isone(g)
            b ÷= g
            x += 1
        else
            break
        end
        v = g
    end
    a = v ^ x
    b = m ÷ a
    a, b
end

"""
    det_Bird(a::Matrix{D}) where D

This is a division-free algorithm to calculate the determinant of `a`.
There are no conditions on `D` except it is a non-empty commutative ring with one.
Derived from "Pearls of Functional Algorithm Design, chap. 22" by Richard Bird - 2010
"""
function det_Bird(a::AbstractMatrix{D}) where D
    n = checksquare(a)
    n == 0 && return one(D)
    x = copy(a)
    for k = n-1:-1:1
        mutx!(x, k, a)
    end
    isodd(n) ? x[1,1] : -x[1,1]
end

function mutx!(x::AbstractMatrix, k::Integer, a::AbstractMatrix{T}) where T
    # assuming size(x) == size(a) == (n, n) && 0 < k < n
    n = size(a, 1)
    for j = 1:n
        for i = j+1:n
            x[i,j] = 0
        end
    end
    s = -x[k+1,k+1]
    x[k+1,k+1] = 0
    for j = k:-1:1
        xjj = x[j,j]
        x[j,j] = s
        s -= xjj
    end
    s = zero(T)
    for l = k:n
        s += x[k,l] * a[l,k]
    end
    x[k,k] = s
    for j = k+1:n
        x[k,j] = 0
    end
    for i = k-1:-1:1
        for j = i:n
            x[j,1] = x[i,j]
        end
        for j = n:-1:1
            s = zero(T)
            for l = i:n
                s += x[l,1] * a[l,j]
            end
            x[i,j] = s
        end
        for j = i+1:n
            x[j,1] = 0
        end
    end
    x
end
"""
    det3(a)

Divisionfree combinatoric algorithm to calculate determinant.
"Determinant: Combinatorics, Algorithms, andComplexity" by Mahajan, Vinay
in Chicago Journal of Theoretical Computer Science1997-5
http://cjtcs.cs.uchicago.edu/articles/1997/5/cj97-05.pdf
"""
function det3_new(a::AbstractMatrix{D}) where D<:Union{Ring,Integer}
    n = checksquare(a)
    A = Vector{Tuple{Bool,Int,Int,D}}(undef, 0)
    B = Vector{Tuple{Bool,Int,Int,D}}(undef, 0)
    C = Dict{Tuple,Int}()

    function addpush!(B, C, t)
        #println("addpush!($t)")
        p, u, v, d = t
        key = (p, u, v)
        nx = length(B) + 1
        ix = get!(C, key, nx)
        if ix == nx
            push!(B, t)
        else
            _, _, _, b = B[ix]
            B[ix] = (p, u, v, b + d)
        end
        nothing
    end

    s = max(min(n * (n+1) ÷ 2, count(!iszero, a)), max(2 * n, 32))
    sizehint!(A, s)
    sizehint!(B, s)
    p = isodd(n)
    for u = 1:n
        push!(A, (p, u, u, one(D)))
    end
    for i = 0:n-2
        empty!(B)
        empty!(C)
        for (p, u, v, Auvp) in A
            avu = a[v,u]
            for w = u+1:n
                avw = a[v,w]
                if !iszero(avw)
                    avw *= Auvp
                    addpush!(B, C, (p, u, w, avw))
                end
            end
            if !iszero(avu)
                avu *= Auvp
                for w = u+1:n
                    addpush!(B, C, (!p, w, w, avu))
                end
            end
        end
        A, B = B, A
        #println("i = $i")
        #display(A)
    end

    d = zero(D)
    for (p, u, v, Auvp) in A
        avu = a[v,u] * (p ? 1 : -1)
        if !iszero(avu)
            d += Auvp * avu
        end
    end
    d
end

function det3(a::AbstractMatrix{D}) where D<:Union{Ring,Integer}
    n = checksquare(a)
    A1 = zeros(D, n, n)
    A2 = zeros(D, n, n)
    B1 = zeros(D, n, n)
    B2 = zeros(D, n, n)
    p = n & 1 + 1
    for u = 1:n
        (p == 1 ? A1 : A2)[u,u] = one(D)
    end
    for i = 0:n-2
        for u = 1:n
            for v = u:n
                Auvp = A1[u,v]
                if !iszero(Auvp)
                    for w = u+1:n
                        #println("addpush(0, $u, $w, $(Auvp*a[v,w]))")
                        B1[u,w] += Auvp * a[v,w]
                        #println("addpush(1, $w, $w, $(Auvp * a[v,u]))")
                        B2[w,w] += Auvp * a[v,u]
                    end
                end
                Auvp = A2[u,v]
                if !iszero(Auvp)
                    for w = u+1:n
                        #println("addpush(1, $u, $w, $(Auvp*a[v,w]))")
                        B2[u,w] += Auvp * a[v,w]
                        #println("addpush(0, $w, $w, $(Auvp * a[v,u]))")
                        B1[w,w] += Auvp * a[v,u]
                    end
                end
            end
        end
        A1, B1 = B1, A1
        A2, B2 = B2, A2
        fill!(B1, zero(D))
        fill!(B2, zero(D))
        #println("i = $i")
        #display([A1 A2])
    end

    for u = 1:n
        for v = u:n
            avu = a[v,u]
            if !iszero(avu)
                B1[1,1] += A1[u,v] * avu
                B2[1,1] += A2[u,v] * avu
            end
        end
    end
    B2[1,1] - B1[1,1]
end

function hamilton_normal_form(a::Matrix{R}) where R<:Union{Ring,Integer}
    m, n = size(a)
    u = Matrix(R.(I(n)))
    hamilton_normal_form!(copy(a), u)
end

function hamilton_normal_form!(a::Matrix, u::Matrix)
    m, n = size(a)

    for i = 1:min(m,n)-1
        for j = 1:i
            ajj = a[j,j]
            aji = a[j,i+1]
            r, p, q = gcdx(ajj, aji)
            pp = -div(aji, r)
            qq = div(ajj, r)
            for k = 1:m
                akj = a[k,j]
                aki = a[k,i+1]
                bkj = akj * p + aki * q
                bki = akj * pp + aki * qq
                a[k,j] = bkj
                a[k,i+1] = bki
            end
            for k = 1:n
                akj = u[k,j]
                aki = u[k,i+1]
                bkj = akj * p + aki * q
                bki = akj * pp + aki * qq
                u[k,j] = bkj
                u[k,i+1] = bki
            end
            reduce_off_diagonal!(a, u, j)
        end
        reduce_off_diagonal!(a, u, i + 1)
    end
    a, u
end

function reduce_off_diagonal!(a, u, k)
    akk = a[k,k]
    if akk < 0
        akk = -akk
        a[:,k] .= -a[:,k]
        u[:,k] .= -u[:,k]
    end
    for z = 1:k-1
        d = div(-a[k,z], akk)
        a[:,z] .+= a[:,k] .* d
        u[:,z] .+= u[:,k] .* d
    end
end


text:
using Test
using DataStructure101
const DS = DataStructure101

@testset "Min Max Heap" begin
    heap = DS.MinMaxHeap()
    @test heap isa DS.MinMaxHeap
    @test heap isa DS.Heap

    for i=30:-1:1
        push!(heap, i)
    end
    DS.build!(heap)

    @test maximum(heap) == 30
    @test minimum(heap) == 1

    @test DS.pop_max!(heap) == 30
    @test DS.pop_min!(heap) == 1

    DS.push_bubble!(heap, 1)
    @test minimum(heap) == 1

    DS.push_bubble!(heap, 30)
    @test maximum(heap) == 30

    @test sort!(heap) == collect(1:30)

    for i=1:30
        push!(heap, i)
    end
    DS.build!(heap)

    @test maximum(heap) == 30
    @test minimum(heap) == 1

    @test DS.pop_max!(heap) == 30
    @test DS.pop_min!(heap) == 1

    DS.push_bubble!(heap, 1)
    @test minimum(heap) == 1

    DS.push_bubble!(heap, 30)
    @test maximum(heap) == 30

    @test sort!(heap) == collect(1:30)

    for i=1:30
        DS.push_bubble!(heap, i)
    end

    @test maximum(heap) == 30
    @test minimum(heap) == 1
    @test sort!(heap) == collect(1:30)

    for i=30:-1:1
        DS.push_bubble!(heap, i)
    end

    @test maximum(heap) == 30
    @test minimum(heap) == 1
    @test sort!(heap) == collect(1:30)
end

ext: jl


text:
using JMPS,LinearAlgebra,OMEinsum
using Test
using CUDA
CUDA.allowscalar(false)

@testset "Constructor" begin
    mps = MPS(5,2,4,FloatType=Float64,ArrayType=Array)
    for i = 1:5
        @test  norm(mps[i]) ≈ 0
    end
    mpo = MPO(5,2,4,FloatType=ComplexF64,ArrayType=Array)
    for i = 1:5
        @test  norm(mpo[i]) ≈ 0
    end
    mps = MPS(5,2,4,FloatType=Float64,ArrayType=Array)
    for i = 1:5
        @test  norm(mps[i]) ≈ 0
    end
    mpo = MPO(5,2,4,FloatType=ComplexF64,ArrayType=Array)
    for i = 1:5
        @test  norm(mpo[i]) ≈ 0
    end
    mps = MPS(5,2,4,FloatType=Float64,ArrayType=CuArray)
    for i = 1:5
        @test  norm(mps[i]) ≈ 0
    end
    mpo = MPO(5,2,4,FloatType=ComplexF64,ArrayType=CuArray)
    for i = 1:5
        @test  norm(mpo[i]) ≈ 0
    end
    mps = MPS(5,2,4,FloatType=Float64,ArrayType=CuArray)
    for i = 1:5
        @test  norm(mps[i]) ≈ 0
    end
    mpo = MPO(5,2,4,FloatType=ComplexF64,ArrayType=CuArray)
    for i = 1:5
        @test  norm(mpo[i]) ≈ 0
    end
end

@testset "get-set index" begin
    mps = randomMPS(5,2,4,FloatType=Float64,ArrayType=Array)
    mps[3] .= 0
    @test norm(mps[3])≈0
    mpo = randomMPO(3,3,3,FloatType=Float64,ArrayType=Array)
    mpo[2] .= 0
    @test norm(mpo[2])≈0
end

@testset "Bdims" begin
    bdim = JMPS.Bdims(4,[2;3;4;5])
    @test bdim[0] == 5
    @test bdim[1] == 2
    @test bdim[2] == 3
    @test bdim[3] == 4
    @test bdim[4] == 5
    @test bdim[5] == 2
end

@testset "single_tensor_multiply" begin
    mpo = rand(4,5,6,7)
    mps = rand(3,7,4)
    @test norm(JMPS.single_tensor_multiply(mpo,mps)-ein"lurd,adb->laurb"(mpo, mps))≈0
end

@testset "Adjoint MPS" begin
    mps1 = randomMPS(4,4,31)
    mps2 = randomMPS(4,4,31)
    @test (mps1*mps2-mps2*mps1)≈0
end

@testset "normalization and overlap" begin
    mps = randomMPS(4,4,51)
    normalization!(mps,LeftNormalization())
    mps1 = deepcopy(mps)
    mps2 = mps
    compress!(mps2,50)
    ovlp = mps2*mps1
    n1 = mps1*mps1
    n2 = mps2*mps2
    @test n1 - n2 ≈ 0 atol = 3E-7
    @test ovlp - n2 ≈ 0 atol = 3E-7
    @test ovlp - n1 ≈ 0 atol = 3E-7
end

@testset "N=2 AKLT chain" begin
    N = 2
    σ⁺ = [0 1;0 0 ]
    σ⁻ = [0 0;1 0 ]
    σᶻ = [1 0;0 -1]
    mps = MPS(N,3,2)
    mps[1][1,1,:] .= ( sqrt(2/3)σ⁺)[1,:]
    mps[1][1,2,:] .= (-sqrt(1/3)σᶻ)[1,:]
    mps[1][1,3,:] .= (-sqrt(2/3)σ⁻)[1,:]
    mps[N][:,1,1] .= ( sqrt(2/3)σ⁺)[:,1]
    mps[N][:,2,1] .= (-sqrt(1/3)σᶻ)[:,1]
    mps[N][:,3,1] .= (-sqrt(2/3)σ⁻)[:,1]
    @test disp(mps,[1;3])≈-2/3
    @test disp(mps,[3;1])≈0
    @test disp(mps,[2;2])≈1/3
    @test mps*mps ≈ 5/9
    x = 1/5
    @test entropy(mps,1) ≈ -(x*log(x)+4x*log(4x)) #will normalize mps
    @test sum(entropy(mps)) ≈ -(x*log(x)+4x*log(4x))
    normalization!(mps,LeftNormalization())
    @test disp(mps,[1;3])≈-2/sqrt(5)
    @test disp(mps,[2;2])≈1/sqrt(5)
end

max_stars_repo_name: qiyang-ustc/JMPS

max_issues_repo_name: qiyang-ustc/JMPS

max_forks_repo_name: qiyang-ustc/JMPS


text:
type ODESolution{uType,uType2,DType,tType,rateType,P,A,IType} <: AbstractODESolution
  u::uType
  u_analytic::uType2
  errors::DType
  t::tType
  k::rateType
  prob::P
  alg::A
  interp::IType
  dense::Bool
  tslocation::Int
  retcode::Symbol
end
(sol::ODESolution)(t,deriv::Type=Val{0};idxs=nothing) = sol.interp(t,idxs,deriv)
(sol::ODESolution)(v,t,deriv::Type=Val{0};idxs=nothing) = sol.interp(v,t,idxs,deriv)

function build_solution{uType,tType,isinplace}(
        prob::Union{AbstractODEProblem{uType,tType,isinplace},AbstractDDEProblem{uType,tType,isinplace}},
        alg,t,u;dense=false,timeseries_errors=true,dense_errors=true,
        calculate_error = true,
        k=[],interp = (tvals) -> nothing, retcode = :Default, kwargs...)

  if has_analytic(prob.f)
    u_analytic = Vector{uType}(0)
    errors = Dict{Symbol,eltype(prob.u0)}()
    sol = ODESolution(u,u_analytic,errors,t,k,prob,alg,interp,dense,0,retcode)
    if calculate_error
      calculate_solution_errors!(sol;timeseries_errors=timeseries_errors,dense_errors=dense_errors)
    end
    return sol
  else
    return ODESolution(u,nothing,nothing,t,k,prob,alg,interp,dense,0,retcode)
  end
end

function calculate_solution_errors!(sol::AbstractODESolution;fill_uanalytic=true,timeseries_errors=true,dense_errors=true)
  if fill_uanalytic
    for i in 1:size(sol.u,1)
      push!(sol.u_analytic,sol.prob.f(Val{:analytic},sol.t[i],sol.prob.u0))
    end
  end

  save_everystep = length(sol.u) > 2
  if !isempty(sol.u_analytic)
    sol.errors[:final] = mean(abs.(sol.u[end]-sol.u_analytic[end]))

    if save_everystep && timeseries_errors
      sol.errors[:l∞] = maximum(vecvecapply((x)->abs.(x),sol.u-sol.u_analytic))
      sol.errors[:l2] = sqrt(mean(vecvecapply((x)->float.(x).^2,sol.u-sol.u_analytic)))
      if sol.dense && dense_errors
        densetimes = collect(linspace(sol.t[1],sol.t[end],100))
        interp_u = sol(densetimes)
        interp_analytic = [sol.prob.f(Val{:analytic},t,sol.u[1]) for t in densetimes]
        sol.errors[:L∞] = maximum(vecvecapply((x)->abs.(x),interp_u-interp_analytic))
        sol.errors[:L2] = sqrt(mean(vecvecapply((x)->float.(x).^2,interp_u-interp_analytic)))
      end
    end
  end
end

function build_solution(sol::AbstractODESolution,u_analytic,errors)
  ODESolution(sol.u,u_analytic,errors,sol.t,sol.k,sol.prob,sol.alg,sol.interp,sol.dense,sol.tslocation,sol.retcode)
end


text:
##############################################################################
##
## get weight
##
##############################################################################

function get_weight(df::AbstractDataFrame, weight::Symbol)
    convert(Vector{Float64}, sqrt(df[weight]))
end
function get_weight(df::AbstractDataFrame, ::Void)
    Ones(size(df, 1))
end

##############################################################################
##
## build model
##
##############################################################################

function reftype(sz)
    sz <= typemax(UInt8)  ? UInt8 :
    sz <= typemax(UInt16) ? UInt16 :
    sz <= typemax(UInt32) ? UInt32 :
    UInt64
end
function simpleModelFrame(df, t, esample)
    df1 = DataFrame(map(x -> df[x], t.eterms))
    names!(df1, convert(Vector{Symbol}, map(string, t.eterms)))
    mf = ModelFrame(df1, t, esample)
end

#  remove observations with negative weights
function isnaorneg{T <: Real}(a::Vector{T})
    BitArray(a .> zero(eltype(a)))
end
function isnaorneg{T <: Real}(a::DataVector{T})
    out = !a.na
    @inbounds @simd for i in 1:length(a)
        if out[i]
            out[i] = a[i] > zero(Float64)
        end
    end
    BitArray(out)
end

# Directly from DataFrames.jl
function remove_response(t::Terms)
    # shallow copy original terms
    t = Terms(t.terms, t.eterms, t.factors, t.order, t.response, t.intercept)
    if t.response
        t.order = t.order[2:end]
        t.eterms = t.eterms[2:end]
        t.factors = t.factors[2:end, 2:end]
        t.response = false
    end
    return t
end

# used when removing certain rows in a dataset
# NA always removed
function dropUnusedLevels!(f::PooledDataVector)
    uu = unique(f.refs)
    length(uu) == length(f.pool) && return f
    sort!(uu)
    T = reftype(length(uu))
    dict = Dict{eltype(uu), T}(zip(uu, collect(1:convert(T, length(uu)))))
    @inbounds @simd for i in 1:length(f.refs)
         f.refs[i] = dict[f.refs[i]]
    end
    f.pool = f.pool[uu]
    f
end

dropUnusedLevels!(f::DataVector) = f

##############################################################################
##
## sum of squares
##
##############################################################################

function compute_tss(y::Vector{Float64}, hasintercept::Bool, ::Ones)
    if hasintercept
        tss = zero(Float64)
        m = mean(y)::Float64
        @inbounds @simd  for i in 1:length(y)
            tss += abs2((y[i] - m))
        end
    else
        tss = sumabs2(y)
    end
    return tss
end

function compute_tss(y::Vector{Float64}, hasintercept::Bool, sqrtw::Vector{Float64})
    if hasintercept
        m = (mean(y) / sum(sqrtw) * length(y))::Float64
        tss = zero(Float64)
        @inbounds @simd  for i in 1:length(y)
         tss += abs2(y[i] - sqrtw[i] * m)
        end
    else
        tss = sumabs2(y)
    end
    return tss
end


text:
##################################################################################
# This file is part of StateSpaceEcon.jl
# BSD 3-Clause License
# Copyright (c) 2020-2022, Bank of Canada
# All rights reserved.
##################################################################################

export shockdecomp
import ..steadystatedata

"""
    shockdecomp(model, plan, exog_data; control, fctype, [options])

Compute the shock decomposition for the given model, plan, exogenous (shocks)
data and control solution.

If `control` option is not specified we use the steady state solution stored in
the model instance. The algorithm assumes that `control` is a solution to the
dynamic model for the given plan range and final condition. We verify the
residual and issue a warning, but do not enforce this. See
[`steadystatedata`](@ref).

As part of the algorithm we run a simulation with the given `plan`,
`exog_data` and `fctype`.  See [`simulate`](@ref) for other options.

!!! note

    For now only the case of `anticipate=true` works. Shock decomp with
    unanticipated shocks is coming soon.

"""
function shockdecomp(m::Model, p::Plan, exog_data::SimData;
    control::SimData = steadystatedata(m, p),
    deviation::Bool = false,
    anticipate::Bool = true,
    verbose::Bool = getoption(m, :verbose, false),
    maxiter::Int = getoption(m, :maxiter, 20),
    tol::Float64 = getoption(m, :tol, 1e-12),
    fctype::FinalCondition = getoption(m, :fctype, fcgiven),
    _debug = false
)

    if !anticipate
        error("Not yet implemented with `anticipate=false`. Try again later.")
    end

    # if m.nauxs > 0
    #     error("Not yet implemented with auxiliary variables.")
    # end

    # we will return result at the end
    result = Workspace()
    result.c = copy(control)

    # apply data transformations of @log variables
    exog_data = transform(exog_data[p.range, :], m)
    control = transform(control[p.range, :], m)

    if deviation
        # with transformed data we simply add the control
        exog_data .+= control
    end

    # prepare the stacked-time
    gdata = StackedTimeSolverData(m, p, fctype)

    # check the residual.           why are we doing this? we know it's 0!
    shocked = copy(exog_data)
    res_shocked = Vector{Float64}(undef, size(gdata.J, 1))
    global_R!(res_shocked, shocked, shocked, gdata)    # Run the "shocked" simulation with the given exogenous data.
    if norm(res_shocked, Inf) > tol
        assign_exog_data!(shocked, exog_data, gdata)
        sim_nr!(shocked, gdata, maxiter, tol, verbose)
    end

    # We need the Jacobian matrix evaluated at the control solution.
    res_control, _ = global_RJ(control, control, gdata)
    if norm(res_control, Inf) > tol
        # What to do if it's not a solution? We just give a warning for now, but maybe we should error()?!
        @warn "Control is not a solution:" norm(res_control, Inf)
    end

    # the difference between shocked and control solutions
    delta = shocked - control

    # now we need to split the unknowns into groups.
    begin
        LI = LinearIndices((p.range, 1:gdata.NU))
        NT = length(p.range)
        # range of initial conditions
        init = 1:m.maxlag
        # range of final conditions
        term = NT-m.maxlead+1:NT
        # simulation range
        sim = m.maxlag+1:NT-m.maxlead

        # indices of endogenous variable-points (the unknowns we were solving for)
        endo_inds = vec(vcat(
            (LI[sim, v.index] for v in m.allvars if !(isexog(v) || isshock(v)))...
        ))

        if length(endo_inds) != sum(gdata.solve_mask) || !all(gdata.solve_mask[endo_inds])
            error("Not yet implemented with non-trivial plan.")
        end

        # indices for the exogenous variable-points split into groups
        exog_inds = Workspace(;
            # contributions of initial conditions
            init = vec(LI[init, :]),
            # contributions of final conditions
            term = vec(LI[term, :]),
            # contributions of shocks
            (v.name => vec(LI[sim, v.index]) for v in m.allvars if isshock(v))...,
            # contributions of @exog variables
            (v.name => vec(LI[sim, v.index]) for v in m.allvars if isexog(v))...
        )
        # Note: we use a Workspace() because it preserves the order in which
        #     members were added to it

        if sum(length, values(exog_inds)) + length(endo_inds) != size(gdata.J, 2)
            error("endogenous and exogenous don't add up!")
        end
    end

    # now do the decomposition
    begin
        # Allocate the shock-decomposition matrix
        SDMAT = zeros(size(gdata.J, 1), length(exog_inds))
        # Fill it with the right-hand side of the system
        for (col, inds) in enumerate(values(exog_inds))
            SDMAT[:, [col]] .-= sum(gdata.J[:, inds] .* delta[inds]', dims = 2)
        end
        # solve the system (in-place)
        ldiv!(gdata.luJ[], SDMAT)
    end

    # now split and decorate the shock-decomposition matrix
    result.sd = Workspace()
    begin
        # in order to split the rows of SDMAT by variable, we need the inverse indexing map
        inv_endo_inds = zeros(Int, size(gdata.solve_mask))
        inv_endo_inds[gdata.solve_mask] .= 1:sum(gdata.solve_mask)
        for v in m.allvars
            v_inds = vec(LI[:, v.index])
            v_endo_mask = gdata.solve_mask[v_inds]
            if !any(v_endo_mask)
                continue
            end
            v_inv_endo_inds = inv_endo_inds[v_inds[v_endo_mask]]
            v_data = MVTSeries(p.range, keys(exog_inds), zeros)
            v_data[v_endo_mask, :] .= SDMAT[v_inv_endo_inds, :]
            # we also append the truncation error due to linear approximation of the decomposition.
            v_data = hcat(v_data; nonlinear = delta[v.name] - sum(v_data, dims = 2))
            push!(result.sd, v.name => v_data)
        end
    end

    result.s = copy(result.c)
    result.s[p.range, :] .= shocked

    if deviation
        logvars = islog.(m.varshks) .| isneglog.(m.varshks)
        result.s[:, logvars] ./= result.c[:, logvars]
        result.s[:, .!logvars] .-= result.c[:, .!logvars]
    end

    if _debug
        return result, gdata
    else
        return result
    end
end


text:
immutable ContinuationFrame
    is_end_headers::Bool
    stream_identifier::UInt32
    fragment::Array{UInt8, 1}
end

==(a::ContinuationFrame, b::ContinuationFrame) =
    a.is_end_headers == b.is_end_headers &&
    a.stream_identifier == b.stream_identifier &&
    a.fragment == b.fragment

function decode_continuation(header, payload)
    is_end_headers = header.flags & 0x4 == 0x4

    return ContinuationFrame(is_end_headers, header.stream_identifier, payload)
end

function encode_continuation(frame)
    flags = 0x0 | (frame.is_end_headers ? 0x4 : 0x0)

    return wrap_payload(frame.fragment, CONTINUATION, flags, frame.stream_identifier)
end


text:
#!/usr/bin/env julia

using StatsBase

"""
This function gives the neighbors of a cell. Takes into
account the diagonal ones.
"""
function neighbors(population,i,j)
    above = [(i-1,k) for k in j-1:j+1 if (k>0 && i-1>0 && k<=size(population)[2])]
    below = [(i+1,k) for k in j-1:j+1 if (k>0 && i<size(population)[1] && k<=size(population)[2])]
    right = [(i,j+1) for i in i if (j<size(population)[2])]
    left  = [(i,j-1) for i in i if (j-1 > 0)]
    neigh = union(above,below,right,left)
    neigh
end

"""
This function initializes the opinion population given
a p probability of 1 and p-1 of -1
"""
function initialize_population(x,y,p)
    sizes = Int(ceil(x*y*(1-p)))
    init = ones(Int64,x,y)
    coord = [tuple(i,j) for i in 1:x for j in 1:y]
    selected = sample(coord,sizes;replace=false)
    map(x -> init[x[1],x[2]]=-init[x[1],x[2]],selected)
    init
end

function friends_opinion(population,i,j)
    n = neighbors(population,i,j)
    op = [population[n[k][1],n[k][2]] for k in 1:length(n)]
    op
end

function avg_friends_opinion(population,i,j)
    op = friends_opinion(population,i,j)
    sum(op)/length(op)
end

function convince(population,i,j)
    my_opinion = population[i,j]
    neighborhood = neighbors(population,i,j)
    map(x-> population[x[1],x[2]]=my_opinion,neighborhood)
    population
end

function convince_2_sites(population,i,j)
    my_opinion = population[i,j]
    neighborhood = neighbors(population,i,j)
    rf = sample(neighborhood,1)
    if population[rf[1][1],rf[1][2]] == my_opinion
        map(x-> population[x[1],x[2]]=my_opinion,neighborhood)
    end
    population
end

function convince_3_sites(population,i,j)
    my_opinion = population[i,j]
    neighborhood = neighbors(population,i,j)
    rf = sample(neighborhood,2)
    if population[rf[1][1],rf[1][2]] === population[rf[2][1],rf[2][2]] == my_opinion
        map(x-> population[x[1],x[2]]=my_opinion,neighborhood)
    end
    population
end

function sznajd_one_site(population,p)
    sizes = Int(ceil(length(population)*p))
    coord = [tuple(i,j) for i in 1:size(population)[1] for j in 1:size(population)[2]]
    selected = sample(coord,sizes;replace=false)
    map(x -> convince(population,x[1],x[2]),selected)
    population
end

function sznajd_two_site(population,p)
    sizes = Int(ceil(length(population)*p))
    coord = [tuple(i,j) for i in 1:size(population)[1] for j in 1:size(population)[2]]
    selected = sample(coord,sizes;replace=false)
    map(x -> convince_2_sites(population,x[1],x[2]),selected)
    population
end

function sznajd_three_site(population,p)
    sizes = Int(ceil(length(population)*p))
    coord = [tuple(i,j) for i in 1:size(population)[1] for j in 1:size(population)[2]]
    selected = sample(coord,sizes;replace=false)
    map(x -> convince_3_sites(population,x[1],x[2]),selected)
    population
end

function stupid_bots(population,p)
    coord = [tuple(i,j) for i in 1:size(population)[1] for j in 1:size(population)[2]]
    hm = Int(ceil(length(population)*p))
    selected = sample(coord,hm;replace=false)
    map(x -> population[x[1],x[2]]=1,selected)
    population
end

function smart_bots(population,p)
    coord = [tuple(i,j) for i in 1:size(population)[1] for j in 1:size(population)[2]
            if population[i,j]!=1]
    hm = Int(ceil(length(coord)*p))
    selected = sample(coord,hm;replace=false)
    map(x -> population[x[1],x[2]]=1,selected)
    population
end

function cambridge_analytica(population,p)
    coord = [tuple(i,j) for i in 1:size(population)[1] for j in 1:size(population)[2] if (population[i,j]!=1 && abs(avg_friends_opinion(population,i,j))!=1)]
    hm = Int(ceil(length(coord)*p))
    selected = sample(coord,hm;replace=true)
    map(x -> convince_2_site(population,x[1],x[2]),selected)
    population
end


text:
using OrthogonalPolynomialsQuasi, ContinuumArrays, DomainSets, FillArrays, Test
import OrthogonalPolynomialsQuasi: jacobimatrix

@testset "Hermite" begin
    H = Hermite()
    @test axes(H) == (Inclusion(ℝ), Base.OneTo(∞))
    x = axes(H,1)
    X = jacobimatrix(H)

    w = HermiteWeight()
    wH = w.*H
    M = H'* ( w.*H)
    S = Diagonal(M.diag .^ (-1/2))
    Si = Diagonal(M.diag .^ (1/2))
    J = Si*X*S

    (J - im*Eye(∞)) \ [1;zeros(∞)]

    @test H[0.1,1] === 1.0 # equivalent to H_0(0.1) == 1.0
    D = Derivative(x)

    h = 0.000001
    @test (D*H)[0.1,1:5] ≈ (H[0.1+h,1:5] - H[0.1,1:5])/h atol=100h
end


text:
using Downloads: download
using OrdinaryDiffEq
using Trixi

###############################################################################
# semidiscretization of the shallow water equations with a continuous
# bottom topography function (set in the initial conditions)

equations = ShallowWaterEquations2D(gravity_constant=1.0, H0=3.0)

# An initial condition with constant total water height and zero velocities to test well-balancedness.
function initial_condition_well_balancedness(x, t, equations::ShallowWaterEquations2D)
  # Set the background values
  H = equations.H0
  v1 = 0.0
  v2 = 0.0
  # bottom topography taken from Pond.control in [HOHQMesh](https://github.com/trixi-framework/HOHQMesh)
  x1, x2 = x
  b = (  1.5 / exp( 0.5 * ((x1 - 1.0)^2 + (x2 - 1.0)^2) )
       + 0.75 / exp( 0.5 * ((x1 + 1.0)^2 + (x2 + 1.0)^2) ) )
  return prim2cons(SVector(H, v1, v2, b), equations)
end

initial_condition = initial_condition_well_balancedness

boundary_condition_constant = BoundaryConditionDirichlet(initial_condition)
boundary_condition = Dict( :OuterCircle => boundary_condition_constant )

###############################################################################
# Get the DG approximation space

volume_flux = (flux_wintermeyer_etal, flux_nonconservative_wintermeyer_etal)
solver = DGSEM(polydeg=4, surface_flux=(flux_hll, flux_nonconservative_fjordholm_etal),
               volume_integral=VolumeIntegralFluxDifferencing(volume_flux))

###############################################################################
# This setup is for the curved, split form well-balancedness testing

# Get the unstructured quad mesh from a file (downloads the file if not available locally)
default_mesh_file = joinpath(@__DIR__, "mesh_outer_circle.mesh")
isfile(default_mesh_file) || download("https://gist.githubusercontent.com/andrewwinters5000/9beddd9cd00e2a0a15865129eeb24928/raw/be71e67fa48bc4e1e97f5f6cd77c3ed34c6ba9be/mesh_outer_circle.mesh",
                                      default_mesh_file)
mesh_file = default_mesh_file

mesh = UnstructuredMesh2D(mesh_file)

# Create the semi discretization object
semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
                                    boundary_conditions=boundary_condition)

###############################################################################
# ODE solvers, callbacks, etc.

tspan = (0.0, 100.0)
ode = semidiscretize(semi, tspan)

summary_callback = SummaryCallback()

analysis_interval = 1000
analysis_callback = AnalysisCallback(semi, interval=analysis_interval,
                                     save_analysis=true,
                                     extra_analysis_integrals=(lake_at_rest_error,))

alive_callback = AliveCallback(analysis_interval=analysis_interval)

save_solution = SaveSolutionCallback(interval=1000,
                                     save_initial_solution=true,
                                     save_final_solution=true)

callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution)

###############################################################################
# run the simulation

# use a Runge-Kutta method with automatic (error based) time step size control
sol = solve(ode, RDPK3SpFSAL49(), abstol=1.0e-11, reltol=1.0e-11,
            save_everystep=false, callback=callbacks);
summary_callback() # print the timer summary


text:
# Autogenerated wrapper script for QuantumEspresso_jll for x86_64-linux-musl-libgfortran3
export bands, carparinello, density_of_states, dynamical_matrix_gamma, dynamical_matrix_generic, hubbardparams, ibrav_to_cell, kpoints, nudged_elastic_band, phonon, plotband, plotrho, pwscf, reciprocal_to_real

using CompilerSupportLibraries_jll
using FFTW_jll
using Libxc_jll
using MPICH_jll
using OpenBLAS32_jll
using SCALAPACK_jll
JLLWrappers.@generate_wrapper_header("QuantumEspresso")
JLLWrappers.@declare_executable_product(bands)
JLLWrappers.@declare_executable_product(carparinello)
JLLWrappers.@declare_executable_product(density_of_states)
JLLWrappers.@declare_executable_product(dynamical_matrix_gamma)
JLLWrappers.@declare_executable_product(dynamical_matrix_generic)
JLLWrappers.@declare_executable_product(hubbardparams)
JLLWrappers.@declare_executable_product(ibrav_to_cell)
JLLWrappers.@declare_executable_product(kpoints)
JLLWrappers.@declare_executable_product(nudged_elastic_band)
JLLWrappers.@declare_executable_product(phonon)
JLLWrappers.@declare_executable_product(plotband)
JLLWrappers.@declare_executable_product(plotrho)
JLLWrappers.@declare_executable_product(pwscf)
JLLWrappers.@declare_executable_product(reciprocal_to_real)
function __init__()
    JLLWrappers.@generate_init_header(CompilerSupportLibraries_jll, FFTW_jll, Libxc_jll, MPICH_jll, OpenBLAS32_jll, SCALAPACK_jll)
    JLLWrappers.@init_executable_product(
        bands,
        "bin/bands.x",
    )

    JLLWrappers.@init_executable_product(
        carparinello,
        "bin/cp.x",
    )

    JLLWrappers.@init_executable_product(
        density_of_states,
        "bin/dos.x",
    )

    JLLWrappers.@init_executable_product(
        dynamical_matrix_gamma,
        "bin/dynmat.x",
    )

    JLLWrappers.@init_executable_product(
        dynamical_matrix_generic,
        "bin/matdyn.x",
    )

    JLLWrappers.@init_executable_product(
        hubbardparams,
        "bin/hp.x",
    )

    JLLWrappers.@init_executable_product(
        ibrav_to_cell,
        "bin/ibrav2cell.x",
    )

    JLLWrappers.@init_executable_product(
        kpoints,
        "bin/kpoints.x",
    )

    JLLWrappers.@init_executable_product(
        nudged_elastic_band,
        "bin/neb.x",
    )

    JLLWrappers.@init_executable_product(
        phonon,
        "bin/ph.x",
    )

    JLLWrappers.@init_executable_product(
        plotband,
        "bin/plotband.x",
    )

    JLLWrappers.@init_executable_product(
        plotrho,
        "bin/plotrho.x",
    )

    JLLWrappers.@init_executable_product(
        pwscf,
        "bin/pw.x",
    )

    JLLWrappers.@init_executable_product(
        reciprocal_to_real,
        "bin/q2r.x",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
##############################################################################
##
## Weight
##
##############################################################################
function get_weights(df::AbstractDataFrame, esample::AbstractVector, weights::Symbol)
    # there are no NA in it. DataVector to Vector
    out = convert(Vector{Float64}, df[esample, weights])
    map!(sqrt, out, out)
    return out
end
get_weights(df::AbstractDataFrame, esample::AbstractVector, ::Nothing) = Ones{Float64}(sum(esample))

#  remove observations with missing or negative weights
isnaorneg(a::AbstractVector) = BitArray(!ismissing(x) & (x > 0) for x in a)


text:
## Copyright (c) 2013 Miguel Bazdresch
##
## This file is distributed under the 2-clause BSD License.

# write commands to gnuplot's pipe
function gnuplot_send(s::String)
    fid = gnuplot_state.fid
    err = ccall(:fputs, Int, (Ptr{Uint8},Ptr{Int}), string(s,"\n"), fid)
    # fputs returns a positive number if everything worked all right
    if err < 0
        println("Something went wrong writing to the gnuplot pipe.")
        return
    end
    err = ccall(:fflush, Int, (Ptr{Int},), fid)
    ## fflush returns 0 if everything worked all right
    if err != 0
        println("Something went wrong writing to the gnuplot pipe.")
        return
    end
end


text:
abstract type GlyphData end

@bitmask_flag SimpleGlyphFlag::UInt8 begin
    ON_CURVE_POINT_BIT =                       0x01
    X_SHORT_VECTOR_BIT =                       0x02
    Y_SHORT_VECTOR_BIT =                       0x04
    REPEAT_FLAG_BIT =                          0x08
    X_IS_SAME_OR_POSITIVE_X_SHORT_VECTOR_BIT = 0x10
    Y_IS_SAME_OR_POSITIVE_Y_SHORT_VECTOR_BIT = 0x20
    OVERLAP_SIMPLE_BIT =                       0x40
    SIMPLE_GLYPH_RESERVED_BIT =                0x80
end

struct GlyphPoint
    coords::Point{2,Int}
    on_curve::Bool
end

"""
Description of a glyph as a series of quadratic bezier patches.

Bezier patches are implicitly defined using a list of `GlyphPoint`s, where two consecutive off-curve points implicitly define an on-curve point halfway.
"""
struct GlyphSimple <: GlyphData
    contour_indices::Vector{Int}
    points::Vector{GlyphPoint}
end

Base.show(io::IO, gdata::GlyphSimple) = print(io, "GlyphSimple(", length(gdata.points), " points, ", length(gdata.contour_indices), " indices)")

struct GlyphHeader
    ncontours::Int16
    xmin::Int16
    ymin::Int16
    xmax::Int16
    ymax::Int16
end

struct Glyph{D<:GlyphData}
    header::GlyphHeader
    data::D
end

"""
Extract all quadratic bezier control points explicitly from a glyph.

In OpenType, control points may be stored using an implicit format, where two consecutive points may be marked as off-curve although the curve is quadratic: in this case, it means that there is an implicit control point halfway.
"""
function uncompress(glyph::Glyph)
    data = glyph.data

    contour_indices = [0; data.contour_indices]
    ranges = map(zip(contour_indices[begin:end-1], contour_indices[begin+1:end])) do (i, j)
        (i+1):j
    end

    curves = Vector{Point{2,Float64}}[]
    for data_points in map(Base.Fix1(getindex, data.points), ranges)
        points = Point{2,Float64}[]

        # make sure data points define a closed contour
        while !(first(data_points).on_curve)
            push!(data_points, popfirst!(data_points))
        end
        if last(data_points) ≠ first(data_points)
            # terminate with a linear segment
            push!(data_points, first(data_points))
        end

        # gather contour points including implicit ones
        on_curve = false
        for point in data_points
            coords = point.coords
            if !on_curve && !point.on_curve || on_curve && point.on_curve
                # there is an implicit on-curve point halfway
                push!(points, (coords + points[end]) / 2)
            end
            push!(points, coords)
            on_curve = point.on_curve
        end

        @assert isodd(length(points)) "Expected an odd number of curve points."
        @assert first(points) == last(points) "Contour is not closed."
        push!(curves, points)
        on_curve = true
    end

    curves
end

function normalize(curves, glyph::Glyph)
    min = Point(glyph.header.xmin, glyph.header.ymin)
    max = Point(glyph.header.xmax, glyph.header.ymax)
    sc = inv(Scaling(max - min))
    transf = sc ∘ Translation(-min)
    map(curves) do points
        @assert all(min[1] ≤ minimum(getindex.(points, 1)))
        @assert all(min[2] ≤ minimum(getindex.(points, 2)))
        @assert all(max[1] ≥ maximum(getindex.(points, 1)))
        @assert all(max[2] ≥ maximum(getindex.(points, 2)))
        res = transf.(points)
        @assert all(p -> all(0 .≤ p .≤ 1), res)
        res
    end
end

function curves(glyph::Glyph)
    patch = Patch(BezierCurve(), 3)
    curves = normalize(uncompress(glyph), glyph)
    [map(Base.Fix2(split, patch), curves)...;]
end

function read_glyphs(io::IO, head::FontHeader, maxp::MaximumProfile, nav::TableNavigationMap, loca::IndexToLocation)
    glyphs = map(glyph_ranges(loca)) do range
        range.stop == range.start && return nothing
        read_table(io, nav["glyf"]; offset = range.start, length = range.stop - range.start) do io
            header = GlyphHeader(
                (read(io, T) for T in fieldtypes(GlyphHeader))...
            )
            data = if header.ncontours ≠ -1
                end_contour_points = [read(io, UInt16) for _ in 1:header.ncontours]

                # convert to 1-based indexing
                end_contour_points .+= 1

                instlength = read(io, UInt16)
                insts = [read(io, UInt8) for _ in 1:instlength]
                end_idx = end_contour_points[end]
                flags = SimpleGlyphFlag[]
                while length(flags) < end_idx
                    flag = SimpleGlyphFlag(read(io, UInt8))
                    push!(flags, flag)
                    if REPEAT_FLAG_BIT in flag
                        repeat_count = read(io, UInt8)
                        append!(flags, (flag for _ in 1:repeat_count))
                    end
                end

                xs = Int[]
                foreach(flags) do flag
                    x = if X_IS_SAME_OR_POSITIVE_X_SHORT_VECTOR_BIT in flag && X_SHORT_VECTOR_BIT ∉ flag
                        val = isempty(xs) ? 0 : last(xs)
                        push!(xs, val)
                        return
                    elseif X_SHORT_VECTOR_BIT in flag
                        val = Int(read(io, UInt8))
                        X_IS_SAME_OR_POSITIVE_X_SHORT_VECTOR_BIT in flag ? val : -val
                    elseif X_IS_SAME_OR_POSITIVE_X_SHORT_VECTOR_BIT ∉ flag && X_SHORT_VECTOR_BIT ∉ flag
                        read(io, Int16)
                    end
                    push!(xs, x + (isempty(xs) ? 0 : last(xs)))
                end

                ys = Int[]
                foreach(flags) do flag
                    y = if Y_IS_SAME_OR_POSITIVE_Y_SHORT_VECTOR_BIT in flag && Y_SHORT_VECTOR_BIT ∉ flag
                        val = isempty(ys) ? 0 : last(ys)
                        push!(ys, val)
                        return
                    elseif Y_SHORT_VECTOR_BIT in flag
                        val = Int(read(io, UInt8))
                        Y_IS_SAME_OR_POSITIVE_Y_SHORT_VECTOR_BIT in flag ? val : -val
                    elseif Y_IS_SAME_OR_POSITIVE_Y_SHORT_VECTOR_BIT ∉ flag && Y_SHORT_VECTOR_BIT ∉ flag
                        read(io, Int16)
                    end
                    push!(ys, y + (isempty(ys) ? 0 : last(ys)))
                end
                GlyphSimple(
                    end_contour_points,
                    GlyphPoint.(collect(zip(xs, ys)), map(Base.Fix1(in, ON_CURVE_POINT_BIT), flags)),
                )
            end
            if position(io) % 4 ≠ 0
                # rest should just be padding zeros
                @assert all(iszero, read(io, UInt8) for _ in 1:4 - position(io) % 4)
            end
            Glyph(header, data)
        end
    end
end


text:
# This file includes RODESystems

import DifferentialEquations: RandomEM, RODEProblem
import UUIDs: uuid4

"""
    @def_rode_system ex

where `ex` is the expression to define to define a new AbstractRODESystem component type. The usage is as follows:
```julia
@def_rode_system mutable struct MyRODESystem{T1,T2,T3,...,TN,OP,RH,RO,ST,IP,OP} <: AbstractRODESystem
    param1::T1 = param1_default                 # optional field
    param2::T2 = param2_default                 # optional field
    param3::T3 = param3_default                 # optional field
        ⋮
    paramN::TN = paramN_default                 # optional field
    righthandside::RH = righthandside_function  # mandatory field
    readout::RO = readout_function              # mandatory field
    state::ST = state_default                   # mandatory field
    input::IP = input_default                   # mandatory field
    output::OP = output_default                 # mandatory field
end
```
Here, `MyRODESystem` has `N` parameters. `MyRODESystem` is represented by the `righthandside` and `readout` function. `state`, `input` and `output` is the initial state, input port and output port of `MyRODESystem`.

!!! warning
    `righthandside` must have the signature
    ```julia
    function righthandside((dx, x, u, t, W, args...; kwargs...)
        dx .= .... # update dx
    end
    ```
    and `readout` must have the signature
    ```julia
    function readout(x, u, t)
        y = ...
        return y
    end
    ```

!!! warning
    New RODE system must be a subtype of `AbstractRODESystem` to function properly.

# Example
```jldoctest
julia> @def_rode_system mutable struct MySystem{RH, RO, IP, OP} <: AbstractRODESystem
           A::Matrix{Float64} = [2. 0.; 0 -2]
           righthandside::RH = (dx, x, u, t, W) -> (dx .= A * x * W)
           readout::RO = (x, u, t) -> x
           state::Vector{Float64} = rand(2)
           input::IP = nothing
           output::OP = Outport(2)
       end

julia> ds = MySystem();
```
"""
macro def_rode_system(ex)
    checksyntax(ex, :AbstractRODESystem)
    fields = quote
        trigger::$(TRIGGER_TYPE_SYMBOL) = Inpin()
        handshake::$(HANDSHAKE_TYPE_SYMBOL) = Outpin{Bool}()
        callbacks::$(CALLBACKS_TYPE_SYMBOL) = nothing
        name::Symbol = Symbol()
        id::$(ID_TYPE_SYMBOL) = Causal.uuid4()
        t::Float64 = 0.
        modelargs::$(MODEL_ARGS_TYPE_SYMBOL) = ()
        modelkwargs::$(MODEL_KWARGS_TYPE_SYMBOL) = NamedTuple()
        solverargs::$(SOLVER_ARGS_TYPE_SYMBOL) = ()
        solverkwargs::$(SOLVER_KWARGS_TYPE_SYMBOL) = (dt=0.01, )
        alg::$(ALG_TYPE_SYMBOL) = Causal.RandomEM()
        integrator::$(INTEGRATOR_TYPE_SYMBOL) = Causal.construct_integrator(Causal.RODEProblem, input, righthandside, state, t, modelargs,
            solverargs; alg=alg, modelkwargs=modelkwargs, solverkwargs=solverkwargs, numtaps=3)
    end, [TRIGGER_TYPE_SYMBOL, HANDSHAKE_TYPE_SYMBOL, CALLBACKS_TYPE_SYMBOL, ID_TYPE_SYMBOL, MODEL_ARGS_TYPE_SYMBOL, MODEL_KWARGS_TYPE_SYMBOL, SOLVER_ARGS_TYPE_SYMBOL, SOLVER_KWARGS_TYPE_SYMBOL, ALG_TYPE_SYMBOL, INTEGRATOR_TYPE_SYMBOL]
    _append_common_fields!(ex, fields...)
    deftype(ex)
end

##### Define RODE sytem library

"""
    RODESystem(; righthandside, readout, state, input, output)

Constructs a generic RODE system
"""
@def_rode_system mutable struct RODESystem{RH, RO, ST, IP, OP} <: AbstractRODESystem
    righthandside::RH
    readout::RO
    state::ST
    input::IP
    output::OP
end

@doc raw"""
    MultiplicativeNoiseLinearSystem()

Constructs a `MultiplicativeNoiseLinearSystem` with the dynamics
```math
\begin{array}{l}
    \dot{x} = A x W
\end{array}
where `W` is the noise process.
```
"""
@def_rode_system mutable struct MultiplicativeNoiseLinearSystem{RH, RO, IP, OP} <: AbstractRODESystem
    A::Matrix{Float64} = [2. 0.; 0 -2]
    righthandside::RH = (dx, x, u, t, W) -> (dx .= A * x * W)
    readout::RO = (x, u, t) -> x
    state::Vector{Float64} = rand(2)
    input::IP = nothing
    output::OP = Outport(2)
end

##### Pretty printing
show(io::IO, ds::RODESystem) = print(io,
    "RODESystem(righthandside:$(ds.righthandside), readout:$(ds.readout), state:$(ds.state), t:$(ds.t), input:$(ds.input), output:$(ds.output))")
show(io::IO, ds::MultiplicativeNoiseLinearSystem) = print(io,
    "MultiplicativeNoiseLinearSystem(A:$(ds.A), state:$(ds.state), t:$(ds.t), input:$(ds.input), output:$(ds.output))")


text:
#!/usr/bin/env julia

using JuMP, Ipopt
using YAML
using JLD, NPZ

println("loading utility functions....")
include("modules/setup.jl")

println("loading configurations ....")
include("obstacle-configuration.jl")
include("vehicle-configuration.jl")
include("modules/mpcPathFollowing.jl")

home            = ENV["HOME"]
param_file 		= string(home,"/barc/workspace/src/barc/launch/control/parking.yaml")
save_path 		= string(home,"/barc/workspace/src/barc/src/control/parking/data")

data = YAML.load(open(param_file))

s0 	= data["initial_state"]
sF 	= data["final_state"]
L 	= data["vehicle_length"]

println("calling hybrid a* to get initial guess for warm start ...")
rx, ry, ryaw = hybrid_a_star.calc_hybrid_astar_path(s0[1], s0[2], s0[3], sF[1], sF[2], sF[3], ox, oy, hybrid_a_star.XY_GRID_RESOLUTION, hybrid_a_star.YAW_GRID_RESOLUTION, hybrid_a_star.OB_MAP_RESOLUTION)

println("smoothing hybrid a* solution ...")
xWS, uWS, N = getSmoothProfile(rx,ry,ryaw,Ts,sampleN,motionStep)

println("computing H-representation for obstacle polytopes ....")
AOb, bOb = obstHrep(nOb, vOb, lOb)

println("solving finite time optimal control problem ...")
sp10, up10, scaleTime10, exitflag10, time10, lp10, np10 = ParkingSignedDist(s0,sF,N,Ts,L,ego,XYbounds,nOb,vObMPC,AOb,bOb,fixTime,xWS,uWS)

println("interpolating solution ....")
tDomOpt = unshift!(cumsum(scaleTime10*Ts),0)
pop!(tDomOpt)
tDomSim 	= 0:0.1:round(tDomOpt[end-1],1)
fsp10 		= intp.interp1d(tDomOpt, sp10)
fup10 		= intp.interp1d(tDomOpt[1:end-1], up10)
sp10Ref 	= fsp10(tDomSim)
up10Ref 	= fup10(tDomSim)
dp10Ref 	=  cumsum( sqrt.( diff(sp10Ref[1,:]).^2 .+ diff(sp10Ref[2,:]).^2 ) )
unshift!(dp10Ref,0)

println("saving solution ...")
save("$save_path/parkPath.jld","sp10",sp10,"up10",up10,"scaleTime10",scaleTime10,"sp10Ref",sp10Ref,"up10Ref",up10Ref,"dp10Ref",dp10Ref)
npzwrite("$save_path/npz/sp10.npz", sp10)
npzwrite("$save_path/npz/up10.npz", up10)
npzwrite("$save_path/npz/scaleTime10.npz", scaleTime10)
npzwrite("$save_path/npz/sp10Ref.npz", sp10Ref)
npzwrite("$save_path/npz/up10Ref.npz", up10Ref)
npzwrite("$save_path/npz/dp10Ref.npz", dp10Ref)

println("plotting solution path")
for j = 1 : nObPlot
	for k = 1 : vObPlot[j]
		plot([lObPlot[j][k][1],lObPlot[j][k+1][1]] , [lObPlot[j][k][2],lObPlot[j][k+1][2]] ,"k")
	end
end
plot(sp10[1,:],sp10[2,:], "-b",  label="H-OBCA")
plot(s0[1],s0[2],"ob")
plot(sF[1],sF[2],"ob")
legend()
axis("equal")

max_stars_repo_name: treyfortmuller/barc

max_issues_repo_name: treyfortmuller/barc