text:
mutable struct FullCovMatrix{T<:AbstractFloat, K <: MatTypes} <: CovMatrix
    C::K
    c_1::T
    c_μ::T
end

function FullCovMatrix(::Type{K}, n, μ_w) where K <: MatTypes
    c_1 = 2/((n+1.3)^2 + μ_w)
    c_μ = min(1-c_1,2*(μ_w-2+1/μ_w)/((n+2)^2 +μ_w))
    C = init_matrix(K,n)

    FullCovMatrix(C,c_1,c_μ)
end

cov_matrix(cov_mat::FullCovMatrix) = cov_mat.C
eigvals(cov_mat::FullCovMatrix) = eigen(cov_mat.C).values #FIXME maybe a bit costly

mv_rand(n,C) = rand(MultivariateNormal(zeros(n),C))
mv_rand(n,C::Diagonal) = rand(MultivariateNormal(zeros(n),diag(C)))

function update_p_σ(::Type{FullCovMatrix{T,K}},c_σ, p_σ, μ_w, C, Δm, σ) where {T,K}
    (1-c_σ)*p_σ  + √(c_σ*(2-c_σ)*μ_w) * C^(-1/2) * Δm/σ
end

function update_p_c(::Type{FullCovMatrix{T,K}},c_c, p_c, H_σ, c_σ, μ_w, Δm, σ) where {T,K}

#    @show (H_σ*√(c_σ*(2 - c_σ)*μ_w), Δm/σ)
#    (1-c_c/3)*p_c + H_σ*√(c_σ*(2 - c_σ)*μ_w) * 2.5Δm/σ
    (1-c_c)*p_c + H_σ*√(c_σ*(2 - c_σ)*μ_w) * Δm/σ
end

function update_C(::Type{FullCovMatrix{T,K}}, C::K, c_1, c_μ, H_σ, c_c, p_c, w, σ, x, m, μ) where {T,K}

