text:
"""
    lectura_uw(red)
Creates a LightGraph's SimpleGraph object from an adjacency list array. Returns the SimpleGraph and the ordered node sequence,
this way node names can be strings.
...
# Arguments
- `red::Array`: Adjacency list.
...
"""
function lectura_uw(red)
    Nodes = union(unique(red[:,1]),unique(red[:,2]))
    g = SimpleGraph()
    last_node = Int64(length(Nodes))
    add_vertices!(g,last_node)
    for n in 1:size(red)[1]
        add_edge!(g,red[n,1],red[n,2])
        add_edge!(g,red[n,2],red[n,1])
    end
    return g, Nodes
end

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

"""
Generates a LFR network
...
# Arguments
- `N::Int`: Number of nodes.
- `k::Float`: Average degree.
- `maxk::Float`: Maximum degree.
- `mu::Float`: Mixing parameter.
- `t1::Float`:: Minus exponent for the degree sequence.
- `t2::Float`:: Minus exponent for the community size distribution.
- `minc::Int`: Minimum for the community sizes.
- `maxc::Int`: Maximum for the community sizes.
...
"""
function lfr_network(the_root, work_path;N = 100, k = 20, maxk = 40, mu = 0.1, t1 = 2, t2 = 1, minc = 10, maxc = 40)
    # the_root = pwd()
    run(`$the_root/benchmark -N $N -k $k -maxk $maxk -mu $mu -t1 $t1 -t2 $t2 -minc $minc -maxc $maxc`)

    out_files_path = joinpath(the_root, "..", "..", work_path)

    # the_net = readdlm(the_root*"/network.dat")
    # true_com = readdlm(the_root*"/community.dat")

    the_net = readdlm(joinpath(out_files_path, "network.dat"))
    true_com = readdlm(joinpath(out_files_path, "community.dat"))

    true_com = Int64.(true_com)
    return the_net, true_com
end

##=================####==============##
#
"""
    init_demographics!(agents; kwargs...)
Initialize agents' demographic attributes and
fraction of infected agents at t = 0
`coop_dist` -> Cooperation Distribution
`meets_dist` -> Distribution of number of meetings per time step
`recovt_dist` -> Agent's recovery time distribution
"""
function init_demographics!(agents; states::Array{String} = ["S","I"], initial::Array{Float64} = [0.5,0.5], coop_dist=[0], meets_dist=[0], recovt_dist=[0])
    for ag in agents
        ag.state = sample(states, Weights(initial))
        push!(ag.previous, ag.state)
        ag.p_cop      = rand(coop_dist)
        ag.num_meets  = 1 + rand(meets_dist)
        ag.recovery_t = ceil(rand(recovt_dist))
        if ag.state == "I"
            ag.counter = ag.counter+1
        end
    end
end

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

"""
    set_fixed_coop_agents!(agents, params)
Set cooperating agents with probability `p_coop_agents`
"""
function set_fixed_coop_agents!(agents; p_coop_agents=0.0)
    for ag in agents
        if rand() <= p_coop_agents
            ag.at_home = true
        end
    end
end

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

"""
    set_coop_agents!(agents, params)
Agent stays at home with probability `Agent.p_cop` at each time step
"""
function set_coop_agents!(agents; p_cop = 0.5)
    Threads.@threads for ag in agents
        if rand() <= p_cop
            ag.at_home = true
            ag.attitude = "ra"
            ag.coop_effect = 1.0
        else
            ag.at_home = false
        end
    end
end

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

"""
    set_adapt_agents!(agents, params)
Agent stays at home with probability `Agent.p_cop` at each time step
"""
function set_adapt_agents!(agents; p_cop = 0.5)
    Threads.@threads for ag in agents
        if rand() <= p_cop
            ag.adapter = true
        else
            ag.adapter = false
        end
    end
end
##=================####==============##

"""
    assign_contacts!(agent, all_agents, adj_mat, row)
Assign contacts to `agent` from adjacency matrix
"""
function assign_contacts!(g, agent)
    agent.contacts_t = neighbors(g, agent.id)
    agent.degree_t = degree(g, agent.id)
end

function get_coop!(agents)
    Threads.@threads for agent in agents
        agent.coopf = findall(x-> agents[x].at_home == true, agent.contacts_t)
        agent.non_coopf = findall(x-> agents[x].at_home == false, agent.contacts_t)
    end
end

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

"""
This function looks at the infected neighbors of each agent to compute wether or not
she will get infected. If both agents adopt a cooperative behavior then the probability
of infection gets reduced by both terms. As the SIS and SIR models are basically the same
excep for reinfections then only a flag for SIR is needed.
"""
function SI_coop!(agent, agents;inf_prob=0.1, rec_prob=0.3, coop_red=0.7, R=false)
    infected_coop = [ag.state for ag in agents[agent.coopf] if ag.state == "I"] |> length
    infected_noncoop = [ag.state for ag in agents[agent.non_coopf] if ag.state == "I"] |> length
    if agent.state == "S"
        if agent.at_home
            inf_prob_co = inf_prob * (1-coop_red)^2
            inf_prob_no = inf_prob * (1-coop_red)
        else
            inf_prob_co = inf_prob * (1-coop_red)
            inf_prob_no = inf_prob
        end
        oddsc = [sample([true,false],Weights([inf_prob_co,1-inf_prob_co])) for i in 1:infected_coop] |> sum
        oddsn = [sample([true,false],Weights([inf_prob_no,1-inf_prob_no])) for i in 1:infected_noncoop] |> sum
        odds = oddsc + oddsn
        if odds > 0
            agent.new_state = "I"
        else
            agent.new_state = "S"
        end
    elseif agent.state == "I"
        if !R
            if rand() <= rec_prob
                agent.new_state = "S"
            else
                agent.new_state = "I"
            end
        else
            if rand() <= rec_prob
                agent.new_state = "R"
            else
                agent.new_state = "I"
            end
        end
    end
end

function SI_attitude!(agent, agents;inf_prob=0.1, rec_prob=0.3, R=false)
    #  sus = findall(x -> x.state == "S",agents)
    if agent.state == "S"
        infec = [ag.coop_effect for ag in agents[agent.contacts_t] if ag.state == "I"]
        the_odds = @. inf_prob * (1-agent.coop_effect) * (1-infec)
        odds = the_odds |> f -> map(x -> sample([true,false], Weights([x,1-x])),f) |> sum
        if odds > 0
            agent.new_state = "I"
        else
            agent.new_state = "S"
        end
    elseif agent.state == "I"
        if !R
            if rand() <= rec_prob
                agent.new_state = "S"
            else
                agent.new_state = "I"
            end
        else
            if rand() <= rec_prob
                agent.new_state = "R"
            else
                agent.new_state = "I"
            end
        end
    end
end

function next_state!(agents;fun=SI_coop!,kwargs...)
    for agent in agents
        fun(agent,agents;kwargs...)
    end
end

function SI_next!(agents; inf_prob=0.1, rec_prob=0.3)
    for agent in agents
        if agent.state == "S"
            infected = [ag.state for ag in agents[agent.contacts_t] if ag.state == "I"] |> length
            odds = [sample([true,false],Weights([inf_prob,1-inf_prob])) for i in 1:infected] |> sum
            if odds > 0
                agent.new_state = "I"
            else
                agent.new_state = "S"
            end
        elseif agent.state == "I"
            if rand() <= rec_prob
                agent.new_state = "S"
            else
                agent.new_state = "I"
            end
        end
    end
end

function SIR_next!(agents; inf_prob=0.1, rec_prob=0.3)
    for agent in agents
        if agent.state == "S"
            infected = [ag.state for ag in agents[agent.contacts_t] if ag.state == "I"] |> length
            odds = [sample([true,false],Weights([inf_prob,1-inf_prob])) for i in 1:infected] |> sum
            if odds > 0
                agent.new_state = "I"
            else
                agent.new_state = "S"
            end
        elseif agent.state == "I"
            if rand() <= rec_prob
                agent.new_state = "R"
            else
                agent.new_state = "I"
            end
        end
    end
end

"""
    update_state!(agents)
Finds agent's next state and updates it, it just packages
`get_next_state!` and `update_state!`
"""
function update_state!(agents)
    Threads.@threads for ag in agents
        push!(ag.previous,ag.state)
        if (ag.state == "S" && ag.new_state == "I")
            ag.counter = ag.counter+1
        end
        ag.state = ag.new_state
    end
end

function update_coop!(agents,threshold;lrt=false)
    Threads.@threads for ag in agents
        if ag.adapter
            if length(ag.coopf)/ag.degree_t > threshold
                ag.at_home = true
            else
                if lrt
                    ag.at_home = false
                end
            end
        end
    end
end

function update_coop_distance!(agents,g,d,threshold;lrt=false)
    Threads.@threads for ag in agents
        if ag.adapter
            status = agents[neighborhood(g,ag.id,d)|> unique] |> f -> map(x -> x.at_home,f)
            thecoop = sum(status)
            tot = length(status)
            if thecoop/tot >= threshold
                ag.at_home = true
            else
                if lrt
                    ag.at_home = false
                end
            end
        end
    end
end

function prev_infected(agent,threshold)
    historic = agent.previous[end-threshold:end]
    Infected = filter(x->x=="I",historic)
    Susceptible = filter(x->x=="S",historic)
    Recovered = filter(x->x=="R",historic)
    Infected, Susceptible, Recovered
