text:
using ADI
using PyCall
using SiriusB

fits = pyimport("astropy.io.fits")
pro = pyimport("proplot")

# rcParams
pro.rc["image.origin"] = "lower"
pro.rc["image.cmap"] = "inferno"
pro.rc["grid"] = false

epochs = [
    "2020feb04",
    "2020nov21",
    "2020nov28"
]
res = [
    fits.getdata(datadir("epoch_2020feb04", "residuals", "2020feb04_sirius-b_residual_median.fits")),
    fits.getdata(datadir("epoch_2020nov21", "residuals", "2020nov21_sirius-b_residual_median.fits")),
    fits.getdata(datadir("epoch_2020nov28", "residuals", "2020nov28_sirius-b_residual_annular_pca-2.fits")),
]
angles = [
    fits.getdata(datadir("epoch_$epoch", "processed", "$(epoch)_sirius-b_pa.fits"))
    for epoch in epochs
]
fwhms = [
    7.985399074538929,
    7.644946278252633,
    8.216162087416189
]

flat_res = collapse.(res, angles)
sigs = detectionmap.(significance, flat_res, fwhms)
stims = stimmap.(res, angles)

tick_locs = range(25, 174, length=5)
tick_labs = @. string(round(SiriusB.auscale * (tick_locs - 99.5), digits=1))

py"""
import proplot as pro
from matplotlib import patches
fig, axs = pro.subplots(ncols=3, wspace="1em", figwidth="7.5in")

axs[0].imshow($(flat_res[1]))
axs[0].text(6, 6, "median", color="w")
axs[0].format(title="Epoch 2020-02-04")

axs[1].imshow($(flat_res[2]))
axs[1].text(6, 6, "median", color="w")
axs[1].format(title="Epoch 2020-11-21")

axs[2].imshow($(flat_res[3]))
axs[2].text(6, 6, "annular PCA(2)", color="w")
axs[2].format(title="Epoch 2020-11-28")

axs.format(
    xticks=$tick_locs,
    xticklabels=$tick_labs,
    yticks=$tick_locs,
    yticklabels=$tick_labs,
    xlabel="x [AU]",
    ylabel="y [AU]",
)

fig.save($(figuredir("residuals.pdf")))
"""

py"""
import proplot as pro
fig, axs = pro.subplots(ncols=3, nrows=2, wspace="1em", figwidth="7.5in")

axs[0, 0].imshow($(sigs[1]), vmin=0, vmax=5)#, colorbar="r", colorbar_kw=dict(space=0))
axs[0, 0].text(6, 6, "median", color="w")
axs[0, 0].format(title="Epoch 2020-02-04")
axs[0, 1].imshow($(sigs[2]), vmin=0, vmax=5)#, colorbar="r", vmin=0, colorbar_kw=dict(space=0))
axs[0, 1].text(6, 6, "median", color="w")
axs[0, 1].format(title="Epoch 2020-11-21")
m= axs[0, 2].imshow($(sigs[3]), vmin=0, vmax=5)#, colorbar="r", vmin=0, colorbar_kw=dict(space=0))
axs[0, 2].text(6, 6, "annular PCA(2)", color="w")
axs[0, 2].format(title="Epoch 2020-11-28")
axs[0, 2].colorbar(m, loc="r", label="significance")

axs[1, 0].imshow($(stims[1]), vmin=0, vmax=1)#, colorbar="r", colorbar_kw=dict(space=0))
axs[1, 0].text(6, 6, "median", color="w")
axs[1, 1].imshow($(stims[2]), vmin=0, vmax=1)#, colorbar="r", vmin=0, colorbar_kw=dict(space=0))
axs[1, 1].text(6, 6, "median", color="w")
m= axs[1, 2].imshow($(stims[3]), vmin=0, vmax=1)#, colorbar="r", vmin=0, colorbar_kw=dict(space=0))
axs[1, 2].text(6, 6, "annular PCA(2)", color="w")
axs[1, 2].colorbar(m, loc="r", label="STIM probability")

axs.format(
    xticks=$tick_locs,
    xticklabels=$tick_labs,
    yticks=$tick_locs,
    yticklabels=$tick_labs,
    xlabel="x [AU]",
    ylabel="y [AU]",
)

fig.save($(figuredir("detection_maps.pdf")))
pro.close(fig)

fig, axs = pro.subplots(ncols=3, wspace="1em", figwidth="7.5in")

axs[0].imshow($(sigs[1]), vmin=3)#, colorbar="r", colorbar_kw=dict(space=0))
axs[0].text(6, 6, "median", color="w")
axs[0].format(title="Epoch 2020-02-04")
axs[1].imshow($(sigs[2]), vmin=3)#, colorbar="r", vmin=0, colorbar_kw=dict(space=0))
axs[1].text(6, 6, "median", color="w")
axs[1].format(title="Epoch 2020-11-21")
m= axs[2].imshow($(sigs[3]), vmin=3)#, colorbar="r", vmin=0, colorbar_kw=dict(space=0))
axs[2].text(6, 6, "annular PCA(2)", color="w")
axs[2].format(title="Epoch 2020-11-28")
fig.colorbar(m, loc="r", label="significance")

axs.format(
    xticks=$tick_locs,
    xticklabels=$tick_labs,
    yticks=$tick_locs,
    yticklabels=$tick_labs,
    xlabel="x [AU]",
    ylabel="y [AU]",
)

fig.save($(figuredir("sig_threshold.pdf")))
pro.close(fig)
"""

py"""
import proplot as pro
# fig, axs = pro.subplots(ncols=3, wspace="1em", figwidth="7.5in")

# axs[0].imshow($(stims[1]), vmin=0, vmax=1)#, colorbar="r", colorbar_kw=dict(space=0))
# axs[0].text(6, 6, "median", color="w")
# axs[0].format(title="Epoch 2020-02-04")
# axs[1].imshow($(stims[2]), vmin=0, vmax=1)#, colorbar="r", vmin=0, colorbar_kw=dict(space=0))
# axs[1].text(6, 6, "median", color="w")
# axs[1].format(title="Epoch 2020-11-21")
# m= axs[2].imshow($(stims[3]), vmin=0, vmax=1)#, colorbar="r", vmin=0, colorbar_kw=dict(space=0))
# axs[2].text(6, 6, "annular PCA(2)", color="w")
# axs[2].format(title="Epoch 2020-11-28")
# fig.colorbar(m, loc="r", label="STIM probability")

# axs.format(
#     xticks=$tick_locs,
#     xticklabels=$tick_labs,
#     yticks=$tick_locs,
#     yticklabels=$tick_labs,
#     xlabel="x [AU]",
#     ylabel="y [AU]",
# )

# fig.save($(figuredir("stim.pdf")))

fig, axs = pro.subplots(ncols=3, wspace="1em", figwidth="7.5in")

axs[0].imshow($(stims[1]), vmin=0.5, vmax=1)#, colorbar="r", colorbar_kw=dict(space=0))
axs[0].text(6, 6, "median", color="w")
axs[0].format(title="Epoch 2020-02-04")
axs[1].imshow($(stims[2]), vmin=0.5, vmax=1)#, colorbar="r", vmin=0, colorbar_kw=dict(space=0))
axs[1].text(6, 6, "median", color="w")
axs[1].format(title="Epoch 2020-11-21")
m= axs[2].imshow($(stims[3]), vmin=0.5, vmax=1)#, colorbar="r", vmin=0, colorbar_kw=dict(space=0))
axs[2].text(6, 6, "annular PCA(2)", color="w")
axs[2].format(title="Epoch 2020-11-28")
fig.colorbar(m, loc="r", label="STIM probability")

axs.format(
    xticks=$tick_locs,
    xticklabels=$tick_labs,
    yticks=$tick_locs,
    yticklabels=$tick_labs,
    xlabel="x [AU]",
    ylabel="y [AU]",
)

fig.save($(figuredir("stim_threshold.pdf")))
pro.close(fig)
"""

max_stars_repo_name: mileslucas/sirius-b

max_issues_repo_name: mileslucas/sirius-b

max_forks_repo_name: mileslucas/sirius-b


text:
using NeXLCore
using Test

@testset "Staging" begin
    stg_coords = Dict(:X=>10.0, :Y=>5.0)
    img_coords = Dict(:X=>0.027777, :Y=>-0.016666 )
    rot = deg2rad(50.0)
    i2s = image2stage(NeXLCore.DefaultStageMapping, stg_coords, img_coords, rot)
    s2i = stage2image(NeXLCore.DefaultStageMapping, stg_coords, i2s, rot)
    @test isapprox(s2i[:X], img_coords[:X], atol=1.0e-6)
    @test isapprox(s2i[:Y], img_coords[:Y], atol=1.0e-6)

    v1 = [ -0.2605, 22.8789, 25.695 ]
    v2 = [ -13.7883, 21.1245, 25.662 ]
    v3 = [-4.3441, -0.1964, 25.820 ]
    (ϕ, θ) = compute_tilt(v1, v2, v3)
    @test isapprox(ϕ, 0.38933, atol=0.00001)
    @test isapprox(θ, 118.2447, atol=0.0001)

    v1 = [ 9.0, 0.0, 0.0 ]
    v2 = [ -8.0, 0.0, 0.0 ]
    v3 = [ 0.0, -10.0*cos(deg2rad(1.0)), 10.0*sin(deg2rad(1.0)) ]
    (ϕ, θ) = compute_tilt(v1, v2, v3)
    @test isapprox(ϕ, 1.0, atol=0.00001)
    @test isapprox(θ, 90.0, atol=0.0001)
