text:
mutable struct LabelIntervalControl <: AbstractControl
    isenabled::Bool
    nested_interval::NestedInterval
    buffer::String
end

function get_start(control::LabelIntervalControl)
    get_start(control.nested_interval)
end

function set_start!(control::LabelIntervalControl, val::Number)
    set_start!(control.nested_interval, val)
end

function get_stop(control::LabelIntervalControl)
    get_stop(control.nested_interval)
end

function set_stop!(control::LabelIntervalControl, val::Number)
    set_stop!(control.nested_interval, val)
end

function (control::LabelIntervalControl)(model::LabelledIntervals, properties::LabelIntervalDisplayProperties, plot_context::AbstractContext)
    id = get_id(properties)
    caption = get_caption(properties)
    col = get_col(properties)
    textcol = get_textcol(properties)
    rectangle = get_layout(properties)
    totalwidth = get_width(rectangle)
    totalheight = get_height(rectangle)
    padding = get_padding(properties)
    width = totalwidth - padding[2]
    height = totalheight

    # Plot the overview
    pos = CImGui.GetCursorScreenPos()
    plot_context()
    pos = ImVec2(pos.x + padding[2], pos.y + padding[1] + 1)

    draw_list = CImGui.GetWindowDrawList()
    yoffset = 3
    # Draw the labelled intervals
    for (key, labelled_interval) in pairs(get_labelled_intervals(model))
        # Map interval position to screen coordinates
        nested_interval = get_nested_interval(labelled_interval)
        interval = get_interval(nested_interval)
        start = get_start(nested_interval)
        stop =  get_stop(nested_interval)
        x₀ = stretch_linearly(start, first(interval), last(interval), pos.x, pos.x + width)
        x₁ = stretch_linearly(stop, first(interval), last(interval), pos.x, pos.x + width)
        # Shaded region
        CImGui.AddRectFilled(draw_list, ImVec2(x₀, pos.y + yoffset), ImVec2(x₁,  pos.y + height + yoffset), Base.convert(ImU32, col))
        # Text Label
        first_nul = findfirst(isequal('\0'), key)
        str = isnothing(first_nul) ? key : key[1:first_nul-1]
        CImGui.AddText(draw_list, ImVec2(x₀ + 5, pos.y + yoffset + 5), Base.convert(ImU32, textcol) , str)
        # Handle event selection
        CImGui.SetCursorScreenPos(x₀, pos.y)
        CImGui.InvisibleButton("###interval button $str", ImVec2(x₁ - x₀, height))
        # TODO Improve tooltip
        #CImGui.IsItemHovered() && CImGui.SetTooltip(string("Hello"))
        if CImGui.IsItemClicked()
            # Update the principal nested interval in accordance with the selection.
            set_start!(control.nested_interval, start)
            set_stop!(control.nested_interval, stop)
            # Update the text for the label in accordance with the selection.
            buffer = key*"\0"^(1)
            # Allow up to 64 characters for the label buffer.
            pad_null = max(0, 255 - length(buffer) + 1)
            control.buffer = buffer*"\0"^(pad_null)
        end
        CImGui.SameLine()
        CImGui.SetCursorScreenPos(pos.x, pos.y)
    end

    # TODO investigate IndentSpacing
    !isempty(get_labelled_intervals(model)) ? CImGui.SetCursorScreenPos(pos.x, pos.y + height + yoffset + 4) : nothing

    # Draw indicators for the currently selected interval
    current_nested_interval = control.nested_interval
    interval = get_interval(current_nested_interval)
    start = get_start(current_nested_interval)
    stop =  get_stop(current_nested_interval)
    x₀ = stretch_linearly(start, first(interval), last(interval), pos.x, pos.x + width)
    x₁ = stretch_linearly(stop, first(interval), last(interval), pos.x, pos.x + width)

    CImGui.AddLine(draw_list, ImVec2(x₀, pos.y + yoffset), ImVec2(x₀ , pos.y + height + yoffset), Base.convert(ImU32, ImVec4(0.0, 0.0, 0.0, 0.9)), Cfloat(2));
    CImGui.AddLine(draw_list, ImVec2(x₁, pos.y + yoffset), ImVec2(x₁ , pos.y + height + yoffset), Base.convert(ImU32, ImVec4(0.0, 0.0, 0.0, 0.9)), Cfloat(2));

    # Handle creation of new event labels
    CImGui.Indent(padding[2])
    CImGui.PushItemWidth(150)
    CImGui.InputText("###event description", control.buffer, length(control.buffer))
    CImGui.PopItemWidth()
    CImGui.SameLine()
    CImGui.Button("Create Interval Label") && add_label!(get_labelled_intervals(model), control.buffer, control.nested_interval)
    CImGui.SameLine()
    CImGui.Button("Remove Interval Label") && remove_label!(get_labelled_intervals(model), control.buffer)
    CImGui.Unindent(padding[2])
end

function add_label!(labelled_intervals::Dict{String, LabelledInterval}, label::String, ni::NestedInterval)
    key = extract_string(label)
    labelled_intervals[key] = LabelledInterval(key, copy(ni))
end

function remove_label!(labelled_intervals::Dict{String, LabelledInterval}, label::String)
    key = extract_string(label)
    if haskey(labelled_intervals, key)
        delete!(labelled_intervals, key)
    end
end


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

using PyPlot, PyCall
clf()

#plots = 1:16
plots = 1:9

for l in plots
    n = rand(1:500)
    y = sort(rand(n), rev=true)
    #plt = semilogy(y, color=ColorMap("tab20c")(l-1))
    col = mycolor(l)
    plt = semilogy(y, color=col)
end
xlabel("function calls")
ylabel("cost")

title("My nice looking figure")
legend(("plot $l" for l in plots))

using Match
function mycolor(l)
    cmap = ColorMap("tab20c")
    @match l begin
        1 => cmap(0)
        2 => cmap(1)
        3 => cmap(2)
        4 => cmap(4)
        5 => cmap(5)
        6 => cmap(6)
        7 => cmap(8)
        8 => cmap(9)
        9 => cmap(10)
        _ => cmap(19)
    end
end

display(plt)
savefig("test.eps")
savefig("test.pdf")
savefig("test.png")
savefig("test.svg")

ext: jl

max_stars_repo_name: briochemc/JuliaPlotting

max_issues_repo_name: briochemc/JuliaPlotting

max_forks_repo_name: briochemc/JuliaPlotting


text:
###
### Mer specific specializations of src/biosequence/indexing.jl
###

@inline _offset(K, i) = 2(K - i)
@inline offset(::Type{T}, i::Integer) where {A,K,T<:AbstractMer{A,K}} = 2(K - i)
@inline offset(x::AbstractMer, i::Integer) = offset(typeof(x), i)

@inline function inbounds_getindex(x::AbstractMer, i::Integer)
    return reinterpret(eltype(x), 0x01 << ((encoded_data(x) >> offset(x, i)) & 0b11))
end


text:
module REPLTreeViews

using TreeViews: numberofnodes, treelabel, treenode, hastreeview, nodelabel

struct REPLTreeViewsDisplay <: AbstractDisplay end

function Base.display(::REPLTreeViewsDisplay, x)
    if hastreeview(x)
        println()
        printMenu(stdout, generateTreeView(x), 999)
        println(stdout)
    else
        # fall through to the Displays lower in the display stack
        throw(MethodError(display, "nope"))
    end
end

terminal = nothing  # The user terminal
function __init__()
    global terminal
    terminal = REPL.Terminals.TTYTerminal(get(ENV, "TERM", Sys.iswindows() ? "" : "dumb"), stdin, stdout, stderr)
    pushdisplay(REPLTreeViewsDisplay())
end

include("renderer.jl")
include("default.jl")

end # module


text:
__precompile__()

module PopGenCore

using CSV, DataFrames, PooledArrays, StaticArrays
using CodecZlib, VariantCallFormat
using StatsBase: countmap, proportionmap
using Random: shuffle, shuffle!

include("PopData.jl")
include("PopDataWrappers.jl")
## Utilities
include("Utils/GeneralUtils.jl")
include("Utils/GenotypeUtils.jl")
include("Utils/ioUtils.jl")
include("Utils/MathUtils.jl")
include("Utils/MissingUtils.jl")
##
include("AlleleFreq.jl")
include("GenoFreq.jl")
include("Conditionals.jl")
include("Permutations.jl")
include("Iterators.jl")
include("Manipulate.jl")
## File IO
include("io/Delimited.jl")
include("io/Genepop.jl")
include("io/Structure.jl")
include("io/VariantCall.jl")
include("io/ReadWrite.jl")
##
include("Datasets.jl")

# precompile some file IO
#include("precompile/precompile.jl") ;

end


text:
using LinearAlgebra, Statistics
using Parameters, QuantEcon, DataFrames, Plots, Random

function ArellanoEconomy(;β = .953,
                          γ = 2.,
                          r = 0.017,
                          ρ = 0.945,
                          η = 0.025,
                          θ = 0.282,
                          ny = 21,
                          nB = 251)

    # create grids
    Bgrid = collect(range(-.4, .4, length = nB))
    mc = tauchen(ny, ρ, η)
    Π = mc.p
    ygrid = exp.(mc.state_values)
    ydefgrid = min.(.969 * mean(ygrid), ygrid)

    # define value functions
    # notice ordered different than Python to take
    # advantage of column major layout of Julia)
    vf = zeros(nB, ny)
    vd = zeros(1, ny)
    vc = zeros(nB, ny)
    policy = zeros(nB, ny)
    q = ones(nB, ny) .* (1 / (1 + r))
    defprob = zeros(nB, ny)

    return (β = β, γ = γ, r = r, ρ = ρ, η = η, θ = θ, ny = ny,
            nB = nB, ygrid = ygrid, ydefgrid = ydefgrid,
            Bgrid = Bgrid, Π = Π, vf = vf, vd = vd, vc = vc,
            policy = policy, q = q, defprob = defprob)
end

u(ae, c) = c^(1 - ae.γ) / (1 - ae.γ)

function one_step_update!(ae,
                          EV,
                          EVd,
                          EVc)

    # unpack stuff
    @unpack β, γ, r, ρ, η, θ, ny, nB = ae
    @unpack ygrid, ydefgrid, Bgrid, Π, vf, vd, vc, policy, q, defprob = ae
    zero_ind = searchsortedfirst(Bgrid, 0.)

    for iy in 1:ny
        y = ae.ygrid[iy]
        ydef = ae.ydefgrid[iy]

        # value of being in default with income y
        defval = u(ae, ydef) + β * (θ * EVc[zero_ind, iy] + (1-θ) * EVd[1, iy])
        ae.vd[1, iy] = defval

        for ib in 1:nB
            B = ae.Bgrid[ib]

            current_max = -1e14
            pol_ind = 0
            for ib_next=1:nB
                c = max(y - ae.q[ib_next, iy]*Bgrid[ib_next] + B, 1e-14)
                m = u(ae, c) + β * EV[ib_next, iy]

                if m > current_max
                current_max = m
                pol_ind = ib_next
                end

            end

            # update value and policy functions
            ae.vc[ib, iy] = current_max
            ae.policy[ib, iy] = pol_ind
            ae.vf[ib, iy] = defval > current_max ? defval : current_max
        end
    end
end

