text:
#script that estimates the dir bin gas network model on emid data and evaluates GAS
# forecasting performances

using Utilities,AReg,StaticNets,JLD,MLBase,StatsBase,DynNets
using PyCall; pygui(:qt); using PyPlot
#estimate and save for half the datase (after LTRO) or whole data?
halfPeriod = false

## Load dataj
fold_Path =  "/home/Domenico/Dropbox/Dynamic_Networks/data/emid_data/juliaFiles/"
loadFilePartialName = "Weekly_eMid_Data_from_"
halfPeriod ? periodEndStr =  "2012_03_12_to_2015_02_27.jld" : periodEndStr =  "2009_06_22_to_2015_02_27.jld"
@load(fold_Path*loadFilePartialName*periodEndStr, AeMidWeekly_T,banksIDs,inactiveBanks,
                YeMidWeekly_T,weekInd,datesONeMid ,degsIO_T,strIO_T)

#density plot
matY_T = YeMidWeekly_T[:,:,3:end]
 matA_T = matY_T.>0
    N = length(matA_T[1,:,1])
    T = size(matA_T)[3]
    Ttrain = 100#round(Int, T/2) #70 106 #
    threshVar = 0.00#5
    degsIO_T = [sumSq(matA_T,2);sumSq(matA_T,1)]
 allFitSS =  StaticNets.estimate( StaticNets.SnapSeqNetDirBin1(degsIO_T);
                        identPost = false,identIter= true )

gasEst_rolling = Array{Vector{Vector{Float64}},1}(T-Ttrain)
#the following variable stores the filtered values of the Parameters
# using at each step the observation t and static gas parameters estimated
# using observations from t-101:t-1
# te values between 1:Ttrain are filtered using the contstant Parameters estimated
# in the first Ttrain observations
gasFiltAndForeFitFromRollEst = zeros(2N,T)
#Estimate gas model on rolling windows of length Ttrain
modGasDirBin1_eMidTrain = DynNets.SdErgmDirBin1(degsIO_T[:,1:Ttrain],"FISHER-DIAG")
gasParEstOnTrain_roll_t,~ = DynNets.estimateTarg(modGasDirBin1_eMidTrain;SSest = allFitSS[:,1:Ttrain] )
gasEst_rolling[1] = gasParEstOnTrain_roll_t
GasforeFit,~ = DynNets.score_driven_filter_or_dgp( DynNets.SdErgmDirBin1(degsIO_T[:,1:Ttrain+1]),
                                    [gasParEstOnTrain_roll_t[1];gasParEstOnTrain_roll_t[2];gasParEstOnTrain_roll_t[3]])
gasforeFit = Float64.(GasforeFit)
gasFiltAndForeFitFromRollEst[:,1:Ttrain+1] = gasforeFit[:,1:Ttrain+1]
#

for t = 2:T-Ttrain
    @show(t)
 modGasDirBin1_eMidTrain = DynNets.SdErgmDirBin1(degsIO_T[:,t:Ttrain+t],"FISHER-DIAG")
 gasParEstOnTrain_roll_t,~ = DynNets.estimateTarg(modGasDirBin1_eMidTrain;SSest = allFitSS[:,1+t:Ttrain+t] )
 gasEst_rolling[t] = gasParEstOnTrain_roll_t
 GasforeFit,~ = DynNets.score_driven_filter_or_dgp( DynNets.SdErgmDirBin1(degsIO_T[:,t:Ttrain+t]),
                                     [gasParEstOnTrain_roll_t[1];gasParEstOnTrain_roll_t[2];gasParEstOnTrain_roll_t[3]])
 gasforeFit = Float64.(GasforeFit)
 # add the last value from this sequence of filtered ones (for which the rolling gas estimates are used) to the sequence of filtered parameters
 gasFiltAndForeFitFromRollEst[:,Ttrain+t] = gasforeFit[:,Ttrain+1]
 end

@save(fold_Path*"_rolling_SD_estimates_"*periodEndStr, gasFiltAndForeFitFromRollEst,allFitSS,gasEst_rolling)

max_stars_repo_name: domenicodigangi/SDERGM

max_issues_repo_name: domenicodigangi/SDERGM

max_forks_repo_name: domenicodigangi/SDERGM


text:
precedence(op::Int) = op < Tokens.end_assignments ?  AssignmentOp :
                       op < Tokens.end_pairarrow ? 2 :
                       op < Tokens.end_conditional ? ConditionalOp :
                       op < Tokens.end_arrow ?       ArrowOp :
                       op < Tokens.end_lazyor ?      LazyOrOp :
                       op < Tokens.end_lazyand ?     LazyAndOp :
                       op < Tokens.end_comparison ?  ComparisonOp :
                       op < Tokens.end_pipe ?        PipeOp :
                       op < Tokens.end_colon ?       ColonOp :
                       op < Tokens.end_plus ?        PlusOp :
                       op < Tokens.end_bitshifts ?   BitShiftOp :
                       op < Tokens.end_times ?       TimesOp :
                       op < Tokens.end_rational ?    RationalOp :
                       op < Tokens.end_power ?       PowerOp :
                       op < Tokens.end_decl ?        DeclarationOp :
                       op < Tokens.end_where ?       WhereOp : DotOp

precedence(kind::Tokens.Kind) = kind === Tokens.DDDOT ? DddotOp :
                        kind < Tokens.begin_assignments ? 0 :
                        kind < Tokens.end_assignments ?   AssignmentOp :
                        kind < Tokens.end_pairarrow ?   2 :
                       kind < Tokens.end_conditional ?    ConditionalOp :
                       kind < Tokens.end_arrow ?          ArrowOp :
                       kind < Tokens.end_lazyor ?         LazyOrOp :
                       kind < Tokens.end_lazyand ?        LazyAndOp :
                       kind < Tokens.end_comparison ?     ComparisonOp :
                       kind < Tokens.end_pipe ?           PipeOp :
                       kind < Tokens.end_colon ?          ColonOp :
                       kind < Tokens.end_plus ?           PlusOp :
                       kind < Tokens.end_bitshifts ?      BitShiftOp :
                       kind < Tokens.end_times ?          TimesOp :
                       kind < Tokens.end_rational ?       RationalOp :
                       kind < Tokens.end_power ?          PowerOp :
                       kind < Tokens.end_decl ?           DeclarationOp :
                       kind < Tokens.end_where ?          WhereOp :
                       kind < Tokens.end_dot ?            DotOp :
                       kind === Tokens.ANON_FUNC ? AnonFuncOp :
                       kind === Tokens.PRIME ?             PrimeOp : 20

precedence(x) = 0
precedence(x::AbstractToken) = precedence(kindof(x))
precedence(x::EXPR) = precedence(kindof(x))

isoperator(kind) = Tokens.begin_ops < kind < Tokens.end_ops
isoperator(t::AbstractToken) = isoperator(kindof(t))
isoperator(x::EXPR) = typof(x) === OPERATOR

isunaryop(op) = false
isunaryop(op::EXPR) = isoperator(op) && isunaryop(kindof(op))
isunaryop(t::AbstractToken) = isunaryop(kindof(t))
@static if VERSION < v"1.2.0"
    isunaryop(kind::Tokens.Kind) = kind === Tokens.ISSUBTYPE ||
                    kind === Tokens.ISSUPERTYPE ||
                    kind === Tokens.PLUS ||
                    kind === Tokens.MINUS ||
                    kind === Tokens.NOT ||
                    kind === Tokens.APPROX ||
                    kind === Tokens.NOT_SIGN ||
                    kind === Tokens.AND ||
                    kind === Tokens.SQUARE_ROOT ||
                    kind === Tokens.CUBE_ROOT ||
                    kind === Tokens.QUAD_ROOT ||
                    kind === Tokens.DECLARATION ||
                    kind === Tokens.EX_OR ||
                    kind === Tokens.COLON
else
    isunaryop(kind::Tokens.Kind) = kind === Tokens.ISSUBTYPE ||
                    kind === Tokens.ISSUPERTYPE ||
                    kind === Tokens.PLUS ||
                    kind === Tokens.MINUS ||
                    kind === Tokens.NOT ||
                    kind === Tokens.APPROX ||
                    kind === Tokens.NOT_SIGN ||
                    kind === Tokens.AND ||
                    kind === Tokens.SQUARE_ROOT ||
                    kind === Tokens.CUBE_ROOT ||
                    kind === Tokens.QUAD_ROOT ||
                    kind === Tokens.DECLARATION ||
                    kind === Tokens.EX_OR ||
                    kind === Tokens.COLON ||
                    kind === Tokens.STAR_OPERATOR
end

isunaryandbinaryop(t) = false
isunaryandbinaryop(t::AbstractToken) = isunaryandbinaryop(kindof(t))
@static if VERSION < v"1.2.0"
    isunaryandbinaryop(kind::Tokens.Kind) = kind === Tokens.PLUS ||
                            kind === Tokens.MINUS ||
                            kind === Tokens.EX_OR ||
                            kind === Tokens.ISSUBTYPE ||
                            kind === Tokens.ISSUPERTYPE ||
                            kind === Tokens.AND ||
                            kind === Tokens.APPROX ||
                            kind === Tokens.DECLARATION ||
                            kind === Tokens.COLON
else
    isunaryandbinaryop(kind::Tokens.Kind) = kind === Tokens.PLUS ||
                            kind === Tokens.MINUS ||
                            kind === Tokens.EX_OR ||
                            kind === Tokens.ISSUBTYPE ||
                            kind === Tokens.ISSUPERTYPE ||
                            kind === Tokens.AND ||
                            kind === Tokens.APPROX ||
                            kind === Tokens.DECLARATION ||
                            kind === Tokens.COLON ||
                            kind === Tokens.STAR_OPERATOR
end

isbinaryop(op) = false
isbinaryop(op::EXPR) = isoperator(op) && isbinaryop(kindof(op))
isbinaryop(t::AbstractToken) = isbinaryop(kindof(t))
isbinaryop(kind::Tokens.Kind) = isoperator(kind) &&
                    !(kind === Tokens.SQUARE_ROOT ||
                    kind === Tokens.CUBE_ROOT ||
                    kind === Tokens.QUAD_ROOT ||
                    kind === Tokens.NOT ||
                    kind === Tokens.NOT_SIGN)

function non_dotted_op(t::AbstractToken)
    k = kindof(t)
    return (k === Tokens.COLON_EQ ||
            k === Tokens.PAIR_ARROW ||
            k === Tokens.EX_OR_EQ ||
            k === Tokens.CONDITIONAL ||
            k === Tokens.LAZY_OR ||
            k === Tokens.LAZY_AND ||
            k === Tokens.ISSUBTYPE ||
            k === Tokens.ISSUPERTYPE ||
            k === Tokens.LPIPE ||
            k === Tokens.RPIPE ||
            k === Tokens.EX_OR ||
            k === Tokens.COLON ||
            k === Tokens.DECLARATION ||
            k === Tokens.IN ||
            k === Tokens.ISA ||
            k === Tokens.WHERE ||
            (isunaryop(k) && !isbinaryop(k) && !(k === Tokens.NOT)))
end

issyntaxcall(op) = false
function issyntaxcall(op::EXPR)
    K = kindof(op)
    P = precedence(K)
    P == AssignmentOp && !(K === Tokens.APPROX || K === Tokens.PAIR_ARROW) ||
    K === Tokens.RIGHT_ARROW ||
    P == LazyOrOp ||
    P == LazyAndOp ||
    K === Tokens.ISSUBTYPE ||
    K === Tokens.ISSUPERTYPE ||
    K === Tokens.COLON ||
    K === Tokens.DECLARATION ||
    K === Tokens.DOT ||
    K === Tokens.DDDOT ||
    K === Tokens.PRIME ||
    K === Tokens.WHERE ||
    K === Tokens.ANON_FUNC
end

issyntaxunarycall(op) = false
function issyntaxunarycall(op::EXPR)
    K = kindof(op)
    !op.dot && (K === Tokens.EX_OR ||
    K === Tokens.AND ||
    K === Tokens.DECLARATION ||
    K === Tokens.DDDOT ||
    K === Tokens.PRIME ||
    K === Tokens.ISSUBTYPE ||
    K === Tokens.ISSUPERTYPE)
end

LtoR(prec::Int) = AssignmentOp ≤ prec ≤ LazyAndOp || prec == PowerOp

"""
    parse_unary(ps)

Having hit a unary operator at the start of an expression return a call.
"""
function parse_unary(ps::ParseState, op::EXPR)
    K, dot = kindof(op), op.dot
    if is_colon(op)
        ret = parse_unary_colon(ps, op)
    elseif should_negate_number_literal(ps, op)
        arg = mLITERAL(next(ps))
        ret = mLITERAL(op.fullspan + arg.fullspan, (op.fullspan + arg.span), string(is_plus(op) ? "+" : "-", val(ps.t, ps)), kindof(ps.t))
    else
        P = precedence(K)
        prec = P == DeclarationOp ? DeclarationOp :
                    K === Tokens.AND ? DeclarationOp :
                    K === Tokens.EX_OR ? 20 : PowerOp
        arg = @closer ps :unary @precedence ps prec parse_expression(ps)
        ret = mUnaryOpCall(op, arg)
    end
    return ret
