text:
module BallisticDeposition

export bdeposition!, bdfull!, bdiso!, bdisofull!, bdrelax!, bdheights

"""
Apply ballistic deposition to a 1D surface with periodic boundary conditions

# Arguments

- `surface::Vector{Integer}`: the 1D array representing the surface
- `time::Integer`: the number of time steps (or the number of deposited particles)

# Returns

the surface heights (Vector) after applying the ballistic deposition
"""
function bdeposition!(surface::Vector{<:Integer}, time::Integer)
    len = length(surface)
    starting_indecies = rand(1:len, time)

    for i in starting_indecies
        # finding the left and right indecies with periodic boundary conditions
        left = i - 1 + ((len + 1 - i) ÷ len) * len
        right = i + 1 - (i ÷ len) * len

        # the particle will "stick" to the largest column
        surface[i] = max(surface[left], surface[i] + 1, surface[right])
    end

    return surface
end

"""
Apply ballistic deposition to a 1D surface with periodic boundary conditions

# Arguments

- `surface::Vector{Integer}`: the 1D array representing the surface
- `time::Integer`: the number of time steps (or the number of deposited particles)

# Returns

boolean matrix containing the state of each cell (full / empty)
"""
function bdfull!(surface::Vector{<:Integer}, time::Integer)
    len = length(surface)
    starting_indecies = rand(1:len, time)
    grid = falses(time + maximum(surface), len)

    for i in starting_indecies
        # finding the left and right indecies with periodic boundary conditions
        left = i - 1 + ((len + 1 - i) ÷ len) * len
        right = i + 1 - (i ÷ len) * len

        # the particle will "stick" to the largest column
        surface[i] = max(surface[left], surface[i] + 1, surface[right])
        grid[surface[i], i] = true
    end

    return grid[1:(maximum(surface)), :]
end

"""
Apply ballistic deposition for an isolated "tree" of particles on a 1D surface

# Arguments

- `surface::Vector{Integer}`: The 1D array representing the surface. Must only contain one
    connected tree of particles.
- `time::Integer`: the number of time steps (or the number of deposited particles)

# Returns

the surface heights (Vector) after applying the ballistic deposition.
"""
function bdiso!(surface::Vector{<:Integer}, time::Integer)
    left, right = findleftright(surface)

    for _ in 1:time
        # The new particles can only be dropped on neighbors of the current points
        index = rand(left:right)

        # the particle will "stick" to the largest column
        surface[index] = max(surface[index - 1], surface[index] + 1, surface[index + 1])

        # expand the domain of dropping the particle
        if index == right
            right += 1
        elseif index == left
            left -= 1
        end
    end

    return surface
end

"""
Apply ballistic deposition for an isolated "tree" of particles on a 1D surface

# Arguments

- `surface::Vector{Integer}`: The 1D array representing the surface. Must only contain one
    connected tree of particles.
- `time::Integer`: the number of time steps (or the number of deposited particles)

# Returns

boolean matrix containing the state of each cell (full / empty)
"""
function bdisofull!(surface::Vector{<:Integer}, time::Integer)
    len = length(surface)
    grid = falses(time + maximum(surface), len)
    left, right = findleftright(surface)

    for _ in 1:time
        # The new particles can only be dropped on neighbors of the current points
        index = rand(left:right)

        # the particle will "stick" to the largest column
        surface[index] = max(surface[index - 1], surface[index] + 1, surface[index + 1])
        grid[surface[index], index] = true

        # expand the domain of dropping the particle
        if index == right
            right += 1
        elseif index == left
            left -= 1
        end
    end

    return grid[1:(maximum(surface)), :]
end

"""
Find left and right of the isolated tree in ballistic deposition (for bdiso functions)
"""
function findleftright(surface)
    leftindex = 1
    while surface[leftindex + 1] == 0
        leftindex += 1
    end

    rightindex = length(surface)
    while surface[rightindex - 1] == 0
        rightindex -= 1
    end

    return leftindex, rightindex
end

"""
Apply ballistic deposition with relaxation to a 1D surface with periodic boundary conditions

# Arguments

- `surface::Vector{Integer}`: the 1D array representing the surface
- `time::Integer`: the number of time steps (or the number of deposited particles)

# Returns

the surface array (vector) after applying the ballistic deposition
"""
function bdrelax!(surface::Vector{<:Integer}, time::Integer)
    len = length(surface)
    starting_indecies = rand(1:len, time)

    for i in starting_indecies
        # finding the left and right indecies with periodic boundary conditions
        left = i - 1 + ((len + 1 - i) ÷ len) * len
        right = i + 1 - (i ÷ len) * len

        # add one particle to the position with the minimum height
        # choose from left, i, and right
        if surface[i] > surface[left]
            if surface[left] > surface[right]
                surface[right] += 1
            else
                surface[left] += 1
            end
        elseif surface[i] > surface[right]
            surface[right] += 1
        else
            surface[i] += 1
        end
    end

    return surface
end

"""
Generate the heights at the given times for the given ballistic deposition

# Arguments

- `bdfunc!::Function`: a function that applies a ballistic deposition to a surface vector
    for a given time and returns the surface (it doesn't matter if the surface itself is
    modified or not)
- `len::Integer`: the length of the 1D surface
- `times::Vector{<:Integer}`: the times of each recorded height

# Returns

A matrix containing the array heights for the given times; Each row corresponds to one time
"""
function bdheights(bdfunc::Function, len::Integer, times::Vector{<:Integer})
    heights = Matrix{Integer}(undef, len, length(times))

    heights[:, 1] = bdfunc(zeros(Integer, len), times[1])
    for i in 2:length(times)
        heights[:, i] = bdfunc(heights[:, i - 1], times[i] - times[i - 1])
    end

    return heights
end

end

ext: jl

max_stars_repo_name: slhshamloo/comp-phys

max_issues_repo_name: slhshamloo/comp-phys

max_forks_repo_name: slhshamloo/comp-phys


text:
#!/usr/bin/env julia
using Gtk, Gtk.ShortNames, Printf
using Random: rand!

row, column, space = 400, 600, 20
mutable struct GOL
    tmpmat::Matrix{Bool}
    now::Matrix{Bool}
    row::Int
    column::Int

    function GOL()
        row_cell = div(row, space)
        column_cell = div(column, space)
        new(falses(row_cell, column_cell), falses(row_cell, column_cell), row_cell, column_cell)
    end
end
gol = GOL()

c = @GtkCanvas(column, row)
win = GtkWindow("Game of Life")
hbox = GtkBox(:h)
vbox = GtkBox(:v)
button_reset = GtkButton("Reset")
button_random = GtkButton("Random")
button_start = GtkButton("Start")
push!(vbox, button_reset)
push!(vbox, button_random)
push!(vbox, button_start)
push!(win, hbox)
push!(hbox, c)
push!(hbox, vbox)
set_gtk_property!(hbox, :expand, c, true)
showall(win)

"white brckgroud"
function makebroad!(c)
    @guarded draw(c) do widget
        ctx = getgc(c)
        rectangle(ctx, 0, 0, column, row)
        set_source_rgb(ctx, 1, 1, 1)
        fill(ctx)
    end
end

"make grid"
function grid!(c)
    @guarded draw(c) do widget
        n = 8
        ctx = getgc(c)
        h = height(c)
        w = width(c)

        for x in 0:space:w
            move_to(ctx, x, 0)
            line_to(ctx, x, h)
        end
        for y in 0:space:h
            move_to(ctx, 0, y)
            line_to(ctx, w, y)
        end
        set_source_rgb(ctx, 0, 0, 0)
        stroke(ctx)
        reveal(widget)
    end
    show(c)
end
# grid!(c)

function initilize!(c)
    makebroad!(c)
    grid!(c)
end
initilize!(c)

function paintcell!(c, x, y, red=0, green=1, blue=1)
    @guarded draw(c) do widget
        ctx = getgc(c)
        set_source_rgb(ctx, red, green, blue)
        rectangle(ctx, x*space, y*space, space, space)
        fill(ctx)
        reveal(widget)
    end
    # show(c)
end

function board_draw!(c, gol)
    row, column = gol.row, gol.column
    # grid!(c)
    for x in 0:column-1
        for y in 0:row-1
            if gol.now[y+1, x+1]
                paintcell!(c, x, y)
            else
                paintcell!(c, x, y, 1, 1, 1)
            end
        end
    end
    # grid!(c)
end

c.mouse.button1press = @guarded (widget, event) -> begin
    x, y = Int(div(event.x, space)), Int(div(event.y, space))
    println(x, " ", y)
    paintcell!(c, x, y)
    grid!(c)
    gol.now[y+1, x+1] = true
end

function board_reset!()
    initilize!(c)
    row, column = gol.row, gol.column
    fill!(gol.now, false)
    fill!(gol.tmpmat, false)
    # show(c)
    return
end
signal_connect(x -> board_reset!(), button_reset, "clicked")

function board_rand!()
    # initilize!(c)
    row, column = gol.row, gol.column
    rand!(gol.now)
    fill!(gol.tmpmat, false)
    board_draw!(c, gol)
    println("rand")
    # show(c)
    return
end
signal_connect(x -> board_rand!(), button_random, "clicked")

function next_gen_state(gol, x, y)
    row, column = gol.row, gol.column
    nowboard = gol.now
    num = 0
    for i in -1:1
        for j in -1:1
            num += nowboard[mod1(x+i, row), mod1(y+j, column)]
        end
    end
    # println(num)

    if !nowboard[x,y] # 死んでいた場合
        if num == 3
            return true
        else
            return false
        end
    else
        if num <= 2 || num >= 5
            return false
        else
            return true
        end
    end
end

function next_gen!(gol)
    row, column = gol.row, gol.column
    for x in 1:row
        for y in 1:column
            gol.tmpmat[x,y] = next_gen_state(gol, x, y)
        end
    end
    copyto!(gol.now, gol.tmpmat)
    return
end

ngen = 0
function board_start!()
    # initilize!(c)
    row, column = gol.row, gol.column
    for i in 1:50
        sleep(0.1)
        next_gen!(gol)
        # fill!(gol.tmpmat, false)
        board_draw!(c, gol)
        global ngen += 1
        @printf("%d\r", ngen)
        # show(c)
    end
    return
end
signal_connect(x -> board_start!(), button_start, "clicked")
showall(win)

if !isinteractive()
    cond = Condition()
    signal_connect(win, :destroy) do widget
        notify(cond)
    end
    wait(cond)
end

max_stars_repo_name: goropikari/JuliaGUIAppExample

max_issues_repo_name: goropikari/JuliaGUIAppExample

max_forks_repo_name: goropikari/JuliaGUIAppExample


text:
using DrWatson
quickactivate(@__DIR__)
@show projectdir()
##
using CovidSurvey
using CSV
using DataFrames
using Dates
# using ArviZ
using PrettyTables
import PyPlot as plt
# using RCall
using Turing
using StatsPlots
using Plots.PlotMeasures
using Random
using BSON
pgfplotsx()
## ============================================================================
# plotting utilities

function plot_effect_size(pp)
    c = group(pp.chain, "grouped_effect")
    effects = DataFrame(c)
    effects = select(effects, Not(["iteration","chain"]))
    rename!(effects, Dict(zip(Symbol.(names(effects)), Symbol.(pp.data.predictors)))...)

    p = plot(xlabel="effect size", legend=:outertopright)
    for (i,(v,n)) in enumerate(zip(eachcol(effects), names(effects)))
        label1 = i == 1 ? "95% HPDI" : nothing
        label2 = i == 1 ? "median" : nothing
        l,u = hpdi(v; alpha=0.05)
        m = median(v)
        plot!([l,u],[n,n], lw=2, c=:midnightblue, lab=label1)
        scatter!([m], [n], mc=:white, lab=label2, msw=2, lc=:black, ms=5)
    end
    return p
end

function last_latent_Rt(region, chain)
    i = Ns[region]
    c = group(chain, "latent_Rts")
    n = Symbol("latent_Rts[$region][$i]")
    j = findfirst(==(n), names(c))
    return Array(c[:,j,:]) |> vec
end

function sampled_effects(effect, region, chain)
    c = group(chain, "effects")
    n = Symbol("effects[$effect,$region]")
    j = findfirst(==(n), names(c))
    return Array(c[:,j,:]) |> vec
end

function sampled_effects(effect, chain)
    c = group(chain, "grouped_effect")
    n = Symbol("grouped_effect[$effect]")
    j = findfirst(==(n), names(c))
    return Array(c[:,j,:]) |> vec
end

# function sampled_effects(region, chain)
#     c = group(chain, "effects")
#     n = Symbol("effects[1,$region]")
#     j = findfirst(==(n), names(c))
#     return Array(c[:,j,:]) |> vec
# end

function relative_change_Rt(x, region, effect, pp)
    # r  = last_latent_Rt(region, pp.chain)
    # R0 = f.(r)
    α    = sampled_effects(effect, region, pp.chain)
    invf = pp.data.turing_data.invlink
    R0   = 1
    r    = invf(1.)
    f    = pp.data.turing_data.link
    RT   = f.(r .+ α * x)
    return @. ( RT / R0 - R0 ) * 100
end

function relative_change_Rt(x, effect, pp)
    # r  = last_latent_Rt(region, pp.chain)
    # R0 = f.(r)
    α    = sampled_effects(effect, pp.chain)
    invf = pp.data.turing_data.invlink
    f    = pp.data.turing_data.link
    R0   = 1
    r    = invf(1.)
    RT   = f.(r .+ α * x)
    return @. ( RT / R0 - R0 ) * 100
end

