text:
"""
Helper function that inserts a new decision node in the search tree,
corresponding to the new time window event that is observed.
success_prob_fn takes two arguments: (i) reference time (ii) interaction event and returns a single floating
point probability of success.
"""
function insert_decision_node!(tree::SearchTree, agent_name::String, ie::InteractionEvent,
                                parent_id::Int, ref_time::Float64, utilval::Float64, success_prob_fn::Function)

    # Start with end index to add nodes
    idx = length(tree.nodes)

    # If ref_time is less than the start of the event, set true ref time to event start
    true_ref_time = (ref_time < ie.timestamps[START]) ? ie.timestamps[START] : ref_time
    @assert true_ref_time <= ie.timestamps[FINISH] "ref_time is $(ref_time) and ie is $(ie.timestamps)"

    # N.B. Decision nodes should be initialized with 0 util
    # Decision node corresponding to attempt
    idx += 1
    attempt_idx = idx
    dnode1 = DecisionNode(agent_name=agent_name, task_name=ie.task_name, attempt=true, timeval=true_ref_time,
                          util=0.0, idx=attempt_idx, utilset=false)
    push!(tree.nodes, dnode1)
    tree.parent_id[attempt_idx] = parent_id

    # Decision node corresponding to not attempt
    # Add node idx to decision leaf
    idx += 1
    no_attempt_idx = idx
    dnode2 = DecisionNode(agent_name=agent_name, task_name=ie.task_name, attempt=false, timeval=true_ref_time,
                          util=0.0, idx=no_attempt_idx, utilset=false)
    push!(tree.nodes, dnode2)
    tree.parent_id[no_attempt_idx] = parent_id
    new_decision_leaf = no_attempt_idx

    # set children of parent id to decision node idxs
    tree.child_ids[parent_id] = [attempt_idx, no_attempt_idx]

    # Use function to get probability of success
    # Arguments of success_prob are the ref_time and the time window
    success_prob = success_prob_fn(ref_time, ie)

    # For now, set util to -Inf because leaf will set it
    idx += 1
    fail_idx = idx
    onode1 = OutcomeNode(agent_name=agent_name, task_name=ie.task_name, outcome=FINISH,
                         timeval=ie.timestamps[FINISH], probability=1.0 - success_prob, util=0.0, idx=fail_idx, utilset=false)
    push!(tree.nodes, onode1)
    tree.parent_id[fail_idx] = attempt_idx

    idx += 1
    succ_idx = idx
    onode2 = OutcomeNode(agent_name=agent_name, task_name=ie.task_name, outcome=SUCCESS,
                         timeval=ie.timestamps[SUCCESS], probability=success_prob, util=utilval, idx=succ_idx, utilset=false)
    push!(tree.nodes, onode2)
    tree.parent_id[succ_idx] = attempt_idx

    # Set child IDs of attempt_idx and new outcome leaves to succ and fail
    tree.child_ids[attempt_idx] = [succ_idx, fail_idx]
    new_outcome_leaves = [succ_idx, fail_idx]

    return new_outcome_leaves, new_decision_leaf
end

"""
Generate the allocation policy tree for a given agent and set of tasks to consider.
The util_val_fn returns the utility of completing an interaction event successfully.
"""
function generate_search_tree!(server::GenericAllocation{A,AP}, agent_name::String,
                               tasks_to_consider::Set{String}, success_prob_fn::Function, util_val_fn::Function,
                               downtime::Float64) where {A,AP}

    # Must be at least one task to consider
    if isempty(tasks_to_consider)
        empty!(server.agent_prop_set[agent_name].tree)
        return
    end

    # Assume that updated time windows have been generated
    # and that agent is currently not occupied with task already
    # Have outcome leaf idxs and outcome decision idxs (when NOT ATTEMPTED)
    # If a new event subsumes all outcome leaf idxs it should NOT be added to not-attempted OUTCOME DECISION IDX
    tree = server.agent_prop_set[agent_name].tree
    interaction_events = server.agent_prop_set[agent_name].interaction_events

    outcome_leaf_idxs = Set{Int}()
    decision_leaf_idxs = Set{Int}()

    first_ie_stamp = Inf

    considered_ie = [ie for ie in interaction_events if ie.task_name in tasks_to_consider]

    num_to_consider = (length(considered_ie) > server.max_tasks_to_consider) ? server.max_tasks_to_consider : length(considered_ie)

    # IMP - The first interaction event timestamp may be shared by multiple
    # events; in that case, should consider all at root
    for ie in considered_ie[1:num_to_consider]

        ie_stamp = ie.timestamps[START]

        if ie_stamp < first_ie_stamp

            utilval = util_val_fn(ie)
            (new_outcome_leaves, new_decision_leaf) = insert_decision_node!(tree, agent_name, ie,
                                                                            0, server.current_time,
                                                                            utilval, success_prob_fn)
            for idx in new_outcome_leaves
                push!(outcome_leaf_idxs, idx)
            end
            push!(decision_leaf_idxs, new_decision_leaf)

            first_ie_stamp = ie_stamp

        else

            # Now iterate through outcome leaf nodes and compare
            new_outcome_leaf_idxs = Set{Int}()
            new_decision_leaf_idxs = Set{Int}()
            outcome_leaves_to_rm = Set{Int}()
            decision_leaves_to_rm = Set{Int}()

            for ol_idx in outcome_leaf_idxs
                onode = tree.nodes[ol_idx]
                @assert typeof(onode) == OutcomeNode # sanity check

                if onode.timeval <= ie.timestamps[START]
                    # Dominated - continue
                    continue
                end

                if onode.timeval >= ie.timestamps[FINISH] - downtime/2.0
                    continue
                end

                # Last case - new IE begins before outcome and ends after outcome
                # now add decision node to try and its two outcomes
                utilval = util_val_fn(ie)
                (temp_outcome_idxs, temp_dec_idx) = insert_decision_node!(tree, agent_name, ie,
                                                                          ol_idx, onode.timeval,
                                                                          utilval, success_prob_fn)
                push!(outcome_leaves_to_rm, ol_idx)

                for to_idx in temp_outcome_idxs
                    push!(new_outcome_leaf_idxs, to_idx)
                end
                push!(new_decision_leaf_idxs, temp_dec_idx)
            end

            # IMP - Not attempting is DIFFERENT from failing immediately, because it affects the dominance principle
            # Current no-attempt decision nodes MUST have had their attempt counterparts added already
            # So if new event is dominated, then it should not follow no-attempt leaf nodes, EVEN if otherwise valid
            @assert isempty(new_decision_leaf_idxs) == isempty(new_outcome_leaf_idxs)

            # Otherwise, loop through decision_leaf_idxs and do the same thing
            for dec_idx in decision_leaf_idxs
                dnode = tree.nodes[dec_idx]

                @assert dnode.attempt == false

                # We know that event is non-dominated by outcomes
                # and we did not jump forward in time, so do not need to do earlier checks
                if dnode.timeval >= ie.timestamps[FINISH] - downtime/2.0
                    continue
                end

                utilval = util_val_fn(ie)
                (temp_outcome_idxs, temp_dec_idx) = insert_decision_node!(tree, agent_name, ie,
                                                                        dec_idx, dnode.timeval,
                                                                        utilval, success_prob_fn)
                push!(decision_leaves_to_rm, dec_idx)

                for to_idx in temp_outcome_idxs
                    push!(new_outcome_leaf_idxs, to_idx)
                end
                push!(new_decision_leaf_idxs, temp_dec_idx)
            end

            # Finally, remove old leaves and add new ones
            setdiff!(outcome_leaf_idxs, outcome_leaves_to_rm)
            union!(outcome_leaf_idxs, new_outcome_leaf_idxs)

            setdiff!(decision_leaf_idxs, decision_leaves_to_rm)
            union!(decision_leaf_idxs, new_decision_leaf_idxs)
        end
    end

    # Now go up from leaves, assigning utilities to parent nodes
    # First assign utilities to actual leaves based on outcomes
    # Then, add parents to the current fringe
    rev_fringe_idxs = Set{Int64}()

    for ol_idx in outcome_leaf_idxs

        # Add to util of parent node
        ol_node = tree.nodes[ol_idx]
        ol_node.utilset = true

        push!(rev_fringe_idxs, ol_idx)
    end

    # Do the same for decision leaves, but now
    # Since there are only positive or zero utilities
    # can just do max with self
    for dec_idx in decision_leaf_idxs

        # this must be a non-attempt decision so util zero
        dec_node = tree.nodes[dec_idx]
        @assert dec_node.attempt == false
        dec_node.util = 0.0 # For good measure
        dec_node.utilset = true
        push!(rev_fringe_idxs, dec_idx)
    end

    # IMP - REV FRINGE IDXS can have their own children in rev_fringe_idxs too!
    # So keep adding and removing until empty
    # N.B - This is for setting PARENT UTIL
    while !isempty(rev_fringe_idxs)

        to_rm = Set{Int64}()
        to_add = Set{Int64}()

        for rfi in rev_fringe_idxs

            # If child of root, should remove (unless still needed)
            if tree.parent_id[rfi] == 0
                push!(to_rm, rfi)
                continue
            end

            node = tree.nodes[rfi]
            @assert node.utilset == true "Node in rev fringe does NOT have util set",node

            # If any sibling not set, continue WHILE keeping in rev_fringe_idxs
            par_id = tree.parent_id[rfi]
            par_node = tree.nodes[par_id]

            siblingset = true
            for child in tree.child_ids[par_id]
                if ~(child in rev_fringe_idxs)
                    siblingset = false
                    break
                end
            end
            if siblingset == false
                continue
            end

            # Do what is necessary for parent and add to rfi
            # Parent is decision node
            if typeof(node) == OutcomeNode
                par_node.util = par_node.util + node.probability*node.util
            else
                # Node type is decision
                # If parent is decision, just compare child decisions and choose max
                # If parent is outcome, add to exp. util if outcome success
                temp_node_util = node.util
                if typeof(par_node) == OutcomeNode && par_node.outcome == SUCCESS
                    temp_node_util = 1.0 + node.util
                end
                par_node.util = (temp_node_util > par_node.util) ? temp_node_util : par_node.util

            end

            push!(to_add, par_id)
            push!(to_rm, rfi)
        end

        for par_id in to_add
            tree.nodes[par_id].utilset = true
        end

        # First remove and then add
        union!(rev_fringe_idxs, to_add)
        setdiff!(rev_fringe_idxs, to_rm)
    end
end

# When there are multiple simultaneous objects, then it is not enough
# that the 'first' object is not attempted
function get_next_attempt_idx(tree::SearchTree)

    if isempty(tree)
        return -1
    end

    # Find the first node at which ATTEMPTING is the decision
    # Until that point, put all objs in entries
    dec_fringe = Deque{Int64}()
    push!(dec_fringe, 0)

    subtree_root_dec = -1

    while !isempty(dec_fringe)

        dec_top = popfirst!(dec_fringe)

        top_children = tree.child_ids[dec_top]

        @assert length(top_children) == 2 # TODO : Check this is valid
        attempt_idx = (tree.nodes[top_children[1]].attempt == true) ? top_children[1] : top_children[2]
        no_attempt_idx = top_children[1] + top_children[2] - attempt_idx

        if tree.nodes[no_attempt_idx].util > tree.nodes[attempt_idx].util
            push!(dec_fringe, no_attempt_idx)
        else
            subtree_root_dec = attempt_idx
        end
    end

    return subtree_root_dec
end


text:
using JSON

args = []
config = JSON.parsefile("config.json")

for option in config["cli"]
    if !isempty(option[2])
        if typeof(option[2])==String
            if option[2]=="true"
                push!(args,"-$(option[1])")
            elseif option[2]!="false"
                push!(args,"-$(option[1])", "$(option[2])")
            end
        else
            push!(args,"-$(option[1])")
            for params in option[2]
                if !isempty(params[2])
                    push!(args,"$(params[1])=$(params[2])")
                end
            end
        end
    end
end

sfpath = "engines/stockfish-11-popcnt/stockfish_20011801_x64_modern.exe"
lc0path = "engines/lc0-cudnn/lc0.exe"

push!(args,"-engine","name=Stockfish-11","cmd=$sfpath")
push!(args,"option.Threads=$(config["Threads"])")
push!(args,"option.SyzygyPath=$(config["tablebases"])")

for net in readdir("nets/")
    push!(args,"-engine","name=Lc0-$net","cmd=$lc0path")
    push!(args,"option.WeightsFile=$(pwd())\\nets\\$net")
    push!(args,"option.Backend=cudnn-auto")
    push!(args,"option.Threads=2")
    push!(args,"option.SyzygyPath=$(config["tablebases"])")

end

# print(args)
run(`cutechess-cli $args`)

ext: jl


text:
"""
    minimumKnapsackProblem(C, w, W)

Solve minimum knapsack problem using vector of costs `C`, weights `w` and overall weight limit `W`.
"""
function minimumKnapsackProblem(C, w, W)
    @assert length(C) == length(w)

    n = length(C)
    model = Model(solver = CplexSolver(CPXPARAM_ScreenOutput = getProperty("minimumKnapsackProblem.cplexLogging")))
    @variable(model, x[1:n], Bin)
    @objective(model, Min, dot(C, x))

    for constraint in getKnapsackConstraints(w, W)
        JuMP.addconstraint(model, constraint(x))
	end

    status = solve(model)
    return map((x) -> Int(x), getvalue(x))
end

"""
    getKnapsackConstraints(w, W)

Return a list of constraints defining a set of feasible solutions of a minimum knapsack problem.
Each constraint is function with one parameter, which is variable of a mathematical programming model.
"""
function getKnapsackConstraints(w, W)
    [x -> @LinearConstraint(vecdot(w, x) >= W)]
end


text:
using Requires

