text:
import Tokenize.Tokens: isoperator
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) = error()#precedence(kindof(x))

isoperator(x) = false
isoperator(t::AbstractToken) = isoperator(kindof(t))

isunaryop(op) = false
isunaryop(op::EXPR) = isoperator(op) && ((valof(op) == "<:" ||
                                          valof(op) == ">:" ||
                                          valof(op) == "+" ||
                                          valof(op) == "-" ||
                                          valof(op) == "!" ||
                                          valof(op) == "~" ||
                                          valof(op) == "¬" ||
                                          valof(op) == "&" ||
                                          valof(op) == "√" ||
                                          valof(op) == "∛"  ||
                                          valof(op) == "∜"  ||
                                          valof(op) == "::" ||
                                          valof(op) == "\$" ||
                                          valof(op) == ":" ||
                                          valof(op) == "⋆") ||
                        (length(valof(op)) == 2 && valof(op)[1] == '.' && (valof(op)[2] == '+' ||
                                                                           valof(op)[2] == '-' ||
                                                                           valof(op)[2] == '!' ||
                                                                           valof(op)[2] == '~' ||
                                                                           valof(op)[2] == '¬' ||
                                                                           valof(op)[2] == '√' ||
                                                                           valof(op)[2] == '∛' ||
                                                                           valof(op)[2] == '∜' ||
                                                                           valof(op)[2] == '⋆')))
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) && !(valof(op) == "√" ||
    valof(op) == "∛" ||
    valof(op) == "∜" ||
    valof(op) == "!" ||
    valof(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
isdotted(x::EXPR) = valof(x) isa String && length(valof(x)) > 1 && valof(x)[1] === '.' && !(valof(x) == "..." || valof(x) == "..")

issyntaxcall(op) = false
function issyntaxcall(op::EXPR)
    v = valof(op)
    if v[1] == '.' && length(v) > 1 && assign_prec(v[2:end]) && v[2] !== '~'
        return true
    end
    assign_prec(v) && !(v == "~" || v == ".~" || v == "=>") ||
    v == "-->" ||
    v == "||" ||
    v == "&&" ||
    v == "<:" ||
    v == ">:" ||
    v == ":" ||
    v == "::" ||
    v == "." ||
    v == "..." ||
    v == "'" ||
    v == "where" ||
    v == "->"
end

issyntaxunarycall(op) = false
function issyntaxunarycall(op::EXPR)
    v = valof(op)
    !isdotted(op) && (v == "\$" || v == "&" || v == "::" || v == "..." || v == "'" || v == "<:" || v == ">:")
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)
    if is_colon(op)
        ret = parse_unary_colon(ps, op)
    elseif should_negate_number_literal(ps, op)
        arg = mLITERAL(next(ps))
        ret = EXPR(literalmap(kindof(ps.t)), op.fullspan + arg.fullspan, (op.fullspan + arg.span), string(is_plus(op) ? "+" : "-", val(ps.t, ps)))
    else
        prec = valof(op) == "::" ? DeclarationOp :
                valof(op) == "&" ? DeclarationOp :
                valof(op) == "\$" ? 20 : PowerOp
        arg = @closer ps :unary @precedence ps prec parse_expression(ps)
        if issyntaxunarycall(op)
            ret = EXPR(op, EXPR[arg], nothing)
        else
            ret = EXPR(:call, EXPR[op, arg], nothing)
        end
    end
    return ret
end

function parse_unary_colon(ps::ParseState, op::EXPR)
    op = requires_no_ws(op, ps)
    if Tokens.iskeyword(kindof(ps.nt)) && op.span == op.fullspan
        ret = EXPR(:quotenode, EXPR[EXPR(:IDENTIFIER, next(ps))], EXPR[op])
    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], nothing)
        end
        ret = EXPR(:quotenode, EXPR[id], EXPR[op])
    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[INSTANCE(next(ps))], EXPR[op])
    elseif closer(ps)
        ret = op
    else
        prev_errored = ps.errored
        arg = @precedence ps 20 @nocloser ps :inref parse_expression(ps)
        if isbracketed(arg) && headof(arg.args[1]) === :errortoken && errorof(arg.args[1]) === UnexpectedAssignmentOp
            ps.errored = prev_errored
            arg.args[1] = arg.args[1].args[1]
            setparent!(arg.args[1], arg)
        end
        # TODO: need special conversion where arg is a n-bracketed terminal (not keywords)
        unwrapped = unwrapbracket(arg)
        if isoperator(unwrapped) || isidentifier(unwrapped) || isliteral(unwrapped)
            ret = EXPR(:quotenode, EXPR[arg], EXPR[op])
        elseif arg.head == :tuple && length(arg.args) == 1 && arg.args[1].head == :parameters && arg.args[1].args !== nothing && length(arg.args[1].args) == 0 && arg.span > 3
            block = EXPR(:BLOCK, EXPR[], EXPR[])
            block.span = arg.span
            block.fullspan = arg.fullspan
            ret = EXPR(:quote, EXPR[block], EXPR[op])
        else
            ret = EXPR(:quote, EXPR[arg], EXPR[op])
        end
    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], nothing)
    end
    if issyntaxcall(op)
        ret = EXPR(op, EXPR[ret, nextarg], nothing)
    else
        ret = EXPR(:call, EXPR[op, ret, nextarg], nothing)
    end
    return ret
end

# Parse conditionals

isconditional(x::EXPR) = headof(x) === :if && hastrivia(x) && isoperator(first(x.trivia))
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 kindof(ps.nt) !== Tokens.COLON
        op2 = mErrorToken(ps, EXPR(:OPERATOR, 0, 0, ":"), MissingColon)
        nextarg2 = mErrorToken(ps, Unknown)
        return EXPR(:if, EXPR[ret, nextarg, nextarg2], EXPR[op, op2])
    else
        op2 = requires_ws(EXPR(:OPERATOR, next(ps)), ps)
    end

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

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

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

    if headof(ret) === :comparison
        push!(ret, op)
        push!(ret, nextarg)
    elseif can_become_comparison(ret)
        if isoperator(headof(ret))
            ret = EXPR(:comparison, EXPR[ret.args[1], ret.head, ret.args[2], op, nextarg], nothing)
        else
            ret = EXPR(:comparison, EXPR[ret.args[2], ret.args[1], ret.args[3], op, nextarg], nothing)
        end
    elseif issyntaxcall(op)
        ret = EXPR(op, EXPR[ret, nextarg], nothing)
    else
        ret = EXPR(:call, EXPR[op, ret, nextarg], nothing)
    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[1])
        ret.trivia = EXPR[]
        pushtotrivia!(ret, op)
        push!(ret, nextarg)
    else
        ret = EXPR(:call, EXPR[op, ret, nextarg], nothing)
    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)
        # TODO: this smells wrong
        nextarg = EXPR(:call, EXPR[op, ret.args[2], nextarg], nothing)
        ret = EXPR(:call, EXPR[ret.args[1], nextarg], nothing)
    else
        ret = EXPR(:call, EXPR[op, ret, nextarg], nothing)
    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 headof(nextarg) === :braces
        pushfirst!(nextarg.args, ret)
        pushfirst!(nextarg.trivia, op)
        ret = EXPR(:where, nextarg.args,nextarg.trivia)
    else
        ret = EXPR(:where, EXPR[ret, nextarg], EXPR[op])
    end
    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(:tuple, sig.args[2:end], sig.trivia)
            if iserred
                nextarg = mErrorToken(ps, nextarg, UnexpectedWhiteSpace)
            end
        else
            sig = @default ps parse_call(ps, ret)
            nextarg = EXPR(:tuple, sig.args[2:end], sig.trivia)
        end
    elseif iskeyword(ps.nt) || both_symbol_and_op(ps.nt)
        nextarg = EXPR(:IDENTIFIER, next(ps))
    elseif kindof(ps.nt) === Tokens.COLON
        op2 = EXPR(:OPERATOR, next(ps))
        if kindof(ps.nt) === Tokens.LPAREN
            nextarg = @closeparen ps @precedence ps DotOp - LtoR(DotOp) parse_expression(ps)
            nextarg = EXPR(:quotenode, EXPR[nextarg], EXPR[op2])
        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 = EXPR(:OPERATOR, next(ps))
        nextarg = parse_call(ps, op2)
    else
        nextarg = @precedence ps DotOp - LtoR(DotOp) parse_expression(ps)
    end

    if isidentifier(nextarg) || isinterpolant(nextarg)
        ret = EXPR(op, EXPR[ret, EXPR(:quotenode, EXPR[nextarg], nothing)], nothing)
    elseif headof(nextarg) === :vect || headof(nextarg) === :braces
        ret = EXPR(op, EXPR[ret, EXPR(:quote, EXPR[nextarg], nothing)], nothing)
    elseif headof(nextarg) === :macrocall
        # TODO : ?
        mname = EXPR(op, EXPR[ret, EXPR(:quotenode, EXPR[nextarg.args[1]], nothing)], nothing)
        ret = EXPR(:macrocall, EXPR[mname], nextarg.trivia)
        for i = 2:length(nextarg.args)
            push!(ret, nextarg.args[i])
        end
    else
        ret = EXPR(op, EXPR[ret, nextarg], nothing)
    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], nothing)
    end
    return EXPR(op, EXPR[ret, arg], nothing)
end

function parse_operator(ps::ParseState, ret::EXPR, op::EXPR)
    P = isdotted(op) ? get_prec(valof(op)[2:end]) : get_prec(valof(op))
    if op.val == "*" && op.fullspan == 0 # implicit multiplication has a very high precedence, but lower than ^
        P = RationalOp
    end
    if headof(ret) === :call && (is_plus(ret.args[1]) || is_star(ret.args[1])) &&
          valof(ret.args[1]) == valof(op) && ret.args[1].span > 0 &&
          !(hastrivia(ret) && headof(ret[end]) === :RPAREN)
        # a + b -> a + b + c
        nextarg = @precedence ps P - LtoR(P) parse_expression(ps)
        !hastrivia(ret) && (ret.trivia = EXPR[])
        pushtotrivia!(ret, op)
        push!(ret, nextarg)
        ret = ret
    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)
        ret = EXPR(op, EXPR[ret], nothing)
    elseif is_prime(op)
        if isidentifier(ret) || isliteral(ret) ||
                headof(ret) in (:call, :tuple, :brackets, :ref, :vect, :vcat, :hcat, :typed_vcat, :typed_hcat, :comprehension, :typed_comprehension, :curly, :braces, :braces_cat) ||
                headof(ret) === :do ||
                is_dot(headof(ret)) ||
                is_prime(headof(ret))
            if valof(op) == "'"
                ret = EXPR(op, EXPR[ret], nothing)
            else
                ret = EXPR(:call, EXPR[op, ret], nothing)
            end
        else
            ret = EXPR(:errortoken, EXPR[ret, op])
            ret.meta = UnexpectedToken
        end
    elseif P == ComparisonOp
        ret = parse_comp_operator(ps, ret, op)
    elseif P == PowerOp
        ret = parse_operator_power(ps, ret, op)
    else
        ltor = valof(op) == "<|" ? true : LtoR(P)
        nextarg = @precedence ps P - ltor parse_expression(ps)
        if issyntaxcall(op)
            ret = EXPR(op, EXPR[ret, nextarg], nothing)
        else
            ret = EXPR(:call, EXPR[op, ret, nextarg], nothing)
        end
    end
    return ret
end

assign_prec(op::String) = get(AllPrecs, op, 0) === AssignmentOp || (length(op) > 1 && op[1] == '.' && assign_prec(op[2:end]))
comp_prec(op::String) = get(AllPrecs, op, 0) === ComparisonOp || (length(op) > 1 && op[1] == '.' && comp_prec(op[2:end]))

# Which ops have different precend when dotted?
# [op for op in Symbol.(keys(CSTParser.AllPrecs)) if Base.isoperator(op) && Base.isoperator(Symbol(".", op)) && Base.operator_precedence(op) !== Base.operator_precedence(Symbol(".", op))]

get_prec(op) = AllPrecs[maybe_strip_suffix(op)]

function maybe_strip_suffix(op::String)
    for (i, c) in enumerate(op)
        if Tokenize.Lexers.isopsuffix(c)
            return op[1:prevind(op, i)]
        end
    end
    return op
end