function plot_Rt_change(x, region, pp)
    predictors = pp.data.predictors
    chain = pp.chain
    title = Regional.regions[region]
    xlabel = region == 5 ? "reduction in Rt [\\%]" : ""
    legend = region == 1 ? :outertopright : nothing
    bottom_margin = region == 5 ? 0mm : 3mm
    top_margin = region == 1 ? 0mm : 3mm

    p = plot(; xlabel, legend, title, bottom_margin, top_margin)
    for (i,n) in enumerate(predictors)
        label1 = i == 1 ? "95% HPDI" : nothing
        label2 = i == 1 ? "median" : nothing
        v = relative_change_Rt(x, region, i, pp)
        l,u = hpdi(v; alpha=0.05)
        m = median(v)
        plot!([l,u],[n,n], lw=2, c=:midnightblue, lab=label1)
        scatter!([m], [n], mc=:white, lab=label2, msw=2, lc=:black, ms=5)
    end
    return p
end

function plot_Rt_change(x, pp)
    ps = Vector{Plots.Plot}()
    for region in 1:5
        p = plot_Rt_reduction(x, region, pp)
        push!(ps, p)
    end
    plot(ps..., layout=(5,1), size=(500, 1000), link=:x, sharex=true)
end

function plot_Rt_reduction_grouped_predictors(x, predictor, pp)
    predictors = pp.data.predictors
    num_predictors = length(predictors)
    title = predictors[predictor]
    chain = pp.chain
    xlabel = predictor == num_predictors ? "change in Rt [\\%]" : ""
    legend = predictor == 1 ? :outertopright : nothing
    bottom_margin = predictor == num_predictors ? 0mm : 3mm
    top_margin = predictor == 1 ? 0mm : 3mm

    p = plot(; xlabel, legend, title, bottom_margin, top_margin)
    for (i,n) in enumerate(Regional.regions)
        label1 = i == 1 ? "95% HPDI" : nothing
        label2 = i == 1 ? "median" : nothing
        v = relative_change_Rt(x, i, predictor, pp)
        l,u = hpdi(v; alpha=0.05)
        m = median(v)
        plot!([l,u],[n,n], lw=1, c=:midnightblue, lab=label1)
        scatter!([m], [n], mc=:white, lab=label2, msw=2, lc=:black, ms=5)
    end
    n = "grouped"
    v = relative_change_Rt(x, predictor, pp)
    l,u = hpdi(v; alpha=0.05)
    m = median(v)
    plot!([l,u],[n,n], lw=4, c=:midnightblue, lab=nothing)
    scatter!([m], [n], mc=:white, msw=3, lc=:black, ms=8, lab=nothing)
    return p
end

function plot_Rt_change_grouped_predictors(x, pp)
    ps = Vector{Plots.Plot}()
    predictors = pp.data.predictors
    N = length(predictors)
    for (i,n) in enumerate(predictors)
        p = plot_Rt_reduction_grouped_predictors(x, i, pp)
        push!(ps, p)
    end
    plot(ps..., layout=(N,1), size=(500, N*250), link=:x, sharex=true)
end

## =============================================================================
# load data
fname = projectdir("reports/regional_effect_estimation/specific_contacts/2021-08-11/CHAIN_chains=10_model=hospit_observ=2021-02-01_preds=CF,CC,CR,CS_regional=true_rwstep=7_std=false_steps=200_warmup=1500.jls")
# fname = projectdir("reports/regional_effect_estimation/specific_contacts/mixed_model/CHAIN_chains=10_model=hospit_observ=2021-02-01_preds=CF,CC,CR,CS_regional=true_rwstep=7_semipar=false_std=false_steps=200_warmup=1000.jls")
fdir, ps, ignores = Regional.parse_fname(fname)
pp = Regional.PostProcessing2(fdir, ps, ignores, [], fname)
pp = Regional.skip_warmup(pp)
# gp = Regional.generate_posterior(pp)
# fname = projectdir("reports/regional_effect_estimation/specific_contacts/mixed_model/GENERATED-QUANTITIES_chains=10_model=hospit_observ=2021-02-01_preds=CF,CC,CR,CS_regional=true_rwstep=7_semipar=false_std=false_steps=200_warmup=1000.bson")
fname = projectdir("reports/regional_effect_estimation/specific_contacts/2021-08-11/GENERATED-QUANTITIES_chains=10_model=hospit_observ=2021-02-01_preds=CF,CC,CR,CS_regional=true_rwstep=7_std=false_steps=200_warmup=1500.bson")
gp = BSON.load(fname) |> NamedTuple;

## =============================================================================
# plot fit result
default(titlefontsize = 20, legendfontsize = 18, labelfontsize = 18, guidefontsize = 18, tickfontsize = 12, framestyle = :zerolines, yminorgrid = true)
# dir = projectdir("figures/")

dir = projectdir("reports/regional_effect_estimation/specific_contacts/mixed_model/")
r = Regional.RegionsOverviewPlottingRecipe(
    pp.data,
    gp,
    :expected_daily_hospits
)
p = plot(r; plot_only_fit=true)
savefig(p, fdir*"overview_hospitalizations.png")
# savefig(p, fdir*"overview_hospitalizations.svg")

r = Regional.RegionsOverviewPlottingRecipe(
    pp.data,
    gp,
    :Rts
)
p = plot(r; plot_only_fit=true)
savefig(p, fdir*"overview_reproduction_number.png")
# savefig(p, fdir*"overview_reproduction_number.svg")

for region in Regional.regions
    r = Regional.RegionPlottingRecipe(pp.data, gp, region)
    p = plot(r; plot_only_fit=true)
    savefig(p, fdir*"overview_region=$(region).png")
    # savefig(p, fdir*"figures/overview_region=$(region).svg")
end

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

p = plot_effect_size(pp)
savefig(p, fdir*"effect_sizes.png")
savefig(p, fdir*"effect_sizes.svg")

x = -0.5
# p = plot_Rt_change(x, pp)
p = plot_Rt_change_grouped_predictors(x, pp)
savefig(p, fdir*"reproduction_number_half_contacts.png")
savefig(p, fdir*"reproduction_number_half_contacts.svg")

x = 2.
# p = plot_Rt_change(x, pp)
p = plot_Rt_change_grouped_predictors(x, pp)
savefig(p, fdir*"reproduction_number_double_contacts.png")
savefig(p, fdir*"reproduction_number_double_contacts.svg")

x = -0.5
d = []
for (i,p) in enumerate(pp.data.predictors)
    for (j,r) in enumerate(Regional.regions)
        v = relative_change_Rt(x, j, i, pp)
        l,u = hpdi(v; alpha=0.05)
        m = median(v)
        push!(d, [m, l, u, p, r])
    end
    r = "grouped"
    v = relative_change_Rt(x, i, pp)
    l,u = hpdi(v; alpha=0.05)
    m = median(v)
    push!(d, [m, l, u, p, r])
end
results = DataFrame(
    (median=m, lower=l, upper=u, predictor=p, region=r)
    for (m, l, u, p, r) in d
)
CSV.write(projectdir("reports/HOPE/effect_sizes_half_contacts.csv"), results)
##
default()
p = let
    c = group(pp.chain, "grouped_effect")
    df = DataFrame(c)
    df = select(df, Not(["iteration","chain"]))
    rename!(df, Dict(zip(Symbol.(names(df)), Symbol.(pp.data.predictors)))...)

    xs = df.colleagues
    p1 = histogram(xs, normalize=true, title="contacts to colleagues", label=nothing, fc=:gray, lw=.1, xlabel="effect size", bottom_margin=3mm)
    plot!(Exponential(0.2), xlims=(0,1), lw=2, label="prior", legend=:topright);
    vline!([mean(xs)], lw=2, label="mean = $(round(mean(xs), sigdigits=2))", c=:black, ls=:dash)
    xs = df.strangers
    p2 = histogram(xs, normalize=true, title="contacts to strangers", label=nothing, fc=:gray, lw=.1, xlabel="effect size", bottom_margin=3mm);
    plot!(Exponential(0.2), xlims=(0,1), lw=2, label="prior", legend=:topright);
    vline!([mean(xs)], lw=2, label="mean = $(round(mean(xs), sigdigits=2))", c=:black, ls=:dash)
    xs = df.family
    p3 = histogram(xs, normalize=true, title="contacts to family", label=nothing, fc=:gray, lw=.1, xlabel="effect size", top_margin=3mm);
    plot!(Exponential(0.2), xlims=(0,1), lw=2, label="prior", legend=:topright);
    vline!([mean(xs)], lw=2, label="mean = $(round(mean(xs), sigdigits=2))", c=:black, ls=:dash)
    xs = df.friends
    p4 = histogram(xs, normalize=true, title="contacts to friends", label=nothing, fc=:gray, lw=.1, xlabel="effect size", top_margin=3mm);
    plot!(Exponential(0.2), xlims=(0,1), lw=2, label="prior", legend=:topright);
    vline!([mean(xs)], lw=2, label="mean = $(round(mean(xs), sigdigits=2))", c=:black, ls=:dash)
    plot(p1,p2,p3,p4, layout=(2,2))
end

savefig(p, projectdir("figures/prior_vs_posterior_effect.png"))
## ==========================================================================
#                         TOTAL CONTACTS
#  ==========================================================================
# load data
fname = projectdir("reports/regional_effect_estimation/total_contacts/2021-08-10/CHAIN_chains=10_model=hospit_observ=2021-02-01_preds=CT_regional=true_rwstep=7_std=false_steps=200_warmup=1000.jls")
fdir, ps, ignores = Regional.parse_fname(fname)
pp = Regional.PostProcessing2(fdir, ps, ignores, [], fname)
pp = Regional.skip_warmup(pp)

fname = projectdir("reports/regional_effect_estimation/total_contacts/2021-08-10/GENERATED-QUANTITIES_chains=10_model=hospit_observ=2021-02-01_preds=CT_regional=true_rwstep=7_std=false_steps=200_warmup=1000.bson")
gp = BSON.load(fname) |> NamedTuple;
## =============================================================================
# plot fit result
default(titlefontsize = 20, legendfontsize = 18, labelfontsize = 18, guidefontsize = 18, tickfontsize = 12, framestyle = :zerolines, yminorgrid = true)

r = Regional.RegionsOverviewPlottingRecipe(
    pp.data,
    gp,
    :expected_daily_hospits
)
p = plot(r; plot_only_fit=true)
savefig(p, projectdir("figures/overview_hospitalizations_TOTAL.png"))
savefig(p, projectdir("figures/overview_hospitalizations_TOTAL.svg"))

r = Regional.RegionsOverviewPlottingRecipe(
    pp.data,
    gp,
    :Rts
)
p = plot(r; plot_only_fit=true)
savefig(p, projectdir("figures/overview_reproduction_number_TOTAL.png"))
savefig(p, projectdir("figures/overview_reproduction_number_TOTAL.svg"))

for region in Regional.regions
    r = Regional.RegionPlottingRecipe(pp.data, gp, region)
    p = plot(r; plot_only_fit=true)
    savefig(p, projectdir("figures/overview_region=$(region)_TOTAL.png"))
end

## ============================================================================
# plot effects

x = -0.5
p = plot_Rt_change_grouped_predictors(x, pp)
savefig(p, projectdir("figures/reproduction_number_half_contacts_TOTAL.png"))
savefig(p, projectdir("figures/reproduction_number_half_contacts_TOTAL.svg"))

x = 2.
p = plot_Rt_change_grouped_predictors(x, pp)
savefig(p, projectdir("figures/reproduction_number_double_contacts_TOTAL.png"))
savefig(p, projectdir("figures/reproduction_number_double_contacts_TOTAL.svg"))

x = 2#-0.5
d = []
for (i,p) in enumerate(pp.data.predictors)
    for (j,r) in enumerate(Regional.regions)
        v = relative_change_Rt(x, j, i, pp)
        l,u = hpdi(v; alpha=0.05)
        m = median(v)
        push!(d, [m, l, u, p, r])
    end
    r = "grouped"
    v = relative_change_Rt(x, i, pp)
    l,u = hpdi(v; alpha=0.05)
    m = median(v)
    push!(d, [m, l, u, p, r])
end
results = DataFrame(
    (median=m, lower=l, upper=u, predictor=p, region=r)
    for (m, l, u, p, r) in d
)

# CSV.write(projectdir("reports/HOPE/effect_sizes_half_contacts_TOTAL.csv"), results)
##
default()
p = let
    xs = df.colleagues
    p1 = histogram(xs, normalize=true, title="contacts to colleagues", label=nothing, fc=:gray, lw=.1, xlabel="effect size", bottom_margin=3mm);
    plot!(Exponential(0.2), xlims=(0,1), lw=2, label="prior", legend=:topright);
    vline!([mean(xs)], lw=2, label="mean = $(round(mean(xs), sigdigits=2))", c=:black, ls=:dash)
    xs = df.strangers
    p2 = histogram(xs, normalize=true, title="contacts to strangers", label=nothing, fc=:gray, lw=.1, xlabel="effect size", bottom_margin=3mm);
    plot!(Exponential(0.2), xlims=(0,1), lw=2, label="prior", legend=:topright);
    vline!([mean(xs)], lw=2, label="mean = $(round(mean(xs), sigdigits=2))", c=:black, ls=:dash)
    xs = df.family
    p3 = histogram(xs, normalize=true, title="contacts to family", label=nothing, fc=:gray, lw=.1, xlabel="effect size", top_margin=3mm);
    plot!(Exponential(0.2), xlims=(0,1), lw=2, label="prior", legend=:topright);
    vline!([mean(xs)], lw=2, label="mean = $(round(mean(xs), sigdigits=2))", c=:black, ls=:dash)
    xs = df.friends
    p4 = histogram(xs, normalize=true, title="contacts to friends", label=nothing, fc=:gray, lw=.1, xlabel="effect size", top_margin=3mm);
    plot!(Exponential(0.2), xlims=(0,1), lw=2, label="prior", legend=:topright);
    vline!([mean(xs)], lw=2, label="mean = $(round(mean(xs), sigdigits=2))", c=:black, ls=:dash)
    plot(p1,p2,p3,p4, layout=(2,2))
