text:
const USE_MIN_RANGE = true

struct AffineEAGO{N}
    c::Float64               # mid-point
    γ::SVector{N,Float64}    # affine terms
    Δ::Float64               # error term
end

function AffineEAGO(x::AffineEAGO{N}, p::Float64, q::Float64, δ::Float64) where N
    c = p * x.c + q
    γ = p .* x.γ
    Δ = p * x.Δ + δ
    AffineEAGO{N}(c, γ, δ)
end
mid(x::Interval{Float64}) = 0.5*(x.lo + x.hi)
AffineEAGO{N}(x::Float64, X::Interval{Float64}, i::Int) where N = AffineEAGO{N}(mid(X), radius(X)*seed_gradient(i, Val(N)), 0.0)

const UNIT_INTERVAL = Interval{Float64}(-1,1)
Interval(x::AffineEAGO{N}) where N = x.c + x.Δ*UNIT_INTERVAL + sum(y -> abs(y)*UNIT_INTERVAL, x.γ)
function bounds(x::AffineEAGO{N}) where N
    z = Interval(x)
    z.lo, z.hi
end

zero(::Type{AffineEAGO{N}}) where N = AffineEAGO{N}(0.0, zero(SVector{N,Float64}), 0.0)
zero(x::AffineEAGO{N}) where N = AffineEAGO{N}(0.0, zero(SVector{N,Float64}), 0.0)

one(::Type{AffineEAGO{N}}) where N = AffineEAGO{N}(1.0, zero(SVector{N,Float64}), 0.0)
one(x::AffineEAGO{N}) where N = AffineEAGO{N}(1.0, zero(SVector{N,Float64}), 0.0)

+(x::AffineEAGO{N}, y::AffineEAGO{N}) where N = AffineEAGO{N}(x.c + y.c, x.γ .+ y.γ, x.Δ + y.Δ)
+(x::AffineEAGO{N}, α::Float64) where N = AffineEAGO{N}(α+x.c, x.γ, x.Δ)
+(α::Float64, x::AffineEAGO{N}) where N = x + α

function *(x::AffineEAGO{N}, y::AffineEAGO{N}) where N
    x0 = x.c
    y0 = y.c
    γ = SVector{N,Float64}(ntuple(i -> x0*y.γ[i] + y0*x.γ[i], Val(N)))
    Δ = abs(x0)*y.Δ + abs(y0)*x.Δ
    sx = abs(x.Δ)
    sy = abs(y.Δ)
    for i = 1:N
        sx += abs(x.γ[i])
        sy += abs(y.γ[i])
    end
    Δ += sx + sy
    return AffineEAGO{N}(x0*y0, γ, Δ)
end
*(x::AffineEAGO{N}, α::Float64) where N = AffineEAGO{N}(α*x.c, α.*x.γ, abs(α)*x.Δ)
*(α::Float64, x::AffineEAGO{N}) where N = x*α

function -(x::AffineEAGO{N}, y::AffineEAGO{N}) where N
    AffineEAGO{N}(x.c - y.c, x.γ .- y.γ, x.Δ + y.Δ)
end
-(x::AffineEAGO{N}) where N = AffineEAGO{N}(-x.c, .-(x.γ), x.Δ)
-(x::AffineEAGO{N}, α::Float64) where N = AffineEAGO{N}(x.c - α, x.γ, x.Δ)
-(α::Float64, x::AffineEAGO{N}) where N = α + (-x)

/(x::AffineEAGO{N}, α::Float64) where N = AffineEAGO{N}(x.c/α, x.γ/α, x.Δ/abs(α))
/(α::Float64, x::AffineEAGO{N}) where N = α*inv(x)

function exp(x::AffineEAGO{N}) where N
    a, b = bounds(x)
    fa = exp(a)
    fb = exp(b)
    if USE_MIN_RANGE
        p = exp(a)
        q = 0.5*(fa + fb - p*(a + b))
        Δ = abs(0.5*(fb - fa - p*(b - a)))
        return AffineEAGO(x, p, q, Δ)
    end
    p = (fb - fa)/(b - a)
    ξ = log(p)
    fξ = p
    q = 0.5*(fa + fξ - p*(a + ξ))
    Δ = abs(0.5*(fξ - fa - p*(ξ - a)))
    return AffineEAGO(x, p, q, Δ)
end

function exp10(x::AffineEAGO{N}) where N
    a, b = bounds(x)
    fa = exp10(a)
    fb = exp10(b)
    if USE_MIN_RANGE
        p = log10(a)
        q = 0.5*(fa + fb - p*(a + b))
        Δ = abs(0.5*(fb - fa - p*(b - a)))
        return AffineEAGO(x, p, q, Δ)
    end
    p = (fb - fa)/(b - a)
    ξ = log10(p)
    fξ = p
    q = 0.5*(fa + fξ - p*(a + ξ))
    Δ = abs(0.5*(fξ - fa - p*(ξ - a)))
    return AffineEAGO(x, p, q, Δ)
end

function log(x::AffineEAGO{N}) where N
    a, b = bounds(x)
    fa = log(a)
    fb = log(b)
    if USE_MIN_RANGE
        p = 1/b
        q = 0.5*(fa + fb - p*(a + b))
        Δ = abs(0.5*(fb - fa - p*(b - a)))
        return AffineEAGO(x, p, q, Δ)
    end
    p = (fb - fa)/(b - a)
    ξ = 1/p
    fξ = log(ξ)
    q = 0.5*(fa + fξ - p*(a + ξ))
    Δ = abs(0.5*(fξ - fa - p*(ξ - a)))
    return AffineEAGO(x, p, q, Δ)
end

function log10(x::AffineEAGO{N}) where N
    a, b = bounds(x)
    fa = log10(a)
    fb = log10(b)
    if USE_MIN_RANGE
        p = 1/(b*log(10))
        q = 0.5*(fa + fb - p*(a + b))
        Δ = abs(0.5*(fb - fa - p*(b - a)))
        return AffineEAGO(x, p, q, Δ)
    end
    p = (fb - fa)/(b - a)
    ξ = 1/p
    fξ = log10(ξ)
    q = 0.5*(fa + fξ - p*(a + ξ))
    Δ = abs(0.5*(fξ - fa - p*(ξ - a)))
    return AffineEAGO(x, p, q, Δ)
end

function pow_1d(x::AffineEAGO{N}, n::Number, p) where N
    a, b = bounds(x)
    fa = a^n
    fb = b^n
    if USE_MIN_RANGE
        q = 0.5*(fa + fb - p*(a + b))
        Δ = abs(0.5*(fb - fa - p*(b - a)))
        return AffineEAGO(x, p, q, Δ)
    end
    p = (fb - fa)/(b - a)
    ξ = (p/n)^(1/(n - 1))
    fξ = ξ^n
    q = 0.5*(fa + fξ - p*(a + ξ))
    Δ = abs(0.5*(fξ - fa - p*(ξ - a)))
    return AffineEAGO(x, p, q, Δ)
end

function pow_even(x::AffineEAGO{N}, n::Int) where N
    a, b = bounds(x)
    fa = a^n
    fb = b^n
    if USE_MIN_RANGE
        m = min(0.0, fa, fb)
        M = max(0.0, fa, fb)
        p = 0.0
        q = 0.5*(m + M)
        Δ = 0.5*(M - m)
        return AffineEAGO(x, p, q, Δ)
    end
    p = (fb - fa)/(b - a)
    ξ = (p/n)^(1/(n - 1))
    fξ = ξ^n
    q = 0.5*(fa + fξ - p*(a + ξ))
    Δ = abs(0.5*(fξ - fa - p*(ξ - a)))
    return AffineEAGO(x, p, q, Δ)
end

function pow_odd(x::AffineEAGO{N}, n::Int) where N
    #println("ran power odd")
    # TODO: DOES THIS HANDLE a <= 0.0 <= b?
    a, b = bounds(x)
    fa = a^n
    fb = b^n
    p = (fb - fa)/(b - a)
    q = 0.0
    ξ = (p/n)^(1/(n-1))
    fξ = ξ^n
    Δ = abs(fξ - p*ξ)
    return AffineEAGO(x, p, q, Δ)
    #y = Base.power_by_squaring(x,n)
   # return y
end

function ^(x::AffineEAGO{N}, n::Int) where N

    iszero(n) && zero(x)
    isone(n)  && one(x)

    xL, xU = bounds(x)
    if (xL > 0.0) || (xU < 0.0)
        return pow_1d(x, n, (n >= 0) ? n*xL^(n-1) : n*xU^(n-1))
    elseif iseven(n)
        return pow_even(x, n)
    end
    return pow_odd(x, n)
end

function ^(x::AffineEAGO{N}, n::Float64) where N
    if isinteger(n)
        return x^Int(n)
    end
    xL, xU = bounds(x)
    (xL < 0.0) && error("Invalid domain...")
    if (n > 1.0) || (n < 0.0)
        return pow_1d(x, n, n*xU^(n-1))
    end
    return pow_1d(x, n, n*xL^(n-1))
end

^(x::AffineEAGO{N}, y::AffineEAGO{N}) where N = exp(y*log(x))
function ^(x::AffineEAGO{N}, n) where N
    if iszero(n)
        return zero(x)
    elseif isone(n)
        return one(x)
    end
    return x^n
end

function inv(x::AffineEAGO{N}) where N
    a, b = bounds(x)
    (a < 0.0 < b) && error("Invalid domain...")
    if b < 0.0
        return -inv(-x)
    end
    if USE_MIN_RANGE
        p = -1/b^2
        q = -(p*(a + b)^2)/(2*a)
        Δ = -(p*(a - b)^2)/(2*a)
        return AffineEAGO(x, p, q, Δ)
    end
    p = -1/(a*b)
    q = -0.5*p*(sqrt(a) + sqrt(b))^2
    Δ = -0.5*p*(sqrt(a) - sqrt(b))^2
    return AffineEAGO(x, p, q, Δ)
end

struct MCAffPnt{N,T}
    v::MC{N,T}
    box::AffineEAGO{N}
end
MC(x::MCAffPnt{N,T}) where {N,T<:RelaxTag} = x.v
MC(x::MC{N,T}) where {N, T<:RelaxTag} = x

relax_info(s::RelaxAA, n::Int, t::T) where {N,T} = MCAffPnt{n,T}

zero(::Type{MCAffPnt{N,T}}) where {N,T} = MCAffPnt{N,T}(zero(MC{N,T}), zero(AffineEAGO{N}))
zero(x::MCAffPnt{N,T}) where {N,T} = MCAffPnt{N,T}(zero(x.v), zero(x.box))

one(::Type{MCAffPnt{N,T}}) where {N,T} = MCAffPnt{N,T}(one(MC{N,T}), one(AffineEAGO{N}))
one(x::MCAffPnt{N,T}) where {N,T} = MCAffPnt{N,T}(one(MC{N,T}), one(AffineEAGO{N}))

+(x::MCAffPnt{N,T}, y::MCAffPnt{N,T}) where {N,T} = MCAffPnt{N,T}(x.v + y.v, x.box + y.box)
+(x::MCAffPnt{N,T}, α::Number) where {N,T} = MCAffPnt{N,T}(x.v + α, x.box + α)
+(α::Number, x::MCAffPnt{N,T}) where {N,T} = x + α

*(x::MCAffPnt{N,T}, y::MCAffPnt{N,T}) where {N,T} = MCAffPnt{N,T}(x.v*y.v, x.box*y.box)
*(x::MCAffPnt{N,T}, α::Number) where {N,T} = MCAffPnt{N,T}(x.v * α, x.box * α)
*(α::Number, x::MCAffPnt{N,T}) where {N,T} = x*α

-(x::MCAffPnt{N,T}, y::MCAffPnt{N,T}) where {N,T} = MCAffPnt{N,T}(x.v-y.v, x.box-y.box)
-(x::MCAffPnt{N,T}) where {N,T} = MCAffPnt{N,T}(-x.v, -x.box)
-(x::MCAffPnt{N,T}, α::Number) where {N,T} = x + (-α)
-(α::Number, x::MCAffPnt{N,T}) where {N,T} = α + (-x)

/(x::MCAffPnt{N,T}, y::MCAffPnt{N,T}) where {N,T} = x*inv(y)
/(x::MCAffPnt{N,T}, α::T) where {N,T} = MCAffPnt{N,T}(x.v/α, x.box/α)
/(α::T, x::MCAffPnt{N,T}) where {N,T} = MCAffPnt{N,T}(α*inv(x.v), α*inv(x.box))

^(x::MCAffPnt{N,T}, n::Integer) where {N,T} = MCAffPnt{N,T}(x.v^n, x.box^n)
^(x::MCAffPnt{N,T}, n::Number) where {N,T} = MCAffPnt{N,T}(x.v^n, x.box^n)
^(x::MCAffPnt{N,T}, n::MCAffPnt{N,T}) where {N,T} = MCAffPnt{N,T}(x.v^n.v, x.box^n.box)

for op in (:inv, :log, :log10, :exp, :exp10)
    @eval ($op)(x::MCAffPnt{N,T}) where {N,T} = MCAffPnt{N,T}(($op)(x.v), ($op)(x.box))
end

Interval(x::MCAffPnt{N,T}) where {N,T<:RelaxTag} = Interval(x.v) ∩ Interval(x.box)

function cut(x::MCAffPnt{N,T}, z::MCAffPnt{N,T}, v::VariableValues, ϵ::Float64, s::Vector{Int}, cflag::Bool, pflag::Bool) where {N,T<:RelaxTag}
    (pflag & cflag)  && (return set_value_post(x ∩ Interval(z), v, s, ϵ))
    (pflag & !cflag) && (return set_value_post(x, v, s, ϵ))
    (pflag & cflag)  && (return x ∩ Interval(z))
    return x
end

function cut(x::MC{N,T}, z::MCAffPnt{N,T}, v::VariableValues, ϵ::Float64, s::Vector{Int}, cflag::Bool, pflag::Bool) where {N,T<:RelaxTag}
    (pflag & cflag)  && (return set_value_post(x ∩ Interval(z), v, s, ϵ))
    (pflag & !cflag) && (return set_value_post(x, v, s, ϵ))
    (pflag & cflag)  && (return x ∩ Interval(z))
    return x
end

function varset(::Type{MCAffPnt{N,T}}, i, x_cv, x_cc, l, u) where {V,N,T<:RelaxTag}
    v = seed_gradient(i, Val(N))
    v_Intv = Interval{Float64}(l, u)
    v_mc = MC{N,T}(x_cv, x_cc, v_Intv, v, v, false)
    v_aff = AffineEAGO{N}(x_cv, v_Intv, i)
    return MCAffPnt{N,T}(v_mc, v_aff)
end

function fprop!(t::RelaxAAInfo, vt::Variable, g::DAT, b::RelaxCache{MCAffPnt{N,T},N,T}, k) where {N,T<:RelaxTag}
    i = first_index(g, k)
    x = val(b, i)
    l = lbd(b, i)
    u = ubd(b, i)
    if l == u
        b[k] = x
        b._is_num[k] = true
    else
        z = varset(MCAffPnt{N,T}, rev_sparsity(g, i, k), x, x, l, u)
        if !first_eval(t, b)
            z = z ∩ interval(b, k)
        end
        b._info[k] = z
        b._is_num[k] = false
    end
    nothing
end