const AllPrecs = Dict(
    "≏" => PlusOp,
    "⤎" => ArrowOp,
    "<:" => ComparisonOp,
    "⫌" => ComparisonOp,
    "⇀" => 0,
    "≁" => ComparisonOp,
    "⧀" => ComparisonOp,
    "⪎" => ComparisonOp,
    "⪱" => ComparisonOp,
    "↬" => 0,
    "⇝" => 0,
    "⩳" => ComparisonOp,
    "⪻" => ComparisonOp,
    "//" => RationalOp,
    "⬿" => ArrowOp,
    "⊅" => ComparisonOp,
    "⋳" => ComparisonOp,
    "≮" => ComparisonOp,
    "\$=" => AssignmentOp,
    "≗" => ComparisonOp,
    "⋗" => ComparisonOp,
    "⩊" => PlusOp,
    "⪵" => ComparisonOp,
    "⋟" => ComparisonOp,
    "⤠" => ArrowOp,
    "≵" => ComparisonOp,
    "⧁" => ComparisonOp,
    "≫" => ComparisonOp,
    "⩫" => ComparisonOp,
    "⪸" => ComparisonOp,
    "⊓" => TimesOp,
    "⭄" => ArrowOp,
    "⇁" => 0,
    "≀" => TimesOp,
    "⪚" => ComparisonOp,
    "∷" => ComparisonOp,
    "⥪" => ArrowOp,
    "⥔" => PowerOp,
    "⊻" => PlusOp,
    "⩯" => ComparisonOp,
    "↷" => 0,
    "⇿" => ArrowOp,
    "⪉" => ComparisonOp,
    "⊰" => ComparisonOp,
    "⨟" => 0,
    "⩻" => ComparisonOp,
    "⭋" => ArrowOp,
    "⩒" => PlusOp,
    "-" => PlusOp,
    "⫸" => ComparisonOp,
    "⤒" => PowerOp,
    "+=" => AssignmentOp,
    "↦" => ArrowOp,
    "⊜" => ComparisonOp,
    "|>" => PipeOp,
    "⩀" => TimesOp,
    "⪫" => ComparisonOp,
    "≡" => ComparisonOp,
    "⋹" => ComparisonOp,
    "⤑" => ArrowOp,
    "⥐" => ArrowOp,
    "≭" => ComparisonOp,
    "≼" => ComparisonOp,
    "⥇" => ArrowOp,
    "⋶" => ComparisonOp,
    "⨸" => TimesOp,
    "⩽" => ComparisonOp,
    "⪐" => ComparisonOp,
    "⨥" => PlusOp,
    "∊" => ComparisonOp,
    "⪙" => ComparisonOp,
    "⥢" => ArrowOp,
    "⅋" => TimesOp,
    "⫘" => ComparisonOp,
    "⊂" => ComparisonOp,
    "⩸" => ComparisonOp,
    "⩄" => TimesOp,
    "⤞" => ArrowOp,
    "⫍" => ComparisonOp,
    "⨢" => PlusOp,
    "⊄" => ComparisonOp,
    "≯" => ComparisonOp,
    "===" => ComparisonOp,
    "⬺" => ArrowOp,
    "⪛" => ComparisonOp,
    "⋡" => ComparisonOp,
    "⫀" => ComparisonOp,
    "⨣" => PlusOp,
    "↞" => 0,
    "⟾" => ArrowOp,
    "⊬" => ComparisonOp,
    "⇌" => 0,
    "≷" => ComparisonOp,
    "∓" => PlusOp,
    "⪶" => ComparisonOp,
    "⊀" => ComparisonOp,
    "⊛" => TimesOp,
    "⋠" => ComparisonOp,
    "↢" => 0,
    "⥦" => ArrowOp,
    "⪪" => ComparisonOp,
    "⋌" => TimesOp,
    "⪅" => ComparisonOp,
    "≿" => ComparisonOp,
    ">>=" => AssignmentOp,
    "\$" => PlusOp,
    "≽" => ComparisonOp,
    "⪟" => ComparisonOp,
    "⋆" => TimesOp,
    "±" => PlusOp,
    "⋴" => ComparisonOp,
    "⥧" => ArrowOp,
    "/" => TimesOp,
    "⥕" => PowerOp,
    "." => DotOp,
    "⩢" => PlusOp,
    "⟖" => TimesOp,
    "⪗" => ComparisonOp,
    "≓" => ComparisonOp,
    "⋨" => ComparisonOp,
    "⥜" => PowerOp,
    "⪊" => ComparisonOp,
    "⨽" => TimesOp,
    "⪦" => ComparisonOp,
    "≱" => ComparisonOp,
    "⧷" => TimesOp,
    ">>" => BitShiftOp,
    "⊏" => ComparisonOp,
    "⬴" => ArrowOp,
    "⪴" => ComparisonOp,
    "⁝" => ColonOp,
    "⤘" => ArrowOp,
    "⩜" => TimesOp,
    "￬" => PowerOp,
    "∤" => TimesOp,
    "⩲" => ComparisonOp,
    "⫙" => ComparisonOp,
    "⇵" => PowerOp,
    "⤗" => ArrowOp,
    "⪀" => ComparisonOp,
    "⩬" => ComparisonOp,
    "↫" => 0,
    "&" => TimesOp,
    "⩺" => ComparisonOp,
    "↚" => ArrowOp,
    "⨇" => TimesOp,
    "⋭" => ComparisonOp,
    "≐" => ComparisonOp,
    "∩" => TimesOp,
    "↓" => PowerOp,
    "⪓" => ComparisonOp,
    "≦" => ComparisonOp,
    "=" => AssignmentOp,
    "∧" => TimesOp,
    "⇇" => 0,
    "≖" => ComparisonOp,
    "⇸" => ArrowOp,
    "⋙" => ComparisonOp,
    "≚" => ComparisonOp,
    "⇴" => ArrowOp,
    "⬲" => ArrowOp,
    "⥙" => PowerOp,
    "⋢" => ComparisonOp,
    "⫹" => ComparisonOp,
    "⩛" => PlusOp,
    "⋿" => ComparisonOp,
    "⪑" => ComparisonOp,
    "⭉" => ArrowOp,
    "⋎" => PlusOp,
    "⋋" => TimesOp,
    "⇋" => 0,
    "≴" => ComparisonOp,
    "…" => ColonOp,
    "⪰" => ComparisonOp,
    "↺" => 0,
    "⪖" => ComparisonOp,
    "⨦" => PlusOp,
    "⬼" => ArrowOp,
    "⤐" => ArrowOp,
    "⊗" => TimesOp,
    "⨮" => PlusOp,
    "⊚" => TimesOp,
    "⥈" => ArrowOp,
    "≣" => ComparisonOp,
    "⋫" => ComparisonOp,
    "≺" => ComparisonOp,
    "⇽" => ArrowOp,
    "⬳" => ArrowOp,
    "⥨" => ArrowOp,
    "⋽" => ComparisonOp,
    "⩟" => TimesOp,
    "⋞" => ComparisonOp,
    "⩖" => PlusOp,
    "⧥" => ComparisonOp,
    "⨝" => TimesOp,
    "⪕" => ComparisonOp,
    "⊟" => PlusOp,
    "⋑" => ComparisonOp,
    "⊲" => ComparisonOp,
    "≟" => ComparisonOp,
    "⫈" => ComparisonOp,
    "⥋" => ArrowOp,
    "⩂" => PlusOp,
    "⥏" => PowerOp,
    "⥩" => ArrowOp,
    "⥑" => PowerOp,
    "⪤" => ComparisonOp,
    "⫖" => ComparisonOp,
    "×" => TimesOp,
    "÷=" => AssignmentOp,
    "⋪" => ComparisonOp,
    "⫓" => ComparisonOp,
    "⫄" => ComparisonOp,
    "⊋" => ComparisonOp,
    "⇆" => 0,
    "⤌" => ArrowOp,
    "⟹" => ArrowOp,
    "⊖" => PlusOp,
    "⫏" => ComparisonOp,
    "⫷" => ComparisonOp,
    "⪜" => ComparisonOp,
    "⊕" => PlusOp,
    "%=" => AssignmentOp,
    "⊍" => TimesOp,
    "⨤" => PlusOp,
    "≃" => ComparisonOp,
    "⩪" => ComparisonOp,
    "⤉" => PowerOp,
    "⪽" => ComparisonOp,
    "⋓" => PlusOp,
    "⤀" => ArrowOp,
    "∺" => ComparisonOp,
    "⤆" => ArrowOp,
    "⊳" => ComparisonOp,
    "//=" => AssignmentOp,
    "≒" => ComparisonOp,
    "⩦" => ComparisonOp,
    "↶" => 0,
    "⊡" => TimesOp,
    "⟈" => ComparisonOp,
    "≩" => ComparisonOp,
    "||" => LazyOrOp,
    "↠" => ArrowOp,
    "⥒" => ArrowOp,
    "⨧" => PlusOp,
    "⪧" => ComparisonOp,
    "⨪" => PlusOp,
    "⇹" => ArrowOp,
    "⥛" => ArrowOp,
    "↣" => ArrowOp,
    "⫁" => ComparisonOp,
    "⭊" => ArrowOp,
    "⋜" => ComparisonOp,
    "⪳" => ComparisonOp,
    "⋼" => ComparisonOp,
    "⥍" => PowerOp,
    "⫎" => ComparisonOp,
    "^=" => AssignmentOp,
    "⪢" => ComparisonOp,
    "∪" => PlusOp,
    "⩅" => PlusOp,
    "⬵" => ArrowOp,
    "⬽" => ArrowOp,
    "≬" => ComparisonOp,
    "⧶" => TimesOp,
    "⪡" => ComparisonOp,
    "≶" => ComparisonOp,
    "⇼" => ArrowOp,
    "^" => PowerOp,
    "⇷" => ArrowOp,
    "⋸" => ComparisonOp,
    "⊩" => ComparisonOp,
    "⭀" => ArrowOp,
    "⧴" => ArrowOp,
    "⊷" => ComparisonOp,
    "⤔" => ArrowOp,
    "~" => AssignmentOp,
    "⫒" => ComparisonOp,
    "⩼" => ComparisonOp,
    "⧺" => PlusOp,
    "⧻" => PlusOp,
    ">>>=" => AssignmentOp,
    "⪂" => ComparisonOp,
    "⨻" => TimesOp,
    "⤏" => ArrowOp,
    "⥘" => PowerOp,
    "≕" => ComparisonOp,
    "⥓" => ArrowOp,
    "⨰" => TimesOp,
    "⨺" => PlusOp,
    "⟵" => ArrowOp,
    "⇜" => 0,
    "<<=" => AssignmentOp,
    "≔" => ComparisonOp,
    "⩋" => TimesOp,
    "⫊" => ComparisonOp,
    "⪠" => ComparisonOp,
    "⥬" => ArrowOp,
    "⥖" => ArrowOp,
    "⪣" => ComparisonOp,
    "↔" => ArrowOp,
    "≥" => ComparisonOp,
    "⧤" => ComparisonOp,
    "⤄" => ArrowOp,
    "⊢" => ComparisonOp,
    "⋛" => ComparisonOp,
    "⩔" => PlusOp,
    "⇶" => ArrowOp,
    "⤇" => ArrowOp,
    "⟽" => ArrowOp,
    "⩵" => ComparisonOp,
    "⬶" => ArrowOp,
    "⩌" => PlusOp,
    "↮" => ArrowOp,
    "⪾" => ComparisonOp,
    "⋕" => ComparisonOp,
    "⥤" => ArrowOp,
    "∋" => ComparisonOp,
    "↽" => 0,
    "⋇" => TimesOp,
    "⪍" => ComparisonOp,
    "⩃" => TimesOp,
    "⨈" => PlusOp,
    "⋚" => ComparisonOp,
    "⩐" => PlusOp,
    "≾" => ComparisonOp,
    "⋤" => ComparisonOp,
    "∔" => PlusOp,
    "⤈" => PowerOp,
    "⇾" => ArrowOp,
    "⊮" => ComparisonOp,
    "⪩" => ComparisonOp,
    "⬷" => ArrowOp,
    "|++|" => PlusOp,
    "↪" => 0,
    "⭇" => ArrowOp,
    "⨨" => PlusOp,
    "⩓" => TimesOp,
    "⤊" => PowerOp,
    "↤" => 0,
    "⊱" => ComparisonOp,
    "⟼" => ArrowOp,
    "⨴" => TimesOp,
    "⤋" => PowerOp,
    "⊴" => ComparisonOp,
    "≊" => ComparisonOp,
    "⧡" => ComparisonOp,
    "⫆" => ComparisonOp,
    "⋯" => ColonOp,
    "≜" => ComparisonOp,
    "∍" => ComparisonOp,
    "::" => DeclarationOp,
    "⥣" => PowerOp,
    "⪃" => ComparisonOp,
    "⪹" => ComparisonOp,
    "⩴" => ComparisonOp,
    "⥫" => ArrowOp,
    "￫" => ArrowOp,
    ">>>" => BitShiftOp,
    "≢" => ComparisonOp,
    "⊈" => ComparisonOp,
    "≤" => ComparisonOp,
    "⟿" => ArrowOp,
    "⪨" => ComparisonOp,
    "⋖" => ComparisonOp,
    "⫉" => ComparisonOp,
    "⩹" => ComparisonOp,
    "￪" => PowerOp,
    "⊆" => ComparisonOp,
    "!==" => ComparisonOp,
    "⪔" => ComparisonOp,
    "⥎" => ArrowOp,
    "≋" => ComparisonOp,
    "⊁" => ComparisonOp,
    "⋥" => ComparisonOp,
    "⋉" => TimesOp,
    "⨱" => TimesOp,
    "?" => ConditionalOp,
    "⥚" => ArrowOp,
    "⋷" => ComparisonOp,
    "⟒" => ComparisonOp,
    "⫅" => ComparisonOp,
    "⪏" => ComparisonOp,
    "⩕" => TimesOp,
    "⪞" => ComparisonOp,
    ">:" => ComparisonOp,
    "⟗" => TimesOp,
    "⪷" => ComparisonOp,
    "⩎" => TimesOp,
    "⬱" => ArrowOp,
    "⇺" => ArrowOp,
    "⩁" => PlusOp,
    "⥆" => ArrowOp,
    "⊘" => TimesOp,
    "⨼" => TimesOp,
    "⊙" => TimesOp,
    "⊇" => ComparisonOp,
    "⋍" => ComparisonOp,
    "≛" => ComparisonOp,
    "⊒" => ComparisonOp,
    "⋏" => TimesOp,
    "⭁" => ArrowOp,
    "⤂" => ArrowOp,
    "⥄" => ArrowOp,
    "⋝" => ComparisonOp,
    "÷" => TimesOp,
    "⧣" => ComparisonOp,
    "⊞" => PlusOp,
    "⬾" => ArrowOp,
    "⇉" => 0,
    "⨶" => TimesOp,
    "in" => ComparisonOp,
    "==" => ComparisonOp,
    "⟶" => ArrowOp,
    "⬰" => ArrowOp,
    ":=" => AssignmentOp,
    "⪋" => ComparisonOp,
    "≂" => PlusOp,
    "<|" => PipeOp,
    "⥊" => ArrowOp,
    "⪺" => ComparisonOp,
    "≎" => ComparisonOp,
    "≈" => ComparisonOp,
    "∌" => ComparisonOp,
    "⨹" => PlusOp,
    "⪘" => ComparisonOp,
    "⨵" => TimesOp,
    "⦾" => TimesOp,
    "⇒" => ArrowOp,
    ".." => ColonOp,
    "⇛" => 0,
    "⪿" => ComparisonOp,
    "≹" => ComparisonOp,
    "⋱" => ColonOp,
    "⇎" => ArrowOp,
    "→" => ArrowOp,
    "⥮" => PowerOp,
    "-=" => AssignmentOp,
    "<" => ComparisonOp,
    "⪇" => ComparisonOp,
    "⥗" => ArrowOp,
    "⪈" => ComparisonOp,
    "⨫" => PlusOp,
    "⥞" => ArrowOp,
    "⭂" => ArrowOp,
    "⩠" => TimesOp,
    "\\" => TimesOp,
    "⩿" => ComparisonOp,
    "↑" => PowerOp,
    "≳" => ComparisonOp,
    "⋧" => ComparisonOp,
    "⩧" => ComparisonOp,
    "--" => 0,
    "⪮" => ComparisonOp,
    "⪄" => ComparisonOp,
    "⊎" => PlusOp,
    "⩾" => ComparisonOp,
    "↝" => 0,
    "⨲" => TimesOp,
    "⊐" => ComparisonOp,
    "⟻" => ArrowOp,
    "⫗" => ComparisonOp,
    "⩞" => TimesOp,
    "⥯" => PowerOp,
    "⪌" => ComparisonOp,
    "≨" => ComparisonOp,
    "⋄" => TimesOp,
    "⪯" => ComparisonOp,
    "⦸" => TimesOp,
    "⫐" => ComparisonOp,
    "↜" => 0,
    "≧" => ComparisonOp,
    "∉" => ComparisonOp,
    "⤟" => ArrowOp,
    "⋾" => ComparisonOp,
    "￩" => ArrowOp,
    "<<" => BitShiftOp,
    "≞" => ComparisonOp,
    "⩍" => TimesOp,
    "⇄" => 0,
    "▷" => TimesOp,
    "⋘" => ComparisonOp,
    "≇" => ComparisonOp,
    "⇏" => ArrowOp,
    "⊶" => ComparisonOp,
    "⫇" => ComparisonOp,
    "*" => TimesOp,
    "⩣" => PlusOp,
    "⬸" => ArrowOp,
    "⥥" => PowerOp,
    "≙" => ComparisonOp,
    "*=" => AssignmentOp,
    "≄" => ComparisonOp,
    "⊑" => ComparisonOp,
    "⫑" => ComparisonOp,
    "⩗" => PlusOp,
    "⟷" => ArrowOp,
    "⟱" => PowerOp,
    "≻" => ComparisonOp,
    "⋲" => ComparisonOp,
    "≑" => ComparisonOp,
    "⭈" => ArrowOp,
    "∾" => ComparisonOp,
    "&=" => AssignmentOp,
    "≪" => ComparisonOp,
    "\\\\=" => TimesOp,
    "↼" => 0,
    "⬹" => ArrowOp,
    "∗" => TimesOp,
    "⊣" => ComparisonOp,
    "⫺" => ComparisonOp,
    "⥌" => PowerOp,
    "⋩" => ComparisonOp,
    "≲" => ComparisonOp,
    "≝" => ComparisonOp,
    "⊽" => PlusOp,
    "/=" => AssignmentOp,
    "⋬" => ComparisonOp,
    "⩶" => ComparisonOp,
    "⫛" => TimesOp,
    "=>" => ConditionalOp,
    "⪁" => ComparisonOp,
    "⨩" => PlusOp,
    "⥝" => PowerOp,
    "⩝" => PlusOp,
    "⟕" => TimesOp,
    "⤖" => ArrowOp,
    "⥉" => PowerOp,
    "⩏" => PlusOp,
    "!=" => ComparisonOp,
    "⨳" => TimesOp,
    "≰" => ComparisonOp,
    "⟂" => ComparisonOp,
    "⋐" => ComparisonOp,
    "\\=" => AssignmentOp,
    "⤝" => ArrowOp,
    "⦷" => ComparisonOp,
    "⊔" => PlusOp,
    ">" => ComparisonOp,
    "⥰" => ArrowOp,
    "∽" => ComparisonOp,
    "⋦" => ComparisonOp,
    "⤅" => ArrowOp,
    "⫕" => ComparisonOp,
    "⋺" => ComparisonOp,
    "∈" => ComparisonOp,
    "∙" => TimesOp,
    "⥅" => ArrowOp,
    "⦼" => TimesOp,
    "<=" => ComparisonOp,
    "←" => ArrowOp,
    "isa" => ComparisonOp,
    "⋻" => ComparisonOp,
    "⫔" => ComparisonOp,
    "⋵" => ComparisonOp,
    "⋅" => TimesOp,
    "⨷" => TimesOp,
    "⟰" => PowerOp,
    "∨" => PlusOp,
    "≆" => ComparisonOp,
    "⪥" => ComparisonOp,
    "≍" => ComparisonOp,
    "⪝" => ComparisonOp,
    "⫂" => ComparisonOp,
    "⊻=" => AssignmentOp,
    "⨭" => PlusOp,
    "⋒" => TimesOp,
    "⦿" => TimesOp,
    "≸" => ComparisonOp,
    "⥡" => PowerOp,
    "+" => PlusOp,
    "⇻" => ArrowOp,
    "⥟" => ArrowOp,
    "⟉" => ComparisonOp,
    "⊵" => ComparisonOp,
    "⩰" => ComparisonOp,
    "%" => TimesOp,
    "⫃" => ComparisonOp,
    "∸" => PlusOp,
    "⋣" => ComparisonOp,
    "-->" => ArrowOp,
    "<--" => ArrowOp,
    "<-->" => ArrowOp,
    "∥" => ComparisonOp,
    "∦" => ComparisonOp,
    "⋮" => ColonOp,
    "⪒" => ComparisonOp,
    "⤃" => ArrowOp,
    "⟑" => TimesOp,
    "⩑" => TimesOp,
    "⋰" => ColonOp,
    "⪬" => ComparisonOp,
    "≠" => ComparisonOp,
    "≅" => ComparisonOp,
    "⊼" => TimesOp,
    "⬻" => ArrowOp,
    "⩘" => TimesOp,
    "⭃" => ArrowOp,
    "⥭" => ArrowOp,
    "⤕" => ArrowOp,
    "≘" => ComparisonOp,
    "⇠" => 0,
    "⟺" => ArrowOp,
    "⪼" => ComparisonOp,
    "⇢" => 0,
    "↩" => 0,
    "⋊" => TimesOp,
    "&&" => LazyAndOp,
    "⩚" => TimesOp,
    "⩮" => ComparisonOp,
    "⤓" => PowerOp,
    "⥠" => PowerOp,
    "⤁" => ArrowOp,
    "≌" => ComparisonOp,
    "⭌" => ArrowOp,
    "⇐" => 0,
    ":" => ColonOp,
    "⇍" => 0,
    "⩭" => ComparisonOp,
    "↛" => ArrowOp,
    "↻" => 0,
    "⪲" => ComparisonOp,
    "⇔" => ArrowOp,
    "⪭" => ComparisonOp,
    "⩱" => ComparisonOp,
    "⤍" => ArrowOp,
    "⊃" => ComparisonOp,
    "∝" => ComparisonOp,
    "⊊" => ComparisonOp,
    "⊠" => TimesOp,
    "∘" => TimesOp,
    "⇚" => 0,
    "⪆" => ComparisonOp,
    "⫋" => ComparisonOp,
    "|" => PlusOp,
    "⊉" => ComparisonOp,
    "∻" => ComparisonOp,
    "≉" => ComparisonOp,
    "⩡" => PlusOp,
    "⨬" => PlusOp,
    "⩷" => ComparisonOp,
    ">=" => ComparisonOp,
    "..." => 0,
    "where" => WhereOp,
    "->" => AnonFuncOp,
    "|=" => AssignmentOp,
    "'" => PrimeOp,
    "!" => 0,
    "√" => 0,
    "∛" => 0,
    "∜" => 0,
    "++" => PlusOp,
    "¬" => 0
)