end

ext: jl


text:
function create_amsg!(sender,receiver,template,model)
    mid = next_amid!(model)
    msg = deepcopy(template)
    msg.id =mid
    msg.sid = sender
    msg.rid = receiver
    msg.ticks = model.ticks
    return msg
end

function send_msg!(receiver::Int64,msg::AGMessage,model)
    rag = getindex(model,receiver)
    #TODO get delay of link in ticks
    g = rag.ctl_graph
    lv = to_local_vertex(g,msg.sid,:aid)
    lva = to_local_vertex(g,rag.id,:aid)
    #need index of nbs
    nbs = neighbors(g,lva)
    i = first(indexin(lv,nbs))
    # log_info(model.ticks,"From $(msg.sid) to  $(receiver) ==> $(msg) ~~~> $(size(rag.msgs_links,1)) --- $i")
    push!(rag.msgs_links[size(rag.msgs_links,1),i],msg)

end

function send_to_nbs!(msg_template::AGMessage,a::Agent,model)::Array{Int64}
    cg = a.ctl_graph
    nbs = neighbors(cg,to_local_vertex(cg,a.id,:aid))
    gid_nbs = [cg[v,:aid] for v in nbs]
    gid_nbs = [gid for gid in gid_nbs if ~(gid in msg_template.body[:trace]) ]
    random_nbs = rand(Binomial(1,model.prob_random_walks),length(gid_nbs))
    random_nbs = [rnb for rnb in random_nbs .* gid_nbs if rnb > 0]
    msgs_sent::Array{Int64} = []

    #disable msgs
    for nb in random_nbs
        if ~(nb in msg_template.body[:trace])
            body = deepcopy(msg_template.body)
            fw_msg = create_amsg!(a.id,nb,msg_template,model)
            send_msg!(nb,fw_msg,model)
            push!(msgs_sent,nb)
        end
    end

    return msgs_sent

end

"""
    Processing for MATCH_PATH msg
"""
function do_match!(msg::AGMessage,a::Agent,model)
    # log_info(model.ticks,a.id,"do_match! -> msg : $(msg)")

    query = msg.body[:query]
    new_path = msg.body[:path]# (tick,confidence,score,path)
    ces = get_controlled_assets(a.id,model)

    #elements this agent controls that are in path
    ces_in_path  = intersect(ces,last(new_path))

    #due diligence: check if path is invalid according to local Knowledge to avoid using/propagating rubbish.
    #however new check required when path query is done locally, as nodes
    # might have dropped since the time the path is received here.
    invalid_path = is_invalid_path(new_path,ces_in_path,model)
    # log_info(model.ticks,a.id,25,"invalidity check: $new_path ==> $invalid_path")
    if !invalid_path

        trace_bk = msg.body[:trace_bk]
        if length(trace_bk) > 1
            save_ctl_path!(a,reverse(trace_bk[1:end-1]))
        end
        # is this agent the original requester of path?
        if first(msg.body[:trace]) == a.id
            for ce in ces_in_path
                spath = last(new_path)[first(indexin(ce,last(new_path))):end]
                # log_info(model.ticks,a.id," procs match for $ce ===> spath: $(spath)")
                #only deals with the exact path, e.g. [7,3,1], not [3,1].
                for i=1:1#length(spath)-1
                    epaths = []
                    if haskey(a.paths,(spath[i],last(spath)))
                        epaths = a.paths[(spath[i],last(spath))]

                        # log_info(model.ticks,a.id," BEFORE: epaths: $epaths")
                        push!(epaths,new_path)
                        # log_info(model.ticks,a.id," AFTER: epaths: $epaths")

                        #sort by score, reverse = false
                        sort!(epaths,lt=isless_paths)

                        #Make sure only model.max_cache_paths are stored
                        if length(epaths) > model.max_cache_paths
                            pop!(epaths)
                        end

                        #TODO cases of older paths
                    else
                        epaths = [new_path]
                    end
                    a.paths[(spath[i],last(spath))] = epaths
                end
            end
            mark_reprocess_of_msg!(a,msg)
            clear_pending_query!(a,query)
            # log_info(model.ticks,a.id,"new path ENDING do_match: => $(a.paths) ++++++++===> $(a.pending)")
        else # This agent is not the original requester of the path
            #continue back propagation of msg
            msg.body[:trace_bk] = trace_bk[1:end-1]
            msg.rid = trace_bk[end-1]
            msg.sid = a.id
            send_msg!(trace_bk[end-1],msg,model)
        end
    else # the path is invalid
        #TODO
    end
end

"""
    Prepare MATCH reply
"""
function do_match!(found_path::Tuple{Int64,Float64,Float64,Array{Int64}},msg::AGMessage,a::Agent,model)
    # log_info(model.ticks,a.id,"*do_match! -> msg : $(msg) -> found: $(found_path)")
    query = msg.body[:query]
    trace = deepcopy(msg.body[:trace])
    append!(trace,a.id)
    trace_bk = deepcopy(trace)
    save_ctl_path!(a,reverse(trace_bk[1:end-1]))
    of_mid = msg.body[:of_mid]

    nbody = Dict(:of_mid=>of_mid,:query=>query,:trace=>trace,:trace_bk=>trace_bk[1:end-1],:path=>found_path)
    # log_info("[$(model.ticks)]{$(a.id)} =-do_match! -> receiver : $(trace_bk)")
    rpy_msg = AGMessage(next_amid!(model),model.ticks,a.id,trace_bk[end-1],MATCH_PATH,nbody)
    send_msg!(trace_bk[end-1],rpy_msg,model)
end

function process_msg!(a::Agent,msg::AGMessage,model)

    @match msg.reason begin
        AG_Protocol(1) =>
                        #Query coming from other agent
                        do_query!(msg,a,model)

        AG_Protocol(2) =>
                        do_match!(msg,a,model)

        AG_Protocol(3) =>
                        do_new_nb!(msg,a,model)
        AG_Protocol(4) =>
                        do_ne_down(a,msg,model)

        # AG_Protocol(5) =>
        #                 do_update_local_graph!(a,msg,model)

        _ => begin
                log_info("[$(model.ticks)]($(a.id)) -> match default!!")
            end
    end

    state = get_state(a)
    state.in_ag_msg+=1.0
    set_state!(a,state)
end

function do_new_nb!(msg::AGMessage,a::Agent,model)

    new_nbs = msg.body[:new_nbs]
    lg = a.ctl_graph
    ls = to_local_vertex(lg,a.id)
    for nb in new_nbs
        ld = to_local_vertex(lg,nb)
        if !has_edge(lg,ls,ld)
            add_edge!(lg,ls,ld)
        end
        if !has_edge(lg,a.id,nb)
            add_edge!(model.ctl_graph,a.id,nb)
        end
    end

end

 """
    It process a msg from a neighbor controller notifying a sne is down
 """
function do_ne_down(a::Agent,msg::AGMessage,model)
    remove_drop_sne!(a,msg.body[:dpid],model.ticks)
    if !isempty(msg.body[:ctl_trace])
        propagate_drop_sne!(a,msg.body[:dpid],msg.body[:ctl_trace],model)
    end
end

"""
 It simulates a NE-controller link down
 In reality this is the routine that checks heartbeats
 from controlled NEs.
"""
function controlled_sne_down!(a::Agent,dpn_id::Int,model)

    set_control_agent!(dpn_id,a.id*-1,model)
    # init_agent!(a,model)
    if isempty(get_live_snes(a,model))
        s = get_state(a)
        s.up = false
        set_state!(a,s)
    end

    lvb = to_local_vertex(a.base_ntw_graph,dpn_id)
    lvc = to_local_vertex(a.ntw_graph,dpn_id)
    a.base_ntw_graph = soft_remove_vertex(a.base_ntw_graph,lvb)
    a.ntw_graph = soft_remove_vertex(a.ntw_graph,lvc)

    #TODO implement when a control agent is down too
    # do_drop!(msg,a,model)
end

"""
 It simulates a NE-controller link up
 In reality this is the routine that checks heartbeats
 from controlled NEs.
"""