function fprop!(t::RelaxAAInfo, ex::Subexpression, g::DAT, b::RelaxCache{MCAffPnt{N,T},N,T}, k) where {N,T<:RelaxTag}
    x =  first_index(g, k)
    if subexpression_is_num(b, x)
        b[k] = subexpression_num(b, x)
        b._is_num[k] = true
    else
        b._info[k] = subexpression_info(b, x)
        b._is_num[k] = false
    end
end

for (F, f) in ((PLUS, :+), (MIN, :min), (MAX, :max), (DIV, :/), (ARH, :arh))
    @eval function fprop_2!(t::RelaxAAInfo, v::Val{$F}, g::DAT, b::RelaxCache{MCAffPnt{N,T},N,T}, k) where {N,T<:RelaxTag}
        x = child(g, 1, k)
        y = child(g, 2, k)
        if !xy_num(b, x, y)
            if xyset(b, x, y)
                z = ($f)(info(b, x), info(b, y))
            elseif xset_ynum(b, x, y)
                z = ($f)(info(b, x), num(b, y))
            else
                z = ($f)(num(b, x), info(b, y))
            end
            b._info[k] = cut(z, info(b,k), b.ic.v, b.ϵ_sg, sparsity(g, k), b.cut, false)
            b._is_num[k] = false
        else
            b[k] = ($f)(num(b, x), num(b, y))
            b._is_num[k] = true
        end
    end
end

function fprop!(t::RelaxAAInfo, v::Val{MINUS}, g::DAT, b::RelaxCache{MCAffPnt{N,T},N,T}, k) where {N,T<:RelaxTag}
    x = child(g, 1, k)
    if is_binary(g, k)
        y = child(g, 2, k)
        if !xy_num(b, x, y)
            if xyset(b, x, y)
                z = info(b, x) - info(b, y)
            elseif xset_ynum(b, x, y)
                z = info(b, x) - num(b, y)
            else
                z = num(b, x) - info(b, y)
            end
            b._info[k] = cut(z, info(b,k), b.ic.v, b.ϵ_sg, sparsity(g, k), b.cut, false)
            b._is_num[k] = false
        else
            b[k] = num(b, x) - num(b, y)
            b._is_num[k] = true
        end
    else
        if is_num(b, x)
            b[k] = -num(b, x)
            b._is_num[k] = true
        else
            b._info[k] = cut(-info(b, x), info(b, k), b.ic.v, b.ϵ_sg, sparsity(g, k), b.cut, false)
            b._is_num[k] = false
        end
    end
end

function fprop_n!(t::RelaxAAInfo, v::Val{PLUS}, g::DAT, b::RelaxCache{MCAffPnt{N,T},N,T}, k::Int) where {N,T<:RelaxTag}
    z = zero(MCAffPnt{N,T})
    znum = 0.0
    numval = true
    for i in children(g, k)
        if is_num(b, i)
            znum += num(b, i)
        else
            numval = false
            z += info(b, i)
        end
    end
    if numval
        b[k] = znum
        b._is_num[k] = true
    else
        z += znum
        b._info[k] = cut(z, info(b, k), b.ic.v, b.ϵ_sg, sparsity(g, k), b.cut, false)
        b._is_num[k] = false
    end
end

function fprop_n!(t::RelaxAAInfo, v::Val{MIN}, g::DAT, b::RelaxCache{MCAffPnt{N,T},N,T}, k::Int) where {N,T<:RelaxTag}
    z = Inf*one(MCAffPnt{N,T})
    znum = Inf
    numval = true
    for i in children(g, k)
        if is_num(b, i)
            znum = min(znum, num(b, i))
        else
            numval = false
            z = min(z, info(b, i))
        end
    end
    if numval
        b[k] = znum
        b._is_num[k] = true
    else
        z = min(z, znum)
        b._info[k] = cut(z, info(b, k), b.ic.v, b.ϵ_sg, sparsity(g, k), b.cut, false)
        b._is_num[k] = false
    end
end

function fprop_n!(t::RelaxAAInfo, v::Val{MAX}, g::DAT, b::RelaxCache{MCAffPnt{N,T},N,T}, k::Int) where {N,T<:RelaxTag}
    z = -Inf*one(MCAffPnt{N,T})
    znum = -Inf
    numval = true
    for i in children(g, k)
        if is_num(b, i)
            znum = max(znum, num(b, i))
        else
            numval = false
            z = max(z, info(b, i))
        end
    end
    if numval
        b[k] = znum
        b._is_num[k] = true
    else
        z = max(z, znum)
        b._info[k] = cut(z, info(b, k), b.ic.v, b.ϵ_sg, sparsity(g, k), b.cut, false)
        b._is_num[k] = false
    end
end

function fprop_2!(t::RelaxAAInfo, v::Val{MULT}, g::DAT, b::RelaxCache{MCAffPnt{N,T},N,T}, k::Int) where {N,T<:RelaxTag}

    x = child(g, 1, k)
    y = child(g, 2, k)

    if !xy_num(b, x, y)
        if xyset(b, x, y)
            xr = info(b, x)
            yr = info(b, y)
            xv = xr.v
            yv = yr.v
            if b.use_apriori_mul
                dp = b.dp
                dP = b.dP
                p_rel = b.p_rel
                p_diam = b.p_diam
                s = sparsity(g, 1)
                u1max, u2max, v1nmax, v2nmax = estimator_extrema(xr, yr, s, dP)
                zv = xv*yv
                wIntv = zv.Intv
                if (u1max < xv.Intv.hi) || (u2max < yv.Intv.hi)
                    u1cv, u2cv, u1cvg, u2cvg = estimator_under(xv, yv, xr, yr, s, dp, dP, p_rel, p_diam)
                    za_l = McCormick.mult_apriori_kernel(xv, yv, wIntv, u1cv, u2cv, u1max, u2max, u1cvg, u2cvg)
                    zv = zv ∩ za_l
                end
                if (v1nmax > -xv.Intv.lo) || (v2nmax > -yv.Intv.lo)
                    v1ccn, v2ccn, v1ccgn, v2ccgn = estimator_over(xv, yv, xr, yr, s, dp, dP, p_rel, p_diam)
                    za_u = McCormick.mult_apriori_kernel(-xv, -yv, wIntv, v1ccn, v2ccn, v1nmax, v2nmax, v1ccgn, v2ccgn)
                    zv = zv ∩ za_u
                end
                z = MCAffPnt{N,T}(zv, xr.box*yr.box)
            else
                z = xv*yv
            end
        elseif xset_ynum(b, x, y)
            z = info(b, x)*num(b, y)
        else
            z = num(b, x)*info(b, y)
        end
        b._info[k] = cut(z, info(b,k), b.ic.v, b.ϵ_sg, sparsity(g, k), b.cut, false)
        b._is_num[k] = true
    else
        b[k] = num(b, x)*num(b, y)
        b._is_num[k] = false
    end
end

function fprop_n!(t::RelaxAAInfo, ::Val{MULT}, g::DAT, b::RelaxCache{MCAffPnt{N,T},N,T}, k::Int) where {N,T<:RelaxTag}

    z = one(MCAffPnt{N,T})
    znum = one(Float64)
    numval = true
    if b.use_apriori_mul
        zr = one(MCAffPnt{N,T})
        dp = b.dp
        dP = b.dP
        s = sparsity(g, 1)
        for (q,i) in enumerate(children(g, k))
            if is_num(b, i)
                znum = znum*num(b, i)
            else
                numval = false
                xi = info(b, i)
                x = xi.v
                xr = info(b, i)
                u1max, u2max, v1nmax, v2nmax = estimator_extrema(zr, xr, s, dP)
                zv = z*x
                wIntv = zv.Intv
                if (u1max < z.Intv.hi) || (u2max < x.Intv.hi)
                    u1cv, u2cv, u1cvg, u2cvg = estimator_under(zr, xr, s, dp, dP)
                    za_l = McCormick.mult_apriori_kernel(z, x, wIntv, u1cv, u2cv, u1max, u2max, u1cvg, u2cvg)
                    zv = zv ∩ za_l
                end
                if (v1nmax > -z.Intv.lo) || (v2nmax > -x.Intv.lo)
                    v1ccn, v2ccn, v1ccgn, v2ccgn = estimator_under(zr, xr, s, dp, dP)
                    za_u = McCormick.mult_apriori_kernel(-z, -x, wIntv, v1ccn, v2ccn, v1nmax, v2nmax, v1ccgn, v2ccgn)
                    zv = zv ∩ za_u
                end
                zr = zr*xr
                zv = cut(zv, zv, b.ic.v, b.ϵ_sg, sparsity(g, i), b.cut, false)
                z = zv
            end
        end
    else
        for (q,i) in enumerate(children(g, k))
            if is_num(b, i)
                znum = znum*num(b, i)
            else
                numval = false
                x = info(b, i)
                z = z*x
                z = cut(z, z, b.ic.v, b.ϵ_sg, sparsity(g, i), b.cut, false)
            end
        end
    end
    if numval
        b[k] = znum
        b._is_num[k] = true
    else
        z = z*znum
        b._info[k] = z
        b._is_num[k] = false
    end
end

for F in (PLUS, MULT, MIN, MAX)
    @eval function fprop!(t::RelaxAAInfo, v::Val{$F}, g::DAT, b::RelaxCache{MCAffPnt{N,T},N,T}, k) where {N,T<:RelaxTag}
        is_binary(g, k) ? fprop_2!(RelaxAAInfo(), Val($F), g, b, k) : fprop_n!(RelaxAAInfo(), Val($F), g, b, k)
    end
end
function fprop!(t::RelaxAAInfo, v::Val{DIV}, g::DAT, b::RelaxCache{MCAffPnt{N,T},N,T}, k) where {N,T<:RelaxTag}
    fprop_2!(RelaxAAInfo(), v, g, b, k)
end

function fprop!(t::RelaxAAInfo, v::Val{POW}, g::DAT, b::RelaxCache{MCAffPnt{N,T},N,T}, k::Int) where {N,T<:RelaxTag}
    x = child(g, 1, k)
    y = child(g, 2, k)
    if is_num(b, y) && isone(num(b, y))
        b._info[k] = info(b, x)
    elseif is_num(b,y) && iszero(num(b, y))
        b._info[k] = one(MCAffPnt{N,T})
    elseif !xy_num(b, x, y)
        if xyset(b, x, y)
            z = info(b, x)^info(b, y)
        elseif xset_ynum(b, x, y)
            z = info(b, x)^num(b, y)
        else
            z = num(b, x)^info(b, y)
        end
        b._info[k] = cut(z, info(b,k), b.ic.v, b.ϵ_sg, sparsity(g, k), b.cut, false)
        b._is_num[k] = false
    else
        b[k] = num(b, x)^num(b, y)
        b._is_num[k] = true
    end
end

function fprop!(t::RelaxAAInfo, v::Val{USER}, g::DAT, b::RelaxCache{MCAffPnt{N,T},N,T}, k::Int) where {N,T<:RelaxTag}
    f = user_univariate_operator(g, first_index(g, k))
    x = child(g, 1, k)
    if is_num(b, x)
        b[k] = f(num(b, x))
        b._is_num[k] = true
    else
        z = f(info(b, x))
        b._info[k] = cut(z, info(b, k), b.ic.v, zero(Float64), sparsity(g, k), b.cut, b.post)
        b._is_num[k] = false
    end
end

function fprop!(t::RelaxAAInfo, v::Val{USERN}, g::DAT, b::RelaxCache{MCAffPnt{N,T},N,T}, k::Int) where {N,T<:RelaxTag}
    mv = user_multivariate_operator(g, first_index(g, k))
    n = arity(g, k)
    set_input = _info_input(b, n)
    num_input = _num_input(b, n)
    anysets = false
    i = 1
    for c in children(g, k)
        if is_num(b, c)
            x = num(b, c)
            if !isinf(x)
                set_input[i] = MCAffPnt{N,T}(x)
                num_input[i] = x
            end
        else
            set_input[i] = info(b, c)
            anysets = true
        end
        i += 1
    end
    if anysets
        z = MOI.eval_objective(mv, set_input)::MCAffPnt{N,T}
        b._info[k] = cut(z, info(b, k), b.ic.v, zero(Float64), sparsity(g,k), b.cut, b.post)
        b._is_num[k] = false
    else
        b[k] = MOI.eval_objective(mv, num_input)
        b._is_num[k] = true
    end
end

for ft in UNIVARIATE_ATOM_TYPES
    f = UNIVARIATE_ATOM_DICT[ft]
    (f == :user || f == :+ || f == :-) && continue
    @eval function fprop!(t::RelaxAAInfo, v::Val{$ft}, g::DAT, b::RelaxCache{MCAffPnt{N,T},N,T}, k) where {N,T<:RelaxTag}
        x = child(g, 1, k)
        if is_num(b, x)
            b._is_num[k] = true
            return b[k] = ($f)(num(b, x))
        else
            b._is_num[k] = false
            z = ($f)(info(b, x))
            b._info[k] = cut(z, info(b, k), b.ic.v, zero(Float64), sparsity(g,k), b.cut, b.post)
        end
    end
end

for (F, f) in ((LOWER_BND, :lower_bnd), (UPPER_BND, :upper_bnd))
    @eval function fprop!(t::RelaxAAInfo, v::Val{$F}, g::DAT, b::RelaxCache{MCAffPnt{N,T},N,T}, k) where {N,T<:RelaxTag}
        x = child(g, 1, k)
        y = child(g, 2, k)
        if is_num(b, y)
            z = info(b, x)
            z = ($f)(z, num(b, y))
            b._info[k] = cut(z, info(b, k), b.ic.v, zero(Float64), sparsity(g, k), b.cut, b.post)
        end
        b._is_num[k] = b._is_num[x]
    end
end

function fprop!(t::RelaxAAInfo, v::Val{BND}, g::DAT, b::RelaxCache{MCAffPnt{N,T},N,T}, k) where {N,T<:RelaxTag}
    x = child(g, 1, k)
    z = info(b, x)
    y = child(g, 2, k)
    r = child(g, 3, k)
    if is_num(b, y) && is_num(b, r)
        z = bnd(z, num(b, y),num(b, r))
    end
    b._is_num[k] = b._is_num[x]
    b._info[k] = cut(z, info(b, k), b.ic.v, zero(Float64), sparsity(g,k), b.cut, b.post)
end

function f_init!(t::RelaxAA, g::DAT, b::RelaxCache)
    tinfo = RelaxAAInfo()
    for k = node_count(g):-1:1
        c = node_class(g, k)
        (c == EXPRESSION)    && fprop!(tinfo, Expression(), g, b, k)
        (c == VARIABLE)      && fprop!(tinfo, Variable(), g, b, k)
        (c == SUBEXPRESSION) && fprop!(tinfo, Subexpression(), g, b, k)
        k_is_num = b._is_num[k]
        b[k] = b._info[k].v
        b._is_num[k] = k_is_num
    end
    fprop!(Relax(), g, b)
    nothing
end

for d in ALL_ATOM_TYPES
    @eval function fprop!(t::RelaxAA, v::Val{$d}, g::DAT, b::RelaxCache{V,N,T}, k::Int) where {V,N,T<:RelaxTag}
        fprop!(Relax(), v, g, b, k)
    end
end

