text:
using Colors
using LinearAlgebra
using Luxor

# Auxiliary drawing functions
function make_drawing(width, height, img_path, bkg_color, origin_p)
    d = Drawing(width, height, img_path)
    background(bkg_color)
    origin(origin_p)
    return d
end

map2luxor(p) = Point(p.x, -p.y)
map2luxor(p) = Point(p...)

# Constants

## Canvas constants
width = 500
height = 500
path = "jarvis_march.png"
color = RGBA(0, 0, 0)
op = [width / 2, height / 2]

# Setting up drawing
my_draw = make_drawing(width, height, path, color, Point(op...))

# Points

p1 = [125, 125]
p2 = [170, -200]
p3 = [-200, 0] #leftmost point
p4 = [-30, -140]
p5 = [0, -5]
p6 = [200, 150]
p7 = [100, 230]
p8 = [220, 0]
p9 = [-100, 80]
p10 = [-90, -200]

sethue(RGBA(1.0, 0.9411764705882353, 0.0))
# pts_list = [p1, p2, p3, p4, p5, p6, p7, p8, p9, p10] .|> x -> Point(x...)
pts_list = [Point(rand(-250:250), rand(-250:250)) for _ in 1:10]
map2luxor.(pts_list) .|> cp -> circle(cp, 15; action = :fill)

sethue(RGBA(1.0, 0.0, 0.0))
leftmost_pt = sort!(pts_list, by = pt -> pt.x) |> x -> x[1]
map2luxor(leftmost_pt) |> x -> circle(x, 15; action = :fill)

sethue(RGBA(0.0, 0.0, 1))
circle(Point(0, 0), 15; action = :fill)

hull = [leftmost_pt]
let
    curr_point = pts_list[2]
    curr_product = crossproduct(leftmost_pt, curr_point) |> abs
    i = 1
    while curr_point != leftmost_pt
        for pt in pts_list
            product = 0
            if i == 1
                if hull[i] != pt
                    product = crossproduct(hull[i], pt)
                end
            else
                if hull[i] != pt && hull[i - 1] != pt
                    product = crossproduct(hull[i - 1], pt)
                end
            end
            if (product > curr_product)
                curr_point = pt
                curr_product = product
            end
        end
        push!(hull, curr_point)
        curr_product = 0
        i += 1
    end

end

sethue(RGBA(1.0, 1.0, 1.0))
poly(hull, action = :stroke; close = true)

finish()

max_stars_repo_name: TheCedarPrince/DataExplorers

max_issues_repo_name: TheCedarPrince/DataExplorers

max_forks_repo_name: TheCedarPrince/DataExplorers


text:
function closer(ps::ParseState)
    (ps.closer.newline && ps.ws.kind == NewLineWS && ps.t.kind != Tokens.COMMA) ||
    (ps.closer.semicolon && ps.ws.kind == SemiColonWS) ||
    (isoperator(ps.nt) && precedence(ps.nt) <= ps.closer.precedence) ||
    (ps.nt.kind == Tokens.WHERE && ps.closer.precedence == 5) ||
    (ps.closer.inwhere && ps.nt.kind == Tokens.WHERE) ||
    (ps.nt.kind == Tokens.LPAREN && ps.closer.precedence > 15) ||
    (ps.nt.kind == Tokens.LBRACE && ps.closer.precedence > 15) ||
    (ps.nt.kind == Tokens.LSQUARE && ps.closer.precedence > 15) ||
    (ps.nt.kind == Tokens.STRING && isemptyws(ps.ws) && ps.closer.precedence > 15) ||
    (ps.closer.precedence > 15 && ps.t.kind == Tokens.RPAREN && ps.nt.kind == Tokens.IDENTIFIER) ||
    (ps.closer.precedence > 15 && ps.t.kind == Tokens.RSQUARE && ps.nt.kind == Tokens.IDENTIFIER) ||
    (ps.nt.kind == Tokens.COMMA && ps.closer.precedence > 0) ||
    ps.nt.kind == Tokens.ENDMARKER ||
    (ps.closer.comma && iscomma(ps.nt)) ||
    (ps.closer.tuple && (iscomma(ps.nt) || isassignment(ps.nt))) ||
    (ps.nt.kind == Tokens.FOR && ps.closer.precedence > -1) ||
    (ps.closer.paren && ps.nt.kind == Tokens.RPAREN) ||
    (ps.closer.brace && ps.nt.kind == Tokens.RBRACE) ||
    (ps.closer.square && ps.nt.kind == Tokens.RSQUARE) ||
    (ps.closer.block && ps.nt.kind == Tokens.END) ||
    (ps.closer.ifelse && ps.nt.kind == Tokens.ELSEIF || ps.nt.kind == Tokens.ELSE) ||
    (ps.closer.ifop && isoperator(ps.nt) && (precedence(ps.nt) <= 0 || ps.nt.kind == Tokens.COLON)) ||
    (ps.closer.trycatch && (ps.nt.kind == Tokens.CATCH || ps.nt.kind == Tokens.FINALLY || ps.nt.kind == Tokens.END)) ||
    (ps.closer.range && ps.nt.kind == Tokens.FOR) ||
    (ps.closer.ws && !isemptyws(ps.ws) &&
        !(ps.nt.kind == Tokens.COMMA) &&
        !(ps.t.kind == Tokens.COMMA) &&
        !(!ps.closer.inmacro && ps.nt.kind == Tokens.FOR) &&
        !(ps.nt.kind == Tokens.DO) &&
        # !((isbinaryop(ps.nt.kind) && (!isemptyws(ps.nws) || !isunaryop(ps.nt))) || (!ps.closer.wsop && isbinaryop(ps.nt.kind)))) ||
        !((isbinaryop(ps.nt) && !isemptyws(ps.nws)) ||
        (isbinaryop(ps.nt) && !isunaryop(ps.nt)) ||
        (isunaryop(ps.t) && ps.ws.kind == WS) ||
        (!ps.closer.wsop && isbinaryop(ps.nt.kind)))) ||
    (ps.nt.startbyte ≥ ps.closer.stop) ||
    ps.errored
end

"""
    @closer ps rule body

Continues parsing closing on `rule`.
"""
macro closer(ps, opt, body)
    quote
        local tmp1 = $(esc(ps)).closer.$opt
        $(esc(ps)).closer.$opt = true
        out = $(esc(body))
        $(esc(ps)).closer.$opt = tmp1
        out
    end
end

"""
    @nocloser ps rule body

Continues parsing not closing on `rule`.
"""
macro nocloser(ps, opt, body)
    quote
        local tmp1 = $(esc(ps)).closer.$opt
        $(esc(ps)).closer.$opt = false
        out = $(esc(body))
        $(esc(ps)).closer.$opt = tmp1
        out
    end
end

"""
    @precedence ps prec body

Continues parsing binary operators until it hits a more loosely binding
operator (with precdence lower than `prec`).
"""
macro precedence(ps, prec, body)
    quote
        local tmp1 = $(esc(ps)).closer.precedence
        $(esc(ps)).closer.precedence = $(esc(prec))
        out = $(esc(body))
        $(esc(ps)).closer.precedence = tmp1
        out
    end
end

"""
    @default ps body

Parses the next expression using default closure rules.
"""
macro default(ps, body)
    quote
        local tmp1 = $(esc(ps)).closer.newline
        local tmp2 = $(esc(ps)).closer.semicolon
        local tmp3 = $(esc(ps)).closer.inmacro
        local tmp4 = $(esc(ps)).closer.tuple
        local tmp5 = $(esc(ps)).closer.comma
        local tmp6 = $(esc(ps)).closer.insquare
        # local tmp7 = $(esc(ps)).closer.brace
        local tmp8 = $(esc(ps)).closer.range
        local tmp9 = $(esc(ps)).closer.block
        local tmp10 = $(esc(ps)).closer.ifelse
        local tmp11 = $(esc(ps)).closer.ifop
        # local tmp12 = $(esc(ps)).closer.trycatch
        local tmp13 = $(esc(ps)).closer.ws
        local tmp14 = $(esc(ps)).closer.wsop
        local tmp15 = $(esc(ps)).closer.precedence
        $(esc(ps)).closer.newline = true
        $(esc(ps)).closer.semicolon = true
        $(esc(ps)).closer.inmacro = false
        $(esc(ps)).closer.tuple = false
        $(esc(ps)).closer.comma = false
        $(esc(ps)).closer.insquare = false
        # $(esc(ps)).closer.brace = false
        # $(esc(ps)).closer.square = false
        $(esc(ps)).closer.range = false
        $(esc(ps)).closer.ifelse = false
        $(esc(ps)).closer.ifop = false
        # $(esc(ps)).closer.trycatch = false
        $(esc(ps)).closer.ws = false
        $(esc(ps)).closer.wsop = false
        $(esc(ps)).closer.precedence = -1

        out = $(esc(body))

        $(esc(ps)).closer.newline = tmp1
        $(esc(ps)).closer.semicolon = tmp2
        $(esc(ps)).closer.inmacro = tmp3
        $(esc(ps)).closer.tuple = tmp4
        $(esc(ps)).closer.comma = tmp5
        $(esc(ps)).closer.insquare = tmp6
        # $(esc(ps)).closer.brace = tmp7
        $(esc(ps)).closer.range = tmp8
        $(esc(ps)).closer.block = tmp9
        $(esc(ps)).closer.ifelse = tmp10
        $(esc(ps)).closer.ifop = tmp11
        # $(esc(ps)).closer.trycatch = tmp12
        $(esc(ps)).closer.ws = tmp13
        $(esc(ps)).closer.wsop = tmp14
        $(esc(ps)).closer.precedence = tmp15
        out
    end
end

"""
    @scope ps scope body

Continues parsing tracking declared variables.
"""
macro scope(ps, new_scope, body)
    quote
        local tmp1 = $(esc(ps)).current_scope
        $(esc(ps)).current_scope = $(esc(new_scope))
        out = $(esc(body))
        $(esc(ps)).current_scope = tmp1
        out
    end
end

"""
    @noscope ps body

Continues parsing not tracking declared variables.
"""
macro noscope(ps, body)
    quote
        local tmp1 = $(esc(ps)).trackscope
        $(esc(ps)).trackscope = false
        out = $(esc(body))
        $(esc(ps)).trackscope = tmp1
        out
    end
end

"""
    @catcherror ps body

Checks for `ps.errored`.
"""
macro catcherror(ps, startbyte, body)
    quote
        $(esc(body))
        if $(esc(ps)).errored
            return EXPR{ERROR}(EXPR[INSTANCE($(esc(ps)))], 0, Variable[], "Unknown error")
        end
    end
end

isidentifier(t::Token) = t.kind == Tokens.IDENTIFIER

isliteral(t::Token) = Tokens.begin_literal < t.kind < Tokens.end_literal

isbool(t::Token) =  Tokens.TRUE ≤ t.kind ≤ Tokens.FALSE
iscomma(t::Token) =  t.kind == Tokens.COMMA

iskw(t::Token) = Tokens.iskeyword(t.kind)

isinstance(t::Token) = isidentifier(t) ||
                       isliteral(t) ||
                       isbool(t) ||
                       iskw(t)

ispunctuation(t::Token) = t.kind == Tokens.COMMA ||
                          t.kind == Tokens.END ||
                          Tokens.LSQUARE ≤ t.kind ≤ Tokens.RPAREN

isstring(x) = false
isstring(x::EXPR{T}) where T <: Union{StringH, LITERAL{Tokens.STRING},LITERAL{Tokens.TRIPLE_STRING}} = true

isajuxtaposition(ps::ParseState, ret) = ((ret isa EXPR{LITERAL{Tokens.INTEGER}} || ret isa EXPR{LITERAL{Tokens.FLOAT}}) && (ps.nt.kind == Tokens.IDENTIFIER || ps.nt.kind == Tokens.LPAREN || ps.nt.kind == Tokens.CMD || ps.nt.kind == Tokens.STRING || ps.nt.kind == Tokens.TRIPLE_STRING)) || (
        (ret isa EXPR{UnarySyntaxOpCall} && ret.args[2] isa EXPR{OPERATOR{16,Tokens.PRIME,false}} && ps.nt.kind == Tokens.IDENTIFIER) ||
        ((ps.t.kind == Tokens.RPAREN || ps.t.kind == Tokens.RSQUARE) && (ps.nt.kind == Tokens.IDENTIFIER || ps.nt.kind == Tokens.CMD)) ||
        ((ps.t.kind == Tokens.STRING || ps.t.kind == Tokens.TRIPLE_STRING) && (ps.nt.kind == Tokens.STRING || ps.nt.kind == Tokens.TRIPLE_STRING))) ||
        (isstring(ret) && ps.nt.kind == Tokens.IDENTIFIER && ps.ws.kind == EmptyWS)

# Testing functions

function test_order(x, out = [])
    if x isa EXPR
        for y in x
            test_order(y, out)
        end
    else
        push!(out, x)
    end
    out
end

function test_find(str)
    x = parse(str, true)
    for i = 1:sizeof(str)
        _find(x, i)
    end
end

"""
    compare(x,y)

Recursively checks whether two Base.Expr are the same. Returns unequal sub-
expressions.
"""
compare(x, y) = x == y ? true : (x, y)

function compare(x::Expr, y::Expr)
    if x == y
        return true
    else
        if x.head != y.head
            return (x, y)
        end
        if length(x.args) != length(y.args)
            return (x.args, y.args)
        end
        for i = 1:length(x.args)
            t = compare(x.args[i], y.args[i])
            if t != true
                return t
            end
        end
    end
end

"""
    span(x, neq = [])

Recursively checks whether the span of an expression equals the sum of the span
of its components. Returns a vector of failing expressions.
"""
function span(x::EXPR{StringH}, neq = []) end
function span(x, neq = [])
    s = 0
    for a in x.args
        span(a, neq)
        s += a.span
    end
    if length(x.args) > 0 && s != x.span
        push!(neq, x)
    end
    neq