function controlled_sne_up!(a::Agent,rjn_id::Int,live_nbs::Array{Int64},model::ABM)

    set_control_agent!(rjn_id,a.id,model)

    if !is_up(a)
        set_up!(a)
    end

    lvb = to_local_vertex(a.base_ntw_graph,rjn_id)
    lvc = to_local_vertex(a.ntw_graph,rjn_id)
    # add_vertex!(a.base_ntw_graph)
    # add_vertex!(a.ntw_graph)

    # set_prop!(a.base_ntw_graph,nv(a.base_ntw_graph),:eid,rjn_id)
    # set_prop!(a.ntw_graph,nv(a.ntw_graph),:eid,rjn_id)

    #just added vertices, last id
    # lvb = nv(a.base_ntw_graph)
    # lvc = nv(a.ntw_graph)

    a.base_ntw_graph = add_edges_gids(a.base_ntw_graph,lvb,live_nbs,:eid)
    a.ntw_graph = add_edges_gids(a.ntw_graph,lvc,live_nbs,:eid)
end

function get_state(a::Agent)::State
    return last(a.state_trj)
end
function set_state!(a::Agent,new_state::ControlAgentState)
    pop!(a.state_trj)
    push!(a.state_trj,new_state)
end

# function get_throughput(a::Agent,model)
#     a.msgs_in
# end

function get_state_trj(a::Agent)::Vector{State}
    return a.state_trj
end

function get_throughput_up(a::Agent,model)
    v_msg_in = [ s.in_ag_msg for s in a.state_trj ]

    # log_info("[$(model.ticks)] ($(a.id)) msg in  ==> $v_msg_in")
    v_up = [ s.up for s in a.state_trj ]
    v_tpt = get_throughput(v_msg_in,model.:interval_tpt)
    return [ v_up[i] ? v_tpt[i] : 0.0   for i=1:length(v_tpt)]
end

"""
 Check if path is invalid according to local knowledge. It does not mean is valid as this agent might not have all info to check validity.
 Only flags invalid if 100% sure.
 cnes:controlled nes by this agent (a)
"""
function is_invalid_path(path,cnes,model)
    result = false

    for cne in cnes
        if last(last(path)) != cne
            sne = getindex(model,cne)
            cne_i = first(indexin(cne,last(path)))
            ports_path = filter(p->p[2]=="s$(last(path)[cne_i+1])",get_port_edge_list(sne))
            result = isempty(ports_path) ? true : result
        end
    end
    return result
end

function save_ctl_path!(a::Agent,path::Array{Int64})
    already_in = false
    for p in a.ctl_paths
        if isless(path,p) || p == path#is path a strict subset or equal to p ?
            already_in = true
        end
    end
    if !already_in
        push!(a.ctl_paths,path)
    end
end

function propagate_drop_sne!(a::Agent,dpid::Int64,ctl_trace::Array{Int64},model)
    sid = a.id
    rid = first(ctl_trace)
    rem_path = length(ctl_trace) > 1 ? ctl_trace[2:end] : []
    dmsg = AGMessage(next_amid!(model),model.ticks,sid,rid,NE_DOWN,Dict(:dpid=>dpid,:ctl_trace=>rem_path))
    send_msg!(rid,dmsg,model)
end

"""
    It removes a dropped sne node from local paths and graph
"""
function remove_drop_sne!(a::Agent,dpid::Int64,drop_time::Int64)
    # log_info(drop_time,a.id,"Existing paths: $(a.paths)")
    new_paths_dict = Dict()
    #delete pre-computed paths containing dropping node
    for path_k in keys(a.paths)
        v_paths = a.paths[path_k]
        new_paths = []
        for path in v_paths
            if !(dpid in last(path))
                push!(new_paths,path)
            end
        end
        if !isempty(new_paths)
            new_paths_dict[path_k] = new_paths
        end
    end
    #update agent pre-computed paths
    a.paths = new_paths_dict

    #update graphs used by the control agent accordingly

   #is dpid in base ntw graph?
   lvb = to_local_vertex(a.base_ntw_graph,dpid)
   lvc = to_local_vertex(a.ntw_graph,dpid)

   log_info(drop_time,a.id," Removing dpid: $dpid from local base: $lvb --- local curr: $lvc")

   if lvb != 0
    a.base_ntw_graph = soft_remove_vertex(a.base_ntw_graph,lvb)
   end
   if lvc != 0
        a.ntw_graph = soft_remove_vertex(a.ntw_graph,lvc)
   end

   state = get_state(a)
   new_ap = Dict()
   for p in pairs(state.active_paths)
        if !(dpid in last(p))
          new_ap[first(p)] = last(p)
        end
   end

   state.active_paths = new_ap
   set_state!(a,state)
end

"""
It deals with prediction of unavailability (for a given time window) of a set of NEs under control.

"""
function do_update_flows_from_changes!(a::Agent,ntw_changes::Vector{Int64},model::ABM)
    #TODO operation for other than centralised agent
    if get_state(a).up

        log_info(model.ticks,a.id,"Pred_Down: ntw_changes!!!!!! $(ntw_changes)")

        joining_nodes = filter(x->x>0,ntw_changes)
        dropping_nodes = -1 * filter(x->x<0,ntw_changes)

        query_graph = deepcopy(a.ntw_graph)

        #add nodes
        for jng_id in joining_nodes
            query_graph = soft_remove_vertex(query_graph,jng_id)

            base_g = a.base_ntw_graph

            lv = to_local_vertex(base_g,jng_id)
            nbs = neighbors(base_g,lv)

            #TODO check if nbs is up?

            for nb_id in nbs
                add_edge!(query_graph,nb_id,jng_id)
                add_edge!(query_graph,jng_id,nb_id)
            end
        end

        #remove nodes
        for dpn_id in dropping_nodes
            query_graph = soft_remove_vertex(query_graph,dpn_id)
        end

        query_time = model.ticks

        #Only trigger queries if maintenance policy is not predictive (where routes/flows come from the optimisation algorithm)
        queries = a.maintenance.policy != PredictiveM  ? model.ntw_services : []

        for query in queries
            query_paths = Dict{Tuple{Int64,Int64},Array{Tuple{Int64,Float64,Float64,Array{Int64}}}}()

            path = do_query(query_time,query,query_graph,query_paths)

            if isempty(path)
                 log_info(model.ticks,a.id,"No path found => $query")
            #     # if !haskey(a.previous_queries,query)
            #     #     sent_to = query_nbs!(a,msg,jg,query,trace,model)
            #     #     a.previous_queries[query] = (model.ticks,sent_to)
            #     # end
            else
            #     # do_match!(path,msg,a,model)
            #     # clear_pending_query!(a,query)
                do_update_flows_from_path!(a,last(path),model)
            end

        end

    end
end

"""
Update flows of the snes controlled by the agent a and for the path given,
"""
function do_update_flows_from_path!(a::Agent,path::Array{Int64,1},model::ABM)
    msg = OFMessage(-1, model.ticks,-1,0,OFPR_ADD_FLOW,[])
    # log_info(model.ticks,a.id,"active path: $path")
    install_flow!(a,path,model,msg)
    if length(path) > 1
        k = (first(path),last(path))
        # log_info(model.ticks,a.id,"key: $(k) ==> Ag Path: $spath")
        get_state(a).active_paths[k] = path
    end
end


text:
# # Summary
#
# * Challenges of parallel computing
#     * Order of execution
#         * execution of out order of Possibility
#         * simultaneous access and mutation
#     * Data access and movement
#     * Code access and movement
#     * Appropriately matching the parallelism strategy to your machine capabilities
#     * Appropriately matching the parallelism strategy with the problem at hand
#
# * Parallelism strategies
#     * SIMD
#     * Multithreading
#     * Tasks
#     * Multi-process
#         * Shared memory
#         * Distributed memory
#     * GPU programming
#
#
# ## Why so many kinds of parallelism?
#
# * Not all problems are created equal
# * Not all computing machines are created equal
# * We want to maximize comuting while minimizing overhead
#     * Chosen solution will depend upon the amount of computing in each inner loop
#       and the amount of syncronization that is required between loops.

max_stars_repo_name: fercarozzi/JuliaAcademyMaterials

max_issues_repo_name: fercarozzi/JuliaAcademyMaterials

max_forks_repo_name: fercarozzi/JuliaAcademyMaterials