@require Spelling begin
    function spell_colorize(s::String, c=1.0)
        words = spellcheck_sentence(lowercase(s), c)
        sp = (w->w.offset).(words)
        ep = (w->w.offset+length(w.match)).(words)
        sar = split(s, "")
        for i in 1:length(sp)
            insert!(sar, sp[i]+2*i-2, Base.text_colors[:light_red])
            insert!(sar, ep[i]+2*i-1, Base.text_colors[:normal])
        end
        join(sar)
    end
end


text:
# module DisplayGadfly

#using GtkUtilities, ..Graphics, Colors

# import Gadfly, Compose, Gtk
import Gadfly: Plot, Aesthetics, Coord.Cartesian
import Gtk: GtkCanvas, showall
# import ..Immerse
# import ..Immerse: absolute_to_data, find_tagged, setproperty!

# export
#     Figure,
#     closefig,
#     closeall,
#     figure,
#     gcf,
#     scf,
#     render_backend,
#     set_limits!

const ppmm = 72/25.4   # pixels per mm FIXME? Get from backend? See dev2data.

struct PanZoomCallbacks
    idpzk::Culong
end

PanZoomCallbacks() = PanZoomCallbacks(0)
initialized(pzc::PanZoomCallbacks) = pzc.idpzk != 0

# Display code was copied & modified from Winston. The original
# contributors to that code included Mike Nolta, Jameson Nash,
# @slangangular, and likely others.

# While figures are by default associated with Windows (each with a
# single Canvas), you can have multiple figures per window.
mutable struct Figure <: Gtk.GtkBox
    handle
    canvas::GtkCanvas
    tb
    prepped         # tuple, a pre-processed Plot to speed rendering
    cc::Compose.Context   # fully-rendered Plot (useful for hit-testing)
    tweaks::Dict{Tuple{Symbol,Symbol},Any}  # post-rendering property changes
    panzoom_cb::PanZoomCallbacks
    figno::Int

    function Figure(c::GtkCanvas, p::Plot)
        prepped = Gadfly.render_prepare(p)
        cc = render_finish(prepped; dynamic=false)
        f = new(nothing, c, nothing, prepped, cc, copy(no_tweaks), PanZoomCallbacks())
    end
    function Figure(p::Plot)
        box, tb, c = createPlotGuiComponents()
        prepped = Gadfly.render_prepare(p)
        cc = render_finish(prepped; dynamic=false)
        f = new(box.handle, c, tb, prepped, cc, copy(no_tweaks), PanZoomCallbacks())
        Gtk.gobject_move_ref(f, box)
    end
    function Figure()
        box, tb, c = createPlotGuiComponents()
        f = new(box.handle, c, tb, nothing, Compose.context(), copy(no_tweaks), PanZoomCallbacks())
        Gtk.gobject_move_ref(f, box)
    end
end

Base.isempty(f::Figure) = f.cc == nothing

# These extract arguments of the arg-tuple that gets passed to
# Gadfly.render_prepared
_plot(prepped::Tuple) = prepped[1]
_plot(fig::Figure) = _plot(fig.prepped)
_coord(prepped::Tuple) = prepped[2]
_coord(fig::Figure) = _coord(fig.prepped)
_aes(prepped::Tuple) = prepped[3]
_aes(fig::Figure) = _aes(fig.prepped)

function setcoord!(fig::Figure, coord)
    fig.prepped = (fig.prepped[1], coord, fig.prepped[3:end]...)
    fig
end

mutable struct GadflyDisplay <: AbstractDisplay
    figs::Dict{Int,Figure}
    fig_order::Vector{Int}
    current_fig::Int
    next_fig::Int

    GadflyDisplay() = new(Dict{Int,Figure}(), Int[], 0, 1)
end

const _display = GadflyDisplay()

Base.display(d::GadflyDisplay, f::Figure) = display(f.canvas, f)

function Base.display(d::GadflyDisplay, p::Plot)
    isempty(d.figs) && figure()
    f = curfig(d)
    # Clear data that might have applied to the previous plot
    clear_hit(f)
    clear_guidata(f.canvas)
    clear_guistate!(f)
    # Supply a background, if not present
    if p.theme.background_color == nothing
        # FIXME: someday one will want to plot transparently.
        # Might be better for Gadfly to default to :auto, and then each
        # renderer can pick how to resolve that (and reserve `nothing`
        # for users who want to turn off backgrounds manually).
        p.theme.background_color = colorant"white"
    end
    # Do most of the time-consuming parts of plotting
    f.prepped = Gadfly.render_prepare(p)
    # Render in the current state
    cc = render_finish(f.prepped; dynamic=false)
    f.cc = apply_tweaks!(cc, f.tweaks)
    # Render the figure
    display(f.canvas, f)
    gcf()
end

function Base.display(d::GadflyDisplay, cc::Compose.Context)
    isempty(d.figs) && figure()
    f = curfig(d)
    # Clear data that might have applied to the previous plot
    clear_hit(f)
    clear_guidata(f.canvas)
    clear_guistate!(f)
    # Render in the current state
    f.cc = apply_tweaks!(cc, f.tweaks)
    # Render the figure
    display(f.canvas, f)
    gcf()
end

function Base.display(c::GtkCanvas, f::Figure)
    c.draw = let bad=false
        function (_)
            bad && return
            if isempty(f)
                fill!(c, RGB(1,1,1))
                return
            end
            xview = get(guidata, (c, :xview), nothing)
            if xview != nothing
                yview = guidata[c, :yview]
                set_limits!(f, xview, yview)
            end
            # Render
            backend = render_backend(c)
            try
                Compose.draw(backend, f.cc)
                guidata[c,:coords], guidata[c,:panelcoords] = backend.coords, backend.panelcoords
            catch e
                bad = true
                rethrow(e)
                # println("Immerse error: ", e.msg)
            end
            nothing
        end
    end
    Gtk.draw(c)
end

const no_tweaks = Dict{Tuple{Symbol,Symbol},Any}()

function apply_tweaks!(cc, tweaks)
    tagged = find_tagged(cc)
    for (k,v) in tweaks
        ctx = tagged[k[1]]
        setproperty!(ctx, v, k[2])
    end
    cc
end

function getproperty(figtag::Tuple{Int,Symbol}, prop::Symbol)
    figno, tag = figtag
    f = Figure(figno)
    getproperty(find_parent(f.cc, tag), prop)
end

function setproperty!(figtag::Tuple{Int,Symbol}, val, prop::Symbol)
    figno, tag = figtag
    f = Figure(figno)
    f.tweaks[(tag,prop)] = val
    apply_tweaks!(f.cc, f.tweaks)
    Gtk.draw(f.canvas)
    val
end

# Co-opt the REPL display
Base.display(::REPL.REPLDisplay, ::MIME"text/html", p::Union{Plot,Compose.Context}) = display(_display, p)

render_backend(c) = Compose.ImmerseBackend(Gtk.cairo_surface(c))

const _hit_data = Dict{Figure,Dict{Symbol,Any}}()

clear_hit(fig::Figure) = delete!(_hit_data, fig)

function addfig(d::GadflyDisplay, i::Int, fig::Figure)
    @assert !haskey(d.figs,i)
    d.figs[i] = fig
    push!(d.fig_order, i)
    while haskey(d.figs,d.next_fig)
        d.next_fig += 1
    end
    d.current_fig = i
    fig.figno = i
end

hasfig(d::GadflyDisplay, i::Int) = haskey(d.figs,i)

function switchfig(d::GadflyDisplay, i::Int)
    haskey(d.figs,i) && (d.current_fig = i)
end

function getfig(d::GadflyDisplay, i::Int)
    haskey(d.figs,i) ? d.figs[i] : error("no figure with index $i")
end
"""
`Figure(3)` gets the underlying Figure object associated with figure #3.
This can be useful if you need to layer on extra drawing on top of what
Gadfly produces.
"""
Figure(i::Int) = getfig(_display, i)

function curfig(d::GadflyDisplay)
    d.figs[d.current_fig]
end

nextfig(d::GadflyDisplay) = d.next_fig

function dropfig(d::GadflyDisplay, i::Int)
    haskey(d.figs,i) || return
    delete!(d.figs, i)
    splice!(d.fig_order, something(findfirst(isequal(i), d.fig_order), 0) )

    d.next_fig = min(d.next_fig, i)
    d.current_fig = isempty(d.fig_order) ? 0 : d.fig_order[end]
end

function initialize_toolbar_callbacks(f::Figure)
    c = f.canvas
    # Work around Gtk #161 & #185
    # Gtk.signal_connect(guidata[c,:save_as], "clicked") do widget
    #     save_as(f)
    # end
    # Gtk.signal_connect(guidata[c,:zoom_button], "clicked") do widget
    #     panzoom_cb(f)
    # end
    # Gtk.signal_connect(guidata[c,:fullview], "clicked") do widget
    #     fullview_cb(f)
    # end
    Gtk.signal_connect(save_as_wrapper, guidata[c,:save_as], "clicked", Nothing, (), false, f)
    Gtk.signal_connect(panzoom_wrapper, guidata[c,:zoom_button], "clicked", Nothing, (), false, f)
    Gtk.signal_connect(fullview_wrapper, guidata[c,:fullview], "clicked", Nothing, (), false, f)
    Immerse.lasso_initialize(f)
end

save_as_wrapper(::Ptr, f) = (save_as(f); nothing)
panzoom_wrapper(::Ptr, f) = (panzoom_cb(f); nothing)
fullview_wrapper(::Ptr, f) = (fullview_cb(f); nothing)

"""
`figure(;name="Figure \$n", width=400, height=400)` creates a new
figure window for displaying plots.

`figure(n)` raises the `n`th figure window and makes it the current
default plotting window, and returns the
"""
function figure(;name::String="Figure $(nextfig(_display))",
                 width::Integer=400,    # TODO: make configurable
                 height::Integer=400)
    i = nextfig(_display)
    f = Figure()

    Gtk.signal_connect(on_figure_destroy, f, "destroy", Nothing, (), false, (i, _display))

    gtkwindow(name, width, height, f)

    initialize_toolbar_callbacks(f)

    addfig(_display, i, f)
    display(f.canvas, f)
    i
end

function on_figure_destroy(widgetptr::Ptr, user_data)
    i, _display = user_data
    dropfig(_display, i)
    nothing
end

function figure(i::Integer; displayfig::Bool = true)
    switchfig(_display, i)
    fig = curfig(_display)
    displayfig && display(_display, fig)
    Gtk.present(Gtk.toplevel(fig.canvas))
    fig
end

# create and initialize a Figure to use an existing GtkCanvas
function figure(c::GtkCanvas)
    i = nextfig(_display)
    f = Figure(c)
    initialize_toolbar_callbacks(f)
    addfig(_display, i, f)
    i
end

"`gcf()` (\"get current figure\") returns the current figure number"
gcf() = _display.current_fig

"`scf()` (\"show current figure\") raises (makes visible) the current figure"
scf() = figure(gcf())

"""
`closefig(n)` closes the `n`th figure window.

`closefig()` closes the current figure window.
"""
closefig() = closefig(_display.current_fig)

closefig(i::Integer) = (fig = getfig(_display,i); clear_hit(fig); gtkdestroy(getfig(_display,i).canvas))

"`closeall()` closes all existing figure windows."
closeall() = (map(closefig, collect(keys(_display.figs))); nothing)

function createPlotGuiComponents()
    box = Gtk.GtkBox(:v)
    tb = Gtk.GtkToolbar()
    push!(box, tb)

    save_as = Gtk.GtkToolButton("gtk-save-as")     # document-save-as
    zb = Gtk.GtkToggleToolButton("gtk-find")       # edit-find
    fullview = Gtk.GtkToolButton("gtk-zoom-100")   # zoom-original
    lasso_button = Gtk.GtkToggleToolButton()
    Gtk.GAccessor.icon_widget(lasso_button, Gtk.GtkImage(joinpath(Immerse.HOME, "images", "lasso_icon_16.png")))
    push!(tb, save_as)
    push!(tb, Gtk.GtkSeparatorToolItem())
    push!(tb, zb)
    push!(tb, fullview)
    push!(tb, Gtk.GtkSeparatorToolItem())
    push!(tb, lasso_button)

    Gtk.G_.icon_size(tb,1)#small icon size (16px)

    c = Gtk.GtkCanvas()
    Gtk.set_gtk_property!(c, :expand, true)
    push!(box, c)

    guidata[c, :save_as] = save_as
    guidata[c, :zoom_button] = zb
    guidata[c, :fullview] = fullview
    guidata[c, :lasso_button] = lasso_button
    box, tb, c
end

function gtkwindow(name, w, h, fig)
    win = Gtk.GtkWindow(fig, name, w, h)
    guidata[win, :toolbar] = fig.tb
    showall(win)
end

function clear_guidata(c)
    gd = guidata[c]
    to_delete = Symbol[]
    for (k,v) in gd
        if !(k in [:save_as, :zoom_button, :fullview, :lasso_button])
            push!(to_delete, k)
        end
    end
    for k in to_delete
        delete!(gd, k)
    end
end

function clear_guistate!(f::Figure)
    if initialized(f.panzoom_cb)
        disconnect(f, f.panzoom_cb)
    end
    Gtk.set_gtk_property!(guidata[f.canvas, :zoom_button], :active, false)
    f.panzoom_cb = PanZoomCallbacks()
end

function gtkdestroy(c::GtkCanvas)
    Gtk.destroy(Gtk.toplevel(c))
    nothing
end

function render_finish(prep; kwargs...)
    p = _plot(prep)
    root_context = Gadfly.render_prepared(prep...; kwargs...)

    ctx =  Compose.pad_inner(root_context, p.theme.plot_padding...)

    if p.theme.background_color != nothing
        Compose.compose!(ctx, (Compose.context(order=-1000000),
                               Compose.fill(p.theme.background_color),
                               Compose.stroke(nothing), Compose.rectangle()))
    end

    return ctx
end