end

function check_base(dir = dirname(Base.find_source_file("base.jl")), display = false)
    N = 0
    neq = 0
    err = 0
    aerr = 0
    fail = 0
    bfail = 0
    ret = []
    oldstderr = STDERR
    redirect_stderr()
    for (rp, d, files) in walkdir(dir)
        for f in files
            file = joinpath(rp, f)
            if endswith(file, ".jl")
                N += 1
                try
                    # print(N)
                    print("\r", rpad(string(N), 5), rpad(string(signif(fail / N * 100, 3)), 8), rpad(string(signif(err / N * 100, 3)), 8), rpad(string(signif(neq / N * 100, 3)), 8))
                    str = readstring(file)
                    ps = ParseState(str)
                    io = IOBuffer(str)
                    x, ps = parse(ps, true)
                    sp = span(x)
                    if length(x.args) > 0 && x.args[1] isa EXPR{LITERAL{nothing}}
                        shift!(x.args)
                    end
                    if length(x.args) > 0 && x.args[end] isa EXPR{LITERAL{nothing}}
                        pop!(x.args)
                    end
                    x0 = Expr(x)
                    x1 = Expr(:file)
                    try
                        while !eof(io)
                            push!(x1.args, Base.parse(io))
                        end
                    catch er
                        bfail += 1
                        continue
                    end
                    if length(x1.args) > 0  && x1.args[end] == nothing
                        pop!(x1.args)
                    end
                    remlineinfo!(x1)
                    print("\r                             ")
                    if !isempty(sp)
                        print_with_color(:blue, file)
                        println()
                        push!(ret, (file, :span))
                    end
                    if ps.errored
                        err += 1
                        print_with_color(:yellow, file)
                        println()
                        push!(ret, (file, :errored))
                    elseif !(x0 == x1)
                        cumfail = 0
                        neq += 1
                        print_with_color(:green, file)
                        println()
                        if display
                            c0, c1 = compare(x0, x1)
                            if !(c0 isa String && c1 isa String)
                                aerr += 1
                                print_with_color(:light_red, string("    ", c0), bold = true)
                                println()
                                print_with_color(:light_green, string("    ", c1), bold = true)
                                println()
                            end

                        end
                        push!(ret, (file, :noteq))
                    end

                catch er
                    fail += 1
                    print_with_color(:red, file)
                    println()
                    push!(ret, (file, :failed))
                end
            end
        end
    end
    redirect_stderr(oldstderr)
    if fail + err + neq > 0
        println("\r$N files")
        print_with_color(:red, "failed")
        println(" : $fail    $(100*fail/N)%")
        print_with_color(:yellow, "errored")
        println(" : $err     $(100*err/N)%")
        print_with_color(:green, "not eq.")
        println(" : $neq    $(100*neq/N)%", "  -  $aerr     $(100*aerr/N)%")
        print_with_color(:magenta, "base failed")
        println(" : $bfail    $(100*bfail/N)%")
    end
    ret
end

function speed_test()
    dir = dirname(Base.find_source_file("base.jl"))
    println("speed test : ", @timed(for i = 1:5
    parse(readstring(joinpath(dir, "inference.jl")), true);
    parse(readstring(joinpath(dir, "random.jl")), true);
    parse(readstring(joinpath(dir, "show.jl")), true);
    parse(readstring(joinpath(dir, "abstractarray.jl")), true);
end)[2])
end

"""
    check_reformat()

Reads and parses all files in current directory, applys formatting fixes and checks that the output AST remains the same.
"""
function check_reformat()
    fs = filter(f -> endswith(f, ".jl"), readdir())
    for (i, f) in enumerate(fs)
        f == "deprecated.jl" && continue
        str = readstring(f)
        x, ps = parse(ParseState(str), true);
        cnt = 0
        for i = 1:length(x) - 1
            y = x[i]
            sstr = str[cnt + (1:y.span)]

            y1 = parse(sstr)
            @assert Expr(y) == Expr(y1)
            cnt += y.span
        end
    end
end

is_func_call(x) = false
is_func_call(x::EXPR) = false
is_func_call(x::EXPR{Call}) = true
is_func_call(x::EXPR{UnaryOpCall}) = true
function is_func_call(x::EXPR{BinarySyntaxOpCall})
    if length(x.args) > 1 && (x.args[2] isa EXPR{OPERATOR{WhereOp,Tokens.WHERE,false}} || x.args[2] isa EXPR{OPERATOR{DeclarationOp,Tokens.DECLARATION,false}})
        return is_func_call(x.args[1])
    else
        return false
    end
end

function get_last_token(x::CSTParser.EXPR)
    if isempty(x.args)
        return x
    else
        return get_last_token(last(x.args))
    end
end

function trailing_ws_length(x::CSTParser.EXPR{CSTParser.IDENTIFIER})
    x.span - sizeof(x.val)
end

function trailing_ws_length(x::CSTParser.EXPR{P}) where P <: CSTParser.PUNCTUATION
    x.span - 1
end

function trailing_ws_length(x::CSTParser.EXPR{K}) where K <: CSTParser.KEYWORD{T} where T
    x.span - sizeof(string(T))
end

function trailing_ws_length(x::CSTParser.EXPR{K}) where K <: CSTParser.LITERAL{T} where T
    x.span - sizeof(x.val)
end

ext: jl


text:
function run_pandamodels_qflex_test(json_path)   # before run_poweramodels
    time_start = time()
    ###############################################################################
    # 0. Initialization
    ###############################################################################
    pm = load_pm_from_json(json_path)
    solver = get_solver(pm["pm_solver"], pm["pm_nl_solver"], pm["pm_mip_solver"],
    pm["pm_log_level"], pm["pm_time_limit"], pm["pm_nl_time_limit"], pm["pm_mip_time_limit"])
    if haskey(pm, "user_defined_params")
        user_defined_params = pm["user_defined_params"]
        data = delete!(pm, "user_defined_params")
    else
        data = pm
    end

    # additional modification
    for (i, branch) in pm["branch"]
        pm["branch"][i]["angmin"] = -6.28
        pm["branch"][i]["angmax"] = 6.28
    end

    # use build_ref to filter out inactive components
    ref = PowerModels.build_ref(data)[:it][:pm][:nw][0]

    ###############################################################################
    # 1. Building the Optimal Power Flow Model
    ###############################################################################
    # Initialize a JuMP Optimization Model
    #-------------------------------------
    model = JuMP.Model(solver)

    # Add voltage angles va for each bus
    JuMP.@variable(model, va[i in keys(ref[:bus])])
    # note: [i in keys(ref[:bus])] adds one `va` variable for each bus in the network

    # Add voltage angles vm for each bus
    JuMP.@variable(model, ref[:bus][i]["vmin"] <= vm[i in keys(ref[:bus])] <= ref[:bus][i]["vmax"], start=1.0)
    # note: this vairable also includes the voltage magnitude limits and a starting value

    # Add active power generation variable pg for each generator (including limits)
    JuMP.@variable(model, ref[:gen][i]["pmin"] <= pg[i in keys(ref[:gen])] <= ref[:gen][i]["pmax"], start=ref[:gen][i]["pg"])
    # Add reactive power generation variable qg for each generator (including limits)
    JuMP.@variable(model, ref[:gen][i]["qmin"] <= qg[i in keys(ref[:gen])] <= ref[:gen][i]["qmax"], start=ref[:gen][i]["qg"])

    # Add power flow variables p to represent the active power flow for each branch
    JuMP.@variable(model, -ref[:branch][l]["rate_a"] <= p[(l,i,j) in ref[:arcs]] <= ref[:branch][l]["rate_a"])
    # Add power flow variables q to represent the reactive power flow for each branch
    JuMP.@variable(model, -ref[:branch][l]["rate_a"] <= q[(l,i,j) in ref[:arcs]] <= ref[:branch][l]["rate_a"])
    # note: ref[:arcs] includes both the from (i,j) and the to (j,i) sides of a branch

    # Add Objective Function
    # ----------------------

    JuMP.@objective(model, Min, sum((q[(content["element_index"],
                                content["f_bus"],
                                content["t_bus"])] - content["value"])^2
                            for (i, content) in user_defined_params["setpoint_q"]))

    # Add Constraints
    for (i,bus) in ref[:ref_buses]
        JuMP.@constraint(model, va[i] == 0)
    end

    # Nodal power balance constraints
    for (i,bus) in ref[:bus]
        # Build a list of the loads and shunt elements connected to the bus i
        bus_loads = [ref[:load][l] for l in ref[:bus_loads][i]]
        bus_shunts = [ref[:shunt][s] for s in ref[:bus_shunts][i]]

        # Active power balance at node i
        JuMP.@constraint(model,
            sum(p[a] for a in ref[:bus_arcs][i]) +                  # sum of active power flow on lines from bus i +
            sum(p_dc[a_dc] for a_dc in ref[:bus_arcs_dc][i]) ==     # sum of active power flow on HVDC lines from bus i =
            sum(pg[g] for g in ref[:bus_gens][i]) -                 # sum of active power generation at bus i -
            sum(load["pd"] for load in bus_loads) -                 # sum of active load consumption at bus i -
            sum(shunt["gs"] for shunt in bus_shunts)*vm[i]^2        # sum of active shunt element injections at bus i
        )

        # Reactive power balance at node i
        JuMP.@constraint(model,
            sum(q[a] for a in ref[:bus_arcs][i]) +                  # sum of reactive power flow on lines from bus i +
            sum(q_dc[a_dc] for a_dc in ref[:bus_arcs_dc][i]) ==     # sum of reactive power flow on HVDC lines from bus i =
            sum(qg[g] for g in ref[:bus_gens][i]) -                 # sum of reactive power generation at bus i -
            sum(load["qd"] for load in bus_loads) +                 # sum of reactive load consumption at bus i -
            sum(shunt["bs"] for shunt in bus_shunts)*vm[i]^2        # sum of reactive shunt element injections at bus i
        )
    end

    # Branch power flow physics and limit constraints
    for (i,branch) in ref[:branch]
        # Build the from variable id of the i-th branch, which is a tuple given by (branch id, from bus, to bus)
        f_idx = (i, branch["f_bus"], branch["t_bus"])
        # Build the to variable id of the i-th branch, which is a tuple given by (branch id, to bus, from bus)
        t_idx = (i, branch["t_bus"], branch["f_bus"])
        # note: it is necessary to distinguish between the from and to sides of a branch due to power losses

        p_fr = p[f_idx]                     # p_fr is a reference to the optimization variable p[f_idx]
        q_fr = q[f_idx]                     # q_fr is a reference to the optimization variable q[f_idx]
        p_to = p[t_idx]                     # p_to is a reference to the optimization variable p[t_idx]
        q_to = q[t_idx]                     # q_to is a reference to the optimization variable q[t_idx]
        # note: adding constraints to p_fr is equivalent to adding constraints to p[f_idx], and so on

        vm_fr = vm[branch["f_bus"]]         # vm_fr is a reference to the optimization variable vm on the from side of the branch
        vm_to = vm[branch["t_bus"]]         # vm_to is a reference to the optimization variable vm on the to side of the branch
        va_fr = va[branch["f_bus"]]         # va_fr is a reference to the optimization variable va on the from side of the branch
        va_to = va[branch["t_bus"]]         # va_fr is a reference to the optimization variable va on the to side of the branch

        # Compute the branch parameters and transformer ratios from the data
        g, b = PowerModels.calc_branch_y(branch)
        tr, ti = PowerModels.calc_branch_t(branch)
        g_fr = branch["g_fr"]
        b_fr = branch["b_fr"]
        g_to = branch["g_to"]
        b_to = branch["b_to"]
        tm = branch["tap"]^2
        # note: tap is assumed to be 1.0 on non-transformer branches

        # AC Power Flow Constraints

        # From side of the branch flow
        JuMP.@NLconstraint(model, p_fr ==  (g+g_fr)/tm*vm_fr^2 + (-g*tr+b*ti)/tm*(vm_fr*vm_to*cos(va_fr-va_to)) + (-b*tr-g*ti)/tm*(vm_fr*vm_to*sin(va_fr-va_to)) )
        JuMP.@NLconstraint(model, q_fr == -(b+b_fr)/tm*vm_fr^2 - (-b*tr-g*ti)/tm*(vm_fr*vm_to*cos(va_fr-va_to)) + (-g*tr+b*ti)/tm*(vm_fr*vm_to*sin(va_fr-va_to)) )

        # To side of the branch flow
        JuMP.@NLconstraint(model, p_to ==  (g+g_to)*vm_to^2 + (-g*tr-b*ti)/tm*(vm_to*vm_fr*cos(va_to-va_fr)) + (-b*tr+g*ti)/tm*(vm_to*vm_fr*sin(va_to-va_fr)) )
        JuMP.@NLconstraint(model, q_to == -(b+b_to)*vm_to^2 - (-b*tr+g*ti)/tm*(vm_to*vm_fr*cos(va_fr-va_to)) + (-g*tr-b*ti)/tm*(vm_to*vm_fr*sin(va_to-va_fr)) )

        # Voltage angle difference limit
        JuMP.@constraint(model, va_fr - va_to <= branch["angmax"])
        JuMP.@constraint(model, va_fr - va_to >= branch["angmin"])

        # Apparent power limit, from side and to side
        JuMP.@constraint(model, p_fr^2 + q_fr^2 <= branch["rate_a"]^2)
        JuMP.@constraint(model, p_to^2 + q_to^2 <= branch["rate_a"]^2)
    end

    ###############################################################################
    # 3. Solve the Optimal Power Flow Model and Review the Results
    ###############################################################################
    JuMP.optimize!(model)

    ###############################################################################
    # 4. Create Result Dictionary such that the PowerModels Results can be used by pandapower
    ###############################################################################
    solution =  Dict{String,Any}()
    push!(solution, "baseMVA" => data["baseMVA"])
    push!(solution, "per_unit" => data["per_unit"])
    push!(solution, "gen" => data["gen"])
    push!(solution, "bus" => data["bus"])
    push!(solution, "branch" => data["branch"])
    push!(solution, "multinetwork" => false)
    push!(solution, "multiinfrastructrue" => false)

    for (i, gen) in solution["gen"]
        index = gen["index"]
        gen["qg"] = value(model[:qg][index])
        gen["pg"] = value(model[:pg][index])
    end

    for (i, bus) in solution["bus"]
        index = bus["index"]
        bus["vm"] = value(model[:vm][index])
        bus["va"] = value(model[:va][index])
    end

    for (i, branch) in solution["branch"]
        index = branch["index"]
        push!(branch, "qf" => value(model[:q][(index, branch["f_bus"], branch["t_bus"])]))
        push!(branch, "qt" => value(model[:q][(index, branch["t_bus"], branch["f_bus"])]))
        push!(branch, "pf" => value(model[:p][(index, branch["f_bus"], branch["t_bus"])]))
        push!(branch, "pt" => value(model[:p][(index, branch["t_bus"], branch["f_bus"])]))
    end
    println("The solver termination status is")

    result = Dict{String,Any}(
            "optimizer" => JuMP.solver_name(model),
            "termination_status" => JuMP.termination_status(model),   # "LOCALLY_SOLVED",
            "primal_status" => JuMP.primal_status(model),
            "dual_status" => JuMP.dual_status(model),
            "objective" => InfrastructureModels._guard_objective_value(model),
            "objective_lb" => InfrastructureModels._guard_objective_bound(model),
            "solve_time" => solve_time,
            "solution" => solution)

    return result