text:
"""
	velo(cmd0::String="", arg1=nothing; kwargs...)

Plot velocity vectors, crosses, and wedges.

Full option list at [`velo`]($(GMTdoc)velo.html)

```julia
    velo(mat2ds([0. -8 0 0 4 6 0.5; -8 5 3 3 0 0 0.5], ["4x6", "3x3"]), pen=(0.6,:red), fill_wedges=:green, outlines=true, Se="0.2/0.39/18", arrow="0.3c+p1p+e+gred", region=(-15,10,-10,10), show=1)
```
"""
function velo(cmd0::String="", arg1=nothing; first=true, kwargs...)

    proggy = (IamModern[1]) ? "velo "  : "psvelo "
	length(kwargs) == 0 && return monolitic(proggy, cmd0, arg1)

	d, K, O = init_module(first, kwargs...)		# Also checks if the user wants ONLY the HELP mode

	cmd, opt_B, opt_J, opt_R = parse_BJR(d, "", "", O, " -JX12cd/0d")
	cmd, = parse_common_opts(d, cmd, [:UVXY :di :e :p :t :params], first)

	if ((val = find_in_dict(d, [:A :arrow])[1]) !== nothing)
		cmd = (isa(val, String)) ? cmd * " -A" * val : cmd * " -A" * vector_attrib(val)
	end

	cmd = add_opt_fill(cmd, d, [:E :fill_wedges], 'E')
	cmd = add_opt_fill(cmd, d, [:G :fill_symbols], 'G')
	cmd = parse_these_opts(cmd, d, [[:D :sigma_scale], [:L :outlines], [:N :no_clip :noclip]])

	if     (haskey(d, :Se) || haskey(d, :vel_NE))       symbs = [:Se :vel_NE]
	elseif (haskey(d, :Sn) || haskey(d, :barscale))     symbs = [:Sn :barscale]
	elseif (haskey(d, :Sr) || haskey(d, :vel_rotated))  symbs = [:Sr :vel_rotated]
	elseif (haskey(d, :Sw) || haskey(d, :wedges))       symbs = [:Sw :wedges]
	elseif (haskey(d, :Sx) || haskey(d, :cross_scale))  symbs = [:Sx :cross_scale]
	elseif (show_kwargs[1])  symbs = [:Se :vel_NE :Sn :barscale :Sr :vel_rotated :Sw :wedges :Sx :cross_scale]
	else  error("Must select one convention (S options. Run gmthelp(velo) to learn about them)")
	end
	cmd = add_opt(d, cmd, string(symbs[1]), symbs)
	cmd *= opt_pen(d, 'W', [:W :pen])

	# If file name sent in, read it and compute a tight -R if it was not provided
	cmd, arg1, opt_R, = read_data(d, cmd0, cmd, arg1, opt_R)

	finish_PS_module(d, proggy * cmd, "", K, O, true, arg1)
end

# ---------------------------------------------------------------------------------------------------
velo!(cmd0::String="", arg1=nothing; kw...) = velo(cmd0, arg1; first=false, kw...)
velo(arg1,  cmd0::String=""; kw...) = velo(cmd0, arg1; first=true, kw...)
velo!(arg1, cmd0::String=""; kw...) = velo(cmd0, arg1; first=false, kw...)

const psvelo  = velo 			# Alias
const psvelo! = velo!			# Alias


text:
#!/usr/bin/julia

using DelimitedFiles

inputs = readdlm("2019/02/input.txt", ',', Int)

inputs[2] = 12
inputs[3] = 2
println(inputs)

for i in eachindex(inputs)
    if i % 4 == 1
        op, i1, i2, i3 = inputs[i:i+3]
        if op == 1
            inputs[i3+1] = inputs[i1+1] + inputs[i2+1]
        elseif op == 2
            inputs[i3+1] = inputs[i1+1] * inputs[i2+1]
        elseif op == 99
            break
        end
    end
end

println(inputs)

max_stars_repo_name: chaserobertson/advent

max_issues_repo_name: chaserobertson/advent

max_forks_repo_name: chaserobertson/advent


text:
using MathematicaParser
using ParserCombinator
using Test

################################################################################
#                                                                              #
#                                Unit Testing                                  #
#                                                                              #
################################################################################

@testset "Number conversion" begin
    @test parsemathematica("1")[1]         ==  1
    @test parsemathematica("-1")[1]        == -1
    @test parsemathematica("1.0000")[1]    ==  1
    @test parsemathematica("1.0000")[1]    !== 1.0
    @test parsemathematica("1.000001")[1]  ==  1.000001
    @test parsemathematica("1.2")[1]       ==  1.2
    @test parsemathematica("1.2e3")[1]     ==  1200
    @test parsemathematica("-1.2")[1]      == -1.2
end

# Testing variable name patterns
vartestdict = Dict(
    "x"         => :x,
    "x01"       => :x01,
    "xAx01xx"   => :xAx01xx,
    "x_1"       => :x_1,
    "x_1 "      => :x_1,
    " x_1 "     => :x_1,
    " x_1 "     => :x_1,
    "Xx01"      => :Xx01,
    "A"         => :A,
)
@testset "Variable Names" begin
    for (t,res) in vartestdict
        @test parsemathematica(t)[1] == res
    end
end

# Testing argument lists
argstestdict = Dict(
    "[x]"               => MathematicaParser.Args([:x]),
    "[a,b,c,z]"         => MathematicaParser.Args([:a, :b, :c, :z]),
    "[a, b,   c  ,z ]"  => MathematicaParser.Args([:a, :b, :c, :z]),
)
@testset "Argument List - No subexpressions" begin
    for (t,res) in argstestdict
        # Note Arguments on their own are not correct syntax
        # Need to import the args ParserCombinator type to parse args alone
        # @info parse_one(t, MathematicaParser.args)[1]
        @test parse_one(t, MathematicaParser.args)[1] == res
    end
end

# Model Function Test (no sub expressions)
funtestdict = Dict(
    "Log[x]"                        => Fun(:Log, [Args([:x])]),
    "Exp[x]"                        => Fun(:Exp, [Args([:x])]),
    "Log2[x]"                       => Fun(:Log2, [Args([:x])]),
    "ArcTan[x,y]"                   => Fun(:ArcTan, [Args([:x, :y])]),
    "Hypergeometric2F1[a,b,c,z]"    => Fun(:Hypergeometric2F1, [Args([:a, :b, :c, :z])]),
    "Derivative[n][f][x]"           => Fun(:Derivative, [Args([:n]), Args([:f]), Args([:x])]),
)
@testset "Function Parsing - No subexpressions" begin
    for (t,res) in funtestdict
        @test parsemathematica(t)[1] == res
    end
end

# Order of operations
@testset "Order of Operations" begin
    @test parsemathematica("a^b")[1] ==  Pow(:a, :b)
    @test parsemathematica("a^b^c")[1] ==  Pow(:a, Pow(:b, :c))
    @test parsemathematica("2+3*4")[1] == Sum([2, Prd([3, 4])])
    @test parsemathematica("4*-6*(3 * 7 + 5) ")[1] == Prd([4, -6, Sum([Prd([3, 7]), 5])])
    @test parsemathematica("7^-3*2")[1] == Prd([Pow(7, -3), 2])
end

parsemathematica("Gamma[x,2*x]/x^3")
parsemathematica("Gamma[x, 2*x]/x^3")

afuntestdict = [
    "f[x]"              => Fun(:f, Args[Args([:x])]),
    "f'[x]"             => Fun(:Derivative, Args[Args(1), Args(:f), Args([:x])]),
    "f''[x]"            => Fun(:Derivative, Args[Args(2), Args(:f), Args([:x])]),
    "f'''[x]"           => Fun(:Derivative, Args[Args(3), Args(:f), Args([:x])]),
    "f'[x]/f[x]"        => Prd([Fun(:Derivative, Args[Args(1), Args(:f), Args([:x])]), Inv(Fun(:f, Args[Args([:x])]))]),
    "Sin[x] + f'[x]"    => Sum([Fun(:Sin, Args[Args([:x])]), Fun(:Derivative, Args[Args(1), Args(:f), Args([:x])])]),
    "f[g'[x]]"          => Fun(:f, Args[Args([Fun(:Derivative, Args[Args(1), Args(:g), Args([:x])])])]),
]
@testset "Derivative shorthand" begin
    for (t,r) in afuntestdict
        @test parsemathematica(t)[1] == r
    end
end

convtestdict = [
    "1*2 + 3"           => :(1 * 2 + 3),
    "3^7^8"             => :(3 ^ (7 ^ 8)),
    "3/8"               => :(3 * (1 / 8)),
    "-5 +9"             => :(-5 + 9),
    "Sin[x]"            => :(Sin(x)),
    "-Sin[x]"           => :(-1 * Sin(x)),
    "Sin[a,b][y][z]"    => :(Sin(a, b, y, z)),
    "4*-6*(3 ^ 7 + 5) + Sin[x]" => :(4 * -6 * (3 ^ 7 + 5) + Sin(x)),
]
@testset "Conversion to Expressions" begin
    for (t,r) in convtestdict
        # println( parsemathematica(t)[1])
        # println(toexpr(parsemathematica(t)[1]))
        @test toexpr(parsemathematica(t)[1]) == r
    end
end


text:
"""
    mixableRaster(raster::Raster)::RasterMix

Convert a raster into a form that is suitable for combining with other rasters.
"""
function mixableRaster(raster::Raster)::RasterMix
    colors = copy(raster[1:3, :, :])
    colors[1, :, :] .*= raster[4, :, :]
    colors[2, :, :] .*= raster[4, :, :]
    colors[3, :, :] .*= raster[4, :, :]
    return (colors, 1 .- raster[4, :, :], raster[4, :, :])