end

savefig(p, projectdir("figures/prior_vs_posterior_effect.png"))

max_stars_repo_name: andreaskoher/CovidSurvey

max_issues_repo_name: andreaskoher/CovidSurvey

max_forks_repo_name: andreaskoher/CovidSurvey


text:
#inspired by nucleus theme for Atom
const dark_palette = [ colorant"#FE4365", # red
orange= colorant"#eca25c", # orange
green = colorant"#3f9778", # green
blue = colorant"#005D7F"] # blue

dark_bg = colorant"#363D46"
dark = PlotTheme(dark_bg,
                    colorant"#30343B",
                    colorant"#ADB2B7",
                    colorant"#FFFFFF",
                    expand_palette(dark_bg, dark_palette; lchoices=linspace(57,57,1),
                                                          cchoices=linspace(100,100,1)),
                    dark_palette[[2,1]])

add_plots_theme(:dark, dark)


text:
module RocketVoidActorTest

using Test
using Rocket

@testset "VoidActor" begin

    println("Testing: actor VoidActor")

    @testset begin
        actor = VoidActor{Int}()

        @test_throws InconsistentSourceActorDataTypesError{Int64,String} next!(actor, "string")
    end

    @testset begin
        factory = void()

        @test Rocket.create_actor(Int, factory) isa VoidActor{Int}
    end

    @testset begin
        @test void(Int) isa VoidActor{Int}
        @test void()    isa Rocket.VoidActorFactory
    end
end

end


text:
refVals = []
refPrecs = []

#! format: off
# SC ========== Test number 1 reference values and precision match template. =======
# SC ========== /home/jmc/cliMa/cliMa_update/test/Ocean/SplitExplicit/simple_box_2dt.jl test reference values ======================================
# BEGIN SCPRINT
# varr - reference values (from reference run)
# parr - digits match precision (hand edit as needed)
#
# [
#  [ MPIStateArray Name, Field Name, Maximum, Minimum, Mean, Standard Deviation ],
#  [         :                :          :        :      :          :           ],
# ]
varr = [
 [  "oce Q_3D",   "u[1]", -1.50801004298594193420e-01,  1.72156284168757833797e-01,  4.68279918998404202923e-03,  1.95837597627477114082e-02 ],
 [  "oce Q_3D",   "u[2]", -1.41680528424183271863e-01,  1.46287507569445834887e-01, -2.51414399339626839017e-03,  1.75694496745556552608e-02 ],
 [  "oce Q_3D",       :η, -7.34574736037479247841e-01,  3.69081182348461100950e-01, -2.03558200029154242510e-04,  2.71392010992079990572e-01 ],
 [  "oce Q_3D",       :θ,  6.65420328206461388076e-03,  9.92226603474087376355e+00,  2.49938050028639668909e+00,  2.17889266711590456893e+00 ],
 [   "oce aux",       :w, -1.93638744994098439866e-04,  1.64381091769879710560e-04,  4.21611820820242154596e-07,  1.43969264053567314515e-05 ],
 [   "oce aux",    :pkin, -8.85156811621749994856e+00,  0.00000000000000000000e+00, -3.26783582829862284669e+00,  2.49771056538873326147e+00 ],
 [   "oce aux",     :wz0, -2.13040831337660739543e-05,  3.03801724923109659664e-05, -1.30760986800648726445e-10,  7.79224859585881761404e-06 ],
 [   "oce aux", "u_d[1]", -1.47512664164697798741e-01,  1.15209893845465857765e-01, -3.40592527368550041888e-05,  1.18072942147697344112e-02 ],
 [   "oce aux", "u_d[2]", -1.41197405478766774278e-01,  1.40391286353623745420e-01, -8.62962475112239205902e-06,  1.16484682883619180410e-02 ],
 [   "oce aux", "ΔGu[1]", -1.83391366147010946068e-06,  1.75788411864240537568e-06, -3.56389896621704670505e-08,  2.63843686150703932422e-07 ],
 [   "oce aux", "ΔGu[2]", -1.20700503402749584594e-06,  2.22174863758920371976e-06,  1.28960840828931523409e-06,  6.84169299520288045132e-07 ],
 [   "oce aux",       :y,  0.00000000000000000000e+00,  4.00000000000000046566e+06,  2.00000000000000000000e+06,  1.15573163901915703900e+06 ],
 [ "baro Q_2D",   "U[1]", -1.28725083574144978371e+01,  6.09709288942672031908e+01,  4.70927375834381667374e+00,  1.49705019656797517769e+01 ],
 [ "baro Q_2D",   "U[2]", -2.82713167453262848028e+01,  6.48523258594305360702e+01, -2.47826878693387353536e+00,  1.34485730829431311406e+01 ],
 [ "baro Q_2D",       :η, -7.34779137552936445488e-01,  3.69122464170923092475e-01, -2.03588902795195280044e-04,  2.71542408145441560574e-01 ],
 [  "baro aux",  "Gᵁ[1]", -1.75788411864240528421e-03,  1.83391366147010953860e-03,  3.56389896621704717621e-05,  2.63856747396035119183e-04 ],
 [  "baro aux",  "Gᵁ[2]", -2.22174863758920360796e-03,  1.20700503402749580867e-03, -1.28960840828931536453e-03,  6.84203168449271982056e-04 ],
 [  "baro aux", "U_c[1]", -1.28505196446931666543e+01,  6.09003088503453255953e+01,  4.71562834799378549633e+00,  1.49564184259850492253e+01 ],
 [  "baro aux", "U_c[2]", -2.82949777356440890230e+01,  6.48066028531587363659e+01, -2.50541914991992387129e+00,  1.34465940520910596234e+01 ],
 [  "baro aux",     :η_c, -7.34574736037479247841e-01,  3.69081182348461100950e-01, -2.03558200029152670417e-04,  2.71405445907609677025e-01 ],
 [  "baro aux", "U_s[1]", -1.28507611526714438099e+01,  6.09013504831384295812e+01,  4.71564903049480754760e+00,  1.49566179983317635305e+01 ],
 [  "baro aux", "U_s[2]", -2.82958333614026145142e+01,  6.48064580265526331004e+01, -2.50548479627681608406e+00,  1.34466238247583795840e+01 ],
 [  "baro aux",     :η_s, -7.34582576718505908531e-01,  3.69085457622194967886e-01, -2.03553625318579639201e-04,  2.71406256014543156851e-01 ],
 [  "baro aux",  "Δu[1]", -4.52663135359685495692e-04,  4.01492679065675853092e-04, -1.17669158892409481397e-05,  1.01693312716598876057e-04 ],
 [  "baro aux",  "Δu[2]", -2.43836577383847778961e-04,  3.86419757650743997315e-04,  5.11479284073749210704e-05,  7.52351009272865279024e-05 ],
 [  "baro aux",  :η_diag, -7.32001048379119167997e-01,  3.67552372776876246085e-01, -2.03164560876945184120e-04,  2.71190595495126973713e-01 ],
 [  "baro aux",      :Δη, -3.56517726875715812440e-03,  4.27489668080521800775e-03, -3.93639152206426446375e-07,  1.12635709821496884930e-03 ],
 [  "baro aux",       :y,  0.00000000000000000000e+00,  4.00000000000000046566e+06,  2.00000000000000000000e+06,  1.15578885204060329124e+06 ],
]
parr = [
 [  "oce Q_3D",   "u[1]",    12,    12,    12,    12 ],
 [  "oce Q_3D",   "u[2]",    12,    12,    12,    12 ],
 [  "oce Q_3D",       :η,    12,    12,     8,    12 ],
 [  "oce Q_3D",       :θ,    12,    12,    12,    12 ],
 [   "oce aux",       :w,    12,    12,     8,    12 ],
 [   "oce aux",    :pkin,    12,    12,    12,    12 ],
 [   "oce aux",     :wz0,    12,    12,     8,    12 ],
 [   "oce aux", "u_d[1]",    12,    12,    12,    12 ],
 [   "oce aux", "u_d[2]",    12,    12,    12,    12 ],
 [   "oce aux", "ΔGu[1]",    12,    12,    12,    12 ],
 [   "oce aux", "ΔGu[2]",    12,    12,    12,    12 ],
 [   "oce aux",       :y,    12,    12,    12,    12 ],
 [ "baro Q_2D",   "U[1]",    12,    12,    12,    12 ],
 [ "baro Q_2D",   "U[2]",    12,    12,    12,    12 ],
 [ "baro Q_2D",       :η,    12,    12,     8,    12 ],
 [  "baro aux",  "Gᵁ[1]",    12,    12,    12,    12 ],
 [  "baro aux",  "Gᵁ[2]",    12,    12,    12,    12 ],
 [  "baro aux", "U_c[1]",    12,    12,    12,    12 ],
 [  "baro aux", "U_c[2]",    12,    12,    12,    12 ],
 [  "baro aux",     :η_c,    12,    12,     8,    12 ],
 [  "baro aux", "U_s[1]",    12,    12,    12,    12 ],
 [  "baro aux", "U_s[2]",    12,    12,    12,    12 ],
 [  "baro aux",     :η_s,    12,    12,     8,    12 ],
 [  "baro aux",  "Δu[1]",    12,    12,    12,    12 ],
 [  "baro aux",  "Δu[2]",    12,    12,    12,    12 ],
 [  "baro aux",  :η_diag,    12,    12,     8,    12 ],
 [  "baro aux",      :Δη,    12,    12,     8,    12 ],
 [  "baro aux",       :y,    12,    12,    12,    12 ],
]
# END SCPRINT
# SC ====================================================================================

    append!(refVals ,[ varr ] )
    append!(refPrecs,[ parr ] )

#! format: on

max_stars_repo_name: mwarusz/CLIMA

max_issues_repo_name: mwarusz/CLIMA

max_forks_repo_name: mwarusz/CLIMA


text:
using PyCall
using Distributions

@pyimport matplotlib.pyplot as plt
@pyimport matplotlib.figure as pyfigure

function part1()

  un = Uniform(-2,2)

  samples = rand(un,500)
  x = 1:500
  oneoverx = map(elem -> 1/elem,x)

  plt.figure(figsize=(20,10))

  # unif
  plt.xlim(xmin=-50,xmax=550)
  plt.ylim(ymin=-5,ymax=5)
  plt.grid()
  plt.plot(x,samples,color="red")

  # sin(x)
  plt.xlim(xmin=-50,xmax=550)
  plt.ylim(ymin=-5,ymax=5)
  plt.grid()
  plt.plot(1:500,sin(x),color="blue")

  # cos(x)
  plt.xlim(xmin=-50,xmax=550)
  plt.ylim(ymin=-5,ymax=5)
  plt.grid()
  plt.plot(1:500,cos(x),color="green")

  # sin(1/x)
  plt.xlim(xmin=-50,xmax=550)
  plt.ylim(ymin=-5,ymax=5)
  plt.grid()
  plt.plot(1:500,sin(oneoverx),color="black")

  # cos(1/x)
  plt.xlim(xmin=-50,xmax=550)
  plt.ylim(ymin=-5,ymax=5)
  plt.grid()
  plt.plot(1:500,cos(oneoverx),color="pink")

  plt.show()

end

function part2()
  un = Uniform(-2,2)

  samples = rand(un,500)
  x = 1:500
  oneoverx = map(elem -> 1/elem,x)

  plt.figure(figsize=(20,10))

  # unif
  plt.subplot(331)
  plt.xlim(xmin=-50,xmax=550)
  plt.ylim(ymin=-5,ymax=5)
  plt.grid()
  plt.plot(x,samples,color="red")

  # sin(x)
  plt.subplot(332)
  plt.xlim(xmin=-50,xmax=550)
  plt.ylim(ymin=-5,ymax=5)
  plt.grid()
  plt.plot(1:500,sin(x),color="blue")

  # cos(x)
  plt.subplot(333)
  plt.xlim(xmin=-50,xmax=550)
  plt.ylim(ymin=-5,ymax=5)
  plt.grid()
  plt.plot(1:500,cos(x),color="green")

  # sin(1/x)
  plt.subplot(334)
  plt.xlim(xmin=-50,xmax=550)
  plt.ylim(ymin=-5,ymax=5)
  plt.grid()
  plt.plot(1:500,sin(oneoverx),color="black")

  # cos(1/x)
  plt.subplot(335)
  plt.xlim(xmin=-50,xmax=550)
  plt.ylim(ymin=-5,ymax=5)
  plt.grid()
  plt.plot(1:500,cos(oneoverx),color="pink")

  plt.show()

end

part2()


text:
#
#  File:    cqo1.jl
#
#  Purpose: Demonstrates how to solve small conic
#           optimization problem using the MOSEK Python API.
##

using Mosek
using Printf

printstream(msg::AbstractString) = print(msg)
callback(where,dinf,iinf,liinf) = 0

# Since the actual value of Infinity is ignores, we define it solely
# for symbolic purposes:

bkc = [ MSK_BK_FX ]
blc = [ 1.0 ]
buc = [ 1.0 ]

c   = [               0.0,              0.0,              0.0,
                      1.0,              1.0,              1.0 ]
bkx = [ MSK_BK_LO,MSK_BK_LO,MSK_BK_LO,
        MSK_BK_FR,MSK_BK_FR,MSK_BK_FR ]
blx = [               0.0,              0.0,              0.0,
                     -Inf,             -Inf,             -Inf ]
bux = [               Inf,              Inf,              Inf,
                      Inf,              Inf,              Inf ]

asub  = [ 1 ,2, 3 ]
aval  = [ 1.0, 1.0, 1.0 ]

numvar = length(bkx)
numcon = length(bkc)