function compute_prices!(ae)
    # unpack parameters
    @unpack β, γ, r, ρ, η, θ, ny, nB = ae

    # create default values with a matching size
    vd_compat = repeat(ae.vd, nB)
    default_states = vd_compat .> ae.vc

    # update default probabilities and prices
    copyto!(ae.defprob, default_states * ae.Π')
    copyto!(ae.q, (1 .- ae.defprob) / (1 + r))
    return
end

function vfi!(ae; tol = 1e-8, maxit = 10000)

    # unpack stuff
    @unpack β, γ, r, ρ, η, θ, ny, nB = ae
    @unpack ygrid, ydefgrid, Bgrid, Π, vf, vd, vc, policy, q, defprob = ae
    Πt = Π'

    # Iteration stuff
    it = 0
    dist = 10.

    # allocate memory for update
    V_upd = similar(ae.vf)

    while dist > tol && it < maxit
        it += 1

        # compute expectations for this iterations
        # (we need Π' because of order value function dimensions)
        copyto!(V_upd, ae.vf)
        EV = ae.vf * Πt
        EVd = ae.vd * Πt
        EVc = ae.vc * Πt

        # update value function
        one_step_update!(ae, EV, EVd, EVc)

        # update prices
        compute_prices!(ae)

        dist = maximum(abs(x - y) for (x, y) in zip(V_upd, ae.vf))

        if it % 25 == 0
            println("Finished iteration $(it) with dist of $(dist)")
        end
    end
end

function QuantEcon.simulate(ae,
                            capT = 5000;
                            y_init = mean(ae.ygrid),
                            B_init = mean(ae.Bgrid),
                            )

    # get initial indices
    zero_index = searchsortedfirst(ae.Bgrid, 0.)
    y_init_ind = searchsortedfirst(ae.ygrid, y_init)
    B_init_ind = searchsortedfirst(ae.Bgrid, B_init)

    # create a QE MarkovChain
    mc = MarkovChain(ae.Π)
    y_sim_indices = simulate(mc, capT + 1; init = y_init_ind)

    # allocate and fill output
    y_sim_val = zeros(capT+1)
    B_sim_val, q_sim_val = similar(y_sim_val), similar(y_sim_val)
    B_sim_indices = fill(0, capT + 1)
    default_status = fill(false, capT + 1)
    B_sim_indices[1], default_status[1] = B_init_ind, false
    y_sim_val[1], B_sim_val[1] = ae.ygrid[y_init_ind], ae.Bgrid[B_init_ind]

    for t in 1:capT
        # get today's indexes
        yi, Bi = y_sim_indices[t], B_sim_indices[t]
        defstat = default_status[t]

        # if you are not in default
        if !defstat
            default_today = ae.vc[Bi, yi] < ae.vd[yi]

            if default_today
                # default values
                default_status[t] = true
                default_status[t + 1] = true
                y_sim_val[t] = ae.ydefgrid[y_sim_indices[t]]
                B_sim_indices[t + 1] = zero_index
                B_sim_val[t+1] = 0.
                q_sim_val[t] = ae.q[zero_index, y_sim_indices[t]]
                else
                default_status[t] = false
                y_sim_val[t] = ae.ygrid[y_sim_indices[t]]
                B_sim_indices[t + 1] = ae.policy[Bi, yi]
                B_sim_val[t + 1] = ae.Bgrid[B_sim_indices[t + 1]]
                q_sim_val[t] = ae.q[B_sim_indices[t + 1], y_sim_indices[t]]
            end

        # if you are in default
        else
        B_sim_indices[t + 1] = zero_index
        B_sim_val[t+1] = 0.
        y_sim_val[t] = ae.ydefgrid[y_sim_indices[t]]
        q_sim_val[t] = ae.q[zero_index, y_sim_indices[t]]

        # with probability θ exit default status
        default_status[t + 1] = rand() ≥ ae.θ
        end
    end

    return (y_sim_val[1:capT], B_sim_val[1:capT], q_sim_val[1:capT],
            default_status[1:capT])
end

## Solutions
using DataFrames, Plots

ae = ArellanoEconomy(β = 0.9832841691248771,     # time discount rate
                     γ = 2.,       # risk aversion
                     r = 0.017,    # international interest rate
                     ρ = .945,     # persistence in output
                     η = 0.025,    # st dev of output shock
                     θ = 0.282,    # prob of regaining access
                     ny = 21,      # number of points in y grid
                     nB = 251)     # number of points in B grid

# now solve the model on the grid.
vfi!(ae)

# create "Y High" and "Y Low" values as 5% devs from mean
high, low = 1.05 * mean(ae.ygrid), 0.95 * mean(ae.ygrid)
iy_high, iy_low = map(x -> searchsortedfirst(ae.ygrid, x), (high, low))

# extract a suitable plot grid
x = zeros(0)
q_low = zeros(0)
q_high = zeros(0)
for i in 1:ae.nB
    b = ae.Bgrid[i]
    if -0.35 ≤ b ≤ 0  # to match fig 3 of Arellano
        push!(x, b)
        push!(q_low, ae.q[i, iy_low])
        push!(q_high, ae.q[i, iy_high])
    end
end

# generate plot
plot(x, q_low, label = "Low")
plot!(x, q_high, label = "High")
plot!(title = "Bond price schedule q(y, B')",
      xlabel = "B'", ylabel = "q", legend_title = "y", legend = :topleft)

plot(ae.Bgrid, ae.vf[:, iy_low], label = "Low")
plot!(ae.Bgrid, ae.vf[:, iy_high], label = "High")
plot!(xlabel = "B", ylabel = "V(y,B)", title = "Value functions",
      legend_title="y", legend = :topleft)

heatmap(ae.Bgrid[1:end-1],
    ae.ygrid[2:end],
    reshape(clamp.(vec(ae.defprob[1:end - 1, 1:end - 1]), 0, 1), 250, 20)')
plot!(xlabel = "B'", ylabel = "y", title = "Probability of default",
      legend = :topleft)
savefig("figs/PofD_Ar.png")

using Random
# set random seed for consistent result
Random.seed!(348938)

# simulate
T = 10000
y_vec, B_vec, q_vec, default_vec = simulate(ae, T)

# find starting and ending periods of recessions
defs = findall(default_vec)
def_breaks = diff(defs) .> 1
def_start = defs[[true; def_breaks]]
def_end = defs[[def_breaks; true]]

y_vals = [y_vec, B_vec, q_vec]
titles = ["Output", "Foreign assets", "Bond price"]

plots = plot(layout = (3, 1), size = (700, 800))

# Plot the three variables, and for each each variable shading the period(s) of default
# in grey
for i in 1:3
    plot!(plots[i], 1:T, y_vals[i], title = titles[i], xlabel = "time", label = "", lw = 2)
    for j in 1:length(def_start)
        plot!(plots[i], [def_start[j], def_end[j]], fill(maximum(y_vals[i]), 2),
            fillrange = [extrema(y_vals[i])...], fcolor = :grey, falpha = 0.3, label = "")
    end
end

plot(plots)

#--------------------------------- Questions -----------------------------------------------#

# Question 1
histogram(y_vec, xlabel = "y", ylabel = "Frequency", label = "")
savefig("figs/endow_hist.png")
histogram(B_vec, xlabel = "B", ylabel = "Frequency", label = "")
savefig("figs/B_hist.png")
histogram(default_vec, xlabel = "Default = 0 No, = 1 Yes", ylabel = "Frequency", label = "")
savefig("figs/default_hist.png")

# Question 2
sum(B_vec)/length(B_vec)
minimum(B_vec)

# Question 3
sum(default_vec .== 1)

# Question 4
D_index = findall(==(1), default_vec)
notD_index = findall(==(0), default_vec)
y_D = y_vec[D_index]
y_notD = y_vec[notD_index]
avg_y_loss = mean(y_notD) - mean(y_D)

# Question 5
function AGEconomy(;β = .953,
    γ = 2.,
    r = 0.017,
    ρ = 0.945,
    η = 0.025,
    θ = 0.282,
    ny = 21,
    nB = 251)

    # create grids
    Bgrid = collect(range(-.4, .4, length = nB))
    mc = tauchen(ny, ρ, η)
    Π = mc.p
    ygrid = exp.(mc.state_values)
    ydefgrid = 0.98 * ygrid

    # define value functions
    # notice ordered different than Python to take
    # advantage of column major layout of Julia)
    vf = zeros(nB, ny)
    vd = zeros(1, ny)
    vc = zeros(nB, ny)
    policy = zeros(nB, ny)
    q = ones(nB, ny) .* (1 / (1 + r))
    defprob = zeros(nB, ny)

    return (β = β, γ = γ, r = r, ρ = ρ, η = η, θ = θ, ny = ny,
    nB = nB, ygrid = ygrid, ydefgrid = ydefgrid,
    Bgrid = Bgrid, Π = Π, vf = vf, vd = vd, vc = vc,
    policy = policy, q = q, defprob = defprob)
end

AGe = AGEconomy(β = .953,     # time discount rate
                     γ = 2.,       # risk aversion
                     r = 0.017,    # international interest rate
                     ρ = .945,     # persistence in output
                     η = 0.025,    # st dev of output shock
                     θ = 0.282,    # prob of regaining access
                     ny = 21,      # number of points in y grid
                     nB = 251)     # number of points in B grid

# now solve the model on the grid.
vfi!(AGe)

# create "Y High" and "Y Low" values as 5% devs from mean
high, low = 1.05 * mean(AGe.ygrid), 0.95 * mean(AGe.ygrid)
iy_high, iy_low = map(x -> searchsortedfirst(AGe.ygrid, x), (high, low))

# extract a suitable plot grid
x = zeros(0)
q_low = zeros(0)
q_high = zeros(0)
for i in 1:AGe.nB
    b = AGe.Bgrid[i]
    if -0.35 ≤ b ≤ 0  # to match fig 3 of Arellano
        push!(x, b)
        push!(q_low, AGe.q[i, iy_low])
        push!(q_high, AGe.q[i, iy_high])
    end
end

# generate plot
plot(x, q_low, label = "Low")
plot!(x, q_high, label = "High")
plot!(title = "Bond price schedule q(y, B')",
      xlabel = "B'", ylabel = "q", legend_title = "y", legend = :topleft)

plot(AGe.Bgrid, ae.vf[:, iy_low], label = "Low")
plot!(AGe.Bgrid, ae.vf[:, iy_high], label = "High")
plot!(xlabel = "B", ylabel = "V(y,B)", title = "Value functions",
      legend_title="y", legend = :topleft)

heatmap(AGe.Bgrid[1:end-1],
AGe.ygrid[2:end],
    reshape(clamp.(vec(AGe.defprob[1:end - 1, 1:end - 1]), 0, 1), 250, 20)')
plot!(xlabel = "B'", ylabel = "y", title = "Probability of default",
      legend = :topleft)

savefig("figs/PofD_AG.png")

using Random
# set random seed for consistent result
Random.seed!(348938)

# simulate
T = 10000
y_vec, B_vec, q_vec, default_vec = simulate(AGe, T)

# find starting and ending periods of recessions
defs = findall(default_vec)
def_breaks = diff(defs) .> 1
def_start = defs[[true; def_breaks]]
def_end = defs[[def_breaks; true]]

y_vals = [y_vec, B_vec, q_vec]
titles = ["Output", "Foreign assets", "Bond price"]

plots = plot(layout = (3, 1), size = (700, 800))

# Plot the three variables, and for each each variable shading the period(s) of default
# in grey
for i in 1:3
    plot!(plots[i], 1:T, y_vals[i], title = titles[i], xlabel = "time", label = "", lw = 2)
    for j in 1:length(def_start)
        plot!(plots[i], [def_start[j], def_end[j]], fill(maximum(y_vals[i]), 2),
            fillrange = [extrema(y_vals[i])...], fcolor = :grey, falpha = 0.3, label = "")
    end
end

plot(plots)

histogram(y_vec, xlabel = "y", ylabel = "Frequency", label = "")
savefig("figs/endow_histQ5.png")
histogram(B_vec, xlabel = "B", ylabel = "Frequency", label = "")
savefig("figs/B_histQ5.png")
histogram(default_vec, yaxis=:log, xlabel = "Default = 0 No, = 1 Yes", ylabel = "Log Frequency", label = "")
savefig("figs/default_histQ5.png")

sum(B_vec)/length(B_vec)
minimum(B_vec)

sum(default_vec .== 1)

D_index = findall(==(1), default_vec)
notD_index = findall(==(0), default_vec)
y_D = y_vec[D_index]
y_notD = y_vec[notD_index]
avg_y_loss = mean(y_notD) - mean(y_D)

# Question 7
histogram(B_vec, xlabel = "B", ylabel = "Frequency", label = "")
savefig("figs/B_histQ7.png")

mean(B_vec)
minimum(B_vec)
maximum(B_vec)

# Question 8
histogram(B_vec, xlabel = "B", ylabel = "Frequency", label = "")
savefig("figs/B_histQ8.png")

mean(B_vec)
minimum(B_vec)
maximum(B_vec)


text:
#=
This file is auto-generated. Do not edit.
=#
"""
    mutable struct AVRTypeII <: AVR
        K0::Float64
        T1::Float64
        T2::Float64
        T3::Float64
        T4::Float64
        Te::Float64
        Tr::Float64
        Vr_max::Float64
        Vr_min::Float64
        Ae::Float64
        Be::Float64
        V_ref::Float64
        ext::Dict{String, Any}
        states::Vector{Symbol}
        n_states::Int64
        internal::InfrastructureSystemsInternal
    end

Parameters of an Automatic Voltage Regulator Type II -  Typical static exciter model

# Arguments
- `K0::Float64`: Regulator Gain, validation range: (0, nothing)
- `T1::Float64`: First Pole in s, validation range: (0, nothing)
- `T2::Float64`: First zero in s, validation range: (0, nothing)
- `T3::Float64`: First Pole in s, validation range: (0, nothing)
- `T4::Float64`: First zero in s, validation range: (0, nothing)
- `Te::Float64`: Field Circuit Time Constant in s, validation range: (0, nothing)
- `Tr::Float64`: Voltage Measurement Time Constant in s, validation range: (0, nothing)
- `Vr_max::Float64`: Maximum regulator voltage in pu, validation range: (0, nothing)
- `Vr_min::Float64`: Minimum regulator voltage in pu, validation range: (0, nothing)
- `Ae::Float64`: 1st ceiling coefficient, validation range: (0, nothing)
- `Be::Float64`: 2nd ceiling coefficient, validation range: (0, nothing)
- `V_ref::Float64`: Reference Voltage Set-point, validation range: (0, nothing)
- `ext::Dict{String, Any}`
- `states::Vector{Symbol}`
- `n_states::Int64`
- `internal::InfrastructureSystemsInternal`: power system internal reference, do not modify
"""
mutable struct AVRTypeII <: AVR
    "Regulator Gain"
    K0::Float64
    "First Pole in s"
    T1::Float64
    "First zero in s"
    T2::Float64
    "First Pole in s"
    T3::Float64
    "First zero in s"
    T4::Float64
    "Field Circuit Time Constant in s"
    Te::Float64
    "Voltage Measurement Time Constant in s"
    Tr::Float64
    "Maximum regulator voltage in pu"
    Vr_max::Float64
    "Minimum regulator voltage in pu"
    Vr_min::Float64
    "1st ceiling coefficient"
    Ae::Float64
    "2nd ceiling coefficient"
    Be::Float64
    "Reference Voltage Set-point"
    V_ref::Float64
    ext::Dict{String, Any}
    states::Vector{Symbol}
    n_states::Int64
    "power system internal reference, do not modify"
    internal::InfrastructureSystemsInternal
end

function AVRTypeII(K0, T1, T2, T3, T4, Te, Tr, Vr_max, Vr_min, Ae, Be, V_ref=1.0, ext=Dict{String, Any}(), )
    AVRTypeII(K0, T1, T2, T3, T4, Te, Tr, Vr_max, Vr_min, Ae, Be, V_ref, ext, [:Vf, :Vr1, :Vr2, :Vm], 4, InfrastructureSystemsInternal(), )
end

function AVRTypeII(; K0, T1, T2, T3, T4, Te, Tr, Vr_max, Vr_min, Ae, Be, V_ref=1.0, ext=Dict{String, Any}(), )
    AVRTypeII(K0, T1, T2, T3, T4, Te, Tr, Vr_max, Vr_min, Ae, Be, V_ref, ext, )
end

# Constructor for demo purposes; non-functional.
function AVRTypeII(::Nothing)
    AVRTypeII(;
        K0=0,
        T1=0,
        T2=0,
        T3=0,
        T4=0,
        Te=0,
        Tr=0,
        Vr_max=0,
        Vr_min=0,
        Ae=0,
        Be=0,
        V_ref=0,
        ext=Dict{String, Any}(),
    )
end

"""Get AVRTypeII K0."""
get_K0(value::AVRTypeII) = value.K0
"""Get AVRTypeII T1."""
get_T1(value::AVRTypeII) = value.T1
"""Get AVRTypeII T2."""
get_T2(value::AVRTypeII) = value.T2
"""Get AVRTypeII T3."""
get_T3(value::AVRTypeII) = value.T3
"""Get AVRTypeII T4."""
get_T4(value::AVRTypeII) = value.T4
"""Get AVRTypeII Te."""
get_Te(value::AVRTypeII) = value.Te
"""Get AVRTypeII Tr."""
get_Tr(value::AVRTypeII) = value.Tr
"""Get AVRTypeII Vr_max."""
get_Vr_max(value::AVRTypeII) = value.Vr_max
"""Get AVRTypeII Vr_min."""
get_Vr_min(value::AVRTypeII) = value.Vr_min
"""Get AVRTypeII Ae."""
get_Ae(value::AVRTypeII) = value.Ae
"""Get AVRTypeII Be."""
get_Be(value::AVRTypeII) = value.Be
"""Get AVRTypeII V_ref."""
get_V_ref(value::AVRTypeII) = value.V_ref
"""Get AVRTypeII ext."""
get_ext(value::AVRTypeII) = value.ext
"""Get AVRTypeII states."""
get_states(value::AVRTypeII) = value.states
"""Get AVRTypeII n_states."""
get_n_states(value::AVRTypeII) = value.n_states
"""Get AVRTypeII internal."""
get_internal(value::AVRTypeII) = value.internal

"""Set AVRTypeII K0."""
set_K0!(value::AVRTypeII, val::Float64) = value.K0 = val
"""Set AVRTypeII T1."""
set_T1!(value::AVRTypeII, val::Float64) = value.T1 = val
"""Set AVRTypeII T2."""
set_T2!(value::AVRTypeII, val::Float64) = value.T2 = val
"""Set AVRTypeII T3."""
set_T3!(value::AVRTypeII, val::Float64) = value.T3 = val
"""Set AVRTypeII T4."""
set_T4!(value::AVRTypeII, val::Float64) = value.T4 = val
"""Set AVRTypeII Te."""
set_Te!(value::AVRTypeII, val::Float64) = value.Te = val
"""Set AVRTypeII Tr."""
set_Tr!(value::AVRTypeII, val::Float64) = value.Tr = val
"""Set AVRTypeII Vr_max."""
set_Vr_max!(value::AVRTypeII, val::Float64) = value.Vr_max = val
"""Set AVRTypeII Vr_min."""
set_Vr_min!(value::AVRTypeII, val::Float64) = value.Vr_min = val
"""Set AVRTypeII Ae."""
set_Ae!(value::AVRTypeII, val::Float64) = value.Ae = val
"""Set AVRTypeII Be."""
set_Be!(value::AVRTypeII, val::Float64) = value.Be = val
"""Set AVRTypeII V_ref."""
set_V_ref!(value::AVRTypeII, val::Float64) = value.V_ref = val
"""Set AVRTypeII ext."""
set_ext!(value::AVRTypeII, val::Dict{String, Any}) = value.ext = val
"""Set AVRTypeII states."""
set_states!(value::AVRTypeII, val::Vector{Symbol}) = value.states = val
"""Set AVRTypeII n_states."""
set_n_states!(value::AVRTypeII, val::Int64) = value.n_states = val
"""Set AVRTypeII internal."""
set_internal!(value::AVRTypeII, val::InfrastructureSystemsInternal) = value.internal = val


text:
mikolovptburl = "http://www.fit.vutbr.cz/~imikolov/rnnlm/simple-examples.tgz"
mikolovptbdir = joinpath(@__DIR__,"mikolovptb")
mikolovptbtrn = "ptb.train.txt"
mikolovptbdev = "ptb.valid.txt"
mikolovptbtst = "ptb.test.txt"

"""

    mikolovptb()

Read [PTB](https://catalog.ldc.upenn.edu/ldc99t42) text from Mikolov's
[RNNLM](http://www.fit.vutbr.cz/~imikolov/rnnlm) toolkit which has
been lowercased and reduced to a 10K vocabulary size.  Return a tuple
(trn,dev,tst,vocab) where

    trn::Vector{Vector{UInt16}}: 42068 sentences, 887521 words
    dev::Vector{Vector{UInt16}}: 3370 sentences, 70390 words
    tst::Vector{Vector{UInt16}}: 3761 sentences, 78669 words
    vocab::Vector{String}: 9999 unique words

"""
function mikolovptb()
    global _mptb_trn, _mptb_dev, _mptb_tst, _mptb_vocab
    if !(@isdefined _mptb_trn)
        isdir(mikolovptbdir) || mkpath(mikolovptbdir)
        if !isfile(joinpath(mikolovptbdir, "ptb.train.txt"))
            @info("Downloading $mikolovptburl")
            tgz = download(mikolovptburl)
            run(`tar --strip-components 3 -C $mikolovptbdir -xzf $tgz ./simple-examples/data/ptb.train.txt ./simple-examples/data/ptb.valid.txt ./simple-examples/data/ptb.test.txt`)
        end
        dict = Dict{String,Int}()
        data = Vector{Vector{UInt16}}[]
        for file in (mikolovptbtrn, mikolovptbdev, mikolovptbtst)
            sentences = Vector{UInt16}[]
            for line in eachline(joinpath(mikolovptbdir,file))
                words = UInt16[]
                for word in split(line)
                    widx = get!(dict, word, 1+length(dict))
                    push!(words, widx)
                end
                push!(sentences, words)
            end
            push!(data, sentences)
        end
        _mptb_trn, _mptb_dev, _mptb_tst = data
        _mptb_vocab = Array{String}(undef, length(dict))
        for (k,v) in dict; _mptb_vocab[v] = k; end
    end
    return _mptb_trn, _mptb_dev, _mptb_tst, _mptb_vocab
end


text:
using FactCheck

inc(x) = x + 1

facts("Succeeding examples") do

    @fact 1 => 1

    context("You can define contexts") do
        @fact 1 => 1
    end

    @fact_throws error("neat")

    context("group facts") do
        @fact 1 => not(2)
        @fact 2 => not(isodd)
    end

    context("generating facts") do
        for i=1:5
            @fact i => i
        end
    end

end

facts("Failing examples") do

    @fact 1 => 2

    context("strings are strings") do
        @fact "bar" => "barr"
        @fact "baz" => "bazz"
    end

    context("some numbers are even") do
        @fact 3 => iseven
    end

    context() do
        x = 10
        @fact inc(inc(inc(0))) => 2
    end

    context() do
        x = 5
        y = 10
        @fact x => y
    end

    context("throws an error") do
        @fact error("foo") => 1
    end

end

facts("Assertion helpers") do

    @fact 1 => not(iseven)
    @fact 1 => truthy
    @fact nothing => falsey
    @fact false => anything
    @fact iseven => exactly(iseven)
    @fact 2.499999 => roughly(2.5)

end


text:
# This file is a part of TypeDBClient.  License is MIT: https://github.com/Humans-of-Julia/TypeDBClient.jl/blob/main/LICENSE

#
# package typedb.client.api.concept.type;
#
# import typedb.client.api.TypeDBTransaction;
# import typedb.client.api.concept.thing.Attribute;
# import typedb.client.common.exception.TypeDBClientException;
# import typedb.protocol.ConceptProto;
#
# import javax.annotation.CheckReturnValue;
# import javax.annotation.Nullable;
# import java.time.LocalDateTime;
# import java.util.stream.Stream;
#
# import static typedb.client.common.exception.ErrorMessage.Concept.BAD_VALUE_TYPE;
#
# public interface AttributeType extends ThingType {
#
#     @CheckReturnValue
#     default ValueType getValueType() {
#         return ValueType.OBJECT;
#     }
#
#     @CheckReturnValue
#     boolean isKeyable();
#
#     @Override
#     @CheckReturnValue
#     default boolean isAttributeType() {
#         return true;
#     }
#
#     @CheckReturnValue
#     default boolean isBoolean() {
#         return false;
#     }
#
#     @CheckReturnValue
#     default boolean isLong() {
#         return false;
#     }
#
#     @CheckReturnValue
#     default boolean isDouble() {
#         return false;
#     }
#
#     @CheckReturnValue
#     default boolean isString() {
#         return false;
#     }
#
#     @CheckReturnValue
#     default boolean isDateTime() {
#         return false;
#     }
#
#     @Override
#     @CheckReturnValue
#     AttributeType.Remote asRemote(TypeDBTransaction transaction);
#
#     @CheckReturnValue
#     AttributeType.Boolean asBoolean();
#
#     @CheckReturnValue
#     AttributeType.Long asLong();
#
#     @CheckReturnValue
#     AttributeType.Double asDouble();
#
#     @CheckReturnValue
#     AttributeType.String asString();
#
#     @CheckReturnValue
#     AttributeType.DateTime asDateTime();
#
#     enum ValueType {
#
#         OBJECT(Object.class, false, false),
#         BOOLEAN(Boolean.class, true, false),
#         LONG(Long.class, true, true),
#         DOUBLE(Double.class, true, false),
#         STRING(String.class, true, true),
#         DATETIME(LocalDateTime.class, true, true);
#
#         private final Class<?> valueClass;
#         private final boolean isWritable;
#         private final boolean isKeyable;
#
#         ValueType(Class<?> valueClass, boolean isWritable, boolean isKeyable) {
#             this.valueClass = valueClass;
#             this.isWritable = isWritable;
#             this.isKeyable = isKeyable;
#         }
#
#         @CheckReturnValue
#         public static ValueType of(Class<?> valueClass) {
#             for (ValueType t : ValueType.values()) {
#                 if (t.valueClass == valueClass) {
#                     return t;
#                 }
#             }
#             throw new TypeDBClientException(BAD_VALUE_TYPE);
#         }
#
#         @CheckReturnValue
#         public Class<?> valueClass() {
#             return valueClass;
#         }
#
#         @CheckReturnValue
#         public boolean isWritable() {
#             return isWritable;
#         }
#
#         @CheckReturnValue
#         public boolean isKeyable() {
#             return isKeyable;
#         }
#
#         @CheckReturnValue
#         public ConceptProto.AttributeType.ValueType proto() {
#             switch (this) {
#                 case OBJECT:
#                     return ConceptProto.AttributeType.ValueType.OBJECT;
#                 case BOOLEAN:
#                     return ConceptProto.AttributeType.ValueType.BOOLEAN;
#                 case LONG:
#                     return ConceptProto.AttributeType.ValueType.LONG;
#                 case DOUBLE:
#                     return ConceptProto.AttributeType.ValueType.DOUBLE;
#                 case STRING:
#                     return ConceptProto.AttributeType.ValueType.STRING;
#                 case DATETIME:
#                     return ConceptProto.AttributeType.ValueType.DATETIME;
#                 default:
#                     return ConceptProto.AttributeType.ValueType.UNRECOGNIZED;
#             }
#         }
#     }
#
#     interface Remote extends ThingType.Remote, AttributeType {
#
#         void setSupertype(AttributeType attributeType);
#
#         @Override
#         @CheckReturnValue
#         Stream<? extends AttributeType> getSubtypes();
#
#         @Override
#         @CheckReturnValue
#         Stream<? extends Attribute<?>> getInstances();
#
#         @CheckReturnValue
#         Stream<? extends ThingType> getOwners();
#
#         @CheckReturnValue
#         Stream<? extends ThingType> getOwners(boolean onlyKey);
#
#         @Override
#         @CheckReturnValue
#         AttributeType.Remote asAttributeType();
#
#         @Override
#         @CheckReturnValue
#         AttributeType.Boolean.Remote asBoolean();
#
#         @Override
#         @CheckReturnValue
#         AttributeType.Long.Remote asLong();
#
#         @Override
#         @CheckReturnValue
#         AttributeType.Double.Remote asDouble();
#
#         @Override
#         @CheckReturnValue
#         AttributeType.String.Remote asString();
#
#         @Override
#         @CheckReturnValue
#         AttributeType.DateTime.Remote asDateTime();
#     }
#
#     interface Boolean extends AttributeType {
#
#         @Override
#         @CheckReturnValue
#         default ValueType getValueType() {
#             return ValueType.BOOLEAN;
#         }
#
#         @Override
#         @CheckReturnValue
#         default boolean isBoolean() {
#             return true;
#         }
#
#         @Override
#         @CheckReturnValue
#         AttributeType.Boolean.Remote asRemote(TypeDBTransaction transaction);
#
#         interface Remote extends AttributeType.Boolean, AttributeType.Remote {
#
#             Attribute.Boolean put(boolean value);
#
#             @Nullable
#             @CheckReturnValue
#             Attribute.Boolean get(boolean value);
#
#             @Override
#             @CheckReturnValue
#             Stream<? extends Attribute.Boolean> getInstances();
#
#             @Override
#             @CheckReturnValue
#             Stream<? extends AttributeType.Boolean> getSubtypes();
#
#             void setSupertype(AttributeType.Boolean booleanAttributeType);
#         }
#     }
#
#     interface Long extends AttributeType {
#
#         @Override
#         @CheckReturnValue
#         default ValueType getValueType() {
#             return ValueType.LONG;
#         }
#
#         @Override
#         @CheckReturnValue
#         default boolean isLong() {
#             return true;
#         }
#
#         @Override
#         @CheckReturnValue
#         AttributeType.Long.Remote asRemote(TypeDBTransaction transaction);
#
#         interface Remote extends AttributeType.Long, AttributeType.Remote {
#
#             Attribute.Long put(long value);
#
#             @Nullable
#             @CheckReturnValue
#             Attribute.Long get(long value);
#
#             @Override
#             @CheckReturnValue
#             Stream<? extends Attribute.Long> getInstances();
#
#             @Override
#             @CheckReturnValue
#             Stream<? extends AttributeType.Long> getSubtypes();
#
#             void setSupertype(AttributeType.Long longAttributeType);
#         }
#     }
#
#     interface Double extends AttributeType {
#
#         @Override
#         @CheckReturnValue
#         default ValueType getValueType() {
#             return ValueType.DOUBLE;
#         }
#
#         @Override
#         @CheckReturnValue
#         default boolean isDouble() {
#             return true;
#         }
#
#         @Override
#         @CheckReturnValue
#         AttributeType.Double.Remote asRemote(TypeDBTransaction transaction);
#
#         interface Remote extends AttributeType.Double, AttributeType.Remote {
#
#             Attribute.Double put(double value);
#
#             @Nullable
#             @CheckReturnValue
#             Attribute.Double get(double value);
#
#             @Override
#             @CheckReturnValue
#             Stream<? extends Attribute.Double> getInstances();
#
#             @Override
#             @CheckReturnValue
#             Stream<? extends AttributeType.Double> getSubtypes();
#
#             void setSupertype(AttributeType.Double doubleAttributeType);
#         }
#     }
#
#     interface String extends AttributeType {
#
#         @Override
#         @CheckReturnValue
#         default ValueType getValueType() {
#             return ValueType.STRING;
#         }
#
#         @Override
#         @CheckReturnValue
#         default boolean isString() {
#             return true;
#         }
#
#         @Override
#         @CheckReturnValue
#         AttributeType.String.Remote asRemote(TypeDBTransaction transaction);
#
#         interface Remote extends AttributeType.String, AttributeType.Remote {
#
#             Attribute.String put(java.lang.String value);
#
#             @Nullable
#             @CheckReturnValue
#             Attribute.String get(java.lang.String value);
#
#             @Override
#             @CheckReturnValue
#             Stream<? extends Attribute.String> getInstances();
#
#             @Nullable
#             @CheckReturnValue
#             java.lang.String getRegex();
#
#             void setRegex(java.lang.String regex);
#
#             @Override
#             @CheckReturnValue
#             Stream<? extends AttributeType.String> getSubtypes();
#
#             void setSupertype(AttributeType.String stringAttributeType);
#         }
#     }
#
#     interface DateTime extends AttributeType {
#
#         @Override
#         @CheckReturnValue
#         default ValueType getValueType() {
#             return ValueType.DATETIME;
#         }
#
#         @Override
#         @CheckReturnValue
#         default boolean isDateTime() {
#             return true;
#         }
#
#         @Override
#         @CheckReturnValue
#         AttributeType.DateTime.Remote asRemote(TypeDBTransaction transaction);
#
#         interface Remote extends AttributeType.DateTime, AttributeType.Remote {
#
#             Attribute.DateTime put(LocalDateTime value);
#
#             @Nullable
#             @CheckReturnValue
#             Attribute.DateTime get(LocalDateTime value);
#
#             @Override
#             @CheckReturnValue
#             Stream<? extends Attribute.DateTime> getInstances();
#
#             @Override
#             @CheckReturnValue
#             Stream<? extends AttributeType.DateTime> getSubtypes();
#
#             void setSupertype(AttributeType.DateTime dateTimeAttributeType);
#         }
#     }
# }


text:
#=
Conway's Game of Life takes place on an infinite two-dimensional board of square cells. Each cell is either dead or alive, and at each tick, the following rules apply:

    Any live cell with less than two live neighbours dies.
    Any live cell with two or three live neighbours remains living.
    Any live cell with more than three live neighbours dies.
    Any dead cell with exactly three live neighbours becomes a live cell.

A cell neighbours another cell if it is horizontally, vertically, or diagonally adjacent.

Implement Conway's Game of Life. It should be able to be initialized with a starting list of live cell coordinates and the number of steps it should run for. Once initialized, it should print out the board state at each step. Since it's an infinite board, print out only the relevant coordinates, i.e. from the top-leftmost live cell to bottom-rightmost live cell.

You can represent a live cell with an asterisk (*) and a dead cell with a dot (.).
=#

max_stars_repo_name: DominiqueCaron/daily-coding-problem

max_issues_repo_name: DominiqueCaron/daily-coding-problem

max_forks_repo_name: DominiqueCaron/daily-coding-problem


text:
list = ["NAME: Rajavee Srivastava","EMAIL: rajvisrivastav18@gmail.com","USED LANGUAGE: JULIA", "BIOSTACK: Genomics","SLACK USERNAME: @Rajvi"]
 for i in list
       println(i)
       end

max_stars_repo_name: jravimicro/team-rosalind-project

max_issues_repo_name: jravimicro/team-rosalind-project

max_forks_repo_name: jravimicro/team-rosalind-project


text:
## Medical Diagnosis Benchmark
## ==========================

using Sigma
using DataFrames
disease_priors = [0.025851035,0.058152077,0.020436158,0.040025596,
                  0.032941836,0.078530717,0.005548213,0.024689777,
                  0.029336553,0.020641227,0.074125078,0.028544798,
                  0.054423213,0.077353640,0.078099579,0.065297382,
                  0.039513058,0.030245903,0.128596589,0.095504584]
diseases = RandArray(map(p->flip(p),disease_priors))

datapath = joinpath(benchdir, "diseasediagnosis","problem-2-edges.csv")
edges = readtable(datapath)
nfindings = size(edges,2)
nfindings = 10
ndiseases = size(edges,1)

findings = RandVar{Bool}[]
for i = 1:nfindings
  ands = RandVar{Bool}[]
  for j = 1:ndiseases
    if edges[j,i] != 0.0
      push!(ands, diseases[j] & flip(edges[j,i]))
    end
  end
  push!(findings, (|)(ands...))
end

findings_randvar = RandArray(findings)

## Benchmarking

samples = cond_sample_mh(diseases,findings_randvar == [randbool() for i = 1:nfindings],10)
@show samples


text:
module ergm

include("spaces.jl")
include("stats.jl")
include("models.jl")
include("sampler.jl")
include("optim.jl")
include("inference.jl")
include("platinum.jl")

end


text:
#=
  de Bruijn sequence in Julia ConstraintSolver

  Implementation of de Bruijn sequences in Comet, both "classical"
  and "arbitrary".

  Compare with the the web based programs:
    http://www.hakank.org/comb/debruijn.cgi
    http://www.hakank.org/comb/debruijn_arb.cgi

  For Base = 2, N = 3, M = 8 there are 2 solutions:
    x : [](0, 1, 3, 7, 6, 5, 2, 4)
    bincode : [0, 0, 0, 1, 1, 1, 0, 1]

    x : [](0, 1, 2, 5, 3, 7, 6, 4)
    bincode : [0, 0, 0, 1, 0, 1, 1, 1]

  There are two symmetry breaking constraints:
  - if possible, we require that the occurrences of number are the same,
    i.e if M mod Base == 0.
  - the sequence X always start with the minumum value (0).

  Model created by Hakan Kjellerstrand, hakank@gmail.com
  See also my Julia page: http://www.hakank.org/julia/

=#
using ConstraintSolver, JuMP
using Cbc, GLPK, Ipopt
const CS = ConstraintSolver
include("constraints_utils.jl")

function debruijn(base=2,n=3,m=base^n,all_solutions=true,print_solutions=true)
    println("base:$base n:$n m:$m")
    cbc_optimizer = optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 0)
    glpk_optimizer = optimizer_with_attributes(GLPK.Optimizer)
    ipopt_optimizer = optimizer_with_attributes(Ipopt.Optimizer)

    model = Model(optimizer_with_attributes(CS.Optimizer,   "all_solutions"=> all_solutions,
                                                            # "all_optimal_solutions"=>true,
                                                            "logging"=>[],

                                                            # "traverse_strategy"=>:BFS,
                                                            # "traverse_strategy"=>:DFS,
                                                            # "traverse_strategy"=>:DBFS,

                                                            # "branch_split"=>:Smallest,
                                                            # "branch_split"=>:Biggest,
                                                            # "branch_split"=>:InHalf,

                                                            # "simplify"=>false,
                                                            # "simplify"=>true, # default

                                                            "time_limit"=>30,

                                                            # "lp_optimizer" => cbc_optimizer,
                                                            # "lp_optimizer" => glpk_optimizer,
                                                            "lp_optimizer" => ipopt_optimizer,
                                        ))

    # x: list of the integers to use which are converted to
    #    'base-ary' numbers below
    @variable(model, 0 <= x[1:m] <= (base^n)-1, Int)

    # binary: the matrix of the "fully expanded" integers in X
    @variable(model, 0 <= binary[1:m,1:n] <= base-1, Int)

    # number of occurrences for each number
    @variable(model, 0 <= gcc[1:base] <= m, Int)

    @constraint(model, x in CS.AllDifferentSet())

    # Connect x and binary
    for i in 1:m
        # convert the integer to base-ary representation
        to_num(model, [binary[i,j] for j in 1:n], base, x[i])
    end

    # The de Bruijn criterion: Connect one element to the next...
    for i in 2:m, j in 2:n
        @constraint(model,binary[i-1, j] == binary[i, j-1])
    end

    # ... and around the corner.
    for j in 2:n
        @constraint(model, binary[m, j] == binary[1, j-1])
    end

    # bin_code: The de Bruijn sequence, i.e. the first
    #           elements in each row in Binary
    #           Note: this is just a slice.
    bin_code = [binary[i,1] for i in 1:m]

    # GCC: Count the number of different elements in the bin code
    # (poor man's global cardinarlity count...)
    for i in 1:base
        count_ctr(model,bin_code, :(==), i-1, gcc[i])
    end

    # If possible, we require that the occurrences of number are the same.
    # Gcc
    if m % base == 0
        for i in 2:base
            @constraint(model,gcc[i] == gcc[i-1])
            @constraint(model,gcc[i] == round(Int,m / base))
        end
    end

    # symmetry breaking
    # my_min is defined in constraints_utils.jl
    my_min(model, x, x[1])

    println("solve")
    # Solve the problem
    optimize!(model)

    status = JuMP.termination_status(model)
    # println("status:$status")
    if status == MOI.OPTIMAL
        num_sols = MOI.get(model, MOI.ResultCount())
        println("num_sols:$num_sols\n")
        if print_solutions
            for sol in 1:num_sols
                println("solution #$sol")
                xx = convert.(Integer,JuMP.value.(x; result=sol))
                binaryx = convert.(Integer,JuMP.value.(binary; result=sol))
                gccx = convert.(Integer,JuMP.value.(gcc; result=sol))
                if base < 10
                    bin_codex = join([binaryx[i,1] for i in 1:m],"")
                else
                    bin_codex = [binaryx[i,1] for i in 1:m]
                end

                println("x: $xx")
                println("binary:")
                for row in eachrow(binaryx)
                    println(row)
                end
                println("bin_code: $bin_codex")
                println("gcc: $gccx")
                println()

            end
            println("num_sols: $num_sols\n")
        end
    else
        println("status: $status")
    end

    return status
end

base=2
n=3
m=base^n
all_solutions= true
print_solutions=true
@time debruijn(base,n,m,all_solutions,print_solutions)

# Test a non standard version: m = 27 (instead of 2^5=32)
# It should be 2000 solutions. And it is, but it takes ≈ 29s!
base = 2
n = 5
m = 27
@time debruijn(base,n,m,all_solutions,false)

max_stars_repo_name: Wikunia/hakank

max_issues_repo_name: Wikunia/hakank

max_forks_repo_name: Wikunia/hakank


text:
const CBF = MathOptFormat.CBF

const CBF_TEST_FILE = "test.cbf"
const MODELS_DIR = joinpath(@__DIR__, "models")

function set_var_and_con_names(model::MOI.ModelLike)
    variable_names = String[]
    for j in MOI.get(model, MOI.ListOfVariableIndices())
        var_name_j = "v" * string(j.value)
        push!(variable_names, var_name_j)
        MOI.set(model, MOI.VariableName(), j, var_name_j)
    end

    idx = 0
    constraint_names = String[]
    for i in MOI.get(model, MOI.ListOfConstraintIndices{MOI.SingleVariable, MOI.Integer}())
        idx += 1
        con_name_i = "c" * string(idx)
        push!(constraint_names, con_name_i)
        MOI.set(model, MOI.ConstraintName(), i, con_name_i)
    end
    for S in [
        MOI.Reals, MOI.Zeros, MOI.Nonnegatives, MOI.Nonpositives,
        MOI.SecondOrderCone, MOI.RotatedSecondOrderCone,
        MOI.PositiveSemidefiniteConeTriangle,
        MOI.ExponentialCone, MOI.DualExponentialCone,
        MOI.PowerCone{Float64}, MOI.DualPowerCone{Float64},
        ],
        F in [MOI.VectorOfVariables, MOI.VectorAffineFunction{Float64}]
        for i in MOI.get(model, MOI.ListOfConstraintIndices{F, S}())
            idx += 1
            con_name_i = "c" * string(idx)
            push!(constraint_names, con_name_i)
            MOI.set(model, MOI.ConstraintName(), i, con_name_i)
        end
    end

    return (variable_names, constraint_names)
end

Base.isapprox(f::MOI.VectorOfVariables, g::MOI.VectorAffineFunction{Float64};
    kwargs...) = isapprox(MOI.VectorAffineFunction{Float64}(f), g)

function test_write_then_read(model_string::String)
    model1 = CBF.Model()
    MOIU.loadfromstring!(model1, model_string)
    (variable_names, constraint_names) = set_var_and_con_names(model1)

    MOI.write_to_file(model1, CBF_TEST_FILE)
    model2 = CBF.Model()
    MOI.read_from_file(model2, CBF_TEST_FILE)
    set_var_and_con_names(model2)

    MOIU.test_models_equal(model1, model2, variable_names, constraint_names)
end

function test_read(filename::String, model_string::String)
    model1 = CBF.Model()
    MOIU.loadfromstring!(model1, model_string)
    (variable_names, constraint_names) = set_var_and_con_names(model1)

    model2 = CBF.Model()
    MOI.read_from_file(model2, filename)
    set_var_and_con_names(model2)

    MOIU.test_models_equal(model1, model2, variable_names, constraint_names)
end

@test sprint(show, CBF.Model()) == "A Conic Benchmark Format (CBF) model"

@testset "Support errors" begin
    @testset "$set variable bound" for set in [
            MOI.EqualTo(1.0),
            MOI.LessThan(1.0),
            MOI.GreaterThan(1.0),
            MOI.Interval(1.0, 2.0),
            MOI.Semiinteger(1.0, 2.0),
            MOI.Semicontinuous(1.0, 2.0),
            MOI.ZeroOne()
        ]
        model_string = """
        variables: x
        minobjective: x
        c: x in $set
        """
        model = CBF.Model()
        err = MOI.UnsupportedConstraint{MOI.SingleVariable, typeof(set)}
        @test_throws err MOIU.loadfromstring!(model, model_string)
    end
end

@testset "Read errors" begin
    @testset "Non-empty model" begin
        model = CBF.Model()
        MOI.add_variable(model)
        @test_throws Exception MOI.read_from_file(model,
            joinpath(MODELS_DIR, "example1.cbf"))
    end

    @testset "Incompatible version" begin
        model = CBF.Model()
        @test_throws Exception MOI.read_from_file(model,
            joinpath(MODELS_DIR, "incompatible_version.cbf"))
    end

    @testset "$filename" for filename in ["bad_cone_string_A.cbf",
        "bad_cone_string_B.cbf", "bad_cone_string_C.cbf",
        "bad_cone_string_D.cbf"]
        model = CBF.Model()
        @test_throws Exception MOI.read_from_file(model,
            joinpath(MODELS_DIR, filename))
    end

    @testset "$filename" for filename in ["bad_power_dim_A.cbf",
        "bad_power_dim_B.cbf"]
        model = CBF.Model()
        @test_throws Exception MOI.read_from_file(model,
            joinpath(MODELS_DIR, filename))
    end

    @testset "$filename" for filename in ["corrupt_line_A.cbf",
        "corrupt_line_B.cbf", "corrupt_line_C.cbf", "corrupt_line_D.cbf",
        "corrupt_line_E.cbf"]
        model = CBF.Model()
        @test_throws Exception MOI.read_from_file(model,
            joinpath(MODELS_DIR, filename))
    end
end

@testset "Write errors" begin
    @testset "Quadratic objective" begin
        model = CBF.Model()
        MOIU.loadfromstring!(model, """
            variables: x
            minobjective: 1 * x * x
        """)
        @test_throws Exception MOI.write_to_file(model, CBF_TEST_FILE)
    end

    # TODO NLP not supported test.
end

write_read_models = [
    ("min SingleVariable", """
        variables: x
        minobjective: x
    """),
    ("min ScalarAffine", """
        variables: x, y
        minobjective: 1.2x + -1y + 1
    """),
    ("max SingleVariable", """
        variables: x
        maxobjective: x
    """),
    ("max ScalarAffine", """
        variables: x, y
        maxobjective: 1.2x + -1y + 1
    """),
    ("SingleVariable in Integer", """
        variables: x, y
        minobjective: 1.2x
        c1: y in Integer()
    """),
    ("VectorOfVariables in Zeros", """
        variables: x, y
        minobjective: x
        c1: [x, y] in Zeros(2)
    """),
    ("VectorOfVariables in Reals", """
        variables: x, y
        minobjective: x
        c1: [x, y] in Reals(2)
    """),
    ("VectorOfVariables in Nonnegatives", """
        variables: x, y
        minobjective: x
        c1: [x, y] in Nonnegatives(2)
    """),
    ("VectorOfVariables in Nonpositives", """
        variables: x, y
        minobjective: x
        c1: [y, x] in Nonpositives(2)
    """),
    ("VectorAffineFunction in Zeros", """
        variables: x, y
        minobjective: 1.2x
        c1: [x + 2y + -1.1, 0] in Zeros(2)
    """),
    ("VectorAffineFunction in Reals", """
        variables: x, y
        minobjective: 1.2x
        c1: [1x, 2y] in Reals(2)
    """),
    ("VectorAffineFunction in Nonnegatives", """
        variables: x, y
        minobjective: 1.2x
        c1: [1.1 * x, y + 1] in Nonnegatives(2)
    """),
    ("VectorAffineFunction in Nonpositives", """
        variables: x
        minobjective: 1.2x
        c1: [-1.1 * x + 1] in Nonpositives(1)
    """),
    ("VectorOfVariables in SecondOrderCone", """
        variables: x, y, z
        minobjective: x
        c1: [x, y, z] in SecondOrderCone(3)
    """),
    ("VectorOfVariables in RotatedSecondOrderCone", """
        variables: x, y, z
        minobjective: x
        c1: [x, y, z] in RotatedSecondOrderCone(3)
    """),
    ("VectorOfVariables in ExponentialCone", """
        variables: x, y, z
        minobjective: x
        c1: [x, y, z] in ExponentialCone()
    """),
    ("VectorOfVariables in DualExponentialCone", """
        variables: x, y, z
        minobjective: x
        c1: [x, y, z] in DualExponentialCone()
    """),
    ("VectorOfVariables in PowerCone", """
        variables: x, y, z
        minobjective: x
        c1: [x, y, z] in PowerCone(2.0)
    """),
    ("VectorOfVariables in DualPowerCone", """
        variables: x, y, z
        minobjective: x
        c1: [x, y, z] in DualPowerCone(2.0)
    """),
    ("VectorOfVariables in PositiveSemidefiniteConeTriangle", """
        variables: x, y, z
        minobjective: x
        c1: [x, y, z] in PositiveSemidefiniteConeTriangle(2)
    """),
    ("VectorAffineFunction in SecondOrderCone", """
        variables: x, y, z
        minobjective: 1.2x
        c1: [1.1x, y + 1, 2x + z] in SecondOrderCone(3)
    """),
    ("VectorAffineFunction in RotatedSecondOrderCone", """
        variables: x, y, z
        minobjective: 1.2x
        c1: [1.1x, y + 1, 2x + z] in RotatedSecondOrderCone(3)
    """),
    ("VectorAffineFunction in ExponentialCone", """
        variables: x, y, z
        minobjective: 1.2x
        c1: [1.1x, y + 1, 2x + z] in ExponentialCone()
    """),
    ("VectorAffineFunction in DualExponentialCone", """
        variables: x, y, z
        minobjective: 1.2x
        c1: [1.1x, y + 1, 2x + z] in DualExponentialCone()
    """),
    ("VectorAffineFunction in PowerCone", """
        variables: x, y, z
        minobjective: 1.2x
        c1: [1.1x, y + 1, 2x + z] in PowerCone(2.0)
    """),
    ("VectorAffineFunction in DualPowerCone", """
        variables: x, y, z
        minobjective: 1.2x
        c1: [1.1x, y + 1, 2x + z] in DualPowerCone(2.0)
    """),
    ("VectorAffineFunction in PositiveSemidefiniteConeTriangle", """
        variables: x, y, z
        minobjective: 1.2x
        c1: [1.1x, y + 1, 2x + z] in PositiveSemidefiniteConeTriangle(2)
    """),
]
@testset "Write/read $model_name" for (model_name, model_string) in
    write_read_models
    test_write_then_read(model_string)
end

example_models = [
    ("example_A.cbf", """
        variables: U, V, W, X, Y, Z, x, y, z
        minobjective: y + 2U + 2V + 2W + 2Y + 2Z
        c1: [U, V, W, X, Y, Z] in PositiveSemidefiniteConeTriangle(3)
        c2: [y + U + W + Z + -1, x + z + U + 2V + W + 2X + 2Y + Z + -0.5] in Zeros(2)
        c3: [y, x, z] in SecondOrderCone(3)
    """),
    ("example_B.cbf", """
        variables: X, Y, Z, x, y
        minobjective: 1 + x + y + X + Z
        c1: [X, Y, Z] in PositiveSemidefiniteConeTriangle(2)
        c2: [2Y + -1x + -1y] in Nonnegatives(1)
        c3: [3y + -1, x + y, 3x + -1] in PositiveSemidefiniteConeTriangle(2)
    """),
    ("example_C.cbf", """
        variables: a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p
        maxobjective: a + b + c + d + e + f + g + h + i + j + k + l + m + n + o + p + -1
        c1: [b] in Zeros(1)
        c2: [c] in Nonnegatives(1)
        c3: [d] in Nonpositives(1)
        c4: [e, f, g] in SecondOrderCone(3)
        c5: [h, i, j] in RotatedSecondOrderCone(3)
        c6: [m, l, k] in ExponentialCone()
        c7: [p, o, n] in DualExponentialCone()
    """),
    ("example_D.cbf", """
        variables: u, v, w, x, y, z
        maxobjective: w + z
        c1: [u, v, w] in PowerCone(0.5)
        c2: [x, y, z] in DualPowerCone(0.1)
        c3: [1, u, u + v] in PowerCone(0.8)
        c4: [1, y, x + y] in DualPowerCone(0.75)
        c5: u in Integer()
        c6: w in Integer()
        c7: y in Integer()
    """),
]
@testset "Read and write/read $model_name" for (model_name, model_string) in example_models
    test_read(joinpath(MODELS_DIR, model_name), model_string)
    test_write_then_read(model_string)
end

# Clean up.
sleep(1.0)  # Allow time for unlink to happen.
rm(CBF_TEST_FILE, force = true)


text:
# # LowLevelParticleFilters
# [![Build Status](https://travis-ci.org/baggepinnen/LowLevelParticleFilters.jl.svg?branch=master)](https://travis-ci.org/baggepinnen/LowLevelParticleFilters.jl)
# [![codecov](https://codecov.io/gh/baggepinnen/LowLevelParticleFilters.jl/branch/master/graph/badge.svg)](https://codecov.io/gh/baggepinnen/LowLevelParticleFilters.jl)

# This readme is auto generated from the file [src/example_lineargaussian.jl](https://github.com/baggepinnen/LowLevelParticleFilters.jl/blob/master/src/example_lineargaussian.jl) using [Literate.jl](https://github.com/fredrikekre/Literate.jl)

# # Types
# We provide a number of filter types
# - `ParticleFilter`: This filter is simple to use and assumes that both dynamics noise and measurement noise are additive.
# - `AuxiliaryParticleFilter`: This filter is identical to `ParticleFilter`, but uses a slightly different proposal mechanism for new particles.
# - `AdvancedParticleFilter`: This filter gives you more flexibility, at the expense of having to define a few more functions. More instructions on this type below.
# - `KalmanFilter`. Is what you would expect. Has the same features as the particle filters, but is restricted to linear dynamics and gaussian noise.
# - `UnscentedKalmanFilter`. Is also what you would expect. Has almost the same features as the kalman filters, but handle nonlinear dynamics and measurement model, still requires an additive Gaussian noise model.

# # Functionality
# - Filtering
# - Smoothing
# - Parameter estimation using ML or PMMH (Particle Marginal Metropolis Hastings)

# # Usage example
# This example demostrates how we set up the filters, both PF and KF, for a simple linear system.

# ## Particle filter
# Defining a particle filter is straightforward, one must define the distribution of the noise `df` in the dynamics function, `dynamics(x,u)` and the noise distribution `dg` in the measurement function `measurement(x)`. The distribution of the initial state `d0` must also be provided. An example for a linear Gaussian system is given below.
using LowLevelParticleFilters, LinearAlgebra, StaticArrays, Distributions, Plots

# Define problem

n = 2   # Dinemsion of state
m = 2   # Dinemsion of input
p = 2   # Dinemsion of measurements
N = 500 # Number of particles

const dg = MvNormal(p,1.0)          # Measurement noise Distribution
const df = MvNormal(n,1.0)          # Dynamics noise Distribution
const d0 = MvNormal(randn(n),2.0)   # Initial state Distribution

# Define random linear state-space system
Tr = randn(n,n)
const A = SMatrix{n,n}(Tr*diagm(0=>LinRange(0.5,0.95,n))/Tr)
const B = @SMatrix randn(n,m)
const C = @SMatrix randn(p,n)

# The following two functions are required by the filter
dynamics(x,u) = A*x .+ B*u
measurement(x) = C*x
vecvec_to_mat(x) = copy(reduce(hcat, x)') # Helper function

# We are now ready to define and use a filter
pf = ParticleFilter(N, dynamics, measurement, df, dg, d0)
xs,u,y = simulate(pf,100,df) # We can simulate the model that the pf represents
pf(u[1], y[1]) # Perform one filtering step using input u and measurement y
particles(pf) # Query the filter for particles, try weights(pf) or expweights(pf) as well
x̂ = weigthed_mean(pf) # using the current state
# If you want to perform filtering using vectors of inputs and measurements, try any of the functions
x,w,we,ll = forward_trajectory(pf, u, y) # Filter whole vectors of signals
x̂,ll = mean_trajectory(pf, u, y)
trajectorydensity(pf,x,w,y,xreal=xs)
# ![window](figs/trajdens.png)

# To see how the performance varies with the number of particles, we simulate several times. The following code simulates the system and performs filtering using the simulated measuerments. We do this for varying number of time steps and varying number of particles.

function run_test()
    particle_count = [10, 20, 50, 100, 200, 500, 1000]
    time_steps = [20, 100, 200]
    RMSE = zeros(length(particle_count),length(time_steps)) # Store the RMS errors
    propagated_particles = 0
    t = @elapsed for (Ti,T) = enumerate(time_steps)
        for (Ni,N) = enumerate(particle_count)
            montecarlo_runs = 2*maximum(particle_count)*maximum(time_steps) ÷ T ÷ N
            E = sum(1:montecarlo_runs) do mc_run
                pf = ParticleFilter(N, dynamics, measurement, df, dg, d0) # Create filter
                u = SVector{2,Float64}(randn(2))
                x = SVector{2,Float64}(rand(d0))
                y = SVector{2,Float64}(sample_measurement(pf,x,1))
                error = 0.0
                @inbounds for t = 1:T-1
                    pf(u, y) # Update the particle filter
                    x = dynamics(x,u) + SVector{2,Float64}(rand(df)) # Simulate the true dynamics and add some noise
                    y = SVector{2,Float64}(sample_measurement(pf,x,t)) # Simulate a measuerment
                    u = @SVector randn(2) # draw a random control input
                    error += sum(abs2,x-weigthed_mean(pf))
                end # t
                √(error/T)
            end # MC
            RMSE[Ni,Ti] = E/montecarlo_runs
            propagated_particles += montecarlo_runs*N*T
            @show N
        end # N
        @show T
    end # T
    println("Propagated $propagated_particles particles in $t seconds for an average of $(propagated_particles/t/1000) particles per millisecond")
    return RMSE
end

@time RMSE = run_test()
# Propagated 8400000 particles in 3.568745383 seconds for an average of 2353.7683691344473 particles per millisecond

# We then plot the results
time_steps     = [20, 100, 200]
particle_count = [10, 20, 50, 100, 200, 500, 1000]
nT             = length(time_steps)
leg            = reshape(["$(time_steps[i]) time steps" for i = 1:nT], 1,:)
plot(particle_count,RMSE,xscale=:log10, ylabel="RMS errors", xlabel=" Number of particles", lab=leg)
gui()
# ![window](figs/rmse.png)

# # Smoothing
# We also provide a particle smoother, based on forward filtering, backward simulation (FFBS)

N     = 2000 # Number of particles
T     = 200 # Number of time steps
M     = 100 # Number of smoothed backwards trajectories
pf    = ParticleFilter(N, dynamics, measurement, df, dg, d0)
du    = MvNormal(2,1) # Control input distribution
x,u,y = simulate(pf,T,du) # Simulate trajectory using the model in the filter
tosvec(y) = reinterpret(SVector{length(y[1]),Float64}, reduce(hcat,y))[:] |> copy
x,u,y = tosvec.((x,u,y))

xb,ll = smooth(pf, M, u, y) # Sample smooting particles
xbm   = smoothed_mean(xb)   # Calculate the mean of smoothing trajectories
xbc   = smoothed_cov(xb)    # And covariance
xbt   = smoothed_trajs(xb)  # Get smoothing trajectories
xbs   = [diag(xbc) for xbc in xbc] |> vecvec_to_mat .|> sqrt
plot(xbm', ribbon=2xbs, lab="PF smooth")
plot!(vecvec_to_mat(x), l=:dash, lab="True")
# ![window](figs/smooth.svg)

# We can plot the particles themselves as well
plot(vecvec_to_mat(x), l=(4,), layout=(2,1), show=false)
scatter!(xbt[1,:,:]', subplot=1, show=false, m=(1,:black, 0.5), lab="")
scatter!(xbt[2,:,:]', subplot=2, m=(1,:black, 0.5), lab="")

# ![window](figs/smooth.png)

# # Kalman filter
# A Kalman filter is easily created using the constructor. Many of the functions defined for particle filters, are defined also for Kalman filters, e.g.:

eye(n) = Matrix{Float64}(I,n,n)
kf     = KalmanFilter(A, B, C, 0, eye(n), eye(p), MvNormal([1.,1.]))
ukf    = UnscentedKalmanFilter(dynamics, measurement, eye(n), eye(p), MvNormal([1.,1.]))
xf,xt,R,Rt,ll = forward_trajectory(kf, u, y) # filtered, prediction, pred cov, filter cov, loglik
xT,R,lls = smooth(kf, u, y) # Smoothed state, smoothed cov, loglik
# It can also be called in a loop like the `pf` above
#md for t = 1:T
#md     kf(u,y) # Performs both correct and predict!!
#md     ## alternatively
#md     ll += correct!(kf, y, t) # Returns loglik
#md     x   = state(kf)
#md     R   = covariance(kf)
#md     predict!(kf, u, t)
#md end

# # Troubleshooting
# Tuning a particle filter can be quite the challenge. To assist with this, we provide som visualization tools
debugplot(pf,u,y, runall=true, xreal=x) # does not work well with gr() as backend, try pyplot()
# ![window](figs/debugplot.png)
#md Time     Surviving    Effective nbr of particles
#md --------------------------------------------------------------
#md t:     1   1.000    1000.0
#md t:     2   1.000     551.0
#md t:     3   1.000     453.0
#md t:     4   1.000     384.3
#md t:     5   1.000     340.9
#md t:     6   1.000     310.5
#md t:     7   1.000     280.0
#md t:     8   1.000     265.9
# The plot displays all states and all measurements. The heatmap in the background represents the weighted particle distributions per time step. For the measurement sequences, the heatmap represent the distibutions of predicted measurements. The blue dots corresponds to measured values. In this case, we simulated the data and we had access to states as well, if we do not have that, just omit `xreal`.
# You can also manually step through the time-series using
# - `commandplot(pf,u,y; kwargs...)`
# For options to the debug plots, see `?pplot`.

# # Parameter estimation
# We provide som basic functionality for maximum likelihood estimation and MAP estimation
# ## ML estimation
# Plot likelihood as function of the variance of the dynamics noise
svec = exp10.(LinRange(-1.5,1.5,60))
llspf = map(svec) do s
    df = MvNormal(n,s)
    pfs = ParticleFilter(2000, dynamics, measurement, df, dg, d0)
    loglik(pfs,u,y)
end
plot(svec, llspf, xscale=:log10, title="Log-likelihood", xlabel="Dynamics noise standard deviation", lab="PF")
vline!([svec[findmax(llspf)[2]]], l=(:dash,:blue), primary=false)
# We can do the same with a Kalman filter
eye(n) = Matrix{Float64}(I,n,n)
llskf = map(svec) do s
    kfs = KalmanFilter(A, B, C, 0, s^2*eye(n), eye(p), d0)
    loglik(kfs,u,y)
end
plot!(svec, llskf, yscale=:identity, xscale=:log10, lab="Kalman", c=:red)
vline!([svec[findmax(llskf)[2]]], l=(:dash,:red), primary=false)
# ![window](figs/svec.png)
# as we can see, the result is quite noisy due to the stochastic nature of particle filtering.

# ### Smoothing using KF
kf = KalmanFilter(A, B, C, 0, eye(n), eye(p), MvNormal(2,1))
xf,xh,R,Rt,ll = forward_trajectory(kf, u, y) # filtered, prediction, pred cov, filter cov, loglik
xT,R,lls = smooth(kf, u, y) # Smoothed state, smoothed cov, loglik

# Plot and compare PF and KF
plot(vecvec_to_mat(xT), lab="Kalman smooth", layout=2)
plot!(xbm', lab="pf smooth")
plot!(vecvec_to_mat(x), lab="true")
# ![window](figs/smoothtrajs.svg)

# ## MAP estiamtion
# To solve a MAP estimation problem, we need to define a function that takes a parameter vector and returns a particle filter

filter_from_parameters(θ,pf=nothing) = ParticleFilter(N, dynamics, measurement, MvNormal(n,exp(θ[1])), MvNormal(p,exp(θ[2])), d0)
# The call to `exp` on the parameters is so that we can define log-normal priors
priors = [Normal(0,2),Normal(0,2)]
# Now we call the function `log_likelihood_fun` that returns a function to be minimized
ll = log_likelihood_fun(filter_from_parameters,priors,u,y)
# Since this is a low-dimensional problem, we can plot the LL on a 2d-grid
function meshgrid(a,b)
    grid_a = [i for i in a, j in b]
    grid_b = [j for i in a, j in b]
    grid_a, grid_b
end
Nv       = 20
v        = LinRange(-0.7,1,Nv)
llxy     = (x,y) -> ll([x;y])
VGx, VGy = meshgrid(v,v)
VGz      = llxy.(VGx, VGy)
heatmap(VGz, xticks=(1:Nv,round.(v,digits=2)),yticks=(1:Nv,round.(v,digits=2)), xlabel="sigma v", ylabel="sigma w") # Yes, labels are reversed

# ![window](figs/heatmap.svg)

# Something seems to be off with this figure as the hottest spot is not really where we would expect it

# Optimization of the log likelihood can be done by, e.g., global/black box methods, see [BlackBoxOptim.jl](https://github.com/robertfeldt/BlackBoxOptim.jl). Standard tricks apply, such as performing the parameter search in log-space etc.

# ## Bayesian inference using PMMH
# This is pretty cool. We procede like we did for MAP above, but when calling the function `metropolis`, we will get the entire posterior distribution of the parameter vector, for the small cost of a massive increase in computational cost.
N = 1000
filter_from_parameters(θ,pf=nothing) = AuxiliaryParticleFilter(N, dynamics, measurement, MvNormal(n,exp(θ[1])), MvNormal(p,exp(θ[2])), d0)
# The call to `exp` on the parameters is so that we can define log-normal priors
priors = [Normal(0,2),Normal(0,2)]
ll     = log_likelihood_fun(filter_from_parameters,priors,u,y)
θ₀     = log.([1.,1.]) # Starting point
# We also need to define a function that suggests a new point from the "proposal distribution". This can be pretty much anything, but it has to be symmetric since I was lazy and simplified an equation.
draw   = θ -> θ .+ rand(MvNormal(0.05ones(2)))
burnin = 200
@info "Starting Metropolis algorithm"
@time theta, lls = metropolis(ll, 2000, θ₀, draw) # Run PMMH for 2000  iterations, takes about half a minute on my laptop
thetam = reduce(hcat, theta)'[burnin+1:end,:] # Build a matrix of the output (was vecofvec)
histogram(exp.(thetam), layout=(3,1)); plot!(lls[burnin+1:end], subplot=3) # Visualize
# ![window](figs/histogram.svg)

# If you are lucky, you can run the above threaded as well. I tried my best to make particle fitlers thread safe with their own rngs etc., but your milage may vary.
#md @time thetalls = LowLevelParticleFilters.metropolis_threaded(burnin, ll, 500, θ₀, draw)
#md histogram(exp.(thetalls[:,1:2]), layout=3)
#md plot!(thetalls[:,3], subplot=3)

# # AdvancedParticleFilter
# The `AdvancedParticleFilter` type requires you to implement the same functions as the regular `ParticleFilter`, but in this case you also need to handle sampling from the noise distributions yourself.
# The function `dynamics` must have a method signature like below. It must provide one method that accepts state vector, control vector, time and `noise::Bool` that indicates whether or not to add noise to the state. If noise should be added, this should be done inside `dynamics` An example is given below
function dynamics(x,u,t,noise=false) # It's important that this defaults to false
    x = A*x .+ B*u # A simple dynamics model
    if noise
        x += rand(df)
    end
    x
end
# The `measurement_likelihood` function must have a method accepting state, measurement and time, and returning the log-likelihood of the measurement given the state, a simple example below:
function measurement_likelihood(x,y,t)
    logpdf(dg, C*x-y) # A simple linear measurement model with normal additive noise
end
# This gives you very high flexibility. The noise model in either function can, for instance, be a function of the state, something that is not possible for the simple `ParticleFilter`
# To be able to simulate the `AdvancedParticleFilter` like we did with the simple filter above, the `measurement` method with the signature `measurement(x,t,noise=false)` must be available and return a sample measurement given state (and possibly time). For our example measurement model above, this would look like this
measurement(x,t,noise=false) = C*x + noise*rand(dg)
# We now create the `AdvancedParticleFilter` and use it in the same way as the other filters:
apf = AdvancedParticleFilter(N, dynamics, measurement, measurement_likelihood, df, d0)
x,w,we,ll = forward_trajectory(apf, u, y)
trajectorydensity(apf, x, we, y, xreal=xs)

# We can even use this type as an AuxiliaryParticleFilter
apfa = AuxiliaryParticleFilter(apf)
x,w,we,ll = forward_trajectory(apfa, u, y)
trajectorydensity(apfa, x, we, y, xreal=xs)
dimensiondensity(apfa, x, we, y, 1, xreal=xs) # Same as above, but only plots a single dimension

# # High performance Distributions
# When `using LowLevelParticleFilters`, a number of methods related to distributions are defined for static arrays, making `logpdf` etc. faster. We also provide a new kind of distribution: `TupleProduct <: MultivariateDistribution` that behaves similarly to the `Product` distribution. The `TupleProduct` however stores the individual distributions in a tuple, has compile-time known length and supports `Mixed <: ValueSupport`, meaning that it can be a product of both `Continuous` and `Discrete` dimensions, somthing not supported by the standard `Product`. Example
#md using BenchmarkTools, LowLevelParticleFilters, Distributions
#md dt = TupleProduct((Normal(0,2), Normal(0,2), Binomial())) # Mixed value support
#md # A small benchmark
#md ```julia
#md sv = @SVector randn(2)
#md d = Product([Normal(0,2), Normal(0,2)])
#md dt = TupleProduct((Normal(0,2), Normal(0,2)))
#md dm = MvNormal(2, 2)
#md @btime logpdf($d,$(Vector(sv))) # 32.449 ns (1 allocation: 32 bytes)
#md @btime logpdf($dt,$(Vector(sv))) # 21.141 ns (0 allocations: 0 bytes)
#md @btime logpdf($dm,$(Vector(sv))) # 48.745 ns (1 allocation: 96 bytes)
#md
#md @btime logpdf($d,$sv) # 22.651 ns (0 allocations: 0 bytes)
#md @btime logpdf($dt,$sv) # 0.021 ns (0 allocations: 0 bytes)
#md @btime logpdf($dm,$sv) # 0.021 ns (0 allocations: 0 bytes)
#md # Without loading `LowLevelParticleFilters`, the timing for the native distributions are the following
#md @btime logpdf($d,$sv) # 32.621 ns (1 allocation: 32 bytes)
#md @btime logpdf($dm,$sv) # 46.415 ns (1 allocation: 96 bytes)
#md ```

#jl # Compile using Literate.markdown("example_lineargaussian.jl", "..", name="README", documenter=false)


text:
using ACTRSimulators, Test, ACTRModels, DataFrames, Random
import ACTRSimulators: start!, press_key!
Random.seed!(8985)
include("task.jl")

function can_attend(actr)
    c1(actr) = !isempty(actr.visual_location.buffer)
    c2(actr) = !actr.visual.state.busy
    return c1, c2
end

function can_encode(actr)
    c1(actr) = !isempty(actr.visual.buffer)
    c2(actr) = !actr.imaginal.state.busy
    return c1, c2
end

function can_respond1(actr)
    c1(actr) = !isempty(actr.imaginal.buffer)
    c2(actr) = !actr.motor.state.busy
    return c1, c2
end

function can_respond2(actr)
    c1(actr) = !isempty(actr.imaginal.buffer)
    c2(actr) = !actr.motor.state.busy
    return c1, c2
end

function attend_action(actr, task)
    buffer = actr.visual_location.buffer
    chunk = deepcopy(buffer[1])
    clear_buffer!(actr.visual_location)
    attending!(actr, chunk)
    return nothing
end

function encode_action(actr, task)
    buffer = actr.visual.buffer
    chunk = deepcopy(buffer[1])
    clear_buffer!(actr.visual)
    encoding!(actr, chunk)
    return nothing
end

function motor_action1(actr, task)
    buffer = actr.imaginal.buffer
    chunk = deepcopy(buffer[1])
    clear_buffer!(actr.imaginal)
    key = "j"
    responding!(actr, task, key)
    return nothing
end

function motor_action2(actr, task)
    buffer = actr.imaginal.buffer
    chunk = deepcopy(buffer[1])
    clear_buffer!(actr.imaginal)
    key = "f"
    responding!(actr, task, key)
    return nothing
end

parms = (
    utility_noise = true,
)

scheduler = ACTRScheduler(;model_trace=false)
task = SimpleTask(;scheduler)
procedural = Procedural()
T = vo_to_chunk() |> typeof
visual_location = VisualLocation(buffer=T[])
visual = Visual(buffer=T[])
visicon = VisualObject[]
motor = Motor()
memory = [Chunk(;animal=:dog), Chunk(;animal=:cat)]
declarative = Declarative(;memory)
actr = ACTR(;
    scheduler,
    procedural,
    visual_location,
    visual,
    motor,
    declarative,
    visicon,
    parms...
)

rule1 = Rule(;conditions=can_attend, action=attend_action, actr, task, name="Attend")
push!(procedural.rules, rule1)

rule2 = Rule(;conditions=can_encode, action=encode_action, actr, task, name="Encode")
push!(procedural.rules, rule2)

rule3 = Rule(;conditions=can_respond1, action=motor_action1, actr, task, name="Respond j")
push!(procedural.rules, rule3)

rule4 = Rule(;conditions=can_respond2, action=motor_action2, actr, task, name="Respond f")
push!(procedural.rules, rule4)

map(_->run!(actr, task), 1:1000)

p_j = mean(task.data.key .== "j")

@test p_j ≈ .5 atol = .01


text:
"""
    precessing_nutating_example([Ωₒ], [Ωₚ], [α], [α̇], [ν], [R₀])

Return an exact quaternionic function of time representing nutating
precessional orbital motion, as described in Sec. 6.2 of [this
paper](https://arxiv.org/abs/1604.08139).  This example is useful because it
provides physically realistic but *very* complicated motion, while still being
simple to code up and differentiate analytically.

The output represents the rotation of the line joining the orbiting bodies and
their angular velocity.  The following are the input parameters, along with
their default values and physical interpretations:

  * `Ωₒ=2π/1_000`: The orbital frequency
  * `Ωₚ=2π/10_000`: The precessional frequency
  * `α=π/8`: The opening angle of the precession cone at ``t=0``
  * `α̇=2α/100_000`: The rate of opening of the precession cone
  * `ν=π/80`: The angle of nutation
  * `R₀=exp(-3α*imx/10)`: Overal rotation of the system

The default values are chosen to be typical of a (potentially) real precessing
binary black hole system shortly before merger.

The returned objects are three functions of time: `R`, `ω⃗`, and `Ṙ`, which
return the orientation as a `Rotor`, followed by the angular velocity as a
`QuatVec`, and the time-derivative of `R` as a `Quaternion`.

# Example
```jldoctest
julia> R, ω⃗, Ṙ = precessing_nutating_example();

julia> R(12.34)
0.9944579779058746 + 0.09804177421238346𝐢 - 0.00084850453525312𝐣 + 0.03795287510453948𝐤
julia> ω⃗(345.67)
0.0 + 0.00046343007342867023𝐢 - 0.0007032818419003181𝐣 + 0.006214814810035088𝐤
julia> ϵ = 1e-6; (R(ϵ) - R(-ϵ)) / 2ϵ  # Approximate derivative at t=0
-3.8491432263754177e-7 + (3.9080960689830135e-6)𝐢 - (6.861695854245619e-5)𝐣 + 0.003076329202503836𝐤
julia> Ṙ(0)
-3.8491240998153413e-7 + (3.9080812828476746e-6)𝐢 - (6.861695854245637e-5)𝐣 + 0.0030763292025038354𝐤
```

"""
function precessing_nutating_example(Ωₒ=2π/1_000, Ωₚ=2π/10_000, α=π/8, α̇=2α/100_000, ν=π/80, R₀=exp(-3α*imx/10))
    R₁(t) = exp(Ωₒ*t*imz/2)
    R₂(t) = exp((α + α̇*t)*imx/2)
    R₃(t) = exp(Ωₚ*t*imz/2)
    R₄(t) = exp(ν*imx/2)
    R(t) = R₀ * R₁(t) * R₄(t) * inv(R₁(t)) * R₃(t) * R₂(t) * inv(R₃(t)) * R₁(t)
    Ṙ₀ = 0imz
    Ṙ₁(t) = R₁(t) * (Ωₒ*imz/2)
    Ṙ₂(t) = R₂(t) * (α̇*imx/2)
    Ṙ₃(t) = R₃(t) * (Ωₚ*imz/2)
    Ṙ₄(t) = 0imz
    Ṙ(t) = (
        Ṙ₀ * R₁(t) * R₄(t) * inv(R₁(t)) * R₃(t) * R₂(t) * inv(R₃(t)) * R₁(t) +
        R₀ * Ṙ₁(t) * R₄(t) * inv(R₁(t)) * R₃(t) * R₂(t) * inv(R₃(t)) * R₁(t) +
        R₀ * R₁(t) * Ṙ₄(t) * inv(R₁(t)) * R₃(t) * R₂(t) * inv(R₃(t)) * R₁(t) +
        -R₀ * R₁(t) * R₄(t) * inv(R₁(t)) * Ṙ₁(t) * inv(R₁(t)) * R₃(t) * R₂(t) * inv(R₃(t)) * R₁(t) +
        R₀ * R₁(t) * R₄(t) * inv(R₁(t)) * Ṙ₃(t) * R₂(t) * inv(R₃(t)) * R₁(t) +
        R₀ * R₁(t) * R₄(t) * inv(R₁(t)) * R₃(t) * Ṙ₂(t) * inv(R₃(t)) * R₁(t) +
        -R₀ * R₁(t) * R₄(t) * inv(R₁(t)) * R₃(t) * R₂(t) * inv(R₃(t)) * Ṙ₃(t) * inv(R₃(t)) * R₁(t) +
        R₀ * R₁(t) * R₄(t) * inv(R₁(t)) * R₃(t) * R₂(t) * inv(R₃(t)) * Ṙ₁(t)
    )
    ω⃗(t) = 2 * QuatVec(Ṙ(t) * inv(R(t)))
    (R, ω⃗, Ṙ)
end


text:
struct DirectPolicyOptimizationProblem
	# problem(s)
	prob

	# objective
	obj

	# constraints
	con

	# number of samples
	Nn # 2n
	Nd # 2d

	# policy
	policy

	# disturbance trajectory
	dist

	# resampling
	sample

	# number of variables
	num_var
	num_con

	# indices
	idx
end

function dpo(
		prob_nom, prob_mean, prob_sample,
		obj,
		policy,
		dist,
		sample;
		Nn = 2 * prob_mean.model.n, Nd = 2 * prob_mean.model.d)

	N = Nn + Nd

	# problem(s)
	prob = DPOProblems(prob_nom, prob_mean, prob_sample)

	# indices
	idx, num_var = dpo_indices(prob, policy.p, Nn, Nd)

	# constraints
	con_dynamics = sample_dynamics_constraints(prob, Nn, Nd)
	con_policy = policy_constraints(prob, Nn, control_dim = policy.output)
	con = [con_dynamics, con_policy]

	num_con = prob.nom.num_con + sum([prob.sample[i].num_con for i = 1:Nn])
	num_con += con_dynamics.n + con_policy.n

	return DirectPolicyOptimizationProblem(
		prob, obj, con,
		Nn, Nd,
		policy,
		dist,
		sample,
		num_var, num_con,
		idx)
end

function dpo_problem(
		prob_nom, prob_mean, prob_sample,
		obj,
		policy,
		dist,
		sample;
		Nn = 2 * prob_mean.model.n, Nd = 2 * prob_mean.model.d)

	prob = dpo(
			prob_nom, prob_mean, prob_sample,
			obj,
			policy,
			dist,
			sample;
			Nn = Nn, Nd = Nd)

	prob_moi = moi_problem(prob)

	return prob_moi
end

struct DPOProblems
	nom
	mean
	sample
end

struct DPOIndices
	nom
	mean
	sample

	xt
	ut

	θ
	policy
end

function dpo_indices(prob::DPOProblems, p, Nn, Nd)
    T = prob.nom.T

    # nominal
    idx_nom = 1:prob.nom.num_var

    # mean
    idx_mean = prob.nom.num_var .+ (1:prob.mean.num_var)

    # sample
    shift = prob.nom.num_var + prob.mean.num_var
    idx_sample = []
    for i = 1:Nn
        push!(idx_sample, shift .+ (1:prob.sample[i].num_var))
        shift += prob.sample[i].num_var
    end

	xt = [vcat([idx_sample[i][prob.sample[i].idx.x[t]] for i = 1:Nn]...) for t = 1:T]
	ut = [vcat([idx_sample[i][prob.sample[i].idx.u[t]] for i = 1:Nn]...) for t = 1:T-1]

    # policy
    idx_θ = [(t - 1) * p .+ (1:p) for t = 1:T-1]
    n_policy = p * (T - 1)
    idx_policy = shift .+ (1:n_policy)
    shift += n_policy

	num_var = shift

    return DPOIndices(idx_nom, idx_mean, idx_sample,
		xt, ut,
		idx_θ, idx_policy), num_var
end

function primal_bounds(dpo::DirectPolicyOptimizationProblem)
    T = dpo.prob.nom.T

    idx_nom = dpo.idx.nom
	idx_mean = dpo.idx.mean
	idx_sample = dpo.idx.sample

    Zl = -Inf * ones(dpo.num_var)
    Zu = Inf * ones(dpo.num_var)

	# nominal
	Zl[idx_nom], Zu[idx_nom] = primal_bounds(dpo.prob.nom)

	# sample
	for i = 1:dpo.Nn
		Zl[idx_sample[i]], Zu[idx_sample[i]] = primal_bounds(dpo.prob.sample[i])
	end

    return Zl, Zu
end

function constraint_bounds(dpo::DirectPolicyOptimizationProblem)
	idx_nom = dpo.idx.nom
	idx_sample = dpo.idx.sample

    cl = zeros(dpo.num_con)
    cu = zeros(dpo.num_con)

	shift = 0

	# nominal
    cu[dpo.prob.nom.con.ineq] .= Inf
	shift += dpo.prob.nom.num_con

	# sample
	for i = 1:dpo.Nn
		cu[(shift .+ (1:dpo.prob.sample[i].num_con))[dpo.prob.sample[i].con.ineq]] .= Inf
		shift += dpo.prob.sample[i].num_con
	end

    return cl, cu
end

function eval_objective(dpo::DirectPolicyOptimizationProblem, Z)
	# nominal
	J = objective(view(Z, dpo.idx.nom), dpo.prob.nom)

	# samples
	for i = 1:dpo.Nn
		J += objective(view(Z, dpo.idx.sample[i]), dpo.prob.sample[i])
	end

	# sample objective
	J += objective(Z, dpo.obj, dpo.prob, dpo.idx, dpo.Nn, dpo.Nd)

	return J
end

function eval_objective_gradient!(∇J, Z, dpo::DirectPolicyOptimizationProblem)
    ∇J .= 0.0

	# nominal
    objective_gradient!(view(∇J, dpo.idx.nom),
		view(Z, dpo.idx.nom), dpo.prob.nom)

	# samples
	for i = 1:dpo.Nn
		objective_gradient!(view(∇J, dpo.idx.sample[i]),
			view(Z, dpo.idx.sample[i]), dpo.prob.sample[i])
	end

	# sample objective
	objective_gradient!(∇J, Z, dpo.obj, dpo.prob, dpo.idx, dpo.Nn, dpo.Nd)

    return nothing
end

function eval_constraint!(c, Z, dpo::DirectPolicyOptimizationProblem)
	shift = 0

	# nominal
    constraints!(view(c, 1:dpo.prob.nom.num_con),
		view(Z, dpo.idx.nom), dpo.prob.nom)
	shift += dpo.prob.nom.num_con

	# sample
	for i = 1:dpo.Nn
		constraints!(view(c, shift .+ (1:dpo.prob.sample[i].num_con)),
			view(Z, dpo.idx.sample[i]), dpo.prob.sample[i])
		shift += dpo.prob.sample[i].num_con
	end

	# sample dynamics
	constraints!(view(c, shift .+ (1:dpo.con[1].n)), Z, dpo.con[1],
		dpo.prob, dpo.idx, dpo.Nn, dpo.Nd, dpo.dist, dpo.sample)
	shift += dpo.con[1].n

	# policy
	constraints!(view(c, shift .+ (1:dpo.con[2].n)), Z, dpo.con[2],
		dpo.policy, dpo.prob, dpo.idx, dpo.Nn)
	shift += dpo.con[2].n

    return nothing
end

function eval_constraint_jacobian!(∇c, Z, dpo::DirectPolicyOptimizationProblem)
	shift = 0

	# nominal
	len_nom = length(sparsity_jacobian(dpo.prob.nom))
    constraints_jacobian!(view(∇c, 1:len_nom), view(Z, dpo.idx.nom), dpo.prob.nom)
	shift += len_nom

	# samples
	for i = 1:dpo.Nn
		len = length(sparsity_jacobian(dpo.prob.sample[i]))
		constraints_jacobian!(view(∇c, shift .+ (1:len)),
			view(Z, dpo.idx.sample[i]), dpo.prob.sample[i])
		shift += len
	end

	# sample dynamics
	len = length(constraints_sparsity(dpo.con[1], dpo.prob, dpo.idx,
		dpo.Nn, dpo.Nd))
	constraints_jacobian!(view(∇c, shift .+ (1:len)), Z, dpo.con[1],
		dpo.prob, dpo.idx,
		dpo.Nn, dpo.Nd,
		dpo.dist, dpo.sample)
	shift += len

	# policy
	len = length(constraints_sparsity(dpo.con[2], dpo.policy, dpo.prob, dpo.idx,
		dpo.Nn))
	constraints_jacobian!(view(∇c, shift .+ (1:len)), Z, dpo.con[2],
		dpo.policy,
		dpo.prob, dpo.idx, dpo.Nn)
	shift += len

    return nothing
end

function sparsity_jacobian(dpo::DirectPolicyOptimizationProblem)

	# nominal
    spar = constraints_sparsity(dpo.prob.nom)
	con_shift = dpo.prob.nom.num_con

	# samples
	shift_col = dpo.prob.nom.num_var + dpo.prob.mean.num_var
	for i = 1:dpo.Nn
		spar = collect([spar...,
			constraints_sparsity(dpo.prob.sample[i],
				shift_row = con_shift,
				shift_col = shift_col)...])
		con_shift += dpo.prob.sample[i].num_con
		shift_col += dpo.prob.sample[i].num_var
	end
	shift_col = 0

	# sample dynamics
	spar = collect([spar...,
					constraints_sparsity(dpo.con[1], dpo.prob, dpo.idx,
						dpo.Nn, dpo.Nd, shift_row = con_shift)...])
	con_shift += dpo.con[1].n

	# policy
	spar = collect([spar...,
					constraints_sparsity(dpo.con[2], dpo.policy,
					dpo.prob, dpo.idx,
					dpo.Nn, shift_row = con_shift)...])
	con_shift += dpo.con[2].n

	return spar
end

"""
	pack nominal trajectory and policy parameters into DPO vector
"""
function pack(z, K, prob::DirectPolicyOptimizationProblem)
	z0 = zeros(prob.num_var)
	z0[prob.idx.nom] = z
	z0[prob.idx.mean] = z

	for i = 1:prob.Nn
		z0[prob.idx.sample[i]] = z
	end

	for t = 1:prob.prob.nom.T-1
		z0[prob.idx.policy[prob.idx.θ[t]]] = vec(copy(K[t]))
	end

	return z0
end

pack(z, K, prob::MOIProblem) = pack(z, K, prob.prob)

"""
	get policy parameters from DPO vector
"""
function get_policy(z, prob::DirectPolicyOptimizationProblem)
	return [reshape(z[prob.idx.policy[prob.idx.θ[t]]],
			prob.policy.output, prob.policy.input) for t = 1:prob.prob.nom.T-1]
end

get_policy(z, prob::MOIProblem) = get_policy(z, prob.prob)


text:
struct Edge
    id::Int64

    from::BaseVertex{Edge}
    to::BaseVertex{Edge}

    cost::Int64
    demand::Int64
end

const Vertex = BaseVertex{Edge}

isRequired(edge::Edge) = edge.demand > 0


text:
using Distributed
# Get one process per vCPU (set less processes if RAM consumption is a concern)
num_processes = Sys.CPU_THREADS
if nprocs() < num_processes
    addprocs(num_processes - nprocs())
end

using ParallelMCWF, QuantumOptics, LinearAlgebra
# Operates the tensor product of all arguments
⨂(x) = length(x) > 1 ? reduce(⊗, x) : x[1];
# Operates the sum of all arguments
∑(x) = sum(x);
##

# SOME DEFINITIONS
# Simulation parameters
nb_trajs = 500;        # Number of trajectories
tspan = Array(range(0,stop=30,length=5));
n_site_max = 5;         # Local Hilbert space cut-off

# Model parameters
L = 3               # Length of the Bose-Hubbard chain
κ1 = 1.             # 1-body dissipation rate (To populate both parities)
κ2 = 1κ1            # 2-body dissipation rate
J  = 50κ1           # Tunneling rate
Δ  = -J             # Cavity detuning
U  = 100κ1          # Nonlinearity strength
F  = 0.             # 1-body driving strength
G  = 50κ1           # 2-photon driving strength

# Hamiltonian
lb = FockBasis(n_site_max); # Local Hilbert space basis
gb = lb⊗lb⊗lb;              # Global Hilbert space basis
Id = one(lb);               # Local Identity operator
la = destroy(lb);           # Local annihilation operator
lN = number(lb);            # Local number operator
# Local Hamiltonian (for a single site)
lH = -Δ*lN + U/2 *dagger(la)*lN*la + F*(la + dagger(la)) + G/2 *(la^2 + dagger(la)^2);
# Tight-binding interaction between the three sites
# This term is computed by summing all possible circular permutations of â†⊗â⊗𝟙̂⊗...⊗𝟙̂
# The Hermitian conjugate is latter added so as to have a Hermitian operator.
Hj = J/2 * ∑([⨂(circshift([dagger(la), la, fill(Id,L-2)...], i)) for i in 1:L]);
# Total Hamiltonian: sum of each site's local Hamiltonian plus the tunneling term and its Hermitian conjugate
# Local Hamiltonians are embedded into the global Hilbert space using the same technique
H = ∑([⨂(circshift([lH, fill(Id,L-1)...], i)) for i in 1:L]) + Hj + dagger(Hj);

# Some observables
lΠ = diagonaloperator(lb,Array(exp.(im*π .* diag(lN.data)))); # Local parity operator
# One could equally have used lΠ = exp(im*π*dense(lN)) but it has very poor performance
gN = ∑([⨂(circshift([lN, fill(Id,L-1)...], i)) for i in 1:L]); # Global number operator
gΠ = diagonaloperator(gb,Array(exp.(im*π .* diag(gN.data))));
##

# MONTE CARLO
# Liouvillians are represented as Arrays of jump operators
# 1-body dissipation jump operators
κ1𝒟 = sqrt(κ1) * [⨂(circshift([la, fill(Id,L-1)...], i)) for i in 1:L]
# 1-body dissipation jump operators
κ2𝒟 = sqrt(κ2) * [⨂(circshift([la^2, fill(Id,L-1)...], i)) for i in 1:L]

# Single ℤ₂-symmetric cavity
println("Single Z2-symmetric cavity:")
ψ0 = fockstate(lb,0) # Start from vacuum
G1  = 100κ1;
H1 = U/2 *dagger(la)*lN*la + F*(la + dagger(la)) + G1/2 *(la^2 + dagger(la)^2);
t, trajs = pmcwf(tspan, ψ0, H1, [sqrt(κ1)*la, sqrt(κ2)*la^2]; Ntrajectories=nb_trajs, parallel_type=:pmap)
# Steady-state kets
kets = [trajs[i][end] for i in 1:nb_trajs]
println("Parity = ",real.(kets_to_obs(lΠ, kets)))
ρ = kets_to_dm(kets)
println("Von Neumann entropy: S = ln($(exp(entropy_vn(ρ))))")

# L=3 global-Z2-symmetric chain
println("L=3 global-Z2-symmetric chain:")
ψ0 = fockstate(lb,0) ⊗ fockstate(lb,0) ⊗ fockstate(lb,0)
t, trajs = pmcwf(tspan, ψ0, H, [κ1𝒟; κ2𝒟]; Ntrajectories=nb_trajs, parallel_type=:pmap)
# Steady-state kets
kets = [trajs[i][end] for i in 1:nb_trajs]
println("Global parity = ",real.(kets_to_obs(gΠ, kets)))
ρ = kets_to_dm(kets)
println("Von Neumann entropy: S = ln($(exp(entropy_vn(ρ))))")

diag(ptrace(ρ, [2,3]).data)

using BenchmarkTools
@btime timeevolution.mcwf(tspan, ψ0, H, [κ1𝒟; κ2𝒟])