end

"""
    mixRasters(
        (c1,ia1,w1)::RasterMix,
        (c2,ia2,w2)::RasterMix)::RasterMix

Combine 2 rasters into mixable form.
"""
function mixRasters((c1, ia1, w1)::RasterMix, (c2, ia2, w2)::RasterMix)::RasterMix
    return (c1 .+ c2, ia1 .* ia2, w1 .+ w2)
end

"""
    mixedRaster((r,ia,w)::RasterMix)::Raster

Convert several (possibly many) mixed rasters back into normal raster.
"""
function mixedRaster((r, ia, ws)::RasterMix)::Raster
    w = copy(reshape(ws, (1, size(ws)...)))
    w[w.==0] .= 1
    cat(
        dims = 1,
        r[1:1, :, :] ./ w,
        r[2:2, :, :] ./ w,
        r[3:3, :, :] ./ w,
        1 .- reshape(ia, (1, size(ia)...)),
    )
end


text:
@show parse(Int, "123459")
@show parse(Int, "0xabcf123")
@show parse(Int, "0o7651")
@show parse(Int, "0b101011001")

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
"""
    AnovaValue <: AnovaEffect

A set of values for an Anova item for which a mean square is not required.

`name` - the name of this value

`ss` - sum of squares

`df` - degrees of freedom
"""
struct AnovaValue <: AnovaEffect
    name::String
    ss::Float64
    df::Float64
end

AnovaValue(ss, df) = AnovaValue("", ss, df)
AnovaValue(name, x::AnovaEffect) = AnovaValue(name, x.ss, x.df)

import Base.isapprox
isapprox(x::AnovaValue, y::AnovaValue) =
    x.name == y.name &&
    x.ss ≈ y.ss &&
    x.df == y.df


text:
println(replace(raw"""
  xxxxx
     x            x
     x            x   x
     x   x  x   x
     x   x  x   x   x    xxx
     x   x  x   x   x   x   x
 x  x   x  x   x   x   x   xx
  xx     xx    x   x    x x xx
 """, "x" => "_/"))


text:
export Presence

"""
A [`User`](@ref)'s presence.
More details [here](https://discordapp.com/developers/docs/topics/gateway#presence-update).
"""
struct Presence
    user::User
    roles::Optional{Vector{Snowflake}}
    game::Nullable{Activity}
    guild_id::Optional{Snowflake}
    status::String
    activities::Vector{Activity}
end
@boilerplate Presence :constructors :docs :lower :merge :mock


text:
import Dolo

path = Dolo.pkg_path

@testset "testing perfect_foresight" begin

    fn = joinpath(path, "examples", "models", "rbc.yaml")

    model = Dolo.yaml_import(fn)

    # we must define the series for exogenous shocks
    n_e = length(model.symbols[:exogenous])

    T_e = 5
    exo = zeros(T_e, n_e)
    exo[1,:] = [0.00]  # this is used to determine initial steady-state of the model
    exo[2,:] = [0.01]
    exo[3,:] = [0.02]
    exo[4,:] = [0.03]
    exo[5,:] = [0.04]  # this is used to determine final steady-state

    @time df = Dolo.perfect_foresight(model, exo, T=20)

    @test true

end

# # Second method: Define a dictionary with shocks
# # Compare the results with the previous method

# path = Dolo.pkg_path
# fn = joinpath(path, "examples", "models", "rbc.yaml");

# model = Dolo.yaml_import(fn);

# n_e = length(model.symbols[:exogenous])
# T_e = 5
# exo = zeros(T_e, n_e)
# exo[1,:] = [0.00]  # this is used to determine initial steady-state of the model
# exo[2,:] = [0.01]
# exo[3,:] = [0.02]
# exo[4,:] = [0.03]
# exo[5,:] = [0.04]
# @time df_1matrix = Dolo.perfect_foresight(model, exo, T=20)

# exo = Dict(:z =>  [0, 0.01, 0.02, 0.03, 0.04])
# @time df_1dict = Dolo.perfect_foresight(model, exo, T=20)

# if df_1dict == df_1matrix
#   println("df is the same")
# end

# # Two shocks
# fn = joinpath(path, "examples", "models", "rbc_catastrophe.yaml");

# model = Dolo.yaml_import(fn)

# n_e = length(model.symbols[:exogenous])
# T_e = 5
# exo = zeros(T_e, n_e)
# exo[1,:] = [0.00, 0.00]  # this is used to determine initial steady-state of the model
# exo[2,:] = [0.01, 0.00]
# exo[3,:] = [0.02, 0.00]
# exo[4,:] = [0.03, 0.00]
# exo[5,:] = [0.04, 0.00]
# @time  df_2matrix = Dolo.perfect_foresight(model, exo, T=20)

# exo = Dict(:z =>  [0.00, 0.01, 0.02, 0.03, 0.04], :xi =>  [0.00, 0.00, 0.00, 0.00, 0.00])
# @time  df_2dict_a = Dolo.perfect_foresight(model, exo, T=20)

# if df_2matrix == df_2dict_a
#   println("df is the same")
# end

# # or

# exo = Dict(:z =>  [0.00, 0.01, 0.02, 0.03, 0.04])
# @time df_2dict_b = Dolo.perfect_foresight(model, exo, T=20)

# if df_2matrix == df_2dict_b
#   println("df is the same")
# end

# #or order changed
# exo = Dict( :xi =>  [0.00, 0.00, 0.00, 0.00, 0.00], :z =>  [0.00, 0.01, 0.02, 0.03, 0.04])
# @time df_2dict_c = Dolo.perfect_foresight(model, exo, T=20)

# if df_2matrix == df_2dict_c
#   println("df is the same")
# end

# #or not defined as Float64
# exo = Dict( :xi =>  [0, 0, 0, 0, 0], :z =>  [0.00, 0.01, 0.02, 0.03, 0.04])

# @time df_2dict_d = Dolo.perfect_foresight(model, exo, T=20)
# if df_2matrix == df_2dict_d
#   println("df is the same")
# end


text:
mutable struct MyCallback{L,F,F2}
  losses::L
  cb::F
  ecb::F2
  nepochs::Int
  nsamples::Int
  iter::Int
  epoch::Int
  print_iter_loss::Bool
  print_epoch_mean_loss::Bool
  plot_iter::Bool
  plot_epoch::Bool
  MyCallback(losses, cb, ecb, nepochs, nsamples, iter, epoch, print_iter_loss, print_epoch_mean_loss, plot_iter, plot_epoch) =
    new{typeof(losses), typeof(cb), typeof(ecb)}(losses, cb, ecb, nepochs, nsamples, iter, epoch, print_iter_loss, print_epoch_mean_loss, plot_iter, plot_epoch)
end

MyCallback(T::DataType=Float32;
              cb=(args...;kwargs...)->false, ecb=(args...;kwargs...)->false,
              nepochs=1, nsamples=1,
              print_iter_loss=false, plot_every_iter=false,
              print_epoch_mean_loss=true, plot_every_epoch=true) =
  MyCallback(T[], cb, ecb, nepochs, nsamples, 0, 0, print_iter_loss, print_epoch_mean_loss, plot_every_iter, plot_every_epoch)

function (mcb::MyCallback)(p,l, args...; kwargs...)
  cbout = invoke_sample_cb!(mcb, p,l, args...; kwargs...)
  mcb.iter % mcb.nsamples == 0 && invoke_epoch_cb!(mcb, p, l, args...; kwargs...)
  return cbout
end

function reset!(mcb::MyCallback)
  mcb.losses = Vector{eltype(mcb.losses)}(undef, 0)
  mcb.iter = 0
  nothing
end

function invoke_sample_cb!(mcb::MyCallback, p,l,args...; kwargs...)
  mcb.print_iter_loss == true && println(l)
  mcb.plot_iter == true && mcb.cb(p,l,args...; kwargs...)
  push!(mcb.losses, l)
  mcb.iter += 1
  return false
end

function invoke_epoch_cb!(mcb::MyCallback, p, l, args...; kwargs...)
  mcb.epoch += 1
  mcb.print_epoch_mean_loss == true && println("Epoch $(mcb.epoch) mean loss: $(mean(mcb.losses))")
  mcb.plot_epoch == true && mcb.ecb(p,l,args...; kwargs...)
  reset!(mcb)
  return false
end


text:
module CompositeStructures

using ..Traits: getcontent, dissolve, rawtype, efficientoperations

import ..Traits: contentnames, dissolve

export CompositeTuple, CompositeNTuple, CompositeVector, CompositeDict
export NamedContainer