end

# json_path = "C:/Users/fmeier/pandapower/pandapower/test/opf/case5_clm_matfile_va.json"
# # #@enter run_powermodels(json_path)
# #
# result = run_powermodels(json_path)
# println(result["termination_status"] == LOCALLY_SOLVED)
# println(isapprox(result["objective"], 17015.5; atol = 1e0))
# mit eingeschränkter slack spannung: 17082.819507648066


text:
#==================

Tests for the LSHFunction() function for constructing hash functions.

==================#

using Test, Random, LSHFunctions

include(joinpath("..", "utils.jl"))

#==================
Tests
==================#

@testset "Test similarity function registration" begin
    Random.seed!(RANDOM_SEED)

    @testset "Register a custom similarity" begin
        ### Test 1: register a function with @register_similarity!
        mysim(x,y) = dot(x,y) / (norm(x) * norm(y))

        @test_throws MethodError LSHFunction(mysim)
        @test_throws MethodError lsh_family(mysim)

        LSHFunctions.@register_similarity!(mysim, SimHash)

        hashfn = LSHFunction(mysim)
        @test isa(hashfn, SimHash)
        @test lsh_family(mysim) == SimHash

        ### Test 2: register a function-like object with @register_similarity!
        struct mytype end

        @test_throws MethodError LSHFunction(mytype())
        @test_throws MethodError lsh_family(mytype())
        @test_throws ErrorException LSHFunctions.@register_similarity!(mytype(), SimHash)

        (::mytype)(x,y) = cossim(x,y)

        LSHFunctions.@register_similarity!(mytype(), SimHash)
        hashfn = LSHFunctions.LSHFunction(mytype())
        @test isa(hashfn, SimHash)
        @test lsh_family(mysim) == SimHash

        ### Test 3: reset the available similarity functions
        LSHFunctions.@reset_similarities!()

        @test_throws MethodError LSHFunction(mysim)
        @test_throws MethodError LSHFunction(mytype())
    end
end

@testset "Test LSHFunction()" begin
    Random.seed!(RANDOM_SEED)

    @testset "Create cosine similarity hash function" begin
        hashfn = LSHFunction(cossim)

        @test similarity(hashfn) == cossim
        @test n_hashes(hashfn) == 1
        @test isa(hashfn, SimHash)
    end

    @testset "Create L^p distance hash function" begin
        hashfn = LSHFunction(ℓ1, 20; scale = 4.0)

        @test similarity(hashfn) == ℓ1
        @test n_hashes(hashfn) == 20
        @test isa(hashfn, LSHFunctions.LpHash)
        @test hashfn.scale == 4.0
    end

    @testset "Create Jaccard similarity hash function" begin
        hashfn = LSHFunction(jaccard, 64)

        @test similarity(hashfn) == jaccard
        @test n_hashes(hashfn) == 64
        @test isa(hashfn, MinHash)
    end

    @testset "Create inner product similarity hash function" begin
        hashfn = LSHFunction(inner_prod, 10; maxnorm=1)

        @test similarity(hashfn) == inner_prod
        @test n_hashes(hashfn) == 10
        @test isa(hashfn, SignALSH)

        # Test that same exceptions are thrown for invalid construction of
        # SignALSH
        @test_throws ErrorException LSHFunction(inner_prod, -1; maxnorm=10)
        @test_throws ErrorException LSHFunction(inner_prod; m=0, maxnorm=10)
        @test_throws ErrorException LSHFunction(inner_prod)
    end

    @testset "Call LSHFunction() with invalid similarity function" begin
        import LinearAlgebra: dot, norm
        my_cossim(x,y) = dot(x,y) / (norm(x) * norm(y))

        @test_throws(MethodError, LSHFunction(my_cossim))
    end
end


text:
# General interface methods for SummationByPartsOperators.jl and Trixi.jl
# TODO: FD. Move to another file
nnodes(D::AbstractDerivativeOperator) = size(D, 1)
eachnode(D::AbstractDerivativeOperator) = Base.OneTo(nnodes(D))
get_nodes(D::AbstractDerivativeOperator) = grid(D)

# For dispatch
# TODO: FD. Move to another file
const FDSBP = DG{Basis} where {Basis<:AbstractDerivativeOperator}

# 2D containers
# TODO: FD. Move to another file
init_mortars(cell_ids, mesh, elements, mortar::Nothing) = nothing

# 2D caches
function create_cache(mesh::TreeMesh{2}, equations,
                      volume_integral::VolumeIntegralStrongForm, dg, uEltype)

  prototype = Array{SVector{nvariables(equations), uEltype}, ndims(mesh)}(
    undef, ntuple(_ -> nnodes(dg), ndims(mesh))...)
  f_threaded = [similar(prototype) for _ in 1:Threads.nthreads()]

  return (; f_threaded,)
end

create_cache(mesh, equations, mortar::Nothing, uEltype) = NamedTuple()
nmortars(mortar::Nothing) = 0

# 2D RHS
function calc_volume_integral!(du, u,
                               mesh::TreeMesh{2},
                               nonconservative_terms::Val{false}, equations,
                               volume_integral::VolumeIntegralStrongForm,
                               dg::FDSBP, cache)
  D = dg.basis # SBP derivative operator
  @unpack f_threaded = cache

  # SBP operators from SummationByPartsOperators.jl implement the basic interface
  # of matrix-vector multiplication. Thus, we pass an "array of structures",
  # packing all variables per node in an `SVector`.
  if nvariables(equations) == 1
    # `reinterpret(reshape, ...)` removes the leading dimension only if more
    # than one variable is used.
    u_vectors  = reshape(reinterpret(SVector{nvariables(equations), eltype(u)}, u),
                         nnodes(dg), nnodes(dg), nelements(dg, cache))
    du_vectors = reshape(reinterpret(SVector{nvariables(equations), eltype(du)}, du),
                         nnodes(dg), nnodes(dg), nelements(dg, cache))
  else
    u_vectors  = reinterpret(reshape, SVector{nvariables(equations), eltype(u)}, u)
    du_vectors = reinterpret(reshape, SVector{nvariables(equations), eltype(du)}, du)
  end

  # Use the tensor product structure to compute the discrete derivatives of
  # the fluxes line-by-line and add them to `du` for each element.
  @threaded for element in eachelement(dg, cache)
    f_element = f_threaded[Threads.threadid()]
    u_element = view(u_vectors,  :, :, element)

    # x direction
    @. f_element = flux(u_element, 1, equations)
    for j in eachnode(dg)
      mul!(view(du_vectors, :, j, element), D, view(f_element, :, j),
           one(eltype(du)), one(eltype(du)))
    end

    # y direction
    @. f_element = flux(u_element, 2, equations)
    for i in eachnode(dg)
      mul!(view(du_vectors, i, :, element), D, view(f_element, i, :),
           one(eltype(du)), one(eltype(du)))
    end
  end

  return nothing
end

function prolong2mortars!(cache, u, mesh, equations, mortar::Nothing,
                          surface_integral, dg::DG)
  @assert isempty(eachmortar(dg, cache))
end

function calc_mortar_flux!(surface_flux_values, mesh,
                           nonconservative_terms, equations,
                           mortar::Nothing,
                           surface_integral, dg::DG, cache)
  @assert isempty(eachmortar(dg, cache))
end

function calc_surface_integral!(du, u, mesh::TreeMesh{2},
                                equations, surface_integral::SurfaceIntegralStrongForm,
                                dg::DG, cache)
  inv_weight_left  = inv(left_boundary_weight(dg.basis))
  inv_weight_right = inv(right_boundary_weight(dg.basis))
  @unpack surface_flux_values = cache.elements

  @threaded for element in eachelement(dg, cache)
    for l in eachnode(dg)
      # surface at -x
      u_node = get_node_vars(u, equations, dg, 1, l, element)
      f_node = flux(u_node, 1, equations)
      f_num  = get_node_vars(surface_flux_values, equations, dg, l, 1, element)
      surface_contribution = -(f_num - f_node) * inv_weight_left
      add_to_node_vars!(du, surface_contribution, equations, dg, 1, l, element)

      # surface at +x
      u_node = get_node_vars(u, equations, dg, nnodes(dg), l, element)
      f_node = flux(u_node, 1, equations)
      f_num  = get_node_vars(surface_flux_values, equations, dg, l, 2, element)
      surface_contribution = +(f_num - f_node) * inv_weight_right
      add_to_node_vars!(du, surface_contribution, equations, dg, nnodes(dg), l, element)

      # surface at -y
      u_node = get_node_vars(u, equations, dg, l, 1, element)
      f_node = flux(u_node, 2, equations)
      f_num  = get_node_vars(surface_flux_values, equations, dg, l, 3, element)
      surface_contribution = -(f_num - f_node) * inv_weight_left
      add_to_node_vars!(du, surface_contribution, equations, dg, l, 1, element)

      # surface at +y
      u_node = get_node_vars(u, equations, dg, l, nnodes(dg), element)
      f_node = flux(u_node, 2, equations)
      f_num  = get_node_vars(surface_flux_values, equations, dg, l, 4, element)
      surface_contribution = +(f_num - f_node) * inv_weight_right
      add_to_node_vars!(du, surface_contribution, equations, dg, l, nnodes(dg), element)
    end
  end

  return nothing
end

# AnalysisCallback
# TODO: FD. Move to another file
SolutionAnalyzer(D::AbstractDerivativeOperator) = D

function integrate_via_indices(func::Func, u,
                               mesh::TreeMesh{2}, equations,
                               dg::FDSBP, cache, args...; normalize=true) where {Func}
  # TODO: FD. This is rather inefficient right now and allocates...
  weights = diag(SummationByPartsOperators.mass_matrix(dg.basis))

  # Initialize integral with zeros of the right shape
  integral = zero(func(u, 1, 1, 1, equations, dg, args...))

  # Use quadrature to numerically integrate over entire domain
  for element in eachelement(dg, cache)
    volume_jacobian_ = volume_jacobian(element, mesh, cache)
    for j in eachnode(dg), i in eachnode(dg)
      integral += volume_jacobian_ * weights[i] * weights[j] * func(u, i, j, element, equations, dg, args...)
    end
  end

  # Normalize with total volume
  if normalize
    integral = integral / total_volume(mesh)
  end

  return integral
end

function calc_error_norms(func, u, t, analyzer,
                          mesh::TreeMesh{2}, equations, initial_condition,
                          dg::FDSBP, cache, cache_analysis)
  # TODO: FD. This is rather inefficient right now and allocates...
  weights = diag(SummationByPartsOperators.mass_matrix(dg.basis))
  @unpack node_coordinates = cache.elements

  # Set up data structures
  l2_error   = zero(func(get_node_vars(u, equations, dg, 1, 1, 1), equations))
  linf_error = copy(l2_error)

  # Iterate over all elements for error calculations
  for element in eachelement(dg, cache)
    # Calculate errors at each node
    volume_jacobian_ = volume_jacobian(element, mesh, cache)

    for j in eachnode(analyzer), i in eachnode(analyzer)
      u_exact = initial_condition(
        get_node_coords(node_coordinates, equations, dg, i, j, element), t, equations)
      diff = func(u_exact, equations) - func(
        get_node_vars(u, equations, dg, i, j, element), equations)
      l2_error += diff.^2 * (weights[i] * weights[j] * volume_jacobian_)
      linf_error = @. max(linf_error, abs(diff))
    end
  end

  # For L2 error, divide by total volume
  total_volume_ = total_volume(mesh)
  l2_error = @. sqrt(l2_error / total_volume_)

  return l2_error, linf_error
end

# TODO: FD. Visualization
#       We need to define better interfaces for all the plotting stuff.
#       Right now, the easiest solution is to use scatter plots such as
#       x = semi.cache.elements.node_coordinates[1, :, :, :] |> vec
#       y = semi.cache.elements.node_coordinates[2, :, :, :] |> vec
#       scatter(x, y, sol.u[end])


text:
using CIM
using Test

@test CIM.dummy_func(2,1) == 5
@test CIM.dummy_func(3,1) == 7


text:
"""
    spliiter(n, k)

Internal utility function, splits 1:n sequence to k chunks of approximately same size.
"""
function splitter(n, k)
    xz = Int.(ceil.(range(0, n, length = k+1)))
    return [xz[i]+1:xz[i+1] for i in 1:k]