text:
# Wishart distribution
#
#   following the Wikipedia parameterization
#
"""
    Wishart(ν, S)
```julia
ν::Real           degrees of freedom (greater than p - 1)
S::AbstractPDMat  p x p scale matrix
```
The [Wishart distribution](http://en.wikipedia.org/wiki/Wishart_distribution)
generalizes the gamma distribution to ``p\\times p`` real, positive definite
matrices ``\\mathbf{H}``. If ``\\mathbf{H}\\sim \\textrm{W}_p(\\nu,\\mathbf{S})``,
then its probability density function is

```math
f(\\mathbf{H};\\nu,\\mathbf{S}) = \\frac{1}{2^{\\nu p/2} \\left|\\mathbf{S}\\right|^{\\nu/2} \\Gamma_p\\left(\\frac {\\nu}{2}\\right ) }{\\left|\\mathbf{H}\\right|}^{(\\nu-p-1)/2} e^{-(1/2)\\operatorname{tr}(\\mathbf{S}^{-1}\\mathbf{H})}.
```

If ``\\nu`` is an integer, then a random matrix ``\\mathbf{H}`` given by

```math
\\mathbf{H} = \\mathbf{X}\\mathbf{X}^{\\rm{T}}, \\quad\\mathbf{X} \\sim \\textrm{MN}_{p,\\nu}(\\mathbf{0}, \\mathbf{S}, \\mathbf{I}_{\\nu})
```

has ``\\mathbf{H}\\sim \\textrm{W}_p(\\nu, \\mathbf{S})``. For non-integer degrees of freedom,
Wishart matrices can be generated via the [Bartlett decomposition](https://en.wikipedia.org/wiki/Wishart_distribution#Bartlett_decomposition).
"""
struct Wishart{T<:Real, ST<:AbstractPDMat} <: ContinuousMatrixDistribution
    df::T     # degree of freedom
    S::ST     # the scale matrix
    logc0::T  # the logarithm of normalizing constant in pdf
end

#  -----------------------------------------------------------------------------
#  Constructors
#  -----------------------------------------------------------------------------

function Wishart(df::T, S::AbstractPDMat{T}) where T<:Real
    p = dim(S)
    df > p - 1 || throw(ArgumentError("df should be greater than dim - 1."))
    logc0 = wishart_logc0(df, S)
    R = Base.promote_eltype(T, logc0)
    prom_S = convert(AbstractArray{T}, S)
    Wishart{R, typeof(prom_S)}(R(df), prom_S, R(logc0))
end

function Wishart(df::Real, S::AbstractPDMat)
    T = Base.promote_eltype(df, S)
    Wishart(T(df), convert(AbstractArray{T}, S))
end

Wishart(df::Real, S::Matrix) = Wishart(df, PDMat(S))

Wishart(df::Real, S::Cholesky) = Wishart(df, PDMat(S))

#  -----------------------------------------------------------------------------
#  REPL display
#  -----------------------------------------------------------------------------

show(io::IO, d::Wishart) = show_multline(io, d, [(:df, d.df), (:S, Matrix(d.S))])

#  -----------------------------------------------------------------------------
#  Conversion
#  -----------------------------------------------------------------------------

function convert(::Type{Wishart{T}}, d::Wishart) where T<:Real
    P = convert(AbstractArray{T}, d.S)
    Wishart{T, typeof(P)}(T(d.df), P, T(d.logc0))
end
function convert(::Type{Wishart{T}}, df, S::AbstractPDMat, logc0) where T<:Real
    P = convert(AbstractArray{T}, S)
    Wishart{T, typeof(P)}(T(df), P, T(logc0))
end

#  -----------------------------------------------------------------------------
#  Properties
#  -----------------------------------------------------------------------------

insupport(::Type{Wishart}, X::Matrix) = isposdef(X)
insupport(d::Wishart, X::Matrix) = size(X) == size(d) && isposdef(X)

dim(d::Wishart) = dim(d.S)
size(d::Wishart) = (p = dim(d); (p, p))
rank(d::Wishart) = dim(d)
params(d::Wishart) = (d.df, d.S)
@inline partype(d::Wishart{T}) where {T<:Real} = T

mean(d::Wishart) = d.df * Matrix(d.S)

function mode(d::Wishart)
    r = d.df - dim(d) - 1.0
    r > 0.0 || throw(ArgumentError("mode is only defined when df > p + 1"))
    return Matrix(d.S) * r
end

function meanlogdet(d::Wishart)
    p = dim(d)
    df = d.df
    v = logdet(d.S) + p * logtwo
    for i = 1:p
        v += digamma(0.5 * (df - (i - 1)))
    end
    return v
end

function entropy(d::Wishart)
    p = dim(d)
    df = d.df
    -d.logc0 - 0.5 * (df - p - 1) * meanlogdet(d) + 0.5 * df * p
end