"""
    CompositeTuple{T<:Tuple}

A composite tuple can be considered as a tuple that is implemented by including an ordinary `Tuple` as its data attribute.
"""
abstract type CompositeTuple{T<:Tuple} end
@inline contentnames(::Type{<:CompositeTuple}) = (:contents,)
@inline dissolve(ct::CompositeTuple, ::Val{:contents}, f::Function, args::Tuple, kwargs::NamedTuple) = f(getcontent(ct, :contents), args...; kwargs...)

@inline Base.length(::CompositeTuple{T}) where T = fieldcount(T)
@inline Base.length(::Type{<:CompositeTuple{T}}) where T = fieldcount(T)
@inline Base.eltype(::CompositeTuple{T}) where T = eltype(T)
@inline Base.eltype(::Type{<:CompositeTuple{T}}) where T = eltype(T)
@inline Base.hash(ct::CompositeTuple, h::UInt) = hash(dissolve(ct), h)
@inline Base.:(==)(ct1::CompositeTuple, ct2::CompositeTuple) = ==(efficientoperations, ct1, ct2)
@inline Base.isequal(ct1::CompositeTuple, ct2::CompositeTuple) = isequal(efficientoperations, ct1, ct2)
@inline Base.getindex(ct::CompositeTuple, i::Union{<:Integer, CartesianIndex}) = getcontent(ct, :contents)[i]
@inline Base.getindex(ct::CompositeTuple, inds) = rawtype(typeof(ct))(dissolve(ct, getindex, (inds,))...)
@inline Base.lastindex(ct::CompositeTuple) = lastindex(getcontent(ct, :contents))
@inline Base.iterate(ct::CompositeTuple) = iterate(getcontent(ct, :contents))
@inline Base.iterate(ct::CompositeTuple, state) = iterate(getcontent(ct, :contents), state)
@inline Base.iterate(rv::Iterators.Reverse{<:CompositeTuple}, state=length(rv.itr)) = (state < 1) ? nothing : (rv.itr[state], state-1)
@inline Base.keys(ct::CompositeTuple) = keys(getcontent(ct, :contents))
@inline Base.values(ct::CompositeTuple) = values(getcontent(ct, :contents))
@inline Base.pairs(ct::CompositeTuple) = pairs(getcontent(ct, :contents))
@inline Base.reverse(ct::CompositeTuple) = rawtype(typeof(ct))(dissolve(ct, reverse)...)
@inline Base.Tuple(ct::CompositeTuple) = getcontent(ct, :contents)

"""
    CompositeNTuple{N, T}

A composite ntuple can be considered as a ntuple that is implemented by including an ordinary `NTuple` as its data attribute.

Alias for `CompositeTuple{NTuple{N, T}}`.
"""
const CompositeNTuple{N, T} = CompositeTuple{NTuple{N, T}}

"""
    CompositeVector{T}

A composite vector can be considered as a vector that is implemented by including a concrete subtype of `AbstractVector` as its data attribute.
"""
abstract type CompositeVector{T} <: AbstractVector{T} end
@inline contentnames(::Type{<:CompositeVector}) = (:contents,)
@inline dissolve(cv::CompositeVector, ::Val{:contents}, f::Function, args::Tuple, kwargs::NamedTuple) = f(getcontent(cv, :contents), args...; kwargs...)

@inline Base.size(cv::CompositeVector) = size(getcontent(cv, :contents))
@inline Base.size(cv::CompositeVector, i) = size(getcontent(cv, :contents), i)
@inline Base.length(cv::CompositeVector) = length(getcontent(cv, :contents))
@inline Base.:(==)(cv1::CompositeVector, cv2::CompositeVector) = ==(efficientoperations, cv1, cv2)
@inline Base.isequal(cv1::CompositeVector, cv2::CompositeVector) = isequal(efficientoperations, cv1, cv2)
@inline Base.getindex(cv::CompositeVector, i::Union{<:Integer, CartesianIndex}) = getcontent(cv, :contents)[i]
@inline Base.getindex(cv::CompositeVector, inds) = rawtype(typeof(cv))(dissolve(cv, getindex, (inds,))...)
@inline Base.lastindex(cv::CompositeVector) = lastindex(getcontent(cv, :contents))
@inline Base.setindex!(cv::CompositeVector, value, inds) = (getcontent(cv, :contents)[inds] = value)
@inline Base.push!(cv::CompositeVector, values...) = (push!(getcontent(cv, :contents), values...); cv)
@inline Base.pushfirst!(cv::CompositeVector, values...) = (pushfirst!(getcontent(cv, :contents), values...); cv)
@inline Base.insert!(cv::CompositeVector, index::Integer, value) = (insert!(getcontent(cv, :contents), index, value); cv)
@inline Base.append!(cv::CompositeVector, values) = (append!(getcontent(cv, :contents), values); cv)
@inline Base.prepend!(cv::CompositeVector, values) = (prepend!(getcontent(cv, :contents), values); cv)
@inline Base.splice!(cv::CompositeVector, index::Integer, replacement=Base._default_splice) = splice!(getcontent(cv, :contents), index, replacement)
@inline Base.splice!(cv::CompositeVector, range::UnitRange{<:Integer}, replacement=Base._default_splice) = rawtype(typeof(cv))(dissolve(cv, splice!, (range, replacement))...)
@inline Base.deleteat!(cv::CompositeVector, indices) = (deleteat!(getcontent(cv, :contents), indices); cv)
@inline Base.pop!(cv::CompositeVector) = pop!(getcontent(cv, :contents))
@inline Base.popfirst!(cv::CompositeVector) = popfirst!(getcontent(cv, :contents))
@inline Base.empty!(cv::CompositeVector) = (empty!(getcontent(cv, :contents)); cv)
@inline Base.empty(cv::CompositeVector) = rawtype(typeof(cv))(dissolve(cv, empty)...)
@inline Base.reverse(cv::CompositeVector) =  rawtype(typeof(cv))(dissolve(cv, reverse)...)
@inline Base.similar(cv::CompositeVector, dtype::Type=eltype(cv), dims::Tuple{Vararg{Int}}=size(cv)) = rawtype(typeof(cv))(dissolve(cv, similar, (dtype, dims))...)
@inline Base.iterate(cv::CompositeVector, state=1) = iterate(getcontent(cv, :contents), state)
@inline Base.keys(cv::CompositeVector) = keys(getcontent(cv, :contents))
@inline Base.values(cv::CompositeVector) = values(getcontent(cv, :contents))
@inline Base.pairs(cv::CompositeVector) = pairs(getcontent(cv, :contents))

"""
    CompositeDict{K, V}

A composite dict can be considered as a dict that is implemented by including a concrete subtype of `AbstractDict` as its data attribute.
"""
abstract type CompositeDict{K, V} <: AbstractDict{K, V} end
@inline contentnames(::Type{<:CompositeDict}) = (:contents,)
@inline dissolve(cd::CompositeDict, ::Val{:contents}, f::Function, args::Tuple, kwargs::NamedTuple) = f(getcontent(cd, :contents), args...; kwargs...)

@inline Base.isempty(cd::CompositeDict) = isempty(getcontent(cd, :contents))
@inline Base.length(cd::CompositeDict) = length(getcontent(cd, :contents))
@inline Base.haskey(cd::CompositeDict, key) = haskey(getcontent(cd, :contents), key)
@inline Base.in(p::Pair, cd::CompositeDict, valcmp=(==)) = in(p, getcontent(cd, :contents), valcmp)
@inline Base.:(==)(cd1::CompositeDict, cd2::CompositeDict) = ==(efficientoperations, cd1, cd2)
@inline Base.isequal(cd1::CompositeDict, cd2::CompositeDict) = isequal(efficientoperations, cd1, cd2)
@inline Base.get(cd::CompositeDict, key, default) = get(getcontent(cd, :contents), key, default)
@inline Base.get(f::Base.Callable, cd::CompositeDict, key) = get(f, getcontent(cd, :contents), key)
@inline Base.getkey(cd::CompositeDict, key, default) = getkey(getcontent(cd, :contents), key, default)
@inline Base.getindex(cd::CompositeDict{K, V}, index::K) where {K, V} = getcontent(cd, :contents)[index]
@inline Base.push!(cd::CompositeDict, ps::Pair...) = (push!(getcontent(cd, :contents), ps...); cd)
@inline Base.get!(cd::CompositeDict, key, default) = get!(getcontent(cd, :contents), key, default)
@inline Base.get!(default::Union{Function, Type}, cd::CompositeDict, key) = get!(default, getcontent(cd, :contents), key)
@inline Base.setindex!(cd::CompositeDict{K, V}, value::V, index::K) where {K, V} = (getcontent(cd, :contents)[index] = value)
@inline Base.pop!(cd::CompositeDict) = pop!(getcontent(cd, :contents))
@inline Base.pop!(cd::CompositeDict, key) = pop!(getcontent(cd, :contents), key)
@inline Base.pop!(cd::CompositeDict, key, default) = pop!(getcontent(cd, :contents), key, default)
@inline Base.delete!(cd::CompositeDict, key) = (delete!(getcontent(cd, :contents), key); cd)
@inline Base.empty!(cd::CompositeDict) = (empty!(getcontent(cd, :contents)); cd)
@inline Base.merge(cd::CD, others::CD...) where CD <: CompositeDict = merge!(empty(cd), cd, others...)
@inline Base.merge(combine::Function, cd::CD, others::CD...) where {CD<:CompositeDict} = merge!(combine, empty(cd), cd, others...)
@inline Base.empty(cd::CompositeDict) = rawtype(typeof(cd))(dissolve(cd, empty)...)
@inline Base.iterate(cd::CompositeDict) = iterate(getcontent(cd, :contents))
@inline Base.iterate(cd::CompositeDict, state) = iterate(getcontent(cd, :contents), state)
@inline Base.keys(cd::CompositeDict) = keys(getcontent(cd, :contents))
@inline Base.values(cd::CompositeDict) = values(getcontent(cd, :contents))
@inline Base.pairs(cd::CompositeDict) = pairs(getcontent(cd, :contents))