end

"""
    colwise!(target, x, y, n_threads)

Internal function for colwise calculations. Let `x` is a matrix `m x n` and `y` is a vector of the length `m`.
Then the `colwise!` function computes distance between each column in `x` and `y` and store result
in `target` array. Argument `n_threads` defines the number of threads used for calculation.
"""
function colwise!(target, x, y, n_threads = Threads.nthreads())
    ncol = size(x, 2)

    # we could have used same algorithm for n_threads == 1
    # but ran into unnecessary allocations and split calculations.
    # Impact is neglible, yet there is no need to do extra calculations
    if n_threads != 1
        ranges = splitter(ncol, n_threads)
        waiting_list = Task[]

        for i in 1:length(ranges) - 1
            push!(waiting_list, @spawn chunk_colwise!(target, x, y, ranges[i]))
        end
        chunk_colwise!(target, x, y, ranges[end])

        for i in 1:length(ranges) - 1
            wait(waiting_list[i])
        end
    else
        chunk_colwise!(target, x, y, axes(x, 2))
    end

    target
end

"""
    chunk_colwise!(target, x, y, r)

Utility function for calculation of the `colwise!(target, x, y, n_threads)` function.
UnitRange argument `r` select subarray of original design matrix `x` that is going
to be processed.
"""
function chunk_colwise!(target, x, y, r)
    @inbounds for j in r
        res = 0.0
        for i in axes(x, 1)
            res += (x[i, j] - y[i])^2
        end
        target[j] = res
    end
end

"""
    smart_init(X, k; init="k-means++")

This function handles the random initialisation of the centroids from the
design matrix (X) and desired groups (k) that a user supplies.

`k-means++` algorithm is used by default with the normal random selection
of centroids from X used if any other string is attempted.

A named tuple representing centroids and indices respecitively is returned.
"""
function smart_init(X, k, n_threads = Threads.nthreads();
        init = "k-means++")

    n_row, n_col = size(X)

    if init == "k-means++"

        # randonmly select the first centroid from the data (X)

        # TODO relax constraints on distances, may be should
        # define `X` as X::AbstractArray{T} where {T <: Number}
        # and use this T for all calculations.
        centroids = zeros(n_row, k)
        rand_indices = Vector{Int}(undef, k)
        rand_idx = rand(1:n_col)
        rand_indices[1] = rand_idx
        centroids[:, 1] .= X[:, rand_idx]
        distances = Vector{Float64}(undef, n_col)
        new_distances = Vector{Float64}(undef, n_col)

        # compute distances from the first centroid chosen to all the other data points
        colwise!(distances, X, centroids[:, 1], n_threads)
        distances[rand_idx] = 0.0

        for i = 2:k
            # choose the next centroid, the probability for each data point to be chosen
            # is directly proportional to its squared distance from the nearest centroid
            r_idx = wsample(1:n_col, vec(distances))
            rand_indices[i] = r_idx
            centroids[:, i] .= X[:, r_idx]

            # no need for final distance update
            i == k && break

            # compute distances from the centroids to all data points
            colwise!(new_distances, X, centroids[:, i], n_threads)

            # and update the squared distance as the minimum distance to all centroid
            for i in 1:n_col
                distances[i] = distances[i] < new_distances[i] ? distances[i] : new_distances[i]
            end
            distances[r_idx] = 0.0
        end

    else
        # randomly select points from the design matrix as the initial centroids
        rand_indices = sample(1:n_col, k, replace = false)
        centroids = X[:, rand_indices]
    end

    return (centroids = centroids, indices = rand_indices)
end


text:
@scenario Simple = begin
    q₁::Float64
    q₂::Float64
    d₁::Float64
    d₂::Float64
end

s₁ = SimpleScenario(-24.0, -28.0, 500.0, 100.0, probability = 0.4)
s₂ = SimpleScenario(-28.0, -32.0, 300.0, 300.0, probability = 0.6)

simple = StochasticProgram([s₁,s₂])

@first_stage simple = begin
    @variable(model, x₁ >= 40)
    @variable(model, x₂ >= 20)
    @objective(model, Min, 100*x₁ + 150*x₂)
    @constraint(model, x₁ + x₂ <= 120)
end

@second_stage simple = begin
    @decision x₁ x₂
    ξ = scenario
    @variable(model, 0 <= y₁ <= ξ.d₁)
    @variable(model, 0 <= y₂ <= ξ.d₂)
    @objective(model, Min, ξ.q₁*y₁ + ξ.q₂*y₂)
    @constraint(model, 6*y₁ + 10*y₂ <= 60*x₁)
    @constraint(model, 8*y₁ + 5*y₂ <= 80*x₂)
end

push!(problems,(simple,"Simple"))


text:
using ReachabilityAnalysis
Dict([
    "opts" => Dict(:T=>7.0, :alg=>TMJets(max_steps=1_000, abs_tol=1e-10, orderT=6, orderQ=2)),
    "X0"   => Hyperrectangle(low=[1.25, 2.25, 1.25, 2.25], high=[1.55, 2.35, 1.55, 2.35]),
    "info" => Dict(
        "name" => "coupled_vanderpol",
        "dim" => 4,
        "linear" => false,
        "hybrid" => false
    )
])


text:
module GenerativeMetrics
    using Random, Statistics, LinearAlgebra

    include("metrics.jl")
end


text:
using ADCME
using PyCall
using LinearAlgebra
using PyPlot
using Random
using DelimitedFiles
Random.seed!(233)

function mpi_tensor_transpose(row,col,ncol,val,n,rank,nt)
    require_mpi()
    mpi_tensor_transpose_ = load_op_and_grad("./build/libMPITensor.so","mpi_tensor_transpose", multiple=true)
    row,col,ncol,val,n,rank,nt = convert_to_tensor(Any[row,col,ncol,val,n,rank,nt], [Int32,Int32,Int32,Float64,Int64,Int64,Int64])
    indices, vals = mpi_tensor_transpose_(row,col,ncol,val,n,rank,nt)
end

mpi_init()

if mpi_rank()==0
    row = [0;1]
    col = [0;1;2;3;0;1;2;3]
    ncol = [4;4]
    val = Float64[1;2;3;4;5;6;7;8]
else
    row = [2,3]
    col = [0;1;2;3;0;1;2;3]
    ncol = [4;4]
    val = Float64[1;2;3;4;5;6;7;8] .+ 8.0
end

n = 2
rank = mpi_rank()
nt = 4

# TODO: specify your input parameters
u = mpi_tensor_transpose(row,col,ncol,val,n,rank,nt)
sess = Session(); init(sess)

indices, vals = run(sess, u)
sleep(mpi_rank())
println("===================================================")
writedlm(stdout, indices)
writedlm(stdout, vals)

mpi_finalize()


text:
using FunctionWrappers; Fn = FunctionWrappers.FunctionWrapper
let
    #= line 1 =#
    var"|main|b_implicitStr"::String = "114514"::String
    var"|main|c_implicitInt"::Function = (x->begin
                    #= /home/tainez/github/MLFS/test/testim.jl:26 =#
                    parse(Int, x)
                end)::Function
    var"|main|d_+"::Function = (((x, y),)->begin
                    #= /home/tainez/github/MLFS/test/testim.jl:28 =#
                    x + y
                end)::Function
    var"|main|e_fun"::Function = function (var"|main|f_x",)
                #= /home/tainez/github/MLFS/test/testim.jl:30 =#
                function (var"|main|g_y",)
                        #= /home/tainez/github/MLFS/test/testim.jl:30 =#
                        (var"|main|d_+"::Function)((var"|main|f_x"::Any, var"|main|g_y"::Any)::Tuple{Any, Any})::Any
                    end::Function
            end::Function
    var"|main|h_var"::Any = (var"|main|e_fun"(var"|main|c_implicitInt"(var"|main|b_implicitStr"::String)::Any)::Function)(1::Any)::Any
    0
end


text:
module DaemonMode

export serve, runfile, runargs, runexpr, sendExitCode

using Sockets

const PORT = 3000

function add_packages(fname::AbstractString)
end

const first_time = Ref(true)
const token_runfile = "DaemonMode::runfile"
const token_runexpr = "DaemonMode::runexpr"
const token_exit = "DaemonMode::exit"
const token_end = "DaemonMode::end"

"""
    serve(port=3000, shared=false)

Run the daemon, running all files and expressions sended by the client function.

# Optionals

- port: port to listen (default=3000).
- shared: Share the environment between calls. If it is false (default) each run
  has its own environment, so the variables/functions are not shared.
"""
function serve(port=PORT, shared=missing)
    server = Sockets.listen(Sockets.localhost, port)
    quit = false
    current = pwd()

    while !quit
        sock = accept(server)
        mode = readline(sock)

        redirect_stdout(sock) do
            redirect_stderr(sock) do

                if mode == token_runfile
                    shared_client = shared

                    if ismissing(shared_client)
                        shared_client = false
                    end
                    serverRunFile(sock, shared_client)
                elseif mode == token_runexpr
                    shared_client = shared

                    if ismissing(shared_client)
                        shared_client = true
                    end
                    serverRunExpr(sock, shared_client)
                elseif mode == token_exit
                    println(sock, token_end)
                    sleep(1)
                    quit = true
                else
                    println(sock, "Error, unrecognised mode, expected (\"runFile()\", \"runExpr()\" or \"exit()\", but received \"$mode\"")
                    quit = true
                end

            end
        end

    end
    close(server)
end
"""
serverRunFile(sock, shared)

Run the source code of the filename push through the socket.

# Parameters

- sock: socket in which is going to receive the dir, the filename, and args to run.
- shared: Share the environment between calls. If it is false (default) each run
has its own environment, so the variables/functions are not shared.
"""
function serverRunFile(sock, shared)
    dir = readline(sock)
    fname = readline(sock)
    args_str = readline(sock)
    args = split(args_str, " ")

    append!(ARGS, args)

    first_time[] = true
    error = ""

    try
        cd(dir) do
            content = join(readlines(fname), "\n")

            if (!shared)
                m = Module()
                include_string(m, content)
            else
                include(fname)
            end
        end
    catch e
        if isa(e, LoadError)
            if :msg in propertynames(e.error)
                error_msg = e.error.msg
            else
                error_msg = "$(e.error)"
            end
            error = "ERROR in line $(e.line): '$(error_msg)'"
        else
            error = "ERROR: could not open file '$fname'"
        end
    end

    try
        !isempty(error) && println(sock, error)
        println(sock, token_end)
    catch e
        if (e isa Base.IOError) && abs(e.code) == abs(Libc.EPIPE)
            # client disconnected early, ignore
        else
            rethrow()
        end
    end
    empty!(ARGS)
end

"""
serverRunExpr(sock, shared)

Run the source code of the filename push through the socket.

# Parameters

- sock: socket in which is going to receive the code to run.
- shared: Share the environment between calls. If it is false (default) each run
has its own environment, so the variables/functions are not shared.
"""
function serverRunExpr(sock, shared)
    dir = readline(sock)
    expr = readuntil(sock, token_end) # Read until token_end to handle multi-line expressions
    error = ""

    try
        cd(dir) do
            if shared
                evaledExpr = Meta.parse(expr)
                Main.eval(evaledExpr)
            else
                evaledExpr = Meta.parse(expr)
                m = Module()
                Base.eval(m, evaledExpr)
            end
        end
    catch e
        if isa(e, LoadError)
            if :msg in propertynames(e.error)
                error_msg = e.error.msg
            else
                error_msg = "$(e.error)"
            end
            error = "ERROR in line $(e.line): '$(error_msg)'"
        else
            error = "ERROR: $e"
        end
    end

    try
        !isempty(error) && println(sock, error)
        println(sock, token_end)
    catch e
        if (e isa Base.IOError) && abs(e.code) == abs(Libc.EPIPE)
            # client disconnected early, ignore
        else
            rethrow()
        end
    end

    empty!(ARGS)
end

"""
    runexpr(expr::AbstractString)

Ask the server to run julia code in a string pass as parameters.

# Parameters

- expr: Julia code to run in the server.

# Optionals

- port: Port (default=3000).
- output: stream in which it is shown the output of the run.
"""
function runexpr(expr::AbstractString ; output = stdout, port = PORT)
    try
        sock = Sockets.connect(port)
        println(sock, token_runexpr)
        println(sock, pwd())
        println(sock, expr)
        println(sock, token_end)

        line = readline(sock)
        while (line != token_end)
            println(output, line)
            line = readline(sock)
        end

    catch e
        println(stderr, "Error, cannot connect with server. Is it running?")
    end
end

"""
    runfile(fname::AbstractString)

Ask the server to run a specific filename.

# Parameters

- fname: Filename to run.

# Optionals

- args: List of arguments (array of String, default=[]).
- port: Port (default=3000)
- output: stream in which it is shown the output of the run.
"""
function runfile(fname::AbstractString; args=String[], port = PORT, output=stdout)
    dir = dirname(fname)

    if isempty(dir)
        fcompletename = joinpath(dir, fname)
    else
        fcompletename = fname
    end
    try
        sock = Sockets.connect(port)
        println(sock, token_runfile)
        println(sock, pwd())
        println(sock, fcompletename)
        println(sock, join(args, " "))
        line = readline(sock)
        while (line != token_end)
            println(output, line)
            line = readline(sock)
        end
    catch e
        println(stderr, "Error, cannot connect with server. Is it running?")
    end
    return
end

"""
    sendExitCode(port)

send the exit code, it closes the server.

# Optionals

- port: port to connect (default=3000).
"""
function sendExitCode(port=PORT)
    sock = Sockets.connect(port)
    println(sock, token_exit)
end

"""
    runargs(port=PORT)

Ask the server to run all files in ARGS.

# Optionals

- port: Port of the server to ask (default=3000).
"""
function runargs(port=PORT)
    if isempty(ARGS)
        println(file=stderr, "Error: missing filename")
    end
    runfile(ARGS[1], args=ARGS[2:end], port=port)