#  Gupta/Nagar (1999) Theorem 3.3.15.i
function cov(d::Wishart, i::Integer, j::Integer, k::Integer, l::Integer)
    S = Matrix(d.S)
    d.df * (S[i, k] * S[j, l] + S[i, l] * S[j, k])
end

function var(d::Wishart, i::Integer, j::Integer)
    S = Matrix(d.S)
    d.df * (S[i, i] * S[j, j] + S[i, j] ^ 2)
end

#  -----------------------------------------------------------------------------
#  Evaluation
#  -----------------------------------------------------------------------------

function wishart_logc0(df::Real, S::AbstractPDMat)
    h_df = df / 2
    p = dim(S)
    -h_df * (logdet(S) + p * typeof(df)(logtwo)) - logmvgamma(p, h_df)
end

function logkernel(d::Wishart, X::AbstractMatrix)
    df = d.df
    p = dim(d)
    Xcf = cholesky(X)
    0.5 * ((df - (p + 1)) * logdet(Xcf) - tr(d.S \ X))
end

#  -----------------------------------------------------------------------------
#  Sampling
#  -----------------------------------------------------------------------------

function _rand!(rng::AbstractRNG, d::Wishart, A::AbstractMatrix)
    _wishart_genA!(rng, dim(d), d.df, A)
    unwhiten!(d.S, A)
    A .= A * A'
end

function _wishart_genA!(rng::AbstractRNG, p::Int, df::Real, A::AbstractMatrix)
    # Generate the matrix A in the Bartlett decomposition
    #
    #   A is a lower triangular matrix, with
    #
    #       A(i, j) ~ sqrt of Chisq(df - i + 1) when i == j
    #               ~ Normal()                  when i > j
    #
    A .= zero(eltype(A))
    for i = 1:p
        @inbounds A[i,i] = rand(rng, Chi(df - i + 1.0))
    end
    for j in 1:p-1, i in j+1:p
        @inbounds A[i,j] = randn(rng)
    end
end

#  -----------------------------------------------------------------------------
#  Test utils
#  -----------------------------------------------------------------------------

function _univariate(d::Wishart)
    check_univariate(d)
    df, S = params(d)
    α = df / 2
    β = 2Matrix(S)[1]
    return Gamma(α, β)
end