end

function parse_unary_colon(ps::ParseState, op::EXPR)
    op = requires_no_ws(op, ps)
    if Tokens.iskeyword(kindof(ps.nt))
        ret = EXPR(Quotenode, EXPR[op, mIDENTIFIER(next(ps))])
    elseif isidentifier(ps.nt)
        id = INSTANCE(next(ps))
        if VERSION > v"1.3.0-" && valof(id) == "var" && isemptyws(ps.ws) && (ps.nt.kind === Tokens.STRING || ps.nt.kind === Tokens.TRIPLE_STRING)
            # Special case for :var"sdf"
            arg = parse_string_or_cmd(next(ps), id)
            id = EXPR(NONSTDIDENTIFIER, EXPR[id, arg])
        end
        ret = EXPR(Quotenode, EXPR[op, id])
    elseif Tokens.begin_literal < kindof(ps.nt) < Tokens.CHAR ||
        isoperator(kindof(ps.nt)) || isidentifier(ps.nt) || kindof(ps.nt) === Tokens.TRUE || kindof(ps.nt) === Tokens.FALSE
        ret = EXPR(Quotenode, EXPR[op, INSTANCE(next(ps))])
    elseif closer(ps)
        ret = op
    else
        prev_errored = ps.errored
        arg = @precedence ps 20 parse_expression(ps)
        if isbracketed(arg)  && typof(arg.args[2]) === ErrorToken && errorof(arg.args[2]) === UnexpectedAssignmentOp
            ps.errored = prev_errored
            arg.args[2] = arg.args[2].args[1]
            setparent!(arg.args[2], arg)
        end
        ret = EXPR(Quote, EXPR[op, arg])
    end
    return ret
end

function parse_operator_eq(ps::ParseState, ret::EXPR, op::EXPR)
    nextarg = @precedence ps AssignmentOp - LtoR(AssignmentOp) parse_expression(ps)

    if is_func_call(ret) && !(isbeginorblock(nextarg))
        nextarg = EXPR(Block, EXPR[nextarg])
    end
    ret = mBinaryOpCall(ret, op, nextarg)
    return ret
end

# Parse conditionals
function parse_operator_cond(ps::ParseState, ret::EXPR, op::EXPR)
    ret = requires_ws(ret, ps)
    op = requires_ws(op, ps)
    nextarg = @closer ps :ifop parse_expression(ps)
    if ps.nt.kind !== Tokens.COLON
        op2 = mErrorToken(ps, mOPERATOR(0, 0, Tokens.COLON, false), MissingColon)
    else
        op2 = requires_ws(mOPERATOR(next(ps)), ps)
    end

    nextarg2 = @closer ps :comma @precedence ps 0 parse_expression(ps)

    return EXPR(ConditionalOpCall, EXPR[ret, op, nextarg, op2, nextarg2])
end

# Parse comparisons
function parse_comp_operator(ps::ParseState, ret::EXPR, op::EXPR)
    nextarg = @precedence ps ComparisonOp - LtoR(ComparisonOp) parse_expression(ps)

    if typof(ret) === Comparison
        push!(ret, op)
        push!(ret, nextarg)
    elseif can_become_comparison(ret)
        ret = EXPR(Comparison, EXPR[ret.args[1], ret.args[2], ret.args[3], op, nextarg])
    else
        ret = mBinaryOpCall(ret, op, nextarg)
    end
    return ret
end

# Parse ranges
function parse_operator_colon(ps::ParseState, ret::EXPR, op::EXPR)
    if isnewlinews(ps.ws) && !ps.closer.paren
        op = mErrorToken(ps, op, UnexpectedNewLine)
    end
    nextarg = @precedence ps ColonOp - LtoR(ColonOp) parse_expression(ps)

    if isbinarycall(ret) && is_colon(ret.args[2])
        ret = EXPR(ColonOpCall, EXPR[ret.args[1], ret.args[2], ret.args[3], op, nextarg])
    else
        ret = mBinaryOpCall(ret, op, nextarg)
    end
    return ret
end

# Parse power (special case for preceding unary ops)
function parse_operator_power(ps::ParseState, ret::EXPR, op::EXPR)
    nextarg = @precedence ps PowerOp - LtoR(PowerOp) @closer ps :inwhere parse_expression(ps)

    if isunarycall(ret)
        nextarg = mBinaryOpCall(ret.args[2], op, nextarg)
        ret = mUnaryOpCall(ret.args[1], nextarg)
    else
        ret = mBinaryOpCall(ret, op, nextarg)
    end
    return ret
end

# parse where
function parse_operator_where(ps::ParseState, ret::EXPR, op::EXPR, setscope=true)
    nextarg = @precedence ps LazyAndOp @closer ps :inwhere parse_expression(ps)

    if typof(nextarg) === Braces
        args = nextarg.args
    else
        args = EXPR[nextarg]
    end
    ret = mWhereOpCall(ret, op, args)
    return ret
end

function parse_operator_dot(ps::ParseState, ret::EXPR, op::EXPR)
    if kindof(ps.nt) === Tokens.LPAREN
        @static if VERSION > v"1.1-"
            iserred = kindof(ps.ws) != Tokens.EMPTY_WS
            sig = @default ps parse_call(ps, ret)
            nextarg = EXPR(TupleH, sig.args[2:end])
            if iserred
                nextarg = mErrorToken(ps, nextarg, UnexpectedWhiteSpace)
            end
        else
            sig = @default ps parse_call(ps, ret)
            nextarg = EXPR(TupleH, sig.args[2:end])
        end
    elseif iskw(ps.nt) || both_symbol_and_op(ps.nt)
        nextarg = mIDENTIFIER(next(ps))
    elseif kindof(ps.nt) === Tokens.COLON
        op2 = mOPERATOR(next(ps))
        if kindof(ps.nt) === Tokens.LPAREN
            nextarg = @closeparen ps @precedence ps DotOp - LtoR(DotOp) parse_expression(ps)
            nextarg = EXPR(Quote, EXPR[op2, nextarg])
        else
            nextarg = @precedence ps DotOp - LtoR(DotOp) parse_unary(ps, op2)
        end
    elseif kindof(ps.nt) === Tokens.EX_OR && kindof(ps.nnt) === Tokens.LPAREN
        op2 = mOPERATOR(next(ps))
        nextarg = parse_call(ps, op2)
    else
        nextarg = @precedence ps DotOp - LtoR(DotOp) parse_expression(ps)
    end

    if isidentifier(nextarg) || isinterpolant(nextarg)
        ret = mBinaryOpCall(ret, op, EXPR(Quotenode, EXPR[nextarg]))
    elseif typof(nextarg) === Vect
        ret = mBinaryOpCall(ret, op, EXPR(Quote, EXPR[nextarg]))
    elseif typof(nextarg) === MacroCall
        mname = mBinaryOpCall(ret, op, EXPR(Quotenode, EXPR[nextarg.args[1]]))
        ret = EXPR(MacroCall, EXPR[mname])
        for i = 2:length(nextarg.args)
            push!(ret, nextarg.args[i])
        end
    else
        ret = mBinaryOpCall(ret, op, nextarg)
    end
    return ret
end

function parse_operator_anon_func(ps::ParseState, ret::EXPR, op::EXPR)
    arg = @closer ps :comma @precedence ps 0 parse_expression(ps)

    if !isbeginorblock(arg)
        arg = EXPR(Block, EXPR[arg])
    end
    return mBinaryOpCall(ret, op, arg)
end

function parse_operator(ps::ParseState, ret::EXPR, op::EXPR)
    K, dot = kindof(op), op.dot
    P = precedence(K)

    if typof(ret) === ChainOpCall && (is_star(op) || is_plus(op)) && kindof(op) == kindof(ret.args[2])
        nextarg = @precedence ps P - LtoR(P) parse_expression(ps)
        push!(ret, op)
        push!(ret, nextarg)
        ret = ret
    elseif can_become_chain(ret, op)
        nextarg = @precedence ps P - LtoR(P) parse_expression(ps)
        ret = EXPR(ChainOpCall, EXPR[ret.args[1], ret.args[2], ret.args[3], op, nextarg])
    elseif is_eq(op)
        ret = parse_operator_eq(ps, ret, op)
    elseif is_cond(op)
        ret = parse_operator_cond(ps, ret, op)
    elseif is_colon(op)
        ret = parse_operator_colon(ps, ret, op)
    elseif is_where(op)
        ret = parse_operator_where(ps, ret, op)
    elseif is_anon_func(op)
        ret = parse_operator_anon_func(ps, ret, op)
    elseif is_dot(op)
        ret = parse_operator_dot(ps, ret, op)
    elseif is_dddot(op) || is_prime(op)
        ret = mUnaryOpCall(ret, op)
    elseif P == ComparisonOp
        ret = parse_comp_operator(ps, ret, op)
    elseif P == PowerOp
        ret = parse_operator_power(ps, ret, op)
    else
        ltor = K === Tokens.LPIPE ? true : LtoR(P)
        nextarg = @precedence ps P - ltor parse_expression(ps)
        ret = mBinaryOpCall(ret, op, nextarg)
    end
    return ret
end

ext: jl


text:
let
    method = ME_MAX_POL_B0
    objider = iJR.BIOMASS_IDER

    # handle cache
    datfile = dat_file(;method)
    check_cache(;method) && return

    # setup
    model = iJR.load_model("max_model")

    lock(WLOCK) do
        @info("Doing... ",
            method, threadid()
        ); println()
    end

    # maxent
    epout = ChEP.maxent_ep(model;
        alpha = Inf, damp = 0.9, epsconv = 1e-4,
        verbose = true, maxiter = 5000
    )

    # storing
    lock(WLOCK) do
        # Storing
        dat = Dict()
        dat[:exp_beta] = 0.0
        dat[:epouts] = Dict(0.0 => epout)
        dat[:model] = model |> ChU.compressed_model

        # caching
        serialize(datfile, dat)

        @info("Finished ", threadid())
        println()
    end
end


text:
using AccurateArithmetic
using Test

using AccurateArithmetic.Summation: accumulate
using AccurateArithmetic.Summation: sumAcc, compSumAcc, mixedSumAcc
using AccurateArithmetic.Summation: dotAcc, compDotAcc, mixedDotAcc
using AccurateArithmetic.Test: generate_sum, generate_dot
using LinearAlgebra
using StableRNGs