end

end # module


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

# types and constructors

# We need a global variable to keep track of gnuplot's state
type GnuplotState
    running::Bool               # true when gnuplot is already running
    current::Int                # current figure
    fid                         # pipe stream id
    tmpdir::String              # where to store data files
    figs::Array                 # storage for all figures

    function GnuplotState(running::Bool,current::Int,fid,tmpdir::String,
        figs::Array)
        # Check to see if tmpdir exists, and create it if not
        # TODO: there has to be a simpler way to do this
        try
            f = open(tmpdir)
            close(f)
        catch
            run(`mkdir $tmpdir`)
        end
        new(running,current,fid,tmpdir,figs)
    end
end

# Structure to keep Gaston's configuration
type GastonConfig
    # default CurveConf values
    legend::String
    plotstyle::String
    color::String
    marker::String
    linewidth::Real
    pointsize::Real
    # default AxesConf values
    title::String
    xlabel::String
    ylabel::String
    zlabel::String
    box::String
    axis::String
    # default grid type
    grid::String
    # default terminal type
    terminal::String
    # for terminals that support filenames
    outputfile::String
    # for printing to file
    print_color::String
    print_fontface::String
    print_fontsize::Real
    print_fontscale::Real
    print_linewidth::Real
    print_size::String
end
GastonConfig() = GastonConfig(
    # CurveConf
    "","lines","","",1,0.5,
    # AxesConf
    "Untitled","x","y","z","inside vertical right top","",
    # grid
    "off",
    # terminal
    "wxt",
    # output file name
    "",
    # print parameters
    "color", "Sans", 12, 0.5, 1, "5in,3in"
    )

# Structs to configure a plot
# Two types of configuration are needed: one to configure a single curve, and
# another to configure a set of curves (the 'axes').
type CurveConf
    legend::String          # legend text
    plotstyle::String
    color::String           # one of gnuplot's builtin colors --
                            # run 'show colornames' inside gnuplot
    marker::String          # point type
    linewidth::Real
    pointsize::Real

end
CurveConf() = CurveConf(
    gaston_config.legend,
    gaston_config.plotstyle,
    gaston_config.color,
    gaston_config.marker,
    gaston_config.linewidth,
    gaston_config.pointsize)

type AxesConf
    title::String      # plot title
    xlabel::String     # xlabel
    ylabel::String     # ylabel
    zlabel::String     # zlabel for 3-d plots
    box::String        # legend box (used with 'set key')
    axis::String       # normal, semilog{x,y}, loglog
end
AxesConf() = AxesConf(
    gaston_config.title,
    gaston_config.xlabel,
    gaston_config.ylabel,
    gaston_config.zlabel,
    gaston_config.box,
    gaston_config.axis)

# coordinates and configuration for a single curve
type CurveData
    x::Vector          # abscissa
    y::Vector          # ordinate
    Z::Array           # 3-d plots and images
    ylow::Vector       # error data
    yhigh::Vector      # error data
    conf::CurveConf    # curve configuration
end
CurveData() = CurveData([],[],[],[],[],CurveConf())
CurveData(x,y,Z,conf::CurveConf) = CurveData(x,y,Z,[],[],conf)

# curves and configuration for a single figure
type Figure
    handle::Int                  # each figure has a unique handle
    curves::Vector{CurveData}    # a vector of curves (x,y,conf)
    conf::AxesConf               # figure configuration
    isempty::Bool                # a flag to indicate if figure is empty
end
Figure(handle) = Figure(handle,[CurveData()],AxesConf(),true)

# coordinate type
Coord = Union(Range1,Range,Matrix,Vector)


text:
module TestReduce

include("preamble.jl")

@testset begin
    n = 3
    C = zeros(n)
    A = rand(-10:10, n, n)
    B = rand(-10:10, n)
    CAB = ndreducible(
        instantiate(broadcasted(Ref, Ref(C), eachindex(C))) => (:i,),
        A => (:i, :j),
        B => (:j,),
    )
    @test CAB.plan.indices == Index.((:j, :i))
    foreach(CAB) do (c, a, b)
        c[] += a * b
    end
    @test C == A * B
end

end  # module


text:
"""
    sun_generators(n)

Returns an array of sparse ``n \\times n`` generators of ``SU(N)``
that are Hermitian and trless.
"""
function sun_generators(n)
    matrices = [spzeros(ComplexF64, n, n) for i = 1 : n^2-1]
    i = 1

    # Diagonal matrices
    for l = 1:n-1
        for j = 1:l
            matrices[i][j, j] = sqrt(2/(l*(l+1)))
        end
        matrices[i][l+1, l+1] = -sqrt(2*l/(l+1))
        i += 1
    end

    # Real matrices
    for k = 2:n
        for j = 1:k-1
            matrices[i][j, k] = matrices[i][k, j] = 1.
            i += 1
        end
    end

    # Imaginary matrices
    for k = 2:n
        for j = 1:k-1
            matrices[i][j, k] -= (matrices[i][k, j] = 1im)
            i += 1
        end
    end

    return matrices
end

"""
    bloch_vector(ρ)

Returns the generalized Bloch vector ``\\mathbf{r}`` for the density operator ρ.

If ρ is a ``N × N`` density operator, ``\\mathbf{r}`` will be a ``N^2 - 1``
vector of real values. We follow the formula
``ρ = \\frac{1}{N} (I + \\sqrt{N} \\mathbf{r} \\cdot \\mathbf{λ})``
"""
function bloch_vector(ρ)
   n = size(ρ)[1]
   λ = sun_generators(n)
   sqrt(n)/2*real([tr(λ[i] * ρ) for i = 1: length(λ)])
end

"""
    density_operator(r)

Returns the density operator ρ corresponding to the Bloch vector `r`.

`vec` must be a vector of ``N^2-1`` real values. We follow the formula
``ρ = \\frac{1}{N} (I + \\sqrt{N} \\mathbf{r} \\cdot \\mathbf{λ})``
"""
function density_operator(vec::Array{Float64})
   n = length(vec)
   Nq = Int(floor(sqrt(n+1))) # Dimension of the corresponding Hilbert space
   λ = sun_generators(Nq) # Obtain the generators of the SU(N) Hilbert space

   1/Nq*(I + sqrt(Nq) * sum([vec[i] * λ[i] for i = 1:n]))
end

# We return im * gen so that the matrix is real (more efficient)

"""
    igen(H)

Returns ``i G(H)``, where G is defined in Eq. (11) of Joynt et al.
"""
function igen(H)
    n = size(H,1)
    λ = sun_generators(n)
    result = spzeros(n^2-1, n^2-1)
    for i = 1:n^2-1
        for j = i+1:n^2-1
            tmp = imag(0.5 * tr(commutator(λ[i],λ[j])*H))
            if tmp != 0
                @inbounds result[i,j] = tmp
                @inbounds result[j,i] = -tmp
            end
        end
    end
    return result
end

"""
    igen(f, λ, H)

Returns ``i G(H)``, where G is defined in Eq. (11) of Joynt et al., with the
structure constants tensor `f` that can be evaluated with `structure_constants(λ)`,
and the array of SU(N) generators λ (`sun_generators(n)`).

When providing `f`, the calculations are much faster (especially if many
generators have to be evaluated).
"""
function igen(f, λ, H)
    n = size(H,1)
    result = spzeros(n^2-1, n^2-1)

    if eltype(H) <: Real
        # If the Hamiltonian is real, the imaginary generators do not contribute
        nimag = Int(floor(n*(n-1)/2))
        ak = zeros(n^2-1)
        ak[1:end-nimag] = real([tr(H*lambda) for lambda in λ[1:end-nimag]])
    else
        ak = real([tr(H*lambda) for lambda in λ])
    end

    for i = 1:n^2-1
        @inbounds result[i, :] .= f[i] * ak
    end
    return result
end

"""
    structure_constants(λ)

Returns the totally antisymmetric tensor ``i f_{ijk}`` of the structure constants
of the list of generators λ, defined as ``f_{ijk} = - 2 i \\text{Tr}([λ_i,λ_j]λ_k)``.

Notice that the true definition has a factor ``-i``, which we avoid to store only
a real value.
"""
function structure_constants(lambda)::Array{SparseArrays.SparseMatrixCSC{Float64,Int64},3}
    Nl = length(lambda)
    f = cat([spzeros(Float64, Nl, Nl) for i = 1:Nl], dims=3)
    for i = 1:Nl
        for j = i+1:Nl
            for k = j+1:Nl
                tmp::Float64 = 1. / 4. * imag(tr(commutator(lambda[i],lambda[j])*lambda[k]))
                if tmp != 0
                    @inbounds begin
                        f[i][j,k] +=  tmp
                        f[i][k,j] -=  tmp
                        f[k][i,j] +=  tmp
                        f[j][i,k] -=  tmp
                        f[j][k,i] +=  tmp
                        f[k][j,i] -=  tmp
                    end
                end
            end
        end
    end
    return f
end


text:
#!/usr/local/bin/julia

using ArgParse

function translate(RNAString)
	CodonTable = Dict("UUU" => 'F', "CUU" => 'L', "AUU" => 'I', "GUU" => 'V',
					"UUC" => 'F', "CUC" => 'L', "AUC" => 'I', "GUC" => 'V',
					"UUA" => 'L', "CUA" => 'L', "AUA" => 'I', "GUA" => 'V',
					"UUG" => 'L', "CUG" => 'L', "AUG" => 'M', "GUG" => 'V',
					"UCU" => 'S', "CCU" => 'P', "ACU" => 'T', "GCU" => 'A',
					"UCC" => 'S', "CCC" => 'P', "ACC" => 'T', "GCC" => 'A',
					"UCA" => 'S', "CCA" => 'P', "ACA" => 'T', "GCA" => 'A',
					"UCG" => 'S', "CCG" => 'P', "ACG" => 'T', "GCG" => 'A',
					"UAU" => 'Y', "CAU" => 'H', "AAU" => 'N', "GAU" => 'D',
					"UAC" => 'Y', "CAC" => 'H', "AAC" => 'N', "GAC" => 'D',
					"UAA" => ' ', "CAA" => 'Q', "AAA" => 'K', "GAA" => 'E',
					"UAG" => ' ', "CAG" => 'Q', "AAG" => 'K', "GAG" => 'E',
					"UGU" => 'C', "CGU" => 'R', "AGU" => 'S', "GGU" => 'G',
					"UGC" => 'C', "CGC" => 'R', "AGC" => 'S', "GGC" => 'G',
					"UGA" => ' ', "CGA" => 'R', "AGA" => 'R', "GGA" => 'G',
					"UGG" => 'W', "CGG" => 'R', "AGG" => 'R', "GGG" => 'G')

	Codons = [RNAString[idx:idx+2] for idx in 1:3:length(RNAString)]
	Protein = strip(join([CodonTable[Codon] for Codon in Codons]))

end

function main(args)
	# Argument parsing
	configs = ArgParseSettings()

	add_arg_group!(configs, "Options:")
	@add_arg_table! configs begin
		"--input", "-i"
			help	 = "Input filename"
			action	 = :store_arg
			arg_type = String
		"--output", "-o"
			help	 = "Output filename"
			action	 = :store_arg
			arg_type = String
		end
	# Store arguments to variables
	parsed_args = parse_args(args, configs)
	infile      = parsed_args["input"]
	outfile     = parsed_args["output"]

	# Read file and compute problem via function
	RNAString = chomp(open(f->read(f, String), infile))
	@show ProteinString = translate(RNAString)
	open(f->write(f, "$ProteinString"), outfile, "w+")

end

main(ARGS)


text:
immutable AffineTransform{T,N}
    scalefwd::Matrix{T}
    scaleinv::Matrix{T}
    offset::Vector{T}
    temp::Vector{T}

    function AffineTransform(S::AbstractMatrix, t::AbstractVector)
        if size(S) != (N,N) || size(t) != (N,)
            throw(DimensionMismatch("Inputs must be $N-dimensional"))
        end
        temp = Array(eltype(t), length(t))
        new(S, inv(S), t, temp)
    end
end
function AffineTransform(S, offset)
    T = promote_type(eltype(S),eltype(offset))
    AffineTransform{T,length(offset)}(convert(Matrix{T},S), convert(Vector{T},offset))
end

ndims{T,N}(tform::AffineTransform{T,N}) = N
eltype{T}(tf::AffineTransform{T}) = T

function show(io::IO, tf::AffineTransform)
    println(io, typeof(tf), ":")
    println(io, "matrix: ", tf.scalefwd)
    println(io, "translation: ", tf.offset)
end

scale(tf::AffineTransform, s::Vector) = AffineTransform(scale(tf.scalefwd, s), tf.offset)
scale(s::Vector, tf::AffineTransform) = AffineTransform(scale(s, tf.scalefwd), tf.offset.*s)

*(a::AffineTransform, v::AbstractVecOrMat) = tformfwd(a, v)
\(a::AffineTransform, x::AbstractVecOrMat) = tforminv(a, x)

tformfwd{T}(a::AffineTransform{T}, x::AbstractVector) = tformfwd!(Array(typeof(one(T)*one(eltype(x))), length(x)), a, x)
tforminv{T}(a::AffineTransform{T}, x::AbstractVector) = tforminv!(Array(typeof(one(T)*one(eltype(x))), length(x)), a, x)
tformfwd{T}(a::AffineTransform{T}, X::AbstractMatrix) = tformfwd!(Array(typeof(one(T)*one(eltype(X))), size(X,1), size(X,2)), a, X)
tforminv{T}(a::AffineTransform{T}, X::AbstractMatrix) = tforminv!(Array(typeof(one(T)*one(eltype(X))), size(X,1), size(X,2)), a, X)