function estimator_extrema(x::MCAffPnt{N,T}, y::MCAffPnt{N,T}, s, dP) where {N,T}

    xIntv = x.box.c + sum(z -> z*UNIT_INTERVAL, x.box.γ)
    xcvU = xIntv.hi - x.box.Δ
    xccL = xIntv.lo + x.box.Δ

    yIntv = y.box.c + sum(z -> z*UNIT_INTERVAL, y.box.γ)
    ycvU = yIntv.hi - y.box.Δ
    yccL = yIntv.lo + y.box.Δ

    return xcvU, ycvU, -xccL, -yccL
end

function estimator_under(xv, yv, x::MCAffPnt{N,T}, y::MCAffPnt{N,T}, s, dP, dp, p_rel, p_diam) where {N,T}
    x_cv = x.box.c - x.box.Δ
    y_cv = y.box.c - y.box.Δ
    for (i,p_i) ∈ enumerate(s)
        rp = 2.0*p_rel[p_i]
        x_cv += x.box.γ[i]*rp
        y_cv += y.box.γ[i]*rp
    end
    x_cv_grad = SVector{N,Float64}(ntuple(i -> 2.0*x.box.γ[i].*p_diam[s[i]], Val(N)))
    y_cv_grad = SVector{N,Float64}(ntuple(i -> 2.0*y.box.γ[i].*p_diam[s[i]], Val(N)))
    x_cv, y_cv, x_cv_grad, y_cv_grad
end

function estimator_over(xv, yv, x::MCAffPnt{N,T}, y::MCAffPnt{N,T}, s, dp, dP, p_rel, p_diam) where {N,T}
    x_cc = x.box.c + x.box.Δ
    y_cc = y.box.c + y.box.Δ
    for (i,p_i) ∈ enumerate(s)
        rp = 2.0*p_rel[p_i]
        x_cc += x.box.γ[i]*rp
        y_cc += x.box.γ[i]*rp
    end
    x_ccn_grad = SVector{N,Float64}(ntuple(i -> -2.0*x.box.γ[i].*p_diam[s[i]], Val(N)))
    y_ccn_grad = SVector{N,Float64}(ntuple(i -> -2.0*y.box.γ[i].*p_diam[s[i]], Val(N)))
    -x_cc, -y_cc, x_ccn_grad, y_ccn_grad
end

#=
function extract_apriori_info(t::Union{RelaxAA,RelaxAAInfo}, v::VariableValues{Float64}, x::AffineEAGO{N}) where {N,T}
    z = 0.0
    for (k,i) in enumerate(t.v)
        l = _bd(v, i)
        u = ubd(v, i)
        z += x.γ[k]*(_val(v, i) - 0.5*(l + u))/(u - l)
    end
    Z = x.c + sum(z -> z*UNIT_INTERVAL, x.γ)
    xcv  = x.c + z - x.Δ
    xcc  = x.c + z + x.Δ
    xcvU = Z.hi - x.Δ
    xccL = Z.lo + x.Δ
    xcvg = x.γ
    xccg = x.γ
    return xcv, xcvU, xcc, xccL, xcvg, xccg
end

function _cut_info(v::VariableValues{Float64}, z::MC{N,T}, x::MCAffPnt{N,T}) where {N,T}
    xcv, xcvU, xcc, xccL, xcvg, xccg = extract_apriori_info(RelaxAA(), v, x.box)
    zaff = MC{N,T}(xcv, xcc, Interval(x.box), xcvg, xccg, false)
    return zaff ∩ z
end
function _cut_info(v::VariableValues{Float64}, z::MCAffPnt{N,T}, x::MCAffPnt{N,T}) where {N,T}
    xcv, xcvU, xcc, xccL, xcvg, xccg = extract_apriori_info(RelaxAA(), v, x.box)
    zaff = MC{N,T}(xcv, xcc, Interval(x.box), xcvg, xccg, false)
    return zaff ∩ z.v
end

function _cut(t::RelaxAAInfo, b, k, x::MCAffPnt{N,T}, z::MCAffPnt{N,T}, v::VariableValues, ϵ::Float64, s::Vector{Int}, c::Bool, p::Bool) where {N,T<:RelaxTag}
    xt = p ? MCAffPnt{N,T}(set_value_post(x.v, v, s, ϵ), x.box) : x
    xtmc = _cut_info(v, xt, xt)
    xtaffp = MCAffPnt{N,T}(xtmc, xt.box)
    b._info[k] = xtaffp
    return
end

function _cut(t::RelaxAA, b, k, x::MC{N,T}, z::MCAffPnt{N,T}, v::VariableValues, ϵ::Float64, s::Vector{Int}, c::Bool, p::Bool) where {N,T<:RelaxTag}
    xMC = set_value_post(c ? x ∩ x.Intv ∩ Interval(z.box) : x, v, s, ϵ)
    xt = p ? xMC : x
    zt = _cut_info(v, xt, info(b, k))
    _store_set!(b, zt, k)
    return
end

function fprop!(t::RelaxAAInfo, vt::Variable, g::DAT, b::RelaxCache{V,N,T}, k::Int) where {V,N,T<:RelaxTag}
    i = first_index(g, k)
    x = val(b, i)
    l = lbd(b, i)
    u = ubd(b, i)
    z = varset(MC{N,T}, rev_sparsity(g, i, k), x, x, l, u)
    zaff = AffineEAGO{N}(x, Interval(l,u), i)
    zinfo = MCAffPnt{N,T}(z, zaff)
    b._info[k] = zinfo
    return
end
function fprop!(t::RelaxAA, vt::Variable, g::DAT, b::RelaxCache{V,N,T}, k::Int) where {V,N,T<:RelaxTag}
    i = first_index(g, k)
    x = val(b, i)
    l = lbd(b, i)
    u = ubd(b, i)
    z = varset(MC{N,T}, rev_sparsity(g, i, k), x, x, l, u)
    z = z ∩ _interval(b, k)
    _store_set!(b, z, k)
end
=#

max_stars_repo_name: MatthewStuber/EAGO

max_issues_repo_name: MatthewStuber/EAGO

max_forks_repo_name: MatthewStuber/EAGO


text:
function extend_seq(points::P; scale=ReflectionScale) where P
    return if isclosed(points)
               extend_closed_seq(points)
           else
               extend_open_seq(points, scale=scale)
           end
end

function unextend_seq(points::P) where P
    pop!(points)
    popfirst!(points)
    return points
end

function extend_closed_seq(points::P) where P
    !isclosed(points) && throw(ErrorException("sequence is not closed"))
    npoints(points) < 3 && throw(ErrorException("cannot extend a sequence with fewer than 3 points"))

    pushfirst!(points, points[end-1])
    push!(points, points[3])
    return points
end

function extend_open_seq(points::P; scale=ReflectionScale) where P
    isclosed(points) && throw(ErrorException("sequence is not open"))
    npoints(points) < 3 && throw(ErrorException("cannot extend a sequence with fewer than 3 points"))

    initialpoint = pointbefore(points[1:4], scale)
    finalpoint   = pointafter(points[end-3:end], scale)
    pushfirst!(points, initialpoint)
    push!(points, finalpoint)
    return points
end

function close_seq!(points::P) where P
    isempty(points) && throw(ErrorException("cannot close an empty sequence"))
    if isapproxclosed(points)
        points[end] = points[1]
    end
    if !isclosed(points)
        push!(points, points[1])
    end
    return points
end

isclosed(firstpoint::P, lastpoint::P) where P = firstpoint == lastpoint
isclosed(points::P) where P = isclosed(first(points), last(points))

isapproxclosed(points::P) where P = isapproxclosed(first(points), last(points))

function isapproxclosed(firstpoint::P, lastpoint::P) where P
    nrm1 = norm(sqrt(eps(norm(firstpoint))), sqrt(eps(norm(lastpoint))))
    nrm2 = norm(firstpoint .- lastpoint)
    return nrm1 >= nrm2
end

npoints(pts::P) where P = length(pts)
npoints(pts::Base.Iterators.Zip) = length(pts)

ncoords(pts::P) where P = eltype(pts) <: NamedTuple ? length(Tuple(pts[1])) : length(pts[1])
ncoords(pts::Base.Iterators.Zip) = eltype(pts) <: NamedTuple ? length(Tuple(first(pts))) : length(first(pts))

coordtype(pts::P) where P = eltype(pts) <: Number ? eltype(pts) : eltype(eltype(pts))

function coordtype(::Type{T}) where {T}
    result = T
    !(result<:Number) ? (result = eltype(result)) : (return result)
    !(result<:Number) ? (result = eltype(result)) : (return result)
    !(result<:Number) ? (result = eltype(result)) : (return result)
    !(result<:Number) ? throw(ErrorException("unable to discern the coordinate type for $T")) : (return result)
end

Base.convert(::Type{Array{T,1}}, x::NTuple{N,T}) where {N,T} = [x...,]
Base.convert(::Type{NTuple{N,T}}, x::Array{T,1}) where {N,T} = (x...,)

function coords_to_cols(coords...)
    n_cols = length(coords)
    n_rows = length(coords[1])
    return reshape(vcat(coords...,), n_rows, n_cols)
end

cols_to_coords(cols::Array{T,2}) where T = [cols[:,i] for i=1:(size(cols)[2])]

coords_to_points(coords...) = collect(zip(coords...,))

function points_to_coords(x)
    n_dims  = length(x[1])
    n_items = length(x)
    typ = eltype(eltype(x[1]))
    res = Array{typ, 2}(undef, (n_items, n_dims))
    for k=1:n_points
       for i=1:n_dims
          res[k,i] = x[k][i]
       end
    end
    return res
end

cols_to_points(x) = coords_to_points(cols_to_coords(x))
points_to_cols(x) = coords_to_cols(points_to_coords(x))

ext: jl


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct IoK8sApiExtensionsV1beta1NetworkPolicyPort <: SwaggerModel
    port::Any # spec type: Union{ Nothing, IoK8sApimachineryPkgUtilIntstrIntOrString } # spec name: port
    protocol::Any # spec type: Union{ Nothing, String } # spec name: protocol

    function IoK8sApiExtensionsV1beta1NetworkPolicyPort(;port=nothing, protocol=nothing)
        o = new()
        validate_property(IoK8sApiExtensionsV1beta1NetworkPolicyPort, Symbol("port"), port)
        setfield!(o, Symbol("port"), port)
        validate_property(IoK8sApiExtensionsV1beta1NetworkPolicyPort, Symbol("protocol"), protocol)
        setfield!(o, Symbol("protocol"), protocol)
        o
    end
end # type IoK8sApiExtensionsV1beta1NetworkPolicyPort

const _property_map_IoK8sApiExtensionsV1beta1NetworkPolicyPort = Dict{Symbol,Symbol}(Symbol("port")=>Symbol("port"), Symbol("protocol")=>Symbol("protocol"))
const _property_types_IoK8sApiExtensionsV1beta1NetworkPolicyPort = Dict{Symbol,String}(Symbol("port")=>"IoK8sApimachineryPkgUtilIntstrIntOrString", Symbol("protocol")=>"String")
Base.propertynames(::Type{ IoK8sApiExtensionsV1beta1NetworkPolicyPort }) = collect(keys(_property_map_IoK8sApiExtensionsV1beta1NetworkPolicyPort))
Swagger.property_type(::Type{ IoK8sApiExtensionsV1beta1NetworkPolicyPort }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiExtensionsV1beta1NetworkPolicyPort[name]))}
Swagger.field_name(::Type{ IoK8sApiExtensionsV1beta1NetworkPolicyPort }, property_name::Symbol) =  _property_map_IoK8sApiExtensionsV1beta1NetworkPolicyPort[property_name]

function check_required(o::IoK8sApiExtensionsV1beta1NetworkPolicyPort)
    true
end

function validate_property(::Type{ IoK8sApiExtensionsV1beta1NetworkPolicyPort }, name::Symbol, val)
end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

using Test, Random

######## Utilities ###########

# Generate an array similar to A[indx1, indx2, ...], but only call
# getindex with scalar-valued indices. This will be safe even if
# `getindex` someday calls `view`.

# The "nodrop" variant does not drop any dimensions (not even trailing ones)
function Agen_nodrop(A::AbstractArray, I...)
    irep = replace_colon(A, I)
    _Agen(A, ensure_iterable(irep)...)
end

# This drops scalar dimensions
function Agen_slice(A::AbstractArray, I...)
    irep = replace_colon(A, I)
    B = _Agen(A, ensure_iterable(irep)...)
    sd = Int[]
    for i = 1:length(I)
        if isa(I[i], Real)
            push!(sd, i)
        end
    end
    squeeze(B, dims=sd)
end

_Agen(A, i1) = [A[j1] for j1 in i1]
_Agen(A, i1, i2) = [A[j1,j2] for j1 in i1, j2 in i2]
_Agen(A, i1, i2, i3) = [A[j1,j2,j3] for j1 in i1, j2 in i2, j3 in i3]
_Agen(A, i1, i2, i3, i4) = [A[j1,j2,j3,j4] for j1 in i1, j2 in i2, j3 in i3, j4 in i4]
_Agen(A, i1, i2, i3, i4, i5) = [A[j1,j2,j3,j4,j5] for j1 in i1, j2 in i2, j3 in i3, j4 in i4, j5 in i5]
_Agen(A, i1, i2, i3, i4, i5, i6) = [A[j1,j2,j3,j4,j5,j6] for j1 in i1, j2 in i2, j3 in i3, j4 in i4, j5 in i5, j6 in i6]

function replace_colon(A::AbstractArray, I)
    Iout = Vector{Any}(undef, length(I))
    I == (:,) && return (1:length(A),)
    for d = 1:length(I)
        Iout[d] = isa(I[d], Colon) ? (1:size(A,d)) : I[d]
    end
    (Iout...,)
end

ensure_iterable(::Tuple{}) = ()
ensure_iterable(t::Tuple{Union{Number, CartesianIndex}, Vararg{Any}}) = ((t[1],), ensure_iterable(Base.tail(t))...)
ensure_iterable(t::Tuple{Any, Vararg{Any}}) = (t[1], ensure_iterable(Base.tail(t))...)

index_ndims(t::Tuple) = tup2val(Base.index_ndims(t))
tup2val(::NTuple{N}) where {N} = Val(N)

# To avoid getting confused by manipulations that are implemented for SubArrays,
# it's good to copy the contents to an Array. This version protects against
# `similar` ever changing its meaning.
function copy_to_array(A::AbstractArray)
    Ac = Array{eltype(A)}(undef, size(A))
    copyto!(Ac, A)
end

# Discover the highest dimension along which the values in A are
# separated by a single increment.  If A was extracted via getindex
# from reshape(1:N, ...), this is equivalent to finding the highest
# dimension of the SubArray consistent with a single stride in the
# parent array.
function single_stride_dim(A::Array)
    ld = 0
    while ld < ndims(A)
        # Collapse all dimensions up to & including ld+1 into the first dimension
        shp = [prod(size(A)[1:ld+1])]
        for j = ld+2:ndims(A)
            push!(shp, size(A,j))
        end
        Ar = reshape(A, shp...)
        # Compute the diff along dimension 1
        if size(Ar, 1) > 1
            indices = map(d->1:size(Ar,d), [1:ndims(Ar);])
            indicesp = copy(indices); indicesp[1] = 2:size(Ar,1)
            indicesm = copy(indices); indicesm[1] = 1:size(Ar,1)-1
            dA = Ar[indicesp...] - Ar[indicesm...]
            ustride = unique(dA[:])
            if length(ustride) == 1  # is it a single stride?
                ld += 1
            else
                break
            end
        else
            ld += 1
        end
    end
    ld