end

function update_coop_infections!(agents,threshold;lrt=false)
    Threads.@threads for ag in agents
        if ag.adapter
            Infected, Susceptible, Recovered = prev_infected(ag, threshold)
            if length(Infected) >= 1
                ag.at_home = true
                #  ag.at_home = !(ag.at_home)
            else
                if lrt
                    ag.at_home = false
                    #  ag.at_home = !(ag.at_home)
                end
            end
        end
    end
end

function update_coop_distance_inf!(agents,g,d,threshold;lrt=false)
    Threads.@threads for ag in agents
        if ag.adapter
            status = agents[neighborhood(g,ag.id,d)|> unique] |> f -> map(x -> prev_infected(x,threshold)[1],f)
            theinfected = sum(map(x -> "I" in x, status))
            tot = length(status)
            if theinfected >= 1
                ag.at_home = true
            else
                if lrt
                    ag.at_home = false
                end
            end
        end
    end
end

##### Check infections and change behavior depending on distance

function update_single_given_distance!(agents,g,v,d,threshold,probs)
    neigh = neighborhood_dists(g,v,d)
    nodes = first.(neigh)
    distances = last.(neigh)
    changed = false
    for dist in 0:d
        current = findall(x->x==dist,distances)
        infected = findall(x->x.counter >= threshold, agents[current])
        for time in 1:length(infected)
            if rand() < probs[dist+1]
                agents[v].at_home = true
                changed = true
                break
            end
        end
        if changed
            break
        end
    end
end

function update_coop_given_distance!(agents,g,d,threshold,probs;lrt=false)
    Threads.@threads for ag in agents
        if ag.adapter
            update_single_given_distance!(agents,g,ag.id,d,threshold,probs)
        end
    end
end

function update_single_effect_distance!(agents,g,d,threshold,step;v)
    neigh = neighborhood_dists(g,v,d)
    nodes = first.(neigh)
    distances = last.(neigh)
    for dist in unique(distances)
        current = findall(x->x==dist,distances)
        current_nodes = nodes[current]
        infected = findall(x->x.counter >= threshold, agents[current_nodes])
        for time in 1:length(infected)
            if agents[v].attitude == "ra"
                agents[v].coop_effect = max(agents[v].coop_effect-step[dist+1],0)
            elseif agents[v].attitude == "rt"
                agents[v].coop_effect = min(agents[v].coop_effect+step[dist+1],1)
            end
        end
    end
end

function update_effect_given_distance!(agents,g,d,threshold,step)
    Threads.@threads for ag in agents
        if ag.adapter
            update_single_effect_distance!(agents,g,d,threshold,step;v=ag.id)
        end
    end
end

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

function update_single_effect_distance_coop!(agents,g,d,threshold,step;v)
    neigh = neighborhood_dists(g,v,d)
    nodes = first.(neigh)[2:end]
    distances = last.(neigh)[2:end]
    for dist in unique(distances)
        current = findall(x->x==dist,distances)
        current_nodes = nodes[current]
        mean_effect = mean([x.coop_effect for x in agents[current_nodes]])
        if agents[v].attitude == "ra"
            agents[v].coop_effect = max(agents[v].coop_effect-mean_effect*step[dist+1],0)
        elseif agents[v].attitude == "rt"
            agents[v].coop_effect = min(agents[v].coop_effect+mean_effect*step[dist+1],1)
        end
    end
end

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

function update_effect_given_distance_coop!(agents,g,d,threshold,step)
    Threads.@threads for ag in agents
        if ag.adapter
            update_single_effect_distance_coop!(agents,g,d,threshold,step;v=ag.id)
        end
    end
end

################################################################################################################################################
################################################################################################################################################
################################################################################################################################################
################################################################################################################################################
################################################################################################################################################
################################################################################################################################################
################################################################################################################################################
################################################################################################################################################
################################################################################################################################################
################################################################################################################################################
################################################################################################################################################
################################################################################################################################################
################################################################################################################################################
################################################################################################################################################
################################################################################################################################################
################################################################################################################################################
################################################################################################################################################
################################################################################################################################################
################################################################################################################################################

"""
    get_next_state!(agent, params)
Finds `agent.new_state` according to its contacts and meetings
"""
function get_next_state!(agent;now_t)

    # IF AGENT IS SUCEPTIBLE AND IS *NOT* AT HOME...
    if agent.state == "S" && agent.at_home == false

        if agent.num_meets <= agent.degree_t
            meets = sample(agent.contacts_t, agent.num_meets, replace=false)
        else
            meets = agent.contacts_t
        end

        # println(agent.id, "|meets:", [x.id for x in meets])

        for c_agent in meets # LOOP THROUGH AGENTS MET THAT TIMESTEP
            coll_state = agent.state * c_agent.state
            # SUCEPTIBLE <- INFECTED
            if coll_state == "SI" && c_agent.at_home == false
                # println(coll_state)
                if rand() < params.attack_rate
                    agent.new_state = "I"
                    agent.infection_t = now_t
                    break
                end
            end
        end
    elseif agent.state == "I" # INFECTED AGENT TO RECOVER

        time_from_infection = now_t - agent.infection_t
        if time_from_infection == agent.recovery_t
            agent.new_state = "R"
            # println(agent.id, "|", time_from_infection, "|RECOVERED!")
        end

    end
end

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

#  """
#      update_state!(agent)
#  Updates `agent.state` with the one computed in `get_next_state!`
#  """
#  function update_state!(agent)
#      agent.state = agent.new_state
#  end

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

"""
    update_all_agents!(agents, params)
Finds agent's next state and updates it, it just packages
`get_next_state!` and `update_state!`
"""
function update_all_agents!(agents, params)
    # FIND AGENTS' NEXT STATE FROM INTERACTIONS
    Threads.@threads for ag in agents
        get_next_state!(ag, params)
    end

    # UPDATE AGENTS' STATE
    Threads.@threads for ag in agents
        update_state!(ag)
    end
end

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

"""
    get_populations(agents, params)
Get number of agents in each state:
    - S: Suceptible
    - I: Infected
    - R: Recovered
    - Q: Quarantained (not implemented)
    - D: Decesead (not implemented)
"""
function get_populations(agents, params)
    # species = num_S, num_I, num_R, num_Q
    species = zeros(4)

    for ag in agents
        if ag.state == "S"
            species[1] += 1.0
        elseif ag.state == "I"
            species[2] += 1.0
        elseif ag.state == "R"
            species[3] += 1.0
        # else
        #     species[4] += 1.0
        end
    end
    return species ./ params.num_agents
end
##=================####==============##

"""
    export_parameters(params, out_path)
Export parameters to file
"""
function export_parameters(params, out_path)

    filename = "parameters_N_$(params.num_agents)_rep_$(params.repetition).txt"

    open(joinpath(out_path, filename), "w") do io
    println(io, "num_agents   :\t", params.num_agents   )
    println(io, "p_link       :\t", params.p_link       )
    println(io, "p_infected_t0:\t", params.p_infected_t0)
    println(io, "attack_rate  :\t", params.attack_rate  )
    println(io, "repetition   :\t", params.repetition   )
    end

end
##=================####==============##


text:
struct ConfusionMatrix
    matrix
    classlabels

    function ConfusionMatrix(tn, fp, fn, tp; classlabels=nothing)
        ConfusionMatrix([tn fp ;fn tp],classlabels=classlabels)
    end

    function ConfusionMatrix(matrix::Matrix; classlabels=nothing)
        if classlabels == nothing
            _, n = size(matrix)
            classlabels = ["class $x" for x in 1:n]
        end
        new(matrix, classlabels)
    end

end