@testset "AccurateArithmetic" begin
    rng = StableRNG(42)

    @testset "Tests" begin
        @testset "generate_sum" begin
            @testset "vector length" begin
                @testset "F64" begin
                    for n in 100:110
                        let (x, s, c_) = generate_sum(n, 1e10; rng=rng)
                            @test eltype(x) === Float64
                            @test s isa Float64
                            @test length(x) == n
                        end
                    end
                end
                @testset "F32" begin
                    for n in 100:110
                        let (x, s, c_) = generate_sum(n, 1f5; rng=rng)
                            @test eltype(x) === Float32
                            @test s isa Float32
                            @test length(x) == n
                        end
                    end
                end
            end

            @testset "condition number" begin
                @testset "F64" begin
                    for c in (1e10, 1e20)
                        let (x, s, c_) = generate_sum(101, c; rng=rng)
                            @test eltype(x) === Float64
                            @test s isa Float64
                            @test c/3 < c_ < 3c
                        end
                    end
                end
                @testset "F32" begin
                    for c in (1f5, 1f10)
                        let (x, s, c_) = generate_sum(101, c; rng=rng)
                            @test eltype(x) === Float32
                            @test s isa Float32
                            @test c/3 < c_ < 3c
                        end
                    end
                end
            end
        end

        @testset "generate_dot" begin
            @testset "vector length" begin
                @testset "F64" begin
                    for n in 100:110
                        let (x, y, s, c_) = generate_dot(n, 1e10; rng=rng)
                            @test eltype(x) === Float64
                            @test s isa Float64
                            @test length(x) == n
                            @test length(y) == n
                        end
                    end
                end
                @testset "F32" begin
                    for n in 100:110
                        let (x, y, s, c_) = generate_dot(n, 1f5; rng=rng)
                            @test eltype(x) === Float32
                            @test s isa Float32
                            @test length(x) == n
                            @test length(y) == n
                        end
                    end
                end
            end

            @testset "condition number" begin
                @testset "F64" begin
                    for c in (1e10, 1e20)
                        let (x, y, s, c_) = generate_dot(101, c; rng=rng)
                            @test s isa Float64
                            @test c/3 < c_ < 3c
                        end
                    end
                end
                @testset "F32" begin
                    for c in (1f5, 1f10)
                        let (x, y, s, c_) = generate_dot(101, c; rng=rng)
                            @test s isa Float32
                            @test c/3 < c_ < 3c
                        end
                    end
                end
            end
        end
    end

    @testset "summation" begin
        @testset "naive" begin
            @testset "F64" begin
                for N in 100:110
                    x = rand(rng, N)
                    ref = sum(x)
                    @test ref isa Float64
                    @test ref ≈ sum_naive(x)

                    acc = sumAcc
                    @test ref ≈ accumulate((x,), acc, Val(:scalar), Val(2))
                    @test ref ≈ accumulate((x,), acc, Val(:mask),   Val(2))
                end
            end
            @testset "F32" begin
                for N in 100:110
                    x = rand(rng, Float32, N)
                    ref = sum(x)
                    @test ref isa Float32
                    @test ref ≈ sum_naive(x)

                    acc = sumAcc
                    @test ref ≈ accumulate((x,), acc, Val(:scalar), Val(2))
                    @test ref ≈ accumulate((x,), acc, Val(:mask),   Val(2))
                end
            end
        end

        @testset "compensated" begin
            @testset "F64" begin
                for N in 100:110
                    x, ref, _ = generate_sum(N, 1e10; rng=rng)
                    @test ref isa Float64
                    @test ref == sum_oro(x)
                    @test ref == sum_kbn(x)

                    acc = compSumAcc(two_sum)
                    @test ref == accumulate((x,), acc, Val(:scalar), Val(2))
                    @test ref == accumulate((x,), acc, Val(:mask),   Val(2))
                end
            end
            @testset "F32" begin
                for N in 100:110
                    x, ref, _ = generate_sum(N, 1f5; rng=rng)
                    @test ref isa Float32
                    @test ref == sum_oro(x)
                    @test ref == sum_kbn(x)

                    acc = compSumAcc(two_sum)
                    @test ref == accumulate((x,), acc, Val(:scalar), Val(2))
                    @test ref == accumulate((x,), acc, Val(:mask),   Val(2))
                end
            end
        end

        @testset "mixed" begin
            for N in 100:110
                # Only test for approximate equality here, since sum_mixed
                # returns a Float64 whereas the reference value is a Float32

                x, ref, _ = generate_sum(N, 1f7; rng=rng)
                @test ref isa Float32
                @test ref ≈ sum_mixed(x)

                acc = mixedSumAcc
                @test ref ≈ accumulate((x,), acc, Val(:scalar), Val(2))
                @test ref ≈ accumulate((x,), acc, Val(:mask),   Val(2))
            end
        end
    end

    @testset "dot product" begin
        @testset "naive" begin
            @testset "F64" begin
                for N in 100:110
                    x = rand(rng, N)
                    y = rand(rng, N)
                    ref = dot(x, y)
                    @test ref isa Float64
                    @test ref ≈ dot_naive(x, y)

                    acc = dotAcc
                    @test ref ≈ accumulate((x,y), acc, Val(:scalar), Val(2))
                    @test ref ≈ accumulate((x,y), acc, Val(:mask),   Val(2))
                end
            end
            @testset "F32" begin
                for N in 100:110
                    x = rand(rng, Float32, N)
                    y = rand(rng, Float32, N)
                    ref = dot(x, y)
                    @test ref isa Float32
                    @test ref ≈ dot_naive(x, y)

                    acc = dotAcc
                    @test ref ≈ accumulate((x,y), acc, Val(:scalar), Val(2))
                    @test ref ≈ accumulate((x,y), acc, Val(:mask),   Val(2))
                end
            end
        end

        @testset "compensated" begin
            @testset "F64" begin
                for N in 100:110
                    x, y, ref, _ = generate_dot(N, 1e10; rng=rng)
                    @test ref isa Float64
                    @test ref == dot_oro(x, y)

                    acc = compDotAcc
                    @test ref == accumulate((x,y), acc, Val(:scalar), Val(2))
                    @test ref == accumulate((x,y), acc, Val(:mask),   Val(2))
                end
            end
            @testset "F32" begin
                for N in 100:110
                    x, y, ref, _ = generate_dot(N, 1f5; rng=rng)
                    @test ref isa Float32
                    @test ref == dot_oro(x, y)

                    acc = compDotAcc
                    @test ref == accumulate((x,y), acc, Val(:scalar), Val(2))
                    @test ref == accumulate((x,y), acc, Val(:mask),   Val(2))
                end
            end
        end

        @testset "mixed" begin
            for N in 100:110
                # Only test for approximate equality here, since dot_mixed
                # returns a Float64 whereas the reference value is a Float32

                x, y, ref, _ = generate_dot(N, 1f7; rng=rng)
                @test ref isa Float32
                @test ref ≈ dot_mixed(x, y)

                acc = mixedDotAcc
                @test ref ≈ accumulate((x,y), acc, Val(:scalar), Val(2))
                @test ref ≈ accumulate((x,y), acc, Val(:mask),   Val(2))
            end
        end
    end
end

using BenchmarkTools
BLAS.set_num_threads(1)
BenchmarkTools.DEFAULT_PARAMETERS.evals = 1000

println("\nsize 32")
x = rand(32)
y = rand(32)
print("  sum_kbn  "); @btime sum_kbn($x)
print("  sum_oro  "); @btime sum_oro($x)
print("  sum_naive"); @btime sum_naive($x)
print("  sum      "); @btime sum($x)
println()
print("  dot_oro  "); @btime dot_oro($x, $y)
print("  dot_naive"); @btime dot_naive($x, $y)
print("  dot      "); @btime dot($x, $y)

println("\nsize 10_000")
x = rand(10_000)
y = rand(10_000)
print("  sum_kbn  "); @btime sum_kbn($x)
print("  sum_oro  "); @btime sum_oro($x)
print("  sum_naive"); @btime sum_naive($x)
print("  sum      "); @btime sum($x)
println()
print("  dot_oro  "); @btime dot_oro($x, $y)
print("  dot_naive"); @btime dot_naive($x, $y)
print("  dot      "); @btime dot($x, $y)

BenchmarkTools.DEFAULT_PARAMETERS.evals = 10
println("\nsize 1_000_000")
x = rand(1_000_000)
y = rand(1_000_000)
print("  sum_kbn  "); @btime sum_kbn($x)
print("  sum_oro  "); @btime sum_oro($x)
print("  sum_naive"); @btime sum_naive($x)
print("  sum      "); @btime sum($x)
println()
print("  dot_oro  "); @btime dot_oro($x, $y)
print("  dot_naive"); @btime dot_naive($x, $y)
print("  dot      "); @btime dot($x, $y)

println("\nsize 100_000_000")
x = rand(100_000_000)
y = rand(100_000_000)
print("  sum_kbn  "); @btime sum_kbn($x)
print("  sum_oro  "); @btime sum_oro($x)
print("  sum_naive"); @btime sum_naive($x)
print("  sum      "); @btime sum($x)
println()
print("  dot_oro  "); @btime dot_oro($x, $y)
print("  dot_naive"); @btime dot_naive($x, $y)
print("  dot      "); @btime dot($x, $y)


text:
module Test_cfid

using DescriptorSystems
using LinearAlgebra
using Polynomials
using Test

println("Test_cfid")
@testset "grcfid and glcfid" begin

@testset "grcfid" begin