end
single_stride_dim(@nospecialize(A)) = single_stride_dim(copy_to_array(A))

# Testing equality of AbstractArrays, using several different methods to access values
function test_cartesian(@nospecialize(A), @nospecialize(B))
    isgood = true
    for (IA, IB) in zip(eachindex(A), eachindex(B))
        if A[IA] != B[IB]
            @show IA IB A[IA] B[IB]
            isgood = false
            break
        end
    end
    if !isgood
        @show A
        @show B
        error("Mismatch")
    end
end

function test_linear(@nospecialize(A), @nospecialize(B))
    length(A) == length(B) || error("length mismatch")
    isgood = true
    for (iA, iB) in zip(1:length(A), 1:length(B))
        if A[iA] != B[iB]
            isgood = false
            break
        end
    end
    if !isgood
        @show A
        @show A.indices
        @show B
        error("Mismatch")
    end
end

# "mixed" means 2 indices even for N-dimensional arrays
test_mixed(::AbstractArray{T,1}, ::Array) where {T} = nothing
test_mixed(::AbstractArray{T,2}, ::Array) where {T} = nothing
test_mixed(A, B::Array) = _test_mixed(A, reshape(B, size(A)))
function _test_mixed(@nospecialize(A), @nospecialize(B))
    m = size(A, 1)
    n = size(A, 2)
    isgood = true
    for J in CartesianIndices(size(A)[2:end]), i in 1:m
        if A[i,J] != B[i,J]
            isgood = false
            break
        end
    end
    if !isgood
        @show A
        @show B
        error("Mismatch")
    end
    nothing
end

function test_bounds(@nospecialize(A))
    @test_throws BoundsError A[0]
    @test_throws BoundsError A[end+1]
    trailing2 = ntuple(x->1, max(ndims(A)-2, 0))
    trailing3 = ntuple(x->1, max(ndims(A)-3, 0))
    @test_throws BoundsError A[1, 0, trailing2...]
    @test_throws BoundsError A[1, end+1, trailing2...]
    @test_throws BoundsError A[1, 1, 0, trailing3...]
    @test_throws BoundsError A[1, 1, end+1, trailing3...]
    @test_throws BoundsError A[0, 1, trailing2...]
    @test_throws BoundsError A[end+1, 1, trailing2...]
    @test_throws BoundsError A[0, 1, 1, trailing3...]
    @test_throws BoundsError A[end+1, 1, 1, trailing3...]
    @test_throws BoundsError A[1, 0, 1, trailing3...]
    @test_throws BoundsError A[1, end+1, 1, trailing3...]
    # TODO: PLI (re-enable after 0.7)
    # @test_throws BoundsError A[1, 0]
    # @test_throws BoundsError A[1, end+1]
    # @test_throws BoundsError A[1, 1, 0]
    # @test_throws BoundsError A[1, 1, end+1]
    # @test_throws BoundsError A[0, 1]
    # @test_throws BoundsError A[end+1, 1]
    # @test_throws BoundsError A[0, 1, 1]
    # @test_throws BoundsError A[end+1, 1, 1]
    # @test_throws BoundsError A[1, 0, 1]
    # @test_throws BoundsError A[1, end+1, 1]
end

function dim_break_linindex(I)
    i = 1
    while i <= length(I) && !isa(I[i], Vector{Int})
        i += 1
    end
    i - 1
end

function runsubarraytests(A::Array, I...)
    # Direct test of linear indexing inference
    C = Agen_nodrop(A, I...)
    ld = min(single_stride_dim(C), dim_break_linindex(I))
    S = view(A, I...)
    if Base.iscontiguous(S)
        @test S.stride1 == 1
    end
    test_linear(S, C)
    test_cartesian(S, C)
    test_mixed(S, C)
end

function runsubarraytests(@nospecialize(A), I...)
    # When A was created with view, we have to check bounds, since some
    # of the "residual" dimensions have size 1. It's possible that we
    # need dedicated tests for view.
    for d = 1:length(I)-1
        if !isa(I[d], Colon) && any(I[d] .> size(A,d))
            return nothing
        end
    end
    if !isa(I[end], Colon) && any(I[end] .> prod(size(A)[length(I):end]))
        return nothing
    end
    AA = copy_to_array(A)
    # Direct test of linear indexing inference
    C = Agen_nodrop(AA, I...)
    Cld = ld = min(single_stride_dim(C), dim_break_linindex(I))
    Cdim = AIindex = 0
    while Cdim <= Cld && AIindex < length(A.indices)
        AIindex += 1
        if isa(A.indices[AIindex], Real)
            ld += 1
        else
            Cdim += 1
        end
    end
    local S
    try
        S = view(A, I...)
    catch err
        @show typeof(A)
        @show A.indices
        @show I
        rethrow(err)
    end
    test_linear(S, C)
    test_cartesian(S, C)
    test_mixed(S, C)
end

# indexN is a cartesian index, indexNN is a linear index for 2 dimensions, and indexNNN is a linear index for 3 dimensions
function runviews(SB::AbstractArray, indexN, indexNN, indexNNN)
    @assert ndims(SB) > 2
    for i3 in indexN, i2 in indexN, i1 in indexN
        runsubarraytests(SB, i1, i2, i3, ntuple(x->1, max(ndims(SB)-3, 0))...)
    end
    for i2 in indexN, i1 in indexN
        runsubarraytests(SB, i1, i2, ntuple(x->1, max(ndims(SB)-2, 0))...)
    end
    for i1 in indexNNN
        runsubarraytests(SB, i1)
    end
end

function runviews(SB::AbstractArray{T, 3} where T, indexN, indexNN, indexNNN)
    @assert ndims(SB) > 2
    for i3 in indexN, i2 in indexN, i1 in indexN
        runsubarraytests(SB, i1, i2, i3)
    end
    for i2 in indexN, i1 in indexN
        runsubarraytests(SB, i1, i2, 1)
    end
    for i1 in indexNNN
        runsubarraytests(SB, i1)
    end
end

function runviews(SB::AbstractArray{T,2}, indexN, indexNN, indexNNN) where T
    for i2 in indexN, i1 in indexN
        runsubarraytests(SB, i1, i2)
    end
    for i1 in indexNN
        runsubarraytests(SB, i1)
    end
end

function runviews(SB::AbstractArray{T,1}, indexN, indexNN, indexNNN) where T
    for i1 in indexN
        runsubarraytests(SB, i1)
    end
end

runviews(SB::AbstractArray{T,0}, indexN, indexNN, indexNNN) where {T} = nothing

######### Tests #########

testfull = Bool(parse(Int,(get(ENV, "JULIA_TESTFULL", "0"))))

### Views from Arrays ###
index5 = (1, :, 2:5, [4,1,5], reshape([2]), view(1:5,[2 3 4 1]))  # all work with at least size 5
index25 = (3, :, 2:11, [19,9,7], reshape([10]), view(1:25,[19 15; 4 24]))
index125 = (113, :, 85:121, [99,14,103], reshape([72]), view(1:125,reshape([25,4,102,67], 1, 2, 2)))

if testfull
    let A = copy(reshape(1:5*7*11, 11, 7, 5))
        runviews(A, index5, index25, index125)
    end
end

### Views from views ###

# "outer" indices create snips that have at least size 5 along each dimension,
# with the exception of Int-slicing
oindex = (:, 6, 3:7, reshape([12]), [8,4,6,12,5,7], [3:7 1:5 2:6 4:8 5:9])

_ndims(::AbstractArray{T,N}) where {T,N} = N
_ndims(x) = 1

if testfull
    let B = copy(reshape(1:13^3, 13, 13, 13))
        for o3 in oindex, o2 in oindex, o1 in oindex
            viewB = view(B, o1, o2, o3)
            runviews(viewB, index5, index25, index125)
        end
    end
end

if !testfull
    let B = copy(reshape(1:13^3, 13, 13, 13))
        for oind in ((:,:,:),
                     (:,:,6),
                     (:,6,:),
                     (6,:,:),
                     (:,3:7,:),
                     (3:7,:,:),
                     (3:7,6,:),
                     (3:7,6,0x6),
                     (6,UInt(3):UInt(7),3:7),
                     (13:-2:1,:,:),
                     ([8,4,6,12,5,7],:,3:7),
                     (6,CartesianIndex.(6,[8,4,6,12,5,7])),
                     (CartesianIndex(13,6),[8,4,6,12,5,7]),
                     (1,:,view(1:13,[9,12,4,13,1])),
                     (view(1:13,[9,12,4,13,1]),2:6,4),
                     ([1:5 2:6 3:7 4:8 5:9], :, 3))
            runsubarraytests(B, oind...)
            viewB = view(B, oind...)
            runviews(viewB, index5, index25, index125)
        end
    end
end

# issue #11289
x11289 = randn(5,5)
@test isempty(view(x11289, Int[], :))
@test isempty(view(x11289, [2,5], Int[]))
@test isempty(view(x11289, Int[], 2))

####### "Classical" tests #######

# Tests where non-trailing dimensions are preserved
A = copy(reshape(1:120, 3, 5, 8))
sA = view(A, 2:2, 1:5, :)
@test @inferred(strides(sA)) == (1, 3, 15)
@test parent(sA) == A
@test parentindices(sA) == (2:2, 1:5, Base.Slice(1:8))
@test Base.parentdims(sA) == [1:3;]
@test size(sA) == (1, 5, 8)
@test axes(sA) === (Base.OneTo(1), Base.OneTo(5), Base.OneTo(8))
@test sA[1, 2, 1:8][:] == [5:15:120;]
sA[2:5:end] .= -1
@test all(sA[2:5:end] .== -1)
@test all(A[5:15:120] .== -1)
@test @inferred(strides(sA)) == (1,3,15)
@test stride(sA,3) == 15
@test stride(sA,4) == 120
test_bounds(sA)
sA = view(A, 1:3, 1:5, 5)
@test Base.parentdims(sA) == [1:2;]
sA[1:3,1:5] .= -2
@test all(A[:,:,5] .== -2)
fill!(sA, -3)
@test all(A[:,:,5] .== -3)
sA[:] .= 4
@test all(A[:,:,5] .== 4)
@test @inferred(strides(sA)) == (1,3)
test_bounds(sA)
sA = view(A, 1:3, 3:3, 2:5)
@test Base.parentdims(sA) == [1:3;]
@test size(sA) == (3,1,4)
@test axes(sA) === (Base.OneTo(3), Base.OneTo(1), Base.OneTo(4))
@test sA == A[1:3,3:3,2:5]
@test sA[:] == A[1:3,3,2:5][:]
test_bounds(sA)
sA = view(A, 1:2:3, 1:3:5, 1:2:8)
@test Base.parentdims(sA) == [1:3;]
@test @inferred(strides(sA)) == (2,9,30)
@test sA[:] == A[1:2:3, 1:3:5, 1:2:8][:]
# issue #8807
@test view(view([1:5;], 1:5), 1:5) == [1:5;]
# Test with mixed types
@test sA[:, Int16[1,2], big(2)] == [31 40; 33 42]
test_bounds(sA)
sA = view(A, 1:1, 1:5, [1 3; 4 2])
@test ndims(sA) == 4
@test axes(sA) === (Base.OneTo(1), Base.OneTo(5), Base.OneTo(2), Base.OneTo(2))
sA = view(A, 1:2, 3, [1 3; 4 2])
@test ndims(sA) == 3
@test axes(sA) === (Base.OneTo(2), Base.OneTo(2), Base.OneTo(2))

# logical indexing #4763
A = view([1:10;], 5:8)
@test A[A.<7] == view(A, A.<7) == [5, 6]
@test Base.unsafe_getindex(A, A.<7) == [5, 6]
B = reshape(1:16, 4, 4)
sB = view(B, 2:3, 2:3)
@test sB[sB.>8] == view(sB, sB.>8) == [10, 11]
@test Base.unsafe_getindex(sB, sB.>8) == [10, 11]

# Tests where dimensions are dropped
A = copy(reshape(1:120, 3, 5, 8))
sA = view(A, 2, :, 1:8)
@test parent(sA) == A
@test parentindices(sA) == (2, Base.Slice(1:5), 1:8)
@test Base.parentdims(sA) == [2:3;]
@test size(sA) == (5, 8)
@test axes(sA) === (Base.OneTo(5), Base.OneTo(8))
@test @inferred(strides(sA)) == (3,15)
@test sA[2, 1:8][:] == [5:15:120;]
@test sA[:,1] == [2:3:14;]
@test sA[2:5:end] == [5:15:110;]
sA[2:5:end] .= -1
@test all(sA[2:5:end] .== -1)
@test all(A[5:15:120] .== -1)
test_bounds(sA)
sA = view(A, 1:3, 1:5, 5)
@test Base.parentdims(sA) == [1:2;]
@test size(sA) == (3,5)
@test axes(sA) === (Base.OneTo(3),Base.OneTo(5))
@test @inferred(strides(sA)) == (1,3)
test_bounds(sA)
sA = view(A, 1:2:3, 3, 1:2:8)
@test Base.parentdims(sA) == [1,3]
@test size(sA) == (2,4)
@test axes(sA) === (Base.OneTo(2), Base.OneTo(4))
@test @inferred(strides(sA)) == (2,30)
@test sA[:] == A[sA.indices...][:]
test_bounds(sA)

let a = [5:8;]
    @test parent(a) == a
    @test parentindices(a) == (1:4,)
end

# issue #6218 - logical indexing
A = rand(2, 2, 3)
msk = fill(true, 2, 2)
msk[2,1] = false
sA = view(A, :, :, 1)
sA[msk] .= 1.0
@test sA[msk] == fill(1, count(msk))

# bounds checking upon construction; see #4044, #10296
@test_throws BoundsError view(1:10, 8:11)
A = reshape(1:20, 5, 4)
sA = view(A, 1:2, 1:3)
@test_throws BoundsError view(sA, 1:3, 1:3)
@test_throws BoundsError view(sA, 1:2, 1:4)
view(sA, 1:2, 1:2)
@test_throws BoundsError view(A, 17:23)
view(A, 17:20)