function set_ticks!(coord::Cartesian, aes::Aesthetics, xview, yview)
    xtick = Gadfly.optimize_ticks(xview.min, xview.max)[1]
    ytick = Gadfly.optimize_ticks(yview.min, yview.max)[1]
    aes.xtick = aes.xgrid = xtick
    aes.ytick = aes.ygrid = ytick
    aes.xtickvisible = fill(true, length(xtick))
    aes.ytickvisible = fill(true, length(ytick))
    aes.xtickscale = ones(length(xtick))
    aes.ytickscale = ones(length(ytick))
    coord = Cartesian(xmin=xview.min,
                      xmax=xview.max,
                      ymin=yview.min,
                      ymax=yview.max,
                      xflip=coord.xflip,
                      yflip=coord.yflip,
                      fixed=coord.fixed,
                      aspect_ratio=coord.aspect_ratio,
                      raster=coord.raster)
    coord, aes
end

function set_limits!(f::Figure, xview, yview)
    coord = _coord(f)
    aes = _aes(f)
    if (coord.xmin, coord.xmax) != (xview.min, xview.max) ||
       (coord.ymin, coord.ymax) != (yview.min, yview.max)
        coord, aes = set_ticks!(coord, aes, xview, yview)
        setcoord!(f, coord)
        cc = render_finish(f.prepped; dynamic=false)
        f.cc = apply_tweaks!(cc, f.tweaks)
    end
    f
end

function dev2data(widget, x, y)
    xmm, ymm = x/ppmm, y/ppmm
    pc = guidata[widget,:panelcoords]
    @assert length(pc) == 1
    parent_box, unit_box, transform = pc[1]
    absolute_to_data(xmm*mm, ymm*mm, transform, unit_box, parent_box)
end

GtkUtilities.panzoom(f::Figure, args...) = panzoom(f.canvas, args...)
function GtkUtilities.panzoom(f::Figure)
    aes = _aes(f)
    xview = (aes.xviewmin, aes.xviewmax)
    yview = (aes.yviewmin, aes.yviewmax)

    if xview == (nothing, nothing)
        xview = (minimum(aes.x),maximum(aes.x))
    end
    if yview == (nothing, nothing)
        yview = (minimum(aes.y),maximum(aes.y))
    end

    panzoom(f.canvas, xview, yview)
end

function GtkUtilities.panzoom_mouse(f::Figure; kwargs...)
    aes = _aes(f)
    xtick = aes.xtick != nothing ? aes.xtick : aes.x
    ytick = aes.ytick != nothing ? aes.ytick : aes.y
    xflip = xtick[end] < xtick[1]
    yflip = ytick[end] > ytick[1]
    panzoom_mouse(f.canvas; xpanflip=xflip, ypanflip=yflip, user_to_data=(c,x,y)->dev2data(c,x,y), kwargs...)
end

function GtkUtilities.panzoom_key(f::Figure; kwargs...)
    aes = _aes(f)
    xtick = aes.xtick != nothing ? aes.xtick : aes.x
    ytick = aes.ytick != nothing ? aes.ytick : aes.y
    xflip = xtick[end] < xtick[1]
    yflip = ytick[end] > ytick[1]
    panzoom_key(f.canvas; xpanflip=xflip, ypanflip=yflip, kwargs...)
end

function block(f::Figure, pcz::PanZoomCallbacks)
    c = f.canvas
    Gtk.signal_handler_block(c, pcz.idpzk)
    pop!((c.mouse,:scroll))
    pop!((c.mouse,:button1press))
end

function disconnect(f::Figure, pcz::PanZoomCallbacks)
    c = f.canvas
    Gtk.signal_handler_disconnect(c, pcz.idpzk)
    pop!((c.mouse,:scroll))
    pop!((c.mouse,:button1press))
end

function unblock(f::Figure, pcz::PanZoomCallbacks)
    Gtk.signal_handler_unblock(f.canvas, pcz.idpzk)
    panzoom_mouse(f)
end

function PanZoomCallbacks(f::Figure)
    panzoom(f)
    panzoom_mouse(f)
    PanZoomCallbacks(panzoom_key(f))
end

const file_backends = Dict(".svg"=>Compose.SVG, ".png"=>Compose.PNG, ".pdf"=>Compose.PDF, ".ps"=>Compose.PS)

function save_as(f::Figure)
    isempty(f) && return nothing
    extensions = (".svg", ".png", ".pdf", ".ps")
    selection = Gtk.save_dialog("Save figure as file", Gtk.toplevel(f.canvas), map(x->string("*",x), extensions))
    isempty(selection) && return nothing
    basename, ext = splitext(selection)
    if !in(ext, extensions)
        Gtk.error_dialog("Extension $ext not recognized: use .svg, .png, .pdf, or .ps")
        return nothing
    end
    w, h = width(f.canvas), height(f.canvas)
    Compose.draw(file_backends[ext](selection, w*Compose.px, h*Compose.px), f.cc)
    nothing
end

function panzoom_cb(f::Figure)
    if isempty(f)
        Gtk.set_gtk_property!(guidata[f.canvas, :zoom_button], :active, false)
        return nothing
    end
    if !initialized(f.panzoom_cb)
        f.panzoom_cb = PanZoomCallbacks(f)
    else
        state = Gtk.get_gtk_property(guidata[f.canvas, :zoom_button], :active, Bool)
        if state
            unblock(f, f.panzoom_cb)
        else
            block(f, f.panzoom_cb)
        end
    end
end

function fullview_cb(f::Figure)
    if initialized(f.panzoom_cb)
        GtkUtilities.zoom_reset(f.canvas)
    end
end

const HOME = splitdir(splitdir(@__FILE__)[1])[1]

# end # module


text:
export HubertAmplifier, HubertAmplifierParams

Base.@kwdef struct HubertAmplifierParams <: DeviceParams
	port::String = "/dev/ttyUSB1"
	baudrate::Integer = 9600
	mode::AmplifierMode = AMP_VOLTAGE_MODE # This should be the safe default
	voltageMode::AmplifierVoltageMode = AMP_LOW_VOLTAGE_MODE # This should be the safe default
	matchingNetwork::Integer = 1
end
HubertAmplifierParams(dict::Dict) = params_from_dict(HubertAmplifierParams, dict)

Base.@kwdef mutable struct HubertAmplifier <: Amplifier
  "Unique device ID for this device as defined in the configuration."
  deviceID::String
  "Parameter struct for this devices read from the configuration."
  params::HubertAmplifierParams
	"Flag if the device is optional."
	optional::Bool = false
	"Flag if the device is present."
  present::Bool = false
  "Vector of dependencies for this device."
  dependencies::Dict{String, Union{Device, Missing}}

  driver::Union{SerialPort, Missing} = missing
end

function init(amp::HubertAmplifier)
  @debug "Initializing Hubert amplifier unit with ID `$(amp.deviceID)`."

	@warn "The code for the Hubert amplifier has not yet been tested within MPIMeasurements!"

	amp.driver = SerialPort(port)
	open(amp.driver)
	set_speed(amp.driver, amp.params.baudrate)

	_hubertSetStartupParameters(amp.driver)

	# Set values given by configuration
	mode(amp, amp.params.mode)
	voltageMode(amp, amp.params.voltageMode)
	matchingNetwork(amp, amp.params.matchingNetwork)

	amp.present = true
end

checkDependencies(amp::HubertAmplifier) = true

Base.close(amp::HubertAmplifier) = close(amp.driver)

# main communication function
function _hubertSerial(sp::SerialPort, input::Array{UInt8}, ack::Array{UInt8})

	write(sp, input)
	sleep(0.2)	#Huberts needs at least 100ms
	answer_utf = readuntil(sp, '\n', 0.1)
	answer_hex = encode(answer_utf, "UTF-8")	#using StringEncodings
	flush(sp)

	if answer_hex == ack
		@debug "Hubert acknowleged: '$(answer_hex)' set."
	else
		@error "Hubert: Serial communication error. Expected '$(ack)', received '$(answer_hex)'"
	end

	return nothing
end

function _hubertSerial(sp::SerialPort, input::Array{UInt8}) #for querys
	write(sp, input)
	sleep(0.2)	#Huberts needs at least 100ms
	answer_utf = readuntil(sp, '\n', 0.1)
	answer_hex = encode(answer_utf, "UTF-8")	#using StringEncodings
	flush(sp)

	return answer_hex
end

function _hubertSetStartupParameters(sp::SerialPort)
	# EINschaltparameter - werde nur beim Einschalten aktualisiert
	input = UInt8[0x0B, 0x2E, #setzten der erw. Einstellung auf einmal
								0x00,		#1. Strommessbereich: high 00, low 01 - use high!
								0x07,		#2. RC_network: depends ... (05 or 07). 07: trimmed
								0x01,		#3. Mode: voltage 00, current 01
								0x00,		#4. 00 (empty)
								0x0f,		#5. Highbyte Limit Control (0x00 - 0x0F)
								0xc3,		#6. Lowbyte Limit Control (0x00 - 0xFF)
								0x01,		#7. Interlock: latching 00, live 01 - needs to be live for RP interlock
								0x00,		#8. Limit Control: current 00, voltage 01
								0x09]		#9. Betriebsspann. mid: 05, high: 09 - do not use auto (01)!
	ack = UInt8[0x2e]
	_hubertSerial(sp, input, ack)
	sleep(0.1)
	input = UInt8[0x03, 0x5D, 0x00] #sensing needs to be off (otherwise DC offset) - but Huberts jumps straight to Overvoltage after Interlock.
	ack = UInt8[0x5D]
	_hubertSerial(sp, input, ack)
	return nothing
end

state(amp::HubertAmplifier) = @error "Getting the current activation state of the amplifier is not yet supported."

function turnOn(amp::HubertAmplifier)
	input = UInt8[0x03, 0x35, 0x01]
	ack   = UInt8[0x01]
	_hubertSerial(amp.driver, input, ack)
	return nothing
end

function turnOff(amp::HubertAmplifier)
	input = UInt8[0x03, 0x35, 0x00]
	ack   = UInt8[0x00]
	_hubertSerial(amp.driver, input, ack)
	return nothing
end

mode(amp::HubertAmplifier) = @error "Getting the current amplifier mode is not yet supported."

function mode(amp::HubertAmplifier, mode::AmplifierMode)
	#Mode: voltage 00, current 01
	if mode == AMP_CURRENT_MODE
		input = UInt8[0x03, 0x2A, 0x01]
	elseif mode == AMP_VOLTAGE_MODE
		input = UInt8[0x03, 0x2A, 0x00]
	else
		throw(ScannerConfigurationError("Unsupported mode for Hubert amplifier: $mode."))
	end

	ack = UInt8[0x2A]
	_hubertSerial(amp.driver, input, ack)

	return nothing
end

voltageMode(amp::HubertAmplifier) = @error "Getting the current amplifier voltage mode is not yet supported."

function voltageMode(amp::HubertAmplifier, mode::AmplifierVoltageMode)
	#high is required, may set to low for below 15mt [??? toDo: after cal]
	#9. Betriebsspann. mid: 05, high: 09 - do not use auto (01)!
	if mode == AMP_HIGH_VOLTAGE_MODE
		input = UInt8[0x03, 0x54, 0x09]
	elseif mode == AMP_LOW_VOLTAGE_MODE
		input = UInt8[0x03, 0x54, 0x05]
	else
		throw(ScannerConfigurationError("Unsupported voltage mode for Hubert amplifier: $mode."))
	end

	ack = UInt8[0x54]
	_hubertSerial(amp.driver, input, ack)

	return nothing
end

matchingNetwork(amp::HubertAmplifier) = @error "Getting the current amplifier network is not yet supported."

function matchingNetwork(amp::HubertAmplifier, network::Integer)
	#2. RC_network: (05 or 07). 07: trimmed //specific for THIS hubert amplifier
	if network >= 1 && network <= 7
		input = UInt8[0x03, 0x29, UInt8(network)]
		ack = UInt8[network]
		_hubertSerial(amp.driver, input, ack)
	else
		throw(ScannerConfigurationError("Hubert: '$(network)' is not a RC-Network {1...7}."))
	end

	return nothing
end

function temperature(amp::HubertAmplifier)
	input = UInt8[0x02, 0x04]
	hex = _hubertSerial(amp.driver, input)
	return Int(hex[1])
end


text:
using Pkg.Artifacts
using SHA

using HTTP

const ARTIFACTS_TOML = joinpath(@__DIR__, "Artifacts.toml")

const CLOUDFRONT_DISTRIB_PREFIX = "https://cdn.huggingface.co"
const S3_BUCKET_PREFIX = "https://s3.amazonaws.com/models.huggingface.co/bert"
const DEFAULT_CONFIG_NAME = "config.json"
const DEFAULT_WEIGHT_NAME = "pytorch_model.bin"

"""
Workflow:
1. Pretrained model is from huggingface (yes=>a. / no=>b.)

a. Already use the model in python before. (yes=>a-1. / no=>a-2.)

a-1. There is a cached pretrain file on the computer. Use [`get_or_download_hgf_file`](@ref) directly. This will copy
the file from cached dir to our Artifacts dir and register on Artifacts.toml.
a-2. No cached files. Also use [`get_or_download_hgf_file`](@ref) directly. This will download the file from
huggingface's model server to our Artifacts dir and register on Artifact.toml.

b. Using the local pretrained files. Use [`find_or_register_hgf_file_hash`](@ref) to register the file to our
Artifacts system. Once the file is registered, you can find the entry appear on Artifacts.toml.

2. Once the model is managed under Julia's Artifacts system. we can use either [`get_registered_file_dir`](@ref) or
[`get_or_download_hgf_file`](@ref) to get the pretrain file or dir.
"""

"""
  `get_registered_config_path(model_name; config=DEFAULT_CONFIG_NAME)`

Get the config path of a given `model_name` which should be already registered on `Artifacts.toml`.
"""
get_registered_config_path(model_name; config=DEFAULT_CONFIG_NAME) = joinpath(get_registered_file_dir(joinpath(model_name, config)), config)