#    c_1 = 0.75*c_1
#    c_µ = 2*c_µ

    C = (1-c_1-c_µ + (1-H_σ)*c_1*c_c*(2-c_c)) * C +
        c_1 * p_c * p_c' +
        c_μ * ∑( w[i]/σ^2*( (x[i]-m)*(x[i]-m)') for i=1:μ)

#    @show (1-c_1-c_µ + (1-H_σ)*c_1*c_c*(2-c_c)) * C
#    @show c_1 * p_c * p_c'
#    @show c_μ * ∑( w[i]/σ^2*( (x[i]-m)*(x[i]-m)') for i=1:μ)
#    @show p_c

    C = (C+C')/2 #keep symmetric part
    C
end

#this is probably wrong
function update_C(::Type{FullCovMatrix{T,K}}, C::K, c_1, c_μ, H_σ, c_c, p_c, w, σ, x, m, μ) where {T,K <: Diagonal}

    C = (1-c_1-c_µ + (1-H_σ)*c_1*c_c*(2-c_c)) * C +
        c_1 * Diagonal(p_c.^2) +
        c_μ * ∑( w[i]/σ^2*( Diagonal((x[i]-m).^2) ) for i=1:μ)
    C
end

function update_state!(d, state::CMAState{T,CM}, method::CMA) where {T,CM <: FullCovMatrix}

    @unpack w,n,μ_w,c_σ,d_σ,c_c,p_σ,p_c,cov_mat,m,chi_D,σ,t,xs = state
    @unpack C,c_1,c_μ = cov_mat

    state.f_x_previous = state.f_x
    copy!(state.x_previous, state.x)

    λ = method.λ
    μ = floor(Int,λ/2)

    x = [m + σ*mv_rand(n,C) for i=1:λ]

    fx = eval_objfun(T, λ, d, x)

    idx = sortperm(fx)
    x, fx = x[idx], fx[idx]

    m_t = m
    m = ∑( w[i] * x[i] for i=1:μ)

    Δm = m-m_t

#    @show C

    p_σ = update_p_σ(CM, c_σ, p_σ, μ_w, C, Δm, σ)
    H_σ = update_H_σ(CM, p_σ, c_σ, t, n)
    p_c = update_p_c(CM, c_c, p_c, H_σ, c_σ, μ_w, Δm, σ)
    C = update_C(CM, C, c_1, c_μ, H_σ, c_c, p_c, w, σ, x, m, μ)

    σ = update_σ(CM, σ, c_σ, d_σ, p_σ, chi_D)

    D = eigvals(cov_mat)

#    @printf("σ: %2.2f H_σ: %2.2f axis-ratio: %2.2e \n",
#            σ, H_σ, maximum(D) / minimum(D) )

    t += 1
    state.x = x[1]
    state.f_x = fx[1]
    xs = x

    @pack! cov_mat = C,c_1,c_μ
    @pack! state = w,n,μ_w,c_σ,d_σ,c_c,p_σ,p_c,cov_mat,m,chi_D,σ,t,xs

    false # should the procedure force quit?
end


text:
# Operator hierarchy
const AssignmentOp  = 1
const ConditionalOp = 2
const ArrowOp       = 3
const LazyOrOp      = 4
const LazyAndOp     = 5
const ComparisonOp  = 6
const PipeOp        = 7
const ColonOp       = 8
const PlusOp        = 9
const BitShiftOp    = 10
const TimesOp       = 11
const RationalOp    = 12
const PowerOp       = 13
const DeclarationOp = 14
const WhereOp       = 15
const DotOp         = 16
const PrimeOp       = 16
const DddotOp       = 7
const AnonFuncOp    = 14

precedence(op::Int) = op < Tokens.end_assignments ?  1 :
                       op < Tokens.end_conditional ? 2 :
                       op < Tokens.end_arrow ?       3 :
                       op < Tokens.end_lazyor ?      4 :
                       op < Tokens.end_lazyand ?     5 :
                       op < Tokens.end_comparison ?  6 :
                       op < Tokens.end_pipe ?        7 :
                       op < Tokens.end_colon ?       8 :
                       op < Tokens.end_plus ?        9 :
                       op < Tokens.end_bitshifts ?   10 :
                       op < Tokens.end_times ?       11 :
                       op < Tfokens.end_rational ?   12 :
                       op < Tokens.end_power ?       13 :
                       op < Tokens.end_decl ?        14 :
                       op < Tokens.end_where ?       15 : 16

precedence(op::Token) = op.kind == Tokens.DDDOT ? DddotOp :
                        op.kind < Tokens.begin_assignments ? 0 :
                        op.kind < Tokens.end_assignments ?   1 :
                       op.kind < Tokens.end_conditional ?    2 :
                       op.kind < Tokens.end_arrow ?          3 :
                       op.kind < Tokens.end_lazyor ?         4 :
                       op.kind < Tokens.end_lazyand ?        5 :
                       op.kind < Tokens.end_comparison ?     6 :
                       op.kind < Tokens.end_pipe ?           7 :
                       op.kind < Tokens.end_colon ?          8 :
                       op.kind < Tokens.end_plus ?           9 :
                       op.kind < Tokens.end_bitshifts ?      10 :
                       op.kind < Tokens.end_times ?          11 :
                       op.kind < Tokens.end_rational ?       12 :
                       op.kind < Tokens.end_power ?          13 :
                       op.kind < Tokens.end_decl ?           14 :
                       op.kind < Tokens.end_where ?          15 :
                       op.kind < Tokens.end_dot ?            16 :
                       op.kind == Tokens.ANON_FUNC ? AnonFuncOp :
                       op.kind == Tokens.PRIME ?             16 : 20

precedence(x) = 0
precedence(x::EXPR{OPERATOR{P,K,dot}}) where {P, K, dot} = P

isoperator(kind) = Tokens.begin_ops < kind < Tokens.end_ops
isoperator(t::Token) = isoperator(t.kind)

isunaryop(op::EXPR{OPERATOR{P,K,D}}) where {P, K, D} = isunaryop(K)
isunaryop(t::Token) = isunaryop(t.kind)

isunaryop(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

isunaryandbinaryop(t::Token) = isunaryandbinaryop(t.kind)
isunaryandbinaryop(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

isbinaryop(op::EXPR{OPERATOR{P,K,D}}) where {P, K, D} = isbinaryop(K)
isbinaryop(t::Token) = isbinaryop(t.kind)
isbinaryop(kind) = isoperator(kind) &&
                    !(kind == Tokens.SQUARE_ROOT ||
                    kind == Tokens.CUBE_ROOT ||
                    kind == Tokens.QUAD_ROOT ||
                    kind == Tokens.NOT ||
                    kind == Tokens.NOT_SIGN)

isassignment(t::Token) = Tokens.begin_assignments < t.kind < Tokens.end_assignments

function non_dotted_op(t::Token)
    k = t.kind
    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{OPERATOR{P,K,dot}}) where {P, K, dot}
    P == 1 && !(K == Tokens.APPROX || K == Tokens.PAIR_ARROW) ||
    K == Tokens.RIGHT_ARROW ||
    P == 4 ||
    P == 5 ||
    K == Tokens.ISSUBTYPE ||
    K == Tokens.ISSUPERTYPE ||
    K == Tokens.COLON ||
    K == Tokens.DECLARATION ||
    K == Tokens.DOT ||
    K == Tokens.DDDOT ||
    K == Tokens.PRIME ||
    K == Tokens.WHERE
end

issyntaxunarycall(op) = false
function issyntaxunarycall(op::EXPR{OPERATOR{P,K,false}}) where {P, K}
    K == Tokens.EX_OR ||
    K == Tokens.AND ||
    K == Tokens.DECLARATION ||
    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{OPERATOR{P,K,dot}}) where {P, K, dot}
    startbyte = ps.nt.startbyte - op.span
    if (op isa EXPR{OPERATOR{PlusOp,Tokens.PLUS,false}} || op isa EXPR{OPERATOR{PlusOp,Tokens.MINUS,false}}) && (ps.nt.kind ==  Tokens.INTEGER || ps.nt.kind == Tokens.FLOAT) && isemptyws(ps.ws)
        next(ps)
        arg = INSTANCE(ps)
        arg.span += op.span
        if op isa EXPR{OPERATOR{PlusOp,Tokens.MINUS,false}}
            arg.val = string("-", arg.val)
        end
        return arg
    end

    # Parsing
    prec = P == DeclarationOp ? DeclarationOp :
                K == Tokens.AND ? 14 :
                K == Tokens.EX_OR ? 20 : 13
    @catcherror ps startbyte arg = @precedence ps prec parse_expression(ps)

    if issyntaxunarycall(op)
        ret = EXPR{UnarySyntaxOpCall}(EXPR[op, arg], op.span + arg.span, Variable[], "")
    else
        ret = EXPR{Call}(EXPR[op, arg], op.span + arg.span, Variable[], "")
    end
    return ret
end

function parse_unary(ps::ParseState, op::EXPR{OPERATOR{ColonOp,Tokens.COLON,false}})
    startbyte = ps.nt.startbyte - op.span
    if Tokens.begin_keywords < ps.nt.kind < Tokens.end_keywords ||
        Tokens.begin_literal < ps.nt.kind < Tokens.end_literal ||
        isoperator(ps.nt.kind) ||
        ps.nt.kind == Tokens.IDENTIFIER
        # Parsing
        next(ps)
        arg = INSTANCE(ps)
        return EXPR{Quotenode}(EXPR[op, arg], op.span + arg.span, Variable[], "")
    elseif closer(ps)
        return op
    else
        # Parsing
        @catcherror ps startbyte arg = @precedence ps 20 parse_expression(ps)
        return EXPR{Quote}(EXPR[op, arg], op.span + arg.span, Variable[], "")
    end
end

# Parse assignments
function parse_operator(ps::ParseState, ret::EXPR, op::EXPR{OPERATOR{AssignmentOp,Tokens.EQ,false}})
    startbyte = ps.nt.startbyte - op.span - ret.span

    is_func_call(ret) && _get_sig_defs!(ret)
    # Parsing
    @catcherror ps startbyte nextarg = @precedence ps AssignmentOp - LtoR(AssignmentOp) parse_expression(ps)

    if is_func_call(ret)
        # Construction
        # NOTE : issue w/ scheme parser
        if length(ret.args) > 1 && !(ret.args[2] isa EXPR{OPERATOR{DeclarationOp,Tokens.DECLARATION,false}}) && ps.closer.precedence != 0
            nextarg = EXPR{Block}(EXPR[nextarg], nextarg.span, Variable[], "")
        end

        ret1 = EXPR{BinarySyntaxOpCall}(EXPR[ret, op, nextarg], op.span + ret.span + nextarg.span, Variable[], "")
        ret1.defs = [Variable(Expr(_get_fname(ret)), :Function, ret1)]
        return ret1
    else
        defs = ps.trackscope ? _track_assignment(ps, ret, nextarg) : Variable[]
        if nextarg isa EXPR{BinarySyntaxOpCall} && nextarg.args[2] isa EXPR{OPERATOR{AssignmentOp,Tokens.EQ,false}}
            append!(defs, nextarg.defs)
            empty!(nextarg.defs)
        end
        ret = EXPR{BinarySyntaxOpCall}(EXPR[ret, op, nextarg], op.span + ret.span + nextarg.span, Variable[], "")
        ret.defs = defs
        return ret
    end
end

# Parse conditionals
function parse_operator(ps::ParseState, ret::EXPR, op::EXPR{OPERATOR{ConditionalOp,Tokens.CONDITIONAL,false}})
    startbyte = ps.nt.startbyte - op.span - ret.span

    # Parsing
    @catcherror ps startbyte nextarg = @closer ps ifop parse_expression(ps)
    @catcherror ps startbyte op2 = INSTANCE(next(ps))
    @catcherror ps startbyte nextarg2 = @closer ps comma @precedence ps 0 parse_expression(ps)

    # Construction
    ret = EXPR{ConditionalOpCall}(EXPR[ret, op, nextarg, op2, nextarg2], ps.nt.startbyte - startbyte, Variable[], "")
    return ret
end

# Parse comparisons
function parse_operator(ps::ParseState, ret::EXPR, op::EXPR{OPERATOR{ComparisonOp,K,dot}}) where {K, dot}
    startbyte = ps.nt.startbyte - op.span - ret.span

    # Parsing
    @catcherror ps startbyte nextarg = @precedence ps ComparisonOp - LtoR(ComparisonOp) parse_expression(ps)

    # Construction
    if ret isa EXPR{Comparison}
        push!(ret.args, op)
        push!(ret.args, nextarg)
        ret.span += op.span + nextarg.span
    elseif ret isa EXPR{BinaryOpCall} && precedence(ret.args[2]) == ComparisonOp
        ret = EXPR{Comparison}(EXPR[ret.args[1], ret.args[2], ret.args[3], op, nextarg], ret.span + op.span + nextarg.span, Variable[], "")
    elseif ret isa EXPR{BinarySyntaxOpCall} && (ret.args[2] isa EXPR{OPERATOR{ComparisonOp,Tokens.ISSUBTYPE,false}} || ret.args[2] isa EXPR{OPERATOR{ComparisonOp,Tokens.ISSUPERTYPE,false}})
        ret = EXPR{Comparison}(EXPR[ret.args[1], ret.args[2], ret.args[3], op, nextarg], ret.span + op.span + nextarg.span, Variable[], "")
    elseif (op isa EXPR{OPERATOR{ComparisonOp,Tokens.ISSUBTYPE,false}} || op isa EXPR{OPERATOR{ComparisonOp,Tokens.ISSUPERTYPE,false}})
        ret = EXPR{BinarySyntaxOpCall}(EXPR[ret, op, nextarg], ret.span + op.span + nextarg.span, Variable[], "")
    else
        ret = EXPR{BinaryOpCall}(EXPR[ret, op, nextarg], op.span + ret.span + nextarg.span, Variable[], "")
    end
    return ret
end

# Parse ranges
function parse_operator(ps::ParseState, ret::EXPR, op::EXPR{OPERATOR{ColonOp,Tokens.COLON,false}})
    startbyte = ps.t.startbyte

    # Parsing
    @catcherror ps startbyte - op.span - ret.span nextarg = @precedence ps ColonOp - LtoR(ColonOp) parse_expression(ps)

    # Construction
    if ret isa EXPR{BinarySyntaxOpCall} && ret.args[2] isa EXPR{OPERATOR{ColonOp,Tokens.COLON,false}} && length(ret.args) == 3
        ret = EXPR{ColonOpCall}(ret.args, ret.span, Variable[], "")
        push!(ret.args, op)
        push!(ret.args, nextarg)
        ret.span += ps.nt.startbyte - startbyte
    else
        ret = EXPR{BinarySyntaxOpCall}(EXPR[ret, op, nextarg], ret.span + ps.nt.startbyte - startbyte, Variable[], "")
    end
    return ret
end

parse_operator(ps::ParseState, ret::EXPR{BinaryOpCall}, op::EXPR{OPERATOR{PlusOp,Tokens.PLUS,false}}) = parse_chain_operator(ps, ret, op)
parse_operator(ps::ParseState, ret::EXPR{ChainOpCall}, op::EXPR{OPERATOR{PlusOp,Tokens.PLUS,false}}) = parse_chain_operator(ps, ret, op)

parse_operator(ps::ParseState, ret::EXPR{BinaryOpCall}, op::EXPR{OPERATOR{TimesOp,Tokens.STAR,false}}) = parse_chain_operator(ps, ret, op)
parse_operator(ps::ParseState, ret::EXPR{ChainOpCall}, op::EXPR{OPERATOR{TimesOp,Tokens.STAR,false}}) = parse_chain_operator(ps, ret, op)

function parse_chain_operator(ps::ParseState, ret::EXPR{ChainOpCall}, op::EXPR{OPERATOR{P,K,false}}) where {P, K}
    startbyte = ps.nt.startbyte - op.span - ret.span

    if ret.args[2] isa EXPR{OPERATOR{P,K,false}}
        # Parsing
        @catcherror ps startbyte nextarg = @precedence ps P - LtoR(P) parse_expression(ps)

        # Construction
        push!(ret.args, op)
        push!(ret.args, nextarg)
        ret.span += nextarg.span + op.span
    else
        ret = invoke(parse_operator, Tuple{ParseState,EXPR,EXPR{OPERATOR{P1,K1,dot}} where {P1, K1, dot}}, ps, ret, op)
    end
    return ret
end

function parse_chain_operator(ps::ParseState, ret::EXPR{BinaryOpCall}, op::EXPR{OPERATOR{P,K,false}}) where {P, K}
    startbyte = ps.nt.startbyte - op.span - ret.span

    if ret.args[2] isa EXPR{OPERATOR{P,K,false}} && ret.args[2].span > 0
        @catcherror ps startbyte nextarg = @precedence ps P - LtoR(P) parse_expression(ps)
        ret = EXPR{ChainOpCall}(ret.args, ret.span, Variable[], "")
        push!(ret.args, op)
        push!(ret.args, nextarg)
        ret.span += nextarg.span + op.span
    else
        ret = invoke(parse_operator, Tuple{ParseState,EXPR,EXPR{OPERATOR{P,K,dot}} where {P1, K1, dot}}, ps, ret, op)
    end
    return ret
end

# Parse power (special case for preceding unary ops)
function parse_operator(ps::ParseState, ret::EXPR, op::EXPR{OPERATOR{PowerOp,K,dot}}) where {K, dot}
    startbyte = ps.nt.startbyte - op.span - ret.span

    # Parsing
    @catcherror ps startbyte nextarg = @precedence ps PowerOp - LtoR(PowerOp) @closer ps inwhere parse_expression(ps)

    # Construction
    # NEEDS FIX
    if ret isa EXPR{UnaryOpCall}
        if false
            xx = EXPR{InvisBrackets}([ret], ret.span, Variable[], "")
            nextarg = EXPR{BinaryOpCall}(EXPR[op, xx, nextarg], op.span + xx.span + nextarg.span, Variable[], "")
        else
            nextarg = EXPR{BinaryOpCall}(EXPR[ret.args[2], op, nextarg], op.span + ret.args[2].span + nextarg.span, Variable[], "")
        end
        ret = EXPR{UnaryOpCall}(EXPR[ret.args[1], nextarg], ret.args[1].span + nextarg.span, Variable[], "")
    else
        ret = EXPR{BinaryOpCall}(EXPR[ret, op, nextarg], op.span + ret.span + nextarg.span, Variable[], "")
    end
    return ret
end

# parse where
function parse_operator(ps::ParseState, ret::EXPR, op::EXPR{OPERATOR{WhereOp,Tokens.WHERE,false}})
    startbyte = ps.nt.startbyte - op.span - ret.span

    # Parsing
    ret = EXPR{BinarySyntaxOpCall}(EXPR[ret, op], 0, Variable[], "")
    # Parsing
    if ps.nt.kind == Tokens.LBRACE
        next(ps)
        push!(ret.args, INSTANCE(ps))
        while ps.nt.kind != Tokens.RBRACE
            @catcherror ps startbyte a = @default ps @nocloser ps newline @closer ps comma @closer ps brace parse_expression(ps)
            push!(ret.args, a)
            if ps.nt.kind == Tokens.COMMA
                next(ps)
                push!(ret.args, INSTANCE(ps))
                format_comma(ps)
            end
        end
        next(ps)
        push!(ret.args, INSTANCE(ps))
    else
        @catcherror ps startbyte nextarg = @precedence ps 5 @closer ps inwhere parse_expression(ps)
        push!(ret.args, nextarg)
    end

    # Construction
    ret.span = ps.nt.startbyte - startbyte
    return ret
end

# parse dot access
function parse_operator(ps::ParseState, ret::EXPR, op::EXPR{OPERATOR{DotOp,Tokens.DOT,false}})
    startbyte = ps.nt.startbyte - op.span - ret.span

    # Parsing
    if ps.nt.kind == Tokens.LPAREN
        startbyte1 = ps.nt.startbyte
        @catcherror ps startbyte sig = @default ps @closer ps paren parse_call(ps, ret)
        args = EXPR{TupleH}(sig.args[2:end], ps.nt.startbyte - startbyte1, Variable[], "")
        ret = EXPR{BinarySyntaxOpCall}(EXPR[ret, op, args], ps.nt.startbyte - startbyte, Variable[], "")
        return ret
    elseif iskw(ps.nt) || ps.nt.kind == Tokens.IN || ps.nt.kind == Tokens.ISA || ps.nt.kind == Tokens.WHERE
        next(ps)
        # nextarg = IDENTIFIER(ps.nt.startbyte - ps.t.startbyte, Symbol(lowercase(string(ps.t.kind))))
        nextarg = EXPR{IDENTIFIER}(EXPR[], ps.nt.startbyte - ps.t.startbyte, Variable[], lowercase(string(ps.t.kind)))
    elseif ps.nt.kind == Tokens.COLON
        next(ps)
        op2 = INSTANCE(ps)
        if ps.nt.kind == Tokens.LPAREN
            @catcherror ps startbyte nextarg = @precedence ps DotOp - LtoR(DotOp) parse_expression(ps)
            nextarg = EXPR{Quote}(EXPR[op2, nextarg], op2.span + nextarg.span, Variable[], "")
        else
            @catcherror ps startbyte nextarg = @precedence ps DotOp - LtoR(DotOp) parse_unary(ps, op2)
        end
    elseif ps.nt.kind == Tokens.EX_OR && ps.nnt.kind == Tokens.LPAREN
        next(ps)
        op2 = OPERATOR(ps)
        @catcherror ps startbyte nextarg = parse_call(ps, op2)
    else
        @catcherror ps startbyte nextarg = @precedence ps DotOp - LtoR(DotOp) parse_expression(ps)
    end

    # Construction
    # NEEDS FIX
    if nextarg isa EXPR{IDENTIFIER} || nextarg isa EXPR{Vect} || nextarg isa EXPR{UnarySyntaxOpCall} && nextarg.args[1] isa EXPR{OPERATOR{PlusOp,Tokens.EX_OR,false}}
        ret = EXPR{BinarySyntaxOpCall}(EXPR[ret, op, Quotenode(nextarg)], op.span + ret.span + nextarg.span, Variable[], "")
    elseif nextarg isa EXPR{MacroCall}
        mname = EXPR{BinarySyntaxOpCall}(EXPR[ret, op, Quotenode(nextarg.args[1])], ret.span + op.span + nextarg.args[1].span, Variable[], "")
        ret = EXPR{MacroCall}(EXPR[mname], ret.span + op.span + nextarg.span, Variable[], "")
        for i = 2:length(nextarg.args)
            push!(ret.args, nextarg.args[i])
        end
    else
        ret = EXPR{BinarySyntaxOpCall}(EXPR[ret, op, nextarg], op.span + ret.span + nextarg.span, Variable[], "")
    end
    return ret
end

function parse_operator(ps::ParseState, ret::EXPR, op::EXPR{OPERATOR{DddotOp,Tokens.DDDOT,false}})
    return EXPR{UnarySyntaxOpCall}(EXPR[ret, op], op.span + ret.span, Variable[], "")
end

function parse_operator(ps::ParseState, ret::EXPR, op::EXPR{OPERATOR{16,Tokens.PRIME,dot}}) where dot
    return EXPR{UnarySyntaxOpCall}(EXPR[ret, op], op.span + ret.span, Variable[], "")
end

function parse_operator(ps::ParseState, ret::EXPR, op::EXPR{OPERATOR{P,Tokens.ANON_FUNC,false}}) where {P}
    startbyte = ps.nt.startbyte - op.span - ret.span
    # Parsing
    @catcherror ps startbyte arg = @closer ps comma @precedence ps 0 parse_expression(ps)

    # Construction
    if ret isa EXPR{TupleH}
        for a in ret.args
            if !(a isa EXPR{P1} where P1 <: PUNCTUATION)
                push!(ret.defs, Variable(Expr(get_id(a)), get_t(a), a))
            end
        end
    else
        push!(ret.defs, Variable(Expr(get_id(ret)), get_t(ret), ret))
    end
    ret = EXPR{BinarySyntaxOpCall}(EXPR[ret, op, EXPR{Block}(EXPR[arg], arg.span, Variable[], "")], ps.nt.startbyte - startbyte, Variable[], "")

    return ret
end

function parse_operator(ps::ParseState, ret::EXPR, op::EXPR{OPERATOR{P,K,dot}}) where {P, K, dot}
    startbyte = ps.nt.startbyte - op.span - ret.span

    # Parsing
    @catcherror ps startbyte nextarg = @precedence ps P - LtoR(P) parse_expression(ps)

    # Construction
    if issyntaxcall(op)
        ret = EXPR{BinarySyntaxOpCall}([ret, op, nextarg], op.span + ret.span + nextarg.span, Variable[], "")
    else
        ret = EXPR{BinaryOpCall}([ret, op, nextarg], op.span + ret.span + nextarg.span, Variable[], "")
    end
    return ret
end

ext: jl


text:
# TemplateName
isNull(x::TemplateName) = clang_TemplateName_isNull(x)

getKind(x::TemplateName) = clang_TemplateName_getKind(x)

getUnderlying(x::TemplateName) = TemplateName(clang_TemplateName_getUnderlying(x))

getNameToSubstitute(x::TemplateName) = TemplateName(clang_TemplateName_getNameToSubstitute(x))

isDependent(x::TemplateName) = clang_TemplateName_isDependent(x)

isInstantiationDependent(x::TemplateName) = clang_TemplateName_isInstantiationDependent(x)

containsUnexpandedParameterPack(x::TemplateName) = clang_TemplateName_containsUnexpandedParameterPack(x)

dump(x::TemplateName) = clang_TemplateName_dump(x)

function getTemplateName(x::TemplateSpecializationType)
    return TemplateName(clang_TemplateSpecializationType_getTemplateName(get_type_ptr(x)))
end


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

JLLWrappers.@generate_main_file_header("OpenAL")
JLLWrappers.@generate_main_file("OpenAL", UUID("c52b6589-6b5d-587d-9bb5-adf8a44d3946"))
end  # module OpenAL_jll


text:
module PCloud
using HTTP
using JSON3
using SHA
using AbstractTrees
import AbstractTrees: children, printnode

export authorize!, PCloudClient, authtoken, update_endpoints!, reset_endpoints!, useglobally!, search, gettree!

# Filesystem tree specific functions
export gettree, fileid, folderid, name, metadata

mutable struct PCloudClient
    apiep::String
    binep::String
    rootep::String
    cookies::Dict{String, String}
end

include("api.jl")
using .API

mutable struct PCloudOpts
    client::PCloudClient
end

include("client.jl")

const OPTS = PCloudOpts(PCloudClient(set_endpoints = false))

include("fs.jl")

end # module


text:
#=
Convert canvas coordinates to voltage vs time coordinates
=#

function coordinate_transform(sc::Single_Channel,event)
    (x1,x2,y1,y2)=coordinate_transform(sc,sc.mi[1],sc.mi[2],event.x,event.y)
end

function coordinate_transform(sc::Single_Channel,xi1::Float64,yi1::Float64,xi2::Float64,yi2::Float64)

    ctx=sc.ctx2

    #This should be optimized
    myx=[1.0;collect(2:(sc.wave_points-1)).*(sc.w2/sc.wave_points)]
    x1=indmin(abs.(myx .- xi1))
    x2=indmin(abs.(myx .- xi2))
    s=sc.s
    o=sc.o
    y1=(yi1-sc.h2/2+o)/s
    y2=(yi2-sc.h2/2+o)/s

    #ensure that left most point is first
    if x1>x2
        x=x1
        x1=x2
        x2=x
        y=y1
        y1=y2
        y2=y
    end
    (x1,x2,y1,y2)
end

#=
Rubber Band functions adopted from GtkUtilities.jl package by Tim Holy 2015
=#

function rubberband_start(sc::Single_Channel, x, y, button_num=1)

    sc.rb = RubberBand(Vec2(x,y), Vec2(x,y), Vec2(x,y), [Vec2(x,y)],false, 2)
    sc.selected=falses(500)
    sc.plotted=falses(500)

    if button_num==1
        push!((sc.c2.mouse, :button1motion),  (c, event) -> rubberband_move(sc,event.x, event.y))
        push!((sc.c2.mouse, :motion), Gtk.default_mouse_cb)
        push!((sc.c2.mouse, :button1release), (c, event) -> rubberband_stop(sc,event.x, event.y,button_num))
    elseif button_num==3
        push!((sc.c2.mouse, :motion),  (c, event) -> rubberband_move(sc,event.x, event.y))
        push!((sc.c2.mouse, :button3release), (c, event) -> rubberband_stop(sc,event.x, event.y,button_num))
    end
    sc.rb_active=true
    nothing
end

function rubberband_move(sc::Single_Channel, x, y)

    sc.rb.moved = true
    sc.rb.pos2 = Vec2(x ,y)
    nothing
end

function rubberband_stop(sc::Single_Channel, x, y,button_num)

    if button_num==1
        pop!((sc.c2.mouse, :button1motion))
        pop!((sc.c2.mouse, :motion))
        pop!((sc.c2.mouse, :button1release))
    elseif button_num==3
        pop!((sc.c2.mouse, :motion))
        pop!((sc.c2.mouse, :button3release))
    end

    sc.rb.moved = false
    sc.rb_active=false
    clear_rb(sc)
    nothing
end

function draw_rb(sc::Single_Channel)

    if sc.rb.moved

        if sc.pause_state == 1

            ctx = sc.ctx2
            clear_rb(sc)

            line(ctx,sc.rb.pos0.x,sc.rb.pos2.x,sc.rb.pos0.y,sc.rb.pos2.y)
            set_line_width(ctx,1.0)
            set_source_rgb(ctx,1.0,1.0,1.0)
            stroke(ctx)

            #Find selected waveforms and plot
            if (sc.buf.selected_clus>0)&((sc.buf.count>0)&(sc.pause))
                get_selected_waveforms(sc,sc.buf.spikes)
                mycolor=1
                if sc.click_button==1
                    mycolor=sc.buf.selected_clus+1
                elseif sc.click_button==3
                    mycolor=1
                end
                plot_selected_waveforms(sc,sc.buf.spikes,mycolor)
            end
            sc.rb.pos1=sc.rb.pos2
        elseif sc.pause_state == 2
            draw_template(sc)
        end

    end

    nothing
end

function clear_rb(sc::Single_Channel)

    line(sc.ctx2,sc.rb.pos0.x,sc.rb.pos1.x,sc.rb.pos0.y,sc.rb.pos1.y)
    set_line_width(sc.ctx2,2.0)
    set_source(sc.ctx2,sc.ctx2s)
    stroke(sc.ctx2)

    nothing
end

#=
Canvas Draw Templates
=#

function draw_start(sc::Single_Channel, x, y, temp,button_num=1)

    vec_list=Array(Vec2,0)

    clus = sc.buf.selected_clus
    for i=1:size(temp.templates,1)
        push!(vec_list,Vec2(temp.templates[i,clus]-temp.sig_min[i,clus],temp.templates[i,clus]+temp.sig_max[i,clus]))
    end

    sc.rb = RubberBand(Vec2(x,y), Vec2(x,y), Vec2(x,y),vec_list,false, 2)
    sc.selected=falses(500)
    sc.plotted=falses(500)

    if button_num==1
        push!((sc.c2.mouse, :button1motion),  (c, event) -> draw_move(sc,event.x, event.y))
        push!((sc.c2.mouse, :motion), Gtk.default_mouse_cb)
        push!((sc.c2.mouse, :button1release), (c, event) -> draw_stop(sc,event.x, event.y,temp,button_num))
    elseif button_num==3
        push!((sc.c2.mouse, :motion),  (c, event) -> draw_move(sc,event.x, event.y))
        push!((sc.c2.mouse, :button3release), (c, event) -> draw_stop(sc,event.x, event.y,temp,button_num))
    end
    sc.rb_active=true
    nothing
end

function draw_move(sc::Single_Channel, x, y)

    sc.rb.moved = true
    sc.rb.pos2 = Vec2(x ,y)
    nothing
end

function draw_stop(sc::Single_Channel, x, y,temp,button_num)

    if button_num==1
        pop!((sc.c2.mouse, :button1motion))
        pop!((sc.c2.mouse, :motion))
        pop!((sc.c2.mouse, :button1release))
    elseif button_num==3
        pop!((sc.c2.mouse, :motion))
        pop!((sc.c2.mouse, :button3release))
    end

    sc.rb.moved = false
    sc.rb_active=false

    clus = sc.buf.selected_clus
    for i=1:length(sc.rb.polygon)
        temp.templates[i,clus] = (sc.rb.polygon[i].x + sc.rb.polygon[i].y) / 2
        temp.sig_min[i,clus] = temp.templates[i,clus] - sc.rb.polygon[i].x
        temp.sig_max[i,clus] = sc.rb.polygon[i].y - temp.templates[i,clus]
    end
    nothing
end

function draw_template(sc::Single_Channel)

    (x1,x2,y1,y2)=coordinate_transform(sc,sc.rb.pos2.x,sc.rb.pos2.y,sc.rb.pos2.x,sc.rb.pos2.y)

    if x1 != sc.rb.pos1.x
        sc.rb.pos1 = Vec2(x1,x1)

        myline=0 #0 for bottom, 1 for top

        mymean=(sc.rb.polygon[x1].x+sc.rb.polygon[x1].y)/2

        if y1 < mymean
            myline = 0
        else
            myline = 1
        end

        s=sc.s
        ctx = sc.ctx2
        wave_points=length(sc.rb.polygon)
        Cairo.translate(ctx,0.0,sc.h2/2)
        Gtk.scale(ctx,sc.w2/wave_points,s)
        clus = sc.buf.selected_clus + 1

        if myline == 0
            move_to(ctx,1,sc.rb.polygon[1].x)
            for i=2:wave_points
                line_to(ctx,i,sc.rb.polygon[i].x)
            end
            sc.rb.polygon[x1] = Vec2(y1,sc.rb.polygon[x1].y)
        else
            move_to(ctx,1,sc.rb.polygon[1].y)
            for i=2:wave_points
                line_to(ctx,i,sc.rb.polygon[i].y)
            end
            sc.rb.polygon[x1] = Vec2(sc.rb.polygon[x1].x,y1)
        end

        identity_matrix(ctx)

        set_line_width(ctx,2.0)
        set_source(ctx,sc.ctx2s)
        stroke(ctx)

        Cairo.translate(ctx,0.0,sc.h2/2)
        Gtk.scale(ctx,sc.w2/wave_points,s)

        if myline == 0
            move_to(ctx,1,sc.rb.polygon[1].x)
            for i=2:wave_points
                line_to(ctx,i,sc.rb.polygon[i].x)
            end
        else
            move_to(ctx,1,sc.rb.polygon[1].y)
            for i=2:wave_points
                line_to(ctx,i,sc.rb.polygon[i].y)
            end
        end

        set_line_width(ctx,1.0)
        select_color(ctx,clus)
        stroke(ctx)

        identity_matrix(ctx)

    end

    nothing
end

#=
Find waveforms that cross line defined by (x1,y1),(x2,y2)
=#
function find_intersected_waveforms(input::Array{T,2},mask,count,x1,y1,x2,y2) where T

    #Bounds check
    x1 = x1 < 2 ? 2 : x1
    x2 = x2 > size(input,2)-2 ? size(input,2)-2 : x2

    for i=1:count
        for j=(x1-1):(x2+1)
            if SpikeSorting.intersect(x1,x2,j,j+1,y1,y2,input[j,i],input[j+1,i])
                mask[i]=false
                break
            end
        end
    end

    nothing
end

#=
Find which of the intersected waveforms are in a different cluster and if that difference has already been plotted
=#
function get_selected_waveforms(sc::Single_Channel,input::Array{T,2}) where T<:Real

    (x1,x2,y1,y2)=coordinate_transform(sc,sc.rb.pos0.x,sc.rb.pos0.y,sc.rb.pos2.x,sc.rb.pos2.y)

    intersection = trues(sc.buf.ind)
    find_intersected_waveforms(sc.buf.spikes,intersection,sc.buf.ind,x1,y1,x2,y2)

    for i=1:sc.buf.ind
        if sc.buf.mask[i]
            if !intersection[i]
                sc.selected[i]=true
            end
            if (sc.plotted[i])&(intersection[i])
                sc.selected[i]=false
            end
        end
    end

    nothing
end

#=
Plot waveforms in incremental way

Selected - true if waveform is captured by incremental capture
Plotted - true if waveform has been replotted in new color since start of incremental capture
=#
function plot_selected_waveforms(sc::Single_Channel,input::Array{T,2},mycolor) where T<:Real

    ctx=sc.ctx2
    s=sc.s
    o=sc.o

    set_line_width(ctx,2.0)
    set_source(ctx,sc.ctx2s)
    Cairo.translate(ctx,0.0,sc.h2/2)
    Gtk.scale(ctx,sc.w2/sc.wave_points,s)

    #=
    Reset waveforms that have changed since the start but are
    no longer selected
    =#
    for j=1:sc.buf.count
        if (!sc.selected[j])&(sc.plotted[j])
            move_to(ctx,1,(input[1,j]-o))
            for jj=2:size(input,1)
                line_to(ctx,jj,input[jj,j]-o)
            end
            sc.plotted[j]=false
        end
    end
    stroke(ctx)

    #=
    Plot selected waveforms in new color that have not
    yet been plotting in new color
    =#
    for i=1:sc.buf.count
        if (sc.selected[i])&(!sc.plotted[i])
            move_to(ctx,1,(input[1,i]-o))
            for jj=2:size(input,1)
                line_to(ctx,jj,input[jj,i]-o)
            end
            sc.plotted[i]=true
        end
    end
    set_line_width(ctx,0.5)
    select_color(ctx,mycolor)
    stroke(ctx)

    identity_matrix(ctx)
    nothing
end

#=
Show threshold on canvas
=#

function plot_thres(sc::Single_Channel)

    ctx = sc.ctx2

    line(ctx,1,sc.w2,sc.h2/2-sc.old_thres,sc.h2/2-sc.old_thres)
    set_line_width(ctx,5.0)
    set_source(ctx,sc.ctx2s)
    stroke(ctx)

    line(ctx,1,sc.w2,sc.h2/2-sc.thres,sc.h2/2-sc.thres)
    set_line_width(ctx,1.0)
    set_source_rgb(ctx,1.0,1.0,1.0)
    stroke(ctx)
    sc.old_thres=sc.thres
    nothing
end

function clear_c2(myc::Gtk.GtkCanvas,num)

    ctx = Gtk.getgc(myc)
    myheight=height(ctx)
    mywidth=width(ctx)

    set_source_rgb(ctx,0.0,0.0,0.0)
    paint(ctx)

    dashes = [10.0,10.0,10.0]
    set_dash(ctx, dashes, 0.0)

    for y = [myheight/6, myheight/3, myheight/2, myheight/6*4, myheight/6*5]
        line(ctx,1,mywidth,y,y)
    end

    for x = [.2*mywidth, .4*mywidth, .6*mywidth, .8*mywidth]
        line(ctx,x,x,1,myheight)
    end

    set_source_rgba(ctx,1.0,1.0,1.0,.5)
    stroke(ctx)

    set_dash(ctx,Float64[])

    line(ctx,1,mywidth,myheight,myheight)
    set_source_rgb(ctx,1.0,1.0,1.0)
    stroke(ctx)

    line(ctx,1,mywidth,myheight/2,myheight/2)
    stroke(ctx)

    move_to(ctx,10,10)
    show_text(ctx,string(num))

    nothing
end

function pause_state_cb(widgetptr::Ptr,user_data::Tuple{Single_Channel,Int64})

    sc, myid = user_data

    sc.pause_state = myid

    nothing
end

#=
Callback for how mouse interacts with canvas
=#
function canvas_press_win(widget::Ptr,param_tuple,user_data::Tuple{Single_Channel})

    sc, = user_data
    event = unsafe_load(param_tuple)

    sc.click_button=event.button

    if event.button == 1 #left click captures window
        sc.mi=(event.x,event.y)
        if sc.pause_state==1
            rubberband_start(sc,event.x,event.y)
        elseif sc.pause_state == 2
            if sc.pause
                draw_start(sc,event.x,event.y,sc.temp)
            end
        end
    elseif event.button == 3 #right click refreshes window
        if !sc.pause
            clear_c2(sc.c2,sc.spike)
            sc.ctx2=Gtk.getgc(sc.c2)
            sc.ctx2s=copy(sc.ctx2)
            sc.buf.ind=1
            sc.buf.count=1
            if sc.sort_cb
                draw_templates(sc)
            end
        else
            sc.mi=(event.x,event.y)
            if sc.pause_state==1
                rubberband_start(sc,event.x,event.y,3)
            elseif sc.pause_state == 2

            end
        end
    end
    nothing
end

function pause_cb(widgetptr::Ptr,user_data::Tuple{Single_Channel})

    sc, = user_data

    widget = convert(ToggleButton, widgetptr)

    if getproperty(widget,:active,Bool)
        sc.pause=true
        change_button_label(widget,"Resume")
        for i=1:length(sc.buf.mask)
            sc.buf.mask[i]=true
        end
    else
        sc.pause=false
        change_button_label(widget,"Pause")
        sc.hold=false
    end

    nothing
end

function clear_button_cb(widget::Ptr,user_data::Tuple{SpikeSorting.Single_Channel})

    sc, = user_data
    SpikeSorting.clear_c2(sc.c2,sc.spike)
    sc.ctx2=Gtk.getgc(sc.c2)
    sc.ctx2s=copy(sc.ctx2)
    #Sort Button
    if sc.sort_cb
        SpikeSorting.draw_templates(sc)
    end

    nothing
end

function restore_button_cb(widget::Ptr,user_data::Tuple{SpikeSorting.Single_Channel})

    sc, = user_data

    if sc.pause
        for i=1:sc.buf.ind
            sc.buf.mask[i]=true
        end

        sc.buf.replot=true
    end

    nothing
end

function canvas_release_template(widget::Ptr,param_tuple,user_data::Tuple{Buffer,Single_Channel})

    buf,sc = user_data
    event = unsafe_load(param_tuple)

    clus=buf.selected_clus

    (x1,x2,y1,y2)=SpikeSorting.coordinate_transform(sc,event)

    if event.button==1

        if sc.pause_state == 1
            #Waveforms selected with a left click rubberband in the paused display will be used to form a cluster template for the currently selected cluster.

            if (clus==0) #do nothing if zeroth cluster
            else

                buf.selected=trues(buf.ind)
                find_intersected_waveforms(buf.spikes,buf.selected,buf.ind,x1,y1,x2,y2)

                buf.c_changed=true
            end

            if (clus>0)&((buf.count>0)&(sc.pause))
                buf.replot=true
            end
        elseif sc.pause_state == 2
            if sc.pause
                buf.c_changed=true
            end
        end
    elseif event.button==3

        if sc.pause_state == 1
            #Waveforms selected with right click in the paused display will be removed from display. If a cluster is selected, these waveforms will also be removed from that selected cluster

            if sc.pause
                find_intersected_waveforms(buf.spikes,buf.mask,buf.count,x1,y1,x2,y2)

                if clus>0
                    buf.selected=.!(buf.clus.==clus)
                    buf.c_changed=true
                end
                buf.replot=true
            end
        elseif sc.pause_state == 2
            if sc.pause

            end
        end
    end

    nothing
end

function win_resize_cb(widget::Ptr,param_tuple,user_data::Tuple{Single_Channel})

    sc, = user_data

    ctx=Gtk.getgc(sc.c2)

    if (height(ctx)!=sc.h2)|(width(ctx)!=sc.w2)

        sc.ctx2=ctx
        sc.ctx2s=copy(sc.ctx2)
        sc.w2=width(sc.ctx2)
        sc.h2=height(sc.ctx2)

        setproperty!(sc.adj_thres,:upper,sc.h2/2)
        setproperty!(sc.adj_thres,:lower,-sc.h2/2)
    end

    nothing
end

#=
Redraw all spikes shown in paused view
=#
function replot_all_spikes(sc::Single_Channel)

    SpikeSorting.clear_c2(sc.c2,sc.spike)
    sc.ctx2=Gtk.getgc(sc.c2)
    sc.ctx2s=copy(sc.ctx2)

    ctx=sc.ctx2s
    s=sc.s
    o=sc.o

    Cairo.translate(ctx,0.0,sc.h2/2)
    Cairo.scale(ctx,sc.w2/sc.wave_points,s)

    for i=1:(sc.total_clus+1)
        for j=1:sc.buf.ind
            if (sc.buf.clus[j]==(i-1))&(sc.buf.mask[j])
                move_to(ctx,1,(sc.buf.spikes[1,j]-o))
                for jj=2:size(sc.buf.spikes,1)
                    line_to(ctx,jj,sc.buf.spikes[jj,j]-o)
                end
            end
        end
        set_line_width(ctx,0.5)
        select_color(ctx,i)
        stroke(ctx)
    end
    identity_matrix(ctx)
    set_source(sc.ctx2,ctx)
    mask_surface(sc.ctx2,ctx,0.0,0.0)
    fill(sc.ctx2)
    reveal(sc.c2)
    nothing
end

function mask_surface(ctx,s,x,y)
    ccall((:cairo_mask_surface,Cairo._jl_libcairo),Void,(Ptr{Void},Ptr{Void},Float64,Float64),ctx.ptr,s.surface.ptr,x,y)
end


text:
#%% cell
#=
Now that we have a model, we will construct one more helper function that takes
the control variables $(S, F, C)$ and shocks $(\delta, R, \eta_G, \eta_a,
\eta_L, \eta_R)$ and applies equilibrium conditions to produce $(\pi, \delta',
Y, L, Y_n, R')$. We will use this function in both the solution and simulation
routines below.
=#

#%% cell
function Base.step(m::Model, S, F, C, δ0, R0, ηG, ηa, ηL, ηR)
    @unpack Θ, ϵ, gbar, ϑ, γ, β, μ, ϕ_π, ϕ_y, πstar = m.p
    πs = m.s.π

    # Compute pie(t) from condition (35) in MM (2015)
    π0 = @. ((1-(1-Θ)*(S/F)^(1-ϵ))/Θ)^(1/(ϵ-1))

    # Compute delta(t) from condition (36) in MM (2015)
    δ1 = @. ((1-Θ)*((1-Θ*π0^(ϵ-1))/(1-Θ))^(ϵ/(ϵ-1))+Θ*π0^ϵ/δ0)^(-1.0)

    # Compute Y(t) from condition (38) in MM (2015)
    Y0 = @. C/(1-gbar/exp(ηG))

    # Compute L(t) from condition (37) in MM (2015)
    L0 = @. Y0/exp(ηa)/δ1

    #  Compute Yn(t) from condition (31) in MM (2015)
    Yn0 = @. (exp(ηa)^(1+ϑ)*(1-gbar/exp(ηG))^(-γ)/exp(ηL))^(1/(ϑ+γ))

    # Compute R(t) from conditions (27), (39) in MM (2015) -- Taylor rule
    R1 = @. πstar / β * (R0*β/πstar)^μ*((π0/πstar)^ϕ_π * (Y0/Yn0)^ϕ_y)^(1-μ)*exp(ηR)

    π0, δ1, Y0, L0, Yn0, R1
end
#%% cell
#=
### Solution routine
=#

#%% cell
# construct an initial guess for the solution
function initial_coefs(m::Model, degree)
    npol = size(m.g.X0_G[degree], 2)
    coefs = fill(1e-5, npol, 3)
    coefs[1, :] = [m.s.S, m.s.F, m.s.C^(-m.p.γ)]
    coefs
end

function solve(m::Model; verbose::Bool=false)
    # simplify notation
    n, n_nodes = size(m.g.ηR, 1), length(m.g.ω_nodes)

    ## allocate memory
    # euler equations
    e = zeros(n, 3)

    # previous iteration S, F, C
    S0_old_G = ones(n)
    F0_old_G = ones(n)
    C0_old_G = ones(n)

    # current iteration S, F, C
    S0_new_G = ones(n)
    F0_new_G = ones(n)
    C0_new_G = ones(n)

    # future S, F, C
    S1 = zeros(n, n_nodes)
    F1 = zeros(n, n_nodes)
    C1 = zeros(n, n_nodes)
    π1 = zeros(n, n_nodes)

    local coefs, start_time

    degs = m.p.deg > 1 ? [1, m.p.deg] : [m.p.deg]

    for deg in degs
        # set up matrices for this degree
        X0_G = m.g.X0_G[deg]

        # future basis matrix for S, F, C
        X1 = Array{Float64}(undef, n, n_complete(8, deg))

        # initialize coefs
        if deg <= 1
            coefs = initial_coefs(m, deg)
        else
            # for higher order, start with degree 2 coefficients
            coefs = X0_G\e
            # old_coefs = copy(coefs)
            # coefs = initial_coefs(m, deg)
            # coefs[1:size(old_coefs, 1), :] = old_coefs
        end

        err = 1.0
        it = 0

        start_time = time()
        # solve at this degree of complete polynomial
        while err > m.p.tol
            it += 1;
            # Current choices (at t)
            # ------------------------------
            S0 = X0_G*coefs[:, 1]                # Compute S(t) using coefs
            F0 = X0_G*coefs[:, 2]                # Compute F(t) using coefs
            C0 = X0_G*coefs[:, 3]
            C0 .^= (-1/m.p.γ)                    # Compute C(t) using coefs

            π0, δ1, Y0, L0, Yn0, R1 = step(m, S0, F0, C0, m.g.δ, m.g.R, m.g.ηG,
                                           m.g.ηa, m.g.ηL, m.g.ηR)

            if m.p.zlb R1 .= max.(R1, 1.0) end

            for u in 1:n_nodes
                # Form complete polynomial of degree "Degree" (at t+1) on future state
                complete_polynomial!(
                    X1,
                    hcat(log.(R1), log.(δ1), m.g.ηR1[:, u], m.g.ηa1[:, u],
                         m.g.ηL1[:, u], m.g.ηu1[:, u], m.g.ηB1[:, u], m.g.ηG1[:, u]),
                    deg
                )

                S1[:, u] = X1*coefs[:, 1]                # Compute S(t+1)
                F1[:, u] = X1*coefs[:, 2]                # Compute F(t+1)
                C1[:, u] = (X1*coefs[:, 3]).^(-1/m.p.γ)  # Compute C(t+1)
            end

            # Compute next-period π using condition
            # (35) in MM (2015)
            @. π1 = ((1-(1-m.p.Θ).*(S1./F1).^(1-m.p.ϵ))/m.p.Θ).^(1/(m.p.ϵ-1))

            # Evaluate conditional expectations in the Euler equations
            #---------------------------------------------------------
            e[:, 1] = @.(exp(m.g.ηu)*exp(m.g.ηL)*L0^m.p.ϑ*Y0/exp(m.g.ηa) + (m.p.β*m.p.Θ*π1^m.p.ϵ*S1))*m.g.ω_nodes
            e[:, 2] = @.(exp(m.g.ηu)*C0^(-m.p.γ)*Y0 + (m.p.β*m.p.Θ*π1^(m.p.ϵ-1)*F1))*m.g.ω_nodes
            e[:, 3] = @.(m.p.β*exp(m.g.ηB)/exp(m.g.ηu)*R1*((exp(m.g.ηu1).*C1^(-m.p.γ)/π1)))*m.g.ω_nodes

            # Variables of the current iteration
            #-----------------------------------
            copy!(S0_new_G, S0)
            copy!(F0_new_G, F0)
            copy!(C0_new_G, C0)

            # Compute and update the coefficients of the decision functions
            # -------------------------------------------------------------
            coefs_hat = X0_G\e   # Compute the new coefficients of the decision
                                 # functions using a backslash operator

            # Update the coefficients using damping
            coefs = m.p.damp*coefs_hat + (1-m.p.damp)*coefs

            # Evaluate the percentage (unit-free) difference between the values
            # on the grid from the previous and current iterations
            # -----------------------------------------------------------------
            # The convergence criterion is adjusted to the damping parameters
            err = mean(abs, 1.0.-S0_new_G./S0_old_G) +
                  mean(abs, 1.0.-F0_new_G./F0_old_G) +
                  mean(abs, 1.0.-C0_new_G./C0_old_G)

            if (it % 20 == 0) && verbose
                @printf "On iteration %d err is %6.7e\n" it err
            end

            # Store the obtained values for S(t), F(t), C(t) on the grid to
            # be used on the subsequent iteration in Section 10.2.6
            #-----------------------------------------------------------------------
            copy!(S0_old_G, S0_new_G)
            copy!(F0_old_G, F0_new_G)
            copy!(C0_old_G, C0_new_G)

        end
    end

    coefs, time() - start_time
end

#%% cell
#=
### Simulation
=#

#%% cell
struct Simulation
    # shocks
    ηR::Vector{Float64}
    ηa::Vector{Float64}
    ηL::Vector{Float64}
    ηu::Vector{Float64}
    ηB::Vector{Float64}
    ηG::Vector{Float64}

    # variables
    δ::Vector{Float64}
    R::Vector{Float64}
    S::Vector{Float64}
    F::Vector{Float64}
    C::Vector{Float64}
    π::Vector{Float64}
    Y::Vector{Float64}
    L::Vector{Float64}
    Yn::Vector{Float64}
    w::Vector{Float64}
end

function Simulation(m::Model, coefs::Matrix)

    # 11. Simualating a time-series solution
    #---------------------------------------
    _path = joinpath(@__DIR__, "..", "epsi_test_NK.mat")
    rands = (matread(_path)["epsi_test_NK"])::Matrix{Float64}
    capT = size(rands, 1)

    # Initialize the values of 6 exogenous shocks
    #--------------------------------------------
    ηR = zeros(capT)
    ηa = zeros(capT)
    ηL = zeros(capT)
    ηu = zeros(capT)
    ηB = zeros(capT)
    ηG = zeros(capT)

    # Generate the series for shocks
    #-------------------------------
    @inbounds for t in 1:capT-1
        ηR[t+1] = m.p.ρηR*ηR[t] + m.p.σηR*rands[t, 1]
        ηa[t+1] = m.p.ρηa*ηa[t] + m.p.σηa*rands[t, 2]
        ηL[t+1] = m.p.ρηL*ηL[t] + m.p.σηL*rands[t, 3]
        ηu[t+1] = m.p.ρηu*ηu[t] + m.p.σηu*rands[t, 4]
        ηB[t+1] = m.p.ρηB*ηB[t] + m.p.σηB*rands[t, 5]
        ηG[t+1] = m.p.ρηG*ηG[t] + m.p.σηG*rands[t, 6]
    end

    δ  = ones(capT+1) # Time series of delta(t)
    R  = ones(capT+1) # Time series of R(t)
    S  = Array{Float64}(undef, capT)   # Time series of S(t)
    F  = Array{Float64}(undef, capT)   # Time series of F(t)
    C  = Array{Float64}(undef, capT)   # Time series of C(t)
    π  = Array{Float64}(undef, capT)   # Time series of π(t)
    Y  = Array{Float64}(undef, capT)   # Time series of Y(t)
    L  = Array{Float64}(undef, capT)   # Time series of L(t)
    Yn = Array{Float64}(undef, capT)   # Time series of Yn(t)
    w  = Array{Float64}(undef, capT)

    pol_bases = Array{Float64}(undef, 1, size(coefs, 1))
    @inbounds for t in 1:capT
        # Construct the matrix of explanatory variables "pol_bases" on the series
        # of state variables; columns of "pol_bases" are given by the basis
        # functions of the polynomial of degree 2
        complete_polynomial!(
            pol_bases,
            hcat(log(R[t]), log(δ[t]), ηR[t], ηa[t], ηL[t], ηu[t], ηB[t], ηG[t]),
            m.p.deg
        )
        S[t], F[t], MU = pol_bases*coefs
        C[t] = (MU).^(-1/m.p.γ)

        π[t], δ[t+1], Y[t], L[t], Yn[t], R[t+1] = step(m, S[t], F[t], C[t], δ[t],
                                                       R[t], ηG[t], ηa[t], ηL[t], ηR[t])

        # Compute real wage
        w[t] = exp(ηL[t])*(L[t]^m.p.ϑ)*(C[t]^m.p.γ)

        # If ZLB is imposed, set R(t)=1 if ZLB binds
        if m.p.zlb; R[t+1] = max(R[t+1],1.0); end
    end

    Simulation(ηR, ηa, ηL, ηu, ηB, ηG, δ, R, S, F, C, π, Y, L, Yn, w)
end

#%% cell
#=
### Accuracy
=#

#%% cell
struct Residuals
    resids::Matrix{Float64}
end

function Residuals(m::Model, coefs::Matrix, sim::Simulation; burn::Int=200)
    capT = length(sim.w)
    resids = zeros(9, capT)

    # Integration method for evaluating accuracy
    # ------------------------------------------
    # Monomial integration rule with 2N^2+1 nodes
    ϵ_nodes, ω_nodes = qnwmonomial2(vcov(m.p))
    n_nodes = length(ω_nodes)

    # Allocate for arrays needed in the loop
    basis_mat = Array{Float64}(undef, n_nodes, 8)
    X1 = Array{Float64}(undef, n_nodes, size(coefs, 1))

    ηR1 = Array{Float64}(undef, n_nodes)
    ηa1 = Array{Float64}(undef, n_nodes)
    ηL1 = Array{Float64}(undef, n_nodes)
    ηu1 = Array{Float64}(undef, n_nodes)
    ηB1 = Array{Float64}(undef, n_nodes)
    ηG1 = Array{Float64}(undef, n_nodes)

    for t in 1:capT                 # For each given point,
        # Take the corresponding value for shocks at t
        #---------------------------------------------
        ηR0 = sim.ηR[t]  # ηR(t)
        ηa0 = sim.ηa[t]  # ηa(t)
        ηL0 = sim.ηL[t]  # ηL(t)
        ηu0 = sim.ηu[t]  # ηu(t)
        ηB0 = sim.ηB[t]  # ηB(t)
        ηG0 = sim.ηG[t]  # ηG(t)

        # Exctract time t values for all other variables (and t+1 for R, δ)
        #------------------------------------------------------------------
        R0  = sim.R[t]   # R(t-1)
        δ0  = sim.δ[t]   # δ(t-1)
        R1  = sim.R[t+1] # R(t)
        δ1  = sim.δ[t+1] # δ(t)

        L0  = sim.L[t]   # L(t)
        Y0  = sim.Y[t]   # Y(t)
        Yn0 = sim.Yn[t]  # Yn(t)
        π0  = sim.π[t]   # π(t)
        S0  = sim.S[t]   # S(t)
        F0  = sim.F[t]   # F(t)
        C0  = sim.C[t]   # C(t)

        # Fill basis matrix with R1, δ1 and shocks
        #-----------------------------------------
        # Note that we do not premultiply by standard deviations as ϵ_nodes
        # already include them. All these variables are vectors of length n_nodes
        copy!(ηR1, ηR0*m.p.ρηR .+ ϵ_nodes[:, 1])
        copy!(ηa1, ηa0*m.p.ρηa .+ ϵ_nodes[:, 2])
        copy!(ηL1, ηL0*m.p.ρηL .+ ϵ_nodes[:, 3])
        copy!(ηu1, ηu0*m.p.ρηu .+ ϵ_nodes[:, 4])
        copy!(ηB1, ηB0*m.p.ρηB .+ ϵ_nodes[:, 5])
        copy!(ηG1, ηG0*m.p.ρηG .+ ϵ_nodes[:, 6])

        basis_mat[:, 1] .= log.(R1)
        basis_mat[:, 2] .= log.(δ1)
        basis_mat[:, 3] = ηR1
        basis_mat[:, 4] = ηa1
        basis_mat[:, 5] = ηL1
        basis_mat[:, 6] = ηu1
        basis_mat[:, 7] = ηB1
        basis_mat[:, 8] = ηG1

        # Future choices at t+1
        #----------------------
        # Form a complete polynomial of degree "Degree" (at t+1) on future state
        # variables; n_nodes-by-npol
        complete_polynomial!(X1, basis_mat, m.p.deg)

        # Compute S(t+1), F(t+1) and C(t+1) in all nodes using coefs
        S1 = X1*coefs[:, 1]
        F1 = X1*coefs[:, 2]
        C1 = (X1*coefs[:, 3]).^(-1/m.p.γ)

        # Compute π(t+1) using condition (35) in MM (2015)
        π1 = @. ((1-(1-m.p.Θ)*(S1./F1).^(1-m.p.ϵ))/m.p.Θ).^(1/(m.p.ϵ-1))

        # Compute residuals for each of the 9 equilibrium conditions
        #-----------------------------------------------------------
        resids[1, t] = 1-dot(ω_nodes,
            @. (exp(ηu0)*exp(ηL0)*L0^m.p.ϑ*Y0/exp(ηa0) + m.p.β*m.p.Θ*π1.^m.p.ϵ.*S1)/S0
        )
        resids[2, t] = 1 - dot(ω_nodes,
            @. (exp(ηu0)*C0^(-m.p.γ)*Y0 + m.p.β*m.p.Θ*π1.^(m.p.ϵ-1).*F1)/F0
        )
        resids[3, t] = 1.0 -dot(ω_nodes,
            @. (m.p.β*exp(ηB0)/exp(ηu0)*R1*exp.(ηu1).*C1.^(-m.p.γ)./π1)/C0^(-m.p.γ)
        )
        resids[4, t] = 1-((1-m.p.Θ*π0^(m.p.ϵ-1))/(1-m.p.Θ))^(1/(1-m.p.ϵ))*F0/S0
        resids[5, t] = 1-((1-m.p.Θ)*((1-m.p.Θ*π0^(m.p.ϵ-1))/(1-m.p.Θ))^(m.p.ϵ/(m.p.ϵ-1)) + m.p.Θ*π0^m.p.ϵ/δ0)^(-1)/δ1
        resids[6, t] = 1-exp(ηa0)*L0*δ1/Y0
        resids[7, t] = 1-(1-m.p.gbar/exp(ηG0))*Y0/C0
        resids[8, t] = 1-(exp(ηa0)^(1+m.p.ϑ)*(1-m.p.gbar/exp(ηG0))^(-m.p.γ)/exp(ηL0))^(1/(m.p.ϑ+m.p.γ))/Yn0
        resids[9, t] = 1-m.s.π/m.p.β*(R0*m.p.β/m.s.π)^m.p.μ*((π0/m.s.π)^m.p.ϕ_π * (Y0/Yn0)^m.p.ϕ_y)^(1-m.p.μ)*exp(ηR0)/R1   # Taylor rule

        # If the ZLB is imposed and R>1, the residuals in the Taylor rule (the
        # 9th equation) are zero
        if m.p.zlb && R1 <= 1; resids[9, t] = 0.0; end

    end
    # discard the first burn observations
    Residuals(resids[:, burn+1:end])
end

Statistics.mean(r::Residuals) = log10(mean(abs, r.resids))
Base.max(r::Residuals) = log10(maximum(abs, r.resids))
max_E(r::Residuals) = log10.(maximum(abs, r.resids, dims=2))[:]

max_stars_repo_name: sglyon/CLMMJuliaPythonMatlab

max_issues_repo_name: sglyon/CLMMJuliaPythonMatlab

max_forks_repo_name: sglyon/CLMMJuliaPythonMatlab


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     formats: jl:light,ipynb
#     text_representation:
#       extension: .jl
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.5.1
#   kernelspec:
#     display_name: Julia 1.4.2
#     language: julia
#     name: julia-1.4
# ---

# # Density Functional Theory: Grid
# ## I. Theoretical Overview
# This tutorial will discuss the basics of DFT and discuss the grid used to evaluate DFT quantities.
# As with HF, DFT aims to solve the generalized eigenvalue problem:
#
# $$\sum_{\nu} F_{\mu\nu}C_{\nu i} = \epsilon_i\sum_{\nu}S_{\mu\nu}C_{\nu i}$$
# $${\bf FC} = {\bf SC\epsilon},$$
#
# While with HF the Fock matrix is constructed as:
#
# $$F^{HF}_{\mu\nu} = H_{\mu\nu} + 2J[D]_{\mu\nu} - K[D]_{\mu\nu}$$
#
# $$D_{\mu\nu} = C_{\mu i} C_{\nu i},$$
#
# with DFT we generalize this construction slightly to:
# $$F^{DFT}_{\mu\nu} = H_{\mu\nu} + 2J[D]_{\mu\nu} - \zeta K[D]_{\mu\nu} + V^{\rm{xc}}_{\mu\nu}.$$
#
# $\zeta$ is an adjustable parameter where we can vary the amount of exact (HF) exchange and $V$ is the DFT potenital which typically attempts to add dynamical correlation in the self-consistent field methodolgy.
#
#

# ## 2. Examining the Grid
# We will discuss the evaluation and manipulation of the grid.

using PyCall: pyimport
psi4 = pyimport("psi4")
np   = pyimport("numpy") # used only to cast to Psi4 arrays
import PyPlot: plt, matplotlib
matplotlib.use(backend="MacOSX")

build_superfunctional = nothing
if VersionNumber(psi4.__version__) >= v"1.3a1"
    build_superfunctional = psi4.driver.dft.build_superfunctional
else
    build_superfunctional = psi4.driver.dft_funcs.build_superfunctional
end

# Set computation options and molecule, any single atom will do.

mol = psi4.geometry("He")
psi4.set_options(Dict("BASIS" => "cc-pvdz",
                      "DFT_SPHERICAL_POINTS" => 50,
                      "DFT_RADIAL_POINTS" => 12))

basis = psi4.core.BasisSet.build(mol, "ORBITAL", "CC-PVDZ")
sup = build_superfunctional("PBE", true)[1]
Vpot = psi4.core.VBase.build(basis, sup, "RV")
Vpot.initialize()

x, y, z, w = Vpot.get_np_xyzw()
R = @. sqrt(x^2 + y^2 + z^2);

fig, ax = plt.subplots()
ax.scatter(x, y, c=w)
#ax.set_xscale("log")

fig = plt.figure()
ax = fig.add_subplot(111, projection="3d")
mask = R .> 8
p = ax.scatter(x[mask], y[mask], z[mask], c=w[mask], marker="o")
plt.colorbar(p)
ax.set_xlabel("X")
ax.set_ylabel("Y")
ax.set_zlabel("Z")

# +
mol = psi4.geometry("""
 O
 H 1 1.1
 H 1 1.1 2 104
""")
mol.update_geometry()
psi4.set_options(Dict("BASIS" => "cc-pvdz",
                      "DFT_SPHERICAL_POINTS" => 26,
                      "DFT_RADIAL_POINTS" => 12))

basis = psi4.core.BasisSet.build(mol, "ORBITAL", "CC-PVDZ")
sup = build_superfunctional("PBE", true)[1]
Vpot = psi4.core.VBase.build(basis, sup, "RV")
Vpot.initialize()
x, y, z, w = Vpot.get_np_xyzw()
R = @. sqrt(x^2 + y^2 + z^2)

fig = plt.figure()
ax = fig.add_subplot(111, projection="3d")
mask = R .> 0
p = ax.scatter(x[mask], y[mask], z[mask], c=w[mask], marker="o")
plt.colorbar(p)
ax.set_xlabel("X")
ax.set_ylabel("Y")
ax.set_zlabel("Z")
# -

# ## Refs:
# - Koch, W. and Holthausen, M.C., **2001**, A Chemist’s Guide to Density Functional Theory, 2nd, Wiley-VCH, Weinheim.
# - Kohn, W. and Sham, L. *J, Phys. Rev.*, **1965**, *140*, A1133- A1138
# - Becke, A.D., *J. Chem. Phys.*, **1988**, *88*, 2547
# - Treutler, O. and Ahlrichs, R., *J. Chem. Phys.*, **1995**, *102*, 346
# - Gill, P.M.W., Johnson, B.G., and Pople, J.A., *Chem. Phys. Lett.*, **1993,209 (5), pp. 506, 16 July 1993.


text:
using BHAPtfem, FinEtools

ProjDir = @__DIR__
cd(ProjDir)

h = 1.0
l = 3.0
nh = 1
nl = 3
nc = 6

fens, fes  = H8block(h, l, 2.0 * pi, nh, nl, nc) # generate a block
# Shape into a cylinder
R = zeros(3, 3)
for i = 1:count(fens)
    x, y, z = fens.xyz[i,:];
    rotmat3!(R, [0, z, 0])
    Q = [cos(psi * pi / 180) sin(psi * pi / 180) 0;
        -sin(psi * pi / 180) cos(psi * pi / 180) 0;
        0 0 1]
    fens.xyz[i,:] = reshape([x + Rmed - h / 2, y - l / 2, 0], 1, 3) * Q * R;
end
# Merge the nodes where the tube  closes up
candidates = selectnode(fens, box = boundingbox([Rmed - h -Inf 0.0; Rmed + h +Inf 0.0]), inflate = tolerance)
fens, fes = mergenodes(fens, fes,  tolerance, candidates);

rin = 0.6*phun("m")
rex = 0.7*phun("m")
nr = 1
nc = 6
Angl = 0.0

fens, fes = Q4annulus(rin::FFlt, rex::FFlt, nr::FInt, nc::FInt, Angl::FFlt)


text:
using Weave

weave(Pkg.dir("Weave","examples","gadfly_md_sample.jmd"),
     informat="markdown", out_path = "build/examples", doctype = "pandoc")

 weave(Pkg.dir("Weave","examples","gadfly_md_sample.jmd"),
      informat="markdown", out_path = "build/examples", doctype = "md2html")

weave(Pkg.dir("Weave","examples","gadfly_md_sample.jmd"),
        informat="markdown", out_path = "build/examples", doctype = "md2pdf")

cp(Pkg.dir("Weave","examples","gadfly_md_sample.jmd"),
    "build/examples/gadfly_md_sample.jmd", remove_destination = true)

cp("build/examples/gadfly_md_sample.md",
      "build/examples/gadfly_md_sample.txt", remove_destination = true)

weave(Pkg.dir("Weave","examples","FIR_design.jl"), out_path = "build/examples")
weave(Pkg.dir("Weave","examples","FIR_design.jl"), doctype = "md2pdf", out_path = "build/examples")
cp(Pkg.dir("Weave","examples","FIR_design.jl"),
    "build/examples/FIR_design.jl", remove_destination = true)


text:
module Bresenham

function line(f, x0::Int, y0::Int, x1::Int, y1::Int)
  dx = 2 * (x1 - x0)
  dy = 2 * (y1 - y0)

  if dx != 0
    if dy != 0
      if abs(dx) >= abs(dy)
        e = dx / 2
        y = y0
        for x = x0:sign(dx):x1
          f(x, y)
          # (a) following code runs even if
          # the loop exits after this iteration
          # but we don't care about that side effect
          e -= sign(dx) * sign(dy) * dy
          if sign(dx) > 0 && e < 0 ||
             sign(dx) < 0 && e >= 0
            y += sign(dy)
            e += dx
          end
        end
      else # abs(dx) < abs(dy) (steep case)
        e = dy / 2
        x = x0
        for y = y0:sign(dy):y1
          f(x, y)
          # see (a)
          e -= sign(dx) * sign(dy) * dx
          if sign(dy) > 0 && e < 0 ||
             sign(dy) < 0 && e >= 0
            x += sign(dx)
            e += dy
          end
        end
      end
    else # dy == 0 (horizontal line)
      for x = x0:sign(dx):x1
        f(x, y0)
      end
    end
  else # dx == 0 (vertical line)
    for y = y0:sign(dy):y1
      f(x0, y)
    end
  end
end # function

function circle(f, x0::Int, y0::Int, r::Int)
  x = r - 1
  y = 0
  dx = 1
  dy = 1
  err = dx - (r << 1)

  while x >= y
    f(x0 + x, y0 + y)
    f(x0 + y, y0 + x)
    f(x0 - y, y0 + x)
    f(x0 - x, y0 + y)
    f(x0 - x, y0 - y)
    f(x0 - y, y0 - x)
    f(x0 + y, y0 - x)
    f(x0 + x, y0 - y)

    if err <= 0
      y += 1
      err += dy
      dy += 2
    end

    if err > 0
      x -= 1
      dx += 2
      err += dx - (r << 1)
    end
  end
end # function

function ordered_circle(f, x0::Int, y0::Int, r::Int)
  pts = Vector{Tuple{Int, Int}}()
  Bresenham.circle(x0, y0, r) do x, y
    push!(pts, (x, y))
  end

  unique!(sort!(pts, by=(X -> atan(X[2], X[1]))))
  for X in pts
    f(X...)
  end
end # function

export line, circle, ordered_circle

end # module


text:
using Test
using IJulia

@testset "stdio" begin

    mktemp() do path, io
        redirect_stdout(IJulia.IJuliaStdio(io, "stdout")) do
            println(Base.stdout, "stdout")
            println("print")
        end
        flush(io)
        seek(io, 0)
        @test read(io, String) == "stdout\nprint\n"
        if VERSION < v"1.7.0-DEV.254"
            @test_throws ArgumentError redirect_stdout(IJulia.IJuliaStdio(io, "stderr"))
            @test_throws ArgumentError redirect_stdout(IJulia.IJuliaStdio(io, "stdin"))
            @test_throws ArgumentError redirect_stderr(IJulia.IJuliaStdio(io, "stdout"))
            @test_throws ArgumentError redirect_stderr(IJulia.IJuliaStdio(io, "stdin"))
            @test_throws ArgumentError redirect_stdin(IJulia.IJuliaStdio(io, "stdout"))
            @test_throws ArgumentError redirect_stdin(IJulia.IJuliaStdio(io, "stderr"))
        end
    end

    mktemp() do path, io
        redirect_stderr(IJulia.IJuliaStdio(io, "stderr")) do
            println(Base.stderr, "stderr")
        end
        flush(io)
        seek(io, 0)
        @test read(io, String) == "stderr\n"
    end

    mktemp() do path, io
        redirect_stdin(IJulia.IJuliaStdio(io, "stdin")) do
            # We can't actually do anything here because `IJuliaexecute_msg` has not
            # yet been initialized, so we just make sure that redirect_stdin does
            # not error.
        end
    end

end


text:
include("ps3000CommonElements.jl")

include("ps3000_close_unit.jl")
include("ps3000_flash_led.jl")
include("ps3000_get_streaming_last_values.jl")
include("ps3000_get_streaming_values.jl")
include("ps3000_get_streaming_values_no_aggregation.jl")
include("ps3000_get_timebase.jl")
include("ps3000_get_times_and_values.jl")
include("ps3000_get_unit_info.jl")
include("ps3000_get_values.jl")
include("ps3000_open_unit.jl")
include("ps3000_open_unit_async.jl")
include("ps3000_open_unit_progress.jl")
include("ps3000_overview_buffer_status.jl")
include("ps3000_ready.jl")
include("ps3000_run_block.jl")
include("ps3000_run_streaming.jl")
include("ps3000_run_streaming_ns.jl")
include("ps3000_save_streaming_data.jl")
include("ps3000_set_channel.jl")
include("ps3000_set_ets.jl")
# include("ps3000_set_siggen.jl")  # I don't have a signal generator in my scope
include("ps3000_set_trigger.jl")  # and ps3000_set_trigger2
include("ps3000_stop.jl")
include("ps3000setadvtriggerchannelconditions.jl")
include("ps3000setadvtriggerchanneldirections.jl")
include("ps3000setadvtriggerchannelproperties.jl")
include("ps3000setadvtriggerdelay.jl")
include("ps3000setpulsewidthqualifier.jl")


text:
# TODO add bangs for functions that modify argument(s)
module SententialDecisionDiagrams

using Parameters

include("sddapi.jl")
include("fnf.jl")

using .SddLibrary

struct VTree
    vtree::Ptr{SddLibrary.VTree_c}
end

struct SddManager
    manager::Ptr{SddLibrary.SddManager_c}
end

struct SddNode
    node::Ptr{SddLibrary.SddNode_c}
    manager::Ptr{SddLibrary.SddManager_c}
end

struct PrimeSub
    prime::SddNode
    sub::SddNode
end

struct WmcManager
    manager::Ptr{SddLibrary.WmcManager_c}
end

function str_to_char(vtree_type::String)::Ptr{UInt8}
    return pointer(vtree_type)
end

# SDD MANAGER FUNCTIONS
function sdd_manager(vtree::VTree)::SddManager
    manager = SddLibrary.sdd_manager_new(vtree.vtree)
    return SddManager(manager)
end
function sdd_manager(var_count::Integer,  auto_gc_and_minimize::Bool)::SddManager
    manager = SddLibrary.sdd_manager_create(convert(SddLibrary.SddLiteral, var_count), auto_gc_and_minimize)
    return SddManager(manager)
end
# TODO function manager_new
function free(manager::SddManager)
    SddLibrary.sdd_manager_free(manager.manager)
end
function Base.print(manager::SddManager)
    SddLibrary.sdd_manager_print(manager.manager)
end
function auto_gc_and_minimize_on(manager::SddManager)
    SddLibrary.sdd_manager_auto_gc_and_minimize_on(manager.manager)
end
function sdd_manager_auto_gc_and_minimize_off(manager::SddManager)
    SddLibrary.sdd_manager_auto_gc_and_minimize_off(manager.manager)
end
function is_auto_gc_and_minimize_on(manager::SddManager)::Bool
    return convert(Bool, SddLibrary.sdd_manager_is_auto_gc_and_minimize_on(manager.manager))
end
# TODO void sdd_manager_set_minimize_function
function unset_minimize_function(manager::SddManager)
    SddLI.sdd_manager_unset_minimize_function(manager.manager)
end
function options(manager::SddManager)
    SddLI.sdd_manager_options(manager.manager)
end
# TODO void sdd_manager_set_options(void* options, SddManager* manager);
function is_var_used(var::Integer, manager::SddManager)::Bool
    return convert(Bool, SddLibrary.sdd_manager_is_var_used(convert(SddLibrary.SddLiteral, var), manager.manager))
end
function vtree_of_var(var::Integer, manager::SddManager)::VTree
    vtree = SddLibrary.sdd_manager_vtree_of_var(convert(SddLibrary.SddLiteral, var), manager.manager)
    return VTree(vtree)
end
# TODO Vtree* sdd_manager_lca_of_literals(int count, SddLiteral* literals, SddManager* manager);
function var_count(manager::SddManager)::SddLibrary.SddLiteral
    return SddLibrary.sdd_manager_var_count(manager.manager)
end
# TODO void sdd_manager_var_order(SddLiteral* var_order, SddManager *manager);
function add_var_before_first(manager::SddManager)
    SddLibrary.sdd_manager_add_var_before_first(manager.manager)
end
function add_var_after_last(manager::SddManager)
    SddLibrary.sdd_manager_add_var_after_last(manager.manager)
end
function add_var_before(manager::SddManager)
    SddLibrary.sdd_manager_add_var_before(manager.manager)
end
function add_var_after(manager::SddManager)
    SddLibrary.sdd_manager_add_var_after(manager.manager)
end

# TERMINAL SDDS
function literal(value::Bool, manager::SddManager)::SddNode
    if value
        node = SddLibrary.sdd_manager_true(manager.manager)
    else
        node = SddLibrary.sdd_manager_false(manager.manager)
    end
    return SddNode(node, manager.manager)
end
function literal(literal::Integer, manager::SddManager)::SddNode
    node = SddLibrary.sdd_manager_literal(convert(SddLibrary.SddLiteral, literal), manager.manager)
    return SddNode(node, manager.manager)
end

# SDD QUERIES AND TRANSFORMATIONS
function apply(node1::SddNode, node2::SddNode, op::Bool, manager::SddManager)::SddNode
    node = SddLibrary.sdd_apply(node1.node, node2.node, convert(SddLibrary.SddBoolOp, op), manager.manager)
    return SddNode(node, node1.manager)
end
function conjoin(node1::SddNode, node2::SddNode, manager::SddManager)::SddNode
    node = SddLibrary.sdd_conjoin(node1.node, node2.node, manager.manager)
    return SddNode(node, manager.manager)
end
function disjoin(node1::SddNode, node2::SddNode, manager::SddManager)::SddNode
    node = SddLibrary.sdd_disjoin(node1.node, node2.node, manager.manager)
    return SddNode(node, manager.manager)
end
function negate(node::SddNode, manager::SddManager)::SddNode
    node = SddLibrary.sdd_negate(node.node, node.manager)
    return SddNode(node, manager.manager)
end
function condition(lit::Integer, node::SddNode, manager::SddManager)::SddNode
    node = SddLibrary.sdd_condition(convert(SddLibrary.SddLiteral, lit), node.node, manager.manager)
    return SddNode(node, manager.manager)
end
function exists(lit::Integer, node::SddNode, manager::SddManager)::SddNode
    node = SddLibrary.sdd_exists(convert(SddLibrary.SddLiteral, lit), node.node, manager.manager)
    return SddNode(node, manager.manager)
end
function exists_multiple(exists_map::Array{Integer,1}, node::SddNode, manager::SddManager; static::Bool=false)::SddNode
    if !static
        node = SddLibrary.sdd_exists(convert(Array{Cint,1}, exists_map), node.node, manager.manager)
    else
        node = SddLibrary.sdd_exists_multiple_static(convert(Array{Cint,1}, exists_map), node.node, manager.manager)
    end
    return SddNode(node, manager.manager)
end
function for_all(lit::Integer, node::SddNode, manager::SddManager)::SddNode
    node = SddLibrary.sdd_forall(convert(SddLibrary.SddLiteral, lit), node.node, manager.manager)
    return SddNode(node, manager.manager)
end
function minimize_cardinality(node::SddNode, manager::SddManager; globally::Bool=false)::SddNode
    if !globally
        node = SddLibrary.sdd_minimize_cardinality(node.node, manager.manager)
    else
        node = SddLibrary.sdd_global_minimize_cardinality(node.node, manager.manager)
    end
    return SddNode(node, manager.manager)
end
function minimum_cardinality(node::SddNode)::SddLibrary.SddLiteral
    return SddLibrary.sdd_minimum_cardinality(node.node)
end
function model_count(node::SddNode, manager::SddManager; globally::Bool=false)::SddLibrary.SddModelCount
    if !globally
        return SddLibrary.sdd_model_count(node.node, manager.manager)
    else
        return SddLibrary.sdd_global_model_count(node.node, manager.manager)
    end
end

# SDD NAVIGATION
function is_true(node::SddNode)::Bool
    res = SddLibrary.sdd_node_is_true(node.node)
    return convert(Bool, res)
end
function is_false(node::SddNode)::Bool
    res = SddLibrary.sdd_node_is_false(node.node)
    return convert(Bool, res)
end
function is_literal(node::SddNode)::Bool
    res = SddLibrary.sdd_node_is_literal(node.node)
    return convert(Bool, res)
end
function is_decision(node::SddNode)::Bool
    res = SddLibrary.sdd_node_is_decision(node.node)
    return convert(Bool, res)
end
function node_size(node::SddNode)::SddLibrary.SddNodeSize
    return SddLibrary.sdd_node_size(node.node)
end
function literal(node::SddNode)::SddLibrary.SddLiteral
    return SddLibrary.sdd_node_literal(node.node)
end
function elements(node::SddNode)::Array{PrimeSub,1}
    # TODO make abstract array for primesubs and avoid copying
    elements_ptr = SddLibrary.sdd_node_elements(node.node)
    m = size(node)
    primesubs = PrimeSub[]
    for i in 0:m - 1
        e = unsafe_load(elements_ptr + i)
        p = SddNode(e, node.manager)
        s = SddNode(e + 1, node.manager)
        push!(primesubs, PrimeSub(p, s))
    end
    return primesubs
end
function set_bit(bit::Int32, node::SddNode)
    SddLibrary.sdd_node_set_bit(bit, node.node)
end
function bit(node::SddNode)::Int32
    return SddLibrary.sdd_node_bit(node.node)
end

#
# # SDD FUNCTIONS

# SDD FILE I/O
# TODO make this safer
function read_sdd(filename::String, manager::SddManager)::SddNode
    node = SddLibrary.sdd_read(str_to_char(filename), manager.manager)
    return SddNode(node, manager.manager)
end
function save(filename::String, node::SddNode)
    SddLibrary.sdd_save(str_to_char(filename), node.node)
end
function save_as_dot(filename::String, node::SddNode)
    SddLibrary.sdd_save_as_dot(str_to_char(filename), node.node)
end
function shared_save_as_dot(filename::String, manager::SddManager)
    SddLibrary.sdd_shared_save_as_dot(str_to_char(filename), manager.manager)
end

# SDD SIZE AND NODE COUNT
# SDD
function count(node::SddNode)::SddLibrary.SddSize
    return SddLibrary.sdd_count(node.node)
end
function size(node::SddNode)::SddLibrary.SddSize
    return SddLibrary.sdd_size(node.node)
end
# SDD OF MANAGER
manager_size_fnames_c = [
    "size", "live_size", "dead_size",
    "count", "live_count", "dead_count"
]
for (fnj, fnc) in zip(manager_size_fnames_c, SddLibrary.manager_size_fnames_c)
    @eval begin
        function $(Symbol(fnj))(manager::SddManager)::SddLibrary.SddSize
            return ((SddLibrary).$(Symbol(fnc)))(manager.manager)
        end
    end
end
# SDD SIZE OF VTREE
vtree_size_fnames_j = [
    "size", "live_size", "dead_size",
    "size_at", "live_size_at", "dead_size_at",
    "size_above", "live_size_above", "dead_size_above",
    "count", "live_count", "dead_count",
    "count_at", "live_count_at", "dead_count_at",
    "count_above", "live_count_above", "dead_count_above"
]
for (fnj, fnc) in zip(vtree_size_fnames_j, SddLibrary.vtree_size_fnames_c)
    @eval begin
        function $(Symbol(fnj))(vtree::VTree)::SddLibrary.SddSize
            return ((SddLibrary).$(Symbol(fnc)))(vtree.vtree)
        end
    end
end

# CREATING VTREES
function vtree(var_count::Integer, vtree_type::String)::VTree
    vtree = SddLibrary.sdd_vtree_new(convert(SddLibrary.SddLiteral, var_count), str_to_char(vtree_type))
    return VTree(vtree)
end
function vtree(var_count::Integer, order::Array{<:Integer,1}, vtree_type::String; order_type::String="var_order")::VTree
    @assert order_type in Set(["var_order", "is_X_var"]) "$order_type not in a valid order type (var_order, is_X_var]"
    if order_type == "var_order"
        vtree = SddLibrary.sdd_vtree_new_with_var_order(convert(SddLibrary.SddLiteral, var_count), convert(Array{SddLibrary.SddLiteral,1}, order), str_to_char(vtree_type))
    elseif order_type == "is_X_var"
        vtree = SddLibrary.sdd_vtree_new_X_constrained(convert(SddLibrary.SddLiteral, var_count), convert(Array{SddLibrary.SddLiteral,1}, order), str_to_char(vtree_type))
    end
    return VTree(vtree)
end
function free(vtree::VTree)
    SddLibrary.sdd_vtree_free(vtree.vtree)
end

# VTREE FILE I/O
# TODO make this safer
function read_vtree(filename::String)::VTree
    vtree = SddLibrary.sdd_vtree_read(str_to_char(filename))
    return VTree(vtree)
end
function save(filename::String, vtree::VTree)
    SddLibrary.sdd_vtree_save(str_to_char(filename), vtree.vtree)
end
function save_as_dot(filename::String, vtree::VTree)
    SddLibrary.sdd_vtree_save_as_dot(str_to_char(filename), vtree.vtree)
end

# SDD MANAGER VTREE
function vtree(manager::SddManager; copy::Bool=false)::VTree
    if !copy
        vtree = SddLibrary.sdd_manager_vtree(manager.manager)
        return VTree(vtree)
    else
        vtree = SddLibrary.sdd_manager_vtree_copy(manager.manager)
        return VTree(vtree)
    end
end

# VTREE NAVIGATION
function left(vtree::VTree)::VTree
    vtree = SddLibrary.sdd_vtree_left(vtree.vtree)
    return VTree(vtree)
end
function right(vtree::VTree)::VTree
    vtree = SddLibrary.sdd_vtree_right(vtree.vtree)
    return VTree(vtree)
end
function parent(vtree::VTree)::VTree
    vtree = SddLibrary.sdd_vtree_parent(vtree.vtree)
    return VTree(vtree)
end

# VTREE FUNCTIONS
function is_leaf(vtree::VTree)::Bool
    return convert(Bool, SddLibrary.sdd_vtree_is_leaf(vtree.vtree))
end
function is_sub(vtree1::VTree, vtree2::VTree)::Bool
    return convert(Bool, SddLibrary.sdd_vtree_is_sub(vtree1.vtree, vtree2.vtree))
end
function lca(vtree1::VTree, vtree2::VTree)::VTree
    vtree =  SddLibrary.sdd_vtree_lca(vtree1.vtree, vtree2.vtree)
    return VTree(vtree)
end
function var_count(vtree::VTree)::SddLibrary.SddLiteral
    return SddLibrary.sdd_vtree_var_count(vtree.vtree)
end
function var(vtree::VTree)::SddLibrary.SddLiteral
    return SddLibrary.sdd_vtree_var(vtree.vtree)
end
function position(vtree::VTree)::SddLibrary.SddLiteral
    return SddLibrary.sdd_vtree_position(vtree.vtree)
end
# Vtree** sdd_vtree_location(Vtree* vtree, SddManager* manager);

# VTREE/SDD EDIT OPERATIONS
function rotate_left(vtree::VTree, manager::SddManager, limited::Union{Bool,Integer})::Bool
    return convert(Bool, SddLibrary.sdd_vtree_rotate_left(vtree.vtree, manager.manager, convert(Cint, limited)))
end
function rotate_right(vtree::VTree, manager::SddManager, limited::Union{Bool,Integer})::Bool
    return convert(Bool, SddLibrary.sdd_vtree_rotate_right(vtree.vtree, manager.manager, convert(Cint, limited)))
end
function swap(vtree::VTree, manager::SddManager, limited::Union{Bool,Integer})::Bool
    return convert(Bool, SddLibrary.sdd_vtree_swap(vtree.vtree, manager.manager, convert(Cint, limited)))
end

# LIMITS FOR VTREE/SDD EDIT OPERATIONS
function init_vtree_size_limit(vtree::VTree, manager::SddManager)
    SddLibrary.sdd_manager_init_vtree_size_limit(vtree.vtree, manager.manager)
end
function update_vtree_size_limit(manager::SddManager)
    SddLibrary.sdd_manager_update_vtree_size_limit(manager.manager)
end

# # VTREE STATE

# GARBAGE COLLECTION
function ref_count(node::SddNode)::SddLibrary.SddRefCount
    return SddLibrary.sdd_ref_count(node.node)
end
function ref(node::SddNode, manager::SddManager)::SddNode
    ref_node = SddLibrary.sdd_ref(node.node, manager.manager)
    return SddNode(ref_node, manager.manager)
end
function deref(node::SddNode, manager::SddManager)::SddNode
    ref_node = SddLibrary.sdd_deref(node.node, manager.manager)
    return SddNode(ref_node, manager.manager)
end
function garbage_collect(manager::SddManager)
    SddLibrary.sdd_manager_garbage_collect(manager.manager)
end
function garbage_collect(vtree::VTree, manager::SddManager)
    SddLibrary.sdd_manager_garbage_collect(vtree.vtree, manager.manager)
end
function garbage_collect_if(dead_node_threshold::Real, manager::SddManager)::Int32
    return SddLibrary.sdd_manager_garbage_collect_if(convert(Float32, dead_node_threshold), manager.manager)
end
function garbage_collect_if(dead_node_threshold::Real, vtree::VTree, manager::SddManager)::Int32
    return SddLibrary.sdd_manager_garbage_collect_if(convert(Float32, dead_node_threshold), vtree.vtree, manager.manager)
end

# MINIMIZATION
function minimize(manager::SddManager; limited::Bool=false)
    if !limited
        SddLibrary.sdd_manager_minimize(manager.manager)
    else
        SddLibrary.sdd_manager_minimize_limited(manager.manager)
    end
end
function minimize(vtree::VTree, manager::SddManager; limited::Bool=false)::VTree
    if !limited
        vtree = SddLibrary.sdd_vtree_minimize(vtree.vtree, manager.manager)
    else
        vtree = SddLibrary.sdd_vtree_minimize_limited(vtree.vtree, manager.manager)
    end
    return VTree(vtree)
end

function set_vtree_search_convergence_threshold(threshold::Real, manager::SddManager)
    SddLibrary.sdd_manager_set_vtree_search_convergence_threshold(convert(Float32, threshold), manager.manager)
end

function set_vtree_search_time_limit(time_limit::Real, manager::SddManager)
    SddLibrary.sdd_manager_set_vtree_search_time_limit(convert(Float32, time_limit), manager.manager)
end
function set_vtree_fragment_time_limit(time_limit::Real, manager::SddManager)
    SddLibrary.sdd_manager_set_vtree_fragment_time_limit(convert(Float32, time_limit), manager.manager)
end
function set_vtree_operation_time_limit(time_limit::Real, manager::SddManager)
    SddLibrary.sdd_manager_set_vtree_operation_time_limit(convert(Float32, time_limit), manager.manager)
end
function set_vtree_apply_time_limit(time_limit::Real, manager::SddManager)
    SddLibrary.sdd_manager_set_vtree_apply_time_limit(convert(Float32, time_limit), manager.manager)
end
function set_vtree_operation_memory_limit(memory_limit::Real, manager::SddManager)
    SddLibrary.sdd_manager_set_vtree_operation_memory_limit(convert(Float32, memory_limit), manager.manager)
end
function set_vtree_operation_size_limit(size_limit::Real, manager::SddManager)
    SddLibrary.sdd_manager_set_vtree_operation_size_limit(convert(Float32, size_limit), manager.manager)
end
function set_vtree_cartesian_product_limit(size_limit::Real, manager::SddManager)
    SddLibrary.sdd_manager_set_vtree_cartesian_product_limit(convert(Float32, size_limit), manager.manager)
end

# WMC
function wmc_manager(node::SddNode, log_mode::Bool, manager::SddManager)::WmcManager
    wmc = SddLibrary.wmc_manager_new(node.node, convert(Cint, log_mode), manager.manager)
    return WmcManager(wmc)
end
function free(manager::WmcManager)
    SddLibrary.wmc_manager_free(manager.manager)
end
function set_literal_weight(node::SddNode, weight::Real, manager::WmcManager)
    literal = SddLibrary.sdd_node_literal(node.node)
    SddLibrary.wmc_set_literal_weight(literal, convert(SddLibrary.SddWmc, weight), manager.manager)
end
function propagate(manager::WmcManager)::SddLibrary.SddWmc
    return SddLibrary.wmc_propagate(manager.manager)
end
function zero(manager::WmcManager)::SddLibrary.SddWmc
    return SddLibrary.wmc_zero_weight(manager.manager)
end
function one(manager::WmcManager)::SddLibrary.SddWmc
    return SddLibrary.wmc_one_weight(manager.manager)
end
function weight(literal::Integer, manager::WmcManager)::SddLibrary.SddWmc
    return SddLibrary.wmc_literal_weight(convert(SddLibrary.SddLiteral, literal), manager.manager)
end
function derivative(literal::Integer, manager::WmcManager)::SddLibrary.SddWmc
    return SddLibrary.wmc_literal_derivative(convert(SddLibrary.SddLiteral, literal), manager.manager)
end
function probability(literal::Integer, manager::WmcManager)::SddLibrary.SddWmc
    return SddLibrary.wmc_literal_pr(convert(SddLibrary.SddLiteral, literal), manager.manager)
end

# CONVENIENCE METHODS
function conjoin(node1::SddNode, node2::SddNode)::SddNode
    node = SddLibrary.sdd_conjoin(node1.node, node2.node, node1.manager)
    return SddNode(node, node1.manager)
end
function disjoin(node1::SddNode, node2::SddNode)::SddNode
    node = SddLibrary.sdd_disjoin(node1.node, node2.node, node1.manager)
    return SddNode(node, node1.manager)
end
function negate(node::SddNode)::SddNode
    nodeptr = SddLibrary.sdd_negate(node.node, node.manager)
    return SddNode(nodeptr, node.manager)
end
function equiv(left::SddNode, right::SddNode)::SddNode
    return (!left | right) & (left | !right)
end
function model_count(node::SddNode; globally::Bool=false)::SddLibrary.SddModelCount
    if !globally
        return SddLibrary.sdd_model_count(node.node, node.manager)
    else
        return SddLibrary.sdd_global_model_count(node.node, node.manager)
    end
end
function ref(node::SddNode)::SddNode
    ref_node = SddLibrary.sdd_ref(node.node, node.manager)
    return SddNode(ref_node, node.manager)
end
function deref(node::SddNode)::SddNode
    ref_node = SddLibrary.sdd_deref(node.node, node.manager)
    return SddNode(ref_node, node.manager)
end

function dot(filename::String, structure::Union{VTree,SddNode})
    save_as_dot(filename, structure)
end
function wmc_manager(node::SddNode; log_mode::Bool=true)::WmcManager
    wmc = SddLibrary.wmc_manager_new(node.node, convert(Cint, log_mode), node.manager)
    return WmcManager(wmc)
end

Base.:&(node1::SddNode, node2::SddNode) = conjoin(node1, node2)
Base.:|(node1::SddNode, node2::SddNode) = disjoin(node1, node2)
Base.:~(node::SddNode) = negate(node)
↔(left::SddNode, right::SddNode) = equiv(left, right)

# FNF methods
@with_kw struct CompilerOptions
    vtree_search_mode::Int32 = -1
    post_search::Bool = false
    verbose::Bool = false
end

function read_cnf(filename::String, manager::SddManager; compiler_options=CompilerOptions())::SddNode
    cnf = read_cnf(filename)
    sdd_node = fnf_to_sdd(cnf, manager.manager, compiler_options)
    return SddNode(sdd_node, manager.manager)
end
function read_dnf(filename::String, manager::SddManager; compiler_options=CompilerOptions())::SddNode
    dnf = read_dnf(filename)
    sdd_node = fnf_to_sdd(dnf, manager.manager, compiler_options)
    return SddNode(sdd_node, manager.manager)
end

export ↔
end


text: include("XtbML.jl")