tformfwd(a::AffineTransform, v::AbstractMatrix) = a.scalefwd*v .+ a.offset
tforminv(a::AffineTransform, x::AbstractMatrix) = a.scaleinv*(x .- a.offset)

# Implementation of a.scalefwd*X .+ a.offset
function tformfwd!(dest, a::AffineTransform, X::AbstractVecOrMat)
    A_mul_B!(dest, a.scalefwd, X)
    @inbounds for j = 1:size(X,2), i = 1:size(X,1)
        dest[i,j] += a.offset[i]
    end
    dest
end

# Implementation of a.scaleinv*(X .- a.offset)
function tforminv!(dest, a::AffineTransform, X::AbstractVecOrMat)
    size(X,1) == ndims(a) || throw(DimensionMismatch("Dimensionality $(size(X,1)) of inputs does not match $(ndims(a)) for the transformation"))
    A_mul_B!(dest, a.scaleinv, X)
    A_mul_B!(a.temp, a.scaleinv, a.offset)
    @inbounds for j = 1:size(X,2), i = 1:size(X,1)
        dest[i,j] -= a.temp[i]
    end
    dest
end

# Faster versions for particular dimensionalities (with vectors)
for VecType in (AbstractVector, CartesianIndex)
    @eval begin
        function tforminv!{T}(v, a::AffineTransform{T,3}, x::$VecType)
            length(x) == 3 || throw(DimensionMismatch("vector must be 3-dimensional"))
            t = a.temp
            o = a.offset
            @inbounds t[1] = x[1]-o[1]
            @inbounds t[2] = x[2]-o[2]
            @inbounds t[3] = x[3]-o[3]
            gemv3!(v, a.scaleinv, t)
        end
        function tformfwd!{T}(x, a::AffineTransform{T,3}, v::$VecType)
            length(x) == 3 || throw(DimensionMismatch("vector must be 3-dimensional"))
            gemv3!(x, a.scalefwd, v)
            o = a.offset
            @inbounds x[1] += o[1]
            @inbounds x[2] += o[2]
            @inbounds x[3] += o[3]
            x
        end
    end
end

function tforminv{T}(a::AffineTransform{T,3}, x1::Number, x2::Number, x3::Number)
    o = a.offset
    @inbounds t1 = x1 - o[1]
    @inbounds t2 = x2 - o[2]
    @inbounds t3 = x3 - o[3]
    gemv3(a.scaleinv, t1, t2, t3)
end

function tformfwd{T}(a::AffineTransform{T,3}, x1::Number, x2::Number, x3::Number)
    t1, t2, t3 = gemv3(a.scalefwd, x1, x2, x3)
    o = a.offset
    @inbounds r1 = t1+o[1]
    @inbounds r2 = t2+o[2]
    @inbounds r3 = t3+o[3]
    r1, r2, r3
end

for VecType in (AbstractVector, CartesianIndex)
    @eval begin
        function tforminv!{T}(v, a::AffineTransform{T,2}, x::$VecType)
            length(x) == 2 || throw(DimensionMismatch("vector must be 2-dimensional"))
            t = a.temp
            o = a.offset
            @inbounds t[1] = x[1]-o[1]
            @inbounds t[2] = x[2]-o[2]
            gemv2!(v, a.scaleinv, t)
        end
        function tformfwd!{T}(x, a::AffineTransform{T,2}, v::$VecType)
            length(x) == 2 || throw(DimensionMismatch("vector must be 2-dimensional"))
            gemv2!(x, a.scalefwd, v)
            t = a.offset
            @inbounds x[1] += t[1]
            @inbounds x[2] += t[2]
            x
        end
    end
end

function tforminv{T}(a::AffineTransform{T,2}, x1::Number, x2::Number)
    o = a.offset
    @inbounds t1 = x1 - o[1]
    @inbounds t2 = x2 - o[2]
    gemv2(a.scaleinv, t1, t2)
end

function tformfwd{T}(a::AffineTransform{T,2}, x1::Number, x2::Number)
    t1, t2 = gemv2(a.scalefwd, x1, x2)
    o = a.offset
    @inbounds r1 = t1+o[1]
    @inbounds r2 = t2+o[2]
    r1, r2
end

# For just the 3rd coordinate (we sometimes only need that one)
function tforminv3rd{T}(a::AffineTransform{T,3}, x::AbstractVector)
    o = a.offset
    A = a.scaleinv
    A[3]*(x[1]-o[1]) + A[6]*(x[2]-o[2]) + A[9]*(x[3]-o[3])
end

function scaled_add3!(x, s::Number, r)
    @inbounds x[1]+=s*r[1]
    @inbounds x[2]+=s*r[2]
    @inbounds x[3]+=s*r[3]
end

function gemv3!(res, A, x)
    @inbounds res[1] = A[1]*x[1] + A[4]*x[2] + A[7]*x[3]
    @inbounds res[2] = A[2]*x[1] + A[5]*x[2] + A[8]*x[3]
    @inbounds res[3] = A[3]*x[1] + A[6]*x[2] + A[9]*x[3]
    res
end

function gemv3(A, x1::Number, x2::Number, x3::Number)
    @inbounds r1 = A[1]*x1 + A[4]*x2 + A[7]*x3
    @inbounds r2 = A[2]*x1 + A[5]*x2 + A[8]*x3
    @inbounds r3 = A[3]*x1 + A[6]*x2 + A[9]*x3
    r1, r2, r3
end

function scaled_add2!(x, s::Number, r)
    @inbounds x[1]+=s*r[1]
    @inbounds x[2]+=s*r[2]
end

function gemv2!(res, A, x)
    @inbounds res[1] = A[1]*x[1] + A[3]*x[2]
    @inbounds res[2] = A[2]*x[1] + A[4]*x[2]
    res
end

function gemv2(A, x1::Number, x2::Number)
    @inbounds r1 = A[1]*x1 + A[3]*x2
    @inbounds r2 = A[2]*x1 + A[4]*x2
    r1, r2
end

### Identity transformation
tformeye(T, nd) = AffineTransform(eye(T, nd), zeros(T, nd))
tformeye(nd) = AffineTransform(eye(nd), zeros(nd))

### Rotations
# Right-hand rotation around an axis (axis points towards viewer)

rotation2(angle) = [cos(angle) -sin(angle); sin(angle) cos(angle)]

function tformrotate(angle)
    A = rotation2(angle)
    AffineTransform{eltype(A), 2}(A, zeros(eltype(A),2))
end

# The following assumes uaxis is normalized
function _rotation3(uaxis::Vector, angle)
    if length(uaxis) != 3
        error("3d rotations only")
    end
    ux, uy, uz = uaxis[1], uaxis[2], uaxis[3]
    c = cos(angle)
    s = sin(angle)
    cm = one(typeof(c)) - c
    R = [c+ux*ux*cm       ux*uy*cm-uz*s      ux*uz*cm+uy*s;
         uy*ux*cm+uz*s      c+uy*uy*cm       uy*uz*cm-ux*s;
         uz*ux*cm-uy*s    uz*uy*cm+ux*s      c+uz*uz*cm]
end

function rotation3{T}(axis::Vector{T}, angle)
    n = norm(axis)
    axisn = n>0 ? axis/n : (tmp = zeros(T,length(axis)); tmp[1] = 1)
    _rotation3(axisn, angle)
end

# Angle/axis representation where the angle is the norm of the vector (so axis is not normalized)
function rotation3{T}(axis::Vector{T})
    n = norm(axis)
    axisn = n>0 ? axis/n : (tmp = zeros(typeof(one(T)/1),length(axis)); tmp[1] = 1; tmp)
    _rotation3(axisn, n)
end

# Convert rotation matrix to axis representation (angle/axis with the norm of the vector coding the angle)
function rotationparameters{T}(R::Matrix{T})
    size(R, 1) == size(R, 2) || error("Matrix must be square")
    if size(R, 1) == 2
        return [atan2(-R[1,2],R[1,1])]
    elseif size(R, 1) == 3
        # See Itzhack Y. Bar-Itzhack.  "New Method for Extracting the Quaternion from a Rotation Matrix",
        # Journal of Guidance, Control, and Dynamics, Vol. 23, No. 6 (2000), pp. 1085-1087.
        # On https://en.wikipedia.org/wiki/Rotation_matrix#Conversions
        # This is said to be robust even if R is not orthogonal
        Q12p, Q13p, Q23p = R[1,2]+R[2,1], R[1,3]+R[3,1], R[2,3]+R[3,2]
        Q12m, Q13m, Q23m = R[1,2]-R[2,1], R[1,3]-R[3,1], R[2,3]-R[3,2]
        K = [R[1,1]-R[2,2]-R[3,3]  Q12p  Q13p  Q23m;
             Q12p  R[2,2]-R[1,1]-R[3,3]  Q23p -Q13m;
             Q13p  Q23p  R[3,3]-R[1,1]-R[2,2]  Q12m;
             Q23m -Q13m  Q12m  R[1,1]+R[2,2]+R[3,3]]
        D, V = eig(K)
        imax = indmax(D)
        q = V[:,imax]
        ax = q[1:3]
        n = norm(ax)
        if n > 0
            ax = ax/n
        end
        th = -2*atan(n/q[4])
        return th*ax
    else
        error("Rotations in $(size(R, 1)) dimensions not supported")
    end
end

tformtranslate(trans::Vector) = AffineTransform(eye(length(trans)), trans)

tformscale(scale::Real, ndims::Int) = AffineTransform(scale*eye(typeof(scale),ndims), zeros(typeof(scale),ndims))

function tformrotate(axis::Vector, angle)
    if length(axis) == 3
        return AffineTransform(rotation3(axis, angle), zeros(eltype(axis),3))
    else
        error("Dimensionality ", length(axis), " not supported")
    end
end

function tformrotate(x::Vector)
    if length(x) == 3
        return AffineTransform(rotation3(x), zeros(eltype(x),3))
    else
        error("Dimensionality ", length(x), " not supported")
    end
end

# The following is useful in optimization routines, which tend to work on parameter vectors
function tformrigid(p::Vector)
    if length(p) == 3
        # 2d: [angle, shiftx, shifty]
        return AffineTransform(rotation2(p[1]), p[2:3])
    elseif length(p) == 6
        # 3d: [axis, shift]
        return AffineTransform(rotation3(p[1:3]), p[4:6])
    else
        error("Rigid transformations with $(length(p)) parameters not supported")
    end
end

function(*)(T2::AffineTransform, T1::AffineTransform)
    A1 = T1.scalefwd
    A2 = T2.scalefwd
    o = T2.offset + A2*T1.offset
    AffineTransform(A2*A1, o)
end

function(\)(T2::AffineTransform, T1::AffineTransform)
    A1 = T1.scalefwd
    A2 = T2.scaleinv
    o = A2*(T1.offset-T2.offset)
    AffineTransform(A2*A1, o)
end


text:
# Autogenerated wrapper script for OpenLibm_jll for aarch64-apple-darwin
export libopenlibm

JLLWrappers.@generate_wrapper_header("OpenLibm")
JLLWrappers.@declare_library_product(libopenlibm, "@rpath/libopenlibm.4.dylib")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libopenlibm,
        "lib/libopenlibm.4.0.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
#=
Functions for constructing initial states
=#
"""
    plus_state(n)

Return the state ``\\ket{+}^{\\otimes n}`` in the computational basis
"""
function plus_state(n::Int)
    @assert n > 0 "n must be a positive integer"
    spinup = Vector{Complex{Float64}}([1., 0.])
    spindown = Vector{Complex{Float64}}([0., 1.])

    return kron([(spinup + spindown) / sqrt(2.) for i in 1:n]...)
end

"""
    ghz_state(n)

Return the GHZ state ``(\\ket{00...0} + \\ket{11...1}) / \\sqrt{2}`` in the computational basis
"""
function ghz_state(n::Int)
    @assert n > 0 "n must be a positive integer"
    r = zeros(Complex{Float64}, 2^n)
    r[1] = r[end] = 1. / sqrt(2.)
    return r
end


text:
module JungleHelperWipeDropdown

using ..Ahorn

# Inject Jungle Helper wipes into the dropdown in map metadata.
Ahorn.MetadataWindow.metaDropdownOptions["Wipe"]["Arrow (Jungle Helper + max480's Helping Hand)"] = "MaxHelpingHand/CustomWipe:JungleHelper/Arrow"
Ahorn.MetadataWindow.metaDropdownOptions["Wipe"]["Dots (Jungle Helper + max480's Helping Hand)"] = "MaxHelpingHand/CustomWipe:JungleHelper/Dots"
Ahorn.MetadataWindow.metaDropdownOptions["Wipe"]["Rocks (Jungle Helper + max480's Helping Hand)"] = "MaxHelpingHand/CustomWipe:JungleHelper/Rocks"
Ahorn.MetadataWindow.metaDropdownOptions["Wipe"]["Spinners (Jungle Helper + max480's Helping Hand)"] = "MaxHelpingHand/CustomWipe:JungleHelper/Speen"
Ahorn.MetadataWindow.metaDropdownOptions["Wipe"]["Vines (Jungle Helper + max480's Helping Hand)"] = "MaxHelpingHand/CustomWipe:JungleHelper/Vines"

end

max_stars_repo_name: Jackal-Celeste/JungleHelper

max_issues_repo_name: Jackal-Celeste/JungleHelper

max_forks_repo_name: Jackal-Celeste/JungleHelper


text:
module TaskHeritableStorageTests #end

using Test

using TaskHeritableStorage

# Examples
struct ShallowCopyKey{sym} end
TaskHeritableStorage.fork_task_local_storage(k::ShallowCopyKey, v) = (k,v)

const key = ShallowCopyKey{:tls}()

fetch(@async begin  # Put each test in its own Task so they don't share a task_heritable_storage
    task_local_storage()[key] = 5
    fetch(@async begin
        @test task_local_storage()[key] == 5
    end)
end)