function _rand_params(::Type{Wishart}, elty, n::Int, p::Int)
    n == p || throw(ArgumentError("dims must be equal for Wishart"))
    ν = elty( n + 1 + abs(10randn()) )
    S = (X = 2rand(elty, n, n) .- 1; X * X')
    return ν, S
end

ext: jl


text:
using VisJS
using Colors
using LightGraphs

label = ["Node $i" for i in 1:100]
g = Graph(100)
gr = [ "$(i%3)" for i in 1:100]
VisJS.drawGraph(g,label=label,group=gr)


text:
# Metagraphs files are simply JLD2 files.

struct MGFormat <: AbstractGraphFormat end

function loadmg(fn::AbstractString)
    @load fn g
    return g
end

function savemg(fn::AbstractString, g::AbstractMetaGraph)
    @save fn g
    return 1
end

loadgraph(fn::AbstractString, gname::String, ::MGFormat) = loadmg(fn)
savegraph(fn::AbstractString, g::AbstractMetaGraph) =  savemg(fn, g)


text:
export Result,
    ResultStatistics,
    statistics,
    seed,
    path_results,
    results,
    nresults,
    solutions,
    real_solutions,
    nonsingular,
    singular,
    at_infinity,
    failed,
    nfinite,
    nsolutions,
    nsingular,
    nat_infinity,
    nexcess_solutions,
    nfailed,
    nnonsingular,
    nreal,
    ntracked

"""
    Result

The result of [`solve`](@ref). This is a wrapper around the results of each single path
([`PathResult`](@ref)) and it contains some additional information like a random seed to
replicate the result.
"""
struct Result
    path_results::Vector{PathResult}
    tracked_paths::Int
    seed::Union{Nothing,UInt32}
    start_system::Union{Nothing,Symbol}
end
Result(
    path_results;
    seed::Union{Nothing,UInt32} = nothing,
    tracked_paths = length(path_results),
    start_system = nothing,
) = Result(path_results, tracked_paths, seed, start_system)

Base.size(r::Result) = (length(r),)
Base.length(r::Result) = length(r.path_results)
Base.getindex(r::Result, I) = getindex(r.path_results, I)

Base.iterate(r::Result) = iterate(r.path_results)
Base.iterate(r::Result, state) = iterate(r.path_results, state)
Base.lastindex(r::Result) = lastindex(r.path_results)
Base.eltype(r::Type{Result}) = PathResult
Base.keys(r::Result) = 1:length(r.path_results)
Base.values(r::Result) = r.path_results

"""
    seed(::Result)

Returns the seed to replicate the result.
"""
seed(R::Result) = R.seed

"""
    path_results(::Result)

Returns the stored [`PathResult`](@ref)s.
"""
path_results(R::Result) = R.path_results

#TODO: is_multiple_result(r::PathResult, r::Result)
is_multiple_result(path_result::PathResult, result) = false

Base.@kwdef struct ResultStatistics
    total::Int
    nonsingular::Int
    singular::Int
    singular_with_multiplicity::Int = singular
    real::Int
    real_nonsingular::Int
    real_singular::Int
    real_singular_with_multiplicity::Int = real_singular
    at_infinity::Int
    excess_solution::Int = 0
    failed::Int = 0
end
Base.show(io::IO, stats::ResultStatistics) = print_fieldnames(io, stats)

function ResultStatistics(
    result::Result;
    real_tol::Float64 = 1e-6,
    singular_tol::Float64 = 1e10,
)
    failed = at_infinity = excess_solution = 0
    nonsingular = singular = real_nonsingular = real_singular = 0
    singular_with_multiplicity = real_singular_with_multiplicity = 0
    for r in result
        is_multiple_result(r, result) && continue
        if is_failed(r)
            failed += 1
        elseif is_at_infinity(r)
            at_infinity += 1
        elseif is_excess_solution(r)
            excess_solution += 1
        elseif is_singular(r, singular_tol)
            if is_real(r, real_tol)
                real_singular += 1
                real_singular_with_multiplicity += unpack(multiplicity(r), 1)
            end
            singular += 1
            singular_with_multiplicity += unpack(multiplicity(r), 1)
        else # finite, nonsingular
            if is_real(r, real_tol)
                real_nonsingular += 1
            end
            nonsingular += 1
        end
    end
    ResultStatistics(
        nonsingular = nonsingular,
        singular = singular,
        singular_with_multiplicity = singular_with_multiplicity,
        real_nonsingular = real_nonsingular,
        real_singular = real_singular,
        real_singular_with_multiplicity = real_singular_with_multiplicity,
        real = real_nonsingular + real_singular,
        at_infinity = at_infinity,
        excess_solution = excess_solution,
        failed = failed,
        total = result.tracked_paths,
    )
end

"""
    statistics(R::Result; real_tol = 1e-6, singular_tol = 1e10)

Statistic about the number of (real) singular and non-singular solutions etc.
"""
statistics(r; kwargs...) = ResultStatistics(r; kwargs...)

const Results = Union{Result,AbstractVector{<:PathResult}}

"""
    results(
        result;
        only_real = false,
        real_tol = 1e-6,
        only_nonsingular = false,
        only_singular = false,
        singular_tol = 1e10,
        only_finite = true,
        multiple_results = false,
    )
    results(f, result; options...)

Return all [`PathResult`](@ref)s for which satisfy the given conditions and apply if provided
the function `f`.
"""
results(R::Results; kwargs...) = results(identity, R; kwargs...)
function results(
    f::Function,
    R::Results;
    only_real::Bool = false,
    real_tol::Float64 = 1e-6,
    only_nonsingular::Bool = false,
    only_singular::Bool = false,
    singular_tol::Float64 = 1e10,
    onlyfinite::Bool = true, # deprecated
    only_finite::Bool = onlyfinite,
    multiple_results::Bool = false,
)
    [
        f(r)
        for
        r in R if
        (!only_real || is_real(r, real_tol)) &&
            (!only_nonsingular || is_nonsingular(r, singular_tol)) &&
            (!only_singular || is_singular(r, singular_tol)) &&
            (!only_finite || is_finite(r)) &&
            (multiple_results || !is_multiple_result(r, R))
    ]
end

"""
    nresults(
        result;
        only_real = false,
        real_tol = 1e-6,
        only_nonsingular = false,
        only_singular = false,
        singular_tol = 1e10,
        only_finite = true,
        multiple_results = false,
    )

Count the number of solutions which satisfy the corresponding conditions. See also
[`results`](@ref).
"""
function nresults(
    R::Results;
    only_real::Bool = false,
    real_tol::Float64 = 1e-6,
    only_nonsingular::Bool = false,
    only_singular::Bool = false,
    singular_tol::Float64 = 1e10,
    onlyfinit::Bool = true, # deprecated
    only_finite::Bool = onlyfinit,
    multiple_results::Bool = false,
)
    count(R) do r
        (!only_real || is_real(r, real_tol)) &&
            (!only_nonsingular || is_nonsingular(r, singular_tol)) &&
            (!only_singular || is_singular(r, singular_tol)) &&
            (!only_finite || isfinite(r)) &&
            (multiple_results || !is_multiple_result(r, R))
    end
end

"""
    solutions(result; conditions...)

Returns all solutions for which the given conditions apply, see [`results`](@ref) for the
possible conditions.

## Example
```julia-repl
julia> @var x y
julia> F = System([(x-2)y, y+x+3]);
julia> solutions(solve(F))
2-element Array{Array{Complex{Float64},1},1}:
 [2.0 + 0.0im, -5.0 + 0.0im]
 [-3.0 + 0.0im, 0.0 + 0.0im]
```
"""
solutions(result::Results; kwargs...) = results(solution, result; kwargs...)

"""
    real_solutions(result; tol=1e-6, conditions...)

Return all real solution for which the given conditions apply.
For the possible `conditions` see [`results`](@ref).
Note that `only_real` is always `true` and `real_tol` is now `tol`.

## Example
```julia-repl
julia> @var x y;
julia> F = System([(x-2)y, y+x+3]);
julia> real_solutions(solve(F))
2-element Array{Array{Float64,1},1}:
 [2.0, -5.0]
 [-3.0, 0.0]
```
"""
function real_solutions(result::Results; tol::Float64 = 1e-6, kwargs...)
    results(real ∘ solution, result; only_real = true, real_tol = tol, kwargs...)
end

"""
    nonsingular(result; conditions...)

Return all [`PathResult`](@ref)s for which the solution is non-singular.
This is just a shorthand for `results(R; only_nonsingular=true, conditions...)`.
For the possible `conditions` see [`results`](@ref).
"""
nonsingular(R::Results; kwargs...) = results(R; only_nonsingular = true, kwargs...)

"""
    singular(result; tol=1e10, multiple_results=false, kwargs...)

Return all [`PathResult`]s for which the solution is singular.
A solution is labeled singular
if the condition number is greater than `singular_tol`, or if the winding number is > 1.
If `multiple_results=false` only one point from each cluster of multiple solutions is returned.
If If `multiple_results=true` all singular solutions in `R` are returned.
For the possible `kwargs` see [`results`](@ref).
"""
function singular(R::Results; tol::Float64 = 1e10, kwargs...)
    results(R; only_singular = true, singular_tol = tol, kwargs...)
end

"""
    real(result, tol=1e-6)

Get all results where the solutions are real with the given tolerance `tol`.
See [`is_real`](@ref) for details regarding the determination of 'realness'.
"""
Base.real(R::Results; tol::Float64 = 1e-6) = filter(r -> is_real(r, tol), path_results(R))

"""
    failed(result)
Get all results where the path tracking failed.
"""
failed(R::Results) = filter(is_failed, path_results(R))

"""
    at_infinity(result)

Get all results where the solutions is at infinity.
"""
at_infinity(R::Results) = filter(is_at_infinity, path_results(R))

# """
#     multiplicities(result; tol =1 e-6)
#
# Returns a `Vector` of `Vector{PathResult}`s grouping the `PathResult`s whose solutions
# appear with multiplicities *greater* than 1 in 'result'.
# Two solutions are regarded as equal, when their pairwise distance is less than 'tol'.
# """
# function multiplicities(results::Results; tol::Float64 = 1e-6)
#     map(i -> results[i], multiplicities(solution, results; tol = tol))
# end

"""
    nsolutions(result)

The number of solutions.
"""
nsolutions(R::Results) = nresults(R)

"""
    nsingular(
        result;
        singular_tol = 1e10,
        counting_multiplicities = false,
        kwargs...,
    )

The number of singular solutions. A solution is considered singular if its winding number is
larger than 1 or the condition number is larger than `tol`.
If `counting_multiplicities=true` the number of singular solutions times their
multiplicities is returned.
"""
function nsingular(
    R::Results;
    singular_tol::Float64 = 1e10,
    counting_multiplicities::Bool = false,
    kwargs...,
)
    S = results(
        R;
        only_singular = true,
        multiple_results = false,
        singular_tol = singular_tol,
        kwargs...,
    )
    isempty(S) && return 0
    counting_multiplicities && return sum(multiplicity, S)
    length(S)
end

"""
    nat_infinity(result)

The number of solutions at infinity.
"""
nat_infinity(R::Results) = count(is_at_infinity, R)

"""
    nexcess_solutions(result)

The number of exess solutions. See also [`excess_solution_check`](@ref).
"""
nexcess_solutions(R::Results) = count(is_excess_solution, R)

"""
    nfailed(result)

The number of failed paths.
"""
nfailed(R::Results) = count(is_failed, R)

"""
    nnonsingular(result; tol=1e-10)

The number of non-singular solutions. See also [`is_singular`](@ref).
"""
nnonsingular(R::Result; tol = 1e10) = count(r -> is_nonsingular(r, tol), R)

"""
    nreal(result; tol=1e-6)

The number of real solutions. See also [`is_real`](@ref).
"""
nreal(R::Results; tol = 1e-6) = count(r -> is_real(r, tol), R)

"""
    ntracked(result)

Returns the total number of paths tracked.
"""
ntracked(R::Result) = R.tracked_paths

###
### Show
####
function Base.show(io::IO, x::Result)
    s = statistics(x)
    total = s.nonsingular + s.singular
    header = "Result with $total $(plural("solution", total))"
    println(io, header)
    println(io, "="^(length(header)))
    println(io, "• $(ntracked(x)) $(plural("path", ntracked(x))) tracked")
    println(
        io,
        "• $(s.nonsingular) non-singular $(plural("solution", s.nonsingular)) ",
        "($(s.real_nonsingular) real)",
    )
    if s.singular > 0
        println(
            io,
            "• $(s.singular) singular $(plural("solution", s.singular)) ",
            "($(s.real_singular) real)",
        )
    end
    s.at_infinity > 0 &&
        println(io, "• $(s.at_infinity) $(plural("solution", s.at_infinity)) at infinity")
    s.excess_solution > 0 && println(
        io,
        "• $(s.excess_solution) excess $(plural("solution", s.excess_solution))",
    )
    println(io, "• random seed: ", sprint(show, seed(x)))
    if !isnothing(x.start_system)
        println(io, "• start_system: :", x.start_system)
    end
    if s.singular > 0
        println(io, "• multiplicity table of singular solutions:")
        #     singular_multiplicities_table(io, x, s)
    end
end

function TreeViews.treelabel(io::IO, x::Result, ::MIME"application/prs.juno.inline")
    s = statistics(x)
    print(
        io,
        "Result with $(s.nonsingular + s.singular) solutions",
    )
end
TreeViews.hastreeview(::Result) = true
TreeViews.numberofnodes(::Result) = 9
function TreeViews.nodelabel(io::IO, x::Result, i::Int, ::MIME"application/prs.juno.inline")
    s = statistics(x)
    if i == 1
        print(io, "Paths tracked")
    elseif i == 2 && s.nonsingular > 0
        print(io, "$(s.nonsingular) non-singular ($(s.real_nonsingular) real)")
    elseif i == 3 && s.singular > 0
        print(io, "$(s.singular) singular ($(s.real_singular) real)")
    elseif i == 4 && (s.real_nonsingular + s.real_singular) > 0
        print(io, "$(s.real_nonsingular+s.real_singular) real")
    elseif i == 5 && s.at_infinity > 0
        print(io, "$(s.at_infinity) at_infinity")
    elseif i == 6 && s.failed > 0
        print(io, "$(s.failed) failed")
    elseif i == 7
        print(io, "Random seed used")
    elseif i == 8
        print(io, "start_system")
    elseif i == 9 && s.singular > 0
        print(io, "  multiplicity table of singular solutions: \n")
        # singular_multiplicities_table(io, x, s)
    end
end

function TreeViews.treenode(r::Result, i::Integer)
    s = statistics(r)
    if i == 1
        return ntracked(r)
    elseif i == 2 && s.nonsingular > 0
        return results(r, only_nonsingular = true)
    elseif i == 3 && s.singular > 0
        return results(r, only_singular = true)
    elseif i == 4 && (s.real_nonsingular + s.real_singular) > 0
        return results(r, only_real = true)
    elseif i == 5 && s.at_infinity > 0
        return at_infinity(r)
    elseif i == 6 && s.failed > 0
        return failed(r)
    elseif i == 7
        return seed(r)
    elseif i == 8
        return something(r.start_system, missing)
    elseif i == 9 && s.singular > 0
        return missing
    end
    missing
end

# function singular_multiplicities_table(io, result::Result, stats = statistics(result))
#     M = result.multiplicity_info[].multiplicities
#     if isempty(M)
#         n_higher_mults_total = 0
#     else
#         n_higher_mults_total = sum(length, values(M))
#     end
#
#     headers = ["mult.", "total", "# real", "# non-real"]
#     mult_1_exists = n_higher_mults_total < stats.singular
#     data = Matrix{Int}(undef, mult_1_exists + length(M), 4)
#
#     n_higher_real = 0
#     i = mult_1_exists + 1
#     for k in sort(collect(keys(M)))
#         data[i, 1] = k
#         n_real_solsᵢ = count(Mᵢ -> is_real(result.pathresults[Mᵢ[1]]), M[k])
#         n_solsᵢ = length(M[k])
#         data[i, 2] = n_solsᵢ
#         data[i, 3] = n_real_solsᵢ
#         data[i, 4] = n_solsᵢ - n_real_solsᵢ
#         n_higher_real += n_real_solsᵢ
#         i += 1
#     end
#
#     if mult_1_exists
#         data[1, 1] = 1
#         data[1, 2] = stats.singular - n_higher_mults_total
#         data[1, 3] = stats.real_singular - n_higher_real
#         data[1, 4] = data[1, 2] - data[1, 3]
#     end
#
#     PrettyTables.pretty_table(
#         io,
#         data,
#         headers;
#         tf = PrettyTables.unicode,
#         alignment = :c,
#         header_crayon = PrettyTables.Crayon(bold = false),
#         border_crayon = PrettyTables.Crayon(faint = true),
#     )
# end


text:
using FewShotAnomalyDetection
using Flux
using MLDataPattern
using FluxExtensions
using Adapt
using DataFrames
using CSV
using Serialization
using IPMeasures
using Random
using Statistics

include("experimentalutils.jl")

outputFolder = mainfolder * "experiments/twostage_vae_corrjac/"
mkpath(outputFolder)

function null_distr_distances(dim, size, k = 500)
    z1 = zeros(Float32, dim, k)
    z2 = copy(z1)
    randn!(z1)
    randn!(z2)
    γs = -10:0.05:2
    cs = [IPMeasures.crit_mmd2_var(IPMeasures.IMQKernel(10.0 ^ γ), z1, z2, IPMeasures.pairwisel2) for γ in γs]
    γ = 10 ^ γs[argmax(cs)]
    null_dst = zeros(Float32, k)
    for i in 1:k
        randn!(z1)
        randn!(z2)
        null_dst[i] = IPMeasures.mmd(IPMeasures.IMQKernel(γ), z1, z2)
    end
    sort!(null_dst)
    return null_dst, γ
end

function get_gamma(m::VAE, x)
    z = Flux.Tracker.data(samplez(m, x))
    zp = randn(size(z))
    γs = -10:0.05:2
    cs = [IPMeasures.crit_mmd2_var(IPMeasures.IMQKernel(10.0 ^ γ), z, zp, 1000, IPMeasures.pairwisel2) for γ in γs]
    γ = 10 ^ γs[argmax(cs)]
end

mmdpval(null_dst, x) = searchsortedfirst(null_dst, x) / length(null_dst)

function runExperiment(datasetName, train, test, inputDim, hiddenDim, hiddenDim2, latentDim, numLayers, nonlinearity, layerType, batchSize = 100, numBatches = 10000, γ_step = 1000, i = 0)

    T = Float32

    println("$datasetName: computing null distribution distances...")
    null_dst, null_γ = null_distr_distances(latentDim, size(train[1], 2))

    println("$datasetName: creating networks...")
    outer_encoder = Adapt.adapt(T, FluxExtensions.layerbuilder(inputDim, hiddenDim, latentDim * 2, numLayers + 1, nonlinearity, "linear", layerType))
    outer_decoder = Adapt.adapt(T, FluxExtensions.layerbuilder(latentDim, hiddenDim, inputDim + 1, numLayers + 1, nonlinearity, "linear", layerType))
    outerVAE = VAE(outer_encoder, outer_decoder, T(1), :scalarsigma)

    inner_encoder = Adapt.adapt(T, FluxExtensions.layerbuilder(latentDim, hiddenDim2, latentDim * 2, numLayers + 1, nonlinearity, "linear", layerType))
    inner_decoder = Adapt.adapt(T, FluxExtensions.layerbuilder(latentDim, hiddenDim2, latentDim + 1, numLayers + 1, nonlinearity, "linear", layerType))
    innerVAE = VAE(inner_encoder, inner_decoder, T(1), :scalarsigma)

    println("$datasetName: training outer VAE...")
    opt = Flux.Optimise.ADAM(1e-4)
    for i in 1:(numBatches ÷ γ_step)
        γ = get_gamma(outerVAE, train[1])
        dist_fun = (x, y) -> IPMeasures.mmd(IPMeasures.IMQKernel(γ), x, y)
        cb = Flux.throttle(() -> println("$datasetName outer VAE: $(printing_wloss(outerVAE, train[1], dist_fun))"), 5)
        Flux.train!(x -> wloss(outerVAE, x, dist_fun), Flux.params(outerVAE), RandomBatches((train[1],), batchSize, γ_step), opt, cb = cb)
    end
    z_train = Flux.Tracker.data(zparams(outerVAE, train[1])[1])
    z_test = Flux.Tracker.data(zparams(outerVAE, test[1])[1])

    println("$datasetName: training inner VAE...")
    opt = Flux.Optimise.ADAM(1e-4)
    for i in 1:(numBatches ÷ γ_step)
        γ = get_gamma(innerVAE, z_train)
        dist_fun = (x, y) -> IPMeasures.mmd(IPMeasures.IMQKernel(γ), x, y)
        cb = Flux.throttle(() -> println("$datasetName inner VAE: $(printing_wloss(innerVAE, z_train, dist_fun))"), 5)
        Flux.train!(x -> wloss(innerVAE, x, dist_fun), Flux.params(innerVAE), RandomBatches((z_train,), batchSize, γ_step), opt, cb = cb)
    end

    # x -> z -> u -> zp -> xp
    println("$datasetName: computing performance metrics...")
    u_train = Flux.Tracker.data(zparams(innerVAE, z_train)[1])
    u_test = Flux.Tracker.data(zparams(innerVAE, z_test)[1])
    u_train_dst = IPMeasures.mmd(IPMeasures.IMQKernel(null_γ), randn(size(u_train)), u_train)
    u_test_dst = IPMeasures.mmd(IPMeasures.IMQKernel(null_γ), randn(size(u_test)), u_test)
    u_train_mmd_pval = mmdpval(null_dst, u_train_dst)
    u_test_mmd_pval = mmdpval(null_dst, u_test_dst)
    zp_train = Flux.Tracker.data(hsplit1softp(innerVAE.g(u_train))[1])
    zp_test = Flux.Tracker.data(hsplit1softp(innerVAE.g(u_test))[1])
    xp_train = Flux.Tracker.data(hsplit1softp(outerVAE.g(zp_train))[1])
    xp_test = Flux.Tracker.data(hsplit1softp(outerVAE.g(zp_test))[1])

    log_pxv_train = vec(collect(log_pxexpectedz(outerVAE, train[1], zp_train)'))
    log_pxv_test = vec(collect(log_pxexpectedz(outerVAE, test[1], zp_test)'))
    log_pzs = vec(collect(log_pz(innerVAE, z_test)'))
    log_det_jac_enc_out = vec(collect(log_det_jacobian_encoder(outerVAE, test[1])'))
    log_det_jac_dec_out = vec(collect(log_det_jacobian_decoder(outerVAE, zp_test)'))
    log_det_jac_enc_inn = vec(collect(log_det_jacobian_encoder(innerVAE, z_test)'))
    log_det_jac_dec_inn = vec(collect(log_det_jacobian_decoder(innerVAE, u_test)'))

    labels = test[2] .- 1
    auc_pxv = computeauc(.-log_pxv_test, labels)
    auc_pz = computeauc(.-log_pzs, labels)
    auc_pz_jacobian_enc = computeauc(.-(log_pzs .+ log_det_jac_enc_out .+ log_det_jac_enc_inn), labels)
    auc_pz_jacobian_dec = computeauc(.-(log_pzs .- log_det_jac_dec_out .- log_det_jac_dec_inn), labels)
    auc_pxv_pz = computeauc(.-(log_pxv_test .+ log_pzs), labels)
    auc_pxv_pz_jacobian_enc = computeauc(.-(log_pxv_test .+ log_pzs .+ log_det_jac_enc_out .+ log_det_jac_enc_inn), labels)
    auc_pxv_pz_jacobian_dec = computeauc(.-(log_pxv_test .+ log_pzs .- log_det_jac_dec_out .- log_det_jac_dec_inn), labels)

    serialize(outputFolder * "twostage_vae-$datasetName-$i-$hiddenDim-$hiddenDim2-$latentDim-innvae.jls", innerVAE)
    serialize(outputFolder * "twostage_vae-$datasetName-$i-$hiddenDim-$hiddenDim2-$latentDim-outvae.jls", outerVAE)
    serialize(outputFolder * "twostage_vae-$datasetName-$i-$hiddenDim-$hiddenDim2-$latentDim-train.jls", train)
    serialize(outputFolder * "twostage_vae-$datasetName-$i-$hiddenDim-$hiddenDim2-$latentDim-test.jls", test)

    df = DataFrame(dataset = datasetName, idim = inputDim, hdim = hiddenDim, hdim2 = hiddenDim2, ldim = latentDim, layers = numLayers, i = i, auc_pxv = auc_pxv, auc_pz = auc_pz, auc_pz_jacobian_enc = auc_pz_jacobian_enc,
                    auc_pz_jacobian_dec = auc_pz_jacobian_dec, auc_pxv_pz = auc_pxv_pz, auc_pxv_pz_jacobian_enc = auc_pxv_pz_jacobian_enc, auc_pxv_pz_jacobian_dec = auc_pxv_pz_jacobian_dec,
                    rec_err_train = .-mean(log_pxv_train), rec_err_test = .-mean(log_pxv_test), u_train_mmd_pval = u_train_mmd_pval, u_test_mmd_pval = u_test_mmd_pval, u_train_dst = u_train_dst,
                    u_test_dst = u_test_dst, meanl2_train = mean(IPMeasures.pairwisel2(train[1], xp_train)), meanl2_test = mean(IPMeasures.pairwisel2(test[1], xp_test)),
                    medianl2_train = median(IPMeasures.pairwisel2(train[1], xp_train)), medianl2_test = median(IPMeasures.pairwisel2(test[1], xp_test)))

    CSV.write(outputFolder * "twostage_vae-$datasetName-$i-$hiddenDim-$hiddenDim2-$latentDim-results.csv", df)
    return df
end

datasets = ["abalone"]
difficulties = ["easy"]
batchSize = 100
iterations = 10000
γ_step = 1000

if length(ARGS) != 0
    datasets = [ARGS[1]]
    difficulties = ["easy"]
end

# datasets = ["abalone"]
# difficulties = ["easy"]

for i in 1:5
    for (dn, df) in zip(datasets, difficulties)
        train, test, clusterdness = loaddata(dn, df)

        evaluateOneConfig = p -> runExperiment(dn, train, test, size(train[1], 1), p..., batchSize, iterations, γ_step, i)
        results = gridsearch(evaluateOneConfig, [64, 32], [32, 16], [8, 2], [3], ["swish"], ["Dense"])

        CSV.write(outputFolder * "twostage_vae-$dn-$i.csv", vcat(results...))
    end
end


text:
"""
    markers(x)

Setup different markers for use in `NumberLinePlot`.
"""
function markers(x)
    m=fill(:rtriangle,size(x))
    m[x.<0] .= :ltriangle
    m[x.==0] .= :dot
    o=fill(-0.5,size(x))
    o[x.<0] .= 0.5
    o[x.==0] .= 0.0
    c=fill(:blue,size(x))
    c[x.<0] .= :red
    c[x.==0] .= :green
    m,o,c
end

"""
    arrows(x)

Setup arrows pointing left (in red) or right (in blue) for use in `NumberLinePlot`.
"""
function arrows(x)
    length(x)>0 ? s=size(x) : s=(1,1)
    a=fill("⊚",s)
    p=fill(:center,s)
    c=fill(:green,s)
    if length(x)>1
        a[x.<0] .= "<"
        p[x.<0] .= :left
        c[x.<0] .= :red
        a[x.>0] .= ">"
        p[x.>0] .= :right
        c[x.>0] .= :blue
    end
    a,p,c
end

"""
    NumberLinePlot(x)

Display the number line sequence of operations, defined by vector `x`, as a graph.
"""
function NumberLinePlot(y)
	x=y[findall((!isnan).(y))]

    nx=length(x)
    nx>1 ? xx=cumsum([0 ;x[2:end]]) : xx=[0.;0.]
    nx>1 ? yy=collect(0:length(x)-1) : yy=[0.;0.]

    x∞=max(maximum(abs.(extrema(xx))),10)
    y∞=maximum(yy)

    xfac=min(x∞,10)/x∞
    yfac=min(y∞+1,10)/(y∞+1)
    f(p::Point) = Point(xfac*200.0*p[1],-yfac*100.0*p[2]+800.0)
	g=10.0
    xdel=Int(ceil(x∞/10))

	@svg begin
		setline(0.5g)
		line(f(Point(-x∞,0)), f(Point(x∞,0)), :stroke)
		[line(f(Point(i,+0.5)), f(Point(i,-0.5)), :stroke) for i in -x∞:x∞]

        [settext(" -$i ",f(Point(-i,-2.0)),
        halign="center",valign="bottom",markup=true,angle=0) for i in xdel:xdel:Int(x∞)]
        settext(" 0 ",f(Point(0,-2.0)),
        halign="center",valign="bottom",markup=true,angle=0)
        [settext(" +$i ",f(Point(i,-2.0)),
        halign="center",valign="bottom",markup=true,angle=0) for i in xdel:xdel:Int(x∞)]

		setline(1.5g)
		[line(f(Point(xx[i],yy[i])), f(Point(xx[i],yy[i+1])),:stroke) for i in 1:nx-1]
  		[line(f(Point(xx[i],yy[i+1])), f(Point(xx[i+1],yy[i+1])),:stroke) for i in 1:nx-1]

		setline(1.5g)
	    setcolor(0,0.75,0)
		#line(f(Point(xx[end],yy[end])), f(Point(xx[end],yy[1])),:stroke)
        a=f(Point(xx[end],yy[end]))
        b=f(Point(xx[end],yy[1]))
        yy[end]!==yy[1] ? arrow(a,b,arrowheadlength=5.0g, arrowheadangle=pi/4, linewidth=1.5g) : nothing

    end 5000 2000 joinpath(tempdir(),"tmp.svg")
end

"""
    NumberLineExpression(x)

Display the number line sequence of operations, defined by vector `x`, as text.
"""
function NumberLineExpression(x)
	tmp1=deepcopy(x)
	xpr=""
	while length(tmp1)>1&&!isnan(tmp1[2])
		tmp2=NumberLineType(pop!(tmp1))
		tmp2<0 ? xpr="$(tmp2)"*xpr : xpr="+$(tmp2)"*xpr
	end
	"Expression : 0"*xpr
end

##

import Base.show

"""
   aSlider(range::AbstractRange,default::Number,show_value::Bool)

Data type used for sliders in Pluto notebooks
"""
struct aSlider
    range::AbstractRange
    default::Number
    show_value::Bool
end

aSlider(range::AbstractRange; default=missing, show_value=true) = aSlider(range, (default === missing) ? first(range) : default, show_value)

function show(io::IO, ::MIME"text/html", slider::aSlider)
    print(io, """<input
        type="range"
        min="$(first(slider.range))"
        step="$(step(slider.range))"
        max="$(last(slider.range))"
        value="$(slider.default)"
        $(slider.show_value ? "oninput=\"this.nextElementSibling.value=this.value\"" : "")
        >""")

    print(io, "  [ x = </output>")

    if slider.show_value
        print(io, """<output>$(slider.default)</output>""")
    end

    print(io, " ] </output>")

end


text:
include("ParserTests.jl")
include("TransformTests.jl")
include("ActionTests.jl")


text:
abstract Simplexer

immutable AffineSimplexer <: Simplexer
    a::Float64
    b::Float64
end
AffineSimplexer(;a = 0.025, b = 0.5) = AffineSimplexer(a, b)

function simplexer{T, N}(S::AffineSimplexer, initial_x::Array{T, N})
    n = length(initial_x)
    initial_simplex = Array{T, N}[copy(initial_x) for i = 1:n+1]
    for j = 1:n
        initial_simplex[j+1][j] = (1+S.b) * initial_simplex[j+1][j] + S.a
    end
    initial_simplex
end

abstract NMParameters

immutable AdaptiveParameters <: NMParameters
    α::Float64
    β::Float64
    γ::Float64
    δ::Float64
end

AdaptiveParameters(;  α = 1.0, β = 1.0, γ = 0.75 , δ = 1.0) = AdaptiveParameters(α, β, γ, δ)
parameters(P::AdaptiveParameters, n::Integer) = (P.α, P.β + 2/n, P.γ - 1/2n, P.δ - 1/n)

immutable FixedParameters <: NMParameters
    α::Float64
    β::Float64
    γ::Float64
    δ::Float64
end

FixedParameters(; α = 1.0, β = 2.0, γ = 0.5, δ = 0.5) = FixedParameters(α, β, γ, δ)
parameters(P::FixedParameters, n::Integer) = (P.α, P.β, P.γ, P.δ)

immutable NelderMead{Ts <: Simplexer, Tp <: NMParameters} <: Optimizer
    initial_simplex::Ts
    parameters::Tp
end

function NelderMead(; kwargs...)
    KW = Dict(kwargs)
    if haskey(KW, :a) || haskey(KW, :g) || haskey(KW, :b)
        a, g, b = 1.0, 2.0, 0.5
        haskey(KW, :a) && (a = KW[:a])
        haskey(KW, :g) && (g = KW[:g])
        haskey(KW, :b) && (b = KW[:b])
        return NelderMead(a, g, b)
    elseif haskey(KW, :initial_simplex) || haskey(KW, :parameters)
        initial_simplex, parameters = AffineSimplexer(), AdaptiveParameters()
        haskey(KW, :initial_simplex) && (initial_simplex = KW[:initial_simplex])
        haskey(KW, :parameters) && (parameters = KW[:parameters])
        return NelderMead(initial_simplex, parameters)
    else
        return NelderMead(AffineSimplexer(), AdaptiveParameters())
    end
end

# centroid except h-th vertex
function centroid!{T}(c::Array{T}, simplex, h=0)
    n = length(c)
    fill!(c, zero(T))
    @inbounds for i in 1:n+1
        if i != h
            xi = simplex[i]
            @simd for j in 1:n
                c[j] += xi[j]
            end
        end
    end
    scale!(c, 1/n)
end

centroid(simplex, h) = centroid!(similar(simplex[1]), simplex, h)

nmobjective(y::Vector, m::Integer, n::Integer) = sqrt(var(y) * (m / n))

function print_header(method::NelderMead)
    @printf "Iter     Function value    √(Σ(yᵢ-ȳ)²)/n \n"
    @printf "------   --------------    --------------\n"
end

function Base.show(io::IO, trace::OptimizationTrace{NelderMead})
    @printf io "Iter     Function value    √(Σ(yᵢ-ȳ)²)/n \n"
    @printf io "------   --------------    --------------\n"
    for state in trace.states
        show(io, state)
    end
    return
end

function Base.show(io::IO, t::OptimizationState{NelderMead})
    @printf io "%6d   %14e    %14e\n" t.iteration t.value t.g_norm
    if !isempty(t.metadata)
        for (key, value) in t.metadata
            @printf io " * %s: %s\n" key value
        end
    end
    return
end

type NelderMeadState{T, N}
    @add_generic_fields()
    m::Int64
    simplex::Vector{Array{T,N}}
    x_centroid::Array{T}
    x_lowest::Array{T}
    x_second_highest::Array{T}
    x_highest::Array{T}
    x_reflect::Array{T}
    x_cache::Array{T}
    f_simplex::Array{T}
    nm_x::T
    f_lowest::T
    i_order::Vector{Int64}
    α::T
    β::T
    γ::T
    δ::T
    step_type::String
end

initial_state(method::NelderMead, options, d, initial_x::Array) = initial_state(method, options, d.f, initial_x)
function initial_state{F<:Function, T}(method::NelderMead, options, f::F, initial_x::Array{T})
    n = length(initial_x)
    m = n + 1
    simplex = simplexer(method.initial_simplex, initial_x)
    f_simplex = zeros(T, m)
    @inbounds for i in 1:length(simplex)
        f_simplex[i] = f(simplex[i])
    end
    # Get the indeces that correspond to the ordering of the f values
    # at the vertices. i_order[1] is the index in the simplex of the vertex
    # with the lowest function value, and i_order[end] is the index in the
    # simplex of the vertex with the highest function value
    i_order = sortperm(f_simplex)

    α, β, γ, δ = parameters(method.parameters, n)

NelderMeadState("Nelder-Mead",
          n, # Dimensionality of the problem
          Array{T}(n), # Variable to hold final minimizer value for MultivariateOptimizationResults
          f_simplex[i_order[1]], # Store Nelder Mead objective in state.f_x = state.f_lowest
          m,
          0,
          0,
          m, # Number of vertices in the simplex
          simplex, # Maintain simplex in state.simplex
          centroid(simplex,  i_order[m]), # Maintain centroid in state.centroid
          Array{T}(n), # Store cache in state.x_lowest
          Array{T}(n), # Store cache in state.x_second_highest
          Array{T}(n), # Store cache in state.x_highest
          Array{T}(n), # Store cache in state.x_reflect
          Array{T}(n), # Store cache in state.x_cache
          f_simplex, # Store objective values at the vertices in state.f_simplex
          T(nmobjective(f_simplex, n, m)), # Store nmobjective in state.nm_x
          f_simplex[i_order[1]], # Store lowest f in state.f_lowest
          i_order, # Store a vector of rankings of objective values
          T(α),
          T(β),
          T(γ),
          T(δ),
          "initial")
end

update_state!(d, state::NelderMeadState, method::NelderMead) = update_state!(d.f, state, method)
function update_state!{F<:Function, T}(f::F, state::NelderMeadState{T}, method::NelderMead)
    # Augment the iteration counter
    shrink = false
    n, m = state.n, state.m

    centroid!(state.x_centroid, state.simplex, state.i_order[m])
    copy!(state.x_lowest, state.simplex[state.i_order[1]])
    copy!(state.x_second_highest, state.simplex[state.i_order[n]])
    copy!(state.x_highest, state.simplex[state.i_order[m]])
    state.f_lowest = state.f_simplex[state.i_order[1]]
    f_second_highest = state.f_simplex[state.i_order[n]]
    f_highest = state.f_simplex[state.i_order[m]]

    # Compute a reflection
    @inbounds for j in 1:n
        state.x_reflect[j] = state.x_centroid[j] + state.α * (state.x_centroid[j]-state.x_highest[j])
    end

    f_reflect = f(state.x_reflect)
    state.f_calls += 1
    if f_reflect < state.f_lowest
        # Compute an expansion
        @inbounds for j in 1:n
            state.x_cache[j] = state.x_centroid[j] + state.β *(state.x_reflect[j] - state.x_centroid[j])
        end
        f_expand = f(state.x_cache)
        state.f_calls += 1

        if f_expand < f_reflect
            copy!(state.simplex[state.i_order[m]], state.x_cache)
            @inbounds state.f_simplex[state.i_order[m]] = f_expand
            state.step_type = "expansion"
        else
            copy!(state.simplex[state.i_order[m]], state.x_reflect)
            @inbounds state.f_simplex[state.i_order[m]] = f_reflect
            state.step_type = "reflection"
        end
        # shift all order indeces, and wrap the last one around to the first
        i_highest = state.i_order[m]
        @inbounds for i = m:-1:2
            state.i_order[i] = state.i_order[i-1]
        end
        state.i_order[1] = i_highest
    elseif f_reflect < f_second_highest
        copy!(state.simplex[state.i_order[m]], state.x_reflect)
        @inbounds state.f_simplex[state.i_order[m]] = f_reflect
        state.step_type = "reflection"
        sortperm!(state.i_order, state.f_simplex)
    else
        if f_reflect < f_highest
            # Outside contraction
            @simd for j in 1:n
                @inbounds state.x_cache[j] = state.x_centroid[j] + state.γ * (state.x_reflect[j]-state.x_centroid[j])
            end
            f_outside_contraction = f(state.x_cache)
            if f_outside_contraction < f_reflect
                copy!(state.simplex[state.i_order[m]], state.x_cache)
                @inbounds state.f_simplex[state.i_order[m]] = f_outside_contraction
                state.step_type = "outside contraction"
                sortperm!(state.i_order, state.f_simplex)

            else
                shrink = true
            end
        else # f_reflect > f_highest
            # Inside constraction
            @simd for j in 1:n
                @inbounds state.x_cache[j] = state.x_centroid[j] - γ *(state.x_reflect[j] - state.x_centroid[j])
            end
            f_inside_contraction = f(state.x_cache)
            if f_inside_contraction < f_highest
                copy!(state.simplex[state.i_order[m]], state.x_cache)
                @inbounds state.f_simplex[state.i_order[m]] = f_inside_contraction
                state.step_type = "inside contraction"
                sortperm!(state.i_order, state.f_simplex)
            else
                shrink = true
            end
        end
    end

    if shrink
        for i = 2:m
            ord = state.i_order[i]
            copy!(state.simplex[ord], state.x_lowest + state.δ*(state.simplex[ord]-state.x_lowest))
            state.f_simplex[ord] = f(state.simplex[ord])
        end
        step_type = "shrink"
        sortperm!(state.i_order, state.f_simplex)
    end

    state.nm_x = nmobjective(state.f_simplex, n, m)
    false
end

after_while!(d, state, method::NelderMead, options) = after_while!(d.f, state, method::NelderMead, options)
function after_while!{F<:Function}(f::F, state, method::NelderMead, options)
    sortperm!(state.i_order, state.f_simplex)
    x_centroid_min = centroid(state.simplex, state.i_order[state.m])
    f_centroid_min = f(x_centroid_min)
    state.f_calls += 1
    f_min, i_f_min = findmin(state.f_simplex)
    x_min = state.simplex[i_f_min]
    if f_centroid_min < f_min
        x_min = x_centroid_min
        f_min = f_centroid_min
    end
    state.f_x = f_min
    state.x[:] = x_min
end


text:
using Test
using TestSetExtensions
using MetaBalls

@testset ExtendedTestSet "MetaBalls Tests" begin
    @includetests ["generators_tests", "edge_lookup_tests"]
end


text:
## add more julia processes
using Distributed
nprocs() < 4 ? addprocs(4-nprocs()) : nothing

# access modules/namespaces
using RoME
@everywhere using RoME
using RoMEPlotting

GTp = Dict{Symbol, Vector{Float64}}()
GTp[:l100] = [0.0;0]
GTp[:l101] = [50.0;0]
GTp[:l102] = [100.0;0]
GTp[:l103] = [100.0;50.0]
GTp[:l104] = [100.0;100.0]
GTp[:l105] = [50.0;100.0]
GTp[:l106] = [0.0;100.0]
GTp[:l107] = [0.0;50.0]
GTp[:l108] = [0.0;-50.0]
GTp[:l109] = [0.0;-100.0]
GTp[:l110] = [50.0;-100.0]
GTp[:l111] = [100.0;-100.0]
GTp[:l112] = [100.0;-50.0]

GTl = Dict{Symbol, Vector{Float64}}()
GTl[:l1] = [10.0;30]
GTl[:l2] = [30.0;-30]
GTl[:l3] = [80.0;40]
GTl[:l4] = [120.0;-50]

# create the factor graph object
fg = initfg()

# first pose with no initial estimate
addVariable!(fg, :l100, Point2)

# add three landmarks
addVariable!(fg, :l1, Point2)
addVariable!(fg, :l2, Point2)
addVariable!(fg, :l3, Point2)

# and put priors on :l101 and :l102
addFactor!(fg, [:l1;], PriorPoint2(MvNormal(GTl[:l1], Matrix{Float64}(LinearAlgebra.I,2,2))) )
addFactor!(fg, [:l2;], PriorPoint2(MvNormal(GTl[:l2], Matrix{Float64}(LinearAlgebra.I,2,2))) )

# first range measurement
rhoZ1 = norm(GTl[:l1]-GTp[:l100])
ppr = Point2Point2Range( Normal(rhoZ1, 2.0) )
addFactor!(fg, [:l100;:l1], ppr)

# second range measurement
rhoZ2 = norm(GTl[:l2]-GTp[:l100])
ppr = Point2Point2Range( Normal(rhoZ2, 3.0) )
addFactor!(fg, [:l100; :l2], ppr)

# second range measurement
rhoZ3 = norm(GTl[:l3]-GTp[:l100])
ppr = Point2Point2Range( Normal(rhoZ3, 3.0) )
addFactor!(fg, [:l100; :l3], ppr)

## solve system
tree, smt, hist = solveTree!(fg)

# plot the first results
plotKDE(fg, [:l1;:l2], dims=[1;2])

plotKDE(fg, :l100, dims=[1;2], levels=6)

drawLandms(fg, from=1, to=101) # |> PDF("/tmp/test.pdf"); #@async run(`evince /tmp/test.pdf`)

function vehicle_drives_to!(fgl::G, pos_sym::Symbol, GTp::Dict, GTl::Dict; measurelimit::R=150.0) where {G <: AbstractDFG, R <: Real}
  currvar = ls(fgl)
  prev_sym = Symbol("l$(maximum(Int[parse(Int,string(currvar[i])[2:end]) for i in 2:length(currvar)]))")
  if !(pos_sym in currvar)
    println("Adding variable vertex $pos_sym, not yet in fgl::FactorGraph.")
    addVariable!(fgl, pos_sym, Point2)
    @show rho = norm(GTp[prev_sym] - GTp[pos_sym])
    ppr = Point2Point2Range( Normal(rho, 3.0) )
    addFactor!(fgl, [prev_sym;pos_sym], ppr)
  else
    @warn "Variable node $pos_sym already in the factor graph."
  end
  beacons = keys(GTl)
  for ll in beacons
    rho = norm(GTl[ll] - GTp[pos_sym])
    # Check for feasible measurements:  vehicle within 150 units from the beacons/landmarks
    if rho < measurelimit
      ppr = Point2Point2Range( Normal(rho, 3.0) )
      if !(ll in currvar)
        println("Adding variable vertex $ll, not yet in fgl::FactorGraph.")
        addVariable!(fgl, ll, Point2)
      end
      addFactor!(fgl, [pos_sym;ll], ppr)
    end
  end
  nothing
end

vehicle_drives_to!(fg, :l101, GTp, GTl)
vehicle_drives_to!(fg, :l102, GTp, GTl)

tree, smt, hist = solveTree!(fg)

pl = plotKDE(fg, [Symbol("l$(100+i)") for i in 0:2], dims=[1;2])

pl = plotKDE(fg, [:l3;:l4], dims=[1;2], levels=4)

vehicle_drives_to!(fg, :l103, GTp, GTl)
vehicle_drives_to!(fg, :l104, GTp, GTl)

getSolverParams(fg).drawtree = true
getSolverParams(fg).showtree = true

tree, smt, hist = solveTree!(fg)

pl = plotKDE(fg, [Symbol("l$(100+i)") for i in 0:4], dims=[1;2]) |> PDF("/tmp/test.pdf")
@async run(`evince /tmp/test.pdf`)

# drive further
vehicle_drives_to!(fg, :l105, GTp, GTl)
vehicle_drives_to!(fg, :l106, GTp, GTl)

# allow potential for incremental updates
tree, smt, hist = solveTree!(fg, tree)

pl = plotKDE(fg, [Symbol("l$(100+i)") for i in 0:4], dims=[1;2]) |> PDF("/tmp/test.pdf")
@async run(`evince /tmp/test.pdf`)

vehicle_drives_to!(fg, :l107, GTp, GTl)

# do a batch solve again
tree, smt, hist = solveTree!(fg)

vehicle_drives_to!(fg, :l108, GTp, GTl)

# and incremental again
tree, smt, hist = solveTree!(fg, tree)

pl = plotKDE(fg, [Symbol("l$(100+i)") for i in 2:8], dims=[1;2], levels=6)
pl |> PDF("/tmp/test.pdf"); @async run(`evince /tmp/test.pdf`)

vehicle_drives_to!(fg, :l109, GTp, GTl)
vehicle_drives_to!(fg, :l110, GTp, GTl)

tree, smt, hist = solveTree!(fg, tree)

vehicle_drives_to!(fg, :l111, GTp, GTl)
vehicle_drives_to!(fg, :l112, GTp, GTl)

tree, smt, hist = solveTree!(fg, tree)

pl = plotKDE(fg, [Symbol("l$(100+i)") for i in 7:12], dims=[1;2])
pl |> PDF("/tmp/test.pdf"); @async run(`evince /tmp/test.pdf`)

pl = plotKDE(fg, [:l1;:l2;:l3;:l4], dims=[1;2], levels=4)
pl |> PDF("/tmp/test.pdf"); @async run(`evince /tmp/test.pdf`)

pl = drawLandms(fg)
pl |> PDF("/tmp/test.pdf"); @async run(`evince /tmp/test.pdf`)

drawTree(tree, show=true, imgs=true)


text:
# internal objects, placed to ASTs in compiling time
module Runtime
using MLStyle
using CamJulia.IdrisList
using CamJulia.IdrisForeignCollections

struct FileHandler
    filename :: String
    mode     :: String
end

simple_open(filename::IdrisString) = FileHandler(string(filename), "r")
simple_read(f::FileHandler) = open(f.filename) do f
    from_text(read(f, String))
end

read_all_text(f::IOStream) = read(f, String)

filesystem_pipe(::Nothing) = (1, 2)

@inline function is(a :: T , b :: T) where T
    a === b
end

is(_, _) = false

@inline function to_list_str(s)
    from_text(string(s))
end

@inline function parse_from_list_str(::Type{T}, s::IdrisString) where T
    parse(T, string(s))
end

function parse_int_from_list_str(s::IdrisString)
    parse_from_list_str(Int, s)
end

function parse_double_from_list_str(s::IdrisString)
    parse_from_list_str(Float64, s)
end

@inline function module_property(m :: Module, s::String)
    getproperty(m, Symbol(s))
end

export rt_support
rt_support = Dict{String, Any}(
    "cam-rt.cmp" => (==),
    "cam-rt.err" => throw,
    "cam-rt.is" => is,

    "prim-plus" => (+),
    "prim-minus" => (-),
    "prim-times" => (*),
    "prim-udiv" => div,
    "prim-sdiv" => div,
    "prim-urem" => rem,
    "prim-srem" => rem,
    "prim-eq" => (==),
    "prim-slt" => (<),
    "prim-sle" => (<=),
    "prim-sgt" => (>),
    "prim-sge" => (>=),
    "prim-and" => (&),
    "prim-or" => (|),

    # str methods1
    "prim-streq" => (==),
    "prim-strlen" => string_len,
    "prim-strindex" => getindex,
    "prim-strlt" => (<=),

    # effect
    "prim-external" => throw,
    "prim-readstr" => readline,
    "prim-writestr" => print,

    # conversion
    "prim-floatstr" => to_list_str,
    "prim-strfloat" => x -> parse(Float64, x) ∘ string,
    "prim-intstr" => to_list_str,
    "prim-strint" => x -> parse(Int, x) ∘ string,
    "prim-intch" => Char,
    "prim-chint" => Int,

    # str method2
    "prim-strhead" => string_head,
    "prim-strtail" => string_tail,
    "prim-strcons" => string_cons,
    "prim-strconcat" => string_concat,
    "prim-crash" => throw,

    # io
    "builtin-println" => println,
    "builtin-simple_open" => simple_open,
    "builtin-simple_read" => simple_read,
    "builtin-filesystem_pipe" => filesystem_pipe,
    "builtin-filesystem_open_file" => open,
    "builtin-filesystem_close_file" => close,
    "builtin-filesystem_read_all_text" => read_all_text,

    # flist
    "builtin-reverse_flist" => reverse,
    "builtin-flist_to_native" => from_flist,
    "builtin-list_to_foreign" => to_flist,

    # fvect
    "builtin-reverse_fvect" => reverse,
    "builtin-fvect_to_native" => from_fvect,
    "builtin-vect_to_foreign" => to_fvect,

    # fhvect
    "builtin-reverse_fhvect" => reverse,
    "builtin-fhvect_to_native" => from_fhvect,
    "builtin-hvect_to_foreign" => to_fhvect,

    # string
    "builtin-to_text" => string,
    "builtin-fstr_to_native" => from_text,
    "builtin-str_to_foreign" => string,

    "builtin-map_hvect" => map_hvect,
    "builtin-module_property" => module_property,

)
end

max_stars_repo_name: thautwarm/idris-cam

max_issues_repo_name: thautwarm/idris-cam

max_forks_repo_name: thautwarm/idris-cam


text:
#=
P12 (**) Decode a run-length encoded list.
Given a run-length code list generated as specified in problem P11. Construct its uncompressed version.
=#

function my_repeat(s, num)
    function _my_repeat(s, num, repeated)
        num == 0 ? repeated : _my_repeat(s, num - 1, vcat(repeated, s))
    end
    _my_repeat(s, num, [])
end

function decode(coll)
    pairs = map(item -> typeof(item) == Symbol ? (1, item) : item, coll)
    repeated = map(pair -> my_repeat(pair[2], pair[1]), pairs)
    reduce(vcat, repeated)
end

@assert decode(((4, :a), :b, (2, :c), (2, :a), :d, (4, :e))) == [:a, :a, :a, :a, :b, :c, :c, :a, :a, :d, :e, :e, :e, :e]

println("Tests passed: JL-12.jl")


text:
function RWS_Example_5_1()
    Al = Alphabet(['a', 'A', 'b', 'B'])
    KnuthBendix.set_inversion!(Al, 'a', 'A')
    KnuthBendix.set_inversion!(Al, 'b', 'B')

    a,A,b,B = Word.([i] for i in 1:4)
    ε = one(a)

    R = RewritingSystem([b*a=>a*b], LenLex(Al))

    return R
end

RWS_ZxZ() = RWS_Example_5_1()

function RWS_Example_5_2()
    Al = Alphabet(['a', 'b', 'B', 'A'])
    KnuthBendix.set_inversion!(Al, 'a', 'A')
    KnuthBendix.set_inversion!(Al, 'b', 'B')

    a,b,B,A = Word.([i] for i in 1:4)
    ε = one(a)

    R = RewritingSystem([b*a=>a*b], LenLex(Al))

    return R
end

RWS_ZxZ_nonterminating() = RWS_Example_5_2()

function RWS_Example_5_3()
    Al = Alphabet(['a', 'b'])

    a,b = Word.([i] for i in 1:2)
    ε = one(a)

    R = RewritingSystem([a^2=>ε, b^3=>ε, (a*b)^3=>ε], LenLex(Al))

    return R
end

function RWS_Example_5_4()
    Al = Alphabet(['a', 'b', 'B'])
    KnuthBendix.set_inversion!(Al, 'b', 'B')

    a,b,B = Word.([i] for i in 1:3)
    ε = one(a)

    R = RewritingSystem([a^2=>ε, b^3=>ε, (a*b)^3=>ε], LenLex(Al))

    return R
end

function RWS_Example_5_5()
    Al = Alphabet(['c', 'C', 'b', 'B', 'a', 'A'])
    KnuthBendix.set_inversion!(Al, 'c', 'C')
    KnuthBendix.set_inversion!(Al, 'b', 'B')
    KnuthBendix.set_inversion!(Al, 'a', 'A')

    c, C, b, B, a, A = Word.([i] for i in 1:6)
    ε = one(a)

    eqns = [c*a =>a*c, c*b=>b*c, b*a=>a*b*c]

    R = RewritingSystem(eqns, WreathOrder(Al))
    return R
end

function RWS_Example_6_4()
    Al = Alphabet(['a', 'b', 'B'])
    KnuthBendix.set_inversion!(Al, 'b', 'B')

    a, b, B = Word.([i] for i in 1:3)
    ε = one(a)
    eqns = [a*a=>ε, b*B=>ε, b^3=>ε, (a*b)^7=>ε, (a*b*a*B)^4=>ε]

    R = RewritingSystem(eqns, LenLex(Al))
    return R
end

function RWS_Example_6_5()
    Al = Alphabet(['a', 'b', 'B'])
    KnuthBendix.set_inversion!(Al, 'b', 'B')

    a, b, B = Word.([i] for i in 1:3)
    ε = one(a)
    eqns = KnuthBendix.Rule.(
        [a*a=>ε, b*B=>ε, b^2=>B, (B*a)^3*B=>(a*b)^3*a, (b*a*B*a)^2=>(a*b*a*B)^2]
    )

    R = RewritingSystem(eqns, LenLex(Al))
    return R
end

function RWS_Example_237_abaB(n)
    Al = Alphabet(['a', 'b', 'B'])
    KnuthBendix.set_inversion!(Al, 'b', 'B')

    a, b, B = Word.([i] for i in 1:3)
    ε = one(a)

    # eqns = [a^2=>ε, b^3=>ε, (a*b)^7=>ε, (a*b*a*B)^n=>ε]
    eqns = [b*B=>ε, B*b=>ε, a^2=>ε, b^3=>ε, (a*b)^7=>ε, (a*b*a*B)^n=>ε]

    R = RewritingSystem(eqns, LenLex(Al), bare=true)
    return R
end

RWS_Example_6_6() = RWS_Example_237_abaB(8)

function RWS_Exercise_6_1(n)
    Al = Alphabet(['a', 'b'])

    a, b = Word.([i] for i in 1:2)
    ε = one(a)

    eqns = [a^2=>ε, b^3=>ε, (a*b)^n=>ε]

    R = RewritingSystem(eqns, LenLex(Al), bare=true)
    return R
end

function RWS_Closed_Orientable_Surface(n)
    ltrs = String[]
    for i in 1:n
        subscript = join('₀'+d for d in reverse(digits(i)))
        append!(ltrs, ["a" * subscript, "A" * subscript, "b" * subscript, "B" * subscript])
    end
    Al = Alphabet(reverse!(ltrs))
    for i in 1:n
        subscript = join('₀'+d for d in reverse(digits(i)))
        KnuthBendix.set_inversion!(Al, "a" * subscript, "A" * subscript)
        KnuthBendix.set_inversion!(Al, "b" * subscript, "B" * subscript)
    end

    ε = one(Word{UInt16})
    rules = Pair{typeof(ε), typeof(ε)}[]
    word = Int[]

    for i in reverse(1:n)
        x = 4 * i
        append!(word, [x, x-2, x-1, x-3])
    end
    push!(rules, Word(word) => ε)
    R = RewritingSystem(rules, RecursivePathOrder(Al))

    return R
end


text:
# ---
# title: 197. Rising Temperature
# id: problem197
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories:
# link: <https://leetcode.com/problems/rising-temperature/description/>
# hidden: true
# ---
#
# Table: `Weather`
#
#
#
#     +---------------+---------+
#     | Column Name   | Type    |
#     +---------------+---------+
#     | id            | int     |
#     | recordDate    | date    |
#     | temperature   | int     |
#     +---------------+---------+
#     id is the primary key for this table.
#     This table contains information about the temperature in a certain day.
#
#
#
#
# Write an SQL query to find all dates' `id` with higher temperature compared to
# its previous dates (yesterday).
#
# Return the result table in **any order**.
#
# The query result format is in the following example:
#
#
#
#     Weather
#     +----+------------+-------------+
#     | id | recordDate | Temperature |
#     +----+------------+-------------+
#     | 1  | 2015-01-01 | 10          |
#     | 2  | 2015-01-02 | 25          |
#     | 3  | 2015-01-03 | 20          |
#     | 4  | 2015-01-04 | 30          |
#     +----+------------+-------------+
#
#     Result table:
#     +----+
#     | id |
#     +----+
#     | 2  |
#     | 4  |
#     +----+
#     In 2015-01-02, temperature was higher than the previous day (10 -> 25).
#     In 2015-01-04, temperature was higher than the previous day (30 -> 20).
#
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
using NLPModels
using LinearOperators
using Plots # graph
Plots.pyplot()
using BenchmarkTools
using BenchmarkProfiles
using MiniLogging

# basic configuration. The root logger level is then INFO
basic_config(MiniLogging.INFO; date_format="%Y-%m-%d %H:%M:%S")
# loggerCRlin = MiniLogging.Logger("CRlin", "CRlin.out", MiniLogging.DEBUG)
# loggerCRlin = MiniLogging.Logger("CRlin", "CRlin.out", MiniLogging.DEBUG)

loggerCRlin = get_logger("CRlin")
loggerCGlin = get_logger("CGlin")
loggerArm = get_logger("Armijo")
loggerLin = get_logger("Linesearch")
loggerlaunch = get_logger("launch_linesearch.jl")
loggerLinCUTEst = get_logger("LinesearchCUTEst")

loggerCRlin.level = MiniLogging.ERROR
loggerCGlin.level = MiniLogging.ERROR
loggerArm.level = MiniLogging.ERROR
loggerLin.level = MiniLogging.ERROR
# loggerlaunch.level = MiniLogging.ERROR
loggerLinCUTEst.level = MiniLogging.ERROR

include("CGlin.jl")
include("CRlin.jl")
include("Armijo.jl")
include("Linesearch.jl")
include("LCG.jl")
include("LCR.jl")
include("LinesearchCUTEst.jl")


text:
include("../examples/rijke.jl")
include("../src/lss.jl")
using LinearAlgebra
using JLD
using OrdinaryDiffEq
function refine_parameter_and_trajectory(n, n_gd_steps)
    #create epsilon orbit.
    n_runup = 400000
    s = [7.0,0.2]
    d = N
    u0 = Rijke_ODE(rand(d), s, n_runup)
    eps = 0.05
    noise = eps*randn(d)
    u_trj = zeros(d,n)
    u_trj[:,1] = u0 + noise
    u_obs = zeros(d,n)
    u_obs[:,1] = u0
    for i = 2:n
        u_trj[:,i] = Rijke_ODE(u_trj[:,i-1], s, 1)
        u_obs[:,i] = Rijke_ODE(u_obs[:,i-1], s, 1)
    end
    z_prd = zeros(n, n_gd_steps)
    z_obs = sum(u_obs[1:2*Ng,:].*u_obs[1:2*Ng,:],dims=1)./4
    z_trj = sum(u_trj[1:2*Ng,:].*u_trj[1:2*Ng,:],dims=1)./4
    z_prd[:,1] = z_trj
    # run gradient descent to minimize observation error.
    gamma = 0.1
    dJds = zeros(n_gd_steps)
    d_u = 2
    # compute_sens assumes many objective functions
    dJ = zeros(1,d,n)
    error = zeros(n_gd_steps)
    n_pert = 100
    flag = 1

    du_trj = zeros(d,d,n)
    f_trj = zeros(d,n)
    X_trj = zeros(d,n)
    @time for i = 1:n_gd_steps
        if flag == 0
            break
        end
        # Set up LSS
        # J is now mean of squared observation error.
        for k = 1:n
            z_prd[k,i] = sum(x->x*x, u_trj[1:2*Ng,k])/4
            error_k = (z_obs[k] .- z_prd[k,i])^2.0/n
            dJ[1,1:2*Ng,k] = -1/n*error_k*u_trj[1:2*Ng,k]
            du_trj[:,:,k] = dRijke(u_trj[:,k], s, 1.e-6)
            X_trj[:,k] = perturbation(u_trj[:,k],s,1.e-6) #ith col in T_{u_{i+1}} M
            f!(view(f_trj,:,k), u_trj[:,k], s, 1.)
            #f = zeros(d,n)
        end
        y, xi = lss(du_trj, X_trj, f_trj, s, d_u)
        dJds[i] = compute_sens(y, xi, dJ, f_trj)[1]
        #println(dJds[i])
        # Make u_trj the computed shadowing trajectory
        ds = gamma*dJds[i]
        s[1] = s[1] - ds
        #u_trj = u_trj - ds*y'
        u2 = u_trj[:,n_pert] .- ds*y[:,n_pert]
        for k = n_pert+1:n
            u_trj[:,k] = Rijke_ODE(u_trj[:,k-1], s, 1)
            error[i] += (z_prd[k,i] - z_obs[k])^2.0/(n-n_pert)
        end
        @show error[i]
        if i==2
            if error[i] > error[i-1]
                flag = 0
                break
            end
        end
    end
    return z_obs, z_prd, error, flag
end
function assimilate_parameter_and_trajectory(filename)
    n_repeat = 1
    n = 2000
    n_gd_steps = 200
    z_obs = zeros(n, n_repeat)
    z_trj = zeros(n, n_gd_steps, n_repeat)
    errors = zeros(n_gd_steps, n_repeat)
    i = 1
    while i <= n_repeat
        @show i
        res1, res2, res3, flag = refine_parameter_and_trajectory(n,
                                            n_gd_steps)
        if flag==1
            z_obs[:,i] = res1
            z_trj[:,:,i] = res2
            errors[:,i] = res3
            i = i + 1
        end
    end
    #println("Errors are ")
    #@show errors[:,1]
    save(filename,
         "z_obs", z_obs, "z_prd", z_trj, "msq_err",
         errors)

    #return z_obs, z_trj, errors
end

max_stars_repo_name: nishaChandramoorthy/energies

max_issues_repo_name: nishaChandramoorthy/energies

max_forks_repo_name: nishaChandramoorthy/energies