"""
  `get_registered_weight_path(model_name; weight=DEFAULT_WEIGHT_NAME)`

Get the weight path of a given `model_name` which should be already registered on `Artifacts.toml`.
"""
get_registered_weight_path(model_name; weight=DEFAULT_WEIGHT_NAME) = joinpath(get_registered_file_dir(joinpath(model_name, weight)), weight)

"""
  `get_registered_file_dir(name)`

Get the dir path of a given `name` which should be already registered on `Artifacts.toml`.

See also: [`get_registered_config_path`](@ref), [`get_registered_weight_path`](@ref)
"""
function get_registered_file_dir(name)
  global ARTIFACTS_TOML
  hash = artifact_hash(name, ARTIFACTS_TOML)
  isnothing(hash) && error("$name not registered.")
  artifact_path(hash)
end

"""
  `get_or_download_hgf_config(model_name; config=DEFAULT_CONFIG_NAME)`

get the config file path of the given `model_name` and `file_name`

See also: [`get_or_download_hgf_file`](@ref)
"""
get_or_download_hgf_config(model_name; config=DEFAULT_CONFIG_NAME) = get_or_download_hgf_file(model_name, config)

"""
  `get_or_download_hgf_weight(model_name; weight=DEFAULT_WEIGHT_NAME)`

get the weight file path of the given `model_name` and `file_name`

See also: [`get_or_download_hgf_file`](@ref)
"""
get_or_download_hgf_weight(model_name; weight=DEFAULT_WEIGHT_NAME) = get_or_download_hgf_file(model_name, weight)

"""
  `get_or_download_hgf_file(model_name, file_name)`

get the file path of the given `model_name` and `file_name`. Automatically download and register from
huggingface server if file not found on `Artifacts.toml`. To use with a local file, register with
[`find_or_register_hgf_file_hash`](@ref) first.
"""
function get_or_download_hgf_file(model_name, file_name)
  global CLOUDFRONT_DISTRIB_PREFIX
  hash = find_or_register_hgf_file_hash(CLOUDFRONT_DISTRIB_PREFIX, model_name, file_name)
  return joinpath(artifact_path(hash), file_name)
end

isregistered(name) = (global ARTIFACTS_TOML; !isnothing(artifact_hash(name, ARTIFACTS_TOML)))

isurl(path::AbstractString) = startswith(path, "http://") || startswith(path, "https://")

"""
  `find_or_register_hgf_config_hash(path, model_name; config=DEFAULT_CONFIG_NAME)`

Get the artifacts hash of config of the give model.

See also: [`find_or_register_hgf_file_hash`](@ref)
"""
find_or_register_hgf_config_hash(path, model_name; config=DEFAULT_CONFIG_NAME) = find_or_register_hgf_file_hash(path, model_name, config)

"""
  `find_or_register_hgf_weight_hash(path, model_name; weight=DEFAULT_WEIGHT_NAME)`

Get the artifacts hash of weight of the give model.

See also: [`find_or_register_hgf_file_hash`](@ref)
"""
find_or_register_hgf_weight_hash(path, model_name; weight=DEFAULT_WEIGHT_NAME) = find_or_register_hgf_file_hash(path, model_name, weight)

"""
  `find_or_register_hgf_file_hash(path, model_name, file_name)`

Get the artifacts hash of the give `<model_name>/<file_name>`. If not found in Artifacts.toml, get the file from `path`
and register on Artifacts.toml. `path` can be either a url or a local path.
"""
function find_or_register_hgf_file_hash(path, model_name, file_name)
  global ARTIFACTS_TOML
  entry_name = joinpath(model_name, file_name)
  file_hash = artifact_hash(entry_name, ARTIFACTS_TOML)

  if isnothing(file_hash) || !artifact_exists(file_hash) # not found, try get model and register
    if isurl(path) # path is a url, download with huggingface url format
      islegacy = !occursin('/', model_name)
      url = islegacy ?
        join((path, model_name, file_name), '/', '-') :
        joinpath(path, entry_name)
      filegetter = Base.Fix1(hgf_file_download, url)
    else # path is a local file, cp to artifact dir
      file_path = joinpath(path, file_name)
      isfile(file_path) || error("Can't register to $entry_name: file $file_path not found.")
      filegetter = function (dest)
        cp(file_path, dest)
      end
    end

    file_hash = create_artifact() do artifact_dir
        filegetter(joinpath(artifact_dir, file_name))
    end

    # register if no in artifact.toml
    if isnothing(artifact_hash(entry_name, ARTIFACTS_TOML))
      bind_artifact!(ARTIFACTS_TOML, entry_name, file_hash)
    end
  end

  file_hash
end

"Get the cache dir used by huggingface/transformers"
function get_th_hgf_cache_dir()
  return get(ENV, "TRANSFORMERS_CACHE",
             get(ENV, "PYTORCH_TRANSFORMERS_CACHE",
                 get(ENV, "PYTORCH_PRETRAINED_BERT_CACHE",
                     joinpath(
                       get(ENV, "TORCH_HOME", get(ENV, "XDG_CACHE_HOME", "~/.cache")),
                       "torch", "transformers")))) |> expanduser
end

"Find the huggingface cached files from their cache dir"
function get_hgf_cached_file(url; cache_dir=get_th_hgf_cache_dir())
  isdir(cache_dir) || return nothing
  global CLOUDFRONT_DISTRIB_PREFIX, S3_BUCKET_PREFIX
  iscdn = startswith(url, CLOUDFRONT_DISTRIB_PREFIX)
  iss3 = startswith(url, S3_BUCKET_PREFIX)

  hashstr(x) = bytes2hex(sha256(x))
  match(hash) = filter(f->startswith(f, hash) && !endswith(f, ".json") && !endswith(f, ".lock"),
                       readdir(cache_dir))

  hash = hashstr(url)
  mf = match(hash)
  !isempty(mf) && return joinpath(cache_dir, last(mf)) # if found then return

  if iscdn # using cdn, retry with s3
    alter_url = replace(url, CLOUDFRONT_DISTRIB_PREFIX=>S3_BUCKET_PREFIX)
  elseif iss3 # using s3, retry with cdn
    alter_url = replace(url, S3_BUCKET_PREFIX=>CLOUDFRONT_DISTRIB_PREFIX)
  else # custom url
    return nothing
  end

  # retry with alternative url
  hash = hashstr(alter_url)
  mf = match(hash)

  if !isempty(mf)
    return joinpath(cache_dir, last(mf))
  else
    return nothing
  end
end

"Get files from huggingface. Copy if found in cache dir, otherwise download from their server."
function hgf_file_download(url, dest)
  cached_file = get_hgf_cached_file(url)
  file_name = basename(dest)
  if isnothing(cached_file) # not found, just download
    @info "No local $file_name found. downloading..."
    HTTP.download(url, dest)
  else # found, copy from cache
    @info "PyTorch cached $file_name found. copying..."
    cp(cached_file, dest)
  end
end


text:
for eps in [0.1, 0.05, 0.02, 0.01]
    n = ceil(1/eps^2)
    println("For epsilon = $eps set n = $n")
end

max_stars_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_issues_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_forks_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook


text:
"""
Generate a regular lattice of points
$(TYPEDSIGNATURES)

# Arguments
- `dx` → element spacing on the x axis
- `dy` → element spacing on the y axis
- `ds` → displacement along x between rows of elements
- `f_cond::Function` → function of two arguments (`f(x,y) = ...`) returning `true` if element at position `x,y` must be kept and `false` otherwise

# Keyord Arguments
- `dx0 = 0` → x coordinate of the origin of the lattice
- `dy0 = 0` → y coordinate of the origin of the lattice
- `M::Int = 70` → Number of elements to generate per row of points before appliying the filtering function `f_cond`
- `N::Int = M` → Number of rows of points to generate before appliying the filtering function `f_cond`
"""
function generate_regular_lattice(dx::T, dy::T, ds::T, f_cond::Function = (x, y)->true;dx0 = T(0), dy0 = T(0), M::Int = 70,N::Int = M) where T<:Real
	# Function to generate x position as function of row,column number m,n
	x(m, n) = m * dx + n * ds + dx0
	# Function to generate y position as function of row,column number m,n
	y(n) = n * dy + dy0
	# Generate the elements. For each row, shift the columns to always have the search domain around x=0
	# out = [(x(m, n), y(n)) for n in -N:N for m in range(-M - round(n * ds / dx);length = 2 * M + 1) if f_cond(x(m, n), y(n))]
	out = [(x(m - round(Int,n * ds / dx), n), y(n)) for n in -N:N,m in -M:M if f_cond(x(m - round(Int,n * ds / dx), n), y(n))]
	return out
end
generate_regular_lattice(dx::Real,dy::Real,ds::Real,args...;kwargs...) = generate_regular_lattice(promote(dx,dy,ds)...,args...;kwargs...)

"""
    generate_rect_lattice(spacing_x::Real,spacing_y::Real[,f_cond];kwargs...)
# Summary
Generate a rectangular lattice of points (with different spacing among x and y directions)
# Arguments
- `spacing_x` → spacing between points on the x axis
- `spacing_y` → spacing between points on the y axis

See [`generate_regular_lattice`](@ref) for a description of `f_cond` and of  the keyword arguments
"""
generate_rect_lattice(spacing_x::Real,spacing_y::Real,args...;kwargs...) = generate_regular_lattice(spacing_x,spacing_y,0,args...;kwargs...)

"""
`generate_square_lattice(spacing::Real[,f_cond];kwargs...)`
# Summary
Generate a square lattice of points (with equal spacing among x and y directions)
# Arguments
- `spacing` → spacing between points on both x and y axis

See [`generate_regular_lattice`](@ref) for a description of `f_cond` and of  the keyword arguments
"""
generate_square_lattice(spacing::Real,args...;kwargs...) = generate_regular_lattice(spacing,spacing,0,args...;kwargs...)

"""
`generate_hex_lattice(spacing::Real[,f_cond];kwargs...)`
# Summary
Generate a hexagonal lattice of points (with equal distance between neighboring points).
The hexagonal lattice generated by this function has distance between points on the same
column √3 times greater than the distance between points on the same row.
# Arguments
- `spacing` → spacing between points

See [`generate_regular_lattice`](@ref) for a description of `f_cond` and of  the keyword arguments
"""
generate_hex_lattice(spacing::Real,args...;kwargs...) = generate_regular_lattice(spacing .* (1,√3/2,.5)...,args...;kwargs...)

# Get the conversion from linear to db and viceversa
"""
$(TYPEDSIGNATURES)
Convert a number from linear to dB
"""
lin2db(x::Real) = 10log10(x)
"""
$(TYPEDSIGNATURES)
Convert a number from dB to linear
"""
db2lin(x::Real) = 10^(x/10)

# Convert between frequency and wavelength
"""
$(TYPEDSIGNATURES)
Get the wavelength (in m) starting from the frequency (in Hz)
"""
f2λ(f::Real) = c₀/f
"""
$(TYPEDSIGNATURES)
Get the frequency (in Hz) starting from the wavelength (in m)
"""
λ2f(λ::Real) = c₀/λ

## Here we have the functions for the coloring computation
function compute_F_cell(max_colours::Int;grid_max::Int=25)
    #=
    This function is used to compute all the possible 2x2 lattice generating matrices for possible coloring schemes up to 'max_colours' colors
    Computation is done with a brute-force approach, generating all possible 2x1 vectors with maximum elements up to grid_max
    =#

    # Find the current length of the pre-computed vector
    current_length = length(F_reuse_matrix.square)
    if current_length >= max_colours
        # We already computed the function for the required number of colors
        return
    end
    n_missing = max_colours - current_length
    append!(F_reuse_matrix.square,Vector{SMatrix{2,2,Float64,4}}(undef,n_missing))
    append!(F_reuse_matrix.triangular,Vector{SMatrix{2,2,Float64,4}}(undef,n_missing))
    # Compute the matrix for the square lattice
    _coloring_inner_bruteforce!(F_reuse_matrix.square,I,grid_max)
    # Compute the matrix for the triangular lattice
    _coloring_inner_bruteforce!(F_reuse_matrix.triangular,@SMatrix([1 0;cosd(60) sind(60)]),grid_max)
end

function _coloring_inner_bruteforce!(T_mat, rot_mat, grid_max)
    max_colours = length(T_mat)
    check_vec = fill((typemax(Int),typemax(Int)),max_colours)
    @inline norm2(x) = sum(abs2.(x))
    @inbounds for x1 = 0:grid_max, y1 = 0:grid_max, x2 = -grid_max:grid_max, y2 = 0:grid_max
        mat = @SMatrix [x1 y1;x2 y2]
        # Compute the determinant
        t_det = Int(det(mat))
        # Skip points which are not likely to give useful results
        if (t_det < 1)  || (t_det > max_colours) || (t_det > grid_max^2) || (maximum(abs.(mat)) > ceil(sqrt(t_det) + 3))
            continue
        end
        # Compute the angle between the basis vectors identified by [x1,y1] and [x2,y2]
        angle = abs(atan(y1,x1) - atan(y2,x2))
        # Skip cases where the angle between vectors is either too acute or too obtuse
        if abs(π/2-angle) > π/4
            continue
        end
        # Create temp variables for computation of the minimum distance
        dmat = mat*rot_mat
        # Compute frobenius norm and minimum squared distance for the candidate lattice generating matrix
        frobe = round(Int,norm2(dmat))
        # display(frobe)
        # Minimum squared distance is either the modulo of one of the vectors or the modulo of sum or difference of them
        dmin = round(Int,minimum((norm2(dmat[1,:]), norm2(dmat[2,:]), norm2(sum(dmat;dims=1)), norm2(diff(dmat;dims=1)))))
        # Check if the current realization is better than the saved one
        if isless((-dmin,frobe),check_vec[t_det])
            # Update the check_vec
            check_vec[t_det] = (-dmin,frobe)
            # Update the vector containing the generating matrices
            T_mat[t_det] = mat'
        end
    end