@testset "assigning doesn't affect parent tasks" begin
    @sync @async begin
        task_local_storage()[key] = 1
        @sync @async begin
            @test task_local_storage()[key] == 1

            # Update it to 2 within this Task
            task_local_storage()[key] = 2
            @test task_local_storage()[key] == 2
        end
        # In the parent task, the value is still 1
        @test task_local_storage()[key] == 1
    end
    # In the outermost task, the value was never set
    @test !haskey(task_local_storage(), key)
end

@testset "callbacks" begin
    # Dummy function that "uses concurrency" when sorting
    function psort(args...; kwargs...)
        fetch(@async sort(args...; kwargs...))
    end

    @testset "PROBLEM: normal task_local_storage isn't composable" begin
        task_local_storage()[:reverse] = false
        reversed() = task_local_storage()[:reverse]

        less_than = (a,b) -> reversed() ? a>b : a<b
        # All good, we can access this task_local_storage
        @test sort(1:10, lt=less_than) == 1:10

        # SURPRISE: the author of sort decides to parallelize it (represented via `psort`)!
        @test_throws Exception psort(1:10, lt=less_than) == 1:10
        # The error is a KeyError, because :reverse is no longer in the task_local_storage:
        e = try                          psort(1:10, lt=less_than) == 1:10; catch e; e end
        @static if VERSION >= v"1.3-"
            @test e.task.exception == KeyError(:reverse)
        else
            @test e == KeyError(:reverse)
        end  # VERSION
    end

    @testset "SOLUTION: Using a `fork_task_local_storage()` override _is_ composable" begin
        reversekey = ShallowCopyKey{:reverse}()

        task_local_storage()[reversekey] = false
        reversed() = task_local_storage()[reversekey]

        less_than = (a,b) -> reversed() ? a>b : a<b
        @test sort(1:10, lt=less_than) == 1:10

        # NOW, if the author of sort wants to parallelize it, everything is *okay*! :)
        @test psort(1:10, lt=less_than) == 1:10
    end
end

@testset "functional interface" begin
    @test haskey(task_local_storage(), key) == false
    task_local_storage(key, 1) do
        # Same Task
        @test task_local_storage()[key] == 1
        # Nested Task
        @test fetch(@async task_local_storage()[key]) == 1
    end
    # Storage removed after the function returns
    @test haskey(task_local_storage(), key) == false
end

end


text:
Transducers.transduce(xf, rf, init, xs, ex::DepthFirstEx) =
    _transduce_depthfirst(xf, rf, init, xs; ex.kwargs...)

function _transduce_depthfirst(
    xform::Transducer,
    step::F,
    init,
    coll0;
    simd::SIMDFlag = Val(false),
    basesize::Union{Integer,Nothing} = nothing,
    stoppable::Union{Bool,Nothing} = nothing,
) where {F}
    rf0 = _reducingfunction(xform, step; init = init)
    rf, coll = retransform(rf0, coll0)
    if stoppable === nothing
        stoppable = _might_return_reduced(rf, init, coll)
    end
    acc = @return_if_reduced _reduce_df(
        stoppable ? CancellableDACContext() : NoopDACContext(),
        DummyTask(),
        maybe_usesimd(rf, simd),
        init,
        SizedReducible(
            coll,
            basesize === nothing ? amount(coll) ÷ Threads.nthreads() : basesize,
        ),
    )
    result = complete(rf, acc)
    if unreduced(result) isa DefaultInitOf
        throw(EmptyResultError(rf))
    end
    return result
end

struct DummyTask end
Base.schedule(::DummyTask) = nothing

function _reduce_df(ctx, next_task, rf::R, init::I, reducible::Reducible) where {R,I}
    if should_abort(ctx)
        # As other tasks may be calling `fetch` on `next_task`, it
        # _must_ be scheduled at some point to avoid dead lock:
        schedule(next_task)
        # Maybe use `error=false`?  Or pass something and get it via `yieldto`?
        return init
    end
    if issmall(reducible)
        schedule(next_task)
        acc = _reduce_basecase(rf, init, reducible)
        if acc isa Reduced
            cancel!(ctx)
        end
        return acc
    else
        left, right = _halve(reducible)
        fg, bg = splitcontext(ctx)
        task = nonsticky!(@task _reduce_df(bg, next_task, rf, init, right))
        a0 = _reduce_df(fg, task, rf, init, left)
        b0 = fetch(task)
        a = @return_if_reduced a0
        should_abort(ctx) && return a  # slight optimization
        b0 isa Reduced && return combine_right_reduced(rf, a, b0)
        return combine(rf, a, b0)
    end
end


text:
@testset "elliptical slice" begin

    mu = [0., 0.]
    cov = [1. 0.; 0. 1.]
    @gen function foo()
        x = @trace(mvnormal(mu, cov), :x)
        y = @trace(mvnormal(x, [0.01 0.; 0. 0.01]), :y)
    end

    # smoke test
    trace, = generate(foo, (), choicemap((:y, [1., 1.])))
    trace = elliptical_slice(trace, :x, mu, cov)
end


text:
###################################
#=
Simple implementations, for dimension 1, of the quadratic estimate, Cvar and Cbias which can be used
to generate estimates for non-local invariant models. Used in `make_figure_a.jl` and `make_figure_b.jl`

=#
######################################

function local_Aℓ(islocvar::Bool, parms)
    ξk, C1k, CZZmobsk = parms.ξk, parms.C1k, parms.CZZmobsk
    deltx, deltk = parms.deltx, parms.deltk
    mlt = islocvar ? -1 : 1
    CXXinvCZZ = squash(1 ./ CZZmobsk)
    Apqx = ifftd( C1k[1] .* conj(C1k[1]) .* CXXinvCZZ, deltk)
    Bx   = ifftd(                    1.0 .* CXXinvCZZ, deltk)
    Cpx  = ifftd(                 C1k[1] .* CXXinvCZZ, deltk)
    Dpx  = ifftd(           conj(C1k[1]) .* CXXinvCZZ, deltk)
    ABCDx = zeros(Complex{Float64}, size(CZZmobsk))
    for ix in eachindex(ABCDx)
        ABCDx[ix] = 2 * Apqx[ix] * Bx[ix] + mlt * Cpx[ix] * Cpx[ix] + mlt * Dpx[ix] * Dpx[ix]
    end
    rtnk = 1./ abs2(ξk[1]) ./ real(fftd(ABCDx, deltx))
    squash!(rtnk)
    return rtnk
end
#= ---- Basic Test -----
include("make_figure_a.jl",1:53)
test1 = (2/√(2π)) .* NonstationaryPhase.local_Aℓ(true, parms)
test2 = (2/√(2π)) ./ invAℓfun(parms.CZZmobsk, parms)
test3 = Cℓvarfun(parms)

test1./test2 # should be [*,1,1 ..., 1]
test1./test3 # should be [*,1,1 ..., 1]
=#

function local_hatθfun(islocvar::Bool, xk, yk, Aell, parms)
    ξk, C1k, CZZmobsk = parms.ξk, parms.C1k, parms.CZZmobsk
    deltx, deltk = parms.deltx, parms.deltk
    mlt = islocvar ? -1 : 1
    xkinvCZZ = xk ./ CZZmobsk
    ykinvCZZ = yk ./ CZZmobsk
    squash!(xkinvCZZ)
    squash!(ykinvCZZ)
    A1x  = ifftd(                 xkinvCZZ, deltk)
    Dpx  = ifftd(      C1k[1]  .* ykinvCZZ, deltk)
    A2x  = ifftd(                 ykinvCZZ, deltk)
    Cpx  = ifftd( conj(C1k[1]) .* xkinvCZZ, deltk)
    rtnk   = conj(ξk[1]) .* fftd(A1x.*Dpx + mlt.*A2x.*Cpx, deltx)
    rtnk .*= Aell
    squash!(rtnk)
    return rtnk
end
#= ---- Test local_hatθfun on locally invariant models -----
include("make_figure_a.jl",1:53)
Aelltest = NonstationaryPhase.local_Aℓ(true, parms)
xk = randn(parms.nside) + im * randn(parms.nside)
yk = randn(parms.nside) + im * randn(parms.nside)
test1 = NonstationaryPhase.local_hatθfun(true, xk, yk, Aelltest, parms)
test2 = let
    rtnk   = unnormalized_estϕkfun(xk, yk, parms)
    rtnk ./= invAℓfun(parms)
    squash!(rtnk)
    rtnk
end

real(test1)./real(test2)  # should be [*,1,1 ..., 1]
imag(test1)./imag(test2)  # should be [*,1,1 ..., 1]
=#

# ---- Test local_hatθfun on nonlocally invariant models -----
function test_local_hatθfun(sims, ϕx, parmsNLI)
    dm     = ndims(parmsNLI.x)
    Aell   = NonstationaryPhase.local_Aℓ(false, parmsNLI)
    Eestϕk = zeros(Complex{Float64}, size(ϕx))
    for iter = 1:sims
        wk     =  fftd(randn(parmsNLI.nside)./√(parmsNLI.deltx ^ dm), parmsNLI.deltx)
        xk     = √(parmsNLI.Ck) .* wk
        yktmp  = √(2π) * (parmsNLI.C1k[1,1] ./ √(parmsNLI.Ck)) .* wk
        yx     = ϕx .* ifftd(yktmp, parmsNLI.deltk)
        yk     = fftd(yx, parmsNLI.deltx)
        estϕk  = NonstationaryPhase.local_hatθfun(false, xk, yk, Aell, parmsNLI)
        estϕk  += NonstationaryPhase.local_hatθfun(false, yk, xk, Aell, parmsNLI)
        Eestϕk += estϕk ./ sims
    end
    Eestϕx = ifftdr(Eestϕk, parmsNLI.deltk)
    return Eestϕx, Eestϕk
end
#= ----- this runs an example of the test using test_local_hatθfun ----
    using NPhaseGRF
    using PyPlot
    const d_const      = 1
    const period_const = 2π
    const nside_const  = nextprod([2,3,5,7], 2_500)
    ν  = 3.0
    ρ  = 0.02
    νΦ = 5.0
    ρΦ = 0.5
    σΦ = 0.002
    parmsNLI               = PhaseParms(d_const, period_const, nside_const)
    parmsNLI.Ck[:]         =   maternk(parmsNLI.k; ν=ν, ρ=ρ)
    parmsNLI.C1k[1,1][:]   = parmsNLI.Ck .* (2 - abs2(parmsNLI.r.*ρ)) ./ (4ν + abs2(parmsNLI.r.*ρ))
    parmsNLI.C1k[1,1][:] .*= (2ν / ρ) / (2*√(2π))
    parmsNLI.ξk[:]         = Array{Complex{Float64},1}[ ones(parmsNLI.k[1]) ]
    parmsNLI.Cϕϕk[:]      =  maternk(parmsNLI.k; ν=νΦ, ρ=ρΦ, σ=σΦ)
    parmsNLI.CZZmk[:]      =  parmsNLI.Ck[:]
    parmsNLI.CZZmobsk[:]   = parmsNLI.CZZmk + parmsNLI.CNNk
    parmsNLI.CZZmobsk[parmsNLI.r .< 1]                 = Inf
    parmsNLI.CZZmobsk[parmsNLI.r .>= 0.9*parmsNLI.nyq] = Inf

    ϕx = NonstationaryPhase.grfsimx(parmsNLI.Cϕϕk, parmsNLI.deltx, parmsNLI.deltk)
    ϕx += -minimum(real(ϕx)) + ρ

    Eestϕx, Eestϕk = NonstationaryPhase.test_local_hatθfun(5000, ϕx, parmsNLI)

    plot(fftshift(Eestϕx), "k", alpha = 0.4)
    plot(fftshift(ϕx), "b")
=#

function local_θ2bias(islocvar::Bool, parms)
    ξk, CZZmobsk, C2k  = parms.ξk,  parms.CZZmobsk, parms.C2k[1,1]
    deltx, deltk, xco, kco = parms.deltx, parms.deltk, parms.x, parms.k
    m1s   = ones(parms.Cϕϕk)
    mlt = islocvar ? -1 : 1
    Aell = local_Aℓ(islocvar, parms)
    hatθC2_kplusℓ  =  local_hatθfun(islocvar, C2k, m1s, Aell, parms)
    rtnk     = zeros(Float64, size(CZZmobsk))
    tmp_rtnk = zeros(Float64, size(CZZmobsk))
    for ωindx in eachindex(kco[1])
        ω = kco[1][ωindx]
        hatθC2_kplusℓminusω =  local_hatθfun(islocvar, shiftfk_by_ω(C2k, -ω, parms), m1s, Aell, parms)
        tmp_rtnk   =  abs2(hatθC2_kplusℓ + mlt.*hatθC2_kplusℓminusω)
        tmp_rtnk .*=  shiftfk_by_ω(ξk[1] .* conj(ξk[1]) .* parms.Cϕϕk, -ω, parms)
        tmp_rtnk .*=  ξk[1][ωindx] * conj(ξk[1][ωindx]) * parms.Cϕϕk[ωindx]
        rtnk     .+=  tmp_rtnk
    end
    rtnk .*= 2 * 4 * (deltk / 2π)
    # the extra 2 in the above constant is because in 1-d, the two terms involving Cθθk.*Cθθk_kminusω are the same.
    return rtnk
end
#= ---- Basic Test -----
include("make_figure_a.jl",1:53)
Cℓbias_test1 = NonstationaryPhase.local_θ2bias(true,  parms)
Cℓbias_test2 = approxCℓbiasfun(parms)
Cℓbias_test3 = Cℓbiasfun(parms)
Cℓbias_test1 ./ Cℓbias_test3   # should be [1,1,1 ..., 1]
Cℓbias_test2 ./ Cℓbias_test3   # should be [1,1,1 ..., 1]

=#


text:
using DataArrays
using DataFrames
using Gadfly