"""
    NamedContainer{T, Names} = NamedTuple{Names, <:Tuple{Vararg{T}}}

NamedContainer is just a wrapper of Julia NamedTuple, but not a composite type.
"""
const NamedContainer{T, Names} = NamedTuple{Names, <:Tuple{Vararg{T}}}
@inline NamedContainer{Names}(contents::Tuple) where Names = NamedTuple{Names}(contents)

end #module


text:
using BinDeps
using Compat

@BinDeps.setup

function cflags_validator(pkg_name)
    return (name, handle) -> begin
        try
            run(`pkg-config --cflags $(pkg_name)`)
            return true
        catch ErrorException
            return false
        end
    end
end

@linux? (
    begin
        deps = [
            python = library_dependency("python", aliases=["libpython2.7.so", "libpython3.2.so", "libpython3.3.so", "libpython3.4.so", "libpython3.5.so", "libpython3.6.so", "libpython3.7.so"], validate=cflags_validator("python"))
            glib = library_dependency("glib", aliases=["libglib-2.0-0", "libglib-2.0", "libglib-2.0.so.0"], depends=[python], validate=cflags_validator("glib-2.0"))
            lcm = library_dependency("lcm", aliases=["liblcm", "liblcm.1"], depends=[glib])

            provides(AptGet, Dict("python-dev" => python, "libglib2.0-dev" => glib))
        ]
    end
    : begin
        deps = [
            glib = library_dependency("glib", aliases = ["libglib-2.0-0", "libglib-2.0", "libglib-2.0.so.0"])
            lcm = library_dependency("lcm", aliases=["liblcm", "liblcm.1"], depends=[glib])
        ]
    end
    )

prefix = joinpath(BinDeps.depsdir(lcm), "usr")
pkg_config_dirs = AbstractString[]
include_dirs = AbstractString[]
@osx_only begin
    if Pkg.installed("Homebrew") === nothing
        error("Homebrew package not installed, please run Pkg.add(\"Homebrew\")")
    end
    using Homebrew
    provides(Homebrew.HB, "glib", glib, os=:Darwin)
    push!(pkg_config_dirs, joinpath(Homebrew.prefix(), "lib", "pkgconfig"))
    push!(include_dirs, joinpath(Homebrew.prefix(), "include"))
end

provides(Yum,
    Dict("glib" => glib))

provides(Sources,
    URI("https://github.com/lcm-proj/lcm/releases/download/v1.3.1/lcm-1.3.1.zip"),
    lcm)

provides(BuildProcess, Dict(Autotools(libtarget="lcm/liblcm.la", include_dirs=include_dirs, pkg_config_dirs=pkg_config_dirs) => lcm), onload="""
using PyCall
sys = pyimport("sys")
unshift!(PyVector(sys["path"]), joinpath("$(prefix)", "lib", "python" * string(sys[:version_info][1]) * "." * string(sys[:version_info][2]), "site-packages"))
"""
)

@BinDeps.install Dict(:lcm => :liblcm)


text:
module CryptoUtils

import Primes: isprime, primes

export legendre, jacobi, sqrt_mod_prime, find_quadratic_non_residue, is_quadratic_residue,
       convergents, continued_fraction, surd

export n2b, b2n

export random_prime, safe_prime, tower_two_prime, get_first_primes, twin_primes

export factor_with_ed, wiener

################################################################
##                   Number theory utilities                  ##
################################################################

"""
    legendre(a::Integer, p::Integer)

Return the Legendre symbol of `(a, p)`.

`p` should be an odd prime number.
"""
function legendre(a::Integer, p::Integer)
    if p % 2 == 0 && !isprime(p)
        throw("Argument p should be an odd prime.")
    end

    b = div(p - 1, 2)
    return powermod(a, b, p)
end

"""
    jacobi(n::Integer, k::Integer)

Return the Jacobi symbol of `n, k`.

`k` should be an odd number.
"""
function jacobi(n::Integer, k::Integer)
    rem(k, 2) != zero(k) || throw("Argument k=$k should be odd.")

    n = mod(n, k)
    t = 1
    while n != zero(n)
        while rem(n, 2) == 0
            n = div(n, 2)
            r = k % 8
            if r == 3 || r == 5
                t = -t
            end
        end
        n, k = k, n
        if n % 4 == 3 && k % 4 == 3
            t = -t
        end
        n = mod(n, k)
    end
    if k == oneunit(k)
        return t
    else
        return 0
    end
end

"""
    find_quadratic_non_residue(p::Integer)

Return a random number `R` which has no square root mod `p`,
i.e., `x^2 == R mod p` has no solutions.
"""
function find_quadratic_non_residue(p::Integer)
    two = big"2"
    qnr = rand(two:p - 1)

    while jacobi(qnr, p) != -1
        qnr = rand(two:p)
    end

    return qnr
end

"""
    tonelli_shanks(a::Integer, p::Integer)

Implements the Tonelli Shanks algorithm for computing
square roots modulo a prime number.

It assumes such square roots exist.
"""
function tonelli_shanks(a::Integer, p::Integer)
    b = find_quadratic_non_residue(p)
    t = p - 1
    s = 0
    while t % 2 == 0
        t = div(t, 2)
        s += 1
    end

    M = s
    Q = t
    c = powermod(b, t, p)
    t = powermod(a, t, p)
    r = powermod(a, div(Q + 1, 2), p)

    while t != 1
        i = 1
        tpow = powermod(t, 2, p)
        while tpow != 1
            tpow = mod(tpow * tpow, p)
            i += 1
        end
        po = M - i - 1
        b = powermod(c, big"2"^po, p)
        M = i
        c = powermod(b, 2, p)
        t = mod(t * powermod(b, 2, p), p)
        r = mod(r * b, p)
    end
    return r
end

"""
    hoc_sqrt(a::Integer, p::Integer)

Algorithm from Handbook of cryptography, Koblitz pp 48-49.
Finds a solution to `x^2 == a mod p`.

It assumes such solution exists.

Running time highly depends on |alpha|, assuming
`p-1 = 2^alpha * s`, for an odd `s`.
"""
function hoc_sqrt(a::Integer, p::Integer)
    n = find_quadratic_non_residue(p)
    s = p - 1
    alpha = 0
    while s % 2 == 0
        s = div(s, 2)
        alpha += 1
    end

    inv_a = invmod(a, p)
    b = powermod(n, s, p)
    r = powermod(a, div(s + 1, 2), p)
    r_sqr = powermod(r, 2, p)
    two = BigInt(2)
    expon = two^(alpha - 1)
    for i in 1:alpha - 1
        expon = div(expon, two)
        d = powermod(inv_a * r_sqr, expon, p)
        if d == p - 1
            r = mod(r * b, p)
            r_sqr = powermod(r, two, p)
        end
        b = powermod(b, two, p)
    end
    return r
end

"""
    is_quadratic_residue(a::Integer, p::Integer) -> Bool

Return true or false depending on wheter `a` is a quadratic residue mod `p`.

That is, it checks if `x^2 == a mod p` has solutions.
"""
function is_quadratic_residue(a::Integer, p::Integer)
    if p == oftype(p, 2)
        return true
    end
    return jacobi(a, p) == 1
end

"""
    sqrt_mod_prime(a::Integer, p::Integer) -> Integer

Solves `x^2 == a mod p` and returns one of the square roots `r`.
The other root is `p - r`.
If there are no solutions, throws an exception.

```julia
julia> sqrt_mod_prime(33^2, 73)
33
```
"""
function sqrt_mod_prime(a::Integer, p::Integer)
    a = mod(a, p)
    is_quadratic_residue(a, p) || throw("$a is not a quadratic residue mod $p.")

    if p % 2 == 0
        return a

    elseif p % 4 == 3
        return powermod(a, div(p + 1, 4), p)

    elseif p % 8 == 5
        d = powermod(a, div(p - 1, 4), p)

        if d == 1
            r = powermod(a, div(p + 3, 8), p)
        elseif d == p - 1
            r = mod(2 * a * powermod(4 * a, div(p - 5, 8), p), p)
        end

        return r

    # If p-1 is of the form 2^k*s for large k, use tonelli-shanks.
    # Here k is large if k > 100
    elseif mod(p - 1, 1267650600228229401496703205376) == 0
        return tonelli_shanks(a, p)

    # depends on size of k
    else
        return hoc_sqrt(a, p)
    end