# Linear indexing by one multidimensional array:
A = reshape(1:120, 3, 5, 8)
sA = view(A, :, :, :)
@test sA[[72 17; 107 117]] == [72 17; 107 117]
@test sA[[99 38 119 14 76 81]] == [99 38 119 14 76 81]
@test sA[[fill(1, (2, 2, 2)); fill(2, (2, 2, 2))]] == [fill(1, (2, 2, 2)); fill(2, (2, 2, 2))]
sA = view(A, 1:2, 2:3, 3:4)
@test sA[(1:8)'] == [34 35 37 38 49 50 52 53]
@test sA[[1 2 4 4; 6 1 1 4]] == [34 35 38 38; 50 34 34 38]

# issue #11871
let a = fill(1., (2,2)),
    b = view(a, 1:2, 1:2)
    b[2] = 2
    @test b[2] === 2.0
end

# issue #15138
let a = [1,2,3],
    b = view(a, UInt(1):UInt(2))
    @test b == view(a, UInt(1):UInt(2)) == view(view(a, :), UInt(1):UInt(2)) == [1,2]
end

let A = reshape(1:4, 2, 2),
    B = view(A, :, :)
    @test parent(B) === A
    @test parent(view(B, 0x1, :)) === parent(view(B, 0x1, :)) === A
end

# issue #15168
let A = rand(10), sA = view(copy(A), :)
    @test sA[Int16(1)] === sA[Int32(1)] === sA[Int64(1)] === A[1]
    permute!(sA, Vector{Int16}(1:10))
    @test A == sA
end

# the following segfaults with LLVM 3.8 on Windows, ref #15417
@test Array(view(view(reshape(1:13^3, 13, 13, 13), 3:7, 6:6, :), 1:2:5, :, 1:2:5)) ==
    cat([68,70,72],[406,408,410],[744,746,748]; dims=3)

# tests @view (and replace_ref_end!)
X = reshape(1:24,2,3,4)
Y = 4:-1:1

@test isa(@view(X[1:3]), SubArray)

@test X[1:end] == @.(@view X[1:end]) # test compatibility of @. and @view
@test X[1:end-3] == @view X[1:end-3]
@test X[1:end,2,2] == @view X[1:end,2,2]
@test X[1,1:end-2,1] == @view X[1,1:end-2,1]
@test X[1,2,1:end-2] == @view X[1,2,1:end-2]
@test X[1,2,Y[2:end]] == @view X[1,2,Y[2:end]]
@test X[1:end,2,Y[2:end]] == @view X[1:end,2,Y[2:end]]

u = (1,2:3)
@test X[u...,2:end] == @view X[u...,2:end]
@test X[(1,)...,(2,)...,2:end] == @view X[(1,)...,(2,)...,2:end]

# test macro hygiene
let size=(x,y)-> error("should not happen"), Base=nothing
    @test X[1:end,2,2] == @view X[1:end,2,2]
end

# test that side effects occur only once
let foo = [X]
    @test X[2:end-1] == @view (push!(foo,X)[1])[2:end-1]
    @test foo == [X, X]
end

# test @views macro
@views let f!(x) = x[1:end-1] .+= x[2:end].^2
    x = [1,2,3,4]
    f!(x)
    @test x == [5,11,19,4]
    @test x[1:3] isa SubArray
    @test x[2] === 11
    @test Dict((1:3) => 4)[1:3] === 4
    x[1:2] .= 0
    @test x == [0,0,19,4]
    x[1:2] .= 5:6
    @test x == [5,6,19,4]
    f!(x[3:end])
    @test x == [5,6,35,4]
    x[Y[2:3]] .= 7:8
    @test x == [5,8,7,4]
    x[(3,)..., ()...] += 3
    @test x == [5,8,10,4]
    i = Int[]
    # test that lhs expressions in update operations are evaluated only once:
    x[push!(i,4)[1]] += 5
    @test x == [5,8,10,9] && i == [4]
    x[push!(i,3)[end]] += 2
    @test x == [5,8,12,9] && i == [4,3]
    @. x[3:end] = 0       # make sure @. works with end expressions in @views
    @test x == [5,8,0,0]
end
@views @test isa(X[1:3], SubArray)
@test X[1:end] == @views X[1:end]
@test X[1:end-3] == @views X[1:end-3]
@test X[1:end,2,2] == @views X[1:end,2,2]
@test X[1,2,1:end-2] == @views X[1,2,1:end-2]
@test X[1,2,Y[2:end]] == @views X[1,2,Y[2:end]]
@test X[1:end,2,Y[2:end]] == @views X[1:end,2,Y[2:end]]
@test X[u...,2:end] == @views X[u...,2:end]
@test X[(1,)...,(2,)...,2:end] == @views X[(1,)...,(2,)...,2:end]
# test macro hygiene
let size=(x,y)-> error("should not happen"), Base=nothing
    @test X[1:end,2,2] == @views X[1:end,2,2]
end

# issue #18034
# ensure that it is possible to create an isbits, IndexLinear view of an immutable Array
let
    struct ImmutableTestArray{T, N} <: Base.DenseArray{T, N}
    end
    Base.size(::Union{ImmutableTestArray, Type{ImmutableTestArray}}) = (0, 0)
    Base.IndexStyle(::Union{ImmutableTestArray, Type{ImmutableTestArray}}) = Base.IndexLinear()
    a = ImmutableTestArray{Float64, 2}()
    @test Base.IndexStyle(view(a, :, :)) == Base.IndexLinear()
    @test isbits(view(a, :, :))
end

# ref issue #17351
@test @inferred(reverse(view([1 2; 3 4], :, 1), dims=1)) == [3, 1]

let
    s = view(reshape(1:6, 2, 3), 1:2, 1:2)
    @test @inferred(s[2,2,1]) === 4
end

# issue #18581: slices with OneTo axes can be linear
let
    A18581 = rand(5, 5)
    B18581 = view(A18581, :, axes(A18581,2))
    @test IndexStyle(B18581) === IndexLinear()
end

@test sizeof(view(zeros(UInt8, 10), 1:4)) == 4
@test sizeof(view(zeros(UInt8, 10), 1:3)) == 3
@test sizeof(view(zeros(Float64, 10, 10), 1:3, 2:6)) == 120

# PR #25321
# checks that issue in type inference is resolved
A = rand(5,5,5,5)
V = view(A, 1:1 ,:, 1:3, :)
@test @inferred(strides(V)) == (1, 5, 25, 125)

# Issue #26263 — ensure that unaliascopy properly trims the array
A = rand(5,5,5,5)
V = view(A, 2:5, :, 2:5, 1:2:5)
@test @inferred(Base.unaliascopy(V)) == V == A[2:5, :, 2:5, 1:2:5]
@test @inferred(sum(Base.unaliascopy(V))) == sum(V) == sum(A[2:5, :, 2:5, 1:2:5])

# issue #27632
function _test_27632(A)
    for J in CartesianIndices(size(A)[2:end])
        A[1, J]
    end
    nothing
end
# check that this doesn't crash
_test_27632(view(ones(Int64, (1, 1, 1)), 1, 1, 1))

max_stars_repo_name: perrutquist/julia

max_issues_repo_name: perrutquist/julia

max_forks_repo_name: perrutquist/julia


text:
using DashBootstrapComponents

input_group = dbc_inputgroup([
    dbc_button("Random name", id = "input-group-button", n_clicks = 0),
    dbc_input(id = "input-group-button-input", placeholder = "name"),
]);

callback!(
    app,
    Output("input-group-button-input", "value"),
    Input("input-group-button", "n_clicks"),
) do n_clicks
    if n_clicks > 0
        names = ["Arthur Dent", "Ford Prefect", "Trillian Astra"]
        which = n_clicks % length(names)
        return names[which+1]
    else
        return ""
    end
end;

max_stars_repo_name: glsdown/dash-bootstrap-components

max_issues_repo_name: glsdown/dash-bootstrap-components

max_forks_repo_name: glsdown/dash-bootstrap-components


text:
function transform(input::AbstractArray)

end


text:
""" Given two sorted vectors x and y of length nx and ny,
    return a sorted vector of length nx+ny containing all elements
    from x and y.
    Equivalent to sort([x; y]) but hopefully a little bit more efficent.
"""
function merge_sorted(x, y)
    nx = length(x)
    ny = length(y)
    merged = Vector{eltype(x)}(undef, nx+ny)
    ix=1
    iy=1
    for ixy in 1:(nx+ny)
        if ix > nx
            # we've run out of x, just fill the rest of the vector with y
            merged[ixy] = y[iy]
            iy +=1
        elseif iy > ny
            # we've run out of y, just fill the rest of the vector with x
            merged[ixy] = x[ix]
            ix +=1
        elseif x[ix] < y[iy]
            merged[ixy] = x[ix]
            ix +=1
        else
            merged[ixy] = y[iy]
            iy +=1
        end
    end
    return merged
end

""" Find all the points where the PDFs of the two distributions in a coupling
    cross. Also return which PDF dominates in each segment.
"""
function get_crossings(coup::MaximalCoupling, nquantiles::Int)
    # Start by comparing the PDFs on a grid of points given by the
    # merged quantiles of the two distributions.
    pquant = quantile.(coup.p, range(0, stop=1, length=nquantiles)[2:end-1])
    qquant = quantile.(coup.q, range(0, stop=1, length=nquantiles)[2:end-1])
    pq_quantiles = merge_sorted(pquant, qquant)

    # function that returns the difference of the two PDFs
    pq_diff(x) = pdf(coup.p, x) - pdf(coup.q, x)

    # p(x) greater than q(x)
    p_gt_q = pq_diff(pq_quantiles[1]) >= 0
    xbefore = -Inf
    crossed_x = Float64[]
    pgreater = Bool[p_gt_q]
    for (idx, x) in enumerate(pq_quantiles)
        x = pq_quantiles[idx]
        if (
            (p_gt_q) # currently expect p(x) > q(x)
            && pq_diff(x) < 0 # but see p(x) < q(x)
           ) || ( # or
            !(p_gt_q)
            && pq_diff(x) > 0
           )
            # p&q have crossed!
            # push!(crossed_indx, idx)
            p_gt_q = !p_gt_q # switch record of dominant PDF

            push!(crossed_x, fzero(pq_diff, xbefore, x))
            push!(pgreater, p_gt_q)
        end
        xbefore = x
    end

    return crossed_x, pgreater
end

""" Probability that X=Y (coupling event) in a maximal coupling
    of two univariate distributions.
"""
function prob_couple(coup::MaximalCoupling; nquantiles::Int=100)
    crossed_x, pgreater = get_crossings(coup, nquantiles)
    xbefore = -Inf # could use minimum support of p and q instead
    pcouple = 0.0
    for (ichange, xcross) in enumerate(crossed_x)
        if pgreater[ichange]
            pcouple += cdf(coup.q, xcross) - cdf(coup.q, xbefore)
        else
            pcouple += cdf(coup.p, xcross) - cdf(coup.p, xbefore)
        end
        xbefore = xcross
    end
    if pgreater[end]
        pcouple += cdf(coup.q, Inf) - cdf(coup.q, xbefore)
    else
        pcouple += cdf(coup.p, Inf) - cdf(coup.p, xbefore)
    end
    @assert 0.0 <= pcouple <= 1.0
    return pcouple
end
function prob_couple(coup::MaximalCoupling{U1,U2}; nquantiles::Int=100) where {U1<:DiscreteDistribution, U2<:DiscreteDistribution}
    p_finite_support = isfinite(minimum(coup.p)) && isfinite(maximum(coup.p))
    q_finite_support = isfinite(minimum(coup.q)) && isfinite(maximum(coup.q))
    local overlap
    if p_finite_support && q_finite_support
        p_support = support(coup.p)
        q_support = support(coup.q)
        overlap = p_support ∩ q_support
    elseif p_finite_support
        overlap = support(coup.p)
    elseif q_finite_support
        overlap = support(coup.q)
    else
        approx_min = min(quantile(coup.p, 1e-10), quantile(coup.q, 1e-10))
        approx_max = max(quantile(coup.p, 1-1e-10), quantile(coup.q, 1-1e-10))
        overlap = approx_min:approx_max
        @assert length(overlap) < 10^5 # otherwise it'll take too long
        # if this is an issue, we should consider using an adaptation
        # of the crossings-detection method for continuous distributions
    end
    # prob_min = [min(pdf(coup.p, x), pdf(coup.q, x)) for x in overlap]
    # pcouple = sum(prob_min)
    pcouple = sum(x -> min(pdf(coup.p, x), pdf(coup.q, x)), overlap)
    return pcouple
end


text:
function variable_converter_filter_voltage(pm::AbstractWModels; kwargs...)
    variable_converter_filter_voltage_magnitude_sqr(pm; kwargs...)
    variable_converter_filter_voltage_cross_products(pm; kwargs...)
end

function variable_converter_internal_voltage(pm::AbstractWModels; kwargs...)
    variable_converter_internal_voltage_magnitude_sqr(pm; kwargs...)
    variable_converter_internal_voltage_cross_products(pm; kwargs...)
end
"""
Creates lossy converter model between AC and DC side
```
pconv_ac[i] + pconv_dc[i] == a + b*I + c*Isq
```
"""
function constraint_converter_losses(pm::AbstractWModels, n::Int, cnd::Int, i::Int, a, b, c, plmax)
    pconv_ac = PowerModels.var(pm, n, cnd, :pconv_ac, i)
    pconv_dc = PowerModels.var(pm, n, cnd, :pconv_dc, i)
    iconv = PowerModels.var(pm, n, cnd, :iconv_ac, i)
    iconv_sq = PowerModels.var(pm, n, cnd, :iconv_ac_sq, i)

    PowerModels.con(pm, n, cnd, :conv_loss)[i] = @constraint(pm.model, pconv_ac + pconv_dc == a + b*iconv + c*iconv_sq)
end
"""
Converter transformer constraints
```
p_tf_fr ==  g/(tm^2)*w_fr + -g/(tm)*wr + -b/(tm)*wi)
q_tf_fr == -b/(tm^2)*w_fr +  b/(tm)*wr + -g/(tm)*wi)
p_tf_to ==  g*w_to + -g/(tm)*wr     + -b/(tm)*(-wi))
q_tf_to == -b*w_to +  b/(tm)*wr     + -g/(tm)*(-wi))
```
"""
function constraint_conv_transformer(pm::AbstractWRModels, n::Int, cnd::Int, i::Int, rtf, xtf, acbus, tm, transformer)
    ptf_fr = PowerModels.var(pm, n, cnd, :pconv_tf_fr, i)
    qtf_fr = PowerModels.var(pm, n, cnd, :qconv_tf_fr, i)
    ptf_to = PowerModels.var(pm, n, cnd, :pconv_tf_to, i)
    qtf_to = PowerModels.var(pm, n, cnd, :qconv_tf_to, i)

    w = PowerModels.var(pm, n, cnd, :w, acbus)  # vm^2
    wf = PowerModels.var(pm, n, cnd, :wf_ac, i)   # vmf * vmf
    wrf = PowerModels.var(pm, n, cnd, :wrf_ac, i) # vm*vmf*cos(va-vaf) =  vmf*vm*cos(vaf-va)
    wif = PowerModels.var(pm, n, cnd, :wif_ac, i) # vm*vmf*sin(va-vaf) = -vmf*vm*sin(vaf-va)

    ztf = rtf + im*xtf
    if transformer
        ytf = 1/(rtf + im*xtf)
        gtf = real(ytf)
        btf = imag(ytf)
        c1, c2, c3, c4 = ac_power_flow_constraints_w(pm, gtf, btf, w, wf, wrf, wif, ptf_fr, ptf_to, qtf_fr, qtf_to, tm)
        PowerModels.con(pm, n, cnd, :conv_tf_p_fr)[i] = c1
        PowerModels.con(pm, n, cnd, :conv_tf_q_fr)[i] = c2
        PowerModels.con(pm, n, cnd, :conv_tf_p_to)[i] = c3
        PowerModels.con(pm, n, cnd, :conv_tf_q_to)[i] = c4

        #@constraint(pm.model, (wrf)^2 + (wif)^2 <= wf*w)
        constraint_voltage_product_converter(pm, wrf, wif, wf, w)
    else
        pcon, qcon = constraint_lossless_section(pm, w, wf, wrf, wif, ptf_fr, ptf_to, qtf_fr, qtf_to)
    end
end
"constraints for a voltage magnitude transformer + series impedance"
function ac_power_flow_constraints_w(pm::AbstractWRModels, g, b, w_fr, w_to, wr, wi, p_fr, p_to, q_fr, q_to, tm)
    c1 = @constraint(pm.model, p_fr ==  g/(tm^2)*w_fr + -g/(tm)*wr + -b/(tm)*wi)
    c2 = @constraint(pm.model, q_fr == -b/(tm^2)*w_fr +  b/(tm)*wr + -g/(tm)*wi)
    c3 = @constraint(pm.model, p_to ==  g*w_to + -g/(tm)*wr     + -b/(tm)*(-wi))
    c4 = @constraint(pm.model, q_to == -b*w_to +  b/(tm)*wr     + -g/(tm)*(-wi))
    return c1, c2, c3, c4
end
"""
Converter reactor constraints
```
p_pr_fr ==  g/(tm^2)*w_fr + -g/(tm)*wr + -b/(tm)*wi)
q_pr_fr == -b/(tm^2)*w_fr +  b/(tm)*wr + -g/(tm)*wi)
p_pr_to ==  g*w_to + -g/(tm)*wr     + -b/(tm)*(-wi))
q_pr_to == -b*w_to +  b/(tm)*wr     + -g/(tm)*(-wi))
```
"""
function constraint_conv_reactor(pm::AbstractWRModels, n::Int, cnd::Int, i::Int, rc, xc, reactor)
    ppr_fr = PowerModels.var(pm, n, cnd, :pconv_pr_fr, i)
    qpr_fr = PowerModels.var(pm, n, cnd, :qconv_pr_fr, i)
    ppr_to = -PowerModels.var(pm, n, cnd, :pconv_ac, i)
    qpr_to = -PowerModels.var(pm, n, cnd, :qconv_ac, i)

    wf = PowerModels.var(pm, n, cnd, :wf_ac, i)
    wc = PowerModels.var(pm, n, cnd, :wc_ac, i)
    wrc = PowerModels.var(pm, n, cnd, :wrc_ac, i)
    wic = PowerModels.var(pm, n, cnd, :wic_ac, i)

    zc  = rc  + im*xc

    if reactor
        yc = 1/(zc)
        gc = real(yc)
        bc = imag(yc)
        constraint_voltage_product_converter(pm, wrc, wic, wf, wc)
        c1, c2, c3, c4 = ac_power_flow_constraints_w(pm, gc, bc, wf, wc, wrc, wic, ppr_fr, ppr_to, qpr_fr, qpr_to, 1)
        PowerModels.con(pm, n, cnd, :conv_pr_p)[i] = c1
        PowerModels.con(pm, n, cnd, :conv_pr_q)[i] = c2
    else
        pcon, qcon = constraint_lossless_section(pm, wf, wc, wrc, wic, ppr_fr, ppr_to, qpr_fr, qpr_to)
        PowerModels.con(pm, n, cnd, :conv_pr_p)[i] = pcon
        PowerModels.con(pm, n, cnd, :conv_pr_q)[i] = qcon
    end
end
"""
Converter filter constraints

```
p_pr_fr + p_tf_to == 0
q_pr_fr + q_tf_to + -bv*filter*wf == 0
```
"""
function constraint_conv_filter(pm::AbstractWModels, n::Int, cnd::Int, i::Int, bv, filter)
    ppr_fr = PowerModels.var(pm, n, cnd, :pconv_pr_fr, i)
    qpr_fr = PowerModels.var(pm, n, cnd, :qconv_pr_fr, i)
    ptf_to = PowerModels.var(pm, n, cnd, :pconv_tf_to, i)
    qtf_to = PowerModels.var(pm, n, cnd, :qconv_tf_to, i)

    wf = PowerModels.var(pm, n, cnd, :wf_ac, i)

    PowerModels.con(pm, n, cnd, :conv_kcl_p)[i] = @constraint(pm.model, ppr_fr + ptf_to == 0 )
    PowerModels.con(pm, n, cnd, :conv_kcl_q)[i] = @constraint(pm.model, qpr_fr + qtf_to + -bv*filter*wf == 0)
end

function constraint_lossless_section(pm::AbstractWModels, w_fr, w_to, wr, wi, p_fr, p_to, q_fr, q_to)
    @constraint(pm.model, w_fr ==  w_to)
    @constraint(pm.model, wr   ==  w_fr)
    @constraint(pm.model, wi   ==  0)

    pcon = @constraint(pm.model, p_fr + p_to == 0)
    qcon = @constraint(pm.model, q_fr + q_to == 0)
    return pcon, qcon
end

function add_converter_voltage_setpoint(sol, pm::AbstractWModels)
    PowerModels.add_setpoint!(sol, pm, "convdc", "vmconv", :wc_ac; scale = (x,item) -> sqrt(x))
    PowerModels.add_setpoint!(sol, pm, "convdc", "vmfilt", :wf_ac; scale = (x,item) -> sqrt(x))
end
"""
LCC firing angle constraints
```
qconv_ac >= Q1 + (pconv_ac-P1) * (Q2-Q1)/(P2-P1)

P1 = cos(0) * Srated
Q1 = sin(0) * Srated
P2 = cos(pi) * Srated
Q2 = sin(pi) * Srated
```
"""
function constraint_conv_firing_angle(pm::AbstractWModels, n::Int, cnd::Int, i::Int, S, P1, Q1, P2, Q2)
    pc = PowerModels.var(pm, n, cnd, :pconv_ac, i)
    qc = PowerModels.var(pm, n, cnd, :qconv_ac, i)
    coeff = (Q2-Q1)/(P2-P1)

    PowerModels.con(pm, n, cnd, :conv_socphi)[i] = @constraint(pm.model, qc >= Q1 + (pc-P1) * coeff )
end


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

name = "SARSOP"
version = v"0.96.0"

# Collection of sources required to complete build
sources = [
    GitSource("https://github.com/JuliaPOMDP/sarsop.git", "903fa55df637d127c1068700377995ee395475a6")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
cd sarsop/src
make -j${nprocs} CC=${CC} CXX=${CXX} LINKER=${CXX} DESTDIR=${prefix}
make install DESTDIR=${prefix}
if [[ "${target}" == *mingw* ]]; then
    for file in ${prefix}/bin/*;
    do
        mv $file $file.exe
    done
fi
install_license ${WORKSPACE}/srcdir/sarsop/license/*
exit
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Linux(:i686, libc=:glibc),
    Linux(:x86_64, libc=:glibc),
    MacOS(:x86_64),
    Windows(:i686),
    Windows(:x86_64)
]
platforms = expand_cxxstring_abis(platforms)

# The products that we will ensure are always built
products = [
    ExecutableProduct("pomdpsim", :pomdpsim),
    ExecutableProduct("pomdpconvert", :pomdpconvert),
    ExecutableProduct("pomdpeval", :pomdpeval),
    ExecutableProduct("pomdpsol", :pomdpsol),
    ExecutableProduct("polgraph", :polgraph)
]

# Dependencies that must be installed before this package can be built
dependencies = Dependency[
]

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


text:
__precompile__()

module MPPI

using mjWrap
using LearningStrategies
import LearningStrategies: setup!, update!, hook
using ExpFunctions
using Common

export MPPIModel, MPPIStrategy, setup!, update!, hook
export MPPIRollout
export MPCRollout

struct MPPIModel
   theta::Matrix{Float64} # parameters to improve
   mjsys::mjWrap.mjSet    # environment
   s0::Vector{Float64}    # initial / current state
   terminalvalue::Function

   Σ::Matrix{Float64}
   λ::Float64

   MPC::Bool

   # scratch stuff
   samples::TrajSamples
   traj::TrajSamples

   trace::Dict{Symbol,Vector{Float64}}

   niter::Int

   function MPPIModel(sigma, lambda, mjsys, mpcT, T, numT;
                      theta=zeros(mjsys.nu, mpcT),
                      s0=zeros(mjsys.nq + mjsys.nv)*NaN, # set NaN init to signal Rollout functions
                      valuefunction=(x)->zeros(numT),
                      niter=T)
      ns = mjsys.ns

      if mpcT == T # not MPC mode
         samples = mjw.allocateTrajSamples(mjsys, T, numT, ns)
         traj    = mjw.allocateTrajSamples(mjsys, T, 1, ns)
         mpc     = false
         warn("MPPI in FULL TRAJECTORY Mode")
      else
         samples = mjw.allocateTrajSamples(mjsys, mpcT, numT, ns)
         traj    = mjw.allocateTrajSamples(mjsys, T, 1, ns)
         mpc     = true
         warn("MPPI in Model Predictive Mode")
      end

      return new(theta, mjsys, s0, valuefunction, sigma, lambda, mpc,
                 samples, traj,
                 Dict(:stocR => Vector{Float64}(),
                      :meanR => Vector{Float64}(),
                      :evalscore => Vector{Float64}() ), niter)
   end
end

struct MPPIStrategy <: LearningStrategy
   costs::Vector{Float64}
   weights::Vector{Float64}

   function MPPIStrategy(model)
      dtype   = Float64
      numT    = model.samples.numT

      costs   = Vector{dtype}(numT)
      weights = Vector{dtype}(numT)

      new(costs, weights)
   end
end
Base.show(io::IO, s::MPPIStrategy) = print(io, "MPPI")

function apply(f::Function, costs::Vector, states::Array{Float64,3})
   for k=1:size(states,3)
      costs[k] -= f(@view states[:,end,k])
   end
   #costs .-= f(states[:,end,:]) # block calculate
end

function update!(model::MPPIModel, s::MPPIStrategy, iter, null)

   BLAS.set_num_threads(Threads.nthreads())

   mjsys   = model.mjsys
   samples = model.samples
   T       = samples.T
   numT    = samples.numT

   s.costs[:] = -1.0.*sum(samples.reward, 1) # rollouts done in different strat

   baseline = minimum(s.costs)

   push!(model.trace[:stocR], mean(s.costs))
   push!(model.trace[:meanR], baseline)

   @. s.costs = exp((-(s.costs-baseline)/model.λ))

   η = sum(s.costs)
   s.weights .= s.costs/η

   #SIGMA = true
   #if SIGMA
   #   #println("before:\n",model.Σ)
   #   sig = reshape(reshape(broadcast(-, samples.ctrl, model.theta),
   #                         mjsys.nu*T, numT)*s.weights,
   #                 mjsys.nu, T)
   #   model.Σ .= 0.0
   #   for t=1:T
   #      BLAS.gemm!('N', 'T', 1.0, sig[:,t], sig[:,t], 1.0, model.Σ)
   #   end
   #   model.Σ ./= T
   #   for u=1:mjsys.nu
   #      model.Σ[u,u] .= max(1e-1, model.Σ[u,u])
   #   end
   #   if mod(iter, 100) == 0
   #      println("Σ:\n",mean(model.Σ))
   #   end
   #end

   model.theta .= reshape(reshape(samples.ctrl,
                                  mjsys.nu*T, numT)*s.weights,
                          mjsys.nu, T)

end

############################################################### ROLLOUTS

struct MPPIRollout <: LearningStrategy
   f::FunctionSet
   c::Vector{Float64}
   o::Vector{Float64}
   function MPPIRollout(f)
      new(f, Array{Float64,1}(), Array{Float64,1}())
   end
end
Base.show(io::IO, s::MPPIRollout) = print(io, "MPPI rollout")

function setup!(s::MPPIRollout, model::MPPIModel)
   samples = model.samples

   push!(s.c, zeros(model.mjsys.nu)...)
   push!(s.o, zeros(model.mjsys.ns)...)

   # traj and samples both start from same place; send inital state to MPPIModel
   if isnan(model.s0[1])
      info("Setting initial state from experiment.")
      s.f.initstate!(model.mjsys, model.traj)
      model.s0 .= model.traj.state[:,1,1]
   else
      model.traj.state[:,1,1] .= model.s0
   end
   samples.state[:,1,:] .= model.s0
end

function update!(model::MPPIModel, s::MPPIRollout, iter, null)
   samples = model.samples
   T       = samples.T
   numT    = samples.numT
   mjsys   = model.mjsys

   noise = reshape( model.Σ * randn(mjsys.nu, T*numT), mjsys.nu, T, numT)
   for k=1:numT
      samples.ctrl[:,1,k] = noise[:,1,k]
      for i=2:T-1
         @. samples.ctrl[:,i,k] = 0.25*noise[:,i-1,k]+0.5*noise[:,i,k]+0.25*noise[:,i+1,k]
      end
      samples.ctrl[:,T,k] = noise[:,T,k]
   end

   samples.ctrl[:,:,1] .= model.theta

   broadcast!(+, samples.ctrl, samples.ctrl, model.theta)
   samples.ctrl[:,:,1] .= model.theta # have at least one thing to do a full theta rollout

   #for u=1:mjsys.nu
   #   samples.ctrl[u,:,:] .= clamp.(samples.ctrl[u,:,:],
   #                                 mjsys.m.actuator_ctrlrange[1,u],
   #                                 mjsys.m.actuator_ctrlrange[2,u])
   #end

   mjWrap.rollout(mjsys, samples, (x...)->nothing, s.f.observe!, s.f.reward)
end

function hook(s::MPPIRollout, model::MPPIModel, iter) # just rollout and store in traj
   if model.MPC # get first controls, mj.step, store in model.traj
      mjsys   = model.mjsys
      traj    = model.traj
      samples = model.samples

      # MPC: apply first controls to system; shift the rist around
      s.c .= model.theta[:,1]
      model.theta[:,1:end-1] .= model.theta[:,2:end]
      model.theta[:,end] .= 0.0

      mjWrap.reset(mjsys, 1, model.s0, s.c, s.o)

      mjWrap.step!(mjsys, 1, s.c, model.s0, s.o) # get next state, next observation
      s.f.observe!(mjsys, 1, model.s0, s.o) # arbitrary observation vector manipulation

      # advance our main trajectory
      traj.state[:,iter,1] .= model.s0
      traj.ctrl[:,iter,1]  .= s.c
      traj.obs[:,iter,1]   .= s.o

      # update the mpc bundle of trajectories to start 'here'
      samples.state[:,1,:] .= model.s0
   elseif mod1(iter, model.niter) == model.niter
      #info("Running rollout for main trajectory.")
      mjsys   = model.mjsys
      traj    = model.traj

      traj.ctrl[:,:,1] = model.theta
      mjWrap.rollout(mjsys, traj, (x...)->nothing, s.f.observe!, s.f.reward)
   end
end

end


text:
module TabularDisplay

using Formatting

export displaytable

"""
    displaytable([io::IO], v::AbstractVector{T}; kwargs...)

Display a vector of things in a compact tabular format.
The `io` argument is defaulted to `stdout`.

# Keyword arguments
- padding: minimum number of spaces between column data. Default value is 2.
- index: prepend cell values with indices. Default value is `false`.
- indexsep: string that separate index and cell values.  Default value is `:`.
- align: either `:left` or `:right`.  Default value is `:left`.
- orientation: either `:column` or `:row`.  Default value is `:column`.
- formatter: custom formatter that takes a value and returns a string.  Default value is `string` function.
- displaywidth: custom display width.  Default value is 0, for which the system will use the terminal's size.

# Examples
```
displaytable([string("randomstr", i) for i in 1:56])

using Formatting
foo = generate_formatter("%7.5f")
displaytable(rand(100); padding=5, align=:right, formatter=foo, index=true, indexsep=" -> ")
```
"""
function displaytable end

function displaytable(v::AbstractVector{T}; kwargs...) where T
    return displaytable(stdout, v; kwargs...)
end

function displaytable(io::IO, v::Vector{T};
            padding = 2,
            index = false,
            indexsep = ":",
            align = :left,
            orientation = :column,
            formatter = string,
            displaywidth = 0,
        ) where T

    # convert vector to strings and prepend with index if needed
    ar = formatter.(v)
    if index
        ar = prependindex(ar, indexsep)
    end

    # Calculate number of rows and columns required
    rowwidth = displaywidth > 0 ? displaywidth : displaysize(io)[2]
    #rowwidth -= 1    # at least one space from right margin
    cellwidth = maximum(length.(ar)) + padding
    rowwidth < cellwidth && error("The value for displaywidth is too small. Please specify a value larger than $cellwidth.")
    columns = div(rowwidth, cellwidth)
    rows = ceil(Int, length(ar) / columns)

    # Generate cell and row format string
    cellfmt = string("{:", align == :left ? "<" : ">", string(cellwidth), "s}")
    rowfmt = repeat(cellfmt, columns)

    # Copy the array and append it with blanks for remaining cells
    # This is done to make the loop below easier
    emptyslots = rows * columns - length(ar)
    ar = vcat(ar, ["" for _ in 1:emptyslots])

    # print each row with stride
    for i in 1:rows
        i > 1 && println(io)
        if orientation == :column
            printfmt(io, rowfmt, ar[i:rows:length(ar)]...)
        else
            startidx = (i - 1) * columns + 1
            endidx = startidx + columns - 1
            printfmt(io, rowfmt, ar[startidx:endidx]...)
        end
    end
end

function prependindex(v::AbstractVector{T}, sep::String) where T
    [string(i, sep, v[i]) for i in 1:length(v)]
end

end # module


text:
module MultiFEOperatorsTests

using Test
using Gridap
import Gridap: ∇

# Define manufactured functions
u1fun(x) = x[1] + x[2]
u2fun(x) = x[1] - x[2]

u1fun_grad(x) = VectorValue(1.0,1.0)
u2fun_grad(x) = VectorValue(1.0,-1.0)

∇(::typeof(u1fun)) = u1fun_grad
∇(::typeof(u2fun)) = u2fun_grad

b1fun(x) = u2fun(x)
b2fun(x) = 0.0

g1fun(x) = 1.0

# Construct the discrete model
model = CartesianDiscreteModel(domain=(0.0,1.0,0.0,1.0), partition=(4,4))

# Construct the FEspace.
# We allow each one of the fields to have different boundary conditions
order = 1
diritags = [1,2,3,4,5,6,7]
fespace1 = H1ConformingFESpace(Float64,model,order,diritags)

diritag = "boundary"
fespace2 = H1ConformingFESpace(Float64,model,order,diritag)

# Define test and trial
V1 = TestFESpace(fespace1)
V2 = TestFESpace(fespace2)
V = [V1, V2]

U1 = TrialFESpace(fespace1,u1fun)
U2 = TrialFESpace(fespace2,u2fun)
U = [U1, U2]

# Define integration mesh and quadrature for volume
trian = Triangulation(model)
quad = CellQuadrature(trian,degree=2)

# Setup integration on Neumann boundary
neumanntags = [8,]
btrian = BoundaryTriangulation(model,neumanntags)
bquad = CellQuadrature(btrian,degree=2)

# Terms in the volume
a(v,u) = inner(∇(v[1]),∇(u[1])) + inner(v[1],u[2]) + inner(∇(v[2]),∇(u[2]))
b(v) = inner(v[1],b1fun) + inner(v[2],b2fun)
t_Ωa = LinearFETerm(a,trian,quad)
t_Ωb = FESource(b,trian,quad)

# Terms on Neumann boundary
# Note that the Neumann BC only applies on the first field
g(v) = inner(v[1],g1fun)
t_Γ = FESource(g,btrian,bquad)

# Define Assembler
assem = SparseMatrixAssembler(V,U)

# Define the FEOperator
op = LinearFEOperator(V,U,assem,t_Ωa,t_Ωb,t_Γ)

# Define the FESolver
ls = LUSolver()
solver = LinearFESolver(ls)

# Solve!
uh = solve(solver,op)

# Define exact solution and error
u1 = CellField(trian,u1fun)
e1 = u1 - uh[1]

u2 = CellField(trian,u2fun)
e2 = u2 - uh[2]

# Define norms to measure the error
l2(u) = inner(u,u)
h1(u) = inner(∇(u),∇(u)) + l2(u)

# Compute errors
e1l2 = sqrt(sum( integrate(l2(e1),trian,quad) ))
e1h1 = sqrt(sum( integrate(h1(e1),trian,quad) ))

e2l2 = sqrt(sum( integrate(l2(e2),trian,quad) ))
e2h1 = sqrt(sum( integrate(h1(e2),trian,quad) ))

@test e1l2 < 1.e-8
@test e1h1 < 1.e-8

@test e2l2 < 1.e-8
@test e2h1 < 1.e-8

# Further tests

t_Ω = AffineFETerm(a,b,trian,quad)
op = LinearFEOperator(a,b,V,U,assem,trian,quad)
op = LinearFEOperator(V,U,t_Ω,t_Γ)
op = LinearFEOperator(V,U,t_Ω)

end


text:
struct SequencingModel{U<:ObjectiveFunction, V<:Int}
    objective               ::U
    length                  ::V
    seqStart                ::Union{Nothing, V}
    seqEnd                  ::Union{Nothing, V}
    precedenceConstraints   ::Bool #true if there are precedence constraints
    precedenceP2F           ::Union{Nothing, Dict{Int, Vector{V}}}#map of priors to followers
    precedenceF2P           ::Union{Nothing, Dict{Int, Vector{V}}} #map of followers to priors
    precedenceNumbers       ::Union{Nothing, Dict{V, V}}#how many things must come before each node
end

Base.length(model::SequencingModel) = model.length

"""
    hasStart(model::SequencingModel)

Check if the sequence has a predefined start

# Arguments
- `model::SequencingModel`: The problem model
"""
function hasStart(model::SequencingModel)
    return !isnothing(model.seqStart)
end

"""
    hasEnd(model::SequencingModel)

Check if the sequence has a predefined end

# Arguments
- `model::SequencingModel`: The problem model
"""
function hasEnd(model::SequencingModel)
    return !isnothing(model.seqEnd)
end

"""
    isBetterSolutionValue(model::SequencingModel, newValue::T, bestKnownValue::T)where{T<:Real}

Check if a solution value is better than the best known value

# Arguments
- `model::SequencingModel`: The problem model
"""
function isBetterSolutionValue(model::SequencingModel, newValue::T, bestKnownValue::T)where{T<:Real}
    if model.objective.type == minimization
        return newValue<bestKnownValue
    else
        return newValue>bestKnownValue
    end
end

"""
    checkValidity(model::SequencingModel, seq::Vector{T})where{T<:Int}

Check if a solution is feasible

# Arguments
- `model::SequencingModel`: The problem model
- `seq::Vector{T}`: The sequence to check
"""
function checkValidity(model::SequencingModel, seq::Vector{T})where{T<:Int}
    #precedence
    if checkPrecedence(model, seq)
        #alldiff
        used = falses(length(model))
        for e in seq
            used[e] = true
        end
        if used == trues(length(model))
            return true
        else
            return false
        end
    else
        return false
    end
end

"""
    checkPrecedence(model::SequencingModel, seq::Vector{T})where{T<:Int}

Check if a solution satisfies all the precedence constraints

# Arguments
- `model::SequencingModel`: The problem model
- `seq::Vector{T}`: The sequence to check
"""
function checkPrecedence(model::SequencingModel, seq::Vector{T})where{T<:Int}
    used = falses(length(seq))
    for e in seq
        used[e] = true
        if haskey(model.precedenceF2P,e)
            for val in values(model.precedenceF2P[e])
                if !used[val]
                    return false
                end
            end
        end
    end
    return true
end


text:
using BinaryBuilder

name = "GnuTLS"
version = v"3.6.13"

# Collection of sources required to build GnuTLS
sources = [
    ArchiveSource("https://www.gnupg.org/ftp/gcrypt/gnutls/v3.6/gnutls-3.6.13.tar.xz",
                  "32041df447d9f4644570cf573c9f60358e865637d69b7e59d1159b7240b52f38"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/gnutls-*/

# Grumble-grumble apple grumble-grumble broken linkers...
#if [[ ${target} == *-apple-* ]]; then
#    export AR=/opt/${target}/bin/ar
#fi

GMP_CFLAGS="-I${prefix}/include" ./configure --prefix=${prefix} --host=${target} \
    --with-included-libtasn1 \
    --with-included-unistring \
    --without-p11-kit

make -j${nproc}
make install
"""

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

# Wimp out of doing FreeBSD since we don't have that for some targets
platforms = [p for p in platforms if !(typeof(p) <: FreeBSD)]

# Disable windows because O_NONBLOCK isn't defined
platforms = [p for p in platforms if !(typeof(p) <: Windows)]

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

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency("Zlib_jll"),
    Dependency("GMP_jll", v"6.1.2"),
    Dependency("Nettle_jll"),
]

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

max_stars_repo_name: mortenpi/Yggdrasil

max_issues_repo_name: mortenpi/Yggdrasil

max_forks_repo_name: mortenpi/Yggdrasil


text:
"Diagonally implicit Runge-Kutta integrator."
immutable IntegratorDIRK{T} <: Integrator{T}
    equation::ODE{T}
    tableau::TableauDIRK{T}
    Δt::T

    x::Array{T,1}
    X::Array{T,2}
    Y::Array{T,2}
    F::Array{T,2}

    function IntegratorDIRK(equation, tableau, Δt)
        D = equation.d
        S = tableau.s
        new(equation, tableau, Δt, zeros(T,D), zeros(T,D,S), zeros(T,D,S), zeros(T,D,S))
    end
end

function IntegratorDIRK(equation::Equation, tableau::TableauDIRK, Δt)
    T = eltype(equation.q₀)
    IntegratorDIRK{T}(equation, tableau, Δt)
end

"Integrate ODE with diagonally implicit Runge-Kutta integrator."
function integrate!(int::IntegratorDIRK, s::SolutionODE)
    # TODO
end

"Integrate partitioned ODE with diagonally implicit Runge-Kutta integrator."
function integrate!(int::IntegratorDIRK, s::SolutionPODE)
    # TODO
end


text:
using Vlasiator, PyPlot, JLD2
using Vlasiator: RE # Earth radius [m]

## Plot velocity space distributions
xloc = range(7.0, 10.0, step=2.35)
yloc = range(-0.0, 5.0, step=2.35)
zloc = 0.0

CIs = CartesianIndices((1:length(xloc), 1:length(yloc)))

fig, axs = plt.subplots(length(xloc), length(yloc), sharex=true, sharey=true)

for i in CIs
   loc = RE .* [xloc[i[1]], yloc[i[2]], zloc]
   vdfslice(meta, loc, axs[i]; verbose=false)
end

for a in axs[end,:] a.set_xlabel("vx [km/s]") end
for a in axs[:  ,1] a.set_ylabel("vy [km/s]") end

## Subtract data of cells with VDF

file = "test/data/bulk.1d.vlsv"

meta = load(file)

cells = getcellwithvdf(meta)

ρ = zeros(Float64, size(cells))
v = zeros(Float64, 3, size(cells)...)
p = zeros(Float64, size(cells))

vsize = meta.meshes["proton"].vblock_size .* meta.meshes["proton"].vblocks
vcellids = [Int64[] for _ in cells]
vcellf = [Float64[] for _ in cells]

for (i, cell) in enumerate(cells)
   ρ[i] = readvariable(meta, "proton/vg_rho", cell)[1]
   v[i] = readvariable(meta, "proton/vg_v", cell)[1]
   p[i] = readvariable(meta, "vg_pressure", cell)[1]

   vcellids[i], vcellf[i] = readvcells(meta, cell)
end

jldopen("example.jld2", "a+") do file
   file["file1/density"] = ρ
   file["file1/velocity"] = v
   file["file1/pressure"] = p
   file["file1/VDF"] = vcellf
end


text:
"""
    Code for a GaussianMixture mask shape for use with CCF
Author: Eric Ford
Created: September 2020
"""

import SpecialFunctions: erf
using StaticArrays

"""   GaussianMixtureCCFMask
A mixture of truncated GaussianMixtures mask with three parameters per mixture component: weight,  standard deviation and offset (in m/s)
and one parameter for the velocity width at which the whole thing is truncated.
Mask weights are stored separately in a line list.
TODO: Test before using
"""
struct GaussianMixtureMixtureCCFMask{NMix::Integer} <: AbstractCCFMaskShape
    weight::SVector{NMix,Float64}
    σ_sqrt2::SVector{NMix,Float64}
    v_offset::SVector{NMix,Float64}
    truncation_Δv::Float64
    normalization::SVector{NMix,Float64}
    cdf_normalization::SVector{NMix,Float64}

    """ GaussianMixtureCCFMask( σ ; half_truncation_width_in_σ=2 ) """
    function GaussianMixtureCCFMask(weight::AbstractVector{Float64}, σ::AbstractVector{Float64}, truncation_Δv::Real; v_offset::AbstractVector{Float64}=zeros(length(weight)) )
        n = legnth(weight)
        @assert 1 <= n <= 10   # Arbitrary upper limit
        @assert length(σ) == n
        @assert length(v_offset) == n
        @assert all(0 .<= weight .<= one(eltype(weight)))   #
        @assert all(0 .< σ .<= 300000)   # 300km/s is arbitrary choice for an upper limit
        @assert all(-1000 .<= v_offset .<= 1000)   # 1km/s is arbitrary choice for an upper limit
        @assert 0 < truncation_Δv <= 300000

        pdf_norm = zeros(n)
        cdf_norm = zeros(n)
        for i in 1:n
            integral = ( erf((truncation_Δv-v_offset[i])/(sqrt(2.0)*σ[i])) - erf((-truncation_Δv-v_offset[i])/(sqrt(2.0*σ[i]))) )
            pdf_norm[i] = weight[i]/(sqrt(2π)*σ[i]*integral)
            cdf_norm[i] = weight[i]/integral
        end
        new(weight, σ.*sqrt(2.0),v_offset, truncation_Δv,pdf_norm,cdf_norm)
    end

end

""" `GaussianMixtureCCFMask( inst ; scale_factor )`
For now just makes a single Gaussian, so for testing purposes only.
"""
function GaussianMixtureCCFMask(inst::AbstractInstrument; σ_scale_factor::Real = 1, truncation_Δv::Real = default_gaussian_mixture_ccf_truncation_Δv )
    σ = σ_scale_factor * default_ccf_mask_v_width(inst)
    GaussianMixtureCCFMask([1.0], [σ],truncation_Δv,v_offset=[0.])
end

λ_min(m::GaussianMixtureCCFMask,λ::Real) = λ/calc_doppler_factor(m.truncation_Δv)
λ_max(m::GaussianMixtureCCFMask,λ::Real) = λ*calc_doppler_factor(m.truncation_Δv)

function cdf(m::GaussianMixtureCCFMask, v_lo::Real,v_hi::Real, i::Integer)
    m.cdf_normalization[i]*( erf((v_hi-v_offset[i])/σ_sqrt2[i]) - erf((-v_lo-v_offset[i])/σ_sqrt2[i]) )
end

function pdf(m::GaussianMixtureCCFMask, Δv::Real, i::Integer)
    m.normalization[i]*exp(-(Δv-m.v_offset[i]/m.σ_sqrt2[i])^2)
end

function integrate(m::GaussianMixtureCCFMask, v_lo::Real,v_hi::Real)
    @assert -m.truncation_Δv <= v_lo < v_hi <= m.truncation_Δv
    #quadgk(m, v_lo, v_hi)[1]
    mapreduce(i->cdf(m,v_lo,v_hi,i),+,1:length(m.weights))
end

""" Functor for returning PSF for Δv <= half_width.  """
function (m::GaussianMixtureCCFMask)(Δv::Real)
    if abs2(Δv) > abs2(m.truncation_Δv)
        return zero(Δv)
    else
        return mapreduce(i->pdf(m,Δv,i),+,1:length(m.weights))
    end
end

function mask_with_increased_fwhm(m::GaussianMixtureCCFMask, Δfwhm::Real )
    fwhm_orig = 2*sqrt(log(2)) .* m.σ_sqrt2
    fwhm_new = sqrt.(fwhm_orig.^2 .+ Δfwhm^2)
    σ_new = 0.5 .* fwhm_new./sqrt(2*log(2))
    return GaussianMixtureCCFMask(m.weight, σ_new, m.truncation_Δv, v_offset=m.v_offset )
end

max_stars_repo_name: RvSpectML/EchelleCCFs

max_issues_repo_name: RvSpectML/EchelleCCFs

max_forks_repo_name: RvSpectML/EchelleCCFs


text:
export iauAe2hd
"""
Horizon to equatorial coordinates:  transform azimuth and altitude
to hour angle and declination.

Given:
   az       double       azimuth
   el       double       altitude (informally, elevation)
   phi      double       site latitude

Returned:
   ha       double       hour angle (local)
   dec      double       declination

Notes:
   1. All the arguments are angles in radians.

   2. The sign convention for azimuth is north zero, east +pi/2.

   3. HA is returned in the range +/-pi.  Declination is returned in
   the range +/-pi/2.

   4. The latitude phi is pi/2 minus the angle between the Earth's
   rotation axis and the adopted zenith.  In many applications it
   will be sufficient to use the published geodetic latitude of the
   site.  In very precise (sub-arcsecond) applications, phi can be
   corrected for polar motion.

   5. The azimuth az must be with respect to the rotational north pole,
   as opposed to the ITRS pole, and an azimuth with respect to north
   on a map of the Earth's surface will need to be adjusted for
   polar motion if sub-arcsecond accuracy is required.

   6. Should the user wish to work with respect to the astronomical
   zenith rather than the geodetic zenith, phi will need to be
   adjusted for deflection of the vertical (often tens of
   arcseconds), and the zero point of ha will also be affected.

   7. The transformation is the same as Ve = Ry(phi-pi/2)*Rz(pi)*Vh,
   where Ve and Vh are lefthanded unit vectors in the (ha,dec) and
   (az,el) systems respectively and Rz and Ry are rotations about
   first the z-axis and then the y-axis.  (n.b. Rz(pi) simply
   reverses the signs of the x and y components.)  For efficiency,
   the algorithm is written out rather than calling other utility
   functions.  For applications that require even greater
   efficiency, additional savings are possible if constant terms
   such as functions of latitude are computed once and for all.

   8. Again for efficiency, no range checking of arguments is carried
   out.

Last revision:   2017 September 12

SOFA release 2018-01-30

Copyright (C) 2018 IAU SOFA Board.  See notes at end.
"""
function iauAe2hd(az::Real, el::Real, phi::Real)
   ha  = zeros(Float64, 1)
   dec = zeros(Float64, 1)

   ccall((:iauAe2hd, libsofa_c), Cvoid,
         (Cdouble, Cdouble, Cdouble, Ptr{Cdouble}, Ptr{Cdouble}),
         convert(Float64, az),
         convert(Float64, el),
         convert(Float64, phi),
         ha, dec)

   return ha[1], dec[1]
end


text:
# This is used in VSCode and others

# If the user displays a AstroImageMat of colors (e.g. one created with imview)
# fal through and display the data as an image
Base.show(io::IO, mime::MIME"image/png", img::AstroImageMat{T}; kwargs...) where {T<:Colorant} =
    show(io, mime, parent(img), kwargs...)

# Otherwise, call imview with the default settings.
Base.show(io::IO, mime::MIME"image/png", img::AstroImageMat{T}; kwargs...) where {T<:Union{Number,Missing}} =
    show(io, mime, imview(img), kwargs...)

# Deprecated
# Lazily reinterpret the AstroImageMat as a Matrix{Color}, upon request.
# By itself, Images.colorview works fine on AstroImages. But
# AstroImages are not normalized to be between [0,1]. So we override
# colorview to first normalize the data using scaleminmax
@deprecate render(img::AstroImageMat) imview(img, clims=extrema, cmap=nothing)


text:
#= This is a code to conduct the bifurcation analysis with all n fixed, one to 1
and two to 0 =#

using SymEngine
using LinearAlgebra
using PolynomialRoots
using Optim
using TensorOperations
using FileIO
include("material_parameters.jl")
include("bifurcation_matrix.jl")
include("element_integration.jl")

#= We start by simplifying and numericising the bifurcation matrix, and applying
the rules constraining the normal vector=#
b_m_s = [symbols("b_ms_s$i$j$t") for i in 1:2*ndim, j in 1:2*ndim, t in 1:n_timesteps+1]
n_rules = [1, 0, 0]
for t in 1:n_timesteps+1
    b_m_s[:, :, t] = coefficient_matrix
    for i in 1:ndim
        for j in 1:ndim
            for k in 1:ndim
                for l in 1:ndim
                    for bm_index1 in 1:2*ndim
                        for bm_index2 in 1:2ndim
                            b_m_s[bm_index1, bm_index2, t] = subs(b_m_s[bm_index1, bm_index2, t], n[i]=>n_rules[i])
                            b_m_s[bm_index1, bm_index2, t] = subs(b_m_s[bm_index1, bm_index2, t], Eepsym[i, j, k, l]=>Eep_res[i, j, k, l, t], Fepsym[i, j, k, l]=>Fep_res[i, j, k, l, t], Kepsym[i, j, k, l]=>Kep_res[i, j, k, l, t], Mepsym[i, j, k, l]=>Mep_res[i, j, k, l, t])
                            b_m_s[bm_index1, bm_index2, t] = subs(b_m_s[bm_index1, bm_index2, t], Irot=>Ires[t])
                        end
                    end
                end
            end
        end
    end
end

# Define a general Levi-Civita function
function generalised_LC(index_vector)::Int64
    value = 1
    for index1 in 1:(length(index_vector)-1)
        for index2 in (index1+1):length(index_vector)
            value = value*sign(index_vector[index2] - index_vector[index1])
        end
    end
    return value
end

#= Get the characteristic polynomial as a function of the time, orientation,
wavelength and value of sq =#
num_mat = [symbols("nm_$t") for t in 1:n_timesteps+1]
fully_subbed = [symbols("fs$i$j") for i in 1:2*ndim, j in 1:2*ndim]
L_mat = [symbols("lm$i$j") for i in 1:2*ndim, j in 1:2*ndim]
U_mat = [symbols("um$i$j") for i in 1:2*ndim, j in 1:2*ndim]
order_6_LC = [generalised_LC([i,j,k,l,m,o]) for i in 1:2*ndim, j in 1:2*ndim, k in 1:2*ndim, l in 1:2*ndim, m in 1:2*ndim, o in 1:2*ndim]
#n_m_t = symbols(:n_m_t)
row_1 = [symbols("rw_$i") for i in 1:2*ndim]
row_2 = [symbols("rw_$i") for i in 1:2*ndim]
row_3 = [symbols("rw_$i") for i in 1:2*ndim]
row_4 = [symbols("rw_$i") for i in 1:2*ndim]
row_5 = [symbols("rw_$i") for i in 1:2*ndim]
row_6 = [symbols("rw_$i") for i in 1:2*ndim]
r1_product = [symbols("r1") for j in 1:2*ndim, k in 1:2*ndim, l in 1:2*ndim, m in 1:2*ndim, o in 1:2*ndim]
r2_product = [symbols("r2") for k in 1:2*ndim, l in 1:2*ndim, m in 1:2*ndim, o in 1:2*ndim]
r3_product = [symbols("r3") for l in 1:2*ndim, m in 1:2*ndim, o in 1:2*ndim]
r4_product = [symbols("r4") for m in 1:2*ndim, o in 1:2*ndim]
r5_product = [symbols("r5") for o in 1:2*ndim]
polycoeff = [symbols("pc$i") for i in 1:2*ndim+1]
function polyΛsq(t)
    fully_subbed = b_m_s[:, :, t]
    #= We use the Levi-Civita tensor formulation of the Leibniz formula to get the
    determinant, but we help it out first by creating some auxiliaries =#
    row_1 = fully_subbed[1, :]
    row_2 = fully_subbed[2, :]
    row_3 = fully_subbed[3, :]
    row_4 = fully_subbed[4, :]
    row_5 = fully_subbed[5, :]
    row_6 = fully_subbed[6, :]
    @tensor begin
        r1_product[j, k, l, m, o] = order_6_LC[i, j , k, l, m, o]*row_1[i]
        r2_product[k, l, m, o] = r1_product[j, k, l, m, o]*row_2[j]
        r3_product[l, m, o] = r2_product[k, l, m, o]*row_3[k]
        r4_product[m, o] = r3_product[l, m, o]*row_4[l]
        r5_product[o] = r4_product[m, o]*row_5[m]
    end
    #= Wrap it up with a loop because TensorOperations doesn't want to give me a
    scalar =#
    n_m_t = symbols(:n_m_t)
    for o = 1:2*ndim
        if o == 1
            n_m_t = r5_product[o]*row_6[o]
        else
            n_m_t = n_m_t + r5_product[o]*row_6[o]
        end
    end
    # Expand it so it's in the form of a + b*s + c*s^2 ...
    expanded_n_m_t = expand(n_m_t)
    return expanded_n_m_t
end

# Run a loop to get all the values of n_m_t
n_m_t = [symbols("nm_$t") for t in 1:n_timesteps+1]
for t in 1:n_timesteps+1
    n_m_t[t] = polyΛsq(t)
end

# Now write the function that returns the roots of the polynomial
function srootfunc(t, Λres)
    # Get the polynomial in clean powers of sq
    subbed_poly = expand(subs(n_m_t[t], Λ=>Λres))
    # Get the polynomial coefficients, using an extremelly inelegant method
    for i in 1:2*ndim+1
        if i == 1
            # Get the s^0 entry
            polycoeff[i] = subs(subbed_poly, sq=>0)
        else
            # Start with all the powers of s
            polycoeff[i] = subbed_poly
            for j in 1:i
                #= Get rid of the entries for the powers lower than what we are
                concerned with =#
                polycoeff[i] = expand(polycoeff[i] - polycoeff[j]*sq^(j-1))
            end
            # Divide out the power that we are concerned with
            polycoeff[i] = expand(polycoeff[i]/(sq^(i-1)))
            # Substitute in sq = 0 to kill the higher powers
            polycoeff[i] = subs(polycoeff[i], sq=>0)
        end
    end
    #= Change the type of the polynomial coefficient from Basic to Complex so it
    can play nicely with the optimisations which require proper numeric types =#
    polycoeffN = lambdify(polycoeff)
    root_vec = roots(polycoeffN)
    return root_vec
end

# Get the largest real part of the roots
function srootmax(t, Λres)
    #= Optim.jl only has minimisation methods, so we need to multiply the maximum
    value by -1 so that we find the minimum of -s, which I think is the same as
    maximising s =#
    max_root = -1.0*maximum(real(srootfunc(t, Λres)))
    return max_root
end

#= Now find the value of Λ that maximises sq (we need to phrase as a
minimisation problem because this is all that's in Optim) =#
function smaxfunc(t, Λmax)
    results = optimize(Λres -> srootmax(t, Λres), dmin, Λmax, Brent(); iterations=optimisation_iterations, abs_tol=optimisation_x_tol)
    minimiser = Optim.minimizer(results)
    return minimiser
end

# Now get the values at each step
function maximisation_run(Λmax)
    Λmax_vec = zeros(Float64, (n_timesteps+1))
    for t in 1:n_timesteps+1
        if t == 1
            Λmax_vec[t] = smaxfunc(t, Λmax)
        else
            Λmax_vec[t] = smaxfunc(t, Λmax)
        end
    end
    return Λmax_vec
end

Λvec = maximisation_run(Λmax_limit)

# Calculate whether we actually have a localisation
function localisation_truth(Λvec)
    truth_result = falses(n_timesteps+1)
    for t in 1:n_timesteps+1
        # Invert the sense because of how our optimisation works
        if srootmax(t, Λvec[:, t]) < 0.0
            truth_result[t] = true
        end
    end
    return truth_result
end

localisation_truth_vec = localisation_truth(nΛvec)

# Save it
save_filenameΛ = save_folder * "Λ_results.jld2"
save(save_filenameΛ, Dict("nΛres"=>nΛvec, "loc_truth_res"=>localisation_truth))


text:
#=
    Example from README.md file
    use heta_update_dev()
=#

using HetaSimulator, Plots

# set the absolute or relative path to the project directory
platform = load_platform("$HetaSimulatorDir/cases/story_0", units_check=true)
# wait for the model compilation...

# get the base Heta model
model = platform.models[:nameless]

# single simulation
results = Scenario(model; tspan = (0., 1200.)) |> sim
#results = Scenario(model; tspan = (0., 1200.), parameters=[:k1=>0.1]) |> sim
plotd = plot(results)
# savefig(plotd,"file.png")

# translate to data frame
df = DataFrame(results)


text:
push!(LOAD_PATH,"../../src/")
import UNSflow

AR = 6.

alphadef = UNSflow.ConstDef(45. *pi/180)

hdef = UNSflow.ConstDef(0.)

udef = UNSflow.ConstDef(1.)

full_kinem = UNSflow.KinemDef(alphadef, hdef, udef)

pvt = 0.25

geometry = "FlatPlate"

lespcrit = [0.1;]

surf = UNSflow.ThreeDSurfSimple(AR, full_kinem, geometry, pvt, lespcrit)

curfield = UNSflow.ThreeDFieldSimple()

dtstar = 0.015

t_tot = 15.

nsteps = Int(round(t_tot/dtstar))+1

startflag = 0

writeflag = 1

writeInterval = t_tot/8.

#delvort = UNSflow.delSpalart(500, 12, 1e-5)
delvort = UNSflow.delNone()

mat, surf, curfield = UNSflow.QSLLT_ldvm(surf, curfield, nsteps, dtstar,startflag, writeflag, writeInterval, delvort)

#UNSflow.makeVortPlots3D()

#UNSflow.makeForcePlots3D()

UNSflow.makeInfoPlots3D()

UNSflow.cleanWrite()