function θestimation_()
    # 1970 batting averages for 18 major league players. first column = batting average for first 45 at bats 2nd column = batting everage for remainder of season
    data = [.400 .346
            .378 .298
            .356 .276
            .333 .222
            .311 .273
            .311 .270
            .289 .263
            .267 .210
            .244 .269
            .244 .230
            .222 .264
            .222 .256
            .222 .303
            .222 .264
            .222 .226
            .200 .285
            .178 .316
            .156 .200]

    # data transformation
    y = data[:,1]
    ytest = data[:,2]
    n = 45
    # arcsin transform
    x = sqrt(n).*asin(2.*y-1)

    # shrinkage estimate
    μBar = mean(x)
    V    = sumabs2(x - μBar)
    d    = length(x)
    s2   = V/d
    # by construction of the arcsin transform
    σ2   = 1
    B    = σ2/(σ2 + max(0, s2-σ2))
    μShrunk = μBar + (1-B).*(x-μBar)

    # back transform
    θShrunk = 0.5*(sin(μShrunk/sqrt(n))+1)
    θMLE    = y
    (ytest, θShrunk, θMLE)
end

"""
Baseball Batting Average Shrinkage Estimates
Reproduce example from "Data Analysis Using Stein's Estimator and its Generalizations" Bradley Efron; Carl Morris JASA Vol. 70, No. 350. (Jun., 1975), pp. 311-319.
"""
function shrinkageDemoBaseball()
    θtrue, θShrunk, θMLE = θestimation_()
    # plot Shrinkage Estimates
    p1 = plot([layer(x=[θMLE[i] θShrunk[i]], y=[1 0], Geom.line) for i=1:d]...,
              layer(x=θMLE, y=ones(1, d), Geom.point),
              layer(x=θShrunk, y=zeros(1, d), Geom.point),
              Theme(default_color=colorant"blue"))

    # histograms
    df = DataFrame(PN=[1:d, 1:d, 1:d], Theta=[ytest, θShrunk, θMLE], Type=[ones(d), 2*ones(d), 3*ones(d)])
    p2 = plot(df, y="Theta", x="PN", color="Type", Geom.bar(position=:dodge))

    [p1, p2]
end

function shrinkageMSEBaseball()
    θtrue, θShrunk, θMLE = θestimation_()
    mseMLE = mean((θtrue-θMLE).^2)
    mseShrink = mean((θtrue-θShrunk).^2)
    [mseMLE, mseShrink, mseMLE/mseShrink]
end


text:
@test isequal(deparse(:(cos(x) + sin(x))), "cos(x) + sin(x)")
@test isequal(deparse(:(cos(x) + sin(x) + exp(-x))), "cos(x) + sin(x) + exp(-x)")
@test isequal(deparse(parse("x+y*z")), "x + y * z")
@test isequal(deparse(parse("(x+y)*z")), "(x + y) * z")
@test isequal(deparse(parse("1/(x/y)")), "1 / (x / y)")
@test isequal(deparse(parse("1/(x*y)")), "1 / (x * y)")
@test isequal(deparse(parse("z^(x+y)")), "z ^ (x + y)")
@test isequal(deparse(parse("z^(x*y)")), "z ^ (x * y)")


text:
# Show AbstractAxis types
Base.show(io::IO, ::MIME"text/plain", ::Axis{IdxMap}) where IdxMap = print(io, "Axis$IdxMap")
Base.show(io::IO, ::Axis{IdxMap}) where IdxMap = print(io, "Axis$IdxMap")

Base.show(io::IO, ::FlatAxis) = print(io, "FlatAxis()")
Base.show(io::IO, ::MIME"text/plain", ::FlatAxis) = print(io, "FlatAxis()")

Base.show(io::IO, ::NullAxis) = print(io, "NullAxis()")

Base.show(io::IO, ::MIME"text/plain", ::PartitionedAxis{PartSz, IdxMap, Ax}) where {PartSz, IdxMap, Ax} =
    print(io, "PartitionedAxis($PartSz, $(Ax()))")
Base.show(io::IO, ::PartitionedAxis{PartSz, IdxMap, Ax}) where {PartSz, IdxMap, Ax} =
    print(io, "PartitionedAxis($PartSz, $(Ax()))")

Base.show(io::IO, ::ShapedAxis{Shape, IdxMap}) where {Shape, IdxMap} =
    print(io, "ShapedAxis($Shape, $IdxMap)")

Base.show(io::IO, ::MIME"text/plain", ::ViewAxis{Inds, IdxMap, Ax}) where {Inds, IdxMap, Ax} =
    print(io, "ViewAxis($Inds, $(Ax()))")
Base.show(io::IO, ::ViewAxis{Inds, IdxMap, <:Ax}) where {Inds, IdxMap, Ax} =
    print(io, "ViewAxis($Inds, $(Ax()))")
Base.show(io::IO, ::ViewAxis{Inds, IdxMap, <:NullorFlatAxis}) where {Inds, IdxMap} =
    print(io, Inds)

Base.show(io::IO, ci::ComponentIndex) = print(io, "ComponentIndex($(ci.idx), $(ci.ax))")

# Show ComponentArrays
_print_type_short(io, ca; color=:normal) = _print_type_short(io, typeof(ca); color=color)
_print_type_short(io, T::Type; color=:normal) = printstyled(io, T; color=color)
_print_type_short(io, ::Type{<:ComponentArray{T,N,<:Array}}; color=:normal) where {T,N} = printstyled(io, "ComponentArray{$T,$N}"; color=color) # do not pollute the stacktrace with verbose type printing
_print_type_short(io, ::Type{<:ComponentArray{T,1,<:Array}}; color=:normal) where {T} = printstyled(io, "ComponentVector{$T}"; color=color)
_print_type_short(io, ::Type{<:ComponentArray{T,2,<:Array}}; color=:normal) where {T} = printstyled(io, "ComponentMatrix{$T}"; color=color)
_print_type_short(io, ::Type{<:ComponentArray{T,N,<:SubArray}}; color=:normal) where {T,N} = printstyled(io, "ComponentArray{$T,$N,SubArray...}"; color=color) # do not pollute the stacktrace with verbose type printing
_print_type_short(io, ::Type{<:ComponentArray{T,1,<:SubArray}}; color=:normal) where {T} = printstyled(io, "ComponentVector{$T,SubArray...}"; color=color)
_print_type_short(io, ::Type{<:ComponentArray{T,2,<:SubArray}}; color=:normal) where {T} = printstyled(io, "ComponentMatrix{$T,SubArray...}"; color=color)

@static if VERSION >= v"1.6"
    Base.print_type_stacktrace(io, CA::Type{<:ComponentArray}; color=:normal) = _print_type_short(io, CA; color=color)
end

function Base.show(io::IO, x::ComponentVector)
    print(io, "(")
    for (i,key) in enumerate(keys(x))
        if i==1
            print(io, "$key = ")
        else
            print(io, ", $key = ")
        end
        show(io, x[key])
    end
    print(io, ")")
    return nothing
end

function Base.show(io::IO, mime::MIME"text/plain", x::ComponentVector)
    len = length(x)
    ax = getaxes(x)[1]
    if last_index(ax) == len
        _print_type_short(io, x)
        show(io, x)
    else
        print(io, "$len-element ")
        _print_type_short(io, x)
        println(io, " with axis $ax:")
        Base.print_array(io, getdata(x))
    end
    return nothing
end

function Base.show(io::IO, ::MIME"text/plain", x::ComponentMatrix{T,A,Axes}) where {T,A,Axes}
    if !haskey(io, :compact) && length(axes(x, 2)) > 1
        io = IOContext(io, :compact => true)
    end
    axs = getaxes(x)
    sz = size(x)
    print(io, "$(sz[1])×$(sz[2]) ")
    _print_type_short(io, x)
    println(io, " with axes $(axs[1]) × $(axs[2])")
    Base.print_matrix(io, getdata(x))
    return nothing
end


text:
const Orbits = [(1,0), (2,0), (0,1), (1,1), (0,2)]

"""
    flatten(pdm::Dict{Int, <:PersistenceDiagram}) -> PersistenceDiagram

Flatten a multi-dimensional pd into a barcode.
"""
flatten(pdm::Dict{Int, PD}) where {T<:Real, PD<:PersistenceDiagram{T}} =
    (collect∘Iterators.flatten∘values)(pdm)

"""
    maxsize(pds) -> Int

Return the number of bars in the largest barcode from the collection of barcodes `pds`
"""
maxsize(pds::AbstractVector{PD}) where {T<:Real, PD<:PersistenceDiagram{T}} =
    maximum(map(length, pds))

tropicdim(n) = n+(div(n*(n+1),2)) # codomain dimension

function tropicstats(pds::AbstractVector{PD}) where {T<:Real, PD<:PersistenceDiagram{T}}
    n = ComputationalHomology.maxsize(pds)
    d = ComputationalHomology.tropicdim(n)
    maxval = filter!(!isinf, map(i->max(birth(i),death(i)), Iterators.flatten(pds))) |> maximum
    m =  map(i->birth(i)/-ComputationalHomology.birthx(i), Iterators.flatten(pds)) |> maximum
    (d, n, ceil(Int, m), maxval)
end

"""
Caluclate 2-symmetric tropical rational polynomial for the persistence diagram
"""
function tropicp(pd::PersistenceDiagram{T},m,l,p; maxval=T(Inf)) where {T<:Real}
    B = Tuple{T,T}[]
    for i in pd
        d = min(maxval, death(i))
        b = birth(i)
        push!(B, (min(m*(d-b), b), d-b))
    end
    # pad zeros
    while length(B) < l+p
        push!(B, (zero(T), zero(T)))
    end
    n = length(B)

    if l == 0 && p == 0
        0
    elseif l == 0
        sort!(map(sum, B), rev=true)[1:p] |> sum
    elseif p == 0
        sort!(map(last,B), rev=true)[1:l] |> sum
    else
        res = 0
        L = hcat(combinations(1:n,l)...)'
        for i in 1:length(L)
            ii = L[i,:]
            e0 = sum(last, B[ii])
            K = hcat(combinations(setdiff(1:n, ii),p)...)'
            for j = 1:length(K)
                jj = K[j,:]
                e1 = sum(map(first, B[jj]) .+ map(last, B[jj]))
                res = max(res, e0 + e1)
            end
        end
        res
    end
end

"""
    tropic(pd)

Calculate a tropical coordinates for a persistent diagram `pd`

Paper: "Tropical Sufficient Statistics for Persistent Homology" by A.Monod et al.
Ref: https://arxiv.org/abs/1709.02647
"""
function tropic(pds::AbstractVector{PD}) where {T<:Real, PD<:PersistenceDiagram{T}}
    d,n,m,mv = ComputationalHomology.tropicstats(pds)
    [tropicp(pd,m,i,j,maxval=mv) for (i,j) in ComputationalHomology.Orbits, pd in pds ]
end
tropic(pd::PersistenceDiagram{T}) where {T<:Real} = tropic([pd]) |> vec
tropic(pdd::Dict{Int, PD}) where {T<:Real, PD<:PersistenceDiagram{T}} = tropic(flatten(pdd))


text:
#example.jl
function echomsg(msg)
    return join([msg, "Thank you !"], ". ")
end

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
using Pkg
Pkg.activate(pwd())

# # Constrained Optimization
# ## Numerical method

using Plots, Random

function create_anim(
    f,
    path,
    xlims,
    ylims,
    file_name = joinpath(pwd(), randstring(12) * ".gif");
    xbounds = xlims,
    ybounds = ylims,
    fps = 15,
)
    xs = range(xlims...; length = 100)
    ys = range(ylims...; length = 100)
    plt = contourf(xs, ys, f; color = :jet)

    # add constraints if provided
    if !(xbounds == xlims && ybounds == ylims)
        x_rect = [xbounds[1]; xbounds[2]; xbounds[2]; xbounds[1]; xbounds[1]]
        y_rect = [ybounds[1]; ybounds[1]; ybounds[2]; ybounds[2]; ybounds[1]]

        plot!(x_rect, y_rect; line = (2, :dash, :red), label="")
    end

    # add an empty plot
    plot!(Float64[], Float64[]; line = (4, :arrow, :black), label = "")

    # extract the last plot series
    plt_path = plt.series_list[end]

    # create the animation and save it
    anim = Animation()
    for x in eachcol(path)
        push!(plt_path, x[1], x[2]) # add a new point
        frame(anim)
    end
    gif(anim, file_name; fps = fps, show_msg = false)
    return nothing
end

f(x) = sin(x[1] + x[2]) + cos(x[1])^2
f(x1,x2) = f([x1;x2])
g(x) = [cos(x[1] + x[2]) - 2*cos(x[1])*sin(x[1]); cos(x[1] + x[2])]

#+

function optim(f, g, P, x, α; max_iter=100)
    xs = zeros(length(x), max_iter+1)
    ys = zeros(length(x), max_iter)
    xs[:,1] = x
    for i in 1:max_iter
        ys[:,i] = xs[:,i] - α*g(xs[:,i])
        xs[:,i+1] = P(ys[:,i])
    end
    return xs, ys
end

#+

P(x, x_min, x_max) = min.(max.(x, x_min), x_max)

x_min = [-1; -1]
x_max = [0; 0]

xs, ys = optim(f, g, x -> P(x,x_min,x_max), [0;-1], 0.1)

xlims = (-3, 1)
ylims = (-2, 1)

#+

create_anim(f, xs, xlims, ylims, joinpath(pwd(), "lecture_08", "anim6.gif");
    xbounds=(x_min[1], x_max[1]),
    ybounds=(x_min[2], x_max[2]),
)

#+

xys = hcat(reshape([xs[:,1:end-1]; ys][:], 2, :), xs[:,end])

create_anim(f, xys, xlims, ylims, joinpath(pwd(), "lecture_08", "anim7.gif");
    xbounds=(x_min[1], x_max[1]),
    ybounds=(x_min[2], x_max[2]),
)