end

"""
    continued_fraction(a::T, b::T) where T <: Integer

Return the continued fraction of the rational `a/b`.

# Example
```julia
julia> continued_fraction(31, 73)
6-element Array{Int64,1}:
 0
 2
 2
 1
 4
 2
```
"""
function continued_fraction(a::T, b::T) where T <: Integer
    r = a // b
    p = Integer(floor(r))
    fraction = typeof(a)[]
    push!(fraction, p)
    f = r - p
    while f != 0
        r = 1 // f
        p = Integer(floor(r))
        push!(fraction, p)
        f = r - p
    end
    return fraction
end

"""
    convergents(a::T, b::T) where T <: Integer

Return the convergents of a rational `a/b`.

# Example
```julia
julia> convergents(31, 73)
6-element Array{Rational,1}:
  0//1
  1//2
  2//5
  3//7
 14//33
 31//73
```
"""
function convergents(a::T, b::T) where T <: Integer
    fraction = continued_fraction(a, b)

    return convergents(fraction)
end

"""
    convergents(cont_fraction::Array)

Return the convergents given the continued fraction of a rational.
"""
function convergents(cont_fraction::Array)

    nums = BigInt[0, 1]
    dens = BigInt[1, 0]

    converg = Rational[]
    for a in cont_fraction
        next_num = a * nums[end] + nums[end - 1]
        next_den = a * dens[end] + dens[end - 1]

        push!(converg, next_num // next_den)
        push!(nums, next_num)
        push!(dens, next_den)
    end
    return converg
end

"""
    solve_quadratic(a::Integer, b::Integer, c::Integer)

Return solutions to the equation `a*x^2 + b*x + c`.

Assumes the solutions are integer.
"""
function solve_quadratic(a::Integer, b::Integer, c::Integer)
    d = isqrt(b^2 - 4 * a * c)
    return div(-b + d, 2 * a), div(-b - d, 2 * a)
end

"""
    surd(n::BigInt, k::Int64)

Return largest integer smaller or equal than the `k`-th root of `n`.
"""
function surd(n::BigInt, k::Int64)
    if k == 1
        return n
    elseif k == 0
        return big"1"
    end

    low = 1
    high = n
    mid = div(low + high, 2)
    while abs(mid^k - n) >= 0.49
        if mid^k < n
            low = mid
        elseif mid^k > n
            high = mid
        end
        mid = div(low + high, 2)
        if high - low == 1
            return low
        end
    end
    return floor(mid)
end

################################################################
##                   Prime number utilities                   ##
################################################################

"""
    random_prime(bitsize::Integer) -> BigInt

Return a random prime with `bitsize` bits.

```
julia> random_prime(42)
2458636110727
```
"""
function random_prime(bitsize::Integer)::BigInt
    n = oneunit(BigInt)
    two = BigInt(2)
    lo = two^(bitsize - 1)
    hi = two^bitsize
    while !isprime(n)
        n = rand(lo:hi)
    end
    return n
end

"""
    tower_two_prime(bitsize::Integer, tower_len::Integer) -> BigInt

Return a random prime of the form `2^towerlen * q + 1`
with `bitsize` bits and where `q` is also a prime.

```
julia> tower_two_prime(22, 6)
2362433
```
"""
function tower_two_prime(bitsize::Integer, tower_len::Integer)::BigInt
    n = oneunit(BigInt)
    tower = big"2"^tower_len
    lo = big"2"^(bitsize - tower_len - 1)
    hi = big"2"^(bitsize - tower_len)
    while !isprime(n)
        n = tower * random_prime(bitsize - tower_len) + 1
    end
    return n
end

"""
    safe_prime(bitsize::Integer) -> BigInt

Return a random safe-prime `q` of the form `q = 2 * p + 1`
where `p` is also a prime number.
The returning prime number has `bitsize` bits.

```
julia> safe_prime(10)
1439
```
"""
function safe_prime(bitsize::Integer)::BigInt
    n = oneunit(BigInt)
    two = BigInt(2)
    while !isprime(n)
        p = random_prime(bitsize - 1)
        n = two * p + 1
    end
    return n
end

"""
    get_first_primes(k::Integer) -> Collection

Output the first `k` prime numbers.

```julia
julia> get_first_primes(10)
10-element Array{Int64,1}:
2
3
5
7
11
13
17
19
23
29
```
"""
function get_first_primes(k::Integer)
    if k <= zero(k)
        return []
    end
    l = log(k + 1)
    hi = max(2, Int(ceil((k + 1) * (l + log(l)))))
    res = primes(hi)
    res = res[1:min(length(res), k)]
    return res
end

"""
    twin_primes(bitsize::Integer)

Return a pair of prime numbers `p, p + 2`
with `bitsize` bits.

This might take a while to run.
"""
function twin_primes(bitsize::Integer)
    q = 1
    p = 1
    while !isprime(q)
        p = random_prime(bitsize)
        q = p + 2
    end
    return p, q
end

################################################################
##                   Conversion utilities                     ##
################################################################

"""
    n2b(n::Integer) -> String

Converts a number to its bytes representation,
effectively writing the number in base 256, and returning
the corresponding bytes.

```julia
julia> n2b(22405534230753963835153736737)
"Hello world!"
```
"""
function n2b(n::Integer)::String
    res::String = ""
    while n != zero(n)
        n, r = divrem(n, 256)
        res = Char(r) * res
    end
    return res
end

"""
    b2n(str::String) -> BigInt

Converts a byte-string to a number,
converting the string from base 256 to base 10.

```julia
julia> b2n("Hello world!")
22405534230753963835153736737
```
"""
function b2n(str::String)::BigInt
    res::BigInt = 0
    for (i, char) in enumerate(reverse(str))
        res += big"256"^(i - 1) * Int(char)
    end
    return res
end

################################################################
##                   Cryptography utilities                   ##
################################################################

"""
    factor_with_ed(n::Integer, e::Integer, d::Integer) -> (Integer, Integer)

Factors `n = p*q` given `(e, d)` such that `e*d = 1 mod phi(n)`
Stinson page 204 - algorithm 5.10
"""
function factor_with_ed(n::Integer, e::Integer, d::Integer)
    k = e * d - 1
    s = 0
    r = k
    while r % 2 == 0
        s += 1
        r = div(r, 2)
    end

    while true
        w = rand(2:n - 1)
        x = gcd(w, n)
        if 1 < x < n
            return x, div(n, x)
        end

        v = powermod(w, r, n)
        v0 = v

        while v != 1
            v0 = v
            v = powermod(v, 2, n)
        end

        if v0 == n - 1
            continue
        end

        x = gcd(v0 + 1, n)
        if 1 < x < n
            return x, div(n, x)
        end
    end
end

"""
    wiener(n::Integer, e::Integer, dujella_bound=20)

Factors the semiprime `n`, assuming Wiener's attack holds:
`d < n^(1/4)`, where `d*e = 1 mod phi(n)`.

Uses Dujella extension attack. Increasing the `dujella_bound` argument
slows the running time but increases chances of finding the correct `d`
in case `d ~ n^(1/4)`.
"""
function wiener(n::Integer, e::Integer, dujella_bound = 20)
    # usual aproximation of k/d
    # convergs = convergents(e, n)

    # better aproximation of k/d
    convergs = convergents(e, n + 1 - 2 * isqrt(n))

    # ignore first convergent 0/1
    deleteat!(convergs, 1)

    old_d = 1

    # ciphertext to test decryption exponent
    test_cipher = powermod(2, e, n)

    # build (r, s) pairs
    # with 1 < s < r < dujella_bound
    # and gcd(r, s) = 1
    # https://bib.irb.hr/datoteka/383127.dujececc.pdf
    pairs_rs = Pair[]
    for r in 1:dujella_bound
        for s in 1:r
            if gcd(r, s) == 1
                push!(pairs_rs, r => s)
            end
        end
    end

    for fraction in convergs
        k = numerator(fraction)
        d = denominator(fraction)

        # regular Wiener attack
        if powermod(test_cipher, d, n) == 2
            phi = div(e * d - 1, k)
            p, q = solve_quadratic(1, phi - n - 1, n)
            @assert p * q == n
            return p, q
        end

        # Dujella extension
        for p in pairs_rs
            dujella = p.first * d + p.second * old_d
            if powermod(test_cipher, dujella, n) == 2
                return factor_with_ed(n, e, dujella)
            end
        end

        old_d = d
    end
end

end # module