end

function generate_F_reuse_matrix(lattice_type::Symbol=:triangular,N_Colours::Int=4;max_colours::Int=max(10,N_Colours))
    compute_F_cell(max_colours)
    return getproperty(F_reuse_matrix,lattice_type)[N_Colours]
end

"""
    generate_colors(BeamCenters::AbstractVector,N_Colours::Int=4;lattice_type::Symbol=:triangular)

Provide a the coloring breakdown for a given set of lattice points.
# Arguments
- `BeamCenters` → Vector of Tuple or StaticVectors expressing the U-V coordinates of each point in the lattice for which the coloring is being computed
- `N_Colours` → Number of colors to divide the lattice in. Defaults to `4`

# keyword Arguments
- `lattice_type` → Symbol that can either be `:triangular` or `:square`, idenifying the type of lattice. Defaults to `:triangular`
"""
function generate_colors(BeamCenters::AbstractVector,N_Colours::Int=4;first_color_coord=nothing,first_color_idx=nothing,precision_digits::Int=7,lattice_type::Symbol=:triangular)
    #=
    **************************************************************************
       Generate frequency colouring file
    **************************************************************************

     References:
     [1]   "On the frequency allocation for mobile radio telephone systems", C. de
           Almeida; R. Palazzo, Proceedings of 6th International Symposium on
           Personal, Indoor and Mobile Radio Communications, Year: 1995, Volume: 1, Pages: 96 - 99 vol.1
     [2]   P. Angeletti, "Simple implementation of vectorial modulo operation based
           on fundamental parallelepiped," in Electronics Letters, vol. 48, no. 3, pp. 159-160,
           February 2 2012. doi: 10.1049/el.2011.3667

     Authors: Alberto Mengali, 2020, European Space Agency

    Input Arguments:
    BeamCenters:         A vector containing the U,V coordinates of the beam centers as tuples or staticvectors
    N_Colours:           An integer number specifying the number of colors to generate in the association
    first_color_coord:   A tuple or static vector containing the U,V coordinates of the beam that will have the first color
    first_color_idx:     The beam index of the beam containing the first color, either this variable of first_order_coord can be specified, not together
    precision_digits:    The number of digits to be used in the flooring operation
    =#

    # Check if either the first color coordinates or first color idx are given
    if first_color_coord === nothing
        if first_color_idx === nothing
            # Initialize to coorinates of the first beam
            first_color_idx = 1;
        end
        first_color_coord = BeamCenters[first_color_idx]
    else
        if first_color_idx !== nothing
            @warn "Both first_color_idx and first_color_coord were given, disregarding the idx variable"
        end
    end
    if lattice_type ∉ (:triangular, :square)
        @error "The specified lattice type ($lattice_type) is not recognized, it should be either :triangular or :square"
    end

    n_beams = length(BeamCenters)

    # If only 1 color is requested, simply return the trivial result
    if N_Colours == 1
        Colours = ones(n_beams)
        nbeams_per_color = n_beams
        idxs = fill(fill(true,nbeams))
    end

    # Find the minimum distance in U and V
    minU_dist = minimum(diff(sort(first.(BeamCenters)) |> x -> unique!(y -> round(y;digits=precision_digits),x)))
    minV_dist = minimum(diff(sort(last.(BeamCenters)) |> x -> unique!(y -> round(y;digits=precision_digits),x)))
    if lattice_type === :triangular
        beamU_spacing = 2minU_dist
        beamV_spacing = 2minV_dist
        # Matrix to normalize the grid points in u-v into a integer grid with y(v) axis not pointing north but north-east with 60° inclination
        D = @SMatrix [1 -1/2;0 1]
    elseif lattice_type === :square
        beamU_spacing = minU_dist
        beamV_spacing = minV_dist
        D = @SMatrix [1 0;0 1]
    end
    # Get the coloring generating matrix
    F_reuse_matrix = generate_F_reuse_matrix(lattice_type,N_Colours)
    # Create the set that will contain the unique results
    unique_colors_vecs = SVector{2,Int}[]
    # Initialize the colors vector
    Colors = similar(BeamCenters, Int)
    @inbounds for (n,p₀) ∈ enumerate(BeamCenters)
        # Compute the integer beam centers
        p₀_normalized = p₀ ./ SA_F64[beamU_spacing, minV_dist]
        # Transform the beam centers from u-v coordinates in radians into integer indexes
        beam_index_vector = round.(Int,D*p₀_normalized);
        # Find the values of the beam indexes modulo F_reuse_matrix
        unique_beam_index = round.(Int,beam_index_vector .- (F_reuse_matrix*floor.(round.(inv(F_reuse_matrix)*beam_index_vector,digits=precision_digits))))
        # Check if this color has already been assigned/registered
        idx = findfirst(x -> x == unique_beam_index,unique_colors_vecs)
        if idx isa Nothing
            push!(unique_colors_vecs, unique_beam_index)
            cidx = length(unique_colors_vecs)
        else
            cidx = idx
        end
        Colors[n] = cidx
    end

    return Colors
end


text:
"""
An aggregator interface for Delay SSA-like algorithms.
### Required Fields
- `next_jump`          # the next jump to execute
- `prev_jump`          # the previous jump that was executed
- `next_jump_time`     # the time of the next jump
- `end_time`           # the time to stop a simulation
- `cur_rates`          # vector of current propensity values
- `sum_rate`           # sum of current propensity values
- `ma_jumps`           # any MassActionJumps for the system (scalar form)
- `rates`              # vector of rate functions for ConstantRateJumps
- `affects!`           # vector of affect functions for ConstantRateJumps
- `save_positions`     # tuple for whether to save the jumps before and/or after event
- `rng`                # random number generator
- `next_delay`         # the index of the delay channel
- `num_next_delay`     # how many times needed for updating the states in the next_delay channel
- `time_to_next_jump`  # the time to the next jump (time gap)
- `dt_delay`           # the time to the next delay reaction
### Optional fields:
- `dep_gr`             # dependency graph, dep_gr[i] = indices of reactions that should
                       # be updated when rx i occurs.
"""
abstract type AbstractDSSAJumpAggregator <: AbstractJumpAggregator end

DiscreteCallback(c::AbstractDSSAJumpAggregator) = DiscreteCallback(c, c, initialize = c, save_positions = c.save_positions)

########### The following routines are templates for all Delay SSAs ###########
########### Generally they should not need to be overloaded.  ###########

## Users will normally define:
# aggregate
# initialize!
# execute_jumps!
# generate_jumps!

# condition for jump to occur
@inline function (p::AbstractDSSAJumpAggregator)(u, t, integrator)
    p.next_jump_time == t
end

# executing jump at the next jump time
function (p::AbstractDSSAJumpAggregator)(integrator)
    execute_jumps!(p, integrator, integrator.u, integrator.p, integrator.t)
    generate_jumps!(p, integrator, integrator.u, integrator.p, integrator.t)
    register_next_jump_time!(integrator, p, integrator.t)
    nothing
end

# setting up a new simulation for initialize
function (p::AbstractDSSAJumpAggregator)(dj, u, t, integrator) # initialize
    initialize!(p, integrator, u, integrator.p, t)
    register_next_jump_time!(integrator, p, integrator.t)
    u_modified!(integrator,false)
    nothing
end

############################## Generic Routines ###############################

"""
    fill_rates_and_sum!(p::AbstractDSSAJumpAggregator, u, params, t)
Reevaluate all rates and their sum.
"""
function fill_rates_and_sum!(p::AbstractDSSAJumpAggregator, u, params, t)
    sum_rate = zero(typeof(p.sum_rate))

    # mass action jumps
    majumps   = p.ma_jumps
    cur_rates = p.cur_rates
    @inbounds for i in 1:get_num_majumps(majumps)
        cur_rates[i] = evalrxrate(u, i, majumps)
        sum_rate    += cur_rates[i]
    end

    # constant rates
    rates = p.rates
    idx   = get_num_majumps(majumps) + 1
    @inbounds for rate in rates
        cur_rates[idx] = rate(u, params, t)
        sum_rate += cur_rates[idx]
        idx += 1
    end

    p.sum_rate = sum_rate
    nothing
end

"""
    calculate_jump_rate(ma_jumps, rates, u, params, t, rx)
Recalculate the rate for the jump with index `rx`.
"""
@inline function calculate_jump_rate(ma_jumps, num_majumps, rates, u, params, t, rx)
    if rx <= num_majumps
        return evalrxrate(u, rx, ma_jumps)
    else
        @inbounds return rates[rx - num_majumps](u, params, t)
    end
end

@inline function dt_delay_generation!(p, integrator)
    next_jump = p.next_jump
    ttnj = p.time_to_next_jump
    @unpack delay_trigger_set, delay_interrupt_set = integrator.delayjumpsets
    if p.next_delay != nothing || any(set->next_jump in set, [delay_interrupt_set, delay_trigger_set])
        find_next_delay_dt!(p, integrator)
    else
        p.dt_delay -= ttnj
    end
end

@inline function update_state_delay!(p::AbstractDSSAJumpAggregator, integrator, u, t)
    @unpack ma_jumps, next_jump, next_delay, num_next_delay, next_jump_time, time_to_next_jump = p
    @unpack delay_interrupt, delay_trigger, delay_trigger_set, delay_interrupt_set, delay_complete = integrator.delayjumpsets

    ttnj = time_to_next_jump
    if next_delay == nothing
        # update_state !
        num_ma_rates = get_num_majumps(ma_jumps)
        if next_jump <= num_ma_rates # is next jump a mass action jump
            if u isa SVector
                integrator.u = executerx(u, next_jump, ma_jumps)
            else
                @inbounds executerx!(u, next_jump, ma_jumps)
            end
        else
            idx = next_jump - num_ma_rates
            @inbounds p.affects![idx](integrator)
        end
        # shift delay channel !
        shift_delay_channel!(integrator.de_chan,ttnj) #shift delays in all channel
        if next_jump in delay_interrupt_set
            update_delay_interrupt!(p, integrator)
        elseif next_jump in delay_trigger_set
            update_delay_trigger!(p, integrator)
        end
    else
        shift_delay_channel!(integrator.de_chan, ttnj)
        update_delay_channel!(integrator.de_chan)
        update_delay_complete!(p, integrator)
    end
    # save jump that was just executed
    p.prev_jump = next_jump
    nothing
end

"""
    function compare_delay!(p::AbstractDSSAJumpAggregator, de_chan, dt_delay, dt_reaction, t)
Compare dt_delay and dt_reaction.
"""
function compare_delay!(p::AbstractDSSAJumpAggregator, de_chan, dt_delay, dt_reaction, t)
    if  dt_reaction < dt_delay
        ttnj = dt_reaction
        next_delay = nothing
        num_next_delay = nothing
    else
        ttnj = dt_delay
        next_delay, num_next_delay = find_next_delay_vec(de_chan, ttnj)
    end
    p.time_to_next_jump = ttnj
    p.next_delay = next_delay
    p.num_next_delay = num_next_delay
    @fastmath p.next_jump_time = t + ttnj
    nothing
end

"""
    function find_next_delay_dt!(p, integrator)
Find the minimal dt_delay in all delay channels.
"""
function find_next_delay_dt!(p, integrator)
    de_chan = integrator.de_chan
    p.dt_delay = find_minimun_dt_delay(de_chan)
    nothing
end

function find_minimun_dt_delay(de_chan::Vector{Vector{T}}) where {T}
    val_vec = Vector{T}(undef,length(de_chan))
    @inbounds for i in eachindex(de_chan)
        val_vec[i] = isempty(de_chan[i]) ? typemax(T) : minimum(de_chan[i])
    end
    minimum(val_vec)
end

function find_next_delay_num!(p, de_chan::Vector{Vector{T}}) where {T}
    @label restart
    val = find_minimun_dt_delay(de_chan)
    if val < eps(T)
        shift_delay_channel!(de_chan, eps(T))
        update_delay_channel!(de_chan)
        @goto restart
    end
    p.next_delay, p.num_next_delay = find_next_delay_vec(de_chan, val)
    p.next_delay_time = val
    nothing
end

@inline function shift_delay_channel!(de_chan::Vector{Vector{T1}},ttnj::T2) where {T1<:Real,T2<:Real}
    @inbounds for idx in eachindex(de_chan)
        for j in eachindex(de_chan[idx])
            de_chan[idx][j] -= ttnj
        end
    end
end

@inline function update_delay_channel!(de_chan::Vector{Vector{T}}) where {T<:Real}
    @inbounds for idx in eachindex(de_chan)
        filter!(x->x>zero(T), de_chan[idx])
    end
end
"""
    function update_delay_interrupt!(p, integrator)
"""
function update_delay_interrupt!(p, integrator)
    delay_interrupt = integrator.delayjumpsets.delay_interrupt
    execute_delay_interrupt!(delay_interrupt[p.next_jump], p, integrator)
end
function execute_delay_interrupt!(delay_interrupt_affect!::Function, p, integrator)
    delay_interrupt_affect!(integrator, p.rng)
end

"""
    function update_delay_trigger!(p, integrator)
"""
function update_delay_trigger!(p, integrator)
    delay_trigger = integrator.delayjumpsets.delay_trigger

    execute_delay_trigger!(delay_trigger[p.next_jump], p, integrator)
end
function execute_delay_trigger!(delay_trigger_affect!::Function, p, integrator)
    delay_trigger_affect!(integrator, p.rng)
end
function execute_delay_trigger!(delay_trigger_affect!::Vector{Pair{Int64,T}}, p, integrator) where {T}
    for (chan_idx, τ) in delay_trigger_affect!
        append!(integrator.de_chan[chan_idx], τ)
    end
end