function tohtml(confusionmatrix::ConfusionMatrix)
    cm = confusionmatrix.matrix
    labels = confusionmatrix.classlabels

    m,n = size(cm)
    stream = IOBuffer()
    write(stream,"""
""")
    write(stream,"
")
    write(stream,"""	""")
    write(stream,"""Predicted class x&#770	""")
    write(stream,"
")
    write(stream,"
")
    for label in labels
        write(stream,"$(label)	")
    end
    write(stream,"Row total	")
    write(stream,"Accuracy  P(x&#770|x)	")
    write(stream,"
")
    write(stream,"
")
    write(stream,"""""")
    write(stream,"""
Actual

class x
""")
    write(stream,"	")
    write(stream,"$(labels[1])	")
    for col in cm[1,:]
        write(stream,"$col	")
    end
    write(stream,"$(sum(cm[1,:]))	") # Row total
    write(stream,"""$(@sprintf("%.2f",(cm[1,1] / (sum(cm[1,:])))))	""") # Accuracy
    write(stream,"
")
    for i in 2:length(labels)
        write(stream,"
")
        write(stream,"$(labels[i])	")
        for col in cm[i,:]
            write(stream,"$col	")
        end
        write(stream,"$(sum(cm[i,:]))	") # Row total
        write(stream,"""$(@sprintf("%.2f",(cm[i,i] / (sum(cm[i,:])))))	""") # Accuracy
        write(stream,"
")
    end
    write(stream,"
")
    write(stream,"""	""")
    write(stream,"Column total	")
    for j in 1:length(labels)
        write(stream,"$(sum(cm[:,j]))	")
    end
    write(stream,"$(sum(cm))	")
    write(stream,"
")
    write(stream,"
")
    write(stream,"""	""")
    write(stream,"Reliability P(x|x&#770)	")
    for i in 1:length(labels)
        write(stream,"""$(@sprintf("%.2f",cm[i,i] / (sum(cm[:,i]))))	""") # Reliability
    end
    write(stream,"""	""")
    write(stream,"""$(@sprintf("%.2f",(sum([cm[i, i] for i = 1:m]) / sum(cm))))	""")
    write(stream,"
")
    write(stream,"
")
    String(stream)
end

Base.show(io::IO, ::MIME"text/html", c::ConfusionMatrix) =
    print(io, tohtml(c))

ext: jl


text:
import Base: length, push!, append!, isempty, collect, first, last, show
using DataStructures
using LightGraphs

mutable struct DoublyLinkedNode{T}
    data::Union{Nothing, T}
    next::DoublyLinkedNode{T}
    prev::DoublyLinkedNode{T}
    function DoublyLinkedNode{T}() where T
        x = new()
        x.data = nothing
        x.next = x
        x.prev = x
        x
    end
    function DoublyLinkedNode{T}(
            data::T,
            next::DoublyLinkedNode{T},
            prev::DoublyLinkedNode{T}) where T
        new(data, next, prev)
    end
end

function show(io::IO, x::DoublyLinkedNode)
    print(io, "LinkedNode($(x.data))")
end

mutable struct DoublyLinkedList{T}
    head::DoublyLinkedNode{T}
    size::Int
    function DoublyLinkedList{T}() where T
        new(DoublyLinkedNode{T}(), 0)
    end
end

function print_next(x::DoublyLinkedList, n::Int)
    v = x.head
    for i in 1:n
        print(v.data)
        v = v.next
        print("-")
    end
    print("\n")
end

function print_prev(x::DoublyLinkedList, n::Int)
    v = x.head
    for i in 1:n
        print(v.data)
        v = v.prev
        print("-")
    end
    print("\n")
end

function show(io::IO, x::DoublyLinkedList)
    print(io, "LinkedList{$(typeof(x.head))} of size $(x.size)")
end

function push!(l::DoublyLinkedList{T}, x::T) where T
    curr_last = l.head.prev
    new_node = DoublyLinkedNode{T}()
    new_node.data = x
    new_node.next = l.head
    new_node.prev = curr_last
    l.head.prev = new_node
    curr_last.next = new_node
    l.size += 1
    return
end

function shift!(l::DoublyLinkedList{T}, x::T) where T
    curr_first = l.head.next
    new_node = DoublyLinkedNode{T}()
    new_node.data = x
    new_node.next = curr_first
    new_node.prev = l.head
    curr_first.prev = new_node
    l.head.next = new_node
    l.size += 1
    return
end

function delete_last!(l::DoublyLinkedList{T}) where T
    l.head.prev.prev.next = l.head
    l.head.prev = l.head.prev.prev
    l.size -= 1
    return
end

function collect(l::DoublyLinkedList{T}) where T
    out = T[]
    v = l.head.next
    for i in 1:l.size
        push!(out, v.data)
        v = v.next
    end
    out
end

function reverse!(l::DoublyLinkedList{T}) where T
    v = l.head
    for i in 1:l.size + 1
        w = v.next
        v.next, v.prev = v.prev, v.next
        v = w
    end
end

function reversed(l::DoublyLinkedList{T}) where T
    rev = DoublyLinkedList{T}()
    rev.size = l.size
    v = rev.head
    w = l.head.prev
    for i in 1:l.size
        v.next = DoublyLinkedNode{T}(w.data, w.prev, w.next)
        v.next.prev = v
        v = v.next
        w = w.prev
    end
    v.next = rev.head
    rev.head.prev = v
    rev
end

function DoublyLinkedList(x::AbstractVector{T}) where T
    l = DoublyLinkedList{T}()
    for xi in x
        push!(l, xi)
    end
    l
end

function append!(
        l::DoublyLinkedList{T},
        l2::DoublyLinkedList{T}) where T
    # let's deal with trivial cases
    if l2.size == 0
        return
    elseif l.size == 0
        l.head = l2.head
        l.size = l2.size
        return
    else
        l.head.prev.next = l2.head.next
        l2.head.next.prev = l.head.prev
        l.head.prev = l2.head.prev
        l2.head.prev.next = l.head
        l.size += l2.size
        return
    end
end

function cat(
        l::DoublyLinkedList{T},
        l2::DoublyLinkedList{T}) where T
    lnew = DoublyLinkedList{T}()
    lnew = deepcopy(l)
    append!(lnew, l2)
end

function first(l::DoublyLinkedList{T}) where T
    l.head.next.data
end

function last(l::DoublyLinkedList{T}) where T
    return l.head.prev.data
end

struct LinkedGraphNode
    id::Int
    nbrs::Vector{LinkedGraphNode}
    inv_nbrs::Vector{Int} # stores the position in the neighbor of the edge to quickly mark as visited
    visited::BitVector
    function LinkedGraphNode(id::Int)
        new(id, LinkedGraphNode[], Int[], falses(0))
    end
end

function show(io::IO, v::LinkedGraphNode)
    msg = "Node($(v.id), $(sum(v.visited))/$(length(v.nbrs)))"
    print(io, msg)
end

struct LinkedGraph
    root::LinkedGraphNode
    num_nodes::Int
    num_edges::Int
end

function show(io::IO, g::LinkedGraph)
    nv = g.num_nodes
    ne = g.num_edges
    msg = "LinkedGraph with $(nv) nodes and $(ne) edges"
    print(io, msg)
end

function find_unvisited(v::LinkedGraphNode)
    n = length(v.nbrs)
    i = 1
    pos_free = n + 1
    num_free = 0
    while i <= n
        if !v.visited[i]
            pos_free = min(pos_free, i)
            num_free += 1
            num_free >= 2 && break
        end
        i += 1
    end
    pos_free, num_free
end

function mark_as_visited!(v::LinkedGraphNode, i::Int)
    v.visited[i] = true
    v.nbrs[i].visited[v.inv_nbrs[i]] = true
end

function LinkedGraph(g::SimpleGraph, root::Int = 1)
    @assert length(connected_components(g)) == 1
    nodes = [LinkedGraphNode(i) for i in 1:nv(g)]
    for e in edges(g)
        i, j = src(e), dst(e)
        κ, η = length(nodes[i].nbrs), length(nodes[j].nbrs)
        push!(nodes[i].nbrs, nodes[j])
        push!(nodes[i].inv_nbrs, η + 1)
        push!(nodes[i].visited, false)
        push!(nodes[j].nbrs, nodes[i])
        push!(nodes[j].inv_nbrs, κ + 1)
        push!(nodes[j].visited, false)
    end
    LinkedGraph(nodes[root], nv(g), ne(g))
end

mutable struct Branch
    path::DoublyLinkedList{LinkedGraphNode}
    deletable::Int
end

function show(io::IO, b::Branch)
    msg = "Branch(size=$(b.path.size), del=$(b.deletable))"
    print(io, msg)
end

function is_cycle(
        v::LinkedGraphNode,
        τ::DoublyLinkedList{LinkedGraphNode})
    v == first(τ)
end

function merge_branches_at_fork!(
    v::LinkedGraphNode,
    branches::Stack{Branch},
    forks::Stack{LinkedGraphNode},
)
    subpaths = DoublyLinkedList{LinkedGraphNode}[]
    deletables = Int[]
    while true
        branch = pop!(branches)
        pop!(forks)
        push!(subpaths, branch.path)
        push!(deletables, branch.deletable)
        # if no more branches or next branch belongs to a different fork
        (isempty(branches) || top(forks) != v) && break
    end

    merged = DoublyLinkedList{LinkedGraphNode}()

    largest = 0
    maxlen = 0
    for (i, τ) in enumerate(subpaths)
        len = τ.size - deletables[i]
        if len > maxlen
            largest = i
            maxlen = len
        end
    end

    τ = subpaths[largest]
    deletable = deletables[largest] + (is_cycle(v, τ) ?  0 : maxlen)
    if !is_cycle(v, τ)
        push!(merged, v)
        append!(merged, reversed(τ))
        delete_last!(merged)
    end
    append!(merged, τ)

    for (i, τ) in enumerate(subpaths)
        if i != largest
            if !is_cycle(v, τ)
                push!(merged, v)
                append!(merged, reversed(τ))
                delete_last!(merged)
            end
            append!(merged, τ)
        end
    end

    branch = Branch(merged, deletable)

    return branch
end

function quicktour(g::LinkedGraph)
    seglens = Stack{Int}()
    branches = Stack{Branch}()
    forks = Stack{LinkedGraphNode}()
    path = Stack{LinkedGraphNode}()

    curr_len = 1
    v = g.root
    push!(path, v)
    push!(seglens, 1)

    while !isempty(seglens)
        push!(path, v)
        i_free, num_free = find_unvisited(v)
        if num_free >= 2
            # case a: commit as fork an advance
            push!(seglens, curr_len)
            curr_len = 1
        elseif num_free == 1
            # case b: proceed to nbr
            curr_len += 1
        else
            # case 3: dead end found, add to branches and check if done with fork
            # flush to branch
            branch_path = DoublyLinkedList{LinkedGraphNode}()
            for i in 1:curr_len
                v = pop!(path)
                push!(branch_path, v)
            end
            push!(branches, Branch(branch_path, 0))
            v = top(path) # now v == curr_fork
            push!(forks, v)

            # if done with fork merge segments, otherwise continue on available
            i_free, num_free = find_unvisited(v)
            if num_free == 0
                # branch merging until finding available path to move
                pop!(path)
                while true
                    #! TODO! DELETE ON THE LEFT AS WELL
                    branch = merge_branches_at_fork!(v, branches, forks)

                    # move to next fork
                    curr_len = pop!(seglens)
                    if isempty(seglens)
                        τ = collect(branch.path)
                        sol = [x.id for (i, x) in enumerate(τ) if i > branch.deletable]
                        return sol
                    end
                    for i = 1:curr_len
                        push!(branch.path, v)
                        v = pop!(path)
                        shift!(branch.path, v)
                        branch.deletable += 1
                    end
                    push!(forks, v)
                    push!(branches, branch)

                    i_free, num_free = find_unvisited(v)
                    num_free > 0 && break
                end
                push!(path, v)
            end

            curr_len = 1
        end

        # advance
        mark_as_visited!(v, i_free)
        v = v.nbrs[i_free]
    end
end

function quicktour(g0::SimpleGraph, k::Int)
    g = LinkedGraph(g0, k)
    quicktour(g)
end

function count_unique_edges(tour::Vector{Int})
    return length(unique([sorted(tour[i], tour[i + 1]) for i in 1:length(tour) - 1]))
end

function quicktour(g0::SimpleGraph)
    v = findfirst(k -> isodd(degree(g0, k)), vertices(g0))
    (v == 0) && (v = 1)
    g = LinkedGraph(g0, v)
    quicktour(g)
end

# """
# Implements pseudotour trails decomposition algorithm from Tansey
# """
function pseudotour_trails(g::AbstractGraph; as_trails::Bool=false)
    odds = [i for i in vertices(g) if isodd(degree(g, i))]
    n_odds = length(odds)
    @assert iseven(n_odds) "pseudotour not possible: number of odd vertices must be even"

    trails = Vector{Vector{Int}}(0)
    if n_odds == 0
        # single tour is possible
        append!(trails, [euler(g, 1)])
    else
        # must add a new vertex to all odd vertices
        n = nv(g)
        g_ = copy(g)
        add_vertex!(g_)
        for v in odds
            add_edge!(g_, n + 1, v)
        end

        # now find a large tour
        tour = euler(g_, 1)
        breaks = findall(tour .== n + 1)

        # split the tour based on the artifial vertex
        if breaks[1] != n + 1
            prepend!(tour, n + 1)
            breaks += 1
            prepend!(breaks, 1)
        end
        if breaks[end] != n + 1
            append!(tour, n + 1)
            append!(breaks, length(tour))
        end
        for i in 1:length(breaks) - 1
            append!(trails, [tour[breaks[i] + 1:breaks[i + 1] - 1]])
        end
    end
    trails = [t for t in trails if length(t) > 1]

    return as_trails ? trailvec_to_trails(trails, nv(g)) : trails
end

function trailvec_to_trails(intvec::AbstractVector{Vector{Int}}, n::Int)
    ptr = vcat(intvec...)
    wts = vcat([[ones(length(τ) - 1); 0.] for τ in intvec]...)
    brks = [1; 1 + cumsum(length.(intvec))]
    return Trails(ptr, brks, wts, n)
end

function euler(g::AbstractGraph, to::Int)
    # check input ---
    @assert has_vertex(g, to) "start is not in g"

    # neighbour map
    g_ = copy(g)
    edge_count = degree(g_)

    # eulerian?
    @assert length(connected_components(g)) == 1 "the graph is not connected"
    n_odds = sum(isodd.(edge_count))
    @assert n_odds ∈ (0, 2) "number of odd vertices must be zero or two"
    @assert n_odds == 0 || isodd(edge_count[to]) "the graph has two odd vertices but start is not one of them"
    # checks ready ---

    # stack to keep vertices
    curr_path = Stack{Int}()

    # stores final circuit
    circuit = Int[]

    # start from start
    curr_v = to
    push!(curr_path, to)

    while !isempty(curr_path)
        # There's a remaining edge
        remaining = edge_count[curr_v]
        if remaining > 0
            # add vertex to current path
            push!(curr_path, curr_v)

            # find the next vertex
            next_v = neighbors(g_, curr_v)[1]

            # reduce number of available edges
            edge_count[curr_v] -= 1
            edge_count[next_v] -= 1
            rem_edge!(g_, curr_v, next_v)

            # update current vertex
            curr_v = next_v
        else
            # add vertex without edges into circuit
            push!(circuit, curr_v)

            # backtrack until a vertex has edges left
            curr_v = pop!(curr_path)
        end
    end

    return circuit
end

struct Trails
    ptr::Vector{Int}
    brks::Vector{Int}
    wts::Vector{Float64}
    num_nodes::Int

    # Validate input
    function Trails(
            ptr::Vector{Int},
            brks::Vector{Int},
            wts::Vector{Float64},
            num_nodes::Int)
        @assert issorted(brks)
        @assert brks[1] == 1
        @assert brks[end] == length(ptr) + 1
        @assert length(wts) == length(ptr)
        new(ptr, brks, wts, num_nodes)
    end

    # Obtain weights and num_nodes from ptr
    function Trails(
            ptr::Vector{Int},
            brks::Vector{Int})
        @assert issorted(brks)
        @assert brks[1] == 1
        @assert brks[end] == length(ptr) + 1
        num_nodes = length(unique(ptr))
        cnts = counter(Tuple{Int, Int})
        for (a, b) in zip(ptr[1:end-1], ptr[2:end])
            x, y = min(a, b), max(a, b)
            cnts[x, y] += 1
        end
        wts = zeros(length(ptr))
        for (i, (a, b)) in enumerate(zip(ptr[1:end-1], ptr[2:end]))
            if i + 1 ∉ brks
                x, y = min(a, b), max(a, b)
                wts[i] += 1.0 / cnts[x, y]
            end
        end
        new(ptr, brks, wts, num_nodes)
    end
end

function show(io::IO, trails::Trails)
    print("Trails object with ", num_trails(trails), " trail(s)")
end

function trail(trails::Trails, i::Int)
    trails.ptr[trails.brks[i]:(trails.brks[i + 1] - 1)]
end

length(trails::Trails) = length(trails.ptr)
num_trails(trails::Trails) = length(trails.brks) - 1
num_nodes(trails::Trails) = trails.num_nodes

sorted(e::Edge) = sorted(src(e), dst(e))

function graph_from_edgelist(
        edges::AbstractArray;
        from_zero::Bool=false)
    @assert size(edges, 2) == 2
    n = maximum(vec(edges))
    g = SimpleGraph(n + Int(from_zero))
    for (i, j) in zip(edges[:, 1], edges[:, 2])
        add_edge!(g, i + Int(from_zero), j + Int(from_zero))
    end
    return g
end

function find_trails(g::SimpleGraph; ntrails::Int = 1)
    @time ptr_chinese = quicktour(g)
    N = length(ptr_chinese) + ntrails - 1
    chunk_size = if N % ntrails == 0
        N / ntrails
    else
        (N + ntrails) ÷ (ntrails)
    end
    ptr = Int[]
    brks = [1]
    for i in 1:ntrails - 1
        chunk_start = (i - 1) * chunk_size + 1
        chunk_end = i * chunk_size
        append!(ptr, ptr_chinese[chunk_start:chunk_end])
        push!(brks, length(ptr) + 1)
        push!(ptr, ptr[end])
    end
    chunk_start = (ntrails - 1) * chunk_size + 1
    append!(ptr, ptr_chinese[chunk_start:end])
    push!(brks, length(ptr) + 1);
    Trails(ptr, brks)
end

function grid_trails(nr::Int, nc::Int)
    ptr = Int[]
    brks = Int[]
    idx = LinearIndices((nr, nc))
    for j=1:nc, i in 1:nr
        push!(ptr, idx[i, j])
        (i == 1) && push!(brks, idx[i, j])
    end
    for i in 1:nr, j=1:nc
        push!(ptr, idx[i, j])
        (i == 1) && push!(brks, nr * nc + idx[i, j])
    end
    push!(brks, 2 * nr * nc + 1)
    ptr, brks
end

# g0 = Graph(7)
# add_edge!(g0, 1, 2)
# add_edge!(g0, 1, 3)
# add_edge!(g0, 2, 4)
# add_edge!(g0, 2, 5)
# add_edge!(g0, 3, 6)
# add_edge!(g0, 3, 7)
# gplot(g0, nodelabel=1:7)

# g = Graph(18)
# add_edge!(g, 1, 2)
# add_edge!(g, 2, 3)
# add_edge!(g, 2, 4)
# add_edge!(g, 2, 6)
# add_edge!(g, 4, 5)
# add_edge!(g, 5, 7)
# add_edge!(g, 6, 7)
# add_edge!(g, 6, 8)
# add_edge!(g, 8, 9)
# add_edge!(g, 9, 10)
# add_edge!(g, 9, 11)
# add_edge!(g, 10, 13)
# add_edge!(g, 11, 12)
# add_edge!(g, 13, 14)
# add_edge!(g, 13, 15)
# add_edge!(g, 15, 16)
# add_edge!(g, 14, 10)
# add_edge!(g, 13, 17)
# add_edge!(g, 11, 18)
# gplot(g, nodelabel=1:nv(g))
# g2 = deepcopy(g)
# tour = quicktour(g)

# srand(321654)
# for i in 1:1000000
#     seed = rand(1:1000000)
#     srand(seed)
#     g = erdos_renyi(10, .25)
#     if length(connected_components(g)) == 1
#         tour = quicktour(g)
#         evisited = count_unique_edges(tour)
#         if ne(g) != evisited
#             println("seed = $seed")
#             for l in 1:length(tour) - 1
#                 v, w = tour[l], tour[l + 1]
#                 if !has_edge(g, v, w)
#                     println("$v-$w")
#                 end
#             end
#             gplot(g, nodelabel=vertices(g))
#             break
#         end
#     end
# end


text:
module INIBuilder

using REPL
using REPL.TerminalMenus
using ..Circuitscape

include("filepicker.jl")
include("run.jl")

export start

end


text:
"Random Variable: a function `Ω -> T` on a probability space"
abstract type RandVar <: NonDetVar end

MaybeRV{T} = Union{T, RandVar} where T

id(rv::RandVar) = rv.id

"Parameters of a random variable"
function params(rv::RandVar) end

"Name of a random variable"
name(rv) = rv

# 1.3 supports abstract fields for function call overload
if VERSION >= v"1.3"
  (rv::RandVar)(ω::Ω) = apl(rv, ω)
end

"""Is `x` a constant random variable, i.e. x(ω) = c ∀ ω ∈ Ω?

Determining constancy is intractable (and likely undecidable) in general.
if `isconstant(x)` is true then `x` is constant, but if `isconstant(x)` is false,
`x` may still be constant, but we have failed to determine it.

```jldoctest
x1 = constant(0.3)
isconstant(x1)
true

x2 = ciid(ω -> 0.3)
isconstant(x2)
true

x3 = ciid(ω -> rand(ω))
isconstant(x3)
false

# False Negative
x3 = ciid(ω -> rand(ω) > 0.5 ? 0.3 : 0.3)
isconstant(x3)
false
```
"""
function isconstant(x::RandVar, ΩT = defΩ(x))
  # This implementation assumes that ΩT is lazy, more general solution would wrap
  # ω of any type and intercept rand(ω)
  ω = ΩT()
  x(ω)
  isempty(ω)
end

isconstant(x) = true

# Printing #
Base.show(io::IO, rv::RandVar) =
  print(io, "$(id(rv)):$(name(rv))($(join(map(name, params(rv)), ", ")))::$(elemtype(rv))")


text:
module INP

export  InpStiffness

using ...TopOptProblems: Metadata, StiffnessTopOptProblem
using JuAFEM

include(joinpath("Parser", "Parser.jl"))
using .Parser

include("inpstiffness.jl")

end


text:
using .tensorboard: NodeDef, AttrValue_ListValue, AttrValue, NameAttrList

"""
	log_graph
"""
function log_graph(logger::TBLogger, g::AbstractGraph; step = nothing, nodelabel::Vector{String} = map(string, vertices(g)), nodeop::Vector{String} = map(string, vertices(g)), nodedevice::Vector{String} = fill("cpu", nv(g)), nodevalue::Vector{Any} = fill(nothing, nv(g)))
	nv(g) == length(nodelabel)  || throw(ArgumentError("length of nodelable must be same as number of vertices"))
	nv(g) == length(nodeop)  || throw(ArgumentError("length of nodeop must be same as number of vertices"))
	nv(g) == length(nodedevice)  || throw(ArgumentError("length of nodedevice must be same as number of vertices"))
	nv(g) == length(nodevalue)  || throw(ArgumentError("length of nodevalue must be same as number of vertices"))
	summ = SummaryCollection(graph_summary(g, nodelabel, nodeop, nodedevice, nodevalue))
    write_event(logger.file, make_event(logger, summ, step=step))
end

function graph_summary(g, nodelabel, nodeop, nodedevice, nodevalue)
	nodes = Vector{NodeDef}()
	for v in vertices(g)
		name = nodelabel[v]
		op = nodeop[v]
		input = [nodelabel[x] for x in inneighbors(g, v)]
		device = nodedevice[v]
		attr = Dict{String, AttrValue}()
		x = nodevalue[v]
		if isa(x, AbstractString)
			attr["value"] = AttrValue(s = Vector{UInt8}(x))
			attr["dtype"] = AttrValue(_type = jltype2tf(typeof(x)))
		elseif isa(x, Integer)
			attr["value"] = AttrValue(i = Int64(x))
			attr["dtype"] = AttrValue(_type = jltype2tf(typeof(x)))
		elseif isa(x, Real)
			attr["value"] = AttrValue(f = Float32(x))
			attr["dtype"] = AttrValue(_type = jltype2tf(typeof(x)))
		elseif isa(x, Bool)
			attr["value"] = AttrValue(b = x)
			attr["dtype"] = AttrValue(_type = jltype2tf(typeof(x)))
		elseif isa(x, AbstractArray)
			shape = TensorShapeProto(dim = [TensorShapeProto_Dim(size = d) for d in (collect(size(x)))])
			t = TensorProto(dtype = jltype2tf(eltype(x)), tensor_shape = shape, tensor_content = serialize_proto(string(x)))
			attr["value"] = AttrValue(tensor = t)
			attr["_output_shapes"] = AttrValue(list = AttrValue_ListValue(shape = [shape]))
		elseif isa(x, Tuple)
			attr["value"] = AttrValue(list = AttrValue_ListValue(s = [Vector{UInt8}(repr(y)) for y in x]))
			shape = TensorShapeProto(dim = [TensorShapeProto_Dim(size = length(x))])
			attr["_output_shapes"] = AttrValue(list = AttrValue_ListValue(shape = [shape]))
		elseif isa(x, Function)
			attr["value"] = AttrValue(func = NameAttrList(name = repr(x)))
		else
			#donothing
		end
		node = NodeDef(name = name, op = op, input = input, device = device, attr = attr)
		push!(nodes, node)
	end
	GraphDef(node = nodes)
end

function jltype2tf(dtype::DataType)
	nodetype =
	dtype == UInt8				? _DataType.DT_UINT8 :
	dtype == UInt16				? _DataType.DT_UINT16 :
	dtype == UInt32				? _DataType.DT_UINT32 :
	dtype == UInt64				? _DataType.DT_UINT64 :
	dtype == Int8				? _DataType.DT_INT8 :
	dtype == Int16				? _DataType.DT_INT16 :
	dtype == Int32				? _DataType.DT_INT32 :
	dtype == Int64				? _DataType.DT_INT64 :
	dtype == Float16			? _DataType.DT_BFLOAT16 :
	dtype == Float32			? _DataType.DT_FLOAT :
	dtype == Float64			? _DataType.DT_DOUBLE :
	dtype <: AbstractString		? _DataType.DT_STRING :
	dtype == Bool				? _DataType.DT_BOOL :
	dtype ∈ [Complex{Float32},
			Complex{Float16},
			Complex{UInt8},
			Complex{UInt16},
			Complex{UInt32},
			Complex{Int8},
			Complex{Int16},
			Complex{Int32}]		? _DataType.DT_COMPLEX64 :
	dtype ∈ [Complex{Float64},
			Complex{UInt64},
			Complex{Int64}]		? _DataType.DT_COMPLEX128 :
	_DataType.DT_INVALID
end


text:
"""
    traits(m::Module, T::Assignable)

Returns a set of traits that the data type `T` exhibits, including
the ones that were assigned to any supertypes of `T`.
See also [`@assign`](@ref).
"""
function traits(m::Module, T::Assignable)
    traits_map = get_traits_map(m)
    base = Set{DataType}()
    for (Tmap, s) in pairs(traits_map)
        if T <: Tmap
            union!(base, s)
        end
    end
    return base
end

"""
    assign(m::Module, T::Assignable, trait::DataType)

Assign data type `T` with the specified trait.
"""
function assign(m::Module, T::Assignable, trait::DataType)
    push_traits_map!(m, T, trait)
    return nothing
end

"""
    @assign <T> with <trait1, trait2, ...>

Assign traits to the data type `T`.  The traits may be
positive or negative e.g. `Can{Fly}` or `Cannot{Swim}`.
"""
macro assign(T::Union{Expr,Symbol}, with::Symbol, traits::Union{Expr,Symbol})
    usage = "Invalid usage: try something like `@assign Duck with Can{Fly}`"
    with === :with || throw(SyntaxError(usage))
    mod = __module__

    assign_impl(mod, T, traits)
end

function assign_impl(mod, T, traits)

    # There are two possible AST cases:
    # 1. @assign T with Can{X}
    # 2. @assign T with Can{X},Can{Y},...
    #
    # The first case would be an Expr with head == :curly
    # The second case would be an Expr with head == :tuple and sub-expr of :curly's
    capabilities = if traits isa Expr && traits.head == :tuple
        [(side = x.args[1], trait = x.args[2]) for x in traits.args]
    elseif traits isa Expr && traits.head == :curly
        x = traits
        [(side = x.args[1], trait = x.args[2])]
    else
        throw(SyntaxError("Must assign types with trait can/cannot-types"))
    end

    # Build up expressions for each capability
    expressions = Expr[]
    for cap in capabilities
        push!(expressions,
                quote
                    BinaryTraits.trait(::Type{$(cap.trait)}, ::Type{<:$T}) =
                        $(cap.side){$(cap.trait)}()
                    BinaryTraits.assign($mod, $T, $(cap.side){$(cap.trait)})
                end)
    end

    expr = quote
        $(expressions...)
    end

    display_expanded_code(expr)
    return esc(expr)
end


text:
#
# Check if structure is in mmCIF format and returns the list of
# desired fields if that is the case
#
function check_mmCIF(file)

    local ifield, mmCIF, line

    mmCIF_fields = Indexes_mmCIF_fields()

    pdbfile = open(file, "r")

    mmCIF = false

    while !eof(pdbfile)
        try
            line = strip(readline(pdbfile))
            if length(line) < 5
                continue
            end
        catch
            break
        end
        if line[1:5] == "loop_"
            mmCIF = true
            ifield = 0
            while !eof(pdbfile)
                try
                    line = strip(readline(pdbfile))
                catch
                    break
                end
                if line[1:1] == "#"
                    continue
                end
                if ifield > 0
                    ifield = ifield + 1
                    if occursin("_atom_site.label_atom_id", line)
                        mmCIF_fields.name = ifield
                    end
                    if occursin("_atom_site.id", line)
                        mmCIF_fields.index = ifield
                    end
                    if occursin("_atom_site.label_comp_id", line)
                        mmCIF_fields.resname = ifield
                    end
                    if occursin("_atom_site.label_asym_id", line)
                        mmCIF_fields.chain = ifield
                    end
                    if occursin("_atom_site.label_seq_id", line)
                        mmCIF_fields.resnum = ifield
                    end
                    if occursin("_atom_site.Cartn_x", line)
                        mmCIF_fields.x = ifield
                    end
                    if occursin("_atom_site.Cartn_y", line)
                        mmCIF_fields.y = ifield
                    end
                    if occursin("_atom_site.Cartn_z", line)
                        mmCIF_fields.z = ifield
                    end
                    if occursin("_atom_site.B_iso_or_equiv", line)
                        mmCIF_fields.beta = ifield
                    end
                    if occursin("_atom_site.occupancy", line)
                        mmCIF_fields.occup = ifield
                    end
                    if occursin("ATOM", line) && occursin("HETATOM", line)
                        break
                    end
                end
                if occursin("_atom_site.group_PDB", line)
                    ifield = 1
                end
            end
        end
    end
    close(pdbfile)

    if mmCIF
        if mmCIF_fields.name == 0 ||
           mmCIF_fields.resname == 0 ||
           mmCIF_fields.chain == 0 ||
           mmCIF_fields.resnum == 0 ||
           mmCIF_fields.x == 0 ||
           mmCIF_fields.y == 0 ||
           mmCIF_fields.z == 0
            error(" ERROR: Could not find all necessary fields in mmCIF file. ")
        end
    end

    return mmCIF, mmCIF_fields

end

max_stars_repo_name: mcubeg/PDBTools

max_issues_repo_name: mcubeg/PDBTools

max_forks_repo_name: mcubeg/PDBTools


text:
push!(LOAD_PATH, "./src/")

using AdaGram
using AdaGram.ArgParse
using Distributed

s = ArgParseSettings()

@add_arg_table s begin
  "train"
    help = "training text data"
    arg_type = AbstractString
    required = true
  "dict"
    help = "dictionary file with word frequencies"
    arg_type = AbstractString
    required = true
  "output"
    help = "file to save the model (in Julia format)"
    arg_type = AbstractString
    required = true
  "--window"
    help = "(max) window size"
    arg_type = Int
    default = 4
  "--workers"
    help = "number of workers for parallel training"
    arg_type = Int
    default = 1
  "--min-freq"
    help = "min. frequency of the word"
    arg_type = Int
    default = 20
  "--remove-top-k"
    help = "remove top K most frequent words"
    arg_type = Int
    default = 0
  "--dim"
    help = "dimensionality of representations"
    arg_type = Int
    default = 100
  "--prototypes"
    help = "number of word prototypes"
    arg_type = Int
    default = 5
  "--alpha"
    help = "prior probability of allocating a new prototype"
    arg_type = Float64
    default = 0.1
  "--d"
    help = "parameter of Pitman-Yor process"
    arg_type = Float64
    default = 0.
  "--subsample"
    help = "subsampling treshold. useful value is 1e-5"
    arg_type = Float64
    default = Inf
  "--context-cut"
    help = "randomly reduce size of the context"
    action = :store_true
  "--epochs"
    help = "number of epochs to train"
    arg_type = Int64
    default = 1
  "--init-count"
    help = "initial weight (count) on first sense for each word"
    arg_type = Float64
    default = -1.
  "--stopwords"
    help = "file with list of stop words"
    arg_type = AbstractString
  "--sense-treshold"
    help = "minimal probability of a meaning to contribute into gradients"
    arg_type = Float64
    default = 1e-10
  "--save-treshold"
    help = "minimal probability of a meaning to save after training"
    arg_type = Float64
    default = 1e-3
  "--regex"
    help = "ignore words not matching provided regex"
    arg_type = AbstractString
    default = ""
end

args = parse_args(ARGS, s)

addprocs(args["workers"])
@everywhere using AdaGram

stopwords = Set{AbstractString}()
if args["stopwords"] != nothing
  stopwords = Set{AbstractString}(readdlm(args["stopwords"]))
end

print("Building dictionary... ")
vm, dict = read_from_file(args["dict"], args["dim"], args["prototypes"],
  args["alpha"], args["d"], args["min-freq"], args["remove-top-k"], stopwords;
  regex=Regex(args["regex"]))
println("Done!")

window = args["window"]

inplace_train_vectors!(vm, dict, args["train"], window;
  threshold=args["subsample"], context_cut=args["context-cut"],
  epochs=args["epochs"], init_count=args["init-count"], sense_treshold=args["sense-treshold"])

save_model(args["output"], vm, dict, args["save-treshold"])


text:
@testset "565" begin
    @test array_nesting([5, 4, 0, 3, 1, 6, 2]) == 4
    @test array_nesting([30, 15, 49, 34, 48, 9, 11, 25, 44, 47, 35, 41, 20, 23, 22, 37, 1,
                         18, 10, 5, 3, 43, 26, 4, 42, 6, 39, 0, 36, 21, 16, 50, 38, 45, 13,
                         46, 17, 40, 28, 7, 24, 12, 27, 14, 33, 31, 2, 8, 29, 19, 32]) == 21
end


text:
@inline soft_th(x, ϵ) = max(x-ϵ,zero(x)) + min(x+ϵ,zero(x))
@inline soft_th(x, ϵ, l) = max(x-ϵ,l) + min(x+ϵ,l) - l
@inline function soft_th(x::Complex, ϵ)
    m,a = abs(x), angle(x)
    m = max(m-ϵ,zero(m)) + min(m+ϵ,zero(m))
    m*cis(a)
end

function soft_hankel!(A, ϵ)
    K,L = size(A)
    N = K+L-1
    for k = 1:N
        ri = min(K,k):-1:max(k-L,1)
        ci = max(1,k-K+1):L
        m = mean(A[r,c] for (r,c) in zip(ri,ci))
        for (r,c) in zip(ri,ci)
            A[r,c] = soft_th(A[r,c], ϵ, m)
        end
    end
    A
end

"""
    unhankel(A)

The inverse of [`hankel`](@ref). Create a 1-D signal by antidiagonal averaging
"""
function unhankel(A)
    K,L = size(A)
    N = L+(K-1)
    y = similar(A, N)
    for k = 1:N
        ri = min(K,k):-1:max(k-L,1)
        ci = max(1,k-K+1):L
        m = mean(A[r,c] for (r,c) in zip(ri,ci))
        y[k] = m
    end
    y
end

"""
    unhankel(A, lag, N, D=1)

The inverse of [`hankel`](@ref). Create a signal of the original dimension by antidiagonal averaging

# Arguments:
- `A`: A Hankel matrix
- `lag`: if lag was used to create `A`, you must provide it to `unhankel`
- `N`: length of the original signal
- `D`: dimension of the original signal
"""
function unhankel(A,lag,N,D=1)
    lag == 1 && D == 1 && (return unhankel(A))
    K      = size(A,1)
    L      = size(A,2)÷D
    y      = zeros(eltype(A), N, D)
    counts = zeros(Int, N, D)
    indmat = CartesianIndex.(1:N, (1:D)')
    inds   = hankel(indmat, L, lag)
    for (Aind,yind) in enumerate(inds)
        y[yind] += A[Aind]
        counts[yind] += 1
    end
    y ./= max.(counts, 1)
    D == 1 && (return vec(y))
    y
end

"""
    X = hankel(x,L,lag=1)

Form a hankel "trajectory matrix" `X` of size KxL, K = N-L+1
x can be a vector or a matrix.
"""
function hankel(x,L,lag=1)
    N = size(x,1)
    D = size(x,2)
    @assert L <= N/2 "L has to be less than N/2 = $(N/2)"
    @assert lag <= L "lag must be <= L"
    K = (N-L)÷lag+1
    X = similar(x, K, L*D)
    for d in 1:D
        inds = 1:L
        colinds = d:D:L*D
        for k = 1:K
            X[k,colinds] = x[inds,d]
            inds = inds .+ lag
        end
    end
    X
end

function ishankel(A)
    K,L = size(A)
    N = K+L-1
    for k = 1:N
        ri = min(K,k):-1:max(k-L,1)
        ci = max(1,k-K+1):L
        val = A[ri[1],ci[1]]
        for (r,c) in zip(ri,ci)
            A[r,c] != val && return false
        end
    end
    true
end

"""
    yf = lowrankfilter(y, n=min(length(y) ÷ 20, 2000); tol=1e-3, kwargs...)

Filter time series `y` by forming a lag-embedding T (a Toeplitz matrix) and using [`rpca`](@ref) to recover a low-rank matrix from which the a filtered signal `yf` can be extracted. The size of the embedding `n` determines the complexity, higher `n` generally gives better filtering at the cost of roughly cubic complexity.

# Arguments:
- `y`: A signal to be filtered, assumed corrupted with sparse noise
- `n`: Embedding size
- `kwargs`: See [`rpca`](@ref) for keyword arguments.
"""
function lowrankfilter(y, n=min(size(y,1)÷20,2000); sv=0, lag=1, tol=1e-3, svd = svd!, kwargs...) where {F <: Function}
    H = hankel(y, n, lag)
    if sv <= 0
        A,E = rpca(H; tol=tol, svd = svd, kwargs...)
    else
        s = svd(H)
        A = s.U[:,1:sv] * Diagonal(s.S[1:sv]) * s.Vt[1:sv,:]
    end
    unhankel(A, lag, size(y,1), size(y,2))
end

"""
    A,E,s,sv = rpca(D::Matrix; λ=1.0 / √(maximum(size(D))), iters=1000, tol=1.0e-7, ρ=1.5, verbose=false, nonnegA=false, nonnegE=false, nukeA=true)

minimize_{A,E} ||A||_* + λ||E||₁ s.t. D = A+E

`s` is the last calculated svd of `A` and `sv` is the estimated rank.

Ref: "The Augmented Lagrange Multiplier Method for Exact Recovery of Corrupted Low-Rank Matrices", Zhouchen Lin, Minming Chen, Leqin Wu, Yi Ma, https://people.eecs.berkeley.edu/~yima/psfile/Lin09-MP.pdf
Significant inspiration taken from an early implementation by Ryuichi Yamamoto in RobustPCA.jl

# Arguments:
- `D`: Design matrix
- `λ`: Sparsity regularization
- `maxrank`: Upper limit on the rank of the estimated matrix. Setting a smaller value makes convergence faster.
- `iters`: Maximum number of iterations
- `tol`: Tolerance
- `ρ`: Algorithm tuning param
- `verbose`: Print status
- `nonnegA`: Hard thresholding on A
- `nonnegE`: Hard thresholding on E
- `nukeA`: Activate the nuclear penalty on `A`, if `false`, then `A` is not assumed to be low rank.
- `hankel`: Indicating whether or not `D` (and thus `A` and `E`) are Hankel matrices (constant anti diagonals). If this fact is known, the expected performance of this alogorithm goes up. If the matrix `D` is Toeplitz (constant diagonals) you may reverse the second dimension, i.e., `Dᵣ = D[:,end:-1:1]`. `hankel=true` should likely be paired with `nukeA=false`.

To speed up convergence you may either increase the tolerance or increase `ρ`. Increasing `tol` is often the best solution.
"""
function rpca(D::AbstractMatrix{T};
                          λ              = real(T)(1.0/sqrt(maximum(size(D)))),
                          maxrank        = typemax(Int),
                          iters::Int     = 1000,
                          tol            = sqrt(eps(real(T))),
                          ρ              = real(T)(1.5),
                          verbose::Bool  = false,
                          nonnegA::Bool  = false,
                          nonnegE::Bool  = false,
                          hankel::Bool   = false,
                          # proxE          = NormL1(λ),
                          nukeA          = true,
                          svd::F1        = svd!,
                          opnorm::F2     = opnorm,
                          kwargs...) where {F1 <: Function, F2 <: Function, T}
    RT        = real(T)
    M, N      = size(D)
    d         = min(M,N)
    A, E      = zeros(T, M, N), zeros(T, M, N)
    Z         = similar(D)
    Y         = copy(D)
    norm²     = opnorm(Y)::RT # can be tuned # norm(Y)/minimum(size(D))^(2/5)
    norm∞     = norm(Y, Inf) / λ
    dual_norm = max(norm², norm∞)
    d_norm    = norm²
    Y       ./= dual_norm
    μ         = RT(1.25 / norm²)
    μ̄         = RT(μ * 1.0e+7)
    sv = svp  = 10
    local s
    for k = 1:iters
        # prox!(E, proxE, D .- A .+ (1/μ) .* Y, 1/μ)
        E .= soft_th.(D .- A .+ (1/μ) .* Y, λ/μ)
        if nonnegE
            E .= max.(E, 0)
        end
        Z .= D .- E .+ (1/μ) .* Y
        if svd ∈ (LinearAlgebra.svd, LinearAlgebra.svd!) || k == 1
            s = LinearAlgebra.svd!(Z) # Z assignment just for storage
        else
            s = svd(Z, sv)
        end
        svp = sum(x-> x >= (1/μ), s.S)::Int
        if svp < sv
            sv = svp#min(svp + 1, N) # the paper says to use these formulas but sv=svp works way better
        else
            sv = svp#min(svp + round(Int, T(0.05) * d), d)
        end
        sv = clamp(sv, 1, maxrank)
        @views if nukeA
            # A .= s.U[:,1:svp] * Diagonal(s.S[1:svp] .- 1/μ) * s.Vt[1:svp,:]
            mul!(Z[:,1:svp], s.U[:,1:svp], Diagonal(s.S[1:svp] .- 1/μ))
            mul!(A, Z[:,1:svp], s.Vt[1:svp,:])
        else
            # A .= s.U[:,1:svp] * Diagonal(s.S[1:svp]) * s.Vt[1:svp,:]
            mul!(Z[:,1:svp], s.U[:,1:svp], Diagonal(s.S[1:svp]))
            mul!(A, Z[:,1:svp], s.Vt[1:svp,:])
        end
        if hankel
            soft_hankel!(A, λ/μ)
        end
        if nonnegA
            A .= max.(A, 0)
        end

        @. Z = D - A - E # Z are the reconstruction errors
        @. Y = Y + μ * Z # Y can not be moved below as it depends on μ which is changed below
        μ = RT(min(μ*ρ, μ̄))

        cost = opnorm(Z) / d_norm
        verbose && println("$(k) cost: $(round(cost, sigdigits=4))")

        if cost < tol
            verbose && println("converged")
            break
        end
        k == iters && @warn "Maximum number of iterations reached, cost: $cost, tol: $tol"
    end
    if hankel
        soft_hankel!(E, λ/μ)
    end

    A, E, s, sv
end

"""
    Q = rpca_ga(X::AbstractMatrix{T}, r=minimum(size(X)), U=similar(X); μ = μ!(s,w,U), verbose=false, kwargs...) where T

"Grassmann Averages for Scalable Robust PCA", Hauberg et al.
http://www2.compute.dtu.dk/~sohau/papers/cvpr2014a/Hauberg_CVPR_2014.pdf

# Arguments:
- `X`: Data matrix
- `r`: Rank (number of components to estimate)
- `U`: Optional pre-allocated buffer
- `verbose`: print status
- `kwargs`: such as `tol=1e-7`, `iters=1000`
- `μ = μ!(s,w,U)` is a function that calculates the spherical average of a all columns in the matrix `U`, weighted by `w` and stores the result in `s`. The default is the standard weighted average. To get a robust estimate, consider using a robust average, such as `entrywise_trimmed_mean` or `entrywise_median` etc.
"""
function rpca_ga(X::AbstractMatrix{T}, r=minimum(size(X)), U = similar(X); verbose = false, kwargs...) where T
    d,N    = size(X)
    X      = copy(X)
    Xs1    = similar(X,1,N)
    Xs2    = similar(X)
    Q      = zeros(T, d, r)
    w      = zeros(T,N)
    Xnorms = zeros(T,N)
    for i = 1:r
        @inbounds @views for n = 1:N
            Xnorms[n] = sqrt(sum(abs2,X[:,n]))
            U[:,n] .=  X[:,n] ./ Xnorms[n]
        end
        q = rpca_ga_1(Xnorms, U, w; verbose = verbose, kwargs...)
        Q[:,i] .= q
        @static if VERSION >= v"1.3"
            mul!(Xs1, q', X)
            mul!(X, q,Xs1, -1, 1)
        else
            X .-= q*(q'X)
        end
    end
    Q
end

"""
Find the first principal component. This is an internal function used by `rpca_ga`.
"""
function rpca_ga_1(Xnorms,U::AbstractMatrix{T},w; tol=1e-7, iters=1000, verbose=false, μ=μ!) where {T}
    d,N = size(U)

    q = randn(d)
    q ./= norm(q)
    qold = copy(q)

    for i = 1:iters
        @inbounds @views for n = eachindex(w,Xnorms)
            w[n] = sign(U[:,n]'q)*Xnorms[n]
        end
        μᵢ = μ(q,w,U)
        q .= μᵢ./norm(μᵢ)
        dq = sqrt(sum(((x,y),)->abs2(x-y), zip(q,qold)))
        verbose && @info "Change at iteration $i: $dq"
        if dq < tol
            verbose && @info "Converged after $i iterations"
            break
        end
        qold .= q
        i == iters && @warn "Reached maximum number of iterations"
    end
    q
end

function μ!(s,w,U)
    ws = zero(eltype(w))
    s .= 0
    @views @inbounds @simd for n = 1:size(U,2)
        ws += w[n]
        s .+= w[n].*U[:,n]
    end
    s ./= ws
end

"""
    entrywise_trimmed_mean(s, w, U, P=0.1)

Remove `P` percent of the data on each side before computing the weighted mean.
"""
function entrywise_trimmed_mean(s,w,U, P=0.1)
    N = size(U,2)
    range = (1+floor(Int, P*N)):floor(Int, (1-P)*N)
    s .= 0
    @views for j = 1:size(U,1)
        I = sortperm(U[j,:])[range]
        s[j] += w[I]'U[j,I] / sum(w[I])
        # s[j] += sum(U[j,I])/length(I)
    end
    s
end
#
#
# function entrywise_trimmed_mean(s,w,U, P=0.05)
#     d,N = size(U)
#     range = (1+floor(Int, P*d)):floor(Int, (1-P)*d)
#     s .= 0
#     ws = zeros(size(s))
#     @views for j = 1:size(U,2)
#         I = sortperm(U[:,j])[range]
#         s[I] .+= w[j] .* U[I,j]
#         ws[I] .+= w[j]
#     end
#     s./ws
# end

function entrywise_median(s,w,U)
    s .= 0
    for j = 1:size(U,1)
        I = sortperm((w).*U[j,:])
        s[j] = sign(w[I[end÷2]])*U[j,I[end÷2]]#, StatsBase.Weights(abs.(w)))
        # s[j] = median(U[j,:])
    end
    s
end


text:
using Statistics
import MLJFlux
import Flux

@testset  "Penalizer" begin
    A = [-1 2; -3 4]
    lambda = 1

    # 100% L2:
    alpha = 0
    penalty = MLJFlux.Penalizer(lambda, alpha)
    @test penalty(A) ≈ 1 + 4 + 9 + 16

    # 100% L1:
    alpha = 1
    penalty = MLJFlux.Penalizer(lambda, alpha)
    @test penalty(A) ≈ 1 + 2 + 3 + 4

    # no strength:
    lambda = 0
    alpha = 42.324
    penalty = MLJFlux.Penalizer(lambda, alpha)
    @test penalty(A) == 0
end

@testset "Penalty" begin
    model = MLJFlux.NeuralNetworkRegressor(lambda=1, alpha=1, loss=Flux.mae)
    chain = Flux.Dense(3, 1, identity)
    w = Flux.params(chain)
    p = MLJFlux.Penalty(model)

    # compare loss by hand and with penalized loss function:
    penalty = (sum(abs.(chain.weight)) + abs(chain.bias[1]))
    @test p(w) ≈ penalty
end


text:
"""
    PyException <: Exception

Represents an exception raised from Python.

It has three fields `tref`, `vref`, `bref` which are all `PyRef`s, and are the type, value and backtrace of the exception.
"""
mutable struct PyException <: Exception
    tref::PyRef
    vref::PyRef
    bref::PyRef
    PyException(
        ::Val{:new},
        t::Ptr = C_NULL,
        v::Ptr = C_NULL,
        b::Ptr = C_NULL,
        borrowed::Bool = false,
    ) = new(
        PyRef(Val(:new), t, borrowed),
        PyRef(Val(:new), v, borrowed),
        PyRef(Val(:new), b, borrowed),
    )
end
export PyException

pythrow() = throw(PyException(Val(:new), C.PyErr_FetchTuple(true)...))

"""
    check(x)

Checks if the Python error indicator is set. If so, throws the current Python exception. Otherwise returns `x`.
"""
check(x) = C.PyErr_IsSet() ? pythrow() : x

"""
    checkm1(x, ambig=false)

Same as `check(x)` but errors are indicated by `x == -1` instead.

If `ambig` is true the error indicator is also checked.
"""
checkm1(x::Number, ambig::Bool = false) = ism1(x) ? ambig ? check(x) : pythrow() : x

"""
    checknull(x, ambig=false)

Same as `check(x)` but errors are indicated by `x == C_NULL` instead.

If `ambig` is true the error indicator is also checked.
"""
checknull(x::Ptr, ambig::Bool = false) = isnull(x) ? ambig ? check(x) : pythrow() : x

"""
    checkerr(x, ambig=false)

Same as `check(x)` but errors are indicated by `x == PYERR()` instead.

If `ambig` is true the error indicator is also checked.
"""
checkerr(x, ambig::Bool = false) = x === PYERR() ? ambig ? check(x) : pythrow() : x

"""
    checknullconvert(T, x, ambig=false) :: T

Same as `checknull(x, ambig)` but steals a reference to `x` and converts the result to a `T`.
"""
checknullconvert(::Type{T}, x::Ptr, ambig::Bool = false) where {T} = begin
    if isnull(x) && (!ambig || C.PyErr_IsSet())
        C.Py_DecRef(x)
        pythrow()
    end
    r = C.PyObject_Convert(x, T)
    C.Py_DecRef(x)
    checkm1(r)
    C.takeresult(T)
end

function Base.showerror(io::IO, e::PyException)
    print(io, "Python: ")

    if isnull(e.tref)
        print(io, "mysterious error (no error was actually set)")
        return
    end

    if CONFIG.sysautolasttraceback
        try
            @py ```
            import sys
            sys.last_type = $(isnull(e.tref) ? nothing : e.tref)
            sys.last_value = $(isnull(e.vref) ? nothing : e.vref)
            sys.last_traceback = $(isnull(e.bref) ? nothing : e.bref)
            ```
        catch err
            print(io, "<error while setting 'sys.last_traceback': $err")
        end
    end

    # if this is a Julia exception then recursively print it and its stacktrace
    if C.PyErr_GivenExceptionMatches(e.tref, C.PyExc_JuliaError()) != 0
        try
            # Extract error value
            vo = @pyv `$(e.vref).args[0]`::Any
            if vo isa Tuple{Exception,Any}
                je, jb = vo
            else
                je = vo
                jb = nothing
            end
            print(io, "Julia: ")
            # Show exception
            if je isa Exception
                showerror(io, je)
            else
                print(io, je)
            end
            # Show backtrace
            if jb === nothing
                println(io)
                printstyled(io, "Stacktrace: none")
            else
                io2 = IOBuffer()
                Base.show_backtrace(
                    IOContext(io2, io),
                    jb,
                )
                printstyled(io, String(take!(io2)))
            end
            # Show Python backtrace
            if isnull(e.bref)
                println(io)
                printstyled(io, "Python stacktrace: none")
                return
            else
                @goto pystacktrace
            end
        catch err
            println("<error while printing Julia exception inside Python exception: $err>")
        end
    end

    # print the type name
    try
        tname = @pyv `$(e.tref).__name__`::String
        print(io, tname)
    catch err
        print(io, "<error while printing type: $err>")
    end

    # print the error message
    if !isnull(e.vref)
        print(io, ": ")
        try
            vstr = @pyv `str($(e.vref))`::String
            print(io, vstr)
        catch err
            print(io, "<error while printing value: $err>")
        end
    end

    # print the stacktrace
    if !isnull(e.bref)
        @label pystacktrace
        println(io)
        printstyled(io, "Python stacktrace:")
        try
            @py ```
            import traceback
            $(fs :: Vector{Tuple{String, String, Int}}) = [(x.name, x.filename, x.lineno) for x in traceback.extract_tb($(e.bref))]
            ```
            for (i, (name, fname, lineno)) in enumerate(reverse(fs))
                println(io)
                printstyled(io, " [", i, "] ")
                printstyled(io, name, bold = true)
                printstyled(io, " at ")
                # if (m=match(r"^(.*):([0-9]+)$", fname)) !== nothing
                #     fname = m.captures[1]
                #     lineno += parse(Int, m.captures[2]) - 1
                # end
                printstyled(io, fname, ":", lineno, bold = true)
            end
        catch err
            print(io, "<error while printing stacktrace: $err>")
        end
    end
end


text:
import Interact.recv_msg

function handle_msg(w::InputWidget, msg)
    if msg.content["data"]["method"] == "backbone"
        recv_msg(w, msg.content["data"]["sync_data"]["value"])
    end
end

function handle_msg{T}(w::Button{T}, msg)
    try
        if msg.content["data"]["method"] == "custom" &&
            msg.content["data"]["content"]["event"] == "click"
            # click event occured
            recv_msg(w, convert(T, w.value))
        end
    catch e
        warn(string("Couldn't handle Button message ", e))
    end
end

function handle_msg{view}(w::Options{view}, msg)
        if msg.content["data"]["method"] == "backbone"
            key = string(msg.content["data"]["sync_data"]["value"])
            if haskey(w.options, key)
                recv_msg(w, w.options[key])
            end
        end
end


text:
module Portinari

export Area, Axis, Context, D3Attr, Direction, Line, Restyle, Scale, Shape

include("./js_base.jl")

include("./d3_abstractions.jl")

include("./context.jl")

include("./area.jl")
include("./line.jl")
include("./shape.jl")
include("./axis.jl")

end # module