sys = rdss(0,0,0);
@time sysni, sysmi = grcfid(sys);
@test gnrank(sys*sysmi-sysni) == 0   &&   #  G(s)*Mi(s)-Ni(s) = 0
      gnrank(sysmi'*sysmi-I) == 0  && # conj(Mi(s))*Mi(s)-I = 0
      isproper(sysmi) & isproper(sysni) && # checking properness of factors
      isstable(sysmi,atol=1.e-7) &&
      isstable(sysni,atol=1.e-7)  &&
      isempty(gzero(gminreal([sysni;sysmi],atol=1.e-7)))

a = -(1. +2im)*ones(1,1); e = (1. -im)*ones(1,1); b = (2. +im)*ones(1,1);
c = (-1. +im)*ones(1,1); d = complex(0)*ones(1,1);
sys = dss(a,e,b,c,d);
@time sysni, sysmi = grcfid(sys,atol=1.e-7,atol3=1.e-7);
@test gnrank(sys*sysmi-sysni,atol=1.e-7) == 0 &&
      gnrank(sysmi'*sysmi-I,atol=1.e-7) == 0  &&
      isproper(sysmi) && isproper(sysni) &&
      isstable(sysmi,atol=1.e-7) &&
      isstable(sysni,atol=1.e-7)  &&
      isempty(gzero(gminreal([sysni;sysmi],atol=1.e-7)))

a = -(1. +2im)*ones(1,1); e = (1. -im)*ones(1,1); b = (2. +im)*ones(1,1);
c = (-1. +im)*ones(1,1); d = complex(0)*ones(1,1);
sys = dss(a,e,b,c,d,Ts=1);
@time sysni, sysmi = grcfid(sys,atol=1.e-7,atol3=1.e-7);
@test gnrank(sys*sysmi-sysni,atol=1.e-7) == 0 &&
      gnrank(sysmi'*sysmi-I,atol=1.e-7) == 0  &&
      isproper(sysmi) && isproper(sysni) &&
      isstable(sysmi,atol=1.e-7) &&
      isstable(sysni,atol=1.e-7)  &&
      isempty(gzero(gminreal([sysni;sysmi],atol=1.e-7)))

s = Polynomial([0, 1],'s'); z = Polynomial([0, 1],'z');     # define the complex variables s and z

# Gd = [z^2 z/(z-2); 0 1/z]     # define the 2-by-2 improper Gd(z)
Nd = [z^2 z; 0 1]; Dd = [1 z-2; 1 z];

# build LTI (nonminimal) descriptor realizations of Gc(s) and Gd(z)
sysd = dss(Nd,Dd,Ts = 1);      # build discrete-time descriptor system realization

sysni, sysmi = grcfid(sysd);
@test gnrank(sysd*sysmi-sysni,atol=1.e-7) == 0  &&    #  Gd(z)*M(z)-N(z) = 0
      iszero(sysmi'-inv(sysmi),atol=1.e-7)  && # conj(M(z))*M(z)-I = 0
      isproper(sysmi) & isproper(sysni)  &&
      isstable(sysmi,atol=1.e-7) &&
      isstable(sysni,atol=1.e-7)  &&
      isempty(gzero(gminreal([sysni;sysmi],atol=1.e-7)))

s = Polynomial([0, 1],'s');
gn = [(s-1) s 1;
    0 (s-2) (s-2);
    (s-1) (s^2+2*s-2) (2*s-1)];
gd = [(s+2) (s+2) (s+2);
    1 (s+1)^2 (s+1)^2;
    (s+2) (s+1)*(s+2) (s+1)*(s+2)];

sys = dss(gn,gd,minimal = true, atol = 1.e-7);
sysni, sysmi = grcfid(sys,mindeg=true);
@test gnrank(sys*sysmi-sysni,atol=1.e-7) == 0  &&    #  Gd(z)*M(z)-N(z) = 0
      iszero(sysmi'-inv(sysmi),atol=1.e-7)  && # conj(M(z))*M(z)-I = 0
      isproper(sysmi) & isproper(sysni)  &&
      isstable(sysmi,atol=1.e-7) &&
      isstable(sysni,atol=1.e-7)  &&
      isempty(gzero(gminreal([sysni;sysmi],atol=1.e-7))) &&
      order(sysmi) == 0

sys = sys'
sysni, sysmi = grcfid(sys);
@test gnrank(sys*sysmi-sysni,atol=1.e-7) == 0  &&    #  Gd(z)*M(z)-N(z) = 0
      iszero(sysmi'-inv(sysmi),atol=1.e-7)  && # conj(M(z))*M(z)-I = 0
      isproper(sysmi) & isproper(sysni)  &&
      isstable(sysmi,atol=1.e-7) &&
      isstable(sysni,atol=1.e-7)  &&
      isempty(gzero(gminreal([sysni;sysmi],atol=1.e-7))) &&
      order(sysmi) == 4

s = Polynomial([0, 1],'s');
gn = [(s-1) s 1;
    0 (s-2) (s-2);
    (s-1) (s^2+2*s-2) (2*s-1)];
gd = [(s-2) (s+2) (s+2);
    1 (s+1)^2 (s+1)^2;
    (s+2) (s-1)*(s+2) (s+1)*(s-2)];

sys = dss(gn,gd,minimal = true, atol = 1.e-7);
sysni, sysmi = grcfid(sys,mindeg = true);
@test gnrank(sys*sysmi-sysni,atol=1.e-7) == 0  &&    #  Gd(z)*M(z)-N(z) = 0
      iszero(sysmi'-inv(sysmi),atol=1.e-7)  && # conj(M(z))*M(z)-I = 0
      isproper(sysmi) & isproper(sysni)  &&
      isstable(sysmi,atol=1.e-7) &&
      isstable(sysni,atol=1.e-7)  &&
      isempty(gzero(gminreal([sysni;sysmi],atol=1.e-7))) &&
      order(sysmi) == 3

z = Polynomial([0, 1],'z')
g = [z^2+z+1 4*z^2+3*z+2 2*z^2-2;
    z 4*z-1 2*z-2;
    z^2 4*z^2-z 2*z^2-2*z];

sys = dss(g,minimal = true, atol = 1.e-7,Ts = 1);
sysni, sysmi = grcfid(sys,mindeg = true,atol=1.e-7,atol3=1.e-7);
@test gnrank(sys*sysmi-sysni,atol=1.e-7) == 0  &&    #  Gd(z)*M(z)-N(z) = 0
      iszero(sysmi'-inv(sysmi),atol=1.e-7)  && # conj(M(z))*M(z)-I = 0
      isproper(sysmi) & isproper(sysni)  &&
      isstable(sysmi,atol=1.e-7) &&
      isstable(sysni,atol=1.e-7)  &&
      isempty(gzero(gminreal([sysni;sysmi],atol=1.e-7)))

z = Polynomial([0, 1],'z')
# Eigenvalue(s) on the unit circle
gn = 1
gd = z-1
sys = dss(gn,gd,Ts = 1)
try
   sysni, sysmi = grcfid(sys)
   @test false
catch
   @test true
end

fast = true; Ty = Complex{Float64}; #Ty = Float64
n = 5; p = 3; m = 2;

for fast in (true, false)
# random examples
for Ty in (Float64, Complex{Float64})

# continuous, standard
sys = rss(n,p,m,T = Ty,disc=false);
@time sysni, sysmi = grcfid(sys, fast = fast,atol=1.e-7,atol3=1.e-7);
@test gnrank(sys*sysmi-sysni,atol=1.e-7) == 0   &&
      iszero(sysmi'-inv(sysmi),atol=1.e-7)   &&
      isproper(sysmi) && isproper(sysni) &&
      isstable(sysmi,atol=1.e-7) &&
      isstable(sysni,atol=1.e-7)  &&
      isempty(gzero(gminreal([sysni;sysmi],atol=1.e-7),atol=1.e-7))

# discrete, standard
sys = rss(n,p,m,T = Ty,disc=true);
@time sysni, sysmi = grcfid(sys, fast = fast,atol=1.e-7,atol3=1.e-7);
@test gnrank(sys*sysmi-sysni,atol=1.e-7) == 0   &&
      iszero(sysmi'-inv(sysmi),atol=1.e-7)   &&
      isproper(sysmi) && isproper(sysni) &&
      isstable(sysmi,atol=1.e-7) &&
      isstable(sysni,atol=1.e-7)  &&
      isempty(gzero(gminreal([sysni;sysmi],atol=1.e-7),atol=1.e-7))

# continuous, standard, uncontrollable
sys = rss(n,p,m,T = Ty,disc=false, nuc = 5);
@time sysni, sysmi = grcfid(sys, fast = fast,atol=1.e-7,atol3=1.e-7);
@test gnrank(sys*sysmi-sysni,atol=1.e-7) == 0   &&
      iszero(sysmi'-inv(sysmi),atol=1.e-7)   &&
      isproper(sysmi) && isproper(sysni) &&
      isstable(sysmi,atol=1.e-7) &&
      isstable(sysni,atol=1.e-7)  &&
      isempty(gzero(gminreal([sysni;sysmi],atol=1.e-7),atol=1.e-7))

# discrete, standard, uncontrollable
sys = rss(n,p,m,T = Ty,disc=true, nuc = 5);
@time sysni, sysmi = grcfid(sys, fast = fast,atol=1.e-7,atol3=1.e-7);
@test gnrank(sys*sysmi-sysni,atol=1.e-7) == 0   &&
      iszero(sysmi'-inv(sysmi),atol=1.e-7)   &&
      isproper(sysmi) && isproper(sysni) &&
      isstable(sysmi,atol=1.e-7) &&
      isstable(sysni,atol=1.e-7)  &&
      isempty(gzero(gminreal([sysni;sysmi],atol=1.e-7),atol=1.e-7))

#failure
# continuous, descriptor, no infinite eigenvalues
sys = rdss(n,p,m,T = Ty,disc=false);
@time sysni, sysmi = grcfid(sys, fast = fast,atol=1.e-7,atol3=1.e-7);
@test gnrank(sys*sysmi-sysni,atol=1.e-7) == 0   &&
      gnrank(sysmi'*sysmi-I,atol=1.e-7) == 0  &&
      isproper(sysmi) && isproper(sysni) &&
      isstable(sysmi,atol=1.e-7) &&
      isstable(sysni,atol=1.e-7)  &&
      isempty(gzero(gminreal([sysni;sysmi],atol=1.e-7)))

# discrete, descriptor, no infinite eigenvalues
sys = rdss(n,p,m,T = Ty, disc=true);
@time sysni, sysmi = grcfid(sys, fast = fast,atol=1.e-7,atol3=1.e-7);
@test gnrank(sys*sysmi-sysni,atol=1.e-7) == 0   &&
      gnrank(sysmi'*sysmi-I,atol=1.e-7) == 0  &&
      isproper(sysmi) && isproper(sysni) &&
      isstable(sysmi,atol=1.e-7) &&
      isstable(sysni,atol=1.e-7)  &&
      isempty(gzero(gminreal([sysni;sysmi],atol=1.e-7)))

# continuous, descriptor, proper
sys = rdss(n,p,m,T = Ty,disc=false,id=ones(Int,3));
@time sysni, sysmi = grcfid(sys, fast = fast,atol=1.e-7,atol3=1.e-7);
@test gnrank(sys*sysmi-sysni,atol=1.e-7) == 0   &&
      gnrank(sysmi'*sysmi-I,atol=1.e-7) == 0  &&
      isproper(sysmi) && isproper(sysni) &&
      isstable(sysmi,atol=1.e-7) &&
      isstable(sysni,atol=1.e-7)  &&
      isempty(gzero(gminreal([sysni;sysmi],atol=1.e-7)))

# discrete, descriptor, proper
sys = rdss(n,p,m,T = Ty, disc=true);
@time sysni, sysmi = grcfid(sys, fast = fast,atol=1.e-7,atol3=1.e-7);
@test gnrank(sys*sysmi-sysni,atol=1.e-7) == 0   &&
      gnrank(sysmi'*sysmi-I,atol=1.e-7) == 0  &&
      isproper(sysmi) && isproper(sysni) &&
      isstable(sysmi,atol=1.e-7) &&
      isstable(sysni,atol=1.e-7)  &&
      isempty(gzero(gminreal([sysni;sysmi],atol=1.e-7)))

# continuous, descriptor, proper
sys = rdss(n,p,m,T = Ty,disc=false,id=ones(Int,3));
@time sysni, sysmi = grcfid(sys, fast = fast, mindeg = true, mininf = true, atol=1.e-7,atol3=1.e-7);
@test gnrank(sys*sysmi-sysni,atol=1.e-7) == 0   &&
      gnrank(sysmi'*sysmi-I,atol=1.e-7) == 0  &&
      isproper(sysmi) && isproper(sysni) &&
      isstable(sysmi,atol=1.e-7) &&
      isstable(sysni,atol=1.e-7)  &&
      isempty(gzero(gminreal([sysni;sysmi],atol=1.e-7)))

# discrete, descriptor, proper
sys = rdss(n,p,m,T = Ty, disc=true,id=ones(Int,3));
@time sysni, sysmi = grcfid(sys, fast = fast,mindeg = true, mininf = true, atol=1.e-7,atol3=1.e-7);
@test gnrank(sys*sysmi-sysni,atol=1.e-7) == 0   &&
      gnrank(sysmi'*sysmi-I,atol=1.e-7) == 0  &&
      isproper(sysmi) && isproper(sysni) &&
      isstable(sysmi,atol=1.e-7) &&
      isstable(sysni,atol=1.e-7)  &&
      isempty(gzero(gminreal([sysni;sysmi],atol=1.e-7)))

# continuous, descriptor, infinite poles -> error
sys = rdss(n,p,m,T = Ty,disc=false,id=[3*ones(Int,1);2*ones(Int,1)]);
try
   @time sysni, sysmi = grcfid(sys, fast = fast, mindeg = true, mininf = true, atol=1.e-7,atol3=1.e-7);
   @test false
catch
   @test true
end

# discrete, descriptor, infinite poles
sys = rdss(n,p,m,T = Ty, disc=true,id=[3*ones(Int,1);2*ones(Int,1)]);
@time sysni, sysmi = grcfid(sys, fast = fast,mindeg = true, mininf = true, atol=1.e-7,atol3=1.e-7);
@test gnrank(sys*sysmi-sysni,atol=1.e-7) == 0   &&
      gnrank(sysmi'*sysmi-I,atol=1.e-7) == 0  &&
      isproper(sysmi) && isproper(sysni) &&
      isstable(sysmi,atol=1.e-7) &&
      isstable(sysni,atol=1.e-7)  &&
      isempty(gzero(gminreal([sysni;sysmi],atol=1.e-7)))

# continuous, descriptor, proper, uncontrollable infinite eigenvalues
sys = rdss(n,p,m,T = Ty,disc=false,iduc=[3*ones(Int,1);2*ones(Int,1)]);
@time sysni, sysmi = grcfid(sys, fast = fast, mindeg = true, mininf = true, atol=1.e-7,atol3=1.e-7);
@test gnrank(sys*sysmi-sysni,atol=1.e-7) == 0   &&
      gnrank(sysmi'*sysmi-I,atol=1.e-7) == 0  &&
      isproper(sysmi) && isproper(sysni) &&
      isstable(sysmi,atol=1.e-7) &&
      isstable(sysni,atol=1.e-7)  &&
      isempty(gzero(gminreal([sysni;sysmi],atol=1.e-7)))

# discrete, descriptor, proper, uncontrollable infinite eigenvalues
sys = rdss(n,p,m,T = Ty,disc=true,iduc=[3*ones(Int,1);2*ones(Int,1)]);
@time sysni, sysmi = grcfid(sys, fast = fast, mindeg = true, mininf = true, atol=1.e-7,atol3=1.e-7);
@test gnrank(sys*sysmi-sysni,atol=1.e-7) == 0   &&
      gnrank(sysmi'*sysmi-I,atol=1.e-7) == 0  &&
      isproper(sysmi) && isproper(sysni) &&
      isstable(sysmi,atol=1.e-7) &&
      isstable(sysni,atol=1.e-7)  &&
      isempty(gzero(gminreal([sysni;sysmi],atol=1.e-7)))

# continuous, descriptor, proper, uncontrollable finite eigenvalues
sys = rdss(n,p,m,T = Ty,disc=false, nfuc = 5);
@time sysni, sysmi = grcfid(sys, fast = fast, mindeg = true, mininf = true, atol=1.e-7,atol3=1.e-7);
@test gnrank(sys*sysmi-sysni,atol=1.e-7) == 0   &&
      gnrank(sysmi'*sysmi-I,atol=1.e-7) == 0  &&
      isproper(sysmi) && isproper(sysni) &&
      isstable(sysmi,atol=1.e-7) &&
      isstable(sysni,atol=1.e-7)  &&
      isempty(gzero(gminreal([sysni;sysmi],atol=1.e-7)))

# discrete, descriptor, proper, uncontrollable finite eigenvalues
sys = rdss(n,p,m,T = Ty,disc=true, nfuc = 5);
@time sysni, sysmi = grcfid(sys, fast = fast, mindeg = true, mininf = true, atol=1.e-7,atol3=1.e-7);
@test gnrank(sys*sysmi-sysni,atol=1.e-7) == 0   &&
      gnrank(sysmi'*sysmi-I,atol=1.e-7) == 0  &&
      isproper(sysmi) && isproper(sysni) &&
      isstable(sysmi,atol=1.e-7) &&
      isstable(sysni,atol=1.e-7)  &&
      isempty(gzero(gminreal([sysni;sysmi],atol=1.e-7)))

end
end

end # grcfid

@testset "glcfid" begin

sys = rdss(0,0,0);
@time sysni, sysmi = glcfid(sys);
@test gnrank(sysmi*sys-sysni) == 0   &&   #  Mi(s)*G(s)-Ni(s) = 0
      gnrank(sysmi*sysmi'-I) == 0  && # Mi(s)*conj(Mi(s))-I = 0
      isproper(sysmi) & isproper(sysni) && # checking properness of factors
      isstable(sysmi,atol=1.e-7) &&
      isstable(sysni,atol=1.e-7)  &&
      isempty(gzero(gminreal([sysni sysmi])))

fast = true; Ty = Complex{Float64}; #Ty = Float64
n = 5; p = 3; m = 2;
for fast in (true, false)
# random examples
for Ty in (Float64, Complex{Float64})

# continuous, standard
sys = rss(n,p,m,T = Ty,disc=false);
@time sysni, sysmi = glcfid(sys, fast = fast,atol=1.e-7,atol3=1.e-7);
@test gnrank(sysmi*sys-sysni,atol=1.e-7) == 0   &&
      iszero(sysmi'-inv(sysmi),atol=1.e-7)   &&
      isproper(sysmi) && isproper(sysni) &&
      isstable(sysmi,atol=1.e-7) &&
      isstable(sysni,atol=1.e-7)  &&
      isempty(gzero(gminreal([sysni sysmi],atol=1.e-7),atol=1.e-7))

# discrete, standard
sys = rss(n,p,m,T = Ty,disc=true);
@time sysni, sysmi = glcfid(sys, fast = fast,atol=1.e-7,atol3=1.e-7);
@test gnrank(sysmi*sys-sysni,atol=1.e-7) == 0   &&
      iszero(sysmi'-inv(sysmi),atol=1.e-7)   &&
      isproper(sysmi) && isproper(sysni) &&
      isstable(sysmi,atol=1.e-7) &&
      isstable(sysni,atol=1.e-7)  &&
      isempty(gzero(gminreal([sysni sysmi],atol=1.e-7),atol=1.e-7))

# continuous, standard, unobservable
sys = rss(n,p,m,T = Ty,disc=false, nuo = 5);
@time sysni, sysmi = glcfid(sys, fast = fast,atol=1.e-7,atol3=1.e-7);
@test gnrank(sysmi*sys-sysni,atol=1.e-7) == 0   &&
      iszero(sysmi'-inv(sysmi),atol=1.e-7)   &&
      isproper(sysmi) && isproper(sysni) &&
      isstable(sysmi,atol=1.e-7) &&
      isstable(sysni,atol=1.e-7)  &&
      isempty(gzero(gminreal([sysni sysmi],atol=1.e-7),atol=1.e-7))

# discrete, standard, unobservable
sys = rss(n,p,m,T = Ty,disc=true, nuo = 5);
@time sysni, sysmi = glcfid(sys, fast = fast,atol=1.e-7,atol3=1.e-7);
@test gnrank(sysmi*sys-sysni,atol=1.e-7) == 0   &&
      iszero(sysmi'-inv(sysmi),atol=1.e-7)   &&
      isproper(sysmi) && isproper(sysni) &&
      isstable(sysmi,atol=1.e-7) &&
      isstable(sysni,atol=1.e-7)  &&
      isempty(gzero(gminreal([sysni sysmi],atol=1.e-7),atol=1.e-7))

# continuous, descriptor, proper, unobservable infinite eigenvalues
sys = rdss(n,p,m,T = Ty,disc=false,iduo=[3*ones(Int,1);2*ones(Int,1)],randlt=true,randrt=false);
#sys = rdss(n,p,m,T = Ty,disc=false,randlt=true,randrt=false);
@time sysni, sysmi = glcfid(sys, fast = fast, mindeg = true, mininf = true, atol=1.e-7,atol3=1.e-7);
#@time sysni, sysmi = glcfid(sys, fast = false, mindeg = false, mininf = false, atol=1.e-7,atol3=1.e-7);
@test gnrank(sysmi*sys-sysni,atol=1.e-7) == 0   &&
      gnrank(sysmi'*sysmi-I,atol=1.e-7) == 0  &&
      isproper(sysmi) && isproper(sysni) &&
      isstable(sysmi,atol=1.e-7) &&
      isstable(sysni,atol=1.e-7)  &&
      isempty(gzero(gminreal([sysni sysmi])))

# discrete, descriptor, proper, unobservable infinite eigenvalues
sys = rdss(n,p,m,T = Ty,disc=true,iduo=[3*ones(Int,1);2*ones(Int,1)]);
@time sysni, sysmi = glcfid(sys, fast = fast, mindeg = true, mininf = true, atol=1.e-7,atol3=1.e-7);
@test gnrank(sysmi*sys-sysni,atol=1.e-7) == 0   &&
      gnrank(sysmi'*sysmi-I,atol=1.e-7) == 0  &&
      isproper(sysmi) && isproper(sysni) &&
      isstable(sysmi,atol=1.e-7) &&
      isstable(sysni,atol=1.e-7)  &&
      isempty(gzero(gminreal([sysni sysmi])))

# continuous, descriptor, proper, unobservable finite eigenvalues
sys = rdss(n,p,m,T = Ty,disc=false, nfuo = 5);
@time sysni, sysmi = glcfid(sys, fast = fast, mindeg = true, mininf = true, atol=1.e-7,atol3=1.e-7);
@test gnrank(sysmi*sys-sysni,atol=1.e-7) == 0   &&
      gnrank(sysmi'*sysmi-I,atol=1.e-7) == 0  &&
      isproper(sysmi) && isproper(sysni) &&
      isstable(sysmi,atol=1.e-7) &&
      isstable(sysni,atol=1.e-7)  &&
      isempty(gzero(gminreal([sysni sysmi])))

# discrete, descriptor, proper, unobservable finite eigenvalues
sys = rdss(n,p,m,T = Ty,disc=true, nfuo = 5);
@time sysni, sysmi = glcfid(sys, fast = fast, mindeg = true, mininf = true, atol=1.e-7,atol3=1.e-7);
@test gnrank(sysmi*sys-sysni,atol=1.e-7) == 0   &&
      gnrank(sysmi'*sysmi-I,atol=1.e-7) == 0  &&
      isproper(sysmi) && isproper(sysni) &&
      isstable(sysmi,atol=1.e-7) &&
      isstable(sysni,atol=1.e-7)  &&
      isempty(gzero(gminreal([sysni sysmi])))

end
end

end # glcfid

end #test

end #module


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule CoinUtils_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("CoinUtils")
JLLWrappers.@generate_main_file("CoinUtils", UUID("be027038-0da8-5614-b30d-e42594cb92df"))
end  # module CoinUtils_jll


text:
"""
    arrowhead(target[, action=:fill];
        shaftangle=0,
        headlength=10,
        headangle=pi/8)

Draw an arrow head. The arrowhead length will be the length of the side of
the arrow's head, and the arrowhead angle is the angle between the sloping
side of the arrowhead and the arrow's shaft.

This doesn't use the current linewidth setting (`setline()`), and defaults to 1,
but you can specify another value.
"""
function arrowhead(target, action=:fill;
        shaftangle=0, headlength=10, headangle=pi/8)
    gsave()
    topangle = shaftangle + headangle
    botangle = shaftangle - headangle

    topx = target.x + cos(topangle) * headlength
    topy = target.y + sin(topangle) * headlength
    botx = target.x + cos(botangle) * headlength
    boty = target.y + sin(botangle) * headlength
    poly([Point(topx, topy), target, Point(botx, boty)], action)
    grestore()
end

"""
    arrow(startpoint::Point, endpoint::Point;
        linewidth         = 1.0,
        arrowheadlength   = 10,
        arrowheadangle    = pi/8,
        decoration        = 0.5 or range(),
        decorate          = nothing,
        arrowheadfunction = nothing)

Draw a line between two points and add an arrowhead at the
end. The arrowhead length will be the length of the side of
the arrow's head, and the arrowhead angle is the angle
between the sloping side of the arrowhead and the arrow's
shaft.

Arrows don't use the current linewidth setting
(`setline()`), and defaults to 1, but you can specify
another value. It doesn't need stroking/filling, the shaft
is stroked and the head filled with the current color.

### Decoration

The `decorate` keyword argument accepts a function with zero
arguments that can execute code at one or more locations on
the arrow's shaft. The inherited graphic environment is
centered at each point on the shaft between 0 and 1 given by
scalar or vector `decoration`, and the x-axis is aligned
with the direction of the curve at that point.

### Arrowheads

A triangular arrowhead is drawn by default. But you can pass
a function to the `arrowheadfunction` keyword argument that
accepts three arguments: the shaft end, the arrow head end,
and the shaft angle. Thsi allows you to draw any shape
arrowhead.

### Example

```
function redbluearrow(shaftendpoint, endpoint, shaftangle)
    @layer begin
        sethue("red")
        sidept1 = shaftendpoint  + polar(10, shaftangle + π/2 )
        sidept2 = shaftendpoint  - polar(10, shaftangle + π/2)
        poly([sidept1, endpoint, sidept2], :fill)
        sethue("blue")
        poly([sidept1, endpoint, sidept2], :stroke, close=false)
    end
end

@drawsvg begin
    background("white")
    arrow(O, O + (120, 120),
        linewidth=4,
        arrowheadlength=40,
        arrowheadangle=π/7,
        arrowheadfunction = redbluearrow)

    arrow(O, 100, 3π/2, π,
        linewidth=4,
        arrowheadlength=20,
        clockwise=false,arrowheadfunction=redbluearrow)
end 800 250
```

"""
function arrow(startpoint::Point, endpoint::Point;
        linewidth         = 1.0,
        arrowheadlength   = 10,
        arrowheadangle    = pi/8,
        decoration        = 0.5,
        decorate          = nothing,
        arrowheadfunction = nothing)

    gsave()
    setlinejoin("butt")
    setline(linewidth)

    isapprox(startpoint, endpoint) && throw(error("can't draw arrow between two identical points"))
    shaftlength = distance(startpoint, endpoint)
    shaftangle = atan(startpoint.y - endpoint.y, startpoint.x - endpoint.x)
    arrowheadtopsideangle = shaftangle + arrowheadangle
    # shorten the length so that lines
    # stop before we get to the arrow
    # thus wide shafts won't stick out through the head of the arrow.
    max_undershoot = shaftlength - ((linewidth/2) / tan(arrowheadangle))
    true_arrowheadlength = arrowheadlength * cos(arrowheadangle)
    if true_arrowheadlength < max_undershoot
        ratio = (shaftlength - true_arrowheadlength)/shaftlength
    else
        ratio = max_undershoot/shaftlength
    end
    tox = startpoint.x + (endpoint.x - startpoint.x) * ratio
    toy = startpoint.y + (endpoint.y - startpoint.y) * ratio
    fromx = startpoint.x
    fromy = startpoint.y

    # draw the shaft of the arrow
    newpath()
    line(Point(fromx, fromy), Point(tox, toy), :stroke)

    #  use a user-supplied arrowhead function?
    toppoint = Point(endpoint.x + cos(arrowheadtopsideangle) * arrowheadlength,
                     endpoint.y + sin(arrowheadtopsideangle) * arrowheadlength)

    arrowheadbottomsideangle = shaftangle - arrowheadangle

    bottompoint = Point(endpoint.x + cos(arrowheadbottomsideangle) * arrowheadlength,
                        endpoint.y + sin(arrowheadbottomsideangle) * arrowheadlength)

    if arrowheadfunction === nothing
        poly([toppoint, endpoint, bottompoint], :fill)
    else
        # shaft end, head end, angle
        arrowheadfunction(Point(tox, toy), endpoint, slope(startpoint, endpoint))
    end

    # prepare to add decorations at point along shaft
    if decorate === nothing
    else
        for decpos in decoration
            decpoint = between(startpoint, endpoint, decpos)
            # slope at this point
            slp = slope(startpoint, endpoint)
            @layer begin
                translate(decpoint)
                rotate(slp)
                decorate()
            end
        end
    end
    grestore()
end

"""
    arrow(centerpos::Point, radius, startangle, endangle;
        linewidth          = 1.0,
        arrowheadlength    = 10,
        arrowheadangle     = π/8,
        decoration         = 0.5,
        decorate           = nothing,
        arrowheadfunction  = nothing,
        clockwise          = true)

Draw a curved arrow, an arc centered at `centerpos` starting
at `startangle` and ending at `endangle` with an arrowhead
at the end. Angles are measured clockwise from the positive
x-axis.

Arrows don't use the current linewidth setting
(`setline()`); you can specify the linewidth.

The `decorate` keyword argument accepts a zero-argument
function that can execute code at one or more locations on
the arrow's shaft. The inherited graphic environment is
centered at points on the shaft between 0 and 1 given by
scalar or vector `decoration`, and the x-axis is aligned
with the direction of the curve at that point.

A triangular arrowhead is drawn by default. But you can pass
a function to the `arrowheadfunction` keyword argument that
accepts three arguments: the shaft end, the arrow head end,
and the shaft angle. Thsi allows you to draw any shape
arrowhead.
"""
function arrow(centerpos::Point, radius, startangle, endangle;
        linewidth          = 1.0,
        arrowheadlength    = 10,
        arrowheadangle     = π/8,
        decoration         = 0.5,
        decorate           = nothing,
        arrowheadfunction  = nothing,
        clockwise          = true)

    # don't bother with them if theyre too small
    if isapprox(startangle, endangle, rtol = 0.01)
        return
    end
    # circular arcs needs swapping to avoid Cairo crash
    if isapprox(startangle - endangle, 2π, rtol = 0.01)
        startangle, endangle = endangle, startangle
    end

    gsave()
    setlinejoin("butt")
    setline(linewidth)
    translate(centerpos)
    θ = mod2pi((2π + endangle) - startangle)
    if clockwise != true
        θ = 2π  - θ
    end
    arclength = radius * θ
    startpoint = Point(radius * cos(startangle), radius * sin(startangle))
    endpoint   = Point(radius * cos(endangle),   radius * sin(endangle))

    # shorten the length so that lines stop before we get to the arrow
    # thus wide shafts won't stick out through the head of the arrow.

    max_undershoot = arclength - ((linewidth/2) / tan(arrowheadangle))
    true_arrowheadlength = arrowheadlength * cos(arrowheadangle)
    if true_arrowheadlength < max_undershoot
        ratio = (arclength - true_arrowheadlength)/arclength
    else
        ratio = max_undershoot/arclength
    end

    # draw the arrow shaft
    newpath()
    move(radius * cos(startangle), radius * sin(startangle))
    newarclength = arclength * ratio
    if clockwise == true
        newendangle = startangle + (newarclength/radius)
        arc(O, radius, startangle, newendangle, :stroke)
    else
        newendangle = startangle - (newarclength/radius)
        carc(O, radius, startangle, newendangle, :stroke)
    end
    closepath()

    # draw the arrowhead
    # rotation of head should be based on end of shaft, not end of arrow
    newendpoint = Point(radius * cos(newendangle), radius * sin(newendangle))
    if clockwise == true
        shaftangle = mod2pi(-π/2 + atan(0 - newendpoint.y, 0 - newendpoint.x))
    else
        shaftangle = mod2pi(-3π/2 + atan(0 - newendpoint.y, 0 - newendpoint.x))
    end

    #  use a user-supplied arrowhead function?
    if arrowheadfunction === nothing
        arrowheadoutersideangle = shaftangle + pi - arrowheadangle
        arrowheadinnersideangle = shaftangle + pi + arrowheadangle
        toppoint =    Point(endpoint.x + cos(arrowheadinnersideangle) * arrowheadlength,
                            endpoint.y + sin(arrowheadinnersideangle) * arrowheadlength)

        bottompoint = Point(endpoint.x + cos(arrowheadoutersideangle) * arrowheadlength,
                            endpoint.y + sin(arrowheadoutersideangle) * arrowheadlength)
        poly([toppoint, endpoint, bottompoint], :fill)
    else
        # shaft end, head end, angle
        arrowheadfunction(newendpoint, endpoint, shaftangle)
    end

    # use-suppplied shaft decoration function?
    if decorate === nothing
    else
        gsave()
        for decpos in decoration
            decorationangle = rescale(decpos, 0, 1, startangle, newendangle)
            decorationpoint = Point(radius * cos(decorationangle), radius * sin(decorationangle))
            ptangle = atan(decorationpoint.y, decorationpoint.x)
            @layer begin
                translate(decorationpoint)
                rotate(π/2 + ptangle)
                decorate()
            end
        end
        grestore()
    end
    grestore()
end

"""
    arrow(start::Point, C1::Point, C2::Point, finish::Point, action=:stroke;
        linewidth       = 1.0,
        arrowheadlength = 10,
        arrowheadangle  = pi/8,
        startarrow      = false,
        finisharrow     = true,
        decoration      = 0.5,
        decorate        = nothing
        arrowheadfunction = nothing)

Draw a Bezier curved arrow, from `start` to `finish`, with
control points `C1` and `C2`. Arrow heads can be
added/hidden by changing `startarrow` and `finisharrow`
options.

The `decorate` keyword argument accepts a function that can
execute code at one or more locations on the arrow's shaft.
The inherited graphic environment is centered at each point
on the shaft given by scalar or vector `decoration`, and the
x-axis is aligned with the direction of the curve at that
point.

### Example

This code draws an arrow head that's filled with orange and outlined in green.

```
function myarrowheadfunction(originalendpoint, newendpoint, shaftangle)
    @layer begin
        setline(5)
        translate(newendpoint)
        rotate(shaftangle)
        sethue("orange")
        ngon(O, 20, 3, 0, :fill)
        sethue("green")
        ngon(O, 20, 3, 0, :stroke)
    end
end

@drawsvg begin
    background("white")
    arrow(O, 220, 0, π,
        linewidth=10,
        arrowheadlength=30,
        arrowheadangle=π/7,
        clockwise=true,
        arrowheadfunction = myarrowheadfunction)
end
```

"""
function arrow(start::Point, C1::Point, C2::Point, finish::Point, action=:stroke;
        # optional
        linewidth=1.0,
        arrowheadlength=10,
        arrowheadangle=pi/8,
        arrowheadfill=true,
        startarrow=false,
        finisharrow=true,
        decoration = 0.5,
        decorate = nothing,
        arrowheadfunction = nothing)
    @layer begin
        setline(linewidth)
        overlap = 0.1 # radians of fudginess
        # TODO rewrite all this
        # arrow heads are a pain
        # length of proposed arrow
        bezlength = polyperimeter(beziertopoly(BezierPathSegment(start, C1, C2, finish)))
        if startarrow && arrowheadfill
            # calculate the shorter version
            actualcurvestart = bezier((arrowheadlength * cos(arrowheadangle + overlap))/bezlength, start, C1, C2, finish)
        else
            actualcurvestart = start
        end
        if finisharrow && arrowheadfill
            # calculate the shorter version
            actualcurvefinish = bezier(1 - (arrowheadlength * cos(arrowheadangle + overlap))/bezlength, start, C1, C2, finish)
        else
            actualcurvefinish = finish
        end

        move(actualcurvestart)
        curve(C1, C2, actualcurvefinish)

        do_action(action)

        start_shaftangle  = slope(start, C1)
        finish_shaftangle = slope(C2, finish)

        if arrowheadfunction === nothing
            finisharrow && arrowhead(finish, arrowheadfill == true ? :fill : :stroke, headlength = arrowheadlength, headangle = arrowheadangle, shaftangle = π + finish_shaftangle)
            startarrow && arrowhead(start, arrowheadfill == true ? :fill : :stroke, headlength = arrowheadlength, headangle = arrowheadangle,   shaftangle = start_shaftangle)
        else
            # shaftangle starts at end and continues outwards
            finisharrow && arrowheadfunction(actualcurvefinish, finish, finish_shaftangle)
            startarrow && arrowheadfunction(actualcurvestart, start, start_shaftangle)
        end

        # prepare to add decorations at points along shaft
        if decorate === nothing
        else
            for decpos in decoration
                decpoint = bezier(decpos, start, C1, C2, finish)
                bp = bezier′(decpos, start, C1, C2, finish)
                slp = atan(bp.y, bp.x)
                @layer begin
                    translate(decpoint)
                    rotate(slp)
                    decorate()
                end
            end
        end
    end # layer
end

"""
    arrow(start::Point, finish::Point, height::Vector, action=:stroke;
        keyword arguments...)

Draw a Bézier arrow between `start` and `finish`, with control points defined to fit in
an imaginary box defined by the two supplied `height` values (see `bezierfrompoints()`). If the height values are different signs, the arrow will change direction on its way.

Keyword arguments are the same as [`arrow(pt1, pt2, pt3, pt4)`](@ref).

### Example

```
arrow(pts[1], pts[end], [15, 15],
    decoration = 0.5,
    decorate = () -> text(string(pts[1])))
```

"""
function arrow(start::Point, finish::Point, height::Vector, action=:stroke;
        # optional kwargs
        linewidth=1.0,
        arrowheadlength=10,
        arrowheadangle=pi/8,
        arrowheadfill=true,
        startarrow=false,
        finisharrow=true,
        decoration = 0.5,
        decorate = nothing,
        arrowheadfunction = nothing)

    @layer begin
        s = slope(start, finish)
        perp1 = start +  polar(height[1], s - π/2)
        perp2 = finish + polar(height[2], s - π/2)
    end

    cpt1 = between(perp1, perp2, 0.33)
    cpt2 = between(perp1, perp2, 0.66)
    pts = bezierfrompoints([start,
        cpt1,
        cpt2,
        finish])

    # i've forgotten how to do this better...
    arrow(pts[1], pts[2], pts[3], pts[4], action,
        linewidth=linewidth,
        arrowheadlength=arrowheadlength,
        arrowheadangle=arrowheadangle,
        arrowheadfill=arrowheadfill,
        startarrow=startarrow,
        finisharrow=finisharrow,
        decoration = decoration,
        decorate = decorate,
        arrowheadfunction = arrowheadfunction)
end

"""
    dimension(p1::Point, p2::Point;
        format::Function   = (d) -> string(d), # process the measured value into a string
        offset             = 0.0,              # left/right, parallel with x axis
        fromextension      = (10.0, 10.0),     # length of extensions lines left and right
        toextension        = (10.0, 10.0),     #
        textverticaloffset = 0.0,              # range 1.0 (top) to -1.0 (bottom)
        texthorizontaloffset = 0.0,            # range 1.0 (top) to -1.0 (bottom)
        textgap            = 5,                # gap between start of each arrow (≈ fontsize?)
        textrotation       = 0.0,
        arrowlinewidth     = 1.0,
        arrowheadlength    = 10,
        arrowheadangle     = π/8)

Calculate and draw dimensioning graphics for the distance between `p1` and `p2`.
The value can be formatted with function `format`.

`p1` is the lower on the page (ie probably the higher y value) point, `p2` is
the higher on the page (ie probably lower y) point.

`offset` is to the left (-x) when negative.

Dimension graphics will be rotated to align with a line between `p1` and `p2`.

In `textverticaloffset`, "vertical" and "horizontal" are best understood by "looking" along the line from the first point to the
second. `textverticaloffset` ranges from -1 to 1, `texthorizontaloffset` in default units.

```
        toextension
        [5  ,  5]
       <---> <--->
                             to
       -----------            +
            ^
            |

           -50

            |
            v
       ----------            +
                            from
       <---> <--->
         [5 , 5]
       fromextension

            <---------------->
                  offset
```

Returns the measured distance and the text.
"""
function dimension(p1::Point, p2::Point;
    format::Function   = (d) -> string(d), # process the measured value into a string
    offset             = 0.0,              # left/right, parallel with x axis
    fromextension      = (10.0, 10.0),     # length of extensions lines left and right
    toextension        = (10.0, 10.0),     #
    textverticaloffset = 0.0,              # range 1.0 (top) to -1.0 (bottom)
    texthorizontaloffset = 0.0,            #
    textgap            = 5,                # gap between start of each arrow (≈ fontsize?)
    textrotation       = 0.0,
    arrowlinewidth     = 1.0,
    arrowheadlength    = 10,
    arrowheadangle     = π/8)

    p1 == p2 && throw(error("dimension(): the two points are the same, can't draw dimensions"))
    d     = distance(p1, p2)
    midpt = midpoint(p1, p2)
    s     = slope(p1, p2)
    t     = format(d)
    @layer begin
        o = perpendicular(midpt, p2, offset)
        translate(o)
        rotate(π/2 + s)

        line(Point(-fromextension[1], d/2), Point(fromextension[2], d/2), :stroke)
        line(Point(-toextension[1], -d/2), Point(toextension[2], -d/2), :stroke)

        # top (lower y)
        apos = O + (0, -textgap - rescale(textverticaloffset, 0, 1, 0, d/2))
        arrow(apos, Point(0, -d/2), linewidth = arrowlinewidth, arrowheadlength = arrowheadlength, arrowheadangle = arrowheadangle)

        # bottom (higher y)
        bpos = O + (0, textgap + rescale(-textverticaloffset, 0, 1, 0, d/2))
        arrow(bpos, Point(0, d/2), linewidth = arrowlinewidth, arrowheadlength = arrowheadlength, arrowheadangle = arrowheadangle)

        te = textextents(t)
        tpos = O + (texthorizontaloffset, -textverticaloffset * d/2)
        @layer begin
            translate(tpos)
            rotate(textrotation)
            text(t, halign=:center, valign=:middle)
        end
    end
    return (d, t)
end


text:
using LinearAlgebraicRepresentation, ViewerGL
GL = ViewerGL; Lar = LinearAlgebraicRepresentation

function randomlines(n=300, t=0.4)
	V = zeros(Float64,2,2*n)
	EV = [zeros(Int64,2) for k=1:n]
	for k=1:n
		v1 = rand(Float64,2)
		v2 = rand(Float64,2)
		vm = (v1+v2)/2
		transl = rand(Float64,2)
		V[:,k] = (v1-vm)*t + transl
		V[:,n+k] = (v2-vm)*t + transl
		EV[k] = [k,n+k]
	end
	return V,EV
end

V,EV = randomlines(300,.4)

meshes = []
for k=1:length(EV)
	color = GL.MayaColors[k%12+1] - (rand(Float64,4)*0.1)
	push!(meshes, GL.GLGrid(V,[EV[k]],color,1) )
end
GL.VIEW(meshes);

meshes = []
model = V,EV
W,EW = Lar.fragmentlines(model);
for k=1:length(EW)
	color = GL.COLORS[k%12+1] - (rand(Float64,4)*0.1)
	push!(meshes, GL.GLGrid(W,[EW[k]],color,1) )
end
GL.VIEW(meshes);


text:
module SIMDPolynomials

export Uninomial, Uniterm, Poly
export PackedMonomial, Monomial, Term, MPoly
export terms, coeffs, content, contprim
export pseudorem

struct Ret{V} <: Function
    value::V
    Ret{V}(value) where {V} = new{V}(value)
    Ret(value) = new{Core.Typeof(value)}(value)
end

(obj::Ret)(args...; kw...) = obj.value

debugmode() = false

include("interface.jl")
include("utils.jl")
include("monomial.jl")
include("packedmonomial.jl")
include("mpoly.jl")
include("poly.jl")

end


text:
###
# Makes the data for series of plots used in the presentation
###
using Distributions, Random, DelimitedFiles
include("../src/JS_SAA_main.jl")

#K = 10,0000 Same set up as above
function genHistogramPlots_multiple(numRuns = 1000)
	Random.seed!(8675309)
	outPath = "../Results/UPDATED_singleKUnifNewsvendor"
	usePoisson=true

	#Gen the data
	K = 10000
	s = .91
	d = 75
	N = 10

	ps = rand(Dirichlet(ones(d)), floor(Int, K/2))
	ps = [ps rand(Dirichlet(.025 * ones(d)), K - floor(Int, K/2)) ]
	@assert size(ps, 2) == K
	@assert size(ps, 1) == d
	p0 = ones(d) ./ d

	alpha_grid = range(0, stop=3N, length=100)
	lams = ones(K)

	cs = JS.getNewsVendorCosts(repeat(collect(1:d), inner=(1, K)), s, K)
	xs = JS.genSSAAtrainers(repeat(collect(1:d), inner=(1, K)), s)

	f = open("$(outPath)_$(K)_$(numRuns).csv", "w")

	#Header
	writedlm(f, ["Run" "Method" "TruePerf" "time" "alpha"], ',')

	#Compute the full-info value once for reference
	t =
	  @elapsed full_info = JS.zstar(xs, cs, ps, lams)
	println(full_info)
	writedlm(f, [1 "FullInfo" full_info t 0.], ',')

	Nhats = fill(N, K)
	for iRun = 1:numRuns
		println("iRun \t:", iRun)
		if usePoisson
			Nhats = rand(Poisson(N), K)
			Nhats[ Nhats.== 0 ] .= 1
		end
		mhats = JS.sim_path(ps, Nhats);

		#for data-driven shrinkage anchor
		phat_avg = JS.get_GM_anchor(mhats)

		#Gen the Oracle const
		t =
		  @elapsed alphaOR, min_indx, or_alpha_curve = JS.oracle_alpha(xs, cs, mhats, ps, lams, p0, alpha_grid)
		writedlm(f, [iRun "Oracle" or_alpha_curve[min_indx] t alphaOR], ',')

		#Gen the LOO cost with the phatAvg shrinkage
		t =
		  @elapsed alphaLOO, min_indx, looUnsc_curve = JS.loo_alpha(xs, cs, mhats, phat_avg, alpha_grid)
		writedlm(f, [iRun "LOO_avg" or_alpha_curve[min_indx] t alphaLOO], ',')

		#The SAA version
		writedlm(f, [iRun "SAA" or_alpha_curve[1] -1 0.0], ',')

		flush(f)
	end #endRun
	close(f)
end

#Generates nice oracle curve used in the presentation slides with a U shape
function genSinglePathCurves(seed=8675309)
	K = 10000
	s = .91
	d = 75
	N = 10
	f = open("../Results/single_path_curves_presentation.csv", "w")

	Random.seed!(seed)

	ps = rand(Dirichlet(ones(d)), floor(Int, K/2))
	ps = [ps rand(Dirichlet(.025 * ones(d)), K - floor(Int, K/2)) ]
	@assert size(ps, 2) == K
	@assert size(ps, 1) == d

	#gen problems
	supp = collect(1:d)
	cs = JS.getNewsVendorCosts(repeat(supp, inner=(1, K)), s, K)
	xs = JS.genSSAAtrainers(repeat(supp, inner=(1, K)), s)
	Nhats = rand(Poisson(N), K)
	Nhats[ Nhats .== 0] .= 1
	lams = ones(K)

	full_info = JS.zstar(xs, cs, ps, lams)

	#gen data
	mhats = JS.sim_path(ps, Nhats)
	p0 = JS.get_GM_anchor(mhats)
	alpha_grid = range(0, stop=3N, length=100)

	#note, this actually computes N * ZLoo
	alphaOR, jstar, outOR =  JS.oracle_alpha(xs, cs, mhats, ps, lams, p0, alpha_grid)
	alphaLOO, jloo, loo_perf = JS.loo_alpha(xs, cs, mhats, p0, alpha_grid)
	loo_perf /= N

	saa_subopt = map(a-> JS.zbar(xs, cs, p0, a, mhats, mhats, lams), alpha_grid) ./ N
	saa_zero = saa_subopt[1]
	saa_subopt = saa_subopt .- saa_zero
	instab = loo_perf - saa_subopt .- saa_zero;

	#Header
	writedlm(f, ["alpha" "OR" "LOO" "SAASubopt" "Instab"], ',')
	writedlm(f, [vec(alpha_grid) outOR loo_perf saa_subopt instab], ',')
	close(f)
end

function genBadSinglePathCurves(seed=8675409)
	K = 10000
	N = 10
	s = .5
	p0 = .3

	#Gen Pks uniformly over interval for now
	ps = rand(K)
	ps = .6 .+ .3 * rand(K)

	f = open("../Results/single_path_curves_bad_presentation.csv", "w")

	Nhats = rand(Poisson(N), K)
	Nhats[Nhats .== 0] .= 1
	mhats = JS.nv_sim_path(ps, Nhats);

	Random.seed!(seed)
	alpha_grid = range(0, stop=2N, length=100)

	#note, this actually computes N * ZLoo
	#loo_perf = map(a->JS.nv_loo(mhats, p0, a, Nhats, s), alpha_grid)
	alphaOR, jstar, outOR = JS.nv_oracle_alpha(mhats, ps, p0, alpha_grid, Nhats, s)
	alphaLOO, jstar, loo_perf = JS.nv_loo_alpha(mhats, p0, alpha_grid, Nhats, s)
	loo_perf /= N

	saa_subopt = map(a-> JS.nv_saa(mhats, p0, a, Nhats, s), alpha_grid) ./ N
	saa_zero = saa_subopt[1]
	saa_subopt = saa_subopt .- saa_zero
	instab = loo_perf - saa_subopt .- saa_zero;

	#Header
	writedlm(f, ["alpha" "OR" "LOO" "SAASubopt" "Instab"], ',')
	writedlm(f, [vec(alpha_grid) outOR loo_perf saa_subopt instab], ',')
	close(f)
end

##
#genHistogramPlots_multiple()
genBadSinglePathCurves()
#genSinglePathCurves()


text:
include("sound/proof_functions.jl")
include("models/acc/acc.jl")
include("models/car/simple_car.jl")
include("models/single_pendulum/single_pend.jl")
include("models/tora/tora.jl")
include("MIP/src/overt_to_mip.jl")
include("MIP/src/mip_utils.jl")
using JLD2
using LaTeXStrings
file_dir = "sound"
file_dir = join(split(@__FILE__, "/")[1:end-1], "/") # get directory of this script
println("Writing results to $file_dir directory.")
ϵ = 0.1
δ = 0.001

function sincos(ϵ, δ; N=-1)
    dyn = :(sin(cos(x)))
    x = [0, π]
    domain = Dict(zip([:x], [x]))
    oa = overapprox(dyn, domain, N=N, ϵ=ϵ)
    t = time()
    result = check_overapprox(oa, domain, [:x], "sincosexample", jobs=2, delta_sat=δ)
    Δt = time() - t
    return result, Δt
end

function sincosysq(ϵ, δ; N=-1, jobs=4)
    dyn = :(sin(cos(x + y^2)))
    x_domain = [0, π]
    y_domain = [-0.25, 1]
    domain = Dict(zip([:x, :y], [x_domain, y_domain]))
    oa = overapprox(dyn, domain, N=N, ϵ=ϵ)
    t = time()
    result = check_overapprox(oa, domain, [:x, :y], "sincos_ysq_example", jobs=jobs, delta_sat=δ)
    Δt = time() - t

    # write to file
    file = open(file_dir*"/sinccosysq_time.txt", "w")
    write(file, "Time to check validity of overapproximation using dreal: ϵ= "*string(ϵ)*", δ="*string(δ)*"\n")
    write(file, "sincosysq check time (sec): "*string(Δt)*"\n")
    write(file, "jobs = "*string(jobs)*"\n")
    write(file, "N=$N\n")
    write(file, "All unsat? $result\n")
    close(file)

    return result, Δt
end

function xy(ϵ, δ; N=-1, jobs=2)
    dyn = :(x*y)
    x = [-1.5, -1.4] # x = [-1.5, 3.5]
    y = [-0.738, -0.6663] # y = [-1.2, 2.2]
    domain = Dict(zip([:x, :y], [x, y]))
    oa = overapprox(dyn, domain, N=N, ϵ=ϵ)
    t = time()
    result = check_overapprox(oa, domain, [:x, :y], "xy_example", jobs=jobs, delta_sat=δ)
    Δt = time() - t
    return result, Δt
end

function xy_SAT(ϵ, δ; N=-1, jobs=2)
    # expect this to produce bad bounds and return SAT
    dyn = :(x*y)
    x = [-1.5, -1.4] # x = [-1.5, 3.5]
    y = [-0.738, -0.6663] # y = [-1.2, 2.2]
    domain = Dict(zip([:x, :y], [x, y]))
    oa = overapprox(dyn, domain, N=N, ϵ=-ϵ) # < note how I am NEGATING ϵ
    t = time()
    result = check_overapprox(oa, domain, [:x, :y], "xy_SAT_example", jobs=jobs, delta_sat=δ)
    Δt = time() - t
    return result, Δt
end

function check_sat_xy(δ; jobs=2)
    # expect result to be sat
    ϕ = [:(y == log(x))]
    ϕ̂ = [:(y == x)]
    defs = []
    x = [1, 2.2]
    domain = Dict(zip([:x], [x]))
    sq = SoundnessQuery(ϕ, # ϕ
                        defs,
                        ϕ̂, # definitions for ϕ̂
                        domain)

    result = check("dreal", sq, "expect_sat_soundness_query.smt2", δ=δ, jobs=jobs) # TODO: pass dreal delta
    println("result is: ", result)
    R = occursin("unsat", result)

    R ? println("all checks pass for check_sat_xy!") : println("Some checks fail :( for check_sat_xy")
    return R
end

function check_xy_whole(ϵ, δ; N=-1, jobs=56)
    # not the right way to do this. expect sat.
    expr = :(x*y)
    x_range = [-1.5, -1.4] # x = [-1.5, 3.5]
    y_range = [-0.738, -0.6663] # y = [-1.2, 2.2]
    domain = Dict(zip([:x, :y], [x_range, y_range]))
    oa = overapprox(expr, domain, N=N, ϵ=ϵ)

    ϕ = [:($(oa.output) == $expr)]
    ####
    all_approx_constraints = [oa.approx_eq..., oa.approx_ineq...]
    all_approx_constraints = [((constraint.args[1] == :≤) || (constraint.args[1] == :≦) ) ? :($(constraint.args[2]) <= $(constraint.args[3])) : constraint for constraint in all_approx_constraints]
    all_approx_constraints = [((constraint.args[1] == :≥) || (constraint.args[1] ==  :≧)) ? :($(constraint.args[2]) >= $(constraint.args[3])) : constraint for constraint in all_approx_constraints]
    ####
    # all_dependencies = all_approx_constraints[get_all_dependecies([oa.output], all_approx_constraints, [:x, :y], [])]
    ## ^ should be all constraints :// but isn't...
    ϕ̂_idx = map( x-> oa.output ∈ Expr.(free_symbols(Basic(x))), all_approx_constraints)
    ϕ̂ = all_approx_constraints[ϕ̂_idx]
    defs = all_approx_constraints[.!ϕ̂_idx]

    sq = SoundnessQuery(ϕ, # ϕ
                        defs,
                        ϕ̂, # definitions for ϕ̂
                        domain)

    result = check("dreal", sq, "mult_soundness_query.smt2", δ=δ, jobs=jobs) # TODO: pass dreal delta
    println("result is: ", result)
    R = occursin("unsat", result)

    R ? println("all checks pass for mult_xy!") : println("Some checks fail :( for mult_xy")
    return R
    # why is this sat?? :(( Because phihat isn't constrained to take on any particular value in it's range at any particular time. the correct way to do this is to propagate bounds. And even then only 1 nesting level of the overapprox can be checked at once without teh use of quantifiers.
end

# slows down...
function xcosy(ϵ, δ; N=-1, jobs=56)
    dyn = :(x*cos(y))
    x = [-1., 1.] #[-1.5, -1.4]
    y = [2., 3.0]#[2.3, 2.4]
    domain = Dict(zip([:x, :y], [x, y]))
    oa = overapprox(dyn, domain, N=N, ϵ=ϵ)
    t = time()
    result = check_overapprox(oa, domain, [:x, :y], "xcosy_example", jobs=jobs, delta_sat=δ)
    Δt = time() - t
    return result, Δt
end

# runs on laptop!
function acc(ϵ, δ; jobs=4)
    # Checking acc
    # In ACC, the derivaties of dimensions 3 and 6 are overapproximated
    x_lead = [90,110]
    v_lead = [32, 32.2]
    gamma_lead = [0,0]
    x_ego = [10,11]
    v_ego = [0, 20]#[30, 30.2]
    gamma_ego = [0, 0]
    input_domains = [x_lead, v_lead, gamma_lead, x_ego, v_ego, gamma_ego]
    input_set = Hyperrectangle(low=[i[1] for i in input_domains], high=[i[2] for i in input_domains])
    bounds = find_controller_bound(file_dir*"/../nnet_files/jair/acc_controller.nnet", input_set, Id()) # returns 1D Hyperrectangle (interval)

    domain = Dict(zip(acc_input_vars, input_domains))
    domain[acc_control_vars[1]] = [low(bounds)..., high(bounds)...]
    # construct overapproximation
    oa_dim3 = overapprox(acc_ẋ₃, domain, N=-1, ϵ=ϵ)
    oa_dim6 = overapprox(acc_ẋ₆, domain, N=-1, ϵ=ϵ)

    println("Using $jobs jobs")
    result_dim3 = check_overapprox(oa_dim3, domain, [acc_input_vars..., acc_control_vars...], "acc_dim3", jobs=jobs, delta_sat=δ)
    result_dim6 = check_overapprox(oa_dim6, domain, [acc_input_vars..., acc_control_vars...], "acc_dim6", jobs=jobs, delta_sat=δ)
    return result_dim3 && result_dim6
end

# first several checks run on laptop, then starts to slow down.
# ran for a month on 25 cores on jodhpur and didn't finish :///
function simple_car(ϵ, δ; N=-1, jobs=4)
    posx = [9.5, 9.55]
    posy = [-4.5, -4.45]
    yaw = [2.1, 2.11]
    vel = [1.5, 1.51]
    input_domains =  [posx, posy, yaw, vel]
    input_set = Hyperrectangle(low=[i[1] for i in input_domains], high=[i[2] for i in input_domains])
    c_bounds = find_controller_bound(file_dir*"/../nnet_files/jair/car_smallest_controller.nnet", input_set, Id()) # returns 2D
    println("Found controller bound.")
    domain = Dict(zip(simple_car_input_vars, input_domains))
    domain[simple_car_control_vars[1]] = [low(c_bounds)[1]..., high(c_bounds)[1]...]  # u is 2d
    domain[simple_car_control_vars[2]] = [low(c_bounds)[2]..., high(c_bounds)[2]...]  # u is 2d

    oa_dim1 = overapprox(simple_car_ẋ, domain, N=N, ϵ=ϵ)
    oa_dim2 = overapprox(simple_car_ẏ, domain, N=N, ϵ=ϵ)

    result_dim1 = check_overapprox(oa_dim1, domain, [simple_car_input_vars..., simple_car_control_vars...], "simple_car_dim1", jobs=56, delta_sat=δ)
    println("Finished checking dim 1")
    result_dim2 = check_overapprox(oa_dim2, domain, [simple_car_input_vars..., simple_car_control_vars...], "simple_car_dim2", jobs=56, delta_sat=δ)
    return result_dim1 && result_dim2
end

# runs on laptop!
function single_pendulum(ϵ, δ; jobs=4)
    θ = [1., 1.2]
    θ_dot = [0., 0.2]
    input_domains = [θ, θ_dot]
    input_set = Hyperrectangle(low=[i[1] for i in input_domains], high=[i[2] for i in input_domains])
    # get controller bounds
    bounds = find_controller_bound(file_dir*"/../nnet_files/jair/single_pendulum_small_controller.nnet", input_set, Id()) # returns 1D
    domain = Dict(zip(single_pend_input_vars, input_domains))
    domain[single_pend_control_vars[1]] = [low(bounds)..., high(bounds)...]  # u is 1d
    oa = overapprox(single_pend_θ_doubledot, domain, N=-1, ϵ=ϵ)

    println("Using $jobs jobs")
    result = check_overapprox(oa, domain, [single_pend_input_vars..., single_pend_control_vars...], "single_pend", jobs=jobs, delta_sat=δ)
    return result
end

# runs on laptop!!!
function tora(ϵ, δ; jobs=4)
    x1 = [0.6, 0.7]
    x2 = [-0.7, -0.6]
    x3 = [-0.4, -0.3]
    x4 = [0.5, 0.6]
    input_domains = [x1, x2, x3, x4]
    input_set = Hyperrectangle(low=[i[1] for i in input_domains], high=[i[2] for i in input_domains])
    # get controller bounds
    bounds = find_controller_bound(file_dir*"/../nnet_files/jair/tora_smallest_controller.nnet", input_set, Id()) # returns 1D
    domain = Dict(zip(tora_input_vars, input_domains))
    domain[tora_control_vars[1]] = [low(bounds)..., high(bounds)...]  # u is 1d
    oa = overapprox(tora_dim2, domain, N=-1, ϵ=ϵ)

    println("Using $jobs jobs")
    result = check_overapprox(oa, domain, [tora_input_vars..., tora_control_vars...], "tora", jobs=jobs, delta_sat=δ)
    return result
end

function run_benchmarks_and_time(ϵ, δ; jobs=4)
    t = time()
    acc(ϵ, δ; jobs=jobs)
    acc_Δt = time() - t
    println("acc done")
    #
    # simple_car(ϵ, δ)
    # println("simple car done")
    #
    t = time()
    single_pendulum(ϵ, δ; jobs=jobs)
    sing_pend_Δt = time() - t
    println("single pendulum done")
    #
    t = time()
    tora(ϵ, δ; jobs=jobs)
    tora_Δt = time() - t
    println("tora done")

    # write to file
    file = open(file_dir*"/benchmark_times.txt", "w")
    write(file, "Time to check validity of overapproximation using dreal: ϵ= "*string(ϵ)*", δ="*string(δ)*"\n")
    write(file, "acc check time (sec): "*string(acc_Δt)*"\n")
    write(file, "single pendulum check time (sec): "*string(sing_pend_Δt)*"\n")
    write(file, "tora check time (sec): "*string(tora_Δt)*"\n")
    write(file, "jobs = "*string(jobs)*"\n")
    close(file)
end

function run_sin_cos(ϵ, δ; N=-1)
    result, Δt = sincos(ϵ, δ, N=N)
    file = open(file_dir*"/sincos_time_Nis"*string(N)*".txt", "w")
    write(file, "Time to check validity of overapproximation using dreal: ϵ= "*string(ϵ)*", δ="*string(δ)*"\n")
    write(file, "Result = "*string(result)*"\n")
    write(file, "sincos time (sec): "*string(Δt)*"\n")
    write(file, "N="*string(N))
    close(file)
end

function run_xy(ϵ, δ; N=-1)
    result, Δt = xy(ϵ, δ, N=N)
    file = open(file_dir*"/xy_time_Nis"*string(N)*".txt", "w")
    write(file, "Time to check validity of overapproximation using dreal: ϵ= "*string(ϵ)*", δ="*string(δ)*"\n")
    write(file, "Result = "*string(result)*"\n")
    write(file, "xy time (sec): "*string(Δt)*"\n")
    write(file, "N="*string(N))
    close(file)
end

function run_simple_car(ϵ, δ; N=-1)
    simple_car(ϵ, δ)
    println("simple car done")
end

# run_sin_cos(ϵ, δ, N=1)
# run_xy(ϵ, δ; N=1)
# xy_SAT(ϵ, δ)
# check_sat_xy(δ)
# xcosy(ϵ, δ)
# run_benchmarks_and_time(ϵ, δ; jobs=1)
#run_simple_car(ϵ, δ, N=1)
sincosysq(ϵ, δ; N=1, jobs=1)


text:
"calculates the entropy of a given state"
entropy(state::InfiniteMPS) = [sum([-j^2*2*log(j) for j in diag(convert(Array,tsvd(state.CR[i])[2]))]) for i in 1:length(state)]

"
given a thermal state, you can map it to an mps by fusing the physical legs together
to prepare a gibbs ensemble, you need to evolve this state with H working on both legs
here we return the 'superhamiltonian' (H*id,id*H)
"
function splitham(ham::MPOHamiltonian)
    fusers = [isomorphism(fuse(p*p'),p*p') for p in ham.pspaces]

    idham = Array{Union{Missing,typeof(ham[1,1,1])},3}(missing,ham.period,ham.odim,ham.odim)
    hamid = Array{Union{Missing,typeof(ham[1,1,1])},3}(missing,ham.period,ham.odim,ham.odim)

    for i in 1:ham.period
        for (k,l) in keys(ham,i)
            idt = isomorphism(ham.pspaces[i],ham.pspaces[i])

            @tensor hamid[i,k,l][-1 -2; -3 -4] :=   ham[i,k,l][-1,12,-3,15]*
                                                    idt[16,13]*
                                                    fusers[i][-2,12,13]*
                                                    conj(fusers[i][-4,15,16])

            @tensor idham[i,k,l][-1 -2; -3 -4] :=   ham[i,k,l][-1,13,-3,16]*
                                                    idt[12,15]*
                                                    fusers[i][-2,12,16]*
                                                    conj(fusers[i][-4,15,13])
        end
    end

    return MPOHamiltonian(hamid),MPOHamiltonian(idham)
end

infinite_temperature(ham::MPOHamiltonian) = [permute(isomorphism(Matrix{eltype(ham[1,1,1])},oneunit(sp)*sp,oneunit(sp)*sp),(1,2,4),(3,)) for sp in ham.pspaces]


text:
using BenchmarkTools
using PkgBenchmark
using LinearAlgebra
using StaticArrays
using TrajectoryOptimization
using TrajOptPlots
using Plots
import TrajectoryOptimization.AbstractSolver
const TO = TrajectoryOptimization

const paramspath = joinpath(@__DIR__,"params.json")
const suite = BenchmarkGroup()

function benchmarkable_solve!(solver; samples=10, evals=10)
    Z0 = deepcopy(get_trajectory(solver))
    solver.opts.verbose = false
    b = @benchmarkable begin
        initial_trajectory!($solver,$Z0)
        solve!($solver)
    end samples=samples evals=evals
    return b
end

# ALTRO
const altro = BenchmarkGroup(["constrained"])
altro["double_int"]    = benchmarkable_solve!(ALTROSolver(Problems.DoubleIntegrator()...))
altro["pendulum"]      = benchmarkable_solve!(ALTROSolver(Problems.Pendulum()...))
altro["cartpole"]      = benchmarkable_solve!(ALTROSolver(Problems.Cartpole()...))
altro["acrobot"]       = benchmarkable_solve!(ALTROSolver(Problems.Acrobot()...))
altro["parallel_park"] = benchmarkable_solve!(ALTROSolver(Problems.DubinsCar(:parallel_park)...))
altro["3obs"]          = benchmarkable_solve!(ALTROSolver(Problems.DubinsCar(:three_obstacles)...))
altro["escape"]        = benchmarkable_solve!(ALTROSolver(Problems.DubinsCar(:escape)...,
    infeasible=true, R_inf=0.1))
altro["quadrotor"]     = benchmarkable_solve!(ALTROSolver(Problems.Quadrotor(:zigzag)...))
altro["airplane"]      = benchmarkable_solve!(ALTROSolver(Problems.YakProblems()...))
suite["ALTRO"] = altro

# iLQR
const ilqr = BenchmarkGroup(["unconstrained"])
ilqr["double_int"]    = benchmarkable_solve!(iLQRSolver(Problems.DoubleIntegrator()...))
ilqr["pendulum"]      = benchmarkable_solve!(iLQRSolver(Problems.Pendulum()...))
ilqr["cartpole"]      = benchmarkable_solve!(iLQRSolver(Problems.Cartpole()...))
ilqr["acrobot"]       = benchmarkable_solve!(iLQRSolver(Problems.Acrobot()...))
ilqr["parallel_park"] = benchmarkable_solve!(iLQRSolver(Problems.DubinsCar(:parallel_park)...))
ilqr["quadrotor"]     = benchmarkable_solve!(iLQRSolver(Problems.Quadrotor(:zigzag)...))
ilqr["airplane"]      = benchmarkable_solve!(iLQRSolver(Problems.YakProblems()...))
suite["iLQR"] = ilqr

# DIRCOL
const dircol = BenchmarkGroup(["constrained"])
dircol["double_int"] = benchmarkable_solve!(DIRCOLSolver(Problems.DoubleIntegrator()..., integration=HermiteSimpson))
dircol["pendulum"]   = benchmarkable_solve!(DIRCOLSolver(Problems.Pendulum()..., integration=HermiteSimpson))
dircol["cartpole"]   = benchmarkable_solve!(DIRCOLSolver(Problems.Cartpole()..., integration=HermiteSimpson))
dircol["acrobot"]    = benchmarkable_solve!(DIRCOLSolver(Problems.Acrobot()..., integration=HermiteSimpson))
dircol["parallel_park"]   = benchmarkable_solve!(DIRCOLSolver(Problems.DubinsCar(:parallel_park)..., integration=HermiteSimpson))
dircol["3obs"]       = benchmarkable_solve!(DIRCOLSolver(Problems.DubinsCar(:three_obstacles)..., integration=HermiteSimpson))
dircol["escape"]     = benchmarkable_solve!(DIRCOLSolver(Problems.DubinsCar(:escape)..., integration=HermiteSimpson))
dircol["quadrotor"]  = benchmarkable_solve!(DIRCOLSolver(Problems.Quadrotor(:zigzag)...,integration=HermiteSimpson))
dircol["airplane"]   = benchmarkable_solve!(DIRCOLSolver(Problems.YakProblems()..., integration=HermiteSimpson))
suite["Ipopt"] = dircol

SUITE = suite