"""
    function update_delay_complete!(p, integrator)
This function modifies `integrator.u` and `integrator.de_chan` upon delay completion
"""
function update_delay_complete!(p, integrator)
    delay_complete = integrator.delayjumpsets.delay_complete
    next_delay_vec, num_next_delay_vec = p.next_delay, p.num_next_delay
    @inbounds for j in eachindex(next_delay_vec)
        next_delay = next_delay_vec[j]
        num_next_delay = num_next_delay_vec[j]
        execute_delay_complete!(delay_complete[next_delay], num_next_delay, integrator, p.rng)
    end
end

function execute_delay_complete!(delay_complete::Vector{Pair{Int64,Int64}}, num_next_delay::Int64, integrator, rng)
    @inbounds for (i, ξ) in delay_complete
        if integrator.u isa SVector
            integrator.u = setindex(integrator.u, integrator.u[i] + num_next_delay*ξ, i)
            # integrator.u = setindex!!(integrator.u, integrator.u[i] + num_next_delay*ξ, i)
        else
            integrator.u[i] += num_next_delay*ξ
        end
    end
end
function execute_delay_complete!(delay_complete::Function, num_next_delay::Int64, integrator, rng)
    for _ in 1:num_next_delay
        delay_complete(integrator, rng)
    end
end

"""
    function dep_gr_delay(p, integrator)
Generate delay dependency graph
input::Int next_delay idx
output::Dict next_delay idx => reactions need to be updated
"""
function dep_gr_delay(p::AbstractDSSAJumpAggregator, integrator)
    # num_delay_chan = length(integrator.de_chan)
    dict_ = Dict{Int,Vector{Int}}()
    dict_complete = integrator.delayjumpsets.delay_complete
    # var_to_jumps = var_to_jumps_map(length(integrator.u),p.ma_jumps)
    @inbounds for key in keys(dict_complete)
        delay_complete_action = dict_complete[key]
        if typeof(delay_complete_action)<:Vector{Pair{Int,Int}}
            vars = first.(delay_complete_action)
        else
            vars = vec(1:length(integrator.u))
        end
        jumps = unique(vars_to_jumps_delay(p, vars))
        push!(dict_, key=>jumps)
    end
    dict_
end
function vars_to_jumps_delay(p, vars)
    jumps = []
    for i in eachindex(p.dep_gr)
        for var in vars
            if var in p.dep_gr[i]
                push!(jumps, i)
            end
        end
    end
    return jumps
end

"""
    find_next_delay_vec(A::Vector{Vector{T}}, x::T)
Returns two vectors. The first is the position indices; `position_indices` of vector `A`, and the second is the vector `num_in_vec` composed of the occurrence of value `x`.
Used in finding the minimal dt_delay in various delay channel.
# Examples
```julia-repl
julia> A =  [[0.09,0.09,0.1],[0.3,0.09,0.1],[0.09]]
3-element Vector{Vector{Float64}}:
 [0.09, 0.09, 0.1]
 [0.3, 0.09, 0.1]
 [0.09]
 julia> x=0.09
 0.09
julia> find_next_delay_vec(A, x)
([1, 2, 3], [2, 1, 1])
```
"""
function find_next_delay_vec(de_chan::Vector{Vector{T}}, ttnj::T) where {T}
    position_indices = Vector{Int64}()
    num_in_vec = Vector{Int64}()
    for idx in eachindex(de_chan)
        if ttnj in de_chan[idx]
            append!(position_indices, idx)
            append!(num_in_vec, count(==(ttnj), de_chan[idx]))
        end
    end
    return position_indices, num_in_vec
end

## decrepit
# """
#     function find_num_in_vec(A::Vector{Vector{T}}, position_indices::Vector{Int64}, x::T)

# Find the occurrence of value `x` in each vector, according to the corresponding position indices specified by the element in the `position_indices` vector in the given vector `A`.

# # Examples
# ```julia-repl
# julia> A =  [[0.09,0.09,0.1],[0.3,0.09,0.1],[0.09]]
# 3-element Vector{Vector{Float64}}:
#  [0.09, 0.09, 0.1]
#  [0.3, 0.09, 0.1]
#  [0.09]
# julia> position_indices =  [1,2]
# 2-element Vector{Int64}:
#  1
#  2
# julia> find_num_in_vec(A, position_indices, 0.09)
# 2-element Vector{Int64}:
# 2
# 1
# ```
# """
# @inline function find_num_in_vec(A::Vector{Vector{T}}, position_indices::Vector{Int64}, x::T) where {T}
#     number_in_vec = Vector{Int64}(undef, length(position_indices))
#     @inbounds for i in eachindex(position_indices)
#         number_in_vec[i] = count(==(x),A[position_indices[i]])
#     end
#     return number_in_vec
# end
# function find_next_delay_vec(de_chan::Vector{Vector{T}}, ttnj::T) where {T}
#     position_indices = findall(de_chan->ttnj in de_chan, de_chan)
#     num_in_vec = find_num_in_vec(de_chan, position_indices, ttnj)
#     return position_indices, num_in_vec
# end


text:
using LaptimeSimulation
using Test

@testset "LaptimeSimulation.jl" begin
    # Write your tests here.
end


text:
using Plots

# Model
include_model("quadruped")
model = free_time_model(model)

include(joinpath(pwd(), "models/visualize.jl"))
vis = Visualizer()
render(vis)

# Horizon
T = 31
T_fix = 5

# Time step
tf = 0.625
h = tf / (T - 1)

# Permutation matrix
perm = @SMatrix [1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0;
                 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0;
				 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0;
				 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0;
				 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0;
				 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0;
				 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0;
				 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0;
				 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0;
				 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0;
				 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0]

function ellipse_traj(x_start, x_goal, z, T)
	dist = x_goal - x_start
	a = 0.5 * dist
	b = z
	z̄ = 0.0
	x = range(x_start, stop = x_goal, length = T)
	z = sqrt.(max.(0.0, (b^2) * (1.0 .- ((x .- (x_start + a)).^2.0) / (a^2.0))))
	return x, z
end

function initial_configuration(model::Quadruped, θ1, θ2, θ3)
    q1 = zeros(model.nq)
    q1[3] = pi / 2.0
    q1[4] = -θ1
    q1[5] = θ2

    q1[8] = -θ1
    q1[9] = θ2

    q1[2] = model.l_thigh1 * cos(q1[4]) + model.l_calf1 * cos(q1[5])

    q1[10] = -θ3
    q1[11] = acos((q1[2] - model.l_thigh2 * cos(q1[10])) / model.l_calf2)

    q1[6] = -θ3
    q1[7] = acos((q1[2] - model.l_thigh2 * cos(q1[6])) / model.l_calf2)

    return q1
end

θ1 = pi / 4.0
θ2 = pi / 4.0
θ3 = pi / 3.0

q1 = initial_configuration(model, θ1, θ2, θ3)
visualize!(vis, model, [q1])

# feet positions
pr1 = kinematics_2(model, q1, body = :calf_1, mode = :ee)
pr2 = kinematics_2(model, q1, body = :calf_2, mode = :ee)

pf1 = kinematics_3(model, q1, body = :calf_3, mode = :ee)
pf2 = kinematics_3(model, q1, body = :calf_4, mode = :ee)

strd = 2 * (pr1 - pr2)[1]
qT = Array(perm) * copy(q1)
qT[1] += 0.5 * strd

# torso height
pt = kinematics_1(model, q1, body = :torso, mode = :com)

zh = 0.05
xr2 = [pr2[1] for t = 1:T]
zr2 = [pr2[2] for t = 1:T]
pr2_ref = [[xr2[t]; zr2[t]] for t = 1:T]

xf1 = [pf1[1] for t = 1:T]
zf1 = [pf1[2] for t = 1:T]
pf1_ref = [[xf1[t]; zf1[t]] for t = 1:T]

xr1_el, zr1_el = ellipse_traj(pr1[1], pr1[1] + strd, zh, T - T_fix)
xr1 = [[xr1_el[1] for t = 1:T_fix]..., xr1_el...]
zr1 = [[zr1_el[1] for t = 1:T_fix]..., zr1_el...]
pr1_ref = [[xr1[t]; zr1[t]] for t = 1:T]

xf2_el, zf2_el = ellipse_traj(pf2[1], pf2[1] + strd, zh, T - T_fix)
xf2 = [[xf2_el[1] for t = 1:T_fix]..., xf2_el...]
zf2 = [[zf2_el[1] for t = 1:T_fix]..., zf2_el...]
pf2_ref = [[xf2[t]; zf2[t]] for t = 1:T]