# Create a task
maketask() do task
    putstreamfunc(task,MSK_STREAM_LOG,printstream)
    putcallbackfunc(task,callback)

    # Append 'numcon' empty constraints.
    # The constraints will initially have no bounds.
    appendcons(task,numcon)

    #Append 'numvar' variables.
    # The variables will initially be fixed at zero (x=0).
    appendvars(task,numvar)

    # Set the linear term c_j in the objective.
    putclist(task,[1:6;],c)

    # Set the bounds on variable j
    # blx[j] <= x_j <= bux[j]
    putvarboundslice(task,1,numvar+1,bkx,blx,bux)

    putarow(task,1,asub,aval)
    putconbound(task,1,bkc[1],blc[1],buc[1])

    # Input the cones
    appendcone(task,MSK_CT_QUAD, 0.0, [ 4, 1, 2 ])
    appendcone(task,MSK_CT_RQUAD, 0.0, [ 5, 6, 3 ])

    # Input the objective sense (minimize/maximize)
    putobjsense(task,MSK_OBJECTIVE_SENSE_MINIMIZE)

    # Optimize the task
    optimize(task,"mosek://solve.mosek.com:30080")
    # Print a summary containing information
    # about the solution for debugging purposes
    solutionsummary(task,MSK_STREAM_MSG)
    prosta = getprosta(task,MSK_SOL_ITR)
    solsta = getsolsta(task,MSK_SOL_ITR)

    if solsta == MSK_SOL_STA_OPTIMAL
        # Output a solution
        xx = getxx(task,MSK_SOL_ITR)
        @printf("Optimal solution: %s\n", xx')
    elseif solsta == MSK_SOL_STA_DUAL_INFEAS_CER
        println("Primal or dual infeasibility.\n")
    elseif solsta == MSK_SOL_STA_PRIM_INFEAS_CER
        println("Primal or dual infeasibility.\n")
    elseif solsta == MSK_SOL_STA_UNKNOWN
        println("Unknown solution status")
    else
        println("Other solution status")
    end
end


text:
mutable struct TriangulateInputMapper
  originalVerticesMarkers::Vector{Cint}
  remappedVerticesMarkers::Vector{Cint}
  originalEdgesList::Vector{Cint}
  remappedEdgesList::Vector{Cint}
  TriangulateInputMapper() = new(Vector{Cint}(), Vector{Cint}(), Vector{Cint}(), Vector{Cint}())
end

function trimap_to_native(verticesMap::Vector{Cint}, edges::Vector{Cint} = Vector{Cint}())
    nativeMap = TriangulateInputMapper()

    nativeMap.originalVerticesMarkers = verticesMap
    nativeMap.remappedVerticesMarkers = collect(1:length(verticesMap))
    dictVertices = Dict{Int64,Int64}()
    for (idx,el) in enumerate(verticesMap)
        dictVertices[el] = idx
    end

    if length(edges) > 0
        nativeMap.originalEdgesList = edges
        nativeMap.remappedEdgesList = Vector{Cint}()
        for el in edges
            push!(nativeMap.remappedEdgesList, dictVertices[el])
        end
    end

    return nativeMap
end

function trimap_from_native(nativeMap::TriangulateInputMapper, trianglelist::Vector{Cint})
    newtrianglelist = Vector{Cint}(undef, length(trianglelist))

    for (idx,el) in enumerate(trianglelist)
        newtrianglelist[idx] = nativeMap.originalVerticesMarkers[el]
    end

    return newtrianglelist
end


text:
module TypeExpression
##=========common head===============
using HMRowUnification
using Sequent
using MLStyle

ImD = Base.ImmutableDict
struct Store{K, O}
    unbox :: ImD{K, O}
end

# given σ ∈ Store{K, O} where K, O
# 'σ[k => v]' add a scoped pair k=>v and returns a new store
# where the old one is not affected.
function Store{K, O}() where {K, O}
    Store(ImD{K, O}())
end
Base.getindex(s::Store{K, O}, x::Pair{K, A}...) where {K, O, A<:O} =
    let im = s.unbox
        for each in x
            im = ImD(im, each)
        end
        Store(im)
    end
Base.getindex(s::Store{K, O}, x::K) where {K, O} = s.unbox[x]
##===========auxiliaries==================

@active T(x) begin
    @match x begin
        App(Nom(:Type), x) => Some(x)
        _ => nothing
    end
end

(::Type{T})(x::HMT) = App(Nom(:Type), x)

(f::HMT)(arg::HMT) = App(f, arg)
(a::HMT) ↦ (r::HMT) = Arrow(a, r)

σ₀ = Store{Symbol, HMT}()[
    :int => T(Nom(:int)),
    :str => T(Nom(:str))
]

##=============semantics===================
tr_type = @sequent TType begin
    σ ⊢ₜ t => t′
    σ ⊢ₜ :(Type[$t]) => T(t′)
end

tr_app = @sequent TAPP begin
    σ ⊢ₜ op => T(op_t)
    σ ⊢ₜ args .=> [T(arg_t) for arg_t in args_t]

    σ ⊢ₜ :($op[$(args...)]) =>
        length(args) == 1 ?
        T(op_t(args_t[1])) :
        let args_t = Tuple(args_t)
            T(op_t(Tup(args_t)))
        end
end

tr_tup = @sequent TTUP begin
    σ ⊢ₜ args .=> [T(arg_t) for arg_t in args_t]
    σ ⊢ₜ Expr(:tuple, args...) =>
        let args_t = Tuple(args_t)
            T(Tup(args_t))
        end
end

tr_sym = @sequent TSYM begin
    σ ⊢ₜ a::Symbol => σ[a]
end

tr_arrow = @sequent TARROW begin
    σ ⊢ₜ f => T(ft)
    σ ⊢ₜ arg => T(argt)
    σ ⊢ₜ :($f ↦ $arg) => T(ft ↦ argt)
end

tr_forall = @sequent TFORALL begin
    bounds = Pair{Symbol, HMT}[s=>T(Fresh(s)) for s in syms]
    σ[bounds...] ⊢ₜ t => T(t′)
    σ ⊢ₜ :(forall[$([s::Symbol for s in syms]...)][$t]) =>
        let freshvars = Tuple(syms)
            T(Forall(freshvars, t′))
        end
end

##=============codegen for rules==========
# the order of rules are always sigficant, because given
#
#      Judge1
#      Judge2
#   ________________
#     Ctx ⊢ Input => Output
#
# we will do pattern matching for (Ctx, Input).
# Only the first rule matching in the given order applies.
# For instance, `tr_forall` should be prior than `tr_app`,
# as the pattern
#    :($op[$(args...)])
# covers
#    :(forall[$([s::Symbol for s in syms]...)][$t])

@semantics [
    tr_forall,
    tr_type,
    tr_app,
    tr_tup,
    tr_sym,
    tr_arrow
]
##==========================================
# after codegen, we have ⊢ₜ

println(σ₀ ⊢ₜ :(forall[a, b][a ↦ (int, b)]))
# Type (forall a b.a->{int, b})

end # module


text:
using LightGraphs
using SimpleWeightedGraphs
using Test

## TODO tests

DIGRAPHS = Dict{String,DiGraph}(
    "complete20" => complete_digraph(20),
    "path40" => path_digraph(40)
)
GRAPHS = Dict{String,Graph}(
    "petersen" => smallgraph(:petersen),
    "tutte" => smallgraph(:tutte),
    "complete20" => complete_graph(20)
)
WEIGHTED_GRAPHS = Dict{String,SimpleWeightedGraph}(
    ""
)

tests = [
    "operators.jl",
    "diversity.jl",
    "sets.jl",
    "fractional_numbers.jl"
]

@testset GraphAlgorithms begin
    for t in tests
        include(t)
    end
end


text:
"""
$SIGNATURES

Take a page (in HTML) and check all `href` on it to see if they lead somewhere.
"""
function verify_links_page(path::AS, online::Bool)
    shortpath = replace(path, PATHS[:folder] => "")
    mdpath    = replace(shortpath, Regex("^pub$(escape_string(PATH_SEP))")=>"pages$(PATH_SEP)")
    mdpath    = splitext(mdpath)[1] * ".md"
    shortpath = replace(shortpath, r"^\/"=>"")
    mdpath    = replace(mdpath, r"^\/"=>"")
    allok     = true
    page      = read(path, String)
    for m in eachmatch(r"\shref=\"(https?://)?(.*?)(?:#(.*?))?\"", page)
        if m.captures[1] === nothing
            # internal link, remove the front / otherwise it fails with joinpath
            link   = replace(m.captures[2], r"^\/"=>"")
            # if it's empty it's `href="/"` which is always ok
            isempty(link) && continue
            anchor = m.captures[3]
            full_link = joinpath(PATHS[:folder], link)
            if !isfile(full_link)
                allok && println("")
                println("- internal link issue on page $mdpath: $link.")
                allok = false
            else
                if !isnothing(anchor)
                    if full_link != path
                        rpage = read(full_link, String)
                    else
                        rpage = page
                    end
                    # look for `id=...` either with or without quotation marks
                    if match(Regex("id=(?:\")?$anchor(?:\")?"), rpage) === nothing
                        allok && println("")
                        println("- internal link issue on page $mdpath: $link.")
                        allok = false
                    end
                end
            end
        else
            online || continue
            # external link
            link = m.captures[1] * m.captures[2]
            try
                HTTP.request("HEAD", link).status == 200 || throw()
            catch
                allok && println("")
                println("- external link issue on page $mdpath: $link")
                allok = false
            end
        end
    end
    return allok
end

"""
$SIGNATURES

Verify all links in generated HTML.
"""
function verify_links()::Nothing
    # check that the user is online (otherwise only verify internal links)
    online = HTTP.request("HEAD", "https://discourse.julialang.org/").status == 200

    print("Verifying links...")
    if online
        print(" [you seem online ✓]")
    else
        print(" [you don't seem online ✗]")
    end

    # go over `index.html` then everything in `pub/`
    overallok = verify_links_page(joinpath(PATHS[:folder], "index.html"), online)

    for (root, _, files) ∈ walkdir(PATHS[:pub])
        for file ∈ files
            splitext(file)[2] == ".html" && continue
            path = joinpath(root, file)
            allok = verify_links_page(path)
            overallok = overallok && allok
        end
    end
    overallok && println("\rAll internal $(online && "and external ")links verified ✓.      ")
    return nothing
end

const JD_PY_MIN_NAME = ".__py_tmp_minscript.py"

"""
$(SIGNATURES)

Does a full pass followed by a pre-rendering and minification step.

* `prerender=true`: whether to pre-render katex and highlight.js (requires `node.js`)
* `minify=true`:    whether to minify output (requires `python3` and `css_html_js_minify`)
* `sig=false`:      whether to return an integer indicating success (see [`publish`](@ref))
* `prepath=""`:     set this to something like "project-name" if it's a project page

Note: if the prerendering is set to `true`, the minification will take longer as the HTML files
will be larger (especially if you have lots of maths on pages).
"""
function optimize(; prerender::Bool=true, minify::Bool=true, sig::Bool=false,
                    prepath::String="", no_fail_prerender::Bool=true,
                    suppress_errors::Bool=true)::Union{Nothing,Bool}
    suppress_errors && (SUPPRESS_ERR[] = true)
    #
    # Prerendering
    #
    if prerender && !JD_CAN_PRERENDER
        @warn "I couldn't find node and so will not be able to pre-render javascript."
        prerender = false
    end
    if prerender && !JD_CAN_HIGHLIGHT
        @warn "I couldn't load 'highlight.js' so will not be able to pre-render code blocks. " *
              "You can install it with `npm install highlight.js`."
    end
    if !isempty(prepath)
        GLOBAL_PAGE_VARS["prepath"] = prepath => (String,)
    end
    # re-do a (silent) full pass
    start = time()
    fmsg = "\r→ Full pass"
    print(fmsg)
    withpre = fmsg * ifelse(prerender, rpad(" (with pre-rendering)", 24), rpad(" (no pre-rendering)", 24))
    print(withpre)
    succ = (serve(single=true, prerender=prerender, nomess=true, isoptim=true, no_fail_prerender=no_fail_prerender) === nothing)
    print_final(withpre, start)

    #
    # Minification
    #
    if minify && (succ || no_fail_prerender)
        if JD_CAN_MINIFY
            start = time()
            mmsg = rpad("→ Minifying *.[html|css] files...", 35)
            print(mmsg)
            # copy the script to the current dir
            cp(joinpath(dirname(pathof(JuDoc)), "scripts", "minify.py"), JD_PY_MIN_NAME; force=true)
            # run it
            succ = success(`$([e for e in split(PY)]) $JD_PY_MIN_NAME`)
            # remove the script file
            rm(JD_PY_MIN_NAME)
            print_final(mmsg, start)
        else
            @warn "I didn't find css_html_js_minify, you can install it via pip."*
                  "The output will not be minified."
        end
    end
    SUPPRESS_ERR[] = false
    return ifelse(sig, succ, nothing)
end

"""
$(SIGNATURES)

This is a simple wrapper doing a git commit and git push without much fanciness. It assumes the
current directory is a git folder.
It also fixes all links if you specify `prepath` (or if it's set in `config.md`).

**Keyword arguments**

* `prerender=true`: prerender javascript before pushing see [`optimize`](@ref)
* `minify=true`:    minify output before pushing see [`optimize`](@ref)
* `nopass=false`:   set this to true if you have already run `optimize` manually.
* `prepath=""`:     set this to something like "project-name" if it's a project page
* `message="jd-update"`: add commit message.
"""
function publish(; prerender::Bool=true, minify::Bool=true, nopass::Bool=false,
                   prepath::String="", message::String="jd-update")::Nothing
    succ = true
    if !isempty(prepath) || !nopass
        succ = optimize(prerender=prerender, minify=minify, sig=true, prepath=prepath)
    end
    if succ
        start = time()
        pubmsg = rpad("→ Pushing updates with git...", 35)
        print(pubmsg)
        try
            run(`git add -A `)
            run(`git commit -m "$message" --quiet`)
            run(`git push --quiet`)
            print_final(pubmsg, start)
        catch e
            println("✘ Could not push updates, verify your connection and try manually.\n")
            @show e
        end
    else
        println("✘ Something went wrong in the optimisation step. Not pushing updates.")
    end
    return nothing
end

"""
$(SIGNATURES)

Cleanpull allows you to pull from your remote git repository after having removed the local
output directory. This will help avoid merge clashes.
"""
function cleanpull()::Nothing
    FOLDER_PATH[] = pwd()
    set_paths!()
    if isdir(PATHS[:pub])
        rmmsg = rpad("→ Removing local output dir...", 35)
        print(rmmsg)
        rm(PATHS[:pub], force=true, recursive=true)
        println("\r" * rmmsg * " [done ✔ ]")
    end
    try
        pmsg = rpad("→ Retrieving updates from the repository...", 35)
        print(pmsg)
        run(`git pull --quiet`)
        println("\r" * pmsg * " [done ✔ ]")
    catch e
        println("✘ Could not pull updates, verify your connection and try manually.\n")
        @show e
    end
    return nothing
end


text:
@testset "import" begin
    sys = pyimport("sys")
    os = pyimport("os")
    @test pyeq(Bool, sys.__name__, "sys")
    @test pyeq(Bool, os.__name__, "os")
    sysos = pyimport("sys", "os")
    @test sysos isa Tuple{Py, Py}
    @test pyis(sysos[1], sys)
    @test pyis(sysos[2], os)
    ver = pyimport("sys" => "version")
    @test pyis(ver, sys.version)
    path = pyimport("sys" => "path")
    @test pyis(path, sys.path)
    verpath = pyimport("sys" => ("version", "path"))
    @test verpath isa Tuple{Py, Py}
    @test pyis(verpath[1], ver)
    @test pyis(verpath[2], path)
end

@testset "consts" begin
    @test pybuiltins.None isa Py
    @test pystr(String, pybuiltins.None) == "None"
end

@testset "str" begin
    @test pyisinstance(pystr("foo"), pybuiltins.str)
    @test pyeq(Bool, pystr(pystr("foo")), pystr("foo"))
    @test pyeq(Bool, pystr(SubString("foobarbaz", 4:6)), pystr("bar"))
    @test pyeq(Bool, pystr('x'), pystr("x"))
    @test pystr(String, pybuiltins.None) === "None"
    @test pystr(String, pyint(123)) === "123"
    @test pystr(String, pystr("foo")) === "foo"
end

@testset "bytes" begin
    @test pyisinstance(pybytes(UInt8[1,2,3]), pybuiltins.bytes)
    @test pyeq(Bool, pybytes(pylist([1,2,3])), pybytes(UInt8[1,2,3]))
    @test pyeq(Bool, pybytes(b"foo"), pystr("foo").encode("ascii"))
    @test pyeq(Bool, pybytes(codeunits(SubString("foobarbaz", 4:6))), pystr("bar").encode("ascii"))
    @test pybytes(Vector, pylist([1,2,3])) == UInt8[1,2,3]
    @test pybytes(Vector{UInt8}, pylist([1,2,3])) == UInt8[1,2,3]
    @test pybytes(Base.CodeUnits, pystr("foo").encode("ascii")) == b"foo"
    @test pybytes(Base.CodeUnits{UInt8,String}, pystr("bar").encode("ascii")) == b"bar"
end

@testset "tuple" begin
    z = pytuple()
    @test pyisinstance(z, pybuiltins.tuple)
    @test pylen(z) == 0
    x = pytuple((1,2,3))
    @test pyisinstance(x, pybuiltins.tuple)
    @test pylen(x) == 3
    @test pyeq(Bool, pygetitem(x, 0), 1)
    @test pyeq(Bool, pygetitem(x, 1), 2)
    @test pyeq(Bool, pygetitem(x, 2), 3)
    @test pyeq(Bool, pytuple([1,2,3]), x)
    @test pyeq(Bool, pytuple(i+1 for i in 0:10 if i<3), x)
    @test pyeq(Bool, pytuple(pytuple((1,2,3))), x)
    @test pyeq(Bool, pytuple(pylist([1,2,3])), x)
end

@testset "list" begin
    z = pylist()
    @test pyisinstance(z, pybuiltins.list)
    @test pylen(z) == 0
    x = pylist((1,2,3))
    @test pyisinstance(x, pybuiltins.list)
    @test pylen(x) == 3
    @test pyeq(Bool, pygetitem(x, 0), 1)
    @test pyeq(Bool, pygetitem(x, 1), 2)
    @test pyeq(Bool, pygetitem(x, 2), 3)
    @test pyeq(Bool, pylist([1,2,3]), x)
    @test pyeq(Bool, pylist(i+1 for i in 0:10 if i<3), x)
    @test pyeq(Bool, pylist(pylist((1,2,3))), x)
    @test pyeq(Bool, pylist(pytuple([1,2,3])), x)
    @test pyeq(Bool, pycollist([1,2,3]), pylist([1,2,3]))
    @test pyeq(Bool, pycollist([1 2; 3 4]), pylist((pylist([1,3]), pylist([2,4]))))
    @test pyeq(Bool, pyrowlist([1,2,3]), pylist([1,2,3]))
    @test pyeq(Bool, pyrowlist([1 2; 3 4]), pylist((pylist([1,2]), pylist([3,4]))))
end

@testset "dict" begin
    z = pydict()
    @test pyisinstance(z, pybuiltins.dict)
    @test pylen(z) == 0
    x = pydict(foo=1, bar=2)
    @test pyisinstance(x, pybuiltins.dict)
    @test pylen(x) == 2
    @test pyeq(Bool, pygetitem(x, "foo"), 1)
    @test pyeq(Bool, pygetitem(x, "bar"), 2)
    @test pyeq(Bool, pydict(["foo"=>1, "bar"=>2]), x)
    @test pyeq(Bool, pydict([("foo"=>1), ("bar"=>2)]), x)
    @test pyeq(Bool, pydict(Dict("foo"=>1, "bar"=>2)), x)
    @test pyeq(Bool, pydict((foo=1, bar=2)), x)
    @test pyeq(Bool, pydict(x), x)
end

@testset "bool" begin
    @test pyis(pybool(), pybuiltins.False)
    @test pyis(pybool(false), pybuiltins.False)
    @test pyis(pybool(true), pybuiltins.True)
    @test pyis(pybool(0.0), pybuiltins.False)
    @test pyis(pybool(-1.2), pybuiltins.True)
    @test pyis(pybool(pybuiltins.None), pybuiltins.False)
    @test pyis(pybool(pylist()), pybuiltins.False)
    @test pyis(pybool(pylist([1,2,3])), pybuiltins.True)
end

@testset "int" begin
    @test pyisinstance(pyint(), pybuiltins.int)
    @test pystr(String, pyint()) == "0"
    x = 123
    y = pyint(x)
    @test pyisinstance(y, pybuiltins.int)
    @test pystr(String, y) == string(x)
    x = BigInt(123) << 200
    y = pyint(x)
    @test pyisinstance(y, pybuiltins.int)
    @test pystr(String, y) == string(x)
    x = UInt(123)
    y = pyint(x)
    @test pyisinstance(y, pybuiltins.int)
    @test pystr(String, y) == string(x)
    x = UInt128(123) << 100
    y = pyint(x)
    @test pyisinstance(y, pybuiltins.int)
    @test pystr(String, y) == string(x)
    @test pyeq(Bool, pyint(pyint(123)), pyint(123))
    @test pyeq(Bool, pyint(pyfloat(12.3)), pyint(12))
end

@testset "float" begin
    y = pyfloat()
    @test pyisinstance(y, pybuiltins.float)
    @test pyeq(Bool, y, pyint(0))
    x = 123
    y = pyfloat(x)
    @test pyisinstance(y, pybuiltins.float)
    @test pyeq(Bool, y, pyint(x))
    x = 0.25
    y = pyfloat(x)
    @test pyisinstance(y, pybuiltins.float)
    @test pyeq(Bool, y, pytruediv(1, 4))
    x = 1//4
    y = pyfloat(x)
    @test pyisinstance(y, pybuiltins.float)
    @test pyeq(Bool, y, pyfloat(float(x)))
    @test pyeq(Bool, pyfloat(pyfloat(12.3)), pyfloat(12.3))
    @test pyeq(Bool, pyfloat(pyint(123)), pyfloat(123))
end

@testset "complex" begin
    y = pycomplex()
    @test pyisinstance(y, pybuiltins.complex)
    @test pyeq(Bool, y, pyint(0))
    x = 12.3
    y = pycomplex(x)
    @test pyisinstance(y, pybuiltins.complex)
    @test pyeq(Bool, y, pyfloat(x))
    xr, xi = 12, 34
    y = pycomplex(xr, xi)
    @test pyisinstance(y, pybuiltins.complex)
    @test pyeq(Bool, y.real, pyfloat(xr))
    @test pyeq(Bool, y.imag, pyfloat(xi))
    x = Complex(12, 34)
    y = pycomplex(x)
    @test pyisinstance(y, pybuiltins.complex)
    @test pyeq(Bool, y.real, pyfloat(real(x)))
    @test pyeq(Bool, y.imag, pyfloat(imag(x)))
    @test pyeq(Bool, pycomplex(y), y)
    @test pyeq(Bool, pycomplex(pyint(12), pyint(34)), y)
end

@testset "set" begin
    y = pyset()
    yf = pyfrozenset()
    @test pyisinstance(y, pybuiltins.set)
    @test pylen(y) == 0
    @test pyisinstance(yf, pybuiltins.frozenset)
    @test pylen(yf) == 0
    @test pyeq(Bool, y, yf)
    x = [1,2,3,2,1]
    y = pyset(x)
    yf = pyfrozenset(x)
    @test pyisinstance(y, pybuiltins.set)
    @test pylen(y) == 3
    @test pycontains(y, 1)
    @test pycontains(y, 2)
    @test pycontains(y, 3)
    @test pyeq(Bool, pyset(y), y)
    @test pyisinstance(yf, pybuiltins.frozenset)
    @test pylen(yf) == 3
    @test pycontains(yf, 1)
    @test pycontains(yf, 2)
    @test pycontains(yf, 3)
    @test pyeq(Bool, pyfrozenset(y), y)
    @test pyeq(Bool, y, yf)
end

@testset "slice" begin
    x = pyslice(12)
    @test pyisinstance(x, pybuiltins.slice)
    @test pyeq(Bool, x.start, pybuiltins.None)
    @test pyeq(Bool, x.stop, 12)
    @test pyeq(Bool, x.step, pybuiltins.None)
    x = pyslice(12, 34)
    @test pyisinstance(x, pybuiltins.slice)
    @test pyeq(Bool, x.start, 12)
    @test pyeq(Bool, x.stop, 34)
    @test pyeq(Bool, x.step, pybuiltins.None)
    x = pyslice(12, 34, 56)
    @test pyisinstance(x, pybuiltins.slice)
    @test pyeq(Bool, x.start, 12)
    @test pyeq(Bool, x.stop, 34)
    @test pyeq(Bool, x.step, 56)
end

@testset "range" begin
    x = pyrange(123)
    @test pyisinstance(x, pybuiltins.range)
    @test pyeq(Bool, x.start, 0)
    @test pyeq(Bool, x.stop, 123)
    @test pyeq(Bool, x.step, 1)
    x = pyrange(12, 123)
    @test pyisinstance(x, pybuiltins.range)
    @test pyeq(Bool, x.start, 12)
    @test pyeq(Bool, x.stop, 123)
    @test pyeq(Bool, x.step, 1)
    x = pyrange(12, 123, 3)
    @test pyisinstance(x, pybuiltins.range)
    @test pyeq(Bool, x.start, 12)
    @test pyeq(Bool, x.stop, 123)
    @test pyeq(Bool, x.step, 3)
end

@testset "none" begin
end

@testset "type" begin
    x = pytype(pyint())
    @test pyisinstance(x, pybuiltins.type)
    @test pyis(x, pybuiltins.int)
    x = pytype(pybuiltins.type)
    @test pyisinstance(x, pybuiltins.type)
    @test pyis(x, pybuiltins.type)
    x = pytype("Foo", (), (foo=1, bar=2))
    @test pyisinstance(x, pybuiltins.type)
    @test pyeq(Bool, x.__name__, "Foo")
    @test pyeq(Bool, x.foo, 1)
    @test pyeq(Bool, x.bar, 2)
    x = pyclass("Bar", (), foo=3, bar=4)
    @test pyisinstance(x, pybuiltins.type)
    @test pyeq(Bool, x.__name__, "Bar")
    @test pyeq(Bool, x.foo, 3)
    @test pyeq(Bool, x.bar, 4)
end

@testset "fraction" begin
    # TODO
end

@testset "method" begin
    # TODO
end

@testset "datetime" begin
    # TODO
end

@testset "code" begin
    # TODO
end


text:
board = [[" ", " ", " "],
         [" ", " ", " "],
         [" ", " ", " "]]

player1 = "X"
player2 = "O"

player1_win = -1
tie = 0
player2_win = 1

function print_board(board)
    println("    1   2   3")
    println(" 1  ", join(board[1], " | "))
    println("   ———┼———┼———")
    println(" 2  ", join(board[2], " | "))
    println("   ———┼———┼———")
    println(" 3  ", join(board[3], " | "))
end

function make_move!(board, player, (x, y))
    board[x][y] = player
end

function check_available(board)
    available_cells = []
    for i in 1:3, j in 1:3
        if board[i][j] == " "
            push!(available_cells, (i, j))
        end
    end
    return available_cells
end

function ai_move!(board)
    pos = rand(check_available(board))
    make_move!(board, player2, pos)
end

print_board(board)

while true
    print_board(board)
    print("Enter: ")
    txt = split(readline())
    make_move!(board, player1, (parse(Int, txt[1]), parse(Int, txt[2])))
    ai_move!(board)
    print_board(board)
end

max_stars_repo_name: soumitradev/youtube-tutorial

max_issues_repo_name: soumitradev/youtube-tutorial

max_forks_repo_name: soumitradev/youtube-tutorial


text:
vis = Visualizer()
open(vis)

# ## model
include("trajopt_model_v2.jl")

# ## horizon
h = 0.05
T = 135
Tm = 15 # mid point for a swing / stance change

# ## centroidal_quadruped
s = get_simulation("centroidal_quadruped", "flat_3D_lc", "flat")
model = s.model
env = s.env

nx = 2 * model.nq
nc = 4 #model.nc
nu = model.nu + nc + 4 * nc + nc + 4 * nc + 1
nθ = 53

# ## model
d1 = DTO.Dynamics((y, x, u, w) -> centroidal_quadruped_dyn1(
	model, env, [h], y, x, u, w), nx + nθ + nx, nx, nu)
dt = DTO.Dynamics((y, x, u, w) -> centroidal_quadruped_dynt(
	model, env, [h], y, x, u, w), nx + nθ + nx, nx + nθ + nx, nu)

dyn = [d1, [dt for t = 2:T-1]...]

# ## initial conditions
body_height = 0.25
foot_x = 0.17
foot_y = 0.17
foot_height = 0.08

q1 = zeros(model.nq)

# body position
q1[1:3] = [0.0; 0.0; body_height]
q1[4:6] = [0.0; 0.0; 0.0]

# foot1
q1[7:9] = [foot_x; foot_y; 0]

# foot2
q1[10:12] = [foot_x; -foot_y; 0]

# foot3
q1[13:15] = [-foot_x; foot_y; 0]

# foot4
q1[16:18] = [-foot_x; -foot_y; 0]

# Terminal Q
qT = copy(q1)

visualize!(vis, model, [qM2], Δt=h);

h_step = 0.20
qM0 = deepcopy(q1)
qM0[0 .+ (1:3)] += [-0.02, 0.0, 0.0] # body

qM1 = deepcopy(q1)
qM1[6 .+ (1:3)] += [+0.06, 0.0, 2h_step] # front left
qM1[0 .+ (1:3)] += [-0.04, 0.0, 0.0] # body

qM2 = deepcopy(q1)
qM2[6 .+ (1:3)] += [+0.16, 0.0, h_step] # front left
qM2[0 .+ (1:3)] += [-0.04, 0.0, 0.0] # body

qM3 = deepcopy(q1)
qM3[6 .+ (1:3)] += [+0.16, 0.0, h_step] # front left
qM3[0 .+ (1:3)] += [+0.00, 0.0, 0.0] # body

qM4 = deepcopy(q1)
qM4[6 .+ (1:3)] += [0.16, 0.0, h_step] # front left
qM4[9 .+ (1:3)] += [0.08, 0.0, 2h_step] # front left
qM4[0 .+ (1:3)] += [0.00, 0.0, 0.0] # body

qM5 = deepcopy(q1)
qM5[6 .+ (1:3)] += [0.16, 0.0, h_step] # front left
qM5[9 .+ (1:3)] += [0.16, 0.0, h_step] # front left
qM5[0 .+ (1:3)] += [0.00, 0.0, 0.0] # body

qT = deepcopy(q1)
qT[6 .+ (1:3)] += [0.16, 0.0, h_step] # front left
qT[9 .+ (1:3)] += [0.16, 0.0, h_step] # front left
qT[0 .+ (1:3)] += [0.08, 0.0, 0.0] # body

set_robot!(vis, model, q1)
set_robot!(vis, model, qM0)
set_robot!(vis, model, qM1)
set_robot!(vis, model, qM2)
set_robot!(vis, model, qM3)
set_robot!(vis, model, qM4)
set_robot!(vis, model, qM5)
set_robot!(vis, model, qT)

q_ref = [q1,
	linear_interpolation(q1, qM0, Tm)...,
	linear_interpolation(qM0, qM1, Tm)...,
	linear_interpolation(qM1, qM2, Tm)...,
	linear_interpolation(qM2, qM2, Tm)...,
	linear_interpolation(qM2, qM3, Tm)...,
	linear_interpolation(qM3, qM4, Tm)...,
	linear_interpolation(qM4, qM5, Tm)...,
	linear_interpolation(qM5, qT, 2Tm)...]

x_ref = [[q_ref[t]; q_ref[t+1]] for t = 1:T]
x1 = x_ref[1]
xM = x_ref[Tm]
xT = x_ref[T]

visualize!(vis, model, q_ref, Δt=h);

plot(hcat(q_ref...)', labels="")

# ## objective
obj = DTO.Cost{Float64}[]
for t = 1:T
    if t == T
        function objT(x, u, w)
            J = 0.0
            v = (x[model.nq .+ (1:model.nq)] - x[1:model.nq]) ./ h
            J += 0.5 * 1.0e-2 * dot(v, v)
            J += 100 * transpose(x[1:nx] - x_ref[t]) * Diagonal(1000.0 * ones(nx)) * (x[1:nx] - x_ref[t])
            return J / T
        end
        push!(obj, DTO.Cost(objT, nx + nθ + nx, 0))
    elseif t == 1
        function obj1(x, u, w)
            J = 0.0
            v = (x[model.nq .+ (1:model.nq)] - x[1:model.nq]) ./ h
            J += 0.5 * 1.0e-2 * dot(v, v)
            J += 100 * transpose(x[1:nx] - x_ref[t]) * Diagonal(1000.0 * ones(nx)) * (x[1:nx] - x_ref[t])
            J += 0.5 * transpose(u[1:model.nu]) * Diagonal(1.0e-3 * ones(model.nu)) * u[1:model.nu]
            J += 1000.0 * u[end] # slack
            return J / T
        end
        push!(obj, DTO.Cost(obj1, nx, nu))
    else
        function objt(x, u, w)
            J = 0.0
            v = (x[model.nq .+ (1:model.nq)] - x[1:model.nq]) ./ h
            J += 0.5 * 1.0e-2 * dot(v, v)
            u_previous = x[nx .+ (1:53)]
            u_control = u
            w = (u_control - u_previous) ./ h
            J += 0.5 * 1.0e-3 * dot(w, w)
            J += 100 * transpose(x[1:nx] - x_ref[t]) * Diagonal(1000.0 * ones(nx)) * (x[1:nx] - x_ref[t])
            J += 0.5 * transpose(u[1:model.nu]) * Diagonal(1.0e-3 * ones(model.nu)) * u[1:model.nu]
            J += 1000.0 * u[end] # slack
            return J / T
        end
        push!(obj, DTO.Cost(objt, nx + nθ + nx, nu))
    end
end

# ## constraints
# initial condition
xl1 = [q1; q1]
xu1 = [q1; q1]
xlt = [-Inf * ones(nx); -Inf * ones(nθ); -Inf * ones(nx)]
xut = [Inf * ones(nx); Inf * ones(nθ); Inf * ones(nx)]

# final condition
xlT = [-Inf * ones(nq); qT; -Inf * ones(nθ); -Inf * ones(nx)]
xuT = [Inf * ones(nq); qT; Inf * ones(nθ); Inf * ones(nx)]

ul = [-Inf * ones(model.nu); zeros(nu - model.nu)]
uu = [Inf * ones(model.nu); Inf * ones(nu - model.nu)]

# bnd1 = DTO.Bound(nx, nu, state_lower=xl1, state_upper=xu1, action_lower=ul, action_upper=uu)
# bndt = DTO.Bound(nx + nθ + nx, nu, state_lower=xlt, state_upper=xut, action_lower=ul, action_upper=uu)
# bndT = DTO.Bound(nx + nθ + nx, 0, state_lower=xlT, state_upper=xuT)
# bnds = [bnd1, [bndt for t = 2:T-1]..., bndT];

bnds = DTO.Bound{Float64}[]
push!(bnds, DTO.Bound(nx, nu, state_lower=xl1, state_upper=xu1, action_lower=ul, action_upper=uu))
for t = 2:T-1
    push!(bnds, DTO.Bound(nx + nθ + nx, nu,
        state_lower=[-Inf * ones(nq); -Inf * ones(nq); -Inf * ones(nθ); -Inf * ones(nx)],
        state_upper=[Inf * ones(nq); Inf * ones(nq); Inf * ones(nθ); Inf * ones(nx)],
        action_lower=ul, action_upper=uu))
end
push!(bnds, DTO.Bound(nx + nθ + nx, 0, state_lower=xlT, state_upper=xuT))

cons = DTO.Constraint{Float64}[]
for t = 1:T
    if t == 1
        function constraints_1(x, u, w)
            [
            # equality (16)
            contact_constraints_equality(model, env, h, x, u, w);
            # inequality (28)
            contact_constraints_inequality_1(model, env, h, x, u, w);

            # body/feet constraints
            x[3] - x_ref[t][3]; # body height
            # x[model.nq + 3] - x_ref[t][model.nq + 3]; # body height
            # x[9:3:18] - x_ref[t][9:3:18];
			x[12 .+ (1:3)] - q1[12 .+ (1:3)] # back_feet
			x[15 .+ (1:3)] - q1[15 .+ (1:3)] # back_feet
            ]
        end
        push!(cons, DTO.Constraint(constraints_1, nx, nu, idx_ineq=collect(16 .+ (1:28))))
    elseif t == T
        function constraints_T(x, u, w)
            [
            # inequality (8)
            contact_constraints_inequality_T(model, env, h, x, u, w);

            # body/feet constraints
            x[3] - x_ref[t][3]; # body height
            # x[model.nq + 3] - x_ref[t][model.nq + 3]; # body height
            # x[9:3:18] - x_ref[t][9:3:18];
			x[12 .+ (1:3)] - q1[12 .+ (1:3)] # back_feet
			x[15 .+ (1:3)] - q1[15 .+ (1:3)] # back_feet
            ]
        end
        push!(cons, DTO.Constraint(constraints_T, nx + nθ + nx, nu, idx_ineq=collect(0 .+ (1:8))))
    else
        function constraints_t(x, u, w)
            [
            # equality (16)
            contact_constraints_equality(model, env, h, x, u, w);
            # inequality (32)
            contact_constraints_inequality_t(model, env, h, x, u, w);

            # body/feet constraints
            x[3] - x_ref[t][3]; # body height
            # x[model.nq + 3] - x_ref[t][model.nq + 3]; # body height
            # x[9:3:18] - x_ref[t][9:3:18];
			x[12 .+ (1:3)] - q1[12 .+ (1:3)] # back_feet
			x[15 .+ (1:3)] - q1[15 .+ (1:3)] # back_feet
            ]
        end
        push!(cons, DTO.Constraint(constraints_t, nx + nθ + nx, nu, idx_ineq=collect(16 .+ (1:32))) )
    end
end

# ## problem
tolerance = 1.0e-3
p = DTO.solver(dyn, obj, cons, bnds,
    options=DTO.Options(
        max_iter=500,
        max_cpu_time=30000.0,
        tol=tolerance,
        constr_viol_tol=tolerance,
        ))

# ## initialize
x_interpolation = [x_ref[1], [[x_ref[t]; zeros(nθ); zeros(nx)] for t = 2:T]...]
u_guess = [1.0e-4 * rand(nu) for t = 1:T-1] # may need to run more than once to get good trajectory
DTO.initialize_states!(p, x_interpolation)
DTO.initialize_controls!(p, u_guess)

# ## solve
@time DTO.solve!(p)

# ## solution
x_sol, u_sol = DTO.get_trajectory(p)
@show x_sol[1]
@show x_sol[T]
maximum([u[end] for u in u_sol[1:end-1]])

# ## visualize
visualize!(vis, model, [x_sol[1][1:nq], [x[nq .+ (1:nq)] for x in x_sol]...], Δt=h);

q_opt = [x_sol[1][1:model.nq], [x[model.nq .+ (1:model.nq)] for x in x_sol]...]
v_opt = [(x[model.nq .+ (1:model.nq)] - x[0 .+ (1:model.nq)]) ./ h for x in x_sol]
u_opt = [u[1:model.nu] for u in u_sol]
λ_opt = [u[model.nu .+ (1:4)] for u in u_sol]
b_opt = [u[model.nu + 4 .+ (1:16)] for u in u_sol]

q_opt = [x_sol[1][1:model.nq], [x[model.nq .+ (1:model.nq)] for x in x_sol]...]
v_opt = [(x[model.nq .+ (1:model.nq)] - x[0 .+ (1:model.nq)]) ./ h for x in x_sol]
u_opt = [u[1:model.nu] for u in u_sol]
γ_opt = [u[model.nu .+ (1:4)] for u in u_sol]
b_opt = [u[model.nu + 4 .+ (1:16)] for u in u_sol]
ψ_opt = [u[model.nu + 4 + 16 .+ (1:4)] for u in u_sol]
η_opt = [u[model.nu + 4 + 16 + 4 .+ (1:16)] for u in u_sol]

qm_ = copy(q_opt)
um_ = copy(u_opt)
γm_ = copy(γ_opt)
bm_ = copy(b_opt)
ψm_ = copy(ψ_opt)
ηm_ = copy(η_opt)
μm = model.μ_world

hm = h
timesteps = range(0.0, stop=(h * (length(qm_) - 2)), length=(length(qm_) - 2))
plot(hcat(qm_...)', labels="")

plot(timesteps, hcat(qm_[2:end-1]...)', labels="")
plot(timesteps, hcat(um_...)', labels="")
plot(timesteps, hcat(γm_...)', labels="")
plot(timesteps, hcat(bm_...)', labels="")
plot(timesteps, hcat(ψm_...)', labels="")
plot(timesteps, hcat(ηm_...)', labels="")

################################################################################
# Trajectory Extension
################################################################################
qm = [qm_; [qm_[end] for i=1:50]]
um = [um_; [um_[end] for i=1:50]]
γm = [γm_; [γm_[end] for i=1:50]]
bm = [bm_; [bm_[end] for i=1:50]]
ψm = [ψm_; [ψm_[end] for i=1:50]]
ηm = [ηm_; [ηm_[end] for i=1:50]]

using JLD2
@save joinpath(@__DIR__, "step_over_box_v3.jld2") qm um γm bm ψm ηm μm hm
@load joinpath(@__DIR__, "step_over_box_v3.jld2") qm um γm bm ψm ηm μm hm


text:
using LinearAlgebra
using QuantumGate

function single_qubit_gate(cite::Int, N::Int, Operator)
    A = σ⁰

    if cite == 1
        A = Operator
    end

    for i in 2:N
        if cite == i
            A = A ⊗ Operator
        else
            A = A ⊗ σ⁰
        end
    end

    return A
end

function Xgate(cite::Int, N::Int)
    single_qubit_gate(cite, N, σˣ)
end

function Ygate(cite::Int, N::Int)
    single_qubit_gate(cite, N, σʸ)
end

function Zgate(cite::Int, N::Int)
    single_qubit_gate(cite, N, σᶻ)
end

function Hadamardgate(cite::Int, N::Int)
    single_qubit_gate(cite, N, Had)
end

function Sgate(cite::Int, N::Int)
    single_qubit_gate(cite, N, Smat)
end

function Tgate(cite::Int, N::Int)
    single_qubit_gate(cite, N, Tmat)
end

function Rˣgate(cite::Int, N::Int, θ::Float64)
    single_qubit_gate(cite, N, Rˣmat(θ))
end

function Rʸgate(cite::Int, N::Int, θ::Float64)
    single_qubit_gate(cite, N, Rʸmat(θ))
end

function Rᶻgate(cite::Int, N::Int, θ::Float64)
    single_qubit_gate(cite, N, Rᶻmat(θ))
end

function controlled_Ugate(control::Int, target::Int, N::Int, U)
    A = single_qubit_gate(control, N, P⁰)

    if control == 1
        B = P¹
    elseif target == 1
        B = U
    else
        B = σ⁰
    end

    for i in 2:N
        if control == i
            B = B ⊗ P¹
        elseif target == i
            B = B ⊗ U
        else
            B = B ⊗ σ⁰
        end
    end
    return A + B
end

function CNOTgate(control::Int, target::Int, N::Int)
    controlled_Ugate(control, target, N, σˣ)
end

function SWAPgate(i, j, N)
    CNOTgate(i, j, N) * CNOTgate(j, i, N) * CNOTgate(i, j, N)
end


text:
#
# This file is a part of MolecularGraph.jl
# Licensed under the MIT License http://opensource.org/licenses/MIT
#

export
    LineGraph, linegraph

struct LineGraphNode <: AbstractNode
    n1::Int
    n2::Int
end

struct LineGraphEdge <: UndirectedEdge
    node::Int
end

struct LineGraph <: OrderedGraph
    neighbormap::Vector{Dict{Int,Int}}
    edges::Vector{Tuple{Int,Int}}
    nodeattrs::Vector{LineGraphNode}
    edgeattrs::Vector{LineGraphEdge}
    cache::Dict{Symbol,Any}
end

"""
    linegraph(G::AbstractGraph) -> LineGraph

Generate line graph.
"""
function linegraph(G::AbstractGraph)
    L = LineGraph([], [], [], [], Dict())
    nmap = Dict{Int,Int}()
    for (i, e) in enumerate(edgeset(G))
        (u, v) = getedge(G, e)
        push!(L.nodeattrs, LineGraphNode(u, v))
        push!(L.neighbormap, Dict())
        nmap[e] = i
    end
    ecnt = 1
    for n in nodeset(G)
        degree(G, n) < 2 && continue
        for (e1, e2) in combinations(incidences(G, n))
            ne1 = nmap[e1]
            ne2 = nmap[e2]
            push!(L.edges, (ne1, ne2))
            push!(L.edgeattrs, LineGraphEdge(n))
            L.neighbormap[ne1][ecnt] = ne2
            L.neighbormap[ne2][ecnt] = ne1
            ecnt += 1
        end
    end
    return L
end

function linegraph(G::OrderedGraph)
    L = LineGraph([], [], [], [], Dict())
    for (i, edge) in enumerate(edgesiter(G))
        push!(L.nodeattrs, LineGraphNode(edge...))
        push!(L.neighbormap, Dict())
    end
    ecnt = 1
    for i in 1:nodecount(G)
        degree(G, i) < 2 && continue
        for (e1, e2) in combinations(incidences(G, i))
            push!(L.edges, (e1, e2))
            push!(L.edgeattrs, LineGraphEdge(i))
            L.neighbormap[e1][ecnt] = e2
            L.neighbormap[e2][ecnt] = e1
            ecnt += 1
        end
    end
    return L
end


text:
#####
##### `frule`/`rrule`
#####

"""
    frule((Δf, Δx...), f, x...)

Expressing the output of `f(x...)` as `Ω`, return the tuple:

    (Ω, ΔΩ)

The second return value is the differential w.r.t. the output.

If no method matching `frule((Δf, Δx...), f, x...)` has been defined, then return `nothing`.

Examples:

unary input, unary output scalar function:

```jldoctest frule
julia> dself = NO_FIELDS;

julia> x = rand()
0.8236475079774124

julia> sinx, Δsinx = frule((dself, 1), sin, x)
(0.7336293678134624, 0.6795498147167869)

julia> sinx == sin(x)
true

julia> Δsinx == cos(x)
true
```

unary input, binary output scalar function:

```jldoctest frule
julia> sincosx, Δsincosx = frule((dself, 1), sincos, x);

julia> sincosx == sincos(x)
true

julia> Δsincosx == (cos(x), -sin(x))
true
```

See also: [`rrule`](@ref), [`@scalar_rule`](@ref)
"""
frule(::Any, ::Vararg{Any}; kwargs...) = nothing

"""
    rrule(f, x...)

Expressing `x` as the tuple `(x₁, x₂, ...)` and the output tuple of `f(x...)`
as `Ω`, return the tuple:

    (Ω, (Ω̄₁, Ω̄₂, ...) -> (s̄elf, x̄₁, x̄₂, ...))

Where the second return value is the the propagation rule or pullback.
It takes in differentials corresponding to the outputs (`x̄₁, x̄₂, ...`),
and `s̄elf`, the internal values of the function itself (for closures)

If no method matching `rrule(f, xs...)` has been defined, then return `nothing`.

Examples:

unary input, unary output scalar function:

```jldoctest
julia> x = rand();

julia> sinx, sin_pullback = rrule(sin, x);

julia> sinx == sin(x)
true

julia> sin_pullback(1) == (NO_FIELDS, cos(x))
true
```

binary input, unary output scalar function:

```jldoctest
julia> x, y = rand(2);

julia> hypotxy, hypot_pullback = rrule(hypot, x, y);

julia> hypotxy == hypot(x, y)
true

julia> hypot_pullback(1) == (NO_FIELDS, (x / hypot(x, y)), (y / hypot(x, y)))
true
```

See also: [`frule`](@ref), [`@scalar_rule`](@ref)
"""
rrule(::Any, ::Vararg{Any}; kwargs...) = nothing


text:
"""
Derivatives-related functions
"""

## Derivative types
abstract type AbstractDiffMethod end

struct ForwardAD <: AbstractDiffMethod end
struct ReverseAD <: AbstractDiffMethod end
struct RevZyg <: AbstractDiffMethod end  # only used for gradients (not jacobians)
struct ForwardFD <: AbstractDiffMethod end
struct CentralFD <: AbstractDiffMethod end
struct ComplexStep <: AbstractDiffMethod end
struct UserDeriv <: AbstractDiffMethod end   # user-specified derivatives

FD = Union{ForwardFD, CentralFD, ComplexStep}

"""
convert to type used in FiniteDiff package
"""
function finitediff_type(dtype)
    if isa(dtype, ForwardFD)
        fdtype = Val{:forward}
    elseif isa(dtype, CentralFD)
        fdtype = Val{:central}
    elseif isa(dtype, ComplexStep)
        fdtype = Val{:complex}
    end
    return fdtype
end

# internally used cache for gradients and jacobians (separately)
struct GradOrJacCache{T1,T2,T3,T4}
    f!::T1
    work::T2
    cache::T3
    dtype::T4
end

function gradientcache(dtype::FD, func!, nx, ng)

    df = zeros(nx)
    x = zeros(nx)
    fdtype = finitediff_type(dtype)
    cache = FiniteDiff.GradientCache(df, x, fdtype)

    return GradOrJacCache(func!, nothing, cache, dtype)
end

function gradient!(df, x, cache::GradOrJacCache{T1,T2,T3,T4}
    where {T1,T2,T3,T4<:FD})

    FiniteDiff.finite_difference_gradient!(df, cache.f!, x, cache.cache)

    return nothing
end

max_stars_repo_name: Yuricst/joptimise

max_issues_repo_name: Yuricst/joptimise

max_forks_repo_name: Yuricst/joptimise


text:
###########################################################################
#                          PowerModels Extensions                         #
###########################################################################

# get rid of angle difference bounds, which will be included in other parts of code
# the valid inequality cuts for voltage products are still kept
function PM.constraint_voltage_angle_difference(pm::PM.AbstractWModels, n::Int, f_idx, angmin, angmax)
    i, f_bus, t_bus = f_idx

    w_fr = PM.var(pm, n, :w, f_bus)
    w_to = PM.var(pm, n, :w, t_bus)
    wr   = PM.var(pm, n, :wr, (f_bus, t_bus))
    wi   = PM.var(pm, n, :wi, (f_bus, t_bus))

    PM.cut_complex_product_and_angle_difference(pm.model, w_fr, w_to, wr, wi, angmin, angmax)
end

# modified from original objective_min_fuel_and_flow_cost to allow for
# the case with no generators
function objective_min_fuel_and_flow_cost_mod(pm::PM.AbstractPowerModel; kwargs...)
    model = PM.check_cost_models(pm)
    if model == 1
        return PM.objective_min_fuel_and_flow_cost_pwl(pm; kwargs...)
    elseif model == 2
        return PM.objective_min_fuel_and_flow_cost_polynomial(pm; kwargs...)
    elseif model === nothing
        return JuMP.@objective(pm.model, Min, 0)
    else
        Memento.error(_LOGGER, "Only cost models of types 1 and 2 are supported at this time, given cost model type of $(model)")
    end
end

# This is essentially copied from original PM.build_opf
function PM.build_opf(pm::NodeWRMPowerModel)
    PM.variable_bus_voltage(pm)
    PM.variable_gen_power(pm)
    PM.variable_branch_power(pm)
    PM.variable_dcline_power(pm)

    objective_min_fuel_and_flow_cost_mod(pm)

    PM.constraint_model_voltage(pm)

    :cut_bus in keys(PM.ref(pm)) ? cut_bus = PM.ids(pm, :cut_bus) : cut_bus = [] # This is needed in order to be compatible with network decomposition later
    for i in setdiff(PM.ids(pm, :bus), cut_bus)
        PM.constraint_power_balance(pm, i)
    end

    for i in PM.ids(pm, :branch)
        PM.constraint_ohms_yt_from(pm, i)
        PM.constraint_ohms_yt_to(pm, i)
        PM.constraint_voltage_angle_difference(pm, i)
        PM.constraint_thermal_limit_from(pm, i)
        PM.constraint_thermal_limit_to(pm, i)
    end

    for i in PM.ids(pm, :dcline)
        PM.constraint_dcline_power_losses(pm, i)
    end
end

max_stars_repo_name: kibaekkim/ComplexBC

max_issues_repo_name: kibaekkim/ComplexBC

max_forks_repo_name: kibaekkim/ComplexBC


text:
# Kaggle tutorial:
# https://www.kaggle.com/c/street-view-getting-started-with-julia/overview/julia-tutorial

# Load packages

using Images
using DataFrames
using CSV
using DecisionTree
using Statistics

;

"""
dataset_type:
    "train" or "test"
labels_info:
    IDs of each image to read
image_size:
    amount of pixels in image
dataset_dir_path:
    path to directory with "trainResized" and "testResized" directories.
"""
function read_image(dataset_type, labels_info, image_size, dataset_dir_path)
    # Initialize X matrix
    X = zeros(size(labels_info, 1), image_size)
    @show size(labels_info)
    for (idx, image_id) in enumerate(labels_info[:ID])
        # Load image
        image_full_path = "$(dataset_dir_path)/$(dataset_type)Resized/$(image_id).Bmp"
        image = load(image_full_path)
        image_greyscale = Gray.(image)
        # Transform image matrix to a vector and store it in data matrix
        X[idx, :] = reshape(image_greyscale, 1, image_size)
    end
    return X
end

;

# The images in the trainResized and testResized directories are 20x20 pixels,
# so image_size should be set to 400.

image_size = 20 * 20
dataset_dir_path = "/home/max/Documents/ai/julia/julia_kaggle/data"

# Load labels (ID and Class) for train and test

# train_set
train_labels_info = CSV.read("$(dataset_dir_path)/trainLabels.csv")
@show size(train_labels_info)
@show names(train_labels_info)

# test_set
test_labels_info = CSV.read("$(dataset_dir_path)/sampleSubmission.csv")
@show size(test_labels_info)
@show names(test_labels_info)

# Show some info about these labels
@show first(train_labels_info, 5)
@show first(test_labels_info, 5)

;

# Read training and testing sets
@time X_train = read_image("train", train_labels_info, image_size, dataset_dir_path)
@time X_test = read_image("test", test_labels_info, image_size, dataset_dir_path)

@show size(X_train)
@show size(X_test)

;

# Create Y vector from each image.
# Target variable == first character of 'Class' column

y_train = map(col_value -> col_value[1], train_labels_info[:Class])
y_train = convert(Array{Int64, 1}, y_train)

@show size(y_train)
@show typeof(y_train)
@show y_train[1:10]

;

# Create a training model

# Learn patterns in the images that identify the character in the label.
# Use Julia version of Random Forest algorithm.
# Parameters to set:
#     number of features to choose at each split
#         sqrt(number_of_features)
#     number of trees
#         larger is better, but it takes more time to train
#     ratio of sampling
#         usually 1.0

num_of_features = sqrt(size(X_train, 2))
num_of_trees = 50
ratio_of_sampling = 1.0

# build_forest(labels::Array{T,1}, features::Array{S,2}, n_subfeatures, n_trees, partial_sampling)

model = build_forest(
    y_train, X_train,  # labels, features
    num_of_features,
    num_of_trees,
    ratio_of_sampling
)

@show model

;

# Use trained model: identify the characters in the test data

test_set_predictions = apply_forest(model, X_test)

;

# Convert predictions back to characters
test_labels_info[:Class] = convert(Array{Char, 1}, test_set_predictions)

;

# Print model accuracy

n_folds = 4

# nfoldCV_forest(
#     labels::Array{T,1}, features::Array{S,2},
#     n_folds::Integer, n_subfeatures::Integer, n_trees::Integer, partial_sampling::Float64,
# )

accuracy = nfoldCV_forest(
    y_train, X_train,  # labels, features
    n_folds,
    convert(Int, num_of_features),
    num_of_trees,
    ratio_of_sampling
)

@show mean(accuracy)

;

# Save predictions

CSV.write(
    "$(dataset_dir_path)/julia_submission.csv",
    test_labels_info,
    writeheader=true
)

;

max_stars_repo_name: honchardev/Fun

max_issues_repo_name: honchardev/Fun

max_forks_repo_name: honchardev/Fun


text:
export PolyVar, @polyvar, @ncpolyvar
export polyvecvar

function polyarrayvar(::Type{PV}, prefix, indices...) where {PV}
    map(i -> PV("$(prefix)[$(join(i, ","))]"), Iterators.product(indices...))
end

function buildpolyvar(::Type{PV}, var) where {PV}
    if isa(var, Symbol)
        var, :($(esc(var)) = $PV($"$var"))
    else
        isa(var, Expr) || error("Expected $var to be a variable name")
        Base.Meta.isexpr(var, :ref) || error("Expected $var to be of the form varname[idxset]")
        (2 ≤ length(var.args)) || error("Expected $var to have at least one index set")
        varname = var.args[1]
        prefix = string(varname)
        varname, :($(esc(varname)) = polyarrayvar($PV, $prefix, $(esc.(var.args[2:end])...)))
    end
end

function buildpolyvars(::Type{PV}, args) where {PV}
    vars = Symbol[]
    exprs = []
    for arg in args
        var, expr = buildpolyvar(PV, arg)
        push!(vars, var)
        push!(exprs, expr)
    end
    vars, exprs
end

# Variable vector x returned garanteed to be sorted so that if p is built with x then vars(p) == x
macro polyvar(args...)
    vars, exprs = buildpolyvars(PolyVar{true}, args)
    :($(foldl((x,y) -> :($x; $y), exprs, init=:())); $(Expr(:tuple, esc.(vars)...)))
end

macro ncpolyvar(args...)
    vars, exprs = buildpolyvars(PolyVar{false}, args)
    :($(foldl((x,y) -> :($x; $y), exprs, init=:())); $(Expr(:tuple, esc.(vars)...)))
end

struct PolyVar{C} <: AbstractVariable
    id::Int
    name::String

    function PolyVar{C}(name::AbstractString) where {C}
        # gensym returns something like Symbol("##42")
        # we first remove "##" and then parse it into an Int
        id = parse(Int, string(gensym())[3:end])
        new(id, convert(String, name))
    end
end

Base.hash(x::PolyVar, u::UInt) = hash(x.id, u)
Base.broadcastable(x::PolyVar) = Ref(x)

MP.name(v::PolyVar) = v.name
function MP.name_base_indices(v::PolyVar)
    splits = split(v.name, r"[\[,\]]\s*", keepempty=false)
    if length(splits) == 1
        return v.name, Int[]
    else
        return splits[1], parse.(Int, splits[2:end])
    end
end

MP.monomial(v::PolyVar) = Monomial(v)
_vars(v::PolyVar) = [v]

iscomm(::Type{PolyVar{C}}) where {C} = C

function mergevars(varsvec::Vector{Vector{PV}}) where {PV<:PolyVar}
    n = length(varsvec)
    is = ones(Int, n)
    maps = zeros.(Int, length.(varsvec))
    nonempty = BitSet(findall(!isempty, varsvec))
    vars = Vector{PV}()
    while !isempty(nonempty)
        imin = 0
        for i in nonempty
            if imin == 0 || varsvec[i][is[i]] > varsvec[imin][is[imin]]
                imin = i
            end
        end
        var = varsvec[imin][is[imin]]
        push!(vars, var)
        for i in nonempty
            if var == varsvec[i][is[i]]
                maps[i][is[i]] = length(vars)
                if is[i] == length(varsvec[i])
                    pop!(nonempty, i)
                else
                    is[i] += 1
                end
            end
        end
    end
    vars, maps
end


text:
using OrdinaryDiffEq, Test

@testset "u_modifed Tests" begin
  T = 1000.0; Ttr = 0.0;
  d0=1e-9; threshold=10^4*d0;
  dt = 0.1; diff_eq_kwargs = Dict()

  @inbounds function eom_lorenz!(du, u, p, t)
    σ = 10.0; ρ = 28.0; β = 8/3
    du[1] = σ*(u[2]-u[1])
    du[2] = u[1]*(ρ-u[3]) - u[2]
    du[3] = u[1]*u[2] - β*u[3]
  end

  prob = OrdinaryDiffEq.ODEProblem(eom_lorenz!, [0., 10.0, 0], (zero(T), T))
  integ1 = init(prob, Tsit5(); diff_eq_kwargs...)
  @. prob.u0 = prob.u0 + d0/sqrt(3)
  integ2 = init(prob, Tsit5(); diff_eq_kwargs...)
  integ1.opts.advance_to_tstop=true
  integ2.opts.advance_to_tstop=true

  dist = d0
  λ = zero(eltype(integ1.u))
  i = 0
  tvector = dt:dt:T

  for τ in tvector
    # evolve until rescaling:
    push!(integ1.opts.tstops, τ); step!(integ1)
    push!(integ2.opts.tstops, τ); step!(integ2)
    dist = norm(integ1.u .- integ2.u)
    # Rescale:
    if dist ≥ threshold
      # Rescale and reset everything:
      integ2.u .= integ1.u #.+ (integ2.u .- integ1.u)./a
      u_modified!(integ2, true)
      set_proposed_dt!(integ2,integ1)
      break
    end
  end

  τ = tvector[end]
  push!(integ1.opts.tstops, τ); step!(integ1)
  push!(integ2.opts.tstops, τ); step!(integ2)
  @test dist = norm(integ1.u .- integ2.u) == 0
end


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

# COPY

copy(corr::Homoscedastic)   = Homoscedastic(corr.adj, corr.method)
copy(corr::Heteroscedastic) = Heteroscedastic(corr.adj)

function copy(corr::CrossCorrelated)
    CrossCorrelated(corr.adj, copy(corr.nonmissing), copy(corr.mat))
end

function copy(corr::Clustered)
    Clustered(corr.adj, copy(corr.nonmissing), copy(corr.mat), copy(corr.ic), copy(corr.nc))
end

# EQUALITY

function Base.isequal(corr₁::C, corr₂::C) where {C <: CorrStructure}

    output = true

    for k in fieldnames(typeof(corr₁))
        output *= isequal(getfield(corr₁, k), getfield(corr₂, k))
    end

    return output
end

Base.isequal(corr₁::CorrStructure, corr₂::CorrStructure) = false

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

# DISTANCE ON THE GLOBE

havd(x::Float64) = 0.5 * (1.0 - cosd(x))

function geodistance(y0::Float64, x0::Float64, y1::Float64, x1::Float64)
    h = havd(y1 - y0) + havd(x1 - x0) * cosd(y0) * cosd(y1)
    return 12742.0 * asin(min(1.0, sqrt(h)))
end

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

# KERNELS

bartlett(x::Real)        = bartlett(float(x))
bartlett(x::Float64)     = (abs(x) >= 1.0) ? 0.0 : (1.0 - abs(x))
truncated(x::Real)       = truncated(float(x))
truncated(x::Float64)    = (abs(x) >= 1.0) ? 0.0 : 1.0
tukeyhanning(x::Real)    = tukeyhanning(float(x))
tukeyhanning(x::Float64) = (abs(x) >= 1.0) ? 0.0 : 0.5 * (1.0 + cospi(x))

parzen(x::Real) = parzen(float(x))

function parzen(x::Float64)
    if abs(x) <= 0.5
        return 1.0 - 6.0 * abs2(x) + 6 * abs(x)^3
    elseif abs(x) <= 1.0
        return 2.0 * (1.0 - abs(x))^3
    else
        return 0.0
    end
end

gallant(x::Real)    = parzen(float(x))
gallant(x::Float64) = parzen(x)


text:
export nmf

using mlpack._Internal.io

import mlpack_jll
const nmfLibrary = mlpack_jll.libmlpack_julia_nmf

# Call the C binding of the mlpack nmf binding.
function nmf_mlpackMain()
  success = ccall((:nmf, nmfLibrary), Bool, ())
  if !success
    # Throw an exception---false means there was a C++ exception.
    throw(ErrorException("mlpack binding error; see output"))
  end
end

" Internal module to hold utility functions. "
module nmf_internal
  import ..nmfLibrary

end # module

"""
    nmf(input, rank; [initial_h, initial_w, max_iterations, min_residue, seed, update_rules, verbose])

This program performs non-negative matrix factorization on the given dataset,
storing the resulting decomposed matrices in the specified files.  For an input
dataset V, NMF decomposes V into two matrices W and H such that

V = W * H

where all elements in W and H are non-negative.  If V is of size (n x m), then W
will be of size (n x r) and H will be of size (r x m), where r is the rank of
the factorization (specified by the `rank` parameter).

Optionally, the desired update rules for each NMF iteration can be chosen from
the following list:

 - multdist: multiplicative distance-based update rules (Lee and Seung 1999)
 - multdiv: multiplicative divergence-based update rules (Lee and Seung 1999)
 - als: alternating least squares update rules (Paatero and Tapper 1994)

The maximum number of iterations is specified with `max_iterations`, and the
minimum residue required for algorithm termination is specified with the
`min_residue` parameter.

For example, to run NMF on the input matrix `V` using the 'multdist' update
rules with a rank-10 decomposition and storing the decomposed matrices into `W`
and `H`, the following command could be used:

```julia
julia> using CSV
julia> V = CSV.read("V.csv")
julia> H, W = nmf(V, 10; update_rules="multdist")
```

# Arguments

 - `input::Array{Float64, 2}`: Input dataset to perform NMF on.
 - `rank::Int`: Rank of the factorization.
 - `initial_h::Array{Float64, 2}`: Initial H matrix.
 - `initial_w::Array{Float64, 2}`: Initial W matrix.
 - `max_iterations::Int`: Number of iterations before NMF terminates (0
      runs until convergence.  Default value `10000`.

 - `min_residue::Float64`: The minimum root mean square residue allowed
      for each iteration, below which the program terminates.  Default value
      `1e-05`.

 - `seed::Int`: Random seed.  If 0, 'std::time(NULL)' is used.  Default
      value `0`.

 - `update_rules::String`: Update rules for each iteration; ( multdist |
      multdiv | als ).  Default value `multdist`.

 - `verbose::Bool`: Display informational messages and the full list of
      parameters and timers at the end of execution.  Default value `false`.

# Return values

 - `h::Array{Float64, 2}`: Matrix to save the calculated H to.
 - `w::Array{Float64, 2}`: Matrix to save the calculated W to.

"""
function nmf(input,
             rank::Int;
             initial_h = missing,
             initial_w = missing,
             max_iterations::Union{Int, Missing} = missing,
             min_residue::Union{Float64, Missing} = missing,
             seed::Union{Int, Missing} = missing,
             update_rules::Union{String, Missing} = missing,
             verbose::Union{Bool, Missing} = missing,
             points_are_rows::Bool = true)
  # Force the symbols to load.
  ccall((:loadSymbols, nmfLibrary), Nothing, ());

  IORestoreSettings("Non-negative Matrix Factorization")

  # Process each input argument before calling mlpackMain().
  IOSetParamMat("input", input, points_are_rows)
  IOSetParam("rank", rank)
  if !ismissing(initial_h)
    IOSetParamMat("initial_h", initial_h, points_are_rows)
  end
  if !ismissing(initial_w)
    IOSetParamMat("initial_w", initial_w, points_are_rows)
  end
  if !ismissing(max_iterations)
    IOSetParam("max_iterations", convert(Int, max_iterations))
  end
  if !ismissing(min_residue)
    IOSetParam("min_residue", convert(Float64, min_residue))
  end
  if !ismissing(seed)
    IOSetParam("seed", convert(Int, seed))
  end
  if !ismissing(update_rules)
    IOSetParam("update_rules", convert(String, update_rules))
  end
  if verbose !== nothing && verbose === true
    IOEnableVerbose()
  else
    IODisableVerbose()
  end

  IOSetPassed("h")
  IOSetPassed("w")
  # Call the program.
  nmf_mlpackMain()

  return IOGetParamMat("h", points_are_rows),
         IOGetParamMat("w", points_are_rows)
end