tr = range(0, stop = tf, length = T)
plot(tr, hcat(pr1_ref...)')
plot!(tr, hcat(pf1_ref...)')

plot(tr, hcat(pr2_ref...)')
plot!(tr, hcat(pf2_ref...)')

# Bounds

# control
# ul <= u <= uu
_uu = Inf * ones(model.m)
_uu[model.idx_u] .= 33.5 * h
_uu[end] = 2.0 * h
_ul = zeros(model.m)
_ul[model.idx_u] .= -33.5 * h
_ul[end] = 0.75 * h
ul, uu = control_bounds(model, T, _ul, _uu)

xl, xu = state_bounds(model, T,
    x1 = [Inf * ones(model.nq); q1],
    xT = [Inf * ones(model.nq); qT[1]; Inf * ones(model.nq - 1)])

# Objective
include_objective(["velocity", "nonlinear_stage", "control_velocity"])
q_ref = linear_interpolation(q1, qT, T+1)
render(vis)
visualize!(vis, model, q_ref, Δt = h)
x0 = configuration_to_state(q_ref)

# penalty on slack variable
obj_penalty = PenaltyObjective(1.0e4, model.m - 1)

# quadratic tracking objective
# Σ (x - xref)' Q (x - x_ref) + (u - u_ref)' R (u - u_ref)
obj_control = quadratic_time_tracking_objective(
    [1.0 * Diagonal(1.0e-5 * ones(model.n)) for t = 1:T],
    [1.0 * Diagonal([1.0e-3 * ones(model.nu)..., 1.0e-3 * ones(model.nc + model.nb)..., zeros(model.m - model.nu - model.nc - model.nb)...]) for t = 1:T-1],
    [[qT; qT] for t = 1:T],
    [zeros(model.m) for t = 1:T],
    1.0)

# quadratic velocity penalty
#Σ v' Q v
v_penalty = 0.0 * ones(model.nq)
obj_velocity = velocity_objective(
    [h * Diagonal(v_penalty) for t = 1:T-1],
    model.nq,
    h = h,
    idx_angle = collect([3, 4, 5, 6, 7, 8, 9, 10, 11]))

obj_ctrl_velocity = control_velocity_objective(Diagonal([1.0e-3 * ones(model.nu)..., 1.0e-3 * ones(model.nc + model.nb)..., zeros(model.m - model.nu - model.nc - model.nb)...]))

function l_stage(x, u, t)
	q1 = view(x, 1:11)
	q2 = view(x, 11 .+ (1:11))
    J = 0.0

	# torso height
    J += 100.0 * (kinematics_1(model, q1, body = :torso, mode = :ee)[2] - kinematics_1(model, view(x0[t], 1:11), body = :torso, mode = :com)[2])^2.0
	J += 100.0 * (kinematics_1(model, q2, body = :torso, mode = :ee)[2] - kinematics_1(model, view(x0[t], 1:11), body = :torso, mode = :com)[2])^2.0

	J += 100.0 * (q1[2] - x0[1][2])^2.0
	J += 100.0 * (q2[2] - x0[1][2])^2.0

    if true
		J += 10000.0 * sum((pr1_ref[t] - kinematics_2(model, q1, body = :calf_1, mode = :ee)).^2.0)
	    J += 10000.0 * sum((pf2_ref[t] - kinematics_3(model, q1, body = :calf_4, mode = :ee)).^2.0)
	end

    return J
end

l_terminal(x) = 0.0
obj_shaping = nonlinear_stage_objective(l_stage, l_terminal)

obj = MultiObjective([obj_penalty,
                      obj_control,
                      obj_velocity,
					  obj_shaping,
					  obj_ctrl_velocity])
# Constraints
include_constraints(["stage", "contact", "free_time", "loop"])
function pinned1!(c, x, u, t)
    q = view(x, 1:11)
    c[1:2] = pr2_ref[t] - kinematics_2(model, q, body = :calf_2, mode = :ee)
    c[3:4] = pf1_ref[t] - kinematics_3(model, q, body = :calf_3, mode = :ee)
    nothing
end

function pinned2!(c, x, u, t)
    q = view(x, 1:11)
	c[1:2] = pr1_ref[t] - kinematics_2(model, q, body = :calf_1, mode = :ee)
    c[3:4] = pf2_ref[t] - kinematics_3(model, q, body = :calf_4, mode = :ee)
    nothing
end

n_stage = 4
t_idx1 = vcat([t for t = 1:T])
t_idx2 = vcat([1:T_fix]...)
con_pinned1 = stage_constraints(pinned1!, n_stage, (1:0), t_idx1)
con_pinned2 = stage_constraints(pinned2!, n_stage, (1:0), t_idx2)

con_contact = contact_constraints(model, T)
con_free_time = free_time_constraints(T)
con_loop = loop_constraints(model, collect([(2:model.nq)...,
	(nq .+ (2:model.nq))...]), 1, T, perm = Array(cat(perm, perm, dims = (1,2))))
con = multiple_constraints([con_contact, con_loop,
    con_free_time, con_pinned1, con_pinned2])

# Problem
prob = trajectory_optimization_problem(model,
               obj,
               T,
               xl = xl,
               xu = xu,
               ul = ul,
               uu = uu,
               con = con)

# trajectory initialization
u0 = [[1.0e-3 * rand(model.m - 1); h] for t = 1:T-1] # random controls

# Pack trajectories into vector
z0 = pack(x0, u0, prob)

# Solve
@time z̄, info = solve(prob, copy(z0),
    nlp = :ipopt,
    tol = 1.0e-3, c_tol = 1.0e-3,
	max_iter = 2000,
    time_limit = 60 * 2, mapl = 5)

check_slack(z̄, prob)
x̄, ū = unpack(z̄, prob)
q̄ = state_to_configuration(x̄)
_tf, _t, h̄ = get_time(ū)
@show h̄[1]

vis = Visualizer()
render(vis)
visualize!(vis, model, state_to_configuration(x̄), Δt = h̄[1])

# check foot trajectories
_pr1_ref = [kinematics_2(model, q, body = :calf_1, mode = :ee) for q in q̄]
_pf1_ref = [kinematics_3(model, q, body = :calf_3, mode = :ee) for q in q̄]

_pr2_ref = [kinematics_2(model, q, body = :calf_2, mode = :ee) for q in q̄]
_pf2_ref = [kinematics_3(model, q, body = :calf_4, mode = :ee) for q in q̄]

plot(hcat(pr1_ref...)', width = 2.0, color = :black)
plot!(hcat(_pr1_ref...)', color = :red)

plot(hcat(pf1_ref...)', width = 2.0, color = :black)
plot!(hcat(_pf1_ref...)', color = :red)

plot(hcat(pr2_ref...)', width = 2.0, color = :black)
plot!(hcat(_pr2_ref...)', color = :red)

plot(hcat(pf2_ref...)', width = 2.0, color = :black)
plot!(hcat(_pf2_ref...)', color = :red)

# check control trajectory
plot(hcat(ū...)[1:model.nu, :]', linetype = :steppost)

# unpack solution
q = state_to_configuration(x̄)
u = [u[model.idx_u] for u in ū]
γ = [u[model.idx_λ] for u in ū]
b = [u[model.idx_b] for u in ū]
ψ = [u[model.idx_ψ] for u in ū]
η = [u[model.idx_η] for u in ū]
hm = mean(h̄)
μm = model.μ
#
# perm4 = [0.0 1.0 0.0 0.0;
#          1.0 0.0 0.0 0.0;
# 		 0.0 0.0 0.0 1.0;
# 		 0.0 0.0 1.0 0.0]
#
# perm8 = [0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0;
#          0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0;
# 		 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0;
# 		 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0;
# 		 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0;
# 		 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0;
# 		 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0;
# 		 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0]
#
# function mirror_gait(q, u, γ, b, ψ, η, T)
# 	qm = [deepcopy(q)...]
# 	um = [deepcopy(u)...]
# 	γm = [deepcopy(γ)...]
# 	bm = [deepcopy(b)...]
# 	ψm = [deepcopy(ψ)...]
# 	ηm = [deepcopy(η)...]
#
# 	stride = zero(qm[1])
# 	@show stride[1] = q[T+1][1] - q[2][1]
# 	@show 0.5 * strd
#
# 	for t = 1:T-1
# 		push!(qm, Array(perm) * q[t+2] + stride)
# 		push!(um, perm8 * u[t])
# 		push!(γm, perm4 * γ[t])
# 		push!(bm, perm8 * b[t])
# 		push!(ψm, perm4 * ψ[t])
# 		push!(ηm, perm8 * η[t])
# 	end
#
# 	return qm, um, γm, bm, ψm, ηm
# end
#
# qm, um, γm, bm, ψm, ηm = mirror_gait(q, u, γ, b, ψ, η, T)
#
# @save joinpath(@__DIR__, "quadruped_mirror_gait.jld2") qm um γm bm ψm ηm μm hm
#
# plot(hcat(q...)', color = :black, width = 2.0, label = "")
# plot!(hcat(qm...)', color = :red, width = 1.0, label = "")
#
# plot(hcat(u...)', color = :black, width = 2.0, label = "", linetype = :steppost)
# plot!(hcat(um...)', color = :red, width = 1.0, label = "", linetype = :steppost)
#
# plot(hcat(γ...)', color = :black, width = 2.0, label = "", linetype = :steppost)
# plot!(hcat(γm...)', color = :red, width = 1.0, label = "", linetype = :steppost)
#
# plot(hcat(b...)', color = :black, width = 2.0, label = "", linetype = :steppost)
# plot!(hcat(bm...)', color = :red, width = 1.0, label = "", linetype = :steppost)
#
# function get_q_viz(q̄)
# 	q_viz = [q̄...]
# 	shift_vec = zeros(model.nq)
# 	shift_vec[1] = q̄[end][1]
# 	for i = 1:5
# 		q_update = [q + shift_vec for q in q̄[2:end]]
# 		push!(q_viz, q_update...)
# 		shift_vec[1] = q_update[end][1]
# 	end
#
# 	return q_viz
# end
#
# vis = Visualizer()
# render(vis)
# visualize!(vis, model,
# 	qm,
# 	Δt = h̄[1])

max_issues_repo_name: thowell/DirectMotionPlanning


text:
# Builtin functions not covered elsewhere

"""
    pyprint(...)

Equivalent to `print(...)` in Python.
"""
pyprint(args...; kwargs...) = (pydel!(pybuiltins.print(args...; kwargs...)); nothing)
export pyprint

"""
    pyhelp([x])

Equivalent to `help(x)` in Python.
"""
pyhelp() = (pydel!(pybuiltins.help()); nothing)
pyhelp(x) = (pydel!(pybuiltins.help(x)); nothing)
export pyhelp

"""
    pyall(x)

Equivalent to `all(x)` in Python.
"""
function pyall(x)
    y = pybuiltins.all(x)
    z = pybool_asbool(y)
    pydel!(y)
    z
end
export pyall

"""
    pyany(x)

Equivalent to `any(x)` in Python.
"""
function pyany(x)
    y = pybuiltins.any(x)
    z = pybool_asbool(y)
    pydel!(y)
    z
end
export pyany

"""
    pycallable(x)

Equivalent to `callable(x)` in Python.
"""
function pycallable(x)
    y = pybuiltins.callable(x)
    z = pybool_asbool(y)
    pydel!(y)
    z
end
export pycallable


text:
export WhithamBoussinesq

"""
    WhithamBoussinesq(param;kwargs)

Define an object of type `AbstractModel` in view of solving the initial-value problem for
a Boussinesq-type model with full-dispersion property.

# Argument
`param` is of type `NamedTuple` and must contain
- dimensionless parameters `ϵ` (nonlinearity) and `μ` (dispersion);
- numerical parameters to construct the mesh of collocation points as `mesh = Mesh(param)`

## Optional keyword arguments
- `Boussinesq`: if `true` (default is `false`), compute the standard Boussinesq system instead (see `Boussinesq(param;kwargs)`);
- a parameter `α` which determines the model solved:
    - If `α = 1` (default), then the model has been introduced in [Dinvay, Dutykh and Kalisch](https://doi.org/10.1016/j.apnum.2018.09.016);
    - If `α = 1/2`, then the model is a quasilinear version;
    - If `α < 1/2`, then expect instabilities stemming from ill-posedness of the model.
- `ktol`: tolerance of the low-pass Krasny filter (default is `0`, i.e. no filtering);
- `dealias`: dealiasing with Orlicz rule `1-dealias/(dealias+2)` (default is `0`, i.e. no dealiasing);
- `label`: a label for future references (default is `"Whitham-Boussinesq"`);
- `verbose`: prints information if `true` (default is `true`).

# Return values
Generate necessary ingredients for solving an initial-value problem via `solve!`:
1. a function `WhithamBoussinesq.f!` to be called in explicit time-integration solvers;
2. a function `WhithamBoussinesq.mapto` which from `(η,v)` of type `InitialData` provides the raw data matrix on which computations are to be executed.
3. a function `WhithamBoussinesq.mapfro` which from such data matrix returns the Tuple of real vectors `(η,v)`, where
    - `η` is the surface deformation;
    - `v` is the derivative of the trace of the velocity potential.

"""
mutable struct WhithamBoussinesq <: AbstractModel

	label   :: String
	f!		:: Function
	mapto	:: Function
	mapfro	:: Function

    function WhithamBoussinesq(param::NamedTuple;Boussinesq=false,
								α = 1, a = -1/3, b = 1/3,
								dealias = 0,
								ktol	= 0,
								label 	= nothing,
								verbose	=true)

		# Set up
		μ 	= param.μ
		ϵ 	= param.ϵ

		if Boussinesq == true
			if label == nothing label = "Boussinesq" end
			info_param = "a=$a, b=$b, c=0 and d=$b"
		else
			if label == nothing label = "Whitham-Boussinesq" end
			info_param = "α=$α"
		end

		if verbose # Print information
			info = "Build the $label model with $info_param.\n"
			info *= "Shallowness parameter μ=$μ, nonlinearity parameter ϵ=$ϵ.\n"
			if dealias == 0
				info *= "No dealiasing. "
			else
				info *= "Dealiasing with Orszag's rule adapted to power $(dealias + 1) nonlinearity. "
			end
			if ktol == 0
				info *= "No Krasny filter. "
			else
				info *= "Krasny filter with tolerance $ktol."
			end
			info *= "\nShut me up with keyword argument `verbose = false`."
			@info info
		end

		# Pre-allocate data
		mesh = Mesh(param)
		x 	= mesh.x
		k = mesh.k
		∂ₓ	=  1im * k
		if Boussinesq==false
			F₁ 	= tanh.(sqrt(μ)*abs.(k))./(sqrt(μ)*abs.(k))
			F₁[1] 	= 1
			F₂ = F₁.^α
		else
			F₂ = 1 ./(1 .+μ*b*abs.(k).^2)
			F₁ 	= (1 .-μ*a*abs.(k).^2).*(F₂.^2)
		end

		if dealias == 0
			Π⅔ 	= ones(size(k)) # no dealiasing (Π⅔=Id)
		else
			K = mesh.kmax * (1-dealias/(2+dealias))
			Π⅔ 	= abs.(k) .<= K # Dealiasing low-pass filter
		end
		η = zeros(Float64, mesh.N)
        v = zeros(Float64, mesh.N)
		fftη = zeros(Complex{Float64}, mesh.N)
        fftv = zeros(Complex{Float64}, mesh.N)

		# Evolution equations are ∂t U = f(U)
		function f!(U)

		    fftv .= U[:,2]
			fftη .= F₂.*U[:,2]
			v .= real(ifft(fftη))
			fftη .= U[:,1]
		   	η .= real(ifft(U[:,1]))

		   	U[:,1] .= -∂ₓ.*(F₁.*fftv.+ϵ*Π⅔.*F₂.*fft(η.*v))
		   	U[:,2] .= -∂ₓ.*(fftη.+ϵ/2*Π⅔.*fft(v.^2))
			U[abs.(U).< ktol ].=0

		end

		# Build raw data from physical data.
		# Discrete Fourier transform with, possibly, dealiasing and Krasny filter.
		function mapto(data::InitialData)
			U = [Π⅔ .* fft(data.η(x)) Π⅔ .*fft(data.v(x))]
			U[abs.(U).< ktol ].=0
			return U
		end

		# Return `(η,v)`, where
		# - `η` is the surface deformation;
		# - `v` is the derivative of the trace of the velocity potential.
		# Inverse Fourier transform and takes the real part.
		function mapfro(U)
			real(ifft(U[:,1])),real(ifft(U[:,2]))
		end

        new(label, f!, mapto, mapfro)
    end
end


text:
using Dash
using DashHtmlComponents
using DashCoreComponents
# using DashCytoscape

app = dash(external_stylesheets = ["https://codepen.io/chriddyp/pen/bWLwgP.css"])

app.layout = html_div() do
    html_h1("Hello Dash"),
        dcc_tabs(id="stage-tabs", value="tab-model") do
            dcc_tab(value="tab-model", label="Model") do
                html_div("Model")
            end,
            dcc_tab(value="tab-identify", label="Identify ") do

            end,
            dcc_tab(value="tab-estimate", label="Estimate") do

            end,
            dcc_tab(value="tab-verify", label="Verify") do

            end
        end,
        html_div(id="tab-content")
    #     cyto_cytoscape(
    #     id="cytoscape-elements-boolean",
    #     layout=(name="preset",),
    #     style=(width="100%",
    #            height="400px"),
    #     elements=[(
    #                 data=(id="one", label="Locked"),
    #                 position=(x=75, y=75),
    #                 locked=true
    #             ),
    #             (
    #                 data=(id="two", label="Selected"),
    #                 position=(x=75, y=200),
    #                 selected=true
    #             ),
    #             (
    #                 data=(id="three", label="Not Selectable"),
    #                 position=(x=200, y=75),
    #                 selectable=false
    #             ),
    #             (
    #                 data=(id="four", label="Not Grabbable"),
    #                 position=(x=200, y=200),
    #                 grabbable=false
    #             ),
    #             (
    #                 data=(source="one", target="two"),
    #             ),
    #             (
    #                 data=(source="two", target="three"),
    #             ),
    #             (
    #                 data=(source="three", target="four"),
    #             ),
    #             (
    #                 data=(source="two", target="four"),
    #             )
    #     ]
    # )
    # dcc_graph(
    #     id = "example-graph",
    #     figure = (
    #         data = [
    #             (x = [1, 2, 3], y = [4, 1, 2], type = "bar", name = "SF"),
    #             (x = [1, 2, 3], y = [2, 4, 5], type = "bar", name = "Montréal"),
    #         ],
    #         layout = (title = "Dash Data Visualization",)
    #     )
    # )
end

callback!(app, Output("tab-content", "children"), Input("stage-tabs", "value")) do tab
    # extra = "-"
    if tab == "tab-model"
        return show_model_creation_content()
    else
        # return tab
        return html_h3("Tab $tab ... ")
    end
    # # isnothing() && throw(PreventUpdate())
end

function show_model_creation_content()
    extra = "uber model"
    return extra
end

run_server(app, "0.0.0.0", 8080, debug=true)

max_stars_repo_name: Magier/Aetia

max_issues_repo_name: Magier/Aetia

max_forks_repo_name: Magier/Aetia


text:
module POPT

import PTChain
@everywhere import PTChain

function minTuple(a, b)
	a[1] < b[1] ? a : b
end

function pSimAnneal(runs, start_x, optimFunction, temps, proposalFunction, itersAtTemp)
	@sync @parallel (minTuple) for i = 1:runs
		PTChain.simulatedAnnealing(start_x, optimFunction, temps, proposalFunction, itersAtTemp)
	end
end

function blockingGreedy(start_x, optimFunction, temps, proposalFunction, iterBetweenSwaps, num_swaps)
	sort!(temps)
	num_temps = length(temps)
	out = SharedArray(repmat(start_x, 1, num_temps))
	func_values = SharedArray(zeros(num_temps))
	func_values[:] = optimFunction(start_x)
	min_val = SharedArray(zeros(num_temps))
	min_val[:] = optimFunction(start_x)
	plot_vals = zeros(num_temps, num_swaps)
	min_x = copy(out)
	for i = 1:num_swaps
	    @sync @parallel for j = 1:num_temps
	        fv, out[:, j], best_val, best_x = PTChain.runChain(out[:,j], optimFunction, temps[j], proposalFunction, iterBetweenSwaps)
			func_values[j] = fv
			if best_val < min_val[j]
	            min_val[j] = copy(best_val)
	            min_x[:,j] = copy(best_x)
	        end
	    end
	    plot_vals[:, i] = func_values
	    out = greedyswapShared(func_values, out)
	end

	return min_val, min_x, plot_vals
end

function epchains(start_x, optimFunction, temps, proposalFunction, iterBetweenSwaps, num_swaps)
	sort!(temps)
	num_temps = length(temps)
	out = SharedArray(repmat(start_x, 1, num_temps))
	func_values = SharedArray(zeros(num_temps))
	func_values[:] = optimFunction(start_x)
	min_val = SharedArray(zeros(num_temps))
	min_val[:] = optimFunction(start_x)
	plot_vals = zeros(num_temps, num_swaps)
	min_x = copy(out)
	@sync @parallel for j = 1:num_temps
		for i = 1:num_swaps
	        func_values[j], out[:, j], best_val, best_x = PTChain.runChain(out[:,j], optimFunction, temps[j], proposalFunction, iterBetweenSwaps)
	        if best_val < min_val[j]
	            min_val[j] = copy(best_val)
	            min_x[:,j] = copy(best_x)
	        end
	    end
	end

	return min_val, min_x, plot_vals
end

function pINSChain(start_x, optimFunction, temps, proposalFunction, iterBetweenSwaps, num_swaps, swap_schemes)
    sort!(temps)
    num_temps = length(temps)
    num_swap_schemes = length(swap_schemes)
    out = SharedArray(repmat(start_x, 1, num_temps))
	func_values = SharedArray(zeros(num_temps))
    func_values[:] = optimFunction(start_x)
	min_val = SharedArray(zeros(num_temps))
	min_val[:] = optimFunction(start_x)
    plot_vals = zeros(num_temps, num_swaps)
    min_x = copy(out)
    current_ss = 1
    for i = 1:num_swaps
    	swap_scheme = swap_schemes[current_ss]
    	num_interactions = length(swap_scheme)
    	@sync @parallel for j = 1:num_interactions
    		@sync @parallel for k = swap_scheme[j]
				func_values[k], out[:, k], best_val, best_x = PTChain.runChain(out[:,k], optimFunction, temps[k], proposalFunction, iterBetweenSwaps)
				if best_val < min_val[k]
					min_val[k] = copy(best_val)
					min_x[:,k] = copy(best_x)
				end
    		end
            func_values[swap_scheme[j]], out[:, swap_scheme[j]] = greedyswap(func_values[swap_scheme[j]], out[:, swap_scheme[j]])
    	end
        if current_ss == num_swap_schemes
            current_ss = 1
        else
            current_ss += 1
        end
    end
	return min_val, min_x, plot_vals
end

function runblock(start_x, optimFunction, temps, proposalFunction, iterBetweenSwaps, num_swaps)
	c_1 = RemoteChannel(1)
	c_2 = RemoteChannel(2)
	res_1 = @spawn blockTwo(start_x, optimFunction, temps[1], proposalFunction, iterBetweenSwaps, num_swaps, c_1, c_2)
	res_2 = @spawn blockTwo(start_x, optimFunction, temps[2], proposalFunction, iterBetweenSwaps, num_swaps, c_2, c_1)
	m_1, b_1, s_1 = fetch(res_1)
	m_2, b_2, s_2 = fetch(res_2)
	return (m_1, m_2, b_1, b_2, s_1, s_2)
end

#goal is to make a non-blocking decentralized implementation
function runnonblock(start_x, optimFunction, temps, proposalFunction, iterBetweenSwaps, num_swaps)
	c_1 = RemoteChannel(1)
	c_2 = RemoteChannel(2)
	res_1 = @spawn nonBlockTwo(start_x, optimFunction, temps[1], proposalFunction, iterBetweenSwaps, num_swaps, c_1, c_2)
	res_2 = @spawn nonBlockTwo(start_x, optimFunction, temps[2], proposalFunction, iterBetweenSwaps, num_swaps, c_2, c_1)
	m_1, b_1, s_1 = fetch(res_1)
	m_2, b_2, s_2 = fetch(res_2)
	return (m_1, m_2, b_1, b_2, s_1, s_2)
end

function nonBlockTwo(start_x, optim, temp, propF, micro_steps, n_swaps, in_c, out_c)
	x = start_x
    min_x = copy(start_x)
    best_val = optim(min_x)
    curr_val = copy(best_val)
    update_pending = false
    out_val = 0
    i = 0
    num_swaps = 0
	while i < n_swaps
	    j = 0
		while j < micro_steps
			#proposal
        	new_x = propF(x)
        	proposal_value = optim(new_x)
			#accept based on temp
        	if ((temp == 0 && curr_val > proposal_value) || (log(rand()) < (curr_val - proposal_value) / temp))
               	x = copy(new_x)
               	curr_val = copy(proposal_value)
          	end
          	if curr_val < best_val
               	best_val = copy(curr_val)
            	min_x = copy(x)
          	end
          	j = j + 1
          	#if waiting for update and in_c has something new available
          	if update_pending && isready(in_c)
          		update_pending = false
          		(swap_temp, in_val, swap_x) = take!(in_c)
          		if (in_val < out_val) && (swap_temp > temp) || (in_val > out_val) && (swap_temp < temp)
          			#println("Swapping $(i): my temp is $(temp) and I'm going from $(out_val) to $(in_val)")
          			#println("throwing out $(j) steps")
          			curr_val = copy(in_val)
          			x = copy(swap_x)
          			j = 0
          			num_swaps = num_swaps + 1
          		end
          	end
		end
		if update_pending
			#println("waiting")
			wait(in_c)
			update_pending = false
          	(swap_temp, in_val, swap_x) = take!(in_c)
          	if (in_val < out_val) && (swap_temp > temp) || (in_val > out_val) && (swap_temp < temp)
          	   	#println("Swapping $(i): my temp is $(temp) and I'm going from $(out_val) to $(in_val)")
          		curr_val = copy(in_val)
          		x = copy(swap_x)
				num_swaps = num_swaps + 1
          		continue
			end
		end
		i = i + 1
		#push current stats to out_c
		update_pending = true
		out_val = copy(curr_val)
		out_v = (temp, out_val, copy(x))
		#println("sending a value of $(out_val)")
		put!(out_c, out_v)
		if update_pending && isready(in_c)
          	update_pending = false
          	(swap_temp, in_val, swap_x) = take!(in_c)
          	if (in_val < out_val) && (swap_temp > temp) || (in_val > out_val) && (swap_temp < temp)
          		#println("Swapping $(i): my temp is $(temp) and I'm going from $(out_val) to $(in_val)")
          		num_swaps = num_swaps + 1
          		curr_val = copy(in_val)
          		x = copy(swap_x)
          		j = 0
          	end
        end
	end
	return best_val, min_x, num_swaps
end

function blockTwo(start_x, optim, temp, propF, micro_steps, n_swaps, in_c, out_c)
	x = start_x
    min_x = copy(start_x)
    best_val = optim(min_x)
    curr_val = copy(best_val)
    out_val = 0
    i = 0
    num_swaps = 0
	while i < n_swaps
	    j = 0
		while j < micro_steps
			#proposal
        	new_x = propF(x)
        	proposal_value = optim(new_x)
			#accept based on temp
        	if (log(rand()) < (curr_val - proposal_value) / temp)
               	x = copy(new_x)
               	curr_val = copy(proposal_value)
          	end
          	if curr_val < best_val
               	best_val = copy(curr_val)
            	min_x = copy(x)
          	end
          	j = j + 1
		end
		i = i + 1
		out_val = copy(curr_val)
		out_v = (temp, out_val, copy(x))
		#println("sending a value of $(out_val)")
		put!(out_c, out_v)
		wait(in_c)
		(swap_temp, in_val, swap_x) = take!(in_c)
        if (in_val < out_val) && (swap_temp > temp) || (in_val > out_val) && (swap_temp < temp)
          	#println("Swapping $(i): my temp is $(temp) and I'm going from $(out_val) to $(in_val)")
          	num_swaps = num_swaps + 1
          	curr_val = copy(in_val)
          	x = copy(swap_x)
        end
	end
	return best_val, min_x, num_swaps
end

function greedyswap(values, x)
    ai = sortperm(values)
    out = x[:, ai]
    values[ai], out
end

function greedyswapShared(values, x)
    ai = sortperm(values)
    out = x[:, ai]
    permute!(values, ai)
    SharedArray(out)
end

end

max_stars_repo_name: jacob-ruth/pt-optimization

max_issues_repo_name: jacob-ruth/pt-optimization

max_forks_repo_name: jacob-ruth/pt-optimization


text:
module LeafOptics

using ClimaCache: HyperspectralAbsorption, HyperspectralRadiation, LeafBiophysics, WaveLengthSet
using DocStringExtensions: METHODLIST
using PkgUtility: H_PLANCK, LIGHT_SPEED, AVOGADRO, numerical∫
using SpecialFunctions: expint
using UnPack: @unpack

# export public types from ClimaCache
export HyperspectralAbsorption, HyperspectralRadiation, LeafBiophysics, WaveLengthSet

# export public functions
export leaf_PAR, leaf_SIF, leaf_spectra!

include("fluorescence.jl" )
include("photon.jl"       )
include("radiation.jl"    )
include("spectra.jl"      )
include("transmittance.jl")

end # module


text:
## State updates ##

function update!(interpreter::ConcreteInterpreter,
                 domain::Domain, state::GenericState, diff::GenericDiff)
    setdiff!(state.facts, diff.del)
    union!(state.facts, diff.add)
    vals = [evaluate(domain, state, v) for v in values(diff.ops)]
    for (term, val) in zip(keys(diff.ops), vals)
        set_fluent!(state, val, term)
    end
    return state
end


text:
using Pkg
Pkg.activate(@__DIR__)

using Documenter, Agents, Distributions, Random, StatsBase, Plots
using EvoDynamics
using Literate

# %% Literate convertion
indir = joinpath(@__DIR__, "..", "examples")
outdir = joinpath(@__DIR__, "src")
for file in ("example1.jl", "example2.jl")
	Literate.markdown(joinpath(indir, file), outdir)
end

# %%
makedocs(modules = [EvoDynamics,],
sitename= "EvoDynamics.jl",
authors = "Ali R. Vahdati and Carlos Melian.",
doctest = false,
format = Documenter.HTML(
    prettyurls = get(ENV, "CI", nothing) == "true",
    ),
pages = [
	"Introduction" => "index.md",
	"Model description" => "model_description.md",
	"Tutorial" => "tutorial.md",
	"Examples" => [
	  "Simple Wright-Fisher" => "example1.md",
	  "Predator prey" => "example2.md"
		],
    ],
)

if get(ENV, "CI", nothing) == "true"
    deploydocs(repo = "github.com/kavir1698/EvoDynamics.jl.git",
               target = "build")
end


text:
# This file is a part of LegendTextIO.jl, licensed under the MIT License (MIT).

module LegendTextIO

using DelimitedFiles

using ArraysOfArrays
using BufferedStreams
using CSV
using LegendDataTypes
using RadiationDetectorSignals
using StaticArrays
using Unitful

using RadiationDetectorSignals: group_by_evtno

include("util.jl")
include("geant4_csv.jl")

## .root.hits files

import Base, Tables

using Mmap: mmap
using Parsers

export DarioHitsFile

include("dario_hits.jl")

end # module


text:
@generated function precomputeB(win, x, n::NTuple{D,Int}, m, M, sigma, T, U1, U2) where D
    quote
        I = zeros(Int64, (2*m+1)^D, M)
        J = zeros(Int64, M+1)
        V = zeros(T, (2*m+1)^D, M)

        J[1] = 1
        @inbounds @simd for k in 1:M
            @nexprs $D d -> xscale_d = x[d,k] * n[d]
            @nexprs $D d -> c_d = floor(Int, xscale_d)

            @nloops $D l d -> (c_d-m):(c_d+m) d->begin
                # preexpr
                gidx_d = rem(l_d+n[d], n[d]) + 1
                Iidx_d = l_d - c_d + m + 1
                idx = abs( (xscale_d - l_d) )
                tmpWin_d =  win(idx / n[d], n[d], m, sigma)

            end begin
                # bodyexpr
                v = 1
                @nexprs $D d -> v *= tmpWin_d
                i1 = 1
                @nexprs $D d -> i1 += (Iidx_d-1) * U1[d]
                i2 = 1
                @nexprs $D d -> i2 += (gidx_d-1) * U2[d]

                I[i1,k] = i2
                V[i1,k] = v
            end
            J[k+1] = J[k] + (2*m+1)^D
        end

        S = SparseMatrixCSC(prod(n), M, J, vec(I), vec(V))
        return S
    end
end